一种基于并行计算的航天器碰撞预警分层快速筛选方法

1.本发明属于航天器碰撞预警筛选领域，具体涉及基于并行计算的航天器碰撞预警分层快速筛选方法。


背景技术：

2.随着航天事业不断发展，各国的星座计划以及太空实验接连部署，越来越多的卫星被送入太空中，这也预示着太空中的空间碎片数目会急剧上升。据不完全统计，截至目前，在轨大小超过10cm的碎片约为34000个，大小在1cm到10cm之间的碎片约为90万，在1mm到1cm之间的碎片超过1.2亿。而目前由人类送入地球轨道的卫星超过10000颗，因此极有可能出现空间目标之间的发生碰撞情况，造成严重的后果，尤其对于近地轨道目标，具有较高的相对运动速度，产生的破坏性更为强烈。为了空间目标的在轨安全，避免造成不可挽回的后果，有必要研究如何快速筛选出对在轨航天器可能造成威胁的空间目标。
3.针对在轨航天器的碰撞预警问题，传统的筛选方法通常是针对绝不可能接近的目标，即轨道高度相差较大或相遇窗口距离相差较远的目标进行筛选，亦或是仅针对单星之间的高精度筛选。但由于实际任务中的空间目标数目庞大，运动情况也各不相同，仅采用几何筛选和一些简单的二体摄动环境下的筛选不能够满足精度需求，而均使用高精度数值方法进行筛选会消耗大量时间，实用性较小。因此，需要一种能够平衡筛选精度与筛选时间的筛选方法。


技术实现要素：

