李小琳,張文俊

(1. 上海大學(xué)通信與信息工程學(xué)院, 上海 200444; 2. 上海機(jī)電工程研究所, 上海 201109)

為了確保波達(dá)方向(direction of arrival, DOA)估計(jì)的精度，陣列相鄰天線間距一般要求不大于入射信號(hào)的半波長(zhǎng)，而當(dāng)工作頻率較高時(shí)，相應(yīng)波長(zhǎng)變短，這就要求陣列天線間距變小, 進(jìn)而導(dǎo)致天線間的互耦效應(yīng)。為了實(shí)現(xiàn)高分辨率，陣列必須具有較大的天線孔徑, 這就需要大量的天線來保證。為了節(jié)約成本并確保陣列性能，稀疏陣列應(yīng)運(yùn)而生。常見的稀疏陣列主要包括最小冗余陣列[1]、Nested(嵌套)陣列[2-4]、互質(zhì)陣列[5]和分布式陣列[6]等。最小冗余陣列是最早出現(xiàn)的一種稀疏陣列形式，但結(jié)構(gòu)設(shè)計(jì)復(fù)雜。2010 年，Pal 等[2]在最小冗余陣列研究的基礎(chǔ)上，提出了一種新穎的稀疏陣列結(jié)構(gòu)，即Nested 陣列。

隨著雷達(dá)與通訊技術(shù)的發(fā)展，調(diào)制類寬帶(linear frequency modulation, LFM)信號(hào)以其承載信息量大、抗干擾能力強(qiáng)等優(yōu)勢(shì)，得到了愈來愈廣泛的應(yīng)用。但由于寬帶信號(hào)的陣列流型與頻率有關(guān)，原有的DOA 估計(jì)經(jīng)典算法已經(jīng)無法直接應(yīng)用于Nested 陣列進(jìn)行寬帶信號(hào)測(cè)角。針對(duì)Nested 陣列高精度DOA 測(cè)角，Pal 等[2]基于二階Nested 陣列提出了類似空間平滑(spatial smoothing) MUSIC 算法，但該算法僅可有效估計(jì)欠定窄帶平穩(wěn)信號(hào)的DOA。張曉鳳等[7]和Huang 等[8]提出了基于矩陣重構(gòu)的ESPRIT 算法，該算法雖不需要求出信號(hào)子空間，但仍需進(jìn)行參數(shù)配對(duì)。在空間平滑算法的基礎(chǔ)上,Coventry 等[9]結(jié)合多項(xiàng)式矩陣，實(shí)現(xiàn)了互質(zhì)或Nested 等稀疏陣列的寬帶信號(hào)空間角度估計(jì)，但在空間平滑這一步驟上，該算法性能有待進(jìn)一步提高。Shi 等[10]通過構(gòu)造參數(shù)化字典，利用聯(lián)合稀疏結(jié)構(gòu)，提出了一種聯(lián)合估計(jì)稀疏信號(hào)和優(yōu)化參數(shù)化字典的迭代最小化方法，解決了Nested 陣列對(duì)寬帶信號(hào)進(jìn)行欠定到達(dá)方向估計(jì)的問題，該方法雖然有效，但未提及是否可以實(shí)現(xiàn)二維測(cè)角。

針對(duì)寬帶信號(hào)測(cè)角難題，一種有效的方法是通過聚焦將各頻率點(diǎn)上的觀測(cè)數(shù)據(jù)在某一子域上對(duì)齊，得到聚焦合成的觀測(cè)量，依此進(jìn)行寬帶LFM 信號(hào)的DOA 估計(jì)[11]。分?jǐn)?shù)階傅里葉變換是一種有效的時(shí)頻域分析工具[12]，而時(shí)頻分析是處理非平穩(wěn)信號(hào)的一種重要方法。作為一種典型的非平穩(wěn)信號(hào)，LFM 信號(hào)在分?jǐn)?shù)階傅里葉域上分別呈現(xiàn)能量聚集性和解線調(diào)等特性，這就使得利用FrFT 旋轉(zhuǎn)時(shí)頻域的特性對(duì)寬帶LFM 信號(hào)進(jìn)行DOA 估計(jì)有了理論依據(jù)。為此，本工作基于分?jǐn)?shù)階傅里葉變換，結(jié)合DOA 矩陣方法[13-14]，提出一種Nested 陣列寬帶LFM 信號(hào)二維空間測(cè)角算法，實(shí)現(xiàn)了對(duì)非平穩(wěn)寬帶LFM 信號(hào)的空間角度精確估計(jì)。

