CN109581275B - 基于非圆信号和三维正交阵的二维水下doa估计方法和装置 - Google Patents

Abstract

本发明公开了一种基于非圆信号和三维正交阵的二维水下DOA估计方法和装置，该估计装置在结构上选取了三维正交均匀阵列，这种阵列通过接收三个维度上子阵列的数据，能够最终得到声速无关的二维DOA估计，有效解决了声速偏差带来的估计精度下降问题，可以达到提升估计性能的目的。另外为了克服水声环境中信号快速衰减的问题，该估计装置通过将非圆信号应用到水下DOA估计中，采用了基于非圆信号的NC‑ESPRIT算法，以提升估计性能。综上，该发明结合非圆信号的优点，同时采用了三维正交阵列，根据水下环境进行了声速无关的优化，估计精度高，具有较强的实用性。

Description

基于非圆信号和三维正交阵的二维水下DOA估计方法和装置

技术领域

本发明涉及目标定位的技术领域，具体涉及一种基于非圆信号和三维正交阵的二维水下DOA估计方法和装置。

背景技术

阵列信号处理技术在众多领域已得到广泛应用，而阵列信号处理的基本问题之一是空间信号波达方向估计(DOA估计)。而水下DOA估计则是指在水面放置传感器阵列利用阵列信号处理技术来对水下目标物进行方位估计的方法。

水下DOA估计采用声波作为传播载体，由于声波信号在水下环境传播时，水声信道中的各种障碍物及崎岖不平的海底造成的声波散射作用，导致了信号的急剧衰减。由此可见，水声信道的信号衰减限制了水声信道在远距离通信中的应用。除了水声环境造成信号的快速衰减，水下DOA估计面临的另一个问题就是声速影响。DOA估计算法的原理是利用阵列接收信号之间的波程差估计空间信号的方向信息。由于河流和海洋等水下环境复杂且不稳定，声波的速度随位置和时间而变化，水下DOA算法的估计精度受到很大影响。目前水下DOA估计方法普遍假定声速为已知的固定量，这将影响波程差的精度。当实际声速偏离预先设定速度，估计精度将因此降低。

为了提高水下DOA估计算法性能，采用基于非圆信号的DOA估计算法成为一个重要研究方向。在非圆信号的伪协方差矩阵不为零的特征用到信号估计中，相当于虚拟扩展阵列，能够显著的提高估计性能。同时非圆信号虚拟阵元的增加也使得算法可以处理多于阵列个数的信源个数，在复杂的水声环境中具有广阔的应用前景。

目前二维DOA估计方法所采用的阵列，多是L型阵列、圆阵和矩形阵列这类平面阵列。平面阵列在常规条件能获得比较好的估计效果，但用于声速变化的水下环境时会损失一定估计精度。

发明内容

本发明的目的是为了解决现有技术中的上述缺陷，提供一种基于非圆信号和三维正交阵的二维水下DOA估计方法和装置，该方法通过对三个维度上均匀子线阵的接收信号进行处理，在二维DOA波达方向估计中消除声速这个因子，从而消除水下声速不确定性对目标定位精度的影响。同时由于选取了非圆信号进行发射和接收，显著的提高了估计性能。

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

一种基于非圆信号和三维正交阵的二维水下DOA估计方法，所述的估计方法包括：

S1、建立十字正交线阵的阵列信号模型，所述的阵列信号模型为一个十字正交线阵和垂直线阵3的组合，所述的十字正交线阵由线阵1和线阵2组成，其中，线阵1排布于坐标系x轴上，线阵2排布于坐标系y轴上，垂直线阵3排布于坐标系z轴上，线阵1和线阵2都各有2M-1个接收阵元，垂直线阵3有M个接收阵元，所有相邻阵元之间的平均间距为d，将中心频率为f，非圆率为ρ,0＜ρ≤1的非圆信号作为发射信号，同时非圆信号满足窄带条件，即当非圆信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移，以坐标系原点为参考点，假设水下目标总个数为K，第k个目标的方位角和仰角表示为θ_k和φ_k，θ_k∈[0,π]，同时目标与坐标系x轴和y轴的夹角分别α_k和β_k，快拍数为L，线阵1和线阵2的接收数据矩阵分别表示为X和Y：

