[SSS14-P14] Reversible-jump MCMCによる平滑化条件を取り入れない断層すべり分布の推定
キーワード：測地インバージョン、平滑化拘束、Reversible-jump MCMC、Trans-dimensional inversion
1. 導入
地殻変動観測データを用いた断層すべり分布の推定には，すべりの平滑化条件を課した最小二乗法が行われ，その度合いを調整するハイパーパラメータをABICやLカーブ等の基準により決定する．しかし，地震時・地震後すべり分布を同時に推定する粘弾性インバージョン [Tomita et al., 2018] では，平滑化度合いを調整するハイパーパラメータ・観測データの重みに対するハイパーパラメータが時間窓の数だけ増大し，計算コストが極めて高くなってしまう．また，解の非負拘束を課した最小二乗法を行う場合，非負拘束を課した場合にはABICが解析的に正しく計算できないこと [Fukuda & Johnson, 2008] に加え，計算コストが通常の最小二乗法に比べて劇的に増大する問題がある．
上記の問題の解決のため，我々は，Reversible-jump MCMC法（Rj-MCMC [Green, 1995]）による粘弾性インバージョン解析の開発に取り組んでいる．本手法では，解の非負拘束を課した場合でも事後確率分布に基づき適切に解を求めることが可能であり，未知パラメータの数も同時に探索することで，最小限の未知パラメータ数で計算可能である．富田・他 [2019，地震学会] では，Simuated Annealingを用いたRj-MCMC手法を粘弾性問題を仮定したSyntheticデータに対して適用し，通常のABICを用いた手法と同様にすべり分布を求められることを示した．しかし，解の探索の効率性と十分性に問題が残っていた．そこで，本研究では効率的な解の探索が可能なParallel Tempering (PT) [e.g., Sambridge, 2013]を用いたRj-MCMCの開発を試みた．
2. 手法・データ
本研究では，Bodin & Sambridge [2009]と同様にVoronoi cellを用いた空間分割と定式化を行った．この空間分割により平滑化のハイパーパラメータは不要となる．PTでは，複数の異なる温度状態のレプリカを作成し，解の探索の過程でその温度状態を入れ替えることで局所解を避け，効率的な解の探索を行う．計算速度を上げるため，MPIによってレプリカを並列させて計算を行った．
Syntheticデータは，弾性媒質 [Okada, 1992] を仮定した地震時すべりに対する応答として求めた．PTを実装したRj-MCMCの解の収束状態・断層すべり分布の再現性の確認するため，粘弾性媒質ではなく，１時間窓の簡易な弾性媒質を仮定している．断層形状は，沈み込み帯プレート境界断層を模した傾斜角15度・幅300 km・長さ500 kmの矩形断層面を仮定し，20 km × 20 kmの小断層を配置した．地震時すべり分布としては，海溝から75 kmの地点をピークとした最大すべり量 8 mの滑らかな滑りを与えた。観測点配置は，海溝から200 kmの距離を海岸線として，①150陸上観測点＋5海底観測点，②150陸上観測点＋150海底観測点，の2パターンを設定した．
3. 結果・議論
上記のSyntheticデータに対し，PTのレプリカ数・最大温度・サンプル数を変えてRj-MCMCインバージョンを実行した．レプリカ数・最大温度は4〜16の間，サンプル数は105〜107の間で検証を行った．結果として，サンプル数が106以上であれば，最終的なすべり分布はレプリカ数・最大温度によらず仮定したすべりを概ね再現できていることを確認した．レプリカ数が多い場合は，すべり分布の収束が早まる傾向が見られた．なお，未知パラメータ数は，観測点配置①では10個程度，②では20個程度に縮約された．これは，観測点分布による解像度の違いを反映していると考えられる．一方で，各小断層でのパラメータの確率密度関数は，観測点配置①では非常に滑らかに得られたのに対し，②では複数の急峻なピークが見られた．観測点配置②では観測データが多いため，特定の観測データへのフィッティングが良くなる場合などの局所解が現れやすくなるためだと考えられる．このような場合でも，複数のレプリカのアンサンブルによって最終的なすべり分布は仮定したすべり分布を偏りなく再現できていた．
本研究では，PTの実装により解の収束性を高め，偏りなくすべり分布を推定可能であることを確かめた．今後はレプリカ数・最大温度のパターンを更に増やすことで，観測点配置②の場合でも局所解に陥りにくくすることが可能かを検証する．また，仮定するすべり分布のパターンを増やし，通常のABICインバージョンの結果との再現性の違いを検証する．
地殻変動観測データを用いた断層すべり分布の推定には，すべりの平滑化条件を課した最小二乗法が行われ，その度合いを調整するハイパーパラメータをABICやLカーブ等の基準により決定する．しかし，地震時・地震後すべり分布を同時に推定する粘弾性インバージョン [Tomita et al., 2018] では，平滑化度合いを調整するハイパーパラメータ・観測データの重みに対するハイパーパラメータが時間窓の数だけ増大し，計算コストが極めて高くなってしまう．また，解の非負拘束を課した最小二乗法を行う場合，非負拘束を課した場合にはABICが解析的に正しく計算できないこと [Fukuda & Johnson, 2008] に加え，計算コストが通常の最小二乗法に比べて劇的に増大する問題がある．
上記の問題の解決のため，我々は，Reversible-jump MCMC法（Rj-MCMC [Green, 1995]）による粘弾性インバージョン解析の開発に取り組んでいる．本手法では，解の非負拘束を課した場合でも事後確率分布に基づき適切に解を求めることが可能であり，未知パラメータの数も同時に探索することで，最小限の未知パラメータ数で計算可能である．富田・他 [2019，地震学会] では，Simuated Annealingを用いたRj-MCMC手法を粘弾性問題を仮定したSyntheticデータに対して適用し，通常のABICを用いた手法と同様にすべり分布を求められることを示した．しかし，解の探索の効率性と十分性に問題が残っていた．そこで，本研究では効率的な解の探索が可能なParallel Tempering (PT) [e.g., Sambridge, 2013]を用いたRj-MCMCの開発を試みた．
2. 手法・データ
本研究では，Bodin & Sambridge [2009]と同様にVoronoi cellを用いた空間分割と定式化を行った．この空間分割により平滑化のハイパーパラメータは不要となる．PTでは，複数の異なる温度状態のレプリカを作成し，解の探索の過程でその温度状態を入れ替えることで局所解を避け，効率的な解の探索を行う．計算速度を上げるため，MPIによってレプリカを並列させて計算を行った．
Syntheticデータは，弾性媒質 [Okada, 1992] を仮定した地震時すべりに対する応答として求めた．PTを実装したRj-MCMCの解の収束状態・断層すべり分布の再現性の確認するため，粘弾性媒質ではなく，１時間窓の簡易な弾性媒質を仮定している．断層形状は，沈み込み帯プレート境界断層を模した傾斜角15度・幅300 km・長さ500 kmの矩形断層面を仮定し，20 km × 20 kmの小断層を配置した．地震時すべり分布としては，海溝から75 kmの地点をピークとした最大すべり量 8 mの滑らかな滑りを与えた。観測点配置は，海溝から200 kmの距離を海岸線として，①150陸上観測点＋5海底観測点，②150陸上観測点＋150海底観測点，の2パターンを設定した．
3. 結果・議論
上記のSyntheticデータに対し，PTのレプリカ数・最大温度・サンプル数を変えてRj-MCMCインバージョンを実行した．レプリカ数・最大温度は4〜16の間，サンプル数は105〜107の間で検証を行った．結果として，サンプル数が106以上であれば，最終的なすべり分布はレプリカ数・最大温度によらず仮定したすべりを概ね再現できていることを確認した．レプリカ数が多い場合は，すべり分布の収束が早まる傾向が見られた．なお，未知パラメータ数は，観測点配置①では10個程度，②では20個程度に縮約された．これは，観測点分布による解像度の違いを反映していると考えられる．一方で，各小断層でのパラメータの確率密度関数は，観測点配置①では非常に滑らかに得られたのに対し，②では複数の急峻なピークが見られた．観測点配置②では観測データが多いため，特定の観測データへのフィッティングが良くなる場合などの局所解が現れやすくなるためだと考えられる．このような場合でも，複数のレプリカのアンサンブルによって最終的なすべり分布は仮定したすべり分布を偏りなく再現できていた．
本研究では，PTの実装により解の収束性を高め，偏りなくすべり分布を推定可能であることを確かめた．今後はレプリカ数・最大温度のパターンを更に増やすことで，観測点配置②の場合でも局所解に陥りにくくすることが可能かを検証する．また，仮定するすべり分布のパターンを増やし，通常のABICインバージョンの結果との再現性の違いを検証する．