基于廣義互質(zhì)雙平行陣列的二維DOA估計方法

王 宏 何培宇 喻偉闖 崔 敖 徐自勵

（1.四川大學(xué)電子信息學(xué)院，四川成都 610065；2.中國民用航空局第二研究所，四川成都 610041）

1 引言

在陣列信號處理領(lǐng)域中，波達角（Direction of Arrival，DOA）估計一直是一個重要的研究方向，在無線通信、語音信號、聲納、電子對抗等許多領(lǐng)域中有著廣泛的研究與應(yīng)用［1-2］。估計DOA時，常需要對信號進行二維的DOA估計［3］，即需要獲取信號入射的俯仰角和方位角。為此，學(xué)者已經(jīng)提出了多種基于面陣結(jié)構(gòu)的DOA估計方法，包括矩形陣［4］、L形陣列［5］和平行線陣［6］等。通常情況下，傳統(tǒng)的平面陣列結(jié)構(gòu)需要滿足空間采樣定理，即相鄰兩陣元之間的間距不能超過入射信號波長的一半，以保證DOA估計不出現(xiàn)模糊［7］。而在當(dāng)前陣列信號處理研究中，與傳統(tǒng)陣列相比，稀疏陣列能在相同物理陣元數(shù)目下能提供更大的陣列孔徑以及更多可用的虛擬陣元，從而擁有更高的分辨力和自由度（Degrees of Freedom，DOF）［8］。

當(dāng)前，稀疏陣列的研究大多集中在一維DOA 估計問題上，其陣列結(jié)構(gòu)主要可分為最小冗余陣［9］、嵌套陣［10］和互質(zhì)陣［11］。其中，最小冗余陣自由度最大，但其陣元位置和陣列自由度無閉合表達式，不易向高維擴展［9］。嵌套陣的陣元位置和自由度雖有閉合表達式，但密集陣元間互耦影響較大［12］。相較于最小冗余陣和嵌套陣，互質(zhì)陣具有陣列結(jié)構(gòu)簡單、陣元位置和自由度有閉合表達式、陣元間互耦影響較小等優(yōu)勢，但該陣列生成的差分虛擬陣列存在“孔洞”缺失，導(dǎo)致其連續(xù)自由度較?。?3］。如何將一維稀疏陣列有效地進行二維推廣，是當(dāng)前學(xué)界仍在研究的重要問題［14］。為此，文獻［15］提出了一種基于嵌套陣的平行線陣結(jié)構(gòu)，該陣列通過以嵌套陣陣列作為平行線陣的子陣，提升了傳統(tǒng)平行線陣的自由度和分辨力，但是也導(dǎo)致其密集陣元間的互耦影響較大。文獻［16］則以互質(zhì)陣為基礎(chǔ)，提出了一種互質(zhì)平行線陣結(jié)構(gòu)，該方法雖有效地避免了嵌套陣互耦影響較大這一缺點，但是由于其互質(zhì)陣列的差分孔洞特性，即差分出的連續(xù)虛擬陣元數(shù)目較少，因此該方法的分辨力和自由度不能得到有效提升。此外，該方法使用傳統(tǒng)的二維MUSIC 估計算法對信號進行估計，計算復(fù)雜度較高。文獻［17］和文獻［18］則針對文獻［16］計算復(fù)雜度較高這一缺點，利用求根算法將二維DOA估計轉(zhuǎn)換成兩個一維DOA估計問題，有效地降低了其算法復(fù)雜度，但是其仍然未解決互質(zhì)陣的差分孔洞特性所造成的性能損失問題。

針對以上問題，本文提出了一種基于廣義互質(zhì)雙平行陣列的二維DOA估計方法。該方法利用了一種廣義互質(zhì)陣列結(jié)構(gòu)，對傳統(tǒng)平行線陣進行稀疏域的推廣，具有更高的分辨力和自由度，又較好地改進了標(biāo)準(zhǔn)互質(zhì)陣的差分虛擬陣列存在的“孔洞”缺失問題，提升了虛擬陣列的連續(xù)陣元個數(shù)。并且，本文方法是在傳統(tǒng)平行線陣估計算法的基礎(chǔ)上，提出的一種改進二維Root-MUSIC估計算法，該方法通過多項式求根，將二維DOA 估計問題轉(zhuǎn)為兩個一維的DOA 估計問題，既提高了算法的精度，又降低了復(fù)雜度。

