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

特集「Hawkes 過程の新展開と応用」［原著論文］

余震誘発効果を考慮した繰り返し地震の予測
野村 俊一1 ・田中 昌之2（受付 2021 年 2 月 4 日；改訂 5 月 24 日；採択 5 月 25 日）

要旨

地震活動の予測に用いられる点過程モデルは，最も単純なポアソン過程を除けば，活断層な
どの同一震源上で周期的に繰り返される地震に対する更新過程と，各地震の余震誘発効果を取り入れた ETAS（Epidemic-type aftershock sequence）モデルとに大別される．しかし，数年周期で繰り返される比較的小さい繰り返し地震は，普段は周期的に発生しながらも，周辺で起こる大地震の影響を受けて発生間隔が急激に狭まることがあり，更新過程と ETAS モデルの両方の特徴を併せ持っている．本稿では，このような繰り返し地震に対して，大地震による余震誘発
効果を相対的な時間進行率の変化として取り入れた非定常更新過程モデルを提案する．提案モ
デルを東北地方太平洋沖に分布する繰り返し地震系列群へと適用し，2011 年東北地方太平洋沖
地震の影響を踏まえた将来の地震発生確率を評価する．
キーワード：繰り返し地震，確率予測，Brownian Passage Time 分布，更新過程，ETAS
モデル．
1.

はじめに

地震は岩石破壊によって引き起こされる断層の急激なずれ滑りがもたらす自然現象である．
同じ断層では，地震が起こるとそれまで断層面に蓄えられた応力がある程度解放されて，そこからプレート運動などによる応力の蓄積が再開するという一連のサイクルによって地震が繰り返されている．したがって，この 1 サイクルにかかる期間すなわち地震発生間隔を確率変数とみなして，それぞれが独立同一分布に従うと仮定された点過程モデルである更新過程がこうした繰り返し地震の解析と予測にしばしば用いられている（例えば，地震調査研究推進本部地震調査
委員会の予測手法，https://www.jishin.go.jp/main/choukihyoka/01b/chouki020326.pdf など）．
国内の主要な活断層や海溝沿いのプレート境界で繰り返される大地震に対しては，政府機関の地震調査研究推進本部（地震本部）
が将来 10〜300 年の地震発生確率を評価し一般公開している（https://www.jishin.go.jp/evaluation/long_term_evaluation/）．地震本部による将来地震発
生確率の評価では，最後の活動時期が推定されている活断層に対して更新過程が適用され，そうでない活断層にはポアソン過程が適用されている．これらの大地震は繰り返される周期がそれぞれ数十年から数万年と非常に長いため，予測が適切であったかを将来の地震活動から検証するには少なくとも数十年の期間を要する．そこで，Okada et al.（2012）は日本海溝から東北
日本太平洋沿岸にかけてのプレート沈み込み帯において同一震源で数年周期に繰り返されているより小規模な繰り返し地震に対して，岡田 他（2007）および岡田（2009）で提案された対数正1

早稲田大学 商学学術院会計研究科：〒 169–0051 東京都新宿区西早稲田 1 丁目 62

気象研究所 地震津波研究部：〒 305–0052 茨城県つくば市長峰 1 丁目 1

統計数理

240

第 69 巻 第 2 号 2021

規分布を用いた更新過程による予測実験を行い，2006 年 7 月から 2010 年の期間において 1 年単位で予測性能を評価した結果を報告している．また，田中他（2018）は同カタログにて繰り返
し回数や予測期間による予測性能の違いについて詳細に議論している．これらの予測実験では概ね良好な予測成績が得られてきたが，2011 年東北地方太平洋沖地震（M9.0）が起こって以降，当該地震の周辺に位置する繰り返し地震の繰り返し間隔が非常に短くなり予測を難しくする状
況が続いてきた（田中・岡田, 2015, 2016; 田中, 2017, 2018, 2019, 2020）．
小規模で発生間隔の短い繰り返し地震に対する予測では，応力蓄積率の短期的な変化に伴って繰り返し地震の発生間隔も変化するという課題がある．小規模な地震活動は，周辺の大地震
による余震活動や，スロースリップと呼ばれる通常の地震よりはるかにゆっくりしたプレートの広範囲でのすべりなどの影響を受けてしばしば変化している．繰り返し地震の予測に従来用いられてきた更新過程は，プレート運動による応力の蓄積が一定の速度で進んでいくことを前提とした定常なモデルであったが，上記の現象によりプレート運動が加速すると応力の蓄積率も加速されてしまうため，小規模な繰り返し地震に対する定常なモデルの当てはまりはときに非常に悪くなる．そのような繰り返し間隔の変化に対して，Nomura et al.（2014）は定常な更新過程に対して時間変化する応力蓄積率に比例した時間変換を施すことで得られる非定常更新過程を提案し，プレート境界に分布する繰り返し地震活動に適用することで相対的な蓄積率の時空間変化を推定した．その際，相対的蓄積率の時空間変化は区分多（relative loading rate）
項式である B - スプライン関数により表現しており，そのため過去の蓄積率の推移をそのまま将来へと補外して予測することは難しい．
そこで本研究では，地震の余震誘発効果を取り入れた点過程である ETAS モデル（Ogata,の考え方をこの非定常更新過程に取り入れたモデルを提案する．マーク付き Hawkes 過1988）
程の一種である ETAS モデルでは，点過程の条件付き強度が定数のバックグラウンド強度と各地震の余震誘発強度の和で表現される．提案モデルでは，上に述べた非定常更新過程の相対的な蓄積率の推移を，定数と大地震後の余震誘発効果の和で表現する．それにより将来の蓄積率の推移の予測が可能となるため，大地震後の余震誘発効果を考慮した繰り返し地震の予測手法も同時に提案する．
提案手法を東北地方太平洋沖に存在する中規模の繰り返し地震群の 2019 年までの発生デー
タへと適用し，2011 年東北地方太平洋沖地震以降の余震発生頻度の時間推移を大森・宇津の法則に従って推定および予測した上で，2020 年における繰り返し地震の発生確率の評価を与える．
2.

