饒偉，桂宇風(fēng)，李旦

（1.南昌工程學(xué)院信息工程學(xué)院，江西 南昌 330099；2.復(fù)旦大學(xué)信息科學(xué)與工程學(xué)院，上海 200433）

1 引言

波達(dá)方向（DoA,direction of arrival）估計作為陣列信號處理領(lǐng)域的一個主要研究方向，廣泛應(yīng)用于通信、雷達(dá)、聲吶等領(lǐng)域[1]。絕大多數(shù)子空間類DoA 估計方法[2-4]最初是針對均勻線性陣列（ULA,uniform linear array）結(jié)構(gòu)提出的，且為避免出現(xiàn)角度模糊問題，相鄰陣元間距需小于或等于載波半波長。在物理陣元總數(shù)和陣列結(jié)構(gòu)確定的情況下，相鄰陣元間距小意味著陣列孔徑小，此時難以獲得更好的DoA 估計性能。為此，人們提出了一些具有高自由度或大孔徑的非均勻線性陣列結(jié)構(gòu)，例如嵌套陣列[5]和互質(zhì)陣列[6-8]。嵌套陣列中包含一個相鄰陣元間距較小的子陣（稱為密集子陣），其存在較嚴(yán)重的互耦效應(yīng)，這給信號參數(shù)估計帶來了一定的負(fù)面影響。為解決這個問題，文獻(xiàn)[6]提出了一種互質(zhì)陣列，它由陣元數(shù)分別為M和N的2 個ULA 組成，且這2 個ULA 中相鄰陣元間距分別為和，其中，λ為載波波長，M和N為互質(zhì)整數(shù)。特別地，文獻(xiàn)[6]證明了僅使用M+N? 1個物理陣元便可獲得O(MN) 的自由度（DoF,degree of freedom）。隨后，文獻(xiàn)[7-8]對互質(zhì)陣列進(jìn)行了更深入的研究，從而進(jìn)一步提升了陣列自由度。

上述均勻和非均勻線性陣列均為一維陣列，一般情況下只能估計信號一維DoA。為實現(xiàn)二維DoA估計，眾多學(xué)者將針對ULA 的傳統(tǒng)子空間類方法推廣到二維領(lǐng)域，提出了二維MUSIC（multiple signal classification）算法、二維ESPRIT（estimation of signal parameters via rotational invariance technique）算法、二維傳播算子（PM,propagator method）[9-12]等，且主要應(yīng)用于均勻矩形陣列（URA,uniform rectangular array）、L 型陣列、圓形陣列等傳統(tǒng)的二維陣列結(jié)構(gòu)。與絕大多數(shù)傳統(tǒng)的一維陣列相似，受相鄰陣元間距不超過載波半波長的限制，這些二維陣列的孔徑也較小。與之形成鮮明對比的是，近年來提出的二維互質(zhì)面陣（CPPA,co-prime planar array）結(jié)構(gòu)[13]具有大陣列孔徑的優(yōu)良特性，受到眾多學(xué)者的關(guān)注[13-16]。互質(zhì)面陣由2 個稀疏的URA 組成，2 個URA 中相鄰陣元間距為載波半波長的互質(zhì)整數(shù)倍，且均大于載波半波長[13]。在該陣列結(jié)構(gòu)下，文獻(xiàn)[13]提出了一種基于二維MUSIC 的陣列信號處理方法，用于信號二維DoA估計。文獻(xiàn)[14]在文獻(xiàn)[13]的基礎(chǔ)上，利用降維轉(zhuǎn)換將頻譜函數(shù)的維度從二維降低到一維，從而減少了算法的計算復(fù)雜度。文獻(xiàn)[15]針對互質(zhì)面陣引入二維PM，使算法的計算復(fù)雜度得到了進(jìn)一步的改善，但在低信噪比或少快拍數(shù)的情況下其信號DoA 估計性能不佳。上述文獻(xiàn)報道的二維互質(zhì)面陣及其相應(yīng)處理方法的優(yōu)勢在于其可通過增加相鄰陣元間距以獲取大陣列孔徑，從而提高信號處理效果，但不足之處在于需要同時借助2 個URA 的獨立工作來消除DoA 估計值中出現(xiàn)的相位模糊。即2 個URA 相互獨立地采用文獻(xiàn)中所提出的方法對相同入射信號進(jìn)行二維DoA 估計。這是因為，2 個URA 的相鄰陣元間距均大于載波半波長，導(dǎo)致它們的DoA 譜估計中均存在“偽峰”，為了去除“偽峰”，就需要借助2 個URA 中相鄰陣元間距為互質(zhì)整數(shù)倍的關(guān)系，對這2 個URA 的估計結(jié)果進(jìn)行比對處理。因此，這些方法雖然具有大陣列孔徑的優(yōu)點，但同時犧牲了陣列一半以上的自由度，即其所能識別處理的入射信號數(shù)小于互質(zhì)面陣陣元總數(shù)的一半。雖然文獻(xiàn)[16]對互質(zhì)面陣結(jié)構(gòu)進(jìn)行了進(jìn)一步優(yōu)化，提出了廣義互質(zhì)面陣結(jié)構(gòu)以保證其2 個URA 的陣元數(shù)相等，但受到陣列信號處理方法的限制，其自由度也只是被提升至陣元總數(shù)的一半。此外，一維互質(zhì)陣列及其處理方法主要用于提高陣列自由度，而二維互質(zhì)面陣及其處理方法主要用于提高陣列孔徑但犧牲了陣列自由度。

