李 濤，李國林2，廖輝榮

(1.海軍航空工程學(xué)院 研究生管理大隊(duì),山東 煙臺(tái) 264001;2.海軍航空工程學(xué)院 七系,山東 煙臺(tái) 264001)

1 引 言

多重信號(hào)分類(MUSIC) DOA估計(jì)算法作為主流的超分辨算法之一，自出現(xiàn)以來就一直受到關(guān)注，由于其運(yùn)算量較大，關(guān)于降低其運(yùn)算量的研究也一直在進(jìn)行。MUSIC算法的運(yùn)算量主要集中在特征值分解和空間譜的搜索兩部分，目前已有的研究也主要集中在這兩方面。采用自適應(yīng)算法來估計(jì)子空間的方法能夠減少計(jì)算量且適合子空間的更新，Oja[1]提出了基于陣列輸出數(shù)據(jù)的協(xié)方差矩陣構(gòu)建代價(jià)函數(shù)，并采用梯度算法遞歸估計(jì)子空間；SarkarT[2]對(duì)其進(jìn)行改進(jìn)，即采用梯度的瞬時(shí)估計(jì)來代替基于協(xié)方差矩陣代價(jià)函數(shù)的梯度，得到最小均方(LMS)型梯度的估計(jì)子空間，但在每一步參數(shù)矢量更新時(shí)需要對(duì)參數(shù)矢量進(jìn)行歸一化處理，引入了多個(gè)乘法和除法運(yùn)算，造成實(shí)時(shí)計(jì)算的瓶頸；Lei Xu[3]對(duì)歸一化過程進(jìn)行改進(jìn)，簡化了運(yùn)算。上述方法均從代價(jià)函數(shù)的最小化設(shè)計(jì)自適應(yīng)過程，得到的運(yùn)算量約為O(M2N)，文獻(xiàn)[4-9]則分別利用接收數(shù)據(jù)近似、陣元降維，以及基于Householoder壓縮理論等方法研究了信號(hào)子空間或噪聲子空間的快速估計(jì)算法，避開了特征值分解問題或簡化了運(yùn)算，得到的運(yùn)算量約為O(MN)；文獻(xiàn)[10-13]分別利用FFT、陣元降維和有限域搜索等方法研究了空間譜的快速搜索問題。這些方法不同程度地降低了MUSIC算法的運(yùn)算量，但同時(shí)多數(shù)存在估計(jì)精度降低的問題。本文從相干信號(hào)的自適應(yīng)對(duì)消的角度考慮噪聲子空間的估計(jì)問題，得到基于陣列輸出自適應(yīng)對(duì)消的噪聲子空間估計(jì)算法。

2 信號(hào)模型與MUSIC算法

圖1 均勻線性陣列Fig.1 Uniformlinear array

考慮一個(gè)由N個(gè)全向陣元組成的均勻線性陣列，陣列間距為d，如圖1所示。假設(shè)M個(gè)遠(yuǎn)場(chǎng)窄帶信號(hào)(M

(1)

式中，si(t)為第i個(gè)信號(hào)，λi為其中心波長，nk(t)為第k個(gè)陣元中的零均值高斯加性白噪聲。則陣列的輸出信號(hào)矢量可表示為

X(t)=[x1(t),x2(t),x3(t),…,xN(t)]T=

A(θ)S(t)+N(t)

(2)

其中：

A(θ)=[a(θ1),a(θ2),…,a(θM)]

(3)

為N×M維陣列流形矩陣，a(θi)為對(duì)應(yīng)的方向向量，且有：

(4)

S(t)=[s1(t),s2(t),s3(t),…,sM(t)]T

(5)

為M個(gè)入射信號(hào)矢量。

N(t)=[n1(t),n2(t),n3(t),…,nN(t)]T

(6)

為陣列噪聲矢量，且滿足：

E[N(t)NH(t)]=σ2I

(7)

由式(2)，陣列輸出信號(hào)自相關(guān)陣為

R=E[X(t)XH(t)]=A(θ)RsAH(θ)+σ2I

(8)

其中：

RS=E[S(t)SH(t)]

(9)

當(dāng)信號(hào)不相關(guān)時(shí)，RS為滿秩的M維對(duì)角陣。

對(duì)R進(jìn)行特征分解，由R的M個(gè)大特征值對(duì)應(yīng)的特征向量E張成M維信號(hào)子空間，而R的N-M個(gè)小特征值對(duì)應(yīng)的特征向量U張成N-M維噪聲子空間，由于：

span[A(θ)]=span[E]

(10)

因而U⊥span[A(θ)]，進(jìn)而U⊥a(θi)。采用MUSIC算法來估計(jì)信號(hào)DOA，空間譜為[14]

(11)

式中，U為噪聲子空間矩陣，a(θ)為方向向量。

由此我們發(fā)現(xiàn)MUSIC算法是由U⊥a(θi)這一本質(zhì)特征來搜索空間譜的，因此對(duì)式(2)左邊乘UH可得：

UHX(t)=UHA(θ)S(t)+UHN(t)=UHN(t)

(12)

式(12)的結(jié)果是N-M個(gè)純?cè)肼曅盘?hào)，且由于U的列相互正交，結(jié)果是N-M個(gè)正交的噪聲信號(hào)。

