CN110109053B - 一种未知声速环境下快速doa估计方法 - Google Patents

Abstract

本发明公开了一种未知声速环境下快速DOA估计方法，步骤如下：采用任意交叉线阵结构进行一维DOA估计；对互协方差矩阵R进行奇异值分解，其奇异值从大到小排列，前K个奇异值对应的左奇异向量张成的空间为信号子空间，这K个奇异向量组成矩阵U_s，将U_s划分为矩阵U₁和矩阵U₂；求出目标信号在x轴上的波达方向角θ_xi以及在y轴上的波达方向角θ_yi的表达式；将目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式；根据两条线阵的几何关系推导出θ_xi和θ_yi的关系式，再最终求解出目标信号的波达方向角。该发明在得到声速无关的DOA估计表达式的同时解决了DOA估计处理过程中的参数配对问题，减小了DOA估计处理的复杂度。

Description

一种未知声速环境下快速DOA估计方法

技术领域

本发明涉及目标定位技术领域，具体涉及一种基于任意夹角的二维均匀线阵在未知声速环境进行的快速水下一维波达方向估计的方法。

背景技术

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

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

发明内容

本发明的目的是为了解决现有技术中的上述缺陷，提供一种未知声速环境下快速DOA估计方法，通过对任意夹角的二维均匀线阵的接收信号进行处理，在DOA波达方向估计中消除声速这个因子，从而消除水下声速不确定性对目标定位精度的影响。同时由于采用了快速配对方法，算法复杂度得到很大的下降，有利于在实际测量中快速估计目标。同时在估计过程中采用双线阵接收信号的互协方差矩阵进行处理，本发明无需对多条线阵的接收信号进行独立估计再进行配对，有效减小了算法复杂度。

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

一种基于未知声速环境进行的快速水下一维波达方向估计的方法，所述的估计方法步骤如下：

S1、用于一维DOA估计的任意交叉线阵结构如图3所示，与常规一维DOA估计方法采用单条均匀线阵不同，本文所提的方法采用两条均匀线阵，每条线阵上各自分布有M个阵元，两条线阵接口处有一个公共阵元，阵元间距为d。两条线阵之间的夹角设为δ且

建立坐标系，设一条线阵所在直线为_x轴，另一条线阵所在直线设为y轴。由于水下DOA估计大部分情况是在水面位置对水面下方的目标进行定位，因此本文只考虑半个平面空间内目标源信号的情况，即_x轴的上半平面空间。目标信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。假设目标信号的个数为K，目标信号的中心频率为f_i，i＝1,2,...,K，各个信号入射路径的声速定义为c_i,i＝1,2,...,K，目标信号在_x轴上的波达方向角设为θ_xi,i＝1,2,...,K，在y轴上的波达方向角设为θ_yi,i＝1,2,...,K。目标信号的待估计的波达方向角设为θ_i且θ_i＝θ_xi,i＝1,2,...,K。第i个信号与x轴线阵的夹角为α_i且α_i∈[0,π]。当有K个远场窄带相互独立的信号入射到图3所示阵列，_x轴和y轴阵列接收到的信号可以写成如下的矢量形式：

X(t)＝A_xS+N_x (1)

Y(t)＝A_yS+N_y (2)

其中S是一个K×1维的源信号矩阵，另外N_x和N_y则是M×1维的噪声矩阵，A_x和A_y分别为_x轴和y轴阵列的M×K阶方向矩阵，写成矢量形式有：

其中a(θ_xi)和a(θ_yi)分别是_x轴和y轴阵列的第i个声源的导向向量，i＝1,2,...,K。

对于两条线阵各自接收的信号X(t)和Y(t)，求它们的互协方差矩阵R_xy＝E{X(t)Y^H(t)}，然后将互协方差矩阵R_xy划分为子矩阵R_xy,1和R_xy,2，最后将子矩阵R_xy,1和R_xy,2组合成新的互协方差矩阵R。

S2、对矩阵R进行奇异值分解，其奇异值从大到小排列，前K个奇异值对应的左奇异向量张成的空间为信号子空间，这K个奇异向量组成矩阵U_s。将U_s划分为矩阵U₁和矩阵U₂。

S3、对矩阵

