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

一、文献阅读

题目信息

• 题目： Physics-Informed Neural Networks of the Saint-Venant Equations for Downscaling a Large-Scale River Model
• 期刊： Water Resources Research
• 作者： Dongyu Feng, Zeli Tan, QiZhi He
• 发表时间： 2023/12/11
• 文章链接： https://arxiv.org/pdf/2210.03240v1

摘要

沿海地区人口增长，人类面临飓风和洪水等自然灾害风险增大。气候变化加剧极端风暴潮、降水及海平面上升，使洪水风险进一步提升，研究潮汐河流动力学对减轻洪灾风险至关重要。按照传统的方法，大规模河流模型是气候变化研究的重要工具，但在模拟局部洪水过程时存在不足。其物理可解释性和网格分辨率低，无法解析洪水泛滥的详细信息；统计和动力降尺度方法，但在河流建模中应用有限。传统线性插值降尺度方法无法解决网格单元内空间变化的水流问题。为解决以上问题，本文提出基于物理信息神经网络的机器学习框架，用于模拟大尺度河流模型在沿海地区的亚网格尺度水流。首先展示PINN能同化各类观测并求解一维圣维南方程，通过多个合成案例研究评估其性能，结果表明PINN在水深模拟上精度良好。针对风暴潮和潮汐引发的洪水波传播，提出基于傅里叶特征嵌入的神经网络架构。此外，研究表明PINN降尺度能通过同化观测数据产生更合理的亚网格解，优于简单线性插值，为改进大尺度模型模拟精细尺度沿海过程能力提供了有前景的途径。

Abstract

Population growth in coastal regions has heightened human exposure to natural disasters such as hurricanes and floods. Climate change exacerbates the risks of extreme storm surges, precipitation, and sea-level rise, further amplifying flood hazards. Investigating tidal river dynamics is critical for mitigating these risks. Traditionally, large-scale river models have served as essential tools in climate change research. However, they exhibit limitations in simulating localized flood processes due to their low physical interpretability and coarse grid resolution, which hinder the detailed resolution of flood inundation dynamics. Statistical and dynamical downscaling methods, while useful, have seen limited application in river modeling. Conventional linear interpolation downscaling approaches fail to address spatial variations in water flow within grid cells. To overcome these challenges, this study proposes a machine learning framework based on Physics-Informed Neural Networks (PINNs) to simulate subgrid-scale water flows in large-scale river models for coastal regions. We first demonstrate that PINNs can assimilate diverse observational data and solve the one-dimensional Saint-Venant equations. Performance evaluation through multiple synthetic case studies reveals that PINNs achieve high accuracy in water depth simulations. For flood wave propagation induced by storm surges and tides, we introduce a neural network architecture incorporating Fourier feature embedding. Furthermore, the study shows that PINN-based downscaling, by assimilating observational data, generates more physically consistent subgrid solutions compared to simple linear interpolation. This approach offers a promising pathway for enhancing the capability of large-scale models to simulate fine-scale coastal processes.

创新点

1. 作者通过基于PINN的数据同化降尺度方法，解决大尺度河流模型在海岸地区子网格尺度河流动力学的变异性问题。该方法能融合多种观测数据，无需修改数值算法或细化网格分辨率，且不受网格限制，为降尺度研究提供新途径。

    降尺度是指从一个粗糙的、低分辨率的数据或模型结果，推导出更精细的、高分辨率的结果的过程。简单来说就是把模糊的大图变成清晰的小图。
   

2. 作者提出了傅里叶特征嵌入的方法，在输入层前将坐标𝑥和𝑡映射到高维傅里叶特征空间。论文特别针对时间坐标𝑡，以适应潮汐的周期性。下面对其进行说明：
   傅里叶变换的核心思想是： 任何周期性的信号都可以分解成一堆不同频率的正弦波和余弦波。因为正弦和余弦是周期函数，可以描述有规律的波动。
   其数学表达如下：
   $$\gamma (\nu )=\left[\begin{array}{c}\cos (2\pi B\nu )\\ \sin (2\pi B\nu )\end{array}\right],\nu =[x,t{]}^T,B\sim N(0,s)$$
   其中：
   $\nu =[x,t{]}^T$输入是空间坐标 𝑥和时间 𝑡，这里是一个二维向量；
   B是一个随机矩阵，里面的值从高斯分布 𝑁(0, 𝑠)中抽取，𝑠 是标准差；
   $\cos (2\pi B\nu )$与$\sin (2\pi B\nu )$则是对输入进行正弦和余弦变换。
   论文中作者通过这个特性将傅里叶特征嵌入到时间中，将来描述潮汐（因为潮汐具有周期性）,具体表示为：
   $$\begin{array}{l}{\gamma }_t^{(i)}(t)=\left[\begin{array}{c}\cos \left(2\pi {\mathbf{B}}_t^{(i)}t\right)\\ \sin \left(2\pi {\mathbf{B}}_t^{(i)}t\right)\end{array}\right],\text{ for }i=1,2,\dots \end{array}$$

网络架构

如图(a)所示，在本片论文中PINN框架由一个深度神经网络与SVE结合的架构。
在输出结果在作者还强调数据同化过程（所谓数据同化，就是让模型预测值不断接近真实值的一个过程，而不是某一个步骤。比如，在算Loss中就算是一个数据同化的步骤）。
下面我们来具体分析一下这个模型的架构：
输入：空间坐标 𝑥 和时间 𝑡 。
结构： 在图示中NN的架构为2个隐藏层和每层4个神经元，其中隐藏层使用激活函数 𝜎
其中𝜎为：$\tanh =\frac{e^x-e^{-x}}{e^x+e^{-x}}$。
输出：预测值 𝑢 流速和ℎ水深，由橙色节点表示。此外，输入数据来源于数值解、现场测站数据、卫星提供的水深快照，这些数据通过损失函数约束神经网络，确保预测接近真实值。
物理约束：

然后通过自动微分计算偏导数并代入SVE检查残差。

损失函数： 损失函数定义为
$$J(\theta )=J_f(\theta )+\sum _j{w}_j{J}_j(\theta )$$
1、 $J_f(\theta )$：SVE残差损失。其中$J_f(\theta )=\frac{1}{N_f}{\sum }_{i=1}^{N_f}\left|r_f^2(x^i,t^i,\theta )\right|$
$N_f$为配置点；$r_f$是SVE的残差表示；$\theta$是神经网络的参数；$x^i,t^i$是第 𝑖个配置点的空间坐标和时间坐标。

