統計数理（2019）
第 67 巻 第 2 号 215–228©2019 統計数理研究所

創立 75 周年記念号［研究ノート］

統計地震学の発展と地震活動予測：
個人的経験と展望
尾形 良彦†（受付 2019 年 1 月 4 日；改訂 3 月 15 日；採択 3 月 26 日）

要旨

地震予測の観点から統計地震学の発展と筆者の研究経験を説明する．点過程モデルによる地震活動の予測と統計的診断解析および地震活動の物理現象との関わりに焦点を当てる．
キーワード：地震活動，点過程モデル，ETAS モデル，階層ベイズ法，確率予測，前
震の識別．
1.

はじめに

筆者の統計科学や地震学の研究者に向けた研究紹介については既に尾形（1993）や Ogata（2013, 2017a）などがある．なので，同様なことを繰り返すのは気が引けたが，本稿に似たタイに倣って筆者も先人の成果を受け継いだ経緯の個人的経験と見トルの論文（Vere-Jones, 2006）
解を述べたい．
筆者の折々の研究活動は，統計数理研究所学術研究レポジトリ https://ismrepo.ism.ac.jp/（最終閲覧日 2019 年 4 月 5 日）で検索でき，とくに統計地震学研究プロジェクトの成果と当時の目標をまとめた外部評価報告書（英文和訳付き 統計数理研究所編 https://www.ism.ac.jp/
も参照できる．
evaluation/index_j/toukei_jisingaku2006.pdf，最終閲覧日 2019 年 4 月 5 日）
本稿では，主に地震の短期・中期予測に関わる点過程解析について解説するが，予測全般については，筆者が編集した地震予測の特集号（統計数理 63 巻 2015 年発行; 統計数理研究所編，
http://www.ism.ac.jp/editsec/toukei/tokeisuri-63j.html，最終閲覧日 2019 年 4 月 16 日）がある．
と略記し，これに論文著者名を付加して引用する．
これを本稿では「統計数理 63」
本稿では割愛するが，長きにわたり膨大化した地震発生様式のデータベースである震源カタ
ログは，観測システムの発展に伴い，時空間的に不均質である．データの欠測推定や系統誤差の補正など，本質的な情報を抽出するために，ベイズモデルと解析法を考慮し応用した．これ
らの概要については Ogata（2017b, Sections 3.5, 3.6）を参照頂きたい．
本稿で述べる内容に関して全て引用するとなると，与えられたページ数を大幅に超えるので本稿では，以上の解説論文に含まれていない引用文献を記すに止め，その他は上記の解説論
文の文献を通して検索されたい．特に，上記の統計学会誌 60 周年号解説論文，尾形（1993）は本稿で割愛した内容を頻繁に引用し，Supplemental Material となるので文献欄に論文の PDF
address を付けた．

†

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

216

2.

統計数理

第 67 巻 第 2 号 2019

統計地震学と地震活動の研究

2.1 地震活動
筆者が，点過程の応用問題を模索し，地震学会に出没して特に興味をもったのは地震活動の分科会であった．当時急速に蓄積された地震データの解析の報告が，相次いだものであった．
（例えば宇津, 1999, 参照）
データに基づく地震の記述と理解のための研究分野としての統計地震学は，明治時代以来の地震学の主分野であった．一回の地震について分かることは「発生時刻，規模，場所」だけであったが，これら 3 つの性質を相手にして，地震発生に関する様々な性質を明らかにしようとした多くの研究がある．後に現代地震学の基礎を築いた安芸（1956）の総合論文から分かる様に，昔からこの分野で議論されていたのは，地震発生の周期性，気圧や海洋潮汐変化など地球物理データとの因果性の研究，震源の移動，地震活動パターン分類，地震の規模（マグニチュード）の分布などであるが，それらの解析結果の妥当性をめぐる最大の争点は，余震や群発地震のような地震の続発性による，統計的検定や仮説の再現性の難しさにあった．
したがってアメリカの地震学界では余震の統計的研究は無視された．地震の規模（マグニ
チュード）を初めて考案したリヒター（Charles Richter）などは，余震は何の役にも立たないデー
タのゴミ
（ノイズ）
と言っていた．地震の予測にはゴミである余震を取り除いた本震の活動を「地震活動見なければいけない，という考えだったようである．この思想は 20 世紀後半でのの静穏化」，「地震発生の季節性や地球潮汐との同期」
や「異常現象の大地震との因果性」などに関する統計的研究に受け継がれる．すなわち，地震活動のデータから
「余震」
を取り除く
「除群
を使った研究である（尾形, 1993, 7.1 節参照）．先ず除群データの帰無仮説（定（de-cluster）法」
常ポアソン過程）を否定してから，焦眉の科学的主張をする．
しかし除群法は様々で，異なった結果を導いたりするだけでなく，余震活動自体の研究に意
味をなさない．事実，地震活動は，時間的には逆ベキ型減衰の自己相関をもつ長記憶性，空間的にはフラクタル性をもつ自己相似性（尾形, 1993, 5.1, 6.5, 8.1 節; Ogata and Katsura, 1991;
がみとめられる．このため，除群は首尾一貫したものでないといえGuo and Ogata, 1997 参照）
る．とくにマグニチュード 7 以上の大地震といえども，履歴依存性や長記憶性がある（Ogata．しかるに，大地震は定常ポアソン過程になるはずだ，という全and Abe, 1991; Ogata, 2017b）
く誤った帰無仮説による検定をし，地震波の原記録に基づいた阿部カタログのマグニチュード
を雑に改竄した，20 世紀の地震カタログ（尾形, 1993, 5.2 節; Engdahl and Villasenor, 2002）が流布し，世界の被害地震の基礎データとして使われているのは深刻である．
2.2 余震の研究と統計地震学の発展
日本では，余震データを重要な情報として研究することで，地震活動の本質的な性質に迫ることになった．大きな地震の後に，より小さな地震が数多く発生するが，古来このような地震
群を余震と呼び，そのきっかけとなった地震は本震と呼ばれている．
余震の発生は，本震直後は極めて頻繁で，その後，時間とともに徐々に低下するが，なかな
か元の状態に戻らない．このような余震減衰の定量的な関係を初めて論じたのは大森房吉である．1894 年の濃尾地震の余震数の頻度の時間経過について，減衰のしかたを調べて「物理現象の減衰だから当然指数関数だろうと考えて当てはめてみたが良く合わない，然るに双曲線だとよく適合する」と述べている．さらに 1957 年，宇津徳治は日本や世界の余震を調べ, 単位時間あたりの余震頻度の減衰率が(2.1)

