CN103605108B - 声矢量阵高精度远程方位估计方法 - Google Patents

Abstract

本发明提供的是一种声矢量阵高精度远程方位估计方法。(1)对二维矢量水听器阵列的接收信号进行窄带滤波，获得待处理频点上的窄带输出信号；(2)将二维矢量水听器阵列相互正交的两个振速分量在复数域进行线性组合，转换为两个新的振速输出分量；(3)采用最大似然估计求出一定数量采样快拍下的声压及复数域双振速分量的互协方差矩阵对；(4)应用矩阵束的ESPRIT算法计算基于广义声能流的声压振速互协方差矩阵对之间的旋转不变因子，从而进行方位估计。本发明能够在任意阵型甚至未知阵型的情况下进行高精度的方位估计，不会出现某些方向的信号被严重削弱甚至完全屏蔽的现象。此外具有更低的可处理信噪比门限。

Description

声矢量阵高精度远程方位估计方法

技术领域

本发明涉及的是一种水声探测及水声定位方法，具体涉及一种声矢量阵的方位估计方法。

背景技术

声场兼具有标量场和矢量场，其中声压提供了标量信息，振速提供了矢量信息，该矢量方向与声波的传播方向一致。矢量水听器由声压水听器和振速水听器组合而成，能够同时拾取声场中的声压与振速信息，相较于传统的声压水听器阵列，矢量水听器阵列获得了更多的信息，因而能够提供更大的处理空间。

国内外学者对基于声矢量阵的方位估计应用进行了大量的研究。文献[1],MalcolmHawkes和Arye Nehorai(Acoustic Vector-sensor Corrections in Ambient Noise.IEEE J.Oceanic Eng,2001,26(3):337-347.)对各向同性噪声场中的声压与振速的相关特性进行了探讨，并给出了声压与振速相关系数的计算公式。文献[2]白兴宇等(一种新的声矢量阵远程ESPRIT方位估计算法.哈尔滨工程大学学报,2006,27(6):891～895.）、文献[3]（基于声压振速联合处理的声矢量阵信源数检测与方位估计.声学学报,2008,33(1):56～61.)提出将声矢量阵的声压振速互协方差矩阵应用于方位估计中，将声矢量阵的抗各向同性噪声能力与高分辨方位估计算法结合在一起，取得了较好的效果。但该方法的振速分量具有空间滤波效果，可能导致感兴趣的信号被屏蔽掉。文献[4]张揽月等(矢量阵的非空间ESPRIT算法.哈尔滨工程大学学报,2009,30(4):406～410.)利用声矢量阵的声压与振速传感器阵列之间的旋转不变性，能够在未知阵列流形的情况下对目标进行方位估计，但该方法未利用声压与振速通道噪声的不相关特性，估计精度还有待提高。

发明内容

本发明的目的在于提出一种能够用于任意阵型甚至未知阵型声矢量阵列，并且结合了声矢量阵抗各向同性噪声能力，精度高的声矢量阵高精度远程方位估计方法。

本发明的目的是这样实现的：

(1)对二维矢量水听器阵列的接收信号进行窄带滤波，获得待处理频点上的窄带输出信号；

(2)将二维矢量水听器阵列相互正交的两个振速分量在复数域进行线性组合，转换为两个新的振速输出分量，即构造复数域双振速分量；

(3)采用最大似然估计求出一定数量采样快拍下的声压及复数域双振速分量的互协方差矩阵对，即构造出基于广义声能流的声压振速互协方差矩阵对；

(4)应用矩阵束的ESPRIT算法计算基于广义声能流的声压振速互协方差矩阵对之间的旋转不变因子，从而进行方位估计。

本发明还可以包括：

1.所述的构造复数域双振速分量具体包括：将y方向的振速输出分量v_y与复数单位j相乘，再与x方向的振速输出分量v_x相加获得新的振速输出分量v₁；将y方向的振速输出分量v_y与负的复数单位-j相乘，再与x方向的振速输出分量v_x相加获得新的振速输出分量v₂。v₁和v₂构成了复数域双振速分量，它们依然携带了振速场的方向信息，并且具有全向的接收能力，是对常规振速分量的一种改进。

2.所述的构造出基于广义声能流的声压振速互协方差矩阵对具体包括：将一定采样快拍数量的声压传感器阵列的输出分别与两个新的复数域振速输出分量的共轭转置相乘，并求其平均值，得到基于广义声能流的互协方差矩阵对R₁和R₂。

