0.   引言

地下水渗流问题涉及到石油工程、采矿工程、岩土工程、水利工程和环境工程等多个工程应用领域，其中，在给定渗流参数的情况下精确求解达西流速场（正向问题）和基于渗流实测数据合理识别渗流参数（反向问题）对于渗流模拟分析具有重要意义^[1-3]。

对于正向问题，鉴于物理试验的局限性，达西渗流流速的求解主要采用数值模拟，具体方式分为两种^[2]：①首先对地下水水头进行分析，再基于达西定律求解渗流流速^{[1-2, 5]}；②建立达西渗流定律和渗流模型的耦合方程，通过联合求解同时得到渗流流速和水头^[6-7]。两种方式都已被证明能够进行准确的流速计算，但却难以依据渗流数据对渗流参数进行反演。关于渗流参数的反演，目前主要基于水头或流速数据，采用最速下降法、高斯-牛顿法、遗传算法等优化算法^[8]进行求解，此类方法虽对均质渗流参数具有较好计算表现，但当地下水流为非均质流情况时（渗流参数呈现空间变异性），却难以获得渗流参数的空间分布特点。因此，发展一种既能计算非均质渗流模型的流速场，又能反演其渗流参数的新算法具有重要意义。

近年来，随着人工智能的快速发展，基于深度学习万能逼近定理^[10]建立的物理信息神经网络（Physics-Informed Neural Networks, PINNs）算法受到了学者的广泛关注，其核心原理：利用自动微分技术(Automatic Differentiation, AD)^[11]将偏微分方程嵌入到神经网络的损失函数当中，再通过优化算法对损失函数进行寻优，从而实现神经网络对方程解的近似，具有不受数据量限制、求解格式稳定等特点。此外，相对于有限元法（FEM）、有限差分法（FDM）等，PINNs算法可通过在损失函数中加入实测数据的损失项，对方程的未知物理参数（如渗流方程中的均质或非均质渗流参数）进行反演，甚至可通过较少数据直接挖掘物理规律并构建出相应的物理方程，故在岩土工程研究中具有一定的应用前景。目前，PINNs算法已在一些领域得到了较广泛的应用，如Haghighat等^[12]将PINNs算法引入到固体力学中，基于Python语言的Keras开发了SciANN框架，对带缺陷的平板拉伸变形行为进行分析，得到了平板内各点的位移。基于SciANN框架，唐明健等^[13]对矩形薄板的正反问题进行了分析，准确给出了矩阵薄板各点的挠度，并利用薄板挠度反演得到了边界的受力情况。Raissi等^[14]通过将PINNs算法应用于流体力学当中，准确得到了Burgers方程、退化Navier-Stokes方程的数值解，并实现了方程中流动参数的识别。Chen等^[15]基于实测数据，利用PINNs算法在大量备选的偏微分算子中，寻找合适的算子，实现了Burgers方程的智能建立。PINNs算法在岩土工程领域也有一些探索，Zhang等^[16]引入了独热编码向量，与PINNs算法相结合，对多孔介质的大变形固结问题进行了分析。Bekele等^[17]采用PINNs研究了多孔介质中的流体流动与固体变形之间的相互作用。兰鹏等^[18]采用PINNs算法对连续排水边界的一维固结问题进行了分析，得到了土体的超孔隙水压力消散规律，并实现了连续排水边界界面参数的准确反演。然而，上述的PINNs算法的应用存在如下限制，一方面，PINNs算法在进行反问题分析时，学者采用的是单物理场神经网络模型，此模型只能对未知参数为定值的情况进行反演，当反演的参数呈现空间变异性时，则无法进行识别。因此，对于非均质地下水渗流问题，采用单物理场神经网络模型PINNs算法无法获得非均质渗流参数分布，此时需引入多物理场神经网络方法^[19-22]；另一方面，上述应用中，PINNs算法在进行求解时，对于边界的处理主要是软边界形式。

事实上，PINNs算法软硬约束边界的设置会对计算结果产生较为显著的影响，原因在于PINNs算法在进行求解时，需要构造方程、初始和边界条件的损失函数。初始和边界条件的损失函数主要是通过初始和边界上随机选取样本点及其取值获得，所以在进行优化时，只能保证在选取的样本点处拟合结果较好，无法保证在初始和边界上均满足条件，故而属于一种软约束形式。针对软约束这一问题，Lu等^[23]和Sun等^[24]引入辅助函数对其进行修正，使PINNs计算的结果自动满足边界或者初始条件，此时边界约束形式被称为硬约束形式，相比于软约束，具有更高的计算精度。目前，PINNs算法对于硬约束的应用较少，主要集中在如下方面，Lagari等^[25]总结了不同类型的边界条件的PINNs硬约束构造方法，随后，Lu等^[26]和陆至彬等^[27]分别将硬约束引入到了超材料设计和平面板导热问题分析中，结果表明，硬约束可以获得更高的求解精度。但是在地下水渗流领域，PINNs硬约束应用暂未形成一定的研究规模。

