CN105403874B

CN105403874B - 非均匀阵列欠定波达方向估计方法

Info

Publication number: CN105403874B
Application number: CN201510990239.2A
Authority: CN
Inventors: 杨明磊; 陈伯孝; 孙磊
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2015-12-25
Filing date: 2015-12-25
Publication date: 2017-11-21
Anticipated expiration: 2035-12-25
Also published as: CN105403874A

Abstract

本发明公开了一种基于非均匀阵列欠定的波达方向估计方法，主要解决现有技术计算复杂度较高的问题。其技术方案是：首先利用接收阵列数据的协方差矩阵构造虚拟阵列；然后利用全部虚拟阵元的数据获得差分合成阵列的观测数据；接着通过构造一个观测数据的Toeplitz矩阵代替传统的空间平滑操作进行解相干，并估计信号和噪声的子空间；最后利用噪声子空间构造一个关于MUSIC谱的多项式并利用其解来估计波达方向。本发明避免了大运算量的角度栅格搜索，以更低的计算复杂度获得比传统SS‑MUSIC算法更高的估计精度，可用于对目标方位的估计。

Description

非均匀阵列欠定波达方向估计方法

技术领域

本发明属于信号处理技术领域，尤其涉及一种波达方向估计方法，用于对目标方位的估计。

背景技术

波达方向DOA估计是雷达、声纳信号处理中的一个重要研究方向。众所周知，对于一个具有N个阵元的均匀线性阵列，采用传统的波达方向估计方法，如MUSIC，ESPRIT等所能分辨的最大目标个数为N-1个。而欠定波达方向估计问题，即目标个数大于阵元数的DOA估计经常出现并引起广泛的研究兴趣。解决该问题的其中一种有效方法是利用一个等效的虚拟阵列来提高波达方向估计的自由度，该虚拟阵列是通过对一个特殊设计的非均匀线性阵列接收信号的协方差矩阵向量化来构造的。最小冗余阵列MRA、嵌套式阵列NA、互质阵列CA及嵌套式最小冗余阵列NMRA等都是针对欠定波达方向估计而设计的非均匀阵列。但是在利用阵列接收数据的协方差矩阵构造虚拟阵列时，虚拟阵列中等效信源被实际信源的功率所替代，因此这些等效信源就表现为完全相干的信号，这样就无法直接利用传统的DOA估计方法来进行波达方向估计。

近来，许多利用非均匀线阵的虚拟阵列来解决欠定波达方向估计问题的新算法被提出。其中一种代表性方法是P.Pal和P.P.Vaidyanathan在《Nested arrays:Anovelapproach to array processing with enhanced degrees of freedom》提出的空间平滑(SS)-MUSIC算法，但该方法在进行DOA估计时没能用到虚拟阵列的所有阵元，从而导致一定的信噪比损失。另一种具有代表性的方法是Y.D.Zhang、M.G.Amin和B.Himed在《Sparsity-based DOAestimationusing co-prime arrays》提出的稀疏信号重构法，该方法利用了信号谱的稀疏性来进行DOA估计，然而实现这种方法的算法通常需要巨大的计算量。

发明内容

本发明的目的在于克服上述已有DOA估计算法的不足，提出一种新的非均匀阵列欠定波达方向估计方法，以充分利用虚拟阵列的数据，减少计算量。

为实现上述目的，本发明的技术思路是:利用接收阵列数据的协方差矩阵来构造虚拟阵列；利用全部虚拟阵元的数据来获得差分合成阵列的观测数据，而不是删除重复的阵元数据；通过构造一个观测数据的Toeplitz矩阵代替传统的空间平滑操作来解相干，并估计信号和噪声的子空间；通过构造一个关于MUSIC谱的多项式并利用其根来估计DOA，从而避免了大运算量的角度栅格搜索。其实现步骤包括如下：

1)根据阵列接收数据x(t)估计非均匀阵列的协方差矩阵并对该该协方差矩阵向量化，得到K²×1的向量

其中t表示采样时刻，t＝1,2,…,N，N表示快拍数，(·)^H表示共轭转置，K表示阵元数；

2)构造降维矩阵R：

R＝(E^TE)^-1E^T

其中是一个在第p_l处为1，其余为0的f_V×1向量，l＝1,2,…,K²，i＝1,2,…,K，f_V表示非均匀阵列的自由度，f_V＝2l_V+1，l_V是非均匀阵列的孔径长度，符号表示复数域，(·)^T表示转置；