进行特征值分解，得到K个特征值λ₁,λ₂,...,λ_K和对应的特征向量α₁,α₂,...,α_K，并由这K个特征向量组成的矩阵估计出非奇异阵T，根据K个特征值求出目标信号在_x轴上的波达方向角θ_xi的表达式，以及求出目标信号在y轴上的波达方向角θ_yi的表达式。

S4、将步骤S3得到的目标信号在y轴上的波达方向角表达式和在_x轴上的波达方向角表达式相除得到一个与声速无关的表达式(33)。

S5、根据两阵列的几何关系推导出θ_xi和θ_yi的关系式，再最终求解出目标信号的波达方向角。

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

根据线阵夹角δ和α_i将目标源信号入射区域分成4部分：当入射信号与x轴的夹角α_i∈[0,δ)时，信号入射区域设为区域①，

时为区域②，

时为区域③，

时为区域④。

设入射信号与x轴法线的夹角为x_ni且

与y轴法线的夹角为y_ni且

以原点处的阵元为参考阵元，对于_x轴阵列，当入射信号最先到达参考阵元时，则波达方向角θ_xi为正值，此时波达方向角等于信号与阵列法线的夹角，即θ_xi＝x_ni。否则当信号最后到达参考阵元时，波达方向角θ_xi为负值，且θ_xi＝-x_ni。对于y轴阵列也有相同的结论。因此有：

sinθ_yi＝sin(θ_xi-δ) (5)

得到X(t)和Y(t)的M×M阶互协方差矩阵为：

由于噪声为零均值白噪声，各噪声之间相互统计独立且都与目标信号统计独立，所以式(6)的后三项都为零，公式(6)改写为：

其中R_s＝E{S(t)S^H(t)}为信源部分的互协方差矩阵，由各信源是统计独立的，且R_s是对角矩阵。

取A_y的前(M-1)行和后(M-1)行分别设为A_y1和A_y2，即：

由A_y的表达式推出：

A_y2＝A_y1Ω^H (9)

其中

接下来将矩阵R_xy的前M-1列划分为子阵R_xy,1，后M-1列划分为子阵R_xy,2，即：

结合式(9)进一步可以得出：

将R_xy,1和R_xy,2按下式组合成新的2M×(M-1)阶矩阵，将此矩阵定义为新的互协方差矩阵R，

进一步地，所述的步骤S2中对矩阵R进行奇异值分解和排列，得到矩阵U₁和矩阵U₂的过程如下：

根据式(12)定义新的方向矩阵为：

则式(12)可以重写为：

对R进行奇异值分解：

其中U和V分别表示R的左右奇异矢量，且它们都是酋矩阵。Σ为对角矩阵，其对角线元素表示奇异值且它们从大到小排列，在无噪声的情况下，前K个奇异值大于零，剩余的奇异值为零。在有噪声影响的情况下，前K个奇异值也远远大于剩余奇异值。因此可以用前K个奇异值对应的奇异矢量构成信号子空间U_s，剩余(M-1-K)个奇异值对应的奇异矢量构成噪声子空间U_n。Σ_s和Σ_n对角线上的元素分别表示信号子空间和噪声子空间的奇异值，V_s和V_n分别为信号子空间和噪声子空间的右奇异矢量。

由奇异值分解性质可知：

RV_n＝Σ_nU_n (16)

两边都取共轭转置并右乘U_n得到：

由式(15)和式(17)可知，由于噪声之间的独立性，矩阵R不受噪声的影响，因此噪声子空间的奇异值都为0，即Σ_n＝0。因此式(17)改写为：

V_n ^HR^HU_n＝0 (18)

又V_n为矩阵R奇异值分解得到的奇异值向量矩阵，所以V_n是满秩矩阵，由满秩矩阵性质可解出式(18)得：R^HU_n＝0，并结合式(14)代入式(18)可得：

由于

是满秩的，因此：

B^HU_n＝0 (20)

由于U＝[U_s,U_n]为酋矩阵，所以有：

结合式(20)和(21)可知，方向矩阵B与矩阵R的信号特征向量组成的子矩阵U_s所张成的子空间相同，因此存在一个非奇异矩阵T使得：

U_s＝BT (22)

