基于差集表遍歷搜索的互素陣列DOA估計器

黃翔東,念天磊,馬 欣

(天津大學電氣自動化與信息工程學院,天津 300072)

0 引 言

陣列波達方向(direction-of-arrival,DOA)[1]估計是信號處理領域的重要研究方向,其應用十分廣泛。傳統方法一般采用陣元間距符合Nyquist采樣定理的陣列進行接收如均勻線性陣列(uniform linear array,ULA),若要增加分辨率,就需通過增加陣元數量來擴大陣列孔徑,但這會導致硬件成本的提高。另外,Nyquist采樣定理要求陣列間距不大于半波長,故隨著工作頻率的升高,信號波長變短,陣元間距也隨之變短,增加陣元數量還會引起陣元間的互耦加大而降低接收性能[2]。因此,在不滿足Nyquist采樣定理情況下,僅依靠稀疏陣列[2-7]估計出多目標的DOA成為急需突破的課題。很顯然,該課題繞不開兩個基本問題:一是選擇稀疏陣列結構的布置,二是DOA估計算法的設計。

在選擇稀疏陣列結構的布置方面,現有的稀疏陣列,例如嵌套陣列[3],最小冗余陣列(minimum redundancy array,MRA)[4-5],最小孔洞陣列[6]等稀疏非均勻陣列均可用于DOA估計,并且在與均勻線陣陣元數目相同的情況下可以獲得更大的陣列孔徑,從而提高DOA分辨率。但這些陣列各有缺陷:最小孔洞陣列和最小冗余陣列的陣列結構無法用具體的閉式表達式描述,這會給優化設計帶來很大困難,其工程應用受到很大限制;而用嵌套陣列做多源目標估計時,雖然可用N個陣元實現O(N2)的自由度[2](即提高了識別分辨率),但其缺陷在于:嵌套結構使得稀疏陣列的觀測協方差矩陣不易轉化成Nyquist協方差矩陣(以方便后續子空間分解算法做源估計),造成目標估計的困難。

在DOA算法設計方面,主要分為壓縮感知算法[5,8-15]和子空間分解算法兩大類。對于第1類,需要根據陣列結構構造壓縮感知系統模型,當為密集陣列時,直接以陣列流型矩陣作為觀測矩陣,以陣元數據做觀測向量,而構造出壓縮感知模型[8-9];當為稀疏陣列時,則需對陣列流型矩陣的各列做克羅內克積構造觀測矩陣,對由陣元數據算出的協方差矩陣做向量化構成觀測向量,而構造出壓縮感知模型[5]。但無論陣列是密集還是稀疏情況,在構造觀測矩陣時,均需要構建候選波達方向足夠豐富的過完備冗余字典(以保證真實DOA在該字典中呈現稀疏分布),進而以該字典元素為變量構造目標函數,通過匹配追蹤(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]、最小絕對值收斂和選擇算子 (least absolute shrinkage and selection operator,LASSO)[15]等手段對目標函數尋最優,而獲得DOA估計。很明顯,冗余字典的構造耗費很大的存儲量,而目標函數的求解涉及復雜的遞歸迭代過程,計算復雜度高,不利于實時應用。對于第2類子空間分解法,通常采用主流的、技術成熟的多信號分類(multiple signal classification,MUSIC)[16]方式,故其計算復雜度比壓縮感知要低,具有更高的工程應用價值,該方法關鍵在于:如何把稀疏陣列的協方差矩陣轉換為Nyquist虛擬陣列的協方差矩陣。

近年來,基于互素欠采樣[7,17-21]的譜分析成為信號處理學科的研究熱點。相比于傳統Nyquist樣本分布,互素欠采樣更為稀疏,這意味著在只需布置同樣數目的陣元即可獲得更大的陣列孔徑,故有利于提升空間譜分辨率;且互素欠采樣有中國余數定理[17]支撐,故易于推導出閉合形式的理論分析結果,這是嵌套陣列、最小冗余陣列、最小孔洞陣列遠不能及的。文獻[22]利用矩陣填充思想對稀疏協方差矩陣進行擴維構造一個更大維度的Toeplitz矩陣(空位補零),然后利用不定點延續算法對零元素進行有效替換,轉化為完整的協方差矩陣。此算法可獲得很大的陣列孔徑,從而獲得很高的自由度,但該算法除了需做MUSIC分解外,在矩陣填充過程中,還消耗了一系列較復雜的最優化操作(以保證零元素替換后矩陣的秩盡可能低)。

