CN113504504B - 一种水下高精度一维doa估计方法 - Google Patents

Abstract

本发明公开了一种水下高精度一维DOA估计方法，步骤如下：采用二维交叉嵌套阵进行一维DOA估计，对两条二级嵌套阵接收数据的协方差矩阵R_xx和R_yy向量化、排序、去冗余，根据矩阵重构算法构建Toeplitz矩阵

和

基于

使用现有的TLS‑ESPRIT算法得到

对

进行特征值分解得到特征值λ_x1，λ_x2，...，λ_xK，基于

使用现有的TLS‑ESPRIT算法得到

对

进行特征值分解得到特征值λ_y1，λ_y2，...，λ_yK，将λ_x1，λ_x2，...，λ_xK从小到大排序，对λ_y1，λ_y2，...，λ_yK从大到小排序，通过排序对两者进行配对；根据二维交叉嵌套阵的几何关系推导出θ_xk、θ_yk和θ_k的关系式，最终求解出目标信号的波达方向角。该发明在得到声速无关的DOA估计表达式的同时，解决了DOA估计过程中参数配对问题，提高了DOA估计精度。

Description

一种水下高精度一维DOA估计方法

技术领域

本发明涉及水下目标定位技术领域，具体涉及一种水下高精度一维DOA估计方法。

背景技术

波达方向估计(Direction of Arrival，DOA)在众多领域已得到广泛应用，而一维水下DOA估计就是指在水面放置传感器阵列利用阵列信号处理技术来对水下声源目标进行一维波达方向估计的方法。与空间目标DOA估计有所不同的是，水下声源目标的DOA估计采用声波作为传播载体，由于声波信号在水下环境传播时，水声信道中的各种障碍物及崎岖不平的海底造成的声波散射作用，导致了信号的衰减。除了水声环境造成信号的快速衰减，水下DOA估计面临的另一个问题就是声速影响。由于河流和海洋等水下环境复杂且不稳定，声波的速度温度、盐度、压力等因素不断变化，水下DOA方法的估计精度受到很大影响，当实际声速偏离预先设定速度，估计精度将因此降低。

现有研究已提出了利用不同构型的二维阵列构造二维角度关系，对声速因子进行消除的一维水下目标DOA估计方法。这类方法提高了未知声速环境下的水下目标DOA估计精度，如专利申请201811241541.8和201811338421.X等。但是在多目标声源的情况下，由于这类方法采用了全组合遍历对多组参数进行配对，使得算法复杂度非常高，装置的估计实时性得不到保证。专利申请201910259714.7中的方法降低水下DOA估计的复杂度，但是估计精度有待提高。因此，如何在消除声速影响的同时提高估计精度成了一个急需解决的问题。

发明内容

本发明的目的是为了解决现有技术中的上述缺陷，提供一种水下高精度一维DOA估计方法，通过对二维交叉嵌套阵的接收信号进行处理，在DOA估计中消除声速这个变量，从而消除水下声速不确定性对目标定位精度的影响。同时由于采用嵌套阵，估计方法的估计精度有了很大的提高，有利于在实际测量中有效估计目标，另外在估计过程中通过排序进行配对，配对复杂度很小。

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

一种水下高精度一维DOA估计方法，所述估计方法包括以下步骤：

S1、采用二维交叉嵌套阵进行一维DOA估计，该用于一维DOA估计的二维交叉嵌套阵的结构如图3所示，与常规一维DOA估计方法采用单条均匀线阵不同，本发明中二维交叉嵌套阵采用两条二级嵌套阵，每条二级嵌套阵上各自分布有(M+N)个阵元，每条二级嵌套阵由两级均匀线阵组成，第一级均匀线阵有M个阵元，阵元间距为d，第二级均匀线阵有N个阵元，阵元间距为(M+1)d，两条嵌套阵连接处有一个公共阵元，两条嵌套阵之间的夹角设为δ，其中

