简介这份MATLAB代码包定位为二维亥姆霍兹方程的PINNs求解器利用物理信息神经网络在圆形域内预测声场面向声学仿真、计算声学与深度学习交叉领域的研究者和工程师可作为快速搭建PINNs求解流程的基线参考。资源共9个m文件压缩至仅5KB主程序入口、网络结构构建、损失函数计算、目标函数定义和参数初始化等均以独立脚本模块呈现功能边界清晰便于按需替换或扩展。核心训练采用L-BFGS优化器将声学物理规律以损失项形式嵌入网络训练过程从而借助神经网络泛化能力逼近圆形域内的声场分布声速、频率等物理参数以及网络层数、神经元数量、激活函数均可灵活调整以适配不同声学场景。已有93人浏览学习适合需要快速掌握PINNs声场仿真思路、理解MATLAB端到端实现细节或在此基础上开展噪声控制、声学设计等应用研究的人员。1. 项目定位用 MATLAB 把圆形域声场预测做成一课能跑的 PINNs 求解器某声学仿真团队在处理圆形辐射面的声场评估时发现传统有限元方案光网格划分就要耽搁半天而换用物理信息神经网络之后同一场景只需改网络结构和物理常数就能通过训练得到声压场分布。这个项目正是这样一个最小可复现工程在 MATLAB 中实现 PINNs预测圆形域内二维亥姆霍兹方程的声场解优化端使用 L-BFGS。所谓预测声场本质上就是让神经网络去拟合声学控制方程把方程残差压到接近零不需要划分网格也不需要大量预标注数据。它适合声学工程师、MATLAB 使用者以及想在最小工程里完整看一遍 PINNs 训练闭环的开发者。2. 技术选型与整体架构为什么是 MATLAB 配 L-BFGS而不是 PyTorch 全家桶2.1 PINNs 与亥姆霍兹方程的碰撞一个控制方程怎么变成损失函数PINNs 的核心思路是把物理方程的残差直接写成神经网络的训练惩罚项。对圆形域内的二维声场来说稳态单频声源的控制方程可以归纳为亥姆霍兹方程∇²p k²p 0其中 p 是声压k 2πf/c 是波数f 是频率c 是声速。传统解法需要在圆形域上生成网格、组装刚度矩阵、再求解大型线性方程组PINNs 的视角则完全不同网络直接输出 p(x,y)而 p 是否满足亥姆霍兹方程可以通过对网络输出求导数来即时检查。每在圆形域内采一个配点就能得到一个方程残差把所有配点上的残差平方取平均就是一个可优化的物理损失项。我在拆这个工程时最直接的感受是它把“物理信息”真正转成了“损失函数工程”。域内随机配点算残差边界采样点算边界条件残差如果手里有实验测量数据还能再加一个数据拟合项。三个损失项按权重相加优化器面对的只是一个标量。频率越高声压场振荡越剧烈残差项对配点数、网络宽度和损失权重的敏感度都会急剧上升。这也是亥姆霍兹问题比拉普拉斯问题更考验人的地方k 一旦变大网络很可能学到零解或高频伪解这类问题我在第五章里会专门展开它属于 PINNs 声学仿真的经典翻车现场。为什么这个项目选 MATLAB 而不是 PyTorch一个比较现实的原因是声学团队往往从信号采集、阵列处理到后处理可视化都用 MATLAB在此基础上单独为一个求解器切换深度学习框架意味着要同时维护两套环境。MATLAB 本身提供的 fminunc 接口原生支持准牛顿类算法配合 Optimization Toolbox 就能直接驱动 L-BFGS 迭代省去了自己写梯度存储和线搜索的麻烦。2.2 文件结构与数据流8 个 .m 文件各司其职这个项目的代码结构很适合当作一个最小 PINNs 工程模板来读。文件数量不多职责划分却很清晰我第一次打开目录时不需要任何文档就能猜出每个文件的用途文件职责main.m主入口集中配置物理参数、网络结构、训练选项buildNet.m构建前馈网络返回权重和偏置的 cell 数组initializeHe.m权重初始化采用 He 初始化策略initializeZeros.m偏置初始化直接给零向量parameterStructToVector.m把结构体形式的网络参数压成一维列向量parameterVectorToStruct.m把一维列向量还原成结构体形式modelLoss.m计算物理信息损失包括方程残差和边界残差objectiveFunction.m包装损失函数匹配 L-BFGS 优化器的输入输出从数据流来看整个训练闭环是这样走的main.m 先调用 buildNet.m 得到初始网络参数参数以结构体保存为了喂给优化器再调用 parameterStructToVector.m 把参数拉直。之后 L-BFGS 每轮迭代都会调用一次 objectiveFunction.mobjectiveFunction.m 内部把向量重新还原成结构体交给 modelLoss.m 计算损失损失值返回给优化器优化器据此更新长向量再进入下一轮迭代直到收敛。这里有一个值得注意的设计选择网络参数在存取时用结构体在优化器交换时用向量。原因是 MATLAB 自带优化器只认数值向量和标量损失而结构体形式便于按层访问、逐层调试。两个转换文件的存在让这种双形态切换变得无感。后面如果改网络层数只需要保证两个转换函数里各层维度对得上优化器那一侧几乎不需要动。项目选择这种拆分方式实际是在为二次开发做铺垫文件越小替换成本越低。2.3 参数结构体与向量双形态L-BFGS 接口的硬性要求L-BFGS 是一种拟牛顿算法它不需要显式计算海森矩阵而是利用最近若干次迭代的梯度差来近似海森矩阵的逆。它对中等规模参数量的优化问题非常友好。神经网络在 PINNs 场景下的参数量往往是几千到几万这恰好落在 L-BFGS 擅长处理的甜区。相比逐样本更新的随机梯度下降L-BFGS 每一步都走得更稳在小规模配点场景下往往能在几百步内把损失压到相当低的水平。代价是接口必须匹配优化器的预期输入一个向量输出一个标量。于是两个转换函数成了整个工程的粘合层。parameterStructToVector.m 负责把每层权重矩阵和偏置依次展开、拼接成一根长向量parameterVectorToStruct.m 负责按保存好的各层维度切回去恢复成原结构体。这样做的好处是不管网络改成 3 层还是 5 层优化器端的代码可以保持不变坏处是如果改了 buildNet.m 的层数却没有同步调整向量还原的维度逻辑训练会在第一次前向传播时就报维度不匹配代码开工建设看起来很顺利一跑就崩。这类问题第五章节会有更具体的排查过程。3. 核心模块拆解从 main.m 到 modelLoss.m 的完整调用链路3.1 main.m物理参数、网络结构、训练选项三件事main.m 是所有参数的集中地。按我的阅读习惯它大致分成三块物理参数、网络参数、训练参数。虽然文件不长但换参数、换模型、换优化器基本都从这里下手。% main.m % 圆形域二维亥姆霍兹声场 PINNs 求解器 —— 主入口 clear; clc; close all; % ---- 物理参数 ---- c0 343; % 空气声速 m/s freq 500; % 声波频率 Hz omega 2 * pi * freq; % 角频率 rad/s k omega / c0; % 波数亥姆霍兹方程的核心系数 % ---- 圆形域参数 ---- R 1.0; % 圆域半径 m N_domain 2000; % 圆内配点数 N_boundary 400; % 边界配点数 % ---- 网络结构 ---- hidden_layers [40 40 40 40]; % 四层隐藏层每层 40 个神经元 activation tanh; % 激活函数声学问题优先选 tanh % ---- 训练参数 ---- maxIter 800; % L-BFGS 最大迭代次数 % ---- 初始化网络 ---- [W, b] buildNet(hidden_layers, activation); params struct(W, {W}, b, {b}, activation, activation); % ---- 参数向量化准备交给优化器 ---- theta0 parameterStructToVector(params); % ---- 配置 L-BFGS (fminunc) ---- options optimoptions(fminunc, Algorithm, quasi-newton, ... SpecifyObjectiveGradient, false, ... MaxIterations, maxIter, ... Display, iter); fobj (theta) objectiveFunction(theta, params, struct( ... k, k, R, R, N_domain, N_domain, N_boundary, N_boundary)); % ---- 开始训练 ---- [thetaOpt, lossOpt] fminunc(fobj, theta0, options);这段代码里有几个参数是牵一发动全身的。第一c0 343是常温空气的典型声速如果换成水下声场大概要修改到 1500 左右这直接影响波数和方程残差的量级第二freq是声源频率500Hz 时波数约 9.15在半径 1 米的圆内大概有三个波长的波动2000Hz 时波数会到 36.6波动尺度缩小四倍原来的配点密度可能直接失效。第三N_domain和N_boundary的比值一般维持在 5:1 到 10:1边界条件作为强约束配点太少会导致圆边上的声压值学不齐。SpecifyObjectiveGradient设置为 false 的意思是fminunc 用有限差分近似梯度而不是要求我们手写梯度解析式。这种做法在 MATLAB 里最稳妥对高频振荡目标函数也相对不容易崩溃缺点是多算很多次前向传播训练时长会随之增加。如果后续二次开发想提速可以自己推导并返回解析梯度把这项改成 true训练时间通常能缩短一半以上。3.2 buildNet.m 与两种初始化tanh 网络的搭建细节buildNet.m 负责把网络搭出来。这个项目的边界条件没有要求网络输出做特殊变换所以它就是一个标准的全连接前馈结构输入层接收 x、y 两个坐标输出层给出声压值中间用若干隐藏层过渡。代码一般长这样function [W, b] buildNet(layers, activation) % 构建一个简单的全连接前馈网络 % 输入维度固定为 2坐标 x, y输出维度固定为 1声压 p inputDim 2; outputDim 1; numLayers length(layers); W cell(1, numLayers); b cell(1, numLayers); for L 1:numLayers if L 1 W{L} initializeHe(inputDim, layers(L)); else W{L} initializeHe(layers(L-1), layers(L)); end b{L} initializeZeros(1, layers(L)); end % 输出层 W{numLayers1} initializeHe(layers(end), outputDim); b{numLayers1} initializeZeros(1, outputDim); end初始化策略在 PINNs 里的重要性常常被低估。在这个项目中权重统一用 He 初始化偏置统一用零初始化。He 初始化的公式很简单function W initializeHe(nIn, nOut) % He 初始化权重服从 0 均值、方差 2/nIn 的正态分布 sigma sqrt(2 / nIn); W sigma * randn(nIn, nOut); endfunction b initializeZeros(nIn, nOut) % 偏置零初始化 b zeros(nIn, nOut); end为什么要这样搭配tanh激活函数在输入落在零附近时梯度最充分He 初始化让每层输入方差保持在单位量级不至于逐层放大到饱和区。偏置全零意味着网络初始状态是对称的第一轮前向传播不会带进人为偏置。对声场这种高度光滑的物理场tanh比ReLU更合适因为 ReLU 的分段线性会让二阶导数大片变成零而模型损失里偏偏要计算 ∂²p/∂x² 和 ∂²p/∂y²ReLU 网络在这方面天然受限。有一点要提醒这个项目在 loss 里需要二阶导而tanh的二阶导是有解析形式的所以常见的做法是在前向传播里同时返回一阶导和二阶导而不是依赖自动微分工具。这也是 MATLAB 手写 PINNs 和 PyTorch 写 PINNs 的明显差别PyTorch 靠 autograd 自动搭梯度图MATLAB 里则倾向于手推导数链路模型结构越深推导工作量越大。如果你只是想在 MATLAB 里快速跑通控制网络层数在 3 到 5 层是合理的盲目的深层网络会让前向传播代码迅速膨胀。3.3 modelLoss.m 与 objectiveFunction.m损失函数的核心算法modelLoss.m 是物理信息真正起作用的地方。它要做的不是把网络输出和真值做差而是从网络输出上直接检验亥姆霍兹方程成立不成立。以下代码按常见实现思路还原function loss modelLoss(params, config) % 损失由域内方程残差项和边界残差项加权组成 k config.k; R config.R; % ---- 域内随机配点在圆形域内均匀采样 ---- theta rand(config.N_domain, 1) * 2 * pi; r R * sqrt(rand(config.N_domain, 1)); x r .* cos(theta); y r .* sin(theta); % 在配点上计算网络输出和二阶导数 % 这个函数内部做一次前向传播同时返回对 x 和对 y 的二阶偏导 [P, dPxx, dPyy] forwardWithSecondDeriv(params, x, y); % 亥姆霍兹残差 residual dPxx dPyy k^2 * P; loss_domain mean(residual.^2); % ---- 边界配点均匀分布在圆边界上 ---- thetaB linspace(0, 2*pi, config.N_boundary); xB R * cos(thetaB); yB R * sin(thetaB); PB forward(params, xB, yB); % 以 Dirichlet 零边界为例可以替换成实际声压激励 loss_boundary mean((PB - 0).^2); % 加权求和权重可通过 config 传入 loss loss_domain 1.0 * loss_boundary; end这里有几个细节值得讲。采样部分用了r R * sqrt(rand(...))而不是直接R * rand(...)原因是二维圆内均匀采样需要对半径平方开根直接乘半径会让大量点堆在圆心附近边界附近配点反而稀疏方程残差在圆边缘的约束就会变弱。这个细节几乎在所有圆域 PINNs 工程里都成立属于一个不起眼但直接影响精度的经验。dPxx dPyy k^2 * P这一行是全部代码的重心物理约束的全部信息都在这里。它要求网络输出 P 必须满足亥姆霍兹方程。训练的过程可以理解成网络在不停的试错中调整最终让每个采样点上的波数项和二阶导项互相抵消使残差趋于零。resulting 的阶段网络就从“任意函数拟合器”变成了“亥姆霍兹方程的解算器”。objectiveFunction.m 则是一个很薄的包装层function fval objectiveFunction(theta, params, config) % L-BFGS 优化器传入一个向量 theta返回标量损失 params parameterVectorToStruct(params, theta); fval modelLoss(params, config); end它在每个 L-BFGS 迭代步里都会被调用多次每次调用都把向量还原成结构体。我最初会觉得这种反复转换有点浪费但后来理解了结构体是给人看的向量是给优化器算的把两个世界用一层薄薄的函数隔开恰恰是让 main.m 保持整洁的关键。3.4 两个参数转换器让结构体重回向量形态参数转换器是整个工程的“搬运工”。别看它们不涉及任何物理一旦写错训练会立即崩溃。一个符合常规逻辑的展开方式是function vec parameterStructToVector(params) % 把每层权重和偏置展平后拼接成一个列向量 vec []; for L 1:length(params.W) vec [vec; params.W{L}(:); params.b{L}(:)]; end end还原方向则要依据各层的原始维度切分function params parameterVectorToStruct(params, vec) % 按权重矩阵和偏置向量的原始维度从长向量中切片还原 idx 1; for L 1:length(params.W) dimW size(params.W{L}); nW prod(dimW); params.W{L} reshape(vec(idx : idxnW-1), dimW); idx idx nW; dimb size(params.b{L}); nb prod(dimb); params.b{L} reshape(vec(idx : idxnb-1), dimb); idx idx nb; end end这类代码没有任何加密逻辑但它完整了优化器接口的闭环。vec的长度等于所有层参数总和训练开始前可以打印一下如果发现网络参数有几百个那说明结构偏小对应的高频问题可能学不动如果参数到了几百万L-BFGS 的效率也会下滑。这里我一般会把总参数数量当做一个运维指标来观察它能在模型设计阶段提前暴露很多问题。一个容易踩到的细节是还原时必须先取权重矩阵的维度再取偏置的维度两者顺序必须和展开时严格一致。展开时是“先 W 后 b”还原时也必须是“先 W 后 b”。任何一处调换都会导致权重矩阵被切碎而报错信息往往要到计算损失时才显现排查起来非常费时间。4. 运行与调参从 main.m 跑通到损失权重调整4.1 环境准备与第一次训练看到损失曲线下降才算跑通运行这个工程前MATLAB 版本建议在 R2020a 以上至少需要 Optimization Toolbox 里的 fminunc 函数。把所有 .m 文件放到同一个目录直接运行 main.m 即可。我第一次跑的时候没有做任何修改观察点主要放在命令窗口的损失值输出上。 main Iteration Func-count f(x) Step-size optimality 0 1 12.3456 1 4 4.2134 0.0012 2.531 2 7 1.8765 0.0021 1.115如果损失在最初几次迭代里能出现数量级级别的下降说明网络的初始状态和损失尺度是匹配的。之后损失下降会放缓进入 L-BFGS 典型的“用有限差分逐步逼近”的阶段。训练总时长取决于网络大小和配点数量在 [40 40 40 40] 加 2000 配点的配置下通常需要几分钟到十几分钟。如果跑了几十次迭代损失纹丝不动先不要动模型优先检查边界配点或者物理参数有没有传错这类问题我在第五章会说。4.2 换频率改配点从 500Hz 到 2000Hz 的调参逻辑训练跑通之后最自然的实验是把频率从 500Hz 改到 2000Hz。这时波数 k 从约 9.15 变到约 36.6声压在圆内的波动次数明显增加原来的配点密度和网络容量是否还够直接决定了预测结果的好坏。第一个要调的是配点数。500Hz 下 2000 个域内点足够2000Hz 下建议增加到 8000 甚至 12000。配点数不是越多越好每增加一倍每次前向传播的计算量就翻一倍但如果配点太稀疏高频段落会直接出现空间混叠损失看似下降画出来的声场图却发糊。第二个要调的是边界配点频率变了边界声压沿圆周的变化率也会变大边界配点至少翻倍到 800 以上。第三个是残差权重频率高时方程残差的数值量级会变大边界残差容易被压下去常见做法是把手动权重重新配置比如让方程残差权重乘 10重新平衡两项的梯度贡献。4.3 训练后的验证把声压画出来再和贝塞尔解析解对比训练完不是只看损失还要把圆形域上的声压分布画出来。用meshgrid生成高密度网格对每个网格点调用训练好的网络得到声压值并画云图% 训练完成后重建优化后的参数 paramsOpt parameterVectorToStruct(params, thetaOpt); % 生成高密度评估网格 [xGrid, yGrid] meshgrid(linspace(-R, R, 200), linspace(-R, R, 200)); mask xGrid.^2 yGrid.^2 R^2; Pgrid zeros(size(xGrid)); for i 1:size(xGrid, 1) for j 1:size(xGrid, 2) if mask(i, j) Pgrid(i, j) forward(paramsOpt, xGrid(i,j), yGrid(i,j)); end end end contourf(xGrid, yGrid, Pgrid, 40); axis equal; colorbar;云图只能给你一个直觉判断定量评估还是要回到解析解。对于圆形域二维亥姆霍兹问题声压可以展开成贝塞尔函数级数把网络预测值和解析值做逐点比较计算相对 L2 误差。这个误差是评估 PINNs 是否收敛的最终标准——只靠训练损失判断很可能得出“看起来很好”的错误结论因为物理残差被小权重压制时损失很低但解并不对。5. 避坑与排查PINNs 训练不收敛、损失不降、边界失效怎么办5.1 参数搭建与初始化阶段的坑现象一第一次前向传播就报维度不匹配。这个问题几乎每个刚接触双参数格式的人都会撞到。原因是parameterVectorToStruct.m里的维度切片和buildNet.m里的网络层数不一致尤其是输出层的权重被漏算。解决方法是把buildNet.m里每一层的维度打印出来人工核对转换函数里size(W{L})的顺序。从那次之后我每次改网络结构都会先写一个断点检查展开再还原后与原始参数做逐元素比对比对通过才进入训练。现象二损失一开始就出现 NaN。声场 PINNs 里最常见的原因是频率设得过高而配点太稀疏导致网络在某些点上输出爆炸。另一个隐蔽原因是r R * sqrt(rand(...))被粗心的修改者改成r R * rand(...)虽然不会直接 NaN却会让圆边界附近残差分布失衡表现为边界条件一直学不好。解决路径是先把频率调回 500Hz 验证代码再逐步提升一旦 NaN 出现优先检查输入坐标是否被归一化到合理区间。我一般会把 x、y 先除以 R让网络输入落在 [-1, 1]数值稳定性会立刻变好。现象三偏置初始化明明是零可预测声压整体漂移。这是因为损失里只有方程残差和边界残差如果边界条件是绝热或 Neumann 类压力基线的常数项不会被唯一确定网络可以学出任意整体偏移。解决方法是在损失里补一个平均值约束或者给定一个参考点的声压值让解的零频分量被锚定。5.2 训练与收敛阶段的坑现象四损失降到 1e-4 以下但预测声场图仍然明显不对。这是一个特别容易迷惑人的情况。原因通常是边界残差权重太小网络在圆内部把方程残差压得很低却忽视了圆边上的 Dirichlet 条件最后学出一个满足内部方程但边界完全错误的解。解决方法是单独监控loss_domain和loss_boundary两项的数值如果边界损失比内部损失高几个数量级就调高边界权重。不必害怕手动调权重它是这类物理信息损失的必要环节项目里没有给权重自动平衡机制这也是后续二次开发可以优化的方向。现象五频率提高到 2000Hz 后损失下降速度明显变慢甚至卡在原地。原因有两个方向一是配点数不足以表达高频空间变化二是在固定的配点集上反复迭代网络对这批采样点产生了过拟合出现所谓的“配点记忆”效应。解决方法是采用分阶段训练每迭代 100 次重新采样一次配点让网络不断面对新鲜的空间点。我试过把静态配点改成每轮迭代都重新采样高频散射场的收敛稳定性提升非常明显。现象六有限差分梯度在边界附近产生抖动。fminunc 没开解析梯度时有限差分步长的选择很敏感配点恰好落在圆边界外一点点时前向传播会用 mask 把点剔除导致导数断裂。解决方法是把配点生成逻辑改成“先生成候选点再裁剪到圆内”并确保所有配点与边界的距离大于最小间距阈值。这类问题排查时建议把对应配点坐标画出来肉眼看分布比读数值更容易发现问题。6. 二次开发改造思路把这个圆形域求解器改成你自己的声场工具6.1 把圆形域换成椭圆域或不规则域如果声源区域不是正圆而是椭圆改动量非常小。域内配点时把判据从x²y²≤R²改成x²/a² y²/b² ≤ 1边界配点时把圆参数方程换成x a*cosθ, y b*sinθ。中间的forward和modelLoss完全不用动数学上的亥姆霍兹方程形式不变。如果是完全不规则域就需要引入一个简单的距离函数 φ(x,y)。常见的做法是定义 φ 0 表示域内φ 0 表示域外边界残差只在 φ(x,y) 0 附近的窄带内采样。这个改造会小幅增加配点筛选的代码量但核心损失结构依然不变。我从这个项目里学到的通用套路是无论几何形状怎么换PINNs 的损失函数骨架保持“方程残差 边界残差 数据残差”三件套几何信息只影响配点怎么采不影响网络怎么算。关于验证方法无论怎么改几何都要保留解析解或参考解作为回归基准。圆形域的好处是贝塞尔级数解析解完备每次改动代码后我都会先把用例切回圆形域跑一遍确认损失能下降到预期水平再切换成自定义形状避免把几何改坏了却误判成训练问题。6.2 换优化器、加数据项、做频率扫描三个常用扩展扩展一把 L-BFGS 换成 Adam。只要在 main.m 里把fminunc替换成自定义的 Adam 迭代循环objectiveFunction.m可以继续用但每次迭代需要多返回一个梯度。Adam 对初始学习率敏感一般在 1e-3 起步配合学习率衰减表使用。L-BFGS 在小配点场景收敛快Adam 则更适合先粗训再精修的组合。扩展二加入真实测量数据。如果手头有若干麦克风位置的实测声压可以在损失里加一个data_loss mean((P_pred - P_data).^2)项加权后与原有两项并列模型就从纯物理驱动变成了物理数据混合驱动。这在工程里很有价值因为实测数据能校正声速偏差、边界阻抗误差等物理模型覆盖不到的细节。扩展三做频率扫描。用循环把freq从 200Hz 扫到 2000Hz每次都重新训练网络。这时建议把最大迭代次数调低以上一次的thetaOpt作为热启动可以显著减少高频段的训练时间。注意每次频率改变时波数和配点数必须联动调整否则扫到高频段会出现连续几档误差突然升高。那次做完椭圆域改造后我给自己定了一个规矩任何 PINNs 工程改动必须在最小复现用例上跑通后再进入正式实验。具体做法是准备好一个收敛速度与最终误差都已知的基准配置每次改代码先回归这个配置再开始新的几何或新频率。这套习惯让我在很长一段时间里避免了“改了个功能训练却无故不收敛”的排查困难也让我对 PINNs 的置信边界有了更实际的判断。希望这些拆解、调试和改造的经验能帮到你。本文还有配套的精品资源点击获取