基于多线谱干扰抑制的水下对空运动声源线谱探测

刘凯悦，彭朝晖，张灵珊，王光旭，李晟昊

(1.中国科学院 声学研究所 声场声信息国家重点实验室，北京 100190；2.中国科学院大学，北京 100049)

0 引言

自20世纪50年代起，国内外学者开展了一系列有关空中运动声源透射入水问题的研究[1-2]，证明了水面附近的空气声源可等效为声源正下方十分接近水面的水下声源。因此，通常的水下目标探测方法理论上可用于水下对空探测。

但是，由于空气、海水介质特性阻抗存在巨大差异，空气中的声能量只有极少部分能够透射入水，加之海洋环境背景噪声的影响，传统的宽带能量检测器从水下探测空中声源相对困难。

空中声源的辐射噪声谱成分包括宽带连续谱和离散线谱，其中线谱强度高且稳定性好[3-5]，因此利用线谱探测是一种可能的途径。但是，由于空中声源的运动速度普遍较高，导致目标方位变化快、接收声场多普勒频偏大，因此水下目标探测技术相继产生的多种对线谱目标的探测方法[6-9]，如基于瞬时频率方差加权、基于线谱方位稳定性的目标检测方法等均难以应用于水下对空运动目标的线谱探测。

陈敬军等[10]提出一种基于人工智能的线谱自动检测算法，充分利用线谱形状特点及其在多个时刻的频率信息，根据最近时刻线谱频率外推得到下一时刻线谱搜索中心，有效解决了存在较大多普勒频移的线谱提取问题。文献[11]借助上述线谱检测算法，初步实现了水下对空中高速运动目标的远程探测。

实际浅海环境中，从水下探测空中目标会不可避免地受到来自多个水下源的干扰，其中大型水面船占主要成分。通常情况下，这些水面船的线谱基频很小[12]，在低频范围内具有较高的线谱能量，其线谱数量远多于空中目标激发的水下线谱，甚至存在同频干扰即目标与干扰的频率相同或十分相近的情况。因此，单纯依靠线谱检测方法难以有效地探测到空中弱线谱目标，还需要进一步进行干扰抑制。

国内外学者在对强干扰抑制的研究方面已总结出许多成熟的理论和技术，包括基于波束形成的干扰抑制算法[13-15]、基于信号分离技术的干扰抑制算法[16]以及基于子空间类的自适应干扰抑制算法[17-18]等。前两种算法分别在干扰与目标位于同一波束内、要求实时处理的情况下性能受损，相比之下，基于子空间类的干扰抑制算法具有更广泛的适用范围。任岁玲等[18]提出了基于特征分析的自适应干扰抑制算法，首先对接收信号的互谱密度矩阵进行特征分析，其次利用目标的方位范围信息剔除非目标信号子空间，重构互谱密度矩阵，实现自适应干扰抑制。该算法具有更稳健的干扰抑制能力，然而需要目标方位的先验信息，且在强干扰与目标轨迹交叉的情况下可能会对交叉后的探测结果产生误导。此外，文献[11]通过提取并剔除强干扰的线谱信息达到了干扰抑制的目的，但对于存在同频干扰的情况，该方法性能受限。

针对上述问题，本文提出一种基于多线谱干扰抑制的水下对空运动声源的线谱探测方法，以实现对同频干扰的抑制。通过2019年5月南中国海北部海域进行的水下探测空中直升机实验结果，验证了所提方法的有效性。

1 基于最小方差无失真响应的线谱探测方法

1.1 最小方差无失真响应

常规波束形成因其计算简单、稳健性高而广泛应用于阵列信号处理中。然而，由于其波束旁瓣较高，可能会导致弱目标信号被强干扰的“能量泄露”所淹没，抗干扰能力较差。基于此，采用最小方差无失真响应(MVDR)波束形成器[19]处理阵列数据。理论上，MVDR是多目标干扰情况下的最佳阵列处理器，具有较强的旁瓣抑制能力。

首先，给出频率分量ωk对应的MVDR波束形成器的加权向量(k=1,2,…,K，K为频点数)：

(1)

式中：ωk表示第k个频率值；a(ωk,θ)表示频率ωk、期望方向θ对应的阵列流形向量；上标H表示共轭转置；R(ωk)是由N个快拍估计得到的频率分量ωk的互谱密度矩阵；R-1(ωk)表示对R(ωk)的求逆运算。

