2025.2.23机器学习笔记：PINN文献阅读

一、文献阅读

题目信息

• 题目： 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的比较以及适用场景

实验

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  # 矩形横截面宽度2m
    A = 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/∂t
    h_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/∂x
    h_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/∂x
    Q_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/∂t
    Q_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)/∂x
    Q2_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=2000m
    t_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)  # 构建网络并移到GPU
        self.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_t
    dh_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 和 Q
    mse_u = torch.mean((h_pred - h_u) ** 2 + (Q_pred - Q_u) ** 2)  # 计算数据误差

    dh_dt_c = torch.zeros_like(x_c)  # 初始化配准点的 Δh_t
    h_c, Q_c = model(x_c, t_c, dh_dt_c)  # 预测配准点的 h 和 Q
    f1, 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/∂t
    h_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_t
    dh_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 和 Q
        mse_val = torch.mean((h_pred - h_true) ** 2 + (Q_pred - Q_true) ** 2)  # 计算MSE
    return 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-BFGS

    mu = 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):  # 检查是否出现NaN
            print(f"NaN detected at epoch {epoch}")
            break
        loss.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)  # 计算验证MSE
        val_mses.append(val_mse)  # 记录验证MSE

        if 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):  # 检查是否出现NaN
            print("NaN detected in L-BFGS")
            return loss
        loss.backward()  # 反向传播
        return loss

    for epoch in range(5000):  # 5000次L-BFGS迭代
        loss = optimizer_lbfgs.step(closure)  # 执行一步优化
        if torch.isnan(loss):  # 检查是否出现NaN
            print(f"NaN detected at epoch {epoch + 2000}")
            break
        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)  # 计算验证MSE
        val_mses.append(val_mse)  # 记录验证MSE

        if (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_t
    dh_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')  # 模型移到CPU
        h_pred, Q_pred = model(x_u, t_u, dh_dt_u)  # 预测 h 和 Q

    h_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')  # 绘制验证MSE
    plt.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)  # 初始化模型并移到GPU
    losses, 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中参数方差的值以提高系统预测能力。