基于差集表遍历搜索的互素阵列DOA估计器

黄翔东,念天磊,马 欣

(天津大学电气自动化与信息工程学院,天津 300072)

0 引 言

阵列波达方向(direction-of-arrival,DOA)[1]估计是信号处理领域的重要研究方向,其应用十分广泛。传统方法一般采用阵元间距符合Nyquist采样定理的阵列进行接收如均匀线性阵列(uniform linear array,ULA),若要增加分辨率,就需通过增加阵元数量来扩大阵列孔径,但这会导致硬件成本的提高。另外,Nyquist采样定理要求阵列间距不大于半波长,故随着工作频率的升高,信号波长变短,阵元间距也随之变短,增加阵元数量还会引起阵元间的互耦加大而降低接收性能[2]。因此,在不满足Nyquist采样定理情况下,仅依靠稀疏阵列[2-7]估计出多目标的DOA成为急需突破的课题。很显然,该课题绕不开两个基本问题:一是选择稀疏阵列结构的布置,二是DOA估计算法的设计。

在选择稀疏阵列结构的布置方面,现有的稀疏阵列,例如嵌套阵列[3],最小冗余阵列(minimum redundancy array,MRA)[4-5],最小孔洞阵列[6]等稀疏非均匀阵列均可用于DOA估计,并且在与均匀线阵阵元数目相同的情况下可以获得更大的阵列孔径,从而提高DOA分辨率。但这些阵列各有缺陷:最小孔洞阵列和最小冗余阵列的阵列结构无法用具体的闭式表达式描述,这会给优化设计带来很大困难,其工程应用受到很大限制;而用嵌套阵列做多源目标估计时,虽然可用N个阵元实现O(N2)的自由度[2](即提高了识别分辨率),但其缺陷在于:嵌套结构使得稀疏阵列的观测协方差矩阵不易转化成Nyquist协方差矩阵(以方便后续子空间分解算法做源估计),造成目标估计的困难。

在DOA算法设计方面,主要分为压缩感知算法[5,8-15]和子空间分解算法两大类。对于第1类,需要根据阵列结构构造压缩感知系统模型,当为密集阵列时,直接以阵列流型矩阵作为观测矩阵,以阵元数据做观测向量,而构造出压缩感知模型[8-9];当为稀疏阵列时,则需对阵列流型矩阵的各列做克罗内克积构造观测矩阵,对由阵元数据算出的协方差矩阵做向量化构成观测向量,而构造出压缩感知模型[5]。但无论阵列是密集还是稀疏情况,在构造观测矩阵时,均需要构建候选波达方向足够丰富的过完备冗余字典(以保证真实DOA在该字典中呈现稀疏分布),进而以该字典元素为变量构造目标函数,通过匹配追踪(matching pursuit,MP)[10]、基追踪(basis pursuit,BP)[11]、正交MP(orthogonal MP,OMP)[12]、压缩采样MP(compressive sampling MP,CoSaMP)[13]等追踪方法或1-奇异值分解(singular value decomposition,SVD)[14]、最小绝对值收敛和选择算子 (least absolute shrinkage and selection operator,LASSO)[15]等手段对目标函数寻最优,而获得DOA估计。很明显,冗余字典的构造耗费很大的存储量,而目标函数的求解涉及复杂的递归迭代过程,计算复杂度高,不利于实时应用。对于第2类子空间分解法,通常采用主流的、技术成熟的多信号分类(multiple signal classification,MUSIC)[16]方式,故其计算复杂度比压缩感知要低,具有更高的工程应用价值,该方法关键在于:如何把稀疏阵列的协方差矩阵转换为Nyquist虚拟阵列的协方差矩阵。

近年来,基于互素欠采样[7,17-21]的谱分析成为信号处理学科的研究热点。相比于传统Nyquist样本分布,互素欠采样更为稀疏,这意味着在只需布置同样数目的阵元即可获得更大的阵列孔径,故有利于提升空间谱分辨率;且互素欠采样有中国余数定理[17]支撑,故易于推导出闭合形式的理论分析结果,这是嵌套阵列、最小冗余阵列、最小孔洞阵列远不能及的。文献[22]利用矩阵填充思想对稀疏协方差矩阵进行扩维构造一个更大维度的Toeplitz矩阵(空位补零),然后利用不定点延续算法对零元素进行有效替换,转化为完整的协方差矩阵。此算法可获得很大的阵列孔径,从而获得很高的自由度,但该算法除了需做MUSIC分解外,在矩阵填充过程中,还消耗了一系列较复杂的最优化操作(以保证零元素替换后矩阵的秩尽可能低)。