X＝A_xS+N_x (1)

Y＝A_yS+N_y (2)

其中S是一个K×L维的源信号矩阵，A_x和A_y是(2M-1)×K维导向向量矩阵，N_x和N_y是(2M-1)×L维的噪声矩阵；

S2、建立垂直线阵3的阵列信号模型，垂直线阵3有M个接收阵元，垂直线阵3的接收数据矩阵可以表示成Z：

Z＝A_zS+N_z (3)

其中A_z是M×K维导向向量矩阵，N_z是M×L维的噪声矩阵；

S3、采用基于非圆信号的DOA算法求出线阵1、线阵2和垂直线阵3对应含声速信息的特征值参数u_k、v_k和w_k,k＝1,2,…,K；

S4、利用子空间投影角度配对法进行三组特征值参数u_k,v_k和w_k,k＝1,2,…,K的配对；

S5、求解目标的二维波达方向估计解，即对于k，k＝1,2,…,K个目标的方位角θ_k和仰角φ_k的估计值。

进一步地，所述的步骤S1中，

根据发射信号的非圆特性有S＝ΦS_R，其中S_R为信源信号的实部，为发射信号的非圆相位，将公式(1)和公式(2)写成

X＝A_xΦS_R+N_x (4)

Y＝A_yΦS_R+N_y (5)

由于阵列接收的非圆信号回波同时也是窄带信号，以旋转点为原点，导向向量矩阵A_x、A_y的表达式写为：

A_x＝[a_x(θ₁,φ₁) a_x(θ₂,φ₂)…a_x(θ_K,φ_K)] (6)

A_y＝[a_y(θ₁,φ₁) a_y(θ₂,φ₂)…a_y(θ_K,φ_K)] (7)

对于第k个目标，则有

a_x(θ_k,φ_k)＝[a_x,-M+1(θ_k,φ_k)…a_x,-1(θ_k,φ_k) a_x,0(θ_k,φ_k) a_x,1(θ_k,φ_k)…a_x,M-1(θ_k,φ_k)]^T (8)

a_y(θ_k,φ_k)＝[a_y,-M+1(θ_k,φ_k)…a_y,-1(θ_k,φ_k) a_y,0(θ_k,φ_k) a_y,1(θ_k,φ_k)…a_y,M-1(θ_k,φ_k)]^T (9)

根据线阵1、线阵2与坐标轴的夹角关系，得到角度α,β关于方位角和仰角的表达式

因此得到

其中λ_k为声波的波长，即均匀线阵两相邻阵元之间的间距d要小于声波信号的半波长，而声波在探测路径上的速度v是未知的，因此取v为其范围中的最小值以确定λ_k的值。

进一步地，所述的步骤S2中，

根据发射信号的非圆特性，将公式(3)写成

Z＝A_zΦS_R+N_z (13)

由于阵列接收的非圆信号回波同时也是窄带信号，以旋转点为原点，A_z的表达式写为：

A_z＝[a_z(θ₁,φ₁) a_z(θ₂,φ₂)…a_z(θ_K,φ_K)] (14)

对于第k个目标，则有

a_z(θ_k,φ_k)＝[a_z，0(θ_k,φ_k)…a_z，M-1(θ_k,φ_k)]^T (15)

进一步地，所述的步骤S3中，

针对线阵1重构出接收信号阵列W_x：

其中J为一个行交换矩阵，

构建W_x的协方差矩阵R_w，其表示为

其中R_s是源信号实部S_R的协方差矩阵，是噪声分量的方差，I_2M为(4M-2)×(4M-2)的单位矩阵。

对协方差矩阵R_w进行特征值分解，得到信号子空间矩阵因为信号子空间U_s与B_x的关系有：span{U_s}＝span{B_x}，因此存在一个满秩矩阵T使得U_sT＝B_x，定义矩阵T₁＝[0_(M-1)×1 I_M-1]、T₂＝[I_M-1 0_(M-1)×1]以及行交换矩阵/>公式中/>为(M-1)×M维的零矩阵，

协方差矩阵R_w的估计值通过采样得到

