揺れ検知から震央を検出してみる
はじめまして。
この記事は 防災アプリ開発 Advent Calendar 2023 13日目です。
12日目は、べにだて さんの 【気象庁HP】推計震度分布図のGeoJSONデータを無料で取得したい！！ でした。
はじめに
地震が起きたとき、地震の発生場所をすぐに知りたくなりませんか？
私はなります(宣言)。
多くの場合、地震の発生から数秒で気象庁から緊急地震速報（予報）が発表されるので、その情報から震央の位置を知ることができます。
しかし、緊急地震速報（予報）には発表条件があります。この条件を満たさない地震では緊急地震速報（予報）が発表されず、震央の位置はわかりません。
いずれかの地震観測点において、P波またはS波の振幅が１００ガル以上となった場合。
地震計で観測された地震波を解析した結果、震源・マグニチュード・各地の予測震度、予測長周期地震動階級が求まり、そのマグニチュードが３．５以上、または最大予測震度が３以上、長周期地震動階級が１以上である場合。
気象庁の発表する緊急地震速報以外でリアルタイムに地震が発生したことを知るには、防災科研の強震モニタがあります。
とくに強震モニタから揺れを検知する機能があるアプリを使用していると、地震が発生して揺れを検知していても、緊急地震速報が発表されていないから震央の位置が分からない！ということが多いのではないでしょうか。
(旧バージョンの画像しかなくてごめんなさい)
ただ、揺れの検知をよく見てみると、地震による揺れの広がりは震央を中心にだいたい同心円状になっています。それなら、揺れの広がりから震央を見つけることができるのでは？と思いつきました
ということで、今回は揺れ検知から震央の位置を検出してみたというお話です。
震央の検出方法、どうしよう
揺れ検知時間と距離の関係は?
地震の揺れはだいたい同心円状に伝播していきます。これは学校でも習った通り、地震波は震源を中心に P波が約7km/s, S波が約5km/s で地中を伝わるからですね。
ということは、震央からの距離が同じ観測点は揺れを検知する時間もほぼ同じになり、距離が遠いほど検知時間は遅くなります。
そこで、揺れ検知の時間と震央からの距離 の関係をグラフにしてみました。
（2023/06/11 苫小牧沖の地震）
これでわかることは、震央の位置が実際の震源に近いほど、時間と距離の関係がある曲線に近づいているということです。
実際の震央に近いと現れるこの曲線は、検知時間と震央距離の関係を表していそうですね。この関係がわかれば、揺れ検知時間をこの曲線に対応させて震央距離が求められそうです。
ということは逆に、
調べたい震央位置を決めて、各観測点の震央距離を計算
この曲線を使って震央距離から到達時間を計算
到達時間が実際の揺れ検知時間にどのくらい近いかを調べる
という感じの流れで、調べたい震央位置が本当の震央に近いかどうかを計算できそうです！
地震波の走時曲線
問題は、この曲線をどうやって求めるかですよね。
気象庁が震源を計算するためにどのような方法を使っているのかを調べてみると、主に JMA2001走時表 から地震波の理論走時を求めて、震源の位置を計算しているそうです。
気象庁では震源から観測点までの理論走時の計算に走時表を使用して震源計算を行っている。
…（中略）…
陸域の観測点ではJMA2001A、日本海溝の陸側の観測点ではJMA2020A、日本海溝の外側の観測点ではJMA2020B、そして南海トラフ領域の観測点ではJMA2020Cという4種類の走時表を震源計算に使用することとした。
※JMA2001A走時表のtjma2001h.+00000はJMA2001走時表のtjma2001と同じ
この走時表にはP波の走時・S波の走時・深さ・震央距離の関係が記述されています。
つまり、検知時刻と震央距離の関係を表していそうなあの曲線の正体(?)と言っても良いと思います。
JMA2001走時表に記述されている距離や深さは飛び飛びの値ですが、間の値を補間すれば震央距離と深さから地震波の走時を求められます。
試しに、先ほどグラフを計算した2023/06/11 苫小牧沖の地震について、震度データベースで公開されている発生時間と深さを使って到達時刻を計算してみました。※揺れ検知はP波によるものであると仮定して、走時はP波のものを使用しました。
走時と比べて検知時間が遅れ気味な気もしますが、だいたい一致していそうですね。
実際の震央にどのくらい近いか計算する方法
これまでの話をまとめると、対象となる震源（緯度・経度・深さ・発生時刻）を決めたとき、震央が実際の震央の位置にどのくらい近いかを調べるには…
緯度・経度から震央距離を計算
震央距離・深さからJMA2001走時表を使って走時を計算
走時・発生時刻から到達時刻を計算
計算した到達時刻が実際の揺れ検知時刻にどのくらい近いかを調べる
という感じになります。
この方法はわりと使えそうなので、いよいよ実際に震央を求めるための作戦を考えましょう。
震央を求めるための手順
ここまできてもうお気づきだと思いますが、震央を検出したいだけだったのに、発生時間や深さも必要になってきて、もはや震源を計算するのと同じですよね。
気象庁の震源計算方法を見てみる
先ほど参照していたJMA2001走時表について書かれていたページに地震月報(カタログ編)の震源計算方法が書いてあったので、ちょっと見てみます。
震源計算手法は、Geiger(1910)の方法(1)を改良した反復法［浜田ほか（1983） (2)]を用いる。
ほう。ガイガーの方法と言われても分からないので調べてみる。
まず大体の震央を定めて, 各観測所からその点まで震央距離Δsを計算する。1つの位相(例えば P)の各所の発現時の観測値を, 縦軸が時間, 横軸が震央距離に取ったグラフの上に記入する。そしてその位相の走時曲線は, それからの各点の時間的誤差の和が零になるような位置にあるのが正しいものと仮定する。
…（中略 すごそうな数式など）…
この同次方程式を解いて x, y, z を求める。かくして求めた値が大きかったら, これを第一近似とし, 上の方法を繰返して第2, 3近似値を逐次求める。
？？？ 内容が高度すぎてわからない…
大体の震央を定めたら走時曲線の傾斜と最小自乗法で各観測点の誤差が0になる場所を計算し新しい震央と決めてそれを何回か繰り返す(？)
みたいなことが書いてありますがよくわかりません。まあ、ちょっと出てきた震央距離と時間のグラフに関しては似たようなものを作ったけれども…
ということで、この流れをちょっと参考にしつつ、(結局)自己流で手順を考えました。
① 揺れ検出時刻を記録する
最初にグラフを作成したときもやっていたものです。
観測点を揺れ検知済み状態に設定するときに、その時点での現在時刻を記録するようにしました。
② 仮の震源を決める
地震波の到達時刻を求めるには、震源（緯度・経度・深さ・発生時刻）が必要でしたね。
ガイガーの方法では、まず大体の震央を定めると書いてあります。
今回は、最初に揺れを検知した観測点の位置を仮の震源として定めました。（緯度・経度は観測点のものを小数第2位で四捨五入した値、深さは10km、発生時刻は検知時刻-2秒）
（震源位置が空の場合 or 最初の検知から10秒以内 の場合に実行）
③ 仮の震源が実際の震源にどのくらい近いかの値を計算する
つぎに、この震源が実際の震源にどのくらい近いか決めるため、誤差レベルという指標を作りました。この値が小さいほど誤差が小さい＝実際の震源に近いことになります。
誤差レベルは以下の方法で求めるようにします。
すべての揺れ検知済み観測点で震央距離と走時を計算する
各観測点で、揺れ検知時刻から走時を引いて発生時刻を計算し、記録しておく
各観測点で、震央距離から重みを計算して記録する。重みは震央距離が近いほど大きくなるように 最初に揺れを検知した観測点の震央距離÷震央距離 とする。ただし震央距離50km以内であれば1に固定する
すべての観測点の発生時刻の平均値を求める
各観測点で、2.で記録した発生時刻から発生時刻の平均値を引いた値を2乗して2乗誤差を求める
2乗誤差に重みをかけたものを合計して、その値を誤差レベルとする
各観測点で地震発生時刻を計算して、そのばらつき(分散)を使って誤差レベルを計算する感じです。
紫の点: 揺れ検知時刻, オレンジの線: 走時(到達時間), オレンジの点: 計算した発生時刻
こんな感じで、実際の震源に近いほど各観測点で計算した発生時刻のばらつきが少なくなり、実際の震源から遠いと発生時刻のばらつきが大きくなります。
そのばらつきに重みをつけてそのまま誤差レベルにしちゃいました。震源パラメータの1つ発生時刻も計算できるので一石二鳥ですね。
でも、このままでは問題がありました。
というのも、揺れ検知済みの観測点が少ないときに、このように明らかに実際の震源から遠い震源で誤差レベルが小さくなってしまいます。
（2023/09/27 大阪府北部の地震）
確かにグラフを見ると、下の方が走時曲線にぴったりフィットしています。でも震央は揺れを検知していない観測点のど真ん中になっていて、明らかに間違いですよね。
この現象をなんとかするために、揺れを検知していない観測点も使用する着未着法を参考にしてみます。
この場合, より遠方に位置する観測点 𝐢 では, 計器の故障でない限り, 後で必ずP波を観測するはずである。このことは, この地震が「ある時刻 𝐓ⁿᵒʷ までに観測点 𝐢 で観測されてはならない」という条件で発生したことを意味している。
まだ揺れを検知していない観測点で走時を計算して、もし到達済みになるようなら誤差レベルを増やすという感じですね。
ということで、揺れ検知から3秒以内、または10秒以内で揺れ検知数30点未満の場合は、以下の方法で誤差レベルをさらに増やします。
揺れを検知した周囲のグリッド内で、最大震央距離+30km以内にある揺れ未検知の観測点を選び、震央距離と走時を計算する
揺れ検知済み観測点の発生時刻の平均値と走時を足して、到達時刻を計算する
到達時刻より現在時刻が大きくなる(計算上は到達済みになるような)観測点では、誤差レベルに1を足す
これにより、揺れを検知していない観測点が到達済み範囲に入るような震源は、誤差レベルが大きくなるようになります。
先ほどの大阪府北部の地震の例でも、より実際の震源に近いほうが誤差レベルが小さくなりました。問題解決です。
上：実際の震源に近い位置 下：誤差レベルが最も小さくなっていた位置
ここまででようやく、実際の震源にどのくらい近いかの値＝誤差レベルを計算できるようになりました！
一応誤差レベルを計算するコードを載せておきます。
謎の変数も登場してますが正体不明です
（文字とか小さすぎてつぶれちゃったかも）
④ 誤差レベルの計算を繰り返して震源位置を特定する
今決めたのはまだ仮の震源なので、より実際の震源に近そうな新たな震源を決める必要があります。
ガイガーの方法では何かの方程式を解いて新たな震央位置を計算して、ある程度小さくなるまで繰り返し計算しているようでした。
…が、数学がまったくわからない私には一体何が行われているか想像もつきません。
ここはもうゴリ押しとして、新たな震源候補すべてに対して誤差レベルを計算し、最も誤差レベルが小さい震源を新たな震源に決定、ある程度小さくなるまで繰り返すようにしたいと思います！！！
ただ、ゴリ押しといっても震源候補を決める必要があります。
誤差レベルは実際の震源に近いほど小さくなり、遠くなるにつれて大きくなります。ということは、誤差レベルは震源を底とした谷になっているはずです。
そこで、仮の震源の周りに震源候補を用意して、その中で最も実際の震源に近そうな候補に仮の震源を移動させます。こうすれば谷底に転がり込むように実際の震源に向かって移動させられそうです。
ということで、具体的には以下のように仮の震源を移動させます。
仮の震源から深さはそのままで、緯度・経度を±0.5度変化させたの4つの震源候補で誤差レベルを求める
仮の震源より誤差レベルが低い震源候補があれば、その位置に仮の震源を移動させる
手順1～2を繰り返して、仮の震源より誤差レベルが低い震源候補がなくなれば、仮の震源の位置を新しい震源と決める
手順1の震源候補を 緯度・経度を±0.1度変化させた4つにして、もう一度手順1～3を行い、新しい震源を決める
手順1の震源候補を、深さを±50km, 緯度・経度を±0.1度変化させた6つに増やして、もう一度手順1～3を行い新しい震源を決める
手順1の震源候補を、深さを±10km, 緯度・経度を±0.1度変化させた6つにして、もう一度手順1～3を行い新しい震源を決める
この移動が終わったとき、だいたい仮の震源の位置は誤差レベルの谷底＝実際の震源に最も近い場所にあるはずです。
ガイガーの方法を改良した反復法にあるとおり、まずは深さを決めずに震源を求め、その後深さも変化させて最終的に震源の位置を決定するようにしました。
また、最初から変化量が0.1度だと、はるか沖や深発などで最初に決める仮震源との実際の震源が遠い場合に、計算回数が非常に多くなってしまいます。
そのため、最初は0.5度,50kmの変化で大まかに震源の位置を求めて、その後0.1度,10kmの変化で細かい震源の位置を決定するようにしました。
こちらも仮震源の移動と震源検出のコードを載せておきます。
（誤差レベル計算の定義を使用）
⑤ 震源位置を決定する
仮の震源を移動させて、誤差レベルの谷で止まった地点を最終的な震源として決定します。
これで、震央の検出は完了です！！めでたしめでたし。結局のところ震源を求めていたわけですが…
最初に観測点の位置を仮震源にする → 誤差レベルを計算する → 仮の震源を移動させる
という流れを実際に可視化してみました。(GIFでは長すぎたので動画です)
仮の震源がだんだん正しい震源に向かってゆく様子がわかりますね。
実際に震央検出を動作させてみる
せっかく震央検出機能ができたので、いろいろ試してみましょう！
とりあえずさんざん実験台にした2023/06/11 苫小牧沖の地震で試しました。震度データベースの内容と比べてみても、これは結構精度がよさそうですね。
最近のEEWが出なかった地震で試す
今回の目的は、EEWが出ない地震で震央を知りたいというものです。最近EEWが出なかった地震で試してみます。
深さはガバガバですが、震央は結構すぐに決まりました。岩手沖でも4秒くらいで震央位置が安定してくれましたね。
ちょっと残念な結果になる地震
続いてあまり震央の精度が良くなさそうな地震です。
検知する観測点が少ないので難しそうな能登半島沖の地震です。これは途中まで良かったのですが、いきなり明後日の方向へ飛んでいき、岩手県沿岸北部を震央とするようになっちゃいました。着未着法リスペクト機能は検知が少ない最初のうちしか動作しないので、この場合は動かないんですよね…
かなり沖な十勝沖の地震です。これはどちらかというと揺れ検知が足を引っ張ってP波の検知が遅れたため、遠くに震央を決めすぎた感じがします。
珍しい地震で試す
最後は単純に試したかっただけの地震です(満足)。
東海道南方沖の地震です。伊豆諸島に観測点があるので、すぐに震央が決まってくれました。ちなみに揺れ検知側の機能なので触れていませんでしたが、S波と思われる検知にはS波の走時を使っています。
日本海中部の深発地震です。計算する観測点の数がかなり多くなってしまうタイプの地震ですが、仮震源を移動させる定義にあった「1フレーム休み判断」というブロック定義のおかげでカクカクにならずに済んでいる例です。
震央検出の精度について
色々な地震で実際に動作させてみた結果、良い感じに震央位置を検出できました！単純な作りにしては結構正しく検出できてそうです。
ただ、実際どのくらい正確に震央を検出できているのか気になりまよね。というか私が気になります。
せっかくなので、震央の検出結果と本当の震源を比較してみたいと思います！
まず、緊急地震速報が発表されないような地震の震央を知りたいという当初の目的に合わせるため、以下のような条件に当てはまる地震を選びました。
現在から半年以内に地震情報を発表した地震
緊急地震速報が発表されないM3.0~M4.0くらいの地震
揺れを検知できなかった地震は除外
同じ震央地名の地震は除外
（手動計測だから数が多いと大変なため）
次に、実際にこれらの地震をリプレイして検知から20秒後の震央検知位置を記録し、地震情報の震央から緯度経度の差を計算しました。
地図に表すとこんな感じです。
(赤い丸が検知した震央、青い×が地震情報の震央)
の実際の震央との距離の分布はこんな感じです。今回は使用する地震に条件を付けているのですべての地震に対する精度はわからないですが、だいたい65%くらいがほぼピッタリになったようです。
0.3度以内なら90%くらいになるので、数字で見ると意外とよく検出できてるように見えますね！
ただ地図をよく見ると、沖で発生する地震が陸寄りになる傾向がありますね。実際は条件を変えたらもっと精度が悪くなるかもしれません。
実際、こんなカオスなことになったパターンもありますし。
ああああ pic.twitter.com/NRvWq1XKa3
— ことほの (@kotoho76) June 17, 2023
おまけ：緊急地震速報の震源推定方法
ここまで震央(震源)を求める方法をひたすら考えて作成してきましたが、ふと思いました。
緊急地震速報で使用されている方法を使えばよいのでは…？
まあ内容が高度すぎて無理だったんですけどね()
ともかく、気象庁の緊急地震速報ではどのような震源推定手法が使われているのでしょうか。
これは今年(2023年9月)に新しくなったので知っている人も多いと思います。
緊急地震速報の震源推定手法について、従来のIPF法（※）を含む複数の手法の併用から、改良を加えたIPF法に一本化する運用を開始します（詳細は別紙参照）。
気象庁の緊急地震速報の震源決定には、改良を加えたIPF法が使われています。
以前はその他の震源推定手法も使用して震源を決め、同一地震判定システムが、発表に用いる震源とマグニチュードを決めていたようです。
しかし、そのシステムが判定を誤り、過大予測により緊急地震速報（警報）を誤って発表する例があったため、2023年9月26日からは改良を加えたIPF法に一本化されたようです。
改良を加えたIPF法、ざっくりこんな雰囲気
私も詳しく解っていないので間違ったことを書いているかもしれませんが、雰囲気だけでも解説してみます。
有識者の方…誰か詳しい解説記事書いてください…！(他力本願)
詳しい資料:
図が見やすい資料:
① トリガ処理
いろいろな観測網から集めた波形データを1秒ごとに受信してトリガ処理を行います。波形にフィルターをかけて、Tᵖᵈ法というものを使ってトリガ判定を行っているようです。
ちなみにここが改良を加えたIPF法の改良ポイントで、従来のIPF法はトリガ機能が最初から搭載されている観測網のデータしか使えなかったみたいです。
② 地震検知
観測点がトリガすると、その観測点や近くの観測点も含めたトリガグループを作成します。
そのグループのうち2~3点以上がトリガしていたり大きな揺れであれば継続地震EQとマークされ、すぐに震源推定が始まります。
そうでなければ保留中の地震EQpとしてマークされ、この状態でしばらくするとトリガが削除されるようです。
トリガグループの作成範囲 (a)近い順に5個 (b)近い順20個+ボロノイ分割の中心から近い順20個
③ 震源推定処理
ここからは従来のIPF法とほぼ同じっぽいです。
まず、水平0.5度、深さに50kmくらいのグリッドを作ります。
これまでの地震の発生状況から求めた確率と、最初にトリガした観測点から最も近いグリッド周辺で高くなる確率を合成して仮定震源を配置するための3D確率分布を決定します。
事前確率の設定方法 (a)が震源の分布による確率 (b)が最初にトリガした最短範囲による確率
次に、この確率分布に従って、大量の仮想震源を配置します。
配置したすべての仮想震源でどのくらい実際の震源に近いかをを示す尤度を計算します。尤度の計算には検測時刻以外にも以下のデータを使っているようです。
検測時刻
地震波の最大振幅
B-Δ法で計算した震央距離
主成分分析法で計算した震源方位
これらの実測値との差を尤度関数というものに入れて、出てきた数字を使って尤度を計算するようです。また、この時最も尤度が高かった震源を実際の震源候補とします。
震源候補が決まって一安心と思いきや、まだ続きます。
尤度から重みを計算し、ベイズの定理というものを使ってもう一度3D確率分布を設定します。
再び設定し直した確率分布に沿って、もう一度大量の仮想震源を設置します。この時、収束を防ぐために緯度経度±0.1度、深さ標準偏差20kmでランダムに震源を動かすようです。
この仮想震源からもう一度尤度を計算し、尤度から確率を計算し、再び仮想震源を設置し… と繰り返してゆくうちに、仮想震源群が実際の震源にどんどん集まってゆきます。
パーティクルフィルタ法を使った流れ
つまり、計算が進むにつれて確率分布がどんどん実際の震源に近くなっていき、計算するべき範囲も絞ることができるみたいです。
④ 複数地震の識別
IPF法は複数の地震が発生した場合の識別が得意です。
新たなトリガ済み観測点が増えた際、そのトリガが同一地震の揺れであるかを判断する同一地震判定を行います。
同一地震判定には、振幅値(M残差)と検測時刻(走時残差)を使った条件式を使います。
この条件式に当てはまる場合は同一地震とみなし、尤度計算に使用されます。
もし条件式に当てはまらない場合は、新たな地震が発生したとみなし、その観測点を最初にトリガした観測点として、③の震源推定処理を別の地震として新たに行います。
(a)で処理中の地震に (b)で新たな地震の観測データが入り、2つの地震を判別して尤度を計算した例
⑤ 震源の収束判定
グループ内で最大理論P波到達時刻から5秒間震源の推定が安定している場合、継続地震EQは収束地震EQcとしてマークされます。この状態では再び新しい地震のトリガを行うことができるようになります。トリガが発生した場合、再び②の地震検知を行います。
ちなみに、M5級の地震では約30秒で収束地震EQcになるようです。
⑥ 処理の終了
システム全体でP波が検出されなくなり、しきい値である300~900秒が経過すると、収束地震EQcは削除され処理が終了します。
と、ここまで書いている自分自身もよくわかっていないですが、ものすごく高度なことをやっていることはなんとなくわかります。
今回作成した自己流の震源検出には揺れ検知時刻しか使っていなかったのに対して、改良を加えたIPF法は検測時刻、最大振幅、震央距離、震源方位の4つのデータをくまなく使っていました。しかもマグニチュードの計算や複数地震の判別まで同時に行ってしまいます。
普段何気なく使っている緊急地震速報が、いかに高度な技術を使っているのかがよくわかりますね。
おわりに
最後まで読んでいただけたでしょうか。ものすごく長い記事になってしまい申し訳ありません。というか、私自身もここまで長くなると思っていませんでした…
防災アプリ開発 Advent Calendar は、自作の地震関連アプリを作成する際にも様々な記事を参考させてもらいました。
今回はあまり防災に関係のない記事でしたが、このアドベントカレンダーをきっかけにもっと防災に関心を持つ人が増えてほしいですね。
防災アプリ開発 Advent Calendar 2023 14日目の担当は、この記事でも登場したJQuakeの作者、フランソワさんです。楽しみですね！
ちなみに、最終的には震源を求めていたのに、なぜ震央にこだわっていたかという話ですが…
実は、震源検出計画を始めた最初の理由は、自作の地震ソフトで揺れ検出範囲を円で囲いたかったからです。
円で囲うにはもちろん震央位置が必要なので、震央を検出できないか試行錯誤していくうちに、いつの間にか震源を求めることになってました。
なお、実際に震央位置を求めるために、この記事に書いてある手法を少し改良した方法を使っていました。
ですが、せっかくの震央検知機能もScratchのコード5MB制限に引っかかってしったうえ、強震モニタの利用条件改定などもあり、今はコードごとごっそり消えてなくなりました(無念)。この記事はその供養も兼ねています。
https://turbowarp.org/636244032/editor
この記事が気に入ったらサポートをしてみませんか？