基于物理信息的极端学习机，结合时间推进算法和随机特征模型，用于非饱和地下水流动的模拟

《Journal of Hydrology》：Physics-informed extreme learning machine with time-marching and random features for unsaturated groundwater flow

　　报告作者：Biao Yuan, Roberto Cudmani, Wei Yan, He Yang, Ana Heitor, Yanjie Song, Xiaohui Chen
所属机构：利兹大学土木工程学院几何建模与人工智能中心，英国利兹LS2 9JT

摘要
机器学习（ML）为工程中的数值建模提供了新的机会，但纯数据驱动的方法需要庞大的数据集且缺乏物理可解释性，而基于物理信息的神经网络（PINNs）由于其基于梯度的优化算法通常计算成本较高。本研究提出了一种基于时间推进策略和随机特征方法（TMS-RFM）的物理信息极端学习机（PIELM），用于高效解决高度非线性的非饱和地下水流动问题，即TR-PIELM。TR-PIELM建立在浅层神经网络和物理信息机器学习（PIML）的基础上，将非饱和渗透的控制方程和初始边界条件编码到损失函数中，实现了软约束学习和硬边界约束的结合。基准测试包括一维稳态/瞬态流动、双层土壤渗透、耦合非饱和渗透与变形、边坡稳定性分析以及二维瞬态流动，结果表明TR-PIELM能够准确捕捉压力头动态。与解析解或有限差分解相比，TR-PIELM在稳态问题上的相对L2误差降低了10^10倍，在瞬态问题上的误差降低了10^2倍，并且其精度可与四阶龙格-库塔（RK4）方案相当甚至更好，同时计算时间显著减少。此外，本文还讨论了影响模型性能的关键因素及潜在的改进方向。该研究强调了该方法在工程中处理高度非线性偏微分方程（PDEs）的适用性，为模拟非饱和渗流提供了一个快速、可靠且可扩展的框架。

1. 引言
非饱和地下水流动的研究在岩土工程和环境工程中的渗流分析及边坡稳定性评估方面具有重要的科学和工程价值（Iverson, 2000; Jiang and Huang, 2016; Wu et al., 2016a），为地质灾害预防和地下水资源管理提供了重要指导（Perani? et al., 2018; Wu et al., 2019）。多孔介质中的非饱和渗透是一个复杂的水土工程过程，通常由理查兹方程（Richards Equation, RE）描述（Richards, 1931）。然而，由于RE的强非线性和复杂的系数，解析解很少能得到（Srivastava & Yeh, 1991）；因此，通常采用有限差分法（Celia et al., 1990; LeVeque, 2007）、有限元法（Huyakorn et al., 1986; Reddy, 1993）、有限体积法（Eymard et al., 2000）和迭代方案（Paniconi & Putti, 1994）来模拟非饱和渗流行为（Chávez-Negrete et al., 2018; Ku et al., 2018; Liu et al., 2015; Su et al., 2022; Zhu et al., 2022b）。当土壤具有变形性时，地下水流动与土壤变形之间的耦合会显著改变非饱和土壤的水力及力学性质（Kakogiannou et al., 2016; Wu et al., 2020）。因此，开发和求解描述渗流与土壤响应之间相互作用的控制方程对于理解和预测非饱和流动行为至关重要。

传统上，上述数值方法在求解控制方程方面取得了相当大的成功。然而，这些方法在处理非线性偏微分方程（PDEs）时仍存在局限性，如离散化误差、高网格成本、数值振荡、昂贵的迭代计算和收敛困难。包括高斯-赛德尔（Gauss–Seidel）、雅可比（Jacobi）、皮卡德（Picard）和牛顿（Newton）方法在内的迭代方案也面临类似挑战（Brenner and Cancès, 2017; Farthing and Ogden, 2017; Paniconi and Putti, 1994; Zhu et al., 2019）。尽管在解决这些问题上取得了一些进展，但仍需要进一步的突破。机器学习（ML）替代模型作为一种有前景的解决方案出现。近年来，ML通过强大的数据拟合能力为岩土工程问题提供了新的解决方案（Ebrahimi et al., 2014; Yuan et al., 2025a）。然而，纯数据驱动的ML模型需要大规模的高质量数据集，且往往缺乏物理可解释性（Sun et al., 2020）。作为突破性的替代方案，物理信息神经网络（PINNs）已被开发并在多个领域成功应用（Guo et al., 2025; Wong et al., xxxx; Yuan et al., 2024）。通过将PDE的残差嵌入损失函数并通过梯度下降进行优化，PINNs可以独立于网格大小或时间步长来近似各种PDE，从而实现灵活的无网格模拟。在岩土工程领域，PINN方法展示了显著的应用潜力（Sheil et al., 2026; Yuan et al., 2025a），在固结分析、基础沉降预测及相关逆问题中取得了有希望的结果（Bekele, 2021; Lan et al., 2024a; Lu and Mei, 2022; Yuan et al., 2025b; Zhang et al., 2022）。此外，PINN能够有效整合实验观测数据或原位监测数据，用于建模和分析复杂的地质力学行为（Wang et al., 2024; Zhang et al., 2023）。该方法在模拟地下水流动和研究多物理耦合过程方面表现出强大的适应性和优势（Ali et al., 2024; Amini et al., 2022; Depina et al., 2022; Tartakovsky et al., 2020; Wang et al., 2023b），例如水力-力学相互作用（Li et al., 2024; Ng et al., 2025）。Lan等人（Lan et al., 2024b）证明了一种改进的PINN变体PECANN是模拟非饱和流动中理查兹方程的有效工具。Bandai和Ghezzehei（Bandai & Ghezzehei, 2022）使用域分解PINN框架来模拟具有不连续水力传导性的异质非饱和土壤中的水流。Kamil等人（Kamil et al., 2024; Kamil et al., 2025）开发了用于模拟非饱和土壤中溶质传输和耦合过程的物理信息深度学习方法。Haruzi和Moreno（Haruzi & Moreno, 2023）将地球电监测数据与PINNs结合，用于模拟非饱和介质中的水流和溶质传输。这些研究凸显了将PINNs应用于地下水流建模的兴趣日益增加，同时也指出了在计算效率和预测精度方面仍存在的挑战。然而，值得注意的是，PINNs通常需要较长的训练时间，并且相比传统数值方案仅带来微小的精度提升。此外，传统数值方法面临的局限性同样适用于PINNs，这大大限制了它们的计算效率和泛化能力。

为了解决这些局限性，极端学习机（ELM）为求解PDEs提供了一种新的范式（Guang-Bin et al., 2004; Huang et al., 2006）。它采用具有固定内部参数的单隐藏层随机神经网络，提供快速训练、强泛化能力和最小的超参数调整需求（Dong & Yang, 2022）。与深度神经网络相比，这种简化的架构大大降低了训练复杂性，同时保留了表示时空依赖性的能力（De Ryck et al., 2025）。物理信息极端学习机（PIELM）通过惩罚项将控制方程和初始边界条件嵌入损失函数（Dwivedi & Srinivasan, 2020）。通过随机特征映射来近似解空间（Chen et al., 2022a; Rahimi and Recht, 2007），PIELM仅需要通过线性或非线性最小二乘法优化输出层系数，从而避免了基于梯度的方法中常见的局部最小值问题（Dong & Li, 2021）。此外，已经严格证明随机特征函数具有通用逼近性质，其凸优化特性确保了训练过程中的数值稳定性（De Ryck et al., 2025）。

为了提高PIELM处理复杂边界条件的性能，一些学者提出了一种将域分解策略与统一分割（PoU）方法相结合的新方法（Chen et al., 2022a; Chen et al., 2024; Dong and Li, 2021）。整个计算域被划分为几个子区域，在每个子区域内构建独立的局部随机神经网络。在子域接口处施加连续性条件以组装全局一致的解。所有局部网络的隐藏层参数保持不变，而输出系数通过最小二乘法进行优化，有效地将原始的高维非线性问题转化为可解的线性系统，从而提高计算效率和收敛稳定性。在此基础上，Chen等人（Chen & Luo, 2023）引入了一种块时间推进策略，以减轻长期模拟过程中的误差累积。在传统数值理论中，通常使用隐式欧拉方法等时间积分方案。例如，Chávez-Negrete等人（Chávez-Negrete et al., 2018）开发了一种结合自适应步长Crank–Nicolson方案的广义有限差分方法，以提高理查兹方程解的精度。Abdellatif等人（Abdellatif et al., 2018）提出了一种高精度离散化技术，将隐式欧拉方法在时间上与空间上的谱方法结合。Radu等人（Radu et al., 2018）进一步将向后欧拉时间离散化与混合有限元空间离散化结合，以更准确地模拟非饱和土壤中的两相流动。此外，可变时间步进策略可以进一步提高皮卡德迭代的计算效率（Crevoisier et al., 2009）。

