2018年9月22日 | ブログ記事

数理計算の講義を活かして震度計を作った話

RLook

本記事は夏休みブログリレー第-3日目です。
こんにちは。数理・計算科学系に所属しているRLookと申します。本記事では表題の通り震度計を作ったので,それについて書いてみたいと思います。

経緯

本学の理系教養科目(1年生向けの科目)の1つに宇宙地球科学基礎ラボ(地球物理)というものがあります。この科目では地球物理に関する実験を3つ行うものです。
そのうちの1つとして,地震発生のメカニズムとそのモデル実験があり,その中でマグニチュードの性質を学びました。
ところで,地震の規模を示す数値としてマグニチュードの他に震度というものがあります。マグニチュードは個人で計測することが困難である一方,震度は加速度センサーさえあれば計測することが可能です。
そこで震度計も作ってみようと思ったのですが,後述する通り震度算出の過程において高速フーリエ変換(FFT)が必要で,なかなか手を出せずにいました。
2年生の第2クォータに入り,私の所属している数理・計算科学系にてアルゴリズムとデータ構造という科目が開講されました。この講義は基本的なアルゴリズムを学んで,それをC言語で実装してみようというもので,その中に高速フーリエ変換があったため,震度計作成に挑むことにしました。

震度とは

ところで震度とはどのように算出されるのでしょうか。現在の震度の算出方法は,平成八年二月十五日気象庁告示第四号にて定められています。
そもそも,以前は震度を人の主観で観測していました。しかし,1996年の4月1日以降に発生した地震は前述の方法によって算出されるようになりました。
一般的にニュース等で報じられる震度は0から7までの10段階ですが,実際にはもっと細かな値が算出されています。これを計測震度とよび,これを四捨五入などすることによって10段階に分けられています。
念のためこの計測震度の算出方法を下記に示します。

  1. 地震動の時刻tにおける直交する3成分の加速度をもとめ,それぞれフーリエ変換する。
  2. 地震波の周期による影響を補正する次の3フィルターを掛ける。
    • 周期効果フィルター Ff(f)F_f(f)

      Ff(f)=1fF_f(f) = \sqrt{ \frac{1}{f} }

    • ハイカットフィルター Fh(y)F_h(y)(ただしy=0.1fy = 0.1f)

      Fh(y)=(1+0.694y2+0.241y4+0.0557y6+0.009664y8+0.00134y10+0.000155y12)12F_h(y) = (1+0.694y^2+0.241y^4+0.0557y^6+0.009664y^8+0.00134y^{10}+0.000155y^{12})^{-\frac{1}{2}}

    • ローカットフィルター Fl(f)F_l(f)

      Fl(f)=1exp((2f)3)F_l(f) = \sqrt{1-\exp(-(2f)^3)}

  3. 逆フーリエ変換し,合成加速度の大きさv(t)v(t)を求める。
  4. 次を満たすa0a_0を求める。ただし集合AAについて,χA(x)\chi_A(x)AAの定義関数である。

    a0=min{aχ{xxa}v(t)dt0.3}a_0 = {\rm min} \left\{a \middle| \int \chi_{\{ x \mid x \geq a\}} \circ v (t)dt \geq 0.3 \right\}

  5. 次の式によって求められるIIが計測震度である。

    I=2log(a0)0.94I = 2\log(a_0)+0.94

準備

とりあえず家に転がっていたRaspberry Pi Type Bを震度計専用機とし,そこにAE-LIS3DHをI2Cモードで接続します。

設計

参考資料によれば,気象庁では基本的に60秒地震動を計測しているそうです。また,主な地震の強震観測データをみると,100Hzでサンプリングを取っているようです。fftの関係上(というか私の知識不足により)2の累乗である213=81922^{13}=8192回サンプリングすることにします。
また,それに合わせてLIS3DHのCTRL_REG1の値は0x57にし,ノーマルモードで100Hzで計測することにします。
ちなみにローパワーモードにすると0x28,0x2A,0x2Cが出力されなくなります。気をつけましょう。

実装

(講義のプログラムをそのまま使うために,)C言語で実装しました。
Github: Rlook32/seismometer

講義からは,fftのほかクイックセレクトのコードも使用しました。
このコードを実行すると80秒ほど計測を行い,その結果をもとに計測震度を算出して出力します。

結果

試しに安定している机の上で動かしてみます。
念の為計測中はキーボードの打鍵など振動を生む行為の一切を慎みます。
----------2018-09-22-1.29.33

ぼくの机は常に震度3の揺れを迎えているようです。お疲れ様でした。

まとめ

加速度センサーの精度には限界がありました。震度計がいかに精密にできているかがわかりました。
また,本コードはマニュアルで作動させるものなので,実際に地震が来ても自動では計測してくれません。
さらに,エイリアシングなどを一切考慮しないので50Hz以上の処理に不正がありそうな気がします。
これらの修正は今後の課題として,今回は終えたいと思います。ありがとうございました。

参考

この記事を書いた人
RLook

アイコンはどこかの生徒会書記です。

この記事をシェア

このエントリーをはてなブックマークに追加

関連する記事

2018年9月24日
夏休みブログリレー エンディング
Yosotsu
2018年9月23日
ゲームを作るものを作っています
Sigma1023
2018年9月20日
私のinit.vimを紹介します(+おまけ)
xecua
2018年9月19日
ゲームの敵の動きを実装した話
mds_boy
2018年9月18日
楽して絵っぽいものを作る
Kejun
2018年9月17日
ロシアでSNS規制されたお話
Aniie

活動の紹介

カテゴリ

タグ