この記事は防災アプリ開発 Advent Calendar 2023 の2日目です。
1日目は、もぐもぐさんの地震速報・監視アプリケーション EQMonitor をアプリストアでstableリリースしましたでした。
3日目は、ingen084 さんの気象庁のコード表を自動で更新するです。

はじめに

おうちで地震観測、ロマンがありますよね。
実際に設置されているような地震計や計測震度計を家に置ければ最高なのですが、金銭面などで色々無理がありますから、現実的な路線としてはマイコンに加速度センサを繋いで揺れを測ってもらうあたりになります。

ここで大事なのは、どの加速度センサを使うかという点です。というのも、加速度センサは価格帯によって分解能や耐ノイズ性といった性能が大きく変化するのです。
今回は、数百円の安価なものから1万円が見えてくるような比較的高価なものまで計5つの加速度センサを用いて揺れを計測し、計測震度の算出結果からそれぞれの性能を比較します。

比較対象

加速度センサを含むセンサは一般的に、加速度データを電圧などのアナログ信号に変換して送るアナログセンサと、SPIやI2Cなどの規格のもとデジタル信号を送るデジタルセンサに大分されます。アナログ信号は配線が受けるノイズの影響を受けやすいほか、結局はどこかでデジタルデータに変換する必要があることを踏まえ、今回は比較対象をデジタルセンサに絞ります。

#	センサ名	メーカー	データビット数	最小分解能	ノイズ密度	販売
1	ADXL345	Analog Devices	10bit	3.9mg	290μg/√Hz（水平軸） 430μg/√Hz（鉛直軸）	秋月電子（450円；終売）樫木総業（396円）スイッチサイエンス（3,669円）
2	LIS3DH	STMicroelectronics	12bit ¹	1mg	220μg/√Hz	秋月電子（600円；終売）マルツオンライン（972円）スイッチサイエンス（1,098円）
3	LSM6DSO²	STMicroelectronics	16bit	0.061mg	70μg/√Hz	スイッチサイエンス（1,998円）³
4	IIS3DHHC	STMicroelectronics	16bit	0.076mg	45μg/√Hz	ストロベリーリナックス（2,420円）
5	ADXL355	Analog Devices	20bit	0.0039mg	25μg/√Hz	ストロベリーリナックス（6,600円）

性能に大きく関わってくるのは最小分解能とノイズ密度です。概ね価格と性能が比例しており、ADXL355に至ってはADXL345に対して最小分解能が1000倍小さく、ノイズ密度が10倍小さくなっています。すごい。

実験

加速度データの収集

前項ではカタログ値を見ましたが、ここでは実際に計測してみましょう。

今回はマイコンにArduino互換の Seeeduino XIAO を使います。秋月では 850円で販売されています。

IIS3DHHCがSPIのみ対応のため、今回は全てのセンサデータをSPIで拾うことにします。配線は以下の通りにします。

Seeeduino XIAO	センサ
3V3	各センサの VCC（あればVSも）
GND	各センサの GND
MISO	各センサの SDO
MOSI	各センサの SDI
SCK	各センサの SCK⁴
D0	ADXL345の CS
D1	LIS3DHの CS
D2	LSM6DSOの CS
D3	IIS3DHHCの CS
D4	ADXL355の CS

Arduino IDEで以下のようなプログラムを作り、Seeeduino XIAOに書き込みます。
（SPIでデータを読み書きする練習も兼ねて作ったプログラムをそのまま載せています。センサのライブラリを使っていないのもあり冗長ですが、ご容赦ください。）

なお、気象庁が運用しているような計測震度計はデータを100Hzでサンプリングしている⁵ことから、同様に100Hzで加速度データを収集します。

#include <TimerTC3.h>
#include <SPI.h>

