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

特集「Hawkes 過程の新展開と応用」［研究詳解］

時空間 ETAS モデルの拡張バージョンとその応用
郭 一村1 ・庄 建倉2,3（受付 2021 年 2 月 15 日；改訂 9 月 14 日；採択 9 月 16 日）

要旨

この論文は，統計地震学で広く使用されている時空間 Epidemic-Type Aftershock Sequence（ETAS）モデルのいくつかの拡張バージョンをまとめたものである．ETAS モデルの拡張バー
ジョンには，2 次元有限震源モデル，3 次元点震源モデル，および 3 次元有限震源モデルが含まれる．有限震源モデルは大地震の断層形状を，3 次元震源モデルは地震の深さを考慮に入れている．この論文は，モデル推定，確率的デクラスタリング，および地震シミュレーションに関連するアルゴリズムを示し，これまでに日本，イタリア，南カリフォルニアで実施されてきた各 ETAS モデルの適用結果をまとめている．フィッティングの結果は，2 次元有限震源および 3 次元有限震源モデルが，点震源モデルよりも大きな α 値をもたらすことにより，本震の余
震生成密度を増大させることを示している．余震の生産性と地震性すべりを比較すると，大きなすべり領域で余震が少なく，断層面上の余震生成率密度の空間パターンは，明らかにすべりの補償を示していることがわかる．
キーワード：確率密度分布，残差解析，地震誘発，余震，確率的再構成．
1.

前書き

過去数十年で，Epidemic-Type Aftershock Sequences（ETAS）モデルは，地震誘発の定量化（Ogata, 1988, 1998; Console and Murru, 2001），常時地震活動の推定（Zhuang et al., 2002;，地震活動異常の抽出（Reverso et al., 2015; Nishikawa
Console et al., 2003; Zhuang et al., 2005），および短期余震の予測（Zhuang, 2011; Tsuruoka et al., 2012）などに有効なツーand Ide, 2017）（時
ルであることが証明された．時空間 ETAS モデルでは，余震の分布は経験的なべき乗法則に従うと想間に関する大森・宇津法則（Utsu, 1969; Omori, 1894）と空間に関する逆べき乗則）
定されている．さらに，このモデルは地震のマグニチュード・頻度分布を 決定するためのグー
テンベルグ・リヒターの法則（GR 則）（Gutenberg and Richter, 1956）を ETAS モデルと組み合わせたシミュレーションを行っている．ETAS モデルは数学的にはホークス過程（Hawkes,
1971a, 1971b）の地震系列への具体化の一例である．ETAS モデルは既に数十年の歴史があり，他分野でのホークス過程の統計的推論の理論や手法の発展にどのような影響を与えた（庄・尾形, 2021）．
ETAS モデルでは，地震活動は 2 つの要素に分けられている：背景の地震活動と誘発された地震活動である．起源が不明な，第 0 世代地震とも呼ばれる背景地震活動の地震は，時間的に
1

中国科学院大学 地球動力学実験室：100049 中国北京市石景山区玉泉路 19 甲2

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

3

総合研究大学院大学 複合科学研究科統計科学専攻：〒 190–8562 東京都立川市緑町 10–3

224

統計数理

第 69 巻 第 2 号 2021

定常的であるが空間的に不均一なポアソン過程によるものである．各イベントが発生すると，それは背景地震または誘発されたイベントに関係なく，確率ルールに従っていくつかの子イベ
ントを誘発する．ETAS モデルは，直接の子孫とその親の間のタイムラグが大森・宇津法則に従って分布し，各イベントの大きさがその発生時間と場所，および親イベントの大きさに依存しないことを前提としている．背景地震（第 0 世代）によって直接誘発されたイベントを第 1 世
代，第 1 世代の直接の子孫を第 2 世代などと順番に命名する．イベントによって間接的に誘発されたイベントは，2 次子孫と呼ばれる．
時空間 ETAS モデルの条件付き強度関数は，次の式で記述できる．
(1.1)

λ(t, x, y | Ht ) = μ(x, y) +κ(mi )g(t − ti )f (x − xi , y − yi ; mi ).

i:ti <t

ここで，Ht はフィッティング期間以前の期間に記録されているイベントであり，解析地域外のイベントも含む．上式において μ(x, y) は位置 (x, y) で推定された常時地震発生率であり，κ(m)はマグニチュード m の余震生産性を示す．
κ(m) = Aeα(m−mc ) .

ここで，Mc はカットオフマグニチュードである．時刻 t における確率密度関数（Probability
Density Function, PDF）g(t) は，余震が次のように減衰する経験的な大森・宇津の法則に従うと仮定する．
p−1
t −pg(t) =
1+, t > 0.
cc
空間位置 PDF は

