聯(lián)合互協(xié)方差矩陣的快速波達(dá)方向估計(jì)

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

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

0　引　言

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

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

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

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

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

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

(1)

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

(2)

圖1　ULA模型

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

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)}表示入射信號(hào)的自協(xié)方差矩陣。因?yàn)樵肼暿噶縩a(t)與nb(t)是統(tǒng)計(jì)獨(dú)立的關(guān)系,所以有基于信號(hào)si(t)(i=1,2,…,K)互不相關(guān)的假設(shè),信號(hào)的自協(xié)方差矩陣RS是對(duì)角陣,且對(duì)角線元素與信號(hào)能量一一對(duì)應(yīng),即

(5)

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

利用前向互協(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和Φ均為對(duì)角陣,因此矩陣Λ也為對(duì)角陣,且可進(jìn)一步表示為

Λ=RSΦH+ΦRS=

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

(7)

式中,ξipi(ejNαi+e-jNαi)(i=1,2,…,K)。因?yàn)樾盘?hào)能量pi為正實(shí)數(shù),所以ξi也為正實(shí)數(shù),進(jìn)一步有Λ=Λ*。

(8)

JNA*(θ)Λ1

(9)

(10)

(11)

(12)

G=CΛD

(13)

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

(15)

(16)

式中

(17)

(18)

實(shí)際應(yīng)用中,協(xié)方差矩陣的計(jì)算公式為

(19)

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

(20)

(21)

(22)

(23)

(24)

(25)

(26)

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

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

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

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

2.3　算法復(fù)雜度分析

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

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

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

3　仿真實(shí)驗(yàn)及分析

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

實(shí)驗(yàn)1DOA估計(jì)的有效性仿真

實(shí)驗(yàn)條件:快拍數(shù)固定為500,信噪比(signal-to-noise ratio,SNR)由-5 dB增加到10 dB。蒙特卡羅實(shí)驗(yàn)次數(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)系

實(shí)驗(yàn)2均方根誤差的仿真

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

(27)

(1) RMSE隨SNR的變化

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

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

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

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

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

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

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

實(shí)驗(yàn)3角度估計(jì)的速度仿真

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

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

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

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

4　結(jié)　論

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

參考文獻(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初值設(shè)置和更新算法[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào),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á)方向估計(jì)[D].哈爾濱：哈爾濱工業(yè)大學(xué),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.