JP6778628B2 - 地下構造推定装置 - Google Patents

Description

本発明は、地下構造推定装置に関する。

従来、地盤・岩盤における透水係数を推定する技術が知られている（例えば、特許文献１を参照）。

特開２００１‐３２４１９号公報

上記特許文献１では、試験用ボーリング孔からの流量計測を必須とせずに、間隙水圧のみの計測から地盤・岩盤における透水係数を推定する。しかし、上記特許文献１の技術では、間隙水圧の計測は必要であり、かつ地下構造としては透水係数のみしか推定することができない。

本発明は上記事実を考慮して、地震の観測データを用いて地下構造を推定することができる地下構造推定装置を提供することを目的とする。

上記目的を達成するために、本発明の地下構造推定装置は、地下構造モデルと、前記地下構造モデルに対応する地点における各時刻の地震の観測データの時系列とを入力とし、前記地下構造モデルのパラメータ及び地震予測モデルから予測される前記時刻の地震の状態から得られる地震データと、前記時刻における地震の観測データとの差分に応じて、前記時刻の地震の状態及び前記地下構造モデルのパラメータを修正することを繰り返す推定部を含んで構成されている。これにより、地震の観測データを用いて地下構造を推定することができる。

本発明の前記地下構造モデルのパラメータは、前記地下構造モデルの地表の各格子点に対する、前記格子点の地下方向に位置する各層に関する前記パラメータを含むようにすることができる。これにより、地下構造を推定する際に計算量を低減させることができる。

本発明の各格子点の前記層に関するパラメータは、前記格子点間で、前記格子点の位置に応じた相関が予め設定されるようにすることができる。これにより、地下構造の連続性を考慮して、地下構造を精度よく推定することができる。

本発明の前記層に関するパラメータは前記層の層厚を含み、前記層厚は、前記格子点間で、前記格子点の位置に応じた相関が予め設定されるようにすることができる。これにより、地下構造の層厚の連続性を考慮して、地下構造を精度よく推定することができる。

本発明では、前記地下構造モデルの前記格子点間の各々について、前記格子点間の距離が近いほど、一方の格子点の地下方向に位置する各層に関する前記パラメータと、他方の格子点の地下方向に位置する各層に関する前記パラメータとが対応するように予め設定されるようにすることができる。これにより、地下構造の格子点間の距離を考慮して、地下構造を精度よく推定することができる。

本発明の地下構造推定装置は、前記地下構造モデルの前記パラメータを表すパラメータベクトルと前記地震による各格子点の速度及び応力を含む状態ベクトルとを有する拡大状態ベクトルと、前記観測データを表す各時刻の前記地震の地震動データとを入力とし、前記推定部は、時刻ｔ−１の前記拡大状態ベクトルと前記地震予測モデルとに基づいて、時刻ｔの前記拡大状態ベクトルを予測する予測部と、前記予測部によって予測された時刻ｔの前記拡大状態ベクトルから生成される前記地震データと、時刻ｔの前記地震の地震動データとに基づいて、時刻ｔの前記拡大状態ベクトルを修正する修正部と、前記観測データの終端時刻Ｔまで、前記予測部による予測及び前記修正部による修正を繰り返し、各時刻において得られた前記修正部による推定結果に含まれる、前記地下構造モデルの前記パラメータを平滑化することにより、前記拡大状態ベクトルに含まれる前記地下構造モデルの前記パラメータを取得する平滑化部と、を含むようにすることができる。これにより、逐次フィルタリング手法を用いて地下構造を精度よく推定することができる。

本発明によれば、地下構造モデルに対応する地点における地震の観測データの時系列を入力とし、地下構造モデルのパラメータ及び地震予測モデルから予測される地震データと、地震の観測データとの差分に応じて、地震の状態及び地下構造モデルのパラメータを修正することを繰り返すことにより、地震の観測データを用いて地下構造を推定することができる、という効果が得られる。

本発明の実施形態に係る地下構造推定装置の概略構成を示すブロック図である。 地下構造モデルを説明するための説明図である。 本発明の実施形態に係る地下構造推定装置の推定部の詳細な構成の一例を示す図である。 逐次データ同化において用いられる分散共分散行列に対して相関が設定された場合の一例を示す図である。 地下構造モデルの層に対して相関を与えない場合のアンサンブルカルマンフィルタの各アンサンブルの一例と、地下構造モデルの層に対して相関を与えた場合のアンサンブルカルマンフィルタの各アンサンブルの一例とを示す図である。 地下構造推定処理ルーチンの一例を示す図である。 データ同化処理ルーチンの一例を示す図である。 シミュレーション実験を説明するための説明図である。 シミュレーション実験の実験結果を示す図である。

以下、本発明の実施形態について詳細に説明する。

＜地下構造推定装置のシステム構成＞

