基于壓縮感知的稀疏重構(gòu)DOA估計(jì)算法

包曉蕾，曲行根，王卓英

(1.上海電子信息職業(yè)技術(shù)學(xué)院 通信與信息工程系，上海 201411；2.哈爾濱工程大學(xué) 信息與通信工程學(xué)院，黑龍江 哈爾濱 150001)

?

基于壓縮感知的稀疏重構(gòu)DOA估計(jì)算法

包曉蕾1，曲行根2，王卓英1

(1.上海電子信息職業(yè)技術(shù)學(xué)院 通信與信息工程系，上海 201411；2.哈爾濱工程大學(xué) 信息與通信工程學(xué)院，黑龍江 哈爾濱 150001)

基于空間目標(biāo)分布的稀疏特性和壓縮感知理論思想，提出一種基于奇異值分解的多測(cè)量梯度投影稀疏重構(gòu)(SVD-MGPSR)算法，將多目標(biāo)DOA估計(jì)轉(zhuǎn)化為一個(gè)稀疏信號(hào)重構(gòu)問(wèn)題。首先利用陣列流形建立的過(guò)完備原子庫(kù)對(duì)信號(hào)進(jìn)行聯(lián)合稀疏表示，然后對(duì)壓縮采樣后的信息矩陣進(jìn)行奇異值分解，可以明顯降低運(yùn)算量，最后基于MGPSR算法對(duì)稀疏信號(hào)進(jìn)行重構(gòu)，從而實(shí)現(xiàn)DOA估計(jì)。相對(duì)于已有算法，該算法不僅在低信噪比、小快拍數(shù)條件下測(cè)向均方誤差較小，而且能夠?qū)ο喔尚盘?hào)進(jìn)行正確估計(jì)，具有較高的測(cè)向精度和角度分辨率。仿真實(shí)驗(yàn)驗(yàn)證了該算法的有效性。

壓縮感知；DOA估計(jì)；奇異值分解；梯度投影

目標(biāo)信號(hào)波達(dá)方向(direction of arrival, DOA)估計(jì)一直是陣列信號(hào)處理領(lǐng)域中的一個(gè)重要研究方向，在通信、雷達(dá)和醫(yī)學(xué)成像等領(lǐng)域有著廣泛的應(yīng)用[1]。目前，國(guó)內(nèi)外學(xué)者已經(jīng)提出了大量有效的DOA估計(jì)算法，其主要分為兩大類：一類是基于波束形成的測(cè)向方法，如由BURG提出的最大熵譜估計(jì)法(maximum entropy spectral estimation method, MEM)[2]和CAPON提出的最小方差無(wú)失真響應(yīng)估計(jì)法(minimum variance distortion response, MVDR)[3]等；另一類基于子空間的測(cè)向方法，如由SCHIMIDT提出的多重信號(hào)分類(multiple signal classification, MUSIC)算法[4]和ROY等提出的旋轉(zhuǎn)不變子空間(estimation of signal parameters via rotational invariance techniques, ESPRIT)算法[5]等。然而，基于波束形成的算法雖然克服了陣列孔徑限制，但估計(jì)分辨率不高且受噪聲影響較大；基于子空間的算法具有高分辨率，但對(duì)于相干信號(hào)由于信號(hào)協(xié)方差矩陣會(huì)出現(xiàn)失秩現(xiàn)象，故無(wú)法有效估計(jì)。

近幾年一種新穎的壓縮感知(compressing sensing, CS)[6]理論已廣泛應(yīng)用于圖像分析、信號(hào)處理等領(lǐng)域，受到了學(xué)者們極大關(guān)注。CS理論是指可以用遠(yuǎn)低于奈奎斯特采樣定理要求的采樣速率對(duì)可壓縮或可稀疏信號(hào)進(jìn)行非自適應(yīng)采樣，然后通過(guò)求解凸優(yōu)化問(wèn)題恢復(fù)信號(hào)。

