統計数理（2015）
第 63 巻 第 1 号 45–64
c 2015 統計数理研究所
特集「地震予測と統計モデル」［研究詳解］

地震活動の異常性とモデリング
熊澤 貴雄†（受付 2014 年 12 月 4 日；改訂 2015 年 3 月 30 日；採択 4 月 1 日）

要

旨点過程 ETAS（epidemic-type aftershock sequences）モデルは地震の発生時系列を表現し解析する統計モデルである．点過程の「条件付き強度関数」（conditional intensity function）は地震発生の切迫度を予測する関数であり，地震の確率予測の計算に必須である．また，通常の地震発生の時系列（地震活動）に何らかの異常が介入した場合に，それらを検出する手法を与える．異常性の特徴を捉え，条件付き強度関数に組み込むことは地震予測の精度向上に必要である．本稿では，ETAS モデルによる地震活動の変化点解析を詳説し，異常地震活動をインバージョン解析する非定常 ETAS モデルを解説する．これらの解析例として東北沖地震によって誘発された地震活動を考察する．
キーワード：点過程，条件付き強度関数，ETAS モデル，異常地震活動，非定常 ETAS
モデル．
1.

はじめに

あらゆる地震には誘発性が伴う．顕著に大きい地震（本震）に対して誘発された地震のことを
伝統的に余震と呼ぶ．この誘発による続発性のため，地震発生の系列を時間軸上に並べてみるするので，何らかの理由で発生パターンに変化が生じたとしてと地震が随所で密集（clustering）
も，それを目視で判別するのは一般に困難を伴う．しかしながら適切な点過程モデルを物差しとすることで地震の発生パターンの変化を統計学的に取り扱うことが可能となる．地震発生パ
ターンの変化は地殻内の応力状態の変化を推測するセンサー（sensor）となり得る．点過程モデ
ルは地震発生の確率予測に於いてさまざまに基礎的な役割を担うが，これは本特集の尾形論文に譲り，本稿では地震発生パターンの変化を計る物差しとしての役割を，筆者らの研究を紹介しながら，基本的な活用から非定常モデルを考慮した応用へと順を追って解説する．
点過程は地震発生などの各種災害発生や機械やシステム等の故障発生等の発生時刻の不規則性の確率構造を研究対象とする数学モデルである．これらには時間とともにデータ値が変動する通常の時系列とは異なる確率構造に則した独特の統計解析が必要である．点過程を記述する
方法は幾つか有るが，条件付き強度関数（conditional intensity function）による表記が予測に直結し一般的かつ応用性に優れている．
（intensity function）
λ(t) は非定常ポアソン過程を特
たとえば時刻 t にのみ依存する強度関数
定し，それぞれの時刻での事象の期待発生率を表す．これの特殊なケースとして λ(t) が時間 tに対して定数となる場合は強度 λ の定常ポアソン過程となる．より一般に強度関数が過去の事（conditional intensity
象の発生履歴にも依存する場合が考えられ，λ(t|Ht ) は条件付き強度関数†

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

統計数理

46

第 63 巻

第 1 号 2015

function）と呼ばれる．履歴 Ht として時刻 t 以前の過去の事象発生の発生時刻やマグニチュー
ド列のデータを考えるのが普通であるが，利用できる他の観測データも含めると一層予測精度が高まることもある．これらのデータが与えられた場合，点過程の発生時刻に関する確率法則
は条件付き強度関数によって完全に特定（Liptzel and Shiryaev, 1978）されるので，点過程のモデリングとは即ち条件付き強度関数のモデリングと言い換えても良い．条件付き強度関数からへと直につなげることができる．
点事象の生成（simulation）
地震発生系列は 1960 年代から 70 年代にかけて，その発生が時間軸上の事象の点過程として初めてモデリングされた（Vere-Jones and Davies, 1967; Vere-Jones, 1970）．その後，条件付き
強度関数（Hawkes, 1970; Hawkes and Adamopoulos, 1973）が出現して以来この分野は目覚まし
く発展しており，特に条件付き強度関数にもとづく最尤法やシミュレーション法が展開された．たとえば 1980 年代に開発された ETAS モデル（Epidemic-Type Aftershock（Ogata, 1978, 1981）
は地震発生時系列（地震活動）の標準的な統計モデル
Sequence model; Ogata, 1985, 1988, 1989）
として地震活動現象の様々な研究仮説検証や予測のために使用されている．
2.

定常 ETAS モデル

広域で長期間に渡る地震発生系列において，群発地震や二次余震（余震の余震）
等の地震間の関係を考慮すると，本震と余震の区別は明確ではない．そこで，全ての地震がそのマグニチュードに応じた規模で付随する誘発地震活動を持つという仮定のもとに展開されたのが ETAS モデ
ルである．このモデルは以下の条件付き強度関数で定義される．すなわち(2.1)

λθ (t|Ht ) = μ +
{i: S<ti <t}K0 eα(Mi −Mc ).
(t − ti + c)p

ここで ti は i 番目の地震が発生した時刻であり，そのマグニチュードを Mi で表す．従って Ht
は時刻 t 以前の地震の発生時刻とそのマグニチュードの履歴である．Mc は基準のマグニチュー
で，一般の地震活動の解析ではモデルに取り込む地震データの様々なド（reference magnitude）
で与えられることが多いが，余震列の解析では本震下限マグニチュード（threshold magnitude）のマグニチュードを採り，異なった下限マグニチュードに従ってどのようにパラメータが変化するかなど，目的に応じて採用する．いずれにしても Mc に依存してパラメータ K0 値の推定値が違ってくる．
（2.1）式では時刻 t 以前に発生した各地震 {i : ti < t} の余震活動の総和として時刻 t 以降における誘発効果が第 2 項に表現される．各々の地震の余震活動は本震からの時間経過 t における発生強度を表す大森・宇津の経験則（改良大森公式，Utsu, 1961）(2.2)

λθ (t) =K(t + c)pに基づいて推移する．時間の逆ベキ関数であるため，余震活動は直後は急激に減るが，その後には長々と尾を引く減衰をし，濃尾地震などは一世紀を過ぎても余震が確認されている（Utsu etal., 1995）．
ETAS モデル（2.1）
は 5 つのパラメータ θ = (μ, K0 , c, α, p) で特徴づけられる．この内で p 値はの減衰パラメータ p 値に対応するが，余震系列によっては互いに少々違った
大森・宇津公式（2.2）．c 値は余震の時間減衰がベキ関数に乗るまでの時間のスケー値をとる（Guo and Ogata, 1997）
ルを表す．c 値は本震直後の余震の欠測に大きく影響されるが，本来の物理的な意味については．パラメータ α は誘発地震のマグニチュード
余震研究での対象になっている（Utsu et al., 1995）に対応する誘発された地震（余震）の期待個数を左右する量である．α 値が十分大きい場合には

地震活動の異常性とモデリング

47

