大规模空间碎片分布演化的数值计算方法及系统与流程

本发明属于航天动力学技术和计算机技术领域，具体涉及一种大规模空间碎片分布演化的数值计算方法及系统。

背景技术：

随着人类空间活动的日益增加，人为发射进入空间的物体不断增多，使得地球外部空间的物体数量迅速增加，在轨航天器受到的碰撞威胁也愈加严重。

据nasa空间监视网(ssn)的统计，地球同步轨道高度以下的大于10cm尺寸的空间物体数量达到15000多个，小于10cm尺寸的空间物体数量更多，粗略估计至少为400000个，致使地球外部空间达到了前所未有的拥挤状态。这些空间物体中，只有很少的一部分是在轨正常工作的航天器，大部分是由于发射任务、在轨航天器爆炸分解、在轨航天器碰撞分解产生的空间碎片。因此，这些在轨道动力学的约束下沿轨道绕地球自由运动的空间碎片，会挤占有限的空间资源，使得诸如近地轨道和静止轨道等区域可用的资源减少。同时，在大气阻力、地球非球形摄动以及太阳、月球三体引力等摄动力的作用下，自由无控的空间碎片轨道会随时间不断变化，当空间碎片轨道与正常工作航天器的轨道存在交点时，就会存在与航天器发生碰撞的可能，对正常工作航天器产生威胁。另外，航天器发生碰撞后，会进一步使空间碎片的数量增加，从而使得产生新的碰撞的可能性增加，进而产生更多碎片，对正常运行航天器产生更大威胁，这种链式的雪崩式效应，最终会导致航天器在轨运行的成本迅速增加，甚至会出现近地或同步轨道范围内航天器均无法生存的情况。

由此可见，大量空间碎片对航天器造成的碰撞威胁，将是人类空间活动面临的主要安全问题，如不加以管理，将导致航天器运行成本大大增加，甚至会使空间资源无法利用。有效的对当前空间碎片的运动分布状态进行预测分析，是航天器在轨运行碰撞规避、空间碎片减缓措施制定以及确保空间资源可持续利用的重要的技术保障。现有技术中，尚未出现有效的对当前空间碎片的运动分布状态进行预测分析的相关内容。

技术实现要素：

针对现有技术存在的缺陷，本发明提供一种大规模空间碎片分布演化的数值计算方法及系统，可有效解决上述问题。

本发明采用的技术方案如下：

本发明提供一种大规模空间碎片分布演化的数值计算方法，包括以下步骤：

s1，空间碎片状态推演子系统获取空间碎片的初始运动状态数据以及空间碎片受到的空间摄动力，并将所述初始运动状态数据以及所述空间摄动力输入到数值积分模型中；所述数值积分模型按照给定的时间步长进行积分推演，得到空间碎片未来运动状态数据；然后，采用所述空间碎片未来运动状态数据更新空间碎片当前运动状态数据；

s2，判断是否达到演化时间，如果未达到，将状态更新后的空间碎片运动状态和属性数据传输给空间碎片碰撞/爆炸解体产生新碎片子系统；

s3，所述空间碎片碰撞/爆炸解体产生新碎片子系统将所述状态更新后的空间碎片运动状态和属性数据输入到碎片解体模型中；所述碎片解体模型进行解体计算，得到空间碎片解体后产生新空间解体碎片的属性数据；并将所述新空间解体碎片的属性数据反馈传输到所述空间碎片状态推演子系统；

s4，所述空间碎片状态推演子系统将所述新空间解体碎片的属性数据以及已有的空间碎片的运动状态数据输入到数值积分模型中，再次按照给定的时间步长进行积分推演，得到空间碎片未来运动状态数据，如此不断循环s2-s4；直到达到演化时间，结束流程。

优选的，s1中，所述空间摄动力包括大气阻力摄动力、地球非球形引力摄动力、太阳光压摄动力以及第三体引力摄动力；所述数值积分模型为4阶阿达姆斯预估矫正积分模型。