ν(t) = K(t + c)−p

統計地震学の発展と地震活動予測：個人的経験と展望217の形になることを示した．ここで t は，本震の発生時刻からの経過時間である．宇津は，余震の頻度 ν(t) と経過時間 t を両対数方眼紙にプロットし，その減衰が漸近的に直線上に乗ることを示し，直線の傾きから指数 p の推定値を得た．そして，濃尾地震の余震は現在に至るまで，式の減衰が続いていることを示した．大森の言う上記の双曲線とは p = 1
一世紀以上でも（2.1）のことである．一方，係数 c の解釈については議論がある（Utsu et al., 1995 参照）．宇津は「改良大森公式」と名付けたが，今日筆者らは「大森・宇津の公式」と呼んでいる．
（2.1）式を
さらに宇津は，多くの余震減衰のプロットで，単一の（2.1）式で説明できない場合があることも示している．大きな余震の後，飛躍的な増加が起こって，再び減衰する二次余震（余震の余震）を発見し，その減衰も大森・宇津の式（2.1）に従っていることを突き止めている．これは，余震が本震のみに誘発されるという当時の常識を覆したものである．このほか宇津は，余震群の空間的広がりと本震の規模（マグニチュード）との関係の定量的研究や，本震と最大余震のマ
グニチュード差の経験分布など，余震の多様な性質に関して多くの事例の統計的解析を行った．
さて，ここまでの研究では，余震の単位時間あたりの数として分析されてきたが，筆者は，（発生強度）と捉え，1 つずつの地震（2.1）式を，極短い時間内に 1 つの地震が起こる確率の微分の発生時刻の記録 {ti ; i = 1, 2, . . .} をそのままデータとして使う，発生強度による推定法を提（さい案した．これは余震列を非定常ポアソン過程として（2.1）式のパラメータ K, c, p を最尤法ゆうほう）で推定する方法である．本震直後の欠測の問題や二次余震を含む場合も，赤池情報．今日，大地震が起き量規準（AIC）でモデル比較の上，偏りの無い推定ができる（Ogata, 1983）式を推定して，余震の予報確率が計算されることになっると，余震発生時刻データを基に（2.1）ている．背景として 1960 年代ごろからの，確率論や統計学の分野で発展した点過程とその統計的解析法の理論的整備が，これらを可能にしたのである．このことは次節で詳しく述べる．
同じく 1960 年代ごろからの地震の物理の飛躍的前進として，地震発生の応力蓄積の根源となるプレートテクトニクス理論が展開され，「地震の揺れの原因は断層面の急激な食い違い運動である」という弾性反発説が定着した．大地震ごとに断層運動モデルが，地震波の解析によって求められる．今世紀では，それは衛星観測の測地的変化から詳しく求められて，地震の発生
様式や地震力学の物理的解明が進んだ．かくて地震カタログは，破壊の始まった時刻と位置，食い違った断層の面積とすべり量を表すデータと理解された．さらに断層面および直交面の方向や断層の重心などのデータも追加されるようになった．今日の地震活動は，これらのデータに基づいて研究されている．
しかし，小さい地震が増々捉えられ，データベースが豊富になるにしたがって，地震の発生
様式の地域性や複雑系が顕著化して，詳細な統計的把握や予測も難しくなってくる．本質的で確かな情報を抽出するために，時間的空間的に非定常または不均質な予測モデルを考慮する必要があり，これらの統計モデルを用いた研究が避けられない様になってきた．時空間モデルなどの逆問題にかかわる大量の未知パラメータを含む大規模モデルを取り扱う階層的なベイズ法の助けが必要であり，この様な手法抜きでは地震予測の研究自体が難しくなっている．
3.

点過程の統計解析と統計地震学の研究の経過

3.1 点過程の地震統計解析の始まり
点過程は，事象（出来事）の発生時間と付置のマーク（ジャンプなどのスカラー量または震源の空間位置など）の時系列の数学的表現として，確率論的に深く研究されてきた．通常の連続確率過程の基礎が正規過程のホワイトノイズ（white noise）であるのに対して，点過程の基礎は
定常ポアソン過程である．各種の応用に対応して点過程の一般化が進んで，1970 年代は点過程の確率論，統計解析，応用の統一的な研究が進んだ時代であった（Lewis ed., 1972 参照）．

218

統計数理

第 67 巻 第 2 号 2019

ヴェア・ジョーンズのトリガーモデル（Trigger model; 詳しくは尾形, 1993, 6.1, 6.4, 6.5 節参照）
は，地震活動の一般化点過程モデルの先駆けである．本震は定常ポアソン過程に従い，余
震は本震のみに誘発され，大森・宇津の減衰関数（2.1）の非定常ポアソン過程で，地震活動は
これらを重ね合わせたものである．これは当時の地震学的常識を反映したものであるが，デー
タのどれが本震であるか余震であるかの同定は明確でない．それらを設定するための組み合
わせ論的複雑さを伴い，尤度（ゆうど）
関数で計算するのが極めて困難である．したがって，
モーメント法に基づいた推論や解析が使われた．特に，ホークス（A. G. Hawkes）と彼の学生
Adamopoulos（1976）
は，トリガーモデルとホークスモデルを二次モーメントのスペクトル尤
で計算したが，精度は良くない．他方，本来のトリガーモデルではない度（尾形, 1993, A3 節）
が，筆者は発震時刻のみのデータから各地震の誘発地震数（クラスターサイズ）などを推定した．
（Ogata, 2001a）
3.2 統計数理研究所における点過程の光臨
1970 年代当時，赤池弘次は，多変量時系列の統計的同定によってフィードバック効果を考慮した予測と制御の実用化を目指し，外力変数を含む自己回帰モデルを状態空間表示し，最小二（FPE）や赤池情報量規
乗法や最尤法による最良予測のモデル選択のため Final Prediction Error
準（AIC）を提案し，世界の統計学界を牽引していた．周囲の若手研究者に対しては，新しい分
野での各種統計モデルの構築を促していた．
ニュージーランドのヴェア・ジョーンズ（David Vere-Jones）は，点過程の理論とモデルで地震データと取り組んだ先駆者であるが，宇津を初めとする日本の地震学者との交流を求め，1976
年，赤池に招待されて統計数理研究所を数か月間訪問した．その間，点過程に関する一連の特別講義をした．このときの経緯が Vere-Jones（2006）や Ogata（2018）によって述べられている．
彼の講義を聞いて直ちに赤池は，時系列解析と同様な予測の展開を，点過程モデルの予測に関する中心概念である「条件付き強度関数」(3.1)

