一种基于平行稀疏阵列虚拟孔洞填充的二维DOA估计算法*

王嘉伟,杨赟秀,陈文东,舒 勤

(1.四川大学 电气工程学院,成都 610065;2.西南技术物理研究所,成都 610045;3.北京遥感设备研究所,北京 100005)

0 引 言

在阵列信号处理领域,波达角(Direction of Arrival,DOA)估计一直是一个重要的技术问题,在无线通信、雷达和声呐等许多领域中有着广泛的研究和应用[1]。通常情况下,阵元数目是影响估计效果的一个重要因素,但是由于成本和运算复杂度的制约,它无法被无限增加。而由稀疏阵列产生的虚拟阵列则在很大程度上突破了这一限制,它发端于压缩感知(Compressed Sensing,CS)理论[2-3],通过计算和处理接收信号的协方差矩阵或伪协方差矩阵,从数学层面实现了对阵列的拓展,从而能以更少的初始阵元数获得更大的阵列孔径与阵元自由度(Degrees-of-freedom,DOF),有着极为广泛的应用前景。为此,研究人员已经提出了许多基于稀疏阵列的高分辨率估计算法,包括SS-MUSIC(Spatial Smoothing MUSIC)[4]、SS-ESPRIT(Spatial Smoothing Estimation of Signal Parameters via Rotational Invariance techniques)[5]等,但它们都需要连续虚拟阵元去进行下一步处理(往往是空间平滑操作),这时随机分布在其上的孔洞就成为制约其连续自由度增长的关键。

当前,对于虚拟孔洞的研究往往集中于一维估计领域,主要分为前处理与后处理两个思路。前处理是指更加合理地设计初始阵列的排布,例如最小冗余阵列[6]、嵌套阵列[7]、改进的互质阵列[8]、CACIS[9]、CADiS[9]等,从而尽可能地减少虚拟孔洞的数目及其带来的影响。而后处理则是指在收到信号后,对于虚拟孔洞位置的空缺信息进行近似化填充,例如核范数最小化[10]、协方差矩阵的半正定结构[11]等,它们通过近似估计孔洞处应有的接收信号,将本该舍弃的虚拟阵元信息更加充分地利用起来。

在二维DOA估计领域,对于虚拟孔洞的研究则较少,往往是直接将孔洞之后的阵元信息舍去。这造成了一定的资源浪费,但若直接将一维孔洞填充思路拓展至二维,则会极大地增加运算复杂度,估计效果也不很理想。针对于此,本文提出一种基于二维DOA估计的虚拟孔洞填充算法,它结合平行稀疏阵列,将繁琐的二维估计问题转化为两个简易的一维估计问题,并通过提前计算其中一子阵列虚拟孔洞位置去设计另一子阵列的排布,从而实现对孔洞的完全填充。本文算法是一种典型的前处理思路,它提升了阵列资源利用率与阵列自由度,能够以更少的初始阵元获得更加精确的估计效果。

本文所用符号说明:(·)T,(·)*,(·)H,(·)-1和(·)†分别表示矩阵的转置、共轭、共轭转置、求逆和求伪逆;diag(·)表示将向量或对角阵对角化;vec(·)表示将矩阵向量化;E(·)表示求矩阵的统计期望;⊗和⊙分别表示Kronecker积和Khatri-Rao积;det(·)表示矩阵行列式;root(·)和arg(·)分别表示求根和求相位角;IK是维度为K×K的单位矩阵;d为基本阵元间距,它被设置为接收信号的半波长λ/2。

1 信号模型

本文的阵列模型采用平行稀疏阵列结构,由两个相互平行的稀疏线性阵列组成,如图1所示,包含子阵1和子阵2,它们分别有L1和L2个阵元。这里,分别用P和Q代表两子阵的阵元位置,表示为

图1 平行稀疏阵列几何模型[12-13]

(1)

(2)

式中:子阵1的方向矩阵A=[a1,a2,…,aK]的维度为L1×K,子阵2的方向矩阵B=[b1,b2,…,bK]的维度为L2×K,它们都只包含角度α信息;响应矢量ak=[efkp1d,efkp2d,…,efkpL1d]T,bk=[efkq1d,efkq2d,…,efkqL2d]T对应第k个入射信号(pld和qld分别为式(1)中表示的两子阵的阵元位置);C=diag([ejπcos β1,ejπcos β2,…,ejπcos βK])为K×K的对角矩阵;n1(t)和n2(t)为两不相关的零均值加性高斯白噪声,且都与信号不相关。

2 基于平行稀疏阵列虚拟孔洞填充的二维DOA估计算法

2.1 稀疏阵列排布结构

为了方便,下面依次介绍互质阵列[15](Co-prime Array,CA)、改进互质阵列[8](Augmented Co-prime Array,ACA)、CACIS[9]和CADiS[9]等4种会在文中用到的稀疏阵列:

PCA={mNd,0≤m≤M-1}∪{nMd,0≤n≤N-1},

