一、文献阅读

题目信息

• 题目： Physics-Informed Neural Networks for Modeling Water Flows in a River Channel
• 期刊： IEEE TRANSACTIONS ON ARTIFICIAL INTELLIGENCE
• 作者： Luis Fernando Nazari，Eduardo Camponogara，Laio Oriel Seman
• 发表时间： 2024/03/03
• 文章链接： https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9863669

摘要

洪水灾害对世界各地的影响重大，造成了严重的社会和经济问题，洪水防控策略依赖于水文预报模型，而这些模型主要通过数据驱动方法获得。所以模型需要达到期望的精度需要大量物理参数，而这些参数往往难以确获取；且传统的数据驱动模型往往需要可靠且真实的数据保证模型的数值模拟准确性。所以在以上背景研究下，本论文提出了一个基于物理信息神经网络（Physics-Informed Neural Networks，PINNs）的河道水流替代模型，作者总共设计了三个PINNs模型架构A、B、C，通过调节模型的输入与输出参数以探究在不同场景中模型的适用性。其中在架构C中作者引入了贝叶斯推断，限制参数偏离合理范围，提升物理可解释性。该论文旨在解决传统方法存在的问题，为流域水流的水文建模领域的研究做出贡献。

Abstract

Flood disasters exert profound impacts on populations worldwide, resulting in severe social and economic challenges. Effective flood prevention strategies hinge on hydrological forecasting models, which are predominantly derived through data-driven approaches. Achieving the desired precision in these models necessitates a multitude of physical parameters, which are often difficult to accurately acquire. Moreover, conventional data-driven models typically require reliable and authentic datasets to ensure the accuracy of numerical simulations. Against this backdrop, this paper proposes a surrogate model for river channel water flow based on Physics-Informed Neural Networks (PINNs). The authors developed three distinct PINN architectures—A, B, and C—by adjusting the input and output parameters to investigate their applicability across varying scenarios. Notably, in Architecture C, Bayesian inference is introduced to constrain parameter deviations within physically plausible ranges, thereby enhancing physical interpretability. This study aims to address the limitations inherent in traditional methodologies, contributing to advancements in the field of hydrological modeling for watershed water flows.

创新点

1. 将PINNs应用于河道水流建模，通过融合圣维南方程的物理约束与实测数据，解决了传统方法中数据稀缺和参数依赖的问题；
2. 通过PINNs直接预测水文动力学的关键参数，如：曼宁系数$n$和河床坡度$S_0$。这种方法减少了传统方法中对地区数据的依赖；
3. 使用了基于贝叶斯推理的PINNs训练新方法，解决参数预测的不确定性问题。

网络架构

本文提出了一种基于物理信息神经网络（PINN）的河流水流建模方法，重点解决了传统水文模型中参数识别困难和数据不足的问题。
作者提出了圣维南方程作为物理模型：

关于 $Q(x,t)$ 与 $P_w(x,t)$的形象表示如下图所示：

下面对论文中三种网络架构进行分析：

架构A

1. 输入：
   $$\psi =[x,t,h(x,t),\Delta h_t(x,t),\Delta h_x(x,t)]\in {\mathbb{R}}^5$$
   其中，参数含义如下：
   $x$：空间；
   $t$：时间；
   $h(x,t)$：水位；
   $\Delta h_t$ ：
   数值微分计算的时间导数；
   $\Delta h_x$：空间导数。
2. 输出：
   模型的输出为 $Q$水流速。
3. 损失函数
   PINN的模型架构的损失函数由两部分组成：数据误差和物理误差，即$\text{Loss}={\text{MSE}}_u+{\text{MSE}}_p$
   $MS{E}_u$为数据误差，公式表达为： ${\text{MSE}}_u=\frac{1}{N_u}{\sum }_{i=1}^{N_u}{\left|{\mathcal{NN}}_Q({\psi }_u^i)-Q_u^i\right|}^2$，表示真实值与预测值之间的误差。具体展开为：${\text{MSE}}_p=\frac{1}{N_c}{\sum }_{i=1}^{N_c}\left[{\left|\frac{\partial A}{\partial t}+\frac{\partial {\mathcal{NN}}_Q}{\partial x}\right|}^2+{\left|\frac{\partial {\mathcal{NN}}_Q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{{\mathcal{NN}}_Q^2}{A}\right)+gA\left(\frac{\partial h}{\partial x}+S_f-S_0\right)\right|}^2\right]$

