互耦效應(yīng)下低復(fù)雜度的二維DOA估計算法

王 凌, 潘 華, 趙 維

(海軍研究院, 上海 200235)

0 引 言

空間譜估計技術(shù)區(qū)別于傳統(tǒng)的測向技術(shù),具有超分辨率測向能力,經(jīng)過近幾十年的的長足發(fā)展,在各領(lǐng)域取得了豐碩成果[1-6]。空間目標(biāo)的二維波達(dá)方向(two-dimension direction of arrival, 2D -DOA)因能更精細(xì)地劃分空間信道,具有更精確的定位目標(biāo),成為空間譜估計技術(shù)中的一個重要分支。目前已經(jīng)形成以二維多重信號分類(two-dimension multiple signal classification algorithm, 2D-MUSIC)算法和二維基于旋轉(zhuǎn)不變技術(shù)的信號參數(shù)估計(two-dimension estimating signal parameter via rotational invariance techniques, 2D -ESPRIT)算法為代表的二維子空間類算法理論體系[7-13]。但此類算法在陣列誤差存在時估計性能明顯下降,甚至失效。2D -MUSIC算法利用關(guān)系式a(α,β)UN來構(gòu)造針狀的二維譜峰,但當(dāng)考慮互耦時,實際的方向向量已經(jīng)變成Ca(α,β),此時仍然使用a(α,β)UN來進(jìn)行譜峰搜索就會出現(xiàn)較大誤差。在實際工程中,可以通過實測或預(yù)估得到互耦系數(shù)矩陣C,文獻(xiàn)[14]證明此時用修正后的方向向量Ca(α,β)來進(jìn)行二維譜峰搜索，同樣也能得到正確的2D -DOA,但此時需要預(yù)先估計互耦系數(shù),無法滿足實時性要求。2D -ESPRIT算法則利用子陣列信號中子空間之間的旋轉(zhuǎn)不變關(guān)系來求得2D -DOA,當(dāng)互耦存在時,這種旋轉(zhuǎn)不變性將被嚴(yán)重破壞,導(dǎo)致基于旋轉(zhuǎn)不變技術(shù)的信號參數(shù)估計(estimating signal parameter via rotational invariance techniques, ESPRIT)算法失效。2D -ESPRIT算法無需二維譜峰搜索,計算量相較于2D -MUSIC算法大幅降低,但對于2D -ESPRIT算法來說,即使已知互耦信息,也無法構(gòu)造出具有旋轉(zhuǎn)不變關(guān)系的信號子空間,因此在互耦應(yīng)用背景下，如何構(gòu)造旋轉(zhuǎn)不變性成為值得研究的課題。

另外,國內(nèi)外學(xué)者針對2D -DOA估計和互耦校正作了大量的研究[15-21]。這類校正算法大致可以分為兩類,一類是有源離線校正方法[15-16],該類方法需要在空間設(shè)置額外的校準(zhǔn)源,增加了操作負(fù)擔(dān),且當(dāng)校正源方位存在誤差時,對校正效果影響明顯。另一類是在線自校正方法[17-21],該類方法將互耦系數(shù)、DOA等未知參數(shù)轉(zhuǎn)化為多參數(shù)的聯(lián)合估計問題,存在多次迭代尋優(yōu)、容易陷入局部最優(yōu)、計算量龐大、不易工程實施等缺點。雖然可利用互耦系數(shù)矩陣的稀疏特性[20]或Toeplitz特性[21]來優(yōu)化和簡化迭代過程,但沒有根本上解決上述問題。

針對互耦效應(yīng)對旋轉(zhuǎn)不變關(guān)系的破壞以及傳統(tǒng)互耦校正算法存在的上述缺點,本文提出了一種新的完全解互耦2D -ESPRIT(mutual coupling 2D-ESPRIT, MC-2D -ESPRIT)算法,該算法通過構(gòu)建互耦效應(yīng)下仍保持旋轉(zhuǎn)不變關(guān)系的子陣,將旋轉(zhuǎn)不變技術(shù)推廣至互耦應(yīng)用背景。由于該算法無需二維譜峰搜索和多次迭代,且互耦抑制過程無需預(yù)知任何互耦信息,相較于傳統(tǒng)算法,該算法計算量小,易于工程實現(xiàn)。仿真結(jié)果驗證了該算法能抑制互耦,2D -DOA估計性能明顯優(yōu)于直接受互耦擾動的2D -ESPRIT算法,且接近于無互耦的標(biāo)準(zhǔn)2D -ESPRIT算法。

