基于L型和差嵌套陣的二維波達(dá)方向估計(jì)

朱大君 何培宇 喻偉闖

（四川大學(xué)電子信息學(xué)院，四川成都 610065）

1 引言

在陣列信號(hào)處理領(lǐng)域，波達(dá)方向估計(jì)（Direction of Arrival，DOA）估計(jì)是至關(guān)重要的部分，在例如雷達(dá)，射電天文學(xué)和無線通信［1-4］等領(lǐng)域有著廣泛的研究與應(yīng)用。傳統(tǒng)的DOA 估計(jì)研究大多集中在均勻線性陣列（Uniform Linear Array，ULA）上，由于ULA的自由度受陣元數(shù)限制，導(dǎo)致具有N個(gè)傳感器的ULA，最多可檢測N-1個(gè)信源。

為了提高檢測能力，有學(xué)者提出稀疏陣列。2010 年Vaidyanathan 和Pal 提出了兩種稀疏陣列，嵌套陣［5］（Nested Array，NA）和互質(zhì)陣［6］（Coprime Array，CA），通過向量化接收到的信號(hào)協(xié)方差矩陣，使用N1+N2個(gè)陣元的NA 就可以達(dá)到2N2(N1+1) -1 的自由度（Degrees of Freedom，DOF）；2N1+N2-1 個(gè)陣元的CA 可以達(dá)到2N1(N2+1) -1 的自由度，遠(yuǎn)大于ULA 的N1+N2的自由度。但是陣元位置間的差分陣（difference co-array，DCA）存在孔洞，導(dǎo)致陣列的DCA 中央連續(xù)部分，才能應(yīng)用空間平滑MUSIC［7-9］（Spatial Smoothing Multiple Signal Classification，SS-MUSIC）算法來估計(jì)入射信號(hào)的DOA。近期，有學(xué)者提出了一些改進(jìn)版本的稀疏陣列，如超級(jí)嵌套陣［10］（Super Nested Array，SNA）、增強(qiáng)嵌套陣［11］（Augmented Nested Array，ANA）和廣義嵌套陣［12］（Generalized Nested Array，GNA）等，來提高陣列自由度。

但這些稀疏陣列僅考慮差分陣列DCA，未考慮求和陣列［13-14］（sum co-array，SCA）。適當(dāng)?shù)乩藐嚵械腟CA，可以填補(bǔ)DCA 中的一些孔洞，得到比僅考慮DCA 的稀疏陣列獲得更高的DOF。2017 年Xinghua Wang 等人提出VCAM［15］（Vectorized Conjugate Augmented MUSIC，VCAM）算法，對物理陣元同時(shí)進(jìn)行DCA 和SCA 操作，稱為DSCA（difference and sum co-array，DSCA）。但是這些稀疏陣列，其DSCA的虛擬陣元間仍然存在很多冗余陣元和孔洞，導(dǎo)致陣列自由度擴(kuò)展有限。為了更好的解決這些問題，本文提出了一種基于求和求差的改進(jìn)嵌套陣，記為INA（Improved Nested Array，INA），使用N個(gè)陣元，就可以實(shí)現(xiàn)O(N2)的自由度。

另外，在DOA 估計(jì)應(yīng)用時(shí)，不僅需要知道入射信號(hào)的方位角，還需要知道俯仰角，即需要進(jìn)行二維的DOA 估計(jì)。在［16］中已經(jīng)證明，L 型陣列不僅有簡單陣列布置結(jié)構(gòu)，由兩個(gè)互相垂直的子陣組成，而且具有更好的估計(jì)性能。但是傳統(tǒng)的L 型陣列的每個(gè)子陣都是ULA，具有低分辨率和物理陣元限制的問題。此外，二維（Two-Dimensional，2D）DOA 估計(jì)時(shí)需要使用2D-MUSIC算法的進(jìn)行二維的角度搜索，其計(jì)算復(fù)雜度高。

