背景

去年のアドベントカレンダーで紹介した地震計は1Hz前後の低周波ノイズを除去しきれなかった。

今年は同じテーマで別のAD変換ボードを使ってみることにした。蓑虫クリップ以外のハンダ付けが無かったので、自作と言えるかは微妙。

ハードウェア

使用したのは千石で手に入る「High-Precision AD HAT」。量子化ビット数が32bitと高分解能で、買ってすぐに使えるのが決め手。

マイコンはAD変換ボードがぴったりサイズの Raspberry Pi Zero W を選択（当初 Raspberry Pi 4 でテストしたが、ノイズが大きく断念した。手待ちの電源に制約があり、本体と電源どちらの問題か切り分けができなかった）。

取得したデータをmseed形式で保存したいと思い、計測から1次処理まで1台で済ませるにも丁度良い組み合わせとなった。

GPIOにHATを接続するだけでハードウェアの準備は完了。ジオフォンをIN1,IN2番端子に繋いで観測する。今回の用途では差動入力で微小電圧の測定を行う。これらはソフトウェア側で設定できる。

ソフトウェア

OSは Raspberry Pi OS Lite (32bit; Kernel Ver. 6.1)。Raspberry Pi はオーバーレイファイルシステムで動かし、データはNASをマウントして保存することにした。

ADS1263の制御には発売元のWikiに従いgithubから入手できるサンプルプログラムの完成度が高く、ほとんどそのまま使える。

AD値の読み出しからmseed形式での保存までをPythonで用意した。サンプリングは100Hzで行い、タイムスタンプを付けるためにインターバルタイムハンドラを使用した。

使用した環境では "gpiomem-bcm2835" が存在せず、ADS1263.py でimportされる公式の config.py ではRaspberryPiZeroがJestonと認識されてしまう。そこで、124行目を下記のように変更した。

config.py

# 変更前：
if os.path.exists('/sys/bus/platform/drivers/gpiomem-bcm2835'):

# 変更後：
if os.path.exists('/sys/bus/platform/drivers/gpiomem-bcm2835') or os.path.exists('/sys/bus/platform/drivers/rpi-gpiomem'):

ソースコード全体はこちら

seismo2.py

#!/usr/bin/python
# -*- coding:utf-8 -*-

import time
import ADS1263
import RPi.GPIO as GPIO

import signal
import sys
import os

from datetime import datetime, timedelta, timezone

from obspy import Stream, Trace
from obspy.core import UTCDateTime
import numpy as np


REF = 5.08          # Modify according to actual voltage
                    # external AVDD and AVSS(Default), or internal 2.5V
ADC = ADS1263.ADS1263()

JST = timezone(timedelta(hours=+9), 'JST')

dataFileLoc = "/mnt/nas/obs_geophone/raw" #location of the mseed file

stnName = "HOM"

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

statsZ['sampling_rate'] = 100.0

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

rawdata = 0
ns = 0
tr = []
fnb="0000000000"

# mseed file length by min
#minutes = 1
minutes = 30
#minutes = 60
#minutes = 180

sdlen = minutes * 60 * statsZ['sampling_rate']

channelList = [0]  # The channel must be less than 10

def writedata():
    global ns
    global statsZ
    global fnb
    global tr

    fnb = dt_now.strftime('%Y%m%d%H%M%S')
    statsZ['starttime'] = UTCDateTime(dt_now)
    tra = np.array(tr)
    tr = []
    print(fnb,tra,ns)

    ns = 0    
    trZ = Trace(data=tra, header = statsZ)
    st = Stream(traces=[trZ])
    st.write(os.path.join(dataFileLoc, f"{fnb}-Z.mseed"), format="MSEED")

    return 0

# インターバルタイマーハンドラー
def intervalHandler(signum,frame):
    global ns
    global tr
    global rawdata
    global dt_now

    ns = ns + 1
    if ns == 1: dt_now = datetime.now(JST)
    val = rawdata
    tr.append(val)

    if ns == sdlen:
        dummy = writedata()

def getVol(raw):
    if(raw>>31 ==1):
        data = -1*(REF*2 - raw * REF / 0x80000000)
    else:
        data = (raw * REF / 0x7fffffff)
   
    return int(data*100000000)


def get30Int(raw):
    if(raw>>31 ==1):
        data = -1 * (0xffffffff - raw)
    else:
        data = raw
   
    return int(data/4)  # STEIM2 can treat only 30bit...


def readData():
    global rawdata
    
    # The faster the rate, the worse the stability
    # and the need to choose a suitable digital filter(REG_MODE1)

#    if (ADC.ADS1263_init_ADC1('ADS1263_400SPS') == -1):
    if (ADC.ADS1263_init_ADC1('ADS1263_100SPS') == -1):
        exit()

    ADC.ADS1263_SetMode(1) # 0 is singleChannel, 1 is diffChannel
   
    while(1):
        ADC_Value = ADC.ADS1263_GetAll(channelList)    # get ADC1 value

#        rawdata = getVol(ADC_Value[0])
        rawdata = get30Int(ADC_Value[0])

    ADC.ADS1263_Exit()

    return rawdata

#=====================================================================

# インターバル起床するハンドラを設定(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)
        fnb = dt_now.strftime('%Y%m%d%H%M%S')
        break

# インターバルタイマー設定
# 最初の開始タイミング1秒後、以後0.01秒間隔
signal.setitimer(signal.ITIMER_REAL, 1, 0.01)


try:
    while True:
        rawdata = readData()

except IOError as e:
    print(e)
   
except KeyboardInterrupt:
    print("ctrl + c:")
    print("Program end")
    ADC.ADS1263_Exit()
    dummy = writedata()
    sys.exit()