由于可以將空間進(jìn)行網(wǎng)格劃分并且實(shí)際目標(biāo)角度僅占較少單元，故目標(biāo)在空間域是稀疏的。基于上述思想，可以將壓縮感知應(yīng)用到DOA估計(jì)中[7-8]，文獻(xiàn)[9]提出一種利用混合l2,0范數(shù)(JLZA-DOA)優(yōu)化算法來(lái)近似求解聯(lián)合稀疏問(wèn)題，從而對(duì)信號(hào)進(jìn)行DOA估計(jì)。文獻(xiàn)[10]采用梯度投影稀疏重構(gòu)(gradient projection for sparse reconstruction, GPSR)算法重構(gòu)信號(hào)，該算法梯度約束條件過(guò)于嚴(yán)格并且只適應(yīng)于單快拍系統(tǒng)，在低信噪比條件下均方誤差較大。針對(duì)上述問(wèn)題，筆者提出一種新的基于奇異值分解(singular value decomposition, SVD)的多測(cè)量梯度投影稀疏重構(gòu)(SVD-MGPSR)算法。該算法首先基于陣列流型矩陣建立的過(guò)完備原子庫(kù)構(gòu)造稀疏基矩陣，構(gòu)造稀疏信號(hào)模型，采用滿足約束等容條件(restricted isometry property, RIP)[11]的隨機(jī)矩陣進(jìn)行時(shí)域壓縮采樣，然后對(duì)壓縮采樣后的信息矩陣進(jìn)行奇異值分解，最后采用新算法重構(gòu)稀疏信號(hào)矩陣，完成DOA估計(jì)。筆者采用的奇異值分解方法進(jìn)行信號(hào)降維處理可以得到低維信號(hào)子空間，運(yùn)算量明顯降低并且能夠提高抗噪聲能力。與已有的MUSIC算法、MVDR算法和文獻(xiàn)[9]中提到的JLZA-DOA算法相比，該算法不僅能夠?qū)ο喔尚盘?hào)正確測(cè)向，而且在低信噪比、小快拍數(shù)條件下測(cè)向均方誤差較小，具有較高的角度分辨力及更優(yōu)的測(cè)向性能。

1 壓縮感知

CS理論原理：信號(hào)若是可壓縮性的或稀疏性的，可以通過(guò)與變換基不相關(guān)的測(cè)量矩陣將高維信號(hào)投影到低維信號(hào)上，再求解凸優(yōu)化問(wèn)題，從少量投影值中重構(gòu)原始信號(hào)。其由3個(gè)核心部分構(gòu)成：信號(hào)稀疏表示、測(cè)量矩陣的設(shè)計(jì)和信號(hào)重構(gòu)算法[12]。

(1)

