近年，新型无人水下潜航器(unmanned underwater vehicle, UUV)作为高技术的水下无人作战平台，在复杂海洋环境中工作时具备低消耗、隐蔽性、机动性等特点，已在水下目标侦察与防御、海洋科学考察与试验等领域发挥了重要作用^[1]，是遂行水下无人攻防作战任务的重要保障之一。但目前水下目标隐身技术迅速发展，使得潜艇和水面舰艇的辐射噪声水平不断降低，从而导致水下航行器对其探测难度不断提高，故现阶段迫切需要大幅度提高UUV的目标探测距离，以保证对远场目标的发现、跟踪以及精确打击。提高UUV探测距离的有效手段之一是降低工作频率，但工作频率的降低受UUV传统头部平面基阵几何尺寸限制，并且会引起空间增益和估计精度的下降，故往往会在其舷侧加装舷侧阵，将若干个声呐纵向排列形成两列长线阵，这种布放形式不占航行器空间位置，也不影响其水动力学性能，同时能够增大水声探测系统基阵有效孔径并尽量降低工作频率^[2]。

由于舷侧阵直接安装在航行器壳体上，平台本身的近场干扰很强，主要包括平台振动引起的噪声、流噪声和螺旋桨辐射噪声等^[3], 这些辐射噪声的频率接近低频声呐的工作频率，故无法用简单的滤波算法来消除，会降低舷侧阵声呐系统的作用距离、目标识别概率及精度。由于舷侧阵基阵的长度大约占艇长1/3~2/3，平台振动引起的噪声就是由艇内安装的众多机械设备运行时对艇体产生的激励带来的，会对目标信号构成本地强干扰，在航行器低速航行时影响最为严重。平台振动引起的噪声主要包含2个部分^[4]: 1)由航行器内部机械振动通过壳体、管道、基阵安装结构传至各水听器，这部分噪声可以通过在声呐与舷侧之间安装减振隔振结构降低其影响；2)由航行器内部机械振动在其周围激起声波，声波直接传至各水听器, 它属于近场辐射噪声，本文针对的就是这部分噪声。

对于平台振动噪声，国内外现有的处理方式大多是采用自适应噪声抵消^[5]、自适应波束形成^[6]或盲信号分离^[7]等基于平面波传播的抑制或抵消的方法；文献[4]则是根据平台噪声特性选用具有更低旁瓣的窗函数进行数据截断来提高系统的信号检测能力。但由于振动干扰传播信道和产生机理十分复杂，上述处理方式的效果都很有限，不能有效提高声呐目标方位估计的性能。空域矩阵滤波技术^[8-9]是一种新兴的对阵列信号进行预处理的技术，通过对空间方位设置通带和阻带，得到一个空域矩阵与阵列接收数据相乘，来保证通带区域内目标信号最大程度无失真通过的同时，抑制阻带区域的强干扰，能获得更高精度的平面波方位估计，为平台近场辐射噪声抑制提供了可行途径。由于空域矩阵滤波器滤波后的数据仍为阵元域，所以此类滤波器可以用在很多目标定向或定位算法之前^[9]。

本文将空域矩阵滤波技术、矢量信号处理技术以及平面波方位估计技术结合，共同设计抑制平台强噪声干扰的空域矩阵滤波器，给出空域矩阵滤波器的最优解及推导过程，并基于矢量舷侧阵的最小方差无失真响应(minimum variance distortionless response，MVDR)波束形成算法对滤波前后的远场弱目标信号进行波达方位(direction of arrival，DOA)估计的仿真，对比分析得到滤波器性能。

1 矢量舷侧阵信号处理模型