因此本文將提出的INA 擴(kuò)展到L 陣，進(jìn)行二維的DOA 估計(jì)。此L 型陣列的兩個(gè)子陣都是INA，通過VCAM 算法對陣元位置之間的差分、求和操作達(dá)到虛擬擴(kuò)展陣元數(shù)目的效果，虛擬陣元數(shù)目具有確定的表達(dá)式。此L 陣能夠提升陣列天線的自由度，更加準(zhǔn)確的估計(jì)出更多的信源波達(dá)方向，突破了物理陣元的限制。為了避免二維角度搜索［17-18］，采用PSCM［19-20］（Pair-matching Signal Covariance Matrices，PSCM）算法，對L 陣兩個(gè)子陣的一維DOA 估計(jì)角度，進(jìn)行角度的匹配。文中第3.2 節(jié)給出自由度的數(shù)學(xué)表達(dá)式和第4 節(jié)實(shí)驗(yàn)仿真來說明此陣列的優(yōu)異性。

本文符號(hào)說明，(.)T、(.)*和(.)H分別表示矩陣的轉(zhuǎn)置、共軛和共軛轉(zhuǎn)置。vec(.)代表向量化操作，⊙和?分別表示Khatri-Rao積和Kronecker積。

2 陣列信號(hào)模型

2.1 INA模型

常規(guī)的均勻線陣（ULA）進(jìn)行DOA 估計(jì)時(shí)需要滿足空間采樣定理，即陣元間間隔不能大于接收信號(hào)的半波長，這限制了陣列的分辨率，對于N個(gè)陣元的ULA最多估計(jì)N-1個(gè)信號(hào)。為了解決上述問題，獲得更高的自由度和分辨率，本文提出了一種基于求和求差的改進(jìn)嵌套陣，記為INA。陣列結(jié)構(gòu)如圖1 所示，陣列由子陣1 和子陣2 兩部分組成，共有N個(gè)陣元，N=N1+N2。子陣1 是一個(gè)密集的均勻線陣，共有N1個(gè)陣元，陣元間距為d=λ/2，λ是入射信號(hào)波長；子陣2 是一個(gè)稀疏線陣，共有N2個(gè)陣元，子陣2的0位置陣元和子陣1的N1-1位置陣元的間距也為d，子陣2前N2-1個(gè)陣元之間的間距為(2N1+2)d，最后兩個(gè)陣元之間的間距為(2N1+1)d。

根據(jù)圖1 所示的INA 陣元位置結(jié)構(gòu)，可以定義SINAd為陣元位置，其中SINA是陣元位置集

陣列孔徑LINA為

2.2 L型INA信號(hào)模型

如圖2 所示的L 型INA，該陣列由位于x軸和z軸上的兩個(gè)子陣組成，每一個(gè)子陣都是由2.1 提出的INA 組成，原點(diǎn)被兩個(gè)子陣共用，總的陣元個(gè)數(shù)為2N-1，N=N1+N2。位于x軸和z軸上陣元的位置可以表示為(xid，0)和(0，zid)，其中xi，zi∈SINA，i=1，2，…，N。

假定有K個(gè)遠(yuǎn)場不相關(guān)的窄帶信號(hào)s(t)以角度(θk，φk)，k=1，2，…，K入射到圖2 所示的陣列，其中θk∈(0，π)是信號(hào)s(t)與x軸的夾角，稱為方位角；φk∈(0，π)是信號(hào)s(t)與z軸的夾角，稱為俯仰角。則對于x軸和z軸上第t個(gè)快拍的陣列輸出可建模為：

其中s(t)=[s1(t)，s2(t)，…，sK(t)]T∈CK×1是信號(hào)源矢量，nx(t)=[nx1(t)，nx2(t)，…，nxN(t)]T∈CN×1和nz(t)=[nz1(t)，nz2(t)，…，nzN(t)]T∈CN×1分別是x軸和z軸上的噪聲矢量，彼此是相互獨(dú)立的均值為零方差為的高斯白噪聲，并且與信號(hào)源不相關(guān)。Ax(θ)和Az(φ)是x軸和z軸上的陣列流形，x(t)和z(t)是x軸和z軸上第t個(gè)快拍的接收信號(hào)，可以表示為

