基于動(dòng)力學(xué)模型的彈道式再入彈道估計(jì)方法*

江曉東，謝京穩(wěn)，郭軍海

(北京跟蹤與通信技術(shù)研究所，北京 100094)

0 引言

飛行器再入大氣飛行的特點(diǎn)是速度快，并由于受到復(fù)雜的空氣動(dòng)力影響而呈現(xiàn)出強(qiáng)的非線性特性[2]。Julier和Uhlmann[3-4]等利用UT(unscented transform)變換，提出unscented Kalman filter(UKF)，由于不需要對(duì)非線性系統(tǒng)進(jìn)行線性化，可有效地應(yīng)用到非線性系統(tǒng)的狀態(tài)估計(jì)中。文獻(xiàn)[5-6]在再入目標(biāo)彈道估計(jì)的背景下，對(duì)3種非線性濾波——擴(kuò)展卡爾曼濾波(extended Kalman filter，EKF)，UKF和粒子濾波(particles filter，PF)的濾波性能進(jìn)行了對(duì)比分析，從計(jì)算量、濾波性能和濾波的穩(wěn)健性上綜合考察，UKF是這3種濾波算法中最優(yōu)的。

1 彈道式再入目標(biāo)的動(dòng)力學(xué)模型

彈道式再入飛行器在再入過程中只受到空氣阻力和重力的作用，在地心慣性坐標(biāo)系[7]下，其總加速度的表達(dá)式為[7]：

a=aD+aG+w,

(1)

式中:a為目標(biāo)的總加速度；aD為空氣阻力引起的加速度；aG為地心引力引起的加速度；w為系統(tǒng)狀態(tài)噪聲，假定服從高斯分布，其方差為模型設(shè)計(jì)參數(shù)。

在地心地固坐標(biāo)系[7]下，加速度表達(dá)式在式(1)的基礎(chǔ)上修正為

(2)

式中:ar為目標(biāo)相對(duì)總加速度;ωe為地球自轉(zhuǎn)角速度向量；p為目標(biāo)的位置向量；v為再入體相對(duì)大氣的速度向量；2ωe×v為科里奧利力引起的加速度；ωe×(ωe×p)為慣性離心力引起的加速度。

空氣阻力的方向與目標(biāo)相對(duì)于空氣的速度矢量v相反，與空氣密度ρ和目標(biāo)速度v的平方成正比，空氣阻力引起的加速度為

(3)

式中：uv=v/|v|為vr的單位向量；ρ(h)為空氣密度，h為高程，本文中采用文獻(xiàn)[8]提出的大氣密度模型。

α(t)為阻力參數(shù)，其倒數(shù)β(t)=1/α(t)也稱彈道系數(shù)。對(duì)彈道式再入飛行器來說，α(t)是一個(gè)非常重要的參數(shù)，它通常是未知變量，因此經(jīng)常也將其作為待估狀態(tài)參量，增加到運(yùn)動(dòng)狀態(tài)向量中，使得狀態(tài)向量從6維變?yōu)?維X=(x,y,z,vx,vy,vz,α)T。α(t)的建模方式有多種，適用于不同的環(huán)境。文獻(xiàn)[7]對(duì)其模型進(jìn)行了系統(tǒng)的總結(jié)。本文選用指數(shù)模型，即對(duì)阻力參數(shù)的比率γln(α/α0)進(jìn)行建模：

α(t)=α0exp(γ(t)),

(4)

(5)

式中:ωγ(t)為高斯白噪聲；α0為預(yù)先知道的典型阻力參數(shù)。這樣，待估狀態(tài)向量變?yōu)閄=(x,y,z,vx,vy,vz,γ)T。這種建模方式需要先驗(yàn)信息α0，但保證了阻力參數(shù)恒正，這對(duì)于保證濾波的穩(wěn)定性至關(guān)重要[3]。

