朱進(jìn)勇，王立冬，孟亞峰

(軍械工程學(xué)院 電子與光學(xué)工程系，石家莊 050003)

虛擬陣列Khatri-Rao積與子空間聯(lián)合稀疏表示的DOA估計方法

朱進(jìn)勇，王立冬，孟亞峰

(軍械工程學(xué)院 電子與光學(xué)工程系，石家莊 050003)

利用目標(biāo)信號在空域分布的稀疏性，該文提出了一種基于虛擬陣列Khatri-Rao(KR)積與信號子空間聯(lián)合稀疏表示的單快拍DOA估計方法；該方法利用單次快拍的采樣數(shù)據(jù)，構(gòu)造出雙向虛擬陣列數(shù)據(jù)，并對虛擬陣列數(shù)據(jù)的協(xié)方差矩陣進(jìn)行KR積變換處理，然后對向量化后的數(shù)據(jù)進(jìn)行順序重構(gòu)，利用重構(gòu)矩陣的大奇異值對應(yīng)的左奇異向量為估計信號子空間；最后，利用凸優(yōu)化工具箱對稀疏模型進(jìn)行二階凸規(guī)劃的優(yōu)化求解，得到高精度的DOA估計值；仿真實驗驗證了算法的有效性，在低信噪比下比傳統(tǒng)MUSIC和OMP算法具有更高的估計精度。

虛擬陣列；Khatri-Rao積；稀疏表示；單快拍；波達(dá)方向估計

0 引言

波達(dá)方向角(direction of arrival, DOA)估計[1-2]作為陣列信號處理領(lǐng)域研究的熱點內(nèi)容，被廣泛應(yīng)用于無線通信、測速、雷達(dá)、地質(zhì)勘探等眾多領(lǐng)域。但傳統(tǒng)MUSIC算法和ESPRIT算法為代表的子空間類算法對空域存在的相關(guān)信源進(jìn)行DOA估計時便失效。近年來，壓縮感知理論的興起，使得DOA估計有了質(zhì)的飛躍性發(fā)展，與傳統(tǒng)的子空間類算法不同，稀疏重構(gòu)類DOA估計算法利用目標(biāo)信號在空域的稀疏特性，可以高概率重構(gòu)原始信號，從而不考慮子空間維數(shù)和秩的虧損問題[3-5]。

基于稀疏分解的DOA估計算法不需要關(guān)注目標(biāo)信號的統(tǒng)計特性、是否相關(guān)以及在較小的快拍數(shù)下就能得到較好的估計精度和分辨精度，因此，國內(nèi)外學(xué)者提出了很多基于稀疏分解的DOA 估計方法。文獻(xiàn)[6]提出了L1SVD算法，采用奇異值分解降低算法復(fù)雜度和噪聲的敏感性，結(jié)合凸優(yōu)化進(jìn)行稀疏重構(gòu)目標(biāo)，解決了壓縮感知在理論在DOA 估計上的基本問題，但是該方法運算量較大；文獻(xiàn)[7]提出了Khatri-Rao子空間分解方法，這種方法能夠處理N≤2M-2個目標(biāo)信源個數(shù)的情況，但該方法是建立在信源的數(shù)目N為已知條件的情況下才能準(zhǔn)確的分離出信號子空間和噪聲子空間；文獻(xiàn)[8]提出了基于KR積變換和空域平滑理論的DOA 估計算法，利用空域平滑處理，使協(xié)方差恢復(fù)滿秩，減少了運算量；文獻(xiàn)[9]提出了基于信號子空間的迭代加權(quán)最小方差算法，利用奇異值分解得到低維的接收信號矩陣，采用迭代加權(quán)最小方差求解欠定的稀疏表示問題，能對任意相關(guān)信號進(jìn)行處理，但是該方法需要已知信源個數(shù)；文獻(xiàn)[10]提出了對多快拍接收信號的協(xié)方差矩陣進(jìn)行l(wèi)1 2混合范數(shù)聯(lián)合稀疏表示，但由于需要對協(xié)方差矩陣求逆，因此在快拍數(shù)較少的情況下，導(dǎo)致算法的性能急劇下降。

本文提出一種建立在KR積理論變換的基礎(chǔ)上，將虛擬陣列技術(shù)和信號子空間相結(jié)合進(jìn)行聯(lián)合稀疏表示的DOA 估計方法。該方法利用陣列接收信號的單次快拍數(shù)據(jù)，進(jìn)行雙向虛擬平移處理，以達(dá)到處理相干信號的目的；在此基礎(chǔ)上，利用KR積變換理論對虛擬平移的數(shù)據(jù)矩陣進(jìn)行處理，得到向量化空間矩陣；然后利用奇異值分解向量化空間矩陣提取信號子空間，得到低維的信號子空間，再采用二階錐規(guī)劃求解的方法以及結(jié)合l1范數(shù)凸優(yōu)化[11]問題進(jìn)行DOA求解。該方法從信號的本質(zhì)上解決了相干問題，減少了運算量。通過仿真實驗證明了本文算法在低信噪比情況下，比傳統(tǒng)的MUSIC算法和壓縮感知DOA 估計算法具有更高的估計精度。

