代正亮崔維嘉 王大鳴 張彥奎

(解放軍信息工程大學(xué)信息系統(tǒng)工程學(xué)院,鄭州 450001)

1 引 言

在雷達(dá)、聲納和移動通信等應(yīng)用領(lǐng)域中,發(fā)射的信號時常會受到復(fù)雜環(huán)境散射和反射等影響,使得天線陣列接收到的信號能量在空間某一范圍呈現(xiàn)一定分布[1,2].此時,基于點源模型假設(shè)的傳統(tǒng)波達(dá)方向(direction of arrival,DOA)估計算法性能會嚴(yán)重惡化,而需要采用參數(shù)化的分布式信源模型進(jìn)行處理[3?5].根據(jù)同一信源不同路徑分量是否相關(guān),分布源主要包括兩類:相干分布源(coherently distributed source,CD)和非相干分布源(incoherently distributed source,ID)[6].對于CD源,同一信源不同路徑分量僅相差一個固定的相位延遲和幅度加權(quán);與之相反,同一ID源的不同路徑分量則完全不相關(guān).

近年來,分布源的DOA估計研究得到了廣泛關(guān)注,相關(guān)學(xué)者陸續(xù)提出了一系列有效的估計算法,比如子空間類算法[7?10]、波束形成類算法[11]、最大似然類算法[12,13]和稀疏重構(gòu)類算法[14,15]等.然而這些算法均只適用于一維分布源,在實際應(yīng)用中,信號源與接收陣列往往不在同一平面上,這種情況下需要將其建立為一個二維分布源模型.一個二維分布源主要包括四個未知角度參數(shù):中心方位角,方位角擴展,中心俯仰角和俯仰角擴展.盡管上述的某些一維估計方法可以直接推廣到二維,如二維協(xié)方差匹配算法[16]等,這些算法由于需要多維參數(shù)搜索復(fù)雜度較高.為解決此類問題,相關(guān)學(xué)者基于一些特定的陣列結(jié)構(gòu)研究了分布源的近似簡化模型,并提出了一些低復(fù)雜度的二維DOA估計算法,如基于三平行線陣的二維傳播算子算法[17,18],基于矩形陣的二維旋轉(zhuǎn)不變子空間(rotation invariant subspace,ESPRIT)算法[19]和基于對稱十字陣的二維廣義ESPRIT算法[20]等.

特別地,均勻圓陣(uniform circular array,UCA)相比于其他陣列結(jié)構(gòu)具有360?的全方位角覆蓋、幾乎不變的方向圖以及額外的俯仰角信息等優(yōu)點[21],基于UCA的分布源二維DOA估計成為了研究熱點[22?25].文獻(xiàn)[22]首次基于空間距離相近的兩個UCA,提出了一種估計二維CD源DOA的一維交替搜索(sequential one-dimensional searching,SOS)算法,該算法首先基于一階泰勒(Taylor)級數(shù)展開得到兩個子陣間的近似旋轉(zhuǎn)不變關(guān)系,接下來利用總體最小二乘旋轉(zhuǎn)不變子空間(total least squares rotation invariant subspace,TLS-ESPRIT)算法得到中心俯仰角估計,最后通過多次一維搜索得到中心方位角估計.算法只需要一維搜索,但俯仰角的初始估計精度對算法性能影響較大,并且靠得很近的兩個子陣容易產(chǎn)生互耦誤差,降低了估計精度.在此基礎(chǔ)上,文獻(xiàn)[23]僅利用單個對稱UCA,提出了一種基于廣義ESPRIT的二維DOA估計算法,該算法無須角信號分布函數(shù)先驗已知,但需要二維搜索,復(fù)雜度較高.文獻(xiàn)[22,23]中的算法都是針對CD源的研究,而對于UCA中的ID源二維DOA估計問題,相關(guān)學(xué)者也取得了一些研究成果.文獻(xiàn)[24]同樣基于靠得很近的兩個UCA,首先根據(jù)子陣間的近似旋轉(zhuǎn)不變關(guān)系得到中心俯仰角估計,進(jìn)而在單個UCA中,通過廣義多重信號分類(generalized multiple signal classif i cation,GMUSIC)方法實現(xiàn)了中心方位角估計.而在文獻(xiàn)[25]中,考慮由多個UCA組成的均勻圓柱陣列(uniform cylindrical array,UCyA),提出了一種低復(fù)雜度的波束空間二維DOA估計算法,該算法估計中心俯仰角的方法與文獻(xiàn)[24]類似,但是在波束域進(jìn)行計算,并且通過波束域GMUSIC(beamspace GMUSIC,BS-GMUSIC)方法得到了中心方位角估計,具有較低的計算復(fù)雜度.盡管如此,在上述文獻(xiàn)[22—25]的算法中,復(fù)雜度較高的譜峰搜索和特征值分解運算依舊無法避免,所以復(fù)雜度仍有待降低.此外,文獻(xiàn)[22,24,25]中的算法均需要多個UCA,陣列結(jié)構(gòu)的硬件復(fù)雜度較高,不利于工程實踐.

