基于压缩感知理论的二维DOA估计

窦慧晶， 梁 霄， 张文倩

(北京工业大学信息学部， 北京 100124)

二维波达方向(direction of arrival，DOA)估计在阵列信号处理领域有着重要的研究意义，与一维DOA估计相比，该估计算法能够更精确描述目标的空间特性，因此DOA估计在二维信号领域更具实际应用价值[1-2]. 二维多重信号分类(two-dimensional multiple signal classification, 2-D MUSIC)算法是目前已有的二维阵列信号DOA估计算法中最为经典的估计算法之一，该算法核心思想是将传统的一维MUSIC估计算法在二维空间进行直接推广，由于该算法需要二维谱峰搜索因而导致计算量巨大，且需要各信源的中心频率已知，因此很难满足实际应用[3]. 为了解决上述缺陷，有学者提出一种无须谱峰搜索的二维旋转不变子空间(two-dimensional estimating signal parameter via rotational invariance techniques, 2-D ESPRIT)算法以及二维传播算子(two-dimensional propagation method, 2-D PM)算法[4]. 这些算法的相继问世使阵列信号的处理性能得到一定的提高，但因其在小快拍数及低信噪比情况下估计性能严重下降而无法推广到实际应用中.

在众多阵列结构中，由于L型阵列具有结构简单、实施容易、估计性能佳等优点而被广泛用于工程领域. 为解决二维信号角度匹配精度不高且计算复杂的问题，文献[5]提出一种基于L型阵列的无须手动配对的二维DOA估计算法，通过引入新的合成角度计算出新的导向矢量，进而获得原信号的俯仰角和方位角. 尽管该方法能够自动完成角度配对，但需要多次谱峰搜索及特征值分解导致计算复杂度过高. 文献[6]提出一种新的二维DOA估计方法，该算法首先将方位角和俯仰角分别估计出来，再通过阵列输出的互相关和信号功率对2个角度进行匹配，由于需要大量的采样信号使得该方法不可有效避免大量的数值计算. 为降低运算量有学者提出利用阵列数据的协方差矩阵进行二维角度估计的算法[7-8]. 文献[9]提出一种利用多相干信号对方位角和俯仰角进行配对的方法，通过利用协方差矩阵最小化构造的代价函数从而实现角度配对，该算法存在的最大弊端是在构造协方差矩阵的过程中可能会引入外界噪声，从而影响其估计性能.

压缩感知(comprehensive sensing，CS)理论的出现为现代信号处理带来一种更高效、更精确的方法，文献[10]提出基于该理论的l1-SVD算法，该算法通过对接收信号进行奇异值分解(singular value decomposition, SVD)来降低算法复杂度和对噪声的敏感性，然后利用二阶锥规划的方法求解相应的优化问题，该算法在小快拍数和低信噪比时有很好的性能，并且可以直接用于相干信号[11]. 该方法摆脱了传统奈奎斯特采样定理带来超大计算量的束缚. 基于此，众多学者将压缩感知理论引入到DOA估计中来，从而达到降低计算量的目的. 文献[12]提出一种基于协方差矩阵联合稀疏重构的降维波达方向估计算法，该算法充分利用阵列孔径，无须预先估计目标数目，参数估计性能在低信噪比及小快拍数据长度下优势明显，但在其他方面尚有改进余地. 本文在基于压缩感知理论的背景下提出一种二维L型阵列信号的DOA估计算法. 该方法在高信噪比、多快拍条件下相较于传统算法具有更高的估计精度和分辨力，且具有较低的运算量.

1 信号模型

本文试验采用L型均匀阵列，该模型中2个子阵互相垂直，成90°夹角，且分别位于x轴和y轴上，如图1所示. 该模型的2个子阵中均含有N个阵元，且子阵X和子阵Y共用L型阵列的原点. 本文规定该模型存在k个非相关的远场窄带信号(k=1,2,…,K)，且每个阵元间的距离为半个波长，入射信号与x轴、y轴的夹角分别为α、β. 定义入射方向角θ、φ分别为信号的方位角、俯仰角. 在某时刻t接收信号的数据模型为

x(t)=AX(θ,φ)s(t)+nx(t)

(1)