基于式(1)和(3)，將加速度的矢量形式寫成分量形式，整理可得

(6)

最后，基于方程(5)和(6)可得彈道式再入飛行器的動(dòng)力學(xué)模型為

(7)

將此模型記為BRV-Exp模型，式(7)為變系數(shù)非線性常微分方程組，難以給出解析解，但可借助歐拉積分法或龍格庫塔法給出數(shù)值解，以便進(jìn)行狀態(tài)參量的預(yù)測(cè)。

2 多臺(tái)雷達(dá)的測(cè)量模型

單臺(tái)雷達(dá)的測(cè)量模型[9]為

(8)

記為

Y1(k)=h1(X(k))+v1(k),

(9)

基于N臺(tái)雷達(dá)的總測(cè)量模型為[10]

Y(k)=h(X(k))+v(k),

(10)

式中：Y(k)=(Y1(k),Y2(k),…,YN(k))T；

h(·)=(h1(·),h2(·),…,hN(·))T；

v(k)=(v1(k),v2(k),…,v3(k))T.

設(shè)v(k)服從正態(tài)分布N(0,R(k))，R(k)=diag(R1(k),R2(k),…,RN(k))，這樣基于狀態(tài)轉(zhuǎn)換模型(7)和測(cè)量模型(10)，就可利用UKF進(jìn)行再入目標(biāo)的彈道估計(jì)了。

3 不敏卡爾曼濾波(UKF)

UKF是由Julier等人首先提出的。UT變換是UKF的基礎(chǔ)，它通過構(gòu)造一個(gè)確定的點(diǎn)集(Sigma點(diǎn))，使得它與當(dāng)前狀態(tài)具有相同的統(tǒng)計(jì)特征，即均值和方差相同。具體來說，UT變換利用當(dāng)前狀態(tài)(均值和方差)構(gòu)造點(diǎn)集{Xi(k)}，然后讓每個(gè)點(diǎn)Xi(k)通過非線性方程，得到非線性變換后的點(diǎn)集{Yi(k)}，預(yù)測(cè)均值和方差可以在{Yi(k)}點(diǎn)集上求得。

設(shè)待估狀態(tài)向量維數(shù)為L，則Sigma點(diǎn)的個(gè)數(shù)為一般為2L+1，為了增加Sigma點(diǎn)的個(gè)數(shù)，以提高估計(jì)精度，通常需要對(duì)狀態(tài)向量進(jìn)行擴(kuò)維，擴(kuò)維后的狀態(tài)向量由目標(biāo)狀態(tài)向量，狀態(tài)噪聲和測(cè)量噪聲組成，即Xa=(XT,vT,nT)。擴(kuò)維UKF(augmented UKF，AUKF)算法流程如下[11]：

初始化：

(11)

(12)

(13)

(14)

計(jì)算Sigma點(diǎn)及其權(quán)值：

(15)

(16)

(17)

(18)

預(yù)測(cè)：

(19)

(20)

(21)

(22)

(23)

更新：

(24)

(25)

(26)

(27)

(28)

4 仿真結(jié)果與分析

圖1 再入目標(biāo)的軌跡及測(cè)量雷達(dá)布局圖Fig.1 Trajectory of re-entry target and measuring-radar layout

圖2 位置分量殘差的均值對(duì)比圖Fig.2 Comparison chart of position residuals error mean

圖3 速度分量殘差的均值對(duì)比圖Fig.3 Comparison chart of velocity residuals error mean

圖4 位置分量的均方根誤差對(duì)比圖Fig.4 Comparison chart of position root-mean-square error

圖5 速度分量的均方根誤差對(duì)比圖Fig.5 Comparison chart of velocity root-mean-square error

表1 均方根誤差的均值Table 1 Mean of velocity root-mean-square error