由式(13)可以得知方向矩阵B为2M×K阶矩阵，由两个M×K的矩阵A_x和A_xΩ组成，易知A_x具有范德蒙德结构。定义两个子阵B₁和B₂，B₁由A_x的前M-1行和A_xΩ的前M-1行组成，B₂由A_x的后M-1行和A_xΩ的后M-1行组成，即：

观察式(13)和(23)很显然可以得到：

B₂＝B₁Ψ (24)

其中Ψ＝diag{ψ₁,ψ₂,...,ψ_K}且

Ψ叫做旋转矩阵其对角线元素为相位旋转算子，易知通过求解旋转算子ψ_i就可以求出θ_xi。将U_s按照相同的方法划分为U₁和U₂，即：

由Us＝BT可得到：

进一步地，所述的步骤S3中对矩阵

进行特征值分解，进行处理后求出目标信号进行处理后求出目标信号在在x轴上的波达方向角θ_xi的表达式和y轴上的波达方向角θ_yi的表达式的过程如下：

由式(24)和(26)可以推出：

U₂＝B₁ΨT＝B₁TT^-1ΨT＝U₁T^-1ΨT (27)

其中T^-1为T的逆矩阵，定义

为矩阵U₁的Moore-Penrose广义逆，公式(27)两边左乘

有：

则矩阵

为矩阵Ψ的相似变换，所以它们有相同的特征值，而Ψ由于是对角矩阵因此其特征值就是其对角线元素

所以可以求解出

的K个特征值λ₁,λ₂,...,λ_K和对应的特征向量α₁,α₂,...,α_K。最后可以求出θ_xi的表达式为：

观察公式(28)可知，矩阵T^-1由

的K个特征向量α₁,α₂,...,α_K组成，则可得到估计值

由(22)可以得出方向矩阵B的估计值为：

观察B的展开式(11)，可以得到如下表达式：

其中

表示矩阵

的第k行第i列的元素。

进一步地，所述的步骤S4中对得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式的过程如下：

此时求出的结果还跟声速和频率有关，还需进一步消除声速和频率的影响。

用公式(32)和(29)相除可以得到：

进一步地，所述的步骤S5中根据两阵列的几何关系推导出θ_xi和θ_yi的关系式，最终求解出目标信号的波达方向角的过程如下：

由于只对矩阵

进行了一次特征值分解，由特征值得到上式的分母，特征值对应的特征向量矩阵得到上式的分子，因此无需进行参数配对。再结合由两线阵之间的几何关系得出的公式(5)可以求出θ_xi的表达式如下：

根据第i个信号的波达方向角θ_i就是其在x轴上的波达方向角θ_xi，得到最终的一维波达方向估计结果为：

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

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

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

3、本发明方法的时间复杂度主要集中在3个部分：一次互协方差矩阵构造、一次奇异值分解和一次特征值分解，其它部分的计算复杂度可以忽略。设采样快拍数为L，则两个线阵的接收信号X(t)和Y(t)的阶数都是M×L，因此构造互协方差矩阵的复杂度为O(M ²L)。矩阵R的阶数为2M×(M-1)，因此对它进行奇异值分解的复杂度为O((M-1)³)。矩阵

的阶数为K×K，对其进行特征值分解的复杂度为O(K ³)，因此总的复杂度为

O(M ²L+(M-1)³+K³)；

相较于其他未知声速环境的DOA估计方法，本发明提出的DOA估计方法时间复杂度更低，而且当目标信源数越多时，本发明所提方法的优势越明显，有利于目标的实时DOA估计。

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

附图说明

图1是本发明装置的硬件结构模块图；

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

图3是本发明所用的任意夹角二维均匀线阵模型示意图；

图4是子线阵3的旋转连接示意图；

图5是x轴均匀线阵的接收信号模型示意图；

图6是本发明公开的未知声速环境下快速DOA估计方法的流程图。

具体实施方式

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

实施例一

本发明采用任意夹角的二维线阵，窄带目标声源为S，中心频率为f。声波入射方向，即方位角和仰角可表示为θ_k和

如附图6所示，本实施例中基于水下DOA估计方法包括以下步骤：

