pyGIMLi

ELSEVIER の Computers & Geosciences で 現在、DL数 ２位の文献です。

pyGIMLi: An open-source library for modelling and inversion in geophysics
https://www.sciencedirect.com/science/article/pii/S0098300417300584

物理探査のforward-inverse modeling が可能です。しかも、input ファイルの形式が簡単。この点で実務向きです。少ない労力で、すぐに答えを得られますので、数値実験にも適したツールです。スパースモデリングについては書かれていません。

まずは 2D を試行。弾性波探査です。
均質な速度層（500m/s）で forward-inverse modeling。観測結果にノイズなしの場合。

ノイズありの場合。結果も問題なし。

トンネルの事前調査では、ray paths が計画高を通るように配置を配慮してきました。が、通っていなくても直下により高速の層がないことは明らかですね。それはそれで有用です（今頃、何を言っているの？という感じですが）。

次は電気探査。均質地盤（100Ω・m）として forward-inverse modeling。ノイズ小の場合。

カラーバーのスケールを合わせ忘れましたが、数値を見る限り問題なさそうです。

次に、ノイズ大。

 

2D上の地形は補正されていますが、評価には注意を要する結果です。ま、河川堤防などで、どのあたりまで信用できそうかのチェックには使えそうかな？

YouTube では地下水を確認した場合などの境界追加について軽く説明されています。
https://www.youtube.com/watch?v=w3pu0H3dXe8
このような設定は未経験。軽い割には結果が劇的に変わります。今まで、電気探査を使えていなかったということでしょうか。

 

 

残念ながら、私の環境では 2D は動くものの、3Dは動きませんでした。

• Civil3D から 書き出した地形の stl が gmsh で読めない。（簡単な形状なら読める）
• iges は読めるが、gmsh から書き出したメッシュを pyGIMLi で読めない。（作り方による）
• stl を pyGIMLi でメッシュとして読めるが、座標を保持できない。（ stl の問題でしょう）
• そもそも、付属の 3D examples がエラーを吐く。

電気探査における側方地形の簡易補正量を試算したかったのですが、できませんでした。

ラッパー？の BERT も web 上で公開されてます。
他にも、多様なコードが無償で多数公開される時代。これから先のある若い方々には羨ましい時代になりました。

物理探査とスパースモデリング その 2

SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications
https://www.sciencedirect.com/science/article/pii/S009830041530056X
https://simpeg.xyz/

文献の掲載は2015年です。ブラックホール撮影のニュースより数年前ですね。
YouTube を見る限り、博論の一環として作られたように見えます。活気があって羨ましいような、それを提供できず身近の若い方々には申し訳ないような。

最適化ルーチンには、最急降下、L-BFGS、Gauss–Newton 法が実装されています。SciPy など外部ソルバーに接続するのも容易とあります。ソースを見ると確かに利用されています。

Regularization では Tikhonov、L2（Small）が適用できるようですね。文献では L1 もカスタマイズで対応可と書かれています。が、現行 Ver.0.15.0 の API 解説では、既にいくつかのスパース正則化が実装されています。

SimPEG.regularization.Sparse(mesh, alpha_s=1.0, alpha_x=1.0, alpha_y=1.0, alpha_z=1.0）

alpha って何？って感じですが、https://giftoolscookbook.readthedocs.io/en/latest/content/fundamentals/Alphas.html と同じ表記でしょうか。そうすると、norm=[[0,1,1,1]]で L1 でしょうか？

norms (Array): Norms used to create the sparse regularization, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*, *), Default: [[2. 2. 2. 2.]]

詳しい解説があれば私でも理解できそうですが、ついていけないですね。時間をかけても理解しましょう。

スパース正則化の利用例もいくつか示されています。
http://docs.simpeg.xyz/content/tutorials/05-dcr/plot_inv_2_dcr2d_irls.html
綺麗ですね。大部分のノイズがそげ落ちています。ま、どちらが適しているかは現場毎で他の情報と比較するしかないでしょう。

物理探査へのスパースモデリングの適用は自然な流れのようです。すぐに機械学習も加わるでしょう。見落とさないようにしないといけませんね。

 

**********************************************
20210604追記

SimPEG.regularization.Sparse の考え方は、以下の文献でした。

