CN104778376A - 一种临近空间高超声速滑翔弹头跳跃弹道预测方法 - Google Patents

Abstract

一种临近空间高超声速滑翔弹头跳跃弹道预测方法，本发明涉及一种飞行器弹道预测方法。本发明是要解决现有方法对机动目标弹道预测精度低的问题，而提供一种临近空间高超声速滑翔弹头跳跃弹道预测方法。它按下述步骤实现：一、建立高超声速滑翔弹头的弹道方程；二、设计实时跟踪高超声速滑翔弹头运动轨迹的卡尔曼滤波器；三、根据跟踪结束时刻的高超声速滑翔弹头的位置、速度和加速度，结合弹道方程估算高超声速滑翔弹头飞行时的攻角和滚转角，在随后的飞行时间内临近空间高超声速滑翔弹头进行等攻角和等滚转角飞行，应用弹道方程向下一时刻循环递推计算，得到一定时间之后的高超声速滑翔弹头的弹道预测值。属于目标跟踪技术领域。

Description

一种临近空间高超声速滑翔弹头跳跃弹道预测方法

技术领域

本发明涉及一种弹道预测方法，特别是一种针对临近空间高超声速滑翔弹头跳跃弹道的预测方法。

背景技术

临近空间高超声速滑翔弹头是一种能够在临近空间无动力滑翔数千至一万余公里，并具有较强的侧向机动能力和突防性能，能携带核弹头或常规弹头实施远距离快速打击,这类飞行器由于具有较高的升阻比，且在大气层内长时间飞行，其运动轨迹往往呈现出“跳跃”特征。

从国内外目前弹道预报的研究情况来看，国外对弹道预测方法的研究仅限于弹道导弹方面，对机动目标预测方法的研究未见公开报道。而国内对弹道预测方法的研究也主要集中在弹道导弹方面，对机动目标弹道预测方法的仅有的研究结果是采用目标机动当前统计模型，利用卡尔曼预报原理向前简单地推算,没有考虑目标的气动力等复杂情况。总之，目前还没有对高超声速滑翔弹头进行高精度弹道预测的方法。

发明内容

本发明是要解决现有方法对机动目标弹道预测精度低的问题，而提供一种临近空间高超声速滑翔弹头跳跃弹道预测方法。

一种临近空间高超声速滑翔弹头跳跃弹道预测方法，它按以下步骤实现：

一、高超声速滑翔弹头的弹道方程的建立；

建立临近空间高超声速滑翔弹头的弹道方程，包括气动力计算方程、加速度计算方程、速度和位置计算方程、速率和弹道角计算方程，以及姿态角计算方程；

综合计算方程，得到高超声速滑翔弹头质心运动速度和位置方程如式(1)～(3)所示：

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_I=v_{\mathrm{xI}}\\ {\stackrel{\·}{v}}_{\mathrm{xI}}=-\frac{\mathrm{\μ}{x}_I}{r_I^3}+a_{\mathrm{axI}}\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}{\stackrel{\·}{y}}_I=v_{\mathrm{yI}}\\ {\stackrel{\·}{v}}_{\mathrm{yI}}=-\frac{\mathrm{\μ}{y}_I}{r_I^3}+a_{\mathrm{ayI}}\end{array}\right.---\left(2\right)$

