CN106709145B - 大规模空间碎片分布状态数值演化的并行计算方法 - Google Patents

Abstract

本发明提供一种大规模空间碎片分布状态数值演化的并行计算方法，包括：将碎片分布的空间，划分成N个高度层；将空间碎片按其地心距大小依次配到每个高度层内，确保每个高度层内碎片数量相等；将高度层i所包含的空间碎片分配给唯一对应的计算进程i；计算进程i对高度层i内的每个空间碎片的运动状态按设定时间步长进行积分推演。针对演化过程碎片的相互碰撞事件，通过将每个高度层进一步划分为N_a个子空间体积元，在子空间体积元及其邻近体积元内进行碎片碰撞概率计算。避免对任意两个碎片进行碰撞判断，可有效降低了碰撞概率判断的计算量。优点为：能够在较短时间内，对碎片环境进行长达数百年的，并行高效的数值演化分析。

Description

大规模空间碎片分布状态数值演化的并行计算方法

技术领域

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

背景技术

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

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

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

发明内容

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

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

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

S1，确定空间碎片的初始运动状态，并将空间碎片按照其地心距从小到大排序；

S2，确定并行计算的进程数量，设共有N个计算进程，依次记为：计算进程1、计算进程2…计算进程N；

S3，将地球表面以上180～40000km范围的空间，划分成N个高度层，即：高度层的数量与计算进程的数量相等；将N个高度层按照距离地球表面由近及远的顺序，依次记为：高度层1、高度层2…高度层N；将步骤S1中排好顺序的空间碎片按设定策略分配到每个高度层内，从而确定每个高度层所包含的空间碎片的空间碎片ID以及空间碎片运动状态；将高度层i所包含的空间碎片分配给唯一对应的计算进程i；其中，i∈(1、2…N)；

S4，计算进程i对高度层i内的每个空间碎片的运动状态按设定时间步长进行积分推演，得到下一时刻空间碎片运动状态数据；

S5，然后，计算进程i根据下一时刻空间碎片运动状态数据，判断下一时刻空间碎片是否发生碰撞解体或坠入大气层；如果未发生，则更新空间碎片运动状态数据，返回S4，再次按照给定的时间步长进行积分推演，如此不断循环S4-S5；直到达到演化时间，结束流程；如果发生，则执行S6；

S6，如果下一时刻空间碎片坠入大气层，则计算进程i将其向量容器中已坠入大气层的空间碎片ID删除；然后，返回S4；

如果下一时刻空间碎片发生碰撞解体，则计算进程i将其向量容器中已解体的空间碎片ID删除；然后，对于已解体的空间碎片，计算得到解体后空间碎片的地心距；判断解体后空间碎片的地心距是否属于本计算进程i对应的高度区间，如果属于，则将解体后空间碎片的空间碎片ID以及运动状态数据增加到自身处理范围向量容器中；然后，返回S4；如果不属于，则计算进程i根据解体后空间碎片的地心距，得到对应的高度区间；根据对应的高度区间，得到对应的计算进程j；其中，j≠i，j∈(1、2…N)；

S7，计算进程i与计算进程j建立通信，计算进程i向计算进程j发送增加新处理对象的通知消息，所述通知消息中携带有解体后空间碎片的ID以及运动状态数据；

S8，所述计算进程j在接收到该通知消息后，将解体后空间碎片的空间碎片ID以及运动状态数据增加到其向量容器中。

优选的，S1中，空间碎片从小到大排序后，顺序存放在计算对象向量容器或链表容器中。

优选的，S2具体为：

根据计算平台计算能力，确定可以使用的计算核心数，并根据可以使用的计算核心数，确定并行计算的进程数量。

优选的，S3中，按照高度划分成N个高度层时，每个高度层所包含的空间碎片的数量相等，从而使各个计算进程的计算任务相等。

优选的，S5中，通过以下方法判断下一时刻空间碎片是否坠入大气层：

如果下一时刻空间碎片距地面高度小于100km时，则该空间碎片将会快速坠入大气层而销毁。

优选的，S5中，通过以下方法判断下一时刻空间碎片是否发生碰撞解体：