4.本发明的目的是为了解决现有筛选方法仅采用几何筛选和一些简单的二体摄动环境下的筛选不能够满足精度需求，而均使用高精度数值方法进行筛选会消耗大量时间，实用性较小的问题，而提出一种基于并行计算的航天器碰撞预警分层快速筛选方法。
5.一种基于并行计算的航天器碰撞预警分层快速筛选方法具体过程为：
6.步骤一、以空间目标的tle根数作为初始数据，首先解算得到所有目标的初始的轨道六根数和近远地点，然后利用预筛选模块从几何角度上筛除不可能接近的空间目标，得到保留的空间目标；
7.步骤二、以目标航天器的历元时刻作为初始时刻，将步骤一中得到的保留的空间目标送入gpu中，在gpu中首先利用线性j2模型在预报周期内进行遍历预报，得到所有目标在相应时间节点的空间位置，然后利用基于误差修正的box方法对所有目标进行筛选，保留存在接近风险的目标以及接近时刻，送回cpu进行储存；
8.步骤三、将保留目标的每个接近时刻前后分别扩充δt得到预报时段，利用sgp4模型对所有目标在预报时段内进行预报，得到最小接近距离，保留最小接近距离小于等于最小距离阈值的目标，利用hpop模型对保留的目标进行递推，得到最小接近距离；
9.步骤四、利用步骤三中保留目标的历史tle数据生成初始时刻的协方差矩阵，并在二体模型下线性外推到接近时刻，得到协方差矩阵，结合最小接近距离利用pc方法的解析
公式得到碰撞概率，与概率阈值比较后，对大于阈值的目标给出碰撞警告。
10.本发明的有益效果：
11.本发明提出了一种基于并行计算的近地轨道碰撞预警快速筛选方法。在本发明中，首先利用最小二乘法得到一种基于迹向误差修正的box方法，相较于传统的box固定空间域的阈值设置，可以更好地降低阈值门限，减少虚警数，然后在此基础上搭建包括预筛选、初筛和精筛在内的完整的筛选模型，并通过gpu对初筛模块进行并行加速。采用本发明方法进行空间目标碰撞预警的筛选问题，只需要给出目标航天器和待筛选目标的初始tle根数以及过去10天的历史tle根数，即可快速筛选得到有碰撞风险的目标，并具有一定的筛选精度。
附图说明
12.图1为本发明基于并行计算的近地轨道碰撞预警快速筛选方法流程图；
13.图2为本发明筛选模块示意图；
14.图3为本发明迹向位置误差修正方法示意图；
15.图4为本发明每个线程块内并行结构示意图；
16.图5为本发明初筛过程示意图；
17.图6为本发明精筛过程示意图；
18.图7为本发明完整筛选流程示意图；
19.图8为修正前迹向位置误差分布情况图；
20.图9为本发明修正后迹向位置误差分布情况图。
具体实施方式
21.具体实施方式一：本实施方式一种基于并行计算的航天器碰撞预警分层快速筛选方法具体过程为：
22.本发明设计了一种基于并行计算的近地轨道碰撞预警快速筛选方法，能够对3天内目标航天器有威胁的目标进行快速筛选。该方法由预筛选、初筛和精筛三个筛选模块组成，适用于近地轨道所有的目标，在保证一定的筛选精度情况下有较高的计算效率，最后能够以碰撞概率的解析形式量化碰撞风险作为预警指标。
23.步骤一、以空间目标的tle根数作为初始数据，首先解算得到所有目标的初始的轨道六根数和近远地点，然后利用预筛选模块从几何角度上筛除不可能接近的空间目标，得到保留的空间目标；
24.步骤二、以目标航天器的历元时刻作为初始时刻，将步骤一中得到的保留的空间目标送入gpu中，在gpu中首先利用线性j2模型在预报周期内(根据任务实际需要设定预报时长，通常为3天，不超过5天)进行遍历预报，得到所有目标在相应时间节点的空间位置，然后利用基于误差修正的box方法对所有目标进行筛选，保留存在接近风险的目标以及接近时刻，送回cpu进行储存；
25.步骤三、将保留目标的每个接近时刻前后分别扩充δt得到预报时段，利用sgp4模型对所有目标在预报时段内进行预报，得到最小接近距离，保留最小接近距离小于等于最小距离阈值的目标，利用hpop模型对保留的目标进行递推，得到最小接近距离；
26.步骤四、利用步骤三中保留目标的历史tle数据生成初始时刻的协方差矩阵，并在二体模型下线性外推到接近时刻，得到协方差矩阵，结合最小接近距离利用pc方法的解析公式得到碰撞概率，与概率阈值比较后，对大于阈值的目标给出碰撞警告。
27.具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤一中以空间目标的tle根数作为初始数据，首先解算得到所有目标的初始的轨道六根数和近远地点，然后利用预筛选模块从几何角度上筛除不可能接近的空间目标，得到保留的空间目标；具体过程为：
28.整个筛选过程使用的惯性系均为j2000地心惯性坐标系；
29.空间中的主目标航天器是卫星，备选目标可以是卫星、空间碎片等，只要其在空间监视网的范围内，具有tle根数即可；
30.tle根数是北美空防空令部发布的人造天体的运行状态数据，不能直接用于轨道预报，需要使用sgp4模型进行解算；
31.首先使用sgp4模型对主目标和备选目标的tle根数进行解算，得到主目标和备选目标历元时刻的初始轨道六根数和近远地点；
32.采用预筛选模块从几何角度上筛除不可能接近的空间目标，几何角度上筛除包含近远地点筛选(公式(1)(2))和倾角筛选(公式(3)(4)(5))两个部分；
33.对于给定主目标和备选目标的初始六根数分别为[a1,e1,i1,ω1,ω1,f1]和[a2,e2,i2,ω2,ω2,f2]，可以分别得到两个轨道的对应的近远地点和
[0034][0035]
a1表示主目标对应的瞬时半长轴，e1表示主目标对应的离心率，i1表示主目标对应的目标轨道的倾角，ω1表示主目标对应的升/降交点赤经，ω1表示主目标对应的近地点角距，f1表示主目标对应的真近点角；
[0036]
a2表示备选目标对应的瞬时半长轴，e2表示备选目标对应的离心率，i2表示备选目标对应的目标轨道的倾角，ω2表示备选目标对应的升/降交点赤经，ω2表示备选目标对应的近地点角距，f2表示备选目标对应的真近点角；
[0037]
对轨道高度相差较大的轨道利用下式便以进行直接筛除，对满足下式的目标删除(比较还有和之间与δh的关系)，保留近远地点半径均在筛选阈值内的目标；
[0038]
对满足下式的目标删除
[0039][0040]
式中，δh为高度阈值，表示主目标的远地点半径，表示备选目标的远地点半径，表示主目标的近地点半径，表示备选目标的近地点半径，||表示或；
[0041]
此外，对于轨道面之间倾角较大的轨道，接近时刻通常只有两个轨道面交线与轨道交点处，通过对两个交点之间高度差和距离阈值比较即可筛除不可能接近目标，记两个轨道面的角动量矢量分别为h1和h2，轨道面交线矢径为：
[0042][0043]
交线有两个共线不同向的矢径和分别对应于式(3)中的正负关系。得出矢
径后，继续计算得到矢径分别和每条轨道交点的真近点角f
±
(正负号分别对应一组交点，一共有两组，)，分别对应不同的矢径：
[0044][0045]
式中，是升/降交点的矢径，是在地心惯性坐标系下，由地心指向轨道升(降)交点的矢量，分别对应和(对应对应对应)；|h3|表示交线矢径的模，通常规一化后取为1；表示升/降交点矢径的模，规一化后取为1；ωi表示轨道近地点幅角；
[0046]
从而可以通过式(5)求出两个轨道在同一交点处的轨道半径，作差后与高度阈值比较，将满足公式(6)的目标进行筛除：
[0047][0048]
|r
1-r2|≥δh
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0049]
式中，r1表示两个轨道在同一交点处时轨道1的轨道半径，r2表示两个轨道在同一交点处时轨道2的轨道半径。
[0050]
其它步骤及参数与具体实施方式一相同。
[0051]
具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤二中以目标航天器的历元时刻作为初始时刻，将步骤一中得到的保留的空间目标送入gpu中，在gpu中首先利用线性j2模型在预报周期内(根据任务实际需要设定预报时长，通常为3天，不超过5天)进行遍历预报，得到所有目标在相应时间节点的空间位置，然后利用基于误差修正的box方法对所有目标进行筛选，保留存在接近风险的目标以及接近时刻，送回cpu进行储存；
[0052]
具体过程为：
[0053]
将步骤一中得到的保留的空间目标，送入gpu中(通过cuda平台送入gpu中进行并行加速)，在gpu中对步骤一保留的空间目标利用线性j2模型在预报周期内进行轨道递推，其中半长轴采用平均根数，得到空间目标各预报时刻的轨道六根数后，基于空间目标各预报时刻的轨道六根数转化得到对应的惯性系下空间位置，继而得到待筛选空间目标在主目标的轨道系下空间位置；
[0054]
然后利用基于迹向误差修正的box方法对所有目标进行初步筛选，具体过程为：
[0055]
该方法利用最小二乘法拟合的多项式对迹向误差进行修正。首先引入一个新的自变量对每个空间目标在轨道系下的位置误差的离散点进行线性化：
[0056][0057]
式中，r
p
表示初始近地点半径，x表示用于线性化的新的自变量，e表示离心率，i表示目标点纬度的倾角，ω表示升/降交点赤经；
[0058]
取n个近地空间目标作打靶仿真，得到3天内沿迹向的误差最大值δr，从而可以得到一组数据序列(xj,δrj)(j＝1,2,...,n)，使用拟合函数如下：
[0059][0060]
式中，表示待拟合的函数，bk表示系数，表示单项函数，n表示拟合函数包含单项函数的总项数，xj表示每个离散点的自变量大小，δrj表示每个离散点位置误差的大小；j表示拟合的数据点共有n个，k表示拟合的多项式是从0到n阶；
[0061]
在本方法中采用多项式进行拟合，即取从而可以解得系数b0,b1,...,bn；具体的拟合流程如图3所示；
[0062]
获取主目标的轨道系下待筛选空间目标中主目标和备选目标的迹向位置、径向位置和法向位置；
[0063]
将主目标和备选目标的径向位置作差得到径向相对位置rr；
[0064]
将主目标和备选目标的法向位置作差得到法向相对位置rw；
[0065]
利用拟合得到的多项式对轨道系下待测空间目标中的主目标和备选目标的迹向位置进行修正，得到修正后的迹向相对位置rs；具体过程为：
[0066]
将待测空间目标中的主目标和备选目标的离散点数据xi带入多项式得到待测空间目标中主目标和备选目标的离散点数据
[0067]
将轨道系下的主目标和备选目标的迹向位置分别减去此时对应得到修正后的主目标和备选目标的迹向位置，将修正后的主目标和备选目标的迹向位置作差得到修正后的迹向相对位置rs；
[0068]
将迹向相对位置rs、径向相对位置rr、法向相对位置rw分别与设置的碰撞阈值进行比较，保留小于等于阈值的目标(满足式(9)的目标)：
[0069][0070]
其中δrr,δrs,δrw分别为轨道系下径向、迹向和法向的碰撞阈值，从而可以得到所有存在接近风险的目标(满足式(9)的目标)及对应的时刻j
′
＝0,1,...,m；i
′
＝0,1,...,nj
′
；
[0071]
其中j
′
表示针对所有保留目标中第j
′
个目标，i
′
表示针对第j
′
个目标所有接近时刻中对应的第i
′
个时刻；m为保留的所有目标数量，nj′
为第j
′
个目标对应的所有的接近时刻数目；
[0072]
并行计算的原理是将数据传入gpu进行计算，基于gpu的结构特点，在面向简单的计算问题时计算效率相较于cpu串行计算可以得到显著提高。在步骤二中由于采用的是对每个目标的遍历预报，需要大量的计算，花费的时间较长，但线性j2模型和box模型复杂度小，且不同的预报节点之间彼此没有关联，十分契合gpu计算的特点。因此考虑将初筛模块通过cuda平台送入gpu中进行并行加速，计算完成后将计算后得到的所有存在接近风险的目标及对应的时刻再传回cpu端进行储存。
[0073]
在gpu端每个线程块(block)内的具体策略如图4所示；
[0074]
从图3可以看出每个线程块内，共有(m 1)
×
(n 1)个线程(thread)同时工作，每个线程分别单独分析某个目标在某一个预报时刻的接近情况。
[0075]
由此可以得到步骤二中的整个初筛方法如图5所示。
[0076]
其它步骤及参数与具体实施方式一或二相同。
[0077]
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤三中将保留目标的每个接近时刻前后分别扩充δt得到预报时段，利用sgp4模型对所有目标在预报时段内进行预报，得到最小接近距离，保留最小接近距离小于等于最小距离阈值的目标，利用hpop模型对保留的目标进行递推，得到最小接近距离；具体过程为：
[0078]
首先对步骤二筛选出的存在接近风险的目标的每个接近时刻进行前后扩充，这是由于线性j2预报存在较大误差，因此对每个时刻进行适当地扩充δt，可以得到一个新的预报时段从而有效地降低虚报、漏报：
[0079][0080]
sgp4模型预报简化了空间目标在运行过程中的复杂摄动项，考虑了大气摄动、四阶位势谐波、半同步和同步的自旋共振以及日月引力影响，是一种精度较高的解析模型，适用于近地目标。
[0081]
利用sgp4模型对式(10)得到的每个时段进行预报，得到每个目标最小接近距离及对应的时刻然后利用：
[0082][0083]
其中r
min
为设置的最小距离阈值；
[0084]
将最小接近距离满足上式的所有目标保留，称为危险目标；
[0085]
此时剩余的目标数目已经极少，对危险目标利用hpop高精度模型从目标航天器的历元时刻进行轨道递推，得到最小的接近距离及对应的接近时刻。
[0086]
其它步骤及参数与具体实施方式一至三之一相同。
[0087]
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述步骤四中利用步骤三中保留目标的历史tle数据生成初始时刻的协方差矩阵，并在二体模型下线性外推到接近时刻，得到协方差矩阵，结合最小接近距离利用pc方法的解析公式得到碰撞概率，与概率阈值比较后，对大于阈值的目标给出碰撞警告；具体过程为：
[0088]
由于本方法针对的是近地轨道目标，运动速度较快，因此考虑目标航天器和待筛选航天器之间为线性相对运动。在具体的计算中将两目标状态分布合并进行考虑，给出在接近时刻相对位置的联合概率分布，且假设空间目标的状态误差分布为正态分布，此时联合误差协方差矩阵等于两目标协方差矩阵之和。
[0089]
在相遇坐标系下，假设空间目标的碰撞形状均为圆形，将两个空间目标的体积集中到原点上，在坐标原点形成一个半径为r的联合球体，同时将两个目标的位置不确定性集中到另一个空间目标上，形成一个联合误差椭球分布p。通过对概率密度函数在联合球体中的体积进行积分就可以得到两个空间目标之间的碰撞概率pc：
[0090][0091]
为了提高计算效率，将两目标的误差椭球向相遇平面进行投影后，再利用rician分布的积分变换进行化简，然后可以得到一个无穷级数形式，取首项作为碰撞概率计算的解析公式进行计算，具体形式为：
[0092][0093]
式中，p0表示无穷级数的首项，σ
x
表示相遇平面坐标系x向标准差，σy表示相遇平面坐标系y向标准差，μ
x
表示最小接近距离在相遇平面坐标系x向投影，μy表示最小接近距离在相遇平面坐标系y向投影，r表示碰撞半径；
[0094]
将碰撞概率与概率阈值比较后，对大于阈值的目标给出碰撞警告；
[0095]
所述标准差σ
x
、σy从协方差矩阵当中得到，标准差是协方差矩阵对角线上的元素，协方差矩阵获取过程为：
[0096]
可以看出上式碰撞概率的解析形式仅与位置标准差、最小接近距离以及碰撞半径相关。下面是协方差矩阵的计算方法。
[0097]
利用步骤三中保留目标的历史tle数据生成初始时刻的协方差矩阵；
[0098]
由于本方法考虑线性相对运动，只需要空间位置的协方差矩阵。考虑利用步骤三中保留空间目标前10天的历史tle数据递推到历元初始时刻形成位置分布椭球，由此可以得到惯性系下的初始协方差矩阵
[0099]
在二体模型下采用线性协方差方法进行外推，可以得到接近时刻的协方差矩阵，表达式为：
[0100][0101]
式中，φ(t,t0)为二体模型下从初始时刻到预报时刻的状态转移矩阵(参见reynolds r g.direct solution of the keplerian state transition matrix[j].journal of guidance control and dynamics,2022.doi:10.2514/1.g006373.)；此时外推得到的协方差矩阵在惯性系下，还需要将接近时刻的协方差矩阵转换到相遇系下进行分析：
[0102][0103]
式中，和分别表示惯性系和相遇系下的协方差矩阵，m为地心惯性系到相遇坐标系的坐标变换矩阵；
[0104]
由此可以得到步骤三中危险目标在最小接近时刻相遇系下的协方差矩阵。
[0105]
由步骤三和步骤四组成的精筛方法如图6所示；
[0106]
综上述，可以得到整个筛选方法的逻辑流程为图7所示。
[0107]
其它步骤及参数与具体实施方式一至四之一相同。
[0108]
采用以下实施例验证本发明的有益效果：
[0109]
实施例一：
[0110]
需要注意，筛选中使用的hpop高精度轨道预报器，其考虑了地球的引力、非球形摄动、大气阻力摄动、太阳光压摄动和日月引力摄动，在不考虑定轨误差下，3天内预报误差很小。具体的参数为：大气阻力系数cd＝2.2，太阳光压系数cr＝1.0，面质比均为0.02m2/kg，选择j2000地心惯性坐标系。
[0111]
(1)首先对步骤二中提出的基于迹向误差修正的box法的有效性进行验证：
[0112]
从celestrak网站申请了450组近地轨道数据，取其中250组数据作为基于迹向误差修正的box方法的原始数据。利用线性j2预报与hpop模型积分结果进行对比，得到各组数据3天内在轨道系下的迹向误差最大值，以步骤二引入的自变量x为横轴，误差最大值为纵轴，结果如图8所示：
[0113]
从图8结果可以看出三天内迹向误差发散最大可达150km，且波动区间较大，不易设置碰撞阈值，利用三次多项式对迹向误差进行修正，可以得到归一化后的多项式为：
[0114][0115]
利用剩下的另外200组数据对修正结果进行验证，通过上式作迹向误差修正后的打靶仿真，可以得到结果如图9所示；
[0116]
从图9结果可以看出修正后三天内迹向位置误差最大值为55km，且大部分误差集中于20km以下。说明基于迹向误差修正的box方法相较于之前，碰撞阈值能够得到明显地降低。
[0117]
(2)然后对并行计算前后初筛模块的计算效率进行对比:
[0118]
从spacetrck网站下载了从2022年4月5日到4月12日之间23000组空间目标数据和700km轨道的surcal 150b作为初始数据进行一星对多星的筛选仿真，仿真参数设置如下：预筛选模型的筛选高度阈值δh＝80km，初筛步长为1s，box模型δrr＝50km，δrs＝130km，δrw＝30km。
[0119]
仿真的环境为：
[0120]
操作系统：arm-linux ubuntu 18.04.4
[0121]
gpu：nvidia tesla v100 pcle 32gb
[0122]
仿真软件：matlab2020a cuda10.1
[0123]
以总步数＝星数
×
步数，在预筛选之后还剩5199个待筛选目标，因此整个初筛过程需要并行的总步数为1.3476
×
109，在两种仿真环境下，分别对比初筛过程在加入gpu并行计算前后，仿真结果如表1所示：
[0124]
表1不同仿真条件下加速效果
[0125][0126]
从上述结果可以看出，并行计算加入后可以对整个初筛过程计算效率的提升为3684.53倍，加速效果明显。
[0127]
(3)最后对筛选方法的有效性进行验证：
[0128]
以2009年2月10日美俄卫星相撞事件为背景，从celestrak网站申请了美国卫星iridium-33(norad编号：24946)和俄罗斯报废卫星cosmos-2251(norad编号：22675)在碰撞之前的历史数据，具体为从2009年1月28日到2009年2月9日的数据，两颗卫星的轨道高度约为790km，以1月28日到2月8日的数据作为历史数据，以2月9日的数据作为初始数据。
[0129]
以cosmos-2251作为目标航天器，将iridium-33的初始tle数据放入1268组tle数据中作为筛选目标，所选数据为和两个碰撞目标同时期在轨，且轨道高度为700km左右的一批卫星。
[0130]
仿真参数设置如下：预筛选模型的筛选高度阈值δh＝80km，初筛步长为1s，box模型δrr＝50km，δrs＝130km，δrw＝30km，精筛步长选为0.02s，扩充时间段δt＝
±
15s，精筛距离阈值r
min
＝2km，假设碰撞体积均为5m。
[0131]
最接近时刻为2009年2月10日16:55:59.82(utc)，与实际碰撞时刻误差δt≤0.1s。通过历史数据递推到历元时刻分别得到位置的初始协方差矩阵如下：
[0132]
iridium-33：
[0133][0134]
cosmos-2251：
[0135][0136]
然后通过初始协方差矩阵外推到最接近时刻可以得到预报时刻的两个航天器的位置协方差矩阵为：
[0137]
iridium-33：
[0138][0139]
cosmos-2251：
[0140][0141]
可以计算得到碰撞概率为：pc＝5.0392
×
10-4
，大于红色碰撞预警阈值，有极大可能发生碰撞，与事实相符。由此也可说明发明提出的筛选方法是正确有效的。
[0142]
本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。