繰り返し地震のモデル

本節では，提案モデルとその予測手法を示す．地震調査研究推進本部（以下，「地震本部」という）が国内の主要活断層における繰り返し地震に適用している更新過程から出発し，それを非定常な更新過程へと拡張し，最後に大地震後の余震誘発効果を取り入れた提案モデルを導入
していく．
2.1 定常更新過程
ランダムなイベント発生時刻の系列 t1 < t2 < · · · < ti < · · · を考え，その発生間隔を
Xi = ti+1 − ti , i = 1, 2, . . . とおく．この発生間隔 {Xi } が独立同一分布に従うとき，{ti } は更新過程に従うという．更新過程は点過程モデルの一種であり，過去の発生履歴が与えられた下での次の条件付き強度関数 h(t|Ht ) によって特徴付けられる：(2.1)

h(t|Ht ) = lim

Δ↓0

P (イベントが [t, t + Δ) で起こる |Ht ).
Δ

余震誘発効果を考慮した繰り返し地震の予測

241

ここで，Ht は時刻 t 以前の過去の発生履歴 {ti ; ti < t} を表す．更新過程では各発生間隔が独立同一分布に従うために，時刻 t 以前の最後の発生時刻を tL とおくと，条件付き強度はh(t|Ht ) = ν(t − tL )(2.2)

の形で表すことができる．この非負値関数 v(·) はハザード関数と呼ばれる．ハザード関数は，発生間隔 X が x をちょうど超えた時点での瞬間的な発生確率として次のように求められる：(2.3)

v(x) = lim

Δ↓0

P (x ≤ X < x + Δ|X > x)
f (x)
=
.
Δ
1 − F (x)た だ し ，F (·) は 発 生 間 隔 の 累 積 分 布 関 数 ，f (·) は 発 生 間 隔 の 確 率 密 度 関 数 で あ る ．未
知 の パ ラ メ ー タ θ に 依 存 す る 条 件 付 き 強 度 hθ (t|Ht ) と ，観 測 期 間 [S, T ) 内 の 発 生 時 刻
S ≤ t1 < t2 < · · · < tn < T が与えられたとき，点過程の尤度は

  T

(2.4)

L[S,T ) (θ|t1 , . . . , tn ) = exp

−


n

hθ (t|Ht )dtShθ (ti |Hti ),i=1と定義される（Daley and Vere-Jones, 1988）．点過程が定常な更新過程であるならば，時刻 t 以
前の最後の発生時刻を tL としたときの条件付き強度は tL ≥ 0 のとき(2.5)

hθ (t|Ht ) =f (t − tL |θ),1 − F (t − tL |θ)と表される．特に，観測期間の初期時刻 S を繰り返し地震の最初の発生時刻 t1 の直後とした
とき，尤度は(2.6)

L[S,T ) (θ|t1 , . . . , tn ) = {1 − F (T − tn |θ)}n
f (ti − ti−1 |θ),

i=2

となる．
観測期間の終了時刻 T から将来の期間 [T, T + Δ) における地震発生確率は，次式により求
まる：
(2.7)

FT (Δ|θ) =F (T + Δ − tn |θ) − F (T − tn |θ).
1 − F (T − tn |θ)

このときパラメータ θ の値が未知であれば，たとえば上の尤度から得られる最尤推定値
θ̂ = arg maxθ L[S,T ) (θ|t1 , . . . , tn ) を代入した FT (Δ|θ̂) により将来の地震発生確率を評価できる．
2.2 Brownian Passage Time 分布
地震本部では国内の主要な活断層における繰り返し地震の予測にこの更新過程を適用している．更新過程における発生間隔の確率分布としては，いくつかの候補となる分布の中から
Brownian Passage Time 分布（以下「BPT 分布」という．
）と呼ばれる確率分布が採用されている．
BPT 分布は，統計学では逆ガウス分布やワルド分布の名で知られているが，地震学においては背景として仮定される物理モデル（ブラウン緩和振動過程）に因んでこの名前で呼ばれている．
Matthews et al.（2002）
により提案されたブラウン緩和振動過程は，断層で地震を起こそうとする応力状態 S(t) の時間推移を次の確率微分方程式により定義している：(2.8)

dS(t) = λdt + σdW (t).

