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

2025.3.2周报

一、文献阅读
- 题目信息
- 摘要
- Abstract
- 创新点
- 网络架构
- 实验
- 结论
- 不足以及展望

一、文献阅读

题目信息

题目： 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.

创新点

作者通过基于PINN的数据同化降尺度方法，解决大尺度河流模型在海岸地区子网格尺度河流动力学的变异性问题。该方法能融合多种观测数据，无需修改数值算法或细化网格分辨率，且不受网格限制，为降尺度研究提供新途径。
```
 降尺度是指从一个粗糙的、低分辨率的数据或模型结果，推导出更精细的、高分辨率的结果的过程。
 简单来说就是把模糊的大图变成清晰的小图。
```
作者提出了傅里叶特征嵌入的方法，在输入层前将坐标𝑥和𝑡映射到高维傅里叶特征空间。论文特别针对时间坐标𝑡，以适应潮汐的周期性。下面对其进行说明：
傅里叶变换的核心思想是： 任何周期性的信号都可以分解成一堆不同频率的正弦波和余弦波。因为正弦和余弦是周期函数，可以描述有规律的波动。
其数学表达如下：
$\gamma(\nu)=\left[\begin{array}{c}\cos (2 \pi B \nu) \\\sin (2 \pi B \nu)\end{array}\right], \quad \nu=[x, t]^{T}, \quad 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 \boldsymbol{B}_{t}^{(i)} t\right) \\\sin \left(2 \pi \boldsymbol{B}_{t}^{(i)} t\right)\end{array}\right], \quad \text { for } i=1,2, \ldots\\\end{array}$

网络架构

如图(a)所示，在本片论文中PINN框架由一个深度神经网络与SVE结合的架构。
在输出结果在作者还强调数据同化过程（所谓数据同化，就是让模型预测值不断接近真实值的一个过程，而不是某一个步骤。比如，在算Loss中就算是一个数据同化的步骤）。
下面我们来具体分析一下这个模型的架构：
输入：空间坐标 𝑥 和时间 𝑡 。
结构： 在图示中NN的架构为2个隐藏层和每层4个神经元，其中隐藏层使用激活函数 𝜎
其中𝜎为： $\tanh =\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}$ 。
输出：预测值 𝑢 流速和ℎ水深，由橙色节点表示。此外，输入数据来源于数值解、现场测站数据、卫星提供的水深快照，这些数据通过损失函数约束神经网络，确保预测接近真实值。
物理约束：
此处为SVE公式
然后通过自动微分计算偏导数并代入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_{BC_h}$ 、 $J_{BC_u}$ ：
边界条件损失，分别衡量流速 𝑢和水深 ℎ在边界上的预测值与已知边界条件𝑢和ℎ的偏差。值越小，说明边界符合要求。
其中：
$J_{BC_u}(\theta) = \frac{1}{N_{BC_u}} \sum_{i=1}^{N_{BC_u}} \left| \hat{u}(x^i, t^i, \theta) - u(x^i, t^i) \right|^2$
$J_{BC_h}(\theta) = \frac{1}{N_{BC_h}} \sum_{i=1}^{N_{BC_h}} \left| \hat{h}(x^i, t^i, \theta) - h(x^i, t^i) \right|^2, \quad x \in \Omega_{BC}, t \in [0, T]$
$N_{BC_u}$ 与 $N_{BC_h}$ 是边界条件点的总数； $\sum_{i=1}^{N_{BC_h}}$ 与 $\sum_{i=1}^{N_{BC_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, \quad x \in \Omega_S$
其中 $N_{S_u}$ 与 $N_{S_h}$ 是初始条件或快照点的总数； $\quad x \in \Omega_S$ 定义了河道距离 $x$ 的范围。

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

观测数据损失 $J_{obs_u}$ 与 $J_{obs_h}$
$J_{obs_u}(\theta) = \frac{1}{N_{obs_u}} \sum_{i=1}^{N_{obs_u}} \left| \hat{u}(x^i, t^i, \theta) - u(x^i, t^i) \right|^2$
$J_{obs_h}(\theta) = \frac{1}{N_{obs_h}} \sum_{i=1}^{N_{obs_h}} \left| \hat{h}(x^i, t^i, \theta) - h(x^i, t^i) \right|^2, \quad x \in \Omega_{obs}, t \in [0, T]$
其中 $N_{obs_u}$ 与 $N_{obs_h}$ 观测点数据总数；在 $\quad 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和2中，PINN在解决简化的圣维南方程（SVE）时展现出较高准确性，能很好地模拟洪水在洪泛平原的传播，其洪水波传播模式和淹没前沿形状与解析解或RK4解接近，相对L2误差和RMSE较小。这表明在有足够观测数据和地形信息时，PINN可用于计算亚网格尺度的洪泛平原淹没情况，提供更详细的淹没地图。
在实验3和4中，PINN解与HEC - RAS参考解在水深方面吻合良好，仅在沿河道剖面的上游端附近有较小模拟偏差，L2误差和RMSE也表明其性能良好。在案例4中，编码傅里叶特征的PINN解能够捕捉到下游潮汐引起的周期性变化，相比标准架构的PINN解，其ϵh和RMSE更小。但由于潮汐和边界条件产生的高振荡行为，案例4的预测准确性比案例3差，意味着PINN训练需要更多物理约束和观测数据。
在实验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 数据类型为 float32

    def __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.0
        self.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))  # 应用激活函数 tanh
        W, b = weights[-1], biases[-1]  # 获取输出层的权重和偏置
        Y = tf.add(tf.matmul(H, W), b)  # 计算输出
        return Y[:, 0:1], Y[:, 1:2]  # 返回流速 u 和水深 h

    def 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):  # 循环运行每个 Case
        print(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'))  # 提取站点 ID
                eles.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 -= 1
                    h_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 -= 1
                    x_ = 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:
                            break
                    h_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],  # 内部点 x
                        model.placeholders['t_f']: X_f_train[:, 1:2]})  # 内部点 t
        model.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())  # 计算 RMSE
        print(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()  # 关闭图形