ETAS モデルは殆ど本震の影響だけを受け，大森・宇津の式（2.2）と同等になる．他方，群発地
震活動においては通常小さな α 値をとる．パラメータ K0 は同じマグニチュード値での余震の生産効率の地域性を示し，本震直後の余震予測に大きな役割を果たす（本特集の近江論文参照）．
μ は常時地震活動率である．定数の場合は μ は，ずっと過去の数々の地震による現在に至る迄
の余震の減衰が重ね合わさった発生率の成分を表すと考えて良い．したがって十分長い先駆期の履歴データが利用できれば推定 μ 値は小さくなる．断層内に応力が蓄積す間（次 2.1 節参照）
る速度と関係するという報告もある（Ide, 2013）．最後の二つのパラメータの時間変化を 4 節で詳しく見ていく．
2.1 パラメータ推定
ETAS モデルなどの点過程の条件付き強度関数のパラメータ値をデータから推定する手法は最尤法が一般的である．尤度関数を最大化する数値計算アルゴリズムでパラメータの最尤推定値が得られる（Ogata, 2006, ソフトウェア SASeis2006 参照）．
期間 [0, T ] 内に n 個の地震が時刻 t1 < t2 < · · · < tn に生起しているとする．対応するマグニ
チュードのデータは {Mi : i = 1, . . . , n} であるとする．n 個の事象が上記の指定時刻のそれぞれに発生するという同時確率密度関数は，点過程の条件付き確率の計算により，
 n

  T


(2.3)
fT (t1 , t2 , . . . , tn ; θ) =λθ (ti |Hti ) exp −λθ (t |Ht ) dt0

i=1の様に条件付き強度関数 λθ (t|Ht ) で表現される．条件は履歴 Ht = {(ti , Mi ); ti < t} である．
λθ (t|Ht ) は θ によってパラメータ化されているので（2.3）の同時確率分布の密度関数も θ によって特徴付けられる．観測された事象データ t1 < t2 < · · · < tn を（2.3）に代入した θ の関数(2.4)

L[0,T ] (θ|HT ) = fT (t1 , t2 , . . . , tn ; θ)を尤度と呼ぶ．対数尤度(2.5)

log L[0,T ] (θ|HT ) =
{i: 0<ti <T }log λθ (ti |Hti ) − T
0log λθ (t |Ht ) dt

を最大にする θ をもって統計モデルを推定する手法が最尤法と呼ばれる．
ETAS モデルの地震列への適合については以下の注意を要する．地震データとして典型的には余震列と広域の地震列の 2 つのケースが考えられる．これらに適用するにあたってパラメー
タ推定の偏りを招きかねない共通する問題がある．先ず，余震列を解析するときに本震直後の余震の欠測を考慮する必要がある．次に広域の地震活動の場合，過去の地震履歴からの誘発率の影
響を無視できない．とくに，当該地域や近辺に大地震が過去にあった場合などである．これらの問題に対処するために，パラメータの推定に当たってはデータを先駆期間（precursory period）[0, S) と推定期間（target period）[S, T ] に分けて考える必要がある．推定期間ではデータの下限
マグニチュード Mc 以上の地震が漏れなくとれている必要があるが，先駆期間では十分大きな地震や余震が含まれていればよい．かくして ETAS モデルのパラメータは対数尤度
 T

(2.6)
log L[S,T ] (θ|HT ) =log λθ (ti |Hti ) −log λθ (t |Ht ) dt
{i: S<ti <T }Sを最大化する．ただし履歴 Ht は先駆期間の地震も含む．先駆期間のデータを全く無視するとパラメータの推定に理論上偏りが出る恐れがある．

48

統計数理

第 63 巻

第 1 号 2015

2.2 診断解析
推定したパラメータのモデルがデータの性質を再現しているかどうかを確認する為には，モデを行ってデータの性質とルから多量の点過程シミュレーション（例えば Ogata, 1981, 1988 参照）
比較することは良く使われる手段の一つであるが，以下に説明する診断解析（diagnostic analysis）
もある．推定モデルによる予測データと元のデータを比較し外れ度合いを予測「残渣」（residual）と呼ぶ．観測値から予測値を算術的に差し引いた通常の時系列解析の「残差」（residual）とは異なり，点過程に於ける残渣は推定した条件付き強度関数に基づく時間伸縮により以下のように定義される．いま，条件付き強度関数 λθ (t|Ht ) が得られたとする．この強度関数の，当てはめ開始時刻 S から期間内の任意の時刻 t までの時間積分
 t
(2.7)
Λθ (t|Ht ) =λθ (u|Hu )duSは与えられたモデルの下での，時刻 S から時刻 t までの事象の累積発生数の期待値を表す．
そこで（2.7）式にもとづいて通常時間 t からの時間変換 τ = Λθ (t|Ht ) を考えると，観測事象列
t1 < t2 < · · · < tn は τ1 < τ2 < · · · < τn へと一対一に変換される．この変換時間上での事象列を残渣過程，このときの変換時間を残渣時間と呼ぶ．観測事象列が当のモデルから生成された（Papangelou,
場合は変換事象列 τ1 < τ2 < · · · < τn は標準ポアソン過程になることが証明される．言い換えるとモデルが真の場合は残渣時間上での観測データの累積関数は傾き 1 の直線
1972）に乗った一様分布になる．データ解析に於いて実際の真のモデルは未知であるが，（2.7）式に最尤推定値を代入して，この時間変換を行ったときに累積関数が傾き 1 の直線に沿った分布をしていれば，モデルは真の過程を良く近似していると言える．以上の変換に関する誤差解析につなどを参照されたい．
いては例えば Ogata（1992, Appendix）
2.3 応用例（長期広域地震活動の推移）
ここで ETAS モデルの実際の地震時系列への適用例を示す．図 1 は 2011 年 3 月 11 日に発生した M9.0 東北沖地震前後の 1997 年 10 月 1 日から 2012 年 6 月 1 日までの期間に，北日本内陸
における下限マグニチュード M3.0 以上の地震系列に対して，東北沖地震直及び沿岸部（領域 1）前（2011 年 3 月 11 日）までを用いて ETAS モデルを推定し，その後の予測として外挿したものを図 2 に示す．
この例ではデータは 1997 年 10 月 1 日を t = 0 として，2012 年 6 月 1 日の t = 5354 日までの
3 以上の地震を含む．パラメータの推定期間（target period）[S, T ] は t = 92
マグニチュード（M）日（1998 年 1 月 1 日に相当）から t = 4906（2011 年 3 月 11 日）までの期間を設定した．t = 0 から t = 92 日までは先駆期間（precursory period）である．ここでは先行する地震から誘発された発生率が推定期間の開始時において定常性を持つ為に先駆期間を設けた．下限マグニチュード
Mc = 3.0 は推定期間においてデータの欠測が顕著に現れないレベルで設定した．データの欠損が有ると見かけ上の誘発効果が本来より低くなり，またマグニチュード列で大きめの地震や余震が続けて並ぶ．これは Gutenberg-Richter 分布（指数分布）の独立同分布からの生成列としては不自然で，珍しいことであるが，事情を知らないと b 値が低くなったり相対的に静穏化した．
りしているように見える（Ogata, 1989）
別の技術的な注意点として，α, p, c などの推定パラメータが定常地震活動としての典型的な推から大きく外れる場合はデータに多くの欠測があっ定値の範囲（例えば宇津, 1999, 10 章, p.485）
たり，地震系列そのものが非定常な為に ETAS モデルの適用が適切ではない場合がある．その他の原因として解が収束していない可能性がある．このような場合，最後に得られたパラメー
タから再度始めるか，改めて初期値として標準的な値を用いることが推奨される．

地震活動の異常性とモデリング49