この式の右辺は，プレート運動による定常的な応力蓄積に対応するドリフト項 λdt と，周辺地（W (t) は標準ブラウン運動）の和として定義
震活動による応力擾乱に対応する拡散項 σdW (t)

統計数理

242

第 69 巻 第 2 号 2021

されている．応力が一定の破壊強度 Sf に達したときに地震が起こり，そのときの応力開放によって応力は一定の水準 S0 まで降下するという一連の物理的プロセスを仮定している．さらに，応力状態 S(t) が断層の破壊強度 Sf に達すると地震が発生し，一定水準 S0 まで応力降下することを仮定することにより，応力降下から再び破壊強度に至るまでの地震発生間隔が理論
的に BPT 分布に従うこととなる．
BPT 分布の累積分布関数は，標準正規分布の累積分布関数 Φ を用いて



(2.9)

F (x|μ, α) = Φ

x−μ
√
α μx





2
+ exp
α2

と表され，確率密度関数は(2.10)

f (x| μ, α) =μexp
2πα2 x3



x+μ
Φ − √
α μx



(x − μ)2−2μα2 x(x > 0),(x > 0),と表される．ここで，x は発生間隔であり，μ と α はその平均および変動係数（標準偏差／平均）を表すパラメータとなっている．これらのパラメータは，上述したブラウン緩和振動過程
におけるパラメータと以下のとおり対応している．
(2.11)

μ=(2.12)

α=Sf − S 0,λσλ(Sf − S0 ).

上式から分かるように，地震活動間隔の平均 μ は応力の蓄積率 λ と反比例し，変動係数 α は拡散項の係数 σ に比例した関係となっている．
2.3 非定常更新過程
ここまで紹介した更新過程は，プレート運動による恒常的な一定の応力蓄積率を前提として発生間隔が独立同一分布に従う定常なモデルであった．ところが，実際のプレート運動は時空間的に不均質であり，大地震やスロースリップなどの影響を受けてしばしば変化する．短期的なプレート運動の加速などの影響で応力の蓄積速度が変化すると，その影響で小規模な繰り返
し地震の発生間隔も変化するため，定常な更新過程では当てはまりが悪くなってしまう．
こうした繰り返し地震活動の変化に対応するため，Nomura et al.（2014）は BPT 分布の更新
過程を時間的に非定常なモデルへと拡張して小繰り返し地震に適用している．繰り返し地震系
列の震源断層における相対的蓄積率の変化を時間関数 ν(t) で表現すると，BPT 分布更新過程
における（2.9）式の応力推移は次式の形へと変わる：(2.13)

dS(t) = λ · ν(t)dt + σ ·ν(t)dW (t).

この応力変分の式では，ドリフト係数が λ · ν(t)，拡散係数（拡散項の係数の二乗）が σ 2 · ν(t) と，それぞれ相対的蓄積率 ν(t) に比例している．これは，拡散項は周辺地震活動による独立増分を表すため，周辺地震の発生数が相対的蓄積率 ν(t) に比例して変化すると，分散の加法性から拡散係数も同様に変化するという考えに基づいている．このようなモデルでは，繰り返し地震の発生間隔はそのまま BPT 分布には従わなくなるが，相対的蓄積率 ν(t) に比例して時間が進行す
t（2.9）
るような時間変換 Λ(t) = S v(s)ds を施すことにより，変換後時間における応力推移式は式のようになり，変換後時間 Λ(t) における発生間隔が BPT 分布に従うようになる．このよ
うに，時間変換によって通常の更新過程に帰着される非定常更新過程は Trend renewal process（Lindqvist et al., 2003）と呼ばれる．相対的蓄積率 ν(t) = ν(t|θ) もパラメータ θ の一部に依存余震誘発効果を考慮した繰り返し地震の予測

243
t

するものとして，発生時刻 t1 , . . . , tn および観測終了時刻 T の Λ(t) = Λ(t|θ) = S v(s|θ)ds による変換後時間 Λ(t1 ), . . . , Λ(tn ), Λ(T ) を用いると，観測期間の初期時刻 S を繰り返し地震の最初の発生時刻 t1 の直後としたときの非定常更新過程の尤度は次の形となる（近江・野村, 2019）：(2.14)

L̃[S,T ) (θ|t1 , . . . , tn ) = (1 − F (Λ(T ) − Λ(tn )|θ))n

f (Λ(ti ) − Λ(ti−1 )|θ)

i=2

n


ν(ti |θ).

i=2

ここで，
（2.6）式の定常な更新過程の尤度と区別するために非定常更新過程の尤度を L̃ と表している．
2.4 大地震の余震誘発効果を考慮した繰り返し地震予測
大地震後の余震発生頻度の推移を記述する大森・宇津の法則（宇津, 1957）を踏まえて，本震発生時刻を tmain とおいて蓄積率の推移を次のようにモデル化する：(2.15)vOU (t|K, c, p) =1,K,
1 + (t−tmain+c)pt ≤ tmain ,t > tmain .

ここで，パラメータ θ には K, c, p が含まれることとなる．このとき(2.16) t
 t − S,vOU (s|K, c, p)ds =ΛOU (t|K, c, p) =K1