其中L为快拍数，

构建并再次进行特征分解，得到正交矩阵：

即得到特征参数u_k,k＝1,2,…K；

同理可得线阵2和垂直线阵3对应的特征参数v_k和w_k：

进一步地，所述的步骤S4过程如下：

首先构造矩阵Q＝[X₁,Y₁,Z]^T＝A·S+N_Q，N_Q为对应的3M×L维噪声矩阵，X₁为X的(M,2M-1)行数据，Y₁为Y的(M,2M-1)行数据，导向向量矩阵A由A_x1，A_y1和A_z构成

A＝[A_x1,A_y1,A_z]^T (24)

其中，A_x1为A_x的(M,2M-1)行，A_y1为A_y的(M,2M-1)行；

获取Q的协方差矩阵R_Q有：

R_Q＝Q·Q^H (25)

对R_Q进行特征值分解，获得对应的噪声子空间U_NQ，由于导向向量矩阵与噪声子空间存在正交关系即：A^H·U_NQ＝0和构造代价函数F

构造导向向量a(θ_i,φ_i)＝[u,v,w]^T，其中u,v,w都是K×1维的向量，分别为三组特征值参数u_k,v_k和w_k的一种排列组合，一共组合出K³组导向向量，其中代价函数最大值的K组所对应的组合即为参数配对成功的组合。

进一步地，所述的步骤S5过程如下：

对于第k个目标，求出方位角θ_k和仰角φ_k的估计值，根据线阵1、线阵2以及垂直线阵3的接收信号X，Y和Z，并利用一维NC-ESPRIT算法分别求出三个子线阵对应的参数u_k，v_k和w_k,k＝1,2,…,K，得到方位角θ_k和仰角φ_k的估计：

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

一种基于非圆信号和三维正交阵的二维水下DOA估计装置，所述的估计装置包括数据处理与控制模块、以及分别与数据处理与控制模块相连的发射模块、接收模块、输出模块和电源模块，其中，所述的数据处理与控制模块包括依次连接的A/D转换器、D/A转换器和处理器，

所述的发射模块包括依次连接的功率放大器、阻抗匹配电路和超声波发射探头，通过D/A转换器与处理器相连，根据处理器发出的指令发射指定的非圆信号；

所述的接收模块采用十字正交线阵的阵列信号模型，所述的阵列信号模型为一个十字正交线阵和垂直线阵3的组合，所述的十字正交线阵由线阵1和线阵2组成，其中，线阵1排布于坐标系x轴上，线阵2排布于坐标系y轴上，垂直线阵3排布于坐标系z轴上，线阵1和线阵2都各有2M-1个接收阵元，垂直线阵3有M个接收阵元，所有相邻阵元之间的平均间距为d；

所述的输出模块包括USB接口和显示器，将数据处理与控制模块中处理完毕的数据通过USB接口输出到外部装置或者显示器进行显示；

所述的电源模块分别与数据处理与控制模块、发射模块、接收模块和输出模块相连并进行供电。

进一步地，所述的线阵1、线阵2和垂直线阵3通过塑料材质的固定支架连接。

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

1、本发明基于非圆信号的NC-ESPRIT算法应用于二维水下波达方向估计，充分利用了信号的非圆特性，可以获得相当于阵列孔径扩展的效果，使得水下DOA估计结果更加精确。不仅如此，虚拟阵元的增加也使得本发明估计方法在相同阵元数的条件下能够估计更多的信源数。

2、与利用传统的水下二维DOA算法相比，本发明应用了三维正交均匀阵列，通过3个子阵列与波达方向角之间的角度关系消去了声速这个变量，使得最后的运算结果与声速无关，消除声速偏差，使得估计的精确度更高，在未知声速环境进行DOA估计时更具优势。

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

附图说明

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

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

图3是本发明装置中三维正交均匀线阵的阵元排列图；

图4是本发明装置中三维正交均匀线阵列与目标模型示意图；

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

图6是本发明公开的基于非圆信号和三维正交阵的二维水下DOA估计方法流程图。

具体实施方式

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

实施例一

