統計数理（2015）
第 63 巻 第 1 号 105–127
c 2015 統計数理研究所
特集「地震予測と統計モデル」［研究詳解］GPS データの逆解析と地震の発生予測
松浦 充宏1 ・野田 朱美2（受付 2014 年 12 月 9 日；改訂 2015 年 2 月 6 日；採択 2 月 9 日）

要旨

地震の原因は岩石の脆性破壊であるから，その発生を決定論的に予測することはできない．そこで，1890 年代以降，過去の地震活動データから未来の地震の発生を確率論的に予測する様々な統計モデルが提案されてきた．一方，1990 年代に入ると地震の発生過程を物理的に理解しようとする研究分野が誕生し，現在では様々な断層構成則（摩擦則）を支配法則とする地震発生の物理モデルが提案されている．これらの物理モデルでは，地震の発生は断層面に働く剪断応力の解放過程として記述される．従って，物理モデルと統計モデルは応力を介して密接に結びついているはずである．しかし，地殻の応力状態を知るのは容易なことではない．現時点では，地震のメカニズム解や CMT（Centroid Moment Tensor）解の逆解析から地殻の応力パターンが，データの逆解析から応力変化が推定できるに過ぎない．
また，GPS（Global Positioning System）
本稿では，主に後者について，その考え方と手法を解説し，実際の解析結果を地震の発生予測という観点から考察する．
キーワード：地震発生，統計モデル，物理モデル，地殻応力，GPS データ，逆解析．
1.

はじめに

地震の発生は地殻に蓄えられた弾性歪みエネルギーの一部が断層の剪断破壊によって一気に解放される過程である（Reid, 1910; Knopoﬀ, 1958; Savage, 1969）．地震破壊は断層強度の急激な低下であり，強度を超えて応力は存在し得ないので，それは応力の低下でもある．破壊面での応力低下によって生じた変形は弾性波速度で周辺域へと伝わり，そこに蓄えられていた歪みエネルギーを解放する．重要なのは，その解放されたエネルギーが破壊の進展を加速し，周辺域での更なるエネルギー解放を引き起こす（Kostrov, 1966; Freund, 1972）という点である．従って，何かのきっかけで開始した破壊が加速し，強い地震波を放射する高速破壊（地震）へと発展するには，震源域に解放すべき充分な弾性歪みエネルギーが蓄えられていなければならない．

或る弾性領域 V の歪みエネルギーは，応力 σij と歪み εij を用いて E ≡ 12 V σij εij dV と定義されるが，弾性体の応力と歪みの線形関係（σij = cijkl εkl 或いは εij = sijkl σkl ）により，

E = 12 V cijkl εij εkl dV 或いは E = 12 V sijkl σij σkl dV と表すことができる．従って，弾性歪み
エネルギーの解放過程である地震破壊は，震源域の応力（或いは弾性歪み）を減少させるように発生する．逆に言えば，地震の CMT（Centroid Moment Tensor）解は，その地震の規模に応じた空間領域の応力場に関する情報を持っていることになる．このような基本的考えに基づいて，Terakawa and Matsu’ura（2008）は，地震の CMT データから地殻の応力パターン（応力テンソ
1
2

統計数理研究所：〒 190–8562 東京都立川市緑町 10–3
構造計画研究所 防災ソリューション部：〒 164–0011 東京都中野区中央 4–5–3

106

統計数理

第 63 巻

第 1 号 2015

ルの配向）を推定する逆解析手法を定式化した．また，Terakawa and Matsu’ura（2010）は，この（1997.1.1〜2007.1.31）
逆解析手法を防災科学技術研究所の F-net モーメントテンソル・カタログの約 15,000 個の中小地震（M 3.5–5.0）に適用し，日本列島域の応力パターンを求めている．個々の中小地震は，地殻内に分布する強度の弱い面を選んで発生するので，必ずしも真の応力場を反．それにも拘わらず，こうして求め映しているわけではない（例えば，Twiss and Unruh, 1998）
た応力パターンは，大地震を発生させ，長期的な地殻変形を引き起こす広域応力場（テクトニック応力場）と調和的であることが確かめられている．
大きな地震が発生すると，震源域とその周辺部では応力と強度の変化が生じ，一連の中小地震が誘発される．こうした地震活動（余震活動）は普遍的特性を持ち，経験的に導かれた大森 - 宇に従って時間と共に減衰していくことが知られている．余津の公式（Omori, 1894; Utsu, 1961）
震活動の一因である本震に伴う静的な応力場の変化は，断層面の幾何学的形状とすべり分布が与えられれば，食い違いの弾性論（elastic dislocation theory）に基づいて理論的に計算できる．
断層面が滑らかなら，そこに働いていた剪断応力は全体的に低下する（高速破壊領域内のすべ
り残しによる局所的応力上昇は物理的に考え難い）．一方，破壊領域のすぐ外側では，剪断応
力の局所的上昇（応力集中）
が起こる．その応力集中を解消するように，破壊領域を囲む浅部の脆性領域では余震が，深部の脆性 - 塑性遷移領域では余効すべりが発生する．この種の余震活
動は，破壊核形成に伴う前震活動と本震の発生に引き続く一連の過程として，合理的に説明できる（Matsu’ura et al., 1992; Shibazaki and Matsu’ura, 1992）．断層面が不規則な形状をしていると，地震時のすべりは応力の局所的不均一を作り出すので，主破壊面から分岐した副次的破
壊を引き起こし，断層面に沿った損傷領域（damage zone）が形成されると考えられる（例えば，Dieterich and Smith, 2009）．また，断層面が滑らかな場合であっても，すべり速度弱化の摩擦則を仮定すれば，地震時すべりの分布が不規則になるので，結果として断層面の応力不均一が．Nakanishi（1992）は，その応力不作り出される（Carlson and Langer, 1989; Carson et al., 1991）
均一が地殻下の粘性流体の応力緩和によるローディングで逐次解消される過程を余震活動と考えたが，中規模以下の地震では粘性応力緩和による効果は無視し得るほど小さいので，これは普遍的な説明とは言えない．
余震は，本震の断層面に限らず，そこから外れた剪断応力の低下領域でも発生する．この種の余震活動は，間隙流体が本震で損傷を受けた領域に拡散していくことで強度低下を引き起こした結果と考えるのが合理的であろう（Nur and Booker, 1972; Bosl and Nur, 2002）．実際，高圧注水によって誘発された地震活動の解析結果（Terakawa, 2014）は，こうした考えを強く支持している．余震を地殻内の既存の弱面での剪断破壊とすると，その起こり易さはクーロンの破壊関で測ることができる．クーロンの破壊関数は対象とする断層面に数（Coulomb Failure Function）
働く剪断応力と断層面の摩擦強度の差として定義され，その本震の発生に伴う変化（ΔCFF ）は，間隙流体が存在する場合，ΔCFF = Δτr − μ(Δσn − Δp) と記述される．ここで，μ は岩石の摩擦係数，Δτr と Δσn は断層面に働く剪断応力と法線（圧縮）応力の変化，そして Δp は間隙流体である圧の変化である．この内，摩擦係数は岩石の種類や環境条件によらずほぼ一定（μ ∼= 0.6）
ことが実験的に確かめられている（Byerlee, 1978）．また，剪断応力と法線応力の変化は，本震のすべり分布が与えられれば，理論的に見積もることができる．しかし，本震後の新たな応力場の下で損傷を受けた地殻内を間隙流体が拡散していく過程は未だ定量的には解明されていな
いので，間隙流体圧 Δp が時間的・空間的にどのように変化するのか不明である．ちなみに，現在広く用いられている Reasenberg and Simpson（1992）や Stein et al.（1992）の ΔCFF の評価は，本震による損傷や間隙流体の拡散は考えず，あまり現実的とは
式（ΔCFF = Δτr − μ Δσn ）．
言えない幾つもの仮定を設けて得られた便宜的なものに過ぎない（例えば，Beeler et al., 2000）
余震発生の物理的メカニズムは未だ完全には解明されていないが，結局，余震活動は本震GPS データの逆解析と地震の発生予測