− (t−t 1 +c)p−1 ,t − S + p−1 cp−1
Smaint ≤ tmain ,t > tmain ,による時間変換を施すことにより，時間変換後の発生間隔が BPT 分布に従うようになる．
t
将来の蓄積率 ν(t) に基づく時間変換を Λ(t) = S v(s)ds とおくとき，観測期間の終了時刻 T
から将来の期間 [T, T + Δ) における地震発生確率は，次式により求まる：(2.17)F̃T (Δ|θ) =F (Λ(T + Δ) − Λ(tn )|θ) − F (Λ(T ) − Λ(tn )|θ).
1 − F (Λ(T ) − Λ(tn )|θ)

上式の分母は地震発生時刻 tn から T までの間に地震が起こらない確率，分子は地震発生時刻
tn から T + Δ までの間に地震が起こらない確率を表している．これらの確率の比を取ること
で，上式は地震発生時刻 tn から T までの間に地震が起こらなかったとき，さらに時刻 T + Δ
までの間に地震が起こらない条件付き確率となっている．このときパラメータ θ の値が未知であれば，たとえば上の尤度から得られる最尤推定値 θ̂ = arg maxθ L̃[S,T ) (θ|t1 , . . . , tn ) を代入した F̃T (Δ|θ̂) により将来の地震発生確率を評価できる．
3.

データ

本解析では，気象庁の震度計（加速度計）
で観測された震度 3 以上の地震波形記録から抽出された 2019 年末までの地震カタログから特定された東北地方太平洋沖の繰り返し地震データを用いた．繰り返し地震の特定にあたっては，3 軸（上下方向と直行する水平 2 方向）それぞれにおいて波形の相関が，2012 年以前は設置数が少ない，東北地方太平洋沖地震により一部使えない地点などがあり 0.8 以上とし，2013 年以降は 0.9 以上であることを基準とした．そうして特定された繰り返し地震系列のうち，2011 年 3 月 11 日の東北地方太平洋沖地震を挟んで地震が繰り返されている 5 系列を選んで解析対象とする．図 1 には各系列の震央位置を番号付きの星印で示している．Group 1 の震源深さは約 36 km，Groups 2-5 は 50 km 前後であり，いずれも
太平洋プレート沈み込み帯のプレート境界付近で繰り返されている地震系列となっている．
図 2 には，系列ごとの発生時期とマグニチュードを示している．同じ系列内でもマグニ

統計数理

244

第 69 巻 第 2 号 2021

図 1．繰り返し地震系列の震央位置（番号付き星印）と 2011 年東北地方太平洋沖地震の震央位置（番号無し星印）．破線は太平洋プレート境界の等深度線（Kita et al., 2010;
Nakajima and Hasegawa, 2006; 弘瀬 他, 2008; Nakajima et al., 2009）を表す．

Fig. 1. Epicenters of repeating earthquake sequences (stars with the group numbers)
and the 2011 Tohoku earthquake (a large star without a number). The dashedcontour lines represent the depth to the upper boundary of the Paciﬁc Plate.

チュードには多少のばらつきがあり，Group 2 のマグニチュードは M4.4 前後，それ以外の系列のマグニチュードは M5 前後となっている．各系列の繰り返し回数は，Group 1 が 4 回，Group2 と 5 が 7 回ずつ，Group 3 が 16 回，Group 4 が 9 回となっている．いずれの系列も，最初の地震の発生時期を各系列の観測開始時刻，2019 年末を観測終了時刻として前節で提案した非定
常更新過程を適用して解析を行い，2020 年の 1 年間および 2020〜2022 年の 3 年間における各系列の地震発生確率を評価する．
4.

解析結果

4.1 単独の系列に対する解析と予測
まず，解析に用いる 5 系列のうち，最も繰り返し回数の多い Group 3 に対して，2011 年 3 月
11 日の東北地方太平洋沖地震の発生時刻を tmain とした（2.15）式の相対的蓄積率の推移をもつ式の尤度を最大化することによりパラメータを最尤推定し非定常更新過程を適用して，（2.14）
た．ただし，更新過程の発生間隔が従う分布には，平均 μ3 ，変動係数 α の BPT 分布を適用している．以上から，推定すべきパラメータは θ = (K, c, p, μ3 , α) の 5 つであり，その最尤推定値を表 1 に示した．ここで，BPT 分布の平均パラメータである μ3 は 2011 年東北地方太平洋沖地震前の期間における平均発生間隔と解釈することができ，東北地方太平洋沖地震後には相対的蓄積率の上昇により発生間隔はそれよりも狭まることとなる．

余震誘発効果を考慮した繰り返し地震の予測245

図 2．繰り返し地震系列ごとの発生時刻とマグニチュード．灰色の縦破線は 2011 年東北地方
太平洋沖地震（M9.0）の発生時刻を表す．

Fig. 2. Magnitude-time diagrams of repeating earthquake sequences. Gray dashed vertical line represents the occurrence time of the 2011 Tohoku earthquake (M9.0).

