JP6460523B2

JP6460523B2 - 地下構造探査システム及び地下構造探査方法

Info

Publication number: JP6460523B2
Application number: JP2015010738A
Authority: JP
Inventors: 宏光水谷; 研志河合; 芳伸狩野
Original assignee: 株式会社セオコンプ
Priority date: 2015-01-22
Filing date: 2015-01-22
Publication date: 2019-01-30
Anticipated expiration: 2035-01-22
Also published as: JP2016133497A

Description

本発明は、地震波を利用して地下の構造を探査する地下構造探査システム及び地下構造探査方法に関する。

従来、地中に含まれる石油，ガス，石炭，鉱物などの資源の探査、あるいは地中の断層などの構造の探査を行うために、地震波を使用した探査が行われている。この場合の探査に使用する地震波としては、実際の地震波でもよいが、通常は人工的に発生させた地震による地震波が使用される。

地震波を利用した解析としては、フォワード問題（順問題）による解析とインバース問題（逆問題）による解析がある。フォワード問題による解析は、震源や地下構造が予め判っているときに、震源からの地震波の地中での拡散状態を計算して、地震波動場を計算するものである。フォワード問題による解析では、例えば図８Ａに示すように、地下構造Ｘａが既知であるとき、震源ｅから地震波を発生させた場合の各地点ａ１，ａ２，ａ３，ａ４での地震波の伝わり状況を計算する。

一方、インバース問題による解析は、観測地点で観測された地震波形を利用して、実際の地震波の動きとは逆の方向に解析を進めて、震源や地下構造のモデルを推定するものである。インバース問題による解析では、例えば図８Ｂに示すように、震源ｅで発生させた地震波を各地点ａ１，ａ２，ａ３，ａ４で観測し、各地点ａ１，ａ２，ａ３，ａ４での観測結果から、未知の地下構造Ｘｂあるいは震源ｅ、もしくはその両方を算出するものである。地震波形を用いたインバース問題は、波形インバージョンと称される。

特許文献１には、先に提案した地震波を使ったシミュレーション手法の一例が開示されている。

特開２０１０−１２３０５６号公報

波形インバージョン法を用いた地下構造の探査では、観測点で収録された波形そのものをデータとして用いることにより波形の持つ情報を最大限引き出すことができる。そのため、既存の地下構造推定手法よりも精度および解像度が高いという優位性を持つ。しかしながら、波形インバージョンを行う際には、理論及び技術上の困難さのみならず、計算機での計算負担が極めて大きいという問題があった。計算機の性能向上および理論及び技術的な進展により、近年ようやく実現できるようになった手法である。

このように波形インバージョンによる地下構造の探査は、従来の他の手法に比べて優位性を持つ一方、計算負担が極めて大きいという問題がある。また、波形インバージョン法で用いる順伝播計算および逆伝播計算手法の精度及び効率に関する進歩により、求まる地下構造の精度および解像度の大幅な向上が望まれていた。

本発明は、少ない計算負担で、精度および解像度の高い地下構造の推定ができる地下構造探査システム及び地下構造探査方法を提供することを目的とする。

本発明の地下構造探査システムは、地球内部を伝播し、観測された地震波形データを解析することで、地球内部の地下空間の構造を探査する地下構造探査システムである。
本発明の地下構造探査システムは、記憶部と、伝播波動場計算部と、モデル修正部と、収束判断部と、出力部とを備える。
記憶部には、地震波の震源と推定される地下構造の初期モデルと観測点との情報とを含む入力データ又はその初期モデルを修正した中間モデルが記憶される。
伝播波動場計算部は、記憶部に記憶された初期モデル又は中間モデルに対して、離散化された時間領域の運動方程式を使用して、過去から未来に向けた順方向の伝播波動場を演算すると共に、未来から過去に向けた逆方向の伝播波動場を演算する。
モデル修正部は、伝播波動場計算部で演算された順方向の伝播波動場と逆方向の伝播波動場を利用して初期モデル又は中間モデルの修正量を算出し、算出された修正量の修正を記憶部に記憶された初期モデル又は中間モデルに対して施す。
収束判断部は、モデル修正部での修正状態の収束を判断する。
出力部は、収束判断部で収束を判断した場合に、モデル修正部で修正されて記憶部に記憶された中間モデルによる地下構造を、探査対象の地下構造として読み出す。
そして、伝播波動場計算部での演算に使用する運動方程式は、時間領域の反復法を用いて行う計算手法であって、密度に関する項の逆演算行列と、弾性定数に関する部分を離散化した行列とを、それぞれ正定値、半正定値の対称行列もしくはエルミート行列として定義して、相反定理及び時間反転対称性が成り立つように定義した運動方程式とした。

