互質(zhì)線陣中一種基于共軛增廣的DOA 估計算法

林新平 張小飛 沈金清

（南京航空航天大學電子信息工程學院，南京，211106）

引 言

近年來，信號波達方向(Direction of arrival,DOA)估計被廣泛應(yīng)用于雷達、無線通信、聲吶和生物醫(yī)學等領(lǐng)域，是陣列信號處理中的一個研究熱點[1-3]。經(jīng)典的超分辨率DOA估計算法，例如多重信號分類(Multiple signal classification,MUSIC)算法[4]、借助旋轉(zhuǎn)不變性進行信號參數(shù)估計(Estimation of signal parameters via rotational invariance technique,ESPRIT)算法[5]等，都是基于特征結(jié)構(gòu)的子空間類方法。子空間類方法最初是針對傳統(tǒng)滿陣提出的，要求陣元間距嚴格小于等于電磁波半波長，以避免出現(xiàn)角度估計模糊。然而，傳統(tǒng)滿陣由于陣元間距過小，存在陣元間互耦嚴重、陣列孔徑不足以及測向精度不高等缺點，已經(jīng)無法滿足許多航空航天以及軍事民用領(lǐng)域的應(yīng)用需求。

文獻[6]提出了互質(zhì)陣列的概念。它是一種陣元間距大于半波長的非均勻陣列，由兩個陣元數(shù)與陣元間距存在互質(zhì)關(guān)系的均勻子陣穿插拓撲構(gòu)成。相比于傳統(tǒng)滿陣，互質(zhì)陣列具有陣元互耦更低，陣列孔徑更大，定位測向精度更高等優(yōu)點。文獻[6]中證明了一個具有M+N-1個陣元的互質(zhì)線陣，能夠獲得O{MN}的空間自由度。因此，基于互質(zhì)陣列的空間譜估計研究成為當下信號處理領(lǐng)域的熱點問題之一。

文獻[7]中提出了一種互質(zhì)線陣下基于矢量化協(xié)方差矩陣的DOA估計方法，稱之為虛擬化方法。該方法通過矢量化接收信號協(xié)方差矩陣進行數(shù)據(jù)重構(gòu)，從而得到一個通過虛擬陣列接收到的單快拍信號。特別地，虛擬陣列能夠提供比物理陣列更大的陣列孔徑和更高的空間自由度。文獻[8]中則提出了基于互質(zhì)特性的聯(lián)合MUSIC方法，該方法利用互質(zhì)線陣的兩個子陣單獨估計DOA，通過比對子陣間的估計結(jié)果消除測向模糊，稱之為解模糊方法。相比于虛擬化方法，解模糊方法實現(xiàn)簡單，且DOA估計性能更優(yōu)[9]。然而文獻[8]中的MUSIC算法需要全局譜搜索，算法復雜度較高。文獻[10]和文獻[11]中分別提出了MUSIC部分空間譜搜索的方法和求根MUSIC算法，復雜度相對于文獻[8]中的方法具有較大改善。相對于MUSIC算法，互質(zhì)線陣下的ESPRIT算法[12]顯然在算法復雜度上要低得多。文獻[13-15]中通過對接收數(shù)據(jù)矩陣重構(gòu)和實值化的方法達到降低互質(zhì)陣列中ESPRTI算法角度估計復雜度的目的。然而，上述DOA估計方法只利用了接收信號的空域信息，并未對時域信息加以利用[16]。

為了能夠充分利用互質(zhì)線陣接收信號的時域與空域信息，本文提出了一種基于共軛增廣(Conjugate augmented)的DOA估計算法。算法首先通過利用不同時長間隔下陣元接收信號間互相關(guān)函數(shù)的共軛對稱特性，構(gòu)造一個共軛增廣虛擬陣列接收信號，虛擬接收信號具有和實際信號相同的信號形式，但相比于原有的互質(zhì)陣列，虛擬陣列中增加了近一倍的鏡像虛擬陣元，陣列孔徑更大，空間自由度更高。然后采用基于互質(zhì)特性的聯(lián)合酉ESPRIT(Unitary ESPRIT)算法[17]對虛擬信號實現(xiàn)DOA估計，稱該算法為基于共軛增廣的酉ESPRIT(CA-UESPRIT)算法。相比于傳統(tǒng)互質(zhì)線陣下的聯(lián)合UESPRIT算法，CA-UESPRIT算法自由度更高，DOA估計性能更優(yōu)。CA-UESPRIT算法通過引入酉變換矩陣降低了ESPRIT算法中的復乘運算次數(shù)，算法復雜度較低。

