張 超，黎仁剛，顧 軍

(中國船舶重工集團(tuán)公司第七二三研究所，江蘇 揚(yáng)州 225001)

基于多級(jí)維納濾波的ESPRIT算法

張 超，黎仁剛，顧 軍

(中國船舶重工集團(tuán)公司第七二三研究所，江蘇 揚(yáng)州 225001)

因奇異值分解引入了龐大的計(jì)算量，子空間分解類的測向算法在工程中難以實(shí)現(xiàn)實(shí)時(shí)處理。多級(jí)維納濾波(MSWF)算法，避免了奇異值分解運(yùn)算，有效地減少了運(yùn)算量。通過對(duì)多級(jí)維納濾波算法進(jìn)行改進(jìn)，提出了一種信源數(shù)目的估計(jì)方法。與傳統(tǒng)多級(jí)維納濾波算法相比，該方法提高了測向精度。通過仿真證明了該算法的有效性。

陣列信號(hào)處理；多級(jí)維納濾波；信源數(shù)目估計(jì)；測向精度

0 引 言

波達(dá)方向估計(jì)(DOA)在通信、雷達(dá)、電子對(duì)抗、導(dǎo)航等領(lǐng)域中有著舉足輕重的地位[1]。隨著信息技術(shù)的不斷發(fā)展，電磁環(huán)境變得越來越復(fù)雜，傳統(tǒng)的測向技術(shù)無法滿足當(dāng)下的需求，以MUSIC[2]和ESPRIT[3]算法為代表的超分辨測向算法應(yīng)運(yùn)而生，并受到廣大研究人員的青睞。但其算法需要進(jìn)行子空間分解，計(jì)算量龐大，不易實(shí)現(xiàn)實(shí)時(shí)性。為了減少算法的運(yùn)算量，提高測向速度，許多學(xué)者也提出了一些方案，文獻(xiàn)[4]通過協(xié)方差矩陣實(shí)值化處理，將復(fù)矩陣分解轉(zhuǎn)換成矩陣分解，文獻(xiàn)[5]通過將陣元空間轉(zhuǎn)換到波束空間，降低矩陣維度。

以上所提出的算法都需要進(jìn)行子空間分解，文獻(xiàn)[6]提到了傳播算子(PA)，此方法利用線性算子從陣列數(shù)據(jù)中提取信號(hào)子空間與噪聲子空間，以矩陣相乘代替特征值分解，減小了運(yùn)算量，但此方法無法對(duì)信源數(shù)目進(jìn)行估計(jì)，需要預(yù)先知道信源數(shù)目才能構(gòu)造分塊矩陣。Goldstein等人在1998年提出了基于多級(jí)維納濾波的方法[7]，此方法在自適應(yīng)波束形成、自適應(yīng)干擾抑制等領(lǐng)域廣泛應(yīng)用。文獻(xiàn)[8]將多級(jí)維納濾波技術(shù)應(yīng)用在子空間分解類算法中，避免了特征值分解，大大降低了運(yùn)算量。文獻(xiàn)[9]通過協(xié)方差矩陣來初始化數(shù)據(jù)，不需要訓(xùn)練信號(hào)就能進(jìn)行測向，但在低信噪比、小快拍數(shù)的情況下，誤差變大甚至可能失效。文獻(xiàn)[10]通過對(duì)初始化數(shù)據(jù)進(jìn)行優(yōu)化，使其在低信噪比、小快拍數(shù)的條件下，估計(jì)性能得到改善。文獻(xiàn)[11]提出了一種信號(hào)源數(shù)目與波達(dá)方向聯(lián)合估計(jì)的方法，但此方法需要進(jìn)行M(陣元數(shù))次迭代，在陣元數(shù)較多的情況下，運(yùn)算量也較為龐大。本文通過改進(jìn)多級(jí)維納濾波算法，優(yōu)化初始矩陣，提出一種高效的信源數(shù)目估計(jì)方案，在陣元數(shù)較多的條件下，與文獻(xiàn)[11]中方法相比，有效降低了計(jì)算量。

1 改進(jìn)算法

1.1 信號(hào)模型

對(duì)于均勻圓陣，N個(gè)全向天線均勻分布在圓周上，如圖1所示，D個(gè)不相關(guān)的遠(yuǎn)場窄帶信號(hào)入射，方位角與俯仰角分別為θi和βi(i=1,2,…,D)。陣元的輸出數(shù)據(jù)模型為：

X=As+n

(1)

A=[α(θ1,β1),…,α(θD,βD)]

(2)

α(θ,β)=[ej2πcosθsinβ,ej2πcos(θ-φ1)sinβ,…, ej2πcos(θ-φN-1)sinβ]T

