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

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

0　引　言

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

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

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

1　陣列結(jié)構(gòu)及數(shù)據(jù)模型

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

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

(1)

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

(2)

圖1　ULA模型

2　聯(lián)合互協(xié)方差矩陣的DOA估計

2.1　算法原理

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

(3)

(4)

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

(5)

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

利用前向互協(xié)方差矩陣和后向互協(xié)方差矩陣構(gòu)造聯(lián)合互協(xié)方差矩陣(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為正實數(shù),所以ξi也為正實數(shù),進一步有Λ=Λ*。

(8)

JNA*(θ)Λ1

(9)

(10)

(11)

(12)

G=CΛD

(13)

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

(15)

(16)

式中

(17)

(18)

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

(19)

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

(20)

(21)

(22)

(23)

(24)

(25)

(26)

2.2　算法實現(xiàn)步驟描述

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

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

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

2.3　算法復雜度分析

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

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

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

3　仿真實驗及分析

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

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

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

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

圖2　成功概率與SNR的關(guān)系

實驗2均方根誤差的仿真

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

(27)

(1) RMSE隨SNR的變化

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

圖3　RMSE與SNR的關(guān)系

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

(2) RMSE隨快拍數(shù)的變化

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

圖4　RMSE與快拍數(shù)的關(guān)系

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

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

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

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

圖5　估計時間與快拍數(shù)的關(guān)系

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

4　結(jié)　論

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

參考文獻：

[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].哈爾濱工業(yè)大學學報,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].哈爾濱：哈爾濱工業(yè)大學,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.