図 1．東北沖地震（M9.0）前後の本州北部の地震活動．マグニチュード 3 以上の地震の震央図．
淡灰色の丸は 1997 年 10 月から 2011 年 3 月 11 日東北沖地震直前まで，黒色丸は東北沖地震直後から 2012 年 6 月までを表す．広領域 1 は北日本内陸部及び沿岸部の地震を（M7.2）の余震域，右下図は領含む．領域 2 は 2008 年 6 月 14 日岩手・宮城内陸部地震
域 3 を拡大したもので，マグニチュード 2.5 以上の震央．黒塗りの星印は 2004 年中越地震（M6.8），白塗りの星印は 2007 年中越沖地震の震央を表す．
Fig. 1. Epicenters of earthquakes of magnitude (M ) ≥ 3.0 in the Northern Honshu region, Japan, for the period from Oct. 1997 to Jun. 2012. The gray and black dots
represent the earthquakes that occurred before and after the M9.0 Tohoku-Oki
earthquake, respectively. The polygonal Region 1 covers the Tohoku inland andits western oﬀshore region. The rectangular Regions 2 and 3 include the aftershocks of the 2008 Iwate-Miyagi Prefectures Inland Earthquake of M7.2 and theswarm near Lake Inawashiro, respectively. The inset panel magniﬁes the epicenter distribution of Region 3 with M ≥ 2.5. The closed black star represents theepicenter of the 2004 Chuetsu earthquake of M6.8, and the open star represents
the 2007 Chuetsu-Oki earthquake of M6.8.

診断解析結果の図 2 から 2 カ所の顕著な相違が浮かび上がる．一つ目は 2007 年から 2008 年の期間（楕円で囲った部分）で，ここでは観測地震数が推定地震数を大きく下回る．これは 2007 年に発生した白抜き星印を震源とする中越沖地震（M6.8）の下限 M3 以上の余震が極端に少なかったことに原因すると考えられる．因みに隣接する中越地震（M6.8）は一か月で 6〜7 倍の余震数でを参照されたい．二つ目の相違は東北沖地震以降の期間（四角で囲っある．詳細は Ogata（2011）た部分）である．モデルの予測を少々上回る数の地震が発生した結果となったが，この活発化は東北地方より震源から遠方である南関東での活発化ほど顕著でない（尾形, 2012; Ishibe et al.,
．この理由については 3.2 節で議論する．この例のように発生パターンが特定時間を境と2015）して変化する場合は，次節で述べる変化点的解析が適している．

50

統計数理

第 63 巻

第 1 号 2015

図 2．図 1 の領域 1 での地震活動（M ≥ 3）に対する ETAS モデルの推定．黒曲線は地震の累積数，灰色の太い曲線は ETAS モデルによる推定累積数を表す．ただし左図では横軸を実時間に，右図では残渣時間にとったものである．時間の単位は 1 日で 1997 年 10
月 1 日から数える．ETAS の推定は，先駆期間（t = 0 日から t = 92 日）の履歴データ（最初の縦点線，1998 年初頭）から t = 4906 日（縦のもとで，推定区間として t = 92 日の点線，東北沖地震の発生時）迄のデータに合わせた．図中の楕円部と四角部については本文参照．
Fig. 2. The ETAS model ﬁtted to the events of M ≥ 3 in Region 1 of Fig. 1. The
gray cumulative curves of the estimated ETAS model are compared to those of
the observed earthquakes in black curves. The horizontal axes show the ordinarytime（left panel）in days from Oct. 1997 and the transformed time（right panel）.
The model is ﬁtted to the earthquakes in the period from t = 92 days to 4906
days（corresponding to the period from 1998 to right before the occurrence of the
Tohoku-Oki earthquake）, with the data for the events history in the preliminary
period of t = 0 to 92. The ellipses and rectangles are explained in the text.

3.

ETAS モデルの変化点解析

3.1 変化点の推定と有意性ある時点を境としてその発生パターンを記述するモデルのパラメータが前後で変化する時点（変化点）を探し，その有意性を取り扱う．ある時点を境に全部の，または一部の，パラメータを別々に推定するとモデルの適合度は向上するのかというモデル選択の問題である．推定区間[S, T ] の地震列に対する，モデルの適合度を比較する赤池情報量規準（赤池, 1976, 1981）は(3.1)

AIC = −2 max log L[S,T ] (θ|HT ) + 2k
θである．ここで，log L[S,T ] (θ|HT ) はパラメータ θ を採用するモデルの対数尤度（2.6）で，k はモデルのパラメータ数である．AIC が小さな値を得るモデルほど適合度が良いと言うことができる．
いま推定期間 [S, T ] が与えられたとき，S < T0 < T となる時刻 T0 を変化点候補としてとるとする．AIC 0 を全期間 [S, T ] で単一のパラメータをとるモデルの AIC 値とし，AIC 1 と AIC 2をそれぞれ分割期間 [S, T0 ) と [T0 , T ] において別々に推定したモデルの期間における AIC 値とすると

地震活動の異常性とモデリング

51

図 3．地震データ数 N に対する変化点パラメータの修正 AIC におけるペナルティ q(N )．右は N = 50 から 600 までを拡大したもの．
Fig. 3. The penalty function of q(N ) in the modiﬁed AIC. The right panel magniﬁes itschanges from N = 50 to 600.

AIC 0 = −2 max log L[S,T ] (θ0 ; HT ) + 2k0
θ0

(3.2)

AIC 1 = −2 max log L[S,T0 ) (θ1 ; HT0 ) + 2k1
θ1

AIC 2 = −2 max log L[T0 ,T ] (θ2 ; HT ) + 2k2
θ2の様に表現される．データから ETAS モデルの全てのパラメータ θ = (μ, K, c, α, p) を求めると
k = 5 であるが，場合によっては一部のパラメータを固定することもあるので，k0 , k1 , k2 はそれぞれで最適化で調節したパラメータの数である．
変化時点 T0 が推定期間 [S, T ] のデータからではなく別の情報で決めたとき，分割期間を合わせた 2 つのモデルの AIC の和AIC 12 = AIC 1 + AIC 2(3.3)と AIC 0 を比較して小さい方がより良い適合である．しかし，T0 を [S,T ] のデータ HT にもとづいて（3.3）が最小になるように決めるとき，T0 は未知のパラメータであり，これの最尤推定値を求めるので，これに対応した自由度としてペナルティ q を付加して
AIC 12 = AIC 1 + AIC 2 + 2q

(3.4)

と AIC 0 を比較して小さい方が良い適合である．ここで変化点パラメータの自由度の q 値はETAS モデルのパラメータと違って，1 でなく推定区間 [S, T ] の地震データの総数 N によって異なった値になる．q(N ) は解析手法的からは得ることができないが，シミュレーション計算により求めることができる（Ogata, 1992, 1999 参照）．
本稿の図 3 に q(N ) の Pade 近似式は ν = N/10 として(3.5)

q(N ) = 1 +15.325ν + 3.9376ν 2 + 0.045644ν 3
1 + 5.0900ν + 0.95595ν 2 + 0.0090963ν 3で与えられる．これは上記の原論文の誤った結果を修正した Kumazawa et al.（2010）によるものである．
3.2 応用例（地震活動の誘発・抑制効果）（図 1 の広領域 1）での地震活動を変化点解析の適用例として 2.3 節で取り上げた北日本内陸部

52

統計数理

第 63 巻

第 1 号 2015

