文献：3D Wave‐Equation Dispersion Inversion of Distributed Acoustic Sensing Data

3D Wave‐Equation Dispersion Inversion of Distributed Acoustic Sensing Data to Reveal the Shallow Near‐Surface Geology - Liu - 2025 - Journal of Geophysical Research: Solid Earth - Wiley Online Library

AI要約

3D波動方程式分散インバージョン（3D-DAS-WD）法の解析手順

3D-DAS-WD法は、DASひずみデータに直接基づく3D波動方程式分散インバージョン手法であり、ひずみ成分を速度に変換する際の数値誤差を回避し、低SNR（信号対雑音比）のDASデータにおけるFull Waveform Inversion（FWI）の収束問題に対処することを目的としています。この手法は、浅い地表付近の地質構造を高分解能で3Dイメージングする能力を示しており、地下の空洞などの特徴を効果的に特定できます。

1. DASの基本原理とデータ取得

• DASの測定原理: 従来の地震計が粒子速度を測定するのに対し、DASは光ファイバーのゲージ長に沿ったひずみ速度を測定します。

2. 初期S波速度モデルの作成

• 解析の出発点として、初期S波速度モデルを準備します。このモデルは、反復的なインバージョンプロセスを通じて改善されていきます。

3. 予測ひずみ速度DASデータの順モデリング

• 作成した初期S波速度モデルを用いて、forward modeling を実行し、予測されるひずみ速度DASデータを生成します。この研究では、スペクトル要素法に基づく弾性モデリングコードSPECFEM3Dが使用されています。

• 計算される速度成分は、式(14)および(15)を使ってDASひずみ速度成分に変換されます。

4. 観測データと予測データの分散曲線の導出

• MASW として分散曲線を作成します。

• 分散曲線は、linear Radon transform を用いて計算された分散スペクトルの最大値に基づいて抽出されます。

• この研究では、各ショットから10度間隔で分散曲線が抽出されます。

5. 目的関数の定義と勾配の計算

• 目的関数（Objective Function）: 3D-DAS-WD法は、観測された波数κobs (θ,ω)と予測されたDASデータ表面波分散曲線κ(θ,ω)の二乗差の合計を最小化します。目的関数Jは以下のように定義されます: 

J = (1/2) ∑θ ∑ω (κ(θ,ω) − κobs (θ,ω))^2

ここで、θは方位角、ωは周波数です。

• S波速度に対する勾配: S波速度s(x)に対する目的関数の勾配δ(x)は、以下の式で計算されます: 

δ(x) = ∂J / ∂s(x) = ∑θ ∑ω Δκ(θ,ω) (∂κ(θ,ω) / ∂s(x))

• Fréchet 微分の導出: Fréchet 微分 (∂κ(θ,ω) / ∂s(x)) は、予測データと観測データのk-ω領域における相関を示す接続関数（connection function）Φを用いて導出されます。これは、随伴状態法（adjoint state method）を使用して計算されます。最終的にS波速度 s を λ (ラメの第一定数) と μ (剛性率) に関連付けることで、Fréchet 微分が導出されます

∂κ(θ,ω) / ∂s(x)

 = − (1/A) R{∫ ∂D(g,ω)/∂s(x) D̂*obs(g,θ,ω) dg} 

 = − (1/A) R{∫⟨GR ∂u(g,ω)/∂s(x), D̂obs(g,θ,ω)⟩ dg}

= − 4ρs (∂κ(m) / ∂λ) + 2ρs (∂κ(m) / ∂μ) 

・∂κ(m) / ∂λ、∂κ(m) / ∂μは、 波場u(m)（合成地震計データ）とバックプロパゲーションされた残差u'（D̂∗obs(g,θ,ω)をソースとする）の内積の積分として表される。

6. 速度モデルの更新

• ステップ長の決定: safeguarded backtracking line search method を用いて、モデル更新のためのステップ長αが決定されます。

• モデルの反復更新: S波速度モデルは、計算された勾配とステップ長αを用いて、steepest-descent formula により反復的に更新されます。 ここで、

7. 収束判定と最終S波速度モデルの出力