SPISettings spiSettings0 = SPISettings(10e6, MSBFIRST, SPI_MODE0);
SPISettings spiSettings3 = SPISettings( 5e6, MSBFIRST, SPI_MODE3);
int adxl345Pin  = 0;
int lis3dhPin   = 1;
int lsm6dsoPin  = 2;
int iis3dhhcPin = 3;
int adxl355Pin  = 4;

void setup(void) {
  Serial.begin(115200);
  SPI.begin();

  // タイマーインターバル10ms（100Hz）
  TimerTc3.initialize(10 * 1000);
  TimerTc3.attachInterrupt(timerIsr);

  pinMode(adxl345Pin, OUTPUT);
  digitalWrite(adxl345Pin, HIGH);
  pinMode(lis3dhPin, OUTPUT);
  digitalWrite(lis3dhPin, HIGH);
  pinMode(lsm6dsoPin, OUTPUT);
  digitalWrite(lsm6dsoPin, HIGH);
  pinMode(iis3dhhcPin, OUTPUT);
  digitalWrite(iis3dhhcPin, HIGH);
  pinMode(adxl355Pin, OUTPUT);
  digitalWrite(adxl355Pin, HIGH);

  // センサの設定
  spiWrite(spiSettings3, adxl345Pin, 0x2D, 3, 0x08); //計測モードオン(ODR100Hz, LPF50Hz)

  spiWrite(spiSettings3, lis3dhPin, 0x20, 3, 0x57); //ODR100Hz(LPF不明), 計測モードオン
  spiWrite(spiSettings3, lis3dhPin, 0x23, 3, 0x08); //12bitモード

  spiWrite(spiSettings0, lsm6dsoPin, 0x10, 1, 0x42); //ODR104Hz, LPF任意, 2g, 計測モードオン
  spiWrite(spiSettings0, lsm6dsoPin, 0x17, 1, 0x00); //LPF26Hz

  spiWrite(spiSettings0, iis3dhhcPin, 0x23, 1, 0x40); //LPF235Hz(ODR1.1kHz cosnt)
  spiWrite(spiSettings0, iis3dhhcPin, 0x20, 1, 0xC0); //計測モードオン
  
  spiWrite(spiSettings0, adxl355Pin, 0x28, 2, 0x05); //ODR125Hz, LPF31.25Hz
  spiWrite(spiSettings0, adxl355Pin, 0x2D, 2, 0x00); //計測モードオン
}