本文选定互素稀疏阵列实现DOA估计,笔者在前期工作中,已经把互素采样应用于时间序列谱分析[17-19]和DOA估计[20-21],但文献[20-21]是通过简单的频谱校正途径来做DOA估计,该方法不能估计多源目标。为实现多目标DOA估计,本文提出差集表遍历搜索方法,借助该方法,可以很容易地将互素稀疏阵列的协方差矩阵转换为Nyquist虚拟阵列的协方差矩阵,其扩张的阵列孔径尺寸虽然不如文献[22],但该方法可直接实现从观测协方差矩阵到虚拟协方差矩阵的解析映射,且无需矩阵填充的优化过程,除MUSIC分解外,其他附加的计算量可忽略不计。故出于计算复杂度与分辨率的权衡,本文估计器具有较高实用价值。

1 系统模型

令Nyquist采样间距为d=λ/2,图1给出的互素阵列由两个均匀线性子阵组成[23],其阵元间距分别为Md和Nd(M、N为互素正整数),阵元个数分别为N和2M,由于两个子阵的第1个阵元的空间坐标重合,故整个互素阵列的阵元个数为L=N+2M-1,其归一化(d为归一化因子)后的坐标集合S可闭式表达为

S= {mM,0≤m≤N-1}∪

{nN,0≤n≤2M-1}

(1)

图1 互素阵列模型Fig.1 Coprime array model

式(1)和图1表明,互素阵列的阵元间距允许远大于Nyquist间距d=λ/2,故具有很高的稀疏度(表现为阵元间存在较多的“空洞”)。

令式(1)的升序排列的坐标集合为

P=sort(S)={p0,pl,…,pN+2M-2}

(2)

以第1个阵元坐标p0=0为参考,假设远场有D个不相关的窄带信号s1,s2,…,sD,分别以到达角θ1,θ2,…,θD入射到接收阵列中,则互素阵列接收到的信号模型为

(3)

其中

a(θi)=[1,e-jπp1sin(θi),…,e-jπpN+2M-2sin(θi)]T

(4)

2 基于差集表遍历搜索的互素阵列DOA估计器

2.1 互素矩阵的差集表构造

对于阵元位置直接按Nyquist间距设定情况,只需把L×1的接收阵元向量x与其转置乘积后,做统计平均(假定K为快拍数量),即可获得子空间分解所需的统计协方差矩阵,即

(5)

当信号源不相干时,有rank(Rxx)=L,可直接做MUSIC分解识别L-1个源。令阵元快拍的空间协方差函数为φ(l),当快拍充足时,Rxx近似为Toeplitz阵。理想情况下,Rxx内部元素表示为

Rxx(u,v)=φ(u-v),0≤u,v≤L-1

(6)

进而,对于某个φ(l)值(l=-L+1,…,0,…,L-1)在矩阵Rxx对应的坐标(u,v)满足

l=u-v

(7)

从式(5)～式(7)的理想定义可看出,经典MUSIC分解要求协方差矩阵的元素行、列坐标的差值正好为阵元归一化坐标的差值,这决定了φ(l)值分布在与矩阵主对角线平行的线上。

然而,互素矩阵与该情况不同,由于阵元稀疏分布,在矩阵Rxx内,与φ(l)所对应的Rxx内的坐标(u,v)不再满足式(6)、式(7)的Toeplitz关系。因而,问题的关键在于:如何将不符合Toeplitz关系的Rxx转换为符合Toeplitz关系的Nyquist虚拟阵列(即将图1互素阵列的“空洞”填满、全部阵元间距为d=λ/2的假想阵列)的协方差矩阵。

具体而言,从图1可知:互素协方差矩阵Rxx内的任一元素Rxx(u,v),一定与互素阵列的两个阵元坐标pu、pv相对应,其对应关系描述为

Rxx(u,v)=φ(pu-pv),0≤u,v≤L-1

(8)

另外,由于pu或pv只能取自子阵1或子阵2,即

(9)

联立式(8)、式(9),可推知Rxx(u,v)与子阵1坐标m,m′、子阵2的坐标n,n′的对应关系,不外乎4种情况,即

(10)

式中,m,m′∈[0,N-1];n,n′∈[0,2M-1]。

结合图1的互素阵列结构,可发现:式(10)的前两种情况为“自差”协方差函数值(由同一子阵上的两阵元的快拍做“自相关”而得),后两种情况为“互差”协方差函数(子阵1的某阵元的快拍与子阵2的某阵元的快拍做“互相关”而得)。根据中国余数定理[17,24-25],任取l∈{-MN,…,0,…,MN},在以上自差和互差的4种情况中,至少有一种情况与φ(l)相对应。即

(11)

联立式(2)、式(8)～式(11)可知,在获得升序坐标集合P后,为把Nyquist协方差函数的空间时延l与观测阵元的协方差矩阵Rxx的行标u、列标v关联起来,不妨构造一张尺寸为L×L的“差集表”T,该差集表的内部元素为