对于宽带信号，经频率非相干处理后的宽带MVDR空间谱可表示如下：

(2)

当进行线谱检测时，需要利用MVDR波束形成器对每条线谱进行空间谱估计，并对所有线谱的估计结果进行综合处理。

假设分析时间长度为L，经线谱检测后得到I条线谱，设第i(i=1,2,…,I)条线谱在l(l=1,2,…,L)时刻共有K个频点，则第i条线谱在l时刻的空间功率谱估计可表示为

(3)

由于不同线谱的能量有所差异，为了提高对低信噪比目标的探测能力，下面对各条线谱能量进行归一化，即对各条线谱的空间功率谱估计结果进行归一化。由此，l时刻综合所有线谱贡献的归一化空间功率谱表示如下：

(4)

1.2 线谱提取

参照文献[10]中的算法进行线谱提取，通过模拟声纳兵在低频分析与记录(LOFAR)图中检测和提取线谱的过程，将线谱识别、检测和跟踪相结合，可以在低信噪比、线谱数量未知的情况下，以较高的检测概率实现线谱的自动提取。方法流程包括背景均衡[20]、线谱识别、线谱跟踪、噪声点剔除和丢失线谱产生。

下面以一次水下探测直升机实验为例，给出上述线谱检测算法的执行结果，实验对应的水文环境在3.1节中给出(但与3.1节中所给实验数据属不同批次)。图1(a)与图1(b)是利用上述线谱检测方法处理得到的线谱识别结果以及最终的线谱提取结果。图1(b)中不同目标的线谱分别用不同颜色的圆点表示，各干扰源线谱对应其他颜色圆点。当选取第6个时间点作为线谱搜索起始时刻时，该方法在前6个时间点内未能检测到直升机谱线。这是因为该方法在进行线谱跟踪时，只选取搜索中心时刻上存在的线谱点作为跟踪起始点，而未在搜索中心时刻处检测到线谱点的谱线可能会被忽略。

图1 文献[10]线谱检测各流程结果Fig.1 Detected results of line spectrum in each process in Ref.[10]

针对上述问题，为了降低线谱的漏检概率，修改线谱跟踪过程如下：以搜索中心时刻及其左右相邻时刻包含的所有线谱点(去除重复元素)作为线谱搜索的起始点，进行线谱跟踪，具体跟踪方法以及后续剔除噪声点、产生丢失线谱等步骤保持不变。

经上述步骤检测得到的线谱中可能会存在两条或两条以上重合度较高的谱线，它们实际对应同一条线谱，此外，在检测过程中对目标线谱频点的识别还有可能会存在偏差。基于此，在方法执行结束后增加线谱分离、线谱加密两步操作：

1) 线谱分离。对线谱跟踪结果进行归类，合并相同谱线。分为如下两步：①计算每条线谱在搜索时间区间内的平均频率值，将平均值之差在一定门限(可参考1 Hz)内的线谱划分至疑似合并线谱集合。②以每个集合内的第1条线谱作为参考线谱，选取其他线谱与参考线谱同时存在线谱点时对应的线谱信息，求取两组数据差值的标准差；若标准差小于一定门限，则与参考线谱视为同一条谱线，并将该线谱信息合并至参考线谱中，依次类推。

2) 线谱加密。对分离出的每条线谱在每个时刻分别进行加密，将线谱点左右一定频率范围内能量大于线谱检测门限的频点一并纳入该线谱。

利用改进后的线谱检测方法提取到的线谱如图2所示。从图2中可以看出，改进后的方法通过借助左右相邻时刻(第5、7个时间点)处的目标线谱点信息进行搜索，可跟踪得到完整的直升机谱线，并通过插值产生第6个时间点处的丢失线谱点。

图2 改进后线谱提取结果Fig.2 Extracted results of improved line spectrum

利用MVDR，分别对图1(b)、图2中的每条线谱做波束形成处理，并综合所有线谱贡献绘制方位历程图，如图3所示。图3(a)、图3(b)分别对应方法改进前后的检测结果。与线谱提取结果相对应，图3(a)中前6个时间点内均未检测到目标方位历程，而图3(b)中第2～11个时间点内直升机目标方位历程明显，与真实方位轨迹较为吻合。

图3 方位历程图Fig.3 Detected target bearings of two methods

