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

12 －９

2023 年 3 月発行

背景地震活動度を用いた内陸地震の長期予測と検証評価
Prediction and validation of long-term earthquake probabilities in inland Japan
using the hierarchical space–time ETAS and space–time Poisson process models
尾形良彦（統計数理研究所）

Yosihiko Ogata (The Institute of Statistical Mathematics)
１．内陸部の大地震の長期確率予測のための４モデル
グーテンベルク・リヒター（G-R）則のｂ 値が 0.9 としてマグニチュード分布を仮定し以下の４種類の空間ポアソン過程モデルを考慮した．

先ず，大地震は「何時何処に起きてもおかしくない」と良く言われているが，この極端なモデルは時間的に一様で内陸部で場所を選ばず同じ確率で起きる一様ポアソン空間モデルである．この発生率密度は内陸部の各位置 (x ,y ) で同一値(a) λ ( x, y ) ∝ λˆinland × 10− b ( M c − 4.0)

で，その外域（海域）では発生率０と定義される．

次に考えられるのは時間的に一様だが非一様なポアソン空間モデルで(b) λ ( x, y ) ∝ λˆ ( x, y ) × 10− b ( M c − 4.0),

ただし λˆ ( x, y ) は局所線形デロネ関数（第 1 図参照）で表現され，それらの係数の推定は学習時間

区間 [S 0 , S ] の震源データ（例えば 1923 - 2018, M ≧ 4.0）にあてはめる（本稿第 1 図参照）．赤池ベ
イズ情報量規準（ABIC）の最小化によって，最適な平滑化の重み（超パラメタ）を求める．その重みで逆問題最適な事後分布関数の最大（MAP 解）を得る１- ３）

本稿では更に，階層ベイズ型的時空間 Epidemic-Type Aftershock Sequence (HIST-ETAS) モデル 1-5)−1t

K 0 ( xi , yi )  ( x − xi , y − yi ) S j ( x − xi , y − yi )+dλθ (t , x, y | H t ) =µ ( x, y ) + ∑p ( xi , yi ) α ( xi , yi )( M i − M 0 )e{i ; t <ti } ( t − ti + c )− q ( xi , yi )の誘発項（右辺第 2 項）を除去した背景（常時）地震活動密度 µˆ ( x, y ) を利用する．ここで誘発項を適切に除去するためには，震央座標 (x j , y j ) や分散行列 S j は準リアルタイム（例えば１時間以内）で検出地震群から AIC で必要に応じてセントロイド型座標に補正し計算・推定する（参考文

献５の第 3 図参照）．ここで，定式の３つのパラメタα, p , q が場所に依存せず定数の場合は HIST-ETAS- µΚ モデルと呼び，その場合の背景（常時）地震活動密度は⒞ λ ( x, y ) ∝ µˆ HIST − ETAS − µ K ( x, y ) × 10− b ( M c − 4.0)

で３番目のモデルとする．

３パラメタα , p , q が，場所 (x , y ) に依存して最適化した場合は HIST-ETAS-5pa モデルと呼び，対応する背景（常時）地震活動密度は(d) λ ( x, y ) ∝ µˆ HIST − ETAS −5 pa ( x, y ) × 10− b ( M c − 4.0)
これらのパラメタは (b) モデルと同様に，局所線形デロネ関数で表現され ABIC の最小化によって，最適なベイズ尤度を求める．但し先駆区間の情報データ 1885 － 1922 の宇津の明治大正カタログ６）― 591 ―

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

2023 年 3 月発行

も使用して，長期履歴に依存する HIST-ETAS モデルの定常性を担保する．

以上の (a), (b), (c), 及び (d) モデルのまとめられた図示は参考文献５の第 4 図参照．

２．予測モデルと結果の特徴について

空間的に非一様なポアソン時空間モデル (b) では λˆ ( x, y ) が学習期間の選び方によって，違った密度パターンになる（第 2 図）．それに対して，モデル (c) の µˆ HIST − ETAS − µ K ( x, y ) と (d) のµˆ HIST − ETAS −5 pa は学習期間の選び方によってロバストで，外見上互いに違わない密度パターンになる（本稿第 2 図および参考文献５の第 1 図参照）．