107が周辺域に作り出した応力と強度の不均一を平準化して新たな平衡状態に移行する必然的な過程と理解すべきであろう．そうすると，余震もまた（程度の差こそあれ）周辺域に応力と強度の不均一を作り出すのだから，大森 - 宇津の余震発生率に関する公式（n(t) = K(t + c)−p ）は，本震に限らず，余震に対しても適用されなければならない．実際，Ogata（1988）が提案（マグニ
した ETAS（Epidemic Type Aftershock Sequence）モデルでは，全ての地震がその規模
チュード）に応じた余震活動を伴うとしている．具体的には，時刻 t = ti に発生した地震のマ
グニチュードを Mi として，時刻 t (> ti ) に於けるマグニチュード Mz 以上の地震の発生率は
λ(t) = μ + ti <t K exp[α(Mi − Mz )](t − ti + c)−p と記述される（但し，μ は常時地震活動度）．
つまり，或る時刻までに発生した余震のデータから 5 つのパラメータ（μ, K, α, c, p）の最適値を推定してやれば，それ以降の余震活動の推移は ETAS モデルで予測可能ということである．このモデルを，余震活動だけでなく，時間的にも空間的にも複雑に変化する一般の地震活動に適用するには，全てのパラメータを地殻の構造特性と応力状態を表す場の量（つまり場所の関数）として扱う時空間 ETAS モデル（Ogata, 1998, 2004）へと発展させる必要があった．現時点では，場の量として推定したパラメータの値が地殻の構造特性や応力状態をどのように反映しているか不明であるが，その時間変化は直接的或いは間接的に地殻応力状態の変化を反映していると考えるのが自然であろう．
1990 年代以降，GPS データの解析を通じて，地殻応力状態の時々刻々の変化を知ることができるようになった．そこで，地震活動データの解析結果と GPS データの解析結果を比較す
れば，地震発生確率予測に不可欠な地殻の構造特性と応力状態に関する新たな知見が得られる
ものと期待される．本稿の以下の章では，最近の一連の研究（Matsu’ura et al., 2007; Noda andに基づき，GPS データから逆解析により地殻応力場の時空
Matsu’ura, 2010; Noda et al., 2013）
間変動を推定する方法について詳しく解説する．
2.

GPS データに含まれる情報

通常，地殻は線形弾性体と見なせるので，その歪み変化 Δεkl が分かれば弾性スティフネス・
テンソル cijkl を介して応力変化 Δσij (= cijkl Δεkl ) を知ることができる．日本では，1880 年
代以降，全国規模の測地測量（測点の緯度・経度を決める三角測量と標高を決める水準測量）
が繰り返し行われるようになった．しかし，三角測量から得られる地殻歪みの情報は水平 3 成分の内の剪断歪みについてだけである（Tsuboi, 1932）．1960 年代に光波測距の技術が確立され，
1970 年代には三角測量に替わって三辺測量が行われるようになった．三辺測量では，それまでの三角測量とは異なり，水平歪み 3 成分全てについての情報が得られる．やがて，1990 年代に入ると GPS による測地測量が行われるようになり，測点の変位が直接測定されるようになった．変位場 ui (x) が与えられれば，歪みは変形勾配テンソル Dij = ∂ui /∂xj の対称部分として
εij = 12 (∂ui /∂xj + ∂uj /∂xi ) と計算される．このことから，GPS 変位データに含まれる情報は三辺測量データに含まれる情報と本質的には変わらないように思える．それでは，変形勾配テ
ンソルの反対称部分として計算される回転 ωij = 12 (∂uj /∂xi − ∂ui /∂xj ) はどんな情報を持っているのだろうか．
いま，鉛直上向きの x3 軸を持つ直交直線座標系 (x1 , x2 , x3 ) をとると，或る地殻ブロック内（1965, Chapter 4）に従い，次のような歪みの経路積分での任意の点 Pk に於ける変位は，Fung表すことができる：
 Pk
0
(2.1)
uki = u0i − (xkj − x0j )ωij+
Uiα (x)dxα (i, j = 1, 2).
P00
ここで，u0i と ωij
はブロック内の起点 P0 に於ける変位と回転を表し，Uiα (x) は歪みを用いて

統計数理

108

第 63 巻

第 1 号 2015

以下のように定義される量である：
Uiα (x) = εiα (x) + (xkj − xj )[∂εiα (x)/∂xj − ∂εjα (x)/∂xi ].

(2.2)

P
0
起点 P0 を u0i = 0 且つ ωij= 0 であるような場所に設定すれば，（2.1）式は uki = P0k Uiα (x)dxαと書けるので，歪みから変位場を一意的に復元できるかのように思える．しかし，実際には，
0はブロックの剛体的な並進成分 ci 及び回転成分 γij を含んでいる．その（2.1）式の u0i 及び ωij
ことを明示的に記せば，
 Pk
0
0
(2.3)
uki = (u i + ci ) − (xkj − x0j )(ω  ij + γij ) +Uiα (x)dxα (i, j = 1, 2).
P0





ここで，u i 及び ω ij はそれぞれ変位 ui 及び回転 ωij の空間変動部分を表す：ui (x) = u i (x) + ci ,00（u i = 0 且つ ω  ij = 0）に設定
ωij (x) = ω  ij (x) + γij ．従って，起点 P0 を変形が無視できる地点
式はしたとしても，
（2.3）
 Pk
(2.4)
uki = ci − (xkj − x0j )γij +Uiα (x)dxα (i, j = 1, 2)
P0となるだけで，剛体並進・回転運動の影響は依然として残ったままである．これは，変位及び回転の空間変動部分は歪みを用いて表せる（∂ω  ij /∂xα = ∂εjα /∂xi − ∂εiα /∂xj ）のに対して，それらの定数部分は歪みと独立である（∂γij /∂xα = 0）ことに起因している．つまり，GPS 観測から得られる変位データには地殻変形の情報に加えてブロックの剛体並進・回転運動の情報が含．GPSまれているが，それらは互いに独立なので分離可能ということである（Noda et al., 2013）
変位データからブロックの剛体並進・回転運動を除去し，地殻変形に関する情報のみを取り出
すには，例えば，隣接する 3 測点の水平変位データをその 3 測点が作る三角形領域の平均水平歪みに変換してやれば良い（具体的変換公式は Noda et al., 2013 の Appendix を参照）．この変換は線形なので，GPS 変位データに含まれる情報が歪められることはない．
3.

地殻変形の物理モデル

日本列島で現在進行中の地殻変形は，今や GPS 観測を通じて直接捉えることができる．しかし，GPS 観測は地表面に沿って 2 次元的に行われるため，変位ベクトルの鉛直勾配 3 成分を得ることができず，従来の幾何学的歪み解析では地表での水平歪み 3 成分しか求めることができない．しかも，こうして求めた幾何学的歪みは応力と関係する弾性歪みと関係しない非弾性歪みを足した全歪みなので，地殻応力を知るには全歪みから弾性歪みだけ取り出す必要がある．
Noda and Matsu’ura（2010）
は，地殻変形の物理モデルを介して GPS データを解析することにより，これらの問題を一挙に解決できることを示した．以下では，その基本的考え方について解説する．
3.1 非弾性歪みの発生と弾性歪みの解放
先ず，地殻を内部に様々な弱面を含む線形弾性体としてモデル化し，その弱面での脆性破壊や塑性変形が地表の変形を引き起こすと考える．脆性破壊や塑性変形などの非弾性的な物理過程に伴って発生する歪み Δεaij (ξ) は応力と関係しない非弾性歪みであり，それに周辺媒質の弾性スティフネス・テンソル cijkl を作用させると，対応するモーメントテンソル密度 Δmij (ξ) が：定義される（Backus and Mulcahy, 1976）(3.1)

Δmij (ξ) = cijkl Δεakl (ξ)(i, j = 1, 2, 3).

GPS データの逆解析と地震の発生予測

109

ちなみに，モーメントテンソル密度の空間導関数 ei = −∂mij /∂xj は，線形弾性を仮定した連続体の運動（平衡）方程式に於いて，周辺領域を変形させる仮想的な力（等価物体力）と見なされる．
地殻内の領域 V に無数の弱面が分布している場合，モーメントテンソル密度 Δmij (ξ) による
地表変位場は，単位のステップ的な力に対する変位応答関数を Gip (x; ξ)，その震源座標 ξq に
（具体的表現式は Yabuki and Matsu’ura, 1992 の Appendix A を参関する偏導関数を Gip,q (x; ξ)照）
として，以下のように計算される：(3.2)
ui (x) =Gip,q (x; ξ)Δmpq (ξ)dξ (i = 1, 2, 3).
V

従って，GPS 変位データから地殻内のモーメントテンソル密度の 3 次元分布 Δmpq (ξ) を推定す
式の逆問題を解いてモーメントテンソル密度分布が得られれば，る逆問題が設定できる．（3.2）
それに弾性コンプライアンス・テンソル sijpq を作用させることにより，地殻の 3 次元非弾性歪
み変化は
(3.3)

Δεaij (x) = sijpq Δmpq (x)(i, j = 1, 2, 3).

と求まる．一方，モーメントテンソル密度分布が与えられれば，地殻の 3 次元弾性歪み変化
Δεeij (x) は，Hijpq (x; ξ) を単位のモーメントテンソルに対する歪み応答関数（具体的表現式はHashima et al., 2014 を参照）
として，

(3.4)
Δεeij (x) =Hijpq (x; ξ)Δmpq (ξ)dξ (i, j = 1, 2, 3)
Vと計算することができる．ちなみに，幾何学的解析から得られる歪みは，上記の弾性歪みと非弾性歪みを足したトータルな歪み Δεij (= Δεeij + Δεaij ) である．
良く発達した活断層のように弱面の位置と形状が分かっている場合，その弱面 Σ(η) での変
位の食い違い Δw(η) に対応する非弾性歪み変化は，デルタ関数 δ(ξ − η) を用いて
(3.5)