1 陣列結(jié)構(gòu)與數(shù)據(jù)模型

考慮如圖1所示互質(zhì)線陣，它由兩個均勻子陣重合第1個陣元拓撲而成，子陣結(jié)構(gòu)如圖2所示，其中子陣1和子陣2陣元數(shù)分別為M和N,陣元間距為d1=Nλ/2和d2=Mλ/2，λ為電磁波波長且M和N互為質(zhì)數(shù)?；ベ|(zhì)線陣總陣元數(shù)為M+N-1。

假設(shè)互質(zhì)線陣置于x軸坐標系上，以圖1虛線位置為坐標原點，那么互質(zhì)線陣陣元位置與坐標原點距離可表示為

圖1 互質(zhì)線陣Fig.1 Coprime linear array

圖2 子陣拓撲圖Fig.2 Topology of subarray

假設(shè)K個遠場窄帶獨立平面波信號以角度θ=[θ1,θ2,…,θK]入射到上述互質(zhì)陣列中，信號快拍數(shù)為J。在t時刻，假設(shè)第k個入射信號表示為sk(t)=Akejωkt，其中Ak為信號幅度，ωk表示信號頻率，k=1,2,…,K。那么互質(zhì)線陣中第p(p=1,2,…,M+N-1)個陣元接收信號可以表示為[18]

式中：np(t)表示加性噪聲，lp表示集合L的第p個元素，根據(jù)式(1)可知l1=0。

陣列接收信號用矩陣形式可以表示為[19]

式中：s(t)=[s1(t),s2(t),…,sK(t)]T為信源矢量，n(t)為均值為0的加性高斯白噪聲，協(xié)方差均為，且噪聲與信號間互不相關(guān)。為互質(zhì)陣列的方向矩陣，第k列表示為

2 基于共軛增廣的DOA估計算法

CA-UESPRIT算法實現(xiàn)DOA估計主要包含兩個步驟，先是利用不同時長間隔下陣元接收信號間互相關(guān)函數(shù)的共軛對稱特性，構(gòu)造一個共軛增廣虛擬陣列接收信號，隨后通過基于互質(zhì)特性的聯(lián)合UESPRIT算法[17]求解虛擬陣列的DOA估計結(jié)果。

2.1 共軛增廣過程

本節(jié)通過構(gòu)造共軛增廣虛擬陣列以擴展互質(zhì)線陣的陣列孔徑及提高空間自由度。根據(jù)式(2)，互質(zhì)線陣中第1個陣元接收信號與第p個陣元接收信號關(guān)于時長τ(τ≠0)的互相關(guān)函數(shù)可以表示為[16]

式中，Rsksk(τ)表示成[16]

(τ)可以表示為[16]

根據(jù)式(7)，可以將式(5)重寫為

構(gòu)造矢量Rs(τ)=[Rs1s1(τ),Rs2s2(τ),…,Rsksk(τ)]T，R(τ)=[Rx1x1(τ),Rx2x1(τ),…,RxM+N-1x1(τ)]T，可以得到

根據(jù)式(6)可知 Rs(τ)=Rs*(-τ)，因此可以得到 (R(-τ))*=A*Rs(τ)。令矩陣 A(2)和矢量 R(2)(τ)分別表示矩陣A和矢量R(τ)的第2～M+N-1行，那么R(2)(τ)=A(2)Rs(τ)，構(gòu)造為

圖3 虛擬互質(zhì)線陣Fig.3 Virtual coprime linear array

設(shè)圖3中虛線位置為坐標原點，可以看出，虛擬陣列關(guān)于原點鏡像對稱，稱之為共軛增廣虛擬陣列。其中，虛擬陣列陣元總數(shù)為2(M+N-1)-1。根據(jù)式(10)，構(gòu)造偽快拍數(shù)據(jù)矩陣

式中：Ts為偽采樣間隔，偽快拍數(shù)為偽信號矩陣。那么y可以看作是虛擬互質(zhì)陣列的虛擬接收信號值。

2.2 DOA估計

本節(jié)利用基于互質(zhì)特性的聯(lián)合UESPRIT算法求解虛擬信號y的DOA估計結(jié)果。將圖3中虛擬陣列拆分為兩個均勻子陣，如圖4所示。