建立坐标系，设一条二级嵌套阵所在直线为x轴，另一条二级嵌套阵所在直线设为y轴，假设只考虑半个平面空间内目标源信号的情况，即x轴的上半平面空间，假设目标信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。假设目标信号的个数为K，目标信号的中心频率为f，实时声速定义为c,目标信号在x轴上的波达方向角为θ_xk,k＝1,2,…,K，在y轴上的波达方向角为θ_yk,k＝1,2,…,K，目标信号的波达方向角设为θ_k,k＝1,2,…,K，且θ_k＝θ_xk，当有K个远场窄带相互独立的信号入射到上述二维交叉嵌套阵，x轴和y轴阵列接收到的信号写成如下的矢量形式：

其中s(t)是一个K×1维的源信号矩阵，n_x(t)和n_y(t)分别为x轴和y轴阵列的噪声矩阵，A_x和A_y分别为x轴和y轴阵列的方向矩阵，写成矢量形式有：

其中a(θ_xk)和a(θ_yk)分别是x轴和y轴阵列的第k(k＝1,2,…,K)个声源的导向向量，

对于两条线阵各自接收的信号x(t)和y(t)，分别求它们的协方差矩阵为R_xx＝E[x(t)x(t)^H]和R_yy＝E[y(t)y(t)^H]；

S2、向量化R_xx和R_yy得到z_xx和z_yy，对z_xx和z_yy进行去冗余并排序，得到

和/>

使用矩阵重构算法基于/>

和/>

构建Toeplitz矩阵/>

和/>

S3、基于

使用TLS-ESPRIT算法得到/>

对/>

进行特征值分解得到特征值λ_x1，λ_x2，...，λ_xK，基于/>

使用TLS-ESPRIT算法得到/>

对/>

进行特征值分解得到特征值λ_y1，λ_y2，...，λ_yK，将λ_x1，λ_x2，...，λ_xK从小到大排序，对λ_y1，λ_y2，...，λ_yK从大到小排序，通过排序对两者进行配对；

S4、根据二维交叉嵌套阵的几何关系推导出θ_xk、θ_yk和θ_k的关系式，最终求解出目标信号的波达方向角。

进一步地，所述步骤S1中对接收信号矩阵进行处理得到协方差矩阵R_xx和R_yy的过程如下：

由于噪声和信号相互独立，得到x轴和y轴阵列接收数据的M×M阶协方差矩阵分别为：

其中R_s＝E[s(t)s^H(t)]为信源部分的协方差矩阵，由于信源相互独立，因此R_s为对角矩阵，I为单位矩阵，σ_n ²为噪声功率，(·)^H表示求共轭转置。

进一步地，所述步骤S2中得到Toeplitz矩阵

和/>

的过程如下：

向量化R_xx和R_yy，得到向量z_xx和z_yy：

其中p为向量化后的信号矩阵，I_n为向量化后的单位矩阵，⊙表示求Khatri-Rao积，()^*表示求共轭；

对向量z_xx和z_yy进行去冗余并排序，得到新的向量

其中

向量化、排序、去冗余后的

和/>

的信号矩阵p的秩为1，根据矩阵重构算法得到

其中

表示/>

的第m个元素，其中/>

同理可以得到

进一步地，所述步骤S3中对λ_x1，λ_x2，…，λ_xK和λ_y1，λ_y2，...，λ_yK进行配对的过程如下：

基于Toeplitz矩阵

和/>

根据现有的TLS-ESPRIT算法得到/>

和/>

形式如下

其中Φ_x和T_x ^-1分别为

的特征值组成的对角矩阵以及对应的特征向量组成的满秩方阵，Φ_y和T_y ^-1分别为/>

的特征值组成的对角矩阵以及对应的特征向量组成的满秩方阵，Φ_x和Φ_y的具体形式为

对

进行特征分解，得到特征值λ_x1，λ_x2，...，λ_xK，对/>

