改正メルカリ震度階級から最大速度を求める式を作り、気象庁震度階級を計算する(式の例付き)

Last updated at 2024-01-09Posted at 2023-03-04

注意

語彙力がないので所々おかしいですが頑張って解読してください。
私は詳しいわけではないのでご了承ください。間違っていたら教えてください。(対数関数を使ってますが習ってないときに書いたのでおかしいかもしれません)
自分の解釈や想像が含まれているので詳しくは各自調べてください。
データは適当に抽出しているので精度は低めです。

前書き

※長いので軽く読み飛ばしてもらってもOKです　本編(導出)
WorldQuakeViewerを作ってる人です。USGSなどでは揺れの大きさを表す指標として改正メルカリ震度階級(MMI)が用いられています。海外で大きな地震が起きた時「気象庁震度階級○～○相当」と聞いたことがあると思いますが違和感を感じたことはありませんか？(特にウェザーニュースの記事等) マグニチュードの割に弱かったり。(MMIが小さく出ることもよくありますが。)ここではUSGSのデータを参考にしてMMIから気象庁震度への変換をやってみたいと思います。

そもそも改正メルカリ震度階級とは

MMIは気象庁震度階級とは異なり、計測データだけではなく過去の地震の被害状況をもとに被害の大きさがに分かるようになっています。 つまり、土地によって変わります。(USGSの解説)

USGSの改正メルカリ震度階級

とはいえ、発生後すぐにわかるわけではないためUSGSでは最大加速度と最大速度で求めています。ほぼM5.5(アメリカの内陸だともっと小さくM3台でも出ることも)くらい以上あると発生から20分~30分くらいでMMIが発表されます。Feedには細かく3.456などまで出ます。

「気象庁震度階級○～○相当」というのはUSGSのデータを対応表に照らし合わせて求めていると思います。そこで注意したいのは、ShakeMapの数値はあくまで推定ということです。またDid You Feel It?は各ユーザーの体感です。気象庁のように観測データの最大というわけではありません。(とはいえ日本も観測された以上に局地的に震度○ということはよくありますが)StationListのものはDownloadsのjsonを見る限り同じPGVでもMMIが違っているのでちゃんとしたやつ・・・だと思います。

住宅地から離れたところで起きた地震、例えば2023/02/23のタジキスタンの地震を見てみると、ShakeMapではVIIに対しStation ListではIVとなっています。PAGERの通り、内陸M6.8という規模でも被害は少ないです。

他にも耐震性や津波の被害などを考えると海外の被害を推測するにはPAGERが適切です。USGSのページによると、M5.5以上の地震が起こると20分以内に人口や過去の被害を考慮して被害を推定するようです。Feed見る限りもうちょっと遅い気がしますが…(基本MMIが出た後1分以内にPAGERが発表されます。MMIが大きいとPAGERはもう少し遅いです。(2023/02/06のトルコの地震1回目は第一報12分、MMI20分、PAGER24分))

改正メルカリ震度階級と最大速度

USGSのShakeMapには↓の画像がありますが、どうやらWorden et al.(2012)だそうです。

~~調べてもよくわかりませんが、~~ 防災科研の解説(必ずしも公式見解ではない)によると最大加速度より最大速度のほうが良い相関があるといわれているらしいのでPGVのほうを使います・・・と言いたいところですがこの表だけでは足りないので、ShakeMapのStation ListのDownloadsのStation Listのjsonデータにある観測点のPGVとMMIを使って式を作っちゃいましょう。

気象庁震度階級と最大速度

J-RISQ地震速報にも使われている、藤本･翠川(2005)の換算式を用います。

I=\left\{\begin{align}
2.165+2.262\log_{10}(PGV)\,\,\,\,\,
(I<4)\approx(PGV<7cm/s)\\
2.002+2.603\log_{10}(PGV)-0.213\log_{10}(PGV)^{2}\,\,\,\,\,
(I\ge4)\approx(PGV\ge7cm/s)\\
\end{align}\right\}

本編(導出)

ShakeMapのStation ListのDownloadsのStation ListのjsonデータのPGVとMMIをcsvにでも抜き出しましょう。場所によって差があるのでいろいろな地震をまとめるなどしましょう。(地域ごとに出してみるのもいいかもしれない)
トルコの地震などPGVが大きいものは少ないのでしっかり入れましょう。(PGV215.34cm/s,MMI9.8観測したところもあります)

C#の例

//using Newtonsoft.Json.Linq;
string JsonText = "";//jsonデータ
JObject json = JObject.Parse(JsonText);
string csv = "PGV,MMI";
foreach (JToken json_ in json.SelectToken("features"))
{
  if ((string)json_.SelectToken("properties.pgv") == "null")
    continue;//null(体感のもの(UTMが付くやつ)?)があると邪魔だしExcelの散布図が作りにくいから飛ばす　たぶんmmiはnullにはならない
  csv += $"\n{json_.SelectToken("properties.pgv")},{json_.SelectToken("properties.intensity")}";
  Console.WriteLine($"{json_.SelectToken("properties.pgv")},{json_.SelectToken("properties.intensity")}");
}
File.WriteAllText($"output.csv", csv);

