特許 6131098 弾性波波群解析方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6131098

(24)【登録日】2017年4月21日

(45)【発行日】2017年5月17日

(54)【発明の名称】弾性波波群解析方法

(51)【国際特許分類】

   G01V 1/28 20060101AFI20170508BHJP

【ＦＩ】

   G01V1/28

【請求項の数】3

【全頁数】12

(21)【出願番号】特願2013-103932(P2013-103932)

(22)【出願日】2013年5月16日

(65)【公開番号】特開2014-224748(P2014-224748A)

(43)【公開日】2014年12月4日

【審査請求日】2016年2月16日

(73)【特許権者】

【識別番号】000206211

【氏名又は名称】大成建設株式会社

(74)【代理人】

【識別番号】110001807

【氏名又は名称】特許業務法人磯野国際特許商標事務所

(74)【代理人】

【識別番号】100064414

【弁理士】

【氏名又は名称】磯野 道造

(74)【代理人】

【識別番号】100111545

【弁理士】

【氏名又は名称】多田 悦夫

(74)【代理人】

【識別番号】100129067

【弁理士】

【氏名又は名称】町田 能章

(72)【発明者】

【氏名】今井 博

【審査官】 田中 秀直

(56)【参考文献】

【文献】 特開２００４−２８６５２７（ＪＰ，Ａ）

【文献】 特開２００９−１０３６７２（ＪＰ，Ａ）

【文献】 特開２０１２−１２７９２９（ＪＰ，Ａ）

【文献】 特開２００１−２０７４３２（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｖ １／００−１３／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、
前記三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の前記波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、
３種類の前記ペアのそれぞれについて最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、
前記３種類の回帰直線を用いて三次元合成位置ベクトルを作成するベクトル作成作業と、を備え、
前記ペアは、前記所定の時間窓幅に含まれる時刻に対応する一の振動方向の成分と前記時刻に対応する他の振動方向の成分との組み合わせたものであって、
前記波形データ抽出作業では、前記一の振動方向をＹ軸、前記他の振動方向をＸ軸とした直交座標に、前記ペアをプロットするとともに、当該プロットを時刻順に線で結んでパーティクル・オービットを形成し、
前記回帰曲線算出作業では、３種類の前記ペアに対してそれぞれ形成した前記パーティクル・オービットについて、それぞれ回帰曲線を算出し、
異なる時間窓幅において前記波形データ抽出作業、前記回帰直線算出作業および前記ベクトル作成作業を繰り返すことで、前記算定時間毎に前記三次元合成位置ベクトルを作成し、前記三次元合成位置ベクトルを時系列で表示することを特徴とする、弾性波波群解析方法。

【請求項2】

振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、
前記三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の前記波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、
３種類の前記ペアのそれぞれに対して最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、
前記３種類の回帰直線のそれぞれについて回帰直線上に設けた回帰ベクトルのベクトル長を算出するベクトル長算出作業と、を備え、
前記ペアは、前記所定の時間窓幅に含まれる時刻に対応する一の振動方向の成分と前記時刻に対応する他の振動方向の成分との組み合わせたものであって、
前記波形データ抽出作業では、前記一の振動方向をＹ軸、前記他の振動方向をＸ軸とした直交座標に、前記ペアをプロットするとともに、当該プロットを時刻順に線で結んでパーティクル・オービットを形成し、
前記回帰曲線算出作業では、３種類の前記ペアに対してそれぞれ形成した前記パーティクル・オービットについて、それぞれ回帰曲線を算出し、
異なる時間窓幅において前記波形データ抽出作業、前記回帰直線算出作業および前記ベクトル長算出作業を繰り返すことで、前記算定時間毎に前記ベクトル長を算出し、前記ベクトル長を時系列で表示することを特徴とする、弾性波波群解析方法。

【請求項3】

前記ベクトル長算出作業は、前記他の振動方向の成分の中からｘ座標の最大値ｘ_maxおよび最小値ｘ_minを算出し、
前記回帰直線に前記ｘ座標の最大値ｘ_maxおよび最小値ｘ_minを代入して、ｙ座標の最大値ｙ_maxおよびｙ座標の最小値ｙ_minを算出し、
前記回帰直線上の第一の点（ｘ_max，ｙ_max）と第二の点（ｘ_min，ｙ_min）との距離を前記ベクトル長とすることを特徴とする、請求項２に記載の弾性波波群解析方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、弾性波波群解析方法に関する。

