距離減衰式で長周期地震動階級を推定してみた

この記事は、[防災アプリ開発 Advent Calendar 2023](https://adventar.org/calendars/9301) 8日目、
[木更津高専 Advent Calendar 2023](https://qiita.com/advent-calendar/2023/nit_kisarazu) 8日目の記事です。

# はじめに
:::note warn
一般の人は、気象庁の手順では計算できないです。(予報業務許可を取れば、後述する補正係数と調整値がもらえるかもしれないので頑張ってください)
:::
初参加のマイクラ電鉄です。
某団体で防災情報関連のシステムの開発と維持管理をしています。

# 概要
2023年02月01日から、緊急地震速報の発表基準に長周期地震動階級に基づく基準が追加されたので、距離減衰式で長周期地震動階級の推定をしてみます。

# 計算式
気象庁ホームページの「[地震動の予報業務許可等の申請の手引き](https://www.jma.go.jp/jma/kishou/minkan/tebiki/jishin_tebiki.pdf)」[^1]に詳細が書かれているので読みたい方はどうぞ。
長周期地震動階級を推定するためには、まず「絶対速度応答スペクトル」というのを距離減衰式で推定します。(震度を推定する際にまず最大速度を推定するのと同じ)

### 距離減衰式
```math
log_{10}\,Sva(T)=c(T) + a(T)M - log_{10}\,R - b(T)R + siteFactor(T)
```
$Sva(T)$: 周期 T の絶対速度応答スペクトル［cm/s］
$c(T)$、$a(T)$、$b(T)$: 定数([地震動の予報業務許可等の申請の手引き](https://www.jma.go.jp/jma/kishou/minkan/tebiki/jishin_tebiki.pdf)の別表第2)
$M$: 気象庁(変位)マグニチュード
$R$: 予想地点から震源までの距離［km］
$siteFactor(T)$: 予想地点における補正係数
$T$: 周期(1.6-7.8秒、0.2秒刻み)

ここで、$(T)$は周期ごとであることを表します。また、震度を推定する際に使用する距離減衰式とは違い、気象庁マグニチュードを使い、断層最短距離ではなく震源までの距離を使います。

### 補正係数
気象庁が提供する補正係数がありますが、(おそらく)提供してもらえないので、地盤情報による補正係数を、防災科学技術研究所が提供する[地震ハザードステーション(J-SHIS)](https://www.j-shis.bosai.go.jp/)のS波速度1.3km/s層下面の深さ $D$［m］と、地下 30mまでの表層における平均S波速度 $AVS30$［m/s］を用いて、補正係数$DSC(T)+\varepsilon'(T)$を計算します。

$DSC(T)$: 周期 T における D を用いて約1kmメッシュ毎に求められる補正係数
$\varepsilon'(T)$: 周期 T における AVS30 を用いて約250mメッシュ毎に求められる補正係数

#### DSC(T)の計算
$D \leqq D_0(T)$の場合は、
```math
DSC(T) = k_1(T)
```
$D > D_0(T)$の場合は、
```math
DSC(T) = k_1(T)+ k_2(T) \ log_{10} \ (D/D_0(T))
```

#### ε’(T)の計算
$AVS30 \leqq V_0(T)$の場合は
```math
\varepsilon'(T) = p_1(T) + p_2(T) \ log_{10} \ AVS30
```
$AVS30 > V_0(T)$の場合は、
```math
\varepsilon'(T) = p_1(T) + p_2(T) \ log_{10} \ V_0(T)
```

$D_0(T)$、$k_1(T)$、$k_2(T)$、$V_0(T)$、$p_1(T)$、$p_2(T)$: 定数([地震動の予報業務許可等の申請の手引き](https://www.jma.go.jp/jma/kishou/minkan/tebiki/jishin_tebiki.pdf)の別表第1)

計算が出来たら、$siteFactor(T)$に$DSC(T)+\varepsilon'(T)$を代入します。

### 計算

準備ができたので、1.6秒から7.8秒まで0.2秒刻みで距離減衰式に代入し、$Sva(T)$を計算します。
そしたら、下の表のとおりに長周期地震動階級を出せばよいのですが、
|長周期地震動階級|絶対速度応答スペクトル|
|---|---|
|1|5 cm/s ～ 15 cm/s|
|2|15 cm/s ～ 50 cm/s|
|3|50 cm/s ～ 100 cm/s|
|4|100 cm/s ～|

1秒台、2秒台・・・の長周期地震動階級や、最大の長周期地震動階級を推定するには$Sva(T)$の最大値を計算しようとすると思います(自分は)。しかし、気象庁の手順では、任意の周期帯の$Sva(T)$の最大値を出すには、その周期帯の$Sva(T)$の最大値を出した上で、それに気象庁が提供する調整値を掛けなければいけないので、一般の人で気象庁の手順通りに推定することは無理に近いと思います。(調整する理由については、最大値の系統的なずれを補正するためと思われる)

なので、文献を探し、最大値の系統的なずれを補正しました。[2]
具体的には下のようになります。
```math
log_{10} \ MaxSva' = log_{10} \ MaxSva + 0.1276
```

$MaxSva$: $Sva(T)$の最大値
$MaxSva'$: $MaxSva$に最大値の系統的なずれの補正をしたもの
ただし、補正値の0.1276は、補正係数が$DSC(T)+\varepsilon'(T)$の場合のものです。

この値を長周期地震動階級に変換すれば完成！

# 裏話
実はこれを書いているときに、いままで使用していた長周期地震動階級を推定するプログラムで最大値の系統的なずれを補正していないことが分かって修正しました。

(開発者向けに)
J-SHISから情報を取得する際に、
1つの特定の場所の情報を取得する場合は表層地盤情報提供API、深部地下構造情報提供APIを使えばよいと思いますが、震度観測点毎に情報を取得する際は、ダウンロードページから表層地盤データと深部地盤データをダウンロードして、震度観測点の緯度経度からメッシュコードを計算して取得する方がAPIに負荷を掛けなくてよいかも？

# おまけ
<img width="85%" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479018/9a0fdfb0-5d6f-2137-c933-942846d9c6d8.png">

↑緊急地震速報の推定データ
↓今回の手法で推定したデータ

<img width="85%" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479018/a7ff1c7d-9b5a-4fcd-aa98-2844ab4d8624.png">


# おわりに
最後まで読んでいただきありがとうございました。
初めて記事を書いたので、わかりにくかったらすみません。
誤ったところがあればコメントなどで指摘してください。

また、予報業務許可を持っていない者が震度や長周期地震動階級を、リアルタイムに推定した結果を第三者に提供することは禁止されているので注意してください。

震度の推定に対して長周期地震動階級を推定するにはめんどうになりますが、やってみてはどうでしょうか？

# 文献
1. 気象庁情報基盤部. “地震動の予報業務許可等の申請の手引き”.
気象庁. 2023-11-30. https://www.jma.go.jp/jma/kishou/minkan/tebiki/jishin_tebiki.pdf, (参照2023-12-02)

2. Yadab P. DHAKAL, Wataru SUZUKI, Takashi KUNUGI, Shin AOI, GROUND MOTION PREDICTION EQUATIONS FOR ABSOLUTE VELOCITY RESPONSE SPECTRA (1-10 S) IN JAPAN FOR EARTHQUAKE EARLY WARNING, 日本地震工学会論文集, 2015, 15 巻, 6 号, p.104, https://doi.org/10.5610/jaee.15.6_91

[^1]: 震度の推定や走時表の作り方についても書かれている。