(3)

φi=2π(i-1)/N，i=1,2,…,N

(4)

式中:A為均勻圓陣的方向矢量；s為入射信號(hào)；n為陣列上的噪聲，均值為0，方差為σ2。

噪聲相互獨(dú)立，且與信號(hào)互不相關(guān)。

對(duì)于均勻圓陣而言，其方向矢量的表達(dá)式相對(duì)較復(fù)雜，且不具備均勻線陣的范得蒙矩陣形式。因此，首先要將其進(jìn)行波束空間轉(zhuǎn)換變成虛擬線陣。

設(shè)空間離散傅里葉變換(DFT)矩陣:

(5)

(6)

(7)

(8)

定義變換矩陣C:

C=diag{j-M,…,j-1,j0,j-1,…,j-M}

(9)

FH=VHCWH

(10)

根據(jù)X=As+n，求出X的協(xié)方差矩陣RX，再將矩陣RX進(jìn)行變換，得到:

RF=FHRXF

(11)

式中：RF為在波束空間下對(duì)應(yīng)的協(xié)方差矩陣，取RF的實(shí)部代替RF。

1.2 多級(jí)維納濾波優(yōu)化

維納濾波技術(shù)是在最小均方誤差(MMSE)準(zhǔn)則下的一種線性濾波，在已知信號(hào)與噪聲相關(guān)函數(shù)的條件下，通過Winner-Holf方程RxWx=r=As+nxd對(duì)信號(hào)進(jìn)行最優(yōu)濾波。Goldstein等人提出的多級(jí)維納濾波技術(shù)的思想是通過對(duì)觀測數(shù)據(jù)進(jìn)行多次正交投影，每次投影得到2個(gè)子空間，其中一個(gè)子空間平行于參考信號(hào)與上次觀測數(shù)據(jù)的互相關(guān)矢量，另一個(gè)子空間正交于這個(gè)子空間，然后再對(duì)這個(gè)正交于互相關(guān)矢量的子空間用同樣的方法繼續(xù)進(jìn)行正交投影。其結(jié)構(gòu)如圖1所示。

多級(jí)維納濾波算法需要參考信號(hào)di(k)與觀測數(shù)據(jù)Xi(k)，而實(shí)際工程中，并不一定能夠獲得，因此采用任意單元的數(shù)據(jù)來代替參考信號(hào)。信號(hào)之間相互獨(dú)立，假設(shè)為高斯白噪聲。觀測信號(hào)每次投影，分別得到與觀測數(shù)據(jù)平行的分量和正交的分量，這些代表的就是信號(hào)子空間與噪聲子空間。具體步驟如下所示：

第1步，初始化觀測數(shù)據(jù)和陣列接收數(shù)據(jù)：

(12)

(13)

式中：R為協(xié)方差矩陣。

第2步，求互相關(guān)矢量hi(i=1,2,3，…,D)：

(14)

(15)

(16)

(17)

第3步，估計(jì)信號(hào)源數(shù)目，確定信號(hào)子空間：

(18)

(19)

式中：S代表信號(hào)子空間；G代表噪聲子空間。

數(shù)據(jù)初始化使用協(xié)方差矩陣來代替，在文獻(xiàn)[10]中提出的初始數(shù)據(jù)優(yōu)化方法，通過將協(xié)方差矩陣行取平均，得到：

(20)

此方法應(yīng)用于均勻線陣，在低信噪比、小快拍數(shù)條件下，適當(dāng)提高了算法的估計(jì)性能。

本文算法是基于均勻圓陣的，經(jīng)過波束空間映射到虛擬均勻線陣，模式數(shù)為2M+1。本文提出的優(yōu)化方法是選用協(xié)方差矩陣第M行初始化數(shù)據(jù)：

(21)

針對(duì)均勻圓陣，本文所提出的優(yōu)化方法比起文獻(xiàn)[10]中提出的方法，在低信噪比、小快拍數(shù)下的精度更高，下文仿真實(shí)驗(yàn)可見。

1.3 信號(hào)源數(shù)目估計(jì)

在常規(guī)的子空間分解測向算法中，信號(hào)源數(shù)目的估計(jì)一般采用Akaike信息論(AIC)準(zhǔn)則和最小描述長度(MDL)準(zhǔn)則，這2個(gè)準(zhǔn)則具有相同的似然函數(shù)，只是懲罰函數(shù)不一樣。但這2個(gè)準(zhǔn)則的使用前提需要對(duì)陣列協(xié)方差矩陣進(jìn)行特征值分解。而在多級(jí)維納濾波算法中，通過正交投影取代了特征值分解。陣列協(xié)方差矩陣分解所得到的特征值分別代表信號(hào)或噪聲功率的大小。與特征值分解類似，通過多級(jí)維納濾波算法得到的正交投影信號(hào)d1(k),d2(k),…,dM(k)，可以用σi來代表其功率：