此外，由于改进后方法对线谱跟踪结果进行了线谱点加密(见图2)，使得图3(b)中经加密处理后的目标方位更接近真实方位历程。因此，改进后线谱检测方法在一定情况下可相对提高线谱检测概率，增强线谱检测方法的探测性能。

2 基于多线谱干扰抑制的线谱探测

2.1 算法背景及思路

线谱检测对宽带背景噪声、宽带强干扰起到了很好的抑制作用，但是，当目标为空中声源时，近水面船因存在线谱成分而成为窄带强干扰。图4给出了直升机噪声入水实验中常见的水下接收声场功率谱密度曲线。从图4中可以看出，分析频段内干扰1、干扰2对应3条线谱，而目标只对应1条线谱；此外，3个干扰源均存在能量高于目标的线谱。即使对每条线谱的能量进行归一化，在综合所有线谱贡献后，仍存在干扰1、干扰2的空间谱能量大于空中目标，从而影响对目标的探测能力。因此，有必要在线谱探测的基础上继续进行干扰抑制。

图4 常见接收声场功率谱密度曲线Fig.4 Common power spectral density curve of receiving sound field

文献[11]与文献[18]所提出的干扰抑制算法分别在抑制同频干扰和同向干扰的情况下性能受损。基于此，本节提出一种针对多线谱干扰的抑制方法，以有效抑制同频干扰，且对于存在同向干扰的情况避免因干扰抑制而造成的目标方位信息损失。

首先，根据线谱提取结果，将各条线谱的顺序按照频率由高到低排列，利用MVDR波束形成器，依次对各条线谱进行空间功率谱估计。其次，对于每一时刻，当存在两条或两条以上线谱的空间谱最大值方位相同或相近时(认为这些线谱对同一信源起主要贡献)，只允许接收来自第1条线谱的全部贡献，抑制其他线谱在该信源方位附近的空间谱估计结果。最后，按照这一准则累积所有线谱贡献，得到方位历程图。

下面仍然以图4为例进行具体说明，图4中干扰2存在3条线谱，按频率轴从右到左依次估计每条线谱的空间谱。由于干扰2第2、3条线谱空间谱的最大值方位与第1条线谱空间谱最大值方位相同或相近，在累积线谱贡献时，只加入第1条线谱的完整空间谱以及第2、3条线谱空间谱经抑制后的结果。因此，最终的空间谱估计中相当于只有第1条线谱对干扰2起主要贡献。

2.2 算法具体流程

步骤1对第1个时刻，计算第1条线谱的空间谱并进行归一化，根据空间谱最大值对应的方位信息确定向量Θ1,1.

步骤2计算第2条线谱的归一化空间谱，若其空间谱最大值所对应的方位与Θ1,1中的元素相同，则找出空间谱中与该最大值相邻的左右两个极小值点(对应方位分别为θ1和θ2)，将介于θ1和θ2之间的空间谱替换为这两点之间的插值结果。

步骤3按照步骤2的方法，依次计算其他线谱的归一化空间谱，通过比较空间谱最大值方位与Θi-1,l中的元素，决定是否更新空间谱。

步骤4将所有线谱在第1个时刻处的归一化空间谱相加，即为第1个时刻干扰抑制处理后的空间谱P_IS1(θ).

步骤5按照上述步骤依次计算第2，3，4，…，L时刻的空间谱。

经过上述干扰抑制，l时刻第i条线谱的归一化空间谱可表示如下：

P_ISi,l(θ)=

(5)

综合所有线谱贡献，l时刻最终的归一化空间功率谱为

(6)

在每个信源所对应的空间谱中，该方法只允许一条线谱起主要贡献，从而抑制了多线谱干扰的能量。由于此处计算的是每条线谱的归一化空间谱，不同线谱之间的能量差异不予考虑。因此，选取哪条线谱作为干扰源的主要贡献者，对结果的影响不大。为了尽量提高线谱空间谱的方位分辨率，此处以同一信源的最高频线谱作为具有主要贡献的线谱。

为了提高方法的稳健性，可根据实际情况采用以下两种方法对参数进行微调。

1) 单条线谱空间谱方位分辨率较差，导致对空间谱最大值方位的估计存在一定误差，因此在判断第i条线谱对应的目标是否为多线谱目标时，可适当放宽条件：若第i条线谱对应空间谱的最大值方位与数组Θi-1,l中某一个元素值之差的绝对值小于一定的角度容限，则将二者视为同一信源，并抑制此线谱在最大值方位处的空间谱能量。