1 陣列信號模型及問題提出

1.1 陣列模型建立

所建立的陣列模型如圖1所示，考慮空間有M個不相關(guān)的窄帶信源(s1(t),s2(t),…,sM(t))分別從矢量角(θ1,θ2,…,θM)入射到由4個相互平行的均勻線陣組成的平面矩形陣列上,其中θk=(αk,βk,γk),αk,βk,γk分別為第k個信源入射方向與x軸、y軸和z軸的夾角,由于其中只有兩個角度獨立,故空間信源入射矢量角可表示為θk=(αk,βk),k=1,2,…,M。矩形陣列各子陣陣元間距為dx,子陣之間間距為dy。

圖1 陣列模型Fig.1 Array model

對于矩形陣列互耦模型的建立,本文考慮陣元周圍8個陣元對其產(chǎn)生的互耦效應(yīng)影響,同一子陣相鄰陣元間互耦系數(shù)定義為cx,相鄰子陣陣元間互耦系數(shù)定義為cy和cxy。此時矩形陣列互耦系數(shù)矩陣C可以表示為

(1)

式中:C1和C2為滿足Toeplitz矩陣形式的N×N維子互耦系數(shù)矩陣,且

{C1=toeplitz{[1,cx,0,…,0]}

C2=toeplitz{[cy,cxy,0,…,0]}

則此時整個矩形陣列的輸出可表示為

X(t)=[z1(t),…,zN(t),x1(t),…,xN(t),

y1(t),…,yN(t),l1(t),…,lN(t)]T=

CAS(t)+N(t)

(2)

式中：S(t)=[s1(t),s2(t),…,sM(t)]T為M個入射信源;N(t)為噪聲矢量,其方差為σ2;A為整個陣列系統(tǒng)的陣列流型矩陣,可表示為

(3)

其中

A1=[a(α1),a(α2),…,a(αM)]

(4)

a(αk)=[1,u(αk),u2(αk),…,uN-1(αk)]T,k=1,2,…,M

(5)

Φ1=diag[v-1(β1),v-1(β2),…,v-1(βM)]

(6)

Φ2=diag[v(β1),v(β2),…,v(βM)]

(7)

Φ3=diag[v2(β1),v2(β2),…,v2(βM)]

(8)

(9)

(10)

1.2 問題提出

2D -ESPRIT是ESPRIT算法的二維推廣形式,用2D -ESPRIT算法實現(xiàn)2D -DOA估計的關(guān)鍵是找到或構(gòu)造出具有選擇不變關(guān)系的子陣。經(jīng)典的2D -ESPRIT算法利用圖1中空心圓代表的平行陣列,構(gòu)造出3個具有選擇不變關(guān)系的子陣列X1、X2和Y。當(dāng)不考慮互耦效應(yīng)時,3個子陣信號的子空間滿足:

(11)

Ψ=diag[u(α1),u(α2),…,u(αM)]

(12)

2 互耦效應(yīng)下的2D -ESPRIT算法

2.1 算法實現(xiàn)

定義如下3個選擇矩陣Q1、Q2和Q3,分別表示為

Q1=[0,J1,0,0]

(13)

Q2=[0,J2,0,0]

(14)

Q3=[0,0,J1,0]

(15)

式中:0、J1和J2為(N-3)×N維矩陣,J1、J2可表示為

(16)

(17)

在圖1中將子陣x中空心圓所示陣列劃分為兩個子陣X1和X2,并定義子陣y中空心圓所示陣列為子陣Y,則子陣X1的接收數(shù)據(jù)可以表示為

X1(t)=[x2,x3,…,xN-2]T=Q1X(t)

(18)

將式(2)代入式(18),進(jìn)一步展開,可以得到如下關(guān)系:

X1(t)=Q1X(t)=Q1CAS(t)+Q1N(t)

(19)

式中:Q1C可以寫作

Q1C=[J1C2,J1C1,J1C2,0]

(20)

式中:J1C1和J1C2可表示為

(21)

(22)

注意到J1C1和J1C2最后一列全為0元素,根據(jù)矩陣?yán)碚?存在下述關(guān)系式:

Bm×(n-1)D(n-1)×p

(23)

則Q1CA表達(dá)式如下:

Q1CA=J1C2A1Φ1+J1C1A1+J1C2A1Φ2=

(24)

X2(t)=[x3,x4,…,xN-1]T=Q2X(t)=

Q2CAS(t)+Q2N(t)

(25)

式中:Q2C可以寫作

Q2C=[J2C2,J2C1,J2C2,0]

(26)

式中:J2C1和J2C2形式如下:

(27)

(28)

注意到J2C1和J2C2第一列全為0元素,根據(jù)矩陣?yán)碚?存在下述關(guān)系式:

Bm×(n-1)D(n-1)×p

(29)

則Q2CA表達(dá)式如下:

Q2CA=J2C2A1Φ1+J2C1A1+J2C2A1Φ2=

(30)

式中:A12為A1的后N-1行。由于A11和A12之間存在關(guān)系式A12=A11Ψ,則式(30)進(jìn)一步可寫為

Q2CA=Q1CAΨ

(31)

子陣Y的接收數(shù)據(jù)可以表示為

Y(t)=[y2,y3,…,yN-2]T=Q3X(t)=

Q3CAS(t)+Q3N(t)

(32)

式中:Q3C可以寫作

Q3C=[0,J1C2,J1C1,J1C2]

(33)

則Q3CA表達(dá)式如下：

Q3CA=J1C2A1+J1C1A1Φ2+J1C2A1Φ3=

(34)

由于Φ1Φ2=I,Φ2Φ2=Φ3,則式(34)進(jìn)一步可寫為

Q3CA=Q1CAΦ2

(35)

將3個子陣接收數(shù)據(jù)進(jìn)行合并,即

(36)

(37)

得到各子陣信號子空間之間的關(guān)系為

US2=US1T-1ΨT=US1φ1

(38)

US3=US1T-1Φ2T=US1φ2

(39)

將式(38)和式(39)求最小二乘解,可得到

(40)

(41)

對φ1和φ2特征分解即可得到旋轉(zhuǎn)不變關(guān)系矩陣Ψ和Φ2,利用式(9)和式(10)對應(yīng)關(guān)系,就可以得到M個入射信源分別與x軸和y軸夾角αk,βj(k,j=1,2，…,M)。式(40)和式(41)的特征分解是獨立進(jìn)行的,設(shè)T1和T2分別為φ1和φ2的特征向量矩陣。由于信源sk(t)與x軸夾角αk對應(yīng)的特征向量和其與y軸夾角βk對應(yīng)的特征向量是強(qiáng)相關(guān)的,因此可以利用式(42)中矩陣G每列的最大值來調(diào)整βj的順序,從而到達(dá)匹配。矩陣G可表示為

(42)

2.2 算法步驟

本文提出的完全解互耦2D -ESPRIT算法可以總結(jié)為以下4個步驟。

步驟 1從圖1所示陣列系統(tǒng)中提取3個子陣的接收數(shù)據(jù)X1(t)、X2(t)和Y(t)。

步驟 3對信號子空間US分塊,并按照式(40)和式(41)求得φ1和φ2,進(jìn)一步求得Ψ和Φ2。

步驟 4根據(jù)式(9)和式(10)對應(yīng)關(guān)系求得αk,βj(k,j=1,2,…,M)并用式(42)來配對。

3 算法性能評估與數(shù)值仿真

前文從理論上驗證了通過平面矩形陣列中對3個子陣的選取,能夠抑制互耦效應(yīng)對旋轉(zhuǎn)不變關(guān)系的擾動,本節(jié)將對提出的完全解互耦2D -ESPRIT算法進(jìn)行數(shù)值仿真驗證,定量分析算法的估計性能。仿真中設(shè)定互耦系數(shù)cx=cy=0.433 1+0.251 2i,cxy=0.141 2+0.141 2i,信源DOA估計精度采用均方根誤差(root mean square error, RMSE)。