对于矢量舷侧阵来说，被探测目标信号源距阵列的距离D≥L²/λ，满足远场条件，辐射的声信号经过远距离传播到达基阵后能量很小，为远场弱信号，而平台辐射噪声与舷侧阵声呐基阵的距离较近，能量相对远场目标来说比较强，应该视为近场强干扰^[10]。为便于分析，将舷侧阵抽象为均匀线列阵，取航行器左右舷中任意一侧的舷侧阵进行信号处理模型的建立。假设P个水下目标、Q个近场辐射干扰源(为点干扰)、基阵位于同一水平面，矢量水听器阵是由M个水声换能器构成，各个阵元等间距布放，间距为d，阵的总长度为(M-1)d，建立如图 1所示二维UUV舷侧阵信号处理模型。

则接收阵列信号模型可用矩阵形式表示为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{X}}(t,f) = {\mathit{\boldsymbol{A}}_s}(\mathit{\boldsymbol{\theta }},f){\mathit{\boldsymbol{S}}_s}(t,f) + {\mathit{\boldsymbol{A}}_n}(\mathit{\boldsymbol{\varphi }},f){\mathit{\boldsymbol{S}}_n}(t,f) + }\\ {\mathit{\boldsymbol{N}}(t,f)} \end{array} $

(1)

式中：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{X}}(t,f) = {{\left[ {{\mathit{\boldsymbol{X}}_1}(t,f) \cdots {\mathit{\boldsymbol{X}}_M}(t,f)} \right]}^{\rm{T}}}}\\ {{\mathit{\boldsymbol{S}}_s}(t,f) = {{\left[ {{\mathit{\boldsymbol{S}}_{s1}}(t,f) \cdots {\mathit{\boldsymbol{S}}_{sP}}(t,f)} \right]}^{\rm{T}}}}\\ {{\mathit{\boldsymbol{S}}_n}(t,f) = {{\left[ {{\mathit{\boldsymbol{S}}_{n1}}(t,f) \cdots {\mathit{\boldsymbol{S}}_{nQ}}(t,f)} \right]}^{\rm{T}}}}\\ {{\mathit{\boldsymbol{A}}_s}(\mathit{\boldsymbol{\theta }},f) = \left[ {{\mathit{\boldsymbol{a}}_s}\left( {{\theta _1},f} \right) \cdots {\mathit{\boldsymbol{a}}_s}\left( {{\theta _p},f} \right) \cdots {\mathit{\boldsymbol{a}}_s}\left( {{\theta _P},f} \right)} \right]}\\ {{\mathit{\boldsymbol{A}}_n}(\mathit{\boldsymbol{\varphi }},f) = \left[ {{\mathit{\boldsymbol{a}}_n}\left( {{\varphi _1},f} \right) \cdots {\mathit{\boldsymbol{a}}_n}\left( {{\varphi _q},f} \right) \cdots {\mathit{\boldsymbol{a}}_n}\left( {{\varphi _Q},f} \right)} \right]}\\ {\mathit{\boldsymbol{N}}(t,f) = {{\left[ {{\mathit{\boldsymbol{n}}_1}(t,f) \cdots {\mathit{\boldsymbol{n}}_M}(t,f)} \right]}^{\rm{T}}}} \end{array}} \right. $

(2)

