news 2026/10/20 4:13:05

用PINN求解基本渗流方程:从原理到工程实践

作者头像

张小明

前端开发工程师

1.2k 24
文章封面图
用PINN求解基本渗流方程:从原理到工程实践

1. 这个标题到底在说什么:PINN 和渗流方程为什么绑在一起

如果你去翻计算力学或者油藏数值模拟相关的资料,最近几年一定会频繁撞见一个缩写:PINN,全称 Physics-Informed Neural Network,中文一般叫物理约束的神经网络。这个标题里写得很明确,它要干的事情不是普通的“数据拟合”,而是用神经网络去解偏微分方程,并且拿一个非常具体的对象开刀——基本渗流方程。

先别被“偏微分方程”这几个字吓退。说白了,渗流方程描述的就是流体在多孔介质里的流动规律,比如地下水往井里流、油藏在注水开发时的水驱前缘推进,都属于这类问题。传统的解法是有限差分、有限元那一套,把求解域剖分成网格,然后迭代求解代数方程组。这个方法成熟可靠,但有个绕不开的痛点:一旦计算域形状复杂、边界条件不规则,或者网格质量差,前处理的工作量就能让人崩溃。

PINN 的思路完全不一样。它把物理方程本身当作一种“约束”,直接塞进神经网络的损失函数里。网络负责输出一个候选解,损失函数负责告诉这个候选解符不符合物理规律。训练完成后,你就得到了一个可以连续输出任意位置解的函数,根本不需要网格。这个思路本质上是用“机器学习的训练过程”替代“传统数值方法的矩阵求解过程”。

那这个标题里的“基本渗流方程”是什么水平?它是渗流力学里最入门的一类控制方程,通常指达西定律和质量守恒联立后得到的压力控制方程。稳态情况下就是一个椭圆型方程,非稳态情况下加一个时间项变成抛物型方程。用 PINN 去解它,恰恰是入门 PINN 的最佳路径:方程形式简单、物理意义明确、解析解或者传统数值解容易拿到做对比,可以让你把注意力集中在方法本身,而不是被复杂的偏微分方程困住。

这篇文章我就围绕这个标题展开,从 PINN 的原理讲起,再给出一套完整的、可以照着复现的渗流方程求解流程,最后聊一聊我实际跑实验时遇到的那些坑。适合谁看?对 PINN 感兴趣但还没动手写过代码的人,以及已经在传统数值模拟里深耕、想看看神经网络这条路靠不靠谱的人。

2. PINN 的核心机制:为什么物理方程能变成神经网络的损失函数

2.1 从“数据拟合”到“物理约束”,这一步跨越了什么

传统神经网络能干什么?给它一堆输入输出样本,它学一个映射关系。比如你给它某个位置的坐标,让它输出该位置的压力值,这本质上是一个回归任务。但这里有个致命问题:如果训练数据不够、采样位置不合理,神经网络学到的映射可能只在样本附近有效,离开样本区域立刻“放飞自我”。

PINN 的出发点就是解决这个泛化问题。它不满足于“只有数据位置的解是对的”,而是要求“整个求解域内,这个解必须满足物理方程”。具体做法是:把偏微分方程改写成一个残差形式,比如方程本来是

[ \frac{\partial p}{\partial t} - \alpha abla^2 p = 0 ]

那残差就定义为网络输出的压力值带入上述方程左边之后剩下的量。如果网络输出是精确解,残差应该是 0;如果偏差很大,残差就很大。这个残差直接加到损失函数里,训练过程就会同时往“匹配已知边界/初始数据”和“满足物理方程”两个方向使劲。

你可以这么理解:传统神经网络的训练样本来自“实验测量”,而 PINN 的训练样本里多了一部分来自“物理定律本身”的虚拟样本。这些虚拟样本不需要真实测量值,只需要你在求解域内撒点,然后计算残差。撒点位置越多,物理约束越强,解就越可靠。

2.2 损失函数的构成:数据项、物理项、边界项如何平衡

一个典型的 PINN 损失函数长这样:

[ L = L_{data} + L_{PDE} + L_{BC/IC} ]