Δεakl (ξ) = Δw(η)νk (η)nl (η)δ(ξ − η)(i, j = 1, 2, 3)と表せる．ここで，n(η) は弱面の法線ベクトル，ν(η) は変位の食い違い方向（弱面のすべり方
式から，対応するモーメントテンソル密度は向）ベクトルを表す．従って，
（3.1）
(3.6)

Δmij (ξ) = cijkl Δw(η)νk (η)nl (η)δ(ξ − η)(i, j = 1, 2, 3),

また，（3.2）式から，地表変位は
(3.7)
ui (x) = cpqkl
Gip,q (x; η)Δw(η)νk (η)nl (η)dη(i = 1, 2, 3)Σとなる．つまり，この場合は，GPS 変位データから既知の弱面 Σ(η) のすべり分布 Δw(η) を推式の逆問題を解いてすべり分布が得られれば，（3.5）
定する馴染みの逆問題が設定できる．（3.7）
式から弱面に局在化した非弾性歪み変化が求まり，また，（3.4）式と（3.6）式から地殻の 3 次元弾
性歪み変化が(3.8)
Δεeij (x) = cpqkl
Hijpq (x; η)Δw(η)νk (η)nl (η)dηΣと計算される．
こうして GPS データから弾性歪み変化 Δεeij (x) が分かれば，それに弾性スティフネス・テン
ソル cijkl を作用させることで，地殻応力変化の 3 次元分布 Δσij (x) = cijkl Δεekl (x) を知ることができる．ここで重要なのは，地殻内の脆性破壊や塑性変形は，基本的に，周辺域の弾性歪み

統計数理

110

第 63 巻

第 1 号 2015

エネルギーを減少させるように発生するということである．つまり，地殻内の脆性破壊や塑性
変形は，局所的な非弾性歪みを生みだすと同時に，周辺域に蓄積されていた弾性歪み（応力）を解放する過程と理解することができる．
3.2 プレート運動が作り出す地殻弾性歪み
地殻内の脆性破壊や塑性変形が地殻応力の解放過程なら，それとは逆の地殻応力を生成する過程が存在するはずである．1 万年オーダーの時間スケールで見れば，プレート境界はその強度に見合う抗力（摩擦抵抗）を受けながら定常すべり運動を続けている．Matsu’ura et al.（2014）は，理論的考察に基づき，この摩擦抵抗を受けながらの定常すべり運動がプレート境界周辺域にテクトニック環境に応じた固有の地殻応力場（テクトニック応力場）を形成すると結論付けている．
一方，百年オーダーの時間スケールで見れば，プレート境界の摩擦強度が相対的に大きい部分（アスペリティ）
では，定常すべりからの遅れによる応力蓄積が徐々に進行し，やがてそれが限界に達すると急激な応力解放を伴う断層運動（地震）が生じてすべり遅れを解消する（Hashimoto．こうした地震発生サイクルに伴う地殻変動は，
and Matsu’ura, 2002; Hashimoto et al., 2014）1990 年代以降，GPS 観測を通じて直接捉えることができるようになった．
プレート境界の地震であっても，地震時の地殻変形に関しては，前節の地殻内活断層での変位の食い違いと同じようにモデル化できるので，地表変位の表現式（3.7）を用いて GPS データ式を用いてプレート境界に局在化した非弾性歪み変化を逆化し，得られたすべり分布から（3.5）
を求め，（3.8）式を用いて地殻の 3 次元弾性歪み変化を計算すれば良い．問題は地震間のすべり遅れに伴う地殻変形で，この場合は，プレート境界の非地震性すべりやアセノスフェアの応力緩和を考慮する必要がある．以下では，Noda et al.（2013）に従って，地震間のすべり遅れに伴う地殻変形のモデル化について概説する．
先ず，弾性的地表層（リソスフェアに相当）と粘弾性的基盤層（アセノスフェアに相当）から成（プレート境界面に相当）を設定する
る半無限領域を考え，弾性的地表層を切断する弱面 Σ(η)．いま，プレート境界面 Σ(η) の法線ベクトルを n(η)，定常すべりベクトルの方向を ν(η)（図 1）とすると，時刻 t = 0 に開始したプレート境界でのすべり運動 w(η, t) による累積地表変位は，履歴積分の形で次のように表現される：
 t

(3.9)
ui (x, t) =dtUi (x, t − t ; η, 0)ẇ(η, t )dη (i = 1, 2, 3).
0Σ

但し，ẇ はすべり w の時間微分（すべり速度）を表す．また，Ui は単位すべりに対する粘弾性変位応答関数，(3.10)Ui (x, t; η, t ) = Gip,q (x, t; η, t )cpqkl νk (η)nl (η),

を表す（具体的表現式は Fukahata and Matsu’ura, 2006 を参照）．
次に，プレート境界 Σ を地震発生領域 Σs とそれ以外の定常すべり領域 Σ − Σs に分け，地震発生領域では w(η, t) = Vpl (η)t + ws (η, t)，それ以外の領域では w(η, t) = Vpl (η)t とする．但
し，Vpl (η) はプレート相対運動速度，ws (η, t) はそのプレート相対運動からの摂動を表す．また，簡単のため，すべりの摂動 ws (η, t) は一定のレート ẇs (η) で一定期間 T だけすべり遅れ続
けた後，地震時のすべりで一気に元に戻るものと仮定する：(3.11)

ws (η, t) = −ẇs (η)t + ẇs (η)Tn
H(t − kT )0 ≤ t < (n + 1)T.

k=1

ここで，H(t) はヘヴィサイドの階段関数を表す．すると，時刻 t = nT（ τe = アセノスフェ

GPS データの逆解析と地震の発生予測

111

図 1．沈み込み帯のプレート境界モデル．黒矢印 n はプレート境界の法線ベクトル，ν はすべ
り方向ベクトルを表す．影を付けた部分 Σs はプレート境界上の地震発生領域．
Fig. 1. Geometry and notation of a plate interface model. The shaded portion of theplate interface Σ represents a potential seismogenic region Σs . The black arrows,n and ν, denote the unit normal and slip-direction vectors, respectively. The
gray arrow shows the horizontal motion of the underlying plate relative to theoverlying plate.

アの実効的緩和時間）に発生した地震の直後から次の地震まで（nT < t < (n + 1)T ）の地表変位
は，（3.9）式から，

(3.12) ui (x, t) = ci (x) + (t − τe )Ui∞ (x; η)Vpl (η)dη − (t − τe )Ui∞ (x; η)ẇs (η)dηΣΣs[Ui (x, t − nT ; η, 0) − Ui (x, 0; η, 0)]ẇs (η)dη+TΣsとなる．但し，ci (x) は初期（0 ≤ t ≤ τe ）の累積変位，Ui∞ (x; η) は t → ∞ での粘弾性変位応答（イメージとしては水に浮かぶ弾性板の応答）を表す．上式を時間 t で微分すると，プレート境
界でのすべり遅れによる地震間（nT < t < (n + 1)T ）の地表変位速度の表現式が次のように得られる：(3.13)u̇i (x, t) =Σ

Ui∞ (x; η)Vpl (η)dη −ΣsUi∞ (x; η)ẇs (η)dη + T

ΣsU̇i (x, t − nT ; η, 0)ẇs (η)dη.

ここで，右辺第 1，第 2，及び第 3 項は，それぞれ，プレート境界全域での定常すべり，地震発
生領域での定常すべり遅れ，及び地震後のアセノスフェアの応力緩和による地表変位速度を表（3.13）
している．第 3 項は時間と共に指数関数的に減衰していくので，地震後しばらく経てば，
式の代わりに以下の近似式を解析に用いることができる：

(3.14)u̇i (x) =Ui∞ (x; η)Vpl (η)dη −Ui∞ (x; η)ẇs (η)dη.
ΣΣs

この場合，右辺第 1 項は理論的に計算可能なので，GPS 変位速度データからプレート境界でのすべり遅れのレート ẇs (η) を推定する逆問題が設定できる．（3.14）式の逆問題を解いてすべり式からプレート境界に局在化した非弾性歪み速度分布
遅れレート分布が得られれば，（3.5）(3.15)

ε̇akl (ξ) = ẇs (η)νk (η)nl (η)δ(ξ − η)が求まり，また，
（3.4）
式と
（3.6）
式から地殻の 3 次元弾性歪み速度分布が(3.16)ε̇eij (x) = cpqkl
Hijpq (x; η)ẇs (η)νk (η)nl (η)dη
Σ

統計数理

112

第 63 巻

第 1 号 2015と計算される．
4.

逆問題の定式化と解法

先ず，前章の
（3.2）
式，
（3.7）
式，或いは（3.14）式をパラメータ化するために，モーメントテン
ソル密度分布 Δmpq (ξ)，弱面のすべり分布 Δw(η)，或いはプレート境界のすべり遅れレート分の基底関数の重ね合わせとして表現する．具体的には，モーメン
布 ẇs (η) を適切な有限個（R）
トテンソル密度分布 Δmpq (ξ) の場合，Φr (ξ) を 3 次元の基底関数として，(4.1)

