低信噪比下非冗余阵列的无网格DOA估计

王 宁, 吕晓德, 李苗苗,3

(1. 中国科学院空天信息创新研究院, 北京 100094; 2. 微波成像技术国家级重点实验室, 北京 100190;3. 中国科学院大学电子电气与通信工程学院, 北京 100049)

0 引 言

波达方向(direction of arrival,DOA)估计是指从阵列中多个接收天线采集的数据中获取电磁波信号方向信息的过程。它是阵列信号处理中的一个重要问题,广泛应用于雷达、声纳、通信、电子对抗、语音信号处理等领域[1-3]。DOA是信号模型中的非线性参数,其均方根误差(root mean square error, RMSE)曲线会产生阈值效应。当信噪比或快拍数量低于某个阈值时,容易出现偏离真实值的全局误差,RMSE迅速增大并偏离Cramer-Rao 下界(Cramer-Rao lower bound, CRLB)。这种现象是许多非线性参数估计问题的共同特征[1]。在实际应用中,通常会选择阈值较低的算法,以适应较低的信噪比和较少的快拍。从估计问题的模型来看,DOA估计方法可以分为3类:参数化法、非参数化法和半参数化法[4]。

参数化方法往往将信号模型建立为参数的已知函数,在一定的最优准则下构造优化问题,通过优化方法直接得到参数的估计值。最经典的参数化方法是非线性最小二乘法(nonlinear least squares, NLS),可以从最大似然(maximum likelihood, ML)估计的角度推导出来[1]。当有足够的快拍时,ML估计渐近有效,渐近无偏并达到CRLB[5],因此NLS具有良好的统计特性。但由于NLS是非线性非凸问题,优化过程不能保证全局最优值,计算量大。多重信号分类(multiple signal classification, MUSIC)是一种基于协方差矩阵的特征值分解(eigenvalue decomposition, EVD)的参数化方法。为了利用导向矢量和噪声子空间之间的正交性来估计DOA,将特征向量分为噪声子空间和信号子空间两个子集。常用的方法有MUSIC、Root-MUSIC和最小模法[6]。当快拍数足够且源不相关时,这些算法可以接近 ML 的统计特征;但当快拍数较少且信噪比较低时,这些算法无法很好地区分噪声子空间和信号子空间,降低了DOA估计的准确性,且Root-MUSIC不能用于非冗余阵列。 基于旋转不变技术的信号参数估计(estimating signal parameter via rotational invariance techiques, ESPRIT)也是一种基于EVD的DOA估计方法,需要两个相同的子阵,利用子阵之间相位的旋转不变性来估计DOA[1,6]。与 MUSIC 一样,估计协方差矩阵也需要足够的快拍。参数化法可以直接得到DOA,而非参数化法需要形成谱,通过找到谱的最大值来间接确定DOA。较常用的非参数方法包括传统的波束形成和Capon方法。传统的波束形成器基于傅里叶分析,性能受到低分辨率和频谱泄漏的限制。Capon方法与 MUSIC 和 ESPRIT方法一样,是基于协方差矩阵的方法,需要足够的快拍来估计协方差矩阵。

