基于加权l 1最小化的低复杂度波达方向估计算法

基于加权l1最小化的低复杂度波达方向估计算法

段素馨1，2张颢1孙秀志2郑春弟3

(1.清华大学电子工程系,北京 100084;2.中国电子设备系统工程公司无线电管理部,北京 100840;3.海军陆战学院,广东 广州 510430)

摘要基于阵列协方差矩阵的稀疏表征和阵列响应矩阵的Khatri-Rao积,提出了一种低运算复杂度的波达方向估计算法. 所提算法在减少未知数个数的同时,通过线性变换降低约束方程的维数,可有效减少优化问题的计算复杂度. 为充分利用阵列协方差矩阵中蕴涵的信息,使用Capon谱的倒数作为权值构建出了加权l1最小化问题,这使得所提算法在降低运算量的同时能够获得较好的估计性能. 仿真实验验证了所提算法的有效性.

关键词波达方向估计;加权l1最小化;稀疏恢复;等距线阵

中图分类号TN911.7

文献标志码A

文章编号1005-0388(2015)04-0640-07

AbstractBased on the sparse representation of the array covariance matrix and the Khatri-Rao product of the array response matrix, a low computational complexity sparse recovery method for direction-of-arrival (DOA) estimation is presented. The proposed algorithm not only lessens the number of unknown variable, but also can cut down the dimension of the constraints, which considerably reduce the computational complexity of the second order cone programming. Moreover, a weighted l1 minimization is designed by using the reciprocal of the Capon spectrum as a weighting vector. As a result, the proposed algorithm can achieve better performance while the computational complexity is reduced. Simulations demonstrate the performance of the proposed method.

收稿日期：2014-03-04

作者简介

A low complexity algorithm based on the weightedl1

inimization for DOA estimation

DUAN Suxin1,2ZHANG Hao1SUN Xiuzhi2ZHENG Chundi3

(1.DepartmentofElectronicEngineering,TsinghuaUniversity,Beijing100084,China;

2.ChinaElectronicSystemEngineeringCompany,Beijing100840,China;

3.NavyMarinesCollege,Guangzhou,Guangdong510430,China)

Key words direction-of-arrival (DOA) estimation; weightedl1minimization; sparse recovery; uniformly-spaced linear array (ULA)

资助项目: 国家自然科学基金(61401496)

联系人： 段素馨 E-mail：007-yx@163.com

引言

作为阵列信号处理领域中的基本问题,波达方向(Direction of Arrival,DOA)估计在电磁学、声学、水声学、地震学和生物医学等领域有着广泛的应用[1-3]. 对于DOA估计算法而言,高分辨性能始终是研究的重点. 近年来,利用观测角度空间中潜在目标具有稀疏性这一信息,众多基于稀疏表征的DOA估计算法被提出[4-10],其中采用基追踪(Basis Pursuit,BP)准则的DOA估计算法由于求解方便且精度高成为最重要的研究分支[4-9]. 相比于传统的多重信号分类(Multiple Signal Classification, MUSIC)、子空间旋转不变技术(Estimating Signal Parameters via Rotational Invariance Techniques, ESPRIT)等子空间类DOA估计算法,基于稀疏表征的DOA估计算法最大的优势之一就是可以获得更高的分辨性能[4].

基追踪可用二阶锥规划(Second Order Cone Programming,SOCP)进行求解. 采用内点法的SOCP所需的运算量为O((N×T)3)[4],而子空间类算法的运算复杂度为O(M3),其中N、T和M分别为稀疏信号长度、快拍数和阵元个数,M≪N,因此采用BP准则的DOA估计算法的运算复杂度一般要远高于传统的DOA估计算法. 为了降低基于稀疏表征的DOA估计算法的运算复杂度,经典的l1-SVD算法[4]采用主分量分析(Principal Component Analysis,PCA)降低稀疏恢复的运算量,该算法使用截尾奇异值分解(Truncated Singular Value Decomposition, TSVD)将未知数的个数从N×T减少到N×K,其中K为信号个数,通常K≤T. 文献[7]在等距线阵(Uniformly-spaced Linear Array, ULA)条件下,将接收数据协方差矩阵进行稀疏表征,提出了一种低复杂度的稀疏恢复算法(本文称之为He’s method),该算法使用Khatri-Rao积在虚拟扩展阵列孔径的同时,将多观测向量(Multiple Measurements Vectors,MMV)稀疏恢复问题转换为单观测向量(Single Measurements Vectors,SMV),从而使得未知数的个数降低为N. 总的来说这两种算法均是从减少未知数个数这个角度来降低稀疏恢复的运算复杂度.

