基于能量分析的航天器LEO圆轨道衰降过程参数预报方法

摘要

本发明公开了基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，涉及航天工程及轨道动力学、航天器速度测定领域，包括：由LEO圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量；根据时间步长推进计算飞行过程中大气阻尼对动能的耗散，从而得到对总能量的耗散；根据耗散后的总能量计算新的圆轨道运行参数；根据轨道高度判断是否满足再入条件，不满足返回上一步骤，满足则继续下一步骤；得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论，包括自初始时刻至再入时刻的时间预报。适用于该类运行模式的航天飞行器的轨道衰降过程参数的中长期分析预报，为非回收航天器的再入预警和地面风险评估提供重要技术支撑。

权利要求

1.基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，包括以下步骤：步骤一、由LEO圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量，其中总能量包括动能和势能；

步骤二、根据时间步长推进计算飞行过程中大气阻尼对动能的耗散，从而得到对总能量的耗散；

步骤三、根据耗散后的总能量计算新的圆轨道运行参数，包括圆轨道半径或高度；

步骤四、根据轨道高度判断是否满足再入条件，不满足返回步骤二，满足则继续下一步骤；

步骤五、得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论，包括自初始时刻至再入时刻的时间预报；

所述步骤二中，大气阻尼对总能量的耗散的计算方法包括：S20、设时间间隔dt内，航天器运动轨迹近似为ds＝Vdt，其中V为航天器运行速度；WA为气动阻力对航天器做功，即大气阻尼对总能量的耗散，其计算公式如下：其中，气动阻力做的功取决于做功距离ds；实际处理时，时间步长内不考虑气动阻力变化，即时间步长内不考虑相关的阻力系数CD、密度ρ、航天器运行速度V及迎风截面积A的变化；

S21、所述S20的公式中，气动阻力做功的微分形式或表示如下：其中，通过换算关系把气动阻力做功的表达式表示为航天器距地心距离H的函数，G为万有引力常数，mE为地球质量；

S22、所述S20的公式中，FA为气动阻力，表达式为：其中，阻力系数CD的参考面积取为单位面积；

所述步骤三中，根据耗散后的总能量计算新的圆轨道运行参数的计算方法包括：S30、设航天器距地心距离H在时间步长dt前后分别是H1和H2，则根据耗散后总能量计算新的轨道运行参数，其计算公式如下：其中，ρ为密度，CD为航天器飞行过程中的阻力系数，A为迎风截面积，m为航天器质量，mE为地球质量，G为万有引力常数；

S31、所述S30的公式中，dH的计算依据即为航天器总能变化dE等于大气阻尼对其耗散，即等于大气阻力对航天器做功WA的负值，具体关系如下：dE＝‑WA

或者

由此即可得出dH的表达式，CD为阻力系数。

2.如权利要求1所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述步骤一中，航天器相对地心准惯性坐标系的总能量计算方法包括：S10、设任意时刻轨道高度为h，则LEO圆轨道至地心距离为H＝h+rE，其中rE为地球平均半径；于是航天器在轨运行时的总能量为E，其计算公式如下：‑11 2 2

其中，mE为地球质量，m为航天器质量，G为万有引力常数，其值为6.67×10 N.m/kg；

S11、所述S10的公式中，航天器动能E动的计算公式如下：其中，V为航天器运行速度，其和航天器运行高度的关系由引力和离心力的平衡关系进行确定；

S12、所述S10的公式中，航天器势能E势的计算公式如下：其中，航天器的势能零点取在无穷远处；

S13、把航天器初始参数即初始时刻的飞行轨道高度代入S10的公式中即可获得初始时刻航天器的总能量。

3.如权利要求1所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述步骤四中，再入条件设为某个临界高度he，判断是否满足再入条件的方法包括：S40、根据S30公式得到的航天器距地心距离得到航天器运行高度h2＝H2‑rE，当满足下列条件则说明满足再入条件：h2≤he