(3)

PACA={mNd,0≤m≤2M-1}∪{nMd,0≤n≤N-1},

(4)

(5)

(6)

2.2 虚拟阵列构造

为了扩充初始阵列为虚拟阵列,需要计算子阵1接收信号的协方差矩阵[16]:

(7)

通过对式(7)进行特征值分解,可从中估计出噪声功率(阵元数目小于信号数目时不再适用)并可依此去除噪声影响[12],可得

(8)

接着,对上式做向量化处理,便可得到基于协方差矩阵的虚拟阵列

(9)

考虑到其中存在重复的元素,需要剔除重复元素,从而可得

(10)

式中:F=diag(Rs)。记差集U={uld,1≤l≤L3}={(pi-pj)d|pi,pj∈P},则有虚拟阵列方向矩阵U=[u1,u2,…,uK]和其响应矢量uk=[efku1d,efku2d,…,efkuL3d]T。

为了进一步处理,传统算法往往需要提取式(10)中的连续阵元,而将虚拟孔洞后的阵元信息直接舍弃,这造成了一定程度的资源浪费。

2.3 虚拟孔洞填充的二维DOA估计算法

首先,需要计算子阵2接收信号的协方差矩阵,并对其做类似于式(8)的去噪处理:

(11)

接着,对式(11)做类似式(9)的向量化处理可得

(12)

将其中重复元素剔除,可得

(13)

记差集G={gld,1≤l≤L4}={(qi-qj)d|qi,qj∈Q},则有虚拟阵列方向矩阵G=[g1,g2,…,gK]和其响应矢量gk=[efkg1d,efkg2d,…,efkgL4d]T。

通过比较、观察式(13)和式(10)可以发现,两者在形式上具有一致性,且都是只包含角度α这一个未知量,只是虚拟阵元位置由于两子阵排布的不同而不同。因此,只要能够合理设计两子阵的排布,便能够对虚拟孔洞进行完全填充,从而能够大幅提升阵列资源利用率。

如图2 所示,为了能够以更少的物理阵元获得更好的虚拟孔洞填充效果,需要提前计算基于子阵1获得的虚拟阵列中正半孔洞的位置H。因此,可将子阵2的阵元位置设定为Q={0,H},易知基于此阵列获得的差集G可以对子阵1得到的虚拟孔洞完全填充,从而提升了虚拟阵列的连续自由度与估计精确度。

图2 传统算法的虚拟阵列、物理子阵2以及新算法的虚拟阵列(阵列参数设置为(4,5)1)

显然,由式(10)可得其虚拟孔洞位置,并据此在式(13)结果中寻找该位置应得的接收信号,之后只需将该值填充到对应位置便可以获得一个拥有更大连续自由度的阵列接收信号:

(14)

式中:填充后的虚拟阵列方向矩阵H=[h1,h2,…,hK]的维度为L5×K,其响应矢量为hk=[efkh1d,efkh2d,…,efkhL5d]T。

由于式(14)是一个单快拍信号,为了获得其协方差矩阵,需采用空间平滑算法,即将其分割成W=(L5+1)/2个子阵列,并求取子阵列协方差矩阵的算术平均值:

(15)

(16)

基于式(16),可得角度α的谱峰搜索表达式[4]:

(17)

通过式(17),将二维估计转化为一维搜索,有效地降低了运算复杂度。为了进一步降低算法复杂度,本文采用求根MUSIC[17]的思路,也即用zg1取代H(αg1)中的ejπcos αg1,并可将上述问题转化为求取下式的K个根:

(18)

从而可得角度α的估计值

(19)

随后,为了获得角度β的估计值,并完成角度自动配对过程,需要计算互协方差矩阵。由于两子阵接收噪声不相关,故而相较于两子阵的协方差矩阵,该式没有噪声项。

(20)

对上式向量化处理可得

(21)

式中:J=diag(CRs)。移除其中的重复元素可得

(22)

式中:方向矩阵D=[d1,d2,…,dK]的估计结果可由式(19)获得,维度记作L6×K,其中包含着角度α的顺序信息,故而可自动实现角度匹配过程。

而J可通过求解如下最小化问题得到[18]:

(23)

式中:‖·‖F表示Frobenius范数。其结果可近似于求解

(24)

从而,可得角度β的估计值

βk=arccos(arg(J)/π),k=1,2,…,K。

(25)

通过联立式(19)和式(25),可得相互匹配的方位角和俯仰角,即

(26)

(27)

2.4 算法步骤

基于以上分析,下面给出本文算法的具体步骤:

Step1 选择合适的稀疏阵列作为子阵1,通过计算其差集G得到正半虚拟孔洞位置H,并据此设置子阵2的阵元位置为Q={0,H}。

Step2 根据实际应用,通过有限的快拍对两子阵列的协方差矩阵和两者的互协方差矩阵进行估计,即

(28)

(29)

(30)

式中:Y为接收信号的快拍数。