优选的，s3中，所述空间碎片碰撞/爆炸解体产生新碎片子系统具体用于：

首先，根据能量定律确定空间碎片解体后产生新空间解体碎片的数量n和新空间解体碎片特征尺寸l的分布；

然后，以新空间解体碎片特征尺寸l作为独立变量，利用概率分布模型确定新空间解体碎片的面质比a/md、新空间解体碎片的迎风截面积a、新空间解体碎片相对于解体前空间碎片的速度增量δv；

最后，根据a/md、a确定碎片质量md。

优选的，s3中，所述空间碎片碰撞/爆炸解体产生新碎片子系统具体用于：

s3.1，确定解体碎片数量：

其中，解体碎片区分为爆炸解体碎片和碰撞解体碎片；

1)对于爆炸解体碎片，根据能量守恒定律，解体碎片特征尺寸l大于被研究解体碎片最小特征尺寸lc的碎片数量nf满足函数关系：

其中：为解体碎片的最小特征尺寸；

cs是修正系数，与爆炸类型有关，对于历史爆炸事件，修正系数cs的表达式：

其中：

其中：h为爆炸解体碎片的高度；

2)对于碰撞解体碎片，碰撞解体碎片数量满足函数关系：

式中：为两碰撞目标的等效质量：

其中：

其中：ep为等效碰撞能量，为灾难性碰撞临界能量；mp为相互碰撞中较小质量；mt为相互碰撞中较大质量；vr是碰撞时两碎片的相对速度，即碰撞速度；

s3.2，确定解体碎片尺寸分布律：

爆炸解体碎片尺寸分布律f：

碰撞解体碎片尺寸分布律f：

s3.3，确定解体碎片面质比参数：

本步骤中，爆炸解体碎片和碰撞解体碎片的碎片面质比参数满足同样的分布函数；

(1)对于尺寸大于11cm解体碎片，定义解体碎片面质比的对数值γ＝log(a/md)，解体碎片特征长度的对数值θ＝log(l)，解体碎片的面质比分布由如下二元正态分布决定：

p(γ,θ)＝ε(θ)p1(γ)+(1-ε(θ))p2(γ)(19)

其中：p1(γ)和p2(γ)均是正态分布概率密度函数：

其中：权重参数ε(θ)∈[0,1]，是θ的饱和函数；正态分布概率密度函数的均值μk和正态分布概率密度函数的方差σk均是θ的函数；

对于航天器解体碎片，其面质比分布函数中的参数由下式决定：

对于火箭上面级解体碎片，其面质比分布函数的参数由下式决定

(2)对于尺寸小于8cm的航天器解体碎片和尺寸小于1.7cm的火箭上面级解体碎片，碎片的面质比满足正态分布律p1(γ)，即权重函数ε＝1，

p(γ,θ)＝p1(γ)(23)

此时，分布函数参数由下式统一决定：

(3)对于尺寸介于8cm～11cm之间的航天器解体碎片和尺寸介于1.7cm～11cm之间的火箭上面级碎片，碎片的面质比通过随机采样的方法确定，即：首先生成服从均匀分布且取值区间在[0.0,1.0]内的随机变量ζ，然后将ζ与利用(25)计算得到的进行对比，若则利用公式(19)计算面质比参数，反之则利用公式(23)计算面质比参数；

步骤3.4，确定解体碎片的速度增量：

解体碎片相对于解体前目标的速度增量满足正态分布规律，定义v＝lg(δv)，γ＝lg(a/md)，则v的分布函数p(ν)满足

其中：v为解体碎片的速度；δv为解体碎片的速度增量；

解体碎片的均值μν和解体碎片的方差σν通过下式计算：

σv＝0.4。

优选的，还包括：

s5，空间碎片运动状态三维显示子系统将所述空间碎片状态推演子系统推演到的空间碎片未来运动状态数据进行三维显示，即：动态的显示所有空间碎片在空间中的位置分布状态；具体包括：

步骤s5.1，接收空间碎片状态推演子系统推演到的空间碎片未来运动状态数据；

