14:15 〜 14:30
[S01-10] アジョイント方程式に基づく、2次元P-SV波動場及び震源推定の試み
１. はじめに
アジョイント方程式に基づき、地震波の観測波形から初期波動場すなわち震源の推定に向けた2次元P-SV問題での数値実験を行った。近年、地震動の予測や震源推定のために、数値モデルに基づいた推定値を地震計記録などの観測値を用いて修正しながら推定精度を高めるデータ同化の活用が進められている（例えばHoshiba & Aoki, 2015; Furumura et al., 2018）。これらの既往研究では最適内挿法を用いて、予測結果と観測の同化はタイムステップ毎に１回行われる。 本研究ではアジョイント法に基づくデータ同化により、予測を繰り返し更新するとともに、波動場の物理的性質を表す基礎方程式を同化に用いることで、波動場の推定精度を高めることを目指す。これまで森田・前田・髙野（2023）は、2次元SH波動場におけるアジョイント法の定式化とそれを用いた初期波動場推定の数値実験を通じて、その有効性を確認した。本研究ではこれをP-SV波動場へと拡張し、地震観測波形の逆伝播による震源推定への有効性を数値実験により検討した。
２. 手法・データ
数値実験は4層（堆積層、上部地殻、下部地殻、マントル）からなる速度構造をもつ400 kmｘ100 kmの2次元領域で行った（図1）。モデルの左よりの位置（100 km、30 km）にパルス幅4秒の点震源を与え、地表にはHi-netと同程度の20 kmの間隔で20点の観測点を置いた。地震波動場の計算はスタガードグリッド差分法（空間４次、時間２次精度）を用いて行い、タイムステップは0.0125 s、空間グリッドサイズは0.2 kmとした。 初めに、発震時から75秒間の地震波の伝播の計算（順計算）により、各観測点での観測波形を求め、これを正解データとして以降の実験に用いた。次に逆計算では、各観測点から時刻t=0まで地震波を逆伝播させ初期波源を推定した。この際、初期波源の変位速度は逆伝播の速度相当値をそのまま用いる一方で、応力には逆伝播の応力相当値に弾性コンプライアンステンソルを乗じたものを用いる必要がある。その後、推定された初期波源から順計算により得た各観測点波形の予測と正解の残差を求め、以降、残差の逆伝播と初期波源の推定及び順計算をイタレーションし残差の収束性を確認した。そして、震源深さや用いる速度構造を変えて、アジョイント法による初期波動場（震源波動場）推定の有効性を検討した。
３. 結果
順計算とアジョイント法による残差の逆計算を1500回イタレーションした結果を図2に示す。変位速度の水平・鉛直成分の二乗和から震動エネルギーを求め、それが震源付近で最大になる時刻（発震時に相当）のエネルギー分布を描画した。また、この初期波動場による各観測点の予測地震波形と正解データを比較した。イタレーション１回（アジョイント法無し）の逆伝播計算では、震源波動場が地表から地中へぼんやりと広がり（図 2 - a ）、この初期波動場から計算された予測波形もイタレーション回数が少ない内では観測波形との違いが大きい（図 2 - c）。しかし十分なイタレーションを重ねると、観測点間隔程度の解像度で正解に近い震源波動場が得られ（図 2 - b）、それにより観測波形とほぼ一致する波形が予測できた（図 2 - d）。また、各観測点の推定残差の全時間の二乗和で定義される目的関数もイタレーションに伴い指数関数的に減少した（図 2 - e）。 次に、震源深さを30 kmから5 km、70 kmに変え推定結果の違いを確認した。深さ5 kmの場合は震動エネルギーのピークが震源位置と一致せず、やや左側に求まった（図 3 - a）。震源が浅いことで表面波が観測記録で卓越し、それが震源推定へ用いられたのが理由かもしれない。実体波が卓越する深い震源（70 km）では、推定された震源の深さ方向の解像度が高くなることがわかった（図 3 - b）。これは、震源が深いほど観測記録において直接波が卓越する範囲が広くなることが関係することを示唆する。 最後に実問題への適用に向け、異なる速度構造モデルを用いた場合の初期波動場推定への影響を調べた。まず、図１の４層モデルとは異なるJMA2001（上野・他, 2002）を用いた場合でも、ほぼ同様の推定ができることが確認できた（図 4 - a）。次に、領域右側に最大深さ5 kmの低速度領域（Vp = 4 km/s、Vs = 2 km/s）を与えた場合には、震動エネルギーのピークが盆地と反対方向にずれ、分布も広がる結果となった（図 4 - b）。よって、水平方向の速度構造の理解が重要であることがわかる。
４. まとめと今後の課題
アジョイント法による計算と観測データの同化に基づき、２次元P-SV波の初期波動場の推定を行い、その有効性を数値実験により確認した。結果、一般的な強震観測点間隔（20 km程度）でも高い解像度で初期波動場（震源）が推定できた。しかし、実問題への適用では、地下構造モデルの解像度に依存するため、既知の速度構造の分解能で推定できる波動場の最大周波数について検討が必要である。
５. 謝辞
数値計算には、東大情報基盤センターのWisteria及び地震研究所のEIC計算機を使用した。
アジョイント方程式に基づき、地震波の観測波形から初期波動場すなわち震源の推定に向けた2次元P-SV問題での数値実験を行った。近年、地震動の予測や震源推定のために、数値モデルに基づいた推定値を地震計記録などの観測値を用いて修正しながら推定精度を高めるデータ同化の活用が進められている（例えばHoshiba & Aoki, 2015; Furumura et al., 2018）。これらの既往研究では最適内挿法を用いて、予測結果と観測の同化はタイムステップ毎に１回行われる。 本研究ではアジョイント法に基づくデータ同化により、予測を繰り返し更新するとともに、波動場の物理的性質を表す基礎方程式を同化に用いることで、波動場の推定精度を高めることを目指す。これまで森田・前田・髙野（2023）は、2次元SH波動場におけるアジョイント法の定式化とそれを用いた初期波動場推定の数値実験を通じて、その有効性を確認した。本研究ではこれをP-SV波動場へと拡張し、地震観測波形の逆伝播による震源推定への有効性を数値実験により検討した。
２. 手法・データ
数値実験は4層（堆積層、上部地殻、下部地殻、マントル）からなる速度構造をもつ400 kmｘ100 kmの2次元領域で行った（図1）。モデルの左よりの位置（100 km、30 km）にパルス幅4秒の点震源を与え、地表にはHi-netと同程度の20 kmの間隔で20点の観測点を置いた。地震波動場の計算はスタガードグリッド差分法（空間４次、時間２次精度）を用いて行い、タイムステップは0.0125 s、空間グリッドサイズは0.2 kmとした。 初めに、発震時から75秒間の地震波の伝播の計算（順計算）により、各観測点での観測波形を求め、これを正解データとして以降の実験に用いた。次に逆計算では、各観測点から時刻t=0まで地震波を逆伝播させ初期波源を推定した。この際、初期波源の変位速度は逆伝播の速度相当値をそのまま用いる一方で、応力には逆伝播の応力相当値に弾性コンプライアンステンソルを乗じたものを用いる必要がある。その後、推定された初期波源から順計算により得た各観測点波形の予測と正解の残差を求め、以降、残差の逆伝播と初期波源の推定及び順計算をイタレーションし残差の収束性を確認した。そして、震源深さや用いる速度構造を変えて、アジョイント法による初期波動場（震源波動場）推定の有効性を検討した。
３. 結果
順計算とアジョイント法による残差の逆計算を1500回イタレーションした結果を図2に示す。変位速度の水平・鉛直成分の二乗和から震動エネルギーを求め、それが震源付近で最大になる時刻（発震時に相当）のエネルギー分布を描画した。また、この初期波動場による各観測点の予測地震波形と正解データを比較した。イタレーション１回（アジョイント法無し）の逆伝播計算では、震源波動場が地表から地中へぼんやりと広がり（図 2 - a ）、この初期波動場から計算された予測波形もイタレーション回数が少ない内では観測波形との違いが大きい（図 2 - c）。しかし十分なイタレーションを重ねると、観測点間隔程度の解像度で正解に近い震源波動場が得られ（図 2 - b）、それにより観測波形とほぼ一致する波形が予測できた（図 2 - d）。また、各観測点の推定残差の全時間の二乗和で定義される目的関数もイタレーションに伴い指数関数的に減少した（図 2 - e）。 次に、震源深さを30 kmから5 km、70 kmに変え推定結果の違いを確認した。深さ5 kmの場合は震動エネルギーのピークが震源位置と一致せず、やや左側に求まった（図 3 - a）。震源が浅いことで表面波が観測記録で卓越し、それが震源推定へ用いられたのが理由かもしれない。実体波が卓越する深い震源（70 km）では、推定された震源の深さ方向の解像度が高くなることがわかった（図 3 - b）。これは、震源が深いほど観測記録において直接波が卓越する範囲が広くなることが関係することを示唆する。 最後に実問題への適用に向け、異なる速度構造モデルを用いた場合の初期波動場推定への影響を調べた。まず、図１の４層モデルとは異なるJMA2001（上野・他, 2002）を用いた場合でも、ほぼ同様の推定ができることが確認できた（図 4 - a）。次に、領域右側に最大深さ5 kmの低速度領域（Vp = 4 km/s、Vs = 2 km/s）を与えた場合には、震動エネルギーのピークが盆地と反対方向にずれ、分布も広がる結果となった（図 4 - b）。よって、水平方向の速度構造の理解が重要であることがわかる。
４. まとめと今後の課題
アジョイント法による計算と観測データの同化に基づき、２次元P-SV波の初期波動場の推定を行い、その有効性を数値実験により確認した。結果、一般的な強震観測点間隔（20 km程度）でも高い解像度で初期波動場（震源）が推定できた。しかし、実問題への適用では、地下構造モデルの解像度に依存するため、既知の速度構造の分解能で推定できる波動場の最大周波数について検討が必要である。
５. 謝辞
数値計算には、東大情報基盤センターのWisteria及び地震研究所のEIC計算機を使用した。