高曉峰， 栗蘋， 李國林， 郝新紅， 賈瑞麗

(北京理工大學(xué) 機電動態(tài)控制重點實驗室， 北京 100081)

0 引言

波達(dá)方向(DOA)估計作為陣列信號處理技術(shù)中的重要領(lǐng)域，在無線通信、水下探測[1-2]、聲探測[3-4]、雷達(dá)探測[5]等領(lǐng)域具有廣泛應(yīng)用。近年來，不同陣列結(jié)構(gòu)的二維DOA估計研究越來越受到關(guān)注，如雙平行陣[6-7]、圓陣[8-9]、十字陣[10-11]、矩形陣[12-13]、L型陣[14-19]等。與其他陣列相比，L型陣因其具有更好的估計性能和更低的克拉美羅界[20]成為二維DOA估計的研究熱點。文獻(xiàn)[14]提出將傳播算子方法應(yīng)用到L型陣列上進(jìn)行二維DOA估計,避免了對協(xié)方差矩陣進(jìn)行奇異值分解和特征值分解。文獻(xiàn)[15]提出利用L型陣列互相關(guān)矩陣的二維DOA估計方法，利用互相關(guān)矩陣消除了噪聲影響。文獻(xiàn)[16]提出一種基于L型陣列的互相關(guān)矩陣聯(lián)合奇異值分解的DOA估計方法，利用互相關(guān)矩陣實現(xiàn)角度自動匹配。文獻(xiàn)[17]提出一種基于L型陣列的二維盲DOA估計方法，將不同方向的陣列分解出滿足子空間旋轉(zhuǎn)不變性的子陣，在不同方向分別利用旋轉(zhuǎn)不變子空間技術(shù)(ESPRIT)算法進(jìn)行DOA估計。文獻(xiàn)[18]提出利用L型陣列不同方向的陣元信號進(jìn)行互協(xié)方差運算，利用互協(xié)方差矩陣構(gòu)建兩個互為共軛的矩陣，對矩陣進(jìn)行奇異值分解，實現(xiàn)二維DOA估計。文獻(xiàn)[19]提出一種基于L型陣列的低復(fù)雜度二維DOA估計方法，將互相關(guān)矩陣劃分為不同子矩陣，利用滿足旋轉(zhuǎn)不變性的傳播算子方法進(jìn)行二維DOA估計。上述基于L型陣列的DOA估計算法大多為針對均勻陣列的子空間類算法，需要對協(xié)方差矩陣采用特征值或奇異值分解等信號處理方法劃分信號子空間和噪聲子空間，估計性能受信噪比的影響較大，在低信噪比環(huán)境下算法估計精度較差。均勻陣列的角度分辨率與陣元孔徑呈反比，可估計的空間信源數(shù)小于陣元數(shù)，因此在陣元數(shù)有限時均勻陣列受限于陣列孔徑和陣元個數(shù)，難以獲得較高的DOA估計精度，可估計的空間信源數(shù)較少。

與均勻陣相比，稀疏陣在陣元數(shù)相同情況下具有更大的孔徑、更小的陣元間互耦效應(yīng)，因此獲得了學(xué)術(shù)界的廣泛關(guān)注，如最小冗余陣列[21]、互質(zhì)陣列[22]和嵌套陣列[23]等。與其他稀疏陣相比，嵌套陣列易于構(gòu)造，產(chǎn)生的虛擬陣列具有固定的解析表達(dá)式，利用接收信號的協(xié)方差矩陣由O(M)個陣元(M為陣元數(shù))能夠獲得O(M2)的自由度，具有較好的DOA估計效果。