步骤s5.2，调用三维显示接口函数，将空间碎片的位置和形状在三维模型中进行渲染；

步骤s5.3，驱动三维模型，对空间碎片运动状态进行动态更新。

优选的，还包括：

s6，在固定时刻或者设定的时刻，空间碎片运动状态输出子系统将所述空间碎片状态推演子系统推演到的空间碎片未来运动状态输出；其中，所述空间碎片未来运动状态包括：空间碎片的位置、速度矢量和面值比、质量、来源属性参数；此外，空间碎片运动状态输出子系统还输出不同高度层内空间碎片状态数据。

优选的，还包括：

s7，空间碎片演化结果处理子系统根据空间碎片状态推演子系统推演到的空间碎片未来运动状态，计算空间碎片状态量，进一步计算到空间碎片密度；根据空间碎片密度得到空间碎片沿轨道高度方向和纬度方向上的分布状态。

优选的，所述空间碎片密度通过以下步骤计算得到：

步骤s7.1，确定空间碎片在高度区间[r,r′]内的平均空间密度分布其中，r和r′分别为高度区间的最小地心距和最大地心距；包括以下步骤：

步骤s7.1.1，当目标轨道为椭圆轨道，近地点地心距和远地点地心距分别为ra和rp，当r′≤ra且r≥rp时，通过下式计算：

其中：a为目标轨道半长轴；

步骤s7.1.2，当目标轨道为椭圆轨道，且r′＞ra或r＜rp时，则令r′＝ra或r＝rp，带入(28)进行计算，得到

步骤s7.1.3，当目标轨道为圆轨或近圆轨道，即rp→ra∈[r,r′]时，则令r′＝ra或r＝rp，公式(28)变为公式(29)，通过公式(29)计算

步骤s7.1.4，其他情况，

步骤s7.2，确定空间碎片在纬度区间[β,β′]内的平均空间密度分布函数其中，β和β′分别是纬度区间的最小值和最大值；

具体的，若β′＞0，则直接执行步骤s7.2.1；

若β′≤0，则令β′＝-β，β＝-β′，均转换为关于原点对称的正值区间后，再执行步骤s7.2.1；

步骤s7.2.1，碎片轨道倾角记为i，当β′≤i且β≥-i时

步骤s7.2.2，当β≤-i且β′≤i时

步骤s7.2.3，当β′≥i且-i≤β≤i时

步骤s7.2.4，当β′≥i且β≤-i时

步骤s7.2.5，当β＞i时，

本发明还提供一种大规模空间碎片分布演化的数值计算系统，包括：

空间碎片状态推演子系统，用于获取空间碎片的初始运动状态数据以及空间碎片受到的空间摄动力，并将所述初始运动状态数据以及所述空间摄动力输入到数值积分模型中；所述数值积分模型按照给定的时间步长进行积分推演，得到空间碎片未来运动状态数据；然后，采用所述空间碎片未来运动状态数据更新空间碎片当前运动状态数据，得到空间碎片新运动状态数据；

空间碎片碰撞/爆炸解体产生新碎片子系统，用于将所述空间碎片新运动状态数据输入到碎片解体模型中；所述碎片解体模型进行解体计算，得到空间碎片解体后产生新空间解体碎片的属性数据；并将所述新空间解体碎片的属性数据反馈传输到所述空间碎片状态推演子系统；

空间碎片运动状态三维显示子系统，用于将所述空间碎片状态推演子系统推演到的空间碎片未来运动状态数据进行三维显示，即：动态的显示所有空间碎片在空间中的位置分布状态；

空间碎片运动状态输出子系统，用于将所述空间碎片状态推演子系统推演到的空间碎片未来运动状态输出；其中，所述空间碎片未来运动状态包括：空间碎片的位置、速度矢量和面值比；此外，空间碎片运动状态输出子系统还输出不同高度层内空间碎片状态数据；

