DASのピッキング

Seismic arrival-time picking on distributed acoustic sensing data using semi-supervised learning | Nature Communications

DASのピッキングに関する文献です。
従来の STA/LTA 等ではうまくいかないDASのピッキングに機械学習を使ったという内容です。

この文献、査読結果が一緒に公開されています。オープンサイエンスもここまで来たのだなと感心しました（他の文献でも公開されていますので、私の気づきが遅かったのでしょう）。

その中の指摘で、近年の文献を見ながら感じていた部分が明示されていました。

I do understand that deep learning is currently a big hype and that it is tempfing, especially for young researchers, to become part of this. However, the days where the mere applicafion of some deep　learning model was science, are over. With numerous easily usable programming tools, the training of a neural network (deep or not; this is just somewhat hollow semanfics) has become an easy task for undergraduate students. After all, it is what it is: fifting the coefficients of some funcfional form to a collecfion of data.
This said, it is not clear how exactly the authors go significantly beyond just training yet another deep neural network? What is the innovafion that goes beyond the obvious?

LSM作成でも、機械学習分野でSOTAを達成した手法の流用は見ていて面白いのですが、それまでですよね。できたLSMが実際にどの程度使えるのかがわかりませんので。GISを使って特徴量を作り、機械学習にかけるというのは、学生でもできます。今年のGCONでは高専生も頑張っています。
第3回高専GIRLS SDGs×Technology Contest(高専GCON2024)

ピッキングも同じなのでしょう。他にもいくつか機械学習を利用している文献を見かけます。が、機械学習を使えば精度が上がった、だけの時代はもう終わったと認識せざるを得ません。

土砂災害と地質のスケール

第43回地質調査総合センターシンポジウム「地質を用いた斜面災害リスク評価－高精度化に必須の地質情報整備－」で千木良先生の発表を拝聴しました。

• 能登地震でのある崩壊の素因として、パミスを含んだ厚さ５㎝程度の層が挙げられる
• しかし、産総研の地質図に描かれるスケールではない

これまで全国の災害データの地質的素因をシームレス地質図を利用して分析してきましたが、がけ崩れ、土石流には大きく影響しませんでした（地すべりには地質が効きます）。純粋な地質のみでなく、複雑性や断層からの距離等の特徴量も作成しましたが、機械学習モデルが重要とは判断しませんでした。これは、崩壊データではなく、災害データだから（人の活動が影響する）というのが主要因ですが、上記も大きく影響しています。地質は影響する。しかし、あなたの地質図の選択が悪い、という事実を突き付けられているのです。

１つの崩壊、すべり毎に素因を抽出し、それを地道に、年月をかけて積み重ねる。数が集まったら、そこから共通性や特異性を見出すのが研究の基本手順なのですが、近年の研究ではその過程をすっ飛ばし、結果を早急に示しがちに思えます。手元にある災害データ、広域で使える均質な地質図、それらから得られた結果はどれほどの信頼性があるのでしょうか。「ないものは仕方がない」ではなく、ないから作る、積み重ねるのが研究であり研究者でしょう。目に見える成果や時短に迫られる立場もわからないではないのですが。

Landslide Sasceptibility Map を広域で作るのが他国で流行っています。「ひとまずやってみた」レベルが多いのですが、手法は着実に進歩しています。それに投入する特徴量の質が上がれば、信頼性もぐっと上がるでしょう。将来の方のために地味な仕事も積み重ねていきましょう。

直線性

地震波の直線性、平面性、方位角、傾斜角について知ったのが今年に入ってから。
プロから文献を紹介していただき、いくつか事例も見ましたが、なかなか文献に紹介されているような特徴は見出せませんでした。

河川流量の推定に振動を利用されている例を見ますが、以下では直線性、方位角も示されていました。
Exploring the relationship between seismic noise signals and modeled river flow data: A case study from Sicily, Italy - ScienceDirect

These analyses show that the scattered azimuth and the low rectilinearity values, observed at the beginning and the end of Helios, coincide with periods of increased wind speed (Supplementary Figs. 5a and b), which can generate a seismic noise as a consequence of the impact of the wind on the trees or buildings. Indeed, the wind seismic noise is usually characterized by spherical waves showing a low value of rectilinearity (<0.5) and a lack of horizontal trend values (Zheng and Stewart, 1992; Panzera et al., 2016).

振動を理解するには環境計測、音波等も含め多面的に捉える必要があるというのはつくづく感じているのですが、なかなかこのような結論に結び付きません。まだまだわからないことばかりです。

 

３次元地質モデルと機械学習 その3

Progressive Geological Modeling and Uncertainty Analysis Using Machine Learning

2段階で機械学習による層序区分、岩石区分を推定し、分布確率を同時エントロピーとして表現されています。1段階目では XYZ 座標を入力として層所区分を予測、2段階目はそれらと物性値2種を入力として岩石区分を予測。手法的には半教師あり学習になりそうです。

1段階目をミスると2段階目も誤りになります。が、地質と岩級の両方を扱いたい、などという場合に使えそうな提案です。

機械学習による確率分布の算出には、これまでの文献のようにボクセル毎で正規化するのではなく、ロジット値を算出しておいて後からモデル全体で正規化すべきだと考えています。が、2段階の場合はこの文献の求め方のほうが好みです。

確率分布の表現も複数あるので、目的に応じて使い分けるのが良いのでしょう。

３次元地質モデルと機械学習 その2

A stacking methodology of machine learning for 3D geological modeling with geological-geophysical datasets, Laochang Sn camp, Gejiu (China) - ScienceDirect

一歩進んで、地質、磁気探査、重力探査結果を組み合わせた3D地質モデルです。

 The input data are coordinate (x, y, z), magnetic susceptibility and residual density

Joint inversion とは異なり、これも物理的関係は考慮されていません。機械学習としてブラックモデルが分布確率を提示します。データ駆動型の手軽さと解釈の難しさを併せ持っています。
手法としてはアンサンブル（スタッキング）で判別性能を高める工夫がなされています。

日本でブラックボックスモデルは受け入れられるでしょうか？
近年の文献では、データ駆動型の方が物理ベースのモデルに比べて性能の良くなる傾向が認められると言われています。が、地質モデル作成の場合はトレーニングに利用できるボーリング数が補間するボクセル数に比べ、圧倒的に少なくなります。その場合は物理ベースモデルを利用し、その結果を地質屋さんが解釈し分布を推定する方が性能が良くなるでしょう。
その「解釈」部分のロジックを、先の文献のように損失関数に組み込んだペナルティとして数式で表現すれば、機械学習でも良い地質モデルを作れるようになると思われます。プログラミングのように、自分が何を考えて推定したのかを整理、伝達することから始める必要があるということです。

３次元地質モデルと機械学習

地質の３次元可視化にも機械学習が使われつつあります。

良い点は、推定地質の分布確率を表示しやすいこと。いえ、統計処理した確率ではありません。モデルが考える確率ということで、その根拠は曖昧です。
これまでは indicator kriging のように地球統計学を利用していました。が、それより手軽に算出できること、機械学習ユーザーの方が圧倒的に多いことから、今後事例が増えると思われます。

まずは、コチラ。2024年です。半教師あり学習を利用しています。

従来と同じくソースは XYZ と地質情報のみですが、線形補間しているようなリサンプリング手法を用いて大きなエラーが生まれないように工夫されているようです。
GMD - GeoPDNN 1.0: a semi-supervised deep learning neural network using pseudo-labels for three-dimensional shallow strata modelling and uncertainty analysis in urban areas from borehole data

次は損失関数を工夫した例。これも2024年。
地層の逆転がない場合、古い地層は必ず新しい地層よりも下位に分布するという単純な法則に対し、違反した場合にペナルティとなるよう損失関数に組み込んでいます。これだけの工夫ですが、うまく予測してくれる場合があるようです。PINNs から着想を得られたのかもしれません。
Research on 3D Geological Modeling Method Based on Deep Neural Networks for Drilling Data

いずれもまだ「やってみた」レベルですが、そのうちより良い推定ができるようなモデルの構築方法へ議論が進むことでしょう。一方で、統計処理や機械学習ができるほどボーリングを掘るサイトは少ないので、身近な土木分野では大きくは変化しないかもしれません。大きなサイトでどうなるか、今後に着目しましょう。

AI-driven rapid landslides mapping

NHESSD - Brief Communication: AI-driven rapid landslides mapping following the 2024 Hualien City Earthquake in Taiwan

地震時崩壊個所を、2手法で迅速に特定したという報告です。
ソースと機械学習手法の組み合わせは以下の通り。

• Sentinel-1 (SAR) - CNN
• PlanetScope（可視画像）‐ ViT

 The ViT model was pre-trained and validated on a multi-source landslide segmentation dataset (Fang et al., 2024), the Globally Distributed Coseismic Landslide Dataset (GDCLD). 

データ入手に日数を要していますが、作業自体は 20分～2時間程度だそうです。 
ViTの事前学習に利用したデータは公開されており、またCNNコードも公開されています。残念ながら私の環境では素直に走らなかったので、どの程度の成果が得られるのかは確認できませんでした。