于是进入步骤五；

S41、根据S30公式得到的航天器距地心距离得到航天器运行高度h2＝H2‑rE，当满足下列条件则说明尚未满足再入条件，需要再返回步骤二推进计算：h2＞he

于是返回步骤二，在给定推进时间步长后继续计算大气阻尼对总能量的耗散，并获得新的轨道参数。

4.如权利要求3所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述步骤五中，在步骤五中，得到再入点参数及统计计算轨道衰降过程的重要特征参数，它们的计算方法包括：S50、根据步骤四设定的再入高度he，即某个临界高度，再入速度Ve的计算公式如下：该公式同时适用于任意高度时的飞行速度计算，即S51、根据步骤四设定的再入高度he，再入时刻当地弹道倾角de的计算公式如下：de＝arctan(dhe/Vedte)

其中，dte、dhe、Ve分别是再入时刻时的时间步长、高度增量及再入速度；该公式同时适用于任意高度时的当地弹道倾角的计算，即d＝arctan(dh/Vdt)；

S52、根据步骤二至步骤四推进计算过程对时间步长的累加，即可得到从初始状态到再入时的时间预报，时间预报ttotal的计算公式如下：ttotal＝∑dt

根据需要，采用不同时间单位，包括“天”、“周”或“年”对时间预报值进行换算。

5.如权利要求4所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述临界高度he的高度值在90km至120km之间。

说明书

基于能量分析的航天器LEO圆轨道衰降过程参数预报方法

技术领域

[0001] 本发明属于航天工程及轨道动力学、航天器速度测定技术领域，更具体地说，本发明涉及基于能量分析的航天器LEO圆轨道衰降过程参数预报方法。

背景技术

[0002] 服役期满非回收航天器，包括各种人造卫星、空间实验室、空间站及运载火箭末级等，一般都要经历轨道衰降过程后再入稠密大气层陨落解体。近年来，空间碎片问题引起普遍关注，各国科研人员一直在对航天器再入的中长期预报进行研究，准确可靠的再入预报能够提醒对地面人员设施的安全预警，及时回应国际关切和舆论。对轨道衰降阶段的预报是落区预报的前提，准确可靠的再入点参数对后续陨落解体及地面风险评估具有重要影响。但是，航天器在轨飞行的空间环境非常复杂且多变，其运动轨迹受到的影响因素包括地球/月球/太阳引力、大气密度模式、重力梯度矩、太阳辐射压力、地磁作用、地球非质点非对称质量分布等等；因此，要精确求解航天器在轨运行轨迹的相关参数需要复杂的轨道动力学理论及求解复杂的运动学和动力学方程组，同时需要对预期的一些影响因素进行研究和尽可能可靠的设定。

[0003] 目前，轨道参数的定解计算常用摄动方法，包括各种坐标摄动法、瞬时椭圆法和正则变换；摄动方程求解方法包括积分法、解析法和半解析法。无论那种摄动方法，都涉及摄动方程建模和求解，模型的精细程度取决于摄动小量的展开级数表达式及其展开阶数；数学上模型的精细化和复杂性将带来极大计算代价。用数值法求解轨道寿命，传统方法诸如考虑大气摄动加速度的科威尔（Cowell）法、恩克（Encke）法或轨道要素变动法，这些方法由于时间步长不能取得太大，当轨道高度较高因而衰减非常缓慢时，计算需要非常长的时间，甚至难以实现。现有的解析法虽然计算简便快捷，但仅适用于人为假设的指数型大气模型及一些特定简化模式，不能普遍使用。

[0004] 航天器轨道衰降问题即空间碎片减缓的主要基本问题。空间碎片减缓的主动措施之一是主动变轨快速再入，其次就是针对一些卫星要求在使命完成后进行离轨机动，进入一条短寿命（通常要求在25年内）的轨道，以免其长期滞留太空，污染太空环境。另一方面也有一些航天器丧失了变轨机动能力。上述情况都需要对其自主衰降过程进行预报。

