ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

波数域是什么?从空间频率到声呐阵列处理的实战解读

波数域是什么?从空间频率到声呐阵列处理的实战解读 拿我们做声呐信号处理的人常被问到的一句话开头吧波数域到底是个啥听着挺唬人其实它一点都不神秘。你要是已经在时域频域里摸爬滚打过几年那波数域就是同一个套路换了个舞台——把随时间变化换成随空间变化把时间频率换成空间频率剩下的活儿还是傅里叶变换那套家底。这篇文章我就尝试用最直白的话把波数域、它的用处、还有实际处理里那些容易栽跟头的地方一次说清楚。适合刚接触阵列处理、合成孔径声呐SAS或近场声全息NAH的同行也适合那些做信号检测、波束形成做到一半发现怎么还有一层窗户纸的朋友。1. 从时间频率到空间频率波数到底在说什么1.1 频率概念的空间版波数与波长的关系先回到基础。我们处理时间信号时习惯把一段声压记录 s(t) 做傅里叶变换得到频谱 S(f)然后说这个信号里有哪些频率成分。这个操作背后的物理图景是任何信号都可以看成一系列不同频率正弦波的叠加频率 f 描述的是单位时间内相位转了多少圈。现在换一个视角。当声波打到一条水听器阵列上每个阵元在不同空间位置 x 上记录到的声压是不一样的。如果你把所有阵元的快拍数据拉出来就得到一条随空间位置变化的曲线 s(x)。既然它是一条随自变量变化的曲线我们当然也可以对它做傅里叶变换得到 S(k)。这个 k就叫波数wavenumber单位是 rad/m。波数的定义式很简单k 2π / λλ 是声波在水中的波长。波长越短波数越大意味着声压在空间里拧得越密。这个关系跟时间频率 f 1/T 完全对偶周期 T 越短频率越高信号在时间里抖得越快。用生活里的类比来说想象水面上的一排波纹。你站在岸边看某个固定点水面上下起伏的频率就是时间频率而你沿着波纹传播方向拍一张快照看到的波纹疏密程度就是波数。一个描述时间上的节奏一个描述空间上的纹理本质上都是周期的另一种表达。1.2 为什么非要再造一个波数域出来可能有人会问既然 f 和 k 通过声速 c 关系在一起ω c·k那我知道频率不就够了吗问得好这正是理解波数域价值的关键。频率和波数的关系只有在沿着声波传播方向看时才成立。但是阵列接收到的信号来自四面八方每个方向入射的声波在阵列这个一维投影上呈现出的空间周期是不一样的。设入射方向与阵列法线的夹角为 θ那么阵列方向看到的等效空间周期是 λ / sin θ对应的波数是k_x k·sinθ (2π/λ)·sinθ同一个频率的信号从不同方向来在空间上的表观波数完全不同。这就是为什么必须引入波数域频率域管不了方向波数域天然把方向信息编码在里面。波束形成、测向、空间滤波这些活儿本质上都是在波数域里划分哪些方向的分量我要哪些我不要。再说深一层。声波在均匀介质中传播时波数还分为传播波分量和倏逝波分量。当某个空间波数分量的模超过总波数 k也就是 k_x² k_y² k²对应的 k_z 变成虚数这个分量在空间上是指数衰减的能量传不远。这种物理现象只有站在波数域才能看清楚在纯时域频域里你根本意识不到它的存在。后面讲到近场声全息时这个性质就是核心技术。所以一句话总结波数域是空间维度的傅里叶对偶域它和频率域平起平坐一个是时间维度的频谱一个是空间维度的空间谱。声呐处理中凡是涉及方向、孔径、成像的问题绕开波数域就绕不开真正的物理本质。2. 波数域处理的第一站阵列波束形成的另一种视角2.1 时延求和与空间傅里叶变换的一体两面常规波束形成Conventional Beamforming是每个声呐工程师的入门课。原理讲起来很简单要对 θ 方向聚焦就给每个阵元加上对应的时延或相移补偿再把各路信号加起来同相叠加的信号增强不同相的自然抵消。这个操作从波数域角度看就更有味道了。设均匀线阵有 N 个阵元间距 d入射信号是来自 θ 方向的平面波。你给第 n 个阵元施加的相移正好等于把入射信号在空间采样序列上掰直求和的过程实际上就是在计算这个空间序列在某一特定波数 k_x k·sinθ 处的傅里叶系数。换句话说对空间维做一次离散傅里叶变换得到的是所有方向的波束输出波束形成就是在这条空间频谱上取值或选通常规波束形成 空间维的 DFT 加权求和所以你看波束形成并不是什么独立于傅里叶分析之外的新东西它就是空间傅里叶变换的工程实现。理解了这一点很多概念就不用死记硬背了。2.2 孔径、阵元间距与波数分辨率在波数域里几个关键参数变得特别直观。波数分辨率由阵列有效孔径 L N·d 决定Δk ≈ 2π / L这是空间傅里叶变换的固有分辨率极限跟时域里频率分辨率 1/观测时长完全对应。孔径越大波数分辨越细能区分的方向间隔越小。可观测波数范围由阵元间距 d 决定遵循空间奈奎斯特采样定理k_max π / d超过这个范围的空间频率分量会被折叠到可见范围内产生所谓的栅瓣grating lobe。换算成方向就是大家熟悉的栅瓣条件sinθ_栅瓣 sinθ_主瓣 ± λ/d我用一个具体例子帮大家找找手感。假设工作频率 100 kHz水中声速 1500 m/s波长 λ 0.015 m。取 8 元均匀线阵阵元间距半波长 d 0.0075 m。可观测波数范围k_max π / 0.0075 ≈ 419 rad/m对应端射方向 90°波数分辨率Δk 2π / (8 × 0.0075) ≈ 104.7 rad/m换算成角度分辨率Δθ ≈ λ / (N·d) 0.015 / 0.06 0.25 rad约 14.3°这个例子说明一件事想要更细的角度分辨就得加长孔径而不是单纯加阵元数想要更宽的观测角范围就得缩小阵元间距。这两个需求天生打架工程上永远在权衡。波数域的价值就是把这种权衡摆到明面上可用波数范围是采样间隔决定的波数分辨率是总孔径决定的一清二楚。3. 声呐里真正吃波数域红利的地方3.1 合成孔径声呐SAS波数域里的一次性聚焦做侧扫声呐或者合成孔径声呐的人对波数域应该是又爱又恨。SAS 的核心思想是用一个小物理孔径沿着航迹方向走出一条等效大孔径来从而提高方位向分辨率。这个等效虚拟阵列可能有几十米长这时候逐点做时域聚焦计算量会非常感人。工程上普遍采用 Omega-K 算法也叫 ω-k 算法本质就是在二维波数域里完成匹配滤波。处理流程大概是原始数据先做距离向 FFT 和方位向 FFT进入二维频域距离频率 方位波数乘以参考函数做匹配滤波补偿距离徙动和方位调制做 Stolt 插值将距离频率轴重新映射消除剩余的空间变异性逆 FFT 回到空间域得到聚焦图像为什么非得在波数域做因为聚焦操作在时域是逐点变化的卷积计算量大而在波数域里匹配滤波变成了逐点乘法Stolt 插值是一次性的频谱重排整体上能用 FFT 加速好几个数量级。实测数据我做过对比同样一条 1000 米长的航迹数据时域逐点算法跑完可能要几十分钟ω-k 算法几秒钟就出图了。当然天下没有免费的午餐。波数域算法对运动误差特别敏感航迹的微小偏移会在波数域相位上累积成线性项导致图像散焦。所以 SAS 处理必须搭配高精度惯导数据还要加自聚焦算法如相位梯度自聚焦 PGA来补偿残余相位误差。这个后面单独说。3.2 近场声全息NAH倏逝波帮我们看见亚波长细节近场声全息是波数域处理另一个典型应用也是我认为最能体现波数域物理价值的地方。NAH 的基本思路是这样的你在目标附近某个测量面上测到复声压 p(x, y)然后做二维空间傅里叶变换得到 P(k_x, k_y)。声波从源面传播到测量面在波数域里的传递函数非常简洁G(k_r, Δz) e^(j·k_z·Δz) k_z sqrt(k² - k_r²), k_r sqrt(k_x² k_y²)把测量面的波数谱除以或乘以这个传递函数再逆变换回去就能把声场回溯到源面得到源面上的声压分布。这一套操作在空间域其实是复杂的卷积积分瑞利积分但在波数域就是一个乘除法这就是红利所在。更精彩的是倏逝波部分。当 k_r k 时k_z sqrt(k² - k_r²) 变成虚数传递函数变成 e^(-|k_z|·Δz)是一个指数衰减因子。这意味着源面上那些高波数小尺度细节传播到测量面时已经被大大衰减如果测量面离源面非常近还能捕捉到这些倏逝波的残余从而恢复出亚波长级别的空间细节如果测量面太远倏逝波衰减到噪声底之下这部分信息就永久丢失了所以 NAH 才强调近场。远场测量为什么恢复不了高频空间细节不是算法不行是物理上倏逝波已经消失了。这个道理在波数域里一眼就看明白了。3.3 运动目标与频率-波数谱分析最后一个应用场景跟目标检测有关。对于静止介质中的声场频率和波数满足色散关系ω² c²(k_x² k_y² k_z²)如果目标是匀速运动的接收到的数据在时间和空间两个维度都会呈现特定的调制规律。对接收阵列数据做二维傅里叶变换时间维 FFT 空间维 FFT运动目标的能量会集中在频率-波数平面上的一条直线上直线的斜率携带目标的速度和运动方向信息。这就是频率-波数谱分析也叫 F-k 分析。这个方法在地震信号处理中用得很多声呐里处理水下运动目标、估计目标航速时也常用。它的好处是不需要先做波束形成再逐帧测多普勒而是在二维谱上一次性把方向-速度信息解耦出来直观且计算量可控。实际处理中要注意的是F-k 谱的解读需要对波数轴做归一化否则很容易把横轴单位搞错。另外宽带信号在 F-k 谱上会形成扇形展布处理时通常要划分频带或者做预白化不然目标轨迹会被淹没在色散展宽里。4. 波数域处理中的几个经典坑4.1 栅瓣空间欠采样的镜像效应栅瓣问题是波数域处理里最经典的坑也是空间采样定理的直接后果。前面算过阵元间距超过 λ/2波数域里就会出现周期性重复的谱峰这些谱峰对应空间上就是栅瓣。我在实际项目里见过一次比较典型的错误。某次阵列方案设计为了在有限通道数下做更大的孔径把阵元间距从半波长放宽到了 0.8λ。仿真时在 30° 方向放了一个目标结果发现 60° 左右出现一个几乎同样强度的假峰。当时负责算法的同事一度以为是算法 bug排查了很久最后翻到阵列参数才发现是栅瓣。这个坑的教训是波数域里出现的东西不一定是空间里真实存在的东西。栅瓣来自采样混叠是数学上不可避免的镜像。要么把间距压回 λ/2 以内要么接受它并在后续处理中通过宽带融合或非均匀布阵来抑制。后者的本质其实就是在波数域里打散混叠的周期性让它变成伪随机噪声而不是尖锐假峰。4.2 倏逝波留还是不留这是个问题前面讲到 NAH 里倏逝波携带高分辨率信息。听起来很美好但实际操作时倏逝波分量同时也在指数放大测量噪声。测量面上稍微一点点噪声回溯到源面时可能被放大几十上百倍图像直接花掉。所以工程上几乎没有谁会在回溯时把所有波数分量原样放大。常见的做法是加一个波数域低通滤波器在高波数段做平滑截止。这个截止波数的选择就是分辨率与噪声的平衡点截止波数设得高分辨率好但对测量精度要求极高噪声敏感截止波数设得低图像干净但亚波长细节被抹掉分辨率退化不同文献里有各种最优滤波器的推导但说实话到我手里的项目经验最可靠的还是根据实测数据的信噪比来试出来。先做一版不加滤波的回溯看噪声底在哪再决定切到哪里。这个艺术性的操作就是波数域处理里最需要经验的地方。4.3 窗函数与泄漏空间加窗就是波数域卷积波束形成和 NAH 处理里窗函数是另一个容易踩的坑。很多人只记得在时间维加窗忘了空间维同样需要加权。阵元幅度加权本质上就是空间窗函数它在波数域的效果与时间窗在频域的效果一模一样主瓣展宽、旁瓣降低。举个数字例子。均匀加权时8 元线阵的理论第一旁瓣电平大约是 -13 dB改用海明窗后旁瓣能压到 -40 dB 出头但主瓣宽度大约增加 40%。这个权衡在波数域看就是窗函数的频谱sinc 或加窗后的频谱与信号空间谱的卷积结果。理解了卷积视角你就能预判想要低旁瓣就得多占主瓣宽度世界上没有白吃的午餐。同理在波数域里直接截断频谱比如 NAH 滤波器的硬截止等价于空间域与一个宽 sinc 函数卷积会在图像边缘产生振铃。所以滤波器设计尽量用平滑过渡避免硬截止振铃会明显减轻。4.4 运动补偿波数域相位对误差零容忍这个问题在 SAS 处理中最突出。波数域算法的前提是完美的匀速直线运动但实际平台受涌浪、水流影响航迹总是有偏差。一个小的航迹误差 Δr在波数域相位上表现为 k_r·Δr 的相位误差而且这个误差随波数增大而放大。具体来说高频段大波数对运动误差极度敏感因为同样的位移在高波数下对应更大的相位翻转。这也是为什么 SAS 系统对姿态传感器精度要求那么高通常需要度数级别的航向精度和厘米级的定位精度。处理层面则要靠自聚焦算法来收拾残局先用惯导做粗补偿再用 PGA 或最小熵准则估计残余相位误差迭代去除。我自己的体会是做仿真时一切都很美好但一换上实测航迹数据问题全冒出来了。所以奉劝各位如果刚开始做 SAS 或合成孔径类处理千万别跳过运动误差的仿真环节——人为给航迹加一点随机扰动看看图像怎么劣化的这个过程比跑十遍理想仿真都有价值。5. 从理解到上手一份可执行的实操路线5.1 用一维均匀线阵把概念吃透如果你想短时间内把波数域从听过变成会用我建议按下面的路径做一遍仿真全程用 Python 的 numpy 就能完成半小时出结果生成一个 N16 元、间距半波长的均匀线阵仿真一个来自 30° 方向的单频平面波入射对快拍数据沿阵元维做 FFT注意是对空间维不要对时间维做画出空间谱观察谱峰对应的波数位置验证 k_x k·sinθ 是否吻合改变入射角到 60°再观察谱峰移动把阵元间距改成 0.8λ重新做看栅瓣怎么出现、出现在哪个角度把阵元数从 16 改成 8看主瓣怎么变宽加海明窗对比旁瓣和主瓣宽度的变化做完这几个实验波数域就不再是抽象名词了。你会发现空间 FFT 和时域 FFT 完全是一个妈生的所有直觉都可以平移。5.2 实操中容易翻车的三个细节坐标轴标错。空间 FFT 的波数轴是 k_x 2π·n/(N·d)其中 n 是 FFT 索引。很多人会顺手写成 n/N 或者漏掉 2π结果谱峰位置怎么都对不上。建议写代码时先把理论值算出来再对谱峰做验证。归一化问题。空间 FFT 做完后如果不除以阵元数 N能量会随阵元数增大而增大多组数据之间没法比较。无论用 numpy.fft.fft 还是 MATLAB 的 fft记得归一化。维度混淆。处理阵列数据时通常拿到的是阵元 × 时间快拍的二维矩阵。做波束形成要对空间维做 FFT做频域处理要对时间维做 FFT两个维度别搞混了。我见过有人把矩阵转置之后直接 FFT出来的东西怎么看怎么怪最后发现是轴搞错了。5.3 学习资料和工具推荐书籍方面阵列信号处理可以看《自适应滤波与阵列处理》这类经典教材合成孔径声呐方向推荐《Synthetic Aperture Sonar》和《合成孔径雷达算法与实现》在声呐上的迁移版本。近场声全息方向可以看《Fourier Acoustics》第 7、8 章那是我见过把波数域声学讲得最透彻的参考书之一。工具上仿真用 Python numpy/scipy 就够了FFT 直接调 numpy.fft。真要处理实测数据MATLAB 的 Signal Processing Toolbox 和 Phased Array Toolbox 也很顺手。关键是先把上面那 7 步仿真跑通再上水池或湖试数据最后才是海试数据。顺序别反了。我在实际工作中还有一个体会波数域不是孤立的知识点它跟时域频域是同一棵树上长出来的枝杈。很多声呐问题——波束形成、栅瓣分析、SAS 聚焦、近场重建、运动目标测速——看着各不相同但只要把它们翻译到波数域里立刻就能看出都是采样、变换、滤波这三个动作的排列组合。建议大家在处理任何阵列类数据之前先问自己一句这个问题在波数域里长什么样很多时候答案会比你想象中清爽得多。
返回列表