09:30 〜 09:45
[SSS07-03] アジョイント方程式に基づく地震波動場のデータ同化と震源メカニズム推定
キーワード：地震学、 データ同化、アジョイント法
1. はじめに
近年の稠密地震観測の環境下において、観測波形を時間を遡って震源に向けて逆伝播させ、発震時の初期波動場を推定する新たな震源推定手法が提案されている（McMechan, 1982; Larmat et al, 2006）。従来の震源推定の課題であるP／S波の読み取り精度に依存する震源推定誤差や、立ち上がりが不明瞭な微動への対応、また同時発生した複数地震の分離などの問題解決に有効なアプローチと考えられる。
本研究では、地震波逆伝播計算にデータ同化手法の一つであるアジョイント法を導入することで、震源のエネルギーを含めた高精度な震源推定に挑戦する。これまで、最適内挿法に基づく地震波伝播シミュレーションと観測値のデータ同化により、未来の時刻の揺れや津波の伝播を予測する研究が行なわれている（Hoshiba & Aoki, 2015; Wang et al., 2017；Furumura et al., 2018）。一方、本研究で用いるアジョイント法によるデータ同化では、順伝播計算で得られたシミュレーション波形と観測波形との残差を、基礎方程式を用いて逆伝播させて波動場を修正する過程を繰り返し実行して実波動場に近づける点が異なる。計算量は多くなるが、反復計算を行わない同化手法よりも高い精度と絶対振幅の評価が期待できる。
データ同化におけるアジョイント法の有効性は、これまで２次元音響場（森田・前田・髙野, 2023）や２次元弾性場（森田・古村・前田, 2023）を対象とした数値実験から確認した。これらの成果に基づき、本研究では稠密地震観測の実波形データを用いて震源の位置とメカニズムの推定を試みた。
２. 手法・データ
実験では、MeSO-netの400点の観測点のうち、山梨県道志村から茨城県潮来市にかけて関東平野を直線上に横切る31点の地震観測点（Fig1）での波形データを使用した。観測点間隔は約5 kmである。観測波形には0.05Hz〜0.5Hzのバンドパスフィルタをかけ、速度波形に変換したものを用いた。また、地震波の逆伝播・順伝播計算は、MeSO-netの側線に沿った水平200 km、鉛直100 kmの2次元領域で行ない、速度構造には全国１次地下構造モデル（JIVSM; Koketsu et al., 2012）を使用した（Fig２）。なお、計算の安定化のために、最小S波速度は2.0 km/sに設定し、浅部の堆積層は省略した。地震波計算はスタガード格子差分法（空間４次精度・時間２次精度）により行い、時間ステップを0.01 s、格子サイズを0.2 kmとした。
解析の流れとして、P波、S波を含む40秒間の観測波形を各観測点から時間を遡って震源の方向に向けて逆伝播させ、発震時における最初の波源を求める。次に、これを初期条件として、地震波動場を差分法計算により順伝播させて、各観測点での計算波形を得る。そして、観測と計算波形の残差を求め、これを再び逆伝播させて初期波動場を修正する操作（アジョイント計算）を、残差が収束するまで繰り返す。そして、最終的に得られた初期波動場に対して、地震波がフォーカスした地点とその時刻を、震源及び発震時とする。ここで、地震波のフォーカスの評価（イメージング）は、各格子点での弾性エネルギー（ひずみと応力の各成分の積和）を用いて行った。
３. 結果
2021年12月29日の東京都23区直下の地震（深さ31 km、Mj 3.5）観測データを用いてアジョイント計算を行った。計算は、東京大学情報基盤センターのWisteriaスパコンの１CPUによる32コアを用いたOpenMP並列計算により行ない、2000回のイタレーションに2時間を要した。
アジョイント計算を繰り返し、最終的に得られた計算波形は観測波形をよく説明し（VR＝85 %; Fig 3）、発震時に近い時刻（t = 0.26秒）に気象庁やF-netの推定震源に近い場所にエネルギーのピークが得られ、適切に震源が推定されたと判断できる（Fig 4）。次に、差分法計算により得られた各格子点の応力値を用いて震源メカニズム （CMT）を計算したところ、傾斜角は60度ほど異なるが、F-net解と同じ逆断層型の解が得られた（Fig 5）。
一方、観測波形をより長時間（120秒）までS波の後続相を含め、アジョイント計算を行ったところ、震源においてより大きな弾性エネルギーのフォーカシングが見られたが、VRは50 %まで大幅に低下した。これは、後続相には関東平野の３次元堆積構造の影響を強く受けており、計算を行った２次元モデルの面外方向から到来する波が多く含まれるためと考えられ、この対処には３次元計算が必要である。
次に、観測点間隔をK-NET／KiK-net、Hi-netと同程度の約20 kmと粗くした場合についても検討した。この結果、弾性エネルギーが観測点近傍の複数地点に集中する不適切な結果となったが、VRは98 %と高い値になった。これは、各観測点での波形記録を説明する特殊な解として、それぞれの観測点近傍に複数の震源があり、これらから各観測点の波形に対応した地震波を放射されたものと推定されたためである。
４. まとめと今後の課題
本研究ではアジョイント法によるデータ同化を用いて、稠密地震観測による実波形データによる関東直下の地震の初期波動場（震源）推定を行ない、気象庁やF-net解に近い震源メカニズムが得られ、計算と観測波形も高い一致度（VR）を示した。さらに、多様な場所で起きる地震への対応や、震源推定の高精度化には３次元計算が必要となるが、仮に２次元モデルを面外方向に100 km広げた場合には計算量が800倍に増大する。この規模は、MPIを用いた並列計算により不可能ではないが、実現に向け、アジョイント計算におけるイタレーション数の低減化や、精度の高いイメージング手法の導入などにより、より計算効率を高めることが先決である。
５. 謝辞
数値計算には東大情報基盤センターのWisteriaを使用した。また、観測データは防災科学技術研究所のMeSO-net観測波形を使用した。
近年の稠密地震観測の環境下において、観測波形を時間を遡って震源に向けて逆伝播させ、発震時の初期波動場を推定する新たな震源推定手法が提案されている（McMechan, 1982; Larmat et al, 2006）。従来の震源推定の課題であるP／S波の読み取り精度に依存する震源推定誤差や、立ち上がりが不明瞭な微動への対応、また同時発生した複数地震の分離などの問題解決に有効なアプローチと考えられる。
本研究では、地震波逆伝播計算にデータ同化手法の一つであるアジョイント法を導入することで、震源のエネルギーを含めた高精度な震源推定に挑戦する。これまで、最適内挿法に基づく地震波伝播シミュレーションと観測値のデータ同化により、未来の時刻の揺れや津波の伝播を予測する研究が行なわれている（Hoshiba & Aoki, 2015; Wang et al., 2017；Furumura et al., 2018）。一方、本研究で用いるアジョイント法によるデータ同化では、順伝播計算で得られたシミュレーション波形と観測波形との残差を、基礎方程式を用いて逆伝播させて波動場を修正する過程を繰り返し実行して実波動場に近づける点が異なる。計算量は多くなるが、反復計算を行わない同化手法よりも高い精度と絶対振幅の評価が期待できる。
データ同化におけるアジョイント法の有効性は、これまで２次元音響場（森田・前田・髙野, 2023）や２次元弾性場（森田・古村・前田, 2023）を対象とした数値実験から確認した。これらの成果に基づき、本研究では稠密地震観測の実波形データを用いて震源の位置とメカニズムの推定を試みた。
２. 手法・データ
実験では、MeSO-netの400点の観測点のうち、山梨県道志村から茨城県潮来市にかけて関東平野を直線上に横切る31点の地震観測点（Fig1）での波形データを使用した。観測点間隔は約5 kmである。観測波形には0.05Hz〜0.5Hzのバンドパスフィルタをかけ、速度波形に変換したものを用いた。また、地震波の逆伝播・順伝播計算は、MeSO-netの側線に沿った水平200 km、鉛直100 kmの2次元領域で行ない、速度構造には全国１次地下構造モデル（JIVSM; Koketsu et al., 2012）を使用した（Fig２）。なお、計算の安定化のために、最小S波速度は2.0 km/sに設定し、浅部の堆積層は省略した。地震波計算はスタガード格子差分法（空間４次精度・時間２次精度）により行い、時間ステップを0.01 s、格子サイズを0.2 kmとした。
解析の流れとして、P波、S波を含む40秒間の観測波形を各観測点から時間を遡って震源の方向に向けて逆伝播させ、発震時における最初の波源を求める。次に、これを初期条件として、地震波動場を差分法計算により順伝播させて、各観測点での計算波形を得る。そして、観測と計算波形の残差を求め、これを再び逆伝播させて初期波動場を修正する操作（アジョイント計算）を、残差が収束するまで繰り返す。そして、最終的に得られた初期波動場に対して、地震波がフォーカスした地点とその時刻を、震源及び発震時とする。ここで、地震波のフォーカスの評価（イメージング）は、各格子点での弾性エネルギー（ひずみと応力の各成分の積和）を用いて行った。
３. 結果
2021年12月29日の東京都23区直下の地震（深さ31 km、Mj 3.5）観測データを用いてアジョイント計算を行った。計算は、東京大学情報基盤センターのWisteriaスパコンの１CPUによる32コアを用いたOpenMP並列計算により行ない、2000回のイタレーションに2時間を要した。
アジョイント計算を繰り返し、最終的に得られた計算波形は観測波形をよく説明し（VR＝85 %; Fig 3）、発震時に近い時刻（t = 0.26秒）に気象庁やF-netの推定震源に近い場所にエネルギーのピークが得られ、適切に震源が推定されたと判断できる（Fig 4）。次に、差分法計算により得られた各格子点の応力値を用いて震源メカニズム （CMT）を計算したところ、傾斜角は60度ほど異なるが、F-net解と同じ逆断層型の解が得られた（Fig 5）。
一方、観測波形をより長時間（120秒）までS波の後続相を含め、アジョイント計算を行ったところ、震源においてより大きな弾性エネルギーのフォーカシングが見られたが、VRは50 %まで大幅に低下した。これは、後続相には関東平野の３次元堆積構造の影響を強く受けており、計算を行った２次元モデルの面外方向から到来する波が多く含まれるためと考えられ、この対処には３次元計算が必要である。
次に、観測点間隔をK-NET／KiK-net、Hi-netと同程度の約20 kmと粗くした場合についても検討した。この結果、弾性エネルギーが観測点近傍の複数地点に集中する不適切な結果となったが、VRは98 %と高い値になった。これは、各観測点での波形記録を説明する特殊な解として、それぞれの観測点近傍に複数の震源があり、これらから各観測点の波形に対応した地震波を放射されたものと推定されたためである。
４. まとめと今後の課題
本研究ではアジョイント法によるデータ同化を用いて、稠密地震観測による実波形データによる関東直下の地震の初期波動場（震源）推定を行ない、気象庁やF-net解に近い震源メカニズムが得られ、計算と観測波形も高い一致度（VR）を示した。さらに、多様な場所で起きる地震への対応や、震源推定の高精度化には３次元計算が必要となるが、仮に２次元モデルを面外方向に100 km広げた場合には計算量が800倍に増大する。この規模は、MPIを用いた並列計算により不可能ではないが、実現に向け、アジョイント計算におけるイタレーション数の低減化や、精度の高いイメージング手法の導入などにより、より計算効率を高めることが先決である。
５. 謝辞
数値計算には東大情報基盤センターのWisteriaを使用した。また、観測データは防災科学技術研究所のMeSO-net観測波形を使用した。