本文所用符號說明：(.)T，(.)*，(.)H，分別表示矩陣的轉(zhuǎn)置、共軛與共軛轉(zhuǎn)置；diag(v)表示以向量v為對角元素的對角矩陣；E{.}表示求矩陣的統(tǒng)計期望，vec(.)表示將矩陣向量化；det{.}表示矩陣行列式；?表示Kronecker 積運算，I表示單位矩陣；arg(.)表示取復(fù)數(shù)矢量的相位角。

2 廣義互質(zhì)陣模型

標(biāo)準(zhǔn)互質(zhì)陣由兩個均勻線陣構(gòu)成，兩個子陣位于同一水平線上，其中均勻線陣1的陣元個數(shù)為M，陣元間距為Nd，d=λ/2；均勻線陣2 的陣元個數(shù)為N，陣元間距為Md，其中M和N互質(zhì)，且M＜N。標(biāo)準(zhǔn)互質(zhì)陣的陣列參數(shù)通?？杀硎緸?M，N)。該陣列的陣元位置集SCA為

標(biāo)準(zhǔn)互質(zhì)陣的陣列孔徑為(MN-M)d［11］。根據(jù)式(1)能得到位置SCA的差集Ld，去掉差集中的冗余元素得到，但是由于互質(zhì)陣對應(yīng)的并非完全連續(xù)，導(dǎo)致互質(zhì)陣的虛擬差分陣元位置中擁有較多的“孔洞”［11］。標(biāo)準(zhǔn)互質(zhì)陣的差集只有在[-M-N-1，M+N-1]范圍內(nèi)是連續(xù)的［7］，因此其自由度DOF為

廣義互質(zhì)陣是在標(biāo)準(zhǔn)互質(zhì)陣列結(jié)構(gòu)上發(fā)展出來的一種新的互質(zhì)陣列結(jié)構(gòu)，該陣列仍由兩個子陣組成，具有M+N-1 個陣元，M和N仍然互質(zhì)，但此時M＜N這一條件不再要求成立［19］。與標(biāo)準(zhǔn)互質(zhì)陣不同的是，廣義互質(zhì)陣引入了一個正整數(shù)的壓縮因子p來對其子陣1 的陣元間距進行壓縮，從而減小了這一均勻線陣之間的間距。兩者的陣列結(jié)構(gòu)如圖1所示。

其中，p是2到M之間的一個整數(shù)，也是一個整數(shù)。根據(jù)式（3）可以確定，和N也滿足互質(zhì)關(guān)系，且其陣列位置集SGCA為

對式（4）集合進行自差分和互差分可得陣列的差分位置集合Ld為［19］

對差分位置集合Ld中的元素分析可知，廣義互質(zhì)陣的差分位置集合Ld在[-MN++1，MN--1]范圍內(nèi)是連續(xù)的［19］。廣義互質(zhì)陣的差分優(yōu)化陣自由度DOF為

根據(jù)式（2）和式（6）可得，廣義互質(zhì)陣的虛擬差分優(yōu)化陣的自由度和陣列孔徑大于傳統(tǒng)的標(biāo)準(zhǔn)互質(zhì)陣模型，其性能隨著的減小而增大，或隨著壓縮因子p的增大而增大，但過大的壓縮因子會增加其陣元間的互耦影響［20］。

3 基于廣義互質(zhì)雙平行陣列的二維DOA估計

3.1 廣義互質(zhì)雙平行陣列模型

在陣列模型方面，本文以單個廣義互質(zhì)陣為基礎(chǔ)，構(gòu)建了一種廣義互質(zhì)平行陣列結(jié)構(gòu)。如圖2 所示，該陣列由兩個完全相同且相互平行的廣義互質(zhì)子陣組成，其中子陣1 位于y軸上，子陣2 和子陣1相距d=λ/2，λ 為入射信號的波長。在該陣列結(jié)構(gòu)中，每個子陣的陣元總數(shù)均為S=M+N-1，陣列參數(shù)均為(M，N)。從2 節(jié)可得，單個廣義互質(zhì)子陣又由兩個均勻的直線陣構(gòu)成。其中，均勻子陣ULA1 的陣元數(shù)為N，陣元間距為d1=，其中=Md/p，，p∈Z+。均勻子陣ULA 2的陣元數(shù)為M，陣元間距d2=Nd。

