統計数理（2009）
第 57 巻 第 1 号 179–193
c 2009 統計数理研究所
［研究詳解］

地震活動のクラスタリングおよび除群化の
統計モデル
庄 建倉

†

（受付 2008 年 7 月 29 日；改訂 2008 年 11 月 5 日）

要

旨

本稿では地震クラスタのモデリングと除群化に関連する手法に関して解説する．地震の群はETAS（Epidemic-Type Aftershocks Sequence）
モデルによってうまく説明される．このモデルでは個々の地震は，それが他の地震と関係なく起きたもの（常時地震）
であれ，あるいはそれ以前の地震によって誘発されたもの（誘発地震）
であれ，ある種の確率的法則にしたがって未来の地
震を誘発すると考える．確率的除群化法はこのモデルの加法的特性を利用して開発された手法である．この方法を用いることにより，カタログに掲載された個々の地震が常時活動のものであるのか，あるいは以前生起した特定の地震により引き起こされたものであるのかを確率的に識別することが可能である．推定されたこれらの確率により地震の群や常時地震活動に関係する性質の仮説を検証することができる．
キーワード： 地震，ETAS モデル，クラスタリング，除群化，二次残差解析．
1.

はじめに

John Milne と James Ewing および Thomas Gray が日本での研究活動を行っていた 1892 年に最初の地震計を開発して以来，すでに 100 年以上を経過した現代においても，地震の予測は人類に残された科学的難問の 1 つであり続けている（Bolt, 1993）．地震計の発明は地震の近代的かつ定量的な研究の端緒となった業績である．世界中に配置された地震観測所から集められた地震データの研究から良く知られているように，世界中の地震の震源の分布は決して一様ではない．
図 1 に示すのは 1963 年から 1998 年までに発生したマグニチュード 3.5 を超える地震の震源地マップである．この図から分かるように大部分の地震は次の 2 つの地震帯に属している：その一つは環太平洋地震帯であり，世界中の地震エネルギーの 75％がここで放出される．他方はユーラシア大陸の南側の端に沿って伸びる山岳地帯であるアルパイド（Alpide）地震帯であり，
地震エネルギーの 23％ がこの地域で放出される．それ以外の地震活動の活発な地域には大洋の中間に横たわる稜線や大陸棚の地震活動帯などが含まれる（全体の 2％未満）．これらの地震
帯をすべて考慮すると，地震に関係する地殻は自然に幾つかのプレートに別れていることが分かる．日本は地形的に 4 つのプレートにまたがっている：東側で太平洋プレートに接し，南側
でフィリピンプレート（太平洋プレートのサブプレート）に接し，西側でユーラシアプレート，
北側で北米プレートに接している．このように複雑なテクトニック構造が日本を世界中でも最†統計数理研究所：〒106–8569 東京都港区南麻布 4–6–7

統計数理

180

図 1.

第 57 巻

第 1 号 2009

1963 年から 1998 年までに世界中で発生したマグニチュード 3.5 以上の地震の分布．
データ提供：国立地球物理データセンター，国立地震情報センター（いずれも米国）．

図 2. （a）本州の中央および西部地域で生起したマグニチュード 4.2 を超える地震の発生位置．
（b），（c）はそれぞれ関東地域および京都 大阪 神戸地域で生起した地震の累積頻度曲線．

も地震活動の活発な地域としている．
図 2 に示すのは本州の西部と中央部およびその近傍で発生した深さ 100 km 未満の浅い震源を持つ地震の震源地マップであり，地震活動が空間的および時間的にクラスタリングされていることを示している．これらのクラスタは地域ごとに様相が異なるが，いずれもが両極端を示
す 2 つの空間的地震クラスタ（余震や群発活動が経時的に漸次減衰を示すクラスタと，同じ場
所で時間的に変化しない活動を示す常時地震活動のクラスタ）の性質を示す．しかし，基本的には常時地震活動が余震活動よりも優位性を示す．図 2（b）に示すのは房総地域で起こった地震の累積回数であり，地震の回数が時間と共にほぼ一定の割合で増加していることが分かる．
一方，図 2（c）では地震の累積回数に段差が生じており，これは神戸大地震直後からの余震活動による．
このような相違点は科学的観点から重要である．ゾーニング（分帯）や地震発生能力評価と地震活動のクラスタリングおよび除群化の統計モデル

181

いった長期的な地震予測のために，研究者達は余震などの一過性のクラスタリングを排除してその背後に一貫して存在する常時地震活動を評価しようと努めている．しかし，その一方で，短期的またはリアルタイム予測のためには地震クラスタを良く理解することは欠かすことができない．
2.

地震活動の経験則

