佘黎煌, 劉平凡, 張 石, 許方晗

(東北大學 計算機科學與工程學院, 遼寧 沈陽 110169)

陣列信號處理技術包括DOA估計技術及波束形成技術,廣泛應用于空間信號探測、雷達及水聲信號探測等諸多領域[1].其中DOA估計技術作為獲得信號波達方向角度的重要技術,如何使DOA估計算法獲得更高的估計精度和更低的時間復雜度一直是科研人員的主要研究方向.

在眾多DOA估計算法中,Root-MUSIC算法以其較高的估計精度、較低的復雜度受到眾多科研人員的重視[2],大量基于Root-MUSIC算法的改進算法被提出.傳統(tǒng)Root-MUSIC算法存在的最大問題在于隨著陣列規(guī)模的增大,算法的運算負擔成倍增加.在較低陣列規(guī)模時,多項式求根使算法運行效率較高,而當陣列規(guī)模增大后,其求根多項式階數(shù)快速升高,算法運行耗時甚至超過譜峰搜索類算法[3].目前眾多基于Root-MUSIC算法的改進算法都著重針對求根多項式進行降維操作,通過較低的數(shù)據(jù)規(guī)模實現(xiàn)DOA估計.Ren等[4]提出了一種快速Root-MUSIC算法,通過將噪聲子空間進行拆分,構造新的數(shù)據(jù)向量,使得求根多項式階數(shù)降低至僅為信號源個數(shù),大大降低了算法運行復雜度.由于噪聲子空間信息利用率較低,算法精度損失嚴重,王新賀等[5]在Ren等[4]拆分噪聲子空間矩陣的基礎上重構數(shù)據(jù)向量,充分利用噪聲子空間,算法精度有了大幅度提升,但運算復雜度也隨之增大,同時算法在陣列規(guī)模較大時完全失效,不具備魯棒性.Yan等[6]將均勻線陣拆分為兩個相等子陣列,利用兩個子陣列采集的數(shù)據(jù)構造3個無噪聲互相關矩陣,用處理后的新數(shù)據(jù)矩陣進行DOA估計,并將此DOA估計數(shù)據(jù)作為參考值,重復以上方法，進行二次估計得到最終DOA估計值[7].算法通過拆分子陣構造新互相關矩陣的方法避免了特征分解運算,降低了求根運算的復雜度,但在低陣列規(guī)模時算法精度損失嚴重.探測陣列必須為偶數(shù),可估計信源個數(shù)較少,魯棒性較差.本文在Toeplitz矩陣重構及求根多項式降階的基礎上,提出了一種改進的Root-MUSIC算法,算法依據(jù)數(shù)據(jù)觀測矩陣首行重構Toeplitz矩陣.特征值分解后得到噪聲子空間,并將噪聲子空間翻轉拆分,重構新的求根多項式向量,進而通過求根方法得到DOA估計值.算法通過重構Toeplitz矩陣使得新自相關矩陣具備Hermitian性質,相比于傳統(tǒng)自相關矩陣特征值分解算法,本文算法擁有更高的DOA估計精度.噪聲子空間的翻轉拆分重構了新的低階求根多項式,有效降低了算法的時間復雜度.

1 陣列信號模型

設存在M個滿足遠場條件的各向同性陣元所組成的均勻線陣模型,同時有q個遠場窄帶目標信號源,其入射角為θk(k=1,2,…,q),并且滿足M>q.選擇首個陣元作為參考陣,陣元間距設為d,則此時陣列響應向量為

(1)

第m個陣元的輸出為

(2)

式中:sk(t)為第k個由天線陣列接收到的信源信號;λ為信號的波長;q為信號源數(shù)目;nm(t)為陣元接收的高斯白噪聲,方差為σ2.此時的陣列接收數(shù)據(jù)可以表示為向量形式:

X(t)=AS(t)+N(t) .

(3)

式中:

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

(4)

S(t)=[s1(t),s2(t),…,sq(t)]T;

