閆鋒剛 劉 帥 金 銘 喬曉林

(哈爾濱工業(yè)大學(xué) 哈爾濱 150001)

1 引言

波達(dá)方向(Direction Of Arrival, DOA)估計(jì)是陣列信號(hào)處理的重要研究?jī)?nèi)容之一，在雷達(dá)、聲納、無線通信、無源定位等領(lǐng)域有著重要應(yīng)用[1-3]。自發(fā)表MUSIC(MUltiple SIgnal Classification)算法[4]以來，DOA估計(jì)進(jìn)入了一個(gè)新的超分辨階段。隨著應(yīng)用的深入，對(duì)2維DOA進(jìn)行估計(jì)的需求日益迫切，各種算法也不斷被提出[5-15]。然而，經(jīng)典MUSIC算法需在參數(shù)空間進(jìn)行極值搜索，計(jì)算量巨大；ESPRIT算法雖然避免了譜峰搜索，但需參數(shù)配對(duì)。因此，降低計(jì)算量的研究成為學(xué)者研究的一個(gè)熱點(diǎn)。文獻(xiàn)[10]提出的求根MUSIC算法用多項(xiàng)式求根代替了譜峰搜索，但僅適用于均勻線陣。文獻(xiàn)[11]對(duì)求根MUSIC法進(jìn)行了改進(jìn)，擴(kuò)展了其應(yīng)用范圍，但算法比較復(fù)雜且精度下降嚴(yán)重。文獻(xiàn)[12,13]將 ESPRIT算法應(yīng)用于2維DOA估計(jì)，但需要較高信噪比和快拍數(shù)，且算法計(jì)算量依然較大。文獻(xiàn)[14]提出的DOA矩陣法以及文獻(xiàn)[15]推廣的DOA時(shí)空矩陣法充分挖掘了相關(guān)矩陣特征向量包含的信息，降低了計(jì)算量，但算法性能有限且對(duì)陣列結(jié)構(gòu)要求十分嚴(yán)格，不利于工程應(yīng)用。

對(duì)于2維DOA估計(jì)，在高精度搜索步長(zhǎng)下，譜峰搜索占總計(jì)算量的比例大。因此，減少譜峰搜索次數(shù)是降低運(yùn)算量的關(guān)鍵。本文將 MUSIC算法原理推廣到共軛子空間，通過對(duì)原噪聲子空間及其共軛的交集進(jìn)行奇異值分解(Singular Value Decomposing, SVD)，實(shí)現(xiàn)了噪聲子空間的降維。接著，利用降維噪聲子空間與導(dǎo)向矢量及其共軛的雙正交性提出了一種新的2維陣列DOA快速估計(jì)算法。該算法能實(shí)現(xiàn)空間譜范圍的2倍壓縮，從而能將DOA估計(jì)的運(yùn)算量降低到傳統(tǒng)方法的50%。

2 MUSIC算法原理

2.1 陣列和數(shù)據(jù)模型

設(shè)M個(gè)坐標(biāo)為(xm,ym,0),m= 1 ,2,… ,M的陣元位于XOY平面，空間有L個(gè)輻射源Sl(t)，定義波達(dá)方向DOA為(θl,φl) ,l= 1 ,2,… ,L。其中，θl為信號(hào)入射方向與Z軸的夾角，φl為信號(hào)入射方向在XOY平面的投影與X軸的夾角，如圖1所示。對(duì)于遠(yuǎn)場(chǎng)窄帶信號(hào)，陣列一次快拍的接收數(shù)據(jù)為

圖1 2維陣列模型

2.2 MUSIC算法

設(shè)陣列接收數(shù)據(jù)協(xié)方差矩陣R，則由其定義可得

對(duì)R進(jìn)行特征值分解，有

式中Σ= diag(λ1,… ,λM)為對(duì)角矩陣，R的特征值為λ1≥ … ≥λM-L+1= … =λM=σ2，其對(duì)應(yīng)的特征向量為ei,i= 1,2,… ,M。W為特征向量矩陣。由S=[e1,e2,… ,eL]張成的子空間為信號(hào)子空間span(S)，而G= [eM-L+1,… ,eM]張成噪聲子空間span(G)且 滿足span(S) ⊥ span(G)及span(A)=span(S)。于是，可得

由此，可構(gòu)造MUSIC空間譜為

根據(jù)子空間正交性原理，在2維空間搜索，可得DOA為

3 基于降維噪聲子空間的DOA估計(jì)算法

3.1 譜函數(shù)構(gòu)造

對(duì)式(6)兩邊同取共軛，得

由此，定義譜函數(shù)為

圖2 虛擬輻射源與真實(shí)輻射源關(guān)系