能登の地震災害では、２日後に地理院で航空写真が公開され、3日後に崩壊個所の判読結果が公開されています。提案手法の場合、精度を保てないと現行手法と勝負にならないのですが、SARでは天候に大きく影響されない点や、広域でも一人で対応できる点が魅力です。そのうち性能は向上するでしょうから、災害直後の道路ネットワークの寸断箇所（孤立集落）の確認や天然ダム形成の確認に使える日が来るかもしれません。

航測会社も頑張って開発されていると思われますので、この分野で良い勝負がなされるよう期待したいところです。

LSM の信頼性

Landslide Susceptibility Assessment in the Japanese Archipelago Based on a Landslide Distribution Map

日本の研究者が日本の Landslide Susceptibility Map を AHP を用いて作成されています。
崩壊データとしては、「地すべり地形分布図」を利用されています。主として30～40年の期間しか扱えない災害データよりも、数百～数千年程度のイベントを反映した地すべり地形分布図の方が、susceptibility を見出す点で適している部分があるでしょう。
また、国外の研究者が作成するよりも地質を細かく取り扱われています。このあたりは国内の研究者の方が有利でしょう。このような文献を見ると、安心します。

この文献もそうですが、LSM を作成している文献では、基本的にはオープンデータが利用されています。検討ツールは異なるものの、どれも同じようなアプローチです。既に LSM の作成方法や考え方はある程度確立されていると言えるでしょう。
問題は、どの作り方、どの特徴量の選び方が良かったのか比較する場がないことです。文献では、当然「自分たちの提案手法は良い性能を出したよ」「今後も改良していくよ」で終わっており、他の文献と同じデータを使って手法を比較する、優位性を示すまでには至っていません。画像分類における ImageNet のように標準的なデータセットがなく、ILSVRC のように比較する場もありません。このようなデータセットや場の提供を学会等が担うべきだと思います。

機械学習を使えば信頼性の高いマップができるわけではありません。研究者が切磋琢磨できるような場ができて初めて、マップや機械学習の信頼性が向上し始めるのです。まだまだこれからです。

タンクモデルのパラメータ最適化 その２

タンクモデルに関する問い合わせがありました。

最近は機械学習を利用する印象が強く、タンクのことを忘れかけていました。
振り返ってみると、EXCEL、Fortran、Python でパラメータ探索を実施しています。なんだか苦労して重み（ペナルティ）を手で合わせた記憶がありますね。

機械学習のノウハウがある程度身についた今だと、タンクモデルでも同じフレームワークを使用できるなあと思いつきました。

畳み込みブロックを何段にするか --> タンクを何段にするか

ハイパーパラメータの最適化 --> パラメータの最適化

損失関数の選択 ＝ 損失関数の選択

機械学習で利用している Optuna を使えば、容易に最適化できそうです。同定するパラメータや計算量が機械学習に比べて非常に少ないため、探索に時間もかからないでしょう。同じようなことを考える人がいるのでは？と調べてみると、やはりいらっしゃいました。
単流域型タンク・モデルとニューラルネットワークの比較

DNN や GBM だと、過学習を起こしやすい印象を持っています。また。外挿も苦手。タンクモデルくらいの少ないパラメータで長期のデータを扱う、交差検証を取り入れるなどの工夫が過学習を起こさないちょうど良いレベルなのかもしれません。ペナルティの配分の仕方も何かしら自動化できそうです。

「丁度良い」最適化ができそうな気がします。

inter-aquifer connectivity

Using geochemical and geophysical data to characterise inter-aquifer connectivity and impacts on shallow aquifers and groundwater dependent ecosystems. - ScienceDirect

データ量が多く全容を理解できませんでしたが、調査法として参考になる文献でした。これだけ調査するのは大変だったでしょう。

• ピットを掘るので地下水を低下させる。
• ８㎞西に文化的にも重要な湧水群がある。
• その周囲には貴重な植生もある。
• ピット周辺の帯水層と湧水群の帯水層は別。
• 間に厚い不透水層がある。

このような調査結果をシミュレーションに与えると「大きな影響なし」という答えが出てきます。私が担当していたとしても、恐らくそのように結論付けたと思います。
が、実際は西方の井戸に地下水位の低下傾向が認められた。
そこで追加調査を実施。

• 不透水層は部分透水だった。（透水係数、若い地下水（年代測定））
• 上部の風化帯を通じて帯水層がつながっている可能性あり（メジャーイオン、放射性同位体組成の重なり、透水係数、若い地下水（年代測定））

空中電探、水質分析、地下水年代測定等を実施することで、地下水流動経路と影響の原因を特定した、という内容でした。

国内では、建設範囲外でこれだけ広域の調査を実施することはレアでしょう。温泉や重要水源への影響検討の場合は実施する場合があるものの、空中電探まではしていないでしょう。が、やるべきなのでしょうね。いろいろな可能性を考えて調査を提案することが、後のリスク低減に繋がります。それが重要である点をこの文献は示唆しています。
実務のお手本になるような文献でした。

VisionTransformer

先日書き残しましたが、LSM 作成に transformer が利用されていました。
https://phreeqc.blogspot.com/2024/11/landslide-susceptibility-map-using-ml-5.html

transformer を利用したいと思っていたのですが、実務では画像分類での ViT を先に触りました。結果は他のアーキテクチャを抑え Best Score。また一つ、外せない選択肢が増えました。

先月、ViT の図書が発売されました。おそらく、国内では2冊目でしょう。先の本もわかりやすいのですが、こちらも good。変化するテンソルの大きさが良くわかります。CNN同様に、これから解説本が増えていくかもしれません。

「Vision Transformer／最新CNNアーキテクチャ 画像分類入門 」

図書には transformer 以降に提案されたネットワークアーキテクチャがいくつか整理されています。2020年末から２～3年で多くのネットワークが提案されているのは、CNNの頃と変わらないでしょうか。

CNN とのハイブリッド手法については、まだ触れる機会がありません。遅れないようについていきましょう。

国内のLSM、ハザードマップ、リスクマップ

日本では、「Landslide Susceptibility Map」に対応する統一的な用語はまだ確立されていません（産総研さんは「感受性マップ」でした）。
土砂災害に関するマップは古くから作成されていますが、様々な表現が用いられています。「ハザードマップ」と「リスクマップ」についても、明確に区別されずに使用されています。

崩壊の危険度に関するマップとして、まずは深層崩壊。国交省から公開されていますが、ここでは隆起量や地質を考慮されているようです。豪雨時のハザードマップの材料として利用可能でしょう。

●深層崩壊推定頻度マップ
過去の発生事例から得られている情報をもとに深層崩壊の推定頻度に関する全国マップを作成しました。（2010/08）
●渓流レベル評価マップ
空中写真判読等による深層崩壊の渓流（小流域）レベルの調査。（2012/09）
https://www.mlit.go.jp/mizukokudo/sabo/deep_landslide.html

地震時の斜面崩壊に対するハザードマップもあります。
この六甲式は評価手法になるのですが、得られるのはマップです。シンプルな式ですが、経験的に良い性能を発揮します。
●六甲式
国総研資料　第 204 号
参考：地震による斜面崩壊危険度評価判別式「六甲式」の改良と実時間運用

さらに一歩進んだリスクマップも公開されています。国内では崩壊データが整備されていませんので、土砂災害データが利用されていたり、民家の有無も考慮されたうえでの評価がなされています。
●土砂災害発生確率マップ（案）
国総研資料　第 1120 号
●ハザードマップポータル
ハザードマップポータルサイト
●キキクル
気象庁 | キキクル（危険度分布）

現状、機械学習を導入するだけで良いマップを作ることができるとは考えていません。おそらく、よく検討された六甲式や土砂災害確率マップの方が実態に近い結果を得られると思います。
将来、過去の崩壊データが整備される、あるいはリアルタイムで整備されるしくみが整い、人の手に負えなくなる量が毎年加えられる、というようなことになるかもしれません。そのような場合に、重回帰に代わってDNNを使う、適中率と捕捉率の同時向上＝F1マクロを評価指標としてトレーニングする、などという取り換えが可能となり、良質なマップができるようになるのでしょう。

国内では、機械学習による LSM 作成を試しながらも、データ収集・提供の仕組みを整える必要があるのでしょう。

Landslide Susceptibility Map using ML その6

天然ダムによる河道閉塞が発生するか？まで考えられた文献です。

A global-scale applicable framework of landslide dam formation susceptibility | Landslides

イタリアと日本全土が検討済みです（著者は中国の方のようです）。LSMを作成した後、川に重ねて河道閉塞リスクを検討されています。道路への影響検討を扱った文献はいくつか見てきましたが、天然ダム形成については初めてです。（３）式がどこまで実態に即しているのかはわかりませんが、一歩進んだ取り組みでしょう。

Fig. 6

Landslide susceptibility evaluation result in Italy (a) and Japan (b)

Fig. 8

