CN112799008B - 一种声速无关的快速二维波达方向估计方法 - Google Patents

Abstract

本发明公开了一种声速无关的快速二维波达方向估计方法，步骤如下：基于三平行线阵得到三个等阵元间距的子阵列接收数据模型；对子阵列接收数据模型处理得自协方差矩阵和互协方差矩阵，联立得协方差扩展矩阵；对扩展矩阵处理得传播算子矩阵，并基于传播算子矩阵与扩展矩阵间的线性关系得三平行线阵对应含声速信息的特征值参数；利用特征值参数间的角度关系消去声速变量，并求出目标源信号的二维波达方向估计值。本方法基于阵列接收数据的协方差信息实现矩阵扩展，同时以线性变换代替特征分解的方式对扩展矩阵处理，大大降低计算量。此外在计算过程中消除信号在介质中的传播速度，使波达方向估计值与速度无关，排除了速度对算法结果的影响。

Description

一种声速无关的快速二维波达方向估计方法

技术领域

本发明涉及目标定位技术领域，具体涉及一种声速无关的快速二维波达方向估计方法。

背景技术

在阵列信号处理技术中，水下波达方向估计(DOA)占有重要地位。它通过对接收信号进行处理，从而获得信号源的位置信息。

针对二维水下波达方向估计(DOA)，许多高分辨率算法被提出来，例如ML、2D-ESPRIT、2D-MUSIC、2D-ROOT-MUSIC等。这些算法都假定目标信号在介质中的传播速度是常量，没有考虑速度对算法结果的影响。所以在水下环境复杂、速度变化大、未知的情况下容易产生较大误差，降低目标定位的精确度。

在考虑去除声速影响的算法中，利用阵列几何关系的DOA估计方法被提出，如专利申请201822100323.4、201822100446.8。该类专利使用特定的阵列结构如三维正交阵、夹角可调三维阵作为接收阵列，实现对ESPRIT、MUSIC的改进。

基于特征值参数间的角度关系，信号在介质中的传播速度被消去，使最终结果与速度无关。但该类专利所提方法利用阵列接收数据的奇异值分解(SVD)和互相关信息的特征值分解(EVD)获得噪声子空间和信号子空间，计算复杂度高，不利于对信源的快速实时定位。

发明内容

本发明的目的是为了解决现有技术中存在的计算复杂度高、估计值受声速影响等不足，提供一种声速无关的快速二维波达方向估计方法。本方法基于阵列接收数据的协方差信息实现矩阵扩展，使阵列接收数据得到充分利用，同时以线性变换代替特征分解的方式对扩展矩阵处理，大大降低计算量。此外在估计过程中消除信号在介质中的传播速度，使波达方向估计值与速度无关，排除了速度对算法结果的影响。

本发明的目的可以通过采取如下技术方案达到：

一种声速无关的快速二维波达方向估计方法，所述的估计方法包括以下步骤：

S1、基于三平行线阵得到三个等阵元间距的子阵列接收数据模型；

S2、对子阵列接收数据模型处理得自协方差矩阵和互协方差矩阵，联立得协方差扩展矩阵；

S3、对协方差扩展矩阵处理得传播算子矩阵，并基于传播算子矩阵与协方差扩展矩阵间的线性关系得三平行线阵对应含声速信息的特征值参数；

S4、利用特征值参数间的角度关系消去声速变量，并求出目标源信号的二维波达方向估计值。

进一步地，所述步骤S1中的三平行线阵由三个均匀线性子阵相互平行排布而成，三个均匀线性子阵列分别记为子阵1、子阵2和子阵3，其中，子阵1排布于坐标系x轴上，子阵2排布于坐标系x-o-y平面上，子阵3排布于坐标系x-o-z平面上，子阵2、子阵3与子阵1的平行间距为d，各子阵均由M个接收阵元组成，阵元间距均为d，以坐标系原点O作为参考阵元，则总阵元数是3M，假定目标源信号的个数为K，其中心频率为f_i，声速为c_i，i＝1,2,…,K，且相互独立，对应入射到阵列的方向角设为θ_i，俯仰角设为

子阵1、子阵2和子阵3的接收数据模型X(l)、Y(l)、Z(l)分别表示为:

X(l)＝A_xS(l)+N_x(l) (1)

Y(l)＝A_xVS(l)+N_y(l) (2)

Z(l)＝A_xUS(l)+N_z(l) (3)