地震クラスタのモデリングに関する考察を始める前に，今日の地震活動研究の主眼点である
2 つの経験則について概要を説明しておく．
大森 宇津（Omori-Utsu）
則：19 世紀末，日本における地震学の先駆者である大森（以下 Omori）
が何種類かの大規模余震系列について行った詳細な研究（Omori, 1894）をもとに余震の発生頻度 λ(t) が近似的に地震発生後の経過時間 t に逆比例することを示した．すなわち
λ(t) ≈ K/(t + c).

ここに，K は本震を特徴付ける定数である．その後に宇津（以下 Utsu）
が行った研究（たとえばUtsu, 1969 参照）
により更に一般化と改良が行われ，現在では Omori-Utsu 公式として広く認められている：
K
λ(t) =.
(t + c)p
ここに，パラメータ K および c と p は個々の余震系列を特徴付ける定数であるが， c は常に小さな値を示し， p は 1 に近い値を持つ． Omori-Utsu 則に関わる研究とその歴史についての
詳細は Utsu et al. による総説（1995）を参照されたい．
Gutenberg-Richter 則：地震のマグニチュードが定義された後で公式化された法則である．
この公式は，地震のカタログの中で，ある特定の下限マグニチュード値 m0 を超えるマグニ
チュードを持つ，欠測のありえない規模の地震について次の式が成立することを示している，
log N (≥ m) = a − b(m − m0 ).

ここに，N (≥ m) はマグニチュード m を超える地震の回数を表し，定数 b は一般に Gutenberg（Gutenberg and Richter, 1954）．確率表現を用いてこの式を変形す
Richter の b 値と呼ばれる
ると，次に示す等値な式が得られる，
Pr{M ≥ m} ∝ 10−bm = e−βm .

ここに M は地震規模を表す確率変数であり，β = b log e 10 である．すなわち，マグニチュードの経験的分布則は指数分布に従う．
3.

モデルおよび方法

3.1 ETAS モデルの定義と仮定
ETAS モデルは下式で表現される時間変化する地震発生頻度（条件付き強度関数と呼ばれる，
Daley and Vere-Jones, 2003 参照）
を使用する：
λ(t, x, y, m) =lim

δt →0+
δm ,δx ,δy →0

E{N ((t, t + δt ] × (x, x + δx ] × (y, y + δy ] × (m, m + δm ]) | Ht }δt δm δx δy= µ(x, y, m) +ξ(t, x, y, m; ti, xi , yi , mi ).

i:ti <tすなわち，この式は時刻 t に至るまでの観測データ (ti , xi , yi , mi ) が与えられたとした場合の空間 時間 マグニチュード座標 (t, x, y, m) で発生が予測される地震の回数を表現している．ここ

182

統計数理

第 57 巻

第 1 号 2009に Ht は時刻 t に至るまでの（ただし，時刻 t は含まれない）
観測履歴を表し，µ(x, y, m) は定常状態において自発的に生起するマグニチュード m の常時地震発生率， ξ(t, x, y, m; ti, xi , yi , mi )はそれ以前に発生した i 番目の地震からの誘発率をそれぞれ表している．実際的には以下の仮定を置くのが普通である：（1）
過程全体をマグニチュードの頻度分布と地震の発生率との積として考えることができる，すなわち，
λ(t, x, y, m) = λ(t, x, y)s(m).

ここに，次式が成立するものとし，(3.1)

λ(t, x, y) = µ(x, y) +ξ(t, x, y; ti , xi , yi , mi )i:ti <tさらに，s(m) = βe−β(m−mc ) ,m ≥ mcは Gutenberg-Richter 則（指数分布）から導かれる確率密度関数（p.d.f.）である；（2）応答関数
ξ(t, x, y; ti, xi , yi , mi ) は時間と空間的な位置関係に応じて次のように分離して考える
ξ(t, x, y; ti, xi , yi , mi ) = κ(mi )g(t − ti )f (x − xi , y − yi ; mi ).

ここに，κ(m) = Aeα(m−mc ) ,m ≥ mcはマグニチュード m の地震から誘発される地震（子孫地震）の期待数（生産性）を表すポアソン
確率変数であり，

−p
p−1
t
g(t) =, t > 0,1+ccは子孫地震の親地震からの発生時間を表す確率密度関数（p.d.f.）である．また，
−qx2 + y 2
q−1
1
+
(3.2)
f (x, y; m) =,πDeγ(m−mc )Deγ(m−mc )は子孫地震の親地震との相対的位置を表現する確率密度関数であり， mc は下限マグニチュー
ド値である．上に示したのは Zhuang et al.（2004, 2005）および Ogata and Zhuang（2006）による定式化であり，Ogata（1998）による定式化の改良版である．
式（3.1）は地震が時刻 t，位置 (x, y) で発生するリスクを表していると考えることも可能であり，常時活動地震の発生率 µ からの寄与と時刻 t 以前に生起した地震からの誘発率の寄与
ξ(t, x, y; ti, xi , yi , mi ) によって構成されている．
3.2 最尤法
推定強度関数 u(x, y) が与えられたものとし，我々は（3.2）における常時活動地震の発生強度を次のように設定する，µ(x, y) = νu(x, y).