S1、用于一维DOA估计的任意交叉线阵结构如图3所示，与常规一维DOA估计方法采用单条均匀线阵不同，本文所提的方法采用两条均匀线阵，每条线阵上各自分布有M个阵元，两线阵接口处有一个公共阵元，阵元间距为d。两条线阵之间的夹角设为δ且

建立坐标系，设一条线阵所在直线为x轴，另一条线阵所在直线设为y轴。由于水下DOA估计大部分情况是在水面位置对水面下方的目标进行定位，因此本文只考虑半个平面空间内目标源信号的情况，即x轴的上半平面空间。目标信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。当有K个远场窄带相互独立的信号入射到图3所示阵列，x轴和y轴阵列接收到的信号可以写成如下的矢量形式：

X(t)＝A_xS+N_x (1)

Y(t)＝A_yS+N_y (2)

其中S是一个K×1维的源信号矩阵，另外N_x和N_y则是M×1维的噪声矩阵，A_x和A_y分别为x轴和y轴阵列的M×K阶方向矩阵，写成矢量形式有：

其中a(θ_xi)和a(θ_yi)分别是x轴和y轴阵列的第i个声源的导向向量，i＝1,2,...,K。

假设目标信号的个数为K，目标信号的中心频率为f_i，i＝1,2,...,K，各个信号入射路径的声速定义为c_i,i＝1,2,...,K，目标信号在x轴上的波达方向角设为θ_xi,i＝1,2,...,K，在y轴上的波达方向角设为θ_yi,i＝1,2,...,K。目标信号的待估计的波达方向角设为θ_i且θ_i＝θ_xi,i＝1,2,...,K。第i个信号与x轴线阵的夹角为α_i且α_i∈[0,π]。根据线阵夹角δ和α_i将目标源信号入射区域分成4部分：当入射信号与x轴的夹角α_i∈[0,δ)时，信号入射区域设为区域①，

时为区域②，

时为区域③，

时为区域④。

设入射信号与x轴法线的夹角为x_ni且

与y轴法线的夹角为y_ni且

当信号从区域①入射时，有

则有：

以原点处的阵元为参考阵元，对于x轴阵列，当入射信号最先到达参考阵元时，则波达方向角θ_xi为正值，此时波达方向角等于信号与阵列法线的夹角，即θ_xi＝x_ni。否则当信号最后到达参考阵元时，波达方向角θ_xi为负值，且θ_xi＝-x_ni。对于y轴阵列也有相同的结论。因此有：

sinθ_yi＝sin(θ_xi-δ) (5)

X(t)和Y(t)的M×M阶互协方差矩阵为：

由于噪声为零均值高斯白噪声，各噪声之间相互统计独立且都与目标信号统计独立所以式(6)的后三项都为零，公式(6)改写为：

其中R_s＝E{S(t)S^H(t)}为信源部分的互协方差矩阵，由各信源是统计独立的，且R_s是对角矩阵。

取A_y的前(M-1)行和后(M-1)行分别设为A_y1和A_y2，即：

由A_y的表达式推出：

A_y2＝A_y1Ω^H (9)

其中

接下来将矩阵R_xy的前M-1列划分为子阵R_xy,1，后M-1列划分为子阵R_xy,2，即：

结合式(9)进一步可以得出：

将R_xy,1和R_xy,2按下式组合成新的2M×(M-1)阶矩阵，将此矩阵定义为新的互协方差矩阵R，

S2、对互协方差矩阵R进行奇异值分解，其奇异值从大到小排列，前K个奇异值对应的左奇异向量张成的空间为信号子空间，这K个奇异向量组成矩阵U_s。将U_s按照式(36)的方式划分为矩阵U₁和矩阵U₂。根据式(12)定义新的方向矩阵为：

则式(12)可以重写为：

对互协方差矩阵R进行奇异值分解：

其中U和V分别表示R的左右奇异矢量，且它们都是酋矩阵。Σ为对角矩阵，其对角线元素表示奇异值且它们从大到小排列，在无噪声的情况下，前K个奇异值大于零，剩余的奇异值为零。在有噪声影响的情况下，前K个奇异值也远远大于剩余奇异值。因此可以用前K个奇异值对应的奇异矢量构成信号子空间U _s，剩余(M-1-K)个奇异值对应的奇异矢量构成噪声子空间U_n。Σ_s和Σ_n对角线上的元素分别表示信号子空间和噪声子空间的奇异值，V_s和V_n分别为信号子空间和噪声子空间的右奇异矢量。