2、${\sum }_j{w}_j{J}_j(\theta )$：边界条件、初始条件损失、观测数据的误差等。
具体如下：
边界条件损失$J_{B{C}_h}$、$J_{B{C}_u}$：
边界条件损失，分别衡量流速 𝑢和水深 ℎ在边界上的预测值与已知边界条件𝑢和ℎ的偏差。值越小，说明边界符合要求。
其中：
$J_{B{C}_u}(\theta )=\frac{1}{N_{B{C}_u}}{\sum }_{i=1}^{N_{B{C}_u}}{\left|\hat{u}(x^i,t^i,\theta )-u(x^i,t^i)\right|}^2$
$J_{B{C}_h}(\theta )=\frac{1}{N_{B{C}_h}}{\sum }_{i=1}^{N_{B{C}_h}}{\left|\hat{h}(x^i,t^i,\theta )-h(x^i,t^i)\right|}^2,x\in {\Omega }_{BC},t\in [0,T]$
$N_{B{C}_u}$与$N_{B{C}_h}$是边界条件点的总数；${\sum }_{i=1}^{N_{B{C}_h}}$与${\sum }_{i=1}^{N_{B{C}_u}}$对边界上的所有点求和。

初始条件损失$J_{S_u}$和快照损失$J_{S_h}$：
初始条件或快照损失，分别衡量$\hat{u}$流速和$\hat{h}$水深在初始时间 t=0 或特定快照时刻的预测值与已知值的偏差。
$J_{S_u}(\theta )=\frac{1}{N_{S_u}}{\sum }_{i=1}^{N_{S_u}}{\left|\hat{u}(x^i,0,\theta )-u(x^i,0)\right|}^2$
$J_{S_h}(\theta )=\frac{1}{N_{S_h}}{\sum }_{i=1}^{N_{S_h}}{\left|\hat{h}(x^i,0,\theta )-h(x^i,0)\right|}^2,x\in {\Omega }_S$
其中$N_{S_u}$与$N_{S_h}$是初始条件或快照点的总数；$x\in {\Omega }_S$定义了河道距离$x$的范围。

	快照条件是指在特定时间点上对系统状态。如，水深 ℎ或流速u的已知值或观测值。此外，快照条件可以是t=0，也可以是后续时间点的观测值，提供额外的时空约束。

观测数据损失$J_{ob{s}_u}$与$J_{ob{s}_h}$
$J_{ob{s}_u}(\theta )=\frac{1}{N_{ob{s}_u}}{\sum }_{i=1}^{N_{ob{s}_u}}{\left|\hat{u}(x^i,t^i,\theta )-u(x^i,t^i)\right|}^2$
$J_{ob{s}_h}(\theta )=\frac{1}{N_{ob{s}_h}}{\sum }_{i=1}^{N_{ob{s}_h}}{\left|\hat{h}(x^i,t^i,\theta )-h(x^i,t^i)\right|}^2,x\in {\Omega }_{obs},t\in [0,T]$
其中$N_{ob{s}_u}$与$N_{ob{s}_h}$观测点数据总数；在$x\in {\Omega }_{obs},t\in [0,T]$中，${\Omega }_{obs}$为河道范围，T为模拟周期。

图(b)为一个大尺度网格单元内的子网格降尺度过程，包括时间步和空间位置，其目的是在粗网格Δ𝑥内生成更详细的水深和流速解。其中：
时间从$t_n$到$t_{n+1}$，Δt为时间步长。
空间从$x_i$到$x_{x+1}$，$x_{obs}^1$与$x_{obs}^2$为检测站获取的信息。

实验

本文通过六个实验，研究基于PINN的数据同化降尺度方法，以解决大型河流模型在沿海地区河流动力学的亚网格变异性问题。其中，实验1和2分析的是洪水平原，测试PINN求解SVE能力；实验3和4分析的是开阔河道，验证动态水流模拟；实验5和6分析的潮汐河流，引入了傅里叶嵌入来评估风暴潮和潮汐影响下的降尺度效果。
实验结果及分析如下：

1. 在实验1和2中，PINN在解决简化的圣维南方程（SVE）时展现出较高准确性，能很好地模拟洪水在洪泛平原的传播，其洪水波传播模式和淹没前沿形状与解析解或RK4解接近，相对L2误差和RMSE较小。这表明在有足够观测数据和地形信息时，PINN可用于计算亚网格尺度的洪泛平原淹没情况，提供更详细的淹没地图。
   
2. 在实验3和4中，PINN解与HEC - RAS参考解在水深方面吻合良好，仅在沿河道剖面的上游端附近有较小模拟偏差，L2误差和RMSE也表明其性能良好。在案例4中，编码傅里叶特征的PINN解能够捕捉到下游潮汐引起的周期性变化，相比标准架构的PINN解，其ϵh和RMSE更小。但由于潮汐和边界条件产生的高振荡行为，案例4的预测准确性比案例3差，意味着PINN训练需要更多物理约束和观测数据。
   
3. 在实验5和6中，与线性插值相比，PINN在解析河道地形和亚网格尺度的可变流态方面表现更优。线性插值在无原位数据区域的降尺度解较差，无法再现时空模式和沿河道剖面，而PINN解在整个亚网格区域具有更清晰的流动动力学，特别是在亚临界流占主导且水深较大的中间部分，误差远小于线性插值解。
   在案例6中，尽管两种降尺度方法随着原位数据增加性能都有所提升，但PINN解更准确且对观测数据依赖性小。线性插值在缺乏数据的上游区域高估潮汐影响，而PINN能合理捕捉潮汐相位和振幅，不过在潮汐传播尖端附近因流量 - 潮汐相互作用存在较大偏差。总体而言，PINN降尺度在可变流态条件下优于线性插值。
   
   此外，实验6结果表明，同化观测数据的可用性对PINN解的准确性至关重要，有限的配置点下，状态变量时空变化增加会降低PINN性能，但PINN在降尺度一维河道流时对观测的依赖小于线性插值。提出的傅里叶特征嵌入时间t足以编码潮汐物理特征和处理边界条件的周期性，但当沿河道流变化频繁时，PINN性能会下降，因为傅里叶变换难以直接编码空间变量，且高频与低频过程的非线性相互作用在观测有限时难以很好再现。不过，随着同化观测数据量增加，PINN准确性会提高。
   
   代码如下：

