王 宏 何培宇 喻伟闯 崔 敖 徐自励
(1.四川大学电子信息学院,四川成都 610065;2.中国民用航空局第二研究所,四川成都 610041)
在阵列信号处理领域中,波达角(Direction of Arrival,DOA)估计一直是一个重要的研究方向,在无线通信、语音信号、声纳、电子对抗等许多领域中有着广泛的研究与应用[1-2]。估计DOA时,常需要对信号进行二维的DOA估计[3],即需要获取信号入射的俯仰角和方位角。为此,学者已经提出了多种基于面阵结构的DOA估计方法,包括矩形阵[4]、L形阵列[5]和平行线阵[6]等。通常情况下,传统的平面阵列结构需要满足空间采样定理,即相邻两阵元之间的间距不能超过入射信号波长的一半,以保证DOA估计不出现模糊[7]。而在当前阵列信号处理研究中,与传统阵列相比,稀疏阵列能在相同物理阵元数目下能提供更大的阵列孔径以及更多可用的虚拟阵元,从而拥有更高的分辨力和自由度(Degrees of Freedom,DOF)[8]。
当前,稀疏阵列的研究大多集中在一维DOA 估计问题上,其阵列结构主要可分为最小冗余阵[9]、嵌套阵[10]和互质阵[11]。其中,最小冗余阵自由度最大,但其阵元位置和阵列自由度无闭合表达式,不易向高维扩展[9]。嵌套阵的阵元位置和自由度虽有闭合表达式,但密集阵元间互耦影响较大[12]。相较于最小冗余阵和嵌套阵,互质阵具有阵列结构简单、阵元位置和自由度有闭合表达式、阵元间互耦影响较小等优势,但该阵列生成的差分虚拟阵列存在“孔洞”缺失,导致其连续自由度较小[13]。如何将一维稀疏阵列有效地进行二维推广,是当前学界仍在研究的重要问题[14]。为此,文献[15]提出了一种基于嵌套阵的平行线阵结构,该阵列通过以嵌套阵阵列作为平行线阵的子阵,提升了传统平行线阵的自由度和分辨力,但是也导致其密集阵元间的互耦影响较大。文献[16]则以互质阵为基础,提出了一种互质平行线阵结构,该方法虽有效地避免了嵌套阵互耦影响较大这一缺点,但是由于其互质阵列的差分孔洞特性,即差分出的连续虚拟阵元数目较少,因此该方法的分辨力和自由度不能得到有效提升。此外,该方法使用传统的二维MUSIC 估计算法对信号进行估计,计算复杂度较高。文献[17]和文献[18]则针对文献[16]计算复杂度较高这一缺点,利用求根算法将二维DOA估计转换成两个一维DOA估计问题,有效地降低了其算法复杂度,但是其仍然未解决互质阵的差分孔洞特性所造成的性能损失问题。
针对以上问题,本文提出了一种基于广义互质双平行阵列的二维DOA估计方法。该方法利用了一种广义互质阵列结构,对传统平行线阵进行稀疏域的推广,具有更高的分辨力和自由度,又较好地改进了标准互质阵的差分虚拟阵列存在的“孔洞”缺失问题,提升了虚拟阵列的连续阵元个数。并且,本文方法是在传统平行线阵估计算法的基础上,提出的一种改进二维Root-MUSIC估计算法,该方法通过多项式求根,将二维DOA 估计问题转为两个一维的DOA 估计问题,既提高了算法的精度,又降低了复杂度。
本文所用符号说明:(.)T,(.)*,(.)H,分别表示矩阵的转置、共轭与共轭转置;diag(v)表示以向量v为对角元素的对角矩阵;E{.}表示求矩阵的统计期望,vec(.)表示将矩阵向量化;det{.}表示矩阵行列式;⊗表示Kronecker 积运算,I表示单位矩阵;arg(.)表示取复数矢量的相位角。
标准互质阵由两个均匀线阵构成,两个子阵位于同一水平线上,其中均匀线阵1的阵元个数为M,阵元间距为Nd,d=λ/2;均匀线阵2 的阵元个数为N,阵元间距为Md,其中M和N互质,且M<N。标准互质阵的阵列参数通常可表示为(M,N)。该阵列的阵元位置集SCA为
标准互质阵的阵列孔径为(MN-M)d[11]。根据式(1)能得到位置SCA的差集Ld,去掉差集中的冗余元素得到,但是由于互质阵对应的并非完全连续,导致互质阵的虚拟差分阵元位置中拥有较多的“孔洞”[11]。标准互质阵的差集只有在[-M-N-1,M+N-1]范围内是连续的[7],因此其自由度DOF为
广义互质阵是在标准互质阵列结构上发展出来的一种新的互质阵列结构,该阵列仍由两个子阵组成,具有M+N-1 个阵元,M和N仍然互质,但此时M<N这一条件不再要求成立[19]。与标准互质阵不同的是,广义互质阵引入了一个正整数的压缩因子p来对其子阵1 的阵元间距进行压缩,从而减小了这一均匀线阵之间的间距。两者的阵列结构如图1所示。
其中,p是2到M之间的一个整数,也是一个整数。根据式(3)可以确定,和N也满足互质关系,且其阵列位置集SGCA为
对式(4)集合进行自差分和互差分可得阵列的差分位置集合Ld为[19]
对差分位置集合Ld中的元素分析可知,广义互质阵的差分位置集合Ld在[-MN++1,MN--1]范围内是连续的[19]。广义互质阵的差分优化阵自由度DOF为
根据式(2)和式(6)可得,广义互质阵的虚拟差分优化阵的自由度和阵列孔径大于传统的标准互质阵模型,其性能随着的减小而增大,或随着压缩因子p的增大而增大,但过大的压缩因子会增加其阵元间的互耦影响[20]。
在阵列模型方面,本文以单个广义互质阵为基础,构建了一种广义互质平行阵列结构。如图2 所示,该阵列由两个完全相同且相互平行的广义互质子阵组成,其中子阵1 位于y轴上,子阵2 和子阵1相距d=λ/2,λ 为入射信号的波长。在该阵列结构中,每个子阵的阵元总数均为S=M+N-1,阵列参数均为(M,N)。从2 节可得,单个广义互质子阵又由两个均匀的直线阵构成。其中,均匀子阵ULA1 的阵元数为N,阵元间距为d1=,其中=Md/p,,p∈Z+。均匀子阵ULA 2的阵元数为M,阵元间距d2=Nd。
广义互质子阵1 的阵元位置集可表示为(0,zid),其中,zi∈S,i=1,2,…,M+N-1,且集和S 的结构如式(7)所示。类似的,广义互质子阵2 的阵元位置可表示为(d,zid)。
设有K个非相关的远场窄带信号sk(t)(k=1,2,…,K)分别从方位角θk和俯仰角φk入射到广义互质双平行线阵,且用αk和βk分别表示第k个入射信号与y轴和x轴的夹角,则有以下关系成立
整个阵列的接收信号模型可以表示为
其中,s(t)=[s1(t),s2(t),…,sK(t)]T∈CK×1是信号源矢量,n1(t)∈C(M+N-1)×1和n2(t)∈C(M+N-1)×1分别表示为两个与信号独立同分布的高斯白噪声。A1∈C(M+N-1)×K和A2∈C(M+N-1)×K分别表示为子阵1和子阵2的方向矢量矩阵,且分别为
根据式(9)可以得到子阵1在时刻t的接收信号协方差矩阵为
其中,Λ=E{s(t)sH(t)}=表示第k个信号的入射功率,表示噪声的功率,IM为M+N-1阶单位阵。
对协方差矩阵R11向量化得到
式中,B1=,⊗表示Kronecker 积运算,p=为信号功率矢量,e=,ei∈R(M+N-1)×1(i=1,…,M+N-1),矢量ei中除了第i个元素为1,其余元素均为0。
将矢量中的冗余元素进行去重,并取连续部分重排后得到
根据式(9),再次得到广义互质双平行阵之间的互相关矩阵为
同理,向量化互相关矩阵R21得到
式中,B2=[a2(αK,βK)],将矢量z2去除冗余,并取连续重排得到,并且将依Toeplitz重排后可以得到
如3.2 节所述,在获得虚拟化之后的互相关矩阵和自相关矩阵之后,按照式(18)构建一个增广协方差矩阵[15]
根据式(18)可以得到一个虚拟化之后的均匀双平行线线阵接收信号的协方差矩阵,相对广义互质双平行阵列模型本身,其虚拟化之后生成的平行差分优化阵列增加了较多的虚拟阵元。对该矩阵进行特征分解以后,即可得到式(19)中的信号子空间Us和噪声子空间Un,利用这两个空间的正交性,即可利用传统二维MUSIC 估计算法来获得空间谱函数,其公式如下
但是,在实际的系统应用中,由于传统算法中二维谱搜索具有较高的计算复杂度,因此在对时延要求较低的系统中该方案不具备可行性[7],可以通过改进的二维Root-MUSIC 算法最小化谱函数f(α,β),将其分解为两个一维的DOA 估计问题,从而对入射信号进行二维估计。
为此,首先将矩阵Un划分成两个子阵,如下所示
式中,ξ=[1ej2πdcos(β)/λ]T,b1(α)为平行线阵虚拟化后的子阵1连续部分所对应的导向矢量,其结构可表示为
式(23)有效地将二维的角度估计问题分解为两个一维的角度估计问题,即Q(α)只包含角度α的信息,矢量ξ只含有角度β的信息。令v=,将式(24)重新表示为b1(v)=[1,v,…,]T,则角度αk(k=1,2,...,K)通过求解多项式方程det{Q(v)}来得到,即
其中
其中,arg(.)表示求相位,vk是多项式det{Q(v)}与角度对应的解。接着将所求得的代入式(25)即可求得方向余弦估计值为
最终可以求得入射信号的方位角θk和俯仰角φk的估计值为
步骤1:对接收信号进行离散化采样,在有限的快拍数下得到广义互质子阵1 的估计协方差矩阵;对其按照式(12)和式(13)进行向量化,去除冗余元素,在进行重排得到矢量,利用所得的矢量按式(14)进行Toeplitz重排得到。
步骤2:同步骤1,得到广义互质子阵1 和广义互质子阵2 的估计互协方差矩阵;对其按式(12)和式(13)进行向量化,去除冗余元素,在进行重排得到矢量,利用所得的矢量按式(17)进行Toeplitz重排得到。
步骤3:根据式(18)构建增广协方差矩阵,对进行矩阵分解得到其噪声子空间Un,然后将噪声子空间Un按式(22)划分成两个子矩阵Un1和Un2,并按式(25)和式(26)计算得到多项式方程的det{Q(v)}的解。
步骤4:将步骤3中所得到的解代入式(26)和式(27)得到估计的,再将其估计结果代入式(29),即可估计出方位角和俯仰角。
本文算法的步骤主要包括协方差矩阵R11和互协方差矩阵R21的估计,Toeplitz矩阵重构,以及矩阵的特征分解,经推导,算法的复杂度约为O{2S2T+(2L)3},其中S=M+N-1,L=M(N+1),T为采样快拍数,M和N互为质数。
设各个阵列物理实际阵元数量S相等,如S=12,并设标准互质阵和广义互质阵参数相同,如M=6,N=7,广义互质阵的压缩因子p分别在2、3、6 的三种情况下,本文以三种双平行线阵的一个子阵为基础,绘制出了传统非稀疏平行线阵、标准互质阵和广义互质阵三种阵列的子阵实际阵元位置分布,从而进行三种阵列的自由度大小对比。
由图3 可见,广义互质阵的自由度明显高于传统非稀疏线阵和标准互质阵,且自由度随着压缩因子的增大而增大。其中,普通均匀平行线阵由于没有虚拟差分特性,阵列自由度与实际的物理阵元个数相等。而标准互质阵虽然利用了稀疏阵列可以生成虚拟阵元这一特性,一定程度上提升了阵列孔径,但是由于生成的差分虚拟阵列存在较多的“孔洞”,因此可以利用的连续阵元数目不多。广义互质阵由于压缩因子p的引入,使子阵1 的间距减小,从而在子阵做位置的差分时,增加了虚拟差分优化阵列中连续阵元的数量,提高了阵列的自由度。且从图中可以看出广义互质阵的压缩因子越大,差分位置集合的连续阵元数量也就越多,当达到最大时,虚拟差分阵列已经完全无“孔洞”,阵列的自由度达到最大。
设有K个非相关远场窄带信号入射到阵列,阵列单个子阵实际物理阵元个数均为12,稀疏阵的阵列参数为M=6,N=7。信噪比SNR=20dB,快拍数T=2000,广义互质平行阵列的压缩因子p=3。传统双平行线阵方法、文献[16]方法和本文方法的估计入射信号的个数K分别为10、13 和20,在此条件下验证三种方法的估计精度和自由度。
从图4 中可以看出,三种方法均可以对入射信号进行较为准确的二维DOA 估计。但是传统双平行线阵算法,受限于空间采样定律,因而其阵列的孔径较小,估计的精度相对也较低[21]。此外,传统双平行线阵算法可用阵元数等于实际的物理阵元数,因而其自由度小,可估计的信源也较少。而文献[16]算法利用了阵列孔径大的稀疏阵列,具有可以生成虚拟阵元这一特性,明显提升了信号估计的精度和增加了可估计信源的个数,但由于该方法生成的虚拟差分阵列存在较多“孔洞”,以及使用了传统的二维MUSIC 方法,在面对较多信源入射时,估计性能出现了明显的下降。而本文方法采用的广义互质阵由于阵列孔径更大,生成的连续虚拟阵元数目也更多,且利用多项式求根的方法替代了传统二维MUSIC 方法对信号进行估计,因而从图中来看,无论是从可估计信号的数量还是估计的精度,明显都优于前两者。
设有两个非相关远场窄带信号分别从二维角度(65°,40°)和(122°,120°)入射到阵列,即K=2,快拍数固定为T=2000,蒙特卡罗实验次数为200,阵列单个子阵实际物理阵元个数均为S=12,稀疏阵的阵列参数为M=6,N=7,压缩因子p=2,3,6。SNR 从-5 dB 到15 dB 均匀变化,其各个阵列的均方根误差RMSE 随信噪比的关系如图5。同时,固定其信噪比SNR=20 dB,快拍数从100到2000,其各个阵列的均方根误差(RMSE)随快拍数的关系如图6。文中,二维均方根误差的定义为
由图5和图6可知,三种算法的RMSE均随着信噪比和快拍数的增加而减小,但本文方法优于传统非稀疏平行线阵方法和文献[16]方法。还可看到,广义稀疏平行阵的RMSE 随着压缩因子p的增加而增大,这与实验1 和2 的结果相吻合。因而,在同等物理阵元的条件下,本文方法的估计性能优于传统非稀疏平行线阵方法和文献[16]方法。
本文针对标准互质平行线阵在进行DOA 估计时,因差分虚拟阵列“孔洞”过多,导致其性能不佳这一缺点,提出了一种基于广义互质平行阵列的二维DOA 估计方法。该方法将广义互质阵列的特性应用于平行阵列二维DOA估计中,即以降低一定的阵元间距为代价,生成了更多可利用的连续虚拟阵元,较好地改进了标准互质平行阵中“孔洞”过多的缺点,提升了其估计性能。本文算法通过一种改进的Root-MUSIC 算法将一个二维估计问题分解为两个一维的估计问题,既提升了估计的精度,又降低了算法复杂度。然而,本文方法性能的提升是以适当增加阵元间的互耦影响为代价的,且该影响随压缩因子的增大而增大,如何选择适当的压缩因子,使其在一定互耦影响范围内,保持较高的估计性能,在未来还有待进一步的研究。