void timerIsr() {
  // ADXL345
  byte bs1[6];
  spiRead(spiSettings3, adxl345Pin, 0x32, 3, &bs1, sizeof(bs1));
  int16_t xRaw1 = bs1[0] | (bs1[1] << 8);
  int16_t yRaw1 = bs1[2] | (bs1[3] << 8);
  int16_t zRaw1 = bs1[4] | (bs1[5] << 8);
  float x1 = xRaw1 * 3.9 * 0.001 * 9.80665;
  float y1 = yRaw1 * 3.9 * 0.001 * 9.80665;
  float z1 = zRaw1 * 3.9 * 0.001 * 9.80665;

  // LIS3DH
  byte bs2[6];
  spiRead(spiSettings3, lis3dhPin, 0x28, 3, &bs2, sizeof(bs2));
  int16_t xRaw2 = (bs2[0] >> 4) | (bs2[1] << 4);
  int16_t yRaw2 = (bs2[2] >> 4) | (bs2[3] << 4);
  int16_t zRaw2 = (bs2[4] >> 4) | (bs2[5] << 4);
  if(xRaw2 >= 0x800) xRaw2 = -(~xRaw2 & 0xFFF) + 1;
  if(yRaw2 >= 0x800) yRaw2 = -(~yRaw2 & 0xFFF) + 1;
  if(zRaw2 >= 0x800) zRaw2 = -(~zRaw2 & 0xFFF) + 1;
  float x2 = xRaw2 * 0.001 * 9.806665;
  float y2 = yRaw2 * 0.001 * 9.806665;
  float z2 = zRaw2 * 0.001 * 9.806665;

  // LSM6DSO
  byte bs3[6];
  spiRead(spiSettings0, lsm6dsoPin, 0x28, 1, &bs3, sizeof(bs3));
  int16_t xRaw3 = bs3[0] | (bs3[1] << 8);
  int16_t yRaw3 = bs3[2] | (bs3[3] << 8);
  int16_t zRaw3 = bs3[4] | (bs3[5] << 8);
  float x3 = xRaw3 * 0.061 * 0.001 * 9.80665;
  float y3 = yRaw3 * 0.061 * 0.001 * 9.80665;
  float z3 = zRaw3 * 0.061 * 0.001 * 9.80665;

  // IIS3DHHC
  byte bs4[6];
  spiRead(spiSettings0, iis3dhhcPin, 0x28, 1, &bs4, sizeof(bs4));
  int16_t xRaw4 = bs4[0] | (bs4[1] << 8);
  int16_t yRaw4 = bs4[2] | (bs4[3] << 8);
  int16_t zRaw4 = bs4[4] | (bs4[5] << 8);
  float x4 = xRaw4 * 0.076 * 0.001 * 9.80665;
  float y4 = yRaw4 * 0.076 * 0.001 * 9.80665;
  float z4 = zRaw4 * 0.076 * 0.001 * 9.80665;

  // ADXL355
  byte bs5[9];
  spiRead(spiSettings0, adxl355Pin, 0x08, 2, &bs5, sizeof(bs5));
  long xRaw5 = (bs5[0] << 12) + (bs5[1] << 4) + (bs5[2] >> 4);
  long yRaw5 = (bs5[3] << 12) + (bs5[4] << 4) + (bs5[5] >> 4);
  long zRaw5 = (bs5[6] << 12) + (bs5[7] << 4) + (bs5[8] >> 4);
  if(xRaw5 >= 0x80000) xRaw5 = -(~xRaw5 & 0xFFFFF) + 1;
  if(yRaw5 >= 0x80000) yRaw5 = -(~yRaw5 & 0xFFFFF) + 1;
  if(zRaw5 >= 0x80000) zRaw5 = -(~zRaw5 & 0xFFFFF) + 1;
  float x5 = xRaw5 / 256.0 * 0.001 * 9.80665;
  float y5 = yRaw5 / 256.0 * 0.001 * 9.80665;
  float z5 = zRaw5 / 256.0 * 0.001 * 9.80665;

  Serial.print(String(x1, 5) + " " + String(y1, 5) + " " + String(z1, 5) + " ");
  Serial.print(String(x2, 5) + " " + String(y2, 5) + " " + String(z2, 5) + " ");
  Serial.print(String(x3, 5) + " " + String(y3, 5) + " " + String(z3, 5) + " ");
  Serial.print(String(x4, 5) + " " + String(y4, 5) + " " + String(z4, 5) + " ");
  Serial.print(String(x5, 5) + " " + String(y5, 5) + " " + String(z5, 5) + " ");
  Serial.println();
}

void loop() {
  
}

void spiWrite(SPISettings settings, int pin, byte address, int mode, byte value) {
  spiWrite(settings, pin, address, mode, &value, 1);
}

// モードは今回違う5センサの特徴がいずれかに当てはまるようになっています
void spiWrite(SPISettings settings, int pin, byte address, int mode, void *data, size_t dataSize) {
  SPI.beginTransaction(settings);
  digitalWrite(pin, LOW);
  switch(mode) {
    default:
    case 1:
      SPI.transfer(address);
      break;
    case 2:
      SPI.transfer(address << 1);
      break;
    case 3:
      if (dataSize == 1) {
        SPI.transfer(address);
      } else {
        SPI.transfer(address | 0x40);
      }
      break;
  }
  SPI.transfer(data, dataSize);
  digitalWrite(pin, HIGH);
  SPI.endTransaction(); 
}

byte spiRead(SPISettings settings, int pin, byte address, int mode) {
  byte value;
  spiRead(settings, pin, address, mode, &value, 1);
  return value;
}

