張 亮, 杜慶磊, 胡 冰, 周必雷, 王安樂

(空軍預(yù)警學(xué)院預(yù)警技術(shù)系,武漢 430000)

0 引言

雷達(dá)作為獲取目標(biāo)位置信息的傳感器,在戰(zhàn)場指揮引導(dǎo)、態(tài)勢感知、武器制導(dǎo)等領(lǐng)域應(yīng)用廣泛。雷達(dá)信號(hào)處理中,經(jīng)常面臨目標(biāo)跨距離門、損失信噪比(SNR)問題[1],特別是目標(biāo)高速機(jī)動(dòng)或雷達(dá)處于長時(shí)間相參積累模式時(shí)。Keystone變換(KT)是一種經(jīng)典的雷達(dá)目標(biāo)距離走動(dòng)校正工具,適用于脈沖多普勒、捷變頻、合成孔徑等多種雷達(dá)體制[2-7]。KT實(shí)現(xiàn)過程包含2個(gè)關(guān)鍵環(huán)節(jié)[8-10],即模糊數(shù)補(bǔ)償和回波慢時(shí)間去耦合。第1個(gè)環(huán)節(jié)涉及真實(shí)模糊數(shù)估計(jì)問題,通常根據(jù)積累后的目標(biāo)峰值,采取循環(huán)搜索取極大值的方法得到;第2個(gè)環(huán)節(jié)涉及“構(gòu)造虛擬慢時(shí)間”,通常采取對(duì)快時(shí)間頻率單元回波進(jìn)行時(shí)間尺度(TS)實(shí)現(xiàn)。所謂TS,即連續(xù)信號(hào)到其尺度版本的映射過程,通俗地講,TS是一種對(duì)信號(hào)進(jìn)行某種程度拉伸或者壓縮的操作,由于TS數(shù)值計(jì)算方法有多種,現(xiàn)有KT實(shí)現(xiàn)方法研究主要圍繞該問題展開。文獻(xiàn)[8]最早提出基于辛格插值的TS計(jì)算方法,由于使用了時(shí)域內(nèi)插函數(shù),計(jì)算復(fù)雜度較高;為降低算法計(jì)算復(fù)雜度,文獻(xiàn)[9]提出基于Chirp濾波的TS計(jì)算方法,但為實(shí)現(xiàn)1個(gè)快時(shí)間頻率單元回波去耦合需要進(jìn)行2次Chirp乘積、2次Chirp卷積,計(jì)算量同樣較大;為降低算法計(jì)算量,文獻(xiàn)[10]進(jìn)一步提出基于“CZT+IFFT”的計(jì)算方法,其中,CZT為Chirp-Z變換(Chirp-Z Transform),IFFT為快速傅里葉逆變換(Inverse Fast Fourier Transform),該方法為實(shí)現(xiàn)1個(gè)快時(shí)間頻率單元回波去耦合需要進(jìn)行2次Chirp乘積、1次Chirp卷積和1次IFFT,計(jì)算量較文獻(xiàn)[9]低,但由于回波距離門數(shù)量通常很大,算法整體計(jì)算量同樣不容樂觀。TS會(huì)導(dǎo)致信號(hào)帶寬的變化,當(dāng)TS后信號(hào)帶寬增大時(shí)需要通過插值提升采樣頻率,而文獻(xiàn)[10]缺乏該操作,去耦過程中容易出現(xiàn)頻譜混疊,影響KT抗噪效能。另外,還有“DFT+IFFT”方法,DFT為離散傅里葉變換(Discrete Fourier Transform),計(jì)算過程與“CZT+IFFT”方法相同[11]。綜上可知,現(xiàn)有KT實(shí)現(xiàn)方法在計(jì)算量和抗噪效能方面還不夠理想,存在一定的改進(jìn)空間?？紤]到,辛格插值只是諸多插值算法中的一種,此外還有線性插值、高斯插值、FFT插值、三次樣條插值和變換域插值等[12],尋找一種低計(jì)算量的插值算法,結(jié)合時(shí)域截取、補(bǔ)零等操作完全可以解決TS的快速計(jì)算問題,降低KT整體計(jì)算量。

