基于加權l(xiāng) 1最小化的低復雜度波達方向估計算法

基于加權l(xiāng)1最小化的低復雜度波達方向估計算法

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

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

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

關鍵詞波達方向估計;加權l(xiāng)1最小化;稀疏恢復;等距線陣

中圖分類號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)

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

引言

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

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

事實上,SOCP的運算復雜度還與二階錐的維數(shù)有關[11],與上述兩種方法不同,本文試圖通過減少二階錐的維數(shù)(二階錐C={x=(x0,x1)∈R×Rk-1:x0-‖x1‖≥0}的維數(shù),即k)來實現(xiàn)降低運算復雜度的目的. 與He’smethod方法相似,我們對接收數(shù)據(jù)協(xié)方差矩陣進行稀疏表征. 但是在使用Khatri-Rao積實現(xiàn)孔徑虛擬擴展之后,采用線性變換將二階錐的維數(shù)從M2降為2M-1. 另外,為了充分利用協(xié)方差矩陣中蘊涵的信息,本文使用Capon譜的倒數(shù)對目標函數(shù)進行加權處理[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個陣元上的接收數(shù)據(jù)向量; 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)的協(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 ;符號“?”、“*”分別代表矩陣的Kronecher積和Khatri-Rao積; p=[p1,…,pK]T.

根據(jù)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)

根據(jù)文獻[13]中的結論,有

(8)

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

(9)

式中,符號AsN表示漸進高斯分布．

1.2基于Khatri-Rao積的加權稀疏恢復算法

本節(jié)首先簡要回顧文獻[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的非零元素所在位置問題.

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

(11)

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

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為行滿秩矩陣,根據(jù)文獻[13]可知

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

(16)

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

(17)

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

(18)

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

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

(19)

(20)

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

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

(21)

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

(22)

因為累積分布函數(shù)為單調遞增函數(shù),所以α與η2之間存在著一一對應關系,給定一個非常高的α值,可以確保正則化參數(shù)大于擬合殘差的概率非常高. 因此,正則化參數(shù)可以定義為

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

(23)

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

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

(24)

稀疏信號x的第n個元素xn上的權值wn定義為

(25)

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

(26)

1.2.3運算復雜度分析

表1　運算量對比

2仿真實驗

為了驗證本文所提KRWl1算法的性能,本節(jié)將從估計精度、分辨性能和運算量等三個方面對其性能進行考察,并將其與經(jīng)典的l1-SVD[4]、WGMF[6]和He的算法[7]進行比較,使用CVX工具包[15]求解上述四種算法. 實驗中如非特別聲明,陣元個數(shù)為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分別給出了信噪比和快拍數(shù)變化情況下的RMSE. 圖1中快拍數(shù)固定在100,圖2中信噪比固定在5 dB. 根據(jù)圖1和圖2可知,相比于同樣使用了Capon倒譜加權策略的稀疏恢復算法——WGMF算法,本文提出的KRWl1算法擁有更好的估計精度,特別是當快拍數(shù)大于20時,因為所提算法虛擬地擴展了陣列的孔徑. 與同樣使用Khatri-Rao積的He’s method[7]相比,KRWl1算法在目標函數(shù)中引入了信號的稀疏分布信息,稀疏分布信息的應用使得稀疏恢復更傾向于在那些真實目標位置處安排非零元素,所以能夠獲得更好的估計精度.

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

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

2.2分辨性能

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

圖3　分辨性能曲線線

2.3運算量對比

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

表2　平均運行時間

3結論

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

參考文獻

[1] 陳顯舟,楊源,韓靜靜,等. 雙基地多入多出雷達收發(fā)方位角聯(lián)合估計算法[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] 程院兵,顧紅,蘇衛(wèi)民. 雙基地MIMO雷達發(fā)射波束形成與多目標定位[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-),男,河南人,中國電子設備系統(tǒng)工程公司無線電管理部工程師,清華大學電子工程系碩士研究生,研究方向包括電磁頻譜管理、輻射源定位、數(shù)字信號處理等.

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

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

鄭春弟(1979-),男,陜西人,海軍陸戰(zhàn)學院副教授,主要研究方向為稀疏恢復、陣列信號處理、雷達成像.

許勇, 黃勇, 周鑄. 雙時間步時域有限體積方法計算時變電磁場[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