• ステップ3から6のプロセスは、目的関数が指定された誤差閾値以下になるか、または最大反復回数に達するまで繰り返されます。

• 収束条件が満たされた後、最終的なS波速度インバージョンモデルが出力されます。

さすがにAI要約のみでは理解できず、数式を追うのに時間がかかりました。

DAS（1軸センサー）を利用して、どのようにレイリー波を取り出すか、取り出した後、どのようにインバージョンをかけると精度が落ちにくいか、という内容でした。地震分野の知識の蓄積を感じる文献です。

3次元解析を目的としてファイバーをクロスで張るなら、通常のピックを使った方が楽なように感じましたが。ま、一つの方法として覚えておきましょう。

論文 DASによる震源決定法

Source location of volcanic earthquakes and subsurface characterization using fiber-optic cable and distributed acoustic sensing system | Scientific Reports

震源決定

・2手法の併用
・両者を使用して震源位置を決定し、最も可能性の高い結果を採用。
・解析ではS/N比が4より大きい波形データを使用。
・グリッドサーチでは、すべてのグリッドポイントの到着時間差を事前に計算しておくことにより、ほぼリアルタイムで震源決定可能。

１．到着時間差法
道路のカーブを利用し、L字型アレーとして解析。
入射地震波の遅延性(伝搬方向と入射角)を推測。
30.6〜71.4mの組み合わせを利用。コヒーレンスの高い波を利用
到着時間差は、クロススペクトルの位相差より計算。
クロススペクトルは主に振幅の大きい波から計算されるため、到着時間の差はS波の到着と関連していると仮定した。
波形に12秒の時間ウィンドウを適用。
グリッドサーチを利用。

2. 振幅ソース位置(ASL)法
光ケーブルに沿った地震波の最大振幅を測定し、火山性地震の発生源を特定。
ASL法では、発生源から発生する地震波は、幾何学的な広がりと固有の減衰によって減衰すると仮定。
サイト増幅係数は、地域性構造地震のコーダ波の解析から決定。
火山性地震の最大振幅を光ファイバーケーブルに沿ったすべての測定点で読み取る。
Qを10〜100とすると、結果は大きく変化しない。20を使用。
グリッドサーチの利用。

到着時間差の測定
光ファイバーケーブルは長距離をカバーするため、到着時間差を計算する測定ポイントのペアの間隔距離には多くの選択肢がある。
最大間隔距離をS波の波長（４Hz）の約4分の1に設定。
ケーブルに沿って90m以内に位置する一対の測定点の波形を解析。
波形の相関を調べるために、コヒーレンスが0.5より大きい波形を使用。
ケーブル方向が異なる測定ポイントでの動的ひずみ信号は極性が異なる場合があるため、ケーブルの方向が互いに30°未満で一致している測定ポイントのペアのみを使用。
ケーブルに沿って91.8m以内にある測定ポイントのペアについて、9つの到着時間差を計算した。

次の手順を適用して、どの間隔が適切であったかを判断。
(1)10.2mから91.8mまでの9区間を10.2m刻みでデータを用意し、
(2)最適な位置を取得し、各区間の分散減少を評価。
(3) 分散の減少が10%未満の区間のデータを削除。
(4)その後、ステップ(2)からこのプロセスを繰り返した。
(2) から (4) までのプロセスは、差異削減が10%を超えたときに終了。
30〜70mの間隔を決定。

部位増幅率の推定
直達S波に続くコーダ波は、不均質な構造が地震波を十分に散乱させると、空間的に均一に分布する。
サイト増幅係数は、基準局に対するターゲットサイトにおけるコーダ波の相対振幅より推定。
約50〜200 kmの距離と約10〜100kmの深さで発生するマグニチュード3以上の地殻変動地震を11回利用。
ケーブルに沿った測定点における二乗平均平方根の振幅を計算し、基準点(地震計が仮的に設置されたMP1027)の振幅と比較。
S波走時の2倍以上の経過時間を持つコーダ波を使用。
5秒ごとに0〜50秒の経過時間の相対振幅を推定し、各測定ポイントでの平均振幅を分析のサイト増幅係数として使用した。
相対振幅の標準偏差は非常に小さく、サイトファクターが確実に評価されていることを示す。
サイト増幅率は、光ファイバケーブルに沿って対数スケールで約1のオーダーで変化するため、ASLではこの補正が極めて重要である。