void spiRead(SPISettings settings, int pin, byte address, int mode, void *data, size_t dataSize) {
  SPI.beginTransaction(settings);
  digitalWrite(pin, LOW);
  switch(mode) {
    default:
    case 1:
      SPI.transfer(address | 0x80);
      break;
    case 2:
      SPI.transfer((address << 1) | 0x01);
      break;
    case 3:
      if (dataSize == 1) {
        SPI.transfer(address);
      } else {
        SPI.transfer(address | 0xC0);
      }
      break;
  }
  SPI.transfer(data, dataSize);
  digitalWrite(pin, HIGH);
  SPI.endTransaction();
}

パソコンとSeeeduino XIAOを繋ぎ、Arduino IDEのシリアルモニタでデータ出力結果を出しておき、静置状態で60秒経ったらExcelなどに全部コピペして保存します。

保存した結果は以下のようなデータになります。100Hzサンプリングなので1分だと6000サンプルありますが、全部貼るとあまりにも長いので以下では10サンプルだけ示しています。

-0.49799 0.03831 9.15543 -0.02942 0.47072 8.04147 -0.08255 0.07657 9.77229 -0.27949 0.12148 9.79033 -0.03677 0.15748 9.72942
-0.49799 -0.03831 9.27035 -0.04903 0.56879 8.0905 -0.08136 0.07478 9.77229 -0.27874 0.11329 9.78735 -0.03865 0.15499 9.73061
-0.5363 0 9.23204 -0.02942 0.5884 8.06108 -0.08016 0.07478 9.75973 -0.27055 0.11925 9.78437 -0.03873 0.154 9.7326
-0.45969 0 9.15543 0 0.51975 8.05127 -0.07657 0.07717 9.77289 -0.27651 0.11776 9.79182 -0.03655 0.15388 9.73111
-0.5363 0.03831 9.27035 0.00981 0.53937 8.01204 -0.08016 0.06341 9.77468 -0.27725 0.11999 9.78884 -0.03923 0.15396 9.72942
-0.49799 0 9.23204 0 0.52956 8.01204 -0.07836 0.08195 9.7675 -0.278 0.11925 9.7881 -0.03865 0.15369 9.72877
-0.57461 0 9.23204 -0.00981 0.51975 8.05127 -0.08016 0.07178 9.78425 -0.26906 0.11552 9.79331 -0.03451 0.15652 9.73226
-0.49799 0.03831 9.27035 0.00981 0.50995 8.07089 -0.08674 0.08136 9.77528 -0.27576 0.11627 9.78437 -0.03609 0.15468 9.73222
-0.5363 0 9.23204 0 0.51975 8.04147 -0.08136 0.07478 9.77348 -0.27502 0.11925 9.7866 -0.03831 0.15457 9.7303
-0.45969 0.03831 9.15543 0.05884 0.52956 8.08069 -0.08495 0.0664 9.77707 -0.26756 0.11403 9.78586 -0.03831 0.15453 9.72977

6000サンプルを標本として標準偏差を出した結果は次の図のようになりました。

3成分を平均した標本標準偏差は以下の表の通りになりました。

	ADXL345	LIS3DH	LSM6DSO	IIS3DHHC	ADXL355
標本標準偏差 [m/s^2]	0.03653	0.02793	0.00531	0.00356	0.00124

図と表どちらからも、ADXL345とLIS3DHの2つはほか3つに比べてかなりノイズが大きい、すなわち弱い揺れがノイズに埋もれやすいことが分かります。

計測震度の算出

先ほど出したセンサ毎の差は計測震度にどう影響するのでしょうか。ここでは、気象庁HPの計測震度の算出方法に基づいて、計測震度を比較します。

先ほどの60秒分（6000サンプル分）の加速度データを aac.dat に保存してある前提で、以下のようなプログラムを実行します。

import numpy as np