(5)

N(t)=[n1(t),n2(t),…,nM(t)]T;

(6)

X(t)=[x1(t),x2(t),…,xM(t)]T.

(7)

通過式(3)可以得到理想觀測數(shù)據(jù)的協(xié)方差矩陣:

R=E[X(t)XH(t)]=ARsAH+σ2IM.

(8)

式中:Rs=E[S(t)SH(t)]為信號的協(xié)方差矩陣;IM為M階的單位矩陣.

2 本文方法描述

2.1 Toeplitz矩陣重構算法

由于硬件設備采樣快拍數(shù)的限制,由數(shù)據(jù)觀測矩陣獲得的自相關矩陣并不能滿足理想的Hermitian性質,僅能保證對角占優(yōu),因此存在估計誤差.Toeplitz矩陣重構算法通過將數(shù)據(jù)觀測矩陣的第一行即首個陣元獲得的數(shù)據(jù)作為參考向量,構造新的自相關矩陣使得新矩陣具備Hermitian性質.相較于原始自相關矩陣更加逼近理想值,DOA估計精度更高,同時不會產(chǎn)生陣列孔徑損失,估計信源個數(shù)不會受到影響[8].

由式(2)可以得到數(shù)據(jù)觀測矩陣第一行即首個陣元的輸出為

(9)

在實際情況下,由于采樣快拍數(shù)有限,因此天線陣列所接收到的信源信號及高斯白噪聲均表現(xiàn)為由若干采樣值所組成的向量形式,向量維度即采樣快拍數(shù),這里將采樣快拍數(shù)設為T,則有

x1=[x1(1),x1(2),…,x1(T)] .

(10)

同理,第k個陣元的接收數(shù)據(jù)向量表示為

xk=[xk(1),xk(2),…,xk(T)] .

(11)

將式(10)和式(11)作為參考向量,定義

(12)

式中,A(1)和A(k)分別代表矩陣A的第一行和第k行.隨著k的變化,可以得到包含全部信源信息的相關矢量[rtp(0),rtp(1),…,rtp(M-1)],由M個相關函數(shù)構成的矩陣:

Rtp=

(13)

2.2 Root-MUSIC算法

譜峰搜索算法通過全范圍角度的遍歷來搜索DOA估計值,在陣元數(shù)較少時會導致算法運行耗時增大.而求根MUSIC算法(Root-MUSIC)將這一過程省去,替代為簡化的多項式求根過程.

通過式(1)可以得到陣列響應向量,這里定義:

(14)

式(1)的陣列響應向量可以改寫為

(15)

式(4)陣列流型矩陣可以表示為

A=[z1,z2,…,zq] .

(16)

根據(jù)式(8)可以得到自相關矩陣R,對R做特征分解,較大的q個特征值對應特征向量構成信號子空間Us,較小的M-q個特征值構成噪聲子空間Un.依據(jù)MUSIC算法相關理論,信號子空間與噪聲子空間相互正交,則噪聲子空間也必與信號子空間的基即陣列流型矩陣A正交,滿足

AHUn=Oq×(M-q).

(17)

式中，O為q×(M-q)階零矩陣.根據(jù)式(17)定義噪聲子空間Un:

Un=[un1,un2,…,un(M-q)] .

(18)

此時式(17)可以表示為

(19)

通過觀察式(19),可將陣列響應方向向量定義為一般形式:

z=[1,z-1,…,z-(M-1)]T.

(20)

根據(jù)式(19)及式(20)構造多項式向量p(z):

p(z)=[p1(z),p2(z),…,pM-q(z)]T.

(21)

式中:

pi(z)=zHuni=β0,i+β1,iz+…+βM-1,izM-1,

1≤i≤M-q.

(22)

式中，uni=[β0,i,β1,i,…,βM-1,i]T為噪聲子空間列向量.

由式(21)可知,M-q個多項式pi(z)的階數(shù)均為M-1,且有q個公共根z1,z2,…,zq.

