CN110598270A - 一种基于编目根数序列的高精度空间目标陨落预报方法 - Google Patents

Abstract

本发明公开了一种基于编目根数序列的高精度空间目标陨落预报方法，本发明具体包括以下步骤：步骤一、解析编目根数系列作为实测轨道值，步骤二、建立轨道摄动模型，步骤三、计算弹道系数，步骤四、预报陨落时间和落点位置，本发明涉及航天测量与控制技术领域。该基于编目根数序列的高精度空间目标陨落预报方法。可实现通过使用半数值法进行轨道外推，积分步长可以取为轨道周期的整数倍，为几十分钟甚至几十个小时，与数值法相比计算效率更高，同时通过使用半数值法进行轨道积分，很好的保证了预报精度的同时提高了计算效率，可同时满足大量空间目标情况下高精度和高效率的计算需求，从而实现既快速又准确地进行陨落预报。

Description

一种基于编目根数序列的高精度空间目标陨落预报方法

技术领域

本发明涉及航天测量与控制技术领域，具体为一种基于编目根数序列的高精度空间目标陨落预报方法。

背景技术

自1957年人类进行第一次太空活动以来，有超过30,000吨的太空垃圾返回地球，随着人类空间活动的日益频繁、空间目标在轨飞行时间不断积累，陨落的空间目标数量也在不断增长，大质量空间目标在陨落过程中并不会被完全烧毁，仍有10-40％的残骸返回地球表面，对地表的生命群体、建筑设施、生态环境等具有很大威胁，因此，对即将陨落的空间目标进行密集监测、精确预报陨落时间和落点，对于掌握太空态势变化、及时规避损害等具有重要意义，进行空间目标陨落预报计算，最好的方法是在获取大量测量数据的基础上确定精确轨道并解算弹道系数，但是只有少数国家具有获取再入目标测量数据的能力,大部分只能依靠公开的TLE(Two Line Element，两行根数)进行预报，目前，对于如何使用TLE进行陨落预报主要从两个方面进行研究，一是如何基于TLE确定空间目标的精密轨道，二是怎样解析出大气阻力系数、光压系数等空间目标的关键参数。

大气阻力是即将陨落空间目标所受的最主要非保守摄动力。进行精确陨落预报的主要难点在于轨道确定和大气阻力建模，最优的方法是在获取大量测量数据的基础上确定精确轨道并解算弹道系数。但是目前只有少数国家具有获取陨落目标测量数据的能力，大部分只能依靠公开的TLE(Two Line Element，两行根数)进行预报。基于TLE进行陨落预报最简单的方法是通过调整空间目标的弹道系数，使较老TLE预报至较新TLE历元时预报的半长轴与较新TLE的实际半长轴一致，这种方法的缺点是预报结果的精度受TLE精度的影响明显。为减小TLE带来的误差，TLE使用前需先进行过滤筛选，剔除异常值后对多组TLE使用最小二乘法拟合弹道系数，或者根据多组TLE生成伪测量数据，重新确定轨道并解算弹道系数。从公开资料看，当前陨落预报的误差基本在10～30％之间，大气环境异常或者空间目标结构发生较大变化的情况下误差可能超过30％。

传统的基于编目根数(包括TLE)的空间目标陨落预报方法一般使用数值法进行计算，由于可以使用高精度的三维大气密度模型，相较于分析法，数值法精度更高，但是计算效率低，难以适用于大批量空间目标的情况。

发明内容

(一)解决的技术问题

针对现有技术的不足，本发明提供了一种基于编目根数序列的高精度空间目标陨落预报方法，解决了现有分析法预报精度不高、数值法计算效率低的问题。

(二)技术方案

为实现以上目的，本发明通过以下技术方案予以实现：一种基于编目根数序列的高精度空间目标陨落预报方法，具体包括以下步骤：

步骤一、解析编目根数系列作为实测轨道值；

步骤二、建立轨道摄动模型：考虑的摄动力模型包括地球非球形摄动(J₂、J₃项)、大气阻力摄动及日月引力摄动，摄动方程为：

步骤三、计算弹道系数：其具体包括以下步骤：

A1、设置初始弹道系数B₀＝0.02m²/kg；

A2、以第5条编目根数和B₀为初始值，1轨道周期为步长对步骤二的摄动方程公式进行反向积分；

A3、若积分结束时刻与第i条根数历元t_i之差小于半轨道周期，则计算半长轴计算值与实测值之差Δa_i以及

A4、第1条编目根数计算结束后，计算弹道系数的修正值ΔB；

A5、若则计算结束，否则取B₀＝B₀+ΔB，返回步骤A2重新计算，ε可取为0.00001；

步骤四、预报陨落时间和落点位置。

优选的，所述步骤一解析编目根数系列作为实测轨道值的具体步骤为：

C1、获取空间目标的5条编目根数，并按根数历元排序，要求相邻两条编目根数历元的间隔不小于一轨道周期；

C2、解析编目根数序列，并转换为历元时刻的平根数，结果为

