2. 一维波动方程的波速辨识

任务目标

给定一维波动方程中若干时空点的位移观测数据，利用物理信息神经网络（PINN）反演未知波速参数 \( c \)。观测数据由已知波速 \( c_{\text{true}} = 2.0 \) 的正问题生成，并叠加不同强度的高斯噪声，考察反演结果随噪声水平的变化。

问题设定与观测数据

考虑一维波动方程初边值问题：

\[ u_{tt} = c^2 u_{xx}, \quad x \in (0, 1), \ t \in (0, 1) \]

边界条件：\( u(0,t) = u(1,t) = 0 \)
初始条件：\( u(x,0) = \sin(\pi x) \), \( u_t(x,0) = 0 \)

该问题的解析解为 \( u(x,t) = \sin(\pi x)\cos(c\pi t) \)。取 \( c_{\text{true}} = 2.0 \)，在 \( x \in [0,1] \), \( t \in [0,1] \) 内均匀采样 \( N = 2000 \) 个时空点作为观测位置，计算解析解并添加高斯噪声 \( \mathcal{N}(0, \sigma^2) \)，其中 \( \sigma \) 分别取 0（无噪声）、0.01、0.05，生成三组观测数据。

PINN结构与损失函数

神经网络 \( \hat{u}(x,t;\theta) \) 以 \( (x,t) \) 为输入，输出位移。总损失为：

\[ \mathcal{L} = \mathcal{L}_{\text{data}} + \lambda \mathcal{L}_{\text{PDE}} \]

其中数据损失：

\[ \mathcal{L}_{\text{data}} = \frac{1}{N}\sum_{i=1}^{N} \left( \hat{u}(x_i,t_i) - u_{\text{obs}}^{(i)} \right)^2 \]

PDE残差损失（在配点 \( (x_r, t_r) \) 上计算，\( M=1000 \)）：

\[ \mathcal{L}_{\text{PDE}} = \frac{1}{M}\sum_{j=1}^{M} \left( \hat{u}_{tt}(x_r^{(j)}, t_r^{(j)}) - c^2 \hat{u}_{xx}(x_r^{(j)}, t_r^{(j)}) \right)^2 \]

波速 \( c \) 作为可训练参数参与优化。取 \( \lambda = 1.0 \)。

关键完整代码

import numpy as np
import torch
import torch.nn as nn

# ---------- 1. 生成观测数据 ----------
def generate_data(c_true, N=2000, sigma=0.0, seed=0):
    rng = np.random.default_rng(seed)
    x = rng.uniform(0, 1, N)
    t = rng.uniform(0, 1, N)
    u_exact = np.sin(np.pi * x) * np.cos(c_true * np.pi * t)
    noise = rng.normal(0, sigma, N)
    u_obs = u_exact + noise
    return x, t, u_obs

# ---------- 2. PINN 网络 ----------
class PINN(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(2, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 1)
        )
        # 待反演参数 c，初始值 1.0
        self.c = nn.Parameter(torch.tensor(1.0, dtype=torch.float32))

    def forward(self, x, t):
        xt = torch.cat([x, t], dim=1)
        return self.net(xt).squeeze(-1)

    def pde_residual(self, x, t):
        x = x.clone().requires_grad_(True)
        t = t.clone().requires_grad_(True)
        u = self.forward(x, t)
        u_t = torch.autograd.grad(u, t, grad_outputs=torch.ones_like(u), create_graph=True)[0]
        u_tt = torch.autograd.grad(u_t, t, grad_outputs=torch.ones_like(u_t), create_graph=True)[0]
        u_x = torch.autograd.grad(u, x, grad_outputs=torch.ones_like(u), create_graph=True)[0]
        u_xx = torch.autograd.grad(u_x, x, grad_outputs=torch.ones_like(u_x), create_graph=True)[0]
        return u_tt - self.c**2 * u_xx

# ---------- 3. 训练 ----------
def train_pinn(x_obs, t_obs, u_obs, sigma, epochs=3000, lr=1e-3):
    model = PINN()
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)

    x_t = torch.tensor(x_obs, dtype=torch.float32).unsqueeze(1)
    t_t = torch.tensor(t_obs, dtype=torch.float32).unsqueeze(1)
    u_t = torch.tensor(u_obs, dtype=torch.float32)

    # 配点（PDE残差位置）
    rng = np.random.default_rng(1)
    x_r = torch.tensor(rng.uniform(0, 1, 1000), dtype=torch.float32).unsqueeze(1)
    t_r = torch.tensor(rng.uniform(0, 1, 1000), dtype=torch.float32).unsqueeze(1)

    for epoch in range(epochs):
        optimizer.zero_grad()
        u_pred = model(x_t, t_t)
        loss_data = torch.mean((u_pred - u_t)**2)
        res = model.pde_residual(x_r, t_r)
        loss_pde = torch.mean(res**2)
        loss = loss_data + loss_pde
        loss.backward()
        optimizer.step()

        if (epoch+1) % 1000 == 0:
            print(f"epoch {epoch+1}, loss_data={loss_data.item():.2e}, "
                  f"loss_pde={loss_pde.item():.2e}, c={model.c.item():.4f}")

    return model.c.item()

# ---------- 4. 执行反演 ----------
if __name__ == "__main__":
    c_true = 2.0
    for sigma in [0.0, 0.01, 0.05]:
        x, t, u = generate_data(c_true, sigma=sigma)
        c_est = train_pinn(x, t, u, sigma)
        print(f"sigma={sigma:.2f}: c_est={c_est:.4f}, 相对误差={abs(c_est-c_true)/c_true*100:.2f}%")

运行方式

需要安装 numpy 和 torch（CPU版本即可）。将上述代码保存为 wave_pinn.py，在终端执行：

python wave_pinn.py

验证方法与判定标准

程序依次对三种噪声水平执行反演，每个噪声水平训练3000轮。判定通过条件：

  • 训练过程中损失函数（尤其是数据损失）持续下降并趋于稳定；
  • 无噪声（σ=0）时，反演波速 \( c \) 接近真值2.0，相对误差小于5%；
  • 随着噪声增大（σ=0.01→0.05），反演结果可能出现偏差增大，但应保持在合理范围（σ=0.05时相对误差通常不超过15%），且训练仍能收敛。

注意：由于网络随机初始化和Adam优化器的随机性，每次运行结果略有差异，属正常现象。若σ=0.05时误差明显偏大，可适当增加训练轮数至5000或调整学习率。

常见问题

  • 梯度计算失败：确保输入 x、t 在计算PDE残差前调用了 requires_grad_(True)，且网络输出为标量（已通过 squeeze(-1) 处理）。
  • 反演不收敛：尝试增大 lr 至 5e-3，或增加网络宽度（如每层80个神经元）。
  • 边界条件未显式加入：本实现通过数据点覆盖边界附近区域，由数据损失隐式约束边界行为。若需严格满足边界条件，可在损失中额外添加边界残差项。
← 上一章 📖 返回目录 下一章 →