本文針對單個分布源入射的情況,提出了一種基于矢量化差分相位的解耦二維DOA快速估計算法.所提算法首先基于UCA不同陣元間的差分相位不受角度擴展影響的特性,實現(xiàn)了中心DOA的解耦合;進(jìn)而通過最小二乘方法直接得到二維DOA的閉式解.該算法僅需要單個均勻圓陣,無須任何的譜峰搜索和特征值分解運算,具有較低的硬件復(fù)雜度和計算復(fù)雜度.仿真實驗表明,相比于現(xiàn)有的低復(fù)雜度算法,所提算法具有較好的估計性能;同時,該算法無須角信號分布的任何先驗信息.

2 數(shù)學(xué)模型

考慮如圖1所示的均勻圓陣,該陣列由分布于x-y平面的M個陣元組成,各陣元均勻分布在半徑為r的圓周上,第一個陣元位于x軸上且坐標(biāo)原點位于圓陣的圓心.假設(shè)有單個遠(yuǎn)場窄帶分布源入射到該陣列,則t時刻該陣列第m個陣元的輸出信號可以表示為[22?25]

其中γm=2π(m?1)/M,λ為入射信號波長;t=1,2,………,T為采樣時刻,T為快拍總數(shù);s(θ,φ,t;μ)為分布源的角信號密度函數(shù);μ=(θ0,σθ0,φ0,σφ0)為分布源的角度參數(shù)矢量,其中各分量分別表示中心方位角θ0,方位角擴展σθ0,中心俯仰角φ0和俯仰角擴展σφ0;nm(t)是均值為0,方差為σ2n的高斯白噪聲.需要注意的是,為了表達(dá)式的簡化,在下面的相關(guān)積分中均省略了上下限的標(biāo)識.

圖1 均勻圓陣Fig.1.Uniform circular array.

進(jìn)一步可將(1)式寫成向量形式為

其中x(t)和n(t)分別為M×1的輸出信號矢量和噪聲矢量;a(θ,φ)為均勻圓陣的陣列方向矢量,

根據(jù)單個分布源內(nèi)部不同路徑分量是否相關(guān),可以將分布源分成兩類:CD源和ID源,分別對應(yīng)于多徑時不變和多徑時變兩種信道狀態(tài).當(dāng)信道的相關(guān)時間大于觀測周期時,此時分布源為CD源;而當(dāng)信道的相關(guān)時間遠(yuǎn)小于觀測周期時,此時分布源為ID源.本文討論單個CD源和單個ID源入射兩種情況,因此接下來需要對上述分布源統(tǒng)一模型進(jìn)行區(qū)別討論.

首先,對于一個二維CD源,由于不同的入射分量僅相差一個固定的相位延遲和幅度加權(quán),因此角信號密度函數(shù)s(θ,φ,t;μ)可進(jìn)一步表示為其中,s(t)是復(fù)隨機信號源;?(θ,φ;μ)為相應(yīng)的確定性角信號密度函數(shù),通常假設(shè)為以中心DOA為對稱中心的單峰對稱函數(shù),服從均勻分布、高斯分布和拉普拉斯分布[22]等.

接下來將(4)式代入(2)式,可以得到

然而,對于二維ID源而言,由于不同的路徑分量是完全不相關(guān)的,則有下式成立:

3 基于矢量化差分相位的二維DOA快速估計算法

