こんにちは、DX Solution事業本部 Devグループの岩永です！
地下の可視化とは、どのようなことかわからない人も多いでしょう。地下の状態は、地上と違い簡単に観察を行うことが難しく、調査データを行うにもコストが発生するため調査・観測も最小限に制限されます。その少ないデータから可視化するために色々な手法が研究されています。地下の可視化は、土木分野において、トンネル、地下鉄、橋梁の基礎、ダムなどの大規模インフラ建設に不可欠です。
その主な目的は、地盤の強度、断層、湧水の状況を正確に把握することです。これにより、設計段階で最も安全かつ経済的な施工ルートや工法を選定することが可能になります。最終的に、工事中の崩落や浸水といった事故を防ぎ、インフラの長期的な安定性を確保します。
今回は、地震波のデータを利用するGeophysical Waveform Inversionを用いた地下の可視化技術を紹介します！
Geophysical Waveform Inversionとは、観測された地震波データ（波動場）と、コンピューターシミュレーションによって計算されるデータとの差（残差）を最小化することで、地下の物理探査モデル（速度構造など）を高精度に推定する最適化手法です。
これまでの技術
S波反射法
地下構造を調査する弾性波探査の一種で、S波（横波）の反射波を利用して地層の境界面や地盤の物性（硬さなど）を推定する手法です。基本的な仕組みは、地表で人工的に振動（地震波）を発生させ、それが地下の地層境界などで反射して地表に戻ってきた波を、計測器で観測します。
波が往復する時間や波形を解析し、逆問題手法を用いて地下の境界面の深さやS波速度分布を推定します。この手法は、地盤の硬さを示すS波速度を高い解像度で直接的に推定できる反面、深い探査が難しく、特別な震源装置が必要なとなる手法です。
トモグラフィー法
地下探査において、地下の三次元的な内部構造や物性分布を「断層撮影」のように画像化する手法です。S波反射法と同様に地下の特定の点から地震波を発生させ、別の多数の地点で波の到達時間を計測します。送受信間の時間のこと走時といいます。 トモグラフィー手法には以下のようにいくつか種類があります。
電気探査トモグラフィー：地下に電流を流し、地上で電位差（電気抵抗）を測定することで、地下の電気抵抗率（比抵抗）の分布を画像化します。
本手法は、波の到達時間のみを利用するため、計算が高速で安定しており、非破壊で広範囲の調査が可能です。これが大きなメリットです。 一方、デメリットとして、走時情報のみに依存するため解像度が低く、波の回折や散乱といった複雑な物理現象を正確に扱えないという限界があります。
Geophysical Waveform Inversion
Geophysical Waveform Inversionは、従来のトモグラフィー法が利用する波の到達時間だけでなく、振幅や位相を含む観測波形の全情報を解析に用いる地下構造推定技術です。波の反射・回折といったすべての物理現象を考慮するため、トモグラフィーよりもはるかに高い解像度で地下の速度や密度の分布を推定でき、土木分野の他に資源探査や地震防災分野でも注目されています。
基本的に以下のステップを反復的に繰り返す最適化プロセスです。
残差計算：計算された波形と、実際に観測された実測波形との差（残差）を計算します。
勾配計算：この残差を小さくするために、地下モデルのどの部分を、どれだけ、どの方向に修正すればよいかを示す勾配を計算します。
モデル更新：勾配情報に基づいて、地下モデルを修正し、反復計算に戻ります。
実行方法
以下にライブラリPyFWIを用いて簡単な解析を示します。 2種類の速度構造があり、中央の速度帯が解析で推定可能かを検討しています。
対象説明
実装
必要なライブラリを定義します。
import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import sys sys.path.append('../src/') sys.path.append("../src/PyFWI/") import PyFWI.wave_propagation as wave import PyFWI.acquisition as acq import PyFWI.seiplot as splt import PyFWI.model_dataset as md import PyFWI.fwi_tools as tools import PyFWI.processing as process from PyFWI.tl_fwi import TimeLapse import copy
今回のモデルを作成します。
Model = md.ModelGenerator('louboutin') model = Model() im = splt.earth_model(model, cmap='coolwarm')
地震波の送受信の設定を行います。
model_shape = ml[[*ml][0]].shape inpa = { 'ns': 4, # 震源の数 'sdo': 4, # 差分法の字数 'fdom': 15, # 震源の中心周波数 'dh': 7, # 空間サンプリング 'dt': 0.004, # 時間サンプリング 'acq_type': 1, # データ取得タイプ (0: 坑井間, 1: 地表, 2: 両方) 't': 0.6, # シミュレーション時間 'npml': 20, # PML境界（反射を防ぐために用いられる層の厚さ） 'pmlR': 1e-5, # PMLの係数（吸収性能を聖女する係数） 'pml_dir': 2, # PMLの方向：境界を適用する方向、2：全方向 } seisout = 0 # 出力タイプ inpa['rec_dis'] = 1 * inpa['dh'] # 受振間距離
指定されたパラメータに基づいて、震源と受振器の位置を決定します。
offsetx = inpa['dh'] * model_shape[1] depth = inpa['dh'] * model_shape[0] src_loc, rec_loc, n_surface_rec, n_well_rec = acq.acq_parameters(inpa['ns'], inpa['rec_dis'], offsetx, depth, inpa['dh'], inpa['sdo'], acq_type=inpa['acq_type']) rec_loc[:, 1] -= 2 * inpa['dh'] # 起振点を設定 src = acq.Source(src_loc, inpa['dh'], inpa['dt']) src.Ricker(inpa['fdom']
フォワードモデリングの実行をします。
# 波動伝播オブジェクトの生成 W = wave.WavePropagator(inpa, src, rec_loc, model_shape, n_well_rec=n_well_rec, components=seisout, chpr=20) # フォワードモデリングの実行 （ベースライン観測データ） db_obs = W.forward_modeling(bl, show=False) # # 波動伝播オブジェクトの生成 W = wave.WavePropagator(inpa, src, rec_loc, model_shape, n_well_rec=n_well_rec, components=seisout, chpr=20) # フォワードモデリングの実行（現行モデル観測データ） dm_obs = W.forward_modeling(ml, show=False) #
地震波探査のシミュレーション結果（合成データ）を表示します。
fig = plt.figure(figsize=(8, 4)) count = 1 ax = fig.add_subplot(122) ax = splt.seismic_section(ax, d_obs['taux'], t_axis=np.linspace(0, inpa['t'], int(1 + inpa['t'] // inpa['dt']))) ax_loc = [1, 2, 5, 6] snapshots = [40, 80, 130, 180] for i in range(len(snapshots)): ax = fig.add_subplot(2, 4, ax_loc[i]) ax.imshow(W.W['taux'][:, :, 0, snapshots[i]], cmap='coolwarm') ax.axis('off') count += 1
起振点と受振点の配置を示します。
m0 = Model(smoothing=1) m0['vs'] *= 0.0 m0['rho'] = np.ones_like(model['rho']) fig = splt.earth_model(m0, ['vp'], cmap='coolwarm') fig.axes[0].plot(src_loc[:,0]//inpa["dh"], src_loc[:,1]//inpa["dh"], "rv", markersize=5) fig.axes[0].plot(rec_loc[:,0]//inpa["dh"], rec_loc[:,1]//inpa["dh"], "b*", markersize=3)
逆解析を実行します。
#弱すぎる反射波の情報を活用しやすくするため、勾配の更新が特定の領域（特にエネルギーが弱い深部や複雑な部分）に偏りすぎないように、エネルギーバランスを取る inpa['energy_balancing'] = True # 逆解析の準備 fwi = FWI(d_obs, inpa, src, rec_loc, model_shape, components=seisout, chpr=20, n_well_rec=n_well_rec) # 逆解析の実行 m_est, _ = fwi(m0, method="lbfgs", freqs=[25, 45], iter=[2, 2], n_params=1, k_0=1, k_end=2)
初期モデル m0 を出発点として、指定されたL-BFGS法（method="lbfgs"）やマルチスケール（freqs：25Hz 〜 45Hz）に基づきモデルを反復的に修正し、最終的な推定モデル （m_est） を求めます。
無事、最初に設定した速度構造まではいきませんが、中央に何かあるような速度構造が計算されました。今回は、解析時間も短時間に抑えましたが、設定次第では、もっと解像度を上げることができます。
最後に
今回は、地下の可視化技術の一つとしてGeophysical Waveform Inversionを紹介しました。 フルウェーブインバージョン（FWI）などの高度な逆問題解析は、その性質上、物理モデルに基づく計算に膨大な時間を要します。現在、この計算コストを大幅に削減するため、機械学習を用いてデータとモデルの関係を近似させる研究が盛んに行われています。
燈では、Geophysical Waveform Inversionのような新しい技術を調査し、効率的に運用する方法も日々追求しています。
さらに、効率的な運用方法を仕組み化して横展開することで、事業部全体としての開発効率向上も図っています。
We’re Hiring!
燈では、技術の調査・導入に興味があるエンジニアを募集しております！
興味がある方は、ぜひカジュアル面談でお話しましょう！