【背景技術】

【0002】

弾性波の波群を適切に解析することができれば、地震時の防災等に活用する他、トンネル施工における切羽前方探査等、建設工事等の地盤調査にも活用することができる。

【0003】

弾性波を用いて地盤内を詳細に調査する場合には、Ｐ波の初動とＳ波の初動を読み取る必要がある。
Ｐ波の初動の自動読み取りは、アルゴリズムがほぼ確立しているものの、Ｓ波の初動の自動読み取りについては、確立しているとはいえない。
そのため、Ｓ波の初動の読み取りは、解析者が目視により判別しているのが一般的である。

【0004】

ところが、目視によるＳ波の初動の読み取りでは、複数の解析者がいる場合、解析者の経験や主観の相違等により、読み取り誤差が生じるおそれがある。
そのため、読み取り誤差の低減化を目的として、Ｐ波とＳ波の読み取りに関する研究が行われてきた。

【0005】

例えば、特許文献１には、過去に発生した多数の波形データと、過去の解析者による波形データのＰ波、Ｓ波到着時刻の読み取りデータ、および、それぞれの波形データについての振幅や周波数が変化する複数の時間区間におけるＳ／Ｎ比、卓越周波数、水平動と上下動の振幅比、区間の時間長、区間内での振幅変化等のデータを用いて算出された評価関数を利用して、地震波からノイズを分離し、Ｐ波とＳ波の到着時刻を決定する方法が開示されている。

【先行技術文献】

【特許文献】

【0006】

【特許文献1】特開２０１２−１２７９２９号公報

【発明の概要】

【発明が解決しようとする課題】

【0007】

特許文献１に記載のＰ波とＳ波の到着時刻を決定する方法は、膨大なデータに基づいて評価関数を作成する必要があり、データの管理および処理に手間がかかる。
また、地山状況は現場毎に異なっているため、過去のデータから導き出すことで誤差が生じるおそれがある。

【0008】

このような観点から、本発明は、弾性波の解析について、読み取り誤差を最小限に抑えることを可能とした弾性波波群解析方法を提案することを課題とする。

【課題を解決するための手段】

【0009】

前記課題を解決するために、本発明の弾性波波群解析方法は、振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、前記三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の前記波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、３種類の前記ペアのそれぞれについて最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、前記３種類の回帰直線を用いて三次元合成位置ベクトルを作成するベクトル作成作業とを備え、前記ペアは、前記所定の時間窓幅に含まれる時刻に対応する一の振動方向の成分と前記時刻に対応する他の振動方向の成分との組み合わせたものであって、前記波形データ抽出作業では、前記一の振動方向をＹ軸、前記他の振動方向をＸ軸とした直交座標に、前記ペアをプロットするとともに、当該プロットを時刻順に線で結んでパーティクル・オービットを形成し、前記回帰曲線算出作業では、３種類の前記ペアに対してそれぞれ形成した前記パーティクル・オービットについて、それぞれ回帰曲線を算出し、異なる時間窓幅において（時間窓幅を一定時間ずらしながら）前記波形データ抽出作業、前記回帰直線算出作業および前記ベクトル作成作業を繰り返すことで、前記算定時間毎に前記三次元合成位置ベクトルを作成し、前記三次元合成位置ベクトルを時系列で表示することを特徴としている。

【0010】

かかる弾性波波群解析方法によれば、三次元合成位置ベクトルの大きさや方向によりＰ波およびＳ波を確認することができるため、Ｐ波、Ｓ波の初動の読み取りの自動化が可能であるとともに、解析者の経験や主観等に影響されることなく、簡易かつ高精度にＰ波、Ｓ波の初動を推定することができる。

【0011】