為了在保留大陣列孔徑優(yōu)勢的前提下提高二維互質(zhì)面陣的自由度，本文借助張量代數(shù)提出了一種新的陣列信號處理方法。首先將互質(zhì)面陣中所包含的2 個URA 均劃分成若干個重疊子面陣，并將這些子面陣的接收信號組合成一個張量；然后對這2 個URA 所對應(yīng)的2 個接收信號張量進(jìn)行互相關(guān)操作，并將結(jié)果處理成一個虛擬陣列的接收信號張量。分析結(jié)果表明，利用所提方法可將一個具有2L2? 1個物理陣元的互質(zhì)面陣轉(zhuǎn)換成一個具有個陣元的虛擬稀疏非均勻面陣，從而大幅提高了陣列自由度。為避免使用二維譜峰搜索，本文給出了采用張量分解從虛擬面陣信號張量中獲取信號二維DoA 的方法。與文獻(xiàn)[13-16]報道的互質(zhì)面陣信號處理方法相比，所提方法將陣列自由度從L2提升至，并具有更好的信號波達(dá)角估計性能及較低的計算復(fù)雜度，仿真結(jié)果證明了所提方法的有效性。

2 張量基礎(chǔ)

張量代數(shù)，也稱為多線性代數(shù)，是經(jīng)典線性代數(shù)（矩陣代數(shù)）的自然拓展，刻畫的是多變量之間線性關(guān)系的數(shù)學(xué)理論[17]。如果把矩陣視為一個具有2 個指針?biāo)饕?r,c)的二維（二階）數(shù)組，其中r和c分別指向矩陣的行和列，那么張量就是一個具有3個或3 個以上指針?biāo)饕?i,j,k,…)的多維（高階）數(shù)組。有3 個索引的三維數(shù)組稱為三階張量。如果數(shù)組索引有N個，那么這個（超體積）數(shù)組稱為N階張量。對于一個高維數(shù)據(jù)，如果仍以矩陣化的形式對其進(jìn)行表示、分析和處理，就會丟失甚至破壞高維數(shù)據(jù)中可能存在的結(jié)構(gòu)信息。與之形成鮮明對比的是，張量的高維特性賦予了其在表示高維數(shù)據(jù)時的自然性和結(jié)構(gòu)上的緊湊性，若配合張量域的高維運算及張量分解對高維數(shù)據(jù)進(jìn)行分析和處理，則可有效發(fā)掘和利用高維數(shù)據(jù)中存在的結(jié)構(gòu)信息，從而提高高維陣列信號處理性能[17-18]。因此本文將采用張量代數(shù)理論對互質(zhì)面陣信號進(jìn)行分析和處理，以期提高信號處理效果。

