[S01P-06] Forecasting Seismic Wave Propagation by Physics-Informed Neural Network (PINN): Numerical Experiments in a 2-D Acoustic Wavefield
1．背景・目的
近年，機械学習を活用した数値モデリング手法である「Physics-Informed Neural Network（以下，PINN）」が注目されている（Raissi et al. 2019）．通常の機械学習とは異なり，PINNはデータだけでなく物理方程式も学習の拘束条件として用いることで，学習と予測がデータと物理法則の双方に従うことが保証され，高精度な予測の実現が期待できる．
本研究では，地震波伝播を学習したPINNの予測能力を検証するため，2次元の音響波動場を対象として数値実験を行った．差分法で計算した波形データを用いてPINNを訓練し，学習済みモデルを用いて任意の地点と時刻の波動場の再構築（予測）を試みた．
2．手法
40 km x 40 kmの2次元平面を数値実験領域に設定し，6秒間の音波伝播を評価する．まず，空間4次・時間2次精度のstaggered-grid差分法を用いて時間ステップと格子サイズを0.01 sと0.1 kmとする音響波動場を計算し，各gridでの速度波形データを合成した．音波速度は一様（3 km/s）とし，破壊継続時間0.5 sの震源から放射される音響波動場の卓越波長は1.5 kmである．
次に，PINNのPythonライブラリ：DeepXDE（Lu et al. 2021）を用いて，隠れ層4層，各ニューロン数50で，tanh関数を活性化関数とする全結合型ニューラルネットワークを作成し（図1），以下の損失関数Lを最小とするよう学習を行った．
L = w_data*L_data + w_PDE*L_PDE
ここで，L_dataはPINNによる予測と波形データの平均二乗誤差（データ拘束），L_PDEは音響波動場を表す音波方程式の平均二乗誤差（物理拘束），w_dataとw_PDEは各拘束の重みである．なお，optimizerにはAdam法を用いた．10万epochの学習に，東大情報基盤センターのWisteria/BDEC-01スパコン（8GPU並列）で5分程度を要した．
3．結果と考察
領域内のデータ拘束点（L_dataを計算する点）におけるt=0 sからt=6 sまでの速度波形を学習した後，学習済みモデルを用いて予測した音響波動場を数値シミュレーション結果と比較した．データ拘束点の数，物理拘束点（L_PDEを計算する点）の数，w_dataとw_PDEの比率を変え，予測された波動場の精度を調べた．予測の精度は，数値シミュレーション結果と予測結果の平均二乗誤差として評価した．また，データ拘束のみを損失関数の計算に用い，物理拘束を用いない学習（w_PDE=0，すなわち通常の機械学習）との比較も行った．
なお，実際の地震観測点の分布は不均一であるが，拘束点の偏りに起因する予測精度のバイアスを避けるため，本数値実験ではデータ・物理拘束点ともに規則的に配置し，両者の密度の違いによる波動場の予測精度を評価した．
まず，データ拘束点の密度を変えた実験では，データ拘束点の平均間隔が音波の卓越波長と同程度の場合は精度よく予測できた（図2，N_data = 81, 49）が，波長に比べて1桁大きくなると精度の悪化が認められた（N_data = 9）．
同様の傾向は物理拘束点の密度を変えた場合についても確認されたが，データ拘束点の数が十分多ければ，物理拘束点の数が少なくとも精度低下は小さいことがわかった．
データ拘束（w_data）と物理拘束（w_PDE）の比率を変えた実験（図3）では，物理拘束の重みを増やすに従って予測精度の向上が確認できた．ただし，物理拘束の重みを非常に大きく設定し，データ拘束の寄与を小さくすると精度は悪化し，2つの重みを等しく設定した場合に最も良い精度が得られることが確認できた（図3，w_data : w_PDE = 1 : 1）．
以上の結果から，予測精度を高める上で物理拘束を加えたPINNは有効であるが，その場合でも，予測にはある程度のデータ拘束が必要であることが示唆された．
4．今後の展望
今回行ったPINNの数値実験は，波形データを学習した時間内の波動場の再構築（内挿予測）であり，さらに未来の波動場の予測（外挿予測）の性能評価を行う予定である．また，PINNの先行研究（Rasht-Behesht et al. 2022）では，震源からの音波伝播が始まった直後の2つの時刻の波動場のスナップショットから，その後の時刻の波動場が予測できることが示されており，現実的な少ない観測点と限られた観測時間の下での，波動場の予測可能性とその要件について，さらなる数値実験を進める予定である．さらに，PINNではPDEの未知パラメータ推定が可能であるため，波動場の推定とともに速度構造インバージョンへの適用可能性も検討する．
近年，機械学習を活用した数値モデリング手法である「Physics-Informed Neural Network（以下，PINN）」が注目されている（Raissi et al. 2019）．通常の機械学習とは異なり，PINNはデータだけでなく物理方程式も学習の拘束条件として用いることで，学習と予測がデータと物理法則の双方に従うことが保証され，高精度な予測の実現が期待できる．
本研究では，地震波伝播を学習したPINNの予測能力を検証するため，2次元の音響波動場を対象として数値実験を行った．差分法で計算した波形データを用いてPINNを訓練し，学習済みモデルを用いて任意の地点と時刻の波動場の再構築（予測）を試みた．
2．手法
40 km x 40 kmの2次元平面を数値実験領域に設定し，6秒間の音波伝播を評価する．まず，空間4次・時間2次精度のstaggered-grid差分法を用いて時間ステップと格子サイズを0.01 sと0.1 kmとする音響波動場を計算し，各gridでの速度波形データを合成した．音波速度は一様（3 km/s）とし，破壊継続時間0.5 sの震源から放射される音響波動場の卓越波長は1.5 kmである．
次に，PINNのPythonライブラリ：DeepXDE（Lu et al. 2021）を用いて，隠れ層4層，各ニューロン数50で，tanh関数を活性化関数とする全結合型ニューラルネットワークを作成し（図1），以下の損失関数Lを最小とするよう学習を行った．
L = w_data*L_data + w_PDE*L_PDE
ここで，L_dataはPINNによる予測と波形データの平均二乗誤差（データ拘束），L_PDEは音響波動場を表す音波方程式の平均二乗誤差（物理拘束），w_dataとw_PDEは各拘束の重みである．なお，optimizerにはAdam法を用いた．10万epochの学習に，東大情報基盤センターのWisteria/BDEC-01スパコン（8GPU並列）で5分程度を要した．
3．結果と考察
領域内のデータ拘束点（L_dataを計算する点）におけるt=0 sからt=6 sまでの速度波形を学習した後，学習済みモデルを用いて予測した音響波動場を数値シミュレーション結果と比較した．データ拘束点の数，物理拘束点（L_PDEを計算する点）の数，w_dataとw_PDEの比率を変え，予測された波動場の精度を調べた．予測の精度は，数値シミュレーション結果と予測結果の平均二乗誤差として評価した．また，データ拘束のみを損失関数の計算に用い，物理拘束を用いない学習（w_PDE=0，すなわち通常の機械学習）との比較も行った．
なお，実際の地震観測点の分布は不均一であるが，拘束点の偏りに起因する予測精度のバイアスを避けるため，本数値実験ではデータ・物理拘束点ともに規則的に配置し，両者の密度の違いによる波動場の予測精度を評価した．
まず，データ拘束点の密度を変えた実験では，データ拘束点の平均間隔が音波の卓越波長と同程度の場合は精度よく予測できた（図2，N_data = 81, 49）が，波長に比べて1桁大きくなると精度の悪化が認められた（N_data = 9）．
同様の傾向は物理拘束点の密度を変えた場合についても確認されたが，データ拘束点の数が十分多ければ，物理拘束点の数が少なくとも精度低下は小さいことがわかった．
データ拘束（w_data）と物理拘束（w_PDE）の比率を変えた実験（図3）では，物理拘束の重みを増やすに従って予測精度の向上が確認できた．ただし，物理拘束の重みを非常に大きく設定し，データ拘束の寄与を小さくすると精度は悪化し，2つの重みを等しく設定した場合に最も良い精度が得られることが確認できた（図3，w_data : w_PDE = 1 : 1）．
以上の結果から，予測精度を高める上で物理拘束を加えたPINNは有効であるが，その場合でも，予測にはある程度のデータ拘束が必要であることが示唆された．
4．今後の展望
今回行ったPINNの数値実験は，波形データを学習した時間内の波動場の再構築（内挿予測）であり，さらに未来の波動場の予測（外挿予測）の性能評価を行う予定である．また，PINNの先行研究（Rasht-Behesht et al. 2022）では，震源からの音波伝播が始まった直後の2つの時刻の波動場のスナップショットから，その後の時刻の波動場が予測できることが示されており，現実的な少ない観測点と限られた観測時間の下での，波動場の予測可能性とその要件について，さらなる数値実験を進める予定である．さらに，PINNではPDEの未知パラメータ推定が可能であるため，波動場の推定とともに速度構造インバージョンへの適用可能性も検討する．