取り上げる．この地震活動（M3.0 以上，1997 年 10 月から 2012 年 6 月まで）は東北沖地震発生
時（2011 年 3 月 11 日，14 時 46 分）を変化点とすると（3.3）の AIC が 131.5 減少し，活動パター
ンの変化はかなり有意である．ただし，この場合は当該領域のデータ外の情報から変化点を設定しているため q = 0 である．
地震活動が活発化あるいは静穏化する原因として様々な要因が考えられるが，隣接する断層
運動からの作用として発生パターンの変化を説明するのが静的応力誘発説（static stress change
である．
triggering model）
地震が起こると周辺の地殻に応力の変化が伝わる．この応力変化はその周辺部における地震活動に強く影響を与える．ある地震の影響を受けて他の地震が発生する誘発現象を議論する為の指標としてクーロンの破壊応力変化（ΔCFF，Coulomb Failure Function increment）が広く用Δτ 及び断層面を押さえる応力いられる．これは断層面をずらせる応力の変化（剪断応力変化）の変化（法線応力変化）Δσ を用いて次のように表される量である．
(3.6)

ΔCFF = Δτ − μ Δσ

μ = 0.4 は見かけの摩擦係数と呼ばれ，断層面の摩擦と断層内間隙流体圧による強度低下を含．断層運動は剪断の破壊なので Δτ が大きくなればめた便宜的な定数である（King et al., 1994）
地震が起こりやすくなる．逆に Δσ が大きいほど断層面上の摩擦が大きくなり地震発生は抑制される．これらの応力変化は本震断層運動から受け手の断層へ伝わる応力テンソルから計算される．ΔCFF が地震前と地震後でどう変化したか，その正負は本震後の周辺地域に於ける地震
活動変化の要因となる．応力変化は距離による減衰があるので遠方への影響は小さいものとなる．それでも潜在する地震数が十分多ければ統計的に意味が出て地震活動に影響する．震源断層近傍では滑りによる歪みが急激に蓄積するため影響が強く現れ，これが大量の余震およびその連鎖を引き起こす原因となる．
図 1 の領域 1 は主要な断層が東西圧縮の南北走行の逆断層で，比較的均一な地震メカニズムの地域であるので，静的応力誘発説に基づく地震活動としては地震活動が静穏化して，ETAS
モデルの活動予測を下回るはずである．しかし図 2 から明らかなように，該当地域全体で東北沖地震以降に活発化している．したがって，この増加に寄与する地震の大部分は，東西圧縮の逆断層型でない別の断層メカニズムの地震の誘発（Toda et al., 2011b），または流体が絡んだ動．事実，領域 1的誘発（dynamic triggering; 4 節参照）によるものと考えられる（Terakawa, 2013）
で増加に寄与する地震の大部分は東北地方内陸部の局所地域において多発した群発型の地震で．
ある（Kumazawa and Ogata, 2013）
ΔCFF と余震活動の変化を ETAS モデルの変化点解析により関連付けた事例研究として Ogata
et al., 2003a; Ogata, 2005, 2006, 2007, 2010a, 2010b 等が挙げられる．また，余震活動が静穏化
すると隣接部に新たな断層破壊を伴う大きな余震が起こる可能性が高いことが議論されている．このような前兆的静穏化の仕組みとして次のよ（Matsu’ura, 1986; Ogata, 1992, 2004a, 2004b）
うに考えることができる．大地震の断層内の一部または深部，隣接部において地震活動として前駆的ゆっくり滑りが起こり，これに伴う周辺部における応力の変化によって，地震活動が抑制されたり活発化されたりすると考えられる．ETAS モデルで地震活動の変化を検出することで，地震計に掛からないゆっくり滑りを，変化点ペナルティ（3.5）を考慮して検出することが可能と．
なり，GPS による地殻変動異常の観測でこれらを補うことができる（Ogata, 2005, 2007, 2011）
特に岩手・宮城内陸部地震の前駆的な地震活動については Kumazawa et al.（2010）を参照されたい．
東北沖地震前後の地震活動に話を戻す．東北沖地震の断層（国土地理院, 2011）を震源とすれば，図 1 の東北沖地震以前に地震活動が高かった領域 1 の多くの断層において ΔCFF は負にな

地震活動の異常性とモデリング53

図 4．岩手・宮城内陸部地震の余震系列（M ≥ 2.0（）図 1 領域 2）に対する ETAS モデルの推定と変化点．黒曲線は地震の累積数，灰色の太曲線は ETAS モデルによる推定累積数を表
し，時間（日）は 1997 年 10 月 1 日から数える．ETAS の推定は先駆期間として t = 0（本震）から t = 1 日までの履歴データを使って，推定区間として t = 1 日から t = 1000日
（縦の点線，東北沖地震の発生時）迄を用いた．差し込み図は東北沖地震以降の予測期間を拡大したもの．
Fig. 4. The ETAS model applied to the aftershock sequence with M ≥ 2.0, following
the 2008 Iwate-Miyagi earthquake of M7.2 in Region 2. The model is ﬁtted tothe sequence for the period from one day after the main shock（t = 1.0）to the
Tohoku-Oki earthquake（March 11, 2011; dashed line）with the preliminary datain one day period from the M7.2 main shock. The thick gray curve indicatesthe estimated ETAS cumulative function, and the thin black curve indicates theobserved counts of the aftershocks. The inset rectangle magniﬁes the cumulative
curve for the extrapolated period.

り（Toda et al., 2011a, b）
その後の地震は抑制される．たとえば図 1 の四方形小領域 2 は岩手・宮城内陸部地震の余震域だが，この余震活動は東北沖地震以降に総じて静穏化している（図 4 参の AIC は 28.5 減少するのでこの変化はかなり有照）．東北沖地震発生時を変化点とすると（3.3）
意である．
しかし，これは図 2 に示す領域 1 での大局的な活発化と相反する様に見える．先述したように活発化は図 1 の領域 3 を含め，それまで極めて低活動または無活動の地域で起こった現象である．そこで，2.3 節で議論を始めたが，図 1 の領域 3 で活発化がどのような経緯で発生したかに興味が有る．この局所地域は東北沖地震以前は殆ど地震活動が無く，東北沖地震後 1 週間
経った 2011 年 3 月 18 日から M2 を超える群発型の地震が多数発生し始めた．2012 年終わりまでのこの地震活動を解析した結果が表 1 及び図 5 である．2011 年 5 月 7 日に M4.6 の最大地震が起こったが，その発生時が変化点の最尤推定値と一致した．
（3.4）の AIC の減少量はペナル
ティ値 q = 5.3 を考慮しても 127.6 となり，地震活動はこの M4.6 地震を境に極めて有意に変化したと言うことができる．しかし変化点によって前後の区間で当てはめても，図 5 の診断解析
図から分かるように，前後どちらの区間でも適合度は満足できるものでない．しかも ETAS モ
については ETAS パラメータ値の通常の値域とはかなり異なったもの
デルの推定値（表 1 参照）である．例えば，変化点前の α 値は極めて小さな値をとっており，また μ が変化点前は変化点

統計数理

54

第 63 巻

第 1 号 2015

表 1．図 1 の領域 3 での ETAS モデルのパラメータ推定値．地震データの下限マグニチュー（a）2011 年 3 月 18 日から 2012 年 6 月までの全期間，（b）変化点前
ドは Mc = 2.5 での期間，（c）変化点後の期間．括弧内は推定誤差を表す．最後の列は各期間における推（t = 49.8 日）が最尤の変化点となり，定モデルの AIC．最大地震である M4.6 発生時
AIC の差分量は ΔAIC = (422.9 − 118.3) − 442.8 = −138.2 となる．
Table 1. The ETAS parameters for Region 3 in Fig. 1. The model is ﬁtted to the earthquakes of Mc = 2.5 and over from（a）the entire period from 18 Mar. 2011 to Jun.
2012,（b）period before the change point, and（c）period after the change point.
Their standard errors are provided in parentheses. The MLE of the change point
is t = 49.8, which coincides with the time at the M4.6 event. The improvementof the two-stage ETAS model relative to the stationary ETAS model is given by
ΔAIC = (422.9 − 118.3) − 442.8 = −138.2.

