做渗流数值模拟的人接触神经网络之后最容易问的一个问题就是能不能让网络不仅拟合数据还老老实实遵守流体力学方程答案就是 PINN也就是 Physics-Informed Neural Networks物理约束的神经网络。我最近用 PINN 对基本渗流方程做了一轮求解从损失函数搭建、配置点采样到训练调参完整走了一遍中间踩了不少坑也有了一些稳定复现的经验。这篇文章就把整个过程和关键细节掰开讲讲给同样想用 PINN 解决偏微分方程问题的朋友一个可以直接参考的路径。不管你是做地下水渗流、石油工程、岩土工程还是单纯想入门 PINN 这个方向这篇内容都适用。我会先讲清楚 PINN 到底在解什么问题再回到渗流方程本身然后给出一套可运行的一维求解思路最后把训练中常见的收敛问题和边界处理问题逐个说清楚。1. 从偏微分方程到神经网络PINN到底在解什么问题1.1 为什么传统数值方法在渗流问题中会遇到瓶颈做渗流计算的人对有限差分和有限元方法都很熟。经典的承压含水层水头分布、基坑降水、土坝渗流基本都是用这些方法在网格上离散求解。传统方法很成熟精度也可以做到很高但有几个场景会让人头疼第一参数反演。实际工程中渗透系数、储水系数往往不是直接测出来的而是通过观测井的水头数据反推。用传统方法做反演每迭代一次就要重新剖分网格、重新求解正问题计算量非常大而且目标函数容易出现局部极值。第二复杂边界和多尺度问题。裂隙介质、非均质土层、随机渗透系数场这些场景下网格剖分和质量守恒处理非常麻烦。网格太粗解不准网格太细计算成本又上去了。第三快速预测需求。比如做实时决策支持系统需要在几秒内预测不同开采方案下的水头变化传统数值模型通常做不到这种响应速度。PINN 的价值恰好在这几个点上体现出来。它不依赖网格用神经网络直接表示解场通过自动微分把偏微分方程的残差变成损失函数的一部分。训练好之后前向预测只是一次网络推理速度快到可以忽略不计而如果要做参数反演渗透系数也可以直接作为网络输入或者可训练变量一起优化避免了传统方法反复求解正问题的开销。1.2 PINN的核心思路把物理方程变成损失函数先说一个最简单的理解方式。平时做监督学习我们给网络一堆输入和标签让它学会映射关系。PINN 不一样它不只有数据标签还把物理定律直接写进训练目标里。假设我们要找的渗流问题解是 H(x,t)网络输出用 H_θ(x,t) 表示其中 θ 是网络权重。我们把 H_θ 代入偏微分方程比如 S ∂H/∂t K ∂²H/∂x²就会产生一个残差r_θ(x,t) S ∂H_θ/∂t - K ∂²H_θ/∂x²如果 H_θ 是精确解这个残差应该处处等于零。所以训练目标就变成了调整 θ让所有采样点上的 r_θ 尽可能接近零。与此同时初始条件和边界条件也通过同样的方式加入损失函数。这就是“物理约束”的含义。传统神经网络只关心输入输出之间的拟合关系而 PINN 用偏微分方程限制了网络输出必须满足物理规律。你可以把它想象成一个学生在做题不仅对答案还要每一步都符合物理定律不允许出现“看起来对但原理上站不住”的结果。实际训练时损失函数通常是好几项均方误差的加权和L λ_ic L_ic λ_bc L_bc λ_pde L_pde λ_data L_data其中 L_ic 是初始条件误差L_bc 是边界条件误差L_pde 是控制方程残差L_data 是实际观测数据误差。这种设计让 PINN 既能利用物理方程又能利用实测数据两者互相补充。1.3 基本渗流方程的数学形式渗流方程听起来复杂其实基本形式非常经典。地下水渗流遵循达西定律结合连续性方程可以得到承压含水层中水头 H 的控制方程S ∂H/∂t ∂/∂x (K ∂H/∂x) Q其中 S 是储水系数K 是渗透系数Q 是源汇项。如果介质均匀、K 为常数方程可以写成∂H/∂t a ∂²H/∂x²这里的 a K/S就是压力传导系数。它就是典型的抛物型偏微分方程和热传导方程在数学上是同一个形式。很多人刚接触 PINN 时觉得偏微分方程很难其实从最简单的扩散型方程入手是最合适的这也是我开始时选择基本渗流方程的原因。定解条件通常包括初始水头分布 H(x,0) 和边界条件比如定水头边界、定流量边界或自由排水边界。在 PINN 中这些条件不需要单独写成约束方程只需要在损失函数里增加对应的误差项就可以。2. 损失函数到底怎么设计物理约束不是加一条公式那么简单2.1 网络输入输出与激活函数选型PINN 的网络结构并不神秘通常就是一个全连接神经网络。输入是空间坐标 x 和时间 t输出是水头 H。我常用的是 3 到 5 层隐藏层每层 50 到 100 个神经元隐藏层激活函数用 tanh。激活函数的选择非常关键。ReLU 在图像分类里很好用但在 PINN 里要谨慎。因为损失函数需要计算 H 对 x 的二阶导数ReLU 的二阶导数几乎处处为零没法给方程残差提供有效梯度。tanh 是光滑函数导数可以一直求下去而且可以很好地处理中等非线性问题。如果遇到梯度消失可以考虑更深的残差连接结构或者用正弦激活函数做变体但基础场景 tanh 足够。输入归一化也是容易被忽略的点。x 和 t 如果原始尺度相差太大比如 x 是几十米、t 是几百秒网络在刚开始训练时很容易梯度失衡。我通常先把 x、t 都归一化到 [-1,1] 或者 [0,1] 区间目标变量 H 也同样压缩到 [0,1]训练稳定很多。归一化不只是为了收敛速度还直接影响损失函数中各项的量纲平衡后面会再提到。2.2 四项损失函数的具体计算方法在代码实现里每一项损失都是均方误差。初始条件损失 L_ic 是在 t0 时刻采一批点计算网络输出与初始水头值的误差。边界条件损失 L_bc 是在空间边界上采点同样计算网络输出与边界水头值的误差。方程残差损失 L_pde 是在整个时空域内部随机采点计算 H_θ 代入方程后得到的残差平方均值。数据损失 L_data 是可选项。如果在研究区域里有观测井的水头时间序列就把这些点也加入损失函数。这样网络解出来的场既符合物理方程又匹配真实观测。这个特性在传统数值方法里实现起来很别手因为传统方法的一般做法是先求解控制方程再把观测数据用数据同化方法融合进去而 PINN 天然把数据同化纳入了训练框架。有一个容易犯的错误把各项损失直接相加。因为初始条件、边界条件和方程残差的量纲可能不同数值范围差异也很大如果权重不调整训练会被量级最大的一项主导。我一般先用小网络试跑看每一项损失的数量级再设定 λ 权重。比如方程残差量级小就需要给一个较大的权重初始和边界误差如果比较顽固也可以动态调大权重。2.3 采样策略残差点不是随便撒PINN 训练依赖配置点也就是我们计算偏微分方程残差的采样点。最简单的做法是在时空区域内均匀随机采样。但实际效果告诉我单纯随机采样在边界和初始条件附近很容易欠采样特别是时间 t 接近 0 的区域物理信息最丰富却经常被忽略。我对比过几种做法均匀随机、拉丁超立方采样和基于残差的自适应加密。拉丁超立方比均匀随机好一些因为它能保证每个维度都有覆盖。自适应加密效果好但实现复杂度高。折中做法是基础配置点用拉丁超立方生成在初始和边界附近额外加密一批点每训练若干轮重新采样一次。配置点数量也要注意。一维问题我通常用 10000 到 30000 个点二维问题就需要更多。如果损失一直降不下去可以先看看是不是配置点数量不够或者采样分布明显失衡。3. 实操全流程用PINN求解一维非稳态渗流方程3.1 问题定义与参考解设置这次我做的案例是典型的一维承压含水层水头扩散问题。假设一个含水层长度为 L100m初始水头 H(x,0)10m左右边界分别固定为 H(0,t)8m 和 H(L,t)5m压力传导系数 a0.5 m²/s。这个问题的物理过程很直观边界水头突然降低后高水头区域的水会向低水头区域流动形成一条随时间变化的过渡曲线。因为边界条件非齐次解析解不太好写我用标准隐式有限差分格式计算了一个参考解网格取 500 个空间节点、2000 个时间步这个结果就作为验证 PINN 精度的基准。选择这个案例有几个原因一维问题便于调试方程线性扩散行为稳定即使 PINN 训练有些波动也不会完全发散。如果你第一次接触 PINN建议先用这种简单问题把流程跑通再碰非线性问题。3.2 代码骨架与关键配置这里给出一个基于 PyTorch 的简化实现方便你自己上手。核心部分包括网络定义、损失函数计算和训练循环。import torch import torch.nn as nn class PINN(nn.Module): def __init__(self): super().__init__() self.net nn.Sequential( nn.Linear(2, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 1) ) def forward(self, x, t): u torch.cat([x, t], dim1) return self.net(u) def pde_residual(model, x, t, a): x.requires_grad_(True) t.requires_grad_(True) H model(x, t) H_t torch.autograd.grad(H, t, grad_outputstorch.ones_like(H), create_graphTrue)[0] H_x torch.autograd.grad(H, x, grad_outputstorch.ones_like(H), create_graphTrue)[0] H_xx torch.autograd.grad(H_x, x, grad_outputstorch.ones_like(H_x), create_graphTrue)[0] return H_t - a * H_xx上面的代码里torch.autograd.grad负责自动微分这就是 PINN 能计算方程残差的关键。训练时我们在区域内采样随机点计算pde_residual在边界采样计算bc_loss在初始时刻采样计算ic_loss三者加权相加后反向传播。优化器方面我推荐两阶段策略。先使用 Adam 优化器训练 30000 到 50000 步学习率从 1e-3 开始配合余弦衰减然后再切换到 L-BFGS 精细化训练。L-BFGS 是二阶方法收敛速度在接近最优解时明显更优但直接用它从头训练很容易陷入局部极值。实际效果是 Adam 把损失降到比较低的平台后L-BFGS 能再往下压一个量级。3.3 训练结果与误差分析我记录了一组比较有代表性的结果。Adam 阶段前 5000 步方程残差损失从初始的 1e-2 量级快速下降到 1e-4 量级到 40000 步左右稳定在 3e-5 左右。切换到 L-BFGS 后再迭代几百步整体损失降到 6e-6。最终预测与有限差分参考解在整个时空域上的相对 L2 误差约为 0.8%。这个误差水平说明在这个简单问题上 PINN 已经可以做到和传统数值方法接近的精度。需要说明的是误差大小严重依赖配置点数量、网络宽度和训练轮数。如果把网络从 3 层 64 神经元改成更深更宽的结构精度还能提升但训练时间也会增加。对于实际工程问题不要把目标定在“损失无限接近于零”而是要把误差控制在实际应用可接受范围内。我顺手测了模型推理速度。训练完成后给定任意一组 (x,t)一次前向计算只需要毫秒级时间。在需要批量场景预测时PINN 的这个优势会非常明显。4. 收敛困难、边界震荡与超参数调节我踩过的坑和解决手法4.1 损失降不下去配置点采样问题首当其冲我第一次跑这个渗流方程时方程残差损失降到 1e-4 就怎么也降不下去了。检查代码没有明显问题后来发现是配置点分布太集中在区域中间边界附近几乎没有采样点。边界条件残差虽然看起来已经很低但方程残差在边界附近依然很大整体损失被这些局部区域拖住。解决办法很简单在边界和初始时刻附近增加配置点密度。我当时的做法是把区域内部点、边界点、初始时刻点按 6:2:2 的比例采样每 2000 步重新生成一次点。损失很快就突破了平台期。另一个常见问题是配置点数量太少。一维问题至少需要 10000 个以上低于这个数量网络很容易过拟合到那些稀疏点上无法捕捉整体物理场。4.2 边界预测偏差大软约束和硬约束怎么选PINN 对边界条件的处理有两种形式。软约束就是把边界误差作为损失函数的一项让网络尽量满足硬约束则是通过修改网络输出的形式让边界条件被自动满足不需要参与损失最小化。软约束实现简单但存在一个问题边界点和内部点之间没有绝对的物理隔离网络在拟合边界的同时可能会扭曲附近区域的解。尤其在初始条件和边界条件在角点不连续时比如本例中初始时刻为 10m而左边界突然降到 8m网络需要在角点同时满足两个互相矛盾的条件边界附近出现震荡是常事。我试过的有效做法有几种。第一种是修改初始条件用一个非常短的时间过渡函数让初始水头从 10m 平滑变化到边界值消除角点不连续。第二种是增大边界条件在损失函数中的权重。第三种是改用硬约束构造输出形式比如把网络输出乘上一个在边界处为零的距离函数再加上边界插值项这样边界条件自动满足。综合来看硬约束更稳定但不够通用软约束灵活但需要调权重。我的个人习惯是先软约束跑通如果发现边界震荡无法接受再针对具体边界类型构造硬约束。4.3 训练不稳定的其他常见原因除了采样和边界问题还有几个因素经常导致 PINN 训练不稳定。学习率过大是最容易发现的。Adam 初始学习率超过 5e-3 时损失曲线经常出现突然上升这是因为梯度更新跨过了最优区域。我建议从 1e-3 起步观察前 1000 步损失变化再调整。损失权重失衡是更隐蔽的问题。前面提到初始条件、边界条件、方程残差的量级可能差好几个数量级。如果某一项权重太小网络会直接忽略它。我有个习惯每隔一段时间打印每一项损失的值肉眼观察哪些项下降慢、哪些项已经降到很低然后手动调整权重。还有一种情况是特征频谱偏差。网络对高频信号的拟合能力天然弱于低频信号而渗流方程的解有时包含快速变化的梯度区域。这个问题可以通过 Fourier Feature Embedding 缓解简单说就是把输入 x、t 先映射到一组不同频率的正弦余弦特征上再送入网络。这个技巧在几何或材料参数变化剧烈的渗流问题中效果很明显。5. 从一维到复杂渗流场景这些扩展方向值得关注5.1 非线性Richards方程与非饱和渗流基本渗流方程假设含水层是饱和的、渗透系数是常数但在实际工程中地下水位以上的非饱和区渗流更常见。这时控制方程变成 Richards 方程C(h) ∂h/∂t ∇·(K(h) ∇(h z))这里的 h 是压力水头C(h) 是容水度K(h) 是非饱和渗透系数都是关于 h 的强非线性函数。传统数值方法处理这类问题时要特别注意迭代收敛和水量平衡误差参数曲线的不光滑还会加剧数值震荡。PINN 处理非线性方程的方式和线性方程没有本质区别只需要把 C(h) 和 K(h) 这些非线性系数在残差计算时求出来代入即可。关键难点在于残差函数变得更加复杂梯度传播路径变长训练时更容易出现梯度消失或发散。我的建议是先用一个简单的 Brooks-Corey 或 van Genuchten 模型测试固定模型参数只求解压力水头场跑稳之后再考虑参数反演。5.2 多相流、裂隙介质与参数反演在石油工程和地下水污染修复中经常遇到油水两相流或污染物迁移问题。这时候需要求解一套耦合方程组包括压力方程和饱和度方程。传统模拟器大多使用隐式压力显式饱和度方法时间步长受限而且网格需要随饱和度锋面调整。PINN 扩展到这个场景并不难因为损失函数可以同时包含多个方程。网络输出不再是单一水头而是多个输出节点比如压力、饱和度和浓度每个输出都对应一组偏微分方程残差。不过多输出会显著增加训练复杂度损失项之间的权重平衡更难。裂缝介质是另一个有潜力的方向。离散裂缝网络在传统有限元中需要非常精细的网格模拟裂缝单元而 PINN 可以用带符号距离函数描述裂缝位置把裂缝处的流量交换作为内部约束加入损失避免了复杂的网格剖分。参数反演这块非常值得聊。传统方法要通过多次正演迭代来配合目标函数优化PINN 则可以让渗透系数场直接成为一个可训练的参数场用神经网络或者通常的分段常数函数来表示 K(x) 和 S(x)和网络的权重一起训练。只要损失函数中加入观测井的水头数据训练完成后就直接得到参数反演结果。我见过一个案例在某场地的示踪试验中用 PINN 联合水头和浓度观测数据反演渗透系数效果比传统配线法更稳定而且还能够输出整个空间上的参数分布。5.3 PINN与传统数值方法的关系我需要强调一点PINN 不是来取代传统数值方法的。成熟商业软件和科研代码里的有限元、有限体积方法精度高、稳定性好经过了大量工程验证这些优点 PINN 暂时比不过。PINN 真正的定位是传统方法的互补工具。在需要快速预测、参数反演、数据同化、缺少网格或网格生成成本过高的场景里PINN 能提供传统方法不容易给出的解决方案。反过来传统数值方法生成的参考解又可以为 PINN 提供训练数据或验证标准。两者完全可以共存而不是谁替代谁的关系。我在这次一维渗流求解过程中也始终保留有限差分参考解作为校验。每训练完一轮都要对比预测解和参考解的误差分布。如果没有参考解我很难判断网络学到的是不是合理的物理场还是仅仅在数字上降低了损失。最后分享一个我个人的体会。PINN 这种物理约束神经网络看起来只是把偏微分方程写进了损失函数但实际工程应用中远远没那么简单。配置点怎么布、损失权重怎么配、边界条件怎么处理这些细节直接决定了模型能不能用。我踩过的坑基本都是这些环节。如果你也想在自己的渗流或偏微分方程问题上试一试我的建议是先抄一个最简单的线性扩散方程跑通全流程再逐步增加非线性和复杂度。别一上来就网撒得太大不然出了问题根本分不清是代码问题、数据问题还是网络结构问题。还有一个很实用的小技巧训练过程中每隔一段时间就把网络在固定测试点上的预测值打印出来直接观察物理曲线变化比单纯盯着损失值盲目调参要靠谱得多。
