郭 拓， 王英民， 張立琛

(西北工業(yè)大學(xué) 航海學(xué)院，西安 710072)

波達(dá)方向估計(jì)(Direction of Arrival，DOA)是雷達(dá)、聲納、通信、無源定位等陣列信號(hào)處理中的一個(gè)重要研究?jī)?nèi)容之一。陣列波達(dá)方向估計(jì)經(jīng)過近幾十年的發(fā)展，已經(jīng)形成了兩類技術(shù)，即基于波束形成的技術(shù)和基于子空間分解的技術(shù)?；诓ㄊ纬杉夹g(shù)的典型方法有常規(guī)波束形成器(Conventional Beam Forming，CBF)和最小方差無畸變響應(yīng)(Minimum Variance Distortionless Response,MVDR)波束形成器；多重信號(hào)分類(Multiple Signal Classification，MUSIC)算法是經(jīng)典的子空間類方法[1]。上述DOA估計(jì)算法都需要估計(jì)基陣接收到數(shù)據(jù)的統(tǒng)計(jì)協(xié)方差矩陣，然而在實(shí)際數(shù)據(jù)處理中都以樣本協(xié)方差矩陣代替統(tǒng)計(jì)協(xié)方差矩陣，在小快拍情況下估計(jì)出的統(tǒng)計(jì)協(xié)方差矩陣誤差很大，將會(huì)導(dǎo)致傳統(tǒng)DOA估計(jì)算法不能準(zhǔn)確估計(jì)目標(biāo)方位。

為了獲得與統(tǒng)計(jì)協(xié)方差矩陣更加相近的樣本協(xié)方差矩陣，通常情況下會(huì)對(duì)樣本協(xié)方差矩陣進(jìn)行修正，主要的修正方法包括：對(duì)角加載方法[2]，收縮因子方法[3]等，通過對(duì)樣本協(xié)方差矩陣的特征值進(jìn)行加權(quán)操作以改進(jìn)其方位估計(jì)的分辨性能；前后向平均技術(shù)[4]，空間光滑技術(shù)[5]，通過移動(dòng)窗口虛擬出更多的快拍以達(dá)到改進(jìn)樣本協(xié)方差矩陣估計(jì)的目的；偽隨機(jī)噪聲再采樣技術(shù)[6-8]則是人為生成偽隨機(jī)噪聲信號(hào)加在接收到的信號(hào)上，從而對(duì)樣本協(xié)方差樣本協(xié)方差矩陣進(jìn)行修正；而子空間泄漏分析方法[9]總共有兩個(gè)步驟：①采用根MUSIC算法估計(jì)目標(biāo)方位；②根據(jù)估計(jì)所得的目標(biāo)方位修正樣本協(xié)方差矩陣。

上述這些方法的思路都是設(shè)法改進(jìn)樣本協(xié)方差矩陣，以達(dá)到其與大快拍數(shù)時(shí)的樣本協(xié)方差矩陣更加相近，而沒有考慮樣本協(xié)方差矩陣與快拍數(shù)的關(guān)系，以及不同快拍數(shù)時(shí)協(xié)方差矩陣本身所具有的不同特性。

本文將不沿著上文所述那樣通過改進(jìn)樣本協(xié)方差矩陣來改進(jìn)方位估計(jì)效果的思路，而是轉(zhuǎn)換研究思路，研究不同快拍數(shù)時(shí)樣本協(xié)方差矩陣本身具有的一些特性，比如譜的特性，具體的就是通過深入分析相干目標(biāo)與獨(dú)立目標(biāo)的樣本協(xié)方差矩陣在不同快拍下所具有的各異的譜分離(Spectral Separation)特性，提出一種基于小快拍條件下基于樣本協(xié)方差矩陣譜分離特性的DOA估計(jì)方法。

1 陣列接收信號(hào)模型

如圖1所示，水平線陣由N個(gè)陣元組成，假設(shè)有D個(gè)相互獨(dú)立的平面波入射，陣列第k次快拍X(k)的輸出如式(1)所示。

圖1 陣列信號(hào)模型Fig.1 Array signal model

(1)

