はじめに

近年，GPSによる地殻変動連続観測データのインバージョン解析によって，地震の余効すべりやスロースリップイベントなどの非地震性すべりにおける断層すべり速度の時空間変化を推定することが可能になった．すべり速度時空間変化推定のための測地インバージョン手法としてSegall and Matthews [1997]によるベイズモデル（状態空間モデル）に基づく方法がある．この手法では，すべり速度が時間的に滑らかであるという事前情報を導入し，その滑らかさを支配するハイパーパラメータを赤池のベイズ型情報量規準（ABIC）の最小化によって決定している．この手法を地震の余効すべりやスロースリップイベントなどすべり速度の時間変化が大きい現象に適用した場合，すべり速度の変化が大きい期間では解が過度に平滑化され，緩やかに変化する期間では解の滑らかさが不足する傾向があり，すべり速度の時空間変化を正確に求めることが困難である．例えば，地震直後の余効すべり速度やスロースリップイベントの立ち上がりにおけるすべり速度は，過度の平滑化のために正確に推定できない可能性が高い．このような問題が起こる原因として，ハイパーパラメータが時間的に一定と仮定されているために，すべての時刻において同程度の時間平滑化が課されていることが考えられる．本研究ではこの問題を解決することを目的として，ハイパーパラメータの時間変化を考慮した新しい測地インバージョン法の開発を行った．

測地インバージョン法の開発

本研究では，非地震性すべりによる地殻変動を，状態空間モデルを用いてモデル化する．状態空間モデルはモデルパラメータと観測データを関係付ける観測モデルと，モデルパラメータの時間発展に関する先験情報を記述するシステムモデルと呼ばれる2つのモデルから成る統計モデルである．

半無限弾性体中の断層面上の点ξにおける時刻tでの断層すべりs（ξ,t）をすべりの空間基底関数Bk（ξ）と時間変化する係数ck（t）を用いて次のように表す．

Okada [1985]のディスロケーション理論を用いて断層すべりs（ξ,t）を観測された地殻変動データと結び付けることにより，次のような観測モデルが定式化できる:

ここで， xnはck（tn），ck（tn）を含んだモデルパラメータベクトル， ynはデータベクトルであり，添字nは時刻tnに対応している．

すべりの時間平滑化に関する事前情報として，本研究では，「すべり速度がほぼ一定」または「すべり加速度がほぼ一定」というモデルを採用する．前者の場合， k番目の空間基底関数に対する事前情報は次のようになる．

ここでαnは時間変化を考慮したハイパーパラメータである．（3）を用いて，システムモデルは次のように定式化できる．

ハイパーパラメータの時間変化を推定するために，本研究ではハイパーパラメータを時間変化する確率変数Anと考える．Anは〓のいずれかの値を取るものとする．Anの時間変化に対しても事前分布を導入することにより，次のような階層ベイズ型の時系列モデルを構築した．

（5）,（6）はそれぞれシステムモデル，観測モデル，（7）はハイパーパラメータの時間変化に対する事前分布である．

（5）,（6）,（7）を用いてAnとxnの時間変化を求めるための計算アルゴリズムとしては，Higuchi and Fukuda （2003）, Fukuda et al. （2004）によって開発されたモンテカルロ混合カルマンフィルタを用いた．このアルゴリズムでは，ハイパーパラメータの予測分布p（A1:n|y1:n-1,β）とフィルタ分布p（A1:n|y1:n-1,β）を多数の実現値（粒子）で近似表現する．予測分布とフィルタ分布を近似する粒子は，「予測」と「フィルタ」と呼ばれるアルゴリズムを繰り返すことにより，逐次的に推定することができ，最終的にAnの事後分布p（A1:Ne|y1:Ne,β）を近似する粒子が得られる．この分布を近似する粒子が得られると，（5）,（6）からカルマンフィルタのアルゴリズムを用いてp（x1:Ne|y1:Ne,A1:Ne,β）を求めることができる．最後にxnの事後分布p（x1:Ne |y1:Ne,β）はp（x1:Ne|y1:Ne,A1:Ne,β）の全粒子に対する平均を取ることによって求められる．

本研究で開発したインバージョン手法の有効性を確かめるために，余効すべりとスロースリップイベントを想定して作成したsynthetic dataを用いた数値実験を行った．従来の手法を用いた場合，地震直後の余効すべり速度やスロースリップイベントの立ち上がりにおけるすべり速度は，過度の平滑化のために真の値を再現できない．一方，本研究で開発した手法を用いた場合，より真の値に近いすべり速度の推定値が得られた．

2003年十勝沖地震の余効すべりへの適用