S5.1，对于每个高度区间，按照经度和纬度间隔，划分为若干个子空间体积元；每个子空间体积元具有唯一的身份标识信息；

S5.2，统计得到每个子空间体积元内包含的空间碎片的空间碎片ID以及空间碎片运动状态；

S5.3，判断每个子空间体积元内的空间碎片是否发生碰撞解体，如果发生，则对应的计算进程将其向量容器中已解体的空间碎片ID删除；并获得与解体后空间碎片的高度层对应的计算进程，然后，将解体后空间碎片的空间碎片ID以及运动状态数据增加到与其高度层对应的计算进程的向量容器中；

其中，通过以下方法判断每个子空间体积元内的空间碎片是否发生碰撞解体：

每个子空间体积元内的空间碎片碰撞，区分为子空间体积元内的空间碎片之间的碰撞，以及邻近子空间体积元的空间碎片之间的相互碰撞两类；

只有当子空间体积元内空间碎片数量大于1时，才判断该子空间体积元内的空间碎片之间是否发生碰撞；

当子空间体积元内存在空间碎片时，需要和邻近26个子空间体积元内的空间碎片进行碰撞概率分析，确定是否发生邻近子空间体积元的空间碎片之间的相互碰撞。

优选的，通过以下公式确定空间碎片碰撞的概率：

p＝4πr²·Δr·s_u(r,t)·s_v(r,t)·σ_u,v·V_i(r)

其中：p为空间碎片碰撞的概率；r为碎片地心距；Δr为高度区间的跨度；t为时间；s_u(r,t)、s_v(r,t)分别为空间碎片u和空间碎片v所属子空间体积元的空间密度；σ_u,v是空间碎片u和空间碎片v之间的平均相互碰撞截面积；V_i(r)是第i个高度区间内空间碎片的平均相对碰撞速度。

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

可以实现对大规模空间碎片的分布状态进行并行高效的数值演化计算，从而可以有效利用超级计算平台的计算资源，能够在短时间内对碎片分布状态进行长达数百年的演化分析，可以为航天器在轨运行空间碎片碰撞规避、空间碎片减缓措施制定以及确保空间资源可持续利用的重要的保障。

附图说明

图1为本发明提供的大规模空间碎片分布状态数值演化的并行计算方法的结构原理图；

图2为本发明提供的高度区间划分示意图；其中，以地球质心为中心，同心球面之间的空间为划分的高区区间。

图3为本发明提供的子空间体积元划分示意图；其中，对于每个高度区间，按照经度和纬度间隔，划分的子空间体积元，每个子空间体积元均由体积元所在的碎片地心距r_i、经度α_k和纬度δ_j唯一确定。

具体实施方式

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

由于空间碎片规模巨大且动态变化，空间碎片在轨运行时受到复杂的摄动作用，致使空间碎片分布数值长期演化面临较大挑战。为了在可承受的时间段内对空间碎片进行上百年的演化分析，利用超级计算机强大的计算能力，如在超级计算机平台上进行演化计算。而为了综合利用超级计算机的计算处理能力，本发明设计高效并行的计算方法，该方法能够处理大规模空间碎片对象，支持规模态变的计算对象，支持跨计算节点、众多线程并发处理。因此，本发明提供的大规模空间碎片分布状态数值演化的并行计算方法，为一种高效的并行计算方法，可实现空间碎片大规模演化。

结合图1-图3，本发明提供一种大规模空间碎片分布状态数值演化的并行计算方法，包括以下步骤：

S1，确定空间碎片的初始运动状态，空间碎片的初始运动状态可以从空间碎片的测量数据中分析得到，如可以利用两行轨道根数(TLE)数据；并将空间碎片按照其地心距从小到大排序；

另外，空间碎片从小到大排序后，顺序存放在计算对象向量容器或链表容器中。

S2，确定并行计算的进程数量，设共有N个计算进程，依次记为：计算进程1、计算进程2…计算进程N；

S2具体为：