式中：X(k)=[x1(k),x2(k), …,xN(k)]T，A(θ)=[a(θ1),a(θ2),…,a(θD)]；a(θi)=[1, ej2πd cos θi/λ, …, ej2π(N-1)dcos θi/λ]為第i個(gè)入射平面波的陣列流形；陣元間距為d，d為波長(zhǎng)λ的一半；θi為第i個(gè)平面波入射方向與水平面的夾角；A(θ)為陣列流形矩陣；S(k)=[s1(k),s2(k), …,sD(k)]T,si(k)為第i個(gè)平面波信號(hào)的第k次快拍；其中T為轉(zhuǎn)置運(yùn)算。N(k)=[n1(k),n2(k), …,nN(k)]為陣列接收到的加性高斯白噪聲。

當(dāng)信號(hào)被假設(shè)具有零均值時(shí)，陣列接收信號(hào)的協(xié)方差矩R=E[X(k)·X(k)H]，即協(xié)方差矩陣需要通過求集總平均而得到。而實(shí)際應(yīng)用中只能得到有限的采樣數(shù)據(jù)，如有M次觀測(cè)樣本X1,X2，…，XM用這M個(gè)樣本來估計(jì)樣本協(xié)方差矩陣如式(2)，H為復(fù)共軛轉(zhuǎn)置。

(2)

2 基于小快拍協(xié)方差矩陣的主特征空間波達(dá)方向估計(jì)方法

本節(jié)將首先深入分析當(dāng)目標(biāo)信號(hào)分別為相干信號(hào)與獨(dú)立信號(hào)時(shí)樣本協(xié)方差矩陣的譜分離特性，然后在此基礎(chǔ)上提出樣本協(xié)方差矩陣主特征空間的概念，最后提出基于小快拍的主特征空間波達(dá)方向估計(jì)方法。

2.1 相干信號(hào)樣本協(xié)方差矩陣特征的譜分離特性分析

從子空間的角度來解釋上述相干源情況下譜分離特性，就是在y較小時(shí)，樣本協(xié)方差矩陣得到較好的估計(jì)，相干信號(hào)的信息全部集中在一個(gè)最大特征值所對(duì)應(yīng)的子空間上；而在y較大時(shí)，如y為1此時(shí)陣元數(shù)等于快拍數(shù)，樣本協(xié)方差矩陣沒得到很好的估計(jì)，此時(shí)信號(hào)子空間發(fā)生了泄漏，使得最大特征值之前出現(xiàn)一些連續(xù)的較大特征值，但是總體而言相干信號(hào)的主要信息還是集中于最大特征值處。本文將這個(gè)最大特征值所對(duì)應(yīng)的子空間稱為主特征空間(Main Feature Space, MFS)，而不是MUSIC算法所假設(shè)的那樣，即信號(hào)子空間就是s個(gè)信源所對(duì)應(yīng)的s個(gè)較大特征值所構(gòu)成的空間，通過我們上述對(duì)相干源情況下不同快拍數(shù)時(shí)譜分離特性的仿真與分析，對(duì)于相干源來說顯然信號(hào)信息主要集中在主特征空間而不是MUSIC算法所使用的s個(gè)較大特征值所構(gòu)成的空間，這就很好的解釋了MUSIC算法對(duì)相干目標(biāo)DOA估計(jì)時(shí)為何失效。

(a) 相干源譜分離性

(b) 相干源譜分離性局部放大圖2 相干源譜分離性分析Fig.2 Analysis of spectra separation for coherent sources

2.2 獨(dú)立信號(hào)樣本協(xié)方差矩陣特征的譜分離特性分析

圖3 獨(dú)立源的譜分離性Fig.3 Spetral separation of independent source

由“2.1”節(jié)分析可知，對(duì)于相干源，主特征空間中總是有個(gè)信源的信息，只是在快拍數(shù)小時(shí)，主特征空間中有一定信號(hào)信息泄漏入非主特征空間中，做DOA估計(jì)時(shí)會(huì)導(dǎo)致旁瓣會(huì)比快拍數(shù)大時(shí)的高。

對(duì)于獨(dú)立源，在小快拍情況下，主特征空間中會(huì)有所有目標(biāo)信源的信息；快拍數(shù)大時(shí)，主特征空間中會(huì)完全不存在一部分源的信息，導(dǎo)致DOA估計(jì)中沒有這些目標(biāo)。