其中，是M×K维导向向量矩阵，/>为第i个源信号入射到子阵1的导向向量，/> S(l)＝[s₁(l),s₂(l),…,s_i(l)…,s_K(l)]^T为K×L维的目标源矩阵，其中，s_i(l)表示第i个目标源信号，[·]^T表示为矩阵的转置，N_x(l)、N_y(l)、N_z(l)均为M×L维的噪声矩阵，l＝1，2，…，L，l为快拍数。

进一步地，所述步骤S2中对子阵列接收数据模型处理得自协方差矩阵R_xx和互协方差矩阵R_xy、R_xz、R_yx、R_zx：

其中，E[·]表示矩阵的期望，[·]^H表示矩阵的共轭转置，[·]^-1表示矩阵的逆，R_s为S(l)的协方差矩阵，为噪声信号的功率，I_M为秩为M的单位矩阵；

对自协方差矩阵R_xx进行特征值分解，得的估计值/>定义无噪声的自协方差矩阵/>联立协方差矩阵得协方差扩展阵列R：

进一步地，所述步骤S3中得到三平行线阵对应含声速信息的特征值参数的过程如下：

以线性变换的方式对协方差扩展阵列R处理，由R得自协方差矩阵R_w＝E[RR^H]，利用R_w可得K×(5M-K)维的传播算子矩阵P，具体如下所示：

其中，R_w＝[R_w1,R_w2],R_w1为5M×K维矩阵，R_w2为5M×(5M-K)维矩阵，定义扩展传播算子矩阵P_e＝[I_K,P]^H,I_K为秩为K的单位矩阵，由P_e与式(9)中A的线性关系得A＝P_e×A₁，其中，A₁为A的前K行，即K×K维矩阵；

定义筛选矩阵I_v1＝[I_v11,O_2M×M,I_v12,O_2M×2M]、I_v2＝[O_2M×2M,I_v11,O_2M×M,I_v12]、I_u1＝[O_2M×M,I_2M,O_2M×2M]、I_u2＝[O_2M×2M,I_2M,O_2M×M]，其中，/>I₅、I_M-1、I_2M分别是秩为5、M-1、2M的单位矩阵，O_(M-1)×1、O_M×M、O_2M×M、O_2M×2M分别为(M-1)×1、M×M、2M×M、2M×2M维的零矩阵，/>为克罗内克积，利用筛选矩阵I_x1、I_x2、I_v1、I_v2、I_u1、I_u2可得关系式：

I_x1P_eA₁φ_x＝I_x2P_eA₁ (11)

I_v1P_eA₁V＝I_v2P_eA₁ (12)

I_u1P_eA₁U＝I_u2P_eA₁ (13)

其中，

对式(11)阵列变换可得：

其中，表示为矩阵的伪逆，对Π_x特征值分解得到对角阵φ_x的估计值/>A₁的估计值/>同样的，基于式(12)、(13)可得V、U的估计值/>

由上式可知，通过使用A₁的估计值使得V、U的估计值/>中的元素一一对应，进而实现方向角和俯仰角的自动匹配。

进一步地，所述步骤S4中得到与声速无关的目标源信号的方向角估计值和俯仰角估计值/>

其中，表示对角阵/>第i行和第i列元素，/>表示对角阵/>第i行和第i列元素，/>表示对角阵/>第i行和第i列元素，angle(·)表示对复数求相位角。

本发明相对于现有技术具有如下的优点及效果：

1、本发明基于平行线阵特性，对阵列接收数据模型处理得到互协方差矩阵和自协方差矩阵，并联立得到扩展的协方差矩阵。该协方差扩展矩阵充分利用了接收阵列数据，从而在保持阵元数不变的条件下，进一步提高DOA估计方法的精确度。

2、本发明相比于使用传统DOA估计方法进行水下波达方向估计的方法更具实用性，并且计算复杂度更低。在复杂的水声环境中，目标声源反射回来的声波速度是未知的，而传统DOA估计方法假定声速是不变，没有考虑声速的不确定性对估计算法的影响，使得在复杂水声环境下传统DOA估计方法不具有实用性。本发明方法利用得到的特征值参数间的角度关系去除声速，从而排除了声速的不确定性对波达方向估计的影响，可在复杂水声环境下进行测量，使得其实用性更高。同时，本发明方法以线性变换代替传统特征值分解的方式对阵列接收数据处理并实现了方向角估计值和俯仰角估计值的自动匹配，有效降低了计算复杂度。

