侯煜博,沈曉峰,李思思
(電子科技大學(xué)電子工程學(xué)院,成都 611731)
彈道導(dǎo)彈(tactical ballistic missile,TBM)作為現(xiàn)代戰(zhàn)爭中威力極強的進攻性武器,具有射程遠、精度高、速度快、殺傷力大、突防能力強等特點,給導(dǎo)彈防御提出了嚴峻的挑戰(zhàn)。穩(wěn)定的跟蹤和精確的制導(dǎo)是導(dǎo)彈攔截成功的關(guān)鍵,是導(dǎo)彈防御雷達的核心任務(wù)。
由于導(dǎo)彈末段高速機動的特性,傳統(tǒng)的勻速(CV)、勻加速(CA)模型在濾波過程中會產(chǎn)生較大誤差。文中針對再入大氣層的來襲導(dǎo)彈目標,通過分析其運動規(guī)律,從解析導(dǎo)彈運動方程入手,在東北天(East-North-Up,ENU)坐標系下建立目標狀態(tài)方程,應(yīng)用擴展卡爾曼(extended Kalman filter,EKF)濾波,對目標進行跟蹤。
在給定來襲目標運動規(guī)律、攔截彈飛行速度規(guī)律和具體導(dǎo)引方法的條件下,可以確定出將兩者視作質(zhì)點時導(dǎo)彈飛向目標的理想運動學(xué)彈道。許多文獻僅在二維情況下或三維情況下目標做勻速直線運動時給出推導(dǎo)[1-3],文中將針對導(dǎo)彈的標準彈道軌跡給出三維制導(dǎo)圖。
TBM再入大氣層后,主要受到重力和空氣阻力作用,阻力影響通常超過重力,因此可以采用簡潔的平方反比重力加速度模型代替橢球地球重力加速度模型,重力產(chǎn)生的加速度為:
空氣阻力引起的加速度為:
標準大氣參數(shù)公式可以參照文獻[5]。為給后續(xù)雅克比矩陣求解帶來方便,文中利用指數(shù)函數(shù)來近似擬合。圖1為在ρ(h)= ρ0e-κh條件下的擬合結(jié)果與標準結(jié)果的對比圖,其中ρ0和κ均為已知常數(shù)。由圖可以看出擬合效果較為準確,最大偏差為0.0269kg/m3。
圖1 大氣密度函數(shù)精確值與擬合值的對比
仿真生成彈道軌跡是進行跟蹤與制導(dǎo)的基礎(chǔ),根據(jù)受力模型,利用龍哥庫塔積分目標動力學(xué)微分方程組的方式可以得到再入段導(dǎo)彈的彈道軌跡。
設(shè)計一條再入角度為39°的標準彈道,圖2為彈道軌跡仿真圖,圖3為彈道參數(shù)。
圖2 導(dǎo)彈再入段的彈道軌跡
擴展卡爾曼濾波是最通用的非線性濾波方法,在確定了目標的狀態(tài)方程與量測方程及其對應(yīng)的雅克比矩陣后,即可應(yīng)用EKF算法進行濾波。EKF算法流程可以參考文獻[6]。根據(jù)受力模型推導(dǎo)出目標的狀態(tài)方程和量測方程及其對應(yīng)的雅克比矩陣是EKF濾波算法的關(guān)鍵。
圖3 彈道參數(shù)
當TBM再入大氣層后(距海平面約91km以內(nèi)),由式(1)、式(2)可知其加速度表達式為:
在采樣間隔T下目標離散化的狀態(tài)方程[7]為:
狀態(tài)噪聲w(k)為零均值白噪聲序列,其方差為Q(k)。
由于量測信息來自雷達站球坐標系,所以量測方程是非線性的,其表達式為:
量測噪聲ω(k)的方差矩陣R(k)由雷達具體的距離、方位角、俯仰角測量誤差確定。
量測方程的雅克比矩陣為:
在確定了上述狀態(tài)方程與量測方程及其相應(yīng)的雅可比矩陣后,即可應(yīng)用EKF進行濾波。
假定在上述濾波過程中已獲得來襲目標的運動規(guī)律且攔截彈的速度變化規(guī)律已知,如何將導(dǎo)彈引向目標便是制導(dǎo)規(guī)律所要完成的任務(wù)。比例導(dǎo)引法是一種最常用的、制導(dǎo)精度較高的制導(dǎo)規(guī)律。追蹤法、前置角法、平行接近法都可以看作是比例導(dǎo)引法的特殊情況。
比例導(dǎo)引法是指攔截彈飛行過程中速度矢量Vm的轉(zhuǎn)動角速度dθm/dt與導(dǎo)彈、目標連線的旋轉(zhuǎn)角速度dq/dt成比例的導(dǎo)引方法[9]。圖4與圖5分別為導(dǎo)彈與目標相對位置和比例導(dǎo)引示意圖。
圖4 導(dǎo)彈、目標相對位置
圖5 比例導(dǎo)引示意圖
圖中 Mk-1、Mk、Tk-1、Tk分別為 k - 1時刻、k時刻導(dǎo)彈位置和k-1時刻、k時刻目標位置。A為Mk-1Mk與 Tk-1Tk交點,MkB 平行于 Mk-1Tk-1。
設(shè) xt、yt、zt表示目標坐標,任意時刻的 xt、yt、zt已知。算法流程如下:
k-1時刻導(dǎo)彈與目標的距離:
k-1時刻目標速度矢量與導(dǎo)彈、目標連線的夾角:
k時刻導(dǎo)彈、目標連線與基準線夾角增量:為k-1時刻導(dǎo)彈與k時刻目標之間距離,st=vt(k-1)·Δt,即在Δt內(nèi)將目標視為做勻速直線運動。
比例導(dǎo)引法的差分方程為:
利用余弦定理求得:
設(shè)xA、yA、zA表示A點坐標,xm、ym、zm表示導(dǎo)彈坐標,由幾何關(guān)系知:
其中sm=vm(k-1)·Δt,即在Δt內(nèi)將導(dǎo)彈視為做勻速直線運動。在獲得k時刻的導(dǎo)彈坐標后即可更新r(k)進行下一次循環(huán)。必要時還要利用c3對dq進行校正。
假定雷達在跟蹤再入段時的量測標準差為:σR=10m,σA= σE=1mrad。跟蹤數(shù)據(jù)率設(shè)定為100Hz。采用兩點法初始化狀態(tài)變量,初始協(xié)方差可以參考文獻[10]。在濾波過程中,狀態(tài)噪聲為設(shè)計參數(shù),文中采用分段噪聲,先利用較大噪聲使算法快速收斂,再采用較小噪聲降低濾波誤差。對目標的跟蹤精度采用均方根誤差(RMSE)作為判別標準[11],設(shè)X^n( )k|k為k時刻第n次蒙特卡羅的狀態(tài)估計值,則:
對量測數(shù)據(jù)進行一次濾波,得到真實軌跡、量測點跡與濾波軌跡的結(jié)果見圖6。蒙特卡羅仿真500次得到均方根誤差結(jié)果如圖7。
圖6 導(dǎo)彈再入段的真實、量測、濾波軌跡
為方便對比,圖7中給出了針對上述相同軌跡應(yīng)用CV模型的濾波結(jié)果。由圖可以看出算法濾波收斂速度較快,誤差較小且收斂后基本穩(wěn)定。較CV模型濾波精度有較大提升。合理改變跟蹤數(shù)據(jù)率或雷達量測誤差,在不同的彈道模型下進行多次仿真比較,結(jié)果類似。
圖7 再入段濾波的位置RMSE誤差
設(shè)定比例系數(shù)為K=5,導(dǎo)彈初始位置為坐標原點,速度恒定為 3000m/s,時間間隔 Δt=0.01s,導(dǎo)彈與目標之間相對距離小于60m時認為制導(dǎo)成功。獲得三維理想攔截彈彈道仿真圖見圖8,圖9為攔截彈的法向過載。
圖8 理想三維彈道仿真
圖9 攔截彈法向過載
文中針對再入大氣層的來襲導(dǎo)彈目標受力情況推導(dǎo)了ENU坐標系下目標運動加速度表達式,建立了混合坐標系下的目標狀態(tài)方程與量測方程,通過受力分析采用基于其特性的彈道模型進行濾波,給出了狀態(tài)估計與協(xié)方差矩陣的初始值計算方法,采用EKF濾波算法對目標進行跟蹤,為制導(dǎo)過程提供了來襲導(dǎo)彈的速度位置信息。利用比例導(dǎo)引制導(dǎo)律進行仿真得到三維制導(dǎo)圖。仿真結(jié)果表明,算法可以有效的對再入類導(dǎo)彈進行穩(wěn)定的跟蹤和精確的制導(dǎo)。
[1]高尚.比例導(dǎo)引理想彈道仿真[J].計算機工程與設(shè)計,2003,24(8):66-68.
[2]陳宇強,趙育善.導(dǎo)彈最優(yōu)制導(dǎo)律彈道仿真分析[J].指揮控制與仿真,2007,29(3):92-105.
[3]周紀元,童幼堂,張磊,等.典型導(dǎo)引規(guī)律三維彈道仿真分析[J].艦船電子工程,2008,28(2):110-112.
[4]張毅,肖龍旭,王順宏.彈道導(dǎo)彈彈道學(xué)[M].長沙:國防科技大學(xué)出版社,2005.
[5]楊炳尉.標準大氣參數(shù)的公式表示[J].宇航學(xué)報,1983(1):83-86.
[6]何子述,夏威.現(xiàn)代數(shù)字信號處理及其應(yīng)用[M].北京:清華大學(xué)出版社,2009.
[7]易令,呂明.使用交互多模型算法的高速高機動目標跟蹤[J].雷達科學(xué)與技術(shù),2006,4(3):143-147.
[8]辛召強,沈曉峰.利用CUDA快速實現(xiàn)IMM目標跟蹤[J].雷達科學(xué)與技術(shù),2012,10(6):656-659.
[9]梁慜.彈道導(dǎo)彈攔截仿真建模技術(shù)研究[J].中國電子科學(xué)研究院學(xué)報,2013,8(1):56-59.
[10]何友,修建娟,張晶煒.雷達數(shù)據(jù)處理及應(yīng)用[M].北京:電子工業(yè)出版社,2006.
[11]張丕旭,石章松,劉忠.一種基于彈道模型的機動目標跟蹤算法[J].彈箭與制導(dǎo)學(xué)報,2009,29(4):55-57.