緊急地震速報のための同時多発地震を識別する震源推定手法
気象庁地震火山部*

溜渕 功史

京都大学防災研究所**

山田 真澄

カリフォルニア工科大学***

Wu Stephen

A new approach to Identify Multiple Concurrent Events for Improvement ofEarthquake Early Warning
Koji Tamaribuchi
Seismological and Volcanological Department, Japan Meteorological Agency,
Otemachi 1-3-4, Chiyoda-ku, Tokyo, 100-8122, Japan

Masumi Yamada
Disaster Prevention Research Institute, Kyoto Univ.,
Gokasho, Uji, Kyoto, 611-0011, Japan

Wu Stephen
California Institute of Technology,
1200 E. California Blvd., Pasadena, California, 91125, USA*〒100-8122

東京都千代田区大手町 1-3-4**〒611-0011

京都府宇治市五ヶ庄***

1200 E. California Blvd., Pasadena, California 91125, USA

1

Abstract2

The extremely active seismicity after the 2011 off the Pacific3

Coast of Tohoku Earthquake (Tohoku Earthquake) caused

4

serious issues for the Earthquake Early Warning (EEW)

5system in Japan. Multiple concurrent earthquakes generated

6

many false alarms due to the failure of correctly determining

7

the source location. In this study, we propose a method to8

distinguish multiple concurrent events for EEW. We used the9

outputs of a B-delta method and principal component analysis,

10

together with P-wave arrival times and amplitudes of11

triggered stations. The information of non-triggered stations

12

is also used to constrain the locations. We constructed a13

likelihood function using the source parameters and solved

14

this optimization problem using the particle filter. The

15

appropriate choice of a likelihood function improves the16

accuracy of the location determination, and the probabilistic

17

framework evaluates the uncertainty of the hypocenters. We

18applied this method to 72 events for which warnings were

19

issued after the Tohoku Earthquake. The performance of

20

EEW was greatly improved, reducing the false alarms from21

22 to 0. The percentage of inaccurate warnings was also22

reduced from 66% to 6%.

23
24

Key words: Earthquake Early Warning, Multiple earthquakes,

25

Bayesian estimation, particle filter

26

1

1

§1.

はじめに

2

緊急地震速報は，震源に近い観測点で P 波を検知して3から直ちに揺れを予測し，強い揺れが来る前に住民に伝4えるリアルタイム警報システムである．このような警報

5システムは，国や地域によって対象エリアや手法は異な

6

るが，日本のほか，メキシコ，台湾，トルコ，イタリア，

7

アメリカなどでも実用化や開発研究がなされている

8［Allen et al.（2009）］．気象庁は日本全国を対象として，9推定した気象庁マグニチュード（M）3.5 以上または予測

10

した震度が 3 以上の場合に発表する高度利用者向けの緊

11

急地震速報の運用を 2006 年 8 月から，予測震度 5 弱以上12の場合に予測震度 4 以上の地域に対して発表する一般向

13

けの緊急地震速報の運用を 2007 年 10 月からそれぞれ開14

始し，2007 年 12 月にはそれらを予報，警報と位置付け

15

た［Hoshiba et al.（2008）］．2011 年（平成 23 年）東北地

16

方太平洋沖地震では，震源に最も近い宮城県石巻市に設17置された地震計で P 波を検知してから 5.4 秒後に予報（第181 報）を，8.6 秒後の予報（第 4 報）の更新時に警報を発19表した．この警報は，テレビ，ラジオ，携帯電話等を通20

じて広く一般に伝達された．このとき，警報を発表して21から実際に震度 5 弱相当以上の揺れになるまでの時間は，

22

震源に最も近い宮城県石巻市でも 15 秒程度であったと23推定される［Hoshiba et al.（2011）］．

24

この緊急地震速報処理は，地震波を検知した観測点数

25

が 1〜数地点という限られた状況から震源を時々刻々と26推定するため，得られた観測点数の増加に応じて複数の27

手法を適用している．気象庁が用いている現行の震源決28

定手法には，Fig.1（a）に示すように，大きく分けて 3

29

つの処理がある．

30

1 つ目は，全国に整備された気象庁の多機能型地震計2

1

約 200 地点を用いるもので，地震波を検知した観測点（ト2

リガ観測点）が地震発生直後の 1〜2 地点しかない場合に3は，B-Δ 法［Odaka et al.（2003），束田・他（2004）］，主

4

成分分析法［横田（1985）］，及びボロノイ分割を利用し5たテリトリー法［Kamigaichi（2004）］を用いて震源を早6期に推定する．その後時間が経過し，トリガ観測点が 37〜5 地点になった場合はグリッドサーチ法［Kamigaichi

8（2004）］に切り替え，震源の推定精度を向上させる．

9

2 つ目は，Hi-net［Okada et al.（2004）］約 800 地点の10

データを用いた着未着法［Horiuchi et al.（2005）］の処理

11である．この手法は，2 地点以上の Hi-net 観測点がトリ12ガした場合に行うもので，トリガしていない観測点を P

13

波がまだ到達していない地点と仮定して，トリガ観測点

14のトリガ時刻とトリガしていない観測点の理論走時の両15方を満足させる震源を探索する．ただし，Hi-net 観測網

16が疎な海域や島嶼部においては，決定が困難な側面もあ

17

る．

18

3 つ目は，緊急地震速報以前から行われている手法で，

19

多機能型地震計と Hi-net 等のデータから P 相，S 相を

20

AR-AIC 法［横田・他（1981）］等を用いて自動検測し，

21

最小二乗法によって震源を決定する処理である．これは，

22

少なくとも 5 地点以上の観測値を用いて計算を行うため，

23

推定精度は高いものの，迅速性の観点からは前述の 2 手

24法には及ばない．

25

それぞれの手法の迅速性と震源推定精度の利点を相補26的に利用するため，上記 3 つの震源決定処理は独立に行

27

い，それぞれ計算結果が得られた場合は，後者の処理結28

果を優先して緊急地震速報に使用している．

29

また M の推定には，多機能型地震計を用いて，以下に30

示す P 波部分の 3 成分合成変位波形の最大振幅から推定
3

2する P 波 M（𝑀௉ ）と，S 波到達以降も含めた全相 M（𝑀ௗ ）を，

3

気象庁（2010），気象庁（2012）］．

1

理論走時に基づいて使い分けている［Kamigaichi（2004），

0.72𝑀௉ = log(𝐴ௗ ) + 1.2 log(𝑅) + 5.0 × 10ିସ 𝑅− 5.0 × 10ିଷ 𝑑𝑒𝑝+ 0.46,

0.87𝑀ௗ = log(𝐴ௗ )+ 1.0 log(𝑅) + 1.9 × 10ିଷ 𝑅− 5.0

（1）

（2）

× 10ିଷ 𝑑𝑒𝑝+ 0.98.

5𝐴ௗ は最大変位振幅（10µm），R は震源距離（km），dep は6

いて計算される M の平均値を用いる方法が一般的であ7

るが，緊急地震速報処理においては他の値から大きく外8れた値（外れ値）の影響を軽減するため，震源近傍の最9

大 5 地点の中央値（1 地点の場合はその値，2 地点の場合

10

はその平均値）を用いる．ただし，最大変位振幅が小さ11く，基準値に満たない場合は，Hi-net による最大速度振12幅を用いて M 推定を行う［Yamamoto et al.（2008）］．

4

震源の深さ（km）である．M の算出には，各観測点にお13

この緊急地震速報処理では，時間的に限られたデータ

14から震源を推定するため，複数の地震がほぼ同時刻に発15

生すると，それらを 1 つの地震と認識して，誤った緊急16

地震速報を発表してしまうことがある．現行の処理では，

17

複数の観測点のデータが同一の地震によるものか，異な18

る複数の地震によるものであるか判定（以下，
「同一地震

19

判定」という）を行う際，主に以下の 2 つの条件の論理

20

和で判定を行っている．一つ目は震源を利用するもので，

21

P 波理論走時と観測された検測時刻を比較し，その差が22

あらかじめ決められた閾値の範囲内であれば同一と判断

23

する．二つ目は震源に依存しないもので，観測点間の見24かけ速度（観測点間距離を検測時間差で除したもの）を25

計算し，P 波速度よりも同等か速ければ同一と判断する．
4

1

このように，現行の処理における同一地震判定は，トリ

