简介：pyGIMLi 例子是一套开源地球物理建模与反演软件的教程示例合集，面向地球物理、地质工程及相关领域的初学者与进阶使用者。压缩包共11个文件，全部为Python脚本，大小仅19KB，便于快速下载与本地运行。内容按 _tutorials_auto_python 目录组织，涵盖基础操作、数值建模、反演与可视化等主题。基础部分演示 gimli 核心对象使用、网格插值与质量检查；建模部分包含有限元法正演（如直流电阻率、二维/一维热传导）；反演部分提供直流电阻率一维平滑与块状反演示例。此外还包括指数拟合、多项式拟合等入门脚本，帮助用户理解从数据到模型的全流程。资源已有1069人学习下载，适合希望快速上手 pyGIMLi、掌握地质建模与反演实战技能的读者。示例脚本均提供完整可运行代码，配合注释便于逐行理解，通过运行这些脚本可以直观学习代码写法与参数设置，为后续开展科研或工程应用提供可直接改造的模板。

1. 拆解 pyGIMLi 例子：从正演到反演的一套完整地球物理工作流

很多做电法、电阻率层析成像的从业者，手里缺的不是数据，而是一套能直接跑通的代码。pyGIMLi 这个库恰好补上了这块短板：它把网格生成、正演计算、反演迭代、结果可视化和模型评估全串在一条流水线上。我最初接触它是因为同事丢来一个高密度电法数据集，用传统脚本处理完还要手动转格式、手动挑反演参数，折腾一下午才出一个像样的剖面。换成 pyGIMLi 之后，从加载数据到拿到反演结果，整个过程被压缩到一个 Python 脚本里，还能顺手把灵敏度、误差分布这些常规工具不好出的东西一起画出来。适合谁用？做地质勘察、水文探测、工程检测的工程师，以及需要批量处理电阻率数据的研究生和研究员。这份资源的核心价值在于：它不止给了能跑的 examples，而是把地球物理反演里最容易被忽略的边界条件、网格密度、正则化参数这些坑，直接在代码注释里标了出来。

2. 为什么选 pyGIMLi：反演框架不是黑匣子，是半透明的工具箱

2.1 模块化架构：正演、反演、网格各自独立又互相咬合

pyGIMLi 的设计逻辑和很多商业软件不一样。商业软件往往把正演和反演焊死在一个图形界面后面，你点几个按钮就出图，但中间发生了什么完全看不到。pyGIMLi 是 Python 库，底层用 C++ 实现核心计算，上层暴露统一 API，正演、反演、网格、可视化各是各的模块，串起来用的时候又能互相传递数据。

从我的使用经验来看，这种设计最实际的好处是排查方便。正算出问题就在正演模块里查，反演不收敛就在反演参数里调，网格质量差就直接看网格对象。举一个具体场景：做二维电阻率剖面时，如果地表地形起伏很大，商业软件处理起来往往要额外买模块，pyGIMLi 里直接在网格生成阶段把地形高程数组传给createGrid或者通过Mesh类的节点修改实现。代码层面大概长这样：

import pygimli as pg import pygimli.meshtools as mt # 定义测线长度和地形高程数组 x = pg.utils.grange(start=0.0, stop=100.0, n=101) topo = pg.Vector(101, 0.0) # 先用平地表占位 for i in range(len(x)): topo[i] = 5.0 * pg.math.sin(x[i] / 20.0) + 2.0 # 构造一个起伏地形 # 创建包含地形的网格 geom = mt.createWorld( x, topo, worldBoundary=[-20, -10, 120, 20], # 四个值：左、下、右、上边界 boundaryMarkers=[1, 2, 3, 4] # 给四周边界打标记 )

这段代码的逻辑是：先用pg.utils.grange生成从 0 到 100 米的等间距测线坐标，然后基于公式构造一个起伏地形数组，最后createWorld用这些地形点生成网格几何模型。worldBoundary四个参数控制网格外边界范围，boundaryMarkers给边界编号，后续加边界条件时直接引用标记号就行。

参数层面有个细节：topo数组的长度必须和x一致，否则createWorld会直接抛异常。如果你的实际测线坐标不是从 0 开始，记得先做平移，不然地形和测线错位，反演结果会非常难看。

2.2 反演引擎：Gauss-Newton 迭代与正则化参数怎么配合

pyGIMLi 反演核心用的是 Gauss-Newton 方法，配合不同类型的一阶或二阶平滑约束。这套东西听起来很数学，实际用起来需要关心的其实就是两个量：正则化参数lam和最大迭代次数maxIter。lam控制模型平滑度和数据拟合度之间的权衡，调大了结果更平滑，调小了更贴数据但容易出现虚假异常体。

