潍坊 网站推广,个人网站赚广告费,手机自己做网站,网站开发电话2025.3.2周报 一、文献阅读题目信息摘要Abstract创新点网络架构实验结论不足以及展望 一、文献阅读
题目信息
题目#xff1a; Physics-Informed Neural Networks of the Saint-Venant Equations for Downscaling a Large-Scale River Model期刊#xff1a; Water Resource… 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的数据同化降尺度方法解决大尺度河流模型在海岸地区子网格尺度河流动力学的变异性问题。该方法能融合多种观测数据无需修改数值算法或细化网格分辨率且不受网格限制为降尺度研究提供新途径。 降尺度是指从一个粗糙的、低分辨率的数据或模型结果推导出更精细的、高分辨率的结果的过程。简单来说就是把模糊的大图变成清晰的小图。作者提出了傅里叶特征嵌入的方法在输入层前将坐标和映射到高维傅里叶特征空间。论文特别针对时间坐标以适应潮汐的周期性。下面对其进行说明 傅里叶变换的核心思想是 任何周期性的信号都可以分解成一堆不同频率的正弦波和余弦波。因为正弦和余弦是周期函数可以描述有规律的波动。 其数学表达如下 γ ( ν ) [ cos ( 2 π B ν ) sin ( 2 π B ν ) ] , ν [ x , t ] T , B ∼ N ( 0 , s ) \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) γ(ν)[cos(2πBν)sin(2πBν)],ν[x,t]T,B∼N(0,s) 其中 ν [ x , t ] T \nu[x, t]^{T} ν[x,t]T输入是空间坐标 和时间 这里是一个二维向量 B是一个随机矩阵里面的值从高斯分布 (0, )中抽取 是标准差 cos ( 2 π B ν ) \cos (2 \pi B \nu) cos(2πBν)与 sin ( 2 π B ν ) \sin (2 \pi B \nu) sin(2πBν)则是对输入进行正弦和余弦变换。 论文中作者通过这个特性将傅里叶特征嵌入到时间中将来描述潮汐因为潮汐具有周期性,具体表示为 γ t ( i ) ( t ) [ cos ( 2 π B t ( i ) t ) sin ( 2 π B t ( i ) t ) ] , for i 1 , 2 , … \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 } i1,2, \ldots\\\end{array} γt(i)(t) cos(2πBt(i)t)sin(2πBt(i)t) , for i1,2,…
网络架构
如图(a)所示在本片论文中PINN框架由一个深度神经网络与SVE结合的架构。 在输出结果在作者还强调数据同化过程所谓数据同化就是让模型预测值不断接近真实值的一个过程而不是某一个步骤。比如在算Loss中就算是一个数据同化的步骤。 下面我们来具体分析一下这个模型的架构 输入空间坐标 和时间 。 结构 在图示中NN的架构为2个隐藏层和每层4个神经元其中隐藏层使用激活函数 其中为 tanh e x − e − x e x e − x \tanh \frac{e^{x}-e^{-x}}{e^{x}e^{-x}} tanhexe−xex−e−x。 输出预测值 流速和ℎ水深由橙色节点表示。此外输入数据来源于数值解、现场测站数据、卫星提供的水深快照这些数据通过损失函数约束神经网络确保预测接近真实值。 物理约束 然后通过自动微分计算偏导数并代入SVE检查残差。
损失函数 损失函数定义为 J ( θ ) J f ( θ ) ∑ j w j J j ( θ ) J(\theta)J_{f}(\theta)\sum_{j} w_{j} J_{j}(\theta) J(θ)Jf(θ)j∑wjJj(θ) 1、 J f ( θ ) J_{f}(\theta) Jf(θ)SVE残差损失。其中 J f ( θ ) 1 N f ∑ i 1 N f ∣ r f 2 ( x i , t i , θ ) ∣ J_f(\theta) \frac{1}{N_f} \sum_{i1}^{N_f} \left| r_f^2(x^i, t^i, \theta) \right| Jf(θ)Nf1∑i1Nf rf2(xi,ti,θ) N f N_f Nf为配置点 r f r_f rf是SVE的残差表示 θ \theta θ是神经网络的参数 x i , t i x^i, t^i xi,ti是第 个配置点的空间坐标和时间坐标。
2、 ∑ j w j J j ( θ ) \sum_{j} w_{j} J_{j}(\theta) ∑jwjJj(θ)边界条件、初始条件损失、观测数据的误差等。 具体如下 边界条件损失 J B C h J_{BC_h} JBCh、 J B C u J_{BC_u} JBCu 边界条件损失分别衡量流速 和水深 ℎ在边界上的预测值与已知边界条件和ℎ的偏差。值越小说明边界符合要求。 其中 J B C u ( θ ) 1 N B C u ∑ i 1 N B C u ∣ u ^ ( x i , t i , θ ) − u ( x i , t i ) ∣ 2 J_{BC_u}(\theta) \frac{1}{N_{BC_u}} \sum_{i1}^{N_{BC_u}} \left| \hat{u}(x^i, t^i, \theta) - u(x^i, t^i) \right|^2 JBCu(θ)NBCu1∑i1NBCu u^(xi,ti,θ)−u(xi,ti) 2 J B C h ( θ ) 1 N B C h ∑ i 1 N B C h ∣ h ^ ( x i , t i , θ ) − h ( x i , t i ) ∣ 2 , x ∈ Ω B C , t ∈ [ 0 , T ] J_{BC_h}(\theta) \frac{1}{N_{BC_h}} \sum_{i1}^{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] JBCh(θ)NBCh1∑i1NBCh h^(xi,ti,θ)−h(xi,ti) 2,x∈ΩBC,t∈[0,T] N B C u N_{BC_u} NBCu与 N B C h N_{BC_h} NBCh是边界条件点的总数 ∑ i 1 N B C h \sum_{i1}^{N_{BC_h}} ∑i1NBCh与 ∑ i 1 N B C u \sum_{i1}^{N_{BC_u}} ∑i1NBCu对边界上的所有点求和。
初始条件损失 J S u J_{S_u} JSu和快照损失 J S h J_{S_h} JSh 初始条件或快照损失分别衡量 u ^ \hat{u} u^流速和 h ^ \hat{h} h^水深在初始时间 t0 或特定快照时刻的预测值与已知值的偏差。 J S u ( θ ) 1 N S u ∑ i 1 N S u ∣ u ^ ( x i , 0 , θ ) − u ( x i , 0 ) ∣ 2 J_{S_u}(\theta) \frac{1}{N_{S_u}} \sum_{i1}^{N_{S_u}} \left| \hat{u}(x^i, 0, \theta) - u(x^i, 0) \right|^2 JSu(θ)NSu1∑i1NSu u^(xi,0,θ)−u(xi,0) 2 J S h ( θ ) 1 N S h ∑ i 1 N S h ∣ h ^ ( x i , 0 , θ ) − h ( x i , 0 ) ∣ 2 , x ∈ Ω S J_{S_h}(\theta) \frac{1}{N_{S_h}} \sum_{i1}^{N_{S_h}} \left| \hat{h}(x^i, 0, \theta) - h(x^i, 0) \right|^2, \quad x \in \Omega_S JSh(θ)NSh1∑i1NSh h^(xi,0,θ)−h(xi,0) 2,x∈ΩS 其中 N S u N_{S_u} NSu与 N S h N_{S_h} NSh是初始条件或快照点的总数 x ∈ Ω S \quad x \in \Omega_S x∈ΩS定义了河道距离 x x x的范围。 快照条件是指在特定时间点上对系统状态。如水深 ℎ或流速u的已知值或观测值。此外快照条件可以是t0也可以是后续时间点的观测值提供额外的时空约束。观测数据损失 J o b s u J_{obs_u} Jobsu与 J o b s h J_{obs_h} Jobsh J o b s u ( θ ) 1 N o b s u ∑ i 1 N o b s u ∣ u ^ ( x i , t i , θ ) − u ( x i , t i ) ∣ 2 J_{obs_u}(\theta) \frac{1}{N_{obs_u}} \sum_{i1}^{N_{obs_u}} \left| \hat{u}(x^i, t^i, \theta) - u(x^i, t^i) \right|^2 Jobsu(θ)Nobsu1∑i1Nobsu u^(xi,ti,θ)−u(xi,ti) 2 J o b s h ( θ ) 1 N o b s h ∑ i 1 N o b s h ∣ h ^ ( x i , t i , θ ) − h ( x i , t i ) ∣ 2 , x ∈ Ω o b s , t ∈ [ 0 , T ] J_{obs_h}(\theta) \frac{1}{N_{obs_h}} \sum_{i1}^{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] Jobsh(θ)Nobsh1∑i1Nobsh h^(xi,ti,θ)−h(xi,ti) 2,x∈Ωobs,t∈[0,T] 其中 N o b s u N_{obs_u} Nobsu与 N o b s h N_{obs_h} Nobsh观测点数据总数在 x ∈ Ω o b s , t ∈ [ 0 , T ] \quad x \in \Omega_{obs}, t \in [0, T] x∈Ωobs,t∈[0,T]中 Ω o b s \Omega_{obs} Ωobs为河道范围T为模拟周期。 图(b)为一个大尺度网格单元内的子网格降尺度过程包括时间步和空间位置其目的是在粗网格Δ内生成更详细的水深和流速解。其中 时间从 t n t_n tn到 t n 1 t_{n1} tn1Δt为时间步长。 空间从 x i x_i xi到 x x 1 x_{x1} xx1 x o b s 1 x_{obs}^1 xobs1与 x o b s 2 x_{obs}^2 xobs2为检测站获取的信息。
实验
本文通过六个实验研究基于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 数据类型为 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, lr5e-4, ExistModel0, uhDir,wDir, wmffDir, useObsTrue, use_ffFalse):初始化 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, dtypeself.DTYPE) if isinstance(S, (int, float)) else S # 将坡度转换为 TensorFlow 张量self.b tf.constant(b, dtypeself.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(configtf.ConfigProto(allow_soft_placementTrue, log_device_placementTrue)) # 初始化会话支持软放置和设备日志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_min1e-4, clip_value_max50) # 限制水深范围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,methodL-BFGS-B,options{maxiter: 50000, # 最大迭代次数maxfun: 50000, # 最大函数评估次数maxcor: 50, # 最大校正向量数maxls: 50, # 最大线搜索步数ftol: 1.0e-10, # 函数值容差gtol: 0.000001}) # 梯度容差self.global_step tf.Variable(0, trainableFalse) # 全局步数训练时递增self.learning_rate tf.train.exponential_decay(lr, self.global_step, 5000, 0.9, staircaseFalse) # 学习率衰减self.optimizer_adam tf.train.AdamOptimizer(learning_rateself.learning_rate) # Adam 优化器self.train_op_adam self.optimizer_adam.minimize(self.loss, global_stepself.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_dicttf_dict) # 执行 Adam 优化一步if it % 10 0: # 每 10 步打印一次elapsed time.time() - start_time # 计算单步耗时loss_value self.sess.run(self.loss, feed_dicttf_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_dicttf_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]], dtypeself.DTYPE), dtypeself.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, stddevxavier_stddev), dtypeself.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 fHEC-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 2RK4 解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), axis1)])) # 流速边界值h_BC np.vstack((h_exact[:, 0:1], h_exact[:, np.argmin(np.abs(x - u * t), axis1)])) # 水深边界值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, useObsTrue, use_ffuse_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(fCase {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, hspace0.08, wspace0.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], cmaprainbow, levelslevels, alpha0.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],markero, cg, s12, clip_onFalse)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],facecolorsnone, edgecolorsk, markero, s15, clip_onFalse)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*, cr, s25, clip_onFalse)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, fontsize16, transformax0.transAxes) # 添加子图标签divider make_axes_locatable(ax0) # 创建颜色条布局cax divider.append_axes(right, size2%, pad0.05) # 添加颜色条轴cb fig.colorbar(cs, caxcax, orientationvertical) # 添加颜色条cb.ax.tick_params(labelsize14) # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14) # 设置偏移文本字体大小cb.set_label(Water depth (m), fontsize14) # 设置颜色条标签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], cmaprainbow, levelslevels, alpha0.8) # 绘制 PINN 水深ax.set_xticklabels([]) # 隐藏 x 轴标签ax.set_yticklabels([] if ax ! ax1 else ax.get_yticklabels()) # 仅第一个子图显示 y 轴标签ax.text(0.05, 0.9, label, fontsize16, transformax.transAxes) # 添加子图标签divider make_axes_locatable(ax4) # 创建颜色条布局cax divider.append_axes(right, size2%, pad0.05) # 添加颜色条轴cb fig.colorbar(cs, caxcax, orientationvertical) # 添加颜色条cb.ax.tick_params(labelsize14) # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14) # 设置偏移文本字体大小cb.set_label(Water depth (m), fontsize14) # 设置颜色条标签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], cmapbwr, levelsnp.linspace(-0.5, 0.5, 11),alpha0.8, extendboth) # 绘制 PINN 误差ax.set_xticklabels([]) # 隐藏 x 轴标签ax.set_yticklabels([] if ax ! ax5 else ax.get_yticklabels()) # 仅第一个子图显示 y 轴标签ax.text(0.05, 0.9, label, fontsize16, transformax.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], cmaprainbow, levelslevels, alpha0.8) # 绘制插值水深ax.set_xticklabels([]) # 隐藏 x 轴标签ax.set_yticklabels([] if ax ! ax9 else ax.get_yticklabels()) # 仅第一个子图显示 y 轴标签ax.text(0.05, 0.9, label, fontsize16, transformax.transAxes) # 添加子图标签divider make_axes_locatable(ax12) # 创建颜色条布局cax divider.append_axes(right, size2%, pad0.05) # 添加颜色条轴cb fig.colorbar(cs, caxcax, orientationvertical) # 添加颜色条cb.ax.tick_params(labelsize14) # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14) # 设置偏移文本字体大小cb.set_label(Water depth (m), fontsize14) # 设置颜色条标签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], cmapbwr, levelsnp.linspace(-0.5, 0.5, 11),alpha0.8, extendboth) # 绘制插值误差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, fontsize16, transformax.transAxes) # 添加子图标签divider make_axes_locatable(ax16) # 创建颜色条布局cax divider.append_axes(right, size2%, pad0.05) # 添加颜色条轴cb fig.colorbar(cs, caxcax, orientationvertical) # 添加颜色条cb.ax.tick_params(labelsize14) # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14) # 设置偏移文本字体大小cb.set_label(Error (m), fontsize14) # 设置颜色条标签else: # Case 1-4 或 5 的标准布局ax1 fig.add_subplot(gs[1, :2]) # 添加 PINN 子图cs ax1.contourf(xx, tt, h_pred[:Nt], cmaprainbow, levelslevels, alpha0.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, fontsize16, transformax1.transAxes) # 添加子图标签ax2 fig.add_subplot(gs[1, 2:]) # 添加参考或插值子图cs ax2.contourf(xx, tt, h_exact[:Nt], cmaprainbow, levelslevels, alpha0.8) # 绘制参考或插值水深ax2.set_xticklabels([]) # 隐藏 x 轴标签ax2.set_yticklabels([] if case_id 4 else ax2.get_yticklabels()) # 隐藏 y 轴标签若非 Case 5-6ax2.text(0.05, 0.9, (c) Interpolation if case_id 2 else (c) Ref, fontsize16, transformax2.transAxes) # 添加子图标签divider make_axes_locatable(ax2) # 创建颜色条布局cax divider.append_axes(right, size2%, pad0.05) # 添加颜色条轴cb fig.colorbar(cs, caxcax, orientationvertical) # 添加颜色条cb.ax.tick_params(labelsize14) # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14) # 设置偏移文本字体大小cb.set_label(Water depth (m), fontsize14) # 设置颜色条标签ax3 fig.add_subplot(gs[2, :2]) # 添加 PINN 误差子图cs ax3.contourf(xx, tt, h_pred[:Nt] - h_exact[:Nt], cmapbwr, levelsnp.linspace(-0.5, 0.5, 11),alpha0.8, extendboth) # 绘制 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, fontsize16, transformax3.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],cmapbwr, levelsnp.linspace(-0.5, 0.5, 11), alpha0.8, extendboth) # 绘制误差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-6ax4.text(0.05, 0.9, (e) Interpolation-Ref if case_id 2 else (e) Ref-Ref, fontsize16, transformax4.transAxes) # 添加子图标签divider make_axes_locatable(ax4) # 创建颜色条布局cax divider.append_axes(right, size2%, pad0.05) # 添加颜色条轴cb fig.colorbar(cs, caxcax, orientationvertical) # 添加颜色条cb.ax.tick_params(labelsize14) # 设置颜色条刻度字体大小cb.ax.yaxis.offsetText.set_fontsize(14) # 设置偏移文本字体大小cb.set_label(Error (m), fontsize14) # 设置颜色条标签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), sharexTrue, shareyTrue) # 创建沿河图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, labelreference) # 绘制参考曲线axes[k].plot(x, h_pred[tlist[k]] eles if case_id 2 else h_pred[tlist[k]], -r, linewidth2, labelPINN) # 绘制 PINN 曲线if case_id 2: # Case 3-6 添加插值曲线axes[k].plot(x, h_exact[tlist[k]] eles, -g, linewidth2, labelInterpolation) # 绘制插值曲线axes[k].fill_between(x.flatten(), eles, color0.7) # 填充地形axes[k].text(0.85, 0.85, ft{tlist[k]} h if case_id 2 else ft{tlist[k] * dt} s,fontsize15, transformaxes[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(loc2 if case_id 2 else 4, prop{size: 14})plt.tight_layout() # 调整布局plt.savefig(ffigures/case{case_id}/along_channel.pdf) # 保存沿河图plt.close() # 关闭图形fig.savefig(ffigures/case{case_id}/contour.pdf) # 保存轮廓图plt.close() # 关闭图形结论
论文作者首次尝试将PINN应用于大尺度河流模型降尺度至子网格中该框架可解一维SVE方程、通过傅里叶变换编码潮汐边界周期性变化、同化多种数据类型并集成到大型河流模型。实验证明PINN在洪水建模中有重要作用为解决洪灾问题应用提供了方法。
不足以及展望
作者在最后提出了论文中的不足
测试案例相对简单仅应用恒定摩擦系数未考虑河道几何形状影响。现实中常见的数据问题如观测不确定性、数据不一致和观测缺失等也未充分探讨且PINN在解决河口非线性相互作用时精度有待提高。
作者对未来工作的展望如下 未来可将PINN降尺度应用于更广泛、更现实的流动条件。考虑可变摩擦系数、河道几何形状影响开发更强大的抗噪声算法处理数据不确定性。探索扩展PINN到二维领域用于创建复杂流动区域的数字孪生以解决当前大尺度模型无法解决的问题。
