The Negligible Lab

Astray in the forest of electrical and electronic circuits. Adrift in the gap between time and frequency domains. 独立独歩を是とする。

Raspberry Piで地震計を作る ① 震度計算編

Raspberry Pi Python ディジタル信号処理

はじめに

世界の地震の10%以上が日本で発生している──。巷間よく言われていることですが，先日も東京23区と埼玉県で10年ぶりに震度5強を観測したばかりであり，地震について考えない日はないと言っても過言ではないでしょう。

地震が発生した際に，気象庁のホームページやNHKなどで各地の震度が発表されます。しかし，これを待たずして自作の地震計で震度が分かったら，さらに，震度が分かるだけでなく加速度の波形も記録できたら面白そうです。昨今は加速度センサが簡単に手に入るので*1，地震計というのは電子工作やプログラミングの練習にもちょうどいい題材になるのではないでしょうか？

今回，次回，次々回を合わせて，筆者がマイ地震計を作ろうと思い至った経緯や，開発にあたっての課題，Raspberry Piでの実現方法について書いていきます。今回はその第一弾として，加速度の時系列データから計測震度を求める計算方法について書いていこうと思います。

加速度センサから定周期でデータを読み出しながら震度計算するといったリアルタイム処理については次回に，加速度センサとそのI2C接続，プログラムの全体構成については次々回に譲ります。

観測例

いきなりですが，実際に起きた地震を捉えた観測例を示します。図1に，試作中のマイ地震計プロトタイプを示します。マルチプロセッシングを活用するので，4コアのCPUを持つRaspberry Pi 3 Model A+を使うことにし，ブレッドボード上にInvenSense社のモーションセンサMPU6050とOLEDディスプレイを置きました。いずれもI2CバスにてRaspberry Piに接続されますので，Raspberry PiのGPIO2 (SDA)とGPIO3 (SCL)に並列に接続しました（もちろん電源も必要です）。ジャンパワイヤの配線はいたってシンプルです。［追記: 2021-11-19］同じI2CバスにモーションセンサMPU6050とOLEDディスプレイを繋げると，加速度の定期的な読み出しのサイクルが乱される可能性があることが判りました。［追記: 2022-03-15］同じI2Cバスに接続してもOLEDディスプレイのためのフレームの隙間にMPU6050のためのフレームがほぼ10 ms周期で挿入されており，ほぼ問題ないという結論に至りました。詳しくは次々回の付録をご覧下さい。

f:id:s-inoue2010:20211115185649j:plain:w480
図1: Raspberry Pi 3 Model A+とモーションセンサMPU6050によるマイ地震計プロトタイプ

試作中のPythonプログラムは10 ms毎に加速度をサンプリングしつつ，加速度データが300点集まった時点で（つまり，3秒毎に）短時間計測震度を計算します。

ある夜，寝る前にプログラムを走らせておいたところ，早朝に発生した地震による加速度が記録されており，また，加速度から震度も（単に「それらしい値」ですが）計算されておりました。図2に結果を示します。この地震でマイ地震計プロトタイプが計算した計測震度の最大値は2.8であり，震度階級では震度3に相当します。

f:id:s-inoue2010:20211116083500p:plain
図2: 実際の地震の観測例（2021年11月8日03時08分頃の地震）

これは，2021年11月8日03時08分頃に発生した茨城県南部を震源とするマグニチュード4.3の地震です*2。気象庁のデータによれば，筆者の住む街では震度2を記録しています。一方，マイ地震計プロトタイプでは震度3を計測しました。恐らくですが，2階に置いてあったため，家屋での共振によって大きめの加速度を観測したのではないかと推測します。

なお，筆者のマイ地震計プロトタイプは，震度1以下の地震では反応しないようにしました。と言うのも，次々回で詳述しますが，モーションセンサMPU6050のノイズが大きいため，全く揺れていない状態と震度1を区別できないためです。より高精度の加速度センサが手に入れば，震度1以下も測れるのではないかと思います。

目次

開発経緯

数年前，Mbedで遊び始めた頃，筆者は地震計の自作を試みました。しかし，加速度から震度を求める計算に離散フーリエ変換(Discrete Fourier Transform, DFT)または高速フーリエ変換(Fast Fourier Transform, FFT)をはじめとする（当時の筆者にとっては）複雑な計算が必要と判り，いったん挫折しておりました。その後，私生活が多忙となり*3，趣味の電子工作はまったくできなくなってしまったため，休眠プロジェクトとなっていました。