[0005] 针对LEO圆轨道衰降参数的计算分析预报，关键的是一个目标参数即飞行高度的变化及衰降计算，并由此获得对再入时间和部分重要参数进行中长期预报的结论。LEO指的是低地球轨道，一般是高度低于2000km的轨道。LEO圆轨道，是若干人造卫星或空间实验室、空间站等航天器绕地球运行飞行时常见的一种轨道模式，其高度在一个运行周期内变化极小。当实际轨道高度在一个运行周期内变化范围不超过50km时，可以近似认为是圆轨道。经验表明，轨道自然衰降再入中长期预报时，即使仅提前几天，要对再入时刻位置（空间位置或星下点）进行有实用意义的确定也几乎是不可能的。由于再入位置对各种影响因素极为敏感，因而可以把分析重点放在对再入日期的高效计算预报方面。因此，专注于对轨道高度直接相关参数的精细分析，简化或忽略轨道倾角、升交点赤经等参数，从能量分析角度构建一种LEO圆轨道衰降参数的中长期预报快速高效方法即成为本发明的主要内容。

发明内容

[0006] 本发明的一个目的是解决上述问题和/或缺陷，并提供后面将说明的优点。

[0007] 为了实现根据本发明的这些目的和优点，提供了基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，包括以下步骤：

[0008] 步骤一、由LEO圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量，其中总能量包括动能和势能；

[0009] 步骤二、根据时间步长推进计算飞行过程中大气阻尼对动能的耗散，从而得到对总能量的耗散；

[0010] 步骤三、根据耗散后的总能量计算新的圆轨道运行参数，包括圆轨道半径或高度；

[0011] 步骤四、根据轨道高度判断是否满足再入条件，不满足返回步骤二，满足则继续下一步骤；

[0012] 步骤五、得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论，包括自初始时刻至再入时刻的时间预报。

[0013] 优选的是，其中，所述步骤一中，航天器相对地心准惯性坐标系的总能量计算方法包括：

[0014] S10、设任意时刻轨道高度为h，则LEO圆轨道至地心距离为H=h+rE，其中rE为地球平均半径；于是航天器在轨运行时的总能量为E，其计算公式如下：

[0015]

[0016] 其中，mE为地球质量，m为航天器质量，G为万有引力常数，其值为6.67×10‑11（N.m2/2
kg）；

[0017] S11、所述S10的公式中，航天器动能E动的计算公式如下：

[0018]

[0019] 其中，V为航天器飞行速度，其和航天器运行高度的关系由引力和离心力的平衡关系进行确定；

[0020] S12、所述S10的公式中，航天器势能E势的计算公式如下：

[0021]

[0022] 其中，航天器的势能零点取在无穷远处；

[0023] S13、把航天器初始参数即初始时刻的飞行轨道高度代入S10的公式中即可获得初始时刻航天器的总能量。

[0024] 优选的是，其中，所述步骤二中，大气阻尼对总能量的耗散的计算方法包括：

[0025] S20、S20、设时间间隔dt内，大气阻尼对总能量的耗散，即气动阻力对航天器做的功WA，其计算公式如下：

[0026]

[0027] 其中，航天器运动轨迹近似为ds=Vdt，其中V为航天器飞行速度；气动阻力做的功主要取决于做功距离ds；实际处理时，时间步长内不考虑气动阻力变化，即时间步长内不考虑相关的阻力系数CD、大气密度ρ、航天器飞行速度V及迎风截面积A的变化；同时，在不同时间过程中大气密度ρ随高度变化采用相关大气模型，航天器飞行速度V随航天器的高度变化根据航天器动能计算，航天器迎风截面积A根据实际情况取合适的量值，即根据航天器实际垂直于风向的净面积而定，必要时可变化，阻力系数CD取定值或根据航天器高度和航天器几何构型和姿态进行精细化标定；

