简介面向近场声全息与声场重构研究者的MATLAB仿真资源聚焦等效源法在实际声场重建中的应用解决从有限测量数据反推整个声场信息的逆问题。该方法适用于声学测试、噪声控制、声学成像等场景通过将复杂声源等效为虚拟源简化建模并借助正则化技术保证求解稳定性对理解声场重构原理与正则化函数选择有直接帮助。资源包共10个文件均为m脚本压缩包仅9KB包含主程序、正则化函数与辅助脚本可支撑从声源建模、441采样点设置、声压数据计算到声场重构的完整流程。目前已有1368人学习使用适合声学相关专业学生、科研人员及工程实践者参考。仿真中包含Tikhonov正则化、GCV、L曲线以及OMP等多种策略读者可通过对比实验观察不同正则化方法对重构精度和噪声抑制的影响并根据实际测量数据调整参数为后续研究或工程应用提供可直接修改的代码基础。1. 等效源法近场声全息为什么值得用 MATLAB 重新实现一遍近场声全息NAH在发动机罩、齿轮箱这类有限尺寸声源的噪声识别里有个老问题按平面波外推做重建测量一靠近复杂曲面重建声场就会出现边界抖动和倏逝波丢失几条亮纹根本不是真实声源。等效源法把问题改了个方向不直接对测量面做空间变换而是假设声源可以用一组假想源代替用测量声压反演源强再从源强正向叠加出任意位置的声场。反演和正向都在线性代数框架里完成恰好是 MATLAB 矩阵运算的主场。这篇文章从等效源法的声学原理说起给出 G 矩阵的构建方式、Tikhonov 正则化参数的选择方法再用可运行的 MATLAB 代码把一套最小闭环跑通。内容适合做 NVH 测试、声学仿真和状态监测的工程师你不需要预先掌握波数域滤波的全部推导但需要会用矩阵左除、能读懂奇异值分解结果并且愿意为正则化参数花一点调参时间。2. 等效源法近场声全息的声学模型与 MATLAB 基础实现2.1 为什么一组单极子能等效一个任意声源等效源法的前提是线性声学里的叠加原理任意稳态声场在某个区域的复声压都能用一组位于声源内部的简单源辐射场叠加来近似。最常用的是单极子它的自由场格林函数是[ g(\mathbf{r}, \mathbf{r}_s) \frac{e^{-jkR}}{4\pi R}, \quad R |\mathbf{r} - \mathbf{r}_s| ]其中 (k2\pi f/c) 是波数(\mathbf{r}_s) 是等效源位置。一个测量点上的声压是所有等效源贡献的线性叠加写成矩阵形式就是[ \mathbf{p} \mathbf{G} \mathbf{q} ]这里 (\mathbf{p}) 是测量面声压列向量(\mathbf{q}) 是等效源源强向量(\mathbf{G}) 是传递矩阵。第 (m) 行第 (n) 列元素就是第 (n) 个源到第 (m) 个测点的格林函数值。用单极子而不是偶极子、四极子是因为单极子的自由度已经足够描述大多数结构辐射声场。结构表面的振动辐射可以用表面法向振速替代而一组分布合理的单极子加上相位差就能生成方向性。偶极子只是两个相邻单极子的差值放在源强反演里会增加条件数却不一定带来精度收益。实际工程中先全部用单极子试算如果重建结果在高频段出现明显方向性误差再考虑引入少量偶极子。2.2 从测量声压反演源强为什么是个病态问题从 (\mathbf{p} \mathbf{G} \mathbf{q}) 求解 (\mathbf{q}) 时如果 (\mathbf{G}) 是方阵且满秩直接求逆就行。但等效源法里测量点数量通常大于等效源数量问题变成超定最小二乘即使不超定(\mathbf{G}) 也严重病态直接 (\mathbf{q} \mathbf{G}^{-1}\mathbf{p}) 会把测量噪声放大上百倍。病态的本质是近场声场中的倏逝波成分离声源越近倏逝波强度越高但它对应 (\mathbf{G}) 的小奇异值方向。数值上这些奇异值接近于零噪声在投影到这些小奇异值方向时会被除零放大。这就是为什么要在反演公式里加正则化。最常用的形式是Tikhonov正则化[ \mathbf{q} (\mathbf{G}^H\mathbf{G} \alpha \mathbf{I})^{-1}\mathbf{G}^H \mathbf{p} ]其中 (\alpha) 是正则化参数控制解的光滑程度。(\alpha) 太小噪声放大明显(\alpha) 太大重建结果退化成正向叠加的平均场丢失近场细节。2.3 构建 G 矩阵的最小 MATLAB 代码下面这个函数实现等效源法近场声全息的核心流程先根据测量点和等效源位置生成传递矩阵 (\mathbf{G})再做 Tikhonov 反演得到源强分布最后通过重建矩阵 (\mathbf{H}) 正向计算任意位置的声压。function [p_recon, q_est] esm_nah_core(p_meas, src_pos, meas_pos, recon_pos, k, alpha) % 等效源法近场声全息重建 % 输入: % p_meas - Mx1 复声压列向量测量面排序后的声压 % src_pos - Nsx3 等效源位置矩阵每行一个 [x,y,z] % meas_pos - Mx3 测量点位置矩阵 % recon_pos - Nrx3 重建点位置矩阵 % k - 波数 2*pi*f/c % alpha - Tikhonov 正则化参数 % 输出: % p_recon - Nrx1 重建声压 % q_est - Nsx1 估计的等效源源强 G green_matrix(meas_pos, src_pos, k); % MxNs 测量面到源面的传递矩阵 H green_matrix(recon_pos, src_pos, k); % NrxNs 重建面到源面的传递矩阵 Ns size(src_pos, 1); q_est (G*G alpha*eye(Ns)) \ G * p_meas; p_recon H * q_est; end function G green_matrix(obs_pos, src_pos, k) Ns size(src_pos, 1); M size(obs_pos, 1); G zeros(M, Ns); for ii 1:Ns R vecnorm(obs_pos - src_pos(ii, :), 2, 2); % 避免观测点与等效源重合时除零 R(R 1e-6) 1e-6; G(:, ii) exp(-1j*k*R) ./ (4*pi*R); end end这段代码的逻辑分两个阶段第一阶段用green_matrix生成测量面到源面的传递矩阵 (\mathbf{G})第二阶段用正则化最小二乘反演源强。注意G*G alpha*eye(Ns)里的 (\mathbf{G}^H\mathbf{G}) 是 Ns×Ns 方阵即使 M 比 Ns 大很多也能正常左除。green_matrix里的循环是逐等效源计算的在等效源数量几百个、测量点几千个时速度仍然可以接受。如果想提速可以改成向量化计算但要注意vecnorm在三维坐标差上的用法和 MATLAB 的隐式扩展规则。R(R 1e-6) 1e-6这行是关键保护实测中经常遇到等效源放在声源表面、测量点恰好有相同坐标的情况。2.4 坐标生成与调用示例实际使用时不直接手写坐标而是用meshgrid生成网格并重整成 N×3 矩阵。一个典型的三维测量面与等效源面配置如下% 测量面: x, y 从 -0.2 到 0.2 m, 阵元间距 0.02 m [x_m, y_m] meshgrid(-0.2:0.02:0.2, -0.2:0.02:0.2); meas_pos [x_m(:), y_m(:), zeros(numel(x_m), 1)]; % 等效源面: 位于 z -0.03 m, 间距取测量面的 1.5 倍 [x_s, y_s] meshgrid(-0.15:0.03:0.15, -0.15:0.03:0.15); src_pos [x_s(:), y_s(:), -0.03*ones(numel(x_s), 1)]; % 重建面: z 0.01 m, 比测量面更靠近声源 [x_r, y_r] meshgrid(-0.18:0.01:0.18, -0.18:0.01:0.18); recon_pos [x_r(:), y_r(:), 0.01*ones(numel(x_r), 1)]; f 2000; c 343; k 2*pi*f/c; % p_meas 由实测或仿真得到 [p_recon, q_est] esm_nah_core(p_meas, src_pos, meas_pos, recon_pos, k, 1e-6);这个示例展示了最常用的三层结构测量面在中间、等效源面偏向声源一侧、重建面在测量面与声源之间。等效源面与声源表面之间留出距离是为了避免 (\mathbf{G}) 矩阵出现近奇异。等效源间距比测量间距大一点可以降低矩阵条件数代价是重建空间分辨率略下降。参数alpha在这里取 (10^{-6}) 只是演示值。它跟声压幅值的绝对值、阵列尺寸、频率都有关系正确值需要通过第 3 章的正则化参数选择方法来确定不能复制粘贴。3. 正则化参数与阵列网格设计等效源法重建精度的两个旋钮3.1 奇异值谱告诉你的三件事先对 (\mathbf{G}) 做奇异值分解代码是[U,S,V] svd(G, econ)。观察diag(S)的衰减曲线可以得到三个信息奇异值下降越快问题病态越严重需要更大的正则化参数奇异值谱出现明显平台说明网格间距太小或等效源面与测量面距离过大前几个奇异值对应的左奇异向量代表声场的整体辐射成分尾部向量则对应倏逝波和噪声敏感方向。工程上我会用semilogy(diag(S), o)把奇异值画出来。如果曲线在前 10 个奇异值就衰减了三个数量级直接把alpha设在最大奇异值的 1e-3 到 1e-2 倍开始试。不要依赖pinv(G, tol)里的固定容差因为tol不能反映声压数据的信噪比。3.2 L 曲线与 GCV两种实用的正则化参数选择函数最稳妥的选择方式是绘制 L 曲线横轴是残差范数 (|\mathbf{G}\mathbf{q}-\mathbf{p}|)纵轴是解范数 (|\mathbf{q}|)一组候选alpha会画出一条 L 形曲线拐角处就是噪声放大和过度平滑的平衡点。function alpha_opt lcurve_alpha(G, p_meas, alpha_list) % 用 L 曲线拐点选择 Tikhonov 正则化参数 % 输入: % G - MxN 传递矩阵 % p_meas - Mx1 测量声压 % alpha_list - 候选正则化参数向量建议对数等间隔排列 % 输出: % alpha_opt - 选出的正则化参数 [U, S, V] svd(G, econ); s diag(S); d U * p_meas; n_alpha numel(alpha_list); eta zeros(n_alpha, 1); % 解范数 rho zeros(n_alpha, 1); % 残差范数 for ii 1:n_alpha a2 alpha_list(ii)^2; q_proj (s ./ (s.^2 a2)) .* d; % 奇异值坐标系下的解 eta(ii) norm(q_proj); rho(ii) norm(d - s .* q_proj); end % 找折线上离首尾连线最远的点作为拐角 dx rho - rho(1); dy eta - eta(end); seg_len sqrt((rho(end)-rho(1))^2 (eta(end)-eta(1))^2); dist abs(dy*(rho(end)-rho(1)) - dx*(eta(end)-eta(1))) / seg_len; [~, idx] max(dist); alpha_opt alpha_list(idx); end这段代码在每次循环里做了标准的 Tikhonov 滤波s / (s^2 alpha^2)是奇异值域的滤波系数d是声压投影到左奇异向量后的坐标。eta和rho分别在衡量解的能量和残差的能量。拐点检测用了点到直线的距离公式代码短且不依赖额外工具箱适合直接嵌入处理流程。GCV广义交叉验证是另一种常用方法它对噪声模型更鲁棒但计算量略大。当测量面点数超过 1000 且希望全自动选参时我优先用 GCV当数据里有明显的少数坏点时L 曲线更可控。两者的alpha结果不应该差太多如果差了一个数量级先检查数据里是否有个别测量通道损坏。3.3 网格间距、测量距离与最高重建频率的工程边界等效源法近场声全息对阵列几何的要求比平面 NAH 宽松但它不是无条件的。表里列的是我在汽车 NVH 台架上常用的起始参数参数推荐范围说明阵元间距 (\Delta l)(\lambda_{min}/3) 到 (\lambda_{min}/2)超过半波长会出现重建伪峰测量面与声源距离小于 (\lambda_{min}/2)尽量小于 5 cm距离越大倏逝波衰减越严重等效源面与声源距离0.05λ 到 0.2λ太近矩阵近奇异太远丢失近场细节等效源间距测量间距的 1 到 2 倍大于 3 倍会降低重建分辨率测量面相对声源的扩展每边外扩 3 个 (\Delta l) 以上降低边界截断效应判断最高重建频率的经验方法是观察奇异值谱改变频率重新算 (\mathbf{G}) 奇异值当第 10 个奇异值小于最大奇异值的 1e-4 时这个频率附近的近场信息已经淹没在噪声容限里。与其硬提频率不如缩小测量距离或加密等效源面。很多初学者把阵元间距当成最高频率的唯一限制实际上测量面到声源的距离在高频段更致命。测量距离增大一倍倏逝波幅值按指数衰减重建信噪比损失远大于加密网格带来的收益。这也是近场声全息必须贴脸测的原因。4. 等效源法在实测数据中的进阶处理局部测量、稀疏源与多源分离4.1 局部测量面的边界截断与加窗处理实测中很难让测量面完全覆盖声源的整个辐射区域。特别是大型箱体结构测量面覆盖一半时重建场边缘会出现明显的半圆形条纹这是截断效应不是真实声源。常见的做法是对测量声压做二维加窗。汉宁窗可以把边缘不连续压下去代价是重建区只取窗口中间部分。我的处理顺序是先用零填充扩展测量面到原来的 1.5 倍尺寸然后加二维汉宁窗重建完成后只保留中心区域结果。% 测量面原始尺寸 Mx x My win_2d hann(Mx) * hann(My); % 二维汉宁窗 p_meas_windowed p_meas_2d .* win_2d; % 先加窗再扩展 % 扩展测量面: 四周补零 Mx_pad round(Mx*1.5); My_pad round(My*1.5); p_ext zeros(Mx_pad, My_pad); x0 floor((Mx_pad-Mx)/2) 1; y0 floor((My_pad-My)/2) 1; p_ext(x0:x0Mx-1, y0:y0My-1) p_meas_windowed;这段代码把窗函数应用与零填充分开。注意hann(Mx) * hann(My)生成的是 Mx×My 矩阵与p_meas_2d形状一致。扩展后的重构结果在中心区域以外基本不可信画云图时要把边缘裁掉。二维汉宁窗的旁瓣抑制能力够用不需要用更高阶的窗否则会过度展宽主瓣。4.2 稀疏等效源重建声源数量少时的压缩感知方案旋转机械的噪声问题里真正的主要声源往往只有三到五个。这种情况下还在规则网格上布置几百个等效源只会让解在非声源位置上出现虚假小值。改用稀疏重建的思路在候选等效源网格上用 OMP正交匹配追踪只选少数源位置参与声场拟合。function q omp_esm(G, p_meas, sparsity) % 用正交匹配追踪求稀疏等效源源强 % G: MxN 传递矩阵, p_meas: Mx1 测量声压, sparsity: 期望稀疏度 [M, N] size(G); G_norm sqrt(sum(abs(G).^2, 1)); % 各列能量用于归一化 G G ./ G_norm; % 避免靠近测量面的源优先被选中 r p_meas; support []; q_sup []; for iter 1:sparsity corr G * r; % 残差与各候选源的相关性 [~, idx] max(abs(corr)); % 找出最相关的源 support union(support, idx); % 加入支撑集 q_sup G(:, support) \ p_meas; % 用当前支撑集做最小二乘 r p_meas - G(:, support) * q_sup; % 更新残差 end q zeros(N, 1); q(support) q_sup; % 未入选位置源强为 0 q q ./ G_norm; % 还原因列归一化造成的幅值缩放 end这个 OMP 实现没有做第二步正交化迭代但对中等大小的传递矩阵足够稳定。G_norm的归一化很重要格林函数幅值随距离衰减如果直接做相关距离测量面近的候选源天然更容易被选中导致所有匹配都集中在一个空间区域。还原幅值时特别注意q q ./ G_norm的顺序少这一步重建声压正常但源强幅值完全失真。稀疏度参数的选择不宜超过真实声源数量的两倍。如果不知道声源数量可以先跑一遍普通 Tikhonov 重建把重建面上声压级最高的几个局部峰值数量作为sparsity的参考。4.3 多源分离与频带分组多个声源在空间上重叠时单一等效源面很难同时描述两个不同辐射特性的声源。我一般会把候选源分两层一层贴近主声源表面另一层放在次声源附近。反演时两层源被同一个 (\mathbf{G}) 矩阵约束重建完成后按源强幅值分布判断主次源贡献。频带处理上高频段需要更靠近声源的等效源面低频段则可以把源面外移减小条件数。如果整段 200 Hz 到 5000 Hz 用一个源面位置低频重建结果会被高频设置拖累表现为重建声场在低频段出现“空洞”。更合理的做法是分三段设置源面距离低频段 0.2λ、中频段 0.1λ、高频段 0.05λ分别重建后再把频响曲线拼接起来。AI 辅助编码工具在写这类分段处理脚本时容易把 Python 的切片语义带进 MATLAB最常见的是把1:end写成0:end-1。等效源位置矩阵的索引错误不会直接报错但重建云图上会出现一排规则的伪源。用这类工具生成代码后逐行检查坐标矩阵行数比检查算法逻辑更紧迫。5. 用仿真数据闭环验证等效源法重建的四个检查点5.1 先造一个已知声源再让算法自己重建把代码跑在实测数据上之前用仿真数据做一次闭环验证能省下大量排错时间。做法是放置几个已知位置和源强的单极子作为真声源正向生成测量面声压加噪后用等效源法重建最后把重建声场与真声源在重建面上产生的声压对比。% 已知真实源位置与源强 true_src [0.1, 0, -0.02; -0.08, 0.05, -0.02]; q_true [1; -0.6]; % 用这些源在测量面生成声压 p_meas green_matrix(meas_pos, true_src, k) * q_true; noise_level 0.02 * max(abs(p_meas)); % 2% 幅值噪声 p_meas_noisy p_meas noise_level * (randn(size(p_meas)) 1j*randn(size(p_meas))) / sqrt(2); % 重建 alpha_opt lcurve_alpha(G, p_meas_noisy, logspace(-8, 0, 40)); [p_recon, q_est] esm_nah_core(p_meas_noisy, src_pos, meas_pos, recon_pos, k, alpha_opt);加噪时要生成复噪声实部虚部独立且除以sqrt(2)保持总功率恒定。用重建完毕的q_est重新计算重建面声压并与真实值对比如果归一化误差小于 10%说明整条链路正常。误差更大时逐段排查先关掉噪声看误差是否来源于模型失配再加回噪声看是否来源于正则化过平滑。5.2 三个容易误判的重建失效现象与调整方向现象可能原因调整手段重建云图出现径向放射条纹正则化参数过小等效源间距较小导致解振荡增大alpha或让等效源面离声源远 20%高频区重建声压整体偏弱测量面距离声源过远倏逝波已衰减到噪声底缩小测量距离后再测或接受分辨率上限低频区云图中心一个圆斑等效源数量不足声场被少数源平均化加密低频段的等效源网格或检查源面距离放射状条纹在线性阵列中尤其明显因为 y 方向的阵列孔径不足会让 G 矩阵沿该方向奇异值衰减更快。判断条纹来源的一个技巧是改变alpha一个量级看条纹是否随之变化变化明显则是正则化问题不变则是网格配置问题。5.3 用残差分布定位阵列坏点重建完成后计算残差 (\mathbf{r} \mathbf{p}{meas} - \mathbf{G}\mathbf{q}{est})把残差按测量点坐标画成云图。正常数据残差是随机分布的如果残差集中在某个测量通道周围说明该通道存在相位偏差或增益异常。逐个剔除坏点后重新标定阵列比在正则化里加权重更可靠。把残差云图保存成png导出时注意imagesc会翻转 y 轴与plot3的坐标方向不一致分析前先统一坐标轴方向习惯避免把坏点位置判断错。本文还有配套的精品资源点击获取