本文選定互素稀疏陣列實現DOA估計,筆者在前期工作中,已經把互素采樣應用于時間序列譜分析[17-19]和DOA估計[20-21],但文獻[20-21]是通過簡單的頻譜校正途徑來做DOA估計,該方法不能估計多源目標。為實現多目標DOA估計,本文提出差集表遍歷搜索方法,借助該方法,可以很容易地將互素稀疏陣列的協方差矩陣轉換為Nyquist虛擬陣列的協方差矩陣,其擴張的陣列孔徑尺寸雖然不如文獻[22],但該方法可直接實現從觀測協方差矩陣到虛擬協方差矩陣的解析映射,且無需矩陣填充的優化過程,除MUSIC分解外,其他附加的計算量可忽略不計。故出于計算復雜度與分辨率的權衡,本文估計器具有較高實用價值。

1 系統模型

令Nyquist采樣間距為d=λ/2,圖1給出的互素陣列由兩個均勻線性子陣組成[23],其陣元間距分別為Md和Nd(M、N為互素正整數),陣元個數分別為N和2M,由于兩個子陣的第1個陣元的空間坐標重合,故整個互素陣列的陣元個數為L=N+2M-1,其歸一化(d為歸一化因子)后的坐標集合S可閉式表達為

S= {mM,0≤m≤N-1}∪

{nN,0≤n≤2M-1}

(1)

圖1 互素陣列模型Fig.1 Coprime array model

式(1)和圖1表明,互素陣列的陣元間距允許遠大于Nyquist間距d=λ/2,故具有很高的稀疏度(表現為陣元間存在較多的“空洞”)。

令式(1)的升序排列的坐標集合為

P=sort(S)={p0,pl,…,pN+2M-2}

(2)

以第1個陣元坐標p0=0為參考,假設遠場有D個不相關的窄帶信號s1,s2,…,sD,分別以到達角θ1,θ2,…,θD入射到接收陣列中,則互素陣列接收到的信號模型為

(3)

其中

a(θi)=[1,e-jπp1sin(θi),…,e-jπpN+2M-2sin(θi)]T

(4)

2 基于差集表遍歷搜索的互素陣列DOA估計器

2.1 互素矩陣的差集表構造

對于陣元位置直接按Nyquist間距設定情況,只需把L×1的接收陣元向量x與其轉置乘積后,做統計平均(假定K為快拍數量),即可獲得子空間分解所需的統計協方差矩陣,即

(5)

當信號源不相干時,有rank(Rxx)=L,可直接做MUSIC分解識別L-1個源。令陣元快拍的空間協方差函數為φ(l),當快拍充足時,Rxx近似為Toeplitz陣。理想情況下,Rxx內部元素表示為

Rxx(u,v)=φ(u-v),0≤u,v≤L-1

(6)

進而,對于某個φ(l)值(l=-L+1,…,0,…,L-1)在矩陣Rxx對應的坐標(u,v)滿足

l=u-v

(7)

從式(5)～式(7)的理想定義可看出,經典MUSIC分解要求協方差矩陣的元素行、列坐標的差值正好為陣元歸一化坐標的差值,這決定了φ(l)值分布在與矩陣主對角線平行的線上。

然而,互素矩陣與該情況不同,由于陣元稀疏分布,在矩陣Rxx內,與φ(l)所對應的Rxx內的坐標(u,v)不再滿足式(6)、式(7)的Toeplitz關系。因而,問題的關鍵在于:如何將不符合Toeplitz關系的Rxx轉換為符合Toeplitz關系的Nyquist虛擬陣列(即將圖1互素陣列的“空洞”填滿、全部陣元間距為d=λ/2的假想陣列)的協方差矩陣。

具體而言,從圖1可知:互素協方差矩陣Rxx內的任一元素Rxx(u,v),一定與互素陣列的兩個陣元坐標pu、pv相對應,其對應關系描述為

Rxx(u,v)=φ(pu-pv),0≤u,v≤L-1

(8)

另外,由于pu或pv只能取自子陣1或子陣2,即

(9)

聯立式(8)、式(9),可推知Rxx(u,v)與子陣1坐標m,m′、子陣2的坐標n,n′的對應關系,不外乎4種情況,即

(10)

式中,m,m′∈[0,N-1];n,n′∈[0,2M-1]。

結合圖1的互素陣列結構,可發現:式(10)的前兩種情況為“自差”協方差函數值(由同一子陣上的兩陣元的快拍做“自相關”而得),后兩種情況為“互差”協方差函數(子陣1的某陣元的快拍與子陣2的某陣元的快拍做“互相關”而得)。根據中國余數定理[17,24-25],任取l∈{-MN,…,0,…,MN},在以上自差和互差的4種情況中,至少有一種情況與φ(l)相對應。即

(11)

聯立式(2)、式(8)～式(11)可知,在獲得升序坐標集合P后,為把Nyquist協方差函數的空間時延l與觀測陣元的協方差矩陣Rxx的行標u、列標v關聯起來,不妨構造一張尺寸為L×L的“差集表”T,該差集表的內部元素為