3 基于相干對(duì)消的噪聲子空間估計(jì)

噪聲子空間估計(jì)的問題是怎么求取U，由矩陣分解得到的噪聲子空間矩陣U為噪聲子空間的一個(gè)完備正交基底，事實(shí)上噪聲子空間的估計(jì)并不需要得到完備的正交基底，一組完備的互不相關(guān)的基底也能滿足要求。由式(12)可發(fā)現(xiàn)，用U的列向量對(duì)陣列輸出進(jìn)行加權(quán)求和，結(jié)果不含任何信號(hào)分量，換句話說，如果一個(gè)列向量對(duì)陣列輸出進(jìn)行加權(quán)求和，結(jié)果不含任何信號(hào)分量，則此列向量必定位于噪聲子空間。從這個(gè)角度出發(fā)，如果將陣列輸出進(jìn)行相互的自適應(yīng)對(duì)消，則由于各個(gè)陣元的輸出中信號(hào)之間是相干的，輸出將能夠完全對(duì)消，下面將證明這一點(diǎn)在信號(hào)源數(shù)目少于陣元數(shù)目是能夠保證的。為此首先我們看圖2的維納濾波器用于干擾對(duì)消的原理。

圖2 維納濾波器用于干擾對(duì)消原理Fig.2 Principle of Winner filter for coherent signal cancellation

(13)

式中，RX1為X1(t)的自相關(guān)陣，rX1d1為X1(t)與d1(t)的互關(guān)。X1(t)能夠?qū)ο鹍1(t)中與之相關(guān)的信號(hào)，輸出ε1(t)是d1(t)中與X1(t)不相關(guān)的期望信號(hào)。如果d1(t)中存在M個(gè)不相關(guān)的干擾，則WX1向量維數(shù)必須大于或等于M才能實(shí)現(xiàn)干擾的完全對(duì)消。

考慮陣列輸出的信號(hào)分量A(θ)S(t)，以第一個(gè)陣元作為維納對(duì)消器的參考信號(hào)為例，則：

d1(t)=[1…1]S(t)+n1(t)

(14)

X1(t)=A′(θ)S(t)+N′(t)

(15)

式中，A′(θ)為A(θ)除去第一行后的矩陣，N′(t)為N(t)除去第一個(gè)元素后的向量。由于噪聲為互不相關(guān)的，因此無法對(duì)消。由式(14)、式(15)知：

[1…1]S(t)=[1…1](A′HA′)-1A′HX0(t)=

bX1(t)

(16)

式中：

b=[1…1](A′HA′)-1A′H

(17)

當(dāng)A′(θ)列滿秩，即信號(hào)源個(gè)數(shù)小于陣元個(gè)數(shù)時(shí)式(17)存在唯一確定解。式(16)表明，d1(t)與X1(t)中的信號(hào)完全相干，因此在滿足式(17)有解的條件下將能實(shí)現(xiàn)完全對(duì)消，即：

[1 -b]A(θ)S(t)=0

(18)

由于S(t)中時(shí)間標(biāo)量的任意性，必有：

[1 -b]A(θ)=0T

(19)

由式(19)知行向量[1 -b]必定位于噪聲子空間。如果對(duì)陣列的各個(gè)陣元均進(jìn)行對(duì)消，如圖3所示，得到的N個(gè)行向量則張成整個(gè)噪聲子空間。

圖3 陣列輸出相干對(duì)消原理Fig.3 Principle of coherent array output signals cancellation

圖3中Xi(t)為陣列輸出除去第i個(gè)陣元的輸出后構(gòu)成的陣列輸出向量。

