佘黎煌, 劉平凡, 張 石, 許方晗
(東北大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院, 遼寧 沈陽(yáng) 110169)
陣列信號(hào)處理技術(shù)包括DOA估計(jì)技術(shù)及波束形成技術(shù),廣泛應(yīng)用于空間信號(hào)探測(cè)、雷達(dá)及水聲信號(hào)探測(cè)等諸多領(lǐng)域[1].其中DOA估計(jì)技術(shù)作為獲得信號(hào)波達(dá)方向角度的重要技術(shù),如何使DOA估計(jì)算法獲得更高的估計(jì)精度和更低的時(shí)間復(fù)雜度一直是科研人員的主要研究方向.
在眾多DOA估計(jì)算法中,Root-MUSIC算法以其較高的估計(jì)精度、較低的復(fù)雜度受到眾多科研人員的重視[2],大量基于Root-MUSIC算法的改進(jìn)算法被提出.傳統(tǒng)Root-MUSIC算法存在的最大問(wèn)題在于隨著陣列規(guī)模的增大,算法的運(yùn)算負(fù)擔(dān)成倍增加.在較低陣列規(guī)模時(shí),多項(xiàng)式求根使算法運(yùn)行效率較高,而當(dāng)陣列規(guī)模增大后,其求根多項(xiàng)式階數(shù)快速升高,算法運(yùn)行耗時(shí)甚至超過(guò)譜峰搜索類算法[3].目前眾多基于Root-MUSIC算法的改進(jìn)算法都著重針對(duì)求根多項(xiàng)式進(jìn)行降維操作,通過(guò)較低的數(shù)據(jù)規(guī)模實(shí)現(xiàn)DOA估計(jì).Ren等[4]提出了一種快速Root-MUSIC算法,通過(guò)將噪聲子空間進(jìn)行拆分,構(gòu)造新的數(shù)據(jù)向量,使得求根多項(xiàng)式階數(shù)降低至僅為信號(hào)源個(gè)數(shù),大大降低了算法運(yùn)行復(fù)雜度.由于噪聲子空間信息利用率較低,算法精度損失嚴(yán)重,王新賀等[5]在Ren等[4]拆分噪聲子空間矩陣的基礎(chǔ)上重構(gòu)數(shù)據(jù)向量,充分利用噪聲子空間,算法精度有了大幅度提升,但運(yùn)算復(fù)雜度也隨之增大,同時(shí)算法在陣列規(guī)模較大時(shí)完全失效,不具備魯棒性.Yan等[6]將均勻線陣拆分為兩個(gè)相等子陣列,利用兩個(gè)子陣列采集的數(shù)據(jù)構(gòu)造3個(gè)無(wú)噪聲互相關(guān)矩陣,用處理后的新數(shù)據(jù)矩陣進(jìn)行DOA估計(jì),并將此DOA估計(jì)數(shù)據(jù)作為參考值,重復(fù)以上方法,進(jìn)行二次估計(jì)得到最終DOA估計(jì)值[7].算法通過(guò)拆分子陣構(gòu)造新互相關(guān)矩陣的方法避免了特征分解運(yùn)算,降低了求根運(yùn)算的復(fù)雜度,但在低陣列規(guī)模時(shí)算法精度損失嚴(yán)重.探測(cè)陣列必須為偶數(shù),可估計(jì)信源個(gè)數(shù)較少,魯棒性較差.本文在Toeplitz矩陣重構(gòu)及求根多項(xiàng)式降階的基礎(chǔ)上,提出了一種改進(jìn)的Root-MUSIC算法,算法依據(jù)數(shù)據(jù)觀測(cè)矩陣首行重構(gòu)Toeplitz矩陣.特征值分解后得到噪聲子空間,并將噪聲子空間翻轉(zhuǎn)拆分,重構(gòu)新的求根多項(xiàng)式向量,進(jìn)而通過(guò)求根方法得到DOA估計(jì)值.算法通過(guò)重構(gòu)Toeplitz矩陣使得新自相關(guān)矩陣具備Hermitian性質(zhì),相比于傳統(tǒng)自相關(guān)矩陣特征值分解算法,本文算法擁有更高的DOA估計(jì)精度.噪聲子空間的翻轉(zhuǎn)拆分重構(gòu)了新的低階求根多項(xiàng)式,有效降低了算法的時(shí)間復(fù)雜度.
設(shè)存在M個(gè)滿足遠(yuǎn)場(chǎng)條件的各向同性陣元所組成的均勻線陣模型,同時(shí)有q個(gè)遠(yuǎn)場(chǎng)窄帶目標(biāo)信號(hào)源,其入射角為θk(k=1,2,…,q),并且滿足M>q.選擇首個(gè)陣元作為參考陣,陣元間距設(shè)為d,則此時(shí)陣列響應(yīng)向量為
(1)
第m個(gè)陣元的輸出為
(2)
式中:sk(t)為第k個(gè)由天線陣列接收到的信源信號(hào);λ為信號(hào)的波長(zhǎng);q為信號(hào)源數(shù)目;nm(t)為陣元接收的高斯白噪聲,方差為σ2.此時(shí)的陣列接收數(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)
通過(guò)式(3)可以得到理想觀測(cè)數(shù)據(jù)的協(xié)方差矩陣:
R=E[X(t)XH(t)]=ARsAH+σ2IM.
(8)
式中:Rs=E[S(t)SH(t)]為信號(hào)的協(xié)方差矩陣;IM為M階的單位矩陣.
由于硬件設(shè)備采樣快拍數(shù)的限制,由數(shù)據(jù)觀測(cè)矩陣獲得的自相關(guān)矩陣并不能滿足理想的Hermitian性質(zhì),僅能保證對(duì)角占優(yōu),因此存在估計(jì)誤差.Toeplitz矩陣重構(gòu)算法通過(guò)將數(shù)據(jù)觀測(cè)矩陣的第一行即首個(gè)陣元獲得的數(shù)據(jù)作為參考向量,構(gòu)造新的自相關(guān)矩陣使得新矩陣具備Hermitian性質(zhì).相較于原始自相關(guān)矩陣更加逼近理想值,DOA估計(jì)精度更高,同時(shí)不會(huì)產(chǎn)生陣列孔徑損失,估計(jì)信源個(gè)數(shù)不會(huì)受到影響[8].
由式(2)可以得到數(shù)據(jù)觀測(cè)矩陣第一行即首個(gè)陣元的輸出為
(9)
在實(shí)際情況下,由于采樣快拍數(shù)有限,因此天線陣列所接收到的信源信號(hào)及高斯白噪聲均表現(xiàn)為由若干采樣值所組成的向量形式,向量維度即采樣快拍數(shù),這里將采樣快拍數(shù)設(shè)為T,則有
x1=[x1(1),x1(2),…,x1(T)] .
(10)
同理,第k個(gè)陣元的接收數(shù)據(jù)向量表示為
xk=[xk(1),xk(2),…,xk(T)] .
(11)
將式(10)和式(11)作為參考向量,定義
(12)
式中,A(1)和A(k)分別代表矩陣A的第一行和第k行.隨著k的變化,可以得到包含全部信源信息的相關(guān)矢量[rtp(0),rtp(1),…,rtp(M-1)],由M個(gè)相關(guān)函數(shù)構(gòu)成的矩陣:
Rtp=
(13)
譜峰搜索算法通過(guò)全范圍角度的遍歷來(lái)搜索DOA估計(jì)值,在陣元數(shù)較少時(shí)會(huì)導(dǎo)致算法運(yùn)行耗時(shí)增大.而求根MUSIC算法(Root-MUSIC)將這一過(guò)程省去,替代為簡(jiǎn)化的多項(xiàng)式求根過(guò)程.
通過(guò)式(1)可以得到陣列響應(yīng)向量,這里定義:
(14)
式(1)的陣列響應(yīng)向量可以改寫為
(15)
式(4)陣列流型矩陣可以表示為
A=[z1,z2,…,zq] .
(16)
根據(jù)式(8)可以得到自相關(guān)矩陣R,對(duì)R做特征分解,較大的q個(gè)特征值對(duì)應(yīng)特征向量構(gòu)成信號(hào)子空間Us,較小的M-q個(gè)特征值構(gòu)成噪聲子空間Un.依據(jù)MUSIC算法相關(guān)理論,信號(hào)子空間與噪聲子空間相互正交,則噪聲子空間也必與信號(hào)子空間的基即陣列流型矩陣A正交,滿足
AHUn=Oq×(M-q).
(17)
式中,O為q×(M-q)階零矩陣.根據(jù)式(17)定義噪聲子空間Un:
Un=[un1,un2,…,un(M-q)] .
(18)
此時(shí)式(17)可以表示為
(19)
通過(guò)觀察式(19),可將陣列響應(yīng)方向向量定義為一般形式:
z=[1,z-1,…,z-(M-1)]T.
(20)
根據(jù)式(19)及式(20)構(gòu)造多項(xiàng)式向量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個(gè)多項(xiàng)式pi(z)的階數(shù)均為M-1,且有q個(gè)公共根z1,z2,…,zq.
根據(jù)文獻(xiàn)[7]構(gòu)造多項(xiàng)式:
froot(z)=pH(z)p(z) .
(23)
求得式(23)的根即可獲得有關(guān)信號(hào)源到達(dá)角的信息.因多項(xiàng)式存在z*項(xiàng),這會(huì)使得求根過(guò)程變得復(fù)雜,因此對(duì)式(23)修正:
froot(z)=zM-1pT(z-1)p(z) .
(24)
在式(24)中,多項(xiàng)式階數(shù)為2M-2,即有M-1對(duì)根,且每對(duì)根為相互共軛關(guān)系.在這些根中有q個(gè)根z1,z2,…,zq剛好分布在單位圓上.
在實(shí)際情況下因采樣快拍數(shù)的限制無(wú)法獲得理想的自相關(guān)矩陣,因此實(shí)際得到的近似協(xié)方差矩陣存在一定程度的誤差,這時(shí)只需要對(duì)式(24)求近似的q個(gè)根即可(最靠近單位圓).在均勻直線陣下的DOA估計(jì)角度計(jì)算式為
(25)
在Root-MUSIC算法中,時(shí)間復(fù)雜度最高的運(yùn)算為多項(xiàng)式求根,求根運(yùn)算的時(shí)間復(fù)雜度為O(2M-1)3,可知,隨陣元規(guī)模的擴(kuò)大,其時(shí)間復(fù)雜度迅速增加,實(shí)驗(yàn)也證實(shí)了這一點(diǎn).當(dāng)陣元規(guī)模達(dá)到32以上時(shí),其耗時(shí)已經(jīng)超過(guò)了譜峰搜索算法,因此本節(jié)將提出一種降低求根多項(xiàng)式階數(shù)的低復(fù)雜度求根運(yùn)算方法.
根據(jù)2.2節(jié)可知,設(shè)存在多項(xiàng)式p1(z),p2(z),…,pM-q(z)且同時(shí)擁有q個(gè)公共根z1,z2,…,zq,則它們之間存在線性組合Q(z):
(26)
接下來(lái)證明式(26)的合理性,若r為多項(xiàng)式p1(z),p2(z),…,pM-q(z)的公共根,則r為多項(xiàng)式p1(z),p2(z),…,pM-q(z)任意線性組合的根.
設(shè)存在c1,c2,…,cM-q使得Q(z)=∑cipi(x),則將公共根r代入Q(z)可得
Q(r)=∑cipi(r)=∑ci×0=0 .
(27)
通過(guò)式(27)證明了式(26)的合理性,構(gòu)造多項(xiàng)式p1(z),p2(z),…,pM-q(z)的線性組合Q(z),對(duì)應(yīng)系數(shù)設(shè)為α1,α2,…,αM-q,令αM-q=1,可以得到如下表達(dá)式:
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對(duì)應(yīng)階系數(shù)為零,得到如下方程組:
(30)
式(30)可以寫成如下的矩陣形式:
(31)
將式(31)簡(jiǎn)化整理為
x=-H-1h.
(32)
通過(guò)式(32)求得x,得到了對(duì)應(yīng)式(26)中的系數(shù)αi,1≤i≤M-q-1.要得到x,就必須得到矩陣H及向量h.
在2.2節(jié)的討論中,對(duì)自相關(guān)矩陣R進(jìn)行特征值分解得到噪聲子空間Un,其列向量uni=[β0,i,β1,i,…,βM-1,i]T,則將噪聲子空間列向量翻轉(zhuǎn)可得:
(33)
(34)
(35)
依據(jù)式(32)及式(35)計(jì)算得到線性組合Q(z)的系數(shù)向量x,并將式(21)代入式(25)得到一個(gè)新的求根多項(xiàng)式:
(36)
式中:αM-q=1.通過(guò)式(36)所構(gòu)造的新求根多項(xiàng)式來(lái)替代式(24)完成求根運(yùn)算,尋找最靠近單位圓的q個(gè)根,最終通過(guò)式(25)求得DOA估計(jì)角度.
本文算法的整體實(shí)現(xiàn)過(guò)程如下:
1) 根據(jù)式(12)及式(13)構(gòu)造近似數(shù)據(jù)觀測(cè)矩陣下的Hermitian Toeplitz矩陣Rtp;
2) 對(duì)矩陣Rtp做特征值分解,得到噪聲子空間矩陣Un;
3) 依據(jù)式(33)~式(35),對(duì)Un翻轉(zhuǎn)拆分獲得矩陣H及向量h;
4) 通過(guò)計(jì)算式(32)得到系數(shù)向量x,并根據(jù)式(36)構(gòu)造新的求根多項(xiàng)式Q(z);
5) 完成求根運(yùn)算,并尋找最靠近單位圓的q個(gè)根,通過(guò)式(25)求得DOA估計(jì)角度.
針對(duì)改進(jìn)算法的DOA估計(jì)性能及算法運(yùn)行耗時(shí)進(jìn)行仿真測(cè)試.為驗(yàn)證本文算法的各項(xiàng)性能,仿真實(shí)驗(yàn)分別對(duì)比了文獻(xiàn)[4]的快速Root-MUSIC算法、文獻(xiàn)[6]的無(wú)需特征分解計(jì)算的兩步Root-MUSIC算法和經(jīng)典Root-MUSIC算法;同時(shí),對(duì)幾種算法的時(shí)間復(fù)雜度進(jìn)行了分析.
均方根誤差計(jì)算式為
(37)
實(shí)驗(yàn)1:兩個(gè)獨(dú)立入射信源角度分別為-10°和10°;陣列為8陣元均勻直線陣列,陣列間距為λ/2;采樣快拍數(shù)為200,信噪比以2 dB間隔從-10 dB變化至10 dB,進(jìn)行1 000次蒙特卡洛實(shí)驗(yàn),計(jì)算均方根誤差.不同信噪比下算法的DOA估計(jì)精度對(duì)比如圖1所示.
圖1 不同信噪比下算法的DOA估計(jì)精度對(duì)比
由圖1可知,隨著信噪比的增加,4種算法的均方根誤差均有所下降,而本文算法的估計(jì)誤差最接近經(jīng)典Root-MUSIC算法,低于文獻(xiàn)[4]及文獻(xiàn)[6]所提算法;同時(shí)可以發(fā)現(xiàn),本文算法及經(jīng)典算法更加接近算法的無(wú)偏估計(jì)量下界(克拉美羅下界)CRLB,擁有更高的估計(jì)精度.
實(shí)驗(yàn)2: 獨(dú)立入射信源數(shù)量為3個(gè),分別為-10°,10°和20°,其余實(shí)驗(yàn)條件與實(shí)驗(yàn)1一致,同樣使信噪比以2 dB間隔從-10 dB變化至10 dB進(jìn)行1 000次蒙特卡洛實(shí)驗(yàn),計(jì)算均方根誤差.
不同信噪比下算法的DOA估計(jì)精度對(duì)比如圖2所示.
由圖2可知,當(dāng)入射信源增加為3個(gè)時(shí),由于陣列規(guī)模的限制,文獻(xiàn)[4]和文獻(xiàn)[6]的算法已經(jīng)基本失效,而本文算法及經(jīng)典Root-MUSIC算法在低信噪比環(huán)境下效果較差,隨信噪比的增加,均方根誤差都有較大幅度的下降.受限于陣列規(guī)模,相較于2個(gè)入射信源的情況本文算法及經(jīng)典Root-MUSIC算法的DOA估計(jì)性能均有所下降.
圖2 不同信噪比下算法的DOA估計(jì)精度對(duì)比
實(shí)驗(yàn)3: 2個(gè)獨(dú)立入射信源角度為-10°和10°,陣列仍采用8陣元均勻直線陣列,陣列間距為λ/2,信噪比固定為5 dB.采樣快拍數(shù)以50為間隔從100變化至600進(jìn)行1 000次蒙特卡洛實(shí)驗(yàn),計(jì)算均方根誤差.不同快拍數(shù)下算法的DOA估計(jì)精度對(duì)比如圖3所示.
圖3 不同快拍數(shù)下算法的DOA估計(jì)精度對(duì)比
由圖3可知,采樣快拍數(shù)對(duì)算法的DOA估計(jì)精度有顯著影響,隨采樣快拍數(shù)量的增加,4種算法的均方根誤差都呈現(xiàn)出減小的趨勢(shì).但從仿真圖中可以清晰看出,本文算法的整體均方根誤差曲線在文獻(xiàn)[4]和文獻(xiàn)[6]之下,其均方根誤差曲線更加接近經(jīng)典Root-MUSIC算法及無(wú)偏估計(jì)量下界CRLB,擁有更高的DOA估計(jì)精度及更強(qiáng)的魯棒性.
實(shí)驗(yàn)4:獨(dú)立入射信源數(shù)量變?yōu)?個(gè),分別為-10°,10°和20°,其余實(shí)驗(yàn)條件與實(shí)驗(yàn)3一致.采樣快拍數(shù)同樣以50為間隔從100變化至600進(jìn)行1 000次蒙特卡洛實(shí)驗(yàn),計(jì)算均方根誤差.不同快拍數(shù)下算法的DOA估計(jì)精度對(duì)比如圖4所示.
圖4 不同快拍數(shù)下算法的DOA估計(jì)精度對(duì)比
由圖4可知,算法的均方根誤差均隨著快拍數(shù)的增加而下降.但在入射信源個(gè)數(shù)增加至3個(gè)時(shí),由于陣列規(guī)模的限制,文獻(xiàn)[6]由于算法自身需要拆分陣元進(jìn)行互相關(guān)矩陣重構(gòu),其可利用陣元數(shù)量嚴(yán)重下降,DOA估計(jì)精度損失嚴(yán)重;文獻(xiàn)[4]則由于算法本身利用噪聲子空間信息較少,DOA估計(jì)精度同樣不足;本文算法及經(jīng)典Root-MUSIC算法的均方根誤差相較于2個(gè)入射信源時(shí)同樣有所上升,但總體更加趨近于無(wú)偏估計(jì)量下界CRLB,算法的DOA估計(jì)精度更高.
本文算法的目的旨在解決目前多數(shù)低復(fù)雜度Root-MUSIC算法的DOA估計(jì)精度不足問(wèn)題,那么所提算法自身的時(shí)間復(fù)雜度也不能有明顯提升,對(duì)3.1節(jié)中的4種算法的時(shí)間復(fù)雜度進(jìn)行分析(其中T為采樣快拍數(shù)).算法各主要運(yùn)算部分時(shí)間復(fù)雜度對(duì)比如表1所示.
表1 算法各主要運(yùn)算部分時(shí)間復(fù)雜度對(duì)比
由表1可知,文獻(xiàn)[4]算法的求根運(yùn)算時(shí)間復(fù)雜度最低;本文算法在求根多項(xiàng)式運(yùn)算部分的時(shí)間復(fù)雜度上高于文獻(xiàn)[4],但低于文獻(xiàn)[6]及經(jīng)典Root-MUSIC算法;本文算法的協(xié)方差矩陣運(yùn)算時(shí)間復(fù)雜度最低;在陣元規(guī)模較小時(shí),本文算法及文獻(xiàn)[4]算法的子空間提取時(shí)間復(fù)雜度略高于文獻(xiàn)[6]算法,但隨著陣元規(guī)模的增大,本文及文獻(xiàn)[4]算法的子空間提取時(shí)間復(fù)雜度將低于文獻(xiàn)[6]算法.整體看,本文算法的綜合時(shí)間復(fù)雜度略高于文獻(xiàn)[4]算法,但低于文獻(xiàn)[6]及經(jīng)典Root-MUSIC算法.下面通過(guò)實(shí)驗(yàn)5來(lái)驗(yàn)證這一結(jié)果.
實(shí)驗(yàn)5:兩個(gè)獨(dú)立入射信源角度分別為-10°和10°,信噪比固定為5 dB,采樣快拍數(shù)為200,陣元間距為λ/2.陣元數(shù)量以2為間隔從6變化至40,每種算法運(yùn)行1 000次計(jì)算耗時(shí).算法運(yùn)行1 000次的耗時(shí)對(duì)比如圖5所示.
圖5 算法運(yùn)行1 000次的時(shí)間耗時(shí)對(duì)比
由圖5可知,隨著陣元數(shù)量的增加,算法的運(yùn)行耗時(shí)均有所增加.其中文獻(xiàn)[4]的算法耗時(shí)始終最低,這也說(shuō)明了多項(xiàng)式求根運(yùn)算是整個(gè)算法運(yùn)行過(guò)程中耗時(shí)占比最高的部分;本文算法耗時(shí)高于文獻(xiàn)[4]算法,但始終低于文獻(xiàn)[6]算法及經(jīng)典Root-MUSIC算法,從表1的時(shí)間復(fù)雜度分析也可以驗(yàn)證這一點(diǎn).多項(xiàng)式求根運(yùn)算的時(shí)間復(fù)雜度決定了算法整體的時(shí)間復(fù)雜度,而本文算法在不提高其他運(yùn)算部分時(shí)間復(fù)雜度的同時(shí),降低了多項(xiàng)式求根運(yùn)算的時(shí)間復(fù)雜度,使得算法整體運(yùn)行耗時(shí)降低.
針對(duì)多數(shù)低復(fù)雜度Root-MUSIC算法的DOA估計(jì)精度損失問(wèn)題,本文提出了一種高精度低復(fù)雜度的改進(jìn)Root-MUSIC算法.算法通過(guò)近似數(shù)據(jù)觀測(cè)矩陣首行重構(gòu)具有 Hermitian 性質(zhì)的Toeplitz矩陣,并將其作為新的自相關(guān)矩陣做特征值分解運(yùn)算得到噪聲子空間;對(duì)噪聲子空間作翻轉(zhuǎn)拆分處理,重構(gòu)新的求根多項(xiàng)式;最終通過(guò)求根運(yùn)算得到DOA估計(jì)值.仿真實(shí)驗(yàn)結(jié)果表明,相較于部分前人算法,本文算法擁有更高的DOA估計(jì)精度,其均方根誤差更接近于經(jīng)典Root-MUSIC算法及無(wú)偏估計(jì)量下界,同時(shí)擁有更高的魯棒性及穩(wěn)定性,其時(shí)間復(fù)雜度則基本不高于前人算法.