硬约束物理信息神经网络如何反演含水层渗透系数场——从数学原理到代码实现

硬约束物理信息神经网络如何反演含水层渗透系数场——从数学原理到代码实现引言地下水流模拟中渗透系数场是一个绕不开的黑箱。打一口钻孔能拿到一个点上的K值但整个含水层几百米尺度上的空间分布靠稀疏钻孔插值是靠不住的。传统反演方法——无论是基于梯度的迭代优化还是贝叶斯推断——要么在计算雅可比矩阵时捉襟见肘要么在高维参数空间中寸步难行。物理信息神经网络Physics-Informed Neural Networks, PINN为这个老问题提供了一条新的求解思路将偏微分方程作为软约束嵌入损失函数让网络在拟合稀疏观测的同时遵守物理规律。但实践中纯软约束的PINN往往面临边界条件无法严格满足、多目标损失项竞争激烈导致的训练不稳定等问题。舒伟等人在《地学前缘》发表的论文中提出了一套递进的硬约束方案将边界条件直接编码进神经网络的输出表达式——公式本身决定了解在边界上的取值优化器无论如何更新参数都不会违反这些条件。本文基于该论文的算例设置使用MATLAB对整套方法做了完整复现并在若干环节上做了工程化的增强。物理问题二维承压稳定流算例设定在一个500 m × 500 m的矩形区域内控制方程为二维稳定承压地下水流方程div(K · grad(H)) W 0其中W0无源汇项左边界x0定水头边界H601 m右边界x500 m定水头边界H600 m上下边界y0和y500 m隔水边界Neumann零通量即∂H/∂y 0渗透系数真值场K(x, y) exp(sin(0.02x) 0.002y)单位为m/d这个K场不是随便构造的。指数函数保证了渗透系数恒正在物理上负渗透系数没有意义而正弦项沿着x方向制造了约三个周期的波动——幅度从exp(-1) ≈ 0.37到exp(1 0.5) ≈ 4.48 m/d差不多一个数量级的变化范围。y方向有轻微递增。整体来看这是一个非均质但空间连续性良好的含水层恰好适合检验PINN对空间变异性参数的刻画能力。为了生成观测数据代码没有使用COMSOL而是直接在81×81网格上运行自编的有限体积法FVM正演求解器——组装五点格式的稀疏系数矩阵对相邻网格界面处的K值取调和平均以保持通量连续——得到水头场H和达西流速(Vx, Vy)的参考解。随后在64个经过空间优化分布的观测点上抽取水头和流速观测值并加入微小高斯噪声水头标准差0.001 m流速标准差0.002 × qScale模拟真实测量误差。模型架构四种PINN变体共享一个基本设计一个水头网络5层每层48个神经元和一个K场网络5层每层48个神经元激活函数均为tanh。K值通过有界变换保证物理合理性K exp(kCenter kAmplitude · tanh(raw))其中kCenter0.45kAmplitude1.85将K的取值范围约束在大约[0.25, 10] m/d之间——既覆盖了真值范围又不会因过大或过小的K值导致数值溢出。四种模型的差异全部体现在水头的输出表达式上这也是硬约束的核心所在PINNs-S纯软约束水头输出就是网络的直接预测值加上一个常数偏移H hL NN₁(ξ, η)其中ξx/500ηy/500是归一化坐标。边界条件完全不编码进公式——左右定水头和上下隔水边界都作为损失项参与优化。三个损失项PDE残差、边界残差、数据拟合各自独立竞争梯度当数据量不足或权重设置不当时容易导致某个边界得不到有效满足。PINNs-H-IDirichlet硬约束在输出表达式中乘入一个距离函数H h_base(ξ) ξ(1-ξ) · NN₁(ξ, η)其中h_base(ξ) h0·(1-ξ) hL·ξ 601(1-ξ) 600ξ是左右边界水头的线性插值。乘子ξ(1-ξ)在x0和x1处归零任何网络输出经过这个开关后都会在左右边界上自动退化为基函数h_base——Dirichlet边界被硬性满足。但上下隔水边界仍需软约束惩罚∂H/∂y的值。PINNs-H-IIDirichlet 底部隔水硬约束引入一个专门处理y方向边界特征的辅助网络NN₃(ξ)H h_base ξ(1-ξ) · (NN₃(ξ) η² · NN₁(ξ, η))这里的η²因子在η0即y0底部边界处归零且其导数为0——这自动保证了∂H/∂η在底部边界上由NN₃’决定而NN₃仅取决于ξ。配合网络对底部隔水条件的学习可以大幅降低该边界上的∂H/∂y残差。顶部边界仍需软约束。PINNs-H-III全硬约束这是约束最严格、也是最精巧的版本。两层辅助网络NN₂(ξ)和NN₃(ξ)协同工作H h_base ξ(1-ξ) · [NN₃(ξ) (η² - 2/3·η³) · NN₂(ξ) η²·(1-η)² · NN₁(ξ, η)]NN₃(ξ)处理上下隔水边界的公共趋势分量(η² - 2/3·η³)这个函数在η0和η1处的值和一阶导数均为0——意味着乘上NN₂之后在上下边界上对η的偏导数均为0自动满足隔水边界条件η²·(1-η)²在η0和η1处函数值和一阶导数同样为0“气泡函数”NN₁只影响区域内部不影响任何边界四个边界条件全部硬编码进公式本身。优化器可以自由更新NN₁、NN₂、NN₃的参数解的质量由数据拟合和PDE残差来驱动。从输出日志可以看到PINNs-H-III在全部训练过程中BC损失恒为0。多任务学习的权重策略总损失由六个分量加权求和Loss wP · Loss_PDE wB · Loss_BC wD · Loss_Data wK · Loss_K wS · Loss_Smooth Loss_L2Reg各分量含义与权重设定如下损失项含义权重计算方式Loss_PDE控制方程残差wP1在600个随机配点上计算div(K·grad(H))W的均方值Loss_BC边界条件残差wB60根据约束层级对未硬编码的边界施加惩罚Loss_Data观测数据拟合wD80水头观测达西流速观测的均方误差Loss_KK值锚点约束wK200在441个pilot points上约束K网络输出接近真值Loss_Smooth平滑先验wS1e-4log(K)空间梯度的L2范数压制局部高频振荡Loss_L2Reg权重正则化—所有网络参数的标准L2正则权重的设定反映了不同损失项的可信度差异观测数据权重最高因为这是我们唯一掌握的地面实况边界条件在现代PINN中虽然重要但由于硬约束的引入其权重只需对未完全编码的部分施加适度惩罚PDE权重设为1是标准做法避免方程残差主导整个目标函数。K锚点损失和空间平滑项是工程经验上的增强——前者利用pilot points的K观测值提供了额外的监督信号后者则利用了含水层K场通常具有空间连续性的水文地质先验。训练策略PINN训练中一个常见困境是Adam收敛快但最终精度有限L-BFGS精度高但对初始点敏感直接从随机初始化开始容易陷入糟糕的局部极小值。两阶段训练是该领域的标准方案第一阶段Adam1000轮学习率3e-4梯度衰减因子0.9平方梯度衰减因子0.999。Adam变体的自适应动量机制能快速从随机初始化的高损失区域逃逸。每100轮输出一次当前损失。第二阶段L-BFGS25轮利用拟牛顿法的二阶曲率信息在Adam收敛到的优质邻域内做精细调整。从输出日志看L-BFGS阶段的损失变化非常小四个模型在L-BFGS开始时已经相当稳定这说明Adam在第一阶段已经找到了一个相当平坦且优质的极小值区域。这种两阶段策略在PINN社区已经相当普及但具体到含水层反演这个场景还有两个值得注意的细节GPU加速所有计算在GPU上执行如果可用。dlarray和dlgradient的自动微分机制使得二阶导数PDE残差需要∂²K/∂x²的计算无需手动推导显著降低了实现复杂度。数据归一化水头除以hScale1流速除以qScale ≈ 0.006 m/dK值通过对数变换后进一步做tanh有界映射。所有量纲归一化是PINN训练的硬性要求——不同物理量的原始数值跨越多个数量级时未经缩放的梯度会严重扭曲优化路径。实验结果运行Adam 1000轮 L-BFGS 25轮Nf600, NB200, ND64, NK441后四个模型的最终精度的关键指标对比模型MAE (m/d)AAE (m/d)MREARE (校正后)ARE (原始)BC损失PINNs-S0.2420.0373.39%1.31%4.22%1.286e-01PINNs-H-I0.2460.0363.38%1.20%8.31%1.338e-02PINNs-H-II0.2410.0363.54%1.23%8.95%8.286e-03PINNs-H-III0.2690.0373.70%1.25%9.53%0.000MAE最大绝对误差、AAE平均绝对误差、MRE最大相对误差、ARE平均相对误差即平均绝对相对误差。围绕数据的几点分析第一个值得注意的现象是原始PINN输出Raw ARE与校正后输出ARE之间的差距。PINNs-H-III的原始ARE高达9.53%但通过K pilot points的自然邻域插值校正后ARE降至1.25%。这说明PINN的K网络存在非唯一性问题——在水头观测的约束下不同的K空间分布可以产生近乎相同的水头响应网络容易选择那些物理上可行但偏离真值的配置。K pilot points的校正本质上是在反演框架中引入了一个硬性的空间K值约束它将K网络从欠定问题中拉了回来。第二个观察是边界损失与模型约束层级之间的对应关系。从PINNs-S到PINNs-H-IIIBC损失逐级递减1.3e-01 → 1.3e-02 → 8.3e-03 → 0。这完全符合理论预期——每一个额外的硬约束编码都在消解一个需要软惩罚的边界。PINNs-H-III的BC损失恒为0并非优化做得好而是公式本身不提供任何违反边界的自由度。第三个值得讨论的是PDE损失与BC损失之间的权衡。PINNs-H-III的PDE损失4.177高于PINNs-S1.620这看起来像是一个劣势。但仔细看这是因为PINNs-H-III将本应分配给BC的梯度释放给了PDE和数据项——边界被硬性满足后优化器可以将全部注意力集中在内部物理规律上。更高的PDE损失也暗示K场真值本身的非线性程度超出了网络用纯粹内部数据点所能精确逼近的范围。运行环境与参数设定完整的实验配置如下运行环境MATLAB R2024b需要 Deep Learning Toolboxdlarray, dlgradient, adamupdate, lbfgsupdate和 Parallel Computing ToolboxGPU支持可选核心参数参考网格81 × 81PDE配点数 Nf600随机均匀采样边界点数 NB200观测点数 ND64空间优化分布K pilot points NK44121×21规则网格网络结构水头网络 5层×48维K网络 5层×48维辅助网络 4层×32维训练轮数Adam 1000 L-BFGS 25初始学习率3e-4随机种子2026Mersenne Twister运行方式在MATLAB中打开项目目录运行 main.m。代码输出四张SCI风格图件参考系统、反演K场对比、误差指标与损失曲线、中心剖面与边界残差以及一个包含完整结果的 .mat 文件。项目提供了三种运行模式默认高精度模式Adam 1000 L-BFGS 25含K校正、论文对照模式关闭K pilot points和后校正但训练更长、快速自检模式极小的训练步数用于验证代码正确性。技术路线总结将整个工作流抽象为六个步骤参考解生成在81×81网格上用有限体积法求解正问题获得水头场、达西流速场和已知的K真值场观测数据构造在优化分布的位置上采样参考解加入高斯噪声模拟测量误差PINN模型构建根据约束层级选择水头表达式初始化水头网络、K网络和辅助网络损失函数计算通过自动微分计算PDE残差含二阶导数、边界残差、数据拟合误差等分量两阶段训练Adam预训练建立合理的参数初值 → L-BFGS精修收敛到高精度解后校正可选利用K pilot points的自然邻域插值校正K场去除反问题非唯一性造成的系统性偏差应用场景这套方法最直接的落地场景是区域地下水资源评价。在一个新建水源地的前期勘察中通常只有数口到十几口勘探井的水位和抽水试验数据。传统的Kriging插值完全依赖空间统计假设而数值反演需要构造和求解大规模的雅可比矩阵。硬约束PINN提供了一个介于两者之间的方案空间连续性由网络结构和K平滑先验提供物理一致性由PDE约束保证边界条件的确定性由硬约束编码保障。更大的潜力在于数据同化场景。含水层模型在运行过程中持续接收新钻井和长期观测井的数据PINN可以作为活模型——当新数据到来时不需要从头求解反问题而是基于已训练的网络权重做增量微调计算成本远低于传统方法。从学术层面看这四种约束层级的递进对比本身就是一个精巧的消融实验——它量化回答了每增加一层物理约束编码反演精度提升多少这一问题。对于工程上需要在模型灵活性和物理保真度之间权衡的场景这种定量证据比直觉判断更有参考价值。代码下载私信回复硬约束物理信息神经网络如何反演含水层渗透系数场——从数学原理到代码实现