黃翔東,念天磊,馬 欣
(天津大學(xué)電氣自動(dòng)化與信息工程學(xué)院,天津 300072)
陣列波達(dá)方向(direction-of-arrival,DOA)[1]估計(jì)是信號(hào)處理領(lǐng)域的重要研究方向,其應(yīng)用十分廣泛。傳統(tǒng)方法一般采用陣元間距符合Nyquist采樣定理的陣列進(jìn)行接收如均勻線性陣列(uniform linear array,ULA),若要增加分辨率,就需通過(guò)增加陣元數(shù)量來(lái)擴(kuò)大陣列孔徑,但這會(huì)導(dǎo)致硬件成本的提高。另外,Nyquist采樣定理要求陣列間距不大于半波長(zhǎng),故隨著工作頻率的升高,信號(hào)波長(zhǎng)變短,陣元間距也隨之變短,增加陣元數(shù)量還會(huì)引起陣元間的互耦加大而降低接收性能[2]。因此,在不滿足Nyquist采樣定理情況下,僅依靠稀疏陣列[2-7]估計(jì)出多目標(biāo)的DOA成為急需突破的課題。很顯然,該課題繞不開兩個(gè)基本問題:一是選擇稀疏陣列結(jié)構(gòu)的布置,二是DOA估計(jì)算法的設(shè)計(jì)。
在選擇稀疏陣列結(jié)構(gòu)的布置方面,現(xiàn)有的稀疏陣列,例如嵌套陣列[3],最小冗余陣列(minimum redundancy array,MRA)[4-5],最小孔洞陣列[6]等稀疏非均勻陣列均可用于DOA估計(jì),并且在與均勻線陣陣元數(shù)目相同的情況下可以獲得更大的陣列孔徑,從而提高DOA分辨率。但這些陣列各有缺陷:最小孔洞陣列和最小冗余陣列的陣列結(jié)構(gòu)無(wú)法用具體的閉式表達(dá)式描述,這會(huì)給優(yōu)化設(shè)計(jì)帶來(lái)很大困難,其工程應(yīng)用受到很大限制;而用嵌套陣列做多源目標(biāo)估計(jì)時(shí),雖然可用N個(gè)陣元實(shí)現(xiàn)O(N2)的自由度[2](即提高了識(shí)別分辨率),但其缺陷在于:嵌套結(jié)構(gòu)使得稀疏陣列的觀測(cè)協(xié)方差矩陣不易轉(zhuǎn)化成Nyquist協(xié)方差矩陣(以方便后續(xù)子空間分解算法做源估計(jì)),造成目標(biāo)估計(jì)的困難。
在DOA算法設(shè)計(jì)方面,主要分為壓縮感知算法[5,8-15]和子空間分解算法兩大類。對(duì)于第1類,需要根據(jù)陣列結(jié)構(gòu)構(gòu)造壓縮感知系統(tǒng)模型,當(dāng)為密集陣列時(shí),直接以陣列流型矩陣作為觀測(cè)矩陣,以陣元數(shù)據(jù)做觀測(cè)向量,而構(gòu)造出壓縮感知模型[8-9];當(dāng)為稀疏陣列時(shí),則需對(duì)陣列流型矩陣的各列做克羅內(nèi)克積構(gòu)造觀測(cè)矩陣,對(duì)由陣元數(shù)據(jù)算出的協(xié)方差矩陣做向量化構(gòu)成觀測(cè)向量,而構(gòu)造出壓縮感知模型[5]。但無(wú)論陣列是密集還是稀疏情況,在構(gòu)造觀測(cè)矩陣時(shí),均需要構(gòu)建候選波達(dá)方向足夠豐富的過(guò)完備冗余字典(以保證真實(shí)DOA在該字典中呈現(xiàn)稀疏分布),進(jìn)而以該字典元素為變量構(gòu)造目標(biāo)函數(shù),通過(guò)匹配追蹤(matching pursuit,MP)[10]、基追蹤(basis pursuit,BP)[11]、正交MP(orthogonal MP,OMP)[12]、壓縮采樣MP(compressive sampling MP,CoSaMP)[13]等追蹤方法或1-奇異值分解(singular value decomposition,SVD)[14]、最小絕對(duì)值收斂和選擇算子 (least absolute shrinkage and selection operator,LASSO)[15]等手段對(duì)目標(biāo)函數(shù)尋最優(yōu),而獲得DOA估計(jì)。很明顯,冗余字典的構(gòu)造耗費(fèi)很大的存儲(chǔ)量,而目標(biāo)函數(shù)的求解涉及復(fù)雜的遞歸迭代過(guò)程,計(jì)算復(fù)雜度高,不利于實(shí)時(shí)應(yīng)用。對(duì)于第2類子空間分解法,通常采用主流的、技術(shù)成熟的多信號(hào)分類(multiple signal classification,MUSIC)[16]方式,故其計(jì)算復(fù)雜度比壓縮感知要低,具有更高的工程應(yīng)用價(jià)值,該方法關(guān)鍵在于:如何把稀疏陣列的協(xié)方差矩陣轉(zhuǎn)換為Nyquist虛擬陣列的協(xié)方差矩陣。
近年來(lái),基于互素欠采樣[7,17-21]的譜分析成為信號(hào)處理學(xué)科的研究熱點(diǎn)。相比于傳統(tǒng)Nyquist樣本分布,互素欠采樣更為稀疏,這意味著在只需布置同樣數(shù)目的陣元即可獲得更大的陣列孔徑,故有利于提升空間譜分辨率;且互素欠采樣有中國(guó)余數(shù)定理[17]支撐,故易于推導(dǎo)出閉合形式的理論分析結(jié)果,這是嵌套陣列、最小冗余陣列、最小孔洞陣列遠(yuǎn)不能及的。文獻(xiàn)[22]利用矩陣填充思想對(duì)稀疏協(xié)方差矩陣進(jìn)行擴(kuò)維構(gòu)造一個(gè)更大維度的Toeplitz矩陣(空位補(bǔ)零),然后利用不定點(diǎn)延續(xù)算法對(duì)零元素進(jìn)行有效替換,轉(zhuǎn)化為完整的協(xié)方差矩陣。此算法可獲得很大的陣列孔徑,從而獲得很高的自由度,但該算法除了需做MUSIC分解外,在矩陣填充過(guò)程中,還消耗了一系列較復(fù)雜的最優(yōu)化操作(以保證零元素替換后矩陣的秩盡可能低)。
本文選定互素稀疏陣列實(shí)現(xiàn)DOA估計(jì),筆者在前期工作中,已經(jīng)把互素采樣應(yīng)用于時(shí)間序列譜分析[17-19]和DOA估計(jì)[20-21],但文獻(xiàn)[20-21]是通過(guò)簡(jiǎn)單的頻譜校正途徑來(lái)做DOA估計(jì),該方法不能估計(jì)多源目標(biāo)。為實(shí)現(xiàn)多目標(biāo)DOA估計(jì),本文提出差集表遍歷搜索方法,借助該方法,可以很容易地將互素稀疏陣列的協(xié)方差矩陣轉(zhuǎn)換為Nyquist虛擬陣列的協(xié)方差矩陣,其擴(kuò)張的陣列孔徑尺寸雖然不如文獻(xiàn)[22],但該方法可直接實(shí)現(xiàn)從觀測(cè)協(xié)方差矩陣到虛擬協(xié)方差矩陣的解析映射,且無(wú)需矩陣填充的優(yōu)化過(guò)程,除MUSIC分解外,其他附加的計(jì)算量可忽略不計(jì)。故出于計(jì)算復(fù)雜度與分辨率的權(quán)衡,本文估計(jì)器具有較高實(shí)用價(jià)值。
令Nyquist采樣間距為d=λ/2,圖1給出的互素陣列由兩個(gè)均勻線性子陣組成[23],其陣元間距分別為Md和Nd(M、N為互素正整數(shù)),陣元個(gè)數(shù)分別為N和2M,由于兩個(gè)子陣的第1個(gè)陣元的空間坐標(biāo)重合,故整個(gè)互素陣列的陣元個(gè)數(shù)為L(zhǎng)=N+2M-1,其歸一化(d為歸一化因子)后的坐標(biāo)集合S可閉式表達(dá)為
S= {mM,0≤m≤N-1}∪
{nN,0≤n≤2M-1}
(1)
圖1 互素陣列模型Fig.1 Coprime array model
式(1)和圖1表明,互素陣列的陣元間距允許遠(yuǎn)大于Nyquist間距d=λ/2,故具有很高的稀疏度(表現(xiàn)為陣元間存在較多的“空洞”)。
令式(1)的升序排列的坐標(biāo)集合為
P=sort(S)={p0,pl,…,pN+2M-2}
(2)
以第1個(gè)陣元坐標(biāo)p0=0為參考,假設(shè)遠(yuǎn)場(chǎng)有D個(gè)不相關(guān)的窄帶信號(hào)s1,s2,…,sD,分別以到達(dá)角θ1,θ2,…,θD入射到接收陣列中,則互素陣列接收到的信號(hào)模型為
(3)
其中
a(θi)=[1,e-jπp1sin(θi),…,e-jπpN+2M-2sin(θi)]T
(4)
對(duì)于陣元位置直接按Nyquist間距設(shè)定情況,只需把L×1的接收陣元向量x與其轉(zhuǎn)置乘積后,做統(tǒng)計(jì)平均(假定K為快拍數(shù)量),即可獲得子空間分解所需的統(tǒng)計(jì)協(xié)方差矩陣,即
(5)
當(dāng)信號(hào)源不相干時(shí),有rank(Rxx)=L,可直接做MUSIC分解識(shí)別L-1個(gè)源。令陣元快拍的空間協(xié)方差函數(shù)為φ(l),當(dāng)快拍充足時(shí),Rxx近似為Toeplitz陣。理想情況下,Rxx內(nèi)部元素表示為
Rxx(u,v)=φ(u-v),0≤u,v≤L-1
(6)
進(jìn)而,對(duì)于某個(gè)φ(l)值(l=-L+1,…,0,…,L-1)在矩陣Rxx對(duì)應(yīng)的坐標(biāo)(u,v)滿足
l=u-v
(7)
從式(5)~式(7)的理想定義可看出,經(jīng)典MUSIC分解要求協(xié)方差矩陣的元素行、列坐標(biāo)的差值正好為陣元?dú)w一化坐標(biāo)的差值,這決定了φ(l)值分布在與矩陣主對(duì)角線平行的線上。
然而,互素矩陣與該情況不同,由于陣元稀疏分布,在矩陣Rxx內(nèi),與φ(l)所對(duì)應(yīng)的Rxx內(nèi)的坐標(biāo)(u,v)不再滿足式(6)、式(7)的Toeplitz關(guān)系。因而,問題的關(guān)鍵在于:如何將不符合Toeplitz關(guān)系的Rxx轉(zhuǎn)換為符合Toeplitz關(guān)系的Nyquist虛擬陣列(即將圖1互素陣列的“空洞”填滿、全部陣元間距為d=λ/2的假想陣列)的協(xié)方差矩陣。
具體而言,從圖1可知:互素協(xié)方差矩陣Rxx內(nèi)的任一元素Rxx(u,v),一定與互素陣列的兩個(gè)陣元坐標(biāo)pu、pv相對(duì)應(yīng),其對(duì)應(yīng)關(guān)系描述為
Rxx(u,v)=φ(pu-pv),0≤u,v≤L-1
(8)
另外,由于pu或pv只能取自子陣1或子陣2,即
(9)
聯(lián)立式(8)、式(9),可推知Rxx(u,v)與子陣1坐標(biāo)m,m′、子陣2的坐標(biāo)n,n′的對(duì)應(yīng)關(guān)系,不外乎4種情況,即
(10)
式中,m,m′∈[0,N-1];n,n′∈[0,2M-1]。
結(jié)合圖1的互素陣列結(jié)構(gòu),可發(fā)現(xiàn):式(10)的前兩種情況為“自差”協(xié)方差函數(shù)值(由同一子陣上的兩陣元的快拍做“自相關(guān)”而得),后兩種情況為“互差”協(xié)方差函數(shù)(子陣1的某陣元的快拍與子陣2的某陣元的快拍做“互相關(guān)”而得)。根據(jù)中國(guó)余數(shù)定理[17,24-25],任取l∈{-MN,…,0,…,MN},在以上自差和互差的4種情況中,至少有一種情況與φ(l)相對(duì)應(yīng)。即
(11)
聯(lián)立式(2)、式(8)~式(11)可知,在獲得升序坐標(biāo)集合P后,為把Nyquist協(xié)方差函數(shù)的空間時(shí)延l與觀測(cè)陣元的協(xié)方差矩陣Rxx的行標(biāo)u、列標(biāo)v關(guān)聯(lián)起來(lái),不妨構(gòu)造一張尺寸為L(zhǎng)×L的“差集表”T,該差集表的內(nèi)部元素為
T(u,v)=pu-pv,u,v∈[0,L-1]
(12)
由以上分析可看出,差集表T的生成僅涉及式排序和式減法操作,故過(guò)程簡(jiǎn)單。
注意在差集表T中,對(duì)于某Nyquist空間時(shí)延l∈{-MN,…,0,…,MN},可能存在有多個(gè)表目T(u,v)值等于l的情況,故需把對(duì)應(yīng)的行標(biāo)rl、列標(biāo)值cl全部搜索出來(lái),即
(13)
(14)
借助式(14)計(jì)算出2MN+1個(gè)協(xié)方差函數(shù)估計(jì)值,即可構(gòu)造如式(15)的Nyquist虛擬陣列的協(xié)方差矩陣
(15)
相比較而言,文獻(xiàn)[22]則需要構(gòu)造(2MN-N+1)×(2MN-N+1)的Toeplitz矩陣,然而由于在延遲范圍l∈{-MN-M+1,…,0,…,MN+M-1}以外的波程差不連續(xù),故在擴(kuò)展的協(xié)方差矩陣中存在零元素(其數(shù)目隨M、N值呈指數(shù)增長(zhǎng)),故需以矩陣秩最小化為目標(biāo)函數(shù)對(duì)這些令元素做填充,該最優(yōu)化雖然可等效為核范數(shù)最小化問題[22],但在矩陣維度(或秩)很大時(shí)求解復(fù)雜度依然較大。
需強(qiáng)調(diào),式(15)構(gòu)造的畢竟是虛擬的協(xié)方差矩陣,與密集陣列構(gòu)造的真實(shí)協(xié)方差矩陣有所不同。具體體現(xiàn)在:當(dāng)各陣元采集的快拍數(shù)受限時(shí),真實(shí)的協(xié)方差矩陣的對(duì)角線元素仍存在差別,而式(15)的對(duì)角線元素是相等的。這種虛擬構(gòu)造方式是以略微犧牲估計(jì)器的抗噪魯棒性為代價(jià)的,后面的仿真實(shí)驗(yàn)也將證明這一點(diǎn)。
根據(jù)以上描述,可將本文提出的互素稀疏陣列DOA估計(jì)算法總結(jié)為表1的流程。
表1 基于差集表遍歷搜索的DOA估計(jì)流程Table 1 DOA estimation flow based on difference-set table traversal searching
例1令目標(biāo)源個(gè)數(shù)D=4,其信號(hào)功率均設(shè)置為1,其入射角θ1,θ2,…,θ4分別為-10°,20°,25°,60°,陣元數(shù)目L=6,每個(gè)陣元上采樣的快拍數(shù)K=2 000,信噪比(signal-to-noise ratio,SNR)設(shè)置為SNR=0 dB(為高斯白噪聲),試分別用互素稀疏陣列和傳統(tǒng)Nyquist間距的均勻線性陣列做信號(hào)接收,并分別用本文提出的差集表遍歷搜索法和傳統(tǒng)直接MUSIC分解法進(jìn)行DOA估計(jì),對(duì)照其空間譜。
(1) 互素稀疏陣列下基于差集表遍歷搜索的DOA估計(jì)結(jié)果
設(shè)定互素整數(shù)對(duì)M=2,N=3(恰好L=N+2M-1=6),根據(jù)式(1)構(gòu)造互素稀疏陣列,根據(jù)式(3)、式(4)的系統(tǒng)模型而得到各陣元快拍數(shù)據(jù),對(duì)這些數(shù)據(jù)處理過(guò)程如下。
根據(jù)表1步驟1,可得到升序排列的坐標(biāo)集P={0,2,3,4,6,9}和6×6的差集表T(見表2)。
根據(jù)表1步驟2對(duì)表2列出的T(u,v)做遍歷搜索,可得到如表3所示的坐標(biāo)集{(rl,cl)}。
表2 基于差集表遍歷搜索的DOA估計(jì)流程Table 2 DOA estimation process based on difference-set table traversal searching
表3 對(duì)差集表做遍歷搜索后的坐標(biāo)集Table 3 Coordinate sets by difference-set table traversal searching
(2) 進(jìn)而根據(jù)表1步驟3、步驟4,可得到如圖2所示的空間譜
圖2 本文方法得到的空間譜圖(D=4)Fig.2 Spatial spectrum by proposed method (D=4)
(3) 傳統(tǒng)Nyquist間距下的均勻線性陣列接收結(jié)果
圖3給出了傳統(tǒng)Nyquist間距的ULA情況,直接用經(jīng)典MUSIC分解法得到的空間譜。
圖3 傳統(tǒng)Nyquist-ULA的空間譜圖(D=4)Fig.3 Spectrum of traditional Nyquist-ULA (D=4)
對(duì)比圖2和圖3,可得出如下結(jié)論。
(1) 雖然耗費(fèi)同樣數(shù)目的接收陣元(L=6)和快拍(K=2 000),用本文提出的基于差集表遍歷搜索算法的互素陣列DOA估計(jì)精度高于基于傳統(tǒng)Nyquist-ULA的估計(jì)精度,表現(xiàn)在前者的譜峰(見圖2)可精確定位在期望的-10°,20°,25°,60°位置上(虛線所示),而后者的譜峰(見圖3)與這些位置有明顯偏離。
(2) 圖2的DOA譜峰形狀明顯比圖3的譜峰形狀尖銳得多,表現(xiàn)在本文方法仍能分辨出兩個(gè)密集入射角θ2=20°,θ3=25°,而圖3的基于Nyquist-ULA的方法則只能將這兩個(gè)入射角識(shí)別為1個(gè)。這是因?yàn)榛ニ仃嚵械臍w一化陣元排列位置為P={0,2,3,4,6,9},而Nyquist-ULA的陣元排列位置為{0,1,2,3,4,5}。因而,在耗費(fèi)同樣陣元數(shù)目情況下,互素陣列由于陣元稀疏放置,具有比傳統(tǒng)Nyquist陣列更大的孔徑,結(jié)合本文提出的差集表遍歷搜索方法,可獲得更高的空間譜區(qū)分度。
(3) 進(jìn)一步分析,兩者性能差異的根本原因在于:如前述,依據(jù)本文提出的差集表遍歷搜索算法可構(gòu)造出(MN+1)×(MN+1)=7×7的協(xié)方差矩陣,從而容許的信號(hào)子空間最高維度為MN=6,而傳統(tǒng)Nyquist-ULA只能構(gòu)造出L×L=6×6的協(xié)方差矩陣,故容許的信號(hào)子空間最高維度為L(zhǎng)-1=5。故對(duì)D=4個(gè)目標(biāo)源情況,本文信號(hào)子空間容許度高于后者,因而獲得了更高的接收性能。
進(jìn)而可推出:M、N值越大,本文方法的性能改善更明顯,下面給出例子說(shuō)明。
例2仍分別用互素稀疏陣列和傳統(tǒng)Nyquist均勻線性陣列做信號(hào)接收,令目標(biāo)源個(gè)數(shù)D=15,其入射角θ1,θ2,…,θ15設(shè)置為[-50°,62°]范圍內(nèi)的15個(gè)均勻間隔值,陣元數(shù)目L=10(對(duì)于互素陣列則要求M=3,N=5),其他參數(shù)設(shè)置與例1相同。圖4(a)、圖4(b)分別給出了本文方法和傳統(tǒng)直接MUSIC分解法的DOA識(shí)別譜圖。
圖4 不同陣列設(shè)置下兩種方法的DOA檢測(cè)譜圖Fig.4 DOA spectrums of different array configurations based on prosed method and traditional method
從圖4的DOA檢測(cè)譜圖可發(fā)現(xiàn):圖4(a)中,基于互素陣列的本文方法仍能區(qū)分[-50°,62°]范圍內(nèi)的15個(gè)譜峰,而對(duì)于Nyquist-ULA的檢測(cè)方法,由于目標(biāo)源數(shù)D=15超出了其允許的信號(hào)子空間最高維度(即L-1=9),導(dǎo)致圖4(b)中的DOA檢測(cè)基本失效。這驗(yàn)證了差集表遍歷搜索法可提升信號(hào)子空間的容許維度的結(jié)論。
例3仍分別用互素稀疏陣列和傳統(tǒng)Nyquist-ULA做信號(hào)接收,令目標(biāo)源數(shù)D=1,其入射角θ1設(shè)置為45°,其他參數(shù)設(shè)置與例2相同。在-20~25 dB的SNR區(qū)間內(nèi)(間隔為1 dB)對(duì)兩種方法做蒙特卡羅測(cè)試(對(duì)于每種SNR情況,做2 000次測(cè)試),并統(tǒng)計(jì)均方根誤差(root-of-mean-square-error,RMSE),其RMSE曲線如圖5所示。
圖5 RMSE隨信噪比變化情況(L=10,M=3,N=5) Fig.5 Change of RMSE with SNR (L=10,M=3,N=5)
從圖5測(cè)試結(jié)果可得出結(jié)論:
(1)整體上,本文方法的RMSE低于傳統(tǒng)Nyquist-ULA的RMSE,因而具有更高的估計(jì)精度,這仍是差集表遍歷搜索方法可提升信號(hào)子空間的容許維度、增強(qiáng)譜峰分辨能力的反映。
(2)在低SNR區(qū)域,兩個(gè)估計(jì)器都會(huì)出現(xiàn)閾值效應(yīng)(即SNR低于該閾值時(shí)估計(jì)器失效),從圖5可看出,本文方法的SNR閾值比傳統(tǒng)Nyquist-ULA的SNR閾值要高1 dB,這說(shuō)明前者相比于后者,抗噪魯棒性略低一些,這可看作是構(gòu)造虛擬協(xié)防差矩陣所付出的代價(jià)。
統(tǒng)計(jì)RMSE隨快拍數(shù)量K的變化,其中SNR設(shè)為10 dB,K取100~2 000間隔為100,結(jié)果如圖6所示。
圖6 波達(dá)角RMSE隨快拍數(shù)量變化情況(L=10,M=3,N=5)Fig.6 Change of DOA RMSE with snapshot numbers (L=10,M=3,N=5)
可以得出結(jié)論:隨著快拍數(shù)量的增加,兩種估計(jì)器的準(zhǔn)確度均會(huì)提高。在快拍數(shù)量一定的條件下基于差集表遍歷搜索的估計(jì)方法誤差更小。
在Matlab R2016b(Windows10)平臺(tái)上對(duì)兩種估計(jì)器的執(zhí)行時(shí)間進(jìn)行比較,處理器為i5。兩種算法分別執(zhí)行300次,MUSIC掃描間隔設(shè)為0.05,陣元數(shù)量為10(N=5,M=3),2 000個(gè)快拍。提出的估計(jì)器用時(shí)30.0 s,Nyquist-ULA估計(jì)器用時(shí)28.5 s。僅僅多消耗了5.2%的時(shí)間就可以獲得O(MN)的自由度,獲得更高的估計(jì)精度。
本文以互素陣列為基礎(chǔ),提出基于差集表遍歷搜索的DOA估計(jì)方法。該差集表僅需依據(jù)互素稀疏陣列的陣元坐標(biāo)即可構(gòu)造出來(lái),以該差集表為指導(dǎo),無(wú)需任何附加的優(yōu)化迭代措施,即可實(shí)現(xiàn)觀測(cè)陣元的協(xié)方差矩陣到Nyquist虛擬陣列的快速轉(zhuǎn)換。理論分析和仿真實(shí)驗(yàn)均證明:借助差集表遍歷搜索方法,可以保證互素陣列只耗費(fèi)L=N+2M-1個(gè)陣元,即可獲得MN個(gè)源目標(biāo)的區(qū)分度,從而相比于Nyquist間距的傳統(tǒng)陣列,大大提升了空間多目標(biāo)分辨能力和DOA估計(jì)精度。
近年來(lái)隨著S波段、C波段、X波段、K波段等高頻段、短波長(zhǎng)雷達(dá)應(yīng)用逐漸展開,陣元間耦合的副效應(yīng)逐漸趨于嚴(yán)重。由于本文提出的基于互素稀疏陣列的差集表遍歷搜索法可放寬傳統(tǒng)Nyquist陣列的密集布置陣元的限制,因而具有較廣泛的工程應(yīng)用前景。