一种波束谱特征加权的水下弱目标检测方法

范文涛 唐 波 李朋伟 王英志

（中国人民解放军92578部队，北京 100161）

1 引言

水下被动声检测核心问题是如何从复杂时变、充满干扰的环境中检测出感兴趣的弱目标信号［1］。常用处理手段是通过波束形成将无源传感器阵列接收到的被动声信号进行时间域和空间域处理，从而挖掘出目标的宽带和窄带信号特征，实现后续目标跟踪和识别。类似于雷达和语音信号处理领域［2］，应用于水下被动声检测的波束形成方法可分为两类：常规和自适应。其中最常用的检测算法是常规波束形成（Conventional Beamforming，CBF）［3］，该方法可以适配各种水下感知系统和海洋环境，因此也非常受到水下信号处理工程师的青睐。但是，随着安静型潜艇的出现，所辐射的噪声级越来越低，导致传统CBF很难检测到目标。一些研究者考虑使用最优阵处理来解决这个问题，尽量压低波束图的旁瓣，使得系统获得更大的阵增益，但是这种方法存在计算量大的缺点。其中一个著名的检测算法叫最小方差无失真响应（Minimum Variance Distortionless Response，MVDR）波束形成器［4］，这种方法主要是基于信号和噪声协方差矩阵来求解最优加权函数，使得波束形成器能够自适应的将波束图凹槽对准干扰方向，将主瓣对准目标。MVDR 波束形成器比起CBF 方法大大增强了在多目标分辨率和干扰抑制方面的性能，但是MVDR 波束形成器在实际系统中应用时不够稳健，尤其是快拍数较少时，需要对接收信号协方差矩阵进行对角加载，不同的加载量会造成一定的检测增益损失，到目前为止，还未出现一种既能保住MVDR 检测性能又能增强其稳健度的对角加载方法。

最近，一些研究者尝试将CBF 或者MVDR 波束形成器与其他的非线性处理方法融合用以提高传统方法在水下被动声检测中的检测性能。盲源分离（Blind Source Separation，BSS）［5］就是一种在语音信号处理领域中较著名的处理方法，该方法主要是通过从接收的阵列混合信号矩阵中提取出感兴趣的信号。小波变换（wavelet transform，WT）［5］和经验模态分解（Empirical Mode Decomposition，EMD）方法［6］也被考虑用以改善听音波束、LOFAR 谱［2］、DEMON 谱［2］的信噪比。还有二阶锥处理［4］和遗传算法［3］也被用来帮助波束形成器找到适合自己的最优权向量。但是遗憾的是，所有这些非线性处理技术都存在两个问题，一是需要对环境和基阵进行假设，二是工程化实现起来非常复杂。

在反潜战领域，大孔径线列阵通常被用来对安静型潜艇目标搜索和跟踪的［2］，例如，舷侧阵、拖曳阵、海底固定阵［4］、高级分布式阵［2］等等。但是这些阵都面临一个相同的问题，就是当目标运动到阵列端首附近时，波束形成器的检测能力就极大地下降。一般称基阵端首方向为盲区（Blind Zone，BZ）［4］。但是，这个问题又不能回避，因为BZ区域太重要了，从战术上讲，鱼雷就喜欢从这个方向对舰船发起攻击。一些学者提出用逆波束形成来抑制盲区的本舰自噪声强干扰［4］，但是该方法是利用语音中谱减的思路来将所有波束中的干扰剪掉。要知道谱减的过程中，信号也会被减掉，这样会引起检测背景的提升。

本文提出一种新的波束输出谱加权方法，步骤为第一步，从波束输出中估计背景噪声级；第二步，根据背景噪声级从波束谱中提取感兴趣的频谱；第三步，基于ECKART 最优检测滤波器估计波束最优权值；最后，权值被加到波束形成器输出波束上完成弱目标检测。通过使用实际拖曳阵和舷侧阵的数据进行验证，从实验结果来看，所提方法在浅海多干扰环境下，对弱目标能够获得4 dB 左右的检测增益。

2 宽带模型

考虑M个水听器组成的均匀线列阵，假设远场有D个宽带声源（D＜M）以平面波的形式到达基阵，波达方向为Θ=[θ1，...，θd，...，θD]。那么第m个水听器所接收到的时域信号为

其中sd(t)表示第d个声源信号，ξn(θd)表示第d个声源相对于第n个水听器的传播延迟，vn(t)表示叠加的高斯白噪声。由于被动声纳基阵接收的信号为宽频带信号，而大部分方位估计算法都是基于窄带假设，所以处理时一般通过短时傅里叶变换（Shorttime Fourier Transform，STFT），将宽带信号分段，并划分为互不重叠的多个窄带，然后再在每个子带应用相应的信号处理算法［7］。那么，基阵输出的M×1维快拍向量数据的频域响应可表示为