基于上述分析,提出一種基于快速傅里葉變換(FFT)插值的KT實(shí)現(xiàn)方法,利用FFT插值算法解決KT去耦合問題。KT中TS要求非整數(shù)倍插值,如果采取“整數(shù)倍插值—抽取”的方式實(shí)現(xiàn),計(jì)算量同樣不容樂觀。針對(duì)該問題,本文分析了FFT插值算法頻域補(bǔ)零個(gè)數(shù)對(duì)插值后序列的影響,給出非整數(shù)倍插值的高效計(jì)算方法,并進(jìn)行了舉例說明。

1 Keystone變換基本原理

設(shè)雷達(dá)載頻為f0,發(fā)射信號(hào)為線性調(diào)頻脈沖信號(hào),脈寬為Tp,帶寬為B,調(diào)頻斜率為ρ=B/Tp,雷達(dá)探測范圍內(nèi)1個(gè)點(diǎn)目標(biāo)向站飛行,初始距離為R0,徑向速度為vt。設(shè)目標(biāo)反射系數(shù)為1,脈壓后的目標(biāo)回波為[8-11]

(1)

式中:t為快時(shí)間;tm=mTr,為慢時(shí)間,m=0,1,2,…,M-1,M為相參積累個(gè)數(shù),Tr為脈沖重復(fù)周期;c為光速。利用KT進(jìn)行目標(biāo)距離走動(dòng)校正,具體包含以下4步。

1) 步驟1。對(duì)脈壓后回波沿快時(shí)間做傅里葉變換,得到

(2)

式中:f為快時(shí)間頻率;fd為不模糊多普勒頻率;F為模糊數(shù);fr=1/Tr,為脈沖重復(fù)頻率。

2) 步驟2。對(duì)ys(f,tm)沿慢時(shí)間進(jìn)行模糊數(shù)補(bǔ)償,得到

(3)

3) 步驟3。構(gòu)造“虛擬”慢時(shí)間tom=[(f+f0)/f0]tm,代入式(3),去除f與tm的耦合關(guān)系,得到

(4)

(5)

對(duì)式(5)沿tom進(jìn)行FFT可實(shí)現(xiàn)相參積累。步驟2涉及模糊數(shù)估計(jì)問題,通常根據(jù)相參積累后的目標(biāo)最大峰值搜索確定。對(duì)于寬帶雷達(dá),目標(biāo)為多散射點(diǎn)模型,KT實(shí)現(xiàn)環(huán)節(jié)相同[13-14]。

2 理想與實(shí)際去耦合過程

KT包含2個(gè)重要環(huán)節(jié),即模糊數(shù)補(bǔ)償和去耦合,其中,第2個(gè)環(huán)節(jié)最為關(guān)鍵,現(xiàn)有研究多強(qiáng)調(diào)構(gòu)造“虛擬”慢時(shí)間,該闡述較為模糊。為深入理解KT去耦合過程,將式(3)表示為

(6)

式中:φ(f)為快時(shí)間頻率項(xiàng);φ(αf,tm)為耦合項(xiàng),αf=(f+f0)/f0,為頻率尺度因子,對(duì)于雷達(dá)方αf是已知的。φ(αf,tm)可視為不同頻率的復(fù)指數(shù)信號(hào),頻率為αffd,而KT去耦合本質(zhì)是利用已知的αf,將αf≠1時(shí)的φ(αf,tm)映射為φ(1,tm)。如果不模糊多普勒頻率fd已知,對(duì)φ(αf,tm)移頻處理,可將αf≠1時(shí)的φ(αf,tm)調(diào)整為φ(1,tm),即

φ(1,tm)=φ(αf,tm)ei2π[(1-αf)fd]tm。

(7)

當(dāng)fd未知時(shí),需采取時(shí)間尺度(TS)的方法去耦合[9]。所謂TS,即連續(xù)信號(hào)x(t)到y(tǒng)(t)的映射過程[15],可表示為

(8)

圖1 TS前、后信號(hào)時(shí)頻分布示意圖

假設(shè)φ(αf,tm)為慢時(shí)間的連續(xù)信號(hào),令α=1/αf,x(t)=φ(αf,tm),代入式(8)可得

(9)