根据计算平台计算能力，确定可以使用的计算核心数，并根据可以使用的计算核心数，确定并行计算的进程数量。当初始空间碎片规模不超过2万个，建议计算进程数量不超过100个。

S3，将地球表面以上180～40000km范围的空间，划分成N个高度层，即：高度层的数量与计算进程的数量相等；将N个高度层按照距离地球表面由近及远的顺序，依次记为：高度层1、高度层2…高度层N；如，对于100个计算进程，则需划分100个高度区间。

将步骤S1中排好顺序的空间碎片按设定策略分配到每个高度层内，从而确定每个高度层所包含的空间碎片的空间碎片ID以及空间碎片运动状态；

将高度层i所包含的空间碎片分配给唯一对应的计算进程i；其中，i∈(1、2…N)；

其中，可采用以下方法将空间碎片按设定策略分配到每个高度层：按照高度划分成N个高度层时，每个高度层所包含的空间碎片的数量相等，从而使各个计算进程的计算任务相等。需要强调的是，高度层所包含的碎片数量是相等的，由于碎片在空间中的分布不均匀，因此，各高度层的高度范围不要求相等。

需要注意的是，实际将碎片分配给每个计算进程时，并不一定能够实现严格等量分配。如，对于初始15169个碎片对象，采用100个计算进程，平均每个计算核分配152个碎片对象，但某个计算进程只能分配151个。可以将每个计算进程内的碎片对象也存储在向量容器中，这样可以方便的将坠入大气层的对象删除，或将碰撞产生的新对象追加到计算对象向量容器中。

S4，计算进程i对高度层i内的每个空间碎片的运动状态按设定时间步长进行积分推演，得到下一时刻空间碎片运动状态数据；

每个计算进程只更新本进程内的碎片状态，各计算进程之间的碎片状态不共享。

S5，然后，计算进程i根据下一时刻空间碎片运动状态数据，判断下一时刻空间碎片是否发生碰撞解体或坠入大气层；如果未发生，则更新空间碎片运动状态数据，返回S4，再次按照给定的时间步长进行积分推演，如此不断循环S4-S5；直到达到演化时间，结束流程；如果发生，则执行S6；

其中，通过以下方法判断下一时刻空间碎片是否坠入大气层：

如果下一时刻空间碎片距地面高度小于100km时，则该空间碎片将会快速坠入大气层而销毁。因此，删除的碎片对象状态将不会在下一次状态更新中得到更新，即不再参与数值演化计算。

通过以下方法判断下一时刻空间碎片是否发生碰撞解体：

S5.1，对于每个高度区间，按照经度和纬度间隔，划分为若干个子空间体积元；每个子空间体积元具有唯一的身份标识信息；

子空间体积元是依据赤经和赤纬进行划分的，子空间体积元划分的精密程度取决于计算精度要求；统计每个子空间体积元内碎片对象时，包括从其它计算进程通信传递过来的碎片对象。

具体的，每个高度层内子空间体积元的划分，可以根据等经度和纬度间隔进行划分。如对于赤经α∈[0^°,360^°)，赤纬δ∈[-90^°,90^°]，按照经度间隔Δα＝10^°，纬度间隔Δδ＝2^°，可以将每个高度层进一步划分为(360/10)×(180/2)＝648个子空间体积元。符号D_i,j,k表示碎片D所处的空间网格：i表示该碎片当前时刻处于第i个高度区间，即被第i个计算进程处理；j表示碎片处于第j个经度区间；k表示碎片处于第k个纬度区间。则每个碎片都可以通过上述三个下标i、j和k来描述其所属子空间体积元。对每个高度区间内空间碎片遍历一遍，即可确定该高度区间内碎片所在的空间体积元。

S5.2，统计得到每个子空间体积元内包含的空间碎片的空间碎片ID以及空间碎片运动状态；

S5.3，判断每个子空间体积元内的空间碎片是否发生碰撞解体，如果发生，则对应的计算进程将其向量容器中已解体的空间碎片ID删除；并获得与解体后空间碎片的高度层对应的计算进程，然后，将解体后空间碎片的空间碎片ID以及运动状态数据增加到与其高度层对应的计算进程的向量容器中；

