金家偉， 阮懷林

(國防科技大學(xué)電子對(duì)抗學(xué)院， 合肥， 230031)

在彈道中段的飛行過程中，為了保持飛行的穩(wěn)定性，彈頭一般會(huì)進(jìn)行姿態(tài)控制，不僅會(huì)繞自身對(duì)稱軸作自旋運(yùn)動(dòng)，還會(huì)受到?jīng)_擊力矩的作用，力矩消失后對(duì)稱軸將在平衡位置做圓錐運(yùn)動(dòng)，即進(jìn)動(dòng)[1]。而誘餌一般不具備姿態(tài)控制裝置，因此，進(jìn)動(dòng)是中段彈頭特有的微動(dòng)方式，與目標(biāo)形狀結(jié)構(gòu)、質(zhì)量分布以及動(dòng)力學(xué)特性等物理屬性密切相關(guān)[2]。由微動(dòng)引起的雷達(dá)回波中多普勒頻率的調(diào)制，一般被稱為微多普勒效應(yīng)[3]。微多普勒特征體現(xiàn)出目標(biāo)獨(dú)特的動(dòng)態(tài)和結(jié)構(gòu)特征，在目標(biāo)識(shí)別和分類方面具有巨大的潛力[4-5]。

對(duì)于非剛性雷達(dá)目標(biāo)，目標(biāo)主體散射的雷達(dá)回波占返回能量的絕大部分，微動(dòng)信號(hào)始終浸沒在主體雷達(dá)回波中，需要將雷達(dá)回波分解為平動(dòng)分量和微動(dòng)分量，用微動(dòng)分量來估計(jì)微多普勒參數(shù)[6-8]。文獻(xiàn)[9]對(duì)信號(hào)進(jìn)行小波變換，用分解出的微動(dòng)分量的自相關(guān)來估計(jì)振動(dòng)/旋轉(zhuǎn)速率；文獻(xiàn)[10]通過信號(hào)分解提取單個(gè)散射點(diǎn)特征。而為了表示出散射點(diǎn)雷達(dá)回波的時(shí)變頻率調(diào)制特性，時(shí)頻分析被廣泛用于分析微多普勒特征，基于時(shí)頻分布的微動(dòng)參數(shù)估計(jì)已成為微多普勒特征提取的重要部分。文獻(xiàn)[11]采用基于加權(quán)迭代自適應(yīng)的時(shí)頻分析方法結(jié)合逆Radon變換，分離重構(gòu)不同散射點(diǎn)的微多普勒分量；文獻(xiàn)[12]利用短時(shí)分?jǐn)?shù)階傅里葉變換分離肢體和軀干的微多普勒信號(hào)；文獻(xiàn)[13]結(jié)合Gabor時(shí)頻分布和變分模態(tài)分解，來估計(jì)目標(biāo)的自旋頻率和錐旋頻率。

上述方法都需要假設(shè)目標(biāo)平動(dòng)已被準(zhǔn)確補(bǔ)償，在此前提下，目標(biāo)時(shí)頻分布中的微多普勒頻率曲線通常應(yīng)該表現(xiàn)為正弦曲線[3]或是多個(gè)正弦分量的合成[2]?，F(xiàn)有的微多普勒特征提取技術(shù)，如經(jīng)驗(yàn)?zāi)Ｊ椒纸?empirical-mode decomposition, EMD)[14]和擴(kuò)展Hough變換(extended hough transform,EHT)[15]，主要依賴于這種正弦調(diào)制假設(shè)。但目標(biāo)運(yùn)動(dòng)由平動(dòng)與微動(dòng)復(fù)合而成，特別是對(duì)于高速運(yùn)動(dòng)的彈道目標(biāo)，周期進(jìn)動(dòng)引起的彈道目標(biāo)微多普勒效應(yīng)并不像通常那樣表現(xiàn)為正常的正弦調(diào)制[16]。

