CN109696657B - 一种基于矢量水听器的相干声源定位方法 - Google Patents

Abstract

本发明公开了一种基于矢量水听器的相干声源定位方法，使用矢量水听器阵列采集声源信号，并获得采集信号的协方差矩阵；利用互不重叠子阵列抽取数据的方法进行信号解相干；利用最小二乘方法估计解相干后协方差矩阵的信号子空间，根据矢量水听器阵列流型旋转不变特性进行求解，最终得到多个相干声源的二维DOA估计。本发明方法估计相干声源方位准确度高，易于实现，还能够有效克服方位模糊问题和相干声源彼此干扰问题。

Description

一种基于矢量水听器的相干声源定位方法

技术领域

本发明涉及相干声源情况下矢量水听器阵列声源定位领域，具体为一种基于矢量水听器的相干声源定位方法。

背景技术

矢量水听器是由在水下声场中正交分布的两到三个振速水听器加上一个声压水听器组成。由于其比传统声压水听器优越的识别能力，矢量水听器阵列在水下信号处理中发挥了重要作用，并且在水下识别、海湾声学反演和水下通信(见文献：Song A,Abdi A,Badiey M,et al.Experimental Demonstration of Underwater AcousticCommunication by Vector Sensors[J].IEEE Journal of Oceanic Engineering,2011,36(3):454-461.)等领域中得到了广泛的应用。在过去的十年中，许多基于矢量水听器的子空间技术被提出，例如MUSIC和ESPRIT，还有使用矢量水听器阵列来估计2D水下信号方位(见文献：He J,Liu Z.Efficient underwater two-dimensional coherent sourcelocalization with linear vector-hydrophone array[J].Signal Processing,2009,89(9):1715-1722.)。

以上所述方法均采用非相干信号，即信号协方差矩阵具有全秩。然而，这种假设在多路径传播或蓄意干扰而存在相干或者高度相关的情况下往往不适用。相干信号可以降低入射信号协方差矩阵的秩，从而严重降低技术的性能，无法正确估计声源的位置。因此，学者们针对声源相干问题进行了大量的研究并提出了如最大似然法、空间平滑技术和Toeplitz方法等解决方法。

为了用矢量水听器阵列处理相干信号，文献(Tao J,Chang W,Shi Y.Direction-finding of coherent sources via'particle-velocity-field smoothing'[J].IetRadar Sonar&Navigation,2008,2(2):127-134.)提出了通过数据相关矩阵的矢量平滑技术来恢复信号子空间的秩。然而，这种矢量平滑技术需要几何平面阵列或2D迭代搜索来估计入射信号的二维方向。文献(Liu S,Yang L,Xie Y,et al.2D DOA Estimation forCoherent Signals with Acoustic Vector-Sensor Array[J].Wireless PersonalCommunications,2017,95(2):1285-1297)提出了一种使用ESPRIT估计入射信号仰角，通过修正的阵列流型匹配得到相干信号的方位角的方法。该方法相对文献(Gu J F,Wei P,TaiH M.2-D direction-of-arrival estimation of coherent signals using cross-correlation matrix[J].Signal Processing,2008,88(1):75-85.)中的计算过程操作简便，定位更加准确。

上述方法都是分步计算，即先得到俯仰角再根据俯仰角求得方位角，进而得到相干声源的方位。由此可知，俯仰角的估计误差会影响方位角的准确估计，定位精度将受到了极大影响。

发明内容

本发明针对现有的矢量水听器对相干声源定位的不足，提出一种基于矢量水听器的相干声源定位方法，以解决谱峰搜索带来的计算量和分步估计带来的误差等问题。本发明适用于任意结构的矢量水听器阵列，能够实现角度估计自动配对。

为解决上述技术问题，本发明提供了一种基于矢量水听器的相干声源定位方法，包括如下步骤：

步骤1：使用M元矢量水听器阵列采集K个声源信号，建立矢量水听器阵列接收数据模型；

步骤2：根据接收信号求得协方差矩阵，利用互不重叠子阵列抽取数据的方法进行信号解相干；