由式(9)可知,TS操作可將φ(αf,tm)映射為φ(1,tm)實(shí)現(xiàn)去耦合。設(shè)αf為0.5,1.0和2,相參積累個(gè)數(shù)為16,圖2為實(shí)際的去耦合過程示意圖。當(dāng)αf為0.5時(shí),尺度因子α為2,TS后信號(hào)長度為TS前信號(hào)的1/2,當(dāng)αf為2時(shí),尺度因子α為0.5,TS后信號(hào)長度為TS前信號(hào)的2倍。

圖2 實(shí)際的去耦合過程示意圖

圖3給出了基于TS的Keystone變換實(shí)現(xiàn)思路,難點(diǎn)在于TS數(shù)值計(jì)算。

為確保去耦后信號(hào)采樣點(diǎn)數(shù)相同,當(dāng)αf>1時(shí),需要對(duì)TS后信號(hào)時(shí)域截取前M個(gè)采樣點(diǎn),當(dāng)αf<1時(shí),需要在TS后信號(hào)時(shí)域補(bǔ)零,補(bǔ)零個(gè)數(shù)為

(10)

ceil(0.05M)。

3 基于快速FFT插值的Keystone變換去耦合算法原理

為實(shí)現(xiàn)KT去耦合,需要解決TS數(shù)值計(jì)算問題。由于TS會(huì)改變信號(hào)時(shí)寬,插值不可避免,尋找一種低計(jì)算量插值算法,可以解決TS快速計(jì)算問題。FFT插值是一種高效的插值算法,但為整數(shù)倍插值,KT中的TS要求非整數(shù)倍插值,本章將針對(duì)該問題給出解決方案。

3.1 FFT插值

FFT插值是一種整數(shù)倍插值算法[16-17],整個(gè)過程僅使用了1次FFT和1次IFFT。設(shè)插值前序列x(n)為實(shí)序列,長度為N1,插值后序列y(m)同樣為實(shí)序列,長度為N2,R為插值倍數(shù),N2=N1R,FFT插值具體步驟如下。

1) 步驟1。計(jì)算x(n)的N1點(diǎn)FFT,得到XN1(k)。

2) 步驟2。創(chuàng)建長度為N2的序列YN2(l),當(dāng)N1為奇數(shù)時(shí),令

(11)

當(dāng)N1為偶數(shù)時(shí),令

(12)

(13)

式中,real[·]表示取實(shí)部。步驟2為FFT插值算法核心,為直觀顯示,設(shè)N1分別為8和9,插值倍數(shù)為2,圖4為頻域補(bǔ)零示意圖,圖中,紅色標(biāo)記為正頻率,綠色標(biāo)記為負(fù)頻率,黑色標(biāo)記為截止頻率,藍(lán)色標(biāo)記為補(bǔ)零頻點(diǎn)。

圖4 頻域補(bǔ)零示意圖

由上述步驟可知,FFT插值具體操作為高頻補(bǔ)零,等效為低通濾波,與常用辛格插值相似,補(bǔ)零個(gè)數(shù)為N2-N1(N1為奇數(shù))和N2-N1-1(N1為偶數(shù))。由步驟1～3還可知,為得到插值后的序列,需要進(jìn)行1次N1點(diǎn)FFT和1次N2點(diǎn)IFFT,總共需要的復(fù)乘次數(shù)為0.5N1lbN1+0.5N2lbN2,計(jì)算復(fù)雜度為O(N2lbN2),低于文獻(xiàn)[8]方法。另外,FFT插值還有一種計(jì)算復(fù)雜度更低的子序列實(shí)現(xiàn)方法[18],基本思想是將1次N2點(diǎn)的IFFT拆分成R-1次的N1點(diǎn)IFFT,總共需要(R-1)(0.5N1lbN1+2N1)次復(fù)乘運(yùn)算,計(jì)算復(fù)雜度為O(N2lbN1),在此不做詳述。YN2(l)實(shí)際上并不是原始信號(hào)以Rfs為采樣頻率計(jì)算得到的頻譜,因此存在插值誤差,可對(duì)補(bǔ)零后頻譜前后幾個(gè)點(diǎn)加窗以減少誤差[16]。

3.2 利用改進(jìn)的FFT插值算法實(shí)現(xiàn)一維信號(hào)時(shí)間尺度操作