廣義互質(zhì)子陣1 的陣元位置集可表示為(0，zid)，其中，zi∈S，i=1，2，…，M+N-1，且集和S 的結(jié)構(gòu)如式(7)所示。類似的，廣義互質(zhì)子陣2 的陣元位置可表示為(d，zid)。

設(shè)有K個非相關(guān)的遠場窄帶信號sk(t)(k=1，2，…，K)分別從方位角θk和俯仰角φk入射到廣義互質(zhì)雙平行線陣，且用αk和βk分別表示第k個入射信號與y軸和x軸的夾角，則有以下關(guān)系成立

整個陣列的接收信號模型可以表示為

其中，s(t)=[s1(t)，s2(t)，…，sK(t)]T∈CK×1是信號源矢量，n1(t)∈C(M+N-1)×1和n2(t)∈C(M+N-1)×1分別表示為兩個與信號獨立同分布的高斯白噪聲。A1∈C(M+N-1)×K和A2∈C(M+N-1)×K分別表示為子陣1和子陣2的方向矢量矩陣，且分別為

3.2 子陣自相關(guān)矩陣和互相關(guān)矩陣的虛擬化

根據(jù)式(9)可以得到子陣1在時刻t的接收信號協(xié)方差矩陣為

其中，Λ=E{s(t)sH(t)}=表示第k個信號的入射功率，表示噪聲的功率，IM為M+N-1階單位陣。

對協(xié)方差矩陣R11向量化得到

式中，B1=，?表示Kronecker 積運算，p=為信號功率矢量，e=，ei∈R(M+N-1)×1(i=1，…，M+N-1)，矢量ei中除了第i個元素為1，其余元素均為0。

將矢量中的冗余元素進行去重，并取連續(xù)部分重排后得到

根據(jù)式(9)，再次得到廣義互質(zhì)雙平行陣之間的互相關(guān)矩陣為

同理，向量化互相關(guān)矩陣R21得到

式中，B2=[a2(αK，βK)]，將矢量z2去除冗余，并取連續(xù)重排得到，并且將依Toeplitz重排后可以得到

3.3 增廣協(xié)方差矩陣構(gòu)建及二維求根MUSIC算法

如3.2 節(jié)所述，在獲得虛擬化之后的互相關(guān)矩陣和自相關(guān)矩陣之后，按照式（18）構(gòu)建一個增廣協(xié)方差矩陣［15］

根據(jù)式（18）可以得到一個虛擬化之后的均勻雙平行線線陣接收信號的協(xié)方差矩陣，相對廣義互質(zhì)雙平行陣列模型本身，其虛擬化之后生成的平行差分優(yōu)化陣列增加了較多的虛擬陣元。對該矩陣進行特征分解以后，即可得到式（19）中的信號子空間Us和噪聲子空間Un，利用這兩個空間的正交性，即可利用傳統(tǒng)二維MUSIC 估計算法來獲得空間譜函數(shù)，其公式如下

但是，在實際的系統(tǒng)應(yīng)用中，由于傳統(tǒng)算法中二維譜搜索具有較高的計算復(fù)雜度，因此在對時延要求較低的系統(tǒng)中該方案不具備可行性［7］，可以通過改進的二維Root-MUSIC 算法最小化譜函數(shù)f(α，β)，將其分解為兩個一維的DOA 估計問題，從而對入射信號進行二維估計。

為此，首先將矩陣Un劃分成兩個子陣，如下所示

式中，ξ=[1ej2πdcos(β)/λ]T，b1(α)為平行線陣虛擬化后的子陣1連續(xù)部分所對應(yīng)的導(dǎo)向矢量，其結(jié)構(gòu)可表示為