ここに，ν は正の値を持つパラメータである．この表式を用いて次に示す対数尤度関数
 T 
log L(θ) =log λθ (tk , xk , yk ) −λθ (t, x, y)dtdxdy{k:tk ∈[0,T ],(xk ,yk )∈S}

0

Sを最大化させることにより最尤推定値（MLE），θ̂ = (ν̂, Â, α̂, ĉ, p̂, D̂, γ̂)，が得られる．ここに，右
辺第一項の和は対象領域 S および時間幅 [0, T ] に発生するすべての地震 k に関してとるものとする．この計算の詳細については Ogata（1998）を参照されたい．

地震活動のクラスタリングおよび除群化の統計モデル

183

3.3 間引き法（Thinning procedure）
確率論的除群法における最も重要なポイントの 1 つは（3.1）が決定した点過程の間引き操作（Lewis and Shedler, 1979; Ogata, 1981）である．この操作により過程全体を複数の部分過程に分割することができる．
この方法により，(tj , xj , yj ) において生起した地震 j の応答関数によって，地震 i からの寄与の割合を i 番目の地震が j 番目の地震を誘発する確率として説明することができる．
⎧⎨ ξ(tj , xj , yj ; ti , xi , yi , mi ) , when j > iλ(tj , xj , yj )
(3.3)
ρij =
⎩
0,
otherwise.

さらに，地震 j が先行地震の被誘発地震である確率は次式で表現され，

ρj =
ρij ,
i

従って同じ地震 j が常時活動地震に属する確率は次式で表現される．
(3.4)

ϕj = 1 − ρj =µ(xj , yj )λ(tj , xj , yj )確率 ϕj , ρj ，または ρij で特徴付けられる地震 j を確率的に選択してゆくことにより，それぞれ，地震 i により誘発された過程，クラスタ過程，または常時活動過程といった新たな過程を作り出すことができる．
3.4 常時活動地震の頻度の推定
こうして常時活動過程が得られたならば，何らかの平滑化手段を用いることで常時活動強度を推定することが可能になる．たとえば，本稿では次に例を示す可変幅のカーネル関数による推定を使用する：
1
µ̂(x, y) =Zhj (x − xj , y − yj ).
T j

ここに，T は過程の持続期間を表し，和をとるサブスクリプト j は常時活動過程のすべての地震である．各カーネルの幅 hj は次式で決定される：
hj = max{ε, inf(r : N [B(xi , yi ; r)] > np )}.

ここに ε は小さな実数正値であり，B(x, y; r) は (x, y) を中心とする半径 r の円盤，np は正の整数である．この式で hj は np 番目の最近接地震への距離を表している．
常時活動の強度の平均推定のために間引き手続きとカーネル関数による推定操作を何度も繰り返す代わりに，ここでは重み付けを行うことによって直接的に平均値の評価を行う．すな
わち，
1
(3.5)
µ̂(x, y) =ϕj Zhj (x − xj , y − yj ).
T j
ここに，和をとる j は全部の過程のすべての地震である．同様に，全地震の発生率とクラスタ
地震発生率はそれぞれ
1
(3.6)
M̂ (x, y) =Zhj (x − xj , y − yj )
T j
および
(3.7)

Ĉ(x, y) =1(1 − ϕj )Zhj (x − xj , y − yj )
T j

184

統計数理

第 57 巻