Δmpq (ξ) =

R


apq
r Φr (ξ),

r=1

弱面のすべり分布 Δw(η) 或いはプレート境界のすべり遅れレート分布 ẇs (η) の場合，Ψr (η) を
2 次元の基底関数として，(4.2)

Δw(η) =Rar Ψr (η)

或いは ẇs (η) =

r=1

R


ar Ψr (η)

r=1と表す．これらの表現式を（3.2）
式，（3.7）式，及び（3.14）式に代入すれば，地表の点 xj での変の線形結合とし
位或いは変位速度は，それぞれ，基底関数の重ね合わせの係数（apqr 或いは ar ）て以下のように記述される：

3 
3 
R

pq
pq
(4.3)
ui (xj ) =Uir(xj )apq
,
但し
U
(x
)
=
Gip,q (xj ; ξ)Φr (ξ)dξ,
j
r
ir
V

p=1 q=1 r=1(4.4)

ui (xj ) =RUir (xj )ar ,但し Uir (xj ) = cpqkl

r=1

(4.5)

u̇i (xj ) − u̇ssi (xj ) = −RΣGip,q (xj ; η)Ψr (η)νk (η)nl (η)dη,

但し Uir (xj ) =Uir (xj )ar ,r=1

ここで，プレート境界の定常すべりによる地表変位速度 u̇ssi (xj ) ≡論的に計算される量なので左辺に移項した．

Σs

Ui∞ (xj ; η)Ψr (η)dη.

Σ

Ui∞ (xj ; η)Vpl (η)dη は理

4.1 観測方程式とモデル誤差
obs
GPS で観測された地表の点 xj での変位或いは変位速度を uobs
i (xj ) 或いは u̇i (xj ) とすると，（4.5）
式の理論モデルを用い，以下の観測方程式を得る：（4.3）式，（4.4）式，及び
(4.6)

uobs
i (xj ) =

3 
3 
R


pq
sys
Uir
(xj )apq+ eresr + eii ,p=1 q=1 r=1

(4.7)

uobs
i (xj ) =R

Uir (xj )ar + esys
+ eres
i ,
i

r=1

(4.8)

ss
u̇obs
i (xj ) − u̇i (xj ) = −RUir (xj )ar + esys
+ eresi .
ir=1

但し，esysは理論モデルには含まれない地殻ブロックの並進・回転運動に起因する系統誤差，eresiiGPS データの逆解析と地震の発生予測

113

はそれ以外のランダム誤差を表す．これらの観測方程式は以下のようなベクトル形式で記述することができる（具体的な対応関係については省略）：(4.9)

uobs = U a + esys + eres .

我々の問題は，データ uobs と係数行列 U が与えられたとき，ランダム誤差 eres の統計的性質を仮定して，モデルパラメータ a の最適値を推定することである．
そのためには，先ず，モデルでは説明できない系統誤差 esys を除去しなければならない．具体的には，ドロネー三角形分割（Delaunay triangulation）の方法で GPS 観測点の三角網を作り，各三角形を構成する観測点の変位或いは変位速度ベクトル uobs に Noda et al.（2013）の Appendixに示される変位から歪みへの変換行列 R を掛けてやれば良い：(4.10)

Ruobs = RU a + Resys + Reres .

そうすると，Resys = 0 なので系統誤差 esys は自動的に除去され，変換後の観測方程式は，
d = Ruobs , H = RU , e = Reres として，(4.11)

d = Ha + eと記述できる．但し，d, H , a, e の次元はそれぞれ n × 1, n × m, m × 1, n × 1 とする．こうして地殻ブロックの並進・回転運動に起因する系統誤差は除去されたので，後は観測方程式（4.11）の誤差 e の統計的性質をどのように設定するかである．変換される前のランダム誤差 eres は GPS
測量に於ける観測誤差 eobs と理論モデルが不完全なために生じるモデル誤差 emod とから成る．
以下では，それらが互いに独立で共に平均ゼロのガウス分布に従うものと仮定し，変換後のデー
タ誤差 e = Reobs + Remod の統計的性質を平均がゼロで分散・共分散行列が σ 2 E のガウス分布(4.12)e ∼ N (0, σ 2 E)で表すことにする．但し，行列 E の構造の具体的設定は，問題によって異なるので，関連論文を参照されたい．
4.2 ベイズモデルと逆化公式
を（4.12）式の誤差分布を仮定して解くには，大別して二つの方法がある．そ観測方程式（4.11）の一つは，誤差の分散・共分散行列が σ 2 I となるように観測方程式を規格化した後，係数行列
し，大きい方から順に q 個（q ≤ p = 係数行列のラを特異値分解（Singular Value Decomposition）
ンク）の固有値とそれに対応する固有ベクトルで疑似逆行列（q = p の場合は Lanczos, 1961, の一般化逆行列）を構成して解を得る方法である（Wiggins, 1972; Jackson, 1972; Hansen, 1987; Xu,．この方法（truncated SVD method）
は，最適な （q 疑似逆行列のランク）を選ぶ規準が曖昧1998）なため，次第に用いられなくなったが，最近，Noda and Matsu’ura（2010）により，ランク q の
疑似逆行列による解は自由パラメータ数 q のモデルの最尤解とみなし得ること，従って，Akaike（1974）の情報量規準（AIC）を用いて最適な q を客観的に決め得ることが示されたので，再び広く用いられるようになるであろう．実際，GPS データから地殻の弾性／非弾性歪み速度分布を推定する問題では，この AIC を用いた打ち切り SVD 法（truncated SVD method using AIC）が有効である．
もう一つは，観測データからの情報に先験的情報を合体させて柔軟な階層構造を持つベイズモデルを構築し，その構造を支配する超パラメータの最適値を Akaike（1980）のベイズ情報量規（尤度）が最大となる解を求める方法である．この方法に準（ABIC）を用いて決め，事後確率密度が Jackson and Matsu’ura（1985）の直接的な先験的情報を取ついては，Matsu’ura et al.（2007）

統計数理

114

第 63 巻

第 1 号 2015

り込んだ手法と Yabuki and Matsu’ura（1992）の間接的な先験的情報を取り込んだ手法を合理的に統合することに成功しているので，それに従って説明する．
先ず，（4.12）式の誤差分布を仮定すれば，観測データ d とモデルパラメータ a を関係付ける
確率モデルは(4.13)

p(d|a; σ 2 ) = (2πσ 2 )−n/2 |E|−1/2 exp[−(d − H a)T E −1 (d − H a)/2σ 2 ]と表現される．ここで，E −1 は n × n 次元の対称行列 E の逆行列，|E| は行列式を表す．地球物理学の多くの問題では，観測データに加え，直接的或いは間接的な先験的情報が存在する．例
えば，GPS データから断層のすべり分布 Δw(η) を推定する問題では，すべりの量や空間的暴れ具合に関する先験的情報に基づき，モデルパラメータ a について直接的或いは間接的拘束条件を課して観測データの解析を行う．
観測データを手にする前のモデルパラメータ a の最も尤もらしい値を ā，その不確かさを fとしたとき，直接的な先験的情報は(4.14)ā = a + f ,f ∼ N (0, ε2 F )と書け，その確率密度関数は(4.15)

r(a; ε2 ) = (2πε2 )−m/2 |F |−1/2 exp[−(a − ā)T F −1 (a − ā)/2ε2 ]と表せる．一方，すべり分布の暴れ具合（roughness）を Δw の 2 次空間導関数の自乗和，(4.16)
R=[(∂ 2 Δw/∂η12 )2 + 2(∂ 2 Δw/∂η1 ∂η2 )2 + (∂ 2 Δw/∂η22 )2 ]dη,Σで測ることにすると，間接的な先験的情報は確率密度関数の形で(4.17)