Landslide dam formation index result in Italy (a) and Japan (b) (each dot with LDam formation index represented the centre of a river reach presented as 90-m resolution grid in MERIT Hydro; Figs. S1–8 in supplementary presents two instances of visualising the LDam formation index on a smaller scale)

内容や精度はともかく、「ひとまず全国やってみた」「使ったデータは公開したよ」という取り組み姿勢は素晴らしいと思います。日本の税金で整備したデータを他国の研究者が使用して成果を公開するという点に引っかかりがありますが、それも日本の研究者の刺激になれば良いでしょう。

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

20241128追記
川幅は日本の方が作られた全球データを利用されています。落としてGISでプロットしてみましたが、粗くて途切れていました。上記のマップはかなり引いていますので、雰囲気を比べられる程度でしょうか。

Landslide Susceptibility Map using ML その5

Improving landslide susceptibility prediction through ensemble Recursive Feature Elimination and meta-learning framework | Research Square

2024年の文献です。
特徴量選択にアンサンブル（多数決）を利用しています。個々のモデルでは、RFECVを利用。さらに、ベースモデル ‐LRによるスタッキングを通して最終的な出力を得ています。

Full article: Near real-time spatial prediction of earthquake-induced landslides: A novel interpretable self-supervised learning method

こちらは transformer を利用されています。
pre-training 後に fine-tuning を実施したところ、他の手法より AUC がよかったよ、という報告です。

transformer ではグローバルな発生データを活用する pre-training が可能であり、データ量の多さを活かして高度な特徴学習ができそうです。それをローカルの発生・非発生データで fine-tuning することで、未知のデータに対する汎化性能を保ちつつ、ローカルな特性を捉えたモデルを構築できます。ローカルのために世界のデータを利用するという報告は"その2"で書き残しました。が、この場合は XGBoost ですので pre-training の概念がないですし、学習時に利用するにしても相応の非発生データが必要になります。

利用する特徴量を決めておいて、世界でデータを整備しておけば、事前学習済みデータとして配布・利用できそうです。幸か不幸か国内はこれからですので、特徴量に使える国内データの整備が進むとありがたいですね。

Landslide Susceptibility Map using ML その4

Full article: An integrated neural network method for landslide susceptibility assessment based on time-series InSAR deformation dynamic features 

時系列 DInSAR を特徴量として使用されています。変動量はSBASでもチェックをされているようです。SARというと、つい地震前後の差分をイメージしてしまいますが、地震前の変動量を利用することは言われて初めて気づいた重要なポイントですね。

24 stages of time-series InSAR cumulative deformation information are taken every 96 days per quarter. 

国内で Landslide Sasceptibility Map を産総研さんの研究以外で聞いたことがありません。それを作成するための機械学習も土木分野では浸透していません。あっても SVM とか Random Forest などクラシカルな手法が使われているように感じます。が、この分野の研究では複雑な特徴に対応するため、または精度向上のために DNN 等が利用されています。

Based on traditional linear statistical analysis, machine learning methods stand out by virtue of their ability to examine large amounts of data independently. Machine learning methods, such as random forest (RF) (Dou et al. Citation2019) and logistic regression (Zhang et al. Citation2019), have been extensively used in LSA. However, under the requirements of complex scenes or high precision, traditional machine learning algorithms cannot meet actual demand (He et al. Citation2021a; Zhao et al. Citation2022). Building on the neural networks present in machine learning, the neural network method effectively predicts complex nonlinear dynamic systems. It has been widely and successfully introduced into the field of LSA, including the convolutional neural network (CNN) (Wang, Fang, and Hong Citation2019; Gao et al. Citation2023a), recurrent neural network (RNN) and deep belief network (DBN) (Chen et al. Citation2020). The convolutional layer of CNN can extract multidimensional features from the input images and has good performance (Hakim et al. Citation2022). Gated recurrent unit (GRU) network of RNN variant has good performance in processing sequence data (Zhao et al. Citation2022). With the complexity of the environment, when faced with a limited sample, the ensemble learning model is also widely used in the LSA. For example, Wang et al. (Citation2022) conducted the LSA based on the XGBoost ensemble learning model. Lv et al. (Citation2022) combined CNN, DBN and ResNet models with the ensemble learning techniques of Stacking, Bagging and Boosting to generate the LSA.

日本は２ステップ遅れている状況です。せめて LSM 作成のためのデータは整備していただきたいものです。そうすると土木に携わっていない機械学習エンジニアが参加しやすく、多様な目的のマップが作成されるようになるでしょう。その時々の、最新のアーキテクチャで。

Landslide Susceptibility Map using ML その3

 昨日、産総研さんの第41回地質調査総合センターシンポジウムを拝聴しました。

地質に関するCIM やデジタル化の話、九州のハザードマップ作成についての報告がありました。

ハザードマップは産総研の方の報告でした。その中には LSM の話もありました。産総研では LSM を「地すべり感受性マップ」と訳されているようです（が、この訳は一般的ではないとおっしゃっていました）。スライドには「崩れやすさマップ」という記載もありましたので、今後はわかりやすい後者に落ち着くのかもしれません。

LSM 作成の話がメインのように思えましたが、表題はハザードマップで、しかも使用されているインベントリは被災を伴う災害データベースであり、作成されたマップはリスクマップに近いようでした。「感受性とハザードはどう違うのか」と質問された方がいらっしゃいましたが、「いつ、どこで、まで答えるのがハザード」というような回答でした。また、災害データを使われた理由については、「被災の有無にかかわらない崩壊データも一部で使っている。ないものは使えない」といった回答でした。
それをどうするか、は研究者ではなく技術者の仕事、でしょうか？Googlig では産総研さんからのリスクマップ作製に関する委託業務がいくつか引っかかりますので、初めからLSMでなくリスクマップを作成する目的だったのか、単に予算の問題だったのかもしれません。
九州地域における斜面災害リスク評価主題図改良業務：SVMを利用
九州北部地域の地質情報解析によるリスクマップ(案)の作成業務：災害特性に着目

国内では崩壊データを整備する仕組みがありません（国が災害データを収集する仕組みはあります）。災害データを使うと、民家に近いところが危ない、危険区域や危険個所に指定されている範囲が危ない、という結果になります。これは当然で、あちこち崩壊しても被災した箇所しか国には報告されない、報告外の箇所は未崩壊として誤って扱われてしまいやすいというのが理由です。そのインベントリを教師データに使用すれば、属性としての地質よりも区域等が効いてきます。機械学習モデルの中で地質の重要性が薄れてしまうのです（地すべりはそれでも地質が効くのでかなり重要なのでしょう）。

では、LSMを作成するのに災害データしか入手できない場合はどうするか？

①航空写真、衛星写真等から作成する。
産総研さんの「一部」というのはこれかもしれません。時間と費用がかかるものの、現状では最も精度の高い方法と言えるでしょう。写真等から抽出するための画像処理や機械学習が進化すれば解決すると思います。が、現状では最後に人手が必要です。もう少し待ちましょう。

②文献で利用されているデータを利用する。
昨日の文献通りです。ローカルの評価にデータが不足すれば、世界中から集めてくればよいのです。文献個々のサプリメントデータを集めるのも有効です。機械学習にかける際には多少の工夫が必要になると思いますが、データ不足には効果的です。

③LSMをあきらめてリスクマップを作成する。
この場合、LSMを作成したうえで、ある地域に台風〇〇号並みの豪雨が来たらどうなるか、その際に被災しやすい場所はどこか？といったような手順は踏めません。被災関係なく豪雨や地震で崩れそうな場所を予測することも正確にはできないのですが、被災しやすい箇所は推定できます。

産総研さんの報告からは、①と③のいずれの方法をどの程度採用しているのかは明確ではありませんでした。聞く限りでは作成されたハザードマップ？リスクマップ？よりは収集された属性データを公開していただく方が、後々利用目的に応じて柔軟なマップ作りができるように感じました。

Landslide Susceptibility Map using ML その2

2019年～2024年まで、機械学習を利用したLSM作成に関する文献を20編程度集めました。

• 文献で使われている手法は、RF や SVM からスタートし、DNN が現れて、2022年ごろから XGB や LightGBM が加わっていました。
• 2022年までは ML の結果のみで LSM を作成していますが、その後は何かと組み合わせるパターンが加わっています。
• imbalance data への対応は under samplimg が主体。1:1～1:2が良い結果を導いているようです。
• 地形スケールの影響を検討している文献もちらほら。30m程度が多いように感じるものの、何が良いかはケースによる、といった文献がありました。
• ブロック全体からポイントを作るか、中心から作るか。滑落崖からつくるか移動体から作るか。微妙な差でしたが滑落崖全体から作るのが良かったようです。

気になった文献のメモ( ..)φ
データ構成のフローチャート Fig. 6 がわかりやすいし、全体構成が参考になります。

Important considerations in machine learning-based landslide susceptibility assessment under future climate conditions | Acta Geotechnica

