地震予知連絡会会報第 107 巻

12 －８

2022 年 3 月発行

階層的時空間 ETAS モデルなどによる短・中・長期の地震確率予測と検証評価
Prediction and validation of short- medium- and long-term earthquake
probabilities using a hierarchical space-time ETAS (HIST-ETAS) models, etc.
尾形 良彦（統計数理研究所）
YOSIHIKO Ogata (The Institute of Statistical Mathematics)

広域の地震活動度予測のための短期・中期・長期モデルの推定

現在進行中の地震活動の震源データから，各地の時空間地震活動の短期予測を実行できるようにするため階層ベイズ型的時空間 Epidemic-Type Aftershock Sequence (HIST-ETAS) モデル 1, 2)−1t
K 0 ( xi , yi )  ( x − xi , y − yi ) S j ( x − xi , y − yi )+λθ (t , x, y | H t ) =µ ( x, y ) + ∑dp ( xi , yi )eα ( xi , yi )( M i − M 0 ){i ; t <ti } ( t − ti + c )− q ( xi , yi )を使う．ここで震央座標 (x j , y j ) や分散行列 S j は準リアルタイム（例えば 1 時間以内）で検出地震群から AIC で必要に応じてセントロイド型座標に補正される 3)．常時活動度の推定されたパラメタµ (x , y ) は永年確率の候補として使用される．パラメタ K 0 (x , y ) は短期予測に有用な余震生産性の位置情報になり，本震直後（たとえば一時間後）からの余震発生の時空間的予測に使用される．

本稿で パラメタ α , p , q が，上式のように場所 (x , y ) に依存する場合は (1) HIST-ETAS5pa モデルと呼び，定数の場合は (2) HIST-ETAS-µΚ モデルと呼ぶ．これらのパラメタは局所線形デロネ関数で表現され，それらの係数の推定はターゲット時間区間 [S 0 , S ] の気象庁震源データ（例えば1923 - 2018, M ≧ 4.0）にあてはめて，赤池ベイズ情報量規準（ABIC）の最小化によって，最適なベイズ尤度を求める 4)．但し先駆区間の情報データで，例えば 1885 － 1922 の宇津カタログで，地震活動度モデルの長期依存性を担保する．しかる後に，逆問題で局所線形デロネ関数の最適な事後分布関数の最大（MAP 解）を得る 4)．その係数の条件付強度 λˆ (t , x, y | H t ) で短期予測を実施する．
両モデルとも，既に CSEP 検証センター（地震研究所）に予測プログラムを提出し，比較検証が継続中であるが，それらの時間と空間の精度を上げた予測動画を今回の予知連会議で上映したが，第3 図にそれらのスナップショットを示した．

同様に，履歴に依存せず時間的に定常だが空間的に非一様な，ポアソン時空間モデルを 4 種類考慮した（第 4 図）．先ず，直下型が場所を選ばず同じ確率で起きる，内陸部のみで同一値の発生率でその外域では発生率 0 の (3) 内陸部一様ポアソン過程モデル，次に ABIC 最小化で求めた局所線形デロネ関数の (4) 空間非一様ポアソン過程モデル，そして (5) HIST-ETAS5pa モデルの常時活動度補正と (6) HIST-ETAS-µΚ モデルの常時活動度補正の，いずれも空間非一様ポアソン過程モデルである．

以上の 6 モデル含む計算法は FORTRAN ソフトウェア集 5) として公開されている．

検証評価法と結果

これまでの著者の計算経験によれば，殆どの広域地震活動データからの推定において ABIC に基づくと，HIST-ETAS-µK モデルより，より一般的な HIST-ETAS5pa モデルの方が当て嵌まり良い結果になっている．

― 547 ―

地震予知連絡会会報第 107 巻

2022 年 3 月発行