也就是说，当解体后空间碎片的地心距超出本计算进程对应的高度层范围时，计算进程与其他计算进程进行通信，传递空间碎片对象数据，且这种通信是基于一定策略和约定实施的。

在具体实现上，进程间通信策略的约定如下：线程之间高效、可靠的数据通信，是确保并行计算结果正确的基础。在设计线程间通信策略的过程中，首先假定每个计算线程仅与邻近的两个计算线程间发生数据通信。这一假定的物理背景是：对碎片的运动状态进行一步积分更新后，碎片的地心距变化量不会超过邻近高度区间的长度。对于在圆轨道上运行的碎片，这一假设通常是成立的。而80％以上的空间碎片均运行在圆轨道上，因此这一假设在大部分情况下是成立的。对于一些特殊的大椭圆轨道，如Molniya轨道，碎片的地心距会在短时间内发生较大的变化，通过设置合适的基本步长，也能确保上述假设成立。在上述假定的基础上，可以设计一个确定的通信策略，即每个高度区间均进行两次数据发送和两次数据接收通信，且这种通信仅与邻近的高度区间进行。

例如，计算进程之间按照一定通信策略进行通信，交换碎片计算对象，具体实现过程为：计算进程根据碎片地心距大小，判断碎片是否仍在当前进程对应的高区区间内。记高度区间邻近H_i的两个高度区间分别为左高度区间H_i-1和右高度区间H_i+1，H_i-1的高度上限与H_i的高度下限重合，同样H_i+1的高度下限与H_i的高度上限重合。若状态更新后碎片的地心距小于高度区间H_i的高度下限值，则认为该碎片进入邻近的左高度区间H_i-1，则将该碎片对象存储在临时碎片对象容器V_i-1中，并将其从高区区间H_i的对象容器中删除；若状态更新后碎片的地心距大于高度区间H_i的高度上限值，则认为该碎片进入邻近的右高度区间H_i+1，则将该碎片对象存储在临时碎片对象容器V_i+1中，同样将其从高区区间H_i的对象容器中删除。进程间通信时，与高度区间H_i对应的计算进程分别将临时碎片对象容器V_i-1和V_i+1中的碎片对象，传递给高度区间H_i-1和H_i+1对应的计算进程。

其中，通过以下方法判断每个子空间体积元内的空间碎片是否发生碰撞解体：

每个子空间体积元内的空间碎片碰撞，区分为子空间体积元内的空间碎片之间的碰撞，以及邻近子空间体积元的空间碎片之间的相互碰撞两类；

只有当子空间体积元内空间碎片数量大于1时，才判断该子空间体积元内的空间碎片之间是否发生碰撞；

当子空间体积元内存在空间碎片时，需要和邻近26个子空间体积元内的空间碎片进行碰撞概率分析，确定是否发生邻近子空间体积元的空间碎片之间的相互碰撞。

在判断碎片是否发生碰撞时，由于空间碎片受到各种摄动力影响，空间碎片的位置矢量存在较大的不确定性，因此在微观尺度上，即在一个空间体积元内，可以近似认为空间碎片的位置是随机的，碎片之间的碰撞概率可以近似用气体运动学理论来确定。即可以利用下述公式确定：

p＝4πr²·Δr·s_u(r,t)·s_v(r,t)·σ_u,v·V_i(r)

其中：p为空间碎片碰撞的概率；r为碎片地心距；Δr为高度区间的跨度；t为时间；s_u(r,t)、s_v(r,t)分别为空间碎片u和空间碎片v所属子空间体积元的空间密度；σ_u,v是空间碎片u和空间碎片v之间的平均相互碰撞截面积；V_i(r)是第i个高度区间内空间碎片的平均相对碰撞速度。

S6，如果下一时刻空间碎片坠入大气层，则计算进程i将其向量容器中已坠入大气层的空间碎片ID删除；然后，返回S4；