架构B

1. input：
   $\psi =[x,t,h(x,t)]\in {\mathbb{R}}^3$
   其中，参数含义如下：
   $x$：空间；
   $t$：时间；
   $h(x,t)$：水位程高；
2. output：
   输出为： h和Q，分别表示水位高度与流速。
3. Loss function
   损失函数扩展为三部分,分别为：数据误差和圣维南方程质量守恒约束以及圣维南方程动量守恒约束
   $$lossfunction=\frac{1}{N_u}\sum _{i=1}^{N_u}{\left|\mathcal{NN}\left({\psi }_u^i\right)-\left(h^i,Q^i\right)\right|}^2+\frac{1}{N_c}\sum _{i=1}^{N_c}{\left|\frac{\partial A}{\partial h}\frac{\partial }{\partial t}\mathcal{N}{\mathcal{N}}_h\left({\psi }_c^i\right)+\frac{\partial }{\partial x}\mathcal{N}{\mathcal{N}}_Q\left({\psi }_c^i\right)\right|}^2+\frac{1}{N_c}\sum _{i=1}^{N_c}|\frac{\partial }{\partial t}{\mathcal{NN}}_Q\left({\psi }_c^i\right)+\frac{\partial }{\partial x}\frac{{\mathcal{NN}}_Q^2\left({\psi }_c^i\right)}{A\left(x_c^i,t_c^i\right)}+{gA\left(x_c^i,t_c^i\right)\left(\frac{\partial }{\partial x}\mathcal{N}{\mathcal{N}}_h\left({\psi }_c^i\right)+S_f\left(x_c^i,t_c^i\right)-S_0\right)|}^2$$
   其中，$N_c$为配准点（collocation points）的数量，这些点是用来强制满足物理方程的时空坐标；$\mathcal{N}{\mathcal{N}}_h\left({\psi }_c^i\right)$:神经网络预测的水位高度h在配准点$\left(x_c^i,t_c^i\right)$上的值；$\mathcal{N}{\mathcal{N}}_Q\left({\psi }_c^i\right)$：神经网络预测的水流量Q在配准$\left(x_c^i,t_c^i\right)$上的值。
   

架构C

1. input：
   $\psi =[x,t,h(x,t),\Delta h_t(x,t)]\in {\mathbb{R}}^4$
   其中，参数${\mathbb{R}}^4$表示四维空间；
   $x$：空间；
   $t$：时间；
   $h$：水位；
   $\Delta h_t$:表示水面程高对时间t的求导。
2. 输出：
   输出与架构B相同都是$h$和$Q$，输出是一个二维向量。
3. 损失函数
   与架构A和架构B不同的是，架构C在损失函数引入贝叶斯推断，增加惩罚项，其整体的损失函数表示为：
   $$\text{Loss Function}=\underset{\text{数据误差}}{\underbrace{\frac{1}{N_u}\sum _{i=1}^{N_u}{\left|\mathcal{NN}({\psi }_u^i)-[h^i,Q^i]\right|}^2}}+\underset{\text{质量守恒}}{\underbrace{\frac{1}{N_c}\sum _{i=1}^{N_c}{\left|\frac{\partial A}{\partial h}\frac{\partial {\mathcal{NN}}_h}{\partial t}+\frac{\partial {\mathcal{NN}}_Q}{\partial x}\right|}^2}}+\underset{\text{动量守恒}}{\underbrace{\frac{1}{N_c}\sum _{i=1}^{N_c}{\left|\frac{\partial {\mathcal{NN}}_Q}{\partial t}+\frac{\partial }{\partial x}\frac{{\mathcal{NN}}_Q^2}{A}+gA\left(\frac{\partial {\mathcal{NN}}_h}{\partial x}+S_f-S_0\right)\right|}^2}}+\underset{\text{贝叶斯正则化}}{\underbrace{\sum _{i=1}^{N_p}\frac{1}{{\sigma }_i^2}{\left({\lambda }_i-{\mu }_i\right)}^2}}+\underset{\text{导数误差}}{\underbrace{\frac{1}{N_c}\sum _{i=1}^{N_c}{\left|\frac{\partial {\mathcal{NN}}_h}{\partial t}-\Delta h_t\right|}^2}}$$
   引入贝叶斯推断可以参数 $S_0$ 和 $n$ 更靠，贝叶斯推断表示为：$\frac{1}{N_p}{\sum }_{j=1}^{N_p}\frac{1}{{\sigma }_i^2}({\lambda }_i-{\mu }_i{)}^2$
   其中，${\lambda }_i$为待估计参数；$N_p$为参数数量；${\mu }_i$参数的均值；${\sigma }_i^2$为参数的方差。
   