図１は、本発明の実施形態に係る地下構造推定装置の構成の一例を示すブロック図である。地下構造推定装置１０は、機能的には、図１に示されるように、受付部１２、コンピュータ１４、及び表示部１６を含んだ構成で表すことができる。

受付部１２は、地下構造推定装置の外部から入力された情報を受け付ける。例えば、受付部１２は、新たに地震が発生した場合、地震の観測データの一例である震源パラメータ及び各観測地点における各時刻の地震の速度データを受け付ける。震源パラメータは、気象庁等からの地震情報に基づいて入力される。そして、受付部１２は、震源パラメータ及び各観測地点における各時刻の地震の速度データを、後述する観測データ記憶部１８へ格納する。なお、地震の速度データは、地震の地震動データの一例である。

コンピュータ１４は、ＣＰＵ（Central Processing Unit）、各処理ルーチンを実現するためのプログラム等を記憶したＲＯＭ（Read Only Memory）、データを一時的に記憶するＲＡＭ（Random Access Memory）、記憶手段としてのメモリ、ネットワークインタフェース等を含んで構成されている。

コンピュータ１４は、機能的には、図１に示されるように、観測データ記憶部１８と、モデル記憶部２０と、情報取得部２２と、推定部２４とを含んで構成される。

観測データ記憶部１８には、地震毎に、震源パラメータ及び各観測地点における各時刻の地震の速度データが格納される。新たな地震が発生する毎に、当該地震の観測データが観測データ記憶部１８へ格納される。

モデル記憶部２０には、地下構造を表す地下構造モデルと、地下構造モデルのパラメータを表すパラメータベクトルとが格納される。

本実施形態では、データ同化により地下構造モデルのパラメータを推定する。そのため、本実施形態では、観測データ記憶部１８に格納された地震の速度データと、モデル記憶部２０に格納された地下構造モデルのパラメータベクトル及び地震予測モデルから予測される地震データとの差分を修正するように、地下構造モデルのパラメータを推定する。

本実施の形態で用いる地下構造モデルとパラメータとについて、以下説明する。

地震動予測において利用されている有限差分法を用いたシミュレーションでは、地下構造モデルの各格子点に物性値（例えば、伝播速度、密度、及び減衰定数等）を表すパラメータを付与することにより、地下構造を表現する。

例えば、図２の地下構造モデル２５のように、地下構造モデル２５の各格子点（ｉ，ｊ，ｋ）には、密度ρ、Ｓ波のせん断波速度Ｖ_Ｓ、Ｐ波のせん断波速度Ｖ_Ｐ、Ｓ波の減衰定数Ｑ_Ｓ、及びＰ波の減衰定数Ｑ_Ｐの物性値がパラメータとして付与される。

しかし、地下構造モデルの各格子点に対してパラメータを付与し、データ同化により各パラメータを推定する場合、推定対象のパラメータが多くなりすぎてしまい、計算量が増加してしまう。また、パラメータの増加により、パラメータの推定精度が低下する可能性がある。

そこで、本実施形態では、データ同化に用いる地下構造モデルにおいて、地表の各格子点の地下方向に位置する各層に対してパラメータを付与する。具体的には、図２の地下構造モデル２６のように、地表の格子点毎に、地下構造モデル２６の各層に対する、層の層厚Ｈをパラメータとして付与する。また、各層に対して物性値（ρ，Ｖ_Ｓ，Ｖ_Ｐ，Ｑ_Ｓ，Ｑ_Ｐ）をパラメータとして付与する。

ここで、ｘ方向の格子の区切り数をＩ、ｙ方向の格子の区切り数をＪ、ｚ方向の格子の区切り数をＫとすると、上記図２の地下構造モデル２５の場合、総パラメータ数は、以下の式（１）のようになる。

総パラメータ数＝１格子点当たりのパラメータ数×格子点数
＝５×Ｉ×Ｊ×Ｋ ［個］ （１）

一方、上記図２に示される本実施形態の地下構造モデル２６の場合、層数をＬとすると、総パラメータ数は、以下の式（２）のようになる。

総パラメータ数
＝１層当たりのパラメータ数×層数＋１格子点当たりの層数×地表の格子点数
＝５×Ｌ＋Ｉ×Ｊ×（Ｌ−１） ［個］ （２）

地下構造モデル２５では、上記式（１）から分かるように、格子点毎に物性値のパラメータ（ρ，Ｖ_Ｓ，Ｖ_Ｐ，Ｑ_Ｓ，Ｑ_Ｐ）が付与されると、推定対象のパラメータ数が膨大な数となる。一方、本実施形態の地下構造モデル２６では、上記式（２）から分かるように、物性値のパラメータ（ρ，Ｖ_Ｓ，Ｖ_Ｐ，Ｑ_Ｓ，Ｑ_Ｐ）は層毎に付与され、層厚のパラメータＨは地表の格子点毎に付与される。そのため、本実施形態の地下構造モデル２６では、推定対象のパラメータ数が減少し、データ同化によるパラメータ推定の際の計算量が抑制される。

