基于矢量化差分相位的單分布源解耦二維波達方向估計?

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

(解放軍信息工程大學信息系統工程學院,鄭州 450001)

1 引 言

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

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

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

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

2 數學模型

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

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

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

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

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

根據單個分布源內部不同路徑分量是否相關,可以將分布源分成兩類:CD源和ID源,分別對應于多徑時不變和多徑時變兩種信道狀態.當信道的相關時間大于觀測周期時,此時分布源為CD源;而當信道的相關時間遠小于觀測周期時,此時分布源為ID源.本文討論單個CD源和單個ID源入射兩種情況,因此接下來需要對上述分布源統一模型進行區別討論.

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

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

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

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

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

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

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

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

首先對于CD源,由于信號與噪聲無關,根據(5)式的廣義陣列流型,可以得到對應于CD源的陣列協方差矩陣RCD的解析表達式為

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

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

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

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

其中,?表示Schur-Hadamard積.

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

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

對于ID源,根據(2)式和(6)式,可以計算對應的協方差矩陣為

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

其中

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

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

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

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

進而有?p,q ∈[?π,π),因此?p,q不會產生相位模糊.應該指出的是當r/λ>1/4時,相關學者也提出了多種方法來解決相位模糊問題[26?28].

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

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

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

其中

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

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

3.3 算法步驟與復雜度分析

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

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

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

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

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

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

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

4 仿真實驗

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

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

其中,

Q

為蒙特卡羅仿真次數,

q

和

q

分別為第

q

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

仿真1參數估計直方圖分布

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

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

仿真2不同算法性能對比

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

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

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

圖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估計均方根誤差與快拍數的關系 (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的估計能力

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

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

5 結 論

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

附錄A

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

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

[1]Liang G L,Ma W,Fan Z,Wang Y L 2013Acta Phys.Sin.62 144302(in Chinese)[梁國龍,馬巍,范展,王逸林2013物理學報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物理學報64 078403]

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

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

[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)[楊學敏,李廣軍,鄭植2014電子與信息學報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