本文針對上述L型均勻陣列在二維DOA估計中的問題，提出一種基于互協(xié)方差的L型嵌套陣列二維DOA估計算法。采用由兩個相互正交的2階嵌套陣構(gòu)成的L型陣列，利用不同子陣信號進(jìn)行互協(xié)方差運算，得到互協(xié)方差矩陣，并基于不同子陣噪聲不相關(guān)、互協(xié)方差矩陣消除了噪聲的影響；對互協(xié)方差矩陣進(jìn)行列向量化、產(chǎn)生無冗余陣元的虛擬陣列，利用虛擬陣列及其共軛矩陣構(gòu)建滿秩的Toeplitz矩陣作為等效協(xié)方差矩陣；利用等效協(xié)方差矩陣之間的子空間旋轉(zhuǎn)不變性，通過ESPRIT算法實現(xiàn)了目標(biāo)的方位估計；基于虛擬陣列等效信源向量的唯一性，通過估計虛擬陣列的等效信源向量實現(xiàn)空間信源的角度匹配。與L型均勻陣列相比，本文算法通過互協(xié)方差矩陣產(chǎn)生了更多的虛擬陣元，獲得了更高的角度自由度，在低信噪比環(huán)境下具有更高的估計性能，但也不可避免地增加了計算復(fù)雜度，在實際應(yīng)用中針對具體應(yīng)用背景需要加以考慮。

符號說明：(·)T、(·)H、(·)-1分別表示矩陣轉(zhuǎn)置、共軛轉(zhuǎn)置、求逆；diag(·)表示對角矩陣；∏N表示N×N階的反對角單位矩陣；⊙表示Khatri-Rao積；vec(·)表示矩陣列向量化；λ表示來波波長。

1 陣列接收信號模型

陣列模型如圖1所示：L型嵌套陣由兩個相互正交的2階嵌套陣構(gòu)成，x軸天線陣列和y軸天線陣列均由陣元數(shù)為N-1、陣元間隔為d=λ/2的均勻線陣，以及陣元數(shù)為N、陣元間隔為Nd的均勻線陣構(gòu)成；x軸天線陣列的陣元位置為p=[ξd,0,0]，y軸天線陣列的陣元位置為q=[0,ξd,0]，ξ={0,1,…,N-2,N-1,…,N2-1}，原點處公共陣元為參考陣元。

圖1 L型嵌套陣列結(jié)構(gòu)圖Fig.1 Array configuration of L-shaped nested array

圖1中：si(t)表示第i個入射信號，i=1，…，k，k為獨立遠(yuǎn)場窄帶信號的個數(shù)；θ=[θ1,θ2,…,θk]表示入射信號的俯仰角；φ=[φ1,φ2,…,φk]表示入射信號的方位角；α=[α1,α2,…,αk]表示入射信號與x軸之間的夾角；β=[β1,β2,…,βk]表示入射信號與y軸之間的夾角。設(shè)空間存在的k個獨立遠(yuǎn)場窄帶信號入射到陣列上，陣列輸出噪聲為0均值、σ2方差的高斯白噪聲。為了方便計算，在如下公式中用(α,β)代替(θ,φ)，角度轉(zhuǎn)換關(guān)系如(1)式所示：

cosα=cosθsinφ，
cosβ=sinθsinφ.

(1)

x軸天線陣列接收到的信號向量可以表示為

X(t)=[x1(t),x2(t),…,x2N-1(t)]T=
Ax(α)S(t)+Nx(t)，

(2)

式中：Ax(α)=[ax(α1),ax(α2),…,αx(αi)，…,ax(αk)]為x軸陣列的方向矩陣，其中

ax(αi)=[1,…,e-jπ(N-1)cos αi,…,e-jπ(N2-1)cos αi]T；

(3)

S(t)=[s1(t),s2(t),…,sk(t)]T為入射信號源向量；Nx(t)為x軸陣列接收到的高斯白噪聲向量。

y軸天線陣列接收到的信號向量可以表示為

Y(t)=[y1(t),y2(t),…,y2N-1(t)]T=
Ay(β)S(t)+Ny(t)，

(4)

式中：Ay(β)=[ay(β1),ay(β2),…,ay(βi)，…,ay(βk)]為y軸陣列的方向矩陣，其中

ay(βi)=[1,…,e-jπ(N-1)cos βi,…,e-jπ(N2-1)cos βi]T；

(5)

Ny(t)為y軸陣列接收到的高斯白噪聲向量。