空间碎片演化结果处理子系统，用于根据空间碎片状态推演子系统推演到的空间碎片未来运动状态，计算空间碎片状态量，进一步计算到空间碎片密度；根据空间碎片密度得到空间碎片沿轨道高度方向和纬度方向上的分布状态。

本发明提供的大规模空间碎片分布演化的数值计算方法及系统具有以下优点：

本发明提供的大规模空间碎片分布演化的数值计算方法及系统，能够对大规模空间碎片分布状态进行长期演化计算，能够动态三维演示碎片分布状态，提供了碎片分布状态数据的输出和处理分析功能，可以为航天器在轨运行空间碎片碰撞规避、空间碎片减缓措施制定以及确保空间资源可持续利用的重要的保障。

附图说明

图1为本发明提供的大规模空间碎片分布演化的数值计算方法的流程示意图；

图2为本发明提供的空间碎片状态推演子系统的结构示意图；

图3是大规模空间碎片分布演化数值计算系统运行时，空间碎片初始分布状态的三维显示图；

其中：在图3中，积分步数为0步；当前空间碎片数量为17000个；已发生碰撞次数为0次；

图4是大规模空间碎片分布演化数值计算系统运行过程中，产生新的解体碎片的分布状态示意图。

具体实施方式

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

本发明提供一种大规模空间碎片分布演化的数值计算方法，应用于由5个子系统构成的大规模空间碎片分布演化的数值计算系统中，大规模空间碎片分布演化的数值计算系统包括：空间碎片状态推演子系统、空间碎片碰撞/爆炸解体产生新碎片子系统、空间碎片运动状态三维显示子系统、空间碎片运动状态输出子系统和空间碎片演化结果处理子系统。

大规模空间碎片分布演化的数值计算方法包括以下步骤：

s1，空间碎片状态推演子系统获取空间碎片的初始运动状态数据以及空间碎片受到的空间摄动力，并将所述初始运动状态数据以及所述空间摄动力输入到数值积分模型中；所述数值积分模型按照给定的时间步长进行积分推演，得到空间碎片未来运动状态数据；然后，采用所述空间碎片未来运动状态数据更新空间碎片当前运动状态数据。

本步骤中，所述空间摄动力包括大气阻力摄动力、地球非球形引力摄动力、太阳光压摄动力以及第三体引力摄动力；所述数值积分模型为4阶阿达姆斯(adams-bashforth)预估矫正积分模型。

空间碎片状态推演子系统的结构如图2所示，具体实现过程如下：

步骤s1.1，计算空间碎片受到的大气阻力摄动力。

大气阻力作为唯一清除在轨空间碎片的因素，其作用随轨道高度变化较为明显。计算表明，250km的典型空间碎片，在大气阻力的作用下，2个月内会再入到大气层内；而600km高度的空间碎片，再入到大气层的时间达到15年，高于850km的空间碎片，再入到大气层的时间将需要数个世纪之久，而地球同步轨道处的空间碎片在通常情况下永远不会再入到大气层中，除非空间环境发生了大规模变化。

在轨空间碎片受到的大气阻力随着大气状态的不同而变化，阻力加速度计算式为

其中：是大气阻力摄动力；cd是空间碎片的阻力系数；a是空间碎片的迎风截面积；md是空间碎片的质量；ρ是空间碎片所在位置的大气密度；是空间碎片与当地大气的相对速度；

假设大气相对地球静止，即大气相对于地球的速度其中为地球自转角速度，为地心矢径；若空间碎片的速度为则相对速度阻力加速度的计算式可进一步改写为

步骤s1.2，计算空间碎片受到的地球非球形引力摄动力。

在地固坐标系中，地球引力位函数作为拉普拉斯方程的解，其非球形部分u为：

其中：和为勒让德多项式相关项，定义为：

其中：

g为引力常数，me是地球质量，r为地心距；

u表示地球引力场位函数；

λ和分别表示单位质点在地固坐标系中的地心经和纬度；

ae表示地球平均半径；

为规格化的缔合勒让德多项式；

和为规格化的地球引力位系数；