[0028] S21、所述S20的公式中，气动阻力做功的微分形式还表示如下：

[0029]

[0030] 其中，通过换算关系把气动阻力做功的表达式表示为航天器距地心距离H的函数，G为万有引力常数，mE为地球质量；

[0031] S22、所述S20的公式中，FA为气动阻力，表达式为：

[0032]

[0033] 其中，阻力系数CD的参考面积取为单位面积。

[0034] 优选的是，其中，所述步骤三中，根据耗散后的总能量计算新的圆轨道运行参数的计算方法包括：

[0035] S30、设航天器距地心距离H在时间步长dt前后分别是H1和H2，则根据耗散后总能量计算新的轨道运行参数，其计算公式如下：

[0036]

[0037] 其中，ρ为大气密度，CD为航天器飞行过程中的阻力系数，A为迎风截面积，m为航天器质量，mE为地球质量，G为万有引力常数；

[0038] S31、所述S30的公式中，dH的计算依据即为航天器总能变化dE等于大气阻尼对其耗散，即等于大气阻力对航天器做功WA的负值，具体关系如下：

[0039] dE=‑WA

[0040] 或者

[0041]

[0042] 由此即可得出dH的表达式，CD为阻力系数。

[0043] 优选的是，其中，所述步骤四中，再入条件设为某个临界高度he，判断是否满足再入条件的方法包括：

[0044] S40、根据S30公式得到的航天器距地心距离得到航天器运行高度h2=H2‑rE，当满足下列条件则说明满足再入条件：

[0045] h2≤he

[0046] 于是进入步骤五；

[0047] S41、根据S30公式得到的航天器距地心距离得到航天器运行高度h2=H2‑rE，当满足下列条件则说明尚未满足再入条件，需要再返回步骤二推进计算：

[0048] h2＞he

[0049] 于是返回步骤二，在给定推进时间步长后继续计算大气阻尼对总能量的耗散，并获得新的轨道参数。

[0050] 优选的是，其中，所述步骤五中，得到再入点参数及统计计算轨道衰降过程的重要特征参数，它们的计算方法包括：

[0051] S50、根据步骤四设定的再入高度he，即某个临界高度，再入速度Ve的计算公式如下：

[0052]

[0053] 该公式同时适用于任意高度时的飞行速度计算，即 ；

[0054] S51、根据步骤四设定的再入高度he，再入时刻当地弹道倾角de的计算公式如下：

[0055] de=arctan（dhe/Vedte）

[0056] 其中，dte、dhe、Ve分别是再入时刻时的时间步长、高度增量及再入速度；该公式同时适用于任意高度时的当地弹道倾角的计算，即d=arctan（dh/Vdt）；

[0057] S52、根据步骤二至步骤四推进计算过程对时间步长的累加，即可得到从初始状态到再入时的时间预报，时间预报ttotal的计算公式如下：

[0058]

[0059] 根据需要，采用不同时间单位，如“天”、“周”或“年”对时间预报值进行换算。

[0060] 优选的是，其中，所述临界高度he的高度值在90km至120km之间。

[0061] 本发明至少包括以下有益效果：

[0062] 其一，本发明通过简单的积分方法，即可快速高效进行航天器LEO圆轨道衰降过程的主要参数进行计算，从而提供了一种简洁的航天器再入陨落时间的中长期预报方法。

[0063] 其二，本发明构造了LEO圆轨道衰降过程高度变化的算法流程，从而在此基础上推算速度、倾角等其它参数，也可以推算出更多的导出参数，包括飞行角速度、每天衰降高度、累计飞行圈数、飞行航迹航程、累计飞行时间等等。

[0064] 其三，本发明的LEO圆轨道衰降过程参数预报方法，仅存在一个难以确定的参数，即航天器飞行时的迎风截面积。对于复杂外形及姿态变化情况下的航天器，该面积参数可能存在较大不确定性，但仍然可以通过一定的观测数据对该迎风面积进行动态标定。