下面簡述本文使用的張量代數(shù)運算算子[17-20]。

其中，?表示Kronecker 乘積運算。

3 陣列信號張量處理方法

3.1 陣列結(jié)構(gòu)

圖1 二維互質(zhì)面陣結(jié)構(gòu)（L=3）

在圖1 所示的互質(zhì)面陣結(jié)構(gòu)下，文獻(xiàn)[13-16]提出了不同的陣列信號處理方法，但這些方法均需令2 個稀疏URA 獨立工作，即2 個稀疏URA 相互獨立地對相同入射信號進(jìn)行DoA 估計。這是因為2個稀疏URA 的陣元間距均大于載波半波長，因此其DoA 估計值中都存在角度模糊，只有通過比對這2 個稀疏URA對相同入射信號的DoA估計值（含模糊角度），再借助2 個稀疏URA 陣元間距互質(zhì)的關(guān)系，才可消除角度模糊。這就意味著，這些處理方法可識別的入射信號數(shù)不到陣元總數(shù)的一半。因此，雖然二維互質(zhì)陣列的大孔徑帶來了DoA 估計性能的提升，但是其現(xiàn)有的信號處理方式卻犧牲了一半的陣列自由度。

為了提高互質(zhì)面陣的自由度并進(jìn)一步提升其信號參數(shù)估計性能，接下來將借助張量代數(shù)對其接收信號進(jìn)行全新的張量建模和處理。

3.2 陣列信號的張量模型

首先，如圖2 所示，將互質(zhì)面陣中的一個稀疏URA 劃分成N x×Ny個大小均為L x×Ly的子面陣，且L=N x+Lx? 1=N y+Ly? 1，并以第(1,1)個子面陣作為參考子面陣。

圖2 稀疏URA 中的參考子面陣及x 和y 方向上的子面陣

至此，本文已經(jīng)將互質(zhì)面陣中的2 個稀疏URA的接收信號數(shù)據(jù)處理成了2 個具有明確物理含義的四階張量。接下來，對其進(jìn)行求互相關(guān)等操作，即利用其高維二階統(tǒng)計量來提高陣列自由度。

3.3 虛擬陣列的生成

受文獻(xiàn)[5-6]的啟發(fā)，本節(jié)將利用張量代數(shù)運算對互質(zhì)面陣信號的“高維二階統(tǒng)計量”進(jìn)行處理，以期獲得一個具有高自由度及大孔徑的虛擬陣列的接收信號張量。對張量 X1和X2求互相關(guān)張量得

令JD為D×D的互換矩陣，其反對角線上元素均為1，其他元素均為0。又因為ULA 存在共軛倒序不變性[20-23]，即

因此，對張量R 的共軛R?的每一維度的元素指針順序進(jìn)行取反操作，可以獲得一個新張量，即

圖3 虛擬稀疏非均勻面陣結(jié)構(gòu)（L=4，L x=Ly=3，Nx=Ny=2，M1=5，M2=4）

利用上述提出的方法對互質(zhì)面陣信號進(jìn)行處理后，原陣列（如圖1 所示）接收信號張量1X 和X2轉(zhuǎn)換成為一個虛擬的稀疏非均勻面陣（如圖3 所示）接收信號張量。接下來，對該虛擬面陣的可辨識性（即可處理的入射信號數(shù)或自由度）進(jìn)行分析，并給出如何利用張量分解從中估計出入射信號二維DoA 的方法。

3.4 虛擬陣列的可辨識能力分析

3.5 基于張量分解的信號DoA 估計

其中，Ψ1y=diag{Φ11,Φ12,…,Φ1K}，Φ1k=ejπM1vk為第k個入射信號入射至第一個稀疏URA 時，其在y方向上相鄰陣元之間的相位差。

同理，可以獲得Ψ2x=diag{Θ21,Θ22,…,Θ2K}，以及Ψ2y=diag{Φ21,Φ22,…,Φ2K}，即得到K個入射信號入射至第二個稀疏URA 時，其在x和y方向上相鄰陣元的相位差。