図3（a）には，灰色線で Group 3 の累積発生回数を階段状に示している．さらに，パラメー式で与えられる変換後時間 ΛOU (t|K, c, p) を平均発生間隔 μ3タの最尤推定値を用いて，（2.16）
で割った値 ΛOU (t|K, c, p)/μ3 の推移を黒破線により示している．提案モデルでは，変換後時間
ΛOU (t|K, c, p) 上での地震発生時刻は平均 μ3 の発生間隔をもつ定常な更新過程に従うため，観測期間の開始時刻 S を繰り返し地震の最初の発生時刻 t1 としたとき，そこから j 個目の地震の発生時刻 tj+1 では ΛOU (tj+1 |K, c, p) ∼= jμ3 となるべきである．よって，モデルが上手く当てはまっていれば，黒破線の平均発生間隔で割った変換後時間 ΛOU (t|K, c, p)/μ3 と灰色線がおお（a）の灰色線と黒破線の推移は概ね合
よそぶつかる時点で次の地震が起こることになる．図 3
致しており，概ね良く当てはまっていると考えられる．
ここで，比較のために以下の 2 つの蓄積率の関数 vlinear (t|d), vexp (t|a, b) を（2.15）式の蓄積率
vOU (t|K, c, p) の代わりに当てはめることを考える：



(4.1)

vlinear (t|d) =


(4.2)

vexp (t|a, b) =1,d,t ≤ tmain ,t > tmain ,1,
1 + a exp{−b(t − tmain )},t ≤ tmain ,t > tmain .

統計数理

246

第 69 巻 第 2 号 2021

表 1．各解析における非定常更新過程の最尤推定値と AIC．

Table 1. Maximum likelihood estimates and AIC of the nonstationary renewal process
in each analysis.

前者は東北地方太平洋沖地震前後で蓄積率が d 倍に増加するものであり，後者は東北地方太平洋沖地震後の蓄積率の増分が指数関数的に減衰するものとした．上式の蓄積率にそれぞれ従う
とき



 t

(4.3)

Λlinear (t|d) =vlinear (s|d)ds =
S



 t
(4.4)

Λexp (t|d) =t − S,
t − S + d(t − tmain ),vexp (s|a, b)ds =St ≤ tmain ,t > tmain ,t − S,
t − S + a(1 − exp{−b(t − tmain )}),t ≤ tmain ,t > tmain ,により時間変換することで，それぞれ定常な更新過程に従うこととなる．これらの蓄積率を当てはめてパラメータを最尤推定して求めた変換後時間 Λlinear (t|d)/μ3 , Λexp (t|a, b)/μ3 を図 3（b），（c）にそれぞれ黒破線で示している．図 3（b）の変換後時間は灰色実線の累積発生回数と大（c）
でも東北地方太平洋沖地震の直後と直近の 5 年間に変換後時間と累きく乖離しており，図 3
積発生回数との乖離が見られる．ここで，各モデルの当てはまりの良さを次式の赤池情報量規
準（AIC）により評価する：(4.5)

AIC = −2 max L̃[S,T ) (θ|t1 , . . . , tn ) + 2 dim θ.
θただし，θ は各モデルのパラメータセットを表し，dim θ はパラメータ θ の次元すなわちパラ（a）〜（c）の左上に各モデルの AIC を示しており，図 3（a）の大森・宇
メータの個数を表す．図 3
津の法則に基づく蓄積率 vOU (t|K, c, p) のモデルが最も AIC が良い結果となった．
そこで，大森・宇津の法則に基づく蓄積率 vOU (t|K, c, p) を用いたモデルと表 1 の 2 列目に示した最尤推定値に基づいて，将来の次の繰り返し地震の発生確率を評価する．図 4 には，変換後時間 Λ(t) = ΛOU (t|K, c, p) における最後の地震から将来の次の地震までの経過時間 Λ(t) − Λ(tn )の累積分布関数 F (Λ(t)–Λ(tn )| μ̂3 , α̂)（最尤推定値 μ̂3 = 6.31, α̂ = 0.29 を代入している）を示し

余震誘発効果を考慮した繰り返し地震の予測

247

図 3．Group 3 の累積発生回数（灰色実線）と平均発生間隔で標準化された変換後時間 Λ(t)/μ3の推定値（黒破線）．変換後時間には（a）Λ(t) = ΛOU (t|K, c, p)（ 大森・宇津の法則），（指数関数）を用いた．各地（b）Λ(t) = Λlinear (t|d)（線形関数），（c）Λ(t) = Λexp (t|a, b)
震の発生直前における，黒破線の変換後時間と灰色実線の累積発生回数の値の近さによって，モデルの当てはまりの良さを判断できる．

Fig. 3. Cumulative count of the events in Group 3 (gray solid line) and estimatedtransformed time Λ(t)/μ3 (black dashed line) normalized by mean recurrence
time. The time is transformed along (a) Λ(t) = ΛOU (t|K, c, p) (Omori-Utsulaw), (b) Λ(t) = Λlinear (t|d) (linear function), (c) Λ(t) = Λexp (t|a, b) (exponential function) in each panel. Goodness-of-ﬁt of the model can be judged by thediﬀerence between the gray and black lines just before each occurrence.

