統計数理（2021）
第 69 巻 第 2 号 209–222©2021 統計数理研究所

特集「Hawkes 過程の新展開と応用」［総合報告］

地震学における非線形 Hawkes 過程：
摩擦構成則に基づく地震活動モデル
岩田 貴樹1,2（受付 2021 年 2 月 18 日；改訂 10 月 6 日；採択 10 月 6 日）

要

旨

本稿では摩擦構成則に基づく地震活動モデルについて紹介する．これは非線形 Hawkes 過程に相当するものであり，よく知られている ETAS モデルが線形 Hawkes モデルであるのと対照的である．これ以外にも両者には対照的な特徴があり，それらの比較について述べる．また，実際の地震発生時系列解析を介し，このモデルの利点や問題点，そして今後の展開について示す．
キーワード：地震活動，非線形 Hawkes 過程，速度・状態依存摩擦構成則，点過程解
析，応力，Dieterich モデル．
1.

はじめに

地震学における点過程モデルとして最も良く知られているものの 1 つは Ogata（1988）に よ る epidemic-type aftershock sequence（ETAS）モ デ ル で あ る ．こ れ は 地 震 の 発 生 時 刻
tn = {t1 , t2 , . . . , tn } と規模（マグニチュード）mn = {m1 , m2 , . . . , mn } で構成されるマーク付き点過程に対し，時刻 t の直前までの地震発生時刻とマークの履歴（Ht で表すこととする）が与えられた下での条件付き強度関数（以降，簡略化のため単に強度関数と呼ぶ）
を(1.1)

λ(t|Ht ) = μ +φ(mi )g(t − ti )i;ti <tとするものである．なお，ここではマーク（マグニチュード）の値に関する確率分布は問わないものとしている．これは自己励起過程（self-exciting process）である Hawkes 過程（Hawkes,
をマーク付き点過程に拡張したもの（例えば近江・野村, 2019, 6.3 節）に相当する．なお，
1971）
地震学では「励起」
を「トリガリング」
や「誘発」と呼ぶことが多く，本稿でも必要に応じてこの表現を用いる．
各々の地震が引き起こす地震（余震）の規模依存性と活動減衰それぞれについて，ETAS モデ
ルにおいては経験則に基づき(1.2)

φ(m) = eαm(1.3)g(t) = K/(t + c)p1

県立広島大学 大学教育実践センター：〒 723–0023 広島県庄原市七塚町 5562

2

統計数理研究所 客員：〒 190–8562 東京都立川市緑町 10–3

210

統計数理

第 69 巻 第 2 号 2021とする．後者は地震学では「大森・宇津公式」（Utsu, 1961）と呼ばれるものである（ETAS モデ．
ルの詳細については，例えば本特集号の 庄・尾形, 2021 を参照されたい）
ETAS モデルは実際の地震活動によく適合しており，時空間へ拡張したもの（Ogata, 1998）と合わせ幅広い研究で用いられ，標準的な地震活動モデルと考えられている．但し，純粋に経験則から導かれたものであるため，モデルに含まれているパラメータの値をどのようにして地球物理学的知見と結びつけるかや，逆に何らかの物理的な情報をモデルにどう反映させるかという点に課題を抱えている．
ETAS モデルと並び，地震学においてよく知られ用いられるものとして Dieterich（1994）にとよる地震活動モデルがある（以降，地震学でしばしば使われる名称に従い「Dieterich モデル」
呼ぶ）．これは室内における岩石すべり実験から得られた摩擦の特性を定式化した「摩擦構成則」を元に，応力変動の値と地震活動度の変化を定量化したものである．実験結果として得られた摩擦特性は経験的なものであるが，そこから導かれた摩擦構成則についてはその物理機構の解明がある程度進められており，Dieterich モデルは経験モデルと物理モデルの中間的なものと言えよう．
これ以外にも Dieterich モデルは ETAS モデルと対照的な点があり，このモデルを紹介する
ことが本稿の目的である．なお，著者は以前 Dieterich モデルと ETAS モデルに関する解説（岩を著したことがある．本稿の内容の一部はそれと重なる点があることを，あらかじめ田, 2009）
お断りする．
2.

摩擦構成則に基づく地震活動モデル（Dieterich モデル）

2.1 Dieterich モデルの概要
Dieterich モデルの基礎を成すのは，速度・状態依存摩擦構成則（Dieterich, 1979, 1986, 1987,
から与えられる次の式である．
1992）



(2.1)

dγ =

1
dt − γdτ + γ
Aσ





τ
− α dσσ