λ(t|Ht ) = P {an event occurs in [t, t + Δ)|Ht }Δ + o(Δ)

に期待した．これは事象（点）が発生する条件付き確率の微分量である．ここで条件 Ht は時刻
t までの，事象発生の履歴および関連情報のことである．単純な点過程なら Ht は時刻 t 直前までの発生時刻の系列 {t1 , . . . , tn } である．条件付き強度関数をモデル化することによって，パラメータの推定や事象時刻のシミュレーションなどで, 予測問題を考えることができる．例えば，条件付き強度関数で，最後の事象の時刻 tn から次の事象が起こるまでの時間の確率分布
F (t − tn |Htn ) 及びその密度関数 f (t − tn |Htn ) との関係
（Hazard 方程式）が得られ(3.2)

λ(t|Ht ) = f (t − tn |Htn )/{1 − F (t − tn |Htn )}となる．これを解くと隣り合う事象の時間間隔の分布が得られる．
3.3 点過程の尤度計算とシミュレーション法
特に，ヴェア・ジョーンズの講義でモデルの具体例として紹介された，点過程の「自己回帰モデル」ともいうべき，ホークス（Hawkes, 1971）の自己励起過程（self-exciting process）が注目された．このモデルの条件付き強度関数(3.3)

λ(t|Ht ) = μ +


ti <t

 t

g(t − ti ) = μ +g(t − s)dNs0

が，事象の過去の発生時間の線形回帰形式になっているからである．時系列研究グループ（ヴェ統計地震学の発展と地震活動予測：個人的経験と展望

219

ア・ジョーンズは Akaike School と呼んでいた）が，その後直ちに取り組んだのは，ホークス点g(·) を AIC と最尤法によって推定することであった．
過程の応答関数（response function）
情報理論の分野では，条件付き強度関数とその尤度関数が理論的に導かれたばかりで，これを数値的に求めることは当時挑戦するに足る課題であった．ヴェ
によって最尤推定値（MLE）
ア・ジョーンズの助言で，赤池の第 5 研究部の最初の部下であった尾崎統（Ozaki, 1979）が，実式でシミュレーションしたホークス過程の標本データ
行したのは g(·) を指数関数とし，（3.3）
から，対数尤度の最大化を Davidon-Fletcher-Powell 法のような効率的な準ニュートン法によって，計算し係数の MLE を得たことである．
このことに触発されて筆者が取り組んだのは 3 課題あった．以下の 2 課題は学術誌に掲載
後，合わせて学位論文として提出した．もう一つの課題は次節で述べる．
（i）上述のような，ホークス型モデルを含む，条件付き強度関数で特徴付けられた点過程のMLE や最大尤度の漸近理論．すなわち MLE の収束や誤差分布，尤度比統計量の漸近分布をマ
ルチンゲールのエルゴード性や中心極限定理によって証明した．
（ii）条件付き強度関数を直裁的に使った点過程のシミュレーション法の提案．尾崎が実行したシミュレーション法は Hazard 方程式（3.2）をニュートン法で反復して解くものであるが，多変量の場合など条件付き強度関数が複雑になると解が安定して求まらない．そこで筆者はLewis and Shedler（1979）の間引き法（thinning method）に着目した．これは非定常ポアソン過程のシミュレーション法で，任意の密度関数から標本を生成する従来の rejection sampling 法と本質的に同じ方法である．筆者は間引き法を一般化し，多変量（マルチ・チャンネル）などにも拡張した．この正当性をマルチンゲール理論で証明し，多数の点過程モデルのシミュレー
ションの実例を示し，尤度比検定で正確さを実証した（Ogata, 1981）．この方法は現在，多分野で採用されている．
3.4 因果関係分析のための点過程モデリング
3 つ目の課題は，除群法に代わるものとして，地震の続発効果を含む，以下の解析モデルの開発であった．すなわち，トレンド，季節性などの周期的な要素の解析，および他の外的デー
タからの誘発効果の存否を確かめる因果関係の解析である．とくに因果関係の量的解析は，多くの科学分野でますます必要となっている．従来から一応，点過程系列間の相互相関統計量を使用することができるが，多くの続発事象を伴う場（cross Palm intensity; 尾形, 1993, 3.3 節）合や非定常性を含む場合から因果関係を定量的に導くことは困難である．
したがって課題はホークス型モデルを一般化して，トレンド，周期性，および因果関係を分析するモデル(3.4)

λ(t|Ht , Ft ) = μ + f (t) + C(t; T0 ) +



ti <t