2

ガ観測点の検測時刻のみを用いているため，ほぼ同時刻

3にトリガした観測点は基本的に同一地震によるものと判

4断されてしまう問題点がある．一方，着未着法はトリガ

5

していない観測点の情報も使用しており，トリガ観測点

6の検測時刻のみを用いる前述の手法に比べるとある程度

7

分離できる利点はあるものの，前述の通り観測網が疎な8領域では決定が困難となる課題が残る．

9

実際，2007 年 10 月の一般向け緊急地震速報運用開始

10からの 6 年間（2013 年 9 月まで）に発表した 138 回の警11報のうち，警報を発表したにも関わらず対象地域の全域

12

で震度 3 未満であった事例は 32 回に上り，その約 8 割（2513回）を占める要因は，複数の地震を分離できなかったた

14めと報告されている［気象庁（2013a），気象庁（2013b）］．

15

特に， 2013 年 8 月 8 日には和歌山県北部で発生した M2.3

16の地震と東南海沖に設置された海底地震計からのノイズ

17

を 1 つの地震と識別したために，地震の規模を過大に評

18

価し，近畿地方を中心に関東地方から九州地方にかけて

19

警報を発表した．このように，ほぼ同時に発生した複数20の地震（同時多発地震）の分離やノイズとの識別の高精

21度化は，緊急地震速報の精度向上に必要不可欠で，急務22の一つである．

23

このような問題の解決策の一つとして，Hoshiba（2013）

24

は波動場そのものを，データ同化を用いて把握し，境界25

積分方程式法を用いて伝播の予測することで，揺れを予26測する時間発展型予測手法を提案した．この手法は，震27源を推定する必要がなく，同時多発地震や震源断層の広28

がりにも対応できるメリットがある．一方で，波動場の29

把握には多数の観測点が必要で，観測網が疎な島嶼部に30は適用が難しい．また，S 波による揺れをとらえてから51予測するため現在の手法と比較して猶予時間が短くなる

2といったデメリットがある．

3

これとは別に，Liu and Yamada（2014）はベイズ推定

4を用いて，同時に発生した複数の地震をより適切に分離5する震源決定手法を提案した．この手法は振幅分布のみ

6から尤度関数を作成することで，走時データに頼らない

7

同時多発地震識別の可能性を示した．ただし，この手法

8は振幅分布しか利用しないので，海域における震源の安9定性や深さの推定精度には課題がある．

10

そこで本研究では，ベイズ推定が異なる複数のデータ11を容易に扱えることを利用し，Liu and Yamada（2014）

12の手法を発展させて震源推定をすることを試みた．尤度

13関数には最大振幅の他，検測時刻，トリガしていない観14測点の情報，B-Δ 法，主成分分析法の推定結果を使用す15ることで，従来気象庁が独立に行っていた処理を統合す16ることとした．また同一地震判定では，検測時刻だけで17はなく振幅も利用し，さらに震源の不確定性も考慮する

18

ことで，より厳密な識別を行う手法を開発した．これら

19を適用し，従来の緊急地震速報処理の迅速性を維持しつ

20

つ，より適切に同時多発地震を識別できる震源推定手法

21を開発したので報告する．

22
23

§2.

使用データ

24

本研究で対象とした事例は，Table 1 に示す 2011 年 325

月 1 日から 4 月 30 日までの間に警報を発表した 71 事例，26

及び 2013 年 8 月 8 日に警報を発表した 1 事例の計 72 事

27

例である．入力として使用したデータは，現在気象庁の28

緊急地震速報処理で利用されている現地観測点処理の結29

果（多機能型地震計約 200 地点［原田（2007），気象庁

30（2013c）］）である．現地観測点処理とは，各観測点の現6

1

地で個別に波形処理を行う仕組みのことで，以下の処理

2を行う．まず，STA（Short Term Average）/LTA（Long Term

3

Average）を常時監視し，STA/LTA が一定の基準未満のと4

きは 1 分に 1 回生存情報を送信する．STA/LTA が一定の5基準を超えると地震波を検知したと判断し，その時刻か

6

ら振幅の絶対値がノイズレベルを越えた時点まで遡って

7

P 波の検測時刻とする．次に，P 波検測時刻から最初の 2

8秒までの間に，B-Δ 法による震央距離，主成分分析法に9よる震源方位を算出し，3 成分合成の最大変位振幅，上10下成分の最大速度振幅等とともに，気象庁の中枢処理シ

11

ステムに処理結果を送信する［気象庁（2008）］．なお，

12

最大振幅は逐次更新されるため，トリガオフになるまで

13毎秒送信を続ける．

14
15

§3.

手法

16

3.1

アルゴリズムの概要

17

本手法は，前節で述べた現地観測点処理の結果を中枢18側で処理するものである．地震波を検知した地点数が 119〜数地点と少ない状況下では，震源を精度よく求めるこ

20とは基本的には困難である．そこで，グリッドサーチの21

ように広範囲に仮想震源を配置して，最も確からしい震22

源を確率論的に探索することを考える．本研究では，震23

源パラメータ𝜃として緯度，経度，深さを仮定する（す24

なわち，𝜃= {𝑙𝑎𝑡,𝑙𝑜𝑛, 𝑑𝑒𝑝}）．発震時刻（T0）及びマグニ

25

チュード（M）は，3.3.2 節と 3.3.3 節で詳しく述べるが，27

𝜃を仮定すれば観測値（検測時刻と最大振幅）から解析

28

タの推定には，異なる物理量のデータを容易に扱うこと29ができるベイズ推定を用いた．ベイズの定理は，時刻 t30までに得られている観測値を𝐃௧ とすると，

26

的に計算できるのでここでは仮定しない．震源パラメー

7

𝑃(𝜃|𝐃௧ ) =𝐿(𝐃௧ |𝜃)𝑃(𝜃),𝑃(𝐃௧ )（3）2で表される．ここで，𝑃(𝜃|𝐃௧ )は事後確率分布，𝐿(𝐃௧ |𝜃)は3

尤度関数，𝑃(𝜃)は事前確率分布である．分母𝑃(𝐃௧ )は観測値に関する確率密度分布であることから定数である．

4

つまりベイズ推定は，事前に得られる情報（事前確率分5布）と，観測値から定義される尤度関数を用いることで，

6

震源パラメータを事後確率分布として確率論的に求める

7

ことが可能である．

1

8

処理のフローを Fig. 1（b）に示す．はじめに震源を 19つと仮定して，ある時点までに得られた観測値の全てを10説明できる震源位置を探索する．観測値の中に推定した

11

震源では説明できないものがある場合には，震源の数を

12

1 つずつ増やして震源位置を推定する．得られた震源の13推定結果から，距離減衰式に基づいて M と予測震度を算14出し，その時点における推定値とする．これを 1 秒毎に15

繰り返し行い，得られた観測値が増加するに従って，震16源の推定精度を向上させる．
事前確率分布の設定

17

3.2

18

事前確率分布𝑃(𝜃)は，情報がない場合は一様分布によ19って定義するが，事前にある程度どのあたりに震源があ20るか分かっていれば，先見的な知見として与えることが21

望ましい．本研究では，過去の震源分布とボロノイ分割

22の考え方を基に事前確率分布を定義した．