Inversion using spatially variable mixed ℓp norms
Dominique Fournier, Douglas W Oldenburg
https://academic.oup.com/gji/article-abstract/218/1/268/5420370

面白いですね。
いくつかの結果をすべて適度に正しいものと考えると、PCA などで特徴をさらに絞り込むようなアイデアは十分に理解できます。これを読むまで、私にはどれか１つを選ぶという発想しかなかったですね。
実装は前半までのようですが、それでも実務には十分です。

Adaptive Lp inversion for simultaneous recovery of both blocky and smooth features in a geophysical model
https://academic.oup.com/gji/article/197/2/882/622542 

However, in regions of complex geological structures, some physical properties may change dramatically across interfaces between distinct geological units. In this scenario, a minimum-structure inversion using sum-of-squares measure of model structure usually cannot reconstruct these sharp geological changes well. 

どのような地盤、対象でも同じインバージョン手法を用いてきた現状の方が乱暴なのでしょう。その視点では、mixed norms で選択肢のできるコードの方が地質屋さん向きです。

 

物理探査とスパースモデリング

ある方から声をかけられました。

何かな？と思って話を伺うと、「物理探査へのスパースモデリングの適用について」でした。恥ずかしながら、 そのような動向を知りませんでした。

調べてみると、ありますね。
珠玖「物理探査におけるスパースモデリングの可能性」
http://soil.en.a.u-tokyo.ac.jp/jsidre/search/PDFs/19/T-6-3.pdf

L1ノルムを利用したLASSO回帰です。
https://phreeqc.blogspot.com/2020/12/blog-post.html

物理探査学会でも平成31年2月のワンデーセミナーでテーマに取り上げられていました。
http://www.segj.org/letter/SEGJ_NEWS_NO42.pdf

全地連さんの「新マーケット創出・提案型事業」でも取り上げられていました。
https://www.zenchiren.or.jp/market/pdf/h27-3.pdf

スパースモデリング自体は昔からありますが、物理探査の最適化手法には組み込まれていなかったのですね。というか、トモグラフィーを自分で組んだことがない（市販品しか使ったことがない）ので、実装の種類について調べたり比較したりすることがありませんでした。
考えてみると、次元圧縮手法なのでノイズは飛ぶかもしれません。インバージョンに最初に適用されれた方は素晴らしいですね。

海外では物理探査の順・逆解析コードが多数公開されています。中には既に実装されているものがあるようです。話題になっていた機械学習手法も既にありそうですが。
日本では新商売を目指されていますのでコードは公開されないでしょう（少なくとも、最後の全地連さんからは公開されていません）。二番煎じなら特許を取れそうにないですし、残差を小さくするアルゴリズムを変更するだけなので（これまでより負担が増えるわけではありませんから）商売にならないのでは？ま、ヒトの世なので宣伝次第でしょうか。

最後の報告書を読みながら、海外のコードを調べてみましょう。

緊急事態宣言とバイアス

緊急事態宣言が各地で出ています。

個人レベルで対策可能なことは皆さん以前から実施されているでしょう。そういう点では宣言が出ても日常生活に変化がないかもしれません。

勤め先では、リモートワークに対応する支社と対面を続ける支社で分かれています。top executive の 考え方の違いです（ほとんどが前者ですが）。 

「マスクは意味がない」と言われていた方も、心境の変化かマスクをつけられる機会が多くなりました。ほとんどの方が常時マスクを外さないからでしょう。

一方、味覚異常、発熱、咳などの症状を有する方が PCR 陰性ということで休まず出社されていたり、それを静観する executives であったり。ここにきて正常性バイアスを振り払うのは大変なようです。

リスク管理や危機管理に対する個人の考え方や能力が露呈しています。よく見ておきましょう。

地下構造のモニタリング

波動を利用した地下構造のモニタリングにかかわる総論的資料２件。導入部分主体ですが基礎知識として有用な個所を ( ..)φ。

中原恒「地震波干渉法 その 2 応用」地震波干渉法による地下構造のモニタリング

