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

2021 年 3 月発行

12 － 10 階層的時空間 ETAS モデルに基づく短期・中期・長期予測および背景率予測― 自動予測の開発に向けて ―	Short-, medium-, long-term and background-rate forecasts based on a hierarchicalspace-time ETAS model: Towards the development of automated forecasts
尾形良彦（統計数理研究所）

Yosihiko Ogata, The Institute of Statistical Mathematics
1．リアルタイム確率予測を目指して地震予測の的中率を高めるために宇津（1979）・安芸（1981）が「多項目確率予測式（multielements prediction formula）」1, 2) を提案し 40 年経った．従来は，いろいろな観測項目について「前兆」現象を捜して，ある項目について明瞭な異常が認められたときには，その地域周辺で観測を強化して他の項目についても調べ，総合的に判断をするという方向であった．結果として，地震の予知率が極めて低かったといえる．

それゆえ，予知率を高めるために，予め永年確率や各異常現象の確率を与える準備や環境を整えておく必要がある．今日では，大量のデータが殆どリアルタイムで利用可能であり，広領域での地震確率のリアルタイムの予測は喫緊の課題である．機会を逃がさず，大地震の見逃しが少なくなるように，とくに短期・中期予測の試行を積み重ねる必要がある．
点過程の条件付強度関数による多項目確率予測式は

λ (t , x, y, M | H t , Ft ) = µ ( x, y, M )λETAS (t , x, y, M | H t ) K λk (t , x, y, M | Ft k )∏µ ( x, y , M )µ ( x, y , M )k =1(1)の如く，時・空間・マグニチュードの表現 3, 4) ができる．ここで µ (x, y, M) は永年確率の空間・マグニチュードの背景率（ポアソン過程； 第 1 図参照）で，(1) 式の各々の項の比は「危険拡大度（riskenhancement factor）5)」と呼ばれる．これらによって諸異常事象が大地震の前兆かもしれない可能性，その時間的切迫度と地域的影響範囲に関する条件つき時空間的強度関数のリアルタイム・モニタリングを実施する．例え的中率が僅かでも，異常事象に基づく条件つき大中規模地震の強度関数のモデル化の開発を進め，検証・評価・改善する根拠になる．
2．広域の地震活動度のモニタリングと診断解析
先ず現在進行中の地震活動の震源データベースから，時空間の地震活動度のリアルタイム予測可視化モニターシステムを構築し，各地の時空間地震活動の短期予測を実行できるようにする．そのため階層ベイズ型的時空間 Epidemic-Type Aftershock Sequence (HIST-ETAS) モデル 2, 3, 20)

 ( x − xi , y − yi ) S −j 1 ( x − xi , y − yi )tλ (t , x, y | H t ) = µ ( x, y ) ++dp ( xi , yi ) α ( xi , yi )( M i − M 0 )e( t − ti + c )
K 0 ( xi , yi )− q ( xi , yi )(2)を使う．ここで，一定の大きさ以上の地震については，震央座標 (x j , y j) や分散行列 S j は準リアルタイムで検出地震群（例えば 1 時間以内）から AIC で判定し補正される．

まず µ (x, y) が時間に依存しない常時活動度を永年確率として使用される．K 0 (x, y) は短期予測に有用な余震生産性のパラメータになり，本震の一時間後からの余震発生の時空間的予測に使用する．

－ 493 －

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

2021 年 3 月発行

本報告では，これまで CSEP 実験で予測した条件付強度の予測動画を示したが，第 2 図はそのスナップショットである．第 3 図は東北地方太平洋沖地震の 1 年程度前と後の条件つき強度関数や各時点後 1 年間の平均活動度の比である．このように地震活動の活発化（狭義・広義の余震活動）や静穏化（余震減衰や相対的静穏化）が可視化できる．

危険拡大度 λ (t , x, y | H t ) µ ( x, y ) は短期・中期の誘発（トリガリング）の時空間的効果を表す．常時活動度 µ (x, y) が高いところ6, 7, 13）に将来の M6.0 以上の大地震や過去の歴史被害地震の大地震の大方が発生している16) が，µ (x, y) が低いところで異常現象が認められた場合，その危険拡大度が高くなることは注目に値する．
3．地震検出率の時空間発展
今後の課題として，危険拡大度や背景率の解像度を上げるために，下限マグニチュードを十分下げてデータを増やしたい．その場合，長期的なデータの不均質性を補うと同時に，大きな地震直後の余震検出率の低下によるバイアスを補正して，正しい HIST-ETAS モデルの推定と予測が出来るようにしなければならない．そのために地震検出率モデル 21）を時空間的に展開する．このため震源から作ったデローネ 4 面体ネット上での時空間検出モデル

ξ (t , x, y, M ) = 10

