基于加權(quán)l(xiāng)1最小化的低復(fù)雜度波達(dá)方向估計(jì)算法

段素馨1，2張顥1孫秀志2鄭春弟3

(1.清華大學(xué)電子工程系,北京 100084;2.中國電子設(shè)備系統(tǒng)工程公司無線電管理部,北京 100840;3.海軍陸戰(zhàn)學(xué)院,廣東 廣州 510430)

摘要基于陣列協(xié)方差矩陣的稀疏表征和陣列響應(yīng)矩陣的Khatri-Rao積,提出了一種低運(yùn)算復(fù)雜度的波達(dá)方向估計(jì)算法. 所提算法在減少未知數(shù)個(gè)數(shù)的同時(shí),通過線性變換降低約束方程的維數(shù),可有效減少優(yōu)化問題的計(jì)算復(fù)雜度. 為充分利用陣列協(xié)方差矩陣中蘊(yùn)涵的信息,使用Capon譜的倒數(shù)作為權(quán)值構(gòu)建出了加權(quán)l(xiāng)1最小化問題,這使得所提算法在降低運(yùn)算量的同時(shí)能夠獲得較好的估計(jì)性能. 仿真實(shí)驗(yàn)驗(yàn)證了所提算法的有效性.

關(guān)鍵詞波達(dá)方向估計(jì);加權(quán)l(xiāng)1最小化;稀疏恢復(fù);等距線陣

中圖分類號(hào)TN911.7

文獻(xiàn)標(biāo)志碼A

文章編號(hào)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)

資助項(xiàng)目: 國家自然科學(xué)基金(61401496)

聯(lián)系人： 段素馨 E-mail：007-yx@163.com

引言

作為陣列信號(hào)處理領(lǐng)域中的基本問題,波達(dá)方向(Direction of Arrival,DOA)估計(jì)在電磁學(xué)、聲學(xué)、水聲學(xué)、地震學(xué)和生物醫(yī)學(xué)等領(lǐng)域有著廣泛的應(yīng)用[1-3]. 對于DOA估計(jì)算法而言,高分辨性能始終是研究的重點(diǎn). 近年來,利用觀測角度空間中潛在目標(biāo)具有稀疏性這一信息,眾多基于稀疏表征的DOA估計(jì)算法被提出[4-10],其中采用基追蹤(Basis Pursuit,BP)準(zhǔn)則的DOA估計(jì)算法由于求解方便且精度高成為最重要的研究分支[4-9]. 相比于傳統(tǒng)的多重信號(hào)分類(Multiple Signal Classification, MUSIC)、子空間旋轉(zhuǎn)不變技術(shù)(Estimating Signal Parameters via Rotational Invariance Techniques, ESPRIT)等子空間類DOA估計(jì)算法,基于稀疏表征的DOA估計(jì)算法最大的優(yōu)勢之一就是可以獲得更高的分辨性能[4].

基追蹤可用二階錐規(guī)劃(Second Order Cone Programming,SOCP)進(jìn)行求解. 采用內(nèi)點(diǎn)法的SOCP所需的運(yùn)算量為O((N×T)3)[4],而子空間類算法的運(yùn)算復(fù)雜度為O(M3),其中N、T和M分別為稀疏信號(hào)長度、快拍數(shù)和陣元個(gè)數(shù),M?N,因此采用BP準(zhǔn)則的DOA估計(jì)算法的運(yùn)算復(fù)雜度一般要遠(yuǎn)高于傳統(tǒng)的DOA估計(jì)算法. 為了降低基于稀疏表征的DOA估計(jì)算法的運(yùn)算復(fù)雜度,經(jīng)典的l1-SVD算法[4]采用主分量分析(Principal Component Analysis,PCA)降低稀疏恢復(fù)的運(yùn)算量,該算法使用截尾奇異值分解(Truncated Singular Value Decomposition, TSVD)將未知數(shù)的個(gè)數(shù)從N×T減少到N×K,其中K為信號(hào)個(gè)數(shù),通常K≤T. 文獻(xiàn)[7]在等距線陣(Uniformly-spaced Linear Array, ULA)條件下,將接收數(shù)據(jù)協(xié)方差矩陣進(jìn)行稀疏表征,提出了一種低復(fù)雜度的稀疏恢復(fù)算法(本文稱之為He’s method),該算法使用Khatri-Rao積在虛擬擴(kuò)展陣列孔徑的同時(shí),將多觀測向量(Multiple Measurements Vectors,MMV)稀疏恢復(fù)問題轉(zhuǎn)換為單觀測向量(Single Measurements Vectors,SMV),從而使得未知數(shù)的個(gè)數(shù)降低為N. 總的來說這兩種算法均是從減少未知數(shù)個(gè)數(shù)這個(gè)角度來降低稀疏恢復(fù)的運(yùn)算復(fù)雜度.

