文獻閱讀

題目信息

• 題目： Physics-Informed Neural Network Approach for Solving the One-Dimensional Unsteady Shallow-Water Equations in Riverine Systems
• 期刊： Journal of Hydraulic Engineering
• 作者： Zeda Yin, S.M.ASCE; Jimeng Shi; Linlong Bian, S.M.ASCE; William H. Campbell; Sumit R. Zanje, S.M.ASCE; Beichao Hu; and Arturo S. Leon, M.ASCE
• 發表時間： 2025
• 文章鏈接： https://ascelibrary.org/doi/epdf/10.1061/JHEND8.HYENG-13572

摘要

數值方法在求解非線性偏微分方程時，實際應用中存在一定困難。傳統機器學習和深度學習模型依賴大量高質量訓練數據，數據成本高且難度大。此外，這些模型多為黑箱模型，計算過程難以解釋。盡管PINN近年來在多個領域取得成功，但在淺水方程及水文學和水力學領域的應用研究仍不充分。現有研究多集中于求解其他偏微分方程，且在考慮地形信息和摩擦的明渠水流問題上，尚無有效的PINN框架。基于以上背景，本文旨在提出一種新的PINN框架，以解決一維非穩態淺水方程，為水系統工程問題提供更有效的解決方案。文中詳細介紹PINN框架，包括前向步驟、損失函數構建和反向步驟，還對其進行改進以解決大規模問題。通過兩個案例驗證，結果表明PINN能準確預測流速、流量和水位，且可進行位置和時間外推，但存在訓練時間長和泛化性不足的局限。

創新點

該論文使用的PINN是無數據方法，不受數據獲取難題限制。且PINN將物理規律數學表達式融入框架，能進行位置和時間外推，提升極端條件下可靠性。

網絡架構

SVE由質量守恒方程和動量守恒方程組成，適用于任意形狀的橫截面，可寫為：
$$\frac{?\mathbf{U}}{?t}+\frac{?\mathbf{F}}{?x}=\mathbf{S}$$
其中，向量變量定義為：
$$\mathbf{U}=\left[\begin{array}{c}A\\ Q\end{array}\right],\mathbf{F}=\left[\begin{array}{c}Q\\ \frac{Q^2}{A}+g{I}_1\end{array}\right],\mathbf{S}=\left[\begin{array}{c}0\\ gA(S_0?S_f)\end{array}\right]$$
A: 橫截面濕周面積；Q: 橫截面流量；$g{I}_1$ : 靜水推力；$S_f$：摩擦坡度；$S_0$：地形高程坡度；t: 時間；x: 空間坐標。
其中：
$g\frac{?I_1}{?x}=gA\frac{?h}{?x}$
摩擦坡度通過Manning方程計算：
$S_f={\left(\frac{nv}{K{R}^{0.667}}\right)}^2$
地形高程坡度為：
$S_0=?\frac{dz}{dx}$
其中
n: Manning粗糙系數；v: 橫截面流速；K: 單位轉換因子；R: 水半徑；z: 河床高程。
以上是用向量表示，拆開來第一行（即質量守恒），第二行（動量守恒）簡化就是如下圖所示：

論文的PINN結構如下：
輸入： 輸入為空間x和時間t

輸出： 為流速𝑣及水位h

結構： 全連接多層感知器，8個隱藏層，每層80個隱藏單元，ReLU激活函數。

由于橫截面形狀復雜，參數（A,P,R,B）無法解析計算。
論文基于DEM數據，采用數值方法計算。

DEM（Digital Elevation Model，數字高程模型）是一種表示地表高程的數字化數據集，通常以網格或點云的形式存儲。
它記錄了地表在特定位置（通常以經緯度或投影坐標表示）的海拔高度
廣泛應用于地理信息系統、水文建模、地形分析、洪水模擬等領域。
DEM數據可以描述地形特征，如河床、坡度、山谷等，對于模擬水流、洪水傳播和地形相關計算至關重要。