n和m为多项式的阶和次，n为取的最高阶数。

步骤s2.3，计算空间碎片受到的太阳光压摄动力。

太阳光压摄动力加速度可以表示为

其中：为太阳光压摄动力加速度；psp＝4.5605×10-6n/m2为地球附近太阳光压常数；cr为太阳光压系数；asr是与太阳垂直的目标横截面积；au＝1.49597870×1011m是iau1976天文常数系统给出的天文单位；为太阳位置矢量。

步骤s2.4，计算空间碎片受到的第三体引力摄动力。

第三体引力摄动力加速度为

其中：为第三体引力摄动力加速度；s和l分别代表太阳和月球；mj为太阳或月球质量；代表太阳或月球的地心矢径。太阳和月球的位置可以通过平根公式计算，也可以采用jpl星历确定。

步骤s2.5，利用4阶阿达姆斯预估矫正方法对空间碎片的运动状态进行积分更新。该积分方法利用4阶龙格-库塔方法进行参数初始化，充分利用历史计算结果，每步积分仅计算一次摄动力，减少积分计算量；采用4阶积分方法，比高阶积分具有跟高的计算稳定性。

4阶龙格-库塔参数初始化方法为：

针对初值问题：

其中：t是时刻；y是状态量，如空间碎片的位置和速度矢量的坐标；是状态变量y的导数；f(t,y)是约束函数，如空间碎片受到大气阻力加速度、地球非球形引力摄动力加速度以及太阳光压摄动力加速度等；tn表示离散化后的不同时刻；yn是第n个时刻状态变量的取值。

龙格-库塔积分公式为：

其中：h＝tn+1-tn，为积分步长；k为积分式所取的阶数，当取值为4的时候即为4阶龙格-库塔初始化方法；ci,ai,bij均为已知的常数项。

设已知4个时刻的函数值分别为fi-3,fi-2,fi-1,fi，则i+1时刻y的近似估计值为

s2，判断是否达到演化时间，如果未达到，将状态更新后的空间碎片运动状态和属性数据传输给空间碎片碰撞/爆炸解体产生新碎片子系统；

s3，s3，所述空间碎片碰撞/爆炸解体产生新碎片子系统将所述状态更新后的空间碎片运动状态和属性数据输入到碎片解体模型中；所述碎片解体模型进行解体计算，得到空间碎片解体后产生新空间解体碎片的属性数据；并将所述新空间解体碎片的属性数据反馈传输到所述空间碎片状态推演子系统；

其中，所述空间碎片碰撞/爆炸解体产生新碎片子系统，具体用于：

首先，根据能量定律确定空间碎片解体后产生新空间解体碎片的数量n和新空间解体碎片特征尺寸l的分布；其中，碎片的特征尺寸l＝(lx+ly+lz)/3，即为碎片3个轴向尺寸的平均。

然后，以新空间解体碎片特征尺寸l作为独立变量，利用概率分布模型确定新空间解体碎片的面质比a/md、新空间解体碎片的迎风截面积a、新空间解体碎片相对于解体前空间碎片的速度增量δv；

最后，根据a/md、a确定碎片质量md。

更具体的，所述空间碎片碰撞/爆炸解体产生新碎片子系统具体用于：

s3.1，确定解体碎片数量：

其中，解体碎片区分为爆炸解体碎片和碰撞解体碎片；

1)对于爆炸解体碎片，根据能量守恒定律，解体碎片特征尺寸l大于被研究解体碎片最小特征尺寸lc的碎片数量nf满足函数关系：

其中：为解体碎片的最小特征尺寸；

cs是修正系数，与爆炸类型有关，对于历史爆炸事件，修正系数cs的表达式：

其中：修正系数cs与目标探测直径门限dssn直接相关，dssn与目标高度h的关系式如下：

其中：h为爆炸解体碎片的高度；

2)对于碰撞解体碎片，解体碎片的数量和分布与相互碰撞的能量有关。碰撞分为灾难性碰撞和非灾难性碰撞，两种碰撞条件下解体碎片的数量和分布规律满足不同的表达式。