本研究で開発した手法を2003年9月25日（UTC）に発生した十勝沖地震の余効すべりに伴う地殻変動に適用した．推定された地震後3日間のすべり速度の時空間変化を図1（a）に示す．比較のために従来の手法による解析も行った（図1（b））．2種類の手法から得られた結果を比較すると，本研究で開発した手法から得られたすべり速度（図1（a））は従来の手法から得られたもの（図1（b））に比べて顕著に大きいことが分かる．これは従来の手法により推定されたすべり速度が過度な時間平滑化のために過小評価されている一方，本手法で推定されたすべり速度にはその影響が少ないためであると考えられる．数値実験の結果を考慮すると，本手法で推定されたすべり速度の方が真の値を良く再現していると考えられる．

このインバージョンで得られたすべり・すべり速度の推定結果を基に余効すべり領域における剪断応力変化τとすべり速度Vの関係∂τ/∂（log V）を推定した結果，∂τ/∂（log V）〜0.07MPaという結果が得られた．速度・状態依存摩擦法則において定常状態を仮定するとこの値は〓n（a.b）に相当する．本研究で得られた推定値は従来の測地インバージョン法に基づきMiyazaki et al. [2004]によって得られた∂τ/∂（log V）〜0.2MPaに比べて顕著に小さい．本研究で求められた値との違いは，Miyazaki et al. [2004]によって求められた値が，すべり速度の過小評価により，過大評価されていることに起因すると考えられる．

2002年房総沖スロースリップイベントへの適用

上記の手法を2002年10月に房総沖で発生したスロースリップイベントに伴う非定常地殻変動に適用した．比較のために，ハイパーパラメータを一定と仮定した従来の手法を用いた解析も行った．結果を図2に示す．従来の手法を用いて推定されたすべり速度は，本研究で開発した手法から得られたすべり速度に比べて，イベント時により平滑化されている一方，イベントの前後では明らかに滑らかさが不足している．一方，本研究で開発した手法を用いて推定されたすべり速度にはこのような問題は見られない．数値実験の結果を考慮すると，本手法で推定されたすべり速度の方が真の値を良く再現していると考えられる．この結果は，本研究で開発した手法を適用することにより，スロースリップイベントの時間発展に関するより詳細な情報が得られる可能性があることを示す．

図1: 2003年十勝沖地震後3日間のすべり速度の時空間変化．（a）:新しい手法から得られたすべり速度．（b）:従来の手法から得られたすべり速度．

図2: （a）,（b）,（c）: （d）に示された断層面上の3点におけるすべり速度の時間変化．赤は新しい手法，青は従来の手法から得られたすべり速度を示す．（d）:すべり速度を計算した点の位置．カラースケールは2002/8/16から11/23までの累積すべり量

本論文は5章からなる．第1章はイントロダクションであり，地殻変動連続データから非地震性すべりにおける断層すべりの時空間変化を推定する手法について，従来のインヴァージョン解析の問題点が記述されている．これまでの手法では，すべり速度が時間的に滑らかであるという事前情報を導入し，その滑らかさを支配するハイパーパラメータを赤池のベイズ型情報量基準（ABIC）の最小化によって決定しているが，この手法を地震の余効すべりやスロースリップイベントなどすべり速度の時間変化が大きい現象に適用した場合，変化が大きい期間では解が過度に平滑化されて時間分解能が低下したり，逆に緩やかに変化する期間では解が不安定になることがあり，すべり速度の時空間変化を正確に求めることが困難であった．第2章では，この問題を解決することを目的として新たに開発を行ったインヴァージョン法について述べられている．新しい手法では，ハイパーパラメータの時間変化を考慮し，階層ベイズ型の時系列モデルを構築することにより，平滑化の強さをデータの変化の度合いに応じて適切に決定できるようにした．第3章では，その方法を2003年9月25日に発生した十勝沖地震の余効すべりに伴う地殻変動に適用し，本震直後のすべり速度が，従来の手法から得られたものに比べて大きいことを示した．このような結果が得られたのは本手法で推定されたすべり速度が過度な時間平滑化を受けにくいためであると考えられる．この結果を基に，余効すべり領域における速度・状態依存摩擦則のパラメータの推定を行った．第4章では，2002年10月に房総沖で発生したスロースリップイベントに適用し，従来の手法に比べて明らかに時間分解能が高くかつ安定した解が得られることを示した．第5章では，結論と今後の課題について述べられている．

本研究は，プレート境界のすべりの物理解明に重要な貢献をするものであり，また，最近さかんに行われている地震発生の数値シミュレーションの信頼性を高めていくものでもある．なぜなら，地震発生の数値シミュレーションでは，プレート境界での速度・状態依存摩擦則のパラメータをそれほど強い根拠がない値に仮定しているが，本研究を進めていけば，パラメータの値に対して強い拘束を与えることができるようになるからである．

なお，本論文第2章の一部は，樋口知之氏，宮崎真一氏，加藤照之氏との共同研究であるが，論文提出者が主体となって開発及び検証を行ったもので，論文提出者の寄与が十分であると判断する．

以上のことから，本研究は博士の学位を受けるのにふさわしい優れた研究と認める．