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

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

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

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

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

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

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

0 引言

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

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

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

1 信號(hào)模型

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

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個(gè)陣元向右平移d，形成一個(gè)虛擬子陣，如此反復(fù)M-1次，與真實(shí)的天線陣列接收數(shù)據(jù)形成一個(gè)方陣接收數(shù)據(jù)。對(duì)于第p次右相平移得到的陣列輸出信號(hào)表示為:

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

(2)

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

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

(3)

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

[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，定義兩個(gè)矩陣的Khatri-Rao積[12]為:

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

(8)

式中，“×”和“?”分別為KR積與Kronecker積的運(yùn)算符號(hào)。

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

(9)

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

(10)

(11)

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

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

(12)

(13)

(14)

YSV=AXSV+NSV

(15)

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

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

(16)

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

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

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

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

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

3 算法性能分析

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

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

4 仿真實(shí)驗(yàn)與分析

4.1 算法的估計(jì)特性比較

實(shí)驗(yàn)中采用8個(gè)陣元的ULA，陣元間距為半波長，兩個(gè)遠(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個(gè)遠(yuǎn)場窄帶信源分別以θ3=-50°，θ4=-10°，θ5=20°和θ6=70°的方位角入射到ULA，其中θ3和θ4為相干信源，θ5和θ6為不相關(guān)信源，其空間譜圖如圖2所示。

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

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

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

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

(17)

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

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

圖3 DOA估計(jì)的RMSE

圖4 陣元變化下的DOA估計(jì)的RMSE

4.3 成功概率比較

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

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

圖5 DOA估計(jì)成功概率

5 結(jié)束語

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

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

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

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

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

[5] 嚴(yán)金山，彭秀艷，王咸鵬.基于酉變換的虛擬陣列DOA估計(jì)算法[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào)，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ì)算法[J].電子與信息學(xué)報(bào),2016,38(1):168-173.

[9] 解 虎，馮大政，魏倩茹.采用信號(hào)子空間稀疏表示的DOA估計(jì)方法[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é)基金資助項(xiàng)目(61372039)。

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

1671-4598(2017)05-0147-03

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

TN911.7

A