(57)【要約】【課題】地震の観測データを用いて地下構造を推定することができる地下構造推定装置を得る。【解決手段】地下構造推定装置１０は、地下構造モデルに対応する地点における各時刻の地震の観測データの時系列を入力とし、地下構造モデルのパラメータ及び地震予測モデルから予測される地震データと、地震の観測データとの差分に応じて、地下構造モデルのパラメータを修正することを繰り返す。【選択図】図１
2018128929.tif 000002
【特許請求の範囲】
【請求項１】
地下構造モデルと、前記地下構造モデルに対応する地点における各時刻の地震の観測データの時系列とを入力とし、
前記地下構造モデルのパラメータ及び地震予測モデルから予測される前記時刻の地震の状態から得られる地震データと、前記時刻における地震の観測データとの差分に応じて、前記時刻の地震の状態及び前記地下構造モデルのパラメータを修正することを繰り返す推定部
を含む地下構造推定装置。
【請求項２】
前記地下構造モデルのパラメータは、前記地下構造モデルの地表の各格子点に対する、前記格子点の地下方向に位置する各層に関する前記パラメータを含む、
請求項１に記載の地下構造推定装置。
【請求項３】
各格子点の前記層に関するパラメータは、前記格子点間で、前記格子点の位置に応じた相関が予め設定される、
請求項２に記載の地下構造推定装置。
【請求項４】
前記層に関するパラメータは前記層の層厚を含み、
前記層厚は、前記格子点間で、前記格子点の位置に応じた相関が予め設定される、
請求項３に記載の地下構造推定装置。
【請求項５】
前記地下構造モデルの前記格子点間の各々について、前記格子点間の距離が近いほど、一方の格子点の地下方向に位置する各層に関する前記パラメータと、他方の格子点の地下方向に位置する各層に関する前記パラメータとが対応するように予め設定される、
請求項３又は請求項４に記載の地下構造推定装置。
【請求項６】
前記地下構造モデルの前記パラメータを表すパラメータベクトルと前記地震による各格子点の速度及び応力を含む状態ベクトルとを有する拡大状態ベクトルと、前記観測データを表す各時刻の前記地震の地震動データとを入力とし、
前記推定部は、
時刻ｔ−１の前記拡大状態ベクトルと前記地震予測モデルとに基づいて、時刻ｔの前記拡大状態ベクトルを予測する予測部と、
前記予測部によって予測された時刻ｔの前記拡大状態ベクトルから生成される前記地震データと、時刻ｔの前記地震の地震動データとに基づいて、時刻ｔの前記拡大状態ベクトルを修正する修正部と、
前記観測データの終端時刻Ｔまで、前記予測部による予測及び前記修正部による修正を繰り返し、各時刻において得られた前記修正部による推定結果に含まれる、前記地下構造モデルの前記パラメータを平滑化することにより、前記拡大状態ベクトルに含まれる前記地下構造モデルの前記パラメータを取得する平滑化部と、
を含む請求項１〜請求項５の何れか１項に記載の地下構造推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、地下構造推定装置に関する。
【背景技術】
【０００２】
従来、地盤・岩盤における透水係数を推定する技術が知られている（例えば、特許文献１を参照）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００１‐３２４１９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
上記特許文献１では、試験用ボーリング孔からの流量計測を必須とせずに、間隙水圧のみの計測から地盤・岩盤における透水係数を推定する。しかし、上記特許文献１の技術では、間隙水圧の計測は必要であり、かつ地下構造としては透水係数のみしか推定することができない。
【０００５】
本発明は上記事実を考慮して、地震の観測データを用いて地下構造を推定することができる地下構造推定装置を提供することを目的とする。
【課題を解決するための手段】
【０００６】
上記目的を達成するために、本発明の地下構造推定装置は、地下構造モデルと、前記地下構造モデルに対応する地点における各時刻の地震の観測データの時系列とを入力とし、前記地下構造モデルのパラメータ及び地震予測モデルから予測される前記時刻の地震の状態から得られる地震データと、前記時刻における地震の観測データとの差分に応じて、前記時刻の地震の状態及び前記地下構造モデルのパラメータを修正することを繰り返す推定部を含んで構成されている。これにより、地震の観測データを用いて地下構造を推定することができる。
【０００７】
本発明の前記地下構造モデルのパラメータは、前記地下構造モデルの地表の各格子点に対する、前記格子点の地下方向に位置する各層に関する前記パラメータを含むようにすることができる。これにより、地下構造を推定する際に計算量を低減させることができる。
【０００８】
本発明の各格子点の前記層に関するパラメータは、前記格子点間で、前記格子点の位置に応じた相関が予め設定されるようにすることができる。これにより、地下構造の連続性を考慮して、地下構造を精度よく推定することができる。
【０００９】
本発明の前記層に関するパラメータは前記層の層厚を含み、前記層厚は、前記格子点間で、前記格子点の位置に応じた相関が予め設定されるようにすることができる。これにより、地下構造の層厚の連続性を考慮して、地下構造を精度よく推定することができる。
【００１０】
本発明では、前記地下構造モデルの前記格子点間の各々について、前記格子点間の距離が近いほど、一方の格子点の地下方向に位置する各層に関する前記パラメータと、他方の格子点の地下方向に位置する各層に関する前記パラメータとが対応するように予め設定されるようにすることができる。これにより、地下構造の格子点間の距離を考慮して、地下構造を精度よく推定することができる。
【００１１】
本発明の地下構造推定装置は、前記地下構造モデルの前記パラメータを表すパラメータベクトルと前記地震による各格子点の速度及び応力を含む状態ベクトルとを有する拡大状態ベクトルと、前記観測データを表す各時刻の前記地震の地震動データとを入力とし、前記推定部は、時刻ｔ−１の前記拡大状態ベクトルと前記地震予測モデルとに基づいて、時刻ｔの前記拡大状態ベクトルを予測する予測部と、前記予測部によって予測された時刻ｔの前記拡大状態ベクトルから生成される前記地震データと、時刻ｔの前記地震の地震動データとに基づいて、時刻ｔの前記拡大状態ベクトルを修正する修正部と、前記観測データの終端時刻Ｔまで、前記予測部による予測及び前記修正部による修正を繰り返し、各時刻において得られた前記修正部による推定結果に含まれる、前記地下構造モデルの前記パラメータを平滑化することにより、前記拡大状態ベクトルに含まれる前記地下構造モデルの前記パラメータを取得する平滑化部と、を含むようにすることができる。これにより、逐次フィルタリング手法を用いて地下構造を精度よく推定することができる。
【発明の効果】
【００１２】
本発明によれば、地下構造モデルに対応する地点における地震の観測データの時系列を入力とし、地下構造モデルのパラメータ及び地震予測モデルから予測される地震データと、地震の観測データとの差分に応じて、地震の状態及び地下構造モデルのパラメータを修正することを繰り返すことにより、地震の観測データを用いて地下構造を推定することができる、という効果が得られる。
【図面の簡単な説明】
【００１３】
【図１】本発明の実施形態に係る地下構造推定装置の概略構成を示すブロック図である。
【図２】地下構造モデルを説明するための説明図である。
【図３】本発明の実施形態に係る地下構造推定装置の推定部の詳細な構成の一例を示す図である。
【図４】逐次データ同化において用いられる分散共分散行列に対して相関が設定された場合の一例を示す図である。
【図５】地下構造モデルの層に対して相関を与えない場合のアンサンブルカルマンフィルタの各アンサンブルの一例と、地下構造モデルの層に対して相関を与えた場合のアンサンブルカルマンフィルタの各アンサンブルの一例とを示す図である。
【図６】地下構造推定処理ルーチンの一例を示す図である。
【図７】データ同化処理ルーチンの一例を示す図である。
【図８】シミュレーション実験を説明するための説明図である。
【図９】シミュレーション実験の実験結果を示す図である。
【発明を実施するための形態】
【００１４】
以下、本発明の実施形態について詳細に説明する。
【００１５】
＜地下構造推定装置のシステム構成＞
【００１６】
図１は、本発明の実施形態に係る地下構造推定装置の構成の一例を示すブロック図である。地下構造推定装置１０は、機能的には、図１に示されるように、受付部１２、コンピュータ１４、及び表示部１６を含んだ構成で表すことができる。
【００１７】
受付部１２は、地下構造推定装置の外部から入力された情報を受け付ける。例えば、受付部１２は、新たに地震が発生した場合、地震の観測データの一例である震源パラメータ及び各観測地点における各時刻の地震の速度データを受け付ける。震源パラメータは、気象庁等からの地震情報に基づいて入力される。そして、受付部１２は、震源パラメータ及び各観測地点における各時刻の地震の速度データを、後述する観測データ記憶部１８へ格納する。なお、地震の速度データは、地震の地震動データの一例である。
【００１８】
コンピュータ１４は、ＣＰＵ（Central Processing Unit）、各処理ルーチンを実現するためのプログラム等を記憶したＲＯＭ（Read Only Memory）、データを一時的に記憶するＲＡＭ（Random Access Memory）、記憶手段としてのメモリ、ネットワークインタフェース等を含んで構成されている。
【００１９】
コンピュータ１４は、機能的には、図１に示されるように、観測データ記憶部１８と、モデル記憶部２０と、情報取得部２２と、推定部２４とを含んで構成される。
【００２０】
観測データ記憶部１８には、地震毎に、震源パラメータ及び各観測地点における各時刻の地震の速度データが格納される。新たな地震が発生する毎に、当該地震の観測データが観測データ記憶部１８へ格納される。
【００２１】
モデル記憶部２０には、地下構造を表す地下構造モデルと、地下構造モデルのパラメータを表すパラメータベクトルとが格納される。
【００２２】
本実施形態では、データ同化により地下構造モデルのパラメータを推定する。そのため、本実施形態では、観測データ記憶部１８に格納された地震の速度データと、モデル記憶部２０に格納された地下構造モデルのパラメータベクトル及び地震予測モデルから予測される地震データとの差分を修正するように、地下構造モデルのパラメータを推定する。
【００２３】
本実施の形態で用いる地下構造モデルとパラメータとについて、以下説明する。
【００２４】
地震動予測において利用されている有限差分法を用いたシミュレーションでは、地下構造モデルの各格子点に物性値（例えば、伝播速度、密度、及び減衰定数等）を表すパラメータを付与することにより、地下構造を表現する。
【００２５】
例えば、図２の地下構造モデル２５のように、地下構造モデル２５の各格子点（ｉ，ｊ，ｋ）には、密度ρ、Ｓ波のせん断波速度ＶＳ、Ｐ波のせん断波速度ＶＰ、Ｓ波の減衰定数ＱＳ、及びＰ波の減衰定数ＱＰの物性値がパラメータとして付与される。
【００２６】
しかし、地下構造モデルの各格子点に対してパラメータを付与し、データ同化により各パラメータを推定する場合、推定対象のパラメータが多くなりすぎてしまい、計算量が増加してしまう。また、パラメータの増加により、パラメータの推定精度が低下する可能性がある。
【００２７】
そこで、本実施形態では、データ同化に用いる地下構造モデルにおいて、地表の各格子点の地下方向に位置する各層に対してパラメータを付与する。具体的には、図２の地下構造モデル２６のように、地表の格子点毎に、地下構造モデル２６の各層に対する、層の層厚Ｈをパラメータとして付与する。また、各層に対して物性値（ρ，ＶＳ，ＶＰ，ＱＳ，ＱＰ）をパラメータとして付与する。
【００２８】
ここで、ｘ方向の格子の区切り数をＩ、ｙ方向の格子の区切り数をＪ、ｚ方向の格子の区切り数をＫとすると、上記図２の地下構造モデル２５の場合、総パラメータ数は、以下の式（１）のようになる。
【００２９】
総パラメータ数＝１格子点当たりのパラメータ数×格子点数
＝５×Ｉ×Ｊ×Ｋ  ［個］                    （１）
【００３０】
一方、上記図２に示される本実施形態の地下構造モデル２６の場合、層数をＬとすると、総パラメータ数は、以下の式（２）のようになる。
【００３１】
総パラメータ数
＝１層当たりのパラメータ数×層数＋１格子点当たりの層数×地表の格子点数
＝５×Ｌ＋Ｉ×Ｊ×（Ｌ−１）  ［個］                      （２）
【００３２】
地下構造モデル２５では、上記式（１）から分かるように、格子点毎に物性値のパラメータ（ρ，ＶＳ，ＶＰ，ＱＳ，ＱＰ）が付与されると、推定対象のパラメータ数が膨大な数となる。一方、本実施形態の地下構造モデル２６では、上記式（２）から分かるように、物性値のパラメータ（ρ，ＶＳ，ＶＰ，ＱＳ，ＱＰ）は層毎に付与され、層厚のパラメータＨは地表の格子点毎に付与される。そのため、本実施形態の地下構造モデル２６では、推定対象のパラメータ数が減少し、データ同化によるパラメータ推定の際の計算量が抑制される。
【００３３】
情報取得部２２は、観測データ記憶部１８に格納された観測データを取得する。そして、情報取得部２２は、観測データの震源パラメータに含まれる震源位置等の情報に基づき、当該観測データの地震の速度データを、地下構造モデルのパラメータの推定に用いるか否かを判定する。
【００３４】
例えば、情報取得部２２は、震源パラメータに含まれる震源位置と地下構造モデルに対応する地点との間の距離が、予め設定した閾値以上である場合には、地下構造モデルの領域内ではないと判定し、地下構造モデルのパラメータの推定に用いないと判定する。
【００３５】
一方、情報取得部２２は、震源パラメータに含まれる震源位置と地下構造モデルに対応する地点との間の距離が、予め設定した閾値未満である場合には、地下構造モデルの領域内であると判定し、地下構造モデルのパラメータの推定に用いると判定する。地下構造モデルに対応する地点と関係ない領域の地震の速度データについては、地下構造のパラメータを推定する際には適切でないため除外される。
【００３６】
また、情報取得部２２は、観測データである地震の速度データのＳＮ比が所定の条件を満たすか否かに応じて、当該観測データの地震の速度データを地下構造モデルのパラメータの推定に用いるか否かを判定する。ＳＮ比が所定の条件を満たすか否かは、地震の速度データから計算されるフーリエ振幅スペクトルの形状に基づき判定される。
【００３７】
情報取得部２２は、観測データである地震の速度データの低周波成分がノイズレベルを下回る場合には、地下構造のパラメータを推定する際には適切でないと判定し、当該地震の速度データを地下構造モデルのパラメータの推定に用いないことを判定する。一方、情報取得部２２は、地震の速度データの低周波成分がノイズレベルを上回る場合には、当該地震の速度データを地下構造モデルのパラメータの推定に用いることを判定する。
【００３８】
また、情報取得部２２は、モデル記憶部２０に格納された地下構造モデルのパラメータを表すパラメータベクトルを取得する。情報取得部２２により取得される地下構造モデルのパラメータベクトルは、前回までの地震の速度データに基づいてデータ同化により得られたものである。
【００３９】
推定部２４は、情報取得部２２によって取得された地下構造モデルのパラメータベクトル及び地震予測モデルから予測される地震データと、情報取得部２２によって取得された地震の速度データとの差分に応じて、地下構造モデルのパラメータベクトルを修正することを繰り返す。
【００４０】
具体的には、地震予測モデルの一例である有限差分法によるシミュレーションによって予測される地震の状態から得られる地震データと、観測データである地震の速度データとの差分に応じて、データ同化により地下構造モデルのパラメータを推定する。
【００４１】
本実施形態では、データ同化の一例として、逐次データ手法であるアンサンブルカルマンフィルタを用いる場合を例に説明する。逐次データ手法では、時系列データが得られる毎に同化が行われる。
【００４２】
データ同化の処理の流れを以下に示す。本実施形態では、データ同化に用いるシステムモデルは以下の式（３）及び式（４）によって表される。
【００４３】
【数１】
2018128929.tif 000003
（３）
2018128929.tif 000004
（４）
【００４４】
上記式（３）及び式（４）において、ｔは時刻を表し、ｘｔは時刻ｔにおける拡大状態ベクトルを表し、ｙｔは時刻ｔにおける観測データを表す。ｆｔはシステムモデルであるシミュレーションを表す。また、観測ノイズｗｔは正規分布Ｎ（０，Ｒｔ）に従い、システムノイズｖｔは任意の分布に従う。Ｈｔは観測行列である。
【００４５】
上記式（３）及び式（４）を条件付き分布で表現すると、以下の式（５）〜式（７）に示されるような一般状態空間モデルとなる。
【００４６】
【数２】
2018128929.tif 000005
（５）
2018128929.tif 000006
（６）
2018128929.tif 000007
（７）
【００４７】
ここで、ｘ０は地震の観測データと同化する前の拡大状態ベクトルを表す。また、地下構造モデルのパラメータベクトルθについては、既往のモデルや構造探査結果や地質情報に応じて設定される。
【００４８】
次に、参考文献（長尾大道、「データ同化の基礎理論とその応用」、自動チューニング研究会、第６回自動チューニング技術の現状と応用に関するシンポジウム（ＡＴＴＡ２０１４）スライド資料）を参照して、アンサンブルカルマンフィルタについて説明する。
【００４９】
［１］一期先予測
アンサンブルカルマンフィルタの一期先予測において、拡大状態ベクトルｘｔに対応する各アンサンブルｘ（ｉ）ｔ−１｜ｔ−１をシミュレーションｆｔによって時間発展させ、一期先予測分布のアンサンブルｘ（ｉ）ｔ｜ｔ−１を得る。
【００５０】
時刻ｔ−１におけるフィルタ分布は、以下の式（９）によって近似表現される。なお、δ（・）はデルタ関数を表す。また、各アンサンブルｘ（ｉ）ｔ−１｜ｔ−１を、システムモデルであるシミュレーションｆｔによって時間発展させた各アンサンブルｘ（ｉ）ｔ｜ｔ−１は、以下の式（１０）によって表される。また、式（１０）の各アンサンブルｘ（ｉ）ｔ｜ｔ−１によって、一期先予測分布は以下の式（１１）によって近似表現される。
【００５１】
【数３】
2018128929.tif 000008
（８）
2018128929.tif 000009
（９）
2018128929.tif 000010
（１０）
2018128929.tif 000011
（１１）
【００５２】
［２］フィルタリング
次に、アンサンブルカルマンフィルタのフィルタリングにおいて、一期先予測分布のアンサンブルｘ（ｉ）ｔ｜ｔ−１からフィルタ分布のアンサンブルｘ（ｉ）ｔ｜ｔを得る。
【００５３】
一期先予測分布のアンサンブルｘ（ｉ）ｔ｜ｔ−１の分散共分散行列Ｖ＾ｔ｜ｔ−１は、以下の式（１２）によって表される。なお、式中の「’」は行列の転置を表す。また、観測ノイズｗ（ｉ）ｔの分散共分散行列Ｒ＾ｔは、以下の式（１３）によって表される。
【００５４】
また、カルマンゲインＫ＾ｔは、以下の式（１４）によって表される。そして、フィルタ分布のアンサンブルｘ（ｉ）ｔ｜ｔは、以下の式（１５）によって表される。なお、ｙｔは時刻ｔにおける観測データを表す。
【００５５】
【数４】
2018128929.tif 000012
（１２）
2018128929.tif 000013
（１３）
2018128929.tif 000014
（１４）
2018128929.tif 000015
（１５）
【００５６】
［３］平滑化
アンサンブルカルマンフィルタの平滑化において、一期先予測及びフィルタリングの繰り返しにより得られたアンサンブルを平滑化する。具体的には、時刻ステップｔ＝Ｔまでの観測データから得られた各アンサンブルに基づいて、平滑化分布ｐ（ｘ０｜ｙ１：Ｔ）を計算することにより、拡大状態ベクトルを推定する。
【００５７】
本実施の形態におけるアンサンブルカルマンフィルタを用いた各処理について、以下、具体的に説明する。
【００５８】
本実施の形態における拡大状態ベクトルｘｔは、シミュレーションによって得られる地下構造モデルの各格子点の速度ｖｉ，ｊ，ｋ及び応力τｉ，ｊ，ｋを要素とする状態ベクトルｘ〜ｔと、地下構造モデルの各層ｌの物性値を表すパラメータ（ρｌ，ＶＳ，ｌ，ＶＰ，ｌ，ＱＳ，ｌ，ＱＰ，ｌ）及び各層ｌの層厚Ｈｉ，ｊ，ｌを要素とするパラメータベクトルθとを含む。
【００５９】
また、本実施形態の状態ベクトルｘ〜ｔは以下の式（１６）によって表される。状態ベクトルｘ〜ｔのｑｉ，ｊ，ｋは、シミュレーションによって得られる各格子点の速度を表す。また、状態ベクトルｘ〜ｔのτｉ，ｊ，ｋは、シミュレーションによって得られる各格子点の応力を表す。状態ベクトルｘ〜ｔは地震の状態の一例である。
【００６０】
本実施形態のパラメータベクトルθは以下の式（１７）によって表される。パラメータベクトルθには、各層ｌの物性値のパラメータ（ρｌ，ＶＳ，ｌ，ＶＰ，ｌ，ＱＳ，ｌ，ＱＰ，ｌ）と、層厚のパラメータＨｉ，ｊ，ｌとが含まれる。
【００６１】
また、地震データの一例である本実施形態の拡大状態ベクトルｘｔは、以下の式（１８）によって表される。また、Ｍ個の観測点についての時刻ｔの地震の速度データｙｔは、以下の式（１９）によって表される。地震の速度データｙｔの各要素ａｍは、地点ｍにおいて観測された地震の速度データを表す。
【００６２】
【数５】
2018128929.tif 000016
【００６３】
推定部２４は、図３に示されるように、予測部２４０と、修正部２４２と、平滑化部２４４とを含んで構成されている。
【００６４】
予測部２４０は、情報取得部２２によって取得された地下構造モデルのパラメータベクトルθを取得する。また、予測部２４０は、情報取得部２２によって取得された地震の速度データｙｔを取得する。
【００６５】
次に、予測部２４０は、各値の初期値を設定する。具体的には、予測部２４０は、上記式（６）の条件付き確率分布に従って、拡大状態ベクトルの初期値ｘ０を設定する。また、予測部２４０は、拡大状態ベクトルの初期値ｘ０を設定する際、パラメータベクトルθのうちの層厚Ｈについて、格子点間で、格子点の位置に応じた相関を設定する。具体的には、予測部２４０は、地下構造モデルの格子点間の各々について、格子点間の距離が近いほど、一方の格子点の地下方向に位置する各層の層厚と、他方の格子点の地下方向に位置する各層の層厚とが近い値になるように設定する。より詳細には、上記式（１２）の分散共分散行列Ｖ＾ｔ｜ｔ−１の層厚Ｈに関する項について、格子点の位置に応じた相関を設定することにより、層厚の相関を設定する。
【００６６】
図４に、層厚Ｈに関する項についての分散共分散行列Ｖ＾ｔ｜ｔ−１の一例を示す。図４（Ａ）の横軸α及び縦軸βは行列の各要素の位置を表し、明るい領域ほど相関値が高く暗い領域ほど相関値が低い。図４（Ｂ）は、図４（Ａ）を３次元で表したものであり、水平面における横軸α及び縦軸βは行列の各要素の位置を表し、垂直方向は相関値を表す。
【００６７】
図４（Ａ）及び（Ｂ）に示されるように、本実施形態では、各格子点について、格子点間の距離が小さいほど相関値が高くなるように、かつ格子点間の距離が大きいほど相関値が低くなるように、層厚Ｈに関しての分散共分散行列Ｖ＾ｔ｜ｔ−１を設定する。そして、相関が設定された分散共分散行列Ｖ＾ｔ｜ｔ−１に応じて各アンサンブルが生成される。これにより、地下構造の連続性が反映される。
【００６８】
地下構造は、断層が存在しない限り層の形状は滑らかに変化する。本実施形態では、この点を考慮し、同一の層に対応する層厚に相関を与える。
【００６９】
図５に、相関を設定しない場合のアンサンブル１〜５の一例と相関を設定した場合のアンサンブル１〜５の一例とを示す。なお、図中のＸgridは上記図２のｘ方向に対応し、Ｚgridは上記図２のｚ方向に対応する。図５（Ａ）に示されるように、相関が設定されない場合のアンサンブル１〜５においては、Ｚgridの値はＸgridの位置に関係なくランダムな値をとる。しかし、地下構造における層は、水平方向において連続していると考えられる。そのため、各アンサンブルのうちの層厚がランダムな値をとりながらデータ同化が行われると、地下構造モデルのパラメータを精度よく推定することができない。
【００７０】
一方、図５（Ｂ）に示されるように、相関が設定された場合のアンサンブル１〜５においては、Ｚgridの値はＸgridの位置に応じた値となり、例えばＸgridが隣接する点同士では近い値をとる。そのため、層厚Ｈに関する項についての分散共分散行列Ｖ＾ｔ｜ｔ−１に対して相関を設定することにより、地下構造の層の連続性が反映される。
【００７１】
次に、予測部２４０は、時刻ｔ−１の拡大状態ベクトルを表す各アンサンブルｘ（ｉ）ｔ−１｜ｔ−１と地震予測モデルの一例である有限差分法のシミュレーションｆｔとに基づいて、上記式（１０）に従って、時刻ｔの拡大状態ベクトルを表す各アンサンブルｘ（ｉ）ｔ｜ｔ−１を予測する。
【００７２】
修正部２４２は、予測部２４０によって予測された時刻ｔの拡大状態ベクトルを表す各アンサンブルｘ（ｉ）ｔ｜ｔ−１から生成される地震データと、時刻ｔの地震の速度データｙｔとに基づいて、上記式（１５）に従って、時刻ｔの拡大状態ベクトルを表す各アンサンブルｘ（ｉ）ｔ｜ｔ−１を修正する。そして、修正部２４２は、拡大状態ベクトルを表す各アンサンブルｘ（ｉ）ｔ｜ｔを得る。
【００７３】
そして、地震の速度データの終端時刻Ｔまで、予測部２４０による一期先予測及び修正部２４２によるフィルタリングが繰り返される。
【００７４】
平滑化部２４４は、各時刻において得られた修正部２４２による推定結果に含まれる、地下構造モデルのパラメータベクトルθを平滑化することにより、地下構造モデルのパラメータベクトルθを取得する。パラメータベクトルθには、各層の物性値のパラメータ（ρ，ＶＳ，ＶＰ，ＱＳ，ＱＰ）と、各層の層厚のパラメータＨとが含まれているため、データ同化により地下構造のパラメータが得られたことになる。
【００７５】
そして、平滑化部２４４は、平滑化によって得られた地下構造モデルのパラメータベクトルを、モデル記憶部２０へ格納する。モデル記憶部２０に格納されたパラメータベクトルは、次回の地震の観測データが得られたときのデータ同化の際に用いられる。
【００７６】
表示部１６は、モデル記憶部２０に記憶されたパラメータベクトル等の推定結果を表示する。地下構造推定装置１０を操作するユーザは、表示部１６によって表示された結果を確認する。
【００７７】
＜地下構造推定装置の作用＞
次に、図６を参照して、地下構造推定装置１０の作用を説明する。地下構造推定装置１０の受付部１２は、新たに地震が発生した場合、新たな地震の震源パラメータ及び各観測地点における各時刻の地震の速度データを受け付け、観測データ記憶部１８へ格納する。そして、地下構造推定装置１０は、図６に示す地下構造推定処理ルーチンを実行する。地下構造推定処理ルーチンは、新たな地震の観測データが観測データ記憶部１８へ格納される毎に実行される。
【００７８】
ステップＳ１００において、情報取得部２２は、観測データ記憶部１８に格納された観測データの震源パラメータを取得する。
【００７９】
ステップＳ１０２において、情報取得部２２は、震源パラメータに含まれる震源位置の情報に基づいて、観測データが地下構造モデルの領域内であるか否かを判定する。観測データが地下構造モデルの領域内である場合には、ステップＳ１０４へ進む。一方、観測データが地下構造モデルの領域内でない場合には、地下構造推定処理ルーチンを終了する。
【００８０】
ステップＳ１０４において、情報取得部２２は、観測データ記憶部１８に格納された観測データの地震の速度データを取得する。
【００８１】
ステップＳ１０６において、情報取得部２２は、上記ステップＳ１０４で取得された地震の速度データのＳＮ比が所定の条件を満たすか否かを判定する。地震の速度データのＳＮ比が所定の条件を満たす場合には、ステップＳ１０８へ進む。一方、地震の速度データのＳＮ比が所定の条件を満たさない場合には、地下構造推定処理ルーチンを終了する。
【００８２】
ステップＳ１０８において、情報取得部２２は、モデル記憶部２０に格納された地下構造モデルのパラメータを表すパラメータベクトルを取得する。
【００８３】
ステップＳ１１０において、推定部２４は、上記ステップＳ１０８で取得された地下構造モデルのパラメータベクトルθ及びシミュレーションに応じて予測される地震データと、上記ステップＳ１０４で取得された地震の速度データｙｔとの差分に応じて、地下構造モデルのパラメータベクトルθを推定する。ステップＳ１１０の処理は、図７に示すデータ同化処理ルーチンによって実現される。
【００８４】
＜データ同化処理ルーチン＞
ステップＳ２００において、時刻ｔに０を設定する。
【００８５】
ステップＳ２０２において、予測部２４０は、上記式（６）の条件付き確率分布に従って、拡大状態ベクトルの初期値ｘ０を設定する。
【００８６】
ステップＳ２０４において、予測部２４０は、格子点間の距離が近いほど、一方の格子点の地下方向の層厚と他方の格子点の地下方向の層厚とが近い値になるように、上記式（１２）の分散共分散行列Ｖ＾ｔ｜ｔ−１の層厚Ｈに関する項に対し、格子点の位置に応じた相関を設定する。
【００８７】
ステップＳ２０６において、予測部２４０は、前時刻ｔ−１の各アンサンブルｘ（ｉ）ｔ−１｜ｔ−１とシミュレーションｆｔとに基づいて、上記式（１０）に従って、現時刻ｔの拡大状態ベクトルを表す各アンサンブルｘ（ｉ）ｔ｜ｔ−１を予測する。
【００８８】
ステップＳ２０８において、修正部２４２は、上記ステップＳ２０６で予測された時刻ｔの各アンサンブルｘ（ｉ）ｔ｜ｔ−１から生成される地震データと、時刻ｔの地震の速度データｙｔとに基づいて、上記式（１５）に従って、時刻ｔの各アンサンブルｘ（ｉ）ｔ｜ｔ−１を修正する。そして、修正部２４２は、拡大状態ベクトルを表す各アンサンブルｘ（ｉ）ｔ｜ｔを得る。
【００８９】
ステップＳ２１０において、修正部２４２は、地震の速度データｙｔの終端時刻Ｔまで、上記ステップＳ２０６の一期先予測及び上記ステップＳ２０８のフィルタリングが繰り返されたか否かを判定する。終端時刻Ｔまで、一期先予測及びフィルタリングが繰り返された場合には、ステップＳ２１４へ進む。一方、終端時刻Ｔまで、一期先予測及びフィルタリングが繰り返されていない場合には、ステップＳ２１２で時刻を１インクリメントしてステップＳ２０６へ戻る。
【００９０】
ステップＳ２１４において、平滑化部２４４は、上記ステップＳ２０８で各時刻において得られた推定結果に含まれる、地下構造モデルのパラメータベクトルθを平滑化することにより、パラメータベクトルθを取得する。
【００９１】
ステップＳ２１６において、平滑化部２４４は、上記ステップＳ２１４で得られたパラメータベクトルθを結果として出力する。
【００９２】
次に、図６のステップＳ１１２において、平滑化部２４４は、上記ステップＳ１１０で出力されたパラメータベクトルθを、モデル記憶部２０へ格納する。モデル記憶部２０に格納されたパラメータベクトルθは、次回の地震の観測データが得られたときのデータ同化の際に用いられる。
【００９３】
ステップＳ１１４において、表示部１６は、モデル記憶部２０に記憶されたパラメータベクトルθ等の推定結果を表示して、地下構造推定処理ルーチンを終了する。
【００９４】
以上詳細に説明したように、本実施形態では、地下構造モデルに対応する地点における地震の観測データの時系列を入力とし、地下構造モデルのパラメータ及び地震予測モデルから予測される地震データと、地震の観測データとの差分に応じて、地震の状態及び地下構造モデルのパラメータを修正することを繰り返すことにより、地震の観測データを用いて地下構造を推定することができる。
【００９５】
地下構造モデルのパラメータ推定に対してデータ同化を適用しようとする場合、例えば、逐次フィルタリング手法の一例であるアンサンブルカルマンフィルタを用いるときには、拡大状態ベクトルの候補を表すアンサンブルが複数生成される。しかし、地下構造モデルのパラメータ数は多いため、パラメータに応じて生成されるアンサンブルに対する計算量も増大してしまうという課題があった。また、各アンサンブルは確率分布に応じてランダムに生成されるため、連続性を有する地下構造モデルのパラメータを推定したとしても、精度よく推定することができないという課題があった。
【００９６】
そこで、本実施形態では、データ同化により地下構造を推定する際に、地下構造モデルの各層に対してパラメータを付与することにより、パラメータ数を減少させる。これにより、データ同化により地下構造を推定する際に、計算量を低減させることができる。
【００９７】
また、本実施形態では、地下構造モデルのパラメータに対して相関を設定する。これにより、地下構造の連続性が考慮され、地下構造を精度よく推定することができる。
【００９８】
また、地震の観測データと整合度の高い地下構造モデルが推定されるため、地震動の予測精度を向上させることができる。また、地下構造モデルの修正のための追加のボーリング調査などが不要となる。
【００９９】
また、本実施形態では逐次データ同化を用いるため、地震の観測データが増加しても、増加分の観測データのみを用いて地下構造モデルのパラメータ修正を行うだけで良く、非逐次データ同化と比較して計算量が少なくて済む。
【０１００】
＜実験例＞
本実施形態の効果の確認のため、２次元ＳＨ波動場における数値実験を行った。なお、逐次データ同化手法の一つであるアンサンブルカルマンフィルタを用いて解析を行った。
【０１０１】
［実験の概要］
図８に示す２次元地下構造モデルを対象として数値実験を行った。まず、真のモデル（真値）に対してシミュレーションを行い、地表観測点（図８中の三角）の模擬地震動を作成する。また、以下の表１に示すように、真の地下構造モデルに１割の誤差を与えた初期モデル（初期値）を設定する。次に初期モデルに対するシミュレーションと模擬地震動とのデータ同化を実施し、初期モデルが真のモデルに近づくことを確認した。
【０１０２】
【表１】
2018128929.tif 000017
【０１０３】
［層に相関を設定した場合と層に相関を設定しない場合との比較］
【０１０４】
層に対して相関を与えた場合の数値実験の結果を図９（Ａ）に示す。図の横軸は時刻ステップ、縦軸は真値に対する推定値の比を示しており、比が１に近づくほど真値に近づくことを表す。層厚ＨとＳ波のせん断波速度ＶＳとついて、逐次修正され真値に近づいており、データ同化を用いた地下構造推定の効果が確認されている。
【０１０５】
一方、層に対して相関を設定しない場合の数値実験の結果を図９（Ｂ）に示す。相関が設定された場合と比較すると、相関を設定しない場合の層厚Ｈの推定精度が悪い。また、Ｓ波のせん断波速度ＶＳについても相関が設定された場合の方が推定精度は良く、層厚に相関を与えることの有効性が確認された。
【０１０６】
なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１０７】
例えば、上記実施形態では、各層の物性値のパラメータ（ρ，ＶＳ，ＶＰ，ＱＳ，ＱＰ）と、層厚のパラメータＨとをデータ同化により同時に推定する場合を例に説明したがこれに限定されるものではない。例えば、各パラメータのうちの所定のパラメータを推定した後、所定のパラメータ以外のパラメータを推定するようにしてもよい。例えば、各層の層厚Ｈと速度ＶＳとを推定した後に、密度ρと減衰定数ＱＳ及びＱＰを推定するようにしてもよい。
【０１０８】
また、上記実施形態では、図４に示されるように格子点間の水平距離（ｘｙ方向）に応じて相関を設定する場合を例に説明したが、これに限定されるものではなく、ｚ方向に対して相関を設定するようにしてもよい。例えば、ｚ方向の位置が深いほど層の固さは硬くなるように相関を設定するようにしてもよい。これにより、例えばアンサンブルカルマンフィルタにおいて生成される各アンサンブルにおいて、特定の層よりも、特定の層の下に位置する層の方が軟らかいことを表すアンサンブルの生成等が抑制される。
【０１０９】
また、上記実施形態では、層に関するパラメータの一例として層厚Ｈに対して相関を設定する場合を例に説明したが、これに限定されるものではない。例えば、層厚Ｈ以外のパラメータに対して相関を設定してもよい。具体的には、例えば、地下構造モデルのｚ方向の位置が深いほど速度ＶＳが大きくなるように、相関を設定してもよい。また、例えば、速度ＶＳと密度ρとの関係を考慮して、速度ＶＳが大きいほど密度ρも大きくなるように相関を設定してもよい。また、地下構造モデルのｚ方向の位置が深いほど減衰係数Ｑが大きくなるように、相関を設定してもよい。また、速度ＶＳと速度ＶＰとの間の関係を考慮して、相関を設定してもよい。また、実際の地盤調査結果、各パラメータ間の関係式、及び既往の研究の少なくとも１つに応じて、相関を設定してもよい。
【０１１０】
また、上記実施形態では、地下構造モデルのパラメータとして、各層の物性値（ρ，ＶＳ，ＶＰ，ＱＳ，ＱＰ）と層厚Ｈとをデータ同化により推定する場合を例に説明したがこれに限定されるものではない。本実施形態の地下構造モデルのパラメータの一部のみを推定してもよく、例えば、層厚Ｈ及び速度ＶＳのみを推定してもよい。また、層厚Ｈ及び速度ＶＳのみを推定した後に、減衰定数ＱＳを推定してもよい。
【０１１１】
また、上記実施形態では、地下構造推定処理ルーチンにおいて、新たな観測データに基づくデータ同化が行われる毎に、パラメータベクトルに対して相関を設定する場合を例に説明したが、これに限定されるものではない。例えば、はじめて得られた地震の観測データに基づきデータ同化を行う際にのみ、パラメータベクトルに対して相関を設定してもよい。
【０１１２】
また、上記実施形態では、データ同化の逐次フィルタリング手法の一例としてアンサンブルカルマンフィルタを用いる場合を例に説明したがこれに限定されるものではなく、他の逐次フィルタリング手法を用いてもよい。例えば、逐次フィルタリング手法として粒子フィルタを用いてもよい。
【０１１３】
また、上記ではプログラムが記憶部（図示省略）に予め記憶（インストール）されている態様を説明したが、プログラムは、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ及びマイクロＳＤカード等の記録媒体の何れかに記録されている形態で提供することも可能である。
【符号の説明】
【０１１４】
１０   地下構造推定装置
１２   受付部
１４   コンピュータ
１６   表示部
１８   観測データ記憶部
２０   モデル記憶部
２２   情報取得部
２４   推定部
２６   地下構造モデル
２４０ 予測部
２４２ 修正部
２４４ 平滑化部
θ     パラメータベクトル
【図１】
2018128929.tif 000018
【図２】
2018128929.tif 000019
【図３】
2018128929.tif 000020
【図４】
2018128929.tif 000021
【図５】
2018128929.tif 000022
【図６】
2018128929.tif 000023
【図７】
2018128929.tif 000024
【図８】
2018128929.tif 000025
【図９】
2018128929.tif 000026