(桂林電子科技大學(xué)信息與通信學(xué)院， 廣西桂林 541004)

0 引言

近年來(lái)，信源定位成為陣列信號(hào)處理中的一個(gè)熱點(diǎn)問(wèn)題并受到廣泛關(guān)注。大量針對(duì)遠(yuǎn)場(chǎng)條件下的波達(dá)方向(Direction of Arrival, DOA)估計(jì)算法被提出[1-6]。當(dāng)信源位于菲涅爾區(qū)時(shí)，信號(hào)不再以平面波傳播，而是以球面波的形式傳播，且信源的信息由角度和距離兩個(gè)參數(shù)共同決定，此時(shí)為近場(chǎng)源定位。為此，國(guó)內(nèi)外學(xué)者提出了大量估計(jì)方法，如2-D 多重信號(hào)分類(Multiple Signal Classification, MUSIC)算法[7]和高階旋轉(zhuǎn)不變子空間(Estimation of Signal Parameters via Rotational Invariance Techniques，ESPRIT)算法[8]都可以有效估計(jì)。

然而在實(shí)際應(yīng)用中近場(chǎng)信源和遠(yuǎn)場(chǎng)信源有可能同時(shí)存在，由于信源的不匹配，上述所提到的近場(chǎng)或遠(yuǎn)場(chǎng)信源DOA估計(jì)算法將不再適用。為了解決混合信源DOA估計(jì)問(wèn)題，文獻(xiàn)[9]通過(guò)構(gòu)造高階累積量矩陣，然后使用MUSIC算法估計(jì)近場(chǎng)信源參數(shù)，計(jì)算量太大；為了減少運(yùn)算量，文獻(xiàn)[10]提出一種基于二階統(tǒng)計(jì)量的混合信源參數(shù)估計(jì)，然而該算法在進(jìn)行近場(chǎng)DOA估計(jì)時(shí)只利用了部分協(xié)方差矩陣的信息，造成較大的陣列孔徑損失，估計(jì)精度較差；為了減少計(jì)算復(fù)雜度和陣列孔徑損失，文獻(xiàn)[11]將對(duì)稱陣均勻線陣分為2分子陣，利用子陣間的關(guān)系使用ESPRIT-like算法進(jìn)行參數(shù)估計(jì)，但估計(jì)結(jié)果仍難以令人滿意。

以上算法都是基于均勻線陣，存在著一定的陣列孔徑損失，在陣元數(shù)一定的情況下，稀疏布陣可以增加陣列孔徑，從而提高信源參數(shù)估計(jì)精度。文獻(xiàn)[12]采用互素對(duì)稱陣列，通過(guò)兩個(gè)陣元數(shù)互為素?cái)?shù)的子陣構(gòu)造一個(gè)四階累積量矩陣，然后利用MUSIC算法估計(jì)信源角度，有效擴(kuò)展了陣列孔徑。文獻(xiàn)[13]使用了一種特殊幾何陣列結(jié)構(gòu)，利用傳統(tǒng)的MUSIC算法和構(gòu)造特殊點(diǎn)的四階累積量聯(lián)合估計(jì)信源的角度和距離。為了進(jìn)一步擴(kuò)展陣列孔徑，本文提出一種稀疏對(duì)稱陣列模型。首先通過(guò)對(duì)不同子陣的接收數(shù)據(jù)進(jìn)行四階累積量運(yùn)算剔除近場(chǎng)信源的距離參數(shù)，構(gòu)造多個(gè)僅與信源角度有關(guān)的四階累積量向量，通過(guò)這些累積量向量構(gòu)造一個(gè)Topelize矩陣，再利用MUSIC算法估計(jì)出信源角度，最后在估計(jì)出角度的基礎(chǔ)上進(jìn)行距離搜索，根據(jù)近場(chǎng)遠(yuǎn)場(chǎng)信源位于不同區(qū)域估計(jì)出近場(chǎng)信源的距離參數(shù)。該算法避免了二維搜索，且稀疏布陣擴(kuò)展了陣列孔徑，提高了參數(shù)估計(jì)精度。