本発明の地下構造探査方法は、地球内部を伝播する地震波を解析することで、地球内部の地下空間の構造を探査する地下構造探査方法である。
本発明の地下構造探査方法は、モデル記憶ステップと、伝播波動場演算ステップと、モデル修正ステップと、収束判断ステップと、出力ステップと、を含む。
モデル記憶ステップは、地震波の震源と推定される地下構造の初期モデルと観測点との情報とを含む入力データ又はその初期モデルを修正した中間モデルを設定して記憶部が記憶する。
伝播波動場演算ステップは、モデル記憶ステップで記憶部に記憶された初期モデル又は中間モデルに対して、伝播波動場計算部により、離散化された時間領域の運動方程式を使用して、過去から未来に向けた順方向の伝播波動場を演算すると共に、未来から過去に向けた逆方向の伝播波動場を演算する。
モデル修正ステップは、モデル修正部により、伝播波動場演算ステップで演算された順方向の伝播波動場と逆方向の伝播波動場とを内積計算して、内積計算した結果から初期モデル又は中間モデルの修正量を算出し、算出された修正量の修正をモデル記憶ステップにより記憶された初期モデル又は中間モデルに対して施す。
収束判断ステップは、モデル修正ステップでの修正状態の収束を、集束判断部が判断する。
出力ステップは、収束判断ステップで収束したと判断された場合に、モデル修正ステップで修正されて記憶部に記憶された中間モデルによる地下構造を、探査対象の地下構造として確定して出力部が出力する。
そして、伝播波動場計算部による前記伝播波動場演算ステップでの演算に使用する運動方程式は、密度に関する項の逆演算行列と、弾性定数に関する部分を離散化した行列とを、それぞれ正定値、半正定値の対称行列もしくはエルミート行列として定義して、相反定理及び時間反転対称性が成り立つ演算子として使用した運動方程式とした。

本発明によると、地球内部を伝搬する地震波をシミュレーションする際に、順伝播と逆伝播の計算結果の相反性が保障され、推定されたモデルの信頼度を向上することができる。すなわち、順伝播と逆伝播が相反性を持たないことに起因する、予測できない誤差が混入するのを防止することができ、高精度の地下構造の探査が効率よくできる効果を有する。

本発明の一実施の形態例による地下構造探査処理の全体を示す構成図である。本発明の一実施の形態例による演算処理装置の構成例を示すブロック図である。本発明の一実施の形態例による地下構造探査処理の流れを示すフローチャートである。本発明の一実施の形態例によるモデル修正量の計算処理の流れを示すフローチャートである。本発明の一実施の形態例によるメッシュ毎の変動と観測波形との関係の例を示す説明図である。本発明の一実施の形態例による処理の概要を示す説明図である。順伝播と逆伝播の計算結果が非対称の場合（図７Ａ）と、対称の場合（図７Ｂ）とを比較した説明図である。フォワード問題による解析（図８Ａ）とインバース問題による解析（図８Ｂ）とを示す説明図である。