T(u,v)=pu-pv,u,v∈[0,L-1]

(12)

由以上分析可看出,差集表T的生成仅涉及式排序和式减法操作,故过程简单。

2.2 基于差集表遍历搜索的Nyquist协方差函数估计

注意在差集表T中,对于某Nyquist空间时延l∈{-MN,…,0,…,MN},可能存在有多个表目T(u,v)值等于l的情况,故需把对应的行标rl、列标值cl全部搜索出来,即

(13)

(14)

借助式(14)计算出2MN+1个协方差函数估计值,即可构造如式(15)的Nyquist虚拟阵列的协方差矩阵

(15)

相比较而言,文献[22]则需要构造(2MN-N+1)×(2MN-N+1)的Toeplitz矩阵,然而由于在延迟范围l∈{-MN-M+1,…,0,…,MN+M-1}以外的波程差不连续,故在扩展的协方差矩阵中存在零元素(其数目随M、N值呈指数增长),故需以矩阵秩最小化为目标函数对这些令元素做填充,该最优化虽然可等效为核范数最小化问题[22],但在矩阵维度(或秩)很大时求解复杂度依然较大。

需强调,式(15)构造的毕竟是虚拟的协方差矩阵,与密集阵列构造的真实协方差矩阵有所不同。具体体现在:当各阵元采集的快拍数受限时,真实的协方差矩阵的对角线元素仍存在差别,而式(15)的对角线元素是相等的。这种虚拟构造方式是以略微牺牲估计器的抗噪鲁棒性为代价的,后面的仿真实验也将证明这一点。

2.3 DOA估计算法流程总结

根据以上描述,可将本文提出的互素稀疏阵列DOA估计算法总结为表1的流程。

表1 基于差集表遍历搜索的DOA估计流程Table 1 DOA estimation flow based on difference-set table traversal searching

3 仿真实验

3.1 本文算法的处理流程解析与性能对比

例1令目标源个数D=4,其信号功率均设置为1,其入射角θ1,θ2,…,θ4分别为-10°,20°,25°,60°,阵元数目L=6,每个阵元上采样的快拍数K=2 000,信噪比(signal-to-noise ratio,SNR)设置为SNR=0 dB(为高斯白噪声),试分别用互素稀疏阵列和传统Nyquist间距的均匀线性阵列做信号接收,并分别用本文提出的差集表遍历搜索法和传统直接MUSIC分解法进行DOA估计,对照其空间谱。

(1) 互素稀疏阵列下基于差集表遍历搜索的DOA估计结果

设定互素整数对M=2,N=3(恰好L=N+2M-1=6),根据式(1)构造互素稀疏阵列,根据式(3)、式(4)的系统模型而得到各阵元快拍数据,对这些数据处理过程如下。

根据表1步骤1,可得到升序排列的坐标集P={0,2,3,4,6,9}和6×6的差集表T(见表2)。

根据表1步骤2对表2列出的T(u,v)做遍历搜索,可得到如表3所示的坐标集{(rl,cl)}。

表2 基于差集表遍历搜索的DOA估计流程Table 2 DOA estimation process based on difference-set table traversal searching

表3 对差集表做遍历搜索后的坐标集Table 3 Coordinate sets by difference-set table traversal searching

(2) 进而根据表1步骤3、步骤4,可得到如图2所示的空间谱

图2 本文方法得到的空间谱图(D=4)Fig.2 Spatial spectrum by proposed method (D=4)

(3) 传统Nyquist间距下的均匀线性阵列接收结果

图3给出了传统Nyquist间距的ULA情况,直接用经典MUSIC分解法得到的空间谱。

图3 传统Nyquist-ULA的空间谱图(D=4)Fig.3 Spectrum of traditional Nyquist-ULA (D=4)

对比图2和图3,可得出如下结论。

(1) 虽然耗费同样数目的接收阵元(L=6)和快拍(K=2 000),用本文提出的基于差集表遍历搜索算法的互素阵列DOA估计精度高于基于传统Nyquist-ULA的估计精度,表现在前者的谱峰(见图2)可精确定位在期望的-10°,20°,25°,60°位置上(虚线所示),而后者的谱峰(见图3)与这些位置有明显偏离。

(2) 图2的DOA谱峰形状明显比图3的谱峰形状尖锐得多,表现在本文方法仍能分辨出两个密集入射角θ2=20°,θ3=25°,而图3的基于Nyquist-ULA的方法则只能将这两个入射角识别为1个。这是因为互素阵列的归一化阵元排列位置为P={0,2,3,4,6,9},而Nyquist-ULA的阵元排列位置为{0,1,2,3,4,5}。因而,在耗费同样阵元数目情况下,互素阵列由于阵元稀疏放置,具有比传统Nyquist阵列更大的孔径,结合本文提出的差集表遍历搜索方法,可获得更高的空间谱区分度。

