BERT 3D examples

先日、BERT （Boundless electrical resistivity tomography） の「3D examples を改良した」と言う旨の連絡が来ました。

この問題を10ヶ月前に Git の issue に投げていたのですが、その後音沙汰がなく、また simpeg へ移行していたため BERT 自体を削除していました。

それではもう一度試してみましょう、と思いインストール。
が、できません。手元の Ubuntu20.04 では以下のコマンドが機能しませんでした。
conda install -c gimli -c conda-forge pybert

では順にということで、まず仮想環境を作成しPython3.8 をインストール。pygimli, pybert, jupyter を順にインストールし、Git のファイルを一式DL。これで準備OK。
で、 testAll.sh を実行！
ダメでした。全てのexamples が failed！となります。jupyter 上では import できたので正しく入ったとは思うのですが。うーん

もう一台の Ubuntu20.04 で試すと、あっさりインストールできました。
こちらも公式指定の場所に examples がなかったので Git のファイルを一式DL。
で、 testAll.sh を実行！
動きました。

と言っても回答のとおり、slagdump は 走るものの 以下はダメなままでした。
2dflat/gallery:galleryGrid.cfg--preparation failed!
3dflat/gallery:galleryGrid.cfg--preparation failed!
3dtopo/acaur:acaur.cfg--preparation failed!

ま、Simpegがあるので、BERTは補助扱いで良いでしょう。

Joint Inversion

物理探査へのスパースモデリングの適用を調べていた際に、Joint Inversion というキーワードが目に留まりました。

いくつかの解き方があるようです。

複数の物理量を個別に最適化すると、それらの結果から統一された解釈を生み難くなります。本来、ある関係式で物理量が結びつけられているはずなので、それを定式化し inversion 時に考慮しよう、という流れが最も理解しやすいと思います。

pyGIMLi では、岩石の空隙率と飽和度を介し、弾性波速度と比抵抗値を結合しています。
Carsten Rücker（2017）pyGIMLi: An open-source library for modelling and inversion in geophysics

The Archie equation relates the bulk electrical resistivity of a medium　ρ to the fluid resistivity ρf and saturation S depending on porosity ϕ　(Archie, 1942):
ρ = a・ρf・ϕ^(-m)・S^(-n) (11)
To create geophysical parameter distributions, we apply common empirical petrophysical models, e.g., Archie's equation (eq. (11)) that provides the resistivity ρ as a function of saturation S. Sonic velocity v (or its reverse, the slowness s) as a function of porosity ϕ and saturation S is given by the time-average equation (Wyllie et al., 1956):
s = 1/v = (1-ϕ)/vm + ϕS/vw + ϕ(1-S)/va　(12)

petrophysical は貯留性の視点からみた（空隙に目を向けた）岩の物理量でしょうか。Archie equation はリザーバーとしての砂岩に対し報告されていた式のようです。結晶質岩ではありません。透水性との関係も図に表現されています。
https://onepetro.org/TRANS/article/146/01/54/161691/The-Electrical-Resistivity-Log-as-an-Aid-in

この関係式の利用は日本でも報告されています。
https://phreeqc.blogspot.com/2012/10/blog-post_23.html
文献の対象は花崗岩。当時、亀裂の影響について引っかかっていますね。ま、堆積岩や堆積物を対象にしたり、地球規模のマクロなスケールだと使えそうです。が、15年経過した現在、土木分野ではインバージョンにまで至っていません。浸透しない理由を挙げることは容易ですが、努力しない理由にはなりません（反省）。

先日の海外のセミナーでは、8つの物理量をインバージョンで扱われていました。いくつか結合できる点は想像できますが、8個全ては想像できません（他の流儀の解き方かもしれません）。
関係性をコツコツ研究されてきた昔の方々にあらためて敬意を表します。

電気探査 3D＋sparse

電気探査（比抵抗）の3次元処理（地形の考慮）をクリアしておこうと、SimPEG を1日触っていました。

といっても、良い例題があるので、入力ファイルとコードを読めばすんなり理解できます（以前に動作確認済みですので、細部を読むだけです）。

ついでに sparse 正則化を導入 。norms が効かなかったので、2.5Dの例題を移植。結果はあまり変わりませんでした。ま、これで電気探査の 3D + sparse はクリアです。

あとは可視化。これに時間を取られました。
HPの絵のように、地形＋電極位置＋断面のような絵を作りたかったのですが、SimPEG からメッシュ＋計算値を吐き出す方法が見つかりません。あーだ、こーだやっていても時間ばかりが過ぎていきましたので、slack の コミュニティー に投げました。