行向量[1 -bi]的求解有賴于式(13)的求解，對(duì)式(13)采用自適應(yīng)算法求解則可得到噪聲子空間的自適應(yīng)求解算法。我們以LMS算法為例來說明，易知求解式(13)的LMS算法權(quán)值更新方程為

WXi(t+1)=WXi(t)+μxi(t)εi(t)

(20)

式中,μ為自適應(yīng)步長，定義由WXi(k)生成的對(duì)角元素為零矩陣

(21)

式中省略了時(shí)間標(biāo)量，且W的第個(gè)i行向量的非零元素對(duì)應(yīng)于用于對(duì)第i個(gè)陣元進(jìn)行對(duì)消的權(quán)矢量WXi(t)，則：

ε(t+1)=X(t)-W(t)X(t)

(22)

式中,ε(t+1)為陣列輸出對(duì)消后的誤差向量，式(20)的系數(shù)更新過程的矩陣形式為

W(t+1)=W(t)+2με(t)XH(t)

(23)

其中對(duì)角元素不予計(jì)算，且設(shè)為0。由W構(gòu)造矩陣

(24)

理想情況D的各行將最終收斂于行向量[1 -bi]構(gòu)成的矩陣，即噪聲子空間。

由于D的N×N維矩陣，當(dāng)已知信號(hào)個(gè)數(shù)時(shí)，式(22)的迭代過程可以只限于W的前M列，M為信號(hào)個(gè)數(shù)。式(22)和式(23)一次迭代過程需要的計(jì)算量為2(N-1)M次乘法和2(N-1)M次加法運(yùn)算，運(yùn)算量為O[(N-1)M]。將D代替式(11)中的U，則可以進(jìn)行空間譜估計(jì)。當(dāng)信號(hào)個(gè)數(shù)未知時(shí)，由于D為N×N矩陣，搜索過程較之式(11)運(yùn)算量變大，尤其是當(dāng)信號(hào)源數(shù)較多時(shí)，此時(shí)可以對(duì)搜索過程進(jìn)行一定的優(yōu)化，例如利用D的第一行與a(θ)的乘積估計(jì)功率譜，保存功率譜峰值，再通過D的其它行來對(duì)第一步的峰值進(jìn)行驗(yàn)證，采用這種方法將大大降低搜索的計(jì)算量?？偨Y(jié)基于相干對(duì)消的噪聲子空間估計(jì)算法如下(基于LMS算法)：

(1)設(shè)置初始系數(shù)矩陣W為全零矩陣；

(2)由式(22)計(jì)算ε(t)；

(3)由式(23)更新系數(shù)矩陣W；

(4)計(jì)算ε(t)功率，判斷W是否收斂，當(dāng)收斂時(shí)進(jìn)入第5步的搜索過程，否則轉(zhuǎn)入第2步；

(5)設(shè)計(jì)搜索步長Δθ，由式(24)根據(jù)W構(gòu)建D；

(6)由D的第一行搜索功率譜，存儲(chǔ)峰值位置；

(7)由D的其它行驗(yàn)證存儲(chǔ)峰值，保留對(duì)于D所有行均形成峰值的位置作為最終的搜索結(jié)果。

如果對(duì)DOA進(jìn)行連續(xù)跟蹤，只需對(duì)W根據(jù)式(23)進(jìn)行連續(xù)更新，并進(jìn)行連續(xù)空間譜搜索即可。

4 仿真試驗(yàn)及結(jié)果分析

仿真試驗(yàn)驗(yàn)證算法4個(gè)方面的性能，試驗(yàn)1驗(yàn)證算法對(duì)信號(hào)DOA估計(jì)的準(zhǔn)確性，試驗(yàn)2驗(yàn)證算法在不同信噪比條件下的DOA估計(jì)性能，試驗(yàn)3驗(yàn)證算法對(duì)信號(hào)DOA估計(jì)的分辨力，試驗(yàn)4驗(yàn)證算法的DOA跟蹤性能，均通過與MUSIC算法的比較來驗(yàn)證。

(1)試驗(yàn)1：DOA估計(jì)