其中，α=[α1,α2,…,αN]T是信號(hào)x在正交基矩陣Ψ的表示系數(shù)，因此α和x是信號(hào)的等價(jià)表示。若系數(shù)矢量α只含有K(K

在CS理論中，信號(hào)稀疏表示之后，將信號(hào)壓縮投影到M×N的測(cè)量矩陣Φ上得到投影值y=[y(1),y(2),…,y(M)]T，其中M

y=Φx=ΦΨα=Θα

(2)

其中,Θ=ΦΨ是M×N的感知矩陣，也稱之為CS的信息算子[13]。CS理論的關(guān)鍵部分就是從這M個(gè)投影值中重構(gòu)出長(zhǎng)度為N的稀疏信號(hào)或信號(hào)的系數(shù)矢量。求解式(2)是一個(gè)欠定問(wèn)題，一般來(lái)講無(wú)確定解，但研究表明：對(duì)于一個(gè)稀疏信號(hào)，當(dāng)Θ滿足RIP條件時(shí)可以通過(guò)求解l0范數(shù)問(wèn)題精確重構(gòu)原始信號(hào)：

‖α‖l0s.t.y=Θα

(3)

‖α‖l1s.t.y=Θα

(4)

2 基于CS理論的信號(hào)模型

2.1 接收信號(hào)模型

設(shè)各向同性的N個(gè)陣元構(gòu)成一個(gè)均勻線陣，有K個(gè)遠(yuǎn)場(chǎng)信號(hào)sk(t),k=1,2,…,K分別以方向θi,i=1,2，…，K入射到該線陣，則第n個(gè)陣元接收信號(hào)為：

n=1,2，…，N

(5)

式中：nn(t)為第i個(gè)信號(hào)在t時(shí)刻的噪聲；τni為第i個(gè)信號(hào)到達(dá)第n個(gè)陣元相對(duì)于參考陣元的時(shí)延；gni為第n個(gè)陣元對(duì)第i個(gè)信號(hào)的增益。各陣元是同向的，并且陣元之間不考慮通道不一致性、互耦等因素，陣元增益可以歸一化為1。這里令第一個(gè)陣元為參考陣元，則相對(duì)于參考陣元的延時(shí)可以描述為τni=(n-1)dsin(θi)/c，式中d為均勻線陣陣元間隔。假設(shè)陣元接收到的是中心頻率為ω0的窄帶信號(hào)，則第n個(gè)陣元接收信號(hào)為：

nn(t) n=1,2,…,N

(6)

將式(6)寫成矢量形式如下:

X(t)=AS(t)+N(t)

(7)

式中：X(t)=[x1(t),x2(t),…,xN(t)]T為N維接收數(shù)據(jù)矢量；S(t)=[s1(t),s2(t),…,sK(t)]T為K維信號(hào)矢量；N(t)=[n1(t),n2(t),…,nN(t)]T為N維噪聲矢量;A=[a(θ1),a(θ2),…,a(θK)]為N×K維陣列流型矩陣。

2.2 過(guò)完備原子庫(kù)建立

根據(jù)陣列流形結(jié)構(gòu)建立過(guò)完備原子庫(kù)。原子庫(kù)結(jié)構(gòu)需要盡可能地逼近信號(hào)結(jié)構(gòu)，并能用盡可能少的原子最佳線性表示信號(hào)，由此可建立如式(8)所示的過(guò)完備原子庫(kù)的原子：

i=1,2,…,P

(8)

X=ΨZ+N

(9)

其中，Z為信號(hào)S在原子庫(kù)的稀疏系數(shù)矩陣。

2.3 聯(lián)合稀疏模型

系統(tǒng)快拍數(shù)為L(zhǎng)，第i次快拍數(shù)接收的信息矢量為yi，稀疏系數(shù)矢量為zi，噪聲矢量為ni，且滿足均值為零方差為σ2的復(fù)正態(tài)分布，根據(jù)CS理論可得：

Y=ΦX=ΦΨZ+ΦN=ΘZ+ΦN

(10)

3 信號(hào)重構(gòu)算法

3.1 聯(lián)合稀疏模型SVD處理

針對(duì)降低運(yùn)算量和復(fù)雜度的要求，基于子空間思想對(duì)接收到的信息矩陣進(jìn)行SVD降維處理，提取信號(hào)子空間。傳統(tǒng)的子空間算法對(duì)協(xié)方差矩陣進(jìn)行特征分解時(shí)要求信號(hào)不相關(guān)，當(dāng)信號(hào)相關(guān)時(shí)，協(xié)方差矩陣會(huì)出現(xiàn)損秩，降低測(cè)向性能。筆者在不要求信號(hào)是否相關(guān)的條件下，對(duì)式(10)的信息矩陣進(jìn)行降維處理：

Y=UΣVH=[USUN]Σ[VSVN]H

(11)

式中:U和V是正交矩陣；矩陣Σ對(duì)角線元素按照由大到小排列，即σ1≥σ2≥…≥σK≥σK+1≥…≥σM。信號(hào)子空間矩陣US由K個(gè)大奇異值對(duì)應(yīng)的左奇異向量組成，而噪聲子空間UN由其余M-K個(gè)小奇異值對(duì)應(yīng)的左奇異值向量組成。在信源數(shù)已知的條件下進(jìn)一步推導(dǎo)可得M×K維的信號(hào)子空間US=YS=UΣΒK=YVΒK,其中,ΒK=[IK0]H，IK為K×K維單位矩陣，0為K×(L-K)維零矩陣。同時(shí)可得ZS=ZΣΒK，NS=NΣΒK，由此代入式(11)可得降維表示式：

YS=ΘZS+ΦΝS

(12)

由式(12)可知，通過(guò)SVD處理信息，矩陣由M×L維降為M×K維，尤其在大快拍數(shù)條件下，求解式(12)可以顯著降低運(yùn)算時(shí)間，并且由初等列變換得到的ZS與系數(shù)矩陣Z具有相同的稀疏特性，故求解式(12)是具有可行性的。另外，從上述分析可知，SVD分解實(shí)質(zhì)是一個(gè)信號(hào)分量累積的過(guò)程，使得在低信噪比的條件下也具有較好的估計(jì)性能。

3.2 SVD-MGPSR算法

傳統(tǒng)的GPSR算法是利用拉格朗日乘法將式(4)改寫成如下無(wú)約束凸優(yōu)化問(wèn)題：

(13)

其中，τ(>0)稱為拉格朗日因子或者正則因子。引入正則因子的目的是在抑制噪聲干擾和稀疏特性上有效均衡。目標(biāo)函數(shù)中l(wèi)2范數(shù)項(xiàng)表征重構(gòu)信號(hào)與原始信號(hào)能量差值，來(lái)抑制噪聲的影響，而l1范數(shù)項(xiàng)表征重構(gòu)信號(hào)的稀疏度。通過(guò)將zS中的值分為同維數(shù)的正部分u和負(fù)部分v，即ui=(zSi)+、vi=(-zSi)+,i=1,2,…,N，可以把式(13)轉(zhuǎn)化成標(biāo)準(zhǔn)的邊界約束二次規(guī)劃(bound-constrained quadratic programming, BCQP)問(wèn)題：

首先，SVD-MGPSR算法的目標(biāo)函數(shù)如式(15)所示：

(15)

(16)

(17)

根據(jù)矩陣求跡運(yùn)算性質(zhì)可將式(17)進(jìn)一步化簡(jiǎn)為：

(18)

(19)

(20)

4 仿真分析

筆者仿真針對(duì)遠(yuǎn)場(chǎng)窄帶信號(hào)模型，與已有的MUSIC算法、MVDR算法和JLZA-DOA算法進(jìn)行比較，來(lái)驗(yàn)證筆者SVD-MGPSR算法的性能。仿真中，空間信號(hào)隨機(jī)采樣均采用16個(gè)陣元的均勻線陣，陣元間距為半波長(zhǎng)，信號(hào)幅度服從零均值、方差為1的復(fù)高斯變量，且與復(fù)高斯白噪聲互不相關(guān)。SVD-MGPSR算法中測(cè)量矩陣采用隨機(jī)高斯矩陣，且參量τ定義為：

利用均方根誤差作為估計(jì)的性能指標(biāo)，DOA估計(jì)均方根誤差定義為：

實(shí)驗(yàn)1 假設(shè)空間3個(gè)非相干信號(hào)分別以θ1=0°、θ2=49°和θ3=61°方向入射，信噪比為3 dB，快拍數(shù)L=100，比較SVD-MGPSR算法與其他3種算法的歸一化空間譜。圖1給出了4種算法的歸一化空間譜圖，從圖1中可以看出：SVD-MGPSR算法比其他3種算法具有較窄的主瓣寬度和較低的旁瓣，其空間譜類似于針狀，具有較高的角度分辨率，且旁瓣沒(méi)有波動(dòng)現(xiàn)象，噪聲對(duì)新算法測(cè)向性能影響較小。

圖1 非相干信號(hào)空間譜

實(shí)驗(yàn)2 假設(shè)空間入射角度分別為θ1=49°和θ2=61°的兩個(gè)非相干信號(hào)。圖2為各算法在快拍數(shù)L=100和信噪比SNR=-5 dB:10 dB條件下，每個(gè)信噪比做50次蒙特卡洛實(shí)驗(yàn)時(shí)DOA估計(jì)均方根誤差變化情況；圖3為在SNR=15 dB和快拍數(shù)L=25:5:50情況下，每個(gè)快拍數(shù)做50次蒙特卡洛實(shí)驗(yàn)時(shí)4種算法的均方根誤差變化情況。從圖2和圖3可以看出：在低信噪比、小快拍數(shù)條件下SVD-MGPSR算法和JLZA-DOA算法的均方根誤差明顯小于MUSIC算法和MVDR算法，具有較高的測(cè)向精度，能夠保證在低信噪比、小快拍數(shù)條件下的測(cè)向性能；當(dāng)信噪比足夠高、快拍數(shù)足夠大時(shí)各算法均方根誤差相差不大具有相近的測(cè)向性能。

實(shí)驗(yàn)3 考慮空間3個(gè)相干信號(hào)，入射角度分別為θ1=0°、θ2=49°和θ3=61°，信噪比為3 dB，快拍數(shù)L=100，由于MUSIC算法和MVDR算法對(duì)相干信號(hào)失效，故只比較JLZA-DOA算法和SVD-MGPSR算法對(duì)相干信號(hào)的測(cè)向能力，圖4給出了相干信號(hào)時(shí)兩種算法的歸一化空間譜圖，從圖4可以看出，JLZA-DOA算法對(duì)于相干信號(hào)不能正確測(cè)向，而SVD-MGPSR算法能夠較強(qiáng)地分辨信號(hào)且信號(hào)相干性對(duì)測(cè)向影響不大，表現(xiàn)出對(duì)相干信號(hào)具有較強(qiáng)的DOA估計(jì)能力。

圖2 均方根誤差隨信噪比變化曲線

圖3 均方根誤差隨快拍數(shù)變化曲線

圖4 相干信號(hào)空間譜

5 結(jié)論

筆者在壓縮感知框架下建立了聯(lián)合稀疏信號(hào)模型，采用奇異值分解進(jìn)行信號(hào)降維處理能明顯降低計(jì)算量；與已有算法相比，新算法不僅能在低信噪比、小快拍數(shù)條件下完成高精度估計(jì)，而且也可以對(duì)任意相干信號(hào)進(jìn)行有效估計(jì)，具有較高的測(cè)向精度和角度分辨率。

[1] DAI J S, XU W C, ZHAO D A. Real-valued DOA estimation for uniform linear array with unknown mutual coupling[J]. IEEE Trans on Signal Processing, 2012,92(9):2056-2065.

[2] BURG J P. Maximum entropy special analysis [C]∥Proc of the 37th Meeting of the Society of Exploration Geophsicists.[S.l.]:[s.n.],1967:1213-1213.

[3] CAPON J. High-resolution frequency-wave number spectrum analysis [J].Proc IEEE, 1969,57(8):1408-1418.

[4] SCHIMIDT R O. A signal subspace approach to multiple emitter location and spectrum estimation [D]. Stanford: Stanford University, 1981.

[5] ROY R，KAILATH T．ESPRIT-estimation of signal parameters via rotational invariance techniques [J]. IEEE Trans on Acoustics Speech and Signal, 1986,37(7):984-995.

[6] DONOHO D L. Compressed sensing [J]. IEEE Trans on Information Theory, 2006,52(4):1289-1306.

[7] YANG Z, XIE L H, ZHANG C S. Off-grid direction of arrival estimation using sparse bayesian inference [J]. IEEE Trans on Signal Processing, 2013,61(1):38-43.

[8] 賀亞鵬,李洪濤,王克讓,等.基于壓縮感知的高分辨DOA估計(jì)[J].宇航學(xué)報(bào),2011,32(6):1344-1349.

[9] MD M H, KAUSHIK M. Direction-of-arrival estimation using a mixedl2,0norm approximation [J]. IEEE Trans on Signal Processing, 2010,58(9):4646-4655.

[10] FIGUEIREDO A T, ROBERT R D, WRIGHT S J. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems [J]. IEEE Journal of Selected Topics in Signal Processing, 2007,1(4):586-597.

[11] STROHMER T. Measure what should be measured: progress and challenge in compressive sensing[J]. IEEE Trans on Signal Processing Letters, 2012,19(12):887-893.

[12] 林波.基于壓縮感知的輻射源DOA估計(jì)[D].長(zhǎng)沙：國(guó)防科技大學(xué),2010.

[13] 石光明,劉丹華,高大化.壓縮感知回顧與展望[J].電子學(xué)報(bào),2009,37(5):1070-1081.

[14] SHISHKIN S L. Fast and robust compressive sensing method using mixed hadamard sensing matrix [J]. IEEE Journal on Emerging and Selected Topics in Circuits and Systems, 2012,2(3):353-361.

[15] BARZILAI J, BORWEIN J. Two point step size gradient methods [J]. IMA Journal of Numerical Analysis, 1988,8(1):141-148.

BAO Xiaolei:Lect.; School of Communication and Information, Shanghai Technical Institute of Electronics & Information, Shanghai 201411, China.

[編輯：王志全]

Sparse Reconstruction DOA Estimation Algorithm Based on Compressive Sensing Theory

BAOXiaolei,QUXinggen,WANGZhuoying

Based on the sparse property of the spatial targets distribution and the idea of compressive sensing (CS) theory, a multi-measurement gradient projection for sparse reconstruction algorithm based on singular value decomposition (SVD) was proposed, in which a multi-targets DOA estimation problem can be translated into a sparse signal reconstruction problem. Firstly, the signal was joint sparse representation by establishing an over-complete atom dictionary according to array manifold matrix. Then, SVD of information matrix of compressive sampling was done to reduce the amount of computation greatly. Finally, the sparse signal was reconstructed based on multi-measurement vectors gradient projection for sparse signal reconstruction algorithm so as to achieve DOA estimation. Compared with existing algorithm, the proposed algorithm not only has a smaller mean square error in the low SNR but also is able to correctly estimate the coherent signal. Moreover, it offers higher direction finding precision and angular resolution. The simulation results verify its effectiveness.

compressive sensing; DOA estimation; SVD; gradient projection

2015-08-25.

包曉蕾(1981-)，女，山東煙臺(tái)人，上海電子信息職業(yè)技術(shù)學(xué)院通信與信息工程系講師.

上海市教育委員會(huì)教育發(fā)展基金資助項(xiàng)目(09CGB06).

2095-3852(2015)06-0827-05

A

TN911

10.3963/j.issn.2095-3852.2015.06.036