由奇异值分解性质可知：

RV_n＝Σ_nU_n (16)

两边都取共轭转置并右乘U_n得到：

由式(6)和式(7)可知，由于噪声之间的独立性，矩阵R不受噪声的影响，因此噪声子空间的奇异值都为0，即Σ_n＝0。因此式(17)改写为：

V_n ^HR^HU_n＝0 (18)

又V_n为互协方差矩阵R奇异值分解得到的奇异值向量矩阵，所以V_n是满秩矩阵，由满秩矩阵性质可解出式(18)得：

R^HU_n＝0 (19)

结合式(13)代入式(19)可得：

由于

是满秩的，因此B^HU_n＝0。由于U＝[U_s,U_n]为酋矩阵，所以有：

结合式(20)和(21)可知，方向矩阵B与矩阵R的信号特征向量组成的子矩阵Us所张成的子空间相同，因此存在一个非奇异矩阵T使得：

U_s＝BT (22)

由式(13)可以得知方向矩阵B为2M×K阶矩阵，由两个M×K的矩阵A_x和A_xΩ组成，易知A_x具有范德蒙德结构。定义两个子阵B₁和B₂，B₁由A_x的前M-1行和A_xΩ的前M-1行组成，B₂由A_x的后M-1行和A_xΩ的后M-1行组成，即：

观察式(22)和(23)很显然可以得到：

B₂＝B₁Ψ (24)

其中Ψ＝diag{ψ₁,ψ₂,...,ψ_K}且

Ψ叫做旋转矩阵其对角线元素为相位旋转算子，易知通过求解旋转算子ψ_i就可以求出θ_xi。将U_s按照相同的方法划分为U₁和U₂，即：

由U_s＝BT可得到：

S3、对矩阵

进行特征值分解，得到K个特征值λ₁,λ₂,...,λ_K和对应的特征向量α₁,α₂,...,α_K，并由这K个特征向量组成的矩阵估计出非奇异阵T，根据K个特征值求出目标信号在x轴上的波达方向角θ_xi的表达式，求出目标信号在y轴上的波达方向角θ_yi的表达式。

由式(24)和(26)可以推出：

U₂＝B₁ΨT＝B₁TT^-1ΨT＝U₁T^-1ΨT (27)

定义

为矩阵U₁的Moore-Penrose广义逆，公式(27)两边左乘

有：

则矩阵

为矩阵Ψ的相似变换，所以它们有相同的特征值，而Ψ由于是对角矩阵因此其特征值就是其对角线元素

所以可以求解出

的K个特征值λ₁,λ₂,...,λ_K和对应的特征向量α₁,α₂,...,α_K。最后可以求出θ_xi的表达式为：

此时求出的结果还跟声速和频率有关，还需进一步消除声速和频率的影响。观察公式(28)易知矩阵T^-1由

的K个特征向量α₁,α₂,...,α_K组成，则得到估计值

由(22)可以得出方向矩阵B的估计值为：

观察B的展开式(13)，可以得到如下表达式：

其中

表示矩阵

的第k行第i列的元素。

S4、将步骤S3得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相比得到一个与声速无关的表达式。用公式(32)和(29)相除可以得到：

S5、根据两阵列的几何关系推导出θ_xi和θ_yi的关系式，再结合式(33)最终求解出目标信号的波达方向角。由于只对矩阵

进行了一次特征值分解，由特征值得到上式的分母，特征值对应的特征向量矩阵得到上式的分子，因此无需进行参数配对。再结合由两线阵之间的几何关系得出的公式(5)可以求出θ_xi的表达式如下：

因为第i个信号的波达方向角θ_i就是其在x轴上的波达方向角θ_xi，所以最终的一维波达方向估计结果为：

实施例二

本实例公开了一种基于未知声速环境的快速水下波达方向估计装置，所述的估计装置包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块，如图1和图2所示。

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

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

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

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

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

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

实施例三

本实例公开了一种基于未知声速环境的快速水下波达方向估计装置，所述的估计装置包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块，如图1和图2所示。

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

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

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

