李小琳,張文俊

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

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

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

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

1 算法理論依據

1.1 Nested 陣列原理

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

圖1 二階Nested 線陣列結構示意圖Fig.1 Structural diagram of 2-level Nested linear array

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

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

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

多階Nested 陣列可以認為是二階Nested 線陣列的擴展。對M階Nested 陣列，每階陣列的天線數(shù)分別為L1,L2,··· ,LM ∈N+，天線位置分布為[2-3]

式中：

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

1.2 分數(shù)階傅里葉變換

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

圖3 分數(shù)階傅里葉變換示意圖Fig.3 Diagram of fractional Fourier transform

信號x(t)的FrFT 定義為

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

如圖3 所示，通過旋轉坐標軸，會出現(xiàn)兩個特殊的分數(shù)階傅里葉域。這兩個分數(shù)階傅里葉域分別定義為LFM 信號的解線調域和能量聚集域[12,16-17]，相應的FrFT 變換角度分別為

根據式FrFT 變換公式，寬帶LFM 信號的FrFT 為

將LFM 信號的解線調域FrFT 變換角度代入上式并化簡，得

式中,是一個常值。

由式(11)可知，一個寬帶線性調頻信號解線調域分數(shù)階傅里葉域變換后，該信號被解調為一個單頻信號，頻率為

由FrFT 的旋轉可加特性，可得

即LFM 信號的解線調域和能量聚集域之間存在一個傅里葉變換的關系[18]。

2 寬帶LFM 信號模型

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

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

每個子陣列包含L個天線，它們被放置在一個二階線性嵌套陣列結構中。線性Nested 子陣列由兩個級聯(lián)的均勻線陣(uniform linear array, ULA)組成(見圖5)，其中第一個ULA 包含間隔為dy的L1個天線，第二個ULA 包含間距為Dy= (L1+1)dy的L2=L-L2個天線。假設p為包含天線y軸坐標的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表示第?個天線與參考天線之間的延時，c為電磁波傳播速度，ξk=sinθksinφk表示第k個信號的y軸方向余弦；n1,?(t)為由第一個子陣列的第?個天線測量到的加性高斯白噪聲。

第二個子陣列中第?個天線的輸出信號表示為

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

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

對寬帶LFM 信號sk(t)進行FrFT，可得

根據FrFT 時移特性，信號sk(t)延時τ后進行FrFT，則有

對第一個子陣列中第?個天線的輸出信號進行FrFT，根據FrFT 的線性特性可得

因為FrFT 不改變噪聲的時域白特性，也不改變噪聲能量，在任何FRF 域，高斯白噪聲均不呈現(xiàn)能量聚集，所以零均值高斯白噪聲經過FrFT 仍為零均值高斯白噪聲[19-20]。

綜上所述，在FRF 域，由第一個子陣列測量的寬帶LFM 信號可以表示為L×1 維向量：

式中：S(u,α) = (S1(u,α),S2(u,α),··· ,SK(u,α))T是K ×1維入射信號向量；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個信號的L×1 維子陣列導向矢量，其表達形式為

類似地，由第二個子陣列測量的L×1 維寬帶LFM 信號向量具有形式，

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

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

3 算法推導

3.1 DOA 矩陣公式

所提算法的第一步是構造具有增強自由度的DOA 矩陣?；谏鲜黾僭O，第一個Nested子陣測量信號的自相關矩陣為

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

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

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

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

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

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

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

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

Φξ是一個K×K維對角矩陣，由于第k個對角線項為

因此，可以將V a表示為

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

類似地，由第一和第二Nested 子陣所測信號的互相關矩陣為

Rb的向量形式為

移除重復項，可得

由于向量與具有相同的大小，可以將向量分解為L2(L1+1)個重疊子向量以構造L2(L1+1)×L2(L1+1)矩陣

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

基于B的范德蒙結構和前文假設，使用B是全列秩的特性，則有

式中：上標?表示偽逆運算符。從式(35)和(36)，可構造如下DOA 矩陣：

3.2 二維空間角度估計

D的特征分解可寫為

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

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

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

4 仿真結果

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

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

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

圖7 空間角實際值與估計結果Fig.7 True values and estimated values of space angles

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

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

5 結束語

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