[0065] 其四，本发明基于能量分析的航天器LEO圆轨道衰降过程参数预报方法是基于物理定律得到的应用手段，其预报结果的准确性还取决于大气密度模式，为此建议采用较为精细的大气模式算法。在这些方面有若干成熟的研究成果能够直接借鉴

[0066] 本发明的优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

[0067] 图1为航天器LEO圆轨道衰降过程参数分析示意图；

[0068] 图2为TG‑1轨道衰降过程参数预报结果和官方发布参数；

[0069] 图3为TG‑1轨道衰降过程参数再入前最后几天高度结果；

[0070] 图4为TG‑1轨道衰降过程飞行速度及高度变化历程；

[0071] 图5为TG‑1轨道衰降过程飞行当地倾角变化历程；

[0072] 图6为TG‑1轨道衰降过程飞行每日高度衰降变化历程。

具体实施方式

[0073] 下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

[0074] 应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不排除一个或多个元件或其组合的存在或添加。

[0075] 如图1所示：本发明的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，包括以下步骤：

[0076] 步骤一、由LEO圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量，其中总能量包括动能和势能；

[0077] 步骤二、根据时间步长推进计算飞行过程中大气阻尼对动能的耗散，从而得到对总能量的耗散；

[0078] 步骤三、根据耗散后的总能量计算新的圆轨道运行参数，包括圆轨道半径或高度；

[0079] 步骤四、根据轨道高度判断是否满足再入条件，不满足返回步骤二，满足则继续下一步骤；

[0080] 步骤五、得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论，包括自初始时刻至再入时刻的时间预报。

[0081] 在上述技术方案中，所述步骤一中，绕地运行航天器相对地心准惯性坐标系的总能量，在仅考虑两体间万有引力作用下包括动能和势能两部分。对于给定的某一个初始参数，或者任意时刻参数，通过圆轨道高度即可确定航天器运行总能量。航天器相对地心准惯性坐标系的总能量的计算方法包括：

[0082] S10、设任意时刻轨道高度为h，则LEO圆轨道至地心距离为H=h+rE，其中rE为地球平均半径；于是航天器在轨运行时的总能量为E，其计算公式如下：

[0083]

[0084] 其中，mE为地球质量，m为航天器质量，G为万有引力常数，其值为6.67×10‑11N.m2/2
kg；

[0085] S11、所述S10的公式中，航天器动能E动的计算公式如下：

[0086]

[0087] 其中，V为航天器飞行速度，其和航天器运行高度的关系由引力和离心力的平衡关系进行确定；

[0088] S12、所述S10的公式中，航天器势能E势的计算公式如下：

[0089]

[0090] 其中，航天器的势能零点取在无穷远处；

[0091] S13、把航天器初始参数即初始时刻的飞行轨道高度代入S10的公式中即可获得初始时刻航天器的总能量。

[0092] 在上述技术方案中，所述步骤二中，大气阻尼对总能量的耗散的计算方法包括：

[0093] S20、设时间间隔dt内，大气阻尼对总能量的耗散，即气动阻力对航天器做的功WA，其计算公式如下：

[0094]

[0095] 其中，航天器运动轨迹近似为ds=Vdt，其中V为航天器飞行速度；气动阻力做的功主要取决于做功距离ds；实际处理时，时间步长内不考虑气动阻力变化，即时间步长内不考虑相关的阻力系数CD、大气密度ρ、航天器飞行速度V及迎风截面积A的变化；同时，在不同时间过程中大气密度ρ随高度变化采用相关大气模型，航天器飞行速度V随航天器的高度变化根据航天器动能计算，航天器迎风截面积A根据实际情况取合适的量值，即根据航天器实际垂直于风向的净面积而定，必要时可变化，阻力系数CD取定值或根据航天器高度和航天器几何构型和姿态进行精细化标定；