步骤3：使用变换矩阵，提取解相干后阵列流型中的声压矢量和振速矢量；

步骤4：利用最小二乘方法估计解相干后协方差矩阵的信号子空间；

步骤5：根据矢量水听器阵列流型旋转不变特性进行求解得到相干声源的二维DOA估计。

步骤1中矢量水听器阵列接收数据模块建立具体步骤为：

假设K个波长为λ的窄带平面波信号s_i(t),i＝1,2,...,K从远场入射到M元矢量水听器阵列上，且声场中各噪声之间、噪声与信号之间互不相关，K个信号中有L_max个相干声源，设声源的俯仰角为

方位角为φ_i,i＝1,2,...,K，则矢量水听器的阵列流型矢量为：

式中，

每个矢量水听器阵列接收数据模型为：

式中，

表示第m个阵元的阵列流型矢量，n_m(t)表示第m个阵元接收到的噪声，m＝1,2,...,M、t＝1,2,...,N、ψ_i＝(2πd/λ)γ_i，i＝1,2,...K，其中M表示阵元数，N表示信号的快拍数，K表示声源个数。

步骤2具体如下：

由式(2)可得阵列接收信号：

x(t)＝As(t)+n(t)，t＝1,2,...,N (3)

式中，

阵列接收信号的相关矩阵为：

式中，

表示信号的相关矩阵；

把相关矩阵先分为L_max个子阵，L_max表示相干声源的个数，每个子阵的维度为4(M-L_max+1)×4M，第l个子阵标记为R_l，l＝1,2,…l_max，由R_x的第4(l-1)+1行到4(M-L_max+1)行组成，通过这些子阵，构造一个新的矩阵R：

式中，矩阵R的维数为4(M-L_max+1)×4ML_max；

根据阵列接收信号的相关矩阵：

带入矩阵R，可以得到：：

式中，

由于

所以当子阵阵元数目≥K时解相干后的数据协方差矩阵R是满秩的。

步骤3具体如下：

使用转换矩阵J将矩阵R阵列流型中的声压矢量和振速矢量提取出来，得到一个新矩阵

可以表示为

式中，

J＝[J₁ J₂ J₃ J₄]，

e_i是第i个分量为1其他都为零的4(M＝L_max+1)×1单位向量；

分段表示转换后的阵列流型矩阵

可以分段写成：

式中，

由此得A_j和A₄的关系：

A_j＝A₄Γ_j，j＝1,2,3 (10)

式中，Γ₁＝diag{α₁,α₂,…α_K}，Γ₂＝diag{β₁,β₂,…β_K}，Γ₃＝diag{γ₁,γ₂,…γ_K}，三个矩阵都为K×K对角矩阵；

通过估计矩阵Г_j中元素的值得到第i个声源的方位

步骤4具体如下：

矩阵

的K个最大特征向量对应的特征矢量与K个声源的阵列流型矢量是线性关系，由此可以得到：

式中

U_j＝A_jT，j＝1,2,3,4 (14)

从U₄＝A₄T可以推导出A₄＝U₄T^-1。由U₁＝A₁T＝A₄Γ₁T可以得到U₁＝U₄T^-1Γ₁T，相似地可以得到，U₂＝U₄T^-1Γ₂T，U₃＝U₄T^-1Γ₃T；

我们定义

Λ_j＝T^-1Γ_jT，j＝1,2,3 (15)

则Λ_j的特征值是矩阵Γ_i的对角元素。所以构造矩阵Λ_i的估计值，计算它的特征值，就可以得到信号的方位

式(14)可以写成

U_j＝U₄Λ_j (16)

从上式可以推出

矩阵

的K个最大特征值对应的特征向量构成信号子空间特征向量；

假设U_j和Λ_j的估计值分别为

和

利用最小二乘方法估计信号子空间：

得到

的估计值为：

步骤5具体为：

先求得

的特征值中可以得到

令：

得到每个声源的方位角和俯仰角：

有益效果：本发明与现有技术相比，具有以下优点：