y(t)=AY(θ,φ)s(t)+ny(t)

(2)

式中：t=1,2,…,T，T表示快拍数；s(t)=[s1(t),s2(t),…,sK(t)]T表示K个信源t时刻信号矢量；导向矢量矩阵为

AX(θ,φ)=[ax(θ1,φ1),ax(θ2,φ2),…,ax(θk,φk)]

(3)

AY(θ,φ)=[ay(θ1,φ1),ay(θ2,φ2),…,ay(θk,φk)]

(4)

其中

ax(θk,φk)=[1,η1,η2,…，ηM-1]T

(5)

ay(θk,φk)=[1,μ1,μ2,…，μM-1]T

(6)

式(5)(6)分别为子阵X、子阵Y的导向矢量，η=e-j2πdcos (θk)cos(φk)/λ，μ=e-j2πdcos (θk)sin(φk)/λ.

由图1空间几何位置关系可知

图1 均匀L型阵列模型Fig.1 L-type array signal model

(7)

将式(7)分别带入式(5)(6)得到新的导向矢量

ax(αk)=[1,e-j2πdcos (αk)/λ,e-j2πd(M-1)cos(αk)/λ]

(8)

ay(βk)=[1,e-j2πdcos(βk)/λ,e-j2πd(M-1)cos(βk)/λ]

(9)

2 算法原理

2.1 压缩感知基本原理

压缩感知理论指出：信号如果满足稀疏的或者可压缩，就可以对信号进行远低于奈奎斯特的采样频率采样，并能实现对原始信号的精确重构[13-14]. 该理论的核心由以下3个部分组成：

1) 对原始信号进行稀疏表示；

2) 针对不同情况构造合适的测量矩阵；

3) 对上述信号进行稀疏重构[15].

假设信号x为离散信号，长度为N，可以通过一组正交基来表示信号，即

(10)

式中：ψ=[ψ1,ψ2,…,ψN]表示标准正交基，α=[α1,α2,…,αN]T表示系数矢量. 对于权系数向量α，可以定义为

α=ΨTx

(11)

假如当0

(12)

式中：‖α‖p表示向量α的p范数，则表明权系数向量α在某变换域上可进行稀疏表示，同时表明信号具有可压缩性[16].

在信号x可进行稀疏表示或可压缩的前提下，利用不同的观测矩阵对采样信号进行投影变换，即可将原始高维信号降低维度，从而得到所需的稀疏信号[17-18].

观测矢量为

y=Φx=Φψα=Θα

(13)

式中：Θ=Φψ，Φ为Q×N维测量矩阵，N≫Q且Q≥O(Klg (N/K))，通过式(13)可以实现信号的稀疏采样.

对于压缩感知矩阵Θ的构造，文献[19]通过数学公式推导的方式证明Θ在满足限制等距原理(restricte isometry principle, RIP)准则时，可选取最小l0范数法、正交匹配追踪等算法求解式(13)，进而高精度地完成基于压缩感知的二维DOA估计过程.

2.2 基于过完备冗余字典的稀疏表示

在基于压缩感知理论下的DOA估计中，为了提高算法估计精度，文献[20]提出一种对阵列的流形矩阵进行扩展的方式构造二维信号的过完备冗余字典D的方法，该方法可以使字典D中尽可能包含所有存在的空间角度.

D=[a(θ1,φ1),a(θ2,φ2),…,a(θk,φk)]

(14)

式(14)是由一维DOA估计的经典算法构造字典的经典方式直接将其推广到二维DOA估计领域得到过完备冗余字典，其中a(θk,φk)称为原子. 在式(14)中，由于信号维度的扩展，致使字典原子数成倍扩增，为后续处理带来大量不必要的计算. 为避免一维经典算法直接推广带来计算困扰，本文采用对第1章所提出空间合成角α和β分别构造过完备冗余字典的估计算法，进而得到2个空间合成角的DOA估计，最后通过参数匹配得方式找出其中的关系，从而高精度地还原出二维信号的方位角和俯仰角.

由式(8)(9)可知，定义的空间合成角α和β已包含二维信号的方位角和俯仰角等信息，同时已经满足空域稀疏性的要求. 接下来分别对α和β构建超完备冗余字典