ここで，γ は「状態変数」と呼ばれ，後述する通り地震活動度，即ち強度関数はこの γ の逆数に
（図 1）．A と α は共に速度・状態依比例する．τ と σ は，それぞれ剪断応力と法線応力である
存摩擦構成則に現れるパラメータである．A は状態変数 γ の法線応力 σ への依存性を定量的に表すために導入されたものである．α は地震を起こすすべり面（断層面）の状態を反映したものと物理的には解釈される．
式（2.1）に基づいて状態変数の時間発展 γ(t) を求める．モデルの簡略化のため，殆んどの既
時間一定の剪断応力変化があった場合と（ii）瞬間的な（step 的な）剪断応
存研究においては，（i）
力の変化の 2 つの状況のみを考える（図 2）．地震が起きる大きな要因の 1 つはプレートの運動により生じる応力の変化（増加）である．プレートはほぼ一定の速度で移動しており，これに伴もほぼ時間一定と考えられている．（i）はこれに対応するものであう応力増加（tectonic loading）
はある地震が起きた時，その断層運動により周辺に及ぼる（図 2 の “secular stress rate”）．（ii）
す応力変化に相当する（同図の “Sudden Stress Increase/Decrease”）．揺れを引き起こすような高速の断層運動に要する時間は，その地震のマグニチュード（m）に依存する．小地震，例えばm が 3 程度のものであれば 1 秒未満（例えば宇津, 2001, 10.5 節），2011 年東北地方太平洋沖地
震のような m が 9 を超える超巨大地震でも数分程度であり（八木, 2012 のレビューに引用されている論文を参照），地震活動解析においては瞬間的なものと見做して問題ない．
（i）の場合について示す．まず時間一定の剪断応力式（2.1）から導かれる具体的な解を，まず変化を τ̇ (= dτ /dt) と表すこととする．この時，計算の簡便化のため法線応力 σ は変化しない，

地震学における非線形 Hawkes 過程：摩擦構成則に基づく地震活動モデル

図 1．剪断応力と法線応力の概念図．（a）剪断応力は断層（太線，実際は紙面鉛直方向に広がっ（あている面）にずれを生じさせるような応力であり，（b）法線応力は断層を押しつける
るいは引っ張る）ような応力である．
Schematic diagram showing shear stress and normal stress. (a) Shear stress causes
dislocation along a fault (bold line), which is stretched along the direction perpendicular to the sheet of paper. (b) Normal stress is compressive (or extensional)
stress acting on the fault.