本实施例提出了一种基于非圆信号和三维正交均匀阵的二维水下波达方向估计方法，通过对三个维度上均匀子线阵的接收信号进行处理，在二维DOA波达方向估计中消除声速这个因子，从而消除水下声速不确定性对目标定位精度的影响。同时由于选取了非圆信号进行发射和接收，显著的提高了估计性能。

如附图6所示，本实施例中基于非圆信号和三维正交均匀线阵在未知声速环境进行水下二维DOA估计方法包括以下步骤：

S1、建立十字正交线阵的阵列信号模型。放置如图4所示的三维正交线阵，可以看作是一个十字正交线阵和垂直线阵3的组合，十字正交线阵的子线阵分别设为线阵1和线阵2，其中线阵1排布于坐标系x轴上，线阵2排布于坐标系y轴上，垂直线阵3排布于坐标系z轴上。线阵1和线阵2都各有2M-1个接收阵元，垂直线阵3有M个接收阵元，所有相邻阵元之间的平均间距为d。中心频率为f，非圆率为ρ(0＜ρ≤1)的非圆信号作为发射信号，同时信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。以坐标系原点为参考点，假设水下目标总个数为K，第k个目标的方位角和仰角可表示为θ_k和φ_k(θ_k∈[0,π]，)，同时目标与坐标系x轴和y轴的夹角分别α_k和β_k。快拍数为L，线阵1和线阵2的接收数据矩阵可以分别表示为X和Y：

X＝A_xS+N_x (1)

Y＝A_yS+N_y (2)

其中S是一个K×L维的源信号矩阵，A_x和A_y是(2M-1)×K维导向向量矩阵，N_x和N_y是(2M-1)×L维的噪声矩阵；

发射信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移；

根据发射信号的非圆特性有S＝ΦS_R，其中S_R为信源信号的实部，为发射信号的非圆相位，将公式(1)和公式(2)写成

X＝A_xΦS_R+N_x (3)

Y＝A_yΦS_R+N_y (4)

由于阵列接收的非圆信号回波同时也是窄带信号，以旋转点为原点，2个导向向量矩阵的表达式写为：

A_x＝[a_x(θ₁,φ₁) a_x(θ₂,φ₂)…a_x(θ_K,φ_K)] (5)

A_y＝[a_y(θ₁,φ₁) a_y(θ₂,φ₂)…a_y(θ_K,φ_K)] (6)

对于第k个目标，则有

a_x(θ_k,φ_k)＝[a_x,-M+1(θ_k,φ_k)…a_x,-1(θ_k,φ_k) a_x,0(θ_k,φ_k) a_x,1(θ_k,φ_k)…a_x,M-1(θ_k,φ_k)]^T

(7)

a_y(θ_k,φ_k)＝[a_y,-M+1(θ_k,φ_k)…a_y,-1(θ_k,φ_k) a_y,0(θ_k,φ_k) a_y,1(θ_k,φ_k)…a_y,M-1(θ_k,φ_k)]^T

(8)

根据线阵1、线阵2与坐标轴的夹角关系，得到角度α,β关于方位角和仰角的表达式

因此可以得到

其中λ_k为声波的波长，即均匀线阵两相邻阵元之间的间距d要小于声波信号的半波长。而声波在探测路径上的速度v是未知的，因此取v为其范围中的最小值以确定λ_k的值。

S2、建立垂直线阵的阵列信号模型。如图4所示，垂直线阵3有M个接收阵元，垂直线阵3的接收数据矩阵可以表示成Z：

Z＝A_zS+N_z (12)

A_z是M×K维导向向量矩阵，N_z是M×L维的噪声矩阵。根据发射信号的非圆特性，将公式(12)写成

Z＝A_zΦS_R+N_z (13)

由于阵列接收的非圆信号回波同时也是窄带信号，以旋转点为原点，A_z的表达式写为：

A_z＝[a_z(θ₁,φ₁) a_z(θ₂,φ₂)…a_z(θ_K,φ_K)] (14)

对于第k个目标，则有

a_z(θ_k,φ_k)＝[a_z，0(θ_k,φ_k)…a_z，M-1(θ_k,φ_k)]^T (15)