T(u,v)=pu-pv,u,v∈[0,L-1]

(12)

由以上分析可看出,差集表T的生成僅涉及式排序和式減法操作,故過程簡單。

2.2 基于差集表遍歷搜索的Nyquist協方差函數估計

注意在差集表T中,對于某Nyquist空間時延l∈{-MN,…,0,…,MN},可能存在有多個表目T(u,v)值等于l的情況,故需把對應的行標rl、列標值cl全部搜索出來,即

(13)

(14)

借助式(14)計算出2MN+1個協方差函數估計值,即可構造如式(15)的Nyquist虛擬陣列的協方差矩陣

(15)

相比較而言,文獻[22]則需要構造(2MN-N+1)×(2MN-N+1)的Toeplitz矩陣,然而由于在延遲范圍l∈{-MN-M+1,…,0,…,MN+M-1}以外的波程差不連續,故在擴展的協方差矩陣中存在零元素(其數目隨M、N值呈指數增長),故需以矩陣秩最小化為目標函數對這些令元素做填充,該最優化雖然可等效為核范數最小化問題[22],但在矩陣維度(或秩)很大時求解復雜度依然較大。

需強調,式(15)構造的畢竟是虛擬的協方差矩陣,與密集陣列構造的真實協方差矩陣有所不同。具體體現在:當各陣元采集的快拍數受限時,真實的協方差矩陣的對角線元素仍存在差別,而式(15)的對角線元素是相等的。這種虛擬構造方式是以略微犧牲估計器的抗噪魯棒性為代價的,后面的仿真實驗也將證明這一點。

2.3 DOA估計算法流程總結

根據以上描述,可將本文提出的互素稀疏陣列DOA估計算法總結為表1的流程。

表1 基于差集表遍歷搜索的DOA估計流程Table 1 DOA estimation flow based on difference-set table traversal searching

3 仿真實驗

3.1 本文算法的處理流程解析與性能對比

例1令目標源個數D=4,其信號功率均設置為1,其入射角θ1,θ2,…,θ4分別為-10°,20°,25°,60°,陣元數目L=6,每個陣元上采樣的快拍數K=2 000,信噪比(signal-to-noise ratio,SNR)設置為SNR=0 dB(為高斯白噪聲),試分別用互素稀疏陣列和傳統Nyquist間距的均勻線性陣列做信號接收,并分別用本文提出的差集表遍歷搜索法和傳統直接MUSIC分解法進行DOA估計,對照其空間譜。

(1) 互素稀疏陣列下基于差集表遍歷搜索的DOA估計結果

設定互素整數對M=2,N=3(恰好L=N+2M-1=6),根據式(1)構造互素稀疏陣列,根據式(3)、式(4)的系統模型而得到各陣元快拍數據,對這些數據處理過程如下。

根據表1步驟1,可得到升序排列的坐標集P={0,2,3,4,6,9}和6×6的差集表T(見表2)。

根據表1步驟2對表2列出的T(u,v)做遍歷搜索,可得到如表3所示的坐標集{(rl,cl)}。

表2 基于差集表遍歷搜索的DOA估計流程Table 2 DOA estimation process based on difference-set table traversal searching

表3 對差集表做遍歷搜索后的坐標集Table 3 Coordinate sets by difference-set table traversal searching

(2) 進而根據表1步驟3、步驟4,可得到如圖2所示的空間譜

圖2 本文方法得到的空間譜圖(D=4)Fig.2 Spatial spectrum by proposed method (D=4)

(3) 傳統Nyquist間距下的均勻線性陣列接收結果

圖3給出了傳統Nyquist間距的ULA情況,直接用經典MUSIC分解法得到的空間譜。

圖3 傳統Nyquist-ULA的空間譜圖(D=4)Fig.3 Spectrum of traditional Nyquist-ULA (D=4)

對比圖2和圖3,可得出如下結論。

(1) 雖然耗費同樣數目的接收陣元(L=6)和快拍(K=2 000),用本文提出的基于差集表遍歷搜索算法的互素陣列DOA估計精度高于基于傳統Nyquist-ULA的估計精度,表現在前者的譜峰(見圖2)可精確定位在期望的-10°,20°,25°,60°位置上(虛線所示),而后者的譜峰(見圖3)與這些位置有明顯偏離。

(2) 圖2的DOA譜峰形狀明顯比圖3的譜峰形狀尖銳得多,表現在本文方法仍能分辨出兩個密集入射角θ2=20°,θ3=25°,而圖3的基于Nyquist-ULA的方法則只能將這兩個入射角識別為1個。這是因為互素陣列的歸一化陣元排列位置為P={0,2,3,4,6,9},而Nyquist-ULA的陣元排列位置為{0,1,2,3,4,5}。因而,在耗費同樣陣元數目情況下,互素陣列由于陣元稀疏放置,具有比傳統Nyquist陣列更大的孔徑,結合本文提出的差集表遍歷搜索方法,可獲得更高的空間譜區分度。