其中fj表示频点（frequency bin），频点数就是STFT对阵元域数据分段后，每个段做FFT 的点数，频点数就是子带数。s(fj)=[s1(fj)，...，sD(fj)]T为D×1维的源信号向量（(·)T表示转置）。A(Θ，fj)=[a(θ1，fj)，...，a(θD，fj)]是M×D维的阵列流形矩阵。导引矢量a(θi，fj)［8］是：

3 所提出的方法

图1 是本文所提方法的框架结构。在图中，阵列接收的时域信号被STFT 变换成频域多子带信号，以保证窄带信号处理算法能够工作。公式（2）中的X(fj，m，k)表示第m时间段下第j个频点下的第k根谱线：

第二步是自适应波束形成，本文中使用MVDR算法［4］。MVDR 方法的权值通过对如下公式的最小化来实现：

其中R(ω)是协方差矩阵［7］。在实际应用中，一般对协方差矩阵采取时间平滑处理措施，则协方差矩阵为［7］：

公式（5）是一个趋向于单位1的约束优化过程；

公式（5）和（6）可以通过拉格朗日乘子法来求解，所得到的结果为［7］：

频域MVDR波束输出为：

因此，频率-方位矩阵（Frequency Azimuth，FRAZ）可表示为：

ZFRAZ(fj)是J×B的矩阵，J是代表频点数和B是波束数。

利用分裂双通窗（Two-Pass Split Window，TPSW）［8］算法来估计每个波束的噪声级（Noise Level，NL）并且以此提取出背景噪声谱。TPSW 是一种常用的估计噪声谱的处理技术，计算量少，实现简单，其窗函数的公式为［8］：

TPSW 将沿着信号滑动来获得噪声平均值。TPSW最后输出背景噪声谱矩阵NFRAZ(fj)：

在宽带检测理论中，ECKART 滤波器［8，10］被认为是一种能够最大化输出信噪比的最优检测器。但是，ECKART 需要对信号和噪声具有一定的先验信息，实际上信号和噪声谱是未知的，需要估计信号和噪声功率谱。实际被动阵列信号的FRAZ 矩阵中，信号谱分量往往表现出一定的空间相关性，而噪声谱分量往往空间不相关，本文使用矩形窗在FRAZ 矩阵中沿着每三个波束频谱分量中表现为“峰值”的值：

公式（13）寻找到的“峰值”zpeak(fj，θb)被认为候选信 号谱分量。反之，被定义为“谷值”：

由此，本文设定规则：当公式（13）估计得到的“峰值”大于对应的公式（12）背景噪声谱分量值，ECKART 滤波器的信号分量即为该“峰值”，噪声分量为常数1；当“峰值”小于对应的背景噪声谱分量值时，ECKART 滤波器的信号分量即为该“峰值”，噪声分量为背景噪声谱分量的平方；当公式（14）估计得到的“谷值”大于对应的公式（12）背景噪声谱分量值，即认为该“谷值”为噪声或者干扰，对该值叠加背景噪声谱值后赋为负值，作为抑制；当“谷值”小于对应的噪声谱分量值时，即直接对该值赋负值。因此，波束谱特征加权值为：

因此，波束形成器的输出波束可以被加权为［11］：

最后，宽带方位历程谱能通过非相干累加Tnew(fj)的方式得到。所提方法处理框架如图1所示。

4 数值试验

仿真采用阵列信号模拟器产生数据，场景设置为：拖曳线列阵为41元水听器均匀线列阵，其阵元间距满足半波长的假设；拖船辐射噪声源级设为SL拖船=120dB，拖速v拖速=14节，拖缆长度R拖缆=400米，阵深为30米；目标A辐射噪声源级设为SL目标A=120dB，距离拖曳阵R目标A=5.4海里，航速v目标A=18节；目标B 辐射噪声源级设为SL目标B=105 dB，距离拖曳阵R目标B=2.7 海里，航速v目标B=8节；目标C辐射噪声源级设为SL目标C=115 dB，距离拖曳阵R目标C=3.7 海里，航速v目标C=14节；目标A、B、C初始状态分别位拖曳阵舷角是95°、115°和135°，向拖曳阵端射方向运动。其中阵列信号数据采样率fs=48 kHz，FFT 长度N=1024，频率分辨率Δf=fs/N=46.875 Hz，所用的数据段总点数为48000 采样点，将该数据段进行分段，每段长度为1024，相邻两段数据重叠率为50%，处理频段为100～600 Hz，相应的频点数为93。MVDR 方法和本文所提方法均采用“白噪声增益”对角加载，加载量取0.092。

图2 是MVDR 与本文所提方法的宽带警戒方位-时间-幅度瀑布显示图。对比（a）、（b）子图可以看出，本文所提方法对三个目标的检测增益要高于MVDR 方法，尤其是过了40度接近拖曳线列阵声纳检测盲区时，本文所提方法对目标C依然保持接触；本文所提方法的检测背景起伏要大于MVDR 方法，原因在于本文所提方法在提取谱特征过程中，会存在将部分噪声谱分量作为信号提取的问题，导致检测背景出现起伏。

