基于循環自相關/平均幅度差函數的彈道目標微動周期估計

金家偉， 阮懷林

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

在彈道中段的飛行過程中，為了保持飛行的穩定性，彈頭一般會進行姿態控制，不僅會繞自身對稱軸作自旋運動，還會受到沖擊力矩的作用，力矩消失后對稱軸將在平衡位置做圓錐運動，即進動[1]。而誘餌一般不具備姿態控制裝置，因此，進動是中段彈頭特有的微動方式，與目標形狀結構、質量分布以及動力學特性等物理屬性密切相關[2]。由微動引起的雷達回波中多普勒頻率的調制，一般被稱為微多普勒效應[3]。微多普勒特征體現出目標獨特的動態和結構特征，在目標識別和分類方面具有巨大的潛力[4-5]。

對于非剛性雷達目標，目標主體散射的雷達回波占返回能量的絕大部分，微動信號始終浸沒在主體雷達回波中，需要將雷達回波分解為平動分量和微動分量，用微動分量來估計微多普勒參數[6-8]。文獻[9]對信號進行小波變換，用分解出的微動分量的自相關來估計振動/旋轉速率；文獻[10]通過信號分解提取單個散射點特征。而為了表示出散射點雷達回波的時變頻率調制特性，時頻分析被廣泛用于分析微多普勒特征，基于時頻分布的微動參數估計已成為微多普勒特征提取的重要部分。文獻[11]采用基于加權迭代自適應的時頻分析方法結合逆Radon變換，分離重構不同散射點的微多普勒分量；文獻[12]利用短時分數階傅里葉變換分離肢體和軀干的微多普勒信號；文獻[13]結合Gabor時頻分布和變分模態分解，來估計目標的自旋頻率和錐旋頻率。

上述方法都需要假設目標平動已被準確補償，在此前提下，目標時頻分布中的微多普勒頻率曲線通常應該表現為正弦曲線[3]或是多個正弦分量的合成[2]。現有的微多普勒特征提取技術，如經驗模式分解(empirical-mode decomposition, EMD)[14]和擴展Hough變換(extended hough transform,EHT)[15]，主要依賴于這種正弦調制假設。但目標運動由平動與微動復合而成，特別是對于高速運動的彈道目標，周期進動引起的彈道目標微多普勒效應并不像通常那樣表現為正常的正弦調制[16]。

針對目標微動伴隨平動時，時頻分布不再表現為正弦調制曲線，本文在雷達回波信號時頻分布的基礎上，結合循環自相關函數(circular autocorrelation function，CACF)和循環平均幅度差函數(circular average amplitude difference function，CAMDF)，以獲取時頻分布的循環系數矩陣，并通過該矩陣的平均循環系數估計出彈道目標的微動周期。

1 模型建立

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

圖1 目標模型

在光學區，雷達目標的整體散射特性通常可以等效為若干個散射中心的疊加。為不失一般性，我們假設目標的微動是周期性進動，目標的平動是任意且未知的。假設目標等效為K個散射中心，各散射中心各向同性，雷達發射的電磁波為連續單頻波，載頻為f0。考慮目標的平動，第i個散射中心與雷達之間的距離為：

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

(1)

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

rMi(t)=rMi(t+TM)

(2)

式中：TM為目標的進動周期。雷達向由i個散射點組成的目標發射電磁波，返回的基帶信號表示為：

sT(t)sM(t)

(3)

其中：

(4)

(5)

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

(6)

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

sM(t)=sM(t+TM)

(7)

第i個散射中心的瞬時多普勒頻率為：

(8)

(9)

(10)

fMi(t)為目標第i個散射中心的微多普勒頻率，fT(t)為目標的平動多普勒頻率，可以被建模為多項式函數[17]如下：

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

(11)

2 雷達微動信號的時頻分析

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

ρ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)

周期進動引起的微動信號的時頻分布與進動具有相同的周期，式(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]，[ ]表示將元素四舍五入為最接近的整數。

當存在平動時，由式(3)可知雷達回波s(t)=sT(t)sM(t)，其STFT為：

(17)

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

|ρ(n+NT,m)|=

(18)

3 微動周期估計

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

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

Cr(k;n1,n2)=

(19)

Ac(k;n1,n2)=

(20)

然后構造新的函數，定義時間切片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)可以看出，FCC(k;n1,n2)是FCC(k;n1,n2+NT)的循環移位，并且CACF/CAMDF保持了ρ(n,m)的循環周期性。

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

MF=

(26)

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

(27)

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

(28)

特別地，當n2=n1+aNT(a=0,1,…)時，Fm(n1,n2)=0，所以很容易可以看出Fm(n1,n2)與Fm(k;n1,n2)具有相同的周期，即：

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

(29)

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

最后，將循環系數矩陣MF第k個對角線的平均值定義為平均循環系數，即：

(30)

其中：

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

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

(31)

圖2 算法實現流程

4 實驗與分析

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

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

4.1 平動敏感性分析

當不存在噪聲時，圖3給出了平動被補償和存在平動兩種情況下的雷達回波STFT時頻分布圖。

圖3 無噪聲條件下的時頻分布

由圖3可知，當平動被完全補償時，雷達回波的時頻圖是多個正弦分量合成的類正弦曲線，如圖3(a)所示；而時頻圖在目標平動的調制下，其圖形不再表現為正常的類正弦曲線，如圖3(b)所示。此時，CACF/CAMDF的估計結果如圖4所示，其中圖4(a)和圖4(b)分別為時頻分布的循環系數矩陣和平均循環系數。

圖4 平動被補償、無噪聲時的CACF/CAMDF估計結果

圖5 存在平動、無噪聲時的CACF/CAMDF估計結果

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

圖6 存在平動、無噪聲時CAMDF、CACF的平均循環系數

4.2 抗噪性能分析

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

圖7 噪聲環境下存在平動時的時頻分布

由圖5和圖8對比可知，本文算法在噪聲環境下具有良好的表現。

圖8 存在平動、SNR=5 dB時的CACF/CAMDF估計結果

為進一步分析算法的抗噪性能，圖9給出了不同信噪比條件下、目標存在平動時的CACF/CAMDF平均循環系數圖。分析圖9可知，該算法具備良好的抗噪性能。即使目標存在平動、信噪比低至-2 dB，仍可準確估計出目標的進動周期。

圖9 存在平動時CACF/CAMDF的平均循環系數

目標存在平動時，圖10和圖11分別給出了不同信噪比條件下CAMDF和CACF的估計結果。對比圖9～圖11可知，噪聲環境下，CAMDF的性能極限是2 dB，CACF的性能極限是0 dB，與CACF/CAMDF的-2 dB都具有較大差距。這進一步體現了CACF/CAMDF在低噪聲環境中的性能優越性。

圖10 存在平動時CAMDF的平均循環系數

圖11 存在平動時CACF的平均循環系數

5 結語

本文分析了微動和平動對雷達回波時頻分布的影響，結合CAMDF和CACF，在時頻分布的基礎上提出了基于CACF/CAMDF的估計算法。仿真實驗表明，該算法在目標存在平動的情況下，無需補償平動分量可以直接估計微動周期，且估計性能優于CAMDF和CACF，這在信噪比較低時體現得更加明顯。