これは地震情報アプリ界隈 Advent Calendar 2021の16日目の記事です。

概要

距離減衰式は震度計算は細かい震度予測ができない反面、計算が非常に簡単なもののため、緊急地震速報など即時性が必要とされる場合に使われます。

これは、自分への記録です。

後半にはTypeScriptによる実装を例を掲載します。

計算式

M変換

距離減衰式では$M_W$を使用しますが、緊急地震速報で算出される$M$は$M_{JMA}$（気象庁マグニチュード）のため、宇津(1982)の経験式により$M_{JMA}$を$M_W$へ変換します。

$M_W = M_{JMA} - 0.171$

最短距離の計算

宇津（1977）の式を用いて$M$から相似則により予測される断層の長さを算出し、その1/2を半径とした球を震源断層と仮定します。その中心を震源に設置し、その球面からの距離を最短距離とします。
下記の式では、断層長$L$を算出します。

$log_{10}L=0.5M-1.85$

ここでの$M$は$M_W$として計算します。

下記の式では断層長$L$の半径と震源距離$x$から最短距離$X$を算出します。単位はいずれもkmです。
計算した地点と震源の距離が3km未満となるについては、一律で最短距離を3kmに設定します。

$X=\max(x - (L・0.5), 3)$

距離減衰式

推定された$M_W$を用いて距離減衰式（司・翠川、1999）により平均S波速度速度600m/sの工学基板上の最大速度、$PGV_{600}$を求めます。

$log_{10}(PGV_{600})=0.58M_W+0.0038D+d-log_{10}(X+0.0028・10^{0.5M_W})-0.002X-1.29$

ここで、 $D$ は震源の深さ、 $d$ は地震タイプ別の係数（内陸地殻地震: 0、プレート間地震:-0.02、プレート内地震:0.12、緊急地震速報:0）、 $X$ は最短距離です。

最大速度の補正

工学的基板上の推定最大速度を補正するために、地盤増幅度を乗算します。
NIEDが表層地盤情報提供APIを公開していますので、これに出現する「工学的基盤（Vs=400m/s）から地表に至る最大速度の増幅率」を使います。
気象庁予報業務許可に適合するためにはV2を使用してください。

推定最大速度と補正値となる増幅率の工学的基盤のS波平均速度が異なるため合わせます。工学的基盤（Vs=600m/s）から工学的基盤（Vs=400m/s）へ変換を行う増幅率(松岡・翠川[1994])を求め、それにより変換を行います。

増幅率とは、工学的基盤上からどれぐらい増幅されて地表面・別の工学的基盤上で揺れるか示す倍率です。

工学的基盤（Vs=600m/s）から工学的基盤（Vs=400m/s）へ変換を行う増幅率は次で求めます。

$ARV=(600 ÷ 400)^{0.66}\cong1.31$

ここからは簡便的な増幅率、1.31を使用します。

$PGV_{400}=PGV_{600}×1.31$

として、

$PGV_S=PGV_{400}×ARV$

で、地表面での推定最大速度$PGV_S$を得ます。

ただし、$ARV$は先ほどAPIで説明した「工学的基盤（Vs=400m/s）から地表に至る最大速度の増幅率」です。

計測震度の算出

推定最大速度が出れば震度は楽勝です。

$I_{INSTR}=2.68+1.72・log_{10}(PGV_S)±0.21$

$I_{INSTR}$は計測震度です。

実装

実装は**"私が"**いつも使うTypeScriptちゃんです。

震度計算の震源初期値は仮で次の通りとします。

// マグニチュード(Mjma)
const magJMA = 7.0;

// 深さ(km)
const depth = 10;

// 震央(ここでは緯度39.5°,経度135°の日本海)
const epicenterLocaltion = [39.5, 135];

Mの変換及び震源断層長を求めます。

// マグニチュード(Mw)
const magW = magJMA - 0.171;

// 断層長(半径)
const long = 10 ** (0.5 * magW - 1.85) / 2;

次は、任意の地点で予測される工学基板上(Vs=600m/s)の最大速度を求めます。

// 予測したい場所(ここでは福島県いわき市役所)
const pointLocaltion = [37.050475, 140.887327];

// 震央と予測したい場所の直線距離を求める(km)
// 関数の内容は省略する
const epicenterDistance = geoDistance(epicenterLocaltion, pointLocaltion);

// 断層長を引いた震源距離を求める
const hypocenterDistance = (depth ** 2 + epicenterDistance ** 2) ** 0.5 - long;

// 計算で使用する震源距離の最短は3kmなので、大きいほうを取る
const x = Math.max(hypocenterDistance, 3);

// 工学基板上(Vs=600m/s)の最大速度
const gpv600 = 10 ** (
      0.58 * magW +
      0.0038 * depth - 1.29 -
      Math.log10(x + 0.0028 * (10 ** (0.5 * magW))) -
      0.002 * x
);

工学基板上(Vs=600m/s)から予測地点ごとに補正をします。

// 予測する地点の工学的基盤（Vs=400m/s）から地表に至る最大速度の増幅率(ここでは簡便的に1.0とする、つまり増幅しない)
const arv = 1.0;

// 最大速度を工学的基盤（Vs=600m/s）から工学的基盤（Vs=400m/s）へ変換を行う
const pgv400 = gpv600 * 1.31;

// 地表面での推定最大速度を求めます
const pgv = pgv400 * arv;

最後に計測震度を求めます

// 予測する地点の地表面推定最大速度から計測震度への変換
const intensity = 2.68 + 1.72 * Math.log10(pgv);

以上、完了です。

法律の話

気象業務法により予報業務許可を取得していない人が、リアルタイムに推定震度など計算した結果を第三者（家族や会社外）に提供すること（アプリケーション・Web・ライブ配信など）を禁止しています。
ですので、身内での利用にとどめるか素直に予報業務許可を取るようにしてください。

・・・

~~個人で地震動予報業務を取得している人はいませんが。~~
→2023年10月に予報業務許可（地震）第237号を取得しました。

まとめ

数学の知識が少しあればすぐ計算できますので、皆さんぜひ試してみてください。

急いで書きましたので、誤字や誤った表記があるかもしれません。遠慮なく指摘してください。

以上、ありがとうございました。