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

0 引言

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

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

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

1 信號模型

本文所采用的陣列模型如圖1所示，陣元總數(shù)為2(N1+N2)+1,由3個子陣列組成。其中子陣1陣元數(shù)和陣元間距分別為2N1+1和d，子陣2和子陣3的陣元數(shù)和陣元間距分別為N2和(2N1+1)d,子陣列1與子陣列2和子陣列3之間的距離分別為(N1+1)d。圖中各陣元坐標為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個陣元到中心陣元的距離。當陣元數(shù)相同時，均勻線陣的陣列孔徑為(2N1+2N2+1)d，本文所采用陣列的陣列孔徑為(2N2+4N1N2)d，可以看出，本文的陣列具有更大的陣列孔徑。

圖1 稀疏對稱陣列模型

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

t=1,…,T

(1)

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

(2)

(3)

式中，λ為信號的波長，rk，θk分別為第k個近場信源的距離和角度參數(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個信源的信號矢量，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維方向矢量，當信源位于近場時，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)

當信源位于遠場時，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)

本文中，對信號作如下假設：

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

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

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

2 基于稀疏對稱陣列的信源定位

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

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

(8)

(9)

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

c1(m+N1+1)=

(10)

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

c2(m+N1+1)=

(11)

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

c3(m+N1+1)=

(12)

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

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

(13)

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

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

(14)

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

(15)

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

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

重構(gòu)方法：

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

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

向量c構(gòu)造一個(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)造了一個僅包含信源角度信息的特殊矩陣C，類似于遠場協(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)

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

2.2 信源DOA估計

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

(21)

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

(22)

2.3 近場信源距離估計

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

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

(23)

(24)

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

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

(25)

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

(26)

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

(27)

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

2.4 計算復雜度分析

3 計算機仿真

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

(28)

(29)

實驗1：信號源為近場遠場混合信源。入射信號為2個4QAM窄帶信號，近場信源與遠場信源位置分別為(θ1=-13°,r1=3.5λ)和(θ2=21°,r2=∞)。

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

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

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

實驗2：考慮信源均為近場時的情況。入射信號為2個4QAM近場窄帶信號，其位置分別為(θ1=-13°,r1=3.5λ)和(θ1=21°,r1=4.5λ)。

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

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

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

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

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

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

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

4 結(jié)束語

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