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

王 凌, 潘 華, 趙 維

(海軍研究院, 上海 200235)

0 引 言

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

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

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

1 陣列信號(hào)模型及問(wèn)題提出

1.1 陣列模型建立

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

圖1 陣列模型Fig.1 Array model

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

(1)

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

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

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

則此時(shí)整個(gè)矩形陣列的輸出可表示為

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個(gè)入射信源;N(t)為噪聲矢量,其方差為σ2;A為整個(gè)陣列系統(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 問(wèn)題提出

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

(11)

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

(12)

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

2.1 算法實(shí)現(xiàn)

定義如下3個(gè)選擇矩陣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中空心圓所示陣列劃分為兩個(gè)子陣X1和X2,并定義子陣y中空心圓所示陣列為子陣Y,則子陣X1的接收數(shù)據(jù)可以表示為

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

(18)

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

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

(19)

式中:Q1C可以寫(xiě)作

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可以寫(xiě)作

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)一步可寫(xiě)為

Q2CA=Q1CAΨ

(31)

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

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

Q3CAS(t)+Q3N(t)

(32)

式中:Q3C可以寫(xiě)作

Q3C=[0,J1C2,J1C1,J1C2]

(33)

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

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

(34)

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

Q3CA=Q1CAΦ2

(35)

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

(36)

(37)

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

US2=US1T-1ΨT=US1φ1

(38)

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

(39)

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

(40)

(41)

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

(42)

2.2 算法步驟

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

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

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

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

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

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

仿真 13個(gè)空間不相關(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àn)中對(duì)比本文算法和不考慮互耦效應(yīng)時(shí)的經(jīng)典2D -ESPRIT算法以及受互耦影響的2D -ESPRIT算法。

估計(jì)誤差隨SNR變化曲線如圖2所示，可以看出本文提出的算法估計(jì)性能相較于直接受互耦擾動(dòng)的2D -ESPRIT算法,估計(jì)性能得到了大幅提高,當(dāng)SNR高于5 dB時(shí)，估計(jì)性能和標(biāo)準(zhǔn)無(wú)互耦誤差的2D -ESPRIT算法接近,仿真結(jié)果證明了本文構(gòu)造的子陣能夠抑制互耦對(duì)旋轉(zhuǎn)不變關(guān)系的影響。從圖中受互耦影響的2D -ESPRIT算法仿真曲線可以看出,互耦對(duì)子陣陣列流型的擾動(dòng)顯著改變了旋轉(zhuǎn)不變關(guān)系,即使SNR增大,估計(jì)誤差也會(huì)收斂在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,得到的估計(jì)誤差隨陣元數(shù)和SNR變化曲線如圖3所示。

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

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

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

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

4 結(jié) 論

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