我先给一个最常规的反演调用框架：

# 以电阻率反演为例 from pygimli.physics import ert # data 是之前加载并处理好的 ERT 测量数据 mgr = ert.ERTManager(data) model = mgr.invert( lam=20.0, maxIter=10, verbose=True, robustData=True, zWeight=0.5 )

lam=20.0是很多标准例子的起始值，如果反演出来的模型太碎、异常体边缘锯齿严重，就调大lam；如果数据拟合差、残差始终降不下来，就适当减小。robustData=True表示用 L1 范数处理数据中的离群点，实测数据有异常跳变时强烈建议开。zWeight=0.5是垂直方向平滑权重，当地质体倾向水平层状时经常设成 0.5 以下。

还有个被很多人忽略的点：ERTManager会自动做数据误差估计。默认按 3% 相对误差算，如果实际仪器精度达不到，反演结果的信噪比会虚高。我会在数据加载后手动重置误差，见后面避坑章节。

2.3 网格质量对反演结果的决定性影响

网格是反演的容器。容器不对，再好的算法也白搭。pyGIMLi 里网格生成用的是无结构三角网格，相比规则矩形网格，它能把地形和异常体边界贴合得更自然。但无结构网格也意味着你要对网格密度分布有概念：电极附近的网格必须加密，远离测线的区域可以放粗。

常见做法是，先用createWorld生成一个大范围背景网格，再在电极位置用createPolygon或者createCircle嵌入局部加密区。下面这段代码演示了如何把加密区域和背景网格合并：

# 创建一个圆形的局部加密区域 circle = mt.createCircle( pos=[30.0, -5.0], # 圆心坐标 radius=3.0, # 半径 marker=5, # 区域标记 area=0.2 # 该区域最大单元格面积 ) # 合并背景网格和加密区域 merged = mt.merge(geom, circle) # 生成最终的计算网格 mesh = mt.createMesh(merged, quality=34.0, area=2.0)

area参数控制加密区的网格单元面积上限，值越小网格越密。这里circle的area=0.2意味着加密区内每个三角单元面积不超过 0.2 平方米，而背景网格area=2.0就疏松得多。quality=34.0是最小三角形质量角，低于这个角度的单元会被重划分，通常 30-34 之间比较安全。

一个经验值：电极间距是 1 米时，电极附近网格单元面积设到 0.1-0.3 比较合理；间距 5 米时，加密区面积可以放到 1.0-2.0。网格过密会让计算时间成倍增加，反演一次从几十秒变成十几分钟，但精度提升很有限，这个平衡点需要自己根据数据量试。

3. 手把手跑通一个二维电阻率成像例子：完整流程与参数说明

3.1 数据准备：从文本文件到 pyGIMLi 数据对象

拿到任何实测数据，第一步永远是统一格式。pyGIMLi 的 ERT 模块支持几种常见格式，最常用的是 Protocol 格式：每行一条测量值，包含电流电极 A、B，电位电极 M、N，视电阻率值或者电流电压对。真实数据往往没这么干净，可能电极编号从 1 开始，坐标原点不在测线起点，甚至电极排列不是等间距的。我先展示怎么手动构造一个数据对象：

import pygimli as pg from pygimli.physics import ert # 手工构造：电极编号从1开始，4个电极组成的测量序列 elecs = [(0.0, 0.0), (5.0, 0.0), (10.0, 0.0), (15.0, 0.0)] # 每个测量四元组: (A, B, M, N) 对应电极数组下标（从0开始） measurements = [ (0, 1, 2, 3), (0, 2, 1, 3), (1, 2, 0, 3) ] # 对应的视电阻率值（单位: 欧姆米） rho_values = [120.5, 135.2, 118.7] # 转换为 pyGIMLi 数据结构 data = ert.createData(elecs, measurements, rho_values)

这段构造代码的关键点在于measurements里的四元组索引是从 0 开始，而实际仪器记录电极号通常从 1 开始，转换时记得减一，这是个一不留神就翻车的细节。createData内部会把电极坐标和测量序列绑定成DataContainer对象，后面所有正演反演都复用这个对象。

如果你的数据量很大，比如高密度电法几百条测线、上千个电极，手工构造不现实。这时直接用ert.ERTManager(datafilename)加载文件更方便，它会自动识别常见格式并完成解析。加载这块我在避坑章节会再补充一个关于单位的问题。

3.2 正演计算：验证网格与测量方案是否匹配

反演之前先做正演，这个习惯帮我挡掉很多数据问题。正演的核心作用是：给一个已知电阻率模型，计算理论上的测量响应，然后和实测数据对比。如果差得离谱，说明网格或电极方案有问题，早发现早处理。

正演调用非常直接：