第 1 号 2009で推定することができる．
この式を用い，以下に説明するアルゴリズム（Zhuang et al., 2002）を反復的に使用することにより常時地震活動とクラスタリング構造に付随するパラメータを同時に推定することが可能になる．
3.5 反復法
3.2 節から 3.4 節に説明した方法を利用することにより，常時活動地震の発生率 µ，モデルのパラメータ θ，さらには常時活動の確率 ϕi を評価することができる．但しこれらの 3 つの内 2
つが解っているという前提であるが．殆どのケースでは，我々が知っているのは地震カタログ
のみであり，これらのパラメータ値は得られていない．しかし，以下のアルゴリズムを使用すればこれらの値を同時に評価することが可能である．
アルゴリズム A：確率論的除群化
A1. 固定された np と ε の値，たとえば 5 と 0.05 度（地球表面上の 5.56 km に相当，これは地震の位置誤差に近い値である）が与えられたものとし，各地震 (tj , xj , yj , mj : j = 1, 2, . . . , N )のカーネル関数の幅 hj を計算する．
A2.  = 0, u() (x, y) = 1 をセットする．
A3. 最尤法の計算手続き（Ogata, 1998 参照）を使用して次の条件付強度関数を持つモデル
λ(t, x, y) = νµ() (x, y) +κ(mi )g(t − ti )f (x − xi , y − yi ; mi )i:ti <tを地震データに適合させる．ここに，κ, g ，および f は（3.1）の中に定義された関数である．
A4. j < i (i = 1, 2, . . . , N ) を満たすインデックスについて ρij , ρi ，および ϕi を計算する．
A5. （3.5）から µ(x, y) を計算し，それを u(+1) (x, y) として記録する．
A6. 与えられた小さな正数 ε に対して収束の条件式 max |u(+1) (x, y) − u() (x, y)| < ε が成立しないときは  =  + 1 と書き換えてステップ A3 へ戻る．条件式が成立したときは νu(+1) (x, y)を常時活動地震の頻度と解釈して ρij , ρi ，および ϕi を出力する．
地震データとして我々が使用したのは気象庁が編纂した震源カタログである．その中でも，期間として 1926 から 1995 年，北緯 34 から 39 度および東経 131 から 140 度の範囲の矩形領域に属し，マグニチュード 4.2 以上，深さ 100 km 以内の地震データを選択して使用した
（図 2
（a））．
この領域の MLE は次の通りである：A = 0.020 event/(deg2 . day)，α = 1.135, c = 0.040 day，
p = 1.15，および d = 0.0010/ deg 2 ．図 3（a）（
– d）に示すのは全地震発生強度（3.6），相対クラスタリング発生率（3.7），常時活動地震（3.5），および相対クラスタリング効果 C(x, y)/M (x, y) それ
ぞれの推定値である．図 3（d）に示す常時地震活動効果を見る限り，標準的な余震によって引き起こされるクラスタリング効果が取り除かれていることが分かる．これらは北丹後（1927），鳥取（1943），三河（1945），福井（1948），新潟（1964），能登半島沖（1993），および神戸地震（1995）
な
どの余震域を示している．これとは対照的に和歌山および茨城地域では常時活動地震の頻度と全強度が殆ど同等であり，これは一過的にクラスタリングされていない地震群が空間的に発生していることを示している．上記のアルゴリズムが確率として出力する ϕi および ρij はカタログに登録された地震を，確率の実現によって，誘発関係を示す異なる系統図に分類するために役立つ．このような分類を行うアルゴリズムは以下に説明する順序にしたがって実行される．

地震活動のクラスタリングおよび除群化の統計モデル

185

図 3. （a）研究対象領域に式（3.6）
を適用して得られる全地震の発生率 M（単位：地震数/(deg2× 年)）．（b）研究対象領域におけるクラスタ（子孫）地震発生率 C（単位：地震数/(deg2× 年)，式（3.7））．（c）研究対象領域に式（3.5）を適用して得られる常時地震活動発生率µ（単位：地震数/(deg2 × 年)）．（d）クラスタ化係数（C/M ）．

3.6

確率的再構成アルゴリズム B：地震クラスタの確率的分類
B1. アルゴリズム A を使用して ϕi および ρij を計算する：ここに i = 1, 2, . . . , N, j = 1, 2, . . . , Nであり，N は地震の総数である．
B2. それぞれの地震 j (j = 1, 2, . . . , N ) について，区間 [0, 1] に一様に分布する確率変数 Ujを生成させる．
B3. それぞれの j について次式を考える
k
Ij = max k − 1 : ϕj +
ρij ≥ Uj and 0 ≤ k < j .
i=1

Ij = 0 ならばその j を常時活動地震または初期地震として選択し，それ以外の場合は j 番目の地震を Ij 番目の地震から直接誘発された地震であると見なす．
図 4 は，3.5 節の中で使用される気象庁データの確率的除群化の実現結果の一つを示す．
各系統樹の初期地震が選択されると，それに続いて除群化されたカタログ（常時活動地震を含む）が得られる．こうして作られた除群化後のカタログは従来の除群化法を使用して得られたものとは異なっている：すなわち，ここで選択される地震は常に系統樹内で最初に起こる地震ではあるがクラスタ中の最大地震であるとは限らない．したがって，従来の文献で「本震」と呼ばれる地震が確率的な除群化後のカタログには含まれていないことがある．もし，各クラ

統計数理

186

図 4.

第 57 巻

第 1 号 2009

本州地域にアルゴリズム 2 を適用した結果．（a）全地震の位置．（b）子孫地震の位置．
（c）常時地震の位置．（d）全地震の空間（緯度）時間プロット．（e）誘発（子孫）地震の空間（緯度）時間プロット．（f）常時地震の空間（緯度）時間プロット．

