聯合互協方差矩陣的快速波達方向估計

閆鋒剛, 榮加加, 劉　帥, 沈　毅, 金　銘

(1. 哈爾濱工業大學(威海)信息與電氣工程學院, 山東 威海 264209; 2. 哈爾濱工業大學電子與信息工程學院, 黑龍江 哈爾濱 150001)

0　引　言

波達方向(direction-of-arrival, DOA)估計是陣列信號處理的一個重要研究方向,其在雷達[1]、通信[2]、聲吶[3]、地震[4]和天文[5]等科技領域得到廣泛應用。以多重信號分類(multiple signal classification, MUSIC)[6-9]和旋轉不變子空間[10-11]為代表的子空間類估計算法的提出,實現了傳統DOA估計向超分辨DOA估計的跨越式發展。子空間類算法的關鍵在于通過對陣列協方差矩陣進行奇異值分解(singular value decomposition, SVD)或特征值分解(eigenvalue decomposition, EVD)來獲得噪聲子空間或信號子空間。通常,矩陣EVD和SVD運算都涉及龐大的計算量,尤其是在陣列陣元數多的情況下這種劣勢更加明顯,這也阻礙了子空間類算法的工程化進度。

為克服子空間類算法計算量大的問題,中外學者提出了眾多性能可靠的替代性算法。文獻[12]提出一種傳播因子算法(propagator method, PM),通過協方差矩陣的線性運算來構造噪聲子空間,有效避開了對協方差矩陣的EVD或SVD運算;文獻[13]則利用協方差矩陣的任意K(K為信號源數目)行來求得一個低維的投影矩陣,并在此基礎上構造一個低維搜索函數來實現角度估計;文獻[14-16]提出一種多級維納濾波(multistage Wiener filtering,MSWF)算法,該方法通過正交投影來獲取子空間,同樣避開了計算量較大的EVD運算,實現了算法計算量的降低,但該方法獲得的子空間精度不是太高。

文獻[17]提出一種基于L型陣列的快速2維高效子空間算法(computationally efficient subspace algorithm, CESA)。該算法分別利用協方差矩陣的第1行、第1列和主對角線元素構造出用于求解方位角、俯仰角和角度配對的3個向量,并通過重組這3個向量的元素,最終構造出等價的信號子空間[18]。雖然CESA算法是基于L陣提出的2維DOA估計,但該算法中的L陣的每個軸都由均勻線陣(uniform linear array, ULA)構成,因此其同樣適用于基于ULA的1維DOA估計。本文基于CESA的基本思想,提出一種新的基于ULA的JCCM估計算法,所提算法在陣列模型上將完整的ULA均勻劃分成2個子陣,然后求取2個子陣的前向和后向互協方差矩陣,再利用這些信息經過簡單的線性運算重構出替代的信號子空間,從而構造多項式并求根估計出角度。本文最后用若干個實驗仿真來驗證算法的有效性、估計精度和估計速度等性能。

1　陣列結構及數據模型

設K個不相干的窄帶遠場信號s1(t),s2(t),…,sK(t)分別以角度θ1,θ2,…,θK入射到xoy平面內由2N個各向同性且均勻分布的陣元所構成的ULA上,如圖1所示。定義DOA為信號來向與陣列法線的夾角。假設陣列各通道噪聲為加性高斯白噪聲,且各通道噪聲相互獨立并與信號不相關。現將陣列劃分成陣元數目均為N的兩組,分別記為子陣a和子陣b,以第1個陣元為參考陣元,則兩個子陣t時刻的接收數據矢量可分別表示為

A(θ)s(t)+na(t)

(1)

A(θ)Φs(t)+nb(t)

(2)

圖1　ULA模型

2　聯合互協方差矩陣的DOA估計

2.1　算法原理

為抑制加性噪聲的干擾,定義矩陣Rab為子陣a數據與子陣b數據的互協方差,稱為前向互協方差,矩陣Rba為子陣b數據與子陣a數據的互協方差,稱為后向互協方差。根據式(1)和式(2),有