本文将基于PINNs算法，通过施加硬约束边界条件，解决渗流模型的渗流流速计算和渗流参数反演问题。在正向问题中，对比了两种PINNs方式求解渗流流速（渗流方程和达西定律耦合，先求解水头再进行自动微分）的计算精度；在反向问题中，分别构造了单物理场或多物理场神经网络模型对均质或非均质渗流参数进行反演。最后，为了验证PINNs算法在地下水渗流问题应用的可行性，对均质和非均质渗流的若干算例进行分析，表明硬约束PINNs算法对水头进行自动微分求解渗流流速计算精度更高，且能够准确反演不同分布方式的渗流参数。

1.   原理和算法

1.1   渗流控制方程

对于计算域为$\mathit{\Omega}\subset {\mathbb{R}^n}$的稳定渗流问题，其控制模型为^[4]

式中，$\nabla$为梯度算子，x=（${x_1}, {x_2}, \cdots , {x_n}$）为空间坐标，K为渗流参数张量（对于非均质渗流可表示为$\boldsymbol{x}$的函数），H为地下水水头值，W为源汇项。

相关的边界条件（BC）可以表示为

式中，$\partial \mathit{\Omega}$为计算域的边界。基于达西定律，可得到V($\boldsymbol{x}$)的表达式：

1.2   PINNs算法求解达西流速分布和渗流参数

（1）神经网络

运用PINNs算法进行渗流流速求解时，首先需构造神经网络（常见的为前馈神经网络）对水头或流速做近似估计。前馈神经网络由神经元连接而成，分为输入层、隐藏层、输出层3部分。各神经元输入间通过相应的权重和偏置进行线性连接，进而传输到非线性激活函数σ中得到其输出，具体计算方式如下所示，

式中，$z = [{z_1}, {z_2}, \cdots , {z_p}]$为神经元的输入，g(z)为神经元输出，w_j为输入z_j的所占的权重，b为对该神经元所引入的偏置。一般地，常见的激活函数σ包括Tanh=(e^x-e^-x)/(e^x+e^-x)、ReLU=max{0, x}和Leaky ReLU=max{ax, x}(a为任意常数)等。本文所求解的控制方程，需要计算二阶偏导数，而ReLU和Leaky ReLU的二阶导数为0，因此选择能够无限次可微Tanh作为激活函数。

（2）自动微分（AD）

运用PINNs算法对渗流流速进行求解时，主要利用AD技术计算网络输出函数对输入变量的导数。其中，AD是一种介于符号微分和数值微分的方法，其计算原理是将复杂解析函数分解成一系列初等运算的组合，通过符号微分进行求导，再代入相关数值，通过计算机存储中间结果，最后利用链式法则得到期望导数值^[11]。因此，与数值微分相比，AD能够替代网格尺度的差分运算，避免了数值微分在求解过程中方程离散带来的计算误差，提高了计算精度，且自动微分的方式主要针对神经网络的输出进行，这使梯度的计算是精确的。

目前常见的深度学习框架TensorFlow^[28]和PyTorch^[29]等都内置自动微分模块，本文主要采用TensorFlow2.3进行PINNs计算。

（3）PINNs算法求解达西流速

PINNs算法进行渗流流速求解时，有两种不同的求解策略。方法一是通过PINNs算法求解达西渗流定律和渗流模型的耦合方程，其求解原理如图 1所示，首先构造输入为x，输出为$\hat H(x;\boldsymbol{\theta} )$和$\hat V(x;\boldsymbol{\theta} )$的FNN分别对水头H和流速V做近似估计，其中，$\boldsymbol{\theta}$为神经网络中所有权重和偏置参数的集合。其次，通过AD对式（1），（3）和BC构造如下损失函数MSE：

图  1  PINNs算法求解渗流流速示意图

Figure  1.  Schematic diagram of PINNs algorithms for deriving Darcy velocity

下载: 全尺寸图片 幻灯片

其中，