可見， ?θ∈ [ 0,π/2],φ∈ [ 0,π]，有

將G寫為列向量G= [g1,g2,… ,gv],v=M-L，則

將式(14)帶入式(11)，得

由式(6)及式(15)，可得

上述分析表明：f(θ,φ)關(guān)于φ=π對(duì)稱，并在輻射源及其鏡像位置同時(shí)產(chǎn)生極值，因而實(shí)現(xiàn)了MUSIC譜的壓縮，因此我們將其稱為MUSIC對(duì)稱壓縮譜(MUSIC Symmetrical Compressed Spectrum, MSCS)。傳統(tǒng)MUSIC譜覆蓋整個(gè)2維空間，而MSCS將DOA估計(jì)的譜范圍壓縮至原來的一半，因而總體運(yùn)算速度也將提高約1倍。

3.2 噪聲子空間降維及譜函數(shù)化簡(jiǎn)

本文在構(gòu)造共軛噪聲空間*G的同時(shí)，等效于在原輻射源E的鏡像位置增加了一個(gè)虛擬輻射源E'。若輻射源個(gè)數(shù)為L(zhǎng)，則MSCS等效地在空間新增加了L個(gè)虛擬輻射源。因而，信號(hào)子空間被升高了L維；相應(yīng)地，噪聲子空間被降低了L維。

設(shè)升維后的信號(hào)子空間為span(?)，降維后的噪聲子空間為span(?)。則span(?)應(yīng)為原噪聲子空間span(G)和新增噪聲子空間 s pan(G*)的交集，而span(?)則由原信號(hào)子空間span(S)與新增信號(hào)子空間的和構(gòu)成。由GGH+SSH=I，得G*GT+S*ST=I。因此，新增信號(hào)子空間即為 s pan(S*)。所以，

為了求解span(?)，我們給出如下定理：

定理令Ψ=I-GGHG*GT，設(shè)Ψ零空間為υ(Ψ)，那么有

證明設(shè)向量γ∈span(?)，則SSHγ=OM×1，從而，

同理γ=G*GTγ，故γ=GGHG*GTγ。從而，

所以γ∈υ(Ψ) ，這表明：

反之，若γ∈υ(Ψ) ，則

由于G*GT為 s pan(G*)的投影矩陣，故

所以，有

由式(21)和式(25)知定理成立。 證畢

上述定理表明，降維噪聲子空間span(?)與矩陣Ψ的零空間相同。前文中，我們實(shí)質(zhì)上是基于矩陣(I-Ψ)的非簡(jiǎn)化列空間(列向量線性相關(guān))構(gòu)造了MSCS，因而存在計(jì)算上的冗余。

對(duì)Ψ進(jìn)行SVD，可得

其中對(duì)角陣Λ=diag(σ1,σ2,… ,σM)。由于：r an k(GGH)=rank(G*GT) =M-L，故rank(Ψ)=2L。因此，Λ的對(duì)角元素滿足：

于是，我們就得到了MSCS的簡(jiǎn)潔形式為

綜上所述，本文提出的基于降維噪聲子空間的DOA快速估計(jì)算法步驟如下：

步驟1 計(jì)算矩陣Ψ并對(duì)其進(jìn)行奇異值分解得到降維噪聲子空間的標(biāo)準(zhǔn)正交基ΦM-2L；

步驟2 利用ΦM-2L構(gòu)造MSCS，并在其半譜內(nèi)搜索，給出DOA信息(θl,φl) 或其鏡像(θl,φl±π) ,l=1,2,…,L；

步驟3 在MUSIC譜(θl,φl±π) ,l= 1 ,2,…,L的鄰近小區(qū)域進(jìn)行極值檢驗(yàn)，存在極值的位置即為真實(shí)DOA。

從上述步驟可見，如果真實(shí)輻射源位置本身對(duì)稱，由于本文算法在“步驟3”中對(duì)MUSIC譜的對(duì)稱位置鄰域進(jìn)行了極值檢驗(yàn)，因此不會(huì)丟失真實(shí)DOA信息。

3.3 算法性能分析

設(shè)特征值分解得信號(hào)、噪聲子空間估計(jì)值分別為和，設(shè)本文所得到的信號(hào)、噪聲子空間為new和new。為了評(píng)價(jià)子空間的估計(jì)性能，引入信號(hào)、噪聲子空間估計(jì)誤差函數(shù)C()和C()，其定義分別為