總之不管是相干源還是獨(dú)立源，小快拍時(shí)主特征空間中包含所有信號(hào)的信息(這里的小快拍是指快拍數(shù)與陣元數(shù)相當(dāng)，具體是指快拍數(shù)是陣元數(shù)的2倍以下，即y>0.5的情況)，所以可以利用小快拍時(shí)樣本協(xié)方差矩陣的主特征空間進(jìn)行方位估計(jì)。下文將根據(jù)上述分析的譜分離特性提出基于主特征空間的DOA估計(jì)方法。

2.3 基于主特征空間的DOA估計(jì)方法

(3)

(4)

根據(jù)子空間理論，信號(hào)子空間與噪聲子空間正交，且導(dǎo)向向量張成的空間與信號(hào)子空間張成的空間是同一空間，即Ena(θ)=0，Esa(θ)，其中En為噪聲子空間，Es為信號(hào)子空間。

主特征空間中包含的主要是信號(hào)信息，即主特征子空間Emain中主要是Es中的信息。根據(jù)分析，在小快拍情況下，不管是相干源還是獨(dú)立源，主特征空間中都包括所有信號(hào)的信息。這樣，根據(jù)導(dǎo)向向量與信號(hào)子空間形成的是同一空間，則Emaina(θ)→1。在極限情形下如果噪聲子空間En的信息沒有泄漏入Emain，則Emaina(θ)=1，那么arccos(Emaina(θ))=0。因此利用導(dǎo)向向量與信號(hào)子空間形成的是同一空間這一特性，定義arccos(Emaina(θ))的倒數(shù)為主特征空間方法的DOA估計(jì)幅度，且使用Emain與a(θ)的l1范數(shù)來歸一化，即

(5)

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

3.1 相干信號(hào)DOA估計(jì)

由于在水池實(shí)驗(yàn)中陣元數(shù)有限采用了被動(dòng)合成孔徑技術(shù)進(jìn)行了孔徑擴(kuò)展，在仿真中采用被動(dòng)合成孔徑算法ETAM[14](Extended Towed Array Measurements)。仿真采用一個(gè)半波長(zhǎng)分布的6元均勻線陣，陣列運(yùn)動(dòng)速度為0.07 m/s，入射的信號(hào)為4個(gè)遠(yuǎn)場(chǎng)窄帶信號(hào)，頻率為6.25 kHz，采樣頻率50 kHz，4個(gè)信號(hào)的入射角分別為-20°，-10°，0°，10°，信噪比SNR皆為0，參考聲速為1 500 m/s，噪聲為復(fù)高斯白噪聲，快拍數(shù)分別取30、12 000，ETAM算法經(jīng)過8次合成，合成24個(gè)虛擬陣元，加6個(gè)物理陣元，則總線陣為30個(gè)陣元。

合成孔徑后30元陣取30個(gè)快拍和12 000個(gè)快拍時(shí)，所提出算法給出的目標(biāo)DOA估計(jì)對(duì)比，如圖4所示。由圖4可知，無論快拍數(shù)大小，主特征空間方法都可以分辨出相干目標(biāo)，只是在大快拍情況下旁瓣會(huì)低一些。圖5給出使用相同的數(shù)據(jù)，采用MUSIC算法獲得的目標(biāo)DOA估計(jì)。由圖5可知，MUSIC算法在小快拍情況下，會(huì)出現(xiàn)偽目標(biāo)，不能正確分辨目標(biāo)。比較可知，本文提出的方法可以在小快拍情況下，對(duì)相干目標(biāo)

圖4 相干源主特征空間DOA估計(jì)Fig.4 DOA of coherent source with main feature space

圖5 相干源MUSIC算法DOA估計(jì)Fig.5 DOA of coherent source with MUSIC

進(jìn)行很好的分辨，且不需要像MUSIC算法那樣需預(yù)先擁有目標(biāo)個(gè)數(shù)的先驗(yàn)信息，這是相干源時(shí)，樣本協(xié)方差矩陣的譜分離特性帶來的好處，即信號(hào)信息全部集中于主特征空間，只有很少的信號(hào)泄漏于非主特征空間中。

3.2 獨(dú)立信號(hào)DOA估計(jì)