(22)

表1 不同信噪比下σi與特征值λi的比較

從表1可以看出，多級(jí)維納濾波算法在低信噪比時(shí)，信號(hào)能量會(huì)泄露到噪聲子空間中，導(dǎo)致無法判斷信號(hào)子空間。文獻(xiàn)[7]中提出了一種信號(hào)源數(shù)目的估計(jì)方法，判斷信號(hào)子空間。但該方法在陣元數(shù)比較多的時(shí)候，計(jì)算量較大。本文提出一種估計(jì)方法，適合多陣元數(shù)條件下的信號(hào)源數(shù)目估計(jì)。

第1步：初始化后進(jìn)行4次投影得到正交投影信號(hào)，求出代表其功率的σi。

第2步：構(gòu)造判決函數(shù)，判斷信號(hào)子空間需要濾波的級(jí)數(shù)：

fori=1:M-4

δi=(σi+1+2σi+2)/3

Δi=(σi-δi)/δi/2

ifΔi<1

ThenD=i

elseifσi<1

ThenD=i

endif

endfor

根據(jù)判決，得到信號(hào)源數(shù)目為i。

第3步：構(gòu)造信號(hào)子空間S=[h1,h2,…,hD]。

1.4ESPRIT算法

經(jīng)過多級(jí)維納濾波正交投影，得到的信號(hào)子空間S=[h1,h2,…,hD]相似于B：

(23)

式中：Jm(ζ)為第一類貝塞爾函數(shù)，m為階數(shù)，ζ為變量。

根據(jù)貝塞爾函數(shù)的性質(zhì)可得：

ΓZ0=Z-1Φ+Z1Φ*

(24)

Γ=(λ/πr)diag{-(M-1),…,0,…,M-1}

(25)

Φ=diag(sinβ1ejθ1,sinβ2ejθ2,…,sinβDejθD)

(26)

式中：Z-1，Z0和Z1分別為信號(hào)子空間Z前2M-1行，中2M-1行和后2M-1行。

最后求解得到Φ，即可估計(jì)出信號(hào)源方位角與俯仰角。

2 仿真與分析

本文采用陣元數(shù)為16的均勻圓陣，模式數(shù)2M+1為9，2個(gè)不相干的信號(hào)源入射角分別為(75°,20°)與(60°,55°)，對(duì)比常規(guī)波束空間ESPRIT算法(算法1)、文獻(xiàn)[10]中所提出的算法(算法2)以及本文算法(算法3)3種算法的性能。定義聯(lián)合均方根誤差(RMSE)為：

(27)

式中：P為蒙特卡洛試驗(yàn)次數(shù)；θi與βi分別為第i次試驗(yàn)的方位角與俯仰角的估計(jì)值；θ與β分別為方位角與俯仰角的真實(shí)值。

試驗(yàn)1：選擇信噪比為4dB與6dB，在不同的快拍數(shù)下，進(jìn)行200次蒙特卡洛試驗(yàn)，仿真結(jié)果如圖3所示。

從圖3可以看出，小快拍數(shù)的情況下，本文算法的估計(jì)性能優(yōu)于常規(guī)的多級(jí)維納濾波算法，接近協(xié)方差矩陣分解算法的估計(jì)性能。隨著快拍數(shù)的增加，各算法的估計(jì)性能也隨之提高。

試驗(yàn)2：選擇快拍數(shù)為40次與80次，在不同的信噪比下，進(jìn)行200次蒙特卡洛試驗(yàn)，仿真結(jié)果如圖4所示。

從圖4可以看出，在低信噪比的情況下，本文算法的估計(jì)性能優(yōu)于常規(guī)的多級(jí)維納濾波算法。隨著信噪比的提高，各算法的估計(jì)性能也隨之提高，并且均方根誤差趨于一致。

試驗(yàn)3：對(duì)于本文提出的信號(hào)源數(shù)目估計(jì)方法，分別在低信噪比、小快拍數(shù)下進(jìn)行200次蒙特卡洛試驗(yàn)，觀察其成功概率。

從圖5、圖6中可以看出，在低信噪比時(shí)隨著快拍數(shù)增加和信噪比的提高，信號(hào)源數(shù)目估計(jì)的成功概率都將提高。本文提出的信號(hào)源數(shù)目估計(jì)方法中引入了σi<1這一判決條件，可避免在大信噪比時(shí)噪聲功率偏小，導(dǎo)致噪聲功率之間的相對(duì)偏差較大，引起估計(jì)錯(cuò)誤。

