simpeg

SimPEG - 地球物理シミュレーションとインバージョン

クイックリファレンス

from discretize import TensorMesh
from simpeg.electromagnetics.static import resistivity as dc
from simpeg import maps, data_misfit, regularization, optimization
from simpeg import inverse_problem, inversion, directives
import numpy as np

# Create mesh
hx, hz = np.ones(100) * 10, np.ones(50) * 5
mesh = TensorMesh([hx, hz], origin='CN')

# Forward model
simulation = dc.Simulation2DNodal(mesh, survey=survey, sigmaMap=maps.ExpMap(mesh))
dpred = simulation.dpred(model)

# Inversion
dmis = data_misfit.L2DataMisfit(data=data, simulation=simulation)
reg = regularization.WeightedLeastSquares(mesh)
opt = optimization.InexactGaussNewton(maxIter=20)
inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)
inv = inversion.BaseInversion(inv_prob, directiveList=[...])
mrec = inv.run(m0)

主なクラス

クラス	用途
TensorMesh, TreeMesh	離散化（規則格子、適応型オクトリー）
Survey	データ取得ジオメトリ
Simulation	フォワードモデリングエンジン
Data	観測データ／予測データのコンテナ
InvProblem	誤差関数、正則化、最適化を統合

基本操作

メッシュの作成

from discretize import TensorMesh

# 2D mesh (x, z) - centered in x, top at z=0
hx, hz = np.ones(100) * 20, np.ones(50) * 10
mesh = TensorMesh([hx, hz], origin='CN')

# 3D mesh
mesh = TensorMesh([np.ones(50)*25, np.ones(50)*25, np.ones(30)*10], origin='CCN')

DC 比抵抗測量

from simpeg.electromagnetics.static import resistivity as dc

elec_locs = np.c_[np.linspace(-95, 95, 20), np.zeros(20)]
source_list = []
for i in range(17):  # dipole-dipole
    rx = dc.receivers.Dipole(elec_locs[[i+2]], elec_locs[[i+3]])
    src = dc.sources.Dipole([rx], elec_locs[i], elec_locs[i+1])
    source_list.append(src)
survey = dc.Survey(source_list)

フォワードモデル

model = np.ones(mesh.nC) * 100  # 100 ohm-m
simulation = dc.Simulation2DNodal(mesh, survey=survey, sigmaMap=maps.ExpMap(mesh))
dpred = simulation.dpred(np.log(1/model))  # input: log(conductivity)

インバージョン

from simpeg import data_misfit, regularization, optimization
from simpeg import inverse_problem, inversion, directives, data

obs_data = data.Data(survey, dobs=dobs, standard_deviation=0.05*np.abs(dobs))
dmis = data_misfit.L2DataMisfit(data=obs_data, simulation=simulation)
reg = regularization.WeightedLeastSquares(mesh, alpha_s=1e-4, alpha_x=1, alpha_z=1)
opt = optimization.InexactGaussNewton(maxIter=20)
inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)
dir_list = [directives.BetaSchedule(coolingFactor=2), directives.TargetMisfit()]
inv = inversion.BaseInversion(inv_prob, directiveList=dir_list)
mrec = inv.run(m0)

よく使うマップ

マップ	説明	使用例
IdentityMap	変換なし	磁化率、密度
ExpMap	exp(m)	対数パラメーター化された導電率
ReciprocalMap	1/m	比抵抗から導電率への変換
Wires	モデルを分割	ジョイントインバージョン

物性値の範囲

物性	一般的な範囲	単位
比抵抗	1 - 10000	ohm-m
導電率	0.0001 - 1	S/m
磁化率	0 - 0.1	SI
密度コントラスト	-1 to 1	g/cc

代替手段との使い分け

シナリオ	推奨
複数手法による地球物理インバージョン（DC、磁気、重力、EM）	SimPEG - 最も幅広い手法に対応
標準電極配置による浅部 ERT	pyGIMLi - よりシンプルな API、電極配置を標準サポート
GUI エクスポートを伴う ERT 中心のインバージョン	pyGIMLi - ERT 専用ツールが充実
柔軟な物理モデルを用いたカスタムフォワードモデリング	SimPEG - モジュール設計で拡張が容易
複数の地球物理データセットのジョイントインバージョン	SimPEG - Wires マップで標準サポート
商用 ERT 処理	Res2DInv / Res3DInv - 業界標準

SimPEG を選ぶ場合：複数の地球物理手法に対応する統一フレームワーク、カスタムフォワード演算子、または研究用途に適した柔軟性が必要な場合です。モジュール設計（メッシュ + 測量 + シミュレーション + インバージョン）は、複雑な問題や標準的でない問題に適しています。

SimPEG を避ける場合：標準的な ERT インバージョンだけが必要な場合（pyGIMLi のほうがセットアップが速い）、またはすぐに使える商用ソリューションが必要な場合です。

よくあるワークフロー

測量データから DC 比抵抗インバージョンを実行する

• 電極の位置を定義し、双極子-双極子（またはその他の）測量ジオメトリを作成する
• 適切なセルサイズで TensorMesh または TreeMesh を作成する
• メッシュ、測量、ExpMap を指定して dc.Simulation2DNodal を設定する
• 標準偏差とともに観測データを data.Data に読み込む
• L2DataMisfit、WeightedLeastSquares 正則化、最適化器を設定する
• 指令を設定する：BetaSchedule、TargetMisfit
• BaseInvProblem と BaseInversion を構築する
• 均質な初期モデルを使い、inv.run(m0) でインバージョンを実行する
• 復元モデルをプロットし、観測データと予測データを比較する
• データ誤差の収束を確認する（目標カイ二乗値は約 1）

ヒント

1. 正の値をとる量（比抵抗、磁化率）には対数パラメーターを使用します
2. 粗いメッシュから始め、初期テスト後に細分化します
3. 観測データと予測データをプロットして、データ適合度を確認します
4. データ適合度とモデルの滑らかさのバランスをとるよう、正則化を調整します
5. 3D 問題では効率を高めるために TreeMesh を使用します

参考資料

• 測量タイプ - 測量設定と受信器タイプ
• メッシュタイプ - メッシュの離散化と細分化

スクリプト

• scripts/dc_inversion.py - DC 比抵抗インバージョンの完全な例