情報取得部２２は、観測データ記憶部１８に格納された観測データを取得する。そして、情報取得部２２は、観測データの震源パラメータに含まれる震源位置等の情報に基づき、当該観測データの地震の速度データを、地下構造モデルのパラメータの推定に用いるか否かを判定する。

例えば、情報取得部２２は、震源パラメータに含まれる震源位置と地下構造モデルに対応する地点との間の距離が、予め設定した閾値以上である場合には、地下構造モデルの領域内ではないと判定し、地下構造モデルのパラメータの推定に用いないと判定する。

一方、情報取得部２２は、震源パラメータに含まれる震源位置と地下構造モデルに対応する地点との間の距離が、予め設定した閾値未満である場合には、地下構造モデルの領域内であると判定し、地下構造モデルのパラメータの推定に用いると判定する。地下構造モデルに対応する地点と関係ない領域の地震の速度データについては、地下構造のパラメータを推定する際には適切でないため除外される。

また、情報取得部２２は、観測データである地震の速度データのＳＮ比が所定の条件を満たすか否かに応じて、当該観測データの地震の速度データを地下構造モデルのパラメータの推定に用いるか否かを判定する。ＳＮ比が所定の条件を満たすか否かは、地震の速度データから計算されるフーリエ振幅スペクトルの形状に基づき判定される。

情報取得部２２は、観測データである地震の速度データの低周波成分がノイズレベルを下回る場合には、地下構造のパラメータを推定する際には適切でないと判定し、当該地震の速度データを地下構造モデルのパラメータの推定に用いないことを判定する。一方、情報取得部２２は、地震の速度データの低周波成分がノイズレベルを上回る場合には、当該地震の速度データを地下構造モデルのパラメータの推定に用いることを判定する。

また、情報取得部２２は、モデル記憶部２０に格納された地下構造モデルのパラメータを表すパラメータベクトルを取得する。情報取得部２２により取得される地下構造モデルのパラメータベクトルは、前回までの地震の速度データに基づいてデータ同化により得られたものである。

推定部２４は、情報取得部２２によって取得された地下構造モデルのパラメータベクトル及び地震予測モデルから予測される地震データと、情報取得部２２によって取得された地震の速度データとの差分に応じて、地下構造モデルのパラメータベクトルを修正することを繰り返す。

具体的には、地震予測モデルの一例である有限差分法によるシミュレーションによって予測される地震の状態から得られる地震データと、観測データである地震の速度データとの差分に応じて、データ同化により地下構造モデルのパラメータを推定する。

本実施形態では、データ同化の一例として、逐次データ手法であるアンサンブルカルマンフィルタを用いる場合を例に説明する。逐次データ手法では、時系列データが得られる毎に同化が行われる。

データ同化の処理の流れを以下に示す。本実施形態では、データ同化に用いるシステムモデルは以下の式（３）及び式（４）によって表される。

（３）

（４）

上記式（３）及び式（４）において、ｔは時刻を表し、ｘ_ｔは時刻ｔにおける拡大状態ベクトルを表し、ｙ_ｔは時刻ｔにおける観測データを表す。ｆ_ｔはシステムモデルであるシミュレーションを表す。また、観測ノイズｗ_ｔは正規分布Ｎ（０，Ｒ_ｔ）に従い、システムノイズｖ_ｔは任意の分布に従う。Ｈ_ｔは観測行列である。

上記式（３）及び式（４）を条件付き分布で表現すると、以下の式（５）〜式（７）に示されるような一般状態空間モデルとなる。

（５）

（６）

（７）

ここで、ｘ_０は地震の観測データと同化する前の拡大状態ベクトルを表す。また、地下構造モデルのパラメータベクトルθについては、既往のモデルや構造探査結果や地質情報に応じて設定される。

次に、参考文献（長尾大道、「データ同化の基礎理論とその応用」、自動チューニング研究会、第６回自動チューニング技術の現状と応用に関するシンポジウム（ＡＴＴＡ２０１４）スライド資料）を参照して、アンサンブルカルマンフィルタについて説明する。

［１］一期先予測
アンサンブルカルマンフィルタの一期先予測において、拡大状態ベクトルｘ_ｔに対応する各アンサンブルｘ^（ｉ） _{ｔ−１｜ｔ−１}をシミュレーションｆ_ｔによって時間発展させ、一期先予測分布のアンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}を得る。

