王 寧, 呂曉德, 李苗苗,3

(1. 中國科學(xué)院空天信息創(chuàng)新研究院, 北京 100094; 2. 微波成像技術(shù)國家級重點實驗室, 北京 100190;3. 中國科學(xué)院大學(xué)電子電氣與通信工程學(xué)院, 北京 100049)

0 引 言

波達(dá)方向(direction of arrival,DOA)估計是指從陣列中多個接收天線采集的數(shù)據(jù)中獲取電磁波信號方向信息的過程。它是陣列信號處理中的一個重要問題,廣泛應(yīng)用于雷達(dá)、聲納、通信、電子對抗、語音信號處理等領(lǐng)域[1-3]。DOA是信號模型中的非線性參數(shù),其均方根誤差(root mean square error, RMSE)曲線會產(chǎn)生閾值效應(yīng)。當(dāng)信噪比或快拍數(shù)量低于某個閾值時,容易出現(xiàn)偏離真實值的全局誤差,RMSE迅速增大并偏離Cramer-Rao 下界(Cramer-Rao lower bound, CRLB)。這種現(xiàn)象是許多非線性參數(shù)估計問題的共同特征[1]。在實際應(yīng)用中,通常會選擇閾值較低的算法,以適應(yīng)較低的信噪比和較少的快拍。從估計問題的模型來看,DOA估計方法可以分為3類:參數(shù)化法、非參數(shù)化法和半?yún)?shù)化法[4]。

參數(shù)化方法往往將信號模型建立為參數(shù)的已知函數(shù),在一定的最優(yōu)準(zhǔn)則下構(gòu)造優(yōu)化問題,通過優(yōu)化方法直接得到參數(shù)的估計值。最經(jīng)典的參數(shù)化方法是非線性最小二乘法(nonlinear least squares, NLS),可以從最大似然(maximum likelihood, ML)估計的角度推導(dǎo)出來[1]。當(dāng)有足夠的快拍時,ML估計漸近有效,漸近無偏并達(dá)到CRLB[5],因此NLS具有良好的統(tǒng)計特性。但由于NLS是非線性非凸問題,優(yōu)化過程不能保證全局最優(yōu)值,計算量大。多重信號分類(multiple signal classification, MUSIC)是一種基于協(xié)方差矩陣的特征值分解(eigenvalue decomposition, EVD)的參數(shù)化方法。為了利用導(dǎo)向矢量和噪聲子空間之間的正交性來估計DOA,將特征向量分為噪聲子空間和信號子空間兩個子集。常用的方法有MUSIC、Root-MUSIC和最小模法[6]。當(dāng)快拍數(shù)足夠且源不相關(guān)時,這些算法可以接近 ML 的統(tǒng)計特征;但當(dāng)快拍數(shù)較少且信噪比較低時,這些算法無法很好地區(qū)分噪聲子空間和信號子空間,降低了DOA估計的準(zhǔn)確性,且Root-MUSIC不能用于非冗余陣列。 基于旋轉(zhuǎn)不變技術(shù)的信號參數(shù)估計(estimating signal parameter via rotational invariance techiques, ESPRIT)也是一種基于EVD的DOA估計方法,需要兩個相同的子陣,利用子陣之間相位的旋轉(zhuǎn)不變性來估計DOA[1,6]。與 MUSIC 一樣,估計協(xié)方差矩陣也需要足夠的快拍。參數(shù)化法可以直接得到DOA,而非參數(shù)化法需要形成譜,通過找到譜的最大值來間接確定DOA。較常用的非參數(shù)方法包括傳統(tǒng)的波束形成和Capon方法。傳統(tǒng)的波束形成器基于傅里葉分析,性能受到低分辨率和頻譜泄漏的限制。Capon方法與 MUSIC 和 ESPRIT方法一樣,是基于協(xié)方差矩陣的方法,需要足夠的快拍來估計協(xié)方差矩陣。

