王嘉偉,楊赟秀,陳文東,舒 勤

(1.四川大學(xué) 電氣工程學(xué)院,成都 610065;2.西南技術(shù)物理研究所,成都 610045;3.北京遙感設(shè)備研究所,北京 100005)

0 引 言

在陣列信號(hào)處理領(lǐng)域,波達(dá)角(Direction of Arrival,DOA)估計(jì)一直是一個(gè)重要的技術(shù)問(wèn)題,在無(wú)線通信、雷達(dá)和聲吶等許多領(lǐng)域中有著廣泛的研究和應(yīng)用[1]。通常情況下,陣元數(shù)目是影響估計(jì)效果的一個(gè)重要因素,但是由于成本和運(yùn)算復(fù)雜度的制約,它無(wú)法被無(wú)限增加。而由稀疏陣列產(chǎn)生的虛擬陣列則在很大程度上突破了這一限制,它發(fā)端于壓縮感知(Compressed Sensing,CS)理論[2-3],通過(guò)計(jì)算和處理接收信號(hào)的協(xié)方差矩陣或偽協(xié)方差矩陣,從數(shù)學(xué)層面實(shí)現(xiàn)了對(duì)陣列的拓展,從而能以更少的初始陣元數(shù)獲得更大的陣列孔徑與陣元自由度(Degrees-of-freedom,DOF),有著極為廣泛的應(yīng)用前景。為此,研究人員已經(jīng)提出了許多基于稀疏陣列的高分辨率估計(jì)算法,包括SS-MUSIC(Spatial Smoothing MUSIC)[4]、SS-ESPRIT(Spatial Smoothing Estimation of Signal Parameters via Rotational Invariance techniques)[5]等,但它們都需要連續(xù)虛擬陣元去進(jìn)行下一步處理(往往是空間平滑操作),這時(shí)隨機(jī)分布在其上的孔洞就成為制約其連續(xù)自由度增長(zhǎng)的關(guān)鍵。

當(dāng)前,對(duì)于虛擬孔洞的研究往往集中于一維估計(jì)領(lǐng)域,主要分為前處理與后處理兩個(gè)思路。前處理是指更加合理地設(shè)計(jì)初始陣列的排布,例如最小冗余陣列[6]、嵌套陣列[7]、改進(jìn)的互質(zhì)陣列[8]、CACIS[9]、CADiS[9]等,從而盡可能地減少虛擬孔洞的數(shù)目及其帶來(lái)的影響。而后處理則是指在收到信號(hào)后,對(duì)于虛擬孔洞位置的空缺信息進(jìn)行近似化填充,例如核范數(shù)最小化[10]、協(xié)方差矩陣的半正定結(jié)構(gòu)[11]等,它們通過(guò)近似估計(jì)孔洞處應(yīng)有的接收信號(hào),將本該舍棄的虛擬陣元信息更加充分地利用起來(lái)。

在二維DOA估計(jì)領(lǐng)域,對(duì)于虛擬孔洞的研究則較少,往往是直接將孔洞之后的陣元信息舍去。這造成了一定的資源浪費(fèi),但若直接將一維孔洞填充思路拓展至二維,則會(huì)極大地增加運(yùn)算復(fù)雜度,估計(jì)效果也不很理想。針對(duì)于此,本文提出一種基于二維DOA估計(jì)的虛擬孔洞填充算法,它結(jié)合平行稀疏陣列,將繁瑣的二維估計(jì)問(wèn)題轉(zhuǎn)化為兩個(gè)簡(jiǎn)易的一維估計(jì)問(wèn)題,并通過(guò)提前計(jì)算其中一子陣列虛擬孔洞位置去設(shè)計(jì)另一子陣列的排布,從而實(shí)現(xiàn)對(duì)孔洞的完全填充。本文算法是一種典型的前處理思路,它提升了陣列資源利用率與陣列自由度,能夠以更少的初始陣元獲得更加精確的估計(jì)效果。