进行特征分解，得到特征值λ_y1，λ_y2，...，λ_yK；

λ_x1，λ_x2，...，λ_xK、λ_y1，λ_y2，...，λ_yK与θ_xk、θ_yk的关系为

此时λ_x1，λ_x2，...，λ_xK、λ_y1，λ_y2，...，λ_yK是没有配对的，考虑二维交叉嵌套阵的阵列结构，通过排序进行配对，对λ_x1，λ_x2，...，λ_xK按照从大到小的顺序排序，对λ_y1，λ_y2，...，λ_yK按照从小到大的顺序排序，排序之后两组特征值的顺序是配对的。

进一步地，所述步骤S4中根据交叉嵌套阵的几何关系推导出θ_xk、θ_yk和θ_k的关系式，最终求解出目标信号的波达方向角的过程如下：

根据二维嵌套阵的阵列结构，得到目标信号相对于x轴阵列的波达角θ_xk、相对于y轴阵列的波达角θ_yk和待估计的波达角θ_k满足：

cosθ_k＝cosθ_xk＝cos(δ-θ_yk) 公式(14)

然后根据公式(13)，得到目标信号相对于x轴阵列的波达方向角θ_xk为

考虑到第k个信号的波达方向角θ_k就是其在x轴上的波达角θ_xk，因此最终的一维DOA估计结果为：

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

1、本发明公开的DOA估计方法适用于其他任意夹角的二维线阵，在阵列构型的选择上具有更大的灵活空间。

2、本发明与利用传统的一维DOA方法进行水下目标波达方向估计的方法相比更具有实用性，估计精确度也更高。传统的一维DOA方法通常假定声速为一个常量，而在实际的复杂水下环境中，声速往往是不断变化的，如果把其当成一个常量来进行计算的话，会导致较大的误差。本发明通过阵列夹角与波达方向角之间的关系消去了声速这个变量，使得最后的运算结果与声速无关，从而提高了估计精度。

3、本发明采用了二级嵌套阵进行DOA估计，由于相同阵元个数下嵌套阵的天线孔径远大于均匀线阵，因此本发明所提方法的估计精度要比现有的水下一维DOA估计方法的估计精度要高。当目标信源数越多时，本发明所提方法的鲁棒性更好，所提方法的优势更明显。

4、本发明方法的实施装置在传统的测量装置上进行了改进，使用夹角可调节的二维交叉嵌套阵，可行性强，安装简单。除此之外，现代处理器计算处理能力的不断提高，这使得本发明所使用的处理器等芯片的集成度高，并且计算能力强，从而保证了本发明的可行性。

附图说明

图1是本发明中公开的水下高精度一维DOA估计方法的实施装置的硬件结构模块图；

图2是本发明中阵列的接收阵元与处理器连接示意图；

图3是本发明中提出的二维交叉嵌套阵阵列模型示意图；

图4是本发明中提出的二级嵌套阵的几何结构示意图；

图5是本发明中提出的二维交叉嵌套阵正面连接示意图；

图6是本发明中提出的二维交叉嵌套阵侧面连接示意图；

图7是本发明中公开的水下高精度一维DOA估计方法的流程图。

具体实施方式

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

实施例一

现有的水下DOA估计方法往往是基于均匀线阵的，这些方法虽然能够得到声速无关的水下DOA估计结果，但是估计精度不够高。考虑到阵元总数相同时，嵌套阵相对于均匀线阵有更大的天线孔径、更高的估计精度，因此本实施例提出基于嵌套阵的一维DOA估计方法。现有的水下DOA估计方法主要有两种思路，一是遍历所有可能的声速值，按照方差最小原则得到声速无关的DOA估计结果；二是构建特殊的接收阵列，利用阵列的几何结构消除声速变量，进而得到声速无关的DOA估计结果。考虑到思路二不需要遍历且计算复杂度较低，本实施例按照思路二设计了一种基于嵌套阵的接收阵列，即二维交叉嵌套阵。基于二维交叉嵌套阵进行水下DOA估计，会得到两组带声速变量的估计值，并且这两组估计值之间存在一种数量关系。接下来需要进行参数配对，现有的配对方法需要遍历所有组合、计算量较大，为了降低计算复杂度，本实施例继续提出利用两组估计值之间的数量关系，通过排序进行配对，得到声速无关的DOA估计结果。本实施例提出的一种水下高精度一维DOA估计方法以二维交叉嵌套阵为接收阵列，根据接收矩阵计算出两组带声速变量的估计值，通过排序对两组估计值进行配对，最终得到高精度的声速无关的DOA估计结果。