g(t − ti |Mi ) +h(t − uj )ξ(uj )uj <tを開発し実用化することだった．ここで Ht は解析対象の点過程の発生時間の履歴で，Ft は外部入力データの時刻 uj とマーク ξ(u) を表す．時刻 uj は等間隔でも不規則間隔でも良い．これで，日本や世界各地の地震活動の季節性や，地域間の地震活動の因果関係を導くことができた（尾形, 1993, 3.1, 3.3 節, 統計数理 63 尾形論文 8, 9 節参照）．（3.4）式右辺の第 4 項は，後ほど次節で述べる ETAS モデルを含む形にもなっている（例えば Kumazawa et al., 2016）．
このアプローチは AIC を採用し，誘発候補データの地震活動への因果性を調べる方法の原型である．今日まで，様々な地震活動の異常に基づいた大地震の警戒（alarm）情報が，電子メールなどで定期的に報知されているが，これらの有効性は議論の余地がある（Jordan et al., 2011）．
．他
そのような警戒情報には，客観的な有意性の評価が可能である（統計数理 63 庄・尾形論文）

220

統計数理

第 67 巻 第 2 号 2019

方，警戒情報が確率予測として採用されるには，対象地域の標準的地震活動度と比べた優位性を示すために，警戒情報との因果関係の有意性および確率利得が，提供されるべきである．ここで「確率利得」とは「大地震の確率が平常時に比べて何倍高くなるのか」という量である．因は，異常なデータや大きな地震との間の前兆性の統計的関係を調べるために果関係モデル（3.4）
使用することができる（統計数理 63 尾形論文 8 節参照）．
3.5 ETAS モデル
このモデルが考案される頃，統計数理研究所は大学共同利用研究所として改組された．筆者
は「数理地震学」の共同研究会や複雑系・フラクタルの学習会を主催し，自己相似性やフラクタ
ル次元（点群の集中度）の最尤法など，地震活動の特徴付けの研究（尾形, 1993, 5.1, 8.1 節; Ogataを行ったが，それらの単なる解釈に終始せず，震源データから直裁的に地and Katsura, 1991）
震活動を予測する実用的な点過程モデルを模索した．手掛りの多くは前節で述べたような，宇津によって探求・確立された余震の経験則であり，これらに基づいて，一定のマグニチュード
Mc 以上の地震発生を予測するために epidemic-type aftershock sequence（ETAS）モデル(3.5)

λ(t|Ht ) = μ +ν(t − ti )eα(Mi −Mc )ti <tを求めた．ここで Ht は過去の地震発生時とマグニチュードの系列の履歴データである．ν(t)は大森・宇津の式（2.1）
で，それに掛け算しているマグニチュードの指数関数の大小によって，大きい地震には多くの余震が誘発され，小さい地震もそれなりに余震を誘発する．時刻 t での条件付き強度は，線形の重ね合わせである．パラメータ α 値はマグニチュード効率で，大きければ典型的な本震・余震型の地震系列になり，小さければ群発地震のようになる．パラメータμ は常時地震活動率と呼び，その地域特有の地震発生の強さを意味する．ただ，常時活動の一
部は昔の地震の余震活動かもしれない．この様に，規模，常時発生率，余震の起き方など，そ（μ, K, c, α, p）は最尤法
れぞれの地域の特徴が定量化される．ETAS モデルの 5 つのパラメータで求められる．ETAS モデルは，シミュレーションで各地の標準的な地震活動を予測するのに良く使われている．
の一定の期間の積分で，平均の地震数が与えられる．このことに基
条件付き強度関数（3.5）
づいて，ある地域において或る期間に地震の数が ETAS の予測に対して有意に少なくなる現象（尾形, 1993, 7.2, 7.3 節; 実例図に注目した．これを「相対的静穏化（relative quiescence）」というとしては http://www.ism.ac.jp/editsec/toukei/pdf/45-1-139.pdf 参照; 最終閲覧日 2019 年 4 月
5 日）．ある地域で地震活動が静穏化しているように見えても，それは単に以前の地震の余震活
動の終息を意味するものかもしれない．ある程度の活動があっても，本来は以前の地震の余震
活動がもっと盛んであるべきなのに，それに比べ，何らかの理由で不活発なのかもしれない．
ETAS モデルを採用すると，このような効果を計算に入れて，減衰する余震活動でも，期待される活動レべルを基準に相対的な静穏化を論じることができる．この点で，従来の除群に基づく静穏化と違う概念である．同様に「相対的活発化」
が定義できる（Ogata, 2005, 2007, 2011a;
Kumazawa and Ogata, 2013; Kumazawa et al., 2010）．このような診断解析的な応用は多数例あり，下記 4.1 節で述べる．
3.6 階層時空間 ETAS モデル「時空間 ETAS モデル」は ETAS モデル（3.5）を一般化して，地震の位置（震央）の変数を考慮したモデルで，リアルタイムで場所も地震の群の広がりも予測する（Ogata, 1998）．これによると今どこで危険度が高く，どこの場所が，これまでどの程度活発だったか一目瞭然である．こ統計地震学の発展と地震活動予測：個人的経験と展望

221

うしたモデルの様々な変型版が「地震予知可能性の研究のための共同研究」（Collaboratory for
the Study of Earthquake Predictability; CSEP; http://www.cseptesting.org/; 最終閲覧日 2019
年 4 月 5 日）の各国のテストセンターに於いて，現在まで 10 年以上予測実験が継続中である．
さらに，応力変化や発震機構など地殻内情報による追加情報をもとに時空間 ETAS モデルより良い適合度を目指す余震メカニズムの研究が盛んに追求されているが，今のところリアルタイ
ム予測に結びついていない．
ETAS（HIST-ETAS）モデルは，地域的特徴を詳細に求
階層的時空間（hierarchical space-time）
めるために，パラメータも位置に依存させるものである（Ogata et al., 2003; Ogata, 2011b）．
これらは隣接する地震を結んだドローネ 3 角形分割上の部分的線形関数で表現され，特に地震集中域では高精度な変化が表現できる．しかし HIST-ETAS モデルの総パラメータ係数は，
データサイズの数倍となる．そこで，経験的ベイズ法と赤池ベイズ情報量規準を用いて，パラメータ係数間の制約の最適な事前分布を求める必要がある．
HIST-ETAS モデルの常時地震活動率 μ(x, y) の高い地域は，GPS で求めたストレス蓄積率
変化の大きい地域と良く対応し，歴史的被害地震の場所もよく説明している．また，今世紀になってからの大地震の発生場所をよく言い当てており，将来の大きな地震の場所の確率予測の基礎となりうる（統計数理 63 尾形論文; 尾形, 2008, 2017）．
4.

