[发明专利]大规模空间碎片分布演化的数值计算方法及系统

摘要

本发明提供一种大规模空间碎片分布演化的数值计算方法及系统，系统包括：空间碎片状态推演子系统、空间碎片碰撞/爆炸解体产生新碎片子系统、空间碎片运动状态三维显示子系统、空间碎片运动状态输出子系统和空间碎片演化结果处理子系统。本发明提供的大规模空间碎片分布演化的数值计算方法及系统具有以下优点：本发明提供的大规模空间碎片分布演化的数值计算方法及系统，能够对大规模空间碎片分布状态进行长期演化计算，能够动态三维演示碎片分布状态，提供了碎片分布状态数据的输出和处理分析功能，可以为航天器在轨运行空间碎片碰撞规避、空间碎片减缓措施制定以及确保空间资源可持续利用的重要的保障。

主权项

1.一种大规模空间碎片分布演化的数值计算方法，其特征在于，包括以下步骤：S1，空间碎片状态推演子系统获取空间碎片的初始运动状态数据以及空间碎片受到的空间摄动力，并将所述初始运动状态数据以及所述空间摄动力输入到数值积分模型中；所述数值积分模型按照给定的时间步长进行积分推演，得到空间碎片未来运动状态数据；然后，采用所述空间碎片未来运动状态数据更新空间碎片当前运动状态数据；空间碎片状态推演子系统具体实现过程如下：步骤S1.1，计算空间碎片受到的大气阻力摄动力；在轨空间碎片受到的大气阻力随着大气状态的不同而变化，阻力加速度计算式为其中：是大气阻力摄动力；CD是空间碎片的阻力系数；A是空间碎片的迎风截面积；md是空间碎片的质量；ρ是空间碎片所在位置的大气密度；是空间碎片与当地大气的相对速度；假设大气相对地球静止，即大气相对于地球的速度其中为地球自转角速度，为地心矢径；若空间碎片的速度为则相对速度阻力加速度的计算式可进一步改写为步骤S1.2，计算空间碎片受到的地球非球形引力摄动力；在地固坐标系中，地球引力位函数作为拉普拉斯方程的解，其非球形部分U为： U = Σ n = 2 N Σ m = 0 n ( C ‾ n m U ‾ n m + S ‾ n m V ‾ n m ) - - - ( 3 ) ]]>其中：和为勒让德多项式相关项，定义为：其中：G为引力常数，Me是地球质量，r为地心距；U表示地球引力场位函数；λ和分别表示单位质点在地固坐标系中的地心经和纬度；ae表示地球平均半径；为规格化的缔合勒让德多项式；和为规格化的地球引力位系数；n和m为多项式的阶和次，N为取的最高阶数；步骤S1.3，计算空间碎片受到的太阳光压摄动力；太阳光压摄动力加速度可以表示为其中：为太阳光压摄动力加速度；pSP＝4.5605×10-6N/m2为地球附近太阳光压常数；CR为太阳光压系数；ASR是与太阳垂直的目标横截面积；AU＝1.49597870×1011m是IAU1976天文常数系统给出的天文单位；为太阳位置矢量；步骤S1.4，计算空间碎片受到的第三体引力摄动力；第三体引力摄动力加速度为其中：为第三体引力摄动力加速度；S和L分别代表太阳和月球；Mj为太阳或月球质量；代表太阳或月球的地心矢径；太阳和月球的位置可以通过平根公式计算，也可以采用JPL星历确定；步骤S1.5，利用4阶阿达姆斯预估矫正方法对空间碎片的运动状态进行积分更新；该积分方法利用4阶龙格-库塔方法进行参数初始化，充分利用历史计算结果，每步积分仅计算一次摄动力，减少积分计算量；采用4阶积分方法，比高阶积分具有跟高的计算稳定性；4阶龙格-库塔参数初始化方法为：针对初值问题： y · ( t ) = f ( t , y ) y ( t ) = y n - - - ( 7 ) ]]>其中：t是时刻；y是状态量，如空间碎片的位置和速度矢量的坐标；是状态变量y的导数；f(t,y)是约束函数，包括空间碎片受到大气阻力加速度、地球非球形引力摄动力加速度以及太阳光压摄动力加速度；tn表示离散化后的不同时刻；yn是第n个时刻状态变量的取值；龙格-库塔积分公式为： y n + 1 = y n + h Σ i = 0 k C i f i - - - ( 8 ) ]]> f 0 = f ( t n , y n ) f 1 = f ( t n + a 1 h , y n + a 1 b 10 hf 0 ) . . . f k = f ( t n + a k h , y n + a k h Σ i = 0 k b k i f i ) - - - ( 9 ) ]]>其中：h＝tn+1-tn，为积分步长；k为积分式所取的阶数，当取值为4的时候即为4阶龙格-库塔初始化方法；Ci,ai,bij均为已知的常数项；设已知4个时刻的函数值分别为fi-3,fi-2,fi-1,fi，则i+1时刻y的近似估计值为 y i + 1 = y i + h 24 ( - 9 f - 3 + 37 f i - 2 - 59 f i - 1 + 55 f i ) - - - ( 10 ) ]]>S2，判断是否达到演化时间，如果未达到，将状态更新后的空间碎片运动状态和属性数据传输给空间碎片碰撞/爆炸解体产生新碎片子系统；S3，所述空间碎片碰撞/爆炸解体产生新碎片子系统将所述状态更新后的空间碎片运动状态和属性数据输入到碎片解体模型中；所述碎片解体模型进行解体计算，得到空间碎片解体后产生新空间解体碎片的属性数据；并将所述新空间解体碎片的属性数据反馈传输到所述空间碎片状态推演子系统；S4，所述空间碎片状态推演子系统将所述新空间解体碎片的属性数据以及已有的空间碎片的运动状态数据输入到数值积分模型中，再次按照给定的时间步长进行积分推演，得到空间碎片未来运动状态数据，如此不断循环S2-S4；直到达到演化时间，结束流程。