1 信號模型

假設(shè)M個各向同性陣元組成的等距均勻線陣(Uniform linear array, ULA)，陣元間距為d=λ/2，λ為信號波長?？沼虼嬖贜個遠(yuǎn)場窄帶信號s1(t)，s2(t),…,sN(t),分別以方位角θ1,θ2,…，θN入射到ULA。陣列的接收信號可寫為:

X(t)=A(θ)S(t)+N(t)

(1)

式中，X(t)=[x1(t),x2(t),…,xM(t)]T為陣列天線接收矢量，A(θ)=[α(θ1),α(θ2),…,α(θN)]為M×N維的陣列流行矩陣，其中，α(θi)=[1,e-j2πdsinθ1/λ,…,e-j2π(i-1)dsinθi/λ]T,1≤i≤N為目標(biāo)信源入射角θi的導(dǎo)向矢量，N(t)為M×1維的陣列接收噪聲矩陣。

2 本文算法

2.1 雙向虛擬陣列

根據(jù)虛擬陣列技術(shù)，每次將M個陣元向右平移d，形成一個虛擬子陣，如此反復(fù)M-1次，與真實的天線陣列接收數(shù)據(jù)形成一個方陣接收數(shù)據(jù)。對于第p次右相平移得到的陣列輸出信號表示為:

Xp(t)=ADrp-1S(t)+Np(t)

(2)

式中，Drp-1為N×N的對角陣的p-1次乘方，平移矩陣Dr表示為:

Dr=diag(e-j2πdsinθ1/λ,…,e-j2πdsinθN/λ)

(3)

利用M-1次的陣列虛擬數(shù)據(jù)構(gòu)造如下接收陣列信號矩陣:

[X(t),X1(t),…,X(M-1)(t)]T

(4)

那么:

Xr=X(t)[1,Dr,…,Dr(M-1)]T

(5)

因此，獲得陣列虛擬處理后的協(xié)方差矩陣為:

(6)

同理，左向虛擬陣列的協(xié)方差矩陣為:

Rl=E[XlXlH]=

(7)

式中，C2=[1,Dl,…，Dl(M-1)]T，Dl=diag(ej2πdsinθ1/λ,…,ej2πdsinθN/λ)。

2.2 協(xié)方差矩陣的KR積處理

假設(shè)矩陣A=[α1,α2,…,αF]∈CI×F和B=[b1,b2,…,bF]∈CJ×F，定義兩個矩陣的Khatri-Rao積[12]為:

A×B=[a1?b1,a2?b2,…，aF?bF]∈CIJ×F

(8)

式中，“×”和“?”分別為KR積與Kronecker積的運算符號。

假設(shè)ai=[ai1,ai2,…,aiI]T∈CI，bi=[bi1,bi2,…，biI]T∈CI(i=1,2,…,F)，因此Kronecker積表示為:

(9)

因此，根據(jù)KR積變換理論，對式(6)和式(7)進(jìn)行向量化運算可以得到:

(10)

(11)

式中，“vec”表示矩陣的向量化處理，Γr=AC1和Γl=AC2為旋轉(zhuǎn)陣列流行矩陣。

利用左右虛擬陣列協(xié)方差的向量化形式，構(gòu)造出左右虛擬陣列的稀疏重構(gòu)觀測向量：

(12)

(13)

(14)

YSV=AXSV+NSV

(15)

通過SVD分解將M×M維的數(shù)據(jù)變成了M×N維，不僅降低了數(shù)據(jù)的維數(shù)還保持了觀測矩陣的稀疏性和稀疏結(jié)構(gòu)，并不影響DOA估計。

根據(jù)壓縮感知理論中的信號恢復(fù)理論，可以容易的得到目標(biāo)函數(shù)為:

(16)

1)獲取陣列觀測的單次快拍數(shù)據(jù)X(t)，分別進(jìn)行左右平移構(gòu)造出虛擬子陣Xr和Xl；

2)對虛擬子陣Xr和Xl進(jìn)行協(xié)方差處理，得到協(xié)方差矩陣Rr和Rl；

3)利用KR積理論對Rr和Rl進(jìn)行向量化運算，并求得平均稀疏觀測向量y；

4)對觀測向量進(jìn)行順序重構(gòu)得到M×M的方陣Y；

6)利用凸優(yōu)化工具箱求解優(yōu)化函數(shù)。

3 算法性能分析