3.所述的基于广义声能流的ESPRIT算法，其思想是：求解上述所得的声压振速协方差矩阵对的非零广义特征值，从而给出目标的方位估计结果。为了避免在求解非零广义特征值的过程中出现潜在的数值困难，所述应用矩阵束的ESPRIT算法计算基于广义声能流的声压振速互协方差矩阵对之间的旋转不变因子中采用奇异值分解(SVD)和总体最小二乘法(TLS)将一个较大维数的病态最小二乘问题转化为一个较小维数的非病态总体最小二乘问题。

本发明的优点体主要体现在：(1)常规声矢量阵的ESPRIT算法采用数据的自相关矩阵进行特征分解，即使在理想情况下，其信号子空间中仍然存在噪声。对于声矢量阵，当阵元间距满足声压与振速的噪声相关半径时，声压与振速分量之间的噪声是不相关的，本发明采用声压与振速的互协方差矩阵，充分利用了声压与振速噪声场的不相关特性，具有良好的抗各向同性噪声能力；(2)常规的振速输出分量以及基于Givess旋转的组合振速输出都具有空间指向性，对不同方向的入射信号具有不同的幅值加权，尤其是从观测方向的正交方向入射的信号会被完全屏蔽掉。本发明对振速输出分量做了改进，对各方向的入射信号都具有相同的幅度响应，能够完整接收到从各方向入射的信号；(3)本发明所构造的振速输出分量之间具有旋转不变性，其对应的旋转因子仅与信号的入射方向有关，而与阵列的空间位置无关，应用ESPRIT算法能够实现在任意阵型甚至未知阵型的情况下进行高精度的方位估计。

附图说明

图1为各向同性噪声场中声压振速的空间相关函数曲线；

图2为给出的一组不规则阵型声矢量阵位置示意图；

图3为本发明的实现框图；

图4(a)-图4(e)为单目标信号从30°方向入射时声压、振速以及复数域双振速通道输出频谱，其中：图4(a)声压通道频谱，图4(b)振速通道v_x频谱，图4(c)振速通道v_y频谱，图4(d)振速通道v₁频谱，图4(e)振速通道v₂频谱；

图5(a)-图5(e)为双目标信号分别从30°、40°方向入射时声压、振速以及复数域双振速通道输出频谱，其中：图5(a)声压通道频谱，图5(b)振速通道v_x频谱，图5(c)振速通道v_y频谱，图5(d)振速通道v₁频谱，图5(e)振速通道v₂频谱；

图6为单目标方位估计均方根误差曲线；

图7为双目标方位估计均方根误差曲线；

图8(a)-图8(e)为信号从y轴方向入射时声压、振速以及复数域双振速通道输出频谱，其中：图8(a)声压通道频谱，图8(b)振速通道v_x频谱，图8(c)振速通道v_y频谱，图8(d)振速通道v₁频谱，图8(e)振速通道v₂频谱；

图9(a)-图9(b)为从y轴附近方向入射的信号的方位估计结果，图9(a)是采用常规振速输出分量的方位估计结果，图9(b)是采用复数域双振速分量的方位估计结果。

具体实施方式

下面结合附图和具体实施例对本发明作进一步具体说明：

1．计算步骤

1.1.阵列输出模型

假设N元二维声矢量阵的第i个阵元位置的极坐标为[r_i,α_i],i=1,2,…,N，D个中心频率为ω的窄带远场平面波入射到基阵上，其中第k个信号的入射方向为θ_k，则阵列的输出信号可记为：