半参数方法也可以称为稀疏方法，其充分利用了协方差矩阵的Toeplitz Hermitian 特性,可用于少快拍、单快拍、信噪比低的情况，还具有高分辨率特性。在过去的10年里,半参数方法受到了学者们的关注并得到了快速的发展。半参数方法可以分为3类[7]:网格、离网和无网格。 网格类方法对 DOA 值域进行网格划分,网格的数量远大于实际的源数量。假设DOA的真实值在某个网格中,从而可以以稀疏表示的形式构建信号模型，然后可以使用稀疏恢复算法来估计 DOA。常用的网格类方法包括平方根最小绝对值收敛和选择算子(square-root least absolute shrinkage and selection operator, square-root LASSO)[8]、迭代最小化稀疏学习(sparse learning via iterative minimization,SLIM)[9]、稀疏贝叶斯学习(sparse Bayesian learning,SBL)[10]和基于协方差的稀疏迭代估计(sparse iterative covariance-based estimation,SPICE)[4,11-13]。由于DOA的真实值可能不在网格上,因此网格类方法存在网格误差,限制了估计的准确性。离网类方法和网格类方法都需要对DOA 值域进行网格划分,但是在构建模型时,离网类方法通过引入网格偏移或使用动态网格来克服网格误差,例如稀疏总体最小二乘(sparse total least-squares,STLS)[14]、离网的稀疏贝叶斯推断(off-grid sparse Bayesian inference,OGSBI)[15]、网格进化[16]。然而,离网类方法往往会引入更多的参数,使算法更加复杂,而且大多数算法都涉及非凸优化,无法保证全局收敛。文献[17]采用“estimate and subtract”的策略提出了一种简洁的离网类方法,可用于非均匀线性阵列(non-uniform linear array,NLA),但在中等信噪比时会出现阈值效应。无网格类方法不需要网格,可以直接在 DOA 的连续域上进行估计,并且大多数算法都是凸优化问题,可以保证全局收敛。无网格类方法通过凸优化获得协方差矩阵的估计,然后使用Root-MUSIC或范德蒙分解计算DOA,但其中许多仅限于均匀线阵。无网格方法可以分为两类:基于原子范数理论的方法[18-21]和基于协方差匹配准则的方法[22-24]。文献[20-21]将无网格类方法扩展到任意一维阵列,不过文献[21]涉及非凸优化问题,不能保证全局收敛。

本文主要对原有的无网格SPICE进行改进,使其在较低信噪比和非冗余阵列下具有更好的性能。为了降低阈值效应的阈值,本文对Root-MUSIC进行了修改。新的 Root-MUSIC 使用了更多的特征向量,获得了比实际 DOA 数量更多的 DOA 估计值,并根据 ML 准则从中选择合适的值作为最终估计结果。同时,本文还在原始无网格SPICE的优化问题中加入了最大熵准则,可以使非冗余阵列情况下的RMSE曲线更接近CRLB,同时保持优化问题的凸性。将所提出方法的性能与 MSUIC、用于NLA的无网格SPICE[23]、交替投影无网格(alternating projections gridless, AP gridless)[21]和 CRLB 进行了比较。蒙特卡罗实验表明,本文提出的方法具有较低的阈值,可以在较低信噪比的环境中使用,并且还具有更好的高分辨率特性。同时,本文所提出的方法可以在相干源和信源数大于阵元数的情况下使用。

本文的后续内容组织如下:第1节介绍了经典的远场窄带模型、阵列信号处理中的冗余和非冗余阵列;第2节介绍了具有最大熵和基于ML-Root-MUSIC 的无网格 SPICE;第3节通过仿真验证了所提出的方法在低信噪比和非冗余阵列的情况下的有效性;第4节进行了总结。

1 信号模型

假设K个远场窄带平面波照射到M元均匀线阵,DOA记为θ=[θ1,θ2,…,θK]T,(θk∈(-90°,90°]),本文假设K是已知的。那么阵列接收到的数据可以写为

A(θ)s[n]+e[n],n=1,2,…,N

(1)

式中:n是快拍的序号;N是快拍数;y[n]∈CM是n时刻快拍数据;a(θk)∈CM是第k个信号的导向矢量;sk[n]∈C是第k个信号在n时刻的数据;e[n]∈CM是噪声;A(θ)=[a(θ1),…,a(θK)]被称为阵列流形。

可以用更紧凑的方式表示为

Y=A(θ)S+E

(2)

式中：Y=[y[1],y[2],…,y[N]]∈CM×N;S∈CK×N;E∈CM×N。本文对信号模型的统计特性做了如下假设:

(1) 源信号均值为零,并且在空域和时域中不相关,即

E[s[i]sH[j]]=diag(p)δi, j

(3)

(2) 噪声向量的均值为零并且在空域和时域中都是白色,即

E[e[i]eH[j]]=σ2Iδi, j

(4)

(3) 噪声向量和源信号不相关,即

E[s[i]eH[j]]=0

(5)

式中:p=[p1,p2,…,pK];pk∈R+表示源信号的功率;σ2表示噪声方差;δi, j是离散冲激函数。接收数据的协方差矩阵可以表示为