S3、采用基于非圆信号的DOA算法求出线阵1、线阵2和垂直线阵3对应含声速信息的特征值参数u_k、v_k和w_k,k＝1,2,…,K；以线阵1为例，可以重构出接收信号阵列W_x：

其中J为一个行交换矩阵，

构建W_x的协方差矩阵R_w，其表示为

其中R_s是源信号实部S_R的协方差矩阵，是噪声分量的方差，I_2M为_{(4M-2)×(4M-2)}的单位矩阵；

对协方差矩阵R_w进行特征值分解，得到信号子空间矩阵因为信号子空间U_s与B_x的关系有：span{U_s}＝span{B_x}，因此存在一个满秩矩阵T使得U_sT＝B_x，定义矩阵T₁＝[0_(M-1)×1 I_M-1]、T₂＝[I_M-1 0_(M-1)×1]以及行交换矩阵/>公式中/>为(M-1)×M维的零矩阵。

实际情况中，协方差矩阵R_w的估计值通过采样得到

其中L为快拍数，

构建并再次进行特征分解，得到正交矩阵：

即得到特征参数u_k,k＝1,2,…K；同理可得线阵2和垂直线阵3对应的特征参数v_k和w_k：

S4、将三组特征值参数(即u_k,v_k和w_k,k＝1,2,…,K)进行参数配对；为了实现三维阵列的三组参数(即u_k,v_k和w_k,k＝1,2,…,K)的成功配对，本实施例采用了一种基于子空间投影角度的3参数配对，适用于本实施例提出的阵列构型；首先构造矩阵Q＝[X₁,Y₁,Z]^T＝A·S+N_Q(X₁为X的(M,2M-1)行数据，Y₁为Y的(M,2M-1)行数据)，其中N_Q为对应的3M×L维噪声矩阵，导向向量矩阵A由A_x1(A_x的(M,2M-1)行)，A_y1(A_y的(M,2M-1)行)和A_z构成

A＝[A_x1,A_y1,A_z]^T (24)

获取Q的协方差矩阵R_Q有：

R_Q＝Q·Q^H (25)

对R_Q进行特征值分解，获得对应的噪声子空间U_NQ，由于导向向量矩阵与噪声子空间存在正交关系即：A^H·U_NQ＝0和可以构造代价函数F

构造导向向量a(θ_i,φ_i)＝[u,v,w]^T，其中u,v,w都是K×1维的向量，分别为三组特征值参数u_k,v_k和w_k的一种排列组合，一共可以组合出K³组导向向量，其中代价函数最大值的K组所对应的组合即为参数配对成功的组合。

S5、求解目标的二维波达方向估计解，即对于k，k＝1,2,…,K个目标的方位角θ_k和仰角φ_k的估计值。对于第k个目标，求出方位角θ_k和仰角φ_k的估计值。根据三个子线阵的接收信号X，Y和Z，并利用一维NC-ESPRIT算法分别求出三个子线阵对应的参数u_k,v_k和w_k,k＝1,2,…,K。因此可以得到方位角θ_k和仰角φ_k的估计：

实施例二

本实施例公开了一种基于非圆信号和三维正交均匀线阵的水下一维DOA估计装置，所述的估计装置包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块，具体结构如图1和图2所示。

数据处理与控制模块由一对多路A/D、D/A转换器和一个处理器组成，是整个装置的核心部分，其它所有模块都与它直接相连。它可以控制发射模块，使发射模块发射指定的信号；同时能够对接收模块传过来的信号进行处理，通过本发明的算法计算出波达方向角，然后将结果传输至输出模块。

接收模块包括3个以均匀间距摆放的超声波探头阵列。附图2为阵列装置连接图，如图所示，线阵1、线阵2和垂直线阵3连接在一起。

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

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

电源模块分别与数据处理与控制模块、发射模块、接收模块和输出模块相连并进行供电。

接收模块采用十字正交线阵的阵列信号模型，所述的阵列信号模型为一个十字正交线阵和垂直线阵3的组合，所述的十字正交线阵由线阵1和线阵2组成，其中，线阵1排布于坐标系x轴上，线阵2排布于坐标系y轴上，垂直线阵3排布于坐标系z轴上，线阵1和线阵2都各有2M-1个接收阵元，垂直线阵3有M个接收阵元，所有相邻阵元之间的平均间距为d。