• undersampling (one-sided selection, k-means clustering, gridded hyperspace even sampling, random sampling)の比較。グリッド型ハイパースペースサンプリング手法が有効
• XGBoostモデルの外挿予測の確認。

 Even though the case study presented in this work is implemented for the state of California, using only California-based datasets for model training may pose several challenges. For instance, the number of positive data samples may be insufficient, leading to a severely imbalanced dataset. In addition, given the changing climate conditions, developing a robust model for future landslide susceptibility assessment requires a diverse dataset not limited to the historical California data to prevent issues with extrapolation due to the impacts of climate change. To address these challenges, a global dataset is used here for training to ensure adequate data samples and reduce extrapolation errors. 

Research indicates that tree-based models, particularly XGBoost [24], are robust to correlated features, and variables do not necessarily need to be removed due to high　correlation. However, to reduce model dimensionality, only one out of each two highly correlated variables is kept in the model.

機械学習をかじっていると転移学習が有効であることを知っているので、この文献の似た考え方は受け入れやすいでしょう。データを利用させてもらえると、日本の LSM は作成しやすくなります。

機械学習による水位予測 その2

 ”water level machine learning” + Googlimg で多数引っかかる内容です。

頭から読んでみたところ、地すべりのSusceptibility Map の作成の文献とは傾向が異なっており、アルゴリズムの比較が多く目につきました。データセットによって予測性能の良いアルゴリズムが異なる点はよく知られていますが、見事にバラバラ。それらのアンサンブルが少ないのは機械学習分野の流れと異なりますが、浸透するにはもう少し年月が必要なのでしょう。

内容的には玉石混合ですが、アルゴリズムの比較という点では有用ですので、引っかかったものから書き残しておきます。

Exploring machine learning algorithms for accurate water level forecasting in Muda river, Malaysia - PubMed

• DNN, LSTM, XGBoost
• １日前の水位を利用する場合はDNNがBEST。
• ７日間の予想ではLSTMが良好。
• データの量と質にに大きく依存する。

Deep Machine Learning-Based Water Level Prediction Model for Colombo Flood Detention Area

• DNN, LSTM
• Daily rainfall, Daily evaporation, minimum daily temperature, maximum daily temperature, daily relative humidity at daytime/nighttime, and daily average wind speed
• LSTMの方が良好

Prediction of Water Level Using Machine Learning and Deep Learning Techniques | Iranian Journal of Science and Technology, Transactions of Civil Engineering

• Random Forest, XGBoost, RNN, BiLSTM, CONV1D-BiLSTM
• XGBoostがBEST

Water level prediction using various machine learning algorithms: a case study of Durian Tunggal river, Malaysia

• 線形回帰 (LR)、相互作用回帰 (IR)、ロバスト回帰 (RR)、ステップワイズ回帰 (SR)、サポートベクター回帰 (SVR)、ブーストツリーアンサンブル回帰 (BOOSTER)、バッグドツリーアンサンブル回帰 (BAGER)、XGBoost、ツリー回帰 (TR)、ガウス過程回帰 (GPR）
• 29年間の日降水量を使用。
• 自己相関関数によるラグタイムを考慮した４つの入力データ（シナリオ）で試行
• 過学習対応のため、Lossをearly stoppingに使用
• 異なるフォールド（3、5、7、9）を使用
• GPRがBEST
• 精度を高めるには長いラグデータが必要
• 過学習の問題は通常、小さなデータセットを使用してモデルを開発するときに発生するが、1990年から2019年までの日次データはそのような問題を回避するのに許容できる長さであった
• 不確実性分析を実施

Modeling the fluctuations of groundwater level by employing ensemble deep learning techniques

• DL, アンサンブルDL
• 4 つの井戸の地下水位から残りの井戸の地下水位の予測
• 最大 20 日前の地下水位記録の時系列を入力として5つの井戸の翌日の地下水位を予測
• 合計レコード数 276、2017 年 10 月 20 日から 2018 年 5 月 1 日までトレーニングセット (70%、194 レコード) 、2018 年 5 月 2 日から 2018 年 7 月 22 日までテストセット(30%、82 レコード)
• EDLの方が良好

機械学習による水位予測

2016年に DNN を利用した水位予測を試行していました。
https://phreeqc.blogspot.com/2016/11/using-deep-learning-4.html

それから8年が経過し、ようやく地下水チームからも機械学習による予測に興味を持つ方が現れ始めました。
8年の間に、"その3"で書き残していたように、いくつかのラグを有する特徴量を利用しないとうまく機能しないことがわかっています。また、決定木系の回帰では、渇水年など外挿となる予測が困難なことも体感しました（コチラの文献を見ると理解しやすいと思います）。

国内では Signate で自治体主催の水位予測コンペが開催されたり、関連文献も見かけるようになりました。が、国外のほうが圧倒的に進んでいると感じています。国内では物理ベースのモデル構築が好まれて、海外よりもより保守的だったのかもしれません。物理モデルの構築にはいくつかの仮定が必要であり、問題を包括的に理解しシンプルに表現する力量が必要です（昔の人は偉かった）。この力量が問われてきた世界では、仮定が必要なく、説明性も低い機械学習を学問として受け入れるのに抵抗があったのかもしれません。が、大量のデータを容易に得られる現代では、そこから新たな知見を発見するというデータ駆動型のアプローチも重要なんですよね。

関連する内容が、以下に書かれています。
Full article: Water level prediction using various machine learning algorithms: a case study of Durian Tunggal river, Malaysia

For example, in Malaysia, a physical-based model was developed to assess one river's floodplain and water level (Mohamad et al., Citation2014). In order to build such a physical model, there was the need to collect massive data and information on top of the cost involved in building such a physical model. Over time, numerical models overcame the limitations of the physics-based models. For instance, a numerical model was developed by (Wu et al., Citation2014) to forecast water levels at the Yangtze River. However, a study conducted by (Guan et al., Citation2013) reported that many errors were encountered during the development of a numerical model. Therefore, despite the noticeable improvements in the numerical models, they still have limitations, such as the need to mimic some of the physical phenomena to improve their accuracy and reliability.

Recently, data-driven techniques were shown to have overcome traditional models’ drawbacks and proved to be more accurate in modeling complex engineering problems. 

ま、データがたくさん必要なのは機械学習も同じです（いえ、それ以上）。が、仮定の必要がなく、ハイパーパラメータの調整さえすればあとは自動でモデルを構築してくれる、ほぼリアルタイムで予測結果を算出できる、しかも予測性能が高くなるとなれば、使わない手はないでしょう。合理性の観点から国外で進んだのかもしれません。

続く。

水位予測に使用する特徴量

地すべり速度、水位の予測性能を Random Forest と重回帰で比較した文献です。

A comparative study of random forests and multiple linear regression in the prediction of landslide velocity | Landslides

RFの回帰の仕組みとしては Fig8 が直感的に理解しやすいと思います。この木が多数あって平均をとるのがRF、重みを含めているのが勾配ブースティングのイメージです。線形補間と異なり階段状の予測モデルになるので、ある程度深い木 or 木の数が必要、それで過学習にならないためには多数のデータや特徴量が必要、ということになります。

sckit-learn の cheat sheet では、50以下でデータを集めなおし、10万以下で線形SVRなど、さらに必要なら非線形やアンサンブル手法（RF、勾配ブースティングなど）という流れになっています。
12. Choosing the right estimator — scikit-learn 1.5.2 documentation

この文献の良いところは入力データを明示しているところです。地すべりの速度を予測する前に水位を予測するのですが、それに使用した75の特徴量をラグ付きで示されています。文献によってはどの程度のラグを使ったかが明確にされていない場合もあるので、表１のように整理して表示されるとありがたいです（理解できない内容もありますが）。

Python から NASへアクセス

Python から NAS へアクセスする必要が出てきました。

Windowsであれば問題にならないのですが、Ubuntu だとマウントが必要です。
CIFSを利用します。コマンド打ちの方が圧倒的にシンプルです。

$ sudo apt install cifs-utils
$ sudo mount -t cifs //[nas_ip]/[nas_share] /mnt/[name] -o user="user_name",pass="nas_password"

どうしても Python スクリプト内で実施する必要がある場合は以下。

import subprocess
＃マウントディレクトリを作成。
mount_point = "/mnt/[name]"
subprocess.run(["sudo", "mkdir", "-p", mount_point],
                input=[sudo_password].encode(),
                stdout=subprocess.PIPE,
                stderr=subprocess.PIPE)
#CIFS共有をマウント
mount_cmd = ["sudo", "mount", "-t", "cifs",
              f"//{nas_ip}/{nas_share}",
              mount_point,
              "-o",
              f"user={user_name}, pass={nas_password}"]
subprocess.run(mount_cmd,
               input=[sudo_password].encode(),
               stdout=subprocess.PIPE,
               stderr=subprocess.PIPE)

Blind Source Separation

ブラインド信号源分離（BBS)は、2つのソースが混在した観測記録（波形）から、元のソースを分離しようという試みです。

音声であれば同じ媒体（空気）を通すので経路中の波形の変化は無視できそうですが、地震波の場合は２つの観測点まで同じ地盤を通るわけではないので、その変化を無視できそうにありません。その場合でも、どの程度使えそうか試してみました。