R=E[y[n]yH[n]]=

A(θ)diag(p)AH(θ)+σ2I

(6)

如果阵列为均匀线阵,则阵元的位置可表示为d=[1,2,…,M]T,距离单位为半波长。导向矢量可以表示为式(7)。fk称为空间频率,1是元素全为1的列向量。因为a(θk)有复正弦的形式,那么A(θ)就是一个范德蒙矩阵。可以进一步证明R是Toeplitz Hermitian矩阵且rank(R)=M。

fk∈(-0.5,0.5]

(7)

当阵列为非均匀线阵时,阵元的位置可表示为d=[d1,d2,…,dL]T(dl∈Z+)。不失一般性,令d1=1,dL=M。用ad(θk)、Ad(θ)、yd[n]、Rd分别表示非均匀阵列情况下的导向矢量、阵列流形矩阵、快拍和协方差矩阵。由于d不再是包含连续正整数的向量,不再具有复正弦的形式,则Ad(θ)不再是范德蒙德矩阵,而Rd只是Hermitian矩阵。我们可以定义一个选择矩阵Γd,只有第i行(i∈1,2,…,L)第di列是1,其余位置是0。然后可以得到:

yd[n]=ΓdAd(θ)s[n]+Γde[n]=Γdy[n]

(8)

Rd和R的关系如下所示:

(9)

为了显示Rd中包含R的多少个元素,可以定义d的合成D[22,25]。如果D={1,2,…,M}, 则Rd包含R中的所有元素,d对应的阵列就是冗余阵列。如果D不包含{1,2,…,M}中的所有元素,则Rd不包含R中的所有元素,并且d对应的阵列是非冗余阵列。例如,当d=[1,2,5,7]T时为冗余阵列:

D={m1-m2+1:m1,m2∈d,m1≥m2}

(10)

此时，Rd和R如下所示：

(11)

当d=[1,3,5,7]T时,该阵列为非冗余阵列,Rd如下所示：

(12)

2 ME-GLS-ML

本节介绍在原始无网格SPICE的基础上提出的方法,称之为联合最大熵和最大似然思想的无网格方法(maximum entropy-gridless SPICE-ML, ME-GLS-ML)。一方面,为了使无网格SPICE可以用于非冗余阵列,在优化问题中引入了最大熵的思想,加入了负熵项。同时,负熵项可以约束协方差矩阵是正定的。 另一方面,为了降低阈值效应的阈值,在Root-MUSIC中,通过ML准则选择根而不是选择最接近单位圆的根。

2.1 Gridless SPICE

Gridless SPICE 是一种基于协方差拟合的无网格DOA 估计方法,使用R的Toeplitz Hermitian属性将其参数化为r=[r0,r1,…,rM-1]T的函数R=T(r)。T(r)表示Toeplitz Hermitian矩阵,其第一行是rT,并且是rT的线性函数,这避免了使用DOA作为参数形成高度非线性函数[7]。以r作为参数的优化问题是一个凸优化问题,很容易求解。

(13)

将R=T(r)代入式(13)后可以得到

s.t.T(r)≥0

(14)

可以看到,式(14)是一个半正定规划问题,是一个凸优化问题[26],可以使用Matlab中的CVX工具箱来解决[27]。 一旦得到最优解,就可以构造R,然后可以使用Root-MUSIC算法计算DOA,从而达到无网格DOA估计的目的。

(15)

(16)

(17)

2.2 加入最大熵准则

从式(16)和式(17)可以看出,在非冗余阵列的情况下,优化问题与r中的某些变量无关。例如d=[1,3,5,7]T, 式(16)和式(17)只能与r0、r2、r4、r6相关,而不能与r1、r3、r5相关,因此只能得到r0、r2、r4、r6的值。

非冗余阵列引起的问题可以理解为协方差矩阵的r0、r2、r4、r6已知,确定r1、r3、r5的方法如下所示。

R=T([r0,?,r2,?,r4,?,r6])

(18)

式中:？表示该处值未知。