從以上的推導(dǎo)過程可以看出，本文算法在l1-SVE算法的基礎(chǔ)上運用虛擬陣列技術(shù)和KR積變換理論對信號進(jìn)行聯(lián)合稀疏表示。該算法的的計算量主要集中在凸優(yōu)化工具箱進(jìn)行二階錐規(guī)劃問題的求解，因此與l1-SVD算法的復(fù)雜度O(N3P3)相差無幾，P為整個空域平均分的份數(shù)。

本文算法通過雙向虛擬陣列技術(shù)將單快拍的接收數(shù)據(jù)構(gòu)成方陣的協(xié)方差矩陣恢復(fù)到滿秩，使該文算法具有了相關(guān)信號的處理能力。需要注意的是，本文算法只采用了一次的快拍數(shù)據(jù)，在進(jìn)行SVD分解估計信號子空間的過程中得不到多次快拍中信號的能量累積，所以并不能有效的估計DOA。

4 仿真實驗與分析

4.1 算法的估計特性比較

實驗中采用8個陣元的ULA，陣元間距為半波長，兩個遠(yuǎn)場窄帶不相關(guān)的信源分別以θ1=20°和θ2=40°的方位角入射到ULA。本文算法采用的快拍數(shù)K1=1，l1-SVD，OMP和MUSIC算法采用的快拍數(shù)K2=50。將空域的-90°～90°范圍以1°的等間隔劃分為181份來構(gòu)造過完備陣列流行矩陣，圖1為空間譜圖。

圖1 空間譜

4個遠(yuǎn)場窄帶信源分別以θ3=-50°，θ4=-10°，θ5=20°和θ6=70°的方位角入射到ULA，其中θ3和θ4為相干信源，θ5和θ6為不相關(guān)信源，其空間譜圖如圖2所示。

圖2 相關(guān)信源空間譜

由圖1和圖2可知，本文算法對相關(guān)信號和不相關(guān)信號均能進(jìn)行有效的DOA估計，傳統(tǒng)的MUSIC算法未能對相關(guān)信號進(jìn)行DOA估計，l1-SVD和OMP能完成對相關(guān)信號的DOA估計。本文算法通過虛擬陣列技術(shù)將接收信號的協(xié)方差恢復(fù)到了滿秩，有了更強(qiáng)的解相干信號的能力，在后期的信號恢復(fù)上是由于基于壓縮感知理論的DOA估計算法是將空域的目標(biāo)信號的DOA估計當(dāng)作一個稀疏向量的重構(gòu)恢復(fù)問題，這與MUSIC算法的子空間分解類思想是完全不同的。

4.2 算法的統(tǒng)計特性比較

在不同的信噪比下，本文算法與，OMP和MUSIC算法對信號角度估計的均方根誤差值(Root Mean Square Error, RMSE)。實驗中采用兩個窄帶遠(yuǎn)場信號θ1=20°和θ2=40°，信噪比變化范圍為-10～30 dB，仿真結(jié)果為50次Monte-Carlo統(tǒng)計結(jié)果，其它條件與實驗一相同。這里定義角度估計的均方根誤差為:

(17)

由圖3可知，在低SNR時，本文算法，l1-SVD和OMP算法的RMSE遠(yuǎn)高于傳統(tǒng)的MUSIC算法，同時，本文算法的RMSE也小于l1-SVD和OMP算法的RMSE；隨著SNR的增大，四種算法的RMSE均減小。

由圖4可知，隨著ULA陣元數(shù)的增加，可以提高本文算法的DOA估計精度。

圖3 DOA估計的RMSE

圖4 陣元變化下的DOA估計的RMSE

4.3 成功概率比較

為研究本文算法對信號的估計成功概率，兩個窄帶遠(yuǎn)場信號θ1=20°和θ2=40°，仿真結(jié)果為50次Monte-Carlo統(tǒng)計結(jié)果，當(dāng)DOA估計結(jié)果在[19.8°,20.2°]和[39.8°,40.2°]的區(qū)間時，認(rèn)定一次成功，成功概率隨SNR的變化情況如圖5所示

可以看出在相同的SNR下，本文算法的DOA估計成功概率高于其它三種算法，在SNR≥9dB時，就能夠?qū)崿F(xiàn)100%的估計成功概率，而OMP算法具有最差的DOA估計成功概率，這從空間譜圖也能看出；這是由于OMP算法在進(jìn)行觀測矩陣選擇時具有比等距約束性(restrictedisometryproperty,RIP)條件[13]更加嚴(yán)格的要求，這就需要確定的迭代次數(shù)和更多的觀測值，即使是這樣也不能保證把任何信號都能夠精確的重構(gòu)出來；l1-SVD算法采用直接對待恢復(fù)信號進(jìn)行l(wèi)1范數(shù)約束，這將導(dǎo)致在恢復(fù)信號的過程中得不到更稀疏的結(jié)果，進(jìn)而在DOA估計的空間譜上出現(xiàn)偽譜峰，降低了DOA的估計精度。