気づけば、30分強で返答がありました。
「この図はPyVistaを使った」でした。 活発ですし、制作側から返答をいただけるのはありがたい。
SimPEG では、メッシュ分割に discretize を利用していますが、そのメソッドの一部、VTK関連を有効にするには PyVista を入れないとダメ、という仕様でした。いや、わからないでしょう。どこかに書いてあるのでしょうけど。

PyVista をインストールして、例題のコードを参考にすれば、できました。断面の 「scalar＝表示したいarray 名」を見つけるのにやや時間がかかりましたが、比較的あっさり図化できました。何度も結果を確認しながら、ということになれば vtu 書き出し ＋ ParaVeiw より役立つかもしれません。ま、その書き出しも PyVista のインストールが必須だったわけですが。

PyVista で地味にうれしいのが、XYZ データから直接 TIN サーフェスを作成できる点です。matplot などではサーフェスを作成する前に Z データをグリッド化しないと描けないのですが、PyVista では座標データそのままでOK。 これは、使えます。

電気探査の3次元地形補正。長年放置してきましたが、ようやく正攻法でクリアできました（開発チームに感謝！）。

物理探査の教科書から

John M. Reynolds 「An Introduction to Applied and Environmental Geophysics　2nd Edition」の気になる箇所の備忘録です。地下水や物理探査、海外には良い教科書がありますね。

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

主な地質の弾性波速度 p147

主な地質の比抵抗 p291

主な地質の比誘電率、電波速度、導電率　p551

表面波

波の伝播イメージ Fig 5.41

分散曲線（高次モード含む）Fig 5.40-41

EM法

1次磁場と測定した2次磁場の比から見かけ導電率を推定　p434 Box 11.1

見かけ導電率と多層地盤モデルの関係　cumulative response function の利用 p434 Box 11.2 Fig 11.8

EMデータから1次元多層地盤モデルの構築　inversion の利用 p437 11.2.4.2 computer analysis

GPR

電磁波は地中で減衰する。

→海水の表皮深さは1cm。wet clay で30cm。 p540-541 

→表皮深さと比抵抗の関係　Fig 13.5

減衰した電磁波で測定・推定された導電率は、理論導電率よりもかなり小さい。

→小さな導電率で表皮深度を過大に見積もることに注意（実際は導電率が高く、表皮深さが小さい）　p541-542 Fig 13.6

表皮深さは地中レーダーの浸透深度とイコールではない。 p541

波長と鉛直方向の解像度の例 Table13.3

地下水位近傍で含水率が高くなると電磁波速度が落ちる→走時から計算した鉛直スケールが間延びする　Fig13.12

主な地質の探査深度と波長の関係（物理探査ハンドブックにも同図あり）p553

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

空中電磁探査で、表皮深度に1/2をかけて探査深度としている方がいらっしゃいましたが、案外、1次磁場の中間でないと2次磁場が減衰しきって返ってこれない、といった内容を考慮されていたのかもしれません。その場合、減衰した2次磁場から過大な比抵抗・表皮深度を求めないための補正というのも1手必要になるでしょう。

いずれにしても、上記の深度問題は単純な式で解けないようです。現状では分散曲線からの深度変換、EM探査からの深度変換ともに、逆解析を用いざるを得ないようです。それにはボーリング結果や検層結果との比較がより求められるのでしょう。p439には以下の通り記載されています。

However, there has yet to be an unequivocal demonstration of true depth discrimination being achieved using a fixed-dipole length, variable-frequency system. As with any modelled data, it is always sensible to test the results by using groundtruth information (e.g. from boreholes).

重心深度と見かけ比抵抗

表皮深さの1/2を空中電磁探査の深度として利用している文献がありました。

長谷川ほか「ヘリコプターを用いた空中物理探査データの再解析」
http://jolissrch-inter.tokai-sc.jaea.go.jp/search/servlet/search?5042189

ここからたどり着いたのが、こちら。
B. Siemon (2001) Improved and new resistivity-depth profiles for helicopter electromagnetic data
http://www.sciencedirect.com/science/article/pii/S0926985100000409

あるモデルでは表皮深さの1/2を深度として利用する手法がよくあいました、といったような内容です（あわないモデルもあるようです）。あう、あわないといっても、机上の話です。現地で確認しました、電気検層と対比しました、といったような内容ではなく、説得力を得られません。
1/2の理由はさらに先があり、その文献を取り寄せても、さらにまだ先がありました。ただ、この文献ではあわせこみパラメータとも取れる内容でしたので、深追いはしていません。