在分布源二維DOA估計中,現(xiàn)有基于均勻圓陣的估計算法均需要譜峰搜索和特征分解運算,復(fù)雜度較高.針對該問題,本文分別考慮單個CD源和單個ID源入射到均勻圓陣這兩種情形,首先基于空間頻率近似模型,證明了各陣元接收信號間的差分相位可實現(xiàn)中心DOA的解耦合,即差分相位中只包含中心DOA信息,而與角度擴展無關(guān);然后基于此特性,通過對協(xié)方差矩陣的嚴(yán)格上三角元素相位進(jìn)行矢量化處理,使得二維DOA估計問題轉(zhuǎn)化為一個最小二乘問題,進(jìn)而直接得到中心方位角和俯仰角的閉式解.

3.1 陣元間差分相位的解耦合特性證明

根據(jù)(2)式的陣列輸出信號矢量x(t),可以得到相應(yīng)的輸出信號協(xié)方差矩陣為R=E[x(t)xH(t)],則R的第(p,q)個元素可表示為

由上式分析可知[R]p,q的相位即為第p個陣元和第q個陣元輸出信號的差分相位,所以可以通過分析協(xié)方差矩陣中的元素相位來實現(xiàn)陣元間差分相位的解耦合特性證明.由于CD源和ID源的陣列模型有所不同,因此需要分別進(jìn)行討論.

首先對于CD源,由于信號與噪聲無關(guān),根據(jù)(5)式的廣義陣列流型,可以得到對應(yīng)于CD源的陣列協(xié)方差矩陣RCD的解析表達(dá)式為

其中,IM為M×M的單位矩陣.由于協(xié)方差矩陣RCD的主對角線元素中不包含中心DOA信息,因此只需要考慮非主對角線上的元素,即當(dāng)p?=q時,[RCD]p,q可進(jìn)一步表示為

可以定義θ=θ0+?θ和φ=φ0+?φ,其中?θ和?φ分別是θ,φ和中心方位角θ0、中心俯仰角φ0之間的角度偏差.在小角度擴展假設(shè)下,基于空間頻率近似模型[7?10],有以下近似關(guān)系成立:

在廣義方向矢量b(μ)中,中心DOA與角度擴展相互耦合.根據(jù)(10),(11)和(12)式,則可以得到廣義方向矢量b(μ)的參數(shù)解耦形式:

其中g(shù)(μ)為確定性角信號分布函數(shù)矢量,且有

其中,?表示Schur-Hadamard積.

將(15)式代入(9)式,可得

分析上式,雖然g(μ)中同時包含中心DOA和角度擴展信息,但由于其為實向量,不會對[RCD]p,q的相位產(chǎn)生影響,因此[RCD]p,q的相位可以表示為

對于ID源,根據(jù)(2)式和(6)式,可以計算對應(yīng)的協(xié)方差矩陣為

在無噪?yún)f(xié)方差矩陣Rs(μ)中,中心DOA與角度擴展參數(shù)相互耦合.在小角度擴展下,基于(10),(11)和(12)式的近似關(guān)系,同樣可以得到Rs(μ)的參數(shù)解耦形式:

其中

由(17)和(22)式可知,當(dāng)p ?=q時,[RCD]p,q和[RID]p,q的相位中均僅包含待定的分布源中心DOA,而不包含任何的角度擴展參數(shù),因此通過提取協(xié)方差矩陣的非對角線元素相位可將中心DOA與角度擴展分離開,即證明了均勻圓陣中不同陣元接收信號的差分相位具有解耦合的特性.

3.2 中心方位角與俯仰角估計

通過上述的陣元間差分相位解耦合特性的證明,可知協(xié)方差矩陣非對角線上的元素相位僅包含待定的中心DOA信息,基于此,可進(jìn)一步得到中心方位角與俯仰角估計.根據(jù)(17)和(22)式可知,具有相同的形式,因此接下來可以一起討論,統(tǒng)一記為?p,q,于是有:

值得注意的是,為了避免產(chǎn)生相位模糊問題,本文假設(shè)r/λ6 1/4,同時

進(jìn)而有?p,q ∈[?π,π),因此?p,q不會產(chǎn)生相位模糊.應(yīng)該指出的是當(dāng)r/λ>1/4時,相關(guān)學(xué)者也提出了多種方法來解決相位模糊問題[26?28].