2 基于L型嵌套陣的二維DOA估計算法

2.1 基于互協(xié)方差的虛擬陣列生成

2.1.1x軸虛擬陣列生成

x軸陣列劃分為兩個子陣,子陣1為x軸陣列前N個陣元組成的陣元間隔為λ/2的均勻線陣，其接收信號為

X1(t)=[x1(t),x2(t),…,xN(t)]T=
Ax1(α)S(t)+Nx1(t)，

(6)

式中：Nx1(t)為子陣1接收到的噪聲向量；Ax1(α)=[ax1(α1),ax1(α2),…,ax1(αi)，…，ax1(αk)]為子陣1的方向矩陣，其中

ax1(αi)=[1,e-jπcos αi，…,e-jπ(N-1)cos αi]T.

(7)

子陣2為x軸上第N個陣元到第2N-1個陣元組成的陣元間隔為Nλ/2的均勻線陣，其接收信號為

X2(t)=[xN(t),xN+1(t),…,x2N-1(t)]T=
Ax2(α)S(t)+Nx2(t)，

(8)

式中：Nx2(t)為子陣2接收到的噪聲向量；Ax2(α)=[ax2(α1),ax2(α2),…,ax2(αi),…,ax2(αk)]為子陣2的方向矩陣，其中

ax2(αi)=
[e-j(N-1)πcos αi,e-j(2N-1)πcos αi,…,e-j(N2-1)πcos αi]T.

(9)

子陣1的接收信號逆向排序，

X1z(t)=∏NX1(t)=
[xN(t),xN-1(t),…,x1(t)]T=
Ax1z(α)S(t)+Nx1z(t)，

(10)

式中：Ax1z(α)=[ax1z(α1),ax1z(α2),…,ax1z(αi),…,ax1z(αk)]，ax1z(αi)=[e-jπ(N-1)cos αi，e-jπ(N-2)cos αi,…,1]T。

逆向排序后的子陣1接收信號和子陣2接收信號進(jìn)行互協(xié)方差運算，得到互協(xié)方差矩陣Rcx：

(11)

互協(xié)方差矩陣Rcx向量化，按列拉伸成1個N2×1階的長向量，表達(dá)式如下：

(12)

(13)

2.1.2y軸虛擬陣列生成

y軸上的陣元劃分為兩個子陣,子陣3為y軸上前N個陣元組成的陣元間隔為λ/2的均勻線陣，子陣4為y軸上第N個陣元到第2N-1個陣元組成的陣元間隔為Nλ/2的均勻線陣，其接收信號為

Y1(t)=[y1(t),y2(t),…,yN(t)]T=
Ay1(β)S(t)+Ny1(t)，

(14)

Y2(t)=[yN(t),yN+1(t),…,y2N-1(t)]T=
Ay2(β)S(t)+Ny2(t)，

(15)

式中：Ny1(t)、Ny2(t)分別為子陣3和子陣4接收到的噪聲向量；Ay1(β)=[ay1(β1),ay1(β2),…,ay1(βi),…,ay1(βk)]為子陣3的方向矩陣，

ay1(βi)=[1,e-jπcos βi，…,e-j(N-1)πcos βi]T；

(16)

接收信號Ay2(β)=[ay2(β1),ay2(β2),…,ay2(βi),…,ay2(βk)]為子陣4的方向矩陣，

ay2(βi)=
[e-j(N-1)πcos βi,e-j(2N-1)πcos βi,…,e-j(N2-1)πcos βi]T.

(17)

子陣3的接收信號逆向排序得到Y(jié)1z(t)，將Y1z(t)與子陣4的接收信號Y2(t)進(jìn)行互協(xié)方差運算，得到互協(xié)方差矩陣Rcy. 由于不同子陣噪聲不相關(guān)，互協(xié)方差矩陣中的噪聲項被消除。

(18)

式中：Ay1z(β)=[ay1z(β1),ay1z(β2),…,ay1z(βi),…,ay1z(βk)]，ay1z(βi)=[e-j(N-1)πcos βi,e-j(N-2)πcos βi,…,1]T.