また、第二発明の弾性波波群解析方法は、振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、前記三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の前記波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、３種類の前記ペアのそれぞれに対して最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、前記３種類の回帰直線のそれぞれについて回帰直線上に設けた回帰ベクトルのベクトル長を算出するベクトル長算出作業とを備え、前記ペアは、前記所定の時間窓幅に含まれる時刻に対応する一の振動方向の成分と前記時刻に対応する他の振動方向の成分との組み合わせたものであって、前記波形データ抽出作業では、前記一の振動方向をＹ軸、前記他の振動方向をＸ軸とした直交座標に、前記ペアをプロットするとともに、当該プロットを時刻順に線で結んでパーティクル・オービットを形成し、前記回帰曲線算出作業では、３種類の前記ペアに対してそれぞれ形成した前記パーティクル・オービットについて、それぞれ回帰曲線を算出し、異なる時間窓幅において前記波形データ抽出作業、前記回帰直線算出作業および前記ベクトル長算出作業を繰り返すことで、前記算定時間毎に前記ベクトル長を算出し、前記ベクトル長を時系列で表示することを特徴としている。

【0012】

かかる弾性波波群解析方法によれば、ベクトル長の大きさによりＰ波およびＳ波を確認することができるため、Ｐ波、Ｓ波の初動の読み取りの自動化に適しているとともに、解析者の経験や主観等に影響されることなく、簡易かつ高精度にＰ波、Ｓ波の初動を推定することができる。

【0013】

前記ベクトル長算出作業では、前記他の振動方向の成分の中からｘ座標の最大値ｘ_maxおよび最小値ｘ_minを算出し、前記回帰直線に前記ｘ座標の最大値ｘ_maxおよび最小値ｘ_minを代入してｙ座標の最大値ｙ_maxおよびｙ座標の最小値ｙ_minを算出し、前記回帰直線上の第一の点（ｘ_max，ｙ_max）と第二の点（ｘ_min，ｙ_min）との距離を前記ベクトル長とすればよい。

【発明の効果】

【0014】

本発明の弾性波波群解析方法によれば、弾性波の波群の解析を、高精度かつリアルタイムに実施することが可能となる。

【図面の簡単な説明】

【0015】

【図1】本発明の実施形態に係る弾性波波群解析方法の作業手順を示すフローチャート図である。

【図2】図１の弾性波波群解析方法の弾性波形記録と算定時間との関係を示す模式図である。

【図3】（ａ）〜（ｃ）は二成分ごとのペアのパーティクル・オービットである。

【図4】（ａ）はパーティクル・オービットと回帰直線との関係を示すグラフ、（ｂ）は回帰直線のベクトル長を示すグラフ、（ｃ）は回帰直線の傾きを示すグラフである。

【図5】（ａ）〜（ｃ）はベクトル長を時刻暦で表示したグラフ、（ｄ）はＮＳ方向の弾性波形記録、（ｅ）はＥＷ方向の弾性波形記録、（ｆ）はＵＤ方向の弾性波形記録である。

【図6】（ａ）〜（ｃ）はベクトルの方向を時刻暦で表示したグラフ、（ｄ）はＮＳ方向の弾性波形記録、（ｅ）はＥＷ方向の弾性波形記録、（ｆ）はＵＤ方向の弾性波形記録である。

【図7】三次元位置ベクトルを時刻暦で表示したグラフである。

【図8】弾性波波群解析システムの概要を示す模式図である。

【図9】図８の弾性波波群解析システムのコンピュータのモニターを示す正面図である。

【発明を実施するための形態】

【0016】

本発明の実施形態に係る弾性波波群解析方法は、図１に示すように、波形データ取得作業Ｓ１と、波形データ抽出作業Ｓ２と、回帰直線算出作業Ｓ３と、ベクトル長算出作業Ｓ４と、ベクトル作成作業Ｓ５とを備えている。
本実施形態では、三成分地震計１、データロガー２およびコンピュータ３を備えた弾性波波群解析システム（図８、９参照）により弾性波の波群の解析を行う。

【0017】

波形データ取得作業Ｓ１では、振動方向が異なる三つの弾性波形記録ＵＤ，ＮＳ，ＥＷを取得する（図２参照）。
なお、弾性波形記録ＵＤは、地震動の上下方向（鉛直方向）の波形データであり、弾性波形記録ＮＳは地震動の南北方向の波形データであり、弾性波形記録ＥＷは地震動の東西方向の波形データである。

【0018】

三つの弾性波形記録ＵＤ，ＮＳ，ＥＷは、三成分地震計１により計測した地震データを、プレトリガー付きデータロガー２を介してコンピュータ３（図８参照）に自動的に取り込むことで取得する。
三成分地震計は、上下および東西南北に合わせた状態で設置されている。

【0019】