优选的，所述步骤二中，(a,e,i,Ω,ω,M)为开普勒轨道根数，均使用平根数，∑S_i、∑T_i、∑W_i分别表示各摄动力的S、T、W方向分量之和。

优选的，所述步骤A3中ρ为大气密度值，V为空间目标相对大气的运动速度，f为真近点角，ω_E为地球运动速率，μ为地球引力常数，r为空间目标的地心向距。

优选的，所述步骤A4中弹道系数修正值ΔB的计算公式为：

优选的，所述步骤四中预报陨落时间和落点位置的方法，具体包括以下步骤：

D1、设置积分步长控制器λ＝1；

D2、以第5条编目根数和弹道系数结果B为初值，λ倍轨道周期为步长，对摄动方程公式进行积分，记积分开始时刻为t_s，则积分结束时刻为t_s+λT；

D3、若t_s+λT时刻空间目标的近地点高度大于80公里，设置积分开始时刻t_s＝t_s+λT，根据公式调整积分步长控制器，以λT为步长继续积分；

D4、若t_s+λT时刻空间目标的近地点高度小于80公里，且λ>1，则取λ＝λ/2，步长为λT，从t_s时刻重新积分；

D5、若t_s+λT时刻空间目标的近地点高度小于80公里，且λ＝1，则积分结束，t_s即为空间目标的陨落时间，对应的大地经度、大地纬度为落点位置。

(三)有益效果

本发明提供了一种基于编目根数序列的高精度空间目标陨落预报方法。与现有技术相比具备以下有益效果：该基于编目根数序列的高精度空间目标陨落预报方法，具体包括以下步骤：步骤一、解析编目根数系列作为实测轨道值，步骤二、建立轨道摄动模型：考虑的摄动力模型包括地球非球形摄动(J₂、J₃项)、大气阻力摄动及日月引力摄动，步骤三、计算弹道系数，步骤四、预报陨落时间和落点位置，可实现通过使用半数值法进行轨道外推，积分步长可以取为轨道周期的整数倍，为几十分钟甚至几十个小时，与数值法相比计算效率更高，同时通过使用半数值法进行轨道积分，很好的保证了预报精度的同时提高了计算效率，可同时满足大量空间目标情况下高精度和高效率的计算需求，从而实现既快速又准确地进行陨落预报。

附图说明

图1为本发明的工作流程图；

图2为本发明38852的编目根数序列图；

图3为本发明38852的平根数结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1-3，本发明实施例提供一种技术方案：一种基于编目根数序列的高精度空间目标陨落预报方法，以NORAD编号为38852的空间目标为例，使用美国战略司令部于space-track网站发布的TLE作为编目根数。实施步骤如下：

步骤一、解析编目根数系列作为实测轨道值。

C1、获取38852目标的5条TLE，并按根数历元排序，如图2所示；

C2、解析出的38852目标的历元、平半长轴、平偏心率，如图3所示。

步骤二、建立轨道摄动模型：考虑的摄动力模型包括地球非球形摄动(J₂、J₃项)、大气阻力摄动及日月引力摄动，摄动方程如下：

其中，(a,e,i,Ω,ω,M)为开普勒轨道根数，均使用平根数。∑S_i、∑T_i、∑W_i分别表示各摄动力的S、T、W方向分量之和。

本发明使用MSIS-90大气密度模型计算大气密度。

步骤三、计算弹道系数：其具体包括以下步骤：

A1、设置初始弹道系数B₀＝0.02m²/kg；

A2、以第5条编目根数和B₀为初始值，1轨道周期为步长对步骤二的摄动方程公式进行反向积分；

A3、若积分结束时刻与第i条根数历元t_i之差小于半轨道周期，则计算半长轴计算值与实测值之差Δa_i以及 其中，ρ为大气密度值，V为空间目标相对大气的运动速度，f为真近点角，ω_E为地球运动速率，μ为地球引力常数，r为空间目标的地心向距；

A4、第1条编目根数计算结束后，计算弹道系数的修正值ΔB，

A5、若则计算结束；否则取B₀＝B₀+ΔB，返回步骤A2重新计算，ε可取为0.00001。

计算共迭代4次，计算出的38852目标的弹道系数结果为0.030148m2/kg。

步骤四：预报陨落时间和落点位置的方法，具体包括以下步骤：

D1、设置积分步长控制器λ＝1；

D2、以第5条编目根数和弹道系数结果B为初值，λ倍轨道周期为步长，对摄动方程公式进行积分；

D3、记积分开始时刻为t_s，则积分结束时刻为t_s+λT。若t_s+λT时刻空间目标的近地点高度大于80公里，设置积分开始时刻t_s＝t_s+λT，根据公式设置积分步长控制器，以λT为步长继续积分；

D4、若t_s+λT时刻空间目标的近地点高度小于80公里，且λ>1，则取λ＝λ/2，步长为λT，从t_s时刻重新积分；

D5、若t_s+λT时刻空间目标的近地点高度小于80公里，且λ＝1，则积分结束，t_s即为空间目标的陨落时间，对应的大地经度、大地纬度为落点位置。

