简介SimPEG.FLOW是基于SimPEG框架的开源Python子模块专门用于地下水流方程的数值模拟与参数反演面向地球物理学家、水文学者以及从事环境地质、水资源管理的工程师。压缩包内提供了一套完整的工程结构共20个文件其中9个.py脚本为核心代码覆盖达西定律实现、Dirichlet/Neumann/Robin三类边界条件设置、网格离散与反演优化调用2个.rst文档与Makefile辅助生成说明LICENSE、.gitignore、.coveragerc等保证项目规范examples目录提供可直接运行的示例docs目录含API文档。整个包仅45KB体量轻便。已有492人学习下载。通过研读此工程可快速掌握SimPEG.FLOW从正演到反演的整体逻辑理解自适应网格构建与多物理场耦合方式并能与NumPy、SciPy、Matplotlib等科学计算生态无缝衔接将示例迁移至污染物扩散、水源保护、储层评价等实际场景是入门SimPEG地下水模拟的一条高效路径。1. 用SimPEG写地下水流方程一次把水文和地球物理参数装进同一张网格的尝试拿电磁法反演框架来写地下水流方程听起来像是用地球物理的锤子砸水文地质的钉子。但真在项目里写过一轮就会发现这一锤正好砸在联合反演最疼的那条缝上。simpegflow 的思路很直接与其让地下水模型和电阻率模型各自用一套网格、一套离散格式最后在两个黑匣子之间来回倒换参数不如从正演第一步就共用 SimPEG 的网格与算子。于是达西定律被写成了和直流电法几乎一样的稀疏矩阵乘法K 模型和电导率模型落在同一个单元中心上后面想做水文-地球物理联合反演算梯度时再也不用头疼坐标映射。这份代码量不大但把“有限体积离散、边界条件、稀疏矩阵装配、参数化”这条正演链路完整走了一遍。适合两类人看一类是刚接触 SimPEG、想搞懂它内部算子怎么拼出物理方程的人另一类是手里有水文观测数据、想往反演方向扩展的开发者。前者能照着跑通后者能看到怎么把 K 模型换成反演变量。2. 方程离散化为什么达西定律和直流电法在 SimPEG 里是同一套矩阵2.1 从达西方程到泊松型方程数学同构稳态承压含水层里的地下水流忽略源汇项时满足达西定律和质量守恒的组合∇·(K ∇h) 0其中 h 是水头mK 是水力传导率m/s。把它和直流电阻率法的控制方程放在一起看∇·(σ ∇φ) Iσ 是电导率φ 是电位。两个方程除了符号不同结构完全一样都是椭圆型方程。这意味着只要有一套可靠的离散化工具能把“梯度算子、物性系数、散度算子”拼成系统矩阵那么水电和地电的正演就是同一套装配流程。SimPEG 正好在电磁正演里把这套流程打磨得很成熟拿来写地下水流方程本质上只是换了一个物性参数的名字。这个同构关系不是纸面上的巧合。在污染羽监测这类场景里含水层的 K 和电导率往往同时受岩性和污染物浓度影响两个参数在空间上有很强相关性。用同一套网格做两种正演参数就能直接逐单元绑定联合反演的目标函数里加一个耦合项就能把两类观测约束到同一组模型上。这个价值在传统水文软件和地球物理软件分开跑的工作流里很难实现。2.2 有限体积离散化的三个算子梯度、界面平均、散度SimPEG 对这类方程的离散走的是有限体积路线。水头 h 定义在单元中心cell center而通量 K∇h 定义在网格面上face。从单元中心到最终的单元平衡方程需要三个算子衔接cell_gradient把单元中心的水头差分到面中心得到每个面上的 ∇h 分量界面平均把单元中心的 K 插值到面中心得到面位置的等效渗透率 K_faceface_divergence把每个面上已经乘过 K 的通量做散度回到单元中心形成质量守恒方程以二维 TensorMesh 为例假设有 nC 个单元、nF 个面三个算子的维度关系如下算子输入位置输出位置矩阵形状cell_gradient单元中心 (nC)面中心 (nF)nF × nC界面平均单元中心 (nC)面中心 (nF)nF × nCface_divergence面中心 (nF)单元中心 (nC)nC × nF值得注意的是SimPEG 里的 cell_gradient 和 face_divergence 在离散符号上恰好满足“散度是梯度的负伴随”这类相容性关系这是有限体积方法保证无虚假振荡的关键。自己手工填差分矩阵很容易破坏这种配对关系而直接用框架的算子至少矩阵层面的兼容性有保障。2.3 系统矩阵到底长什么样A D M K G把三个算子串起来系统矩阵的装配就是一行稀疏矩阵乘法A Div · diag(K_face) · Grad这个 A 的物理含义是每个单元方程等于它所有面上的通量之和为零。对角线上的值代表流出该单元的总能力非对角线上的值代表相邻两个单元之间的水力联系。和有限元法相比这里没有“质量矩阵”加权也不需要对每个单元单独推导单元刚度矩阵本质上是把有限体积的逐个单元组装压缩成了三个全局稀疏矩阵的乘积。由于水流问题是线性的边界水头确定之后求解 A h b 得到的就是整个网格的水头分布。这和传统地下水软件内部的做法在数学上是等价的差别在于K 被明确地保留成一个可以在单元级别随意赋值的向量后面想把它替换成 exp(m)、想加各向异性张量、想做灵敏度计算都是在矩阵装配层动手而不用爬进某个模拟器的源码里改。3. 核心模块拆解一次从 K 到水头场的完整装配3.1 工程结构单文件正演也能拆成三层拿到 simpegflow 这类代码我习惯先把它整理成最小可用的三层结构方便后面加功能时不动主流程。实际跑通只需要正演入口、模型参数和可视化三个文件我一般会拆成这样simpegflow/ ├─ requirements.txt # SimPEG、numpy、scipy、matplotlib ├─ models.py # K 模型构建均质、透镜体、层状含水层 ├─ flow_forward.py # 网格生成、算子装配、边界条件、求解 └─ viz.py # 水头等值线绘制、VTK 导出flow_forward.py 是这个资源的核心其余文件都围绕它工作。下面按正演链路从网格到出图逐个过一遍代码可以直接保存成文件跑。3.2 TensorMesh 网格让 K 和电阻率落在同一个单元上第一步是生成网格。simmpegflow 选择 TensorMesh 而不是非结构网格主要是为了和电磁正演的网格习惯对齐。TensorMesh 在规则区域上可以用一个向量描述一个方向上的剖分写起来非常紧凑而且 face_divergence、cell_gradient 这些算子都是预先编译好的跑起来比逐个单元循环快几个数量级。import numpy as np import scipy.sparse as sp from scipy.sparse.linalg import spsolve import matplotlib.pyplot as plt from SimPEG import TensorMesh # 100 x 100 个单元每个单元 2 m模拟 200 m x 200 m 区域 hx np.ones(100) * 2.0 hy np.ones(100) * 2.0 mesh TensorMesh([hx, hy], x0[0.0, 0.0]) nC mesh.nC # 水力传导率背景值 1e-4 m/s中部放一个半径 30 m 的透镜体K 1e-2 m/s K_model np.ones(nC) * 1e-4 cx mesh.cell_centers[:, 0] cy mesh.cell_centers[:, 1] inside ((cx - 100.0) ** 2 (cy - 100.0) ** 2) 30.0 ** 2 K_model[inside] 1e-2给 hx 一个等距向量就够了如果想在井点附近加密网格把 hx 换成不等距数组比如hx np.r_[np.ones(20)*2, np.ones(30)*1, ...]完全不需要改后续代码。这也是 TensorMesh 在工程里最好用的点加密区域和粗化区域用一个向量表达清楚算子维度自动跟着变。3.3 组装系统矩阵G、MfK、D 三段乘法系统矩阵的组装放在一个独立函数里这样后面做反演、做瞬态扩展时只需要调用它返回 A不用关心内部的 FV 细节。界面上的 K 一定用谐波平均不能直接在两个相邻单元里取算术平均。def assemble_A(mesh, K_model): # 1. 梯度单元中心水头 - 面中心水头梯度 Grad mesh.cell_gradient # 2. 散度面中心通量 - 单元中心通量散度 Div mesh.face_divergence # 3. 界面上 K 取谐波平均先取倒数平均再倒数 K_face 1.0 / mesh.average_cell_to_face(1.0 / K_model) # 4. 系统矩阵A Div * diag(K_face) * Grad A Div sp.diags(K_face) Grad return A.tocsr(), Grad, Div, K_face谐波平均的物理原因很直接水流串联穿过两个相邻单元等效渗透率应该由低渗透率那侧主导算术平均会高估流量。SimPEG 的average_cell_to_face默认是算术平均所以先在单元中心取 1/K 的倒数做平均保证界面位置得到的是正确的串联等效值。这在含水层里有透镜体、有细砂夹层时影响特别明显算出的水头等值线会整体偏移。3.4 边界条件与求解固定水头怎么塞进稀疏矩阵这个项目里左边界固定水头 10 m右边界固定水头 0 m上下边界默认零通量。上下边界不用特殊处理因为有限体积离散里没有在边界面上强行施加通量的地方自然就是无通量左右边界则需要把对应单元的行替换成单位行固定水头值。def solve_flow(mesh, K_model, h_left10.0, h_right0.0): A, Grad, Div, K_face assemble_A(mesh, K_model) rhs np.zeros(mesh.nC) # 找到左右两列单元 cc mesh.cell_centers left cc[:, 0] cc[:, 0].min() right cc[:, 0] cc[:, 0].max() bc left | right # 固定水头用单位行替换原方程行 A A.tolil() A[bc, :] 0.0 A[bc, bc] 1.0 rhs[bc] h_left # 左侧赋值 10 rhs[bc (right)] h_right # 右侧赋值 0 A A.tocsc() h spsolve(A, rhs) return h, A, Div, K_face这里容易踩的一个小坑是修改稀疏矩阵的行之前先转成 LIL 格式否则直接在 CSR 上按行赋值会非常慢而且容易破坏内部索引。另一个注意点是 rhs 初始化为全零只对边界单元写水头值。如果漏了这一步固定水头条件实际上没有生效解出来的是全零场但程序不会报任何错。3.5 结果输出水头等值线和 VTK 持久化求解完成后把 h 按网格形状写回二维数组直接画等值线。如果需要进 ParaView 查看三维场或做后续后处理一次性导出 VTK 文件是最省事的。h, A, Div, K_face solve_flow(mesh, K_model) # 画等值线 h2d h.reshape((mesh.nCy, mesh.nCx), orderC) plt.contourf(cx.reshape(mesh.nCy, mesh.nCx), cy.reshape(mesh.nCy, mesh.nCx), h2d, levels20) plt.colorbar(labelhydraulic head (m)) plt.savefig(head_field.png, dpi150) # 导出 VTK方便在 ParaView 里加载 mesh.write_UBC? # 如果不熟悉直接用下面的方式导出 VTK 我一般用mesh.write_vtk(head.vtk, {head: h, logK: np.log10(K_model)})。这个 API 在老版本和新版本里基本稳定参数就是文件名加字典字典里每个键对应一个单元中心的标量场。导出来后水头场和 logK 场叠加看能直观确认透镜体对水头分布的影响是否符合预期。4. 常见问题排查五个能把水头算飞的坑4.1 高渗透镜体附近流量虚高界面 K 平均方式不对现象同样一组 K 模型透镜体轮廓画出来总是“胖”一圈流量比手算解析解大一倍。原因在把单元中心 K 插值到面中心时用了算术平均。算术平均对两个相邻单元的渗透率一视同仁但实际水流是串联穿过界面的等效值应该由低渗透率单元主导。高渗透镜体附近一旦用了算术平均界面等效 K 被明显抬高水流全部从透镜体附近“抄近路”水头等值线会整体扭曲。解决统一改成谐波平均先对 1/K 做算术平均再取倒数也就是代码里的1.0 / mesh.average_cell_to_face(1.0 / K_model)。改完之后和解析解对比误差从两位数降到个位数。4.2 固定水头边界没生效改了 A 却忘了 rhs现象左右边界明明设了 10 m 和 0 m画出来的水头场却是全局接近 0 的平坦场。原因只把 A 的边界行替换成单位行但 rhs 对应位置没有写入水头值。稀疏矩阵求解时右侧向量全是零边界行求解得到的就是零。解决改 A 的同时同步rhs[bc] 10.0右侧另写一行rhs[right] 0.0。养成把 A 和 rhs 合起来当“一个系统”来改的习惯两条赋值语句挨着放中间不插别的逻辑。4.3 解出来水头方向反了梯度算子的符号约定现象左边界 10 m、右边界 0 m结果左边水头反而更低整个场整体反号。原因SimPEG 的 cell_gradient 在离散时会按某个方向约定做差分拼出来的系统矩阵不一定是我们手写推导里想要的 M-matrix符号方向取决于网格坐标和算子实现。这一点在不同版本里表现不一致遇到反号不用慌。解决组装后先打印一个小网格的 A 对角元和相邻非对角元符号判断方向。方向反了整体给 A 和 rhs 各乘 -1 即可不会改变解的形态只是把符号拨正。关键是别只在 A 上取负而忘了 rhs否则解会真的翻转。4.4 运行一切正常但结果像“打印错误”单位不统一现象把渗透系数 cm/s 直接代进去水头场没问题但算流速或者和实测抽水数据对比时数值差了几个数量级。原因网格长度用 m水头用 m而 K 用的是 cm/s通量单位 mix 了。SimPEG 不检查单位所有数值都是裸数字单位错误只在后处理阶段才暴露。解决进入代码前先把所有量纲统一写在一行注释里长度 m、时间 s、水头 m、K m/s。K 从 cm/s 换算到 m/s 时统一除以 100不要在代码里散着写换算系数否则换个人维护一定会漏。4.5 直接把 K 当反演参数梯度更新一步就崩现象想拿水头观测数据反演 K用 Gauss-Newton 迭代第一次更新后模型出现负值第二次直接 NaN。原因K 跨了几个数量级最小二乘目标函数对 logK 的敏感度比对 K 平稳得多直接在 K 空间做线性更新步长稍大就冲进负值区。解决反演前把模型参数强制改成 θ log(K)正演函数内部用 exp(θ) 还原 K。这一步在 simpegflow 的正演入口里预留很自然输入 logK输出水头场其余代码一行不用动。后面接 SimPEG 的优化模块时梯度也天然在 log 空间里计算非负约束自动满足。5. 二次改造把它做成水文-地球物理联合反演的共用正演器5.1 封装成 FlowForward 类直接改脚本里的全局变量容易把稳定版本改坏我拿到手后第一件事是把正演链路封装成一个类让网格、边界水头、观测井位置都变成构造参数。class FlowForward: def __init__(self, mesh, h_left10.0, h_right0.0): self.mesh mesh self.h_left h_left self.h_right h_right self.Grad mesh.cell_gradient self.Div mesh.face_divergence def solve(self, logK): K_model np.exp(logK) A, rhs self.assemble(K_model) return spsolve(A, rhs) def assemble(self, K_model): K_face 1.0 / self.mesh.average_cell_to_face(1.0 / K_model) A self.Div sp.diags(K_face) self.Grad # ... 边界行替换同上略 return A.tocsr(), rhs封装完的好处是后面的反演脚本可以完全不知道有限体积细节只需要把 logK 传进去再拿残差。5.2 接上反演同一套网格上跑两种正演联合反演的接法也很自然。把直流电法正演和地下水流正演各自封装成类共享同一个 TensorMesh 对象模型的物理参数映射只发生在一个函数里。def residual(logK): K_model np.exp(logK) h_pred flow_fwd.solve(logK) # 水文观测 phi_pred dc_fwd.solve(K_model_to_sigma(K_model)) # 电法观测 return np.r_[h_obs - h_pred[obs_cells], phi_obs - phi_pred[dc_obs_cells]]这个残差函数可以直接丢进 scipy 的 least_squares也可以接 SimPEG 的 optimization 模块做带正则化的反演。因为两种观测的网格完全一致参数梯度没有任何空间插值误差。后来我在做某场地污染羽调查时就是把水文观测井的水头和地面电法数据拼在一个目标函数里两次正演用的是同一个缝出来的 K 场收敛速度比分开反演再叠图快得多。从那以后我每次把新物理方程塞进 SimPEG都强制自己先做一步一维解析解对比再上二维网格跑完整的装配和边界条件确认无误后再往反演方向扩。这套最小验证流程救过我不少次希望帮到你。本文还有配套的精品资源点击获取