半?yún)?shù)方法也可以稱為稀疏方法，其充分利用了協(xié)方差矩陣的Toeplitz Hermitian 特性,可用于少快拍、單快拍、信噪比低的情況，還具有高分辨率特性。在過去的10年里,半?yún)?shù)方法受到了學(xué)者們的關(guān)注并得到了快速的發(fā)展。半?yún)?shù)方法可以分為3類[7]:網(wǎng)格、離網(wǎng)和無網(wǎng)格。 網(wǎng)格類方法對 DOA 值域進(jìn)行網(wǎng)格劃分,網(wǎng)格的數(shù)量遠(yuǎn)大于實際的源數(shù)量。假設(shè)DOA的真實值在某個網(wǎng)格中,從而可以以稀疏表示的形式構(gòu)建信號模型，然后可以使用稀疏恢復(fù)算法來估計 DOA。常用的網(wǎng)格類方法包括平方根最小絕對值收斂和選擇算子(square-root least absolute shrinkage and selection operator, square-root LASSO)[8]、迭代最小化稀疏學(xué)習(xí)(sparse learning via iterative minimization,SLIM)[9]、稀疏貝葉斯學(xué)習(xí)(sparse Bayesian learning,SBL)[10]和基于協(xié)方差的稀疏迭代估計(sparse iterative covariance-based estimation,SPICE)[4,11-13]。由于DOA的真實值可能不在網(wǎng)格上,因此網(wǎng)格類方法存在網(wǎng)格誤差,限制了估計的準(zhǔn)確性。離網(wǎng)類方法和網(wǎng)格類方法都需要對DOA 值域進(jìn)行網(wǎng)格劃分,但是在構(gòu)建模型時,離網(wǎng)類方法通過引入網(wǎng)格偏移或使用動態(tài)網(wǎng)格來克服網(wǎng)格誤差,例如稀疏總體最小二乘(sparse total least-squares,STLS)[14]、離網(wǎng)的稀疏貝葉斯推斷(off-grid sparse Bayesian inference,OGSBI)[15]、網(wǎng)格進(jìn)化[16]。然而,離網(wǎng)類方法往往會引入更多的參數(shù),使算法更加復(fù)雜,而且大多數(shù)算法都涉及非凸優(yōu)化,無法保證全局收斂。文獻(xiàn)[17]采用“estimate and subtract”的策略提出了一種簡潔的離網(wǎng)類方法,可用于非均勻線性陣列(non-uniform linear array,NLA),但在中等信噪比時會出現(xiàn)閾值效應(yīng)。無網(wǎng)格類方法不需要網(wǎng)格,可以直接在 DOA 的連續(xù)域上進(jìn)行估計,并且大多數(shù)算法都是凸優(yōu)化問題,可以保證全局收斂。無網(wǎng)格類方法通過凸優(yōu)化獲得協(xié)方差矩陣的估計,然后使用Root-MUSIC或范德蒙分解計算DOA,但其中許多僅限于均勻線陣。無網(wǎng)格方法可以分為兩類:基于原子范數(shù)理論的方法[18-21]和基于協(xié)方差匹配準(zhǔn)則的方法[22-24]。文獻(xiàn)[20-21]將無網(wǎng)格類方法擴(kuò)展到任意一維陣列,不過文獻(xiàn)[21]涉及非凸優(yōu)化問題,不能保證全局收斂。

本文主要對原有的無網(wǎng)格SPICE進(jìn)行改進(jìn),使其在較低信噪比和非冗余陣列下具有更好的性能。為了降低閾值效應(yīng)的閾值,本文對Root-MUSIC進(jìn)行了修改。新的 Root-MUSIC 使用了更多的特征向量,獲得了比實際 DOA 數(shù)量更多的 DOA 估計值,并根據(jù) ML 準(zhǔn)則從中選擇合適的值作為最終估計結(jié)果。同時,本文還在原始無網(wǎng)格SPICE的優(yōu)化問題中加入了最大熵準(zhǔn)則,可以使非冗余陣列情況下的RMSE曲線更接近CRLB,同時保持優(yōu)化問題的凸性。將所提出方法的性能與 MSUIC、用于NLA的無網(wǎng)格SPICE[23]、交替投影無網(wǎng)格(alternating projections gridless, AP gridless)[21]和 CRLB 進(jìn)行了比較。蒙特卡羅實驗表明,本文提出的方法具有較低的閾值,可以在較低信噪比的環(huán)境中使用,并且還具有更好的高分辨率特性。同時,本文所提出的方法可以在相干源和信源數(shù)大于陣元數(shù)的情況下使用。