以下、本発明の一実施の形態例（以下、「本例」と称する。）を、図１〜図７を参照して説明する。
［１．全体の概要］
図１は、本例の全体構成例を示す図である。
本例では、震源ｋで発生した地震による地震波を複数の観測点Ｐ１，Ｐ２，Ｐ３で観測する。そして、各観測点Ｐ１，Ｐ２，Ｐ３で観測された地震波形データから、地下構造ｘが算出される。
震源データ収集部１１は、震源ｋで発生した地震の震源データを収集する。地震波形データ収集部１２は、各観測点Ｐ１，Ｐ２，Ｐ３で観測した地震波形データを収集する。初期モデル作成部１３は、地下構造ｘを推定したデータである、初期モデルを保持する。そして、震源データ収集部１１で収集した震源データと、地震波形データ収集部１２で収集した各観測点Ｐ１，Ｐ２，Ｐ３の地震波形データと、初期モデル作成部１３が作成した初期モデルとが、入力データ結合部１４に供給され、震源データや地震波形データを結合したデータセットと初期モデルとが得られる。

入力データ結合部１４で得られたセット（データセットと初期モデル）が、演算処理装置２０に供給される。演算処理装置２０は、演算処理により初期モデルを修正した中間モデルを生成し、さらにその中間モデルを修正する処理を繰り返し行う。演算処理装置２０での演算処理の詳細については後述する。
演算処理装置２０で生成した中間モデルが適正な状態となったとき、出力装置３０が、その中間モデルが持つデータを出力する。出力装置３０は、例えばディスプレイやプリンタを備え、中間モデルが持つデータの表示又は印刷を行う。あるいは、出力装置３０から別の装置に解析結果のデータが転送される。

［２．演算処理装置の構成］
図２は、演算処理装置２０の構成例を示す。
演算処理装置２０は、データ入力部２１を備える。データ入力部２１が受け付けたデータセットと初期モデルのデータがモデル記憶部２２に供給される。モデル記憶部２２が記憶した初期モデルのデータは、モデル修正部２６から修正された中間モデルのデータが供給される毎に更新される。

モデル記憶部２２に記憶されたモデルのデータは、順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４に供給される。
順方向伝播波動場計算部２３では、震源から観測点の方向への波動場の伝播状態が計算される。また、逆方向伝播波動場計算部２４では、実際の地震波の伝播方向とは逆の、震源から観測点の方向の伝播状態が計算される。

順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４が計算を行う際には、地下構造を一定の距離ごとに区切ったメッシュ集合を設定する。そして、順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４は、地下構造（密度と弾性定数）と波動場の状態をメッシュごとに離散化すると共に、時間刻みによって離散化し、それらを離散化した運動方程式に代入することで、各メッシュの地下構造を伝播する波動場を計算する。
順方向伝播波動場計算部２３は、過去と現在の２ステップの波動場の状態を用いて、１ステップ先の未来の波動場を計算することで、順々に時間をすすめていく順伝播の計算を行う。逆方向伝播波動場計算部２４は、未来と現在の２ステップの波動場の状態を用いて、過去のステップの波動場を計算することで、時間を順々に戻していく逆伝播の計算を行う。
なお、順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４での計算に使用する運動方程式については後述する。

そして、順方向伝播波動場計算部２３で得られた順方向の伝播波動場のデータと、逆方向伝播波動場計算部２４で得られた逆方向の伝播波動場のデータとが、内積計算部２５に供給される。
内積計算部２５は、順方向の伝播波動場のデータと逆方向の伝播波動場のデータとを使った内積計算を行う。そして、内積計算部２５で得られた内積計算結果が、モデル修正部２６に供給される。

モデル修正部２６は、内積計算部２５で得られた内積計算結果から、モデル記憶部２２が記憶したモデル（初期モデル又は中間モデル）の修正量を算出する。ここでは、計算された内積と、観測波形データと初期モデルまたは中間モデルに対して生成された波形データの残差とを用いて、モデルの修正量の計算が行われる。モデル修正量の計算を行う際には、波形の残差を震源として逆伝播させて内積処理を行う手法や、データを時間窓で切り出して、残差の一部を用いてモデル変動を計算する手法など、様々な手法がある。
そして、モデル修正部２６は、モデル記憶部２２が記憶したモデルに対して、算出した修正量の修正を施す。