$\left\{\begin{array}{c}p\left(t\right)=A_{p}s\left(t\right)+n_p\left(t\right)\\ v_x\left(t\right)=A_{x}s\left(t\right)+n_x\left(t\right)\\ v_y\left(t\right)=A_{y}s\left(t\right)+n_y\left(t\right)\end{array}\right.---\left(1\right)$

式中p(t)、v_x(t)、v_y(t)分别是声压通道和x、y振速通道的N×1维接收快拍向量；s(t)是D个远场信号源在参考位置的D×1维声压向量；n_p(t)、n_x(t)、n_y(t)是各通道对应的N×1维本地噪声向量。

A_p=[a_p1,…,a_pD]是信号向量对应的声压通道流形矩阵，第k个信号的流形矢量可写为：

a_pk=[exp(jωr_Ncos(θ_k-α₁)/c),…,exp(jωr_Ncos(θ_k-α_N)/c)]^T   (2)

运算符(·)^T表示矩阵转置，A_x=[a_x1,…,a_xD]是信号向量对应的x方向振速通道流形矩阵，第k个信号的流形矢量可写为：

a_xk=[exp(jωr_Ncos(θ_k-α₁)/c),…,exp(jωr_Ncos(θ_k-α_N)/c)]^Tcosθ_k   (3)

A_y=[a_y1,…,a_yD]是信号向量对应的y方向振速通道流形矩阵，第k个信号的流形矢量可写为：

a_yk=[exp(jωr_Ncos(θ_k-α₁)/c),…,exp(jωr_Ncos(θ_k-α_N)/c)]^Tsinθ_k   (4)

由式(2)-(4)可知，A_p、A_x、A_y之间存在如下关系：

A_x=A_pΦ_x,A_y=A_pΦ_y,A_y=A_xΦ_xy   (5)

式中，Φ_x、Φ_y、Φ_xy均是对角矩阵，并且有：

$\left\{\begin{array}{c}{\mathrm{\Φ}}_x=\mathrm{diag}[{\cos \mathrm{\θ}}_1,\·\·\·,{\cos \mathrm{\θ}}_D]\\ {\mathrm{\Φ}}_y=\mathrm{diag}[{\sin \mathrm{\θ}}_1,\·\·\·,{\sin \mathrm{\θ}}_D]\\ {\mathrm{\Φ}}_{\mathrm{xy}}=\mathrm{diag}[{\tan \mathrm{\θ}}_1,\·\·\·,{\tan \mathrm{\θ}}_D]\end{array}\right.---\left(6\right)$

从式(5)可以看出，对于远场平面波入射信号，声矢量阵的声压传感器阵列以及各振速分量传感器阵列之间都存在旋转不变关系，并且这种旋转因子仅与信号入射方向有关，而与阵列的空间位置无关，这是声矢量阵能够在任意阵型甚至未知阵型的情况下应用ESPRIT算法进行方位估计的物理基础。

1.2.新的复数域双振速输出分量

在各向同性噪声场中，声压振速的空间相关函数曲线如图1所示。曲线①显示了声压的空间自相关函数曲线，在d=0(d是空间距离)处相关性最强，而在d=0.5λ(λ是信号的波长)处声压的自相关函数值为0，因此基于声压阵的算法通常采用半波长间距。曲线②显示了声压与振速的空间互相关函数，从图中来看，空间共点的声压与振速是不相关的，而在d=0.715λ处，其相关函数值也为0，可知声压与振速的噪声相关半径约为0.715λ。与传统声压阵不同的是，若考虑声矢量阵的阵元间距满足声压与振速的噪声相关半径时，此时有：

$E\left(n_p{n}_x^H\right)=E\left(n_p{n}_y^H\right)=0---\left(7\right)$

式中，运算符(·)^H表示矩阵的共轭转置。由式(7)可知，采用声压与振速的互协方差矩阵，将具有良好的抗各向同性噪声的能力。值得注意的是，由于振速传感器自身的偶极子指向性，常规的声压振速互协方差矩阵具有空间滤波的效果，有可能会导致感兴趣方向的信号被削弱甚至被完全屏蔽掉。以声压与x方向的振速分量的互协方差矩阵为例，从y轴方向入射的信号将会被完全屏蔽掉。为了解决这个问题，本发明提出在复数域对振速分量进行线性组合，构造出新的复数域双振速输出分量v₁、v₂：

$\left\{\begin{array}{c}{v}_1=v_x+{\mathrm{jv}}_y=A_p{\mathrm{\Φ}}_{1}s+n_1\\ v_2=v_x-{\mathrm{jv}}_y=A_p{\mathrm{\Φ}}_{2}s+n_2\end{array}\right.---\left(8\right)$

式中，n₁、n₂是对应的背景噪声向量；Φ_x、Φ_y均是对角矩阵，并且有：

$\left\{\begin{array}{c}{\mathrm{\Φ}}_1=\mathrm{diag}[e^{{\mathrm{j\θ}}_1},\·\·\·,e^{{\mathrm{j\θ}}_D}]\\ {\mathrm{\Φ}}_2=\mathrm{diag}[e^{{-\mathrm{j\θ}}_1},\·\·\·,e^{{-\mathrm{j\θ}}_D}]\end{array}\right.---\left(9\right)$

1.3.基于广义声能流的声压振速互协方差矩阵

新的复数域双振速输出分量对各方向的入射信号的幅度响应都是相同的，其矢量特点仅体现在相位上，因此能够完整接收到从任何方向入射的信号。另一个特点是新的复数域双振速输出分量的背景噪声与声压通道的背景噪声依然是不相关的，可以从下式看出：

$\left\{\begin{array}{c}E\left(n_p{n}_1^H\right)=E\left(n_p{n}_x^H\right)+\mathrm{jE}\left(n_p{n}_y^H\right)=0\\ E\left(n_p{n}_2^H\right)=E\left(n_p{n}_x^H\right)-\mathrm{jE}\left(n_p{n}_y^H\right)=0\end{array}\right.---\left(10\right)$

考虑声压与新的复数域双振速分量的协方差矩阵如下：

$\left\{\begin{array}{c}{R}_1=E\left({\mathrm{pv}}_1^H\right)=A_p{R}_s{\mathrm{\Φ}}_1^H{A}_p^H\\ R_2=E\left(p_2^H\right)=A_p{R}_s{\mathrm{\Φ}}_2^H{A}_p^H\end{array}---\left(11\right)\right.$

式中：R_s=E(ss^H)是D×D维的入射信号协方差矩阵。当入射信号互不相关时，R_s是一个对角阵，对角元素为各信号的功率。由于新的复数域双振速输出分量是对原声矢量阵振速分量的扩展，故将新的声压振速互协方差矩阵称为基于广义声能流的互协方差矩阵。

实际应用中R₁、R₂采用最大似然估计值代替：

$\left\{\begin{array}{c}{R}_1=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{pv}}_1^H\\ R_2=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{pv}}_2^H\end{array}\right.---\left(12\right)$

式中，M是采样快拍数。

1.4.基于矩阵束ESPRIT算法的方位估计

令γ为一常数，考查矩阵束(R₁,R₂)：

$R_1-\mathrm{\γ}{R}_2=A_p{R}_s({\mathrm{\Φ}}_1^H-\mathrm{\γ}{\mathrm{\Φ}}_2^H)A_p^H---\left(13\right)$

由于A_p是列满秩的，而R_s非奇异，对(13)式两边求矩阵秩：

$\mathrm{rank}(R_1-\mathrm{\γ}{R}_2)=\mathrm{rank}({\mathrm{\Φ}}_1^H-\mathrm{\γ}{\mathrm{\Φ}}_2^H)---\left(14\right)$

当时，矩阵是非奇异阵，而当时，矩阵秩亏缺成为奇异阵。即均是矩阵束(R₁,R₂)的非零广义特征值。

通常情况下，协方差矩阵R₁、R₂是不满秩的，在求解非零广义特征值的过程中，有可能导致潜在的数值困难。本发明采用奇异值分解(SVD)和总体最小二乘法(TLS)，可以避免这个问题。

对R₁奇异值分解：

$R_1={\mathrm{U\ΣV}}^H=\left[\begin{array}{cc}{U}_1& U_2\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Σ}}_1& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V_1}^H\\ {V_2}^H\end{array}\right]---\left(15\right)$

如前文所述，式中Σ₁是由R₁的非零奇异值构成的对角阵；U₁、V₁分别是非零奇异值对应的左、右奇异向量；U₂、V₂则是零奇异值对应的左、右奇异向量。于是可得：

$U_1^H(R_1-\mathrm{\γ}{R}_2)V_1={\mathrm{\Σ}}_1-\mathrm{\γ}{U}_1^H{R}_2{V}_1---\left(16\right)$

可见，N×N维矩阵束(R₁,R₂)的D个非零广义特征值即是D×D维矩阵束的全部广义特征值。

本发明的这种方位估计方法，实现方法如图3所示，具体计算步骤总结如下：

(1)对阵列的接收信号进行窄带滤波，获得待处理频点上的窄带输出信号。

(2)依据式(8)计算新的复数域双振速输出分量。

(3)依据式(12)估计基于广义声能流的声压振速互协方差矩阵。

(4)对协方差矩阵进行奇异值分解，并估计出非零奇异值矩阵Σ₁和对应的左、右奇异向量矩阵U₁、V₁。

(5)对矩阵束进行广义特征值分解，从而得到入射信号的方位估计值。

2.方位估计的具体计算实例

5元声矢量阵阵元位置如图2所示，入射远场窄带平面波的中心频率为1kHz。

2.1.目标方位估计的均方根误差仿真实例

同时采用本发明所提出的方位估计方法和文献[4]中的方法进行了500次仿真实验，分别进行了单目标和双目标方位估计，每次实验的采样快拍数为100，单目标情况下，目标水平方位角为30°；双目标情况下，目标水平方位角分别为30°、40°。

首先，对阵列的接收数据进行窄带滤波，并据式(8)计算新的复数域双振速输出分量v₁、v₂，图4、图5分别显示了单目标和双目标情况时，信噪比(SNR)为0dB时各通道的输出频谱。由图可以看出，由于自身的滤波特性，常规的振速输出分量的频谱峰值要比声压通道小；而复数域双振速分量具有全向的接收能力，其频谱峰值与声压通道基本一致。

然后，根据式(12)估计出基于广义声能流的协方差矩阵对R₁、R₂，最后应用矩阵束ESPRIT方法计算出矩阵对R₁、R₂之间的旋转不变因子，获得方位估计值。图6、图7分别给出了两种情况下目标水平方位角估计的均方根(RMS)误差和SNR的关系曲线。图中横坐标为SNR，图6的纵坐标为单目标估计的RMS误差，图7的纵坐标为双目标估计的RMS误差之和。

由图6和图7可以看出，在对数坐标下，两种估计方法的均方根误差与信噪比大致成线性关系，本发明所提出的方法相较于文献[4]的方法始终具有更低的估计误差。在图6所示的单目标情况下，当SNR=5dB时，文献[4]的方法均方根误差降至1°左右，而本发明所提出的方法在SNR=0dB时，就已降至1°；在图7所示的双目标情况下，当SNR=17dB时，文献[4]的方法均方根误差降至1°左右，而本发明所提出的方法在SNR=10dB时，就能降至1°。可见，本发明所提出的方法有效的结合了矢量水听器抗各向同性噪声的能力，因而具有更低的可处理信噪比门限。

2.2.y轴附近方向入射信号的方位估计仿真实例

考虑单目标入射的情况，信号的入射方向以0.5°的步长从87°变化到93°，信噪比SNR=0dB，采样快拍数为100。分别采用常规的振速输出分量和复数域双振速输出分量构造声压振速互协方差矩阵，并进行方位估计。图8显示了目标从y轴方向入射时，各通道的输出信号频谱，图9显示了这两种方案的方位估计结果。

由图8和图9可以看出，当信号从y轴方向入射时，振速通道v_x已彻底接收不到信号，因此采用声压与常规振速分量的ESPRIT方位估计结果失效，而复数域双振速通道v₁、v₂能完整接收到信号，依然能给出高精度的方位估计结果。

Claims (3)

1.一种声矢量阵高精度远程方位估计方法，其特征是：

(1)对二维矢量水听器阵列的接收信号进行窄带滤波，获得待处理频点上的窄带输出信号；

(2)将二维矢量水听器阵列相互正交的两个振速分量在复数域进行线性组合，转换为两个新的振速输出分量，即构造复数域双振速分量,具体包括：将y方向的振速输出分量v_y与复数单位j相乘，再与x方向的振速输出分量v_x相加获得新的振速输出分量v₁；将y方向的振速输出分量v_y与负的复数单位-j相乘，再与x方向的振速输出分量v_x相加获得新的振速输出分量v₂；

(3)采用最大似然估计求出一定数量采样快拍下的声压及复数域双振速分量的互协方差矩阵对，即构造出基于广义声能流的声压振速互协方差矩阵对；

(4)应用矩阵束的ESPRIT算法计算基于广义声能流的声压振速互协方差矩阵对之间的旋转不变因子，从而进行方位估计。

2.根据权利要求1所述的声矢量阵高精度远程方位估计方法，其特征是所述的构造出基于广义声能流的声压振速互协方差矩阵对具体包括：将一定采样快拍数量的声压传感器阵列的输出分别与两个新的复数域振速输出分量的共轭转置相乘，并求其平均值，得到基于广义声能流的互协方差矩阵对R₁和R₂。

3.根据权利要求1或2所述的声矢量阵高精度远程方位估计方法，其特征是：所述应用矩阵束的ESPRIT算法计算基于广义声能流的声压振速互协方差矩阵对之间的旋转不变因子中采用奇异值分解和总体最小二乘法将一个较大维数的病态最小二乘问题转化为一个较小维数的非病态总体最小二乘问题。