実際には常時微動源の分布が等方的ではないことが多く，その場合は常時微動の相互相関はグリーン関数に対応しない．しかし，その場合でも，常時微動源が時間的に安定しているならば，常時微動の相互相関を用いて地下構造のモニタリングを行うことが可能であることが指摘されている［Hadziioannouetal.(2009)］．
波形から地下構造の変化を検出する手法は，すでに相似地震の解析のために提案されていたコーダ波干渉法［Poupinetetal. (1984)，Sniederetal. (2002)，Snieder (2006)］がよく使われる.
常時微動だけではなく，自然地震を用いることも可能であり，こちらの場合は常時というわけにはいかないが，自然地震が発生するたびにイメージングが可能となり，様々な時期の波動場の相互相関関数や自己相関関数を比較することにより，その期間内に地下構造の変化があったのかどうかを監視できる．
大・中地震の震源域では，2000年鳥取県西部地震に伴う変化を相互相関関数により検出したSawazakietal. (2009) や2011年東北地方太平洋沖地震に伴う変化を自己相関関数により検出したNakahara (2015) がある．

澤崎郁「大地震に伴う地下構造変化とその回復過程について」

地下構造変化に関する研究は，高密度かつ高精度な地震観測網の整備［例えば，Okadaetal. (2004)］と，コーダ波干渉法や地震波干渉法に代表される解析手法の高度化［例えば，Sniederetal. (2002)，Curtisetal. (2006)］に伴い，この十数年の間に飛躍的な進歩を遂げた．
大地震に伴う地震波速度の変化は，古くは断層運動に伴う断層面近傍での静的応力変化を反映するものとして解釈されることが多かった［例えば，早川(1951)，Poupinetetal. (1984)，Lietal. (1998)］が，近年は強震動に伴う地盤浅部の損傷が速度変化の主要因として有力視されている［例えば，Rubinstein andBeroza (2005)，PengandBen-Zion(2006)，Sawazakietal. (2006)，Wuetal. (2009a)，Takagietal.(2012)］．
強震動により地盤の物性が変化して弱震時とは異なる地盤増幅特性を示す現象は，地盤の非線形応答として知られているが，強震動が終息した後も地盤の物性が元に戻らない場合がある［例えば，AguirreandIrikura(1997)，PavlenkoandIrikura(2002)］．
強震動によりいったん低下した地震波速度は，数か月から数年をかけて経過時間の対数に比例しながら回復する［例えば，Rubinstein andBeroza (2004a)，Sawazakietal. (2009)］．
この対数比例型の速度回復は，岩石試料に振動を加えた後の固有振動数の変化を測る実験においても報告されており，岩石の物性が加振前の状態に向かいゆっくりと回復していくこの現象は，スローダイナミクス(slow dynamics)として知られている［例えば，TenCateetal. (2000)］．

２点間の伝播速度を（Active なり Passive なり）何らかの手法で求める、位相 or 群速度をモニタリングする、といった２段階の流れ。ambient noise だと長期間のデータが必要になるので時間の解像度が落ちます。ちなみに、Hi-net のデータを１年分DLしようとすると、１か所あたり丸１日かかります。そういう意味でも、余震のコーダ波を利用するのが現実的でしょうか。

passive 手法

passive 手法のテスト。

*** CC-FJpy ******************************
https://github.com/ColinLii/CC-FJpy
CC-FJpy: A Python Package for seismic ambient noise cross-correlation and the frequency-Bessel transform method
・複数の観測点間の相互相関関数をプロット。
・このコードに限りませんが、基本、長期データが必要です。F-netだと申し込み時に制限（１GB以下、リクエスト数（容量）制限）があるため、難しい。

*** noisepy ******************************
https://noise-python.readthedocs.io/en/latest/
・複数の観測点間の相互相関関数をプロット
・SAC を ASDF に変換して読み込み。
・noise_module が import できないエラー（S0B_to_ASDF.py）。

*** AANT *********************************
http://penglipk.github.io/AANT/
Adaptive ambient noise tomography and its application to the Garlock Fault, southern California
・sac 形式のデータを読み込み、速度を推定し、速度場を逆算してくれます（トモグラフィー）。
・cutsac のソースが含まれておらず中止。
・理想的なツールなのですが、更新が止まっています。非常に惜しい。