モデル修正部２６で得られた修正量のデータは、収束判断部２７に供給される。収束判断部２７は、モデル記憶部２２に記憶されたモデルと、モデル修正部２６で得られた修正量とを比較して、修正状態が収束したか否かの判断を行う。収束判断部２７では、修正状態が収束したと判断したとき、順方向伝播波動場計算部２３と逆方向伝播波動場計算部２４での波動場計算を終了し、出力部２８に対して、モデル記憶部２２に記憶されたモデルの出力を指示する。

出力部２８では、モデル修正部２６からの指示に基づいて、モデル記憶部２２が記憶した最後に更新されたモデルを出力装置３０に出力する。出力装置３０は、出力部２８から供給されたモデルを、地下構造ｘのモデルとして出力する。

［３．演算処理装置での演算処理］
図３のフローチャートは、演算処理装置２０での演算処理の流れを示す。
まず、初期モデル設定部１３は、震源データ及び地震波形データと、地下構造ｘを推測したデータとを持つ初期モデルを設定する（ステップＳ１１）。初期モデル設定部１３が設定した初期モデルは、モデル記憶部２２に記憶される。
そして、モデル記憶部２２に記憶される初期モデルを使って、順方向伝播波動場計算部２３が順方向の伝播波動場を計算し、逆方向伝播波動場計算部２４が逆方向の伝播波動場を計算する。これらの順方向及び逆方向の伝播波動場の内積計算部２５での内積計算結果から、モデル修正部２６が修正量を算出する（ステップＳ１２）。

モデル修正部２６は、入力した初期モデルと修正量を結合して修正モデルを計算する（ステップＳ１３）。モデル修正部２６では、この修正されたモデルの計算結果と観測波形データとの残差、および推定モデル自体への拘束条件などから算出される統計的指標と、修正前のモデルから算出される統計的指標が比較され、その比較結果から修正状態が収束したか否かが判断される（ステップＳ１４）。ここで、収束していないと判断した場合には、ステップＳ１２の処理に戻り、修正したモデルでの順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４での伝播波動場の計算が繰り返し行われる。なお、ここで述べる統計的手法とは、モデル推定の良さを統計的に評価するための、最適化問題で用いられる指標のことである。たとえば、残差二乗和や拘束条件付き重み付き残差二乗和など、最適化問題には様々な指標が存在し、本アルゴリズムでは、どの形式の指標を用いてもよい。

そして、ステップＳ１４の判断で、収束したと判断したとき、モデル修正部２６は、最後に修正されてモデル記憶部２２が記憶した中間モデルを、最終的な推定モデルに確定する（ステップＳ１５）。出力装置３０は、この最終的な推定モデルを、地下構造ｘの推測結果として出力する。

［４．モデル修正量の計算処理］
図４のフローチャートは、演算処理装置２０でのモデル修正量の計算処理の流れを示す。
順方向伝播波動場計算部２３は、震源の地震波形データと初期モデル（又は中間モデル）の地下構造ｘとを使用して、順方向の伝播波動場を計算する（ステップＳ２１，Ｓ２２）。ここでは、離散化した運動方程式を使用して、過去から未来の順方向への地震波形の伝わり状態が計算される。
また、逆方向伝播波動場計算部２４は、各観測点の地震波形データと初期モデル（又は中間モデル）の地下構造ｘとを使用して、逆方向の伝播波動場を計算する（ステップＳ２３，Ｓ２４）。ここでは、離散化した運動方程式を使用して、実際の地震波形の伝わり状態とは逆の未来から過去の逆方向への地震波形の伝わり状態が計算される。

そして、内積計算部２５は、順方向の伝播波動場と逆方向の伝播波動場との内積を計算する（ステップＳ２５）。この内積の計算結果から、モデル修正部２６では、初期モデル又は中間モデルの修正量が算出される（ステップＳ２６）。