本文的后續(xù)內(nèi)容組織如下:第1節(jié)介紹了經(jīng)典的遠(yuǎn)場窄帶模型、陣列信號處理中的冗余和非冗余陣列;第2節(jié)介紹了具有最大熵和基于ML-Root-MUSIC 的無網(wǎng)格 SPICE;第3節(jié)通過仿真驗證了所提出的方法在低信噪比和非冗余陣列的情況下的有效性;第4節(jié)進(jìn)行了總結(jié)。

1 信號模型

假設(shè)K個遠(yuǎn)場窄帶平面波照射到M元均勻線陣,DOA記為θ=[θ1,θ2,…,θK]T,(θk∈(-90°,90°]),本文假設(shè)K是已知的。那么陣列接收到的數(shù)據(jù)可以寫為

A(θ)s[n]+e[n],n=1,2,…,N

(1)

式中:n是快拍的序號;N是快拍數(shù);y[n]∈CM是n時刻快拍數(shù)據(jù);a(θk)∈CM是第k個信號的導(dǎo)向矢量;sk[n]∈C是第k個信號在n時刻的數(shù)據(jù);e[n]∈CM是噪聲;A(θ)=[a(θ1),…,a(θK)]被稱為陣列流形。

可以用更緊湊的方式表示為

Y=A(θ)S+E

(2)

式中：Y=[y[1],y[2],…,y[N]]∈CM×N;S∈CK×N;E∈CM×N。本文對信號模型的統(tǒng)計特性做了如下假設(shè):

(1) 源信號均值為零,并且在空域和時域中不相關(guān),即

E[s[i]sH[j]]=diag(p)δi, j

(3)

(2) 噪聲向量的均值為零并且在空域和時域中都是白色,即

E[e[i]eH[j]]=σ2Iδi, j

(4)

(3) 噪聲向量和源信號不相關(guān),即

E[s[i]eH[j]]=0

(5)

式中:p=[p1,p2,…,pK];pk∈R+表示源信號的功率;σ2表示噪聲方差;δi, j是離散沖激函數(shù)。接收數(shù)據(jù)的協(xié)方差矩陣可以表示為

R=E[y[n]yH[n]]=

A(θ)diag(p)AH(θ)+σ2I

(6)

如果陣列為均勻線陣,則陣元的位置可表示為d=[1,2,…,M]T,距離單位為半波長。導(dǎo)向矢量可以表示為式(7)。fk稱為空間頻率,1是元素全為1的列向量。因為a(θk)有復(fù)正弦的形式,那么A(θ)就是一個范德蒙矩陣?？梢赃M(jìn)一步證明R是Toeplitz Hermitian矩陣且rank(R)=M。

fk∈(-0.5,0.5]

(7)

當(dāng)陣列為非均勻線陣時,陣元的位置可表示為d=[d1,d2,…,dL]T(dl∈Z+)。不失一般性,令d1=1,dL=M。用ad(θk)、Ad(θ)、yd[n]、Rd分別表示非均勻陣列情況下的導(dǎo)向矢量、陣列流形矩陣、快拍和協(xié)方差矩陣。由于d不再是包含連續(xù)正整數(shù)的向量,不再具有復(fù)正弦的形式,則Ad(θ)不再是范德蒙德矩陣,而Rd只是Hermitian矩陣。我們可以定義一個選擇矩陣Γd,只有第i行(i∈1,2,…,L)第di列是1,其余位置是0。然后可以得到:

yd[n]=ΓdAd(θ)s[n]+Γde[n]=Γdy[n]

(8)

Rd和R的關(guān)系如下所示:

(9)

為了顯示Rd中包含R的多少個元素,可以定義d的合成D[22,25]。如果D={1,2,…,M}, 則Rd包含R中的所有元素,d對應(yīng)的陣列就是冗余陣列。如果D不包含{1,2,…,M}中的所有元素,則Rd不包含R中的所有元素,并且d對應(yīng)的陣列是非冗余陣列。例如,當(dāng)d=[1,2,5,7]T時為冗余陣列:

D={m1-m2+1:m1,m2∈d,m1≥m2}

(10)

此時，Rd和R如下所示：

(11)