針對(duì)目標(biāo)微動(dòng)伴隨平動(dòng)時(shí)，時(shí)頻分布不再表現(xiàn)為正弦調(diào)制曲線，本文在雷達(dá)回波信號(hào)時(shí)頻分布的基礎(chǔ)上，結(jié)合循環(huán)自相關(guān)函數(shù)(circular autocorrelation function，CACF)和循環(huán)平均幅度差函數(shù)(circular average amplitude difference function，CAMDF)，以獲取時(shí)頻分布的循環(huán)系數(shù)矩陣，并通過該矩陣的平均循環(huán)系數(shù)估計(jì)出彈道目標(biāo)的微動(dòng)周期。

1 模型建立

空間進(jìn)動(dòng)錐體目標(biāo)模型如圖1所示，坐標(biāo)系(U,V,W)為雷達(dá)坐標(biāo)系，雷達(dá)靜止于原點(diǎn)Q。O為目標(biāo)質(zhì)心，以O(shè)為原點(diǎn)、目標(biāo)對(duì)稱軸Oz為z軸建立目標(biāo)本體坐標(biāo)系O-xyz。以O(shè)為原點(diǎn)建立參考坐標(biāo)系O-XYZ，以初始時(shí)刻與目標(biāo)對(duì)稱軸Oz、進(jìn)動(dòng)軸OZ共面且垂直于OZ的方向?yàn)閅軸，X軸根據(jù)右手準(zhǔn)則確定。目標(biāo)在平動(dòng)的同時(shí)，以角速度ωs繞對(duì)稱軸z軸做自旋運(yùn)動(dòng)，同時(shí)以角速度ωc繞軸OZ做錐旋運(yùn)動(dòng)(ωs和ωc均采用參考坐標(biāo)系中的表達(dá)式)，自旋軸和錐旋軸之間的夾角為進(jìn)動(dòng)角θ。

圖1 目標(biāo)模型

在光學(xué)區(qū)，雷達(dá)目標(biāo)的整體散射特性通?？梢缘刃槿舾蓚€(gè)散射中心的疊加。為不失一般性，我們假設(shè)目標(biāo)的微動(dòng)是周期性進(jìn)動(dòng)，目標(biāo)的平動(dòng)是任意且未知的。假設(shè)目標(biāo)等效為K個(gè)散射中心，各散射中心各向同性，雷達(dá)發(fā)射的電磁波為連續(xù)單頻波，載頻為f0?？紤]目標(biāo)的平動(dòng)，第i個(gè)散射中心與雷達(dá)之間的距離為：

ri(t)≈rT(t)+rMi(t)

(1)

式中：i=1,2，…，K，為散射中心序號(hào)；rT(t)為目標(biāo)平動(dòng)對(duì)應(yīng)的距離；rMi(t)為第i個(gè)散射中心微動(dòng)對(duì)應(yīng)的距離。從式(1)可以看出，散射中心的總體運(yùn)動(dòng)由錐體目標(biāo)平動(dòng)和散射中心微動(dòng)合成，表現(xiàn)在距離上是周期性的正弦運(yùn)動(dòng)疊加一個(gè)平動(dòng)趨勢項(xiàng)。由文獻(xiàn)[3]可知：

rMi(t)=rMi(t+TM)

(2)

式中：TM為目標(biāo)的進(jìn)動(dòng)周期。雷達(dá)向由i個(gè)散射點(diǎn)組成的目標(biāo)發(fā)射電磁波，返回的基帶信號(hào)表示為：

sT(t)sM(t)

(3)

其中：

(4)

(5)

σi(t)=σi(t+TM)

(6)

由式(2)和式(6)可得：

sM(t)=sM(t+TM)

(7)

第i個(gè)散射中心的瞬時(shí)多普勒頻率為：

(8)

(9)

(10)