1923 年から 2018 年までの (b), (c), (d) モデルは一見，似たパターンであるが非一様 Poisson モデル (b) は HIST-ETAS のバックグラウンド密度モデル (c), (d) より，かなり集中的である．すなわち対数スケールで前者の等高線の数が，かなり多い（第 3 図）．

例えば M ≧ 4.0 の 学 習 期 間（ 推 定 期間 ） デ ー タ 1926 ~1995 年 で 得 ら れ た 予 測 モ デ ル (d) のµˆ HIST − ETAS −5 pa に対して 2000-2016 年の M6 以上の大地震（本稿第 4 図左パネル）は概ね予測リスクに良く乗っているように見える．特に従前から熊本地域は高リスクであった．他方，低リスクの予測に拘わらず，福島県浜通り，能登半島，福岡県西方沖が起きているが，これらは一世紀以上，カタログから地震活動が確認できなかった場所である．また 500-1884 年の歴史被害地震（本稿第 4図左パネル）も概ね予測リスクに良く乗っているように見える．ただし，東北地方北部では低リス
クの予測にたいして起きている．

以下，これらの予測結果を以下の定量的なスコアで評価する．

３．長期予測の検証評価法

内陸部地震活動の十分なデータが揃えば HIST-ETAS-µK モデルより，より柔軟な HISTETAS-5paの方が一般的に ABIC に基づく当て嵌まりが良い．しかし，最適な事後分布（MAP）解の予測モデルが良い予測結果になるとは限らない．そこで各種 MAP 解モデルで地震発生予測の結果から計算した以下の対数尤度スコアを評価基準として採用した．また全国一律にｂ ＝ 0.9 と仮定して，中規模から大規模の地震の予測の結果も比較した．

M4.0 以上の地震の予測後の発生結果に対して，G-R 則のｂ 値が全国一律 0.9 と単純に仮定した中規模以上（ M ≥ M c ）の地震予測成績を比較評価する．ただし，十分な数の大地震に対する遡及予測を十分な精度で評価するために，学習期間を 1923-1995 年に短縮し，1996-2021 年の長期予測を実施した．本稿では M ≥ M c の内陸の地震平均数を一定とする条件で空間配置の「対数尤度」
評価スコア５, ７, ８）λˆ ( xi , yi , M i )ˆlogf ( xi , yi , M i ) = log  ∏∏∞{i ; ( xi , yi )∈Inland , M i ≥ M c }{i; M i ≥ M c } ∫M ∫∫Inland λˆ ( x, y, M ) dx dydM cを用いる．このスコアが大きい程，予測が優れていたことを示す．ただし G-R 則の b = 0.9 を仮定
して，各サイズの大地震の結果を比べるため，モデルとして，

λˆ=( x, y, M ) λ ( x, y, 4.0) ×10−0.9 ( M −4.0)を評価する．

参考までに，確率予測の評価規準としての対数尤度の意味は文献８の第４節に記載されている．

― 592 ―

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

2023 年 3 月発行

すなわち対数尤度は元来「ボルツマンの相対エントロピー」と呼ばれるものであり，「尤度」が予測データの出現率を表す量であり，競合するモデル全部の尤度で正規化するとスコアは各モデルによる出現データの確率が得られる．以下の表では，モデル (a) の対数尤度を基準とした他のモデル(b) - (d) との尤度比が予測の相対的な良さの確度になる．
４．結果

HIST-ETAS モデルを気象庁震源データ 1923-1995 年にあてはめる場合，先行情報として宇津の明治大正カタログ６）を使う．1996 年以降の地震に対して各モデルについて，26 年間の遡及的予測と結果の評価を与えた（第 5 図参照）．

1996-2021 年の 26 年間の予測期間において，大地震（例えば Mc6.0 以上）の遡及的長期予測にはモデル (b)~(d) の活動度が高いところに多く発生している．評価スコアでも (d) HIST-ETAS-5pa の背景密度が最も優れた性能を示している．すると，M5.0 までの範囲では最適な非一様ポアソン過程モデルの予測が二つの HIST-ETAS 背景強度より優れており，空間的一様なポアソン過程モデルよりはるかに優れていることが示された．しかし M5.5 以上では HIST-ETAS-5pa モデルの背景空間強度 (d) が非一様ポアソン空間モデルよりも優れており，最高の性能を示した．しかし M7.0 以上については，データが少ないが，一様モデルの 10~20 倍の実現性が見られる．