しかし，ベイズ適合性の ABIC 最小法で求めた事後分布モデルの最大値解（MAP 解）のパラメタ値を挿入した予測モデル (MAP 解 plug-in model) が，予測スキルで勝るとは限らない．そこで時刻S までのデータから求めた各種モデルの MAP 解（ θˆ ）モデルで時間区間 [S , T ] での発生予測とその結果から計算した対数尤度スコアを評価基準として採用する．これはlog L(λˆ; S , T , M c )ただし∑{i ; S <ti <T , M i ≥ M c }log λ (ti , xi , yi , M i | H ti ) − ∫∞McT∫ ∫∫SAreaλ (t , x, y, M | H t ) dx dy dt dM| H t ) λθˆ (t , x, y | H t ) ×10− b ( M − 4.0) andλˆ (t , x, y, M==
H t {(t j , x j , y j , M j ); t j < t , M j ≥ 4.0}cである．ここで履歴情報 H t は M4.0 以上の地震の時刻 t までの発生情報である．また G-R 則のｂ値が全国一律 0.9 と仮定した中規模以上（ M ≥ M c ）の地震予測の結果に対しても比較評価できる．

ポアソン過程モデルの場合，履歴情報 H t は不要である．

気象庁震源データ（先行情報として宇津カタログを含む）に基づき，2019 年以降の地震について予測と結果の評価を与えた（第 5 図）．表の評価では，予測確率が変わらない内陸部一様ポアソン過程モデル（3）を基準とした対数尤度比スコア 6) で短期予測の結果を評価したが，中期予測ではポアソンモデル同士のみの比較になる．

最近 2 年 9 ヶ月間の短期予測結果は (1) HIST-ETAS-5pa モデルが M5 クラスまでの地震に対する予測結果が最も優れており (2) HIST-ETAS-µK モデルがそれに次ぐ．中期予測に限れば，(4) 定常非一様ポアソンモデルが他より優れており，特に (3) 内陸部均一ポアソンモデルより遥かに優れている．

大地震（例えば Mc6.0 以上）のレトロスペクティブ長期予測には常時活動度が高いところに多く発生しており 7)，見かけでは既にモデル（5）や（6）が有望であるが，今回の 3 年足らずの検証では評価すべき予測の地震数が足りなく，数年足らずの予測では長期予測は難しい．

その代わりに，歴史地震に対して「逆予測」を試みる．すなわちポアソン過程の特性から，時間軸の因果性を無視できる．そこで空間分布モデルとしての逆予測はNlog L(λˆ; Area, M c ) = log ∏i =1∞∫ ∫∫McAreaλˆ ( xi , yi , M i )λˆ ( x, y, M ) dx dydMをスコアとして評価できる．ただし， λ̂ は気象庁カタログ（1923 - 2018）から求めたモデルを意{}

味し，( x j , y j , M j ); M j ≥ M c は宇津 9) の直下型被害歴史地震（古代 ~ 1884）の配置（第 6 図）である．
こうして当てはまりの性能を比較したスコアの表も第 6 図にある．結果，HIST-ETAS-5pa 時空間モ

デルの常時地震活動の空間 µ 5pa (x ,y ) 分布が，非一様ポアソン空間モデル λˆ ( x, y ) を大きく離して，最も良い結果を導いている．M7.5 クラス以上の歴史地震では，小標本の為かスコアの差が小さく，地域的な特徴が見られない．本報告では内陸部直下型の地震予測を議論したが，海域を含む全日本の予測 8) でも同様な検証が可能である．

本解析に関して気象庁震源カタログおよび地震活動可視化システム TSEIS を使用した．
（尾形良彦）

OGATA Yosihiko― 548 ―

地震予知連絡会会報第 107 巻

2022 年 3 月発行

文献1) 統計数理研究所 ( 尾形良彦 ), 予知連会報 96 (12-23), 2016.