3)根据降维矩阵R和向量计算非均匀阵列的差分合成阵列的观测数据z：

4)由观测数据z，构造Toeplitz矩阵Y：

其中z_m是观测数据z的第m个元素，m＝1,2,…,f_V；

5)对Toeplitz矩阵Y进行特征值分解，即：

其中U_S为信号子空间，Λ_S表示信号的特征值，σ²表示噪声功率，U_N表示噪声子空间；

6)根据5)所获得的噪声子空间U_N，构造root-MUSIC多项式f_MUSIC(r)：

其中r是关于目标到达角θ的未知量，L(r)表示r的l_V次多项式所构成的向量；

7)求解6)所构造的多项式f_MUSIC(r)，得到多项式的解r_n，n＝1,2,…,Q，Q为目标数，根据r_n，进而获得目标到达角θ_n：

其中arg(·)表示取相角运算，d为阵元间隔。

本发明与现有DOA估计算法相比具有以下优点：

1.本发明综合利用了所有虚拟阵元接收到的数据，避免了接收信号的信噪比SNR的损失；

2.本发明通过构造观测数据的Toeplitz矩阵代替传统的空间平滑操作来解相干，节省空间平滑解相干的运算量；

3.本发明通过构造一个关于MUSIC谱的多项式并利用其根来估计DOA，从而避免了普通谱MUSIC算法复杂的角度栅格搜索，在减少计算复杂度的同时还能获得更好地DOA估计性能。

附图说明

图1是本发明的实现流程图；

图2是本发明仿真使用的12阵元嵌套式最小冗余阵列NMRA的几何结构示意图；

图3是用本发明和SS-MUSIC算法对波达方向DOA的估计均方根误差与信噪比关系的比较图。

具体实施方式

下面结合附图详细说明本发明的内容和效果。

参照图1，本发明的实现步骤如下：

步骤1：根据阵列接收数据x(t)估计非均匀阵列的协方差矩阵

1a)对于给定K个阵元的非均匀线性天线阵列，设阵元位置为：

v＝[v₁,v₂,…,v_i,…,v_K]d,

其中v_i为第i个阵元位置系数，i＝1,2,…,K，d为阵元间隔，取值为半个信号波长；

1b)将目标n从角度θ_n入射到阵列的导向矢量表示为：

其中表示第n个信号在第i个阵元处的阵因子，n＝1,2,…,Q，Q为信号数目，κ为半波数，(·)^T表示转置；

根据入射信号的导向矢量a(θ_n)，构造导向矢量矩阵A：

A＝[a(θ₁),a(θ₂),…,a(θ_n),…,a(θ_Q)]，

1d)设Q个不相关的窄带信号从角度θ_n入射到阵列，将阵列的接收数据表示为：

x(t)＝As(t)+n(t)，

其中s(t)表示Q个不相关的窄带信号，s(t)＝[s₁(t),s₂(t),…,s_n(t),…,s_Q(t)]^T，s_n(t)表示第n个入射信号，各信号在时间上相互独立且非相关，且服从复高斯分布表示第n个信号s_n(t)的功率；n(t)表示均值为0，方差为σ²高斯白噪声，其满足独立同分布，并且与入射信号不相关；t表示采样时刻，t＝1,2,…,N，N表示快拍数；

1e)根据阵列接收数据x(t)，估计非均匀阵列的协方差矩阵

其中，I_K是K阶单位方阵，(·)^H表示共轭转置。

步骤2，向量化阵列协方差矩阵计算虚拟阵列的接收数据

其中B表示虚拟的差分合成阵列的方位矢量矩阵，

c表示等效信源矢量，1_n为单位矩阵，表示第i个元素为1，其余元素为0的单位列向量；符号vec(·)表示矩阵的向量化操作，符号代表Kronecker积，(·)^*表示共轭。

步骤3：构造降维矩阵R。

3a)规定单位向量

其中在第p_l处为1，其余为0的单位向量，p_l表示虚拟阵元位置，p_l＝1,2,…,f_V，l＝1,2,…,K²，f_V表示非均匀阵列的自由度，f_V＝2l_V+1，l_V是非均匀阵列的孔径长度；

3b)根据3a)所规定的向量生成单位矩阵E：

其中E是K²×f_V的单位矩阵；

3c)由3b)所生成的单位矩阵E，得到降维矩阵R：

R＝(E^TE)^-1E^T；

其中R是f_V×K²的矩阵。

步骤4：根据接收数据和降维矩阵R，计算非均匀阵列的差分合成阵列观测的数据z：