Step3 向量化处理子阵1和2的协方差矩阵,获得虚拟阵列1和2,并将后者的结果填充于前者的虚拟孔洞处,从而获得一个拥有更多连续自由度的虚拟阵列3。

3 实验分析

3.1 阵列自由度分析

采用上文提到的四种稀疏阵列,假定初始阵元数目相同,比较传统算法与本文算法可获得的连续自由度。经过分析,可得严格的数目通式如表1所示。

表1 阵列的连续自由度[9,19]

由表1可知,虚拟孔洞填充后的阵列连续自由度有了较大提升,大多数情况都能有1倍甚至更多的提升。分析可知,新算法下获得的虚拟阵列的连续自由度值已经达到了单条物理阵列能达到的理论极限值,因此本算法能够极大地提升阵元利用率,并由此获得了更为精确的二维DOA估计结果。

3.2 二维DOA估计精度分析

将本文算法与文献[12]和文献[16]中算法进行比较,两算法都采用了两个相互平行且相同的子阵,前者为ACA,后者为CACIS。而为了能够更加清晰地比较算法间的性能,本文算法分别采用上述两种阵列进行实验,且三种算法的子阵1阵元数目设置相同。

首先,选取多个特殊值,直观比较三种算法的连续自由度,结果如表2和表3所示。分析可知,新算法以更少的初始阵元获得了更高的连续自由度值,从而节约了阵元成本。

表2 本文算法和文献[12]算法的连续自由度比较

表3 本文算法和文献[16]算法的连续自由度比较

为比较三种算法的高分辨率性能,进行了单次估计实验。假设入射信号数为2,快拍数为500,信噪比为10 dB,假设文献[12]算法的阵列参数为(3,5)2,文献[16]算法的阵列参数为(6,5,3)3,并将本文算法的初始阵列分别设定为上述两阵列参数,实验结果如图3所示。可以发现,虽然有一定误差,但三种算法都能清晰地分辨出两个较为接近的角度,且此时新算法用到的初始阵元数目最少,相较另外两算法有着明显优势。

图3 DOA高分辨率性能比较

考虑到单次实验存在随机性,为了进一步比较三种算法的真实性能,下面设置蒙特卡洛实验,采用均方误差(Root Mean Square Error,RMSE)准则,定义为

(31)

图4 不同信噪比算法性能比较

接着,分析快拍数对两种算法的影响。假设文献[12]算法的阵列参数为(4,3)2,文献[16]算法的阵列参数为(6,5,2)3,显然两者的初始阵元数目都为20,信噪比为0 dB,信源数及其角度设置和上文信噪比实验中完全一致,结果如图5所示。观察可知,本文算法的估计精度在快拍数较小时,相比另外两算法有了较大的提升,这论证了本文算法在增加连续虚拟阵元数目上的显著性能。

图5 不同快拍数算法性能比较

最后,比较两种算法的运算复杂度。由文献[12]可知其算法复杂度为O{2LY+2K3+(2MN+2M)3},其中L为两子阵列阵元数目。由文献[16]可知其算法复杂度为O{2(M+N-1)2Y+(2M(N+1))3}。而由于初始阵列排布的不同,新算法下的子阵2的排布也不同,从而将直接影响本文算法复杂度,故这里不易于直接比较。为此,进行500次实验,快拍数为500,信源与上文蒙特卡洛实验一致,比较三种算法的运行时间,结果如表4所示。

表4 500次实验运行时间比较

可以发现,本文算法在所用运算时间上也明显比另外两算法少。而且由表1可知随着子阵1阵元数目的增加,本文算法在总阵元数目上相比另外两算法也会减少得更为显著,它显然也会直接反映为运行时间上的减少,也就是说随着阵元数目的增加,本文算法在节约运算成本上的能力会体现得更为显著。

4 结束语

本文针对采用稀疏阵列进行二维DOA估计时由于出现孔洞而需要舍弃一些信息,导致其性能不佳这一缺点,提出了一种低复杂度基于平行稀疏阵列虚拟孔洞填充的二维DOA估计算法。该算法最直观地优势在于能以更少的物理阵元数目获得更大的连续虚拟阵元数目,从而为后续的估计工作提供更多的有用信息,极大提升了阵元资源的利用率。同时,阵元数目的减少也在一定程度上降低了运算复杂度,节约了运算成本。另外,本文舍弃传统的二维角度谱峰搜索思路,将二维问题转化为两个一维问题,并借用Root-MUSIC算法思想进一步压缩了算法的运算复杂度。关键的是,本文算法推广性较好,对于几乎所有会产生虚拟孔洞的稀疏线性阵列,它都能在不同程度上提升对于物理阵元的资源利用率。

但是,其中一子阵阵元数目的减少,会影响用于角度β估计的两子阵互协方差矩阵在有限快拍下的估计精度,若单论该角的估计效果是不如同算法下设置两子阵阵元数目相同时的情况,这也是本算法的局限性,并将会在未来的研究中着重探讨。