CN116384599A - 基于能量分析的航天器leo圆轨道衰降过程参数预报方法 - Google Patents

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

Info

Publication number
CN116384599A
CN116384599A CN202310660927.7A CN202310660927A CN116384599A CN 116384599 A CN116384599 A CN 116384599A CN 202310660927 A CN202310660927 A CN 202310660927A CN 116384599 A CN116384599 A CN 116384599A
Authority
CN
China
Prior art keywords
spacecraft
reentry
parameters
orbit
formula
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202310660927.7A
Other languages
English (en)
Other versions
CN116384599B (zh
Inventor
唐小伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Original Assignee
Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center filed Critical Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Priority to CN202310660927.7A priority Critical patent/CN116384599B/zh
Publication of CN116384599A publication Critical patent/CN116384599A/zh
Application granted granted Critical
Publication of CN116384599B publication Critical patent/CN116384599B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F7/00—Methods or arrangements for processing data by operating upon the order or content of the data handled
    • G06F7/38—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation
    • G06F7/48—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation using non-contact-making devices, e.g. tube, solid state device; using unspecified devices
    • G06F7/544—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation using non-contact-making devices, e.g. tube, solid state device; using unspecified devices for evaluating functions by calculation
    • G06F7/548—Trigonometric functions; Co-ordinate transformations
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Software Systems (AREA)
  • Strategic Management (AREA)
  • Economics (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Operations Research (AREA)
  • Human Resources & Organizations (AREA)
  • Marketing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Business, Economics & Management (AREA)
  • Game Theory and Decision Science (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • Development Economics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

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

Description

基于能量分析的航天器LEO圆轨道衰降过程参数预报方法
技术领域
本发明属于航天工程及轨道动力学、航天器速度测定技术领域，更具体地说，本发明涉及基于能量分析的航天器LEO圆轨道衰降过程参数预报方法。
背景技术
服役期满非回收航天器，包括各种人造卫星、空间实验室、空间站及运载火箭末级等，一般都要经历轨道衰降过程后再入稠密大气层陨落解体。近年来，空间碎片问题引起普遍关注，各国科研人员一直在对航天器再入的中长期预报进行研究，准确可靠的再入预报能够提醒对地面人员设施的安全预警，及时回应国际关切和舆论。对轨道衰降阶段的预报是落区预报的前提，准确可靠的再入点参数对后续陨落解体及地面风险评估具有重要影响。但是，航天器在轨飞行的空间环境非常复杂且多变，其运动轨迹受到的影响因素包括地球/月球/太阳引力、大气密度模式、重力梯度矩、太阳辐射压力、地磁作用、地球非质点非对称质量分布等等；因此，要精确求解航天器在轨运行轨迹的相关参数需要复杂的轨道动力学理论及求解复杂的运动学和动力学方程组，同时需要对预期的一些影响因素进行研究和尽可能可靠的设定。
目前，轨道参数的定解计算常用摄动方法，包括各种坐标摄动法、瞬时椭圆法和正则变换；摄动方程求解方法包括积分法、解析法和半解析法。无论那种摄动方法，都涉及摄动方程建模和求解，模型的精细程度取决于摄动小量的展开级数表达式及其展开阶数；数学上模型的精细化和复杂性将带来极大计算代价。用数值法求解轨道寿命，传统方法诸如考虑大气摄动加速度的科威尔（Cowell）法、恩克（Encke）法或轨道要素变动法，这些方法由于时间步长不能取得太大，当轨道高度较高因而衰减非常缓慢时，计算需要非常长的时间，甚至难以实现。现有的解析法虽然计算简便快捷，但仅适用于人为假设的指数型大气模型及一些特定简化模式，不能普遍使用。
航天器轨道衰降问题即空间碎片减缓的主要基本问题。空间碎片减缓的主动措施之一是主动变轨快速再入，其次就是针对一些卫星要求在使命完成后进行离轨机动，进入一条短寿命（通常要求在25年内）的轨道，以免其长期滞留太空，污染太空环境。另一方面也有一些航天器丧失了变轨机动能力。上述情况都需要对其自主衰降过程进行预报。
针对LEO圆轨道衰降参数的计算分析预报，关键的是一个目标参数即飞行高度的变化及衰降计算，并由此获得对再入时间和部分重要参数进行中长期预报的结论。LEO指的是低地球轨道，一般是高度低于2000km的轨道。LEO圆轨道，是若干人造卫星或空间实验室、空间站等航天器绕地球运行飞行时常见的一种轨道模式，其高度在一个运行周期内变化极小。当实际轨道高度在一个运行周期内变化范围不超过50km时，可以近似认为是圆轨道。经验表明，轨道自然衰降再入中长期预报时，即使仅提前几天，要对再入时刻位置（空间位置或星下点）进行有实用意义的确定也几乎是不可能的。由于再入位置对各种影响因素极为敏感，因而可以把分析重点放在对再入日期的高效计算预报方面。因此，专注于对轨道高度直接相关参数的精细分析，简化或忽略轨道倾角、升交点赤经等参数，从能量分析角度构建一种LEO圆轨道衰降参数的中长期预报快速高效方法即成为本发明的主要内容。
发明内容
本发明的一个目的是解决上述问题和/或缺陷，并提供后面将说明的优点。
为了实现根据本发明的这些目的和优点，提供了基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，包括以下步骤：
步骤一、由LEO圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量，其中总能量包括动能和势能；
步骤二、根据时间步长推进计算飞行过程中大气阻尼对动能的耗散，从而得到对总能量的耗散；
步骤三、根据耗散后的总能量计算新的圆轨道运行参数，包括圆轨道半径或高度；
步骤四、根据轨道高度判断是否满足再入条件，不满足返回步骤二，满足则继续下一步骤；
步骤五、得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论，包括自初始时刻至再入时刻的时间预报。
优选的是，其中，所述步骤一中，航天器相对地心准惯性坐标系的总能量计算方法包括：
S10、设任意时刻轨道高度为h，则LEO圆轨道至地心距离为H=h+r E ，其中r E 为地球平均半径；于是航天器在轨运行时的总能量为E，其计算公式如下：
Figure SMS_1
其中，m E 为地球质量，m为航天器质量，G为万有引力常数，其值为6.67×10-11（N.m2/kg2）；
S11、所述S10的公式中，航天器动能E 动的计算公式如下：
Figure SMS_2
其中，V为航天器飞行速度，其和航天器运行高度的关系由引力和离心力的平衡关系进行确定；
S12、所述S10的公式中，航天器势能E 势的计算公式如下：
Figure SMS_3
其中，航天器的势能零点取在无穷远处；
S13、把航天器初始参数即初始时刻的飞行轨道高度代入S10的公式中即可获得初始时刻航天器的总能量。
优选的是，其中，所述步骤二中，大气阻尼对总能量的耗散的计算方法包括：
S20、S20、设时间间隔dt内，大气阻尼对总能量的耗散，即气动阻力对航天器做的功W A ，其计算公式如下：
Figure SMS_4
其中，航天器运动轨迹近似为ds=Vdt，其中V为航天器飞行速度；气动阻力做的功主要取决于做功距离ds；实际处理时，时间步长内不考虑气动阻力变化，即时间步长内不考虑相关的阻力系数C D 、大气密度ρ、航天器飞行速度V及迎风截面积A的变化；同时，在不同时间过程中大气密度ρ随高度变化采用相关大气模型，航天器飞行速度V随航天器的高度变化根据航天器动能计算，航天器迎风截面积A根据实际情况取合适的量值，即根据航天器实际垂直于风向的净面积而定，必要时可变化，阻力系数C D 取定值或根据航天器高度和航天器几何构型和姿态进行精细化标定；
S21、所述S20的公式中，气动阻力做功的微分形式还表示如下：
Figure SMS_5
其中，通过换算关系把气动阻力做功的表达式表示为航天器距地心距离H的函数，G为万有引力常数，m E 为地球质量；
S22、所述S20的公式中，F A 为气动阻力，表达式为：
Figure SMS_6
其中，阻力系数C D 的参考面积取为单位面积。
优选的是，其中，所述步骤三中，根据耗散后的总能量计算新的圆轨道运行参数的计算方法包括：
S30、设航天器距地心距离H在时间步长dt前后分别是H 1和H 2，则根据耗散后总能量计算新的轨道运行参数，其计算公式如下：
Figure SMS_7
其中，ρ为大气密度，C D 为航天器飞行过程中的阻力系数，A为迎风截面积，m为航天器质量，m E 为地球质量，G为万有引力常数；
S31、所述S30的公式中，dH的计算依据即为航天器总能变化dE等于大气阻尼对其耗散，即等于大气阻力对航天器做功W A 的负值，具体关系如下：
dE=-W A
或者
Figure SMS_8
由此即可得出dH的表达式，C D 为阻力系数。
优选的是，其中，所述步骤四中，再入条件设为某个临界高度h e ，判断是否满足再入条件的方法包括：
S40、根据S30公式得到的航天器距地心距离得到航天器运行高度h 2=H 2-r E，当满足下列条件则说明满足再入条件：
h 2≤h e
于是进入步骤五；
S41、根据S30公式得到的航天器距地心距离得到航天器运行高度h 2=H 2-r E，当满足下列条件则说明尚未满足再入条件，需要再返回步骤二推进计算：
h 2＞h e
于是返回步骤二，在给定推进时间步长后继续计算大气阻尼对总能量的耗散，并获得新的轨道参数。
优选的是，其中，所述步骤五中，得到再入点参数及统计计算轨道衰降过程的重要特征参数，它们的计算方法包括：
S50、根据步骤四设定的再入高度h e ，即某个临界高度，再入速度V e 的计算公式如下：
Figure SMS_9
该公式同时适用于任意高度时的飞行速度计算，即
Figure SMS_10
；
S51、根据步骤四设定的再入高度h e ，再入时刻当地弹道倾角d e 的计算公式如下：
d e =arctan（dh e /V e dt e ）
其中，dt e 、dh e 、V e 分别是再入时刻时的时间步长、高度增量及再入速度；该公式同时适用于任意高度时的当地弹道倾角的计算，即d=arctan（dh/Vdt）；
S52、根据步骤二至步骤四推进计算过程对时间步长的累加，即可得到从初始状态到再入时的时间预报，时间预报t total 的计算公式如下：
Figure SMS_11
根据需要，采用不同时间单位，如“天”、“周”或“年”对时间预报值进行换算。
优选的是，其中，所述临界高度h e 的高度值在90km至120km之间。
本发明至少包括以下有益效果：
其一，本发明通过简单的积分方法，即可快速高效进行航天器LEO圆轨道衰降过程的主要参数进行计算，从而提供了一种简洁的航天器再入陨落时间的中长期预报方法。
其二，本发明构造了LEO圆轨道衰降过程高度变化的算法流程，从而在此基础上推算速度、倾角等其它参数，也可以推算出更多的导出参数，包括飞行角速度、每天衰降高度、累计飞行圈数、飞行航迹航程、累计飞行时间等等。
其三，本发明的LEO圆轨道衰降过程参数预报方法，仅存在一个难以确定的参数，即航天器飞行时的迎风截面积。对于复杂外形及姿态变化情况下的航天器，该面积参数可能存在较大不确定性，但仍然可以通过一定的观测数据对该迎风面积进行动态标定。
其四，本发明基于能量分析的航天器LEO圆轨道衰降过程参数预报方法是基于物理定律得到的应用手段，其预报结果的准确性还取决于大气密度模式，为此建议采用较为精细的大气模式算法。在这些方面有若干成熟的研究成果能够直接借鉴
本发明的优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。
附图说明
图1为航天器LEO圆轨道衰降过程参数分析示意图；
图2为TG-1轨道衰降过程参数预报结果和官方发布参数；
图3为TG-1轨道衰降过程参数再入前最后几天高度结果；
图4为TG-1轨道衰降过程飞行速度及高度变化历程；
图5为TG-1轨道衰降过程飞行当地倾角变化历程；
图6为TG-1轨道衰降过程飞行每日高度衰降变化历程。
具体实施方式
下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。
应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不排除一个或多个元件或其组合的存在或添加。
如图1所示：本发明的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，包括以下步骤：
步骤一、由LEO圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量，其中总能量包括动能和势能；
步骤二、根据时间步长推进计算飞行过程中大气阻尼对动能的耗散，从而得到对总能量的耗散；
步骤三、根据耗散后的总能量计算新的圆轨道运行参数，包括圆轨道半径或高度；
步骤四、根据轨道高度判断是否满足再入条件，不满足返回步骤二，满足则继续下一步骤；
步骤五、得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论，包括自初始时刻至再入时刻的时间预报。
在上述技术方案中，所述步骤一中，绕地运行航天器相对地心准惯性坐标系的总能量，在仅考虑两体间万有引力作用下包括动能和势能两部分。对于给定的某一个初始参数，或者任意时刻参数，通过圆轨道高度即可确定航天器运行总能量。航天器相对地心准惯性坐标系的总能量的计算方法包括：
S10、设任意时刻轨道高度为h，则LEO圆轨道至地心距离为H=h+r E ，其中r E 为地球平均半径；于是航天器在轨运行时的总能量为E，其计算公式如下：
Figure SMS_12
其中，m E 为地球质量，m为航天器质量，G为万有引力常数，其值为6.67×10-11N.m2/kg2；
S11、所述S10的公式中，航天器动能E 动的计算公式如下：
Figure SMS_13
其中，V为航天器飞行速度，其和航天器运行高度的关系由引力和离心力的平衡关系进行确定；
S12、所述S10的公式中，航天器势能E 势的计算公式如下：
Figure SMS_14
其中，航天器的势能零点取在无穷远处；
S13、把航天器初始参数即初始时刻的飞行轨道高度代入S10的公式中即可获得初始时刻航天器的总能量。
在上述技术方案中，所述步骤二中，大气阻尼对总能量的耗散的计算方法包括：
S20、设时间间隔dt内，大气阻尼对总能量的耗散，即气动阻力对航天器做的功W A ，其计算公式如下：
Figure SMS_15
其中，航天器运动轨迹近似为ds=Vdt，其中V为航天器飞行速度；气动阻力做的功主要取决于做功距离ds；实际处理时，时间步长内不考虑气动阻力变化，即时间步长内不考虑相关的阻力系数C D 、大气密度ρ、航天器飞行速度V及迎风截面积A的变化；同时，在不同时间过程中大气密度ρ随高度变化采用相关大气模型，航天器飞行速度V随航天器的高度变化根据航天器动能计算，航天器迎风截面积A根据实际情况取合适的量值，即根据航天器实际垂直于风向的净面积而定，必要时可变化，阻力系数C D 取定值或根据航天器高度和航天器几何构型和姿态进行精细化标定；
S21、所述S20的公式中，气动阻力做功的微分形式还表示如下：
Figure SMS_16
其中，通过换算关系把气动阻力做功的表达式表示为航天器距地心距离H的函数，G为万有引力常数，m E 为地球质量；
S22、所述S20的公式中，F A 为气动阻力，表达式为：
Figure SMS_17
其中，阻力系数C D 的参考面积取为单位面积。
在上述技术方案中，所述步骤三中，根据耗散后的总能量计算新的圆轨道运行参数的计算方法包括：
S30、设航天器距地心距离H在时间步长dt前后分别是H 1和H 2，则根据耗散后总能量计算新的轨道运行参数，其计算公式如下：
Figure SMS_18
其中，ρ为大气密度，C D 为航天器飞行过程中的阻力系数，A为迎风截面积，m为航天器质量，m E 为地球质量，G为万有引力常数；
S31、所述S30的公式中，dH的计算依据即为航天器总能变化dE等于大气阻尼对其耗散，即等于大气阻力对航天器做功W A 的负值，具体关系如下：
dE=-W A
或者
Figure SMS_19
由此即可得出dH的表达式，C D 为阻力系数。
在上述技术方案中，所述步骤四中，再入条件设为某个临界高度h e ，判断是否满足再入条件的方法包括：
S40、根据S30公式得到的航天器距地心距离得到航天器运行高度h 2=H 2-r E，当满足下列条件则说明满足再入条件：
h 2≤h e
于是进入步骤五；
S41、根据S30公式得到的航天器距地心距离得到航天器运行高度h 2=H 2-r E，当满足下列条件则说明尚未满足再入条件，需要再返回步骤二推进计算：
h 2＞h e
于是返回步骤二，在给定推进时间步长后继续计算大气阻尼对总能量的耗散，并获得新的轨道参数。
在上述技术方案中，所述步骤五中，得到再入点参数及统计计算轨道衰降过程的重要特征参数，它们的计算方法包括：
S50、根据步骤四设定的再入高度h e ，即某个临界高度，再入速度V e 的计算公式如下：
Figure SMS_20
该公式同时适用于任意高度时的飞行速度计算，即
Figure SMS_21
；
S51、根据步骤四设定的再入高度h e ，再入时刻当地弹道倾角d e 的计算公式如下：
d e =arctan（dh e /V e dt e ）
其中，dt e 、dh e 、V e 分别是再入时刻时的时间步长、高度增量及再入速度；该公式同时适用于任意高度时的当地弹道倾角的计算，即d=arctan（dh/Vdt）；
S52、根据步骤二至步骤四推进计算过程对时间步长的累加，即可得到从初始状态到再入时的时间预报，时间预报t total 的计算公式如下：
Figure SMS_22
根据需要，采用不同时间单位，如“天”、“周”或“年”对时间预报值进行换算。
在上述技术方案中，所述临界高度h e 的高度值在90km至120km之间。
实施例：
本实施例基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，为了更清楚地说明本发明的技术方案，通过实施例进行说明。实施例针对某大型航天器（天宫一号，简称TG-1）自主衰降过程，采用本发明描述的方法进行中长期再入参数预报，和实际观测结果进行对比分析，验证了本发明的有效性。
表1所示为TG-1再入之前一段时间，中国载人航天官网发布的部分在轨运行参数。本实施例以此作为初始状态和中间状态对本发明的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法进行实际应用。
TG-1离轨及再入时的整器质量取为7660kg。
表1官网发布的TG-1在轨运行参数
Figure SMS_23
由圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量。
以表1中2017年3月13日TG-1的运行高度作为初始状态，计算得到航天器相对地心准惯性坐标系的总能量以及其它若干关联参数如表2所示。
表2 TG-1在轨运行初始参数
参数名称 数值 单位
初始日期 2017/03/13 —
初始高度 348.3 km
初始速度 7698 m/s
初始总能量 -2.2696×1011 J
初始动能 2.2696×1011 J
初始势能 -4.5392×1011 J
初始角动量 3.9624×1014 N.m.s
初始轨道下降速率 115.035 m/天
根据时间步长推进计算飞行过程中大气阻尼对总能量的耗散，根据耗散后总能量计算新的轨道运行参数。
根据计算经验，轨道衰降过程参数推进计算时间步长可以从几秒到几百秒不等，选择的主要依据是根据轨道高度动态变化，高度较高时步长较大，高度较低时步长较小。
由于航天器轨道衰降运行过程中自由分子流稀薄大气阻尼作用下的航迹计算基于瞬时密切圆轨道参数的弧段解析参数，因此在作用区间内阻尼相对变化极小的情况下，计算精度是能够得到充分保障的。
图2和图3所示分别为利用本发明方法计算得到的TG-1运行高度随时间的变化历程。图2从整体上显示了计算结果和官网发布的实测数据具有良好一致性，图3给出了本发明计算获得的最终再入日期（2017年3月13日起始+382.4天）的预报数据，和实际再入日期（2018年4月2日）基本吻合。
需要说明的是，实际计算过程中的航天器迎风截面积根据实际数据进行了反演锚定及修正，这是对于迫近最终再入不断融合测试信息从而不断提高预报精度的必然过程。
判断是否满足再入条件，不满足返回步骤二，满足则继续
本实施例中，选择高度120km作为进入稠密大气层再入高度。当上述计算结果高度大于120km时则不断推进时间步进行计算；当计算结果小于等于120km时则终止计算，进入下一步骤。
表3所示为采用本发明方法计算得到的TG-1起始高度与衰降至120km高度时所需时间关系。
表3 TG-1起始高度与衰降至120km高度所需时间关系
Figure SMS_24
图4-图6分别给出了根据本发明方法计算得到的轨道衰降过程中航天器的飞行速度、飞行高度、当地弹道倾角、每日轨道高度衰降值随时间的变化历程，这些参数同时也是随轨道高度相关而变化的。从图4可知，航天器轨道衰降过程中，其飞行速度是逐渐有所增大的（在再入稠密大气层之前），其原理在于部分势能转化为动能的结果，且大气阻尼对总能的耗散减量有限，还没有达到能够同时减小动能的程度。从图5可知，航天器在飞行高度约150km以上时，其当地弹道倾角和变化极小且接近0度，当高度低于120km时才有明显减小变化；从图6及计算详细数据结果可知，航天器轨道衰降每日降低高度随高度降低而增大，在飞行高度约150km以上时，轨道每日降低高度较小，在百米量级，高度低于150km时每日降低高度明显增大，当高度低于120km时每日高度降低急剧增大。
得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论
根据再入高度120km对应的时间推算计算结果，获得再入点参数并统计计算轨道衰降过程重要特征参数如表4所示。
表4 TG-1再入参数
参数名称 数值 单位
再入日期 2018/04/02 —
再入高度 120.0 km
再入速度 7832 m/s
再入弹道倾角 -0.031 deg.
再入总能量 -2.3494×1011 J
再入动能 2.3494×1011 J
再入势能 -4.6988×1011 J
再入角动量 3.8945×1014 N.m.s
上述参数中，对再入日期/时间的中长期快速预报是最重要、最根本的一个方面。
这里说明的处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。
尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (7)

1.基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，包括以下步骤：
步骤一、由LEO圆轨道初始参数计算航天器相对地心准惯性坐标系的总能量，其中总能量包括动能和势能；
步骤二、根据时间步长推进计算飞行过程中大气阻尼对动能的耗散，从而得到对总能量的耗散；
步骤三、根据耗散后的总能量计算新的圆轨道运行参数，包括圆轨道半径或高度；
步骤四、根据轨道高度判断是否满足再入条件，不满足返回步骤二，满足则继续下一步骤；
步骤五、得到再入点参数并统计计算轨道衰降过程重要特征参数，给出预报结论，包括自初始时刻至再入时刻的时间预报。
2.如权利要求1所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述步骤一中，航天器相对地心准惯性坐标系的总能量计算方法包括：
S10、设任意时刻轨道高度为h，则LEO圆轨道至地心距离为H=h+r E ，其中r E 为地球平均半径；于是航天器在轨运行时的总能量为E，其计算公式如下：
Figure QLYQS_1
其中，m E 为地球质量，m为航天器质量，G为万有引力常数，其值为6.67×10-11 N.m2/kg2；
S11、所述S10的公式中，航天器动能E 动的计算公式如下：
Figure QLYQS_2
其中，V为航天器飞行速度，其和航天器运行高度的关系由引力和离心力的平衡关系进行确定；
S12、所述S10的公式中，航天器势能E 势的计算公式如下：
Figure QLYQS_3
其中，航天器的势能零点取在无穷远处；
S13、把航天器初始参数即初始时刻的飞行轨道高度代入S10的公式中获得初始时刻航天器的总能量。
3.如权利要求1所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述步骤二中，大气阻尼对总能量的耗散的计算方法包括：
S20、设时间间隔dt内，航天器运动轨迹近似为ds=Vdt，其中V为航天器飞行速度；W A 为气动阻力对航天器做功，即大气阻尼对总能量的耗散，其计算公式如下：
Figure QLYQS_4
其中，气动阻力做的功主要取决于做功距离ds；实际处理时，时间步长内不考虑气动阻力变化，即时间步长内不考虑相关的阻力系数C D 、大气密度ρ、航天器飞行速度V及迎风截面积A的变化；
S21、所述S20的公式中，气动阻力做功的微分形式或表示如下：
Figure QLYQS_5
其中，通过换算关系把气动阻力做功的表达式表示为航天器距地心距离H的函数，G为万有引力常数，m E 为地球质量；
S22、所述S20的公式中，F A 为气动阻力，表达式为：
Figure QLYQS_6
其中，阻力系数C D 的参考面积取为单位面积。
4.如权利要求1所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述步骤三中，根据耗散后的总能量计算新的圆轨道运行参数的计算方法包括：
S30、设航天器距地心距离H在时间步长dt前后分别是H 1和H 2，则根据耗散后总能量计算新的轨道运行参数，其计算公式如下：
Figure QLYQS_7
其中，ρ为大气密度，C D 为航天器飞行过程中的阻力系数，A为迎风截面积，m为航天器质量，m E 为地球质量，G为万有引力常数；
S31、所述S30的公式中，dH的计算依据即为航天器总能变化dE等于大气阻尼对其耗散，即等于大气阻力对航天器做功W A 的负值，具体关系如下：
dE=-W A
或者
Figure QLYQS_8
由此得出dH的表达式，C D 为阻力系数。
5.如权利要求4所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述步骤四中，再入条件设为某个临界高度h e ，判断是否满足再入条件的方法包括：
S40、根据S30公式得到的航天器距地心距离得到航天器运行高度h 2=H 2-r E，当满足下列条件则说明满足再入条件：
h 2≤h e
于是进入步骤五；
S41、根据S30公式得到的航天器距地心距离得到航天器运行高度h 2=H 2-r E，当满足下列条件则说明尚未满足再入条件，需要再返回步骤二推进计算：
h 2＞h e
于是返回步骤二，在给定推进时间步长后继续计算大气阻尼对总能量的耗散，并获得新的轨道参数。
6.如权利要求5所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述步骤五中，得到再入点参数及统计计算轨道衰降过程的重要特征参数，它们的计算方法包括：
S50、根据步骤四设定的再入高度h e ，即某个临界高度，再入速度V e 的计算公式如下：
Figure QLYQS_9
该公式同时适用于任意高度时的飞行速度计算，即
Figure QLYQS_10
；
S51、根据步骤四设定的再入高度h e ，再入时刻当地弹道倾角d e 的计算公式如下：
d e =arctan（dh e /V e dt e ）
其中，dt e 、dh e 、V e 分别是再入时刻时的时间步长、高度增量及再入速度；该公式同时适用于任意高度时的当地弹道倾角的计算，即d=arctan（dh/Vdt）；
S52、根据步骤二至步骤四推进计算过程对时间步长的累加，得到从初始状态到再入时的时间预报，时间预报t total 的计算公式如下：
Figure QLYQS_11
根据需要，采用不同时间单位，包括“天”、“周”或“年”对时间预报值进行换算。
7.如权利要求5所述的基于能量分析的航天器LEO圆轨道衰降过程参数预报方法，其特征在于，所述临界高度h e 的高度值在90km至120km之间。
CN202310660927.7A 2023-06-06 2023-06-06 基于能量分析的航天器leo圆轨道衰降过程参数预报方法 Active CN116384599B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310660927.7A CN116384599B (zh) 2023-06-06 2023-06-06 基于能量分析的航天器leo圆轨道衰降过程参数预报方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310660927.7A CN116384599B (zh) 2023-06-06 2023-06-06 基于能量分析的航天器leo圆轨道衰降过程参数预报方法

Publications (2)

Publication Number Publication Date
CN116384599A true CN116384599A (zh) 2023-07-04
CN116384599B CN116384599B (zh) 2023-08-22

Family

ID=86979103

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310660927.7A Active CN116384599B (zh) 2023-06-06 2023-06-06 基于能量分析的航天器leo圆轨道衰降过程参数预报方法

Country Status (1)

Country Link
CN (1) CN116384599B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101888339B1 (ko) * 2018-06-01 2018-09-20 한국 천문 연구원 통제 불가능한 인공우주물체 재진입 예측 방법
CN109558660A (zh) * 2018-11-21 2019-04-02 中国航天空气动力技术研究院 一种航天器碎片陨落落区预报方法
CN109992927A (zh) * 2019-04-27 2019-07-09 中国人民解放军32035部队 稀疏数据情况下小椭圆目标的再入预报方法
CN111241634A (zh) * 2019-11-19 2020-06-05 中国空气动力研究与发展中心超高速空气动力研究所 一种航天器再入陨落的分析预报方法
CN111353121A (zh) * 2020-03-31 2020-06-30 中国空气动力研究与发展中心超高速空气动力研究所 一种用于航天器解体碎片不确定性参数的分布方法
CN113343369A (zh) * 2021-08-06 2021-09-03 中国空气动力研究与发展中心设备设计与测试技术研究所 一种航天器气动融合轨道摄动分析方法
US11312512B1 (en) * 2019-03-04 2022-04-26 United States Of America As Represented By The Secretary Of The Air Force Early warning reentry system comprising high efficiency module for determining spacecraft reentry time
CN114580224A (zh) * 2022-05-09 2022-06-03 中国空气动力研究与发展中心设备设计与测试技术研究所 一种分布式气动融合轨道耦合姿态摄动分析方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101888339B1 (ko) * 2018-06-01 2018-09-20 한국 천문 연구원 통제 불가능한 인공우주물체 재진입 예측 방법
CN109558660A (zh) * 2018-11-21 2019-04-02 中国航天空气动力技术研究院 一种航天器碎片陨落落区预报方法
US11312512B1 (en) * 2019-03-04 2022-04-26 United States Of America As Represented By The Secretary Of The Air Force Early warning reentry system comprising high efficiency module for determining spacecraft reentry time
CN109992927A (zh) * 2019-04-27 2019-07-09 中国人民解放军32035部队 稀疏数据情况下小椭圆目标的再入预报方法
CN111241634A (zh) * 2019-11-19 2020-06-05 中国空气动力研究与发展中心超高速空气动力研究所 一种航天器再入陨落的分析预报方法
CN111353121A (zh) * 2020-03-31 2020-06-30 中国空气动力研究与发展中心超高速空气动力研究所 一种用于航天器解体碎片不确定性参数的分布方法
CN113343369A (zh) * 2021-08-06 2021-09-03 中国空气动力研究与发展中心设备设计与测试技术研究所 一种航天器气动融合轨道摄动分析方法
CN114580224A (zh) * 2022-05-09 2022-06-03 中国空气动力研究与发展中心设备设计与测试技术研究所 一种分布式气动融合轨道耦合姿态摄动分析方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
GAO XINGLONG: "Prediction of Orbit Decay for Large-Scale Spacecraft considering Rarefied Aerodynamic Perturbation Effects", 《INTERNATIONAL JOURNAL OF《 AEROSPACE ENGINEERING》, vol. 2022, pages 1 - 14 *
党雷宁等: "基于等效迎角的气动融合轨道直接积分计算无控航天器轨道衰降研究", 《载人航天》, vol. 26, no. 4, pages 452 - 458 *
刘劲宏等: "基于雷诺数的大气阻力模型在飞行器再入预报中的应用", 《空间科学学报》, vol. 42, no. 2, pages 227 - 283 *

Also Published As

Publication number Publication date
CN116384599B (zh) 2023-08-22

Similar Documents

Publication Publication Date Title
CN106697333B (zh) 一种航天器轨道控制策略的鲁棒性分析方法
CN110595485B (zh) 基于两行根数的低轨卫星长期轨道预报方法
CN109255096B (zh) 一种基于微分代数的地球同步卫星轨道不确定演化方法
CN109323698B (zh) 一种空间目标陨落多模型跟踪引导方法
CN109657256B (zh) 一种高精度弹道式再入标称返回轨道仿真方法
CN114936471B (zh) 一种基于并行计算的航天器碰撞预警分层快速筛选方法
CN109992927A (zh) 稀疏数据情况下小椭圆目标的再入预报方法
CN111428339B (zh) 基于空间密度模型的空间物体长期碰撞风险分析方法
CN110489879B (zh) 一种适用于空间环境扰动情况下的空间目标陨落预报方法
CN109752005A (zh) 一种基于精确轨道模型的航天器初轨确定方法
CN116384600B (zh) 基于能量分析的航天器leo椭圆轨道衰降过程参数预报方法
CN110059285B (zh) 考虑j2项影响的导弹自由段弹道偏差解析预报方法
CN116384599B (zh) 基于能量分析的航天器leo圆轨道衰降过程参数预报方法
CN109117543B (zh) 载人航天器对近地小行星探测并返回的轨道设计方法
CN109765141A (zh) 一种基于swarm-c卫星提取大气密度的方法
CN111125874A (zh) 一种动平台高精度测轨预报方法
Grenestedt et al. Dynamic soaring in hurricanes
DeRose Trim attitude, lift and drag of the Apollo command module with offset center-of-gravity positions at Mach numbers to 29
Patera Hazard analysis for uncontrolled space vehicle reentry
CN111382514A (zh) 一种基于监督学习的火星探测飞行轨道精确计算方法及系统
Fattig Comparison of precision orbit derived density estimates for CHAMP and GRACE satellites
Yoon et al. Problems of controlling the motion of small satellite using inflatable thin-film shells for removal space objects from orbit
CN109211225A (zh) 获取大椭圆轨道空间物体剩余轨道寿命的方法、系统及设备
Patera Risk to commercial aircraft from reentering space debris
Boone et al. Artificial debris collision risk following a catastrophic spacecraft mishap in lunar orbit

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant