简介这份资源面向声学成像、近场声全息与声场重建方向的学习者与研究人员聚焦声干涉获取全息图、相位全息图解析及声场复振幅重建这一完整技术链路适合具备信号处理与MATLAB基础、希望动手复现声全息算法的中高级读者。压缩包共2个文件以1个m脚本和1个txt说明文档为主包体约2KB脚本对应相位全息图逆傅里叶变换等重建流程文本则梳理声全息与声干涉的理论要点便于对照理解算法实现。目前已有542人学习下载可作为入门近场声全息的轻量实践素材。读者可借此掌握从声波振幅与相位信息采集、相位全息图频率域转换到空间域声场重构的关键思路并延伸至无损检测、声学成像、噪声控制等应用场景为后续搭建麦克风阵列实验或编写重建程序提供参考。1. 从一张相位全息图到可听声场nah.zip 到底在重建什么如果你手头只有一只麦克风阵列采到的一帧复声压数据想把它变成「任意位置都能算出声压」的三维声场最直接的路子就是近场声全息NAH。而nah.zip_全息 重建_声场_声干涉_相位全息图_重建声场这个标题把整条链路的关键词都摆出来了全息、声干涉、相位全息图、声场重建。它讲的不是某个玄学算法而是一套从「记录干涉条纹」到「反演声场」的工程流程——用相位全息图把声波的相位信息编码下来再通过重建算法把目标面上的声压、振速还原出来。这套东西解决的核心问题是你没法在声源表面贴满传感器但你又想知道声源表面的声场分布。NAH 的价值就在于「由外推内」——在离声源一定距离的测量面上采数据反推回声源面。适合谁做扬声器阵列测试、电机异响定位、水下声呐成像、MEMS 声学器件表征的工程师以及想把声全息从论文搬到代码里的研究生。相位全息图在这里不是光学里的那张底片而是声场复振幅的相位分布记录声干涉则是它形成的物理机制。2. 声干涉与相位全息图重建声场的物理底座2.1 声干涉条纹是怎么被记录成相位全息图的声全息和光学全息共享同一个底层逻辑参考波与物波干涉把相位信息转成强度信息记录下来。在声学里参考波通常是一路已知相位的同频信号物波是声源辐射后到达全息面的声压。两路叠加后全息面上的复声压分布就带上了干涉条纹条纹的疏密和走向直接对应物波的相位梯度。工程上记录相位全息图有两种常见做法。第一种是扫描法用一只参考麦克风固定在已知位置另一只测量麦克风在全息面上逐点扫描每点同时采参考通道和测量通道做互谱得到该点的复声压。第二种是阵列法用麦克风阵列一次性采完整个全息面参考通道用阵列中某一只或外部参考。扫描法设备便宜但慢阵列法快但对通道一致性要求高。记录下来的复声压写成极坐标形式就是幅度加相位相位那一项就是相位全息图。注意这里说的「全息图」不是一张图片文件而是一个复数矩阵实部虚部或者幅度相位都行。很多初学者以为要存成图像格式其实没必要.npy或.mat存复数矩阵最省事。提示参考通道的信噪比直接决定相位全息图的质量。参考麦克风要放在声源近场但不要被声源直接冲击的位置否则参考本身被污染整张相位图都会歪。2.2 为什么重建必须回到波数域做角谱外推有了全息面上的复声压下一步是反演。NAH 最经典的算法是波数域角谱法K-space 或 angular spectrum method。它的思路是把全息面上的声压做二维傅里叶变换从空间域变到波数域然后在波数域乘上一个传播算子把场从全息面外推到目标面再逆变换回空间域。传播算子长这样对于向声源方向反推evanescent 波放大传播因子是exp(j*kz*d)其中kz sqrt(k^2 - kx^2 - ky^2)d是外推距离。当kx^2 ky^2 k^2时kz变成虚数传播因子变成指数衰减或指数增长。反推时这些倏逝波被指数放大这就是 NAH 分辨率能突破半波长限制的原因也是它最脆弱的地方——噪声也被一起放大。所以重建流程不是简单乘个算子就完事必须配正则化。常见做法是加 Tikhonov 正则化或者截断奇异值分解把那些放大倍数过高的波数分量砍掉。这一步做不好重建出来的声场就是一片雪花。2.3 从复声压矩阵到目标面声场的最小可跑流程下面这段 Python 代码演示了从全息面复声压到目标面重建的最小闭环。假设你已经有了全息面复声压矩阵p_holo形状(Ny, Nx)以及采样间距dx, dy声速c频率f外推距离d。import numpy as np def nah_reconstruct(p_holo, dx, dy, c, f, d, reg_alpha1e-3): p_holo: 全息面复声压, shape (Ny, Nx) dx, dy: 采样间距 (m) c: 声速 (m/s) f: 频率 (Hz) d: 外推距离 (m), 正值表示向声源方向反推 reg_alpha: Tikhonov 正则化系数 返回: 目标面复声压, shape (Ny, Nx) Ny, Nx p_holo.shape k 2 * np.pi * f / c # 波数域坐标 kx 2 * np.pi * np.fft.fftfreq(Nx, ddx) ky 2 * np.pi * np.fft.fftfreq(Ny, ddy) KX, KY np.meshgrid(kx, ky) # 传播波数 kz kz_sq k**2 - KX**2 - KY**2 kz np.sqrt(kz_sq.astype(complex)) # 反推传播算子, 倏逝波用指数衰减/增长描述 propagator np.exp(1j * kz * d) # Tikhonov 正则化: 抑制倏逝波放大 propagator_reg propagator / (1 reg_alpha * np.abs(propagator)**2) # 波数域重建 P_holo np.fft.fft2(p_holo) P_target P_holo * propagator_reg p_target np.fft.ifft2(P_target) return p_target逻辑说明先把全息面声压做二维 FFT 变到波数域然后对每个波数分量乘上反推传播因子。kz用复数开方保证倏逝波区域自动进入指数分支。正则化项reg_alpha控制放大抑制强度值越大抑制越狠但分辨率越低。参数说明dx, dy必须和实际扫描或阵列间距一致填错会导致波数轴整体偏移重建结果完全错位。d的符号约定要统一这里正值表示从全息面向声源方向反推。reg_alpha一般从1e-3起调噪声大就加大到1e-2但别超过1e-1否则传播波也被压掉。f是分析频率窄带分析时取峰值频率宽带分析要逐频点循环。跑完这一步p_target就是目标面上的复声压。你可以取np.abs(p_target)看幅度分布取np.angle(p_target)看相位分布。如果目标面就是声源面那这个结果就是声源表面的声压重建。3. 相位全息图重建声场的参数怎么设五个必调项3.1 全息面采样间距与阵列孔径的匹配采样间距dx, dy决定了波数域的奈奎斯特上限k_nyq pi/dx。如果声波频率对应的波数k接近或超过k_nyq就会出现空间混叠重建结果里出现虚假声源。经验规则是每个波长至少采 4 到 6 个点也就是dx lambda/4。对于 1 kHz 的空气声波长约 0.34 mdx要小于 8.5 cm。对于 20 kHz波长 1.7 cmdx要小于 4 mm这时候就得用密集阵列或者扫描。阵列孔径L Nx * dx决定了波数域的分辨率dk 2*pi/L。孔径越大波数分辨率越高能区分的声源细节越多。但孔径大了边缘效应也明显通常要做加窗或者补零。补零能改善波数域插值但不能真正提高分辨率。注意采样间距和孔径是一对矛盾。间距小意味着阵元多孔径大意味着阵列物理尺寸大。实际选型时先确定最高分析频率反推dx再根据需要的空间分辨率确定Nx。3.2 正则化系数对倏逝波放大的抑制边界正则化系数reg_alpha是 NAH 里最需要手调的参数。它的作用是在传播算子上加一个阻尼项形式可以写成H_reg H / (1 alpha * |H|^2)其中H是原始传播算子。当|H|很大倏逝波区域H_reg趋近于1/(alpha * |H|)放大被限制住。alpha太小倏逝波放大没被压住重建结果里高频噪声爆炸。alpha太大倏逝波信息被过度压制分辨率退化到和普通波束形成差不多。判断标准是看重建声压的动态范围如果目标面上出现大量幅度远超全息面最大值的点说明alpha太小如果重建结果过于平滑、细节全丢说明alpha太大。一个实用的调参方法是画 L 曲线横轴是正则化残差纵轴是解的范数拐点对应的alpha就是比较平衡的值。工程上懒得画曲线的话从1e-3开始每次乘 3 或除以 3看重建结果什么时候从「雪花」变成「能看」。3.3 外推距离的符号约定与数值稳定性外推距离d的符号是最容易翻车的地方。不同文献的约定不一样有的定义d为从声源到全息面的距离反推时用exp(-j*kz*d)有的定义d为从全息面到目标面的距离反推时用exp(j*kz*d)。代码里必须统一否则重建出来的声场要么发散要么完全对不上。数值稳定性方面d越大倏逝波的放大倍数越高。反推距离超过 2 到 3 个波长后正则化很难同时保住分辨率和稳定性。所以 NAH 的实用反推距离一般控制在 1 到 2 个波长以内。如果声源面离全息面太远建议先用波束形成做声源定位再用 NAH 做局部精细重建。3.4 频率分辨率与窄带分析的取舍NAH 是窄带算法每个频点单独重建。如果声源是宽带的需要对每个频点循环执行重建然后叠加或者分别看。频率分辨率df由采样时长决定df 1/T。df太大峰值频率可能落在两个频点之间重建幅度偏低df太小每个频点的数据信噪比下降相位全息图质量变差。常见做法是先用 Welch 法做互谱取峰值频率附近的几个频点分别重建看哪个频点的重建结果最干净。对于旋转机械的阶次噪声直接锁定阶次频率df取转速的 1/10 左右就够了。3.5 参考通道相位对齐的校准步骤相位全息图的质量极度依赖参考通道和测量通道之间的相位一致性。多通道采集系统里每个通道的模拟前端都有不同的相位延迟不校准的话相位全息图里会混入通道间的相位偏差重建出来的声源位置会整体偏移。校准步骤用一个已知相位的信号源同时馈入所有通道采一段数据算每个通道相对于参考通道的互谱相位得到相位修正因子。正式测量时每个测量通道的复声压乘以对应的相位修正因子。这一步在阵列法里尤其重要扫描法里如果参考和测量用同一套前端通道间偏差可以忽略但参考麦克风的物理位置必须精确记录。4. 避坑与排查相位全息图重建声场时最容易翻车的五件事4.1 重建结果出现镜像声源现象目标面上除了真实声源位置对称位置还出现一个虚假声源幅度和真实声源接近。原因二维 FFT 默认把空间域当成周期延拓全息面边缘不连续时波数域出现高频泄漏逆变换后表现为镜像。另外如果全息面数据没有去均值零波数分量过强也会在重建面上产生背景。解决做 FFT 之前先对全息面复声压加二维 Tukey 窗或者 Hann 窗窗宽覆盖有效数据区域。同时去掉全息面的均值。如果镜像仍然存在检查dx, dy是否填错间距填错会导致波数轴镜像。4.2 倏逝波放大导致高频雪花现象重建声压幅度分布里出现大量细密的高频斑点幅度远超全息面最大值看起来像雪花屏。原因反推传播算子在倏逝波区域是指数增长噪声被同等放大。正则化系数太小或者根本没加正则化。解决加大reg_alpha从1e-3逐步加到1e-2甚至1e-1。同时检查全息面数据的信噪比如果全息面本身噪声就大先做空间平均或者多次平均降噪。另一个办法是限制波数域的有效范围把kx^2 ky^2 (k_max)^2的分量直接置零k_max取1.2*k到1.5*k。4.3 重建声源位置整体偏移现象重建出来的声源位置和实际声源位置差了一个固定距离方向一致。原因参考通道和测量通道之间的相位偏差没有校准或者外推距离d的符号约定和实际测量几何不一致。解决先做通道相位校准用已知相位信号测出每个通道的相位修正因子。然后核对d的定义如果全息面在声源面前方z0处反推距离应该是z0符号根据传播算子约定确定。建议用一个已知位置的单极子声源做验证确认重建位置和实际位置重合后再测实际对象。4.4 低频重建结果模糊不清现象高频重建结果还行低频重建出来的声源又大又模糊分辨率很差。原因低频波长长倏逝波衰减慢但传播波的波数分辨率dk由阵列孔径决定。孔径不够大时低频的波数分辨率不足以区分声源细节。另外低频段全息面数据的信噪比通常较低相位估计误差大。解决低频重建需要更大的阵列孔径。如果孔径固定只能接受较低的空间分辨率或者改用近场声全息的统计优化方法如 SONAH它在小孔径下比波数域法表现更好。检查低频段互谱的相干函数相干低于 0.8 的频点直接丢弃。4.5 多频点重建结果无法叠加现象每个频点单独重建都正常但把多个频点的重建声压直接相加后总声场完全不对。原因不同频点的复声压相位基准不一致。每个频点单独做 FFT 和重建时相位参考是各自频点的互谱相位频点之间没有统一的绝对相位。直接相加会互相抵消。解决如果要合成总声场必须在互谱阶段保留统一的参考通道相位或者对每个频点的重建结果取幅度后按能量叠加而不是复数相加。对于声压级合成用sqrt(sum(|p_i|^2))更合理。如果确实需要复数叠加确保所有频点共享同一个时间零点。5. 用单极子仿真验证重建链路一个可复现的闭环5.1 构造单极子声源的全息面复声压验证 NAH 链路最可靠的办法是自己造一个已知答案。单极子声源在自由场中的声压解析解是p(r) A * exp(-j*k*r) / r其中r是场点到声源的距离。用这个解析解生成全息面复声压然后跑重建看能不能把声源位置和幅度还原出来。import numpy as np def monopole_pressure(x, y, z, xs, ys, zs, k, A1.0): 单极子声源在自由场中的复声压 r np.sqrt((x - xs)**2 (y - ys)**2 (z - zs)**2) r np.maximum(r, 1e-6) # 避免除零 return A * np.exp(-1j * k * r) / r # 参数设置 c 343.0 # 声速 m/s f 2000.0 # 频率 Hz k 2 * np.pi * f / c lam c / f # 波长 # 全息面网格 Nx, Ny 64, 64 dx lam / 6 # 采样间距, 每波长6点 x np.arange(Nx) * dx y np.arange(Ny) * dx X, Y np.meshgrid(x, y) # 声源位置: 全息面中心正下方 0.5 波长处 xs (Nx - 1) * dx / 2 ys (Ny - 1) * dx / 2 zs -0.5 * lam z_holo 0.0 # 生成全息面复声压 p_holo monopole_pressure(X, Y, z_holo, xs, ys, zs, k) # 加少量噪声模拟实测 noise_level 0.01 p_holo noise_level * np.std(np.abs(p_holo)) * (np.random.randn(Ny, Nx) 1j * np.random.randn(Ny, Nx))逻辑说明monopole_pressure直接实现单极子解析解X, Y是全息面网格zs是声源深度负值表示在全息面下方。加噪声是为了模拟实测中的信噪比噪声水平取信号标准差的 1%。参数说明dx lam/6满足每波长 6 点避免空间混叠。zs -0.5*lam表示声源离全息面半个波长反推距离就是0.5*lam。噪声水平0.01对应约 40 dB 信噪比比较接近实际测量条件。5.2 跑重建并对比解析解把p_holo喂给第 2 章的nah_reconstruct反推距离d 0.5*lam然后和声源面上的解析解对比。# 重建到声源面 d 0.5 * lam p_recon nah_reconstruct(p_holo, dx, dx, c, f, d, reg_alpha1e-3) # 声源面解析解 z_src zs p_exact monopole_pressure(X, Y, z_src, xs, ys, zs, k) # 对比幅度分布 amp_recon np.abs(p_recon) amp_exact np.abs(p_exact) # 归一化后比较 amp_recon_norm amp_recon / np.max(amp_recon) amp_exact_norm amp_exact / np.max(amp_exact) # 找峰值位置 idx_recon np.unravel_index(np.argmax(amp_recon_norm), amp_recon_norm.shape) idx_exact np.unravel_index(np.argmax(amp_exact_norm), amp_exact_norm.shape) print(f重建峰值位置: {idx_recon}, 解析峰值位置: {idx_exact}) print(f重建峰值幅度: {amp_recon_norm[idx_recon]:.4f}, 解析峰值幅度: {amp_exact_norm[idx_exact]:.4f}) # 计算相对误差 error np.linalg.norm(amp_recon_norm - amp_exact_norm) / np.linalg.norm(amp_exact_norm) print(f归一化幅度相对误差: {error:.4f})逻辑说明重建结果和解析解都做归一化后比较避免绝对幅度标定的影响。峰值位置对比验证声源定位精度相对误差验证整体重建质量。参数说明reg_alpha1e-3是起点如果误差偏大可以调整。d必须和zs的绝对值一致。如果重建峰值位置和解析峰值位置差超过一个网格说明相位校准或者d的符号有问题。5.3 从误差曲线反推正则化系数单极子仿真的另一个用途是标定reg_alpha。固定噪声水平扫一遍reg_alpha画相对误差曲线误差最低点就是当前信噪比下的最优正则化系数。alphas np.logspace(-4, -1, 20) errors [] for alpha in alphas: p_r nah_reconstruct(p_holo, dx, dx, c, f, d, reg_alphaalpha) amp_r np.abs(p_r) amp_r_norm amp_r / np.max(amp_r) err np.linalg.norm(amp_r_norm - amp_exact_norm) / np.linalg.norm(amp_exact_norm) errors.append(err) best_alpha alphas[np.argmin(errors)] print(f最优正则化系数: {best_alpha:.2e}, 对应误差: {min(errors):.4f})逻辑说明np.logspace(-4, -1, 20)在四个数量级上取 20 个点覆盖从欠正则化到过正则化的范围。每个alpha跑一次重建算归一化误差取最小值对应的alpha。参数说明这个标定只对当前噪声水平和频率有效。换频率或者换信噪比后需要重新标定。实际测量中如果不知道噪声水平可以用全息面数据的相干函数估计相干越低噪声越大alpha相应加大。5.4 把验证链路固化成回归测试单极子仿真最大的价值是它可以变成一个回归测试。每次改了重建代码跑一遍单极子看误差有没有变大。误差突然跳变说明改动引入了 bug。我一般会把monopole_pressure、nah_reconstruct和误差计算封装成一个测试函数参数化频率和噪声水平跑一组用例。频率覆盖 500 Hz 到 10 kHz噪声水平覆盖 1% 到 10%。所有用例的误差都在阈值以内才认为重建链路是健康的。这个习惯帮我省了很多后悔药。有一次改波数域坐标生成方式单极子测试立刻报误差从 0.02 跳到 0.3查了半天发现是fftfreq的d参数单位从米错写成毫米。没有这个测试这种错误要到实际测量数据上才会暴露那时候排查成本就高多了。希望帮到你。本文还有配套的精品资源点击获取 近场声全息NAH：相位全息图与声干涉重建声场实战
简介这份资源面向声学成像、近场声全息与声场重建方向的学习者与研究人员聚焦声干涉获取全息图、相位全息图解析及声场复振幅重建这一完整技术链路适合具备信号处理与MATLAB基础、希望动手复现声全息算法的中高级读者。压缩包共2个文件以1个m脚本和1个txt说明文档为主包体约2KB脚本对应相位全息图逆傅里叶变换等重建流程文本则梳理声全息与声干涉的理论要点便于对照理解算法实现。目前已有542人学习下载可作为入门近场声全息的轻量实践素材。读者可借此掌握从声波振幅与相位信息采集、相位全息图频率域转换到空间域声场重构的关键思路并延伸至无损检测、声学成像、噪声控制等应用场景为后续搭建麦克风阵列实验或编写重建程序提供参考。1. 从一张相位全息图到可听声场nah.zip 到底在重建什么如果你手头只有一只麦克风阵列采到的一帧复声压数据想把它变成「任意位置都能算出声压」的三维声场最直接的路子就是近场声全息NAH。而nah.zip_全息 重建_声场_声干涉_相位全息图_重建声场这个标题把整条链路的关键词都摆出来了全息、声干涉、相位全息图、声场重建。它讲的不是某个玄学算法而是一套从「记录干涉条纹」到「反演声场」的工程流程——用相位全息图把声波的相位信息编码下来再通过重建算法把目标面上的声压、振速还原出来。这套东西解决的核心问题是你没法在声源表面贴满传感器但你又想知道声源表面的声场分布。NAH 的价值就在于「由外推内」——在离声源一定距离的测量面上采数据反推回声源面。适合谁做扬声器阵列测试、电机异响定位、水下声呐成像、MEMS 声学器件表征的工程师以及想把声全息从论文搬到代码里的研究生。相位全息图在这里不是光学里的那张底片而是声场复振幅的相位分布记录声干涉则是它形成的物理机制。2. 声干涉与相位全息图重建声场的物理底座2.1 声干涉条纹是怎么被记录成相位全息图的声全息和光学全息共享同一个底层逻辑参考波与物波干涉把相位信息转成强度信息记录下来。在声学里参考波通常是一路已知相位的同频信号物波是声源辐射后到达全息面的声压。两路叠加后全息面上的复声压分布就带上了干涉条纹条纹的疏密和走向直接对应物波的相位梯度。工程上记录相位全息图有两种常见做法。第一种是扫描法用一只参考麦克风固定在已知位置另一只测量麦克风在全息面上逐点扫描每点同时采参考通道和测量通道做互谱得到该点的复声压。第二种是阵列法用麦克风阵列一次性采完整个全息面参考通道用阵列中某一只或外部参考。扫描法设备便宜但慢阵列法快但对通道一致性要求高。记录下来的复声压写成极坐标形式就是幅度加相位相位那一项就是相位全息图。注意这里说的「全息图」不是一张图片文件而是一个复数矩阵实部虚部或者幅度相位都行。很多初学者以为要存成图像格式其实没必要.npy或.mat存复数矩阵最省事。提示参考通道的信噪比直接决定相位全息图的质量。参考麦克风要放在声源近场但不要被声源直接冲击的位置否则参考本身被污染整张相位图都会歪。2.2 为什么重建必须回到波数域做角谱外推有了全息面上的复声压下一步是反演。NAH 最经典的算法是波数域角谱法K-space 或 angular spectrum method。它的思路是把全息面上的声压做二维傅里叶变换从空间域变到波数域然后在波数域乘上一个传播算子把场从全息面外推到目标面再逆变换回空间域。传播算子长这样对于向声源方向反推evanescent 波放大传播因子是exp(j*kz*d)其中kz sqrt(k^2 - kx^2 - ky^2)d是外推距离。当kx^2 ky^2 k^2时kz变成虚数传播因子变成指数衰减或指数增长。反推时这些倏逝波被指数放大这就是 NAH 分辨率能突破半波长限制的原因也是它最脆弱的地方——噪声也被一起放大。所以重建流程不是简单乘个算子就完事必须配正则化。常见做法是加 Tikhonov 正则化或者截断奇异值分解把那些放大倍数过高的波数分量砍掉。这一步做不好重建出来的声场就是一片雪花。2.3 从复声压矩阵到目标面声场的最小可跑流程下面这段 Python 代码演示了从全息面复声压到目标面重建的最小闭环。假设你已经有了全息面复声压矩阵p_holo形状(Ny, Nx)以及采样间距dx, dy声速c频率f外推距离d。import numpy as np def nah_reconstruct(p_holo, dx, dy, c, f, d, reg_alpha1e-3): p_holo: 全息面复声压, shape (Ny, Nx) dx, dy: 采样间距 (m) c: 声速 (m/s) f: 频率 (Hz) d: 外推距离 (m), 正值表示向声源方向反推 reg_alpha: Tikhonov 正则化系数 返回: 目标面复声压, shape (Ny, Nx) Ny, Nx p_holo.shape k 2 * np.pi * f / c # 波数域坐标 kx 2 * np.pi * np.fft.fftfreq(Nx, ddx) ky 2 * np.pi * np.fft.fftfreq(Ny, ddy) KX, KY np.meshgrid(kx, ky) # 传播波数 kz kz_sq k**2 - KX**2 - KY**2 kz np.sqrt(kz_sq.astype(complex)) # 反推传播算子, 倏逝波用指数衰减/增长描述 propagator np.exp(1j * kz * d) # Tikhonov 正则化: 抑制倏逝波放大 propagator_reg propagator / (1 reg_alpha * np.abs(propagator)**2) # 波数域重建 P_holo np.fft.fft2(p_holo) P_target P_holo * propagator_reg p_target np.fft.ifft2(P_target) return p_target逻辑说明先把全息面声压做二维 FFT 变到波数域然后对每个波数分量乘上反推传播因子。kz用复数开方保证倏逝波区域自动进入指数分支。正则化项reg_alpha控制放大抑制强度值越大抑制越狠但分辨率越低。参数说明dx, dy必须和实际扫描或阵列间距一致填错会导致波数轴整体偏移重建结果完全错位。d的符号约定要统一这里正值表示从全息面向声源方向反推。reg_alpha一般从1e-3起调噪声大就加大到1e-2但别超过1e-1否则传播波也被压掉。f是分析频率窄带分析时取峰值频率宽带分析要逐频点循环。跑完这一步p_target就是目标面上的复声压。你可以取np.abs(p_target)看幅度分布取np.angle(p_target)看相位分布。如果目标面就是声源面那这个结果就是声源表面的声压重建。3. 相位全息图重建声场的参数怎么设五个必调项3.1 全息面采样间距与阵列孔径的匹配采样间距dx, dy决定了波数域的奈奎斯特上限k_nyq pi/dx。如果声波频率对应的波数k接近或超过k_nyq就会出现空间混叠重建结果里出现虚假声源。经验规则是每个波长至少采 4 到 6 个点也就是dx lambda/4。对于 1 kHz 的空气声波长约 0.34 mdx要小于 8.5 cm。对于 20 kHz波长 1.7 cmdx要小于 4 mm这时候就得用密集阵列或者扫描。阵列孔径L Nx * dx决定了波数域的分辨率dk 2*pi/L。孔径越大波数分辨率越高能区分的声源细节越多。但孔径大了边缘效应也明显通常要做加窗或者补零。补零能改善波数域插值但不能真正提高分辨率。注意采样间距和孔径是一对矛盾。间距小意味着阵元多孔径大意味着阵列物理尺寸大。实际选型时先确定最高分析频率反推dx再根据需要的空间分辨率确定Nx。3.2 正则化系数对倏逝波放大的抑制边界正则化系数reg_alpha是 NAH 里最需要手调的参数。它的作用是在传播算子上加一个阻尼项形式可以写成H_reg H / (1 alpha * |H|^2)其中H是原始传播算子。当|H|很大倏逝波区域H_reg趋近于1/(alpha * |H|)放大被限制住。alpha太小倏逝波放大没被压住重建结果里高频噪声爆炸。alpha太大倏逝波信息被过度压制分辨率退化到和普通波束形成差不多。判断标准是看重建声压的动态范围如果目标面上出现大量幅度远超全息面最大值的点说明alpha太小如果重建结果过于平滑、细节全丢说明alpha太大。一个实用的调参方法是画 L 曲线横轴是正则化残差纵轴是解的范数拐点对应的alpha就是比较平衡的值。工程上懒得画曲线的话从1e-3开始每次乘 3 或除以 3看重建结果什么时候从「雪花」变成「能看」。3.3 外推距离的符号约定与数值稳定性外推距离d的符号是最容易翻车的地方。不同文献的约定不一样有的定义d为从声源到全息面的距离反推时用exp(-j*kz*d)有的定义d为从全息面到目标面的距离反推时用exp(j*kz*d)。代码里必须统一否则重建出来的声场要么发散要么完全对不上。数值稳定性方面d越大倏逝波的放大倍数越高。反推距离超过 2 到 3 个波长后正则化很难同时保住分辨率和稳定性。所以 NAH 的实用反推距离一般控制在 1 到 2 个波长以内。如果声源面离全息面太远建议先用波束形成做声源定位再用 NAH 做局部精细重建。3.4 频率分辨率与窄带分析的取舍NAH 是窄带算法每个频点单独重建。如果声源是宽带的需要对每个频点循环执行重建然后叠加或者分别看。频率分辨率df由采样时长决定df 1/T。df太大峰值频率可能落在两个频点之间重建幅度偏低df太小每个频点的数据信噪比下降相位全息图质量变差。常见做法是先用 Welch 法做互谱取峰值频率附近的几个频点分别重建看哪个频点的重建结果最干净。对于旋转机械的阶次噪声直接锁定阶次频率df取转速的 1/10 左右就够了。3.5 参考通道相位对齐的校准步骤相位全息图的质量极度依赖参考通道和测量通道之间的相位一致性。多通道采集系统里每个通道的模拟前端都有不同的相位延迟不校准的话相位全息图里会混入通道间的相位偏差重建出来的声源位置会整体偏移。校准步骤用一个已知相位的信号源同时馈入所有通道采一段数据算每个通道相对于参考通道的互谱相位得到相位修正因子。正式测量时每个测量通道的复声压乘以对应的相位修正因子。这一步在阵列法里尤其重要扫描法里如果参考和测量用同一套前端通道间偏差可以忽略但参考麦克风的物理位置必须精确记录。4. 避坑与排查相位全息图重建声场时最容易翻车的五件事4.1 重建结果出现镜像声源现象目标面上除了真实声源位置对称位置还出现一个虚假声源幅度和真实声源接近。原因二维 FFT 默认把空间域当成周期延拓全息面边缘不连续时波数域出现高频泄漏逆变换后表现为镜像。另外如果全息面数据没有去均值零波数分量过强也会在重建面上产生背景。解决做 FFT 之前先对全息面复声压加二维 Tukey 窗或者 Hann 窗窗宽覆盖有效数据区域。同时去掉全息面的均值。如果镜像仍然存在检查dx, dy是否填错间距填错会导致波数轴镜像。4.2 倏逝波放大导致高频雪花现象重建声压幅度分布里出现大量细密的高频斑点幅度远超全息面最大值看起来像雪花屏。原因反推传播算子在倏逝波区域是指数增长噪声被同等放大。正则化系数太小或者根本没加正则化。解决加大reg_alpha从1e-3逐步加到1e-2甚至1e-1。同时检查全息面数据的信噪比如果全息面本身噪声就大先做空间平均或者多次平均降噪。另一个办法是限制波数域的有效范围把kx^2 ky^2 (k_max)^2的分量直接置零k_max取1.2*k到1.5*k。4.3 重建声源位置整体偏移现象重建出来的声源位置和实际声源位置差了一个固定距离方向一致。原因参考通道和测量通道之间的相位偏差没有校准或者外推距离d的符号约定和实际测量几何不一致。解决先做通道相位校准用已知相位信号测出每个通道的相位修正因子。然后核对d的定义如果全息面在声源面前方z0处反推距离应该是z0符号根据传播算子约定确定。建议用一个已知位置的单极子声源做验证确认重建位置和实际位置重合后再测实际对象。4.4 低频重建结果模糊不清现象高频重建结果还行低频重建出来的声源又大又模糊分辨率很差。原因低频波长长倏逝波衰减慢但传播波的波数分辨率dk由阵列孔径决定。孔径不够大时低频的波数分辨率不足以区分声源细节。另外低频段全息面数据的信噪比通常较低相位估计误差大。解决低频重建需要更大的阵列孔径。如果孔径固定只能接受较低的空间分辨率或者改用近场声全息的统计优化方法如 SONAH它在小孔径下比波数域法表现更好。检查低频段互谱的相干函数相干低于 0.8 的频点直接丢弃。4.5 多频点重建结果无法叠加现象每个频点单独重建都正常但把多个频点的重建声压直接相加后总声场完全不对。原因不同频点的复声压相位基准不一致。每个频点单独做 FFT 和重建时相位参考是各自频点的互谱相位频点之间没有统一的绝对相位。直接相加会互相抵消。解决如果要合成总声场必须在互谱阶段保留统一的参考通道相位或者对每个频点的重建结果取幅度后按能量叠加而不是复数相加。对于声压级合成用sqrt(sum(|p_i|^2))更合理。如果确实需要复数叠加确保所有频点共享同一个时间零点。5. 用单极子仿真验证重建链路一个可复现的闭环5.1 构造单极子声源的全息面复声压验证 NAH 链路最可靠的办法是自己造一个已知答案。单极子声源在自由场中的声压解析解是p(r) A * exp(-j*k*r) / r其中r是场点到声源的距离。用这个解析解生成全息面复声压然后跑重建看能不能把声源位置和幅度还原出来。import numpy as np def monopole_pressure(x, y, z, xs, ys, zs, k, A1.0): 单极子声源在自由场中的复声压 r np.sqrt((x - xs)**2 (y - ys)**2 (z - zs)**2) r np.maximum(r, 1e-6) # 避免除零 return A * np.exp(-1j * k * r) / r # 参数设置 c 343.0 # 声速 m/s f 2000.0 # 频率 Hz k 2 * np.pi * f / c lam c / f # 波长 # 全息面网格 Nx, Ny 64, 64 dx lam / 6 # 采样间距, 每波长6点 x np.arange(Nx) * dx y np.arange(Ny) * dx X, Y np.meshgrid(x, y) # 声源位置: 全息面中心正下方 0.5 波长处 xs (Nx - 1) * dx / 2 ys (Ny - 1) * dx / 2 zs -0.5 * lam z_holo 0.0 # 生成全息面复声压 p_holo monopole_pressure(X, Y, z_holo, xs, ys, zs, k) # 加少量噪声模拟实测 noise_level 0.01 p_holo noise_level * np.std(np.abs(p_holo)) * (np.random.randn(Ny, Nx) 1j * np.random.randn(Ny, Nx))逻辑说明monopole_pressure直接实现单极子解析解X, Y是全息面网格zs是声源深度负值表示在全息面下方。加噪声是为了模拟实测中的信噪比噪声水平取信号标准差的 1%。参数说明dx lam/6满足每波长 6 点避免空间混叠。zs -0.5*lam表示声源离全息面半个波长反推距离就是0.5*lam。噪声水平0.01对应约 40 dB 信噪比比较接近实际测量条件。5.2 跑重建并对比解析解把p_holo喂给第 2 章的nah_reconstruct反推距离d 0.5*lam然后和声源面上的解析解对比。# 重建到声源面 d 0.5 * lam p_recon nah_reconstruct(p_holo, dx, dx, c, f, d, reg_alpha1e-3) # 声源面解析解 z_src zs p_exact monopole_pressure(X, Y, z_src, xs, ys, zs, k) # 对比幅度分布 amp_recon np.abs(p_recon) amp_exact np.abs(p_exact) # 归一化后比较 amp_recon_norm amp_recon / np.max(amp_recon) amp_exact_norm amp_exact / np.max(amp_exact) # 找峰值位置 idx_recon np.unravel_index(np.argmax(amp_recon_norm), amp_recon_norm.shape) idx_exact np.unravel_index(np.argmax(amp_exact_norm), amp_exact_norm.shape) print(f重建峰值位置: {idx_recon}, 解析峰值位置: {idx_exact}) print(f重建峰值幅度: {amp_recon_norm[idx_recon]:.4f}, 解析峰值幅度: {amp_exact_norm[idx_exact]:.4f}) # 计算相对误差 error np.linalg.norm(amp_recon_norm - amp_exact_norm) / np.linalg.norm(amp_exact_norm) print(f归一化幅度相对误差: {error:.4f})逻辑说明重建结果和解析解都做归一化后比较避免绝对幅度标定的影响。峰值位置对比验证声源定位精度相对误差验证整体重建质量。参数说明reg_alpha1e-3是起点如果误差偏大可以调整。d必须和zs的绝对值一致。如果重建峰值位置和解析峰值位置差超过一个网格说明相位校准或者d的符号有问题。5.3 从误差曲线反推正则化系数单极子仿真的另一个用途是标定reg_alpha。固定噪声水平扫一遍reg_alpha画相对误差曲线误差最低点就是当前信噪比下的最优正则化系数。alphas np.logspace(-4, -1, 20) errors [] for alpha in alphas: p_r nah_reconstruct(p_holo, dx, dx, c, f, d, reg_alphaalpha) amp_r np.abs(p_r) amp_r_norm amp_r / np.max(amp_r) err np.linalg.norm(amp_r_norm - amp_exact_norm) / np.linalg.norm(amp_exact_norm) errors.append(err) best_alpha alphas[np.argmin(errors)] print(f最优正则化系数: {best_alpha:.2e}, 对应误差: {min(errors):.4f})逻辑说明np.logspace(-4, -1, 20)在四个数量级上取 20 个点覆盖从欠正则化到过正则化的范围。每个alpha跑一次重建算归一化误差取最小值对应的alpha。参数说明这个标定只对当前噪声水平和频率有效。换频率或者换信噪比后需要重新标定。实际测量中如果不知道噪声水平可以用全息面数据的相干函数估计相干越低噪声越大alpha相应加大。5.4 把验证链路固化成回归测试单极子仿真最大的价值是它可以变成一个回归测试。每次改了重建代码跑一遍单极子看误差有没有变大。误差突然跳变说明改动引入了 bug。我一般会把monopole_pressure、nah_reconstruct和误差计算封装成一个测试函数参数化频率和噪声水平跑一组用例。频率覆盖 500 Hz 到 10 kHz噪声水平覆盖 1% 到 10%。所有用例的误差都在阈值以内才认为重建链路是健康的。这个习惯帮我省了很多后悔药。有一次改波数域坐标生成方式单极子测试立刻报误差从 0.02 跳到 0.3查了半天发现是fftfreq的d参数单位从米错写成毫米。没有这个测试这种错误要到实际测量数据上才会暴露那时候排查成本就高多了。希望帮到你。本文还有配套的精品资源点击获取