事實(shí)上,SOCP的運(yùn)算復(fù)雜度還與二階錐的維數(shù)有關(guān)[11],與上述兩種方法不同,本文試圖通過減少二階錐的維數(shù)(二階錐C={x=(x0,x1)∈R×Rk-1:x0-‖x1‖≥0}的維數(shù),即k)來實(shí)現(xiàn)降低運(yùn)算復(fù)雜度的目的. 與He’smethod方法相似,我們對接收數(shù)據(jù)協(xié)方差矩陣進(jìn)行稀疏表征. 但是在使用Khatri-Rao積實(shí)現(xiàn)孔徑虛擬擴(kuò)展之后,采用線性變換將二階錐的維數(shù)從M2降為2M-1. 另外,為了充分利用協(xié)方差矩陣中蘊(yùn)涵的信息,本文使用Capon譜的倒數(shù)對目標(biāo)函數(shù)進(jìn)行加權(quán)處理[6,9-10],通過加權(quán)處理之后,恢復(fù)算法更傾向于在那些更有可能是真實(shí)DOA的位置安排非零元素,從而可以有效提升算法的性能.

1理論分析

1.1信號(hào)模型

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

(1)

式中: y(t)=[y1(t),…,yM(t)]T∈CM×1為M個(gè)陣元上的接收數(shù)據(jù)向量; n(t)=[n1(t),…,nM(t)]T∈CM×1為M個(gè)陣元上的噪聲向量; A=[a(θ1),…,a(θK)]∈CM×K為陣列響應(yīng)矩陣,a(θk)=[1,ej2πdsin(θk)/λ ,

…,ej(M-1)2πdsin(θk)/λ]T∈CM×1為陣列的導(dǎo)向矢量. 假設(shè)噪聲n(t)為均值為0、方差為σ2的高斯白噪聲. 不失一般性,假設(shè)n(t)與s(t)不相關(guān). 假設(shè)s(t)為不相關(guān)信號(hào),則式(1)的協(xié)方差矩陣可以表示為

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 ;符號(hào)“?”、“*”分別代表矩陣的Kronecher積和Khatri-Rao積; p=[p1,…,pK]T.

根據(jù)Khatri-Rao積的特點(diǎn)和ULA的假設(shè),容易得到[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為虛擬的陣列響應(yīng)矩陣,相應(yīng)的虛擬導(dǎo)向矢量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. 相比初始的導(dǎo)向矢量a(θk),陣列孔徑虛擬地從M擴(kuò)展為2M-1,這是使用Khatri-Rao積帶來的最大好處,陣列孔徑的擴(kuò)展將使得DOA估計(jì)算法的性能獲得極大的改善. 將式(4)代入式(3)可得

z=GBp+σ2vec(IM).

(5)

(6)

令

(7)

根據(jù)文獻(xiàn)[13]中的結(jié)論,有

(8)

文獻(xiàn)[13]中指出,協(xié)方差矩陣的估計(jì)誤差服從如下的高斯分布

(9)

式中,符號(hào)AsN表示漸進(jìn)高斯分布．

1.2基于Khatri-Rao積的加權(quán)稀疏恢復(fù)算法

本節(jié)首先簡要回顧文獻(xiàn)[7]中提出的低復(fù)雜度稀疏恢復(fù)算法,在分析該算法運(yùn)算量的基礎(chǔ)上,尋求進(jìn)一步降低運(yùn)算量的途徑,并將其應(yīng)用于稀疏恢復(fù),提出基于Khatri-Rao積的加權(quán)稀疏恢復(fù)算法.

1.2.1基于Khatri-Rao積的稀疏恢復(fù)算法[7]

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

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

(10)

式中: Φ=[b(φ1),…,b(φN)]∈CM×N是超完備基矩陣,即有M?N,集合{φ1,…,φN}表示觀測角度空間[0,π]被均勻劃分為N個(gè)格點(diǎn)所獲得的角度集合; x∈RN×1表示稀疏信號(hào),當(dāng)且僅當(dāng)φn=θk時(shí),xn≠0,否則xn=0,n∈{1,2,…,N}. 于是在稀疏表征框架下,DOA估計(jì)問題被轉(zhuǎn)化為確定x的非零元素所在位置問題.

基于協(xié)方差矩陣的稀疏表征公式(10),文獻(xiàn)[7]提出了一種低復(fù)雜度的DOA估計(jì)算法:

(11)

(p,M2)表示置信度為p、自由度為M2的卡方分布的逆函數(shù).

1.2.2基于Khatri-Rao積的加權(quán)稀疏恢復(fù)算法

在式(5)等號(hào)兩邊同乘以矩陣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為行滿秩矩陣,根據(jù)文獻(xiàn)[13]可知

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

(16)

同樣地,根據(jù)文獻(xiàn)[13]容易得到,γ的估計(jì)誤差Δγ服從如下的漸進(jìn)高斯分布

(17)

根據(jù)式(17)容易獲得

(18)

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

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

(19)

(20)

式中,η2為正則化參數(shù),它折中著擬合誤差和解的稀疏性.

為了順利求解公式(20),需要獲得η2的一個(gè)合適的估計(jì)值. 根據(jù)式(18)可知

(21)

根據(jù)累積分布函數(shù)的定義,進(jìn)一步可得

(22)

因?yàn)槔鄯e分布函數(shù)為單調(diào)遞增函數(shù),所以α與η2之間存在著一一對應(yīng)關(guān)系,給定一個(gè)非常高的α值,可以確保正則化參數(shù)大于擬合殘差的概率非常高. 因此,正則化參數(shù)可以定義為

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