2) 尾形良彦 . 統計数理 63 (1), https://www.ism.ac.jp/editsec/toukei/pdf/63-1-029.pdf, 20153) Ogata, Y. Earth, Planets and Space 63 (3), doi:10.5047/eps.2010.09.001, 2011.
4) Ogata, Y., et al. Applied Statistics (JRSSC) 52 (4), 499-509, 2003.
5) Ogata et al. ISM Computer Science Monographs 35, https://www.ism.ac.jp/editsec/csm/pdf/csm_035.pdf,
2021.
6) Ogata, Y., et al. Bull. Seismol. Soc. Am. , 103 (3), doi:10.1785/0120120063, 2013.
7) 統計数理研究所 ( 尾形良彦 ) 予知連会報 97 (1-3), 2017.

8) 統計数理研究所 ( 尾形良彦 ) 予知連会報 79 (12-1), 2008.

9) 宇津徳治，世界の被害地震の表 , http://iisee.kenken.go.jp/utsu/index.html, 1990.

― 549 ―

地震予知連絡会会報第 107 巻

第1図Fig. 1

2022 年 3 月発行

Ogata (1998) の時空間的 ETAS モデルの特徴は，震源カタログの震央をセントロイド型に自動補正することである．震央座標（x j , y j ）と分散行列 S j は，大地震（例えば M5 以上）の後に例えば 1 時間スパンに検出された全ての地震の震央座標を基に，AIC を用いて同定する．上図の赤丸はカタログに掲載されている震央，星印は推定重心座標，点線は行列 S j に対応する推定楕円である．
A characteristic feature of the spatio-temporal ETAS model of Ogata (1998) is the automatic centroidal correctionof epicenters in the hypocenter source catalog. The source coordinates (x j , y j ) and variance matrix S j are identiﬁedusing AIC based on the source coordinates of all detected earthquakes in a short span, say one hour, after a major
earthquake (e.g., M5 or greater). The red circles in the above panels are the epicenters listed in the catalog, the stars
are the estimated centroidal coordinates, and the dotted curves are the estimated ellipses corresponding to the matrix
Sj .

― 550 ―

地震予知連絡会会報第 107 巻

第2図Fig. 22022 年 3 月発行

各地域の時空間的な地震活動の短期的な予測を可能にするために，階層型ベイズ時空間 (HIST-ETAS) モデルを採用した．まず，パラメタ µ (x , y ) は時間に依存しない一定の活動の通年確率の密度関数とする．また，パラメタ K 0 (x , y ) は短期予測に有用な余震生産性の位置情報を提供する．パラメタ α , p , q が上式のように位置に依存する場合を (1) HIST-ETAS5pa モデル，定数の場合を (2) HIST-ETAS-µK モデルと呼ぶ．これ
らのモデルの係数の推定を，十分に長い目標時間間隔 [S 0 , S ] の気象庁元データ (1923 - 2018, M ≧ 4.0) に適用し，最適なベイズ尤度を赤池ベイズ情報量規準 (ABIC) で決定し，逆問題により局所線形ドローネ関数の最適事後分布の最大パラメタ（MAP）解を求める．
Hierarchical space-time epidemic-type aftershock sequence (HIST-ETAS) models are adopted to enable efficient
short-term forecasts of spatiotemporal seismic activity in each region. First, the parameter µ (x , y ) does not dependon time is used as a candidate for the perennial earthquake probability. The parameter K 0 (x , y ) provides position
information for aftershock productivity useful for short-term prediction. When the parameters α , p , and q depend on
the location as in the above equation, it is called (1) HIST-ETAS5pa model; and when it is a constant, it is called (2)HIST-ETAS-µK model. The estimation of the coeﬃcients of these models is applied to the JMA hypocenter data (1923- 2018, M ≧ 4.0) in a suﬃciently long target time interval [S 0 , S ], and the optimal Bayesian likelihood is determinedby the Akaike Bayesian Information Criterion (ABIC). After that, the maximum a posteriori parameter (MAP)
solution of the optimal posterior distribution of the local-linear Delaunay function is obtained by the inverse problem.

― 551 ―

地震予知連絡会会報第 107 巻

