李海鹏 孙大军 郑翠娥

(哈尔滨工程大学水声技术重点实验室 哈尔滨 150001)

(海洋信息获取与安全工信部重点实验室(哈尔滨工程大学) 哈尔滨 150001)

(哈尔滨工程大学水声工程学院 哈尔滨 150001)

1 引言

海洋作为地球最大的生态系统影响着全球能量流动、气候变化与生态安全，将地球连结为一个命运共同体。我国是一个海洋大国，海洋面积辽阔，认识海洋、经略海洋、建设海洋强国具有重要的战略地位。深海面积超过海洋总面积的90%，走向深海是海洋强国的必由之路。对深海环境特性的精确认知和对深海资源的科学开发利用是建设海洋强国的基础[1–3]。随着对海洋特别是深海探索和开发的深入，对各类水下潜器、平台的高精度定位导航需求越来越强烈。水声定位系统是现代深海作业必备的高精度水下定位装备，针对复杂多变的海洋环境开展深海高精度水下声学定位技术研究，将成为推动海洋强国建设不断取得新成就的必要手段[3]。

基于多元传感器阵列的水声定位系统通过水面声学基阵与水下声学应答器进行声波交互，通过估计声波从应答器到达基阵中心的传播时延和声波在水中的传播速度估计基阵与目标的距离，通过声波到达各基元的时延差(或相位差)估计目标方位，从而获得目标相对基阵中心的位置，再结合罗经、全球定位系统(Global Positioning System, GPS)等外接辅助设备转换得到目标的绝对位置。目前，高精度水声定位面临着定位信号长距离传播导致能量衰减，水面作业船、海洋生物以及海洋环境等因素产生的干扰(噪声)，导致接收数据的信干噪比(Signal to Interference and Noise Radio, SINR)较低，而SINR是决定时延估计精度的重要因素，从而制约定位精度进一步提高，因此抑制干扰(噪声)的影响是高精水声定位中不可避免且亟待解决的问题之一。

目前的干扰抑制算法主要分为两类，一类是针对单通道的降噪方法，如最小均方误差(Least Mean Square, LMS)自适应滤波器法[4–6]、短时傅里叶变换(Short Time Fourier Transform, STFT)[7,8]等。LMS 算法基于最小均方误差准则，通过输入量与期望响应的差值对权值进行迭代更新，以获取最优权值。LMS 算法具有计算量小、稳健性强、易于实现等优点，但该方法的收敛过程慢，而且对于随机干扰的适应性较差，而在复杂的海洋环境中干扰的统计特性往往是复杂且随机的。基于短时傅里叶变换的干扰抑制算法通过对接收数据进行短时傅里叶变换，根据期望信号和干扰在时频域的能量分布对期望信号进行重构，从而达到抑制干扰的效果，但该方法对于SINR要求较高。与此同时，LMS和STFT算法的共同问题是会影响期望信号的相位，这会严重影响定位系统的定位精度。另一类算法是针对阵列信号的算法，包括波束形成类方法和子空间类算法。波束形成类算法[9–12]通过对各基元的接收数据进行加权，从而在期望方向形成波束，抑制非期望方向的信号，可以视为一种空域滤波器。但是在进行波束设计时往往需要已知阵列流型，从而针对性地设计波束，而在水声定位作业中目标的方位通常是未知的。子空间类算法是依靠对数据矩阵的奇异值分解或对协方差矩阵的特征值分解估计信号子空间和噪声子空间，子空间类算法有3个重要理论基础：(1)信号子空间与噪声子空间垂直；(2)对于窄带信号模型，阵列流型所张成的子空间与信号子空间相等；(3)对于窄带信号模型，信号子空间维度等于信源个数。子空间类算法的典型应用为多信号分类(MUltiple SIgnal Classification, MUSIC)算法，利用信号子空间与噪声子空间垂直特性进行频率估计、方位估计等。基于子空间的特性，还有学者提出了基于子空间理论的干扰抑制方法。Bose等人[13]提出了通过子空间类算法进行语音降噪，通过对信号协方差矩阵进行特征值分解估计信号子空间，并利用信号子空间与噪声子空间的垂直特性将带噪数据线性投影到信号子空间中以实现数据降噪，该算法只考虑了信号中仅包含高斯噪声的情况，当干扰存在时，信号子空间与干扰子空间将发生空间纠缠，导致无法抑制干扰。针对干扰和噪声同时存在的情况，目前的研究内容主要集中在对窄带信号的抑制。张春海等人[14]提出基于子空间跟踪的直接序列扩频 (Direct-Sequence Spread-Spectrum, DSSS)通信系统抗窄带干扰研究，通过跟踪接收信号自相关矩阵大特征值对应特征矢量构成的干扰子空间，实现对窄带干扰的有效抑制。周峰等人[15]提出了一种用于合成孔径雷达的基于回波数据特征子空间滤波的干扰抑制方法，首先在频域对窄带干扰进行识别，然后在时域对干扰进行抑制处理。张小飞等人[16]提出一种基于斜投影的波束形成算法，算法通过构造斜投影矩阵先对接收数据进行斜投影抑制干扰和噪声的影响，然后进行波束设计，进而提高了波束形成的鲁棒性，但算法只适用于窄带模型，且要求阵列流型已知，但在水声定位过程中，信号和干扰均为宽带且目标的方位是未知的。

