劉魯濤，司錫才，王立國

(哈爾濱工程大學 信息與通信工程學院，黑龍江 哈爾濱 150001)

在高技術條件下的現(xiàn)代戰(zhàn)爭中，被動制導系統(tǒng)必須能在各種干擾背景下對雷達目標進行正確的捕捉與跟蹤.而現(xiàn)代雷達為提高自身工作性能和抗干擾能力和躲避被動偵測系統(tǒng)的跟蹤，在電子戰(zhàn)雙方奪取制電磁權，采取的主要技術措施包括擴展雷達工作頻域范圍、頻率捷變、頻率分集、超寬帶等等，并配備雷達誘餌.針對雷達技術的變化，跟蹤多個新體制雷達目標，使得傳統(tǒng)的比幅/比相等被動偵測技術難以應付，因此對被動制導系統(tǒng)提出更高的要求.基于陣列信號處理的空間譜估計技術更多的應用在被動制導系統(tǒng)當中，用以跟蹤多部雷達信號［1］.但傳統(tǒng)的譜估計方法，如MUSIC［2］等等，要求多目標信號工作頻率一致，并通過角度搜索估計目標方向.而對于同時到達目標信號工作頻率不同的情況，傳統(tǒng)譜估計算法的角度搜索不再適用.只能采用波達方向-頻率聯(lián)合估計的參數(shù)估計方法.如文獻［3-4］是基于狀態(tài)空間實現(xiàn)并結合ESPRIT方法，將波達方向-頻率二維參數(shù)估計問題轉化為2個一維參數(shù)估計.文獻［5］基于PM(propagator method)對相關矩陣進行降維處理，給出波達方向-頻率估計進而得出偵測目標方向.上述算法的局限性在于算法都是以ESPRIT為基礎，要求特定的陣列結構且相對譜峰搜索技術精度較低;而且對于處理脈沖間捷變頻的信號來講，相鄰脈沖頻率跳變使信號子空間不停變化，需要在各個信號脈沖內重復的對信號波達方向－頻率進行估計，工作量巨大，工程難以實現(xiàn).文章提出了一種跟蹤多部捷變頻雷達目標的算法，其基本思想將多個捷變頻雷達信號，映射成為單一參考頻率的信號，從而構造一個對應單一頻率的組合信號子空間，采用空間譜峰搜索的方法，得到多部捷變頻雷達信號的方向估計.

1 信號模型

假設D個遠場窄帶信號入射到由M個陣元構成的天線陣(M＞D).在信號源不相關情況下，設陣列接收數(shù)據(jù)為x(t)，第i個陣元的接收數(shù)據(jù)為

式中:ni(t)為第i個天線接收噪聲;假定各陣元的噪聲是均值為零的平穩(wěn)空間白噪聲過程，方差為σ2，并且噪聲之間不相關，且與信號不相關;sd表示第d個目標信號;θd表示第d個信號源的到達角(DOA);τi(θd)相對于參考陣元，第d個信源傳播到第i個陣元的延遲時間.在均勻線陣情況下，相鄰陣元間的距離為l，τi(θd)=(i－1)lsin θd/c(c為波的傳播速度).

如果信號源是窄帶信號，可以將信號時間延遲近似等效為信號相位延遲，則信號可以寫成:

M個陣元在特定時刻接收的信號可以寫成:

式中:x(t)=［x1(t)x2(t)… xM(t)］T為M維的接收數(shù)據(jù)矢量，s(t)=［s1(t) s2(t) …sD(t)］T為D維信號矢量，A=A(Ω，F(xiàn))=［a(θ1，f1) a(θ2，f2)…a(θD，fD)］為M×D維的陣列流型矩陣，a(θd，fd)=［1 e－j2πfdτ1(θd)… e－j2πfdτM－1(θd)］T為M維的導向矢量，Ω=［θ1θ2…θD］目標信號到達角集合，F(xiàn)=［f1f2… fD］為目標信號頻率集合，n(t)=［n1(t)n2(t) … nM(t)］T為M維的噪聲矢量.