式（23）有效地將二維的角度估計問題分解為兩個一維的角度估計問題，即Q(α)只包含角度α的信息，矢量ξ只含有角度β的信息。令v=，將式（24）重新表示為b1(v)=[1，v，…，]T，則角度αk(k=1，2，...，K)通過求解多項式方程det{Q(v)}來得到，即

其中

其中，arg(.)表示求相位，vk是多項式det{Q(v)}與角度對應(yīng)的解。接著將所求得的代入式（25）即可求得方向余弦估計值為

最終可以求得入射信號的方位角θk和俯仰角φk的估計值為

3.4 算法步驟總結(jié)

步驟1：對接收信號進行離散化采樣，在有限的快拍數(shù)下得到廣義互質(zhì)子陣1 的估計協(xié)方差矩陣；對其按照式（12）和式（13）進行向量化，去除冗余元素，在進行重排得到矢量，利用所得的矢量按式（14）進行Toeplitz重排得到。

步驟2：同步驟1，得到廣義互質(zhì)子陣1 和廣義互質(zhì)子陣2 的估計互協(xié)方差矩陣；對其按式（12）和式（13）進行向量化，去除冗余元素，在進行重排得到矢量，利用所得的矢量按式（17）進行Toeplitz重排得到。

步驟3：根據(jù)式（18）構(gòu)建增廣協(xié)方差矩陣，對進行矩陣分解得到其噪聲子空間Un，然后將噪聲子空間Un按式（22）劃分成兩個子矩陣Un1和Un2，并按式（25）和式（26）計算得到多項式方程的det{Q(v)}的解。

步驟4：將步驟3中所得到的解代入式（26）和式（27）得到估計的，再將其估計結(jié)果代入式（29），即可估計出方位角和俯仰角。

本文算法的步驟主要包括協(xié)方差矩陣R11和互協(xié)方差矩陣R21的估計，Toeplitz矩陣重構(gòu)，以及矩陣的特征分解，經(jīng)推導(dǎo)，算法的復(fù)雜度約為O{2S2T+(2L)3}，其中S=M+N-1，L=M(N+1)，T為采樣快拍數(shù)，M和N互為質(zhì)數(shù)。

4 陣列自由度分析與仿真實驗

4.1 陣列自由度分析

設(shè)各個陣列物理實際陣元數(shù)量S相等，如S=12，并設(shè)標(biāo)準(zhǔn)互質(zhì)陣和廣義互質(zhì)陣參數(shù)相同，如M=6，N=7，廣義互質(zhì)陣的壓縮因子p分別在2、3、6 的三種情況下，本文以三種雙平行線陣的一個子陣為基礎(chǔ)，繪制出了傳統(tǒng)非稀疏平行線陣、標(biāo)準(zhǔn)互質(zhì)陣和廣義互質(zhì)陣三種陣列的子陣實際陣元位置分布，從而進行三種陣列的自由度大小對比。

由圖3 可見，廣義互質(zhì)陣的自由度明顯高于傳統(tǒng)非稀疏線陣和標(biāo)準(zhǔn)互質(zhì)陣，且自由度隨著壓縮因子的增大而增大。其中，普通均勻平行線陣由于沒有虛擬差分特性，陣列自由度與實際的物理陣元個數(shù)相等。而標(biāo)準(zhǔn)互質(zhì)陣雖然利用了稀疏陣列可以生成虛擬陣元這一特性，一定程度上提升了陣列孔徑，但是由于生成的差分虛擬陣列存在較多的“孔洞”，因此可以利用的連續(xù)陣元數(shù)目不多。廣義互質(zhì)陣由于壓縮因子p的引入，使子陣1 的間距減小，從而在子陣做位置的差分時，增加了虛擬差分優(yōu)化陣列中連續(xù)陣元的數(shù)量，提高了陣列的自由度。且從圖中可以看出廣義互質(zhì)陣的壓縮因子越大，差分位置集合的連續(xù)陣元數(shù)量也就越多，當(dāng)達到最大時，虛擬差分陣列已經(jīng)完全無“孔洞”，陣列的自由度達到最大。

4.2 多信源條件下DOA估計