在谱估计中,周期图法无法估计具有高延迟的相关函数的值并将其设置为0,但突然转变为0会导致频谱中出现虚假的频谱分量。为了解决这个问题,Burg方法将相关函数的高延迟设置为对数据做出最少假设的值,即最大化随机过程的熵[28]。对于式(18)所示的问题,也可以利用最大熵的思想来获得协方差矩阵的未知值[25]。

假设y[n]是圆复高斯的,则其熵为ln[(πe)MdetR][29-30]。将熵代入式(16)和式(17)得到：

(19)

λln[detT(r)]

(20)

最大化熵相当于最小化负熵,所以熵项前面有一个负号。忽略与优化问题无关且大于零的常数项。由于负熵是一个凸函数,所以新生成的优化问题(式(19)和式(20))仍然是凸优化问题,可以使用CVX工具箱解决。负熵可以看作是约束T(r)>0的障碍函数。同时,其进一步在最大熵的统计意义下限制了T(r)。

2.3 ML-Root-MUSIC

(21)

(22)

PA(θ)=I-A(θ)(AH(θ)A(θ))-1AH(θ)

(23)

在本文中,使用ML从pk个候选根中选择K个根,如下所示：

(24)

称这种新方法为ML-Root-MUSIC。在随后的仿真中,发现ML-Root-MUSIC可以降低阈值,从而可以在较低的信噪比下应用无网格SPICE。

本文所提出的方法的处理过程如算法1所示。

算法1 ME-GLS-ML输入:观测数据Y,信号子空间维数ps,备选根数pk,源数目K,常数λ。输出:K个DOA。(1) 计算相关阵R^=(1/N)YYH;(2) 求解优化问题式(19)和式(20),得到R～;(3) 对R～进行EVD, 得到ps个最大特征值对应的特征向量,构成U～s, 剩余特征向量构成U～n;(4) 用U～n构造多项式并求解pk个最接近单位圆的根;(5) 用式(24)从pk个根中选择K个,记为z^i(i=1,2,…,k);(6) 输出θ^i=arcsin(argzi·π-1)(i=1,2,…,k)。

3 数值仿真

3.1 不同信噪比下的RMSE

本节通过蒙特卡罗仿真实验比较了所提方法与其他方法在非冗余阵列情况下的估计精度。 在仿真中,设置快拍数为200,阵元的位置为[1,3,6,8,11,14]T,该阵列是非冗余的。空间频率设置为[0.2,0.3,0.4]T,蒙特卡罗试验次数为1 000。源和噪声均使用圆复高斯白噪声,信噪比定义为10lg(p/σ2),p为源的方差,σ2为噪声方差,并且每个源的方差是相同的。CRLB的计算公式[1,6]为

(PAH(θ)R-1A(θ)P)T]}-1

(25)

(1) MUSIC:搜索谱峰时的网格数为4 096。

(2) NLA-GLS:文献[23]中用于NLA的无网格SPICE,在本文称为NLA-GLS。

(3) AP gridless:文献[21]中的方法。

(4) NLA-GLS+:使用了ML-Root-MUSIC的NLA-GLS,设置信号子空间的维数ps=3,设置备选根数pk=8。

(5) ME-GLS-ML:本文提出的方法,设置信号子空间的维数ps=3,设置备选根数pk=8。

5种方法非冗余阵列下RMSE随信噪比的变化如图1所示。可以看出,这5种方法在信噪比小于某个阈值后都会产生阈值效应。MUSIC和AP gridless的阈值为-6 dB,NLA-GLS的阈值为-2 dB,NLA-GLS+和ME-GLS-ML的阈值为-8 dB。NLA-GLS+和ME-GLS-ML的阈值最低。与NLA-GLS相比,由于NLA-GLS+增加了ML-Root-MUSIC,阈值从-2 dB降低到-8 dB。与NLA-GLS+相比,ME-GLS-ML增加了负熵项,RMSE曲线更接近CRLB,说明在非冗余阵列情况下增加负熵项有利于提高估计精度。在这5种方法中,ME-GLS-ML的阈值最低,最贴近CRLB 曲线。因此,本文提出的方法在非冗余阵列的情况下更具优势。

图1 非冗余阵列下RMSE随信噪比的变化Fig.1 Change of RMSE with signal to noise ratio in non-redundant array