文獻(xiàn)[13]指出，當(dāng)陣元間距大于載波半波長時，陣列孔徑較大，接收信號估計性能較好，但會出現(xiàn)角度模糊問題。而二維互質(zhì)陣列中的2 個稀疏URA 的相鄰陣元間距均大于載波半波長，因此從1xΨ和Ψ1y中估計的入射信號DoA 值均存在角度模糊，從Ψ2x和Ψ2y估計的入射信號DoA 值也都存在角度模糊。

但文獻(xiàn)[13]同時還指出，在使用二維互質(zhì)陣列進(jìn)行DoA 估計時，其包含的2 個URA 陣元間距之間的互質(zhì)特性可用來消除角度模糊。因此，也可以利用該特性消除Ψ1x、Ψ1y、Ψ2x、Ψ2y中存在的角度模糊，具體方法如下。

首先，分別取出Ψ1x、Ψ1y、Ψ2x、Ψ2y對角元素的相位α1k、β1k、α2k、β2k，且表示為

其中，k=1,…,K，α1k、β1k、α2k、β2k均為含有角度模糊的相位觀測值，且均在[?π,π]范圍內(nèi)。

特別地，當(dāng)M1=M2=1時，即相鄰陣元間距為載波半波長時，k1x=k1y=k2x=k2y=0，即不包含角度模糊，此時入射信號的俯仰角和方位角可通過式(22)求得。

當(dāng)M1=M2＞ 1時，方位角和俯仰角估計值可能出現(xiàn)“模糊”現(xiàn)象[13]，即k1x、k1y、k2x、k2y的值不為零，但此時文獻(xiàn)[13-16]中的解模糊方法也適用于式(21)，即

4 仿真結(jié)果分析

將所提方法與文獻(xiàn)[13-16]中的方法進(jìn)行仿真對比研究，且將互質(zhì)陣列的克拉美勞界（CRB,Cramér-Rao bound）[14]作為參考。需要注意的是，文獻(xiàn)[13]和文獻(xiàn)[16]的區(qū)別在于陣列結(jié)構(gòu)，而它們所使用的陣列信號處理方法是相同的，且本文采用的是文獻(xiàn)[16]的陣列結(jié)構(gòu)。

入射信號DoA 均方誤差（RMSE,root mean square error）表示為

4.1 所提方法陣列自由度驗證

圖4 所提方法的DoA 估計結(jié)果

從圖4 中可以看出，與理論分析一致，本文所提的具有37 個自由度的方法可以成功辨識出這36個信號（物理陣元僅為31）。

在該仿真條件下，互質(zhì)面陣共包含2 個大小均為4 ×4 的稀疏URA。由于文獻(xiàn)[13-16]中方法均需令2 個稀疏URA 獨立工作。因此，在該陣列結(jié)構(gòu)參數(shù)情況下，文獻(xiàn)[13-16]中方法無法完成這36 個入射信號的DoA 估計。具體原因分析如下。

文獻(xiàn)[13-14,16]中均采用MUSIC 類方法，在該仿真條件下，其每個URA 大小均為4 ×4，共16 個陣元。按照空間譜估計理論可知，每個URA接收信號協(xié)方差矩陣大小為16 ×16 。對其進(jìn)行特征值分解可得到16 個特征值，因此其最多只能處理15 個入射信號。相似地，文獻(xiàn)[15]采用FM 對互質(zhì)面陣中的2 個稀疏URA 接收信號分別進(jìn)行DoA 估計，因此需要對每個稀疏URA（大小為4 ×4）陣列流型矩陣進(jìn)行分塊處理。而該陣列流型矩陣大小為16 ×36 （對應(yīng)16 個陣元和36 個入射信號），但FM 要求陣列流型矩陣的行向量數(shù)多于列向量數(shù)，故此時算法失效，即該算法只能處理15 個入射信號。

可見，在自由度提升方面，所提方法顯著優(yōu)于文獻(xiàn)[13-16]算法。

4.2 信號DoA 估計性能對比