(3) 进一步分析,两者性能差异的根本原因在于:如前述,依据本文提出的差集表遍历搜索算法可构造出(MN+1)×(MN+1)=7×7的协方差矩阵,从而容许的信号子空间最高维度为MN=6,而传统Nyquist-ULA只能构造出L×L=6×6的协方差矩阵,故容许的信号子空间最高维度为L-1=5。故对D=4个目标源情况,本文信号子空间容许度高于后者,因而获得了更高的接收性能。

进而可推出:M、N值越大,本文方法的性能改善更明显,下面给出例子说明。

3.2 提高阵元数目情况的多目标检测性能对比

例2仍分别用互素稀疏阵列和传统Nyquist均匀线性阵列做信号接收,令目标源个数D=15,其入射角θ1,θ2,…,θ15设置为[-50°,62°]范围内的15个均匀间隔值,阵元数目L=10(对于互素阵列则要求M=3,N=5),其他参数设置与例1相同。图4(a)、图4(b)分别给出了本文方法和传统直接MUSIC分解法的DOA识别谱图。

图4 不同阵列设置下两种方法的DOA检测谱图Fig.4 DOA spectrums of different array configurations based on prosed method and traditional method

从图4的DOA检测谱图可发现:图4(a)中,基于互素阵列的本文方法仍能区分[-50°,62°]范围内的15个谱峰,而对于Nyquist-ULA的检测方法,由于目标源数D=15超出了其允许的信号子空间最高维度(即L-1=9),导致图4(b)中的DOA检测基本失效。这验证了差集表遍历搜索法可提升信号子空间的容许维度的结论。

3.3 检测精度和鲁棒性对比

例3仍分别用互素稀疏阵列和传统Nyquist-ULA做信号接收,令目标源数D=1,其入射角θ1设置为45°,其他参数设置与例2相同。在-20～25 dB的SNR区间内(间隔为1 dB)对两种方法做蒙特卡罗测试(对于每种SNR情况,做2 000次测试),并统计均方根误差(root-of-mean-square-error,RMSE),其RMSE曲线如图5所示。

图5 RMSE随信噪比变化情况(L=10,M=3,N=5) Fig.5 Change of RMSE with SNR (L=10,M=3,N=5)

从图5测试结果可得出结论:

(1)整体上,本文方法的RMSE低于传统Nyquist-ULA的RMSE,因而具有更高的估计精度,这仍是差集表遍历搜索方法可提升信号子空间的容许维度、增强谱峰分辨能力的反映。

(2)在低SNR区域,两个估计器都会出现阈值效应(即SNR低于该阈值时估计器失效),从图5可看出,本文方法的SNR阈值比传统Nyquist-ULA的SNR阈值要高1 dB,这说明前者相比于后者,抗噪鲁棒性略低一些,这可看作是构造虚拟协防差矩阵所付出的代价。

统计RMSE随快拍数量K的变化,其中SNR设为10 dB,K取100～2 000间隔为100,结果如图6所示。

图6 波达角RMSE随快拍数量变化情况(L=10,M=3,N=5)Fig.6 Change of DOA RMSE with snapshot numbers (L=10,M=3,N=5)

可以得出结论:随着快拍数量的增加,两种估计器的准确度均会提高。在快拍数量一定的条件下基于差集表遍历搜索的估计方法误差更小。

3.4 计算时间对比

在Matlab R2016b(Windows10)平台上对两种估计器的执行时间进行比较,处理器为i5。两种算法分别执行300次,MUSIC扫描间隔设为0.05,阵元数量为10(N=5,M=3),2 000个快拍。提出的估计器用时30.0 s,Nyquist-ULA估计器用时28.5 s。仅仅多消耗了5.2%的时间就可以获得O(MN)的自由度,获得更高的估计精度。

4 结 论

本文以互素阵列为基础,提出基于差集表遍历搜索的DOA估计方法。该差集表仅需依据互素稀疏阵列的阵元坐标即可构造出来,以该差集表为指导,无需任何附加的优化迭代措施,即可实现观测阵元的协方差矩阵到Nyquist虚拟阵列的快速转换。理论分析和仿真实验均证明:借助差集表遍历搜索方法,可以保证互素阵列只耗费L=N+2M-1个阵元,即可获得MN个源目标的区分度,从而相比于Nyquist间距的传统阵列,大大提升了空间多目标分辨能力和DOA估计精度。

近年来随着S波段、C波段、X波段、K波段等高频段、短波长雷达应用逐渐展开,阵元间耦合的副效应逐渐趋于严重。由于本文提出的基于互素稀疏阵列的差集表遍历搜索法可放宽传统Nyquist阵列的密集布置阵元的限制,因而具有较广泛的工程应用前景。