許旭光，馮存前,2，蘇于童，張文強

(1.空軍工程大學防空反導學院，陜西 西安 710051；2.信息感知技術協(xié)同創(chuàng)新中心，陜西 西安 710077；3.西安交通大學數(shù)學與統(tǒng)計學院，陜西 西安 710049)

0 引言

彈道導彈作為現(xiàn)代高科技戰(zhàn)爭條件下的主要攻擊武器之一，具有強大的戰(zhàn)略威懾能力和戰(zhàn)役戰(zhàn)術打擊能力[1]。與此同時，對彈道導彈的識別與分類技術也得到迅速發(fā)展。當前，可供彈道導彈識別的特征包括微動特征、極化特征、一維距離像以及二維ISAR像等[2-3]。其中，微動作為彈道目標在中段特有的運動方式，能夠精確描述彈道目標物理特征和運動特征，對于彈道目標的特征提取和分類有重要的意義。然而，除了進行各種形式的微動之外，彈道目標在中段本身還具有高速平動。高速平動會導致微多普勒發(fā)生模糊或折疊，影響微動特征的提取[4]。

目前對于中段彈道目標平動補償方面的研究已經(jīng)取得一定成果。文獻[5]針對高速機動目標提出一種基于EMD(經(jīng)驗模態(tài)分解)的平動補償算法，該算法利用平動多普勒和微多普勒之間的頻率差異實現(xiàn)平動分量和微動分量的分離，而EMD算法存在端點效應、包絡插值誤差等缺點，會影響平動參數(shù)的估計。文獻[6]針對寬帶雷達回波，采用簡化分數(shù)階傅里葉變換實現(xiàn)對目標平動參數(shù)的估計，該算法將平動等效為勻速運動，所建立平動模型過于簡單，補償精度不夠高。文獻[7]針對多普勒曲線極值點中僅包含平動信息這一特點，提取多普勒曲線極值點并進行最小二乘估計，實現(xiàn)對平動參數(shù)的估計，該算法僅在單散射點能夠?qū)崿F(xiàn)較好的補償，在多散射點條件下則存在算法復雜，估計精度不足的問題。文獻[8]針對平動對進動目標微多普勒帶來的影響，采用基于最小熵值的平動補償算法對加速度進行補償，進而消除加速度對微多普勒的影響。該算法僅考慮加速度對時頻曲線的影響，對于速度對微多普勒的影響考慮不足。

針對上述平動補償方法存在的補償精度不足、算法復雜的問題，本文提出一種基于高階模糊函數(shù)的彈道目標平動參數(shù)估計方法。

1 復合運動模型

1.1 微動回波模型

建立錐體彈道目標進動模型如圖1所示。設進動軸為z′軸，將平底錐形錐彈道目標彈體對稱軸與進動軸交點O′作為坐標原點，建立參考坐標系O′x′y′z′，雷達視線LOS在參考坐標系中的方位角和俯仰角分別為α和β，ψ為雷達視線與彈體對稱軸夾角，ξ為彈體對稱軸與z′軸夾角，A為錐頂散射中心，B，C為兩個滑動散射中心，O′′為圓錐彈頭底面圓心，設底面圓半徑為r，O′A長度為h1，O′O′′長度為h2?？紤]到平底錐型彈道目標為旋轉(zhuǎn)對稱目標，自旋對目標微距離影響不會在回波中體現(xiàn)出來，因此，僅考慮錐旋對目標的影響，設錐旋頻率為ωc。

圖1 進動目標滑動散射模型Fig.1 The sliding scattering model of precession target

根據(jù)文獻[9]可以得出各散射點的微距離表示為：

(1)

式(1)中，

(2)

若雷達發(fā)射單頻信號，微動分量對應的基頻回波為：

(3)

式(3)中，f0為雷達發(fā)射電磁波頻率。

1.2 復合運動回波模型

當考慮高速平動對目標回波的影響時，目標的距離函數(shù)可以修正為：

r(t)=rm(t)+rtr(t)

(4)

式(4)中，平動引起的距離變化可以表示為

(5)

式(5)中，r0為初始時刻目標到雷達的距離，P為平動引起的距離變化的等效多項式的階數(shù)。

若考慮平動，目標回波由式(3)修正為：

(6)

2 基于高階模糊函數(shù)的平動補償原理

2.1 高階模糊函數(shù)理論

加州大學的Shimon Peleg提出高階模糊函數(shù)(High-order Ambiguity Function, HAF),并將其應用到多項式相位信號(Polynomial Phase Signal,PPS)參數(shù)估計中[10]。與分數(shù)階傅里葉變換、匹配傅里葉變換等信號參數(shù)估計算法相比，由于其可估計階數(shù)高、運算簡單且運算量相對較小，因而得到廣泛應用于PPS參數(shù)估計。下面介紹基于高階模糊函數(shù)的多項式相位信號估計原理：