−q
q−1
x2 + y 2f (x, y; m) =1+
πDeγ(m−mc )Deγ(m−mc )である．
ETAS モデルが連続型ホークス過程の特殊な形式であることは明らかである．この別の説明として，位置 (x, y) と時間 t での地震リスクが，背景の地震活動と前の各イベントによって引き起こされた地震活動からなるということが挙げられる．イベントが大きく，空間と時間が近いほど，トリガー効果は大きくなっている．
時空間 ETAS モデルによる大地震のモデル化を改善することは，常に重要な目標となっている．これまでの研究によると，断層の幾何形状を無視して等方的な余震応答を仮定すると，パラメータ推定に大きな偏りが生じる可能性があり（Hainzl et al., 2008），余震分布の異方性に．
向けた改良が行われてきた（Ogata, 1998; Marsan and Lengliné, 2010; Bach and Hainzl, 2012）
この問題を解決するには, 断層形状が重要であり，このことは余震のほとんどが本震の破壊断層に沿って発生することからわかる．
大規模地震の直接余震の偏った推定を改善するために，Guo et al.（2015a）は，本震の断層形
状を組み込んだ 2 次元有限震源 ETAS モデルを提案した．このモデルは，本震の生産性を向上．
させ，本震破断部における余震の異方性を再現することに成功している（Zhuang et al., 2018）
は 3 次元震源 ETAS モデルを提案しており，このモデルでは，ベーさらに，Guo et al.（2015b）
タ分布に従って震源の深さが組み込まれている．さらに，Guo et al.（2019）は，有限震源 ETAS
モデルの 3 次元バージョンを提案し，熊本地震について検証した．それに応じて，本来の時空
間 ETAS モデルは 2 次元点震源 ETAS モデルと呼ばれる．
この論文は，時空間 ETAS モデルの拡張バージョンをまとめ，これまでのさまざまな地域への適用事例を示すことにより，余震の生成における大地震の断層形状の重要性を説明すること

時空間 ETAS モデルの拡張バージョンとその応用225を目的としている．
2.

方法論

2.1 モデルの定式化
3 次元点震源と 2 次元有限震源 ETAS モデルの条件付き強度関数は，それぞれ(2.1)

λ(t, x, y, z | Ht ) = μ(x, y, z) +κ(mi )g(t − ti )f (x − xi , y − yi ; mi )h(z; zi )i:ti <tと(2.2)λ(t, x, y | Ht ) = μ(x, y) +κ(mi )g(t − ti )f2D (x − xi , y − yi ; Si , mi )i:ti <tで記述する．ここで，z は深さを示し，Ht と μ と κ と g および f は式（1.1）と同様に定義される．3 次元点震源モデルにおける誘発された地震の震源深さがベータ分布に従うと仮定し，(2.3)

h(z, zi ) = z η zZi η(1− zZi )1 − Zz
Z
 zi



zi

ZB η Z + 1, η 1 − Z + 1

が示す PDF がある．式（2.3）
で最大の深さ Z を与えることにより，地震が限られた地震発生層に分布しているという仮定を簡単に満たすことができる．2 次元有限震源 ETAS モデルにおける誘発された地震の空間位置 PDF は

⎧

−q
q−1
x2 + y 2⎪⎪, 地震 i は点震源である場合,⎪⎨ πDeγ(mi −mc ) 1 + Deγ(mi −mc )f2D (x, y; Si , mi ) =

⎪
⎪
⎪
⎩

Si

f ∗ (x − u, y − v)τi (u, v) dudvSiτi (u, v) dudv

q−1
f (x, y) =
πD
∗



, 地震 i は有限震源である場合,

x2 + y 2
1+
D

−q

である．ここで，τi (u, v) はイベント i の地震断層の位置 (u, v) について子地震の生成率密度で
ある．
（2.2）
を組み込んで，Guo et al.（2019）は 3 次元有限震源 ETAS モデルを提案した．
式（2.1）とこのモデルの条件付き強度関数は次の式で記述する．
(2.4)

λ(t, x, y, z | Ht ) = μ(x, y, z) +κ(mi )g(t − ti )f3D (x − xi , y − yi , z; zi , Si , mi ).

i:ti <t

ここで，地震 i は点震源
（Si = {(xi , yi , zi )}）の場合，f3D (x, y, z; Si , mi ) =q−1
πDeγ(mi −mc )1+x2 + y 2
Deγ(mi −mc )−qh(z, zi ).

断層が数十 km 以上の大地震では，断層形状をモデルに組み込む必要がある．地震 i は有限震
源の場合，
(2.5)

f3D (x, y, z; zi , Si , mi ) =Sif ∗ (x − u, y − v)h(z, w)τi (u, v, w) dudvdw
Si

τi (u, v, w) dudvdw.

ここで，τi (u, v, w) は震源 Si 上の任意の 3 次元位置 (u, v, w) の子地震生成率密度である．

統計数理

226

第 69 巻 第 2 号 2021

したがって，2 次元点震源，2 次元有限震源，3 次元点震源，3 次元有限震源 ETAS モデルで推定（θ = (A, α, c, p, D, q, γ)），8 個（θ = (A, α, c, p, D, D , q, γ)），
されるパラメータはそれぞれ合計 7 個
8個（θ = (A, α, c, p, D, q, γ, η)），9 個（θ = (A, α, c, p, D, D , q, γ, η)）となる．
2.2 地震誘発の確率論的な説明
ETAS モデルでは，カタログから常時地震を抽出するために，以下のような定義を用いて，本震と余震の誘発関係を確率的に記述している．3 次元有限震源 ETAS モデルを例として，イ
ベント j が常時地震である確率 ϕj と，イベント j がそれ以前のイベント i によって引き起こされる確率 ρij を以下のように定義する．
(2.6)

ϕj =

μ(xj , yj , zj ).
λ(tj , xj , yj | Htj )

(2.7)

ρij =

κ(mi )g(tj − ti )f3D (xj − xi , yj − yi , zj ; Si , mi ).
λ(tj , xj , yj | Htj )有限震源 Si が多くのパッチ，Siι , ι = 1, 2, . . . , ni に分割されていると仮定する．同様に，Si上のパッチ ι によって引き起こされるイベント j の確率，つまり ρiιj は，次のように定義することができる．
(2.8)

