昨年度やり残した短期目標「優先度高：機械学習のスキル増強」「優先度中：DAS」は完了です。我乍ら、うまくいきました。
中期目標も順調な滑り出し。非常に優秀な後輩君が走り出し、安全管理や解析スキル等をスポンジのように吸収しました。GPTがすぐに使える環境があったこと、その精度がこの1年でかなり高くなったことも効率良く彼が育った要因の一つでしょう。時代とともに、育成方法も変化するのですね。さらに育ってもらうためには、相応の場の提供が必須です。プレッシャーです。
やり残し事項は優先度低のみ。
優先度低：流体＋個体（不連続体＋連続体）＋振動
優先度低：Dtransu の MPI/GPU 対応
優先度低：地表流＋地下水＋移流拡散
地すべりのすべり面を3次元で表現する場合、ボーリングで決定した深度と平面ブロックの位置を通すようにサーフェスを作成します。
ボーリングがない場合、あるいは地震後など多数のすべりが発生した場合に、尤もらしいすべり面を作成したいということもあるでしょう。精度が悪くとも、GIS等で広域かつ一括処理できれば便利です。
SLBLという手法を最近知ったのですが、これが手軽でしょう。
作成される深度に地質的な根拠はなく、処理も周辺4～８方向の高さを使用するだけで非常に単純なのですが、経験的におかしくない面を作成できそうです。公開されているArcGISのプラグインは不安定でしたが、うまく動いた時の結果だけを見ると、それなりに使えそうな考え方でした。
近年ではInSARと組み合わせてキャリブレーションしようとする方もいらっしゃるようです。国内では、そこまで大きなすべりは扱わないでしょうかね。
GPTでも組めそうな程度の内容です。活用してみましょう。
**************************************
20250105 追記
SLBL なら csv のみで処理できるでしょう、と自分で組んでみました。
ポリラインで地すべり範囲を囲まず、 csv 全体を作成対象とすると、接谷面図のような結果になります。山頂の尾根部にも面ができてしまうので、囲まないとダメなのでしょう。
ま、これはこれで使えそうですが。
PHREEQC を用いてバッチ試験を繰り返すことで概念モデルの妥当性を検証した内容です。水‐熱のシミュレーションはすでに実施済みのようで、残された反応部分の報告です。あえて kinetic model を使用せず、扱いやすいバッチ試験を利用されている点では、実務上の利便性が高いと言えます。
この報告、水質計算に目が行きがちですが、地質屋さんとして難しいと感じるのは岩石側、鉱物構成の設定でしょう。地下深い場所の鉱物構成をいかに尤もらしく設定できるかがポイントです。計算で水質を再現できると、設定の妥当性が高くなると考えることもできます。このトライ＆エラーは Geochemist しかできないでしょう。
固液比については、感度分析のように複数実施した結果から、尤もらしい設定を採用するだけですので難しくはないでしょう。ただ、浸透流の結果との整合性を考えるのは難しいでしょうね。バッチ試験ですから余計に。
実務では源泉への影響調査などに適用できそうです。一つのベンチマークになりそうな文献です。
PS と SBAS の自動解析スクリプトを作成したよ、という内容です。コードは公開されており、私の環境でも少しの修正で動きました。
文献にも書かれていますが、時系列変位の整理までは実施してくれません。データもローカルにダウンロードしてから計算されますので、時間がかかります。必要となるストレージサイズも大きくなるなため、個人的にはLiCSBASなどの他のコードの方が魅力的でした。
SBASベースの手法には、多くの種類があります。近年ではMSBAS（多軌道・多視点解析）による立体的な変動把握や、インフラ監視・建物傾斜などに対してP-SBASなどの高精度型の需要が強いようです。このあたりは勉強していなかったので、いずれフォローしましょう。
欧州ではデータも自動解析結果もweb上で無償公開されていますが、国内ではまだまだ。これは日本も頑張ってほしい。
一方で、自動解析ではローカル用にチューニングされていないそうですので、プロの技が必要になるとも。素人ユーザーとしては、様々な種類の自動解析結果を得られた場合でも、その精度や誤差に関して留意すべきなのでしょう。
地すべり調査に DAS を利用する方法です。
文献では平面計測のみですが、ひずみ計のかわりに途中でボーリング孔を経由する（埋め込む）ことも可能です。地すべり調査と DAS をご存じの方なら誰もが思いつくレベルの調査法なのですが、国内の実務での実施例は見たことがありません（研究で平面のみ、ボーリング孔のみは見たことがあります）。
問題は費用。千万円単位で大きくなるので発注するにも容易ではないのでしょう。精度の高い技術があることは分かっていても高価すぎる、今まで通りでも良い、といったところでしょうか。
DAS、安くなりませんかね。
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).
振動を理解するには環境計測、音波等も含め多面的に捉える必要があるというのはつくづく感じているのですが、なかなかこのような結論に結び付きません。まだまだわからないことばかりです。
Progressive Geological Modeling and Uncertainty Analysis Using Machine Learning
2段階で機械学習による層序区分、岩石区分を推定し、分布確率を同時エントロピーとして表現されています。1段階目では XYZ 座標を入力として層所区分を予測、2段階目はそれらと物性値2種を入力として岩石区分を予測。手法的には半教師あり学習になりそうです。
1段階目をミスると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 とは異なり、これも物理的関係は考慮されていません。機械学習としてブラックモデルが分布確率を提示します。データ駆動型の手軽さと解釈の難しさを併せ持っています。
手法としてはアンサンブル（スタッキング）で判別性能を高める工夫がなされています。
日本でブラックボックスモデルは受け入れられるでしょうか？
近年の文献では、データ駆動型の方が物理ベースのモデルに比べて性能の良くなる傾向が認められると言われています。が、地質モデル作成の場合はトレーニングに利用できるボーリング数が補間するボクセル数に比べ、圧倒的に少なくなります。その場合は物理ベースモデルを利用し、その結果を地質屋さんが解釈し分布を推定する方が性能が良くなるでしょう。
その「解釈」部分のロジックを、先の文献のように損失関数に組み込んだペナルティとして数式で表現すれば、機械学習でも良い地質モデルを作れるようになると思われます。プログラミングのように、自分が何を考えて推定したのかを整理、伝達することから始める必要があるということです。
地震時崩壊個所を、2手法で迅速に特定したという報告です。
ソースと機械学習手法の組み合わせは以下の通り。
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では天候に大きく影響されない点や、広域でも一人で対応できる点が魅力です。そのうち性能は向上するでしょうから、災害直後の道路ネットワークの寸断箇所（孤立集落）の確認や天然ダム形成の確認に使える日が来るかもしれません。
航測会社も頑張って開発されていると思われますので、この分野で良い勝負がなされるよう期待したいところです。
日本の研究者が日本の Landslide Susceptibility Map を AHP を用いて作成されています。
崩壊データとしては、「地すべり地形分布図」を利用されています。主として30～40年の期間しか扱えない災害データよりも、数百～数千年程度のイベントを反映した地すべり地形分布図の方が、susceptibility を見出す点で適している部分があるでしょう。
また、国外の研究者が作成するよりも地質を細かく取り扱われています。このあたりは国内の研究者の方が有利でしょう。このような文献を見ると、安心します。
この文献もそうですが、LSM を作成している文献では、基本的にはオープンデータが利用されています。検討ツールは異なるものの、どれも同じようなアプローチです。既に LSM の作成方法や考え方はある程度確立されていると言えるでしょう。
問題は、どの作り方、どの特徴量の選び方が良かったのか比較する場がないことです。文献では、当然「自分たちの提案手法は良い性能を出したよ」「今後も改良していくよ」で終わっており、他の文献と同じデータを使って手法を比較する、優位性を示すまでには至っていません。画像分類における ImageNet のように標準的なデータセットがなく、ILSVRC のように比較する場もありません。このようなデータセットや場の提供を学会等が担うべきだと思います。
機械学習を使えば信頼性の高いマップができるわけではありません。研究者が切磋琢磨できるような場ができて初めて、マップや機械学習の信頼性が向上し始めるのです。まだまだこれからです。
**************************************
20250118追記
「地すべり地形分布図」を正解として使用すると、すべり前の地形データはないのですべり後の地形（DEM）に対し特徴量を抽出することになります。すべり前の地形に対し危険度を予測するというよりは、地すべりの可能性が高い地形を危険度として示す、という解釈になるのでしょう。
データ量が多く全容を理解できませんでしたが、調査法として参考になる文献でした。これだけ調査するのは大変だったでしょう。
空中電探、水質分析、地下水年代測定等を実施することで、地下水流動経路と影響の原因を特定した、という内容でした。
国内では、建設範囲外でこれだけ広域の調査を実施することはレアでしょう。温泉や重要水源への影響検討の場合は実施する場合があるものの、空中電探まではしていないでしょう。が、やるべきなのでしょうね。いろいろな可能性を考えて調査を提案することが、後のリスク低減に繋がります。それが重要である点をこの文献は示唆しています。
実務のお手本になるような文献でした。
先日書き残しましたが、LSM 作成に transformer が利用されていました。
https://phreeqc.blogspot.com/2024/11/landslide-susceptibility-map-using-ml-5.html
transformer を利用したいと思っていたのですが、実務では画像分類での ViT を先に触りました。結果は他のアーキテクチャを抑え Best Score。また一つ、外せない選択肢が増えました。
日本では、「Landslide Susceptibility Map」に対応する統一的な用語はまだ確立されていません（産総研さんは「感受性マップ」でした）。
土砂災害に関するマップは古くから作成されていますが、様々な表現が用いられています。「ハザードマップ」と「リスクマップ」についても、明確に区別されずに使用されています。
崩壊の危険度に関するマップとして、まずは深層崩壊。国交省から公開されていますが、ここでは隆起量や地質を考慮されているようです。豪雨時のハザードマップの材料として利用可能でしょう。
●深層崩壊推定頻度マップ
過去の発生事例から得られている情報をもとに深層崩壊の推定頻度に関する全国マップを作成しました。（2010/08）
●渓流レベル評価マップ
空中写真判読等による深層崩壊の渓流（小流域）レベルの調査。（2012/09）
https://www.mlit.go.jp/mizukokudo/sabo/deep_landslide.html
地震時の斜面崩壊に対するハザードマップもあります。
この六甲式は評価手法になるのですが、得られるのはマップです。シンプルな式ですが、経験的に良い性能を発揮します。
●六甲式
国総研資料 第 204 号
参考：地震による斜面崩壊危険度評価判別式「六甲式」の改良と実時間運用
さらに一歩進んだリスクマップも公開されています。国内では崩壊データが整備されていませんので、土砂災害データが利用されていたり、民家の有無も考慮されたうえでの評価がなされています。
●土砂災害発生確率マップ（案）
国総研資料 第 1120 号
●ハザードマップポータル
ハザードマップポータルサイト
●キキクル
気象庁 | キキクル（危険度分布）
現状、機械学習を導入するだけで良いマップを作ることができるとは考えていません。おそらく、よく検討された六甲式や土砂災害確率マップの方が実態に近い結果を得られると思います。
将来、過去の崩壊データが整備される、あるいはリアルタイムで整備されるしくみが整い、人の手に負えなくなる量が毎年加えられる、というようなことになるかもしれません。そのような場合に、重回帰に代わってDNNを使う、適中率と捕捉率の同時向上＝F1マクロを評価指標としてトレーニングする、などという取り換えが可能となり、良質なマップができるようになるのでしょう。
国内では、機械学習による LSM 作成を試しながらも、データ収集・提供の仕組みを整える必要があるのでしょう。
Landslide susceptibility evaluation result in Italy (a) and Japan (b)
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でプロットしてみましたが、粗くて途切れていました。上記のマップはかなり引いていますので、雰囲気を比べられる程度でしょうか。
2024年の文献です。
特徴量選択にアンサンブル（多数決）を利用しています。個々のモデルでは、RFECVを利用。さらに、ベースモデル ‐LRによるスタッキングを通して最終的な出力を得ています。
こちらは transformer を利用されています。
pre-training 後に fine-tuning を実施したところ、他の手法より AUC がよかったよ、という報告です。
transformer ではグローバルな発生データを活用する pre-training が可能であり、データ量の多さを活かして高度な特徴学習ができそうです。それをローカルの発生・非発生データで fine-tuning することで、未知のデータに対する汎化性能を保ちつつ、ローカルな特性を捉えたモデルを構築できます。ローカルのために世界のデータを利用するという報告は"その2"で書き残しました。が、この場合は XGBoost ですので pre-training の概念がないですし、学習時に利用するにしても相応の非発生データが必要になります。
利用する特徴量を決めておいて、世界でデータを整備しておけば、事前学習済みデータとして配布・利用できそうです。幸か不幸か国内はこれからですので、特徴量に使える国内データの整備が進むとありがたいですね。
時系列 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 作成のためのデータは整備していただきたいものです。そうすると土木に携わっていない機械学習エンジニアが参加しやすく、多様な目的のマップが作成されるようになるでしょう。その時々の、最新のアーキテクチャで。
昨日、産総研さんの第41回地質調査総合センターシンポジウムを拝聴しました。
地質に関するCIM やデジタル化の話、九州のハザードマップ作成についての報告がありました。
ハザードマップは産総研の方の報告でした。その中には LSM の話もありました。産総研では LSM を「地すべり感受性マップ」と訳されているようです（が、この訳は一般的ではないとおっしゃっていました）。スライドには「崩れやすさマップ」という記載もありましたので、今後はわかりやすい後者に落ち着くのかもしれません。
LSM 作成の話がメインのように思えましたが、表題はハザードマップで、しかも使用されているインベントリは被災を伴う災害データベースであり、作成されたマップはリスクマップに近いようでした。「感受性とハザードはどう違うのか」と質問された方がいらっしゃいましたが、「いつ、どこで、まで答えるのがハザード」というような回答でした。また、災害データを使われた理由については、「被災の有無にかかわらない崩壊データも一部で使っている。ないものは使えない」といった回答でした。
それをどうするか、は研究者ではなく技術者の仕事、でしょうか？Googlig では産総研さんからのリスクマップ作製に関する委託業務がいくつか引っかかりますので、初めからLSMでなくリスクマップを作成する目的だったのか、単に予算の問題だったのかもしれません。
九州地域における斜面災害リスク評価主題図改良業務：SVMを利用
九州北部地域の地質情報解析によるリスクマップ(案)の作成業務：災害特性に着目
国内では崩壊データを整備する仕組みがありません（国が災害データを収集する仕組みはあります）。災害データを使うと、民家に近いところが危ない、危険区域や危険個所に指定されている範囲が危ない、という結果になります。これは当然で、あちこち崩壊しても被災した箇所しか国には報告されない、報告外の箇所は未崩壊として誤って扱われてしまいやすいというのが理由です。そのインベントリを教師データに使用すれば、属性としての地質よりも区域等が効いてきます。機械学習モデルの中で地質の重要性が薄れてしまうのです（地すべりはそれでも地質が効くのでかなり重要なのでしょう）。
では、LSMを作成するのに災害データしか入手できない場合はどうするか？
①航空写真、衛星写真等から作成する。
産総研さんの「一部」というのはこれかもしれません。時間と費用がかかるものの、現状では最も精度の高い方法と言えるでしょう。写真等から抽出するための画像処理や機械学習が進化すれば解決すると思います。が、現状では最後に人手が必要です。もう少し待ちましょう。
②文献で利用されているデータを利用する。
昨日の文献通りです。ローカルの評価にデータが不足すれば、世界中から集めてくればよいのです。文献個々のサプリメントデータを集めるのも有効です。機械学習にかける際には多少の工夫が必要になると思いますが、データ不足には効果的です。
③LSMをあきらめてリスクマップを作成する。
この場合、LSMを作成したうえで、ある地域に台風〇〇号並みの豪雨が来たらどうなるか、その際に被災しやすい場所はどこか？といったような手順は踏めません。被災関係なく豪雨や地震で崩れそうな場所を予測することも正確にはできないのですが、被災しやすい箇所は推定できます。
産総研さんの報告からは、①と③のいずれの方法をどの程度採用しているのかは明確ではありませんでした。聞く限りでは作成されたハザードマップ？リスクマップ？よりは収集された属性データを公開していただく方が、後々利用目的に応じて柔軟なマップ作りができるように感じました。
2019年～2024年まで、機械学習を利用したLSM作成に関する文献を20編程度集めました。
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 は作成しやすくなります。
”water level machine learning” + Googlimg で多数引っかかる内容です。
頭から読んでみたところ、地すべりのSusceptibility Map の作成の文献とは傾向が異なっており、アルゴリズムの比較が多く目につきました。データセットによって予測性能の良いアルゴリズムが異なる点はよく知られていますが、見事にバラバラ。それらのアンサンブルが少ないのは機械学習分野の流れと異なりますが、浸透するにはもう少し年月が必要なのでしょう。
内容的には玉石混合ですが、アルゴリズムの比較という点では有用ですので、引っかかったものから書き残しておきます。
Deep Machine Learning-Based Water Level Prediction Model for Colombo Flood Detention Area
Modeling the fluctuations of groundwater level by employing ensemble deep learning techniques
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.
ま、データがたくさん必要なのは機械学習も同じです（いえ、それ以上）。が、仮定の必要がなく、ハイパーパラメータの調整さえすればあとは自動でモデルを構築してくれる、ほぼリアルタイムで予測結果を算出できる、しかも予測性能が高くなるとなれば、使わない手はないでしょう。合理性の観点から国外で進んだのかもしれません。
続く。
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)
ブラインド信号源分離（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()
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)
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が簡単です。
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と何を組み合わせるか、アイデアはたくさんあると思います。基本ツールの結果を組み合わせ、本来解くべき問題に早々に挑めると、少しは前進するのではないでしょうか。
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に応じて地すべりの大きさを扱える点が魅力です。日本は完全に置いてけぼりです。
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 について、大まかな流れは理解できました（細かい点はわからないことばかりで先は長い）。
以下、流れの抜粋です。( ..)φ
狭帯域レーザーの出力は、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．選択した距離間隔（ゲージ長）にて位相差を空間微分。
※ゲージ長全体にわたる位相の変化をレーザーショット数の関数として記録することも可能。これにより、ひずみの時間微分（ひずみ速度）が得られる。これらの処理後、ファイバーに沿った距離と遅い時間（つまり、プローブパルスの連続に見合った時間スケール）の関数として、アンラップされた位相の画像が利用可能になる。
昨年購入し、何度か読むにつれ徐々に頭に入ってきました。基本原理についても最初に記載があります。
１章
<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はジオフォンからのデータよりも本質的に低くなる。しかし、高密度で取得可能であるため大幅に改善可。
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)]/Lwhere 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.
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 SPHの実装が概ね完了しました。
実装を通し、理論と実際を理解できました。
現状、2点の問題があるようです。
１つ目は seepage force (drag force) の計算式。
これまでの文献にも、DualSPHysicsでも同じ式が使われています。が、透水係数の大きな場合にしか実用的ではないようです。1e-4㎝/sオーダーだと、γ/kが大きくなりすぎて、少しの速度差が大きな拘束力を生みます。土塊がすべろうとしても、地下水がその動きに追随するまで拘束されるような挙動も見受けられました。
２つ目は密度変化。
土の間隙内の水が地表に浸出した場合、水相の密度は大きく変化しますが、この表現が困難。いえ、pressureの式で参照している基準密度を変更すれば容易に反映できるのですが、粒子が重なり合い、挙動がやや不安定でした。初期配置も面倒で、実用的ではないですね。水収支はともかく、土の動きだけを正しく反映したい、ということであれば計算の軽い single-layer two-phase formulation でよいと思います。
土相側でも層毎に密度を変えている場合、その境界で密度が訛るのは現状では仕方がないと考えています。が、この密度変化の問題は将来的に何とか解決されそうな気がします。たとえば、土層別にフラグを付けて、個別に計算するとか。それが実現象を正しく反映するようになれば、より実務に浸透していくのでしょう。
今年前半より取り組んできた SPH コードの修正と実装は、これで一旦終了です。
結果的には単相での計算はOK、2相は「定性的評価」まで、土木分野の実務で求められる諸問題には適用が難しい場合が多いと判断されます。が、もう少しだと思います。
将来に期待しましょう。
振り返ると、double-layer での実装に取り掛かってから１か月以上かかっています。
以下を参考に始めたのですが、試行錯誤の結果、最終的にはPersianSPHの文献と同じような形に落ち着きました。
A coupled fluid-solid SPH approach to modelling flow through deformable porous media - ScienceDirect
片道４時間、日帰りで、観測システムを修復してきました。
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/*
ひとまず周辺環境を含め目についた問題は全てつぶしておきました。今回は停電により、隠れていた問題も同時に顕在化しました。最初に見逃さないようにしないと、このように不要な対応が発生します。気を付けましょう。
前回から約1か月かかりましたが、２相２層での水圧計算に目途が立ちました。
https://www.sciencedirect.com/science/article/pii/S002076831730286X
上記の78-82式を実装することで解けるようにはなったのですが、まだ土層を出た後の水の密度変化までは安定して解けません。最初の含水率のまま解くことはできるのですが、空中でそれは誤りですし、勝手に密度を上昇させることもできません。圧力が下がって自然に密度が上昇する、というように計算させるとすぐに発散します。そのあたり、この論文では明記されていません。土層の計算だけ追うのであれば、１層での計算が良いと思いますが、まだ何か見落としているようです。
ひとまず、次の drag force を安定して解けるか見てみましょう。
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
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ならではでしょう。
弾塑性（一部流体）の挙動はコチラ。端から土圧により崩壊し始めます。
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”.
当たり前すぎて一般的には解説されていないのでしょうね。初心者には思いつかないノウハウです。
密度でなく圧力を変更。数行の追加で計算時間もかからず、簡単に解決できました。
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