2) 当存在两个或两个以上强干扰源对应的线谱频率相同或相近时，线谱空间谱将同时出现多个较大极值。如果只抑制能量最大的干扰源，即只将最大值方位写入数组Θi,l而忽略其他极大值方位，则其他干扰源在其他频率处仍会贡献线谱，干扰抑制不完全。因此，记录入数组Θi,l中的每条线谱空间谱极大值的个数可根据实际情况在1～3之间进行选择。

综上所述，给出基于多线谱干扰抑制的水下对空运动声源线谱探测流程图，如图5所示。

图5 方法流程图Fig.5 Flowchart of the proposed method

3 海试数据验证

将本文多线谱干扰抑制方法应用于水下对空探测的海试数据中，并与文献[11]、文献[18]提出的两种干扰抑制方法相比较，验证该方法在实际工程应用中的有效性。分别选取两组来自不同实验的海试数据进行处理分析，两组数据中空中声源均以较高速度飞行，采用布放于海底的水平阵接收声信号，接收声场中均伴随有较多的水面船只干扰。

3.1 存在同频干扰，且目标与干扰的方位历程存在交叉

3.1.1 实验概况

2019年5月，在中国南海北部海域进行了一次海底水平阵探测空中直升机的实验。海底水平阵由56个阵元组成。实验期间直升机绕海底阵运动，运动轨迹如图6所示，运动范围始终保持在与阵相距2.5 km以内，声源的速度与高度随时间波动(此处直升机的位置及工况与1.2节中的实验数据有所差别)。实验海区海深约80 m，声速剖面如图7所示，水深5 m以内、5～50 m之间存在两个温跃层，呈现出不同的负梯度，水深50 m以下接近等声速分布。

图6 空中声源相对参考阵元运动轨迹Fig.6 Trajectory of airborne sound source relative to reference array element

图7 实验海区声速剖面Fig.7 Sound velocity profile of experimental sea area

3.1.2 线谱检测

选取低频范围作为分析频段，归一化后的频率范围为[6×10-3，12.4×10-3]。接收信号总时长为440 s，其中1次空间谱估计的积分时间为40 s，将其分成互不重叠的10个频域快拍，每个快拍对应的傅里叶变换点数为65 536.相邻时段空间谱估计的滑动步长为20 s，整段信号被分为21个时间段。

图8给出了参考阵元接收声场时频图，横坐标第1～21个时间点，分别对应1～21个分析时段的中间时刻。观察图8可以看出，分析频段内存在多条线谱，不同线谱之间的能量相差较大。按照1.2节所给方法分别进行背景均衡、线谱识别与跟踪、线谱分离，最终分离出的10条线谱如图9所示，相邻两条谱线用不同颜色的圆点进行区分。图9中黑色实线为根据声源全球定位系统(GPS)信息计算得到的接收声场瞬时频率曲线，由于声源的运动速度与运动方向都在不断变化，因此得到的瞬时频率曲线具有明显的波动。

图8 参考阵元时频图Fig.8 Spectrogram of reference array element

图9 线谱分离结果Fig.9 Result of line spectrum separation

对每条线谱进行加密，并利用MVDR波束形成器进行线谱检测，综合所有线谱的贡献，累积时间得到方位历程图，如图10所示。由图10可以看出，在190°和230°方向附近均存在较强的干扰，干扰的方位历程较为清晰、连续，从而掩盖了第7～11个时间点内的目标轨迹。此外，在第15～18个时间点处目标方位历程出现间断。总体而言，由于多个强干扰的存在，仅依靠线谱探测无法得到较清晰完整的目标方位历程。

图10 线谱探测结果Fig.10 Detection result of line spectrum

为了明确多个水面干扰源的存在对水下探测空中目标的影响，图11给出了各个时间点处空中目标与各干扰源之间的信干比随时间的变化曲线。选取图10中3个强度较大的干扰源，对应方位190°、227°和236°附近，分别计算目标与3个干扰源之间的信干比。由图11可以看出，目标与单个干扰源的信干比基本维持在0 dB以下，最低可达到-9 dB.因此，为了提高对空中目标的探测能力，需要寻求有效的干扰抑制方法。

图11 空中目标与水面干扰的信干比估计Fig.11 Estimation of signal-to-interference ratio between airborne target and surface interferences

3.1.3 干扰抑制