综上，基于子空间的干扰抑制算法目前还面临如下挑战：(1)接收数据中同时包含期望信号、干扰和噪声；(2)期望信号和干扰均为宽带；(3)期望信号入射方位未知。基于上述挑战，本文提出一种基于子空间理论的宽带强干扰抑制方法，首先通过贝叶斯信息量准则(Bayesian Information Criterion,BIC)估计信号子空间和干扰子空间的维度，然后推导不同信号假设下的概率密度函数，求解未知参数的最大似然估计，构造广义似然比并通过最优匹配广义似然比检测法估计与期望信号最匹配的子空间，然后以此构造空间投影算子对接收数据进行线性投影，实现对干扰和噪声的抑制。本方法的优点在于，适用于宽带阵列信号，且无需已知阵列流型和信源个数，同时不影响期望信号的相位信息。仿真结果表明本文所提方法能够在低SINR条件下有效抑制干扰的影响，提高水声定位系统时延估计精度。

2 阵列信号模型与子空间理论

2.1 子空间框架下的阵列信号模型

假设定位系统接收基阵由 N 个无指向性的基元构成，每个基元进行时域均匀采样接收 Lx个快拍(即采样点个数)，位于远场的目标发射快拍数为Ls的波形已 知 的 定 位 信 号s0(t) , t ∈Ωs≡{K+1,K+2,···,K+Ls}且 K+Ls

即在 H0假设下，接收信号只包含干扰信号和噪声；在 H1假设下，接收信号包含期望信号、干扰信号和噪声。其中，期望信号可以表示为

其中，[ ·]T表示转置运算。τk表示期望信号到达第k个基元和参考点之间的时延差，该时延差由入射信号方位和基元位置关系决定。

因此，包含 Lx个快拍的接收数据矩阵X ∈ℜN×Lx可以表示为

在子空间理论框架下，假设期望信号和干扰分别位于两个独立的子空间< Q >和< F>中。因此，期望信号可以表示为一个列满秩矩阵 Q ∈ℜN×r所有列的线性组合[17,18]，r为期望信号子空间的维度。于是s (t)可以表示为

其中， A(t)∈ℜr×1为 期望信号在子空间中的坐标。

同理，干扰信号可以表示为一个列满秩矩阵F ∈ℜN×q所有列的线性组合，q为干扰子空间的维度，于是i (t)可以表示为

其中，B (t)∈ℜq×1为 干扰信号在子空间< F>中的坐标。

将式(4)和式(5)代入式(1)可得子空间理论框架下基阵在t 时刻的接收数据

进一步假设干扰为随机干扰且协方差为Ri，即

其 中，σ2I为噪声的方差。

2.2 线性子空间投影

根据2.1节所述的接收信号模型，假设期望信号、干扰信号和噪声相互独立，则接收数据的协方差矩阵可表示为

其中，RA=EA(t)AT(t)。