本文所用符號(hào)說(shuō)明:(·)T,(·)*,(·)H,(·)-1和(·)?分別表示矩陣的轉(zhuǎn)置、共軛、共軛轉(zhuǎn)置、求逆和求偽逆;diag(·)表示將向量或?qū)顷噷?duì)角化;vec(·)表示將矩陣向量化;E(·)表示求矩陣的統(tǒng)計(jì)期望;?和⊙分別表示Kronecker積和Khatri-Rao積;det(·)表示矩陣行列式;root(·)和arg(·)分別表示求根和求相位角;IK是維度為K×K的單位矩陣;d為基本陣元間距,它被設(shè)置為接收信號(hào)的半波長(zhǎng)λ/2。

1 信號(hào)模型

本文的陣列模型采用平行稀疏陣列結(jié)構(gòu),由兩個(gè)相互平行的稀疏線性陣列組成,如圖1所示,包含子陣1和子陣2,它們分別有L1和L2個(gè)陣元。這里,分別用P和Q代表兩子陣的陣元位置,表示為

圖1 平行稀疏陣列幾何模型[12-13]

(1)

(2)

式中:子陣1的方向矩陣A=[a1,a2,…,aK]的維度為L(zhǎng)1×K,子陣2的方向矩陣B=[b1,b2,…,bK]的維度為L(zhǎng)2×K,它們都只包含角度α信息;響應(yīng)矢量ak=[efkp1d,efkp2d,…,efkpL1d]T,bk=[efkq1d,efkq2d,…,efkqL2d]T對(duì)應(yīng)第k個(gè)入射信號(hào)(pld和qld分別為式(1)中表示的兩子陣的陣元位置);C=diag([ejπcos β1,ejπcos β2,…,ejπcos βK])為K×K的對(duì)角矩陣;n1(t)和n2(t)為兩不相關(guān)的零均值加性高斯白噪聲,且都與信號(hào)不相關(guān)。

2 基于平行稀疏陣列虛擬孔洞填充的二維DOA估計(jì)算法

2.1 稀疏陣列排布結(jié)構(gòu)

為了方便,下面依次介紹互質(zhì)陣列[15](Co-prime Array,CA)、改進(jìn)互質(zhì)陣列[8](Augmented Co-prime Array,ACA)、CACIS[9]和CADiS[9]等4種會(huì)在文中用到的稀疏陣列:

PCA={mNd,0≤m≤M-1}∪{nMd,0≤n≤N-1},

(3)

PACA={mNd,0≤m≤2M-1}∪{nMd,0≤n≤N-1},

(4)

(5)

(6)

2.2 虛擬陣列構(gòu)造

為了擴(kuò)充初始陣列為虛擬陣列,需要計(jì)算子陣1接收信號(hào)的協(xié)方差矩陣[16]:

(7)

通過(guò)對(duì)式(7)進(jìn)行特征值分解,可從中估計(jì)出噪聲功率(陣元數(shù)目小于信號(hào)數(shù)目時(shí)不再適用)并可依此去除噪聲影響[12],可得

(8)

接著,對(duì)上式做向量化處理,便可得到基于協(xié)方差矩陣的虛擬陣列

(9)

考慮到其中存在重復(fù)的元素,需要剔除重復(fù)元素,從而可得

(10)

式中:F=diag(Rs)。記差集U={uld,1≤l≤L3}={(pi-pj)d|pi,pj∈P},則有虛擬陣列方向矩陣U=[u1,u2,…,uK]和其響應(yīng)矢量uk=[efku1d,efku2d,…,efkuL3d]T。

為了進(jìn)一步處理,傳統(tǒng)算法往往需要提取式(10)中的連續(xù)陣元,而將虛擬孔洞后的陣元信息直接舍棄,這造成了一定程度的資源浪費(fèi)。

2.3 虛擬孔洞填充的二維DOA估計(jì)算法

首先,需要計(jì)算子陣2接收信號(hào)的協(xié)方差矩陣,并對(duì)其做類似于式(8)的去噪處理:

(11)

接著,對(duì)式(11)做類似式(9)的向量化處理可得

(12)

將其中重復(fù)元素剔除,可得

(13)

記差集G={gld,1≤l≤L4}={(qi-qj)d|qi,qj∈Q},則有虛擬陣列方向矩陣G=[g1,g2,…,gK]和其響應(yīng)矢量gk=[efkg1d,efkg2d,…,efkgL4d]T。

通過(guò)比較、觀察式(13)和式(10)可以發(fā)現(xiàn),兩者在形式上具有一致性,且都是只包含角度α這一個(gè)未知量,只是虛擬陣元位置由于兩子陣排布的不同而不同。因此,只要能夠合理設(shè)計(jì)兩子陣的排布,便能夠?qū)μ摂M孔洞進(jìn)行完全填充,從而能夠大幅提升陣列資源利用率。

如圖2 所示,為了能夠以更少的物理陣元獲得更好的虛擬孔洞填充效果,需要提前計(jì)算基于子陣1獲得的虛擬陣列中正半孔洞的位置H。因此,可將子陣2的陣元位置設(shè)定為Q={0,H},易知基于此陣列獲得的差集G可以對(duì)子陣1得到的虛擬孔洞完全填充,從而提升了虛擬陣列的連續(xù)自由度與估計(jì)精確度。

圖2 傳統(tǒng)算法的虛擬陣列、物理子陣2以及新算法的虛擬陣列(陣列參數(shù)設(shè)置為(4,5)1)

顯然,由式(10)可得其虛擬孔洞位置,并據(jù)此在式(13)結(jié)果中尋找該位置應(yīng)得的接收信號(hào),之后只需將該值填充到對(duì)應(yīng)位置便可以獲得一個(gè)擁有更大連續(xù)自由度的陣列接收信號(hào):

(14)

式中:填充后的虛擬陣列方向矩陣H=[h1,h2,…,hK]的維度為L(zhǎng)5×K,其響應(yīng)矢量為hk=[efkh1d,efkh2d,…,efkhL5d]T。

由于式(14)是一個(gè)單快拍信號(hào),為了獲得其協(xié)方差矩陣,需采用空間平滑算法,即將其分割成W=(L5+1)/2個(gè)子陣列,并求取子陣列協(xié)方差矩陣的算術(shù)平均值:

(15)

(16)

基于式(16),可得角度α的譜峰搜索表達(dá)式[4]:

(17)

通過(guò)式(17),將二維估計(jì)轉(zhuǎn)化為一維搜索,有效地降低了運(yùn)算復(fù)雜度。為了進(jìn)一步降低算法復(fù)雜度,本文采用求根MUSIC[17]的思路,也即用zg1取代H(αg1)中的ejπcos αg1,并可將上述問(wèn)題轉(zhuǎn)化為求取下式的K個(gè)根:

(18)

從而可得角度α的估計(jì)值

(19)

隨后,為了獲得角度β的估計(jì)值,并完成角度自動(dòng)配對(duì)過(guò)程,需要計(jì)算互協(xié)方差矩陣。由于兩子陣接收噪聲不相關(guān),故而相較于兩子陣的協(xié)方差矩陣,該式?jīng)]有噪聲項(xiàng)。

(20)

對(duì)上式向量化處理可得

(21)

式中:J=diag(CRs)。移除其中的重復(fù)元素可得

(22)

式中:方向矩陣D=[d1,d2,…,dK]的估計(jì)結(jié)果可由式(19)獲得,維度記作L6×K,其中包含著角度α的順序信息,故而可自動(dòng)實(shí)現(xiàn)角度匹配過(guò)程。

而J可通過(guò)求解如下最小化問(wèn)題得到[18]:

(23)

式中:‖·‖F(xiàn)表示Frobenius范數(shù)。其結(jié)果可近似于求解

(24)

從而,可得角度β的估計(jì)值

βk=arccos(arg(J)/π),k=1,2,…,K。