時刻ｔ−１におけるフィルタ分布は、以下の式（９）によって近似表現される。なお、δ（・）はデルタ関数を表す。また、各アンサンブルｘ^（ｉ） _{ｔ−１｜ｔ−１}を、システムモデルであるシミュレーションｆ_ｔによって時間発展させた各アンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}は、以下の式（１０）によって表される。また、式（１０）の各アンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}によって、一期先予測分布は以下の式（１１）によって近似表現される。

（８）

（９）

（１０）

（１１）

［２］フィルタリング
次に、アンサンブルカルマンフィルタのフィルタリングにおいて、一期先予測分布のアンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}からフィルタ分布のアンサンブルｘ^（ｉ） _ｔ｜ｔを得る。

一期先予測分布のアンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}の分散共分散行列Ｖ＾_{ｔ｜ｔ−１}は、以下の式（１２）によって表される。なお、式中の「’」は行列の転置を表す。また、観測ノイズｗ^（ｉ） _ｔの分散共分散行列Ｒ＾_ｔは、以下の式（１３）によって表される。

また、カルマンゲインＫ＾_ｔは、以下の式（１４）によって表される。そして、フィルタ分布のアンサンブルｘ^（ｉ） _ｔ｜ｔは、以下の式（１５）によって表される。なお、ｙ_ｔは時刻ｔにおける観測データを表す。

（１２）

（１３）

（１４）

（１５）

［３］平滑化
アンサンブルカルマンフィルタの平滑化において、一期先予測及びフィルタリングの繰り返しにより得られたアンサンブルを平滑化する。具体的には、時刻ステップｔ＝Ｔまでの観測データから得られた各アンサンブルに基づいて、平滑化分布ｐ（ｘ_０｜ｙ_１：Ｔ）を計算することにより、拡大状態ベクトルを推定する。

本実施の形態におけるアンサンブルカルマンフィルタを用いた各処理について、以下、具体的に説明する。

本実施の形態における拡大状態ベクトルｘ_ｔは、シミュレーションによって得られる地下構造モデルの各格子点の速度ｖ_{ｉ，ｊ，ｋ}及び応力τ_{ｉ，ｊ，ｋ}を要素とする状態ベクトルｘ^〜 _ｔと、地下構造モデルの各層ｌの物性値を表すパラメータ（ρ_ｌ，Ｖ_Ｓ，ｌ，Ｖ_Ｐ，ｌ，Ｑ_Ｓ，ｌ，Ｑ_Ｐ，ｌ）及び各層ｌの層厚Ｈ_{ｉ，ｊ，ｌ}を要素とするパラメータベクトルθとを含む。

また、本実施形態の状態ベクトルｘ^〜 _ｔは以下の式（１６）によって表される。状態ベクトルｘ^〜 _ｔのｑ_{ｉ，ｊ，ｋ}は、シミュレーションによって得られる各格子点の速度を表す。また、状態ベクトルｘ^〜 _ｔのτ_{ｉ，ｊ，ｋ}は、シミュレーションによって得られる各格子点の応力を表す。状態ベクトルｘ^〜 _ｔは地震の状態の一例である。

本実施形態のパラメータベクトルθは以下の式（１７）によって表される。パラメータベクトルθには、各層ｌの物性値のパラメータ（ρ_ｌ，Ｖ_Ｓ，ｌ，Ｖ_Ｐ，ｌ，Ｑ_Ｓ，ｌ，Ｑ_Ｐ，ｌ）と、層厚のパラメータＨ_{ｉ，ｊ，ｌ}とが含まれる。

また、地震データの一例である本実施形態の拡大状態ベクトルｘ_ｔは、以下の式（１８）によって表される。また、Ｍ個の観測点についての時刻ｔの地震の速度データｙ_ｔは、以下の式（１９）によって表される。地震の速度データｙ_ｔの各要素ａ_ｍは、地点ｍにおいて観測された地震の速度データを表す。

推定部２４は、図３に示されるように、予測部２４０と、修正部２４２と、平滑化部２４４とを含んで構成されている。

予測部２４０は、情報取得部２２によって取得された地下構造モデルのパラメータベクトルθを取得する。また、予測部２４０は、情報取得部２２によって取得された地震の速度データｙ_ｔを取得する。

次に、予測部２４０は、各値の初期値を設定する。具体的には、予測部２４０は、上記式（６）の条件付き確率分布に従って、拡大状態ベクトルの初期値ｘ_０を設定する。また、予測部２４０は、拡大状態ベクトルの初期値ｘ_０を設定する際、パラメータベクトルθのうちの層厚Ｈについて、格子点間で、格子点の位置に応じた相関を設定する。具体的には、予測部２４０は、地下構造モデルの格子点間の各々について、格子点間の距離が近いほど、一方の格子点の地下方向に位置する各層の層厚と、他方の格子点の地下方向に位置する各層の層厚とが近い値になるように設定する。より詳細には、上記式（１２）の分散共分散行列Ｖ＾_{ｔ｜ｔ−１}の層厚Ｈに関する項について、格子点の位置に応じた相関を設定することにより、層厚の相関を設定する。