本節(jié)仿真條件與“3.1”節(jié)相同，快拍數(shù)分別取30， 12 000，使用ETAM算法經(jīng)過8次合成，合成24個(gè)虛擬陣元，加6個(gè)物理陣元，則總線陣為30個(gè)陣元，信號(hào)源改為獨(dú)立源，頻率為6.24 kHz，6.245 kHz，6.25 kHz和6.254 kHz ，入射角還是-20°，-10°，0°，10°。主特征空間方法DOA估計(jì)，如圖6所示。

由圖6可知，當(dāng)目標(biāo)為獨(dú)立信源，小快拍(30快拍)情況下，主特征空間方法可以很好的估計(jì)出目標(biāo)方位；而在大快拍(12 000快拍)情況下，會(huì)發(fā)現(xiàn)不能正確分辨出所有目標(biāo)，有一個(gè)目標(biāo)(0°方向的目標(biāo))沒能估計(jì)出。此仿真驗(yàn)證了“2.2”節(jié)所分析的結(jié)論，即小快拍時(shí)，主特征空間中會(huì)有所有目標(biāo)信源的信息，大快拍時(shí)，主特征空間中會(huì)完全不存在某些源的信息，導(dǎo)致DOA估計(jì)中沒有得到這些目標(biāo)的方位。圖7為與圖6使用相同數(shù)據(jù)情況下，MUSIC算法在不同快拍情形下的DOA估計(jì)，發(fā)現(xiàn)與圖6的情況恰好相反，是在小快拍(30快拍)時(shí)不能準(zhǔn)確估計(jì)出4個(gè)目標(biāo)方位，而在多快拍時(shí)能準(zhǔn)確估計(jì)出目標(biāo)方位。

圖6 獨(dú)立源主特征空間DOA估計(jì)Fig.6 DOA of independent source with main feature space

圖7 獨(dú)立源MUSIC算法DOA估計(jì)Fig.7 DOA of independent source with MUSIC

4 水池驗(yàn)證

為了驗(yàn)證本文所提算法，在消聲水池(長(zhǎng)20 m,寬8 m,深7 m，6面消聲)進(jìn)行了相關(guān)實(shí)驗(yàn)。由于可用陣元數(shù)有限，故使用合成孔徑方法進(jìn)行孔徑的擴(kuò)展，然后進(jìn)行DOA估計(jì)。

實(shí)驗(yàn)采用由6個(gè)水聽器組成的均勻線陣，設(shè)計(jì)頻率為6.25 kHz，水平布放于水下2.5 m深處，陣元間距為半波長(zhǎng)0.12 m，陣列勻速運(yùn)動(dòng)速度為0.07 m/s，采用電機(jī)勻速拉動(dòng)。相干源時(shí)，3個(gè)發(fā)射源發(fā)射頻率皆為6.25 kHz的連續(xù)CW(單頻)信號(hào)；獨(dú)立源時(shí)，3個(gè)源頻率分別為6.24 kHz，6.25 kHz和6.254 kHz；與接收陣在同一深度，3個(gè)源按遠(yuǎn)場(chǎng)條件布放，與陣形成的角度分別為24°，33°，44° ；采樣頻率為50 kHz；參考聲速為1 500 m/s。ETAM算法合成6次，得到18個(gè)虛擬陣元，加6個(gè)物理陣元，共得到24個(gè)陣元，快拍數(shù)分別取24,12 000。實(shí)驗(yàn)系統(tǒng)如圖8所示。

圖8 水池實(shí)驗(yàn)系統(tǒng)Fig.8 Tank experiment system

4.1 相干信號(hào)DOA估計(jì)

圖9為水池實(shí)驗(yàn)對(duì)相干信號(hào)使用主特征空間DOA估計(jì)的結(jié)果，正如“2.1”節(jié)相干信號(hào)樣本協(xié)方差矩陣譜分離特性的分析以及“3.1”節(jié)仿真驗(yàn)證，無論快拍數(shù)大小，主特征空間方法都可以分辨出相干目標(biāo)，只是在大快拍情況下旁瓣會(huì)低一些，圖10中快拍數(shù)為12 000時(shí)比24個(gè)快拍時(shí)DOA估計(jì)旁瓣低。圖10為水池實(shí)驗(yàn)對(duì)相干信號(hào)使用MUSIC算法DOA估計(jì)的結(jié)果，由圖10可知，在小快拍情況下，MUSIC算法有偽目標(biāo)出現(xiàn)；在大快拍情況下，才可以勉強(qiáng)分開相干目標(biāo)。