スタ内の最大地震を常時地震活動の代表例として選択することが好ましいのであれば，そのよ
うにすることも可能である．
地震のクラスタリングでは常時活動地震と異種地震クラスタが空間的および時間的に複雑に入り組むことによる特徴的な性格が現れるが，このような問題のクラスタに関する仮説の検証に取り組もうとするのが上記のアルゴリズムである．地震クラスタを分離して異なる確率論タ
イプを得るためにアルゴリズム B を多数回繰り返すことも可能である．得られる結果が決して一意ではないということ自体が地震クラスタ決定における不確実性の存在を如実に示しており，地震活動クラスタリングパターンのある種の特性が持つ重要性を評価する上でアルゴリズ
ムの反復適用が非常に役立つ．これらの検証は確率 ϕj と ρij を直接の対象として実施することも可能である．以下の節ではこれらの確率値を使用して地震クラスタリングに付随する特性を再構築する方法について説明する．そこで利用するのは日本気象庁のカタログに記録されている地震である．また，比較のためにシミュレーションにより生成されたカタログを対象として同様の再構築操作を行うものとする．
4.

応用例：地震クラスタに関連する仮説の検証（3.2）は空間応答関数として使用され
位置分布：空間 時間 ETAS モデルの初期の研究では式ておらず，Ogata（1998）
や Zhuang et al.（2002, 2004）の研究ではその代わりに以下の 2 種類の形式が使用されている．すなわち，モデル I−qx2 + y 2
q−1
exp
−
,
(4.1)
f (x, y; m) =πDeα(m−mc )Deα(m−mc )地震活動のクラスタリングおよび除群化の統計モデル

およびモデル II(4.2)

f (x, y; m) =q−1
πDeα(m−mc )1+x2 + y 2Deα(m−mc )

187

−q
,

である．この 2 つの関数のどちらが地震活動のモデリングにより適しているか，すなわち余震
活動はモデル I が示すように短い空間的範囲内で減衰してゆくのか，それともモデル II のように長い距離で減衰していくのかが重要な問題である． Ogata（1998）によれば，AIC モデル選択
手順を使用した場合は 2 番目のモデルの方が第 1 のモデルよりもデータへの適合性が優れている．我々も，以下のような単純なヒストグラム法を用いることにより同じ結論に達した．
誘発された地震 j とその直接の先祖（インデックス i）との間の変換距離を次式によって定義
する，(xj − xi )2 + (yj − yi )2rij =.
Deα(mi −mc )
式
（4.1）
，
（4.2）
から，モデル I およびモデル II の rij はそれぞれ次の密度関数にしたがう：すな
わち，
2fR (r) = 2re−r ,(4.3)

r ≥ 0,

および

2r(q − 1), r ≥ 0,(1 + r 2 )qである．式（4.3）で定義される密度を持つ分布は Rayleigh 分布と呼ばれる．他方，fR (r) のノ
ンパラメトリックな推定量として次式のような再構成が可能である：
∆ri,j ρij I |rij − r| <2.
fˆR (r) =∆r i,j ρijfR (r) =ここに，∆r は小さな値を持つ正数である． 2 つのモデルそれぞれについて fˆR (r) と fR を比較した結果を図 5 に示す．この結果から，モデル I を使用した場合であっても先祖と直接の子孫との間の変換距離を使用して再構築した確率密度は理論確率密度と非常に異なっていることが分かる．一方，モデル II を使用した場合の再構成確率密度は理論値に非常に近い値を示した．

図 5.

規格化した誘発距離 fR (r) の経験分布（円）にモデル （I 左図）およびモデル II（右図）を適用して分布を再構築した結果を示す．（左図）と（右図）において実線で示すのが fR (r)の理論曲線である．

統計数理

188

第 57 巻

第 1 号 2009

これらの結果は余震が，空間的に狭い範囲というよりも，空間的に長い距離にわたって減衰するという事実を確認するものである（Ogata, 1998; Console et al., 2003）．同時に，これらの結果は上に説明した再構成手法の堅牢性を示唆するものでもあり，たとえモデル I のような不適切と思われるモデルを使用したとしても，モデル II を使用した場合の関数に類似した再構成
確率密度関数が得られることを示している．
モデル II の方がモデル I よりもはるかに優れた地震活動データへのフィッティングが得られたことから，以下の節における再構成の考察では専らモデル II を使用するものとする．
常時活動地震と誘発地震の特徴比較
空間 時間 ETAS モデルにおける最も重要な仮定のひとつは常時活動地震と誘発地震間に実質的な差異はないという考え方である．すなわち，一旦ある地震が生起したならば，それが常時活動地震により生起したか，あるいは別な時間的に先行する地震によって生起したかを問わず，マグニチュードはすべての地震に共通なマグニチュード分布から導かれ，他のすべての地震と同様に子孫地震を誘発すると考える．以下のセクションで我々は常時活動地震と誘発地震との間にマグニチュード分布や誘発能力，子孫地震の生起時間と空間的分布などに差異があるか否かの検証を行う．
気象庁カタログに登録されている地震とシミュレーションで作成したカタログを対象として
全地震，常時活動地震と誘発地震の誘発能力を再構築（図 6 参照）
することができるが，そのときに使用するのは次の 3 つのノンパラメトリックな推定量である：
∆mi
j ρij I |mi − m| <2
κ̂(m) =
(4.4)
,
∆m