最後に，1923 年から 2018 年の期間のデータによって推定された対象データの候補モデルについて，歴史的な被害地震９）（599 ～ 1884 年）の対数尤度スコアの検討を行った．歴史記録の精度を考慮して期間を 1585 年から 1884 年に限定し，北海道と東北地方北部を除いた，北緯 38 度以南の地域を限定した（第 7 図）と，HIST-ETAS-5pa モデルの背景地震活動の空間強度は全てのマグニチュードの範囲において最もスコアが高く，他のポアソン空間モデルよりもはるかに優れていることが分かった．M8 クラス以上の地震については，歴史的な被害地震の記録数が少ないため，評価スコアの差は小さく，地域的な有意差は大きくなかった．

第 7 図は確率をセルごとに示している．特に M6 以上の直下型地震の 30 年確率は凡そ 0.5% ~10% の地域性が見られる．
５．追記

最近 2 年 9 ヶ月の短期予測結果は HIST-ETAS-5pa がマグニチュード 4.0 から 5.0 まで一番優れており HIST-ETAS-µK がそれに次ぐ．中期予測に限れば ( ｂ ) 非一様ポアソンモデルが最も優れている５）．

HIST-ETAS-5pa モデルを気象庁震源データ 1923-2010.5 年にあてはめる場合，先行情報として宇津の明治大正カタログ６）に加えて，宇津の歴史地震９）も考慮すると，当てはまりの ABIC は -325.8ほど改善するが背景（常時）地震活動密度 ˆ ( x, y ) は殆ど変わらなかった．

本解析に関して気象庁震源カタログおよび地震活動可視化システム TSEIS を使用した．
参考文献
1) Ogata, Y et al. Applied Statistics (JRSSC) 52 (4), 499-509, 2003.
2) Ogata, Y. Earth, Planets and Space 63 (3), doi:10.5047/eps.2010.09.001, 2011.
3) Ogata Y et al. ISM Computer Science Monographs 35, https://www.ism.ac.jp/editsec/csm/pdf/csm_035.
pdf, 2021.
― 593 ―

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

2023 年 3 月発行

4) 統 計 数 理 研 究 所 ( 尾 形 良 彦 ) 予 知 連 会 報 97 (1-3), 2017. https://cais.gsi.go.jp/YOCHIREN/report/
kaihou97/01_03.pdf

5) 統計数理研究所 ( 尾形良彦 ) 予知連会報 107 (12-8), 2022. https://cais.gsi.go.jp/YOCHIREN/report/
kaihou107/12_08.pdf

6) 宇津徳治，地震研究所彙報 ，57, 401-463, 1982. https://repository.dl.itc.u-tokyo.ac.jp/records/33001#.
Y2NyU-TP0uU
7) Ogata Y et al. Bull. Seismol. Soc. Am. , 103 (3), doi:10.1785/0120120063, 2013.

8) 庄 建倉，尾形良彦 . 統計数理 63 (1), https://www.ism.ac.jp/editsec/toukei/pdf/63-1-029.pdf, 20159) 宇津徳治，世界の被害地震の表 , http://iisee.kenken.go.jp/utsu/index.html, 1990.

― 594 ―

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

2023 年 3 月発行

Delaunay 型
部分線形関数:

第1図

左図の点は学習期間（1923 年〜 2018 年）の日本内陸部で発生した M ≧ 4.0 および宇津カタログ（明治・大正）の地震の震央位置を示している．右図では，震央（左図の点）と内陸域外の追加点を結ぶデロネ分割を示している．デロネ三角形は補間の為の局所的な線形関数（デロネ関数；右パネルの挿入図）を定義する．最適な平滑化は赤池ベイズ情報量規準（ABIC）によって決定され，デロネ関数の係数は最大事後分布解を逆問題により求める 1-3)．

Fig. 1The points in the left ﬁgure show the epicenter locations of M ≥ 4.0 earthquakes that occurred in inland Japan during
the learning period (1923-2018) in addition to Utsu’s data for 1885-1922. The right figure shows the Delaunay
partition connecting the epicenters (the points in the left ﬁgure) and additional points outside the inland interior. The
Delaunay triangle deﬁnes a local linear function of interpolation. The optimal smoothing is determined by the Akaike
Bayesian Information Criterion (ABIC) and the coefficients of the Delaunay function are obtained by an inverse
problem1-3) , which is called the maximum a posteriori (MAP) estimate.