後より二桁大きい値をとっている．
群発地震活動は例えば火山地帯においてマグマや熱水等の流体の貫入と関連付けられる．これはプレート運動に由来する蓄積応力の解放というテクトニックな地震活動とは本質的に違う
非定常な地震活動である．この例でも図 5 から見て取れるように定常 ETAS モデルの理論累積関数は実際のデータの累積関数を大きく外している．これに対して変化点を増やした局所定常
ETAS モデルを用いて群発地震を含む地震活動に対応するのも一つの手段であるが，変化点の自由度は大きすぎて一般的には現実的な方法でない．そこで以下ではパラメータの連続的な時
間変化を考慮した非定常 ETAS モデルを拡張する方法を紹介する．
4.

非定常 ETAS モデル

4.1 ベイズ的平滑化
地震活動の性質はモデルのパラメータで特徴付けられるので，非定常性はパラメータの時間（2.1）のパラメータのうち，以下のにおける変化として表現される．ここでは定常 ETAS モデル
様に，常時地震活動の強度を表す μ と余震の生産性効率を表す K0 を時間依存するパラメータとする．
K0 (ti )eα(Mi −Mc )(4.1)
λθ (t|Ht ) = μ(t) +(t − ti + c)p{i: S<ti <t}

これについて安定的に推定可能なパラメータ化を考える．先ず基準となる定常 ETAS モデル（2.1）のパラメータ θ = (μ, K0 , c, α, p) を決め固定する．次に μ(t) と K0 (t) を以下のようにモデ
ル化する．
(4.2)

μ(t) = μ × qμ (t)K0 (t) = K0 × qK (t)地震活動の異常性とモデリング

55

図 5．図 1 の領域 3 での群発地震活動（M ≥ 2.5）に対する定常 ETAS モデルの推定結果．左（単位は 1 日）を
図は横軸に群発地震の発生開始時間（2011 年 3 月 18 日）からの実時間取っているが，右図は起きた順に地震の番号をとる．黒曲線は観測累積数，灰色線は定I 青）の実線部分は変化点（t = 49.8 日，縦常 ETAS モデルによる累積曲線を表す．線 （（赤）の実線部分は変化点後の期間での推点線）前の期間での推定された累積曲線，線 II定された累積曲線を表し，それぞれ外延したものを同色の破線で表す．灰色線 III は全期間で推定された累積曲線を表す．
Fig. 5. Stationary ETAS models in（2.1）ﬁtted to the datasets in diﬀerent periods from
Region 3. The left panel plots cumulative functions of the ETAS models against
the ordinary time, whereas the right panel plots those against the number orderof the earthquakes. The thin black curves indicate the observed cumulative counts
against the order of earthquakes. The blue colored cumulative curve I shows the
ETAS model ﬁtted to the earthquakes before the M4.6 event（t = 49.8 days）, thered curve II represents that of the ETAS model ﬁtted to the earthquakes after
the M4.6 event. The dashed portions of the cumulative curves represent their extrapolations. The gray curve III shows the model ﬁtted to the earthquakes fromthe entire period.

ただし

(4.3)



 
N
qμ,i+1 − qμ,i(t − ti ) + qμ,i =
qμ,i Fi (t)ti+1 − ti
i=1
i=1

 
N
N

qK,i+1 − qK,i
I(ti ,ti+1 ) (t)(t − ti ) + qK,i =qK,i Fi (t)qK (t) =ti+1 − ti
i=1
i=1qμ (t) =NI(ti ,ti+1 ) (t)

ここで，I(a,b) (t) は a < t < b で 1，それ以外で 0 の値をとる指標関数をあらわし，関数 Fi (t) は以下の形をとる．
t − ti−1
ti+1 − t(4.4)
Fi (t) =I(t ,t ) (t) +I(t ,t ) (t)ti − ti−1 i−1 i
ti+1 − ti i i+1
すなわち，qμ (t) と qK (t) は観測値に対してモデルがより良く適合するように μ と K0 を修正する為の補正量であり，各地震発生時刻 {ti ; i = 1, . . . , N } における補正値 {qμ,i } と {qK,i } をそれ
ぞれ線形に繋いだものである．本節では，この補正値の組 q = {(qμ,i , qK,i ), i = 1, 2, . . . , N } が推定されるべきパラメータとなる．このモデルの対数尤度は

統計数理

56

(4.5)

第 63 巻log L(q) =第 1 号 2015log λq (ti |Hti ) −{i; S<ti <T } TSλq (t|Ht ) dtと表記される．ただし履歴 Ht は先駆期間からの地震の履歴である．
パラメータ数が地震の数とほぼ同じであるため，これらの推定の過剰変化（overﬁtting）を防ぐ（4.6）を課し，以下の罰則付き対数尤度（penalized
為に（4.3）の関数の起伏に対して以下の罰則項
log likelihood; Good and Gaskins, 1971）
を採用する．
Q(q |wμ , wK ) = log L(q) − wμ Φμ − wK ΦK

(4.6)

ただし
(4.7)

Φμ =

N

i=0

qμ,i+1 − qμ,i
ti+1 − ti2(ti+1 − ti ) and ΦK =

N

i=0

qK,i+1 − qK,i
ti+1 − ti2(ti+1 − ti )である．ここで w = (wμ , wK ) はそれぞれの罰則関数に対する重みを表し，これが大きいほどパラメータの折れ線関数は水平線（定数）に近くなる．補正値の組 q = {(qμ,i , qK,i ), i = 1, 2, . . . , N }が条件付き強度関数（4.1）
に関して線形であるので対数尤度（4.5）は高次元に拘らずパラメータに関して単峰である（Ogata, 1978 参照）．したがって罰則付き対数尤度関数（4.7）も単峰であるので安定的に推定可能なパラメータ化である．
最適な重み w を得る為に，罰則付き対数尤度（4.7）を客観的ベイズ統計法の枠組みで取り扱う（4.6）式の第 2 項と 3 項（Akaike, 1980）．すなわち事後分布 exp{Q(q|w)} に対応する事前分布はであるが，事前分布 π(q ∗ |w) を仮定する際に注意を要する．すなわち q ∗ は q = {(qμ,i , qK,i ), i =1, 2, . . . , N } から最後の項 qμ,N および qK,N を取り除いたものである．これは縮退した分布からq ∗ に関する確率事前分布を定義するためのものである．そして，事前分布 π(q ∗ |w) を仮定したときに，q ∗ に対する事後分布は
(4.8)
Ψ(w) = L(q ∗ )π(q ∗ |w )dq ∗で与えられる．
このとき赤池ベイズ情報量基準（Akaike Bayesian Information Criterion; Akaike, 1980）は(4.9)