如果下一时刻空间碎片发生碰撞解体，则计算进程i将其向量容器中已解体的空间碎片ID删除；然后，对于已解体的空间碎片，计算得到解体后空间碎片的地心距；判断解体后空间碎片的地心距是否属于本计算进程i对应的高度区间，如果属于，则将解体后空间碎片的空间碎片ID以及运动状态数据增加到自身处理范围向量容器中；然后，返回S4；如果不属于，则计算进程i根据解体后空间碎片的地心距，得到对应的高度区间；根据对应的高度区间，得到对应的计算进程j；其中，j≠i，j∈(1、2…N)；

S7，计算进程i与计算进程j建立通信，计算进程i向计算进程j发送增加新处理对象的通知消息，所述通知消息中携带有解体后空间碎片的ID以及运动状态数据；

S8，所述计算进程j在接收到该通知消息后，将解体后空间碎片的空间碎片ID以及运动状态数据增加到其向量容器中。

另外，对于每个计算进程，可包括空间碎片状态推演子模块和空间碎片碰撞/爆炸解体产生新碎片子模块，这两个子模块可具体采用以下方法进行积分推演，需要强调的是，下面描述的仅为一种优选实施方式，本发明对此并不限制：

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

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

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

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

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

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

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

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

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

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

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

其中：

G为引力常数，M_e是地球质量，r为地心距；

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

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

a_e表示地球平均半径；

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

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

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

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

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

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

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

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

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

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

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

针对初值问题：

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

龙格-库塔积分公式为：

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

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

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

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

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

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

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

最后，根据A/m_d、A确定碎片质量m_d。

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

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

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

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

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

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

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

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

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

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

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

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

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

其中：

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

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

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

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

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

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

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

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

p(γ,θ)＝ε(θ)p₁(γ)+(1-ε(θ))p₂(γ) (19)

其中：p₁(γ)和p₂(γ)均是正态分布概率密度函数：

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

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

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

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

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

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

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

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

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

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

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

。

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

综上所述，本发明提供的大规模空间碎片分布状态数值演化的并行计算方法，具有以下优点：

(1)本发明在进行大规模空间碎片分布状态数值演化计算过程中，能够综合利用多个计算进程，可以有效利用超级计算平台的计算资源，实现在段时间内对空间碎片分布状态进行计算分析，具有演化过程效率高的优点；

(2)本发明解决了大规模空间碎片数值演化计算过程中碎片之间相互碰撞判断的效率问题。即：通过划分子空间体积元，当子空间体积元内碎片数量大于1时，才判断子空间体积元内碎片是否碰撞；当子空间体积元内碎片数量大于0时，分析与邻近空间体积元内碎片的碰撞情况。这种综合利用空间网格判断碎片之间是否发生碰撞的方法，可以避免碎片之间两两对比，有效降低了碰撞概率判断的计算量，从而更有利于大规模碎片的长期演化计算。

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

Claims (6)

1.一种大规模空间碎片分布状态数值演化的并行计算方法，其特征在于，包括以下步骤：

S1，确定空间碎片的初始运动状态，并将空间碎片按照其地心距从小到大排序；

S1中，空间碎片从小到大排序后，顺序存放在计算对象向量容器或链表容器中；

S2，确定并行计算的进程数量，设共有N个计算进程，依次记为：计算进程1、计算进程2…计算进程N；

S3，将地球表面以上180～40000km范围的空间，划分成N个高度层，即：高度层的数量与计算进程的数量相等；将N个高度层按照距离地球表面由近及远的顺序，依次记为：高度层1、高度层2…高度层N；将步骤S1中排好顺序的空间碎片按设定策略分配到每个高度层内，从而确定每个高度层所包含的空间碎片的空间碎片ID以及空间碎片运动状态；将高度层i所包含的空间碎片分配给唯一对应的计算进程i；其中，i∈(1、2…N)；

S4，计算进程i对高度层i内的每个空间碎片的运动状态按设定时间步长进行积分推演，得到下一时刻空间碎片运动状态数据；

S5，然后，计算进程i根据下一时刻空间碎片运动状态数据，判断下一时刻空间碎片是否发生碰撞解体或坠入大气层；如果未发生，则更新空间碎片运动状态数据，返回S4，再次按照给定的时间步长进行积分推演，如此不断循环S4-S5；直到达到演化时间，结束流程；如果发生，则执行S6；