如图7所示，本实施例公开了一种水下高精度一维DOA估计方法，该估计方法包括以下步骤：

S1、采用二维交叉嵌套阵进行一维DOA估计，该用于一维DOA估计的二维交叉嵌套阵的结构如图3所示，与常规一维DOA估计方法采用单条均匀线阵不同，本发明中二维交叉嵌套阵采用两条二级嵌套阵，每条二级嵌套阵上各自分布有(M+N)个阵元，每条二级嵌套阵由两级均匀线阵组成，第一级均匀线阵有M个阵元，阵元间距为d，第二级均匀线阵有N个阵元，阵元间距为(M+1)d，两条嵌套阵连接处有一个公共阵元，两条嵌套阵之间的夹角设为δ，其中

建立坐标系，设一条二级嵌套阵所在直线为x轴，另一条二级嵌套阵所在直线设为y轴，假设只考虑半个平面空间内目标源信号的情况，即x轴的上半平面空间，假设目标信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。假设目标信号的个数为K，目标信号的中心频率为f，实时声速定义为c,目标信号在x轴上的波达方向角为θ_xk,k＝1,2,…,K，在y轴上的波达方向角为θ_yk,k＝1,2,…,K，目标信号的波达方向角设为θ_k,k＝1,2,…,K，且θ_k＝θ_xk，当有K个远场窄带相互独立的信号入射到上述二维交叉嵌套阵，x轴和y轴阵列接收到的信号写成如下的矢量形式：

其中s(t)是一个K×1维的源信号矩阵，n_x(t)和n_y(t)分别为x轴和y轴阵列的噪声矩阵，A_x和A_y分别为x轴和y轴阵列的方向矩阵，写成矢量形式有：

其中a(θ_xk)和a(θ_yk)分别是x轴和y轴阵列的第k(k＝1,2,…,K)个声源的导向向量，

对于两条线阵各自接收的信号x(t)和y(t)，分别求它们的协方差矩阵为R_xx＝E[x(t)x(t)^H]和R_yy＝E[y(t)y(t)^H]；

本实施例的步骤S1中对接收信号矩阵进行处理得到协方差矩阵R_xx和R_yy的过程如下：

由于噪声和信号相互独立，得到x轴和y轴阵列接收数据的M×M阶协方差矩阵分别为：

其中R_s＝E[s(t)s^H(t)]为信源部分的协方差矩阵，由于信源相互独立，因此R_s为对角矩阵，I为单位矩阵，σ_n ²为噪声功率，(·)^H表示求共轭转置。

S2、向量化R_xx和R_yy得到z_xx和z_yy，对z_xx和z_yy进行去冗余并排序，得到

和/>

使用矩阵重构算法基于/>

和/>

构建Toeplitz矩阵/>

和/>

本实施例的步骤S2中得到Toeplitz矩阵

和/>

的过程如下：

向量化R_xx和R_yy，得到向量z_xx和z_yy：

其中p为向量化后的信号矩阵，I_n为向量化后的单位矩阵，⊙表示求Khatri-Rao积，()^*表示求共轭；

对向量z_xx和z_yy进行去冗余并排序，得到新的向量

其中

向量化、排序、去冗余后的

和/>

的信号矩阵p的秩为1，根据矩阵重构算法得到

其中

表示/>