根據(jù)文獻[7]構造多項式:

froot(z)=pH(z)p(z) .

(23)

求得式(23)的根即可獲得有關信號源到達角的信息.因多項式存在z*項,這會使得求根過程變得復雜,因此對式(23)修正:

froot(z)=zM-1pT(z-1)p(z) .

(24)

在式(24)中,多項式階數(shù)為2M-2,即有M-1對根,且每對根為相互共軛關系.在這些根中有q個根z1,z2,…,zq剛好分布在單位圓上.

在實際情況下因采樣快拍數(shù)的限制無法獲得理想的自相關矩陣,因此實際得到的近似協(xié)方差矩陣存在一定程度的誤差,這時只需要對式(24)求近似的q個根即可(最靠近單位圓).在均勻直線陣下的DOA估計角度計算式為

(25)

2.3 本文提出的低復雜度求根運算方法

在Root-MUSIC算法中,時間復雜度最高的運算為多項式求根,求根運算的時間復雜度為O(2M-1)3,可知,隨陣元規(guī)模的擴大,其時間復雜度迅速增加,實驗也證實了這一點.當陣元規(guī)模達到32以上時,其耗時已經(jīng)超過了譜峰搜索算法,因此本節(jié)將提出一種降低求根多項式階數(shù)的低復雜度求根運算方法.

根據(jù)2.2節(jié)可知,設存在多項式p1(z),p2(z),…,pM-q(z)且同時擁有q個公共根z1,z2,…,zq,則它們之間存在線性組合Q(z):

(26)

接下來證明式(26)的合理性,若r為多項式p1(z),p2(z),…,pM-q(z)的公共根,則r為多項式p1(z),p2(z),…,pM-q(z)任意線性組合的根.

設存在c1,c2,…,cM-q使得Q(z)=∑cipi(x),則將公共根r代入Q(z)可得

Q(r)=∑cipi(r)=∑ci×0=0 .

(27)

通過式(27)證明了式(26)的合理性,構造多項式p1(z),p2(z),…,pM-q(z)的線性組合Q(z),對應系數(shù)設為α1,α2,…,αM-q,令αM-q=1,可以得到如下表達式:

Q(z)=α1p1(z)+…+αM-q-1pM-q-1(z)+

pM-q(z) .

(28)

由式(22)可知：

αipi(z)=αi(βM-1,izM-1+βM-2,izM-2+…+

β0,i).

因此,Q(z)可以被重寫為

Q(z)=(α1βM-1,1+…+αM-q-1βM-1,M-q-1+

βM-1,M-q)zM-1+(α1βM-2,1+…+αM-q-1×

βM-2,M-q-1+βM-2,M-q)zM-2+(α1β1,1+…+

αM-q-1β1,M-q-1+β1,M-q)z+(α1β0,1+…+

αM-q-1β0,M-q-1+β0,M-q) .

(29)

由式(27)可知,將公共根r代入Q(z)的前M-q-1階,即令式(29)中zM-1,zM-2,…,zq+1對應階系數(shù)為零,得到如下方程組:

(30)

式(30)可以寫成如下的矩陣形式:

(31)

將式(31)簡化整理為

x=-H-1h.

(32)

通過式(32)求得x,得到了對應式(26)中的系數(shù)αi,1≤i≤M-q-1.要得到x,就必須得到矩陣H及向量h.

在2.2節(jié)的討論中,對自相關矩陣R進行特征值分解得到噪聲子空間Un,其列向量uni=[β0,i,β1,i,…,βM-1,i]T,則將噪聲子空間列向量翻轉可得：

(33)

(34)

(35)

依據(jù)式(32)及式(35)計算得到線性組合Q(z)的系數(shù)向量x,并將式(21)代入式(25)得到一個新的求根多項式:

(36)

式中:αM-q=1.通過式(36)所構造的新求根多項式來替代式(24)完成求根運算,尋找最靠近單位圓的q個根,最終通過式(25)求得DOA估計角度.