為了得到中心方位角和俯仰角估計,可以將(23)式進(jìn)一步展開,改寫為以下形式:

由于協(xié)方差矩陣是共軛對稱的,并且主對角線的元素相位不含方位信息,因此只需要考慮協(xié)方差矩陣嚴(yán)格上三角部分的元素,將這部分元素相位按行提取,并進(jìn)行矢量化處理,進(jìn)而可以得到差分相位矢量W:

接下來可將(24)式改寫為如下矩陣形式:

其中

因此v的最小二乘估計由下式給出:

最后,分別得到單分布源的中心方位角和俯仰角估計:

3.3 算法步驟與復(fù)雜度分析

根據(jù)上述分析可以將本文估計單分布源二維DOA的方法歸納為以下步驟:

1)根據(jù)(2)式的分布源陣列輸出信號矢量,計算相應(yīng)的協(xié)方差矩陣R=E[x(t)xH(t)];

2)提取協(xié)方差矩陣R的嚴(yán)格上三角元素相位,并進(jìn)行矢量化處理,得到差分相位矢量W;

3)根據(jù)(27)和(28)式分別構(gòu)造矩陣A和v,進(jìn)而得到如(26)式所示的等式關(guān)系,接下來根據(jù)(29)式直接得到v的最小二乘估計;

4)由(30)和(31)式計算中心方位角和俯仰角估計值.

接下來將本文算法與文獻(xiàn)[22]中的SOS算法、文獻(xiàn)[23]中的NAM算法、文獻(xiàn)[24]中的GMUSIC算法和文獻(xiàn)[25]中的BS-GMUSIC算法進(jìn)行復(fù)雜度對比分析.假設(shè)陣列陣元數(shù)為M,快拍數(shù)為T,則本文算法的復(fù)雜度主要在于估計一個M階的陣列協(xié)方差矩陣以及最小二乘計算過程,相應(yīng)的計算量可以表示為O(M2T+(M2?M)/2);SOS算法的復(fù)雜度主要在于估計一個M階的陣列協(xié)方差矩陣和特征分解運算以及一維譜峰搜索過程,可以表示為O(M3+9M2T+L(1+M));NAM算法的復(fù)雜度主要在于估計一個M階的陣列協(xié)方差矩陣和特征分解運算以及二維譜峰搜索過程,可以表示為O(M2T+M3+L2(M+1));GMUSIC算法的主要復(fù)雜度同樣在于估計一個M階的陣列協(xié)方差矩陣和特征分解運算以及一維譜峰搜索過程,可以表示為O(M2T+M3+M2L);BS-GMUSIC算法實現(xiàn)了波束域計算,其復(fù)雜度主要包括波束空間轉(zhuǎn)換過程、估計一個M階協(xié)方差矩陣和特征分解運算以及一維波束域譜峰搜索過程,可以表示為O(MPT+P2T+P3+P2L).注意在上述復(fù)雜度分析中,L表示搜索點數(shù),P為輸出信號矩陣波束空間轉(zhuǎn)換后的維數(shù).

由上述的復(fù)雜度分析可知,當(dāng)搜索點數(shù)L較大時,本文算法的復(fù)雜度明顯低于其他算法,這是因為本文算法無須任何的譜峰搜索和特征分解運算.

4 仿真實驗

本文研究的是單個分布源DOA估計算法,擬分別采用單個CD源與單個ID源作為發(fā)射信號.仿真實驗采用如圖1所示的均勻圓陣,圓陣中包含20個陣元,半徑為r=λ/4.實驗中假設(shè)噪聲為高斯白噪聲.為了驗證本文算法的實用性與魯棒性,采用蒙特卡羅實驗將本文算法與文獻(xiàn)[22]中的SOS算法、文獻(xiàn)[23]中的NAM算法、文獻(xiàn)[24]中的GMUSIC算法、文獻(xiàn)[25]中的BS-GMUSIC算法和相應(yīng)的CRLB進(jìn)行對比分析.

為衡量算法性能,首先定義中心方位角和俯仰角估計的二維均方根誤差2D-RMSE為

其中,

Q

為蒙特卡羅仿真次數(shù),

q

和

q

分別為第

q

次蒙特卡羅實驗中分布源的中心方位角和俯仰角估計值.

仿真1參數(shù)估計直方圖分布