當(dāng)d=[1,3,5,7]T時,該陣列為非冗余陣列,Rd如下所示：

(12)

2 ME-GLS-ML

本節(jié)介紹在原始無網(wǎng)格SPICE的基礎(chǔ)上提出的方法,稱之為聯(lián)合最大熵和最大似然思想的無網(wǎng)格方法(maximum entropy-gridless SPICE-ML, ME-GLS-ML)。一方面,為了使無網(wǎng)格SPICE可以用于非冗余陣列,在優(yōu)化問題中引入了最大熵的思想,加入了負(fù)熵項。同時,負(fù)熵項可以約束協(xié)方差矩陣是正定的。 另一方面,為了降低閾值效應(yīng)的閾值,在Root-MUSIC中,通過ML準(zhǔn)則選擇根而不是選擇最接近單位圓的根。

2.1 Gridless SPICE

Gridless SPICE 是一種基于協(xié)方差擬合的無網(wǎng)格DOA 估計方法,使用R的Toeplitz Hermitian屬性將其參數(shù)化為r=[r0,r1,…,rM-1]T的函數(shù)R=T(r)。T(r)表示Toeplitz Hermitian矩陣,其第一行是rT,并且是rT的線性函數(shù),這避免了使用DOA作為參數(shù)形成高度非線性函數(shù)[7]。以r作為參數(shù)的優(yōu)化問題是一個凸優(yōu)化問題,很容易求解。

(13)

將R=T(r)代入式(13)后可以得到

s.t.T(r)≥0

(14)

可以看到,式(14)是一個半正定規(guī)劃問題,是一個凸優(yōu)化問題[26],可以使用Matlab中的CVX工具箱來解決[27]。 一旦得到最優(yōu)解,就可以構(gòu)造R,然后可以使用Root-MUSIC算法計算DOA,從而達(dá)到無網(wǎng)格DOA估計的目的。

(15)

(16)

(17)

2.2 加入最大熵準(zhǔn)則

從式(16)和式(17)可以看出,在非冗余陣列的情況下,優(yōu)化問題與r中的某些變量無關(guān)。例如d=[1,3,5,7]T, 式(16)和式(17)只能與r0、r2、r4、r6相關(guān),而不能與r1、r3、r5相關(guān),因此只能得到r0、r2、r4、r6的值。

非冗余陣列引起的問題可以理解為協(xié)方差矩陣的r0、r2、r4、r6已知,確定r1、r3、r5的方法如下所示。

R=T([r0,?,r2,?,r4,?,r6])

(18)

式中:？表示該處值未知。

在譜估計中,周期圖法無法估計具有高延遲的相關(guān)函數(shù)的值并將其設(shè)置為0,但突然轉(zhuǎn)變?yōu)?會導(dǎo)致頻譜中出現(xiàn)虛假的頻譜分量。為了解決這個問題,Burg方法將相關(guān)函數(shù)的高延遲設(shè)置為對數(shù)據(jù)做出最少假設(shè)的值,即最大化隨機過程的熵[28]。對于式(18)所示的問題,也可以利用最大熵的思想來獲得協(xié)方差矩陣的未知值[25]。

假設(shè)y[n]是圓復(fù)高斯的,則其熵為ln[(πe)MdetR][29-30]。將熵代入式(16)和式(17)得到：

(19)

λln[detT(r)]

(20)

最大化熵相當(dāng)于最小化負(fù)熵,所以熵項前面有一個負(fù)號。忽略與優(yōu)化問題無關(guān)且大于零的常數(shù)項。由于負(fù)熵是一個凸函數(shù),所以新生成的優(yōu)化問題(式(19)和式(20))仍然是凸優(yōu)化問題,可以使用CVX工具箱解決。負(fù)熵可以看作是約束T(r)>0的障礙函數(shù)。同時,其進(jìn)一步在最大熵的統(tǒng)計意義下限制了T(r)。

2.3 ML-Root-MUSIC

(21)

(22)

PA(θ)=I-A(θ)(AH(θ)A(θ))-1AH(θ)

(23)

在本文中,使用ML從pk個候選根中選擇K個根,如下所示：

(24)