因为接收模块会放置在水中，所以固定支架采用塑料材质以增大浮力。

本发明装置的主要工作流程如下：在实测过程中根据想要发射的信号参数，通过数据处理与控制模块输入对应的参数，使处理器产生相应的数字信号，然后通过D/A转换后传给发射模块，超声波发射探头就能产生需要的信号并进行发射。接收模块中的接收阵列收到从目标声源反射回来的信号后将其通过A/D转换成数字信号后发送给处理器，然后处理器根据水下DOA估计方法计算出结果。最后数据处理与控制模块将计算结果传给输出模块，输出模块将结果通过USB接口传给外部设备或者通过显示器显示出来。电源模块为所有其它模块供电。

实施例三

本实施例具体公开一种基于非圆信号和三维正交均匀线阵的水下一维DOA估计装置，包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块。数据处理与控制模块可以用DSP芯片实现(如：TI公司TMS320VC5509A型号的DSP芯片)，此DSP芯片可实现A/D转换和D/A转换的功能，并能够实现最终波达方向的计算；3个均匀阵列按图3所示组装；发射模块使用一个超声波发射探头；输出模块使用一个USB接口和一个LCD显示屏。图1即为本发明所述装置的硬件结构模块图。

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

步骤T1、按图2连接好具体装置，其中接收模块中的线阵1和线阵2中的阵元个数定为15，垂直线阵3中的阵元个数定为8。线阵1、线阵2和垂直线阵3均是均匀线阵。利用数据处理与控制模块发送指令，控制超声发射探头发射超声信号s(t)，发射信号为初相位为20°，非圆率ρ＝1的BPSK信号，信号的频率为f_s＝10kHz，脉冲长度5ms；海水中声速范围大致为1430m/s-1550m/s，则取最小声速为1430m/s，可以求出最小半波长为7.15cm。任意两相邻线阵之间的距离必须小于7.15cm，在满足此限制条件下可以任意选取阵元间距，所以设置子线阵中2个相邻阵元的平均间距为5cm，即第一个阵元和最后一个阵元相隔35cm。在水下放置一个目标声源，入射到三维正交均匀阵列的二维入射角度为(60°,45°)。

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

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

1)根据接收到的信号分别得出3个均匀线阵的接收信号矩阵X、Y和Z，随后用基于非圆信号的一维NC-ESPRIT算法求出对应的参数u_k、v_k和w_k,k＝1,2,…,K。

2)利用求出的参数u_k、v_k和w_k,k＝1,2,…,K，进行3组参数的配对。根据代价函数最大值原则，对参数u_k、v_k和w_k，k＝1,2,…,K进行全组合遍历，每一种组合都根据代价函数进行比较，对于最大的K个值所对应的组合，即为配对成功的组合。

3)求解目标的二维波达方向估计解，即对于第k个目标，求出目标角度θ_k的估计值。利用已经配对成功的参数，根据公式(27)分别求出共K个二维波达方向角度(方位角θ和仰角φ)，并将本次估计的结果进行存储。

步骤T4、将求得平均值后的最终估计结果并传送给输出模块，使其通过USB接口输出给外部装置或者显示在LCD显示屏上。根据本专利提出的方法和装置，目标估计结果为(58.86°,45.21°)，达到了预期精度，说明估计结果正确，本发明方法及装置可行。

综上所述，上述实施例为了克服水声环境中信号快速衰减的问题，通过将非圆信号应用到水下DOA估计中，采用基于非圆信号的NC-ESPRIT算法，以提升估计性能。同时上述实施例为了消除声速影响造成的估计偏差，使用一个三维正交均匀线阵作为接收阵列，根据线阵之间的角度关系，采用了与声速无关的二维DOA估计表达式，进一步提升了水下DOA估计精度。本发明结合了非圆信号的优点，结合水下环境利用三维正交均匀线阵进行了声速无关的优化，估计精度高，具有较强的实用性。

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

Claims (5)

1.一种基于非圆信号和三维正交阵的二维水下DOA估计方法，其特征在于，所述的估计方法包括：