ている．東北地方太平洋沖地震後の蓄積率の高まりが未だ続いているために，図 4 上横軸の実時間上では最後の地震から 3 年も経過していない 2022 年の初めには変換後時間における平均発生間隔に達していることが分かる．ここで，観測終了時刻（2019 年末）から 1 年以内（2020年末まで）に次の繰り返し地震が起こる確率を（2.17）式により評価する．図 4 の 2019 年末時
での累積分布関数はほぼ 0 であるため（2.17）式の分母はほぼ 1 となるが，点（2020 年の縦点線）（2021 年の縦点線）での累積分布関数は 2.8%
（0.028）となる（2.17）式の分子にくる 2020 年末時点

248

統計数理

第 69 巻 第 2 号 2021

図 4．Group 3 の最後の地震（2019/9/25）から将来の次の地震までの変換後時間における発生間隔 Λ(t) − Λ(tn ) に対する累積分布関数 F (Λ(t) − Λ(tn )|μ̂3 , α̂)（μ̂3 = 6.31, α̂ = 0.29は最尤推定値）．上横軸と縦点線は実時間（年）を示している．

Fig. 4. Cumulative distribution function F (Λ(t) − Λ(tn )|μ̂3 , α̂) (μ̂3 = 6.31, α̂ = 0.29
are the maximum likelihood estimates) of lapse time from the latest event on
2019/9/25 to the next event. Upper horizontal axis and vertical dotted lines
show the real time (year) of the next event.

ため，（2.17）式に基づく 1 年以内の地震発生確率は 2.8% となる．同様に，図 4 の 2022 年末時
での累積分布関数は 89.1%（0.891）となるため，（2.17）式に基づく 3 年以内
点（2023 年の縦点線）に次の繰り返し地震が起こる確率は 89.1% となる．これらの確率の評価結果（2022 年末まで）
は表 2 の 2 列目にも示している．
4.2 複数の系列に対する解析と予測
以上のように Group 3 については単独の繰り返し系列で解析ができたものの，他の系列については繰り返し回数が多くないため，単独の系列で 5 つものパラメータを推定するのは現実的でない．そこで，震源が近い位置にある複数の系列に対して，一部のパラメータを共有して推定することを考える．
繰り返し地震系列が J 系列あり，j 番目の系列の繰り返し地震数を nj ，それぞれの地震発
生時刻を tj,1 , . . . , tj,nj とおく．この J 系列に対して，（2.15）式で表される共通の相対的蓄積率
ν(t|K, c, p) の推移を仮定して非定常更新過程を適用する．さらに，変換後時間 Λ(t|K, c, p) 上の発生間隔が従う BPT 分布のパラメータについて，発生間隔の変動係数 α は系列間で共通の値をもつものと仮定し，一方で発生間隔の平均は系列間で異なるパラメータ μ1 , . . . , μJ を用意する．
このとき，J 系列を合算した尤度が次式のように得られる：(4.6)

L̃(K, c, p, μ1 , . . . , μJ , α) =JL̃[Sj ,T ) (K, c, p, μj , α|tj,1 , . . . , tj,nj )j=1

ただし，Sj = tj,1 とおいた．このとき，推定すべきパラメータは θ = (K, c, p, μ1 , . . . , μJ , α) のJ + 4 個となる．
まず，Groups 1-3 の 3 系列に対して以上のように共通のパラメータを置いて推定した結果を
表 1 の 3 列目に示す．BPT 分布の平均パラメータ（東北地方太平洋沖地震前の平均発生間隔）

余震誘発効果を考慮した繰り返し地震の予測249

表 2．各解析における次の繰り返し地震の発生確率の評価結果．

Table 2. Probabilistic forecasts of the next repeating earthquakes in each analysis.

は系列ごとの大きく異なるため別々の値として推定しているが，それ以外のパラメータについては系列間で共通のものとして推定している．
図 5 には，Groups 1-3 の系列ごとに，累積発生回数を灰色線で，推定された変換後時間を平均発生間隔で割ったものの推移を黒破線で示しており，いずれの系列も両者はおおまかに合致している．さらに，表 2 の 3 列目には系列ごとの観測終了時刻の 2019 年末から 1 年以内（2020年末まで）および 3 年以内（2022 年末まで）の次の繰り返し地震の発生確率の評価結果を示している．Group 1 は最後の地震からおよそ 3 年が経過しており，1 年以内の発生確率は 6.3% で，3 年以内の発生確率は 24.1% となった．Group 2 は最後の地震から 1 年も経過していないた
め，1 年以内の発生確率は 0.2% と低いが，3 年以内の発生確率を評価すると 20.3% となった．
Group 3 の地震発生確率の評価値は前節の Group 3 単独で解析した結果と大きくは変わらない
が，1 年以内の発生確率は若干下がって 2.4%，3 年以内の発生確率はさらに上がって 94.8% となった．
次に，Groups 4-5 の 2 系列に対して平均パラメータ以外は共通のパラメータを仮定して最尤推定した結果を表 1 の 4 列目に示す．BPT 分布の変動係数パラメータ α の推定値は 0.76 となり，表 2 に示した Groups 1-3 に対する α の推定値 0.25 よりも遥かに大きい値となった．
図 6 には，Groups 4-5 の系列ごとに，累積発生回数を灰色線で，推定された変換後時間を平均発生間隔で割ったものの推移を黒破線で示している．2011 年東北地方太平洋沖地震後の黒破線の勾配は図 5 と比べると緩やかであり，Groups 4-5 における東北地方太平洋沖地震による相対的蓄積率の変化は Groups 1-3 よりも小さかったと考えられる．表 2 の 4 列目には系列ごと（2020 年末まで）
および 3 年以内（2022 年末まで）の次の繰り返し地震の 2019 年末から 1 年以内の発生確率の評価結果を示している．いずれの系列も最後の地震からある程度の時間が経過しているため，Group 4 の 1 年以内の発生確率は 60.1%，3 年以内の発生確率は 93.8% となった．
Group 5 の 1 年以内の発生確率は 40.0%，3 年以内の発生確率は 78.9% となった．