a +b (t , x , y )( M − M 0 )q {M | µ '(t , x, y ), σ (t , x, y )}(3)を準備する．ここで b (t, x, y) は b 値の時空間変化，q {M | ⋅, ⋅} は正規分布の分布関数で，平均パラメー

タ µ ' (t, x, y) は 50% 検出のマグニチュードに対応し，σ (t, x, y) は部分的に検出されるマグニチュードの範囲を示す．ゆえに µ ' (t, x, y)＋3σ (t, x, y) が 99.5% 検出されているマグニチュードを与える．

本報告ではベイズ解 µ ' (t, x, y) の動画を示したが，東北沖地震直後のスナップショットを第 3 図に示した．これらを基に，低い下限マグニチュードで HIST-ETAS モデルの再推定が期待できる．
4．広域の地震活動のための短期・中期・長期モデル

主にマグニチュード系列データに基づいた前震の確率予測が，時空間 ETAS モデルによる誘発予測より確率利得や情報利得が高いという報告 8, 9 10）の根拠を解説した 11, 23)．また大地震発生直後（特に一日以内）は，余震検出率の低下によるバイアスの補正で，大きな余震のみならず，前震確率予測も可能 11, 12）で，地震検出率の時空間発展を組み込めば，時空間危険拡大度を評価できる．

時間または時空間 ETAS モデルの診断解析の異常に基づく中期予測に関して，危険拡大度の作成の試みを説明した．特に典型的な余震活動の解析の静穏化については時空間的な危険拡大度関数を確保することが可能である22)（第 4 図）．時空間 ETAS モデルで，検出率効果を補正して下限マグニチュードを下げて時空間の相対的静穏化 20）を詳細に可視化することも期待される．

長期予測の時空間的な危険拡大度は地震調査委員会の活断層群データのモンテカルロ法による地域別地震発生確率の算出方法 13）が有用である．また HIST-ETAS モデルの予測で宇津カタログ（明治大正＋歴史被害）を先駆情報として採用すると情報利得が有意に優れていることがわかる．

測地学的歪み変化の地域的ゆらぎ 17) は測定誤差だけでなく周辺地域の中小地震や遠地の大規模地震の地震時やその前後のすべり変動が混ざり合っている．これらの中から地震活動変化につながるモデル化・解析法 14) を開発したい．特に内陸部の地震活動について，稠密に網羅されているGEONET システムから追求したい．

－ 494 －

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

2021 年 3 月発行

その他，特定の物理量（歪み変化量 14)，地下電位など 18））と地震の発生強度との因果関係についての統計モデルは ETAS 効果と入力成分の和として情報利得や有意性を議論できるが，入力成分を積の形 19）で考慮できれば，尤度計算は数値的に複雑になるかも知れないが，時空間危険拡大度の多項目確率予測の枠組みに適している．

各危険拡大度の変化と多項目確率予測の適中率はたいして高くないことも見込まれる 15）が，何時何処でどのように低下・向上するか，変動経過を準リアルタイムでモニターして，大・中地震に至った場合の適中率や情報利得の地域性などを求め，今後の展開の礎としたい．
まとめ① 予め永年確率や各異常現象の確率を準リアルタイムで与える準備や環境を整えることが必要．
② 広領域での多項目確率予測式は点過程の条件付強度関数を使って，HIST-ETAS モデルの時空間背景率 µ (x, y, M) を永年強度として危険拡大度 (1) で時空間表現する．

③ 階層時空間 ETAS（HIST-ETAS) モデルは，標準的な短期予測に使用する．