这里根据观测数据z得到了所有虚拟阵元的数据，但并不删除其中重复的量。

步骤5：由观测数据，构造Toeplitz矩阵Y。

5a)展开观测数据z的各个元素z_m，m＝1,2,…,2l_V+1；

5b)根据5a)中的z_m，从m＝1开始，选取连续的l_V+1个元素，生成第一个向量m逐次后移，生成第m个向量直至m＝l_V+1，最后生成第l_V+1个向量

5c)根据5b)中生成的l_V+1个列向量y_m，得到Toeplitz矩阵Y：

这里矩阵Y是个秩为Q的矩阵，如果Q＜l_V+1，则等效信源的相干性就可以被解掉，矩阵Y就能被用来正确的估计信号和噪声子空间。由于矩阵Y，能够代替传统的空间平滑操作来解相干，而且由于矩阵Y是对观测数据重新排列来构造的，所以不需要乘法，因而能够节省空间平滑解相干的运算量。

步骤6：对Toeplitz矩阵Y进行特征值分解，获得噪声子空间U_N。

定义矩阵Y的特征值分解为：

其中分别为信号和噪声子空间，Λ_S表示信号的特征值，λ_n是第n个信号的特征值，符号Diag{·}表示对角矩阵。

步骤7：根据步骤6所获得的噪声子空间U_N，构造root-MUSIC多项式f_MUSIC(r)。

7a)设关于目标波达方向θ的未知量r：

r＝e^{j(2π/λ)dsinθ}，

其中r＝[r₁,r₂,…,r_n,…,r_Q]^T，r_n表示关于目标n的未知量，θ＝[θ₁,θ₂,…,θ_n,…,θ_Q]^T，θ_n表示目标n的波达方向，λ为信号波长；

7b)生成r的l_V次多项式L(r)：

7c)结合L(r)和噪声子空间U_N，得到root-MUSIC多项式f_MUSIC(r)：

步骤8：求解步骤7构造的关于MUSIC谱的多项式f_MUSIC(r)，利用其根来估计DOA。

8a)求解步骤7所构造的多项式f_MUSIC(r)，得到其根r_n；

8b)选取离单位圆最近的Q个根r_n，估计目标n的到达角θ_n：

其中arg(·)表示取相角运算。

通过求解关于MUSIC谱的多项式并利用其根来估计DOA，避免了普通谱MUSIC算法复杂的角度栅格搜索，在减少计算复杂度的同时还能获得更好地DOA估计性能。

本发明的效果通过以下计算仿真进一步说明：

仿真1：对本发明和SS-MUSIC算法的波达方向DOA估计的均方根误差RMSE与信噪比之间的关系进行仿真。

1.1)仿真条件：图2是本发明仿真使用的12阵元嵌套式最小冗余阵列NMRA的几何结构示意图，阵元位置为[0,1,4,6,13,14,17,19,39,40,43,45]d，阵元个数为12，阵元间隔d＝λ/2，目标个数Q＝16，快拍数N＝500，500次Monte-Carlo实验。

1.2)仿真内容与结果：

用现有的SS-MUSIC算法和本发明的波达方向DOA估计性能进行比较，在上述1.1)仿真条件下，对本发明和SS-MUSIC算法的波达方向DOA估计均方根误差与信噪比关系做仿真，仿真结果如图3所示，其中横轴表示信噪比从-20dB至10dB变化，纵轴表示均方根误差。

从图3可以看出，本发明DOA估计性能随着信噪比的提高而提高，而且本发明角度估计的均方根误差小于传统SS-MUSIC，说明本发明具有更好的DOA估计精度。

仿真2：计算复杂度比较。

列举非均匀阵列的阵元数分别为12和24时，将本发明和现有的SS-MUSIC算法的计算复杂度进行比较，结果如表1。

表1本发明和SS-MUSIC算法的计算复杂度

注：1)J表示非均匀阵列波速宽度中角度栅格的个数；

2)在12阵元嵌套式最小冗余阵列NMRA中，l_V＝45,J＝1588。

3)在24阵元嵌套式最小冗余阵列NMRA中，l_V＝157,J＝6176。

从表1中可以看出，在阵元数为12的情况下，本发明的计算复杂度不足SS-MUSIC算法的一半。而且阵元数越多，复杂度降低越多。在阵元数为24的情况下，本发明的计算复杂度仅为SS-MUSIC算法的计算复杂度的1/3，可见本发明能明显降低信号处理的运算量，提高目标波达方向估计的计算速度。

