本发明涉及残骸分离后的落地预测技术，具体是一种基于统计原理的实时积分外推落点预报方法。
背景技术：
一般卫星发射不同于载人航天回收任务，各子级残骸分离后是没有外测测量设备继续跟踪的，因此，在进行落点计算时，残骸初始运动信息是无法通过直接测量得出的，而是通过分离时刻运载器弹道信息估算得出的。其原因是运载器分离时刻箭体运动矢量过于复杂(除了在射向方向加速运动外，箭体自身相对运动方向还有方位、俯仰和滚动等矢量)，同时分离时的相对作用力大小、方向和作用时间都是一个模糊的控制量。既然精确计算比较困难，任务实施中就只能利用经验值的模式进行修正。
在落点计算中，航天器被动段的运动可视为质点运动，我们仅从物体运动学原理方面理解，航天器被动段飞行主要受地球引力、哥氏力和空气阻力的影响。被动段飞行从受力不同上可分为自由飞行段和再入段，自由飞行段航天器只在地球引力和自转引起的哥氏力的作用下运动，再入段与自由飞行段的差别仅仅增加了空气阻力的影响。运动方程如下所示：
上式运动方程中，空气阻力项：xd
其中sm为残骸迎风面的面积，v为残骸综合速度，g0为地面重力加速度，ρ/ρ0为跟高度有关的空气密度，可通过查落区大气密度表得到，cd为空气阻力系数，其值与残骸速度有关，可通过查马赫数对应的气阻系数表得到。从公式可以看出，在任务实际执行中，残骸迎风面的面积的精确值不容易获取，随高度变化的空气密度获取也不够精准，这些都是导致计算大气阻力不精确的因素。
落点计算是由运动方程构成的一个微分方程，在进行数值计算中需设计一定的步长进行迭代计算，选择适当的步长是非常重要的，特别是实时任务中，步长过小会导致计算量过大而影响实时性需求，过大则造成计算误差过大的情况发生。
技术实现要素：
为解决上述问题，本发明的目的在于提供一种基于历史数据统计反演的实时落点预报方法。具体涉及残骸分离信息计算、大气阻滞补偿计算、无动力陨落加速度计算、数值迭代计算，以及基于历史真实数据的大气阻滞综合系数、修正参数δv和δk的数值反演。
本发明具体包括以下，(1)在分离点弹道计算时采用单向高阶最小二乘非线性拟合来消除分离时刻弹道为拐点的影响，同时为进一步提高分离时刻弹道的精度，对残骸初始运动速度信息进行了修正；(2)在常规大气阻力模型的基础之上，通过分析大气密度随高度h的分布规律，用最小二乘拟合了大气密度随高度变化的指数，结合火箭残骸的截面、重量、马赫数的影响，建立了含大气阻滞综合系数、马赫数影响系数以及大气密度指数的大气阻滞模型，同时为进一步提高大气阻滞综合系数精度，对大气阻滞综合系数进行补偿修正；(3)综合考虑地球引力、哥氏力和空气阻力等因素的影响，计算无动力陨落加速度；(4)进行数值迭代计算落点，实现落点的预报，采用5ms步长的“龙格-库塔”模型进行迭代计算，迭代计算终止标识一般设置为落点相对地球椭球体表面的高程h小于500米；(5)结合分离作用力影响和历史经验，设计了修正参数的取值范围，δv取值范围为[2，15]和δk取值范围为[0.5,4.5]；(6)基于历史数据统计反演大气阻滞系数、大气阻滞修正系数δk、修正速度δv。
鉴于此，本发明采用的技术方案是，一种基于统计的重残骸实时落点预报方法，包括以下步骤：
基于历史数据进行反演，获取最佳修正速度δv和最佳修正系数δk。
从遥测时间中获取最佳残骸分离时间信息。
从运载器众多飞行弹道中挑选最优飞行弹道。
计算残骸分离时刻弹道信息，包括弹道速度信息和弹道位置信息。
根据最佳修正速度δv修正残骸分离模型中残骸分离时刻弹道速度信息。
根据最佳修正系数δk修正大气阻滞补偿模型。
计算大气阻力。
根据无动力陨落加速度模型计算无动力加速度。
按迭代算法和时间迭代步长进行数值迭代。
判断迭代是否满足终止条件，若满足则输出残骸预示落点，若不满足则执行计算大气阻力步骤。
采用上述技术的原因是落点计算时无法准确获取残骸初始运动信息和大气阻力。
本发明的有益技术效果是：
1、充分考虑了运载器分离时刻运动过程复杂性和弹道特性，采用单向高阶最小二乘非线性拟合来消除分离时刻弹道为拐点的影响。
2、对残骸初始运动速度信息进行修正，进一步提高残骸初始运动速度信息精度。
3、综合考虑多方面因素，建立大气阻滞模型，基于历史数据统计反演大气阻滞系数。
4、对大气阻滞系数进行修正，进一步提高大气阻滞系数精度。
5、结合分离作用力影响和历史经验，设计了修正参数(修正系数和修正速度)的取值范围。
6、基于历史数据统计反演大气阻滞修正系数和修正速度。
7、结合计算精度和实时计算需求，设计了合适的迭代方法。
8、本发明编程易实现，提高了实时计算落点的效率和精度。
附图说明
图1为本发明的处理流程图。
具体实施方式
为使本发明的目的、技术方案及优点得到更好的理解，现以具体实施例并结合附图，具体说明如下：
基于统计原理的实时落点预报方法是以常规运动模型为基础建立，根据综合误差模糊计算理论，在基础计算模型中增加两个模糊修正量，一个用来补偿计算与实际运动的速度误差；一个用来调节模型计算与实际风阻的关系，同时根据重残骸运动特性，设定值域减少迭代计算次数，包括：
(1)建立残骸分离信息计算模型
进行残骸落点计算第一步是需要知道残骸分离时刻的弹道信息，然后才能根据加速度情况进行数值计算。这里之所以采用高阶拟合计算，主要原因是因为在残骸分离前后，其主体火箭加速度变化较大，如果利用线性拟合计算将导致计算结果误差较大。模型利用分离前、后弹道信息分别拟合，最后利用两组结果求均值得出。
设yj,j＝1,2,……n为(n为多项式回归拟合段落点数)弹道参数(包括位置信息和速度信息)时间序列t的第j个采样点测量值；bj为yj序列的多项式回归系数；p为拟合多项式的阶数(取为2,3,4……；)，解算模型为：
b＝(dtd)-1dty(1)
其中：
bj,j＝1,2,……p为待求多项式回归系数，p为待求多项式阶数；ti,i＝1,2,……n为弹道序列对应的时间，yi,i＝1,2,……n为为弹道信息中时间序列t的第j个采样点测量值。
y＝b0 b1x … bpxp,p＝2,3…(3)
公式(3)即为最小二乘非线性拟合方程，将分离点时刻带入上式即可求出拟合弹道信息，本方法仿真计算中拟合多项式阶数取4。x表示数据拟合时间，y表示x时刻的拟合值。
假设经模型公式(3)解算的残骸分离的速度、位置信息分别为：v0＝(vx0vy0vz0)，p0＝(x0y0z0)。vx0vy0vz0分别表示速度在x方向的分量、速度在y方向的分量、速度在z方向的分量，x0y0z0分别表示位置在x方向的分量、位置在y方向的分量、位置在z方向的分量。加入补偿修正量(δv)后，其分离信息(根据分离特性，只做速度补偿，位置补偿不做修正)为：
式中
xd＝x0；yd＝y0；zd＝z0。
(2)建立大气阻滞补偿计算模型
大气阻滞加速度矢量模型推导过程较为复杂，严格意义上需要考虑残骸陨落过程空域的大气密度、风向等指标，实际任务中，往往采用地区通用大气密度数据进行计算，本发明中模型是在常规大气阻力模型的基础之上，通过分析大气密度随高度h的分布规律，用最小二乘拟合了大气密度随高度变化的指数，结合火箭残骸的截面、重量、马赫数的影响，建立了含大气阻滞综合系数、马赫数影响系数以及大气密度指数的大气阻滞模型：
gbx表示阻力加速度在x方向的分量，gby表示阻力加速度在y方向的分量，gbz表示阻力加速度在z方向的分量，vx表示速度在x方向的分量，vy表示速度在y方向的分量，vz表示速度在z方向的分量，v表示综合速度。
指数id是一个跟高度(公里)有关的一个高阶拟合值，拟合方程如下：
a1……a8，b1……b6表示指数系数，如表1所示。hk表示高度。
参数kv是与残骸飞行的马赫数有关的一个拟合值。
cv表示马赫数；k为大气阻滞综合系数。
表1指数系数表
补偿系数δk带入公式(5)，从而形成带补偿系数的大气阻滞计算模型。
(3)建立无动力陨落加速度模型
残骸陨落过程中，自身已无动力系统工作，下坠过程中除受大气阻力影响外，就是受引力影响，这里的引力一般指地球引力与因地球自转而产生的哥氏力，此力是残骸下坠过程中影响最大，计算模型如下：
其中gxd表示重力加速度在x方向的分量，gyd表示重力加速度在y方向的分量，gzd表示重力加速度在z方向的分量。
上式中里斯常数j＝0.00108263，u＝3.986004415e14为地球引力常数，we＝0.000072921151467为地球自转角速率，地球半径re＝6378245.0。
残骸无动力陨落过程中，其运动加速度由两部分构成，一个是引力加速度，一个是大气阻滞加速度，因此残骸无动力陨落加速度模型可由公式(6)和公式(7)得出：
(4)进行落点数值迭代计算方法
残骸陨落轨迹和落点的计算，根据物体运动学原理，已知物体运动初始速度和位置，在加速度可计算的条件下，即可计算出下一时刻的速度和位置，在进行积分迭代计算中，时间步长不宜取的过大，一般不超过10毫秒，实际计算中所取步长越小，其迭代计算精度越高，但同时计算所用资源也随着步长的减少而成倍增加，实际应用中需结合计算需要和计算机运算速度来综合确定(仿真计算t取5ms)。迭代计算公式见公式(9)。在结合分离作用力影响和历史经验，在进行计算中δv取值范围为[2，15]，步长取0.24，δk取值范围为[0.5,4.5],步长取1/50.0。
参见图1，基于历史数据统计反演修正系数流程：
1、历史数据收集整理，包括残骸实际落点、残骸分离时间、运载器飞行弹道等。
2、结合分离作用力影响和历史经验，设计δv取值范围[2，15]，迭代步长0.24，δk取值范围[0.5,4.5],迭代步长1/50.0。
3、设计迭代算法、时间迭代步长、迭代终止条件等。
4、计算残骸分离时刻弹道信息。
5、对残骸分离时刻弹道速度信息按修正步长进行修正，并判断修正量是否超出取值范围，若超出取值范围执行步骤11，若未超出取值范围执行步骤6。
6、对大气阻滞系数按修正步长进行修正，并判断修正量是否超出取值范围，若超出取值范围执行步骤5，若未超出取值范围执行步骤7。
7、计算大气阻力。
8、计算无动力加速度。
9、按迭代算法和时间迭代步长进行数值迭代。
10、判断迭代是否满足终止条件，若满足则执行步骤6，若不满足执行步骤7。
11、将计算结果与历史实际落点进行对比分析。
12、选取最优修正系数。
基于统计原理的实时落点预报：
1、设计迭代算法、时间迭代步长、迭代终止条件等。
2、装填最优修正系数。
3、从遥测时间中获取最佳残骸分离时间信息。
4、从运载器众多飞行弹道中挑选最优飞行弹道。
5、计算残骸分离时刻弹道信息。
6、修正残骸分离时刻弹道速度信息。
7、修正大气阻滞系数。
8、计算大气阻力。
9、计算无动力加速度。
10、按迭代算法和时间迭代步长进行数值迭代。
11、判断迭代是否满足终止条件，若满足则执行步骤12，若不满足执行步骤7。
12、输出残骸预示落点。
1.一种基于统计的重残骸实时落点预报方法，其特征在于，包括以下步骤：
基于历史数据进行反演，获取最佳修正速度δv和最佳修正系数δk；
从遥测时间中获取最佳残骸分离时间信息；
从运载器众多飞行弹道中挑选最优飞行弹道；
计算残骸分离时刻弹道信息，包括弹道速度信息和弹道位置信息；
根据最佳修正速度δv修正残骸分离模型中残骸分离时刻弹道速度信息；
根据最佳修正系数δk修正大气阻滞补偿模型；
计算大气阻力；
根据无动力陨落加速度模型计算无动力加速度；
按迭代算法和时间迭代步长进行数值迭代；
判断迭代是否满足终止条件，若满足则输出残骸预示落点，若不满足则执行计算大气阻力步骤。
2.根据权利要求1所述一种基于统计的重残骸实时落点预报方法，其特征在于：所述残骸分离时刻弹道信息利用分离前、后弹道信息分别拟合，最后利用两组结果求均值得出。
3.根据权利要求2所述一种基于统计的重残骸实时落点预报方法，其特征在于：所述弹道信息拟合采用单向高阶最小二乘非线性拟合。
4.根据权利要求1所述一种基于统计的重残骸实时落点预报方法，其特征在于：所述修正后残骸分离模型为
式中：残骸分离时刻的速度信息v0＝(vx0vy0vz0)，vx0vy0vz0分别表示速度在x方向的分量、速度在y方向的分量、速度在z方向的分量，残骸分离时刻的位置信息p0＝(x0y0z0)，x0y0z0分别表示位置在x方向的分量、位置在y方向的分量、位置在z方向的分量，修正后的速度信息和位置信息分别为vd和pd；
xd＝x0；yd＝y0；zd＝z0。
5.根据权利要求1所述一种基于统计的重残骸实时落点预报方法，其特征在于：所述修正后的大气阻滞补偿模型为
式中，gbx表示阻力加速度在x方向的分量，gby表示阻力加速度在y方向的分量，gbz表示阻力加速度在z方向的分量，vx表示速度在x方向的分量，vy表示速度在y方向的分量，vz表示速度在z方向的分量，v表示综合速度；
指数id的拟合方程如下：
a1……a8，b1……b6表示指数系数；
参数kv是与残骸飞行的马赫数有关的一个拟合值；
cv为马赫数；k为大气阻滞综合系数。
6.根据权利要求1所述一种基于统计的重残骸实时落点预报方法，其特征在于：残骸无动力陨落过程中，其运动加速度由两部分构成，一个是引力加速度，一个是大气阻滞加速度，因此所述无动力陨落加速度模型为
gz＝(gzx,gzy,gzz)
gxd表示重力加速度在x方向的分量，gyd表示重力加速度在y方向的分量，gzd表示重力加速度在z方向的分量，gbx表示阻力加速度在x方向的分量，gby表示阻力加速度在y方向的分量，gbz表示阻力加速度在z方向的分量。
7.根据权利要求6所述一种基于统计的重残骸实时落点预报方法，其特征在于：所述引力加速度的计算模型为
其中
上式中j为里斯常数，u为地球引力常数，we为地球自转角速率，re为地球半径。
8.根据权利要求1-7任一项所述一种基于统计的重残骸实时落点预报方法，其特征在于：所述基于历史数据进行反演，包括以下步骤：
1)历史数据收集整理，包括残骸实际落点、残骸分离时间、运载器飞行弹道；
2)结合分离作用力影响和历史经验，设计修正速度和修正系数；
3)设计迭代算法、时间迭代步长、迭代终止条件；
4)计算残骸分离时刻弹道信息；
5)对残骸分离时刻弹道速度信息按修正速度进行修正，并判断修正量是否超出取值范围，若超出取值范围执行步骤11)，若未超出取值范围执行步骤6)；
6)对大气阻滞系数按修正系数进行修正，并判断修正量是否超出取值范围，若超出取值范围执行步骤5)，若未超出取值范围执行步骤7)；
7)计算大气阻力；
8)计算无动力加速度；
9)按迭代算法和时间迭代步长进行数值迭代；
10)判断迭代是否满足终止条件，若满足则执行步骤6)，若不满足执行步骤7)；
11)将计算结果与历史实际落点进行对比分析；
12)选取最佳修正速度δv和最佳修正系数δk。
9.根据权利要求8所述一种基于统计的重残骸实时落点预报方法，其特征在于：结合分离作用力影响和历史经验，设计的修正速度δv和修正系数δk，δv取值范围为[2，15]和δk取值范围为[0.5,4.5]。
技术总结