捷變頻雷達的瞬時工作頻率在頻率捷變帶寬中的多個頻率點隨機跳動，其信號的特點是相鄰2個發(fā)射脈沖的載頻在一個很大的范圍內以很高的速度躍變，也就是說在相鄰脈沖內雷達信號的載頻是不相等的.如圖1(a)和(b)為第d部雷達信號脈沖間頻率跳變時域和頻譜示意圖.捷變頻雷達原理性能見文獻［6，10］.從式(2)、(3)可以看出由于不同脈沖間每部雷達信號頻率fd的跳變，導致每部雷達信號對應導向矢量a(θd，fd)在不同的脈沖間也是不斷變化的.假設陣列觀測時間足夠長，采樣數(shù)據(jù)來自K個回波脈沖，每個脈沖包含J個采樣點(每個脈沖采樣點數(shù)).第k(k=1，2，…，K)個脈沖內的采樣數(shù)據(jù):

這里S=［s1s2… sD］T為D×J維信號矩陣，不同脈沖采樣數(shù)據(jù)Xk(M×J維)對應的不同的陣列流型Ak.所以簡單組合所有采樣數(shù)據(jù)進行測向不能實現(xiàn).而對于文獻［3-5］等波達方向—頻率參數(shù)聯(lián)合估計的方法，需要在各個脈沖間重復使用，工作量巨大;而且由于每個脈沖內采樣數(shù)目J有限，難以達到滿意的測向精度.

圖1 捷變頻雷達信號時頻域示意圖Fig.1 Scheme of radar signal of frequency agility

2 參考頻率子空間估計及算法描述

2.1 參考頻率子空間估計

引理 在一定條件下M×D維矩陣A(Ω，F(xiàn))有rank［A(Ω，F(xiàn))］=D，則一定存在一個非奇異的M× M維變換矩陣T(F)，使得

式中:A(Ω，f0)為M×D矩陣.

證明 因為A(Ω，F(xiàn))和A(Ω，f0)都是秩為D的M×D維矩陣，則存在M×(M－D)維矩陣B(F)和B(f0)，使得M×M維矩陣［A(Ω，F(xiàn))|B(F))］和［A(Ω，f0)|B(f0)］是非奇異陣.

T(F)的一個明顯的選擇為

這里變換矩陣T(F)不是唯一的.通常情況下可以選擇B(f0)=BT(F)=［0D×(M－D)|IM－D］，這樣選擇的好處可以減少計算復雜度.通過分塊矩陣的逆變換公式，把M×M維矩陣的逆變換轉化為D× D維矩陣的逆變換的問題.可以得到變換矩陣:

式中:

其中，Al和Au分別為矩陣A的上側D×D維塊陣和下側的(M－D)×D維塊陣.

使用變換矩陣T可以將不同頻率觀測矩陣變化成同一頻率的觀測矩陣，第k個采樣脈沖內的觀測數(shù)據(jù)矩陣變換為

則由式(4)和(7)可以將采樣協(xié)方差矩陣變?yōu)?/p>

式中:

以上變換將導向矢量陣A(Ω，F(xiàn))變?yōu)锳(Ω，f0)，重新構造了信號子空間，使得不同的脈沖內的采樣數(shù)據(jù)都對應相同的頻率f0，通過這個方法使Rk(k=1，2，…，K)中的跳頻信號能量都聚集參考頻率f0內，并對應的信號子空間A(Ω，f0)內，因此可以對所有脈沖內的采樣數(shù)據(jù)組合在一起，進行組合子空間測向處理.從而改善了對多個捷變頻信號方向的估計性能.但是也產生了一個問題，這種變換可能使白噪聲產生變化，改變輸出信號處理的SNR.文獻［7］分析了矩陣變換對白噪聲矩陣的影響.

2.2 變換矩陣構造與算法描述

由式(5)可以看出，對于每個脈沖內數(shù)據(jù)的變換矩陣構造都需要求得該信號的頻率和到達角度的預估計.對于雷達信號來講，由于脈沖的間隔時間相對較短，目標信號的到達角在各個脈沖基本不變，則角度的預估計在任意一個快拍內完成即可(預估計的目標到達角度可以存在較大誤差).而對與信號的頻率是在不停的跳變的，所以對信號頻率需要對每個快拍分別估計.

由此可以對第1個脈沖內的采樣數(shù)據(jù)采用角度－頻率聯(lián)合估計［4，8］，得到信號的頻率和角度預估計，構造T1.對于接下來的脈沖內數(shù)據(jù)只估計信號頻率即可.通過觀測數(shù)據(jù)的特征矢量對預估計的角度和頻率進行參數(shù)配對.例如第k個脈沖的數(shù)據(jù)包括兩部雷達信號，信號頻率與角度的參數(shù)為f1、f2、θ1、θ2.則可能的信號陣列流行有以下2種可能:

