基于机动模式识别的临近空间高超声速飞行器弹道预测方法

1.本发明属于飞行器弹道预测技术领域，具体涉及一种基于机动模式识别的临近空间高超声速飞行器弹道预测方法。


背景技术：

2.临近空间高超声速飞行器是指在临近空间中以不小于5ma的飞行速度做巡航或滑翔飞行的一类空间飞行器，它们响应快速、突防能力超强，是一种具有战略威慑能力的武器。高超声速飞行器助推到预定高度后，高超声速滑翔体或高超声速巡航体与助推器分离，继续攻击预定目标，末段攻击速度可能达马赫数5以上。在临近空间远程滑翔和机动飞行能力较强。
3.现有的针对临近空间飞行器的弹道预测方法是使用singer模型通过卡尔曼滤波进行飞行器状态估计，得到飞行器位置、速度以及加速度的信息，且假设临近空间高超声速飞行器保持姿态角恒定进行后续的弹道预测。但是现有方法并未考虑不同的机动类型对飞行器弹道的影响，导致采用现有方法对飞行器弹道进行预测时，所获得的弹道预测结果的误差较大。


技术实现要素：

4.本发明的目的是为解决现有方法并未考虑不同的机动类型对飞行器弹道的影响，导致采用现有方法对飞行器弹道进行预测时获得的弹道预测结果的误差大的问题，而提出的一种基于机动模式识别的临近空间高超声速飞行器弹道预测方法。
5.本发明为解决上述技术问题所采取的技术方案是：
6.一种基于机动模式识别的临近空间高超声速飞行器弹道预测方法，所述方法在跟踪阶段，利用卡尔曼滤波器对高超声速飞行器的状态信息进行跟踪，将跟踪阶段结束时刻的状态估计值记为：
[0007][0008]
其中，kt代表跟踪阶段结束时刻，x
p
(kt)、y
p
(kt)和z
p
(kt)为跟踪结束kt时刻的高超声速飞行器在观测惯性坐标系的x轴、y轴和z轴向的位置，v
xp
(kt)、v
yp
(kt)和v
zp
(kt)分别为跟踪结束kt时刻的高超声速飞行器在观测惯性坐标系的x轴、y轴和z轴向的速度；和分别为跟踪结束kt时刻对z
x
(kt)、zy(kt)和zz(kt)的估计值，z
x
(kt)、zy(kt)和zz(kt)为跟踪结束kt时刻的状态变量，将x
p
(kt)、y
p
(kt)、z
p
(kt)、v
xp
(kt)、v
yp
(kt)和v
zp
(kt)作为弹道预测的初值；
[0009]
根据跟踪过程中状态变量z
x
、zy和zz的变化特点判断高超声速飞行器的机动类型；
[0010]
若机动类型为正弦加速度机动，则使用模型f(t)＝a1sin(a2t+a3)分别对z
x
、zy、zz进行拟合，得到函数z
x
(t)、zy(t)、zz(t)，a1、a2和a3为拟合系数；
[0011]
若机动类型为常值加速度机动，则使用模型f(t)＝a4t+a5分别对z
x
、zy、zz进行拟合得到函数z
x
(t)、zy(t)、zz(t)，a4和a5为拟合系数；
[0012]
根据拟合得到的函数z
x
(t)、zy(t)和zz(t)以及弹道预测的初值，进行临近空间高超声速飞行器的弹道预测。
[0013]
进一步地，所述观测惯性坐标系的定义为：
[0014]
观测惯性坐标系以观测基站为原点o，以地心到观测基站的连线指向天方向为y轴正方向；在观测基站、地心和目标所构成的平面内，由原点指向目标方向为x轴正方向；再根据右手定则确定z轴方向。
[0015]
进一步地，所述根据跟踪过程中状态变量z
x
、zy和zz的变化特点判断高超声速飞行器的机动类型时，所采用的是聚类方法。
[0016]
进一步地，所述根据拟合得到的函数z
x
(t)、zy(t)和zz(t)以及弹道预测的初值，进行临近空间高超声速飞行器的弹道预测；其具体过程为：
[0017]
根据拟合得到的函数z
x
(t)、zy(t)和zz(t)，计算(k+1)t时刻高超声速飞行器在观测惯性坐标系的x轴、y轴和z轴向的加速度a
(k+1)t,x
、a
(k+1)t,y
和a
(k+1)t,z
：
[0018][0019]
其中，ρ为空气密度，v
t(k+1)t,x
、v
t(k+1)t,y
和v
t(k+1)t,z
分别为(k+1)t时刻高超声速飞行器在观测惯性坐标系的x轴、y轴和z轴向的速度；为由弹道坐标系到观测惯性坐标系的变换矩阵；z
x
((k+1)t)为根据函数z
x
(t)得到的(k+1)t时刻的状态变量，zy((k+1)t)为根据函数zy(t)得到的(k+1)t时刻的状态变量，zz((k+1)t)为根据函数zz(t)得到的(k+1)t时刻的状态变量；
[0020]
根据弹道预测的初值以及a
(k+1)t,x
、a
(k+1)t,y
和a
(k+1)t,z
，得到(k+1)t时刻高超声速飞行器的位置和速度在观测惯性坐标系下的预测值；
[0021]
再根据(k+1)t时刻高超声速飞行器的位置和速度在观测惯性坐标系下的预测值以及(k+2)t时刻高超声速飞行器在观测惯性坐标系的x轴、y轴和z轴向的加速度，得到(k+2)t时刻高超声速飞行器的位置和速度在观测惯性坐标系下的预测值，以此类推，直至获得了每个时刻高超声速飞行器的位置和速度的预测值，完成临近空间高超声速飞行器的弹道预测。
[0022]
进一步地，所述弹道坐标系的定义为：
[0023]
弹道坐标系以观测基站为原点o
′
，以目标速度方向为x
′
轴正方向，y
′
轴在观测基站、地心和目标所构成的平面内，y
′
轴与x
′
轴垂直，且y
′
轴以指向目标方向为正；再根据右手定则确定z
′
轴方向。
[0024]
进一步地，所述根据弹道预测的初值以及a
(k+1)t,x
、a
(k+1)t,y
和a
(k+1)t,z
，得到(k+1)t时刻高超声速飞行器的位置和速度在观测惯性坐标系下的预测值；其具体过程为：
[0025]
建立微分方程如下：
[0026][0027][0028][0029][0030][0031][0032]
其中，μ为引力系数；re为地球半径；v
(k+1)t,x
、v
(k+1)t,y
和v
(k+1)t,z
分别为(k+1)t时刻高超声速飞行器在观测惯性坐标系下x轴、y轴和z轴向的速度；x
(k+1)t
、y
(k+1)t
和z
(k+1)t
分别为(k+1)t时刻高超声速飞行器在观测惯性坐标系下x轴、y轴和z轴向的位置坐标；为x
(k+1)t
的一阶导数，为v
(k+1)t,y
的一阶导数，为z
(k+1)t
的一阶导数，为v
(k+1)t,x
的一阶导数，为v
(k+1)t,y
的一阶导数，为v
(k+1)t,z
的一阶导数；
[0033]
根据弹道预测的初值以及a
(k+1)t,x
、a
(k+1)t,y
和a
(k+1)t,z
求解微分方程，得到(k+1)t时刻高超声速飞行器的位置和速度在观测惯性坐标系下的预测值。
[0034]
更进一步地，所述根据弹道预测的初值以及a
(k+1)t,x
、a
(k+1)t,y
和a
(k+1)t,z
求解微分方程，所采用的是龙格库塔法。
[0035]
本发明的有益效果是：
[0036]
本发明针对临近空间高超声速飞行器，考虑目标受气动力等复杂情况的影响，在飞行器质量、参考面积、推力大小，气动力参数等敌方飞行器参数未知的情况下，基于当前时刻对于位置和速度的预测，通过对观测出来的状态量进行区分，辨识出其机动方式，然后确定数据拟合方案，利用拟合得到的函数，通过求解微分方程，对下一时刻的位置及速度进行预测，直至完成弹道预测。相比于传统方法，本发明方法通过辨识出机动模式提高了弹道预测精度，减小了弹道预测误差。
附图说明
[0037]
目标机动模式为正弦加速度机动情况下的仿真：
[0038]
图1为仿真实验中高超声速飞行器在观测惯性坐标系中纵向平面飞行弹道图；
[0039]
图2为仿真实验中高超声速飞行器在观测惯性坐标系中侧向平面飞行弹道图；
[0040]
图3为仿真实验中高超声速飞行器的x轴的加速度估计误差图；
[0041]
图中：t代表时间，error in a
x
代表x轴的加速度估计误差；
[0042]
图4为仿真实验中高超声速飞行器的y轴的加速度估计误差图；
[0043]
error in ay代表y轴的加速度估计误差；
[0044]
图5为仿真实验中高超声速飞行器的z轴的加速度估计误差图；
[0045]
error in az代表z轴的加速度估计误差；
[0046]
图6为仿真实验中高超声速飞行器的x轴的速度估计误差图；
[0047]
error in v
x
代表x轴的速度估计误差；
[0048]
图7为仿真实验中高超声速飞行器的y轴的速度估计误差图；
[0049]
error in vy代表y轴的速度估计误差；
[0050]
图8为仿真实验中高超声速飞行器的z轴的速度估计误差图；
[0051]
error in vz代表z轴的速度估计误差；
[0052]
图9为仿真实验中高超声速飞行器在观测坐标系的x轴向位置估计误差图；
[0053]
error in x代表x轴的位置估计误差；
[0054]
图10为仿真实验中高超声速飞行器在观测坐标系的y轴向位置估计误差图；
[0055]
error in y代表y轴的位置估计误差；
[0056]
图11为仿真实验中高超声速飞行器在观测坐标系的z轴向位置估计误差图；
[0057]
error in z代表z轴的位置估计误差；
[0058]
图12为仿真实验中高超声速飞行器在观测坐标系的x轴向加速度估计、预报与实际值对比图；
[0059]
图13为仿真实验中高超声速飞行器在观测坐标系的y轴向加速度估计、预报与实际值对比图；
[0060]
图14为仿真实验中高超声速飞行器在观测坐标系的z轴向加速度估计、预报与实际值对比图；
[0061]
图15为仿真实验中高超声速飞行器在观测坐标系的x轴向速度估计、预报与实际值对比图；
[0062]
图16为仿真实验中高超声速飞行器在观测坐标系的y轴向速度估计、预报与实际值对比图；
[0063]
图17为仿真实验中高超声速飞行器在观测坐标系的z轴向速度估计、预报与实际值对比图；
[0064]
图18为仿真实验中高超声速飞行器在观测坐标系的x轴向位置估计、预报与实际值对比图；
[0065]
图19为仿真实验中高超声速飞行器在观测坐标系的y轴向位置估计、预报与实际值对比图；
[0066]
图20为仿真实验中高超声速飞行器在观测坐标系的z轴向位置估计、预报与实际值对比图；
[0067]
图21为仿真实验中高超声速飞行器的位置估计、预报与实际值的误差曲线图；
[0068]
error in r代表估计、预报与实际值的位置误差；
[0069]
图22为仿真实验中高超声速飞行器终端位置预报误差统计结果图；
[0070]
p
rerror
代表终端位置预报误差；
[0071]
图23为仿真实验中高超声速飞行器状态量z
x
拟合曲线结果图；
[0072]
图24为仿真实验中高超声速飞行器状态量zy拟合曲线结果图；
[0073]
图25为仿真实验中高超声速飞行器状态量zz拟合曲线结果图；
[0074]
目标机动模式为常值加速度机动情况下的仿真：
[0075]
图26为仿真实验中高超声速飞行器在观测惯性坐标系中纵向平面飞行弹道图；
[0076]
图27为仿真实验中高超声速飞行器在观测惯性坐标系中侧向平面飞行弹道图；
[0077]
图28为仿真实验中高超声速飞行器的x轴的加速度估计误差图；
[0078]
图中：t代表时间，error in ay代表x轴的加速度估计误差；
[0079]
图29为仿真实验中高超声速飞行器的y轴的加速度估计误差图；
[0080]
error in ay代表y轴的加速度估计误差；
[0081]
图30为仿真实验中高超声速飞行器的z轴的加速度估计误差图；
[0082]
error in az代表z轴的加速度估计误差；
[0083]
图31为仿真实验中高超声速飞行器的x轴的速度估计误差图；
[0084]
error in v
x
代表x轴的速度估计误差；
[0085]
图32为仿真实验中高超声速飞行器的y轴的速度估计误差图；
[0086]
error in vy代表y轴的速度估计误差；
[0087]
图33为仿真实验中高超声速飞行器的z轴的速度估计误差图；
[0088]
error in vz代表z轴的速度估计误差；
[0089]
图34为仿真实验中高超声速飞行器在观测坐标系的x轴向位置估计误差图；
[0090]
error in x代表x轴的位置估计误差；
[0091]
图35为仿真实验中高超声速飞行器在观测坐标系的y轴向位置估计误差图；
[0092]
error in y代表y轴的位置估计误差；
[0093]
图36为仿真实验中高超声速飞行器在观测坐标系的z轴向位置估计误差图；
[0094]
error in z代表z轴的位置估计误差；
[0095]
图37为仿真实验中高超声速飞行器在观测坐标系的x轴向加速度估计、预报与实际值对比图；
[0096]
图38为仿真实验中高超声速飞行器在观测坐标系的y轴向加速度估计、预报与实际值对比图；
[0097]
图39为仿真实验中高超声速飞行器在观测坐标系的z轴向加速度估计、预报与实际值对比图；
[0098]
图40为仿真实验中高超声速飞行器在观测坐标系的x轴向速度估计、预报与实际值对比图；
[0099]
图41为仿真实验中高超声速飞行器在观测坐标系的y轴向速度估计、预报与实际值对比图；
[0100]
图42为仿真实验中高超声速飞行器在观测坐标系的z轴向速度估计、预报与实际值对比图；
[0101]
图43为仿真实验中高超声速飞行器在观测坐标系的x轴向位置估计、预报与实际值对比图；
[0102]
图44为仿真实验中高超声速飞行器在观测坐标系的y轴向位置估计、预报与实际值对比图；
[0103]
图45为仿真实验中高超声速飞行器在观测坐标系的z轴向位置估计、预报与实际
值对比图；
[0104]
图46为仿真实验中高超声速飞行器的位置估计、预报与实际值的误差曲线图；
[0105]
error in r代表估计、预报与实际值的位置误差；
[0106]
图47为仿真实验中高超声速飞行器终端位置预报误差统计结果图；
[0107]
p
rerror
代表终端位置预报误差；
[0108]
图48为仿真实验中高超声速飞行器状态量z
x
拟合曲线结果图；
[0109]
图49为仿真实验中高超声速飞行器状态量zy拟合曲线结果图；
[0110]
图50为仿真实验中高超声速飞行器状态量zz拟合曲线结果图；
[0111]
图51为飞行器可能落点范围及飞行走廊的示意图。
具体实施方式
[0112]
下面结合附图和现有技术对本发明方法进行进一步地详细说明：
[0113]
一种基于机动模式识别的临近空间高超声速飞行器弹道预测方法，预测方法的具体执行步骤如下：
[0114]
步骤一、建立临近空间高超声速飞行器的状态方程；所述步骤一的具体过程为：
[0115]
引入3个状态变量z
x
、zy和zz，且状态变量z
x
、zy和zz的表达式分别为：
[0116][0117]
其中：c
x0
、c
y0
和c
z0
分别代表高超声速飞行器在速度坐标系下的x0轴、y0轴和z0轴向的空气动力系数；c
p0
为推力系数；m
t
为高超声速飞行器的质量；s为高超声速飞行器的等效参考面积；
[0118]
选取状态变量x为：
[0119]
x＝[x y z v
x v
y v
z z
x z
y zz]
t
ꢀꢀꢀ
(2)
[0120]
其中：x、y和z分别为高超声速飞行器在观测惯性坐标系下x轴、y轴和z轴向的位置坐标；v
x
、vy和vz分别为高超声速飞行器在观测惯性坐标系下x轴、y轴和z轴向的速度；
[0121]
则建立临近空间高超声速飞行器的状态方程为：
[0122]
所述临近空间高超声速飞行器的状态方程包含加速度计算方程、速度计算方程、位置计算方程；
[0123][0124]
其中：为x的一阶导数，为y的一阶导数，为z的一阶导数，为v
x
的一阶导数，为vy的一阶导数，为vz的一阶导数，为z
x
的一阶导数，为zy的一阶导数，为zz的一阶导数；a
x
、ay和az分别为高超声速飞行器在观测惯性坐标系下x轴、y轴和z轴向的加速度；w1、w2和w3均为高斯白噪声，代表了慢时变的随机游走过程。
[0125]
所述中间变量a
tx
、a
ty
和a
tz
组成的中间变量矩阵以及中间变量g
x
、gy和gz组成的中间变量矩阵的表达式如下：
[0126][0127][0128]
其中：ρ为空气密度；代表高超声速飞行器的速度大小；利用弹道坐标系对速度坐标系进行替换，则为由弹道坐标系到观测惯性坐标系的变换矩阵(为由速度坐标系到观测惯性坐标系的变换矩阵，对于轴对称式飞行器，侧向机动一般由偏航
通道通过侧滑角进行控制，实际飞行过程中飞行器滚转角一般较小，因而速度坐标系可以用弹道坐标系进行替换，所以，代表由弹道坐标系到观测惯性坐标系的变换矩阵)；μ为引力系数；re为地球半径；
[0129]
弹道坐标系到观测惯性坐标系的变换矩阵的表达式为：
[0130][0131]
其中：
[0132][0133]
步骤二、设计用于临近空间高超声速飞行器状态跟踪的卡尔曼滤波器，利用设计的卡尔曼滤波器实时跟踪高超声速飞行器的状态信息；所述步骤二的具体过程为：
[0134]
所述用于临近空间高超声速飞行器状态跟踪的卡尔曼滤波器的预报方程为：
[0135][0136]
其中，为卡尔曼滤波器对状态变量x的估计值，表示使用状态变量的估计值下的高超声速飞行器系统状态方程；g是系统噪声矩阵；q为高超声速飞行器系统状态模型误差协方差矩阵；rn为测量误差协方差矩阵；上标-1
代表矩阵求逆运算；上标
t
代表矩阵转置运算；h为卡尔曼滤波器的测量矩阵；p为卡尔曼滤波器的协方差矩阵，为p对时间的一阶导数；f为雅克比矩阵；
[0137][0138][0139]
采用龙格库塔法求解公式(8)，得到中间变量和p(k+1/k)的值；
[0140]
所述用于临近空间高超声速飞行器状态跟踪的卡尔曼滤波器的测量修正方程为：
[0141][0142]
其中：i代表单位矩阵，k(k+1)代表卡尔曼滤波器第k+1步的滤波增益矩阵，η(k+1)
代表卡尔曼滤波器第k+1步的测量信息，为卡尔曼滤波器第k+1步的估计值，p(k+1)为卡尔曼滤波器第k+1步的协方差矩阵的值；
[0143]
η(k+1)＝hx+rnꢀꢀꢀ
(12)
[0144]
综合公式(8)和公式(11)，获得用于临近空间高超声速飞行器状态跟踪的卡尔曼滤波器，利用设计的卡尔曼滤波器实时跟踪高超声速飞行器的状态信息。
[0145]
步骤三、根据步骤二设计的卡尔曼滤波器获得跟踪开始时刻0至结束kt时刻的高超声速飞行器的状态信息，并根据跟踪得到的状态信息获得跟踪阶段的高超声速飞行器的状态估计值；
[0146]
其中，k为大于0的整数，t为计算周期，一般可取t＝0.01s；将跟踪结束kt时刻的高超声速飞行器的状态估计值作为弹道预测的初值；
[0147]
所述步骤三的具体过程为：
[0148]
根据步骤二设计的卡尔曼滤波器获得开始时刻0至跟踪结束kt时刻的高超声速飞行器的状态信息，并根据跟踪得到的的状态信息获得跟踪结束kt时刻的高超声速飞行器的状态估计值如下：
[0149][0150]
其中：x
p
(kt)、y
p
(kt)和z
p
(kt)为跟踪结束kt时刻的高超声速飞行器在观测惯性坐标系下x轴、y轴和z轴向的位置，v
xp
(kt)、v
yp
(kt)和v
zp
(kt)分别为跟踪结束kt时刻的高超声速飞行器在观测惯性坐标系下x轴、y轴和z轴向的速度；和分别为跟踪结束kt时刻对z
x
、zy和zz估计值；
[0151]
将x
p
(kt)、y
p
(kt)、z
p
(kt)、v
xp
(kt)、v
yp
(kt)和v
zp
(kt)作为弹道预测的初值。
[0152]
由式(3)至式(5)可知，对于临近空间高超声速飞行器，其在飞行末段的加速度由空气动力和地球引力产生，在高超声速飞行器进入最后下压之前，需要进行机动，由于飞行器的姿态角通常会发生变化；考虑到高超声速飞行器机动大的情况下，空气动力系数c
x
、cy、cz变化相对较大大。
[0153]
由式(1)可知在机动条件下，z
x
、zy、zz变化较大，且机动模式未知。因而根据跟踪阶段的估计值，利用聚类分析的方法辨识出高超声速飞行器的机动类型，进而确定预测方案，对状态量进行拟合，使用拟合的结果作为预报的根据，z
x
、zy、zz的值由拟合的结果确定。
[0154]
式(4)中ρ、v、和地球引力对加速度变化影响较大，空气密度ρ只与高超声速飞行器所在高度有关，而高超声速飞行器的高度与地球引力大小仅与状态量x、y、z有关，速度v及坐标变换矩阵只与状态v
x
、vy、vz有关。
[0155]
即由kt时刻观测坐标系下x轴、y轴和z轴向的位置和速度的预测值就可以获得(k+1)t时刻加速度的预测值，而由(k+1)t时刻加速度的预测值，根据式(3)又可以获得(k+1)t时刻的观测坐标系下x轴、y轴和z轴向的速度和位置的预测值，从而实现迭代计算，进行全程的预测。该过程是一个含有6个参数的微分方程，可使用龙格库塔法求解此微分方程。
[0156]
步骤四、根据跟踪阶段的估计值，利用聚类分析的方法辨识出的机动类型，进而确定预测方案，对状态量进行拟合，使用拟合的结果作为预报的根据，利用弹道预测的初值和对状态量的拟合的结果求解(k+1)t时刻高超声速飞行器的位置和速度的预测值，并继续利
用(k+1)t时刻的位置和速度的预测值以及对状态量的拟合的结果来求解(k+2)t时刻高超声速飞行器的位置和速度的预测值；所述步骤四的具体过程为：
[0157]
在对目标进行建模的过程中，状态变量z
x
、zy、zz是与目标加速度有直接的关系，使用估计过程中进入稳态的z
x
、zy、zz的估计值，根据z
x
、zy、zz的变化特点可将x-51的机动类型分为两种，即正弦加速度机动和常值加速度机动。首先根据聚类的方法判断x-51的机动类型，如果机动类型为正弦加速度机动，则使用模型f(t)＝a
1 sin(a2t+a3)对z
x
、zy、zz进行拟合得到函数z
x
(t)、zy(t)、zz(t)；如果机动类型为常值加速度机动，则使用模型f(t)＝a1t+a2对z
x
、zy、zz进行拟合得到函数z
x
(t)、zy(t)、zz(t)。将z
x
(t)、zy(t)、zz(t)代入公式(4)可以得到
[0158]
综合公式(3)至公式(7)获得如下的微分方程：
[0159][0160]
根据x
p
(kt)、y
p
(kt)、z
p
(kt)、v
xp
(kt)、v
yp
(kt)、v
zp
(kt)和采用龙格库塔法求解微分方程(14)，得到(k+1)t时刻高超声速飞行器的位置和速度在观测坐标系下的预测值；其中：x
p
((k+1)t)、y
p
((k+1)t)和z
p
((k+1)t)代表(k+1)t时刻高超声速飞行器的位置在观测坐标系下x轴、y轴和z轴向的预测值，v
xp
((k+1)t)、v
yp
((k+1)t)和v
zp
((k+1)t)代表(k+1)t时刻高超声速飞行器的速度在观测坐标系下x轴、y轴和z轴向的预测值；
[0161]
同理，根据x
p
((k+1)t)、y
p
((k+1)t)、z
p
((k+1)t)、v
xp
((k+1)t)、v
yp
((k+1)t)、v
zp
((k+1)t)和采用龙格库塔法求解微分方程(14)，得到(k+2)t时刻高超声速飞行器的位置和速度在观测坐标系下x轴、y轴和z轴向的预测值。
[0162]
聚类方法的具体过程为：
[0163]
步骤1、首先是要获取一定的样本
[0164]
将目标飞行器的机动模式设置为正弦形式，分别将加速度设置为不同频率的正弦
信号，只进行跟踪过程，将跟踪过程估计出来的数据存储起来记为x＝{x
1t
,x
2t
,......,x
nt
}
t
，其中n为跟踪过程重复的次数；将目标飞行器的机动模式设置为常值形式，分别将加速度设置为不同频率的常值信号，只进行跟踪过程，将跟踪过程估计出来的数据存储起来记为y＝{y
1t
,y
2t
,......,y
mt
}
t
，其中m为跟踪过程重复的次数。将已经分类的m+n个数据作为行向量，得到一个矩阵a＝{x
t
,y
t
}
t
，记a＝{a
1t
,a
2t
,......,a
m+nt
}
t
。
[0165]
步骤2、其次是将待测样本输入，与步骤1中获得的数据样本一起进行分类，分类采取线性判别分析(linear discriminant analysis,lda)的方法，可以看成一种线性分类器。lda希望投影后的数据，不同类别之间的距离尽可能大，同类别内部尽量紧凑。假设a∈a，通过线性投影变换表示为一维向量b，变换关系为：
[0166]
b＝wa
t
ꢀꢀꢀ
(15)
[0167]
然后进一步设定一个阈值b0，若b≥b0时判定属于第一类，否则属于第二类，于是lda的问题就可以理解为如何调整投影方向w，以保证在一维空间上b尽可能地区分两个类别。
[0168]
lda同时考虑类间和类内的关系。一方面，希望两个类别的数据投影后的均值之间的距离尽可能大。假设类别一的均值为类别二的均值为则投影后两类均值的距离为
[0169]
另一方面，希望同类别的数据尽可能集中。假设类别i中的数据表示为i＝1,2，则原始数据的类内协方差表示为
[0170][0171]
投影后数据的类内方差表示为：
[0172][0173]
结合类间距最大与类内方差最小，lda求解的优化问题为；
[0174][0175]
解得
[0176]
步骤五、不断重复步骤四的过程，直至获得每个时刻高超声速飞行器的位置和速度的预测值，完成临近空间高超声速飞行器的弹道预测。
[0177]
仿真部分
[0178]
本次仿真共进行400s，82～282s时刻对临近空间飞行器进行状态估计；
[0179]
如图1至图51所示，对于攻击方来讲，本发明所设计的临近空间飞行器，不管是巡航式高超声速飞行器还是滑翔式高超声速飞行器，在攻击目标点时会有末速度和攻击角度的要求，通常要求末速度和攻击角度达到某一最小值。仿真结果表明，若要求末速度大于2mach，攻击角度大于70
°
的情况下，在进入最终下压点时速度需要在2800m/s左右。在这种实际情况的制约下，临近空间高超声速飞行器在末段的机动能力已经受到了约束限制。若在距目标点1000km处飞行器速度大小及方向已知的情况下，飞行器可能的落点范围约为以
此目标点为圆心半径50km的圆形区域。在此情况下，若飞行器最终落点确定，则末段飞行走廊宽度为仅为20km左右(z轴正负方向各10km左右)。飞行器可能落点范围、飞行走廊宽度等关系如图51所示。在目标机动模式未知的条件下，根据滤波器跟踪得到的状态量可以识别目标的机动模式，并以此作为状态拟合模型的选取依据。选用适当的拟合模型对状态量进行拟合，所得到的结果如图23、图24、图25和图48、图49、图50所示。
[0180]
在本发明的仿真中，直接将飞行器最终攻击点定为了观测坐标系的原点，因此飞行器最大横向机动范围为20km左右。若飞行器攻击点在可能落点范围内任意取值而与观测坐标系原点不重合，并不会对状态估计与预测的结果产生影响。此外，对于临近空间飞行器的弹道预测，知道该飞行器的大致目标点是必要的，雷达等探测装置通常部署在防御目标附近。同时，进行弹道预测的目的是对临近空间飞行器进行后续的拦截做准备，而只有预测出下压点之前的弹道才具有拦截的意义。如果敌方飞行器的大致攻击范围不确定，飞行器随时可能进入下压段进行攻击，这种情况是难以预测拦截的。由图51可知，在飞行器可能攻击范围限定为50km半径的圆后，最终下压点横向取值范围约为120km，本发明研究的意义在于经过200s时间的状态估计，可通过弹道预测将400s时刻下压点之前120km宽度的走廊范围缩小到小于10km的范围。
[0181]
在上述条件下，在目标机动模式为正弦加速度机动和常值加速度机动时分别做100次monte-carlo仿真，经过100次monte-carlo仿真得到的临近空间高超声速飞行器最终位置预报误差的统计结果如图22和图47所示，目标机动模式为正弦加速度机动时统计计算得出终端位置预报误差均小于10km，目标机动模式为常值加速度机动时统计计算得出终端位置预报误差均小于3km。
[0182]
本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。