简介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(start0.0, stop100.0, n101) 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( lam20.0, maxIter10, verboseTrue, robustDataTrue, zWeight0.5 )lam20.0是很多标准例子的起始值如果反演出来的模型太碎、异常体边缘锯齿严重就调大lam如果数据拟合差、残差始终降不下来就适当减小。robustDataTrue表示用 L1 范数处理数据中的离群点实测数据有异常跳变时强烈建议开。zWeight0.5是垂直方向平滑权重当地质体倾向水平层状时经常设成 0.5 以下。还有个被很多人忽略的点ERTManager会自动做数据误差估计。默认按 3% 相对误差算如果实际仪器精度达不到反演结果的信噪比会虚高。我会在数据加载后手动重置误差见后面避坑章节。2.3 网格质量对反演结果的决定性影响网格是反演的容器。容器不对再好的算法也白搭。pyGIMLi 里网格生成用的是无结构三角网格相比规则矩形网格它能把地形和异常体边界贴合得更自然。但无结构网格也意味着你要对网格密度分布有概念电极附近的网格必须加密远离测线的区域可以放粗。常见做法是先用createWorld生成一个大范围背景网格再在电极位置用createPolygon或者createCircle嵌入局部加密区。下面这段代码演示了如何把加密区域和背景网格合并# 创建一个圆形的局部加密区域 circle mt.createCircle( pos[30.0, -5.0], # 圆心坐标 radius3.0, # 半径 marker5, # 区域标记 area0.2 # 该区域最大单元格面积 ) # 合并背景网格和加密区域 merged mt.merge(geom, circle) # 生成最终的计算网格 mesh mt.createMesh(merged, quality34.0, area2.0)area参数控制加密区的网格单元面积上限值越小网格越密。这里circle的area0.2意味着加密区内每个三角单元面积不超过 0.2 平方米而背景网格area2.0就疏松得多。quality34.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, schemeforward_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( lam20.0, maxIter15, startModelpg.Vector(1, 100.0), # 初始模型均匀100欧姆米 lamTolerance0.1, # 正则化参数下降阈值 fitTolerance0.01, # 数据拟合相对变化阈值 verboseTrue ) # 查看反演结果摘要 print(mgr.inv.invModel) print(均方根残差:, mgr.inv.absRMS)参数含义拆开看startModelpg.Vector(1, 100.0)构造一个只有 1 个元素的向量初始值 100 欧姆米pyGIMLi 会把它自动扩展到所有网格单元。lamTolerance0.1表示当正则化参数的变化率小于 10% 时认为达到收敛fitTolerance0.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 迭代被这个异常值持续拉扯。另一个可能原因是初始模型远离真实模型迭代在局部极小值附近抖动。解决打开robustDataTrue这个参数会切换为 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, coveragesens_norm, labelResistivity (Ohm m))coverage参数控制模型显示时被灵敏度阈值遮罩的程度。coverage0.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( lam20.0, maxIter12, robustDataTrue, verboseFalse ) # 保存反演结果和灵敏度 mgr.saveResults(f.replace(.dat, _inverted.png)) mgr.inv.save(f.replace(.dat, _invModel.npy)) print(f第{idx1}/{len(dat_files)}条测线完成残差: {mgr.inv.absRMS:.2f}%)saveResults会生成带色标和地形显示的反演剖面图save把模型向量直接存成 numpy 格式方便后续批量读取做统计对比。一个建议批处理时verboseFalse关掉逐次迭代输出只在每个文件结尾打印残差摘要日志会干净很多。absRMS的单位是百分比表示数据拟合的均方根相对误差这个值也是你判断每条测线数据质量最直接的指标。5.3 与地质解释衔接剖面图叠加钻孔信息反演剖面最终要能跟地质分层对上话。我经常在 pyGIMLi 生成的图上叠加钻孔岩性分界点用来校验反演边界和实际地质界线是否一致。import matplotlib.pyplot as plt # 假设一个钻孔在 x45 米处岩性分界在 -12 米 ax pg.show(mesh, model, labelResistivity (Ohm m)) ax[0].plot([45.0, 45.0], [0.0, -25.0], k--, linewidth1.5, labelborehole) ax[0].plot([45.0], [-12.0], rv, markersize10, labelinterface) ax[0].legend() plt.show()坐标对应关系要注意pyGIMLi 的 y 轴向下为正所以深度 -12 米在垂直方向上正好显示在下方。钻孔位置如果落在测线范围之外plot还是会把线画出来但已经超出网格覆盖区需要在图上留意别误解。叠加图的意义在于快速验证反演结果的可靠性为最终地质解释提供一层独立证据。从那以后我每次跑完测线数据都会强制走一遍「正演验证 → 反演 → 灵敏度评估 → 钻孔对比」这个流程少选条反演直接交图的情况就很少发生了。希望帮到你。本文还有配套的精品资源点击获取 pyGIMLi 电阻率层析成像全流程解析：正演、反演与灵敏度评估
简介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(start0.0, stop100.0, n101) 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( lam20.0, maxIter10, verboseTrue, robustDataTrue, zWeight0.5 )lam20.0是很多标准例子的起始值如果反演出来的模型太碎、异常体边缘锯齿严重就调大lam如果数据拟合差、残差始终降不下来就适当减小。robustDataTrue表示用 L1 范数处理数据中的离群点实测数据有异常跳变时强烈建议开。zWeight0.5是垂直方向平滑权重当地质体倾向水平层状时经常设成 0.5 以下。还有个被很多人忽略的点ERTManager会自动做数据误差估计。默认按 3% 相对误差算如果实际仪器精度达不到反演结果的信噪比会虚高。我会在数据加载后手动重置误差见后面避坑章节。2.3 网格质量对反演结果的决定性影响网格是反演的容器。容器不对再好的算法也白搭。pyGIMLi 里网格生成用的是无结构三角网格相比规则矩形网格它能把地形和异常体边界贴合得更自然。但无结构网格也意味着你要对网格密度分布有概念电极附近的网格必须加密远离测线的区域可以放粗。常见做法是先用createWorld生成一个大范围背景网格再在电极位置用createPolygon或者createCircle嵌入局部加密区。下面这段代码演示了如何把加密区域和背景网格合并# 创建一个圆形的局部加密区域 circle mt.createCircle( pos[30.0, -5.0], # 圆心坐标 radius3.0, # 半径 marker5, # 区域标记 area0.2 # 该区域最大单元格面积 ) # 合并背景网格和加密区域 merged mt.merge(geom, circle) # 生成最终的计算网格 mesh mt.createMesh(merged, quality34.0, area2.0)area参数控制加密区的网格单元面积上限值越小网格越密。这里circle的area0.2意味着加密区内每个三角单元面积不超过 0.2 平方米而背景网格area2.0就疏松得多。quality34.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, schemeforward_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( lam20.0, maxIter15, startModelpg.Vector(1, 100.0), # 初始模型均匀100欧姆米 lamTolerance0.1, # 正则化参数下降阈值 fitTolerance0.01, # 数据拟合相对变化阈值 verboseTrue ) # 查看反演结果摘要 print(mgr.inv.invModel) print(均方根残差:, mgr.inv.absRMS)参数含义拆开看startModelpg.Vector(1, 100.0)构造一个只有 1 个元素的向量初始值 100 欧姆米pyGIMLi 会把它自动扩展到所有网格单元。lamTolerance0.1表示当正则化参数的变化率小于 10% 时认为达到收敛fitTolerance0.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 迭代被这个异常值持续拉扯。另一个可能原因是初始模型远离真实模型迭代在局部极小值附近抖动。解决打开robustDataTrue这个参数会切换为 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, coveragesens_norm, labelResistivity (Ohm m))coverage参数控制模型显示时被灵敏度阈值遮罩的程度。coverage0.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( lam20.0, maxIter12, robustDataTrue, verboseFalse ) # 保存反演结果和灵敏度 mgr.saveResults(f.replace(.dat, _inverted.png)) mgr.inv.save(f.replace(.dat, _invModel.npy)) print(f第{idx1}/{len(dat_files)}条测线完成残差: {mgr.inv.absRMS:.2f}%)saveResults会生成带色标和地形显示的反演剖面图save把模型向量直接存成 numpy 格式方便后续批量读取做统计对比。一个建议批处理时verboseFalse关掉逐次迭代输出只在每个文件结尾打印残差摘要日志会干净很多。absRMS的单位是百分比表示数据拟合的均方根相对误差这个值也是你判断每条测线数据质量最直接的指标。5.3 与地质解释衔接剖面图叠加钻孔信息反演剖面最终要能跟地质分层对上话。我经常在 pyGIMLi 生成的图上叠加钻孔岩性分界点用来校验反演边界和实际地质界线是否一致。import matplotlib.pyplot as plt # 假设一个钻孔在 x45 米处岩性分界在 -12 米 ax pg.show(mesh, model, labelResistivity (Ohm m)) ax[0].plot([45.0, 45.0], [0.0, -25.0], k--, linewidth1.5, labelborehole) ax[0].plot([45.0], [-12.0], rv, markersize10, labelinterface) ax[0].legend() plt.show()坐标对应关系要注意pyGIMLi 的 y 轴向下为正所以深度 -12 米在垂直方向上正好显示在下方。钻孔位置如果落在测线范围之外plot还是会把线画出来但已经超出网格覆盖区需要在图上留意别误解。叠加图的意义在于快速验证反演结果的可靠性为最终地质解释提供一层独立证据。从那以后我每次跑完测线数据都会强制走一遍「正演验证 → 反演 → 灵敏度评估 → 钻孔对比」这个流程少选条反演直接交图的情况就很少发生了。希望帮到你。本文还有配套的精品资源点击获取