其中:

  • (L_{data}) 是监督项,用有真实标签的位置计算误差,比如某些观测井的压力数据;
  • (L_{PDE}) 是物理残差项,在求解域内部随机采样点,把网络输出代入控制方程,计算残差的均方值;
  • (L_{BC/IC}) 是边界条件和初始条件项,在边界或者初始时刻采样,让网络输出符合给定的约束。

每一项前面都可以加权重系数。实际调试中,我发现权重比选择对收敛速度影响巨大。物理残差项的权重如果太小,网络会偏向拟合数据但忽略物理规律;如果太大,又可能导致边界条件被“稀释”,最终解整体偏移。我常用的做法是让三项的初始量级差不多,然后在训练过程中根据损失变化动态调整。

2.3 自动微分:PINN 能够落地的技术基石

PINN 能跑起来,完全依赖深度学习框架里的自动微分。传统数值方法用差分公式近似导数,PINN 则直接用链式法则精确求出网络输出对输入坐标的导数,精度可以达到机器精度。

具体来说,你定义网络输出为 (p(x, t)),然后框架会自动计算 (\partial p / \partial x)、(\partial^2 p / \partial x^2)、(\partial p / \partial t) 这些导数,你只需要把它们组合成残差表达式就行。这意味着什么?意味着你不用手写复杂的离散格式,也不需要担心差分步长选多少。这个优势在求解高阶偏微分方程或者方程含有交叉导数项时特别明显。

注意:自动微分求导是精确的,但前提是你的网络输出对输入坐标是光滑可导的。这也是为什么 PINN 里激活函数通常选 (\tanh) 或 (\text{Swish}),而不是 ReLU——ReLU 在零点不可导,虽然实际中还没出过大问题,但理论上有隐患。

3. 基本渗流方程的数学形式与物理背景

3.1 达西定律和质量守恒:渗流方程从哪来

要解渗流方程,先得明白它是怎么推导出来的。流体在多孔介质中的流动遵循达西定律,它说的是流速和压力梯度成正比:

[ \mathbf{v} = -\frac{k}{\mu} abla p ]

其中 (k) 是渗透率,(\mu) 是流体黏度,(p) 是压力。这个公式等价于电学里的欧姆定律:压力差相当于电压,流速相当于电流,渗透率和黏度的比值相当于电导。

再结合质量守恒方程:

[ \frac{\partial (\phi \rho)}{\partial t} + abla \cdot (\rho \mathbf{v}) = q ]

假设流体微可压缩、密度变化不大,把达西定律代入质量守恒,就得到了压力控制方程。最常见的形式是:

[ \phi \mu c_t \frac{\partial p}{\partial t} = abla \cdot (k abla p) + q ]

如果进一步假设渗透率在空间均匀,各项参数合并成一个系数 (\eta = k / (\phi \mu c_t)),方程简化为:

[ \frac{\partial p}{\partial t} = \eta abla^2 p + \frac{q}{\phi c_t} ]

这个形式就是经典的抛物型方程,也就是标题里“基本渗流方程”的典型代表。稳态情况下时间项为零,方程退化成拉普拉斯方程或泊松方程:

[ abla^2 p = -\frac{q}{k} ]

3.2 边界条件和初始条件:控制方程必须有“配套约束”

偏微分方程本身只是描述“场内属性如何相互关联”,真正决定一个具体问题解的唯一性的,是边界条件和初始条件。渗流问题里最常见的有三类:

  • 第一类边界条件(狄利克雷条件):直接给定边界上的压力值,比如定压边界,井底压力固定;
  • 第二类边界条件(诺伊曼条件):给定边界上的压力梯度,也就是流速,比如封闭边界,法向流速为零;
  • 初始条件:给定初始时刻整个求解域内的压力分布。

在 PINN 里,这些条件不需要像传统方法那样单独建立离散方程,而是放进损失函数的 (L_{BC/IC}) 项中,通过采样点让网络输出逼近给定值。

3.3 为什么说“基本”二字:从最简模型入手的好处

