声呐阵列信号处理：波数域、空间FFT与波束形成的本质

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°，那 u=0，g=0，所有阵元同相，波数谱的峰会出现在整个阵列孔径的正中央，对应波数0。如果目标偏到端射方向 θ=0°，u=1，g=0.5，波数谱的峰就跑到FFT横轴的最右端。如果目标在 θ=180°，u=-1，g=-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 次复乘加。如果 M=256，L=512，那就是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, lw=1.2) plt.axvline(g_true[0], ls='--', color='orange', label='50deg target') plt.axvline(g_true[1], ls='--', color='green', label='120deg target') plt.xlabel('g = (d/lambda) * cos(theta)') plt.ylabel('Power (dB)') plt.grid(True, alpha=0.3) plt.legend() plt.tight_layout() plt.show()

这段代码里有几个地方特别值得注意。

第一，阵列流形 A 的构造用了 np.outer(x, g_true)，把矩阵一次性算出来。这比一层层循环遍历方向干净得多，运行速度也快。第二，噪声功率控制是通过 10^{-SNR/20} 乘一个标准复高斯实现的，记得复噪声要除以√2，否则实部和虚部叠加后功率会偏大。第三，fftshift 一定要用来把零波数放到中间，输出顺序才是我们脑子里的正负频顺序，不 shift 的话，横坐标从低频到高频排列，极容易看反。

跑完这段代码，你应该能看到波数谱在两个 g_true 的位置有明显凸起。50°对应 g=0.5·cos50°≈0.3214，120° 对应 g=0.5·cos120°=-0.25。一个在正半轴，一个在负半轴，完全分开。

4.3 从波数谱回到角度：目标怎么读出来

谱里只是一个一个的峰，要从波数峰还原成工程上习惯的“角度”，步骤很简单：找到峰值的横坐标 g_peak，除回 d/λ 得到 u_peak，再做 arccos 得到 θ。比如上面代码里，如果程序输出 g_peak = 0.3214，那 u = 0.3214/0.5 = 0.6428，arccos(0.6428) = 50°，完美对应仿真设定。

要注意，arccos 这个函数在 u 接近±1时非常敏感，在 u 接近0时则比较平缓。这带来的实际后果是：波数谱在端射方向附近，角度分辨能力会被拉伸得很差；在正横方向附近，角度分辨能力相对较好。这个性质不是 bug，而是极坐标和笛卡尔坐标映射的固有几何特性。工程上做目标显示时，我一般直接在波数域找峰，不在角度域找峰，找完再换算成角度，这样最稳。

5. 声呐工程师踩过的坑：常见问题与排查

5.1 栅瓣：波数谱里的“假目标”

栅瓣是阵列信号处理里最经典的坑，波数域里尤其明显。当阵元间距 d 超过 λ/2 时，归一化波数 g = (d/λ)u 的取值会超出 FFT 的主周期范围。一个真实的 g，会和 g+n（n=±1,±2...）出现在相邻周期里，看起来就像多了几个假峰。

比如 d=λ 时，g 的范围是 [-1,1]，而 FFT 的主周期是 [-0.5,0.5]。于是 u 对应的真实峰在 g=0.6 处，同时还会在 g=-0.4 处出现一个等高的假峰，看起来就像是另一个方向的来波。这就是栅瓣。解决思路有两个：一是在硬件设计上保证 d≤λ/2，让全空间可见；二是实在需要大孔径间距时，把阵元排成非均匀阵，破坏周期结构，让栅瓣变成高旁瓣而不是等高峰。

5.2 横坐标标错了，方向直接差90°

我刚用波数域的时候犯过一个无语的错误：把 FFT 的横轴 g 当成 u 直接用，结果算出来的角度全都偏得离谱。比如 g=0.3，d/λ=0.5，我用 arccos(0.3) 算出 72.5°，但正确的 u=0.3/0.5=0.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 实操建议速查表

最后分享一点个人体会。刚转用波数域那阵子，我总觉得这是“理论派”才用的工具，实际调声呐还是看波束图顺手。直到有一次排一个很隐蔽的阵元故障，波束图和角度谱都看不出名堂，换成波数谱后周期性波纹清清楚楚，才彻底改观。现在我做阵列调试，第一时间总是先把原始快拍拉成波数谱看一轮，再决定往哪个方向追。波数域处理本质上就是“从空间维度做FFT”这个朴素操作，它替代不了自适应、超分辨那些更高级的算法，但它是最直观的第一道检视窗口，也是理解水声阵列信号处理的一块关键拼图。后续如果你在学MVDR、MUSIC这类频域高分辨算法，会发现它们很多前置推导都是站在波数域框架上的，这也算是一条自然延伸的学习路径。