M5AtomとRaspberryPiで自宅の震度測定システムを作る

本記事で説明する震度測定システムは、M5AtomLite (MatrixやEcho等でも可)にIMU Unitを取り付けたものと、Raspberry Piとから構成されます。本システムを用いることで、以下のようなことができます。

加速度の生値の記録
計算した震度(80秒毎)の記録
計測した震度のAmbientへの送信(外出先等でも自宅の震度を確認できる)

以前より筆者が住んでいる地域は、周辺の地域に比べ震度が小さく、体感に比べても低めの値が出ているのではないかという気がしていました。
そこで、定量的に自宅の震度を測定できるようなシステムを作ることとしました。

作っているうちに機能を追加したり、カッコいい見た目やきれいなソースコードにしようと色々やっていたため中々出来上がりませんでした。が、先日大きな地震があったので、とりあえず版でもいいのでとにかく早く震度が測定できるシステムを作ることとしました。

気象庁のHPによると、震度は下記の手順で求められるようです。

水平2方向+上下方向、合わせて3方向の加速度を測定する
→気象庁のHPに上がっている加速度のデータは、加速度は100Hzで記録されています。
加速度データをフーリエ変換し、フィルタをかけ、逆フーリエ変換をかける
フィルター処理済みの３成分の波形をベクトル的に合成する
合成した加速度をソートし、上から0.3秒目のデータの加速度をaとしたとき、以下の数式で求められる値Iの小数第3位を四捨五入し、小数第2位を切り捨てる

$I=2log_{10}a+0.94$

例えば、I=5.195だったとき、小数第3位を四捨五入すると5.20となり、計測震度は5.2となる。一方、I=5.194ならば小数第3位を四捨五入すると5.19となり、計測震度は5.1となることに注意が必要です。

本システムの処理の中で難しいポイントは、「100Hzで加速度を測定する」「フーリエ変換・逆変換を行う(それなりに計算処理能力が必要)」ことです。

計算した震度を出力する手段としては、簡単に確認できるようAmbientを使用したいと思います。また作成するシステムは震度を出すだけでなく、後から出力された震度の妥当性を検証できるよう、加速度の測定値も記録したいと思います。

以上を踏まえ、手元にあるデバイスで簡単に実装できる仕組みとして、以下の図の通り実装することとしました。

加速度の取得方法は先日記事にした通りですが、MPU6886が出力する値の単位は[G]であり、震度を測定する際に使用するのは[gal]です。よって以下のように変換する必要があります。

$1G=9.8m/s^2=980cm/s^2=980gal$

以上を踏まえ、M5 Atom Liteに書き込むプログラムは下記の通りです。なお、この記事に記載の「IMU_Unit.h」及び「IMU_Unit.cpp」を追加する必要があります。


main.ino
#include "IMU_Unit.h"
#define SAMPLE_PERIOD 10    // サンプリング間隔(ミリ秒)、100Hzで取得したいので1000/100=10[ms]
float accX = 0, accY = 0, accZ = 0;
IMU_Unit unit1;
void setup()
{
    M5.begin(true, false, true);
    if (unit1.Init() == 0){
        unit1.SetAccelFsr(unit1.AFS_2G);
    }
    else{   
        Serial.printf("Error!\r\n");
    }
}
void loop()
{
    long t = micros();  // 処理開始時刻を記録
    unit1.getAccelData(&accX, &accY, &accZ);
    //MPU6886が出力する値の単位は[G]のため、1G=980galに直して送信する
    Serial.printf("%.2f,%.2f,%.2f\r\n", accX * 980, accY * 980, accZ * 980);
    delayMicroseconds(SAMPLE_PERIOD * 1000 - (micros() - t));//この処理でピッタリ10ms待つ
}

Raspberry Pi側では、以下のプログラムを使用します。なお、Ambientにデータを送信する処理内で、「★各自のチャネルIdに変える★」「★各自のライトキーに変える★」と記載されている場所は、各自の環境に合わせて修正する必要があります。また、Ambientを使用するためには、ライブラリをインストールする必要があり、Raspberry Piのターミナル上で以下のコマンドを実行する必要があります。
「sudo pip install git+https://github.com/AmbientDataInc/ambient-python-lib.git」
このプログラムは、例えば以下のようなコマンドで実行できます(Python3じゃないと動かない？です)。

cd /home/pi・・・/★data_collector.pyのあるディレクトリ名
→プログラムの格納されたディレクトリまで移動する
sudo python3 data_collector.py

Raspberry Piの起動に合わせて自動でこのプログラムを実行するには、以下のようにして「rc.local」でプログラムを実行するよう記載するとよいです。