步骤T1、按图2连接好具体装置，其中接收模块中的每个均匀线阵中的阵元个数M统一定为8。利用数据处理与控制模块发送指令，控制超声发射探头发射超声信号s(t)；海水中声速范围大致为1430m/s-1550m/s，则取最小声速为1430m/s，可以求出最小半波长为7.15cm。所以设置两个均匀线阵的平均间距为5cm，即第一个阵元和最后一个阵元相隔35cm。任意两相邻线阵之间的距离必须小于7.15cm，在满足此限制条件下可以任意选取阵元间距，这里2个均匀线阵的间距都取为4cm。线阵1与线阵2之间的夹角

设置线阵夹角值为72°。在数据处理与控制模块设定线阵夹角值，首先将均匀线阵夹角δ转为72°。在水下放置K＝3个目标声源，目标信号的中心频率分别为f＝10,15,20kHz，各个信号入射路径的声速分别为c＝1450,1500,1550m/s，入射的波达方向角分别为(40°,60°，70°)。

步骤T2、对超声接收探头线阵接收到的目标声源信号进行采样；均匀线阵1接收到的信号为x₁(t),x₂(t),…,x₈(t)，均匀线阵2接收的信号为y1(t),y2(t),…,y₈(t)，共采样接收200次，并将接收到的信号传递给数据处理与控制模块进行分析处理。

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

1)对于两条线阵各自接收的信号X(t)和Y(t)，求它们的互协方差矩阵R_xy＝E{X(t)Y^H(t)}，然后将互协方差矩阵R_xy划分为子矩阵R_xy,1和R_xy,2，最后将R_xy,1和R_xy,2按组合成新的矩阵R。

2)对矩阵R进行奇异值分解，其奇异值从大到小排列，前K个奇异值对应的左奇异向量张成的空间为信号子空间，这K个奇异向量组成矩阵U_s。将U_s划分为矩阵U₁和矩阵U₂。

3)对矩阵

进行特征值分解，得到K个特征值λ₁,λ₂,...,λ_K和对应的特征向量α₁,α₂,...,α_K，并由这K个特征向量组成的矩阵估计出非奇异阵T，根据K个特征值按式(29)求出目标信号在x轴上的波达方向角θ_xi的表达式，按式(32)求出目标信号在y轴上的波达方向角θ_yi的表达式。

4)将以上步骤3)得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相比得到一个与声速无关的表达式(33)。

5)根据两阵列的几何关系推导出θ_xi和θ_yi的关系式，再结合式(33)最终求解出目标信号的波达方向角。

根据本发明中快速DOA估计方法，估计出的二维波达方向角(40.15°,59.89°,70.23°)，对目标估计达到了预期精度，说明估计结果正确，本发明方法及装置可行。

综上所述，上述实施例提出一种基于任意交叉线阵的无需配对的声速无关快速一维DOA估计方法与装置，利用对交叉线阵的接收信号进行联合处理，在得到声速无关的DOA估计表达式的同时解决了DOA估计处理过程中的参数配对问题，减小了DOA估计处理的复杂度。

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

Claims (6)

1.一种未知声速环境下快速DOA估计方法，其特征在于，所述的估计方法包括以下步骤：

S1、采用任意交叉线阵结构进行一维DOA估计，该任意交叉线阵结构采用两条均匀线阵，每条线阵上各自分布有M个阵元，两条线阵接口处有一个公共阵元，阵元间距为d，两条线阵之间的夹角设为δ且

建立坐标系，设一条线阵所在直线为x轴，另一条线阵所在直线设为y轴，假设只考虑半个平面空间内目标源信号的情况，即x轴的上半平面空间，假设目标信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。假设目标信号的个数为K，目标信号的中心频率为f_i,i＝1,2,...,K，各个信号入射路径的声速定义为c_i,i＝1,2,...,K，目标信号在x轴上的波达方向角设为θ_xi,i＝1,2,...,K，在y轴上的波达方向角设为θ_yi,i＝1,2,...,K，目标信号的待估计的波达方向角设为θ_i且θ_i＝θ_xi,i＝1,2,...,K，第i个信号与x轴线阵的夹角为α_i且α_i∈[0,π]，当有K个远场窄带相互独立的信号入射到上述任意交叉线阵结构，x轴和y轴阵列接收到的信号写成如下的矢量形式：