式中：${\left\|\cdot\right\|_2}$为欧几里得范数，${\tau _\mathrm{P}}$，${\tau _\mathrm{B}}$分别为计算区域内部和边界区域中用于训练神经网络的残差点集合，一般可以通过在定义域内随机抽样^[30]或进行指定位置设置；${\omega _{\text{p}}}$，${\omega _\mathrm{B}}$分别为修正损失项MSE_PDE，MSE_BC的权重。最后，通过Adam^[31]、L-BFGS-B^[32]等优化算法对损失函数进行寻优即可得到最优的神经网络参数$\hat{\boldsymbol{\theta}}$：

同理，方法二则是通过PINNs算法构造输入为x，输出为$\hat H(x;\boldsymbol{\theta} )$的神经网络计算地下水头分布$\hat H(x;\hat {\boldsymbol{\theta}})$，再基于达西定律，利用AD计算求解$\hat V(x;\boldsymbol{\theta} )$。值得一提的是，由于方法二求解的控制方程不同，此时MSE_PDE变为

（4）PINNs算法识别渗流参数

PINNs算法不仅能够实现达西流速模拟，还能根据水头的实测或模拟数据对渗流参数进行反演，其原理如图 2所示。

图  2  PINNs算法反演渗流参数示意图

Figure  2.  Schematic diagram of PINNs algorithms for identifying seepage parameters

下载: 全尺寸图片 幻灯片

在进行均质渗流参数反演时，基于正向水头求解的损失函数及实测点${\boldsymbol{x}_D}$的H值（H(${\boldsymbol{x}_D}$)），加入衡量实测数据与预测值间偏差的MSE_Data损失项及其权重${\omega _D}$：

式中：${\tau _D}$为实测数据组成的集合。同理，通过优化算法对总损失函数进行优化可得到最终的神经网络参数$\hat {\boldsymbol{\theta}}$和均质渗流参数$\hat {\boldsymbol{K}}$。

PINNs算法进行非均质渗流参数反演时，由于K随x变化，需要构造多物理场神经网络对渗流参数进行反演，即构造两个神经网络$\hat H(x;{\boldsymbol{\theta} _{\text{1}}})$和$\hat {\boldsymbol{K}}(x;{\boldsymbol{\theta }_{\text{2}}})$分别对水头$H(\boldsymbol{x})$和$K(\boldsymbol{x})$做近似估计，此时${\boldsymbol{\theta} _{\text{1}}}$和${\boldsymbol{\theta} _{\text{2}}}$为各自神经网络的待优化参数，总的神经网络参数可表示为$\boldsymbol{\theta}= \{ {\boldsymbol{\theta }_{\text{1}}}, {\boldsymbol{\theta} _{\text{2}}}\}$。相应地，MSE_PDE、MSE_BC和MSE_Data分别变为

通过优化算法对MSE优化可得到最终的参数$\hat {\boldsymbol{\theta}}$：

（5）硬约束PINNs算法

PINNs算法进行偏微分方程求解或参数反演时，边界的约束主要是软约束的形式，导致计算精度会受到一定的影响。

为了解决上述问题，可以引入合适的函数A(x)和Z(x)，使PINNs计算的结果

能够自动满足边界条件^[25]，式中，${a^L}(\boldsymbol{x};\boldsymbol{\theta} )$为神经网络的输出值。以第三类边界条件为例，假定偏微分方程满足如下边界条件：

式中：f(x)为偏微分方程的解，x_j为x第j个空间坐标；$\hat B_j^L$，$\hat B_j^R$分别为x_j=a_j和x_j=b_j的左、右边界条件算子；$\lambda$和μ为已知常数。A(x)和Z(x)对于第三类边界条件，可以进行如下表示^[25]：

式中，${\mathcal{M}_j}$为作用在f(x)的算子，其定义为

其中，

2.   数值试验

在本节中，将PINNs算法应用于二维（空间坐标为x，y）均质和非均质无量纲渗流模型（参数量纲均为1）求解中，通过与解析解和实际解进行对比，验证硬约束PINNs算法对于渗流流速的求解（以x方向的渗流流速V_x为例）和渗流参数反演的准确性。采用如下缩写方式：PINNs-H-I、PINNs-H-II分别表示方法一、二的硬约束PINNs算法，PINNs-S-I、PINNs-S-II分别表示方法一、二的软约束PINNs算法。所有的PINNs算法都使用Python编写，修正的损失项权重设置为1，优化算法采用学习率为0.001，迭代步数为20000的Adam算法和L-BFGS-B算法，即使用Adam算法迭代至指定步数后，再通过L-BFGS-B进行优化，这样能够提高计算精度和效率^[18]，其他超参数(PINNs算法运行之前设置的参数值)选择如表 1所示。

