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

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

引言

地下水流模拟中，渗透系数场是一个绕不开的"黑箱"。打一口钻孔能拿到一个点上的K值，但整个含水层几百米尺度上的空间分布，靠稀疏钻孔插值是靠不住的。传统反演方法——无论是基于梯度的迭代优化还是贝叶斯推断——要么在计算雅可比矩阵时捉襟见肘，要么在高维参数空间中寸步难行。

物理信息神经网络（Physics-Informed Neural Networks, PINN）为这个老问题提供了一条新的求解思路：将偏微分方程作为软约束嵌入损失函数，让网络在拟合稀疏观测的同时"遵守"物理规律。但实践中，纯软约束的PINN往往面临边界条件无法严格满足、多目标损失项竞争激烈导致的训练不稳定等问题。

舒伟等人在《地学前缘》发表的论文中提出了一套递进的硬约束方案，将边界条件直接编码进神经网络的输出表达式——公式本身决定了解在边界上的取值，优化器无论如何更新参数都不会违反这些条件。本文基于该论文的算例设置，使用MATLAB对整套方法做了完整复现，并在若干环节上做了工程化的增强。

物理问题：二维承压稳定流

算例设定在一个500 m × 500 m的矩形区域内：

• 控制方程为二维稳定承压地下水流方程：div(K · grad(H)) + W = 0，其中W=0（无源汇项）
• 左边界x=0：定水头边界，H=601 m
• 右边界x=500 m：定水头边界，H=600 m
• 上下边界y=0和y=500 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))

其中kCenter=0.45，kAmplitude=1.85，将K的取值范围约束在大约[0.25, 10] m/d之间——既覆盖了真值范围，又不会因过大或过小的K值导致数值溢出。

四种模型的差异全部体现在水头的输出表达式上，这也是"硬约束"的核心所在：

PINNs-S：纯软约束

水头输出就是网络的直接预测值加上一个常数偏移：

H = hL + NN₁(ξ, η)，其中ξ=x/500，η=y/500是归一化坐标。

边界条件完全不编码进公式——左右定水头和上下隔水边界都作为损失项参与优化。三个损失项（PDE残差、边界残差、数据拟合）各自独立竞争梯度，当数据量不足或权重设置不当时，容易导致某个边界得不到有效满足。

PINNs-H-I：Dirichlet硬约束

在输出表达式中乘入一个"距离函数"：

H = h_base(ξ) + ξ(1-ξ) · NN₁(ξ, η)

其中h_base(ξ) = h0·(1-ξ) + hL·ξ = 601(1-ξ) + 600ξ，是左右边界水头的线性插值。乘子ξ(1-ξ)在x=0和x=1处归零，任何网络输出经过这个"开关"后都会在左右边界上自动退化为基函数h_base——Dirichlet边界被硬性满足。但上下隔水边界仍需软约束惩罚∂H/∂y的值。

PINNs-H-II：Dirichlet + 底部隔水硬约束

引入一个专门处理y方向边界特征的辅助网络NN₃(ξ)：

H = h_base + ξ(1-ξ) · (NN₃(ξ) + η² · NN₁(ξ, η))

这里的η²因子在η=0（即y=0，底部边界）处归零，且其导数为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

各分量含义与权重设定如下：

权重的设定反映了不同损失项的"可信度"差异：观测数据权重最高，因为这是我们唯一掌握的地面实况；边界条件在现代PINN中虽然重要，但由于硬约束的引入，其权重只需对未完全编码的部分施加适度惩罚；PDE权重设为1是标准做法，避免方程残差主导整个目标函数。

K锚点损失和空间平滑项是工程经验上的增强——前者利用pilot points的K观测值提供了额外的监督信号，后者则利用了含水层K场通常具有空间连续性的水文地质先验。

训练策略

PINN训练中一个常见困境是：Adam收敛快但最终精度有限；L-BFGS精度高但对初始点敏感，直接从随机初始化开始容易陷入糟糕的局部极小值。两阶段训练是该领域的标准方案：

第一阶段（Adam，1000轮）：学习率3e-4，梯度衰减因子0.9，平方梯度衰减因子0.999。Adam变体的自适应动量机制能快速从随机初始化的高损失区域逃逸。每100轮输出一次当前损失。

第二阶段（L-BFGS，25轮）：利用拟牛顿法的二阶曲率信息在Adam收敛到的优质邻域内做精细调整。从输出日志看，L-BFGS阶段的损失变化非常小（四个模型在L-BFGS开始时已经相当稳定），这说明Adam在第一阶段已经找到了一个相当平坦且优质的极小值区域。

这种两阶段策略在PINN社区已经相当普及，但具体到含水层反演这个场景，还有两个值得注意的细节：