先求得第k個脈沖數(shù)據(jù)Xk的協(xié)方差陣，對其進行特征分解，得到信號和噪聲的特征子空間Us和UN.由于信號流型與噪聲的特征子空間UN正交，比較‖AH1UNUHNA1‖和‖AH2UNUHNA2‖大小，即可完成頻率與角度參數(shù)配對，從而得到A(Ω，F(xiàn)).選定參考頻率f0(參考頻率選擇在探測目標信號頻率頻段附近為宜)后，根據(jù)式(5)構造變換矩陣Tk.

算法實現(xiàn)步驟:

1)對第1個脈沖采樣數(shù)據(jù)進行頻率-角度預估計，使用文獻［4，8］方法，得到第1個脈沖內各個信號的波達方向-頻率預估計，作為初始值;

2)對余下不同脈沖的采樣數(shù)據(jù)分別進行FFT或采用ESPRIT方法［9］，得到各脈沖內同時到達信號的頻率估計;

3)頻率-角度參數(shù)配對;

4)選定參考頻率f0，構造各個脈沖內的數(shù)據(jù)變換矩陣Tk(F)，k=1，2，…，K;

5)通過式(7)形成變換后的觀測數(shù)據(jù)Yk，使其陣列流行對應于同一頻率f0，構造新的觀測數(shù)據(jù)組合Y=［Y1Y2…YK］;

6)構造Y的協(xié)方差陣R，然后估計雷達目標波達方向(如使用MUSIC);

步驟4)～6)迭代可以改善估計性能.

3 計算復雜度分析

算法主要涉及的計算包括頻率估計、協(xié)方差矩陣計算、特征值分解、變換矩陣構造以及譜峰搜索.步驟2)頻率預估計需要的復數(shù)乘法數(shù)O(KJlbJ)，步驟3)中參數(shù)配對及步驟6)共K次特征分解需要O (M2DK)+O(M3K)次復數(shù)乘法;步驟4)K個變換矩陣構造見式(5)(6)，需要復數(shù)乘法數(shù)為O(D3K) +O(D2MK);步驟6)構造Yk需要O(M2JK)個復數(shù)乘法;構造協(xié)方差矩陣O(KJM(M+1)/2)復數(shù)乘法數(shù);步驟7)采用MUSIC譜峰搜索得到目標到達角，假設采用線陣(一維搜索步長為1°，搜索范圍［－ 90°，+ 90°］)需要的乘法數(shù)為O(181M(M－D)2).

由于M＞D，KJ＞J＞M.計算量主要集中在譜峰搜索和步驟5)的Yk構造.算法需要復數(shù)乘法數(shù)約為O(181M(M－D)2)+O(M2N)(N=JK為觀測數(shù)據(jù)總的采樣數(shù)).在實時處理系統(tǒng)中由于采樣點數(shù)目有限，一般滿足181M＞N時，計算量主要集中在譜峰搜索步驟.所以本算法與一般MUSIC算法計算量基本相同.同時相比于參數(shù)估計方法來講，保證了譜峰搜索的高精度優(yōu)點.

4 實驗結果與分析

為了驗證所提出算法的有效性及測向精度，對該算法進行計算機仿真試驗.在下面仿真中，假設各向同性陣元數(shù)為M=8的的等距線陣.相鄰陣元間距為d=c/2f'，f'信號可達到的最大頻率.各個信號相互獨立，陣列輸入噪聲為平穩(wěn)零均值白高斯過程.采樣頻率符合奈奎斯特要求，觀測脈沖數(shù)為K=4，每個脈沖內采樣J=100.SNR定義為單個陣元信號能量與噪聲能量的比值.

4.1 對2個捷變頻信號的到達角估計

2信號來波方向分別為 θ=10°和 θ=30°，SNR=10 dB，2個信號載頻在［500 MHz，600 MHz］任意跳變，參考頻率選擇為f0=600 MHz.圖1給出了文中組合子空間算法2個信號到達角譜估計情況.由圖可以看出該組合子空間算法可以正確估計出多個捷變頻信號的DOA.

圖2 組合子空間DOA估計(10次)Fig.2 DOA estimation of the combined subspace algorithm(10 times)

4.2 預估計對捷變頻信號到達角估計性能影響

該方法在對付脈沖捷變頻信號時，由于需要對某一個脈沖內的信號進行預估計(通常選擇第一個到達脈沖)，所以有必要討論預估計測向誤差對該方法的測角性能進行討論.