④ 前震の統計的識別法は 独立 G-R 分布の ETAS モデルより優れた短期予測を与える．
⑤ 本震直後 ( 例えば 1 日以内）の余震予報は前震の確率予測も与える．

⑥ 相対的静穏化などの異常事象に基づき大中規模地震を予測する条件つき強度関数で時間的切迫度や影響地域範囲のスケーリングを示すモデルを備える．

⑦ 地震活動の異常強度の精度を上げるために，データのマグニチュード下限を下げ，欠測によるバイアスの補正ために時空間的検出率を準備する．

謝辞

気象庁震源カタログを使用しました．本研究は科研費＃17H00727 の補助を受けている．

文献1） 統計数理研究所（尾形）（2014），予知連会報 , 91 (12-6).

2） 統計数理研究所（尾形）（2016），予知連会報 , 96 (12-23).

3） 尾形 (2015). 統計数理 , 63 (1), https://www.ism.ac.jp/editsec/toukei/pdf/63-1-029.pdf
4） Ogata (2017), Ann. Rev. Earth Planet. Sci., 45 doi:10.1146/annurev-earth-063016-015918.
5） Vere-Jones (1978), J. Phys. Earth, 26, 126-146.

6） 統計数理研究所（尾形）（2008），予知連会報 , 79 (12-1).

7） 統計数理研究所（尾形）（2017），予知連会報 , 97 (1-3).

8） Ogata and Katsura (2014), J. Gophys. Res., 119 (11), 8457-8477.

9） 気象研究所（前田，弘瀬，溜渕）（2019），予知連会報 , 103 (12-7).

10）Lippiello et al. (2012), Scientific Reports, 2 (846).

11） 統計数理研究所（野村，尾形）（2019），予知連会報 , 103 (12-8).
12）Ogata and Omi (2020), Bull. Seism. Soc. Am., 110 (4), 1781–1798.
13）地震調査委員会 (2013), https://www.jishin.go.jp/main/chousa/13feb_chi_kyushu/kyushu_gaiyo.pdf
14）統計数理研究所（熊澤，尾形）(2019), 予知連会報 , 103 (12-12).
15）Cao and Aki (1983), J. Gophys. Res., 88 (B3), 2185-2190.

16）統計数理研究所（尾形）（2020），予知連会報 , 103 (12-13).

－ 495 －地震予知連絡会会報第 105 巻

2021 年 3 月発行

17）Wang et al. (2013), Geophys. Res. Lett., 40, doi:10.1002/grl.50554.
18）統計数理研究所（尾形，庄）（2001），予知連会報 , 66 (10-4).

19）岩田（2015），統計数理，63 (1), https://www.ism.ac.jp/editsec/toukei/pdf/63-1-129.pdf.
20）Ogata et al. (2003), Appl. Statist. (J. Roy. Statist. Soc. C), 52 (4), 499–509.
21）Ogata and Katsura (1996), Geophys. J. Int., 113, 727-738.
22）Ogata (2001), J. Geophys. Res. 106 (B5), 8729-8744.
23）Ogata and Katsura (2014), J. Geophys. Res., 119 (11), 8457-8477.

－ 496 －地震予知連絡会会報第 105 巻(a) GCMT M≧5.4
 , (b) JMA M≧5.0
 , 2021 年 3 月発行(c) JMA 内陸 M≧4.0
 , 地震

M≧4/deg^2 / day

時空間の µ パラメータ（通年確率の微分値）．
(a) Global CMT カタログの世界の µ 画像は，日本付近に限定しています．対象期間は 1973 年から 1995 年で，M ≥ 5.4．(b) 気象庁データからの推定イメージ，対象期間は 1926 年～ 1995 年で M ≥ 5.0，宇津のカタログを用いた前駆期間は 1885 年～ 1926 年．(c) 内陸部に限定した推計．対象期間は 1926 年～ 1995 年で M ≥ 4.0，宇津のカタログを用いて前駆期間は 1885 年～ 1926 年である．
Fig. 1 The µ parameter in space-time that is the derivatives of the perennial probability.	(a) World µ image from the Global CMT catalog is restricted in around Japan. Target period is 1973 - 1995 with M ≥5.4. (b) The estimate image from JMA data, Target period is 1926 - 1995 with M ≥ 5.0, and precursory period is 1885 1926 using Utsu’s catalog. (c) Estimates restricted to inland areas. Target period is 1926 - 1995 with M ≥ 4.0; precursory
period is 1885 - 1926 using Utsu’s catalog.
第1図Time : 2011.03.14 10:02:38.00

第2図HIST-ETAS モデルの予測動画の記載日時におけるスナップショット．一平方度，一日当たりの M ≥ 4 地震の発生率予測の対数が下記の色表で示されている．
Fig. 2 A snapshot of the HIST-ETAS model prediction movie at the indicated date and time, which is being tested in CSEPJapan. The logarithm of the predicted incidence of M≥4 earthquakes per square degree and per day is shown in the colortable below.

－ 497 －地震予知連絡会会報第 105 巻

2021 年 3 月発行

10(a)M ≥ 4.0

2

Ratio ofone-year
probabilities

3.

地震活動の増加・減少比

λ (t2 , x, y | H t )λ (t1 , x, y | H t )(b)

M ≥ 5.0

地震活動の増加・減少比

Ratio ofdaily
forecasts

113∫
∫

1

t2 + 365

t2
t1 + 365

t2 = 2012.05.01
t1 = 2009.11.01

t1

λ (t , x, y | H t ) dt
λ (t , x, y | H t ) dt

1 10
第3図

東北地方太平洋沖地震前後の強度の比．
(a) 2009 年 11 月 1 日 12 時時点の予測と 2012 年 5 月 1 日 12 時時点の予測の強度の予測強度の比．東北地方
太平洋沖地震の震源域とその周辺では，かなりの上昇率が見られる．なお，赤と青の小さな点は必ずしもノイズではなく，詳細なシグナルが多い．赤色は活発化で，東北沖地震の余震域と，その周辺部に，たくさん点々としている広義余震が見られる．紺色の地域は，以前に起きた大地震の余震活動の自然な減衰を示しているが，岩手宮城県境の余震活動や M9 地震の前震活動の，濃い青スポットは，いわゆる相対的静穏化を示している．（b）下限マグニチュードを M5 に高くとった，各時点後 1 年間の平均活動度の比．自然な活発化や静穏化が安定して見られる．
Fig. 3 The ratio of the forecasts before and after the Tohoku-Oki earthquake.	(a) The ratio of the forecasts at the indicated two time points, with considerable increase in the epicenter area of theMarch 11, 2011 M9 Tohoku-Oki earthquake and its vicinity. The small dots in red and blue are not always noises, but
possibly some signals. The red color is activation, and there are many red spots of wide-sense aftershocks in and aroundthe aftershock area of the Tohoku-Oki earthquake. Dark blue spots of aftershock activity on the Iwate-Miyagi border andforeshock activity of the M9 earthquake indicate so-called relative quietness, while darker areas indicate natural decay ofaftershock activity from previous large earthquakes. (b) Ratio of average activity for one year after each time point, where
the lower limit of magnitude was taken higher to M5. Regions of natural activation or the quiescence are seen stably.

Detection rates of M≧≧2 earthquakes

2011年3月11日本震直後

ｂ値

M≧2.0余震の検出率

第4図

時空間データの不均質性
等高線は，東北地方太平洋沖地震の直後のスナップショットから求めた b 値分布と，例えば M2 の地震が99.5％検出できている割合を示している．
Fig. 4 Heterogeneity of space-time data	The contours show the b-value distribution obtained from the immediate snapshot of the M9 Tohoku earthquake and thepercentage of, say, M2 earthquakes that have been detected 99.5% of the time.

－ 498 －

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

2021 年 3 月発行(c) Distance (in degrees) from the main shockto the subsequent earthquake(d) Risk Enhancement Factor of Relative Quiescence ∝λ0 (t − ti , | r − ri |),µ ( x, y )ri = ( xi , yi )where λ0 (t , r ) is smoothed Palm intensity on [0,8+)deg × [0, 20+) yrs .
調査した余震列の本震から後続する全ての M6.5 以上の地震との時空間の距離関係 22)．
(a) 余震列が変化点を持たない 43 個の本震から M ≧ 6.5 の大規模事象の時間差と距離のスタッキング図，横軸は与えられた距離を半径とする円盤の面積に比例するように，距離の二乗でスケーリングされている．
(b) 余震列が有意に相対的に静穏化 している 34 個の本震についての同様なスタッキング図．最初の 6 年程度の間に，距離の集中が顕著に見られる．(c) 調査した余震列の本震からの距離に対する単位面積 (1 平方度 )あたりの M ≧ 6.5 の地震の平均回数で，青と黄色のヒストグラムは最初の 6 年間，赤と紫色のヒストグラ
ムは 7 年から 20 年までの間の後続 M ≧ 6.5 地震平均数をそれぞれ表している．静穏化ありの場合，6 年以内かつ近傍範囲で，他のケースより 2 倍ほど後続地震が多く，特に 200km 以内で数倍の確率利得が見込まれる．(d) 相対的静穏化の後続地震についてのスケーリングを示す時空間危険度拡大関数は（b）の点パターンを非定常ポアソン時空間過程としてベイズ平滑化したものを使う．
Fig. 5 Distance in time and space from the main shock of the studied aftershock sequence to all subsequent earthquakes of M ≥6.522).	(a) Stacked diagram of time differences and distances between the 43 main shocks with no aftershock sequences withno change points and large events with M≥6.5, where the horizontal axis scaled by the square of the distance so thatthe horizontal axis is proportional to the area of the disk with the radius of the given distance. (b) Similar stacking plots
for 34 main shocks for which the aftershock sequence is significantly quiet relative to the ETAS prediction. There isa significant concentration of distances during the first 6 years or so. (c) Average number of M ≥ 6.5 earthquakes perunit area (1 square degree) against the distance from the main shocks of the studied aftershock sequence, is shown by
blue and yellow color histograms representing the first 6 years period, and red and purple color histograms representing
aftershock activity for the next 7 to 20 years, respectively. In the case of quieting, there are approximately twice as many
subsequent earthquakes within 6 years and in the vicinity as in the other cases, and the probability gain is expected to be
several times greater, especially within 200 km. (d) Spatio-temporal risk magnification function showing the scaling ofrelative quieting for subsequent earthquakes. This risk enhancement factor in time and space is obtained by the Bayesiansmoothing of the point pattern in (b) assuming Poisson space-time process.
第5図－ 499 －