2022 年 3 月発行

1923−2018, M≧4

2019−2021.9.30, M≧4第3図Fig. 3

HIST-ETAS モデルで短期予測は，既に CSEP 日本（地震研究所）に予測プログラムを提出し，検証が継続中であるが，それらの時間と空間の精度を上げた予測の動画を今回の予知連会議で上映し，そのスナップショットを右側に示す．左側の図は示した時間区間での震央図である．
Short-term predictions and veriﬁcations the HIST-ETAS models have been ongoing in the CSEP Japan Testing Centerat the Earthquake Research Institute, University of Tokyo. Here are snapshots that are taken from the screening
predictive videos with improved accuracy of predictions in times and spaces that were carried out at the present CCEP
meeting. The ﬁgures on the left are the epicentral map for the corresponding time interval.

― 552 ―

地震予知連絡会会報第 107 巻

第4図Fig. 42022 年 3 月発行

履歴に依存せず , 時間的に定常で , 空間的に非一様なポアソン時空間モデルを 4 種類考慮した．先ず，「(3) 内陸部一様ポアソン過程モデル」は内陸部のみで一定発生率でその外域は 0 の発生率，「(4) 空間非一様ポアソン過程モデル」は全ての M ≧ 4 地震データから ABIC 最小化で求めた最適 MAP 解，そして「(5) HIST-ETAS-µK モデルの常時活動度」を定数倍補正した空間非一様ポアソン過程，最後に「(6) HISTETAS5pa モデル」の常時活動度を定数倍補正した空間非一様ポアソン過程である．右側の図 (3) ～ (6) は，日本の内陸部における 4 つの対応する空間ポアソン強度率を対数で表したものである．
Four nonhomogeneous Poisson space-time models that are time-independent but spatially non-uniform. Those are:(3) inland uniform Poisson process model; (4) non-uniform spatial Poisson process model that minimize the Akaike
Bayesian Information Criterion (ABIC); (5) proportionally corrected background µ (x ,y ) intensities of the HISTETAS-µK model and (6) HIST-ETAS -5pa model, respectively. The ﬁgures (3) - (6) on the right are contours of the 4corresponding spatial Poisson intensity rates on Japan inland area in logarithmic scale.

― 553 ―

地震予知連絡会会報第 107 巻

2022 年 3 月発行

対数尤度スコア: ln L(S,T | M≧
≧4.0 の履歴 )Mc4.0 Mc4.5

時空間点過程モデル
地震数(2019-2021.sept)

126

42123

(1) HIST-ETAS-5pa

638.7

276.6

55.8-3.9(2) HIS-TETAS-k551.2

252.5

49.9-3.2

inland : 内陸部 ⼀様 Poisson

0.0

0.0

0.0

0.0

(x,y) : ⾮⼀様 Poisson

157.9

105.6

41.7

9.3

(x,y)-Poisson (HIST-ETAS-5pa)

64.7

85.2

36.3

8.3

 (x,y)-Poisson (HIST-ETAS-)

111.5

91.5

35.3

7.9



1885

0第5図Fig. 5Mc5.0 Mc5.5

M≧~６2021.9.30
2019

1923

RM≧4S

T

