機械学習による水位予測

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