(1)本发明采用线性矢量水听器阵列，相比常见L型阵列，阵列结构更为简便。

(2)本发明采用互不重叠的子阵提取方法解相干，在保留较大阵列孔径的同时，减少相干声源对定位估计的影响。

(3)本发明充分利用矢量水听器的阵列流型结构及其旋转不变特性，确保定位估计效果更稳定。

(4)本发明的实施案例表明，本发明比传统方法定位效果更好。

附图说明

图1是本发明的流程图。

图2是本发明的阵列结构图。

图3是本发明在三个不同方位的相干声源定位结果仿真图。

图4是本发明与PM、ESPRIT-AMM方法不同SNR与RMSE的对比关系。

图5是本发明与PM、ESPRIT-AMM方法不同快拍数与RMSE的对比关系。

图6是本发明在不同阵元数时不同SNR与RMSE的关系。

具体实施方式

下面结合附图和实施例进一步阐述本发明。

如图1所示，本发明包括以下步骤：

步骤1：使用M元矢量水听器阵列采集K个声源信号，建立矢量水听器阵列接收数据模型；

假设K个波长为λ的窄带平面波信号s_i(t),i＝1,2,...,K从远场入射到M元矢量水听器阵列上，且声场中各噪声之间、噪声与信号之间互不相关，如图2所示。K个信号中有L_max个相干信号，设信号的俯仰角为

方位角为φ_i,i＝1,2,...,K，则矢量水听器的阵列流型矢量为

式中，

每个矢量水听器的接收数据模型为：

式中，

表示第m个阵元的阵列流型矢量，n_m(t)表示第m个阵元接收到的噪声，m＝1,2,...,M、t＝1,2,...,N、ψ_i＝(2πd/λ)γ_i，i＝1,2,...K，其中M表示阵元数，N表示信号的快拍数，K表示声源个数。

步骤2：根据接收信号求得协方差矩阵，利用互不重叠子阵列抽取数据的方法进行信号解相干；

由式(2)可得阵列接收信号：

x(t)＝As(t)+n(t)t＝1,2,...,N (3)

式中，

阵列接收信号的相关矩阵为：

式中，

表示信号的相关矩阵；

为了能够将相干信号解相干，把相关矩阵先分为L_max个子阵，每个子阵的维度为4(M-L_max+1)×4M。第l个子阵标记为R_l,l＝1,2,…l_max，由R_x的第4(l-1)+1行到4(M-L_max+1)行组成。通过这些子阵，我们构造一个新的矩阵R：

式中，矩阵R的维数为4(M-L_max+1)×4ML_max。

使用式(4)和式(5)矩阵R可以写成分段形式，如下：

式中，

由于

所以当子阵阵元数目≥K时解相干后的数据协方差矩阵R是满秩的。

步骤3：使用转换矩阵，提取解相干后阵列流型中的声压矢量和振速矢量；

使用转换矩阵J将矩阵R阵列流型中的声压矢量和振速矢量提取出来，得到一个新矩阵

可以表示为

式中，

J＝[J₁ J₂ J₃ J₄]，

e_i是第i个分量为1其他都为零的4(M-L_max+1)×1单位向量。由于矩阵J的每一列与其他列正交，所以矩阵J的秩等于4(M-L_max+1)，根据式(7)可以直接得出矩阵

的秩等于K。通过对矩阵

的特征分解，可以得到K个正交向量构成信号子空间，即

列向量的线性空间。

阵列流型矩阵

可以分段写成

式中，

由此可以得到A_j和A₄的关系：

A_j＝A₄Γ_j，j＝1,2,3 (10)

式中，Γ₁＝diag{α₁,α₂,…α_K}，Γ₂＝diag{β₁,β₂,…β_K}，Γ₃＝diag{γ₁,γ₂,…γ_K}，三个矩阵都为K×K对角矩阵。式(10)表明了每个矩阵对(A_j,A₄)的关系，因此，可以通过估计矩阵Γ_i中元素的值得到第i个声源的方位

步骤4：利用最小二乘方法估计解相干后协方差矩阵的信号子空间；