Siι

ρiιj =f ∗ (x − u, y − v)h(z, w)τi (u, v, w) dudvdw

λ(tj , xj , yj Htj )Siτi (u, v, w) dudvdw

.

j−1

なお，イベント j が誘発されたイベント ρj = i=1 ρij となる確率は，ϕj + ρj = 1 である．
j = 1 の場合には，ρj = 0, ϕj = 1 である．
（2.7）
は，直接余震と二次余震を区別するための基準を与える．本震 i に対する，
式（2.6）と式
イベント j, j = i + 1, . . . , N, が直接余震である確率 ρij とすることで，その派生の過程を表すこによって提案された確率的デクラスタリング手法は，このよとができる．Zhuang et al.（2002）
うな考え方に基づいており，カタログ内のイベントをそれぞれ樹形図に分類することができる．
2.3 μ(x, y, z) と τ (u, v, w) の推定
3 次元 ETAS モデルでは，総地震発生率 M (x, y, z) をカーネル法で推定している．
(2.9)

 z dz zZi dz (1− zZi )N
1 − Zz
1 
Z



.
M̂ (x, y, z) =Kdj (x − xj , y − yj )T
ZB dz zZi + 1, dz 1 − zZi + 1
j=1

ここで，dz は深さ成分のベータカーネルのバンド幅であり，dj は経度・緯度成分のガウシア
ンカーネルのバンド幅である．



Kdj (x, y) =1x2 + y 2
exp −
2
2πdj
2d2j



.

なお，dj は，イベント j とその n 番目に近いイベントとの距離で表される異なるイベントに対して，最小のカットオフ値（一般的に 0.02〜0.05 度）を用いてそれぞれ異なる値を取るが，これは，震源決定誤差に相当すると考えられる．式（2.6）と式（2.9）を組み合わせることで，以下のようにして常時地震発生率を推定することができる．
(2.10)

 z dz zZi dz (1− zZi )N
1 − Zz
1 
Z



.
μ̂(x, y, z) =ϕj Kdj (x − xj , y − yj )T
ZB dz zZi + 1, dz 1 − zZi + 1
j=1

時空間 ETAS モデルの拡張バージョンとその応用

227

τ (u, v, w) を計算するには，式（2.5）
を離散化したものを使用する．
(2.11)

f3D (x, y, z; Si , mi ) =

ni

τiι
ι=1

Pi

f ∗ (x − uiι , y − viι )h(z; wiι ).

ni
なお，τiι は (uiι , viι , wiι ) を中心とする Si 上のパッチ ι の生産率であり，Pi = ι=1τiι はイベント i の総生産性，ni はパッチ数である．式（2.8）を用いれば，τiι = j:j>i ρiιj により τiι を容易に算出することができる．しかし，我々の経験によれば，このような単純な処理は余震が高度にクラスター化されているパッチに大きな重みを与え，最終的にはパラメーター推定の数値を不安定にする可能性がある．これを避けるために，τiι の計算に平滑化係数を導入する．

(2.12)

τ̂iι =

1
Hd1 ,d2 (ι)

niKd∗1 ,d2 (uj − uι , vj − vι , wj |wι )j=1

ρiιk .

k:k>i

ここで，平滑カーネル zd は次の形をとる．
(2.13)

Kd∗1 ,d2 (u, v, w|w ) = Kd1 (u, v)

 w 
 w d2 wZ 
w d2 1− Z
1
−
Z
Z
 w


．
w

ZB d2 Z + 1, d2 1 − Z+1

ここで，d1 および d2 は，それぞれ水平成分のガウスカーネルおよび深さ成分のベータカーネ
では，j は Si 上のすべてのパッチに適用され，次の式で正規ルの帯域幅の値である．式（2.12）
化されていることに注意する．
Hd1 ,d2 (ι) =niKd∗1 ,d2 (uj − uι , vj − vι , wj |wι ).

j=1

3.

アルゴリズム3.1 モデルの最尤推定
常時地震活動 μ(x, y, z) と本震の生産性 τ (u, v, w) が与えられると，最尤法を用いてモデルパ
ラメータを推定することができる．対数尤度関数は次のように書くことができる（Daley and．
Vere-Jones, 2003）
log L(θ) =jlog λ(tj , xi , yj , zj ) − Z   T2λ(t, x, y, z)dt dxdy dz.
0

R

T1

なお，ここで θ = (A, α, c, p, D, q, γ, η, D ) において 3 次元有限震源 ETAS モデルの場合，j は，
時間窓 [T1 , T2 ]，解析領域 R，および深度範囲 [0, Z] 内のカタログ内のすべてのイベントを対象
は，2 次元有限震源 ETAS モデルでモデルパラメータ，常時地としている．Guo et al.（2017）
震活動 μ(x, y)，および本震の余震生成率密度 τ (u, v) を同時に推定するための反復アルゴリズ
ムを提供している．3 次元有限震源の場合，以下のような同様のアルゴリズムを使用してモデル推定を行う．
S1. イベント j = 1, 2, . . . , N ごとに異なる変数ガウス帯域幅 dj と全ての j に共通な深さ成分のためのベータバンド幅 dz を設定する．
S2. = 0，μ() (x, y, z) = τ () (u, v, w) = 1，初期パラメータ θ = (A, α, c, p, D, q, γ, η, D  )，および式
（2.13）の平滑カーネルのバンド幅 d1 , d2 を設定する．
S3. 式（2.11）
および（2.12）
によって τ () (u, v, w) を数回更新して，安定した τ () (u, v, w) を取得
する．