地すべりブロック把握のための DAS 適用

地すべり調査に DAS を利用する方法です。

Previously hidden landslide processes revealed using distributed acoustic sensing with nanostrain-rate sensitivity | Nature Communications

文献では平面計測のみですが、ひずみ計のかわりに途中でボーリング孔を経由する（埋め込む）ことも可能です。地すべり調査と DAS をご存じの方なら誰もが思いつくレベルの調査法なのですが、国内の実務での実施例は見たことがありません（研究で平面のみ、ボーリング孔のみは見たことがあります）。
問題は費用。千万円単位で大きくなるので発注するにも容易ではないのでしょう。精度の高い技術があることは分かっていても高価すぎる、今まで通りでも良い、といったところでしょうか。

DAS、安くなりませんかね。

DASのピッキング

Seismic arrival-time picking on distributed acoustic sensing data using semi-supervised learning | Nature Communications

DASのピッキングに関する文献です。
従来の STA/LTA 等ではうまくいかないDASのピッキングに機械学習を使ったという内容です。

この文献、査読結果が一緒に公開されています。オープンサイエンスもここまで来たのだなと感心しました（他の文献でも公開されていますので、私の気づきが遅かったのでしょう）。

その中の指摘で、近年の文献を見ながら感じていた部分が明示されていました。

I do understand that deep learning is currently a big hype and that it is tempfing, especially for young researchers, to become part of this. However, the days where the mere applicafion of some deep　learning model was science, are over. With numerous easily usable programming tools, the training of a neural network (deep or not; this is just somewhat hollow semanfics) has become an easy task for undergraduate students. After all, it is what it is: fifting the coefficients of some funcfional form to a collecfion of data.
This said, it is not clear how exactly the authors go significantly beyond just training yet another deep neural network? What is the innovafion that goes beyond the obvious?

LSM作成でも、機械学習分野でSOTAを達成した手法の流用は見ていて面白いのですが、それまでですよね。できたLSMが実際にどの程度使えるのかがわかりませんので。GISを使って特徴量を作り、機械学習にかけるというのは、学生でもできます。今年のGCONでは高専生も頑張っています。
第3回高専GIRLS SDGs×Technology Contest(高専GCON2024)

ピッキングも同じなのでしょう。他にもいくつか機械学習を利用している文献を見かけます。が、機械学習を使えば精度が上がった、だけの時代はもう終わったと認識せざるを得ません。

ヘテロダイン検波

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が簡単です。

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. 

DAS

18～20日まで、鹿児島で火山学会の秋季大会がありました。

私は参加していませんが、同僚が参加しており発表要旨を見せてもらいました。

と、俄然、気になるモノが。

DAS (Distributed Acoustic Sensing) で火山性地震や土石流を検知した事例です。
これ、最近の流行ですよね。機材が欲しいのですが、高価なので指をくわえて眺めています。

DASは既設の光ファイバーを利用できるので、接続する機材のイニシャルコストに目を瞑れば非常に手軽です。既に全国に受信器が設置済みと言えますので。そのためか、国内でも地震や火山以外に、いくつかの用途に向けて研究されているようです。例えば、河川。

国総研：河川堤防の変状検知システム実験結果nilim.go.jp/lab/fbg/gijyutsukoubo/jikkennkekka/report/10.enaa.pdf

中国地勢：分布型光ファイバセンサによる河川堤防の 観測方法の検討について
https://www.cgr.mlit.go.jp/tyokugi/74/d4/05.pdf

数百ch、大きいと 1000Hz でデータを取られているので、処理はヘビーです。空間方向を時間方向になおして波形を作る前処理も必要になることがあります。が、Pythonでも処理は可能です。

地すべりや物理探査、実務でも使えるのですが、市場が大きくなるまでは皆さん様子見でしょうね。先行投資できる体力のある会社が羨ましい。