sudo vi /etc/rc.local
→「/etc」にある「rc.local」を「vi」というツールを使って書き換えます
→書き換えには権限が必要なため、頭に「sudo」をつける必要があります
→「vi」は使い方がちょっと難しいので、詳細は調べてください。
※「i」を入力すると文字の挿入モードになり、次に示すコマンドを入力後、「Escキー」で戻り、「:wq」を入力して「Enterキー」で上書き保存できます。
sudo python3 /home/pi/Desktop/Shindokei/data_collector.py
→上記のコマンドは、「data_collector.py」をデスクトップのフォルダ「Shindokei」に置いた場合のコマンドです。各自の環境に合わせて修正してください。


data_collector.py
# -*- coding: utf-8 -*-
import serial
import numpy as np
import time
import threading
import queue
import datetime
import sys
import os
import shutil
from math import log10 as log
import ambient

#記録するデータの個数をここで指定する
#FFTでの処理がしやすいよう、2のn乗にするとよい
DATA_SIZE=2**13#約80秒ごとにデータを記録する

SAMPLING_FREQ=100#サンプリング周波数は100Hz
D_SAVE=1#加速度データを保存する期間[日]
q = queue.Queue() #キューの作成

##フーリエ変換後に使用するフィルタの係数の計算
freq=np.linspace(0.0,float(SAMPLING_FREQ),DATA_SIZE)#フーリエ変換後の周波数の値を計算
#フィルターの係数を計算
X=freq/10
FH=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))
FL=np.sqrt(1-np.exp(-(freq/0.5)**3))
#直流成分の値でzero divisionが発生してしまうため、nan_to_numを使う
FC=np.nan_to_num(1/np.sqrt(freq))
F=FH*FL*FC
#ナイキスト周波数以降は無視する。代わりに係数を2倍にする。
F[freq>SAMPLING_FREQ/2]=0
F=F*2

#シリアル通信を管理するスレッド
def serial_receiver(port,baud):
   ser = serial.Serial(port, baud, timeout=0.1)
   data = np.zeros((DATA_SIZE, 3)) #データはXYZの3系列
   counter=0
   while True:
   	mes='error'
   	try:
   		mes=ser.readline().decode('utf-8')#シリアル通信でデータを受信
   		#TODO#
   		#誤り検知などの仕組みを検討
   		result = np.array(mes.split(',')).astype(float) #「○,○,○\r\n」の形式でデータを受信するので、numpyのarrayにする
   		data[counter]=result #値を格納
   	except:
   		print(mes)
   	counter+=1
   	if counter==DATA_SIZE: #データがそろったらキューに追加
   		counter=0
   		q.put(data) 
   ser.close()

#データを保存するスレッド
def save_data():
   while(True):
   	try:
   		data_save = q.get()
   		dt_now=datetime.datetime.now()
   		os.makedirs(dt_now.strftime('data/%Y%m%d'), exist_ok=True) #日ごとにフォルダを作成
   		np.savetxt(dt_now.strftime('data/%Y%m%d/%H%M%S_%f.csv'), data_save,delimiter=',', fmt ='%.3f') #小数第3位まで保存
   		shindo=calc_shindo(data_save[:,0],data_save[:,1],data_save[:,2])
   		print(shindo)
   		f=open(dt_now.strftime('data/%Y%m%d/shindo.csv'),'a')
   		f.write(dt_now.strftime('%H%M%S_%f')+',{0:.2f}\n'.format(shindo))
   		f.close()
   		am = ambient.Ambient(★各自のチャネルIdに変える★,'★各自のライトキーに変える★')
   		r = am.send({'d1':shindo})
   	except Exception as e:
   		print(e)
   	time.sleep(0.1)

#古いフォルダを削除する
def delete_folder():
   while True:
   	try:
   		files = os.listdir('data')
   		dt_now=datetime.datetime.now()
   		TIME_NOW=os.path.getctime(dt_now.strftime('data/%Y%m%d')) #最新の日付のフォルダの作成日時をエポック秒で取得
   		
   		for file in files: 
   			t = os.path.getctime('data/'+file)
   			if(abs(TIME_NOW-t)>60*60*24*D_SAVE):#エポック秒でD_SAVE日以上古いフォルダを削除
   				print(abs(TIME_NOW-t),t,TIME_NOW,file)
   				shutil.rmtree('data/'+file)#フォルダごとファイルを削除
   	except:
   		pass
   	time.sleep(10)

#加速度データをFFTしてフィルタかけてIFFTする関数
def calc_filter(data):
   return np.fft.ifft(np.fft.fft(data)*F).real