228

統計数理

第 69 巻 第 2 号 2021

図 1．確率的除群法の図．各セグメントの長さは，対応する確率の値を表す．i つ目のイベン
トには，ni パッチを含む有限震源がある．U は，[0, 1] に一様分布する乱数である．この図では，U はセグメント ρkj に分類され，イベント k が j の親として指定されていることを示している（Zhuang et al., 2019）．

S4. 最尤推定（Maximum Likelihood Estimation, MLE）手順を使用して λ(t, x, y, z | Ht ) を(+1)として記録する．
フィットさせ，モデルパラメーター θ̂ を更新し，θ̂S5. 式（2.6）の ϕj と式（2.8）の ρiιj を，j = 1, 2, . . . , N, ι = 1, 2, . . . , ni および i = 1, 2, . . . , Nsub
について計算する．ここで，Nsub は有限震源の総数を表す．
S6. 式（2.10）により μ̂(x, y, z) を再計算し，μ(+1) (x, y, z) として記録し，式（2.12）により,(+1)(u, v, w) として記録する．
i = 1, 2, . . . , Nsub ，に対する τ̂i (u, v, w) を再計算し，τ̂i
S7. もし max |μ(+1) (x, y, z) − μ() (x, y, z)| > であると，その が収束を保証するのに十分に小さい正の数で， を + 1 と更新し，ステップ S3 に進む．それ以外の場合は，μ̂(+1) (x, y, z)と τ̂i(+1) (u, v, w) を出力し，停止する．
3.2 確率的デクラスタリング
2.2 節の式
（2.6）
と式
（2.7）
で示されているように，地震が常時イベントおよびトリガイベントになる確率は，ETAS モデルによって計算できる．トリガイベントは，それぞれのサブプロセ
スに分類され，それぞれのサブプロセスは特定のマザーイベントによってトリガされる．
これらの確率に基づいて，Zhuang et al.（2002）は地震カタログから常時地震活動を客観的に与えるための確率的デクラスタリング法を提案している．この節では，ETAS モデルの有限震源バージョンのデクラスタリングアルゴリズムをまとめる．カタログから系図を確率的に与えるアルゴリズムは以下の通りである．
S1. j = 1 を設定する．
S2. 区間 [0, 1] の一様乱数 Uj を生成する．
S3. Uj < ϕj ならば，イベント j を常時イベントとみなす．
I
S4. それ以外の場合は，Uj < ϕj + i=1 ρij を満たす最小の I を選択し，イベント j をイベ
ント I によってトリガされたものと見なす．
I−1
JS5. イベント I の震源が有限である場合，Uj < ϕj + i=1 ρij + ι=1 ρiιj を満たす最小の Jを見つける．そして，イベント j は震源 Si 上のパッチ ι によってトリガされたものとみなす．
それ以外の場合，ステップ 6 に進む．
S6. j = N の場合，停止する．そうでなければ，j = j + 1 に設定し，手順 S2 に進む．

このような，手順に固有の不確実性は，選択に使用される乱数のシードによって説明され，選択を何度も繰り返すことによって定量化される．
この手順は図 1 に説明された．
3.3 シミュレーションアルゴリズム[0, t] の観測カタログがあれば，Zhuang（2011）の 2 次元 ETAS モデルのアルゴリズムに基づいて，[t, t + Δt] の地震活動をシミュレーションすることができる．マグニチュード度数分布（GR則）で b 値を 1.0 と仮定することにより，我々のアルゴリズムを 3 次元の場合に一般化する．

時空間 ETAS モデルの拡張バージョンとその応用

229

S1. 常時地震カタログを生成する．カタログ内の各イベント i について，[0, 1] に一様分布するランダム変数 Ui を生成し，Ui < νϕi Δt/(t − t0 ) の場合受け入れる．ここで，ν はデータのフィッティングの収束を速めるために μ(x, y, z) = νμ0 (x, y, z) として導入され，t0 はフィッティングの開始時間である．選択されたイベントについて，その発生時刻を [t, t + Δt] の一様分布，震源を帯域幅 di の 2 次元ガウス分布，深さをデータフィッティングから得た η のベー(0)
タ分布でそれぞれ与える．これらの新しいイベントを Generate 0，すなわち G として記録
する．
(0)
S2. G(0) を G の集合とし，カタログ中のすべてのイベントを t 以前のものとする．
S3. ← 0 を設定する．
S4. カタログ G() の各イベント (ti , xi , yi , zi , mi ) について，
()

Oi

(i)

(i)

(i)

(i)

(i)

← {(tk , xk , yk , zk , mk ) : k = 1, 2, . . . , N (i) }として記録された N (i) 派生地震のシミュレーションを行う．N (i) は，平均値が κ(mi ) であ
(i)
(i)
(i)
るポアソン分布に従う乱数である．t(i)k , xk , yk および zk はそれぞれ，確率密度 g(· − ti )，(i)()f (· − xi , · − yi ) および h(· | zi ) から生成される．tk ∈ [t, t + Δt] で Oi のイベントを選択し，()O i として記録する．
()S5. G(+1) ← i∈G() O i を設定する．
(+1)
S6. G
が空でない場合は ← + 1 と設定して，ステップ S4 に戻り，そうでなければ(0)G ∪ j=1 G(j) を返す．
4.