図 2．簡略化された応力の時間変化の概念図．時間一定の剪断応力変化（“secular stressrate”）と，ある地震が起きた際にその断層運動により周辺に及ぼす応力変化（“Sudden
Stress Increase/Decrease”）の 2 つの要因により応力（“Coulomb Stress”）が時間変化
し，それがある値（“FAILURE STRESS”）に達すると地震が発生する．この図ではある 1 つの地震の発生のみを考えているため，瞬間的な応力増加あるいは減少によりその地震の発生が早まったり遅くなったりする（“Time Advance/Delay”）．これを多数の地震に対して考えることで，強度関数の増加あるいは減少が導かれる（Toda et al.,1998 より引用）．
Schematic diagram showing the simpliﬁed temporal evolution of shear stress(“Coulomb Stress”). Two factors are considered: the one is stress change of whichrate is constant in time (“secular stress rate”), and the other is static stress change
caused by co-seismic dislocations (“Sudden Stress Increase/Decrease”）. In thisﬁgure, only the occurrence of a particular earthquake is considered and it is advanced or delayed by “Sudden Stress Increase/Decrease.” The increase or decreaseof an intensity function is derived from the application of this concept to a set ofmany earthquakes. This ﬁgure is taken from Toda et al. (1998) and is producedwith the permission of American Geophysical Union. Copyright 1998 American
Geophysical Union.

即ち dσ/dt = 0 を合わせて仮定すると，式（2.1）は(2.2)となる．これよりdγ =1[1 − γ τ̇ ]dt
Aσ

211

212

統計数理

(2.3)

γ(t) = γt0 −

第 69 巻 第 2 号 2021







1
τ̇
1
exp −
(t − t0 ) +
τ̇
Aσ
τ̇

を得る．ここで，γt0 は t = t0 における γ の値
（即ち γ(t0 )）である．
次に（ii）についてである．ここでは，t = tb に瞬間的な剪断応力の変化 Δτ があったとする．
（i）の場合と同様，この場合も法線応力 σ は変化し「瞬間的」であるので dt = 0 であり，加えてはないと仮定する．すると式（2.1）
1dγ = −(2.4)γdτAσであり，これより
(2.5)


γ(t+
b ) = γ(tb ) exp −

Δτ
Aσ



となる．ここで γ(t+b ) は時刻 t = tb に起きた剪断応力変化の直後における γ の値である．
1 節に記した ETAS モデルの場合と同様，地震の発生時刻 tn = {t1 , t2 , . . . , tn } とし，マグニチュードに換えて各地震が引き起こす剪断応力の変化 Δτn = {Δτ1 , Δτ2 , . . . , Δτn } をマークとするマーク付き点過程を考える．また，解析対象期間を [S, T ] とし，簡便化のため S = t0 ，T = tn+1 とする．
（2.5）
をまとめることで ti−1 < t ≤ ti (i = 1, 2, . . . , n + 1) に対して
式（2.3）と式

(2.6)γ(t) = γ(ti−1 ) exp −

Δτi−1
Aσ



−





1
τ̇
1
exp −
(t − ti−1 ) +
τ̇
Aσ
τ̇

として，解析対象期間全体に亘る γ(t) を得ることが出来る．
この γ(t) に対して，条件付き強度関数は(2.7)

λ(t|Ht ) = r/γ(t)として与えられる．ここで r は定常状態の地震活動度（長期に亘る地震活動の，単位時間あたりの平均的な地震数）である．
なお，上では「マークとして剪断応力変化を使う」旨を記した．しかし，現実の破壊現象を支配するものとして剪断応力と法線応力を組み合わせたクーロン破壊関数（Coulomb failure（例えば Jeager et al., 2007）がある．特に King et al.（1994）
以降，本震による CFS
stress; CFS）の変化量の空間分布と余震の空間分布との比較が盛んになされるようになったこともあって，
Dieterich モデルを用いる際も剪断応力に換えて CFS を用いることが伝統的に行なわれている．
以上に示したように，Dieterich モデルにおけるパラメータは何らかの物理的な意味を持ったものであり，ETAS モデルの場合に比べ，得られたパラメータの値の解釈を行ない易い．また，マークとして用いられるものは剪断応力変化であり，このような直接的な物理量をモデルに取り込むことが可能である．こういった利点を生かして，応力変化量から期待される地震活動と実際のものとの定量的比較や，逆に地震活動から応力変化量を推定する試みが 1990 年代後半．
より盛んに行なわれている（Dieterich, 2015 のレビューに詳しい）「負の励起」（図 2（b）参照）を生じさせることも容易であさらに Δτ に負の値を与えることでり，実際に大地震の発生後，地震活動が抑制（静穏化）された事例への解析が試みられている（例えば Toda et al., 2012）．ETAS モデルにおいても負の励起を作り出すことは原理的には可能であるが，現状のものを拡張する必要がある（例えば式（1.3）のパラメータ K を正の励起と負の励起が生じた場合で別々の値とする）．また，強度関数の値が負にならないよう何らかの制約を与える必要もあり，現実的な解析においては簡単ではない．

地震学における非線形 Hawkes 過程：摩擦構成則に基づく地震活動モデル

213

2.2 Hawkes 過程としての Dieterich モデル
より一般化された Hawkes 過程は，強度関数が次のような式で表されるものである（例えば，
Daley and Vere-Jones, 2003, Section 7.3）
．



(2.8)



t

h(t − s)N (ds)λ(t|Ht ) = Φ μ +
S

=Φ

μ+



h(t − ti )ti <t

ここで Φ(·) は非負の値を取る関数である．また，Φ(·) が非線形である時，これを（線形の場合と区別するために）非線形 Hawkes 過程と呼ぶことがある．
では，再帰的に用いることにより Dieterich モデルの状態変数 γ(t) を表した．
さて，式
（2.6）
これは Γ(t) = τ̇ γ(t) として，以下のように書き換えることが出来る．
(2.9)

Γ(t) =j

[exp(−Δτi ) − 1] exp

i=0

−



j


Δτk

exp −

k=i+1

t − tita

ここで j は t < ti を満たす最小の i の値である．
（2.8）におまた式（2.7）から λ(t|Ht ) = rτ̇ /Γ(t) である．以上を踏まえると，この強度関数は式
いてφ(Δτi ) = [exp(−Δτi ) − 1] exp

−

j


Δτk

k=i+1

(2.10)

h(t) = exp(−t/ta )Φ(x) = rτ̇ /xとした場合の非線形 Hawkes 過程になっている．これに対し ETAS モデルは Φ(x) = x であり，
線形 Hawkes 過程である．即ち，非線形／線形という点でも ETAS モデルと Dieterich モデルは対照的な関係にある．
2.3 Dieterich モデルの問題点）に示した通り，余震活動の減衰は，時間経過に対するベキ乗則に従大森・宇津公式（式（1.3）
うことが経験的に知られている．この式におけるベキ指数 p は 1 前後の値，具体的には 0.9 から 1.5 程度と余震系列によって多様である（宇津, 1999, 7.3.2 節）．
等に示した Dieterich モデルによる地震活動の減衰は，漸近的に大森・宇津公
一方，式
（2.3）
（図 3）．即ち，Dieterich モデ
式における p = 1 の場合にのみ対応していることが分かっている
ルは p = 1 の場合には対応しておらず，余震減衰の多様性を表現しきれないという問題がある．
この問題は Dieterich（1994）自身が既に指摘しており，同論文では 1 つの解決として複数のstep 的な応力変化によるトリガリング効果により，p の値を見かけ上 1 以外のもの（特に p > 1の場合）とし得ることを指摘している（図 4）．他に Dieterich（1994）
は式（2.3）の導出時に仮定さに換えて「時間変化する τ̇ 」を仮定することで p = 1 となる地震活動のれている「時間一定の τ̇ 」
．さらに地震発生時の減衰を作り得ることも述べている（これについては 4 節で改めて触れる）
における Δτ ）の空間的な不均質性を考慮することも指摘しているが，この場応力変化（式（2.5）
合は p < 1 となる減衰は作り出せる一方，p > 1 となる減衰には対応出来ない．Borovkov andや Helmstetter and Shaw（2009）をはじめとする Dieterich モデル拡張の試み
Bebbington（2003）

214

統計数理

第 69 巻 第 2 号 2021

図 3．式（2.3）から得られる Dieterich モデルにおける余震活動の時間減衰の様子．減衰の時（強度間変動に影響する τ̇ /Aσ の値を変えた例について示す．なお，縦軸は地震活動度関数）を定常状態の地震活動度 r で割ったものである．
Decays of aftershock activity derived from eq. (2.3) of the Dieterich model. Decay
patterns with various values of τ̇ /Aσ, which dominates the pattern, are shown.
The vertical axis indicates the seismicity rate (intensity function) normalized by
the steady (reference) seismicity rate r.

図 4．Dieterich モデルにおいて，step 的な応力変動によるトリガリング効果により，大森・
宇津公式の p の値を変化させる様子を示す概念図．地震活動が減衰する途中で step 的な応力変動がなかった場合（破線．p = 1 に相当）に対し，途中で step 的な応力変動が（実線）．Dieterich（1994）に加筆．
あると見掛け上 p の値を大きく出来る
A schematic diagram showing that the triggering eﬀect caused by stepwise stresschange makes the value of p of the Omori-Utsu formula variable. Without any
stepwise stress change, the decay of the seismic activity (earthquake rate) approximately follows the Omori-Utsu formula with p = 1 (dashed line). Contrastively,
several stepwise stress changes can produce the decay with p > 1 (solid line). Thisﬁgure is modiﬁed from Dieterich (1994) and is produced with the permission ofAmerican Geophysical Union. Copyright 1994 American Geophysical Union.

も同様の問題を抱えている．
「全てまた，ここまでに示した通り，ETAS モデルと同じように Dieterich モデルにおいてもの地震が自己励起（誘発）能力を持つ」と仮定されている．そして，2.1 節に記したように，各々地震学における非線形 Hawkes 過程：摩擦構成則に基づく地震活動モデル

215の地震の自己励起能力，即ち余震を引き起こす能力の多寡を定量化する際に必要となるマークに，ETAS モデルはマグニチュードを用いるのに対し，Dieterich モデルは各地震が周囲に与える瞬間的な剪断応力変化量を用いている．この点は物理量を直接取り込めるという利点となる
一方で，実データに適用する際には問題となる．これは以下に述べる通り，この応力変化量の計算が地震断層同士の位置関係などに依存し複雑であるため，容易ではないからである．
地震学で通常用いられる仮定の範囲内で，断層運動などに伴う応力変化量を計算する公開プ
ログラムは以前から存在し（Okada, 1992），計算そのものは可能である．但し応力は 2 階のテンソルであるため，その成分を適切に座標変換することで剪断応力変化量という特定の応力成分を求めることが出来る．この座標変換の際には，余震を引き起こす側（応力を与える側）・引き起こされる側（応力を与えられる側）
双方の地震断層面の相対的な位置と向き，そして地震時に生じるすべりの向きが必要となる．即ち，これらが不明であれば剪断応力変化量の値を求めることが出来ない．
こういった情報は地震波から求められる「メカニズム解」
より得ることが出来る．但し，地震観測点が密に設置されている好条件下にあっても，特に小地震に対しては，精度良くメカニズ
ム解を求めることは必ずしも容易ではない（例えば Yukutake et al., 2020）．こういった理由により，Dieterich モデルを実際の地震活動に適用する際には，本震・余震系列であれば本震に相当する最大の m を持つ地震のみ，あるいは本震および本震と同程度の大きな m である限られた個数の地震による応力変化の影響だけを考慮してモデリングする（例えば Toda and Stein,
2003）
ことが長年行なわれてきた．即ち，全ての地震による誘発性を考慮できず，現実的なモ
デリングの面で ETAS モデルに劣るものとなっていた．
3.

Dieterich モデルの実データへの適用例：Iwata（2016）

ここでは 2.3 節で示した Dieterich モデルの問題点に対する改善を試み，実際の地震活動へ適を紹介する．
用した一例として，Iwata（2016）（本震・余震系列であれば余震）のトリガリング効果しか取り入れ問題の 1 つは限られた地震
られていないという点にあった．これに対して，Iwata（2016）では，ある地震による応力変化とその地震の m との関係が Δτi = a exp(dmi ) として表されるという仮定を置き，ETAS モデルのように全ての地震によるトリガリング効果の影響を取り入れる工夫を施した．
Iwata（2016）
で扱ったデータは 2004 年新潟県中越地震と 1995 年兵庫県南部地震の，本震発
生後 0.5 日から 4 年までの余震系列である．本震発生直後のデータを除いてあるのは，一般
的に本震発生直後の地震検知能力が低く，データの質に問題があるためである（例えば Nanjo．これらのデータに対し，全ての地震のトリガリング効果の影響を入et al., 2007; Iwata, 2008）と本震のみのトリガリング効果の影響を入れたモデル（モデル B）の 2 つれたモデル（モデル A）
を適用した．モデルパラメータは最尤法により求め，赤池情報量規準（AIC）（Akaike, 1974）によるモデル比較を行なった．表 1 に示す通り，モデル A の AIC の値はモデル B のものに比べて数百あるいは数十小さく，本震だけでなく全ての地震のトリガリング効果を取り入れることでモデルが大きく改善されたことが分かる．
また最尤モデルから得られた強度関数を図 5 に示す．同図には余震活動の大局的な時間変動を示すため大森・宇津公式（本震のみのトリガリング効果を考慮した場合に相当する）
による強度関数も合わせて示してある．大森・宇津公式のベキ指数 p の最尤推定値は，2004 年新潟県中越地震と 1995 年兵庫県南部地震それぞれに対して 1.35 と 1.22 である．
これに対し，モデル B の強度関数は傾きが緩やかであり t−1 とほぼ平行である．これは 2.3
節に示した通り，Dieterich モデルによる地震活動の減衰が大森・宇津公式における p = 1 の場

216

統計数理

第 69 巻 第 2 号 2021

表 1．複数のタイプの Dieterich モデル（モデル A–C）および ETAS モデルを 2 つの余震系列に当てはめて得られた赤池情報量規準（AIC）の値．
Akaike’s Information Criterion (AIC) values of the several types of the Dieterichmodel and ETAS model applied to the two aftershock sequences.

図 5．複数のタイプの Dieterich モデルと大森・宇津公式に関して推定された強度関数．
（b）1995 年兵庫県南部地震に対するもの．Iwata（a）2004 年新潟県中越地震および（2016）を改変．
Estimated intensity functions of the several types of the Dieterich model andOmori-Utsu formula for (a) the 2004 Mid-Niigata (Chuuetsu) earthquake and(b) the 1995 Kobe earthquake. This ﬁgure is modiﬁed from Iwata (2016).

合にのみ対応していることに依るものであり，実際の余震活動の時間減衰に対応し切れていない様子を示している．
（図 4），全ての余震によるトリガ
一方，モデル A では Dieterich（1994）が示唆していた通り
リング効果を取り入れることで，（見かけ上）p の値を変化させることが出来るようになった．
このことも AIC の値を改善させた一因となっている．但し，新潟県中越地震については，モデ
ル A の強度関数は大森・宇津公式のものに追随している一方，兵庫県南部地震についてはよく見ると 10〜100 日目以前は，モデル A の強度関数は大森・宇津公式のものより全般的に値が小さく，それ以後は値が大きくなっている．即ち，図 5 のような両対数グラフで言えば全体的な傾きが緩やかになっており，p の値が 1 よりは大きいものの大森・宇津公式による p = 1.22 ほ
どまでには大きく出来なかったことを示唆している．
新潟県中越地震の余震系列においては，解析期間中，本震の m との差が 1 以下の地震が 3 つ．言い換えると，余震系列全起きている（本震の m = 6.8 に対し，発生した順に 5.8，6.1，5.9）
体の中で本震以外の地震によって誘発された地震が占める割合が比較的高い．これに対し，兵地震学における非線形 Hawkes 過程：摩擦構成則に基づく地震活動モデル

217

庫県南部地震においては，本震の m = 7.3 に対し，解析期間中での最大余震 m = 5.1 と，m の差は 2 以上である．このため，本震による誘発の影響が大きい．この違いが，図 4 に示したよ
うな状況の実現し易さに現れ，新潟県中越地震では十分に p を大きくすることが出来た一方，
兵庫県南部地震ではいくらかは p を大きく出来たものの十分と言えるほどには至らなかったと考えられる．
また，こうした改善を行なってもやはり表 1 に示されている通り，Dieterich モデルの適合性
は，ETAS モデルのものに比べてまだ劣っている．このことは，前段落に述べた通り大森・宇津公式から得られた p の値に対応し切れていない兵庫県南部地震に対しては当然の結果と言える．しかし，対応している（ように見える）
新潟県中越地震に対しても同様であり，これには別の要因を考える必要がある．この要因の 1 つとしては，Dieterich モデルと ETAS モデルにおける各々の地震によるトリガリングの効果の現れ方の違いが挙げられる．
「和」，
式（1.1）に示した通り，ETAS モデルにおいては各々の地震によるトリガリング効果が即ち加法的に現れる．これに対して，Dieterich モデルにおいては式（2.5）あるいは（2.9）が示すように各々の地震トリガリング効果は「積」，即ち乗法的に現れる．この違いは 2.2 節で触れた
両モデルの非線形／線形モデルの違いでもある．
さて，地震が短期間に続発し，そのうちの最初の 1 つ目が本震に相当するような大きな mを持つ場合を考えてみる．ETAS モデルにおいては一番最初の地震により強度関数を十分に大きくしても問題がない．一方，Dieterich モデルにおいては，一番最初の地震で強度関数を大きくしてしまうと，乗法的であるがためにその後に続く地震のトリガリング効果でさらに強度関
数が雪だるま式に上がっていき，強度関数の値が非現実に高くなり過ぎてしまうことになる．
これを避けるため，実際の解析においては，最初の地震直後の強度関数を抑制し，その後の地
震のトリガリング効果で強度関数を高くすることになる．言い換えると，Dieterich モデルは最初の地震発生から，いくらか遅れて強度関数が最大値を取ることになる．一方，ETAS モデ
ルにおいては最初の地震の直後に，強度関数の値が最大になることが可能である．現実の地震活動に見られる典型的な本震・余震系列においては，最初の地震が本震であり，それが最大のm を取ることが殆んどである．そして，その本震発生直後に最も地震活動が激しくなる．線形Hawkes 過程である ETAS モデルの加法性はこれに合致した強度関数を生み出せる一方，非線形 Hawkes モデルである Dieterich モデルの持つ乗法性は矛盾した強度関数となってしまう．このことが新潟県中越地震に対しても Dieterich モデルの適合性が ETAS モデルのそれに劣った
大きな要因である．
4.

Dieterich モデルに関する今後の展開

2.3 節に記した通り，大森・宇津公式におけるベキ指数 p = 1 の場合に対応し切れていない点を改善することは，Dieterich モデルにおける大きな課題である．これについては，Dieterich（1994）が示唆し，また 3 節で実データ解析の結果が示したようにトリガリング効果を重ね合わせることで，ある程度は解決された．但し，兵庫県南部地震の余震系列に対する解析結果が示すように十分とは言えない場合もあり，さらなる工夫が必要である．
2.3 節で簡単に触れたが，Dieterich（1994）は時間変化する応力変化を導入することで p = 1とし得ることを指摘している．具体的には，(4.1)

τ = τ0 + u ln(wt + 1)というものを考え，これにより特に p > 1 となる地震活動の減衰を作り出すことが出来ることを示した．大地震の余震域では「余効変動」と呼ばれる，本震発生の後，長期間に亘ってゆっく

218

統計数理

第 69 巻 第 2 号 2021

りとしたすべり（スロースリップ）が続く現象が観測されている（例えば Heki et al., 1997）．時
定数が異なるなど，様々な種類のスロースリップが見つかっており（Obara and Kato, 2016 のレビューが詳しい），こういった現象によって非定常的な応力変化が起きている可能性はある．
という仮定に換えて式（2.1）に代入して解
実際に上記の τ の式を，「剪断応力変化 τ̇ が一定」
き，式（2.5）とまとめると，ti−1 < t ≤ ti (i = 1, 2, . . . , n + 1) に対して



(4.2)



Δτi−1
1
γ(t) = γ(ti−1 ) exp −−[w(t − ti−1 ) + 1]−mAσAσw(m + 1)
1+
[w(t − ti−1 ) + 1] （但し m = u/Aσ ）Aσw(m + 1)となる（Dieterich, 1994 の式
（B21）
を参照）
．即ち，これを式（2.6）に換えて用いることで，この場合の強度関数を得ることが出来る．
式（4.2）の示すところは，−m > 1 即ち m < −1 の時，第 1 項が支配的となり，近似的にγ(t) ∝ (t − ti−1 )−m であることである．強度関数は式（2.7）に示した通り γ(t) の逆数であるゆmえ，これは (t − ti−1 ) で減衰することになる．こうして，大森・宇津公式の p = 1 の場合を，
u と Aσ の比を調整することで，Dieterich モデルからも作り出すことが出来る．
では，式（4.2）に基づく強度関数を用いれば，更にモデルが改善されるだろうか．Iwata（2010）（以下，モデル C）を適用しに基づき，試みとして 3 節で扱った 2 つの余震系列に，このモデル
た結果を示す．すると，兵庫県南部地震に対しては AIC の値が小さくなる一方，新潟県中越地．
震に対しては大きくなった（表 1）
合わせて強度関数がどのようになったかを見てみると，図 5 にある通り，新潟県中越地震に対しては殆んどの凸凹が消滅し，ごく限られたものによるトリガリングの影響が見られるだけである．またその影響もごくわずかである．兵庫県南部地震に至っては，（強度関数の全体的な減衰を生み出した本震によるものを除いて）全くトリガリングの影響が現れていない．
このようになった理由は，応力変化が時間一定であった場合と式（4.1）のようなものであった
場合の強度関数の特徴の違いにある．図 6 に両者の応力変化に 2 回の step 的な応力変化を加えた場合の地震活動度（強度関数）の時間変化の例を示す．図には合わせて step 的な応力変化を与えなかった場合や，2 回のうちの 1 回だけを与えた場合のものも示してある．「時間一定」の場合は step 的な応力変化の有無に関わらず最終的には同じ状態へと戻っていく．一方，式図 6．複数の step 的な応力変化があった場合（黒実線）となかった場合（灰色点線）とで強度関数の時間変化がどのように異なるかを示す概念図．（a）応力変化一定の場合および（4.1）に従う応力変化の場合．
（b）式
A schematic diagram showing the temporal evolution of the intensity function(seismicity rate) with multiple stepwise stress changes (solid black curves) andwith none of one stepwise stress change (dotted gray curves) for (a) constantstress rate and (b) stress change as shown in eq. (4.1).

地震学における非線形 Hawkes 過程：摩擦構成則に基づく地震活動モデル

219（4.1）のような応力変化の場合は，step 的な応力変化で地震活動度が跳ね上がってしまうとそ「各々の余震の影響が永劫的に残り，step 的な応力変化がなかった状態には戻らない．これは活動はいつしか定常的な活動に飲み込まれ（少なくとも見かけ上は）消滅する」というこれまでの観測事実（例えば Ogata and Shimazaki, 1984）と矛盾している．そのため，本震以外の地震によるトリガリングの効果を極力抑えようとし，殆んどあるいは全く凸凹のない強度関数が選ばれることとなった．つまりトリガリング効果が殆んど抑えられた結果，新潟県中越地震に対しては AIC の値が悪くなった．他方，元々トリガリング効果がさほど大きくなく大森・宇津公
式の p = 1.22 に対応できていなかった兵庫県南部地震については，式（4.1）の応力変化を導入（実際，強度関数は大森・宇津公式によすることで，この p の値には対応できるようになったるものと殆んど一致しており図の上では見分けが付かない）．これによるモデル改善が大きく，
AIC の値はよくなったと考えられる．
以上の結果が示すことは，応力変化が時間変化するようにして地震活動の減衰の多様性と，全ての地震のトリガリング効果の影響を Dieterich モデルで両立させることは簡単ではなく，適切な定式化が必要であるという点である．また，ここでは Dieterich（1994）の提案に基づき
式（4.1）の応力変化を扱ったが，この微分は t → ∞ で漸近的に 0 に近付くため，地震活動度も0 に近付く（図 6（b）参照）．これについては式（4.1）に換えて，例えば τ = τ0 + u ln(wt + 1) + ct
（c は定数）のようなものを考えればよいが，トリガリング効果の影響はより複雑となり，やはり定式化における何らかの工夫が必要となるであろう．
また 3 節の最後にある通り，特に本震・余震系列に Dieterich モデルを適用した場合，最も地震活動が盛んである筈の本震発生直後に強度関数が最大になることが出来ないという問題もある．これは Dieterich モデルが非線形 Hawkes 過程であるという本質的な点に起因するものであるが，ともあれ現実の地震活動とどう整合させるかが，今後の展開に向けて重要である．
この他，2.1 節に記した通り，元々の摩擦構成則で考えられていたものは剪断応力であった
のに，実際のモデリングにおいては長年に亘って CFS を用いており，モデルの精密化においてはこの点についての再考の必要があろう．さらに，式（2.3）や（2.5）の導出の際には「法線応力σ が変化しない」という仮定が置かれているが，実際のところは非現実的なものである．そもが現実の余震活動を最もよく説明する応力成分かという問題もあそも剪断応力（あるいは CFS）
る（DeVries et al., 2018）．即ち，摩擦構成則の取り扱い方や，摩擦構成則で入力とされる物理量をそのまま現実の地震活動に適用していいのかという点について，まだまだ考察する必要があると言える．
冒頭に述べた通り，現時点では ETAS モデルは現実の地震活動に最もよく適合するモデルである．とはいえ，ETAS モデルとてあらゆる地震活動に完璧に合うものということではなく，ETAS モデルにも様々な改善の余地があろう．ここまで見てきたように Dieterich モデルは様々な問題を抱えたままであるが，ETAS モデルと対照的な特徴を持つものである．こういった対照的なモデルを扱うことは，データ生成の背景構造（ここでは地震発生のメカニズム）への理解を深めるのに役立つし，今後，ETAS モデルや Dieterich モデルを超える地震活動モデルの誕生につながるであろう．
，
なお，本稿では触れなかったが ETAS モデルが時空間へと拡張されたように（Ogata, 1998）
Dieterich モデルの時空間版への拡張も行なわれつつある（Cattania et al., 2015, 2018; Mancini
et al., 2019）．但し，応力変動の与え方に空間的な差異を考慮する必要が生じるため，当然ながらより複雑なモデリングとなる．これをどう与えるかに関してはあまり踏み込んだ解析が行な
われていないため，この方面についても今後の発展が期待される．

統計数理

220

謝

第 69 巻 第 2 号 2021

辞

統計数理研究所の庄建倉氏には，本稿を書く機会を提供して下さると共に，非線形 Hawkes
過程に関する議論をして頂きました．同研究所の尾形良彦氏には，出版前の解説に関する情報を提供して頂きました．東北大学の遠田晋次氏，Temblor 社の Ross S. Stein 氏からは図の原画を提供して頂きました．匿名の査読者からは数々の適切なご指摘を頂き，それにより本稿は大きく改善されました．本稿の作成にあたっては，公益財団法人サタケ技術振興財団・大学研究助成金および日本学術振興会・科学研究費補助金（基盤研究 (B)，課題番号 21H01191）を使用しました．ここに記して謝意を表します．

参

考

文

献

Akaike, H. (1974). A new look at the statistical model identiﬁcation, IEEE Transactions on AutomaticControl, 19(6), 716–723.
Borovkov, K. and Bebbington, M. S. (2003). A stochastic two-node stress transfer model reproducing
Omori’s Law, Pure and Applied Geophysics, 168(8), 1429–1445.
Cattania, C., Hainzl, S., Wang, L., Enescu, B. and Roth, F. (2015). Aftershock triggering by postseismic stresses: A study based on Coulomb rate-and-state models, Journal of Geophysical Research:
Solid Earth, 120(4), 2388–2407.
Cattania, C., Werner, M. J., Marzocchi, W., Hainzl, S., Rhoades, D., Gerstenberger, M., Liukis, M.,
Savran, W., Christophersen, A., Helmstetter, A., Jimenez, A., Steacy, S. and Jordan, T. H.
(2018). The forecasting skill of physics-based seismicity models during the 2010-2012 Canterbury,
New Zealand, earthquake sequence, Seismological Research Letters, 89(4), 1238–1250.
Daley, D. and Vere-Jones, D. (2003). An Introduction to the Theory of Point Processes, I, 2nd ed.,Springer-Verlag, New York.
DeVries, P. M. R., Viégas, F., Wattenberg, M. and Meade, B. J. (2018). Deep learning of aftershock
patterns following large earthquakes, Nature, 560, 632–634.
Dieterich, J. (1994). A constitutive law for rate of earthquake production and its application to earthquake clustering, Journal of Geophysical Research: Solid Earth, 99(B2), 2601–2618.
Dieterich, J. H. (1979). Modeling of rock friction: 1. Experimental results and constitutive equations,
Journal of Geophysical Research: Solid Earth, 84(B5), 2161–2168.
Dieterich, J. H. (1986). A Model for the Nucleation of Earthquake Slip, Earthquake Source Mechanics(eds. S. Das, J. Boatwright and C. H. Scholz), 37–47, American Geophysical Union, Washington
D.C.
Dieterich, J. H. (1987). Nucleation and triggering of earthquake slip: Eﬀect of periodic stresses,
Tectonophysics, 144(1), 127–139.
Dieterich, J. H. (1992). Earthquake nucleation on faults with rate-and state-dependent strength,
Tectonophysics, 211(1), 115–134.
Dieterich, J. H. (2015). Applications of rate- and state-dependent friction to models of fault-slip andearthquake occurrence, Treatise on Geophysics (ed. G. Schubert), 2nd ed., 93–110, Elsevier, Oxford.
Hawkes, A. G. (1971). Spectra of some self-exciting and mutually exciting point processes, Biometrika,
58(1), 83–90.
Heki, K., Miyazaki, S. and Tsuji, H. (1997). Silent fault slip following an interplate thrust earthquakeat the Japan Trench, Nature, 386, 595–598.
Helmstetter, A. and Shaw, B. E. (2009). Afterslip and aftershocks in the rate-and-state friction law,
Journal of Geophysical Research: Solid Earth, 114(B1), B01308.
Iwata, T. (2008). Low detection capability of global earthquakes after the occurrence of large earth-地震学における非線形 Hawkes 過程：摩擦構成則に基づく地震活動モデル

221

quakes: Investigation of the Harvard CMT catalogue, Geophysical Journal International, 174(3),849–856.
岩田貴樹 (2009). 地震活動の時空間分布を表すための定量的モデルの展開, 地震 第 2 輯, 61（特集号）,
133–141.
Iwata, T. (2010). An application of rate- and state-friction model to observed aftershock sequenceswith logarithmical stress evolution in time, EOS Transactions, American Geophysical Union, 91,
Fall Meeting Supplement, Abstract S33B–2095.
Iwata, T. (2016). A variety of aftershock decays in the rate- and state-friction model due to the eﬀectof secondary aftershocks: Implications derived from an analysis of real aftershock sequences, Pureand Applied Geophysics, 173(1), 21–33.
Jeager, J., Cock, N. and Zimmerman, R. (2007). Fundamentals of Rock Mechanics, 4th ed., Blackwell,
Oxford.
King, G. C. P., Stein, R. S. and Lin, J. (1994). Static stress changes and the triggering of earthquakes,
Bulletin of the Seismological Society of America, 84(3), 935–953.
Mancini, S., Segou, M., Werner, M. J. and Cattania, C. (2019). Improving physics-based aftershock
forecasts during the 2016-2017 Central Italy Earthquake Cascade, Journal of Geophysical Research: Solid Earth, 124(8), 8626–8643.
Nanjo, K. Z., Enescu, B., Shcherbakov, R., Turcotte, D. L., Iwata, T. and Ogata, Y. (2007). Decayof aftershock activity for Japanese earthquakes, Journal of Geophysical Research: Solid Earth,112(B8), B08309.
Obara, K. and Kato, A. (2016). Connecting slow earthquakes to huge earthquakes, Science, 353(6296),
253–257.
Ogata, Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes, Journal of the American Statistical Association, 83(401), 9–27.
Ogata, Y. (1998). Space-time point-process models for earthquake occurrences, Annals of the Instituteof Statistical Mathematics, 50(2), 379–402.
Ogata, Y. and Shimazaki, K. (1984). Transition from aftershock to normal activity: The 1965 Rat
Islands earthquake aftershock sequence, Bulletin of the Seismological Society of America, 74(5),1757–1765.
Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space, Bulletin of the
Seismological Society of America, 82(2), 1018–1040.
近江崇宏, 野村俊一 (2019).『点過程の時系列解析』, 共立出版, 東京.
Toda, S. and Stein, R. (2003). Toggling of seismicity by the 1997 Kagoshima earthquake couplet: Ademonstration of time-dependent stress transfer, Journal of Geophysical Research: Solid Earth,
108(B12), 2567.
Toda, S., Stein, R. S., Reasenberg, P. A., Dieterich, J. H. and Yoshida, A. (1998). Stress transferred by
the 1995 Mw = 6.9 Kobe, Japan, shock: Eﬀect on aftershocks and future earthquake probabilities,
Journal of Geophysical Research: Solid Earth, 103(B10), 24543–24565.
Toda, S., Stein, R. S., Beroza, G. C. and Marsan, D. (2012). Aftershocks halted by static stress shadows,
Nature Geoscience, 5(6), 410–413.
Utsu, T. (1961). A statistical study on the occurrence of aftershocks, Geophysical Magazine, 30(4),521–605.
宇津徳治 (1999).『地震活動総説』, 東京大学出版会, 東京.
宇津徳治 (2001).『地震学』, 第 3 版, 共立出版, 東京.
八木勇治 (2012). 2011 年東北地方太平洋沖地震の震源過程, 地震 第 2 輯, 64(3), 143–153.
Yukutake, Y., Iwata, T. and Iio, Y. (2020). Estimation of the heterogeneity of stress ﬁelds using misﬁt
angles in focal mechanisms, Tectonophysics, 790, 228553.
庄建倉, 尾形良彦 (2021). 統計地震学における ETAS モデル — その進展とホークス型モデル, 統計数理,
69(2), 145–163.

222

Proceedings of the Institute of Statistical Mathematics Vol. 69, No. 2, 209–222 (2021)

A Nonlinear Hawkes Process in Seismology: A Seismicity Model Based on
Rate- and State-dependent Friction Law
Takaki Iwata1,2
1 Comprehensive Education Center, Prefectural University of Hiroshima
2 The Institute of Statistical Mathematics

This review introduces the Dieterich model, which is a seismicity model based on
rate- and state-dependent friction law. Herein the Dieterich model is compared to theepidemic-type aftershock sequence (ETAS) model, which is the best-known model forseismicity analysis. The Dieterich model is a nonlinear Hawkes process, while ETAS is alinear Hawkes process. This review also presents other contrastive characteristics between
the two models. Additionally, the advantages and disadvantages of the Dieterich model
along with future directions are discussed.

Key words: Seismicity, nonlinear Hawkes process, rate- and state-dependent friction law, point-process analysis, stress, Dieterich model.