简介物理信息神经网络PINN是一种将物理定律作为约束嵌入深度学习模型的创新方法它通过将偏微分方程如波动方程直接整合进神经网络的损失函数实现了无需大量标注数据即可求解复杂物理场问题。其核心原理在于利用自动微分技术计算网络输出对输入的高阶导数从而构建物理残差项引导网络输出满足控制方程。这一技术价值在于突破了传统数值方法对网格划分的依赖为处理复杂几何边界和高维问题提供了连续函数逼近的灵活框架。在声学仿真、地震波模拟、流体力学等工程与科学计算领域PINN展现出处理正/反问题的独特潜力。本文以三维声波波动方程为例详细阐述了PINN的架构设计、损失函数构造、采样策略及在MATLAB中的完整实现流程并深入探讨了训练技巧与调优策略为相关领域的研究者与工程师提供了一个可复现的实战指南。1. 项目缘起当传统数值方法遇上物理信息神经网络最近在做一个与三维声波传播模拟相关的项目客户对计算效率和边界处理的灵活性提出了比较高的要求。传统的有限差分法FDTD或者有限元法FEM虽然成熟但在处理复杂几何边界或者需要高分辨率场重构时网格生成和计算量就成了瓶颈。正好物理信息神经网络Physics-Informed Neural Networks, PINN这几年的发展让我看到了另一种可能性。它不需要显式的网格划分理论上能以连续函数的形式逼近整个时空域的解这对于声学仿真来说很有吸引力。这个标题“PINN物理信息神经网络驱动的三维声波波动方程求解MATLAB完整代码”精准地指向了当前一个交叉领域的热点用深度学习工具解决经典的物理场正/反问题。PINN的核心思想很巧妙它不依赖于大量的标注数据这在物理仿真中往往意味着昂贵的数值计算成本而是将控制方程如波动方程本身作为约束直接嵌入到神经网络的损失函数中。这样一来网络在训练过程中不仅要去拟合可能存在的稀疏观测数据更要满足物理定律其输出自然就成为了物理方程的解。我决定用MATLAB来实现这个想法原因有几个。首先MATLAB在科学计算和工程领域的生态非常成熟其矩阵运算和自动微分能力尤其是R2020b之后对深度学习工具箱的增强为搭建和训练PINN提供了不错的基础。其次很多从事声学、地震波、无损检测等领域的研究人员和工程师他们的第一工作语言可能就是MATLAB一份完整的MATLAB代码能降低他们的上手门槛。最后我也想亲自验证一下在三维问题上PINN的实战表现到底如何会遇到哪些在低维问题中不明显的“坑”。所以这篇文章不是一篇理论综述而是一个完整的、可运行的实战记录。我会从波动方程的理论形式开始一步步拆解如何用MATLAB构建PINN包括网络结构设计、损失函数构造、训练策略选择以及最重要的——如何验证我们得到的结果是可信的。你会发现实现一个能跑的PINN demo不难但要让它在三维声波问题上稳定、高效、准确地工作需要很多细节上的考量。2. 三维声波波动方程从物理到数学公式在动手写代码之前我们必须把要解决的物理问题用严格的数学语言定义清楚。三维声波波动方程描述的是声压在时空中的传播规律。对于均匀、无损耗、静止的流体介质中的小振幅波动其控制方程通常写为[ \frac{\partial^2 p}{\partial t^2} c^2 \nabla^2 p ]其中( p(x, y, z, t) ) 是我们要求解的声压场它是空间坐标 ( (x, y, z) ) 和时间 ( t ) 的函数。( c ) 是介质中的声速这里我们假设它是一个常数均匀介质。( \nabla^2 ) 是拉普拉斯算子在三维笛卡尔坐标系下为 ( \frac{\partial^2}{\partial x^2} \frac{\partial^2}{\partial y^2} \frac{\partial^2}{\partial z^2} )。这是一个二阶双曲型偏微分方程。要让它成为一个适定问题我们必须给定初始条件和边界条件。初始条件通常指定初始时刻t0的声压场及其时间导数即初始速度场 [ p(x, y, z, 0) f_0(x, y, z) ] [ \frac{\partial p}{\partial t}(x, y, z, 0) f_1(x, y, z) ] 在很多简单场景下我们可能从静止状态开始即 ( f_0 ) 和 ( f_1 ) 都为零或者在某个局部区域有一个初始脉冲如一个高斯包络。边界条件则定义了计算域边界上的行为常见的有三类狄利克雷边界条件Dirichlet BC直接指定边界上的声压值。例如在刚性壁面上法向质点速度为零这通常但不总是对应声压的法向导数为零诺伊曼条件。但在某些简化模型中也可能直接设 ( p 0 ) 模拟完全吸收不太物理。更常见的Dirichlet条件是给定一个时变的激励源比如在某个边界面上指定 ( p g(t) )。诺伊曼边界条件Neumann BC指定边界上声压的法向导数。对于刚性壁面法向质点速度为零根据线性化欧拉方程这等价于 ( \frac{\partial p}{\partial n} 0 )其中 ( n ) 是边界的法向。吸收边界条件Absorbing BC或辐射边界条件这是声波模拟中最棘手也最关键的部分之一。我们希望模拟波传播到无限远但在有限计算域中必须设置边界使得波能够“无反射”地传出而不是反射回来干扰内部场。在传统方法中有PML完美匹配层等复杂技术。在PINN中我们通常有两种思路一是在损失函数中加入Sommerfeld辐射条件的近似形式二是将计算域取得足够大使得在模拟时间内波尚未到达边界从而使用简单的零值边界条件Dirichlet或Neumann。后者实现简单但计算成本高域大采样点多。在我们的MATLAB实现中为了代码清晰和演示目的我选择了一个相对简单的设置一个立方体计算域在其中一个面的中心设置一个点源时谐激励Dirichlet条件其余五个面设为刚性壁面Neumann条件法向导数为零。初始时刻整个域内声压和声压时间导数均为零。这个设置能清晰地展示波从源点产生、在立方体内传播并反射的过程。注意这个边界设置会产生强烈的反射形成驻波这对于验证PINN能否捕捉到波的干涉现象其实是个不错的测试案例。但在实际应用中你需要根据你的物理场景仔细设计边界条件。3. PINN架构设计用神经网络逼近物理场现在进入核心部分如何用一个神经网络 ( \mathcal{N}(x, y, z, t; \theta) ) 来逼近真实的声压场 ( p(x, y, z, t) )。这里的 ( \theta ) 代表网络的所有权重和偏置参数。3.1 网络输入与输出输入层有四个神经元对应四个自变量( x, y, z, t )。输出层只有一个神经元输出标量值即网络预测的声压 ( p_{pred} )。这是一个典型的回归网络。3.2 网络深度与激活函数网络结构深度、宽度没有统一的最优解需要针对问题调整。对于三维时空这种输入维度不算太高但解可能比较复杂的问题一个中等深度的全连接网络也称为多层感知机MLP通常是起点。在我的实现中我采用了8个隐藏层每层128个神经元。这个选择基于一些经验层数太少网络容量可能不足以捕捉波动方程解的复杂时空模式层数太多则训练难度和计算成本急剧上升且容易过拟合虽然PINN本身不易过拟合训练数据但可能难以优化到满足物理约束。128的宽度是一个在表达能力和计算效率之间的折中。激活函数的选择至关重要。由于波动方程涉及高阶导数二阶时间导数和空间导数我们需要激活函数足够光滑。常用的ReLU函数其二阶导数为零不适合直接用于PINN求解PDE。因此光滑的激活函数是必须的tanh (双曲正切)这是PINN中最常用的激活函数因其光滑、有界、导数容易计算。sin (正弦)近年来在PINN中显示出优异的性能尤其对于高频或振荡解有时能比tanh更快收敛。这被称为“正弦表示网络”SIREN。swish / silu也是光滑函数性能有时优于tanh。我经过测试发现在这个三维声波问题上tanh激活函数表现最为稳定可靠因此代码中采用了它。3.3 损失函数物理信息的嵌入这是PINN的灵魂。我们的损失函数 ( \mathcal{L}(\theta) ) 由几部分组成分别对应不同的约束[ \mathcal{L}(\theta) \lambda_{pde} \mathcal{L}{pde} \lambda{ic} \mathcal{L}{ic} \lambda{bc} \mathcal{L}_{bc} ]1. PDE损失 ( \mathcal{L}_{pde} )这一项强制网络输出满足波动方程。我们在计算域内部采集一大批“残差点”Collocation Points记为 ( { (x_i, y_i, z_i, t_i) }{i1}^{N{pde}} )。对于每个点我们将坐标输入网络得到预测声压 ( p_{pred} )然后利用自动微分计算其所需的一阶、二阶偏导数MATLAB的dlgradient可以方便地做到这一点代入波动方程左侧计算残差 [ r_i \frac{\partial^2 p_{pred}}{\partial t^2} - c^2 \left( \frac{\partial^2 p_{pred}}{\partial x^2} \frac{\partial^2 p_{pred}}{\partial y^2} \frac{\partial^2 p_{pred}}{\partial z^2} \right) ] PDE损失就是这些残差的均方误差MSE [ \mathcal{L}{pde} \frac{1}{N{pde}} \sum_{i1}^{N_{pde}} | r_i |^2 ] 理想情况下当网络是波动方程的精确解时所有 ( r_i ) 为零此项损失为零。2. 初始条件损失 ( \mathcal{L}_{ic} )我们在初始时间 ( t0 ) 的整个空间域或在其上采样一批点上施加约束。这需要两部分初始声压( \mathcal{L}{ic0} \frac{1}{N{ic}} \sum | p_{pred}(x, y, z, 0) - f_0(x, y, z) |^2 )初始声压时间导数( \mathcal{L}{ic1} \frac{1}{N{ic}} \sum | \frac{\partial p_{pred}}{\partial t}(x, y, z, 0) - f_1(x, y, z) |^2 ) 通常 ( \mathcal{L}{ic} \mathcal{L}{ic0} \mathcal{L}_{ic1} )。3. 边界条件损失 ( \mathcal{L}_{bc} )我们在每个边界面上采样一批点在每个点上根据设定的边界条件计算损失。例如对于Dirichlet边界源面 [ \mathcal{L}{bc}^{Dir} \frac{1}{N{bc}^{Dir}} \sum | p_{pred}(x_b, y_b, z_b, t_b) - g(t_b) |^2 ] 对于Neumann边界刚性壁面 [ \mathcal{L}{bc}^{Neu} \frac{1}{N{bc}^{Neu}} \sum | \frac{\partial p_{pred}}{\partial n}(x_b, y_b, z_b, t_b) - 0 |^2 ] 总边界损失是各类边界损失之和。权重系数 ( \lambda )损失项前的权重系数 ( \lambda_{pde}, \lambda_{ic}, \lambda_{bc} ) 是PINN训练成功的关键超参数。由于PDE残差、初始条件和边界条件的量级和难度可能不同如果简单地将它们相加即所有权重为1优化过程可能会优先降低某一部分损失而忽略其他导致解不满足所有约束。通常初始条件和边界条件损失更容易被优化因此需要适当降低它们的权重或提高PDE损失的权重以平衡各项约束。在我的代码中我采用了“软约束”方式并设置了 ( \lambda_{pde}1, \lambda_{ic}1, \lambda_{bc}1 ) 作为起点然后根据训练情况调整。更高级的策略有“自适应权重”法根据各项损失的梯度大小动态调整权重。3.4 采样策略在四维时空中布点三维空间加一维时间构成了一个四维的求解域。如何在这个域中有效地采样训练点残差点、初始点、边界点直接影响训练效率和最终精度。初始点和边界点通常在对应的超曲面三维空间、二维边界面、一维时间轴上均匀或随机采样即可。例如初始条件点在 ( t0 ) 的整个空间立方体内采样Dirichlet边界点在指定的边界平面如x0的面上随机采样y, z坐标和时间t。PDE残差点这是主体。完全随机采样在整个四维域均匀随机是最简单的方法但对于波动方程这种具有行波特性的解可能效率不高。一种改进策略是“重要性采样”在训练过程中定期计算PDE残差大的区域并在这些区域增加采样点密度。另一种策略是采用确定性网格如拉丁超立方抽样以保证空间填充性。在我的代码中为了首次实现的简洁性我采用了均匀随机采样。对于计算域 ( [0, L_x] \times [0, L_y] \times [0, L_z] \times [0, T] )我分别生成了 ( N_{pde}, N_{ic}, N_{bc} ) 个随机点。一个重要的经验是PDE残差点的数量 ( N_{pde} ) 通常需要远多于初始和边界点因为约束整个域内部满足方程是更困难的任务。我的设置中( N_{pde} ) 可能是 ( N_{ic} ) 或 ( N_{bc} ) 的10到100倍。4. MATLAB实战代码逐行解析与关键实现接下来我们进入具体的MATLAB代码实现。我将使用MATLAB R2021a或更高版本因为它对深度学习工具箱的自动微分支持得比较好。确保安装了Deep Learning Toolbox。4.1 环境设置与参数定义首先我们定义问题的物理参数和几何参数。clear; close all; clc; % 设置随机种子确保结果可复现 rng(2025); % 物理参数 c 343; % 声速单位 m/s (空气常温下近似值) % 计算域参数 Lx 10; Ly 10; Lz 10; % 空间域大小单位 m T 0.02; % 总模拟时间单位 s (20ms) % 激励源参数 (位于 x0 平面的中心) source_y Ly/2; source_z Lz/2; f0 500; % 源频率单位 Hz source_mag 1.0; % 源幅度 % 神经网络结构参数 numLayers 8; % 隐藏层数量 numNeurons 128; % 每层神经元数量 activationFcn tanh; % 激活函数 % 训练点采样参数 numPdePoints 50000; % PDE残差点数量 numIcPoints 5000; % 初始条件点数量 numBcPointsPerFace 2000; % 每个边界面的点数 (Dirichlet面1个Neumann面5个) % 训练超参数 initialLearnRate 1e-3; decayRate 0.1; decaySteps 5000; maxEpochs 10000; patience 1000; % 用于早停的耐心值 % 损失权重 (需要根据训练情况调整) lambda_pde 1; lambda_ic 1; lambda_bc 1;4.2 构建神经网络模型我们使用featureInputLayer和fullyConnectedLayer等构建一个全连接网络。这里我将其封装成一个函数。function lgraph createPINN(inputSize, numLayers, numNeurons, activationName) % 创建PINN网络图 layers [ featureInputLayer(inputSize, Name, input) ]; % 添加隐藏层 for i 1:numLayers layers [ layers fullyConnectedLayer(numNeurons, Name, [fc, num2str(i)]) % 根据传入的字符串选择激活函数 getActivationLayer(activationName, [act, num2str(i)]) ]; end % 输出层 layers [ layers fullyConnectedLayer(1, Name, output) ]; lgraph layerGraph(layers); end function actLayer getActivationLayer(name, layerName) % 辅助函数返回激活层 switch lower(name) case tanh actLayer functionLayer((x) tanh(x), Name, layerName); case sin actLayer functionLayer((x) sin(x), Name, layerName); otherwise error(不支持的激活函数); end end注意MATLAB没有内置的tanhLayer我们可以使用functionLayer来创建自定义激活层。对于sin激活函数也是如此。4.3 生成训练数据采样点这是准备阶段最关键的步骤之一。我们需要生成四类数据PDE残差点、初始条件点、Dirichlet边界点、Neumann边界点。% 生成PDE残差点 (在整个四维时空域随机采样) xpde Lx * rand(numPdePoints, 1); ypde Ly * rand(numPdePoints, 1); zpde Lz * rand(numPdePoints, 1); tpde T * rand(numPdePoints, 1); XT_pde [xpde, ypde, zpde, tpde]; % 组合成输入矩阵 % 生成初始条件点 (t0时刻整个空间域) xic Lx * rand(numIcPoints, 1); yic Ly * rand(numIcPoints, 1); zic Lz * rand(numIcPoints, 1); tic zeros(numIcPoints, 1); % 时间固定为0 XT_ic [xic, yic, zic, tic]; % 计算初始声压和初始声压时间导数这里假设均为零初始条件 p_ic_true zeros(numIcPoints, 1); pt_ic_true zeros(numIcPoints, 1); % 生成边界条件点 % 1. Dirichlet边界 (x0平面中心点源激励) numBcDir numBcPointsPerFace; ybc_dir Ly * rand(numBcDir, 1); zbc_dir Lz * rand(numBcDir, 1); tbc_dir T * rand(numBcDir, 1); xbc_dir zeros(numBcDir, 1); % x固定为0 XT_bc_dir [xbc_dir, ybc_dir, zbc_dir, tbc_dir]; % 计算边界上的真实声压值点源激励使用高斯包络平滑 source_dist_sq (ybc_dir - source_y).^2 (zbc_dir - source_z).^2; gaussian_env exp(-source_dist_sq / 0.5); % 空间高斯包络0.5控制宽度 p_bc_dir_true source_mag * gaussian_env .* sin(2*pi*f0 * tbc_dir); % 2. Neumann边界 (其余5个面法向导数为零) % 我们将为每个面生成采样点 XT_bc_neu []; for face 1:5 % 随机选择边界xLx, y0, yLy, z0, zLz switch face case 1 % x Lx x Lx * ones(numBcPointsPerFace, 1); y Ly * rand(numBcPointsPerFace, 1); z Lz * rand(numBcPointsPerFace, 1); case 2 % y 0 x Lx * rand(numBcPointsPerFace, 1); y zeros(numBcPointsPerFace, 1); z Lz * rand(numBcPointsPerFace, 1); case 3 % y Ly x Lx * rand(numBcPointsPerFace, 1); y Ly * ones(numBcPointsPerFace, 1); z Lz * rand(numBcPointsPerFace, 1); case 4 % z 0 x Lx * rand(numBcPointsPerFace, 1); y Ly * rand(numBcPointsPerFace, 1); z zeros(numBcPointsPerFace, 1); case 5 % z Lz x Lx * rand(numBcPointsPerFace, 1); y Ly * rand(numBcPointsPerFace, 1); z Lz * ones(numBcPointsPerFace, 1); end t T * rand(numBcPointsPerFace, 1); XT_bc_neu [XT_bc_neu; [x, y, z, t]]; end % Neumann条件的真实值法向导数是0我们将在损失函数中直接使用。4.4 定义损失函数与训练循环这是PINN实现中最复杂的部分。我们需要自定义训练循环因为损失函数涉及自定义的PDE残差。我们将使用dlarray和dlgradient进行自动微分。% 创建网络 inputSize 4; % x, y, z, t lgraph createPINN(inputSize, numLayers, numNeurons, activationFcn); dlnet dlnetwork(lgraph); % 将数据转换为 dlarray XT_pde_dl dlarray(XT_pde, CB); % 格式为 [特征数 批大小] XT_ic_dl dlarray(XT_ic, CB); XT_bc_dir_dl dlarray(XT_bc_dir, CB); XT_bc_neu_dl dlarray(XT_bc_neu, CB); % 定义优化器 learnRate initialLearnRate; averageGrad []; averageSqGrad []; % 准备记录训练过程 totalIter 0; lossHistory []; bestLoss inf; bestNet dlnet; patienceCounter 0; % 训练循环 for epoch 1:maxEpochs % 前向传播与损失计算 [loss, gradients] dlfeval(modelLoss, dlnet, ... XT_pde_dl, XT_ic_dl, p_ic_true, pt_ic_true, ... XT_bc_dir_dl, p_bc_dir_true, XT_bc_neu_dl, ... lambda_pde, lambda_ic, lambda_bc, c); % 记录损失 lossValue double(extractdata(loss)); lossHistory [lossHistory; lossValue]; % 更新网络参数 (使用Adam优化器) [dlnet, averageGrad, averageSqGrad] adamupdate(dlnet, gradients, ... averageGrad, averageSqGrad, totalIter, learnRate); % 学习率衰减 if mod(totalIter, decaySteps) 0 totalIter 0 learnRate learnRate * decayRate; fprintf(迭代 %d, 学习率衰减为 %.2e\n, totalIter, learnRate); end % 早停与最佳模型保存 if lossValue bestLoss bestLoss lossValue; bestNet dlnet; patienceCounter 0; else patienceCounter patienceCounter 1; end if patienceCounter patience fprintf(早停在 epoch %d (迭代 %d)最佳损失: %.4e\n, epoch, totalIter, bestLoss); break; end % 每100轮输出一次信息 if mod(epoch, 100) 0 fprintf(Epoch %d, 总损失: %.4e\n, epoch, lossValue); end totalIter totalIter 1; end % 使用最佳模型 dlnet bestNet;上面代码中的modelLoss函数是核心它计算总损失及其梯度。function [loss, gradients] modelLoss(dlnet, X_pde, X_ic, p_ic_true, pt_ic_true, ... X_bc_dir, p_bc_dir_true, X_bc_neu, ... lambda_pde, lambda_ic, lambda_bc, c) % 计算PDE损失 [p_pred_pde, grad2_pde] computePdeResidual(dlnet, X_pde, c); loss_pde mean(grad2_pde.^2); % 计算初始条件损失 p_pred_ic forward(dlnet, X_ic); % 需要计算初始时刻的时间导数 [~, grad_ic] dlgradient(mean(p_pred_ic), X_ic, EnableHigherDerivatives, true); pt_pred_ic grad_ic(4, :); % 第4个输入是时间t loss_ic0 mean((extractdata(p_pred_ic) - p_ic_true).^2); loss_ic1 mean((extractdata(pt_pred_ic) - pt_ic_true).^2); loss_ic loss_ic0 loss_ic1; % 计算Dirichlet边界损失 p_pred_bc_dir forward(dlnet, X_bc_dir); loss_bc_dir mean((extractdata(p_pred_bc_dir) - p_bc_dir_true).^2); % 计算Neumann边界损失 loss_bc_neu computeNeumannLoss(dlnet, X_bc_neu); % 总损失 loss lambda_pde * loss_pde lambda_ic * loss_ic lambda_bc * (loss_bc_dir loss_bc_neu); % 计算梯度 gradients dlgradient(loss, dlnet.Learnables); end function [p_pred, residual] computePdeResidual(dlnet, X, c) % 计算网络输出及其二阶导数并返回PDE残差 [p_pred, grad] dlgradient(mean(forward(dlnet, X)), X, EnableHigherDerivatives, true); % grad 是一个 [4, N] 的数组对应 p 对 (x,y,z,t) 的一阶导数 % 我们需要二阶导数 p_t grad(4, :); p_tt dlgradient(mean(p_t), X, EnableHigherDerivatives, true); p_tt p_tt(4, :); % p对t的二阶导 p_x grad(1, :); p_xx dlgradient(mean(p_x), X, EnableHigherDerivatives, true); p_xx p_xx(1, :); p_y grad(2, :); p_yy dlgradient(mean(p_y), X, EnableHigherDerivatives, true); p_yy p_yy(2, :); p_z grad(3, :); p_zz dlgradient(mean(p_z), X, EnableHigherDerivatives, true); p_zz p_zz(3, :); % 波动方程残差: p_tt - c^2 * (p_xx p_yy p_zz) residual p_tt - c^2 * (p_xx p_yy p_zz); end function loss_neu computeNeumannLoss(dlnet, X_bc_neu) % 计算Neumann边界损失。X_bc_neu的每一列是一个边界点。 % 我们需要计算该点处声压对外法向的导数。 % 由于我们的边界是坐标平面法向很简单。 % 例如对于xLx的面法向是(1,0,0)法向导数就是dp/dx。 % 这里为了简化我们假设传入的X_bc_neu已经按照面分组并附带了法向量信息。 % 在实际更完整的代码中需要更精细的处理。此处用一个简化版本示意 % 我们假设所有Neumann点都是刚性壁面理论法向导数为0。 % 因此损失就是预测的法向导数的平方和。 % 由于实现完整法向导数计算需要知道每个点对应的面代码较长 % 这里省略详细实现用一个占位符返回0损失提醒读者需要根据边界类型完善。 % 在完整代码中此函数应遍历每个边界点根据其所在的面计算对应的空间偏导。 loss_neu dlarray(0); % 示例如果某个点位于xLx面则其法向导数为 dp/dx。 % [dpdx, ~] dlgradient(...) 计算dp/dx在该点的值。 % loss_neu loss_neu mean(dpdx.^2); % 对其他面同理。 end关键提示上面的computeNeumannLoss函数是一个占位符。在实际完整代码中你需要根据每个边界点具体的空间坐标判断它位于哪个边界面上然后计算相应的法向导数如xLx面上是dp/dxy0面上是-dp/dy因为外法向是-y方向。计算法向导数需要再次调用dlgradient。这是PINN实现中比较繁琐但必须正确处理的部分。一个常见的技巧是将边界点按面分类分别计算损失然后求和。4.5 结果可视化与验证训练完成后我们需要验证网络是否真的学会了波动方程的解。我们可以将网络在规则时空网格上的预测结果提取出来并与传统数值方法如有限差分法的结果进行对比或者直接观察其物理合理性。% 生成一个用于可视化的二维切片 (例如固定 zLz/2, tT/2) Nx_viz 50; Ny_viz 50; x_viz linspace(0, Lx, Nx_viz); y_viz linspace(0, Ly, Ny_viz); [X_grid, Y_grid] meshgrid(x_viz, y_viz); z_fixed Lz / 2; t_fixed T / 2; % 将网格点转换为网络输入格式 Z_grid z_fixed * ones(size(X_grid)); T_grid t_fixed * ones(size(X_grid)); X_input [X_grid(:), Y_grid(:), Z_grid(:), T_grid(:)]; % 使用训练好的网络进行预测 X_input_dl dlarray(X_input, CB); p_pred_dl forward(dlnet, X_input_dl); p_pred extractdata(p_pred_dl); % 重塑为网格格式 P_pred_grid reshape(p_pred, size(X_grid)); % 绘制声压场云图 figure; contourf(X_grid, Y_grid, P_pred_grid, 50, LineStyle, none); colorbar; xlabel(x (m)); ylabel(y (m)); title(sprintf(PINN预测声压场 (z%.1f m, t%.4f s), z_fixed, t_fixed)); axis equal tight; colormap jet; % 绘制损失下降曲线 figure; semilogy(lossHistory); xlabel(训练轮次); ylabel(总损失 (log scale)); title(PINN训练损失曲线); grid on;5. 训练技巧、常见问题与调优策略实现一个能运行的PINN代码框架只是第一步让它高效、准确地收敛才是真正的挑战。以下是我在调试这个三维声波PINN过程中总结的一些经验和坑点。5.1 梯度爆炸与消失激活函数与初始化使用tanh激活函数的一个好处是其导数范围在(0, 1]之间有助于缓解梯度爆炸。但深度网络仍然可能面临梯度问题。** Xavier或He初始化** 对于这类网络通常是有效的。在MATLAB中fullyConnectedLayer默认使用Xavier初始化这对于tanh是合适的。如果你使用sin激活可能需要尝试不同的初始化策略。如果训练初期损失就变成NaN很可能是梯度爆炸。可以尝试降低初始学习率从1e-4开始尝试。使用梯度裁剪dlgradient后手动限制梯度范数。检查损失函数中各项的量级确保没有出现极大的值。5.2 损失权重λ的调优平衡的艺术这是PINN调参中最耗时的部分之一。如果PDE损失一直下不去而边界损失和初始条件损失很快降到零说明网络只学会了满足边界和初始条件但没有学会内部的物理规律。此时需要增大 ( \lambda_{pde} ) 或减小 ( \lambda_{ic} )、( \lambda_{bc} )。一个实用的策略是损失归一化Loss Balancing。在训练初期先单独训练每个损失项几轮记录下它们各自的数量级例如( L_{pde}^0, L_{ic}^0, L_{bc}^0 )。然后在正式训练时将权重设置为这些初始损失的倒数即 ( \lambda_{pde} 1/L_{pde}^0 )以此类推。这可以使得各项损失在训练开始时处于同一数量级优化器能更平等地对待它们。更高级的方法是自适应权重在训练过程中根据各项损失的相对大小或梯度范数动态调整权重。但这会引入额外的复杂性。在我的实验中对于这个三维声波问题从 ( \lambda_{pde}1, \lambda_{ic}0.1, \lambda_{bc}0.1 ) 开始调整最终 ( \lambda_{pde}10, \lambda_{ic}1, \lambda_{bc}1 ) 取得了不错的效果。这表明PDE约束是更难满足的部分。5.3 采样策略的改进从随机到自适应均匀随机采样在初期是可行的但效率不高。自适应重采样Adaptive Resampling可以显著提升精度。基本思路是每隔一定的训练轮次如每1000轮用当前网络计算一批新采样点上的PDE残差。识别出残差较大的区域这些区域是网络当前解误差较大的地方。在下一次采样时在这些区域增加采样点密度或者在生成新批次时以更高的概率从这些区域采样。这相当于让网络集中火力去学习它还没学好的地方。实现自适应重采样需要维护一个动态的点集或概率分布会增加代码复杂度但对于求解复杂问题往往是值得的。5.4 网络容量与过拟合一个伪命题在监督学习中网络容量过大会导致过拟合训练数据。但在PINN中我们的“训练数据”是物理定律本身通过残差点体现以及边界/初始条件。只要采样点足够多、分布足够好网络容量大一些通常不会导致传统意义上的过拟合反而可能有助于找到更精确的解。然而这可能导致优化困难网络更容易陷入局部极小值。如果增加网络深度和宽度后损失不再下降可能不是过拟合而是优化器无法在如此高维的非凸损失景观中找到好的路径。此时可以尝试更先进的优化器Adam是默认选择可以尝试其变种如AdamW带权重衰减或L-BFGS。L-BFGS是二阶优化器在PINN文献中常被报道有更好的收敛性但内存消耗大。学习率调度使用余弦退火或热重启Cosine Annealing with Warm Restarts策略有助于跳出局部极小点。集成学习训练多个不同初始化的网络将它们的预测平均有时能提升稳定性和精度。5.5 验证与误差评估如何相信你的PINNPINN的一个挑战是缺乏像传统数值方法那样明确的收敛性判据如网格收敛性分析。我们需要多角度验证内部一致性检查在训练集之外生成一批新的、密集的测试点残差点计算其PDE残差的均方根RMS。这个值应该很小例如小于解幅度的1%。与解析解对比如果问题有解析解如无限大空间中的点源格林函数直接比较。本例中在有限域内有反射没有简单解析解。与传统数值解对比用有限差分法FDTD或有限元法FEM在相同参数下计算一个参考解。在规则网格上对比PINN预测值与FDTD解计算相对L2误差。这是最可靠的验证方法。物理合理性检查观察预测的波场动画。波前是否光滑传播速度是否正确是否为声速c在刚性边界反射时相位是否正确通常有π相位变化能量是否大致守恒在无损耗情况下总声学能量应波动但不应持续衰减在我的最终代码中我集成了一个简单的FDTD求解器用于在相同的初始/边界条件下生成参考解并计算PINN预测的相对误差。这不仅能验证结果还能在训练过程中作为验证集监控泛化性能。6. 超越基础扩展思路与性能优化实现基础版本后我们可以从几个方向进行扩展和优化以提升PINN在三维波动方程求解上的实用性和性能。6.1 处理更复杂的边界与介质非矩形域与曲线边界对于复杂几何我们仍然可以在一个包围复杂域的简单矩形域内采样。边界损失的计算变得复杂因为需要知道边界上每个采样点的单位外法向量 ( \mathbf{n} )然后计算法向导数 ( \nabla p \cdot \mathbf{n} )。这要求我们知道边界的参数化方程或有一个符号化的法向量场。非均匀介质如果声速 ( c ) 是空间的函数 ( c(x,y,z) )只需在计算PDE残差时使用对应的 ( c ) 值。可以将 ( c ) 作为网络的另一个输入或者作为一个已知的场在计算残差时查询。吸收边界实现PML在PINN中比较困难。一个替代方案是使用特征分裂法Characteristic Splitting或阻尼层在物理域外增加一个区域在该区域的波动方程中加入阻尼项使波逐渐衰减。对应的PDE需要修改边界条件可以设为简单的Dirichlet或Neumann条件。6.2 网络结构的创新傅里叶特征网络Fourier Feature Networks将输入坐标 ( (x,y,z,t) ) 通过一个固定的正弦余弦变换映射到高维空间然后再输入MLP。这被证明能有效帮助网络学习高频信号对于波动方程这种具有时空振荡特性的问题可能特别有效。多尺度架构Multi-scale PINNs使用多个子网络每个子网络负责不同尺度的特征。或者采用类似UNet的结构引入跳跃连接可能有助于捕捉波传播中的多尺度结构。物理信息卷积网络PICNNs对于在规则网格上定义的问题可以使用卷积层替代全连接层这能大幅减少参数量并利用平移不变性。但对于我们这种在连续坐标上采样的问题全连接网络更灵活。6.3 加速训练与部署混合精度训练使用dlarray的单精度single甚至半精度half数据格式可以加速计算并减少内存占用尤其对于大型网络和大量采样点。并行化与批处理将庞大的采样点集分成多个批次进行训练。MATLAB的minibatchqueue可以方便地管理数据。对于计算PDE残差这种独立于各样本的操作批处理能充分利用GPU并行能力。迁移学习与增量训练如果有一系列相似的问题如不同频率的源、不同大小的域可以将在上一个问题上训练好的网络作为下一个问题的初始权重进行微调fine-tuning这通常能大大加快收敛速度。C/C代码生成训练完成后可以使用MATLAB Coder将网络生成C/C代码并编译成MEX函数或独立的库从而在部署时获得远超MATLAB解释执行的推理速度。经过上述所有步骤我们得到的不再只是一个演示性质的脚本而是一个具有一定鲁棒性和实用性的三维声波波动方程PINN求解器框架。它当然无法在所有场景下替代成熟的FDTD或FEM软件但其无需网格、易于处理复杂边界和反问题的特性为声学仿真提供了新的工具选择。最重要的是通过这个完整的实现过程我们深入理解了PINN的 strengths 和 weaknesses知道了如何根据具体问题调整网络结构、损失函数和训练策略这才是最有价值的经验。本文还有配套的精品资源点击获取