X(t)＝A_xS+N_x (1)

Y(t)＝A_yS+N_y (2)

其中S是一个K×1维的源信号矩阵，N_x和N_y是M×1维的噪声矩阵，A_x和A_y分别为x轴和y轴阵列的M×K阶方向矩阵，写成矢量形式有：

其中a(θ_xi)和a(θ_yi)分别是x轴和y轴阵列的第i个声源的导向向量，i＝1,2,...,K，

对于两条线阵各自接收的信号X(t)和Y(t)，求它们的互协方差矩阵R_xy＝E{X(t)Y^H(t)}，然后将互协方差矩阵R_xy划分为子矩阵R_xy,1和R_xy,2，最后将子矩阵R_xy,1和R_xy,2组合成新的互协方差矩阵R；

S2、对互协方差矩阵R进行奇异值分解，其奇异值从大到小排列，前K个奇异值对应的左奇异向量张成的空间为信号子空间，这K个奇异向量组成矩阵U_s，将U_s划分为矩阵U₁和矩阵U₂；

S3、对矩阵U₁ ⁺U₂进行特征值分解，得到K个特征值λ₁,λ₂,...,λ_K和对应的特征向量α₁,α₂,...,α_K，并由这K个特征向量组成的矩阵估计出非奇异阵T，根据K个特征值求出目标信号在x轴上的波达方向角θ_xi的表达式，以及求出目标信号在y轴上的波达方向角θ_yi的表达式；

S4、将步骤S3得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式；

S5、根据两条线阵的几何关系推导出θ_xi和θ_yi的关系式，再最终求解出目标信号的波达方向角。

2.根据权利要求1所述的一种未知声速环境下快速DOA估计方法，其特征在于，所述的步骤S1中对接收信号矩阵进行处理得到新的互协方差矩阵R的过程如下：

根据线阵夹角δ和α_i将目标源信号入射区域分成4部分：当入射信号与x轴的夹角α_i∈[0,δ)时，信号入射区域设为区域①，

时为区域②，

时为区域③，

时为区域④；

设入射信号与x轴法线的夹角为x_ni且

与y轴法线的夹角为y_ni且

以原点处的阵元为参考阵元，对于x轴阵列，当入射信号最先到达参考阵元时，则波达方向角θ_xi为正值，此时波达方向角等于信号与阵列法线的夹角，即θ_xi＝x_ni，否则当信号最后到达参考阵元时，波达方向角θ_xi为负值，且θ_xi＝-x_ni，对于y轴阵列也有相同的结论，因此有：

sinθ_yi＝sin(θ_xi-δ) (5)

得到X(t)和Y(t)的M×M阶互协方差矩阵为：

由于噪声为零均值白噪声，各噪声之间相互统计独立且都与目标信号统计独立，将式(6)的后三项都为零，式(6)改写为：

其中R_s＝E{S(t)S^H(t)}为信源部分的互协方差矩阵，由各信源是统计独立的，且R_s是对角矩阵；

取A_y的前(M-1)行和后(M-1)行分别设为A_y1和A_y2，即：

由A_y的表达式推出：

A_y2＝A_y1Ω^H (9)

其中

接下来将矩阵R_xy的前M-1列划分为子阵R_xy,1，后M-1列划分为子阵R_xy,2，即：

结合式(9)得出：

将R_xy,1和R_xy,2按下式组合成新的2M×(M-1)阶矩阵，将此矩阵定义为新的互协方差矩阵R，

3.根据权利要求2所述的一种未知声速环境下快速DOA估计方法，其特征在于，所述的步骤S2中对互协方差矩阵R进行奇异值分解和排列，得到矩阵U₁和矩阵U₂的过程如下：

根据式(12)重构出新的方向矩阵B为：

则式(12)重写为：

对R进行奇异值分解：