1 信號(hào)模型

本文所采用的陣列模型如圖1所示，陣元總數(shù)為2(N1+N2)+1,由3個(gè)子陣列組成。其中子陣1陣元數(shù)和陣元間距分別為2N1+1和d，子陣2和子陣3的陣元數(shù)和陣元間距分別為N2和(2N1+1)d,子陣列1與子陣列2和子陣列3之間的距離分別為(N1+1)d。圖中各陣元坐標(biāo)為pi，pi=-N2(2N1+1),-(N2-1)(2N1+1),…,-(2N1+1),-N1,…,0,…,N1,(2N1+1),…,(N2-1)(2N1+1),N2(2N1+1),pid為第i個(gè)陣元到中心陣元的距離。當(dāng)陣元數(shù)相同時(shí)，均勻線陣的陣列孔徑為(2N1+2N2+1)d，本文所采用陣列的陣列孔徑為(2N2+4N1N2)d，可以看出，本文的陣列具有更大的陣列孔徑。

圖1 稀疏對(duì)稱陣列模型

假設(shè)空間存在K個(gè)獨(dú)立窄帶信號(hào)入射到陣列，以中心陣元為相位參考點(diǎn)，則第i個(gè)陣元接收數(shù)據(jù)為

t=1,…,T

(1)

式中，sk(t)為第k個(gè)信號(hào)源，T為快拍數(shù)，ni(t)為第i個(gè)陣元接收到的噪聲，μk和φk分別為[14]

(2)

(3)

式中，λ為信號(hào)的波長(zhǎng)，rk，θk分別為第k個(gè)近場(chǎng)信源的距離和角度參數(shù)。

將陣元接收數(shù)據(jù)寫為矢量形式為

x(t)=As(t)+n(t)

(4)

式中,x(t)=[x-N1-N2(t),…,x-N1(t),…,xN1(t),…,xN1+N2(t)]T為陣元接收數(shù)據(jù)，s(t)=[s1(t),s2(t),…,sK(t)]T為K個(gè)信源的信號(hào)矢量，n(t)=[n-N1-N2(t),…,nN1+N2(t)]T為陣元接收的噪聲矢量，A為(2(N1+N2)+1)×K維陣列流型矢量，表示為

A=[a(θ1,r1),…,a(θk,rk),…a(θK,rK)]

(5)

式中，a(θk,rk)為(2(N1+N2)+1)×1維方向矢量，當(dāng)信源位于近場(chǎng)時(shí)，a(θk,rk)表示為

a(θk,rk)=[ej[-(2N1N2+N2)μk+(-(2N1N2+N2))2φk],…,

ej[-(2N1+1)μk+(-(2N1+1))2φk],ej[-(N1)μk+(-N1)2φk],…,1,…,

ej[(2N1N2+N2)μk+(2N1N2+N2)2φk]]

(6)

當(dāng)信源位于遠(yuǎn)場(chǎng)時(shí)，a(θk,rk)表示為

a(θk,rk)=[ej[-(2N1N2+N2)μk],…,

ej[-(2N1+1)μk],ej[-(N1)μk],…,

1,…,ej[(N1)μk],ej[(2N1+1)μk],…,

ej[(2N1N2+N2)μk]]T

(7)

本文中，對(duì)信號(hào)作如下假設(shè)：

1) 信號(hào)為零均值、非高斯的窄帶平穩(wěn)隨機(jī)過(guò)程，且具有非零峰度，信號(hào)之間不相關(guān)。

2) 陣元接收的噪聲為零均值的高斯白噪聲，并且與信號(hào)相互獨(dú)立。

3) 為了確保信源參數(shù)估計(jì)的唯一性，陣列1陣元最小間距d≤λ/4。