一次碰撞是否为灾难性碰撞可通过下式判断：

其中：ep为等效碰撞能量，为灾难性碰撞临界能量。mp为相互碰撞中较小质量；mt为相互碰撞中较大质量；vr是碰撞时两碎片的相对速度，即碰撞速度；当上式成立时，碰撞为灾难性碰撞，反之则为非灾难性碰撞。

因此，给出下面的碰撞解体碎片数量满足函数关系：

式中：为两碰撞目标的等效质量：

其中：

s3.2，确定解体碎片尺寸分布律：

爆炸解体碎片尺寸分布律f：

碰撞解体碎片尺寸分布律f：

计算时，首先确定要研究解体碎片的最小尺寸lc，然后根据关系式(11)和(15)计算爆炸或解体产生大于lc尺寸碎片的总数目，再根据分布律(17)或(18)确定这些碎片具有的尺寸。利用分布律确定解体碎片尺寸时，需要对随机变量抽样，可采用反函数的方法进行随机变量的抽样。

s3.3，确定解体碎片面质比参数：

得到解体碎片的特征尺寸以后，碎片的面质比参数满足确定的分布函数，利用随机变量采样法可以得到每一个碎片具有的面质比。

本步骤中，爆炸解体碎片和碰撞解体碎片的碎片面质比参数满足同样的分布函数；

(1)对于尺寸大于11cm解体碎片，定义解体碎片面质比的对数值γ＝log(a/md)，解体碎片特征长度的对数值θ＝log(l)，解体碎片的面质比分布由如下二元正态分布决定：

p(γ,θ)＝ε(θ)p1(γ)+(1-ε(θ))p2(γ)(19)

其中：p1(γ)和p2(γ)均是正态分布概率密度函数：

其中：权重参数ε(θ)∈[0,1]，是θ的饱和函数；正态分布概率密度函数的均值μi和正态分布概率密度函数的方差σi均是θ的函数；

对于航天器解体碎片，其面质比分布函数中的参数由下式决定：

对于火箭上面级解体碎片，其面质比分布函数的参数由下式决定

(2)对于尺寸小于8cm的航天器解体碎片和尺寸小于1.7cm的火箭上面级解体碎片，碎片的面质比满足正态分布律p1(γ)，即权重函数ε＝1，

p(γ,θ)＝p1(γ)(23)

此时，分布函数参数由下式统一决定：

(3)对于尺寸介于8cm～11cm之间的航天器解体碎片和尺寸介于1.7cm～11cm之间的火箭上面级碎片，碎片的面质比通过随机采样的方法确定，即：首先生成服从均匀分布且取值区间在[0.0,1.0]内的随机变量ζ，然后将ζ与利用(25)计算得到的进行对比，若则利用公式(19)计算面质比参数，反之则利用公式(23)计算面质比参数；

步骤3.4，确定解体碎片的速度增量：

解体碎片相对于解体前目标的速度增量满足正态分布规律，定义v＝lg(δv)，γ＝lg(a/md)，则v的分布函数p(ν)满足

其中：v为解体碎片的速度；δv为解体碎片的速度增量；

解体碎片的均值μν和解体碎片的方差σν通过下式计算：

σv＝0.4。

s4，所述空间碎片状态推演子系统将所述新空间解体碎片的属性数据以及已有的空间碎片的运动状态数据输入到数值积分模型中，再次按照给定的时间步长进行积分推演，得到空间碎片未来运动状态数据，如此不断循环s2-s4；直到达到演化时间，结束流程。

还包括：

s5，空间碎片运动状态三维显示子系统将所述空间碎片状态推演子系统推演到的空间碎片未来运动状态数据进行三维显示，即：利用directx等三维显示技术，动态的显示所有空间碎片在空间中的位置分布状态；

具体包括：

步骤s5.1，接收空间碎片状态推演子系统推演到的空间碎片未来运动状态数据；