空中電磁探査での見かけ比抵抗等の求め方を詳しく知らなかったのですが、複数の方法があるようです。そのうち、表皮深さを用いた方法は Siemon による以下の2通りがメジャーなようでした。それ以外にも、√2で割る方法がありました。

1-1. R,Q→Da→da→Zsim

A = √(R^2+Q^2)
γ ＝ s/h
A'^（1/3)=A^（1/3)/γ
Da = s(A'/A)^（1/3)
Zsim = da + p/2

s: sensor horizontal distance
h: sensor height
p: skin depth

1-2. h,p→ρsim


2-1. calculating phase ratioε

ε=|Q/R|

R （real） ppm
Q (quadrant) ppm

2-2. log（ε）→log（A'^（1/3)）→Da→da→Zsim

2-3. log（ε）→log（δ）→ρsim

R,Q,h を 測定し、ρsim と Zsim を求める際に p （と係数）をパラメータとして導入しているように見えます。たとえ1/2や表皮深さの導入に理論的背景があったとしても、測定対象と整合しなければ他の方法を選択できるというのは、理論がまだ確立されていない、推定するのは難しいということなのでしょう。
それに、以前にも書き残していますが、実際に磁場が到達しているのかどうかは確認しないとわかりません。また、1つ目の文献のように、得られた信号がノイズレベル以下になっていないかどうかもチェックが必要でしょう。

そうなると、現段階では電気検層など他の手法と空中電磁探査結果の対比、その補正が必須となります（物理探査ですので、当然なのですが）。見た目の綺麗な断面図に惑わされがちですが、別の物性値と対比がなされているか、推定式の妥当性は確認されているか？といったステップを踏まないといけないのでしょう。

新旧調査手法

砂防学会の発表予稿を通して見ていたのですが、以下の3点が印象に残りました。

・空中電磁探査

・酸素同位体比

・UAV

上記調査が増えたように感じましたね。深層崩壊を通じて砂防分野に興味を持ち始めたのが最近ですので、ようやく気付いたのかもしれません。

酸素同位体比は、地下水調査で他の同位体と併せて実施されています。UAVについても比較的安く導入できます。空中電磁探査は前回記載した通りです。これだけは何とかしたいところです。

一方、最も興味を惹かれたのは、古典的な堰堤の老朽化調査でした。

数点ほど発表されていましたが、いずれも表面だけでなく、トモグラフィー的な弾性波測定を計画・実施されていました。これ、地味ですが解像度はかなり上がるでしょうね。

新しい調査法を導入する必要もありますが、古いものをより深く使えるようにすべきと、あらためて感じた予稿集でした。

空中探査

ヘリが何か運んでいたので見ていると、空中探査でした。電磁探査でしょうね

今年の砂防学会でも、いくつかの発表があり、気になっていました。
実際に見たのは初めて。昨今、UAV も実務に適用され始めているのですが、思わず、周りの方と一緒に写真を撮っていました。

さすがに、早いですね。山をなめるように数往復し、見えなくなりました。往復する密度も決まりがあるのでしょうね。あのスピードなら、きっと周辺の山々も一緒に測って、今後の営業材料に使うのでしょうね。

トンネルの事前調査でよく実施される2次元比抵抗探査ですが、最近は3次元地形で補正しないとダメだと思っています（比抵抗は当たらないといわれる一因でしょう）。3次元の測定ができればBESTですが、時間がかかるのが難点。また、崖など人が入れない場所は測れません。その点、空中探査なら問題なさそうです。深度はどの程度測れるのでしょうか？
出遅れ感、満載ですが、基礎知識だけはつけておく必要がありそうです。

テスターによる抵抗測定

踏査で電気伝導度計を持ち歩くことがあります。

地表水や湧水を対象とするのですが、今回の目的は汚染の有無でした。

現場道具の中にテスターも常備しているのですが、これなら軽く、思いついた時に測れて便利だなあと思い、とりあえず電気伝導度計と一緒に抵抗を測っておきました。

帰社してから、さあ、どのように電流の通過した面積や体積を考えれば良いか？などと考えて資料を探していると、同じような事を考える方がいらっしゃいました。以下の資料ですが、積雪の電気抵抗をテスターで測定されていました。測定中の不安定な挙動も同じですね。
松岡,春樹 (1969) テスターによる積雪の電気抵抗測定, 福井大学工学部研究報告 3月

幸い、テスターをまったく同じように使っていましたので、掲載されていた式で抵抗値を電気伝導度にそのまま変換できました。結果は以下の通り。

データ数が少ないので、相関性が有意かどうかは分かりませんし、テスターでは正確な値は分からないでしょう。しかし、高低の差は出ます。精度は落ちても数が必要な時には良いかも知れません。覚えておきましょう。

比抵抗と地山分類

以前、お客様から、「比抵抗探査結果は地山分類に使えるの？」と直球を投げられました。

2次元比抵抗探査は古い手法だと思うのですが、屈折法弾性波探査のように、地山等級に直結して評価された例は少ないようです。

鹿島出版会「地盤の可視化技術と評価法」には、その例が若干載っています。「トンネル毎に差があるなー」などと思ってみていますと、それとは別に、ちょっと変わった評価法が載っていました。
間隙率、飽和度を弾性波速度、比抵抗値と関連付け、それらを連立させて同定するというものです。面白いなーと思って出典を見てみました。
https://www.jstage.jst.go.jp/article/jscejf/62/4/62_4_603/_article/-char/ja/

コアの試験値との関連づけであるため、式中の亀裂の影響の取り込み方が明確でない点のみ引っ掛かりました。ただ、施工時の支保実績との間に良い相関が認められるなら、実務的には無視しても差し支えないでしょう。参考書に載っていますので、ある程度実証済みなのか、データを集めている段階なのか分かりません。ま、この手法に限らず、比抵抗の利用について今後の展開に期待しましょう。

電気探査

先日扱っていた HandyViewer miniOHM では、いくつかのアダプターが接続できます。
そのアダプターによって比抵抗値の計算方法は変わります。が、基本的な考え方は同じでしょう。

電極を多数配置して測定する電気探査の場合も、電流が放射状に流れ、等電位面が球面を形成するというイメージではないでしょうか。地表から電流を流し地表で測定する場合は、放射状電流の重ね合わせを考え、平均的な比抵抗値を算出します。地下水で言う虚像井と同じ組み立て方です。上記の壁面アダプターも同じ考えでしょう。

電気探査の場合は電極間隔（球面直径）を大きく取り、深部の情報まで得ようとします。それを、解析的にフィッティングし、（平坦であろうが斜面であろうが）鉛直2次元として作画するわけです。少しおかしいのですが、平坦面でも鉛直、斜面でも鉛直で作画します。構造物や技術者に都合の良いように、解釈、投影するわけです。ま、その程度の精度であるということです。

時々見るのですが、複数の直線で曲線を分割した測線計画。探査深度と電極間隔による球面のイメージ内であれば、多少の曲がりは個人的にOKなのですが、直線で計画される方はそうではないようです。それなら鉛直かどうかにこだわった方が良いと思うのですが。
では、どこまでのRであればOKなのかというと、基準はありませんし、個人的にも感覚でしかありませんでした。

先日、このような話をダラダラとプロにしていると、ある指標を出してこられました。Rなりの電極間隔（探査深度）と、それを直線で結んだ場合の距離の比率を、解析に対する影響度として定義されていました。ま、何もないよりはマシといった程度ですが、「感覚」よりは説明しやすいかも知れません。機会があれば利用させてもらいましょう。

比抵抗値の整理

先日測ってきた岩盤露頭や河川水の比抵抗値を整理していました。

使用した機械は OYO さんの HandyViewer miniOHM Model-2121 です。15年前の機械。古いので内部電池は切れています。また、現場で計測していると、接触が悪いのか時々電源が落ちます。

帰ってから分解して振動を与えてみると、どうもコネクタ類ではなく、液晶周りか電池ボックス。電池ボックスを見ると、端子の一部が押し出されておらず、電池にスキができる状態。これが原因でした。端子をハンダで盛り付けし、スキをなくすことで解消しました。
一部は押し出されていたので、この端子、不良品ですね。15年間、使われた方はおかしいと思われなかったのでしょうか？

次に、結果を読み込もうとしましたが、整理ソフトがありません。15年間、誰も整理しようとしなかったようです。ま、比抵抗値さえ読めれば良いので、問題はなかったのかもしれません。
今回は数ヶ月前に SEISMOGRAPH Model-1816 のデータ整理用EXCELブックを作っていましたので、それを流用することにしました。データ形式は似たようなものなので、読み込み部分はすぐに加工できました。
が、読み込んだ結果を見ると、電流・電位ともに波形の値とヘッダーの値があっていません。波形からの算出法はマニュにも書かれていません。仕方ないのでOYOさんに聞くことに。

帰ってきた答え、マニュにまったく書かれていない演算内容でした。
その演算内容を取り込むと、確かにヘッダーの値になります（丸め誤差で微妙に違う場合もあります）。マニュアルとして、値として微妙ですが、ま、良しとしましょう。

整理用BookをUPしました。コチラ↓
https://sites.google.com/site/geochemist001/resources/handyviewer_mo