データ

このセクションと次のセクションでは，新しい形式の ETAS モデルを日本，イタリア，南カ
リフォルニアのデータに適用した例を示す（Guo et al., 2017, 2021; Zhuang et al., 2018）．
2 次元点震源と 2 次元有限震源の ETAS モデルを日本地域に適用し，気象庁のカタログ（1964
年 1 月 1 日から 2014 年 7 月 31 日まで）を用いて，北緯 32〜46 度，東経 130〜148 度を調査地域とした．1983 年以降に発生した M7.5 よりも大きな 6 つの主要な地震を，有限震源とし，マ（表 1）．データとモデル設定の詳細については，Guo et
グニチュードの閾値は 4.0 に設定した
al.（2017）
を参照されたい．
イタリア地域については，INGV が編集した 2005 年から 2017 年までの ISIDE 地震カタログ
表 1．日本，イタリア，南カリフォルニアの主要な地震の一覧．SC は南カリフォルニアの
略号．

統計数理

230

第 69 巻 第 2 号 2021が使用されている．調査領域は北緯 35〜48 度，東経 6〜19 度とした．マグニチュード閾値は2.9 に設定されている．有限震源として調査地域における 6 つの M5.5 以上の本震断層を地表に投影したものを用いた（表 1）．データと研究地域の詳細は Zhuang et al.（2018）に記載されている．
南カリフォルニアへの適用には，1980 年から 2019 年までの南カリフォルニア地震データセンターのカタログデータを用いた．調査地域は，北緯 33〜37 度，西経 122〜114 度である．
1990 年以降に発生した M6.5 以上のマグニチュードを持つ 4 つの主要な地震の破断形状が，有限震源 ETAS モデルに組み込まれている（表 1）．この研究のより詳細な設定については，Guoet al.（2021）
を参照されたい．
5.

結果

5.1 モデルフィッティング（1）点震源モデルと有限震源モデルの比較．
表 2 の 1 行目と 2 行目は，JMA データへの 2 次元点震源および 2 次元有限震源の ETAS モ
デルのフィッティング結果を示している．イタリア地域へのフィッティング結果は，3 行目と5 行目に示されている．表 2 の 6〜7 行目には南カリフォルニア地域の 3 次元点震源と 3 次元
有限震源の ETAS モデルのフィッティング結果を示す．ここで，最後の列の尤度を比較すると，有限震源モデルが対応の点震源モデルより優れていることがわかる．また，余震生産性パ
ラメータ α は，点震源モデルにおいて有限震源モデルよりも小さい．これは α 値が高いほど，
小さなイベントと大きなイベントの余震生産性の差が大きくなるためである．表 2 の最後の 2
行については，マグニチュード 7.0 の場合，3 次元有限震源 ETAS モデルでは想定される余震の生産性が 34.3 から 224.5 に増加しており（Mc = 3.0），破断形状を取り入れることで大規模地震の生産性が飛躍的に向上することを示唆している．ここで α 値が大きいことは，余震がカタ
ログ内の最大のイベントらによって直接発生する可能性が高いことを意味する．α 値の違いは等方性空間応答カーネルを使用して，非等方性地震破壊源に対する余震発生率の異方性応答を，Guo et al.（2015b）おモデル化することによって引き起こされる．これも Hainzl et al.（2008）
によって結論付けられた．また，A と α は反相関関係にあり，α がよび Zhuang et al.（2019）
増加するとトレードオフとして各モデルのフィッティングで A が減少する．
2 次元モデルと 3 次元モデルの比較．
（2）2 次元 ETAS モデルと 3 次元 ETAS モデルの唯一の比較は，3 行目と 5 行目である．尤度を比較すると，3 次元モデルは 2 次元モデルよりもイタリアのデータにはるかによく適合し，尤
表 2．日本，イタリア，南カリフォルニアにおける ETAS モデルの各バージョンのフィッ（PDF: h(z) = 1/Z, 0 ≤ z ≤ Z ）テ ィ ン グ 結 果 ．*イ ベ ン ト の 深 さ が 一 様 分 布（PDF:に従う 2 次元モデルの対数尤度値．**イベントの深さがヒストグラム分布
I(int(z)<z<int(z)
+
1)
）
に従う
2D
モデルの対数尤度値．
h(z) ∝
ii

時空間 ETAS モデルの拡張バージョンとその応用

231

度の増分が 2 次元点震源モデルと 2 次有限点震源モデルの尤度差よりもはるかに大きいことがわかる．これは，深さの相関関係が地震活動の誘発に重要な役割を果たしていることを意味する．しかし，驚くべきことに，余震の生産性パラメータ α は，3 次元点震源モデルにおいて，2次元点震源モデルよりも小さくなっている．同じ結果が Guo et al.（2015a）によって発見されは，これは，余震イベントの場所が地表の 2 次元経度緯度空間に投影た．Zhuang et al.（2019）
されると，3 次元経度 - 緯度 - 深度空間のこのような異方性が減少するためであると説明した．
このように α 値が低いため，Zhuang et al.（2019）は余震を分類する可能性が高い 3 次元点震源は第 2 世代以上の本震であると断定した．当時，3 次元有限震源 ETAS モデルは実装されていは 3 次元有限震源 ETAS モデルが地震活動のモデリングにおなかったが，Zhuang et al.（2019）
いて理想的であることを示した．上記の結果は，3 次元有限震源 ETAS モデルがこれらすべてのモデルの中で最も適切であるという結論を裏付けている．
マグニチュード 7.0 の場合，3 次元有限震源 ETAS モデルでは，想定される余震の生産性が34.3 から 224.5 に増加しており（mc = 3.0），破断形状を取り入れることで大規模地震の生産性が飛躍的に向上することを示唆している．
5.2 モデルの検証
（2）と
（3）
と（5）
のなかで使用した定式が確率的再構築法で評価できた（Zhuang et al.,（1）と．たとえば，余震生成率と時間差分布の再構成に対応する方程式はそれぞれ，
2004）

 
ρ I(mi ∈ [m − Δm, m + Δm])ij ij,
κ̂(m) = ρij I(mi ∈ [m − Δm, m + Δm])