ABIC = (−2) max log Ψ(w) + 2 × 2
wとなる．ABIC を最小化する最適な超パラメータ w が求まったときに（4.6）を最小にする推定量は最適最大事後分布推定量（optimal maximum a posteriori estimates; OMAP estimates）と呼ぶ．
（4.6）と（4.8）を交互にそれぞれ最実際の数値計算では超パラメータ w とパラメータ q について式
大化を繰り返すことによって最適解は最終的に同時推定される．詳細は Kumazawa and Ogata（2013, 2014）を参照して頂きたい．
非定常 ETAS モデルはパラメータの折れ線の起伏への罰則（4.6）を課すため，OMAP の μ(t) とK0 (t) の時間変化は比較的滑らかな変動を持つ．しかし大地震発生時の急激な ΔCFF の変化には起こりうる．この急激な変化に対処する為に，事前に定基づく地震活動の急激な変化（jump）
jump)
まっている変化点の直前の間隔においてのみ別の重み（超パラメータ）wjump = (wμjump , wK
で小さな値に設定する処置をとる．この処置が適切かどうかは，その超パラメータ（重み）wjump
を十分小さく取って(4.10)ABIC jump = (−2) max log Ψ(w, wjump ) + 2 × 4
w

を計算して jump 無しの OMAP 解と適合度を比べることになる．しかし，（4.8）式が近似計算

地震活動の異常性とモデリング57

式（たとえば Ogata et al., 2003b 参照）で与えられているので，ABIC の定数項の厳密な計算評（4.10）の ABIC 値を直接的に比べることはできない．そこでそれぞれのモデ
価が難しく（4.9）とルの超パラメータ w を十分大きくとって μ(t) と K0 (t) の MAP を定数解にした状態の ABIC 0に関して基準化した差分(4.11)ΔABIC = ABIC − ABIC 0で計算評価の難しい定数項をキャンセルする．これを使って相対的に比べる必要がある．たと
えば，ΔABIC 値が正なら非定常 ETAS モデルより定常 ETAS モデル（2.1）が良い適合ということになる．これで他のベイズモデルも同様に基準化した差分 ΔABIC 値で較べる．通常の非ベイズのパラメトリックモデルでも定数解を基準にした差分 ΔAIC の大小比較でベイズモデルとの適合性の比較が可能である．これらの値が小さいモデルほど良い適合度を持つといえる．
4.2 応用例（群発地震活動）を図 1 の領域 3 での群発地震活動に適用してみる．領域 3 では
この非定常 ETAS モデル（4.1）
東北沖地震以前には地震活動が殆どなく，また東北沖地震後のこの地震活動の際立って異常な活動パターンのために同領域の地震時系列からは定常的な基準モデルを求めることができない．
代案として領域 3 を含む広領域 1 の東北沖地震前の地震時系列を用いる．3.2 節で述べたように
東北地方内陸部の地震は東北沖地震まで東西圧縮で南北走行の方向性を持つ比較的均質な逆断
層メカニズムで起こっており，同領域の大方の地震発生地域での共通の地震活動特性を持っているということができる．また頑健性のある広域の地震時系列のパラメータを採用することに式の補正量で適切に反映される．具体的に基準パラメータは 2より局所的な異常があれば（4.2）
節の図 2 で推定された c, α, p を固定して μ と K0 を領域 3 の地震時系列で推定し直したものを．
採用した（表 2 参照）
3 節の図 5 で求めたように，変化点を群発の最大地震である M4.6
（この地震時系列発生時から49.8 日後）に変化点を置いたモデルは，変化点を採用しないモデルより ΔABIC が 10.7 大きくなる．その最適な最大事後分布推定値（OMAP）の結果を図 6 に載せる．縦軸は対数スケールである．群発地震の発生開始より変化時点まで，黒の実線で表される μ(t) は当初大きい値をとりつつ徐々に減衰し，M4.6 の変化点で急激に減少して以降，比較的緩やかに減衰する．スパイク状に立ち上がり収束を繰り返す細い灰色曲線は全体の地震活動の条件付き強度 λ(t|Ht ) を表す．
誤差幅の線分を付けた黒点で表される K0 (t) の推定値は，変化点までは小さく，M4.6 の変化点
で急激に増大して以降緩やかに変化する．変化点前の，K0 (t) に比べて μ(t) の値が数桁大きい
点は注目すべきである．全体の条件付き強度 λ(t|Ht ) が変化点以前では μ(t) とほぼ一致している．この初期の 50 日間に於いて地震活動が非定常ポアソン過程に近い状況で発生し，M4.6 の表 2．図 1 の領域 3 の基準パラメータ．図 1 領域 1 で 1997 年 10 月から 2011 年 3 月 10 日
まで M3 以上の地震時系列から推定した c, α, p を固定し，μ と K0 のみ図 1 領域 3 のM ≥ 2.5 の地震時系列から推定した．
Table 2. The reference parameters for Region 3. The parameters c, α, p are estimatedfrom the events above the threshold magnitude Mc = 2.95, in the period from
Oct. 1997 to 10 Mar. 2011 in Region 1. Fixing those parameters, μ and K0 are
estimated from the events in Region 3 with M ≥ 2.5.

58

統計数理

第 63 巻

第 1 号 2015

図 6．図 1 の領域 3 の群発地震活動に対する非定常 ETAS モデルの推定．上段図の先端の尖った灰色細曲線は強度関数 λ(t|Ht ) の変化，黒曲線と付随する点曲線は常時地震活動度
μ(t) の変化とその標準誤差，そして小黒丸と縦線分は誘発効率係数 K0 (t) の変化とその標準誤差を表す．水平実直線と水平点線はそれぞれの基準値 μ と K0 である．滑ら（大森・宇かに減衰する灰色の太曲線は式（4.11）の第 3 項の東北沖地震による誘発効果（t = 49.8 日）を示す．横軸は群津の減衰式）を表す．縦の点線は M4.6 地震時の変化点
発地震発生開始（2011 年 3 月 18 日）からの日数である．中段図は地震震央の経度対発
生時間，下段図はマグニチュード対発生時間を表す．
Fig. 6. The non-stationary ETAS model in（4.1）ﬁtted to the datasets from Region 3.
In the top panel, gray spiky curve represents the variation of the conditional intensity rates λ(t|Ht ). The background rate μ(t) is shown in solid black, with thestandard error bounds in black dashed lines. K0 (t) is shown in black dots with thestandard error bars. The solid and dashed gray horizontal lines represent the reference μ and K0 values of the stationary ETAS model, respectively. The smoothgray curve represents the third term of the equation（4.11）by the M9.0. All of theabove estimates are plotted in a logarithmic scale. The horizontal axis represents
the elapsed times in days from March 18, 2011. The vertical dashed line shows
the change point (t = 49.8）at the M4.6 event. The middle panel displays thelongitudes versus the elapsed times of the earthquakes. The bottom panel shows
magnitudes of earthquakes versus the elapsed times.

地震活動の異常性とモデリング59

地震発生以降に地震活動は誘発地震（余震）を伴う地震系列に戻って減衰していると考えられる．
領域 3 は東北沖地震の震源域から近い地域であるので，本震からの静的応力変化も大きい．
この誘発効果が主に影響しているか否かを調べる．そこで Dieterich（1994）の摩擦構成則に基づく余震誘発の式（大森公式と同値）
を大森・宇津公式に拡張して東北沖地震からの直接的な誘発項を B 地域内での地震同士誘発を示す ETAS モデルに付け加えたモデル(4.12)λθ (t|Ht ) = μ + K0

{i: ti <t}

K1 · I{tM 9 <t} (t)K0 eα(Mi −Mc )+(t − ti + c)p(t − tM 9 + c1 )p1を導入する．第 3 項の tM 9 は東北沖地震の発生時刻，K1 , c1 , p1 が大森・宇津減衰式のパラメータを表し，この項を時刻 t > tM 9 の区間にのみ追加する．全てのパラメータは最尤法で求めた．
結果，このモデルの ΔAIC 値は上記の最適 OMAP の ΔABIC 値に 134.0 も及ばない適合度である．さらに最尤推定値によると，この第 3 項目は図 6 の滑らかな灰色曲線で表される．初期の50 日間において（4.1）
式で推定した μ(t) のサイズに遥かに及ばない変化となっている．この項は東北沖地震断層のずれから及んだ静的応力変化に基づく誘発効果というよりは，非定常 ETAS
モデル
（4.1）式で推定された μ(t) の初期の変化を説明できる別の何らかの誘発効果ということができる．
このように非定常活動度 μ(t) に支配される群発性地震の，本震から数日遅れての発生の要因としてどういった現象が考えられるのか．そもそも東西圧縮の逆断層系が支配的な東北地方では ΔCFF の減少をもたらし，全体として静穏化すべきであるが，領域 3 では逆に活発化となっは東北沖地震の地震動によって断層群内の間隙に流体た．これについて Terakawa et al.（2013）の見かけの摩擦 μ の減少をもたらして通常では起こり得ないメカニズが貫入し，これが式（3.6）
ムの地震を誘発したと議論している．図 7 のパネル（1）〜（3）は初期の約 50 日での地震時系列の時空間分布をプロットしたものであるが，震央が広がって行く様子，パネル（1）と全活動区間のパネル（4）では細部に染みいく様子が見て取れ，流体が拡散する様子と関連付けられる（Hainzl．Okada et al.（2014）では震源の移動と流体拡散の関係がand Ogata, 2005; Llenos et al., 2009）
調査されている．K0 (t) の時間変化は図 7 のパネル（4）の空間上の分布を見ることで理解することができる．K0 (t) 値の高い点と低い点がそれぞれパッチ状に分布しており，これは誘発地震（位置）に依存していると考えられる．
の生産効率 K0 が断層系
5.

まとめと今後の展望

通常の地震発生の時系列を有効に取り扱う点過程 ETAS モデルとその非定常版拡張モデルを概観してきた．定常の ETAS モデルは本震・余震の区別を取り払って地震の誘発関係を組み入れた条件付き強度関数で記述される．それは地震活動の確率予測を与える一方で，地震活動異常を検出する為の物差しとなる．より精度の高い予測を得る為には後者としての役割は無視できない．定常 ETAS モデルによる地震活動の異常変化を解析する為に変化点を挟む 2 段階の局
所定常 ETAS モデルと，ETAS モデルのパラメータが時間に依存して漸次変化する，群発地震の解析に有用な非定常 ETAS モデルを紹介した．ただし，非定常 ETAS モデルは地震系列からの推定に基づく非定常性のインバージョンを目的としているため，遡及的でありこのままでは予測に向かない．地震活動の予測として目指すべきは，この解析結果にもとづいて地震活動以外の何らかの観測時系列などの入力を含んだモデルの開発であろう（熊澤 他, 2015 参照）．
謝辞本研究は尾形良彦博士のサポートや議論が無ければ遂行することができませんでした．この場

60

統計数理

第 63 巻

第 1 号 2015

図 7．非定常 ETAS モデルの qK の空間分布．値によって色付けした qK の空間分布を表す．
（3）はそれぞれ変化点（t = 49.8）までの前期間での経度 - 緯度，時間 - 緯度，経度 （1）から
時間分布．（4）は全期間での経度 - 緯度分布．
Fig. 7. The qK value distributions on:（1）epicenters in the ﬁrst 150 days period till
the change-point,（2）hypocenters of latitude against time in the ﬁrst 150 days
period,（3）hypocenters of time against longitude in the ﬁrst 150 days period,and（4）epicenters in the whole period.

を借りて厚く御礼申し上げます．またこの原稿を査読し，改良点を指摘してくださった査読者に感謝致します．本研究では気象庁一元化震源カタログを使用しました．ここに感謝を記します．

参

考

文

献

赤池弘次
（1976）. 情報量規準 AIC とは何か, 数理科学, 14（3）, 5–11.
Akaike, H. （1980）. On the use of the predictive likelihood of a Gaussian model, Annals of theInstitute of Statistical Mathematics, 32, 311–324 （http://www.ism.ac.jp/editsec/aism/pdf/
032 3 0311.pdf）.
赤池弘次（1981）. モデルによってデータを測る, 数理科学, 19（3）, 7–10.
Dieterich, J.（1994）. A constitutive law for rate of earthquake production and its application to earthquake clustering, Journal of Geophysical Research, 99, 2601–2618.
Good, I. J. and Gaskins, R. A.（1971）. Nonparametric roughness penalties for probability densities,
Biometrika, 58, 255–277.

地震活動の異常性とモデリング

61

Guo, Z. and Ogata, Y.（1997）. Statistical relations between the parameters of aftershocks in time, spaceand magnitude, Journal Geophysical Research, 102, 2857–2873.
Hainzl, S. and Ogata, Y.（2005）. Detecting ﬂuid signals in seismicity data through statistical earthquake
modeling, Journal of Geophysical Research, 110（B5）, B05S07, doi:10.1029/2004JB003247.
Hawkes, A. G.（1970）. Spectra of some self-exciting and mutually exciting point processes, Biometrika,
58, 83–90, doi:10.1093/biomet/58.1.83.
Hawkes, A. G. and Adamopoulos, L.（1973）. Cluster models for earthquakes–regional comparisons,
Bulletin of International Statistics Institute, 45, 454–461.
Ide, S.（2013）. The proportionality between relative plate velocity and seismicity in subduction zones,
Nature Geoscience, 10.1038/ngeo1901.
Ishibe, T., Satake, K., Sakai, S., Shimazaki, K., Tsuruoka, H., Yokota, Y., Nakagawa, S. and Hirata,
N.（2015）. Correlation between Coulomb stress imparted by the 2011 Tohoku-Oki earthquakeand seismicity rate change in Kanto, Japan, Geophysical Journal International (accepted).
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, 935–953.

国土地理院（2011）. 東北地方の地殻変動，地震予知連絡会会報，86（3-34）, 184–272, http://cais.gsi.go.jp/
YOCHIREN/report/index86.html.
Kumazawa, T. and Ogata, Y.（2013）. Quantitative description of induced seismic activity before andafter the 2011 Tohoku-Oki Earthquake by nonstationary ETAS models, Journal of GeophysicalResearch, 118（12）, 6165–6182, doi:10.1002/2013JB010259.
Kumazawa, T. and Ogata, Y.（2014）. Nonstationary ETAS models for nonstandard earthquakes, Annals of Applied Statistics, 8, 1825–1852, doi:10.1214/14-AOAS759, http://projecteuclid.org/
euclid.aoas/1414091236.
Kumazawa, T., Ogata, Y. and Toda, S.（2010）. Precursory seismic anomalies and transient crustal
deformation prior to the 2008 Mw = 6.9 Iwate-Miyagi Nairiku, Japan, earthquake, Journal ofGeophysical Research, 115, B10312, doi:10.1029/2010JB007567.
熊澤貴雄, 尾形良彦, 木村一洋, 前田憲二, 小林昭夫（2015）. 伊豆東部の群発地震活動の推移の短期予測について, 第 206 回地震予知連絡会議, 2015 年 2 月 16 日, 地震予知連絡会, 九段第二合同庁舎,
東京, 地震予知連絡会会報, 94（掲載予定）.
Liptzel, R. S. and Shiryaev, A. N.（1978）. Statistics of Random Process, Part II, Springer, Berlin.
Llenos, A. L., McGuire, J. J. and Ogata, Y.（2009）. Modeling seismic swarms triggered by aseismic transients, Earth and Planetary Science Letters, 281, Issues 1-2, 59–69, doi:10.1016/
j.epsl.2009.02.011.
Matsu’ura, R. S.（1986）. Precursory quiescence and recovery of aftershock activities before some aftershocks, Bulletin of the Earthquake Research Institute, University of Tokyo, 61, 1–65.
Ogata, Y.（1978）. The asymptotic behavior of maximum likelihood estimators for stationary point
processes, Annals of the Institute of Statistical Mathematics, 30, 243–261.
Ogata, Y.（1981）. On Lewis’ simulation method for point process, IEEE Transactions on Information
Theory, IT-30, 23–31.
Ogata, Y.（1985）. Statistical models for earthquake occurrences and residual analysis for point processes, Research Memo., No.288, The Institute of Statistical Mathematics, Tokyo.
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.（1989）. Statistical model for standard seismicity and detection of anomalies by residual
analysis, Tectonophysics, 169, 159–174.
Ogata, Y.（1992）. Detection of precursory relative quiescence before great earthquakes through a statistical model, Journal of Geophysical Research, 97, 19845–19871.

62

統計数理

第 63 巻

第 1 号 2015

Ogata, Y.（1999）. Seismicity analysis through point-process modeling: A review, Pure and AppliedGeophysics, 155, 471–507.
Ogata, Y.（2004a）. Space-time model for regional seismicity and detection of crustal stress changes,
Journal of Geophysical Research, 109, B03308, doi:10.1029/2003JB002621.
Ogata, Y.（2004b）. Seismicity quiescence and activation in western Japan associated with the 1944 and1946 great earthquakes near the Nankai trough, Journal of Geophysical Research, 109（B4）,B04305.
Ogata, Y.（2005）. Detection of anomalous seismicity as a stress change sensor, Journal of GeophysicalResearch, 110, B05S06, doi:10.1029/2004JB003245.
Ogata, Y.（2006）. Fortran programs statistical analysis of seismicity — Updated version（SASeis2006）,
Computer Science Monograph, No.33, The Institute of Statistical Mathematics, Tokyo, Japan（http://www.ism.ac.jp/editsec/csm/index j.html）.
Ogata, Y.（2007）. Seismicity and geodetic anomalies in a wide preceding the Niigata-Ken-Chuetsu
earthquake of 23 October 2004, central Japan, Journal of Geophysical Research, 112, B10301,
doi:10.1029/2006JB004697.
Ogata, Y.（2010a）. Anomalies of seismic activity and transient crustal deformations preceding the2005 M7.0 earthquake west of Fukuoka, Pure and Applied Geophysics, 167, 1115–1127,
doi:10.1007/s00024-010-0096-y.
Ogata, Y.（2010b）. Space-time heterogeneity in aftershock activity, Geophysical Journal International,
181, Issue 3, 1575–1592, doi:10.1111/j.1365-246X.2010.04542.x.
Ogata, Y.（2011）. Pre-seismic anomalies in seismicity and crustal deformation: case studies of the2007 Noto Hanto earthquake of M6.9 and the 2007 Chuetsu-oki earthquake of M6.8 after the 2004 Chuetsu earthquake of M6.8, Geophysical Journal International, 186, 331–348,
doi:10.1111/j.1365-246X.2011.05033.x.
尾形良彦（2012）. 東北地方太平洋沖地震によって誘発された地域の地震活動の長期確率予測, 地震予知連（3-2）, 92–99, http://cais.gsi.go.jp/YOCHIREN/report/kaihou88/03 02.pdf.
絡会会報, 88
Ogata, Y., Jones, L. M. and Toda, S.（2003a）. When and where the aftershock activity was depressed:
Contrasting decay patterns of the proximate large earthquakes in southern California, Journalof Geophysical Research, 108（B6）, 2318, doi:10.1029/2002JB002009（1-12）.
Ogata, Y., Katsura, K. and Tanemura, M.（2003b）. Modelling heterogeneous space-time occurrencesof earthquakes and its residual analysis, Applied Statistics（JRSSC）, 52, Part4, 499–509.
Okada, T., Matsuzawa, T., Umino, N., Takahashi, H., Yamada, T., Kosuga, M., Takeda, T., Kato,A., Igarashi, T., Obara, K., Sakai, S., Saiga, A., Iidaka, T., Iwasaki, T., Hirata, N., Tsumura,N., Yamanaka, Y., Terakawa, T., Nakamichi, H., Okuda, T., Horikawa, S., Katao, H., Miura,T., Kubo, A., Matsushima, T., Goto, K. and Miyamachi, H.（2014）. Hypocenter migrationand crustal seismic velocity distribution observed for the inland earthquake swarms induced by
the 2011 Tohoku-Oki earthquake in NE Japan: Implications for crustal ﬂuid distribution andcrustal permeability, Geofluids, 15（1-2）, 293–309, doi:10.1111/gﬂ.12112.
Papangelou, F.（1972）. Integrability of expected increments of point processes and related randomchange of scale, Transactions of the American Mathematical Society, 165, 483–506.
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.
Toda, S., Lian, L. and Ross, S.（2011a）. Using the 2011 M=9.0 Tohoku earthquake to test the Coulomb
stress triggering hypothesis and to calculate faults brought closer to failure, Earth Planets andSpace, 63, 725–730.
Toda, S., Stein, R. S. and Jian, L.（2011b）. Widespread seismicity excitation throughout central Japan

地震活動の異常性とモデリング63

following the 2011 M = 9.0 Tohoku earthquake and its interpretation in terms of Coulomb
stress transfer, Geophysical Research Letters, 38, L00G03, doi:10.1029/2011GL047834.
Utsu, T.（1961）. Statistical study on the occurrence of aftershocks, Geophysical Magazine, 30, 521–605.
宇津徳治（1999）.『地震活動総説』, 東京大学出版会, 東京.
Utsu, T., Ogata, Y. and Matsu’ura, R. S.（1995）. The centenary of the Omori formula for a decay lawof aftershock activity, Journal of Physics of the Earth, 43（1）, 1–33.
Vere-Jones, D.（1970）. Stochastic models for earthquake occurrences（with discussion）, Journal of theRoyal Statistical Society, B32, 1–62.
Vere-Jones, D. and Davies, R. B.（1967）. A statistical survey of earthquakes in the main seismic region of
New Zealand, Part 2—Time series analyses, New Zealand Journal of Geology and Geophysics,
9, 251–284, doi:10.1080/00288306.1966.10422815.

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

Modeling Seismicity Anomalies
Takao Kumazawa
The Institute of Statistical Mathematics

The epidemic-type aftershock sequences（ETAS）model is useful for statistical analysis of time sequences of earthquake occurrences. The conditional intensity function of thepoint process model is defined for the occurrence rate of an earthquake in the immediate
future, and is indispensable for probability forecasting of earthquakes in a given future
period. Also, it provides a standard seismicity model that is useful for detecting anomalies
from a series of earthquake occurrences. Then remodeling of the intensity function can
improve the earthquake forecasting. In this manuscript, application methods of the stationary ETAS model are explained in detail for analyzing seismicity changes. Then thismodel is proposed for inversion of abnormal seismicity. This is applied to a dataset ofswarm activity triggered by the 2011 Tohoku-Oki mega-earthquake.

Key words: Point process, conditional intensity function, stationary ETAS model, anomaly in seismicity, non-stationary ETAS model.