設(shè)陣元個(gè)數(shù)為16，陣元間距d=λ/2，λ為信號(hào)波長。信號(hào)采樣率為2.5 GHz，信號(hào)載頻為500 MHz。4個(gè)等功率窄帶信號(hào)分別從-45°、-30°、30°和22.5° 4個(gè)方向入射到陣列，圖4(a)給出了信噪比為10 dB的仿真結(jié)果，圖4(b)給出了信噪比為0 dB的仿真結(jié)果，仿真采用256次快拍，μ=0.001。

圖4共進(jìn)行30次蒙特卡羅仿真，為便于觀察對(duì)MUSIC算法的結(jié)果分別下移了30 dB和15 dB。結(jié)果顯示在10 dB信噪比條件下本文算法能夠很好逼近MUSIC算法性能，信噪比為0 dB時(shí)，算法性能較之MUSIC算法下降，但依然保持良好的高分辨率性能。

(a)信噪比為10 dB的結(jié)果

(b)信噪比為0 dB 的結(jié)果圖4 DOA估計(jì)的蒙特卡羅仿真結(jié)果Fig.4 Simulation results of DOA estimation by Monte Carlo

(2)試驗(yàn)2：不同比信噪比條件下的DOA估計(jì)

設(shè)陣元個(gè)數(shù)為12，2個(gè)等功率窄帶信號(hào)分別10°、25°方向入射到陣列，信噪比為-20～5 dB，其它參數(shù)不變，在每個(gè)信噪比條件下做200次蒙特卡羅仿真，用兩個(gè)信號(hào)DOA估計(jì)方差的平均作為最終的方差估計(jì)，圖5給出了不同信噪比時(shí)的DOA估計(jì)方差。

圖5 不同信噪比的DOA估計(jì)均方差Fig.5 RMSE of DOA estimation by different SNR

圖5的結(jié)果顯示，在低信噪比的情形下，采用本文算法估計(jì)的DOA的方差略大于MUSIC算法，這主要是因?yàn)長MS算法是一種次優(yōu)估計(jì)技術(shù)，存在固有的估計(jì)誤差，即額外最小均方誤差(MSE)。

(3)試驗(yàn)3：DOA估計(jì)的分辨力

設(shè)陣元個(gè)數(shù)為12，2個(gè)等功率窄帶信號(hào)分別從20°-Δθ和20°+Δθ方向入射到陣列，Δθ取0.2°到6°范圍，步進(jìn)為0.2°，信噪比為5 dB，其它參數(shù)不變，對(duì)每個(gè)Δθ采用128次快拍做200次蒙特卡羅仿真，以能夠分辨兩個(gè)信號(hào)DOA為成功標(biāo)準(zhǔn)，圖6給出了不同Δθ時(shí)的成功率。圖6的結(jié)果顯示采用本文算法估計(jì)的DOA的分辨力略小于MUSIC算法，同樣是因?yàn)楣逃泄烙?jì)誤差的結(jié)果。

圖6 角度分辨力Fig.6 Discrimination performance

(4)試驗(yàn)4：DOA跟蹤性能

(25)

圖7 DOA跟蹤性能Fig.7 Performance of DOA Tracking

圖7中的兩條曲線是空間譜峰值的軌跡。在信號(hào)源方向間隔較大時(shí)，算法能夠起到良好的跟蹤效果。其跟蹤效果與基于陣列輸出信號(hào)自相關(guān)陣更新的MUSIC算法效果相當(dāng)。在信號(hào)源交會(huì)段DOA的軌跡出現(xiàn)波動(dòng)甚至重疊，但依然保持良好的平滑性，較之MUSIC算法不會(huì)出現(xiàn)錯(cuò)誤的跟蹤結(jié)果，這主要?dú)w功于基于LMS系數(shù)矩陣W的更新是連續(xù)的，體現(xiàn)了式(23)連續(xù)迭代運(yùn)算的優(yōu)越性。

5 結(jié) 論

綜上所述，采用陣列輸出自適應(yīng)對(duì)消的噪聲子空間估計(jì)方法，能夠在不需要進(jìn)行特征分解的前提下估計(jì)噪聲子空間，大大降低了計(jì)算量，且算法的計(jì)算過程適合并行計(jì)算，便于工程實(shí)現(xiàn)。由于自適應(yīng)算法迭代計(jì)算的本質(zhì)，適合用來進(jìn)行噪聲子空間的更新估計(jì)，對(duì)于運(yùn)動(dòng)信號(hào)源來說，方位的跟蹤估計(jì)的計(jì)算量較之特征分解的方法降低是非常可觀的。工程實(shí)用時(shí)有以下幾點(diǎn)考慮：

