CN106054123B - 一种稀疏l阵及其二维doa估计方法 - Google Patents

Abstract

本发明公开了一种稀疏L阵及其二维DOA估计方法，属于无线移动通信技术领域。本发明的稀疏L阵列包括由阵元间距等于波长的稀疏均匀线阵和一个辅助阵元构成第一子阵、由最小阵元间距小于或等于半倍波长的任意稀疏线阵构成的第二子阵，两个线阵的共有阵元为参考阵元，辅助阵元到参考阵元的距离为半倍波长。在二维DOA估计时，首先基于第二子阵的接收数据计算其自相关矩阵，并对其进行特征分解后估计对应的第二角度，再基于其计算信源自相关矩阵；基于两个子阵接收数据的互相关矩阵、信源自相关矩阵得到第一子阵的阵流行矩阵，从而完成第一子阵所对应的第一角度的估计处理，得到二维DOA。本发明的复杂程度低、DOA估计的精确度高。

Description

一种稀疏L阵及其二维DOA估计方法

技术领域

本发明属于无线移动通信技术领域，特别是涉及一种基于线阵构造的L阵列及其二维波达方向(DOA)估计方法。

背景技术

随着近年来无线通信技术的迅速发展，人们对通信业务量和通信质量的需求也越来越大，以阵列信号处理技术为核心的空分多址技术已成为下一代移动通信的关键。

现有的二维DOA估计大多是基于由阵元间距等于半倍波长的简化面阵构成的。其中，L阵由于有更大的有效孔径、更小的运算量、更易实现、更强的方法适用性等优点得到了广泛的关注和应用。近几十年来，人们已经做了很多的利用L阵估计2-D DOA的研究并提出了大量的算法。主要分为两大类：需要额外配对的算法，如文献“Nizar Tayem and Hyuck MKwon,L-shape 2dimensional arrival angle estimation with propagator method,Antennas and Propagation,IEEE Transactions on,vol.53,no.5,pp.1622–1630,2005.”和可以自动配对的算法,如文献“Jian-Feng Gu and Ping Wei,Joint svd of twocrosscorrelation matrices to achieve automatic pairing in 2-D angleestimation problems,Antennas and Wireless Propagation Letters,IEEE,vol.6,pp.553–556,2007.，简称JSVD”和”J.Gu,P.Wei,and H.-M.Tai,“DOAestimation usingcross-correlation matrix,”in Phased Array Systems and Technology(ARRAY),2010IEEE International Symposium on.IEEE,2010,pp.593–598.”，简称CCM-based。但现有的L阵大多是由常规ULA(均匀线阵)构成，存在测向精度和系统成本之间的矛盾，为了缓解该矛盾，2015年，文献“Jian-Feng Gu,Wei-Ping Zhu,and MNS Swamy,Joint 2-d doaestimation via sparse l-shaped array,Signal Processing,IEEE Transactions on,vol.63,no.5,pp.1171–1182,2015.”提出了基于SLA(稀疏线阵)和ULA的稀疏L阵，该算法和传统由常规ULA构成的L阵相比虽然在性能上有所提升，但是该阵列并没有充分利用稀疏阵列的优势，此外在利用SLA求解方位角时的性能不稳定，不能达到理想所需要的结果。

发明内容

本发明的发明目的在于：为克服现有的稀疏L阵和传统L阵在估计算法复杂度，估计精度方面的不足，提出一种结构简单的稀疏L阵及其相应的二维波达角(DOA)估计算法，以达到降低计算复杂度、系统成本、简化处理程序、有效提高估计精度等目的。

本发明为了利用稀疏阵有更大的阵列孔径，均匀线阵的平移不变性等特性而提出一种基于SLA和SULA的L阵结构，然后利用互相关矩阵不受噪声影响的特性、ULA的平移不变性、递归思想和最小二乘(LS)技术求解二维DOA。

本发明的稀疏L阵列，包括阵元数不同的线性第一子阵(M1个阵元)、第二子阵(M2个阵元)构成的L形阵列，两个均匀线阵的共有阵元定义为参考阵元，第一子阵由阵元间距等于波长的稀疏均匀线阵和一个辅助阵元构成，辅助阵元到参考阵元的距离为半倍波长，第二子阵由最小阵元间距小于或等于半倍波长的任意稀疏线阵构成，且M₁≥4，M₂≥3。本发明的稀疏L阵列(以下简称L阵)可以是位于x-z平面或y-z平面，第一子阵、第二子阵对应二维DOA的两个角：方位角、俯仰角，第一子阵可对应方位角或俯仰角，第二子阵可对应俯仰角或方位角，取决于L阵的放置方式。