2 基于稀疏對(duì)稱陣列的信源定位

2.1 構(gòu)造特殊四階累積量矩陣

高階累積量具有抑制高斯噪聲，同時(shí)還可以擴(kuò)展陣列孔徑等優(yōu)點(diǎn)，因此基于高階累積量的空間譜估計(jì)得到廣泛關(guān)注。本文所采用的四階累積量公式為[15]

(8)

(9)

通過(guò)四階累積量構(gòu)造一個(gè)僅與信源角度有關(guān)的特殊矩陣，使其稀疏對(duì)稱陣等效于陣元間距為λ/4的均勻線陣列。令m∈[-N1,…,N1]，n為0，首先，將子陣1陣元的接收數(shù)據(jù)與中心處的陣元的接收數(shù)據(jù)進(jìn)行四階累積運(yùn)算得到一個(gè)(2N1+1)×1維四階累積量向量c1，其第m個(gè)元素為

c1(m+N1+1)=

(10)

同理，將子陣1的陣元接收到的數(shù)據(jù)與子陣3的第一個(gè)陣元接收的數(shù)據(jù)進(jìn)行四階累積量的運(yùn)算，構(gòu)造出(2N1+1)×1維四階累積量向量c2，其第m個(gè)元素分別為

c2(m+N1+1)=

(11)

將子陣1的陣元接收到的數(shù)據(jù)與子陣2的第一個(gè)陣元接收的數(shù)據(jù)進(jìn)行四階累積量的運(yùn)算，構(gòu)造出(2N1+1)×1維四階累積量向量c3，其第m個(gè)元素分別為

c3(m+N1+1)=

(12)

令n∈[0,…,N1]，m∈[N1+2,…,N1+N2]，將子陣3陣元接收到的數(shù)據(jù)與子陣1的參考陣元右半部分陣元接收到的數(shù)據(jù)進(jìn)行四階累積量運(yùn)算，構(gòu)造出(N1+1)(N2-1)×1維向量c4，其第 ((m-N1-2)(N1+1)+n+1)個(gè)元素為

c4((m-N1-2)(N1+1)+n+1)=

(13)

同理，將子陣1陣元接收到的數(shù)據(jù)與子陣1的參考陣元左半部分陣元接收到的數(shù)據(jù)進(jìn)行四階累積量運(yùn)算，構(gòu)造出(N1+1)(N2-1)×1維向量c5，其第((m-N1-2)(N1+1)+n+1)個(gè)元素為

c5((m-N1-2)(N1+1)+n+1)=

(14)

將向量c5，c3，c1，c2，c4壘成一個(gè)(2(N1N2+2N1+N2)+1)×1維的長(zhǎng)向量c,表示為

(15)

可以發(fā)現(xiàn)，向量c中元素的相位項(xiàng)m-n的值跟均勻線陣的相位等效，分別為：-(N1N2+2N1+N2),-(N1N2+2N1+N2-1),…,0,…,(N1N2+2N1+N2-1),(N1N2+ 2N1+N2)。通過(guò)四階累積量運(yùn)算使稀疏對(duì)稱陣列的接收數(shù)據(jù)等效為一個(gè)均勻線陣，如圖2所示。

圖2 重構(gòu)均勻線陣模型

重構(gòu)方法：

矢量c1是由子陣1內(nèi)各陣元的接收數(shù)據(jù)分別與子陣1中的中心參考陣元的接收數(shù)據(jù)做累積量運(yùn)算的結(jié)果值組合而得;矢量c2是由子陣1內(nèi)各陣元的接收數(shù)據(jù)分別與子陣3中第一個(gè)陣元的接收數(shù)據(jù)做累積量運(yùn)算的結(jié)果值組合而得，矢量c3是由子陣1內(nèi)各陣元的接收數(shù)據(jù)分別與子陣2中第一個(gè)陣元的接收數(shù)據(jù)做累積量運(yùn)算的結(jié)果值組合而得，矢量c4是由子陣1內(nèi)各陣元的接收數(shù)據(jù)分別與子陣3中除第一個(gè)陣元外的所有陣元接收數(shù)據(jù)做累積量運(yùn)算的結(jié)果值組合而得，矢量c5是由子陣1內(nèi)各陣元的接收數(shù)據(jù)分別與子陣2中除第一個(gè)陣元外的所有陣元接收數(shù)據(jù)做累積量運(yùn)算的結(jié)果值組合而得。

