简介:面向声全息与近场声成像研究的MATLAB实现包,围绕近场声全息算法展开,适用于声学检测、无损检测及噪声控制等场景的科研人员与工程师。压缩包共10个文件,含6个m脚本与4个txt说明文档,脚本包括主程序Main.m及f2_5.m、f1_9.m、shuju.m等子程序,txt文档则提供模块注释与补充说明,整体仅12KB,轻量便捷。目前已有270人学习下载。内容覆盖从麦克风阵列数据采集、预处理、干涉计算、成像重建到后处理的完整流程,通过运行示例可直观理解近场声全息算法如何将近场声压数据转化为可视化的声场图像。资源虽小巧,但各脚本结构清晰、参数易改,适合学习MATLAB声学仿真或快速搭建近场声全息原型系统,为实际工程问题提供可参照的实现思路。
1. 声全息程序包动手之前,先搞清它到底解决什么问题
声全息(Acoustic Holography)这名字听着玄,落到现场其实就一句话:用一块麦克风阵列贴近噪声源表面采集复声压场,再通过逆传播算法还原出声源面上的声压和振速分布,直接指出“噪声到底是从哪块区域辐射出来的”。你手里这个 holoGRAPHY_GitHub.rar 压缩包里,全息成像、声全息程序、声成像、近场成像四个词指向的是同一套技术路线,organized7ah 更像是整理者手动盖的版本标记,不代表功能分级,也不代表源码质量。
这套程序适合三类人:做噪声源识别的测试工程师,想知道电机端盖、齿轮箱、轮胎接地面到底是哪个局部在叫;做结构辐射声学仿真的研发,想把仿真频响和实测重建结果互相校核;以及一直在用远场波束形成、但被低频分辨率卡住喉咙、想换近场方案的技术人员。它值不值得跑,取决于你要不要分辨小于半个波长的声源细节——如果只是看个大致方位,远场波束形成更省事;如果要看到“哪一圈螺栓在漏声”,近场声全息是少数能干到亚波长分辨的实测手段。
下文从原理、最小复现、参数到验收,按一线落地顺序拆开讲。你不需要提前懂太多数学,但每个参数的物理含义我会交代清楚,因为这套算法的翻车点几乎全藏在参数的物理直觉里。
2. 近场声全息到底在算什么:从测量声压到声源面重建
2.1 为什么是近场而不是远场:倏逝波里才有亚波长细节
常规声学测量在远场做,麦克风收到的声压近似看成平面波叠加,声源的方向信息由各传声器之间的相位差体现。远场模型的极限是瑞利判据:分辨率被波长按阵列孔径和测量距离锁死。低频时一个 1 kHz 的声波波长 0.34 m,你要在 0.5 m 外分辨 0.2 m 尺度的两个声源,理论上就不可能——两个源在远场生成的波前几乎一样。
近场声全息换了个思路:把测量面贴到声源附近,离声源表面 0.1 到 1 个波长的距离。在这个区域里,声场除了传播波,还携带倏逝波(evanescent waves)。倏逝波的特点是波数分量大于自由场波数 k=2πf/c,沿声源表面传播时幅度随离面距离指数衰减,远场完全测不到,但在近场它还活着,而且恰恰是它携带了声源表面小于波长的结构信息。所以近场测量等于把“高分辨率信息”先抓住,再用算法从测量面数据里把声源面反算回去。
这也是为什么这套程序名字里“近场成像”跟在“声全息”后面。近场两个字不是测量习惯,是物理前提:离远了,信息已经衰减没了,再好的算法也补不回来。
2.2 从测量面角谱到声源面角谱:一步逆传播数学
近场声全息的经典实现是空间傅里叶变换法(角谱法),数学骨架三步走。设声源面为 z=0,测量面为 z=z_h,声源面上的声压 p_s(x,y,0) 做二维空间傅里叶变换,得到角谱 P_s(kx,ky)。声场向测量面正向传播时,角谱只发生相位变化:
P_h(kx,ky) = P_s(kx,ky) · exp(j·kz·z_h)
其中 kz = sqrt(k² - kx² - ky²),k=2πf/c。当 kx²+ky² > k² 时,kz 变成虚数,取负虚根,指数项变成 exp(-|kz|·z_h),这就是倏逝波的指数衰减。
反问题是把这个过程倒过来:把测量面角谱 P_h 乘上逆传播算子 exp(-j·kz·z_h),得到声源面角谱 P_s,再做逆傅里叶变换得到声源面声压。在倏逝波区域,这个逆算子会变成 exp(+|kz|·z_h),也就是把已经衰减掉的成分按指数放大回去——算法能不能用,关键就在于这步指数放大有没有把测量噪声一起放大成灾难。
除声压外,法向振速也能重建。由欧拉方程可得声压角谱与法向振速角谱的关系 V_s = P_s · kz / (ρ·c·k0),再做逆傅里叶变换即在声源面的法向振速分布。对结构辐射问题,振速云图比声压云图更能反映“哪块表面在振动发声”。
2.3 平面/柱面/球面坐标选型:rar 包默认平面的判断依据
声全息按坐标系统分平面、柱面、球面三种。平面 NAH 要求测量面与声源面平行,适合电机端盖、平板结构、大型板壳的辐射面重建;柱面 NAH 适合管道、轴流风扇外壁;球面 NAH 适合小型整机包络测量。这个 rar 包标题未说明坐标系,但带“近场成像”字样的开源实现绝大多数默认平面 NAH,原因很实际:平面角谱法只需要两次二维 FFT,计算量低,不依赖迭代,工程现场最快能落地。
如果你的被测对象是管道外壁或球状外壳,平面程序也能跑,但重建图会有几何畸变,这时优先找包里是否有多坐标入口。判断方法很简单:看数据读取函数里有没有 r、theta、z 圆柱参数,或者看示例数据文件是矩阵还是极坐标网格。没有的话就按平面处理,把管壁展开成矩形面测量,只在轴向上保留近场距离。
提示:平面 NAH 的反传播公式量纲和符号约定在不同代码里容易打架。下文代码统一采用“测量面 z 为正,声源面 z=0,kz 的倏逝波取负虚支”的约定,你拿到其他工程时先对齐符号,再对拍结果。3. 把 rar 里的声全息工程跑通:解压、数据格式与最小复现
3.1 解压前先做三件事:读说明、看目录结构、确认运行环境
从 GitHub 直接拿到 rar 打包的工程,比一帧一帧签出源码省事,不用操心子模块漏拉、文件缺失这类问题。但解压不能双击完就开跑,我一般先做三件事:第一,打开压缩包看路径——确认没有把一堆文件直接摊在根目录,也没有奇怪的长路径嵌套,这决定了解压后工程结构是否完整;第二,找 readme 或 version 说明,重点看运行环境,是全 MATLAB 工程还是 Python 为主,有没有注明依赖的 toolbox;第三,确认数据文件在不在包里,很多声全息程序把示例数据单独放在 data 目录,rar 打包时容易漏。
解压命令用 unrar 或 7-Zip 都可以,先在不解压的情况下看清单:
# 先看压缩包内容, 确认目录结构和有没有密码 unrar l hologaphy_GitHub.rar # 确认无误后解压到独立目录, 别解压到当前目录摊一地 unrar x hologaphy_GitHub.rar ./holography_src/ # 没有 unrar 时, 7z 也支持 rar 格式, 7z 的 -o 参数指定输出目录 7z x hologaphy_GitHub.rar -oholography_srcunrar l列出压缩包内文件清单,先看一层目录,能判断工程是源码加数据还是只有源码。unrar x保留包内完整目录结构解压,-o参数在 7z 里指定输出目录,养成解压到独立目录的习惯,避免脚本里的相对路径把文件写到别处。如果解压时提示密码,密码一般写在仓库介绍页或随包文本里,先找说明,别急着用恢复工具。
解压完成后,目录里大概率是三类东西:算法核心文件(.py或.m)、示例数据文件、说明文档。先跑说明文档里给的示例数据,再换自己的数据,这是最快的验证路径。
3.2 最小复现骨架:用 NumPy 把逆传播写成能跑的代码
不管 rar 里原工程用什么语言写的,NAH 的数学骨架就那几步。这里给一个最小 Python 实现,你把它和包里的算法做交叉验证,也能快速判断原代码的坐标系、符号约定和你手头数据是否一致:
import numpy as np from numpy.fft import fft2, ifft2, fftshift, ifftshift def nah_backprop(p_m, dx, f, z, c=343.0, rho=1.2): """平面近场声全息逆传播重建 p_m : 测量面复声压, 形状 (Ny, Nx), 单位 Pa dx : 网格步长, 单位 m, 要求 dx <= lambda/2 f : 分析频率, 单位 Hz z : 测量面到声源面的距离, 单位 m """ Ny, Nx = p_m.shape p_c = p_m - p_m.mean() # 去直流偏置 Pk = fftshift(fft2(p_c)) # 中心化角谱 # 角波数网格, 注意 fftfreq 默认给的是 cycle/unit, 要乘 2*pi kx = fftshift(np.fft.fftfreq(Nx, d=dx)) * 2 * np.pi ky = fftshift(np.fft.fftfreq(Ny, d=dx)) * 2 * np.pi KX, KY = np.meshgrid(kx, ky) kr2 = KX**2 + KY**2 k0 = 2 * np.pi * f / c # 轴向波数: 传播波取正实根, 倏逝波取负虚支 kz = np.zeros_like(kr2) kz[kr2 <= k0**2] = np.sqrt(k0**2 - kr2[kr2 <= k0**2]) kz[kr2 > k0**2] = -1j * np.sqrt(kr2[kr2 > k0**2] - k0**2) G = np.exp(-1j * kz * z) # 逆传播算子 Ps = Pk * G # 声源面角谱 p_s = np.real(ifft2(ifftshift(Ps))) # 声源面声压 Vs = Ps * kz / (rho * c * k0) # 声压角谱 -> 法向振速角谱 v_s = np.real(ifft2(ifftshift(Vs))) # 声源面法向振速 return p_s, v_s逻辑上,fft2做二维空间傅里叶变换,fftshift把零频移到矩阵中心,保证后续波数网格的 kx、ky 与角谱矩阵的排列一致。kz分两块计算,传播波区域是实数,对应常规相位传播;倏逝波区域取负虚支,让正向传播时指数衰减、反传播时指数放大,这是 NAH 区别于普通声场外推的核心。G是逆传播算子,乘到测量面角谱上再逆傅里叶变换,就得到声源面声压。
Vs的计算用到了法向振速与声压角谱的阻抗关系,注意空气密度rho和声速c要按实测环境微调——室温 20°C 时 c 取 343,10°C 时降到 337 左右。若把c差 10 个点,结果云图只会发生轻微尺度偏移,但dx超过半波长会直接出现空间混叠,云图出现重复副本,这是硬条件。
3.3 输入数据三种格式:.mat、.dat、.npy 的读取要点
声全息程序的数据入口各有不同,rar 包里常见三种:MATLAB 的.mat、文本/二进制的.dat、Python 的.npy。.mat文件用 scipy.io 读取,注意变量名在不同版本 MATLAB 里可能是p_meas、pressure或P_mic,读出来后先用shape确认维度顺序——矩阵是(Ny, Nx)还是(Nx, Ny),这直接影响重建图是否转置:
from scipy.io import loadmat mat = loadmat('data/measurement.mat') print([k for k in mat.keys() if not k.startswith('_')]) # 先看变量名 # 假设取到的是 p_mic, 统一reshape成 (Ny, Nx) 复声压 p_m = mat['p_mic'] if p_m.shape[0] < p_m.shape[1]: p_m = p_m.T # 把长边放到x方向, 和kx网格对齐.dat文件要区分是文本还是二进制:文本文件第一行一般是注释,用np.loadtxt加skiprows跳过;二进制分块存储时,常见格式是每帧数据前有一个文件头记录采样率和通道数,先用np.fromfile按块读取,再 reshape。.npy最简单,np.load直接读,但仍要检查数据类型是complex64还是float64。如果存的是实声压,说明这份数据只保留了幅值,相位已经丢了——遇到这种情况,NAH 程序跑出来也是白跑,原因放到第 5 章讲。
提示:接任何数据的第一步是画出测量面声压的幅值和相位两张图。幅值分布看阵列有没有坏通道,相位分布看声场是否连续。相位图如果是一团乱麻,先修采集,不要急着调算法。4. 四个必调参数:分析频率、截止波数、正则化与测量距离
4.1 分析频率与网格步长:先算半波长再动参数
NAH 对网格的要求比波束形成苛刻得多,空间采样定理要求网格步长 dx 至少小于半波长。这句话的实际含义是:你的测量面网格密度限定了可重建的最高频率。若 dx=0.05 m,半波长 0.05 m 对应频率约为 343/(2×0.05)=3430 Hz,超过这个频率的重建结果会出现空间混叠,云图上出现声源的“镜像副本”,看着像多个源,实则全是假的。
所以我拿到程序的第一件事不是调算法参数,而是用 dx 反推有效频段,再决定分析频率 f。工程现场常用的经验是让 dx 落在 λ/4 到 λ/2 之间,保证每个波长至少 4 个采样点。测量面孔径尺寸决定最低可分辨波数,孔径越小,低频空间分辨率越差,这和远场阵列的物理规律一致,近场改变了频率上限,改不了孔径的低频极限。
4.2 截止波数与正则化:一组能上手的初始参数表
逆传播在倏逝波区域是指数放大,测量噪声会被同等放大,所以必须做波数域滤波。最朴素的做法是硬截断:只保留 kx²+ky² ≤ (k_factor·k0)² 的分量,其余置零。k_factor 取 1.0 到 1.5 之间,超过 1.5 后放大噪声的风险急剧升高。
同样效果的另一种方法是 Tikhonov 正则化,在逆传播算子后面加一个抑制项,让倏逝波区域的放大倍数从指数级降为有界放大。给一组能直接上手的初始参数:
| 参数 | 初始值 | 调节依据 |
|---|---|---|
| 分析频率 f | ≤ 343/(2·dx) | 超过则出现空间混叠 |
| 截止波数 k_factor | 1.2 | 云图噪点过多时降为 1.05 |
| 正则化 alpha | 1e-4 | 重建面能量发散时增大到 1e-2 |
| 测量距离 z | ≤ λ/4 | 低频细节模糊时减小 z |
| 声速 c | 343 m/s | 按环境温度修正,每降 10°C 约减 6 m/s |
k_factor 和 alpha 是两个可以互换的“降压阀”:一个在波数域硬切,一个在逆算子上软抑制。实操经验是先用 k_factor=1.2、alpha=1e-4 跑通,再看重建云图的噪声水平——如果声源区域外全是细小颗粒状伪影,先降 k_factor;如果云图整体偏糊、峰值不锐利,再微增 alpha。这两个参数很少需要同时大改,一次动一个,否则无法判断是谁导致的失真。
4.3 窗函数与孔径外推:消除边缘 Gibbs 环的常见做法
测量面是有限孔径,空间傅里叶变换时边缘截断会产生 Gibbs 现象,表现为重建图边缘出现一整套平行波纹或环形伪影。最简单有效的缓解办法是空间域加窗,但加汉宁窗会连同孔径内的有效数据一起衰减,导致重建峰值变钝。工程上我更推荐镜像外推:在测量数据四周扩展 20% 的镜像副本,再做 FFT,重建后裁掉外扩区域。
# 四边做对称外扩, 降低孔径截断效应 n_ext = int(0.2 * p_m.shape[0]) # 外扩20%, 按行数估算 p_ext = np.pad(p_m, pad_width=n_ext, mode='symmetric') # 外扩后按新尺寸重新生成波数网格并重建 p_s_ext, v_s_ext = nah_backprop(p_ext, dx, f, z) p_s = p_s_ext[n_ext:-n_ext, n_ext:-n_ext] # 裁掉外扩, 保留原孔径symmetric模式把边缘数据镜像翻转进行延拓,等效于让 FFT 认为孔径更大,频谱泄漏明显降低。代价是计算量增大约 1.4 倍,对测量面 64×64 的阵列来说可忽略。注意外扩比例不宜超过 50%,延拓过多会改变有效孔径内的窗效应,反而让主瓣变宽。
4.4 测量距离 z:分辨率与信噪比的折中
z 是 NAH 里最需要现场经验的参数。z 越小,倏逝波携带的高频信息衰减越少,重建分辨率越高,但同时测量面会落入声源的近场反应区,传声器位置的声场对阵列定位误差极度敏感,几毫米的位置偏差就会让重建图局部扭曲。z 越大,测量越稳,但亚波长细节随指数衰减消失,重建结果退化成相当于波束形成的分辨率。
我一般把 z 控制在 0.1λ 到 1.0λ 之间。被测对象表面不平整时取上限,平整光滑表面可取下限。判断是否合适的土办法:把 z 设为两个值跑结果,比如 0.3λ 和 0.6λ,如果两幅云图的声源峰位置一致、只是锐度不同,说明 z 的选择在合理区间;如果峰的位置都变了,说明测量面离得太近,阵列局部散射已经污染了测量声场。
提示:NAH 的重建分辨率理论上可到波长的十分之一以下,但那是理想网格、理想信噪比下的结论。现场数据信噪比每下降 10 dB,有效分辨能力大约倒退一个档位,参数怎么调都救不回来。5. 声全息程序常见的五个坑:现象、原因与排查步骤
5.1 重建云图左右颠倒:坐标轴方向与网格顺序不匹配
现象:明明测量的是一个中心对称声源,重建图却是左右不对称,或者声源峰位置整体偏移几十厘米。原因:数据矩阵的行列顺序与波数网格 kx、ky 的定义没有对齐。常见于 MATLAB 转 Python 的工程,MATLAB 的矩阵按列存储,Python 的 NumPy 按行存储,数组转置后没有同步更新网格方向。解决:在跑真实数据前,先构造一个已知位置的点源解析解做全流程验证;左右颠倒就把数据矩阵转置后重跑,上下颠倒则是 ky 的符号取反,而不是盲目调滤波参数。
5.2 满屏高频噪点:倏逝波放大失控,先降截止波数
现象:重建云图在声源区域外全是密集的颗粒状亮斑,看起来像是噪声被“放大”了,而不是声源本身的形状。原因:逆传播对倏逝波区做了指数放大,而测量数据里的通道噪声、量化误差、阵元位置误差在这个区域同样被指数放大,信噪比低的频率完全被噪声主导。解决:第一步把 k_factor 降到 1.0 甚至 0.9,只保留传播波分量,确认噪点消失;若噪点仍在,检查各测量通道的幅值一致性,单个坏通道的数据会在重建图里形成一条亮线,这种坏通道问题调滤波参数永远解决不了,要回采集端补测或剔除该通道。
5.3 低频定位糊成一片:测量面离声源太远
现象:低频段(比如 800 Hz 以下)重建结果没有清晰的声源峰,功率云图像一个均匀发光的大包,完全看不出辐射细节。原因:z 已经大于几个波长,倏逝波在到达测量面前全部衰减完,测量面数据里只剩传播波信息,NAH 的理论优势不存在了,分辨率退化为远场阵列量级。解决:把测量面降到 z ≤ λ/2 以内重新测量;如果机械结构不允许贴太近,就放弃低频段的亚波长分辨诉求,改用波束形成做辅助定位,不要在算法参数上硬耗。
5.4 整幅图存在直流偏置:通道校准与均值处理
现象:重建云的背景颜色整体偏亮或偏暗,像蒙了一层雾,声源峰被背景淹没。原因:测量阵列各通道的直流偏置没有校准,或者数据里的人为偏置没有在 FFT 前去干净。代码里p_m - p_m.mean()只去掉了全局均值,如果各通道偏置不同,这一步处理不彻底。解决:采集前记录一段无声环境的数据,逐通道减去本底;再用代码里的去均值逻辑处理;处理后看频域零频分量是否为接近零的小值。零频分量就是直流,它会在逆传播后变成一个均匀分布在重建面上的常数背景。
5.5 有幅值没相位:实声压跑 NAH 必然失败
现象:程序能跑完,但重建云图没有像样的聚焦效果,怎么调参数都“糊”。原因:NAH 的输入必须是复声压,声源面的相位重建依赖测量面的相位信息;如果数据采集时只保存了各通道的幅值(比如用了普通声级计逐个位置测的 RMS 值),相位关系不存在,逆传播数学上就没有意义。解决:换多通道同步采集设备重新测;如果只有历史 RMS 数据,就别上 NAH,转用波束形成或声强法做定性分析。这是最容易被忽略、也最没有参数可救的坑,我在这上面浪费过一个下午。
6. 重建结果怎么验证:点源解析解与两个土办法
6.1 用点源解析解验证重建精度
跑通程序的标志不是能出图,而是出图可信。最可靠的验证不是拿别人的数据碰运气,而是构造一个自由场点源的理论声场生成仿真测量面,喂给程序看重建结果能不能还原点源位置。自由场点源声压有解析表达式:
p(R) = A · exp(j·k·R) / R
其中 R 是声源到测量点的距离。用这个解析解生成一个与真实数据同网格尺寸的复声压矩阵,跑一遍重建设定流程,重建图上应该在预设坐标位置出现一个尖锐单峰,云图背景应接近均匀、无环状伪影。这个验证同时校对了坐标方向、FFT 象限和符号约定——如果点源重建位置不对,先修坐标映射,再处理真实数据。
# 生成理论点源测量面: 网格 64x64, dx=0.05m, 源位于中心 nx, ny = 64, 64 x = (np.arange(nx) - nx/2) * dx y = (np.arange(ny) - ny/2) * dx X, Y = np.meshgrid(x, y) R = np.sqrt(X**2 + Y**2 + z**2) # z 为测量面到源面的距离 p_theory = 1.0 * np.exp(1j * k0 * R) / R # 单位幅值点源跑完后看重建峰位置。峰位置差一个网格以内,说明坐标映射和 FFT 方向全部正确;差得远则回头查网格生成顺序。
6.2 保留运行参数日志,让结果可追溯
声全息重建的每一步都对参数敏感,最恼火的情况是同一份数据跑出了两版明显不同的云图,却说不清是哪个参数改动导致的。我现在每次重建都会在输出目录落一份参数日志,连同云图一并存储,养成习惯后排查问题省一半时间:
# config.txt 与输出图放同一目录 cat > config.txt << EOF f = 1800 dx = 0.05 z = 0.08 k_factor = 1.20 alpha = 1e-4 c = 343 data = measurement_05.mat note = 电机端盖第2次测量, 阵列中心对准轴承座 EOF这份日志同时是数据交接的最小档案,换人接手项目时,有参数日志的数据可以直接复现,没有的直接作废重测。
做声全息这些年,我最深的体会是这算法有点“宁缺毋滥”:参数表能给的只是起点,真正的功力在于判断哪个参数不匹配现场物理条件。希望这套验证思路能帮你少走弯路,也希望帮到你——下次拿到类似的近场成像工程,先验相位、再对坐标、后调滤波,这个顺序别乱。
本文还有配套的精品资源,点击获取