(3)

(4)

式中,RsE{s(t)sH(t)}表示入射信號的自協方差矩陣。因為噪聲矢量na(t)與nb(t)是統計獨立的關系,所以有基于信號si(t)(i=1,2,…,K)互不相關的假設,信號的自協方差矩陣RS是對角陣,且對角線元素與信號能量一一對應,即

(5)

式中,pi(i=1,2,…,K)表示第i個信號的能量。

利用前向互協方差矩陣和后向互協方差矩陣構造聯合互協方差矩陣(joint cross-covariance matrix,JCCM),記為RJCCM,具體表示為

RJCCMRab+Rba=A(θ)RSΦHAH(θ)+A(θ)ΦRSAH(θ)=

A(θ)(RSΦH+ΦRS)AH(θ)=A(θ)ΛAH(θ)

(6)

式中,ΛRSΦH+ΦRS。

由于Rs和Φ均為對角陣,因此矩陣Λ也為對角陣,且可進一步表示為

Λ=RSΦH+ΦRS=

RSΦH+RSΦ=RS(ΦH+Φ)=

(7)

式中,ξipi(ejNαi+e-jNαi)(i=1,2,…,K)。因為信號能量pi為正實數,所以ξi也為正實數,進一步有Λ=Λ*。

(8)

JNA*(θ)Λ1

(9)

(10)

(11)

(12)

G=CΛD

(13)

式中,C是(2N-K)×K維矩陣;D是K×K維矩陣,且

(15)

(16)

式中

(17)

(18)

實際應用中,協方差矩陣的計算公式為

(19)

所以式(3)、式(4)的實現公式為

(20)

(21)

(22)

(23)

(24)

(25)

(26)

2.2　算法實現步驟描述

綜上所述,本文的JCCM算法具體實現步驟如下:

步驟6由式(16)構造求根多項式fJCCM(z),求其K個最接近于單位圓的根,由式(18)求得信號入射方向。

從上述步驟可見,本文提出的JCCM算法利用子陣劃分和矩陣重構思想,通過對低維(N×1維)聯合互協方差矩陣RJCCM進行線性操作得到等效的高維((2N-K)×K維)信號子空間,避免了直接對全陣列2N×1維協方差矩陣進行EVD或SVD操作,實現了計算量的簡化。

2.3　算法復雜度分析

對陣元數為2N的均勻線陣分別運用本文的JCCM算法和root-MUSIC算法進行DOA估計,設快拍數為L,入射信號源數目為K。

常規root-MUSIC算法的計算量主要在于協方差矩陣計算、特征值分解和多項式求根,這3步操作的計算量分別為O(4N2L)、O(8N3)和O(8(2N-1)3),所以root-MUSIC算法總的計算量為O(4N2L+8N3+8(2N-1)3)。

由于通常情況下快拍數L、陣元數N及信號源數K滿足L?N>K,所以JCCM算法相對于root-MUSIC算法計算量顯著降低。

3　仿真實驗及分析

進行仿真實驗以驗證本文提出的JCCM算法的性能。仿真中,陣列結構采用圖1所示ULA模型,總陣元數2N=16,即各子陣包含8個陣元,入射信號源數K=2,入射角度分別為θ1=-15°,θ2=0°,陣元間距為半波長。實驗還作出了root-MUSIC算法及CESA算法的性能曲線作為參考。

實驗1DOA估計的有效性仿真

實驗條件:快拍數固定為500,信噪比(signal-to-noise ratio,SNR)由-5 dB增加到10 dB。蒙特卡羅實驗次數為100次。仿真得出的成功概率與SNR的關系如圖2所示。

圖2表明,在SNR由-5 dB增加到10 dB的過程中,JCCM與root-MUSIC的成功概率相近,而CESA的成功概率始終不如JCCM及root-MUSIC。說明本文的JCCM算法提高了算法的穩定性。