1 算法理論依據(jù)

1.1 Nested 陣列原理

Nested 陣列是一種天線布置在規(guī)則柵格上的稀疏陣列，一般由兩個(gè)或多個(gè)陣元間距不同的均勻線陣嵌套而成。相對(duì)于相同天線數(shù)的傳統(tǒng)均勻陣列，Nested 陣列具有更大的陣列孔徑和更多的自由度，因而具有更高的測(cè)角精度與分辨率。下面以常見的二階Nested 線陣列為例，說明Nested 陣列的結(jié)構(gòu)和特性。二階Nested 線陣列通常是由兩個(gè)相鄰的均勻線陣嵌套組成：第一級(jí)(內(nèi)層)均勻線陣有L1個(gè)天線, 天線間距為d；第二級(jí)(外層)均勻線陣有L2=L-L1個(gè)天線,其天線間距為D=(L1+1)d。圖1 為天線總數(shù)為L(zhǎng)的二階Nested 線陣列的結(jié)構(gòu)示意圖。

圖1 二階Nested 線陣列結(jié)構(gòu)示意圖Fig.1 Structural diagram of 2-level Nested linear array

對(duì)于一個(gè)有L=L1+L2個(gè)天線的二階Nested 線陣列，其差分協(xié)同陣列可以生成2L2(L1+1)-1 個(gè)虛擬天線[2]，此虛擬陣為均勻線陣，天線位置為

式中，K=L2(L1+1)-1。虛擬陣列的陣元分布如圖2 所示。一個(gè)有L=L1+L2個(gè)天線的二階Nested 線陣列可以獲得2L2(L1+1)-1個(gè)自由度。

圖2 虛擬陣列的陣元分布Fig.2 Elements distribution of virtual array

多階Nested 陣列可以認(rèn)為是二階Nested 線陣列的擴(kuò)展。對(duì)M階Nested 陣列，每階陣列的天線數(shù)分別為L(zhǎng)1,L2,··· ,LM ∈N+，天線位置分布為[2-3]

式中：

M階Nested 陣列的總物理天線個(gè)數(shù)為，對(duì)應(yīng)的自由度(degree of freedom,DoF)為

1.2 分?jǐn)?shù)階傅里葉變換

作為普通Fourier 變換的擴(kuò)展，分?jǐn)?shù)階傅里葉變換(fractional Fourier transform,FrFT)非常適于處理瞬時(shí)頻率特性呈線性變化的LFM 類信號(hào). 圖3 為FrFT 示意圖?？梢钥闯觯喝绻盘?hào)的傅里葉變換可看成將其由時(shí)間t軸上逆時(shí)針旋轉(zhuǎn)π/2 后到頻率f軸上的表示，則FrFT可以看成是將信號(hào)在時(shí)間t軸上逆時(shí)針旋轉(zhuǎn)α角度到u軸上的射線表示[15]。

圖3 分?jǐn)?shù)階傅里葉變換示意圖Fig.3 Diagram of fractional Fourier transform

信號(hào)x(t)的FrFT 定義為

式中：p為FrFT 的階數(shù)，周期為4；u=tanα=tan(pπ)；Kp(t,u)為FrFT 的變換核。根據(jù)定義可知，F(xiàn)rFT 的逆變換為角度α=-pπ/2 的FrFT，即

如圖3 所示，通過旋轉(zhuǎn)坐標(biāo)軸，會(huì)出現(xiàn)兩個(gè)特殊的分?jǐn)?shù)階傅里葉域。這兩個(gè)分?jǐn)?shù)階傅里葉域分別定義為L(zhǎng)FM 信號(hào)的解線調(diào)域和能量聚集域[12,16-17]，相應(yīng)的FrFT 變換角度分別為

根據(jù)式FrFT 變換公式，寬帶LFM 信號(hào)的FrFT 為

將LFM 信號(hào)的解線調(diào)域FrFT 變換角度代入上式并化簡(jiǎn)，得

式中,是一個(gè)常值。

由式(11)可知，一個(gè)寬帶線性調(diào)頻信號(hào)解線調(diào)域分?jǐn)?shù)階傅里葉域變換后，該信號(hào)被解調(diào)為一個(gè)單頻信號(hào)，頻率為

由FrFT 的旋轉(zhuǎn)可加特性，可得

即LFM 信號(hào)的解線調(diào)域和能量聚集域之間存在一個(gè)傅里葉變換的關(guān)系[18]。