これでできたcsvをExcelなどに読み込ませて散布図を作りましょう。

トルコの地震の例。対数近似はずれるので手動で調整しましょう！
別のところに適当に式を入れて一定区間ごとに計算し表示してみましょう。
(もっといいやり方があるかもしれない。VBAでFunctionに式を定義し、Subを実行し適当なセルの値を変えると楽)

VBAの例

Function MMI(PGV As Double) As Double
  If PGV < 4 Then
    MMI = 3.5 + 0.65 * WorksheetFunction.Ln(PGV)'LnかLog10かはご自由に
  Else
    MMI = 2.5 + 1.37 * WorksheetFunction.Ln(PGV)
End Function

Sub change()
  Cells(1, 28).Value = "" '(任意のセル)　これ実行すると再実行されるので楽
End Sub

こんな感じ。
ここで、1つの式に表すとずれるので藤本･翠川(2005)の換算式のように2つに分けるといいでしょう。
1.一定以下のPGVを無視して近似対数関数を作る
2.無視したところの近似対数関数を作る
という風にやりましょう。

PGVが小さいところは難しいですが頑張りましょう。

式の例

MMI=\left\{\begin{align}
3.5+0.65\ln(PGV)\,\,\,\,\,
(PGV<4cm/s)\approx(MMI<4.4)\\
2.5+1.37\ln(PGV)\,\,\,\,\,
(PGV\ge4cm/s)\approx(MMI\ge4.4)\\
\end{align}\right\}

PGV=\left\{\begin{align}
e^{((MMI-3.5)/0.65)}\,\,\,\,\,
(MMI<4.4)\approx(PGV<4cm/s)\\
e^{((MMI-2.5)/1.37)}\,\,\,\,\,
(MMI\ge4.4)\approx(PGV\ge4cm/s)\\
\end{align}\right\}

常用対数にしてみると、

MMI=\left\{\begin{align}
3.5+1.5\log_{10}(PGV)\,\,\,\,\,
(PGV<4cm/s)\approx(MMI<4.4)\\
2.5+3.155\log_{10}(PGV)\,\,\,\,\,
(PGV\ge4cm/s)\approx(MMI\ge4.4)\\
\end{align}\right\}

PGV=\left\{\begin{align}
10^{((MMI-3.5)/1.5)}\,\,\,\,\,
(MMI<4.4)\approx(PGV<4cm/s)\\
10^{((MMI-2.5)/3.155)}\,\,\,\,\,
(MMI\ge4.4)\approx(PGV\ge4cm/s)\\
\end{align}\right\}

となります。目視で求めているのでlogへの変換はこれくらいで十分でしょう。
この式は他の地震も適当に混ぜたものでもいい感じに一致しました。
ただし、この画像から逆算するとMMIが少しずれます。(この表がよくわかっていないのでとりあえずこのまま)

青:各データ　灰:PGV<4cm/s　黄色:PGV≧4cm/s(少し範囲外も描画しています)
順番に計算してみるとこのようになります。

(ln)はMMI→PGAの計算に自然対数、(log)はMMI→PGAの計算に常用対数を使ったものです。
もう少し細かくすると、

ここで、MMIと気象庁震度階級の揺れの大きさの対応表を比べて見ると思います。
4<MMI<6のあたりは式が変わるところなのでこの辺りは曲線になると良い…のかな？
一応多項式近似を出してみましたが大きな違いはなさそうですね。どこが誤差が多いのかわからないので直接の式を出すのはやめときます。

コードの例

private void Main()
{
  double MMI = 5.678;
  double PGV = MMI2PGV(MMI);
  double JI = PGV2JI(PGV);
  Console.WriteLine("MMI:" + MMI);
  Console.WriteLine("PGV:" + PGV);
  Console.WriteLine("震度:" + JI);
}

public double MMI2PGV(double MMI)//MMI→PGV
{
  if (MMI < 4.4)
    return Math.Exp((MMI - 3.5) / 0.65);
  else
    return Math.Exp((MMI - 2.5) / 1.37);
}

public double PGV2JI(double PGV)//PGV→JmaIntensity
{
  if (PGV < 7)
    return 2.165 + 2.262 * Math.Log10(PGV);
  else
    return 2.002 + 2.603 * Math.Log10(PGV) - 0.213 * Math.Pow(Math.Log10(PGV), 2);
}

VBA

Function MMI2PGV(MMI As Double) As Double
  If MMI < 4.4 Then
      MMI2PGV = Exp((MMI - 3.5) / 0.65)
  Else
      MMI2PGV = Exp((MMI - 2.5) / 1.37)
  End If
End Function