圖2　成功概率與SNR的關系

實驗2均方根誤差的仿真

算法的均方根誤差(root-mean-square error, RMSE)隨SNR和快拍數的變化情況可用來反映不同情況下的算法估計精度。定義RMSE為

(27)

(1) RMSE隨SNR的變化

仿真條件:快拍數固定為500,SNR由0 dB增加到15 dB。蒙特卡羅實驗次數為500次。仿真結果如圖3所示,其中圖3(a)為3種算法共同比較的結果,為突出顯示JCCM與root-MUSIC的曲線差異,圖3(b)給出了JCCM與root-MUSIC單獨比較的結果。

圖3　RMSE與SNR的關系

由圖3可以看出,在整個SNR變化范圍內,JCCM的RMSE明顯低于CESA,但是始終高于root-MUSIC。從理論的角度分析,JCCM的估計精度高于CESA主要得益于JCCM算法同時利用了前向、后向互協方差,而CESA算法只求解了前向互協方差;然而,JCCM只取聯合互協方差矩陣的第1列構造信號子空間,相對于root-MUSIC利用了整個協方差矩陣獲取子空間而言,還是丟失了一定的有用信息,所以估計精度也相應降低。事實上,從圖3(b)可以看出,當SNR大于3 dB時,JCCM算法可以達到RMSE低于0.1°的估計精度,這樣的估計精度在多數工程應用中是可被接受的。

(2) RMSE隨快拍數的變化

仿真條件:信噪比固定為10 dB,快拍數從100增加到1 000。蒙特卡羅實驗次數為500次。仿真結果如圖4所示,其中圖4(a)為3種算法共同比較的結果,與圖3(b)同理,圖4(b)給出了JCCM與root-MUSIC單獨比較的結果。

圖4　RMSE與快拍數的關系

圖4反映的RMSE隨快拍數變化的趨勢與圖3反映的變化趨勢基本一致,即在整個快拍數變化范圍內,JCCM算法的RMSE介于CESA和root-MUSIC之間,且要明顯低于CESA,雖然高于root-MUSIC,但在快拍數大于200時便已小于0.1°,同樣已經可以達到多數工程應用的精度要求了。

實驗3角度估計的速度仿真

角度估計速度使用算法在不同快拍數下的單次估計時間來衡量。實驗中取蒙特卡羅實驗運行時間的平均值作為算法的單次估計時間。

仿真條件:SNR固定為10 dB,快拍數由100變化到2 000,蒙特卡羅實驗次數為500。估計時間與快拍數關系的仿真結果如圖5所示。

圖5　估計時間與快拍數的關系

由圖5可以看出,在整個快拍數變化范圍內,JCCM與CESA算法的單次平均估計時間始終小于root-MUSIC算法,仿真結果驗證了算法估計速度上的優勢。需要說明的是,圖5顯示出小快拍數下估計時間明顯小于大快拍數下的估計時間,從理論角度分析,這是由于快拍數較小時求得的協方差矩陣并不滿秩,從而特征值分解運算量小且求根多項式階次低,所以總體運算量較小。當快拍數較大時,協方差矩陣近似為滿秩,此時特征值分解運算量及求根多項式的階次均與陣元數成正比,而不受快拍數影響,快拍數僅影響求解協方差矩陣的運算量,因此3種算法隨快拍數變化的趨勢相一致。

4　結　論

本文提出一種基于ULA的聯合互協方差矩陣的DOA估計算法。算法避開繁瑣的子空間分解過程,利用子陣接收數據的前向和后向互協方差矩陣構造了聯合互協方差矩陣,并利用其重構出了等價的信號子空間,最后構造多項式通過求根實現了入射角度的估計。結果表明,本文的JCCM算法在保證估計精度可接受的同時,有效降低了計算量,實現了估計速度的提高。

參考文獻：