i

 
i

ĝ(t) =

ρij I(ti − ti ∈ [t − Δt, t + Δt])

j

2Δt



i,j

ρij

である．ここで，ρij は式（8）で定義され，I(X) は，X が真であれば 1，そうでなければ 0 の値を持つ指示関数である．
点震源の場合，次の式によって標準化された距離を定義し,



rij =

(xj − xi )2 + (yj − yi )2,Deγ(m−mc )さらに，空間応答の再構成では，式は，

 ifˆR (r) =j

ρij I(rij ∈ [r − Δr, r + Δr])
2Δr



i,j

ρij

で表され，有限震源の場合，余震は有限パッチの中心からの距離に応じて減衰し，再構成式は
次のようになる．

(5.1)

fˆ R (r) =i,ι,jρiιj I(riιj ∈ [r − Δr, r + Δr])
2Δr



i,ι,j

ρij

.

rιj は，トリガとなったイベント j からパッチ ι までの標準化された距離であり，以下の式で定
義される．



(5.2)

riιj =(xj − uiι )2 + (yj − viι )2.
D図2（a）の ETAS の 2 次元点震源モデルと 2 次元有限震源モデルは日本地域の生産性分布，時間差分布，空間分布の再構築結果を示している．したがって，再構成された生産性分布は，理

232

統計数理

第 69 巻 第 2 号 2021

図 2．（a）生産性分布，（b）タイムラグ（本震後の経過時間）分布，（c）震源間距離分布の再構成
結果．κ(M ) の時間分布と空間分布の確率密度関数の理論曲線を点線と破線でそれぞれプロットしたものである．（c）の丸は，有限震源上のパッチまでの距離から再構成した結果である．この図は Guo et al.（2017）の図 3 を修正したものである．

論直線よりも低い A 値と大きい α 値を持つ傾向がある．有限震源モデルでは，この偏差は有（b）では，時間差分布は 2 つのモデルで限震源を持つ大きな地震に対して補正されている．図 2（c）から，空間的なほぼ同じであり，また確率密度は理論曲線と非常によく一致している．図 2
余震の集中性は有限震源モデルの方が大きく，また余震の空間的な減衰が速くなることがわかる．特に，断層形状を取り入れることにより，有限震源モデルの式（3.1）の結果は理論曲線と完全に一致していることが示唆される．
5.3 断層面で生産性密度分布
有限震源 ETAS モデルによって明らかにされたもう 1 つの重要な結果は，余震の生産率密度パターンと本震によって引き起こされた地震時のすべりとの相関関係である．図 3 は，3 つの
データセットからの 3 つの主要な地震の結果を示している．その 3 つの地震は東日本大地震（Guo et al., 2017 を参照），ノルチャ地震（Zhuang et al., 2018），リッジクレスト地震（Guo etである．その中で，最初の 2 つのパターンは 2 次元バージョンの有限震源 ETAS モ
al., 2021）
デルに基づいており，3 番目のパターンは 3 次元バージョンの有限震源 FS モデルに基づいている．
図 3(a) は東北地方太平洋沖地震の結果は，すべりパターンの詳細は大きく異なるが，大きなすべり領域の張り出しや全モーメントの解放などの主要な特徴は類似している．それらを我々の結果と比較すると，断層の浅い部分では生産性が非常に低く，すべりが非常に大きいことがわかった．これは，プレート界面に蓄積された弾性ひずみが既に地震時に解放されているため．したがって，子地震生成率密度の高い場所は，主にと考えられる（Yagi and Fukahata 2011）
すべり量の大きい領域から小さい領域への移行部や，2 つのアスペリティの中間に位置すると（c）の考えられる．結果を図 2 に示す本震のすべり分布は，Dreger（1994）から引用した．図 3
リッジクレスト地震では生産性の高い値が断層面右上に存在している．高い余震生産性を示す
領域は主に 5km から 15km の間に分布しており，中央部では 20km の深さまですべりが広がっている．
図 3 では，生産性の高いパッチが，全体的に高スリップ領域と相補的なパターンを形成している．例えば，図 3(c) に示したリッジクレスト地震のすべりと余震生産性のパターンは，全体的にやや相補的なものとなっている．生産性が最も高い領域は断層の北端であり，ここでは地震時のすべりが比較的低い．

時空間 ETAS モデルの拡張バージョンとその応用

233