給定信號x(t),定義其高階矩函數(shù)為：

(7)

將式(7)中信號高階矩函數(shù)進行拓展，得到信號高階矩函數(shù)的統(tǒng)一表達式為：

(8)

式(8)中，

(9)

PM(t,τ)=A2M-1ej(2πfct+φ)

(10)

式(10)中，

(11)

式(11)中，!表示階乘。

對式(10)作傅里葉變換，得到高階模糊函數(shù)(High-order Ambiguity Function,HAF)為：

(12)

對于式(10)中形成的單頻信號,通過式(12)將信號展現(xiàn)在頻域，通過對信號頻域峰值的搜索,利用式(11)可實現(xiàn)對PPS相位參數(shù)的估計。

2.2 平動參數(shù)估計

1.2節(jié)中指出中段目標回波包括平動項和微動項，其中平動項為多項式相位信號。因此，本文擬將高階模糊函數(shù)引入平動參數(shù)的估計中。當觀測時間為幾個周期時，三階多項式可以較好地描述目標平動[11]。根據(jù)式(7)，回波的三階矩可以表示為

(13)

式(13)中，

(14)

式(14)中，F(xiàn)3(t,τ)為由平動分量相乘得到，為一單頻信號，其相位?3對于參數(shù)估計沒有影響，因此不詳細寫出其表達式；G3(t,τ)表示回波三階矩中同一散射中心微動部分延遲共軛相乘結(jié)果之和，為微動自項，當延遲時間為微動周期時，該項為常數(shù)；H3(t,τ)表示回波三階矩中不同散射中心微動信號的耦合，為微動互項，當延遲時間為微動周期時，該項在頻域中能量相對較小，不會影響平動參數(shù)的估計。綜上分析可得，當延遲時間為微動周期時，其對應階數(shù)的模糊函數(shù)相應的頻率處會產(chǎn)生峰值。

對式(13)進行傅里葉變換,得到三階模糊函數(shù)X3(f,τ),進行二維峰值搜索：

(15)

當對微動信號高階模糊函數(shù)進行峰值搜索，通過峰值對應的頻率和延遲時間的位置，既可以實現(xiàn)對目標微動周期的估計，又可以實現(xiàn)對相應階數(shù)平動項系數(shù)的估計。

(16)

式(16)中，

(17)

(18)

(19)

(20)

通過式(20)可以實現(xiàn)對速度的估計，進而實現(xiàn)對平動項的補償。

通過上述算法分析，可以看出本文所提算法計算簡單，對具有周期性的微動形式均可以實現(xiàn)有效的平動補償。由于對參數(shù)的估計采用的是峰值搜索的思想，因此估計性能也相對穩(wěn)定。

綜合以上分析，總結(jié)得到彈道目標平動補償流程為：

步驟3 對補償后的回波作半周期延遲相乘處理，對得到信號作傅里葉變換，對a1進行估計，并對回波進行補償。

3 仿真分析

3.1 彈道仿真

本文采用Matlab 2014a仿真平臺進行仿真實驗。假設導彈關機點高度250 km，關機點速度4.5 km/s，按照最佳速度傾角飛行[14]，地球半徑6 337 km，開普勒常數(shù)3.98×1014，得到彈道仿真如圖2所示。

針對平動仿真的研究，多將平動等效為二階多項式或三階多項式。為進一步分析兩種不同階數(shù)多項式等效效果，取導彈飛行時間為0～6 s，將距離參數(shù)分別進行二階多項式擬合和三階多項式擬合，具體擬合方法為采用Matlab中的“polyfit”函數(shù)得到進行擬合，擬合結(jié)果如圖3所示。

圖2 彈道仿真Fig.2 The trajectory simulation

圖3 不同階數(shù)多項式擬合效果Fig.3 Performance of different order polynomial fitting

對比圖3(a)和圖3(b)可以發(fā)現(xiàn)，二階多項式的擬合誤差在10-2m，三階多項式的擬合誤差在10-6m，說明采用三階多項式對平動進行等效的精度相對更高。

3.2 平動補償仿真

假設雷達發(fā)射單載脈沖信號，載頻為f0=10 GHz，脈沖重復頻率為fr=2 000 Hz，觀測時間為T=5 s，目標為圓錐彈頭，h1=1 m，h2=0.6 m，底面半徑r=0.8 m,半錐角γ=26.6°，雷達視線的方位角和高低角(α,β)=(20°,70°)；彈頭錐旋角ξ=8°，進動頻率ωc=5π rad/s。不同散射點的散射系數(shù)不同，假設A，B，C散射點散射系數(shù)之比為1.5∶1∶1。設經(jīng)粗補償后平動速度V=5 m/s，一階加速度a1=2m/s2,二階加速度a2=0.5m/s3?？紤]到在觀測時間內(nèi)，目標存在姿態(tài)上的起伏，這會導致單個脈沖的信噪比難以被量化。因此，在仿真過程中設定觀測時間內(nèi)的平均信噪比為3 dB。后文關于信噪比的敘述與此處意義相同。采用短時傅里葉變換時頻分析方法對回波進行時頻分析,窗函數(shù)選擇為hamming窗，對應表達式為