在基于本发明的L阵进行二维DOA求解时，先基于第二子阵的接收数据计算其自相关矩阵，并对该自相关矩阵进行特征分解后估计对应的方位角或俯仰角，再基于其得到信源自相关矩阵；基于两个子阵接收数据的互相关矩阵、信源自相关矩阵得到第一子阵的阵流行矩阵，从而完成第一子阵所对应的俯仰角或方位角的估计处理，得到二维DOA，从而大幅度降低处理量及处理的复杂程度、有效提高DOA估计的精确度，从而实现本发明目的。用于本发明的稀疏L阵列的二维波达方向的估计方法具体包括下列步骤：

步骤1：设置天线阵列并建立系统模型：

设置L阵列的两个子阵与二维波达方向DOA的两个角的对应关系，对应第一子阵的定义为第一角度，对应第二子阵的定义为第二角度，例如第一子阵放于x轴(定义方位角)、第二子阵放于z轴(定义为俯仰角)，则第一角度对应方位角，第二子阵对应俯仰角。为了便于描述，以下用x-z平面的稀疏L阵对本发明进行说明。

在t时刻第一子阵和第二子阵的接收数据分别为：其中，分别为第一子阵和第二子阵的接收数据矢量，s(t)＝[s₁(t),…,s_K(t)]^T为信号矢量，即信号源，A_x＝[a_x(φ₁),…,a_x(φ_K)]为第一子阵的阵列流型矩阵，表示A_x的第k列导向矢量，φ_k表示第k个信源的第一角度，在x-z平面中，因x轴对应方位角，z轴对应俯仰角，φ_k表示第k个信源的方位角，λ表示信源波长，辅助阵元与参考阵元的间距d＝0.5λ，其他阵元间距d_x＝λ，e表示自然底数，j表示虚数单位。A_z＝[a_z(θ₁),…,a_z(θ_K)]为第二子阵的阵列流型矩阵，表示A_z的第k列导向矢量，k＝1,…,K，θ_k表示第k个信源的第二角度，在x-z平面中，θ_k表示第k个信源的俯仰角，d_i表示第二子阵的第i个阵元与第(i-1)个阵元之间的间距，i＝0,…,M2-1。

步骤2：计算第二子阵(SLA)的所有阵元的接收数据的自相关矩阵利用时间平均代替统计平均，可以求得第二子阵的接收数据的自相关矩阵：

步骤3：确定的噪声子空间U_zn：对步骤2所得协方差矩阵进行特征值分解，取的前K个最大特征值对应的特征向量为列构建特征向量矩阵U_zs作为信号子空间，剩余的特征值(M₂-K个)对应的特征向量为列构建特征向量矩阵U_zn作为噪声子空间；取的前K个最大特征值组成对角矩阵D_s，剩余的特征值组成对角矩阵D_n，即

步骤4：求第二角度：利用信号子空间和噪声子空间的正交性确定MUSIC算法的空间谱函数为：其中a(θ)表示关于搜索角度θ的方向矢量，θ∈[0°,180°]，即例如采用递归网格划分的方式来进行谱峰搜索：

1)第一次谱峰搜索时，先对空间谱[0°,180]°进行一个粗略的划分，θ以步长L₁从0°增长到180°，遍历搜索得到P_MUSIC的前K个最大峰值，这K个最大峰值对应的θ值就是所求的K个信号的俯仰角的大概的估计值θ′_k；

2)在上一步所得K个θ′_k附近用一个加密的网格即更小θ的增长步长，进行谱峰搜索，得到K个更精确的俯仰角估计值；

3)重复第2步，直到空间搜索的网格足够精细。

步骤5：求信源自相关矩阵首先利用步骤3所得的U_zs及D_s得到然后利用步骤4所得的俯仰角的估计值求其相应的阵列流型矩阵最后由和求得信源协方差矩阵的估计值从而建立俯仰角和方位角之间一一对应的关系。

步骤6：求互相关矩阵利用时间平均代替统计平均，求x轴全部阵元的接收数据x(t)和z轴除参考阵元以外的阵元的接收数据z′(t)之间的互相关矩阵即：