3.1節(jié)介紹了FFT插值基本原理,經(jīng)典FFT插值設(shè)定插值前序列為實(shí)序列、插值倍數(shù)為整數(shù),而KT耦合項(xiàng)為解析信號(hào)、尺度因子為1附近的小數(shù),因此利用FFT插值計(jì)算TS需要解決解析信號(hào)適用性及非整數(shù)倍插值問題。對(duì)于第1個(gè)問題,易知解析信號(hào)DFT為單邊譜,而實(shí)序列DFT為雙邊譜,利用式(11)、式(12)進(jìn)行補(bǔ)零同樣可行,僅是當(dāng)N1為偶數(shù)時(shí),將YN2(N2-N1/2)置零即可,因?yàn)槭?12)中對(duì)l=N2-N1/2的賦值是為了滿足實(shí)值序列頻譜的共軛對(duì)稱性;對(duì)于第2個(gè)問題,需要知道插值后序列表示式,計(jì)算式(11)、式(12)的離散傅里葉逆變換。當(dāng)N1為奇數(shù)時(shí),可得

(14)

令m=N2r/N1,得到

(15)

將式(15)代入式(13)可得插值后的序列表示式為

y(rN2/N1)=y(Rr)=x(r)。

(16)

同理,當(dāng)N1為偶數(shù)時(shí),可得

(17)

同樣,令m=N2r/N1,進(jìn)一步得到

(18)

將式(18)代入式(13)可得插值后的序列表示式為

y(rN2/N1)=y(Rr)≈x(r)。

(19)

綜合式(16)、式(19)可知,當(dāng)插值倍數(shù)為整數(shù)時(shí),插值后的第Rr值與插值前的第r值近似一致;當(dāng)插值倍數(shù)非整數(shù)時(shí),插值后的第N2r/N1值與插值前的第r值近似一致,FFT插值可理解為TS的離散形式,尺度因子α為N1/N2(對(duì)應(yīng)αf=N2/N1),通過改變頻域補(bǔ)零個(gè)數(shù)可實(shí)現(xiàn)TS。當(dāng)αf>1時(shí),頻域補(bǔ)零個(gè)數(shù)取

Nz1(αf)=ceil(αfN1-N1)。

(20)

當(dāng)0<αf<1時(shí),利用式(20)補(bǔ)零顯然是不合理的,因?yàn)镹z1(αf)為負(fù)整數(shù),可采取“頻域補(bǔ)零+時(shí)域抽取”的方法實(shí)現(xiàn),補(bǔ)零個(gè)數(shù)為

Nz2(αf)=ceil(LαfN1-N1)

(21)

式中,L為抽取倍數(shù),L=ceil(1/αf),對(duì)于窄帶雷達(dá),L=2。為直觀顯示,圖5給出了基于FFT插值的TS實(shí)現(xiàn)流程,結(jié)合圖3可得完整的KT實(shí)現(xiàn)過程。

圖5 基于FFT插值的TS實(shí)現(xiàn)流程

圖5涉及頻域補(bǔ)零、L倍抽取問題,而圖3又涉及時(shí)域截取、補(bǔ)零,在此舉例說明。設(shè)TS前信號(hào)即耦合項(xiàng)φ(αf,tm)長度為256,256同樣為相參積累個(gè)數(shù),頻率尺度因子αf為1.05,根據(jù)式(20)可得頻域補(bǔ)零個(gè)數(shù)為13,進(jìn)而得到TS后信號(hào)φ(1,tm)長度為269,由于理想KT去耦合過程要求去耦后信號(hào)數(shù)字采樣點(diǎn)數(shù)相同,截取TS后信號(hào)前256個(gè)點(diǎn);同理,設(shè)頻率尺度因子αf為0.95,根據(jù)式(21)可得抽取倍數(shù)L=ceil(1/αf)=2,頻域補(bǔ)零個(gè)數(shù)為231,進(jìn)而得到TS后信號(hào)φ(1,tm)長度為487,再進(jìn)行2倍抽取(長度為243),由于理想KT去耦合過程要求去耦后信號(hào)數(shù)字采樣點(diǎn)數(shù)相同,對(duì)抽取后的信號(hào)補(bǔ)13個(gè)零。

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