在矢量的重構(gòu)過(guò)程中，將矢量c1的第m個(gè)元素對(duì)應(yīng)于重構(gòu)矢量c的第m個(gè)元素，矢量c2的第m個(gè)元素對(duì)應(yīng)于重構(gòu)矢量c的第m+2N1+1個(gè)元素，矢量c3的第m個(gè)元素對(duì)應(yīng)于重構(gòu)矢量c的第-(m+2N1+1)個(gè)元素，矢量c4的第m個(gè)元素對(duì)應(yīng)于重構(gòu)矢量c的第(m-N1)(N1+1)+N1+n個(gè)元素，矢量c5的第m個(gè)元素對(duì)應(yīng)于重構(gòu)矢量c的第-((m-N1)(N1+1)+N1+n)個(gè)元素。

向量c構(gòu)造一個(gè)(N1N2+2N1+N2+1)×(N1N2+2N1+N2+1)維的Topelize矩陣C，C的第m列可以表示為

C(:,m)=c(N1N2+2N1+N2+2-

m:2(N1N2+2N1+N2)+2-m)

(16)

因此，我們構(gòu)造了一個(gè)僅包含信源角度信息的特殊矩陣C，類似于遠(yuǎn)場(chǎng)協(xié)方差矩陣，可以表示為

C=A(θ)C4,SAH(θ)

(17)

式中，

C4,S=diag[c4,s1,…,c4,sP]

(18)

A(θ)=[a(θ1),…,a(θp)]

(19)

a(θp)=[1,ej2μp,…,ej2(N1N2+2N1+N2)μp]T

(20)

針對(duì)接收端對(duì)回波數(shù)據(jù)的接收方式，與已有的將稀疏對(duì)稱陣列劃分為4個(gè)子陣來(lái)分別接收數(shù)據(jù)的方式不同[13]，本文中將稀疏對(duì)稱陣列劃分為3個(gè)子陣來(lái)分別接收數(shù)據(jù)。

2.2 信源DOA估計(jì)

由式(17)～式(20)可以看出，矩陣A(θ)類似于K個(gè)遠(yuǎn)場(chǎng)信號(hào)入射到陣元數(shù)為(N1N2+2N1+N2+1)的均勻線陣所產(chǎn)生的陣列流型矩陣，a(θk)為第k個(gè)信號(hào)的類遠(yuǎn)場(chǎng)陣列流型矢量。由于信源sk的四階累積量非零，且A(θ)和C4,S均為滿秩，因此可以運(yùn)用MUSIC算法來(lái)估計(jì)信源角度，對(duì)構(gòu)造的四階累積量矩陣C進(jìn)行特征值分解為

(21)

Vs=R(N1N2+2N1+N2+1)×K為較大K個(gè)特征值構(gòu)成的信號(hào)子空間，Vn=R(N1N2+2N1+N2+1)×((N1N2+2N1+N2+1)-K)為較小(N1N2+2N1+N2+1-K)個(gè)特征值構(gòu)成的噪聲子空間。由式(22)所示的MUSIC譜可估計(jì)出信源的方位角：

(22)

2.3 近場(chǎng)信源距離估計(jì)

陣列接收數(shù)據(jù)的協(xié)方差矩陣為

R=E[x(t)xH(t)]

(23)

(24)