圖3和圖4分別給出了M=16時(shí)，不同L下MUSIC和MSCS子空間準(zhǔn)確度對(duì)比關(guān)系。由圖可見：當(dāng)L較小時(shí)，MSCS噪聲子空間和信號(hào)子空間準(zhǔn)確度均比MUSIC略差。隨著L增大，MSCS噪聲子空間與 MUSIC噪聲子空間差異變小；而在L的整個(gè)變化過程中，MSCS信號(hào)子空間準(zhǔn)確度都較MUSIC略差，且隨L增加而更甚。這是容易理解的，因?yàn)镸SCS比MUSIC多了L維“鏡像”信號(hào)，從而累積了子空間的估計(jì)誤差，這使得其總體誤差較大。

但是，這里需要強(qiáng)調(diào)的是：本文算法由于在得到“鏡像”信號(hào)后，又在最后一步對(duì) MUSIC譜對(duì)稱位置進(jìn)行了峰值檢驗(yàn)，故其估計(jì)精度與 MUSIC算法一致。因此，本文算法在提高 DOA估計(jì)速度的同時(shí)，保持了估計(jì)精度未下降。

采用式(7)和式(30)計(jì)算一個(gè)譜值點(diǎn)分別需要3M2-2ML和 3M2-4ML次復(fù)數(shù)乘法。若DOA搜索步長(zhǎng)為δ，則 MUSIC算法的計(jì)算量為π(3M2- 2ML)/δ。本文構(gòu)造 MSCS時(shí)，矩陣Ψ奇異值分解運(yùn)算量[16]為M(M+ 1 )2+ 1 7(M+ 1 )3/3 ≈ 6M3，故本文算法總運(yùn)算量為 6M3+π(3M2- 4ML)/2δ。通常δ?π，故π(3M2- 4ML) /δ? 6M3。因而，相比于MUSIC算法，本文算法將DOA估計(jì)速度提高為原來的2倍左右。

4 仿真及分析

實(shí)驗(yàn)設(shè)置陣元數(shù)M=16，陣元間距為半個(gè)波長(zhǎng)，快拍數(shù)N= 2 00, Monte-Carlo實(shí)驗(yàn)次數(shù)均為200次。

為了對(duì)比MSCS與MUSIC譜的差異，選取二者空間譜的公共部分進(jìn)行DOA估計(jì)實(shí)驗(yàn)。實(shí)驗(yàn)中，L增加時(shí)，增加的輻射源均非相干。

圖7和圖8給出了以L(L＜M/ 2)為參變量，采用MUSIC和MSCS進(jìn)行DOA估計(jì)的均方根誤差(Root Mean Square Error, RMSE)和成功概率與SNR的關(guān)系。其中，RMSE定義為

圖3 噪聲子空間準(zhǔn)確度

圖4 信號(hào)子空間準(zhǔn)確度

圖5 MUSIC 空間譜

圖6 MSCS空間譜

圖7 DOA估計(jì)均方根誤差

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

由圖可見：當(dāng)L=2時(shí)，MSCS與MUSIC的估計(jì)誤差相當(dāng)，隨著L增大，MSCS的估計(jì)誤差較MUSIC略差。這與性能分析部分MSCS子空間估計(jì)精度的現(xiàn)象也保持一致。然而，在整個(gè)L的變化過程中，MSCS的 DOA估計(jì)成功概率均略優(yōu)于MUSIC，這是因?yàn)镸SCS噪聲子空間的維度下降而使得其空間譜變得更為“尖銳”了的緣故。

表1給出了不同SNR下，采用MUSIC算法和本文算法進(jìn)行 DOA估計(jì)所需時(shí)間的對(duì)比關(guān)系。實(shí)驗(yàn)中，采用Matlab7.0自帶的“cputime”命令記錄程序運(yùn)行時(shí)間，DOA估計(jì)的范圍為40°≤θ≤60°,10°≤φ≤220°。由表 1可見：MUSIC算法 DOA估計(jì)的時(shí)間約為 0.1205 s；而本文算法則需要約0.0601 s。這說明本文算法能將DOA估計(jì)的速度提高約2倍。

5 結(jié)論

波達(dá)方向估計(jì)是空間譜估計(jì)的重要研究?jī)?nèi)容。傳統(tǒng)MUSIC算法需在2維空間進(jìn)行峰值搜索，計(jì)算量巨大。本文從等效添加虛擬輻射源的角度入手，提出了一種基于降維噪聲子空間的 2維陣列快速DOA估計(jì)算法，在保持DOA估計(jì)精度不下降的同時(shí)將計(jì)算量降低了約50%，并經(jīng)過適當(dāng)變換能適用于任意陣型，具有較大的理論和應(yīng)用價(jià)值。

表1 DOA估計(jì)用時(shí)比較 (s)