表  1  算例所用的计算超参数

Table  1.  Hyperparameters of PINNs algorithms used in all examples of Section 2

求解方法	${\tau _{\text{P}}}$	${\tau _{\text{B}}}$	${\tau _{\text{D}}}$	NN(H或H，Vx)	NN(K)
正向求解	2000	200	—	50×3	—
反向求解	400	200	400	100×6	60×6
注：NN表示神经网络隐藏层结构。

下载: 导出CSV 

| 显示表格

为了对计算结果的精度进行分析，采用L₂相对误差（$\varepsilon$）作为评价指标：

式中，$({x_i}, {y_j})$为计算域内随机分布的点，$u({x_i}, {y_j})$为解析解，$\hat u({x_i}, {y_j})$为PINNs数值计算结果。

2.1   渗流参数为均质的二维稳定流

地下水二维稳定流的水头和V_x的控制方程可由式（1）和式（2）得出，令K_x=K_y=1，求解的计算区域$\mathit{\Omega}$=[0, 1]×[0, 1]，控制方程的解析解H= xy(1-x)(1-y)，控制方程满足的Dirichlet条件和源汇项W均由解析解反推得到。

为了进行对比，本例分析了在不同y条件下（y=0.3, 0.5, 0.9），PINNs-H-II、PINNs-S-II、PINNs-H-I、PINNs-S-I和FEM的渗流流速求解精度，其中FEM采用计算域为1024个单元的三角形剖分方法求解H，再基于达西定律求解V_x，如图 3所示。可以看出，4种PINNs计算方式都具有较高的计算稳定性，这表明在不同分布的${\tau _{\text{P}}}$、${\tau _{\text{B}}}$下（随机抽样设置），利用优化算法进行神经网络参数寻优时，能够迭代至较优的计算结果；FEM的计算误差（主要来自于划分网格时的离散误差，以及达西渗透流速在节点处不连续导致的误差^[2]）高于PINNs-H-II产生的误差（由近似误差、泛化误差和优化误差共同决定^[23]），低于其他3种PINNs算法所产生的误差，表明PINNs-H-II可替代FEM作为精确求解地下水达西渗流流速的计算方法；同时，还注意到方法一施加硬约束对计算性能提升不明显，且计算精度低于PINNs-H-II，其原因在于：方法一的MSE_PDE是通过对水头和流速控制方程同时定义的，硬约束只对水头施加，因此无法影响流速计算结果，而方法二的MSE_PDE仅由水头控制方程构造，优化过程硬约束能够对使优化算法更易迭代至全局最优，并且H进行自动微分运算得到V_x的过程是数值准确的，因此优化算法对损失函数进行优化时，方法一相较于方法二会产生更大的优化误差，且施加硬约束对计算性能提升不明显。图 4即为采用PINNs-H-II求解V_x的计算结果，其中图 4（a）为V_x的求解结果，图 4（b）为求解结果与解析解的绝对误差，可以看出，采用自动微分技术能够得到连续光滑且数值准确的数值计算结果，且绝对误差最大为6.8×10^-6，具有较高的计算精度。

图  3  不同位置y下，PINNs-H-II、PINNs-S-II、PINNs-H-I、PINNs-S-I和FEM所计算V_x的L₂相对误差图

Figure  3.  Relative errors of L₂ given by PINNs-H-II, PINNs-S-II, PINNs-H-I, PINNs-S-I and FEM methods under different values of y

下载: 全尺寸图片 幻灯片

图  4  PINNs-H-II求解均质渗流流速V_x的计算结果图

Figure  4.  Values of Darcy velocity V_x of homogeneous seepage and absolute errors given by PINNs-H-II method

下载: 全尺寸图片 幻灯片

为了验证PINNs算法在进行渗流参数反演的准确性，假定K_x=K_y=1是需要反演的量，并设定残差点${\tau _{\text{P}}} = {\tau _{\text{D}}} = 400$（由x，y方向上均匀选取20个点生成），实测值则根据解析解得到。在进行计算时，在使用方法二进行正向求解所构造损失函数的基础上增加一项关于实测值的损失函数，利用优化算法对损失函数进行全局寻优，最终反演结果如图 5所示。可以看出，随着优化迭代步数的增加，PINNs算法施加硬约束后能够更快收敛至准确值，而最终20000步后软约束反演结果为0.97，需要更多的迭代步数才能收敛至准确值，这表明，通过施加硬约束，PINNs算法能够提高对均质二维稳定流的渗流参数反演精度。