S1、建立十字正交线阵的阵列信号模型，所述的阵列信号模型为一个十字正交线阵和垂直线阵3的组合，所述的十字正交线阵由线阵1和线阵2组成，其中，线阵1排布于坐标系x轴上，线阵2排布于坐标系y轴上，垂直线阵3排布于坐标系z轴上，线阵1和线阵2都各有2M-1个接收阵元，垂直线阵3有M个接收阵元，所有相邻阵元之间的平均间距为d，将中心频率为f，非圆率为ρ,0＜ρ≤1的非圆信号作为发射信号，同时非圆信号满足窄带条件，即当非圆信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移，以坐标系原点为参考点，假设水下目标总个数为K，第k个目标的方位角和仰角表示为θ_k和φ_k，θ_k∈[0,π]，同时目标与坐标系x轴和y轴的夹角分别α_k和β_k，快拍数为L，线阵1和线阵2的接收数据矩阵分别表示为X和Y：

X＝A_xS+N_x (1)

Y＝A_yS+N_y (2)

其中S是一个K×L维的源信号矩阵，A_x和A_y是(2M-1)×K维导向向量矩阵，N_x和N_y是(2M-1)×L维的噪声矩阵；

根据发射信号的非圆特性有S＝ΦS_R，其中S_R为信源信号的实部， 为发射信号的非圆相位，将公式(1)和公式(2)写成

X＝A_xΦS_R+N_x (4)

Y＝A_yΦS_R+N_y (5)

由于阵列接收的非圆信号回波同时也是窄带信号，以旋转点为原点，导向向量矩阵A_x、A_y的表达式写为：

A_x＝[a_x(θ₁,φ₁) a_x(θ₂,φ₂)…a_x(θ_K,φ_K)] (6)

A_y＝[a_y(θ₁,φ₁) a_y(θ₂,φ₂)…a_y(θ_K,φ_K)] (7)

对于第k个目标，则有

a_x(θ_k,φ_k)＝[a_x,-M+1(θ_k,φ_k)…a_x,-1(θ_k,φ_k) a_x,0(θ_k,φ_k) a_x,1(θ_k,φ_k)…a_x,M-1(θ_k,φ_k)]^T (8)

a_y(θ_k,φ_k)＝[a_y,-M+1(θ_k,φ_k)…a_y,-1(θ_k,φ_k) a_y,0(θ_k,φ_k) a_y,1(θ_k,φ_k)…a_y,M-1(θ_k,φ_k)]^T (9)

根据线阵1、线阵2与坐标轴的夹角关系，得到角度α，β关于方位角和仰角的表达式

因此得到

其中λ_k为声波的波长，即均匀线阵两相邻阵元之间的间距d要小于声波信号的半波长，而声波在探测路径上的速度v是未知的，因此取v为其范围中的最小值以确定λ_k的值；

S2、建立垂直线阵3的阵列信号模型，垂直线阵3有M个接收阵元，垂直线阵3的接收数据矩阵表示成Z：

Z＝A_zS+N_z (3)

A_z是M×K维导向向量矩阵，N_z是M×L维的噪声矩阵；

根据发射信号的非圆特性，将公式(3)写成

Z＝A_zΦS_R+N_z (13)

由于阵列接收的非圆信号回波同时也是窄带信号，以旋转点为原点，A_z的表达式写为：

A_z＝[a_z(θ₁,φ₁) a_z(θ₂,φ₂)…a_z(θ_K,φ_K)] (14)

对于第k个目标，则有

a_z(θ_k,φ_k)＝[a_z,0(θ_k,φ_k)…a_z,M-1(θ_k,φ_k)]^T (15)

S3、采用基于非圆信号的DOA算法求出线阵1、线阵2和垂直线阵3对应含声速信息的特征值参数u_k、v_k和w_k，k＝1,2,…,K；

针对线阵1重构出接收信号阵列W_x：

其中J为一个行交换矩阵，

构建W_x的协方差矩阵R_w，其表示为

其中R_s是源信号实部S_R的协方差矩阵，是噪声分量的方差，I_2M为(4M-2)×(4M-2)的单位矩阵；