的第m个元素，其中/>

同理可以得到

S3、基于

使用TLS-ESPRIT算法得到/>

对/>

进行特征值分解得到特征值λ_x1，λ_x2，...，λ_xK，基于/>

使用TLS-ESPRIT算法得到/>

对/>

进行特征值分解得到特征值λ_y1，λ_y2，...，λ_yK，将λ_x1，λ_x2，...，λ_xK从小到大排序，对λ_y1，λ_y2，...，λ_yK从大到小排序，通过排序对两者进行配对；

本实施例的步骤S3中对λ_x1，λ_x2，...，λ_xK和λ_y1，λ_y2，...，λ_yK进行配对的过程如下：

基于Toeplitz矩阵

和/>

根据现有的TLS-ESPRIT算法得到/>

和/>

形式如下

其中Φ_x和T_x ^-1分别为

的特征值组成的对角矩阵以及对应的特征向量组成的满秩方阵，Φ_y和T_y ^-1分别为/>

的特征值组成的对角矩阵以及对应的特征向量组成的满秩方阵，Φ_x和Φ_y的具体形式为

对

进行特征分解，得到特征值λ_x1，λ_x2，...，λ_xK，对/>

进行特征分解，得到特征值λ_y1，λ_y2，...，λ_yK；

λ_x1，λ_x2，...，λ_xK、λ_y1，λ_y2，...，λ_yK与θ_xk、θ_yk的关系为

此时λ_x1，λ_x2，...，λ_xK、λ_y1，λ_y2，...，λ_yK是没有配对的，考虑二维交叉嵌套阵的阵列结构，通过排序进行配对，对λ_x1，λ_x2，...，λ_xK按照从大到小的顺序排序，对λ_y1，λ_y2，...，λ_yK按照从小到大的顺序排序，排序之后两组特征值的顺序是配对的。

S4、根据二维交叉嵌套阵的几何关系推导出θ_xk、θ_yk和θ_k的关系式，最终求解出目标信号的波达方向角。

本实施例的步骤S4中根据交叉嵌套阵的几何关系推导出θ_xk、θ_yk和θ_k的关系式，最终求解出目标信号的波达方向角的过程如下：

根据二维嵌套阵的阵列结构，得到目标信号相对于x轴阵列的波达角θ_xk、相对于y轴阵列的波达角θ_yk和待估计的波达角θ_k满足：

cosθ_k＝cosθ_xk＝cos(δ-θ_yk) 公式(14)

然后根据公式(13)，得到目标信号相对于x轴阵列的波达方向角θ_xk为

考虑到第k个信号的波达方向角θ_k就是其在x轴上的波达角θ_xk，因此最终的一维DOA估计结果为：

实施例二

本实施例公开了一种水下高精度一维DOA估计方法的具体实施装置，水下高精度一维DOA估计装置包括：数据处理与控制模块、发射模块、接收模块、输出模块和电源模块，如图1和图2所示。

数据处理与控制模块由一对A/D、D/A转换器和一个处理器组成，是整个装置的核心部分，其它所有模块都与它直接相连。它可以控制发射模块，使发射模块发射我们指定的信号；可以控制接收模块的夹角可调线阵，使嵌套阵1保持固定，嵌套阵2以连接点为中心进行自由旋转，从而调节阵列夹角到需要的角度；同时能够对接收模块传过来的信号进行处理，通过本发明的方法计算出波达方向角，然后将结果传输至输出模块。