import numpy as np
# 导入插值函数
from scipy.interpolate import interp1d
# 导入 Matplotlib 用于绘图
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1 import make_axes_locatable
import h5py
from datetime import datetime
import time
# 导入拉丁超立方采样工具
from pyDOE import lhs
import tensorflow as tf
import pickle
import os
import warnings# 忽略警告信息以保持输出清洁
warnings.filterwarnings("ignore")
# 设置 TensorFlow 随机种子为 1234 确保可重复性
tf.set_random_seed(1234)
# 设置 NumPy 随机种子为 1234 确保可重复性
np.random.seed(1234)# 定义 SVE 模型类，替代缺失模块
class SVE:DTYPE = tf.float32  # 定义 TensorFlow 数据类型为 float32def __init__(self, X_h_IC, X_u_BC, X_h_BC, X_u_obs, X_h_obs, X_f, h_IC, u_BC, h_BC, u_obs, h_obs,layers, lb, ub, S, b, X_star, u_star, h_star, lr=5e-4, ExistModel=0, uhDir='',wDir='', wmffDir='', useObs=True, use_ff=False):"""初始化 PINN 模型用于 SVE。参数:X_h_IC, X_u_BC, X_h_BC, X_u_obs, X_h_obs, X_f: 输入坐标 (x, t)h_IC, u_BC, h_BC, u_obs, h_obs: 初始条件、边界条件和观测值layers: 神经网络架构lb, ub: 域边界S: 渠道坡度b: 渠道宽度X_star, u_star, h_star: 测试数据lr: 学习率ExistModel: 加载已有模型标志 (0: 新建, 1/2: 加载)uhDir, wDir, wmffDir: 权重保存/加载路径useObs: 是否使用观测数据use_ff: 是否使用傅里叶特征嵌入"""self.count = 0  # 初始化回调计数器，用于跟踪训练迭代次数self.lb, self.ub = lb, ub  # 设置输入数据的下界和上界self.S = tf.constant(S, dtype=self.DTYPE) if isinstance(S, (int, float)) else S  # 将坡度转换为 TensorFlow 张量self.b = tf.constant(b, dtype=self.DTYPE)  # 将渠道宽度转换为 TensorFlow 张量self.useObs, self.use_ff = useObs, use_ff  # 设置是否使用观测数据和傅里叶特征的标志self.X_star, self.u_star, self.h_star = X_star, u_star, h_star  # 存储测试数据 (x, t, u, h)self.beta = 0.9  # 自适应权重更新因子self.adaptive_constants = {  # 初始化自适应权重字典'bcs_u': np.array(1.0), 'bcs_h': np.array(1.0),  # 边界条件权重'ics_h': np.array(1.0), 'obs_u': np.array(1.0),  # 初始条件和观测权重'obs_h': np.array(1.0)}# 准备输入数据self.x_h_IC, self.t_h_IC = X_h_IC[:, 0:1], X_h_IC[:, 1:2]  # 提取初始条件 x 和 t 坐标self.x_u_BC, self.t_u_BC = X_u_BC[:, 0:1], X_u_BC[:, 1:2]  # 提取流速边界条件 x 和 t 坐标self.x_h_BC, self.t_h_BC = X_h_BC[:, 0:1], X_h_BC[:, 1:2]  # 提取水深边界条件 x 和 t 坐标self.x_u_obs, self.t_u_obs = X_u_obs[:, 0:1], X_u_obs[:, 1:2] if useObs else (None, None)  # 提取流速观测数据坐标self.x_h_obs, self.t_h_obs = X_h_obs[:, 0:1], X_h_obs[:, 1:2] if useObs else (None, None)  # 提取水深观测数据坐标self.x_f, self.t_f = X_f[:, 0:1], X_f[:, 1:2]  # 提取内部点 x 和 t 坐标self.h_IC, self.u_BC, self.h_BC = h_IC, u_BC, h_BC  # 存储初始和边界条件值self.u_obs, self.h_obs = u_obs, h_obs if useObs else (None, None)  # 存储观测值self.layers = layers  # 存储神经网络层结构# 初始化权重和偏置if ExistModel == 0:self.weights, self.biases = self.initialize_NN(layers, use_ff)  # 新建网络参数else:self.weights, self.biases = self.load_NN(uhDir, layers)  # 加载已有模型参数# 创建 TensorFlow 会话self.sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True, log_device_placement=True))  # 初始化会话，支持软放置和设备日志self.placeholders = self._define_placeholders()  # 定义输入数据的占位符self._build_model()  # 构建神经网络模型self._define_loss_and_optimizer()  # 定义损失函数和优化器def _define_placeholders(self):"""定义 TensorFlow 占位符用于输入数据。"""placeholders = {'x_h_IC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：初始条件 x 坐标't_h_IC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：初始条件 t 坐标'h_IC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：初始水深值'x_u_BC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：流速边界 x 坐标't_u_BC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：流速边界 t 坐标'u_BC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：边界流速值'x_h_BC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：水深边界 x 坐标't_h_BC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：水深边界 t 坐标'h_BC': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：边界水深值'x_u_obs': tf.placeholder(self.DTYPE, [None, 1]) if self.useObs else None,  # 占位符：流速观测 x 坐标（若使用观测）'t_u_obs': tf.placeholder(self.DTYPE, [None, 1]) if self.useObs else None,  # 占位符：流速观测 t 坐标（若使用观测）'u_obs': tf.placeholder(self.DTYPE, [None, 1]) if self.useObs else None,  # 占位符：观测流速值（若使用观测）'x_h_obs': tf.placeholder(self.DTYPE, [None, 1]) if self.useObs else None,  # 占位符：水深观测 x 坐标（若使用观测）'t_h_obs': tf.placeholder(self.DTYPE, [None, 1]) if self.useObs else None,  # 占位符：水深观测 t 坐标（若使用观测）'h_obs': tf.placeholder(self.DTYPE, [None, 1]) if self.useObs else None,  # 占位符：观测水深值（若使用观测）'x_f': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：内部点 x 坐标't_f': tf.placeholder(self.DTYPE, [None, 1]),  # 占位符：内部点 t 坐标'adaptive_bcs_u': tf.placeholder(tf.float32, [1]) if ExistModel == 0 else None,  # 占位符：流速边界自适应权重（新建模型时使用）'adaptive_bcs_h': tf.placeholder(tf.float32, [1]) if ExistModel == 0 else None,  # 占位符：水深边界自适应权重（新建模型时使用）'adaptive_ics_h': tf.placeholder(tf.float32, [1]) if ExistModel == 0 else None,  # 占位符：初始条件自适应权重（新建模型时使用）'adaptive_obs_u': tf.placeholder(tf.float32, [1]) if ExistModel == 0 and self.useObs else None,  # 占位符：流速观测自适应权重（新建模型且使用观测时使用）'adaptive_obs_h': tf.placeholder(tf.float32, [1]) if ExistModel == 0 and self.useObs else None  # 占位符：水深观测自适应权重（新建模型且使用观测时使用）}return placeholders  # 返回占位符字典def _build_model(self):"""构建神经网络模型，包括预测和残差计算。"""# 归一化输入数据到 [-1, 1]X = 2.0 * (tf.concat([self.placeholders['x_f'], self.placeholders['t_f']], 1) - self.lb) / (self.ub - self.lb) - 1.0self.u_pred, self.h_pred = self.neural_net(X, self.weights, self.biases)  # 预测流速和水深self.u_ic_pred, self.h_ic_pred = self.neural_net(2.0 * (tf.concat([self.placeholders['x_h_IC'], self.placeholders['t_h_IC']], 1) - self.lb) / (self.ub - self.lb) - 1.0,self.weights, self.biases)  # 初始条件预测self.u_bc_pred, self.h_bc_pred = self.neural_net(2.0 * (tf.concat([self.placeholders['x_u_BC'], self.placeholders['t_u_BC']], 1) - self.lb) / (self.ub - self.lb) - 1.0,self.weights, self.biases)  # 边界条件预测if self.useObs:self.u_obs_pred, self.h_obs_pred = self.neural_net(2.0 * (tf.concat([self.placeholders['x_u_obs'], self.placeholders['t_u_obs']], 1) - self.lb) / (self.ub - self.lb) - 1.0,self.weights, self.biases)  # 观测数据预测self.eq1_pred, self.eq2_pred = self._compute_residuals()  # 计算残差def _compute_residuals(self):"""计算 SVE 的连续性和动量残差。"""u_t = tf.gradients(self.u_pred, self.placeholders['t_f'])[0]  # 流速时间导数u_x = tf.gradients(self.u_pred, self.placeholders['x_f'])[0]  # 流速空间导数h_t = tf.gradients(self.h_pred, self.placeholders['t_f'])[0]  # 水深时间导数h_x = tf.gradients(self.h_pred, self.placeholders['x_f'])[0]  # 水深空间导数eq1 = self.fun_r_mass(self.u_pred, h_t, h_x)  # 连续性残差eq2 = self.fun_r_momentum(self.u_pred, self.h_pred, u_t, u_x, h_x)  # 动量残差return eq1, eq2  # 返回残差def fun_r_mass(self, u, h_t, h_x):"""计算连续性方程残差：∂h/∂t + u ∂h/∂x = 0。"""return h_t + u * h_x  # 返回连续性残差def fun_r_momentum(self, u, h, u_t, u_x, h_x):"""计算动量方程残差：∂u/∂t + u ∂u/∂x + g ∂h/∂x + g (n²|u|u / R^(4/3) - S) = 0。"""n = 0.015  # Manning 粗糙度系数h = tf.clip_by_value(h, clip_value_min=1e-4, clip_value_max=50)  # 限制水深范围R = self.b * h / (2 * h + self.b)  # 液压半径return u_t + u * u_x + 9.81 * h_x + 9.81 * (n * n * tf.abs(u) * u / tf.pow(tf.square(R), 2.0 / 3.0) - self.S)  # 返回动量残差def _define_loss_and_optimizer(self):"""定义损失函数和优化器。"""self.loss_f_c = tf.reduce_mean(tf.square(self.eq1_pred))  # 计算连续性损失的均方误差self.loss_f_m = tf.reduce_mean(tf.square(self.eq2_pred))  # 计算动量损失的均方误差self.loss_f = self.loss_f_c + self.loss_f_m  # 总物理损失self.loss_bc_u = tf.reduce_mean(tf.square(self.placeholders['u_BC'] - self.u_bc_pred))  # 流速边界损失self.loss_bc_h = tf.reduce_mean(tf.square(self.placeholders['h_BC'] - self.h_bc_pred))  # 水深边界损失self.loss_bcs = self.loss_bc_u + self.loss_bc_h  # 总边界损失self.loss_ic_h = tf.reduce_mean(tf.square(self.placeholders['h_IC'] - self.h_ic_pred))  # 初始条件损失self.loss = self.loss_f + self.loss_bcs + self.loss_ic_h  # 总损失if self.useObs:self.loss_obs_u = tf.reduce_mean(tf.square(self.placeholders['u_obs'] - self.u_obs_pred))  # 流速观测损失self.loss_obs_h = tf.reduce_mean(tf.square(self.placeholders['h_obs'] - self.h_obs_pred))  # 水深观测损失self.loss += self.loss_obs_u + self.loss_obs_h  # 加入观测损失# 定义 L-BFGS-B 优化器self.optimizer = tf.contrib.opt.ScipyOptimizerInterface(self.loss,method='L-BFGS-B',options={'maxiter': 50000,  # 最大迭代次数'maxfun': 50000,  # 最大函数评估次数'maxcor': 50,  # 最大校正向量数'maxls': 50,  # 最大线搜索步数'ftol': 1.0e-10,  # 函数值容差'gtol': 0.000001})  # 梯度容差self.global_step = tf.Variable(0, trainable=False)  # 全局步数，训练时递增self.learning_rate = tf.train.exponential_decay(lr, self.global_step, 5000, 0.9, staircase=False)  # 学习率衰减self.optimizer_adam = tf.train.AdamOptimizer(learning_rate=self.learning_rate)  # Adam 优化器self.train_op_adam = self.optimizer_adam.minimize(self.loss, global_step=self.global_step)  # Adam 优化操作def train(self, num_epochs, tf_dict):"""训练模型，使用 Adam 优化器。"""for it in range(num_epochs):  # 循环指定次数start_time = time.time()  # 记录训练开始时间self.sess.run(self.train_op_adam, feed_dict=tf_dict)  # 执行 Adam 优化一步if it % 10 == 0:  # 每 10 步打印一次elapsed = time.time() - start_time  # 计算单步耗时loss_value = self.sess.run(self.loss, feed_dict=tf_dict)  # 获取当前损失learning_rate = self.sess.run(self.learning_rate)  # 获取当前学习率print('Iteration: %d, Loss: %.3e, Time: %.2f, Learning Rate: %.3e' % (it, loss_value, elapsed, learning_rate))  # 打印信息u_pred, h_pred = self.predict(self.X_star[:, 0:1], self.X_star[:, 1:2])  # 预测结果error_u = np.linalg.norm(self.u_star - u_pred, 2) / np.linalg.norm(self.u_star, 2)  # 流速 L2 误差error_h = np.linalg.norm(self.h_star - h_pred, 2) / np.linalg.norm(self.h_star, 2)  # 水深 L2 误差print('Error u: %.3e, Error h: %.3e' % (error_u, error_h))  # 打印误差def predict(self, x_star, t_star):"""预测流速和水深。"""tf_dict = {self.placeholders['x_f']: x_star,  # 输入 x 坐标self.placeholders['t_f']: t_star,  # 输入 t 坐标self.placeholders['x_u_tf']: x_star,  # 流速预测 x 坐标self.placeholders['t_u_tf']: t_star,  # 流速预测 t 坐标self.placeholders['x_h_tf']: x_star,  # 水深预测 x 坐标self.placeholders['t_h_tf']: t_star  # 水深预测 t 坐标}u_pred, h_pred = self.sess.run([self.u_pred, self.h_pred], feed_dict=tf_dict)  # 获取预测值return u_pred, h_pred  # 返回流速和水深预测def initialize_NN(self, layers, use_ff):"""初始化神经网络权重和偏置，使用 Xavier 初始化。"""weights, biases = [], []  # 初始化权重和偏置列表num_layers = len(layers)  # 获取层数if use_ff:  # 如果使用傅里叶特征，调整输入层layers = [2 * layers[0]] + layers[1:]for l in range(num_layers - 1):  # 遍历每一层W = self.xavier_init([layers[l], layers[l + 1]])  # 初始化权重b = tf.Variable(tf.zeros([1, layers[l + 1]], dtype=self.DTYPE), dtype=self.DTYPE)  # 初始化偏置weights.append(W)  # 添加权重biases.append(b)  # 添加偏置return weights, biases  # 返回权重和偏置def xavier_init(self, size):"""Xavier 初始化权重。"""in_dim, out_dim = size  # 获取输入和输出维度xavier_stddev = np.sqrt(2 / (in_dim + out_dim))  # 计算 Xavier 标准差return tf.Variable(tf.truncated_normal(size, stddev=xavier_stddev), dtype=self.DTYPE)  # 返回初始化权重def neural_net(self, X, weights, biases):"""定义前馈神经网络。"""H = X  # 输入数据for l in range(len(weights) - 1):  # 遍历隐藏层W, b = weights[l], biases[l]  # 获取当前层的权重和偏置H = tf.tanh(tf.add(tf.matmul(H, W), b))  # 应用激活函数 tanhW, b = weights[-1], biases[-1]  # 获取输出层的权重和偏置Y = tf.add(tf.matmul(H, W), b)  # 计算输出return Y[:, 0:1], Y[:, 1:2]  # 返回流速 u 和水深 hdef add_noise(self, signal):"""添加高斯噪声到信号。"""target_snr_db = 500  # 目标信噪比 (dB)sig_avg_db = 10 * np.log10(np.mean(signal))  # 信号平均功率 (dB)noise_avg_db = sig_avg_db - target_snr_db  # 噪声平均功率 (dB)noise_avg = 10 ** (noise_avg_db / 10)  # 噪声平均功率 (W)noise = np.random.normal(0, np.sqrt(noise_avg), len(signal))  # 生成白噪声return signal + noise  # 返回含噪声的信号# 主程序
if __name__ == "__main__":# 定义英尺到米的转换因子FACTOR = 0.3048# 定义每个 Case 的配置字典CASES = {1: {'file': 'case1/MixedFlow.p02.hdf', 'dx': 30, 'dt': 30, 'n': 0.005, 'u': 1.0},2: {'file': 'case2/MixedFlow.p02.hdf', 'dx': 20, 'dt': 30, 'n': 0.025, 'u': 1.0, 'S': 1e-3},3: {'file': 'case3/MixedFlow.p02.hdf', 'obs': [16, 32]},4: {'file': 'case4/MixedFlow.p02.hdf', 'obs': [16, 32], 'use_ff': True},5: {'file': 'case5/MixedFlow.p02.hdf', 'obs': [118, 134]},6: {'file': 'case6/MixedFlow.p02.hdf', 'obs': [[118, 134], [25, 118, 134], [25, 50, 118, 134], [25, 50, 100, 118, 134]]}}# 定义默认网络层结构LAYERS = [2] + 4 * [1 * 64] + [2]for case_id in range(1, 7):  # 循环运行每个 Caseprint(f"\nRunning Case {case_id}...")  # 打印当前 Case 编号case = CASES[case_id]  # 获取当前 Case 配置hdf_file = f'HEC-RAS/{case["file"]}' if case_id > 2 else None  # 根据 Case 类型设置 HDF 文件路径use_ff = case.get('use_ff', False)  # 获取是否使用傅里叶特征obs_indices = case.get('obs', [[int(1200 / case.get('dx', 30)), int(2400 / case.get('dx', 30))]])  # 获取观测点索引# 数据准备if hdf_file:  # 如果是 Case 3-6，使用 HDF5 数据hf = h5py.File(hdf_file, 'r')  # 打开 HDF5 文件attrs = hf['Geometry']['Cross Sections']['Attributes'][:]  # 读取几何属性staid, eles, reach_len = [], [], []  # 初始化列表for attr in attrs:  # 遍历属性staid.append(attr[2].decode('utf-8'))  # 提取站点 IDeles.append(attr[14])  # 提取高程reach_len.append(attr[6])  # 提取河段长度coor = np.cumsum(np.array(reach_len[:-1]))  # 计算累积坐标coor = [0] + coor.tolist()  # 添加起始点coor = coor[::-1]  # 反转坐标eles = np.array(eles)  # 转换为数组slope = np.gradient(eles, coor)  # 计算坡度water_surface = hf['Results']['Unsteady']['Output']["/Results/Unsteady/Output"]['Output Blocks']['Base Output']["Unsteady Time Series"]["Cross Sections"]['Water Surface'][1:]  # 读取水面高程velocity_total = hf['Results']['Unsteady']['Output']["/Results/Unsteady/Output"]['Output Blocks']['Base Output']["Unsteady Time Series"]["Cross Sections"]['Velocity Total'][1:]  # 读取总流速timestamp = hf['Results']['Unsteady']['Output']["/Results/Unsteady/Output"]['Output Blocks']['Base Output']["Unsteady Time Series"]['Time Date Stamp'][1:]  # 读取时间戳time_model = [datetime.strptime(t.decode('utf-8'), '%d%b%Y %H:%M:%S') for t in timestamp]  # 转换时间格式water_depth = water_surface - eles[None, :]  # 计算水深Nt, Nx = water_depth.shape  # 获取时间和空间维度t = np.arange(Nt)[:, None]  # 生成时间数组x = np.array(coor[::-1])[:, None]  # 生成空间数组u_exact, h_exact = velocity_total, water_depth  # 提取真实流速和水深else:  # Case 1-2，使用模拟数据dx, dt = case['dx'], case['dt']  # 获取空间和时间步长n, u, S = case['n'], case['u'], case.get('S', 0)  # 获取 Manning 系数、流速和坡度t = np.arange(0, 3600 + dt, dt)  # 生成时间序列x = np.arange(0, 3600 + dx, dx)  # 生成空间序列Nt, Nx = len(t), len(x)  # 获取维度if case_id == 1:  # Case 1：解析解h_exact = np.zeros((Nt, Nx))  # 初始化水深数组for i in range(Nt):  # 遍历时间xb = u * t[i]  # 计算边界位置indb = np.argwhere(np.abs(x - xb) == np.abs(x - xb).min())[0][0]  # 找到最近点if x[indb] > xb:  # 调整索引indb -= 1h_exact[i, :indb + 1] = np.power(-7 / 3 * n * n * u * u * (x[:indb + 1] - u * t[i]), 3 / 7)  # 计算解析解elif case_id == 2:  # Case 2：RK4 解h_exact = np.zeros((Nt, Nx))  # 初始化水深数组h_exact[0, 0] = 4 * np.sin(2 * np.pi / (4 * 3600) * t[0])  # 初始边界条件for i in range(1, Nt):  # 遍历时间xb = u * t[i]  # 计算边界位置indb = np.argwhere(np.abs(x - xb) == np.abs(x - xb).min())[0][0]  # 找到最近点if x[indb] > xb:  # 调整索引indb -= 1x_ = x[:indb + 1]  # 截取空间范围y0 = 4 * np.sin(2 * np.pi / (4 * 3600) * t[i])  # 边界水深h_ = [y0]  # 初始值for j in range(1, len(x_)):  # RK4 迭代k1 = dx * (-S - n * n * u * u / np.power(np.power(h_[-1], 4), 1 / 3.))k2 = dx * (-S - n * n * u * u / np.power(np.power(h_[-1] + 0.5 * k1, 4), 1 / 3.))k3 = dx * (-S - n * n * u * u / np.power(np.power(h_[-1] + 0.5 * k2, 4), 1 / 3.))k4 = dx * (-S - n * n * u * u / np.power(np.power(h_[-1] + k3, 4), 1 / 3.))y = h_[-1] + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)  # 更新水深if y >= 0:  # 确保水深非负h_.append(y)else:breakh_exact[i, :len(h_)] = h_  # 存储 RK4 结果u_exact = np.zeros_like(h_exact) + u  # 假设恒定流速X, T = np.meshgrid(x, t)  # 生成网格X_star = np.hstack((X.flatten()[:, None], T.flatten()[:, None]))  # 展平坐标u_star, h_star = u_exact.flatten()[:, None], h_exact.flatten()[:, None]  # 展平目标值lb, ub = X_star.min(0), X_star.max(0)  # 获取域边界# 准备训练数据tsteps = [0, Nt - 1]  # 初始和结束时间步X_h_IC, h_IC = None, None  # 初始化初始条件for i, tstep in enumerate(tsteps):  # 遍历时间步xx1_ = np.hstack((X[tstep:tstep + 1, :].T, T[tstep:tstep + 1, :].T))  # 提取坐标hh1_ = self.add_noise(h_exact[tstep:tstep + 1, :].T)  # 添加噪声if i == 0:  # 首次赋值X_h_IC, h_IC = xx1_, hh1_else:  # 追加数据X_h_IC, h_IC = np.vstack((X_h_IC, xx1_)), np.vstack((h_IC, hh1_))X_u_BC = np.vstack((np.hstack((X[:, 0:1], T[:, 0:1])), np.hstack((u * t, t))))  # 上下游边界坐标X_h_BC = X_u_BC  # 水深边界坐标u_BC = np.vstack((u_exact[:, 0:1], h_exact[:, np.argmin(np.abs(x - u * t), axis=1)]))  # 流速边界值h_BC = np.vstack((h_exact[:, 0:1], h_exact[:, np.argmin(np.abs(x - u * t), axis=1)]))  # 水深边界值X_u_obs, X_h_obs, u_obs, h_obs = None, None, None, None  # 初始化观测数据if case_id > 2 or useObs:  # 如果使用观测数据t_obs_u, x_obs_u, u_obs = None, None, None  # 初始化流速观测t_obs_h, x_obs_h, h_obs = None, None, None  # 初始化水深观测for obs_idx in obs_indices[0] if case_id == 6 else obs_indices:  # 遍历观测点t_obs_u = np.append(t_obs_u, t.flatten()) if 't_obs_u' in locals() else t.flatten()  # 追加时间x_obs_u = np.append(x_obs_u, np.ones(Nt) * x[obs_idx]) if 'x_obs_u' in locals() else np.ones(Nt) * x[obs_idx]  # 追加 x 坐标u_obs = np.append(u_obs, self.add_noise(u_exact[:, obs_idx])) if 'u_obs' in locals() else self.add_noise(u_exact[:, obs_idx])  # 追加流速t_obs_h = np.append(t_obs_h, t.flatten()) if 't_obs_h' in locals() else t.flatten()  # 追加时间x_obs_h = np.append(x_obs_h, np.ones(Nt) * x[obs_idx]) if 'x_obs_h' in locals() else np.ones(Nt) * x[obs_idx]  # 追加 x 坐标h_obs = np.append(h_obs, self.add_noise(h_exact[:, obs_idx])) if 'h_obs' in locals() else self.add_noise(h_exact[:, obs_idx])  # 追加水深X_u_obs, X_h_obs = np.vstack((x_obs_u, t_obs_u)).T, np.vstack((x_obs_h, t_obs_h)).T  # 组合坐标u_obs, h_obs = u_obs[:, None], h_obs[:, None]  # 转换为列向量X_f_train = X_star  # 内部训练点slope = np.hstack([np.array(slope) for _ in range(Nt)])[:, None] if hdf_file else 0  # 扩展坡度数组# 模型训练model = SVE(X_h_IC, X_u_BC, X_h_BC, X_u_obs, X_h_obs, X_f_train, h_IC, u_BC, h_BC, u_obs, h_obs,LAYERS if case_id > 2 else [2] + 3 * [1 * 32] + [1], lb, ub, slope, 10,X_star, u_star, h_star, useObs=True, use_ff=use_ff)  # 创建模型实例tf_dict = {k: v for k, v in model.placeholders.items() if v is not None}  # 创建字典，仅包含非空占位符tf_dict.update({  # 更新字典model.placeholders['x_h_IC']: X_h_IC,model.placeholders['t_h_IC']: X_h_IC[:, 1:2],model.placeholders['h_IC']: h_IC,model.placeholders['x_u_BC']: X_u_BC[:, 0:1],model.placeholders['t_u_BC']: X_u_BC[:, 1:2],model.placeholders['u_BC']: u_BC,model.placeholders['x_h_BC']: X_h_BC[:, 0:1],model.placeholders['t_h_BC']: X_h_BC[:, 1:2],model.placeholders['h_BC']: h_BC})if useObs:  # 如果使用观测数据tf_dict.update({model.placeholders['x_u_obs']: X_u_obs[:, 0:1],model.placeholders['t_u_obs']: X_u_obs[:, 1:2],model.placeholders['u_obs']: u_obs,model.placeholders['x_h_obs']: X_h_obs[:, 0:1],model.placeholders['t_h_obs']: X_h_obs[:, 1:2],model.placeholders['h_obs']: h_obs})tf_dict.update({model.placeholders['x_f']: X_f_train[:, 0:1],  # 内部点 xmodel.placeholders['t_f']: X_f_train[:, 1:2]})  # 内部点 tmodel.train(10000, tf_dict)  # 训练模型 10000 轮# 预测和评估x_test, t_test = X_star[:Nt * Nx, 0:1], X_star[:Nt * Nx, 1:2]  # 准备测试数据u_pred, h_pred = model.predict(x_test, t_test)  # 预测结果u_pred, h_pred = u_pred.reshape(Nt, Nx), h_pred.reshape(Nt, Nx)  # 重塑为原始形状u_exact, h_exact = u_exact[:Nt], h_exact[:Nt]  # 截取测试范围error_h = np.linalg.norm(h_exact - h_pred, 2) / np.linalg.norm(h_exact, 2)  # 计算 L2 误差rmse_h = np.sqrt(((h_exact - h_pred) ** 2).mean())  # 计算 RMSEprint(f'Case {case_id} - Error h: {error_h:.3e}, RMSE h: {rmse_h:.3f} m')  # 打印结果# 可视化hours = np.arange(Nt) / 120. if case_id > 2 else np.arange(Nt) * dt / 3600  # 转换为小时x = x[::-1] * FACTOR if case_id > 2 else x  # 反转并转换单位xx, tt = np.meshgrid(x, hours)  # 生成网格h_pred *= FACTOR if case_id > 2 else 1  # 转换单位h_exact *= FACTOR if case_id > 2 else 1  # 转换单位eles = np.zeros_like(x) if case_id <= 2 else eles * FACTOR  # 地形高程fig = plt.figure(figsize=(16.5, 15 if case_id == 6 else 8))  # 创建图形对象gs = gridspec.GridSpec(3 if case_id <= 4 else 5, 4 if case_id == 6 else 3, hspace=0.08, wspace=0.15)  # 设置网格levels = np.linspace(0, 5.4 if case_id > 2 else 0.6, 9 if case_id > 2 else 10)  # 设置颜色条范围ax0 = fig.add_subplot(gs[0, 1:3 if case_id <= 4 else 3:5])  # 添加第一个子图cs = ax0.contourf(xx, tt, h_exact[:Nt], cmap='rainbow', levels=levels, alpha=0.8)  # 绘制参考水深图ax0.scatter(X_u_BC[:, 0] * FACTOR if case_id > 2 else X[:, 0:1][::3],  # 绘制边界条件点X_u_BC[:, 1] / 120. if case_id > 2 else T[:, 0:1][::3],marker='o', c='g', s=12, clip_on=False)ax0.scatter(x_obs_h * FACTOR if case_id > 2 else x_obs[::2],  # 绘制观测点t_obs_h / 120. if case_id > 2 else t_obs[::2],facecolors='none', edgecolors='k', marker='o', s=15, clip_on=False)ax0.scatter(X_h_IC[:, 0] * FACTOR if case_id > 2 else xx1_[:, 0][::3],  # 绘制快照点X_h_IC[:, 1] / 120. if case_id > 2 else xx1_[:, 1][::3],marker='*', c='r', s=25, clip_on=False)ax0.set_ylabel('Time (h)')  # 设置 y 轴标签ax0.set_xticklabels([] if case_id <= 4 else ax0.get_xticklabels())  # 隐藏 x 轴标签（若非最后一行）ax0.text(0.05, 0.9, f'(a) Ref', fontsize=16, transform=ax0.transAxes)  # 添加子图标签divider = make_axes_locatable(ax0)  # 创建颜色条布局cax = divider.append_axes("right", size="2%", pad=0.05)  # 添加颜色条轴cb = fig.colorbar(cs, cax=cax, orientation='vertical')  # 添加颜色条cb.ax.tick_params(labelsize=14)  # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14)  # 设置偏移文本字体大小cb.set_label('Water depth (m)', fontsize=14)  # 设置颜色条标签if case_id == 6:  # Case 6 特殊处理，4 列 PINN 结果ax1, ax2, ax3, ax4 = [fig.add_subplot(gs[1, i * 2:(i + 1) * 2]) for i in range(4)]  # 添加 4 个子图for ax, h_pred_, label in zip([ax1, ax2, ax3, ax4], [h_pred2obs, h_pred3obs, h_pred4obs, h_pred5obs],['(b) PINN', '(c) PINN', '(d) PINN', '(e) PINN']):cs = ax.contourf(xx, tt, h_pred_[:Nt], cmap='rainbow', levels=levels, alpha=0.8)  # 绘制 PINN 水深ax.set_xticklabels([])  # 隐藏 x 轴标签ax.set_yticklabels([] if ax != ax1 else ax.get_yticklabels())  # 仅第一个子图显示 y 轴标签ax.text(0.05, 0.9, label, fontsize=16, transform=ax.transAxes)  # 添加子图标签divider = make_axes_locatable(ax4)  # 创建颜色条布局cax = divider.append_axes("right", size="2%", pad=0.05)  # 添加颜色条轴cb = fig.colorbar(cs, cax=cax, orientation='vertical')  # 添加颜色条cb.ax.tick_params(labelsize=14)  # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14)  # 设置偏移文本字体大小cb.set_label('Water depth (m)', fontsize=14)  # 设置颜色条标签ax5, ax6, ax7, ax8 = [fig.add_subplot(gs[2, i * 2:(i + 1) * 2]) for i in range(4)]  # 添加 4 个误差子图for ax, h_pred_, label in zip([ax5, ax6, ax7, ax8], [h_pred2obs, h_pred3obs, h_pred4obs, h_pred5obs],['(f) PINN-Ref', '(g) PINN-Ref', '(h) PINN-Ref', '(i) PINN-Ref']):cs = ax.contourf(xx, tt, h_pred_[:Nt] - h_exact[:Nt], cmap='bwr', levels=np.linspace(-0.5, 0.5, 11),alpha=0.8, extend='both')  # 绘制 PINN 误差ax.set_xticklabels([])  # 隐藏 x 轴标签ax.set_yticklabels([] if ax != ax5 else ax.get_yticklabels())  # 仅第一个子图显示 y 轴标签ax.text(0.05, 0.9, label, fontsize=16, transform=ax.transAxes)  # 添加子图标签ax9, ax10, ax11, ax12 = [fig.add_subplot(gs[3, i * 2:(i + 1) * 2]) for i in range(4)]  # 添加 4 个插值子图for ax, h_interp_, label in zip([ax9, ax10, ax11, ax12], [h_interp2obs, h_interp3obs, h_interp4obs, h_interp5obs],['(j) Interpolation', '(k) Interpolation', '(l) Interpolation', '(m) Interpolation']):cs = ax.contourf(xx, tt, h_interp_[:Nt], cmap='rainbow', levels=levels, alpha=0.8)  # 绘制插值水深ax.set_xticklabels([])  # 隐藏 x 轴标签ax.set_yticklabels([] if ax != ax9 else ax.get_yticklabels())  # 仅第一个子图显示 y 轴标签ax.text(0.05, 0.9, label, fontsize=16, transform=ax.transAxes)  # 添加子图标签divider = make_axes_locatable(ax12)  # 创建颜色条布局cax = divider.append_axes("right", size="2%", pad=0.05)  # 添加颜色条轴cb = fig.colorbar(cs, cax=cax, orientation='vertical')  # 添加颜色条cb.ax.tick_params(labelsize=14)  # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14)  # 设置偏移文本字体大小cb.set_label('Water depth (m)', fontsize=14)  # 设置颜色条标签ax13, ax14, ax15, ax16 = [fig.add_subplot(gs[4, i * 2:(i + 1) * 2]) for i in range(4)]  # 添加 4 个插值误差子图for ax, h_interp_, label in zip([ax13, ax14, ax15, ax16], [h_interp2obs, h_interp3obs, h_interp4obs, h_interp5obs],['(n) Interpolation-Ref', '(o) Interpolation-Ref', '(p) Interpolation-Ref', '(q) Interpolation-Ref']):cs = ax.contourf(xx, tt, h_interp_[:Nt] - h_exact[:Nt], cmap='bwr', levels=np.linspace(-0.5, 0.5, 11),alpha=0.8, extend='both')  # 绘制插值误差ax.set_xlabel('Distance upstream (m)' if case_id > 2 else 'Distance (m)')  # 设置 x 轴标签ax.set_yticklabels([] if ax != ax13 else ax.get_yticklabels())  # 仅第一个子图显示 y 轴标签ax.text(0.05, 0.9, label, fontsize=16, transform=ax.transAxes)  # 添加子图标签divider = make_axes_locatable(ax16)  # 创建颜色条布局cax = divider.append_axes("right", size="2%", pad=0.05)  # 添加颜色条轴cb = fig.colorbar(cs, cax=cax, orientation='vertical')  # 添加颜色条cb.ax.tick_params(labelsize=14)  # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14)  # 设置偏移文本字体大小cb.set_label('Error (m)', fontsize=14)  # 设置颜色条标签else:  # Case 1-4 或 5 的标准布局ax1 = fig.add_subplot(gs[1, :2])  # 添加 PINN 子图cs = ax1.contourf(xx, tt, h_pred[:Nt], cmap='rainbow', levels=levels, alpha=0.8)  # 绘制 PINN 水深ax1.set_ylabel('Time (h)' if case_id > 2 else 'Time (s)')  # 设置 y 轴标签ax1.set_xticklabels([])  # 隐藏 x 轴标签ax1.text(0.05, 0.9, '(b) PINN', fontsize=16, transform=ax1.transAxes)  # 添加子图标签ax2 = fig.add_subplot(gs[1, 2:])  # 添加参考或插值子图cs = ax2.contourf(xx, tt, h_exact[:Nt], cmap='rainbow', levels=levels, alpha=0.8)  # 绘制参考或插值水深ax2.set_xticklabels([])  # 隐藏 x 轴标签ax2.set_yticklabels([] if case_id <= 4 else ax2.get_yticklabels())  # 隐藏 y 轴标签（若非 Case 5-6）ax2.text(0.05, 0.9, '(c) Interpolation' if case_id > 2 else '(c) Ref', fontsize=16, transform=ax2.transAxes)  # 添加子图标签divider = make_axes_locatable(ax2)  # 创建颜色条布局cax = divider.append_axes("right", size="2%", pad=0.05)  # 添加颜色条轴cb = fig.colorbar(cs, cax=cax, orientation='vertical')  # 添加颜色条cb.ax.tick_params(labelsize=14)  # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14)  # 设置偏移文本字体大小cb.set_label('Water depth (m)', fontsize=14)  # 设置颜色条标签ax3 = fig.add_subplot(gs[2, :2])  # 添加 PINN 误差子图cs = ax3.contourf(xx, tt, h_pred[:Nt] - h_exact[:Nt], cmap='bwr', levels=np.linspace(-0.5, 0.5, 11),alpha=0.8, extend='both')  # 绘制 PINN 误差ax3.set_xlabel('Distance upstream (m)' if case_id > 2 else 'Distance (m)')  # 设置 x 轴标签ax3.set_ylabel('Time (h)' if case_id > 2 else 'Time (s)')  # 设置 y 轴标签ax3.text(0.05, 0.9, '(d) PINN-Ref', fontsize=16, transform=ax3.transAxes)  # 添加子图标签ax4 = fig.add_subplot(gs[2, 2:])  # 添加参考或插值误差子图cs = ax4.contourf(xx, tt, h_exact[:Nt] - h_exact[:Nt] if case_id <= 2 else h_exact[:Nt],cmap='bwr', levels=np.linspace(-0.5, 0.5, 11), alpha=0.8, extend='both')  # 绘制误差ax4.set_xlabel('Distance upstream (m)' if case_id > 2 else 'Distance (m)')  # 设置 x 轴标签ax4.set_yticklabels([] if case_id <= 4 else ax4.get_yticklabels())  # 隐藏 y 轴标签（若非 Case 5-6）ax4.text(0.05, 0.9, '(e) Interpolation-Ref' if case_id > 2 else '(e) Ref-Ref', fontsize=16, transform=ax4.transAxes)  # 添加子图标签divider = make_axes_locatable(ax4)  # 创建颜色条布局cax = divider.append_axes("right", size="2%", pad=0.05)  # 添加颜色条轴cb = fig.colorbar(cs, cax=cax, orientation='vertical')  # 添加颜色条cb.ax.tick_params(labelsize=14)  # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14)  # 设置偏移文本字体大小cb.set_label('Error (m)', fontsize=14)  # 设置颜色条标签tlist = [20, 50, 100, 200] if case_id > 2 else [int(i * 60 / dt) for i in [20, 30, 40, 60]]  # 选择时间点fig, axes = plt.subplots(2, 2, figsize=(15, 6 if case_id <= 4 else 12), sharex=True, sharey=True)  # 创建沿河图axes = axes.ravel()  # 展平轴数组for k in range(4):  # 遍历四个时间点axes[k].plot(x, h_exact[tlist[k]] + eles if case_id > 2 else h_exact[tlist[k]], 'ok', label='reference')  # 绘制参考曲线axes[k].plot(x, h_pred[tlist[k]] + eles if case_id > 2 else h_pred[tlist[k]], '-r', linewidth=2, label='PINN')  # 绘制 PINN 曲线if case_id > 2:  # Case 3-6 添加插值曲线axes[k].plot(x, h_exact[tlist[k]] + eles, '-g', linewidth=2, label='Interpolation')  # 绘制插值曲线axes[k].fill_between(x.flatten(), eles, color='0.7')  # 填充地形axes[k].text(0.85, 0.85, f't={tlist[k]} h' if case_id > 2 else f't={tlist[k] * dt} s',fontsize=15, transform=axes[k].transAxes)  # 添加时间标签if k in [0, 2]:  # 设置 y 轴标签axes[k].set_ylabel('Water stage (m)')axes[k].set_ylim([0, 5 if case_id > 2 else 0.6])  # 设置 y 轴范围axes[k].grid()  # 添加网格if k in [2, 3]:  # 设置 x 轴标签axes[k].set_xlabel('Distance upstream (m)' if case_id > 2 else 'Distance (m)')if k == 0:  # 添加图例axes[k].legend(loc=2 if case_id > 2 else 4, prop={'size': 14})plt.tight_layout()  # 调整布局plt.savefig(f'figures/case{case_id}/along_channel.pdf')  # 保存沿河图plt.close()  # 关闭图形fig.savefig(f'figures/case{case_id}/contour.pdf')  # 保存轮廓图plt.close()  # 关闭图形

结论

论文作者首次尝试将PINN应用于大尺度河流模型降尺度至子网格中，该框架可解一维SVE方程、通过傅里叶变换编码潮汐边界周期性变化、同化多种数据类型并集成到大型河流模型。实验证明PINN在洪水建模中有重要作用，为解决洪灾问题应用提供了方法。

不足以及展望

作者在最后提出了论文中的不足：

1. 测试案例相对简单，仅应用恒定摩擦系数，未考虑河道几何形状影响。
2. 现实中常见的数据问题如观测不确定性、数据不一致和观测缺失等也未充分探讨，且PINN在解决河口非线性相互作用时精度有待提高。

作者对未来工作的展望如下：
未来可将PINN降尺度应用于更广泛、更现实的流动条件。考虑可变摩擦系数、河道几何形状影响，开发更强大的抗噪声算法处理数据不确定性。探索扩展PINN到二维领域，用于创建复杂流动区域的数字孪生，以解决当前大尺度模型无法解决的问题。