(1)當(dāng)信號(hào)源數(shù)較少時(shí)，對(duì)消系數(shù)矩陣W可以進(jìn)行簡化，可以僅保留W的一定量的斜對(duì)角元素，而將其它元素設(shè)為零，這樣用來對(duì)消的陣列輸出減少，從而減少計(jì)算量；

(2)關(guān)于陣元數(shù)。仿真過程發(fā)現(xiàn)，當(dāng)陣元數(shù)較少時(shí)，例如8個(gè)陣元，算法性能在信號(hào)源數(shù)增加時(shí)下降較快，但當(dāng)陣元數(shù)較大時(shí)，能很好逼近MUSIC算法的分辨率性能；

(3)關(guān)于自適應(yīng)過程。論文的仿真采用的是LMS算法，也可以采用RLS算法，可以預(yù)見RLS會(huì)有更快的收斂速度。

參考文獻(xiàn)：

[1] Erkki Oja. On Stochastic Approximation of the Eigenvectors and Eigenvalues of the Expectation of a Random Matrix[J]. Journal of Mathematic Analysis and Applications，1985,33(1):69-84.

[2] Sarkar T K,Dianat S A,Chen H,et al. Adaptive Spectral Estimation by the Conjugate Gradient Method[J]. IEEE Transactions on ASSP, 1986,34(2):272-284.

[3] LEI Xu, Erkki Oja, Ching Y S. Modified Hebbian Learning for Curve and Surface Fitting[J].Neural Networks,1992,5(3):441-457.

[4] Gierull C H. Angle estimation for small sample size with fast eigenvector-free subspace method[J]. IEEE Proceedings-Radar Sonar Navigation,1999,146(3):126-132.

[5] Mohammed A Hasan, Jawad A K Hasan. Fast algorithms for signal subspace estimation with applications to DOA estimation[C]//Proceedings of the 1999 IEEE International Symposium on Circuits and Systems.Orlando, FL, USA:IEEE,1999:223-226.

[6] Mohammed A Hasan. Fast rational approximation algorithms of signal and noise[C]//Proceedings of 2001 Sixth International Symposium on Signal Processing and its Applications. Kuala Lumpur,Malaysia:IEEE,2001: 124-127.

[7] Strobach Peter. The Householder compressor theorem and its application in subspace tracking[J].Signal Processing, 2009, 89(5): 857-875.

[8] 黃磊,吳順君,張林讓,等.快速子空間分解方法及其維數(shù)的快速估計(jì)[J]. 電子學(xué)報(bào),2005,33(6):977-981.

HUANG Lei, WU Shun-jun, ZHANG Lin-rang,et al. A fast method for subspace decomposition and its dimension estimation[J]. ACTA ELECTRONICA SINICA,2005,33(6):977-981.(in Chinese)

[9] Stoica P.Maximum-Likelihood DOA Estimation by Data-Supported Grid Search[J]. IEEE Signal Processing Letters, 1999,10(6):273-275.

[10] Paine A S, Qineti Q Malvern. Fast MUSIC for large 2-D element digitized phased array radar[C]//Proceedings of the International Radar Conference.Adelaide Australia:IEEE,2003: 200-205.

[11] CUI Wei-wei, CAO Zhi-gang. Fast Source Location Method Using Anti-reverberant Searching Space Pre-estimation[C]//Proceedings of 2006 IEEE Region 10 Conference(TENCON 2006). Hongkong:IEEE,2006:1-4.

[12] 蔣毅,古天祥. 基于有限域搜索的MUSIC法快速頻率估計(jì)[J].儀器儀表學(xué)報(bào),2006,27(11):1526-1528.

JIANG Yi, GU Tian-xiang. Quick frequencyestimation based on MUSIC algorithm[J]. Chinese Journal of Scientific Instrument, 2006,27(11):1525-1528.(in Chinese)

[13] 曾浩,張迎輝,馮文江.快速子空間譜峰搜索方法[J].計(jì)算機(jī)應(yīng)用,2009,29(9):2546-2547.

ZENG Hao,ZHANG Ying-hui,FENG Wen-jiang.New peak searching method in subspace spectrum[J].Journal of Computer Applications,2009,29(9):2564-2547.(in Chinese)

[14] Schmidt R O. Multiple emitter location and signal parameter estimation[J]. IEEE Transactions on Antennas and Propagation, 1986,34(3):276-280.