23)
まず，過去の震源データを基にした事前確率分布𝑃
௦ (𝜃

24

は，1997 年 10 月から 2012 年 12 月までの一元化震源カ25

タログ（例えば，気象庁［2013d］）を用いて，緯度，経26度方向にそれぞれ 0.5°，鉛直方向に 50km のブロック内

27で震源の数をカウントし，震源パラメータ𝜃が含まれる

28

ブロックに𝑛௦ 個の震源が存在するとき，
8

0.01
(𝜃
)
𝑃
∝
൝
0.01𝑛
௦
௦
1

𝑖
𝑓0 = 𝑛௦𝑖𝑓0 < 𝑛௦ ≤ 100.
𝑖𝑓100 < 𝑛௦（4）1と定義した．Fig.2（a）に宮城県周辺の例を示す．群発2

地震や余震等の影響で局所的に確率が上昇することがあ34

るので，𝑛௦ が 100 以上の場合には 1 に固定した．なお，

5合から，東経 90 度〜180 度，北緯 0 度〜70 度，深さ 0km

6〜900km とした．

本手法における仮想震源の最大探索範囲は，計算上の都

7次に，ボロノイ分割は，ある地点がどの観測点に最も8近いかによって領域分けされる分割方法である．具体的

9には，最も早くトリガした地点と，他の観測点との垂直10

2 等分線で囲まれた領域を指し，例を Fig. 2（b）に破線11で示す．しかし，特に東北地方太平洋沖地震以降は，余12

震活動によってノイズレベルが上昇したために，必ずし13も震源に最も近い観測点がトリガしないことがあったの

14)を以下のようで，ボロノイ分割を拡張して事前確率𝑃
௩ (𝜃

15

に定義した．

) ∝ exp ቆ−𝑃௩ (𝜃(𝑘− 1)ଶ
ଶ
2𝑐
ଵ

ቇ.

（5）

16

ここで，c1 はボロノイ分割をどのくらい拡張するかを規17定するパラメータで，本研究では c1 = 10 とし，Fig.2（b）

18に例を示す．

19（4）式と（5）式の積により，事前確率分布𝑃(𝜃)は𝑃(𝜃) =

)𝑃
)
𝑃
௦ (𝜃
௩ (𝜃

)𝑃
)𝑑𝜃
∫ఏ 𝑃
௦ (𝜃
௩ (𝜃

.

（6）

20

で表わされる．（6）式によって定義された事前確率分布21の例を Fig. 2（c）に示す．
尤度関数の設定

22

3.3

23

本研究では，震源の推定に B-Δ 法による震央距離，主24

成分分析法による震源方位，検測時刻，最大振幅の各観91

測値を尤度関数に用いた．

2

まず，それぞれの観測値について尤度関数を定義する．

3

それぞれの観測値と震源パラメータ𝜃における期待値と4の差が正規分布で近似できると仮定して，尤度関数を設5定した．観測値を𝑥，期待値を𝜇，モデルの推定誤差を𝜎で6表すとき，正規分布 𝑓(𝑥|𝜇, 𝜎)は以下のように表わされる．
𝑓(𝑥|𝜇, 𝜎) = exp ቆ−(𝑥− 𝜇)ଶ
ቇ.
2𝜎ଶ

（7）

7

なお，（6）式において総和が 1 になるように規格化する

8

ため，（7）式以降では定数である規格化因子は省略して9いる．このとき，他の地震やノイズによるデータが混入

10

すると，そのデータによって尤度が大きく影響されて，

11

適切な震源パラメータの推定が困難となるため，実用上12問題となる．この影響を軽減するために，モデルから大13

きくかけ離れた観測値（いわゆる「異常値」）に対して下14

限値（𝑔଴ ）を設定した．遠方の観測点が尤度に与える影

15
16
17

18

響を小さくするため，下限値は震源から遠くなるほど大きくなるように定義し，下限値（𝑔଴ ）を設定した正規分

布を変形した確率密度関数𝑔(𝑥|𝜇, 𝜎, 𝑘௜ )を以下の式で表す．
𝑔(𝑥|𝜇, 𝜎, 𝑘௜ ) = ൫1 − 𝑔଴ (𝑘௜ )൯𝑓(𝑥|𝜇, 𝜎) + 𝑔଴ (𝑘௜ ).

（8）

ただし，

(𝑘௜ − 1)ଶ1
(𝑘)𝑔଴ ௜ = 1 − exp ቆ−
ቇ.
ଶ
2
2𝑐
ଶ

（9）

20

𝑘௜ は，観測点 i が震源から𝑘௜ 番目に近い観測点であること

21

影響を及ぼすかを規定するパラメータである．正確な震22源決定には複数の観測値が必要であるが，震源から遠く23離れた地点の観測値が震源推定に影響を与えないように24するため，本研究では c2=10 とした．Fig. 3 に𝑘௜ を変えた

1925を表す．c2 はどの程度遠方の観測点まで確率密度関数にときの下限値𝑔଴ と，確率密度関数𝑔の例を示す．ここで10

1

定義した確率密度関数を用いて，各観測値に対する尤度

23

関数𝐿௕ ，𝐿௧ ，𝐿ெ を定義する．𝐿௕ ，𝐿௧ ，𝐿ெ はそれぞれ震央

4

ある．

5
6

距離と震源方位，検測時刻，及び最大振幅の尤度関数で本研究では簡単のため，各観測値を独立とみなし，求める尤度関数𝐿(𝐃௧ |𝜃)を以下の式で表す．

𝐿(𝐃௧ |𝜃) = 𝐿௕ (𝐃௧ |𝜃) ∙ 𝐿௧ (𝐃௧ |𝜃) ∙ 𝐿ெ (𝐃௧ |𝜃).

（10）

7

Fig. 4 にそれぞれの尤度分布の例を示す．以降の節では，

8それぞれの観測値の尤度関数について具体的に述べる．
震央距離及び震源方位

9

3.3.1

10

観測点 i において B-Δ 法及び主成分分析から得られる

11෨௱௜ ， 𝑅
෨ ௔௜ とす
震央距離及び震源方位の推定値を，それぞれ𝑅

12
13
14

る．また，震源パラメータ𝜃について，観測点 i との震央距離及び震源方位を𝑅௱௜ ，𝑅௔௜ とすると，震央距離及び震源方位の尤度関数𝐿௕ は，
ே

15
16
17

෨௱୧ |𝑅௱௜ , 𝜎෨ ௔௜ |𝑅௔௜ , 𝜎௔ , 𝑘௜ ൯,𝐿௕ (𝐃௧ |𝜃) = ෑ 𝑔൫𝑅൫𝑅௱, 𝑘௜൯ ∙ 𝑔（11）௜ୀଵで定義した．ここで𝜎௱ ，𝜎௔ はそれぞれ B-Δ 法，主成分分析法の推定誤差で，気象庁（2011）を参考に𝜎௱ =100km，𝜎௔ =20°を与えた．N は全観測点の数で，データが得られ19ていない観測点については𝑔(𝑥|𝜇, 𝜎, 𝑘௜ ) = 𝑔଴ (𝑘௜ )を与えた．

20

検測時刻の尤度関数を定義するために，以下のように21
22

して発震時刻𝑇を用い
଴ を求めた．まず，震源パラメータ𝜃

23

観測点の検測時刻からこれを引くことで，発震時刻の初24期値を得る．得られた発震時刻の初期値と𝜃から各観測

25

点の P 波到達予想時刻を求め，各観測点で予想時刻の±26

3 秒の範囲内にある検測時刻を選択した．選択した検測27時刻から P 波理論走時を引くことで，各観測点で発震時

18

3.3.2

検測時刻

て，1st トリガ観測点の P 波理論走時を求める．1st トリガ

11

1
2
3
4

5

刻の推定値が得られ，その平均値を発震時刻𝑇
଴ とした．

次に検測時刻の尤度関数を定義する．𝑇
から観測点
଴ と𝜃

௜
i における P 波到達予想時刻𝑇
௖௔௟ を得られるので，検測時
௜
刻𝑇
௧ は，
௧௥௚ の尤度関数𝐿

ただし，

ே

𝐿௧ (𝐃௧ |𝜃) = ෑ 𝐿௧ ൫𝐷௧௜ |𝜃൯.

（12）

௜ୀଵ

௜
௜
௜
𝑔൫𝑇
௧, 𝑘
௜൯ 𝑖
𝑓𝑇
𝑥𝑖
𝑠
𝑡
𝑠
,
௧௥௚ |𝑇௖௔௟ , 𝜎௧௥௚ 𝑒௜𝐿௧ ൫𝐷௧௜ |𝜃൯ = ቐ𝑔൫𝑇(∗)௡௢௪ |𝑇௧, 𝑘௜൯
௖௔௟ , 𝜎
𝑜
𝑡
ℎ𝑒
𝑟
𝑤𝑖
𝑠
𝑒,
𝑔଴ (𝑘௜ )

（13）

௜
௜
(∗) 𝑇
𝑜𝑒𝑠𝑛𝑜𝑡𝑒𝑥𝑖𝑠𝑡,𝑇𝑛𝑑𝑘௜ ≤ 5.
௡௢௪ − 𝑇௧௥௚ 𝑑௖௔௟ ≥ 1, 𝑎7と定義した． 𝜎௧ は走時残差の標準偏差で経験的に𝜎௧ =1.5

8

ム上の時刻（現在時刻）である．（13）式の（*）は，着69を与えた．𝑇௡௢௪ は処理を行っている時点におけるシステ
未着法の考え方を基に，震源近傍 5 地点以内で，かつ𝑇
௡௢௪

11

௜
が P 波到達予想時刻𝑇
௖௔௟ を 1 秒以上経過してもトリガし12௜
ていない観測点については，その残差（𝑇௡௢௪ − 𝑇௖௔௟ ）をペナルティとして尤度関数に与えることを意味する．そ13

れ以遠の観測点で検測時刻が得られていない場合は，下14

限値𝑔଴ を与えた．

10

15

3.3.3

最大振幅

16

最大振幅の尤度関数は，各観測点において算出された

17

M の中央値からの偏差を用いて定義する．震源パラメー

1819

タ𝜃のとき，観測点 i における速度マグニチュード𝑀௩௜ は20

最大速度振幅𝐴௜௩ を用いて以下のように求められる［舟崎・地震予知情報課（2004）］．

21

𝑀௩௜ = 𝛼log൫𝐴௜௩ ൯ + 𝛽௩ ൫𝑅௱௜ , 𝑑𝑒𝑝൯ + 𝐶௩ .

（14）

22

𝛼は定数，𝛽௩ は距離減衰項，𝐶௩ は観測点補正項である．

23

る M のばらつきを示すように，以下の式で定義した．

このとき，最大振幅の尤度関数𝐿ெ は，各観測点におけ
12ே෡௩ , 𝜎ெ , 𝑘𝐿ெ (𝐃௧ |𝜃) = ෑ 𝑔൫𝑀௩௜ |𝑀
௝ ൯.

（15）

௜ୀଵ

2

ここで，𝜎ெ は M の標準偏差で，経験的に𝜎ெ =0.7 を与え

3

෡௩ は近傍 5 地点における𝑀௩௜ の中央値である．この尤た．𝑀
度関数の M 計算では，各観測点における振幅の条件を揃4えるため，検測時刻から 5 秒以内に得られた最大振幅を5用いて，速度 M を計算した．速度 M を用いる理由は，

6

対象とする地震が M4 未満であっても SN 比がよく，ま7た変位振幅に比べて時間的に減衰が早く，異なる地震の8

後続波の影響を比較的受けにくい点である．速度センサ

9は振り切れる可能性があるが，情報発表に用いる M とは

10異なり，尤度計算では各観測点における M の偏差のみを11計算できればよく，さらに検測時刻から 5 秒以内に得ら12れた最大振幅を用いるので振り切れの問題はほぼ無視で13きる．ただし，速度センサが併設されていない観測点な

14

どで，最大速度振幅が得られなかった場合は，最大変位15

振幅から（2）式を用いて観測点 M を推定した．

1

パーティクルフィルタ法

16

3.4

17

リアルタイムに震源の推定を行うためには，処理の効18率化も大きな課題の一つである．グリッドサーチを用い

19ても前節までで述べた事後確率を推定することは可能で20あるが，本研究では処理をより効率的に行うことができ21

るパーティクルフィルタ法［例えば，Arulampalam et al.

22（2002）］を用いた．パーティクルフィルタ法は，モンテ

23

カルロ法の一種で，事前確率分布に従ってランダムに配24置した仮想震源について，観測値に基づいて尤度を算出

25

し，尤度の高い仮想震源の周辺に重点的に仮想震源を再26

配置する手法である．具体的には，以下の 4 つのステッ

27

プの繰り返すことで，震源の推定を行う（ただし Step 1

28

は 1 回目の推定のみ）．
13

1

Step 1. 初期化：事前確率分布𝑃(𝜃)（（6）式）に従って

2

仮想震源を配置する．本研究では，仮想震源の数は 1000

3

とした．

4

Step 2. 尤度計算：それぞれの仮想震源において，現時56点までに得られた観測値を用いて尤度𝐿(𝐃௧ |𝜃)（（10）式）

7

Step 3. 最尤値推定：最も尤度が高い仮想震源を時刻 tを計算する．

9

෠とする．
における震源の推定値𝜃

10の尤度を規格化し重み w を求める．つまり，配置した震11

源分布を𝜃(௝) （𝑗= 1, … , 𝑚; 𝑚 = 1000）とすると𝑤(௝) は，

8Step 4. リサンプリング：Step 2 で得られた各仮想震源(௝)

𝑤

𝐿൫𝐃௧ |𝜃(௝) ൯= ௠,∑௝ୀଵ 𝐿(𝐃௧ |𝜃(௝) )

（16）

12

である．得られた重み分布を基に仮想震源を復元抽出し13た．ここで抽出された仮想震源は，次の時間ステップに14

おいて利用されるが，復元抽出であるため，ある特定の15仮想震源に収束してしまう危険性がある．そのため，抽16出された仮想震源に，わずかに摂動を与えた．具体的に

17は，水平方向に 0.1 度，深さ方向に 20km の標準偏差を18持つ正規乱数を摂動として与えた．次の時間ステップに

19

おいては，ここで求めた仮想震源を利用し，Step 2 以降

20を繰り返す．

21

෠ を，時間の経過ととこのようにして，震源の推定値𝜃

22

もに繰り返し更新する．Step 3 での震源の推定値として，

23

重み付き平均ではなく最も尤度が高い仮想震源を用いて24

いるのは，観測値が少ない初期の段階で，事前確率に大25きく左右されることを防ぐためである．深さについては，

26

トリガした観測点数が 3 地点以下の場合，または最も近27

いトリガ観測点からの震央距離が 150km 以上離れてい28る場合は，深さの精度がないため，過去の震源分布を考14

1

慮して 10km に固定した．

2

情報発表に利用する M は，現行の緊急地震速報処理と3同様にトリガ時刻以降の全ての観測値の最大変位振幅か4らそれぞれ M を算出し，近傍 5 地点の中央値を採用した．

5

また，M の算出においては，他の地震やノイズの影響を6軽減するため，震央距離が「トリガした観測点と震源の7

最短距離＋150km」よりも遠い地点の M は用いないこと

8

とした．

9

3.5

10

複数地震を識別するために，時刻 t，観測点 i における

11

観測値𝐷௧௜ がパーティクルフィルタ処理によって得られ12
13

複数地震の識別（同一地震判定）

た重み分布によって説明できるか判定する．震源パラメータ𝜃の下で観測値が𝐷௧௜ である確率𝑝൫𝐷௧௜ |𝜃൯は，
𝑝൫𝐷௧௜ |𝜃൯=

141

ඥ2𝜋𝜎௧ଶ

exp ൝−ଶ௜௜൫𝑇
௧௥௚ − 𝑇
௖௔௟ ൯

2𝜎௧ଶ

ଶ

෡௩ ൯
൫𝑀௩௜ − 𝑀exp ൝−
ൡ
ൡ
ଶ
ଶ
2𝜎ெ
ඥ2𝜋𝜎ெ
1

ଶ

௜
௜
෡௩ ൯ଶ
1
1 ൫𝑇
൫𝑀௩௜ − 𝑀௧௥௚ − 𝑇௖௔௟ ൯=
exp ൥− ൝+ൡ൩ , (17)ଶ
22𝜋𝜎௧ 𝜎ெ𝜎௧ଶ
𝜎ெで表わされる．ここで，同一の地震と判断する条件を，
1
𝐴ଶ
exp ቆ− ቇ,
2𝜋𝜎௧ 𝜎ெ2𝑝൫𝐷௧௜ |𝜃൯ ≥

（18）

15

と設定する．A は閾値で，本研究では A=3 を与えた．つ16

まり，（17）式，（18）式から，震源パラメータ𝜃と観測

17

値𝐷௧௜ が同一の地震と判断する条件は，
ଶ

௜
௜
൫𝑇
௧௥௚ − 𝑇
௖௔௟ ൯

𝜎௧ଶ

+

ଶ

෡௩ ൯
൫𝑀௩௜ − 𝑀≤ 𝐴ଶ ,
ଶ
𝜎ெ

（19）

18

を満たせばよい．これは±3𝜎程度の誤差楕円の範囲内で19あれば，同一地震による観測値とみなすことを意味する．

20

このとき，重み分布と観測値𝐷௧௜ との同一地震判定は，そ21

れぞれの震源パラメータの重み付き和によって表わされ，15௠෍ 𝑤(௝) 𝑝൫𝐷௧௜ |𝜃(௝) ൯ ≥

1
2

௝ୀଵ

1
𝐴ଶ
exp ቆ− ቇ,
2𝜋𝜎௧ 𝜎ெ2（20）である．𝜃(௝) はパーティクルフィルタ処理によって配置

した震源パラメータ，𝑤(௝) は𝜃(௝) における重みである．従3

来手法における同一地震判定は，主に各観測点の検測時

4

刻のみに頼っていたが，本手法では振幅値（M 残差）と5

検測時刻（走時残差）が含まれるため，従来手法に比べ

6て複数イベントをより適切に分離することが可能となる．

7

これに加え，重み分布を利用することで，震源の不確定

8性も考慮できる．

9

複数の地震が発生した際に識別する方法について，以10下に示す．まず，ある観測点𝑠ଵ が最初に地震波を検知し11stたとき，観測点𝑠
ଵ を 1 トリガ観測点とする．そこで観測

13

点𝑠
ଵ のデータを満たす 1 つ目の地震について 3.4 節で述

14

点が地震波を検知したとき，1 つ目の地震の重み分布と15

検測時刻及び最大振幅を比較し，1 つ目の地震による揺16

れとみなせるか判定を行う．この判定を時間ステップご

17とに繰り返し，判定の結果同一とみなせない観測点𝑠ᇱ が12

18

べたパーティクルフィルタ処理を行う．次に，他の観測存在すると，2 つ目の地震が発生したと判断して，その20

ᇱ
観測点を 1st トリガ観測点𝑠ଵ として新たにパーティクルフ

21の識別を行う．新しい観測値が複数の地震のどちらでも

22説明できる場合は，両方の尤度関数に観測値を利用する．

23

また逆に，1 つの地震をノイズ等の影響で誤って複数24に分離してしまった場合でも，後から 1 つの地震に統合

25

できるように考慮しておく必要がある．そこで，まず分26離された 2 つの地震で観測値が（20）式の条件を満たす

27

かどうか確認する．それぞれの地震で（20）式の条件を28満たす観測値のうち，半数以上が両方の地震で（20）式

19ィルタ処理を行う．以下同様にして，3 つ目以降の地震

161を満たした場合は，2 つの地震を 1 つの地震に統合させ2ることとした．例えば，1 つ目の地震に対して（20）式3の条件を満たす観測値𝐷௧ଵ ，𝐷௧ଶ ，𝐷௧ଷ ，𝐷௧ସ ，𝐷௧ହ があり，25

つ目の地震に対して条件を満たす観測値𝐷௧ସ ，𝐷௧ହ ，𝐷௧଺ があ

6式の条件を満たしており，かつ 2 つ目の地震は 3 観測値4ったとする．ここで，観測値𝐷௧ସ ，𝐷௧ହ は両方の地震で（20）8（𝐷௧ସ ，𝐷௧ହ ，𝐷௧଺ ）のうち 2 観測値が重複しているので統合

9

最適な仮想震源の尤度を比較し，大きい尤度を持つ地震

10（つまり（20）式を満たす観測値が多い方の地震）を採11用する．小さい尤度を持つ地震については，パーティク

12

ルフィルタ処理をこれ以上継続させないこととした．

7

する対象となる．なお，統合させる際には，それぞれの

13
14

§4.

結果と議論

15

4.1

単独地震の事例（平成 23 年（2011 年）東北地方

16

太平洋沖地震）

17

まず，地震が単独で発生した場合の事例（2011 年 3 月

1811 日 14 時 46 分に発生した平成 23 年（2011 年）東北地

19

方太平洋沖地震）を示す．従来処理では，震源に最も近20い石巻大瓜（OURI）観測点で最初にトリガ検知してから21

8.6 秒後に M7.2（最大震度 5 弱）を予測し，警報を発表

22

した．その後，時間の経過とともに徐々に M の推定値が23大きくなり，最終的にはトリガ検知から 105.0 秒後に

24

M8.1（最大震度 6 弱）を予測した．

25

本手法の処理では，震源に最も近い観測点で最初にト26リガ検知してから 2.1 秒後に事前確率分布に従って仮想

27震源を配置した．Fig. 5（a）に，このときの仮想震源分

28

布を示す．仮想震源の色は重みを示す．その後，石巻大29

瓜観測点でトリガ検知してから 9.1 秒後までに受信した

30

データを用いて，M7.4（最大震度 5 弱）を推定した（Fig.
1715（b））．これは，従来処理と比較しても，警報発表時間

2の遅延はほとんどない．また，M，震源の位置，及び深3さの比較を Fig. 5（c），（d），（e）に示すが，従来の手法

4に比べて震源決定精度に大きな差は確認されなかった．
複数地震の事例（2011 年 3 月 19 日 6 時 18 分の

5

4.2

6

地震）

7

次に，複数の地震がほぼ同時刻に発生した場合の事例

8を示す．2011 年 3 月 19 日の事例では，6 時 18 分 48 秒に

9

栃木県北部で M2.5（最大震度 1）の地震が発生し，5 秒

10

後の 6 時 18 分 53 秒に茨城県沖で M4.8（最大震度 2）の11

地震が発生した．従来処理では，栃木県北部の地震の理

12

論走時と，茨城県沖の地震による検測時刻が重なったた13

め，これらの地震を 1 つの地震と識別した．その結果，214

つの地震の振幅を混在して M を過大に算出（M5.8，最大15

震度 5 弱）したために，警報を発表した［気象庁（2012）］．

16

本手法では，最初の地震波を検知してから 2.1 秒後に

17

1 つ目の地震の処理を開始した（Fig. 6（a））．その後，

18

最 初 の ト リ ガ 検 知 か ら 14.2 秒 後 に ， 銚 子 天 王 台

19（TENNOD）観測点で地震波を検知した．ここで同一地

20

震判定の結果，このデータが 1 つ目の地震によるもので21はないと判断し， 2 つの地震に適切に分離した（Fig. 622（b））．2 つ目の地震である茨城県沖の地震（M4.8）の推23

定結果の比較を Fig. 6（c）〜（e）に示す．本手法では初24期の段階から震源を 2 つに分離できたため，茨城県沖の25

地震の M の推定値は M4.6〜M4.8 程度と，誤差が 0.2 以

26

下に減少した（Fig. 6（c））．一元化震源の位置と比較す

27

る．Fig. 6（d）に示すように，従来の手法では 1 つ目の28

地震波を検知してから 30 秒後に，一元化震源と比べて水29平距離で 150km 以上離れた位置に震源を決定した．しか30

し本手法ではその水平距離は 50km 程度と大きく改善さ
18

1

れた．深さについても，本手法では 30km 程度で，一元2化震源の深さ 44km と概ね整合している（Fig. 6（e））．
地震とノイズが重なった事例（2013 年 8 月 8 日

3

4.3

4

の地震）

5

2013 年 8 月 8 日 16 時 56 分に和歌山県北部で発生した

6

M2.3（震度 1 以上の観測はなし）の地震の事例では，ト7リガ検知から 17.4 秒後に東南海 OBS のノイズを 1 つの8

地震と識別し，OBS のノイズを用いて M を計算した結9果，M7.8（最大震度 7）を予測し，警報を発表した．こ

10のノイズの直接的な原因は陸上中継局の不具合によるも11のであったが，雷ノイズのように，地震計は誤ったトリ

12ガを送信する可能性を常に持っている．

13

本手法では，Fig. 7（a）に示すように，和歌山県北部

14の地震によるトリガを検知した直後から震源の推定を開15始した．その後，最初のトリガ検知から 17.4 秒後に東南

16

海 OBS3（TT3OBS）観測点でトリガを検知した．このと17

き，東南海 OBS3 の観測値は検測時刻及び最大振幅から

18

和歌山県北部の地震によるものではないと判断し，2 つ19目の地震として処理を開始した（Fig. 7（b））．この 2 つ

20目の地震は，それ以降しばらく時間が経過しても同一と21判定される他観測点のデータが得られなかったため，約

22

60 秒後に処理を終了した．2 つ目の地震について，仮に

23

過大な振幅を用いて M を計算しても，警報を発表するに

24

は 2 地点以上での検知が必要なので，警報発表には至ら

25

ない．和歌山県北部の地震における M は，初期の段階か

26

ら安定して M2 程度（ただし速度 M）を推定し，警報を27

発表することはない（Fig. 7（c））．震源の位置及び深さ

28についても，ノイズを適切に分離したために一元化震源29と比較して大きく誤って決定することはなかった（Fig. 730（d），（e））．このように，ノイズが混在した場合でも，
19

1

最大振幅値を活用することで適切に分離することが可能

2

である．

3

4.4

4

本研究で対象とした 72 事例について，本処理で推定さ5れた警報対象地域と実際に観測された震度 3 以上の地域

6を比較し，警報対象地域の全域で，震度 3 以上が観測さ7れた場合を適切（Accurate），警報対象地域の 1 つでも震8度 2 以下であった場合を不適切（Inaccurate）と分類した．

9

さらに，不適切な事例のうち，警報対象地域の全てで最10

大震度 2 以下であった事例を誤報（False）と分類した．

11

これらの分類結果を Fig. 8 に示す．従来の処理で誤報に12分類された事例は 22 事例であったが，本手法では 0 事例

13であった．また，不適切な警報を発表する割合は，従来14の処理では 66%（48/72 事例）であった一方，本手法で15は 6%（4/72 事例）であった．従来の手法と本手法によ16

る最大予測震度と観測された最大震度の階級差を Fig. 917に示す．従来の手法において過大に予測した事例が軽減18されたことが分かる．また，過少に予測した事例もわず

19

かながら改善した．なお，震度階級差は，震度 5 弱と 5

20

強，6 弱と 6 強でも 1 階級の差と定義しているため，震21

度 0（無感）の地震を震度 7 と予測した場合は，9 階級異22なることになる．

全事例のまとめ

23

本手法において，警報対象地域の一部で震度 2 以下だ

24ったため不適切と分類された 4 事例のうち 2 事例（201125

年 4 月 14 日 20 時 24 分の福島県浜通りの地震と 2011 年

26

4 月 19 日 4 時 14 分の秋田県南部の地震）は，震源直上27の観測点で M を過大に推定し，M 推定の誤差が大きくな

28

った．震源近傍における距離減衰式や，震度推定の誤差

29の改善は今後の課題である．その他の 2 事例のうち，201130

年 3 月 12 日 3 時 11 分の福島県沖の地震では，数秒前に20

1

ほぼ同じ場所で別の小さな地震が発生し，地動レベルが2上昇したため，震源近傍の観測点でトリガしなかった．

3

本手法では，トリガしていない観測点を「まだ揺れてい4ない」と仮定して処理を行っているため，このような場

5合は不適切な警報を発表する可能性がある．また，20116

年 3 月 19 日 18 時 56 分の茨城県北部（M4.0，最大震度 4）7の事例では，約 30 秒後にほぼ同じ場所で M6.1 の地震が8発生した．これらの震源は適切に分離できたものの，前9者の地震の M 推定において後者の地震による振幅が混

10

入し，M を過大に推定した．このように，震源近傍の観11測点であるにも関わらずトリガしない場合や，時空間的12に非常に近接して地震が発生した場合には，震源や M を13

適切に推定することは現時点では困難であり，引き続き

14

改善に向けた検討を進めたい．

15

逆に，適用した事例のうち，震度 5 弱以上を観測した

16にも関わらず警報発表基準に満たなかった事例（見逃し17事例）は 4 事例だった．そのうち震度 2 階級以上の差が18あった 2 事例について，詳細を以下に述べる．1 つ目の19

事例は 2011 年 3 月 12 日 5 時 42 分の新潟県の地震（M5.3，

20

最大震度 6 弱）である．この事例では，震源直上の観測

21点で M の割に比較的大きい震度（6 弱）を観測した．本22

手法による推定 M は M5.1（予測最大震度 4）で，マグニ

23

チュードは適切に推定できていたが，震度予測式の誤差24が大きく，適切な震度推定ができなかった．もう 1 つの25

事例は，2011 年 3 月 23 日 7 時 36 分の福島県の地震（M5.8，

26

最大震度 5 強）である．この事例は，2 分前の 7 時 34 分

27に福島県で発生した M5.5（最大震度 5 強，現行処理では

28

警報未発表）の地震が発生し，7 時 34 分の地震について

29

は本処理では適切に警報を発表したが，地動レベルが上30昇したことによって，7 時 36 分の地震では 1 地点しかト
21

1

リガデータが得られなかった．前述と同様に，大規模な

2

地震発生直後に近傍の別地震でトリガ検知することは現

3行のシステムでは困難で，今後の課題である．ただし，

4

このような見逃し事例は，Hoshiba（2013）のような時間

5

発展型予測を併用することで，改善が期待できると考え

6

られる．

7

本手法で特筆すべき点は，既存の緊急地震速報処理で8使用している観測値と全く同じものを使用している点で9ある．つまり，既存の枠組みの範囲内でも緊急地震速報

10の震源推定精度を大きく向上させることが可能である．

11

また，本処理の処理にかかる時間は，簡易的な並列化処12理である OpenMP を用いた 12 スレッド並列処理で，1 イ

13

ベントあたり 0.3~0.5 秒程度である（動作環境は，CPU:

14

Xeon 2.26GHz，OS: Linux 2.6.32，メモリ:32 GB，コンパ

15

イラ:gcc ver.4.4.7）．すなわち，2 つ地震が発生した場合

16でも，約 1 秒ごとに震源情報を更新でき，緊急地震速報

17の用途として十分に実時間で処理を行うことが可能であ

18

る．

19

なお，今後の展望として，緊急地震速報の迅速性を更20に向上させるため，今回利用した気象庁の観測網だけで21はなく，海底地震観測網や Hi-net，KiK-net のような地中

22に設置された高密度観測網を統合した処理の検討を進め

23

たい．その際，今回経験的に与えた尤度関数等のパラメ

24ータを観測点，観測網ごとに評価することや，観測点の25増加に伴う処理の更なる効率化に向けた検討が今後必要

26になるであろう．

27
28

§5.

まとめ

29

本研究では，先験的な情報と物理量が異なる複数の観30測値を統合的に扱うことのできるベイズ推定を用いて，
22

1

震源をリアルタイムで推定する手法を開発した．また，

2

振幅値やトリガしていない観測点の情報を利用すること3で，ほぼ同時刻に発生した複数の地震を適切に分離する

4

ことが可能となった．開発した手法を 2011 年 3 月 1 日か

5

ら 4 月 30 日に緊急地震速報を発表した 71 事例，及び 20136

年 8 月 8 日に発表した 1 事例に対して適用した．その結7果，従来処理で 22 回あった誤報回数は，本手法では 08回となった．また，不適切な警報を発表する割合も，66%

9（48/72 事例）から 6%（4/72 事例）に大幅に改善された．

10

本手法は，従来の緊急地震速報処理と同じ観測値を入力11として用いているため，既存の枠組みの範囲内でも，緊

12

急地震速報の精度を改善できることが示された．

13
14

謝辞

15

本論文をまとめるにあたり，気象庁の土井恵治氏，束16

田進也氏，中村雅基氏，平野和幸氏，山田安之氏には有17

益な議論をしていただきました．また，東京大学地震研

18究所の長尾大道氏と 1 名の匿名の査読者及び担当編集委

19

員の独立行政法人海洋研究開発機構の有吉慶介氏からの20

ご意見は，本論文の改善に大変役立ちました．本研究で21は，公益財団法人鉄道総合技術研究所と気象庁による共22

同技術開発の成果である現地観測点処理結果を利用しま

23

した．過去の震源分布及び処理結果の比較には，北海道

24

大学，弘前大学，東北大学，東京大学，名古屋大学，京25

都大学，高知大学，九州大学， 鹿児島大学，独立行政法

26

人防災科学技術研究所，独立行政法人産業技術総合研究

27

所，国土地理院，青森県，東京都，静岡県， 神奈川県温28

泉地学研究所，横浜市，独立行政法人海洋研究開発機構

29及び気象庁の地震波形を気象庁と文部科学省が協力して

30

データを処理した一元化震源を使用しました．図の作成
23

1には Wessel and Smith（1998）の Generic Mapping Tool を2使用しました．本研究は，最先端・次世代研究開発支援

3プログラムの助成を受けています．記して感謝いたしま

4

す．

5
6

文献

7

Allen, M. R., P. Gasparini, O. Kamigaichi and M. Böse, 2009,8

The status of Earthquake Early Warning around the9world: An introductory overview, Seismol. Res. Let., 80,

10682-693.

11

Arulampalam, M. S., S. Maskell, N. Gordon and T. Clapp,

12

2002,

13

nonlinear/non-Gaussian Bayesian tracking. IEEE Trans.

14

Signal Proc., 50, 174-188.

15
16
17
18

A tutorial

onparticlefiltersforonline原田智史, 2007, 多機能型地震観測装置の概要, 験震時
報, 70, 73-81.
舟崎淳・地震予知情報課, 2004, 気象庁速度マグニチュー
ドの改訂について, 験震時報, 67, 11-20.

19

Horiuchi, S., H. Negishi, K. Abe, A. Kaminuma and Y.

20

Fujinawa, 2005, An automatic processing system for21

broadcasting earthquake alarms, Bull. Seism. Soc. Am.,22

95, 708-718.

23

Hoshiba, M., O. Kamigaichi, M. Saito, S. Tsukada and N.

24

Hamada,2008,

Earthquake

25

nationwide in Japan, EOS Trans. AGU, 89, 73-74.

26

Hoshiba, M., K. Iwakiri, N. Hayashimoto, T. Shimoyama, K.

27

Hirano, Y. Yamada, Y. Ishigaki and H. Kikuta, 2011,28

Outline of the 2011 off the Pacific coast of Tohoku

29

Earthquake (Mw 9.0) – Earthquake Early Warning and30

observed seismic intensity –, Earth Planets Space, 63,
24

Early Warning starts

1547-551.

2

Hoshiba, M., 2013, Real-time prediction of ground motion by

3

Kirchhoff-Fresnel boundary integral equation method:

4

Extended front detection method for Earthquake Early

5

Warning, J. Geophys. Res. Solid Earth, 118, 1-13.

6

Kamigaichi, O., 2004, JMA Earthquake Early Warning,

7

Journal of Japan Association for Earthquake Engineering,8

4, 134-137.

9気象庁, 2008, 緊急地震速報の概要や処理手法に関する10

技術的参考資料, <http://www.seisvol.kishou.go.jp/eq/

11

EEW/kaisetsu/Whats_EEW/reference.pdf>（参照 2013-128-30）.

13

気象庁, 2010, 緊急地震速報評価・改善検討会 技術部会

14

（第 2 回）, <http://www.seisvol.kishou.go.jp/eq/EEW/

15

Meeting_HYOUKA/t02/index.html>（参照 2013-9-17）.

16

気象庁, 2011, 緊急地震速報評価・改善検討会 技術部会

17

（第 3 回）, <http://www.seisvol.kishou.go.jp/eq/EEW/

18

Meeting_HYOUKA/t03/index.html>（参照 2013-9-17）.

19

気象庁, 2012, 緊急地震速報評価・改善検討会 技術部会

20

（第 4 回）, <http://www.seisvol.kishou.go.jp/eq/EEW/

21

Meeting_HYOUKA/t04/index.html>（参照 2013-9-17）.

22

気象庁, 2013a, 緊急地震速報評価・改善検討会（第 4 回）

23

資料, <http://www.seisvol.kishou.go.jp/eq/EEW/Meetin24

g_HYOUKA/04/index.html> （参照 2013-8-30）.

25

気象庁, 2013b, 緊急地震速報の発表状況, <http://www.se26

isvol.kishou.go.jp/eq/EEW/kaisetsu/joho/joho.html>（参27

照 2013-10-16）.

28

気象庁, 2013c, 緊急地震速報のしくみ, <http://www.seis

29

vol.kishou.go.jp/eq/EEW/portal/shikumi/eew_shikumi.ht

30

ml>（参照 2013-10-17）.
25

1
2

気象庁, 2013d,地震火山月報（カタログ編）平成 25 年 3
月.

3Odaka, T., K. Ashiya, S. Tsukada, S. Sato, K. Ohtake and D.

4

Nozaka, 2003, A new method of quickly estimating

5

epicentral distance and magnitude from a single seismic

6

record, Bull. Seism. Soc. Am., 93, 526-532.

7

Okada, Y., K. Kasahara, S. Hori, K. Obara, S. Sekiguchi, H.

8

Fujiwara and A. Yamamoto, 2004, Recent progress of

9seismic observation networks in Japan –Hi-net, F-net,

10

K-NET and KiK-net–, Earth Planets Space, 56,

11

xv-xxviii.

12

束田進也・小高俊一・芦谷公稔・大竹和生・野坂大輔, 2004,

13

P 波エンベロープ形状を用いた早期地震諸元推定法,14

地震 2, 56, 351-361.

15

Wessel, P. and W. H. F. Smith, 1998, New, improved version

16of Generic Mapping Tools released. EOS Trans. AGU,

1779, 579.

18

Liu, A. and M. Yamada, 2014, Bayesian Approach for

19

Identification of Multiple Events in an Early Warning

20

System, Bull. Seism. Soc. Am., 104-3, (accepted).

21

Yamamoto, S., P. Rydelek, S. Horiuchi, C. Wu and H.

22

Nakamura, 2008, On the estimation of seismic intensity

23

in earthquake early warning systems, Geophys. Res.

24

Lett., 35, L07302.

25
26

横田崇, 1985, 自動検測手法の研究, 気象研究所技術報
告, No. 16, 56-100.

27

横田崇・周勝奎・溝上恵・中村功, 1981, 地震波データの28

自動検測方式とオンライン処理システムにおける

29

稼動実験, 東京大学地震研究所彙報, 56, 449-484.

30
26

1
2

Figure Captions

3
4
5

Fig. 1. Flow charts of (a) the current JMA EEW system and(b) this study.

6
7

Fig. 2. Example of prior probability density at the depth of

8

10km. (a) Past seismicity, (b) Voronoi cell, and (c)9product of (a) and (b). Open triangles show location of

10

stations. A filled triangle shows the first triggered station.

11

Bottom figures show the section at the triggered station.

12Black dots show the seismicity between 1997 and 2012.

13

Note that each plot is normalized. Dashed line in (b)14shows the Voronoi cell of the first triggered station.

15
16

Fig. 3. (a) A function of lower bound g0 in Eq. 9. (b)17Probability density function f with μ=0 and σ=1 in Eq.

18

10. It is normal distribution with the lower bound of g0.

19
20

Fig. 4. Examples of likelihood for (a) B-Δ, (b) travel time,21

and (c) amplitude. Filled and open triangles show

22

triggered and non-triggered station, respectively. A star

23

shows the epicenter location. Note that each plot is24normalized.

25
26

Fig. 5. Example of the Tohoku Earthquake. Distribution of

27

samples at (a) t=2.1s and (b) t=9.1s. Color indicates the28weight of each sample. A cross shows an optimal29

෠ at the time. Open star shows theparameter 𝜃30hypocenter location in the JMA catalog. Time series of27

1(c) estimated magnitude, (d) distance error of estimated2location, (e) estimated depth. Dashed line shows the3value in the catalog. Horizontal axis shows the elapsed

4time from the pick time of the first triggered station.

5
6
7

Fig. 6. Example of two small earthquakes on Mar. 19, 2011.
See Fig. 5 for the caption.

8
9
10

Fig. 7. Example of a small earthquake and noise on Aug. 8,
2013. See Fig. 5 for the caption.

11
12
13

Fig. 8. Performance of EEW issued by current JMA systemand this study.

14
15

Fig. 9. Histogram of the residual between the expected and16

observed maximum seismic intensity scale. (a) Current

17

JMA system and (b) this study.

18
19

Table Captions

20
21

Table 1. List of EEW from Mar. 1, 2011 to Apr. 30, 2011, and22Aug. 8, 2013. Magnitude (M) and maximum observed

23

seismic intensity (Obs. I.) in the JMA catalog, and24

maximum expected seismic intensity (Exp. I.) of current

25

JMA system and this study are listed.

28

Fig. 1(a) JMA EEW

PCA: Principal Component Analysis
LSM: Least Square Method of Amplitudes

B-Δ, PCA,
Voronoi cell

JMA station
Traveltime grid search

LSM

Estimate
Mag. &
Intensity

EEW

Estimate
Mag. &
Intensity

EEW

Not yet arrival determination

Hi-net
1

2

3

4

5~

No. of stations(b) This study
JMA station(RΔ,Ra,Ttrg,Ad,Av)

Past seismicity (Ps)
Prior
Voronoi cell concept (Pv)

B-Δ & PCA (Lb)×Traveltime (Lt)×Amplitude (LM)Particle Filter

Likelihood

Any
other
EQ?
yes

no

Fig. 2

(a)

(b)

(c)
×10-3

Ps

Pv

P=Ps Pv

Fig. 3(a)(b)

Fig. 4

(a)

(b)
×10-3

×10-3

Lb

Lt

(c)

(d)
×10-3

LM

×10-3

L=Lb Lt LM

Fig. 5

(a)

(c)

(b)
(d) Location error
×10-3

OURI
(e)

Fig. 6(a)(c)

TENNOD(b)×10-3(d) Location error(e)Fig. 7(a)(c)

TT3OBS(b)×10-3(d) Location error(e)Fig. 80
JMA EEW

20
24

40

60

26

22

4
This study
Accurate

26
Inaccurate

False

No warningFig. 93020(a) JMA EEW

100-4 -3 -2 -1 0 1 2 3 4 5 6 7 8 93020(b) This study

10
0-4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9Exp. – Obs. Intensity ScaleTable 1
Date
2011/3/11

Catalog
Region
14:46 Off Miyagi Pref.
17:40 Fukushima Pref.
19:35 Off Fukushima Pref.
3:11 Off Fukushima Pref.
3:59 Niigata Pref.
4:08 Off Ibaraki Pref.
4:16 Off Fukushima Pref.
4:31 Niigata Pref.
5:11 Off Miyagi Pref.
5:42 Niigata Pref.
6:19 Nagano Pref.
6:34 Off Fukushima Pref.
6:48 E off Chiba Pref.
22:15 Off Fukushima Pref.
22:24 Off Iwate Pref.
22:26 Off Iwate Pref.
23:34 Niigata Pref.
23:43 Off Iwate Pref.
8:25 Off Miyagi Pref.
10:26 Off Ibaraki Pref.
10:02 Off Ibaraki Pref.
15:12 Off Fukushima Pref.
16:25 Off Ibaraki Pref.
1:36 Tokyo Bay
5:33 E off Chiba Pref.
7:29 Fukushima Pref.
22:31 Yamanashi Pref.
2:40 Chiba Pref.
12:23 Off Fukushima Pref.
12:52 E off Chiba Pref.
21:32 E off Chiba Pref.
6:18 Off Ibaraki Pref.
8:32 Off Iwate Pref.
18:56 Ibaraki Pref.

1234
2011/3/12
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
2011/3/13
20
21
2011/3/14
22
23
24
2011/3/15
25
26
27
28
2011/3/16
29
30
31
2011/3/17
32
2011/3/19
33
34
35
18:57 Ibaraki Pref.
36
37
2011/3/20
14:19 Fukushima Pref.
38
2011/3/22
12:38 E off Chiba Pref.
39
2011/3/23
1:12 Off Ibaraki Pref.
40
7:12 Fukushima Pref.
41
7:36 Fukushima Pref.
42
8:46 E off Chiba Pref.
43
2011/3/25
20:36 Off Miyagi Pref.
44
2011/3/27
19:23 E off Chiba Pref.
45
2011/3/28
7:23 Off Miyagi Pref.
46
2011/4/1
19:49 Akita Pref.
47
2011/4/3
16:38 Off Fukushima Pref.
48
2011/4/4
18:29 Off Fukushima Pref.
49
2011/4/7
23:33 Off Miyagi Pref.
50
2011/4/11
17:16 Fukushima Pref.
51
17:26 Fukushima Pref.
52
18:05 Fukushima Pref.
53
20:42 Off Fukushima Pref.
54
2011/4/12
8:08 E off Chiba Pref.
55
56
10:23 Chiba Pref.
57
12:20 E off Chiba Pref.
58
14:07 Fukushima Pref.
59
16:14 Nagano Pref.
60
2011/4/13
10:07 Fukushima Pref.
61
2011/4/14
6:43 Fukushima Pref.
62
12:08 Fukushima Pref.
63
20:24 Fukushima Pref.
64
21:24 Fukushima Pref.
65
2011/4/15
23:34 Off Iwate Pref.
66
2011/4/16
11:19 Ibaraki Pref.
67
2011/4/19
4:14 Akita Pref.
68
6:33 Ibaraki Pref.
69
2011/4/21
22:37 Chiba Pref.
70
2011/4/24
20:50 Fukushima Pref.
71
2011/4/30
2:04 E off Chiba Pref.
72
2013/8/8
16:56 Wakayama Pref.
*1 Observation Blackout*2 Simultaneous Occurrence*3 The other*4 Two events make three warnings*5 Underestimate*6 Only single station

M
9.0
6.0
5.1
6.0
6.7
5.2
4.0
5.9
6.4
5.3
3.7
4.8
4.6
6.2
5.0
5.4
3.7
5.9
6.2
6.6
6.2
5.2
5.0
3.3
3.6
4.3
6.4
4.0
4.7
6.1
5.7
4.8
5.7
4.0

Obs. I.
7
5+
4
3
6+
4
3
63
63
3
3
53
2
54
54
54
3
2
1
3
6+
2
2
54
2
4
4

6.1

5+

4.7
5.9
5.4
6.0
5.8
5.0
6.3
5.0
6.5
5.0
5.4
4.0
7.2
7.0
5.4
5.1
5.9

3
4
3
5+
5+
2
4
2
55+
4
2
6+
654
5-

6.4

5-

4.2
3.8
6.4
3.1
5.7
4.1
5.4
4.4
3.9
5.0
5.9
4.9
4.8
6.0
3.1
4.7
2.3

2
2
62
53
4
3
3
3
5+
53
53
3
0

Issued warning
Evaluation
Exp. I.
6Accurate
5+
Accurate
5+
Inaccurate*1
5Accurate
6Accurate
5Accurate5+
Inaccurate*25+Inaccurate*2
5+
False*2
5Accurate
6False*2
6+Inaccurate*2
5Inaccurate*3
5Inaccurate*1
5Inaccurate*3
5False*2
6+Inaccurate*2
5Accurate
5Inaccurate*1
5Inaccurate*1
5Accurate
6False*2
6False*2
5False*25+False*2
6+
Inaccurate*2
5Accurate
5+False*2
5False*2
6Accurate
5Inaccurate*3
5False*25Inaccurate*2
5+
Accurate
5False*45+False*4
6Inaccurate*2
6+Inaccurate*2
6+Inaccurate*2
5+
Accurate
5+
Accurate
5False*36Inaccurate*3
5+False*2
6Inaccurate*2
5Accurate5+
Inaccurate*2
5+
False*2
6Accurate
6+
Accurate
5Accurate
5Accurate
6Accurate
7
Inaccurate*2
5Inaccurate*2
5+
False*2
5False*2
6Accurate
6False*36+Inaccurate*2
6False*2
5Accurate
6Inaccurate*2
5Accurate
6+False*2
5Accurate
5Inaccurate*3
5Inaccurate*3
5Accurate
6Inaccurate*2
6+
False*2
7
False*3

This study
Exp. I.
Evaluation
7
Accurate
5Accurate
4(No warning)
5Inaccurate*3
6+
Accurate
5Accurate
4(No warning)
6Accurate
4(No warning)
4(No warning)*5
4(No warning)
4(No warning)
4(No warning)
5+
Accurate(No warning)
4(No warning)
4(No warning)
5Accurate
5Accurate5(No warning)*6
5+
Accurate
4(No warning)3(No warning)3(No warning)3(No warning)
4(No warning)
6Accurate
4(No warning)3(No warning)
6Accurate
5Accurate3(No warning)
4(No warning)
6Inaccurate*3
6-

Accurate

4
4
3
65+
3
53
554
3
67
5+
54

(No warning)(No warning)(No warning)
Accurate(No warning)*6(No warning)
Accurate(No warning)
Accurate
Accurate(No warning)(No warning)
Accurate
Accurate
Accurate
Accurate(No warning)

5-

Accurate

2
6+
4
53
4
54
3
6+
63
54
4
-

(No warning)(No warning)
Accurate(No warning)
Accurate(No warning)(No warning)
Inaccurate*3(No warning)(No warning)
Accurate
Inaccurate*3(No warning)
Accurate(No warning)(No warning)(No warning)