∆m
ϕ
ρ
I
|m
−
m|
<
i
i
j i ij2(4.5),
κ̂b (m) =∆m i ϕi

図 6.

左図：再構成された生産性関数（同じマグニチュードの事象が誘発する子事象の平均数を表す）．式（4.4），（4.5）および（4.6）に対応してすべての事象（円），常時地震（緑色の“X”）および誘発地震（赤色 + ）をそれぞれ表している．直線はモデル（3.1）を 気象庁
カタログに適用して最尤推定を行った結果を示す．右図：ETAS モデルを使用してシミュレーションを行った合成カタログを対象として同じ処理を行った．

地震活動のクラスタリングおよび除群化の統計モデル

および，
i

(4.6)

κ̂c (m) =

189



∆m
(1
−
ϕ
)ρ
I
|m
−
m|
<
i ij
i
j
2
.
∆m i (1 − ϕi )常時活動地震と誘発地震は両者ともに近似的に指数法則（指数係数は異なるが）にしたがって子孫地震を生起させる．先祖地震のマグニチュードが同じであるとすれば，誘発地震は常時活動
地震よりも多くの子孫地震を生起させる．地震のマグニチュードが大きくなるほど，常時活動
地震と誘発地震における誘発能力の差異は小さくなってゆく．
これらの差異は数値的な誤差により生ずるものではない．他方シミュレーションで作成したカタログからの結果によれば 2 種類の地震間で誘発能力の差は認められない．
5.

理論的基礎

二次残差解析

この節では，まず最初に予測可能過程という概念を導入するとともに，過程の条件付強度関数に関係する 2 つの定理を導入する．
点過程 N に付随する予測可能な関数 f は時間を変数とする関数であり，その値は時刻 t よりも過去の N の発生履歴により決定される．すなわち，N = {ti : i = 1, . . . , n} とすれば f (t, N )は t および過去の発生履歴 {tk : tk < t} のみの関数であり，将来の発生時刻 {tk : tk ≥ t} には依存しない．たとえば，時間的点過程の条件付強度 λ(t) は予測可能であるが， N (t) = N [0, t] は予測可能でない．N (t) が予測可能でない理由は，N (t) が t 以前ばかりではなく時刻 t における N の生起に関する情報にも依存するためである．混同の恐れがないことから，f (t, N ) を簡単のために f (t) と表記する．
以下の定理は点過程の条件付強度がマルチンゲール性を持つことを示している：
定理 1. 単純時間的点過程 N が条件付強度 λ(t) で表現されるとすれば，任意の予測可能な関数 f および実数の正則（ルベーグ測度が存在する）部分集合について次の関係が成立する（ただし，下式のどちらか一方の平均が存在しなければならない）
：
⎡
⎤



E⎣
f (ti )⎦ = E
f (t)λ(t)dt .
ti ∈N∩S

S

上記の特性は二次モーメント特性への一般化が可能である．
単純時間的点過程 N が条件付強度 λ(t) で表現され，かつ h1 (s)h2 (t)定理 2. （Zhuang, 2006）の様式の単調限界を持つ関数の線形結合として関数 f (s, t) が表現されるとすれば（h1 , h2 はい
ずれも予測可能な関数とする），ユークリッド平面上の正則（ルベーグ測度が存在する）部分集
合 D について次の関係式が成立する：





E
f (ti , tj )I(i = j ∧ (ti , tj ) ∈ N × N ∩ D) = E
f (s, t)λ(s)λ(t)ds dt .
i,jD

ただし，上式において少なくとも一方の平均が存在するものとする．
上記の 2 つの定理は ETAS モデルのような時空的点過程についても成立する．
ここで，論点を再び常時地震活動と誘発地震活動の差異検出の問題に戻す．次式に示す条件付強度によって表現される，より複雑なモデルを考えるものとする，(5.1)

λ1 (t, x, s) = λ1 (t, x, s, ω)I(ω = 0) + λ1 (t, x, s, ω)I(ω = 1).

統計数理

190

ただし

第 57 巻

第 1 号 2009⎧⎪
⎨ u(x)s0 (m),
λ1 (t, x, m, ω) = s (m)ξ(t, x; ti , xi , mi , ωi ),
1
⎪
⎩

(5.2)

if ω = 0,
if ω = 1,

