1. 先搞懂“波数”声呐里的空间频率1.1 我为什么想专门聊聊这个名词早几年调试一部多波束声呐的时候我最怕听到三个字波数域。那会儿日常工作已经习惯了画波束图在角度域里调阵列总觉得所谓“波数域处理”是另一套高深理论得翻开一堆数学书才敢碰。后来被一位老工程师拉着看了一次波数谱他指着屏上几个峰说“你看这个峰是目标旁边这个是栅瓣那边是干扰多清楚。”当场我就悟了——波数域不是另一套信号处理体系它其实就是我们天天用的波束形成只是换了个坐标系来观察问题。这篇文章面向的读者是那些已经在做声呐、水声信号处理或者正在入门阵列信号处理的朋友。哪怕你现在还分不清时域和频域只要耐心看完也能大概明白波数域处理是怎么一回事以及为什么它在实际声呐系统里这么常用。我不堆公式尽量用大白话和你能亲手跑起来的代码把这件事讲透。1.2 波数就是空间里的“频率”先做个类比。时间信号里我们常说“频率”指的是信号每秒变化了多少个周期单位是Hz。一个20kHz的正弦波就是每秒钟来回振荡两万次。这个“振荡得快慢”只跟时间有关跟位置无关。声波在水里传播时还有另外一个维度的振荡——空间维。你可以想象一列声波沿着某个方向往前推如果在某一瞬间拍张快照水里的声压沿着传播方向也是高低起伏的像一道一道的波纹。那么问题来了这道波纹在空间上是“挤得密”还是“拉得疏”这就要用到“波数”这个概念。波数的定义很干脆k 2π / λ其中 λ 是声波波长。波长越短k 越大说明声压在空间上变化越快。所以波数 k 本质上就是“空间里的频率”它描述的是声学量沿着空间方向变化的快慢。比如说同样是20kHz的声波在海水里大概以1500m/s左右的速度传播波长大约7.5cm算出来的波数就是 k ≈ 83.8 rad/m。这个数意味着你沿着声传播方向走1米相位要转83.8弧度也就是转了十几圈。如果写成平面波公式会更直观。一个沿x方向传播的声压可以写成p(t, x) A·cos(ωt - kx)这里的 ωt 管的是时间振荡kx 管的是空间振荡。你把时间 t 固定住只看 x 方向的变化那它就是一个个余弦波k 就是它每单位长度转过的相位。这个角度理解到位了波数域的门就推开了一半。1.3 阵列上的“花纹”不同方向的声波长得不一样声呐用到的阵列是一组在空间上按一定位置布放的换能器阵元。每个阵元其实就是个“空间采样点”跟时间采样一个道理。时间采样是等间隔地读取信号变化空间采样则是等间隔或不均匀地在不同位置感知声场。那阵列到底怎么“看见”目标的靠的是相位差。假设有一个平面波从某个角度 θ 传来约定0°是端射方向90°是正横方向。当它掠过一条均匀线列阵时每个阵元收到的波形其实是同一个波但到达时间不一样相位也就不一样。相邻两个阵元的相位差是多少这取决于波传播方向在阵列方向上的投影。如果把方向写成 u cosθ那相邻阵元的相位差就是Δφ 2π·(d/λ)·u其中 d 是阵元间距λ 是波长。这句话可以换个更漂亮的说法入射方向不同阵列上看到的“空间振荡快慢”就不同。u 越大目标越靠近端射方向波在阵列表面留下的“花纹”越密表现在相位差上就是相邻阵元走得越快u 越小目标越靠近正横方向花纹越疏相位差趋近于0。这个“空间振荡快慢”说白了就是某个方向上的空间频率也就是波数。只不过入射波本身的波数 k 2π/λ 是固定的阵列上看到的却是它在阵列方向上的投影即 kx k·u。所以阵列输出信号里天然就带着“空间频率”的信息——不同方向的信号对应不同的空间频率。我们要做的就是把混合在一起的空间频率分开。这件事就是波数域处理。2. 波数域处理跟波束形成是一回事对只是换了个坐标系2.1 从“时延求和”到“空间FFT”的三条路先想一个基础问题传统波束形成是怎么把某个方向的信号“捞”出来的最朴素的做法是时延求和。既然目标方向来的声波到每个阵元有到达时差那我就把每个阵元的输出往前补上这个时延再把所有阵元对齐后的信号加起来。对齐之后这个方向的信号同相叠加幅度增强其他方向的信号没有对齐叠加时互相抵消幅度变弱。这就是最经典的时延求和波束形成物理上最直观。到了窄带场景比如只关注某个频率点时延可以换算成相位旋转。比如想观察方向 u就乘上 e^{-j2π(d/λ)nu} 这样的加权系数再求和。每个方向对应一组加权系数于是可以用不同的 u 去扫描得到角度的功率谱。这是频域波束形成的视角也是大多数人先接触到的。但如果把目光从“逐个方向扫描”换成“把所有阵元数据一起做一次空间傅里叶变换”那就进入了波数域。注意空间FFT和窄带波束形成在数学上是完全等价的。对 M 个阵元的复包络 x_n 做X(g) Σ x_n · e^{-j2π n g}得到的 X(g) 就是波数谱峰值出现的位置 g 对应噪声源的空间频率。如果令 g (d/λ)u那么每个 g 直接对应一个入射方向。波数谱里哪个 g 有峰就说明哪个方向上来波能量强。三个视角一个内核。区别只在于实现路径一个是逐角度扫描的循环一个是一次性把整个空间频带拆开。2.2 波数谱的横坐标到底怎么读刚接触波数谱的人最容易懵的就在这横坐标到底是什么为什么有时候写 kx有时候写 u有时候又写 sinθ这里有一个约定问题。我在文里统一用 u cosθ0°为端射方向90°为正横方向。这时入射波在阵列方向上的波数分量为kx k·u (2π/λ)·cosθ而阵元间距 d 通常用波长归一化所以实际更常用的无量纲量是g (d/λ)·cosθ这个 g 就是空间FFT直接输出的横轴坐标。它为什么好用因为阵元接收信号的相位差写成 e^{j2π n g}n 是阵元序号g 相当于“每个阵元序号对应的周期数”正好是FFT输出的自然频率轴。举个例子。阵元间距 dλ/2入射方向如果是正横 θ90°那 u0g0所有阵元同相波数谱的峰会出现在整个阵列孔径的正中央对应波数0。如果目标偏到端射方向 θ0°u1g0.5波数谱的峰就跑到FFT横轴的最右端。如果目标在 θ180°u-1g-0.5峰会出现在最左端。换句话说波数谱的横轴范围对半波长间距的阵列来说正好是整个可观察的“空间频率带”。超出这个范围的 g 值对应的方向不存在或者说已经进入了“不可见区”。2.3 关键词归一化波数 g工程里的通用货币很多资料上又会写另一个量空间频率 fx sinθ/λ甚至会用“cycles/m”这种单位。听着很乱其实一回事只是归一化方式不同。只要记住一个关键换算物理波数kx 2π·fx (2π/λ)·cosθ归一化空间频率fx·λ cosθFFT自然轴g fx·d (d/λ)·cosθ在实际系统里我强烈建议你在代码和文档里统一用 g 这个无量纲量也就是“以阵元间距为单位的空间频率”。原因很简单不管频率变到多少只要 d/λ 确定g 的范围就是 [-d/λ, d/λ]如果你只看单边就是 [0, d/λ]。程序里画图、找峰、判断栅瓣全都用 g最后需要显示角度了再用 arccos 反算回来。这样能少踩很多坐标混乱的坑。3. 声呐系统为什么爱用波数域处理3.1 一次空间FFT等于同时扫出所有方向的波束第一个理由是效率。常规波束形成要形成 L 个波束每个波束对 M 个阵元做加权求和总共要做 M×L 次复乘加。如果 M256L512那就是13万次乘加还不算中间临时变量。而空间FFT呢256点的FFT约 M×log2(M) 2048 次复数蝶形运算一对复数乘法加加法大约算4次实数乘加也就在8000次左右。比常规波束形成少了一个数量级以上。更重要的是空间FFT天然把所有可能的波束方向都“算”了一遍。虽然 FFT 输出的 bin 对应的是均匀分布的 g 值不是均匀分布的角度但方向分辨率本来也接近均匀分布在小角度附近的。对于多波束测深声呐这种需要同时形成几十上百个波束的设备空间FFT几乎是必然选择。早年DSP性能紧张的时候很多系统就是靠一条 FFT 流水线把多波束撑起来的。3.2 波数谱是阵列的“体检报告”第二个理由是排查问题方便。在角度域看波束图你得逐个方向取最大值再拼成一张方位谱。目标一多、干扰一多图上就是一坨坨鼓包很难分清谁是旁瓣谁是栅瓣谁是真目标。但波数谱不一样。它直接把空间频率铺开所有成分都是竖线一样的分立峰。哪个峰是信号哪个峰是栅瓣位置在不在可见区内一眼就能判断。比如看到 g±0.5 的位置上有异常大峰而目标入射方向根本不可能是端射那就是阵列流形错误或者某几个阵元接反了。去年我调试一条64元线列阵时有一路水密连接器进水导致第17号阵元输出几乎为0常规波束图上看只是旁瓣稍高波数谱上却出现了一个非常规整的周期波纹顺着这个线索几分钟就定位到了问题阵元。3.3 宽带信号处理波数域是天然接口水下目标辐射噪声不是单频的而是宽带信号。工程上标准的做法是先把时域信号做FFT拆成多个窄带频点每个频点上的窄带阵列数据再做一遍空间处理。这种“频域空间域”两步走的架构本质上就是二维傅里叶变换的分离实现。波数域处理在这里的优势很明显不同频点上的波长不同直接放在一起比较没有意义但只要在频率轴上除以 λ或者用 g 做归一化不同频点就能对应同一套空间频率轴。于是你可以把多个频点的波数谱做非相干累加或者频域平滑得到更稳健的目标方位估计。这在窄带波束形成里是很难直接做到的因为角度谱的峰值随频率会漂移而波数谱经过归一化后不会。4. Python仿真从零跑通一次波数域处理4.1 仿真场景设计我先设定一个最标准的场景均匀线列阵16个阵元阵元间距取半波长这样全空间可见且不会有栅瓣。两个目标一个在50°方向一个在120°方向信噪比都不高20dB左右模拟真实环境里两个分得开但又不算太远的声源。采样上做窄带假设也就是只观察一个频点每个阵元输出一个复数快拍。为了看出“空间FFT直接出波数谱”的效果我故意加了一个小技巧FFT点数补到1024。数据只有16个点补零不会提高分辨率但能把谱线画得更圆滑峰值位置也更容易肉眼判断。这也是工程上常用的小手段想看谱形补零不亏。4.2 核心代码与逐段讲解import numpy as np M 16 # 阵元数量 d_over_lambda 0.5 # 阵元间距以波长为单位 thetas np.array([50.0, 120.0]) # 目标方向度0度为端射 snr 20.0 # 信噪比单位dB Nfft 1024 # 补零后的FFT点数 # 阵元位置以波长为单位 x np.arange(M) * d_over_lambda # 目标来波方向对应的 u cos(theta) u np.cos(np.deg2rad(thetas)) g_true d_over_lambda * u # 真正的归一化波数位置 # 随机复数幅度 rng np.random.default_rng(42) s rng.standard_normal(len(thetas)) 1j * rng.standard_normal(len(thetas)) # 阵列流形每个阵元对应每个目标的相位 A np.exp(1j * 2 * np.pi * np.outer(x, g_true)) # 理想接收数据窄带复数快拍 x_data A s # 加高斯白噪声 noise (rng.standard_normal(M) 1j * rng.standard_normal(M)) / np.sqrt(2) x_noisy x_data noise * (10 ** (-snr / 20)) # 空间FFT Xg np.fft.fft(x_noisy, Nfft) Xg np.fft.fftshift(Xg) # 横轴FFT自然频率轴单位是“每阵元序号多少周期” g_axis np.fft.fftshift(np.fft.fftfreq(Nfft, 1.0)) # 根据 d/λ 还原成 u再从 u 还原成角度 u_axis g_axis / d_over_lambda theta_axis np.rad2deg(np.arccos(u_axis)) # 画波数谱 import matplotlib.pyplot as plt power 20 * np.log10(np.abs(Xg) 1e-12) plt.figure(figsize(10, 4)) plt.plot(g_axis, power, lw1.2) plt.axvline(g_true[0], ls--, colororange, label50deg target) plt.axvline(g_true[1], ls--, colorgreen, label120deg target) plt.xlabel(g (d/lambda) * cos(theta)) plt.ylabel(Power (dB)) plt.grid(True, alpha0.3) plt.legend() plt.tight_layout() plt.show()这段代码里有几个地方特别值得注意。第一阵列流形 A 的构造用了 np.outer(x, g_true)把矩阵一次性算出来。这比一层层循环遍历方向干净得多运行速度也快。第二噪声功率控制是通过 10^{-SNR/20} 乘一个标准复高斯实现的记得复噪声要除以√2否则实部和虚部叠加后功率会偏大。第三fftshift 一定要用来把零波数放到中间输出顺序才是我们脑子里的正负频顺序不 shift 的话横坐标从低频到高频排列极容易看反。跑完这段代码你应该能看到波数谱在两个 g_true 的位置有明显凸起。50°对应 g0.5·cos50°≈0.3214120° 对应 g0.5·cos120°-0.25。一个在正半轴一个在负半轴完全分开。4.3 从波数谱回到角度目标怎么读出来谱里只是一个一个的峰要从波数峰还原成工程上习惯的“角度”步骤很简单找到峰值的横坐标 g_peak除回 d/λ 得到 u_peak再做 arccos 得到 θ。比如上面代码里如果程序输出 g_peak 0.3214那 u 0.3214/0.5 0.6428arccos(0.6428) 50°完美对应仿真设定。要注意arccos 这个函数在 u 接近±1时非常敏感在 u 接近0时则比较平缓。这带来的实际后果是波数谱在端射方向附近角度分辨能力会被拉伸得很差在正横方向附近角度分辨能力相对较好。这个性质不是 bug而是极坐标和笛卡尔坐标映射的固有几何特性。工程上做目标显示时我一般直接在波数域找峰不在角度域找峰找完再换算成角度这样最稳。5. 声呐工程师踩过的坑常见问题与排查5.1 栅瓣波数谱里的“假目标”栅瓣是阵列信号处理里最经典的坑波数域里尤其明显。当阵元间距 d 超过 λ/2 时归一化波数 g (d/λ)u 的取值会超出 FFT 的主周期范围。一个真实的 g会和 gnn±1,±2...出现在相邻周期里看起来就像多了几个假峰。比如 dλ 时g 的范围是 [-1,1]而 FFT 的主周期是 [-0.5,0.5]。于是 u 对应的真实峰在 g0.6 处同时还会在 g-0.4 处出现一个等高的假峰看起来就像是另一个方向的来波。这就是栅瓣。解决思路有两个一是在硬件设计上保证 d≤λ/2让全空间可见二是实在需要大孔径间距时把阵元排成非均匀阵破坏周期结构让栅瓣变成高旁瓣而不是等高峰。5.2 横坐标标错了方向直接差90°我刚用波数域的时候犯过一个无语的错误把 FFT 的横轴 g 当成 u 直接用结果算出来的角度全都偏得离谱。比如 g0.3d/λ0.5我用 arccos(0.3) 算出 72.5°但正确的 u0.3/0.50.6对应 53.1°。在正横附近误差看着不大靠近端射区直接就歪了十几度。这个问题的根源是 g 本身还含着阵元间距和波长的比。换算一定要做完两步先除 d/λ再做反余弦。不要跳步不要省那 0.5。5.3 波数谱“糊了”泄漏、窗函数与分辨率阵列孔径有限相当于在无限空间上截了一段。这个截断会带来频谱泄漏波数谱上原本细尖的峰变得“糊”旁瓣也抬高。跟时域加窗道理完全一致。想压低旁瓣就在空间FFT之前给阵元数据加窗比如汉宁窗、海明窗。加了窗旁瓣降了但主瓣变宽角度分辨率变差。这是不可调和的矛盾只能按场景取舍。另外要记得补零不会提高物理分辨率。16个阵元、孔径8λ波数主瓣宽度大约在 1/8 0.125按g为单位换算成角度在正横附近约7°。补到1024点只是把谱画细了两个距离小于物理分辨率的信号依然是一个鼓包。所以看到两个很窄的峰挨在一起别急着说“两个目标”先算一下分辨极限。5.4 近场不能用平面波套路波数域处理默认入射波是平面波也就是目标在无穷远处。可现实中目标往往就在近场比如港口水域的噪声源测距、浅海调查船拖曳阵离目标很近时波前是弯曲的。这时候同样的入射方向中间阵元和边缘阵元看到的相位关系就不再是简单的线性递增波数谱上的峰会展宽、偏位。解决办法是要么把数据先做距离聚焦再进入波数域要么用近场聚焦波束形成这时候波数域那套不加权的简单FFT就不灵了。具体工程上近场判据就是距离 R 是否远大于 L²/λL是阵列孔径。我习惯先算这个数不够大就放弃平面波假设。5.5 实操建议速查表现象可能原因排查/处理建议波数谱出现等间隔假峰阵元间距大于λ/2栅瓣检查d/λ改用非均匀阵目标角度整体偏移g没除d/λ就arccos先算ug/(d/λ)再反余弦谱峰很宽、分不清双目标阵孔径太小物理分辨率不足算分辨极限别靠补零改善谱峰旁边旁瓣抬头矩形窗的频谱泄漏空间FFT前加汉宁窗近场目标峰飘忽不定波前弯曲平面波假设失效确认远场判据改用聚焦算法波数谱像周期波纹个别阵元失效或接线错误检查该阵元输出观察空间FFT单点值最后分享一点个人体会。刚转用波数域那阵子我总觉得这是“理论派”才用的工具实际调声呐还是看波束图顺手。直到有一次排一个很隐蔽的阵元故障波束图和角度谱都看不出名堂换成波数谱后周期性波纹清清楚楚才彻底改观。现在我做阵列调试第一时间总是先把原始快拍拉成波数谱看一轮再决定往哪个方向追。波数域处理本质上就是“从空间维度做FFT”这个朴素操作它替代不了自适应、超分辨那些更高级的算法但它是最直观的第一道检视窗口也是理解水声阵列信号处理的一块关键拼图。后续如果你在学MVDR、MUSIC这类频域高分辨算法会发现它们很多前置推导都是站在波数域框架上的这也算是一条自然延伸的学习路径。