復活のきっかけとして背中を押してくれたのは話題の地震計キット「PiDAS」です。

k-s.booth.pm

PiDASはRaspberry Pi Picoを使ってリアルタイムに震度を計算・表示してくれるようです。これは面白い🎵 しかしやはり，技術者の端くれとして自分で作りたいではありませんか…❗

前々回，前回でDFTやFFTのアルゴリズムに触れてみました。Wikipediaの震度計算に関する記述や気象庁の資料を今になって改めて眺めてみると，確かに複雑ではありますが，プログラムとして作れないこともなさそうだという印象を持ちました*4。そこで，まずは震度計算のアルゴリズムをJupyterLabで確認しながら実装しました。気象庁が過去の大きな地震のデータをCSV形式で公開していますので，これを使って加速度と震度の関係を照らし合わせることで，震度を正しく計算できているか確認できます。

ふと部屋に転がっているRaspberry Piを見ると，これはPython (+ NumPy)が動作し，そして加速度センサも繋げられるLinuxコンピュータです。つまり，加速度センサからデータを取得して，それをJupyterLabで作った震度計算プログラムに食わせれば，よもやよもや，意外と簡単に地震計を作れるのではないかと考えました*5。

開発課題

マイ地震計を作るにあたり，地震がひと通り終わってから震度を計算するのではなく，防災科学技術研究所（防災技研/NIED）の「強震モニタ」のように，ある程度のリアルタイム性をもって短時間毎の震度を表示したいと考えました（図2の通り，3秒毎に計算することにしました）。前掲のPiDASも強震モニタと同様の仕組みを実現しているそうです。

www.kmoni.bosai.go.jp

この場合，地震が継続している，つまり，加速度センサによる加速度のサンプリングが続いている最中に，それまでに取得した数秒間の加速度データをもとに震度計算のルーチンを走らせなければいけません。気象庁のデータを見ると，加速度のサンプリング周期は10 msです。10 ms周期のサンプリングという定周期性を守りながら，震度計算も同時に行う──。これを実現するには，マルチプロセッシングを活用し，加速度をサンプリングするプロセスと震度計算するプロセスに分け，かつ，前者ではLinuxのシステムコールを使ったタイマ割り込みを作れば良いのではないかと考えました。

筆者としては経験のないことが多いので，問題を切り分け，

  1. 震度計算プログラムの実装
  2. システムコールsetitimerを活用したタイマ割り込み
  3. マルチプロセッシングとキュー
  4. タイマ割り込みとマルチプロセッシングの組み合わせ
  5. 加速度センサからの加速度取得
  6. OLEDディスプレイへの表示
  7. マイ地震計全体のプログラムの組み立て

という順序で組み上げていくことにしました。なお，今回は上記1についてのみ取り扱い，2 ~ 7については次回，次々回に譲ります。

計測震度の計算方法

観測された加速度から震度を計算する手順は，（私のような凡人にとっては）意外なほどに複雑です。ここでは，震度計算について詳述した後，Pythonで実装していきます。

用語について

普段，無意識に「震度」という語を使っていますが，詳しくは表1のように「計測震度」と「震度階級」に分けられます。なお，表中の定義は筆者の理解であり，気象庁の公式の定義ではありませんのでご注意ください。

表1: 計測震度と震度階級

用語筆者の理解
計測震度
(JMA instrumental seismic intensity)		 計測した加速度から一定の手順によって算出される数値であり，小数点以下第1位までを有する。
震度階級
(JMA seismic intensity scale)		 計測震度の範囲に割り当てられたレベルであり，震度0，1，2，3，4，5弱，5強，6弱，6強，7の10段階がある。

本記事では，基本姿勢として計測震度と震度階級を総称して単に「震度」と呼ぶことにしますが，混乱の恐れがなければどちらも単に「震度」と呼ぶことも許容します。以下では加速度から計測震度を求める手順について説明します。

加速度から計測震度を求める手順

気象庁の下記ページにて，加速度から計測震度を求める手順が概説されています。

www.data.jma.go.jp

ごく簡単に書けば，

  1. 南北，東西，上下方向の各加速度成分（単位[gal] = [cm/s2]）をフーリエ変換し，周波数領域に移す
  2. 各成分に対して地震波の周期による影響を補正するフィルタ（後述）を掛けた後，逆フーリエ変換で時間領域に戻す
  3. 時間領域に戻した3つの成分の合成加速度の絶対値がある値a以上となっている時間がちょうど0.3秒となるようなaを求める
  4. 上記で求めたaから，I = 2 log10 a + 0.94としてIという値を求める
  5. Iの小数点以下第3位を四捨五入，小数点以下第2位を切り捨てる