点過程の診断解析と地震の物理

4.1 相対的静穏化と後続の大地震との関係
近傍で大地震が連発するのは珍しくない．日本に於ける各々の大地震について，その後に起きた全ての大地震との時間差と距離を重ね合わせプロットすると，直後に近傍で別の大地震が起きる頻度は平時より数倍高くなり，離れるにしたがって逆比例することが分かる（Ogata,．これには理由がある．断層が急激にずれるので周辺の応力が増加し，近傍の断層がず
2017b）
れ易くなるからである（統計数理 63 尾形 5 節参照）．
また著者は日本に於ける過去の大地震の余震列を 76 例調べた（Ogata, 2001b; 詳細記録はOgata, 2001; http://bemlar.ism.ac.jp/ogata/JGR01supplement/ 参照, 最終閲覧日 2019 年 4 月
9 日）．そのうち，相対的静穏化は 45% ほどの余震列で認められ特に珍しくないが，近辺で別の大きな地震が起こる確率は，前述の場合より更に 3〜4 倍高くなる（Ogata, 2017b）．相対的「ゆっ
静穏化の原因の一つとして，余震域で破壊応力の低下（stress shadow）をもたらすようなが近傍の断層内で起きる場合が考えられる．その際，そのような slow
くりすべり（slow slip）」
slip がきっかけとなって，新たな断層の破壊が誘発される確率があると考える．この様な先行
slow slip を数例，GPS 測地データの異常変化の解析で確認することができた（Kumazawa etが，更なる事例を数多く解析し，リアルタイムでモニター
al., 2010; Ogata, 2007, 2010b, 2011a）
できる様になることが望まれる．
1990 年頃から現在まで，大地震と余震活動や広域地震活動の前駆的な，相対的静穏化や相対的活性化について，筆者は地震予知連絡会で報告した（地震予知連絡会報，国土地理院;
http://cais.gsi.go.jp/YOCHIREN/report.html; 最終閲覧日 2019 年 4 月 5 日; 報告内容は「統計数理研究所」
として検索）．さらにそれらの大半の数十に亘る論文を地震関連学術誌（統計数理
63 尾形 7 節参照）に掲載したが，それらは一例を除いて事後報告である．したがって，今後ETAS で相対的静穏化・活発化現象や，slow slip を地震の予測に繋げる為には，地震活動や地殻変動の異常検出を事前に探索し，効率的に検出できる統計的モデルや方法の研究を勧めたい．そのヒントの一部を以下に記す．

統計数理

222

第 67 巻 第 2 号 2019

4.2 群発地震
プレート境界における群発地震と slow slip との関係については，数例の報告がある（たとえば Nishikawa and Ide, 2018）．従来から，群発的な地震活動は多様であることが指摘されている（宇津, 1970）．小さい α 値の ETAS モデルで表現できる群発地震（第 2 種群発地震）もある（3.5）式が，火山性の群発地震のように ETAS モデルで良く表現できないものもある．これはの μ のパラメータが時間とともに変化する非定常 ETAS モデルが良い適合を示す（Kumazawa．パラメータ μ の時間変化は，マグet al., 2010, 2016, 2017; Kumazawa and Ogata, 2013, 2014）
マや熱水などの流体の断層系貫入に起因する固着断層の弱化と地震発生率との間に，量的な物
理的関係を提供する可能性がある（統計数理 63 熊沢論文 4 節，尾形論文 9 節参照）．
4.3 余震の時空間診断解析の積分曲線で時間を伸縮変換し，余震群内の時空間的分布を詳
余震活動の大森・宇津式（2.1）
しく見ると，余震域全域で一様な経緯を示す場合もあるが，各部分領域で異なったパターンのこともある（Ogata, 2010a）．特に，大きな余震が発生する前には局所の静穏化パターンを示すことが多い．このようなパターン異常は局所的な slow slip などが余震域の中でも起きている可能性がある．
5.

前震

5.1 マグニチュード系列について
地震のマグニチュード系列に関する研究は未だに発展途上である．伝統的な標準モデル（Gutenberg-Richter 則; Gutenberg and Richter, 1944）は，一定のマグニチュード Mc 以上の指
数分布
(5.1)

F (M |b) = 1 − 10−b(M −Mc ) = 1 − e−β(M −Mc ) ; M ≥ Mc , β = b log10 eに則り，独立同分布である．係数 b は最尤法で決めると Mi − Mc の標本平均の逆数に比例したものとなり，日本地域基準値では凡そ b = 0.9 である．ただし，気象庁などの多くの地震カ
タログの場合，マグニチュード値は四捨五入して 0.1 刻みで与えられているので，不偏推定の為には，下限のマグニチュードは 0.05 を差し引いたものにする．CSEP の予測実験に提出されを仮定している．
ているモデルは殆んど（5.1）
しかし，マグニチュードの時系列は，地震の発生パターンなどの履歴データに依存している可能性がある．また，マグニチュード分布の裾の部分を高めたものも考えられる．例えばアメの第三次カルフォルニア地震予測計画 UCERF3- ETAS（Field
リカ合衆国地質調査所（USGS）
では，サンアンドレアス断層などの近傍活断層への誘発可能性を含む，固有地震型
et al., 2017）
マグニチュード分布を提案している．
様々な観測異常現象と大地震との因果関係をマグニチュード予測に組み込むことが考えられる．たとえば，次節で議論する「前震」の事前識別によって有効な予測を構成できる（Ogata,2017b; Ogata et al., 2018）．その予測の優劣は標準 Gutenberg-Richter 分布モデル（5.1）の予測結果に比べた対数尤度比基準（情報利得）
で検証できる．このような研究は現時点で他に未だ余り見られない．
5.2 前震の統計的特徴
本節では，先ず本震を基準にした前震の相対的な研究（尾形, 1993, 8.3 節）を紹介する．これ
まで様々な前震の定義があるが，次節で予測を考慮するために，以下の定義で考える．