计算出的38852目标的陨落时间为UTC时2013年8年5月14时4分，落点位于(98.28°E，47.26°N)。

1、计算精度

使用陨落时间误差百分比δ进行误差评估，δ用下式计算：其中，t_real为目标的真实陨落时间，t_pred为预报的陨落时间，t₀为进行陨落预报时使用轨道的历元。

对38852目标而言，根据space-track网站发布的结果，38852的陨落时间为2013年8年6月3时39分，经分析，此结果较为接近实际陨落时间，可作为t_real；预报的陨落时间t_pred为2013年8年5月14时4分；预报时从第5条根数开始分析，即t₀为2013年7年26月13时14分，根据公式计算出陨落时间误差百分比为5.3％，低于10％至30％的正常水平，预报精度很高。

同样以日本FITSAT 1卫星为例进行精度验证，按轨道寿命长短进行6次计算，精度结果如表1所示，可以看出误差基本小于10％，低于10％至30％的正常水平。

表1 FITSAT 1卫星的陨落预报精度

2、计算效率

随机选取100个近地空间目标作为样本，在Intel酷睿2处理器(2.7G主频，1G内存)，Windows XP操作系统计算环境下，计算弹道系数时平均迭代次数约为5.6次，平均耗时小于3秒，进行陨落时间预报时预报100天耗时约5秒，满足大批量目标的实际应用要求。

综上所述

本发明可实现通过使用半数值法进行轨道外推，积分步长可以取为轨道周期的整数倍，为几十分钟甚至几十个小时，与数值法相比计算效率更高，同时通过使用半数值法进行轨道积分，很好的保证了预报精度的同时提高了计算效率，可同时满足大量空间目标情况下高精度和高效率的计算需求，从而实现既快速又准确地进行陨落预报。

需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (6)

1.一种基于编目根数序列的高精度空间目标陨落预报方法，其特征在于：具体包括以下步骤：

步骤一、解析编目根数系列作为实测轨道值；

步骤二、建立轨道摄动模型：考虑的摄动力模型包括地球非球形摄动(J₂、J₃项)、大气阻力摄动及日月引力摄动，摄动方程为：

步骤三、计算弹道系数：其具体包括以下步骤：

A1、设置初始弹道系数B₀＝0.02m²/kg；

A2、以第5条编目根数和B₀为初始值，1轨道周期为步长对步骤二的摄动方程公式进行反向积分；

A3、若积分结束时刻与第i条根数历元t_i之差小于半轨道周期，则计算半长轴计算值与实测值之差Δa_i以及

A4、第1条编目根数计算结束后，计算弹道系数的修正值ΔB；

A5、若则计算结束，否则取B₀＝B₀+ΔB，返回步骤A2重新计算，ε可取为0.00001；

步骤四、预报陨落时间和落点位置。

2.根据权利要求1所述的一种基于编目根数序列的高精度空间目标陨落预报方法，其特征在于：所述步骤一解析编目根数系列作为实测轨道值的具体步骤为：

C1、获取空间目标的5条编目根数，并按根数历元排序，要求相邻两条编目根数历元的间隔不小于一轨道周期；

C2、解析编目根数序列，并转换为历元时刻的平根数，结果为

3.根据权利要求1所述的一种基于编目根数序列的高精度空间目标陨落预报方法，其特征在于：所述步骤二中，(a,e,i,Ω,ω,M)为开普勒轨道根数，均使用平根数，∑S_i、∑T_i、∑W_i分别表示各摄动力的S、T、W方向分量之和。

4.根据权利要求1所述的一种基于编目根数序列的高精度空间目标陨落预报方法，其特征在于：所述步骤A3中ρ为大气密度值，V为空间目标相对大气的运动速度，f为真近点角，ω_E为地球运动速率，μ为地球引力常数，r为空间目标的地心向距。

5.根据权利要求1所述的一种基于编目根数序列的高精度空间目标陨落预报方法，其特征在于：所述步骤A4中弹道系数修正值ΔB的计算公式为：

6.根据权利要求1所述的一种基于编目根数序列的高精度空间目标陨落预报方法，其特征在于：所述步骤四中预报陨落时间和落点位置的方法，具体包括以下步骤：

D1、设置积分步长控制器λ＝1；

D2、以第5条编目根数和弹道系数结果B为初值，λ倍轨道周期为步长，对摄动方程公式进行积分，记积分开始时刻为t_s，则积分结束时刻为t_s+λT；

D3、若t_s+λT时刻空间目标的近地点高度大于80公里，设置积分开始时刻t_s＝t_s+λT，根据公式调整积分步长控制器，以λT为步长继续积分；

D4、若t_s+λT时刻空间目标的近地点高度小于80公里，且λ>1，则取λ＝λ/2，步长为λT，从t_s时刻重新积分；

D5、若t_s+λT时刻空间目标的近地点高度小于80公里，且λ＝1，则积分结束，t_s即为空间目标的陨落时间，对应的大地经度、大地纬度为落点位置。