事实上,SOCP的运算复杂度还与二阶锥的维数有关[11],与上述两种方法不同,本文试图通过减少二阶锥的维数(二阶锥C={x=(x0,x1)∈R×Rk-1:x0-‖x1‖≥0}的维数,即k)来实现降低运算复杂度的目的. 与He’smethod方法相似,我们对接收数据协方差矩阵进行稀疏表征. 但是在使用Khatri-Rao积实现孔径虚拟扩展之后,采用线性变换将二阶锥的维数从M2降为2M-1. 另外,为了充分利用协方差矩阵中蕴涵的信息,本文使用Capon谱的倒数对目标函数进行加权处理[6,9-10],通过加权处理之后,恢复算法更倾向于在那些更有可能是真实DOA的位置安排非零元素,从而可以有效提升算法的性能.

1理论分析

1.1信号模型

y(t)=As(t)+n(t),t=1,2,…,T.

(1)

式中: y(t)=[y1(t),…,yM(t)]T∈CM×1为M个阵元上的接收数据向量; n(t)=[n1(t),…,nM(t)]T∈CM×1为M个阵元上的噪声向量; A=[a(θ1),…,a(θK)]∈CM×K为阵列响应矩阵,a(θk)=[1,ej2πdsin(θk)/λ ,

…,ej(M-1)2πdsin(θk)/λ]T∈CM×1为阵列的导向矢量. 假设噪声n(t)为均值为0、方差为σ2的高斯白噪声. 不失一般性,假设n(t)与s(t)不相关. 假设s(t)为不相关信号,则式(1)的协方差矩阵可以表示为

R=E{y(t)yH(t)}=ARsAH+σ2IM.

(2)

z≜vec(R) =(A*⊗A)vec(Rs)

=(A**A)p+σ2vec(IM).

(3)

式中: A**A=[a*(θ1)⊗a(θ1),…,a*(θK)⊗a(θK)]∈CM2×K ;符号“⊗”、“*”分别代表矩阵的Kronecher积和Khatri-Rao积; p=[p1,…,pK]T.

根据Khatri-Rao积的特点和ULA的假设,容易得到[12]

A**A=GB.

(4)

式中: G(i×M+1∶(i+1)×M,∶)=[0M×(M-1-i),IM,0M×i],i=0,1,…,M-1; B=[b(θ1),…,b(θK)]∈C(2M-1)×K为虚拟的阵列响应矩阵,相应的虚拟导向矢量b(θk)=[e-j(M-1)2πdsin(θk)/λ ，…，e-j2πdsin(θk)/λ ,1,ej2πdsin(θk)/λ ,…,ej(M-1)2πdsin(θk)/λ]T∈C(2M-1)×1. 相比初始的导向矢量a(θk),阵列孔径虚拟地从M扩展为2M-1,这是使用Khatri-Rao积带来的最大好处,阵列孔径的扩展将使得DOA估计算法的性能获得极大的改善. 将式(4)代入式(3)可得

z=GBp+σ2vec(IM).

(5)

(6)

令

(7)

根据文献[13]中的结论,有

(8)

文献[13]中指出,协方差矩阵的估计误差服从如下的高斯分布

(9)

式中,符号AsN表示渐进高斯分布．

1.2基于Khatri-Rao积的加权稀疏恢复算法

本节首先简要回顾文献[7]中提出的低复杂度稀疏恢复算法,在分析该算法运算量的基础上,寻求进一步降低运算量的途径,并将其应用于稀疏恢复,提出基于Khatri-Rao积的加权稀疏恢复算法.

1.2.1基于Khatri-Rao积的稀疏恢复算法[7]