如果rk∈[0.62(D3/λ),2D2/λ]，其中D為陣列孔徑，則與rk對(duì)應(yīng)的信源位于菲涅耳區(qū)，屬于近場(chǎng)信源。若rk>2D2/λ，則與rp對(duì)應(yīng)的信源為遠(yuǎn)場(chǎng)信源。由此可以估計(jì)出近場(chǎng)信源參數(shù)，無(wú)需二維搜索，且近場(chǎng)信源的角度和距離參數(shù)自動(dòng)配對(duì)。

由于采用稀疏對(duì)稱陣列接收數(shù)據(jù)的協(xié)方差矩陣來(lái)估計(jì)近場(chǎng)信源的距離參數(shù)，直接用稀疏陣列估計(jì)信源參數(shù)是需要考慮空間模糊問(wèn)題。信源θk方向上產(chǎn)生距離模糊的條件為

(25)

式中，pk為第k個(gè)信源的位置，l為整數(shù)，|l|≥1。由式(23)可得

(26)

要使式(26)右邊取得最小值，令cos2θk=1，rk=0.62(D3/λ)1/2,r′k=∞，化簡(jiǎn)得

(27)

假設(shè)陣列孔徑D=λ,最小陣元間距d=λ/4代入式(25)可得pk≥4.45。因此pk≥5或pk≤-5時(shí)陣列流型向量的第k個(gè)元素會(huì)產(chǎn)生模糊，然而當(dāng)D=λ，d=λ/4時(shí)，-4≤pk≤4，所以整個(gè)陣列流型向量都不會(huì)產(chǎn)生模糊現(xiàn)象，估計(jì)出的每一個(gè)近場(chǎng)信源的距離參數(shù)都是唯一的。

2.4 計(jì)算復(fù)雜度分析

3 計(jì)算機(jī)仿真

為了驗(yàn)證所提算法的優(yōu)良性能，將本文所提算法與文獻(xiàn)[11]中基于二階統(tǒng)計(jì)量的ESPRIT算法、文獻(xiàn)[12]中基于互素對(duì)稱陣的近場(chǎng)源定位和文獻(xiàn)[13]中基于特殊的幾何陣列結(jié)構(gòu)的遠(yuǎn)近場(chǎng)混合定位進(jìn)行對(duì)比。假設(shè)陣元總數(shù)7(N1=1，N2=2)，最小陣元間距為λ/4，采用角度和距離的均方根誤差(RMSE)作為衡量標(biāo)準(zhǔn)，分別定義為

(28)

(29)

實(shí)驗(yàn)1：信號(hào)源為近場(chǎng)遠(yuǎn)場(chǎng)混合信源。入射信號(hào)為2個(gè)4QAM窄帶信號(hào)，近場(chǎng)信源與遠(yuǎn)場(chǎng)信源位置分別為(θ1=-13°,r1=3.5λ)和(θ2=21°,r2=∞)。

在本實(shí)驗(yàn)中，驗(yàn)證本文算法DOA估計(jì)精度隨信噪比和快拍數(shù)變化的情況。第一，信噪比從-5 dB到10 dB不斷變化，步長(zhǎng)為2 dB，快拍數(shù)為700，蒙特卡洛次數(shù)為500；第二，信噪比為7 dB，快拍數(shù)從100到1 100不斷變化，步長(zhǎng)為200。圖3為信源角度均方根誤差隨信噪比的變化情況，圖4為近場(chǎng)信源距離的均方根誤差隨信噪比的變化。圖5和圖6分別是方位角和近場(chǎng)信源距離隨快拍數(shù)的變化曲線圖。從圖3可以看出，信源角度的估計(jì)精度隨信噪比的增加而提高。由于本文算法采用稀疏對(duì)稱陣列，其陣列孔徑比均勻線陣[11]、互質(zhì)對(duì)稱陣列[12]以及特殊的幾何陣列結(jié)構(gòu)[13]大，因此本文算法角度估計(jì)精度要高于對(duì)比的算法。從圖5和圖6可以看出，4種算法在快拍數(shù)較小時(shí)參數(shù)估計(jì)精度較差，均隨著快拍數(shù)的增加而提高。從圖4和圖6可以看出，本文算法中的近場(chǎng)信源距離估計(jì)精度在隨著信噪比和快拍數(shù)的變化中都高于對(duì)比算法。本文算法和對(duì)比算法都是采用MUSIC算法來(lái)估計(jì)近場(chǎng)信源的距離，由于距離估計(jì)是基于角度估計(jì)，所以所提算法也具有更高的距離估計(jì)精度。