[1]Wan F, Zhu W P, and Swamy M N S. Spacial extrapolation based blind DOA estimation approach for closely spaced sources[J].IEEE Transactions on Aerospace and Electronic Systems, 2010, 46(2): 569-582.

[2]Chen F J, Kwong S, and Kok C W. Esprit like two dimensional DOA estimation for coherent signals[J].IEEE Transactions on Aerospace and Electronic Systems, 2010,46(3): 1477-1484.

[3]Yo Han Ko,et al.. 2-D DOA estimation with cell searching for a mobile relay station with uniform circular array[J].IEEE Transactions on Communications, 2010, 58(10):2805-2809.

[4]Schmidt R O. Multiple emitter location and signal parameter estimation[J].IEEE Transactions on Antennas and Propagation, 1986, 34(3): 276-280.

[5]Mathews C and Zohtowski M D. Eigenstructure techniques for 2-D angle estimation with uniform circular array[J].IEEE Transactions on Signal Processing, 1994, 42(9): 2395-2404.

[6]Wang Guang-min, Xin Jin-ming, Zheng Nan-ning,et al..Computationally efficient subspace-based method for two-dimensional direction estimation with L-shaped array[J].IEEE Transactions on Signal Processing,2011, 59(7):3197-3212.

[7]Tayem N and Kwon H. L-shape 2-dimensional arrival angle estimation with propagator method[J].IEEE Transactions on Antennas Propagation, 2005, 53(5): 1622-1630.

[8]顧陳, 何勁, 李彧晟, 等. 基于擴(kuò)展孔徑波達(dá)方向矩陣法的高精度二維DOA估計(jì)[J]. 電子與信息學(xué)報(bào), 2010, 32(12):2879-2883.

Gu Chen, He Jin, Li Yu-sheng,et al.. Extended aperture DOA matrix method for high accurate 2D direction estimation[J].Journal of Electronics&Information Technology, 2010, 32(12): 2879-2883.

[9]艾名舜, 馬紅光, 劉剛. 基于噪聲子空間解析形式的快速DOA估計(jì)算法[J]. 電子與信息學(xué)報(bào), 2010, 32(5): 1071-1075.

Ai Ming-shun, Ma Hong-guang, and Liu Gang. DOA estimation algorithm base on analytical solution of noise subspace[J].Journal of Electronics&Information Technology, 2010, 32(5): 1071-1075.

[10]Barabell A J. Improving the resolution performance of eigenstructure based direction finding algorithms [C]. IEEE International Conference on Acoustics, Speech, and Signal Processing, Boston, 1983: 336-339.

[11]Rubsamen M and Gershman A B. Direction-of-arrival estimation for nonunifom sensor arrays: from manifold separation to Fourier domain music methods[J].IEEE Transactions on Signal Processing, 2009, 57(2): 588-599.

[12]Roy R and Kailath T. ESPRIT-Estimation of signal parameters via rotational invariance techniques[J].IEEE Transactions on Acoustics,Speech and Signal Processing,1989, 37(7): 984-995.

[13]Jian C, Wang S, and Lin L. Two-dimensional DOA estimation of coherent signals based on 2D unitary ESPRIT method [C]. The 8th International Conference on Signal Processing, Beijing, 2006: 16-20.

[14]季飛, 余華, 謝澤明, 等. 基于DOA矩陣法的矢量傳感器陣列二維波達(dá)方向估計(jì)[J]. 電子與信息學(xué)報(bào), 2008, 30(8):1886-1889.

Ji Fei, Yu Hua, Xie Ze-ming,et al.. A DOA matrix algorithm for 2-D direction finding estimation with electromagnetic vector sensor arrays[J].Journal of Electronics&Information Technology, 2008, 30(8): 1886-1889.

[15]Du Rui-yan, Wang Jin-kuan, and Liu Fu-lai. Space-time matrix method for joint frequency and 2-D DOA estimation[C]. The 1st International Conference on Information Science and Engineering, Qinhuangdao, China, 2009: 671-674.

[16]張賢達(dá). 矩陣分析與應(yīng)用[M]. 北京: 清華大學(xué)出版社, 2004, 9:611-613.

[17]張立杰, 黃建國(guó), 史文濤, 等. 正交投影陣列信號(hào)子空間估計(jì)方法[J]. 系統(tǒng)工程與電子技術(shù), 2009, 31(9): 2063-2066.

Zhang Li-jie, Huang Jian-guo, Shi Wen-tao,et al..Orthogonal projection method for array signal subspace estimation[J].Systems Engineering and Electronic, 2009,31(9): 2063-2066.

[18]Huang Lei, Wu Shunjun, Feng Dazhang,et al.. Low complexity method for signal subspace fitting [J].Electronics Letters, 2004, 40(14): 847-848.