一般に HIST-ETAS-µK モデルより，より柔軟な HISTETAS-5pa の方が ABIC に基づく当て嵌まりが良い．
しかし，MAP パラメタを挿入した予測モデル (plug-in model) が良い結果になるとは限らない．そこで時刻S までのデータから求めた各種 MAP 解モデルで時間区間 [S , T ] 間での地震発生予測とその結果から計算した以下の対数尤度スコアの大きさを評価基準として採用した．また全国一律ｂ ＝ 0.9 と仮定して，中規模以上の地震の予測の結果も比較した．最近 2 年 9 ヶ月の短期予測結果は HIST-ETAS-5pa がマグニチュード 4.0から 5.0 まで一番優れており HIST-ETAS-µK がそれに次ぐ．中期予測に限れば「(4) 定常非一様ポアソンモ
デル」が最も優れている．
In general, the more ﬂexible HISTETAS-5pa ﬁts better based on ABIC than the HIST-ETAS-µK model. However, thepredictive model (plug-in model) with the MAP parameter inserted does not always give best forecast. Therefore, thelog-likelihood score calculated from the prediction of earthquake occurrence during the time interval [S , T ] and theresult of various MAP solution models obtained from the data up to time S is adopted as the evaluation criterion. Also,assuming b = 0.9 nationwide, the forecast evaluation results of the higher score were also compared based on theoccurrence information of earthquakes of magnitude M4.0 or higher. As for the short-term forecast results for the last
2 years and 9 months, HIST-ETAS-5pa is the best for forecasting the lower limit magnitude from 4.0 to 5.0, followed
by HIST-ETAS-µK . As far as the medium-term forecast is concerned, the MAP of non-uniform Poisson model issuperior to the other models including the inland uniform Poisson model.

― 554 ―

地震予知連絡会会報第 107 巻

2022 年 3 月発行○: M≧6.０ 歴史被害地震（宇津）, 599 ‐ 1884

Spatial probability density; log-likelihood scoreˆ ( xi , yi , M i )Nlog L(ˆ; Area, M c )  log 
ただし ˆ は下記モデル

i 1



 McArea

ˆ ( x, y, M ) dx dydM

( x y , M ); M  M  は被害歴史地震 (古代~1884) の配置
j,jjjcUtsu Historical data;
Mc 4.0Mc 6.0Mc 6.5

Mc 7.0

number of events

206

185

101

50

8

non-homogeneous (4)-911.4-828.0-411.4-233.3-36.4

magnitude thresholds

第6図Fig. 6○: M≧6.8 歴史被害地震（宇津）, 599 ‐ 1884

Mc 7.5

HIST-ETAS-K Background (5)-826.4-746.8-412.2-207.2-33.7

HIST-ETAS-5pa Background (6)-820.4-739.4-408.1-205.8-34.8

uniform in inland Japan(3)-881.6-791.7-432.2-214.0-34.2

uniform in all Japan region-1210.8-1087.4-593.7-293.9-47.0

被害歴史地震に対する気象庁カタログ（1923 - 2018）による長期的逆予測の性能評価．右上パネルの赤丸
印は M6 クラス以上の直下型歴史被害地震 9) による．左上の評価式は気象庁データ (1923-2018, M ≧ 4.0)に基づいて推定した各種非一様ポアソン空間分布モデルによる，歴史地震の配置データに対する当てはまりを比較した評価スコア式．左下の表は各々の結果スコア．HIST-ETAS5pa モデルの µ 値の空間変化が最も良い結果となる．M7.5 クラス以上の歴史地震では小標本数の為かスコアの差が小さく，発生率の地域差が明確でない．右下パネルは，M6.8 以上の歴史被害地震の発生場所が µ 値（色標は M ≧ 4 の地震 /km2/day）の高い地域で発生し易いことを示している．
Performance evaluation of long-term‘reverse prediction’of historical earthquakes based on the Japan Meteorological
Agency hypocenter catalog (1923 – 2018, M ≧ 4.0). The red circles on the upper right panel are shallow historicaldamage earthquakes of M6 class or higher 9). The upper left panel is an evaluation score statistic to compare the ﬁtof the spatial Poisson models to the coordinate data of historical earthquakes. The lower left panel is the table of theresults. The µ value change of the HIST-ETAS5pa model gives the best results, but the earthquakes of M7.5 class orabove, the diﬀerence in scores is small possibly due to a small number of such events; namely, regional diﬀerences
cannot be seen. The lower right panel shows the location of historical damaging earthquakes of M6.8 or higher andthe µ value distribution of the HIST-ETAS5pa model (color table indicates expected number of M ≧ 4 earthquakes /km2 / day).

― 555 ―