統計地震学の発展と地震活動予測：個人的経験と展望

223

全ての地震
（M ≥ 4.0）に関して，或る時空間的な距離より近いものを全て繋いで，孤立した地
震や地震の群れに分離する（Single-link 法）．各群れの中の最大の地震を本震と呼ぶ．本震と，群内で先行する最大の地震とのマグニチュード差が，0.45 以上の場合，先行する地震を「前震」と呼ぶ．このマグニチュードの差は常識的に，なるべく大きくしたいが，そうすると全体の地震群に占める前震群の割合が少なくなるので，統計的な議論のために便宜上決めた．
余震に比べて前震は希少で，事例ごとに起こり方の相違が著しいので，全ての前震の時刻や位置座標を，本震が原点になるように移動して重ね合わせた統計的な特徴を述べると以下のようになる．（i）その時間頻度は，過去に向かって大森・宇津の（2.1）式に則っている．（ii）本震と（ii）は言い換えると本震に集中することを示すの位置から離れるとベキ則で頻度が減る．
（i）が，それぞれの座標軸での周辺分布が一様になるように座標変換をすると，（iii）あたかもドー
ナツ型が本震に向かって収束する様相になる．これは前震の本震へのマイグレーションを意味．
する（Ogata et al., 1995）
他方，Gutenberg-Richter 独立分布でのマグニチュード系列から時空間 ETAS モデルでシミュ（- iii）が再現できる
レーションされた合成時空間データの前震からでも，重ね合わせの性質（i）．時空間 ETAS モデルは，単なる余震活動（Helmstetter et al., 2003; Ogata and Katsura, 2014）の重ね合わせなのに，なぜ前震の性質を再現するのか，この結果は，前震を物理学的に追究している研究者には意外に思われている．筆者は前震の物理的知見に基づく情報を持ち合わせていないが，この議論を深める鍵となるデータはマグニチュードの時系列である．その周辺分布であるが，その違いは，自然界の地震（気象庁カタロは何れも指数分布（Gutenberg-Richter 則）グ）には自己相関があるのに，前述の合成カタログのものは独立である．
そこで，筆者らは，一方で，気象庁カタログのマグニチュード時系列を逐次 ETAS に入力し，発生時刻と震央座標を作った第 1 種の合成カタログを作製し，他方では，気象庁カタログのマ
グニチュードをブートストラップ（boot strapping）して独立性を確保し，それらを ETAS に入力することよって第 2 種の合成カタログを作製した．それぞれの合成カタログは，いずれも時（- iii）の性質を満たすが，定量的には有意な違いがある．さらに，伝統的に知ら
空間的特徴
（i）
れている前震と余震の b 値の大小の関係については，本震の大きさなどの条件によって異なった結果を導くことを示した（Ogata and Katsura, 2014 参照）．
5.3 前震の事前識別ある地域で新規の地震群が始まったとき，それらが「前震」である確率をリアルタイムで求めたい．すなわち，群れの中で，それまでの最大マグニチュードより 0.45 以上飛びぬけた，格段に大きな地震が起こりうる確率である．先ず，群れの先頭の地震の位置座標，群内の各々の地
震の時間間隔，震央間距離，およびマグニチュード差を計算する．それらの履歴データを説明変数とし，最適な logit 関数のモデルによって，現在進行中の地震群が前震系列である確率，すなわち 1 か月以内にマグニチュードが，これまでより 0.45 以上の大きな地震が出現する確率を計算し予測する（統計数理 63 尾形論文 10 節参照）．
筆者らがそのような計算式を提案してから経過した 15 年間の予測の結果を検証した（Ogata．すなわち，経験的な一定確率での前震予測に対する対数尤度比基準（情報and Katsura, 2012）
利得）によって比較した結果，予測成績が上回るということになった．特に，本震が大きい地に限ると，識別予測の結果は極めて明瞭であった．
震の群（例えば M0 ≥ 6.5）
6.

多重確率予測式による大地震の予測

無情報の基での大地震の確率（永年確率）
は極めて小さいが，何らかの異常現象が現れると，

統計数理

224

第 67 巻 第 2 号 2019

その予測確率（確率利得）
は相対的に増えることが期待される．しかし，前兆性や切迫性の識別には不確定さが伴う．なので，我々は経時観測データの異常を明瞭に定義し，それらが前兆と）や確率利得を求める作業をする必要がある．また地して大地震に至る確率（「適中率，hit rate」
震の見逃し（不意打ち）を避ける為には，低い適中率でもできるだけ多くの異常現象を発掘する必要がある．
大地震に至る高い適中率や確率利得を出すことは困難であろうが，それでも独立な幾つかの異常の予測因子が重なれば，「多重確率予測公式」（統計数理 63 尾形論文）によって，大地震予測の確率が実用的なものまで高まることがある．独立性が保証できない場合は，多重確率予測関数の展開モデルで算出し，AIC 比較でモデルの期待予測を公式を一般化したロジット（logit）
高めることができる．
現状では，独立性を担保するような短期・中期・長期の予測確率を組み合わせるのが有望な策であり，先行研究では，宇津による 1978 年の伊豆大島近海地震，および安芸らによる中国の海城地震や唐山地震などに関して計算されている．最近著者も，M7.3 熊本地震に関して多重確率予測公式で試算した結果と予測確率のバラツキを示した（Ogata, 2017b）．
7.

おわりに

