简介:这份资源面向声学成像、近场声全息(NAH)方向的学习者与研究人员,聚焦声干涉获取声全息图、相位全息图解析与声场重建这一完整技术链路,适合具备一定信号处理与MATLAB基础、希望动手复现声场重建流程的读者。压缩包共2个文件,包含1个m脚本与1个txt说明文档,整体仅2KB,体积轻量,便于快速下载与本地运行。其中脚本对应相位全息图经逆傅里叶变换重构空间域声场的核心计算过程,文本文件则用于补充声全息原理与实验背景,二者配合可帮助读者理解振幅与相位信息如何共同决定复振幅分布。目前已有542人学习下载,说明该方向具备一定关注度。资源虽小,但覆盖了从声干涉记录到数字重建的关键环节,可作为无损检测、声学成像、噪声控制等场景的入门实践素材,也便于在此基础上扩展麦克风阵列数据处理与算法验证。
1. 从 nah.zip 说起:声全息重建到底在算什么
如果你手头有一个叫nah.zip的压缩包,里面大概率躺着近场声全息(Near-field Acoustic Holography,NAH)的实测数据、相位全息图或者一套重建脚本。这个方向的核心问题很具体:用麦克风阵列在靠近声源的一个平面上测到声压,怎么反推出声源表面的声场分布,或者预测另一个平面上的声场。声干涉和相位全息图是绕不开的两个关键词——前者决定了你测到的数据里哪些是真实辐射、哪些是倏逝波在捣乱,后者决定了你重建时相位参考从哪来。
我最早接触这套东西是为了定位一块 PCB 上几个电容的异常振动。远场测出来一团糊,近场全息重建之后能直接看到声源面上的热点。适合做这件事的人很明确:做 NVH 的、做声源定位的、做扬声器阵列校准的,以及被“声场重建”四个字吸引过来但不知道从哪下手的人。下面按我实际跑通的路径,从数据长什么样一路讲到参数怎么调、坑在哪。
2. 近场声全息的数据长什么样:相位全息图与声压复矩阵
2.1 全息面复声压矩阵的构成
NAH 的输入不是一张图,而是一个复数矩阵。假设你用 K 个麦克风组成一个平面阵列,在距离声源面z_h的全息面上扫描 M 个点,每个点测到的是复声压p(x, y, z_h),包含幅值和相位。这个矩阵的维度通常是M × K或者按扫描网格重排成Ny × Nx的二维复数数组。
相位全息图这个名字容易让人误解。它不是光学里那种干涉条纹照片,而是指全息面上的相位分布∠p(x, y, z_h)。如果你用的是参考麦克风法,相位是相对于参考通道的;如果用声光调制或者相位步进干涉,相位来自多帧相移。声干涉在这里的作用是:全息面上的声压是声源各点辐射的球面波叠加结果,干涉条纹的疏密直接对应声源频率和阵列孔径的关系。
一个典型的全息面复声压矩阵在 Python 里长这样:
import numpy as np # 假设扫描网格 32x32,测量频率 2000 Hz Nx, Ny = 32, 32 freq = 2000.0 c0 = 343.0 # 声速 m/s k = 2 * np.pi * freq / c0 # 波数 # 全息面坐标(单位:米),假设间距 0.02 m dx = dy = 0.02 x = np.arange(Nx) * dx y = np.arange(Ny) * dy X, Y = np.meshgrid(x, y) # 模拟一个点声源在全息面上产生的复声压 # 声源位于 (0.32, 0.32, 0),全息面 z_h = 0.05 m zs = 0.0 zh = 0.05 xs, ys = 0.32, 0.32 r = np.sqrt((X - xs)**2 + (Y - ys)**2 + (zh - zs)**2) p_hologram = np.exp(-1j * k * r) / r # 复声压,含幅值和相位 print("全息面复声压矩阵形状:", p_hologram.shape) print("相位范围: {:.2f} 到 {:.2f} rad".format(np.angle(p_hologram).min(), np.angle(p_hologram).max()))这段代码构造了一个理想点声源在全息面上的复声压。实际测量中p_hologram来自采集设备,但维度、物理意义完全一致。关键参数是k和zh:k决定倏逝波衰减速度,zh决定你离声源多近。一般要求zh < λ/2,否则倏逝波衰减到噪声里,重建分辨率会崩。
2.2 为什么必须用近场:倏逝波与分辨率的关系
远场测量只能拿到传播波,波数分量满足k_r < k。近场测量的价值在于捕捉到k_r > k的倏逝波分量,这些分量携带亚波长信息,是超分辨率重建的物理基础。但倏逝波随距离指数衰减,衰减因子是exp(-sqrt(k_r^2 - k^2) * zh)。这意味着zh每增加一点点,高频空间分量就掉一个数量级。
我一般会先算一下最关心的空间频率对应的衰减量。如果衰减超过 60 dB,基本可以认为这个分量被噪声淹没了,重建时强行放大只会得到一堆虚假热点。这一步没有代码,但可以用一个简单表格判断:
| 空间频率分量 | 波数比 k_r/k | 衰减因子 (zh=0.02m, f=2kHz) | 是否可用 |
|---|---|---|---|
| 传播波 | 0.5 | 1.0 | 是 |
| 传播波 | 0.9 | 1.0 | 是 |
| 倏逝波 | 1.2 | 0.12 | 勉强 |
| 倏逝波 | 2.0 | 0.001 | 否 |
| 倏逝波 | 3.0 | 1e-6 | 否 |
这张表解释了一个常见翻车场景:有人拿zh = 0.1 m的数据做 5 kHz 重建,结果全是噪声。不是算法不行,是物理上那些高频分量根本没传到全息面。
3. 从全息面反推声源面:NAH 重建的三种实现路径
3.1 空间傅里叶变换法:最直接但边界最敏感
空间傅里叶变换法(SFT-NAH)的思路很干净:对全息面复声压做二维 FFT,得到波数域谱P(kx, ky, zh),然后乘以反向传播算子exp(j * kz * (zh - zs)),再逆变换回空间域。其中kz = sqrt(k^2 - kx^2 - ky^2),当kx^2 + ky^2 > k^2时kz是虚数,对应倏逝波的指数衰减/放大。
import numpy as np def nah_reconstruct_sft(p_hologram, dx, dy, freq, zh, zs, c0=343.0): """ 空间傅里叶变换法近场声全息重建 p_hologram: 全息面复声压矩阵 (Ny, Nx) dx, dy: 网格间距 (m) freq: 频率 (Hz) zh: 全息面 z 坐标 (m) zs: 重建面 z 坐标 (m) """ Ny, Nx = p_hologram.shape k = 2 * np.pi * freq / c0 # 波数域坐标 kx = 2 * np.pi * np.fft.fftfreq(Nx, d=dx) ky = 2 * np.pi * np.fft.fftfreq(Ny, d=dy) KX, KY = np.meshgrid(kx, ky) # 传播算子 kz_sq = k**2 - KX**2 - KY**2 kz = np.sqrt(kz_sq.astype(complex)) # 反向传播:从 zh 到 zs # 注意:zs < zh 时,传播因子为 exp(1j * kz * (zh - zs)) # 倏逝波对应 kz 为虚数,exp(1j * 1j * |kz| * dz) = exp(-|kz| * dz),是放大 dz = zh - zs propagator = np.exp(1j * kz * dz) # 正则化:限制倏逝波放大倍数 max_gain = 100.0 # 最大放大倍数,根据信噪比调整 gain = np.abs(propagator) propagator[gain > max_gain] *= max_gain / gain[gain > max_gain] # 波数域滤波与重建 P_hologram = np.fft.fft2(p_hologram) P_source = P_hologram * propagator p_source = np.fft.ifft2(P_source) return p_source # 使用示例 p_source = nah_reconstruct_sft(p_hologram, dx, dy, freq, zh, zs=0.0) print("重建面声压幅值范围: {:.4f} 到 {:.4f}".format(np.abs(p_source).min(), np.abs(p_source).max()))这段代码里最关键的是max_gain这个正则化参数。理论上反向传播对倏逝波是无限放大的,但实测数据有噪声,不限制放大倍数的话,重建结果会被噪声主导。我一般从 100 开始试,如果重建面出现明显不合理的孤立尖峰,就降到 30 或 10。另一个参数是dz = zh - zs,重建面越靠近声源面,倏逝波放大越剧烈,对正则化越敏感。
SFT 法的边界问题很突出:FFT 默认数据是周期的,但实际测量阵列有限,边界处声压不连续,会在波数域产生泄漏。常见做法是加 Tukey 窗或者补零。补零能改善表观分辨率,但不增加真实信息,补零到 2 倍尺寸通常够用。
3.2 边界元法:适合任意形状声源但计算量大
当声源不是平面或者重建面形状不规则时,SFT 法就不够用了。边界元法(BEM-NAH)把声源表面离散成网格,建立全息面声压与表面声压的传递矩阵,然后求逆。传递矩阵的维度是M × N,M 是测量点数,N 是表面网格节点数。通常 M < N,所以是欠定问题,需要 Tikhonov 正则化或者 L1 稀疏约束。
import numpy as np from scipy.linalg import pinv def nah_reconstruct_bem(G, p_hologram, lambda_reg=1e-3): """ 边界元法重建(简化版) G: 传递矩阵 (M, N),由格林函数计算 p_hologram: 全息面复声压 (M,) lambda_reg: Tikhonov 正则化参数 """ # Tikhonov 正则化:min ||G p - p_h||^2 + lambda ||p||^2 # 解为 p = (G^H G + lambda I)^(-1) G^H p_h M, N = G.shape GH = G.conj().T A = GH @ G + lambda_reg * np.eye(N) b = GH @ p_hologram p_source = np.linalg.solve(A, b) return p_source # 传递矩阵构造示例(自由场格林函数) def build_transfer_matrix(source_nodes, hologram_nodes, freq, c0=343.0): k = 2 * np.pi * freq / c0 M = hologram_nodes.shape[0] N = source_nodes.shape[0] G = np.zeros((M, N), dtype=complex) for i in range(M): for j in range(N): r = np.linalg.norm(hologram_nodes[i] - source_nodes[j]) if r < 1e-12: r = 1e-12 G[i, j] = np.exp(-1j * k * r) / (4 * np.pi * r) return GBEM 法的参数核心是lambda_reg。太小则解不稳定,太大则过度平滑。我一般用 L 曲线法找拐点:把||G p - p_h||和||p||都算出来,画在双对数坐标上,拐点对应的lambda_reg就是比较平衡的值。BEM 的另一个坑是传递矩阵的条件数。如果全息面离声源太远,矩阵接近奇异,求逆会炸。这时候要么拉近测量距离,要么增加测量点数。
3.3 等效源法:在工程中最容易落地
等效源法(ESM-NAH)是我最常用的。它不直接求表面声压,而是在声源内部布置一堆等效点源,用这些点源的辐射来拟合全息面测量值,然后拿拟合好的等效源去预测任意位置的声场。好处是传递矩阵构造简单,不需要处理边界积分奇异点,而且等效源位置可以灵活布置。
import numpy as np def nah_reconstruct_esm(p_hologram, hologram_coords, equiv_coords, freq, c0=343.0, lambda_reg=1e-4): """ 等效源法近场声全息重建 p_hologram: 全息面复声压 (M,) hologram_coords: 全息面坐标 (M, 3) equiv_coords: 等效源坐标 (N, 3) freq: 频率 lambda_reg: 正则化参数 """ k = 2 * np.pi * freq / c0 M = hologram_coords.shape[0] N = equiv_coords.shape[0] # 构造传递矩阵 G = np.zeros((M, N), dtype=complex) for i in range(M): for j in range(N): r = np.linalg.norm(hologram_coords[i] - equiv_coords[j]) G[i, j] = np.exp(-1j * k * r) / (4 * np.pi * r) # Tikhonov 正则化求解等效源强度 GH = G.conj().T A = GH @ G + lambda_reg * np.eye(N) b = GH @ p_hologram q = np.linalg.solve(A, b) return q, G # 预测重建面声场 def predict_field(q, G_pred): """q: 等效源强度, G_pred: 预测点与等效源之间的传递矩阵""" return G_pred @ qESM 的参数比 BEM 多一个:等效源的位置和数量。我一般把等效源布置在声源面后方0.5 * dx到1 * dx的位置,数量取测量点数的 1/2 到 1/3。等效源太靠近声源面会导致传递矩阵病态,太远则拟合精度下降。lambda_reg同样用 L 曲线选。
三种方法对比:
| 方法 | 适用场景 | 参数敏感度 | 计算量 | 我的使用频率 |
|---|---|---|---|---|
| SFT | 平面阵列、规则网格 | 中,边界窗和 max_gain | 低 | 高,快速验证 |
| BEM | 任意形状声源 | 高,lambda 和网格质量 | 高 | 中,复杂结构 |
| ESM | 一般工程问题 | 中,等效源位置和 lambda | 中 | 最高,落地首选 |
4. 相位全息图获取与声干涉处理的避坑清单
4.1 相位参考丢失导致重建面出现镜像声源
现象:重建出来的声源面出现两个对称的热点,一个是真的,一个是镜像。
原因:相位全息图的参考相位选错了。如果用双麦克风法,参考麦克风放在声源另一侧,参考信号本身包含了反向传播的相位,导致全息面相位多了一个共轭分量。
解决:参考麦克风必须放在声源和全息面之间,或者用声源上的加速度计作为参考。如果已经测完了,可以尝试对全息面相位取共轭再重建,看镜像是否消失。我一般会在测量前用一个小喇叭在已知位置验证相位参考方向。
4.2 倏逝波放大倍数设太大导致重建面全是噪点
现象:重建面声压幅值范围异常大,比如比全息面大 1000 倍,而且空间分布像随机噪声。
原因:max_gain或lambda_reg设得太松,噪声的高空间频率分量被当成倏逝波放大了。
解决:先看全息面相位图。如果相位在相邻点之间跳变超过 π,说明信噪比不够,那些高频分量不可信。把max_gain降到 10 到 30,或者对全息面先做低通滤波。我习惯先画全息面幅值图,如果幅值本身就有明显噪声纹理,重建前必须滤波。
4.3 阵列网格间距大于半波长导致空间混叠
现象:重建结果出现栅瓣,声源位置偏移,或者出现多个虚假声源。
原因:空间采样定理要求dx < λ/2。如果测量频率 4 kHz,λ ≈ 0.086 m,dx必须小于 0.043 m。很多扫描架的最小步进是 0.05 m,刚好不满足。
解决:要么提高扫描密度,要么降低分析频率上限。如果数据已经采了,可以在波数域加一个低通滤波器,把|k| > π/dx的分量砍掉,代价是损失部分高频信息。我一般会在测量前算好dx和最高频率的对应关系,避免白跑一趟。
4.4 全息面尺寸不够导致低频重建失败
现象:低频重建时声源面出现明显的边缘振荡,声源中心幅值也不对。
原因:全息面尺寸必须大于声源尺寸加上几个波长的余量。低频波长长,如果全息面只比声源大一点点,边缘的声压截断会在波数域产生严重泄漏。
解决:全息面每边至少比声源大λ/2。如果做不到,用补零加窗,但补零不能替代真实测量。我做过一个 500 Hz 的案例,声源 0.3 m,全息面只做了 0.4 m,重建结果完全不可用,后来扩到 0.8 m 才正常。
4.5 声干涉条纹被误判为声源分布
现象:重建面出现周期性条纹,看起来像多个声源,但实际只有一个。
原因:全息面上的声干涉条纹被反向传播算子放大后,在重建面形成了虚假的周期性结构。这不是噪声,是真实的干涉图样,但它对应的是全息面的测量位置,不是声源面。
解决:检查重建距离dz。如果dz太小,干涉条纹还没充分发散,就会在重建面保留。适当增大dz或者对波数域加一个传播波滤波,只保留k_r < k的分量,可以抑制这种假象。但注意,这样也会损失倏逝波信息,分辨率会下降。
5. 验证重建质量的三个硬指标与一个实战技巧
5.1 用全息面回代残差判断正则化是否合理
重建完成后,把等效源或者表面声压正向传播回全息面,和实测全息面声压比较。残差定义为||p_h - G q|| / ||p_h||。如果残差小于 5%,说明拟合没问题;如果大于 20%,要么正则化太重,要么传递矩阵模型不对。我一般会扫一遍lambda_reg,画残差曲线,选拐点。
def compute_residual(p_hologram, G, q): """计算全息面回代残差""" p_pred = G @ q residual = np.linalg.norm(p_hologram - p_pred) / np.linalg.norm(p_hologram) return residual # 扫描 lambda_reg lambdas = np.logspace(-6, -1, 20) residuals = [] for lam in lambdas: q, G = nah_reconstruct_esm(p_hologram, hologram_coords, equiv_coords, freq, lambda_reg=lam) residuals.append(compute_residual(p_hologram, G, q)) # 找拐点:残差开始快速上升的位置 # 实际使用时可以画图,这里打印几个关键值 for lam, res in zip(lambdas, residuals): print(f"lambda={lam:.2e}, residual={res:.4f}")5.2 重建面声压与直接测量对比
如果条件允许,在重建面位置再测一遍声压,和重建结果对比。这是最直接的验证。我一般会在声源面附近选几个点,用探针麦克风测,然后和重建值比幅值和相位。幅值误差在 2 dB 以内、相位误差在 20 度以内,就算合格。如果相位差很大,多半是相位参考或者传播算子符号搞反了。
5.3 空间分辨率用两个靠近的点声源测试
想知道你的 NAH 系统实际分辨率,可以放两个点声源,间距从λ/2开始逐渐缩小,看重建面能否分开。能分开的最小间距就是实际分辨率。我实测过,zh = 0.02 m、dx = 0.01 m、2 kHz 的条件下,分辨率大概在λ/4左右。如果间距小于这个值,重建面会合并成一个热点。
5.4 一个实战技巧:先做传播波重建,再加倏逝波
很多人一上来就全波数重建,结果被倏逝波噪声搞崩。我的习惯是分两步:第一步只保留传播波分量k_r < k做重建,得到一个稳定的低频结果;第二步逐步加入倏逝波分量,每加一档看重建面是否出现不合理尖峰。这样能清楚知道哪些空间频率是可信的,哪些是噪声。具体做法是在波数域加一个平滑过渡窗,而不是硬截断。
def apply_kfilter(P_k, kx, ky, k, k_cut_ratio=1.0, transition=0.1): """ 波数域滤波:保留 k_r < k_cut_ratio * k 的分量 transition: 过渡带宽度比例 """ KX, KY = np.meshgrid(kx, ky) kr = np.sqrt(KX**2 + KY**2) k_cut = k_cut_ratio * k # 平滑窗:1 在通带,0 在阻带,余弦过渡 H = np.ones_like(kr) H[kr > k_cut * (1 + transition)] = 0 mask = (kr > k_cut * (1 - transition)) & (kr <= k_cut * (1 + transition)) H[mask] = 0.5 * (1 + np.cos(np.pi * (kr[mask] - k_cut * (1 - transition)) / (2 * k_cut * transition))) return P_k * H这个滤波器的好处是不会在波数域产生振铃。k_cut_ratio从 1.0 开始,只保留传播波,然后逐步加到 1.5、2.0,观察重建面变化。如果加到 1.5 就开始出现噪点,说明你的信噪比只支持到 1.5 倍波数。这个上限直接决定了你的重建分辨率极限。
我做了这么多年声全息,最大的教训是:不要迷信算法,先看物理。全息面离声源多远、阵列间距多大、信噪比多少,这三个数基本决定了你能重建出什么。算法只是在这个边界内尽量逼近。每次拿到新数据,我第一件事是画全息面幅值和相位图,第二件事是算倏逝波衰减表,第三件事才是选重建方法。希望帮到你。
本文还有配套的精品资源,点击获取