ti <t

および，



(5.3)

ξ(t, x; ti , xi , mi , ωi ) =κ0 (mi )g0 (t − ti )f0 (x − xi , mi ),if ωi = 0,

κ1 (mi )g1 (t − ti )f1 (x − xi , mi ),if ωi = 1.

である．ここで，ω = 0 は常時活動地震，ω = 1 は誘発地震を表す．この混合型モデルでは常時活動地震 (t, x, m) が確率密度関数 s0 (m) から導かれるマグニチュード s を持ち，平均として
κ0 (s) 個の子孫地震を誘発する．これら子孫地震のマグニチュード，時間および位置分布はそれぞれ γ1 , g0 および f0 によって特徴付けられる．誘発地震 (t, x, m) のマグニチュードは密度関
数 s1 (m) から導かれ，平均として κ1 (m) 個の子孫地震を誘発する．これらの子孫地震のマグニチュード，時間，および位置密度はそれぞれ s1 ，g1 ，f1 である．
κ0 値の「比率が不偏」な推定をするために次の 2 つの式を導入する：
s1 (m )ξ(t , x ; t, x, m, ω)λ1(t, x, m, ω)I(ω = 0)
λ1 (t , x , m )λ1 (t, x, m)
s1 (m )κ0 (m)g0 (t − t)f0 (x − x, m)u(x)s0 (m)(t, x, m; t , x , m ) =λ1 (t , x , m )λ1 (t, x, m)

H(t, x, m; t, x , m ) =および
λ1 (t, x, m, ω)I(ω = 0) u(x)s0 (m)=.
λ1 (t, x, m)
λ1 (t, x, m)
次に f1 ，f2 ，g1 ，および g2 を積分すると 1 が得られることに着目すると次の 2 式が得られる：
(5.4) EH(ti , xi , mi ; tj , xj , mj )I(mi ∈ [m0 − δ, m0 + δ])h(t, x, m) =i,j 
=E
H(t, x, m; t, x , m )I(m ∈ [m0 − δ, m0 + δ])dt dx dm dt dx dm= κ0 (m)s0 (m)I(m ∈ [m0 − δ, m0 + δ])dm ×
u(x)dt dx
u(x)dt dx≈ 2δκ0 (m0 )s0 (m0 ) ×

および
(5.5)


E


i


h(ti , xi , mi )I(mi ∈ [m0 − δ, m0 + δ])h(t, x, m)I(m ∈ [m0 − δ, m0 + δ])dt dx dm=E= s0 (m)I(m ∈ [m0 − δ, m0 + δ])dm ×
u(x)dt dx
u(x)dt dx.
≈ 2δs0 (m0 ) ×

式
（5.4）
および
（5.5）
から得られる「不偏比率」（ratio-unbiased）κ0 (m) 推定値は次式で与えられる：κ̂0 (m0 ) =i,j H(ti , xi , mi ; tj , xj , mj )I(mi ∈ [m0 − δ, m0 + δ])i h(ti , xi , mi )I(mi ∈ [m0 − δ, m0 + δ]).

地震活動のクラスタリングおよび除群化の統計モデル

191

(1)
(0)
(1)
最初のステップにおいて γ0(0) = γ0(1) = γ̂ ，κ(0)= f0 = fˆ0 の条件が成立するも0 = κ0 = κ̂，f0のとし，γ̂ ，κ̂，fˆ，ĝ が ETAS モデルのフィッティング後の条件付強度から導かれるものとす
れば次式が得られる：(1)

κ̂0 (m0 ) =i,j ϕ̂i ρ̂ij I(mi ∈ [m0 − δ, m0 + δ])i ϕ̂i I(mi ∈ [m0 − δ, m0 + δ]).

ここに，ϕ̂i および ρ̂ij は式（3.3），（3.4）で定義される ϕi , ρij の推定値である．同様に次式が成立することが分かる：
κ̂1 (m0 ) =
i,j H (ti , xi , mi ; tj , xj , mj )I(mi ∈ [m0 − δ, m0 + δ])i [1 − h(ti , xi , mi )]I(mi ∈ [m0 − δ, m0 + δ]).

ここに，次の 2 つの式が成立する
s1 (m )ξ(t , x ; t, x, m, ω)λ1(t, x, m, ω)I(ω = 1)
λ1 (t , x , m )λ1 (t, x, m)
u(x)s0 (m)s1 (m )κ1 (m)g1 (t − t)f1 (x − x, m)
1−(t, x, m; t, x , m ) =λ1 (t , x , m )λ1 (t, x, m)

H  (t, x, m; t, x , m ) =

および
(1)

κ̂1 (m0 ) =i,j (1 − ϕ̂i )ρ̂ij I(mi ∈ [m0 − δ, m0 + δ])i (1 − ϕ̂i )I(mi ∈ [m0 − δ, m0 + δ]).

