近场声全息复现：从阵列数据到声源成像的完整流程

简介面向全息成像与近场声全息分析这套MATLAB程序包借助声波干涉实现对三维声场的记录与重建适宜声学检测、无损检测、近场声成像等方向的研究者、工程师和高年级学生参考。资源共10个文件其中6个.m源码脚本、4个.txt辅助文档压缩包整体仅12KB结构轻巧便于快速通读这些脚本覆盖数据读取、预处理、干涉相位计算、成像重建等关键环节txt文档用于存放数据或补充说明适合作为教学与练习样例。目前已有270人学习/下载是小而实用的声全息入门素材。脚本围绕近场声全息的典型流程展开主程序与分步函数齐全读者可对照数据采集、预处理、干涉计算、图像重建和后处理这几个关键阶段理解算法实现通过调整输入参数也能观察近场与远场声场分布的差异进而体会近场声全息在微小结构或复杂声环境中的优势。对需要从事声学检测、噪声控制或结构无损检测的读者这份程序包提供了可运行的MATLAB基础框架方便在此基础上扩展改造。1. 全息成像资料太多九成是 Demo先判断这份 organized7ah 值不值得复现从 GitHub 上扒下来的声学资源十有八九是论文附带代码解压密码、依赖缺失、数据残缺三个坑轮着踩。我拿到这份标题写着“全息成像、声全息程序、声成像、近场成像”的 rar 包时第一反应不是急着找 main.py而是先判断它是拿来演示的玩具还是能把近场声全息NAH从数据读取一路跑到出图的真程序。拆包后确认这是一份偏声成像方向的可复现项目organized7ah 是解压后的顶层目录里面包含了传声器阵列数据、波数域重建脚本和可视化模块。它解决的是声学工程师最常碰到的诉求拿一块阵列在近场扫一圈把时域压力变成某个频率下的声压幅值分布再推到声源表面识别噪声位置。适合两类人用一类是做噪声源定位、异音分析的工程师另一类是刚接触 NAH、想拿真实数据跑通全流程的学生。我按自己的习惯把压缩包重新整理了一遍把直接跑不通的地方都标了出来。下面先从 NAH 的原理和文件结构讲再给可复现的 Python 步骤最后是参数边界和踩坑清单。这张包里真正能复用的不是某一个函数而是整套“时域采样 → 频域提取 → 波数域传播 → 空间域成像”的流程。2. 近场声全息为什么能还原亚波长细节原理和文件结构先对齐2.1 近场测量与倏逝波声全息区别于波束成形的关键常规波束成形把声源假设成远场平面波测点离声源好几米远只能给出角度信息横向分辨率在低频段差得离谱。近场声全息把测量面挪到离声源很近的位置通常是十分之一到二分之一个波长量级。此时声场里还没有被空气衰减掉的倏逝波仍能被测到这就是 NAH 能还原小于一个波长的声源结构的原因。空气里 1 kHz 的波长大约是 343 mmλ/2 就是 171 mm。如果一个测量面贴在离声源 20 mm 的位置测到的声压里就同时含有传播波和倏逝波。传播波能一路传到远场倏逝波却在几个波长内快速衰减距离稍远就测不到。远场方法丢掉的信息恰恰是近场声全息最看重的部分。波数域里的传播关系可以写成一行判断逻辑KZ np.sqrt(k**2 - KX**2 - KY**2 0j)当 KX² KY² k² 时KZ 是实数对应传播波当 KX² KY² k² 时KZ 是纯虚数对应倏逝波。回推时乘以 exp(1j * KZ * z)虚数指数会把衰减掉的成分指数放大这也是 NAH 对噪声极度敏感的根本原因。理解这一点后面调截止波数就有依据了。所谓近场不是随便离近一点就行工程上可以按 z λ/2 做粗判实际项目里我一般控制在 λ/6 以内也就是 1 kHz 时测距不超过 55 mm。这份包里测量面与声源面的距离通常写在配置文件或脚本顶部复现前先把它找出来别等结果出来再猜。2.2 重建算法主流程从复声压到二维空间 FFT整个 NAH 重建可以拆成五步每一步都有对应的数据结构对阵列每个传声器做时间 FFT把时域压力变成某个频率下的复声压 P(x, y)。对复声压做二维空间 FFT得到波数域复声压 P(kx, ky)。在波数域里乘传播算子 exp(1j * KZ * Δz)把声场从测量面回推或外推到目标面。加低通或正则化压掉放大后的高频噪声。对处理后的波数谱做二维 IFFT得到目标面上的重建声压。这里核心的数学技巧是“空间卷积变波数域相乘”。测量面到重建面的传播在空间域是一个二维卷积直接用矩阵算64×64 的网格还能撑住到了 256×256 就会慢到让人怀疑人生在波数域里卷积变成逐点相乘速度完全是另一回事。这份 rar 包里如果是按传统思路写的脚本主体结构基本就是这几步。注意回推和外推的方向约定。从测量面向声源表面回推时传播距离用 z_h - z_s反过来从声源表面外推到测量面距离写 z_s - z_h。很多包会在注释里写清楚约定但现实中经常出现符号反了、结果整体镜像的情况后面避坑部分会专门讲。2.3 解压后先认识目录organized7ah 里代码和数据的约定GitHub 拉下来的资源第一件事不是急着装环境而是先把压缩包解开确认目录结构和数据格式。下面是我重新整理后习惯使用的解压方式# 解压到当前目录保持 organized7ah 这个顶层文件夹 unrar x hologaphy_GitHub.rar # 如果系统没装 unrar用 7-Zip 也可以 # 7z x hologaphy_GitHub.rar解压后按 unpack 的习惯先列一下顶层结构。这份包整理后的典型结构大致是这样organized7ah/ ├── README.md ├── main.py ├── nah_lib.py ├── config.yaml ├── data/ │ ├── holo_64x64_fs48000.csv │ ├── mic_positions.csv │ └── source_plane_meas.csv └── results/我并不保证你拿到手的文件名一字不差以包内 README 为准但这种布局是声全息资源里最常见的一类主脚本负责读数据和调流程库文件放重建函数data 目录放某次实测或仿真的阵列数据results 目录留给输出图像与结果数据。解压时有一个建议路径里不要带中文、空格和特殊符号。Windows 下解压到“C:\Users\你的用户名\Desktop\organized7ah”这类位置Python 脚本读相对路径时省掉很多坑。如果解压出来发现 data 目录为空八成是 GitHub 的 Git LFS 没拉全需要单独下载大文件跟代码本身没关系。3. 从原始阵列数据到二维声压图复现近场声全息的可执行步骤3.1 环境准备Python 版本、NumPy 和 SciPy 依赖复现这份材料不需要重型依赖NumPy、SciPy、Matplotlib 三件套足够。Python 版本建议 3.9 以上太老的版本对复数数组运算和文件路径处理都会有额外麻烦。我一般会新建一个虚拟环境避免跟系统 Python 打架python3 -m venv venv source venv/bin/activate pip install numpy scipy matplotlib pandaspandas 不是必须的但读 CSV 数据时它能少写好几行解析代码。如果这个包里的实验数据是 .mat 格式还需要补装 scipy 的 io 模块这个通常随 SciPy 一起装好了。环境装完后用下面的命令验证点关键依赖版本复数 FFT 的接口在较新版本里没有破坏性变化但老版本的 NumPy 在 fftfreq 返回类型上有差异。python -c import numpy, scipy, matplotlib; print(numpy.__version__, scipy.__version__)3.2 读取传声器阵列数据并提取目标频率的复声压阵列数据通常是一张二维表每一行是一个传声器通道每一列是时间采样点。读取后先不要急着做 FFT打印形状确认维度再决定怎么 reshape。下面这段代码把 CSV 数据读到内存并按 64×64 的阵列排布重组成三维数组import numpy as np import pandas as pd raw pd.read_csv(data/holo_64x64_fs48000.csv, headerNone).values print(原始数据形状:, raw.shape) # 期望 (N_sensors, N_time) Nx Ny 64 fs 48000 # 采样率, 单位 Hz f_target 1000.0 # 按行对应传声器来理解: 前 64 行是第一排, 后 64 行是第二排 p_time raw.reshape(Ny, Nx, -1) print(重组后形状:, p_time.shape) # (64, 64, 时间点数)逻辑说明reshape(Ny, Nx, -1)假设原始数据按“先空间后时间”存放且一行是一个通道。如果你的数据格式是每一列一个通道先对raw做一次转置再 reshape。频率索引用rfftfreq计算不要手写除法采样点数不同会导致索引偏一位n_time p_time.shape[2] freqs np.fft.rfftfreq(n_time, 1 / fs) idx int(np.argmin(np.abs(freqs - f_target))) # 加汉宁窗减少时域泄漏 window np.hanning(n_time)[np.newaxis, np.newaxis, :] p_windowed p_time * window # 沿时间轴做 FFT, 取出目标频率的复声压 P_cplx np.fft.rfft(p_windowed, axis2)[..., idx] print(复声压形状:, P_cplx.shape) # (64, 64)参数说明汉宁窗会压低旁瓣代价是主瓣变宽。如果目标是精确定位一个纯音噪声源加窗是划算的如果扫频信号本身谱线很宽可以跳过窗函数。rfft返回的频率点数是n_time // 2 1低频索引靠前高频索引靠后argmin这一步是为了防止你手动算索引算错。3.3 空间 FFT 与传播算子NAH 主循环拿到某个频率下的复声压分布后就进入了 NAH 的核心。下面这段代码实现了测量面到重建面的波数域传播dx dy 0.01 # 传感器间距, 单位 m z_dist 0.02 # 测量面到重建面的距离, 单位 m c0 343.0 # 空气中声速 # 二维空间 FFT, 注意前后各一次 shift Pk np.fft.fftshift(np.fft.fft2(P_cplx)) # 波数坐标 kx np.fft.fftshift(np.fft.fftfreq(Nx, dx)) ky np.fft.fftshift(np.fft.fftfreq(Ny, dy)) KX, KY np.meshgrid(kx, ky) # 圆波数与传播常数 k 2.0 * np.pi * f_target / c0 KZ np.sqrt(k**2 - KX**2 - KY**2 0j) # 回推到声源面: 传播算子 Pk_recon Pk * np.exp(1j * KZ * z_dist) # 硬截止滤波, 防止高波数噪声放大 k_c 1.2 * k mask (KX**2 KY**2) k_c**2 Pk_recon Pk_recon * mask # 回到空间域 p_recon np.fft.ifft2(np.fft.ifftshift(Pk_recon)) print(重建声压形状:, p_recon.shape)逻辑说明fft2之后的波数坐标顺序和 FFT 输出是错位的必须先fftshift把零波数挪到中心最后ifftshift还原这个顺序错了图像会四角分裂。传播算子里的0j是刻意的它强制sqrt返回复数负数的平方根才会变成纯虚数的倏逝波项否则会报 RuntimeWarning。参数说明dx是阵列相邻传感器的间距直接决定最高可重建频率z_dist是回推距离越大倏逝波放大越剧烈噪声越明显k_c是截止波数我用的是 1.2 倍圆波数也就是允许少量倏逝波通过。这个系数不是拍脑袋定的它跟测量面的信噪比强相关信噪比差就降到 0.8k信噪比很好可以放到 1.5k。3.4 可视化与 dB 标定别把声压当声压级重建出来的p_recon是复数数组值域可以到正负几十帕直接用imshow画虚部或实部会得到一张只有黑白两色的图因为数值动态范围太大。常见做法是把幅度转成 dB 再画import matplotlib.pyplot as plt p_db 20 * np.log10(np.abs(p_recon) / 2e-5 1e-12) vmax np.percentile(p_db, 99) vmin vmax - 40 # 显示 40 dB 动态范围 plt.imshow(p_db, extent[0, Nx*dx, 0, Ny*dy], originlower, cmapjet, vminvmin, vmaxvmax) plt.colorbar(labeldB SPL (re 20 μPa)) plt.xlabel(x / m) plt.ylabel(y / m) plt.title(fNAH reconstruction at {f_target:.0f} Hz) plt.savefig(results/nah_recon_1k.png, dpi150)参数说明参考声压 2e-5 是空气声学的标准参考值如果这份包里的数据本来就用 Pa 表示声压幅值画对比图时统一用 20 μPa 没问题如果数据已经做过归一化建议先确认量纲再写参考值。vmax取 99 分位数、vmin取vmax - 40是为了避免顶部一个特别亮的点把整张图压成暗色这个技巧在处理重建声压时比固定上下限安全得多。如果画相位图直接plt.imshow(np.angle(p_recon))不要对相位做 dB 转换相位是角度量不是幅度量。4. 参数调优与误差控制传感器间距、阵列孔径和截止波数怎么配合4.1 传感器间距决定空间奈奎斯特频率先检查硬上限阵列传感器间距是最先要确认的参数它直接划定了这道题能做多高的频率。空间采样同样要满足奈奎斯特条件最高频率可以这样估算f_max c0 / (2 * max(dx, dy)) print(f最大可重建频率约 {f_max:.0f} Hz)用 10 mm 间距的空气阵列最高可重建频率约为 17 kHz20 mm 间距只有 8.6 kHz。想做 1 kHz 到 10 kHz 的噪声定位10 mm 间距是合理的起点如果只需要低频段20 mm 间距可以减少通道数、降低成本。传感器间距 (mm)空气中最高可重建频率 (Hz)适合的频段534300全频段到超声边缘10171501~10 kHz 噪声定位208575低频机械噪声503430纯低频、大结构声源这个表按 c0343 m/s 计算实际温度不同会有几个百分点的差别。间距确定后如果目标频率超过上限降采样也救不回来只能重新设计测量阵列这是 NAH 里的硬约束。4.2 阵列孔径与重建范围低频分辨率不够时怎么办阵列的物理孔径 L N × dx对应波数域的最小分辨间隔 Δkx 2π / L。孔径越大波数域分辨率越高重建图像对声源的空间分辨越好。64×64 的阵列、10 mm 间距孔径是 640 mm在 1 kHz 时能分辨的声源细节大约在一个波长量级够用但不富裕。常见的做法是零填充。把复声压矩阵四周补零到 128×128 甚至 256×256可以让重建图更平滑看起来更细腻但要注意零填充不提高真实分辨率它只做插值。对 FFT 来说补零等于在波数域做 sinc 插值不会产生新的物理信息。另一个方向是减小间距但保持孔径不变这需要增加通道数成本随之上升。工程上如果低频分辨率不够优先考虑加大孔径而不是减小间距因为小间距只解决高频上限不解决低频空间分辨率。4.3 截止波数与正则化抑制高波数噪声的实际手法NAH 对噪声敏感的点在于倏逝波放大。回推距离越远这个放大越夸张。假设测点信噪比是 30 dB回推 50 mm 后高波数分量的噪声可能直接盖过真实信号这就是为什么必须在波数域加低通。硬截止就是我前面代码里的矩形窗(KX**2 KY**2) k_c**2实现简单但会在截止边界引入振铃。软截止更常用在截止区做一个余弦过渡k_rho np.sqrt(KX**2 KY**2) k_lo 0.9 * k # 开始衰减 k_hi 1.5 * k # 完全截止 window np.ones_like(k_rho) band (k_rho k_lo) (k_rho k_hi) window[band] 0.5 * (1 - np.cos(np.pi * (k_rho[band] - k_lo) / (k_hi - k_lo))) window[k_rho k_hi] 0 Pk_recon Pk * np.exp(1j * KZ * z_dist) * window逻辑说明这个窗从 0.9k 开始让幅度逐渐下降到 1.5k 处归零相当于把倏逝波部分“软着陆”地抑制掉比硬截止的振铃伪影小很多。k_lo和k_hi是调参的重心信噪比好可以调大到 1.2k 和 2.0k 保住更多细节信噪比差调到 0.6k 和 1.0k 优先保图像干净。Tikhonov 正则是另一种思路用1 / (|G|² λ)替代直接乘传播算子其中 λ 是正则化系数。它不会像截止窗那样一刀切而是按信噪比自适应压低低贡献波数分量。实际效果对 λ 的选择很敏感我一般先用截止窗调出基线图再用 Tikhonov 做对比观感差不多时选图像更平滑的那个。4.4 单源自检用已知位置点声源验证重建参数参数设得对不对空对空看理论没用。我拿到一个新包时第一件事是生成一个已知位置的点声源用代码生成“伪测量数据”跑一遍完整重建流程跟理论值对比误差。这比直接拿实验数据试错要快得多N 64 xx, yy np.meshgrid(np.arange(N)*dx, np.arange(N)*dy) # 假设声源在 (0.3m, 0.3m), 距离测量面 20mm src_x, src_y, src_z 0.3, 0.3, 0.02 dist np.sqrt((xx - src_x)**2 (yy - src_y)**2 src_z**2) # 点源在近场近似为球面波 p_meas np.exp(1j * k * dist) / dist # 加入一点噪声模拟真实测量 rng np.random.default_rng(42) p_noisy p_meas 0.01 * rng.standard_normal(p_meas.shape) * np.abs(p_meas).max() # 用 3.3 节的重建函数处理 p_recon reconstruct_plane(p_noisy, dx, dy, z_dist, f_target, k_c1.2*k) # 误差评估 p_true np.exp(1j * k * np.sqrt((xx - src_x)**2 (yy - src_y)**2 0.0**2)) / 1.0 err np.linalg.norm(p_recon - p_true) / np.linalg.norm(p_true) print(f相对重建误差: {err:.2%})参数说明伪测量数据里加了 1% 量的噪声这是模拟真实传声器阵列的最低噪声水平。err如果超过 30%先怀疑截止波数设得偏高降k_c重跑如果误差能压到 10% 以内这个参数组合就可以套用到实测数据上。这个自检脚本值得留在包里换一个测量面距离或换一批数据时先跑它做“参数指纹”能省一整天的盲调时间。5. 复现近场声全息的避坑清单五个我真实翻过车的点5.1 重建图出现棋盘格伪影现象重建声压图里出现规则的高频格子像像素级棋盘声源区域边缘尤其明显。原因空间混叠。传感器间距太大目标频率超出了空间奈奎斯特极限或者数据 reshape 时行列顺序搞反把 64×64 的网格硬拆成 32×128也会出类似格子。解决先用f_max c0 / (2 * dx)算硬上限确认目标频率没超。再检查 reshape 顺序打印P_cplx.shape和网格坐标对齐情况。我遇到过最气人的一次是 CSV 里每个通道的数据按列存放reshape 前没转置导致整个空间排列错了那种情况 c0 算出来再快也是白搭。5.2 图像中心出现同心环状振铃现象声源周围出现一圈一圈的同心圆条纹距离越远越密像水面波纹。原因测量面边缘的声压不是平滑趋近于零直接做空间 FFT 相当于对有限孔径内的信号做矩形截断频谱两侧产生 sinc 旁瓣。回推距离越大旁瓣越明显。解决在做二维 FFT 前对空间数据加 Tukey 窗。下面这行代码可以插在读取数据之后、FFT 之前from scipy.signal.windows import tukey w2d np.outer(tukey(Nx, 0.25), tukey(Ny, 0.25)) P_cplx P_cplx * w2d0.25表示边缘 25% 的区域做锥削太小压不住振铃太大又会把边缘真实信息削掉。这个数是我实测比较平衡的值你可以从 0.2 到 0.4 试一圈看图像边界干净为止。5.3 FFT shift 对不齐重建结果上下左右翻转现象重建出的声源位置跟实际测量位置镜像对称或者在图像四角各出现半个声源中间却什么都没有。原因只对输入做了fftshift没有在 IFFT 前做ifftshift或者两者顺序反了。FFT 输出的是从零频开始的半边谱fftshift把它移到中心做完处理后必须先把中心谱移回原始位置再做 IFFT否则时域空间坐标错位。解决严格按这个模板走不要凭感觉省略哪一步。Pk np.fft.fftshift(np.fft.fft2(P_cplx)) # ... 中间处理 ... p_recon np.fft.ifft2(np.fft.ifftshift(Pk_recon))注意这不是玄学fftshift和ifftshift对偶数长度序列行为一致对奇数长度序列有区别。阵列尺寸尽量用偶数64、128、256 都行奇数会让坐标对齐更麻烦。5.4 重建图几乎全黑或全白只有零星亮点现象保存的 PNG 图要么一片漆黑只有几个亮点要么白茫茫什么都看不见调节vmin/vmax也没用。原因复数数组里有 0 值或接近 0 的值log10(0)得到负无穷dB 数组变成 NaN。imshow遇到 NaN 就直接渲染成底色大范围 NaN 区域就会全黑或全白另外也可能是vmin/vmax设置离实际数据分布太远。解决log10前加一个极小偏移量 1e-12颜色范围用分位数而非绝对值先np.percentile(p_db, 1)和99扫一眼数据分布再定绘制范围。我习惯写一个简单的渲染函数内部固定用vmax p99, vmin p99 - 40这样每次输出的对比度是一致的不会因为某次数据里有异常尖峰导致整图失真。5.5 计算慢到像死机FFT 网格和复数数组类型问题现象64×64 的数据很快出结果换成 256×256 之后单次重建卡了几分钟还伴随内存占用暴涨。原因FFT 对尺寸敏感256 虽然本身是 2 的幂但如果数组之前被 pandas 或解析逻辑弄成了 object 类型NumPy 的 FFT 会退化成慢速路径另外如果数据里含复数却用 Python 原生 complex 列表来存内存占用会多出好几倍。解决进入 FFT 前强制转类型P_cplx np.asarray(P_cplx, dtypenp.complex128) N P_cplx.shape[0] # 如果 N 有大量质因子, 先补零到最近的 2 的幂 M 2 ** int(np.ceil(np.log2(N))) if M ! N: P_cplx np.pad(P_cplx, ((0, M-N), (0, M-N)), modeconstant)逻辑说明complex128是 NumPy 复数 FFT 的标准输入类型asarray强制转换可以避免 Python 对象逐元素访问的额外开销。补零到 2 的幂能显著提升 FFT 速度代价是重建图分辨率略微变化不影响相对声源分布。6. 进阶用法把重建流程封装成带参数的命令行入口直接批量出图6.1 命令行入口与常用默认值跑通单频点重建后下一步就是把它变成能复用的工具。最常见的方式是把主流程包进一个带argparse的函数让频率、距离、输出目录都变成显式参数避免每次都改脚本里的硬编码import argparse def main(): p argparse.ArgumentParser(descriptionnear-field acoustic holography reconstruction) p.add_argument(--freq, typefloat, default1000.0, helptarget frequency in Hz) p.add_argument(--z-dist, typefloat, default0.02, helpmeasurement-to-source distance in m) p.add_argument(--k-scale, typefloat, default1.2, helpcutoff wavenumber as k scale) p.add_argument(--data, defaultdata/holo_64x64_fs48000.csv) p.add_argument(--out, defaultresults) args p.parse_args() # 把 args 传给重建函数, save_figure 内部完成 dB 标定 reconstruct_and_save(args) if __name__ __main__: main()这样每次换数据或换参数只需要改命令行不用反复打开脚本编辑器。我一般会把默认参数设成上次标定好的最优值下次拿到新数据先跑一遍默认值再按结果微调。6.2 批量重建并输出声强向量实际项目往往不是看单一频率而是要看一段频带的噪声分布。批量循环几十个频点时输出声压图之外声强向量能更直观地指示声源方向。近场中粒子的振动速度可以通过声压的横向梯度来近似rho0 1.2 # 空气密度, kg/m^3 omega 2 * np.pi * f_target # 重建面上按中心差分求横向压力梯度 dpx np.gradient(p_recon, dx, axis1) dpy np.gradient(p_recon, dy, axis0) # 粒子速度近似: u grad(p) / (j * omega * rho0) ux dpx / (1j * omega * rho0) uy dpy / (1j * omega * rho0) # 声强向量: I 0.5 * Re(p * conj(u)) Ix 0.5 * np.real(p_recon * np.conj(ux)) Iy 0.5 * np.real(p_recon * np.conj(uy))逻辑说明这个近似假设声场局部满足欧拉方程重建网格足够细时误差可控。每个频率重建后保存一张声压 dB 图和一组 Ix/Iy 分量后续可以用np.load汇总成矢量图声源位置会直接体现在声强向量指向汇聚的地方。批量跑 20 个频点时加一个tqdm进度条会比干等舒服很多这也算是我这个工作流里最后一块拼图。从那以后我每次拿到一个新的 NAH 资源包不管代码写得多么眼花缭乱都会先强制走一遍“单点源自检 → 频点扫描 → 声强校验”这条固定管线。参数对不对、数据格式有没有理解偏在这套流程里都会现出原形。希望这份拆解能帮你在 organized7ah 上少花几个下午的盲调时间。本文还有配套的精品资源点击获取