圖9 水池實(shí)驗(yàn)相干源主特征空間DOA估計(jì)Fig.9 Tank experiment DOA of coherent source with main feature space

圖10 水池實(shí)驗(yàn)相干源MUSIC算法DOA估計(jì)Fig.10 Tank experiment DOA of coherent source with MUSIC

4.2 獨(dú)立信號(hào)DOA估計(jì)

水池實(shí)驗(yàn)采用主特征方法對(duì)獨(dú)立信號(hào)DOA估計(jì)結(jié)果，從圖11可知，在小快拍情況下，3個(gè)目標(biāo)都可以分辨開，這是由于在小快拍情況下主特征空間擁有所有目標(biāo)信號(hào)的信息；而在大快拍情況下，主特征空間信號(hào)信息泄漏，從圖11虛線可知，3個(gè)目標(biāo)有2個(gè)已很難分辨出。這與“2.2”節(jié)及“3.2”節(jié)仿真結(jié)論一致，進(jìn)一步驗(yàn)證了“2.2”節(jié)分析的結(jié)果。

圖11 水池實(shí)驗(yàn)獨(dú)立源主特征空間DOA估計(jì)Fig.11 Tank experiment DOA of independent source with main feature space

5 海試驗(yàn)證

海試數(shù)據(jù)采用1993-10-26/10-27，SACLANT研究中心在意大利西海岸北厄爾巴島附近進(jìn)行2天海試的數(shù)據(jù)，海況等級(jí)為2～3級(jí)，該實(shí)驗(yàn)使用的接收陣是48元的均勻線陣，采樣頻率1 kHz。參照文獻(xiàn)[15]，在335 Hz處有2個(gè)目標(biāo)，本文對(duì)時(shí)域信號(hào)分段，每1 000數(shù)據(jù)點(diǎn)做為一個(gè)快拍，傅里葉變換后，取335 Hz處的數(shù)據(jù)構(gòu)造樣本協(xié)方差矩陣。

圖12為使用本文提出方法得到的目標(biāo)DOA估計(jì)，使用了48個(gè)快拍(快拍數(shù)等于陣元數(shù))，與Qian等獲得的目標(biāo)方位結(jié)果一致，即準(zhǔn)確估計(jì)出目標(biāo)的方位。圖13為與圖12使用相同的樣本協(xié)方差矩陣，通過MUSIC算法DOA估計(jì)所得到的結(jié)果，由圖13可以發(fā)現(xiàn)其分辨出3個(gè)目標(biāo)，出現(xiàn)了偽目標(biāo)，無法準(zhǔn)確得到目標(biāo)方位。

圖12 海試主特征空間DOA估計(jì)Fig.12 Sea trial DOA with main feature space

圖13 海試MUSIC算法DOA估計(jì)Fig.13 Sea trial DOA with MUSIC

6 結(jié) 論

本文針對(duì)MUSIC算法對(duì)相干目標(biāo)分辨的失效以及在小快拍情況下傳統(tǒng)波達(dá)方向估計(jì)方法目標(biāo)分辨能力的下降，通過分析相干目標(biāo)與獨(dú)立目標(biāo)的樣本協(xié)方差矩陣的譜分離特性，提出了基于主特征空間的波達(dá)方向估計(jì)方法。

(1)該方法在深入分析相干目標(biāo)與獨(dú)立目標(biāo)信號(hào)信息在特征空間分布的基礎(chǔ)上，發(fā)現(xiàn)在小快拍情況下，目標(biāo)信息主要分布在樣本協(xié)方差矩陣的最大特征值及所對(duì)應(yīng)的特征向量上，將其定義為主特征空間；利用導(dǎo)向向量與噪聲子空間正交，導(dǎo)向向量張成的空間與信號(hào)子空間為同一空間，定義了基于MFS的目標(biāo)方位估計(jì)方法。

