一种稀疏l阵及其二维doa估计方法
【专利摘要】本发明公开了一种稀疏L阵及其二维DOA估计方法，属于无线移动通信技术领域。本发明的稀疏L阵列包括由阵元间距等于波长的稀疏均匀线阵和一个辅助阵元构成第一子阵、由最小阵元间距小于或等于半倍波长的任意稀疏线阵构成的第二子阵，两个线阵的共有阵元为参考阵元，辅助阵元到参考阵元的距离为半倍波长。在二维DOA估计时，首先基于第二子阵的接收数据计算其自相关矩阵，并对其进行特征分解后估计对应的第二角度，再基于其计算信源自相关矩阵；基于两个子阵接收数据的互相关矩阵、信源自相关矩阵得到第一子阵的阵流行矩阵，从而完成第一子阵所对应的第一角度的估计处理，得到二维DOA。本发明的复杂程度低、DOA估计的精确度高。
【专利说明】
一种稀疏L阵及其二维DOA估计方法
技术领域
[0001] 本发明属于无线移动通信技术领域，特别是涉及一种基于线阵构造的L阵列及其 二维波达方向（D0A)估计方法。
【背景技术】
[0002] 随着近年来无线通信技术的迅速发展，人们对通信业务量和通信质量的需求也越 来越大，以阵列信号处理技术为核心的空分多址技术已成为下一代移动通信的关键。
[0003] 现有的二维D0A估计大多是基于由阵元间距等于半倍波长的简化面阵构成的。其 中，L阵由于有更大的有效孔径、更小的运算量、更易实现、更强的方法适用性等优点得到了 广泛的关注和应用。近几十年来，人们已经做了很多的利用L阵估计2-D D0A的研究并提出 了大量的算法。主要分为两大类:需要额外配对的算法，如文献"Nizar Tayem and Hyuck M Kwon,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 two crosscorrelation matrices to achieve automatic pairing in 2-D angle estimation problems,Antennas and Wireless Propagation Letters,IEEE,vol.6, pp.553-556,2007.，简称JSVD"和"J.Gu，P.Wei,and H.-M.Tai, "DOA estimation using cross-correlation matrix，''in Phased Array Systems and Technology(ARRAY),2010 IEEE International Symposium on. IEEE,2〇10，pp.593_598."，简称CCM-based。但现有的 L阵大多是由常规ULA(均匀线阵)构成，存在测向精度和系统成本之间的矛盾，为了缓解该 矛盾，2015 年，文献"Jian-Feng Gu，Wei_Ping Zhu，and MNS Swamy，Joint 2_d doa estimation via sparse 1-shaped array,Signal Processing,IEEE Transactions on, vol. 63，no. 5，pp. 1171 -1182，2015提出了基于SLA(稀疏线阵)和ULA的稀疏L阵，该算法和 传统由常规ULA构成的L阵相比虽然在性能上有所提升，但是该阵列并没有充分利用稀疏阵 列的优势，此外在利用SLA求解方位角时的性能不稳定，不能达到理想所需要的结果。
【发明内容】
[0004] 本发明的发明目的在于:为克服现有的稀疏L阵和传统L阵在估计算法复杂度，估 计精度方面的不足，提出一种结构简单的稀疏L阵及其相应的二维波达角(D0A)估计算法， 以达到降低计算复杂度、系统成本、简化处理程序、有效提高估计精度等目的。
[0005] 本发明为了利用稀疏阵有更大的阵列孔径，均匀线阵的平移不变性等特性而提出 一种基于SLA和SULA的L阵结构，然后利用互相关矩阵不受噪声影响的特性、ULA的平移不变 性、递归思想和最小二乘(LS)技术求解二维D0A。
[0006] 本发明的稀疏L阵列，包括阵元数不同的线性第一子阵(Ml个阵元）、第二子阵(M2 个阵元)构成的L形阵列，两个均匀线阵的共有阵元定义为参考阵元，第一子阵由阵元间距 等于波长的稀疏均匀线阵和一个辅助阵元构成，辅助阵元到参考阵元的距离为半倍波长， 第二子阵由最小阵元间距小于或等于半倍波长的任意稀疏线阵构成，且施多4，112多3。本发 明的稀疏L阵列（以下简称L阵)可以是位于x-z平面或y-z平面，第一子阵、第二子阵对应二 维DOA的两个角:方位角、俯仰角，第一子阵可对应方位角或俯仰角，第二子阵可对应俯仰角 或方位角，取决于L阵的放置方式。
[0007] 在基于本发明的L阵进行二维D0A求解时，先基于第二子阵的接收数据计算其自相 关矩阵，并对该自相关矩阵进行特征分解后估计对应的方位角或俯仰角，再基于其得到信 源自相关矩阵;基于两个子阵接收数据的互相关矩阵、信源自相关矩阵得到第一子阵的阵 流行矩阵，从而完成第一子阵所对应的俯仰角或方位角的估计处理，得到二维D0A，从而大 幅度降低处理量及处理的复杂程度、有效提高D0A估计的精确度，从而实现本发明目的。用 于本发明的稀疏L阵列的二维波达方向的估计方法具体包括下列步骤：
[0008] 步骤1:设置天线阵列并建立系统模型：
[0009] 设置L阵列的两个子阵与二维波达方向D0A的两个角的对应关系，对应第一子阵的 定义为第一角度，对应第二子阵的定义为第二角度，例如第一子阵放于x轴(定义方位角）、 第二子阵放于z轴(定义为俯仰角），则第一角度对应方位角，第二子阵对应俯仰角。为了便 于描述，以下用x-z平面的稀疏L阵对本发明进行说明。
[0010] 在t时刻第一子阵和第二子阵的接收数据分别为：
，其中， x(〇[-M0,…>.rw (or，…分别为第一子阵和第二子阵的接收数据矢 量，s⑴= [S1(t)，…，SK(t)]T为信号矢量，即信号源，AfEaxMO,…，a x((i)K)]为第一子阵 的阵列流型矩阵，
表示Ax的第k列导向矢 量，巾k表示第k个信源的第一角度，在x-z平面中，因 x轴对应方位角，z轴对应俯仰角，巾!^表 示第k个信源的方位角，A表示信源波长，辅助阵元与参考阵元的间距d = 0.5A，其他阵元间 距dx =入，e表示自然底数，j表示虚数单位。^=[82(0:)，…，az(0K)]为第二子阵的阵列流型 矩阵
.表示Az的第k列导向矢量，k = 1，…，K，0k表示第k个信源的第二角度，在x-z平面中，0k表示第k个信源的俯仰角，cU表示第 二子阵的第i个阵元与第(i_l)个阵元之间的间距，i = 〇,…，M2-1。
[0011] 步骤2:计算第二子阵(SLA)的所有阵元的接收数据玢)唯,(〇,⑴]y的自相关 矩阵利用时间平均代替统计平均，可以求得第二子阵的接收数据的自相关矩阵：
[0012] 步骤3:确定：^=的噪声子空间Uzn:对步骤2所得协方差矩阵氧2.进行特征值分解，取 的前K个最大特征值对应的特征向量为列构建特征向量矩阵U zs作为信号子空间，剩余的特 征值(Ms-K个)对应的特征向量为列构建特征向量矩阵Uzn作为噪声子空间;取的前K个最大特 征值组成对角矩阵Ds，剩余的特征值组成对角矩阵Dn，
[0013]步骤4:求第二角度:利用信号子空间和噪声子空间的正交性确定MUSIC算法的空
，其中a(0)表示关于搜索角度0的方向矢量，0G[0°，
180° ] 例如采用递归网格划分的方式 〇 来进行谱峰搜索：
[0014] 1)第一次谱峰搜索时，先对空间谱[0°，180]°进行一个粗略的划分，0以步长1^从 0°增长到180°，遍历搜索得到Pmusk的前K个最大峰值，这K个最大峰值对应的0值就是所求的 K个信号的俯仰角的大概的估计值0i/ ;
[0015] 2)在上一步所得K个0!/附近用一个加密的网格即更小0的增长步长，进行谱峰搜 索，得到K个更精确的俯仰角估计值；
[0016] 3)重复第2步，直到空间搜索的网格足够精细。
[0017] 步骤5:求信源自相关矩阵A :首先利用步骤3所得的Uzs及Ds得到UzsDsllf，然后利用步 骤4所得的俯仰角的估计值4,4，…,4求其相应的阵列流型矩阵I = a(之)]，最后 由UZSDSUZHS和求得信源协方差矩阵的估计值1\，从而建立俯仰角和方位角之间一一对应 的关系。
[0018] 步骤6:求互相关矩阵众、.__.:利用时间平均代替统计平均，求x轴全部阵元的接收数 据x(t)和z轴除参考阵元以外的阵元的接收数据^ (t)之间的互相关矩阵或z，，即：
[0019]步骤7:求x轴阵列流型矩阵:利用互相关矩阵Rw不受噪声影响的特性，将最大似 然(ML)估计转化为最小二乘(LS)问题，得到AXR^估计值，然后利用步骤4所得的z轴阵列流 型矩阵和步骤5所得的信源自相关矩阵束求x轴阵列流型矩阵的估计值I，即 人T二XXk1，其中AXRS的估计值IX = RT,.(人f (2:6,〇)+，符号"A(a: b，：）"表示对应矩阵A的 第a至第b行，符号（?）+表示M-P广义逆。
[0020] 步骤8:求方位角的粗估计值:利用步骤7所得的x轴第一子阵的阵列流型矩阵 的前三行数据组成的矩阵^4及1]1^的平移不变性求得方位角的粗估计值，其中第k个方位角 的粗估计值巾i/可根据公式
，K求得， 和 2分别表示矩阵的第k列的前两行和最后两行的数据。
[0021] 步骤9:求方位角的细估计值:将第一子阵的阵列流型矩阵人,的第二行数据删除 得到矩阵人：：;将180°均分为3个不重叠区间：最小区间、中间区间、最大区间；
[0022] 若巾1/在最小区间，
[0025] 其中，,4分别表示矩阵Ij：的第k列的前M1-2行和最后M1-2行的数据。
[0026] 本发明首先利用第二子阵的接收数据求得相应的自相关矩阵，接着利用一次特征 值分解，得到相应的噪声子空间Uzn，信号子空间U zs及信号的能量矩阵Ds，再基于MUSIC算法 得到高精度的俯仰角（或方位角）的估计值，从而很好的利用了 SLA的稀疏性带来的阵列孔 径的优势，且有效的降低了 MUSIC谱峰搜索的计算量;然后利用所得的俯仰角(或方位角)得 到相应的高精度的阵列流型矩阵又2，为后面估计方位角(或俯仰角）提供更精确的参数，再 利用特征值分解得到的信号子空间和信号能量矩阵及1得到信号自相关矩阵的估计值 食,，从而很大程度上提高了 的准确性；再然后利用不受加性噪声影响的互相关矩阵 食V，及焱，求的第一子阵的阵列流型矩阵的估计值人，；最后利用辅助阵元及ULA的平移 不变性先后估计方位角的粗略值和精确值，既解决了 SULA带来的周期性模糊问题，提高了 估计精度，又降低了计算复杂度，提高了测向的效率。因而，本发明具有方法简单、可大幅度 提高估计精度、测向效率，降低系统成本及计算复杂度的特点，可应用于雷达、声呐、无线通 信系统及智能天线系统等领域中。
【附图说明】
[0027]图1是本发明提出的基于SULA和SLA构造的稀疏L阵的阵列结构；
[0028] 图2是本发明所提稀疏L阵的x轴(子阵1)和z轴(子阵2)的阵元配置；
[0029] 图3是利用本发明估计二维D0A的俯仰角和方位角的角度散布图；
[0030] 图4是本发明的二维D0A估计与现有方式的性能随信噪比变化的对比图；
[0031 ]图5是本发明的二维D0A估计与现有方式的性能随采样快拍数变化的对比图。
【具体实施方式】
[0032]为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发 明作进一步地详细描述。
[0033]步骤1:设置天线阵列
[0034] 设置一个如图1所示的x-z平面的L形天线阵列，子阵l(subarrayl)是一个阵元数 为M = 6的SULA加上一个辅助阵元构成的稀疏线阵，因此Mi=M+l = 7,子阵2(subarray2)是 一个阵元数为M2 = M=6的SLA，由于原点处参考阵元是两个子阵的共有阵元，因此本案例所 提稀疏L阵共含MdMrl = 12个阵元。本技术方案中设信号波长为入=0.8m，则子阵1中的SULA 的阵元间距为dx = A = 〇 .8m，辅助阵元和原点处的参考阵元间距为(1 =人/2 = 0.4111;子阵2的 阵元间距分别为：di =人=0 ? 8m、cb = V2 = 0 ? 4m、cb = 3入/2 = 1 ? 2m、cU = 2入=1 ? 6m、d5 =入= 0.8m。本实施方式有K = 2个窄带非相干信号以不同方向入射到此阵列，各阵元上的噪声为 加性高斯白噪声，且噪声与信号不相关。
[0035] 步骤2:求子阵2(SLA)的所有阵元的接收数据z(t)的自相关矩阵歲^在实际工作 中，用z(t)的N = 200次采样数据{z(l)，z(2)，. . .，z(200)}建立如下协方差矩阵：
.式中，t为采样的时间序号。
[0036] 步骤3:确定之的噪声子空间Uzn:
[0037] 对步骤2所得协方差矩阵U故特征值分解，得到特征值和其相应的特征向量，并 利用其中最大的K = 2个特征值所对应的特征向量uzl，uz2构建矩阵Uzs={u zl，uz2}作为信号 子空间，利用其余M-K = 4个小特征值对应的特征向量UZ3，UZ4，UZ5，UZ6构建矩阵U zn= {UZ3，UZ4， uZ5，uz6}作为噪声子空间，2个大特征值组成对角矩阵Ds，其余4个小特征值组成对角矩阵Dn， 如下式所示
[0038] 步骤4:求俯仰角
[0039] 利用信号子空间和噪声子空间的正交性得到MU SI C算法的空间谱函数为
，并进行谱峰搜索处理：
[0040] 第一次谱峰搜索时，搜索角9的变化范围为:从0°以步长Lo = l°增长到180°，PmlisiC 可以得到2个最大峰值，将这两个最大峰值对应的角度值分别记为9/，如'；然后设置第二次 搜索的搜索步长为1^ = 0.01°，对应的两个搜索区间分别为：[01/-2°，01/+2°]和[0 2/-2°，02/ +2° ]，每个搜索区间得到一个最大峰值，即为最终的估计值$和堯。
[0041 ]步骤5:求信源自相关矩阵
[0042]根据估计值4和4可得相应的阵列流型矩阵的精确估计值人，，其 中：
+= 1，2，基于J^的特征分解所得到 的uzs、ds可得信源自相关矩阵的近似值为我=(Af)+，其中人：(人f A1 Af表 示又2.的M-P伪逆。
[0043] 步骤6:求互相关矩阵I,
[0044] 计算x(t) = [xi(t) ,X2(t), ? ? ? ,X7(t)]T的N = 200次采样数据{x(l) ,x(2), ? ? ?，x (200)}和子阵2除参考阵元外的其余阵元的接收数据z' (t) = [z2(t)，z3(t)，. . .，z6(t)]T的 200次采样数据{，（ 1)(2)，. . .，Z' ( 200 ) }之间的互相关矩阵I，即：
[0045] 步骤7:求x轴(子阵1)阵列流型矩阵：
[0046] 由于步骤6中所得的互相关矩阵是不受加性高斯白噪声影响的矩阵，且理论上 =￡^(/)*2"/(/)}=八，，( /(2:6,〇，因此，关于厶几的最大似然(1^)估计问题可以转化 为如下最小二乘（LS)问题：
，其中Af (2:6，:)表示 Az的第2到M2=M = 6行的共辄转置。因此，利用步骤4、5得到的&及步骤6得到的互相关矩 阵，AXRS的估计值可以表示为:乂瓦=R。.(又f (2: 6，:))+;再结合步骤5得到的信源自相关矩阵， 可得方位角的阵列流型矩阵为 :又,=.4)飞、1^。由于&是由I和良求得的，因此&中方位 角和又z中的俯仰角是一一对应的，从而后面可以对按列求取相应的方位角而不需要任 何子空间分解。
[0047]步骤8:求方位角的粗估计值
[0048] 考虑到子阵1是由一个阵元间距dx = A的SULA和一个距原点处参考阵元间距d =入/ 2辅助阵元组成，且若单独用SULA来求解方位角时存在周期性的模糊，因此，先用辅助阵元 和SULA的前两个阵元构成一个阵元间距等于半倍波长的传统ULA，然后利用ULA的平移不变 性，利用步骤7所得的的前三行数据组成的矩阵Al，按列求取相应的方位角的估计值，则 第k个方位角的粗估计值巾!/可表示为
> & = 1，2， 其中和<&2分别表示欠的第k列的前两行和最后两行的数据。
[0049]步骤9:求方位角的细估计值
[0050]由于所选SULA的阵元间距dx = A，因此将〇~180度的角度空间划分为3段：
;在得到步骤8得到方位角的粗估计值巾!/后，判断巾i/属于上述三段中的 明卜段，然后利用SULA的平移不变性，求相应的方位角的精确估计值，记步骤7所得的又、的除第 二行以外的数据组成的矩阵1〖42(<%-，&(成.)].其中，3(4〃) = [1，4-'"、,以::^^7表 示的第k列，求取方位角的细估计值& :
[0052] 其中，和心.,.4分别表不a( <H〃）的前5行和后5行的数据。
[0053] 本方案中，当取信源方向（01;，55°)，（02,伞2) = (80°，65°)，M = 6，N = 200,独立重复实验次数P= 1000及信噪比SNR = 5dB时，可以得到相应的二维波达方向的估 计值。图3是上述仿真环境下的角度散布图，由该图可知，俯仰角和方位角比较集中的分布 于真实值附近，说明了本发明所提方案能够实现角度自动配对，且估计角度精度较高。图4 (a)和4(b)分别是本方案所选的两个信号的俯仰角和方位角的均方根误差(RMSE)随信噪比 SNR 变化的仿真结果图，其中，（0:，，55°)，（02,巾2) = (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， ~=2〇〇~37〇〇,3冊=5(18，（0 1，巾1) = (5〇,55)，（02，巾2) = (8〇,65)，由图5可知，采用本专利 提出的稀疏L阵和相应的二维D0A估计算法之后，俯仰角和方位角的性能均比对比算法提升 至少5dB，且即使在少快拍数下，本发明所提算法也能获得较为精确的估计值。
[0054]因此，本发明所提的新的L阵列及其相应的二维D0A估计算法能够很好的提高二维 D0A估计的测向精度，降低系统成本，并在一定程度上降低了计算复杂度。
[0055]以上所述，仅为本发明的【具体实施方式】，本说明书中所公开的任一特征，除非特别 叙述，均可被其他等效或具有类似目的的替代特征加以替换;所公开的所有特征、或所有方 法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。
【主权项】
1. 一种稀疏L阵列，包括阵元数不同的线性第一子阵、第二子阵构成的L形阵列，两个线 阵的共有阵元定义为参考阵元，其特征在于，第一子阵由阵元间距等于波长的稀疏均匀线 阵和一个辅助阵元构成，所述辅助阵元到参考阵元的距离为半倍波长;第二子阵由最小阵 元间距小于或等于半倍波长的任意稀疏线阵构成;且第一子阵的阵元数大于或等于4,第二 子阵的阵元数大于或等于3。2. 如权利要求1所述的L阵列，将所述L形阵列的第一子阵、第二子阵分别放置于x轴和z 轴，x轴对应方位角，z轴对应俯仰角。3. 如权利要求1所述的L阵列，将所述L形阵列的第一子阵、第二子阵分别放置于y轴和z 轴，y轴对应方位角，z轴对应俯仰角。4. 一种用于权利要求1所述的稀疏L阵列的二维波达方向的估计方法，其特征在于，包 括下列步骤： 步骤1:设置L阵列的两个子阵与二维波达方向DOA的两个角的对应关系，对应第一子阵 的定义为第一角度，对应第二子阵的定义为第二角度； 步骤2:L阵列接收K个不相关信源的入射信号，得到第一子阵、第二子阵的接收数据，其 中K小于第二子阵的陈元数M2; 步骤3:计算第二子阵上所有阵元的接收数据在N次采样下的自相关矩阵食s ,并对：^2做 特征值分解，取的前K个最大特征值对应的特征向量为列构建特征向量矩阵Uzs作为信号 子空间，剩余的特征值对应的特征向量为列构建特征向量矩阵U zn作为噪声子空间；取&的 前K个最大特征值组成对角矩阵Ds，剩余的特征值组成对角矩阵D n; 步骤4:基于信号子空间Uzs和噪声子空间Uzn得到MUSIC算法的空间普函数PMUSIC::，a(0)表不关于搜索角度0的方向矢量，其中0 G [〇°，180° ];对Pmusic 进行普峰搜索，取前K个最大峰值所对应的搜索角度0作为K个信源的第二角度估计值； 步骤5:由K个信源的第二角度估计值的方向矢量得到第二子阵的阵列流型矩阵，再 根据公式兔=人(人,)+得到信源协方差矩阵的估计值％，其中符号（?）+表示M-P 广义逆； 步骤6:计算第一子阵的所有阵元的接收数据和第二子阵除参考陈元外的所有阵元的 接收数据在N次采样下的互相关矩阵食; 步骤7 :根据公式I = RTy (人f (2: M2,:))+ #计算第一子阵的阵列流型矩阵4，其中 if (2:M2,〇表示第二子阵的阵列流型矩阵I的第2到M2行的共辄转置； 步骤8:基于第一子阵的阵列流型矩阵的前三行数据组成矩阵又:；计算各信源的第一角度粗估 计值Vk，其中<4.;1和（^2分别表示&的第k列的前两行和最后两行的数据； 步骤9:对各第一角度粗估计值巾'k进行角度调整，得到第一角度精估计值& : 将第一子阵的阵列流型矩阵的第二行数据删除得到矩阵又::; 将180°均分为3个不重叠区间:最小区间、中间区间、最大区间；其中，仏_和％,4分别表示矩阵乂的第k列的前M1-2和最后M1-2行的数据，Ml为第一 子阵的陈元数。5.如权利要求4所述的方法，其特征在于，步骤4中，采用递归网格划分的方式对PMUSIC进 行普峰搜索： 设置迭代搜索次数n，每次的搜索步长U，其中i = 0，1，…，n-1; 对于初始搜索，基于搜索步长Lo,在搜索范围[0°，180° ]内，取前K个最大峰值所对应的 9,得到K个0\，其中k=l，…，K; 对于第1~n-1次搜索，基于当前搜索步长U和当前的取值，在搜索范围0/k±2L1- 1 内，取最大峰值所对应的角度作为当前9\的更新值，直到n次迭代搜索结束并输出当前 作为K个信源的第二角度估计值。
【文档编号】G01S3/14GK106054123SQ201610404072
【公开日】2016年10月26日
【申请日】2016年6月6日
【发明人】杨雨轩, 郑植, 杨姣, 杨海芬, 闫波, 孟会鹏
【申请人】电子科技大学