fMi(t)為目標(biāo)第i個(gè)散射中心的微多普勒頻率，fT(t)為目標(biāo)的平動(dòng)多普勒頻率，可以被建模為多項(xiàng)式函數(shù)[17]如下：

fT(t)=c0+c1t+c2t2+…+cLtK

(11)

2 雷達(dá)微動(dòng)信號(hào)的時(shí)頻分析

由式(7)可知，微動(dòng)信號(hào)有著與進(jìn)動(dòng)相同的周期，目標(biāo)周期性的進(jìn)動(dòng)對(duì)雷達(dá)回波進(jìn)行了周期性的時(shí)變頻率調(diào)制。時(shí)頻分析可以描述信號(hào)的頻率隨著時(shí)間變化的規(guī)律，因此可用來體現(xiàn)進(jìn)動(dòng)目標(biāo)微多普勒效應(yīng)的周期性。短時(shí)傅里葉變換(STFT)是時(shí)變信號(hào)的線性時(shí)頻表示，為簡單起見，這里使用STFT來證明微動(dòng)信號(hào)時(shí)頻分布的周期性：

ρM(t,f)=

(12)

下面分析式(12)的周期性，令t=t+TM，由于sM(t)=sM(t+TM)，可得：

ρm(t+TM,f)=

(13)

令τ=τ1+TM，可得：

ρM(t+TM,f)=

exp(-j2πfTM)ρM(t,f)

(14)

因此，可得：

|ρM(t+TM,f)|=|ρM(t,f)|

(15)

周期進(jìn)動(dòng)引起的微動(dòng)信號(hào)的時(shí)頻分布與進(jìn)動(dòng)具有相同的周期，式(15)的離散形式為：

|ρM(n+NT,m)|=|ρM(n,m)|

(16)

式中：t=nΔt,n=0,1,…，N-1；f=mΔf,m=0,1,…，M-1,Δf=fs/M，Δt是采樣間隔，fs是采樣頻率；NT=[TM/Δt]，[ ]表示將元素四舍五入為最接近的整數(shù)。

當(dāng)存在平動(dòng)時(shí)，由式(3)可知雷達(dá)回波s(t)=sT(t)sM(t)，其STFT為：

(17)

ρ(t,f)的離散形式為ρ(n,m)，則根據(jù)文獻(xiàn)[18]可知：

|ρ(n+NT,m)|=

(18)

3 微動(dòng)周期估計(jì)

CACF和CAMDF都能獨(dú)立進(jìn)行周期估計(jì)[17]，但如果將兩種方法組合，可更為充分地體現(xiàn)出微動(dòng)信號(hào)的周期性，并且會(huì)有更強(qiáng)的穩(wěn)定性，文獻(xiàn)[19]通過合成CACF和CAMDF來估計(jì)彈道中段目標(biāo)的RCS周期。當(dāng)雷達(dá)目標(biāo)存在平動(dòng)時(shí)，雖然雷達(dá)回波信號(hào)的時(shí)頻分布并不表現(xiàn)為正常的正弦調(diào)制形式，但雷達(dá)回波信號(hào)的循環(huán)周期性仍然存在，本文正是利用這一特性，將時(shí)頻分布的CAMDF的倒數(shù)和CACF相乘，來估計(jì)彈道目標(biāo)的進(jìn)動(dòng)周期。

首先，分別定義ρ(n,m)2個(gè)時(shí)間切片n1和n2的循環(huán)自相關(guān)函數(shù)CACF[18]和平均幅度差函數(shù)CAMDF[20]如下：

Cr(k;n1,n2)=

(19)

Ac(k;n1,n2)=

(20)

然后構(gòu)造新的函數(shù)，定義時(shí)間切片n1和n2的CACF/CAMDF如下:

(21)

令n2=n2+NT，可得：

FCC(k;n1,n2+NT)=

(22)

將式(18)代入式(22)得：

FCC(k;n1,n2+NT)=

(23)

n1,n2+NT)=

(24)