(21)

式(21)中，窗長N=55。

在進行微動周期搜索時，設置搜索步長為0.2 s。

圖4 平動未補償時回波分析Fig.4 The analysis of echo withour translatioinal compensation

圖4(a)為未補償時回波波形，圖4(b)為平動未補償?shù)幕夭〞r頻圖。由圖4(b)可以看出多普勒曲線出現(xiàn)嚴重變形，無法通過現(xiàn)有曲線實現(xiàn)微動參數(shù)估計。

圖5(a)為對信號的三階模糊函數(shù)進行延遲時間搜索。圖中出現(xiàn)一系列峰值，通過采用式(15)對峰值進行處理，實現(xiàn)對二階加速度的估計，進而實現(xiàn)對二階加速度的補償，圖5(b)為對二階加速度補償后的結(jié)果。同理，圖5(c)和圖5(d)為加速度補償結(jié)果。圖5(e)和圖5(f)分別為采用本文的半周期延遲相乘和文獻[13]中所提算法對速度進行估計，可見本文所提算法能夠頻域中有效抑制微動分量的影響，圖5(g)為最終補償后得到的時頻圖。

表1 參數(shù)估計性能分析

表1中反映了本文方法對平動參數(shù)及微動周期的估計效果。由表中可知，各參數(shù)的估計誤差均保持在5%以內(nèi)，因此可以認為文中所提方法具有較高的平動參數(shù)估計精度。需要說明的是由于本文設定微動周期搜索步長為0.2 s，而微動周期為0.4 s，恰好可以實現(xiàn)微動周期的0誤差估計。實際過程中無法實現(xiàn)如此高精度的估計，但是通過合理調(diào)整搜索步長，可以提高微動周期估計精度。由于文中采取的基于高階模糊函數(shù)的參數(shù)估計實質(zhì)是對回波進行延遲共軛相乘處理，而處理之后得到的回波的信噪比一定會下降[15]，因此當回波信噪比降低到某臨界值時，峰值搜索過程中出現(xiàn)多個虛假峰值，導致a1、a2、a3無法被正確估計。仿真發(fā)現(xiàn)該臨界信噪比值為-4 dB。

圖5 平動參數(shù)估計及補償Fig.5 Parameter estimation and compensation of translational motion

分別采用文獻[5]中提出的平動補償算法和本文提出的平動補償算法對復合運動信號進行處理，對補償后的微多普勒作誤差分析如圖6所示。

圖6 平動補償后微多普勒誤差分析Fig.6 The error anslysis of micro-Doppler withranslational motion compensated

從圖6可以看出，相對于文獻[5]所提的算法，經(jīng)過本文所提算法補償后的微多普勒具有相對誤差低，補償效果好。

文獻[7]中基于最強散射點的平動參數(shù)估計僅適用于微距離為正弦形式的情況。而本文提出的方法對于具有周期性的微動形式均可以實現(xiàn)有效的平動補償，具有較強的實際應用性，并且本文提出的算法在對平動參數(shù)進行估計的同時，還可以實現(xiàn)微動周期的有效估計。文獻[8]提出的基于最小時頻熵的平動補償方法雖然操作簡單且具有一定的補償能力，該算法僅能實現(xiàn)去除徑向加速度對微多普勒分析帶來的影響，對于速度對微多普勒的影響則依賴于雷達自身的測速精度，具有一定的局限性。與文獻[12]中的平動補償方法相比，平動參數(shù)的估計誤差基本為一個量級。仿真發(fā)現(xiàn)在信噪當大于-4 dB時，各次估計結(jié)果完全相同，因此估計均方誤差明顯小于文獻[12]中所提方法，所以估計效果相對穩(wěn)定。

4 結(jié)論

本文提出了基于高階模糊函數(shù)的進動目標平動補償算法。該算法將平動引起的距離變化等效為三階多項式，采用基于高階模糊函數(shù)的思想實現(xiàn)了對微動周期和平動參數(shù)的估計，從而實現(xiàn)對中段彈道目標平動分量的有效補償。仿真實驗證明，該算法受噪聲影響相對較小，參數(shù)估計精度高，且適用于任何具有周期性的微動形式?？紤]到本文對于回波高階矩的求解需要相對較長的觀測時間，后續(xù)將研究如何在較短的觀測時間內(nèi)實現(xiàn)更為精確的平動補償。