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