式中：S_s(t, f)、S_n(t, f)分别为目标信号源矩阵和平台噪声干扰源矩阵；A_s(θ, f)、A_n(θ, f)为远场平面波阵列流型矩阵和近场球面波阵列流型矩阵; N(t, f)为加性高斯环境噪声矩阵；θ=[θ₁…θ_P]、φ=[φ₁φ₂…φ_Q]分别为目标信号与噪声的入射角度。信号源和干扰源到接收阵的方向向量：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{a}}_s}\left( {{\theta _p},f} \right) = \left[ {1,\exp \left( { - {\rm{j}}2{\rm{ \mathsf{ π} }}fd\sin {\theta _p}/{c_z}} \right) \cdots } \right.\\ {\left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \exp \left( { - {\rm{j}}2{\rm{ \mathsf{ π} }}f(M - 1)d\sin {\theta _p}/{c_z}} \right)} \right]^{\rm{T}}}\\ {\mathit{\boldsymbol{a}}_n}\left( {{\varphi _q},f} \right) = \left[ {1,\exp \left( { - {\rm{j}}2{\rm{ \mathsf{ π} }}f\left( {{r_{q,1}} - {r_{q,0}}} \right)/{c_n}} \right) \cdots } \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \exp {\left( { - {\rm{j}}2{\rm{ \mathsf{ π} }}f\left( {{r_{q,M - 1}} - {r_{q,0}}} \right)/{c_n}} \right]^{\rm{T}}} \end{array} \right. $

其中：

$ \left\{ {\begin{array}{*{20}{l}} {{r_{q,0}} = \left| {h/\cos {\varphi _q}} \right|}\\ {{r_{q,1}} = \sqrt {r_{q,0}^2 + {d^2} - 2{r_{q,0}}d\sin {\varphi _q}} }\\ {{r_{q,M - 1}} = \sqrt {r_{q,0}^2 + {{[(M - 1)d]}^2} - 2{r_{q,0}}(M - 1)d\sin {\varphi _q}} } \end{array}} \right. $

式中c_z代表水下声速，c_z=1 500 m/s；c_n代表空气中声速，c_n=340 m/s。

由于声压与振速的关系：

$ \mathit{\boldsymbol{v}}(t) = \left[ {\begin{array}{*{20}{l}} {{v_y}(t)}\\ {{v_y}(t)} \end{array}} \right] = p(t)\left[ {\begin{array}{*{20}{l}} {\sin \theta }\\ {\cos \theta } \end{array}} \right] $

此时声矢量传感器的阵列流形为：

$ \mathit{\boldsymbol{u}}(\theta ) = \left[ {\begin{array}{*{20}{c}} 1\\ {\sin \theta }\\ {\cos \theta } \end{array}} \right] $

代入式(1)中得到矢量阵接收阵列信号模型：

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}_v}(t,f) = {\mathit{\boldsymbol{A}}_{sv}}(\theta ,f){\mathit{\boldsymbol{S}}_s}(t,f) + }\\ {{\mathit{\boldsymbol{A}}_{nv}}(\mathit{\boldsymbol{\varphi }},f){\mathit{\boldsymbol{S}}_n}(t,f) + {\mathit{\boldsymbol{N}}_v}(t,f)} \end{array} $

(3)

其中的阵列流型矩阵和对应的方向向量变为：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{A}}_{sv}}(\mathit{\boldsymbol{\theta }},f) = \left[ {{\mathit{\boldsymbol{a}}_{sv}}\left( {{\theta _1},f} \right) \cdots {\mathit{\boldsymbol{a}}_{sv}}\left( {{\theta _p},f} \right) \cdots {\mathit{\boldsymbol{a}}_{sv}}\left( {{\theta _P},f} \right)} \right]}\\ {{\mathit{\boldsymbol{A}}_{nv}}(\mathit{\boldsymbol{\varphi }},f) = \left[ {{\mathit{\boldsymbol{a}}_{nv}}\left( {{\varphi _1},f} \right) \cdots {\mathit{\boldsymbol{a}}_{nv}}\left( {{\varphi _q},f} \right) \cdots {\mathit{\boldsymbol{a}}_{nv}}\left( {{\varphi _Q},f} \right)} \right]}\\ {{\mathit{\boldsymbol{a}}_{sv}}{{\left( {{\theta _p},f} \right)}_{3M \times 1}} = {\mathit{\boldsymbol{a}}_s}\left( {{\theta _p},f} \right) \otimes \mathit{\boldsymbol{u}}\left( {{\theta _p}} \right)}\\ {{\mathit{\boldsymbol{a}}_{nv}}{{\left( {{\varphi _q},f} \right)}_{3M \times 1}} = {\mathit{\boldsymbol{a}}_n}\left( {{\varphi _q},f} \right) \otimes \mathit{\boldsymbol{u}}\left( {{\varphi _q}} \right)} \end{array}} \right. $

式中⊗为Kronecker积的运算符。

2 空域矩阵滤波器设计