圖4中子陣1和子陣2陣元數(shù)分別為2M-1和2N-1,陣元間距分別為d1=Nλ/2和d2=Mλ/2。設(shè)y1和y2分別為子陣1和子陣2的接收信號值，那么y1和y2可以表示成

式中：Ae,1=[ae,1(θ1),ae,1(θ2),…,ae,1(θK)]∈ C(2M-1)×K和 Ae,2=[ae,2(θ1),ae,2(θ2),…,ae,2(θK)]∈C(2N-1)×K分別表示子陣1和子陣2的方向矩陣，Ae,1和Ae,2第k列分別表示為

圖4 虛擬陣列子陣拓撲圖Fig.4 Subarray topology of virtual array

式中：k=1,2,…,K?？梢缘玫教摂M子陣1和子陣2的接收信號協(xié)方差矩陣分別為R1=E{y1yH1}和

虛擬子陣中陣元間距均大于半波長，因此先利用UESPRIT算法分別得到兩個子陣的模糊DOA估計結(jié)果，再利用互質(zhì)特性聯(lián)合子陣DOA估計值消除測向模糊。

2.2.1 UESPRIT算法

兩個虛擬子陣具有相似的陣列結(jié)構(gòu)，因此，以子陣1為例詳細闡述UESPRIT算法過程，然后推廣至子陣2同理得到子陣2的DOA估計結(jié)果。

對于子陣1，引入酉矩陣為

式中：ΠM-1為反對角矩陣。酉矩陣滿足對協(xié)方差矩陣R1進行酉變換

式中：Re(·)表示取實值操作。酉變換后虛擬子陣1協(xié)方差矩陣變?yōu)閷崝?shù)矩陣，方向矩陣由Ae,1變?yōu)?/p>

可知方向矩陣Ae,1滿足旋轉(zhuǎn)不變關(guān)系式

式中：φx=diag(e-j2πd1sin(θ1)/λ,e-j2πd1sin(θ2)/λ,…,I(2M-2)×(2M-2)]。J1Ae,1和J2Ae,1分別表示選取方向矩陣Ae,1的前(2M-2)行和后(2M-2)行。根據(jù)文獻[20]，酉變換后，方向矩陣由 Ae,1變?yōu)?Ω1，因此式(18)變?yōu)?/p>

對式(17)中實數(shù)協(xié)方差矩陣進行特征分解，取最大K個特征值對應(yīng)的特征矢量構(gòu)成信號子空間Es1，可知Es1與實數(shù)方向矩陣Ω1間滿足關(guān)系式Es1=Ω1T，其中T是一個維度為K×K的非奇異矩陣，將Es1=Ω1T代入式(19)可得

對式(20)中T-1φuT=(K1Es1)+K2Es1進行特征分解，可得特征值為φu的對角元素，記第k個特征值為αk，k=1,2,…,K，那么虛擬子陣1的角度估計值可以表示為

同理，對于虛擬子陣2，對協(xié)方差矩陣R2進行酉變換后得到R2,u，參考式(18—21)，可以得到虛擬子陣2的角度估計結(jié)果，k=1,2,…,K。注意到虛擬子陣1和虛擬子陣2陣元間距均大于半波長，因此得到的角度估計值均為模糊值，需要解模糊處理。

2.2.2 角度解模糊

對于第k個入射信源θk，假設(shè)虛擬子陣1和子陣2中分別存在估計模糊角度值θ'k和θ″k，并有以下關(guān)系式[8]