q(a; ρ2 ) = (2πρ2 )−m/2 |G|1/2 exp[−(a − ā)T G(a − ā)/2ρ2 ]と表せる（行列 G の具体的表現式は問題によって異なるので関連論文を参照）．これらの先験的に従って適切に統合すれば，トータルの先験的情報の確率密度関情報を Fukahata et al.（2004）
数として(4.18) p(a; ρ2 , ε2 ) = (2π)−m/2 |ρ−2 G + ε−2 F −1 |1/2 exp[−(a − ā)T (ρ−2 G + ε−2 F −1 )(a − ā)/2]が得られる．
式の先験的情報をベイズの規則によって (4.13）式の観測データからの情報とそこで，（4.18）
合体させると，以下のような階層構造を持つベイズモデルが構築される：(4.19)

p(a; σ 2 , α2 , β 2 |d) = c(2πσ 2 )−(m+n)/2 |E|−1/2 |α2 G + β 2 F −1 |1/2 exp[−s(a)/2σ 2 ].

但し，
(4.20)

s(a) = (d − H a)T E −1 (d − H a) + (a − ā)T (α2 G + β 2 F −1 )(a − ā).

ここで，α2 = σ 2 /ρ2 と β 2 = σ 2 /ε2 は，ρ2 と ε2 の代わりに導入した超パラメータで，それぞれ，観測データに対する直接的及び間接的な先験的情報の重みを支配する．結局，我々の問題
は，観測データ d が与えられた時のベイズモデルの事後確率密度（尤度）を最大にするモデルパを求めることに帰着する．
ラメータと超パラメータの組（a, σ 2 , α2 , β 2 ）2いま，超パラメータ α と β 2 の値を特定すれば，その時の（4.19）式の事後確率密度を最大にするモデルパラメータの値（最尤解）は次式で与えられる：GPS データの逆解析と地震の発生予測(4.21)

115a∗ = ā + (H T E −1 H + α2 G + β 2 F −1 )−1 H T E −1 (d − H ā).

実際には超パラメータ α2 と β 2 は可変なので，その最適値を求めるには，上記の最尤解を用のベイズ情報量規準（ABIC = −2 log（marginal likelihood）+2（number of adいて Akaike（1980）
を
justable hyper-parameters）(4.22) ABIC(α2 , β 2 |d) = n log s(a∗ )−log |α2 G + β 2 F −1 |+log |H T E −1 H + α2 G + β 2 F −1 |+C と評価し，ABIC を最小にする α2 と β 2 を数値的に探してやれば良い（詳細は Matsu’ura et al.,．
2007 を参照）
こうして求めた超パラメータ α2 及び β 2 の最適値を α̂2 及び β̂ 2 と表記すると，最適モデルとその推定誤差の分散・共分散行列は，σ̂ 2 = s(â)/n として，次式で与えられる：(4.23)â = ā + (H T E −1 H + α̂2 G + β̂ 2 F −1 )−1 H T E −1 (d − H ā),(4.24)C(â) = σ̂ 2 (H T E −1 H + α̂2 G + β̂ 2 F −1 )−1 .

この一般的な逆化公式は，間接的な先験的情報がない場合（α̂2 G → O ），Jackson and Matsu’ura(1985) の公式，
(4.25) â = ā +(H T E −1 H + β̂ 2 F −1 )−1 H T E −1 (d−H ā),C(â) = σ̂ 2 (H T E −1 H + β̂ 2 F −1 )−1 ,に帰着する．また，直接的な先験的情報がない場合（β̂ 2 F −1 → O ）は，Yabuki and Matsu’ura（1992）の公式，(4.26) â = ā + (H T E −1 H + α̂2 G)−1 H T E −1 (d − H ā),C(â) = σ̂ 2 (H T E −1 H + α̂2 G)−1 ,に帰着する．更に，直接的な先験的情報も間接的な先験的情報もない場合（α̂2 G → O 且つ
2 −1
→ O）
，もし解が存在するなら，良く知られた最小自乗解，β̂ F(4.27)â = (H T E −1 H )−1 H T E −1 d,C(â) = σ̂ 2 (H T E −1 H )−1 ,に帰着する．但し，この場合は σ̂ 2 = (d − H â)T E −1 (d − H â)/(n − m) となることに注意．
5.

実際の解析例と考察

最後に，これまで述べてきた逆化手法を実際の GPS データに適用した例を幾つか紹介し，その解析結果について地震の発生予測という観点から考察する．
5.1 新潟 - 神戸変形集中帯の弾性／非弾性歪み解析
Sagiya et al.（2000）
は，Shen et al.（1996）の方法を用いて GEONET（国土地理院の高密度 GPS
観測網）による水平変位速度データ（1997–1999）の 2 次元幾何学的歪み解析を行い，新潟から神戸が存在することを明らかにしている．Nodaに延びる顕著な変形集中帯（新潟 - 神戸変形集中帯）
は，これとほぼ同じ GPS 水平変位速度データに 4.2 節の「モデル選択にand Matsu’ura（2010）
AIC を用いた打ち切り SVD 法」
を適用し，先ず，新潟 - 神戸変形集中帯の 3 次元モーメントレー
ト・テンソル密度分布を求めた（図 2）．次に，その結果を用いて非弾性歪み速度場と全歪み速度
場を独立に計算し，弾性歪みと非弾性歪みの分離が可能であることを示した．彼らの解析結果
は，地表の全歪みのパターン（領域中央部の顕著な体積収縮と剪断変形）
を見る限り，Sagiya etal.（2000）の解析結果と調和的である．しかし，弾性歪みと非弾性歪みに分離して 3 次元的な分布を見ると，領域中央部の体積収縮は弾性的なもので地表付近に限られていること，また，剪

116

統計数理

第 63 巻

第 1 号 2015

図 2．物理モデルを介した逆解析により GPS データから求めたモーメントレート・テンソル
密度分布（Noda and Matsu’ura, 2010）．（a）等方部分．（b）偏差部分．
Fig. 2. The distribution of moment rate tensor density in Central Japan, estimated fromthe 3-D physical inversion analysis of interseismic（1996–2000）GPS data（modiﬁed after Noda and Matsu’ura, 2010）.（a）The distribution of isotopic part at thedepth of 5 km. The red and blue of the color scale indicate expansion and compression, respectively.（b）The distribution of deviatoric part at the depth of 5km. The color scale indicates the magnitude of maximum shear. The red and blue
bars indicate the directions of maximum tension and compression, respectively.

断変形は非弾性的なもので上部地殻全体に及んでいることが分かる（図 3）．このことから，新潟 - 神戸変形集中帯は，地殻内の継続的な非弾性変形運動とそれに伴う地殻応力の再配分が進行している場所と理解することができる．
が GPS データから求めた最大剪断歪み速度の分布
（図 4
（a）
）と尾形
実際，Sagiya et al.（2000）（background（2014）が時空間 ETAS モデルを用いて地震活動データから求めた常時地震活動度
の分布
（図 4
（b）
）を比較すると，両者の間には強い正の相関が見てとれる．これとはseismicity）
独立に地震の CMT データから日本列島域の 3 次元広域応力場のパターンが求められているので，GPS データの弾性／非弾性歪み解析から得られた応力変化が広域応力場と同センスの領域
では大地震の発生確率が高まる（逆センスの領域では低下する）という検証可能な作業仮説を立てることができるであろう．
5.2 プレート間すべり遅れと 2011 年東北沖地震
2011 年 3 月 11 日に東北沖の北米 - 太平洋プレート境界で Mw 9.0 の地震（東北沖地震）が発生した．この地域の地殻変形運動は，1996 年以降現在に至るまで，GEONET によって詳細にモニタリングされてきた．東北沖地震が発生する前の GEONET のデータには，プレート境界の摩擦強度が相対的に大きい部分（アスペリティ）
でのすべり遅れが引き起こす地表変形運動が記録されている．また，東北沖地震の発生を跨ぐ期間の GEONET のデータには，地震時のすべりに起因する地表変形が記録されている．
Hashimoto et al.（2009, 2012）
は，東北沖地震が発生する前の静穏期（1996–2000）の GPS デー
を適用し，タに 4.2 節の「階層的ベイズモデルの構造パラメータの選択に ABIC を用いた最尤法」

GPS データの逆解析と地震の発生予測

117

図 3．深さ 10 km の非弾性，弾性，及び全歪み速度分布（Noda and Matsu’ura, 2010）．上（c）全体積膨張歪み速度分布，下段左から順に段左から順に（a）非弾性，（b）弾性，及び（f）全最大剪断歪み速度分布を示す．
（d）非弾性，（e）弾性，及びFig. 3. Inelastic, elastic and total strain rates in Central Japan（modiﬁed after Nodaand Matsu’ura, 2010）. In the ﬁrst row, we show the distribution of（a）inelastic,（b）elastic and（c）total dilatation rates at the depth of 10 km. The red and blueof the color scale indicate dilatation and contraction, respectively. In the secondrow, we show the distribution of（d）inelastic,（e）elastic and（f）total maximum
shear strain rates at the depth of 10 km. The color scale represents the magnitudeof maximum shear strain rate. The red and blue bars indicate the directions ofmaximum extension and contraction, respectively.

図5（a）に示すような東北・北海道沖のプレート境界のすべり遅れレートの分布を推定している．
の 3 次元プレート境界モデルが用いられている．こちなみに，解析には Hashimoto et al.（2004）の図から，南千島 - 日本海溝に沿って幾つもの顕著なすべり遅れ領域が分布していることが見てとれる．これらのすべり遅れ領域は，20 世紀に発生した M 7.5 以上のプレート境界地震の震源域とほぼ完全に一致する．このことは，そこにアスペリティが存在して，少なくとも 2000 年の時点では，次の大地震に向けた応力蓄積が着々と進行していたことを意味している．実際，十勝沖のすべり遅れ領域では，1952 年の十勝沖地震から約 50 年経った 2003 年に Mw 8.1 の地震が発生している．
東北沖地震の発生によって東日本の地殻は大きく変形した．GEONET が捉えた東北沖地震に伴う地表変位データを，地震間のすべり遅れレートを推定したのと同じプレート境界モデル，同じ手法を用いて逆解析すると，図 5（b）に示すような地震時のすべり分布が得られる（Hashimoto

118

統計数理

第 63 巻

第 1 号 2015

図 4．地殻の歪み速度と常時地震活動度の比較．（a）GPS データの幾何学的解析で求めた最大．（b）時空間 ETAS モデルで求めた常時地震活動度
剪断歪み速度（Sagiya et al., 2000）（尾形, 2014）．
Fig. 4. Comparison of background seismicity with crustal strain rates in Japan.
（a）Maximum shear strain rates estimated from the 2-D geometrical analysis ofinterseismic（1997–1999）GPS data（modiﬁed after Sagiya et al., 2000）. The blackbars indicate the direction of maximum right-lateral shear strain.（b）Backgroundseismicity（1997–2011）estimated with a space-time ETAS model（modiﬁed after
尾形, 2014）. The color scale indicates the number of events（M ≥ 4.0）per 1 × 1
deg2 per day.

et al., 2012）．地震時の顕著なすべりは，図 5（a）の宮城沖及び福島沖の二つのすべり遅れ領域を併せた，南北 400 km，東西 200 km の範囲に及んでいる．すべり量は宮城沖の主破壊領域で最大 25 (±3) m に達する．地震前のすべり遅れのパターンと地震時のすべりのパターンの良い一致は，アスペリティのすべり遅れによって蓄積された応力が地震時のすべりによって解放されるという基本的考え方が，観測的にも支持されることを意味している．問題は宮城沖のすべり遅れ領域で発生した主破壊の規模が大きすぎるという点である．
宮城沖のすべり遅れ領域では，過去 200 年以上に亘り，M 7.5 クラスの地震がほぼ 40 年間隔で繰り返し発生してきた．その同じ領域で今回は M 9 クラスの超巨大地震が発生した．アスペリ
ティを物理的実体とみなす従来の考えでは，プレート境界には地震発生単位としての大小のア
スペリティが互いに重なることなく分布しているとされてきた．この場合，隣接するアスペリ
ティの相互作用による連動破壊は説明できても，同じ場所で異なる規模の地震が異なる時間間隔で繰り返し発生することは説明できない．しかし，アスペリティを断層摩擦特性（臨界強度と臨界変位量）の空間的不均一をスペクトル分解したときの各モードに対応する概念（松浦, 2012）と理解すれば，異なるサイズのアスペリティが同じ場所に階層的に存在しても不思議ではない．
そして，異なる臨界変位量を持つ異なるサイズのアスペリティは異なる時間間隔で異なる規模の地震発生を繰り返すので，全体として見ればスケールに依存するマルチ地震サイクルの様相を呈する．このように考えると，40 年で臨界状態に達する 30 km サイズの局所的アスペリティの下に隠れていた 300 km サイズの基盤アスペリティが，何百年もかけて応力を蓄積し，一挙に破壊したのが今回の超巨大地震ということになる．

GPS データの逆解析と地震の発生予測

119

図 5．GPS データの逆解析で求めた東日本の地震間すべり遅れレート分布と東北沖地震のす．（a）北米 - 太平洋プレート境界のすべり遅れレート
べり分布（Hashimoto et al., 2012）
分布．
（b）東北沖地震のすべり分布．赤丸はプレート境界で発生した M 6.0 以上の余震．
紫色のコンターは余効すべり分布．
Fig. 5. Interseismic slip-deﬁcit rate distribution in East Japan and the slip distributionof the 2011 Tohoku-oki earthquake, estimated from GPS data（modiﬁed after
Hashimoto et al., 2012）.（a）Slip-deﬁcit rate distribution（blue contours）on the
North American-Paciﬁc plate interface, estimated from interseismic（1996–2000）
GPS data.（b）Coseismic slip distribution（red contours）of the 2011 Tohoku-oki
earthquake, estimated from coseismic GPS data. Afterslip distribution is also
shown by the purple contours. The yellow star, white squares, and red circles indicate the epicenters of the main shock, foreshocks, and interplate aftershocks of
the Tohoku-oki earthquake, respectively. The gray contours represent the depth
to the plate interface.

Terakawa and Matsu’ura（2010）
が CMT データの逆解析で求めた東北地域の地殻応力場のパ（a））は，基本的に，プレート境界の摩擦強度（プレート境界に働く剪断応力）に支配
ターン（図 6
されている．2011 年東北沖地震は，プレート境界に働く剪断応力を低下させ，震源周辺域に大（b））．東北沖地震の発生に伴う地震活動の時空間変化を合理
きな応力変化を引き起こした（図 6が指摘しているように，地震時の応力変化だけでは的に説明するには，Terakawa et al.（2013）
不充分で，地殻応力の絶対レベル（プレート境界の摩擦強度）と地殻流体の分布と移動についての新たな知見が不可欠である．ここでは，参考として，Hashimoto et al.（2012）の地震時すべりが計算した，地震時及び地震後の歪みと応力の変化（図 7）
分布に基づいて Hashima et al.（2014）を示すにとどめる．
5.3 大正及び元禄関東地震の震源域と再来間隔
関東地域では，北米プレートの下にフィリピン海プレートが沈み込み，そのフィリピン海
プレートの下に太平洋プレートが沈み込んでいる．Noda et al.（2013）は，この地域の地震間

120

統計数理

第 63 巻

第 1 号 2015

図 6．東北沖地震発生前のテクトニック応力場と地震時の応力変化（Terakawa et al., 2013）．
（a）CMT データの逆解析で求めた東北沖地震発生前のテクトニック応力パターン．
（b）地震時すべり分布（図 5（b））から計算した応力変化．
Fig. 6. Background tectonic stress and coseismic stress changes in East Japan（Terakawa
et al., 2013）.（a）The tectonic stress pattern（at the depth of 10 km）estimatedfrom CMT data of seismic events before the 2011 Tohoku-oki earthquake is represented with the lower hemisphere projection of focal mechanisms of potential
seismic events.（b）The stress changes（at the depth of 10 km）due to the 2011
Tohoku-oki earthquake are computed for the coseismic slip distribution（darkblue contours）in Fig. 5（b）. The pattern of stress changes is represented in thesame way as in（a）. The color scale indicates the magnitude of maximum shear
stress. The star indicates the epicenter of the 2011 Tohoku-oki earthquake.

（1996–2000）の GPS 速度データを歪みに変換し（図 8），それに 4.2 節の「階層的ベイズモデルの構造パラメータの選択に ABIC を用いた最尤法」を適用して，現在のプレート間すべり遅れレー．一方，Matsu’ura et al.（2007）は，1923 年（大正）関東地震に伴う地殻トを推定している（図 9）変動データ（三角測量による水平変位及び水準測量による垂直変位）を逆解析し，地震時のすべり分布を推定している．いずれの場合も解析には同じ 3 次元プレート境界モデル（Hashimoto etを用いているので，推定した地震時すべり分布と地震間すべり遅れレート分布の比較
al., 2004）（関東）地震の震源域と再来間隔について以下のように推論することができる．
から，1703 年元禄（a）に示すように，房総半島南東沖から相模トラフ
北米 - フィリピン海プレート境界では，図 9を経て伊豆半島北端部に至る帯状の領域で顕著なすべり遅れが進行している．このすべり遅れ（a）
レートの分布図に大正関東地震時のすべり分布を黄色のコンターで書き加えたものが図 10である．先ず，伊豆半島北端部は衝突境界なので，そこでのすべり遅れがプレート境界地震によって解放されることはない．一方，大正関東地震の震源域は相模トラフに沿ったすべり遅れ
領域と重なる．また，そこでの平均すべり遅れレート（2 mm/yr）と地震時の平均すべり量（6 m）
から，大正型関東地震の再来間隔は約 300 年と見積もられる．これは宍倉（2003）が岩井低地の離水浜堤群の C14 年代から見積もった大正型関東地震の再来間隔（約 400 年）と概ね調和的であGPS データの逆解析と地震の発生予測

121

図 7．東北沖地震に伴う地震時及び地震後の（a）変位，（b）膨張歪み，（c）最大剪断歪み，及び（Hashima et al., 2014）．いずれも図 5（b）の地震時すべり分布（d）最大剪断応力の変化（b）の A-A に沿った鉛直断面が示されている．
から計算したもので，図 6
Fig. 7. Coseismic and postseismic changes in（a）displacement,（b）dilatation,（c）maximum shear strain, and（d）maximum shear stress of the 2011 Tohoku-oki earthquake（modiﬁed after Hashima et al., 2014）. The vertical sections of displacement,strain and stress ﬁelds along the line A-A in Fig. 6（b）are computed for the coseismic slip distribution in Fig. 5（b）. The thick solid line indicates the plate interface
shallower than 50 km. The thin horizontal line indicates the elastic-viscoelastic
layer interface. The black and white bars in（c）indicate the directions of extension
and contraction, respectively.

る．従って，房総半島南端部に 4 m 以上の隆起を引き起こした 1703 年元禄（関東）地震は，房総半島南東沖のすべり遅れ領域（灰色の楕円）での地震破壊が大正型関東地震の震源域にまで及ん（b））．問題はその再来間隔である．
だと考えるのが合理的であろう（図 10
房総半島南端部には沼 I〜IV 面と呼ばれる完新世海岸段丘が発達している．最低位の沼 IV 面
は元禄地震の際に離水したことが確かめられている．この元禄段丘面と高位の沼 I〜III 面の高度分布パターンが良く似ていることから，従来，沼 I〜III 面も昔の元禄型地震によって離水したと考えられてきた（Matsuda et al., 1978; Shimazaki and Nakata, 1980; 宍倉, 2003）．ちなみに，
沼 I〜IV 面の離水年代から元禄型関東地震の再来間隔を見積もると，平均で約 2000 年となる．
しかし，プレート境界地震の場合，地震時にすべった領域（震源域）
はやがて再固着するが，それ以外のところでは地震間に非地震性すべりが進行するので，地震時の隆起・沈降パターンは
時間と共に徐々に失われていき，最終的に残るのはプレートの定常沈み込みによる変動だけで．従って，沼 I〜III 面の形成は，太平洋プレートとフィリピン
ある（Matsu’ura and Sato, 1989）
海プレートの定常沈み込みによる房総半島南端部の急速な隆起運動（例えば，Hashimoto et al.,と完新世の平均海面変動（例えば，Siddall et al., 2003）に帰すべきもので，元禄型地震の
2004）

122

統計数理

第 63 巻

第 1 号 2015

図 8．関東地域の地震間（1996–2000）の地殻変形速度（Noda et al., 2013）．（a）GPS 水平速度データ．
（b）GPS 水平速度データからプレート定常沈み込みによる分を差し引いて変換した水平歪み速度．
Fig. 8. Interseismic crustal deformation rates in the Kanto region（modiﬁed after Noda et al.,2013）.（a）Horizontal velocities obtained from GPS measurements for 1996–2000（Sagiya,2004）. The gray arrows represent the relative velocity vectors to a reference point YST（solid circle）.（b）The horizontal strain rates transformed from the GPS velocity data after
subtracting theoretical velocities due to steady plate subduction. The gray solid lines show
the optimum triangular mesh constructed from GPS stations with Delaunay triangulation.
The black bar in each triangle indicates the direction of maximum horizontal contraction.

図 9．GPS 水平歪み速度データの逆解析で求めた関東地域の地震間すべり遅れレート（Noda et al.,2013）．（a）北米 - フィリピン海プレート境界のすべり遅れレート分布．（b）フィリピン海 - 太平洋
プレート境界のすべり遅れレート分布．
Fig. 9. Interseismic slip-deﬁcit rate distribution in the Kanto region, inverted from GPS strain
rate data（Noda et al., 2013）.（a）The North American-Philippine Sea plate interface.
（b）The Philippine Sea-Paciﬁc plate interface. The blue and red contours show slip-deﬁcit
and -excess rates at intervals of 10 mm/yr, respectively. The thick arrows indicate slip-rate
vectors on each plate interface.

GPS データの逆解析と地震の発生予測

123

図 10．1923 年関東地震及び 1703 年元禄地震の震源域．
（a）図 9（a）のすべり遅れレート分布（黄色のコンター）を重ねたもに Matsu’ura et al.（2007）による関東地震のすべり分布（黒のコンター）に関東地の．（b）Matsuda et al.（1978）による元禄地震の地表隆起分布
震のすべり分布を重ねたもの．影を付けた部分が元禄地震の推定震源域．
Fig. 10. The source regions of the 1923 Kanto earthquake and the 1703 Genroku earthquake.（a）The coseismic slip distribution of the 1923 Kanto earthquake（yellow
contours）estimated by Matsu’ura et al.（2007）is shown together with the interseismic slip-deﬁcit rate distribution（blue contours）estimated by Noda et al.
（2013）. The 1923 Kanto earthquake has not released the slip deﬁcits southeast
oﬀ the Boso peninsula.（b）A plausible source region of the 1703 Genroku earthquake（shaded region）is shown together with the coseismic surface uplifts（black
contours）estimated by Matsuda et al.（1978）.

発生とは無関係である．そうすると，元禄型地震の再来間隔を見積もるには，既に失われつつある 1703 年元禄地震時の海岸隆起データに頼るしかない．以下は概算であるが，元禄地震時の（宍倉, 2003）とすると，GPS
海岸隆起データを説明する房総半島南東沖の平均すべり量を 12 m
データから推定したそこでの平均すべり遅れレートは 2 mm/yr なので，元禄型地震の再来間隔
は約 600 年と見積もられる．この見積が正しければ，関東地域では，大正型と元禄型のプレー
ト境界地震が 300 年置きに交互に発生することになる．
謝辞

本論文原稿を丁寧にお読み頂き，有益な助言を下さった匿名の査読者に感謝いたします．

参

考

文

献

Akaike, H.（1974）. A new look at the statistical model identiﬁcation, IEEE（Institute of Electrical andElectronics Engineers ）Transactions on Automatic Control, 19, 716–723.
Akaike, H.（1980）. Likelihood and the Bayes procedure, Bayesian Statistics（eds. J. M. Bernardo, M.
H. DeGroot, D. V. Lindley and A. F. M. Smith）, 143–166, University Press, Valencia, Spain.
Backus, G. and Mulcahy, M.（1976）. Moment tensors and other phenomenological descriptions of seismic sources — I. Continuous displacements, Geophysical Journal of the Royal Astronomical

統計数理

124

第 63 巻

第 1 号 2015

Society, 46, 341–361.
Beeler, N. M., Simpson, R. W., Hickman, S. H. and Lockner, D. A.（2000）. Pore ﬂuid pressure, apparent
friction, and Coulomb failure, Journal of Geophysical Research, 105, 25533–25542.
Bosl, W. J. and Nur, A.（2002）. Aftershocks and pore ﬂuid diﬀusion following the 1992 Landers earthquake, Journal of Geophysical Research, 107, doi:10.1029/2001JB000155.
Byerlee, J.（1978）. Friction of rocks, Pure and Applied Geophysics, 116, 615–626.
Carlson, J. M. and Langer, J. S.（1989）. Mechanical model of an earthquake fault, Physical Review A,
40, 6470–6484.
Carlson, J. M., Langer, J. S., Shaw, B. E. and Tang, C.（1991）. Intrinsic properties of Burridge-Knopoﬀ
model of an earthquake fault, Physical Review A, 44, 884–897.
Dieterich, J. H. and Smith, D. E.（2009）. Nonplanar faults: Mechanics of slip and oﬀ-fault damage,
Pure and Applied Geophysics, 166, 1799–1815.
Freund, L. B.（1972）. Energy ﬂux into the tip of an extending crack in an elastic solid, Journal of
Elasticity, 2, 341–349.
Fukahata, Y. and Matsu’ura, M.（2006）. Quasi-static internal deformation due to a dislocation source ina multilayered elastic/viscoelastic half-space and an equivalence theorem, Geophysical JournalInternational, 166, 418–434.
Fukahata, Y., Nishitani, A. and Matsu’ura, M.（2004）. Geodetic data inversion using ABIC to estimate
slip history during one earthquake cycle with viscoelastic slip-response functions, GeophysicalJournal International, 156, 140–153.
Fung, Y. C.（1965）. Foundations of Solid Mechanics, Prentice-Hall, Englewood Cliﬀs, New Jersey.
Hansen, P. C.（1987）. The truncated SVD as a method for regularization, BIT Numerical Mathematics,
27, 534–553.
Hashima, A., Fukahata, Y., Hashimoto, C. and Matsu’ura, M.（2014）. Quasi-static strain and stressﬁelds due to a moment tensor in elastic-viscoelastic layered half-space, Pure and Applied Geophysics, 171, 1669–1693.
Hashimoto, C. and Matsu’ura, M.（2002）. 3-D simulation of earthquake generation cycles and evolutionof fault constitutive properties, Pure and Applied Geophysics, 159, 2175–2199.
Hashimoto, C., Fukui, K. and Matsu’ura, M.（2004）. 3-D modelling of plate interfaces and numericalsimulation of long-term crustal deformation in and around Japan, Pure and Applied Geophysics,
161, 2053–2068.
Hashimoto, C., Noda, A., Sagiya, T. and Matsu’ura, M.（2009）. Interplate seismogenic zones along the
Kuril-Japan trench inferred from GPS data inversion, Nature Geoscience, 2, 141–144.
Hashimoto, C., Noda, A. and Matsu’ura, M.（2012）. The Mw9.0 northeast Japan earthquake: Total
rupture of a basement asperity, Geophysical Journal International, 189, 1–5.
Hashimoto, C., Fukuyama, E. and Matsu’ura, M.（2014）. Physics-based 3-D simulation of earthquake
generation cycles at plate interfaces in subduction zones, Pure and Applied Geophysics, 171,1705–1728.
Jackson, D. D.（1972）. Interpretation of inaccurate, insuﬃcient and inconsistent data, GeophysicalJournal of the Royal Astronomical Society, 28, 97–109.
Jackson, D. D. and Matsu’ura, M.（1985）. A Bayesian approach to nonlinear inversion, Journal ofGeophysical Research, 90, 581–591.
Knopoﬀ, L.（1958）. Energy release in earthquakes, Geophysical Journal of the Royal Astronomical
Society, 1, 44–52.
Kostrov, B. V.（1966）. Unsteady propagation of longitudinal shear cracks, Journal of Applied Mathematics and Mechanics, 30, 1241–1248.
Lanczos, C.（1961）. Linear Diﬀerential Operators, Van Nostrand, London.

GPS データの逆解析と地震の発生予測

125

Matsuda, T., Ota, Y., Ando, M. and Yonekura, N.（1978）. Fault mechanism and recurrence time ofmajor earthquakes in southern Kanto district, Japan, as deduced from coastal terrace data,
Geological Society of America Bulletin, 89, 1610–1618.
松浦充宏（2012）. 東北沖超巨大地震とプレート沈み込み帯のマルチ地震サイクル, 地質学雑誌, 118,
313–322.
Matsu’ura, M. and Sato, T.（1989）. A dislocation model for the earthquake cycle at convergent plateboundaries, Geophysical Journal International, 96, 23–32.
Matsu’ura, M., Kataoka, H. and Shibazaki, B.（1992）. Slip-dependent friction law and nucleation processes in earthquake rupture, Tectonophysics, 211, 135–148.
Matsu’ura, M., Noda, A. and Fukahata, Y.（2007）. Geodetic data inversion based on Bayesian formulation with direct and indirect prior information, Geophysical Journal International, 171,1342–1351.
Matsu’ura, M., Noda, A., Terakawa, T. and Fukahata, Y.（2014）. Tectonic stress ﬁelds in subduction
zones governed by frictional strength of plate interfaces, Abstracts of the Japan Geoscience
Union 2014 Annual Meeting, SCG66–04.
Nakanishi, H.（1992）. Earthquake dynamics driven by a viscous ﬂuid, Physical Review A, 46, 4689–
4692.
Noda, A. and Matsu’ura, M.（2010）. Physics-based GPS data inversion to estimate 3-D elastic andinelastic strain ﬁelds, Geophysical Journal International, 182, 513–530.
Noda, A., Hashimoto, C., Fukahata, Y. and Matsu’ura, M.（2013）. Interseismic GPS strain data inversion to estimate slip-deﬁcit rates at plate interfaces: Application to the Kanto region, central
Japan, Geophysical Journal International, 193, 61–77.
Nur, A. and Booker, J. R.（1972）. Aftershocks caused by pore ﬂuid ﬂow?, Science, 175, 885–887.
Ogata, Y.（1988）. Statistical models for earthquake occurrences and residual analysis for point processes, Journal of the American Statistical Association, 83, 9–27.
Ogata, Y.（1998）. Space-time point-process models for earthquake occurrences, Annals of the Instituteof Statistical Mathematics, 50, 379–402.
Ogata, Y.（2004）. Space-time model for regional seismicity and detection of crustal stress changes,
Journal of Geophysical Research, 109, B03308, doi:10.1029/2003JB002621.

尾形良彦（2014）. 大地震の確率予測と統計モデル, 地震予知連絡会会報, 91, 387–395.
Omori, F.（1894）. On the aftershocks of earthquakes, Journal of College Science, Imperial Universityof Tokyo, 7, 111–200.
Reasenberg, P. A. and Simpson, R. W.（1992）. Response of regional seismicity to the static stresschange produced by the Lorna Prieta earthquake, Science, 255, 1687–1690.
Reid, H. F.（1910）. The Mechanics of the Earthquake, The California Earthquake of April 18, 1906
（Report of the State Earthquake Investigation Commission ）, Vol. 2, Carnegie Institution of
Washington, Washington, D. C.
Sagiya, T.（2004）. Interplate coupling in the Kanto district, central Japan, and the Boso peninsulasilent earthquake in May 1996, Pure and Applied Geophysics, 161, 2327–2342.
Sagiya, T., Miyazaki, S. and Tada, T.（2000）. Continuous GPS array and present-day crustal deformation of Japan, Pure and Applied Geophysics, 157, 2303–2322.
Savage, J. C.（1969）. Steketee’s paradox, Bulletin of the Seismological Society of America, 59, 381–384.
Shen, Z.-K., Jackson, D. D. and Ge, B. X.（1996）. Crustal deformation across and beyond the Los Angeles basin from geodetic measurements, Journal of Geophysical Research, 101, 27957–27980.
Shibazaki, B. and Matsu’ura, M.（1992）. Spontaneous processes for nucleation, dynamic propagation,and stop of earthquake rupture, Geophysical Research Letters, 19, 1189–1192.
Shimazaki, K. and Nakata, T.（1980）. Time-predictable recurrence model for large earthquakes, Geo-

126

統計数理

第 63 巻

第 1 号 2015

physical Research Letters, 7, 279–282.
宍倉正展（2003）. 変動地形から見た相模トラフにおけるプレート間地震サイクル, 地震研究所彙報, 78,
245–254.
Siddall, M., Rohling, E. J., Almogi-Labin, A., Hemleben, Ch., Meischner, D., Schmelzer, I. and Smeed,
D. A.（2003）. Sea-level ﬂuctuations during the last glacial cycle, Nature, 423, 853–858.
Stein, R. S., King, G. C. P. and Lin, J.（1992）. Change in failure stress on the southern San Andreas
fault system caused by the 1992 magnitude = 7.4 Landers earthquake, Science, 258, 1328–1332.
Terakawa, T.（2014）. Evolution of pore ﬂuid pressures in a stimulated geothermal reservoir inferredfrom earthquake focal mechanisms, Geophysical Research Letters, 41, 7468–7476.
Terakawa, T. and Matsu’ura, M.（2008）. CMT data inversion using a Bayesian information criterion
to estimate seismogenic stress ﬁelds, Geophysical Journal International, 172, 674–685.
Terakawa, T. and Matsu’ura, M. （2010）. The 3-D tectonic stress ﬁelds in and around Japan
inverted from centroid moment tensor data of seismic events, Tectonics, 29, TC6008,
doi:10.1029/2009TC002626.
Terakawa, T., Hashimoto, C. and Matsu’ura, M.（2013）. Changes in seismic activity following the 2011
Tohoku-oki earthquake: Eﬀects of pore ﬂuid pressure, Earth and Planetary Science Letters,
365, 17–24.
Tsuboi, C.（1932）. Investigation on the deformation of the Earth’s crust in Idu Peninsula connected with
the Idu earthquake of Nov. 26, 1930, Bulletin of the Earthquake Research Institute, Universityof Tokyo, 10, 435–448.
Twiss, R. J. and Unruh, J. R.（1998）. Analysis of fault slip inversions: Do they constrain stress or strain
rate?, Journal of Geophysical Research, 103, 12205–12222.
Utsu, T.（1961）. A statistical study on the occurrence of aftershocks, Geophysical Magazine, 30, 521–
605.
Wiggins, R. A.（1972）. The general linear inverse problem: Implication of surface waves and free oscillations for Earth structure, Reviews of Geophysics and Space Physics, 10, 251–285.
Xu, P.（1998）. Truncated SVD methods for discrete linear ill-posed problems, Geophysical JournalInternational, 135, 505–514.
Yabuki, T. and Matsu’ura, M.（1992）. Geodetic data inversion using a Bayesian information criterion
for spatial distribution of fault slip, Geophysical Journal International, 109, 363–375.

Proceedings of the Institute of Statistical Mathematics Vol. 63, No. 1, 105–127 (2015)

127

Inversion Analysis of GPS Data and Forecast of Earthquake Occurrence
Mitsuhiro Matsu’ura1 and Akemi Noda2
1 Institute of Statistical Mathematics
2 Disaster Reduction Consulting Department, Kozo Keikaku Engineering Inc.

The occurrence of earthquakes is brittle fracture of rocks, and so its deterministic
prediction is diﬃcult. Since the 1890s, various statistical models have been proposed forprobabilistic prediction of future earthquakes from past seismicity data. In contrast, in the
1990s, a new research ﬁeld that aims to understand earthquake generation on the basisof physics was born, and nowadays physics-based earthquake generation models governed
by various fault constitutive（friction）laws have been proposed. In these physical models,the occurrence of earthquakes is described as the release process of shear stress acting
on faults. Therefore, the physical model and the statistical model must be closely relatedwith each other. However, it is not easy to determine the actual stress state of the earth’s
crust. At present, we can only estimate stress patterns from inversion analysis of focal
mechanism solutions or CMT（Centroid Moment Tensor）solutions of seismic events andstress changes from inversion analysis of GPS（Global Positioning System）data. This article introduces the basic idea and method of GPS data inversion, and examine inversion
results from the viewpoint of earthquake prediction.

Key words: Earthquake generation, statistical model, physical model, crustal stress, GPS data, inversion analysis.