彈道導(dǎo)彈再入段跟蹤制導(dǎo)算法研究*

侯煜博，沈曉峰，李思思

（電子科技大學(xué)電子工程學(xué)院，成都 611731）

0 引言

彈道導(dǎo)彈（tactical ballistic missile，TBM）作為現(xiàn)代戰(zhàn)爭(zhēng)中威力極強(qiáng)的進(jìn)攻性武器，具有射程遠(yuǎn)、精度高、速度快、殺傷力大、突防能力強(qiáng)等特點(diǎn)，給導(dǎo)彈防御提出了嚴(yán)峻的挑戰(zhàn)。穩(wěn)定的跟蹤和精確的制導(dǎo)是導(dǎo)彈攔截成功的關(guān)鍵，是導(dǎo)彈防御雷達(dá)的核心任務(wù)。

由于導(dǎo)彈末段高速機(jī)動(dòng)的特性，傳統(tǒng)的勻速（CV）、勻加速（CA）模型在濾波過(guò)程中會(huì)產(chǎn)生較大誤差。文中針對(duì)再入大氣層的來(lái)襲導(dǎo)彈目標(biāo)，通過(guò)分析其運(yùn)動(dòng)規(guī)律，從解析導(dǎo)彈運(yùn)動(dòng)方程入手，在東北天（East-North-Up，ENU）坐標(biāo)系下建立目標(biāo)狀態(tài)方程，應(yīng)用擴(kuò)展卡爾曼（extended Kalman filter，EKF）濾波，對(duì)目標(biāo)進(jìn)行跟蹤。

在給定來(lái)襲目標(biāo)運(yùn)動(dòng)規(guī)律、攔截彈飛行速度規(guī)律和具體導(dǎo)引方法的條件下，可以確定出將兩者視作質(zhì)點(diǎn)時(shí)導(dǎo)彈飛向目標(biāo)的理想運(yùn)動(dòng)學(xué)彈道。許多文獻(xiàn)僅在二維情況下或三維情況下目標(biāo)做勻速直線運(yùn)動(dòng)時(shí)給出推導(dǎo)［1－3］，文中將針對(duì)導(dǎo)彈的標(biāo)準(zhǔn)彈道軌跡給出三維制導(dǎo)圖。

1 彈道軌跡生成

1.1 目標(biāo)受力模型

TBM再入大氣層后，主要受到重力和空氣阻力作用，阻力影響通常超過(guò)重力，因此可以采用簡(jiǎn)潔的平方反比重力加速度模型代替橢球地球重力加速度模型，重力產(chǎn)生的加速度為:

空氣阻力引起的加速度為:

標(biāo)準(zhǔn)大氣參數(shù)公式可以參照文獻(xiàn)［5］。為給后續(xù)雅克比矩陣求解帶來(lái)方便，文中利用指數(shù)函數(shù)來(lái)近似擬合。圖1為在ρ（h）= ρ0e－κh條件下的擬合結(jié)果與標(biāo)準(zhǔn)結(jié)果的對(duì)比圖，其中ρ0和κ均為已知常數(shù)。由圖可以看出擬合效果較為準(zhǔn)確，最大偏差為0.0269kg/m3。

圖1 大氣密度函數(shù)精確值與擬合值的對(duì)比

1.2 彈道生成

仿真生成彈道軌跡是進(jìn)行跟蹤與制導(dǎo)的基礎(chǔ)，根據(jù)受力模型，利用龍哥庫(kù)塔積分目標(biāo)動(dòng)力學(xué)微分方程組的方式可以得到再入段導(dǎo)彈的彈道軌跡。

設(shè)計(jì)一條再入角度為39°的標(biāo)準(zhǔn)彈道，圖2為彈道軌跡仿真圖，圖3為彈道參數(shù)。

圖2 導(dǎo)彈再入段的彈道軌跡

2 跟蹤算法

擴(kuò)展卡爾曼濾波是最通用的非線性濾波方法，在確定了目標(biāo)的狀態(tài)方程與量測(cè)方程及其對(duì)應(yīng)的雅克比矩陣后，即可應(yīng)用EKF算法進(jìn)行濾波。EKF算法流程可以參考文獻(xiàn)［6］。根據(jù)受力模型推導(dǎo)出目標(biāo)的狀態(tài)方程和量測(cè)方程及其對(duì)應(yīng)的雅克比矩陣是EKF濾波算法的關(guān)鍵。

圖3 彈道參數(shù)

2.1 目標(biāo)狀態(tài)方程

當(dāng)TBM再入大氣層后（距海平面約91km以內(nèi)），由式（1）、式（2）可知其加速度表達(dá)式為:

在采樣間隔T下目標(biāo)離散化的狀態(tài)方程［7］為:

狀態(tài)噪聲w（k）為零均值白噪聲序列，其方差為Q（k）。

2.2 目標(biāo)量測(cè)方程

由于量測(cè)信息來(lái)自雷達(dá)站球坐標(biāo)系，所以量測(cè)方程是非線性的，其表達(dá)式為:

量測(cè)噪聲ω（k）的方差矩陣R（k）由雷達(dá)具體的距離、方位角、俯仰角測(cè)量誤差確定。

量測(cè)方程的雅克比矩陣為:

在確定了上述狀態(tài)方程與量測(cè)方程及其相應(yīng)的雅可比矩陣后，即可應(yīng)用EKF進(jìn)行濾波。

3 比例導(dǎo)引制導(dǎo)

假定在上述濾波過(guò)程中已獲得來(lái)襲目標(biāo)的運(yùn)動(dòng)規(guī)律且攔截彈的速度變化規(guī)律已知，如何將導(dǎo)彈引向目標(biāo)便是制導(dǎo)規(guī)律所要完成的任務(wù)。比例導(dǎo)引法是一種最常用的、制導(dǎo)精度較高的制導(dǎo)規(guī)律。追蹤法、前置角法、平行接近法都可以看作是比例導(dǎo)引法的特殊情況。

比例導(dǎo)引法是指攔截彈飛行過(guò)程中速度矢量Vm的轉(zhuǎn)動(dòng)角速度dθm/dt與導(dǎo)彈、目標(biāo)連線的旋轉(zhuǎn)角速度dq/dt成比例的導(dǎo)引方法［9］。圖4與圖5分別為導(dǎo)彈與目標(biāo)相對(duì)位置和比例導(dǎo)引示意圖。

圖4 導(dǎo)彈、目標(biāo)相對(duì)位置

圖5 比例導(dǎo)引示意圖

圖中 Mk－1、Mk、Tk－1、Tk分別為 k － 1時(shí)刻、k時(shí)刻導(dǎo)彈位置和k－1時(shí)刻、k時(shí)刻目標(biāo)位置。A為Mk－1Mk與 Tk－1Tk交點(diǎn)，MkB 平行于 Mk－1Tk－1。

設(shè) xt、yt、zt表示目標(biāo)坐標(biāo)，任意時(shí)刻的 xt、yt、zt已知。算法流程如下:

k－1時(shí)刻導(dǎo)彈與目標(biāo)的距離:

k－1時(shí)刻目標(biāo)速度矢量與導(dǎo)彈、目標(biāo)連線的夾角:

k時(shí)刻導(dǎo)彈、目標(biāo)連線與基準(zhǔn)線夾角增量:為k－1時(shí)刻導(dǎo)彈與k時(shí)刻目標(biāo)之間距離，st=vt（k－1）·Δt，即在Δt內(nèi)將目標(biāo)視為做勻速直線運(yùn)動(dòng)。

比例導(dǎo)引法的差分方程為:

利用余弦定理求得:

設(shè)xA、yA、zA表示A點(diǎn)坐標(biāo)，xm、ym、zm表示導(dǎo)彈坐標(biāo)，由幾何關(guān)系知:

其中sm=vm（k－1）·Δt，即在Δt內(nèi)將導(dǎo)彈視為做勻速直線運(yùn)動(dòng)。在獲得k時(shí)刻的導(dǎo)彈坐標(biāo)后即可更新r（k）進(jìn)行下一次循環(huán)。必要時(shí)還要利用c3對(duì)dq進(jìn)行校正。

4 仿真分析

4.1 跟蹤濾波結(jié)果

假定雷達(dá)在跟蹤再入段時(shí)的量測(cè)標(biāo)準(zhǔn)差為:σR=10m，σA= σE=1mrad。跟蹤數(shù)據(jù)率設(shè)定為100Hz。采用兩點(diǎn)法初始化狀態(tài)變量，初始協(xié)方差可以參考文獻(xiàn)［10］。在濾波過(guò)程中，狀態(tài)噪聲為設(shè)計(jì)參數(shù)，文中采用分段噪聲，先利用較大噪聲使算法快速收斂，再采用較小噪聲降低濾波誤差。對(duì)目標(biāo)的跟蹤精度采用均方根誤差（RMSE）作為判別標(biāo)準(zhǔn)［11］，設(shè)X^n（ ）k|k為k時(shí)刻第n次蒙特卡羅的狀態(tài)估計(jì)值，則:

對(duì)量測(cè)數(shù)據(jù)進(jìn)行一次濾波，得到真實(shí)軌跡、量測(cè)點(diǎn)跡與濾波軌跡的結(jié)果見圖6。蒙特卡羅仿真500次得到均方根誤差結(jié)果如圖7。

圖6 導(dǎo)彈再入段的真實(shí)、量測(cè)、濾波軌跡

為方便對(duì)比，圖7中給出了針對(duì)上述相同軌跡應(yīng)用CV模型的濾波結(jié)果。由圖可以看出算法濾波收斂速度較快，誤差較小且收斂后基本穩(wěn)定。較CV模型濾波精度有較大提升。合理改變跟蹤數(shù)據(jù)率或雷達(dá)量測(cè)誤差，在不同的彈道模型下進(jìn)行多次仿真比較，結(jié)果類似。

圖7 再入段濾波的位置RMSE誤差

4.2 制導(dǎo)結(jié)果

設(shè)定比例系數(shù)為K=5，導(dǎo)彈初始位置為坐標(biāo)原點(diǎn)，速度恒定為 3000m/s，時(shí)間間隔 Δt=0.01s，導(dǎo)彈與目標(biāo)之間相對(duì)距離小于60m時(shí)認(rèn)為制導(dǎo)成功。獲得三維理想攔截彈彈道仿真圖見圖8，圖9為攔截彈的法向過(guò)載。

圖8 理想三維彈道仿真

圖9 攔截彈法向過(guò)載

5 結(jié)束語(yǔ)

文中針對(duì)再入大氣層的來(lái)襲導(dǎo)彈目標(biāo)受力情況推導(dǎo)了ENU坐標(biāo)系下目標(biāo)運(yùn)動(dòng)加速度表達(dá)式，建立了混合坐標(biāo)系下的目標(biāo)狀態(tài)方程與量測(cè)方程，通過(guò)受力分析采用基于其特性的彈道模型進(jìn)行濾波，給出了狀態(tài)估計(jì)與協(xié)方差矩陣的初始值計(jì)算方法，采用EKF濾波算法對(duì)目標(biāo)進(jìn)行跟蹤，為制導(dǎo)過(guò)程提供了來(lái)襲導(dǎo)彈的速度位置信息。利用比例導(dǎo)引制導(dǎo)律進(jìn)行仿真得到三維制導(dǎo)圖。仿真結(jié)果表明，算法可以有效的對(duì)再入類導(dǎo)彈進(jìn)行穩(wěn)定的跟蹤和精確的制導(dǎo)。

［1］高尚.比例導(dǎo)引理想彈道仿真［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2003，24（8）:66－68.

［2］陳宇強(qiáng)，趙育善.導(dǎo)彈最優(yōu)制導(dǎo)律彈道仿真分析［J］.指揮控制與仿真，2007，29（3）:92－105.

［3］周紀(jì)元，童幼堂，張磊，等.典型導(dǎo)引規(guī)律三維彈道仿真分析［J］.艦船電子工程，2008，28（2）:110－112.

［4］張毅，肖龍旭，王順宏.彈道導(dǎo)彈彈道學(xué)［M］.長(zhǎng)沙:國(guó)防科技大學(xué)出版社，2005.

［5］楊炳尉.標(biāo)準(zhǔn)大氣參數(shù)的公式表示［J］.宇航學(xué)報(bào)，1983（1）:83－86.

［6］何子述，夏威.現(xiàn)代數(shù)字信號(hào)處理及其應(yīng)用［M］.北京:清華大學(xué)出版社，2009.

［7］易令，呂明.使用交互多模型算法的高速高機(jī)動(dòng)目標(biāo)跟蹤［J］.雷達(dá)科學(xué)與技術(shù)，2006，4（3）:143－147.

［8］辛召?gòu)?qiáng)，沈曉峰.利用CUDA快速實(shí)現(xiàn)IMM目標(biāo)跟蹤［J］.雷達(dá)科學(xué)與技術(shù)，2012，10（6）:656－659.

［9］梁慜.彈道導(dǎo)彈攔截仿真建模技術(shù)研究［J］.中國(guó)電子科學(xué)研究院學(xué)報(bào)，2013，8（1）:56－59.

［10］何友，修建娟，張晶煒.雷達(dá)數(shù)據(jù)處理及應(yīng)用［M］.北京:電子工業(yè)出版社，2006.

［11］張丕旭，石章松，劉忠.一種基于彈道模型的機(jī)動(dòng)目標(biāo)跟蹤算法［J］.彈箭與制導(dǎo)學(xué)報(bào)，2009，29（4）:55－57.