其中U和V分别表示R的左右奇异矢量，且它们都是酋矩阵，Σ为对角矩阵，其对角线元素表示奇异值且它们从大到小排列，在无噪声的情况下，前K个奇异值大于零，剩余的奇异值为零，在有噪声影响的情况下，前K个奇异值也远远大于剩余奇异值，因此用前K个奇异值对应的奇异矢量构成信号子空间U_s，剩余(M-1-K)个奇异值对应的奇异矢量构成噪声子空间U_n，Σ_s和Σ_n对角线上的元素分别表示信号子空间和噪声子空间的奇异值，V_s和V_n分别为信号子空间和噪声子空间的右奇异矢量；

由奇异值分解性质可知：

RV_n＝Σ_nU_n (16)

两边都取共轭转置并右乘U_n得到：

由式(15)和式(17)可知，由于噪声之间的独立性，互协方差矩阵R不受噪声的影响，因此噪声子空间的奇异值都为0，即Σ_n＝0，式(17)改写为：

V_n ^HR^HU_n＝0 (18)

又V_n为矩阵R奇异值分解得到的奇异值向量矩阵，所以V_n是满秩矩阵，由满秩矩阵性质可解出式(18)得：R^HU_n＝0，并结合式(14)代入式(18)可得：

由于

是满秩的，因此：

B^HU_n＝0 (20)

由于U＝[U_s,U_n]为酋矩阵，所以有：

结合式(20)和(21)可知，方向矩阵B与互协方差矩阵R的信号特征向量组成的子矩阵U_s所张成的子空间相同，因此存在一个非奇异矩阵T使得：

U_s＝BT (22)

由式(13)可以得知方向矩阵B为2M×K阶矩阵，由两个M×K的矩阵A_x和A_xΩ组成，得出A_x具有范德蒙德结构，定义两个子阵B₁和B₂，B₁由A_x的前M-1行和A_xΩ的前M-1行组成，B₂由A_x的后M-1行和A_xΩ的后M-1行组成，即：

由式(13)和(23)得到：

B₂＝B₁Ψ (24)

其中Ψ＝diag{ψ₁,ψ₂,...,ψ_K}且

Ψ叫做旋转矩阵其对角线元素为相位旋转算子，通过求解旋转算子ψ_i求出θ_xi，将U_s按照相同的方法划分为U₁和U₂，即：

由U_s＝BT可得到：

4.根据权利要求3所述的一种未知声速环境下快速DOA估计方法，其特征在于，所述的步骤S3中对矩阵U₁ ⁺U₂进行特征值分解，进行处理后求出目标信号在在x轴上的波达方向角θ_xi的表达式和y轴上的波达方向角θ_yi的表达式的过程如下：

由式(24)和(26)推出：

U₂＝B₁ΨT＝B₁TT^-1ΨT＝U₁T^-1ΨT (27)

其中T^-1为T的逆矩阵，定义U₁ ⁺为矩阵U₁的Moore-Penrose广义逆，公式(27)两边左乘U₁ ⁺有：

U₁ ⁺U₂＝T^-1ΨT (28)

则矩阵U₁ ⁺U₂为矩阵Ψ的相似变换，所以它们有相同的特征值，而Ψ由于是对角矩阵因此其特征值就是其对角线元素

求解出U₁ ⁺U₂的K个特征值λ₁,λ₂,...,λ_K和对应的特征向量α₁,α₂,...,α_K，最后求出θ_xi的表达式为：

根据式(28)可知，矩阵T^-1由U₁ ⁺U₂的K个特征向量α₁,α₂,...,α_K组成，则可得到估计值

由(22)可以得出方向矩阵B的估计值为：

由方向矩阵B的展开式(11)得到如下表达式：

其中

表示矩阵

的第k行第i列的元素。

5.根据权利要求4所述的一种未知声速环境下快速DOA估计方法，其特征在于，所述的步骤S4中对得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式的过程如下：

用公式(32)和公式(29)相除得到：

6.根据权利要求5所述的一种未知声速环境下快速DOA估计方法，其特征在于，所述的步骤S5中根据两条线阵的几何关系推导出θ_xi和θ_yi的关系式，最终求解出目标信号的波达方向角的过程如下：

通过结合由两条线阵之间的几何关系得出的公式(5)求出θ_xi的表达式如下：

根据第i个信号的波达方向角θ_i就是其在x轴上的波达方向角θ_xi，得到最终的一维波达方向估计结果为：

Images