對各方法的DoA 估計RMSE 性能隨信噪比及快拍數(shù)變化情況進(jìn)行仿真實驗對比。假設(shè)有2 個遠(yuǎn)場窄帶不相關(guān)信號入射至陣列，其DoA 分別為(θ1,φ1)=(28?,31?)，(θ2,φ2)=(58?,64?)，固定快拍數(shù)為200，蒙特卡洛數(shù)為500，信噪比從?10 dB變化到20 dB，其他條件與4.1 節(jié)相同。各方法的RMSE 如圖5 所示?？梢钥闯觯岱椒ň哂凶罴训腞MSE 性能，且隨信噪比的變化較平穩(wěn)。此外，其他方法在較低信噪比下的RMSE 性能并不理想，其中文獻(xiàn)[15]采用FM，性能最差，這是因為低信噪比可能造成解模糊失?。磦畏迮c偽峰的距離過近，被誤認(rèn)為是2 個真實的譜峰[13]），而所提方法由于具有較高的陣列自由度因此表現(xiàn)較優(yōu)。

圖5 RMSE 隨信噪比變化的結(jié)果

固定SNR 為5 dB，對比各方法RMSE 性能隨快怕數(shù)變化的情況?？炫聰?shù)從100 變化到800，則各方法的RMSE 性能如圖6 所示。由圖6 可知，所提方法同樣具有最優(yōu)的RMSE 性能。

圖6 RMSE 隨快怕數(shù)變化的結(jié)果

由實驗結(jié)果可知，和文獻(xiàn)[13-16]的方法相比，所提方法具有最優(yōu)的信號DoA 估計性能，且在信噪比較低或快拍數(shù)較少的情況下，其優(yōu)勢更加明顯。

4.3 計算復(fù)雜度對比

在陣列物理陣元總數(shù)變化的情況下，利用各方法的平均執(zhí)行時間對它們的計算復(fù)雜度進(jìn)行仿真對比。仿真的硬件條件為Intel(R) Core(TM)2 Quad CPU Q9550 @2.83 GHz，6 GB RAM。仿真軟件為MATLAB2016a。其他仿真條件設(shè)置如下?？炫臄?shù)為500，信噪比為?5 dB，文獻(xiàn)[13-14,16]中網(wǎng)格搜索步長均設(shè)置為0.000 2，蒙特卡洛數(shù)為100；互質(zhì)面陣結(jié)構(gòu)設(shè)置為M1=6，M2=5，L從3 變化到7，對應(yīng)的陣列物理陣元總數(shù)分別為17、31、49、71、97。假設(shè)有2 個遠(yuǎn)場窄帶不相關(guān)入射信號，其DoA分別為(θ1,φ1)=(20?,30?)，(θ2,φ2)=(50?,60?)。各方法平均每次運行所需的時間隨物理陣元數(shù)變化結(jié)果如圖7 所示，各方法RMSE 隨物理陣元數(shù)變化的結(jié)果如圖8 所示。

由圖7 可知，由于所提方法不需要二維譜峰搜索，因此表現(xiàn)出僅次于文獻(xiàn)[15]方法的良好的系統(tǒng)實時性，但需注意的是文獻(xiàn)[15]方法在低信噪比或少快拍數(shù)的情況下信號估計性能較差（如前實驗結(jié)果所示）。此外，由圖8 可知，在不同物理陣元數(shù)的情況下，所提方法同樣具有最好的RMSE 性能。

圖7 運行時間隨陣元總數(shù)變化的結(jié)果

圖8 RMSE 隨物理陣元總數(shù)變化的結(jié)果

5 結(jié)束語

本文提出了一種具有高自由度的互質(zhì)面陣信號的張量處理方法。該方法可將一個由2L2?1 個物理陣元構(gòu)成的互質(zhì)面陣轉(zhuǎn)換成一個由個陣元構(gòu)成的虛擬稀疏非均勻面陣，從而大幅提高了陣列自由度。此外，本文還給出了利用張量分解從該虛擬陣列接收信號張量中獲取信號二維DoA的方法，從而避免了二維譜峰搜索，降低了計算復(fù)雜度。理論分析和仿真實驗驗證了所提方法的有效性。