步骤7：求x轴阵列流型矩阵：利用互相关矩阵不受噪声影响的特性，将最大似然(ML)估计转化为最小二乘(LS)问题，得到A_xR_s的估计值，然后利用步骤4所得的z轴阵列流型矩阵和步骤5所得的信源自相关矩阵求x轴阵列流型矩阵的估计值即其中A_xR_s的估计值符号“A(a:b,:)”表示对应矩阵A的第a至第b行，符号(·)⁺表示M-P广义逆。

步骤8：求方位角的粗估计值：利用步骤7所得的x轴第一子阵的阵列流型矩阵的前三行数据组成的矩阵及ULA的平移不变性求得方位角的粗估计值，其中第k个方位角的粗估计值φ_k′可根据公式求得，和分别表示矩阵的第k列的前两行和最后两行的数据。

步骤9：求方位角的细估计值：将第一子阵的阵列流型矩阵的第二行数据删除得到矩阵将180°均分为3个不重叠区间：最小区间、中间区间、最大区间；

若φ_k′在最小区间，则

若φ_k′在中间区间，则

若φ_k′在最大区间，则

其中，和分别表示矩阵的第k列的前M1-2行和最后M1-2行的数据。

本发明首先利用第二子阵的接收数据求得相应的自相关矩阵，接着利用一次特征值分解，得到相应的噪声子空间U_zn，信号子空间U_zs及信号的能量矩阵D_s，再基于MUSIC算法得到高精度的俯仰角(或方位角)的估计值，从而很好的利用了SLA的稀疏性带来的阵列孔径的优势，且有效的降低了MUSIC谱峰搜索的计算量；然后利用所得的俯仰角(或方位角)得到相应的高精度的阵列流型矩阵为后面估计方位角(或俯仰角)提供更精确的参数，再利用特征值分解得到的信号子空间和信号能量矩阵及得到信号自相关矩阵的估计值从而很大程度上提高了的准确性；再然后利用不受加性噪声影响的互相关矩阵及求的第一子阵的阵列流型矩阵的估计值最后利用辅助阵元及ULA的平移不变性先后估计方位角的粗略值和精确值，既解决了SULA带来的周期性模糊问题，提高了估计精度，又降低了计算复杂度，提高了测向的效率。因而，本发明具有方法简单、可大幅度提高估计精度、测向效率，降低系统成本及计算复杂度的特点，可应用于雷达、声呐、无线通信系统及智能天线系统等领域中。

附图说明

图1是本发明提出的基于SULA和SLA构造的稀疏L阵的阵列结构；

图2是本发明所提稀疏L阵的x轴(子阵1)和z轴(子阵2)的阵元配置；

图3是利用本发明估计二维DOA的俯仰角和方位角的角度散布图；

图4是本发明的二维DOA估计与现有方式的性能随信噪比变化的对比图；

图5是本发明的二维DOA估计与现有方式的性能随采样快拍数变化的对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

步骤1：设置天线阵列

设置一个如图1所示的x-z平面的L形天线阵列，子阵1(subarray1)是一个阵元数为M＝6的SULA加上一个辅助阵元构成的稀疏线阵，因此M₁＝M+1＝7，子阵2(subarray2)是一个阵元数为M₂＝M＝6的SLA，由于原点处参考阵元是两个子阵的共有阵元，因此本案例所提稀疏L阵共含M₁+M₂-1＝12个阵元。本技术方案中设信号波长为λ＝0.8m，则子阵1中的SULA的阵元间距为d_x＝λ＝0.8m，辅助阵元和原点处的参考阵元间距为d＝λ/2＝0.4m；子阵2的阵元间距分别为：d₁＝λ＝0.8m、d₂＝λ/2＝0.4m、d₃＝3λ/2＝1.2m、d₄＝2λ＝1.6m、d₅＝λ＝0.8m。本实施方式有K＝2个窄带非相干信号以不同方向入射到此阵列，各阵元上的噪声为加性高斯白噪声，且噪声与信号不相关。

步骤2：求子阵2(SLA)的所有阵元的接收数据z(t)的自相关矩阵在实际工作中，用z(t)的N＝200次采样数据{z(1),z(2),...,z(200)}建立如下协方差矩阵：式中，t为采样的时间序号。

步骤3：确定的噪声子空间U_zn：

对步骤2所得协方差矩阵做特征值分解，得到特征值和其相应的特征向量，并利用其中最大的K＝2个特征值所对应的特征向量u_z1,u_z2构建矩阵U_zs＝{u_z1,u_z2}作为信号子空间，利用其余M-K＝4个小特征值对应的特征向量u_z3,u_z4,u_z5,u_z6构建矩阵U_zn＝{u_z3,u_z4,u_z5,u_z6}作为噪声子空间，2个大特征值组成对角矩阵D_s，其余4个小特征值组成对角矩阵D_n，如下式所示：