架构A、B、C的比较以及适用场景

架构	对输入参数的依赖车程度	物理可解释性	适用场景
A	高	高	密集监测场景
B	低	低	监测站稀疏的实际场景
C	中	最高	需平衡数据与物理定律的场景

实验

1. PINNs用于河道水流建模的初步实验结果
   在模拟环境下，以圣维南方程配置构建的系统中进行实验，目的是在评估PINNs预测河道水流的有效性。实验分两个阶段，第一阶段假设PDE参数和水位函数已知，验证PINNs的数据同化能力；第二阶段聚焦于识别定义PDE的向量参数，以描述水流动力学并捕捉河道内在特征。
   通过数值解圣维南方程得到训练数据和配置点集，采用如架构A进行实验，结果表明PINN方法在系统识别方面，使用相对少的数据点，在所有实验中达到相当的精度和强大的数据同化能力，验证的均方误差在0.00648% - 0.0466%之间。
   

2. 不同架构下的建模结果
   架构A
   在将$S_0$和$n$作为优化变量加入损失函数后，经过一定训练轮次，架构A在所有实验中的预测误差（MSE）低于0.305%，最佳情况MSE为0.0263%。
   
   基于PINNs发现的摩擦斜率和曼宁系数近似于圣维南方程中的实际参数，表明PINNs通过利用微分方程模型知识实现有效正则化。尽管没有理论保证神经网络模型和参数发现的联合学习收敛到全局最优，但在一定条件下可达到较好预测精度。
   
   架构B
   为解决数据点缺乏问题提出架构B，其输入变量仅包含距离、时间和河流水位，输出除水流率外还尝试学习水位函数用于自动微分计算物理正则化中的偏导数。
   
   实验结果显示架构B发现真实PDE参数的能力不如架构A，但能较好地近似系统行为，MSE验证指数较低。
   
   架构C
   为改进参数识别提出架构C，采用新结构并受益于新的损失函数，引入贝叶斯推理到训练过程。
   架构C添加$\Delta h_t$作为输入变量，新的损失函数包含对参数偏离均值的惩罚项和对输出自动微分与数值导数差异的惩罚项。
   实验结果表明，架构C在系统行为的泛化方面，MSE验证指数在0.0762% - 0.298%之间，相比架构B有所改进，在参数估计方面也有显著提升，能以较好的初始先验值逼近真实参数并保持物理特性。对比固定方差向量的情况，虽然参数有显著改进，但会导致系统识别能力的损失。
   

3. 伊塔雅伊米林河中的实验结果
   以伊塔雅伊河的伊塔雅伊米林河河段为研究对象，该区域有水位和流量测量数据，以及政府提供的地理研究数据。
   
   利用2020年1月10 - 14日一次洪水事件的960个样本数据训练PINN模型，采用架构C并加入贝叶斯推理正则化的实验结果显示，架构为八层隐藏层、每层30个神经元、80000个配置点时，数据集同化效果最佳（MSE = 4.126×10⁻³），发现的床坡和曼宁系数参数，其中曼宁系数收敛到接近文献中的值，床坡参数收敛到区域研究中的极限值内。
   
   
   代码如下：