圖2 分布源二維DOA估計直方圖 (a)CD源;(b)ID源Fig.2.Histogram of two dimensional DOA estimation for distributed source:(a)CD source;(b)ID source.

假設(shè)入射的單個分布源角度參數(shù)矢量為μ=(30? ,1? ,50? ,1?),固定信噪比(SNR)為10 dB,快拍數(shù)為500,分別在CD源和ID源兩種情況下,進(jìn)行300次蒙特卡羅實驗,相應(yīng)的二維DOA估計值的分布直方圖如圖2所示.從圖中可以看出,中心方位角和中心俯仰角估計的誤差范圍均在1?以內(nèi),并且越接近真值,估計次數(shù)越多,因此本文算法能夠?qū)崿F(xiàn)對單個分布源二維DOA的有效估計.

仿真2不同算法性能對比

在CD源的條件下,將本文算法與SOS算法和NAM算法性能進(jìn)行對比,而對于ID源,則將本文算法與GMUSIC算法和BS-GMUSIC算法性能進(jìn)行對比.入射的CD源或ID源的角度參數(shù)矢量均為μ=(30? ,2? ,60? ,3?),SNR從0 dB變化到30 dB,快拍數(shù)為500,分別繪制這些算法的2DRMSE曲線與信噪比SNR的關(guān)系曲線如圖3所示.從圖中可以看出,在CD源角度估計中,本文算法的估計性能明顯優(yōu)于SOS算法,而與NAM算法近乎相同;在ID源角度估計,本文算法的性能在低信噪比時明顯優(yōu)于GMUSIC算法和BS-GMUSIC算法.因此結(jié)合上述的復(fù)雜度分析可得出結(jié)論:所提算法在具有較低復(fù)雜度的同時也具有一個較好的估計性能.

仿真3驗證快拍數(shù)對算法性能的影響

假設(shè)快拍數(shù)由100到1500之間變化,SNR固定為10 dB,其他參數(shù)與仿真二相同,分別繪制不同算法的2D-RMSE曲線與快拍數(shù)的關(guān)系,如圖4所示.從圖中可以得出與仿真二類似的結(jié)論,并且隨著快拍數(shù)的增多,算法的估計性能也越來越好.

圖3 不同算法二維DOA均方根誤差對比 (a)CD源;(b)ID源Fig.3.2D-RMSE contrast of two-dimensional DOA of different algorithms:(a)CD source;(b)ID source.

圖4 不同算法二維DOA估計均方根誤差與快拍數(shù)的關(guān)系 (a)CD源;(b)ID源Fig.4.The relationship between 2D-RMSE of two-dimensional DOA and the number of snapshots of different algorithms:(a)CD source;(b)ID source.

仿真4驗證算法對不同中心DOA的估計能力

選取分布源的角度擴展參數(shù)為σθ0=σφ0=1.5?, 首先固定中心俯仰角φ0=60?,SNR=10 dB,快拍數(shù),研究中心方位角的變化對算法估計性能的影響,如圖5(a)所示.從圖中可以看出,在全方位角范圍內(nèi)[?180? ,180?),所提算法均具有近乎相同的估計性能與分辨力,因此本文算法可以實現(xiàn)全方位測角;類似地,固定中心方位角θ0=30?,則中心俯仰角的變化對算法估計性能的影響如圖5(b)所示,分析可知當(dāng)中心俯仰角接近角度域下界φ0=0?或者上界φ0=90?時,估計的均方根誤差顯著增大但仍然小于1?,這說明算法在整個角度區(qū)域(θ0∈[?180? ,180?),φ0∈(0? ,90?))內(nèi)都可進(jìn)行有效估計.

圖5 算法對不同中心DOA的估計能力 (a)中心方位角;(b)中心俯仰角Fig.5.Estimation performance of the algorithm for different center DOAs:(a)Center azimuth;(b)center elevation.

5 結(jié) 論

