起震車の揺れをiPhoneで計測して震度を計算

iPhoneは3軸の加速度計を内蔵しています。iPhoneの本体を傾けると画面が縦モードになったり横モードになったりしますが、あれは加速度計で重力方向を検出して、向きを検知しているわけです。しかしせっかく加速度計があるので、もっといろいろな応用、例えば震度計としても使えるのではないでしょうか。そこで、夏休みの自由研究として、iPhoneの加速度計で地震の震度を計ってみました。そして、気象庁がどうやって震度を決めているのかも調べました。

とは言うものの、地震が来るまで待ってはいられません。運が良いことに、防災行事で起震車に乗る機会があったので、iPhoneのSensorLogというアプリで地震発生から終了までの加速度を計測してみました。震度7の揺れです。

2回に分けて大きな揺れが発生したのがわかります。このデータを元にして震度を求めてみます。震度と言えば、昔は職員の体感で震度を決めていたのですが、最近は機械で判定するようになったようで、そのアルゴリズムが気象庁のホームページで公開されています。「最大加速度値を何かの関数に放り込 んで震度を出すのだろう」ぐらいに思っていたのですが、意外に込み入ったアルゴリズムが使われていました。その手順通りに計算を進めてみます。

１　ディジタル加速度記録３成分（水平動２成分、上下動１成分）のそれぞれのフーリエ変換を求める。

FFTにかけてみました。X,Y軸とZ軸で最大成分の周波数が違うのが興味深いです。

２　地震波の周期による影響を補正するフィルターを掛ける。

ハイパス・ローパスに加えて、長周期の影響を弱めるようなフィルターを掛けます。フィルターの具体的な式は気象庁のページに載っていなかったので、こちらの『第１３講　計測震度』を参考にしました。

３　逆フーリエ変換を行い、時刻歴の波形にもどす。

フィルターを掛ける前（青）と後（オレンジ）のグラフを、揺れが激しかった時間帯について表示しました。特にローパスフィルターが効いていて、平滑化されているのがわかります。

４　得られたフィルター処理済みの３成分の波形をベクトル的に合成をする。

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

この意味は、フィルターを掛けた後の波形の最大加速度をaとするということです。ただし、外れ値（極端値）を拾わないように、全体のうち高い方から0.3秒分は除いているのかと思います。

波形の絶対値についてソートした結果はこうなります。0.3秒（赤線）のところを見ると582.2Galで、これがaの値です。

６　５で求めたaを、I = 2 log a + 0.94　により計測震度 Ｉ を計算する。計算されたIの小数第３位を四捨五入し、小数第２位を切り捨てたものを計測震度とする。

計測震度ⅠはI=2log 582.2+0.94=6.47014です。

これを四捨五入すればそのまま震度階級になるかと言えばそうではなく、微妙なトラップがあります。「Iの小数第３位を四捨五入し、小数第２位を切り捨てたものを計測震度とする」なので、例えば計測震度Ⅰが3.497の場合は計測震度は3.5、これを四捨五入して震度階級は4となります。表にまとめるとこんな感じかなと思います。

おそらく、発表されるのは計測震度と震度階級なのでしょう。計測震度が同じなのに違う震度階級が発表されるという事態を避けるために、このような手順が決められているのかと思います。

そういうわけで、今回iPhoneで計測した加速度値による計測震度Ⅰは6.47014、計測震度は6.4、震度階級は6となります。7には届きませんでした。計測震度Ⅰがあと0.02486だけ大きければ震度階級7だったのでで、惜しかったとも言えます。そもそも起震車が正しく震度7の揺れを起こしていたのかという問題もあるため、結論はよくわからないということにしておきたいと思います。

ところで、最大加速度を元に計算するという想像は、あながち外れてはいませんでした。ただし、

• 波形をそのままではなく、周波数フィルターに通して使う。
• 外れ値は除く。

という点に注意は必要です。

計測震度Ⅰは2log a+0.94なので、a^2が10倍になる毎に1増えます。加速度の2乗はその地点の地面が持つエネルギーに比例するので、計測震度Ⅰはエネルギーの相対的なベル値（デシベル値ではなく）を計ったことになります。おおざっぱに言えば、震度とは局所的なエネルギーの『10を底とした対数値』であるということのようです。