图  5  均质情况下反演渗流参数随迭代步数变化情况

Figure  5.  Variation of homogeneous seepage parameter with iteration step numbers in backward problem

下载: 全尺寸图片 幻灯片

2.2   渗流参数渐变的二维非均质稳定流

本例控制方程可由式（1）和式（2）得出，令K_x=K_y=$\frac{{20}}{{2(1 - x)(1 - y) + 1}}$，求解的计算区域$\mathit{\Omega}$=[0, 1]×[0, 1]，控制方程的解析解H=xy(1-x)(1-y)，控制方程满足的Dirichlet条件和源汇项W均由解析解反推得到。

图 6展示了基于两种硬约束PINNs算法的渗流流速的求解结果，其中图 6（a），（b）分别为PINNs-H-I、PINNs-H-II与解析解之间的绝对误差，与2.1节中得到的结果类似，PINNs-H-II求解结果的绝对误差比PINNs-H-I低了约3个数量级，表明在非均质渗流条件下，PINNs-H-II的计算精度高于PINNs-H-I。

图  6  PINNs算法求解非均质渗流流速V_x与解析解的绝对误差云图

Figure  6.  Absolute errors of V_x of non-homogeneous seepage given by PINNs-H-I and PINNs-H-II methods

下载: 全尺寸图片 幻灯片

为了进一步说明软硬约束在边界的作用情况，对比了在y=0下，PINNs-H-II和PINNs-S-II的x与H的计算关系，如图 7所示。可以看出，硬约束PINNs算法的水头计算结果在边界处恒等于0，这与实际水头相吻合，而软约束PINNs算法在边界计算结果存在振荡现象，因此可说明硬约束PINNs算法能够对软约束边界条件进行修正，具有更高的计算精度。

图  7  y=0时，PINNs-H-II和PINNs-S-II所计算的x与H的关系

Figure  7.  Relationship between x and H via PINNs-H-II and PINNs-S-II methods at y=0

下载: 全尺寸图片 幻灯片

对于渗流参数的反演，在方法二求解水头问题的基础上，假定K_x=K_y=$\frac{{20}}{{2(1 - x)(1 - y) + 1}}$为待反演函数，构建多物理场神经网络模型，求解结果如图 8所示，图 8（a）为反演的渗流参数，图 8（b）为反演的渗流参数与实际解的绝对误差云图。可以看出，最大绝对误差在4.5×10^-2附近，因此可证明硬约束PINNs算法反演具有较高的计算精度，并且从图 8（c）可知，在y=0.5处，硬约束PINNs算法能够准确反演出渗流参数随空间位置x的变化关系，而软约束PINNs算法却无法获得准确的求解结果，说明通过施加硬约束，能够提高反演的准确性。

图  8  非均质情况下渗流参数反演结果

Figure  8.  Identified seepage parameters and corresponding absolute errors via soft- and hard-constraint PINNs in non-homogeneous seepage

下载: 全尺寸图片 幻灯片

3.   结论

传统的数值算法主要是通过划分网格的方式进行地下水模型的渗流流速计算，较难实现非均质模型的渗流参数反演，本文将PINNs算法应用于地下渗流模型的正向和反向求解当中，通过施加硬约束提高了算法的精度，主要得到以下3点结论。

（1）基于硬约束PINNs算法构造了PINNs-H-I和PINNs-H-II两种渗流流速计算方法，分析表明，两种计算方法都具有较高的求解稳定性，并通过对比可知，PINNs-H-I在进行损失函数构造时，MSE_PDE需要增加求解V的损失项，导致损失函数更为复杂，因此优化算法更难收敛至最优点，而PINNs-H-II仅由渗流方程确定，进行寻优时具有更高的计算精度。

（2）对PINNs算法施加硬约束可自动满足偏微分方程的边界要求，并能够提高对渗流模型求解和渗流参数反演的精度。

（3）针对渗流参数的分布特点，采用硬约束PINNs算法分布构造了单物理场或多物理场的神经网络模型，通过最小化神经网络的损失函数，对均质或非均质渗流参数进行了反演。在算例中，通过与实际的渗流参数进行对比，验证了反演参数值的准确性。