(3) 進一步分析,兩者性能差異的根本原因在于:如前述,依據本文提出的差集表遍歷搜索算法可構造出(MN+1)×(MN+1)=7×7的協方差矩陣,從而容許的信號子空間最高維度為MN=6,而傳統Nyquist-ULA只能構造出L×L=6×6的協方差矩陣,故容許的信號子空間最高維度為L-1=5。故對D=4個目標源情況,本文信號子空間容許度高于后者,因而獲得了更高的接收性能。

進而可推出:M、N值越大,本文方法的性能改善更明顯,下面給出例子說明。

3.2 提高陣元數目情況的多目標檢測性能對比

例2仍分別用互素稀疏陣列和傳統Nyquist均勻線性陣列做信號接收,令目標源個數D=15,其入射角θ1,θ2,…,θ15設置為[-50°,62°]范圍內的15個均勻間隔值,陣元數目L=10(對于互素陣列則要求M=3,N=5),其他參數設置與例1相同。圖4(a)、圖4(b)分別給出了本文方法和傳統直接MUSIC分解法的DOA識別譜圖。

圖4 不同陣列設置下兩種方法的DOA檢測譜圖Fig.4 DOA spectrums of different array configurations based on prosed method and traditional method

從圖4的DOA檢測譜圖可發現:圖4(a)中,基于互素陣列的本文方法仍能區分[-50°,62°]范圍內的15個譜峰,而對于Nyquist-ULA的檢測方法,由于目標源數D=15超出了其允許的信號子空間最高維度(即L-1=9),導致圖4(b)中的DOA檢測基本失效。這驗證了差集表遍歷搜索法可提升信號子空間的容許維度的結論。

3.3 檢測精度和魯棒性對比

例3仍分別用互素稀疏陣列和傳統Nyquist-ULA做信號接收,令目標源數D=1,其入射角θ1設置為45°,其他參數設置與例2相同。在-20～25 dB的SNR區間內(間隔為1 dB)對兩種方法做蒙特卡羅測試(對于每種SNR情況,做2 000次測試),并統計均方根誤差(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測試結果可得出結論:

(1)整體上,本文方法的RMSE低于傳統Nyquist-ULA的RMSE,因而具有更高的估計精度,這仍是差集表遍歷搜索方法可提升信號子空間的容許維度、增強譜峰分辨能力的反映。

(2)在低SNR區域,兩個估計器都會出現閾值效應(即SNR低于該閾值時估計器失效),從圖5可看出,本文方法的SNR閾值比傳統Nyquist-ULA的SNR閾值要高1 dB,這說明前者相比于后者,抗噪魯棒性略低一些,這可看作是構造虛擬協防差矩陣所付出的代價。

統計RMSE隨快拍數量K的變化,其中SNR設為10 dB,K取100～2 000間隔為100,結果如圖6所示。

圖6 波達角RMSE隨快拍數量變化情況(L=10,M=3,N=5)Fig.6 Change of DOA RMSE with snapshot numbers (L=10,M=3,N=5)

可以得出結論:隨著快拍數量的增加,兩種估計器的準確度均會提高。在快拍數量一定的條件下基于差集表遍歷搜索的估計方法誤差更小。

3.4 計算時間對比

在Matlab R2016b(Windows10)平臺上對兩種估計器的執行時間進行比較,處理器為i5。兩種算法分別執行300次,MUSIC掃描間隔設為0.05,陣元數量為10(N=5,M=3),2 000個快拍。提出的估計器用時30.0 s,Nyquist-ULA估計器用時28.5 s。僅僅多消耗了5.2%的時間就可以獲得O(MN)的自由度,獲得更高的估計精度。

4 結 論

本文以互素陣列為基礎,提出基于差集表遍歷搜索的DOA估計方法。該差集表僅需依據互素稀疏陣列的陣元坐標即可構造出來,以該差集表為指導,無需任何附加的優化迭代措施,即可實現觀測陣元的協方差矩陣到Nyquist虛擬陣列的快速轉換。理論分析和仿真實驗均證明:借助差集表遍歷搜索方法,可以保證互素陣列只耗費L=N+2M-1個陣元,即可獲得MN個源目標的區分度,從而相比于Nyquist間距的傳統陣列,大大提升了空間多目標分辨能力和DOA估計精度。

近年來隨著S波段、C波段、X波段、K波段等高頻段、短波長雷達應用逐漸展開,陣元間耦合的副效應逐漸趨于嚴重。由于本文提出的基于互素稀疏陣列的差集表遍歷搜索法可放寬傳統Nyquist陣列的密集布置陣元的限制,因而具有較廣泛的工程應用前景。