附图说明

图1是本发明实施例中三平行均匀线阵的结构示意图；

图2是本发明中公开的声速无关的快速二维波达方向估计方法流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

本实施例提供了一种声速无关的快速二维波达方向估计方法。本方法基于阵列接收数据的协方差信息实现矩阵扩展，使阵列接收数据得到充分利用，同时以线性变换代替特征分解的方式对扩展矩阵处理，大大降低计算量。此外在估计方法中消除信号在介质中的传播速度，使波达方向估计值与速度无关，排除了速度对估计方法结果的影响。

如图2所示，本实施例中声速无关的快速二维波达方向估计方法包括以下步骤：

S1、基于三平行线阵得到三个等阵元间距的子阵列接收数据模型。所使用的三平行均匀线阵如图1所示，三个均匀线性子阵列分别记为子阵1、子阵2和子阵3，其中，线阵1排布于坐标系x轴上，线阵2排布于坐标系x-o-y平面上，线阵3排布于坐标系x-o-z平面上，线阵2、3与线阵1的平行间距为d，各子阵均由M个接收阵元组成，阵元间距均为d，以坐标系原点O作为参考阵元，则总阵元数是3M，假定目标源信号的个数为K，其中心频率为f_i，声速为c_i，i＝1,2,…,K，且相互独立，对应入射到阵列的方向角设为θ_i，俯仰角设为

子阵1、子阵2和子阵3的接收数据模型X(l)、Y(l)、Z(l)分别表示为:

X(l)＝A_xS(l)+N_x(l) (1)

Y(l)＝A_xVS(l)+N_y(l) (2)

Z(l)＝A_xUS(l)+N_z(l) (3)

其中，是M×K维导向向量矩阵，/>为第i个源信号入射到子阵1的导向向量，/> S(l)＝[s₁(l),s₂(l),…,s_i(l)…,s_K(l)]^T为K×L维的目标源矩阵，其中，s_i(l)表示第i个目标源信号，[·]^T表示为矩阵的转置，N_x(l)、N_y(l)、N_z(l)均为M×L维的噪声矩阵，l＝1，2，…，L，l为快拍数。

S2、对子阵列接收数据模型处理得互协方差矩阵和自协方差矩阵，联立得协方差扩展矩阵。

对子阵列接收数据模型处理得自协方差矩阵R_xx和互协方差矩阵R_xy、R_xz、R_yx、R_zx：

其中，E[·]表示矩阵的期望，[·]^H表示矩阵的共轭转置，[·]^-1表示矩阵的逆，R_s为S(l)的协方差矩阵，为噪声信号的功率，I_M为秩为M的单位矩阵。

对自协方差矩阵R_xx进行特征值分解，得的估计值/>定义无噪声的自协方差矩阵/>联立协方差矩阵得协方差扩展阵列R：

S3、对扩展矩阵处理得传播算子矩阵，并基于传播算子矩阵与扩展矩阵间的线性关系得三平行线阵对应含声速信息的特征值参数。

以线性变换代替特征值分解的方式对协方差扩展阵列R处理，本实施例使用传播算子算法对其进行说明。

由R得自协方差矩阵R_w＝E[RR^H]，利用R_w可得K×(5M-K)维的传播算子矩阵P，具体如下所示：

其中，R_w＝[R_w1,R_w2],R_w1为5M×K维矩阵，R_w2为5M×(5M-K)维矩阵。定义扩展传播算子矩阵P_e＝[I_K,P]^H,I_K为秩为K的单位矩阵，由P_e与式(9)中A的线性关系得A＝P_e×A₁，其中，A₁为A的前K行，即K×K维矩阵。

定义筛选矩阵I_v1＝[I_v11,O_2M×M,I_v12,O_2M×2M]、I_v2＝[O_2M×2M,I_v11,O_2M×M,I_v12]、I_u1＝[O_2M×M,I_2M,O_2M×2M]、I_u2＝[O_2M×2M,I_2M,O_2M×M]，其中，/>I₅、I_M-1、I_2M分别为秩为5、M-1、2M的单位矩阵，O_(M-1)×1、O_M×M、O_2M×M、O_2M×2M分别为(M-1)×1、M×M、2M×M、2M×2M维的零矩阵，/>为克罗内克积，利用筛选矩阵I_x1、I_x2、I_v1、I_v2、I_u1、I_u2可得关系式：