(25)

通過(guò)聯(lián)立式(19)和式(25),可得相互匹配的方位角和俯仰角,即

(26)

(27)

2.4 算法步驟

基于以上分析,下面給出本文算法的具體步驟:

Step1 選擇合適的稀疏陣列作為子陣1,通過(guò)計(jì)算其差集G得到正半虛擬孔洞位置H,并據(jù)此設(shè)置子陣2的陣元位置為Q={0,H}。

Step2 根據(jù)實(shí)際應(yīng)用,通過(guò)有限的快拍對(duì)兩子陣列的協(xié)方差矩陣和兩者的互協(xié)方差矩陣進(jìn)行估計(jì),即

(28)

(29)

(30)

式中:Y為接收信號(hào)的快拍數(shù)。

Step3 向量化處理子陣1和2的協(xié)方差矩陣,獲得虛擬陣列1和2,并將后者的結(jié)果填充于前者的虛擬孔洞處,從而獲得一個(gè)擁有更多連續(xù)自由度的虛擬陣列3。

3 實(shí)驗(yàn)分析

3.1 陣列自由度分析

采用上文提到的四種稀疏陣列,假定初始陣元數(shù)目相同,比較傳統(tǒng)算法與本文算法可獲得的連續(xù)自由度。經(jīng)過(guò)分析,可得嚴(yán)格的數(shù)目通式如表1所示。

表1 陣列的連續(xù)自由度[9,19]

由表1可知,虛擬孔洞填充后的陣列連續(xù)自由度有了較大提升,大多數(shù)情況都能有1倍甚至更多的提升。分析可知,新算法下獲得的虛擬陣列的連續(xù)自由度值已經(jīng)達(dá)到了單條物理陣列能達(dá)到的理論極限值,因此本算法能夠極大地提升陣元利用率,并由此獲得了更為精確的二維DOA估計(jì)結(jié)果。

3.2 二維DOA估計(jì)精度分析

將本文算法與文獻(xiàn)[12]和文獻(xiàn)[16]中算法進(jìn)行比較,兩算法都采用了兩個(gè)相互平行且相同的子陣,前者為ACA,后者為CACIS。而為了能夠更加清晰地比較算法間的性能,本文算法分別采用上述兩種陣列進(jìn)行實(shí)驗(yàn),且三種算法的子陣1陣元數(shù)目設(shè)置相同。

首先,選取多個(gè)特殊值,直觀比較三種算法的連續(xù)自由度,結(jié)果如表2和表3所示。分析可知,新算法以更少的初始陣元獲得了更高的連續(xù)自由度值,從而節(jié)約了陣元成本。

表2 本文算法和文獻(xiàn)[12]算法的連續(xù)自由度比較

表3 本文算法和文獻(xiàn)[16]算法的連續(xù)自由度比較

為比較三種算法的高分辨率性能,進(jìn)行了單次估計(jì)實(shí)驗(yàn)。假設(shè)入射信號(hào)數(shù)為2,快拍數(shù)為500,信噪比為10 dB,假設(shè)文獻(xiàn)[12]算法的陣列參數(shù)為(3,5)2,文獻(xiàn)[16]算法的陣列參數(shù)為(6,5,3)3,并將本文算法的初始陣列分別設(shè)定為上述兩陣列參數(shù),實(shí)驗(yàn)結(jié)果如圖3所示。可以發(fā)現(xiàn),雖然有一定誤差,但三種算法都能清晰地分辨出兩個(gè)較為接近的角度,且此時(shí)新算法用到的初始陣元數(shù)目最少,相較另外兩算法有著明顯優(yōu)勢(shì)。

圖3 DOA高分辨率性能比較

考慮到單次實(shí)驗(yàn)存在隨機(jī)性,為了進(jìn)一步比較三種算法的真實(shí)性能,下面設(shè)置蒙特卡洛實(shí)驗(yàn),采用均方誤差(Root Mean Square Error,RMSE)準(zhǔn)則,定義為