稱這種新方法為ML-Root-MUSIC。在隨后的仿真中,發(fā)現(xiàn)ML-Root-MUSIC可以降低閾值,從而可以在較低的信噪比下應(yīng)用無網(wǎng)格SPICE。

本文所提出的方法的處理過程如算法1所示。

算法1 ME-GLS-ML輸入:觀測數(shù)據(jù)Y,信號子空間維數(shù)ps,備選根數(shù)pk,源數(shù)目K,常數(shù)λ。輸出:K個DOA。(1) 計算相關(guān)陣R^=(1/N)YYH;(2) 求解優(yōu)化問題式(19)和式(20),得到R～;(3) 對R～進(jìn)行EVD, 得到ps個最大特征值對應(yīng)的特征向量,構(gòu)成U～s, 剩余特征向量構(gòu)成U～n;(4) 用U～n構(gòu)造多項式并求解pk個最接近單位圓的根;(5) 用式(24)從pk個根中選擇K個,記為z^i(i=1,2,…,k);(6) 輸出θ^i=arcsin(argzi·π-1)(i=1,2,…,k)。

3 數(shù)值仿真

3.1 不同信噪比下的RMSE

本節(jié)通過蒙特卡羅仿真實驗比較了所提方法與其他方法在非冗余陣列情況下的估計精度。 在仿真中,設(shè)置快拍數(shù)為200,陣元的位置為[1,3,6,8,11,14]T,該陣列是非冗余的?？臻g頻率設(shè)置為[0.2,0.3,0.4]T,蒙特卡羅試驗次數(shù)為1 000。源和噪聲均使用圓復(fù)高斯白噪聲,信噪比定義為10lg(p/σ2),p為源的方差,σ2為噪聲方差,并且每個源的方差是相同的。CRLB的計算公式[1,6]為

(PAH(θ)R-1A(θ)P)T]}-1

(25)

(1) MUSIC:搜索譜峰時的網(wǎng)格數(shù)為4 096。

(2) NLA-GLS:文獻(xiàn)[23]中用于NLA的無網(wǎng)格SPICE,在本文稱為NLA-GLS。

(3) AP gridless:文獻(xiàn)[21]中的方法。

(4) NLA-GLS+:使用了ML-Root-MUSIC的NLA-GLS,設(shè)置信號子空間的維數(shù)ps=3,設(shè)置備選根數(shù)pk=8。

(5) ME-GLS-ML:本文提出的方法,設(shè)置信號子空間的維數(shù)ps=3,設(shè)置備選根數(shù)pk=8。

5種方法非冗余陣列下RMSE隨信噪比的變化如圖1所示?？梢钥闯?這5種方法在信噪比小于某個閾值后都會產(chǎn)生閾值效應(yīng)。MUSIC和AP gridless的閾值為-6 dB,NLA-GLS的閾值為-2 dB,NLA-GLS+和ME-GLS-ML的閾值為-8 dB。NLA-GLS+和ME-GLS-ML的閾值最低。與NLA-GLS相比,由于NLA-GLS+增加了ML-Root-MUSIC,閾值從-2 dB降低到-8 dB。與NLA-GLS+相比,ME-GLS-ML增加了負(fù)熵項,RMSE曲線更接近CRLB,說明在非冗余陣列情況下增加負(fù)熵項有利于提高估計精度。在這5種方法中,ME-GLS-ML的閾值最低,最貼近CRLB 曲線。因此,本文提出的方法在非冗余陣列的情況下更具優(yōu)勢。

圖1 非冗余陣列下RMSE隨信噪比的變化Fig.1 Change of RMSE with signal to noise ratio in non-redundant array

圖2是NLA-GLS和ME-GLS-ML在-8 dB信噪比下1 000次蒙特卡羅結(jié)果的直方圖。可以看出ME-GLS-ML的估計結(jié)果比NLA-GLS更集中在真實值附近,NLA-GLS存在偏離真實值的異常值。

圖2 1 000 次蒙特卡羅試驗直方圖Fig.2 Histogram of 1 000 Monte Carlo

3.2 分辨率

在本節(jié)增加如下兩個方法的仿真。

(1) NLA-GLS++:設(shè)置NLA-GLS+信號子空間的維數(shù)ps=8,設(shè)置備選根數(shù)pk=8。