在稀疏表征框架下,式(5)可以表示为

z=GΦx+σ2vec(IM) .

(10)

式中: Φ=[b(φ1),…,b(φN)]∈CM×N是超完备基矩阵,即有M≪N,集合{φ1,…,φN}表示观测角度空间[0,π]被均匀划分为N个格点所获得的角度集合; x∈RN×1表示稀疏信号,当且仅当φn=θk时,xn≠0,否则xn=0,n∈{1,2,…,N}. 于是在稀疏表征框架下,DOA估计问题被转化为确定x的非零元素所在位置问题.

基于协方差矩阵的稀疏表征公式(10),文献[7]提出了一种低复杂度的DOA估计算法:

(11)

(p,M2)表示置信度为p、自由度为M2的卡方分布的逆函数.

1.2.2基于Khatri-Rao积的加权稀疏恢复算法

在式(5)等号两边同乘以矩阵GT可得

GTz=GTGBp+σ2GTvec(IM).

(12)

注意到矩阵G为列正交矩阵,即有

W≜ GTG=diag{1,2,…,M-1,M,

M-1,…,2,1}.

(13)

由于矩阵W可逆,所以容易得到

W-1GTz=Bp+σ2W-1GTvec(IM).

(14)

注意到

W-1GTvec(IM)=[01×(M-1),1,01×(M-1)]T≜e0,式(14)可以表示为

W-1GTz=Bp+σ2e0.

(15)

令D=W-1GT,注意到D为行满秩矩阵,根据文献[13]可知

≜∈C(2M-1)×(2M-1).

(16)

同样地,根据文献[13]容易得到,γ的估计误差Δγ服从如下的渐进高斯分布

(17)

根据式(17)容易获得

(18)

在稀疏表征框架下,式(15)可以表示为

γ=Φx+σ2e0∈C(2M-1)×1.

(19)

(20)

式中,η2为正则化参数,它折中着拟合误差和解的稀疏性.

为了顺利求解公式(20),需要获得η2的一个合适的估计值. 根据式(18)可知

(21)

根据累积分布函数的定义,进一步可得

(22)

因为累积分布函数为单调递增函数,所以α与η2之间存在着一一对应关系,给定一个非常高的α值,可以确保正则化参数大于拟合残差的概率非常高. 因此,正则化参数可以定义为

η2≜chi2inv(α,2M-1).

(23)

值得注意的是,优化公式(20)中只在约束方程中应用到了采样协方差矩阵. 事实上,还可以通过在目标函数中融入采样协方差矩阵所获得信息的方式来提高算法的估计精度和分辨性能[6,9-10,14],例如基于Capon谱估计技术的加权稀疏恢复算法[6,10,14]. 这些算法将Capon算法视作“粗处理”用以获取信号的稀疏分布信息和相应的权值,即在稀疏恢复进行之前,先进行数据学习,以便于利用学习获得的信息来改善稀疏恢复算法的性能.

为充分利用采样协方差矩阵中蕴涵的信息,在这里我们借鉴基于Capon谱加权稀疏恢复算法的思想,使用Capon谱的倒数作为权值用以引入信号的稀疏分布信息,最终达到提升算法性能之目的.

(24)

稀疏信号x的第n个元素xn上的权值wn定义为

(25)

结合优化公式(20)和上述的加权策略,可以得到基于Khatri-Rao积的加权l1稀疏恢复算法(KRWl1算法)为

(26)

1.2.3运算复杂度分析

表1　运算量对比

2仿真实验

为了验证本文所提KRWl1算法的性能,本节将从估计精度、分辨性能和运算量等三个方面对其性能进行考察,并将其与经典的l1-SVD[4]、WGMF[6]和He的算法[7]进行比较,使用CVX工具包[15]求解上述四种算法. 实验中如非特别声明,阵元个数为M=8,阵元之间的间隔d=0.5λ,对观察角度空间[0,π]以0.1°的间隔进行均匀采样,即N=1 801. 在计算均方根误差(Root Mean Squared Error,RMSE)时,使用了500次蒙特卡洛实验,故RMSE的计算公式为

(27)

2.1RMSE