図４に、層厚Ｈに関する項についての分散共分散行列Ｖ＾_{ｔ｜ｔ−１}の一例を示す。図４（Ａ）の横軸α及び縦軸βは行列の各要素の位置を表し、明るい領域ほど相関値が高く暗い領域ほど相関値が低い。図４（Ｂ）は、図４（Ａ）を３次元で表したものであり、水平面における横軸α及び縦軸βは行列の各要素の位置を表し、垂直方向は相関値を表す。

図４（Ａ）及び（Ｂ）に示されるように、本実施形態では、各格子点について、格子点間の距離が小さいほど相関値が高くなるように、かつ格子点間の距離が大きいほど相関値が低くなるように、層厚Ｈに関しての分散共分散行列Ｖ＾_{ｔ｜ｔ−１}を設定する。そして、相関が設定された分散共分散行列Ｖ＾_{ｔ｜ｔ−１}に応じて各アンサンブルが生成される。これにより、地下構造の連続性が反映される。

地下構造は、断層が存在しない限り層の形状は滑らかに変化する。本実施形態では、この点を考慮し、同一の層に対応する層厚に相関を与える。

図５に、相関を設定しない場合のアンサンブル１〜５の一例と相関を設定した場合のアンサンブル１〜５の一例とを示す。なお、図中のＸgridは上記図２のｘ方向に対応し、Ｚgridは上記図２のｚ方向に対応する。図５（Ａ）に示されるように、相関が設定されない場合のアンサンブル１〜５においては、Ｚgridの値はＸgridの位置に関係なくランダムな値をとる。しかし、地下構造における層は、水平方向において連続していると考えられる。そのため、各アンサンブルのうちの層厚がランダムな値をとりながらデータ同化が行われると、地下構造モデルのパラメータを精度よく推定することができない。

一方、図５（Ｂ）に示されるように、相関が設定された場合のアンサンブル１〜５においては、Ｚgridの値はＸgridの位置に応じた値となり、例えばＸgridが隣接する点同士では近い値をとる。そのため、層厚Ｈに関する項についての分散共分散行列Ｖ＾_{ｔ｜ｔ−１}に対して相関を設定することにより、地下構造の層の連続性が反映される。

次に、予測部２４０は、時刻ｔ−１の拡大状態ベクトルを表す各アンサンブルｘ^（ｉ） _{ｔ−１｜ｔ−１}と地震予測モデルの一例である有限差分法のシミュレーションｆ_ｔとに基づいて、上記式（１０）に従って、時刻ｔの拡大状態ベクトルを表す各アンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}を予測する。

修正部２４２は、予測部２４０によって予測された時刻ｔの拡大状態ベクトルを表す各アンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}から生成される地震データと、時刻ｔの地震の速度データｙ_ｔとに基づいて、上記式（１５）に従って、時刻ｔの拡大状態ベクトルを表す各アンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}を修正する。そして、修正部２４２は、拡大状態ベクトルを表す各アンサンブルｘ^（ｉ） _ｔ｜ｔを得る。

そして、地震の速度データの終端時刻Ｔまで、予測部２４０による一期先予測及び修正部２４２によるフィルタリングが繰り返される。

平滑化部２４４は、各時刻において得られた修正部２４２による推定結果に含まれる、地下構造モデルのパラメータベクトルθを平滑化することにより、地下構造モデルのパラメータベクトルθを取得する。パラメータベクトルθには、各層の物性値のパラメータ（ρ，Ｖ_Ｓ，Ｖ_Ｐ，Ｑ_Ｓ，Ｑ_Ｐ）と、各層の層厚のパラメータＨとが含まれているため、データ同化により地下構造のパラメータが得られたことになる。

そして、平滑化部２４４は、平滑化によって得られた地下構造モデルのパラメータベクトルを、モデル記憶部２０へ格納する。モデル記憶部２０に格納されたパラメータベクトルは、次回の地震の観測データが得られたときのデータ同化の際に用いられる。

表示部１６は、モデル記憶部２０に記憶されたパラメータベクトル等の推定結果を表示する。地下構造推定装置１０を操作するユーザは、表示部１６によって表示された結果を確認する。

＜地下構造推定装置の作用＞
次に、図６を参照して、地下構造推定装置１０の作用を説明する。地下構造推定装置１０の受付部１２は、新たに地震が発生した場合、新たな地震の震源パラメータ及び各観測地点における各時刻の地震の速度データを受け付け、観測データ記憶部１８へ格納する。そして、地下構造推定装置１０は、図６に示す地下構造推定処理ルーチンを実行する。地下構造推定処理ルーチンは、新たな地震の観測データが観測データ記憶部１８へ格納される毎に実行される。