首先针对声压阵接收数据设计中心频率为f的空域矩阵滤波器H(f)_M×M，再推广到矢量阵。滤波后的输出为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Y}}(t,f) = \mathit{\boldsymbol{H}}(f)\mathit{\boldsymbol{X}}(t,f) = \mathit{\boldsymbol{H}}(f){\mathit{\boldsymbol{A}}_s}({\bf{ \pmb{\mathsf{ θ}} }},f){\mathit{\boldsymbol{S}}_s}(t,f) + }\\ {\mathit{\boldsymbol{H}}(f){\mathit{\boldsymbol{A}}_n}(\mathit{\boldsymbol{\varphi }},f){\mathit{\boldsymbol{S}}_n}(t,f) + \mathit{\boldsymbol{H}}(f)\mathit{\boldsymbol{N}}(t,f)} \end{array} $

(4)

空域矩阵滤波器是通过作用于信号方向向量来对其产生通过或抑制效果。定义‖H(f)a(θ_i, f)‖_F²为向量的响应，‖H(f)a(θ_i, f)-a(θ_i, f)‖_F²为向量的响应误差，‖·‖_F代表Frobenius范数。

滤波器的设计需将全空间划分成通带区域Θ_P和阻带区域Θ_S，其分别对应的阵列流型矩阵为：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{V}}_P}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_P}} \right) = \left[ {\mathit{\boldsymbol{a}}\left( {{\theta _1}} \right) \cdots \mathit{\boldsymbol{a}}\left( {{\theta _p}} \right) \cdots \mathit{\boldsymbol{a}}\left( {{\theta _P}} \right)} \right],\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _p} \in {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_P},p = 1, \cdots ,P\\ {\mathit{\boldsymbol{V}}_S}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_S}} \right) = \left[ {\mathit{\boldsymbol{a}}\left( {{\theta _1}} \right) \cdots \mathit{\boldsymbol{a}}\left( {{\theta _s}} \right) \cdots \mathit{\boldsymbol{a}}\left( {{\theta _S}} \right)} \right],\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _s} \in {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_S},s = 1,2, \cdots ,S \end{array} \right. $

式中：a(θ_p)和a(θ_s)分别代表通、阻带离散化后的方向向量；P和S分别为通、阻带离散化的点数。

由于平台强干扰来自于某确定的小范围方位内，因此可以选择对窄阻带设置零响应约束，求解阻带零点约束的空域矩阵滤波器来消除近场干扰的影响。故建立如下的最优化问题^[11]：

$ \left\{ \begin{array}{l} \mathop {\min }\limits_H \left\| {\mathit{\boldsymbol{H}}(f){\mathit{\boldsymbol{V}}_P}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_P}} \right) - {\mathit{\boldsymbol{V}}_P}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_P}} \right)} \right\|_{\rm{F}}^2\\ {\rm{subject}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{to}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{H}}(f){\mathit{\boldsymbol{V}}_S}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_S}} \right) = {{\bf{0}}_{M \times S}} \end{array} \right. $

其中约束条件H(f)V_S(Θ_S)=0_M×S可以表示为：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathop{\rm Re}\nolimits} \left[ {\mathit{\boldsymbol{H}}(f)\mathit{\boldsymbol{a}}\left( {{\theta _s}} \right)} \right] = {{\bf{0}}_{M \times 1}},s = 1,2, \cdots ,S}\\ {{\mathop{\rm Im}\nolimits} \left[ {\mathit{\boldsymbol{H}}(f)\mathit{\boldsymbol{a}}\left( {{\theta _s}} \right)} \right] = {{\bf{0}}_{M \times 1}},s = 1,2, \cdots ,S} \end{array}} \right. $

因此最优化问题转化为：