接收模块包括2个按二级嵌套阵进行排列的超声波探头阵列，步进电机和步进电机驱动电路。步进电机是将电脉冲信号转变为角位移或线位移的开环控制电机，当步进电机驱动电路收到一个脉冲信号，它就驱动步进电机按设定的方向转动固定的角度，称为步距角。所以可以通过使数据处理与控制模块发射一定数量的脉冲信号来达到期望的角度值。如图3所示，嵌套阵1和嵌套阵2是结构相似的二级嵌套阵，每条二级嵌套阵上分布有(M+N)个阵元，每条二级嵌套阵由两级均匀线阵组成，第一级均匀线阵有M个阵元，阵元间距为d，第二级均匀线阵有N个阵元，阵元间距为(M+1)d，两条嵌套阵连接处有一个公共阵元。使嵌套阵1排布于坐标系x轴上保持固定，因为接收模块会放置在水中，所以固定支架采用塑料材质以增大浮力。嵌套阵2安装到步进电机上，可由步进电机带动旋转，从而达到夹角调节的目的，图6为嵌套阵2与步进电机的连接旋转示意图，如图所示，步进电机通过旋转转子与嵌套阵2连接在一起，以控制嵌套阵2的旋转。

发射模块由一个阻抗匹配电路和一个超声波发射探头组成，通过D/A转换器与处理器相连，能够根据处理器发出的指令发射指定的信号。

输出模块由一个USB接口和一个显示器组成，并且与数据处理与控制模块和电源模块相连。它能够提供人机交互，将数据处理与控制模块中处理好的数据通过USB接口输出到外部装置或者在显示器上显示出来。

电源模块由一个电源组成，并且与数据处理与控制模块、发射模块、接收模块和输出模块相连。它能够为这些模块供电。

本装置的主要工作流程如下：在实测过程中根据用户想要发射的信号参数，通过数据处理与控制模块输入对应的参数，使处理器产生相应的数字信号，然后通过D/A转换后传给发射模块，超声波发射探头就能产生我们需要的信号并进行发射。嵌套阵1和嵌套阵2之间的夹角值δ可以通过数据处理与控制模块进行设定，处理器发送特定的脉冲信号到步进电机驱动电路，然后驱动步进电机转动至我们需要的角度。接收模块中的接收阵列收到从目标声源反射回来的信号后将其通过A/D转换成数字信号后发送给处理器，然后处理器根据本发明提供的算法计算出结果。最后数据处理与控制模块将计算结果传给输出模块，输出模块将结果通过USB接口传给外部设备或者通过显示器显示出来。电源模块为所有其它模块供电。

实施例三

本实施例公开了一种水下高精度一维DOA估计装置，所述的估计装置包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块，如图1和图2所示。

数据处理与控制模块可以用DSP芯片实现(如：TI公司TMS320VC5509A型号的DSP芯片)，此DSP芯片可实现A/D转换和D/A转换的功能，并能够实现二维嵌套阵矩阵重构和最终波达方向的计算；

接收模块中的步进电机采用东芝公司的23HY6606-CP型号的电机，此步进电机的步距角为1.8度，步进电机驱动电路采用东芝公司的TC78S600FTG型芯片。此外接收模块还使用1个固定嵌套阵和1个自由旋转的嵌套阵，其中每个阵列包括多个超声接收探头，并且数量相同，2个嵌套阵按图5所示组装；发射模块使用一个超声波发射探头；输出模块使用一个USB接口和一个LCD显示屏。图1即为本发明所述装置的硬件结构模块图。

本发明的主要工作步骤具体如下：

步骤T1、按图2连接好具体装置，其中接收模块中每个嵌套阵中设置M为3，N为4，单条嵌套阵一共有7个阵元。利用数据处理与控制模块发送指令，控制超声发射探头发射超声信号s(t)；海水中声速范围大致为1430m/s-1550m/s，则取最小声速为1430m/s，可以求出最小半波长为7.15cm。考虑到DOA估计中最小阵元间距应该小于半波长，设置每条嵌套阵的第一级均匀线阵的阵元间距为5cm，第二级均匀线阵的阵元间距为20cm，设置嵌套阵1和嵌套阵2之间的夹角为150°。在数据处理与控制模块设定线阵夹角值，并控制步进电机使嵌套阵2转到设定的位置。在水下放置K＝3个目标信源，信源的中心频率均为15kHz，入射的波达方向角分别为(40°，50°，60°)，假设实时声速为c，其中c∈(1430m/s，1550m/s)。