という手順で計算するようです。

上記2における「地震波の周期による影響を補正するフィルタ」については，同じく気象庁の下記資料のII-35ページ以降，また，Wikipediaの「気象庁震度階級」の項目にも記載されております。

https://www.data.jma.go.jp/svd/eqev/data/study-panel/shindo-kentokai/hensen.pdf

ja.wikipedia.org

上記のフィルタは，周期効果フィルタ，ハイカットフィルタ，ローカットフィルタの3つから構成されます。なお，「フィルタを掛ける」，またWikipediaの記事では「畳み込む」と大仰な言葉を使っていますが，周波数領域では単に掛け算するだけです。しかもこのフィルタはゲイン特性のみを持っており，位相特性はありません（伝達関数が実数のみ）。以下，それぞれのフィルタについて述べます。

3つのフィルタ

周期効果フィルタ

周波数領域にて，加速度のスペクトルに次式の重みを乗算します。fは周波数です。

W_{PE} (f) = \displaystyle \sqrt{\frac{1}{f}} \tag{1}
ハイカットフィルタ

周波数領域にて，加速度のスペクトルに次式の重みを乗算します。

W_{HC} (f) = \displaystyle \frac{1}{\sqrt{P(x)}} \tag{2}

ただし，

\begin{align} P(x) &= 1 + 0.694 x^{2} + 0.241 x^{4} + 0.0557 x^{6} \\ &+ 0.009664 x^{8} + 0.00134 x^{10} + 0.000155 x^{12} \end{align} \tag{3}

また，

x = \displaystyle \frac{f}{10} \tag{4}

です*6。

ローカットフィルタ

周波数領域にて，加速度のスペクトルに次式の重みを乗算します。

W_{LC}(f) = \displaystyle \sqrt{1 - \exp \left \{- \left ( \frac{f}{0.5} \right )^{3} \right \}} \tag{5}
3つのフィルタの特性

図3に，上記3つのフィルタの周波数特性をプロットします。

f:id:s-inoue2010:20211115165752p:plain
図3: 周期効果，ハイカット，ローカットの各フィルタの周波数特性

3つのフィルタを合成した

W(f) = W_{PE}(f) \cdot W_{HC}(f) \cdot W_{LC}(f) \tag{6}

は，図3の黒色のカーブとなり，低域と高域がカットされている他，0.7 ~ 7 Hz辺りにおいても，緩やかに高域が低下していく特性となっております。

筆者としてはこのフィルタの物理的意味はまったく理解していませんが，Python-NumPyでこれを実現することは容易に思えてきました。

aという値を求める

南北，東西，上下方向の加速度に対して，周波数領域で(6)式のフィルタを掛けた後，逆フーリエ変換にて時間領域に戻します。その後，3つの成分をベクトルとして合成した後，気象庁のホームページに記載の通り，

ベクトル波形の絶対値がある値a以上となる時間の合計を計算したとき，これがちょうど0.3秒となるようなaを求める。

という処理を行います。少し分かり難いので図4に概略図を描いてみました。

f:id:s-inoue2010:20211115184615p:plain:w600
図4: 加速度がa以上となる時間が0.3秒となるようなaを求める

ベクトルとして合成した加速度の絶対値は必ず正になります。この加速度がある値a [gal]以上となっている期間は，図4の例ではT1，T2，T3，T4の4か所となっています。ある値aを正しく見つけた場合，その和TはT = T1 + T2 + T3 + T4 = 0.3 sとなります。

図4のように，合成加速度の絶対値の波形からこのaをという値を見つける必要があります。aをものすごく大きな値にすると，合成加速度の絶対値がa以上となる時間帯はまったくない（ゼロ秒）でしょう。aを少しずつ小さくしていくと，どこかで合成加速度の絶対値の波形のピークにぶつかり，それ以降はaを小さくしていくと，合成加速度の絶対値がa以上となる時間は少しずつ長くなっていくことでしょう（直観的には単純増加であるはず…）。したがって，ソート済みのデータに対して目的の値を探索する二分探索（バイナリサーチ）を応用してaの値を見つけることができます。

計測震度と震度階級

aの値が求まったら，次式でIという値を求めます。

I = 2 \log_{10} a + 0.94 \tag{7}

Iの小数点以下第3位を四捨五入，小数点以下第2位を切り捨てることによって，計測震度を得ます。混乱の恐れがなければ，四捨五入や切り捨ての後の値もIと呼ぶことにします。

