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()