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')

直流电阻率测量

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

何时使用及替代方案

场景	建议
多方法地球物理反演（直流电、磁法、重力、电磁法）	SimPEG - 方法覆盖面最广
使用标准阵列的近地表 ERT	pyGIMLi - API 更简单，内置阵列支持
专注于 ERT 且需要 GUI 导出功能的反演	pyGIMLi - ERT 专用工具更完善
使用灵活物理模型的自定义正演建模	SimPEG - 模块化设计，易于扩展
多种地球物理数据集的联合反演	SimPEG - 通过 Wires 映射内置支持
商业 ERT 处理	Res2DInv / Res3DInv - 行业标准

选择 SimPEG 的情况：需要适用于多种地球物理方法的统一框架、自定义正演算子或研究级灵活性。其模块化设计（网格 + 测量 + 模拟 + 反演）适合复杂和非标准问题。

避免使用 SimPEG 的情况：只需要标准 ERT 反演（pyGIMLi 的配置速度更快），或需要开箱即用的商业解决方案。

常见工作流

根据测量数据进行直流电阻率反演

• 定义电极位置，并构建偶极－偶极（或其他）测量几何结构
• 使用适当的单元尺寸创建 TensorMesh 或 TreeMesh
• 使用网格、测量和 ExpMap 设置 dc.Simulation2DNodal
• 将观测数据及标准差载入 data.Data
• 配置 L2DataMisfit、WeightedLeastSquares 正则化和优化器
• 设置指令：BetaSchedule、TargetMisfit
• 构建 BaseInvProblem 和 BaseInversion
• 使用均匀初始模型运行反演：inv.run(m0)
• 绘制恢复模型，并比较观测数据与预测数据
• 检查数据失配的收敛情况（目标卡方值约为 1）

提示

1. 使用对数参数表示正值物性（电阻率、磁化率）
2. 从粗网格开始，并在初步测试后细化
3. 检查数据拟合情况：绘制观测数据与预测数据
4. 调整正则化，平衡数据拟合和模型平滑度
5. 在三维问题中使用 TreeMesh，以提高效率

参考资料

• 测量类型 - 测量配置和接收器类型
• 网格类型 - 网格离散化和细化

脚本

• scripts/dc_inversion.py - 完整的直流电阻率反演示例