步骤T2、对超声接收探头线阵接收到的目标声源信号进行采样；嵌套阵1接收到的信号为x₁(t)，x₂(t)，...，x₇(t)，嵌套阵2接收到的信号为y₁(t)，y₂(t)，...，y₇(t)，采样200次，将采样后的数据传递给数据处理与控制模块进行分析处理。

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

1)将两条嵌套阵200次采样后的数据以矢量形式表示为X，Y，分别计算两者的协方差矩阵

2)对R_xx和R_yy向量化得到z_xx和z_yy，对z_xx和z_yy进行去冗余并排序，得到

和/>

使用矩阵重构算法基于/>

和/>

构建Toeplitz矩阵/>

和/>

/>

3)基于

使用TLS-ESPRIT算法得到/>

对/>

进行特征值分解得到特征值λ_x1，λ_x2，...，λ_xK，基于/>

使用TLS-ESPRIT算法得到/>

对/>

进行特征值分解得到特征值λ_y1，λ_y2，...，λ_yK，将λ_x1，λ_x2，...，λ_xK从小到大排序，对λ_y1，λ_y2，...，λ_yK从大到小排序，通过排序对两者进行配对；

4)根据交叉嵌套阵的几何关系推导出θ_xk、θ_yk和θ_k的关系式，结合公式公式(16)最终求解出目标信号的波达方向角。

根据本发明中水下高精度一维DOA估计方法，实时声速为1480m/s时估计出的波达方向角为(40.15°，50.11°，59.89°)，仅使用13个阵元得到了声速无关的高精度的估计结果，验证了本发明方法及装置的可行性。

综上所述，上述实施例提出一种水下高精度一维DOA估计方法与装置，基于两条嵌套阵接收数据得到两组角度估计值，通过排序对两组估计值进行配对，最终得到声速无关的一维DOA估计结果。所提方法的配对复杂度很小，估计精度较高。

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

Claims (4)

1.一种水下高精度一维DOA估计方法，其特征在于，所述估计方法包括以下步骤：

S1、采用二维交叉嵌套阵进行一维DOA估计，该二维交叉嵌套阵采用两条二级嵌套阵，每条二级嵌套阵上各自分布有(M+N)个阵元，每条二级嵌套阵由两级均匀线阵组成，第一级均匀线阵有M个阵元，阵元间距为d，第二级均匀线阵有N个阵元，阵元间距为(M+1)d，两条二级嵌套阵的连接处有一个公共阵元，两条二级嵌套阵之间的夹角设为δ，其中

建立坐标系，设一条二级嵌套阵所在直线为x轴，另一条二级嵌套阵所在直线设为y轴，假设只考虑半个平面空间内目标源信号的情况，即x轴的上半平面空间，假设目标信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移，假设目标信号的个数为K，目标信号的中心频率为f，实时声速定义为c，目标信号在x轴上的波达方向角为θ_xk,k＝1,2,…,K，在y轴上的波达方向角为θ_yk，目标信号的波达方向角设为θ_k，且θ_k＝θ_xk，当有K个远场窄带相互独立的信号入射到上述二维交叉嵌套阵，x轴和y轴阵列接收到的信号写成如下的矢量形式：

其中s(t)是一个K×1维的源信号矩阵，n_x(t)和n_y(t)分别为x轴和y轴阵列的噪声矩阵，A_x和A_y分别为x轴和y轴阵列的方向矩阵，写成矢量形式有：

其中a(θ_xk)和a(θ_yk)分别是x轴和y轴阵列的第k个声源的导向向量，

对于两条线阵各自接收的信号x(t)和y(t)，分别求它们的协方差矩阵为R_xx＝E[x(t)x(t)^H]和R_yy＝E[y(t)y(t)^H]；

S2、向量化R_xx和R_yy得到z_xx和z_yy，对z_xx和z_yy进行去冗余并排序，得到