1. GPU加速：所有计算在GPU上执行（如果可用）。dlarray和dlgradient的自动微分机制使得二阶导数（PDE残差需要∂²K/∂x²）的计算无需手动推导，显著降低了实现复杂度。

2. 数据归一化：水头除以hScale=1，流速除以qScale ≈ 0.006 m/d，K值通过对数变换后进一步做tanh有界映射。所有量纲归一化是PINN训练的硬性要求——不同物理量的原始数值跨越多个数量级时，未经缩放的梯度会严重扭曲优化路径。

实验结果

运行Adam 1000轮 + L-BFGS 25轮（Nf=600, NB=200, ND=64, NK=441）后，四个模型的最终精度的关键指标对比：

MAE（最大绝对误差）、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-III，BC损失逐级递减：1.3e-01 → 1.3e-02 → 8.3e-03 → 0。这完全符合理论预期——每一个额外的硬约束编码都在消解一个需要软惩罚的边界。PINNs-H-III的BC损失恒为0并非优化"做得好"，而是公式本身不提供任何违反边界的自由度。

第三个值得讨论的是PDE损失与BC损失之间的权衡。PINNs-H-III的PDE损失（4.177）高于PINNs-S（1.620），这看起来像是一个劣势。但仔细看，这是因为PINNs-H-III将本应分配给BC的梯度"释放"给了PDE和数据项——边界被硬性满足后，优化器可以将全部注意力集中在内部物理规律上。更高的PDE损失也暗示K场真值本身的非线性程度超出了网络用纯粹内部数据点所能精确逼近的范围。

运行环境与参数设定

完整的实验配置如下：

运行环境：MATLAB R2024b+，需要 Deep Learning Toolbox（dlarray, dlgradient, adamupdate, lbfgsupdate）和 Parallel Computing Toolbox（GPU支持，可选）

核心参数：

• 参考网格：81 × 81
• PDE配点数 Nf=600（随机均匀采样）
• 边界点数 NB=200
• 观测点数 ND=64（空间优化分布）
• K pilot points NK=441（21×21规则网格）
• 网络结构：水头网络 5层×48维，K网络 5层×48维，辅助网络 4层×32维
• 训练轮数：Adam 1000 + L-BFGS 25
• 初始学习率：3e-4
• 随机种子：2026（Mersenne Twister）

运行方式：在MATLAB中打开项目目录，运行 main.m。代码输出四张SCI风格图件（参考系统、反演K场对比、误差指标与损失曲线、中心剖面与边界残差），以及一个包含完整结果的 .mat 文件。

项目提供了三种运行模式：默认高精度模式（Adam 1000 + L-BFGS 25，含K校正）、论文对照模式（关闭K pilot points和后校正，但训练更长）、快速自检模式（极小的训练步数，用于验证代码正确性）。

技术路线总结

将整个工作流抽象为六个步骤：

1. 参考解生成：在81×81网格上用有限体积法求解正问题，获得水头场、达西流速场和已知的K真值场
2. 观测数据构造：在优化分布的位置上采样参考解，加入高斯噪声模拟测量误差
3. PINN模型构建：根据约束层级选择水头表达式，初始化水头网络、K网络和辅助网络
4. 损失函数计算：通过自动微分计算PDE残差（含二阶导数）、边界残差、数据拟合误差等分量
5. 两阶段训练：Adam预训练建立合理的参数初值 → L-BFGS精修收敛到高精度解
6. 后校正（可选）：利用K pilot points的自然邻域插值校正K场，去除反问题非唯一性造成的系统性偏差

应用场景

这套方法最直接的落地场景是区域地下水资源评价。在一个新建水源地的前期勘察中，通常只有数口到十几口勘探井的水位和抽水试验数据。传统的Kriging插值完全依赖空间统计假设，而数值反演需要构造和求解大规模的雅可比矩阵。硬约束PINN提供了一个介于两者之间的方案：空间连续性由网络结构和K平滑先验提供，物理一致性由PDE约束保证，边界条件的确定性由硬约束编码保障。

更大的潜力在于数据同化场景。含水层模型在运行过程中持续接收新钻井和长期观测井的数据，PINN可以作为"活模型"——当新数据到来时，不需要从头求解反问题，而是基于已训练的网络权重做增量微调，计算成本远低于传统方法。

从学术层面看，这四种约束层级的递进对比本身就是一个精巧的消融实验——它量化回答了"每增加一层物理约束编码，反演精度提升多少"这一问题。对于工程上需要在"模型灵活性"和"物理保真度"之间权衡的场景，这种定量证据比直觉判断更有参考价值。

代码下载

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