DEM提供離散采樣點，每個點包含橫截面距離a和高程b。
然后將高程值轉換為相對于預測水深的坐標，如下圖所示：

濕周面積A通過梯形規則積分高程和距離：
$A=?{\sum }_{i=1}^{n?1}\frac{(b_{i+1}+b_i)?(a_{i+1}?a_i)}{2}$
濕周周長P累加相鄰點間的歐幾里得距離：
$P={\sum }_{i=1}^{n?1}\sqrt{(a_{i+1}?a_i{)}^2+(b_{i+1}?b_i{)}^2}$
水面半徑R和頂部寬度B直接從A,P和邊界點計算：
$R=\frac{A}{P}$
$B=a_m?a_k$
損失函數： 損失函數包括上游邊界、下游邊界、質量方程和動量方程損失組成。
由于各部分量級差異，比如，水位h為10²而偏導數$\frac{?h}{?t}$為10^-3-10^-8，直接求和會導致優化偏向某些分量。每個損失分量乘以權重（$W_1?W_4$），使其量級接近，權重通過邊界條件的量級和網格劃分，如下圖所示：

$${\mathcal{L}}_{\text{total}}=W_1{\mathcal{L}}_{{\text{up}}_{bc}}+W_2{\mathcal{L}}_{{\text{down}}_{bc}}+W_3{\mathcal{L}}_{{\text{physics}}_{mass}}+W_4{\mathcal{L}}_{{\text{physics}}_{momentum}}$$
PDE殘差項基于SVE：
${\mathcal{L}}_{\text{physics}}=\frac{1}{N_{\text{physics}}}{\sum }_{i=1}^{N_{\text{physics}}}{\left|\frac{?\mathbf{U}}{?t}+\frac{?\mathbf{F}}{?x}?\mathbf{S}\right|}^2$
上游邊界條件損失：
${\mathcal{L}}_{{\text{up}}_{bc}}=\frac{1}{N_{{\text{up}}_{bc}}}{\sum }_{i=1}^{N_{{\text{up}}_{bc}}}{\left|{\hat{v}}_{bc}?v_{bc}\right|}^2$
下游邊界條件損失：
${\mathcal{L}}_{{\text{down}}_{bc}}=\frac{1}{N_{{\text{down}}_{bc}}}{\sum }_{i=1}^{N_{{\text{down}}_{bc}}}{\left|{\hat{h}}_{bc}?h_{bc}\right|}^2$
其中：
$N_{\text{physics}}$為配點總數；$N_{{\text{up}}_{bc}}$與$N_{{\text{down}}_{bc}}$上下游邊界點數；$W_1?W_4$為權重；

實驗

實驗采用多階段訓練策略，如下圖所示：
分階段優化邊界和物理方程，模擬了先確定邊界再求解內部場的邏輯，加速收斂。論文將訓練分為三個階段優化。
第一階段為放大邊界條件權重，優先優化邊界收斂。
第二階段則平衡權重，降低學習率，優化SVE殘差。
第三階段略微增加SWEs權重，進一步提高精度。

論文通過假設和實際場景研究展示了PINN框架求解一維非定常淺水方程的性能，具體結果如下：

1. 假設的均勻梯形渠道中流量和水位突然變化的場景
   PINN框架和HEC - RAS輸出的速度和水位剖面吻合良好，所有橫截面速度和水位的平均絕對誤差分別為0.002743 m/s和0.001219 m。
   
   PINN能高精度求解流量基一維淺水方程，預測水位的平均絕對誤差與速度基PINN相近，約為0.0012 m，但流量的平均絕對誤差比速度大，這是因為流量尺度更大且未進行歸一化。
   