S6，如果下一时刻空间碎片坠入大气层，则计算进程i将其向量容器中已坠入大气层的空间碎片ID删除；然后，返回S4；

如果下一时刻空间碎片发生碰撞解体，则计算进程i将其向量容器中已解体的空间碎片ID删除；然后，对于已解体的空间碎片，计算得到解体后空间碎片的地心距；判断解体后空间碎片的地心距是否属于本计算进程i对应的高度区间，如果属于，则将解体后空间碎片的空间碎片ID以及运动状态数据增加到自身处理范围向量容器中；然后，返回S4；如果不属于，则计算进程i根据解体后空间碎片的地心距，得到对应的高度区间；根据对应的高度区间，得到对应的计算进程j；其中，j≠i，j∈(1、2…N)；

S7，计算进程i与计算进程j建立通信，计算进程i向计算进程j发送增加新处理对象的通知消息，所述通知消息中携带有解体后空间碎片的ID以及运动状态数据；

S8，所述计算进程j在接收到该通知消息后，将解体后空间碎片的空间碎片ID以及运动状态数据增加到其向量容器中。

2.根据权利要求1所述的大规模空间碎片分布状态数值演化的并行计算方法，其特征在于，S2具体为：

根据计算平台计算能力，确定可以使用的计算核心数，并根据可以使用的计算核心数，确定并行计算的进程数量。

3.根据权利要求1所述的大规模空间碎片分布状态数值演化的并行计算方法，其特征在于，S3中，按照高度划分成N个高度层时，每个高度层所包含的空间碎片的数量相等，从而使各个计算进程的计算任务相等。

4.根据权利要求1所述的大规模空间碎片分布状态数值演化的并行计算方法，其特征在于，S5中，通过以下方法判断下一时刻空间碎片是否坠入大气层：

如果下一时刻空间碎片距地面高度小于100km时，则该空间碎片将会快速坠入大气层而销毁。

5.根据权利要求1所述的大规模空间碎片分布状态数值演化的并行计算方法，其特征在于，S5中，通过以下方法判断下一时刻空间碎片是否发生碰撞解体：

S5.1，对于每个高度区间，按照经度和纬度间隔，划分为若干个子空间体积元；每个子空间体积元具有唯一的身份标识信息；

S5.2，统计得到每个子空间体积元内包含的空间碎片的空间碎片ID以及空间碎片运动状态；

S5.3，判断每个子空间体积元内的空间碎片是否发生碰撞解体，如果发生，则对应的计算进程将其向量容器中已解体的空间碎片ID删除；并获得与解体后空间碎片的高度层对应的计算进程，然后，将解体后空间碎片的空间碎片ID以及运动状态数据增加到与其高度层对应的计算进程的向量容器中；

其中，通过以下方法判断每个子空间体积元内的空间碎片是否发生碰撞解体：

每个子空间体积元内的空间碎片碰撞，区分为子空间体积元内的空间碎片之间的碰撞，以及邻近子空间体积元的空间碎片之间的相互碰撞两类；

只有当子空间体积元内空间碎片数量大于1时，才判断该子空间体积元内的空间碎片之间是否发生碰撞；

当子空间体积元内存在空间碎片时，需要和邻近26个子空间体积元内的空间碎片进行碰撞概率分析，确定是否发生邻近子空间体积元的空间碎片之间的相互碰撞。

6.根据权利要求5所述的大规模空间碎片分布状态数值演化的并行计算方法，其特征在于，通过以下公式确定空间碎片碰撞的概率：

p＝4πr²·△r·s_u(r,t)·s_v(r,t)·σ_u,v·V_i(r)

其中：p为空间碎片碰撞的概率；r为碎片地心距；△r为高度区间的跨度；t为时间；s_u(r,t)、s_v(r,t)分别为空间碎片u和空间碎片v所属子空间体积元的空间密度；σ_u,v是空间碎片u和空间碎片v之间的平均相互碰撞截面积；V_i(r)是第i个高度区间内空间碎片的平均相对碰撞速度。