Dx=[ax(α1),ax(α2),…,ax(αN)]

(15)

Dy=[ay(β1),ay(β2),…,ay(βN)]

(16)

在构建过完备冗余字典后，由于多数文献采用的空间等角度网格划分算法存在恢复精度不高的缺陷，本文进行一定改进，利用一种空间等余弦网络划分的方法构建过完备冗余字典，及通过取cosα=-1+(n-1)/N,n=1,2,…,2N+1，将式(1)(2)转换为稀疏表示形式

(17)

式中hx和hy分别代表空间合成角α和β的系数矢量.

2.3 测量矩阵的构造

如何设计测量矩阵Φ是信号稀疏表示后的最重要一步. 压缩感知理论指出：在信号可以进行稀疏表示的情况下，通常利用一个远小于信号维度的观测矩阵对该信号进行压缩采样，从而利用远低于奈奎斯特采样的方式实现信号重构，但测量矩阵的构造需满足RIP准则.

高斯随机观测矩阵是满足RIP准则最常用测量矩阵构造方法之一，该矩阵元素在(0,1)区间随机独立正态分布性质使它几乎与任何稀疏信号都不相关. 但该矩阵的非结构化性质为计算带来较高的复杂度. 将高斯矩阵正交化可有效增强它的结构化的特性，为后稀疏分解过程降低计算复杂度，因此本文选择正交化的高斯矩阵作为观测矩阵. 具体方案如下：

1) 定义高斯测量矩阵Φ的RIP参数δk为满足公式

(18)

的最小值δ，式中x为K阶稀疏信号. 若δk<1，则称测量矩阵Φ满足K阶RIP.

2) 若PTP=I，式中I为单位测量矩阵，则称P为正交矩阵. 设高斯随机矩阵Φ为M×N维，P0和P2为N×N维的正交矩阵，P1为M×M的正交矩阵，可以通过公式

Φ1=ΦP0

(19)

Φ2=P2ΦP1

(20)

得到2个新的观测矩阵.

3) 因为Φ1=ΦP0，所以有

(21)

又有

(22)

其中Φ满足δk的RIP，所以

(23)

由式(19)(20)得

(24)

即Φ1也满足RIP准则，同理可得Φ2也满足RIP准则，因此Φ1、Φ2均可用于构造单位测量矩阵.

2.4 稀疏表示的求解

由上述观测矩阵和过完备冗余字典可以得到观测矢量

xcs=ΦDXhx+Φnx=ΘXhx+vx

(25)

ycs=ΦDYhy+Φny=ΘYhy+vy

(26)

观测模型为

Xcs=ΘXHX+VX

(27)

Ycs=ΘYHY+VY

(28)

式中：ΘX=ΦDX,ΘY=ΦDY,vX=Φnx,vY=Φny，VX=ΘXNX,VY=ΘYNY.

由式(20)可知，压缩采样可将原始信号数据从N×1维降低到M×1维，其中M≪N，再利用每次快拍获得的M×1维观测数据求解得到稀疏系数，即得到相应的空间频率，而稀疏系数的求解可以转化成式(27)(28)的优化问题. 大多数研究方向在于如何高精度低成本重构原始信号，考虑实际应用中信号的压缩采样过程是在模拟域进行，考虑到模拟器件的非理想性，因此建模时将压缩感知矩阵中含有噪声E的条件加入模型，并假设压缩感知矩阵中存在的噪声是和观测矩阵中存在的噪声相互独立且均为高斯白噪声，因此该问题的含噪模型可表示为

Xcs=(ΘX+EX)HX+VX

(29)

Ycs=(ΘY+EY)HY+VY

(30)

式中：ΘX、ΘY为分解后的字典；EX、EY分别为与压缩感知矩阵ΘX、ΘY同维度的噪声矩阵.

针对该模型通过改进RM- FOCUSS算法进行求解. 利用循环迭代的方法求解最优矩阵H(其中H为HX、HY总称)，RM- FOCUSS针对不含噪声的E时优化问题可描述为

(31)

其代价函数可表示为

(32)

式中

表示稀疏度的度量，γ=2λ/|p|≥0，p∈[0,1].

该算法的求解过程可以简单描述为