图3 为对弱目标利用蒙特卡洛仿真方法，比较不同输入信干噪比、距离和信噪比条件下，本文所提方法与MVDR 方法的阵增益性能，其中干噪比为10 dB，设定-3 dB为检测阈值。从子图（a）中可以看出，本文所提方法在输入信干噪比为-34 dB 时就能检测目标，比MVDR 方法小了9 dB；子图（b）可以看出，本文所提方法的检测距离比MVDR 方法要远5 公里左右；子图（c）可以看出，MVDR 方法阵增益信噪比门限比本文所提方法要低3 dB，意味着本文所提方法的抗噪性能不如MVDR方法。

可以得到以下结论：

（1）从方位历程图的检测结果来看，本文所提方法增大了背景的起伏，但是对于弱目标尤其是盲区附近目标的检测性能要好于MVDR方法。

（2）表1 为总结图3 得出的有关两种方法的弱目标最小可检测性能，包括最小输入信干噪比、最远检测距离和最小输入信噪比三个度量指标。从表中可以看出，本文所提方法比MVDR 方法的最小输入信干噪比低了9 dB；本文所提方法比MVDR 方法的最远检测距离多了4 公里；本文所提方法比MVDR方法的信噪比门限高了3 dB。

表1 弱目标最小可检测性能表Tab.1 Weak target minimum detectable performance table

5 实际海试结果

本文利用舷侧阵数据验证算法性能。水听器阵为32 阵元等间隔直线阵，间距1 米。试验中目标以6 节航速定深航行，有合作目标。根据多次测量中比较典型的不同阵元数和不同频率上的噪声能量分布，设定检测处理频段为310～630 Hz，采样频率为6000 Hz。处理数据时间长度为3000 s，该段时间内共有6 批目标。将该数据段进行分段，每段长度为2048，相邻两段数据重叠率为75%。MVDR 方法和本文所提方法均采用“白噪声增益”对角加载，加载量取0.092，积分时间为2 s。本次试验中没有合作目标，通过目力观测和导航雷达记录，观察到了6批次水面渔船和商船目标，其中起始于170度舷角的是大型商船，其余为渔船。

图4 子图（a）是MVDR 方法的宽带警戒方位-时间-幅度瀑布显示二维、三维对比图。从子图（a）检测结果来看，目标A、目标B、目标C、目标D、目标E和目标F度附近有存在目标，其中目标C 能量较强，目标A、目标B、目标D 和目标E 的前1500 s 历程非常弱，目标F 从舷侧阵端首往另一侧端首（0 度）运动，靠近端首盲区的历程部分MVDR 检测效果不佳，出现了方位模糊现象。

图4 子图（b）为本文所提方法得到的三维和二维方位历程显示，可以看出本文所提方法比MVDR方法对六批目标的历程跟踪较为清晰，其中对于起始于目标A和目标B的两个渔船目标检测增益明显提升，全程历程显示清晰；对于起始于目标D 和目标E 的两个渔船目标，也较为明显的将其前1500 s的历程呈现出来；对于170 度目标靠近两端端首部分的历程显示最清晰完整。本文所提方法在100度至120度方位扇区内的检测背景起伏要高于MVDR方法，原因在于受到正横目标C 附近强干扰能量泄露的影响，本文所提方法在提取信号谱特征分量过程中，将泄露至该方位扇区内的干扰噪声分量作为信号分量提取，导致最后在该方位扇区内的背景起伏大于MVDR方法。

表2为舷侧阵数据阵增益的统计结果。由于该实录数据无法预知准确的目标态势、信干噪比、干噪比和信信比的信息，因此本文考虑选取整个数据分析历程中的四秒（10 s、500 s、1000 s 和2000 s）时间，针对目标A 附近的弱目标方位谱峰值与平均背景的对数比值作为检测阵增益，经过比较本文所提方法比MVDR方法平均阵增益提高了2.785 dB。

表2 舷侧阵数据阵增益统计表Tab.2 Flank array data array gain statistics table

由于本文所提方法，为了更加直观地说明本文所提方法的复杂度，将本文所提方法与MVDR 方法处理同一批试验数据，统计观察时间2 s内方位谱输出所需的计算时间。计算机硬件采用Intel i7-8500CPU，主频3.1 Hz，内存8.0 GB。计算软件采用MATLAB2019a版本。如表3所示，MVDR，本文所提方法所需计算时间分别为0.83 s和0.97 s。对比可以看出本文所提方法基于MVDR方法波束输出谱增加了特征提取和ECKART 滤波两个步骤，相比MVDR方法慢了0.14 s，计算复杂度没有明显提升。

表3 方法性能比较表Tab.3 Method performance comparison table

6 结论

本文提出一种波束谱特征加权的水下弱目标检测方法，通过对波束输出信号进行背景级估计，从而利用ECKART最优滤波器从波束输出谱中提取感兴趣的目标谱信号，实现对弱目标的检测。计算机仿真和实际海试数据证明该方法对于水下被动弱目标检测具有一定的效果，同时所提方法基于MVDR 方法波束输出谱进行特征处理，计算复杂度未明显提升，利于在实际水下传感器处理系统上工程实现。