Function PGV2JI(PGV As Double) As Double
  If PGV < 7 Then
    PGV2JI = 2.165 + 2.262 * WorksheetFunction.Log10(PGV)
  Else
    PGV2JI = 2.002 + 2.603 * WorksheetFunction.Log10(PGV) - 0.213 * WorksheetFunction.Log10(PGV) ^ 2
  End If
End Function

Function MMI2JI(MMI As Double) As Double
  MMI2JI = PGV2JI(MMI2PGV(MMI))
End Function

日本の地震で比較

福島県沖地震(2022年)

気象庁 USGS

右は計算した震度をマップに表示するものです。いつか公開予定です。配色はKiwi Monitor カラースキーム第2版を少し変更したものです。震度0.5未満も描画しているので注意してください。
↓配色(基本背景(ここでは不使用)、震度0、震度1....)

震度が低いところは特に1階級くらい差ありますがほとんど同じですね。

東北地方太平洋沖地震

気象庁 USGS

同じ感じですが、震度4が多くなっています。

熊本地震(16日)

気象庁 USGS

震度7を観測したところがUSGSにはないですね(緑丸のあたり？)

北海道胆振東部地震

気象庁 USGS

計算震度7はUSGSによるとおそらくHKD126 鵡川ですね。気象庁観測点ではない？
MMI高いと近いMMIでもPGVの範囲が広いから誤差が大きくなりやすいですね。(その分気象庁震度階級も変わりにくいですが)
全体的に計算PGVが高いですが、震度にはほとんど影響してないように見えます。StationListのMMIの計算方法が世界共通かMMIが計算で出したものでないなら、被害が外国より小さくなりやすい日本ではMMIが小さくなるはずなので、その場合なら妥当ではないでしょうか。

やりたかった検証

USGSのStationListは大きな地震の場合防災科研のデータがあります(要検証)が、実際の詳しい震度(震度5.3、震度5.7など)が分かれば詳しく比較できますがどこにあるかわからないのでパス。(気象庁は防災科研のデータは詳しく公開していない) そういえばわかるやん。まあ気が向いたら()
DYFIはアンケートを基に推定(メルカリ)震度が出ますが、ついでにPGAとPGVも出ます。この値が震度から計算されたのかアンケートも考慮しているのかただ推定震度分布的なものから出しているのかはわかりませんが、これも何かに生かせそうな...?

実践(USGSの推定最大MMIから)

USGSのShakeMapにはIIIとかVとかVIIとかしか表示されませんが、Feedにはもっと詳しく書かれています。 こっちを使いましょう。
WorldQuakeViewerには詳細なMMIを表示しているので発生して数時間以内だったらどうぞ。
↓は面倒なので簡単に見れるページ作っておきます。

json(geojson)を見ることになるので、JsonFormatterなどで検索して拡張機能を入れたりサイトで見やすくすることをお勧めします。私はこれを使用しています。↓では説明用に[]を用いていますがその部分は実際の値です。
USGSの地震の詳細ページのURLとかからIDを探してコピーしておきましょう。(URLだったらhttps://earthquake.usgs.gov/earthquakes/eventpage/[ここがID]/executive)

Feedのフォーマットはここで見れます。

1.地震ごとのFeedを見る(こっちのほうが簡単)

https://earthquake.usgs.gov/earthquakes/feed/v1.0/detail/[ID].geojson を開く。
"mmi":[ここがMMI、nullだったら未発表]のところを探す。最初のほうにあります。
計算する

2.総合Feedを見る

https://earthquake.usgs.gov/earthquakes/feed/v1.0/geojson.php の右側にあるFeedの、発生から1日以内ならPast Day、7日以内ならPast 7 DaysのSignificant Earthquakes(重大な地震(M6クラス、MMI5クラス、アメリカならM3~4でもでも入る))かM4.5+ Earthquakes(M4.5以上)を開く。
該当の地震を頑張って探すか、IDで検索する
{"type": "Feature",～"id": "[ここにID]"}, が1つの地震です。
上の塊の中で、"mmi":[ここがMMI、nullだったら未発表]のところを探す。
計算する。

USGSのShakeMapのMMI-0.5~+0.5での計算震度と観測震度の差を比べてみるとShakeMapの方は小さいことがあります。このページで近年のMMI5以上の地震が調べられます。(検索条件を確認してください)

応用

StationListのデータを利用してこれのように気象庁震度マップを作ることができます。(現在作成中です。)

結論・感想

一応求めることができます。
USGSのShakeMapは推定値なのでマグニチュードと深さ的に明らかに小さければ(日本の地震は特に小さくなりやすい)マグニチュードと深さで震度予想したほうが正確なこともあります。
被害を推測するならPAGERを参考にしましょう。
私が出した式はある程度の精度はありますがまだまだ改良の余地があるのでぜひ暇があれば求めてみてください。
~~4日もかかった。疲れた。~~

参考・データ使用

気象庁
USGS
各先行研究?

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up