引言

矢量阵不仅可测量声场中的声压信息, 还能同步共点测量质点振速信息^[1,2]。在相同孔径和阵元数情况下, 利用声压和振速信息, 矢量阵比声压阵具有更高的分辨率和测向精度^[3-5], 并能解决线阵的左右舷模糊问题^[6-9]。然而, 矢量阵中的声压和振速是两个不同的物理量。已有研究证明, 即使不考虑矢量阵的幅相误差, 在各向同性噪声场中, 矢量水听器声压通道和振速通道接收的噪声功率也并不一致^[10]。声压通道和噪声通道的噪声不一致会引发虚源^[11], 从而导致常规高分辨测向方法性能急剧下降。

针对矢量阵高分辨测向问题, 一些学者将标量传感器阵列的波达方向(DOA)估计方法扩展到矢量阵列, 如最小方差无失真响应(MVDR)波束形成^[12]、子空间类方法^[13]和稀疏类方法^[14-15]等。这些方法均假设声压通道和振速通道的噪声功率是一致的。然而, 直接采用均匀噪声模型与实际矢量阵模型并不匹配。对于矢量阵通道噪声功率不一致问题, 一些针对性方法相继被提出。文献[10]建立了声压通道与矢量通道噪声不一致的矢量阵模型, 并将波束形成和Capon方向估计器扩展到矢量阵。文献[16]提出了一种噪声功率估计方法, 要求阵元数不小于信源数的3倍。然而, 上述两种方法不能满足阵元数较少的矢量阵测向的分辨率需求。文献[11]提出了一种基于子空间理论的增广子空间多重信号分类(MUSIC)方法(简称AS-MUSIC方法)。此方法将声压通道和振速通道中不一致的噪声归于信号子空间, 取得优于传统MUSIC方法的方位分辨率和估计精度, 但此方法仍然建立在MUSIC方法基础上, 理论上它的高分辨性能弱于稀疏恢复类方法。

近年来, 利用信号在空间域的稀疏性, 大量基于稀疏恢复的测向方法相继提出。与子空间类方法相比, 稀疏恢复类方法在快拍数不足、信噪比低、少阵元数和相关源等方面具有广阔的应用前景。其中, 基于奇异值分解(SVD)的L₁-SVD算法^[17]采用SVD技术降低阵列输出维数, 从而减少计算量, 提高DOA估计性能。其他稀疏恢复类方法, 如基于协方差的稀疏迭代估计(SPICE)算法^[18]和基于阵列协方差向量的稀疏表示(SRACV)算法^[19]也采用了L₁范数惩罚, 但不同于L₁-SVD算法, 它们利用协方差矩阵的统计特性, 通过对协方差矩阵进行约束来实现小快拍下的高分辨测向。文献[20]提出了一种加权L_2,1最小化方法, 该方法利用噪声子空间和过完备基之间的关系设计权值来增强解的稀疏性。相比较而言, 稀疏表示类的方法复杂度要高于子空间类方法, 但是这类方法的多目标方位分辨性能更优。

受子空间理论^[11]以及重加权L₁范数最小化思想^[20,21]的启发, 本文提出了一种应用加权子空间稀疏恢复的水声矢量阵测向方法(UAVSA-WSSR)。通过矢量阵协方差矩阵和声压通道协方差矩阵的特性, 分别估计声压通道和振速通道的噪声功率。然后, 利用估计的噪声功率构建加权矩阵, 对矢量阵接收数据进行加权, 克服矢量阵声压通道和振速通道噪声功率差异对高分辨方法的负面影响。此外, 利用加权子空间拟合和重加权L₁范数最小化将方位估计问题转化为概率参数约束的稀疏恢复问题, 以实现水声矢量阵的高分辨测向。仿真试验和消声水池试验数据处理结果表明UAVSA-WSSR方法有效地提高了多目标方位分辨性能和测向精度。

1.   矢量阵信号模型及通道噪声功率的不一致性

1.1   矢量阵频域快拍模型

假设$N$个远场窄带信号入射到一个任意$M$元矢量阵, 矢量阵窄带频域第k个频域快拍可表示为

其中, ${{\boldsymbol{y}}_p}(k) = {[{y_{p,1}}(k), \cdots ,{y_{p,M}}(k)]^{\text{T}}}$表示矢量阵声压通道接收数据, 上标T表示转置。${{\boldsymbol{y}}_{vx}}(k) = [{y_{vx,1}}(k), \cdots , {y_{vx,M}}(k)]^{\text{T}}$和${{\boldsymbol{y}}_{vy}}(k) = {[{y_{vy,1}}(k), \cdots ,{y_{vy,M}}(k)]^{\text{T}}}$表示两个振速通道的数据信息。${\boldsymbol{S}}(k) = {[{s_1}(k), \cdots ,{s_N}(k)]^{\text{T}}}$是$N \times 1$维的信号向量。${{\boldsymbol{n}}_p}(k)$, ${{\boldsymbol{n}}_x}(k)$, ${{\boldsymbol{n}}_y}(k)$分别是声压和两个振速通道接收的环境噪声, 本文中假设不同通道之间的噪声互不相关, 且噪声与信号相互独立。${\boldsymbol{\varPhi}} = \{ {\theta _1}, \cdots ,{\theta _N}\}$表示信号的入射方位集合, ${{\boldsymbol{A}}_p}({\boldsymbol{\varPhi}} )$, ${{\boldsymbol{A}}_x}({\boldsymbol{\varPhi}} )$, ${{\boldsymbol{A}}_y}({\boldsymbol{\varPhi}} )$分别为3个通道信号的阵列流形矩阵:

其中, ${{\boldsymbol{a}}_p}({\theta _i})$, ${{\boldsymbol{a}}_x}({\theta _i}) = \cos {\theta _i}{{\boldsymbol{a}}_p}({\theta _i})$和${{\boldsymbol{a}}_y}({\theta _i}) = \sin {\theta _i}{{\boldsymbol{a}}_p}({\theta _i})$表示第i个信号声压和振速通道的导向向量。将矢量阵声压与振速通道输出组成一个向量, 式(1)可简化为

其中, ${\boldsymbol{A}}({\boldsymbol{\varPhi}} ) \,=\, [\, {{\boldsymbol{A}}_p}({\boldsymbol{\varPhi}} ); \; {{\boldsymbol{A}}_x}({\boldsymbol{\varPhi}} ); \; {{\boldsymbol{A}}_y}({\boldsymbol{\varPhi}} ) \,]$, ${\boldsymbol{n}}(k) \,=\, [\,{{\boldsymbol{n}}_p}(k); {{\boldsymbol{n}}_x}(k); \; {{\boldsymbol{n}}_y}(k)]$。

矢量阵接收数据的协方差矩阵可表示为

其中,$\text{E}[ \cdot ]$表示期望, ${{\boldsymbol{R}}_s} = \text{E}[{\boldsymbol{s}}(k){{\boldsymbol{s}}^{\text{H}}}(k)]$表示信号协方差矩阵, 噪声协方差矩阵为

其中,${{\boldsymbol{R}}_v} = {\text{diag}}(\sigma _{np}^2,\sigma _{nx}^2,\sigma _{ny}^2)$表示协方差矩阵, ${\text{diag}}( \cdot )$表示对角矩阵, $\otimes$表示克罗内克积运算。$\sigma _{np}^2$, $\sigma _{nx}^2$, $\sigma _{ny}^2$分别表示声压和两个振速通道的噪声功率。

1.2   矢量阵通道噪声功率的不一致性

由于声压与振速是两个不同的物理量, 一般情况下, 各向同性噪声场中声压通道与振速通道的噪声功率满足^[10]:

在矢量阵噪声模型中, 一般情况下认为声压通道的噪声功率大于振速通道的噪声功率。根据文献[10], 假设${\rho _0}c$是已知的(${\rho _0}$是水的密度, $c$是水中声速), 通过对振速测量值乘以${\rho _0}c$, 则在三维各向同性噪声场中, 振速通道噪声功率是声压通道噪声功率的1/3^[10,22]; 对于二维水平各向同性噪声场, 振速通道噪声功率是声压通道噪声功率的1/2^[23]。文献[22,24]的海试数据也表明, 在大部分频带范围内声压通道噪声功率明显大于振速通道噪声功率。为了便于分析, 下文将统一采用$\sigma _{nv}^2$表示振速通道的噪声功率, 模型假设声压通道噪声功率大于振速通道噪声功率。由于声压通道与振速通道噪声功率差异的存在, 矢量阵的噪声协方差矩阵并不是单位阵。根据式(5)—式(6), 噪声协方差矩阵可以重写为

其中, ${\boldsymbol{h}}_i^{}$是第i个元素为1并且其余元素为0的$3{{M}} \times 1$维向量。将${\boldsymbol{h}}_i^{}$视为信号向量, 那么剩余部分噪声即为均匀噪声。令${\boldsymbol{H}} = [{\boldsymbol{h}}_1^{}, \cdots ,{\boldsymbol{h}}_i^{}, \cdots ,{\boldsymbol{h}}_M^{}]$, 式(4)重新表示为

此时, 接收到的数据的协方差矩阵可改写为

其中, ${\boldsymbol{ C}}$和${{\boldsymbol{R}}_c}$矩阵表示为

根据式(11)可知, ${{\boldsymbol{R}}_c}$是一个$(M + N) \times (M + N)$维的矩阵。当信号不相关时, ${{\boldsymbol{R}}_c}$是一个满秩矩阵, 满足$\text{rank}({{\boldsymbol{R}}_c}) = M + N$。根据文献[11], 由于声压和振速传感器接收到的噪声功率不一致, 产生M个虚拟声源。虚源的存在会导致常规高分辨测向方法的性能下降, 影响矢量阵的多目标方位分辨能力。