［５．メッシュの説明］
順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４で伝播波動場を計算する際には、地震波が伝わる地下を、一定の距離ごとに区切ったメッシュ集合として、それぞれのメッシュごとの波動場を計算する。順方向、逆方向それぞれの波動場を適切に内積計算することによって、メッシュの構造（密度や弾性定数など）に修正があった時に、観測点Pで観測される波形に生じる変動を計算することができる。
すなわち、図５に示すように、順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４は、震源ｋと観測点Ｐの間を、一定の体積のメッシュｍ_１，ｍ_２，・・・ｍ_ｉ（ｉは任意の整数）に分割する。
ここで、観測点Ｐでは、震源ｋにより生じた波動場が観測波形Ｗｐとして観測され、また初期モデルに対して、初期波形Ｗｋが計算される。
このとき、観測点Ｐでの観測波形Ｗｐは、初期波形Ｗｋに各メッシュｍ_１〜ｍ_ｉの変更が観測点Ｐでの観測波形に与える変動Ｗm_１〜Ｗm_ｉの成分を加算したものになる。すなわち、各メッシュに与える変化量が十分に小さいとすれば、摂動論にしたがって、全体の変化量は、各部分の変動の影響の和で示すことができる。初期モデルとして適切なモデルを設定して、図３のフローチャートで説明した修正を繰り返すことで、各メッシュの地下構造は正確なものに収束すると期待される。

順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４では、それぞれのメッシュｍ_１〜ｍ_ｉごとの変動が計算される。ここで、順方向伝播波動場計算部２３は、震源ｋに近い側のメッシュから順に伝播波動場の伝わり状態を計算し、逆方向伝播波動場計算部２４は、観測点Ｐに近い側のメッシュから順に逆方向に伝播波動場の伝わり状態を計算する。

図６は、それぞれのメッシュで内積が行われる状態を示す。
内積計算部２５は、ｋ番目の震源から順伝播させた波動場ｃと、ｐ番目の観測点から逆伝播させた波動場ｚとについて、ｌ番目のメッシュｍ_ｌの地点で、内積計算を行う。
次式は、内積計算の例を示す。

［数１］
ｃ^（ｋ）＊（ω^２Ｔ^（ｌ）−Ｈ^（ｌ））Ｚ_ｉ ^（ｐ）

この［数１］式においては内積計算の概念を示すために、時間領域ではなく、周波数領域で記した。ωは周波数であり、Ｔは密度に関する部分を離散化した行列（質量行列）であり、Ｈは弾性定数に関する部分を離散化した行列（剛性行列）である。上付きの(l)は１番目のメッシュに対する変動に関する部分を取り出したものであることを意味する。
このように、地下構造の各地点のメッシュごとに内積計算が行われる。
この［数１］式は、後述する運動方程式（［数２］式）の一次摂動を計算することで導出される。

ここで、本例の場合には、順伝播と逆伝播の相反性が保障されるため、内積演算により全ての情報を取り出すことができる。
図７は、順伝播と逆伝播とが相反性を満たさない場合（図７Ａ）と、順伝播と逆伝播とが相反性を満たす場合（図７Ｂ）とを比較した模式図である。
順伝播と逆伝播とが相反性を満たさない場合には、内積の結果から取り出される情報が減り、推定における誤差源となる。一方、順伝播と逆伝播とが相反定理を満たす場合には、地下構造推測のための情報を、内積演算で全て取り出すことができるようになる。

［６．具体的な計算手法］
次に、モデル修正量を計算する具体的な例について説明する。
本例での計算処理の手法は、時間領域差分法と呼ばれる時間領域の既存手法を改良して、順伝播と逆伝播の相反性を保ちつつ、高精度でかつ効率のよい計算が行えるようにしたものである。特に、順伝播と逆伝播の相反性が保たれ、時間反転が可能であるという点が重要である。本例の場合には、安定性は自動的に保たれる。

順方向伝播波動場計算部２３及び逆方向伝播波動場計算部２４が波動場を計算する際には、地球内部を弾性体として仮定し、弾性体の運動方程式が設定される。この運動方程式では、地下の密度や弾性定数などをパラメータとし、地震の揺れの状態がメッシュ状に離散化して計算される。
波動場計算手法において、時刻方向の離散化手法には、周波数領域の手法と、時間領域の手法とがある。本例では、時間領域の手法を採用して波動場を計算する。
順方向の伝播波動場を計算する場合と、逆方向の伝播波動場を計算する場合とでは、時間刻みの正負を入れ替えて計算する。

