简介一套PyGIMLi示例脚本聚焦地质建模与反演面向地球物理、地质工程领域的初学者和研究者。压缩包内共11个Python文件大小约19KB全部为可直接运行的.py脚本并按基础、建模、反演三大主题分类组织。内容涵盖网格创建与插值、网格质量检查、有限元正演热传导问题以及直流电阻率一维反演等典型场景每个脚本对应一个独立任务如指数拟合、多项式拟合、一维块状反演等便于边运行边理解PyGIMLi的调用流程和参数设置。示例还展示了软件在几何建模、多种数值方法有限元、有限差分与反演框架上的核心用法适合作为系统学习的起始模板。已有1069人学习尤其适合希望通过实例快速掌握几何建模、数值求解、反演算法与可视化技巧的读者也可作为实际项目中的参考模板。1. 遇到 pyGIMLi 例子从一头雾水到能跑通第一行反演如果你刚接触 geophysical 反演大概率会在某个深夜搜到 pyGIMLi 这个词然后被一堆 example 脚本和术语劝退。pyGIMLi 是一个用于二维/三维电阻率、地震波、电磁法等地球物理正演与反演的开源库它的例子不是在教你怎么调用某个函数而是在展示一条完整的处理链路网格生成、正演计算、合成数据加噪、反演参数设置、结果可视化最后才是你真正要的那张地电断面图。看懂这些例子等于拿到了从数据到模型的钥匙。这篇文章把 pyGIMLi 的例子拆开讲按「例子在解决什么问题 → 怎么把例子跑起来 → 参数怎么调 → 反演结果怎么判断 → 坑在哪 → 还能怎么玩」的顺序走一遍全程用可复现的命令和代码不绕弯子。适合刚入行想快速上手的学生也适合接到勘探任务需要短期内用 pyGIMLi 出图的一线从业者。读完你能独立改例子、换数据、调反演参数并且知道哪些地方容易翻车。2. pyGIMLi 例子到底在讲什么核心对象和一条完整链路2.1 例子背后的三个核心对象网格、数据和反演框架pyGIMLi 的例子里几乎都会同时出现三个东西网格Mesh、数据容器以及反演算法入口。理解这三个对象比背函数名重要得多。网格负责描述模型空间数据容器装着你在野外测到的视电阻率或地震走时反演框架则负责把数据映射到模型上。以最常见的电阻率例子为例你会看到类似pg.meshtools.createGrid创建网格然后用pg.DataContainer承载电极编号和测量值再用pg.RobustInversion或pg.Inversion做反演。这三个对象不是孤立存在的它们通过一个「正演算子」串起来pyGIMLi 里叫pg.ERTModelling它负责根据当前模型算出理论观测数据也就是正演反演过程就是在不断调整网格单元的电导率让正演结果逼近实测数据。我个人的经验是先别急着跑完整例子把这三个对象在例子里各自干了什么标记出来然后你会发现 pyGIMLi 所有例子的骨架几乎都一样。这个共性就是它最值钱的地方学会一个等于学会了这一类。2.2 例子的标准链路从网格生成到反演出图pyGIMLi 官方仓库里的例子无论领域怎么变链路都是固定的五步第一步生成网格。大多数例子用pg.meshtools.createWorld或者直接读入一个已有的网格文件。第二步加载或仿真数据。电气法例子里常用pg.ERTManager里的simulate方法地震走时例子则用pg.SeismicManager。第三步设置反演参数包括正则化权重、Z 权重、迭代次数等。第四步调用inv.run()启动反演。第五步用pg.show可视化模型。这里要提醒的是很多例子里的「数据」是为了演示生成的合成数据并不是你野外的实测数据。合成数据的好处是真实模型已知可以用来验证反演代码对不对但代价是它太干净了野外数据里的电极接触电阻、地形起伏、坏点这些问题它都没有。所以当你把例子跑通以后一定要意识到例子只是模板替换成自己的数据才是开始。2.3 最小可运行例子一个二维电阻率反演的完整代码这里直接给一个我在学习和给新人讲解时最常用的最小例子它基于 pyGIMLi 的 ERT电阻率层析成像模块代码很短但足够说明整条链路怎么闭合。import pygimli as pg from pygimli.physics import ert # 1. 创建一个简单的两层地电模型网格 grid pg.meshtools.createWorld(start[-10, 0], end[10, -10], layers[-2], world[False]) # 2. 在网格上定义真实电阻率值上层 100 欧姆米下层 10 欧姆米 model_true pg.solver.parseMapToCellArray([[1, 100.0], [2, 10.0]], grid) grid[res] model_true # 3. 生成一个施伦贝格装置的数据容器 data ert.createData(elecspg.utils.grange(-10, 10, dx1.0), schemeNameschlumberger) # 4. 用正演模拟生成合成视电阻率数据并加 3% 高斯噪声 data ert.simulate(meshgrid, resmodel_true, datadata, noiseLevel0.03, noiseAbs1e-2) # 5. 创建反演管理器并运行反演 mgr ert.ERTManager() model mgr.invert(data, verboseFalse) # 6. 输出反演结果并绘图 print(model) pg.show(meshmodel.parent(), datamodel, showMeshTrue)这段代码的逻辑说明createWorld生成一个带有地下一层的矩形网格parseMapToCellArray给不同区域赋电阻率值createData用线性排列的 21 个电极生成标准施伦贝格测量方案simulate调用有限元正演得到视电阻率并添加噪声ERTManager.invert启动反演优化最后model对象里就是反演出的每个网格单元的电阻率。参数说明里最重要的是noiseLevel它直接影响反演难度建议从 0.01 开始逐步加大太大反演会失真太小结果会过拟合。这段代码我建议你逐行读过一遍再跑因为它是后续所有电阻率例子的骨架。如果你跑的是地震或直流电法例子把ert换成对应模块流程是一样的差异只在于正演算子和数据容器。3. 选对例子并改造成自己的数据从读例子到替换实测数据3.1 官方例子的目录结构怎么判断哪个例子适合自己pyGIMLi 的 examples 目录一般按物探方法分类常见的有ert电阻率、traveltime地震走时、tdem时间域电磁、inversion纯反演框架等。第一次接触的人容易直接挑一个名字看起来顺眼的文件然后发现里面用到了自己没装好的第三方库或者数据文件路径对不上。我的建议是先看文件名里有没有simulate或synthetic如果有说明它自带数据生成逻辑不依赖外部数据文件最适合先跑通。再挑一个自己不熟悉的装置类型比如从wenner换到schlumberger对比同一条反演链路的差异会比重复跑同类型例子收获大得多。官方例子的作用是给你提供「最小可用框架」千万别把它当成标准答案。3.2 把合成数据换成实测数据的关键改动点当你拿到自己的野外数据格式往往是protocol文件或者csv里面有电极编号、电极坐标和测量的视电阻率值。pyGIMLi 提供了ert.load方法可以直接读入多种格式但实测数据里最常见的坑是坐标单位不一致、电极编号从 1 开始而你导入后续索引从 0 开始还有缺失的电极位置信息。我一般在替换数据时会做三件事第一把数据文件转成 pyGIMLi 的DataContainer如果原数据是txt的「ABMN R」格式我写一个很小的解析脚本第二检查数据范围删除视电阻率为负值或明显偏离中位数几个数量级的点第三给数据容器附带正确的电极坐标这一步决定了网格能不能正确覆盖测线。import pygimli as pg from pygimli.physics import ert # 假设你的实测数据文件是 csv列依次为A, B, M, N, rho # 先读成 numpy 数组再构建 pyGIMLi 数据容器 import numpy as np raw np.loadtxt(field_data.csv, delimiter,, skiprows1) # 提取电极坐标这里假设电极沿 x 方向从 0 到 40 均匀分布 elecs pg.utils.grange(0, 40, nlen(set(raw[:, :4].ravel()))) # 创建容器并设置四极装置 data ert.createData(elecselecs, schemeNamecustom) # 把实测数据填入容器注意 pyGIMLi 的 data[rhoa] 对应视电阻率 data[rhoa] raw[:, 4] data[valid] np.ones(len(data), dtypebool) # 过滤掉明显异常值视电阻率小于 1 或大于 10000 欧姆米 bad (data[rhoa] 1) | (data[rhoa] 10000) data.markInvalid(bad) data.removeInvalid() # 打印数据概况检查 print(data)这段代码里有个关键的参数是schemeNamecustom当你没有标准装置方案时用它最省事pyGIMLi 不会强加电极配置约束。markInvalid和removeInvalid是处理坏点的标准动作实测数据里经常有接地电阻不好导致的飞点不处理会让反演结果出现大片假异常。注意我这里的电极坐标用grange均匀生成如果你的实测电极间距不等这个假设就错了需要从数据文件里真正读取每个电极的坐标。3.3 为什么要复用反演管理器而不是每次重建很多例子最后都会出现mgr ert.ERTManager()然后mgr.invert(data)。有的新手会问为什么不直接用pg.Inversion类非要套一个管理器。管理器的价值在于它封装了正演算子、数据变换、误差加权、默认边界条件和可视化等一堆反演外围逻辑你只改数据就可以迭代反演不需要重复配置。实际使用中我更推荐这种模式数据预处理一次生成多个测线数据文件然后循环喂给同一个ERTManager每次invert会复用上一次的网格初始猜测对多测线联合成像特别有用。但注意复用管理器有一个副作用上一次反演得到的模型会成为下一次的初始模型。如果不是连续测线记得调用mgr.reset()或者在创建管理器时指定mgr ert.ERTManager(verboseTrue)避免上一测线的残差影响本次结果。4. 反演参数逐个拆解从默认值开始调而不是冲动改4.1 正则化参数 lambda 和 Z 权重控制平滑度的两个旋钮pyGIMLi 的大部分例子默认使用二阶平滑约束正则化通过pg.meshtools生成的网格配合盒状约束使反演结果在地层边界处仍能保持清晰。反演控制器里常用的参数是lam正则化权重和zWeight垂直与水平平滑比例。如果你直接用mgr.invert(data)库会启动一个自动选择 lambda 的机制从较大值开始逐渐降低直到拟合残差达到目标。这看起来省心但在实际数据上自动选项往往偏保守结果过于平滑分辨率不够。我的做法是先把自动反演跑一遍记下它最终选择的 lambda 值然后把它乘以 0.3 到 0.5 作为手动初始值再看模型细节是否增多且不发散。zWeight 的作用是调整垂直和水平方向平滑的权重默认值在 pyGIMLi 里通常是 1意味着各向同性平滑。对测线方向拉长的二维结构我会把 zWeight 设为 0.5 到 0.8让垂直方向有更多变化空间。inv mgr.inversion # 获取反演对象 inv.setLambda(20.0) # 手动指定正则化权重而不是自动选择 inv.setZWeight(0.8) # 垂直方向平滑比水平方向略弱 # 设置反演停止条件相对残差变化小于 0.01 或达到 15 次迭代 inv.setMaxIter(15) inv.setChiLimit(1.01) mgr.invert(data, startModelpg.meshtools.createGrid(...))参数说明setLambda(20.0)在我的同类型数据集上是从自动值 54 调下来的模型细节明显增加但仍在稳定范围setChiLimit控制卡方目标1.01 对含 3% 噪声的数据比较合理如果噪声更大可以放宽到 1.05。这里要注意startModel必须是一个网格或者与待求模型长度一致的数组否则反演会直接报错。4.2 数据误差模型比正则化更值得先调的东西pyGIMLi 默认认为数据是高斯独立同分布的误差由绝对误差和相对误差组成。在 ERT 模块里simulate函数生成合成数据时加的noiseLevel0.03和noiseAbs1e-2到了反演阶段pyGIMLi 会默认用数据本身的中位数来估计误差比例。这个假设在大部分野外数据上是脆弱的。如果你用的是实测数据强烈建议在反演前显式设置数据误差。pyGIMLi 中data[err]是一个可选字段如果没有反演内部会用经验公式生成。我踩过最直观的坑是不同电极排列的测量数据误差差异很大比如偶极-偶极排列在长极距时信噪比低如果统一用一个固定相对误差长极距的噪点会被放大成浅部假异常。正确做法是给每个测量值设定独立误差一般做法是data.setError(np.abs(data[rhoa]) * 0.02 0.005)或者更稳妥地从重复测量数据里估算标准差。pyGIMLi 反演框架在加权最小二乘里会把误差作为分母误差越小权重越高所以低误差数据会对反演产生主导作用。如果你发现自己的反演结果被某个数据点控制去查一下那个点的误差是否过小。4.3 网格密度与地形两个容易被忽视的边界pyGIMLi 例子里常用createWorld生成规则网格边界处网格无限向外延伸避免边界反射。当你把例子改成实测数据时必须重新检查网格是否覆盖所有电极坐标并考虑地形。pyGIMLi 支持带地形的网格常见做法是用pg.meshtools.createSurveyMesh根据电极坐标生成贴合地表的网格。网格太疏会导致分辨率不够太密会导致反演矩阵巨大而且浅层容易出现锯齿。我的经验是在电极间距的 1/5 到 1/10 之间设置网格细化级别。你不需要手工建网格用pg.meshtools.createSurveyMesh( data, quality34)会自动在电极附近加密网格。地形处理方面如果测线高差超过电极间距的 5%必须把电极坐标的 z 值传入网格生成函数否则反演模型会把人造地形起伏当成浅部电性异常。5. pyGIMLi 例子避坑我在实际项目里踩过的五个真实问题5.1 现象反演迭代发散每次迭代残差反而增大原因最常见的是数据里有负值或零值。pyGIMLi 反演通常采用对数电阻率参数零或负值取对数直接报错即使没报错也会导致正演矩阵奇异。另一种可能是起始模型给了一个极端值比如startModel给了 1e-6 欧姆米正演计算浮点溢出。解决先做数据清洗data[rhoa] 0的点一律删掉再看startModel用中位数初始化。代码里顺手加一行data.removeInvalid()能避免 80% 的发散。5.2 现象反演结果出现带状条纹沿测线方向规则排列原因电极位置坐标有误最常见是电极间距不等但用了均匀网格或者 x 和 z 坐标填反。另一个原因是测量数据存在系统偏差比如电极接触电阻大导致所有相邻极距的数据整体偏高反演不得不以条带状模型补偿。解决先用pg.show(data)查看数据的伪剖面图如果伪剖面里有明显的横向突变优先怀疑数据质量再检查网格与电极坐标是否严格匹配使用mgr.model与电极位置比对。5.3 现象反演结果与真实地质分层不符深部电阻率平均化严重原因这不是 pyGIMLi 的问题而是观测装置对深部的灵敏度不足。偶极-偶极装置比温纳装置有更好的横向分辨率但深部信号弱且误差大。如果沿用例子里的施伦贝格装置深部约束本来就弱。解决换装置类型重新采集数据或者在反演中加入先验约束。pyGIMLi 支持pg.Region约束某个区域电阻率固定比如你知道表层是冲积层就把浅层区域约束在一定范围内减轻深部平均化。5.4 现象运行反演时内存爆炸电脑直接卡死原因网格太密或者反演求解器默认用了直接法。三维反演时尤其明显二维问题一般还好。在高分辨率地形起伏大的测线上网格数量上了百万直接法求一次正演就吃掉几十 G 内存。解决降低网格质量参数从quality34调到quality30或者限制网格最大面积如果 pyGIMLi 版本支持指定求解器为迭代法常见做法是创建pg.solver.Solver并设置solvercg然后传给反演管理器。不过这会让反演时间变长需要权衡。5.5 现象pg.show显示模型时颜色范围永远被少量极端值霸占断面图看不出结构原因这个是可视化常见误用pyGIMLi 默认的色标范围是数据的最小最大值非鲁棒。解决手动指定cMin和cMax一般用数据的 5% 和 95% 分位数或者直接先对模型取对数再显示。pg.show(mesh, model, cMin10, cMax200)就比自动范围实用得多。这个小技巧在实际报告里非常重要审阅者第一眼看到的就是图颜色范围不对等于模型白反演了。6. 把例子进化成自己的研究工具批量反演与结果质量评估6.1 用脚本批量处理多条测线从手动跑例子到流程自动化你手里的野外数据一般不止一条测线。手动改文件名跑例子太慢我通常会写一个循环批量读入测线、反演、保存结果图。这一步能最大限度发挥 pyGIMLi 的价值。import glob from pygimli.physics import ert # 获取所有测线数据文件 files glob.glob(survey_*.csv) for f in files: data ert.load(f) # 假定是 pyGIMLi 支持格式 mgr ert.ERTManager() model mgr.invert(data, verboseFalse) # 保存反演结果模型为 vtk 格式方便后续在 paraview 里查看 pg.save(f.replace(.csv, _model.vtk), model) # 保存反演拟合图 ax pg.show(mgr.fop.mesh(), model, cMin10, cMax200) pg.show(data, axax, cMin10, cMax200)逻辑说明glob匹配所有测线文件ert.load自动识别常见格式然后每条测线单独创建ERTManager避免初始模型串扰。保存为 vtk 是因为后续我想在专业可视化软件里剖切联合显示pyGIMLi 自带的绘图适合快速验证但不好排版。运行这段之前务必在每条测线上跑一次单条数据确认没有坏点。自动化里最容易出问题的就是某条测线数据有问题导致反演中断所以在循环外层加try...except记录失败文件避免整个批处理崩掉。6.2 用 L-curve 和模型粗糙度判断反演结果可信度pyGIMLi 反演完成后除了看图还需要量化说明结果可信度。mgr.invert的历史记录里保存了每次迭代的残差和模型粗糙度我常用mgr.inversion.chi2看最终卡方值如果它远大于 1说明数据误差估计偏低或者模型不满足数据如果小于 0.8可能过拟合了。一种更直观的验证方法是做「正演验证」把反演得到的模型作为真实模型用同样的测量方案重新生成模拟数据再对比模拟数据和实测数据的拟合程度也就是交叉验证思路。我经常这么做用一套网格做反演用另一套更密或更疏的网格做正演看结果是否稳健。如果模型对网格变化特别敏感说明它主要由正则化控制而不是由数据驱动解释时要格外小心。import pygimli as pg from pygimli.physics import ert # 已知反演模型 model_res 和对应的网格 mesh mgr_check ert.ERTManager() new_data mgr_check.simulate(meshmesh, resmodel_res, schemedata.scheme()) # 对比模拟数据与实测数据的归一化 RMS rms pg.norm((new_data[rhoa] - data[rhoa]) / data[rhoa]) / pg.norm(data[rhoa]) print(Normalized RMS:, rms)这段代码的用途不是重新反演而是单纯做正演模拟来验证模型。如果模拟数据和实测数据的 RMS 在 0.05 以内说明这个模型对当前数据是合理的但如果 RMS 很大即使反演收敛了也要怀疑网格或误差设置有问题。6.3 我最常用的一个技巧反演前把数据做对数变换pyGIMLi 内部在反演时默认对电阻率做对数变换但很多例子在预处理阶段不会主动告诉你。野外视电阻率数据跨度经常从个位数到几千欧姆米如果不做对数变换高阻值的数据点会在目标函数中占据主导低阻地区的细节会被掩盖。我有个习惯在把实测数据塞进反演之前先把rhoa取对数并对对数数据设置误差这样反演更稳定。pyGIMLi 里可以通过data.log()复制一份对数数据容器或者直接在读取后设置data[rhoa] np.log(data[rhoa])然后反演管理器内部会自动检测。最终显示的时候你再指数还原注意pg.show的色标标签标成欧姆米即可。这个技巧对电阻率跨几个数量级的项目帮助很大。我印象最深的是一次在第四系覆盖区做工程勘察表层干砂电阻率上千下部含水黏土只有几十不做对数变换反演含水层边界基本糊成一片做了对数变换后分层清晰多了。反演不仅仅是调参数数据变换本身也是反演的一部分这个观点在 pyGIMLi 的例子里体现得不明显需要你自己在实际项目中体会。希望我讲的这些能让你少走一点弯路跑通例子只是一小步用好参数和数据判断才是在这个方向投入时间的真正回报。本文还有配套的精品资源点击获取