2.   加权子空间稀疏恢复测向方法

2.1   声压通道和振速通道的噪声功率

在实际应用中, 只能得到协方差矩阵$\boldsymbol{R}$的估计值:

其中, $K$表示采样的快拍数。结合式(8)—式(12), 矩阵$\widehat {\boldsymbol{R}}$的特征分解表示为

其中, ${\eta _i}$表示第i个特征值, ${{\boldsymbol{u}}_i}$是与之对应的特征向量。特征值由大到小进行排列: ${\eta _1} \geqslant \cdots \geqslant {\eta _{N + M}} \geqslant {\eta _{N + M + 1}} \simeq \cdots \simeq {\eta _{3M}} \simeq \sigma _{nv}^2$。利用声压协方差矩阵的特征值大小关系, 信源的个数N可以通过矢量阵信源数估计方法^[25]求得。在本文中, 假设信源的个数已知。$\sigma _{nv}^2$的估计值可以通过取后面$2M - N$个特征值的平均值获得。

声压通道的协方差矩表示为

其中, 协方差矩阵可以通过${\widehat {\boldsymbol{R}}_p} = \sum\nolimits_{k = 1}^K {{{\boldsymbol{y}}_p}(k){\boldsymbol{y}}_p^{\text{H}}(k)} /K$估计。对协方差矩阵${\widehat{\boldsymbol{ R}}_p}$进行特征分解:

其中, ${\mu _i}$表示第i个特征值, ${{\boldsymbol{v}}_i}$是与其对应的特征向量。特征值由大到小进行排列: ${\mu _1} \geqslant \cdots \geqslant {\mu _N} \geqslant {\mu _{N + 1}} \simeq \cdots \simeq {\mu _M} \simeq \sigma _{np}^2$, $\sigma _{np}^2$的估计值可以通过取后面$M - N$个特征值的平均值获得。在本文中, 假设$N < M$。

2.2   幅度加权

根据2.1节估计的声压和振速通道的噪声功率, 可以构建加权协方差矩阵:

利用${{\boldsymbol{W}}_1}$, 加权后的阵列数据为

其中, $\widetilde {\boldsymbol{A}}({\boldsymbol{\varPhi}} ) = {{\boldsymbol{W}}_1}{\boldsymbol{A}}({\boldsymbol{\varPhi}} )$表示加权后的阵列流形, $\widetilde{\boldsymbol{ n}}(k) = {{\boldsymbol{W}}_1}{\boldsymbol{n}}(k)$是加权后的噪声。显然, 在${{\boldsymbol{W}}_1}$准确估计情况下, $\widetilde {\boldsymbol{n}}(k)$变为功率为$\sigma _{np}^2$均匀噪声。此时, 新的协方差矩阵为

由于加权补偿后的矢量阵声压通道和振速通道的噪声变为均匀噪声, ${\boldsymbol{R}}'$中的噪声协方差矩阵是单位阵乘以P通道的噪声功率。${\boldsymbol{R}}'$的特征值分解为

其中, ${\alpha _i}$表示第i个特征值, ${\boldsymbol{b}_i}$是与其对应的特征向量。特征值由大到小进行排列${\alpha _1} \geqslant \cdots \geqslant {\alpha _N} \geqslant {\alpha _{N + 1}} \simeq \cdots \simeq {\alpha _{3M}} \simeq \sigma _{np}^2$。${{\boldsymbol{U}}_s} = [{{\boldsymbol{b}}_1}, \cdots ,{{\boldsymbol{b}}_N}]$和${{\boldsymbol{U}}_n} = [{{\boldsymbol{b}}_{N + 1}}, \cdots ,{{\boldsymbol{b}}_{3M}}]$分别表示信号子空间和噪声子空间。${{\boldsymbol{\varLambda}} _s}$和${{\boldsymbol{\varLambda}} _n}$分别表示包含信号和噪声特征值的对角矩阵。

2.3   基于信号子空间加权的稀疏恢复

现有文献[26, 27]揭示了大部分测向方法都可以转化为一个信号子空间拟合问题, 表示为