[0096] S21、所述S20的公式中，气动阻力做功的微分形式还表示如下：

[0097]

[0098] 其中，通过换算关系把气动阻力做功的表达式表示为航天器距地心距离H的函数，G为万有引力常数，mE为地球质量；

[0099] S22、所述S20的公式中，FA为气动阻力，表达式为：

[0100]

[0101] 其中，阻力系数CD的参考面积取为单位面积。

[0102] 在上述技术方案中，所述步骤三中，根据耗散后的总能量计算新的圆轨道运行参数的计算方法包括：

[0103] S30、设航天器距地心距离H在时间步长dt前后分别是H1和H2，则根据耗散后总能量计算新的轨道运行参数，其计算公式如下：

[0104]

[0105] 其中，ρ为大气密度，CD为航天器飞行过程中的阻力系数，A为迎风截面积，m为航天器质量，mE为地球质量，G为万有引力常数；

[0106] S31、所述S30的公式中，dH的计算依据即为航天器总能变化dE等于大气阻尼对其耗散，即等于大气阻力对航天器做功WA的负值，具体关系如下：

[0107] dE=‑WA

[0108] 或者

[0109]

[0110] 由此即可得出dH的表达式，CD为阻力系数。

[0111] 在上述技术方案中，所述步骤四中，再入条件设为某个临界高度he，判断是否满足再入条件的方法包括：

[0112] S40、根据S30公式得到的航天器距地心距离得到航天器运行高度h2=H2‑rE，当满足下列条件则说明满足再入条件：

[0113] h2≤he

[0114] 于是进入步骤五；

[0115] S41、根据S30公式得到的航天器距地心距离得到航天器运行高度h2=H2‑rE，当满足下列条件则说明尚未满足再入条件，需要再返回步骤二推进计算：

[0116] h2＞he

[0117] 于是返回步骤二，在给定推进时间步长后继续计算大气阻尼对总能量的耗散，并获得新的轨道参数。

[0118] 在上述技术方案中，所述步骤五中，得到再入点参数及统计计算轨道衰降过程的重要特征参数，它们的计算方法包括：

[0119] S50、根据步骤四设定的再入高度he，即某个临界高度，再入速度Ve的计算公式如下：

[0120]

[0121] 该公式同时适用于任意高度时的飞行速度计算，即 ；

[0122] S51、根据步骤四设定的再入高度he，再入时刻当地弹道倾角de的计算公式如下：

[0123] de=arctan（dhe/Vedte）

[0124] 其中，dte、dhe、Ve分别是再入时刻时的时间步长、高度增量及再入速度；该公式同时适用于任意高度时的当地弹道倾角的计算，即d=arctan（dh/Vdt）；

[0125] S52、根据步骤二至步骤四推进计算过程对时间步长的累加，即可得到从初始状态到再入时的时间预报，时间预报ttotal的计算公式如下：

[0126]

[0127] 根据需要，采用不同时间单位，如“天”、“周”或“年”对时间预报值进行换算。

[0128] 在上述技术方案中，所述临界高度he的高度值在90km至120km之间。

[0129] 实施例：

[0130] 本实施例基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，为了更清楚地说明本发明的技术方案，通过实施例进行说明。实施例针对某大型航天器（天宫一号，简称TG‑1）自主衰降过程，采用本发明描述的方法进行中长期再入参数预报，和实际观测结果进行对比分析，验证了本发明的有效性。

[0131] 表1所示为TG‑1再入之前一段时间，中国载人航天官网发布的部分在轨运行参数。本实施例以此作为初始状态和中间状态对本发明的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法进行实际应用。

[0132] TG‑1离轨及再入时的整器质量取为7660kg。

[0133] 表1官网发布的TG‑1在轨运行参数

[0134]

[0135] 由圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量。

[0136] 以表1中2017年3月13日TG‑1的运行高度作为初始状态，计算得到航天器相对地心准惯性坐标系的总能量以及其它若干关联参数如表2所示。