在计算RMSE时,3个不相关信号源分别位于{12°,23°,35°}. 图1和图2分别给出了信噪比和快拍数变化情况下的RMSE. 图1中快拍数固定在100,图2中信噪比固定在5 dB. 根据图1和图2可知,相比于同样使用了Capon倒谱加权策略的稀疏恢复算法——WGMF算法,本文提出的KRWl1算法拥有更好的估计精度,特别是当快拍数大于20时,因为所提算法虚拟地扩展了阵列的孔径. 与同样使用Khatri-Rao积的He’s method[7]相比,KRWl1算法在目标函数中引入了信号的稀疏分布信息,稀疏分布信息的应用使得稀疏恢复更倾向于在那些真实目标位置处安排非零元素,所以能够获得更好的估计精度.

图1　信噪比变化情况下的RMSE曲线

图2　快拍数变化情况下的RMSE曲线

2.2分辨性能

图3给出了算法的分辨性能曲线,该图中信噪比为5 dB,快拍数为100. 根据图3可知,在本文给定条件下,所提KRWl1算法的分辨极限为0.9倍的瑞利限,He’s method的分辨力与瑞利限相当,而其他两种算法在给定条件下均无法突破瑞利限. KRWl1算法和He’s method具有较好分辨力的原因在于它们均采用了Khatri-Rao积扩展了阵列孔径. 前者分辨力更好的原因是采用加权l1最小化之后能够更好地逼近理想的l0最小化问题[10,14]. Candes等人[18]指出,l0最小化有唯一精确解所需的观测数要小于l1最小化有唯一精确解对观测数的要求,对于DOA估计问题,观测数意指阵元个数. 因此,根据Candes等人的结论可以得出当阵元个数相同时,加权l1最小化的分辨力更高.

图3　分辨性能曲线线

2.3运算量对比

本例主要考察提出算法的运算复杂度．仿真中,PC采用英特尔CPU,主频3.4 GHz,运行软件为Matlab2012b,运行时间来自CVX 2.0 beta的求解结果．仿真条件如下:M=12,N=1 081,K=3,信噪比为10 dB,T=100,实验次数为500．如表2所示,KRWl1算法的平均总运行时间大幅减少,这与理论分析一致,原因在于KRWl1算法有效地减少了二阶锥的维数．

表2　平均运行时间

3结论

本文在等距线阵及信号不相关的条件下,提出了一种低复杂度的稀疏恢复算法. 基于Khatri-Rao积获得了协方差矩阵的稀疏表征形式,该表征形式可以将MMV问题转化为SMV问题,这会有效地减少未知数的个数. 在此基础上,所提算法使用线性变换来进一步降低二阶锥规划的维数,从而导致所提KRWl1算法具有较低的运算复杂度. 此外算法通过使用加权l1最小化策略来提升算法的估计性能. 仿真实验表明,所提KRWl1算法不仅具有较低的运算量,而且具有较好的估计精度和分辨性能.

参考文献

[1] 陈显舟,杨源,韩静静,等. 双基地多入多出雷达收发方位角联合估计算法[J]. 电波科学学报, 2013, 28(1):176-182.

CHEN Xianzhou, YANG Yuan, HAN Jingjing, et al. Joint DOD and DOA estimation using polynomial rooting for bistatic MIMO radar [J]. Chinese Journal of Radio Science, 2013, 28(1):176-182. (in Chinese)

[2] 程院兵,顾红,苏卫民. 双基地MIMO雷达发射波束形成与多目标定位[J]. 电波科学学报,2012, 27(2): 275-281.

CHENG Yuanbing, GU Hong, SU Weimin. Transmit beamforming and multi-target localization in bistatic MIMO radar [J]. Chinese Journal of Radio Science, 2012, 27(2): 275-281. (in Chinese)

[3] 王永良,陈辉,彭应宁,等.空间谱估计理论与算法[M]. 北京:清华大学出版社,2004.

WANG Yongliang, CHEN Hui, PENG Yingning, et al. Spatial Spectrum Estimation Theory and Algorithm [M]. Beijing: Tsinghua University Press, 2004. (in Chinese)