1923 ~2018

1926 ~1995

Log10( M≧4 / deg^2 / day)Log10(M≧4 / deg^2 / day)

第2図

パネルに示した異なる期間の内陸部のデータに対する非均質ポアソン空間モデルの MAP 推定値．色と等高線は M ≥ 4 地震 /deg2/day の予想数の対数スケールである．

Fig. 2MAP estimates of the non-homogeneous Poisson space model for the inland data for the diﬀerent periods shown in
the panel. Colors and contours are on the log scale of the expected number of M ≥ 4 earthquakes/deg2/day.

― 595 ―

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

非一様 Poisson モデル

HIST-ETAS-Kモデル

1923 ~2018

1923 ~2018

Log10(M≧4 地震 / deg^2 / day)

HIST-ETAS-5paモデル

1923 ~2018log10 ˆ  K ( x, y )log10 ˆ ( x, y )2023 年 3 月発行

Log10(M≧4 地震 / deg^2 / day)log10 ˆ 5 pa ( x, y )Log10(M≧4 地震 / deg^2 / day)

左図の非均質空間ポアソンモデル（本文 b モデル），中央図の HIST-ETAS-µK モデル（本文 c モデル），右図の HIST-ETAS-5pa モデル（本文 d モデル）の期間中の同じデータに対する MAP 推定値．色と等高線は
図 2 のものと同じである．

Fig. 3

MAP estimates of the non-homogeneous spatial Poisson model (b model in text) in the left panel, the HIST-ETAS-µKmodel (c model in text) in the middle and the HIST-ETAS-5pa model (d model in text) in the right, for the same datain the period. Colors and contours are the same as those in Fig. 2.

M≧4.0／２）／日／km^２

第3図≧6.8 地震○: M≧
発生期間 599 - 1884

≧6.0 地震○: M≧
発生期間 2000 - 2016

第4図

1926-1995 年の M ≥ 4 の地震データに対する HIST-ETAS-5pa モデルの推定背景強度画像．画像カラーはカラーバーに示すように 10-8 ~ 10-6/km2/day の範囲．画像上の円は (a) 2000 年から 2016 年の M ≥ 6 の浅い地震，(b) 599 年から 1884 年の M ≥ 6.8 の歴史的被害地震 9) の位置を示している．

Fig. 4

Estimated background intensity image of the HIST-ETAS-5pa model for M ≥ 4 earthquake data for the period of interest
1926-1995 where image colors range 10-8 ~ 10-6 /km2/day as shown in the color bar. The circles on the image show (a)shallow M ≥ 6 earthquakes for period 2000 -2016 and (b) the locations of historical damaging earthquakes of M ≥ 6.8
for the period from 599-1884.

― 596 ―

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

long-term
forecast

学習・推定

準備区間

長期予測(26y)

1885

1926

予測・検証
1995 1996≧4M≧

M≧≧~６６ R

2023 年 3 月発行

S2021

T

1926 ~1995

1926 ~1995(d) HISTETAS-5pa(ｂｂ) Non-uniform
PoissonLog (M≧4 / deg^2 / day)Log10(M≧4 / deg^2 / day)10

M ≥ 4.0

M ≥ 4.5

number of events276599030510343185(a) inland uniform

0.00.0

0.0

0.0

0.0

0.00.0(b) non-uniform

2835.01024.9

264.8

60.1

12.3

9.8

2.5⒞ HIST-ETAS-m
mK

2270.7

843.5

230.0

53.7

16.9

9.2

3.0(d) HIST-ETAS-5pa

2576.7

954.1

263.7

61.6

19.9

10.6

2.5

Forecast magnituderange

M ≥ 5.0 M ≥ 5.5 M ≥ 6.0 M ≥ 6.5 M ≥ 7.0

第5図

学習期間を 1995 年までとして，1996 年から 2021 年までを予測・検証した．上段図： 一，二位を競う非一様ポアソンモデル（左）と HIST-ETAS-5pa の背景強度密度モデル（右），それぞれの画像と等高線．下段の表：各マグニチュード範囲の地震予測の対数尤度スコア．スコアは内陸一様ポアソンモデル（a）のスコアとの相対値を示している．赤い数字は最も良い予測，青い数字は 2 番目の予測である．