[0137] 表2 TG‑1在轨运行初始参数

[0138] 参数名称 数值 单位初始日期 2017/03/13 —
初始高度 348.3 km
初始速度 7698 m/s
11
初始总能量 ‑2.2696×10 J
11
初始动能 2.2696×10 J
11
初始势能 ‑4.5392×10 J
14
初始角动量 3.9624×10 N.m.s
初始轨道下降速率 115.035 m/天

[0139] 根据时间步长推进计算飞行过程中大气阻尼对总能量的耗散，根据耗散后总能量计算新的轨道运行参数。

[0140] 根据计算经验，轨道衰降过程参数推进计算时间步长可以从几秒到几百秒不等，选择的主要依据是根据轨道高度动态变化，高度较高时步长较大，高度较低时步长较小。

[0141] 由于航天器轨道衰降运行过程中自由分子流稀薄大气阻尼作用下的航迹计算基于瞬时密切圆轨道参数的弧段解析参数，因此在作用区间内阻尼相对变化极小的情况下，计算精度是能够得到充分保障的。

[0142] 图2和图3所示分别为利用本发明方法计算得到的TG‑1运行高度随时间的变化历程。图2从整体上显示了计算结果和官网发布的实测数据具有良好一致性，图3给出了本发明计算获得的最终再入日期（2017年3月13日起始+382.4天）的预报数据，和实际再入日期（2018年4月2日）基本吻合。

[0143] 需要说明的是，实际计算过程中的航天器迎风截面积根据实际数据进行了反演锚定及修正，这是对于迫近最终再入不断融合测试信息从而不断提高预报精度的必然过程。

[0144] 判断是否满足再入条件，不满足返回步骤二，满足则继续

[0145] 本实施例中，选择高度120km作为进入稠密大气层再入高度。当上述计算结果高度大于120km时则不断推进时间步进行计算；当计算结果小于等于120km时则终止计算，进入下一步骤。

[0146] 表3所示为采用本发明方法计算得到的TG‑1起始高度与衰降至120km高度时所需时间关系。

[0147] 表3 TG‑1起始高度与衰降至120km高度所需时间关系

[0148]

[0149] 图4‑图6分别给出了根据本发明方法计算得到的轨道衰降过程中航天器的飞行速度、飞行高度、当地弹道倾角、每日轨道高度衰降值随时间的变化历程，这些参数同时也是随轨道高度相关而变化的。从图4可知，航天器轨道衰降过程中，其飞行速度是逐渐有所增大的（在再入稠密大气层之前），其原理在于部分势能转化为动能的结果，且大气阻尼对总能的耗散减量有限，还没有达到能够同时减小动能的程度。从图5可知，航天器在飞行高度约150km以上时，其当地弹道倾角和变化极小且接近0度，当高度低于120km时才有明显减小变化；从图6及计算详细数据结果可知，航天器轨道衰降每日降低高度随高度降低而增大，在飞行高度约150km以上时，轨道每日降低高度较小，在百米量级，高度低于150km时每日降低高度明显增大，当高度低于120km时每日高度降低急剧增大。

[0150] 得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论[0151] 根据再入高度120km对应的时间推算计算结果，获得再入点参数并统计计算轨道衰降过程重要特征参数如表4所示。

[0152] 表4 TG‑1再入参数

[0153] 参数名称 数值 单位再入日期 2018/04/02 —
再入高度 120.0 km
再入速度 7832 m/s
再入弹道倾角 ‑0.031 deg.
11
再入总能量 ‑2.3494×10 J
11
再入动能 2.3494×10 J
11
再入势能 ‑4.6988×10 J
14
再入角动量 3.8945×10 N.m.s

[0154] 上述参数中，对再入日期/时间的中长期快速预报是最重要、最根本的一个方面。

[0155] 这里说明的处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。

[0156] 尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。