互協(xié)方差矩陣Rcy向量化，按列拉伸成1個N2×1階的長向量，表達(dá)式如下：

(19)

(20)

2.2 基于共軛矩陣的協(xié)方差矩陣構(gòu)造

定義x軸的虛擬陣列接收信號為

(21)

(22)

(23)

(24)

針對X1和X2，利用ESPRIT算法求解角度α，構(gòu)造如下矩陣：

(25)

對矩陣Cx進(jìn)行奇異值分解，得到k個大特征值對應(yīng)的特征向量構(gòu)成的矩陣Ex：

(26)

式中：T為k×k階的滿秩矩陣。利用虛擬陣列信號的旋轉(zhuǎn)不變性可得Ωx(α)：

(27)

由(27)式可知，對Ωx(α)進(jìn)行特征值分解可得特征值u，即可得到角度，

=arccos(arg(u)/π).

(28)

同理，定義y軸的虛擬陣列接收信號為ry:

(29)

(30)

(31)

同理，針對Y1和Y2，利用ESPRIT算法求解角度β，構(gòu)造如下矩陣：

(32)

由(26)式、(27)式，同理對矩陣Cy進(jìn)行奇異值分解，可以求得Ωy(β)，對Ωy(β)進(jìn)行特征值分解可得特征值v，即可得到角度：

=arccos(arg(v)/π).

(33)

2.3 角度匹配

由(28)式、(33)式可知，本文所提算法的兩個角度估計值和是通過兩次一維DOA估計得到的，角度和擔(dān)的估計值是失配的。本文提出利用虛擬陣列等效信源向量的唯一性進(jìn)行角度匹配。采用(28)式得到的角度重構(gòu)x軸方向虛擬陣列的方向矩陣由(12)式可知信源向量p可由陣列方向矩陣和虛擬陣列信號rx估計得到。通過最小二乘法可得信源向量的估計值1：

(34)

(35)

(36)

(37)

(38)

=Ц.

(39)

本文所提算法具體計算步驟如下：

1)x軸和y軸上的陣元分別劃分為兩個子陣，其接收信號分別為X1、X2和Y1、Y2，將X1和Y1逆向排序得到X1z、Y1z，將X1z和X2進(jìn)行互相關(guān)運算得到Rcx，將Y1z和Y2進(jìn)行互相關(guān)運算得到Rcy；

2) 由(12)式、(19)式，利用互協(xié)方差矩陣Rcx和Rcy向量化產(chǎn)生虛擬陣列信號rx和ry，等效信號源為全相干實包絡(luò)信號；

3) 利用虛擬陣列信號及其共軛矩陣構(gòu)造滿秩的Toeplitz矩陣，作為等效協(xié)方差矩陣實現(xiàn)解相干，基于ESPRIT算法得到目標(biāo)角度的估計值和；

2.4 算法分析

常規(guī)嵌套陣列利用自協(xié)方差矩陣產(chǎn)生的虛擬陣列存在冗余陣元，需要對產(chǎn)生的虛擬陣列進(jìn)行去冗余處理。去冗余后的虛擬陣列信號為相干等價信號，需要進(jìn)行解相干處理。對于M個陣元的二級嵌套陣列，采用空間平滑[23]、矩陣重構(gòu)等方法進(jìn)行解相干處理后，可估計的最大信源數(shù)Kmax為

(40)

根據(jù)(13)式、(20)式可知，本文所提算法利用M個陣元二級嵌套陣不同子陣信號的互協(xié)方差矩陣生成了無冗余陣元的虛擬陣列，消除了噪聲。針對虛擬陣列等效信號轉(zhuǎn)化為相干信號，本文利用虛擬陣列及其共軛矩陣構(gòu)造滿秩的等效協(xié)方差矩陣實現(xiàn)了解相干，由[23]式、[24]式、[30]式、[31]式可知，可估計的最大信源數(shù)為M2/4+M/2-3/4，與利用自協(xié)方差矩陣的嵌套陣列最大可估計信源數(shù)相同，優(yōu)于文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]中可估計的最大信源數(shù)M-1.