圖5 DOA估計成功概率

5 結(jié)束語

本文在KR積變換理論的基礎(chǔ)上，提出了一種在單快拍條件下對順序重構(gòu)虛擬陣列的提取信號子空間稀疏表示的DOA估計算法，經(jīng)過仿真驗證本文算法對相關(guān)信號具有更強(qiáng)的處理能力，尤其是在低信噪比情況下相對于l1-SVD，OMP和MUSIC算法估計精度也大大提高。然而，本文算法利用凸優(yōu)化工具箱處理，在運算量方面并沒有太大的優(yōu)勢，這方面還需要更進(jìn)一步的研究。

[1] GRANT M and BOYD S. CVX: Matlab software for disciplined convex programming[OL].http://cvxr.com/cvx,2012.

[2] 王 凌，李國林，劉堅強(qiáng)，等.一種基于數(shù)據(jù)矩陣重構(gòu)的相干信源二維測向新方法[J].西安電子科技大學(xué)學(xué)報，2013，40(2):130-137.

[3] 杜新鵬.聯(lián)合稀疏恢復(fù)新型算法及其應(yīng)用研究[D].長沙：國防科技大學(xué),2013.

[4] 劉慶華，歐陽繕，何振清.準(zhǔn)平穩(wěn)信號的Khatri-Rao積聯(lián)合稀疏分解DOA估計方法[J].系統(tǒng)工程與電子技術(shù)，2012,34(9):1753-1757.

[5] 嚴(yán)金山，彭秀艷，王咸鵬.基于酉變換的虛擬陣列DOA估計算法[J].哈爾濱工業(yè)大學(xué)學(xué)報，2012,44(4):136-140.

[6] Malioutov D, Mujdat C, Willsky A. A sparse signal reconstruction perspective for source localization with sensor arrays[J].IEEE Transaction on Signal Processing,2005,53(8):3010-3022.

[7] Ma W K, Hsieh T H, Chi C Y. DOA estimation of quasi-stationary signals with less sensors than sources and unknown spatial noise covariance: A Khatri-Rao subspace approach[J].IEEE Transaction. On Signal Processing, 2010,58(4):2168-2180.

[8] 蔡晶晶，宗 汝，蔡 輝.基于空域平滑稀疏重構(gòu)的DOA估計算法[J].電子與信息學(xué)報,2016,38(1):168-173.

[9] 解 虎，馮大政，魏倩茹.采用信號子空間稀疏表示的DOA估計方法[J].系統(tǒng)工程與電子技術(shù)，2015,37(8):1717-1722.

[10] Yin J H, Chen T Q. Direction-of-arrival estimation using a sparse representation of array covariance vectors[J].IEEE Transaction. on Signal Processing letter, 2011,59(9):4489-4493.

[11] Tibshirani R. Regression shrinkage and selection via the LASSO[J].Journal of the Royal Statistical Society: Series B,1996,58(1):267-288.

[12] Ma W K, Hsieh T H, Chi C Y. DOA estimation of quasi-stationary signals with less sensors than sources and unknown spatial noise covariance: A Khatri-Rao subspace approach[J].IEEE Trans. On signal Processing, 2010,58(4):2168-2180.

[13] Candes E J, Tao T. Decoding by linear programming[J]. IEEE Transactions on Information Theory,2005,51(12):4203-4215.

Virtual Array Khatri-Rao Product and Subspace Joint Sparse Representation Method of DOA Estimation

Zhu Jinyong, Wang Lidong, Meng Yafeng

(Department of Electronic and Optical Engineering, Ordnance Engineering College, Shijiazhuang 050003,China)

Using the target signal in the spatial distribution of sparse, this paper puts forward a Khatri-Rao(KR) product based on virtual array and signal subspace joint sparse representation of single snapshot DOA estimation method. The method uses a single snapshot sampling data, constructs the two-way virtual array data, and the covariance matrix of the virtual array data for KR product transformation process, and then to reconstruct the order of data after vectorization, by using the large singular values of reconstruction matrix left singular vectors of the corresponding to estimate the signal subspace; Finally, using convex optimization toolbox for sparse matrix model of quadratic convex programming optimization solution, get high accuracy DOA estimate. Simulation experiments verify the effectiveness of the algorithm, under the low SNR has higher estimation accuracy than traditional MUSIC,SVD and OMP algorithm.

virtual array; Khatri-Rao product; sparse representation; single snapshot; DOA

2016-10-11；

2016-12-20。

國家自然科學(xué)基金資助項目(61372039)。

朱進(jìn)勇(1992-)，男，碩士研究生，主要從事信息與通信方向的研究。

1671-4598(2017)05-0147-03

10.16526/j.cnki.11-4762/tp.2017.05.041

TN911.7

A