基于降維噪聲子空間的二維陣列DOA估計算法

閆鋒剛 劉 帥 金 銘 喬曉林

(哈爾濱工業大學 哈爾濱 150001)

1 引言

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

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

2 MUSIC算法原理

2.1 陣列和數據模型

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

圖1 2維陣列模型

2.2 MUSIC算法

設陣列接收數據協方差矩陣R，則由其定義可得

對R進行特征值分解，有

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

由此，可構造MUSIC空間譜為

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

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

3.1 譜函數構造

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

由此，定義譜函數為

圖2 虛擬輻射源與真實輻射源關系

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

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

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

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

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

3.2 噪聲子空間降維及譜函數化簡

本文在構造共軛噪聲空間*G的同時，等效于在原輻射源E的鏡像位置增加了一個虛擬輻射源E'。若輻射源個數為L，則MSCS等效地在空間新增加了L個虛擬輻射源。因而，信號子空間被升高了L維；相應地，噪聲子空間被降低了L維。

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

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

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

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

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

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

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

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

所以，有

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

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

對Ψ進行SVD，可得

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

于是，我們就得到了MSCS的簡潔形式為

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

步驟1 計算矩陣Ψ并對其進行奇異值分解得到降維噪聲子空間的標準正交基ΦM-2L；

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

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

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

3.3 算法性能分析

設特征值分解得信號、噪聲子空間估計值分別為和，設本文所得到的信號、噪聲子空間為new和new。為了評價子空間的估計性能，引入信號、噪聲子空間估計誤差函數C()和C()，其定義分別為

圖3和圖4分別給出了M=16時，不同L下MUSIC和MSCS子空間準確度對比關系。由圖可見：當L較小時，MSCS噪聲子空間和信號子空間準確度均比MUSIC略差。隨著L增大，MSCS噪聲子空間與 MUSIC噪聲子空間差異變小；而在L的整個變化過程中，MSCS信號子空間準確度都較MUSIC略差，且隨L增加而更甚。這是容易理解的，因為MSCS比MUSIC多了L維“鏡像”信號，從而累積了子空間的估計誤差，這使得其總體誤差較大。

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

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

4 仿真及分析

實驗設置陣元數M=16，陣元間距為半個波長，快拍數N= 2 00, Monte-Carlo實驗次數均為200次。

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

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

圖3 噪聲子空間準確度

圖4 信號子空間準確度

圖5 MUSIC 空間譜

圖6 MSCS空間譜

圖7 DOA估計均方根誤差

圖8 DOA估計成功概率

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

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

5 結論

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

表1 DOA估計用時比較 (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]顧陳, 何勁, 李彧晟, 等. 基于擴展孔徑波達方向矩陣法的高精度二維DOA估計[J]. 電子與信息學報, 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]. 電子與信息學報, 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矩陣法的矢量傳感器陣列二維波達方向估計[J]. 電子與信息學報, 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]張賢達. 矩陣分析與應用[M]. 北京: 清華大學出版社, 2004, 9:611-613.

[17]張立杰, 黃建國, 史文濤, 等. 正交投影陣列信號子空間估計方法[J]. 系統工程與電子技術, 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.