図 3．有限震源 ETAS モデルによる推測 1 次余震生産率密度と破裂領域の地震すべりのパターンの比較．（a）2011 年東北地方太平洋沖地震の相対的な生産性分布．等値線は，Yagi and Fukahata（2011）のすべり分布であり，左のカラーバーはすべり量，右側では，生産性を示す．この図は Guo et al.（2017）の図 5 を修正したものである．（b）ノ
ルチャ地震．ゼロから最大までの地震時すべりの値は，虹色の緑から赤までの等高線で示されている．赤い星は本地震の震央位置を表し，青い点はその地域の町や都市の場所を表している．黒い点は，本震の直後に発生した小地震の位置を示している．また，活断層の痕跡もプロットされている．この図は Zhuang et al.（2018）の図 10 から選択さ（左パネル）と余震生産率密度分布（右パネル）．
れた．（c）リッジクレスト地震性すべり
すべりモデルは，Jin と Fialko（2020）を参考にしている．この図は Guo et al.（2021）の図 S1 を修正したものである．

6.

まとめ

いずれの地震に対しても，有限震源モデルでは，震源域の生産性のスケーリングパラメータα が大きくなり，より大きなマグニチュードの地震の誘発能力が向上している．すなわち，FS
モデルでは，本震の破壊域での余震は直接余震として捉えられ，PS モデルでは二次余震として認識される可能性が高い．
日本，イタリア，南カリフォルニアのデータに ETAS モデルを適用し，余震生産性の再構築されたパターンにより，本震の断層面上に生産性の値を補間することが可能となった．全体として，すべりと余震の生産性は相補的なパターンを形成する．このことは，余震が地震後の応力緩和の役割を果たし，応力の変化に対応していることを示す． また，大規模なすべり領

統計数理

234

第 69 巻 第 2 号 2021

域の内部では余震はほとんど観測されておらず，本震の発生時には，このような大規模なす
べり領域は破壊されていることから，本震発生後に発生する余震が限られていることが示唆される．その結果，その後に発生する余震は，大規模なアスペリティの周辺付近で発生する傾向がある．この結果は，これまでの余震分布に関する研究（Das and Henry, 2003; Ebel andと一致している．
Chambers, 2016; Wetzler et al., 2018; Woessner et al., 2006）
新しいモデルの地震活動を予測する際のパフォーマンスについて簡単に説明する．1-day の地震予測に 3 次元点震源 ETAS モデルを適用することは簡単であり，2 次元点震源 ETAS モデルに対して確実で優れたパフォーマンスを提供する．ただし，有限震源 ETAS モデルを適用す
ることは容易ではない．これは，生産率密度の推定に地震系列全体の情報が必要であり，本震の発生直後には利用できないためである．Guo et al.（2021）は，初期余震の位置の情報を使用することにより，生産率密度の空間的変動の大まかな推定が初期余震の場所の予測を改善できることを示した．
謝辞

著者は有益なコメントを寄せてくれた査読者に感謝する．この研究は科研費 19H04073 に
よってサポートされたものである.

参

考

文

献

Bach, C. and Hainzl, S. (2012). Improving empirical aftershock modeling based on additional source
information, Journal of Geophysical Research: Solid Earth (1978–2012), 117(B4).
Console, R. and Murru, M. (2001). A simple and testable model for earthquake clustering, Journal ofGeophysical Research, 106(B5), 8699–8711.
Console, R., Murru, M. and Lombardi, A. M. (2003). Refining earthquake clustering models, Journalof Geophysical Research, 108(B10), 2468.
Daley, D. D. and Vere-Jones, D. (2003). An Introduction to Theory of Point Processes—Volume 1 :Elementrary Theory and Methods, 2nd ed., Springer, New York.
Das, S. and Henry, C. (2003). Spatial relation between main earthquake slip and its aftershock distribution, Reviews of Geophysics, 41(3), 1013.
Dreger, D. (1994). Empirical Green’s function study of the January 17, 1994 Northridge, California
earthquake, Geophysical Research Letters, 21(24), 2633–2636.
Ebel, J. E. and Chambers, D. W. (2016). Using the locations of M>= 4 earthquakes to delineate theextents of the ruptures of past major earthquakes, Geophysical Journal International, 207,
862–875.
Guo, Y., Zhuang, J. and Zhou, S. (2015a). A hypocentral version of the space–time ETAS model,
Geophysical Journal International, 203(1), 366–372.
Guo, Y., Zhuang, J. and Zhou, S. (2015b). An improved space-time ETAS model for inverting the rupture geometry from seismicity triggering, Journal of Geophysical Research, 120, 3309–3323.
Guo, Y., Zhuang, J., Hirata, N. and Zhou, S. (2017). Heterogeneity of direct aftershock productivityof the main shock rupture, Journal of Geophysical Research: Solid Earth, 122(7), 5288–5305.
Guo, Y., Zhuang, J. and Ogata, Y. (2019). Modeling and forecasting aftershocks can be improvedby incorporating rupture geometry in the ETAS model, Geophysical Research Letters, 46(22),12881–12889.

時空間 ETAS モデルの拡張バージョンとその応用

235