Wl+1=diag(cl[i])p/2

(33)

式中

(34)

式中：Θl+1=ΘWl+1(Θ为ΘX、ΘY的统称)；λ≥0.

Hl+1=Wl+1Ql+1

(35)

通过J(p)(H)的结果可知，矩阵HX、HY的列在lp(p≤1)范数意义下产生了稀疏效果，而行元素是由2范数结合的，所以并没有产生稀疏效果. 式(35)中l表示迭代次数，式(31)(32)中的参数γ表示平衡逼近误差与稀疏性的参数，对于不同的信噪比及解的稀疏度采用修正l曲线法可预先选取最优的值. 通过查阅相关文献得知0.8≤p≤1.0时，针对信号稀疏性与计算复杂度综合考虑可知：该算法与原算法相比，具有计算量相对较小且估计精度更高的优势.

当字典Θ中含有噪声E的情况下，本文将RM- FOCUSS算法进行改进，提出一种通过循环迭代从而求解最优矩阵H的方法，其代价函数可表示为

(36)

式中：Y代表Xcs、Ycs的总称；V代表VX、VY的总称；E代表EX、EY的总称. 在每次迭代中，H和E的所有元素依次循环更新. 换句话说，已知E，需要通过式(31)(32)求解代价函数Hl，求解出Hl后，可以在最小F-范数的限制下估计El+1. 也就是说：每次迭代过程中，需要在H和E之间轮流计算最优值.

该算法具体过程为：首先初始化噪声E0为与字典Θ维数相同的零矩阵，在l次迭代时，已知E的估计值为El，可推导出估计值Hl的表达式为

(37)

求得Hl后，通过求解优化问题可得估计值

(38)

对式(38)求一阶导数后置零，其最优解可表示为

(39)

多次重复以上步骤，直到收敛条件可以满足‖Hl-Hl‖F/‖Hl‖F<δ(本文δ=0.001)时，停止迭代过程，得到满足要求的稀疏解H=Hl+1.

通过上述方法即可求出满足条件的稀疏解HX、HY.

2.5 方位角和俯仰角的配对

相较于传统DOA估计算法需要分别求解2个子阵的观测矢量，以及接收数据的协方差矩阵进行特征值分解导致特征矢量的排列顺序产生变换的弊端，本文采用的方法可有效避免特征值分解所带来的冗余计算，最后利用频率或相关性对空间角进行匹配，再根据式

(40)

给出的空间角定义反推出所求的方位角和俯仰角.

3 仿真结果分析

3.1 实验一：3种算法有效性的比较

假设存在3个非相关窄带信号的入射角度分别是(60°，80°)、(100°，100°)、(140°，120°)，设定仿真实验的快拍数等于200，阵列中阵元总个数为16个，信噪比是-5 dB，且保持恒定不变. 在此条件下，分别对经典2-D MUSIC算法，参考文献[11]算法及本文算法进行仿真试验. 得到仿真结果如图2～4所示.

图2 经典2-D MUSIC算法仿真结果Fig.2 Traditional 2-D MUSIC algorithm simulation results

图3 文献[11]算法仿真结果Fig.3 Literature [11] algorithm simulation results

图4 本文算法仿真结果Fig.4 Proposed algorithm simulation results

通过图2～4的散点图可以直观看出，在低信噪比及小快拍数下的条件下，传统2-D MUSIC算法，文献[11]对信号的二维角度估计虽能看出俯仰角与方位角的大致位置，但估计效果与本文算法的估计效果相比而言还是有较明显的差距，本文的估计效果更佳.

3.2 实验二：通过均方根误差衡量估计性能

在DOA估计实验中，衡量估计性能的方法有很多种，最重要的指标为均方根误差的大小. 该定义为

(41)

1) 单一信源信号均方根误差随信噪比变化的对比实验. 假设该单信源位置是(60°,80°). 设定仿真实验快拍数为200且保持恒定，信噪比从-10 dB到30 dB以5 dB间隔进行变化. 比较此时3种算法在信噪比变化时的均方根误差曲线. 为了保证实验的准确性，每个信噪比下做200次蒙特卡罗实验，结果如图5所示.