(31)

圖4 不同信噪比算法性能比較

接著,分析快拍數(shù)對(duì)兩種算法的影響。假設(shè)文獻(xiàn)[12]算法的陣列參數(shù)為(4,3)2,文獻(xiàn)[16]算法的陣列參數(shù)為(6,5,2)3,顯然兩者的初始陣元數(shù)目都為20,信噪比為0 dB,信源數(shù)及其角度設(shè)置和上文信噪比實(shí)驗(yàn)中完全一致,結(jié)果如圖5所示。觀察可知,本文算法的估計(jì)精度在快拍數(shù)較小時(shí),相比另外兩算法有了較大的提升,這論證了本文算法在增加連續(xù)虛擬陣元數(shù)目上的顯著性能。

圖5 不同快拍數(shù)算法性能比較

最后,比較兩種算法的運(yùn)算復(fù)雜度。由文獻(xiàn)[12]可知其算法復(fù)雜度為O{2LY+2K3+(2MN+2M)3},其中L為兩子陣列陣元數(shù)目。由文獻(xiàn)[16]可知其算法復(fù)雜度為O{2(M+N-1)2Y+(2M(N+1))3}。而由于初始陣列排布的不同,新算法下的子陣2的排布也不同,從而將直接影響本文算法復(fù)雜度,故這里不易于直接比較。為此,進(jìn)行500次實(shí)驗(yàn),快拍數(shù)為500,信源與上文蒙特卡洛實(shí)驗(yàn)一致,比較三種算法的運(yùn)行時(shí)間,結(jié)果如表4所示。

表4 500次實(shí)驗(yàn)運(yùn)行時(shí)間比較

可以發(fā)現(xiàn),本文算法在所用運(yùn)算時(shí)間上也明顯比另外兩算法少。而且由表1可知隨著子陣1陣元數(shù)目的增加,本文算法在總陣元數(shù)目上相比另外兩算法也會(huì)減少得更為顯著,它顯然也會(huì)直接反映為運(yùn)行時(shí)間上的減少,也就是說(shuō)隨著陣元數(shù)目的增加,本文算法在節(jié)約運(yùn)算成本上的能力會(huì)體現(xiàn)得更為顯著。

4 結(jié)束語(yǔ)

本文針對(duì)采用稀疏陣列進(jìn)行二維DOA估計(jì)時(shí)由于出現(xiàn)孔洞而需要舍棄一些信息,導(dǎo)致其性能不佳這一缺點(diǎn),提出了一種低復(fù)雜度基于平行稀疏陣列虛擬孔洞填充的二維DOA估計(jì)算法。該算法最直觀地優(yōu)勢(shì)在于能以更少的物理陣元數(shù)目獲得更大的連續(xù)虛擬陣元數(shù)目,從而為后續(xù)的估計(jì)工作提供更多的有用信息,極大提升了陣元資源的利用率。同時(shí),陣元數(shù)目的減少也在一定程度上降低了運(yùn)算復(fù)雜度,節(jié)約了運(yùn)算成本。另外,本文舍棄傳統(tǒng)的二維角度譜峰搜索思路,將二維問(wèn)題轉(zhuǎn)化為兩個(gè)一維問(wèn)題,并借用Root-MUSIC算法思想進(jìn)一步壓縮了算法的運(yùn)算復(fù)雜度。關(guān)鍵的是,本文算法推廣性較好,對(duì)于幾乎所有會(huì)產(chǎn)生虛擬孔洞的稀疏線性陣列,它都能在不同程度上提升對(duì)于物理陣元的資源利用率。

但是,其中一子陣陣元數(shù)目的減少,會(huì)影響用于角度β估計(jì)的兩子陣互協(xié)方差矩陣在有限快拍下的估計(jì)精度,若單論該角的估計(jì)效果是不如同算法下設(shè)置兩子陣陣元數(shù)目相同時(shí)的情況,這也是本算法的局限性,并將會(huì)在未來(lái)的研究中著重探討。