*** BIDO *********************************
https://staff.aist.go.jp/ikuo-chou/BIDO/2.0/bidodl.html
・時間領域と周波数領域の結果を比較する目的で実施。
・2.0 では 水平動も利用されていましたが、 3.0 では鉛直動のみ扱うとのこと。
・Ubuntu の binary が入っていましたが、動かなかったのでコンパイルしました。
・15分程度のデータであれば動作。速度が異常。
・1時間程度のデータは読み込み時にエラー。ディメンジョン修正箇所がわからず。

*** FMST *********************************
http://rses.anu.edu.au/~nick/surftomo.html
・トモグラフィーのみのようでした。分散曲線を出した後のステップです。
・速度を事前に用意することで、それらを満たす速度場を逆算。
・マニュアルを読みましたが、今回は手を出さず。

*** seismic-noise-tomography *************
https://github.com/bgoutorbe/seismic-noise-tomography
Rayleigh wave group velocities at periods of 6–23 s across Brazil from ambient noise tomography
・Python2.7でした。ライブラリがそろっていた 3.8 環境で動かし始めましたが、変換が面倒になり断念。

*** amb_noise_processing *****************
http://hestia.lgs.jussieu.fr/~boschil/downloads.html
Two-receiver measurements of phase velocity: cross-validationof ambient-noise and earthquake-based observations
・これもPython2.7でした。変換が print 文のみでしたので、3.8 にて実施。

この他、MATLAB を利用した GUI 対応コードもいくつかありました。さすがに購入しようとは思いませんでしたが、仕事で必要になれば選択肢に入るでしょう。

私の環境では、最後の amb_noise_processing のみ分散曲線まで確認できました。詳細を追いましょう。

************************************
20210516追加
「Two-receiver measurements of phase velocity: cross-validationof ambient-noise and earthquake-based observations」より、passive 手法で共通と思われる部分。

• Microseisms are most energetic in the period band between 5 and 30 s. At the shorter period end of this range, surface waves from teleseismic earthquakes are difficult to observe. Local to regional seismicity is needed to obtain phase velocities from earthquake data in the frequency range between about 3 and 10 s (Endrun et al. 2004).
• 1年のデータが必要（１年以上は不要）
• Records are filtered in the frequency band of interest (2–200s). サンプリング周波数1Hz以下⇒容量削減、計算負荷軽減。
• The data is then cut into overlapping time windows (e.g. half an hour with 50 per cent overlap), cosine-tapered at both ends.

 

CPS 3.30

 先週入れた  CPS 3.30 のコマンドです。

 ひとまず、動作確認ができたモノから記載。

*** 震源（dist）追加 *********************************************
$ gsac
GSAC> r *.U
GSAC> ch evla -4.52 evlo 101.374 evdp 34
GSAC> wh
GSAC> quit          

*** plot *********************************************************
GSAC> sort up dist
GSAC> p

*** mft *********************************************************
$ sacmft96 -f 'filename'.sac -PMIN 4.0 -PMAX 100.0 -a0 100.0 -A -VMIN 2.0 -VMAX 5.0 -U cm -R -S
$ plotnps -K < MFT96.PLT > temp.ps

*** pom ********************************************************
$ ls *.sac > cmdfil
$ sacpom96 -C cmdfil -PMIN 4.0 -PMAX 100.0 -nray 100 -A -VMIN 2.0 -VMAX 5.0 -R -S
$ rm cmdfil
$ plotnps -K -R < POM96.PLT > temp.ps

****************************************************************

（理解できていませんが）active 手法はひとまずクリア。あとは passive。

表面波の伝播 その2

先日の文献で紹介されていたツールです。

Computer Programs in Seismology 3.30
http://www.eas.slu.edu/eqc/eqccps.html

MFT が提供されていましたので、試しに使ってみました。

Cygwin 環境を作るのはイマサラなので、Linux 環境を選択。SSH や WSL2 だと GUI がダメなので、手元の laptop を Win10 との dual boot に変更（途中、BIOS を触って Win の boot loader が立ち上がるようになってしまいましたが、grub を選択して F10 で書き込めばOKでした（機種によって書き込み方が違うので、また迷わないように( ..)φ））。

これで Ubuntu 20.04 が入ったので、次はCPS。インストール手順は以下の通り。

$ sudo apt install make gcc gfortran

隠しファイル表示
/home/.profile にパス追加
PATH=$PATH:/home/CPS/PROGRAMS.330/bin
export PATH