基于均勻圓陣的分布源二維DOA估計可以實現(xiàn)全方位測角,并且具有較高的分辨率,然而現(xiàn)有的估計算法均需要譜峰搜索和特征值分解,復(fù)雜度較高.針對該問題,本文算法在證明任意單分布源入射時,均勻圓陣不同陣元接收信號間的差分相位具有解耦合特性的基礎(chǔ)上,通過提取協(xié)方差矩陣的嚴(yán)格上三角元素相位,并進(jìn)行矢量化處理,最終將波達(dá)方向估計問題轉(zhuǎn)化為一個最小二乘問題,從而可以直接得到閉式解.通過仿真實驗和復(fù)雜度分析可以看出,所提算法與現(xiàn)有的算法相比,在大幅度降低復(fù)雜度的同時,依舊能夠獲得較好的估計性能;并且僅需要單個圓陣,具有較低的硬件復(fù)雜度,同時也無須角度分布的任何先驗信息.

附錄A

因此g(μ)為實值矢量.

因此B(μ)為實對稱矩陣.

[1]Liang G L,Ma W,Fan Z,Wang Y L 2013Acta Phys.Sin.62 144302(in Chinese)[梁國龍,馬巍,范展,王逸林2013物理學(xué)報62 144302]

[2]Xiong W,Picheral J,Marcss S 2017Digital Signal Process.63 155

[3]Ba B,Liu G C,Fan Z,Li T,Fan Z,Lin Y C,Wang Y 2015Acta Phys.Sin.64 078403(in Chinese)[巴斌, 劉國春,李韜,范展,林禹丞,王瑜2015物理學(xué)報64 078403]

[4]Jiang H,Zhou J,Hisakazu K,Shao G F 2014Acta Phys.Sin.63 048702(in Chinese)[江浩,周杰,菊池久和,邵根富2014物理學(xué)報63 048702]

[5]Zheng Z 2011Ph.D.Dissertation(Chengdu:University of Electronic Science and Technology)(in Chinese)[鄭植2011博士學(xué)位論文(成都:電子科技大學(xué))]

[6]Valaee S,Champagne B,Kabal P 1995IEEE Trans.Signal Process.43 2144

[7]Shahbazpanahi S,Valaee S,Bastani M H 2001IEEE Trans.Signal Process.49 2169

[8]Zheng Z,Li G 2013Multi.Sys.Signal Process.24 573

[9]Yang X M,Li G J,Chi C K,Zheng Z,Yeo T S 2015Circ.Sys.Signal Process.34 3697

[10]Cao R Z,Gao F,Zhang X 2016IEEE Trans.Signal Process.64 1

[11]Hassanien A,Shahbazpanahi S,Gershman A B 2004IEEE Trans.Signal Process.52 280

[12]Shahbazpanahi S,Valaee S,Gershman A B 2004IEEE Trans.Signal Process.52 592

[13]Sieskul B T 2010IEEE Trans.Vehicul.Technol.59 1534

[14]Yang X,Li G J,Zheng Z 2014J.Electron.Informat.Technol.36 164(in Chinese)[楊學(xué)敏,李廣軍,鄭植2014電子與信息學(xué)報36 164]

[15]Yang X M,Zheng Z,Hu B 2016Electro.Lett.52 262

[16]Boujemaa H 2005European Trans.Telecommun.16 557

[17]Zheng Z,Li G,Teng Y 2012Circ.Sys.Signal Process.31 255

[18]Dai Z L,Ba B,Cui W J,Sun Y M 2017IEEE Acce.99 1

[19]Hu A,Lü T,Gao H,Zhang Z,Yang S 2014IEEE J.Sel.Topics in Signal Process.8 996

[20]Dai Z L,Cui W J,Ba B,Wang D M,Sun Y M 2017Sens.17 1300

[21]Cao M Y,Huang L,Qian C,Xue J Y,So H C 2015Signal Process.106 41

[22]Lee J,Song I,Kwon H,Lee S R 2003Signal Process.83 1789

[23]Nam J G,Lee S H,Lee K K 2014IEEE Anten.Wire.Propag.Lett.13 415

[24]Guo X,Wan Q,Shen X,Dou H 2011Tur.J.Elec.Enging.Com.Sci.19 445

[25]Lü T,Tan F,Gao H,Yang G 2016Signal Process.121 30

[26]Sundaram K R,MalliK R K,Murthy U M S 2000IEEE Trans.on Aero.Electro.Sys.36 1391

[27]Ballal T,Blealley C J 2008IEEE Signal Process.Lett.15 853

[28]Chen X,Liu Z,Wei X 2016IEEE Anten.Wire.Propag.Lett.99 1