步骤4：求俯仰角

利用信号子空间和噪声子空间的正交性得到MUSIC算法的空间谱函数为并进行谱峰搜索处理：

第一次谱峰搜索时，搜索角θ的变化范围为：从0°以步长L₀＝1°增长到180°，P_MUSIC可以得到2个最大峰值，将这两个最大峰值对应的角度值分别记为θ′₁，θ′₂；然后设置第二次搜索的搜索步长为L₁＝0.01°，对应的两个搜索区间分别为：[θ′₁-2°,θ′₁+2°]和[θ′₂-2°,θ′₂+2°]，每个搜索区间得到一个最大峰值，即为最终的估计值和

步骤5：求信源自相关矩阵

根据估计值和可得相应的阵列流型矩阵的精确估计值其中，基于的特征分解所得到的U_zs、D_s可得信源自相关矩阵的近似值为：其中表示的M-P伪逆。

步骤6：求互相关矩阵

计算x(t)＝[x₁(t),x₂(t),...,x₇(t)]^T的N＝200次采样数据{x(1),x(2),...,x(200)}和子阵2除参考阵元外的其余阵元的接收数据z′(t)＝[z₂(t),z₃(t),...,z₆(t)]^T的200次采样数据{z′(1),z′(2),...,z′(200)}之间的互相关矩阵即：

步骤7：求x轴(子阵1)阵列流型矩阵：

由于步骤6中所得的互相关矩阵是不受加性高斯白噪声影响的矩阵，且理论上因此，关于A_xR_s的最大似然(ML)估计问题可以转化为如下最小二乘(LS)问题：其中表示的第2到M₂＝M＝6行的共轭转置。因此，利用步骤4、5得到的及步骤6得到的互相关矩阵，A_xR_s的估计值可以表示为：再结合步骤5得到的信源自相关矩阵，可得方位角的阵列流型矩阵为：由于是由和求得的，因此中方位角和中的俯仰角是一一对应的，从而后面可以对按列求取相应的方位角而不需要任何子空间分解。

步骤8：求方位角的粗估计值

考虑到子阵1是由一个阵元间距d_x＝λ的SULA和一个距原点处参考阵元间距d＝λ/2辅助阵元组成，且若单独用SULA来求解方位角时存在周期性的模糊，因此，先用辅助阵元和SULA的前两个阵元构成一个阵元间距等于半倍波长的传统ULA，然后利用ULA的平移不变性，利用步骤7所得的的前三行数据组成的矩阵按列求取相应的方位角的估计值，则第k个方位角的粗估计值φ_k′可表示为：其中和分别表示的第k列的前两行和最后两行的数据。

步骤9：求方位角的细估计值

由于所选SULA的阵元间距d_x＝λ，因此将0～180度的角度空间划分为3段： 在得到步骤8得到方位角的粗估计值φ_k′后，判断φ_k′属于上述三段中的哪一段，然后利用SULA的平移不变性，求相应的方位角的精确估计值，记步骤7所得的的除第二行以外的数据组成的矩阵其中，

表示的第k列，求取方位角的细估计值

其中，和分别表示a(φ_k″)的前5行和后5行的数据。

本方案中，当取信源方向(θ₁,φ₁)＝(50°,55°)，(θ₂,φ₂)＝(80°,65°)，M＝6，N＝200，独立重复实验次数P＝1000及信噪比SNR＝5dB时，可以得到相应的二维波达方向的估计值。图3是上述仿真环境下的角度散布图，由该图可知，俯仰角和方位角比较集中的分布于真实值附近，说明了本发明所提方案能够实现角度自动配对，且估计角度精度较高。图4(a)和4(b)分别是本方案所选的两个信号的俯仰角和方位角的均方根误差(RMSE)随信噪比SNR变化的仿真结果图，其中，(θ₁,φ₁)＝(50°,55°)，(θ₂,φ₂)＝(80°,65°)，M＝6，N＝200，P＝1000，SNR＝(0～25)dB，由图4可知，采用本专利提出的稀疏L阵和相应的二维DOA估计算法之后，俯仰角和方位角的性能都有很大的提升，比本方案所提对比算法(CCM-based和JSVD)至少有10dB的提升，且即使在低信噪比下，本发明所提算法也能得到一个较好的估计性能。图5(a)和5(b)分别是俯仰角和方位角的RMSE随采样快拍数N的变化曲线，其中，M＝6，N＝200～3700，SNR＝5dB，(θ₁,φ₁)＝(50,55)，(θ₂,φ₂)＝(80,65)，由图5可知，采用本专利提出的稀疏L阵和相应的二维DOA估计算法之后，俯仰角和方位角的性能均比对比算法提升至少5dB，且即使在少快拍数下，本发明所提算法也能获得较为精确的估计值。

