本公开涉及空间物体碰撞概率计算领域，尤其涉及一种基于cube模型的碰撞概率计算方法。
背景技术：
空间碎片的日益增多引发了国际上关于空间环境长期稳定性的研究。而计算机技术的发展使得对n体系统的长期演化仿真成为可能。在空间碎片系统演化中，需要考虑碎片的增加机制和减少机制，与之相应的，需要一种能够利用不断更新的碎片轨道数据进行碰撞概率计算的算法。美国国家宇航局研制了一种称为“cube”模型的碰撞概率算法。cube模型的核心思想是，在空间碎片系统长期演化中，每隔一定时间dt对系统采样一次，在采样时刻，将整个近地空间划分为长度h的小立方体，利用更新的轨道数据计算出该时刻每个空间物体所在的位置，然后找出处于同一个立方体的空间物体，计算其两两之间的碰撞概率，若一个立方体中只有一个物体，则不考虑其与其他立方体中物体的碰撞概率。cube模型的优点在于通过对空间划分立方体可以对n体碎片系统进行快速配对，大大减少了计算量；其缺点则是遗漏了对处于相邻立方体中的空间物体之间碰撞概率的计算，这种碰撞概率甚至会大于同一立方体中两个物体之间的碰撞概率。
技术实现要素：
(一)要解决的技术问题
本公开旨在弥补cube模型中没有考虑相邻立方体中空间物体之间碰撞概率的缺陷，提出一种基于cube模型的碰撞概率计算方法。
(二)技术方案
本公开提供了一种基于cube模型的碰撞概率计算方法，包括：将近地空间划分立方体网格，并对立方体进行编号；计算空间物体的位置及其所处的立方体编号；对空间物体进行配对，将处于相邻和相同立方体内的空间物体作为配对目标；以及计算配对目标之间的碰撞概率。
在本公开的一些实施例中，所述将近地空间划分立方体网格，并对立方体进行编号包括：以地心为坐标系原点，建立笛卡尔坐标系，将近地空间划分为立方体网格，并对每个立方体以边长为单位进行编号。
在本公开的一些实施例中，所述编号表示为(x，y，z)，x、y、z分别表示立方体在x轴方向、y轴方向和z轴方向的编号。
在本公开的一些实施例中，所述计算空间物体的位置及其所处的立方体编号包括：根据每个空间物体更新后的轨道根数计算其所在的位置，查找所述位置所处的立方体，并记录所处立方体的编号。
在本公开的一些实施例中，所述对空间物体进行配对，将处于相邻和相同立方体内的空间物体作为配对目标包括：分别在三个坐标轴方向，按序号顺序将空间物体所处立方体的编号进行排序，得到与目标空间物体所处立方体的编号相邻和相同的编号；与所述编号对应的立方体内的空间物体序号构成三个坐标轴方向的空间物体序号集；计算三个坐标轴方向的空间物体序号集的交集，得到与目标空间物体所处立方体相邻和相同的立方体的空间物体序号集；遍历空间物体序号集中的空间物体与目标空间物体的欧氏距离，保留欧式距离在小于阈值的空间物体作为配对目标。
在本公开的一些实施例中，利用二分法按序号顺序将空间物体所处立方体的编号进行排序。
在本公开的一些实施例中，所述阈值取立方体边长乘以
在本公开的一些实施例中，所述计算配对目标之间的碰撞概率包括：对处于相邻和相同立方体内的空间物体，空间物体i和j的平均碰撞数为：
c＝sisjvimpacdudt
碰撞概率为：
pij＝1-exp(-c)
其中，du为体元；dt为时间采样间隔；si、sj分别表示空间物体i和j在体元du内的分布密度；vimp是空间物体i和j的相对碰撞速度；ac为空间物体i和j的碰撞截面。
在本公开的一些实施例中，所述体元du是以空间物体i和j之间欧式距离的阈值为半径的球体体积。
在本公开的一些实施例中，其中，
空间物体在体元du内均匀分布，
si＝sj＝1/du
其中，dc为空间物体i和j之间欧式距离的阈值；h为立方体边长。
(三)有益效果
从上述技术方案可以看出，本公开基于cube模型的碰撞概率计算方法具有以下有益效果：
(1)本公开通过增加对相邻cube空间物体的碰撞可能性的考虑，弥补了原有cube模型可能遗漏碰撞情况的缺陷，使得改进后的cube模型更符合物理真实，分析结果更加合理可信。
(2)本公开适用于空间碎片系统在长期演化过程中对碰撞概率的评估，有利于提高空间碎片长期演化模型的合理性与可信度。而对空间碎片长期演化模型的研究能够对未来我国进行主动碎片清除提供策略分析和模型支撑。
附图说明
图1是本公开实施例的基于cube模型的碰撞概率计算方法的流程图。
具体实施方式
本公开通过对空间物体快速配对，计算空间物体之间的碰撞概率。该方法适用于空间碎片系统在长期演化过程中对碰撞概率的评估。在本公开中，考虑空间目标之间的碰撞，可将其视为空气分子间的碰撞，依据气体动力学来评估空间目标之间的碰撞概率。空气分子间的碰撞概率服从泊松分布。
为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。
本公开实施例提供了一种基于cube模型的碰撞概率计算方法，参见图1，包括：
步骤s1：将近地空间划分立方体网格，并对立方体进行编号。
该步骤具体包括：在采样时刻，以地心为坐标系原点，建立笛卡尔坐标系，将地球周围的近地空间划分为边长为h的立方体形成的网格，并对每个立方体按边长h的倍数进行编号，即以边长h为单位进行编号。
立方体的编号表示为(x，y，z)，x、y、z分别表示立方体在x轴方向、y轴方向、z轴方向的编号；x轴方向编号表示该立方体处于x轴的第x个单位的位置，同理，y轴方向编号表示该立方体处于y轴的第y个单位的位置，y轴方向编号表示该立方体处于z轴的第z个单位的位置。
步骤s2：计算空间物体的位置及其所处立方体的编号。
该步骤具体包括：在采样时刻，根据每个空间物体更新后的轨道根数计算其所在的位置，查找该位置所处的立方体，并记录所处立方体的编号。在该步骤中可先对空间物体进行排序，对于序号为i的空间物体，其位置坐标表示为(xi，yi，zi)，该位置所处立方体的编号表示为(xi，yi，zi)。
步骤s3：对空间物体进行配对，将处于相邻和相同立方体内的空间物体作为配对目标。
该步骤是在采样时刻，通过对所有空间物体所处立方体的编号进行比对，查找出处于相邻或者相同立方体内的空间目标。
具体来说，
首先，分别在三个坐标轴方向，按序号顺序将空间物体所处立方体的编号进行二分法排序，得到与目标空间物体所处立方体的编号相邻和相同的编号。
与上述编号对应的立方体内的空间物体序号构成三个坐标轴方向的空间物体序号集。
在一个示例中，可以先按空间物体的序号顺序，将各个空间物体所处立方体的x轴方向编号进行二分法排序，得到空间物体所处立方体的x轴方向编号集，记录x轴方向编号集中与空间物体i所处立方体的x轴方向编号相邻和相同的x轴方向编号，得到与这些x轴方向编号对应的立方体的集合，该集合中的立方体内的空间物体序号构成x轴方向空间物体序号集mx。
此处的二分法排序具体包括：在对第k个空间物体所处立方体的x轴方向编号进行排序时，之前的k-1个空间物体所处立方体的x轴方向编号已经按照x轴排好序并形成数组nk-1，此时首先利用二分法查找第k个空间物体所处立方体的x轴方向编号在数组nk-1位置，然后将第k个空间物体所处立方体的x轴方向编号插入nk-1的对应位置得到数组nk。
同样，按序号顺序将空间物体所处立方体的y轴方向编号进行二分法排序，得到空间物体所处立方体的y轴方向编号集，记录y轴方向编号集中与空间物体i所处立方体的y轴方向编号相邻和相同的y轴方向编号，得到与这些y轴方向编号对应的立方体的集合，该集合中的立方体内的空间物体序号构成y轴方向空间物体序号集my。按序号顺序将空间物体所处立方体的z轴方向编号进行二分法排序，得到空间物体所处立方体的z轴方向编号集，记录z轴方向编号集中与空间物体i所处立方体的z轴方向编号相邻和相同的z轴方向编号，得到与这些z轴方向编号对应的立方体的集合，该集合中的立方体内的空间物体序号构成z轴方向空间物体序号集mz。
在其他示例中，也可以先进行y轴方向或z轴方向的排序。
然后，计算三个坐标轴方向的空间物体序号集的交集mi＝mx∩my∩mz，得到与空间物体i所处立方体相邻和相同立方体的空间物体序号集mi。
最后，遍历空间物体序号集mi中的空间物体与空间物体i的欧氏距离，保留欧式距离在小于阈值dc的空间物体作为配对目标。
在一个示例中，阈值dc取
步骤s4：计算配对目标之间的碰撞概率。
对处于相邻和相同立方体内的空间物体，通过以下方法计算其两两之间的碰撞概率：
根据气体动力学，在体元du内，在dt时间内，空间物体i和j的平均碰撞数为：
c＝sisjvimpacdudt(1)
碰撞次数服从泊松统计，故碰撞概率为：
pij＝1-exp(-c)(2)
其中，si、sj分别表示两个空间物体在体元du内的分布密度，vimp是两个空间物体的相对碰撞速度，ac为两个空间物体的碰撞截面。碰撞截面
本方法中考虑的配对目标是欧氏距离在dc之内的空间物体，故du不再是立方体的体积，而是以dc为半径的球体体积，即
同时空间物体在体元du内均匀分布，故
si＝sj＝1/du(4)
时间采样间隔dt设为5天，可计算得到空间物体i和j之间的碰撞概率pij。
综上所述，本公开基于cube模型而开发的碰撞概率方法，通过增加对相邻cube空间物体的碰撞可能性的考虑，弥补了原有cube模型可能遗漏碰撞情况的缺陷，使得改进后的cube模型更符合物理真实，分析结果更加合理可信。本公开适用于空间碎片系统在长期演化过程中对碰撞概率的评估，有利于提高空间碎片长期演化模型的合理性与可信度。而对空间碎片长期演化模型的研究能够对未来我国进行主动碎片清除提供策略分析和模型支撑。
至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本公开有了清楚的认识。
需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换，例如：
(1)实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向，并非用来限制本公开的保护范围；
(2)上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。
以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。