・ウェーブレット変換
・独立成分分析（ICA）
・非負値行列因子分解(NMF)

ウェーブレット変換は、簡単なsin波の合成でも振幅が大きく異なるとダメでした。
ICAはOKです。ソースの1つとして同じ地震波形を利用した場合、2つ目のソースの振幅が大きく異なっても奇麗に抽出できました。が、自然界にそのような都合の良いモデル波形はなく、実際の2カ所で観測された波形を使うとダメでした。NMFも同じくダメでした。

なかなか難しいようです。

ヘテロダイン検波

DAS（ヘテロダインOTDR）を理解するため、順に可視化してみました（Thanks! GPT）。
誤っていたらいずれ気づくでしょう。その際に修正しましょう。

まずは import。

import numpy as np
from scipy.signal import hilbert, butter, filtfilt
import matplotlib.pyplot as plt
 

レーザーを作成します。

# サンプリング周波数
fs = 5000  # サンプリング周波数
t = np.linspace(0, 1, fs, endpoint=False)  # 1秒間の時間軸

# 信号の周波数
f_signal = 1000  # 1000Hzの振動信号
signal = np.sin(2 * np.pi * f_signal * t)

# 周波数シフト
carrier_freq = 1010
carrier_wave = np.sin(2 * np.pi * carrier_freq * t)

# 後方散乱光
# 光ファイバーから得られた2ショット
# 位相が少しずれている
scattered_signal1 = np.sin(2*np.pi*(carrier_freq)*t
　　　　　　　　　　　　　　　　　　　+ 2*np.pi*0.1)
scattered_signal2 = np.sin(2*np.pi*(carrier_freq)*t
　　　　　　　　　　　　　　　　　　　+ 2*np.pi*0.2)
                   
plt.plot(t[:100],signal[:100],label="signal")
plt.plot(t[:100],carrier_wave[:100],label="carrier",c="green")
plt.plot(t[:100],scattered_signal1[:100],label="scatteres1")
plt.plot(t[:100],scattered_signal2[:100],label="scatteres2")
plt.legend()
 

戻ってきた光を干渉させます。
sin波の積和の公式より、高周波成分と低周波成分が生まれます。

 # 干渉信号の生成
interfered_signal1 = scattered_signal1 * signal
interfered_signal2 = scattered_signal2 * signal
plt.plot(t[:100],interfered_signal1[:100])
plt.plot(t[:100],interfered_signal2[:100])
 

光の高周波数は検出できませんが、干渉後の低周波成分は検出可能です。

 # Butterworthフィルタ(ローパス)
cutoff = 50
# 正規化周波数
Wn = cutoff / (f_signal / 2)
# Butterworthフィルタの設計
b, a = butter(4, Wn, 'low')
filtered_signal1 = filtfilt(b, a, interfered_signal1)
filtered_signal2 = filtfilt(b, a, interfered_signal2)
plt.plot(t,filtered_signal1)
plt.plot(t,filtered_signal2)
 

低周波側の信号から位相を取り出し、２ショット間の位相差を算出します。
ここはプロのノウハウがあるそうです。図書にも見当たらないので推定です。

 # ヒルベルト変換を用いた瞬時位相の抽出
analytic_signal1 = hilbert(filtered_signal1)
instantaneous_phase1 = np.unwrap(np.angle(analytic_signal1))
analytic_signal2 = hilbert(filtered_signal2)
instantaneous_phase2 = np.unwrap(np.angle(analytic_signal2))
# 瞬時周波数の計算
instantaneous_frequency1 = np.diff(instantaneous_phase1)
                                   / (2.0 * np.pi) * fs
instantaneous_frequency2 = np.diff(instantaneous_phase2)
                                   / (2.0 * np.pi) * fs

plt.figure(figsize=(10, 10))
# プロット：位相
plt.subplot(311)
plt.plot(t, instantaneous_phase1)
plt.plot(t, instantaneous_phase2)
plt.title('Extracted Phase from Interfered Signal')
plt.xlabel('Time [s]')
plt.ylabel('Phase [radians]')
plt.grid(True)

plt.subplot(312)
plt.plot(t, instantaneous_phase2-instantaneous_phase1)
plt.title('Extracted Phase from Interfered Signal')
plt.xlabel('Time [s]')
plt.ylabel('Phase [radians]')
plt.grid(True)

# プロット：周波数
plt.subplot(313)
# t変数の長さをinstantaneous_frequencyに合わせる
plt.plot(t[:-1], instantaneous_frequency1)
plt.plot(t[:-1], instantaneous_frequency2)
plt.plot(t[:-1], instantaneous_frequency2
                  -instantaneous_frequency1)
plt.title('Instantaneous Frequency')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.ylim(-20,20)
plt.grid(True)
plt.tight_layout()
 

Complex Hilbert PCA

DASでの震源決定に複素主成分分析を利用されている事例を拝聴しました。
確かに、密な観測点ですので相互相関は得られそうですし、ノイズも落としやすいと思います。

複素主成分分析についてはよく知らなかったので調べてみました。
時系列データに対し、ヒルベルト変換をかませてからPCAにかけるだけのようですね。手順としては以下の通りです。 (by GPTさん)

1．時間領域の波形読み込み

2. 信号の標準化
標準化（mean = 0, standard deviation = 1）することで相互相関の計算が比較的均一になります。
signals_mod = [(s - np.mean(s)) / np.std(s) for s in signals]

3. ヒルベルト変換
ヒルベルト変換を使って各信号の解析信号（analytic signals）を計算。ヒルベルト変換は信号を複素数表現に変換し、瞬時位相や包絡線を得るために使用されます。
analytic_signals = [scipy.signal.hilbert(signal) for signal in signals_mod]

4.エルミート共役行列の計算
解析信号のエルミート共役（Hermitian transpose）を計算。エルミート共役は、行列の転置を取ってから各要素の複素共役を取る操作です。
analytic_signals_h =np.conjugate(analytic_signals.T)

5. 相互相関行列の計算
解析信号の相互相関行列を計算。複素数表現を使用することで、信号の振幅変化だけでなく、位相変化も考慮した相互相関を計算することができます。
corr_matrix = np.dot(analytic_signals, analytic_signals_h) * (1/ｆs)
1/fs倍することで正規化しています

6.主成分分析
相互相関行列の固有値（eigenvalues）と固有ベクトル（eigenvectors）を計算。
重要な情報を保持しつつデータの次元を削減します。小さな固有値に対応する成分は、ノイズである可能性が高く、これらの成分を無視することでデータのノイズ除去が可能です。
eigenvalues, eigenvectors = np.linalg.eigh(corr_matrix)
principal_components = np.dot(eigenvectors.T, analytic_signals)

ここから主成分をいくつか抜き出して、ノイズを落とした信号を再構築すれば、前処理は終わりです。再度、相互相関でΔtを推定するなり、位相差を計算するなりして利用するのでしょう。

***********************************************
20241007追記

複素行列Aのエルミート転置イロイロ。
A.H
np.swapaxes(np.conjugate(A),1,0)
np.conjugate(A).T
np.conjugate(A.T)

.Hが簡単です。

Landslide Susceptibility Map using ML

2020年頃から機械学習が Landslide Susceptibility Map の作成に用いられるようになりました。機械学習は安定計算と同じような基本ツールの一つになっています。日本の土木分野は、かなり出遅れています。

初期の文献では、ある事象を単一モデルで学習し、予測マッピングするだけでした。その後、インバランスデータの取り扱いやアルゴリズムの選択などの文献が現れています。参考になるのは以下。1：2がBESTだった例です。
Comparative study on landslide susceptibility mapping based on unbalanced sample ratio | Scientific Reports (nature.com)

今日では既知の地震に対してのみならず、新たな地震に対しても利用できるような汎用性を持たせる作り方が報告されています。新たな地震に対する危険度や損害の程度、それがどこで発生するかを予測したいからでしょう。過去の地震をトレース、予測するだけだと予防や対策につながりませんので。

また、MLと別の分析結果を組み合わせる例が増えているようです。予測された危険度毎に建築物の数を集計するなど、家屋への被害の程度を簡易的に表現する単純な方法でも現段階では有用だと思われます。（建築物の数を考慮すると Risk Map になります。）
Enhancing co-seismic landslide susceptibility, building exposure, and risk analysis through machine learning | Scientific Reports (nature.com)

MLとKDEとのアダマール積、InSARの結果などを組み合わせて評価している例もあります。残念ながら、これらは汎用性に劣るため流行らないでしょう。先のReg3Dは後で組み合わせるより先に計算しておいてMLに投入する方が良いでしょう。どの範囲、というのはDEMだけだと難しいですから、確率をメッシュに与えた方が当たりそうな気がします。

MLと何を組み合わせるか、アイデアはたくさんあると思います。基本ツールの結果を組み合わせ、本来解くべき問題に早々に挑めると、少しは前進するのではないでしょうか。

Landslide Susceptibility Map using 3D LEM