图5 单信源不同算法随信噪比变化的均方根误差曲线Fig.5 RMS error curves of different algorithms for single

2) 相干信源均方根误差随信噪比变化的对比进实验. 假设空间存在位置为(60°,120°)和(100°,100°)的2个相干信源. 设定仿真试验快拍数为200且保持恒定，信噪比从-10 dB到30 dB以5 dB间隔进行变化. 比较此时3种算法在信噪比变化时的均方根误差曲线. 为了保证实验的准确性，每个信噪比下做200次蒙特卡罗实验，结果如图6所示.

图6 相干信源不同算法随信噪比变换的均方根误差曲线Fig.6 RMS error curves of different algorithms for coherent sources varying with SNR

由图5、6可以看出，随着信噪比的增大，3种算法的均方根误差均有不同程度降低. 对于相干信源信号角度的估计：在相同信噪比的条件下其余2种算法均方误差明显高于本文所提出的算法，并且在单信源的情况下，上述2种算法的估计精度也同本算法具有一定差距. 这是由于本文算法利用压缩感知理论得到更多信号信息，提高了阵列信号的利用率，进而提升了估计精度.

3) 均方根误差随快拍数变化的试验

比较3种算法在信噪比为10 dB且保持恒定，快拍数从20到200以10为间隔进行变化的均方根误差曲线. 为了保证实验的准确性，每个快拍数下做200次蒙特卡罗实验，实验结果如图7所示.

图7 不同算法随快拍数变化的均方根误差比较Fig.7 Comparisons of root mean square errors of different algorithms varying with the number of snapshots

通过图7可以看出，3种算法随着快拍数的增加均方根误差都有下降的趋势，当快拍数相同的时候，本文算法的均方根误差相较前2种算法更低，表明本文算法性能更优.

3.3 实验三：3种算法角度估计成功率比较

首先约定实验成功的定义为：估计角度与真实角度相差在±3°以内视即可视本次估计成功. 假设3个非相关窄带信号的入射角度分别是(60°,120°)、(100°,100°)、(140°,80°).

1) 比较3种不同算法随信噪比变化时的成功率变化，实验中将快拍数设置为50且保持恒定，信噪比从-10 dB到30 dB以2 dB间隔进行变化. 为了保证实验的准确性，每个快拍数下做200次蒙特卡罗实验，实验结果如图8所示.

图8 不同算法随信噪比变化的成功率比较Fig.8 Comparison of success rates of different algorithms with signal-to-noise ratio

2) 比较3种算法随快拍数变化时的成功率变化，实验中将信噪比设置为10 dB且保持恒定不变，快拍数从20到200以10间隔进行变化. 为了保证实验的准确性，每个快拍数下做200次蒙特卡罗实验，实验结果如图9所示.

图9 不同算法随快拍数变化的成功率比较Fig.9 Comparison of success rates of different algorithms with the change of snapshot number

从图8、9三种算法随信噪比和快拍数变化成功率的对比中可以较为直观地看出：随着信噪比或快拍数的增大，3种算法的成功率均有所提升，但在相同快拍数或信噪比条件下，本文算法与其他2种算法相比性能更优.

4 结论

本文提出的基于压缩感知理论下二维DOA估计方法对压缩感知理论与DOA估计进行结合，首先将俯仰角与方位角进行分解，分别构建过完备冗余字典；接着采用正交化的方法构造高斯随机矩阵；最后利用改进的RM- FOCUSS算法和反三角余弦求解出所求的俯仰角与方位角. 并通过MATLAB仿真试验验证了算法的可行性. 在多快拍、低信噪比情况下分别研究单信源和相干信源估计精度. 本文所提算法和传统2-D MUSIC算法以及文献[11]算法相比具有更高的估计性能. 本文算法具有如下优点：

1) 利用空间合成角的方法将二维信号有效转换为2个一维信号，降低了估计的难度.

2) 利用压缩感知的方法，通过少量信号即可高精度地还原出信号的俯仰角及方位角，可有效降低计算量.

3) 通过采用正交化的方法构造高斯矩阵可以降低矩阵恢复误差，提升恢复精度.

4) 考虑到构造字典可能存在噪声的情况采用改进的RM- FOCUSS算法得到更加高精度的角度估计结果.