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.

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

hDVS

An Introduction to Distributed Optical Fibre Sensors

先の洋書よりもハードウェアや原理について詳しく書かれています。先日使用した 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．選択した距離間隔（ゲージ長）にて位相差を空間微分。
※ゲージ長全体にわたる位相の変化をレーザーショット数の関数として記録することも可能。これにより、ひずみの時間微分（ひずみ速度）が得られる。

これらの処理後、ファイバーに沿った距離と遅い時間（つまり、プローブパルスの連続に見合った時間スケール）の関数として、アンラップされた位相の画像が利用可能になる。

Distributed Acoustic Sensing in Geophysics

Distributed Acoustic Sensing in Geophysics: Methods and Applications 

昨年購入し、何度か読むにつれ徐々に頭に入ってきました。基本原理についても最初に記載があります。

１章

<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はジオフォンからのデータよりも本質的に低くなる。しかし、高密度で取得可能であるため大幅に改善可。

DAS その2

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)]/L

where 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. 

２段階認証

 google の2段階認証アプリに異常が発生し、他のアプリに切り替えました。

その間、約１週間。ひとまず会社とGoogleの２アカウントは復帰できましたが、それ以外は不具合の発生しているアプリの中で眠っています。ここにもようやくたどり着きました。

GitHub がまだ眠っています。ひとまずローカルからアップはできるのでしばらく待つつもりです。

怖いですね。