傳統(tǒng)的多級(jí)維納濾波算法的運(yùn)算量為O(M3)，而本文所提出的信號(hào)源數(shù)目估計(jì)方法的運(yùn)算量為O[M2(K+3)]，其中M為陣元數(shù)目，K為信號(hào)源數(shù)目。相比文獻(xiàn)[11]中的信號(hào)源數(shù)目估計(jì)方法，本文方法在陣元數(shù)較多而信號(hào)源數(shù)目較少的情況下，能有效減少運(yùn)算量。

3 結(jié)束語

根據(jù)多級(jí)維納濾波的特點(diǎn)，本文提出了一種基于多級(jí)維納濾波的改進(jìn)算法，用協(xié)方差矩陣代替參考信號(hào)，并通過優(yōu)化初始數(shù)據(jù)來提高估計(jì)精度。與常規(guī)的子空間分解算法相比，其避免了特征值分解，大大減小了運(yùn)算量。此外，在多級(jí)維納濾波投影過程中，低信噪比條件下存在子空間泄露問題。本文提出了一種信號(hào)源數(shù)目的估計(jì)方法，能夠在低信噪比時(shí)有效估計(jì)出信號(hào)源的數(shù)目。尤其是在大陣元數(shù)、信號(hào)源數(shù)目不多的情況下，本文方法比常規(guī)多級(jí)維納濾波算法運(yùn)算量更低，更適用于工程實(shí)現(xiàn)。

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

[2]SCHMIDTRO.Multipleemitterlocationandsignalparameterestimation[J].IEEETransactionsonAntennasandPropagation,1986,34(3):276-280.

[3]ROYR,PAULRAJA,KAILATHT.ESPRIT——asubspacerotationapprochtoeatimationofParameterofcisoidsinnoise[J].IEEETransactionsonAcoustics,SpeechandSignalProcessing,1986,34(5):1340- 1342.

[4] 孫繼紅,楊科,張瑾.一種實(shí)值MUSIC算法[J].空間電子技術(shù),2006(3):36-38.

[5] 李紹濱,趙宜楠,胡航.均勻圓陣波束空間二維DOA算法[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào),2004,36(6):796- 798.

[6]MARCOSS,BEBIDIRM.Onahighresolutionarrayprocessingmethodnon-basedoneigenanalysisApproach[C]//ETTERDM.InternationalConferenceonAcoustics,speechandSignalProcessing.NewMexico,USA:IEEEPress,1990:2955-2958.

[7]GOLDSTAINJS,REEDIS,ScharfLL.Amultistagerepresentionofthewinnerfilterbasedonorthogonalprojection[J].IEEETransactionsonInformationTheory,1998,44(7):2943-2959.

[8] 黃志忠,李萬玉,強(qiáng)勇.多級(jí)維納濾波在子空間類測向算法中的應(yīng)用[J].火控雷達(dá)技術(shù),2009(1):11-14.

[9] 于紅旗,黃知濤,周一宇,徐欣.一種不需要特征值分解的ESPRIT方法[J].信號(hào)處理,2008(3):357-360.

[10]李堰,宋愛民,劉劍.多級(jí)維納濾波測向算法的參考信號(hào)優(yōu)化方法[J].電訊技術(shù),2010(12):75-80.

[11]包志強(qiáng),吳順君,張林讓.一種信源個(gè)數(shù)與波達(dá)方向聯(lián)合估計(jì)的新算法[J].電子學(xué)報(bào),2006(12):2170- 2174.

ESPRIT Algorithm Baesd on Multi-stage Winner Filtering

ZHANG Chao,LI Ren-gang,GU Jun

(The 723 Institute of CSIC,Yangzhou 225001,China)

Because of the enormous calculating quantity introduced by singular value decomposition,the direction finding algorithm based on subspace decomposition is difficult to be realized in engineering in real time.The multi-stage Wiener filtering (MSWF) algorithm avoids the singular value decomposition calculation,reduces the calculating quantity effectively.This paper puts forward an estimation method of source number by improving MSWF algorithm.The algorithm improves the direction-finding accuracy compared with traditional MSWF algorithms.The validity of this algorithm is validated through simulation.

array signal processing;multi-stage Wiener filtering;estimation of source number;direction-finding accuracy

2016-09-01

TN971.1

A

CN32-1413(2017)01-0074-05

10.16426/j.cnki.jcdzdk.2017.01.016