图2是NLA-GLS和ME-GLS-ML在-8 dB信噪比下1 000次蒙特卡罗结果的直方图。可以看出ME-GLS-ML的估计结果比NLA-GLS更集中在真实值附近,NLA-GLS存在偏离真实值的异常值。

图2 1 000 次蒙特卡罗试验直方图Fig.2 Histogram of 1 000 Monte Carlo

3.2 分辨率

在本节增加如下两个方法的仿真。

(1) NLA-GLS++:设置NLA-GLS+信号子空间的维数ps=8,设置备选根数pk=8。

(2) ME-GLS-ML+:设置ME-GLS-ML信号子空间的维数ps=8,设置备选根数pk=8。

在仿真中,设置快拍数为200,阵元的位置为[1,3,6,8,11,14]T,该阵列是非冗余的。基于傅里叶分析的DOA估计方法的空间频率分辨率约为1/M[1,6]。一般能够使分辨率小于1/M的方法称为超分辨率方法。在仿真中生成两个源信号,其中一个具有固定空间频率0.2,另一个具有空间频率0.2+Δf/M,Δf为频率间隔(Δf∈{0.1,0.2,…,1})。两个源信号功率一样,信噪比固定为-4 dB。每个Δf分别进行200次蒙特卡罗模试验。当满足下式时,认为是成功分辨的[1]：

(26)

仿真结果如图3所示,所有方法均可实现超分辨。NLA-GLS+和ME-GLS-ML的超分辨性能并不比NLA-GLS好,但是当信号子空间的维数ps增加到8,大于真实DOA数目时,超分辨性能得到了改善,如NLA-GLS++和ME-GLS-ML+的曲线所示。这一结果表明,适当增加信号子空间维数有利于提高超分辨率性能。NLA-GLS++在所有比较方法中分辨性能最好,说明使用ML-Root-MUSIC可以提高分辨性能。

图3 分辨成功率Fig.3 Success rate of resolution

3.3 相干源

本小节通过仿真说明ME-GLS-ML在源信号相干的情况下依旧可以进行DOA估计,信噪比固定为0 dB,各源信号使用相同的圆复高斯白噪声来模拟相干源的场景,其余仿真条件与第3.1节一致。仿真结果如图4所示,经典的MUSIC方法在相干源的情况下会出现很多伪峰,无法得到正确得DOA估计,而本文所提出的ME-GLS-ML可以正确得到DOA的估计。

图4 相干源情况下的仿真Fig.4 Simulation with coherent sources

3.4 源数目大于阵元数

本节考虑源数目大于阵元数的情况,源数目设置为8,相应的空间频率分别设置为-0.4、-0.3、-0.2、-0.1、0.1、0.2、0.3、0.4,信噪比为0 dB,阵元位置的设置与第3.1节相同, 此时阵元数为6。仿真结果如图5所示,ME-GLS-ML可以成功得到所有源的DOA估计。由于ME-GLS-ML通过优化的方式得到了虚拟的14阵元均匀线阵所对应的协方差矩阵,因此可以对多余真实阵元数的源进行DOA估计,理论上最多可以对13个源进行DOA估计。

图5 源数目大于阵元数情况下的仿真Fig.5 Simulation with the number of sources greater than the number of array elements

4 结 论

本文主要关注DOA估计中的阈值效应,提出了一种适用于低信噪比和非冗余阵列的无网格DOA估计方法。该方法在原有无网格SPICE的基础上,引入最大熵的思想,构造了一个新的凸优化问题。该凸优化问题可以用来拟合非冗余阵列的协方差矩阵,然后通过ML-Root-MUSIC完成DOA估计。根据蒙特卡罗仿真实验,该方法具有更好的性能。该方法在相同数量的快拍下具有较低的阈值,可应用于信噪比较低的环境。同时,在相同的信噪比和相同快拍数量下,本文提出的方法具有更好的高分辨率特性。在 ML-Root-MUSIC中,使用了比实际信源数更多的特征向量,并且使用ML准则进行DOA选择,这个想法也可以扩展到其他基于EVD的方法。