ステップＳ１００において、情報取得部２２は、観測データ記憶部１８に格納された観測データの震源パラメータを取得する。

ステップＳ１０２において、情報取得部２２は、震源パラメータに含まれる震源位置の情報に基づいて、観測データが地下構造モデルの領域内であるか否かを判定する。観測データが地下構造モデルの領域内である場合には、ステップＳ１０４へ進む。一方、観測データが地下構造モデルの領域内でない場合には、地下構造推定処理ルーチンを終了する。

ステップＳ１０４において、情報取得部２２は、観測データ記憶部１８に格納された観測データの地震の速度データを取得する。

ステップＳ１０６において、情報取得部２２は、上記ステップＳ１０４で取得された地震の速度データのＳＮ比が所定の条件を満たすか否かを判定する。地震の速度データのＳＮ比が所定の条件を満たす場合には、ステップＳ１０８へ進む。一方、地震の速度データのＳＮ比が所定の条件を満たさない場合には、地下構造推定処理ルーチンを終了する。

ステップＳ１０８において、情報取得部２２は、モデル記憶部２０に格納された地下構造モデルのパラメータを表すパラメータベクトルを取得する。

ステップＳ１１０において、推定部２４は、上記ステップＳ１０８で取得された地下構造モデルのパラメータベクトルθ及びシミュレーションに応じて予測される地震データと、上記ステップＳ１０４で取得された地震の速度データｙ_ｔとの差分に応じて、地下構造モデルのパラメータベクトルθを推定する。ステップＳ１１０の処理は、図７に示すデータ同化処理ルーチンによって実現される。

＜データ同化処理ルーチン＞
ステップＳ２００において、時刻ｔに０を設定する。

ステップＳ２０２において、予測部２４０は、上記式（６）の条件付き確率分布に従って、拡大状態ベクトルの初期値ｘ_０を設定する。

ステップＳ２０４において、予測部２４０は、格子点間の距離が近いほど、一方の格子点の地下方向の層厚と他方の格子点の地下方向の層厚とが近い値になるように、上記式（１２）の分散共分散行列Ｖ＾_{ｔ｜ｔ−１}の層厚Ｈに関する項に対し、格子点の位置に応じた相関を設定する。

ステップＳ２０６において、予測部２４０は、前時刻ｔ−１の各アンサンブルｘ^（ｉ） _{ｔ−１｜ｔ−１}とシミュレーションｆ_ｔとに基づいて、上記式（１０）に従って、現時刻ｔの拡大状態ベクトルを表す各アンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}を予測する。

ステップＳ２０８において、修正部２４２は、上記ステップＳ２０６で予測された時刻ｔの各アンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}から生成される地震データと、時刻ｔの地震の速度データｙ_ｔとに基づいて、上記式（１５）に従って、時刻ｔの各アンサンブルｘ^（ｉ） _{ｔ｜ｔ−１}を修正する。そして、修正部２４２は、拡大状態ベクトルを表す各アンサンブルｘ^（ｉ） _ｔ｜ｔを得る。

ステップＳ２１０において、修正部２４２は、地震の速度データｙ_ｔの終端時刻Ｔまで、上記ステップＳ２０６の一期先予測及び上記ステップＳ２０８のフィルタリングが繰り返されたか否かを判定する。終端時刻Ｔまで、一期先予測及びフィルタリングが繰り返された場合には、ステップＳ２１４へ進む。一方、終端時刻Ｔまで、一期先予測及びフィルタリングが繰り返されていない場合には、ステップＳ２１２で時刻を１インクリメントしてステップＳ２０６へ戻る。

ステップＳ２１４において、平滑化部２４４は、上記ステップＳ２０８で各時刻において得られた推定結果に含まれる、地下構造モデルのパラメータベクトルθを平滑化することにより、パラメータベクトルθを取得する。

ステップＳ２１６において、平滑化部２４４は、上記ステップＳ２１４で得られたパラメータベクトルθを結果として出力する。

次に、図６のステップＳ１１２において、平滑化部２４４は、上記ステップＳ１１０で出力されたパラメータベクトルθを、モデル記憶部２０へ格納する。モデル記憶部２０に格納されたパラメータベクトルθは、次回の地震の観測データが得られたときのデータ同化の際に用いられる。

ステップＳ１１４において、表示部１６は、モデル記憶部２０に記憶されたパラメータベクトルθ等の推定結果を表示して、地下構造推定処理ルーチンを終了する。

以上詳細に説明したように、本実施形態では、地下構造モデルに対応する地点における地震の観測データの時系列を入力とし、地下構造モデルのパラメータ及び地震予測モデルから予測される地震データと、地震の観測データとの差分に応じて、地震の状態及び地下構造モデルのパラメータを修正することを繰り返すことにより、地震の観測データを用いて地下構造を推定することができる。