[1] WANG X, WANG L, LI X, et al. An efficient sparse representation algorithm for DOA estimation in MIMO radar system[C]∥Proc.of the IEEE International Workshop on Signal Processing Advances in Wireless Communications, 2016:1-4.

[2] LIU L, LIU H. Joint estimation of DOA and TDOA of multiple reflections in mobile communications[J]. IEEE Access, 2016, 4:3815-3823.

[3] SAUCAN A A, CHONAVEL T, SINTES C, et al. CPHD-DOA tracking of multiple extended sonar targets in impulsive environments[J].IEEE Trans.on Signal Processing,2016,64(5):1147-1160.

[4] YAN G, LI G, NING L. Method on fast DOA estimation of moving nodes in ad-hoc network[C]∥Proc.of the IEEE International Symposium on Communications and Information Technology, 2006:1169-1172.

[5] LEVANDA R, LESHEM A. Adaptive selective sidelobe canceller beamformer with applications to interference mitigation in radio astronomy[J].IEEE Trans.on Signal Processing,2013, 61(20):5063-5074.

[6] SANTOSH S, SHARMA K. A review on multiple emitter location and signal parameter estimation[J]. International Journal of Engineering Research, 2013, 2(3):276-280.

[7] YAN F G, JIN M, LIU S, et al. Real-valued music for efficient direction estimation with arbitrary array geometries[J]. IEEE Trans.on Signal Processing, 2014, 62(6):1548-1560.

[8] BASIKOLO T, ARAI H. APRD-MUSIC algorithm DOA estimation for reactance based uniform circular array[J]. IEEE Trans.on Antennas & Propagation,2016,64(10):4415-4422.

[9] 閆鋒剛,張薇,金銘.求根MUSIC初值設置和更新算法[J].哈爾濱工業大學學報,2015, 47(3): 88-92.

YAN F G, ZHANG W, JIN M. A new method for setting and updating the initiation of root-MUSIC[J]. Journal of Harbin Institute of Technology, 2015, 47(3): 88-92.

[10] ROY R, PAULRAJ A, KAILATH T. ESPRIT: a subspace rotation approach to estimation of parameters of cissoids in noise[J]. IEEE Trans.on Acoustics,Speech,and Signal Processing，1986, 34(5): 1340-1342.

[11] LIN J, MA X, YAN S, et al. Time-frequency multi-invariance ESPRIT for DOA estimation[J]. IEEE Antennas & Wireless Propagation Letters, 2016, 15(1):770-773.

[12] MARCOS S, MARSAL A, BENIDIR M. Performances analysis of the propagator method for source bearing estimation[C]∥Proc.of the IEEE International Conference on Acoustics, Speech, & Signal Processing,1994:237-240.

[13] YEH C C. Simple computation of projection matrix for bearing estimations[J]. IEE Proceedings F-Communications, Radar and Signal Processing, 1987, 134(2):146-150.

[14] GOLDSTEIN J S, REED I S, SCHARF L L. A multistage representation of the Wiener filter based on orthogonal projections[J].IEEE Trans.on Information Theory, 1998, 44(7):2943-2959.

[15] HUANG L, WU S, FENG D, et al. Low complexity method for signal subspace fitting[J]. Electronics Letters, 2004, 40(14): 847-848.

[16] HUANG L, WU S. Low-complexity MDL method for accurate source enumeration[J]. IEEE Signal Processing Letters, 2007, 14(9): 581-584.

[17] XI N, LI L. A computationally efficient subspace algorithm for 2-D DOA estimation with l-shaped array[J]. IEEE Signal Processing Letters, 2014, 21(8):971-974.

[18] 榮加加.特殊陣列下的快速波達方向估計[D].哈爾濱：哈爾濱工業大學,2017:19-32.

RONG J J. Fast direction of arrival estimation based on special array configurations[D].Harbin:Harbin Institute of Technology, 2017:19-32.