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.

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