基于FFT插值的Keystone變換實現(xiàn)方法

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

(空軍預警學院預警技術系,武漢 430000)

0 引言

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

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

1 Keystone變換基本原理

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

(1)

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

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

(2)

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

2) 步驟2。對ys(f,tm)沿慢時間進行模糊數(shù)補償,得到

(3)

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

(4)

(5)

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

2 理想與實際去耦合過程

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

(6)

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

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

(7)

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

(8)

圖1 TS前、后信號時頻分布示意圖

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

(9)

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

圖2 實際的去耦合過程示意圖

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

為確保去耦后信號采樣點數(shù)相同,當αf>1時,需要對TS后信號時域截取前M個采樣點,當αf<1時,需要在TS后信號時域補零,補零個數(shù)為

(10)

ceil(0.05M)。

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

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

3.1 FFT插值

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

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

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

(11)

當N1為偶數(shù)時,令

(12)

(13)

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

圖4 頻域補零示意圖

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

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

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

(14)

令m=N2r/N1,得到

(15)

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

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

(16)

同理,當N1為偶數(shù)時,可得

(17)

同樣,令m=N2r/N1,進一步得到

(18)

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

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

(19)

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

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

(20)

當0<αf<1時,利用式(20)補零顯然是不合理的,因為Nz1(αf)為負整數(shù),可采取“頻域補零+時域抽取”的方法實現(xiàn),補零個數(shù)為

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

(21)

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

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

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

3.3 計算復雜度分析

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

(22)

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

(23)

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

(24)

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

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

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

4 仿真結果與分析

4.1 參數(shù)設置

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

4.2 可行性分析

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

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

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

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

圖9 時間尺度后回波

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

圖10 回波相參積累結果

4.3 抗噪效能分析

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

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

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

5 結束語

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