[SVC45-26] 火山性地震のエンベロープ幅の逆問題解析を用いたタール火山の3次元散乱構造の推定
火山の内部構造の推定は、多様な火山活動の理解や火山監視において重要である。高周波数(5-10 Hz)で卓越する散乱S波は、数十m程度の短波長の構造の影響を強く受けるため、マグマや熱水の挙動を捉えやすいというメリットがある。そこで本研究では、高周波数の地震波形記録から火山内部散乱構造についての逆問題解析を用いた推定手法の開発を行い、フィリピンのタール火山に適用した。
散乱構造は平均自由行程l0と内部減衰Qiによって特徴付けられる。Kumagai et al. (JGR, 2018)では、火山性地震の5-10 Hz帯の観測波形のエンベロープ幅が、l0とQiに依存する性質を利用して1次元散乱構造モデルを推定し、深さ1 kmまでの表層の散乱と減衰が強い(l0 = 690 mとQi = 110)ことを示した。この先行研究では、各層の変数であるl0、Qi、深さについてのグリッドサーチを用いて、観測エンベロープ幅と計算エンベロープ幅のフィッティングを行った。しかしエンベロープ幅の計算に用いるモンテカルロ法(Yoshimoto, JGR, 2000)では、1つの条件での計算に時間がかかり、かつグリッドサーチでは推定散乱構造モデルの複雑化に伴い急激に変数の組み合わせの数が増大するため、より複雑な散乱構造モデルを現実的な時間で推定するのが難しいという問題があった。この問題を解決するために本研究では逆問題解析を用いて3次元散乱構造の推定を試みた。
逆問題解析では1次元モデルに設定した異常領域中のl0とQiの1次元モデルからの変化量（dl0 と dQi）を求めることで3次元モデルを推定する。まず各異常領域のl0とQiを微小量変化させて、モンテカルロ法でエンベロープ波形を合成し計算エンベロープ幅を求める。これにより異常領域中の各パラメータに対するエンベロープ幅の傾き（偏微分）を計算し、この傾きを利用して計算したエンベロープ幅と観測エンベロープ幅の残差が小さくなるように最小二乗法を用いてモデルを推定する。これにより、モンテカルロ法で計算するパラメータの組み合わせの数をグリッドサーチの場合に対して大幅に減らすことで、計算の効率化を図ることができる。本研究では2011年11月から2013年5月にタール火山の8観測点で観測された火山構造性地震の5-10 Hz 帯のエンベロープ幅を用いた。
観測エンベロープ幅と1次元モデルによる計算エンベロープ幅の残差から、タール火山についての先行研究を参考に2つの異常領域を先験的に設定した。異常領域1はKumagai et al. (GRL, 2014)によって推定された火山島の東斜面地下に位置する地震波減衰領域であり、異常領域2はYamaya et al. (BV, 2013)が推定した中央の火口湖（MCL）に位置する電気比抵抗の周囲より高い領域である。異常領域1については、Niino and Kumagai (AGU Fall Meeting, 2018)により、ほとんど同じ震源位置であるが高周波数帯で異なるエンベロープ波形を示す3イベント（2012年5月15日、6月25日、11月25日にそれぞれ発生）のVT地震のエンベロープ波形のフィッティングを行い、散乱構造の月単位での時間変化を推定している。
設定した異常領域の影響の強いVT地震の観測データを用いて、2つの異常領域のl0とQiを制約付き最小二乗法による逆問題解析で同時に求めた結果、異常領域1では1次元モデルの表層と比べて散乱が弱く(l0 = 980 m)、異常領域2では散乱が強い(l0 = 400 m)ことがわかった。Qiについては両異常領域で1次元モデルと大きく変わらない結果となった。さらに時間変化を示す3イベントを用いて異常領域1についての逆問題解析を行い、2012年5月には散乱が弱く(l0 = 1300 m)、6月には比較的強く(l0 = 660 m)、11月には弱い(l0 = 1320 m)という結果を得た。
異常領域2の強い散乱は、過去の噴火の火道の割れ目などによる強い構造不均質を反映していると考えられる。異常領域1のl0の月単位の時間変化については、Niino and Kumagai (2018)で解釈された通り、この領域に貫入したマグマの初期含水率の違いによって発泡度が変化したことが原因と考えられる。今回行った逆問題解析の結果はタール火山の様々な先行研究の解釈と整合的であり、エンベロープ幅の逆問題解析を用いた手法は3次元散乱構造の定量的な推定に有効であることを示している。
散乱構造は平均自由行程l0と内部減衰Qiによって特徴付けられる。Kumagai et al. (JGR, 2018)では、火山性地震の5-10 Hz帯の観測波形のエンベロープ幅が、l0とQiに依存する性質を利用して1次元散乱構造モデルを推定し、深さ1 kmまでの表層の散乱と減衰が強い(l0 = 690 mとQi = 110)ことを示した。この先行研究では、各層の変数であるl0、Qi、深さについてのグリッドサーチを用いて、観測エンベロープ幅と計算エンベロープ幅のフィッティングを行った。しかしエンベロープ幅の計算に用いるモンテカルロ法(Yoshimoto, JGR, 2000)では、1つの条件での計算に時間がかかり、かつグリッドサーチでは推定散乱構造モデルの複雑化に伴い急激に変数の組み合わせの数が増大するため、より複雑な散乱構造モデルを現実的な時間で推定するのが難しいという問題があった。この問題を解決するために本研究では逆問題解析を用いて3次元散乱構造の推定を試みた。
逆問題解析では1次元モデルに設定した異常領域中のl0とQiの1次元モデルからの変化量（dl0 と dQi）を求めることで3次元モデルを推定する。まず各異常領域のl0とQiを微小量変化させて、モンテカルロ法でエンベロープ波形を合成し計算エンベロープ幅を求める。これにより異常領域中の各パラメータに対するエンベロープ幅の傾き（偏微分）を計算し、この傾きを利用して計算したエンベロープ幅と観測エンベロープ幅の残差が小さくなるように最小二乗法を用いてモデルを推定する。これにより、モンテカルロ法で計算するパラメータの組み合わせの数をグリッドサーチの場合に対して大幅に減らすことで、計算の効率化を図ることができる。本研究では2011年11月から2013年5月にタール火山の8観測点で観測された火山構造性地震の5-10 Hz 帯のエンベロープ幅を用いた。
観測エンベロープ幅と1次元モデルによる計算エンベロープ幅の残差から、タール火山についての先行研究を参考に2つの異常領域を先験的に設定した。異常領域1はKumagai et al. (GRL, 2014)によって推定された火山島の東斜面地下に位置する地震波減衰領域であり、異常領域2はYamaya et al. (BV, 2013)が推定した中央の火口湖（MCL）に位置する電気比抵抗の周囲より高い領域である。異常領域1については、Niino and Kumagai (AGU Fall Meeting, 2018)により、ほとんど同じ震源位置であるが高周波数帯で異なるエンベロープ波形を示す3イベント（2012年5月15日、6月25日、11月25日にそれぞれ発生）のVT地震のエンベロープ波形のフィッティングを行い、散乱構造の月単位での時間変化を推定している。
設定した異常領域の影響の強いVT地震の観測データを用いて、2つの異常領域のl0とQiを制約付き最小二乗法による逆問題解析で同時に求めた結果、異常領域1では1次元モデルの表層と比べて散乱が弱く(l0 = 980 m)、異常領域2では散乱が強い(l0 = 400 m)ことがわかった。Qiについては両異常領域で1次元モデルと大きく変わらない結果となった。さらに時間変化を示す3イベントを用いて異常領域1についての逆問題解析を行い、2012年5月には散乱が弱く(l0 = 1300 m)、6月には比較的強く(l0 = 660 m)、11月には弱い(l0 = 1320 m)という結果を得た。
異常領域2の強い散乱は、過去の噴火の火道の割れ目などによる強い構造不均質を反映していると考えられる。異常領域1のl0の月単位の時間変化については、Niino and Kumagai (2018)で解釈された通り、この領域に貫入したマグマの初期含水率の違いによって発泡度が変化したことが原因と考えられる。今回行った逆問題解析の結果はタール火山の様々な先行研究の解釈と整合的であり、エンベロープ幅の逆問題解析を用いた手法は3次元散乱構造の定量的な推定に有効であることを示している。