2 寬帶LFM 信號(hào)模型

假定有K個(gè)寬帶LFM 信號(hào)入射至非均勻平行陣列，信源與平行Nested 陣列的幾何關(guān)系如圖4 所示。非均勻平行陣列由兩個(gè)子陣列組成，每個(gè)子陣列是一個(gè)二階Nested 線陣列。第一個(gè)子陣列沿y軸展開，第二個(gè)子陣列沿x軸平移距離dx展開。第k個(gè)寬帶LFM 信號(hào)sk(t)的方位角和俯仰角分別為φk和θk，且φk ∈[0,]，θk=[0,/2]。

圖4 信源與平行Nested 陣列幾何示意圖Fig.4 Geometric sketch of sources and parallel Nested array

每個(gè)子陣列包含L個(gè)天線，它們被放置在一個(gè)二階線性嵌套陣列結(jié)構(gòu)中。線性Nested 子陣列由兩個(gè)級(jí)聯(lián)的均勻線陣(uniform linear array, ULA)組成(見圖5)，其中第一個(gè)ULA 包含間隔為dy的L1個(gè)天線，第二個(gè)ULA 包含間距為Dy= (L1+1)dy的L2=L-L2個(gè)天線。假設(shè)p為包含天線y軸坐標(biāo)的L×1維向量，則有p=(0,1,L,L1-1,L1,2L1,L,L2L1)dy.

圖5 Nested 子陣列天線陣元分布圖Fig.5 Elements distribution map of Nested sub-array

式中：τ1,?k=p?ξk/c表示第?個(gè)天線與參考天線之間的延時(shí)，c為電磁波傳播速度，ξk=sinθksinφk表示第k個(gè)信號(hào)的y軸方向余弦；n1,?(t)為由第一個(gè)子陣列的第?個(gè)天線測(cè)量到的加性高斯白噪聲。

第二個(gè)子陣列中第?個(gè)天線的輸出信號(hào)表示為

式中：τ2,?k= (p?ξk+dxγk)/c表示第二個(gè)子陣列上第?個(gè)天線與參考天線之間的延時(shí)，γk=sinθkcosφk是第k個(gè)信號(hào)的x軸方向余弦；n2,?(t)為由第二個(gè)子陣列的第?個(gè)天線測(cè)量到的加性高斯白噪聲。由式(14)和(15)可知，兩個(gè)子陣列的流形矩陣都是時(shí)變的，傳統(tǒng)DOA 估計(jì)算法不能直接用于寬帶LFM 信號(hào)。

本工作利用分?jǐn)?shù)階傅里葉變換工具，首先把時(shí)頻域內(nèi)的寬帶問題轉(zhuǎn)換成為頻率響應(yīng)函數(shù)(frequency response function, FRF)域內(nèi)的窄帶問題，進(jìn)而提出新的算法實(shí)現(xiàn)Nested 平行陣列的寬帶LFM 信號(hào)二維測(cè)角。

對(duì)寬帶LFM 信號(hào)sk(t)進(jìn)行FrFT，可得

根據(jù)FrFT 時(shí)移特性，信號(hào)sk(t)延時(shí)τ后進(jìn)行FrFT，則有

對(duì)第一個(gè)子陣列中第?個(gè)天線的輸出信號(hào)進(jìn)行FrFT，根據(jù)FrFT 的線性特性可得

因?yàn)镕rFT 不改變?cè)肼暤臅r(shí)域白特性，也不改變?cè)肼暷芰浚谌魏蜦RF 域，高斯白噪聲均不呈現(xiàn)能量聚集，所以零均值高斯白噪聲經(jīng)過FrFT 仍為零均值高斯白噪聲[19-20]。

綜上所述，在FRF 域，由第一個(gè)子陣列測(cè)量的寬帶LFM 信號(hào)可以表示為L(zhǎng)×1 維向量：

式中：S(u,α) = (S1(u,α),S2(u,α),··· ,SK(u,α))T是K ×1維入射信號(hào)向量；N1(u,α) =(N1,1(u,α),N1,2(u,α),··· ,N1,L(u,α))T是L×1 維噪聲向量；A(u,α)=(a1(u,α),a2(u,α),··· ,aK(u,α))是L×K維子陣列流形矩陣，ak(u,α)是第k個(gè)信號(hào)的L×1 維子陣列導(dǎo)向矢量，其表達(dá)形式為

