[S01P-10] Data assimilation of seismic wavefields and estimation of source mechanism based on adjoint equations: Application to MeSO-net data
1. はじめに
本研究では、波形逆伝播手法にデータ同化手法の一つであるアジョイント法を用いた新たな震源イメージング手法を開発した（森田・古村・前田、2023SSJ）。本手法では、観測波形のP、S走時だけでなく、観測波形全体を時間に遡って地中に戻し、震源イメージングを行うため、P、Sの読み取りが難しい波形でも高精度な推定が期待できる。従来の最適内挿法（Furumura & Maeda, 2021）による（１回の）データ同化に対し、アジョイント法は波動伝播計算を繰り返し実行して同化を繰り返し行うことで、高精度の震源イメージングとその絶対値の評価が可能である。本発表では、これまでの数値実験に基づく成果を踏まえ、MeSO-netによる首都圏の小規模地震（M3〜5）の実観測データを用いて震源イメージングを行い、本手法の有効性を検討した。
２. 手法・データ
解析では、はじめに発震時刻の10秒以上前から50秒間の観測波形を各観測点から地中に向けて逆伝播させ、発震時付近の波動場（初期波動場）を求める。次に、これを初期値として波動伝播の順計算を行い、各観測点での計算波形を得る。そして、計算波形と観測波形の残差を再び逆伝播させ初期波動場の推定を修正する操作を、残差が収束するまで繰り返す（Fig. 1）。
最終的に推定した初期波動場の順計算結果により得られる各時刻の波動場に対して、弾性エネルギー（ひずみと応力の各成分の積和）が最大となる地点と時刻をそれぞれ震源と発震時として震源イメージングを行う。震源イメージングでは、波動場の震幅、運動エネルギー、P波とS波エネルギーの相関係数にもとづき行われることが多いが（Gajewski & Tessmer, 2005; Artman et al., 2010等）が、本研究ではS/Nの低い波形データにおいて、弾性エネルギーに基づくイメージングが適していることを数値実験から確認した。アジョイント計算の対象領域をFig. 2に示す。本研究ではMeSO-net観測網のうち、山梨県東部から茨城県東部にかけて関東平野を直線上に横切る31点の観測点（間隔約5 km）の波形記録を使用した。観測波形は速度波形に変換後、バンドパスフィルタ（0.1〜1Hz）をかけて使用した。また、計算領域は観測側線に沿った水平200 km、鉛直100 kmの2次元領域とし、速度構造はJIVSM（Koketsu et al., 2012）を使用した。なお、JIVSMの最小S波速度（0.35 km/s）に対して、本研究では計算安定化のために最小速度を2.0 km/sに設定した。波動計算はスタガード格子差分法（空間４次・時間２次精度）により、時間ステップを0.01 s、格子サイズを0.2 kmとして行った。計算は、東京大学情報基盤センターのWisteria-aスパコンの1CPU（64スレッド）による並列計算により行い、1000回のイタレーションに約1時間半を要した。
３. 結果
本研究で対象とする地震をFig. 2に示す；① 2021/12/29東京都23区の地震（深さ31 km、Mj 3.5）、② 2016/2/5神奈川県東部の地震（深さ26 km、Mj4.6）、③ 2012/5/29千葉県北西部の地震（深さ64 km、Mj5.2）を対象に震源イメージングを行った。結果をFig. 3に示す。まず、①のフィリピン海プレート上面の地震では、最終的に得られた計算波形が観測波形とV.R=56 %で適合し、気象庁震源に近い位置と時刻にイメージングがされた。地震②についても同様の結果を得た。しかし、初期波動場の応力分布から推定され震源メカニズム（モーメント・テンソル）は、F-netと傾斜角は類似するが、滑り角の符号が反転した結果となった。こうしたメカニズム解の符号の反転は、地震波伝播における振幅の正・負震動に対応したものであり、震源メカニズムの推定への適用に向けさらなる検討が必要である。次に、太平洋プレート上面の地震③では、震源位置とは離れた場所 (特に地表)にイメージング結果が得られた。原因の一つとして、震源が深く伝播経路が長いため速度構造モデルの違いの影響の累積が考えられる。ところが、この場合でも得られた計算と観測波形の残差は小さく（V.R.=81%）、見かけ上は良く一致する結果となった。これは、アジョイント計算の過程で、それぞれの観測波形を説明するための「偽」の震源が観測点直下に推定されたためであり、V,R.を推定の正しさの指標とする際は注意が必要である。
４. まとめ・課題
本研究ではMeSO-net側線に沿った２次元断面でアジョイント計算を行なった。しかし、現実問題への適用には３次元化が必須である。本研究で用いた２次元モデルを面外方向に100 km広げると計算量は数百倍に増大する。そのため、計算コードの並列化とともに、適切な初期条件の設定等によるイタレーション数を削減など、現行の計算効率を１〜２桁高める必要がある。
５. 謝辞
数値計算には東大情報基盤センターのWisteriaを使用した。また、観測データは防災科学技術研究所のMeSO-net観測波形を使用した。
本研究では、波形逆伝播手法にデータ同化手法の一つであるアジョイント法を用いた新たな震源イメージング手法を開発した（森田・古村・前田、2023SSJ）。本手法では、観測波形のP、S走時だけでなく、観測波形全体を時間に遡って地中に戻し、震源イメージングを行うため、P、Sの読み取りが難しい波形でも高精度な推定が期待できる。従来の最適内挿法（Furumura & Maeda, 2021）による（１回の）データ同化に対し、アジョイント法は波動伝播計算を繰り返し実行して同化を繰り返し行うことで、高精度の震源イメージングとその絶対値の評価が可能である。本発表では、これまでの数値実験に基づく成果を踏まえ、MeSO-netによる首都圏の小規模地震（M3〜5）の実観測データを用いて震源イメージングを行い、本手法の有効性を検討した。
２. 手法・データ
解析では、はじめに発震時刻の10秒以上前から50秒間の観測波形を各観測点から地中に向けて逆伝播させ、発震時付近の波動場（初期波動場）を求める。次に、これを初期値として波動伝播の順計算を行い、各観測点での計算波形を得る。そして、計算波形と観測波形の残差を再び逆伝播させ初期波動場の推定を修正する操作を、残差が収束するまで繰り返す（Fig. 1）。
最終的に推定した初期波動場の順計算結果により得られる各時刻の波動場に対して、弾性エネルギー（ひずみと応力の各成分の積和）が最大となる地点と時刻をそれぞれ震源と発震時として震源イメージングを行う。震源イメージングでは、波動場の震幅、運動エネルギー、P波とS波エネルギーの相関係数にもとづき行われることが多いが（Gajewski & Tessmer, 2005; Artman et al., 2010等）が、本研究ではS/Nの低い波形データにおいて、弾性エネルギーに基づくイメージングが適していることを数値実験から確認した。アジョイント計算の対象領域をFig. 2に示す。本研究ではMeSO-net観測網のうち、山梨県東部から茨城県東部にかけて関東平野を直線上に横切る31点の観測点（間隔約5 km）の波形記録を使用した。観測波形は速度波形に変換後、バンドパスフィルタ（0.1〜1Hz）をかけて使用した。また、計算領域は観測側線に沿った水平200 km、鉛直100 kmの2次元領域とし、速度構造はJIVSM（Koketsu et al., 2012）を使用した。なお、JIVSMの最小S波速度（0.35 km/s）に対して、本研究では計算安定化のために最小速度を2.0 km/sに設定した。波動計算はスタガード格子差分法（空間４次・時間２次精度）により、時間ステップを0.01 s、格子サイズを0.2 kmとして行った。計算は、東京大学情報基盤センターのWisteria-aスパコンの1CPU（64スレッド）による並列計算により行い、1000回のイタレーションに約1時間半を要した。
３. 結果
本研究で対象とする地震をFig. 2に示す；① 2021/12/29東京都23区の地震（深さ31 km、Mj 3.5）、② 2016/2/5神奈川県東部の地震（深さ26 km、Mj4.6）、③ 2012/5/29千葉県北西部の地震（深さ64 km、Mj5.2）を対象に震源イメージングを行った。結果をFig. 3に示す。まず、①のフィリピン海プレート上面の地震では、最終的に得られた計算波形が観測波形とV.R=56 %で適合し、気象庁震源に近い位置と時刻にイメージングがされた。地震②についても同様の結果を得た。しかし、初期波動場の応力分布から推定され震源メカニズム（モーメント・テンソル）は、F-netと傾斜角は類似するが、滑り角の符号が反転した結果となった。こうしたメカニズム解の符号の反転は、地震波伝播における振幅の正・負震動に対応したものであり、震源メカニズムの推定への適用に向けさらなる検討が必要である。次に、太平洋プレート上面の地震③では、震源位置とは離れた場所 (特に地表)にイメージング結果が得られた。原因の一つとして、震源が深く伝播経路が長いため速度構造モデルの違いの影響の累積が考えられる。ところが、この場合でも得られた計算と観測波形の残差は小さく（V.R.=81%）、見かけ上は良く一致する結果となった。これは、アジョイント計算の過程で、それぞれの観測波形を説明するための「偽」の震源が観測点直下に推定されたためであり、V,R.を推定の正しさの指標とする際は注意が必要である。
４. まとめ・課題
本研究ではMeSO-net側線に沿った２次元断面でアジョイント計算を行なった。しかし、現実問題への適用には３次元化が必須である。本研究で用いた２次元モデルを面外方向に100 km広げると計算量は数百倍に増大する。そのため、計算コードの並列化とともに、適切な初期条件の設定等によるイタレーション数を削減など、現行の計算効率を１〜２桁高める必要がある。
５. 謝辞
数値計算には東大情報基盤センターのWisteriaを使用した。また、観測データは防災科学技術研究所のMeSO-net観測波形を使用した。