(23)

值得注意的是,優(yōu)化公式(20)中只在約束方程中應(yīng)用到了采樣協(xié)方差矩陣. 事實(shí)上,還可以通過在目標(biāo)函數(shù)中融入采樣協(xié)方差矩陣所獲得信息的方式來提高算法的估計(jì)精度和分辨性能[6,9-10,14],例如基于Capon譜估計(jì)技術(shù)的加權(quán)稀疏恢復(fù)算法[6,10,14]. 這些算法將Capon算法視作“粗處理”用以獲取信號(hào)的稀疏分布信息和相應(yīng)的權(quán)值,即在稀疏恢復(fù)進(jìn)行之前,先進(jìn)行數(shù)據(jù)學(xué)習(xí),以便于利用學(xué)習(xí)獲得的信息來改善稀疏恢復(fù)算法的性能.

為充分利用采樣協(xié)方差矩陣中蘊(yùn)涵的信息,在這里我們借鑒基于Capon譜加權(quán)稀疏恢復(fù)算法的思想,使用Capon譜的倒數(shù)作為權(quán)值用以引入信號(hào)的稀疏分布信息,最終達(dá)到提升算法性能之目的.

(24)

稀疏信號(hào)x的第n個(gè)元素xn上的權(quán)值wn定義為

(25)

結(jié)合優(yōu)化公式(20)和上述的加權(quán)策略,可以得到基于Khatri-Rao積的加權(quán)l(xiāng)1稀疏恢復(fù)算法(KRWl1算法)為

(26)

1.2.3運(yùn)算復(fù)雜度分析

表1　運(yùn)算量對比

2仿真實(shí)驗(yàn)

為了驗(yàn)證本文所提KRWl1算法的性能,本節(jié)將從估計(jì)精度、分辨性能和運(yùn)算量等三個(gè)方面對其性能進(jìn)行考察,并將其與經(jīng)典的l1-SVD[4]、WGMF[6]和He的算法[7]進(jìn)行比較,使用CVX工具包[15]求解上述四種算法. 實(shí)驗(yàn)中如非特別聲明,陣元個(gè)數(shù)為M=8,陣元之間的間隔d=0.5λ,對觀察角度空間[0,π]以0.1°的間隔進(jìn)行均勻采樣,即N=1 801. 在計(jì)算均方根誤差(Root Mean Squared Error,RMSE)時(shí),使用了500次蒙特卡洛實(shí)驗(yàn),故RMSE的計(jì)算公式為

(27)

2.1RMSE

在計(jì)算RMSE時(shí),3個(gè)不相關(guān)信號(hào)源分別位于{12°,23°,35°}. 圖1和圖2分別給出了信噪比和快拍數(shù)變化情況下的RMSE. 圖1中快拍數(shù)固定在100,圖2中信噪比固定在5 dB. 根據(jù)圖1和圖2可知,相比于同樣使用了Capon倒譜加權(quán)策略的稀疏恢復(fù)算法——WGMF算法,本文提出的KRWl1算法擁有更好的估計(jì)精度,特別是當(dāng)快拍數(shù)大于20時(shí),因?yàn)樗崴惴ㄌ摂M地?cái)U(kuò)展了陣列的孔徑. 與同樣使用Khatri-Rao積的He’s method[7]相比,KRWl1算法在目標(biāo)函數(shù)中引入了信號(hào)的稀疏分布信息,稀疏分布信息的應(yīng)用使得稀疏恢復(fù)更傾向于在那些真實(shí)目標(biāo)位置處安排非零元素,所以能夠獲得更好的估計(jì)精度.