(2)通過仿真、水池實(shí)驗(yàn)以及海試數(shù)據(jù)驗(yàn)證表明，對(duì)于相干目標(biāo)，不管快拍數(shù)大小，本文方法都可以實(shí)現(xiàn)目標(biāo)方位的準(zhǔn)確分辨，只是快拍數(shù)較大時(shí)旁瓣較低，彌補(bǔ)了MUSIC算法對(duì)相干目標(biāo)方位估計(jì)性能下降的缺陷，且不需要先驗(yàn)信息信源個(gè)數(shù)；對(duì)于獨(dú)立目標(biāo)，在小快拍條件下也可以準(zhǔn)確分辨目標(biāo)，改善了在小快拍情況下，樣本協(xié)方差矩陣估計(jì)誤差導(dǎo)致的傳統(tǒng)DOA估計(jì)方法無法準(zhǔn)確估計(jì)目標(biāo)，甚至出現(xiàn)偽目標(biāo)的情況?？傊?，本文所提方法在不需要信源個(gè)數(shù)等先驗(yàn)信息情況下，使用小快拍樣本協(xié)方差矩陣，就能實(shí)現(xiàn)相干目標(biāo)與獨(dú)立目標(biāo)的估計(jì)，將會(huì)在UUV(Unmanned Underwater Vehicle)、魚雷等水下目標(biāo)方位估計(jì)的小快拍應(yīng)用場(chǎng)景中有廣闊的前景。

參 考 文 獻(xiàn)

[ 1 ] 王永良, 陳輝, 彭應(yīng)寧, 等. 空間譜估計(jì)理論與算法[M]. 北京： 清華大學(xué)出版社, 2004.

[ 2 ] CARLSON B D. Covariance matrix estimation errors and diagonal loading in adaptive arrays[J]. IEEE Transactions on Aerospace Electronic Systems, 1988, 24(4): 397-401.

[ 3 ] CHEN Y, WIESEL A, ELDAR Y C, et al. Shrinkage algorithms for MMSE covariance estimation[J]. IEEE Transactions on Signal Processing, 2010, 58(10): 5016-

5029.

[ 4 ] PILLAI S U, KWON B H. Forward/backward spatial smoothing techniques for coherent signal identification[J]. IEEE Transactions on Acoustics Speech & Signal Processing, 1989, 37(1): 8-15.

[ 5 ] EVANS J E, SUN D F, JOHNSON J R. Application of advanced signal processing techniques to angle of arrival estimation in ATC navigation and surveillance systems[R]. Lexington MA: MIT Lincoln Laboratory, 1982.

[ 6 ] GERSHMAN A B, BOHME J F. Improved DOA estimation via pseudorandom resampling of spatial spectrum[J]. IEEE Signal Processing Letters, 1997, 4(2): 54 - 57.

[ 7 ] VASYLYSHYN V. Removing the outliers in root-MUSIC via pseudo-noise resampling and conventional beamformer[J]. Signal Processing, 2013, 93(12): 3423-3429.

[ 8 ] QIAN C, HUANG L, SO H C. Improved unitary root-MUSIC for DOA estimation based on pseudo-noise resampling[J]. IEEE Signal Processing Letters, 2014, 21(2): 140-144.

[ 9 ] SHAGHAGHI M, VOROBYOV S A. Subspace leakage analysis and improved DOA estimation with small sample size[J]. Signal Processing IEEE Transactions on, 2015, 63(12): 3251-3265.

[10] BAI Z D, SILVERSTEIN J W. No eigenvalues outside the support of the limiting spectral distribution of large-dimensional sample covariance matrices[J]. Annals of Probability, 1998, 26(1): 316-345.

[11] BAI Z D, SILVERSTEIN J W. Exact separation of eigenvalues of large dimensional sample covariance matrices[J]. Annals of Probability, 1999, 27(3): 1536-1555.

[12] 李華, 白志東, 肖玉山. 大維隨機(jī)矩陣的漸進(jìn)特征[J]. 東北師大學(xué)報(bào)(自然科學(xué)版),2014,46 (4): 1-8.

LI Hua, BAI Zhidong, XIAO Yushan. The asymptotic properties of the large dimension random matrix[J]. Journal of Northeast Normal University (Natural Science), 2014 ,46(4): 1-8.

[14] STERGIOPOULOS S, SULLIVAN E J. Extended towed array processing by an overlap correlator[J]. The Journal of the Acoustical Society of America, 1989, 86(1): 158-171.

[15] CHAVALI V, WAGE K E, BUCK J R. Coprime processing for the Elba Island sonar data set[C]∥ 48th Asilomar Conference on Signals, Systems and Computers. California: IEEE, 2014.