Raspberry Pi PicoとA/D変換ボードで作る地震計

はじめに

アドベントカレンダーに参加したことが無く、かつなかなか合うテーマが無くて尻込みしていました。
今回、なんとなく投稿できそうなカレンダーに空きがあったので初投稿です。

今回はこれらの投稿のまとめです。

ハードウェア

メインのコントローラには　RaspberryPiPicoを使いました。今回のプロダクトはI2Cが使えれば良いので、主な理由は「使ってみたかったから」になるかもしれません。

(これまでRaspberry Pi や Pi Zero系以外には Arduino系や ESP** を使ってIoT的な遊びをやってきました。Picoは買ったもののLチカで満足して出番を待っている状態でした)

信号の増幅とデジタイズにはSX8725Cを採用して、ほぼ満足する結果を得ました。

SX8725C
分解能	16 bit
アンプ	3段(最大1000倍)
動作電圧	2.5~5.5 V
内蔵基準電圧	1.22 V
変換レート	0.483 ksps
I2C転送速度	最大400kHz
消費電流(5V時)	250 uA ~

配線はとてもシンプルで、PicoとADCは電源, GND, SDA, SDCの計4本、センサ側は差動入力の2線でOK。今回はセンサ部分にジオフォンを繋いでシステムを「地震計」として使いますが、出力が微小電圧のセンサなら汎用的に使える電圧測定器になっています。今回の設定では測定電圧範囲が±0.61[V]で、分解能は16bitです。

ソフトウェア

ハードウェア側

慣れ親しんだPythonでの開発ができるところはPicoの利点と思っていて、今後の利用頻度が上がりそうです。今回初めてMicroPythonを使ってみましたが、特に違和感なく開発できました。

Raspberry Pi Pico の main.py はこちら↓↓

ソースコード（長いので折りたたみ表示）

SX8725C.py

import utime
import machine

led = machine.Pin(25, machine.Pin.OUT)
RT6150_PS = machine.Pin(23, machine.Pin.OUT)
timer = machine.Timer()

led.value(1)
utime.sleep(0.5)
led.value(0)
utime.sleep(0.5)

i2c_bus_number  = 0
i2c_sda=machine.Pin(0)
i2c_scl=machine.Pin(1)
i2c=machine.I2C(i2c_bus_number,sda=i2c_sda, scl=i2c_scl, freq=100000)

#if we use only one i2c device, this code is useful to find the address.
addrlist = i2c.scan()
print( "address is :" + str(hex(addrlist[0])) )

addr = addrlist[0]
#addr = 0x48  # Default address of SX8725C
i2c_address = addr

Vref = 1.22  # for internal refernce
#Vref = 5.00   # for Vbatt reference
conversion_factor = Vref / (65535)
zero_shift = -3.7
#zero_shift = 0


# Define registers
RegACOutLsb = 0x50
RegACOutMsb = 0x51
RegACCfg0 = 0x52
RegACCfg1 = 0x53
RegACCfg2 = 0x54
RegACCfg3 = 0x55
RegACCfg4 = 0x56
RegACCfg5 = 0x57
RegMode   = 0x70

# Define values
ValACCfg0 = b'\x1E' #30   #0b00011110
ValACCfg1 = b'\xFF' #255  #0b11111111
ValACCfg2 = b'\x00' #0    #0b00000000 #PGA2x1
ValACCfg3 = b'\x8C' #140  #0b10001100 #PGA1x10 PGA3 12/12
ValACCfg4 = b'\x00' #0    #0b00000000 #PGA3OFF 0 
#ValACCfg5 = b'\x02' #2    #0b00000010 #diff, VBATT
ValACCfg5 = b'\x03' #3    #0b00000011 #diff, Vref
ValMode   = b'\x88' #132  #0b10000100 #pump off

reglist = []
reglist = list(range(0x50, 0x50+8))
reglist.append(0x70)

# --- Check default values (not required) ---
data=[]
for reg in reglist:
    readdata = i2c.readfrom_mem(addr,reg,1)
    data.append(int.from_bytes(readdata,'big'))
print(data)
#--------------------------------------------

