馮維婷, 梁 青, 谷 靜
(西安郵電大學(xué)電子工程學(xué)院, 陜西 西安 710121)
在雷達目標(biāo)中有一類是旋轉(zhuǎn)目標(biāo),比如無人機旋翼的旋轉(zhuǎn),彈道中段目標(biāo)的自旋、錐旋,這一類目標(biāo)除質(zhì)心平動以外內(nèi)部旋轉(zhuǎn)的局部運動稱為微動[1-3]。無人機處于懸停狀態(tài)時的微動特征是識別目標(biāo)或干擾的重要依據(jù),彈道中段目標(biāo)的自旋、錐旋特征也是區(qū)別于誘餌和欺騙干擾的重要特征。因此,對旋轉(zhuǎn)目標(biāo)微動特征的提取分析一直是雷達目標(biāo)識別領(lǐng)域的研究熱點[4-5]。
雷達目標(biāo)的微動會調(diào)制雷達回波相位,使得目標(biāo)主體平動引起的多普勒頻率周圍產(chǎn)生調(diào)制邊帶,這是微多普勒效應(yīng)[6]。對微多普勒參數(shù)的提取方法中,文獻[7-8]采用參數(shù)化提取方法,基于目標(biāo)的微動特征構(gòu)造包含若干個參數(shù)的解調(diào)算子對信號進行解調(diào),若干個參數(shù)在多維平面上進行搜索,當(dāng)合適的參數(shù)構(gòu)成的解調(diào)算子與信號匹配時解調(diào)信號能量最大,從而通過搜索峰值進行目標(biāo)微動參數(shù)的估計。參數(shù)化方法是多維搜索方法,參數(shù)估計精度取決于搜索步長,高精度的參數(shù)估計值必然有巨大的運算量。
微多普勒頻率具有時變特征,大量的文獻利用時頻分析方法處理這類非平穩(wěn)信號。文獻[9]對時頻分布結(jié)果采用經(jīng)驗?zāi)J椒纸夥椒▽r頻曲線分解成若干獨立曲線來估計目標(biāo)的微動參數(shù),但經(jīng)驗?zāi)J椒纸夥椒ù嬖谀J交殳B問題,無法有效處理復(fù)雜的雷達信號。文獻[10]計算信號的時頻分布,基于時頻脊線分離出各散射點的時頻曲線來取得目標(biāo)微動參數(shù),然而存在相互交疊處的時頻脊線提取錯誤問題。為了解決多條時頻曲線相互交疊難以提取的問題,文獻[11]結(jié)合時頻分布與逆Radon變換進行微動參數(shù)估計,逆Radon變換需要提前知道時頻曲線的周期,文中采用自相關(guān)法進行周期估計,但自相關(guān)法適合單分量信號,無法有效估計多分量的信號周期。
針對以上問題,以旋轉(zhuǎn)目標(biāo)微動參數(shù)提取為研究對象,利用德州儀器公司的60 GHz頻段調(diào)頻連續(xù)波(frequency modulated continuous wave, FMCW)毫米波小型雷達系統(tǒng)IWR6843搭建實驗平臺,進行旋轉(zhuǎn)目標(biāo)雷達實測數(shù)據(jù)的采集。從FMCW雷達接收到的回波信號出發(fā),分析了旋轉(zhuǎn)目標(biāo)微多普勒信號特征、利用奇異值比譜、時頻分析和逆Radon變換相結(jié)合的方法進行雷達采集數(shù)據(jù)的分析,完成了旋轉(zhuǎn)目標(biāo)微動參數(shù)的有效估計。
雷達發(fā)射寬帶線性調(diào)頻信號,在一個調(diào)頻周期內(nèi),發(fā)射信號的解析表達式為
(1)
式中:t表示一個調(diào)頻周期TC內(nèi)的快時間;AT表示信號振幅;fc表示信號載波頻率;K表示信號調(diào)頻斜率。經(jīng)時延τ后,雷達接收的回波信號表示為
(2)
式中:AR表示接收信號振幅;時延τ=2d/c,其中d表示雷達與目標(biāo)之間距離,c表示光速。
sR(t)與sT(t)混頻處理,再經(jīng)過低通濾波器,得到的中頻信號表示為
sIF(t)=AIFexp[j2π(fcτ+fIFt)]
(3)
式中:AIF表示中頻信號振幅;fIF表示中頻信號頻率;fIF=Kτ。fIF與距離d的關(guān)系為
(4)
式中:B表示信號帶寬。由于TC很小,一個TC內(nèi)的距離近似為常量,因此對式(3)進行頻域分析估得fIF,再利用式(4)可以得到目標(biāo)距離。
FMCW雷達連續(xù)發(fā)射多幀信號,每幀有多個線性調(diào)頻信號,對多幀中頻信號需要考慮目標(biāo)運動引起的多普勒現(xiàn)象,離散化后的多幀中頻信號可以表示為
sIF(n,m)=AIF(n,m)·exp{j2π[fIFn+fd(m)n)]}
(5)
式中:n表示快時間離散序列,長度為N;m表示慢時間離散序列,長度為M;fd是距離變化引起的多普勒頻率。
為了獲取旋轉(zhuǎn)目標(biāo)的微動參數(shù),對式(5)在每個TC進行脈沖壓縮,脈沖壓縮結(jié)果中目標(biāo)所在距離單元對應(yīng)的M點信號就是包含多普勒頻率的雷達多普勒信號。
圖1 雷達與旋轉(zhuǎn)目標(biāo)幾何關(guān)系示意圖Fig.1 Diagram of geometrical relationship between radar and rotating target
當(dāng)旋轉(zhuǎn)目標(biāo)滿足遠場條件,即(r/D)2→0,雷達與旋轉(zhuǎn)散射點P的距離可近似為
dP(t)≈D+rcos(2πfrt+θ0)
(6)
對式(6)求導(dǎo)得目標(biāo)速度為v(t)=d/dt[dP(t)]=-2πrfrsin(2πfrt+θ0)。根據(jù)fd的定義,旋轉(zhuǎn)目標(biāo)的微多普勒頻率可以表示為
(7)
式中:λ為載波波長。
由式(7)可見,旋轉(zhuǎn)散射點P的微多普勒頻率是隨時間按正弦規(guī)律變化的曲線。正弦曲線的頻率就是目標(biāo)的旋轉(zhuǎn)頻率fr,初始相位就是散射點P在觀測初始時刻的相位θ0,振幅即是最大微多普勒頻率,該值由旋轉(zhuǎn)半徑、旋轉(zhuǎn)頻率和載波波長共同決定,包含了目標(biāo)的結(jié)構(gòu)信息。一般情況下繞固定中心旋轉(zhuǎn)的散射點有多個,它們具有相同的旋轉(zhuǎn)頻率fr、不同的旋轉(zhuǎn)半徑和初始相位,Q個旋轉(zhuǎn)散射點的多普勒頻率可以表示為
(8)
可以看出,同一旋轉(zhuǎn)中心的多個散射點,其微動特征表現(xiàn)為微多普勒頻率是一組頻率相同的正弦信號的疊加,各正弦信號的振幅和初始相位不同。因此,旋轉(zhuǎn)目標(biāo)微多普勒頻率中各正弦信號振幅、頻率和初始相位3個參數(shù)的提取是檢測和識別不同旋轉(zhuǎn)目標(biāo)的基本依據(jù)。
根據(jù)第1節(jié)結(jié)論,旋轉(zhuǎn)目標(biāo)微動參數(shù)包括旋轉(zhuǎn)頻率fr、初始相位θ0和包含在最大微多普勒頻率中的旋轉(zhuǎn)半徑r。為了提取微動參數(shù),采用的方法主要有3部分:基于奇異值比譜的旋轉(zhuǎn)頻率估計、時頻分析技術(shù)和利用逆Radon變換的振幅和初始相位估計。
雷達脈沖壓縮后的多普勒信號表示為s(m),m=1,2,…,M。用信號s(m)構(gòu)造矩陣A,其表達式表示為
(9)
式中:Mj、Mi是正整數(shù),Mj·Mi≤M。對A進行奇異值分解[13]
A=UDVT
(10)
式中:U為Mj階矩陣,U的列為A·AT的正交特征向量;V為Mi階矩陣,V的列為AT·A的正交特征向量;D為Mj×Mi的對角矩陣,對角線上前z個非零元素稱為奇異值,表示為σi(i=1,2,…,z)。
實際上基于式(9)構(gòu)造的矩陣A存在假設(shè)周期與實際周期之間的截斷誤差不斷積累的問題,常常無法有效估計頻率。受到周期信號奇異值分解特點的啟發(fā),改進矩陣A的構(gòu)造方法。重新構(gòu)造的矩陣A為
(11)
式中:Mk=[k·fPRF/fri],符號[·]表示取整,fri是旋轉(zhuǎn)頻率在一定估計范圍內(nèi)進行搜索的第i次值,搜索范圍表示為[frmin,frmax];L=[fPRF/fri]。
對每一個fri,進行式(11)構(gòu)造的矩陣A的奇異值分解,計算奇異值比譜SVR(fri),旋轉(zhuǎn)頻率估計范圍[frmin,frmax]的峰值處對應(yīng)頻率就是fr的估計值。
時頻分析技術(shù)是分析非平穩(wěn)信號的有利工具,可以得到微多普勒頻率隨時間變化規(guī)律。這里采用同步壓縮變換(synchro-squeezing transform,SST)的時頻分析技術(shù)。SST以短時傅里葉變換(short time Fourier transform,STFT)為基礎(chǔ)計算信號的線性時頻譜,再進行壓縮,最后輸出分辨率高于STFT的時頻譜[15-16]。信號s(t)的STFT定義為
(12)
式中:g(τ-t)代表隨t移動的窗函數(shù)。用一定長度的窗函數(shù)g截取信號s(t),對截取信號進行傅里葉變換得到t時刻信號的局部頻譜信息。隨著窗函數(shù)g沿整個時間軸的移動,就得到了信號在整個時間-頻率軸上的完整頻譜分布。
對于諧波信號sk(t)=ej2πfkt,經(jīng)STFT后的時頻譜能量集中在以頻率值fk為中心的一定頻帶寬度內(nèi),并且其兩側(cè)的時頻系數(shù)隨著與fk的距離增大而減小。由于窗函數(shù)的截斷效應(yīng),信號的時頻譜存在模糊,時頻分辨率較低。
為了提高時頻分辨率,SST基于STFT進行瞬時頻率fI(t,f)的估計[17],表達式為
(13)
(14)
對旋轉(zhuǎn)目標(biāo)的多普勒信號進行時頻分析,微多普勒頻率表現(xiàn)為時頻譜圖中多條同頻率的正弦曲線的疊加,各曲線之間相互交疊。傳統(tǒng)的微多普勒參數(shù)提取方法是通過檢測時頻譜的能量峰值,提取出各條時頻脊線完成參數(shù)估計,但是相互交疊的曲線在交叉點處容易檢測到錯誤曲線,導(dǎo)致方法失效。為了解決時頻曲線難以提取的問題,根據(jù)旋轉(zhuǎn)目標(biāo)微多普勒頻率的正弦形式特點,采用逆Radon變換提取正弦信號的振幅和初始相位兩個參數(shù)。
逆Radon變換能對時頻譜圖中的正弦曲線進行能量積累、聚焦,將微動參數(shù)估計問題轉(zhuǎn)化為逆Radon變換后參數(shù)空間中的峰值檢測問題。
旋轉(zhuǎn)散射點P的微多普勒頻率重新寫成如下形式:
fd(l,θ)=δ(l-lPsin(θ+θ0))
(15)
式中:lP是正弦信號的振幅,與式(7)的對應(yīng)關(guān)系為lP=4πrfr/λ,即最大微多普勒頻率;θ是隨時間變化的旋轉(zhuǎn)相位,θ=2πfrt。
時頻平面上fd(l,θ) 的逆Radon變換[19]為
g(x,y)=
(16)
式中:kx=ηcosθ;ky=ηsinθ。
Q個散射點逆Radon變換結(jié)果為
(17)
由式(17)可知,時頻平面上的一條正弦曲線經(jīng)逆Radon變換后是x-y參數(shù)空間中的一個點(xi,yi),該點對應(yīng)于正弦曲線的振幅lPi和初始相位θ0i。lPi與θ0i估計值與點坐標(biāo)(xi,yi)的關(guān)系為
(18)
進行逆Radon變換前需要已知各正弦曲線的周期,其值根據(jù)上文中改進的奇異值比譜方法得到。
綜上,多個散射點的微動參數(shù)估計具體算法過程如下:
步驟 1將接收到的雷達中頻信號存放為二維數(shù)據(jù)矩陣形式,每行是一個調(diào)頻周期的N點快時間信號,每列是M個調(diào)頻周期構(gòu)成的慢時間信號。
步驟 2對快時間維進行脈沖壓縮得到距離維信息,抑制靜止目標(biāo)、干擾和噪聲,取出旋轉(zhuǎn)目標(biāo)所在的若干個距離單元并進行相加,得到目標(biāo)慢時間維的M點多普勒信號。
步驟 3用多普勒信號構(gòu)造矩陣A,并進行奇異值分解,見式(10)和式(11),基于奇異值比譜并搜索峰值估計旋轉(zhuǎn)頻率fr。
步驟 4對多普勒信號計算SST時頻譜,對時頻譜圖進行逆Radon變換,估計最強目標(biāo)分量的旋轉(zhuǎn)半徑和初始相位。
步驟 5利用估計值構(gòu)造最強分量信號,將其從原多普勒信號中剔除。
步驟 6重復(fù)步驟4和步驟5,直到估計出所有分量信號參數(shù)為止。
算法的處理流程如圖2所示。
圖2 算法流程圖Fig.2 Flow chart of the proposed algorithm
選取3個旋轉(zhuǎn)散射點進行仿真,雷達波長為0.005 m;多普勒信號的振幅為1;散射點的旋轉(zhuǎn)頻率為4.5 Hz;各散射點的旋轉(zhuǎn)半徑為r1=0.10 m,r2=0.15 m,r3=0.20 m;初始相位為θ01=0.780 rad,θ02=2.870 rad,θ03=4.970 rad。設(shè)置高斯白噪聲下多普勒信號的信噪比(signal to noise ratio,SNR)為0 dB。
圖3 奇異值比譜Fig.3 Singular value ratio spectrum
多普勒信號的STFT和SST時頻譜如圖4所示。
由圖4可以看出,時頻譜中3個旋轉(zhuǎn)散射點微多普勒頻率為標(biāo)準(zhǔn)的正弦曲線形式,各曲線相互交疊,曲線的頻率相同,振幅和起始相位不同。圖4(a)的STFT時頻譜時頻聚集性較差,頻率分辨率較低。圖4(b)的SST時頻譜譜線清晰,時頻分辨率較高。
圖5 逆Radon變換結(jié)果Fig.5 Inverse Radon transform
采用所提方法在SNR分別為-6 dB、-3 dB、0 dB、3 dB、6 dB時,各進行100點蒙特卡羅實驗,得到旋轉(zhuǎn)頻率的估值均為4.5 Hz, 3個散射點的旋轉(zhuǎn)半徑和起始相位平均估計值見表1和表2。由表1和表2可以看出,各參數(shù)估計誤差隨著SNR的增加而減小,并且在強噪聲環(huán)境下所提方法仍能得到旋轉(zhuǎn)半徑和起始相位的高精度估計值。
表1 旋轉(zhuǎn)半徑估計值
表2 起始相位估計值
實測實驗使用的FMCW雷達是德州儀器公司的IWR6843毫米波雷達系統(tǒng),配置了DCA1000高速數(shù)據(jù)采集卡,采用1發(fā)4收模式采集信號。試驗中雷達參數(shù)設(shè)置如下:載波頻率fc=62.5 GHz,采樣頻率FS=2 MHz,調(diào)頻斜率K=25.998 MHz/μs,一幀數(shù)據(jù)包含128個調(diào)頻周期數(shù),調(diào)頻周期TC=240 μs,一個調(diào)頻周期內(nèi)的采樣點數(shù)為256點,實際帶寬B=3.328 GHz,脈沖重復(fù)頻率為3 200 Hz。
實測實驗裝置如圖6所示,設(shè)置的轉(zhuǎn)臺中心與雷達之間的距離約1.2 m, 旋轉(zhuǎn)半徑約0.14 m,兩個旋轉(zhuǎn)目標(biāo)的相位差為π rad。
圖6 實驗裝置Fig.6 Experimental devices
接收到的雷達中頻信號中取一幀數(shù)據(jù),首先對每個調(diào)頻周期的256點快時間信號進行脈沖壓縮得到一維距離像,再對一幀數(shù)據(jù)的128個調(diào)頻周期數(shù)據(jù)在慢時間維進行傅里葉變換得到距離-多普勒像。一幀實測數(shù)據(jù)的距離-多普勒像如圖7所示。
圖7 一幀數(shù)據(jù)距離-多普勒像Fig.7 Range-Doppler imaging of one frame’s data
圖8 奇異值比譜Fig.8 Singular value ratio spectrum
對實測實驗數(shù)據(jù)進行SST變換,時頻分析結(jié)果如圖9所示。
圖9 旋轉(zhuǎn)目標(biāo)SST時頻譜Fig.9 SST time-frequency spectrum of rotating targets
圖9的時頻譜中是兩條相位相反的正弦曲線,這與圖6的實驗設(shè)置相符。圖中零多普勒頻率對應(yīng)的轉(zhuǎn)盤主體能量大于兩個散射點的能量,不利于微動參數(shù)提取。對多普勒信號求導(dǎo)處理,消除強零頻干擾,而兩個散射點的微多普勒頻率波形基本不受影響,求導(dǎo)處理后的時頻譜如圖10所示。
圖10 求導(dǎo)處理后的SST時頻譜Fig.10 SST time-frequency spectrum after calculation of derivation
實測實驗數(shù)據(jù)的逆Radon變換結(jié)果如圖11所示。
圖11 逆Radon變換結(jié)果Fig.11 Inverse Radon transform
圖11中逆Radon變換結(jié)果存在展寬現(xiàn)象,這是因為圖10中的時頻曲線并非理想等幅的標(biāo)準(zhǔn)正弦曲線的緣故。
2個旋轉(zhuǎn)目標(biāo)的微動參數(shù)平均估計值見表3。實際旋轉(zhuǎn)目標(biāo)的初始相位和旋轉(zhuǎn)頻率未知,但從表3的旋轉(zhuǎn)半徑估值與真值的對比可以看出旋轉(zhuǎn)半徑的估值與真值基本符合;雖然初始相位未知,但是2個旋轉(zhuǎn)目標(biāo)的相位差是確定的,相位差的估值與真值也基本符合;旋轉(zhuǎn)頻率的估值是其他微動參數(shù)估計的基礎(chǔ),其他微動參數(shù)的有效估計間接表明了旋轉(zhuǎn)頻率估值的有效性,因此所提方法能有效估計實際雷達旋轉(zhuǎn)目標(biāo)的微動參數(shù)。
表3 實測數(shù)據(jù)微動參數(shù)估計值
以旋轉(zhuǎn)目標(biāo)為研究對象,基于FMCW雷達推導(dǎo)了其多普勒信號的解析形式,并根據(jù)信號時頻譜的正弦曲線特征,提出了綜合利用奇異值比譜、SST時頻分析技術(shù)和逆Radon變換以提取旋轉(zhuǎn)目標(biāo)微動參數(shù)的方法。仿真分析和雷達實測數(shù)據(jù)的實驗結(jié)果表明,所提方法能有效提取旋轉(zhuǎn)目標(biāo)旋轉(zhuǎn)頻率、旋轉(zhuǎn)半徑和初始相位值,算法抗干擾能力強,參數(shù)估計精度較高,所得的微動參數(shù)為后繼的雷達目標(biāo)分類、識別、監(jiān)測打下了基礎(chǔ),也為復(fù)雜背景中復(fù)雜雷達目標(biāo)微動特征分析提取提供了一種途徑。