# 导入必要的库
import torch  # 用于张量计算和神经网络搭建
import torch.nn as nn  # 提供神经网络层和激活函数
import numpy as np  # 用于数据转换和绘图
import matplotlib.pyplot as plt  # 用于可视化训练结果
from time import time  # 用于测量训练时间
from pytorch_lbfgs import LBFGS  # 导入支持GPU的L-BFGS优化器# 设置随机种子，确保结果可重复
torch.manual_seed(42)# 检查并设置计算设备，确保使用GPU
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
if not torch.cuda.is_available():raise RuntimeError("GPU not available. This code requires CUDA support.")
print(f"Using device: {device} ({torch.cuda.get_device_name(0)})")# 定义Saint-Venant方程
def saint_venant_equations(h, Q, x, t, S0, n, g=9.81, eps=1e-6):"""计算Saint-Venant方程的残差，用于物理约束"""width = 2.0  # 矩形横截面宽度2mA = h * width + eps  # 横截面积，添加偏移避免除零P_w = 2 * h + width + eps  # 湿周长，添加偏移R = A / P_w  # 水力半径# 自动微分计算偏导数h_t = torch.autograd.grad(h.sum(), t, create_graph=True)  # 计算 ∂h/∂th_t = h_t[0] if h_t is not None else torch.zeros_like(h)  # 若梯度为空，返回零张量h_x = torch.autograd.grad(h.sum(), x, create_graph=True)  # 计算 ∂h/∂xh_x = h_x[0] if h_x is not None else torch.zeros_like(h)  # 若梯度为空，返回零张量Q_x = torch.autograd.grad(Q.sum(), x, create_graph=True)  # 计算 ∂Q/∂xQ_x = Q_x[0] if Q_x is not None else torch.zeros_like(Q)  # 若梯度为空，返回零张量Q_t = torch.autograd.grad(Q.sum(), t, create_graph=True)  # 计算 ∂Q/∂tQ_t = Q_t[0] if Q_t is not None else torch.zeros_like(Q)  # 若梯度为空，返回零张量Q2_A = Q ** 2 / (A + eps)  # 计算 Q^2/A，添加偏移避免除零Q2_A_x = torch.autograd.grad(Q2_A.sum(), x, create_graph=True)  # 计算 ∂(Q^2/A)/∂xQ2_A_x = Q2_A_x[0] if Q2_A_x is not None else torch.zeros_like(Q)  # 若梯度为空，返回零张量S_f = (n ** 2 * Q * torch.abs(Q)) / (A ** 2 * R ** (4 / 3) + eps)  # 摩擦坡度，添加偏移f1 = (width * h_t) + Q_x  # 质量守恒方程f2 = Q_t + Q2_A_x + g * A * (h_x + S_f - S0)  # 动量守恒方程return f1, f2  # 返回残差# 数据生成函数
def generate_data(N_u=1200, N_c=150000, N_val=100):"""生成训练数据、配准点和验证点，模拟2km河道洪水波"""t = torch.linspace(0, 12000, N_u // 2, device=device, dtype=torch.float32).reshape(-1, 1)  # 生成时间序列x_train = torch.cat([torch.zeros_like(t), 2000 * torch.ones_like(t)], dim=0)  # 设置x=0m和x=2000mt_train = torch.cat([t, t], dim=0)  # 拼接完整时间序列h_train = torch.zeros_like(t_train)  # 初始化水深数组mask1 = t_train <= 6000  # 上升阶段mask2 = (t_train > 6000) & (t_train <= 9000)  # 下降阶段mask3 = t_train > 9000  # 平稳阶段h_train[mask1] = (1.000735) ** (t_train[mask1] / 10) * 0.4  # 上升阶段水深h_ta = (1.000735) ** (6000 / 10) * 0.4  # t=6000s参考值h_train[mask2] = (0.9985) ** ((t_train[mask2] - 6000) / 10) * h_ta  # 下降阶段水深h_train[mask3] = 0.4  # 平稳阶段水深A0 = h_train * 2.0 + 1e-6  # 初始横截面积P0 = 2 * h_train + 2.0 + 1e-6  # 初始湿周长Q_train = ((0.0025) ** 0.5 / 0.02) * (A0 ** (5 / 3)) / (P0 ** (2 / 3))  # 稳态流量x_c = torch.rand(N_c, 1, device=device, dtype=torch.float32) * 2000  # 配准点空间t_c = torch.rand(N_c, 1, device=device, dtype=torch.float32) * 12000  # 配准点时间x_val = torch.tensor([500, 1000, 1500], device=device, dtype=torch.float32).repeat(N_val).reshape(-1, 1)  # 验证点位置t_val = torch.linspace(0, 12000, N_val, device=device, dtype=torch.float32).repeat(3).reshape(-1, 1)  # 验证点时间return (x_train.requires_grad_(True), t_train.requires_grad_(True), h_train, Q_train), \(x_c.requires_grad_(True), t_c.requires_grad_(True)), \(x_val.requires_grad_(True), t_val.requires_grad_(True))# 定义PINN模型（架构C）
class PINN_C(nn.Module):def __init__(self, hidden_layers=8, neurons=20):"""定义架构C的神经网络：输入[x, t, Δh_t]，输出[h, Q]"""super(PINN_C, self).__init__()layers = [nn.Linear(3, neurons), nn.Tanh()]  # 输入层：3维 -> neurons，使用tanh激活for _ in range(hidden_layers - 1):layers += [nn.Linear(neurons, neurons), nn.Tanh()]  # 添加隐层layers += [nn.Linear(neurons, 2)]  # 输出层：neurons -> 2 (h, Q)self.net = nn.Sequential(*layers).to(device)  # 构建网络并移到GPUself.S0 = nn.Parameter(torch.tensor([0.002], device=device, dtype=torch.float32))  # 初始化床坡self.n = nn.Parameter(torch.tensor([0.03], device=device, dtype=torch.float32))  # 初始化曼宁系数def forward(self, x, t, dh_dt):"""前向传播，输入[x, t, Δh_t]，输出[h, Q]"""inputs = torch.cat([x, t, dh_dt], dim=1)  # 拼接输入向量out = self.net(inputs)  # 通过神经网络h = torch.relu(out[:, 0:1]) + 0.01  # 输出水深，确保非负并加偏移Q = torch.relu(out[:, 1:2]) + 0.01  # 输出流量，确保非负并加偏移return h, Q# 损失函数
def loss_fn_C(model, x_u, t_u, h_u, Q_u, x_c, t_c, mu, sigma):"""计算损失函数，包括数据项、物理项、正则化项和导数惩罚项"""dh_dt_u = torch.zeros_like(h_u)  # 初始化训练数据的 Δh_tdh_dt_u[1:] = (h_u[1:] - h_u[:-1]) / 10  # 计算数值导数dh_dt_u[0] = dh_dt_u[1]  # 边界填充h_pred, Q_pred = model(x_u, t_u, dh_dt_u)  # 预测训练点的 h 和 Qmse_u = torch.mean((h_pred - h_u) ** 2 + (Q_pred - Q_u) ** 2)  # 计算数据误差dh_dt_c = torch.zeros_like(x_c)  # 初始化配准点的 Δh_th_c, Q_c = model(x_c, t_c, dh_dt_c)  # 预测配准点的 h 和 Qf1, f2 = saint_venant_equations(h_c, Q_c, x_c, t_c, model.S0, model.n)  # 计算PDE残差mse_p = torch.mean(f1 ** 2) + torch.mean(f2 ** 2)  # 计算物理误差reg = (1 / (sigma[0] ** 2 + 1e-6)) * (model.S0 - mu[0]) ** 2 + \(1 / (sigma[1] ** 2 + 1e-6)) * (model.n - mu[1]) ** 2  # 参数正则化项h_t = torch.autograd.grad(h_c.sum(), t_c, create_graph=True)  # 计算配准点 ∂h/∂th_t = h_t[0] if h_t is not None else torch.zeros_like(h_c)  # 若梯度为空，返回零张量penalty = torch.mean((h_t - dh_dt_c) ** 2)  # 计算导数惩罚项return mse_u + mse_p + reg + penalty  # 返回总损失# 验证MSE计算
def compute_validation_mse(model, x_val, t_val, h_true, Q_true):"""计算验证点的MSE"""dh_dt_val = torch.zeros_like(h_true)  # 初始化验证点 Δh_tdh_dt_val[1:] = (h_true[1:] - h_true[:-1]) / 10  # 计算数值导数dh_dt_val[0] = dh_dt_val[1]  # 边界填充with torch.no_grad():  # 无梯度计算h_pred, Q_pred = model(x_val, t_val, dh_dt_val)  # 预测验证点 h 和 Qmse_val = torch.mean((h_pred - h_true) ** 2 + (Q_pred - Q_true) ** 2)  # 计算MSEreturn mse_val.item()  # 返回标量值# 训练函数
def train_model(model, x_u, t_u, h_u, Q_u, x_c, t_c, x_val, t_val, h_val, Q_val):"""训练模型：2000次ADAM + 5000次L-BFGS，全部在GPU上"""optimizer_adam = torch.optim.Adam(model.parameters(), lr=0.0001)  # 初始化ADAM优化器optimizer_lbfgs = LBFGS(model.parameters(), lr=0.1, max_iter=20)  # 初始化GPU支持的L-BFGSmu = torch.tensor([0.002, 0.03], device=device, dtype=torch.float32)  # 初始化先验均值sigma = torch.tensor([4.17e-4, 5e-3], device=device, dtype=torch.float32)  # 初始化先验方差S0_history, n_history, losses, val_mses = [], [], [], []  # 初始化历史记录列表print("Training with ADAM on GPU (2000 epochs)...")start_time = time()  # 记录ADAM阶段开始时间for epoch in range(2000):  # 2000次ADAM迭代optimizer_adam.zero_grad()  # 清零梯度loss = loss_fn_C(model, x_u, t_u, h_u, Q_u, x_c, t_c, mu, sigma)  # 计算总损失if torch.isnan(loss):  # 检查是否出现NaNprint(f"NaN detected at epoch {epoch}")breakloss.backward()  # 反向传播torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)  # 梯度裁剪optimizer_adam.step()  # 更新参数S0_history.append(model.S0.item())  # 记录床坡n_history.append(model.n.item())  # 记录曼宁系数losses.append(loss.item())  # 记录损失mu[0] = torch.tensor(S0_history, device=device, dtype=torch.float32).mean()  # 更新S0均值mu[1] = torch.tensor(n_history, device=device, dtype=torch.float32).mean()  # 更新n均值sigma[0] = max((max(S0_history) - min(S0_history)) / 6, 1e-6)  # 更新S0方差sigma[1] = max((max(n_history) - min(n_history)) / 6, 1e-6)  # 更新n方差val_mse = compute_validation_mse(model, x_val, t_val, h_val, Q_val)  # 计算验证MSEval_mses.append(val_mse)  # 记录验证MSEif epoch % 200 == 0:  # 每200次打印结果print(f"Epoch {epoch}, Loss: {loss.item():.6f}, Val MSE: {val_mse:.6f}, "f"S0: {model.S0.item():.6f}, n: {model.n.item():.6f}")print(f"ADAM Training Time: {time() - start_time:.2f}s")  # 打印ADAM阶段时间print("Training with L-BFGS on GPU (5000 epochs)...")start_time = time()  # 记录L-BFGS阶段开始时间def closure():  # 定义L-BFGS闭包函数optimizer_lbfgs.zero_grad()  # 清零梯度loss = loss_fn_C(model, x_u, t_u, h_u, Q_u, x_c, t_c, mu, sigma)  # 计算总损失if torch.isnan(loss):  # 检查是否出现NaNprint("NaN detected in L-BFGS")return lossloss.backward()  # 反向传播return lossfor epoch in range(5000):  # 5000次L-BFGS迭代loss = optimizer_lbfgs.step(closure)  # 执行一步优化if torch.isnan(loss):  # 检查是否出现NaNprint(f"NaN detected at epoch {epoch + 2000}")breakS0_history.append(model.S0.item())  # 记录床坡n_history.append(model.n.item())  # 记录曼宁系数losses.append(loss.item())  # 记录损失mu[0] = torch.tensor(S0_history, device=device, dtype=torch.float32).mean()  # 更新S0均值mu[1] = torch.tensor(n_history, device=device, dtype=torch.float32).mean()  # 更新n均值sigma[0] = max((max(S0_history) - min(S0_history)) / 6, 1e-6)  # 更新S0方差sigma[1] = max((max(n_history) - min(n_history)) / 6, 1e-6)  # 更新n方差val_mse = compute_validation_mse(model, x_val, t_val, h_val, Q_val)  # 计算验证MSEval_mses.append(val_mse)  # 记录验证MSEif (epoch + 2000) % 500 == 0:  # 每500次打印结果print(f"Epoch {epoch + 2000}, Loss: {loss.item():.6f}, Val MSE: {val_mse:.6f}, "f"S0: {model.S0.item():.6f}, n: {model.n.item():.6f}")print(f"L-BFGS Training Time: {time() - start_time:.2f}s")  # 打印L-BFGS阶段时间return losses, val_mses, S0_history, n_history  # 返回训练历史# 可视化结果
def plot_results(model, x_u, t_u, h_u, Q_u, x_val, t_val, losses, val_mses, S0_history, n_history):"""绘制训练结果"""x_u, t_u, h_u = x_u.to('cpu'), t_u.to('cpu'), h_u.to('cpu')  # 数据移到CPU用于绘图dh_dt_u = torch.zeros_like(h_u)  # 初始化 Δh_tdh_dt_u[1:] = (h_u[1:] - h_u[:-1]) / 10  # 计算数值导数dh_dt_u[0] = dh_dt_u[1]  # 边界填充with torch.no_grad():  # 无梯度计算model = model.to('cpu')  # 模型移到CPUh_pred, Q_pred = model(x_u, t_u, dh_dt_u)  # 预测 h 和 Qh_u, Q_u = h_u.numpy(), Q_u.numpy()  # 转换为NumPy数组h_pred, Q_pred = h_pred.numpy(), Q_pred.numpy()  # 同上t_u = t_u.numpy()  # 同上plt.figure(figsize=(15, 10))  # 设置图形大小plt.subplot(2, 2, 1)  # 子图1：损失曲线plt.plot(losses, label='Loss')  # 绘制总损失plt.plot(val_mses, label='Validation MSE')  # 绘制验证MSEplt.xlabel('Epoch')  # x轴标签plt.ylabel('Value')  # y轴标签plt.yscale('log')  # 对数刻度plt.legend()  # 添加图例plt.subplot(2, 2, 2)  # 子图2：流量预测plt.plot(t_u[:len(t_u) // 2], Q_u[:len(t_u) // 2], label='True Q (x=0)')  # 真实流量plt.plot(t_u[:len(t_u) // 2], Q_pred[:len(t_u) // 2], '--', label='Pred Q (x=0)')  # 预测流量plt.xlabel('Time (s)')  # x轴标签plt.ylabel('Flow Rate (m^3/s)')  # y轴标签plt.legend()  # 添加图例plt.subplot(2, 2, 3)  # 子图3：水深预测plt.plot(t_u[:len(t_u) // 2], h_u[:len(t_u) // 2], label='True h (x=0)')  # 真实水深plt.plot(t_u[:len(t_u) // 2], h_pred[:len(t_u) // 2], '--', label='Pred h (x=0)')  # 预测水深plt.xlabel('Time (s)')  # x轴标签plt.ylabel('Water Depth (m)')  # y轴标签plt.legend()  # 添加图例plt.subplot(2, 2, 4)  # 子图4：参数估计plt.plot(S0_history, label='S0')  # 绘制床坡估计plt.plot(n_history, label='n')  # 绘制曼宁系数估计plt.axhline(y=0.0025, color='r', linestyle='--', label='True S0')  # 真实床坡plt.axhline(y=0.02, color='g', linestyle='--', label='True n')  # 真实曼宁系数plt.xlabel('Epoch')  # x轴标签plt.ylabel('Parameter Value')  # y轴标签plt.legend()  # 添加图例plt.tight_layout()  # 调整布局plt.show()  # 显示图形# 主程序
if __name__ == "__main__":(x_u, t_u, h_u, Q_u), (x_c, t_c), (x_val, t_val) = generate_data()  # 生成数据h_val = h_u[:300]  # 设置验证集水深Q_val = Q_u[:300]  # 设置验证集流量model = PINN_C().to(device)  # 初始化模型并移到GPUlosses, val_mses, S0_history, n_history = train_model(model, x_u, t_u, h_u, Q_u, x_c, t_c, x_val, t_val, h_val,Q_val)  # 训练模型plot_results(model, x_u, t_u, h_u, Q_u, x_val, t_val, losses, val_mses, S0_history, n_history)  # 可视化结果print(f"Final S0: {model.S0.item():.6f}, True S0: 0.0025")  # 打印最终床坡估计print(f"Final n: {model.n.item():.6f}, True n: 0.02")  # 打印最终曼宁系数估计print(f"Final Validation MSE: {val_mses[-1]:.6f}")  # 打印最终验证MSE

结论

本论文将PINNs作为河流水流率的替代水文模型，研究其在模拟环境中的性能。结果表明PINNs可同化数据并发现控制系统动力学的微分方程参数，能规避概念模型中参数估计的劣势和经验方法中数据缺乏的问题。且在不同条件下，架构A和架构C的PINNs模型分别有较好表现。

后续展望

作者认为未来可将PINNs应用于解决流域的实时系统。此外，进一步探索PINNs在不同真实世界场景中的应用和深入研究物理和数据学习平衡任务的理论。之后可以优化架构C中参数方差的值以提高系统预测能力。