综上，本发明能够以更低的计算量获得更好的DOA估计精度。

Claims

1.一种非均匀阵列欠定波达方向估计方法，包括以下步骤：

1)根据阵列接收数据x(t)估计非均匀阵列的协方差矩阵并对该协方差矩阵向量化，得到K²×1的向量

2)构造降维矩阵R：

R＝(E^TE)^-1E^T

其中是一个在第p_l处为1，其余为0的f_V×1向量，l＝1,2,…,K²，f_V表示非均匀阵列的自由度，f_V＝2l_V+1，l_V是非均匀阵列的孔径长度，(·)^T表示转置；

3)根据降维矩阵R和向量计算非均匀阵列的差分合成阵列观测的数据z：

4)由观测数据z，构造Toeplitz矩阵Y：

其中z_m是观测数据z的第m个元素，m＝1,2,…,2l_V+1；

5)对Toeplitz矩阵Y进行特征值分解，即：

<mrow> <mi>Y</mi> <mo>=</mo> <msub> <mi>U</mi> <mi>S</mi> </msub> <msub> <mi>&Lambda;</mi> <mi>S</mi> </msub> <msubsup> <mi>U</mi> <mi>S</mi> <mi>H</mi> </msubsup> <mo>+</mo> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> <msub> <mi>U</mi> <mi>N</mi> </msub> <msubsup> <mi>U</mi> <mi>N</mi> <mi>H</mi> </msubsup> </mrow>

其中r是关于目标到达角θ的未知量，L(r)表示r的l_V次多项式所构成的向量，r＝[r₁,r₂,…,r_n,…,r_Q]^T，r_n表示关于目标n的未知量，θ＝[θ₁,θ₂,…,θ_n,…,θ_Q]^T，θ_n表示第n个目标的到达角，n＝1,2,…,Q，Q为目标数；

7)求解6)所构造的多项式f_MUSIC(r)，得到多项式的解r_n，根据r_n进而获得第n个目标的到达角θ_n：

<mrow> <msub> <mi>&theta;</mi> <mi>n</mi> </msub> <mo>=</mo> <mi>a</mi> <mi>r</mi> <mi>c</mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>arg</mi> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mi>&lambda;</mi> </mrow> <mrow> <mn>2</mn> <mi>&pi;</mi> <mi>d</mi> </mrow> </mfrac> <mo>)</mo> </mrow> </mrow>

其中arg(·)表示取相角运算，d为阵元间隔。

2.根据权利要求1所述的非均匀阵列欠定波达方向估计方法，其中步骤2)中构造的降维矩阵R，按如下步骤进行:

2a)规定单位向量

其中在第p_l处为1，其余为0，p_l表示虚拟阵元位置，p_l＝1,2,…,f_V，f_V表示非均匀阵列的自由度，f_V＝2l_V+1，l＝1,2,…,K²；

2b)根据2a)所规定的生成单位矩阵E：

2c)由2b)所生成的单位矩阵E，得到降维矩阵R：

R＝(E^TE)^-1E^T。

3.根据权利要求1所述的非均匀阵列欠定波达方向估计方法，其中步骤4)中由观测数据z，构造Toeplitz矩阵Y，按如下步骤进行:

3a)展开观测数据z的各个元素z_m，m＝1,2,…,2l_V+1；

3b)根据3a)中的z_m，从m＝1开始，选取连续的l_V+1个元素，生成第一个向量m逐次后移，生成第m个向量直至m＝l_V+1，最后生成第l_V+1个向量

3c)根据3b)中生成的l_V+1个列向量y_m，得到Toeplitz矩阵Y：

4.根据权利要求1所述的非均匀阵列欠定波达方向估计方法，其中步骤6)中构造root-MUSIC多项式f_MUSIC(r)，按如下步骤进行：

4a)设定关于目标波达方向θ的未知量r：

r＝e^{j(2n/λ)dsinθ}

其中λ为信号波长，d为阵元间隔；

4b)生成r的l_V次多项式L(r)：

4c)结合L(r)和噪声子空间U_N，得到root-MUSIC多项式f_MUSIC(r)：

<mrow> <msub> <mi>f</mi> <mrow> <mi>M</mi> <mi>U</mi> <mi>S</mi> <mi>I</mi> <mi>C</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>L</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>/</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>U</mi> <mi>N</mi> </msub> <msubsup> <mi>U</mi> <mi>N</mi> <mi>H</mi> </msubsup> <mi>L</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>.</mo> </mrow> 2