和/>

使用矩阵重构算法基于/>

和/>

构建Toeplitz矩阵/>

和/>

S3、基于

使用TLS-ESPRIT算法得到/>

对/>

进行特征值分解得到特征值λ_x1，λ_x2，...，λ_xK，基于/>

使用TLS-ESPRIT算法得到/>

对/>

进行特征值分解得到特征值λ_y1，λ_y2，...，λ_yK，将λ_x1，λ_x2，...，λ_xK从小到大排序，对λ_y1，λ_y2，...，λ_yK从大到小排序，通过排序对两者进行配对；

S4、根据二维交叉嵌套阵的几何关系推导出θ_xk、θ_yk和θ_k的关系式，最终求解出目标信号的波达方向角，过程如下：

根据二维嵌套阵的阵列结构，得到目标信号相对于x轴阵列的波达方向角θ_xk、相对于y轴阵列的波达方向角θ_yk和待估计的波达方向角θ_k满足：

cosθk＝cosθ_xk＝cos(δ-θ_yk) 公式(14)

然后根据公式(13)，得到目标信号相对于x轴阵列的波达方向角θ_xk为

考虑到第k个信号的波达方向角θ_k就是其在x轴上的波达方向角θ_xk，因此最终的一维DOA估计结果为：

2.根据权利要求1所述的一种水下高精度一维DOA估计方法，其特征在于，所述步骤S1中对接收信号矩阵进行处理得到协方差矩阵R_xx和R_yy的过程如下：

由于噪声和信号相互独立，得到x轴和y轴阵列接收数据的M×M阶协方差矩阵分别为：

其中R_s＝E[s(t)s^H(t)]为信源部分的协方差矩阵，由于信源相互独立，因此R_s为对角矩阵，I为单位矩阵，σ_n ²为噪声功率，()^H表示求共轭转置。

3.根据权利要求2所述的一种水下高精度一维DOA估计方法，其特征在于，所述步骤S2中得到Toeplitz矩阵

和/>

的过程如下：

向量化R_xx和R_yy，得到向量z_xx和z_yy：

其中p为向量化后的信号矩阵，I_n为向量化后的单位矩阵，⊙表示求Khatri-Rao积，()^*表示求共轭；

对向量z_xx和z_yy进行去冗余并排序，得到新的向量

其中

向量化、排序、去冗余后的

和/>

的信号矩阵p的秩为1，根据矩阵重构算法得到

其中

表示/>

的第m个元素，其中/>

同理可以得到

4.根据权利要求3所述的一种水下高精度一维DOA估计方法，其特征在于，所述步骤S3中对λ_x1，λ_x2，...，λ_xK和λ_y1，λ_y2，...，λ_yk进行配对的过程如下：

基于Toeplitz矩阵

和/>

根据现有的TLS-ESPRIT算法得到/>

和/>

形式如下

其中Φ_x和T_x ^-1分别为

的特征值组成的对角矩阵以及对应的特征向量组成的满秩方阵，Φ_y和T_y ^-1分别为/>

的特征值组成的对角矩阵以及对应的特征向量组成的满秩方阵，Φ_x和Φ_y的具体形式为

对

进行特征分解，得到特征值λ_x1，λ_x2，...，λ_xK，对/>

进行特征分解，得到特征值λ_y1，λ_y2，...，λ_yK；

λ_x1，λ_x2，...，λ_xK、λ_y1，λ_y2，...，λ_yK与θ_xk、θ_yk的关系为

此时λ_x1，λ_x2，...，λ_xK、λ_y1，λ_y2，...，λ_yK是没有配对的，考虑二维交叉嵌套阵的阵列结构，通过排序进行配对，对λ_x1，λ_x2，...，λ_xK按照从大到小的顺序排序，对λ_y1，λ_y2，...，λ_yK按照从小到大的顺序排序，排序之后两组特征值的顺序是配对的。

Images