ここで、順方向の伝播波動場と、逆方向の伝播波動場とは、相反定理を満たす必要がある。すなわち、震源に置いた点震源が作る波動場に対する観測点に置いた点震源による作用と、観測点に置いた点震源が作る波動場に対する震源に置いた点震源による作用とが等しいことを示す。これは、簡単に言うと、震源と観測点の関係を入れ替えても、同じ現象を記述できることである。

相反定理の条件を満たすために、運動方程式中に含まれる、最終的に計算に用いる演算子を対称（もしくはエルミート）行列の形とし、離散空間で相反定理が成立するように設定する。このように離散空間で相反定理が成立する運動方程式とすることで、自動的に時間反転対称となる。
この条件を満たす運動方程式の一例を示すと、ここでは、次の［数２］式の運動方程式を採用する。

［数２］

この運動方程式において、Δｔは時間刻みの幅、ｃ_＋，ｃ_０，ｃ₋はそれぞれ未来，現在，過去の波動場を離散かした量をベクトルで示すものである。また、Ｂは弾性体運動方程式の密度に関する項の逆演算行列、Ｈは弾性定数に関する部分を離散化した行列（剛性行列）であり、ｈは震源における外力項である。
この［数２］式において、ＢとＨを、それぞれ正定値，半正定値の対称行列もしくはエルミート行列として定義することで、予測子修正子法を用いて計算する際にも、自動的に相反定理及び時間反転対称性が成り立つ演算子が生成される。

ここで、［数２］式が相反定理を満たすことについて説明する。
［数２］式は、一次の反復法を用いて、以下のように展開することができる。

［数３］

このような手順で、未来のステップの解ｃ_＋を陽的解法で求めることができる。この手順は次に示すようい、１ステップにまとめることができる。

［数４］

ここで、両辺にＢの逆行列（Ｂが対称なので対称行列となる）を形式的にかけると、１ステップで計算する陽的解法は以下のようになる。

［数５］

この［数５］式において、下線を引いた行列が対称もしくはエルミートなことは共役転置をとれば明らかである。

［７．本例の計算手法による効果］
以上説明した本例の計算手法により地下構造を推定することで、非常に精度の高い地下構造の推定が可能になる。例えば、地下に存在する含油層の厚さを推定する際の精度として、従来よりも約１０％精度が向上した。このため、本例によると、石油やシェールオイルなどの含油層や断層などの地下構造の推定が、従来よりも高精度に行えるようになる。

［８．変形例］
なお、上述した実施の形態では、地下構造の推定を行う手法について説明したが、同じ計算手法を震源情報の推定のための波形インバージョンに利用することもできる。また、上述の実施形態の例では単一の地震の例を図示したが、複数の地震の震源情報及び複数の観測点で観測された波形データを入力データに含めることができる。
また、上述した実施の形態では、相反定理及び時間反転対称性が成り立つ演算子として使用した運動方程式として、［数２］式とし、時間領域差分法の反復法を用いて計算した。これに対して、密度に関する項の逆演算行列と、弾性定数に関する部分を離散化した行列とを、正定値又は半正定値の対称行列もしくはエルミート行列として定義して、相反定理及び時間反転対称性が成り立つように定義した運動方程式であれば、その他の離散化手法を使用してもよい。

また、図１及び図２で説明した演算処理の構成や、図３及び図４で説明した演算処理の手順は、好適な一例を示したものであり、図示の構成や手順に限定されるものではない。

１１…震源データ収集部、１２…地震波形データ収集部、１３…初期モデル作成部、１４…入力データ結合部、２０…演算処理装置、２１…データ入力部、２２…モデル記憶部、２３…順方向伝播波動場計算部、２４…逆方向伝播波動場計算部、２５…内積計算部、２６…モデル修正部、２７…収束判断部、２８…出力部、３０…出力装置

Claims