#init
i2c.writeto_mem(addr,RegACCfg0,ValACCfg0)
i2c.writeto_mem(addr,RegACCfg1,ValACCfg1)
i2c.writeto_mem(addr,RegACCfg2,ValACCfg2)
i2c.writeto_mem(addr,RegACCfg3,ValACCfg3)
i2c.writeto_mem(addr,RegACCfg4,ValACCfg4)
i2c.writeto_mem(addr,RegACCfg5,ValACCfg5)
i2c.writeto_mem(addr,RegMode,ValMode) 

# --- Check current values (not required) ---
data=[]
for reg in reglist:
    readdata = i2c.readfrom_mem(addr,reg,1)
    data.append(int.from_bytes(readdata,'big'))
print(data)
#--------------------------------------------


def sign16(x):
    return ( -(x & 32768) | (x & 32767) );

#main
def readdata(timer):
    data = []
    for reg in (RegACOutLsb, RegACOutMsb):
        readdata = i2c.readfrom_mem(addr,reg,1)
        data.append(int.from_bytes(readdata,'big'))
    data16 = (data[1]<<8 | data[0])

    dataDiff = sign16(int(hex(data16),16))
    mVol = 1000 * dataDiff * conversion_factor
    var = "{:0=+12.4f}".format(mVol + zero_shift)
    print (var)

led.value(1)
RT6150_PS.value(1)
timer.init(freq=100, mode=machine.Timer.PERIODIC, callback=readdata)

PC側

データの後処理を考慮して、3本目の記事公開後に収録用のコードを少し改変しました。具体的には、きちんと秒の変わり目(00秒)から100Hzで記録できるように割り込み処理でシリアルのデータを回収するようにしました。

シリアル入力にタイムスタンプを付与して保存するコード

ソースコード（折りたたみ表示）

serialread.py

import time,signal
import sys
from datetime import datetime
import serial

dt_now=datetime.now()
fnb = dt_now.strftime('%Y%m%d%H%M%S')
fn = fnb + ".csv"

ser = serial.Serial("/dev/cu.usbmodem14401", 9600)
f = open(fn, 'a')
line = ser.readline().decode('utf-8').rstrip()

# インターバルタイマーハンドラー
def intervalHandler(signum,frame):
    dtnow = datetime.now()
    char=str(dtnow)
    out = char[0:22] + "+09:00, " + line + "\n"
    f.write(out)

# インターバル起床するハンドラを設定(SIGALRMハンドラー)
signal.signal(signal.SIGALRM,intervalHandler)

# 秒の変わり目をポーリング
while True:
    dt_now = datetime.now()
    if dt_now.microsecond < 1000: # 秒の変わり目から1ms未満ならOK
        print('microsec = %d' % dt_now.microsecond)
        break

# インターバルタイマー設定
signal.setitimer(signal.ITIMER_REAL, 1, 0.01)

try:
    while True:
       # インターバルタイマー待ちの間はシリアルの読み込み
       # ここで拾っておかないと、ハンドラのタイミング次第で間に合わない事がある
       # 場合によっては取りこぼすかもしれないが、処理落ちするより良い事にした。
        line = ser.readline().decode('utf-8').rstrip()

except KeyboardInterrupt:
    print('\nCTRL-C pressed!!')
    f.close()
    sys.exit()

得られるのは タイムスタンプ, 電圧値, ADカウント の生値がひたすら続くCSVです。地震の観測データではデータの劣化を防ぐためカウント値を使うのが一般的で、最後に物理量（加速度、速度や変位など）に変換します。

YYmmddHHMMSS.csv

...
2022-12-04 16:13:32.40+09:00, -000000.9916,+000124
2022-12-04 16:13:32.41+09:00, +000001.8566,+000277
2022-12-04 16:13:32.42+09:00, -000001.4942,+000097
2022-12-04 16:13:32.43+09:00, -000001.4942,+000097
2022-12-04 16:13:32.44+09:00, +000000.2370,+000190
2022-12-04 16:13:32.45+09:00, -000000.5262,+000149
...

データ処理

地震計として使うために、CSVからmseedやsac形式に変換します。ついでにデータの概要が分かるようにプロットも作成してみました。

このコードはまだスクラッチ状態ですが、最低限の機能は果たします。

ソースコード（折りたたみ表示）

obspytest.py

#memo obs sgy PasteUp