[4] MALIOUTOV D, CETIN M, WILLSKY A. A sparse signal reconstruction perspective for source localization with sensor arrays [J]. IEEE Transactions on Signal Processing, 2005, 53(8):3010-3022.

[5] YIN J, CHEN T. Direction-of-arrival estimation using a sparse representation of array covariance vectors [J]. IEEE Transactions on Signal Processing, 2011, 59(9):4489-4493.

[6] XU X, WEI X, YE Z. DOA estimation based on sparse signal recovery utilizing weightedl1norm penalty [J]. IEEE Signal Processing Letters, 2012, 19(3):155-158.

[7] HE Z Q, LIU Q H, JIN L N, et al. Low complexity method for DOA estimation using array covariance matrix sparse representation [J]. Electronics Letters, 2013, 49(3): 228-230.

[8] ZHENG C, LI G, LIU Y, et al. Subspace weightedl2,1minimization for sparse signal recovery [J]. EURASIP Journal on Advances in Signal Processing, 2012, 2012:98.

[9] ZHENG Chundi, LI Gang, WANG Xiqin. Combination of weightedl2,1minimization with unitary transformation for DOA estimation [J]. Signal Processing, 2013, 93(12): 3430-3434.

[10]STOICA P, BABU P, LI J. New method of sparse parameter estimation in separable models and its use for spectral analysis of irregularly sampled data [J]. IEEE Transactions on Signal Processing, 2011, 59(1):35-47.

[11]NESTEROV Y, NEMIROVSKII A, YE Y. Interior-point polynomial algorithms in convex programming [M]. SIAM Philadelphia, 1994.

[12]MA W K, HSIEH T H, CHI C Y. Direction-of-arrival estimation of quasi-stationary signals with less sensors than sources and unknown spatial noise covariance: a Khatri-Rao subspace approach [J]. IEEE Trans. Signal Processing, 2010, 58, (4): 2168-2180.

[13]OTTERSTEN B, STOICA P, ROY R. Covariance matching estimation techniques for array signal processing applications [J]. Digital Signal Processing, 1998, 8(3):185-210.

[14]ZHENG C, LI G, XIA X, et al. Weighted l2,1minimization for high resolution range profile with stepped frequency radar [J]. Electronics Letters, 2012, 48(18):1155-1156.

[15]GRANT M, BOYD S. CVX: Matlab software for disciplined convex programming, version 1.22[EB/OL]2012[2012-03-04]. http://cvxr.com/cvx/.

[16]FUCHS J. DOA estimation in the presence of modeling errors, the global matched filter approach[C]//Proceedings of 17th European Signal Processing Conference.Glasgow, 2009: 1963-1967.

[17]JOHNSON D, DUDGEON D. Array Signal Processing: Concepts and Techniques [M]. Englewood Cliffs: Prentice-Hall, 1993.

[18]CANDES E, WAKIN M, BOYD S. Enhancing sparsity by reweightedl1minimization [J]. Journal of Fourier Analysis and Applications, 2008, 14(5):877-905.

段素馨(1983-),男,河南人,中国电子设备系统工程公司无线电管理部工程师,清华大学电子工程系硕士研究生,研究方向包括电磁频谱管理、辐射源定位、数字信号处理等.

张颢(1972-),男,河北人,清华大学电子工程系副教授,研究方向包括雷达信号处理、稀疏信号分析、阵列信号处理等.

孙秀志(1965-),男,山东人,中国电子设备系统工程公司无线电管理部高级工程师,研究方向包括电磁频谱管理、数字信号处理等.

郑春弟(1979-),男,陕西人,海军陆战学院副教授,主要研究方向为稀疏恢复、阵列信号处理、雷达成像.

许勇, 黄勇, 周铸. 双时间步时域有限体积方法计算时变电磁场[J]. 电波科学学报，2015，30(4)：647-652. doi:10.13443/j.cjors. 2014062601

XU Yong, HUANG Yong, ZHOU Zhu. Dual time stepping FVTD method for computation of time-dependent electromagnetic fields [J]. Chinese Journal of Radio Science，2015，30(4)：647-652. (in Chinese). doi:10.13443/j.cjors. 2014062601