(2) ME-GLS-ML+:設(shè)置ME-GLS-ML信號子空間的維數(shù)ps=8,設(shè)置備選根數(shù)pk=8。

在仿真中,設(shè)置快拍數(shù)為200,陣元的位置為[1,3,6,8,11,14]T,該陣列是非冗余的?；诟道锶~分析的DOA估計方法的空間頻率分辨率約為1/M[1,6]。一般能夠使分辨率小于1/M的方法稱為超分辨率方法。在仿真中生成兩個源信號,其中一個具有固定空間頻率0.2,另一個具有空間頻率0.2+Δf/M,Δf為頻率間隔(Δf∈{0.1,0.2,…,1})。兩個源信號功率一樣,信噪比固定為-4 dB。每個Δf分別進(jìn)行200次蒙特卡羅模試驗。當(dāng)滿足下式時,認(rèn)為是成功分辨的[1]：

(26)

仿真結(jié)果如圖3所示,所有方法均可實現(xiàn)超分辨。NLA-GLS+和ME-GLS-ML的超分辨性能并不比NLA-GLS好,但是當(dāng)信號子空間的維數(shù)ps增加到8,大于真實DOA數(shù)目時,超分辨性能得到了改善,如NLA-GLS++和ME-GLS-ML+的曲線所示。這一結(jié)果表明,適當(dāng)增加信號子空間維數(shù)有利于提高超分辨率性能。NLA-GLS++在所有比較方法中分辨性能最好,說明使用ML-Root-MUSIC可以提高分辨性能。

圖3 分辨成功率Fig.3 Success rate of resolution

3.3 相干源

本小節(jié)通過仿真說明ME-GLS-ML在源信號相干的情況下依舊可以進(jìn)行DOA估計,信噪比固定為0 dB,各源信號使用相同的圓復(fù)高斯白噪聲來模擬相干源的場景,其余仿真條件與第3.1節(jié)一致。仿真結(jié)果如圖4所示,經(jīng)典的MUSIC方法在相干源的情況下會出現(xiàn)很多偽峰,無法得到正確得DOA估計,而本文所提出的ME-GLS-ML可以正確得到DOA的估計。

圖4 相干源情況下的仿真Fig.4 Simulation with coherent sources

3.4 源數(shù)目大于陣元數(shù)

本節(jié)考慮源數(shù)目大于陣元數(shù)的情況,源數(shù)目設(shè)置為8,相應(yīng)的空間頻率分別設(shè)置為-0.4、-0.3、-0.2、-0.1、0.1、0.2、0.3、0.4,信噪比為0 dB,陣元位置的設(shè)置與第3.1節(jié)相同, 此時陣元數(shù)為6。仿真結(jié)果如圖5所示,ME-GLS-ML可以成功得到所有源的DOA估計。由于ME-GLS-ML通過優(yōu)化的方式得到了虛擬的14陣元均勻線陣所對應(yīng)的協(xié)方差矩陣,因此可以對多余真實陣元數(shù)的源進(jìn)行DOA估計,理論上最多可以對13個源進(jìn)行DOA估計。

圖5 源數(shù)目大于陣元數(shù)情況下的仿真Fig.5 Simulation with the number of sources greater than the number of array elements

4 結(jié) 論

本文主要關(guān)注DOA估計中的閾值效應(yīng),提出了一種適用于低信噪比和非冗余陣列的無網(wǎng)格DOA估計方法。該方法在原有無網(wǎng)格SPICE的基礎(chǔ)上,引入最大熵的思想,構(gòu)造了一個新的凸優(yōu)化問題。該凸優(yōu)化問題可以用來擬合非冗余陣列的協(xié)方差矩陣,然后通過ML-Root-MUSIC完成DOA估計。根據(jù)蒙特卡羅仿真實驗,該方法具有更好的性能。該方法在相同數(shù)量的快拍下具有較低的閾值,可應(yīng)用于信噪比較低的環(huán)境。同時,在相同的信噪比和相同快拍數(shù)量下,本文提出的方法具有更好的高分辨率特性。在 ML-Root-MUSIC中,使用了比實際信源數(shù)更多的特征向量,并且使用ML準(zhǔn)則進(jìn)行DOA選擇,這個想法也可以擴(kuò)展到其他基于EVD的方法。