類似地，由第二個(gè)子陣列測(cè)量的L×1 維寬帶LFM 信號(hào)向量具有形式，

式中：Φ(u,α)是K × K維對(duì)角矩陣，它的第k個(gè)對(duì)角線元素為；N2(u,α) = [N2,1(u,α),N2,2(u,α),··· ,N2,L(u,α)]是第二個(gè)子陣的L×1 維加性噪聲向量。

本工作討論的問題是對(duì)不同時(shí)刻tn(n= 1,2,··· ,N)寬帶LFM 入射信號(hào)sk(t)的方向(θk,fk) (k= 1,2,··· ,K)進(jìn)行估計(jì)并提出一種Nested 陣列測(cè)量寬帶信號(hào)二維空間角的新算法，為此進(jìn)行以下假設(shè)：①并行Nested 陣列天線特性一致，為全相寬帶接收天線；②入射信號(hào)為點(diǎn)源信號(hào)，各點(diǎn)源相對(duì)于Nested 陣列方向唯一確定；③入射信號(hào)的數(shù)量是預(yù)先知道或正確估計(jì)的；④入射信號(hào)的功率和帶寬相等, 且彼此不相干；⑤噪聲是均值為0、方差為σ2寬帶高斯白噪聲，與入射信號(hào)不相關(guān)。

3 算法推導(dǎo)

3.1 DOA 矩陣公式

所提算法的第一步是構(gòu)造具有增強(qiáng)自由度的DOA 矩陣。基于上述假設(shè)，第一個(gè)Nested子陣測(cè)量信號(hào)的自相關(guān)矩陣為

式中：上標(biāo)H 表示共軛轉(zhuǎn)置；I表示單位矩陣；對(duì)角矩陣Rs= diag(ρ1,ρ2,··· ,ρK)是K×K階入射信號(hào)相關(guān)矩陣，其中ρk是第k個(gè)信源的“聚集”功率。由于Rs是對(duì)角矩陣，Ra的向量形式可簡(jiǎn)化為

式中：上標(biāo)*代表復(fù)共軛；⊙代表Khatri-Rao乘積；ρ=(ρ1,ρ2,··· ,ρK)T；E=(eT1,eT2,··· ,eTL)T，除了第?位置的1 外，e?是全零的向量，上標(biāo)T 是轉(zhuǎn)置算子[21]。由于用特征值估計(jì)法可以很容易地估計(jì)出式(23)中的噪聲功率，因此可以去掉式(23)中的噪聲項(xiàng)，得到

接下來，利用Nested 陣列的差分共陣屬性，很容易驗(yàn)證(A*⊙A)的每列只有=2L2(L1+1)-1個(gè)不同的項(xiàng)。首先，從(A*⊙A)中刪除重復(fù)行，并對(duì)其余行進(jìn)行排序，使第?行與y軸坐標(biāo)相關(guān)聯(lián)；然后，可以通過選擇va的非重復(fù)項(xiàng)來形成一個(gè)新的向量va

式中：B(u,α)=(b1(u,α),b2(u,α),...,bK(u,α))是一個(gè)階矩陣，第k列是一個(gè)維向量，

很明顯，式(25)用流形矩陣B(u,α)和信號(hào)向量ρ表示單快拍測(cè)向問題。為了估計(jì)信號(hào)的方向，向量被劃分為L(zhǎng)2(L1+1)個(gè)重疊子向量，每個(gè)子向量包含L2(L1+1)個(gè)元素。那么，第m個(gè)子向量可以表示為

式中：Bm(u,α)是矩陣B(u,α)的第m到(m+L2(L1+1)-1)行。

對(duì)于所有m=1,2,··· ,L2(L1+1)，可以形成一個(gè)L2(L1+1)×L2(L1+1)維矩陣

由于Bm(u,α)與B1(u,α)相關(guān)，即

Φξ是一個(gè)K×K維對(duì)角矩陣，由于第k個(gè)對(duì)角線項(xiàng)為

因此，可以將V a表示為

式中：D(·)表示通過將矢量放在主對(duì)角線上的括號(hào)中來生成對(duì)角矩陣的運(yùn)算符。

類似地，由第一和第二Nested 子陣所測(cè)信號(hào)的互相關(guān)矩陣為

Rb的向量形式為

移除重復(fù)項(xiàng)，可得

由于向量與具有相同的大小，可以將向量分解為L(zhǎng)2(L1+1)個(gè)重疊子向量以構(gòu)造L2(L1+1)×L2(L1+1)矩陣