現(xiàn)有KT實(shí)現(xiàn)方法主要包括辛格插值法[8]、Chirp濾波法[9]和“CZT+IFFT”方法[10],其中,“CZT+IFFT”方法計(jì)算量最低。本節(jié)對(duì)比“CZT+IFFT”方法,對(duì)所提方法計(jì)算復(fù)雜度進(jìn)行分析。設(shè)回波信號(hào)經(jīng)下變頻、脈沖壓縮、快時(shí)間FFT,得到M×P的二維矩陣,M為積累脈沖數(shù),P為快時(shí)間頻率點(diǎn)數(shù)(奇數(shù))?！癈ZT+IFFT”方法去耦合過程包含CZT和IFFT兩步,共需(P-1)[M1+2M+1.5M1lbM1+0.5MlbM]次復(fù)乘[9],M1為大于2M-1的正整數(shù),計(jì)算復(fù)雜度為O(PM1lbM1)。所提方法包含3步,分別為FFT、頻域補(bǔ)零和IFFT,對(duì)于1個(gè)快時(shí)間頻率單元回波,計(jì)算FFT均需要進(jìn)行0.5MlbM次復(fù)數(shù)乘法,而計(jì)算不同快時(shí)間頻率單元回波IFFT所需的復(fù)乘次數(shù),與補(bǔ)零后的頻譜長度有關(guān)。由式(20),(21)可知,當(dāng)αf>1時(shí),補(bǔ)零后的頻譜長度為Nz1(αf)+N1,當(dāng)0<αf<1時(shí),補(bǔ)零后的頻譜長度為Nz2(αf)+N1。由于αf=(f+f0)/f0是關(guān)于快時(shí)間頻率f的連續(xù)信號(hào),將其表示成離散形式為

(22)

式中,p=-(P-1)/2,…,0,…,(P-1)/2。所提算法總共需要的復(fù)乘次數(shù)Nnum為

(23)

Nnum<0.5(P-1)MlbM+(P-1)Mlb (2M)。

(24)

所提方法計(jì)算復(fù)雜度為O[PMlb (2M)],低于“CZT+IFFT”方法。設(shè)P分別為2001,4001,αf取值0.95～1.05,相參積累個(gè)數(shù)為16～1000,圖6給出了復(fù)乘次數(shù)隨積累脈沖數(shù)變化曲線,相同快時(shí)間頻率點(diǎn)數(shù)下所提方法需要的復(fù)乘次數(shù)更少。

圖6 復(fù)乘次數(shù)隨積累脈沖數(shù)變化曲線

直觀上看,“CZT+IFFT”方法為實(shí)現(xiàn)1個(gè)快時(shí)間頻率單元回波去耦合需要進(jìn)行2次Chirp乘積、1次Chirp卷積和1次IFFT,所提方法僅需1次FFT和1次IFFT(不考慮頻域補(bǔ)零、時(shí)域截取、時(shí)域補(bǔ)零產(chǎn)生的運(yùn)算量),計(jì)算量明顯更低。

4 仿真結(jié)果與分析

4.1 參數(shù)設(shè)置

雷達(dá)載頻為0.5 GHz,脈沖重復(fù)頻率為10 kHz,相參積累個(gè)數(shù)為256,發(fā)射信號(hào)為LFM脈沖信號(hào),脈寬為15 μs,帶寬為20 MHz,采樣頻率為40 MHz,雷達(dá)探測范圍內(nèi)1個(gè)高速點(diǎn)目標(biāo)向站勻速飛行,初始距離為8 km,不模糊速度為450 m/s,模糊數(shù)為41。

4.2 可行性分析

對(duì)所提基于FFT插值的KT可行性進(jìn)行驗(yàn)證,為顯示細(xì)節(jié),仿真回波中不考慮噪聲,下一節(jié)將進(jìn)行算法抗噪效能分析。首先,對(duì)回波快時(shí)間脈沖壓縮結(jié)果如圖7所示,目標(biāo)存在明顯距離走動(dòng)。

其次,對(duì)脈壓后回波快時(shí)間做FFT,沿慢時(shí)間進(jìn)行模糊數(shù)補(bǔ)償(真實(shí)模糊數(shù)),再沿快時(shí)間頻率做IFFT,觀察“模糊數(shù)補(bǔ)償”環(huán)節(jié)對(duì)目標(biāo)距離走動(dòng)校正的影響,結(jié)果如圖8所示,目標(biāo)慢時(shí)間上仍未完全對(duì)齊。