波形データ抽出作業Ｓ２では、三つの弾性波形記録ＵＤ，ＮＳ，ＥＷのそれぞれについて、弾性波形記録の収録時間（地震動の継続時間）よりも短い、所定の時間窓幅Δｔの波形データを抽出し、三成分の波形データから二成分ごとのペアを３種類設定する。

【0020】

波形データの抽出は、所定の時間窓幅Δｔを設定したうえで、図２に示すように、弾性波形記録に対して、この時間窓幅Δｔ内の波形データを抽出することにより行う。
なお、時間窓幅Δｔの大きさは限定されるものではないが、例えば、データ（加速度、速度データ等）のサンプリングが０．０２ｍｓｅｃの場合、時間窓幅Δｔ＝０．２秒とする場合、パーティクル・オービットの回帰直線解析のデータ数は１１個となる。

【0021】

波形データを抽出したら、二成分ごとのペア（ＮＳ−ＥＷ，ＵＤ−ＮＳ，ＵＤ−ＥＷ）を設定する。このとき、設定された各ペアについて、図３に示すようなパーティクル・オービットを作成する。

【0022】

なお、図３の（ａ）は、時間窓幅Δｔに含まれる時刻ｔ_１，ｔ_２，…ｔ_ｎに対応するＥＷ方向のデータｘ_１，ｘ_２，…ｘ_ｎとＮＳ方向のデータｙ_１，ｙ_２，…ｙ_ｎとを用いて作成した第一のペア（ＮＳ−ＥＷ）のパーティクル・オービットである。このパーティクル・オービットは、「ＥＷ」を横軸、「ＮＳ」を縦軸とした直交座標に、時刻ｔ_１，ｔ_２，…ｔ_ｎにそれぞれ取得されたＥＷ成分およびＮＳ成分のペア（ｘ_１，ｙ_１）、（ｘ_２，ｙ_２）、…（ｘ_ｎ，ｙ_ｎ） をプロットし、各プロットを時刻順に線で結んで形成したものである。
同様に、図３の（ｂ）は、「ＮＳ」を横軸、「ＵＤ」を縦軸とした直交座標に、時刻ｔ_１，ｔ_２，…ｔ_ｎにそれぞれ取得されたＮＳ成分およびＵＤ成分のペア（ｙ_１，ｚ_１）、（ｙ_２，ｚ_２）、…（ｙ_ｎ，ｚ_ｎ） をプロットし、各プロットを時刻順に線で結んで形成した第二のペア（ＵＤ−ＮＳ）のパーティクル・オービットであり、図３の（ｃ）は、「ＥＷ」を横軸、「ＵＤ」を縦軸とした直交座標に、時刻ｔ_１，ｔ_２，…ｔ_ｎにそれぞれ取得されたＥＷ成分およびＵＤ成分のペア（ｘ_１，ｚ_１）、（ｘ_２，ｚ_２）、…（ｘ_ｎ，ｚ_ｎ） をプロットし、各プロットを時刻順に線で結んで形成した第三のペア（ＵＤ−ＥＷ）のパーティクル・オービットである。

【0023】

回帰直線算出作業Ｓ３では、波形データ抽出作業Ｓ２において設定された３種類のペア（ＮＳ−ＥＷ，ＵＤ−ＮＳ，ＵＤ−ＥＷ）のそれぞれについて３種類の回帰直線を算出する。

【0024】

回帰直線の算出は、時間窓幅Δｔ分の３成分のうちの２成分のペアに対して、最小二乗法を適用することにより行う。
すなわち、時間窓幅Δｔ内の時刻ｔ_１，ｔ_２，…ｔ_ｎにおけるＥＷ成分およびＮＳ成分のペア（ｘ_１，ｙ_１）、（ｘ_２，ｙ_２）、…（ｘ_ｎ，ｙ_ｎ）に対して最小二乗法を適用して回帰直線を求めるとともに、ＮＳ成分およびＵＤ成分のペア（ｙ_１，ｚ_１）、（ｙ_２，ｚ_２）、…（ｙ_ｎ，ｚ_ｎ）に対して最小二乗法を適用して回帰直線を求め、さらに、ＥＷ成分およびＵＤ成分のペア（ｘ_１，ｚ_１）、（ｘ_２，ｚ_２）、…（ｘ_ｎ，ｚ_ｎ）に対して最小二乗法を適用して回帰直線を求める。