使用特征分解可以分为两个正交子空间，一个是K维信号子空间，包含K个最大特征值对应的特征向量；另一个是[4(M-L_max+1)-K]维噪声子空间。通过对矩阵

的特征分解，得到噪声子空间和信号子空间。设信号子空间为U_s，则有

所以，信号子空间U_s可以表示为

T为K×K的非奇异矩阵。从式(12)可以看出，矩阵

的K个最大特征向量对应的特征矢量与K个声源的阵列流型矢量是线性关系。由式(9)和(12)可以得到

式中

U_j＝A_jT，j＝1,2,3,4 (14)

从U₄＝A₄T可以推导出A₄＝U₄T^-1。由U₁＝A₁T＝A₄Γ₁T可以得到U₁＝U₄T^-1Γ₁T，相似地可以得到，U₂＝U₄T^-1Γ₂T，U₃＝U₄T^-1Γ₃T。

我们定义

Λ_j＝T^-1Γ_jT，j＝1,2,3 (15)

则Λ_j的特征值是矩阵Γ_i的对角元素。所以构造矩阵Λ_i的估计值，计算它的特征值，就可以得到信号的方位

式(14)可以写成

U_j＝U₄Λ_j (16)

从上式可以推出

矩阵

的K个最大特征值对应的特征向量构成信号子空间特征向量。假设U_j和Λ_j的估计值分别为

和

通过信号子空间特征向量估计得到的U₄并不精确，所以采用LS-ESPRIT对式(16)进行估计

得到

的估计值为

步骤5：根据矢量水听器阵列流型旋转不变特性进行求解得到相干声源的二维DOA估计。

从

的特征值中可以得到

因此，α_i,β_i,γ_i的值可以表示为

最后，得到每个声源的方位角和俯仰角

实施结果，如图3-6所示：

以下实施例所用参数如下：

图3是三个不同方位的相干声源定位结果仿真图。矢量水听器阵列个数为10，阵元间隔为d＝λ/2，其中λ表示声波的波长，3个相干信源分别来自

分别在SNR＝0dB和SNR＝20dB情况下的100次仿真结果。

图4是本发明与PM、ESPRIT-AMM方法不同SNR与RMSE的对比关系。矢量水听器阵列个数为10，阵元间隔为d＝λ/2，其中λ表示声波的波长。角度估计的均方根误差为

J表示蒙特卡罗试验次数，K表示声源的个数，

和

表示j次试验中第k个目标的到达角估计值。取快拍数为500，信噪比范围为0dB～20dB，蒙特卡罗实验次数为200次时，各个方法定位结果的均方误差，验证方法的定位估计性能。

图5是本发明与PM、ESPRIT-AMM方法不同快拍数与RMSE的对比关系。矢量水听器阵列个数为10，阵元间隔为d＝λ/2，其中λ表示声波的波长。取信噪比为20dB，快拍数为100～1000，蒙特卡罗实验次数为200次时，各个方法定位结果的均方误差，验证各个方法的定位估计性能。

图6是本发明在不同阵元数时不同SNR与RMSE的关系。矢量水听器阵列个数为10，阵元间隔为d＝λ/2，其中λ表示声波的波长。在快拍数为500，信噪比为20dB，蒙特卡罗实验次数为200次，阵元数分别为5、8、10时，对两个相干声源定位结果进行仿真对比。

Claims (5)

1.一种基于矢量水听器的相干声源定位方法，其特征在于：包括如下步骤：

步骤1：使用M元矢量水听器阵列采集K个声源信号，建立矢量水听器阵列接收数据模型；

步骤2：根据接收信号求得协方差矩阵，利用互不重叠子阵列抽取数据的方法进行信号解相干；

步骤3：使用变换矩阵，提取解相干后阵列流型中的声压矢量和振速矢量；

步骤4：利用最小二乘方法估计解相干后协方差矩阵的信号子空间；

步骤5：根据矢量水听器阵列流型旋转不变特性进行求解得到相干声源的二维DOA估计；