$ \left\{ {\begin{array}{*{20}{l}} {\mathop {\min }\limits_H \left\| {\mathit{\boldsymbol{H}}(f){\mathit{\boldsymbol{V}}_P}\left( {{\mathit{\boldsymbol{\theta }}_P}} \right) - {\mathit{\boldsymbol{V}}_P}\left( {{\mathit{\boldsymbol{\theta }}_P}} \right)} \right\|_{\rm{F}}^2}\\ {{\rm{subject}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{to }}\left\{ {\begin{array}{*{20}{l}} {{\mathop{\rm Re}\nolimits} \left[ {\mathit{\boldsymbol{H}}(f)\mathit{\boldsymbol{a}}\left( {{\theta _s}} \right)} \right] = {\bf{0}},s = 1,2, \cdots ,S}\\ {{\mathop{\rm Im}\nolimits} \left[ {\mathit{\boldsymbol{H}}(f)\mathit{\boldsymbol{a}}\left( {{\theta _s}} \right)} \right] = {\bf{0}},s = 1,2, \cdots ,S} \end{array}} \right.} \end{array}} \right. $

(5)

为方便求解，将H(f)、V_P(Θ_P)、V_S(Θ_S)简记为H、V_P、V_S，通过构造Lagrange函数的方法解上述最优化问题，式(5)对应的Lagrange函数为：

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} L\left( {\mathit{\boldsymbol{H}},{\mathit{\boldsymbol{\lambda }}_1}, \cdots ,{\mathit{\boldsymbol{\lambda }}_s},{\mathit{\boldsymbol{\delta }}_1}, \cdots ,{\mathit{\boldsymbol{\delta }}_S}} \right) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\| {\mathit{\boldsymbol{H}}{\mathit{\boldsymbol{V}}_P} - {V_P}} \right\|_{\rm{F}}^2 - \left\{ {\mathit{\boldsymbol{\lambda }}_1^{\rm{T}}{\mathop{\rm Re}\nolimits} \left[ {\mathit{\boldsymbol{Ha}}\left( {{\theta _1}} \right]} \right) + \cdots + } \right.\\ \left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{\lambda }}_{\rm{S}}^{\rm{T}}{\mathop{\rm Re}\nolimits} \left[ {\mathit{\boldsymbol{Ha}}\left( {{\theta _S}} \right)} \right]} \right\} - \left\{ {\mathit{\boldsymbol{\delta }}_1^{\rm{T}}{\mathop{\rm Im}\nolimits} \left[ {\mathit{\boldsymbol{Ha}}\left( {{\theta _1}} \right)} \right] + \cdots + } \right.\\ \left. {\mathit{\boldsymbol{\delta }}_{\rm{S}}^{\rm{T}}{\mathop{\rm Im}\nolimits} \left[ {\mathit{\boldsymbol{H}}\left( {{\theta _S}} \right)} \right]} \right\} = {\mathop{\rm tr}\nolimits} \left( {\mathit{\boldsymbol{H}}{\mathit{\boldsymbol{V}}_P}\mathit{\boldsymbol{V}}_{\rm{P}}^{\rm{H}}{\mathit{\boldsymbol{H}}^{\rm{H}}}} \right) - {\mathop{\rm tr}\nolimits} \left( {\mathit{\boldsymbol{H}}{\mathit{\boldsymbol{V}}_{\rm{P}}}\mathit{\boldsymbol{V}}_{\rm{P}}^{\rm{H}}} \right) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{c}} {{\mathop{\rm tr}\nolimits} \left( {{\mathit{\boldsymbol{V}}_{\rm{P}}}\mathit{\boldsymbol{V}}_{\rm{P}}^{\rm{H}}{\mathit{\boldsymbol{H}}^{\rm{H}}}} \right) + {\mathop{\rm tr}\nolimits} \left( {{\mathit{\boldsymbol{V}}_{\rm{P}}}\mathit{\boldsymbol{V}}_{\rm{P}}^{\rm{H}}} \right) - }\\ \sum\limits_{s = 1}^S {\left( {\frac{{\mathit{\boldsymbol{\lambda }}_{\rm{s}}^{\rm{T}}}}{2} + \frac{{\mathit{\boldsymbol{\delta }}_{\rm{s}}^{\rm{T}}}}{{2{\rm{j}}}}} \right)\mathit{\boldsymbol{Ha}}\left( {{\theta _s}} \right) - \sum\limits_{s = 1}^S {\left( {\frac{{\mathit{\boldsymbol{\lambda }}_{\rm{s}}^{\rm{T}}}}{2} + \frac{{\mathit{\boldsymbol{\delta }}_{\rm{s}}^{\rm{T}}}}{{2{\rm{j}}}}} \right)} {\mathit{\boldsymbol{H}}^*}\mathit{\boldsymbol{a}}{{\left( {{\theta _s}} \right)}^*}} \end{array} \end{array} $