【0025】

ベクトル長算出作業Ｓ４（図１）では、回帰直線算出作業Ｓ３（図１）で算出された３種類の回帰直線のそれぞれについて、回帰直線上に設けた回帰ベクトルのベクトル長を算出する。

【0026】

ベクトル長の算出方法は、図４を参照して、以下に示すように行う。なお、以下では、時間窓幅Δｔにおける第一のペア（ＵＤ−ＥＷ）を例に、ベクトル長の算出方法を説明するが、他の時間幅および他のペアにおけるベクトル長の算出方法も同様である。また、図４を参照してベクトル長の算出方法を説明するが、パーティクル・オービットの図示は必須の手順ではない。
第一のペア（ＵＤ−ＥＷ）におけるベクトル長を算出する場合には、先ず、時間窓幅ΔｔにおけるＥＷ成分の波形データｘ_１，ｘ_２…ｘ_ｎの中から最大値ｘ_maxおよび最小値ｘ_minを抽出する（図４の（ａ）参照）。

【0027】

次に、回帰直線（式１）にｘ座標の最大値ｘ_maxおよび最小値ｘ_minを代入して、ＵＤ成分の最大値ｙ_maxおよびｙ座標の最小値ｙ_minを算出する。本実施形態では、回帰直線上の第一の点（ｘ_max，ｙ_max）から第二の点（ｘ_min，ｙ_min）に向かうベクトルを回帰ベクトルとする。

ｙ＝ａ_ｘｙｘ＋ｂ_ｘｙ ・・・ 式１

【0028】

そして、式２を利用して、回帰直線上の第一の点（ｘ_max，ｙ_max）と第二の点（ｘ_min，ｙ_min）との距離Ｌ_ｘｙを算出し、この算出された距離Ｌ_ｘｙをベクトル長Ｌとする（図４の（ｂ）参照）。

【0029】

【数1】

【0030】

ベクトル作成作業Ｓ５では、３種類の回帰直線を用いて三次元合成位置ベクトルを作成する。

【0031】

３つのペア（ＮＳ−ＥＷ，ＵＤ−ＮＳ，ＵＤ−ＥＷ）のｉ番目の時間窓幅Δｔ_ｉにおける回帰ベクトルの空間ベクトルを、それぞれⁱｐ_NSEW，ⁱｐ_UDNS，ⁱｐ_UDEWとした場合、各回帰ベクトルは以下のようになる。

【0032】

ⁱｐ_NSEW（ⁱｘ_max−ⁱｘ_min，ⁱｙ_max−ⁱｙ_min，０）
ⁱｐ_UDNS（０，ⁱｙ_max−ⁱｙ_min， ⁱｚ_max−ⁱｚ_min）
ⁱｐ_UDEW（ⁱｘ_max−ⁱｘ_min，０，ⁱｚ_max−ⁱｚ_min）

【0033】

そして、式３に示すように、三次元の各回帰ベクトルを成分合成して、ｉ番目の時間窓幅Δｔ_ｉにおける三次元合成位置ベクトル（波群到来方向ベクトル）を算出する。

ⁱＰ＝ⁱｐ_NSEW＋ⁱｐ_UDNS＋ⁱｐ_UDEW ・・・ 式３

【0034】

三次元合成位置ベクトルを算出したら、時間窓幅Δｔとは異なる時間窓幅において、波形データ抽出作業Ｓ２、回帰直線算出作業Ｓ３、ベクトル長算出作業Ｓ４およびベクトル作成作業Ｓ５を繰り返すことで、時間窓幅Δｔ毎のベクトル長Ｌ、方向ベクトルおよび三次元合成位置ベクトルⁱＰを算出する。

【0035】

次の時間窓幅Δｔ_ｎ＋１の設定方法は、前回算出した時間窓幅Δｔ_ｎからシフト量ｄｔだけずらした時間窓幅とすればよい。
なお、シフト量ｄｔの大きさは限定されるものではないが、例えば、時間窓幅Δｔ＝０．２秒と設定した場合には、解析のオーバーラップを考えて、例えば、シフト量ｄｔ＝０．１秒にすればよい。この場合、前回の時間窓幅の一部（後半部分）と今回の時間窓幅の一部（前半部分）とがラップすることになる。

【0036】