式（5.1）から（5.3）に現れる他の関数（g1 ，g2 ，f1 ，f2 など）についても同様の方法で導き出すことができる．
6.

まとめ

本稿では地震クラスタの理解に必要な ETAS モデルに付随する何種類かの有力な技法について解説した．これらの技法の基本となる考え方は，複数の地震クラスタ区分に付随する不確実性を確率で表現し，しかる後にクラスタの再構築を行い，モデルによるシミュレーションデー
タの比較を行うことにより，実データとモデルの差異が明瞭になるような方法である．後半の説明で示したとおり，この方法の理論的基礎となるのは残差解析である．確率的な除群化法は，前震を含む地震のクラスタに関連するあらゆる仮説の検証に関して，従来の主観的な除群法にとって代わって，広く使用されていくであろう．
謝

辞

本研究所の尾形良彦教授から，この論文を書いた間，私に大きな支援をいただきました。審
査員の先生方ならびに編集委員の中村隆教授には，有益なご助言と修正のご指摘を賜りました．
記して深い感謝の意を表します．

参 考 文 献Bolt, B. A.（1993）. Earthquakes and Geological Discovery, Scientific American Library from W. H.
Freeman and Co., Salt Lake City, Utah.
Console, R., Murru, M. and Lombardi, A. M.（2003）. Refining earthquake clustering models, Journal（B10），2468, doi:10.1029/2002JB002130.
of Geophysical Research, 108
Daley, D. and Vere-Jones, D.（2003）. An Introduction to the Theory of Point Process, Vol. 1, SpringerVerlag, New York.

統計数理

192

第 57 巻

第 1 号 2009

. Seismicity of the Earth and Associated Phenomena, 2ndGutenberg, B. and Richter, C. F.（1954）
ed., Princeton University Press, Princeton, New Jersey.
. Simulation of non-homogeneous Poisson processes by thinLewis, P. A. W. and Shedler, E.（1979）
ning, Naval Research Logistics Quarterly, 26, 403–413.
Ogata, Y.（1981）. On Lewis’ simulation method for point processes, IEEE Transactions on Information Theory, IT-27, 23–31.
. Space-time point-process models for earthquake occurrences, Annals of the InstiOgata, Y.（1998）
tute of Statistical Mathematics, 50, 379–402.
. Space-time ETAS models and an improved extension, TectonoOgata, Y. and Zhuang, J.（2006）
physics, 413, 13–23.
. On after-shocks of earthquakes, Journal of the Faculty of Science, University of
Omori, F.（1894）
Tokyo, 7, 111–200.
Utsu, T.（1969）. Aftershock and earthquake statistics（I）: Some parameters which characterize anaftershock sequence and their interrelations, Journal of the Faculty of Science, Hokkaido Uni, 3, 129–195.
versity, Series. VII（Geophysics ）Utsu, T., Ogata, Y. and Matsu’ura, R. S.（1995）. The centenary of the Omori formula for a decay
law of aftershock activity, Journal of Physics of the Earth, 43, 1–33.
. Second-order residual analysis of spatiotemporal point processes and applications
Zhuang, J.（2006）
in model evaluation, Journal of the Royal Statistical Society, Series B（Statistical Methodology ）, 68, 635–653. doi: 10.1111/j.1467-9868.2006.00559.x.
Zhuang, J., Ogata, Y. and Vere-Jones, D.（2002）. Stochastic declustering of space-time earthquake
occurrences, Journal of the American Statistical Association, 97, 369–380.
. Analyzing earthquake clustering features by using
Zhuang, J., Ogata Y. and Vere-Jones D.（2004）
stochastic reconstruction, Journal of Geophysical Research, 109, B05301, doi:10.1029/2003JB
002879.
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 a point process model, Journal of GeophysicalResearch, 110, B05S18, doi:10.1029/2004JB003157.

Proceedings of the Institute of Statistical Mathematics Vol. 57, No. 1, 179–193 (2009)

193

Statistical Models for Earthquake Clustering and Declustering
Jiancang Zhuang
The Institute of Statistical Mathematics

This review paper summarizes the statistical methods associated with modeling earthquake clusters and declustering. Earthquake clusters are well described by the epidemictype aftershocks sequence (ETAS) model. In this model, each earthquake, whether it isfrom the background or triggered by a previous earthquake, triggers its own child events
independently according to some probability rules. The stochastic declustering method isdeveloped by making use of the additive property of this model. With this method, each
earthquake in the catalog can be identified to be a background event or be triggered froma particular previous event in probabilities. These estimated probabilities can be used totest hypotheses associated with seismicity clustering or background.

Key words: Earthquake, ETAS model, cluster, declustering, second-order residual analysis.