通過以上分析可知，本文所提算法的主要計算復(fù)雜度包括互協(xié)方差矩陣運算、矩陣奇異值分解、最小二乘法，所提算法涉及到的復(fù)數(shù)乘法次數(shù)為O[(M+1)2L/2+4((M+1)2/4-1)3+3k2(M+1)+3k3，其中k表示信源數(shù)，L表示快拍數(shù)。文獻(xiàn)[16]利用互協(xié)方差矩陣進(jìn)行聯(lián)合奇異值分解(JSVD)算法需要波束形成進(jìn)行角度搜索，假設(shè)搜索范圍為[0°,180°]，搜索步長為0.01°，其算法復(fù)雜度為O[M2L+8M3+18 000M2]；文獻(xiàn)[18]中所需的復(fù)數(shù)乘法次數(shù)為O[M2L+3M3+2M4]；文獻(xiàn)[19]中所需的復(fù)數(shù)乘法次數(shù)為O[M2L+2k3+(7M-4)k2+k(M-1)(2M-k)]。

圖2所示為本文所提算法與文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]所提算法的復(fù)雜度對比圖。從圖2中可以看出，本文所提算法無需進(jìn)行角度搜索，計算復(fù)雜度低于文獻(xiàn)[16]。但由于產(chǎn)生了更多的虛擬陣元，增大了角度自由度，本文所提算法的計算復(fù)雜度高于文獻(xiàn)[18-19]。

圖2 算法復(fù)雜度與陣元個數(shù)關(guān)系Fig.2 Algorithm complexity versus number of array elements

3 仿真實驗及性能分析

為了驗證本文方法的有效性，本文進(jìn)行系列仿真實驗，分析比較本文所提算法以及文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]中算法的估計性能，仿真實驗中陣元間隔為半波長。采用均方根誤差、角度分辨概率衡量DOA估計性能，如(41)式、(42)式所示：

(41)

(42)

為了驗證所提算法的有效性，下面通過仿真實驗分別分析所提算法的DOA估計性能、DOA估計性能與信噪比的關(guān)系，以及DOA估計性能與入射信源數(shù)的關(guān)系。

1) 仿真1.x軸方向陣元和y軸方向陣元個數(shù)均為9，子陣1～子陣4的陣元個數(shù)均為5，3個獨立窄帶遠(yuǎn)場的信號入射到陣列上，其入射信號角度(α,β)分別為(20°,10°)、(50°,40°)、(60°,50°)，快拍數(shù)為100，信噪比SNR分別為0 dB、5 dB、10 dB，進(jìn)行200次Monte Carlo仿真實驗。

圖3、圖4、圖5分別表示信噪比SNR為0 dB、5 dB、10 dB時本文所提算法的估計結(jié)果。從圖中可以看出，在信噪比分別為0 dB、5 dB、10 dB時，本文所提方法均成功實現(xiàn)了空間信源的方向估計，表明該算法可以在低信噪比環(huán)境下有效地進(jìn)行DOA估計。

圖3 SNR=0 dB的估計結(jié)果Fig.3 DOA estimated results for SNR=0 dB

圖4 SNR=5 dB的估計結(jié)果Fig.4 DOA estimated results for SNR=5 dB

圖5 SNR=10 dB的估計結(jié)果Fig.5 DOA estimated results for SNR=10 dB

2) 仿真2.x軸方向和y軸方向陣元個數(shù)均為9，子陣1～子陣4的陣元個數(shù)均為5，3個獨立窄帶遠(yuǎn)場信號入射到陣列上，其入射信號角度(α,β)分別為(25°,10°)、(50°,45°)、(60°,55°)，信噪比SNR為0～20 dB，進(jìn)行200次Monte Carlo仿真實驗。圖6所示為本文所提算法以及文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法的DOA估計均方根誤差與信噪比的關(guān)系圖。圖7所示為本文所提算法以及文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法的角度分辨概率與信噪比的關(guān)系圖。