你可能觉得,既然 PINN 这么厉害,为什么不直接上复杂方程?我觉得恰恰相反,从基本渗流方程入手是最聪明的选择。原因有几点:

第一,基本方程的解析解或半解析解在不少简化几何条件下是存在的,比如一维无限大介质中的压力传播解。你可以拿 PINN 的结果和解析解对比,判断方法本身是否实现正确。如果一上来就解一个复杂的耦合方程组,出了偏差根本分不清是 PINN 的问题还是方程模型的问题。

第二,基本渗流方程是抛物型方程,这个类型在 PINN 的文献里研究得最透彻,训练稳定性和收敛性相对可控,适合新手建立信心。

第三,渗流方程的系数物理意义明确,便于在调试时通过改变参数观察解的变化趋势,帮助建立直觉。比如说渗透率增大,压力扩散应该加快,这个直观判断可以用来检验 PINN 是否学对了物理。

4. 完整实操:用 PINN 求解一维瞬态渗流方程

4.1 问题设定:一维孔隙压力扩散模型

我实际跑的实验,选的是最简单的一维瞬态渗流问题,也就是孔隙压力扩散方程。控制方程:

[ \frac{\partial p}{\partial t} = \eta \frac{\partial^2 p}{\partial x^2} ]

计算域取 (x \in [0, 1]),时间取 (t \in [0, 1]),假设初始条件:

[ p(x, 0) = 0 ]

边界条件:

[ p(0, t) = 1, \quad \frac{\partial p}{\partial x}\bigg|_{x=1} = 0 ]

这个问题的物理场景可以想象成一根一维岩心,左端突然施加一个恒定压力,右端封闭。压力波会从左边向右传播,最终全场的压力都趋于 1。

扩散系数 (\eta) 取 0.5,这个值纯粹是为了让现象在时间窗口内展示得比较明显。系数太大,压力传播过快,后期全场都是稳态,看不出动态过程;系数太小,压力还没传播到右端计算就结束了,误差集中在传播前沿,也不利于观察。

4.2 网络结构与采样策略

网络结构方面,我用的不是特别复杂的网络,就是 4 层全连接网络,每层 50 个神经元,激活函数用 (\tanh)。输入是二维向量 ([x, t]),输出是压力 (p)。

为什么这么设计?因为对于这个简单问题,层数过多会导致训练变慢,甚至出现梯度消失;神经元数过多则容易过拟合到采样点上的噪声。4 层 50 个神经元是文献里比较常见的配置,也是我实测在速度和精度之间比较平衡的配置。

采样策略非常关键。我把整个计算域 ([0,1] \times [0,1]) 当成一个二维区域来处理:

  • 内部残差点:用拉丁超立方抽样取 5000 个点,计算 PDE 残差;
  • 边界点:左边界 (x=0) 和时间轴 (t \in [0,1]) 上取 500 个点,约束压力等于 1;
  • 右边界 (x=1) 上取 500 个点,约束压力梯度为零;
  • 初始时刻 (t=0) 上取 500 个点,约束压力等于 0。

这里有个我后来优化过的细节:初始时刻和左边界交汇处的点会同时具有两个约束,如果处理不好,训练后期容易出现震荡。我的做法是让边界点集合包含这个角点,但它的权重在训练初期调低一点,让网络先学习内部物理规律,再慢慢收紧角点约束。

4.3 损失函数的具体实现与权重设置

用 Python 和深度学习框架实现时,关键是写好 PDE 残差的计算逻辑。网络输出 p_pred 之后,用框架的自动微分分别对 x 和 t 求导,构造残差:

p_x = tf.gradients(p_pred, x)[0] p_xx = tf.gradients(p_x, x)[0] p_t = tf.gradients(p_pred, t)[0] residual = p_t - eta * p_xx

然后整个损失函数就是三部分均方误差的和:

loss_pde = tf.reduce_mean(tf.square(residual)) loss_bc = tf.reduce_mean(tf.square(p_left_bc - 1.0)) + tf.reduce_mean(tf.square(p_right_grad)) loss_ic = tf.reduce_mean(tf.square(p_initial)) loss = loss_pde + 1.0 * loss_bc + 1.0 * loss_ic