Iは小数点以下第1位までの値を持つ数になります。ここから表2に従っていわゆる震度階級を決定します。

表2: 計測震度と震度

計測震度震度階級
0.5未満0
0.5以上1.5未満1
1.5以上2.5未満2
2.5以上3.5未満3
3.5以上4.5未満4
4.5以上5.0未満5弱
5.0以上5.5未満5強
5.5以上6.0未満6弱
6.0以上6.5未満6強
6.5以上7

以上の手順で，加速度から計測震度および震度階級を求めることができます。

Pythonでの実装

shindoモジュール

以上で述べた手順を元に，計測震度を求めるgetShindo関数と，そこから表2に基づいて震度階級を得るgetShindoName関数を作り，これをshindo.pyというファイルに格納しました。他の.pyファイルからモジュールとしてimportできます。GitHubに置きましたので，ご興味のある方は試してみて下さい*7。

github.com

全文は長いので，getShindo関数に関わる部分を以下に抜粋します*8。

import numpy as np

def _filter(A: np.ndarray, Ts: float) -> None:
    """
    @brief Apply filter to the accleration spectra
    @param A 3-D acceleration frequency-domain spectra, N-S, E-W, and U-D
    @param Ts Sampling period
    """
    N = len(A)
    k = np.arange(N)
    f = k / (N * Ts * 2)

    # Periodic-effect filter
    epsilon = 0.0001  # To prevent division by zero
    W_pe = np.sqrt(1 / (f + epsilon))

    # High-cut filter
    x = f / 10
    W_hc = 1 / np.sqrt( \
             1 + 0.694 * x**2 + 0.241 * x**4 + 0.0557 * x**6 \
             + 0.009664 * x**8 + 0.00134 * x**10 + 0.000155 * x**12 \
    )

    # Low-cut filter
    W_lc = np.sqrt(1 - np.exp(-(f / 0.5)**3))

    # Apply filter
    A[:,0] *= (W_pe * W_hc * W_lc)
    A[:,1] *= (W_pe * W_hc * W_lc)
    A[:,2] *= (W_pe * W_hc * W_lc)

def _search_aval(a: np.ndarray, Ts: float) -> float:
    """
    @brief Search for the a value
    @param a 3-D acceleration time-domain data in [gal], N-S, E-W, and U-D
    @param Ts Sampling period
    @return The a value found
    """
    aval = 2000.0            # Initial value of search [gal]
    T_ref = 0.3              # Time where acceleration is above the a value
    epsilon = T_ref * 0.001  # Acceptable error
    while True:
        T_above_aval = np.count_nonzero(a >= aval) * Ts

        # Too high
        if T_above_aval < T_ref - epsilon:
            aval -= aval / 2
            continue

        # Too low
        if T_above_aval > T_ref + epsilon:
            aval += aval / 2
            continue

        # The a value found
        break
    return aval

def getShindo(a: np.ndarray, Ts: float) -> float:
    """
    @brief Calculates JMA shindo scale from acceleration data as ndarray 
    @param a 3-D acceleration time-domain data in [gal], N-S, E-W, and U-D
    @param Ts Sampling period
    @return Calculated instrumental shindo scale
    """

    # Perform FFT
    A = np.fft.rfft(a, axis = 0)

    # Apply filter defined by JMA
    _filter(A, Ts)

    # Perform inverse FFT
    afil = np.fft.irfft(A, axis = 0)
    afil_total = np.sqrt(np.sum(afil**2, axis = 1))

    # Search for the a value
    aval = _search_aval(afil_total, Ts)

    # Calculate JMA instrumental seismic intensity
    I_raw = 2 * np.log10(aval) + 0.94
    I = np.floor(np.round(I_raw, decimals = 2) * 10) / 10

    return I

N点の3軸加速度（単位[gal]）のデータがNumPyのN行3列の配列aに格納されているものとします。配列aをサンプリング周期TsとともにgetShindo関数に渡せば，戻り値としてfloat型の計測震度を返します。また，getShindo関数から呼ばれることを意図した_filter関数があり，これは加速度に前述の3つのフィルタを掛ける働きをします。同様に_serarch_aval関数があり，これは前述のaの値を二分探索（バイナリサーチ）で求めます。

テストベンチ

shindoモジュールでは，いわゆるif __name__ == '__main__':以下の部分にテストベンチを作っております。前述のように気象庁では過去の大きな地震について，加速度や震度のデータを波形画像およびCSVファイルとして公開しております。このデータを活用すれば，震度計算プログラムの計算結果を検証できます。