(6)

式中，λ和δ均为Lagrange乘子，λ=[λ₁ … λ_S], δ=[δ₁ … δ_S]。

由于目标函数和约束条件都是严格的凸函数，故对于任意λ和δ，对式(6)Lagrange函数求偏导等于零的解即为最优化问题的全局最优解，即：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{\partial L}}{{\partial \mathit{\boldsymbol{H}}}} = {\bf{0}} \to {{\mathit{\boldsymbol{\hat H}}}^*}{\mathit{\boldsymbol{V}}_P}\mathit{\boldsymbol{V}}_{\rm{P}}^{\rm{T}} - \mathit{\boldsymbol{V}}_{\rm{P}}^*\mathit{\boldsymbol{V}}_{\rm{P}}^{\rm{T}} - \mathit{\boldsymbol{\hat \gamma V}}_{\rm{S}}^{\rm{T}} = {\bf{0}}}\\ {\frac{{\partial L}}{{\partial {\mathit{\boldsymbol{\lambda }}_s}}} = {\bf{0}} \to {\mathop{\rm Re}\nolimits} \left[ {\mathit{\boldsymbol{Ha}}\left( {{\theta _s}} \right)} \right] = {\bf{0}},s = 1,2, \cdots ,S}\\ {\frac{{\partial L}}{{\partial {\mathit{\boldsymbol{\delta }}_s}}} = {\bf{0}} \to {\mathop{\rm Im}\nolimits} \left[ {\mathit{\boldsymbol{Ha}}\left( {{\theta _s}} \right)} \right] = {\bf{0}},s = 1,2, \cdots ,S} \end{array}} \right. $

(7)

式中$\hat{\boldsymbol{\gamma}}=\hat{\boldsymbol{\lambda}} / 2+\hat{\boldsymbol{\delta}} / 2 \rm{j}$，对式(7)进行求解得:

$ \mathit{\boldsymbol{\hat H}} = {\mathit{\boldsymbol{I}}_{M \times M}} - {\mathit{\boldsymbol{V}}_S}{\left[ {\mathit{\boldsymbol{V}}_S^{\rm{H}}{{\left( {{\mathit{\boldsymbol{V}}_P}\mathit{\boldsymbol{V}}_P^{\rm{H}}} \right)}^{ - 1}}{\mathit{\boldsymbol{V}}_S}} \right]^{ - 1}}\mathit{\boldsymbol{V}}_S^{\rm{H}}{\left( {{\mathit{\boldsymbol{V}}_P}\mathit{\boldsymbol{V}}_P^{\rm{H}}} \right)^{ - 1}} $

(8)

式中I为M×M维单位向量。

对于矢量阵，结合矢量信号处理的方式可以得到通阻带阵列流型矩阵的矢量表达形式V_Pv和V_Sv：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{V}}_{Pv}} = {{\left[ {{\mathit{\boldsymbol{V}}_{pv}}\left( {{\theta _1}} \right) \cdots {\mathit{\boldsymbol{V}}_{pv}}\left( {{\theta _p}} \right) \cdots {\mathit{\boldsymbol{V}}_{pv}}\left( {{\theta _P}} \right)} \right]}_{3M \times P}}}\\ {{\mathit{\boldsymbol{V}}_{Sv}} = {{\left[ {{\mathit{\boldsymbol{V}}_{sv}}\left( {{\theta _1}} \right) \cdots {\mathit{\boldsymbol{V}}_{sv}}\left( {{\theta _s}} \right) \cdots {\mathit{\boldsymbol{V}}_{sv}}\left( {{\theta _S}} \right)} \right]}_{3M \times S}}} \end{array}} \right. $

