
简介这是一份基于开源软件PyGIMLi的Python教程示例包面向地球物理、地质工程领域的科研人员、工程师及入门学习者旨在帮助读者快速掌握地质建模与反演的基本流程。压缩包内共有11个纯Python脚本压缩后体积仅19KB却完整覆盖了基础、建模、反演三大主题基础模块演示模型创建、网格插值与质量检查建模模块基于有限元法实现二维直流电法勘探与热传导问题的正演模拟反演模块则包含光滑约束与块体约束的一维直流电阻率数据反演案例并配有结果可视化脚本方便对比不同反演策略的效果。这些示例均源自官方教程的自动整理版本目录结构化强每个脚本可直接运行读者可边运行边观察输出从而将抽象的理论概念与具体代码实现一一对应。目前已有1069人参与学习下载无论初学者还是一线研究人员都能通过这套轻量级示例快速上手PyGIMLi的几何建模、数值求解与反演算法为后续实际数据处理和科研工作打下扎实基础。1. 从 pyGIMLi 例子开始先弄懂这个库到底解决什么问题一次野外电阻率测深数据拿回来丢进常规成图软件里只能画视电阻率等值线可你真正想知道的是地下几十米内的电阻率分布模型。想过把测深曲线转成深度剖面、把异常体边界固定下来却发现手里的工具做不到。这个问题我在某工地的水文勘察里撞过多次最后是靠一套 Python 地球物理反演框架解决的它叫 pyGIMLi。pyGIMLi 的定位很直接用统一的数据结构和反演引擎把直流电法、电磁法、地震、重力等正反演问题封装成几个类让你在 Jupyter Notebook 里就能从原始野外观测数据走到带深度刻度的模型剖面。本文不打算复述官方手册而是按“最小可复现例子”整理成一条能直接搬走的路径适合水文地质、工程勘察、环境评估、考古物探方向的一线从业者。2. 环境与数据准备装好 pyGIMLi把第一份 ERT 数据读进来2.1 安装与环境pip 之外的三个注意点pyGIMLi 的安装看起来只是一条命令但实际环境里翻车大多出在依赖和 Python 版本上。最常见的做法是直接 pip 安装命令就一行pip install pygimli这条命令会把核心库和 pygimli.physics 子模块一起装上。装完后验证是否可用比急着跑例子更重要我一般先执行一段最小导入import pygimli as pg from pygimli.physics import ert print(pg.__version__)如果上面这段能输出版本号说明核心库和 ERT 模块都正常。装环境时有三个注意点值得单独说。第一pyGIMLi 依赖 numpy、scipy、matplotlib如果这些基础库版本太旧导入时会直接报符号错误建议先升级它们再装 pyGIMLi。第二低版本 Python 在某些操作系统上会遇到编译问题我在某台只能装 Python 3.9 的旧机器上编译失败过换到 conda 环境装预编译包就稳定了所以遇到 pip 安装失败时优先用 conda 建一个干净环境。第三如果你之前装过旧版 pyGIMLi重新安装前最好先卸载干净我见过一次新旧版本混装导致pg.createGrid行为异常的情况排查了很久才意识到是路径冲突。环境这块的坑不是每个项目都踩但踩一次代价很高建议第一步就做导入验证。2.2 数据容器用 ert.load 读 dat 文件并检查数据结构pyGIMLi 对野外数据的抽象集中在一个叫 DataContainer 的对象里ERT 数据就是它的一个实例。读一个标准的 dat 文件只需要两行import pygimli as pg from pygimli.physics import ert data ert.load(field_data.dat) print(data)执行print(data)后会列出传感器数量、四极测量数量以及数据里存了哪些列。一个典型的二维电阻率测深 dat 文件内部大致分两块先是传感器位置块记录每个电极的 x、y 坐标然后是数据块按“a b m n r”的字段顺序记录每次四极测量的正负电流电极、正负电位电极和电阻值。读进来之后a、b、m、n 是电极索引从 0 开始计数r 是对应的测量值。这一步非常重要因为后边反演网格的边界、电极位置映射全依赖这份索引关系索引错一位反演出来的异常位置能偏出好几米。读进来后别急着反演先检查数据质量。一个我做过的快速检查是看测量值范围import numpy as np r np.array(data[r]) print(测量值范围, r.min(), r.max()) print(负值个数, (r 0).sum())电阻率测量值理论上都应该是正数出现负值或零值通常意味着电极接反、接触不良或数据文件里混入了坏点。pyGIMLi 的ert.load不会自动剔除这些点它把剔除的权力留给使用者所以这类检查能帮你把明显坏数据挡在反演之前。别小看这一步我遇到过一份 500 组数据里混了 20 个负值的数据默认反演后模型出现了一个孤立的极高阻异常最后排查才发现是数据问题不是地下真有高阻体。2.3 创建你自己的合成数据从电极坐标到四极测量没有现成 dat 文件或者想验证流程是否可靠时可以先用 pyGIMLi 生成一套合成数据。这个做法贯穿整个学习过程因为你知道真实模型是什么就能检验反演流程有没有把模型找回来。下面这段代码建立一条测线、设定电极间隔并生成 Wenner 装置的测量序列import numpy as np import pygimli as pg from pygimli.physics import ert # 布置 21 个电极间距 1 米沿地表排列 x np.arange(0, 21, 1.0) y np.zeros(len(x)) sensors pg.PosVector() for xi, yi in zip(x, y): sensors.push_back(pg.Pos(xi, yi)) # 建立一个二维网格横向 0~20 米纵向 0~-10 米 mesh pg.createGrid(np.linspace(0, 20, 41), np.linspace(0, -10, 21)) # 给每个网格单元赋电阻率浅层 100 欧姆米深层 10 欧姆米 res pg.RVector(mesh.cellCount(), 100.0) for cell in mesh.cells(): if cell.center().y() -5.0: res[cell.id()] 10.0 # 由模型正演得到视电阻率数据 mgr ert.ERTManager() data mgr.simulate(meshmesh, resres, schemewenner, electrodePositionssensors) print(data)这里解释几个关键参数。pg.createGrid的两个数组分别是 x 方向和 y 方向的网格边界y 用负值代表深度方向这是 pyGIMLi 的习惯很多新人在这里把正负搞反结果模型跑到地表以上去了。pg.RVector是 pyGIMLi 内部的向量类型传入单元数和一个初始值生成与网格单元一一对应的电阻率数组。simulate方法的schemewenner表示用温纳装置生成测量序列electrodePositions传定位好的电极列表。整个流程相当于你先构造了一个已知的地下模型再让 pyGIMLi 替你“做一次野外测量”得到的data就可以当作后续反演的输入了。这比拿实测数据练手更安心因为答案就在你手里。3. 跑通第一次二维 ERT 反演参数选型与结果解读3.1 最小反演命令ERTManager 的默认参数能直接出结果数据在手下一步就是反演。pyGIMLi 把反演过程封装得很干净最小反演命令只有三行import pygimli as pg from pygimli.physics import ert data ert.load(field_data.dat) # 或直接用上一节生成的 data mgr ert.ERTManager(data) model mgr.invert(maxIter8) print(反演完成模型单元数, len(model))这一段代码能跑通是因为ERTManager内部自动完成了网格生成、初始模型估计、正演算子和反演迭代的控制。model返回的是一个向量长度与反演网格的单元数一致每个元素是该单元的电阻率。跑完后调用mgr.showResult()会在当前绘图窗口画出一条带深度刻度的电阻率剖面。需要注意默认网格和默认参数只是给了一个“能出结果”的下限。我第一次用默认参数跑一块已知模型时得到的剖面形状大致对但异常体边界模糊、幅度被压低这说明默认的平滑约束起了作用但也把真实对比度磨掉了。所以最小反演命令适合用来确认整条流程通不通真正要用于解释的结果必须进入参数调整这一步。在反演圈里有句俗话反演出来的不是真实模型而是满足约束条件的一个平滑版本。这句话不是泼冷水而是提醒你参数选型决定了结果往哪个方向偏。3.2 反演参数选型lam、zWeight、maxIter、blocky 到底调什么ERT 反演本质上是在拟合实测数据和保持模型平滑之间找一个平衡点ERTManager.invert()提供的那几个参数就是这个平衡点的控制器。下面是我常用的参数表按重要程度排参数默认值作用调整方向lam自动估计正则化权重越大模型越平滑过拟合时增大欠拟合时减小zWeight1.0垂直方向与水平方向平滑的权重比地层水平层状时设 0.5~0.8maxIter默认 10 次左右迭代次数上限收敛慢时调大但别超过 30blockyFalse使用 L1 范数模型约束异常体边界锐利时设 TruerobustFalse对数据中的离群点降权数据干扰大时设 True实际调用时这样传参import pygimli as pg from pygimli.physics import ert mgr ert.ERTManager(data) model mgr.invert(lam20, zWeight0.8, maxIter15, blockyFalse)参数不是越多越好我建议先动zWeight和blocky这两个因为它们直接对应地下水文或工程勘察里最常见的地质假设地层近似水平层状还是存在倾斜的、边界清晰的岩体接触带。某次在某岩溶区做探测我默认平滑约束反演出的溶洞边界糊成一团把blocky设为 True 之后溶洞顶底界面才显形。而lam这种参数调节起来有点像玄学不同数据规模、不同电极数量下最优值差异很大我的习惯是先固定其他参数从默认值开始按 10、20、50、100 几档试看模型形态发生突变的那一档就是敏感区。3.3 结果解读视电阻率伪截面与真实模型的差别在哪反演完成不等于拿到了最终答案读图这一步才是决定解释方向的关口。mgr.showResult()画出来的是反演模型剖面横坐标是测线距离纵坐标是深度颜色代表电阻率。但很多人容易忽略这个剖面里的网格单元尺寸并不是均匀的pyGIMLi 会在电极附近自动加密所以近地表的分辨率高深部单元大且平滑读图时不要对深部的细碎异常过度解释那很可能是网格粗化带来的伪影。与视电阻率伪截面相比反演剖面的优势在于它近似还原了真实模型的深度位置但代价是引入了平滑约束和网格剖分带来的不确定性。判断结果靠不靠谱一个直接办法是看拟合差pyGIMLi 反演过程中会在日志里打印 RMS 值如果 RMS 在 3% 以内说明正演响应和实测数据拟合得不错如果超过 10%就说明模型结构和数据结构存在系统性问题先回头检查数据质量而不是继续调参数。我一般还会把mgr.chi2()打出来看一眼这个值接近 1 说明数据误差估计合理远大于 1 说明模型太简单或误差被低估。4. 网格、约束与数据误差把反演结果做可靠的关键4.1 网格生成用 createGrid 控制反演网格密度反演网格是 pyGIMLi 里最终决定计算量和精度的因素很多人在这上面栽跟头。pyGIMLi 能自动生成网格但自动生成的网格不一定适合所有数据尤其是测线长、电极密的情况默认网格会在远离电极的区域堆太多单元白白消耗内存。手动控制网格用pg.createGrid下面是一个典型做法import numpy as np import pygimli as pg x np.linspace(0, 20, 51) y np.linspace(0, -12, 31) mesh pg.createGrid(x, y) print(网格单元数, mesh.cellCount())这个网格横向 51 个节点、纵向 31 个节点单元数大约在 1500 到 2000 之间跑二维 ERT 反演完全够用。原则是电极附近网格要密没有电极覆盖的深部和侧向区域要疏。pg.createGrid生成的网格是均匀的做不到局部加密但它足够演示流程。实际项目中更常见的做法是先生成一个均匀粗网格反演稳定后再细化网格反演一次用粗网格的结果做初始模型。这里有个经验值网格最浅层单元大小不要超过最小电极间距的 1/2否则电极附近的电位梯度拟合不准。4.2 数据误差模型为什么 relativeError 要设成 0.03~0.05反演算法需要知道测量数据的可信程度误差模型设得不对反演就会为了拟合一个 5% 的噪声点而把模型扯变形。pyGIMLi 里最常用的误差设置是相对误差加绝对误差通过ERTManager构造参数传入from pygimli.physics import ert mgr ert.ERTManager(data, relativeError0.05, absoluteError0.001)relativeError0.05表示认为每条测量值都带 5% 的相对噪声absoluteError表示有一个最小绝对噪声下限防止测量值很小时相对噪声变成无穷大。常规野外数据设 3% 到 5% 是合理的起点高质量实验室数据可以压低到 1%。如果你不做任何设置pyGIMLi 会用一个默认估计但默认值不一定贴合你的仪器和数据质量。数据误差的影响在实践中很容易被忽视。某次我用一组高精度数据反演RMS 一直压不下来后来发现是绝对误差设成了 0导致小电阻率测量值被算法当成了绝对精确反演结果在高阻区出现锯齿状抖动。把这个参数加上后结果立刻稳定。误差模型这块本质上是反演的“信任度声明”交给黑匣子自动处理省心但想控制结果质量必须自己动手设。4.3 结果导出把模型和网格交给其他工具反演结果不能只留在内存里导出去给制图软件或三维可视化工具是常态。pyGIMLi 提供了一套保存机制几行命令就能把模型、网格、响应数据落盘import pygimli as pg from pygimli.physics import ert mgr ert.ERTManager(data) model mgr.invert(maxIter10) mgr.saveResult(outdirresult_dir)执行后会在result_dir目录下保存模型文件、网格文件、正演响应和反演日志模型矢量可以继续加载到 numpy 里做二次计算import numpy as np from pygimli.physics import ert model np.asarray(mgr.model) # 模型电阻率向量 data_response np.asarray(mgr.response) # 模型正演响应 r np.asarray(data[r]) misfit (r - data_response) / r # 逐点相对误差 print(平均相对误差, np.mean(np.abs(misfit)))这套导出逻辑的价值在于它可以接进你自己的报告流程先把反演模型导出成网格和向量再用脚本生成剖面图或转成其他格式给制图软件。我习惯把相对误差逐点算出来画成误差分布图贴在剖面图旁边给看不懂反演机制的同事一个直观的信任依据。5. pyGIMLi 反演常见问题与避坑清单3 个让你翻车的真实场景5.1 现象一反演出来的模型看着平滑但对不上实测曲线你跑了十来次反演RMS 一直在 10% 附近降不下去模型剖面却平滑得很漂亮。这种情况单独看剖面容易误判为成功实际上正演响应和实测数据严重不吻合剖面是平滑约束硬推出来的不是数据驱动的。原因多半是数据里有系统误差或者误差模型设得不够反演算法只能通过拉平模型来妥协。解决分两步先检查实测数据与正演响应逐点比较找出偏离最大的测量点确认是不是野值再确认你的相对误差设置是否符合实际仪器精度。如果数据本身没问题就把robustTrue打开让它对离群点自动降权往往能压住 RMS。5.2 现象二一调网格就内存爆炸杀进程才罢休把网格从 2000 个单元加到 8000 个单元后程序跑到一半内存被吃光机器卡到鼠标都动不了。二维反演的计算量随网格单元数增长但实际翻车通常不是因为单元总数而是网格里混进了大量细碎的畸形单元稀疏矩阵求解时填充项爆炸。解决的办法是不要盲目加密整个网格。我现在的做法是先用粗网格跑通反演确认模型大致形态需要精细结果时只对电极附近的浅层区域加密远离电极的深部和侧向区域保持粗网格。pyGIMLi 的自动网格已经考虑了电极加密你手动传入反演的网格时也要保留这个逻辑。另一个经验是先把测线范围缩短到你要解释的目标区段再决定网格边界别让一条几百米的测线带着一台 16G 内存的机器硬跑那纯粹是和自己过不去。5.3 现象三异常体位置对但边界扭成波浪状反演出的异常体中心位置准但边界像波浪一样起伏与实际地质条件完全不符。这个问题在埋深大、电极距大的目标上尤其常见。原因有两层一是反演网格在深度方向太粗模型对边界变化的表达能力不足二是深层数据覆盖角度有限横向约束弱模型在深层可以自由波动而不显著影响拟合差。解决思路是先查网格在目标深度处的垂向单元尺寸如果比异常体的估计厚度大好几倍先细化这一带的网格。如果网格已经够密问题就在约束权重把zWeight从 1.0 调到 0.5 左右强化水平方向平滑波浪状边界会明显收敛。要提醒的是波浪边界不一定全是网格问题也可能是地表地形起伏没有正确叠进网格pyGIMLi 支持高程数据的叠加地形没进网格时反演剖面会把浅层异常拉歪。每次做地形起伏大的测线时检查一下网格节点的高程是否正确能省下很多无效调参时间。6. 用合成数据做盲测给反演流程上个保险6.1 合成盲测的三步法合成数据不仅能用来学习更是检验反演流程可靠性的标准手段。每到一个新工区我会先花一小时做一个盲测流程固定为三步。第一步根据已知的地质背景设定一个合理的二到三层电阻率模型把异常体的位置、大小、电阻率对比度都写死第二步用这个模型正演生成合成数据人为加入 3% 到 5% 的高斯噪声模拟野外数据的真实噪声水平第三步把合成数据当成实测数据跑一遍你计划使用的反演参数把反演模型和真实模型放在一起对比。如果异常体位置偏差超过一个网格单元电阻率值偏差超过 30%就说明这套参数吃不住这个工区的数据特点得回头调参数或检查流程。这个流程看起来简单但它能筛掉大量操作层面的错误。某次在某山前地段准备做电法勘察我先按 80 米埋深的低阻透镜体做了盲测发现反演结果把透镜体顶界抬高了约 10 米初始以为参数问题反复检查后才发现是数据文件中电极编号顺序写反了。如果没有盲测这个错误会被当成真实地质现象写进报告后果严重得多。盲测还有一个附加好处它能帮你确定这个工区用 ERT 到底能分辨多大尺寸的异常提前知道哪些目标属于“超出分辨率能力”的范围别等出了野外报告再后悔。6.2 一个我坚持的习惯现在每接手一个新测区、一套新参数、甚至一次软件版本升级我都会先跑一遍小规模的合成数据确认反演流程没有因环境变化而失效。这个习惯源自一次血泪经验某次在某矿区做高密度电法升级了 pyGIMLi 版本后没有重跑验证直接用旧参数处理新数据结果反演剖面上出现了一条与实际地质完全不符的横向高阻带耗费了整整两天排查最后发现是版本变化导致默认误差模型计算方式变了参数虽然看起来一样实际行为已经不同。那次之后我定了个规矩凡是长期项目环境版本固定下来凡是新项目或新数据先合成后实测两步走完再上正式反演。版本升级带来的行为变化不是 pyGIMLi 独有的问题而是所有数值计算工具都存在的隐藏风险盲测就是给这种风险上的保险。这套流程虽然多花半小时但能让你解释每一条异常时都底气足一点。希望帮到你少替我交这些学费。本文还有配套的精品资源点击获取