由式(21)和式(24)可得:

FCC(k;n1,n2)=

n1,n2+NT)

(25)

由式(25)可以看出，F(xiàn)CC(k;n1,n2)是FCC(k;n1,n2+NT)的循環(huán)移位，并且CACF/CAMDF保持了ρ(n,m)的循環(huán)周期性。

然后，基于CACF/CAMDF定義ρ(n,m)的循環(huán)系數(shù)矩陣為:

MF=

(26)

式中：Cm(n1,n2)是Cr(k;n1,n2)的歸一化最大值，Am(n1,n2)是Ac(k;n1,n2)的歸一化最小值，即：

(27)

由Fm(n1,n2)的定義，可得:

(28)

特別地，當(dāng)n2=n1+aNT(a=0,1,…)時(shí)，F(xiàn)m(n1,n2)=0，所以很容易可以看出Fm(n1,n2)與Fm(k;n1,n2)具有相同的周期，即：

Fm(n1,n2)=Fm(n1,n2+NT)

(29)

因此，ρ(n,m)的循環(huán)系數(shù)矩陣保留了行和列維的周期性，并且在同一對(duì)角線上的滯后是常數(shù)。

最后，將循環(huán)系數(shù)矩陣MF第k個(gè)對(duì)角線的平均值定義為平均循環(huán)系數(shù)，即：

(30)

其中：

diag(MF,k)={MF(n1,n2),

n2=n1+k,1≤n1≤N,1≤n2≤N}

(31)

圖2 算法實(shí)現(xiàn)流程

4 實(shí)驗(yàn)與分析

下面通過仿真驗(yàn)證本文提出的CACF/CAMDF的估計(jì)效果和抗噪性能，并將其與CACF、CAMDF進(jìn)行對(duì)比。首先分別在目標(biāo)的平動(dòng)被補(bǔ)償、目標(biāo)存在平動(dòng)條件下，對(duì)比3種算法的平動(dòng)敏感性；然后加入不同信噪比的噪聲，進(jìn)一步對(duì)比3種算法在低信噪比條件下的估計(jì)性能。

仿真條件：假設(shè)雷達(dá)工作在10 GHz，且雷達(dá)坐標(biāo)系中本體坐標(biāo)系原點(diǎn)O的坐標(biāo)為(400,500,100)km，本地坐標(biāo)系和參考坐標(biāo)系之間的初始?xì)W拉角(x-y-z序列)為(30°,60°,45°)。假設(shè)目標(biāo)繞z軸旋轉(zhuǎn)，目標(biāo)上有2個(gè)散射點(diǎn)：第一散射點(diǎn)A位于錐頂，在本體坐標(biāo)系中的坐標(biāo)為(0,0,1)m；第二個(gè)散射點(diǎn)B位于錐底的尾翼，在本體坐標(biāo)系中的坐標(biāo)是(0.5,0,-0.5)m。雷達(dá)照射時(shí)間為8 s。旋轉(zhuǎn)頻率為fs=1 Hz，圓錐運(yùn)動(dòng)頻率為fc=0.5 Hz。當(dāng)目標(biāo)存在平動(dòng)時(shí)，其平動(dòng)模型為rT(t)=-vt+0.4t2+0.5t3。

4.1 平動(dòng)敏感性分析

當(dāng)不存在噪聲時(shí)，圖3給出了平動(dòng)被補(bǔ)償和存在平動(dòng)兩種情況下的雷達(dá)回波STFT時(shí)頻分布圖。

圖3 無噪聲條件下的時(shí)頻分布

由圖3可知，當(dāng)平動(dòng)被完全補(bǔ)償時(shí)，雷達(dá)回波的時(shí)頻圖是多個(gè)正弦分量合成的類正弦曲線，如圖3(a)所示；而時(shí)頻圖在目標(biāo)平動(dòng)的調(diào)制下，其圖形不再表現(xiàn)為正常的類正弦曲線，如圖3(b)所示。此時(shí)，CACF/CAMDF的估計(jì)結(jié)果如圖4所示，其中圖4(a)和圖4(b)分別為時(shí)頻分布的循環(huán)系數(shù)矩陣和平均循環(huán)系數(shù)。