www.data.jma.go.jp

気象庁の前掲のページでは，2000年10月6日に発生した鳥取県西部地震における米子市での観測記録（計測震度5.1）を例に，加速度から震度を求める方法を説明しています。この地震では米子市で計測震度5.1を観測しています。

https://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/001006_tottori-seibu/png/AA06EA01.png
図5: 鳥取県西部地震における米子市での観測記録（出典: 気象庁ホームページ）

図5に，気象庁の「強震波形（平成12年(2000年)鳥取県西部地震）」から引用した波形を示します。CSVとして数値データも公開されており，南北，東西，上下方向の加速度が10 ms毎に記録されています。データ点数は30,000点（5分間）です。

NumPyはCSVファイルを読み取ることもできます。これを使って，気象庁のCSVファイルから3軸の加速度を読み込み，前述のgetShindo関数に食わせてみます。図6のように，Raspberry Piにてこのテストベンチを実行してみると，気象庁の値と同じく計測震度5.1（震度5強）を得ることができました。

f:id:s-inoue2010:20211116094015p:plain
図6: shindo.pyのテストベンチをRaspberry Piで実行した様子

う～む，30,000点のデータに163 msも掛かっていますね*9。

短時間計測震度

以上で計算した計測震度と震度階級は，加速度の波形全体に対するものでした。一方，図5の左下に計測震度の時系列グラフが参考値としてプロットされております。これは，過去5秒間の加速度から短時間の震度を計算したものです。

マイ地震計プロトタイプでも，このように，地震が継続している最中での震度の変化を計算・表示したいと考えております。しかし，10 ms毎のそれぞれのサンプリングの後で，過去数秒間の加速度から移動平均のように計測震度を計算するのは恐らく計算パワーの観点から難しいと思われるので，数秒毎に飛び飛びのタイミングで短時間の計測震度を計算することにしました。先ほどの鳥取県米子市での波形を例とすれば，図7のようなイメージになります。

f:id:s-inoue2010:20211116100049p:plain
図7: 短時間計測震度（5秒毎）の計算例

これは，気象庁のCSVファイルに記録されている加速度を500点毎にスライスしてgetShindo関数に与えたものです。図5の左下のグラフのように滑らかではありませんが，はじめの10秒ほどで地震が強くなり，そして3つの小さな山を作りながら，徐々に収束していく様子が見て取れます。ただし，短時間計測震度の最大値が5.0である一方，波形全体から計算した計測震度は5.1です。若干の誤差は許容する必要がありそうです。

この短時間震度の「窓」は自由に決められますが，短かすぎると波形全体から計算した計測震度との誤差が大きくなり，長すぎると，短時間計測震度がまばらなタイミングしか算出されません。マイ地震計プロトタイプでは，短時間計測震度の窓を前述の通り3秒間としました。

まとめ

以上で，加速度の時系列データから計測震度を計算する方法について述べ，これをRaspberry Pi上のPythonプログラムとして実装した結果についてまとめました。また，数秒間のデータをスライスとして取り出して短時間計測震度を計算した例を示しました。

次回，次々回は，一定時間(10 ms)毎に加速度をサンプリングしながら短時間震度計算を実行するリアルタイム処理について書き，マイ地震計プロトタイプに実装していきます。

*1:後から分かりましたが，地震計として本格的に使うためには，1 galを下回るような小さな加速度を精度よく測れる必要があり，趣味の電子工作で手に入るような加速度センサでは，本来は対応できないと考えます。

*2:https://www.data.jma.go.jp/multi/quake/quake_detail.html?eventID=20211108031241&lang=jp

*3:お察し下さい。

*4:どうしてこんなに複雑な計算になったのか（あるいはこの程度の複雑さで済んでいるのか）について，残念ながら筆者は理解していません。

*5:筆者はまだRaspberry Pi Picoの使い方が分かっていません。1枚所有してはおりますが…。

*6:なんでこんなまどろっこしい式なのかは不明です。

*7:型ヒントとかDoxygenのコメントのマネゴトを書いていますが，きちんと機能するのか分かりません…💦

*8:FFTですか？ 自分では作らず，numpy.fft.rfft関数（実数信号に特化したFFT）に頼ります😅

*9:次節の短時間計測震度の場合はどうだろうと考え，データ点数を最初の300点に絞ってみると，2.5 msでした。一応，毎回のサンプリングの後に移動平均のように短時間計測震度を計算することはできそうですね…。