圖1　信噪比變化情況下的RMSE曲線

圖2　快拍數(shù)變化情況下的RMSE曲線

2.2分辨性能

圖3給出了算法的分辨性能曲線,該圖中信噪比為5 dB,快拍數(shù)為100. 根據(jù)圖3可知,在本文給定條件下,所提KRWl1算法的分辨極限為0.9倍的瑞利限,He’s method的分辨力與瑞利限相當(dāng),而其他兩種算法在給定條件下均無法突破瑞利限. KRWl1算法和He’s method具有較好分辨力的原因在于它們均采用了Khatri-Rao積擴(kuò)展了陣列孔徑. 前者分辨力更好的原因是采用加權(quán)l(xiāng)1最小化之后能夠更好地逼近理想的l0最小化問題[10,14]. Candes等人[18]指出,l0最小化有唯一精確解所需的觀測數(shù)要小于l1最小化有唯一精確解對觀測數(shù)的要求,對于DOA估計(jì)問題,觀測數(shù)意指陣元個(gè)數(shù). 因此,根據(jù)Candes等人的結(jié)論可以得出當(dāng)陣元個(gè)數(shù)相同時(shí),加權(quán)l(xiāng)1最小化的分辨力更高.

圖3　分辨性能曲線線

2.3運(yùn)算量對比

本例主要考察提出算法的運(yùn)算復(fù)雜度．仿真中,PC采用英特爾CPU,主頻3.4 GHz,運(yùn)行軟件為Matlab2012b,運(yùn)行時(shí)間來自CVX 2.0 beta的求解結(jié)果．仿真條件如下:M=12,N=1 081,K=3,信噪比為10 dB,T=100,實(shí)驗(yàn)次數(shù)為500．如表2所示,KRWl1算法的平均總運(yùn)行時(shí)間大幅減少,這與理論分析一致,原因在于KRWl1算法有效地減少了二階錐的維數(shù)．

表2　平均運(yùn)行時(shí)間

3結(jié)論

本文在等距線陣及信號(hào)不相關(guān)的條件下,提出了一種低復(fù)雜度的稀疏恢復(fù)算法. 基于Khatri-Rao積獲得了協(xié)方差矩陣的稀疏表征形式,該表征形式可以將MMV問題轉(zhuǎn)化為SMV問題,這會(huì)有效地減少未知數(shù)的個(gè)數(shù). 在此基礎(chǔ)上,所提算法使用線性變換來進(jìn)一步降低二階錐規(guī)劃的維數(shù),從而導(dǎo)致所提KRWl1算法具有較低的運(yùn)算復(fù)雜度. 此外算法通過使用加權(quán)l(xiāng)1最小化策略來提升算法的估計(jì)性能. 仿真實(shí)驗(yàn)表明,所提KRWl1算法不僅具有較低的運(yùn)算量,而且具有較好的估計(jì)精度和分辨性能.

參考文獻(xiàn)

[1] 陳顯舟,楊源,韓靜靜,等. 雙基地多入多出雷達(dá)收發(fā)方位角聯(lián)合估計(jì)算法[J]. 電波科學(xué)學(xué)報(bào), 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] 程院兵,顧紅,蘇衛(wèi)民. 雙基地MIMO雷達(dá)發(fā)射波束形成與多目標(biāo)定位[J]. 電波科學(xué)學(xué)報(bào),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] 王永良,陳輝,彭應(yīng)寧,等.空間譜估計(jì)理論與算法[M]. 北京:清華大學(xué)出版社,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-),男,河南人,中國電子設(shè)備系統(tǒng)工程公司無線電管理部工程師,清華大學(xué)電子工程系碩士研究生,研究方向包括電磁頻譜管理、輻射源定位、數(shù)字信號(hào)處理等.

張顥(1972-),男,河北人,清華大學(xué)電子工程系副教授,研究方向包括雷達(dá)信號(hào)處理、稀疏信號(hào)分析、陣列信號(hào)處理等.

孫秀志(1965-),男,山東人,中國電子設(shè)備系統(tǒng)工程公司無線電管理部高級(jí)工程師,研究方向包括電磁頻譜管理、數(shù)字信號(hào)處理等.

鄭春弟(1979-),男,陜西人,海軍陸戰(zhàn)學(xué)院副教授,主要研究方向?yàn)橄∈杌謴?fù)、陣列信號(hào)處理、雷達(dá)成像.

許勇, 黃勇, 周鑄. 雙時(shí)間步時(shí)域有限體積方法計(jì)算時(shí)變電磁場[J]. 電波科學(xué)學(xué)報(bào)，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