圖8 模糊數(shù)補(bǔ)償后回波

然后,對(duì)圖8回波快時(shí)間做FFT,沿慢時(shí)間進(jìn)行TS,再沿快時(shí)間頻率做IFFT,完成目標(biāo)距離走動(dòng)校正,結(jié)果如圖9所示,目標(biāo)在慢時(shí)間上完全對(duì)齊。上述涉及模糊數(shù)估計(jì)問題,可以采取循環(huán)搜索最大峰值的方式得到。

圖9 時(shí)間尺度后回波

最后,對(duì)圖7、圖9分別沿慢時(shí)間做FFT完成相參積累結(jié)果,如圖10所示,校正前無法檢測目標(biāo),校正后目標(biāo)得到充分聚焦,峰值搜索可得目標(biāo)初始距離為8 km,不模糊速度為453 m/s,與仿真參數(shù)基本一致。由于圖10已對(duì)回波進(jìn)行了增益歸一化,目標(biāo)理論幅值應(yīng)接近1 V,而實(shí)際得到的目標(biāo)幅值為0.82 V,說明KT對(duì)目標(biāo)存在一定的處理損失(現(xiàn)有方法均如此,且不可避免),下面將針對(duì)該問題,對(duì)不同的KT實(shí)現(xiàn)方法處理損失及能量聚焦性進(jìn)行評(píng)估。

圖10 回波相參積累結(jié)果

4.3 抗噪效能分析

對(duì)基于FFT插值的KT實(shí)現(xiàn)方法抗噪效能進(jìn)行分析,使用的評(píng)估指標(biāo)為輸出信噪比(SNR-out)和瑞利熵(RE)[19],前者用于評(píng)估不同KT實(shí)現(xiàn)方法對(duì)目標(biāo)的處理損失程度,后者用于評(píng)估相參積累后回波能量聚焦程度。相同輸入信噪比(SNR-in)條件下,SNR-out越大,熵值越小,說明KT對(duì)目標(biāo)的處理損失越小,目標(biāo)聚焦性越好。另外,仿真計(jì)算SNR-out時(shí),已對(duì)回波脈壓增益、相參積累增益進(jìn)行了歸一化,且下文中的參數(shù)設(shè)置能夠確保目標(biāo)被有效檢測。設(shè)相參積累個(gè)數(shù)為128,SNR-in取值-15～5 dB,間隔2 dB,運(yùn)行蒙特卡羅仿真500次,KT抗噪效能曲線如圖11所示。設(shè)目標(biāo)不模糊徑向速度為900 m/s,KT抗噪效能曲線如圖12所示,進(jìn)一步設(shè)雷達(dá)載頻為1 GHz,KT抗噪效能曲線如圖13所示?？梢钥闯?圖11(a),12(a)和13(a)中的黑色線條(對(duì)應(yīng)所提方法)多位于其他線條的上方,說明所提方法處理后的目標(biāo)峰值高于文獻(xiàn)[8-10]方法;圖11(b),12(b)和13(b)中的黑色線條多位于其他線條的下方,說明所提方法處理后的回波能量聚焦性更強(qiáng)。

圖11 KT抗噪效能曲線(回波脈沖數(shù)為128)

圖12 KT抗噪效能曲線(目標(biāo)不模糊徑向速度為900 m/s)

5 結(jié)束語

針對(duì)現(xiàn)有Keystone變換(KT)實(shí)現(xiàn)方法計(jì)算量和抗噪效能不夠理想的問題,提出基于FFT插值的KT實(shí)現(xiàn)方法。仿真結(jié)果表明,所提方法能夠有效校正目標(biāo)距離走動(dòng),抗噪效能優(yōu)于現(xiàn)有方法,且計(jì)算復(fù)雜度更低。另外,所提算法使用了FFT,IFFT,頻域補(bǔ)零及時(shí)域抽取等基礎(chǔ)操作,存在工程實(shí)現(xiàn)可能。本文利用FFT插值實(shí)現(xiàn)時(shí)間尺度操作,解決了KT去耦合問題,考慮到呂分布、寬帶模糊函數(shù)及合成孔徑雷達(dá)成像等同樣面臨時(shí)間尺度數(shù)值計(jì)算問題,所提方法均可應(yīng)用。