式中：θk為第 k個入射信源的真實值，k=1,2,…,K。根據(jù) d1=Nλ/2和d2=Mλ/2，將式(14)和(15)分別代入式(22)和(23)，可以得到式中：k1，k2為整數(shù)，因為|sin(θk)-sin(θ'k)|＜2，|sin(θk)-sin(θ″k)|＜2，所以k1和k2的取值范圍分別為k1=-(N-1),…,0,…,(N-1)和k2=-(M-1),…,0,…,(M-1)?？紤]式(24)和(25)中θ'k，θ″k和θ″k可互換，因此虛擬子陣1和子陣2關(guān)于θk的模糊值個數(shù)分別為N和M(包含真實值)。其中，聯(lián)立式(24)和(25)可以得到

因為M和N互為質(zhì)數(shù)，當且僅當k1=k2=0時上式成立，也就是說，當且僅當θk=θ'k=θ″k時，式(26)成立。可以得出結(jié)論，虛擬子陣1和子陣2中關(guān)于θk的所有模糊值中只有一個是相等的，也就是真實估計值。

2.3 所提DOA估計算法步驟

本文提出的CA-UESPRIT算法步驟如下：

(1)得到互質(zhì)線陣中第p(p=1,2,…,M+N-1)個陣元接收信號xp(t)，利用不同時長間隔τ求陣元接收信號間的互相關(guān)函數(shù)Rxpx1(τ)。

(2)根據(jù)式(8—11)得到共軛增廣虛擬陣列接收信號y。

(3)利用UESPRIT算法分別對虛擬子陣1和子陣2接收信號y1和y2進行DOA估計，步驟如式(16—21)，得到兩個子陣的模糊角估計值

3 性能分析

3.1 最大可辨識信源數(shù)

根據(jù)文獻[1,2],一個具有N個傳感器的陣列最多能辨識N-1信源。因此，互質(zhì)陣列解模糊方法下最大可辨識信源數(shù)取決于傳感器數(shù)較小的子陣。令T=Min(M,N)，其中Min(M,N)表示取M和N中的最小值。根據(jù)文獻[17]，傳統(tǒng)互質(zhì)線陣下基于互質(zhì)特性的聯(lián)合UESPRIT算法，最大可辨識信源數(shù)為T-1。然而，根據(jù)2.2節(jié)描述，本文提出的CA-UESPRIT算法通過共軛增廣特性使得虛擬子陣陣元數(shù)分別達到2M-1和2N-1，因此CA-UESPRIT算法的最大可辨識信源數(shù)為Min(2M-1,2N-1)-1，也就是2T-2。

3.2 復雜度分析

以復乘次數(shù)來評估DOA估計算法的復雜度，CA-UESPRIT算法主要運算復雜度包括：構(gòu)造虛擬陣列接收信號y所需復雜度為O{Jw(2(M+N-1)-1)}，其中Jw為偽快拍數(shù)；得到虛擬子陣1和子陣2協(xié)方差矩陣R1和R2的復雜度分別為O{Jw(2M-1)2}和O{Jw(2N-1)2}；得到實數(shù)協(xié)方差矩陣R1,u和R2,u的復雜度分別為O{2(2M-1)3}和O{2(2N-1)3}。因此，總運算復雜度為O{Jw(2M+2N-3+(2M-1)2+(2N-1)2)+2(2M-1)3+2(2N-1)3}。本文利用UESPRIT算法實現(xiàn)DOA估計，不需要計算協(xié)方差矩陣特征分解以及后續(xù)總體最小二乘法的復雜度，因此復雜度相比于ESPRIT算法要低。

3.3 優(yōu)點總結(jié)

本文提出的CA-UESPRIT算法具有如下優(yōu)點：

(1)相比于傳統(tǒng)互質(zhì)線陣下基于互質(zhì)特性的聯(lián)合UESPRIT算法，CA-UESPRIT算法空間自由度更高。

(2)CA-UESPRIT算法具有比聯(lián)合UESPRIT算法更優(yōu)的DOA估計性能。

(3)CA-UESPRIT算法通過引入酉變換矩陣降低了ESPRIT算法中的復乘運算次數(shù)，算法復雜度較低。

4 仿真結(jié)果

本文采用1000次蒙特卡洛仿真的求根均方誤差(Root mean squares error,RMSE)評估DOA估計算法性能，定義RMSE為

圖5為本文DOA估計算法在信噪比SNR=0dB條件下角度估計圖。圖5中橫坐標為仿真次數(shù)，縱坐標為估計出的角度值。仿真參數(shù)設(shè)置為陣元數(shù)M=3,N=5，快拍數(shù)J=200。其中圖5(a)入射信源數(shù) K=3, θ =[25°,30°,35°],圖 5(b)入射信源數(shù) K=4, θ=[25°,30°,35°,40°]。從圖 5中可以看出，本文提出的CA-UESPRIT算法可以準確估計出信源角度值，且仿真結(jié)果驗證了3.1節(jié)中提出的本文算法可辨識信源數(shù)能達到2T-2=4，其中T=Min(M,N)。而根據(jù)文獻[17]，傳統(tǒng)互質(zhì)線陣下基于互質(zhì)特性的聯(lián)合UESPRIT算法，可辨識信源數(shù)僅能達到T-1=2。此外，對比圖5(a)和圖5(b),可以看出隨著信源數(shù)的增加，算法估計角度精確性變差，圖5(a)中RMSE值對應(yīng)于不同仿真次數(shù)分別為0.0842,0.0831,0.0838,0.0866,0.0827，如圖 5(b)所示,RMSE 值分別為 0.1330,0.1295,0.1372，0.1276,0.1372。可以得出信源數(shù)增加使得RMSE性能降低。

圖6為互質(zhì)線陣下不同DOA估計算法隨信噪比SNR變化下的RMSE性能對比圖，仿真參數(shù)均設(shè)置為陣元數(shù)M=3,N=5和K=2，θ=[20°,40°]。算法包括基于共軛增廣的CA-UESPRIT算法、CA-ESPRIT算法和CA-PM算法以及傳統(tǒng)互質(zhì)線陣下基于互質(zhì)特性的聯(lián)合UESPRIT算法[17]、ESPRIT算法[9]和PM算法。從圖6可以看出，隨著信噪比的增加，所提CA-UESPRIT算法的RMSE性能逐步變優(yōu)，也就是說，信噪比的增加帶來角度估計精度的提升。對比不同的DOA估計算法性能，可以看出，所提的CA-UESPRIT算法DOA估計性能與CA-ESPRIT算法十分相近，但CA-UESPRIT算法將協(xié)方差矩陣從復數(shù)域轉(zhuǎn)化到實數(shù)域，所以算法復雜度更低，且CA-UESPRIT算法DOA估計性能要優(yōu)于CA-PM算法。此外，相比于互質(zhì)線陣下基于互質(zhì)特性的聯(lián)合DOA估計算法[9，17],本文提出基于共軛增廣特性DOA估計算法由于構(gòu)造的虛擬陣列信號具有更大的空間自由度，因此在RMSE性能全面占優(yōu)。

圖5 DOA估計結(jié)果散點圖Fig.5 Scatter of DOA estimation results

圖7 為互質(zhì)線陣下不同DOA估計算法隨快拍數(shù)J變化下的RMSE性能對比圖，仿真參數(shù)均設(shè)置為陣元數(shù)M=3,N=5和K=2，θ=[20°,40°]，算法類型與圖6中相同。從圖7可以看出，隨著快拍數(shù)的增加，所提的CA-UESPRIT算法DOA估計性能提升。此外，對比不同的DOA估計算法，可以看出基于共軛增廣類的DOA估計算法具有更優(yōu)的RMSE性能。此外，所提的CA-UESPRIT算法DOA估計性能優(yōu)于CA-PM算法且接近于CA-ESPRIT算法，但CA-UESPRIT算法復雜度較低(具體分析參考圖6分析結(jié)果)。

圖6 互質(zhì)陣列下不同算法RMSE性能對比圖(J=200)Fig.6 Comparison of RMSE performance with different algorithms for coprime array(J=200)

圖7 互質(zhì)陣列下不同算法RMSE性能對比圖(SNR=0dB)Fig.7 Comparison of RMSE performance with different algorithms for coprime array(SNR=0dB)

圖8 為本文提出的CA-UESPRIT算法下RMSE隨信噪比SNR和快拍數(shù)J變化的性能對比圖，仿真參數(shù)均設(shè)置為陣元數(shù)M=3,N=5和K=2，θ=[20°,40°]。從圖8中可以看出，隨著信噪比SNR和快拍數(shù)J的增大，所提的DOA估計算法角度估計性能變優(yōu)。

5 結(jié)束語

本文提出了互質(zhì)線陣中基于共軛增廣的DOA估計算法。該算法首先利用不同時長間隔下物理陣元接收信號間互相關(guān)函數(shù)的共軛對稱特性，構(gòu)造虛擬陣列接收信號以擴展陣列孔徑和提高空間自由度。然后采用基于互質(zhì)特性的聯(lián)合UESPRIT算法實現(xiàn)DOA估計。相比于傳統(tǒng)互質(zhì)線陣下聯(lián)合的UESPRIT算法，本文提出的基于共軛增廣的UESPRIT(CA-UESPRIT)算法角度估計性能更優(yōu)，空間自由度更高。此外，相比于ESPRIT算法，CA-UESPRIT算法通過引入酉變換矩陣降低復乘運算次數(shù)，算法復雜度較低。探索更高效的DOA估計算法將會是下一步的研究方向。

圖8 不同快拍數(shù)條件下本文DOA估計算法RMSE性能對比圖Fig.8 Comparison of RMSE performance of DOA estimation algorithm under different snapshots