from obspy import read
from obspy import Stream, Trace
from obspy.core import UTCDateTime
import pandas as pd
import os
import matplotlib.pyplot as plt
import datetime as dt
import sys

outputdir='./'

args = sys.argv
inputfile = args[1]
outputdir = args[2]

dataFileLoc = "./" #location of the ascii file

basename = os.path.splitext(os.path.basename(inputfile))[0]

df1 = pd.read_csv(inputfile, header=None, names=('datetime','mV','Z'),parse_dates=['datetime'])
df1 = df1.set_index('datetime', drop=True)

stnName = "HOM"

print(df1)

statsZ = {}
statsZ['network'] = "HOME" #fake network name
statsZ['station'] = stnName
statsZ['component'] = 'Z'

## Station location
statsZ['stla'] = 34
statsZ['stlo'] = 135

## Event location
statsZ['evla'] = 34
statsZ['evlo'] = 135
statsZ['evdp'] = 0

dt = UTCDateTime(df1.index[0])
sdt = dt

statsZ['starttime'] = sdt
statsZ['sampling_rate'] = 100.0
trZ = Trace(data=df1['Z'].values, header = statsZ)

st = Stream(traces=[trZ])

print(st)

## Write the stream to file
st.write(os.path.join(dataFileLoc, f"{basename}-Z-obspy.mseed"), format="MSEED")
st.write(os.path.join(dataFileLoc, f"{basename}-Z-obspy.sac"), format="SAC")

## Plot the stream
stFig = st.plot(show=False,
        size=(1500,600), number_of_ticks=6,
        type='relative', tick_rotation=60, handle=True,
        linewidth = 0.3, starttime=sdt)
plt.savefig(os.path.join(outputdir, f"{basename}-Z-obspy.png"), dpi=300)

trZ.plot(type='dayplot', interval=60, color=['k'], vertical_scaling_range=1e5, outfile=(os.path.join(outputdir, f"{basename}_Z.png")))

st_et = UTCDateTime(df1.index[-1])
st_st = st_et - 60
st_r = st.trim(st_st, st_et)
print(st_st)
strFig = st_r.plot(show=False,
        size=(1500,600), number_of_ticks=6,
        type='relative', tick_rotation=60, handle=True,
        linewidth = 0.3, starttime=st_st)
plt.savefig(os.path.join(outputdir, f"{basename}-Z-recent.png"), dpi=300)
st_r.spectrogram(log=True, per_lap=0.8, outfile=(os.path.join(outputdir, f"{basename}_Z_spec.png")))

出力例：

課題

ジオフォンを繋ぐとどうしても低周波ノイズが現れてしまい、周波数が1Hz前後と観測ターゲットの帯域に含まれているのが課題です。ジオフォンと言っても中身はただのコイルなので、共振してしまうのか？とも思ったのですが、アンプの1/fノイズなど様々な原因を切り分けることが出来ていません。

センサーを接続して周囲の振動ノイズが少ない時に見てみると、看過できない綺麗な1Hzの波が見えます。高調波の2Hzに加えてセンサーの固有周期と思われる15Hz程度のピーク、60Hzの折り返しである40Hz(および高調波の20Hz)などが混ざっているようです。都会など振動ノイズが大きい場所ではノイズの10%程度の振幅なので埋もれてしまいますが、微小地震観測を目的とした場合には失格ですね。

センサーを外して、ダンピング用の負荷抵抗のみの時は、殆どホワイトノイズで良さそうに見えます。しかしよく見ると0.9Hzくらいに他より卓越したピークが認められます。センサー接続時に増長しなければ問題無い振幅なのですが、問題は回路側にあるようです。

まとめ

安価な材料で微小電圧測定器を作成してみました。結果はやはり厳しいものでしたが、趣味として楽しんだり実演用のガジェットとしては十分使用できるものができました。

研究用途で使う低ノイズ・低消費電力の本格的なロガーは軽く数百万円ですので、仕方ありませんね。もう少し安価でホビーと実用の間くらいの地震計としては、RaspberryShakeという製品があるので、こちらもいつか試してみたいです。センサー込みで424ドルくらい（から）なので個人でも何とか手が出るかな、、という感じです。

地震計ネットワークに加入するためには設置場所が良くないといけないので、持ち家など地面（あるいは地下室）に直接置ける環境じゃないと厳しいかもしれません。