其中, $\left\| \cdot \right\|_{\text{F}}^{}$表示Frobenius范数。$\overline {\boldsymbol{\varPhi}} = \{ {\overline \theta _1}, \cdots ,{\overline \theta _{N'}}\}$是空间域中所有潜在方向的采样网格。一般情况下, $N \ll N'$。注意$\widetilde {\boldsymbol{A}}(\overline {\boldsymbol{\varPhi}} )$是已知的, 不依赖于实际信号方向。这意味着只要找到非零值的位置, 就可以估计信号源的方向, 即将稀疏恢复问题转化为DOA估计问题。${{\boldsymbol{W}}_2}$是一个加权矩阵, 根据文献[26], 最小渐近方差的最优加权为${{\boldsymbol{W}}_2} = {({{\boldsymbol{\varLambda}} _s} - \sigma _{np}^2{{\boldsymbol{I}}_N})^2}{\boldsymbol{\varLambda}} _s^{ - 1}$。${\boldsymbol{Q}}$是一个$N' \times N$的满秩矩阵。矩阵${\boldsymbol{Q}}$的行元素与信号幅值相对应, 在非目标方位处的值是一个近似为0的小量。这意味着, 可以利用式(20)求解矩阵${\boldsymbol{Q}}$, 然后利用${\boldsymbol{Q}}$在空域水平方位平面的分布情况实现目标方位的估计。

式(20)的稀疏问题可以表示为如下联合稀疏优化形式:

其中, ${{\boldsymbol Q}^{{l_{2,0}}}} = {[{\boldsymbol Q}_1^{({l_2})},{\boldsymbol Q}_2^{({l_2})}, \cdots ,{\boldsymbol Q}_{N'}^{({l_2})}]^{({l_0})}}$。${[ \cdot ]^{({l_0})}}$表示向量的L₀范数, 并且${\boldsymbol Q}_i^{({l_2})}$ 表示矩阵${\boldsymbol{Q}}$的第i行元素的L₂范数。$\beta$是信号子空间拟合的误差上界, 它表示稀疏子空间的拟合误差。假设矩阵${\boldsymbol{Q}}$被准确恢复, 拟合误差为$\sqrt {{\boldsymbol{Q}}(\overline {\boldsymbol{\varPhi}} )}$。根据文献[26]的结论, 如果$\overline {\boldsymbol{\varPhi}}$等于正确的方位, 则函数$(2K/\sigma _{np}^2){\boldsymbol{Q}}(\overline {\boldsymbol{\varPhi}} )$在自由度为$2N \cdot (3M - N)$的卡方分布下渐近成立。因而, $\sqrt {Q(\overline {\boldsymbol{\varPhi}} )} \leqslant \beta$在概率为$q$下成立。在本文中, 设定$q = 0.995$。

注意, 式(21)是非凸的, 难以求解。一种替代方法是使用L₁范数代替L₀范数约束稀疏性^[17]:

其中, ${{\boldsymbol Q}^{{l_{2,1}}}} = {[{\boldsymbol{Q}}_1^{({l_2})},{\boldsymbol{Q}}_2^{({l_2})}, \cdots ,{\boldsymbol{Q}}_{N'}^{({l_2})}]^{({l_1})}}$, ${[ \cdot ]^{({l_1})}}$表示向量的L₁范数。然而, L₁范数约束并没有考虑信号幅度的影响, 并不能保证${\boldsymbol Q}_i^{({l_2})} \leqslant 1$, 这意味着L₁范数并不是L₀范数的凸包。因此, 式(22)并不是式(21)良好的凸优化逼近。

2.4   重加权L₁范数最小化设计

当${\boldsymbol{Q }}$中的${\boldsymbol Q}_i^{({l_2})}$是0或者1时, ${{\boldsymbol Q}^{{l_{2,0}}}} = {{\boldsymbol Q}^{{l_{2,1}}}}$。然而, 在大部分情况下, ${\boldsymbol Q}_i^{({l_2})}$是一个0和1以外的非负值。文献[21]表明, 通过加权可以增强真实信号位置处的稀疏性, 大的权值可以阻止恢复信号中的非零项, 而小的权值可以鼓励非零项。根据上述原则, 加权系数需要满足如下特性:

其中, ${w_i}$表示${\boldsymbol{Q}}$中${\boldsymbol Q}_i^{({l_2})}$的加权系数。根据式(23), 当${\boldsymbol Q}_i^{({l_2})} \ne 0$时, ${w_i}{\boldsymbol Q}_i^{({l_2})} = 1$; 当${\boldsymbol Q}_i^{({l_2})} = 0$时, ${w_i}{\boldsymbol Q}_i^{({l_2})} = 0$。因而, 式(23)的加权系数会增强信号稀疏性, 使得${{\boldsymbol Q}^{{l_{2,1}}}}$逼近${{\boldsymbol Q}^{{l_{2,0}}}}$。当然, ${\boldsymbol{Q}}$是一个未知量, 直接构造式(23)的权值是很困难的。

根据式(15), 信号子空间和噪声子空间对应的特征值都是大于0的正数。在没有信号情况下, ${\boldsymbol Q}_i^{({l_2})}$对应各向同性噪声幅度并且是一个大于0的小值, 与式(23)中${\boldsymbol Q}_i^{({l_2})} = 0$情况近似, ${w_i}$是一个较大的权值; 在有信号情况下, ${\boldsymbol Q}_i^{({l_2})} \ne 0$并且${w_i}$是一个相对较小的权值。共轭转置后的搜索导向矢量与噪声子空间乘积的范数恰好符合上述规律^[20], 用其来近似代替式(23)中的权值:

其中, ${\overline \theta _i} \in \overline {\boldsymbol{\varPhi}}$。显然, $\widetilde {\boldsymbol{A}}(\overline {\boldsymbol{\varPhi}} )$中真实信号对应的共轭转置后导向矢量与噪声子空间乘积的L₂范数数值较小; 而$\widetilde {\boldsymbol{A}}(\overline {\boldsymbol{\varPhi}} )$中其余共轭转置后导向矢量与噪声子空间乘积的L₂范数数值较大, 详细推导参考附录A。

根据附录A,一般情况下, 真实信号以外的导向矢量对应的权值远大于真实信号的导向矢量对应的权值。权值特性符合文献[21]对加权系数的要求, 可以起到一定增强信号稀疏性的作用。令${{\boldsymbol{W}}_3} = \text{diag} {\text{([}}{\overline w_1}{\text{,}}{\overline w_2}, \cdots ,{\overline w_{N'}}{\text{])}}$, 式(21)可以转化为重加权L₁范数最小问题:

式(25)的优化问题可以在二阶锥优化框架^[17]下进行求解, 它的二阶锥优化规范形式如下:

其中, $\boldsymbol{1}$是$3M \times 1$维的1向量。对于凸优化问题的数值求解, 式(26)可通过CVX工具箱等编程软件包计算得到。

2.5   算法流程

根据上述内容, UAVSA-WSSR算法的流程图如图1所示, 具体实施步骤如下:

图  1  流程框图

下载: 全尺寸图片 幻灯片

(1) 根据矢量阵接收数据, 分别估计声压通道和振速通道的噪声功率: $\widehat \sigma _{np}^2$和$\widehat \sigma _{nv}^2$, 进而得到加权矩阵${{\boldsymbol{W}}_1}$;

(2) 利用${{\boldsymbol{W}}_1}$对接收数据进行加权, 获得新的协方差矩阵${\boldsymbol{R}}'$。通过对协方差矩阵${\boldsymbol{R}}'$进行特征分解, 得到信号子空间${{\boldsymbol{U}}_s}$、噪声子空间${{\boldsymbol{U}}_n}$和特征值;

(3) 根据步骤(2)中的特征值, 计算信号子空间的加权矩阵${{\boldsymbol{W}}_2}$;

(4) 利用噪声子空间和过完备基矩阵之间的关系, 设计L₁范数最小化的权值${{\boldsymbol{W}}_3}$;

(5) 根据搜索的角度范围设定$\widetilde {\boldsymbol{A}}(\overline {\boldsymbol{\varPhi}} )$, 并且设定加权信号子空间稀疏重构的误差上界$\beta$, 采用CVX工具箱求解式(26), 得出空间谱。

3.   仿真结果及分析

通过仿真分析验证UAVSA-WSSR方法的性能, 并且与MUSIC^[28]、AS-MUSIC^[11]和 L₁-SVD^[17]比较, 其中L₁-SVD的正则化系数选用2, UAVSA-WSSR的概率系数$q$选用0.995。主要研究以下性能: 方位估计的有效性, 方位估计的分辨概率和均方根误差随信噪比变化的情况, 方位估计的分辨概率和均方根误差随角度间隔变化的情况。考虑一个四元均匀线列矢量阵, 阵元间距0.5 m, 其中每个矢量传感器由声压分量和二维质点振速分量组成, 两个振速分量的方向分别沿x轴和y轴, 矢量传感器沿x轴布放。假设两个远场信号频率为1.5 kHz, 信号源互不相关并且相互独立, 背景噪声是各向同性高斯噪声, 频域快拍数为5000。此部分直接仿真频域快拍模型^[29], 采用零均值高斯随机过程的采样函数模拟信号在1.5 kHz下的5000个频域快拍, 同样采用零均值高斯随机过程的另一个采样函数模拟噪声在1.5 kHz下的5000个频域快拍。振速通道和声压通道的窄带噪声功率采用常用的理想矢量阵各向同性噪声模型, 分别设定为$\sigma _{nx}^2 = \sigma _{ny}^2 = 0.5$和$\sigma _{np}^2 = 1$。需要注意的是, 在实际工程应用中, 不一定会对振速测量值乘以${\rho _0}c$, 因此, 振速通道噪声功率也可能是其他值。第i个信源的信噪比定义为$10{\lg _{10}}(\sigma _{{s_i}}^2/\sigma _{np}^2)$, $i = 1,2$, 其中, $\sigma _{{s_i}}^2$和$\sigma _{np}^2$分别为声压通道第i个信源和噪声的窄带功率。

首先, 分析信噪比为10 dB, 5 dB, 0 dB, −5 dB, −10 dB时几种方法的方位估计结果。假设两个信号的方位分别为90°和100°, 信号功率相等。空域角度间隔为1°, 归一化空间谱对比结果如图2所示。在信噪比10 dB情况下, MUSIC方法、AS-MUSIC方法、L₁-SVD方法和UAVSA-WSSR方法均可以分辨两个信源, 但L₁-SVD方法的估计误差明显大于其他3种方法。在信噪比5 dB情况下, L₁-SVD方法无法分辨两个信源, 而MUSIC方法、AS-MUSIC方法和UAVSA-WSSR方法可以分辨出两个信源。在信噪比为−5 dB和0 dB情况下, MUSIC方法和 L₁-SVD方法无法分辨出两个信源, 而AS-MUSIC方法和UAVSA-WSSR方法可以分辨出两个信源; 与AS-MUSIC方法相比, UAVSA-WSSR方法有着更尖锐的谱峰。在信噪比−10 dB情况下, MUSIC方法、AS-MUSIC方法和 L₁-SVD方法完全无法分辨两个信源, 而UAVSA-WSSR方法仍然可以分辨出两个信源, 并且有效估计两个目标的方位。图2的结果表明, 与其他3种方法相比, UAVSA-WSSR方法的多目标分辨能力更优, 并且对旁瓣抑制能力更强。

图  2  不同信噪比下归一化空间谱 (a) SNR = 10 dB; (b) SNR = 5 dB; (c) SNR = 0 dB; (d) SNR = −5 dB; (e) SNR = −10 dB

下载: 全尺寸图片 幻灯片

接着, 分析几种方法DOA估计的分辨概率随信噪比变化的情况。为了获得准确的统计结果, 空域角度间隔设为0.1°, 每个信噪比条件下均进行200次独立的蒙特卡罗实验。在每次实验中, 若$| {{{\widehat \theta }_1} - {\theta _1}}|$和$| {{{\widehat \theta }_2} - {\theta _2}} |$都小于${{\left| {{\theta _1} - {\theta _2}} \right|} \mathord{\left/ {\vphantom {{\left| {{\theta _1} - {\theta _2}} \right|} 2}} \right. } 2}$, 则定义该次实验中两个方位被成功分辨; 否则, 定义其分辨失败。假设两个等功率的远场信号从${90^ \circ }$和${100^ \circ }$方向入射, 设置信噪比从−16 dB变化到18 dB, DOA估计的分辨概率随信噪比的变化如图3所示。图3的仿真结果表明, UAVSA-WSSR方法在−10 dB的分辨概率就稳定达到1, 远高于同等条件下其他3种方法。AS-MUSIC在−4 dB达到1, MUSIC和 L₁-SVD分别在2 dB和10 dB达到1。显然, 在信噪比为−6 dB到0 dB, MUSIC和 L₁-SVD的分辨概率远差于AS MUSIC和UAVSA-WSSR方法, 这主要是因为MUSIC和 L₁-SVD受矢量阵声压和振速通道噪声不一致的影响较大。

图  3  DOA估计的分辨概率随信噪比的变化

下载: 全尺寸图片 幻灯片

随后, 考察几种方法在不同信噪比下对方位近邻目标进行估计的精度, 实验条件和参数与图3一致。DOA的估计精度由均方根误差来衡量, 其中均方根误差定义为

考虑到分辨概率太低情况下, 只能估计到一个目标方位, 对统计两个方位均方根误差的意义不大。根据图2和图3结果, L₁-SVD方法的多目标分辨性能远差于其他3种方法。因而, 此部分主要统计MUSIC方法、AS-MUSIC方法和UAVSA-WSSR方法在完全能分辨出两个目标方位情况下(即信噪比从2 dB到18 dB)的均方根误差, 结果如图4所示。随着信噪比的增大, 3种方法的均方根误差都在降低。在信噪比较低时, UAVSA-WSSR方法的均方根误差最小, AS-MUSIC方法次之, MUSIC方法误差最大且性能差距明显。在信噪比较高时, 3种方法的均方根误差的差距减小。因而, UAVSA-WSSR方法的估计精度优于另两种方法。

图  4  DOA估计的均方根误差随信噪比的变化

下载: 全尺寸图片 幻灯片

接下来, 考察几种方法在不同角度间隔近邻目标下的分辨能力。假设两个信噪比为0 dB的远场信号入射到线列矢量阵, 两个信号的角度间隔从2°变化到32°, 其他条件与图3一致, 图5为DOA估计的分辨概率随两个信源角度间隔的变化情况。两目标方位间隔为4°时, UAVSA-WSSR方法有0.995的概率将其成功分辨, 而其他3种方法完全不能将其分辨。这说明UAVSA-WSSR方法在小角度间隔情况下的多目标方位分辨性能优于其他方法。随着角度间隔的增大, 其他3种方法的分辨率都在不断提高, 其中AS-MUSIC和MUSIC方法比L₁-SVD方法的分辨率先达到1。这说明矢量阵声压和振速通道的差异对L₁-SVD方法的分辨率影响很大。

图  5  DOA估计的分辨概率随角度间隔的变化

下载: 全尺寸图片 幻灯片

然后, 考察几种方法在不同角度间隔近邻目标下的估计精度, 实验条件和参数与图5一致。与图4类似, 考虑到图5中L₁-SVD方法的多目标分辨性能远差于其他方法, 此部分主要统计MUSIC方法、AS-MUSIC方法和UAVSA-WSSR方法在完全能分辨出两个目标方位情况下(即两个信源角度间隔从12°到32°)的均方根误差, 结果如图6所示。随着两个信源角度间隔的增大, 3种方法的均方根误差都在降低。综上所述, 相比于对比方法, UAVSA-WSSR方法在小角度间隔条件下具有更优的估计精度。

图  6  DOA估计的均方根误差随角度间隔的变化

下载: 全尺寸图片 幻灯片

最后, 考察几种方法对一强一弱目标的分辨能力。两个信号的方位分别为90°和100°, 信噪比分别为0 dB和20 dB, 实验条件的其他参数与图2一致。图7为4种方法的方位历程图, 图8是图7方位历程中截取第1秒单帧的结果对比。UAVSA-WSSR方法和AS MUSIC方法可以分辨出两个目标, 而MUSIC 方法和L₁-SVD方法只能估计出强目标方位。虽然AS MUSIC方法能分辨两个目标方位, 但是主瓣宽度要大于UAVSA-WSSR方法。另外, UAVSA-WSSR方法对噪声抑制能力强于其他3种方法。

图  7  方位历程图 (a) MUSIC; (b) AS MUSIC; (c) L₁-SVD; (d) UAVSA-WSSR

下载: 全尺寸图片 幻灯片

图  8  第1秒单帧的空间谱图

下载: 全尺寸图片 幻灯片

4.   水池试验结果

本节利用消声水池实测数据验证UAVSA-WSSR方法的性能。本文选用3个矢量水听器组成的直线阵, 阵元间距为0.1875 m。此外, 单个矢量水听器由1个声压通道和2个振速通道组成。两个目标声源和矢量基阵入水深度相同, 均为4 m。两个声源距离矢量基阵均为11 m左右, 分别发射带宽为2.8~4.1 kHz的连续噪声, 两个信号互相独立并且不相关。文献[22,30,31]的研究表明, 两个传感器间距为半波长时, 其声压通道噪声之间的相关系数为0, 而振速通道噪声之间的相关系数也是一个不大的值。实验采用直线矢量阵的半波长频率是4 kHz, 相邻两个矢量水听器的声压通道噪声在3980~4020 Hz的相关系数为0.0064, Vx振速通道噪声在3980~4020 Hz的相关系数为0.0061, Vy振速通道噪声在3980~4020 Hz的相关系数为0.0060。相邻矢量水听器的3个通道相关系数很小而且数值相近, 此窄带频率下的不同水听器的噪声可近似为互相独立的噪声。因而, 选用3980~4020 Hz的窄带频率验证方法性能。通过快速傅里叶变换(FFT), 选取4 kHz处带宽为40 Hz的窄带信号作为阵列输出矢量, 频域快拍数为41, 分别采用MUSIC、AS-MUSIC、L₁-SVD和UAVSA-WSSR方法进行处理。

图9是1号矢量水听器接收的信号波形结果, 图10是1号矢量水听器接收的噪声结果。声压通道的噪声时域波形幅度大于振速通道, 而且声压通道噪声的频谱明显也高于振速通道噪声的频谱。表1为4 kHz下三个矢量水听器的声压P通道与两个振速V_x和V_y通道的实测噪声关系, 噪声功率选取3980~4020 Hz的频率进行计算, 其中$\sigma _{n\overline p}^2$是声压通道噪声功率参考值(选取1号矢量水听器)。4 kHz下振速通道噪声幅值大约是声压通道噪声幅值的0.5倍。而UAVSA-WSSR方法估计的$\sqrt { {\sigma _{ny}^2} / {\sigma _{np}^2} }$为0.49, 与实际的通道噪声关系近似一致。值得注意的是, 振速通道噪声幅值大约是声压通道噪声幅值的0.5倍, 那么振速通道的噪声功率是声压通道噪声功率的0.25倍, 与理想矢量阵各向同性噪声模型不一致。这可能是未测量${\rho _0}c$进行振速归一化和矢量水听器灵敏度差异造成的。表2为4 kHz下两个信号源在3个矢量水听器声压通道的实测信噪比, 信噪比通过计算信号和噪声在频率为3980~4020 Hz下功率之比所得。对于同一个信号源, 不同矢量水听器声压通道的信噪比存在一些差异, 推测是由水听器灵敏度和电路噪声等差异引起; 两个信号源的信噪比大概差20 dB。由于本实验采用的搜索导向向量是水池校准后的导向向量, 矢量水听器信噪比的差异对方法性能的验证影响较小。

图  9  1号矢量水听器接收信号波形结果 (a) 时域; (b) 频域

下载: 全尺寸图片 幻灯片

图  10  1号矢量水听器接收噪声结果 (a) 时域; (b) 频域

下载: 全尺寸图片 幻灯片

表  1  P通道噪声与V_x通道和V_y通道在4 kHz处噪声关系

水听器序号	$\sqrt {{{\sigma _{nx}^2} \mathord{\left/ {\vphantom {{\sigma _{nx}^2} {\sigma _{np}^2}}} \right. } {\sigma _{np}^2}}}$	$\sqrt {{{\sigma _{ny}^2} \mathord{\left/ {\vphantom {{\sigma _{ny}^2} {\sigma _{np}^2}}} \right. } {\sigma _{np}^2}}}$	$\sqrt {{{\sigma _{np}^2} \mathord{\left/ {\vphantom {{\sigma _{np}^2} {\sigma _{n\overline p}^2}}} \right. } {\sigma _{n\overline p}^2}}}$
1	0.49	0.48	1
2	0.52	0.52	0.95
3	0.54	0.52	0.94

下载: 导出CSV 

| 显示表格

表  2  两个信号源在4 kHz处声压通道的信噪比

水听器序号	信号1的SNR (dB)	信号2的SNR (dB)
1	16.44	36.69
2	8.35	28.39
3	20.13	40.19

下载: 导出CSV 

| 显示表格

图11为4种方法在4 s数据下的方位历程图, 从上到下分别是MUSIC、AS-MUSIC、L₁-SVD和UAVSA-WSSR的结果。图12是4种方法对第1 秒单帧数据估计的归一化空间谱图, 其中红色圆圈表示两个目标的方位, 分别位于74°和86°。考虑到MUSIC、AS-MUSIC的空间谱与L₁-SVD、UAVSA-WSSR的空间谱旁瓣抑制能力差异很大, 为便于观察分两子图给出, 图12(a)为MUSIC和AS-MUSIC的归一化空间谱估计结果, 图12(b)为L₁-SVD和UAVSA-WSSR的归一化空间谱估计结果。由图11可见, MUSIC和AS-MUSIC的谱峰更宽, 而L₁-SVD和UAVSA-WSSR的谱峰更尖锐。这主要是由方法特性决定的, L₁-SVD和UAVSA-WSSR属于稀疏恢复类方法, MUSIC和AS-MUSIC属于子空间类方法, 而稀疏恢复类方法的谱峰尖锐程度本就优于子空间类方法^[17]。由图12可见, MUSIC、AS-MUSIC和L₁-SVD只能估计出一个目标, 而UAVSA-WSSR可以较为准确估计出两个目标。

图  11  方位历程图(上至下: MUSIC, AS-MUSIC, L₁-SVD, UAVSA-WSSR)

下载: 全尺寸图片 幻灯片

图  12  第1秒单帧数据的空间谱图 (a) MUSIC和AS-MUSIC; (b) L₁-SVD和UAVSA-WSSR

下载: 全尺寸图片 幻灯片

5.   结论

本文提出一种应用加权子空间稀疏恢复的水声矢量阵测向方法。该方法通过噪声归一化, 克服了矢量阵声压通道与振速通道噪声功率不一致对高分辨方法性能的负面影响; 利用加权子空间拟合和重加权L₁范数最小化技术, 将方位估计问题转化为加权子空间稀疏恢复问题, 提高了矢量阵的多目标分辨能力。仿真结果和消声水池试验结果表明, 所提方法的空间谱在旁瓣抑制、多目标方位分辨及测向精度方面明显优于其他对比方法, 尤其是在信噪比较低和小角度间隔情况下, 本文方法可以有效地提高水声矢量阵的高分辨测向能力。

附录A

文献[28]证明了真实信号对应导向矢量与噪声子空间满足:

其中, 有足够的快拍数时, ${\left( {{{\boldsymbol{E}}_1}} \right)_{i,j}} \approx 0$。${\left( \cdot \right)_{i,j}}$表示矩阵的第i行和第j列元素。$\widetilde {\boldsymbol{A}}(\overline {\boldsymbol{\varPhi}} )$与噪声子空间之间关系可以表示为

其中, ${{\boldsymbol{E}}_2} = \widetilde {\boldsymbol{A}}_{c,{\boldsymbol{\varPhi}} }^{\text{H}}(\overline {\boldsymbol{\varPhi}} ){{\boldsymbol{U}}_n}$。$\widetilde {\boldsymbol{A}}_{c,{\boldsymbol{\varPhi}} }^{\text{H}}(\overline {\boldsymbol{\varPhi}} )$表示$\overline {\boldsymbol{\varPhi}}$中除${\boldsymbol{\varPhi}}$之外其他方位角对应的导向向量。当有足够的快拍数时, ${\boldsymbol{E}}$第i行的L₂范数满足

其中, ${0^ + }$表示大于0但是非常接近0。${c_i}$是一个大于0的值, 一般情况下${c_i} \gg {0^ + }$。这与式(23)的特性是一致, 即信号部分被赋予较小的权值, 没有信号部分赋予较大权值。在实际应用中使用有限快拍时, 式(A3)依然适用。结合式(24), $\left( \boldsymbol{E} \right)_{i,:}^{({l_2})}$与${\overline w_i}$等同。