ax(θk)和az(θk)，k=1，2，…，K分別是第k個(gè)信號(hào)在x軸和z軸上的導(dǎo)向矢量，表示為

3 陣列孔徑虛擬擴(kuò)展

3.1 VCAM算法

利用VCAM 算法可以生成2.1 節(jié)所提出的INA的虛擬擴(kuò)展陣列，擴(kuò)展后的虛擬陣列為均勻線陣，虛擬陣元之間的間距為d=λ/2。由N1+N2個(gè)物理陣元可虛擬出4N1N2+4N2-3個(gè)連續(xù)的虛擬陣元，即陣列的自由度由N1+N2提升到了4N1N2+4N2-3。下面介紹VCAM算法。

對于2.2 節(jié)中所提出的x軸上的陣列，在位置(xmd，0)和(xnd，0)，m，n=1，2，…，N上陣元的輸出可以表示為xm(t)和xn(t)。則xm(t)和xn(t)的時(shí)間自相關(guān)函數(shù)可以表示為

當(dāng)xn=0 時(shí)，即取(0，0)位置陣元作為參考陣元，向量化可以得

聯(lián)合式子（12）和式子（13）可以得

用將Ts和Ns表示為偽采樣周期和偽快照的數(shù)量，r(τ)的偽數(shù)據(jù)矩陣可以表示為

R的協(xié)方差矩陣表示為

同理對于z軸上的和差嵌套陣也可以通過VCAM算法生成其虛擬擴(kuò)展陣列。

3.2 INA的虛擬擴(kuò)展

利用VCAM 算法可以生成INA 的DSCA，陣元間 距d=λ/2，由N1+N2個(gè)物理陣元可虛擬出4N1N2+4N2-3 個(gè)連續(xù)的虛擬陣元，虛擬陣元的位置集SVINA在內(nèi)連續(xù)

當(dāng)N1=N2=4 時(shí)，即一共有8 個(gè)陣元的INA 的物理陣元結(jié)構(gòu)，和其經(jīng)過DSCA 虛擬后陣元結(jié)構(gòu)，記為INA（DSCA）。如圖3 所示，連續(xù)陣元的區(qū)間為［-38，38］，即陣列的自由度由8提升到了77。

為了說明本陣列的優(yōu)越性，與同樣是8 個(gè)陣元的嵌套陣NA、互質(zhì)陣CA 和增強(qiáng)嵌套陣ANA 來進(jìn)行比較。對NA 和CA 同時(shí)進(jìn)行DCA 和DSCA 操作，對ANA 僅進(jìn)行DCA 操作，為了方便，將其相對應(yīng)的虛擬陣列分別記為NA（DCA）、NA（DSCA）、CA（DCA）、CA（DSCA）、ANA（DCA）。上述陣列的結(jié)果如圖4 所示，可以發(fā)現(xiàn)同樣是8 個(gè)陣元的情況下，下面三種陣列的最大連續(xù)范圍（圖（d））是［-28，28］，自由度為57，比本文提出的INA 陣列的自由度要少20，說明了本陣列的具有更高的自由度。

除此之外，同時(shí)對比圖（a）和圖（b）、或者圖（c）和圖（d），可以發(fā)現(xiàn)同一陣列經(jīng)過DSCA 后的虛擬陣元數(shù)是要大于只進(jìn)行DCA 的虛擬陣元數(shù)，說明了陣列進(jìn)行DSCA后能夠獲得更高的自由度。

當(dāng)NA（DSCA）、CA（DSCA）、ANA（DCA）以及INA（DSCA）的物理總陣元數(shù)N不變時(shí)，合理分配N1、N2的值可得到虛擬陣元數(shù)的最優(yōu)表達(dá)式如表1所示。注意下面只列舉了當(dāng)N為偶數(shù)時(shí)的表達(dá)式，N為奇數(shù)的表達(dá)式與N為偶數(shù)的數(shù)量級(jí)是一致的。