def main():
    # 全センサの加速度データ
    waveforms = np.loadtxt("acc.dat")

    # m/s^2からgal(cm/^2)に変換するため100倍
    waveforms = waveforms * 100

    # 成分ごとに取り扱えるように配列の軸を入れ替え
    waveforms = waveforms.T

    # 成分ごとにオフセット除去
    waveforms = waveforms - np.tile(waveforms.mean(axis=1), (waveforms.shape[1], 1)).T

    # 成分ごとに周波数補正フィルタ適用
    filtered_waveforms = np.empty_like(waveforms)
    for waveform_idx, waveform in enumerate(waveforms):
        # フーリエ変換
        sp = np.fft.fft(waveform)

        # 周波数補正フィルタ適用
        sample_num = len(waveform) # サンプル数
        sample_spc = 0.01 # サンプリング周期（0.01秒=100Hz）
        freqs = np.fft.fftfreq(sample_num, sample_spc)
        for i in range(sample_num):
            if freqs[i] == 0:
                # 直流成分は除去
                sp[i] = 0
            elif i < np.ceil(sample_num / 2):
                # ナイキスト周波数までは通常通りに適用
                sp[i] = sp[i] * get_filter_coef(abs(freqs[i]))
            else:
                # ナイキスト周波数よりも高い周波数では折り返して適用
                sp[i] = sp[i] * get_filter_coef(abs(freqs[sample_num - i]))

        # 逆フーリエ変換
        filtered_waveform = np.fft.ifft(sp).real
        filtered_waveforms[waveform_idx] = filtered_waveform

    # センサごとに計測震度を出力
    sensor_names = ["ADXL345", "LIS3DH", "LSM6DSO", "IIS3DHHC", "ADXL355"]
    sensor_waveforms = filtered_waveforms.reshape(filtered_waveforms.shape[0] // 3, 3, filtered_waveforms.shape[1])
    for sensor_name, sensor_waveform in zip(sensor_names, sensor_waveforms):
        # 加速度を3成分合成
        resultant_accelerations = np.sqrt(sensor_waveform[0] ** 2 + sensor_waveform[1] ** 2 + sensor_waveform[2] ** 2)

        # 合成加速度の絶対値がある値a以上となる時間の合計を計算したとき、これがちょうど0.3秒となるようなa
        # 100Hzサンプリングの場合は降順ソートで30番目の値
        a = sorted(resultant_accelerations, reverse=True)[29]

        # 計測震度を算出
        s = 2 * np.log10(a) + 0.94

        # 出力
        print(f"{sensor_name:8}: {s:8.5f}")


# フィルター特性
def get_filter_coef(freq: float) -> float:
    y = freq / 10

    lowpass_filt  = (1 - np.exp(-(freq / 0.5) ** 3)) ** (1 / 2)

    highpass_filt = (1 + 0.694 * y**2 + 0.241 * y**4 + 0.0557 * y**6 +
                     0.009664 * y**8 + 0.00134 * y**10 + 0.000155 * y**12) ** (-1 / 2)

    period_filt   = (1 / freq) ** (1 / 2)

    return lowpass_filt * highpass_filt * period_filt


if __name__ == "__main__":
    main()

実行結果は以下のようになりました⁶。

ADXL345 :  1.93675
LIS3DH  :  2.17672
LSM6DSO :  0.35153
IIS3DHHC: -0.01853
ADXL355 : -0.91114

計測震度は小数点以下第3位を四捨五入し、第2位を切り捨てます。その結果を標本標準偏差と併せて表にすると以下のようになりました。

	ADXL345	LIS3DH	LSM6DSO	IIS3DHHC	ADXL355
標本標準偏差 [m/s^2]	0.03653	0.02793	0.00531	0.00356	0.00124
計測震度	1.9	2.1	0.3	0.0	-0.9

所感

ADXL345は、静置状態で計測震度1.9程度となりました。つまり、震度2程度以下の揺れはノイズに埋もれて検知できないということになります。
データシートを振り返ると分解能は3.9mg=4gal程度であり、0galから4galまでの間は技術的に検知できないのです。震度3以上の揺れを検知したい！という目的では価格面からもかなり優秀ですが、震度1程度の小さな揺れも知りたい、という目的には適しません。

LIS3DHは、静置状態で計測震度2.1程度となりました。ADXL345と同様に、震度2程度以下の揺れはノイズに埋もれて検知できないということになります。
データシート上はADXL345よりも分解能もノイズ密度も良いのですが、計測震度がこちらのほうが高くなってしまっています。考えられる原因として、LIS3DHは計測震度に与える影響が大きい周期帯（0.1秒-10秒）のノイズが比較的大きいのかもしれません。

LSM6DSOは、静置状態で計測震度0.3程度となりました。
有感（震度1以上）の基準は計測震度0.5以上ですから、人間が「あれ地震かも」と思ったときにきちんと波形に収まるくらいのセンサがほしい、という目的にはぴったりです。ただ、ノイズの乗り方によっては時々計測震度0.5を超えるときがある可能性もあり、ノイズを抑える工夫は別途必要かもしれません。
国内ではあまり使いやすい形での流通はないのですが、アリエクとかDigiKeyでは1個1000円弱～1500円くらいで取り扱いがあるので、コスパが良いです。

IIS3DHHCは、静置状態で計測震度0.0程度となりました。
こちらも静置状態では有感の基準を下回っています。計測震度の計算式に基づくと、計測震度0.5に到達するにはノイズが約2倍まで膨らむ必要がありますから、静置状態で時々有感判定になるリスクは低いと言えると思います。
ストロベリーリナックスで取り扱いがあるため手軽に入手できるのも心強いです。

ADXL355は、静置状態で計測震度-0.9程度となりました。
筆者はこのセンサを含め5センサを定常的に動かしていますが、緊急地震速報（予報）の受信時にADXL355の波形を見ていると、体感では揺れがないのにセンサに地震っぽい揺れが記録されていることが度々あります。つまり、体感で揺れたか揺れてないかくらいの揺れはこのセンサはもうお手の物で、さらに体感できない揺れもこのセンサでキャッチできる、ということです。
正直、1万円以下でここまで性能の良い地震観測ができるのはかなりすごいと思います。STマイクロさんありがとうございます。

おわりに

震度3以上の揺れを押さえたい → ADXL345、LIS3DH が便利
震度1以上の揺れを押さえたい → LSM6DSO、IIS3DHHC が便利
体に感じない揺れも押さえたい → ADXL355 が便利

といった感じになりました。ちゃんと比べてみると目的別に色々見えていて面白いですね。

宣伝

今回は色々なプログラムを作りましたが、少し手間もあるものです。
そんな、おうちで地震観測を簡単に叶えられる品 EQIS-1 が、今年8月に発売されました。こちらはLSM6DSOを使用しているので、震度0と震度1の判別も出来ちゃいます。
ささやかながら私も企画に参加させていただいております。ぜひ。

高分解能モード時。標準モードでは10bit、省電力モードでは8bitです。 ↩
ジャイロ（傾斜）センサ付きの6軸センサです。今回は加速度センサだけ使います。 ↩
Qwiicシステム用に販売されているので、素で使うのにはあまり適していないかもしれません。 ↩
SCKではなく、SCL、SCKL、SCLKと表記されていることがあります。 ↩
https://www.data.jma.go.jp/eqev/data/kyoshin/jishin/collection.html より ↩
計測震度は加速度の常用対数の2倍に比例します。つまり、性質として加速度が10倍になれば計測震度は2だけ上がり、同様に加速度が31.6倍になれば3だけ上がります。ADXL345とADXL355を比べると、加速度の標準偏差はちょうど30倍近く違いますので、計測震度が3ずれるというPythonプログラムの計算は合っている…はずです。間違ってたらごめんなさい。 ↩