このまま make すると ヘッダーファイルがないといわれます。追加。
Xlib.hの追加
$ sudo apt install libx11-dev
curses.hの追加
$ sudo apt-get install libncurses5-dev

$ ./Setup LINUX6440
$ ./C

これでOK。で、起動。
MFT は GUI を備えていましたので、まずこれから触ってみました（do_mft）。
残念ながら、Hi-net は読み込めず。F-net は読めました。同じ防災科研さんのツールで SAC 形式に変換していたのですが、前者についてはヘッダーを読めなかったようです。gsac では Hi-net も読み込み & プロット可能でしたので、do_mft の問題でしょうか。

続いて同じく GUI の do_pom。こちらは以下が出典。
Analysis of dispersive waves by wave field transformation
https://pubs.geoscienceworld.org/geophysics/article-abstract/46/6/869/68551/analysis-of-dispersive-by-wave-field
相変わらず理屈はついていけないのですが、結果は MFT よりマシ。先日の位相速度と同程度でした。うーん。

なかなか前に進みません。やはり、素人がすぐに結果を出せるほど甘くない内容です。

**************************************
20210511追記 その１

追加すべきライブラリは TOP ページに書かれていました。
http://www.eas.slu.edu/eqc/eqccps.html

Ubuntu Linux

You can use the program download page, or open a terminal window and enter the following commands:

sudo apt-get install make
sudo apt-get install gcc
sudo apt-get install gfortran
sudo apt-get install xorg-dev
sudo apt-get install libncurses5-dev
sudo apt-get install gawk

the following programs are useful:

sudo apt-get install gv
sudo apt-get install ImageMagick

追記 その2
F-net も読めないモノがありました。長いと読めないのでしょうか？

20210514追記
do_mft のバックで動いている sacmft96 などでは読めました。do_mft、do_pom のファイル認識に不具合があるようです。

 
20210626追記
Kimiyuki Asano et al.(2017)Surface wave group velocity in the Osaka sedimentary basin, Japan, estimated using ambient noise cross-correlation functions

The group delay time of the surface wave at each period was estimated using the multiple flter analysis technique (Dziewonski et al. 1969). 

MFTから求められる群速度。まだ理論は理解できていません。
http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/LessonA/index.html
http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/EMPIRICAL_GREEN/MFT.pdf

表面波の伝播

Hi-net、F-net から表面波を取り出せないか試行錯誤していました。

先日見た時間領域フィルタの文献では、遠地地震で綺麗に表面波が現れていました。
これはS波より遅れて到達する点、継続時間が分～時間単位になる点で、見た目だけでも同定は容易でしょう。が、震源が観測局に近接している場合や微小地震では判定困難。地震分野では、この表面波速度をどのように同定しているのでしょうか？

調べてみると、ありました。

まず、passive な手法。
微動＋相互相関を利用した表面波位相速度の推定です。この手法は図書にも載っていました。

・林田ほか「中京堆積盆地における表面波群速度の推定-Hi-net連続地震観測記録を用いた地震干渉法に基づく検証-」
https://www.jstage.jst.go.jp/article/zisin/66/4/66_127/_article/-char/ja/
・Bruno Goutorbe　et al.(2015) Rayleigh wave group velocities at periods of 6–23 s across Brazil from ambient noise tomography
https://academic.oup.com/gji/article/203/2/869/578005 

passive な手法は、土木分野でもハイブリッド表面波探査として利用されています。アレー探査から着想されたのでしょう。また、地震後の剛性回復チェックにも利用されているようです。
このような微動でレーリー波位相速度を求めておき、地盤構成を推定しておけば、地震が起こった際に各周波数で目立つ波がレーリ波かどうかチェックできるかもしれません。

次に、 active な手法。地震動を利用します。土木分野でよく利用されるアクティブ探査と同様の考え方です。違いは１イベントかつランダムな受振点配置。震央距離を利用するようですが、どのように群速度を計算するのでしょうか？

・Corchete V. et al.(2007) Shear wave velocity structure of the Sinai Peninsula from Rayleigh wave analysis
https://link.springer.com/article/10.1007/s10712-007-9027-6
・Dziewonski A. et al.(1969) A technique for the analysis of transient seismic signals.
https://pubs.geoscienceworld.org/ssa/bssa/article-abstract/59/1/427/101646/A-technique-for-the-analysis-of-transient-seismic