設2個信號來波方向分別為θ1=10°和θ2=30°， SNR=10 dB，2個信號載頻在［500 MHz，600 MHz］任意跳變，參考頻率選擇為f=600 MHz.對預估計到達角分別為)θ1=7°和)θ2=27°、)θ1=5°和)θ2=25°以及)θ1=4°和)θ2=24°的3種情況該方法的測向準確性進行仿真.50次蒙特卡洛實驗(如圖3所示)結果表明在預估計偏離目標角度較大的情況，該方法可以較精確估計多部跳頻信號的到達角.只要預估計誤差在目標達到角(±5°)的波束范圍內，該方法可以有效估計到達角;如果預估計偏差加大，會增大測向誤差甚至測向錯誤.

圖3 預估計對測向方法的影響(50次)Fig.3 Performance of DOA estimation in different pre-estimation(50 times)

4.3 多捷變頻信號到達角估計精度的統(tǒng)計分析

兩信號來波方向分別為θ1=10°和θ2=30°，2個信號載頻在［500MHz，600MHz］任意跳變，參考頻率選擇為f=600 MHz，譜峰搜索步進角度為0.1°.均方根誤差RMSE(root mean squared error)隨信噪比變化如圖4所示.

圖4 DOA均方根誤差隨信噪比變化曲線Fig.4 RMSE of DOA estimation versus SNR

從圖4可以看出信噪比較高情況下，該算法的均方根誤差與經典MUSIC算法很接近，SNR在10 dB以上時基本保證0.1°以下，而SNR在0dB時也只有1°左右，符合被動測向精度要求.

5 結束語

基于陣列的目標波達方向估計都以估計導向矢量信息為基礎.捷變頻信號導向矢量隨頻率變化使得估計導向矢量變得困難.文章通過選擇參考頻率構造對應于同頻率導向矢量，簡化了導向矢量結構;利用同頻導向矢量的組合信號子空間結構實現(xiàn)了多個捷變頻信號的同時估計.相比于頻率-DOA聯(lián)合估計方法該方法在保證測向精度的同時，具有較小的計算復雜度適于實時處理的特點，適合于被動偵測系統(tǒng)對復雜雷達信號入射方向實現(xiàn)快速估計.

［1］謝機靈.超分辨測向算法在被動雷達導引頭系統(tǒng)中的應用［D］.哈爾濱:哈爾濱工程大學，2008:78-91.

［2］SCHMIT R O.Multiple emitter location and signal parameter estimation［J］.IEEE Trans on AP，1986，34(3):276-280.

［3］LI F，VACCARO R.On frequency-wavenumber estimation by state-space realization［J］.IEEE Trans on CAS，1991，38(3):800-804.

［4］LEMMA A N.Joint angle-frequency estimation using multiresolution ESPRIT［C］//Proceedings of the 1998 IEEE International Conference on Acoustics，Speech and Signal Processing.Seattle，USA，1998:1957-1960.

［5］吳云韜，廖桂生，田孝華.一種波達方向—頻率聯(lián)合估計快速算法［J］.電波科學學報，2003，18(4):380-384.

WU Yuntao，LIAO Guisheng，TIAN Xiaohua.Fast algorithm for joint direction-of-arrival and frequency estimation［J］.Chinese Journal of Radio Science，2003，18(4):380-384.

［6］茅以海.頻率捷變雷達［M］.北京:國防工業(yè)出版社，1981:27-31.

［7］HUNG H，KAVEH M.Focusing matrix for coherent signalsubspace processing［J］.IEEE Trans Accost，Speech，Signal Processing，1988(36):1272-1281.

［8］黃可生，黃知濤，周一宇.信號頻率和DOA聯(lián)合估計方法［J］.電子對抗，2003，6(6):1-4.

HUANG Kesheng，HUANG Zhitao，ZHOU Yiyu.An estimation method of the frequency and DOA of signals［J］.E-lectronic Warfare，2003，6(6):1-4.

［9］ZOLTOWSKI M D，MATHEWS C P.Real-time frequency and 2-D angle estimation with sub-Nyquist spatio-temporal sampling［J］.IEEE Trans Signals Processing，1994，42: 2781-2794.

［10］侯寶林，李陟，宋才水.捷變頻雷達信號頻域偵察截獲技術研究［J］.現(xiàn)代防御技術，2007，35(1):87-91.

HOU Baolin，LI Zhi，SONG Caishui.Study on frequencydomain intercept technique for frequency agility radar［J］.Modern Defence Technology，2007，35(1):87-91.