Guo, Y., Zhuang, J. and Zhang, H. (2021). Heterogeneity of aftershock productivity along the mainshock ruptures and its advantage in improving short-term aftershock forecast, Journal of Geophysical Research: Solid Earth, 126, e2020JB020494.
Gutenberg, B. and Richter, C. F. (1956). Magnitude and energy of earthquakes, Annals of Geophysics,
9(1), 1–15.
Hainzl, S., Christophersen, A. and Enescu, B. (2008). Impact of earthquake rupture extensions on parameter estimations of point-process models, Bulletin of the Seismological Society of America,
98(4), 2066–2072.
Hawkes, A. G. (1971a). Point spectra of some mutually exciting point processes, Journal of the Royal
Statistical Society: Series B (Statistical Methodology), 33(3), 438–443.
Hawkes, A. G. (1971b). Spectra of some self-exciting and mutually exciting point processes, Biometrika,
58(1), 83–90.
Jin, Z. and Fialko, Y. (2020). Finite slip models of the 2019 Ridgecrest earthquake sequence constrained by space geodetic data and aftershock locations, Bulletin of the Seismological Societyof America, 110(4), 1660–1679.
Marsan, D. and Lengliné, O. (2010). A new estimation of the decay of aftershock density with distanceto the mainshock, Journal of Geophysical Research: Solid Earth (1978–2012), 115(B9).
Nishikawa, T. and Ide, S. (2017). Detection of earthquake swarms at subduction zones globally: Insights into tectonic controls on swarm activity, Journal of Geophysical Research: Solid Earth,
122(7), 5325–5343.
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.
Omori, F. (1894). On the aftershocks of earthquakes, Journal of the College of Science, Imperial
University of Tokyo, 7, 111–200.
Reverso, T., Marsan, D. and Helmstetter, A. (2015). Detection and characterization of transient forcing
episodes aecting earthquake activity in the Aleutian Arc system, Earth and Planetary Science
Letters, 412, 25–34.
Tsuruoka, H., Hirata, N., Schorlemmer, D., Euchner, F., Nanjo, K. Z. and Jordan, T. H. (2012). CSEP
testing center and the first results of the earthquake forecast testing experiment in Japan,
Earth, Planets and Space, 64(8), 661–671.
Utsu, T. (1969). Aftershock and earthquake statistics (I), Some parameters which characterize an aftershock sequence and their interrelations, Journal of the Faculty of Science, Hokkaido University,
Japan, Series, 3(3), 129–195.
Wetzler, N., Lay, T., Brodsky, E. E. and Kanamori, H. (2018). Systematic deficiency of aftershocks
in areas of high coseismic slip for large subduction zone earthquakes, Science Advances, 4(2),eaao3225.
Woessner, J., Schorlemmer, D., Wiemer, S. and Mai, P. (2006). Spatial correlation of after-shock
locations and on-fault main shock properties, Journal of Geophysical Research: Solid Earth,111(B8), B08301.
Yagi, Y. and Fukahata, Y. (2011). Rupture process of the 2011 Tohoku-oki earthquake and absoluteelastic strain release, Geophysical Research Letters, 38(19), L19307.
Zhuang, J. (2011). Next-day earthquake forecasts by using the ETAS model, Earth, Planet, and Space,63, 207–216.
庄 建倉, 尾形 良彦 (2021). 統計地震学における ETAS モデル - その進展とホークス型モデル, 統計数理,
69(2), 145–163.

236

統計数理

第 69 巻 第 2 号 2021

Zhuang, J., Ogata, Y. and Vere-Jones, D. (2002). Stochastic declustering of space-time earthquake
occurrences, Journal of the American Statistical Association, 97(3), 369–380.
Zhuang, J., Ogata, Y. and Vere-Jones, D. (2004). Analyzing earthquake clustering features by using
stochastic reconstruction, Journal of Geophysical Research, 109(3), B05301.
Zhuang, J., Chang, C.-P., Ogata, Y. and Chen, Y.-I. (2005). A study on the background and clustering seismicity in the Taiwan region by using point process models, Journal of GeophysicalResearch: Solid Earth, 110(B5), 1978–2012.
Zhuang, J., Murru, M., Falcone, G. and Guo, Y. (2018). An extensive study of clustering features ofseismicity in Italy from 2005 to 2016, Geophysical Journal International, 216(1), 302–318.
Zhuang, J., Murru, M., Falcone, G. and Guo, Y. (2019). An extensive study of clustering features ofseismicity in Italy from 2005 to 2016, Geophysical Journal International, 216, 302–318.

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

237

Extended Versions of the Space-time ETAS Model and Their Applications
Yicun Guo1 and Jiancang Zhuang2,3
1 Key Laboratory of Computational Geodynamics, University of Chinese Academy of Sciences
2 The Institute of Statistical Mathematics

3 Department of Statistical Science, School of Multidisciplinary Sciences,

The Graduate University for Advanced Studies (SOKENDAI)

This paper summarizes several extended versions of the space-time Epidemic-type
Aftershock Sequence (ETAS) model. The ETAS model is commonly accepted as a standard model in statistical seismology, and its extended versions include the 2D finite source(FS) model, the 3D hypocentral model, and the 3D-FS model. The finite source model
incorporates the rupture geometries of large earthquakes, while the 3D hypocentral model
considers focal depths of earthquakes. The relevant algorithms for model estimation,
stochastic declustering, and earthquake simulations are given. Additionally, the resultsfrom previous applications of the finite source and 3D ETAS models to Japan, Italy, andSouthern California are summarized. The 2D-FS and 3D-FS models enhance the productivity of mainshocks by yielding a larger α value than the point source model. Comparing
the aftershock productivity with the coseismic slip suggests that large slip areas are depleted of aftershocks, and the trajectory of the productivity pattern on the fault planedemonstrates the apparent compensation to the slip.

Key words: ETAS model, probability density, residual analysis, earthquake triggering, aftershock, stochastic
reconstruction.