地球内部を伝播し、観測された地震波形データを解析することで、地球内部の地下空間の構造を探査する地下構造探査システムであって、
地震波の震源と推定される地下構造の初期モデルと観測点との情報とを含む入力データ又はその初期モデルを修正した中間モデルが記憶される記憶部と、
前記記憶部に記憶された初期モデル又は中間モデルに対して、離散化された時間領域の運動方程式を使用して、過去から未来に向けた順方向の伝播波動場を演算すると共に、未来から過去に向けた逆方向の伝播波動場を演算する伝播波動場計算部と、
前記伝播波動場計算部で演算された順方向の伝播波動場と逆方向の伝播波動場を利用して初期モデル又は中間モデルの修正量を算出し、算出された修正量の修正を前記記憶部に記憶された初期モデル又は中間モデルに対して施すモデル修正部と、
前記モデル修正部での修正状態の収束を判断する収束判断部と、
前記収束判断部で収束を判断した場合に、前記モデル修正部で修正されて前記記憶部に記憶された中間モデルによる地下構造を、探査対象の地下構造として読み出す出力部とを備え、
前記伝播波動場計算部での演算に使用する前記運動方程式は、時間領域の反復法を用いて行う計算手法であって、密度に関する項の逆演算行列と、弾性定数に関する部分を離散化した行列とを、それぞれ正定値、半正定値の対称行列もしくはエルミート行列として定義して、相反定理及び時間反転対称性が成り立つ演算子として使用した運動方程式とした
地下構造探査システム。
前記運動方程式は、

（ここで、Ｃ_＋，Ｃ_０，Ｃ₋は未来，現在，過去の波動場を離散化した量のベクトル、Ｂは密度に関する項の逆演算行列、Ｈは弾性定数に関する部分を離散化した行列、ｈは震源（もしくは観測点）における外力項、Δｔは時間刻みの幅）とした
請求項１に記載の地下構造探査システム。
地球内部を伝播する地震波を解析することで、地球内部の地下空間の構造を探査する地下構造探査方法であって、
地震波の震源と推定される地下構造の初期モデルと観測点との情報とを含む入力データ又はその初期モデルを修正した中間モデルを設定して記憶部に記憶されるモデル記憶ステップと、
前記モデル記憶ステップで前記記憶部に記憶された初期モデル又は中間モデルに対して、伝播波動場計算部により、離散化された時間領域の運動方程式を使用して、過去から未来に向けた順方向の伝播波動場を演算すると共に、未来から過去に向けた逆方向の伝播波動場を演算する伝播波動場演算ステップと、
モデル修正部により、前記伝播波動場演算ステップで演算された順方向の伝播波動場と逆方向の伝播波動場とを内積計算して、内積計算した結果から初期モデル又は中間モデルの修正量を算出し、算出された修正量の修正を前記モデル記憶ステップにより記憶された初期モデル又は中間モデルに対して施すモデル修正ステップと、
前記モデル修正ステップでの修正状態の収束を、集束判断部が判断する収束判断ステップと、
前記収束判断ステップで収束したと判断された場合に、前記モデル修正ステップで修正されて前記記憶部に記憶された中間モデルによる地下構造を、探査対象の地下構造として確定して出力部が出力する出力ステップと、を含む地下構造探査方法であり、
前記伝播波動場計算部による前記伝播波動場演算ステップでの演算に使用する前記運動方程式は、密度に関する項の逆演算行列と、弾性定数に関する部分を離散化した行列とを、正定値又は半正定値の対称行列もしくはエルミート行列として定義して、相反定理及び時間反転対称性が成り立つ演算子として使用した運動方程式とした
地下構造探査方法。
前記運動方程式は、

（ここで、Ｃ_＋，Ｃ_０，Ｃ₋は未来，現在，過去の波動場を離散化した量のベクトル、Ｂは密度に関する項の逆演算行列、Ｈは弾性定数に関する部分を離散化した行列、ｈは震源（もしくは観測点）における外力項、Δｔは時間刻みの幅）とした
請求項３に記載の地下構造探査方法。