統計数理

250

第 69 巻 第 2 号 2021

図 5．Groups 1-3 の累積発生回数（灰色実線）と平均発生間隔で標準化された変換後時間の推定値（黒破線）．各地震の発生直前における，黒破線の変換後時間と灰色実線の累積発生
回数の値の近さによって，モデルの当てはまりの良さを判断できる．

Fig. 5. Cumulative count of the events in Groups 1-3 (gray solid line) and estimatedtransformed time (black dashed line) normalized by mean recurrence time.
Goodness-of-ﬁt of the model can be judged by the diﬀerence between the grayand black lines just before each occurrence.

最後に，Groups 1-5 の全ての系列に対して平均発生間隔以外は共通のパラメータを仮定して最尤推定した結果を表 1 の 5 列目に示す．このとき AIC は 93.36 となる．一方，表 1 の 3，4 列目に示した Groups 1-3，Groups 4-5 の解析結果における AIC はそれぞれ 41.35，38.41 とあり，それらを合計すると 79.76 となる．したがって，前者の全ての系列の同時解析よりも，後者のGroups 1-3 と Groups 4-5 を別々に解析した場合の方が良い AIC が得られる結果となった．
5.

終わりに

本研究では，数年周期で繰り返される中規模の繰り返し地震について，周辺域の大地震による余震誘発効果を考慮した非定常更新過程モデルによる予測手法を提案した．定常な更新過程

余震誘発効果を考慮した繰り返し地震の予測

251

図 6．Groups 4-5 の累積発生回数（灰色実線）と平均発生間隔で標準化された変換後時間の推定値（黒破線）．各地震の発生直前における，黒破線の変換後時間と灰色実線の累積発生
回数の値の近さによって，モデルの当てはまりの良さを判断できる．

Fig. 6. Cumulative count of the events in Groups 4-5 (gray solid line) and estimatedtransformed time (black dashed line) normalized by mean recurrence time.
Goodness-of-ﬁt of the model can be judged by the diﬀerence between the grayand black lines just before each occurrence.

を時間変換することで得られる非定常更新過程をベースとし，その時間変化率（相対的蓄積率）の推移を，定数項と大地震による余震誘発効果の和として表現することにより，将来にわたる
時間変化率を予測して地震発生確率を評価した．東北地方太平洋沖のプレート境界付近に分（BPT 分布）のパラメータと，布する 5 つの繰り返し地震系列を扱い，系列ごとの発生間隔分布2011 年東北地方太平洋沖地震の影響による相対的蓄積率の変化をそれぞれ推定した．ETAS モ
デルの考え方では全ての地震の影響を取り入れるべきであるが，そのいくらかは発生間隔のばらつきを表すパラメータ α に吸収されるものとして簡便的なモデルを立てた．しかし，より精が繰り返し地震の発生間隔からプレート間準静的
密な予測を行うには，Nomura et al. （2017）
滑りの時空間変化を推定したように，相対的蓄積率も東北地方太平洋沖地震のような巨大地震
だけではなく近隣の大地震やスロースリップの影響を考慮すべきである．スロースリップによる相対的蓄積率の変化を将来的に予測することはより難しい課題となるため，今後の取り組み
としたい．
謝

辞

本研究の一部は，JSPS 科研費 JP20K11722 の助成を受けたものである．また，本論文の執筆にあたり，有益な助言を下さった査読者の方々に心より感謝する．

252

統計数理

参

第 69 巻 第 2 号 2021

考

文

献

Daley, D. and Vere-Jones, D. (1988). An Introduction to the Theory of Point Processes, Springer-Verlag,
New York.