因此，本发明所提的新的L阵列及其相应的二维DOA估计算法能够很好的提高二维DOA估计的测向精度，降低系统成本，并在一定程度上降低了计算复杂度。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。

Claims (2)

1.一种稀疏L阵列的二维波达方向的估计方法，其特征在于，包括下列步骤：

步骤1：设置稀疏L阵列：阵元数不同的线性第一子阵、第二子阵构成的L形阵列，两个线阵的共有阵元定义为参考阵元，第一子阵由阵元间距等于波长的稀疏均匀线阵和一个辅助阵元构成，所述辅助阵元到参考阵元的距离为半倍波长；第二子阵由最小阵元间距小于或等于半倍波长的任意稀疏线阵构成；且第一子阵的阵元数大于或等于4，第二子阵的阵元数大于或等于3；

设置L阵列的两个子阵与二维波达方向DOA的两个角的对应关系，对应第一子阵的定义为第一角度，对应第二子阵的定义为第二角度；

步骤2：L阵列接收K个不相关信源的入射信号，得到第一子阵、第二子阵的接收数据，其中K小于第二子阵的阵元数M2；

步骤3：计算第二子阵上所有阵元的接收数据在N次采样下的自相关矩阵并对做特征值分解，取的前K个最大特征值对应的特征向量为列构建特征向量矩阵U_zs作为信号子空间，剩余的特征值对应的特征向量为列构建特征向量矩阵U_zn作为噪声子空间；取的前K个最大特征值组成对角矩阵D_s，剩余的特征值组成对角矩阵D_n；

步骤4：基于信号子空间U_zs和噪声子空间U_zn得到MUSIC算法的空间普函数P_MUSIC：a(θ)表示关于搜索角度θ的方向矢量，其中θ∈[0°,180°]；对P_MUSIC进行普峰搜索，取前K个最大峰值所对应的搜索角度θ作为K个信源的第二角度估计值；

步骤5：由K个信源的第二角度估计值的方向矢量得到第二子阵的阵列流型矩阵再根据公式得到信源协方差矩阵的估计值其中符号(·)⁺表示M-P广义逆；

步骤6：计算第一子阵的所有阵元的接收数据和第二子阵除参考阵元外的所有阵元的接收数据在N次采样下的互相关矩阵

步骤7：根据公式计算第一子阵的阵列流型矩阵其中表示第二子阵的阵列流型矩阵的第2到M2行的共轭转置；

步骤8：基于第一子阵的阵列流型矩阵的前三行数据组成矩阵

根据公式计算各信源的第一角度粗估计值φ_k′，其中和分别表示的第k列的前两行和最后两行的数据；

步骤9：对各第一角度粗估计值φ_k′进行角度调整，得到第一角度精估计值

将第一子阵的阵列流型矩阵的第二行数据删除得到矩阵

将180°均分为3个不重叠区间：最小区间、中间区间、最大区间；

若φ_k′在最小区间，则

若φ_k′在中间区间，则

若φ_k′在最大区间，则

其中，和分别表示矩阵的第k列的前M1-2和最后M1-2行的数据，M1为第一子阵的阵元数。

2.如权利要求1所述的方法，其特征在于，步骤4中，采用递归网格划分的方式对P_MUSIC进行普峰搜索：

设置迭代搜索次数n,每次的搜索步长L_i，其中i＝0,1,…,n-1；

对于初始搜索，基于搜索步长L₀，在搜索范围[0°,180°]内，取前K个最大峰值所对应的θ，得到K个θ′_k，其中k＝1,…,K；

对于第1～n-1次搜索，基于当前搜索步长L_i和当前θ′_k的取值，在搜索范围θ′_k±2L_i-1内，取最大峰值所对应的角度作为当前θ′_k的更新值，直到n次迭代搜索结束并输出当前θ′_k作为K个信源的第二角度估计值。