地下構造モデルのパラメータ推定に対してデータ同化を適用しようとする場合、例えば、逐次フィルタリング手法の一例であるアンサンブルカルマンフィルタを用いるときには、拡大状態ベクトルの候補を表すアンサンブルが複数生成される。しかし、地下構造モデルのパラメータ数は多いため、パラメータに応じて生成されるアンサンブルに対する計算量も増大してしまうという課題があった。また、各アンサンブルは確率分布に応じてランダムに生成されるため、連続性を有する地下構造モデルのパラメータを推定したとしても、精度よく推定することができないという課題があった。

そこで、本実施形態では、データ同化により地下構造を推定する際に、地下構造モデルの各層に対してパラメータを付与することにより、パラメータ数を減少させる。これにより、データ同化により地下構造を推定する際に、計算量を低減させることができる。

また、本実施形態では、地下構造モデルのパラメータに対して相関を設定する。これにより、地下構造の連続性が考慮され、地下構造を精度よく推定することができる。

また、地震の観測データと整合度の高い地下構造モデルが推定されるため、地震動の予測精度を向上させることができる。また、地下構造モデルの修正のための追加のボーリング調査などが不要となる。

また、本実施形態では逐次データ同化を用いるため、地震の観測データが増加しても、増加分の観測データのみを用いて地下構造モデルのパラメータ修正を行うだけで良く、非逐次データ同化と比較して計算量が少なくて済む。

＜実験例＞
本実施形態の効果の確認のため、２次元ＳＨ波動場における数値実験を行った。なお、逐次データ同化手法の一つであるアンサンブルカルマンフィルタを用いて解析を行った。

［実験の概要］
図８に示す２次元地下構造モデルを対象として数値実験を行った。まず、真のモデル（真値）に対してシミュレーションを行い、地表観測点（図８中の三角）の模擬地震動を作成する。また、以下の表１に示すように、真の地下構造モデルに１割の誤差を与えた初期モデル（初期値）を設定する。次に初期モデルに対するシミュレーションと模擬地震動とのデータ同化を実施し、初期モデルが真のモデルに近づくことを確認した。

［層に相関を設定した場合と層に相関を設定しない場合との比較］

層に対して相関を与えた場合の数値実験の結果を図９（Ａ）に示す。図の横軸は時刻ステップ、縦軸は真値に対する推定値の比を示しており、比が１に近づくほど真値に近づくことを表す。層厚ＨとＳ波のせん断波速度Ｖ_Ｓとついて、逐次修正され真値に近づいており、データ同化を用いた地下構造推定の効果が確認されている。

一方、層に対して相関を設定しない場合の数値実験の結果を図９（Ｂ）に示す。相関が設定された場合と比較すると、相関を設定しない場合の層厚Ｈの推定精度が悪い。また、Ｓ波のせん断波速度Ｖ_Ｓについても相関が設定された場合の方が推定精度は良く、層厚に相関を与えることの有効性が確認された。

なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、上記実施形態では、各層の物性値のパラメータ（ρ，Ｖ_Ｓ，Ｖ_Ｐ，Ｑ_Ｓ，Ｑ_Ｐ）と、層厚のパラメータＨとをデータ同化により同時に推定する場合を例に説明したがこれに限定されるものではない。例えば、各パラメータのうちの所定のパラメータを推定した後、所定のパラメータ以外のパラメータを推定するようにしてもよい。例えば、各層の層厚Ｈと速度Ｖ_Ｓとを推定した後に、密度ρと減衰定数Ｑ_Ｓ及びＱ_Ｐを推定するようにしてもよい。

また、上記実施形態では、図４に示されるように格子点間の水平距離（ｘｙ方向）に応じて相関を設定する場合を例に説明したが、これに限定されるものではなく、ｚ方向に対して相関を設定するようにしてもよい。例えば、ｚ方向の位置が深いほど層の固さは硬くなるように相関を設定するようにしてもよい。これにより、例えばアンサンブルカルマンフィルタにおいて生成される各アンサンブルにおいて、特定の層よりも、特定の層の下に位置する層の方が軟らかいことを表すアンサンブルの生成等が抑制される。