为了研究基于时间推进策略和随机特征方法（TMS-RFM）的物理信息极端学习机（PIELM，即TR-PIELM）在岩土工程中模拟非饱和土壤地下水渗透的适用性，本研究重点关注了几种代表性的非饱和渗透模型，包括：一维（1D）稳态渗透、1D瞬态渗透、通过双层土壤系统的一维渗透、带变形的一维渗透、一维瞬态坡面流动以及二维（2D）瞬态渗透。对解决方案过程进行了全面分析，并通过与参考解的比较验证了模型的性能。比较结果证实了所提出方法在解决岩土工程中的非饱和渗透问题方面的有效性和可行性。本文结构如下：第2节概述了非饱和流动的理论基础及结合时间推进和随机特征方法的TR-PIELM框架。第3节展示了六个基准非饱和流动案例的结果。第4节提供了评估该方法性能的敏感性分析。第5节讨论了时间推进策略的变异性。最后，第6节总结了讨论和结论。

2. 理论和方法
2.1. 控制方程
理查兹方程可用于描述土壤中水的运动。因此，非饱和土壤中的一维地下水渗透过程可以表示如下（Liu et al., 2015; Richards, 1931; Zha et al., 2017），受给定初始和边界条件的约束：
(1)
??zKzψ
?ψ
?z + ?Kzψ
?ψ
?z = ?θ
?t,
z ∈ [0, L],
t ∈ [0, T]
Iψ,
z, t = 0 = 0,
z ∈ [0, L]
Bψ,
z, t = 0,
z ∈ ?Ω,
t ∈ [0, T]
其中，ψ表示压力头，z表示垂直和流动方向，t表示时间。Kzψ表示z方向上的非饱和水力传导率，θ表示体积含水量，L表示土壤层厚度。I和B分别表示初始和边界条件，可以采用狄利克雷（Dirichlet）、诺伊曼（Neumann）、罗宾（Robin）或周期性条件。假设水力传导率Kzψ和含水量θ随压力头变化，它们的相互依赖性可以通过本构关系来描述。一般来说，θ和Kzψ的变化可以用经验指数模型表示（Gardner, 1958）：
(2)
θψ = θr + (θs - θr)e^αψ
(3)
Kzψ = Ks e^αψ
其中ψ表示压力头，θs表示饱和含水量，θr表示残余含水量，Ks表示饱和水力传导率。参数α是一个与土壤性质密切相关的拟合系数，用于表征土壤的孔隙大小分布（Gardner, 1958）。为了便于求解方程（1），通过引入新参数ψ*来线性化控制方程：
(4)
ψ* = e^αψ - e^αψd
其中ψd表示干土的压力头。将方程（2）、（3）和（4）代入控制方程（1）得到：
(5)
Kα?2ψ*
?z2 + Kθ?ψ*
?z = ?ψ*
?t
(6)
ψ = 1/α ln(ψ*) + e^αψd
其中Kθ = Ks / (θs - θr)，Kα = Kθ / α。

2.2. 从PINN到PIELM
2.2.1. 物理信息神经网络
物理信息神经网络（PINNs）是一种基于神经网络的PDE求解方法。它们遵循基本的物理定律并给出控制PDE的近似解（Guo et al., 2025）（见图1）。

