王 宏 何培宇 喻偉闖 崔 敖 徐自勵(lì)

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

1 引言

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

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

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

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

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

標(biāo)準(zhǔn)互質(zhì)陣由兩個(gè)均勻線陣構(gòu)成，兩個(gè)子陣位于同一水平線上，其中均勻線陣1的陣元個(gè)數(shù)為M，陣元間距為Nd，d=λ/2；均勻線陣2 的陣元個(gè)數(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ì)陣對(duì)應(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ā)展出來(lái)的一種新的互質(zhì)陣列結(jié)構(gòu)，該陣列仍由兩個(gè)子陣組成，具有M+N-1 個(gè)陣元，M和N仍然互質(zhì)，但此時(shí)M＜N這一條件不再要求成立［19］。與標(biāo)準(zhǔn)互質(zhì)陣不同的是，廣義互質(zhì)陣引入了一個(gè)正整數(shù)的壓縮因子p來(lái)對(duì)其子陣1 的陣元間距進(jìn)行壓縮，從而減小了這一均勻線陣之間的間距。兩者的陣列結(jié)構(gòu)如圖1所示。

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

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

對(duì)差分位置集合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的增大而增大，但過大的壓縮因子會(huì)增加其陣元間的互耦影響［20］。

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

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

在陣列模型方面，本文以單個(gè)廣義互質(zhì)陣為基礎(chǔ)，構(gòu)建了一種廣義互質(zhì)平行陣列結(jié)構(gòu)。如圖2 所示，該陣列由兩個(gè)完全相同且相互平行的廣義互質(zhì)子陣組成，其中子陣1 位于y軸上，子陣2 和子陣1相距d=λ/2，λ 為入射信號(hào)的波長(zhǎng)。在該陣列結(jié)構(gòu)中，每個(gè)子陣的陣元總數(shù)均為S=M+N-1，陣列參數(shù)均為(M，N)。從2 節(jié)可得，單個(gè)廣義互質(zhì)子陣又由兩個(gè)均勻的直線陣構(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個(gè)非相關(guān)的遠(yuǎn)場(chǎng)窄帶信號(hào)sk(t)(k=1，2，…，K)分別從方位角θk和俯仰角φk入射到廣義互質(zhì)雙平行線陣，且用αk和βk分別表示第k個(gè)入射信號(hào)與y軸和x軸的夾角，則有以下關(guān)系成立

整個(gè)陣列的接收信號(hào)模型可以表示為

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

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

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

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

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

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

將矢量中的冗余元素進(jìn)行去重，并取連續(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)建一個(gè)增廣協(xié)方差矩陣［15］

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

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

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

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

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

其中

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

最終可以求得入射信號(hào)的方位角θk和俯仰角φk的估計(jì)值為

3.4 算法步驟總結(jié)

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

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

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

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

本文算法的步驟主要包括協(xié)方差矩陣R11和互協(xié)方差矩陣R21的估計(jì)，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 陣列自由度分析與仿真實(shí)驗(yàn)

4.1 陣列自由度分析

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

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

4.2 多信源條件下DOA估計(jì)

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

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

4.3 均方根誤差

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

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

5 結(jié)論

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