弘瀬冬樹, 中島淳一, 長谷川昭 (2008). Double-Diﬀerence Tomography 法による関東地方の 3 次元地震
波速度構造およびフィリピン海プレートの形状の推定, 地震第 2 輯, 60, 123–138.
Kita, S., Okada, T., Hasegawa, A., Nakajima, J. and Matsuzawa, T. (2010). Anomalous deepeningof a seismic belt in the upper-plane of the double seismic zone in the Paciﬁc slab beneath the
Hokkaido corner: Possible evidence for thermal shielding caused by subducted forearc crustmaterials, Earth and Planetary Science Letters, 290, 415–426.
Lindqvist, B. H., Elvebakk, G. and Heggland, K. (2003). The trend-renewal process for statisticalanalysis of repairable systems, Technometrics, 45, 31–44.
Matthews, M. V., Ellsworth, W. L. and Reasenberg, P. A. (2002). A Brownian model for recurrent
earthquakes, Bulletin of the Seismological Society of America, 92, 2233–2250.
Nakajima, J. and Hasegawa, A. (2006). Anomalous low-velocity zone and linear alignment of seismicity along it in the subducted Paciﬁc slab beneath Kanto, Japan: Reactivation of subducted
fracture zone?, Geophysical Research Letters, 33, L16309, doi: 10.1029/2006GL026773.
Nakajima, J., Hirose, F. and Hasegawa, A. (2009). Seismotectonics beneath the Tokyo metropolitan area, Japan: Eﬀect of slab-slab contact and overlap on seismicity, Journal of GeophysicalResearch, 114, B08309, doi:10.1029/2008JB006101.
Nomura, S., Ogata. Y. and Nadeau, R. M. (2014). Space-time model for repeating earthquakes andanalysis of recurrence intervals on the San Andreas Fault near Parkﬁeld, California, Journal ofGeophysical Research, 119, 7092–7122.
Nomura, S., Ogata, Y., Uchida, N. and Matsu’ura, M. (2017). Spatiotemporal variations of interplate slip rates in northeast Japan inverted from recurrence intervals of repeating earthquakes,
Geophysical Journal International, 208, 468–481.
Ogata. Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes, Journal of American Statistical Association, 83, 9–27.
近江崇宏, 野村俊一 (2019).『点過程の時系列解析』, 共立出版, 東京.
岡田正実 (2009). 繰り返し地震および余震の確率予測, 地震第 2 輯, 61, S143-S153.
岡田正実, 高山博之, 弘瀬冬樹, 内田直希 (2007). 地震長期確率予測に使用する更新過程対数正規分布モ
デルのパラメータ事前分布, 地震第 2 輯, 60, 85–100.
Okada, M., Uchida, N. and Aoki, S. (2012). Statistical forecasts and tests for small interplate repeating
earthquakes along the Japan Trench, Earth, Planets, and Space, 64, 703–715.
田中昌之 (2017). 中規模繰り返し相似地震の発生状況と発生確率（2017）Probability of moderate repeating earthquakes within one and three years from 1 February 2017, 地震予知連絡会会報,
98, 437–442.
田中昌之 (2018). 中規模繰り返し相似地震の発生状況と発生確率（2018）Probability of moderate repeating earthquakes within one and three years from 1 February 2018, 地震予知連絡会会報, 100,356–359.
田中昌之 (2019). 中規模繰り返し相似地震の発生状況と発生確率（2019）Probability of moderate repeating earthquakes within one and three years from 1 February 2019, 地震予知連絡会会報, 102,
380–384.
田中昌之 (2020). 中規模繰り返し相似地震の発生状況と発生確率（2020）Probability of moderate repeating earthquakes within one and three years from 1 January 2020, 地震予知連絡会会報, 104,
448–452.
田中昌之, 岡田正実 (2015). 中規模繰り返し相似地震の発生状況と発生確率（2015）Probability of moderate repeating earthquakes within one and three years from 1 April 2015, 地震予知連絡会会

余震誘発効果を考慮した繰り返し地震の予測

253

報, 94, 372–376.
田中昌之, 岡田正実 (2016). 中規模繰り返し相似地震の発生状況と発生確率（2016）Probability of moderate repeating earthquakes within one and three years from 1 February 2016, 地震予知連絡会
会報, 96, 456–461.
田中昌之, 岡田正実, 内田直希 (2018). 繰り返し回数の少ない繰り返し地震系列に対する長期的地震発生
確率予測の成績と検証, 地震第 2 輯, 70, 195–213.
宇津徳治 (1957). 地震のマグニチュードと余震の起こりかた, 地震第 2 輯, 10, 35–45.

254

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

Forecasting Repeating Earthquakes Considering Aftershock-triggering
EﬀectsShunichi Nomura1 and Masayuki Tanaka2
1 Faculty of Commerce, Graduate School of Accountancy, Waseda University

2 Department of Seismology and Tsunami Research, Meteorological Research Institute

Point process models to predict earthquake occurrences (other than the simplest
Poisson processes) can be roughly classiﬁed into two types: the renewal process for earthquakes repeating periodically on the same hypocenter, such as active faults; and the ETAS
(epidemic-type aftershock sequence) model, which takes the aftershock-triggering eﬀect
of every earthquake into account. However, relatively small repeating earthquakes havethe characteristics of both models: they usually recur periodically, but have very shortrecurrence intervals after large earthquakes nearby. In this paper, we propose a nonstationary renewal process model for such repeating earthquakes that incorporates theaftershock-triggering eﬀect of nearby large earthquakes as a relative change in the loadingrate. We apply the proposed model to the repeating earthquake catalog on the Paciﬁc
Plate subduction zone in the northeastern Japan and evaluate probabilistic forecasts of
the next repeating events, considering the aftershock-triggering eﬀect of the 2011 Tohoku
earthquake.

Key words: Recurrent earthquakes, long-term forecast, Brownian Passage Time distribution, renewal processes, ETAS model.