I_x1P_eA₁φ_x＝I_x2P_eA₁ (11)

I_v1P_eA₁V＝I_v2P_eA₁ (12)

I_u1P_eA₁U＝I_u2P_eA₁ (13)

其中，

对式(11)阵列变换可得：

其中，表示为矩阵的伪逆，对Π_x特征值分解得到对角阵φ_x的估计值/>A₁的估计值/>同样的，基于式(12)、(13)可得V、U的估计值/>

由上式可知，通过使用A₁的估计值使得V、U的估计值/>中的元素一一对应，进而实现方向角和俯仰角的自动匹配。

S4、利用特征值参数间的角度关系消去声速变量，并求出目标源信号的二维波达方向估计值。

利用特征值参数求得与声速无关的目标源信号的方向角估计值和俯仰角估计值/>

其中，表示为对角阵/>第i行和第i列元素，/>表示为对角阵/>第i行和第i列元素，/>表示为对角阵/>第i行和第i列元素，angle(·)表示为对复数求相位角。

实施例二

本实施例具体公开一种声速无关的快速二维波达方向估计方法，该方法的实施步骤具体如下：

T1、目标信源个数设为K＝5，分别表示为s₁(l),s₂(l),s₃(l),s₄(l),s₅(l)，其中心频率f分别为9.8kHz,9.9kHz,10kHz,10.1kHz,10.2kHz，脉冲长度为5ms。

海水中声速范围大致为1430m/s–1550m/s，则取最小声速1430m/s，最大中心频率10.2kHz，可以求出最小半波长λ/2为7cm。接收模块中的三平行均匀线阵的总阵元个数3M定为18，三个子阵列分别记为子阵1、子阵2和子阵3，且阵元数均为M＝6。两相邻阵元之间的距离必须小于半波长，在满足此限制条件下取阵元间距d＝5cm。

T2、对三平行均匀线阵中的子阵列接收到的目标源信号进行采样；由阵元个数可知，

子阵1接收到的信号为X(l)＝{x₁(l),x₂(l),x₃(l),x₄(l),x₅(l),x₆(l)}，子阵2接收到的信号为Y(l)＝{y₁(l),y₂(l),y₃(l),y₄(l),y₅(l),y₆(l)}，子阵3接收到的信号为Z(l)＝{z₁(l),z₂(l),z₃(l),z₄(l),z₅(l),z₆(l)}。共采样接收200次，并将接收到的信号按照本专利所提供的方法进行处理。

T3、信号的处理分析步骤具体如下：

T31、根据接收阵列的阵列平行特性，得到子阵1、子阵2和子阵3的信号接收数据模型。

T32、对子接收数据模型处理得自协方差矩阵R_xx和互协方差矩阵R_xy、R_xz、R_yx、R_zx，联立得协方差扩展阵列R。

T33、对协方差扩展阵列R处理得传播算子矩阵P，并利用传播算子矩阵与扩展矩阵间的线性关系得三平行线阵对应含声速信息的特征值参数。

T34、利用特征值参数间的角度关系消去声速变量，并求出全部K个目标源信号的方向角估计值和俯仰角估计值。

综上所述，上述实施例为了解决现有技术中存在的计算复杂度高、估计值受声速影响等不足，提供了一种声速无关的快速二维波达方向估计方法。本方法基于阵列接收数据的协方差信息实现矩阵扩展，对扩展矩阵处理时，以线性变换代替特征分解，使阵列接收数据得到充分利用的同时大大降低计算量。此外在算法中消除信号在介质中的传播速度，使波达方向估计值与速度无关，排除了速度对算法结果的影响。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (2)

1.一种声速无关的快速二维波达方向估计方法，其特征在于，所述的估计方法包括以下步骤：

S1、基于三平行线阵得到三个等阵元间距的子阵列接收数据模型；

所述步骤S1中的三平行线阵由三个均匀线性子阵相互平行排布而成，三个均匀线性子阵列分别记为子阵1、子阵2和子阵3，其中，子阵1排布于坐标系x轴上，子阵2排布于坐标系x-o-y平面上，子阵3排布于坐标系x-o-z平面上，子阵2、子阵3与子阵1的平行间距为d，各子阵均由M个接收阵元组成，阵元间距均为d，以坐标系原点O作为参考阵元，则总阵元数是3M，假定目标源信号的个数为K，其中心频率为f_i，声速为c_i，i＝1,2,…,K，且相互独立，对应入射到阵列的方向角设为θ_i，俯仰角设为