步骤s5.2，调用三维显示接口函数(api)，将空间碎片的位置和形状在三维模型中进行渲染；

步骤s5.3，驱动三维模型，对空间碎片运动状态进行动态更新。

空间碎片的分布状态如图3和图4所示。

还包括：

s6，在固定时刻或者设定的时刻，空间碎片运动状态输出子系统将所述空间碎片状态推演子系统推演到的空间碎片未来运动状态输出；其中，所述空间碎片未来运动状态包括：空间碎片的位置、速度矢量和面值比、质量、来源属性参数；此外，空间碎片运动状态输出子系统还输出不同高度层内空间碎片状态数据。

还包括：

s7，空间碎片演化结果处理子系统根据空间碎片状态推演子系统推演到的空间碎片未来运动状态，计算空间碎片状态量，进一步计算到空间碎片密度；根据空间碎片密度得到空间碎片沿轨道高度方向和纬度方向上的分布状态。

所述空间碎片密度通过以下步骤计算得到：

步骤s7.1，确定空间碎片在高度区间[r,r′]内的平均空间密度分布其中，r和r′分别为高度区间的最小地心距和最大地心距；包括以下步骤：

步骤s7.1.1，当目标轨道为椭圆轨道，近地点地心距和远地点地心距分别为ra和rp，当r′≤ra且r≥rp时，通过下式计算：

其中：a为目标轨道半长轴；

步骤s7.1.2，当目标轨道为椭圆轨道，且r′＞ra或r＜rp时，则令r′＝ra或r＝rp，带入(28)进行计算，得到

步骤s7.1.3，当目标轨道为圆轨或近圆轨道，即rp→ra∈[r,r′]时，则令r′＝ra或r＝rp，公式(28)变为公式(29)，通过公式(29)计算

步骤s7.1.4，其他情况，

步骤s7.2，确定空间碎片在纬度区间[β,β′]内的平均空间密度分布函数其中，β和β′分别是纬度区间的最小值和最大值；

具体的，若β′＞0，则直接执行步骤s7.2.1；

若β′≤0，则令β′＝-β，β＝-β′，均转换为关于原点对称的正值区间后，再执行步骤s7.2.1；

步骤s7.2.1，碎片轨道倾角记为i，当β′≤i且β≥-i时

步骤s7.2.2，当β≤-i且β′≤i时

步骤s7.2.3，当β′≥i且-i≤β≤i时

步骤s7.2.4，当β′≥i且β≤-i时

步骤s7.2.5，当β＞i时，

本发明还提供一种大规模空间碎片分布演化的数值计算系统，包括空间碎片状态推演子系统、空间碎片碰撞/爆炸解体产生新碎片子系统、空间碎片运动状态三维显示子系统、空间碎片运动状态输出子系统和空间碎片演化结果处理子系统。下在对这5个子系统详细介绍：

(1)空间碎片状态推演子系统

空间碎片状态推演子系统，用于获取空间碎片的初始运动状态数据以及空间碎片受到的空间摄动力，并将所述初始运动状态数据以及所述空间摄动力输入到数值积分模型中；所述数值积分模型按照给定的时间步长进行积分推演，得到空间碎片未来运动状态数据；然后，采用所述空间碎片未来运动状态数据更新空间碎片当前运动状态数据，得到空间碎片新运动状态数据。

也就是说，空间碎片状态推演子系统以空间碎片的初始运动状态为输入，综合考虑空间碎片受到的大气阻力、非球形引力、太阳光压、三体引力等摄动力的作用，采用高效稳定的积分方法，积分推演碎片的未来的运动状态。

其中，在计算大气阻力的过程中，要首先确定高层大气密度。太阳辐射作用使大气密度呈现周日变化、季节变化等周期性变化规律；太阳活动会导致大气密度发生剧烈变化，严重影响大气阻力对目标的衰减作用。计算系统选取了harris-priester大气密度模型，该模型能够反映太阳活动对大气密度的影响，模型复杂性低，对计算能力要求不高，适宜进行大规模的长期演化计算。