Fig. 5The study period was set to 1995, and predictions were made and validated from 1996 to 2021. Top panels: Images
and contours of the non-uniform Poisson model and the HIST-ETAS-5pa background intensity density model,respectively, competing for ﬁrst and second place. Bottom table: Log-likelihood scores for earthquake predictions ineach magnitude range. Scores are relative to the score of the "inland uniform Poisson" model (a). Red numbers are thebest predictions; blue numbers are the second best predictions.

― 597 ―

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

CUMLATIVE NUMBERS OF EARTHQUAKES

予測・検証
1585

599

M8M7M6M5

Forecast of magnitude range All events M ≥ 5.5
number of events
159
156
0.00.0(a) uniform in inland Japan
-51.2
-53.3
(b) non-homogeneous
9.45.0(c) HIST-ETAS-m
mK BackgroundYEAR(d) HIST-ETAS-5pa Background

推定
1884

R599 ~ 1884

13.5

9.1

2023 年 3 月発行

1923 2019

STM ≥ 6.0 M ≥ 6.5 M ≥ 7.0 M ≥ 7.5
1380.0-54.52.3
7.8

700.0
-38.1
-3.7-0.5

310.0
-24.9
-4.6

40.0
-5.6
-2.4

-2.6

-1.9

1585 ~ 1884 北緯38°以南（北海道，北東北を除く）〇 599-1584, M≥6〇 1585-1884, M≥6
歴史被害地震（宇津）,○ 599 – 1499; ○ 1500 - 1884

Forecast of magnitude range All events M ≥ 5.5 M ≥ 6.0 M ≥ 6.5 M ≥ 7.0 M ≥ 7.5
114
57
25
3
number of events
131
129
0.0
0.0
0.0
0.0
0.0
0.0
(a) uniform in inland Japan
-7.8
-11.3
-15.7
-16.3
-12.4
-1.2
(b) non-homogeneous
32.9
28.1
24.0
9.84.01.1(c) HIST-ETAS-m
mK Background
36.7
31.8
29.4
12.5
5.5
1.6(d) HIST-ETAS-5pa Background

第6図

配色と等高線は 1926 年から 2019 年の M ≥ 4 地震データに対する HIST-ETAS-5pa モデルの推定背景密度．
挿入図は検出された歴史被害地震の時代対地震累積数と M-T 図．画像上の青丸は 599 年から 1584 年，赤丸は 1585 年から 1884 年の歴史的な被害地震の位置と規模をそれぞれ示している．上のスコア表は全ての歴史的被害地震記録に対する各モデルの各マグニチュード範囲の対数尤度評価，下の表は 1585 年から1884 年までの北緯 38 度以南に制限された地震に対する各モデルの対数尤度評価である．

Fig. 6

Color scheme and contours are the estimated background density of the HIST-ETAS-5pa model for the 1926-2019 M≥ 4 earthquake data. The inset shows the cumulative number of earthquakes versus the age of the detected historicaldamage earthquakes and the M-T diagram. Blue and red circles on the image indicate the location and size of theearthquakes from 599 to 1584 and 1585 to 1884, respectively. The top score table is the log-likelihood for eachmagnitude range rating of each model for all historical damage earthquake records, and the bottom table is the loglikelihood rating of each model for earthquakes restricted south of 38° N from 1585 to 1884.

M４以上の１日確率10��10��M６以上の30年確率(%) 0.510��1.05.010.

50.

第7図

全国一律の b 値を 0.9 とすると，左図は，日本内陸部の 0.1°
× 0.1°（約 100km2）の各セルで 1 日に M ≥ 4 の地震が発生する確率の常用対数を示している．右図は，今後 30 年間に 0.2°× 0.2°（約 400km2）の各セルで M ≥ 6 の地震が発生する予測確率 (%) を示す．

Fig. 7

Assuming a uniform nationwide b value of 0.9, the left ﬁgure shows the ordinary logarithm of the probability of an M ≥ 4earthquake occurring in each 0.1°
× 0.1°(approx. 100 km2) cell in inland Japan on a single day. The ﬁgure on the rightshows the probability of an M ≥ 6 earthquake occurring in each 0.2°
× 0.2°(about 400 km2) cell over the next 30 years.

― 598 ―