子阵1、子阵2和子阵3的接收数据模型X(l)、Y(l)、Z(l)分别表示为:

X(l)＝A_xS(l)+N_x(l) (1)

Y(l)＝A_xVS(l)+N_y(l) (2)

Z(l)＝A_xUS(l)+N_z(l) (3)

其中，是M×K维导向向量矩阵，/>为第i个源信号入射到子阵1的导向向量，/> S(l)＝[s₁(l),s₂(l),…,s_i(l),…，s_K(l)]^T为K×L维的目标源矩阵，其中，s_i(l)表示第i个目标源信号，[·]^T表示为矩阵的转置，N_x(l)、N_y(l)、N_z(l)均为M×L维的噪声矩阵，l＝1，2，…，L,l为快拍数；

S2、对子阵列接收数据模型处理得自协方差矩阵和互协方差矩阵，联立得协方差扩展矩阵；

所述步骤S2中对子阵列接收数据模型处理得自协方差矩阵R_xx和互协方差矩阵R_xy、R_xz、R_yx、R_zx：

其中，E[·]表示矩阵的期望，[·]^H表示矩阵的共轭转置，[·]^-1表示矩阵的逆，R_s为S(l)的协方差矩阵，为噪声信号的功率，I_M为秩为M的单位矩阵；

对自协方差矩阵R_xx进行特征值分解，得的估计值/>定义无噪声的自协方差矩阵联立协方差矩阵得协方差扩展阵列R：

S3、对协方差扩展矩阵处理得传播算子矩阵，并基于传播算子矩阵与协方差扩展矩阵间的线性关系得三平行线阵对应含声速信息的特征值参数；

所述步骤S3中得到三平行线阵对应含声速信息的特征值参数的过程如下：

以线性变换的方式对协方差扩展阵列R处理，由R得自协方差矩阵R_w＝E[RR^H]，利用R_w可得K×(5M-K)维的传播算子矩阵P，具体如下所示：

其中，R_w＝[R_w1,R_w2],R_w1为5M×K维矩阵，R_w2为5M×(5M-K)维矩阵，定义扩展传播算子矩阵P_e＝[I_K,P]^H,I_K为秩为K的单位矩阵，由P_e与式(9)中A的线性关系得A＝P_e×A₁，其中，A₁为A的前K行，即K×K维矩阵；

定义筛选矩阵I_v1＝[I_v11,O_2M×M,I_v12,O_2M×2M]、I_v2＝[O_2M×2M,I_v11,O_2M×M,I_v12]、I_u1＝[O_2M×M,I_2M,O_2M×2M]、I_u2＝[O_2M×2M,I_2M,O_2M×M]，其中，/>I₅、I_M-1、I_2M分别是秩为5、M-1、2M的单位矩阵，O_(M-1)×1、O_M×M、O_2M×M、O_2M×2M分别为(M-1)×1、M×M、2M×M、2M×2M维的零矩阵，/>为克罗内克积，利用筛选矩阵I_x1、I_x2、I_v1、I_v2、I_u1、I_u2可得关系式：

I_x1P_eA₁φ_x＝I_x2P_eA₁ (11)

I_v1P_eA₁V＝I_v2P_eA₁ (12)

I_u1P_eA₁U＝I_u2P_eA₁ (13)

其中，

对式(11)阵列变换可得：

其中，表示为矩阵的伪逆，对Π_x特征值分解得到对角阵φ_x的估计值/>A₁的估计值/>同样的，基于式(12)、(13)可得V、U的估计值/>

由式(15)、(16)可知，通过使用A₁的估计值使得V、U的估计值/>中的元素一一对应，进而实现方向角和俯仰角的自动匹配；

S4、利用特征值参数间的角度关系消去声速变量，并求出目标源信号的二维波达方向估计值。

2.根据权利要求1所述的一种声速无关的快速二维波达方向估计方法，其特征在于，所述步骤S4中得到与声速无关的目标源信号的方向角估计值和俯仰角估计值/>

其中，表示对角阵/>第i行和第i列元素，/>表示对角阵/>第i行和第i列元素，/>表示对角阵/>第i行和第i列元素，angle(·)表示对复数求相位角。