2. 休斯頓賽普拉斯溪下游的實際案例
   PINN框架和HEC - RAS輸出對比，流量預測的平均絕對誤差在0.0833 - 0.10083 m3/s之間，水位預測的平均絕對誤差在0.0152 - 0.079 m之間，PINN輸出趨勢對流量和水位略有低估。
   
   通過追蹤內部橫截面在每次迭代中的損失函數值，發現存在全局最小值，證明PINN框架理論上可獲得一維淺水方程的小殘差解。
   
   PINN在計算域內的位置外推結果準確，平均絕對誤差甚至小于非外推預測結果。
   
   在計算域外進行外推測試，與邊界位置的結果無顯著差異，表明PINN可在域外進行合理外推。
   
   以水位站記錄的歷史數據為參考，對未來30小時進行外推，外推結果的平均絕對誤差與非外推預測相似，趨勢與歷史曲線擬合良好。
   

代碼如下：

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
from uuid import uuid4# 設置隨機種子，確保可重復性
tf.random.set_seed(42)
np.random.seed(42)# 物理參數
g = 9.81  # 重力加速度 (m/s^2)
n = 0.013  # Manning粗糙系數
K = 1.0  # 單位轉換因子（SI單位）# DEM數據
def compute_hydraulic_params(h, x):# 假設梯形橫截面：底寬2m，邊坡1:1b0 = 2.0  # 底寬z = 1.0   # 邊坡A = h * (b0 + z * h)  # 濕周面積P = b0 + 2 * h * tf.sqrt(1 + z**2)  # 濕周周長R = A / P  # 液壓半徑B = b0 + 2 * z * h  # 頂部寬度S0 = -0.001  # 假設地形坡度（簡化）return A, P, R, B, S0# 定義神經網絡模型
class PINN(tf.keras.Model):def __init__(self):super(PINN, self).__init__()# 8層全連接網絡，每層80個單元，ReLU激活self.layers_list = [tf.keras.layers.Dense(80, activation='relu') for _ in range(8)]self.output_layer = tf.keras.layers.Dense(2)  # 輸出：v（流速），h（水位）def call(self, inputs):x, t = inputs[:, 0:1], inputs[:, 1:2]X = tf.concat([x, t], axis=1)for layer in self.layers_list:X = layer(X)return self.output_layer(X)  # [v, h]# 計算物理殘差
def compute_physics_loss(model, x, t):with tf.GradientTape(persistent=True) as tape:tape.watch([x, t])inputs = tf.concat([x, t], axis=1)vh = model(inputs)  # 預測 [v, h]v, h = vh[:, 0:1], vh[:, 1:2]# 計算一階偏導數dv_dx = tape.gradient(v, x)dh_dx = tape.gradient(h, x)dh_dt = tape.gradient(h, t)# 計算水文參數A, P, R, B, S0 = compute_hydraulic_params(h, x)# 摩擦坡度Sf = (n * v / (K * tf.pow(R, 2/3)))**2# 質量守恒方程殘差（簡化形式，公式18）res_mass = dh_dt + (A / B) * dv_dx + v * dh_dx# 動量守恒方程殘差（忽略Q^2/A項，簡化）res_momentum = g * A * (dh_dx + Sf - S0)# L2范數平方損失loss_mass = tf.reduce_mean(tf.square(res_mass))loss_momentum = tf.reduce_mean(tf.square(res_momentum))return loss_mass, loss_momentum# 計算邊界條件損失
def compute_bc_loss(model, x_bc, t_bc, v_bc, h_bc):inputs = tf.concat([x_bc, t_bc], axis=1)vh = model(inputs)v_hat, h_hat = vh[:, 0:1], vh[:, 1:2]# 上游流速損失loss_up_bc = tf.reduce_mean(tf.square(v_hat - v_bc)) if v_bc is not None else 0.0# 下游水位損失loss_down_bc = tf.reduce_mean(tf.square(h_hat - h_bc)) if h_bc is not None else 0.0return loss_up_bc, loss_down_bc# 訓練函數
def train_pinn(model, x_physics, t_physics, x_bc_up, t_bc_up, v_bc_up, x_bc_down, t_bc_down, h_bc_down, epochs, lr, weights):optimizer = tf.keras.optimizers.Adam(learning_rate=lr)for epoch in range(epochs):with tf.GradientTape() as tape:# 物理殘差loss_mass, loss_momentum = compute_physics_loss(model, x_physics, t_physics)# 邊界條件loss_up_bc, loss_down_bc = compute_bc_loss(model, x_bc_up, t_bc_up, v_bc_up, x_bc_down, t_bc_down, h_bc_down)# 總損失loss = (weights[0] * loss_up_bc + weights[1] * loss_down_bc +weights[2] * loss_mass + weights[3] * loss_momentum)gradients = tape.gradient(loss, model.trainable_variables)optimizer.apply_gradients(zip(gradients, model.trainable_variables))if epoch % 100 == 0:print(f"Epoch {epoch}, Loss: {loss.numpy():.6f}, "f"BC_up: {loss_up_bc.numpy():.6f}, BC_down: {loss_down_bc.numpy():.6f}, "f"Mass: {loss_mass.numpy():.6f}, Momentum: {loss_momentum.numpy():.6f}")# 主函數
def main():L = 79.25  # 渠道長度T = 600.0  # 時間范圍N_physics = 1000  # 物理配點數N_bc = 100  # 邊界點數x_physics = tf.random.uniform((N_physics, 1), 0, L)t_physics = tf.random.uniform((N_physics, 1), 0, T)# 生成邊界點t_bc = tf.random.uniform((N_bc, 1), 0, T)x_bc_up = tf.zeros_like(t_bc)  # 上游 x=0x_bc_down = L * tf.ones_like(t_bc)  # 下游 x=L# 模擬邊界條件v_bc_up = tf.sin(2 * np.pi * t_bc / T) * 0.5 + 0.5  # 上游流速尖峰h_bc_down = tf.sin(2 * np.pi * t_bc / T) * 0.2 + 0.5  # 下游水位尖峰# 初始化模型model = PINN()# 多階段訓練stages = [# 第一階段：優先優化邊界條件{"epochs": 5000, "lr": 1e-4, "weights": [1.0, 1.0, 0.1, 0.1]},# 第二階段：平衡物理和邊界{"epochs": 4000, "lr": 1e-5, "weights": [1.0, 1.0, 1.0, 1.0]}]for stage in stages:print(f"訓練階段：學習率={stage['lr']}, 權重={stage['weights']}")train_pinn(model, x_physics, t_physics,x_bc_up, t_bc, v_bc_up,x_bc_down, t_bc, h_bc_down,stage["epochs"], stage["lr"], stage["weights"])x_test = tf.linspace(0, L, 100)[:, None]t_test = tf.ones_like(x_test) * 300.0  # t=300sinputs = tf.concat([x_test, t_test], axis=1)vh_pred = model(inputs)v_pred, h_pred = vh_pred[:, 0], vh_pred[:, 1]plt.figure(figsize=(10, 5))plt.plot(x_test, v_pred, label="預測流速 (v)")plt.plot(x_test, h_pred, label="預測水位 (h)")plt.xlabel("x (m)")plt.ylabel("值")plt.legend()plt.savefig("pinn_result.png")plt.close()if __name__ == "__main__":main()

結論

本文提出并測試了用于求解一維非定常淺水方程的物理信息神經網絡。PINN框架能準確預測假設場景和歷史洪水場景結果，誤差小。可求解基于速度和流量的淺水方程，能對下游賽普拉斯溪案例在大流量下準確預測流量和水位。響應面表明PINN理論上可獲小殘差解，還能進行位置和時間外推，與參考數據高度吻合。

不足以及展望

PINN訓練時間長，因優化器找最優解難，且本文方法集成數值計算，GPU處理表現差。且PDE在邊界和初始條件不確定時有無限解，PINN只能在特定邊界條件下求近似解，不同條件需重新訓練。后續希望探索優化訓練算法或硬件加速方式，減少PINN訓練時間，如改進GPU對數值計算的處理能力。將邊界條件作為輸入變量，擴大淺水方程近似形式，使其能在一定邊界條件范圍內回歸方程。完善PINN模型，使其能處理涉及各種水工結構的復雜河流系統，減少對額外理論或經驗方程的依賴。