圖6 均方根誤差與信噪比關(guān)系圖(快拍數(shù)為100)Fig.6 RMSE versus SNR (Snapshots=100)

圖7 分辨概率與信噪比關(guān)系圖(快拍數(shù)為100)Fig.7 Angular resolution probability versus SNR (Snapshots=100)

從圖6、圖7中可以看出，本文所提算法的DOA估計均方根誤差遠(yuǎn)小于文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法，而角度分辨概率分別優(yōu)于文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法，表明本文所提算法的DOA估計性能優(yōu)于文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法。當(dāng)信噪比SNR=0時，本文所提算法的均方根誤差小于1°，角度分辨概率接近90%，表明本文所提算法在低信噪比環(huán)境下具有較高的估計性能。隨著信噪比SNR的升高，本文所提算法與文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法相比，其均方根誤差變化較為平緩，表明本文所提算法對空間噪聲具有更強的魯棒性。

3) 仿真3.x軸方向和y軸方向陣元個數(shù)均為9，子陣1～子陣4的陣元個數(shù)均為5，入射到陣列上的信源數(shù)為1～10，入射角度(α,β)依次為(25°,30°)、(35°,40°)、(45°,50°)、(55°,60°)、(65°,70°)、(75°,80°)、(85°,90°)、(95°,100°)、(105°,110°)、(115°,120°)，信噪比SNR為10 dB，快拍數(shù)為100，進(jìn)行500次Monte Carlo仿真實驗。圖8所示為本文所提算法以及文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法的DOA估計均方根誤差與入射信源數(shù)的關(guān)系。

圖8 均方根誤差與入射信源數(shù)的關(guān)系圖 (快拍數(shù)為100)Fig.8 RMSE versus number of incident signals (Snapshots=100)

從圖8中可以看出，隨著入射信源數(shù)從1增加到4，文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法與本文所提算法的DOA估計均方根誤差均增大。在入射信源數(shù)大于5、小于8時，文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法的可估計最大信源數(shù)為M-1=8，此時入射信源數(shù)接近可估計最大信源數(shù)，算法的DOA估計均方根誤差不隨入射信源數(shù)的增大而增大。而本文所提算法的可估計最大信源數(shù)為M2/4+M/2-3/4=24，遠(yuǎn)大于仿真中入射信源數(shù)的變化范圍，其DOA估計均方根誤差隨著入射信源數(shù)的增多緩慢增大，且本文所提算法的DOA估計均方根誤差明顯小于文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法，表明隨著信源數(shù)的增多，本文所提算法的DOA估計性能優(yōu)于文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]中的算法。當(dāng)信源數(shù)為9、10時，入射信源數(shù)已經(jīng)大于文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法的可估計最大信源數(shù)M-1=8，此時它們均已失效，而本文所提算法依然能夠正常工作。

4 結(jié)論

本文針對L型均勻陣列角度分辨率較低、可估計空間信源數(shù)受限于陣元數(shù)及估計精度易受到信噪比影響的問題，提出一種新的L型嵌套陣二維DOA估計算法。該方法利用不同子陣信號的互協(xié)方差矩陣產(chǎn)生虛擬陣列并消除了噪聲的干擾，利用虛擬陣列輸出信號及其共軛矩陣構(gòu)建了等效協(xié)方差矩陣，通過ESPRIT算法分別獲得了角度α和β的估計值，利用虛擬陣列等效信源的唯一性實現(xiàn)了角度匹配。通過仿真實驗，將本文所提算法與文獻(xiàn)[16]、文獻(xiàn)[18]、文獻(xiàn)[19]算法進(jìn)行了對比，仿真結(jié)果表明本文所提算法的估計性能明顯優(yōu)于其他算法，具有更大的角度自由度，可以辨識更多的空間信源，在低信噪比環(huán)境下具有較高的DOA估計精度。本文所提算法在獲得更大角度自由度和更高估計精度的同時，不可避免地提高了算法的計算復(fù)雜度，在實際應(yīng)用過程中需加以考慮。