简介这份资源是面向地球物理与地质工程方向学习者的 PyGIMLi 入门示例合集适合希望用 Python 开展地质建模与反演分析的初学者和科研人员。压缩包共 11 个文件全部为 py 脚本体积约 19KB按基础、建模、反演三个层次组织覆盖几何建模、有限元数值求解、一维与二维反演、网格插值与质量检查等典型环节并包含热传导与直流电法等场景的示例代码。每个脚本均可直接运行便于读者对照理解建模流程、数值方法调用方式以及反演参数设置思路从而把抽象的地球物理理论落到可执行的代码上。目前已有 1069 人学习说明该示例集在同类资源中具有一定参考价值适合作为熟悉 PyGIMLi 工作流、积累建模与反演实践经验的起点。1. pyGIMLi 例子从一条电阻率剖面到可复现反演的最小闭环很多人第一次搜 pyGIMLi 例子是因为手里已经有一份野外电法数据或者一段想验证的地下水文模型却卡在“装完之后不知道下一步点哪里”。pyGIMLi 是一个用 Python 做地球物理建模与反演的开源库覆盖电阻率、探地雷达、地震走时、重力等多种方法。它真正有用的地方不是 API 有多全而是能把“正演—反演—可视化”串成一条脚本化的流水线让一次剖面处理从手工点软件变成可版本管理的代码。这篇笔记面向已经会写基础 Python、想用 pyGIMLi 跑通第一个完整例子的从业者从环境、数据组织、正演、反演一路写到参数怎么调、哪里最容易翻车。读完你应该能独立复现一条二维电阻率反演并知道每个数字背后的含义。2. 环境与数据把 pyGIMLi 例子跑起来前的三件事2.1 安装方式的选择与依赖边界pyGIMLi 的安装是新手第一个坎。它底层有 C 核心纯 pip 安装在某些平台会缺预编译轮子常见做法是优先用 conda 环境因为 conda-forge 上有维护好的二进制包能避开编译工具链的坑。我一般会单独建一个环境不和其它科学计算库混在一起避免 BLAS、OpenMP 版本冲突导致反演中途崩溃。# 建一个独立环境Python 版本不要追最新3.10 附近最稳 conda create -n gimli python3.10 -y conda activate gimli # 从 conda-forge 装pygimli 和它的绘图依赖 pyvista 一起进来 conda install -c conda-forge pygimli pyvista -y # 验证核心是否可导入能打印版本说明二进制没问题 python -c import pygimli as pg; print(pg.__version__)这段命令的逻辑是先隔离环境再从 conda-forge 拉预编译包最后用一次导入验证 C 核心是否加载成功。参数上python3.10是经验值太新的版本有时轮子还没跟上-c conda-forge必须显式指定否则默认频道可能找不到包。如果导入时报ImportError且提到共享库八成是环境里混进了 pip 装的同名依赖删掉环境重来比修依赖快。提示不要在同一环境里既 conda 装又 pip 装 pyGIMLi两套二进制混用是玄学崩溃的常见来源。2.2 数据格式与最小数据集的组织pyGIMLi 例子通常围绕两类数据一是电极排列的几何信息二是每个测点的视电阻率。常见做法是用一个文本文件存四列——A、B、M、N 四个电极的位置和对应的视电阻率值这是电法里最通用的交换格式。你也可以直接在脚本里用数组构造但真实项目里从文件读更贴近落地。import numpy as np import pygimli as pg from pygimli.physics import ert # 假设数据文件是空格分隔的五列a b m n rhoa # 这里用 numpy 读入再转成 pyGIMLi 的 DataContainer raw np.loadtxt(profile_01.dat) data pg.DataContainerERT() # 逐行添加四个电极位置和视电阻率 for a, b, m, n, rhoa in raw: data.createSensor([a, 0.0]) # 电极坐标二维剖面 y 固定为 0 data.createSensor([b, 0.0]) data.createSensor([m, 0.0]) data.createSensor([n, 0.0]) data[rhoa].append(rhoa) # 关键一步告诉 pyGIMLi 每个测点用的是哪四个电极编号 # 这里假设电极按添加顺序编号实际项目要按你的排列索引填 data[a] pg.Vector(range(0, len(raw))) data[b] pg.Vector(range(1, len(raw) 1)) data[m] pg.Vector(range(2, len(raw) 2)) data[n] pg.Vector(range(3, len(raw) 3)) data.checkData()逻辑说明DataContainerERT是电法专用的数据容器它要求电极坐标和测量配置分开存。createSensor负责几何data[a]这类字段负责索引。参数上checkData()会检查索引越界、重复电极、负视电阻率等问题跑反演前一定要调用否则错误会推迟到反演内部才爆排查成本翻倍。真实数据里电极编号往往不是顺序的需要你按实际排列建立映射表这一步错了后面全错。2.3 用内置例子做冒烟测试在碰自己的数据前先用 pyGIMLi 自带的例子确认整条链路是通的。库本身附带了一些合成数据生成函数可以造一个已知模型再反演回去看结果能不能对上。这是判断环境是否正常的标准动作。import pygimli as pg from pygimli.physics import ert # 生成一个合成 ERT 数据内部是一个低阻异常体 data ert.simulate( ert.createGeometricFactors, # 占位实际用下面的标准调用 )上面这段是示意正确的最小冒烟测试应该用库里的合成模型接口先建一个网格给定电阻率分布用ert.simulate正演出视电阻率再把这个视电阻率当输入做反演比较反演结果和真实模型。判断标准很简单——反演出的异常体位置和真实模型偏差在一个网格单元以内就说明环境和算法都正常。这一步过了再上真实数据心里有底。3. 正演与网格pyGIMLi 例子里的几何怎么搭才不出错3.1 网格类型的选择结构化还是非结构化pyGIMLi 支持结构化网格和非结构化网格这是新手最容易选错的地方。结构化网格是规则矩形生成快、索引简单适合地形平坦、地下结构规整的场景。非结构化网格用三角形或四面体能贴合起伏地形和复杂边界但生成和反演都更慢。我一般的原则是地形起伏小于剖面长度的百分之五用结构化超过就上非结构化否则地形会被强行拉平反演出的异常深度系统性偏移。import pygimli as pg import pygimli.meshtools as mt # 构造一个二维世界边界宽度 100 米深度 50 米 world mt.createWorld( start[0, -50], # 左下角坐标 end[100, 0], # 右上角坐标 marker1, # 区域标记反演时用来区分背景 ) # 在中间放一个低阻异常体模拟含水破碎带 body mt.createRectangle( start[40, -20], end[60, -10], marker2, # 异常体标记 ) # 合并几何生成非结构化三角形网格 geom world body mesh mt.createMesh(geom, quality33, area1.0) print(节点数:, mesh.nodeCount(), 单元数:, mesh.cellCount())逻辑说明createWorld定义计算域createRectangle定义异常体两者相加得到带标记的几何。createMesh把几何离散成网格quality33控制三角形最小角值越大网格质量越高但单元数也越多area1.0是最大单元面积直接决定分辨率。参数上area要和你关心的最小异常尺度挂钩异常体宽 20 米单元面积 1 平方米意味着横向约 20 个单元够用如果设成 100异常体只剩两三个单元反演根本分辨不出来。3.2 正演算子的建立与边界条件正演是把已知电阻率模型算成视电阻率的过程。pyGIMLi 里用ERTModelling建立正演算子它内部会组装有限元矩阵并求解。边界条件默认是混合边界对大多数剖面够用但如果你把计算域设得太小电流会从边界“漏出去”导致视电阻率整体偏低。from pygimli.physics import ert # 基于网格建立正演管理器 fop ert.ERTModelling() fop.setData(data) # 绑定测量配置 fop.setMesh(mesh) # 绑定网格 # 给每个区域赋电阻率背景 100 欧姆米异常体 10 欧姆米 rho pg.Vector(mesh.cellCount(), 100.0) for cell in mesh.cells(): if cell.marker() 2: rho[cell.id()] 10.0 # 执行正演得到合成视电阻率 rhoa_syn fop.response(rho) print(正演完成视电阻率范围:, min(rhoa_syn), max(rhoa_syn))逻辑说明setData和setMesh把测量几何和网格绑到算子上response接收一个按单元排列的电阻率向量返回视电阻率。参数上rho的长度必须等于mesh.cellCount()顺序和单元编号一致搞错顺序结果就是乱的。边界条件如果默认不满足需求可以在ERTModelling上调整但多数剖面用默认即可。正演跑通后把rhoa_syn和真实数据对比能提前发现数据里的系统性错误。3.3 网格质量对反演结果的影响网格不是越密越好。单元太密反演矩阵规模暴涨一次迭代从几秒变成几分钟而且过密网格会让反演欠定问题更严重出现虚假的小尺度异常。常见做法是先粗网格跑通看异常大体位置再局部加密验证。判断网格是否够用的方法把网格加密一倍如果反演结果的主异常位置和幅值变化小于百分之十说明当前网格已经收敛。4. 反演参数pyGIMLi 例子中最容易调崩的几个旋钮4.1 正则化强度与平滑约束反演的本质是在数据拟合和模型平滑之间找平衡。pyGIMLi 用正则化参数 lambda 控制这个平衡lambda 大则模型平滑、异常模糊lambda 小则拟合好但可能出现剧烈振荡的假结构。默认值往往不适合你的数据必须手动调。# 建立反演管理器 inv ert.ERTManager(data) inv.setMesh(mesh) # 关键参数lambda 初始值、是否用对数变换、最大迭代次数 inv.inv.setLambda(20.0) # 正则化强度先给一个中等值 inv.inv.setMaxIter(15) # 最多迭代 15 次 inv.inv.setMarquardtScheme(0.9) # 马夸特阻尼0.9 表示较保守 # 执行反演 model inv.invert() print(反演完成模型电阻率范围:, min(model), max(model))逻辑说明setLambda是核心旋钮20 是二维 ERT 的常见起点数据噪声大就往上加异常需要锐利就往下减。setMarquardtScheme控制每次迭代的步长衰减0.9 偏保守收敛慢但稳设成 0.5 会激进容易在噪声数据上发散。setMaxIter不是越大越好通常 10 到 20 次就够再往后拟合改善很小但假结构会累积。判断收敛看数据拟合残差曲线残差降到噪声水平附近就该停。4.2 数据误差估计与加权反演里每个数据点都有权重权重来自误差估计。pyGIMLi 允许你给相对误差和绝对误差这两个数直接决定反演信谁。相对误差设太小反演会拼命拟合噪声设太大真实异常被当噪声压掉。# 给数据设置误差2% 相对误差 1e-4 绝对误差 data[err] ert.estimateError( data, relativeError0.02, # 相对误差 2% absoluteUError1e-4, # 电位绝对误差 ) # 也可以手动指定一个统一的相对误差 # data[err] data[rhoa] * 0.03逻辑说明estimateError根据电位测量值估算误差relativeError反映仪器和耦合的综合不确定性野外数据一般 2% 到 5%。参数上absoluteUError是电位测量的本底噪声室内标定可以测出来。如果数据里有个别点明显是坏点不要靠调大误差去掩盖直接剔除或降权更干净。误差设完可以画误差分布图如果误差和视电阻率强相关说明模型设错了。4.3 反演结果的可视化与剖面解释反演完不画图等于没做。pyGIMLi 的可视化基于 pyvista能画二维剖面和三维切片。二维剖面是电法最常用的展示方式横轴是水平距离纵轴是深度颜色是电阻率对数。import matplotlib.pyplot as plt from pygimli.viewer.mpl import drawModel # 画反演模型剖面 fig, ax plt.subplots(figsize(10, 4)) drawModel( ax, mesh, model, logScaleTrue, # 电阻率用对数色标 cMapSpectral_r, # 色标低阻偏蓝高阻偏红 labelResistivity (Ohm-m), ) ax.set_xlabel(Distance (m)) ax.set_ylabel(Depth (m)) plt.tight_layout() plt.savefig(inversion_profile.png, dpi150)逻辑说明drawModel把单元电阻率映射到网格上着色logScaleTrue是电法惯例因为电阻率跨几个数量级。参数上cMap选发散色标能突出异常但解释时要注意色标范围范围设太宽异常会被压平。保存图片时dpi150够用于报告再高文件会很大。解释剖面时先看整体趋势再看局部异常最后结合地质资料判断异常是真实结构还是反演假象。5. 避坑与排查pyGIMLi 例子跑不通时的五条血泪经验5.1 现象反演结果全是背景值看不到任何异常原因通常是数据误差设得过大或者 lambda 太大反演直接把模型压成均匀半空间。解决方法是先把relativeError降到 1% 到 2%再把lambda从 20 降到 5 附近重新跑一次。如果还是平的检查数据本身是不是视电阻率范围太窄窄到没有异常信息。5.2 现象反演中途报矩阵奇异或内存溢出原因多是网格单元数太多或者电极索引有重复导致正演矩阵病态。解决方法是先用粗网格跑通确认索引无误后再加密同时检查data.checkData()的输出把重复电极和零长度排列剔除。内存溢出时把area调大单元数降一个数量级再试。5.3 现象反演出的异常深度和已知地质对不上原因通常是计算域深度不够电流被底部边界截断或者地形没有正确建模。解决方法是把createWorld的深度加大到剖面长度的三分之一以上有地形就用非结构化网格贴合。深度偏移是系统性的不要靠调 lambda 去修那是治标不治本。5.4 现象可视化时颜色全糊在一起看不出结构原因是色标范围设得太宽或者用了线性色标。解决方法是改用对数色标并把色标范围卡在模型电阻率的 5% 到 95% 分位数之间把极端值截掉。这样主异常会清晰很多代价是极端高阻或低阻区域被饱和显示解释时心里有数即可。5.5 现象换一台机器或换一个环境同样的脚本结果不一样原因是底层 BLAS 库或随机种子不同反演这类迭代算法对浮点运算顺序敏感。解决方法是固定随机种子并在脚本里记录 pyGIMLi 版本和依赖版本。如果结果差异在百分之几以内属于正常浮点误差差异很大说明环境里有冲突的数学库需要重建干净环境。6. 进阶技巧把 pyGIMLi 例子变成可复用的反演流水线跑通单个例子只是起点真正省时间的是把流程脚本化、参数化。我习惯把一次反演拆成三个文件一个配置文件存网格参数、反演参数和数据路径一个主脚本读配置跑正演反演一个后处理脚本专门画图和导出。这样换一条剖面只改配置不用动代码。一个具体技巧是用 pyGIMLi 的ert.ERTManager做交叉验证把数据分成奇偶两组各自反演比较两个模型的一致性。如果两组反演出的主异常位置和幅值接近说明数据质量够、反演稳定如果差异大说明数据里有噪声主导的区域解释时要保守。这个验证花不了几分钟但能避免把假异常写进报告。# 交叉验证奇偶测点分组反演 mask_odd np.arange(len(data[rhoa])) % 2 0 data_odd data.copy() data_odd[rhoa] data[rhoa][mask_odd] # 误差字段也要同步切片否则长度不匹配会报错 data_odd[err] data[err][mask_odd] inv_odd ert.ERTManager(data_odd) inv_odd.setMesh(mesh) model_odd inv_odd.invert() # 比较 model 和 model_odd 在主异常区的均值差异 # 差异小于 15% 认为反演稳定逻辑说明copy复制数据容器后rhoa和err必须同步切片只切一个会触发长度校验失败。参数上分组比例可以按测点编号奇偶也可以按空间位置分左右两半后者对横向不均匀更敏感。判断标准是主异常区电阻率均值差异15% 是经验阈值超过就要回头查数据质量。我自己的习惯是每跑完一条剖面先把反演参数、网格参数、数据版本写进一个日志文件再画图。这样三个月后回头看能准确复现当时的每一步而不是对着一张图猜当时用了什么参数。反演这东西后悔药就是日志。希望帮到你。本文还有配套的精品资源点击获取