你可能会问,为什么权重都取 1?因为我先尝试让不同损失项的量级匹配。如果训练过程中发现边界条件收敛太慢,单独把边界权重调到 5 或 10 就可以解决。调整权重没有固定的万能公式,需要观察损失曲线的相对大小来定。

4.4 训练策略:从 Adam 到 L-BFGS 的两阶段法

训练过程我强烈推荐两阶段法。第一阶段用 Adam 优化器,学习率设为 (10^{-3}),训练大概 2 万步。这个阶段的目的不是精确收敛,而是让网络参数快速进入一个“合理的盆地”,整体压力分布的大致形态是对的,但细节可能还有偏差。

第二阶段换成 L-BFGS 优化器继续训练。L-BFGS 是拟牛顿方法,收敛速度在光滑问题上远快于 Adam,而且不需要手动调学习率。我一般让它跑几百步就能看到损失显著下降。这里有个小技巧:L-BFGS 对初始值敏感,如果直接从头开始用 L-BFGS,很容易收敛到局部极小值;先让 Adam 跑一段,把网络参数带到一个比较合理的区域,再切换,效果会好很多。

关于训练轮数,我的经验值是总计 3 万步以内可以收敛到不错的结果,损失函数可以从初始的 0.1 量级降到 (10^{-5}) 到 (10^{-6}) 量级。再进一步降低就需要调整网络结构或者增加采样密度了。

4.5 结果验证:除了看损失,还要看剖面和渐近行为

训练结束之后,不能只看损失函数数值,必须要做物理层面的验证。我的验证流程有三步:

第一,画出不同时间点的压力剖面。(t=0.1) 时压力波应该刚传播到 (x\approx0.3) 附近,(t=0.5) 时整个左半场应该接近 1,右端压力开始上升,(t=1.0) 时全场应该接近 1。如果这些渐近行为不对,说明网络学到的不是正确物理过程。

第二,和解析解对比。一维瞬态扩散问题在封闭右端条件下有分离变量解,我可以解析地算出精确压力分布,然后和 PINN 预测结果在每个采样点对比。我实测下来,在不靠近传播前沿的位置,相对误差可以控制在 1% 以内;在传播前沿附近误差稍大,因为压力梯度很陡,网络需要更多采样点才能刻画。

第三,做一次网格无关性和采样密度对比。把内部采样点从 5000 增加到 20000,观察结果是否变化明显。如果变化很小,说明结果已经收敛;如果变化很大,说明采样还没足够。

5. 从一维到二维:渗流方程 PINN 的进阶玩法

5.1 二维渗流场和复杂边界,PINN 的优势才真正释放

一维问题只能算热身,PINN 真正的优势要等到二维、三维问题才显现。我后来尝试了二维矩形求解域内的稳态渗流问题,控制方程为:

[ \frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2} = -\frac{q(x,y)}{k} ]

域内有两个“井点”,一个注水井一个采油井。注水井附近压力高,采油井附近压力低,流体从注水井流向采油井。我用 PINN 求解这个压力分布,边界全部封闭,也就是法向流速为零。

这里就体现出传统方法的痛苦了:如果要在这个矩形域里布网格,遇到井点这种奇异性强的源汇项,网格需要局部加密。而 PINN 的做法只是在井点附近增加采样点密度,网络自然会在高梯度区域花更多“注意力”去匹配,代码量几乎不用变。

5.2 源汇项的处理:用点源函数替代加载项

二维问题里有个非常实际的细节:井点在实际油藏里直径很小,相对于整个求解域可以近似看作点源。这个点源在数学上是一个狄拉克函数,直接在残差里实现会有奇异性。

我的处理方式是:不用严格的狄拉克函数,而是在井点周围的高斯函数做平滑近似。具体来说,源项表达式为:

[ q(x,y) = \sum_{i} Q_i \frac{1}{2\pi \sigma^2} \exp\left(-\frac{(x-x_i)^2 + (y-y_i)^2}{2\sigma^2}\right) ]