# 创建一个均匀半空间模型，电阻率固定为100欧姆米 rhomap = [[1.0, 100.0]] # 网格标记1对应电阻率100 forward_scheme = mgr.createForwardScheme(mesh) # 执行正演，返回模拟数据对象 sim_data = mgr.forward(mesh, rhomap, scheme=forward_scheme) # 对比模拟数据和实测数据的视电阻率范围 print("实测数据范围:", min(rho_values), max(rho_values)) print("模拟数据范围:", min(sim_data["rhoa"]), max(sim_data["rhoa"]))

这里rhomap的写法是[网格标记, 电阻率值]的列表，网格标记在createWorld和createCircle阶段定义。forward方法返回的sim_data["rhoa"]是模拟视电阻率数组。如果模拟值和实测值差一个数量级以上，先检查单位：pyGIMLi 默认长度单位是米，电阻率单位是欧姆米，如果你的电极坐标或电阻率值用的是厘米或毫西门子每米，结果必然对不上。

还有一个常见坑：电极位置必须落在网格节点上。用createWorld生成背景网格后，如果电极坐标超出了网格覆盖的 x 范围，正演计算会跳过这些电极，导致模拟数据不完整。检查方法很简单，打印mesh.nodes()的最大最小值，跟电极坐标对照一下。

3.3 反演执行：参数矩阵与收敛判据

正演验证通过后，进入反演阶段。反演的本质是反复正演、计算残差、更新模型，直到残差不再下降。pyGIMLi 的ERTManager.invert把这一整套迭代过程封装了起来，但理解内部逻辑对调参非常关键。

一个完整的反演执行片段：

