試作品を公開しました。
https://github.com/T40O0/ADXL355_SPI_M5_SD_FIR.git
振り返ると、1年かかっています（ラベル：MEMS Seismometer）。仕事が終わってからコツコツ積み上げて、ようやく形にしました。
MEMSセンサーの取り扱いができるようになったという点よりも、FIRフィルタの働きを理解できた点のほうが勉強になりました。仕事では完成した地震計を扱うので、周波数領域で正しくデータを扱えるようにするためのハードウェア側の知識が不足していました。やはり、手は動かすものです。
来年は何にチャレンジしましょうか。
a.txtをb.txtにリネーム
$ mv a.txt b.txt
*.txtをdir1の下位に移動
$ mv *.txt dir1/
dir1をdir0の下位に移動
$ mv dir1/ dir0/
-vで経過表示
$ mv −v dir1/ dir0/
ftp サイトからカレントディレクトリにコピー
$ wget ftp://ID:PASS@address/a/b/c/*
DAS の noise reduction が目に留まりました。
2016-06-RECORDER-Simultaneous_Acquisition.pdf (mit.edu)
Panel B shows the result of applying the noise-reduction processing as a running weighted average over 11 adjacent channels (20m).
このような方法で良いのか？奇麗な結果だから良いか。
研究内容を社内アドバイザーに説明し、承認をいただく必要があります。
が、アドバイザーとは名ばかりの executives が多く、逆に技術的な内容を嚙み砕いた資料を作成しなければなりません。
そこで GPT-4 の出番です。聞かれそうな内容を代わりに聞きます。
「○○について技術的内容を教えて」
「○○の海外での適用事例を紹介して」
「○○の国内での適用事例を紹介して」
これを Word データに書き出し、絵を入れて体裁を整えたら７割完了です。味も深みもありませんが、ソツのない内容の資料ができました。この程度なら、御自身で聞いて理解していただくと、少しは有意義な議論ができると思うのですが。
社内向けの資料にGPTの回答レベルがちょうど良いのは皮肉ですが、無駄な作業に対する時間短縮と割り切り、楽になったと考えましょう。
Raspberry Pi 4 を定期的に再起動する必要が出てきました。
調べてみるといくつかありますね。
少し古そうですが、簡単な以下で書き込みました。
introduce tasks to be run by cron.
$ sudo crontab -e
minute (m)
hour (h)
day of month (dom)
month (mon)
day of week (dow)
'*' for 'any'
m h dom mon dow command
5 0 * * * /sbin/reboot
ESurf - Seismic monitoring of torrential and fluvial processes (copernicus.org)
地震計を使って何を測っている？に対して総論のようにまとめられている文献です。
河川
Gimbert et al.（2014）
乱流によって誘発されるノイズは、床荷重によって誘発されるノイズよりも低い周波数で励起される。
発生源と受信源の距離にも敏感で、河川に近い（遠い）観測点ほど、底質（乱流）によって誘発されるノイズの割合が大きい。
観測されたヒステリシスは、8Hz以下の乱流の影響を大きく受けている。
図7
アルプスの「サンピエールの激流」放水量の比較。2-5Hzの周波数帯域は流出変動を最もよく説明し、20-30Hzの周波数帯域はこの関係を失っている。
落石、岩すべり、岩なだれ
地震と同様に、数十秒から数分という継続時間を持つ。
多くの落石における地震エネルギーの長い立ち上がり時間は、地震における急激な上昇とは対照的であり、識別に利用することができる。
岩石崩壊
Dammeier et al.（2011）
スイス地震学研究所で記録された岩石崩落の地震信号を用いて、
など、崩壊と地震の特徴との関係を定義。
地すべり
地震信号に含まれる超低周波（<0.1 Hz）は、あまり減衰せずに長距離を伝播する。
ほとんどが高周波信号（>1Hz）を生成する落石や岩なだれのような他の斜面事象との大きな違いである。
土石流
広い高周波数帯域（1-60Hz）で地震エネルギーの急激な増加を示す（図8）。
風
局所的な気圧場の変動により、30秒を超える周期の振動を発生させることがある（例えば、De Angelis and Bodin, 2012）。
雨
Burtin et al. (2011)
70Hz以上の周波数の雨信号を確認したが、これは大きな岩の隣に設置されたステーションにのみに存在。
大きな岩に近接していない場所に設置されたステーションでは、この信号は見られなかった。
柔らかい地面や土の上に落ちた雨滴は、効果的な信号源にはならないかもしれない。
雷、稲妻
20Hz以下の低周波信号を発生（Scarpetta et al., 2005; Assink et al., 2008 ）
音響領域と地震領域の結合が効果的であれば、地震計で記録可（Walker et al. , 2011）
人為的な騒音
明確な日周期が生じる。
多くの場合持続時間が短く、主に1-20Hzの周波数帯域に影響を与える（図8）。
振幅は通常比較的小さく、近傍の観測点でのみで記録される。
人為的な信号はアレイ全体にわたってコヒーレントであることはまれであり、分離に役立つ。
火山性地震と落石の Automated identification に利用する特徴量の候補５つ。
SDカードにデータを記録していましたが、トラブル発生。
なぜか一定期間を過ぎると記録が停止していました。
調べてみると、ルートのファイル数が21,844。どの機種も記録日時は異なるのに同じ数。この数字、記憶にありました。
FAT32の制限です。
ファイルシステムに応じた総数とディレクトリ毎のファイル数に制限がありました。
https://qiita.com/ryoma-jp/items/f9124a057ac8c646ae61
対策は簡単で、異なるファイルシステムでフォーマットする、ある程度の期間でファイルをまとめる、1日毎にフォルダを分ける、ファイル名を短くするなど。
あまり平たい階層にしても整理時のアクセスに問題が生じるので、フォルダ分けの方針で進めましょう。
https://www.t3.gsic.titech.ac.jp/node/441
先の論文に出てきたゼロクロス法の説明は以下がわかりやすいと思います。
Determination of surface‐wave phase velocities across USArray from noise and Aki's spectral formulation - Ekström - 2009 - Geophysical Research Letters - Wiley Online Library
Because the amplitude of the real part of the spectrum depends on both the background noise spectrum and non-linear effects of the data processing, dispersion information cannot readily be deciphered from the detailed shape of the spectrum. The locations of the zero crossings in the spectrum should, however, be insensitive to variations in the spectral power of the background noise, and we choose to use the locations of these zero crossings as the dispersion observables.
前提として、ランダムに全周囲から表面波が到来する場合、２点間の相互相関係数は第一種0次のベッセル関数J0を介して表現できる（それは円周上の１点と中心点とのそれを円周上の多点で平均した空間自己相関係数（SPAC係数）と一致する）。
J0（ωr/c(ω)）=0 の時の角周波数ωがゼロクロス点でわかり、２点間の距離rも既知のため、位相速度 c(ω) が求まる流れです。０の時を選ぶのは上記理由です。
以前、passive 手法等をまとめた際に含めていたかな？と思い見返しましたが、入れていませんでした。
https://phreeqc.blogspot.com/2021/05/passive.html
クロスコリレーションのゼロクロッシング位置を用いて観測曲間の分散曲線（位相速度）を取得しています。ゼロクロス法については後日。
文献の対象は米国でしたが、日本でも観測局の長期データを基に速度構造を推定している例をいくつか見たことがあります。ま、日本の場合はもっと細かいデータが公開されていますので、積極的な利用はないでしょう。
物理探査にも使えそうです。
フローレンスとローマの間に、チヴィタという町があります。溶結凝灰岩の崖の上に孤立する町で、世界遺産候補です。観光客からは「死にゆく町」として人気だそうです。
今回の出張で見たかった場所の一つだったのですが、見学が中止になりました。が、それでも個別に見に行かれた地質屋さんがいらっしゃいました。後日、写真を見せてもらうと、香川県の屋島の上に中世の村が作られているような印象でした。
中には、鋼管杭？からアンカーをとって崖を引き留めているような写真もありました。このような対策は日本にもあるそうです。以下の文献には、おおまかな構造が示されています。
(PDF) The dying town of Civita di Bagnoregio and the killer landslide (researchgate.net)
当然ながら、古くから「dying」だった訳で、似たような姿は絵画の題材にもされています。多くはデフォルメされて描かれているようです。いずれにしても、古くから人を惹きつける景観であったことは間違いありません。
対策は観光地としての延命治療のようにも感じますが、住民の方にとっては確かに死活問題です。厳しい土地で工夫しながら長期にわたり歩んでこられた歴史。対策も含め世界遺産候補なのでしょう。
機会があれば、じっくり見てみたいですね。
reactive transport を扱えるOpenFOAM11 用のソルバーが実装されています。二重空隙モデルを扱える点が特徴でしょうか。反応はPhreeqcRM 利用とのこと。
H/Vをすべり面深度決定に利用している方がいらっしゃいました。
まずは、H/Vの基本。
中村 (2008) H/V スペクトル比の基本構造
H/V の地すべり調査への適用例です。
F. Panzera et al. (2012) Seismic site response of unstable steep slope using noise measurements: the case study of Xemxija Bay area, Malta
V. Pazzi et al. (2016) H/V measurements as an effective tool for the reliable detection of landslide slip surfaces: Case studies of Castagnola (La Spezia, Italy) and Roccalbegna (Grosseto, Italy)
つまり、岩盤が崩壊したような基盤とすべり土塊のインピーダンスに明確なコントラストがある場合に H/V を地すべり調査に利用できるということです。地すべり土塊と岩盤の間に速度変化がない場合や、深度とともに速度が徐々に増加する場合、この手法を適用することはできません。これは、地質屋さんが現場を見たらある程度判断できますし、ボーリングを掘れば明らかになります。
深度を決定するには Vs が必要です。単点に加えアレーを加えたら良いのでしょう。が、大雑把な手法なので、やりすぎると矛盾が出るでしょうね。ボーリングで確認した深度に一致するよう Vs を調整し、それが妥当な場合に「適用可」と判断するほうが良いのでしょう。
いずれにしても非常に手軽な補間手法です。積極的に利用したいですね。
桜島 19日の噴火。
フェリー乗り場から迫力のある噴煙が見えました。
これでレベル４にならないのか地元の年配の方に伺うと、これは日常とのことでした。笑ってらっしゃったので、何ともないのでしょう。いや、すごい。
さらに話を伺うと、昭和30年ごろまでは北側から御岳の広場まで遠足で登っていらしたそうです。噴火による死傷者が出てから立入り禁止になったそうですが、確かに、その頃の地図には山頂へのルートが書かれています。
今昔マップ：https://ktgis.net/kjmapw/kjmapw.html?lat=31.584930&lng=130.652478&zoom=15&dataset=kagoshima&age=2&screen=2&scr1tile=k_cj4&scr2tile=k_cj4&scr3tile=k_cj4&scr4tile=k_cj4&mapOpacity=10&overGSItile=no&altitudeOpacity=2
登山道があったなら、火口の写真が残っているだろうと探しましたが、当時のものは見つかりませんでした（ヤマレコには201〇年に顔出し写真を投稿されている方はいらっしゃいましたが）。
UAVを飛ばすのはどうなのでしょうね。火口を間近で見てみたいものです。
18～20日まで、鹿児島で火山学会の秋季大会がありました。
私は参加していませんが、同僚が参加しており発表要旨を見せてもらいました。
と、俄然、気になるモノが。
DAS (Distributed Acoustic Sensing) で火山性地震や土石流を検知した事例です。
これ、最近の流行ですよね。機材が欲しいのですが、高価なので指をくわえて眺めています。
DASは既設の光ファイバーを利用できるので、接続する機材のイニシャルコストに目を瞑れば非常に手軽です。既に全国に受信器が設置済みと言えますので。そのためか、国内でも地震や火山以外に、いくつかの用途に向けて研究されているようです。例えば、河川。
国総研：河川堤防の変状検知システム実験結果nilim.go.jp/lab/fbg/gijyutsukoubo/jikkennkekka/report/10.enaa.pdf
中国地勢：分布型光ファイバセンサによる河川堤防の 観測方法の検討について
https://www.cgr.mlit.go.jp/tyokugi/74/d4/05.pdf
数百ch、大きいと 1000Hz でデータを取られているので、処理はヘビーです。空間方向を時間方向になおして波形を作る前処理も必要になることがあります。が、Pythonでも処理は可能です。
地すべりや物理探査、実務でも使えるのですが、市場が大きくなるまでは皆さん様子見でしょうね。先行投資できる体力のある会社が羨ましい。
DualSPHysics を用いた土砂移動の新たなモデルかなと期待して読み始めましたが、非ニュートン流体実装を用いた2Dの計算のようで途中でやめてしまいました。モデル自体はこちらがわかりやすいでしょうか。210629.pptx (jsce.or.jp)
不安定化に対する補正を加えられていますね。SPHではこれらの補正がないと色々と不具合が出るので仕方ないのですが、多くの補正が加わると本当に正しい計算なのか疑問に思えてくることがあります（いえ、私が作ったコードのことです）。
ま、それでも補正して動くようにしたものを公開される研究者には頭が下がります。
余裕ができたので、以前購入していた GNSS モジュール L76K を引っ張り出してきました。
ボードは眠っていた M5 ATOM S3。以前、マイクをつけて、avatar を動かして遊んでいたボードです。
https://github.com/meganetaaan/m5stack-avatar/blob/master/README_ja.md
いくつかのサンプルコードが公開されていますので、その中からESP32用を利用。Serial2 の PIN を設定してコンパイルするだけでした。
が、なぜか PC側にデータが送られてこない。
数時間、試行錯誤しました。原因は接続するタイミングでした。マイコンを起動した後にPC側のシリアルモニタを開いて接続する手順ならOK。先に開いて待ち受けにしているとダメ。これは初めてのケースでした。
肝心の位置は遠く離れた場所にプロットされていましたし、時刻も数秒遅れて表示。うーん。先が長そう。
移動中に、引き続き MEMS 加速度センサーのノイズを低減する方法を考えていました。
参考書をパラパラ見ていて、目に留まったのがエイリアシング。
ナイキスト周波数よりも高い側に信号が含まれていると、FFTの結果にノイズとなって現れます。これ、大事です。
https://www.analog.com/jp/technical-articles/elusive-tones-aliasing-effects-in-digital-mems-accelerometers-in-condition-monitoring.html
では、どうすれば良いか？
答えは簡単で、事前に LPF を通せば良いわけです。が、カットする高周波を求めるために FFT をかけてしまうとエイリアシングが生じます。
そこで、時刻歴データの段階で LPF をかけてしまうという発想。これがデジタルフィルタ。FIR と IIR が有名ですが、実装上、計算の軽い前者がよく利用されているようです。PolyPhase フィルタも同様です。
FIR フィルタを今まで使ってきたものの、自分で作ったことはありませんでした。今では Python でも簡単に作成できるようです。
帰社後に探してみると、すぐに見つかりました。scipy.signal.firwin で作成できます。https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firwin.html
FIR filter design using the window method.
This function computes the coefficients of a finite impulse response filter. The filter will have linear phase; it will be Type I if numtaps is odd and Type II if numtaps is even.
これを手元の信号にかけると、位相が遅れます。ま、タップ数の半分の時間を補正してやればOKです。直線位相から最小位相にしたい場合は scipy.signal.minimum_phase で良いみたい。https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.minimum_phase.html
Convert a linear-phase FIR filter to minimum phase
MEMSセンサーから500Hzで吐き出したサンプルに対し作成した係数を乗じ、その後に100Hzへ間引き。Tap 数 201 で 2msec 以内に回るか心配でしたが、全く問題なし。しかも効果は抜群で、ノイズを大きく取り除くことができました。成功です。
解決はしませんでしたが、ヒントはこちら。
中川ほか「WIN フォーマットデータを ObsPy で読み込む 新しいモジュール」https://www.eri.u-tokyo.ac.jp/GIHOU/archive/26_031-036.pdf
Windows 環境でしたので、パッチを充てるだけに留まりましたが、ダメでした。ま、既にツールはあるので支障はありません。寝かしましょう。
MEMS加速度センサーのノイズを低減する方法を考えていました。
いくつかの製品を比較してみたところ、同じセンサーを使用していてもプログラムの組み方でノイズが3倍ほど異なります。尊敬すべきプロのこだわりでしょうね。でも、どうするのが正解なのでしょうか？
数日悩みましたが、解決せず。やはりプロの領域です。
素人でもできそうなのはスタック。物理探査ではおなじみの方法です。最後に試してみました。
まずは1台。（ノイズに無関係ですが、40Hzで‐6dBのハイカットフィルターを通しています。）
若干、小さくなりました。
以前、16個 スタックしてノイズを1/4にした記事も見かけましたが、そこまですると安価ではなくなります。
やはり、スマートな解決策が欲しいところです。
