简介:这是一份基于开源软件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(mesh=mesh, res=res, scheme="wenner", electrodePositions=sensors) print(data)这里解释几个关键参数。pg.createGrid的两个数组分别是 x 方向和 y 方向的网格边界,y 用负值代表深度方向,这是 pyGIMLi 的习惯,很多新人在这里把正负搞反,结果模型跑到地表以上去了。pg.RVector是 pyGIMLi 内部的向量类型,传入单元数和一个初始值,生成与网格单元一一对应的电阻率数组。simulate方法的scheme="wenner"表示用温纳装置生成测量序列,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(maxIter=8) print("反演完成,模型单元数:", len(model))这一段代码能跑通,是因为ERTManager内部自动完成了网格生成、初始模型估计、正演算子和反演迭代的控制。model返回的是一个向量,长度与反演网格的单元数一致,每个元素是该单元的电阻率。跑完后调用mgr.showResult()会在当前绘图窗口画出一条带深度刻度的电阻率剖面。
需要注意,默认网格和默认参数只是给了一个“能出结果”的下限。我第一次用默认参数跑一块已知模型时,得到的剖面形状大致对,但异常体边界模糊、幅度被压低,这说明默认的平滑约束起了作用,但也把真实对比度磨掉了。所以最小反演命令适合用来确认整条流程通不通,真正要用于解释的结果,必须进入参数调整这一步。在反演圈里有句俗话:反演出来的不是真实模型,而是满足约束条件的一个平滑版本。这句话不是泼冷水,而是提醒你,参数选型决定了结果往哪个方向偏。
3.2 反演参数选型:lam、zWeight、maxIter、blocky 到底调什么
ERT 反演本质上是在拟合实测数据和保持模型平滑之间找一个平衡点,ERTManager.invert()提供的那几个参数就是这个平衡点的控制器。下面是我常用的参数表,按重要程度排:
| 参数 | 默认值 | 作用 | 调整方向 |
|---|---|---|---|
| lam | 自动估计 | 正则化权重,越大模型越平滑 | 过拟合时增大,欠拟合时减小 |
| zWeight | 1.0 | 垂直方向与水平方向平滑的权重比 | 地层水平层状时设 0.5~0.8 |
| maxIter | 默认 10 次左右 | 迭代次数上限 | 收敛慢时调大,但别超过 30 |
| blocky | False | 使用 L1 范数模型约束 | 异常体边界锐利时设 True |
| robust | False | 对数据中的离群点降权 | 数据干扰大时设 True |
实际调用时这样传参:
import pygimli as pg from pygimli.physics import ert mgr = ert.ERTManager(data) model = mgr.invert(lam=20, zWeight=0.8, maxIter=15, blocky=False)参数不是越多越好,我建议先动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, relativeError=0.05, absoluteError=0.001)relativeError=0.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(maxIter=10) mgr.saveResult(outdir="result_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% 附近降不下去,模型剖面却平滑得很漂亮。这种情况单独看剖面容易误判为成功,实际上正演响应和实测数据严重不吻合,剖面是平滑约束硬推出来的,不是数据驱动的。
原因多半是数据里有系统误差,或者误差模型设得不够,反演算法只能通过拉平模型来妥协。解决分两步:先检查实测数据与正演响应逐点比较,找出偏离最大的测量点,确认是不是野值;再确认你的相对误差设置是否符合实际仪器精度。如果数据本身没问题,就把robust=True打开,让它对离群点自动降权,往往能压住 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 独有的问题,而是所有数值计算工具都存在的隐藏风险,盲测就是给这种风险上的保险。这套流程虽然多花半小时,但能让你解释每一条异常时都底气足一点。希望帮到你,少替我交这些学费。
本文还有配套的精品资源,点击获取