令子空间由列满秩矩阵[ Q F]的列张成，即

其中， span{M} 表示由矩阵 M的列张成子空间。因此

其中，对角阵 P=diag{RARB}表示复合向量[A(t) B(t)] 的协方差矩阵。因此Rx可以简化为

由于 M为列满秩矩阵且期望信号和干扰信号是相互独立的，所以矩阵 P为非奇异的埃尔米特矩阵，因此可以推断 Rx一定是非负正定埃尔米特矩阵[19]。根据埃尔米特矩阵的特性，可以对 Rx进行特征值分解，则有

其 中， W=diag(λ1λ2··· λN)为 特 征 值 矩 阵，U=[u1u2··· uN]为 特征向量矩阵，当r+q

当i r+q时，{λiui} 对应噪声的特征值和特征向量。因此，Rx的特征分解可以进一步表示为

需要注意的是， UM和 M的列分别为子空间< M>的一组基，在数值上并不一定相等，这就意味着无法直接从UM中将Q 和F 分离。

如果能正确估计期望信号和干扰信号的子空间，可通过线性投影的方式实现抑制干扰和噪声的影响。一般地，子空间和< F>的正交投影算子可以表示为

综上，若已知张成信号子空间和干扰子空间的矩阵 Q和 F，可以构造对应的正交投影算子或斜投影算子，通过线性投影运算实现干扰和噪声抑制，主要有两种实现方式：

图1 子空间及其投影关系示意图

方式1 将接收数据投影在信号子空间上从而在保留期望信号的前提下，抑制部分干扰和噪声的影响

方式2 将接收数据投影在斜投影空间上从而在保留期望信号的前提下，抑制全部干扰和部分噪声的影响

为了便于讨论，下文分别将两种投影方式命名为 SP −I和 SP −II。

3 基于最大似然估计的信号子空间估计

根据式(13)—式(15)可知，决定子空间的有两个参数，即子空间维度和接收数据协方差矩阵。当干扰存在时，信号子空间和干扰子空间会出现纠缠现象，通过对接收数据协方差矩阵进行特征值分解仅能得到信号和干扰的联合子空间，为了抑制干扰的影响则需要单独估计信号子空间，本节的目的是分别估计子空间维度和信号协方差矩阵，然后构造投影算子抑制干扰和噪声。

3.1 子空间维度估计

根据文献[21]可知，入射信号确定时，信号子空间的维度与入射角度无关。因此对于定位信号已知的水声定位系统而言，可以假定任意方位的入射信号，然后直接估计信号子空间的维度。假设长度为Ls接收数据xs(t),t ∈{1,2,···,Ls}中只包含期望信号和高斯白噪声，即

其中，‖ ·‖表示矩阵(向量)的二范数。对式(29)取对数，则

其中， Tr(·)表 示矩阵的迹。因此，通过遍历r =1,2,···,N并使BIC最小即可获得信号子空间的维度。同理，在水声定位过程中，通常可以在每次接收期望信号前，采集一段仅包含干扰和噪声的测试信号xc(t)，同样可以通过上述方法估计干扰信号子空间维 度。

3.2 期望信号协方差矩阵估计

广义似然比作为最大的不变统计量，广泛用于信号检测、方位估计等领域[23]。本节推导了不同信号假设下的概率密度函数，求解出未知参数的最大似然估计，构造广义似然比并用匹配广义似然比检测法估计与期望信号最匹配的子空间。广义似然比定义为有约束条件下的概率密度函数最大值与无约束条件下概率密度函数最大值之比，首先假设信号和干扰的子空间均已知，根据式(6)所定义的信号模型，构造广义似然比检测的表达式为

根据第2节所述内容可知，在信号波形确定的情况下，信号子空间只与期望信号的入射方位有关，因此可以将Q 写成Qα,β，其中α 和 β 分别表示信号入射的俯仰角和方位角。期望信号的协方差可以表示为

4 仿真验证