複数の時間窓幅毎にベクトル長Ｌ、方向ベクトルおよび三次元合成位置ベクトルⁱＰを算出したら、必要に応じて、図５〜７に示すように、それぞれ時系列で表示する。

【0037】

ベクトル長Ｌを時系列で表示すると、図５の（ａ）〜（ｃ）に示すように、ベクトル長Ｌに大きな変化点が生じる場合がある。
このベクトル長Ｌの変化点をＰ波とＳ波の変化点として推定することができる。
また、方向ベクトルを時系列で表示すると、図６の（ａ）〜（ｃ）のようになる。

【0038】

また、三次元合成位置ベクトルを時系列で表示すると、図７のようになる。図７では、ＥＷ方向を時間軸としている。^ｉＰは、ｉ番目の時間窓幅における三次元合成ベクトルである。三次元合成ベクトル^１Ｐ〜^５Ｐのベクトル長が小さいことから、未だ弾性波が到来していないことが分かる。また、三次元合成ベクトル^６Ｐは、そのベクトル長が三次元合成ベクトル^１Ｐ〜^５Ｐのベクトル長に比べて明らかに大きくなっていることから、６番目の時間窓幅にＰ波が到来したことが分かる。さらに、１１番目の時間窓幅の三次元合成ベクトル^１１Ｐは、１０番目の時間窓幅以前の三次元合成ベクトル^１Ｐ〜^１０Ｐとは大きく異なる方向に向かっているので、１１番目の時間窓幅にＳ波が到来したことが分かり、さらに、その後の三次元合成ベクトルのベクトル長が小さくなっていることから、Ｓ波が一旦収束したことが分かる。また、三次元合成ベクトル^２１Ｐは、２０番目の時間窓幅の三次元合成ベクトルとは大きく異なる方向に向かっているので、２１番目の時間窓幅に表面波などが到来したことが分かる。

【0039】

そのため、Ｐ波の到着からＳ波の到着までの時間差である「ＰＳ時間」の推定、震源距離の推定、震央の推定等が可能となり、地震の震源位置の速報に使用することができる。
また、トンネルの弾性波探査に採用することで、切羽前方の地山状況の把握に役立てることもできる。

【0040】

なお、震源距離は、以下に示すように、簡易的には、大森公式により算出することができる。
震源距離をＤ_ｓ、Ｐ波の速度をｖ_ｐ、Ｓ波の速度をｖ_ｓとすると、Ｐ波の到着にかかった時間ｔ_ｐおよびＳ波の到着にかかった時間ｔ_ｓは、それぞれ式４および式５のようになる。

【0041】

ｔ_ｐ＝Ｄ_ｓ／ｖ_ｐ ・・・ 式４
ｔ_ｓ＝Ｄ_ｓ／ｖ_ｓ ・・・ 式５

【0042】

よって、ＰＳ時間（初期微動継続時間）ｔ_ｐｓは、式６のようになる。

【0043】

【数2】

【0044】

式６より、震源距離Ｄ_ｓは、式７により示すことができる。

【0045】

【数3】

【0046】

このように、Ｐ波の方向がＵＤ−ＥＷ、ＵＤ−ＮＳで決定すれば、方向と震源距離Ｄ_ｓからおおよその震源の場所を推定することができる。

【0047】

以上、本発明に係る実施形態について説明した。しかし、本発明は、前述の実施形態に限られず、前記の各構成要素については、本発明の趣旨を逸脱しない範囲で、適宜変更が可能である。

【0048】

なお、ｉ番目の時間窓幅における三次元合成位置ベクトル（波群到来方向ベクトル）は、式８により算出することができる。
三次元合成位置ベクトル^ｉＰの成分を（^ｉｐ_ｘ，^ｉｐ_ｙ，^ｉｐ_ｚ）とすると、三次元合成位置ベクトルのベクトル長Ｌは式８で算出でき、三次元合成位置ベクトル^ｉｐの方向余弦（^ｉｌ，^ｉｍ，^ｉｎ）は式９で算出することができる。

【0049】

【数4】

【符号の説明】

【0050】

Ｓ１ 波形データ取得作業
Ｓ２ 波形データ抽出作業
Ｓ３ 回帰直線算出作業
Ｓ４ ベクトル長算出作業
Ｓ５ ベクトル作成作業

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】