(57)【要約】【課題】弾性波の波群の解析について、高精度かつリアルタイムに解析することを可能とした弾性波波群解析方法を提案する。【解決手段】振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、３種類のペアのそれぞれについて最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、３種類の回帰直線を用いて３種類の三次元合成位置ベクトルを作成するベクトル作成作業とを備え、異なる時間窓幅において波形データ抽出作業、回帰直線算出作業およびベクトル作成作業を繰り返すことで、算定時間毎に三次元合成位置ベクトルを作成し、三次元合成位置ベクトルを時系列で表示する弾性波波群解析方法。【選択図】図４
2014224748.tif 000002
【特許請求の範囲】
【請求項１】
振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、
前記三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の前記波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、
前記３種類のペアのそれぞれについて最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、
前記３種類の回帰直線を用いて３種類の三次元合成位置ベクトルを作成するベクトル作成作業と、を備え、
異なる時間窓幅において前記波形データ抽出作業、前記回帰直線算出作業および前記ベクトル作成作業を繰り返すことで、前記算定時間毎に前記三次元合成位置ベクトルを作成し、前記三次元合成位置ベクトルを時系列で表示することを特徴とする、弾性波波群解析方法。
【請求項２】
振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、
前記三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の前記波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、
前記３種類のペアのそれぞれに対して最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、
前記３種類の回帰直線のそれぞれについて回帰直線上に設けた回帰ベクトルのベクトル長を算出するベクトル長算出作業と、を備え、
異なる時間窓幅において前記波形データ抽出作業、前記回帰直線算出作業および前記ベクトル長算出作業を繰り返すことで、前記算定時間毎に前記ベクトル長を算出し、前記ベクトル長を時系列で表示することを特徴とする、弾性波波群解析方法。
【請求項３】
前記ベクトル長算出作業は、前記波形データのペアのｘ座標の最大値ｘmaxおよび最小値ｘminを算出し、
前記回帰直線に前記ｘ座標の最大値ｘmaxおよび最小値ｘminを代入して、ｙ座標の最大値ｙmaxおよびｙ座標の最小値ｙminを算出し、
前記回帰直線上の第一の点（ｘmax，ｙmax）と第二の点（ｘmin，ｙmin）との距離を前記ベクトル長とすることを特徴とする、請求項２に記載の弾性波波群解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、弾性波波群解析方法に関する。
【背景技術】
【０００２】
弾性波の波群を適切に解析することができれば、地震時の防災等に活用する他、トンネル施工における切羽前方探査等、建設工事等の地盤調査にも活用することができる。
【０００３】
弾性波を用いて地盤内を詳細に調査する場合には、Ｐ波の初動とＳ波の初動を読み取る必要がある。
Ｐ波の初動の自動読み取りは、アルゴリズムがほぼ確立しているものの、Ｓ波の初動の自動読み取りについては、確立しているとはいえない。
そのため、Ｓ波の初動の読み取りは、解析者が目視により判別しているのが一般的である。
【０００４】
ところが、目視によるＳ波の初動の読み取りでは、複数の解析者がいる場合、解析者の経験や主観の相違等により、読み取り誤差が生じるおそれがある。
そのため、読み取り誤差の低減化を目的として、Ｐ波とＳ波の読み取りに関する研究が行われてきた。
【０００５】
例えば、特許文献１には、過去に発生した多数の波形データと、過去の解析者による波形データのＰ波、Ｓ波到着時刻の読み取りデータ、および、それぞれの波形データについての振幅や周波数が変化する複数の時間区間におけるＳ／Ｎ比、卓越周波数、水平動と上下動の振幅比、区間の時間長、区間内での振幅変化等のデータを用いて算出された評価関数を利用して、地震波からノイズを分離し、Ｐ波とＳ波の到着時刻を決定する方法が開示されている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１２−１２７９２９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
特許文献１に記載のＰ波とＳ波の到着時刻を決定する方法は、膨大なデータに基づいて評価関数を作成する必要があり、データの管理および処理に手間がかかる。
また、地山状況は現場毎に異なっているため、過去のデータから導き出すことで誤差が生じるおそれがある。
【０００８】
このような観点から、本発明は、弾性波の解析について、読み取り誤差を最小限に抑えることを可能とした弾性波波群解析方法を提案することを課題とする。
【課題を解決するための手段】
【０００９】
前記課題を解決するために、本発明の弾性波波群解析方法は、振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、前記三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の前記波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、前記３種類のペアのそれぞれについて最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、前記３種類の回帰直線を用いて３種類の三次元合成位置ベクトルを作成するベクトル作成作業とを備え、異なる時間窓幅において（時間窓幅を一定時間ずらしながら）前記波形データ抽出作業、前記回帰直線算出作業および前記ベクトル作成作業を繰り返すことで、前記算定時間毎に前記三次元合成位置ベクトルを作成し、前記三次元合成位置ベクトルを時系列で表示することを特徴としている。
【００１０】
かかる弾性波波群解析方法によれば、三次元合成位置ベクトルの大きさや方向によりＰ波およびＳ波を確認することができるため、Ｐ波、Ｓ波の初動の読み取りの自動化が可能であるとともに、解析者の経験や主観等に影響されることなく、簡易かつ高精度にＰ波、Ｓ波の初動を推定することができる。
【００１１】
また、第二発明の弾性波波群解析方法は、振動方向が異なる三つの弾性波形記録を取得する波形データ取得作業と、前記三つの弾性波形記録のそれぞれについて所定の時間窓幅の波形データを抽出し、三成分の前記波形データから二成分ごとのペアを３種類設定する波形データ抽出作業と、前記３種類のペアのそれぞれに対して最小二乗法を適用し、３種類の回帰直線を算出する回帰直線算出作業と、前記３種類の回帰直線のそれぞれについて回帰直線上に設けた回帰ベクトルのベクトル長を算出するベクトル長算出作業とを備え、異なる時間窓幅において前記波形データ抽出作業、前記回帰直線算出作業および前記ベクトル長算出作業を繰り返すことで、前記算定時間毎に前記ベクトル長を算出し、前記ベクトル長を時系列で表示することを特徴としている。
【００１２】
かかる弾性波波群解析方法によれば、ベクトル長の大きさによりＰ波およびＳ波を確認することができるため、Ｐ波、Ｓ波の初動の読み取りの自動化に適しているとともに、解析者の経験や主観等に影響されることなく、簡易かつ高精度にＰ波、Ｓ波の初動を推定することができる。
【００１３】
前記ベクトル長算出作業では、前記波形データのペアのｘ座標の最大値ｘmaxおよび最小値ｘminを算出し、前記回帰直線に前記ｘ座標の最大値ｘmaxおよび最小値ｘminを代入してｙ座標の最大値ｙmaxおよびｙ座標の最小値ｙminを算出し、前記回帰直線上の第一の点（ｘmax，ｙmax）と第二の点（ｘmin，ｙmin）との距離を前記ベクトル長とすればよい。
【発明の効果】
【００１４】
本発明の弾性波波群解析方法によれば、弾性波の波群の解析を、高精度かつリアルタイムに実施することが可能となる。
【図面の簡単な説明】
【００１５】
【図１】本発明の実施形態に係る弾性波波群解析方法の作業手順を示すフローチャート図である。
【図２】図１の弾性波波群解析方法の弾性波形記録と算定時間との関係を示す模式図である。
【図３】（ａ）〜（ｃ）は二成分ごとのペアのパーティクル・オービットである。
【図４】（ａ）はパーティクル・オービットと回帰直線との関係を示すグラフ、（ｂ）は回帰直線のベクトル長を示すグラフ、（ｃ）は回帰直線の傾きを示すグラフである。
【図５】（ａ）〜（ｃ）はベクトル長を時刻暦で表示したグラフ、（ｄ）はＮＳ方向の弾性波形記録、（ｅ）はＥＷ方向の弾性波形記録、（ｆ）はＵＤ方向の弾性波形記録である。
【図６】（ａ）〜（ｃ）はベクトルの方向を時刻暦で表示したグラフ、（ｄ）はＮＳ方向の弾性波形記録、（ｅ）はＥＷ方向の弾性波形記録、（ｆ）はＵＤ方向の弾性波形記録である。
【図７】三次元位置ベクトルを時刻暦で表示したグラフである。
【図８】弾性波波群解析システムの概要を示す模式図である。
【図９】図８の弾性波波群解析システムのコンピュータのモニターを示す正面図である。
【発明を実施するための形態】
【００１６】
本発明の実施形態に係る弾性波波群解析方法は、図１に示すように、波形データ取得作業Ｓ１と、波形データ抽出作業Ｓ２と、回帰直線算出作業Ｓ３と、ベクトル長算出作業Ｓ４と、ベクトル作成作業Ｓ５とを備えている。
本実施形態では、三成分地震計１、データロガー２およびコンピュータ３を備えた弾性波波群解析システム（図８、９参照）により弾性波の波群の解析を行う。
【００１７】
波形データ取得作業Ｓ１では、振動方向が異なる三つの弾性波形記録ＵＤ，ＮＳ，ＥＷを取得する（図２参照）。
なお、弾性波形記録ＵＤは、地震動の上下方向（鉛直方向）の波形データであり、弾性波形記録ＮＳは地震動の南北方向の波形データであり、弾性波形記録ＥＷは地震動の東西方向の波形データである。
【００１８】
三つの弾性波形記録ＵＤ，ＮＳ，ＥＷは、三成分地震計１により計測した地震データを、プレトリガー付きデータロガー２を介してコンピュータ３（図８参照）に自動的に取り込むことで取得する。
三成分地震計は、上下および東西南北に合わせた状態で設置されている。
【００１９】
波形データ抽出作業Ｓ２では、三つの弾性波形記録ＵＤ，ＮＳ，ＥＷのそれぞれについて、弾性波形記録の収録時間（地震動の継続時間）よりも短い、所定の時間窓幅Δｔの波形データを抽出し、三成分の波形データから二成分ごとのペアを３種類設定する。
【００２０】
波形データの抽出は、所定の時間窓幅Δｔを設定したうえで、図２に示すように、弾性波形記録に対して、この時間窓幅Δｔ内の波形データを抽出することにより行う。
なお、時間窓幅Δｔの大きさは限定されるものではないが、例えば、データ（加速度、速度データ等）のサンプリングが０．０２ｍｓｅｃの場合、時間窓幅Δｔ＝０．２秒とする場合、パーティクル・オービットの回帰直線解析のデータ数は１１個となる。
【００２１】
波形データを抽出したら、二成分ごとのペア（ＮＳ−ＥＷ，ＵＤ−ＮＳ，ＵＤ−ＥＷ）を設定する。このとき、設定された各ペアについて、図３に示すようなパーティクル・オービットを作成する。
【００２２】
なお、図３の（ａ）は、時間窓幅Δｔに含まれる時刻ｔ１，ｔ２，…ｔｎに対応するＥＷ方向のデータｘ１，ｘ２，…ｘｎとＮＳ方向のデータｙ１，ｙ２，…ｙｎとを用いて作成した第一のペア（ＮＳ−ＥＷ）のパーティクル・オービットである。このパーティクル・オービットは、「ＥＷ」を横軸、「ＮＳ」を縦軸とした直交座標に、時刻ｔ１，ｔ２，…ｔｎにそれぞれ取得されたＥＷ成分およびＮＳ成分のペア（ｘ１，ｙ１）、（ｘ２，ｙ２）、…（ｘｎ，ｙｎ） をプロットし、各プロットを時刻順に線で結んで形成したものである。
同様に、図３の（ｂ）は、「ＮＳ」を横軸、「ＵＤ」を縦軸とした直交座標に、時刻ｔ１，ｔ２，…ｔｎにそれぞれ取得されたＮＳ成分およびＵＤ成分のペア（ｙ１，ｚ１）、（ｙ２，ｚ２）、…（ｙｎ，ｚｎ） をプロットし、各プロットを時刻順に線で結んで形成した第二のペア（ＵＤ−ＮＳ）のパーティクル・オービットであり、図３の（ｃ）は、「ＥＷ」を横軸、「ＵＤ」を縦軸とした直交座標に、時刻ｔ１，ｔ２，…ｔｎにそれぞれ取得されたＥＷ成分およびＵＤ成分のペア（ｘ１，ｚ１）、（ｘ２，ｚ２）、…（ｘｎ，ｚｎ） をプロットし、各プロットを時刻順に線で結んで形成した第三のペア（ＵＤ−ＥＷ）のパーティクル・オービットである。
【００２３】
回帰直線算出作業Ｓ３では、波形データ抽出作業Ｓ２において設定された３種類のペア（ＮＳ−ＥＷ，ＵＤ−ＮＳ，ＵＤ−ＥＷ）のそれぞれについて３種類の回帰直線を算出する。
【００２４】
回帰直線の算出は、時間窓幅Δｔ分の３成分のうちの２成分のペアに対して、最小二乗法を適用することにより行う。
すなわち、時間窓幅Δｔ内の時刻ｔ１，ｔ２，…ｔｎにおけるＥＷ成分およびＮＳ成分のペア（ｘ１，ｙ１）、（ｘ２，ｙ２）、…（ｘｎ，ｙｎ）に対して最小二乗法を適用して回帰直線を求めるとともに、ＮＳ成分およびＵＤ成分のペア（ｙ１，ｚ１）、（ｙ２，ｚ２）、…（ｙｎ，ｚｎ）に対して最小二乗法を適用して回帰直線を求め、さらに、ＥＷ成分およびＵＤ成分のペア（ｘ１，ｚ１）、（ｘ２，ｚ２）、…（ｘｎ，ｚｎ）に対して最小二乗法を適用して回帰直線を求める。
【００２５】
ベクトル長算出作業Ｓ４（図１）では、回帰直線算出作業Ｓ３（図１）で算出された３種類の回帰直線のそれぞれについて、回帰直線上に設けた回帰ベクトルのベクトル長を算出する。
【００２６】
ベクトル長の算出方法は、図４を参照して、以下に示すように行う。なお、以下では、時間窓幅Δｔにおける第一のペア（ＵＤ−ＥＷ）を例に、ベクトル長の算出方法を説明するが、他の時間幅および他のペアにおけるベクトル長の算出方法も同様である。また、図４を参照してベクトル長の算出方法を説明するが、パーティクル・オービットの図示は必須の手順ではない。
第一のペア（ＵＤ−ＥＷ）におけるベクトル長を算出する場合には、先ず、時間窓幅ΔｔにおけるＥＷ成分の波形データｘ１，ｘ２…ｘｎの中から最大値ｘmaxおよび最小値ｘminを抽出する（図４の（ａ）参照）。
【００２７】
次に、回帰直線（式１）にｘ座標の最大値ｘmaxおよび最小値ｘminを代入して、ＵＤ成分の最大値ｙmaxおよびｙ座標の最小値ｙminを算出する。本実施形態では、回帰直線上の第一の点（ｘmax，ｙmax）から第二の点（ｘmin，ｙmin）に向かうベクトルを回帰ベクトルとする。
ｙ＝ａｘｙｘ＋ｂｘｙ  ・・・  式１
【００２８】
そして、式２を利用して、回帰直線上の第一の点（ｘmax，ｙmax）と第二の点（ｘmin，ｙmin）との距離Ｌｘｙを算出し、この算出された距離Ｌｘｙをベクトル長Ｌとする（図４の（ｂ）参照）。
【００２９】
【数１】
2014224748.tif 000003
【００３０】
ベクトル作成作業Ｓ５では、３種類の回帰直線を用いて３種類の三次元合成位置ベクトルを作成する。
【００３１】
３つのペア（ＮＳ−ＥＷ，ＵＤ−ＮＳ，ＵＤ−ＥＷ）のｉ番目の時間窓幅Δｔｉにおける回帰ベクトルの空間ベクトルを、それぞれiｐNSEW，iｐUDNS，iｐUDEWとした場合、各回帰ベクトルは以下のようになる。
【００３２】
iｐNSEW（iｘmax−iｘmin，iｙmax−iｙmin，０）
iｐUDNS（０，iｙmax−iｙmin， iｚmax−iｚmin）
iｐUDEW（iｘmax−iｘmin，０，iｚmax−iｚmin）
【００３３】
そして、式３に示すように、三次元の各回帰ベクトルを成分合成して、ｉ番目の時間窓幅Δｔｉにおける三次元合成位置ベクトル（波群到来方向ベクトル）を算出する。
iＰ＝iｐNSEW＋iｐUDNS＋iｐUDEW  ・・・  式３
【００３４】
三次元合成位置ベクトルを算出したら、時間窓幅Δｔとは異なる時間窓幅において、波形データ抽出作業Ｓ２、回帰直線算出作業Ｓ３、ベクトル長算出作業Ｓ４およびベクトル作成作業Ｓ５を繰り返すことで、時間窓幅Δｔ毎のベクトル長Ｌ、方向ベクトルおよび三次元合成位置ベクトルiＰを算出する。
【００３５】
次の時間窓幅Δｔｎ＋１の設定方法は、前回算出した時間窓幅Δｔｎからシフト量ｄｔだけずらした時間窓幅とすればよい。
なお、シフト量ｄｔの大きさは限定されるものではないが、例えば、時間窓幅Δｔ＝０．２秒と設定した場合には、解析のオーバーラップを考えて、例えば、シフト量ｄｔ＝０．１秒にすればよい。この場合、前回の時間窓幅の一部（後半部分）と今回の時間窓幅の一部（前半部分）とがラップすることになる。
【００３６】
複数の時間窓幅毎にベクトル長Ｌ、方向ベクトルおよび三次元合成位置ベクトルiＰを算出したら、必要に応じて、図５〜７に示すように、それぞれ時系列で表示する。
【００３７】
ベクトル長Ｌを時系列で表示すると、図５の（ａ）〜（ｃ）に示すように、ベクトル長Ｌに大きな変化点が生じる場合がある。
このベクトル長Ｌの変化点をＰ波とＳ波の変化点として推定することができる。
また、方向ベクトルを時系列で表示すると、図６の（ａ）〜（ｃ）のようになる。
【００３８】
また、三次元合成位置ベクトルを時系列で表示すると、図７のようになる。図７では、ＥＷ方向を時間軸としている。ｉＰは、ｉ番目の時間窓幅における三次元合成ベクトルである。三次元合成ベクトル１Ｐ〜５Ｐのベクトル長が小さいことから、未だ弾性波が到来していないことが分かる。また、三次元合成ベクトル６Ｐは、そのベクトル長が三次元合成ベクトル１Ｐ〜５Ｐのベクトル長に比べて明らかに大きくなっていることから、６番目の時間窓幅にＰ波が到来したことが分かる。さらに、１１番目の時間窓幅の三次元合成ベクトル１１Ｐは、１０番目の時間窓幅以前の三次元合成ベクトル１Ｐ〜１０Ｐとは大きく異なる方向に向かっているので、１１番目の時間窓幅にＳ波が到来したことが分かり、さらに、その後の三次元合成ベクトルのベクトル長が小さくなっていることから、Ｓ波が一旦収束したことが分かる。また、三次元合成ベクトル２１Ｐは、２０番目の時間窓幅の三次元合成ベクトルとは大きく異なる方向に向かっているので、２１番目の時間窓幅に表面波などが到来したことが分かる。
【００３９】
そのため、Ｐ波の到着からＳ波の到着までの時間差である「ＰＳ時間」の推定、震源距離の推定、震央の推定等が可能となり、地震の震源位置の速報に使用することができる。
また、トンネルの弾性波探査に採用することで、切羽前方の地山状況の把握に役立てることもできる。
【００４０】
なお、震源距離は、以下に示すように、簡易的には、大森公式により算出することができる。
震源距離をＤｓ、Ｐ波の速度をｖｐ、Ｓ波の速度をｖｓとすると、Ｐ波の到着にかかった時間ｔｐおよびＳ波の到着にかかった時間ｔｓは、それぞれ式４および式５のようになる。
【００４１】
ｔｐ＝Ｄｓ／ｖｐ  ・・・  式４
ｔｓ＝Ｄｓ／ｖｓ  ・・・  式５
【００４２】
よって、ＰＳ時間（初期微動継続時間）ｔｐｓは、式６のようになる。
【００４３】
【数２】
2014224748.tif 000004
【００４４】
式６より、震源距離Ｄｓは、式７により示すことができる。
【００４５】
【数３】
2014224748.tif 000005
【００４６】
このように、Ｐ波の方向がＵＤ−ＥＷ、ＵＤ−ＮＳで決定すれば、方向と震源距離Ｄｓからおおよその震源の場所を推定することができる。
【００４７】
以上、本発明に係る実施形態について説明した。しかし、本発明は、前述の実施形態に限られず、前記の各構成要素については、本発明の趣旨を逸脱しない範囲で、適宜変更が可能である。
【００４８】
なお、ｉ番目の時間窓幅における三次元合成位置ベクトル（波群到来方向ベクトル）は、式８により算出することができる。
三次元合成位置ベクトルｉＰの成分を（ｉｐｘ，ｉｐｙ，ｉｐｚ）とすると、三次元合成位置ベクトルのベクトル長Ｌは式８で算出でき、三次元合成位置ベクトルｉｐの方向余弦（ｉｌ，ｉｍ，ｉｎ）は式９で算出することができる。
【００４９】
【数４】
2014224748.tif 000006
【符号の説明】
【００５０】
Ｓ１  波形データ取得作業
Ｓ２  波形データ抽出作業
Ｓ３  回帰直線算出作業
Ｓ４  ベクトル長算出作業
Ｓ５  ベクトル作成作業
【図１】
2014224748.tif 000007
【図２】
2014224748.tif 000008
【図３】
2014224748.tif 000009
【図４】
2014224748.tif 000010
【図５】
2014224748.tif 000011
【図６】
2014224748.tif 000012
【図７】
2014224748.tif 000013
【図８】
2014224748.tif 000014
【図９】
2014224748.tif 000015