仿真 13個空間不相關(guān)信源分別從波達(dá)方向(30°,90°)、(60°,60°)和(85°,75°)入射至陣列系統(tǒng),設(shè)置N=10,快拍數(shù)為1 000次。設(shè)置信噪比(signal to noise ratio, SNR)為0～20 dB,試驗中對比本文算法和不考慮互耦效應(yīng)時的經(jīng)典2D -ESPRIT算法以及受互耦影響的2D -ESPRIT算法。

估計誤差隨SNR變化曲線如圖2所示，可以看出本文提出的算法估計性能相較于直接受互耦擾動的2D -ESPRIT算法,估計性能得到了大幅提高,當(dāng)SNR高于5 dB時，估計性能和標(biāo)準(zhǔn)無互耦誤差的2D -ESPRIT算法接近,仿真結(jié)果證明了本文構(gòu)造的子陣能夠抑制互耦對旋轉(zhuǎn)不變關(guān)系的影響。從圖中受互耦影響的2D -ESPRIT算法仿真曲線可以看出,互耦對子陣陣列流型的擾動顯著改變了旋轉(zhuǎn)不變關(guān)系,即使SNR增大,估計誤差也會收斂在12°左右,因此互耦導(dǎo)致經(jīng)典的ESPRIT算法失效。

圖2 RMSE隨SNR變化曲線Fig.2 RMSE curve with SNR

仿真 2固定快拍數(shù)為1 000次,改變子陣陣元數(shù),使N分別等于10、11、12、13和15,設(shè)置SNR為0～20 dB,得到的估計誤差隨陣元數(shù)和SNR變化曲線如圖3所示。

圖3 RMSE隨陣元數(shù)和SNR變化曲線Fig.3 RMSE curve with array number and SNR

從圖3所示結(jié)果易知,本文算法在子陣陣元數(shù)小于12時,低SNR的估計性能出現(xiàn)較大誤差,但只要SNR高于4 dB,2D -DOA的估計誤差能控制在1°以內(nèi)。當(dāng)子陣陣元數(shù)大于12時,明顯改善了低SNR時算法估計性能,圖3也進(jìn)一步驗證了本文算法對互耦的抑制能力。

仿真 3固定陣元數(shù)為15,改變陣列系統(tǒng)快拍數(shù),使快拍數(shù)分別等于700、800、1 000、1 200和1 500,設(shè)置SNR為0～20 dB,得到估計誤差隨SNR和快拍數(shù)變化曲線如圖4所示。從結(jié)果可以看出,當(dāng)快拍數(shù)高于1 000次時,算法估計誤差能控制在2°以內(nèi),當(dāng)SNR高于10 dB時算法估計誤差已經(jīng)小于0.5°。

圖4 RMSE隨SNR變化曲線Fig.4 RMSE curve with SNR

4 結(jié) 論

傳統(tǒng)實現(xiàn)互耦抑制的思路多是將互耦系數(shù)和2D -DOA作為整體進(jìn)行估計,并轉(zhuǎn)化為多參量非線性優(yōu)化問題,在估計得到互耦參數(shù)后再利用2D -MUSIC算法進(jìn)行譜峰搜索,這樣無疑會導(dǎo)致計算量龐大等問題。而本文算法計算量僅對協(xié)方差矩陣進(jìn)行一次特征分解。本文通過選取矩形陣列中3個在互耦效應(yīng)影響下仍保持旋轉(zhuǎn)不變關(guān)系的子陣列,無需任何互耦信息,將旋轉(zhuǎn)不變思想推廣至互耦擾動下的2D -DOA估計,提出了一種互耦效應(yīng)下的低計算復(fù)雜度MC-2D -ESPRIT,數(shù)值仿真結(jié)果驗證了該算法對互耦的抑制能力,估計精度能夠接近標(biāo)準(zhǔn)無互耦誤差的2D -ESPRIT算法。