方法名稱Δx/mΔy/mΔz/mΔvx/(m·s-1)Δvy/(m·s-1)Δvz/(m·s-1)基線方向余弦法0.870 20.967 11.965 80.055 40.094 40.156 4動(dòng)力學(xué)模型法 0.133 70.135 20.253 60.029 60.036 70.070 1

5 結(jié)束語

參考文獻(xiàn)：

[1] 蔡洪．Unscented Kalman濾波用于再入飛行器跟蹤[J]．飛行器測(cè)控學(xué)報(bào)，2003，22(3)：12-16．

CAI Hong．Unscented Kalman Filtering for Reentry Vehicle Tracking[J]．Journal of Spacecraft TT&C Technology，2003，22(3)：12-16．

[2] 蔡洪，張金槐．導(dǎo)彈飛行試驗(yàn)多站測(cè)量定軌的自適應(yīng)濾波方法[J]．航天控制，1997(2)：17-22．

CAI Hong，ZHANG Jin-huai．Adaptive Filtering Approaeh for Missile Flight Test Trajectory Deterrulnation Under Multiple Station Measuring[J]．Aerospace Control，1997(2)：17-22．

[3] JULIER S J，UHLMANN J K．Unscented Filtering and Nonlinear Estimation[J]．Proc．IEEE，March 2004，92(3)：401-422．

[4] JULIER S J，UHLMANN J K，DURRANT-WHYTE H F．A New Method for Nonlinear Transformation of Means and Covariance in Filters and Estimators[J]．IEEE Trans．Automatic Control，Mar．1996，AC-45(3)：477-482．

[5] FARINA A，RISTIC B，BENVENUTI D． Tracking a Ballistic Target：Comparison of Several Nonlinear Filters[J].IEEE Transactions On Aerospace and Electronic Systems July 2002，38(3)：854-867．

[6] RISTIC B，F(xiàn)ARINA A，BENVENUTI D． Tracking a Ballistic-Entry Object：Performance Bounds and Comparison of Non-Linear Filters[C]∥IDC-2002，Adelaide，Australia，F(xiàn)ebruary 2002：259-264．

[7] LI X R,JILKOV V P．A Survey of Maneuvering Target Tracking-Part II：Ballistic Target Models[C]∥2001 SPIE Conf on Signal and Data Processing of Small Targets，San Diego,CA,USA,July August 2001,Vol．4473：559-581．

[8] 楊炳尉. 標(biāo)準(zhǔn)大氣參數(shù)的公式表示[J]. 宇航學(xué)報(bào)，1983,1(1)：83-86．

YANG Bing-wei．Formulization of Standard Atmospheric Parameters[J]. Journal of Astronautics，1983,1(1)：83-86．

[9] LI X R,JILKOV V P．A Survey of Maneuvering Target Tracking-Part Ⅲ：Measurement Models[C]∥2001 SPIE Conf on Signal and Data Processing of Small Targets，San Diego,CA,USA,July August 2001,4473： 423-446．

[10] 丁維福．多傳感器信息融合理論及其在機(jī)動(dòng)目標(biāo)跟蹤中的應(yīng)用[D]．西安：西北工業(yè)大學(xué)，2007：19．

DING Wei-fu．Multi-Sensor Information Fusion Theory and Its Applications to Maneuvering Target Tracking[D]．Xi’an：Northwestern Polytechnical University，2007：19．

[11] WAN E A，MERWE R M．The Unscented Kalman Filter [M/OL]．(2001-08-01)[2009-04-03]http:∥stomach.v2.n1/docs-/TechPubs/Tracking_and_AR/wan01/unscented.pdf．

[12] 胡紹林，許愛華，郭小紅．脈沖雷達(dá)跟蹤測(cè)量數(shù)據(jù)處理技術(shù)[M]．北京：國防工業(yè)出版社，2007：120-135．

HU Shao-lin，XU Ai-hua，GUO Xiao-hong ．Pulse Radar Tracking Measurement Data Processing Technology[M]．Beijing：National Defence Industry Press，2007：120-135．