(\sigma) 取网格特征尺寸的 0.1 倍左右。这样做的好处是避免了残差无穷大带来的数值问题,同时又能在宏观尺度上保持注入/采出的总通量守恒。这个处理在传统有限元里叫“等效源项”,在 PINN 里同样适用。

5.3 井间压力干扰的物理验证

二维算例完成后,我会检查井间中点处的压力是否低于注水井压力且高于采油井压力,同时检查压力等值线是否垂直于封闭边界。如果等值线在边界处是斜着穿过去的,说明诺伊曼边界条件没有被很好满足,需要检查边界采样点密度或者损失权重。

这套验证逻辑非常朴素,但确实能抓住绝大部分实现错误。我记得有一次训练结果看着损失很低,结果一画等值线图,发现封闭边界附近的等值线明显异常,最后查下来是右边界采样点坐标写错了,把 (x=1) 写成了 (x=0)。这种错误在传统数值方法中一眼就能看出来,在 PINN 中如果不画图验证,很容易被低损失值欺骗。

6. 实测避坑指南:我踩过的最典型的六个坑

6.1 采样点不是越多越好,要均匀性与关键区域加密并重

我一开始犯的错误是单纯增加采样点数量,以为点数够了精度自然就高。后来发现,均匀撒 20000 个点,不如在边界和源汇项附近加密采样 10000 个点加 5000 个均匀点。因为 PINN 的残差损失是所有采样点的均值,如果大部分点都位于解比较平滑的区域,网络会倾向于优先拟合这些区域,而忽略高梯度区域。这就是所谓的“损失被平滑区域主导”现象。

6.2 tanh 激活函数的初始化尺度影响巨大

同样一个网络结构,初始化权重标准差从 0.1 换成 0.05,训练结果可能差别很大。(\tanh) 在零附近近似线性,如果初始化权重过大,输入稍微偏移就进入饱和区,梯度接近零,训练停滞。我后来固定使用 Xavier 初始化,并在网络输入做标准化,让 (x) 和 (t) 的取值都落在 ([-1, 1]) 范围附近,训练稳定性明显提升。

6.3 损失函数三项权重不是死的,要配合训练进程调整

我见过很多初学者问“PDE 损失、边界损失、数据损失的权重到底怎么设”。我的回答是:设一个初始值,然后盯着训练曲线调整。如果 PDE 损失下降得很快但边界损失还在高位震荡,说明边界权重太小,网络觉得牺牲边界条件能换来物理方程满足;这时候把边界权重调大 5 倍,一般就能纠正回来。反过来,如果边界损失早就降到很低但 PDE 损失停滞,就说明物理项权重不足。

6.4 靠近传播前沿的误差是常态,不要被单个点的大误差吓到

在一维瞬态案例里,压力传播前沿附近曲率最大,网络预测误差天然偏高。我用解析解算过,前沿附近最大误差可以达到 5%,而其他区域误差不到 1%。如果你把整体 L2 误差当成唯一指标,会误以为网络没有收敛。正确的做法是分区统计误差,或者画出误差分布云图,看大误差集中在哪些区域。只要大误差集中在高梯度区域,且随着采样密度增加逐步下降,就没有大问题。

6.5 时间域的处理有两种方案,不要混用

处理含时间的问题时,有两种主流做法:一种是把时间当成网络的第二个输入,一次训练得到整个时空域的解;另一种是把问题看作某种“参数化初值问题逐步推进”,网络只是一个时间步内空间解的映射。我建议新手从第一种开始,因为实现简单,思路直观。第二种虽然在某些问题上效率更高,但容易引入累计误差,调试复杂度更高,等熟悉了 PINN 的脾性再碰也不迟。

6.6 损失降到多少算“收敛”,要用物理指标说话

最后一条经验:别只看损失函数绝对值。(10^{-6}) 的损失听起来很棒,但你得看这个损失对应的物理量级。压力本身的量级是 0 到 1,损失 (10^{-6}) 相当于均方误差在 (10^{-3}) 左右,也就是说预测压力平均误差大约是千分之一,这已经很好了。但如果你预测的量是渗透率,取值范围从 (10^{-3}) 到 (10^3),同样的损失可能意味着相对误差巨大。所以判断收敛一定要结合你所求解问题的物理尺度。