前者は以下の2種を組み合わせているようですが、頭がついていきません。
・Multiple Filter Technique (MFT)
・Time Variable Filtering (TVF)

GitHub 上で両者を扱うソースを順に見てみましたが、どれも一長一短。記載通り Python2.7 の環境を作っても動かなかったり、fortran だったり（最初はインタラクティブ希望）、面倒だったり。

ひとまず、弾性波探査のP波のピックアップと同じように、各周波数帯の波形から位相速度を取り出してみることにしました。
が、観測点が100、対象周波数が20あれば、2000個の波形に対してマニュアルピッキングになります。これは非効率。自動でできないかと obspy の Trigger/Picker をいくつか試したものの、どれもイマイチ（そういえば、表面波探査の解析では OYO さんのソフトで自動処理していましたが、どのようにピックアップしていたのでしょう？）。
結局、record section にして直線を引くことで平均位相速度を見出し、それを中心周波数に対してプロットしました。エンベロープを書けば、群速度も得られると思いますが、ひとまずここまで。

手を動かすほど、理解していなかった点、学ぶべき点が見えてきました。
機会があればプロに指導を仰ぎたいですね。

*************************************
20210508追記

現代地球科学入門シリーズ6「地震学」
p95-97 「6.6 表面波分散曲線の基づく推定」
・最終的に実施した手法
・鉛直成分でレーリー波、トラバース成分でラブ波を測定
p273-274 「21.2 常時微動の相互相関関数解析に基づく地震波伝播速度の推定
・passive
・地震直後の速度低下とその後の回復

C言語の基礎

GW 中に、C言語を眺めておこうと考えていました。

英語同様に、辞書があれば読めるものの全く書けない、という状況です。必要に迫られることが少なく、積極的に利用する機会はありませんでした。が、書ける能力があれば楽だったかも、ということは数回ありました。将来？に備えた弱点補強です。

取り急ぎ作りたいモノもないので、いくつかのサイトを見て１週間程度で完了できそうな以下を実施。

一週間で身につくC言語の基本
http://c-lang.sevendays-study.com/index.html

基礎編はほぼ復習になりましたが、応用編は知らないことばかり。あわせてメモしておきましょう。

基礎編
・void とは、関数の戻り値が無いことを意味する。
・Visaul Studio では、セキュリティの関係で scanf を使えない（scanf_s　ならOK）。
 プロジェクトを右クリック → プロパティ → c/c++ → SDL チェックを「いいえ」に変更
　※画面に直接入力する仕様は少ない（が、たまにあると忘れる）のでメモ。
・７日目、ファイル分割の御作法。
　#ifndef、#define、#endif マクロの意味、２重インクルードの防止、global 変数に対する extern修飾子。

応用編
２日目　変数のアドレス
・通常の変数：int a    値a    アドレス&a
・ポインタ変数：int* p    値*p    アドレスp
    p = &a:    pに、aのアドレスを代入    *pはaと同じものになる
    *p = 300:    *pに300を代入
3日目　ポインタと配列
・ポインタへのアドレス代入：p1 = &ar1[0]
    配列変数    アドレス    ポインタ    ポインタ変数の値
    ar1[0]    &ar1[0]    p1        *p1
    ar1[1]    &ar1[1]    p1+1    *(p1+1)
・ポインタ変数ｐ1に配列d[3] = { 0.2 , 0.4 , 0.6 }を代入：配列変数dの先頭アドレス&d[0]の値を代入するのと同じ処理
・動的メモリ確保 malloc()～解放 free()を使う場合、stdlib.hをインクルード
5日目　構造体
・構造体の成分の変数「メンバ」：(構造体変数名).(メンバ)
・データ渡し⇒スタック領域の圧迫、データコピー（無駄な処理）発生。⇒ポインタ渡しが通常

初心者がまず躓くといわれているポインタ。shallow・deep copy の関係と似ているようです。概念は理解できましたが、慣れるまでに時間がかかるでしょう。

丁寧な御説明で分かりやすく、2日で終わりました。感謝。
今回はココまで。