RegionGrow3D: A Deterministic Analysis for Characterizing Discrete Three‐Dimensional Landslide Source Areas on a Regional Scale - Mathews - 2024 - Journal of Geophysical Research: Earth Surface - Wiley Online Library

In this study, the RegionGrow3D (RG3D) model is developed to broadly simulate the area, volume, and location of landslides on a regional scale (≥1,000 km2) using 3D, limit-equilibrium (LE)-based slope stability modeling. Furthermore, RG3D is incorporated into a susceptibility framework that quantifies landsliding uncertainty using a distribution of soil shear strengths and their associated probabilities, back-calculated from inventoried landslides using 3D LE-based landslide forensics.

動画もあります。
Landslide susceptibility map generated using RegionGrow3D (youtube.com) 

インベントリがあれば、そこから逆算してΦを出す、その確率密度の分布から累積分布を書いて susceptibility とみなす、その確率でどの大きさの地すべりがどこに分布するかを推定する。このような土質力学ベースの計算が可能な Matlab ツールをUSGSが提供しています。

文献では、すべり層厚を決めるのに隆起量などをパラメータとした別のシミュレーションを利用していますが、この辺りはインベントリの平均層厚でも構わないでしょう。

すべる・すべらないのバイナリを拡張し確率を持たせた点、thresholdに応じて地すべりの大きさを扱える点が魅力です。日本は完全に置いてけぼりです。

Susceptibility, Hazard, Risk

4～5年前から、海外の文献で “Landslide Susceptibility Map” をよく見かけるようになりました。

”感受性マップ”と直訳するとわかりにくいのですが、危険性や事象が起こりそうな確率を表す地図を意味しています。日本で使われるハザードマップが近いと思います。

が、Susceptibility と Hazard は使い方が異なります。前者は空間に対し、後者は時空間に適用されます。さらに、Hazard が発生した場合の影響度を考慮するとRiskになります（KY活動からリスクアセスメントに移った頃を思い出せば、容易に理解できると思います）。

このあたり、検索してみるとすぐに出てきます。例えば以下。

Landslide Inventory, Susceptibility, Hazard and Risk Mapping | IntechOpen

DOI: 10.5772/intechopen.100504

2. Definition and concepts

・Susceptibility refers to the probability of occurrence of an event within a selected type during a given location.

・Hazard refers to the probability of occurrence of an event within a selected type and magnitude during a given location within a reference period.

・Risk refers to the expected losses or damage by events during a given region, which are the products of susceptibility, hazard, and elements in peril.

・Vulnerability means the degree of loss to a given element of the set of elements in peril resulting from the occurrence of natural phenomena of a given magnitude. It is expressed on a scale from 0 (no damage) to 1 (total damage). 

・Elements at risk is potentially vulnerable of properties, population, and economic activities including public services in peril during a given area.

用語を正しく用いないと伝わりません。忘れないようにしましょう。

hDVS

An Introduction to Distributed Optical Fibre Sensors

先の洋書よりもハードウェアや原理について詳しく書かれています。先日使用した hDVS について、大まかな流れは理解できました（細かい点はわからないことばかりで先は長い）。

以下、流れの抜粋です。( ..)φ

狭帯域レーザーの出力は、Optical local oscilator （OLO）パスとパルス形成パスに分割される。

前者は角周波数 ω0 で放射され、ELo・exp(jω0t) の形の時間依存電界を持つ。

後者はAcoustic optical module （AOM） によってω0 + Δω のプローブパルスに周波数シフトされる。

サーキュレータ等を使用して、プローブパルスをセンシングファイバーに発射し、後方散乱光を OLO と混合する別の結合器を介してバランス型受信機に導く。

検出器に降りかかる電界は次のように表される。

Eωt = EL0・exp(jω0t)+Eb(t)・exp{j[(ω0 + Δω)t + Φ(t)]},  (3.29)

光検出器は表面に到達する光パワーに応じて電流を供給するので、光電流は次のように表される。