分别利用强干扰线谱剔除、基于特征分析的自适应干扰抑制以及多线谱干扰抑制3种方法进行数据处理，比较3种方法的探测效果。

1) 强干扰线谱剔除：由于对230°附近2个干扰的方位分辨力较差，此处选取190°干扰作为抑制对象，借助文献[11]中的方法对干扰进行线谱提取并予以剔除。

2) 基于特征分析的自适应干扰抑制：从图10中大致确定目标在初始时刻的方位为105°，后续时刻的目标方位由前一时刻经干扰抑制后的方位估计结果给出。空中运动目标方位变化较快，此处设置目标的方位搜索范围为前一时刻方位估计结果的左右15°区间内。

3) 多线谱干扰抑制：基于每个信源所在方位的空间谱主要接收一条线谱贡献的准则，依次更新每条线谱的归一化空间谱，并综合所有线谱贡献得到最终方位历程图。其中，角度容限取5°，记录入数组Θi,l中每条线谱空间谱极大值的个数(取为2)。

图12(a)、图12(b)、图12(c)分别给出了3种方法处理得到的方位历程图。从图12(a)中可以看出，190°方位附近的干扰已被抑制，目标方位历程在第1～8个时间点内清晰可见，但在第15～18个时间点处仍未显现。图13为190°干扰的低频率分析与记录(LOFAR)谱图，由图13可以发现，在第15～18个时间点内，190°干扰线谱与目标线谱的频率十分接近，因此在对190°干扰进行线谱剔除的同时损失了目标线谱信息，导致目标在该段时间内方位缺失。

图12 干扰抑制结果Fig.12 Results of interference suppression

图13 190°干扰LOFAR谱图Fig.13 LOFAR spectrogram of 190° interference

观察图12(b)发现，在第1～6个时间点内，由于算法对非目标信号子空间进行了抑制，使得190°和230°附近的干扰并未出现在方位历程图中。然而在第7个时间点处，目标与190°干扰发生轨迹交叉，由于干扰能量远大于目标，下一时刻得到的方位估计值为干扰所对应的方位，以至于自第8个时刻起，对目标的方位估计发生偏差，转变为对190°干扰的方位跟踪。此外，由于所给的目标方位搜索范围较大，对非干扰子空间的选取不够精确，导致方位历程图中的峰值旁瓣比较低，背景相对模糊。

相比图12(a)与图12(b)，图12(c)中的结果既显示出较完整的目标、干扰运动轨迹，又保持了较高的峰值旁瓣比，便于同时对目标与干扰进行分析。图12对应的海试数据处理结果表明，本文所提方法对于解决存在同频干扰或目标与强干扰轨迹交叉的线谱探测问题具有明显的优势。

3.2 不存在同频干扰

为了证明本文所提方法同样适用于抑制非同频干扰，下面用该方法处理文献[11]中直升机探测实验的海试数据。本次实验与3.1节中所述实验毫无关联，具体实验细节可参考文献[11]。图14(a)是利用所提方法处理得到的方位历程结果，为了进行对比，在图14(b)中给出了文献[11]中的处理结果。对比发现，图14(a)中目标的方位历程在前11个时间点内更为连续，能量更加集中；在第12～13个时间点内，由于150°干扰与目标的运动轨迹存在交叉，且干扰的能量相对较强，使得目标的方位历程在后续时间段内衔接性较差。总体而言，利用所提方法可有效处理文献[11]中的实验数据，实现非同频干扰下水下对空运动声源的线谱探测。

图14 目标探测结果对比Fig.14 Results of target detection

4 结论

为了解决水下对空中高速运动声源的探测问题，本文在线谱探测的基础上提出了一种对多线谱干扰的抑制方法。针对分析频段内干扰源线谱数量远大于目标线谱且干扰线谱能量较强的特征，限定每个信源的空间谱主要接收来自一条能量归一化的线谱的贡献，最后综合所有贡献线谱的空间谱估计得到最终的方位估计结果。由于分析频段内目标与干扰的空间谱各自只对应一条线谱的贡献，且对线谱的能量进行了归一化，因此目标的方位历程得到增强。本文方法可同时有效抗多个干扰源，且在干扰与目标同频的情况下依旧保持较好的性能。海试实验结果表明，在多个近水面船干扰情况下，基于多线谱干扰抑制的线谱探测方法可成功实现水下对空中运动声源的探测，相比其他方法具有更广泛的适用范围。