圖3 信源角度均方根誤差隨信噪比的變化

圖4 信源距離均方根誤差隨信噪比的變化

實(shí)驗(yàn)2：考慮信源均為近場(chǎng)時(shí)的情況。入射信號(hào)為2個(gè)4QAM近場(chǎng)窄帶信號(hào)，其位置分別為(θ1=-13°,r1=3.5λ)和(θ1=21°,r1=4.5λ)。

圖5 信源角度均方根誤差隨快拍數(shù)的變化

圖6 信源距離均方根誤差隨快拍數(shù)的變化

在本實(shí)驗(yàn)中，驗(yàn)證本文算法DOA估計(jì)精度隨信噪比和快拍數(shù)變化的情況。第一，信噪比從-5 dB到10 dB不斷變化，步長(zhǎng)為2 dB，快拍數(shù)為700，蒙特卡洛次數(shù)為500；第二，信噪比為7 dB，快拍數(shù)從100到1 100不斷變化，步長(zhǎng)為200。圖7為DOA估計(jì)的角度均方根誤差隨信噪比的變化情況，圖8為近場(chǎng)信源距離的均方根誤差隨信噪比的變化。圖9和圖10分別是方位角和距離隨快拍數(shù)的變化曲線圖。從圖7可以看出，本文算法的信源角度估計(jì)精度高于對(duì)比算法。由于本文算法所采用稀疏對(duì)稱陣列，其陣列孔徑要略大于互質(zhì)對(duì)稱陣列[11]和特殊的幾何陣列結(jié)構(gòu)[13]，且遠(yuǎn)大于均勻線陣[12]，因此本文算法具有更高的角度估計(jì)精度。從圖9和圖10可以看出，4種算法在快拍數(shù)較小時(shí)估計(jì)精度較差，隨著快拍數(shù)的增大性能有所改善并趨于穩(wěn)定。從圖8和圖10可以看出，本文算法的近場(chǎng)信源距離估計(jì)精度在隨著信噪比和快拍數(shù)變化下都高于對(duì)比算法。由于距離估計(jì)是基于角度估計(jì)，所提算法角度估計(jì)性能最好，所以也具有更高的距離估計(jì)精度。

圖7 信源角度均方根誤差隨信噪比的變化

圖8 信源距離均方根誤差隨信噪比的變化

圖9 信源角度均方根誤差隨快拍數(shù)的變化

圖10 信源距離均方根誤差隨快拍數(shù)的變化

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

本文在稀疏對(duì)稱陣列模型的基礎(chǔ)上，提出了一種混合信源DOA估計(jì)算法。通過(guò)不同子陣列接收數(shù)據(jù)的四階累積量運(yùn)算，構(gòu)造一個(gè)僅與信源角度有關(guān)的特殊四階累積量矩陣，然后使用MUSIC算法得到所有信源的角度估計(jì)；在每個(gè)角度估計(jì)的基礎(chǔ)上進(jìn)行距離維的搜索，進(jìn)而估計(jì)出近場(chǎng)信源的距離參數(shù)。本文算法避免了二維搜索和參數(shù)配對(duì)，稀疏對(duì)稱陣列擴(kuò)展了陣列孔徑，仿真實(shí)驗(yàn)結(jié)果表明，在相同的陣元情況下，本文算法具有更高的估計(jì)精度。