圖4 平動(dòng)被補(bǔ)償、無噪聲時(shí)的CACF/CAMDF估計(jì)結(jié)果

圖5 存在平動(dòng)、無噪聲時(shí)的CACF/CAMDF估計(jì)結(jié)果

當(dāng)目標(biāo)存在平動(dòng)時(shí)，其他條件相同，圖6分別給出了CAMDF和CACF的估計(jì)結(jié)果。將圖6與圖5對(duì)比可看出，在平動(dòng)條件下，CAMDF和CACF兩種算法都具有一定的估計(jì)能力，且CACF性能更好，但性能都劣于CACF/CAMDF。這說明了CACF/CAMDF具有更低的平動(dòng)敏感性，能夠在目標(biāo)存在平動(dòng)時(shí)，具有更好的估計(jì)性能。

圖6 存在平動(dòng)、無噪聲時(shí)CAMDF、CACF的平均循環(huán)系數(shù)

4.2 抗噪性能分析

為進(jìn)一步驗(yàn)證CACF/CAMDF在噪聲條件下的表現(xiàn)，分別向目標(biāo)信號(hào)中添加不同信噪比的高斯白噪聲進(jìn)行仿真。由圖7可知，即使信噪比低至-5 dB，STFT仍然能夠較好地獲取雷達(dá)回波的時(shí)頻分布。

圖7 噪聲環(huán)境下存在平動(dòng)時(shí)的時(shí)頻分布

由圖5和圖8對(duì)比可知，本文算法在噪聲環(huán)境下具有良好的表現(xiàn)。

圖8 存在平動(dòng)、SNR=5 dB時(shí)的CACF/CAMDF估計(jì)結(jié)果

為進(jìn)一步分析算法的抗噪性能，圖9給出了不同信噪比條件下、目標(biāo)存在平動(dòng)時(shí)的CACF/CAMDF平均循環(huán)系數(shù)圖。分析圖9可知，該算法具備良好的抗噪性能。即使目標(biāo)存在平動(dòng)、信噪比低至-2 dB，仍可準(zhǔn)確估計(jì)出目標(biāo)的進(jìn)動(dòng)周期。

圖9 存在平動(dòng)時(shí)CACF/CAMDF的平均循環(huán)系數(shù)

目標(biāo)存在平動(dòng)時(shí)，圖10和圖11分別給出了不同信噪比條件下CAMDF和CACF的估計(jì)結(jié)果。對(duì)比圖9～圖11可知，噪聲環(huán)境下，CAMDF的性能極限是2 dB，CACF的性能極限是0 dB，與CACF/CAMDF的-2 dB都具有較大差距。這進(jìn)一步體現(xiàn)了CACF/CAMDF在低噪聲環(huán)境中的性能優(yōu)越性。

圖10 存在平動(dòng)時(shí)CAMDF的平均循環(huán)系數(shù)

圖11 存在平動(dòng)時(shí)CACF的平均循環(huán)系數(shù)

5 結(jié)語

本文分析了微動(dòng)和平動(dòng)對(duì)雷達(dá)回波時(shí)頻分布的影響，結(jié)合CAMDF和CACF，在時(shí)頻分布的基礎(chǔ)上提出了基于CACF/CAMDF的估計(jì)算法。仿真實(shí)驗(yàn)表明，該算法在目標(biāo)存在平動(dòng)的情況下，無需補(bǔ)償平動(dòng)分量可以直接估計(jì)微動(dòng)周期，且估計(jì)性能優(yōu)于CAMDF和CACF，這在信噪比較低時(shí)體現(xiàn)得更加明顯。