また、上記実施形態では、層に関するパラメータの一例として層厚Ｈに対して相関を設定する場合を例に説明したが、これに限定されるものではない。例えば、層厚Ｈ以外のパラメータに対して相関を設定してもよい。具体的には、例えば、地下構造モデルのｚ方向の位置が深いほど速度Ｖ_Ｓが大きくなるように、相関を設定してもよい。また、例えば、速度Ｖ_Ｓと密度ρとの関係を考慮して、速度Ｖ_Ｓが大きいほど密度ρも大きくなるように相関を設定してもよい。また、地下構造モデルのｚ方向の位置が深いほど減衰係数Ｑが大きくなるように、相関を設定してもよい。また、速度Ｖ_Ｓと速度Ｖ_Ｐとの間の関係を考慮して、相関を設定してもよい。また、実際の地盤調査結果、各パラメータ間の関係式、及び既往の研究の少なくとも１つに応じて、相関を設定してもよい。

また、上記実施形態では、地下構造モデルのパラメータとして、各層の物性値（ρ，Ｖ_Ｓ，Ｖ_Ｐ，Ｑ_Ｓ，Ｑ_Ｐ）と層厚Ｈとをデータ同化により推定する場合を例に説明したがこれに限定されるものではない。本実施形態の地下構造モデルのパラメータの一部のみを推定してもよく、例えば、層厚Ｈ及び速度Ｖ_Ｓのみを推定してもよい。また、層厚Ｈ及び速度Ｖ_Ｓのみを推定した後に、減衰定数Ｑ_Ｓを推定してもよい。

また、上記実施形態では、地下構造推定処理ルーチンにおいて、新たな観測データに基づくデータ同化が行われる毎に、パラメータベクトルに対して相関を設定する場合を例に説明したが、これに限定されるものではない。例えば、はじめて得られた地震の観測データに基づきデータ同化を行う際にのみ、パラメータベクトルに対して相関を設定してもよい。

また、上記実施形態では、データ同化の逐次フィルタリング手法の一例としてアンサンブルカルマンフィルタを用いる場合を例に説明したがこれに限定されるものではなく、他の逐次フィルタリング手法を用いてもよい。例えば、逐次フィルタリング手法として粒子フィルタを用いてもよい。

また、上記ではプログラムが記憶部（図示省略）に予め記憶（インストール）されている態様を説明したが、プログラムは、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ及びマイクロＳＤカード等の記録媒体の何れかに記録されている形態で提供することも可能である。

１０ 地下構造推定装置
１２ 受付部
１４ コンピュータ
１６ 表示部
１８ 観測データ記憶部
２０ モデル記憶部
２２ 情報取得部
２４ 推定部
２６ 地下構造モデル
２４０ 予測部
２４２ 修正部
２４４ 平滑化部
θ パラメータベクトル

Claims (5)

1. 地下構造モデルと、前記地下構造モデルに対応する地点における各時刻の地震の観測データの時系列とを入力とし、
   前記地下構造モデルのパラメータ及び地震予測モデルから予測される前記時刻の地震の状態から得られる地震データと、前記時刻における地震の観測データとの差分に応じて、前記時刻の地震の状態及び前記地下構造モデルのパラメータを修正することを繰り返す推定部
   を含み、
   前記地下構造モデルの前記パラメータを表すパラメータベクトルと前記地震による各格子点の速度及び応力を含む状態ベクトルとを有する拡大状態ベクトルと、前記観測データを表す各時刻の前記地震の地震動データとを入力とし、
   前記推定部は、
   時刻ｔ−１の前記拡大状態ベクトルと前記地震予測モデルとに基づいて、時刻ｔの前記拡大状態ベクトルを予測する予測部と、
   前記予測部によって予測された時刻ｔの前記拡大状態ベクトルから生成される前記地震データと、時刻ｔの前記地震の地震動データとに基づいて、時刻ｔの前記拡大状態ベクトルを修正する修正部と、
   前記観測データの終端時刻Ｔまで、前記予測部による予測及び前記修正部による修正を繰り返し、各時刻において得られた前記修正部による推定結果に含まれる、前記地下構造モデルの前記パラメータを平滑化することにより、前記拡大状態ベクトルに含まれる前記地下構造モデルの前記パラメータを取得する平滑化部と、
   を含む地下構造推定装置。
2. 前記地下構造モデルのパラメータは、前記地下構造モデルの地表の各格子点に対する、前記格子点の地下方向に位置する各層に関する前記パラメータを含む、
   請求項１に記載の地下構造推定装置。
3. 各格子点の前記層に関するパラメータは、前記格子点間で、前記格子点の位置に応じた相関が予め設定される、
   請求項２に記載の地下構造推定装置。
4. 前記層に関するパラメータは前記層の層厚を含み、
   前記層厚は、前記格子点間で、前記格子点の位置に応じた相関が予め設定される、
   請求項３に記載の地下構造推定装置。
5. 前記地下構造モデルの前記格子点間の各々について、前記格子点間の距離が近いほど、一方の格子点の地下方向に位置する各層に関する前記パラメータと、他方の格子点の地下方向に位置する各層に関する前記パラメータとが対応するように予め設定される、
   請求項３又は請求項４に記載の地下構造推定装置。