步骤1中矢量水听器阵列接收数据模块建立具体步骤为：

假设K个波长为λ的窄带平面波信号s_i(t),i＝1,2,...,K从远场入射到M元矢量水听器阵列上，且声场中各噪声之间、噪声与信号之间互不相关，K个信号中有L_max个相干声源，设声源的俯仰角为

方位角为φ_i,i＝1,2,...,K，则矢量水听器的阵列流型矢量为：

式中，

每个矢量水听器阵列接收数据模型为：

式中，

表示第m个阵元的阵列流型矢量，n_m(t)表示第m个阵元接收到的噪声，m＝1,2,...,M、t＝1,2,...,N、ψ_i＝(2πd/λ)γ_i，i＝1,2,...K，其中M表示阵元数，N表示信号的快拍数，K表示声源个数。

2.根据权利要求1所述的一种基于矢量水听器的相干声源定位方法，其特征在于，所述步骤2具体如下：

由式(2)可得阵列接收信号：

x(t)＝As(t)+n(t)，t＝1,2,...,N (3)

式中，

阵列接收信号的相关矩阵为：

式中，

表示信号的相关矩阵；

把相关矩阵先分为L_max个子阵，L_max表示相干声源的个数，每个子阵的维度为4(M-L_max+1)×4M，第l个子阵标记为R_l,l＝1,2,…l_max，由R_x的第4(l-1)+1行到4(M-L_max+1)行组成，通过这些子阵，构造一个新的矩阵R：

式中，矩阵R的维数为4(M-L_max+1)×4ML_max；

根据阵列接收信号的相关矩阵：

带入矩阵R，可以得到：

式中，

由于

所以当子阵阵元数目≥K时解相干后的数据协方差矩阵R是满秩的。

3.根据权利要求2所述的一种基于矢量水听器的相干声源定位方法，其特征在于，所述步骤3具体如下：

使用转换矩阵J将矩阵R阵列流型中的声压矢量和振速矢量提取出来，得到一个新矩阵

可以表示为

式中，

J＝[J₁ J₂ J₃ J₄]，

e_i是第i个分量为1其他都为零的4(M-L_max+1)×1单位向量；

分段表示转换后的阵列流型矩阵

可以分段写成：

式中，

由此可以得到A_j和A₄的关系：

A_j＝A₄Γ_j，j＝1,2,3 (10)

式中，Γ₁＝diag{α₁,α₂,…α_K}，Γ₂＝diag{β₁,β₂,…β_K}，Γ₃＝diag{γ₁,γ₂,…γ_K}，三个矩阵都为K×K对角矩阵；

通过估计矩阵Γ_j中元素的值得到第i个声源的方位

4.根据权利要求3所述的一种基于矢量水听器的相干声源定位方法，其特征在于，所述步骤4具体如下：

矩阵

的K个最大特征向量对应的特征矢量与K个声源的阵列流型矢量是线性关系，由此可以得到：

式中

U_j＝A_jT，j＝1,2,3,4 (14)

从U₄＝A₄T可以推导出A₄＝U₄T^-1，由U₁＝A₁T＝A₄Γ₁T可以得到U₁＝U₄T^-1Γ₁T，相似地可以得到，U₂＝U₄T^-1Γ₂T，U₃＝U₄T^-1Γ₃T；

我们定义

Λ_j＝T^-1Γ_jT，j＝1,2,3 (15)

则Λ_j的特征值是矩阵Γ_i的对角元素，所以构造矩阵Λ_i的估计值，计算它的特征值，就可以得到信号的方位

式(14)可以写成

U_j＝U₄Λ_j (16)

从上式可以推出

矩阵

的K个最大特征值对应的特征向量构成信号子空间特征向量；

假设U_j和Λ_j的估计值分别为U_j和Λ_j，利用最小二乘方法估计信号子空间：

得到

的估计值为：

5.根据权利要求4所述的一种基于矢量水听器的相干声源定位方法，其特征在于：

所述步骤5具体为：

先求得

的特征值中可以得到

令：

得到每个声源的方位角和俯仰角：

Images