是的第m到(m+L2(L1+1)-1)行。矩陣可以表示為

基于B的范德蒙結(jié)構(gòu)和前文假設(shè)，使用B是全列秩的特性，則有

式中：上標(biāo)?表示偽逆運(yùn)算符。從式(35)和(36)，可構(gòu)造如下DOA 矩陣：

3.2 二維空間角度估計(jì)

D的特征分解可寫為

由式(37)和(38)可以很容易地證明，Φγ的第k個(gè)對(duì)角元素等于Ψγ的第k個(gè)非零對(duì)角元素。令[Ψγ]k,k為D的非零特征值，并且設(shè)fk(k= 1,2,··· ,K)是與這些特征值相關(guān)的特征向量，則x軸方向余弦估計(jì)可以從fk估計(jì)為

式中：fk,m表示fk的第m個(gè)元素。y軸方向余弦估計(jì)可通過下式計(jì)算得到：

本工作提出的算法從矩陣的特征值估計(jì)x軸方向余弦，并且從D的特征向量估計(jì)y軸方向余弦ξk。由于特征值及其特征向量是成對(duì)的，因此估計(jì)的方向余弦和將自動(dòng)匹配而無需任何額外的配對(duì)匹配處理?？梢宰R(shí)別的信源的最大數(shù)量由D的信號(hào)子空間維度確定。由于D的維數(shù)為L(zhǎng)2(L1+1)×L2(L1+1)，因此該算法可以利用2L物理傳感器陣列解析K≤L2(L1+1)-1個(gè)寬帶LFM 信號(hào)。

4 仿真結(jié)果

為了驗(yàn)證所提算法的分辨性能，使用L1= 3,L2= 2的并行嵌套陣列測(cè)量7 個(gè)寬帶LFM入射信號(hào)。設(shè)LFM 信號(hào)的初始頻率分別為0, 50, 100, 150, 200, 230 和280 MHz, 帶寬同為300 MHz。本算法首先需要對(duì)寬帶LFM 信號(hào)進(jìn)行FrFT。圖6 為疊加信噪比(signal noise ratio,SNR)為20 dB 白噪聲后，初始頻率為100 MHz 寬帶LFM 信號(hào)FrFT 前后的頻譜。

圖6 初始頻率100 MHz 的LFM 信號(hào)FrFT 前后頻譜Fig.6 Frequency spectrum of LFM signal before and after FrFT with initial frequency of 100 MHz

若7 個(gè)信源的方位、高低角度分別為(30°,10°)，(30°,50°)，(50°,60°)，(80°,20°)，(100°,80°)，(110°,50°)和(160°,30°)，信噪比和快拍數(shù)分別設(shè)置為20 dB 和1 000，共進(jìn)行了500 次獨(dú)立蒙特卡羅試驗(yàn)。圖7 為10 次獨(dú)立運(yùn)行的空間角估計(jì)。從圖中可以看出，所提算法可以精確地分辨出這7 個(gè)寬帶LFM 信號(hào)。

圖7 空間角實(shí)際值與估計(jì)結(jié)果Fig.7 True values and estimated values of space angles

應(yīng)用本工作所提算法，將L1=L2=3 的并行嵌套陣列與L=8，L=12 的平行ULA 陣列進(jìn)行性能對(duì)比。對(duì)于角度為(10°,50°)的信號(hào)入射，圖8 為不同陣列結(jié)構(gòu)算法所估計(jì)角度的均方根誤差隨SNR 的變化關(guān)系。從圖中可以看出，所提算法的性能接近L= 12的平行ULA陣列，而明顯優(yōu)于L=8的平行ULA 陣列。

圖8 隨信噪比變化的均方根誤差對(duì)比曲線Fig.8 Root mean square error (RMSE) contrast curves as the function of SNR

5 結(jié)束語

本工作提出了一種使用平行Nested 陣列來估計(jì)寬帶LFM 信號(hào)二維方向的新算法。該算法的關(guān)鍵思想是首先基于分?jǐn)?shù)階傅里葉變換，把寬帶LFM 信號(hào)轉(zhuǎn)換為FRF 域的窄帶信號(hào)，然后利用嵌入在并Nested 陣列的自相關(guān)矩陣和互相關(guān)矩陣中的差分共陣列屬性，形成具有增強(qiáng)自由度的DOA 矩陣，在不涉及二維非線性搜索和參數(shù)匹配的情況下，可實(shí)現(xiàn)自動(dòng)配對(duì)的二維角度估計(jì)。