觀察表1可知，對于本文提出來的INA 的DSCA的自由度是O(N2)級(jí)別的，而NA 和CA 的DSCA 自由度僅為O(N2/2)級(jí)別，雖然ANA 的DCA 自由度為O(2N2/3)，但是還是小于O(N2)。表1 在數(shù)學(xué)表達(dá)式上說明了本陣INA具有更高的自由度。

表1 N陣元不同陣列的最大自由度Tab.1 The maximum DOF of different arrays of N elements

3.3 基于L型INA的二維DOA估計(jì)

通過VCAM 可以生成INA 的和差嵌套陣，其計(jì)算復(fù)雜度約為O(T(4N1N2+4N2-3)4)，其中T代表快拍數(shù)，N1、N2分別代表子陣1、2 的陣元數(shù)，式子（17）是向量化接收信號(hào)的協(xié)方差矩陣，矢量z等價(jià)于經(jīng)過DSCA 虛擬化后的和差嵌套陣的接收信號(hào)。但是由于式子中矩陣p的秩為1，無法直接用MUSIC算法進(jìn)行正確的DOA 估計(jì)，需要先進(jìn)行秩的恢復(fù)，可以采用文獻(xiàn)［9］提出的Toeplitz 矩陣重構(gòu)的方法來恢復(fù)矩陣的秩。

通過對物理陣元位置的DSCA 操作，會(huì)得到一些冗余陣元，矢量z中的每個(gè)元素都對應(yīng)一個(gè)虛擬陣元，因此z中也會(huì)存在冗余元素，對z中的元素進(jìn)行去冗余，并按陣元位置從小到大排序，取中間連續(xù)的部分為得下式。

采用Toeplitz 矩陣重構(gòu)的方法會(huì)損失一半的自由度，因此基于和差嵌套陣的Toeplitz 矩陣重構(gòu)法最多可以估計(jì)-1個(gè)信號(hào)的入射角度。

基于L型INA的二維DOA估計(jì)步驟總結(jié)如下：

4 實(shí)驗(yàn)仿真

實(shí)驗(yàn)1 將對比INA（DSCA）與NA（DSCA）、CA（DSCA）、ANA（DCA）分別擴(kuò)展至L 型陣上的二維DOA估計(jì)效果，L陣x軸和z軸的兩個(gè)子陣相同，其物理陣元排列方式同3.2節(jié)，每個(gè)子陣的陣元設(shè)為N=8，則L陣總陣元數(shù)為2N-1=15，四種陣列采用Toeplitz矩陣重構(gòu)法時(shí)的理論自由度分別為38、23、28、22。先對L陣的兩個(gè)子陣分別采用MUSIC算法進(jìn)行一維的DOA 估計(jì)，然后再采用PSCM 算法進(jìn)行一維角度配對，進(jìn)而完成L陣的二維DOA估計(jì)。

進(jìn)行蒙特卡洛實(shí)驗(yàn)，引入均方根誤差（RMSE）來評(píng)估上述四種陣的DOA 估計(jì)性能。定義RMSE為

其中W代表進(jìn)行獨(dú)立重復(fù)實(shí)驗(yàn)的次數(shù)，K代表總的入射信號(hào)個(gè)數(shù)，而分別代表第w次獨(dú)立實(shí)驗(yàn)中對于第k信號(hào)的方位角θk和俯仰角φk的估計(jì)值。

設(shè)有K=3 個(gè)遠(yuǎn)場不相關(guān)的窄帶信號(hào)其入射方向?yàn)?45°，70°)、(80°，40°)和(130°，90°)。