$\left\{\begin{array}{c}{\stackrel{\·}{z}}_I={\stackrel{\·}{v}}_{\mathrm{zI}}\\ {\stackrel{\·}{v}}_{\mathrm{zI}}=-\frac{{\mathrm{\μz}}_I}{r_I^3}+a_{\mathrm{azI}}\end{array}\right.---\left(3\right)$

其中，x_I、y_I和z_I代表高超声速滑翔弹头在地心惯性坐标系中的位置，代表高超声速滑翔弹头的位置到地球中心的距离,和则代表高超声速滑翔弹头在地心惯性坐标系中的速度，a_axI、a_ayI和a_azI代表气动力产生的高超声速滑翔弹头机动加速度在地心惯性坐标系中的分量，μ为引力参数；

二、高超声速滑翔弹头的跟踪滤波器设计；

设计高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器，实时地跟踪并估计高超声速滑翔弹头的位置、速度和加速度；

其中，所述高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器为

$\left\{\begin{array}{c}{\stackrel{\‾}{X}}_i(k+1)={\mathrm{\Φ}}_i{\hat{X}}_i\left(k\right)\\ P_i(k+1/k){=\mathrm{\Φ}}_i{P}_i\left(k\right){\mathrm{\Φ}}_i^T+Q_i\end{array}\right.,i=x,y,z---\left(4\right)$

上式中，代表滤波器i第k步的状态估计值，代表滤波器i从第k步到第k+1步的状态预报值，P_i(k)代表滤波器i第k步的状态估计误差协方差矩阵，P_i(k+1/k)代表滤波器i从第k步到第k+1步的状态预报误差协方差矩阵，Q_i为滤波器i模型预报误差协方差矩阵；

${\mathrm{\Φ}}_i=\exp \left(F_i\mathrm{\Δt}\right)=\left[\begin{array}{ccc}1& \mathrm{\Δt}& 1/{\mathrm{\λ}}_i^2(e^{-{\mathrm{\λ}}_i\mathrm{\Δt}}+{\mathrm{\λ}}_i\mathrm{\Δt}-1)\\ 0& 1& 1/{\mathrm{\λ}}_i(1-e^{-{\mathrm{\λ}}_i\mathrm{\Δt}})\\ 0& 0& e^{-{\mathrm{\λ}}_i\mathrm{\Δt}}\end{array}\right]---\left(5\right)$

所述运动轨迹跟踪卡尔曼滤波器的测量修正方程为

$\left\{\begin{array}{c}{K}_i(k+1)=P_i(k+1/k)H_i^T{[J_i{P}_i(k+1/k)H_i^T+R_i]}^{-1}\\ {\hat{X}}_i(k+1)={\stackrel{\‾}{X}}_i(k+1)+K_i(k+1)[{\mathrm{\η}}_i(k+1)-H_i{\stackrel{\‾}{X}}_i(k+1)],i=x,y,z\\ P_i(k+1)=[I-K_i(k+1)H_i]P_i(k+1/k)\end{array}\right.---\left(6\right)$

其中，上标^T代表矩阵转置运算，上标^-1代表矩阵求逆运算，滤波器i的测量矩阵

H_i＝[1 0 0],i＝x,y,z          (7)

R_i为滤波器i测量误差协方差矩阵，I代表单位矩，K_i(k+1)代表滤波器i第k+1步的滤波增益矩阵，η_i(k+1)代表滤波器i第k+1步的测量信息，P_i(k+1)代表滤波器i第k+1步的状态估计误差协方差矩阵；

三、高超声速滑翔弹头的弹道预测；

(三一)根据步骤二中运动轨迹跟踪卡尔曼滤波器跟踪结束kT时刻的高超声速滑翔弹头的加速度，其中，k为大于0的整数，T为计算周期，一般可取T＝0.01s，应用步骤一中的临近空间高超声速滑翔弹头的加速度计算方程求出跟踪结束kT时刻临近空间高超声速滑翔弹头的气动力产生的加速度，并根据跟踪结束kT时刻高超声速滑翔弹头的速度，应用步骤一中的临近空间高超声速滑翔弹头的弹道角计算方程计算其弹道角；

(三二)根据计算出来的kT时刻临近空间高超声速滑翔弹头的气动力产生的加速度和弹道角,结合步骤一中的临近空间高超声速滑翔弹头的气动力计算方程和姿态角计算方程估算跟踪结束kT时刻高超声速滑翔弹头飞行时的攻角和姿态角；

(三三)设在随后的kT到(k+1)T，(k+2)T,...,(k+N)T时间内，N为远大于2的正整数，临近空间高超声速滑翔弹头进行等攻角和等滚转角飞行，则纵向平面内呈现出上下跳跃式的飞行弹道；同时，高超声速滑翔弹头采用BTT控制方式，侧滑角β等于零，通过控制滚转角在其体系侧向输出一个机动加速度进行转弯，即滚转角指令γ_c保持常值；

(三四)根据等攻角、零侧滑角和等滚转角飞行条件，计算kT时刻气动力产生的加速度、kT时刻的位置x_I、y_I和z_I、kT时刻的速度v_xI、v_yI和v_zI代入式(1)～(3)，经过一步数值积分运算后得到高超声速滑翔弹头在(k+1)T时刻在地心惯性坐标系中的速度和位置；

计算高超声速滑翔弹头在(k+1)T时刻受到的气动力，计算高超声速滑翔弹头在(k+1)T时刻受到的气动力加速度，再把(k+1)T时刻的气动力加速度以及(k+1)T时刻的位置和速度代入(1)～(3)，计算出(k+2)T时刻的位置和速度，以此类推，计算出直到(k+N)T时刻的位置和速度。

发明效果：

在上述条件下，经过100次Monte-Carlo仿真得到的高超声速滑翔弹头最终位置预报的统计结果如图13所示，统计计算得出终端位置预报误差小于50km的概率为95％。

针对临近空间高超声速滑翔弹头，本发明考虑了目标受气动力等复杂情况，提出了一种新型跳跃弹道预测方法，相比于传统方法提高了弹道预测精度，减小了弹道预测误差，填补了临近空间高超声速飞行器弹道预报领域的空白。

附图说明

图1为具体实施方式二中地心惯性坐标系和探测点惯性坐标系图；

图2为仿真实验中高超声速滑翔弹头飞行弹道图，上图为纵向平面飞行弹道,下图为侧向平面飞行弹道；

图3为仿真实验中高超声速滑翔弹头X轴加速度估计图；

图4为仿真实验中高超声速滑翔弹头Y轴加速度估计图；

图5为仿真实验中高超声速滑翔弹头Z轴加速度估计图；

图6为仿真实验中高超声速滑翔弹头X轴速度估计图；

图7为仿真实验中高超声速滑翔弹头Y轴速度估计图；

图8为仿真实验中高超声速滑翔弹头Z轴速度估计图；

图9为仿真实验中高超声速滑翔弹头X轴位置估计图，上图为X轴位置的真实值和估计值,下图为X轴位置的真实值与估计值之差；

图10为仿真实验中高超声速滑翔弹头Y轴位置估计图，上图为Y轴位置的真实值和估计值,下图为Y轴位置的真实值与估计值之差；

图11为仿真实验中高超声速滑翔弹头Z轴位置估计图，上图为Z轴位置的真实值和估计值,下图为Z轴位置的真实值与估计值之差；

图12为仿真实验中高超声速滑翔弹头位置估计误差图；

图13为仿真实验中高超声速滑翔弹头终端位置预报误差统计结果。

具体实施方式

具体实施方式一：本实施方式的一种临近空间高超声速滑翔弹头跳跃弹道预测方法，它按以下步骤实现：

一、高超声速滑翔弹头的弹道方程的建立；

建立临近空间高超声速滑翔弹头的弹道方程，包括气动力计算方程、加速度计算方程、速度和位置计算方程、速率和弹道角计算方程，以及姿态角计算方程；

综合计算方程，得到高超声速滑翔弹头质心运动速度和位置方程如式(1)～(3)。

二、高超声速滑翔弹头的跟踪滤波器设计；

设计高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器，实时地跟踪并估计高超声速滑翔弹头的位置、速度和加速度；

其中，所述高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器的状态一步预报方程为(4)，所述运动轨迹跟踪卡尔曼滤波器的测量修正方程为(6)；

三、高超声速滑翔弹头的弹道预测；

(三一)根据步骤二中运动轨迹跟踪卡尔曼滤波器跟踪结束kT时刻的高超声速滑翔弹头的加速度，其中，k为大于0的整数，T为计算周期，一般可取T＝0.01s，应用步骤一中的临近空间高超声速滑翔弹头的加速度计算方程求出跟踪结束kT时刻临近空间高超声速滑翔弹头的气动力产生的加速度，并根据跟踪结束kT时刻高超声速滑翔弹头的速度，应用步骤一中的临近空间高超声速滑翔弹头的弹道角计算方程计算其弹道角；

(三二)根据计算出来的kT时刻临近空间高超声速滑翔弹头的气动力产生的加速度和弹道角,结合步骤一中的临近空间高超声速滑翔弹头的气动力计算方程和姿态角计算方程估算跟踪结束kT时刻高超声速滑翔弹头飞行时的攻角和姿态角；

(三三)设在随后的kT到(k+1)T，(k+2)T,...,(k+N)T时间内，N为远大于2的正整数，临近空间高超声速滑翔弹头进行等攻角和等滚转角飞行，则纵向平面内呈现出上下跳跃式的飞行弹道；同时，高超声速滑翔弹头采用BTT控制方式，侧滑角β等于零，通过控制滚转角在其体系侧向输出一个机动加速度进行转弯，即滚转角指令γ_c保持常值；

(三四)根据等攻角、零侧滑角和等滚转角飞行条件，计算kT时刻气动力产生的加速度、kT时刻的位置x_I、y_I和z_I、kT时刻的速度v_xI、v_yI和v_zI代入式(1)～(3)，经过一步数值积分运算后得到高超声速滑翔弹头在(k+1)T时刻在地心惯性坐标系中的速度和位置；

计算高超声速滑翔弹头在(k+1)T时刻受到的气动力，计算高超声速滑翔弹头在(k+1)T时刻受到的气动力加速度，再把(k+1)T时刻的气动力加速度以及(k+1)T时刻的位置和速度代入式(1)～(3)，计算出(k+2)T时刻的位置和速度，以此类推，计算出直到(k+N)T时刻的位置和速度。

具体实施方式二：本实施方式与具体实施方式一不同的是：

第一步坐标系的定义

为了得出临近空间高超声速滑翔弹头跳跃弹道的预测方法，需要建立一组方程来描述临近空间高超声速滑翔弹头的运动，为此，首先需要定义几个有关的坐标系。

(一)地心惯性坐标系(o_Ix_Iy_Iz_I)

如图1所示，原点o_I位于地心上，o_Iy_I指向地面雷达探测点，o_Ix_I轴位于目标来袭平面内，与o_Iy_I轴垂直，指向目标为正，o_Iz_I按右手定则确定。

(二)探测点惯性坐标系(o_Fxyz)

原点为地面雷达所在点o_F，o_Fy轴铅锤向上，o_Fx轴在目标来袭平面内，垂直于o_Fy轴，指向目标为正，o_Fz轴按右手定则确定。此坐标系在导弹发射瞬间固定于惯性空间。它与地心惯性坐标系的关系如图1所示，即o_Fx、o_Fy和o_Fz轴分别平行于o_Ix_I、o_Iy_I和o_Iz_I轴。

(三)目标惯性坐标系(o_tFx_ty_tz_t)

原点为目标发射点o_tF，将导弹发射点惯性坐标系绕着o_Fy轴旋转180°，即得到此坐标系。此坐标系用于确定目标的飞行弹道角。

(四)目标体坐标系(ox₁y₁z₁)

原点o位于目标质心上。ox₁轴与其纵轴重合，指向头部为正，oy₁轴在其纵向对称平面内，垂直ox₁轴，指向上方为正，oz₁轴方向按右手定则确定。

(五)几个坐标系之间的转换关系

(六)由探测点惯性坐标系到目标惯性坐标系的变换矩阵

$C_{0\→t0}=\left[\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$

(七)由目标惯性坐标系到目标体坐标系的变换矩阵

采用231转序推导出来的坐标转换矩阵为

其中，、ψ和γ分别为目标的俯仰角，偏航角和滚动角。

步骤一具体为：

(一)计算高超声速滑翔弹头气动力

在滑翔阶段，高超声速滑翔弹头受到的气动力在体系内的投影按照下式计算：

$\left\{\begin{array}{c}{R}_{x1}=-c_x\mathrm{qS}\\ R_{y1}=c_y\mathrm{sign}\left(\mathrm{\α}\right)\mathrm{qS}\\ R_{z1}=0\end{array}\right.---\left(8\right)$

其中，α为攻角，

$\mathrm{sign}\left(\mathrm{\&alpha;}\right)=\left\{\begin{array}{c}1,\mathrm{\&alpha;}>0\\ -1,\mathrm{\&alpha;}<0\end{array}\right.$

c_y为升力系数，c_x为阻力系数，它们与高超声速滑翔弹头的速度相关，分别需要查表1和表2求得。表1和表2中的Ma代表马赫数，即高超声速滑翔弹头的速率v除以音速，高超声速滑翔弹头在临近空间滑翔时的速率大约为15Ma；S为特征参考面积，设高超声速滑翔弹头的特征参考面积为S＝0.4837m²；

q＝ρv²/2        (9)

其中，q为动压头，ρ为空气密度，它随着飞行高度h变化，可查标准的大气密度表得到；高超声速滑翔弹头飞行高度h用下式计算:

$h=\sqrt{x_I^2+y_I^2+z_I^2}-R_e---\left(10\right)$

其中，x_I、y_I和z_I代表高超声速滑翔弹头在地心惯性坐标系中的位置，R_e＝6378km为地球平均半径；

表1 高超声速滑翔弹头的升力系数c_y

表2 高超声速滑翔弹头的阻力系数c_x

设在飞行中侧滑角β＝0，因此气动力侧向力R_z1＝0；

(二)计算高超声速滑翔弹头的速率

高超声速滑翔弹头的速率为

$v=\sqrt{v_{\mathrm{xI}}^2+v_{\mathrm{yI}}^2+v_{\mathrm{zI}}^2}---\left(11\right)$

其中，v_xI、v_yI和v_zI代表高超声速滑翔弹头在地心惯性坐标系中的速度；

将高超声速滑翔弹头的速度矢量投影到目标惯性坐标系v_xt0、v_yt0与v_zt0，即

$\left[\begin{array}{c}{v}_{\mathrm{xt}0}\\ v_{\mathrm{yt}0}\\ v_{\mathrm{zt}0}\end{array}\right]=C_{0\&RightArrow;t0}\left[\begin{array}{c}{v}_{\mathrm{xI}}\\ v_{\mathrm{yI}}\\ v_{\mathrm{zI}}\end{array}\right]---\left(12\right)$

其中，

$C_{0\&RightArrow;t0}=\left[\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}\right];$

(三)计算高超声速滑翔弹头的弹道仰角θ和弹道偏角ψ_v：

$\left\{\begin{array}{c}\mathrm{\&theta;}={\sin }^{-1}(v_{\mathrm{yt}0}/v)\\ {\mathrm{\&psi;}}_v=-{\tan }^{-1}(v_{\mathrm{zt}0}/v_{\mathrm{xt}0})\end{array}\right.---\left(13\right)$

高超声速滑翔弹头的姿态角、ψ和γ按照下式计算：

其中，为俯仰角，ψ为偏航角，γ为滚转角，γ_c为滚转角指令；

(四)计算高超声速滑翔弹头机动加速度在地心惯性坐标系中的投影

$\left[\begin{array}{c}{a}_{\mathrm{axI}}\\ a_{\mathrm{ayI}}\\ a_{\mathrm{azI}}\end{array}\right]=C_{0\&RightArrow;t0}^T{C}_{t0\&RightArrow;1}^T\left[\begin{array}{c}{R}_{x1}/m\\ R_{y1}/m\\ R_{z1}/m\end{array}\right]---\left(15\right)$

其中，

m为高超声速滑翔弹头质量，假设m＝907.2kg；

(五)高超声速滑翔弹头质心运动速度和位置方程如式(1)～(3)所示，其中包含了高超声速滑翔弹头受到的重力加速度，即

$\left\{\begin{array}{c}{g}_{\mathrm{xI}}=-\frac{{\mathrm{\&mu;x}}_I}{r_I^3}\\ g_{\mathrm{yI}}=-\frac{\mathrm{\&mu;}{y}_I}{r_I^3}\\ g_{\mathrm{zI}}=-\frac{\mathrm{\&mu;}{z}_I}{r_I^3}\end{array}\right.---\left(16\right)$

其中，g_xI、g_yI与g_zI代表重力加速度在地心惯性系中的投影。

当目标高超声速滑翔弹头的气动力产生的机动加速度已经充分表现出来，跟踪滤波器已经进入稳态，充分地估计出来目标的机动加速度之后，转入采用弹道预报算法进行预测。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤二具体为：

(一)在探测点惯性坐标系的三个轴上，均用Singer模型来分别描述高超声速滑翔弹头的加速度分量的变化，即

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{a}}_x=-{\mathrm{\&lambda;}}_x{a}_x+w_x\\ {\stackrel{\&CenterDot;}{a}}_y=-{\mathrm{\&lambda;}}_y{a}_y+w_y\\ {\stackrel{\&CenterDot;}{a}}_z=-{\mathrm{\&lambda;}}_z{a}_z+w_z\end{array}\right.---\left(17\right)$

其中，λ_x、λ_y和λ_z分别代表在探测点惯性坐标系下o_Fx轴、o_Fy轴和o_Fz轴方向上目标机动时间常数的倒数，w_x、w_y和w_z分别代表探测点惯性坐标系下o_Fx轴、o_Fy轴和o_Fz轴方向上的零均值高斯白噪声；加速度a_x、a_y和a_z中既包含气动力产生的加速度a_axI、a_ayI和a_azI，又包含重力产生的加速度；探测点惯性坐标系的o_Fx轴、o_Fy轴和o_Fz轴与地心惯性系的三个轴o_Ix_I、o_Iy_I和o_Iz_I分别平行，气动力产生的加速度在地心惯性系中的投影a_axI、a_ayI和a_azI就是此加速度在探测点惯性系中的投影；

(二)在探测点惯性坐标系下建立起如下三组运动方程

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=v_x\\ {\stackrel{\&CenterDot;}{v}}_x=a_x\\ {\stackrel{\&CenterDot;}{a}}_x=-{\mathrm{\&lambda;}}_x{a}_x+w_x\end{array}\right.---\left(18\right)$

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{y}=v_y\\ {\stackrel{\&CenterDot;}{v}}_y=a_y\\ {\stackrel{\&CenterDot;}{a}}_y=-{\mathrm{\&lambda;}}_y{a}_y+w_y\end{array}\right.---\left(19\right)$

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{z}=v_z\\ {\stackrel{\&CenterDot;}{v}}_z=a_z\\ {\stackrel{\&CenterDot;}{a}}_z=-{\mathrm{\&lambda;}}_z{a}_z+w_z\end{array}\right.---\left(20\right)$

其中，x、y和z代表高超声速滑翔弹头在探测点惯性坐标系中的位置，和则代表高超声速滑翔弹头速度在探测点惯性坐标系中的投影。因为探测点惯性坐标系的o_Fx轴、o_Fy轴和o_Fz轴与地心惯性系的三个轴o_Ix_I、o_Iy_I和o_Iz_I分别平行，所以

$\left\{\begin{array}{c}{v}_x=v_{\mathrm{xI}}\\ v_y=v_{\mathrm{yI}}\\ v_z=v_{\mathrm{zI}}\end{array}\right.---\left(21\right)$

而由于探测点惯性坐标系只是相对于地心惯性系探测点惯性坐标系在o_Iy_I轴方向上作了平移，所以有

$\left\{\begin{array}{c}x=x_I\\ y=y_I-R_e\\ z=z_I\end{array}\right.---\left(22\right)$

对于***(18)，定义x轴子***状态变量X_x＝[x v_x a_x]^T，则得到如下状态方程

${\stackrel{\&CenterDot;}{X}}_x=F_x{X}_x+G_x{w}_x---\left(23\right)$

其中，

$F_x=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -{\mathrm{\&lambda;}}_x\end{array}\right],G_x=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]---\left(24\right)$

对式(23)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_x(k+1)＝Φ_xX_x(k)+ω_x(k)       (25)

其中，

${\mathrm{\&Phi;}}_x=\exp \left(F_x\mathrm{\&Delta;t}\right)=\left[\begin{array}{ccc}1& \mathrm{\&Delta;t}& 1/{\mathrm{\&lambda;}}_x^2(e^{-{\mathrm{\&lambda;}}_x\mathrm{\&Delta;t}}+{\mathrm{\&lambda;}}_x\mathrm{\&Delta;t}-1)\\ 0& 1& 1/{\mathrm{\&lambda;}}_x(1-e^{-{\mathrm{\&lambda;}}_x\mathrm{\&Delta;t}})\\ 0& 0& e^{-{\mathrm{\&lambda;}}_x\mathrm{\&Delta;t}}\end{array}\right]---\left(26\right)$

ω_x(k)代表x轴子***状态噪声向量，为零均值高斯白噪声向量。

高超声速滑翔弹头的位置由地面目标跟踪装置测量得到，则测量方程为

η_x＝x+υ_x        (27)

其中,υ_x代表x轴子***测量噪声向量，为零均值高斯白噪声向量。

同理，对于***(19)，定义y轴子***状态变量X_y＝[y v_y a_y]^T，则得到如下状态方程

${\stackrel{\&CenterDot;}{X}}_y=F_y{X}_y+G_y{w}_y---\left(28\right)$

其中，

$F_y=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -{\mathrm{\&lambda;}}_y\end{array}\right],G_y=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]---\left(29\right)$

对式(28)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_y(k+1)＝Φ_yX_y(k)+ω_y(k)      (30)

其中，

${\mathrm{\&Phi;}}_y=\exp \left(F_y\mathrm{\&Delta;t}\right)=\left[\begin{array}{ccc}1& \mathrm{\&Delta;t}& 1/{\mathrm{\&lambda;}}_y^2(e^{-{\mathrm{\&lambda;}}_y\mathrm{\&Delta;t}}+{\mathrm{\&lambda;}}_y\mathrm{\&Delta;t}-1)\\ 0& 1& 1/{\mathrm{\&lambda;}}_y(1-e^{-{\mathrm{\&lambda;}}_y\mathrm{\&Delta;t}})\\ 0& 0& e^{-{\mathrm{\&lambda;}}_y\mathrm{\&Delta;t}}\end{array}\right]---\left(31\right)$

ω_y(k)代表y轴子***状态噪声向量，为零均值高斯白噪声向量。

测量方程为

η_y＝y+υ_y        (32)

其中,υ_y代表y轴子***测量噪声向量，为零均值高斯白噪声向量。

对于***(20)，定义z轴子***状态变量X_z＝[z v_z a_z]^T，则得到如下状态方程

${\stackrel{\&CenterDot;}{X}}_z=F_z{X}_z+G_z{w}_z---\left(33\right)$

其中，

$F_z=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -{\mathrm{\&lambda;}}_z\end{array}\right],G_z=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]---\left(34\right)$

对式(28)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_z(k+1)＝Φ_zX_z(k)+ω_z(k)       (35)

其中，

${\mathrm{\&Phi;}}_z=\exp \left(F_z\mathrm{\&Delta;t}\right)=\left[\begin{array}{ccc}1& \mathrm{\&Delta;t}& 1/{\mathrm{\&lambda;}}_z^2(e^{-{\mathrm{\&lambda;}}_z\mathrm{\&Delta;t}}+{\mathrm{\&lambda;}}_z\mathrm{\&Delta;t}-1)\\ 0& 1& 1/{\mathrm{\&lambda;}}_z(1-e^{-{\mathrm{\&lambda;}}_z\mathrm{\&Delta;t}})\\ 0& 0& e^{-{\mathrm{\&lambda;}}_z\mathrm{\&Delta;t}}\end{array}\right]---\left(36\right)$

ω_z(k)代表z轴子***状态噪声向量，为零均值高斯白噪声向量。

测量方程为

η_z＝z+υ_z           (37)

其中,υ_z代表z轴子***测量噪声向量，为零均值高斯白噪声向量。

把探测点惯性坐标系三个轴上的子***，即式(23)和(27)构成的x轴子***，式(28)和(32)构成的y轴子***，以及式(33)和(37)构成的z轴子***，分别设计Kalman滤波器，其预报方程为公式(4),即

$\left\{\begin{array}{c}{\stackrel{\&OverBar;}{X}}_i(k+1)={\mathrm{\&Phi;}}_i{\hat{X}}_i\left(k\right)\\ P_i(k+1/k)={\mathrm{\&Phi;}}_i{P}_i\left(k\right){\mathrm{\&Phi;}}_i^T+Q_i\end{array}\right.,i=x,y,z$

上式中，代表滤波器i第k步的状态估计值，代表滤波器i从第k步到第k+1步的状态预报值，P_i(k)代表滤波器i第k步的状态估计误差协方差矩阵，P_i(k+1/k)代表滤波器i从第k步到第k+1步的状态预报误差协方差矩阵，Q_i为滤波器i模型预报误差协方差矩阵,

滤波器的测量修正方程为公式(6),即

$\left\{\begin{array}{c}{K}_i(k+1)=P_i(k+1/k)H_i^T{[J_i{P}_i(k+1/k)H_i^T+R_i]}^{-1}\\ {\hat{X}}_i(k+1)={\stackrel{\&OverBar;}{X}}_i(k+1)+K_i(k+1)[{\mathrm{\&eta;}}_i(k+1)-H_i{\stackrel{\&OverBar;}{X}}_i(k+1)],i=x,y,z\\ P_i(k+1)=[I-K_i(k+1)H_i]P_i(k+1/k)\end{array}\right.$

其中，上标^T代表矩阵转置运算，上标^-1代表矩阵求逆运算，滤波器i的测量矩阵为式(7),即

H_i＝[1 0 0],i＝x,y,z

R_i为滤波器i测量误差协方差矩阵，I代表单位矩，K_i(k+1)代表滤波器i第k+1步的滤波增益矩阵，P_i(k+1)代表滤波器i第k+1步的状态估计误差协方差矩阵,η_i(k+1)代表滤波器i第k+1步的测量信息，即η_x(k+1)＝x(k+1)+υ_x(k+1)，η_y(k+1)＝y(k+1)+υ_y(k+1),η_z(k+1)＝z(k+1)+υ_z(k+1)。

当没有地面目标跟踪装置测量信息时，只运行预报方程(4)，有地面目标跟踪装置测量信息时，则同时运行预报方程(4)和测量修正方程(6)。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：当高超声速滑翔弹头气动力产生的机动加速度已经充分表现出来，跟踪滤波器已经进入稳态，充分地估计出来目标的机动加速度之后，转入弹道预报算法。

由图1可知，步骤(三一)具体如下：

首先，根据三个子滤波器的估计结果得到跟踪结束时刻kT时空间高超声速滑翔弹头在探测点惯性系中的位置[x y z]^T，探测点惯性坐标系中的位置与地心惯性系中的位置关系为

[x_I y_I z_I]^T＝[x y z]^T+[0 R_e 0]^T      (38)

其中，R_e＝6378km为地球平均半径；

根据跟踪结束时刻的高超声速滑翔弹头的加速度估计值[a_x a_y a_z]^T，由于探测点惯性坐标系的三个与地心惯性系的三个轴平行，可得

[a_xI a_yI a_zI]^T＝[a_x a_y a_z]^T      (39)

其中，a_xI、a_yI和a_zI代表高超声速滑翔弹头的加速度在地心惯性系中的投影

用下式计算出当前高超声速滑翔弹头受到的气动力：

$R=m(\left[\begin{array}{c}{a}_{\mathrm{xI}}\\ a_{\mathrm{yI}}\\ a_{\mathrm{zI}}\end{array}\right]-\left[\begin{array}{c}{g}_{\mathrm{xI}}\\ g_{\mathrm{yI}}\\ g_{\mathrm{zI}}\end{array}\right])---\left(40\right)$

在15Mach左右的时候，升阻比大约为3:1，由于R_z1＝0，可得

$R_{y1}=\frac{3}{\sqrt{10}}\left|\right|R\left|\right|---\left(41\right)$

根据跟踪结束时刻的高超声速滑翔弹头的速度估计值[v_x v_y v_z]^T，由于探测点惯性坐标系的三个与地心惯性系的三个轴平行，可得

[v_xI v_yI v_zI]^T＝[v_x v_y v_z]^T     (42)

根据跟踪结束时高超声速滑翔弹头的速度[v_xI v_yI v_zI]^T，用式(12)和(13)计算其弹道角θ和ψ_v。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤(三二)具体如下：

根据跟踪结束时高超声速滑翔弹头的速度[v_xI v_yI v_zI]^T，用式(11)计算出其速率v，再利用式(9)求得当前的动压头q；

用下式计算出当前的升力系数为

$c_y=\frac{R_{y1}}{\mathrm{qS}}---\left(43\right)$

攻角为10°时，升力系数为0.37，所以用下式估算当前的攻角：

$\mathrm{\&alpha;}=\frac{c_y}{0.37}\frac{10}{57.3}---\left(44\right)$

高超声速滑翔弹头采用倾斜转弯控制方式，即侧滑角β＝0，然后由式(14)中的前两式可以计算出俯仰角和偏航角ψ；

估算出当前的滚转角γ：由式(40)已经计算出气动力在惯性系的投影，把它投影到目标体坐标系，即

$\left[\begin{array}{c}{R}_{x1}\\ R_{y1}\\ R_{z1}\end{array}\right]=C_{t0\&RightArrow;1}{C}_{0\&RightArrow;t0}R$

上式进一步可写作

令

则根据式(45)可得

$\left[\begin{array}{cc}{b}_3& b_2\\ {-b}_2& b_3\end{array}\right]\left[\begin{array}{c}\sin \mathrm{\&gamma;}\\ \cos \mathrm{\&gamma;}\end{array}\right]=\left[\begin{array}{c}{R}_{y1}\\ R_{z1}\end{array}\right]---\left(46\right)$

其中，R_y1由式(41)估算出来，R_z1＝0，用式(46)求得sinγ和cosγ，然后用下式唯一地确定当前的滚转角

γ＝tan2^-1(sinγ,cosγ)      (47)

其它步骤及参数与具体实施方式一至四之一相同。

仿真实验：

针对临近空间高超声速滑翔弹头，本发明考虑了目标受气动力等复杂情况，提出了一种新型跳跃弹道预测方法，相比于传统方法提高了弹道预测精度，减小了弹道预测误差，填补了临近空间高超声速飞行器弹道预报领域的空白。下面通过数值仿真验证本发明的效果。

设仿真初始时刻t₀＝82s，高超声速滑翔弹头在探测点惯性坐标系中的初始位置为x_I(0)＝1161km，y_I(0)＝-76km，z_I(0)＝-1.1km，初始速度在探测点惯性坐标系中的分量为v_xI(0)＝-3024m/s，v_yI(0)＝558m/s，v_zI(0)＝0m/s。高超声速滑翔弹头保持等攻角和等滚转角飞行，攻角为α＝10°，滚转角为γ＝45°。仿真结束时刻为t_f＝400s。

高超声速滑翔弹头跟踪Kalman滤波器的初始估计值取为

X_x(0)＝[x_I(0)+300 v_xI(0)+20 0]^T，X_y(0)＝[y_I(0)+300 v_yI(0)+20 0]^T，

X_z(0)＝[z_I(0)+300 v_zI(0)+20 0]^T

状态估计方差阵的初值取为

$P_x\left(0\right)=P_y\left(0\right)=P_z\left(0\right)=\left[\begin{array}{ccc}10& 0& 0\\ 0& 10& 0\\ 0& 0& 10\end{array}\right]$

地面测距误差取为300m，测角误差取为俯仰方向0.5mrad，偏航方向0.3mrad。滤波器一步预报周期取为0.01s，滤波器测量修正周期即地面测量信息的数据更新率为0.2s。在前118s利用测量信息进行修正。后200s进行弹道预报。

以一次仿真结果为例，高超声速滑翔弹头的飞行弹道如图2所示，加速度估计结果如图3—图5所示，速度估计结果如图6—图8所示，位置估计如图9—图11所示，总位置误差如图12所示。

在上述条件下，经过100次Monte-Carlo仿真得到的高超声速滑翔弹头最终位置预报的统计结果如图13所示，统计计算得出终端位置预报误差小于50km的概率为95％。

Claims (5)

1.一种临近空间高超声速滑翔弹头跳跃弹道预测方法，其特征在于它按以下步骤实现：

一、高超声速滑翔弹头的弹道方程的建立；

建立临近空间高超声速滑翔弹头的弹道方程，包括气动力计算方程、加速度计算方程、速度和位置计算方程、速率和弹道角计算方程，以及姿态角计算方程；

综合计算方程，得到高超声速滑翔弹头质心运动速度和位置方程如式(1)～(3)所示：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_I=v_{\mathrm{xI}}\\ {\stackrel{\&CenterDot;}{v}}_{\mathrm{xI}}=-\frac{\mathrm{\&mu;}{x}_I}{r_I^3}+a_{\mathrm{axI}}\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}\stackrel{}{{\stackrel{\&CenterDot;}{y}}_I}{=v}_{\mathrm{yI}}\\ {\stackrel{\&CenterDot;}{v}}_{\mathrm{yI}}=-\frac{\mathrm{\&mu;}{y}_I}{r_I^3}+a_{\mathrm{ayI}}\end{array}\right.---\left(2\right)$

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{z}}_I=v_{\mathrm{zI}}\\ {\stackrel{\&CenterDot;}{v}}_{\mathrm{zI}}=-\frac{\mathrm{\&mu;}{z}_I}{r_I^3}+a_{\mathrm{azI}}\end{array}\right.---\left(3\right)$

其中，x_I、y_I和z_I代表高超声速滑翔弹头在地心惯性坐标系中的位置，代表高超声速滑翔弹头的位置到地球中心的距离,和则代表高超声速滑翔弹头在地心惯性坐标系中的速度，a_axI、a_ayI和a_azI代表气动力产生的高超声速滑翔弹头机动加速度在地心惯性坐标系中的分量，μ为引力参数；

二、高超声速滑翔弹头的跟踪滤波器设计；

设计高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器，实时地跟踪并估计高超声速滑翔弹头的位置、速度和加速度；

其中，所述高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器的状态一步预报方程为

$\left\{\begin{array}{c}{\stackrel{\&OverBar;}{X}}_i(k+1)={\mathrm{\&Phi;}}_i{\hat{X}}_i\left(k\right)\\ P_i(k+1/k)={\mathrm{\&Phi;}}_i{P}_i\left(k\right){\mathrm{\&Phi;}}_i^T+Q_i\end{array}\right.,i=x,y,z---\left(4\right)$

上式中，代表滤波器i第k步的状态估计值，代表滤波器i从第k步到第k+1步的状态预报值，P_i(k)代表滤波器i第k步的状态估计误差协方差矩阵，P_i(k+1/k)代表滤波器i从第k步到第k+1步的状态预报误差协方差矩阵，Q_i为滤波器i模型预报误差协方差矩阵,

${\mathrm{\&Phi;}}_i=\exp \left(F_i\mathrm{\&Delta;t}\right)=\left[\begin{array}{ccc}1& \mathrm{\&Delta;t}& 1/{\mathrm{\&lambda;}}_i^2(e^{-{\mathrm{\&lambda;}}_{i}t}+{\mathrm{\&lambda;}}_i\mathrm{\&Delta;t}-1)\\ 0& 1& 1/{\mathrm{\&lambda;}}_{i}1-e^{-{\mathrm{\&lambda;}}_i\mathrm{\&Delta;t}}\\ 0& 0& e^{-{\mathrm{\&lambda;}}_i\mathrm{\&Delta;t}}\end{array}\right]---\left(5\right)$

所述运动轨迹跟踪卡尔曼滤波器的测量修正方程为

$\left\{\begin{array}{c}{K}_i(k+1)=P_i(k+1/k)H_i^T[H_i{P}_i(k+1/k)H_i^T+R_i{]}^{-1}\\ {\hat{X}}_i(k+1)={\stackrel{\&OverBar;}{X}}_i(k+1)+K_i(k+1)[{\mathrm{\&eta;}}_i(k+1)-H_i{\stackrel{\&OverBar;}{X}}_i(k+1)]\\ P_i(k+1)=[I-K_i(k+1)H_i]P_i(k+1/k)\end{array}\right.,i=x,y,z---\left(6\right)$

其中，上标^T代表矩阵转置运算，上标^-1代表矩阵求逆运算，滤波器i的测量矩阵

H_i＝[1 0 0],i＝x,y,z           (7)

R_i为滤波器i测量误差协方差矩阵，I代表单位矩阵，K_i(k+1)代表滤波器i第k+1步的滤波增益矩阵，η_i(k+1)代表滤波器i第k+1步的测量信息，P_i(k+1)代表滤波器i第k+1步的状态估计误差协方差矩阵；

三、高超声速滑翔弹头的弹道预测；

(三一)根据步骤二中运动轨迹跟踪卡尔曼滤波器跟踪结束kT时刻的高超声速滑翔弹头的加速度，其中，k为大于0的整数，T为计算周期，一般可取T＝0.01s，应用步骤一中的临近空间高超声速滑翔弹头的加速度计算方程求出跟踪结束kT时刻临近空间高超声速滑翔弹头的气动力产生的加速度，并根据跟踪结束kT时刻高超声速滑翔弹头的速度，应用步骤一中的临近空间高超声速滑翔弹头的弹道角计算方程计算其弹道角；

(三二)根据计算出来的kT时刻临近空间高超声速滑翔弹头的气动力产生的加速度和弹道角,结合步骤一中的临近空间高超声速滑翔弹头的气动力计算方程和姿态角计算方程估算跟踪结束kT时刻高超声速滑翔弹头飞行时的攻角和姿态角；

(三三)设在随后的kT到(k+1)T，(k+2)T,...,(k+N)T时间内，N为大于2的正整数，临近空间高超声速滑翔弹头进行等攻角和等滚转角飞行，则纵向平面内呈现出上下跳跃式的飞行弹道；同时，高超声速滑翔弹头采用BTT控制方式，侧滑角β等于零，通过控制滚转角在其体系侧向输出一个机动加速度进行转弯，即滚转角指令γ_c保持常值；

(三四)根据等攻角、零侧滑角和等滚转角飞行条件，计算kT时刻气动力产生的加速度、kT时刻的位置x_I、y_I和z_I、kT时刻的速度v_xI、v_yI和v_zI代入式(1)～(3)，经过一步数值积分运算后得到高超声速滑翔弹头在(k+1)T时刻在地心惯性坐标系中的速度和位置；

计算高超声速滑翔弹头在(k+1)T时刻受到的气动力，计算高超声速滑翔弹头在(k+1)T时刻受到的气动力加速度，再把(k+1)T时刻的气动力加速度以及(k+1)T时刻的位置和速度代入式(1)～(3),计算出(k+2)T时刻的位置和速度，以此类推，计算出直到(k+N)T时刻的位置和速度。

2.根据权利要求1所述的一种临近空间高超声速滑翔弹头跳跃弹道预测方法，其特征在于步骤一具体为：

(一)计算高超声速滑翔弹头气动力

在滑翔阶段，高超声速滑翔弹头受到的气动力在体系内的投影按照下式计算：

$\left\{\begin{array}{c}{R}_{x1}=-c_x\mathrm{qS}\\ R_{y1}=c_y\mathrm{sign}\left(\mathrm{\&alpha;}\right)\mathrm{qS}\\ R_{z1}=0\end{array}\right.---\left(8\right)$

其中，α为攻角，

$\mathrm{sign}\left(\mathrm{\&alpha;}\right)=\left\{\begin{array}{c}1,\mathrm{\&alpha;}>0\\ -1,\mathrm{\&alpha;}<0\end{array}\right.$

c_y为升力系数，c_x为阻力系数，S为特征参考面积；

q＝ρv²/2        (9)

其中，q为动压头，v为高超声速滑翔弹头的速率，ρ为空气密度，它随着飞行高度h变化，可查标准的大气密度表得到；高超声速滑翔弹头飞行高度h用下式计算:

$h=\sqrt{x_I^2+y_I^2+z_I^2}-R_e---\left(10\right)$

其中，x_I、y_I和z_I代表高超声速滑翔弹头在地心惯性坐标系中的位置，R_e＝6378km为地球平均半径；

设在飞行中侧滑角β＝0，因此气动力侧向力R_z1＝0；

(二)计算高超声速滑翔弹头的速率

高超声速滑翔弹头的速率为

$v=\sqrt{v_{\mathrm{xI}}^2+v_{\mathrm{yI}}^2+v_{\mathrm{zI}}^2}---\left(11\right)$

其中，v_xI、v_yI和v_zI代表高超声速滑翔弹头在地心惯性坐标系中的速度；

将高超声速滑翔弹头的速度矢量投影到目标惯性坐标系v_xt0、v_yt0与v_zt0，即

$\left[\begin{array}{c}{v}_{\mathrm{xt}0}\\ v_{\mathrm{yt}0}\\ v_{\mathrm{zt}0}\end{array}\right]=C_{0\&RightArrow;t0}\left[\begin{array}{c}{v}_{\mathrm{xI}}\\ v_{\mathrm{yI}}\\ v_{\mathrm{zI}}\end{array}\right]---\left(12\right)$

其中，

$C_{0\&RightArrow;t0}=\left[\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}\right]$

(三)计算高超声速滑翔弹头的弹道仰角θ和弹道偏角ψ_v：

$\left\{\begin{array}{c}\mathrm{\&theta;}={\sin }^{-1}(v_{\mathrm{yt}0}/v)\\ {\mathrm{\&psi;}}_v=-{\tan }^{-1}(v_{\mathrm{zt}0}/v_{\mathrm{xt}0})\end{array}\right.---\left(13\right)$

高超声速滑翔弹头的姿态角ψ和γ按照下式计算：

其中，为俯仰角，ψ为偏航角，γ为滚转角，γ_c为滚转角指令；

(四)计算高超声速滑翔弹头机动加速度在地心惯性坐标系中的投影

$\left[\begin{array}{c}{a}_{\mathrm{axI}}\\ a_{\mathrm{ayI}}\\ a_{\mathrm{azI}}\end{array}\right]=C_{0\&RightArrow;t0}^T{C}_{t0\&RightArrow;1}^T\left[\begin{array}{c}{R}_{x1}/m\\ R_{y1}/m\\ R_{z1}\text{/m}\end{array}\right]---\left(15\right)$

其中，

m为高超声速滑翔弹头质量；

(五)计算高超声速滑翔弹头的速度和位置

高超声速滑翔弹头质心运动速度和位置方程如式(1)～(3)所示,其中包含了高超声速滑翔弹头受到的重力加速度，即

$\left\{\begin{array}{c}{g}_{\mathrm{xI}}=-\frac{\mathrm{\&mu;}{x}_I}{r_I^3}\\ g_{\mathrm{yI}}=-\frac{\mathrm{\&mu;}{y}_I}{r_I^3}\\ g_{\mathrm{zI}}=-\frac{\mathrm{\&mu;}{z}_I}{r_I^3}\end{array}\right.---\left(16\right)$

其中，g_xI、g_yI与g_zI代表重力加速度在地心惯性系中的投影。

3.根据权利要求2所述的一种临近空间高超声速滑翔弹头跳跃弹道预测方法，其特征在于步骤二具体为：

(一)在探测点惯性坐标系的三个轴上，均用Singer模型来分别描述高超声速滑翔弹头的加速度分量的变化，即

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{a}}_x=-{\mathrm{\&lambda;}}_x{a}_x+w_x\\ {\stackrel{\&CenterDot;}{a}}_y=-{\mathrm{\&lambda;}}_y{a}_y+w_y\\ {\stackrel{\&CenterDot;}{a}}_z=-{\mathrm{\&lambda;}}_z{a}_z+w_z\end{array}\right.---\left(17\right)$

其中，λ_x、λ_y和λ_z分别代表在探测点惯性坐标系下o_Fx轴、o_Fy轴和o_Fz轴方向上目标机动时间常数的倒数，w_x、w_y和w_z分别代表探测点惯性坐标系下o_Fx轴、o_Fy轴和o_Fz轴方向上的零均值高斯白噪声；加速度a_x、a_y和a_z中既包含气动力产生的加速度a_axI、a_ayI和a_azI，又包含重力产生的加速度g_xI、g_yI和g_zI；探测点惯性坐标系的o_Fx轴、o_Fy轴和o_Fz轴与地心惯性系的三个轴o_Ix_I、o_Iy_I和o_Iz_I分别平行，气动力产生的加速度在地心惯性系中的投影a_axI、a_ayI和a_azI就是此加速度在探测点惯性系中的投影；

(二)在探测点惯性坐标系下建立起如下三组运动方程

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=v_x\\ {\stackrel{\&CenterDot;}{v}}_x=a_x\\ {\stackrel{\&CenterDot;}{a}}_x=-{\mathrm{\&lambda;}}_x{a}_x+w_x\end{array}\right.---\left(18\right)$

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{y}=v_y\\ {\stackrel{\&CenterDot;}{v}}_y=a_y\\ {\stackrel{\&CenterDot;}{a}}_y=-{\mathrm{\&lambda;}}_y{a}_y+w_y\end{array}\right.---\left(19\right)$

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{z}=v_z\\ {\stackrel{\&CenterDot;}{v}}_z=a_z\\ {\stackrel{\&CenterDot;}{a}}_z=-{\mathrm{\&lambda;}}_z{a}_z+w_z\end{array}\right.---\left(20\right)$

其中，x、y和z代表高超声速滑翔弹头在探测点惯性坐标系中的位置，和则代表高超声速滑翔弹头速度在探测点惯性坐标系中的投影。因为探测点惯性坐标系的o_Fx轴、o_Fy轴和o_Fz轴与地心惯性系的三个轴o_Ix_I、o_Iy_I和o_Iz_I分别平行，所以

$\left\{\begin{array}{c}{v}_x=v_{\mathrm{xI}}\\ v_y=v_{\mathrm{yI}}\\ v_z=v_{\mathrm{zI}}\end{array}\right.---\left(21\right)$

而由于探测点惯性坐标系只是相对于地心惯性系探测点惯性坐标系在o_Iy_I轴方向上作了平移，所以有

$\left\{\begin{array}{c}x=x_I\\ y=y_I-R_e\\ z=z_I\end{array}\right.---\left(22\right)$

对于***(18)，定义x轴子***状态变量X_x＝[x v_x a_x]^T，则得到如下状态方程

${\stackrel{\&CenterDot;}{X}}_x=F_x{X}_x+G_x{w}_x---\left(23\right)$

其中，

$F_x=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -{\mathrm{\&lambda;}}_x\end{array}\right],G_x=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]---\left(24\right)$

对式(23)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_x(k+1)＝Φ_xX_x(k)+ω_x(k)    (25)

其中，

${\mathrm{\&Phi;}}_x=\exp \left(F_x\mathrm{\&Delta;t}\right)=\left[\begin{array}{ccc}1& \mathrm{\&Delta;t}& 1/{\mathrm{\&lambda;}}_x^2(e^{-{\mathrm{\&lambda;}}_{x}t}+{\mathrm{\&lambda;}}_x\mathrm{\&Delta;t}-1)\\ 0& 1& 1/{\mathrm{\&lambda;}}_{x}1-e^{-{\mathrm{\&lambda;}}_x\mathrm{\&Delta;t}}\\ 0& 0& e^{-{\mathrm{\&lambda;}}_x\mathrm{\&Delta;t}}\end{array}\right]---\left(26\right)$

ω_x(k)代表x轴子***状态噪声向量，为零均值高斯白噪声向量。

高超声速滑翔弹头的位置由地面目标跟踪装置测量得到，则测量方程为

η_x＝x+υ_x       (27)

其中,υ_x代表x轴子***测量噪声向量，为零均值高斯白噪声向量。

同理，对于***(19)，定义y轴子***状态变量X_y＝[y v_y a_y]^T，则得到如下状态方程

${\stackrel{\&CenterDot;}{X}}_y=F_y{X}_y+G_y{w}_y---\left(28\right)$

其中，

$F_y=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -{\mathrm{\&lambda;}}_y\end{array}\right],G_y=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]---\left(29\right)$

对式(28)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_y(k+1)＝Φ_yX_y(k)+ω_y(k)    (30)

其中，

${\mathrm{\&Phi;}}_y=\exp \left(F_y\mathrm{\&Delta;t}\right)=\left[\begin{array}{ccc}1& \mathrm{\&Delta;t}& 1/{\mathrm{\&lambda;}}_y^2(e^{-{\mathrm{\&lambda;}}_{y}t}+{\mathrm{\&lambda;}}_y\mathrm{\&Delta;t}-1)\\ 0& 1& 1/{\mathrm{\&lambda;}}_{x}1-e^{-{\mathrm{\&lambda;}}_y\mathrm{\&Delta;t}}\\ 0& 0& e^{-{\mathrm{\&lambda;}}_y\mathrm{\&Delta;t}}\end{array}\right]---\left(31\right)$

ω_y(k)代表y轴子***状态噪声向量，为零均值高斯白噪声向量。

测量方程为

η_y＝y+υ_y   (32)

其中,υ_y代表y轴子***测量噪声向量，为零均值高斯白噪声向量。

对于***(20)，定义z轴子***状态变量X_z＝[z v_z a_z]^T，则得到如下状态方程

${\stackrel{\&CenterDot;}{X}}_z=F_z{X}_z+G_z{w}_z---\left(33\right)$

其中，

$F_z=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -{\mathrm{\&lambda;}}_z\end{array}\right],G_z=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]---\left(34\right)$

对式(28)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_z(k+1)＝Φ_zX_z(k)+ω_z(k)        (35)

其中，

${\mathrm{\&Phi;}}_z=\exp \left(F_z\mathrm{\&Delta;t}\right)=\left[\begin{array}{ccc}1& \mathrm{\&Delta;t}& 1/{\mathrm{\&lambda;}}_z^2(e^{-{\mathrm{\&lambda;}}_{z}t}+{\mathrm{\&lambda;}}_z\mathrm{\&Delta;t}-1)\\ 0& 1& 1/{\mathrm{\&lambda;}}_{z}1-e^{-{\mathrm{\&lambda;}}_z\mathrm{\&Delta;t}}\\ 0& 0& e^{-{\mathrm{\&lambda;}}_z\mathrm{\&Delta;t}}\end{array}\right]---\left(36\right)$

ω_z(k)代表z轴子***状态噪声向量，为零均值高斯白噪声向量。

测量方程为

η_z＝z+υ_z         (37)

其中,υ_z代表z轴子***测量噪声向量，为零均值高斯白噪声向量。

把探测点惯性坐标系三个轴上的子***，即式(23)和(27)构成的x轴子***，式(28)和(32)构成的y轴子***，以及式(33)和(37)构成的z轴子***，分别设计Kalman滤波器，其预报方程为公式(4),即

$\left\{\begin{array}{c}{\stackrel{\&OverBar;}{X}}_i(k+1)={\mathrm{\&Phi;}}_i{\hat{X}}_i\left(k\right)\\ P_i(k+1/k)={\mathrm{\&Phi;}}_i{P}_i\left(k\right){\mathrm{\&Phi;}}_i^T+Q_i\end{array}\right.,i=x,y,z$

上式中，代表滤波器i第k步的状态估计值，代表滤波器i从第k步到第k+1步的状态预报值，P_i(k)代表滤波器i第k步的状态估计误差协方差矩阵，P_i(k+1/k)代表滤波器i从第k步到第k+1步的状态预报误差协方差矩阵，Q_i为滤波器i模型预报误差协方差矩阵,

滤波器的测量修正方程为公式(6),即

$\left\{\begin{array}{c}{K}_i(k+1)=P_i(k+1/k)H_i^T[H_i{P}_i(k+1/k)H_i^T+R_i{]}^{-1}\\ {\hat{X}}_i(k+1)={\stackrel{\&OverBar;}{X}}_i(k+1)+K_i(k+1)[{\mathrm{\&eta;}}_i(k+1)-H_i{\stackrel{\&OverBar;}{X}}_i(k+1)]\\ P_i(k+1)=[I-K_i(k+1)H_i]P_i(k+1/k)\end{array}\right.,i=x,y,z$

其中，上标^T代表矩阵转置运算，上标^-1代表矩阵求逆运算，滤波器i的测量矩阵为式(7),即

H_i＝[1 0 0],i＝x,y,z

R_i为滤波器i测量误差协方差矩阵，I代表单位矩，K_i(k+1)代表滤波器i第k+1步的滤波增益矩阵，P_i(k+1)代表滤波器i第k+1步的状态估计误差协方差矩阵,η_i(k+1)代表滤波器i第k+1步的测量信息，即η_x(k+1)＝x(k+1)+υ_x(k+1)，η_y(k+1)＝y(k+1)+υ_y(k+1),η_z(k+1)＝z(k+1)+υ_z(k+1)。

当没有地面目标跟踪装置测量信息时，只运行预报方程(4)，有地面目标跟踪装置测量信息时，则同时运行预报方程(4)和测量修正方程(6)。

4.根据权利要求3所述的一种临近空间高超声速滑翔弹头跳跃弹道预测方法，其特征在于步骤(三一)具体如下：

首先，根据三个子滤波器的估计结果得到跟踪结束时刻kT时空间高超声速滑翔弹头在探测点惯性系中的位置[x y z]^T，探测点惯性坐标系中的位置与地心惯性系中的位置关系为

[x_I y_I z_I]^T＝[x y z]^T+[0 R_e 0]^T      (38)

其中，R_e＝6378km为地球平均半径；

根据跟踪结束时刻的高超声速滑翔弹头的加速度估计值[a_x a_y a_z]^T，由于探测点惯性坐标系的三个与地心惯性系的三个轴平行，可得

[a_xI a_yI a_zI]^T＝[a_x a_y a_z]^T     (39)

其中，a_xI、a_yI和a_zI代表高超声速滑翔弹头的加速度在地心惯性系中的投影

用下式计算出当前高超声速滑翔弹头受到的气动力：

$R=m(\left[\begin{array}{c}{a}_{\mathrm{xI}}\\ a_{\mathrm{yI}}\\ a_{\mathrm{zI}}\end{array}\right]-\left[\begin{array}{c}{g}_{\mathrm{xI}}\\ g_{\mathrm{yI}}\\ g_{\mathrm{zI}}\end{array}\right])---\left(40\right)$

在15Mach左右的时候，升阻比大约为3:1，由于R_z1＝0，可得

$R_{y1}=\frac{3}{\sqrt{10}}\left|\right|R\left|\right|---\left(41\right)$

根据跟踪结束时刻的高超声速滑翔弹头的速度估计值[v_x v_y v_z]^T，由于探测点惯性坐标系的三个与地心惯性系的三个轴平行，可得

[v_xI v_yI v_zI]^T＝[v_x v_y v_z] ^T    (42)

根据跟踪结束时高超声速滑翔弹头的速度[v_xI v_yI v_zI]^T，用式(12)和(13)计算其弹道角θ和ψ_v。

5.根据权利要求4所述的一种临近空间高超声速滑翔弹头跳跃弹道预测方法，其特征在于步骤(三二)具体如下：

根据跟踪结束时高超声速滑翔弹头的速度[v_xI v_yI v_zI]^T，用式(11)计算出其速率v，再利用式(9)求得当前的动压头q；

用下式计算出当前的升力系数为

$c_y=\frac{R_{y1}}{\mathrm{qS}}---\left(43\right)$

攻角为10°时，升力系数为0.37，所以用下式估算当前的攻角：

$\mathrm{\&alpha;}=\frac{c_y}{0.37}\frac{10}{57.3}---\left(44\right)$

高超声速滑翔弹头采用倾斜转弯控制方式，即侧滑角β＝0，然后由式(14)中的前两式可以计算出俯仰角和偏航角ψ；

估算出当前的滚转角γ：由式(40)已经计算出气动力在惯性系的投影，把它投影到目标体坐标系，即

$\left[\begin{array}{c}{R}_{x1}\\ R_{y1}\\ R_{z1}\end{array}\right]=C_{t0\&RightArrow;1}{C}_{0\&RightArrow;t0}R$

上式进一步可写作

令

则根据式(45)可得

$\left[\begin{array}{cc}{b}_3& b_2\\ -b_2& b_3\end{array}\right]\left[\begin{array}{c}\sin \mathrm{\&gamma;}\\ \cos \mathrm{\&gamma;}\end{array}\right]=\left[\begin{array}{c}{R}_{y1}\\ R_{z1}\end{array}\right]---\left(46\right)$

其中，R_y1由式(41)估算出来，R_z1＝0，用式(46)求得sinγ和cosγ，用下式确定当前的滚转角

γ＝tan2^-1(sinγ,cosγ)       (47)。