本文算法的整體實現(xiàn)過程如下:

1) 根據(jù)式(12)及式(13)構造近似數(shù)據(jù)觀測矩陣下的Hermitian Toeplitz矩陣Rtp;

2) 對矩陣Rtp做特征值分解,得到噪聲子空間矩陣Un;

3) 依據(jù)式(33)～式(35),對Un翻轉拆分獲得矩陣H及向量h;

4) 通過計算式(32)得到系數(shù)向量x,并根據(jù)式(36)構造新的求根多項式Q(z);

5) 完成求根運算,并尋找最靠近單位圓的q個根,通過式(25)求得DOA估計角度.

3 實驗驗證與比較

針對改進算法的DOA估計性能及算法運行耗時進行仿真測試.為驗證本文算法的各項性能,仿真實驗分別對比了文獻[4]的快速Root-MUSIC算法、文獻[6]的無需特征分解計算的兩步Root-MUSIC算法和經(jīng)典Root-MUSIC算法;同時,對幾種算法的時間復雜度進行了分析.

均方根誤差計算式為

(37)

3.1 DOA估計精度實驗

實驗1:兩個獨立入射信源角度分別為-10°和10°;陣列為8陣元均勻直線陣列,陣列間距為λ/2;采樣快拍數(shù)為200,信噪比以2 dB間隔從-10 dB變化至10 dB，進行1 000次蒙特卡洛實驗，計算均方根誤差.不同信噪比下算法的DOA估計精度對比如圖1所示.

圖1 不同信噪比下算法的DOA估計精度對比

由圖1可知,隨著信噪比的增加,4種算法的均方根誤差均有所下降,而本文算法的估計誤差最接近經(jīng)典Root-MUSIC算法,低于文獻[4]及文獻[6]所提算法;同時可以發(fā)現(xiàn),本文算法及經(jīng)典算法更加接近算法的無偏估計量下界(克拉美羅下界)CRLB,擁有更高的估計精度.

實驗2: 獨立入射信源數(shù)量為3個,分別為-10°,10°和20°,其余實驗條件與實驗1一致,同樣使信噪比以2 dB間隔從-10 dB變化至10 dB進行1 000次蒙特卡洛實驗,計算均方根誤差.

不同信噪比下算法的DOA估計精度對比如圖2所示.

由圖2可知,當入射信源增加為3個時,由于陣列規(guī)模的限制,文獻[4]和文獻[6]的算法已經(jīng)基本失效,而本文算法及經(jīng)典Root-MUSIC算法在低信噪比環(huán)境下效果較差,隨信噪比的增加,均方根誤差都有較大幅度的下降.受限于陣列規(guī)模,相較于2個入射信源的情況本文算法及經(jīng)典Root-MUSIC算法的DOA估計性能均有所下降.

圖2 不同信噪比下算法的DOA估計精度對比

實驗3: 2個獨立入射信源角度為-10°和10°,陣列仍采用8陣元均勻直線陣列,陣列間距為λ/2,信噪比固定為5 dB.采樣快拍數(shù)以50為間隔從100變化至600進行1 000次蒙特卡洛實驗,計算均方根誤差.不同快拍數(shù)下算法的DOA估計精度對比如圖3所示.

圖3 不同快拍數(shù)下算法的DOA估計精度對比

由圖3可知,采樣快拍數(shù)對算法的DOA估計精度有顯著影響,隨采樣快拍數(shù)量的增加,4種算法的均方根誤差都呈現(xiàn)出減小的趨勢.但從仿真圖中可以清晰看出,本文算法的整體均方根誤差曲線在文獻[4]和文獻[6]之下,其均方根誤差曲線更加接近經(jīng)典Root-MUSIC算法及無偏估計量下界CRLB,擁有更高的DOA估計精度及更強的魯棒性.