mgr = ert.ERTManager(data) model = mgr.invert( lam=20.0, maxIter=15, startModel=pg.Vector(1, 100.0), # 初始模型均匀100欧姆米 lamTolerance=0.1, # 正则化参数下降阈值 fitTolerance=0.01, # 数据拟合相对变化阈值 verbose=True ) # 查看反演结果摘要 print(mgr.inv.invModel) print("均方根残差:", mgr.inv.absRMS)

参数含义拆开看：startModel=pg.Vector(1, 100.0)构造一个只有 1 个元素的向量，初始值 100 欧姆米，pyGIMLi 会把它自动扩展到所有网格单元。lamTolerance=0.1表示当正则化参数的变化率小于 10% 时认为达到收敛；fitTolerance=0.01表示数据拟合度的相对改善小于 1% 时停止迭代。这两个容差是常见的收敛判据，如果反演一直跑不完，可以把fitTolerance放宽到 0.03。

注意：startModel的选值最好接近测区背景电阻率。如果测区是含水地层，背景可能只有几十欧姆米，你却给了 1000，反演的迭代次数会明显增加，甚至陷入局部极小值。没有先验信息时，用所有视电阻率的对数平均值作为初始模型是最稳妥的。

4. 避坑指南：反演跑不通、结果跳变、网格错位这几个真实踩坑记录

4.1 电极坐标单位不匹配，正演数据全是 NaN

现象：加载实测数据后，正演计算出的模拟数据全是nan，或者反演进行到一半直接报错退出。

原因：这是我一个同事踩过的真坑。他手里的电极坐标以公里为单位，而 pyGIMLi 默认按米处理，测线长 1.2 公里，网格范围却只有几十米，电极全部跑到网格外面去了。正演在电极位置找不到对应的网格单元，计算自然失败。

解决：加载数据后，先检查坐标范围。如果坐标数值很大或很小，主动缩放单位。

data.setSensorPositions([(p[0] * 1000.0, p[1] * 1000.0) for p in data.sensorPositions()])

把坐标统一到米以后，重新打印范围确认一下，再进入正演环节。

4.2 反演结果出现条纹状虚假异常

现象：反演完成的剖面图上，地表浅部出现交替排列的高阻-低阻条纹，看起来像斑马线。这些条纹在实测数据里没有明确对应。

原因：网格加密区域和背景网格疏密差异过大，导致反演时浅部单元权重失衡。更常见的原因是lam值设得太小，模型过度拟合了数据中的噪声。

解决：遇到条纹先别急着改网格，优先把lam调大试试。我会按 20 → 50 → 100 的顺序递进测试，同时观察均方根残差变化。如果残差变化不大，而条纹明显消失，说明是正则化约束不够。如果残差显著上升，才考虑网格问题。

4.3 反演不收敛：残差停在 10% 上下反复横跳

现象：迭代到第 5 到第 8 次时，均方根残差不再单调下降，而是在 10% 附近小幅震荡，verbose 输出里能看到数值反复跳。

原因：数据中存在离群点，单个测量值的误差远大于其他值，Gauss-Newton 迭代被这个异常值持续拉扯。另一个可能原因是初始模型远离真实模型，迭代在局部极小值附近抖动。

解决：打开robustData=True，这个参数会切换为 L1 范数约束，显著降低离群点的影响。同时把maxIter调大到 20 以上，让迭代有足够空间跳出抖动区。

4.4 网格边界条件错误导致深部模型失真

现象：反演结果的深部区域出现明显的高阻异常，即使测线下方实际是低阻含水层，也会被反演成高阻。

原因：createWorld的worldBoundary设置得太窄，底部边界离测线深度太近，边界条件对深部区域产生了强烈约束。

解决：把底部边界放到目标探测深度的 3 倍以上。例如测线长度 100 米、最大探测深度按 30 米算，worldBoundary的底部至少放到 -100 米。

geom = mt.createWorld( x, topo, worldBoundary=[-20, -100, 120, 20], boundaryMarkers=[1, 2, 3, 4] )

这个修改会略微增加网格单元总数，但深部可靠性显著提升。

5. 进阶技巧：用 pyGIMLi 做模型质量评估与可视化输出

5.1 灵敏度图：识别反演结果哪些区域可信，哪些是陪跑

反演结果不能拿到手就画图交差。我建议每次反演完都额外输出一张灵敏度图。灵敏度表示模型参数变化对观测数据的响应程度，灵敏度高的区域反演结果可信度高，灵敏度低的区域基本属于数值填充，没有实际地质意义。

pyGIMLi 里计算灵敏度并不复杂：

sens = mgr.inv.sensitivity() # 归一化到 0-1 范围 sens_norm = (sens - sens.min()) / (sens.max() - sens.min()) # 和反演结果一起可视化 ax1 = pg.show(mesh, model, coverage=sens_norm, label="Resistivity (Ohm m)")

coverage参数控制模型显示时被灵敏度阈值遮罩的程度。coverage=0.1表示灵敏度低于 10% 的区域将以半透明方式显示，这样剖面上不可信的区域一目了然。实际交付报告时，我通常把coverage设到 0.05 到 0.15 之间，太低遮不住噪声区，太高会把浅部有效信息也糊掉。

5.2 多测线批量反演：脚本化处理的一个闭环范例

当你有几十条平行测线要处理时，逐条手工配置和等待就是不现实的。pyGIMLi 的脚本化特性在这里能体现最大价值。我一般把一个批处理流程组织成三步：加载配置、循环反演、汇总输出。

import os import glob # 假设测线数据文件按 std_01.dat, std_02.dat 排列 dat_files = sorted(glob.glob("std_*.dat")) for idx, f in enumerate(dat_files): mgr = ert.ERTManager(f) model = mgr.invert( lam=20.0, maxIter=12, robustData=True, verbose=False ) # 保存反演结果和灵敏度 mgr.saveResults(f.replace(".dat", "_inverted.png")) mgr.inv.save(f.replace(".dat", "_invModel.npy")) print(f"第{idx+1}/{len(dat_files)}条测线完成，残差: {mgr.inv.absRMS:.2f}%")

saveResults会生成带色标和地形显示的反演剖面图，save把模型向量直接存成 numpy 格式，方便后续批量读取做统计对比。一个建议：批处理时verbose=False关掉逐次迭代输出，只在每个文件结尾打印残差摘要，日志会干净很多。absRMS的单位是百分比，表示数据拟合的均方根相对误差，这个值也是你判断每条测线数据质量最直接的指标。

5.3 与地质解释衔接：剖面图叠加钻孔信息

反演剖面最终要能跟地质分层对上话。我经常在 pyGIMLi 生成的图上叠加钻孔岩性分界点，用来校验反演边界和实际地质界线是否一致。

import matplotlib.pyplot as plt # 假设一个钻孔在 x=45 米处，岩性分界在 -12 米 ax = pg.show(mesh, model, label="Resistivity (Ohm m)") ax[0].plot([45.0, 45.0], [0.0, -25.0], "k--", linewidth=1.5, label="borehole") ax[0].plot([45.0], [-12.0], "rv", markersize=10, label="interface") ax[0].legend() plt.show()

坐标对应关系要注意：pyGIMLi 的 y 轴向下为正，所以深度 -12 米在垂直方向上正好显示在下方。钻孔位置如果落在测线范围之外，plot还是会把线画出来，但已经超出网格覆盖区，需要在图上留意别误解。叠加图的意义在于快速验证反演结果的可靠性，为最终地质解释提供一层独立证据。从那以后，我每次跑完测线数据，都会强制走一遍「正演验证 → 反演 → 灵敏度评估 → 钻孔对比」这个流程，少选条反演直接交图的情况就很少发生了。希望帮到你。

本文还有配套的精品资源，点击获取