何時だったか思い出せないが「統計モデルはデータ解析の望遠鏡や顕微鏡である」という赤池の言がある．古来，科学的発見を進める原動力は，適切な問題意識の上に，望遠鏡や顕微鏡などの発見や改良のように，方法論的革新が大きく絡んでいることはよく知られている．また，科学的仮説の実証は予測の結果によって決着がつく．同様に，統計科学の研究対象は誠に複雑系ではあるが，適切な統計モデルは予測をそれなりに全うするだけでなく，モデルを作るにあたって考慮されていなかった
「異常」
や
「想定外」
とも呼ばれる科学的新事実を，現実のデータから露出することもある．統計モデルや，それに基づくグラフや画像による表現は，辛うじて見えるものや見えないものをはっきり見せ，新知見を導く科学的方法としての役割を果たしうる．
著者は地震統計にもとづく多くの経験則と，地震学の物理的仮説を，統計的点過程モデルとして表現して，統計的方法の有用性を示すように心がけてきた．たとえば ETAS モデルは地震のデータベースから短期発生率を予測するために，余震活動の経験則に基づいて構成されたが，地震活動の微妙な異常を検出する物差しとして使える可能性も提供する．地震活動研究と密接に関係する統計地震学の展開と，それらの予測における意義をお伝えできたなら幸いである．

参

考

文

献

Adamopoulos, L. (1976). Cluster models for earthquakes: Regional comparisons, Mathematical Geology, 8, 463-475.
安芸敬一 (1956). 統計地震学の現状, 地震, 8(4), 205-228, doi: 10.4294/zisin1948.8.4_205.
Engdahl, E. R. and Villasenor, A. (2002). Global seismicity: 1900-1999, International Handbookof Earthquake and Engineering Seismology (eds. H. K. Lee, H. Kanamori, P. C. Jennings,
et al.), Part A, Academic Press, Amsterdam, https://earthquake.usgs.gov/data/centennial/
centennial.pdf.
Field, E. H., Milner, K. R., Hardebeck, J. L., Page, M., van der Elst, N., Jordan, T. H., Michael, A. J.,
Shaw, B. E. and Werner, M. J. (2017). A spatiotemporal clustering model for the third uniform
California earthquake rupture forecast (UCERF3-ETAS): Toward an operational earthquake
forecast, Bulletin of the Seismological Society of America, 107(3) 1049-1081, doi:10.1785/

統計地震学の発展と地震活動予測：個人的経験と展望

225

0120160173.
Guo, Z. and Ogata, Y. (1997). Statistical relations between the parameters of aftershocks in time,
space and magnitude, Journal of Geophysical Research, 102(B2), 2857-2873.
Gutenberg, B. and Richter, C. F. (1944). Frequency of earthquakes in California, Bulletin of the
Seismological Society of America, 34, 185–188.
Hawkes, A. G. (1971). Point spectra of some mutually exciting point processes, Journal of the Royal
Statistical Society, Series B, 33, 438–443.
Helmstetter, A., Sornette, D. and Grasso, J. R. (2003). Mainshocks are aftershocks of conditional
foreshocks: How do foreshock statistical properties emerge from aftershock laws, Journal ofGeophysical Research, 108(B1), 2046, doi:10.1029/2002JB001991.
Jordan, T. H., Chen, Y.-T., Gasparini, P., Madariaga, R., Main, I., Marzocchi, W., Papadopoulos, G.,
Sobolev, G., Yamaoka, K. and Zschau, J. (2011). Operational earthquake forecasting: Stateof knowledge and guidelines for implementation, Final report of the International Commission on Earthquake Forecasting for Civil Protection, Annals of Geophysics, 54(4), 315–391,
doi:10.4401/ag-5350, http://www.annalsofgeophysics.eu/index.php/annals/article/view/5350/
5371.
熊澤貴雄 (2015). 地震活動の異常性とモデリング, 統計数理, 63(1), 45-64.
Kumazawa, T. and Ogata, Y. (2013). Quantitative description of induced seismic activity before andafter the 2011 Tohoku-Oki earthquake by non-stationary 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(3), 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.
Kumazawa, T., Ogata, Y., Kimura, K., Maeda, K. and Kobayashi, A. (2016). Background rates ofswarm earthquakes that are synchronized with volumetric strain changes, Earth and Planetary
Science Letters, 442, 51-60, doi:10.1016/j.epsl.2016.02.049.
Kumazawa, T., Ogata, Y. and Tsuruoka, H. (2017). Measuring seismicity diversity and anomalies
by point process models: Case studies before and after the 2016 Kumamoto earthquakes in
Kyushu, Japan, Earth, Planets and Space, 69, Article 169.
Lewis, P. A. W. (ed.) (1972). Stochastic Point Processes: Statistical Analysis, Theory and Applications,Wiley, New York.
Lewis, P. A. W. and Shedler, G. S. (1979). Simulation of nonhomogeneous Poisson processes bythinning, Naval Research Logistics Quartary, 26(3), 403-413, https://apps.dtic.mil/dtic/tr/
fulltext/u2/a059904.pdf.
Nishikawa, T. and Ide, S. (2018). Recurring slow slip events and earthquake nucleation in the source
region of the M7 Ibaraki-Oki earthquakes revealed by earthquake swarm and foreshock activity,
Journal of Geophysical Research, doi:10.1029/2018JB015642, https://drive.google.com/file/d/
16Ah30kGNuq3za038b2r_I-LNLSOvkGu5/view.
Ogata, Y. (1981). On Lewis’ simulation method for point processes, IEEE Transaction of Information
Theory, IT-27, 23-31.
Ogata, Y. (1983). Estimation of the parameters in the modified Omori formula for aftershock frequencies by the maximum likelihood procedure, Journal of Physics of the Earth, 31, 115-124,https://www.jstage.jst.go.jp/article/jpe1952/31/2/31_2_115/_article/-char/en.
尾形良彦 (1993). 地震学とその周辺の地球学分野に於ける統計モデルと統計的手法, 日本統計学会誌,

226

統計数理

第 67 巻 第 2 号 201922(3), 413-463, doi: 10.11329/jjss1970.22.413.
Ogata, Y. (1998). Space-time point-process models for earthquake occurrences, Annals of Institute of
Statistical Mathematics, 50(2), 379-402.
Ogata, Y. (2001a). Exploratory analysis of earthquake clusters by likelihood-based trigger models,Festschrift volume for Professor Vere-Jones, Journal of Applied Probability, 38A, 202-212.
Ogata, Y. (2001b). Increased probability of large earthquakes near aftershock regions with relative quiescence, Journal of Geophysical Research, 106(B5), 8729-8744, doi: 10.1029/2000JB900400.
Ogata, Y. (2005). Detection of anomalous seismicity as a stress change sensor, Journal of GeophysicalResearch, 110(B5), B05S06, doi:10.1029/2004JB003245.
Ogata, Y. (2007). Seismicity and geodetic anomalies in a wide area preceding the Niigata-Ken-Chuetsu
earthquake of 23 October 2004, central Japan, Journal of Geophysical Research, 112, B10301,
doi:10.1029/2006JB004697.