在本研究中，考虑了一个没有标记数据的PINN模型。根据通用逼近定理，基于深度神经网络的PINN可以逼近任何连续函数，如下所示：
(7)
ψ^(z, t; W, b) = ∑?=1^m W?·σ·…·σW?·σW?·z, t
T + b? + b?·…·+ b?
其中，z, t表示输入，σ是激活函数，W和b表示网络的权重和偏置。PINN由m层组成。为了获得近似解ψ^，PINN将方程（1）中的三个方程嵌入损失函数。具体来说，根据非饱和土壤中的一维地下水渗透模型的物理模型，初始条件损失、边界条件损失以及基于控制方程残差的PDE损失可以通过均方误差（MSE）损失来定义：
(8)
Lic = 1/Nic ∑?=1^Nc(Iψ^(z, t; W, b), z, t = 0
(9)
Lbc = 1/Nbc ∑?=1^Nbc(Bψ^(z, t; W, b), z ∈ ?Ω, t)
(10)
Lpde = 1/Npde ∑?=1^Npde(??z Kzψ^(z, t; W, b)
?ψ^(z, t; W, b)
?z + ?Kzψ^(z, t; W, b)
?z - ?θ(ψ^(z, t; W, b))
?t
在本研究中，ψ^和ψ分别表示预测值和参考值。Nic、Nbc 和 Npde 分别表示从初始条件、边界条件和计算域中的控制偏微分方程（PDE）中随机抽取的样本点数量。最后，所有三个损失项都使用均方误差（MSE）进行计算，并分别表示为 Lic、Lbc 和 Lpde。每个组成部分都被分配了一个相应的权重，总损失函数被定义为这三个项的加权和：(11)L=ωic·Lic+ωbc·Lbc+ωpde·Lpd，其中 ωic、ωbc 和 ωpde 分别代表三个损失项的权重。根据损失函数方程 (11)，PINN 的最优参数为：(12)(W,b)?=argmin(W,b)L(W,b) PINN 的目标是通过最小化损失函数 L(W,b) 来确定可学习的网络参数 (W,b)，从而获得 PDE 的近似解。参数 (W,b) 通常使用 Xavier 初始化等方法随机初始化。在优化（训练）过程中，当损失函数 L(W,b) 变得足够小或达到预定义的训练迭代次数时，迭代终止。在许多研究中，最小化损失函数的一种常见方法是在每个训练阶段使用梯度下降法迭代更新神经网络参数：(13)(W,b)k+1=(W,b)k-η·?W,bLW,b=(W,b)k-∑j=ic,bc,pdeηωj·?W,bLjW,b，其中 η 表示学习率，是一个可调节的正超参数，?W,bLjW,b 表示损失函数相对于参数 (W,b)k 的梯度，可以通过链式法则计算得出。

2.2.2. 物理信息极端学习机（Physics-Informed Extreme Learning Machine, PIELM）是一种基于神经网络的新兴方法，用于求解一般的偏微分方程（PDE）。它整合了 PINN 和 ELM 的关键概念（Dwivedi & Srinivasan, 2020）。（见图 2）。图 2 展示了用于求解偏微分方程的 PIELM 网络架构示意图。PIELM 通常采用单隐藏层前馈神经网络，隐藏层包含 n 个神经元/节点。因此，域 Ω 中的近似解 ψ^z,t;C 可以表示为线性组合：(14)ψ^z,t;C=∑i=1nci·σiwi·z,tT+bi=ΦC，其中 Ω 是 z∈[0,L], t∈[0,T]，根据控制方程 (1)，激活函数 σi 可以选择为 Tanh、Sine 或高斯函数。权重参数 wi 和偏置参数 bi 从均匀分布中随机初始化并在 PIELM 中保持不变，而 ci 表示待确定的输出层组合系数，C=(c1,c2,?,cn)T。为了训练单隐藏层神经网络并确保最小的训练误差，我们的目标是确定：(15)minC‖Φz,t,wi,bici-ψz,t‖，这里 i=1,2,?,n，‖?‖ 表示欧几里得范数。这等同于最小化以下损失函数：(16)argminCL(w,b,C)=argminC(∑i=1nciσwi·z,tT+bi-ψi(z,t))2。像方程 (16) 这样的优化问题是非凸的，因此解决起来具有挑战性。可以应用传统的基于梯度的算法，如 PINN 中的算法，但它们需要迭代调整所有参数，导致计算时间长、收敛慢，并且由于学习率选择不当可能会陷入局部最小值。相比之下，在 ELM 中，一旦隐藏层权重 w 和偏置 b 被随机初始化并固定，输出层组合系数矩阵 C 就唯一确定了。训练这样的单隐藏层网络简化为求解线性代数系统，输出层系数 C 可以通过最小二乘法唯一获得：(17)ΦC=Ψ，(18)C=Φ?Ψ，这里 Φ? 表示矩阵 Φ 的摩尔-彭若斯伪逆。可以证明，从方程 (18) 获得的解 C 具有最小范数且是唯一的。Ψ 表示目标解矩阵，并受到 PIELM 中物理定律的约束。考虑到方程 (1) 中的控制方程形式，将方程 (14) 代入 ψ 可得到残差函数：(19)Ric(z,t;C)Rbc(z,t;C)Rpde(z,t;C)=Iψ^(z,t;C)，z,t=0Bψ^(z,t;C)，z∈?Ω,t??zKzψ^(z,t;C)?ψ^(z,t;C)?z+?Kzψ^(z,t;C)?z-?θ(ψ^(z,t;C))?t，其中 Ric、Rbc 和 Rpde 分别表示对应于初始条件、边界条件和控制方程的残差项。Nic、Nbc、Npde 是从初始条件、边界条件和解域 Ω 中随机抽取的样本点。因此，配置点的总数为 N=Nic+Nbc+Npde。通过强制方程 (19) 中的残差函数在所有配置点处消失，可以获得解的线性系统。与 PINN 框架类似，也可以为每个损失项分配适当的权重：(20)Rz,t;C=ωic·Ric(z,t;C)ωbc·Rbc(z,t;C)ωpde·Rpde(z,t;C)=ωic·IΦicCωbc·BΦbcCωpde·NΦpdeC=0，其中 I(?)、B(?) 和 N(?) 分别表示初始条件、边界条件和 PDE 算子。Rz,t;C∈RN×1，Ric(z,t;C)∈RNic×1，Rbc(z,t;C)∈RNbc×1，Rpde(z,t;C)∈RNpde×1，Φic∈RNic×n，Φbc∈RNbc×n，Φpde∈RNpde×n。方程 (20) 中的系统是关于 C 的代数系统，包含 N 个方程和 n 个未知数。对于线性偏微分方程，代数系统方程 (20) 变为关于 C 的线性系统，可以根据 (17) 和 (18) 重写为：(21)Φic,Φbc,ΦpdeTC=Ψ，在这种情况下，线性组合系数 C 直接使用线性最小二乘法确定，如果系数矩阵是秩亏的，则采用最小范数解：(22)C=[Φic,Φbc,ΦpdeT]?Ψ，这里 ? 表示相应矩阵的摩尔-彭若斯广义逆。对于非线性偏微分方程，代数系统方程 (20) 关于 C 是非线性的。在这种情况下，通常使用非线性最小二乘（NLSQ）方法来求解 C。NLSQ 算法通常表现良好，并表现出平滑的收敛行为。

2.3. 带有时间推进和随机特征的物理信息极端学习机（Physics-Informed Extreme Learning Machine with Time-Marching and Random Features）由于初始条件和边界条件的变化，瞬态渗透问题可能表现出尖锐的空间梯度和局部变化。对于由 Gardner 水力模型描述的土壤，控制 Richards 方程可以转换为线性形式；然而，瞬态渗透过程仍然可能涉及压力头分布的显著梯度。此外，瞬态状态和准稳态之间的交替主导进一步复杂化了流动动力学，使得传统的全局方法难以准确捕捉解的细尺度变化。因此，在时间演化过程中有效地模拟多个尺度之间的相互作用同时保持数值稳定性和准确性仍然是一个重大挑战。为了解决这个问题，本研究引入了一种结合了时间推进策略（Time-Marching Strategy, TMS）和随机特征方法（Random Feature Method, RFM）的 PIELM，以高效模拟瞬态非饱和流动（Chen 等人，2022a, 2024, Chen 和 Luo, 2023, Dong 和 Li, 2021）。

2.3.1. 时间推进策略（Time-Marching Strategy）通过时间推进策略，将整个时间域划分为一系列较小的子区间，在每个子区间内，独立的 ELM 构建一个局部近似。每个局部解在 PIELM 框架内定义在一组全局一致的随机特征基函数上。通过顺序连接这些局部解，该方法可以在整个时间域内重构全局解。这种方法的核心思想是将时间域分解，以便分阶段解决复杂的瞬态行为，同时提高计算效率并保持高精度。首先定义时间步长 Δt，并将全局域 Ω 划分为 Nt 个不重叠的子域：(23)Ω=Ω1∪Ω2∪?∪ΩNt，其中 Ωi 表示第 i 个子域。如果 Ωi 和 Ωk (1≤i,k≤Nt) 相邻并共享一个公共边界，则这个共享界面定义为 Γik。需要注意的是，必须使用上一个时间步长网络在共享边界 Γik 上的预测作为当前网络的初始条件，因为最初规定的条件仅在 t=0 时有效。时间推进策略的详细实现如图 3 所示。为了增强相邻时间子域之间的物理耦合，可以结合相应的 PDE 残差。因此，根据方程 (5)，首先构建以下中间函数：(24)fψ^?,z=Kα?2ψ^??z2+Kθ?ψ^??z，其中 ψ^? 表示模型输出。然后，可以如下构建时间域插值项：(25)k1=fψ^i?,z;k2=fψ^i?+dt2k1,z;k3=fψ^i?+dt2k2,z;k4=fψ^i?+dt?k3,z。可以考虑几种方法来建立相邻时间子域之间的耦合。基于 Euler 方案的 PDE 残差表示为：(26)Rpde=ψ^i?+dt?k1-ψ^i+1?，基于 Crank–Nicolson 方案的 PDE 残差表示为：(27)Rpde=ψ^i?+dt2?fψ^i?,z+fψ^i+1?,z-ψ^i+1?，基于四阶 Runge-Kutta 方案的 PDE 残差表示为：(28)Rpde=ψ^i?+dt?·k1+2k2+2k3+k?-ψ^i+1?。在每个时间步长，神经网络在当前时间区间 [t,+Δt] 内构建一个局部近似函数。这个函数同时满足控制方程残差以及子域内的相应初始和边界条件。因此，网络训练实际上是在局部子域上进行的，在该子域内通过最小化 PDE 残差和约束误差来获得时间区间内的解分布。随后，使用当前区间结束时的预测解 t+Δt 作为下一个时间段的初始条件，从而逐步推进整个时间域。这种策略本质上是将时间域分解与顺序时间推进解框架相结合。尽管本研究提出的时间步进方案预测了 t+Δt 时的解，但其原始动机是解决尖锐梯度现象。为了提高适应性，第 4 节进一步研究了一种可变时间推进策略，在其中只有在出现尖锐梯度的区域才选择性地应用时间步进。一般来说，时间步长的选择直接关系到模拟精度和计算效率。过小的时间步长可能导致计算负担过重，导致模拟时间过长，不切实际。相反，过大的时间步长可能会导致忽略渗流过程的重要动态特征，从而影响模拟结果的可靠性和准确性。

2.3.2. 随机特征方法（Random Feature Method, RFM）随机特征方法（RFM）是指 PIELM 在模型中使用全局定义的随机函数作为基函数来近似 Richards 方程的解（Chen 等人，2022b）。具体来说，第 2.2.2 节中的 Φ 表示多变量随机特征函数的集合，形成随机特征矩阵。因此，基于随机特征方法的通用逼近定理可以表示为（Wang & Dong, 2024）：设 Id 表示 d 维空间 [0,1]d，C(Id) 表示在 Id 上定义的连续函数集。对于任何函数 u∈C(Id) 和任何 ε>0，存在一个属于 C(Id) 的极端 ELM 网络，它具有 m 个随机特征基函数，其权重参数 wi 和偏置参数 bi 根据连续概率分布随机初始化并固定，使得：(29)‖∑n=0Ntfn(x)∑j=1majg(wjx+bj)-u(x)‖<ε，x∈id，其中 ‖?‖ 表示 l2 范数，g 是一个绝对可积的激活函数，g(wjx+bj) 表示用于在每个 pielm 中近似 pde 解的随机特征基函数，fn(x) 表示每个子域中的指示函数，如本研究中是时间推进的时间子域。方程 (29) 表明，理论上，这种方法具有解决高维非饱和流动问题的潜力。网络激活函数和随机初始化方法的选择在很大程度上决定了随机函数家族的特性。确保随机特征函数空间足够丰富以准确近似 pde 的解是至关重要的。图 4 展示了不同家族的随机特征基函数的生成。

3. 结果（results）在本节中，进行了六项数值实验来评估和验证所提出方法在求解 richards 方程方面的有效性。所提出的方法，此后称为 tr-pielm，应用于一维稳态、一维瞬态和二维瞬态非饱和地下水流动问题。将获得的结果与解析解、使用四阶 runge–kutta 方案（rk4）的传统有限差分方法（fdm）和物理信息神经网络（pinn）进行比较，以验证所提出方法的准确性和效率。具体来说，为了进一步研究复杂的地质技术场景，还检验了三个额外的模型：一个双层非饱和土壤地下水渗透模型、一个一维非饱和渗透和土壤变形的耦合模型，以及一个考虑 van genuchten-mualem 模型的一维瞬态非饱和渗透模型用于边坡稳定性分析。所提出方法的结果与基于 rk4 的数值解进行了比较，表明 tr-pielm 达到了可比的计算精度。在数值实验中，tr-pielm 的架构定义为：(30)tr-pielmarch=[nin,n=200,nout]，其中 nin 和 nout 分别表示输入节点和输出节点的数量，在本研究中 nout=1。作为参考，pinn 采用更深的网络架构，定义为：(31)pinnarch=[nin,40,40,40,40,40,nout]，这种配置确保了 pinn 和 pielm 中的神经元总数保持一致。pinn 在每次测试中采用与 tr-pielm 相同的初始化和激活函数。 ‖?‖ 表示 l2 范数，g 是一个绝对可积的激活函数，g(wjx+bj) 表示用于在每个 pielm 中近似 pde 解的随机特征基函数，fn(x) 表示每个子域中的指示函数，如本研究中是时间推进的时间子域。方程 (29) 表明，理论上，这种方法具有解决高维非饱和流动问题的潜力。网络激活函数和随机初始化方法的选择在很大程度上决定了随机函数家族的特性。确保随机特征函数空间足够丰富以准确近似 pde 的解是至关重要的。图 4 展示了不同家族的随机特征基函数的生成。 3. 结果（results）在本节中，进行了六项数值实验来评估和验证所提出方法在求解 richards 方程方面的有效性。所提出的方法，此后称为 tr-pielm，应用于一维稳态、一维瞬态和二维瞬态非饱和地下水流动问题。将获得的结果与解析解、使用四阶 runge–kutta 方案（rk4）的传统有限差分方法（fdm）和物理信息神经网络（pinn）进行比较，以验证所提出方法的准确性和效率。具体来说，为了进一步研究复杂的地质技术场景，还检验了三个额外的模型：一个双层非饱和土壤地下水渗透模型、一个一维非饱和渗透和土壤变形的耦合模型，以及一个考虑 van genuchten-mualem 模型的一维瞬态非饱和渗透模型用于边坡稳定性分析。所提出方法的结果与基于 rk4 的数值解进行了比较，表明 tr-pielm 达到了可比的计算精度。在数值实验中，tr-pielm 的架构定义为：(30)tr-pielmarch=[Nin,n=200,Nout]，其中 nin 和 nout 分别表示输入节点和输出节点的数量，在本研究中 nout=1。作为参考，PINN 采用更深的网络架构，定义为：(31)pinnarch=[Nin,40,40,40,40,40,Nout]，这种配置确保了 pinn 和 pielm 中的神经元总数保持一致。pinn 在每次测试中采用与 tr-pielm>
3. 结果（results）在本节中，进行了六项数值实验来评估和验证所提出方法在求解 richards 方程方面的有效性。所提出的方法，此后称为 tr-pielm，应用于一维稳态、一维瞬态和二维瞬态非饱和地下水流动问题。将获得的结果与解析解、使用四阶 runge–kutta 方案（rk4）的传统有限差分方法（fdm）和物理信息神经网络（pinn）进行比较，以验证所提出方法的准确性和效率。具体来说，为了进一步研究复杂的地质技术场景，还检验了三个额外的模型：一个双层非饱和土壤地下水渗透模型、一个一维非饱和渗透和土壤变形的耦合模型，以及一个考虑 van genuchten-mualem 模型的一维瞬态非饱和渗透模型用于边坡稳定性分析。所提出方法的结果与基于 rk4 的数值解进行了比较，表明 tr-pielm 达到了可比的计算精度。在数值实验中，tr-pielm 的架构定义为：(30)tr-pielmarch=[nin,n=200,nout]，其中 nin 和 nout 分别表示输入节点和输出节点的数量，在本研究中 nout=1。作为参考，pinn 采用更深的网络架构，定义为：(31)pinnarch=[nin,40,40,40,40,40,nout]，这种配置确保了 pinn 和 pielm 中的神经元总数保持一致。pinn 在每次测试中采用与 tr-pielm 相同的初始化和激活函数。>初始残差和偏微分方程（PDE）残差的惩罚系数统一设置为1。同时，在每种情况下，PINN都采用与TR-PIELM相同的边界条件处理方法，例如硬编码，并采用了一种改进的MLP架构，结合了随机傅里叶特征嵌入（Wang等人，2023a）。对于一维稳态问题，PINN使用100个随机采样的配置点。对于一维瞬态问题，使用拉丁超立方抽样（LHS）生成10,000个配置点。对于二维情况，使用100,000个LHS配置点。由于初始条件对瞬态问题至关重要，因此包括了初始条件边界上的所有点。PINN解在256个空间网格点和100个时间网格点上进行评估。训练过程总共进行了10,000次Adam优化器的迭代以实现全局收敛，随后使用L-BFGS算法进行微调。为了定量评估所提出方法的准确性，进行了详细的误差分析。引入了相对L2误差作为评估指标，以衡量预测值和参考值之间的差异，较小的值表示更高的预测准确性，其定义为：
(32) RelativeL2Error = ∑i=1N‖ψ^-ψ‖L2 / ∑i=1N‖ψ‖L2
其中ψ^表示从框架获得的预测值，而ψ表示真实值或参考值。‖?‖L2表示L2范数，N是整个计算域内的时空点总数。在本研究中，所有数值实验都是在PyTorch中实现的，并在配备NVIDIA Tesla T4 GPU和第12代Intel(R) Core(TM) i7-12700H（2.30 GHz）32.0 GB CPU的计算机上执行的。

3.1 一维稳态非饱和流动
该测试模拟了一维稳态地下水在非饱和均匀土壤中的流动，如图5所示。土壤属性由参数α=8.0×10^-3m^-1、θs=0.35、θr=0.14和Ks=1×10^-4m/h表征（Liu等人，2015年；Zhu等人，2022a年）。土壤层厚度为L=10m。根据第2.1节，一维稳态流动的控制方程如下：
(33) ?2ψ*?z2 + α?ψ*?z = 0
该问题的边界条件定义为：
(34) ψz=0 对于z=0；ψz=L 对于z=L
其中ψd被视为-1000 m。
下载：下载高分辨率图像（48KB）
下载：下载全尺寸图像
图5. 均匀非饱和土壤的一维稳态渗透问题。
该问题的解析解（Tracy，2006年）表示为：
(36) ψz = 1/αln(1-e^αψd) / (1-e^α)^2L - zsin(αz/2) + e^αψd
为了数值计算的方便，控制方程经过空间-时间归一化和无量纲化处理，其中z′=z/L。关于边界条件的施加，可以应用第2.2.2节中描述的软约束方法和硬约束方法。在这个测试中，构建了一个试验解来将边界条件作为硬约束施加：
(37) ψ^z′ = gz′(1-z′) + z′(1-z′)ΦC
其中gz′=ψd用于施加边界条件，ΦC是PIELM的最终输出。一维稳态问题可以被视为具有单个时间不变子域的物理系统。因此，只需要一个PIELM替代模型进行模拟。在这个测试中，PIELM使用Swish激活函数和统一的权重及偏置初始化。在计算域内分布了总共100个配置点。
图6展示了所提出方法获得的模拟水头及其相应的误差，显示出与解析解几乎相同的准确性。
下载：下载高分辨率图像（198KB）
下载：下载全尺寸图像
图6. (a) 均匀多孔非饱和介质中一维稳态地下水流动的PIELM解与解析解的比较；(b) PIELM解与解析解的绝对误差。
表1展示了不同方法计算水头得到的准确性结果。TR-PIELM的结果为1.250e-16±3.393e-18，这是基于10次独立运行（使用不同的随机种子）的相对L2误差的平均值±标准差，最小相对L2误差为1.195e-16。相对较小的标准差表明所提出的方法对随机初始化具有鲁棒性。与现有文献中的结果（Zhu等人，2022a年；Zhu等人，2019年）相比，所提出的方法在模拟一维稳态流动方面表现出压倒性的优势。值得注意的是，PIELM对于某些线性或轻微非线性的PDE可以实现极高的准确性，因为物理信息损失、固定的随机特征和最小二乘输出训练可以以非常小的优化和离散化误差近似平滑的解流形。然而，这并不意味着PIELM可以普遍替代数值求解器。其有效性取决于问题的复杂性、维度、条件以及物理编码的质量。因此，PIELM应被视为一种补充的高效替代模型或求解器加速器，而不是经典数值算法的通用替代品。

表1. 不同方法计算水头的相对L2误差。
方法 计算时间（秒） 相对L2误差
PI (Δz=0.025) 23.26 2.267e-2
GS (Δz=0.025) 18.43 2.422e-4
SOR (Δz=0.025) 3.34 5.682e-6
MP(4)-GS (Δz=0.025) 0.27 4.651e-7
FDM (Δz=0.025) 0.097 3.409e-9
TS-IJIM (Δz=0.02) 0.41 9.449e-7
FDM (Δz=0.02) 0.17 9.232e-9
PINN 29 1.64 3.713e-6
PIELM 0.008 4.912e-16

3.2 一维瞬态非饱和流动
本节研究了一维瞬态非饱和流动问题，以进一步评估所提出的方法。对于均匀多孔介质中的地下水流动，可以通过求解控制方程获得瞬态解。根据第2.1节，本研究中的控制方程可以表示为：
(38) ?2ψ*?z2 + α?ψ*?z = c?ψ*?t
其中c=αθs-θrKs。土壤层厚度和总模拟时间分别设为L=10 m和T=5 h。该问题的初始和边界条件如下：
(39) ψz,t=0 对于z=0；ψz=L,t=0
其中ψd被视为-10 m。该问题的解析解（Tracy，2006年）表示为：
(42) ψz,t = 1/αln(ψt*z,t + ψs*z + e^αψd)
其中：
(43) ψs*z = (1-e^αψd) × (1-e^(-αz/2)sin(λmz) / e^(-μmz)
其中λm=kπL，μm=α^2/4+λm^2c。
本研究为了测试目的，研究了两种类型的土壤，即粉质壤土和沙质土壤。每种土壤类型的相应参数列在表2中（Ku等人，2018年；Zhu等人，2022a年）。
表2. 沙质土壤和粉质壤土的参数。
土壤类型 α(m^-1) θs θr Ks(m/h)
粉质壤土 8.0×10^-3 0.46 0.14 1.5×10^-3
沙质土壤 1.6×10^-2 0.50 0.11 6.0×10^-2
为了计算方便，在模拟之前，控制方程根据z′=z/L和t′=t/T进行了时间和空间归一化和无量纲化。采用类似于第3.1节中的方程(37)的试验解来硬编码边界条件。在本研究中，TR-PIELM采用了Cos激活函数和统一初始化。时间推进方案使用了四阶龙格-库塔迭代，每个子域有400个等分的时间子域和100个配置点。
图7显示TR-PIELM的结果与非饱和土壤中一维瞬态地下水流动的解析解一致。结果表明，所提出的方法能够有效模拟非饱和均匀多孔介质中的一维瞬态地下水流动。
下载：下载高分辨率图像（406KB）
下载：下载全尺寸图像
图7. TR-PIELM解及其与均匀多孔非饱和介质中一维瞬态地下水流动的解析解的比较。(a) 粉质壤土；(b) 沙质土壤。
此外，为了评估所提出方法的准确性，进行了误差分析，并与其他解决方法进行了比较，结果总结在表3中。对于粉质壤土，TR-PIELM在10次独立运行中的相对L2误差为3.338e-4±1.682e-5，最小误差为3.014e-4。对于沙质土壤，相应的误差为4.296e-4±9.806e-5，最小误差为3.654e-4。研究发现，即使在相对较少的子域分解下，所提出的方法也能达到与高分辨率网格划分的RK4方案相当的准确性。因此，所提出的方法在相同的准确性下显著提高了计算效率，而随机特征函数的近似在一定程度上缓解了传统数值方法固有的收敛限制。在本研究中，PINN方法在计算时间和准确性方面都没有显著优势。

表3. 不同方法计算水头的相对L2误差。
土壤类型 方法 计算时间（秒） 相对L2误差
粉质壤土 FDM (Δz=0.01) 17.74 1.921e-3
RK4 (Δz=0.01) 10.58 7.025e-4
PINN 87 2.445.075e-2
TR-PIELM 0.70 3.338e-4
沙质土壤 FDM (Δz=0.01) 17.36 2.491e-3
RK4 (Δz=0.01) 10.55 7.959e-4
PINN 85 9.54 1.737e-2
TR-PIELM 0.72 4.296e-4

3.3 两层土壤中的一维非饱和渗透
在这个案例研究中，所提出的方法用于模拟通过两层土壤系统的一维非饱和地下水流动，如图8所示。该模型考虑了具有不同物理性质的分层土壤。土壤层1和土壤层2的厚度分别为L1=5 m和L2=5 m，它们的物理参数列在表4中。此外，本研究考虑了两种不同的两层土壤系统，以评估所提出方法对具有不同土壤层的异质土壤系统的适用性。
表5. 两种不同两层土壤系统的模型参数。
表4. 两种不同两层土壤系统的参数。
土壤类型 α(m^-1) θs θr Ks(m/h)
土壤1 第一层 8.0×10^-3 0.46 0.14 1.5×10^-3
土壤1 第二层 1.6×10^-2 0.50 0.11 6.0×10^-2
土壤2 第一层 3.2×10^-4 0.46 0.14 1.0×10^-6
土壤2 第二层 8.0×10^-4 0.50 0.11 1.0×10^-2
表5. 不同方法计算水头的相对L2误差。
土壤类型 方法 计算时间（秒） 相对L2误差
未耦合 FDM (Δz=0.01) 12.94 --
RK4 (Δz=0.01) 29.83 --
PINN 88 1.61 1.311e-2
TR-PIELM 0.83 1.862e-4
耦合 F>0 FDM (Δz=0.01) 13.27 --
RK4 (Δz=0.01) 30.27 --
PINN 89 0.74 2.721e-2
TR-PIELM 0.78 2.486e-4
耦合 F<0 FDM (Δz=0.01) 12.96 --
RK4 (Δz=0.01) 30.06 --
PINN 89 1.32 1.259e-2
TR-PIELM 0.81 1.714e-4
控制方程、初始条件、边界条件和总时间与第3.2节中定义的一致。两层土壤系统中的界面处理特别关键，因为它必须满足物理守恒和连续性条件。根据第2.1节中的线性化控制方程，在使用有限差分方法生成参考解时，两层土壤之间的界面处的有限差分方程修改如下：
(45) Kα1ψi-1?j-2ψi*j + 2α2ψi+1?jΔz2 + Kθ1ψi*j-ψi-1?j + β2ψi+1?j-ψi*jΔz = ψi*j-ψi*j-1Δ
其中，i、j、Δz和Δt分别表示空间节点索引、时间节点索引、沿z轴的空间离散间隔和时间间隔。α1=Kα1/(Kα1+Kα2)，α2=Kα2/(Kα1+Kα2)，β1=Kθ1/Kθ2和β2=Kθ2/Kθ1是一组新引入的参数，用于处理两层土壤之间的界面。根据上述推导，两层土壤之间的土壤属性可以表示为：
(46) Kαi = 2Kα1/Kα2
(47) Kθi = 2Kθ1/Kθ2
(48) Kθi = 2Kθ1/Kθ2
在这个测试中，每层土壤都由单独的TR-PIELM替代模型模拟。两个TR-PIELM通过上述界面连续性条件耦合在一起，以模拟两层土壤系统，这在概念上类似于空间域分解。每个TR-PIELM使用Tanh激活函数，时间推进方案使用四阶龙格-库塔迭代，每个子域有100个配置点。每个TR-PIELM的控制方程都经过空间-时间归一化和无量纲化，其中z1′=z1/L1，z2′=z2/L2，t′=tT。如第2.2.2节所讨论的，也可以应用边界软约束以获得满意的结果。在本研究中，层1的下边界和层2的上边界通过以下试验解进行硬约束：
(48) ψ^1z1′=ψ^z1′=0 + z1′-0 / 1-0Φ1C1
(49) ψ^2z2′=ψ^z2′=2 + 2-z2′ / 2-1Φ2C2
其中1和2分别表示代表土壤层1和层2的两个网络。模拟结果如图9所示。
下载：下载高分辨率图像（374KB）
下载：下载全尺寸图像
图9. 两层多孔非饱和土壤中一维瞬态地下水流动的TR-PIELM解。实线代表RK4解，虚线代表TR-PIELM预测解。(a) 土壤1；(b) 土壤2。
TR-PIELM的模拟结果进一步与RK4参考解进行了比较，如下所示。
图10显示了每种方法的计算时间。由于这种情况没有提供解析解，图表报告了PINN和所提出方法相对于RK4参考解的相对L2误差，Δz=0.01。TR-PIELM获得的解与RK4参考解吻合良好，计算效率高。对于土壤1，TR-PIELM在10次独立运行中的相对L2误差为2.665e-3±1.473e-4，最小误差为2.406e-3。对于土壤2，相应的误差为6.220e-3±5.320e-4，最小误差为5.488e-3。此外，RK4方法在模拟两层土壤渗流问题时需要极其细的网格，时间步数为3990000，导致计算成本很高。相比之下，所提出的方法在较少的时间子域下实现了相当的准确性，显示出更高的计算效率。
下载：下载高分辨率图像（154KB）
下载：下载全尺寸图像
图10. 不同方法计算水头的计算时间和相对L2误差。

3.4 一维耦合非饱和渗透和变形
为了进一步评估所提出方法的有效性，本节考虑了一个地质工程中的经典物理问题：非饱和土壤中的一维耦合渗透和变形过程。将所提出方法的结果和计算效率与RK4方法进行了比较。Wu等人（此外，本研究考虑了以下初始条件：底部采用狄利克雷边界条件以保持恒定的水头压力，顶部采用罗宾边界条件以保持恒定的流量通量。(52)uz,t=0=e-αz(53)uz=0,t=1(54)?uz=L,t?z+uz=L,t=qKs，其中q=3.6×10-3m/h代表降雨通量。上述方程中的系数α=0.1m-1，M=θsα+γω/F，θs=0.3表示饱和含水量，γω=10kN/m3表示水的单位重量，F表示部分饱和土壤在变化土壤吸力下的弹性模量，这是应力的函数。当我们将F设置为-1000kPa时，土壤会发生坍塌；当我们将F设置为1000kPa时，土壤会发生膨胀。当不考虑耦合关系时，M=θsα。Ks=3.6×10-3m/h是饱和水力传导系数（Wu等人，2016b）。土壤厚度为L=5m，总演化时间设为T=10h。与前面的章节一样，在计算之前，控制方程在空间和时间上都被归一化和无量纲化。在TR-PIELM框架中，采用了Tanh激活函数和均匀初始化。时间推进方案采用了包含100个配置点的400个时间子域的欧拉方案。边界条件通过软约束施加，而狄利克雷边界条件可以作为硬约束强制执行。为了说明目的，图11展示了三种情况的结果：未耦合情况、F>0的耦合情况和F<0的耦合情况。下载：下载高分辨率图像（704KB）下载：下载全尺寸图像图11. 在未耦合和耦合条件下水头压力的结果及其与RK4的比较。(a) 未耦合；(b) F>0的耦合；(c) F<0的耦合。结果表明，所提出的方法具有与RK4模拟相当的高预测精度。此外，随着渗透时间的增加，坍塌土壤情况（F<0）的水头压力相对较高，而膨胀土壤情况（F>0）的水头压力相对于未耦合条件较低。上述比较总结了PINN方法和所提出方法与RK4方案在计算性能方面的优势。可以观察到，所提出的方法在准确性和计算效率方面都表现出相当的竞争优势。对于未耦合土壤情况，TR-PIELM在10次独立运行中的相对L2误差为1.862e-4±2.767e-5，最小值为1.450e-4。对于膨胀土壤情况（F>0），相应的误差为2.486e-4±5.948e-5，最小值为1.803e-4。对于坍塌土壤情况（F<0），误差为1.714e-4±3.251e-5，最小值为1.190e-4。此外，RK4方法受到离散化方案稳定性要求的限制，通常需要精细的空间和时间网格，这显著降低了其计算效率。相比之下，所提出的方法有效地缓解了这一限制，使其特别适合快速预测土壤水分分布。3.5. 1D瞬态非饱和流用于边坡稳定性分析本节研究了与边坡稳定性分析相关的1D瞬态非饱和流问题，以进一步评估所提出方法的性能。对于具有均匀多孔介质的边坡中的地下水流动，可以通过求解控制方程来获得瞬态响应。根据第2.1节，本研究中考虑的控制方程如下（Wu等人，2020）：(55)??zKzψ?ψ?z+cosβ?Kzψ?z=?θ?t，其中β=34°。本问题中的初始和边界条件、土壤层厚度、总模拟时间和土壤物理性质与第3.2节中的设置相同。然而，为了进一步评估所提出方法在更复杂条件下的适用性，本节考虑了van Genuchten-Mualem模型来表示土壤本构关系θψ和Kψ（Wu, Liu等人，2016）：(56)θψ=θr+θs-θr[1+(αψn)]m，m=1-1n和：(57)Kψ=KsSe12[1-(1-Se1m)m]2，Se=θ-θrθs-θr，其中n=1.41适用于粉质壤土，n=1.89适用于沙质土壤。Se表示有效饱和度。安全系数Fs定义如下，用于分析无限长边坡的稳定性（Wu等人，2020）：(58)Fs=c′Lγtcosβsinβ+tanφ′tanβ-ψγw[(θψ-θr)/(θs-θr)]tanφ′Lγtcosβsinβ，其中c′和φ′表示有效粘聚力和摩擦角，γt和γw分别表示土壤和水的单位重量。这些参数在表6中总结。表6. 粉质壤土和沙质土壤的边坡安全系数参数土壤质地c′(kPa)φ′γt(kN/m3)γw(kN/m3)粉质壤土1020°19.09.8沙质土壤032°19.09.8表7. 不同方法在水头压力模拟中的相对L2误差。方法计算时间（秒）相对L2误差FDM (Δz=0.01)88.765.860e-3RK4 (Δz=0.01)321.401.814e-3PINN27761.485e-2TR-PIELM28.151.576e-3在本研究中，边界条件仍然通过硬编码强制执行。TR-PIELM采用了Tanh激活函数和随机正态初始化。采用了欧拉时间推进方案，时间域被均匀划分为400个子域，每个子域包含100个空间配置点。由于van Genuchten-Mualem模型引入的强非线性，这种情况比其他测试案例更具挑战性。为了减少计算负担，并基于第4节中的敏感性分析，本例中使用了相对较少的神经元数量60。此外，在TR-PIELM的每个时间子域内还使用了通常用于数值方法的Picard迭代来稳定计算。图12展示了TR-PIELM预测的边坡内非饱和土壤中一维瞬态地下水流动的安全系数Fs。此情况的参考解是使用有限差分方法生成的，并对结果进行了比较。粉质壤土和沙质土壤的预测Fs值与参考值之间的相对L2误差分别为4.730e-4±4.952e-5和4.687e-4±4.501e-5，最小误差分别为4.069e-4和3.881e-4，相应的TR-PIELM计算时间分别为3.597秒和1.458秒。这些结果表明，所提出的方法可以准确估计边坡安全系数，并可以作为有效的替代模型。相比之下，PINN方法在计算效率或准确性方面没有明显优势。其相对于参考解的相对L2误差分别为5.448e-2和2.969e-2，计算时间分别为1318.182秒和1310.385秒（经过10,000个时期后）。因此，与PINNs相比，所提出的方法显著提高了计算效率。此外，与传统的数值方法相比，随机特征近似为替代建模提供了更大的灵活性。下载：下载高分辨率图像（402KB）下载：下载全尺寸图像图12. TR-PIELM对具有均匀多孔非饱和介质的边坡中的一维瞬态地下水流动的安全系数Fs的解。(a)(b) 粉质壤土；(c)(d) 沙质土壤。3.6. 2D瞬态非饱和流在前面的章节中，所提出的方法被用来模拟1D非饱和地下水流动，并通过与解析解或参考解的比较进行了验证。在本研究中，进一步调查了所提出方法在2D瞬态非饱和流中的适用性。考虑了一个深度为Z=1m、宽度为W=1m的均匀2D土壤域，以及参数α=9.0×10-3m-1、θs=0.46、θr=0.14和Ks=8.0×10-4m/h。总演化时间设为T=5小时。二维非饱和流动的控制方程表示如下：(59)??xKxψ?ψ?x+??zKzψ?ψ?z+?Kzψ?z=?θ?t根据第2.1节，控制方程可以重新表述为：(60)?2ψ??x2+?2ψ??z2+α?ψ??z=c?ψ??t本研究的初始条件设置为ψx,z,t=0=ψd=-1m。根据Green-Ampt边界条件，在土壤表面施加渗透，产生相应的水头（Heber Green & Ampt，1911）。在土壤域的中点，水头设为零，并逐渐向两侧边界减小。这些条件作为边界约束施加在土壤表面。此外，土壤的左侧、右侧和底部边界保持在干燥条件下。边界条件总结如下：(61)ψx=0,z,t=ψd(62)ψW,z,t=ψd(63)ψx,z=0,t=ψd(64)ψx,z=L,t=1αlneαψd+1-eαψdsin(πxW)这个问题的解析解（Tracy，2006）可以通过以下方程推导出来：(65)ψx,z,t=1αln(ψt?x,z,t+ψs?x,z,t+eαψd)(66)ψt?x,z,t=2(1-eαψd)Lcsin(πxW)eα(L-z)2∑k=1∞(-1)k(λkμk)sin(λkz)e-μkt(67)ψs?x,z,t=(1-eαψd)eα(L-z)2sin(πxW)sinh(β1z)sinh(β1L)(68)β1=α24+(πW)2在模拟之前，控制方程像前面的案例一样进行了时空归一化。在本研究中，TR-PIELM采用了Tanh激活函数和均匀初始化，400个等分的时间子域和100个配置点。所有边界条件都是狄利克雷类型，可以使用以下试验解构造在替代模型中硬编码。(69)ψ^z′=gz′+x′(1-x′)z′(1-z′)ΦC，其中gz′=1-eαψdsin(πxW)。图13展示了使用所提出方法和参考解获得的水头压力结果之间的比较。下载：下载高分辨率图像（571KB）下载：下载全尺寸图像图13. TR-PIELM预测与参考解之间的2D瞬态非饱和流动水头ψ的比较。接下来，将所提出方法预测的水头压力结果与RK4方法的结果进行了比较，如下所示。TR-PIELM的相对L2误差为1.576e-3±3.950e-04，这是基于10次独立运行（使用不同的随机种子）的平均值±标准差；观察到的最小误差为1.164e-3。这些结果表明，所提出方法得到的解决方案在不同时间点与解析解吻合得很好。此外，所提出的方法消除了PINN模型中对权重调整的依赖。边界条件的硬编码使替代模型能够更好地捕捉物理过程的演变，实现了与传统方法相当的精度，同时提高了计算效率，从而能够更高效和可靠地表示土壤中水流的动态行为。4. 敏感性分析在前面章节中，所提出的TR-PIELM方法被用来解决非饱和土壤中的某些地下水渗透问题，验证了其可行性和有效性。然而，作为一种基于浅层神经网络的方法，TR-PIELM需要反复实验来确定合适的超参数，包括激活函数类型、初始化方案、神经元数量、配置点数量、软边界约束的损失权重和时间推进策略。其中一些因素直接或间接影响随机特征基函数的功能族，而其他因素则影响计算过程中的误差累积。为了研究这些超参数如何影响所提出方法的计算性能，并为未来的研究提供见解和指导，本节以第3.2节中的一维瞬态非饱和流问题的粉质壤土案例为例。使用第3.2节中呈现的解析解作为参考。需要注意的是，在每次分析中，只改变一个影响因素，而其他所有因素保持不变，以确保结果的可靠性和可比性。以下讨论总结并分析了这些因素对结果的影响。4.1. 激活函数和初始化本节研究了不同激活函数和初始化方案对所提出方法计算性能的影响。实际上，这两个因素直接影响用于近似偏微分方程解的随机特征基函数的功能属性。图14展示了TR-PIELM在不同激活函数和初始化策略下获得的预测水头压力的相对L2误差。图中显示的缩写表示测试的配置，其中常用的激活函数包括T（Tanh）、S（Sigmoid）、R（ReLU）、GE（GELU）、Swish（GaussianSwish）、G（Gaussian-RBF）、W（Wavelet）、M（Mish）、L（LeakyReLU）和C（Cos），初始化方法标记为XU（Xavier_Uniform）、XN（Xavier_Normal）、KU（Kaiming_Uniform）、KN（Kaiming_Normal）、LN（Lecun_Normal）、U（Uniform）、N（Normal）、SR（Sparse）和O（Orthogonal）。这些超参数可以从现有文献中参考，这里不再详细说明，因为它们超出了本研究的主要范围。下载：下载高分辨率图像（163KB）下载：下载全尺寸图像图14. 粉质壤土中1D瞬态非饱和流动的水头压力模拟的相对L2误差，分别在(a)不同激活函数；(b)不同初始化方法下。结果表明，ReLU和LeakyReLU激活函数不适合本研究中研究的问题，而Cos激活函数被证明是TR-PIELM模拟一维瞬态非饱和流动的最佳选择。这归因于非饱和流的指数时间演变和非线性特性，例如在t=0时初始条件和边界条件之间的冲突导致的解决方案中的尖锐梯度。此外，如第2.3节图4所示，在所提出的单层神经网络中，ReLU激活函数无法为随机特征基函数赋予平滑的非线性特征。相比之下，具有傅里叶特性的激活函数在近似尖锐梯度方面表现出更好的效果。均匀初始化方法取得了最佳性能，其次是Normal和Kaiming_Uniform方案。在均匀初始化中，网络权重和偏置在[-R,R]的均匀分布区间内随机初始化，确保参数的足够随机性，并防止所有神经元从过于相似的状态开始。进一步研究了不同的R值，相应的结果显示在图15中。下载：下载高分辨率图像（87KB）下载：下载全尺寸图像图15。在均匀初始化条件下，使用不同R值对粉壤中一维瞬态非饱和流的压力头模拟进行相对L2误差分析。研究表明，R值对结果有显著影响。随着R值的增加，相对L2误差先减小后增大，表明存在一个近似最小点。对于当前问题，R值在1.0左右时性能最佳。实际上，最优R值会因不同问题而异。关于初始化问题，Peng等人（Peng et al., 2025）报告称，PIELM在求解偏微分方程（PDEs）时确实对初始化敏感。为了解决这个问题，他们引入了基于秩的神经网络（RINN）框架来增强网络输出基函数的正交性，有效降低了初始化敏感性，并确保PIELM在各种初始化设置下保持高精度。

4.2 神经元数量和配置点数量
对于TR-PIELM来说，神经元数量和配置点数量对模型性能有显著影响。神经元数量决定了随机特征基函数的数量，而配置点数量则影响模型的表示能力。为了评估所提方法对这些因素的敏感性，进行了一系列研究，相应的结果展示在图16中。
下载：下载高分辨率图像（215KB）
下载：下载全尺寸图像
图16. 在不同神经元数量（a）和不同配置点数量（b）下，粉壤中一维瞬态非饱和流的压力头模拟的相对L2误差。
结果表明，当神经元数量过少时，模型的表示能力受限，导致精度较低。一旦神经元数量达到20个或更多，模型就能达到稳定的精度，增加神经元数量不再有显著提升。在低维问题中，Chen等人（Chen et al., 2024）通过完整的SVD分析了随机特征矩阵的秩，发现其奇异值迅速衰减，无论神经元数量或配置点数量如何增加，秩都表现出明显的上限。此外，计算时间随神经元数量略有增加，但增长幅度不大，额外增加2000个神经元仅使计算时间增加约1.2秒。同样，当配置点数量少于40个时，模型捕捉非饱和流的空间特性的能力有限，导致精度较低。随着配置点数量从40增加到2000个，相对L2误差在同一个数量级内波动，先略有减小后略有增加。计算时间也相应增加，但保持在大约0.4到1.4秒的狭窄范围内。

4.3 损失权重
作为一种基于物理的机器学习方法，TR-PIELM需要强制执行边界条件（如2.2.2节所述），以唯一确定PDE解。因此，本节研究了不同边界条件执行策略对模型性能的影响。具体来说，采用Euler方案作为时间推进策略，在软边界约束下，边界条件残差权重表示为λ，而PDE残差权重固定为1。测试的λ值包括0.01、0.1、1、5、10、20、50、100和1000，H表示硬边界条件执行情况。
图17展示了在不同边界损失权重λ下，一维瞬态非饱和流压力头预测的相对L2误差。结果表明，当λ过小时，模型无法有效执行边界条件，导致解无效。具有软边界约束的TR-PIELM模型表现相对较差。有趣的是，当λ超过某个水平时，精度会降低，因为过大的边界权重会导致模型忽略内部PDE残差，这种行为与通常在PINNs中观察到的行为类似。对于当前问题，最优的软边界权重为λ=1，而硬边界执行可以获得最高的精度。
下载：下载高分辨率图像（88KB）
下载：下载全尺寸图像
图17. 在使用Euler方案时间推进策略和软约束边界条件（BCs）的情况下，不同BC损失权重下粉壤中一维瞬态非饱和流压力头模拟的相对L2误差。H表示硬编码的BC。

5. 时间推进策略的变异性
除了前述影响随机特征基函数特性和网络结构表达能力的超参数外，时间推进方案的选择也对所提方法的性能起着关键作用。本节重点分析不同时间推进策略对模型性能的影响。如2.3节所述，整个时空域中的时间推进策略包括Euler（E）、Crank–Nicolson（C）和四阶Runge–Kutta（R）方案。这些方案的比较结果展示在图18中。
下载：下载高分辨率图像（116KB）
下载：下载全尺寸图像
图18. 在不同时间推进积分方案下，粉壤中一维瞬态非饱和流压力头模拟的相对L2误差。深蓝色条‘E’：Euler方案；红色条‘C’：Crank–Nicolson方案；浅蓝色条‘R’：四阶Runge-Kutta方案。（关于图中颜色参考的解释，请参阅本文的网络版本。）
Crank–Nicolson和四阶Runge–Kutta方案都比Euler方案具有更高的精度。然而，无论选择哪种时间推进策略，过长的时间域都无法捕捉到非饱和流中压力头的急剧梯度。此外，随着时间子域数量的增加，计算时间呈指数级增长，而整体精度没有显著提高。这种限制是由于与额外子域相关的浮点误差累积限制了可实现的精度。
因此，本节进一步探讨了通过考虑非均匀时间子域或可变时间推进策略来减少计算误差同时控制模拟时间的策略。为了准确捕捉非饱和流中压力头梯度的变化，采用时间推进策略来解决这些陡峭的梯度。一旦陡峭梯度过渡到更平滑的状态，可以使用更大的时间域进行模拟。基于这一想法，提出了一种两阶段的“早期+后期”可变时间推进策略。具体来说，第一阶段使用细粒度的时间子域划分来解决陡峭梯度，而第二阶段将剩余时间段视为一个没有子域划分的单一全局时间域。图19总结了不同阶段划分方案的结果，以使用均匀4k子域时间推进策略的RK4方案作为参考。
下载：下载高分辨率图像（724KB）
下载：下载全尺寸图像
图19. 在可变时间推进策略下，粉壤中一维瞬态非饱和流压力头模拟的相对L2误差。蓝色水平虚线代表参考情况的相对L2误差，橙色水平虚线表示其计算时间。(a) 0.2T + 0.8T；(b) 0.3T + 0.7T；(c) 0.6T + 0.4T；(d) 0.7T + 0.3T；(e) 0.8T + 0.2T；(f) 0.9T + 0.1T。（关于图中颜色参考的解释，请参阅本文的网络版本。）
结果表明，所有可变时间推进策略所需的计算时间都少于参考情况。此外，对于“0.7T + 0.3T”之后的方案，始终存在计算误差低于参考误差的结果，表明可变时间推进策略不仅提高了精度，还减少了计算时间，从而显著提高了效率。特别是，“0.9T + 0.1T”方案实现了最低的L2误差2.628e-04。

6. 结论
在地下水流动的数值模拟中，传统数值方法通常需要高度精细的离散化网格、较长的计算时间，并受到稳定性和收敛条件的限制。为了寻找替代方案，本研究提出了一种结合物理信息的极端学习机（Physics-Informed Extreme Learning Machine）与时间推进随机特征方法（TR-PIELM），旨在高效模拟非饱和地下水流动过程。所提方法采用浅层神经网络框架，将非饱和流动的控制方程和初始边界条件嵌入损失函数中。通过结合软约束和边界硬编码，模型有效增强了其表示物理约束的能力。构建了六个基准案例来系统评估所提方法：（1）一维稳态渗透，（2）一维瞬态渗透，（3）两层多孔介质中的一维非饱和流动，（4）与土壤变形耦合的一维渗透，（5）用于坡度稳定性分析的一维瞬态非饱和流动，以及（6）二维瞬态渗透。与解析解和有限差分解的比较表明，所提方法在稳态和瞬态条件下都能准确再现压力头分布。与传统PINNs相比，它具有更高的精度和计算效率，并且预测精度可与四阶Runge–Kutta（RK4）方案相媲美。值得注意的是，对于稳态问题，所提方法实现了接近解析的精度，并显著提高了效率。此外，本研究还探讨了影响模型性能的关键因素，包括激活函数、初始化方案、随机特征函数的数量和权重、配置点数量、软边界约束的损失权重以及时间推进策略，并为模型改进提供了实用指导。结果证实，所提出的TR-PIELM为快速准确模拟非饱和流动问题提供了一个稳健且高效的框架。
Vanilla ELM确实存在速度-精度权衡，其固定的随机隐藏特征和封闭形式的输出使得训练非常快速，但可能会降低对高度复杂映射的表示灵活性。因此，ELM最适合需要快速近似、重复评估或替代建模的问题。尽管Richards方程本身并不需要极端效率，但所提出的TR-PIELM因其在许多需要重复求解的瞬态情况下能够实现快速准确的替代预测而变得合理。在这种设置中，时间推进设计有助于在保持ELM的计算优势的同时保持精度。
尽管取得了这些有希望的结果，但仍有一些未来的工作方向：（1）将观测数据或实验数据与Richards方程结合进行数据驱动或逆向渗透分析，在实际应用中具有巨大潜力。（2）可以进一步扩展和优化时间推进策略，以确保PDE求解的计算效率和绝对收敛性。（3）未来的研究将探索改进的网络架构和超参数调整策略，以丰富随机特征基空间，增强方法捕捉非线性行为的能力。（4）将框架扩展到耦合的PDE系统也将是重要的下一步。（5）在本研究中，只考虑了简单的异质条件，即两层土壤系统；更复杂的异质性是未来研究的方向。虽然所提出的TR-PIELM框架预计会保持其效率优势，但其与高度异质介质中数值算法的相对精度仍有待系统评估。因此，解决这些挑战将进一步确立所提方法作为在广泛科学和工程领域中模拟复杂物理过程的强大且通用工具的地位。

CRediT作者贡献声明
Biao Yuan：撰写——原始草稿、可视化、验证、软件、方法论、调查、形式分析、数据整理、概念化。
Roberto Cudmani：撰写——审阅与编辑、资源、项目管理、方法论、调查、概念化。
Wei Yan：撰写——审阅与编辑、资源、项目管理、方法论、调查、概念化。
He Yang：撰写——审阅与编辑、方法论、调查、概念化。
Ana Heitor：撰写——审阅与编辑、监督、概念化。
Yanjie Song：撰写——审阅与编辑、方法论、调查、概念化。
Xiaohui Chen：撰写——审阅与编辑、监督、资源、项目管理、方法论、调查、资金获取、概念化。