计算空间碎片受到的非球形摄动时，计算系统对其量级进行了分析，并进行了合理的取舍。地球扁率j2项是主要的摄动项，导致空间碎片半长轴和升交点赤经周期性变化，影响着物体沿经度方向上的分布规律；在扁率j2项和j2,2项影响下，geo轨道带物体在地球赤道短轴上空，即东经75°和西经105°附近，呈现聚集分布规律；j3项和j4项的量级与j2,2项量级相当，均为10^-6，计算空间碎片受到的非球形摄动时也加以考虑。

计算系统中的积分方法是经过优选确定的。对空间碎片的运动进行动力学积分时，要求积分算法精度高且具有有一定稳定性，避免长期计算时结果发散；由于空间碎片数量规模大，需要开发高效的并行计算方法。系统采用4阶阿达姆斯预估矫正方法(adams-bashforth)进行积分，该积分方法利用4阶龙格-库塔方法进行参数初始化，充分利用历史计算结果，每步积分仅计算一次摄动力，减少积分计算量；采用4阶积分方法，比高阶积分具有跟高的计算稳定性；模型基于mpi的并行计算框架，利用多进程和多核处理技术，设计了高效的并行演化系统。

(2)空间碎片碰撞/爆炸解体产生新碎片子系统

空间碎片碰撞/爆炸解体产生新碎片子系统，用于将所述空间碎片新运动状态数据输入到碎片解体模型中；所述碎片解体模型进行解体计算，得到空间碎片解体后产生新空间解体碎片的属性数据；并将所述新空间解体碎片的属性数据反馈传输到所述空间碎片状态推演子系统。

空间碎片碰撞/爆炸解体产生新碎片子系统的功能包括两个方面：首先确定空间碎片之间发生碰撞的概率，然后根据碎片解体模型生成大量解体碎片。然后，根据能量定律确定解体后产生碎片的数量n和碎片特征尺寸l的分布；然后以特征尺寸l作为独立变量，利用概率分布模型确定碎片的面质比a/md、有效横截面积a和速度增量δv；最后根据a/md、a确定碎片质量md。

(3)空间碎片运动状态三维显示子系统

空间碎片运动状态三维显示子系统，用于利用directx等三维显示技术，将所述空间碎片状态推演子系统推演到的空间碎片未来运动状态数据进行三维显示，子系统能够根据物体的推演状态，动态的显示所有物体在空间中的位置分布状态。

(4)空间碎片运动状态输出子系统

空间碎片运动状态输出子系统，用于在固定时刻或者设定的时刻，将所述空间碎片状态推演子系统推演到的空间碎片未来运动状态输出，具体的，输出包块解体碎片在内的当前所有碎片的位置、速度矢量和物体的面值比等属性参数；此外，空间碎片运动状态输出子系统还输出不同高度层内空间碎片状态数据；

(5)空间碎片演化结果处理子系统

空间碎片演化结果处理子系统，用于根据空间碎片状态推演子系统推演到的空间碎片未来运动状态，分析碎片沿轨道高度方向、纬度方向等空间维度方向上的变化情况，分析碎片规模随时间的变化规律。计算空间碎片状态量，进一步计算到空间碎片密度；根据空间碎片密度得到空间碎片沿轨道高度方向和纬度方向上的分布状态。

本发明提供的大规模空间碎片分布演化数值计算方法及系统，能够根据碎片的初始运动状态，综合考虑大气阻力、地球非球形摄动力以及太阳/月球三体引力等复杂摄动力的作用，碎片相互碰撞产生新碎片的作用，对碎片的分布状态进行长期演化计算。系统能够三维动态的显示空间碎片在轨分布状态随时间的变化，能够对碎片分布状态的演化计算结果进行统计分析，得到空间碎片在不同轨道高度、不同纬度等方向上的分布规律，从而为航天器在轨运行碰撞规避、空间碎片减缓措施制定以及确保空间资源可持续利用的重要的技术保障。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。