a）研究均方根誤差與信噪比的關(guān)系，信噪比SNR 在［-5 dB，15 dB］內(nèi)均勻變化，快拍數(shù)為500，譜峰搜索步長為0.01°，進(jìn)行300 次蒙特卡洛實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果如圖5 所示，隨著SNR 增大，這四種陣列的DOA 估計(jì)性能有明顯的提高，且INA（DSCA）相比其他三個(gè)虛擬陣列一直有著更低RMSE 這是因?yàn)镮NA（DSCA）具有更多的連續(xù)陣元，更高的自由度。

b）研究均方根誤差與快拍數(shù)的關(guān)系，快拍在[100，1000]內(nèi)均勻變化，信噪比SNR=5 dB，譜峰搜索步長為0.01°，進(jìn)行300次蒙特卡洛實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果如圖6 所示，隨著快拍數(shù)增加，這四種陣列的DOA估計(jì)性能有明顯的提高，且INA（DSCA）相比其他三個(gè)虛擬陣列一直有著更低RMSE，這是因?yàn)镮NA（DSCA）具有更高的自由度。

實(shí)驗(yàn)2 研究INA（DSCA）的L陣可估計(jì)信源數(shù)，并與NA（DSCA）、CA（DSCA）和ANA（DCA）做對比，其他L 陣的結(jié)構(gòu)同實(shí)驗(yàn)1。設(shè)有11個(gè)遠(yuǎn)場不相關(guān)的窄帶信號(hào)其入射方向?yàn)椋?0°，75°）、（45°，80°）（50°，140°）、（60°，55°）、（65°，85°）、（70°，120°）、（80°，60°）、（90°，90°）、（120°，95°）、（135°，35°）和（140°，50°），信噪比SNR=10 dB，快拍數(shù)為500，譜峰搜索步長為0.01°。進(jìn)行DOA 估計(jì)的結(jié)果如圖7 所示，結(jié)果表明每個(gè)子陣是8陣元INA（DSCA）組成的L陣可以準(zhǔn)確估計(jì)出11個(gè)入射信號(hào)的二維角度，突破了物理陣元個(gè)數(shù)的限制。且對比NA（DSCA）、CA（DSCA）和ANA（DCA）有更高的估計(jì)的估計(jì)精度，這是因?yàn)镮NA（DSCA）具有更高的自由度。8陣元本陣列理論自由度38，這里只估計(jì)11 個(gè)入射信號(hào)的原因是，PSCM 這類的匹配算法，盡管是針對稀疏陣列的角度配對，可以減少計(jì)算復(fù)雜度，在信源數(shù)不多時(shí)可以實(shí)現(xiàn)角度的準(zhǔn)確配對，但是當(dāng)信源數(shù)遠(yuǎn)遠(yuǎn)大于陣元數(shù)會(huì)出現(xiàn)匹配失敗的現(xiàn)象。

5 結(jié)論

本文針對傳統(tǒng)L型均勻陣列二維波達(dá)方向估計(jì)中可估計(jì)信源數(shù)目受限于陣元數(shù)、分辨率低等問題，提出了一種L 型INA 的陣列結(jié)構(gòu)。INA 由一個(gè)密集線陣和一個(gè)稀疏線陣兩部分組成，通過VCAM算法對INA 的物理陣元位置進(jìn)行求和求差（DSCA）的虛擬擴(kuò)展操作，由N1+N2個(gè)物理陣元可得到4N1N2+4N2-3 的自由度，相比NA（DSCA）、CA（DSCA）和ANA（DCA）具有更高的自由度。并與這三種陣列進(jìn)行二維DOA 估計(jì)的實(shí)驗(yàn)仿真對比，結(jié)果顯示此陣具有更高的估計(jì)精度和分辨率。以及利用此陣去估計(jì)超過子陣元數(shù)的信源數(shù)，實(shí)驗(yàn)表明此陣具有高估計(jì)精度和能突破物理陣元限制的結(jié)論。此外，文中采用的VSCM算法，是將二階統(tǒng)計(jì)量轉(zhuǎn)換為四階統(tǒng)計(jì)量來進(jìn)行運(yùn)算，即陣列自由度和估計(jì)精度提高是以增加一定的計(jì)算復(fù)雜度為代價(jià)的。