本节通过仿真分析所提方法的性能并与现有方法进行对比。如图2所示，仿真采用一个均匀分布的平面阵列，基元个数为30个，相邻阵元间距为5 cm。每个定位周期的采样信号长度为25 ms，采样频率为100 kHz。定位信号为线性调频信号，频带宽度为9～kHz，信号长度为5 ms。信号入射的俯仰角和方位角分别为75°和30°。干扰为随机干扰，入射的俯仰角和方位角分别为10°和60°。噪声为高斯白噪声。

互相关时延估计法是水声定位中最常用且最有效的时延估计算法，然而互相关时延估计法的精度受到干扰和噪声的影响，本节首先通过仿真SINR对互相关时延估计精度的影响。图3展示了不同的干扰噪声比(Interference to Noise Radio, INR)条件下互相关时延估计精度随SINR的变化。

根据图3可以发现，互相关时延估计精度会随着SINR的减小而降低；在相同的SINR条件下，INR越小时延估计误差越大，原因是干扰与期望信号是非相干的，当干扰为主要分量时，互相关时延估计精度要优于噪声为主要分量的情况。接下来验证本文所提出的方法对于时延估计精度的提升效果。按上述仿真条件，得到SINR为0 dB，INR为0 dB的基阵接收数据，图4展示了1号基元的接收数据。

根据图4所示的接收数据估计期望信号子空间，图5展示了在角度空间内的搜索结果，其中颜色对应广义似然比的大小。为了便于分析，分别提取峰值位置对应的俯仰角和方位角方向的切面，如图6所示，结果表明本文所提方法能够精确估计信号子空间。

接下来，根据式(47)的估计结果，分别构造正交投影算子和斜投影算子，将接收数据进行对应的线性投影，然后通过互相关法估计各基元接收信号的时延，并与LMS和STFT算法的结果进行对比。当INR=0 dB时，经过1000次蒙特卡罗仿真，各种方法输出数据的时延估计均方根误差(Root Mean Squared Error, RMSE)随SINR的变化如图7所示。

为了进一步验证算法在不同INR情况下的性能，分别对INR=20 dB和INR=50 dB两种情况进行仿真，经过1000次蒙特卡罗仿真，上述两种情况的时延估计的RMSE随SINR的变化分别如图8和图9所示。

图2 定位系统接收基阵阵型图

图3 相关时延估计随SINR的变化

图4 基阵1号基元接收信号

图5 广义似然比在角度空间内的搜索结果

图6 广义似然比随入射角度变化图

图7 时延估计误差随SINR变化图，INR=0 dB

图8 时延估计误差随SINR变化图，INR=20 dB

图9 时延估计误差随SINR变化图，INR=50 dB

根据图7—图9的仿真结果可以发现，当INR较大即非期望信号中干扰信号为主要分量时，基于斜投影的降噪方法( SP −II)的性能明显优于其他方法；当INR较小即非期望信号中高斯噪声为主要分量时，基于正交投影的降噪方法(S P −I)明显优于其他方法；随着INR的降低， SP −I 和S P −II的性能均会下降，但是 SP −I对于噪声的敏感程度低于SP −II。上述现象产生的原因是正交投影和斜投影的投影方式导致的，斜投影沿着干扰子空间的方向将期望信号投影到信号子空间，因此能够最大限度地消除干扰，而正交投影仅将期望信号进行正交投影，当干扰子空间与信号子空间非正交时，会有部分干扰信号分量投影到信号子空间中。当噪声分量较高时，会对干扰子空间的估计产生影响，进而影响斜投影的性能，而正交投影不受干扰子空间的影响。

5 结论

本文针对强干扰降低水声定位系统时延估计精度的问题，提出一种基于子空间理论的宽带强干扰抑制方法，通过估计期望信号子空间和干扰子空间，构造投影算子并对接收数据进行线性投影，从而抑制干扰和噪声对定位系统时延估计精度的影响。相比传统方法，本文所提方法可适用于宽带阵列信号，且无需已知阵列流型和信源个数，同时不影响期望信号的相位信息。仿真结果显示，本文所述的方法能有效抑制宽带强干扰的影响，提高系统时延估计精度。时延估计误差的仿真结果显示，当非期望信号中干扰信号为主要分量时，基于斜投影的降噪方法性能最优；当非期望信号中高斯噪声为主要分量时，基于正交投影的降噪方法最优。