7. 工程实践中的一点体会

跑了一圈 PINN 解渗流方程之后,我对这个方法的看法是:它不是来取代传统数值模拟的,更像是给传统方法打补丁、扩边界。复杂几何、反问题、参数识别、多物理场耦合这些场景下,PINN 的灵活性和无网格特性确实是亮点。但要说它能在常规工程计算里完全替代成熟软件,目前还早。

以我自己的项目为例,真正有用的场景反而是反问题:已知部分观测压力数据,反推渗透率场分布。传统方法做反演要反复求解正问题,每次都要重新剖网格;PINN 把正问题求解和参数反演放进同一个训练过程中,一次训练同时得到压力场和渗透率场的估计,这个思路在很多场景下效率优势非常明显。

如果有人问我“新手该怎么入门 PINN”,我的建议很直接:拿本篇这个一维渗流方程开始,一步步把代码跑通,然后对比解析解验证。等你确认自己的网络确实“学到了物理”,再逐步增加维度、增加源汇项、换用不规则几何。这条路虽然不算短,但每一步都有明确的验证手段,不容易走偏。

最后再分享一个小技巧:训练结束后,把网络权重保存下来,然后在小范围扰动边界条件的情况下重新微调训练。你会发现 PINN 模型对边界条件变化有一定的“记忆效应”,微调收敛非常快,这在实际工程中对多方案比选很有用。这个特性是传统数值模拟完全不具备的。

版权声明: 本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若内容造成侵权/违法违规/事实不符,请联系邮箱:809451989@qq.com进行投诉反馈,一经查实,立即删除!
网站建设 2026/10/20 4:12:14

C++初阶学习之模板进阶

前言 模板的入门内容通常是「写一个 template <class T> 的函数&#xff0c;让同一个算法适配多种类型」。真正写库的时候才会发现&#xff0c;入门那点知识不够用&#xff1a;参数可以不是类型&#xff0c;类可以针对某类类型单独实现&#xff0c;名字查找分两个阶段&am…

作者头像 李华
网站建设 2026/10/20 4:12:06

多智能体分布式RBF标量场估计:从原理到MATLAB实操

1. 多智能体分布式RBF标量场估计的核心思路拆解1.1 为什么需要多智能体来做标量场估计单智能体做标量场估计这件事&#xff0c;做过的人都知道痛点在哪。假设你有一片区域需要监测某种物理量——温度、浓度、信号强度都行——一个传感器跑一圈&#xff0c;效率低不说&#xff0…

作者头像 李华
网站建设 2026/10/20 4:11:16

四种深度学习模型实战:TensorFlow实现CNN、Transformer、GAN与胶囊网络

上一篇我们完成了TensorFlow和Keras的基础入门&#xff0c;聊了张量、自动微分&#xff0c;以及用Sequential搭一个线性模型的完整流程。但说实话&#xff0c;线性模型和真正能上任务的深度网络之间还隔着一大段路&#xff0c;所以这篇我想直接把四种真正能打的架构拿出来动手过…

作者头像 李华
网站建设 2026/10/20 4:10:47

纯Socket实现西门子S7-1200/1500 PLC数据采集C#源码

简介&#xff1a;本资源是一套面向工业自动化开发者的西门子S7系列PLC数据采集C#源代码&#xff0c;适用于.NET平台下构建轻量级SCADA通信模块或产线数据监控系统&#xff0c;尤其适合具备C#基础与工业通信入门经验的工程师快速集成PLC实时数据读取功能。压缩包共5个文件&#…

作者头像 李华
网站建设 2026/10/20 4:08:57

运算符避坑指南:从算术到位运算,掌握核心原理与常见陷阱

最近在代码评审里看到好几处与运算符相关的低级错误&#xff0c;而且不是只有新人会犯。恰好经常有读者问我&#xff1a;“算术、关系、逻辑、位运算……这些运算符我都认识&#xff0c;可为什么写出来就是不对&#xff1f;”今天就以“运算符大全”为线索&#xff0c;把四类运…

作者头像 李华