設(shè)有K個非相關(guān)遠場窄帶信號入射到陣列，陣列單個子陣實際物理陣元個數(shù)均為12，稀疏陣的陣列參數(shù)為M=6，N=7。信噪比SNR=20dB，快拍數(shù)T=2000，廣義互質(zhì)平行陣列的壓縮因子p=3。傳統(tǒng)雙平行線陣方法、文獻［16］方法和本文方法的估計入射信號的個數(shù)K分別為10、13 和20，在此條件下驗證三種方法的估計精度和自由度。

從圖4 中可以看出，三種方法均可以對入射信號進行較為準(zhǔn)確的二維DOA 估計。但是傳統(tǒng)雙平行線陣算法，受限于空間采樣定律，因而其陣列的孔徑較小，估計的精度相對也較低［21］。此外，傳統(tǒng)雙平行線陣算法可用陣元數(shù)等于實際的物理陣元數(shù)，因而其自由度小，可估計的信源也較少。而文獻［16］算法利用了陣列孔徑大的稀疏陣列，具有可以生成虛擬陣元這一特性，明顯提升了信號估計的精度和增加了可估計信源的個數(shù)，但由于該方法生成的虛擬差分陣列存在較多“孔洞”，以及使用了傳統(tǒng)的二維MUSIC 方法，在面對較多信源入射時，估計性能出現(xiàn)了明顯的下降。而本文方法采用的廣義互質(zhì)陣由于陣列孔徑更大，生成的連續(xù)虛擬陣元數(shù)目也更多，且利用多項式求根的方法替代了傳統(tǒng)二維MUSIC 方法對信號進行估計，因而從圖中來看，無論是從可估計信號的數(shù)量還是估計的精度，明顯都優(yōu)于前兩者。

4.3 均方根誤差

設(shè)有兩個非相關(guān)遠場窄帶信號分別從二維角度(65°，40°)和(122°，120°)入射到陣列，即K=2，快拍數(shù)固定為T=2000，蒙特卡羅實驗次數(shù)為200，陣列單個子陣實際物理陣元個數(shù)均為S=12，稀疏陣的陣列參數(shù)為M=6，N=7，壓縮因子p=2，3，6。SNR 從-5 dB 到15 dB 均勻變化，其各個陣列的均方根誤差RMSE 隨信噪比的關(guān)系如圖5。同時，固定其信噪比SNR=20 dB，快拍數(shù)從100到2000，其各個陣列的均方根誤差（RMSE）隨快拍數(shù)的關(guān)系如圖6。文中，二維均方根誤差的定義為

由圖5和圖6可知，三種算法的RMSE均隨著信噪比和快拍數(shù)的增加而減小，但本文方法優(yōu)于傳統(tǒng)非稀疏平行線陣方法和文獻［16］方法。還可看到，廣義稀疏平行陣的RMSE 隨著壓縮因子p的增加而增大，這與實驗1 和2 的結(jié)果相吻合。因而，在同等物理陣元的條件下，本文方法的估計性能優(yōu)于傳統(tǒng)非稀疏平行線陣方法和文獻［16］方法。

5 結(jié)論

本文針對標(biāo)準(zhǔn)互質(zhì)平行線陣在進行DOA 估計時，因差分虛擬陣列“孔洞”過多，導(dǎo)致其性能不佳這一缺點，提出了一種基于廣義互質(zhì)平行陣列的二維DOA 估計方法。該方法將廣義互質(zhì)陣列的特性應(yīng)用于平行陣列二維DOA估計中，即以降低一定的陣元間距為代價，生成了更多可利用的連續(xù)虛擬陣元，較好地改進了標(biāo)準(zhǔn)互質(zhì)平行陣中“孔洞”過多的缺點，提升了其估計性能。本文算法通過一種改進的Root-MUSIC 算法將一個二維估計問題分解為兩個一維的估計問題，既提升了估計的精度，又降低了算法復(fù)雜度。然而，本文方法性能的提升是以適當(dāng)增加陣元間的互耦影響為代價的，且該影響隨壓縮因子的增大而增大，如何選擇適當(dāng)?shù)膲嚎s因子，使其在一定互耦影響范圍內(nèi)，保持較高的估計性能，在未來還有待進一步的研究。