實驗4:獨立入射信源數(shù)量變?yōu)?個,分別為-10°,10°和20°,其余實驗條件與實驗3一致.采樣快拍數(shù)同樣以50為間隔從100變化至600進行1 000次蒙特卡洛實驗,計算均方根誤差.不同快拍數(shù)下算法的DOA估計精度對比如圖4所示.

圖4 不同快拍數(shù)下算法的DOA估計精度對比

由圖4可知,算法的均方根誤差均隨著快拍數(shù)的增加而下降.但在入射信源個數(shù)增加至3個時,由于陣列規(guī)模的限制,文獻[6]由于算法自身需要拆分陣元進行互相關矩陣重構,其可利用陣元數(shù)量嚴重下降,DOA估計精度損失嚴重;文獻[4]則由于算法本身利用噪聲子空間信息較少,DOA估計精度同樣不足;本文算法及經(jīng)典Root-MUSIC算法的均方根誤差相較于2個入射信源時同樣有所上升,但總體更加趨近于無偏估計量下界CRLB,算法的DOA估計精度更高.

3.2 時間復雜度分析

本文算法的目的旨在解決目前多數(shù)低復雜度Root-MUSIC算法的DOA估計精度不足問題,那么所提算法自身的時間復雜度也不能有明顯提升,對3.1節(jié)中的4種算法的時間復雜度進行分析(其中T為采樣快拍數(shù)).算法各主要運算部分時間復雜度對比如表1所示.

表1 算法各主要運算部分時間復雜度對比

由表1可知,文獻[4]算法的求根運算時間復雜度最低;本文算法在求根多項式運算部分的時間復雜度上高于文獻[4],但低于文獻[6]及經(jīng)典Root-MUSIC算法;本文算法的協(xié)方差矩陣運算時間復雜度最低;在陣元規(guī)模較小時,本文算法及文獻[4]算法的子空間提取時間復雜度略高于文獻[6]算法,但隨著陣元規(guī)模的增大,本文及文獻[4]算法的子空間提取時間復雜度將低于文獻[6]算法.整體看,本文算法的綜合時間復雜度略高于文獻[4]算法,但低于文獻[6]及經(jīng)典Root-MUSIC算法.下面通過實驗5來驗證這一結果.

實驗5:兩個獨立入射信源角度分別為-10°和10°,信噪比固定為5 dB,采樣快拍數(shù)為200,陣元間距為λ/2.陣元數(shù)量以2為間隔從6變化至40,每種算法運行1 000次計算耗時.算法運行1 000次的耗時對比如圖5所示.

圖5 算法運行1 000次的時間耗時對比

由圖5可知,隨著陣元數(shù)量的增加,算法的運行耗時均有所增加.其中文獻[4]的算法耗時始終最低,這也說明了多項式求根運算是整個算法運行過程中耗時占比最高的部分;本文算法耗時高于文獻[4]算法,但始終低于文獻[6]算法及經(jīng)典Root-MUSIC算法,從表1的時間復雜度分析也可以驗證這一點.多項式求根運算的時間復雜度決定了算法整體的時間復雜度,而本文算法在不提高其他運算部分時間復雜度的同時,降低了多項式求根運算的時間復雜度,使得算法整體運行耗時降低.

4 結 論

針對多數(shù)低復雜度Root-MUSIC算法的DOA估計精度損失問題,本文提出了一種高精度低復雜度的改進Root-MUSIC算法.算法通過近似數(shù)據(jù)觀測矩陣首行重構具有 Hermitian 性質的Toeplitz矩陣,并將其作為新的自相關矩陣做特征值分解運算得到噪聲子空間;對噪聲子空間作翻轉拆分處理,重構新的求根多項式;最終通過求根運算得到DOA估計值.仿真實驗結果表明,相較于部分前人算法,本文算法擁有更高的DOA估計精度,其均方根誤差更接近于經(jīng)典Root-MUSIC算法及無偏估計量下界,同時擁有更高的魯棒性及穩(wěn)定性，其時間復雜度則基本不高于前人算法.