式中：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{V}}_{pv}}\left( {{\theta _p}} \right) = \mathit{\boldsymbol{a}}\left( {{\theta _p}} \right) \otimes \mathit{\boldsymbol{u}}\left( {{\theta _p}} \right),{\theta _p} \in {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_P}}\\ {{\mathit{\boldsymbol{V}}_{sv}}\left( {{\theta _s}} \right) = \mathit{\boldsymbol{a}}\left( {{\theta _s}} \right) \otimes \mathit{\boldsymbol{u}}\left( {{\theta _s}} \right),{\theta _s} \in {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_S}}\\ {\mathit{\boldsymbol{u}}\left( {{\theta _i}} \right) = {{\left[ {\begin{array}{*{20}{l}} 1&{\sin {\theta _i}}&{\cos {\theta _i}} \end{array}} \right]}^{\rm{T}}}} \end{array}} \right. $

因此对于式(8)中的滤波器的矢量表达式为：

$ {\mathit{\boldsymbol{\hat H}}_v} = {\mathit{\boldsymbol{I}}_{3M \times 3M}} - {\mathit{\boldsymbol{V}}_{Sv}}{\left[ {\mathit{\boldsymbol{V}}_{Sv}^{\rm{H}}{{\left( {{\mathit{\boldsymbol{V}}_{Pv}}\mathit{\boldsymbol{V}}_{Pv}^{\rm{H}}} \right)}^{ - 1}}{\mathit{\boldsymbol{V}}_{Sv}}} \right]^{ - 1}}\mathit{\boldsymbol{V}}_{Sv}^{\rm{H}}{\left( {{\mathit{\boldsymbol{V}}_{Pv}}\mathit{\boldsymbol{V}}_{Pv}^{\rm{H}}} \right)^{ - 1}} $

(9)

则最终滤波输出为：

$ \begin{array}{l} \mathit{\boldsymbol{Y}}(t,f) = {\mathit{\boldsymbol{H}}_v}(f){\mathit{\boldsymbol{X}}_v}(t,f) = {\mathit{\boldsymbol{H}}_v}(f){\mathit{\boldsymbol{A}}_{sv}}(\theta ,f){\mathit{\boldsymbol{S}}_s}(t,f) + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{H}}_v}(f){\mathit{\boldsymbol{A}}_{nv}}(\varphi ,f){\mathit{\boldsymbol{S}}_n}(t,f) + {\mathit{\boldsymbol{H}}_v}(f){\mathit{\boldsymbol{N}}_v}(t,f) \end{array} $

(10)

3 平台自噪声抑制算法仿真

首先分别给出空域矩阵滤波器作用前后基阵近场和远场方向向量的归一化指向性图，通过指向性对比图来验证滤波器的作用。

由于UUV尺度有限，限制了其舷侧阵的孔径。现在使用的UUV长度一般不会超过15 m。故考虑一个由16个水听器组成的舷侧阵基阵，信号为1 200 Hz的窄带信号，阵元间距为半波长，阵总长9.375 m。平台噪声干扰源距阵列1 m，为部分连续谱加线谱信号；目标信号源为满足远场条件的平面波信号。滤波器设计通带Θ_P=[-180°, 120°]∪[140°, 180°]，阻带Θ_S=[120°, 140°]。

从图 2中可以看出，经过滤波器处理之后，基阵近场方向向量的指向性幅度在设定的阻带Θ_S=[120°, 140°]内大幅下降，即空域矩阵滤波器对设定阻带内的近场干扰噪声有明显的抑制作用；而对于基阵远场方向向量的指向性图，从图 3中可以看出在设定的通带范围内是几乎没有变化的，即可以保证通带范围内的远场信号无失真通过。这里的指向性图仿真采用的是声压阵，故存在左右舷模糊，即关于90°/-90°的指向性对称，后续采用矢量阵即可消除左右舷模糊。