#震度を計算する関数
def calc_shindo(NS,EW,UD):
   #3軸成分のベクトル合成
   vec=np.sqrt(calc_filter(NS)**2+calc_filter(EW)**2+calc_filter(UD)**2)
   #降順にソート
   s_vec=np.sort(vec)[::-1]
   #震度を計算/加速度の大きいデータから並べて0.3秒(0.3[s]/サンプリング周期：0.01[s]=30番目)のデータで計算
   shindo=2*log(s_vec[int(0.3*SAMPLING_FREQ)])+0.94
   #小数第3位を四捨五入して小数第2位以下を切り捨てる
   return int(round(shindo*1000,-1)/100)/10

if __name__ == '__main__':

   #スレッドを立ち上げる
   thread_1 = threading.Thread(target=serial_receiver, args=('/dev/ttyUSB0', 115200))#raspberry pi上のポート名が違ったら修正する
   thread_2 = threading.Thread(target=save_data)
   thread_3 = threading.Thread(target=delete_folder)
   thread_1.start()
   thread_2.start()
   thread_3.start()

Raspberry Piのプログラムは3つのスレッドから構成され、以下のような動作をします。

シリアル通信を管理するスレッド：シリアル通信で送られてくる加速度データを、numpy array形式の変数に格納し、データが8192個たまったら(81.92秒毎に)震度計算スレッドに送る
震度計算スレッド：加速度を保存/震度を計算して保存/Ambientに震度を送信する
→日付に応じてフォルダを作成し、そこに加速度データのcsvファイルと1日分の震度(「時刻,震度\n」)をまとめたcsvファイルを保存する
ファイル削除スレッド：最終更新時刻から24時間以上経過したフォルダを削除
→フォルダに格納されたcsvデータごと削除します

まず最初に、システムを設置します。本当は加速度センサの取り付け方もしっかり検討したいところですが、とりあえず壁にマスキングテープで固定しました。(本当は床とか壁に、ねじ止めでしっかりと固定したい)

さっそくAmbientに上がってきたデータを見てみます。約80秒ごとにデータが上がってきていることが確認できました。また震度は常に1.7～1.8程度のようです。

次に、記録された生の加速度データを確認します。どうやらY軸が重力方向と一致するよう加速度センサが取り付けられているようです(ちょっとずれていることもわかります・・・)

これをフーリエ変換してフィルタをかけて、逆フーリエ変換したものが下記の通りです。こころなしか高周波成分が減っているように見えます。また、直流成分が消えて平均値が約0galとなっています。なお、XYZ成分を合成した加速度ソートして、大きい方からちょうど30番目(0.3秒目)の値は「2.5gal」で、これをもとに震度を計算すると確かに「震度1.7」となります。
※：自宅が揺れているわけではなく、加速度センサーのノイズに起因するものだと思います。取り付け方の工夫やノイズの少ないセンサーを選定する必要があるのでしょうか？いずれにせよ、気象庁などの地震計測システムで使用されている加速度センサーはきっととてもいいものを使用しているのだろうな、と思いました。

これだけだと何となく信用ならないので、震度の計算処理が正しそうか、気象庁HPの強震観測データをもとに確認してみました。なお、ここではこのページで例として使用されている2000年10月６日に発生した鳥取県西部地震の米子市（計測震度＝5.1）のデータを使用しており、上述の気象庁のHPに示された波形と似た傾向が確かめられました。またこのデータから計算した震度は5.1であり、どうやら震度計算処理はそれなりに正しそうです。

とりあえず版の最低限の機能として、以下の機能を実現できたのではないかと思います。

加速度の記録
震度の計算・記録
震度の時間推移の確認

今後は長時間稼働して問題がないか確認したり、大きな地震が来た時に正確な震度が出せるかなど確認したいと思います。

また、とりあえず版ではない理想版として、今後以下のようなポイントに着手できればと思っています。

M5 AtomとRaspberry Pi間のシリアル通信のロバスト性向上
もっというと無線通信(BLE？)
さらに言うとRaspberry PiなしでM5 Atom単体で震度計算＋加速度記録
※：M5 Atomの処理能力で加速度取得＋フーリエ変換＋震度計算が要求される処理周期で処理できるのか？検討中
※：M5 Atom用のTDカードキット持ってない・・・
※：WiFiの接続関係でトラブったときの処理とか、長時間稼働時の安定性・・・？
加速度センサーの取り付け・固定方法の検討
震度や加速度の確認用インターフェースの工夫
※：Ambientはめちゃくちゃ便利だけど、可能なら加速度の波形も見たい