iphot = Rd{EL0^2+Eb^2+2EL0Eb・exp[j(Δωt+Φ(t))]}  (3.30）

バンドパスフィルタを使用して中間周波数IF=Δω のみを選択することで、OLO の電界と後方散乱信号の積であるヘテロダイン項ihet(つまり、式 3.30 の 3 番目の項) が得られる。これは、受信器の出力を支配する唯一の光学的寄与である。

ihet = Rd{2EL0Eb・exp[j(Δωt+Φ(t))]}  (3.31）

この信号を二乗することにより、後方散乱信号と直接比較できる波形が得られる。つまり、自由空間の各ポイントから返される後方散乱パワーに比例する。

ヘテロダイン OTDR では、バンドパスフィルタ適用時に後方散乱信号の振幅が検出され、位相情報が破棄される。対照的に hDVS では、後方散乱の位相は、位相情報を保持しながらも電子機器で処理できる周波数であるダウンシフト IF を介して測定される。(プローブパルスの周波数が50MHzシフトアップされている場合、ヘテロダイン後方散乱信号の周期は、約2mのファイバーに相当する。）

処理の手順
1．ファイバーに沿った距離として位相を検出。
2．レーザーのショット間の位相差を算出。
3．選択した距離間隔（ゲージ長）にて位相差を空間微分。
※ゲージ長全体にわたる位相の変化をレーザーショット数の関数として記録することも可能。これにより、ひずみの時間微分（ひずみ速度）が得られる。

これらの処理後、ファイバーに沿った距離と遅い時間（つまり、プローブパルスの連続に見合った時間スケール）の関数として、アンラップされた位相の画像が利用可能になる。

Distributed Acoustic Sensing in Geophysics

Distributed Acoustic Sensing in Geophysics: Methods and Applications 

昨年購入し、何度か読むにつれ徐々に頭に入ってきました。基本原理についても最初に記載があります。

１章

<A(z)> = 1/(A0·Fs)·τ(z)⊗[v(z+L/2)−v(z−L/2)] = 1/(A0·Fs)·τ(z)⊗ε˙(z)·L
A0 = λ/(4π·neff·Kε) = 115nm

The elongation corresponding to ΔΦ = 1rad is A0 = 115nm, calculated for λ = 1550, neff = 1.468  and Kε = 0.73, which has been measured for conventional fiber. 

The DAS signal is a convolution of pulse shape  with a measured field which is the spatial difference in fiber elongation speed of points separated by a gauge length.

2章

Φ（z,t） = atan(Q/I)
ε（z,t） = Φ(z,t)·λ/(4π·neff·Kε)/L

Gauge　Length
・ゲージ大 → 帯域小, パルス大 → SNR大
・ゲージ小 → 地震信号の忠実度は高まるが、SNRが低下する。
・必要なデータ帯域幅のゲージ長を確認する必要がある。
・良好な信号を得るために、最適なゲージ長、パルス幅を選択する。

Fading
・Phase unwrapping algorithm
・フェーディングは DAS 取得の自然な特徴。常にハードウェアとソフトウェアの両方で対処する必要あり。
・複数の光周波数を使用する IU は、本質的にフェーディングの低減に優れている。
・取得後の処理でフェーディングに対処可。

Common-Mode Noise
・Stack
・周辺の環境音によって発生。
・簡単な取得後の処理で、そのほとんどを除去可能。

Spacial Calibration of Channels
・Tap test
・チャネルの位置を決定することは、結果として得られる精度にとって重要。
・既知の位置をキャリブレーションポイントとして使用し、インターリーブチャネルをそこから補間する。

その他
・可能であれば、マルチモード ファイバーではなくシングルモードファイバーが望ましい。
・接着されたファイバーケーブルを使用すると、より優れたSNRデータが得られる。回収可能なファイバー配置方法を使用して目的に適したデータを取得することも可能。
・光ファイバー接続を清潔に保つ必要あり。
・DAS データでタイミング情報が確実に保持されるよう、GPS への接続が必要。
・ファイバーの入射角応答がジオフォンの応答と異なる場合は、まずその影響を判断し、次にジオフォンのような応答を除去する。
・DASのSNRはジオフォンからのデータよりも本質的に低くなる。しかし、高密度で取得可能であるため大幅に改善可。

DAS その2

1年ほど前にこんなことを書いていたのですが、先週、DAS (Distributed Acoustic Sensing) に触れる機会がありました。

と言っても、私は作業員で路上を叩くだけ（Tap Test）。後輩君がメインで進めるお手伝いした形です。

まだ機器の構造と理論を理解できていないのですが、理解したいと思わせる結果を見せてくれます。橋の揺れ、ダンプの走行など数㎞の揺れを表示してくれるのが面白い。測定画面を見ているだけでも時間が過ぎるくらい、久々に興味を惹かれる機器でした。

今年の物理探査学会講演会でも１セッション設けられるほど注目されています。早く理解しましょう。以下、文献より備忘録です。

Seismic Monitoring With Distributed Acoustic Sensing From the Near-Surface to the Deep Oceans

φ = k·n·2z
Δφ = 4π/λ·Δn·Δz + 4π/λ·n·ε·Δz = 4π/λ·Δneff·Δz
Δneff = Δn + n·ε
ε˙(z) = [v(z+L/2)-v(z-L/2)]/L

where k is the wavenumber (k = 2π/λ, with λ　the light central wavelength), n is the　refractive index and 2z is the two-way length　that the backscattered light travels. 

２段階認証

 google の2段階認証アプリに異常が発生し、他のアプリに切り替えました。

その間、約１週間。ひとまず会社とGoogleの２アカウントは復帰できましたが、それ以外は不具合の発生しているアプリの中で眠っています。ここにもようやくたどり着きました。

GitHub がまだ眠っています。ひとまずローカルからアップはできるのでしばらく待つつもりです。

怖いですね。

core MP135

core MP135 を購入。
Python コードを動かすまでの手順です。

balena-eacharを管理者モードで起動
ｍicroSDにイメージを焼く

USBケーブル接続後、シリアル通信
ID: root
Pass: root
$ passwd root

ユーザー追加
$ adduser "user name"
sudo group に追加
$ gpasswd -a "user name" sudo

root パーティション拡張
$ cd /usr/local/m5stack
$ bash resize_mmc.sh 

LANに接続後、SSH接続
IP確認
$ ifconfig

・自動起動させるプログラムをWinSCPでコピー。
・python3 で稼働確認。
・vimでrc.local を変更し、自動起動設定。

double-layer two-phase formulation その3

double-layer two-phase SPHの実装が概ね完了しました。

実装を通し、理論と実際を理解できました。
現状、2点の問題があるようです。

１つ目は seepage force (drag force) の計算式。

これまでの文献にも、DualSPHysicsでも同じ式が使われています。が、透水係数の大きな場合にしか実用的ではないようです。1e-4㎝/sオーダーだと、γ/kが大きくなりすぎて、少しの速度差が大きな拘束力を生みます。土塊がすべろうとしても、地下水がその動きに追随するまで拘束されるような挙動も見受けられました。

２つ目は密度変化。

土の間隙内の水が地表に浸出した場合、水相の密度は大きく変化しますが、この表現が困難。いえ、pressureの式で参照している基準密度を変更すれば容易に反映できるのですが、粒子が重なり合い、挙動がやや不安定でした。初期配置も面倒で、実用的ではないですね。水収支はともかく、土の動きだけを正しく反映したい、ということであれば計算の軽い single-layer two-phase formulation でよいと思います。

土相側でも層毎に密度を変えている場合、その境界で密度が訛るのは現状では仕方がないと考えています。が、この密度変化の問題は将来的に何とか解決されそうな気がします。たとえば、土層別にフラグを付けて、個別に計算するとか。それが実現象を正しく反映するようになれば、より実務に浸透していくのでしょう。

今年前半より取り組んできた SPH コードの修正と実装は、これで一旦終了です。
結果的には単相での計算はOK、2相は「定性的評価」まで、土木分野の実務で求められる諸問題には適用が難しい場合が多いと判断されます。が、もう少しだと思います。

将来に期待しましょう。

double-layer two-phase formulation その2

振り返ると、double-layer での実装に取り掛かってから１か月以上かかっています。

以下を参考に始めたのですが、試行錯誤の結果、最終的にはPersianSPHの文献と同じような形に落ち着きました。
A coupled fluid-solid SPH approach to modelling flow through deformable porous media - ScienceDirect

SPH approach for simulating hydro-mechanical processes with large deformations and variable permeabilities | Acta Geotechnica (springer.com)

２次元、３次元共に土、水の相互作用までは動くのですが、地下水が地上に出た場合の密度変化をうまく表現できません。相互作用の検証も終わっていませんので、まだ時間がかかりそうです。

SPHでの２層2相表現は定性的にとどまると伺ったことがあります。それほど新しいものでもないので、国内でも実装されている方はいらっしゃいます。その先生がおっしゃるのですから、そうなのかもしれません。

まだ１か月。もう少し追いかけてみましょう。

システム修復

片道４時間、日帰りで、観測システムを修復してきました。

Linuxに不慣れだった後輩君が設定していたため細部の設定に漏れがあったこと、私の管理ミスでこれを見逃していたこと、夏場の高温でハードの一部に不具合が発生していたこと、最終的には停電がトリガーとなったこと等、いくつかの要因が重なり問題となりました。で、現場まで赴くことに。

原因を突き止めてしまえばその場で対応できる、比較的簡単なソフト側の対策で済んだのですが、熱だけはハード対策が必要そうでした。が、これは高価な製品を購入しないと対応できません。うーん。

システム修復後、ついでに DNS を設定。GUIからなぜか反映されていなかったので、コマンドで変更。

$ sudo vi /etc/systemd/resolved.conf

[Resolve]
DNS=8.8.8.8

$ sudo systemctl restart systemd-resolved.service

apt update は成功しますが、upgrade で失敗します。古いリストが残ったままなのでしょう。リストを全消ししてから再度 update & upgrade。upgrade には時間がかかりました。

$ sudo rm -r /var/lib/apt/lists/*

ひとまず周辺環境を含め目についた問題は全てつぶしておきました。今回は停電により、隠れていた問題も同時に顕在化しました。最初に見逃さないようにしないと、このように不要な対応が発生します。気を付けましょう。

MEMS Rain Gauge

IOTデバイスを用いて、振動から雨量を推定するプログラムを公開しました。https://github.com/T40O0/ADXL355_SPI_M5_SD_RainGauge/tree/main

約１年間観測して、改良を施し、安定して測定できるようになりました。

ポイントはデジタルフィルターを利用しないことでした。雨滴による衝撃を観測しているので、高周波まで扱えるほうが精度よく雨量を推定できることがわかりました。周波数領域は扱えませんが、ただただ時間軸で衝撃をカウントするという単純な方法が推定には良いようです。

転倒マスとは異なり、構造的に火山灰等で詰まることがなく、メンテナンスの楽な点が長所でしょう。火山灰を被っても、雨で洗い流されます。

安価な点も魅力です。個人で作るなら材料費2万円弱、従来の1/5～1/10の価格ですので、多点に設置できます。近い場所でも雨量は異なりますので、各地すべりブロックに1台、などという計画も可能です。
さらに、1滴の振動から観測できますから、時間方向の解像度が転倒マス型に比べ格段に向上します。

欠点は消費電力。電池のみで数か月持てば理想的なのですが、ソーラー＋バッテリーか商用電源が必要です。モバイルバッテリーでも１週間しか持ちません。

また、振動のノイズに弱いことも短所でしょう。雨が降っていなくても、加速度を観測すると微小な雨量としてカウントします。設置場所毎にノイズを除去するための閾値の設定が必要です。

まだまだ改良の余地がありますが、ひとまずここで公開。良いアイデアが浮かべば追加していきましょう。

SPHの水圧計算

前回から約1か月かかりましたが、２相２層での水圧計算に目途が立ちました。

https://www.sciencedirect.com/science/article/pii/S002076831730286X

上記の78-82式を実装することで解けるようにはなったのですが、まだ土層を出た後の水の密度変化までは安定して解けません。最初の含水率のまま解くことはできるのですが、空中でそれは誤りですし、勝手に密度を上昇させることもできません。圧力が下がって自然に密度が上昇する、というように計算させるとすぐに発散します。そのあたり、この論文では明記されていません。土層の計算だけ追うのであれば、１層での計算が良いと思いますが、まだ何か見落としているようです。

ひとまず、次の drag force を安定して解けるか見てみましょう。

double-layer two-phase formulation

SPHにおいて、two-phase を single-layer で実装するか、double-layerにするか？

前者は計算が速い！
でも、水の流れは表現できません。（工夫すればできるのかもしれませんが、single ではなくなります。）

後者は計算が遅い！しかも圧倒的に。
でも、水の流れは表現できます。

考えた結果、もともと計算してみたかった「崩土から水の抜ける表現」が可能な後者を選択しました。既存のプログラムを容易に変更できそうな点も good です。

2相 を double-layer で解く文献はたくさんあります。が、なぜか式が統一されていません。以下が最も詳しく書かれていましたので、これを参考に実装します。
A coupled fluid-solid SPH approach to modelling flow through deformable porous media - ScienceDirect

式を整理しながら不飽和の表現を single-layer の文献からいただきました。これで、飽和・不飽和浸透と力の関係を解けるはずです。

コツコツ、組んでみましょう。

**********************************
20240707追記
地下水の上下で飽和度（水の単体）を変更すると、境界で爆発してしまいます。計算過程を追えば当然なのですが、土層境界よりも挙動が激しいようでした。飽和度を低くするには文献のように水粒子をまばらに配置するしかなさそうです。
似た現象ですが、地表流がある直下で飽和度が過剰に高くなります。透水係数が実際よりも高く計算され、水が入りやすくなります。なかなか難しい。
SPHでは密度変化がないか小さい現象（弱圧縮）を扱うのが前提なので、単体の異なる土層とか飽和‐不飽和‐地表流などを扱うのが難しいのでしょうね。

single-layer two-phase formulation

SPH で流体とソリッドをカップリングさせる方法は比較的古くから提案されています。

以前、PersianSPHを触っていた際に文献でおおむね理解していたのですが、実装目的であらためて詳細を確認しました。

基本は土質力学レベルです。それを丁寧に組み合わせると two-phase formulation の出来上がりです。SPH の場合はこれを2レイヤーで計算します。ここで倍の計算量。さらに相互作用を反映するため、土粒子の位置での水粒子の動き、水粒子位置での土粒子の動きを計算します。単純に4倍ではありませんが、計算量は飛躍的に増えます。

これを解決した、という文献がコチラ⇓
Two-phase fully-coupled smoothed particle hydrodynamics (SPH) model for unsaturated soils and its application to rainfall-induced slope collapse - ScienceDirect
A general smoothed particle hydrodynamics (SPH) formulation for coupled liquid flow and solid deformation in porous media - ScienceDirect
1レイヤーで2相の計算が可能に工夫されています。不飽和浸透も取り扱えるので一通りの計算ができそうです。が、境界条件が面倒。いえ、実装されたものを使用するのであればDtransu を利用するのと変わらないでしょうが、汎用性を持たせて実装しようとすると気が遠くなります。

著者には DualSPHysics の開発者が含まれています。GPU対応で開発されているのでしょう。実装されるまで待ちましょうか。

火山地帯の観測機器

先週は火山につけた観測機器のメンテナンスでした。

いくつか新しい機器を設置し、古い機器を入れ替えて帰社したのですが、あらゆるモノが錆びていました。火山ガスと雨でしょうね。防食、大事です。

PCは火山灰を吸い込み異音が発生。分解、清掃しても、治るかどうか。寿命を短くしたのは間違いなさそうです。

こういった環境では、安価な機器を使い捨て＆ばらまきにするか、お金をかけて防食し長期運用＆集中させるかのどちらかでしょう。費用が１～2桁違いますので、重要度に応じて使い分けるのが良いのでしょう。

SPHのベンチマーク

組み終わったSPHで静水圧のチェックは終わりました。残る弾塑性のチェックをどうしようか、と調べてみましたが、適当なベンチマークはなさそうでした。

流体のベンチマークは多く見かけます。が、弾塑性のものは少ないようです。FEMと同様に、室内３軸圧縮試験であったり、支持力公式との比較をいくつかの文献で見かけますが、ベンチマークとして共有されるには至っていないようです。

３軸は少し面倒ですので、支持力を試してみました。
が、文献に多くあるように、解析値のほうがやや大きく出る結果になります。うーん。ま、FEMと同程度に大きな値になっているので、あっているともいえるのですが。

新しい機能も組み込みつつ、もう少し探してみましょう。

EOS

equation of state

WCSPH で圧力を計算する際に用いられる式です。Monahan らが提案した Tait 型だそうですが、これ、簡単に見えて扱いが難しい。

係数を決めるのに最大流速を考慮しないといけないのですが、流速は計算しないとわかりません。ま、感覚で大体この程度、圧縮率は〇%程度をネラッテ音速ハコノクライ、みたいな決め方をします。DualSPHysics に使われている multi-phase （liquid-sediment）の文献では、音速はドメイン内の最大流速の10倍以上が利用されています。これは、マッハ数が0.1なので、概ね0.5%までの圧縮を許容するということでしょう。
https://research.manchester.ac.uk/en/publications/modelling-multi-phase-liquid-sediment-scour-and-resuspension-indu

では、個体が弾性から塑性になり、さらに非ニュートン流体として流動するような現象を追う場合はどうするか？

弾性、塑性では変形係数とポアソン比から音速を求められますので、それぞれ別の値を利用するという方法が考えられます。流体は水のように扱うことも可能です。
が、それぞれ別の音速を使うと計算が破綻しやすくなります。試算では塑性の計算で圧力が低くなり、塑性領域が広くなって流動化しやすくなるような結果を得ました。あと、2Dでは動くが3Dだと進まなくなるとか（2Dの計算って、問題やバグが内在しても露見しにくいことを、最近強く感じています）。

結論としては、弾性として求めた音速をそのまま利用するか、徐々に切り替えるか、なのでしょう。文献には出てこないのでノウハウの部類なのでしょうが、皆さんどうされているのか知りたいところです。

弾塑性SPH その2

新たに頂いたコードを参考に、以下の組み合わせで計算ができるように改造しました。

・弾性
・弾塑性
・弾性‐流体
・弾塑性‐流体
・流体

これらのほとんどは日本語の図書でも紹介されています。個体と流体を同じように実装できるのはSPHならではでしょう。

弾塑性（一部流体）の挙動はコチラ。端から土圧により崩壊し始めます。

少しc・φをあげると、崩れにくくなります。その場合、まず角が楔のように崩れ、次に長辺の2/3が崩れ、最後に残りが崩れるような時間差が見られました。ちなみに、前回のアップしているイメージは「弾性‐流体」の計算結果で、全体がペシャっと潰れるような挙動を示しています。これだけ自由な表現ができると、多様な自然現象を再現することが可能でしょう（予測は相変わらず困難ですが）。

まだ、いくつか疑問点が残っています。資料を探して早く完成させましょう。

壁粒子の圧力

SPHで壁際の粒子が壁に張り付く現象を、後回しにしていました。
ある程度動くようになったので取り組んでみることに。

まずは mDBC で様子見。
効果はありますが、残念ながら完全ではありません。いくつかの粒子は壁に残り、計算時間も余計にかかりました。

負圧が効いているのは想像がつくのですが、他の修正方法を思い付きません。
そこで googling。
すると、「まさにコレ！」といった回答を見つけました。条件も全く同じです。
Trying my chances here; anyone knows why particles stick to the wall in my SPH simulation? - Specific Domains / Modelling & Simulations - Julia Programming Language (julialang.org)

Here the simplest fix is simply to put the boundary particle density to reference density, if it ever gets below that reference value.
In this way a wall can never produce suction, but only “push”.

当たり前すぎて一般的には解説されていないのでしょうね。初心者には思いつかないノウハウです。
密度でなく圧力を変更。数行の追加で計算時間もかからず、簡単に解決できました。

弾塑性 SPH

2020年の秋から始めた弾塑性SPHコードの修正が完了しました。

基本コードを頂いてから寝かせていたため実質1年程度ですが、SPHを知ってからは14年目です。長い。
何とかSPHの離散化や各種補正方法、境界条件について基本的なしくみを理解したと言える証になるでしょう。いや、よかった。 

いくつかのハマり処があったのですが、式の実装時の静水圧Pと静水圧応力σの符号の整理が効いたと思われます。多くの文献を引用する際に、これらの書き方の流儀を統一してから実装に落とし込む必要があったのでしょう。

粒子が凝集する現象を hourglass,  zero energy mode だと勘違いしたり、他にも似たような勘違いで多くの寄り道をしましたが、振り返るとそれらも良い勉強になったと思います。
https://phreeqc.blogspot.com/2022/05/blog-post.html

早くVerification を済ませたいですね。ベンチマークは何が良いのでしょうか？
それが済んだら高速化。倍精度を多く使っていますので、GPUではなくスパコン利用（MPI並列）のほうが現実的でしょうね。

まだまだ先は長い。

相互相関関数の計算

西田究のホームページ (u-tokyo.ac.jp)
解説 (u-tokyo.ac.jp)

地震動に関する様々なシミュレーション結果をお手軽に見ることができる、素晴らしいサイトです。特に、地震波干渉法で相互相関関数が出来上がっていく過程を見たことがなかったので、非常に興味深く見ていました。

　　初期

数十分後

相互相関関数を求める際は、時間領域（直接法）か周波数領域（FFT）を利用します。配列の要素数が多ければ周波数領域での計算が早いでしょう。Python の scipy.signal.correlate では手法を指定できますし、指定しなくとも早い方を自動で選択してくれます（Ver.1.13.0）。 numpy.correlate は直接法のみです（Ver.1.26）。

以下、計算例です。設定次第なのですが、以下だと答えは①≠②=③
①直接法

sig1 = np.std(data1)
sig2 = np.std(data2)
data1 = data1-data1.mean()
data2 = data2-data2.mean()
cross_correlation = np.correlate(data1, data2, mode='same'
                                ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(data1)
f2 = np.fft.fft(data2)
cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum)
                            ) /sig1/sig2/N

③FFT2

f1 = np.fft.rfft(data1)
f2 = np.fft.rfft(data2)
cross_spectrum2 = f1 * np.conj(f2)
cross_correlation = np.fft.irfft(cross_spectrum2
                                ) /sig1/sig2/N

0をパディング + シフトで①＝②となります。

①直接法

padded_data1 = np.pad(data1, (0, N), 'constant')
padded_data2 = np.pad(data2, (0, N), 'constant')
cross_correlation = np.correlate(padded_data1,
　　　　　　　　　　　　　　　　 padded_data2,
　　　　　　　　　　　　　　　　 mode='same'
                                ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(padded_data1)
f2 = np.fft.fft(padded_data2)
cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum))
cross_correlation = np.fft.fftshift(cross_correlation
                                   ) /sig1/sig2/N