目前的声呐系统大多数情况下都采用波束形成算法进行信号检测，对于阵列和阵元布放形式固定的舷侧阵，最小方差无失真响应(minimum variance distortionless response，MVDR)波束形成算法较为常用且可靠。由于空域矩阵滤波器滤波后的数据仍为阵元域，故可直接用于后续MVDR算法。MVDR波束形成器的输出功率表达式为：

$ {P_{{\rm{MVDR}}}} = \frac{1}{{{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ){\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{a}}(\theta )}} $

(11)

信号的波达方位可通过搜索上式的峰值位置来确定。R=E[YY^H]为接收数据协方差矩阵，当采样快拍数N→∞时，可采用极大似然估计下的输出数据外积的统计平均值来代替：

$ \mathit{\boldsymbol{\hat R}} = \frac{1}{N}\sum\limits_i^N \mathit{\boldsymbol{Y}} {\mathit{\boldsymbol{Y}}^{\rm{H}}} $

(12)

在前述基阵指向性仿真条件的基础上，采用矢量水听器，设置远场目标信号、近场平台干扰噪声以及环境背景噪声，进行MVDR算法DOA估计的仿真。设目标信号的入射方向为25°，信噪比SNR=10 dB，平台干扰入射方向130°，干信比SIR分别为-100、-140、-150、-170 dB。图 4给出了仿真结果。

从图 4可以看出，当SIR=-140 dB时，滤波前的MVDR算法勉强能够定位出信号的来波方向，但由于近场平台噪声的强干扰，会在非目标来波方位出现幅度较高的假目标峰，而滤波后再进行MVDR算法处理的结果仅在目标方位出现单峰，消除了平台噪声强干扰带来的影响；之后随着平台噪声的增强，SIR逐渐降低，滤波前的MVDR算法受噪声干扰的影响愈加严重，显然已丢失目标方位，无法检测，而滤波后的MVDR算法依然具有非常好的信号检测和估计能力，但性能逐渐有所下降；当SIR降低至-170 dB时，可以看出空域矩阵滤波器的能力达到下限，滤波前后MVDR算法均失效。

在此利用最高峰的幅值(dB^-1)与K个次高峰平均幅值的差值来得到输出信干比，等价于：

$ {\rm{SI}}{{\rm{R}}_{{\rm{out}}}} = 10\lg \left[ {P\left( {{\theta _i}} \right)/\left( {\frac{1}{K}\sum\limits_k P \left( {{\theta _k}} \right)} \right)} \right] $

(13)

式中:P(θ_i)为波峰最高的疑似目标方位的波束输出功率；P(θ_k)为波峰高度仅次于θ_i的K个干扰目标方位的波束输出功率；[∑P(θ_k)/K]为平均干扰功率，当SIR_out大于某一阈值时，便认为检测到目标信号。在此设阈值为7 dB，K取为3，图 5给出了通过200次蒙特卡洛仿真计算得到的输入与输出信干比之间的关系曲线。

为了得到空域矩阵滤波后MVDR算法的方位分辨能力，考虑临近双目标信号入射的情况。仿真设置SIR=-140 dB，目标信号方位θ₁=25°, θ₂=28°，30°，35°，40°。如图 6给出了仿真结果。

从图 6可以看出，空域滤波处理后不仅能消除平台噪声强干扰带来的影响从而检测出正确的目标信号方位，而且具有较强的方位分辨能力，能够分辨出2个方位很相近的声目标，分辨极限为5°。

4 结论

1) 设计的空域矩阵滤波器通过作用于信号方向向量对通阻带产生了预期的无失真通过或抑制的效果。

2) 空域矩阵滤波器作用于阵列数据后再用于MVDR波束形成算法，平台自噪声带来的强干扰能够被很好的抑制，输出信噪比以及信号方位分辨能力得到了有效的提高。

3) 该方法为后置信号处理创造了条件，进而应该能够提高水下航行器舷侧阵声呐系统的工作性能及作用距离。