对协方差矩阵R_w进行特征值分解，得到信号子空间矩阵因为信号子空间U_s与B_x的关系有：span{U_s}＝span{B_x}，因此存在一个满秩矩阵T使得U_sT＝B_x，定义矩阵T₁＝[0_(M-1)×1I_M-1]、T₂＝[I_M-1 0_(M-1)×1]以及行交换矩阵/>公式中/>为(M-1)×M维的零矩阵，

协方差矩阵R_w的估计值通过采样得到

其中L为快拍数，

构建并再次进行特征分解，得到正交矩阵：

即得到特征参数u_k,k＝1,2,…K；

同理得线阵2和垂直线阵3对应的特征参数v_k和w_k：

S4、利用子空间投影角度配对法进行三组特征值参数u_k、v_k和w_k，k＝1,2,…,K的配对；

S5、求解目标的二维波达方向估计解，即对于第k,k＝1,2,…,K个目标的方位角θ_k和仰角φ_k的估计值。

2.根据权利要求1所述的基于非圆信号和三维正交阵的二维水下DOA估计方法，其特征在于，所述的步骤S4过程如下：

首先构造矩阵Q＝[X₁,Y₁,Z]^T＝A·S+N_Q，N_Q为对应的3M×L维噪声矩阵，X₁为X的(M,2M-1)行数据，Y₁为Y的(M,2M-1)行数据，导向向量矩阵A由A_x1、A_y1和A_z构成

A＝[A_x1,A_y1,A_z]^T (24)

其中，A_x1为A_x的(M,2M-1)行，A_y1为A_y的(M,2M-1)行；

获取Q的协方差矩阵R_Q有：

R_Q＝Q·Q^H (25)

对R_Q进行特征值分解，获得对应的噪声子空间U_NQ，由于导向向量矩阵与噪声子空间存在正交关系即：A^H·U_NQ＝0和构造代价函数F

构造导向向量a(θ_i,φ_i)＝[u,v,w]^T，其中u,v,w都是K×1维的向量，分别为三组特征值参数u_k，v_k和w_k的一种排列组合，一共组合出K³组导向向量，其中代价函数最大值的K组所对应的组合即为参数配对成功的组合。

3.根据权利要求2所述的基于非圆信号和三维正交阵的二维水下DOA估计方法，其特征在于，所述的步骤S5过程如下：

对于第k个目标，求出方位角θ_k和仰角φ_k的估计值，根据线阵1、线阵2以及垂直线阵3的接收信号X、Y和Z，并利用一维NC-ESPRIT算法分别求出三个子线阵对应的参数u_k，v_k和w_k，k＝1,2,…,K，得到方位角θ_k和仰角φ_k的估计：

4.一种根据权利要求1至3任一所述的基于非圆信号和三维正交阵的二维水下DOA估计方法的估计装置，其特征在于，所述的估计装置包括数据处理与控制模块、以及分别与数据处理与控制模块相连的发射模块、接收模块、输出模块和电源模块，其中，所述的数据处理与控制模块包括依次连接的A/D转换器、D/A转换器和处理器，

所述的发射模块包括依次连接的功率放大器、阻抗匹配电路和超声波发射探头，通过D/A转换器与处理器相连，根据处理器发出的指令发射指定的非圆信号；

所述的接收模块采用十字正交线阵的阵列信号模型，所述的阵列信号模型为一个十字正交线阵和垂直线阵3的组合，所述的十字正交线阵由线阵1和线阵2组成，其中，线阵1排布于坐标系x轴上，线阵2排布于坐标系y轴上，垂直线阵3排布于坐标系z轴上，线阵1和线阵2都各有2M-1个接收阵元，垂直线阵3有M个接收阵元，所有相邻阵元之间的平均间距为d；

所述的输出模块包括USB接口和显示器，将数据处理与控制模块中处理完毕的数据通过USB接口输出到外部装置或者显示器进行显示；

所述的电源模块分别与数据处理与控制模块、发射模块、接收模块和输出模块相连并进行供电。

5.根据权利要求4所述的基于非圆信号和三维正交阵的二维水下DOA估计装置，其特征在于，所述的线阵1、线阵2和垂直线阵3通过塑料材质的固定支架连接。