尾形良彦 (2008). 最近 30 年の大地震発生と指定地域について, 地震予知連絡会会報, 79(12-1), 623-625,
http://cais.gsi.go.jp/YOCHIREN/report/kaihou79/12_01.pdf.
Ogata, Y. (2010a). Space-time heterogeneity in aftershock activity, Geophysical Journal International,
181(3), 1575-1592, doi:10.1111/j.1365-246X.2010.04542.x.
Ogata, Y. (2010b). Anomalies of seismic activity and transient crustal deformations preceding the2005 M7.0 earthquake west of Fukuoka, Pure and Applied Geophysics, 167(8-9), doi:10.1007/
s00024-010-0096-y.
Ogata, Y. (2011a). 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 the2004 Chuetsu earthquake of M6.8, Geophysical Journal International, 186, doi:10.1111/j.1365246X.2011.05033.x.
Ogata, Y. (2011b). Significant improvements of the space-time ETAS model for forecasting of accuratebaseline seismicity, Earth, Planets and Space, 63(3), 217-229, doi:10.5047/eps.2010.09.001.
Ogata, Y. (2013). A prospect of earthquake prediction research, Statistical Science, 28, 521–541,
doi:10.1214/13-STS439.
尾形良彦 (2015). 地震の確率予測の研究 — その展望, 統計数理, 63(1), 3-27.
Ogata, Y. (2017a). Statistics of earthquake activity: Models and methods for earthquake predictability
studies, Annual Review of Earth and Planetary Sciences, 45, 497-527, doi:10.1146/annurevearth-063016-015918.
Ogata, Y. (2017b). Forecasting of a large earthquake: An outlook of the research, Seismological Research Letters, 88(4), 1117-1126, doi:10.1785/0220170006.
尾形良彦 (2017). 日本列島内陸部の常時地震活動度について, 地震予知連絡会会報, 97(9-12), http://cais.
gsi.go.jp/YOCHIREN/report/kaihou97/1_3.pdf.
Ogata, Y. (2018). Comment on “A review of self-exciting spatiotemporal point process and theirapplications” by Alex Reinhart, Statistical Science, 33(3), 319-322, doi:10.1214/18-STS650,
https://projecteuclid.org/euclid.ss/1534147222.
Ogata, Y. and Abe, K. (1991). Some statistical features of the long-term variation of the global andregional seismic activity, International Statistical Reviews, 59, 139-161, https://rmgsc.cr.usgs.
gov/outgoing/threshold_articles/Ogata_Abe1991.pdf.
Ogata, Y. and Katsura, K. (1991). Maximum likelihood estimates of the fractal dimension for random
spatial patterns, Biometrika, 78(3), 463-474.
Ogata, Y. and Katsura, K. (2012). Prospective foreshock forecast experiment during the last 17 years,
Geophysical Journal International, 191(3), 1237-1244, doi:10.1111/j.1365-246X.2012.05645.x.
Ogata, Y. and Katsura, K. (2014). Comparing foreshock characteristics and foreshock forecasting inobserved and simulated earthquake catalogs, Journal of Geophysical Research, 119(11), 8457-

統計地震学の発展と地震活動予測：個人的経験と展望

227

8477, doi:10.1002/2014JB011250.
Ogata, Y., Utsu, T. and Katsura, K. (1995). Statistical features of foreshocks in comparison with otherearthquake clusters, Geophysical Journal International, 121, 233-254.
Ogata, Y., Katsura, K. and Tanemura, M. (2003). Modelling heterogeneous space-time occurrencesof earthquakes and its residual analysis, Applied Statistics (Journal of the Royal StatisticalSociety, Series C ), 52(4), 499-509.
Ogata, Y., Katsura, K., Tsuruoka, H. and Hirata, N. (2018). Exploring magnitude forecasting of thenext earthquake, Seismological Research Letters, doi:10.1785/0220180034.
Ozaki, T. (1979). Maximum likelihood estimation of Hawkes’ self-exciting point processes, Annals of
the Institute of Statistical Mathematics, 31, 145–155.
庄建倉, 尾形良彦 (2015). 地震予測の評価法について, 統計数理, 63(1), 29-44.
宇津徳治 (1970). 5. 地震の時間的分布に関連する諸問題（その 3）：余震, 前震, 群発地震の時間的性質,
北海道大学地球物理学研究報告, 23, 49-71, doi: 10.14943/gbhu.23.49, http://hdl.handle.net/
2115/13973.
宇津徳治 (1999).『地震活動総覧』, 東大出版会, 東京.
Utsu, T., Ogata, Y. and Matsu’ura, R. S. (1995). The centenary of the Omori formula for a decay law of
the aftershock activity, Journal of Physics of the Earth, 43, 1-33, https://www.jstage.jst.go.jp/
article/jpe1952/43/1/43_1_1/_article/-char/en.
Vere-Jones, D. (2006). The development of statistical seismology: A personal experience, Tectonophysics, 413, 5-12.

228

Proceedings of the Institute of Statistical Mathematics Vol. 67, No. 2, 215–228 (2019)

The Development of Statistical Seismology:
A Personal Experience and View
Yosihiko Ogata
Professor Emeritus, The Institute of Statistical Mathematics

I provide an overview of the development of statistical seismology in Japan and myresearch experience. Some focuses are placed on prediction of seismic activity by point
process models, and statistical diagnostic analysis of anomalous seismic activities searching relations to physical phenomena.

Key words: Seismic activity, point processes, ETAS model, hierarchical Bayesian method, probability forecasts, foreshock discrimination.