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

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

Info

Publication number
CN114936471B
CN114936471B CN202210671501.7A CN202210671501A CN114936471B CN 114936471 B CN114936471 B CN 114936471B CN 202210671501 A CN202210671501 A CN 202210671501A CN 114936471 B CN114936471 B CN 114936471B
Authority
CN
China
Prior art keywords
target
targets
representing
screening
collision
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210671501.7A
Other languages
English (en)
Other versions
CN114936471A (zh
Inventor
张刚
卿晨
叶东
吴宝林
李化义
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202210671501.7A priority Critical patent/CN114936471B/zh
Publication of CN114936471A publication Critical patent/CN114936471A/zh
Application granted granted Critical
Publication of CN114936471B publication Critical patent/CN114936471B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Navigation (AREA)

Abstract

一种基于并行计算的航天器碰撞预警分层快速筛选方法，本发明涉及基于并行计算的航天器碰撞预警分层快速筛选方法。本发明的目的是为了解决现有筛选方法仅采用几何筛选和一些简单的二体摄动环境下的筛选不能够满足精度需求，而均使用高精度数值方法进行筛选会消耗大量时间，实用性较小的问题。过程为：一、利用预筛选模块从几何角度上筛除不可能接近的空间目标，得到保留的空间目标；二、保留存在接近风险的目标以及接近时刻；三、得到最小接近距离；四、得到碰撞概率，与概率阈值比较后，对大于阈值的目标给出碰撞警告。本发明用于航天器碰撞预警筛选领域。

Description

一种基于并行计算的航天器碰撞预警分层快速筛选方法
技术领域
本发明属于航天器碰撞预警筛选领域，具体涉及基于并行计算的航天器碰撞预警分层快速筛选方法。
背景技术
随着航天事业不断发展，各国的星座计划以及太空实验接连部署，越来越多的卫星被送入太空中，这也预示着太空中的空间碎片数目会急剧上升。据不完全统计，截至目前，在轨大小超过10cm的碎片约为34000个，大小在1cm到10cm之间的碎片约为90万，在1mm到1cm之间的碎片超过1.2亿。而目前由人类送入地球轨道的卫星超过10000颗，因此极有可能出现空间目标之间的发生碰撞情况，造成严重的后果，尤其对于近地轨道目标，具有较高的相对运动速度，产生的破坏性更为强烈。为了空间目标的在轨安全，避免造成不可挽回的后果，有必要研究如何快速筛选出对在轨航天器可能造成威胁的空间目标。
针对在轨航天器的碰撞预警问题，传统的筛选方法通常是针对绝不可能接近的目标，即轨道高度相差较大或相遇窗口距离相差较远的目标进行筛选，亦或是仅针对单星之间的高精度筛选。但由于实际任务中的空间目标数目庞大，运动情况也各不相同，仅采用几何筛选和一些简单的二体摄动环境下的筛选不能够满足精度需求，而均使用高精度数值方法进行筛选会消耗大量时间，实用性较小。因此，需要一种能够平衡筛选精度与筛选时间的筛选方法。
发明内容
本发明的目的是为了解决现有筛选方法仅采用几何筛选和一些简单的二体摄动环境下的筛选不能够满足精度需求，而均使用高精度数值方法进行筛选会消耗大量时间，实用性较小的问题，而提出一种基于并行计算的航天器碰撞预警分层快速筛选方法。
一种基于并行计算的航天器碰撞预警分层快速筛选方法具体过程为：
步骤一、以空间目标的TLE根数作为初始数据，首先解算得到所有目标的初始的轨道六根数和近远地点，然后利用预筛选模块从几何角度上筛除不可能接近的空间目标，得到保留的空间目标；
步骤二、以目标航天器的历元时刻作为初始时刻，将步骤一中得到的保留的空间目标送入GPU中，在GPU中首先利用线性J2模型在预报周期内进行遍历预报，得到所有目标在相应时间节点的空间位置，然后利用基于误差修正的BOX方法对所有目标进行筛选，保留存在接近风险的目标以及接近时刻，送回CPU进行储存；
步骤三、将保留目标的每个接近时刻前后分别扩充Δt得到预报时段，利用SGP4模型对所有目标在预报时段内进行预报，得到最小接近距离，保留最小接近距离小于等于最小距离阈值的目标，利用HPOP模型对保留的目标进行递推，得到最小接近距离；
步骤四、利用步骤三中保留目标的历史TLE数据生成初始时刻的协方差矩阵，并在二体模型下线性外推到接近时刻，得到协方差矩阵，结合最小接近距离利用Pc方法的解析公式得到碰撞概率，与概率阈值比较后，对大于阈值的目标给出碰撞警告。
本发明的有益效果：
本发明提出了一种基于并行计算的近地轨道碰撞预警快速筛选方法。在本发明中，首先利用最小二乘法得到一种基于迹向误差修正的BOX方法，相较于传统的BOX固定空间域的阈值设置，可以更好地降低阈值门限，减少虚警数，然后在此基础上搭建包括预筛选、初筛和精筛在内的完整的筛选模型，并通过GPU对初筛模块进行并行加速。采用本发明方法进行空间目标碰撞预警的筛选问题，只需要给出目标航天器和待筛选目标的初始TLE根数以及过去10天的历史TLE根数，即可快速筛选得到有碰撞风险的目标，并具有一定的筛选精度。
附图说明
图1为本发明基于并行计算的近地轨道碰撞预警快速筛选方法流程图；
图2为本发明筛选模块示意图；
图3为本发明迹向位置误差修正方法示意图；
图4为本发明每个线程块内并行结构示意图；
图5为本发明初筛过程示意图；
图6为本发明精筛过程示意图；
图7为本发明完整筛选流程示意图；
图8为修正前迹向位置误差分布情况图；
图9为本发明修正后迹向位置误差分布情况图。
具体实施方式
具体实施方式一：本实施方式一种基于并行计算的航天器碰撞预警分层快速筛选方法具体过程为：
本发明设计了一种基于并行计算的近地轨道碰撞预警快速筛选方法，能够对3天内目标航天器有威胁的目标进行快速筛选。该方法由预筛选、初筛和精筛三个筛选模块组成，适用于近地轨道所有的目标，在保证一定的筛选精度情况下有较高的计算效率，最后能够以碰撞概率的解析形式量化碰撞风险作为预警指标。
步骤一、以空间目标的TLE根数作为初始数据，首先解算得到所有目标的初始的轨道六根数和近远地点，然后利用预筛选模块从几何角度上筛除不可能接近的空间目标，得到保留的空间目标；
步骤二、以目标航天器的历元时刻作为初始时刻，将步骤一中得到的保留的空间目标送入GPU中，在GPU中首先利用线性J2模型在预报周期内(根据任务实际需要设定预报时长，通常为3天，不超过5天)进行遍历预报，得到所有目标在相应时间节点的空间位置，然后利用基于误差修正的BOX方法对所有目标进行筛选，保留存在接近风险的目标以及接近时刻，送回CPU进行储存；
步骤三、将保留目标的每个接近时刻前后分别扩充Δt得到预报时段，利用SGP4模型对所有目标在预报时段内进行预报，得到最小接近距离，保留最小接近距离小于等于最小距离阈值的目标，利用HPOP模型对保留的目标进行递推，得到最小接近距离；
步骤四、利用步骤三中保留目标的历史TLE数据生成初始时刻的协方差矩阵，并在二体模型下线性外推到接近时刻，得到协方差矩阵，结合最小接近距离利用Pc方法的解析公式得到碰撞概率，与概率阈值比较后，对大于阈值的目标给出碰撞警告。
具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤一中以空间目标的TLE根数作为初始数据，首先解算得到所有目标的初始的轨道六根数和近远地点，然后利用预筛选模块从几何角度上筛除不可能接近的空间目标，得到保留的空间目标；具体过程为：
整个筛选过程使用的惯性系均为J2000地心惯性坐标系；
空间中的主目标航天器是卫星，备选目标可以是卫星、空间碎片等，只要其在空间监视网的范围内，具有TLE根数即可；
TLE根数是北美空防空令部发布的人造天体的运行状态数据，不能直接用于轨道预报，需要使用SGP4模型进行解算；
首先使用SGP4模型对主目标和备选目标的TLE根数进行解算，得到主目标和备选目标历元时刻的初始轨道六根数和近远地点；
采用预筛选模块从几何角度上筛除不可能接近的空间目标，几何角度上筛除包含近远地点筛选(公式(1)(2))和倾角筛选(公式(3)(4)(5))两个部分；
对于给定主目标和备选目标的初始六根数分别为[a1,e1,I1,Ω1,ω1,f1]和[a2,e2,I2,Ω2,ω2,f2]，可以分别得到两个轨道的对应的近远地点
Figure BDA0003693355030000041
和
Figure BDA0003693355030000042
Figure BDA0003693355030000043
a1表示主目标对应的瞬时半长轴，e1表示主目标对应的离心率，I1表示主目标对应的目标轨道的倾角，Ω1表示主目标对应的升/降交点赤经，ω1表示主目标对应的近地点角距，f1表示主目标对应的真近点角；
a2表示备选目标对应的瞬时半长轴，e2表示备选目标对应的离心率，I2表示备选目标对应的目标轨道的倾角，Ω2表示备选目标对应的升/降交点赤经，ω2表示备选目标对应的近地点角距，f2表示备选目标对应的真近点角；
对轨道高度相差较大的轨道利用下式便以进行直接筛除，对满足下式的目标删除(比较
Figure BDA0003693355030000044
还有
Figure BDA0003693355030000045
和
Figure BDA0003693355030000046
之间与Δh的关系)，保留近远地点半径均在筛选阈值内的目标；
对满足下式的目标删除
Figure BDA0003693355030000047
式中，Δh为高度阈值，
Figure BDA0003693355030000048
表示主目标的远地点半径，
Figure BDA0003693355030000049
表示备选目标的远地点半径，
Figure BDA00036933550300000410
表示主目标的近地点半径，
Figure BDA00036933550300000411
表示备选目标的近地点半径，||表示或；
此外，对于轨道面之间倾角较大的轨道，接近时刻通常只有两个轨道面交线与轨道交点处，通过对两个交点之间高度差和距离阈值比较即可筛除不可能接近目标，记两个轨道面的角动量矢量分别为h1和h2，轨道面交线矢径为：
Figure BDA00036933550300000412
交线有两个共线不同向的矢径
Figure BDA00036933550300000413
和
Figure BDA00036933550300000414
分别对应于式(3)中的正负关系。得出矢径后，继续计算得到矢径分别和每条轨道交点的真近点角f±(正负号分别对应一组交点，一共有两组，)，分别对应不同的矢径：
Figure BDA0003693355030000051
式中，
Figure BDA0003693355030000052
是升/降交点的矢径，是在地心惯性坐标系下，由地心指向轨道升(降)交点的矢量，分别对应
Figure BDA0003693355030000053
和
Figure BDA0003693355030000054
(
Figure BDA0003693355030000055
对应
Figure BDA0003693355030000056
Figure BDA0003693355030000057
对应
Figure BDA0003693355030000058
)；|h3|表示交线矢径的模，通常规一化后取为1；
Figure BDA0003693355030000059
表示升/降交点矢径的模，规一化后取为1；ωi表示轨道近地点幅角；
从而可以通过式(5)求出两个轨道在同一交点处的轨道半径，作差后与高度阈值比较，将满足公式(6)的目标进行筛除：
Figure BDA00036933550300000510
|r1-r2|≥Δh (6)
式中，r1表示两个轨道在同一交点处时轨道1的轨道半径，r2表示两个轨道在同一交点处时轨道2的轨道半径。
其它步骤及参数与具体实施方式一相同。
具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤二中以目标航天器的历元时刻作为初始时刻，将步骤一中得到的保留的空间目标送入GPU中，在GPU中首先利用线性J2模型在预报周期内(根据任务实际需要设定预报时长，通常为3天，不超过5天)进行遍历预报，得到所有目标在相应时间节点的空间位置，然后利用基于误差修正的BOX方法对所有目标进行筛选，保留存在接近风险的目标以及接近时刻，送回CPU进行储存；
具体过程为：
将步骤一中得到的保留的空间目标，送入GPU中(通过CUDA平台送入GPU中进行并行加速)，在GPU中对步骤一保留的空间目标利用线性J2模型在预报周期内进行轨道递推，其中半长轴采用平均根数，得到空间目标各预报时刻的轨道六根数后，基于空间目标各预报时刻的轨道六根数转化得到对应的惯性系下空间位置，继而得到待筛选空间目标在主目标的轨道系下空间位置；
然后利用基于迹向误差修正的BOX方法对所有目标进行初步筛选，具体过程为：
该方法利用最小二乘法拟合的多项式对迹向误差进行修正。首先引入一个新的自变量对每个空间目标在轨道系下的位置误差的离散点进行线性化：
Figure BDA0003693355030000061
式中，rp表示初始近地点半径，x表示用于线性化的新的自变量，e表示离心率，I表示目标点纬度的倾角，Ω表示升/降交点赤经；
取N个近地空间目标作打靶仿真，得到3天内沿迹向的误差最大值Δr，从而可以得到一组数据序列(xj,Δrj)(j＝1,2,...,N)，使用拟合函数如下：
Figure BDA0003693355030000062
式中，
Figure BDA0003693355030000063
表示待拟合的函数，bk表示系数，
Figure BDA0003693355030000064
表示单项函数，n表示拟合函数包含单项函数的总项数，xj表示每个离散点的自变量大小，Δrj表示每个离散点位置误差的大小；j表示拟合的数据点共有N个，k表示拟合的多项式是从0到n阶；
在本方法中采用多项式进行拟合，即取
Figure BDA0003693355030000065
从而可以解得系数b0,b1,...,bn；具体的拟合流程如图3所示；
获取主目标的轨道系下待筛选空间目标中主目标和备选目标的迹向位置、径向位置和法向位置；
将主目标和备选目标的径向位置作差得到径向相对位置rR；
将主目标和备选目标的法向位置作差得到法向相对位置rW；
利用拟合得到的多项式
Figure BDA0003693355030000066
对轨道系下待测空间目标中的主目标和备选目标的迹向位置进行修正，得到修正后的迹向相对位置rS；具体过程为：
将待测空间目标中的主目标和备选目标的离散点数据xi带入多项式
Figure BDA0003693355030000067
得到待测空间目标中主目标和备选目标的离散点数据
Figure BDA0003693355030000068
将轨道系下的主目标和备选目标的迹向位置分别减去此时对应
Figure BDA0003693355030000069
得到修正后的主目标和备选目标的迹向位置，将修正后的主目标和备选目标的迹向位置作差得到修正后的迹向相对位置rS；
将迹向相对位置rS、径向相对位置rR、法向相对位置rW分别与设置的碰撞阈值进行比较，保留小于等于阈值的目标(满足式(9)的目标)：
Figure BDA0003693355030000071
其中ΔrR,ΔrS,ΔrW分别为轨道系下径向、迹向和法向的碰撞阈值，从而可以得到所有存在接近风险的目标(满足式(9)的目标)及对应的时刻
Figure BDA0003693355030000072
j′＝0,1,...,m；i′＝0,1,...,nj′；
其中j′表示针对所有保留目标中第j′个目标，i′表示针对第j′个目标所有接近时刻中对应的第i′个时刻；m为保留的所有目标数量，nj′为第j′个目标对应的所有的接近时刻数目；
并行计算的原理是将数据传入GPU进行计算，基于GPU的结构特点，在面向简单的计算问题时计算效率相较于CPU串行计算可以得到显著提高。在步骤二中由于采用的是对每个目标的遍历预报，需要大量的计算，花费的时间较长，但线性J2模型和BOX模型复杂度小，且不同的预报节点之间彼此没有关联，十分契合GPU计算的特点。因此考虑将初筛模块通过CUDA平台送入GPU中进行并行加速，计算完成后将计算后得到的所有存在接近风险的目标及对应的时刻
Figure BDA0003693355030000073
再传回CPU端进行储存。
在GPU端每个线程块(Block)内的具体策略如图4所示；
从图3可以看出每个线程块内，共有(m+1)×(n+1)个线程(Thread)同时工作，每个线程分别单独分析某个目标在某一个预报时刻的接近情况。
由此可以得到步骤二中的整个初筛方法如图5所示。
其它步骤及参数与具体实施方式一或二相同。
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤三中将保留目标的每个接近时刻前后分别扩充Δt得到预报时段，利用SGP4模型对所有目标在预报时段内进行预报，得到最小接近距离，保留最小接近距离小于等于最小距离阈值的目标，利用HPOP模型对保留的目标进行递推，得到最小接近距离；具体过程为：
首先对步骤二筛选出的存在接近风险的目标的每个接近时刻进行前后扩充，这是由于线性J2预报存在较大误差，因此对每个时刻进行适当地扩充Δt，可以得到一个新的预报时段
Figure BDA0003693355030000081
从而有效地降低虚报、漏报：
Figure BDA0003693355030000082
SGP4模型预报简化了空间目标在运行过程中的复杂摄动项，考虑了大气摄动、四阶位势谐波、半同步和同步的自旋共振以及日月引力影响，是一种精度较高的解析模型，适用于近地目标。
利用SGP4模型对式(10)得到的每个时段进行预报，得到每个目标最小接近距离
Figure BDA0003693355030000083
及对应的时刻
Figure BDA0003693355030000084
然后利用：
Figure BDA0003693355030000085
其中Rmin为设置的最小距离阈值；
将最小接近距离满足上式的所有目标保留，称为危险目标；
此时剩余的目标数目已经极少，对危险目标利用HPOP高精度模型从目标航天器的历元时刻进行轨道递推，得到最小的接近距离及对应的接近时刻。
其它步骤及参数与具体实施方式一至三之一相同。
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述步骤四中利用步骤三中保留目标的历史TLE数据生成初始时刻的协方差矩阵，并在二体模型下线性外推到接近时刻，得到协方差矩阵，结合最小接近距离利用Pc方法的解析公式得到碰撞概率，与概率阈值比较后，对大于阈值的目标给出碰撞警告；具体过程为：
由于本方法针对的是近地轨道目标，运动速度较快，因此考虑目标航天器和待筛选航天器之间为线性相对运动。在具体的计算中将两目标状态分布合并进行考虑，给出在接近时刻相对位置的联合概率分布，且假设空间目标的状态误差分布为正态分布，此时联合误差协方差矩阵等于两目标协方差矩阵之和。
在相遇坐标系下，假设空间目标的碰撞形状均为圆形，将两个空间目标的体积集中到原点上，在坐标原点形成一个半径为R的联合球体，同时将两个目标的位置不确定性集中到另一个空间目标上，形成一个联合误差椭球分布P。通过对概率密度函数在联合球体中的体积进行积分就可以得到两个空间目标之间的碰撞概率Pc：
Figure BDA0003693355030000086
为了提高计算效率，将两目标的误差椭球向相遇平面进行投影后，再利用Rician分布的积分变换进行化简，然后可以得到一个无穷级数形式，取首项作为碰撞概率计算的解析公式进行计算，具体形式为：
Figure BDA0003693355030000091
式中，P0表示无穷级数的首项，σx表示相遇平面坐标系x向标准差，σy表示相遇平面坐标系y向标准差，μx表示最小接近距离在相遇平面坐标系x向投影，μy表示最小接近距离在相遇平面坐标系y向投影，R表示碰撞半径；
将碰撞概率与概率阈值比较后，对大于阈值的目标给出碰撞警告；
所述标准差σx、σy从协方差矩阵当中得到，标准差是协方差矩阵对角线上的元素，协方差矩阵获取过程为：
可以看出上式碰撞概率的解析形式仅与位置标准差、最小接近距离以及碰撞半径相关。下面是协方差矩阵的计算方法。
利用步骤三中保留目标的历史TLE数据生成初始时刻的协方差矩阵；
由于本方法考虑线性相对运动，只需要空间位置的协方差矩阵。考虑利用步骤三中保留空间目标前10天的历史TLE数据递推到历元初始时刻形成位置分布椭球，由此可以得到惯性系下的初始协方差矩阵
Figure BDA0003693355030000092
在二体模型下采用线性协方差方法进行外推，可以得到接近时刻的协方差矩阵，表达式为：
Figure BDA0003693355030000093
式中，Φ(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.)；此时外推得到的协方差矩阵在惯性系下，还需要将接近时刻的协方差矩阵
Figure BDA0003693355030000094
转换到相遇系下进行分析：
Figure BDA0003693355030000095
式中，
Figure BDA0003693355030000096
和
Figure BDA0003693355030000097
分别表示惯性系和相遇系下的协方差矩阵，M为地心惯性系到相遇坐标系的坐标变换矩阵；
由此可以得到步骤三中危险目标在最小接近时刻相遇系下的协方差矩阵。
由步骤三和步骤四组成的精筛方法如图6所示；
综上述，可以得到整个筛选方法的逻辑流程为图7所示。
其它步骤及参数与具体实施方式一至四之一相同。
采用以下实施例验证本发明的有益效果：
实施例一：
需要注意，筛选中使用的HPOP高精度轨道预报器，其考虑了地球的引力、非球形摄动、大气阻力摄动、太阳光压摄动和日月引力摄动，在不考虑定轨误差下，3天内预报误差很小。具体的参数为：大气阻力系数Cd＝2.2，太阳光压系数Cr＝1.0，面质比均为0.02m2/kg，选择J2000地心惯性坐标系。
(1)首先对步骤二中提出的基于迹向误差修正的BOX法的有效性进行验证：
从CelesTrak网站申请了450组近地轨道数据，取其中250组数据作为基于迹向误差修正的BOX方法的原始数据。利用线性J2预报与HPOP模型积分结果进行对比，得到各组数据3天内在轨道系下的迹向误差最大值，以步骤二引入的自变量x为横轴，误差最大值为纵轴，结果如图8所示：
从图8结果可以看出三天内迹向误差发散最大可达150km，且波动区间较大，不易设置碰撞阈值，利用三次多项式对迹向误差进行修正，可以得到归一化后的多项式为：
Figure BDA0003693355030000101
利用剩下的另外200组数据对修正结果进行验证，通过上式作迹向误差修正后的打靶仿真，可以得到结果如图9所示；
从图9结果可以看出修正后三天内迹向位置误差最大值为55km，且大部分误差集中于20km以下。说明基于迹向误差修正的BOX方法相较于之前，碰撞阈值能够得到明显地降低。
(2)然后对并行计算前后初筛模块的计算效率进行对比:
从SpaceTrck网站下载了从2022年4月5日到4月12日之间23000组空间目标数据和700km轨道的SURCAL 150B作为初始数据进行一星对多星的筛选仿真，仿真参数设置如下：预筛选模型的筛选高度阈值Δh＝80km，初筛步长为1s，BOX模型ΔrR＝50km，ΔrS＝130km，ΔrW＝30km。
仿真的环境为：
操作系统：ARM-LINUX Ubuntu 18.04.4
GPU：Nvidia Tesla V100 Pcle 32GB
仿真软件：Matlab2020a+CUDA10.1
以总步数＝星数×步数，在预筛选之后还剩5199个待筛选目标，因此整个初筛过程需要并行的总步数为1.3476×109，在两种仿真环境下，分别对比初筛过程在加入GPU并行计算前后，仿真结果如表1所示：
表1不同仿真条件下加速效果
Figure BDA0003693355030000111
从上述结果可以看出，并行计算加入后可以对整个初筛过程计算效率的提升为3684.53倍，加速效果明显。
(3)最后对筛选方法的有效性进行验证：
以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日的数据作为初始数据。
以Cosmos-2251作为目标航天器，将Iridium-33的初始TLE数据放入1268组TLE数据中作为筛选目标，所选数据为和两个碰撞目标同时期在轨，且轨道高度为700km左右的一批卫星。
仿真参数设置如下：预筛选模型的筛选高度阈值Δh＝80km，初筛步长为1s，BOX模型ΔrR＝50km，ΔrS＝130km，ΔrW＝30km，精筛步长选为0.02s，扩充时间段Δt＝±15s，精筛距离阈值Rmin＝2km，假设碰撞体积均为5m。
最接近时刻为2009年2月10日16:55:59.82(UTC)，与实际碰撞时刻误差Δt≤0.1s。通过历史数据递推到历元时刻分别得到位置的初始协方差矩阵如下：
Iridium-33：
Figure BDA0003693355030000121
Cosmos-2251：
Figure BDA0003693355030000122
然后通过初始协方差矩阵外推到最接近时刻可以得到预报时刻的两个航天器的位置协方差矩阵为：
Iridium-33：
Figure BDA0003693355030000123
Cosmos-2251：
Figure BDA0003693355030000124
可以计算得到碰撞概率为：Pc＝5.0392×10-4，大于红色碰撞预警阈值，有极大可能发生碰撞，与事实相符。由此也可说明发明提出的筛选方法是正确有效的。
本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (3)

1.一种基于并行计算的航天器碰撞预警分层筛选方法，其特征在于：所述方法具体过程为：
步骤一、以空间目标的TLE根数作为初始数据，首先解算得到所有空间目标的初始的轨道六根数和近远地点，然后利用预筛选模块从几何角度上筛除不可能接近的空间目标，得到保留的空间目标；具体过程为：
使用的惯性系均为J2000地心惯性坐标系；
空间中的主目标航天器是卫星，备选目标是卫星、空间碎片；
使用SGP4模型对主目标和备选目标的TLE根数进行解算，得到主目标和备选目标历元时刻的初始轨道六根数和近远地点；
采用预筛选模块从几何角度上筛除不可能接近的空间目标，几何角度上筛除包含近远地点筛选和倾角筛选两个部分；
对于给定主目标和备选目标的初始轨道六根数分别为[a1,e1,I1,Ω1,ω1,f1]和[a2,e2,I2,Ω2,ω2,f2]，可以分别得到两个轨道的对应的近远地点
Figure FDA0004053921240000011
和
Figure FDA0004053921240000012
Figure FDA0004053921240000013
a1表示主目标对应的瞬时半长轴，e1表示主目标对应的离心率，I1表示主目标对应的目标轨道的倾角，Ω1表示主目标对应的升/降交点赤经，ω1表示主目标对应的近地点角距，f1表示主目标对应的真近点角；
a2表示备选目标对应的瞬时半长轴，e2表示备选目标对应的离心率，I2表示备选目标对应的目标轨道的倾角，Ω2表示备选目标对应的升/降交点赤经，ω2表示备选目标对应的近地点角距，f2表示备选目标对应的真近点角；
对满足下式的目标删除
Figure FDA0004053921240000014
式中，Δh为高度阈值，
Figure FDA0004053921240000015
表示主目标的远地点半径，
Figure FDA0004053921240000016
表示备选目标的远地点半径，
Figure FDA0004053921240000017
表示主目标的近地点半径，
Figure FDA0004053921240000018
表示备选目标的近地点半径，||表示或；
记两个轨道面的角动量矢量分别为h1和h2，轨道面交线矢径为：
Figure FDA0004053921240000021
继续计算真近点角f±：
Figure FDA0004053921240000022
式中，
Figure FDA0004053921240000023
是升/降交点的矢径分别对应
Figure FDA0004053921240000024
和
Figure FDA0004053921240000025
|h3|表示交线矢径的模；
Figure FDA0004053921240000026
表示升/降交点矢径的模；ωi表示轨道近地点幅角；
从而可以通过式(5)求出两个轨道在同一交点处的轨道半径，作差后与高度阈值比较，将满足公式(6)的目标进行筛除：
Figure FDA0004053921240000027
|r1-r2|≥Δh (6)
式中，r1表示两个轨道在同一交点处时轨道1的轨道半径，r2表示两个轨道在同一交点处时轨道2的轨道半径；
步骤二、以目标航天器的历元时刻作为初始时刻，将步骤一中得到的保留的空间目标送入GPU中，在GPU中首先利用线性J2模型在预报周期内进行遍历预报，得到所有空间目标在相应时间节点的空间位置，然后利用基于误差修正的BOX方法对所有空间目标进行筛选，保留存在接近风险的目标以及接近时刻，送回CPU端进行储存；具体过程为：
将步骤一中得到的保留的空间目标，送入GPU中，在GPU中对步骤一保留的空间目标利用线性J2模型在预报周期内进行轨道递推，其中半长轴采用平均根数，得到空间目标各预报时刻的轨道六根数后，基于空间目标各预报时刻的轨道六根数转化得到对应的惯性系下空间位置，继而得到待筛选空间目标在主目标的轨道系下空间位置；
然后利用基于误差修正的BOX方法对所有目标进行初步筛选，具体过程为：
首先引入一个新的自变量对每个空间目标在轨道系下的位置误差的离散点进行线性化：
Figure FDA0004053921240000028
式中，rp表示初始近地点半径，
Figure FDA0004053921240000029
表示用于线性化的新的自变量，e表示离心率，I表示目标点纬度的倾角，Ω表示升/降交点赤经；
取N个近地空间目标，得到3天内沿迹向的误差最大值Δr，从而可以得到一组数据序列
Figure FDA0004053921240000031
使用拟合函数如下：
Figure FDA0004053921240000032
式中，
Figure FDA0004053921240000033
表示多项式，bk表示系数，
Figure FDA0004053921240000034
表示单项函数，n表示拟合函数包含单项函数的总项数，
Figure FDA0004053921240000035
表示每个离散点的自变量大小，Δrj表示每个离散点位置误差的大小；j表示拟合的数据点，共有N个，k表示拟合的多项式是从0到n阶；
采用多项式进行拟合，即取
Figure FDA0004053921240000036
从而可以解得系数b0,b1,...,bn；
获取主目标的轨道系下待筛选空间目标中主目标和备选目标的迹向位置、径向位置和法向位置；
将主目标和备选目标的径向位置作差得到径向相对位置rR；
将主目标和备选目标的法向位置作差得到法向相对位置rW；
利用多项式
Figure FDA0004053921240000037
对轨道系下待测空间目标中的主目标和备选目标的迹向位置进行修正，得到修正后的迹向相对位置rS；具体过程为：
将待测空间目标中的主目标和备选目标的离散点数据
Figure FDA0004053921240000038
带入多项式
Figure FDA0004053921240000039
得到待测空间目标中主目标和备选目标的离散点数据
Figure FDA00040539212400000310
将轨道系下的主目标和备选目标的迹向位置分别减去此时对应
Figure FDA00040539212400000311
得到修正后的主目标和备选目标的迹向位置，将修正后的主目标和备选目标的迹向位置作差得到修正后的迹向相对位置rS；
将迹向相对位置rS、径向相对位置r R、法向相对位置rW分别与设置的碰撞阈值进行比较，保留小于等于阈值的目标：
Figure FDA00040539212400000312
其中ΔrR,ΔrS,ΔrW分别为轨道系下径向、迹向和法向的碰撞阈值，从而可以得到所有存在接近风险的目标及对应的时刻
Figure FDA0004053921240000041
其中j′表示针对所有保留目标中第j′个目标，i′表示针对第j′个目标所有接近时刻中对应的第i′个时刻；m为保留的所有目标数量，nj′为第j′个目标对应的所有的接近时刻数目；
计算完成后将计算后得到的所有存在接近风险的目标及对应的时刻
Figure FDA0004053921240000042
再传回CPU端进行储存；
步骤三、将保留目标的每个接近时刻前后分别扩充Δt得到预报时段，利用SGP4模型对所有目标在预报时段内进行预报，得到最小接近距离，保留最小接近距离小于等于最小距离阈值的目标，利用HPOP模型对保留的目标进行递推，得到最小接近距离；
Δt为扩充时间段；
步骤四、利用步骤三中保留目标的历史TLE根数生成初始时刻的协方差矩阵，并在二体模型下线性外推到接近时刻，得到协方差矩阵，结合最小接近距离利用Pc方法的解析公式得到碰撞概率，与概率阈值比较后，对大于概率阈值的目标给出碰撞警告。
2.根据权利要求1所述一种基于并行计算的航天器碰撞预警分层筛选方法，其特征在于：所述步骤三中将保留目标的每个接近时刻前后分别扩充Δt得到预报时段，利用SGP4模型对所有目标在预报时段内进行预报，得到最小接近距离，保留最小接近距离小于等于最小距离阈值的目标，利用HPOP模型对保留的目标进行递推，得到最小接近距离；具体过程为：
首先对步骤二筛选出的存在接近风险的目标的每个接近时刻进行前后扩充，得到一个新的预报时段
Figure FDA0004053921240000043
Figure FDA0004053921240000044
利用SGP4模型对式(10)得到的每个时段进行预报，得到每个目标最小接近距离
Figure FDA0004053921240000045
及对应的时刻
Figure FDA0004053921240000046
然后利用：
Figure FDA0004053921240000047
其中Rmin为设置的最小距离阈值；
将最小接近距离满足上式的所有目标保留，称为危险目标；
对危险目标利用HPOP模型从目标航天器的历元时刻进行轨道递推，得到最小的接近距离及对应的接近时刻。
3.根据权利要求2所述一种基于并行计算的航天器碰撞预警分层筛选方法，其特征在于：所述步骤四中利用步骤三中保留目标的历史TLE根数生成初始时刻的协方差矩阵，并在二体模型下线性外推到接近时刻，得到协方差矩阵，结合最小接近距离利用Pc方法的解析公式得到碰撞概率，与概率阈值比较后，对大于概率阈值的目标给出碰撞警告；具体过程为：
碰撞概率具体形式为：
Figure FDA0004053921240000051
式中，P0表示无穷级数的首项，σx表示相遇平面坐标系x向标准差，σy表示相遇平面坐标系y向标准差，μx表示最小接近距离在相遇平面坐标系x向投影，μy表示最小接近距离在相遇平面坐标系y向投影，R表示碰撞半径；
将碰撞概率与概率阈值比较后，对大于概率阈值的目标给出碰撞警告；
所述标准差σx、σy从协方差矩阵当中得到，标准差是协方差矩阵对角线上的元素，协方差矩阵获取过程为：
利用步骤三中保留空间目标前10天的历史TLE根数递推到历元初始时刻形成位置分布椭球，由此可以得到惯性系下的初始协方差矩阵
Figure FDA0004053921240000052
在二体模型下采用线性协方差方法进行外推，可以得到接近时刻的协方差矩阵，表达式为：
Figure FDA0004053921240000053
式中，Φ(t,t0)为二体模型下从初始时刻到预报时刻的状态转移矩阵；此时外推得到的协方差矩阵在惯性系下，还需要将接近时刻的协方差矩阵
Figure FDA0004053921240000054
转换到相遇系下进行分析：
Figure FDA0004053921240000055
式中，
Figure FDA0004053921240000056
和
Figure FDA0004053921240000057
分别表示惯性系和相遇系下的协方差矩阵，M为地心惯性系到相遇坐标系的坐标变换矩阵；
由此可以得到步骤三中危险目标在最小接近时刻相遇系下的协方差矩阵。
CN202210671501.7A 2022-06-14 2022-06-14 一种基于并行计算的航天器碰撞预警分层快速筛选方法 Active CN114936471B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210671501.7A CN114936471B (zh) 2022-06-14 2022-06-14 一种基于并行计算的航天器碰撞预警分层快速筛选方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210671501.7A CN114936471B (zh) 2022-06-14 2022-06-14 一种基于并行计算的航天器碰撞预警分层快速筛选方法

Publications (2)

Publication Number Publication Date
CN114936471A CN114936471A (zh) 2022-08-23
CN114936471B true CN114936471B (zh) 2023-03-10

Family

ID=82865987

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210671501.7A Active CN114936471B (zh) 2022-06-14 2022-06-14 一种基于并行计算的航天器碰撞预警分层快速筛选方法

Country Status (1)

Country Link
CN (1) CN114936471B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115186521B (zh) * 2022-09-14 2023-01-13 中国科学院空天信息创新研究院 空间目标轨道根数模拟仿真方法、装置、设备及介质
CN115394126B (zh) * 2022-10-25 2023-01-20 北京开运联合信息技术集团股份有限公司 一种空间目标的碰撞预警方法及装置
CN117271127B (zh) * 2023-09-26 2024-03-26 中科星图测控技术股份有限公司 一种实时大规模空间目标碰撞预警系统及方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111127295A (zh) * 2019-11-22 2020-05-08 哈尔滨工业大学 基于gpu的sgp4轨道模型集成式并行方法
CN114327919A (zh) * 2022-03-14 2022-04-12 北京航天驭星科技有限公司 空间目标碰撞预警方法及系统
US11312512B1 (en) * 2019-03-04 2022-04-26 United States Of America As Represented By The Secretary Of The Air Force Early warning reentry system comprising high efficiency module for determining spacecraft reentry time

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114326774B (zh) * 2022-03-14 2022-05-24 北京航天驭星科技有限公司 航天器碰撞规避策略生成的方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11312512B1 (en) * 2019-03-04 2022-04-26 United States Of America As Represented By The Secretary Of The Air Force Early warning reentry system comprising high efficiency module for determining spacecraft reentry time
CN111127295A (zh) * 2019-11-22 2020-05-08 哈尔滨工业大学 基于gpu的sgp4轨道模型集成式并行方法
CN114327919A (zh) * 2022-03-14 2022-04-12 北京航天驭星科技有限公司 空间目标碰撞预警方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Midair collision risk when executing an incorrect approach during established on required navigation performance operations;Cody T. Nichols 等;《IEEE》;20171109;全文 *
卫星轨道递推的GPU集成式并行加速方法;孔繁泽 等;《哈尔滨工业大学学报 》;20210519;全文 *
基于TLE数据的空间目标大气再入预报;刘劲宏;《中国博士学位论文全文数据库工程科技Ⅱ辑 》;20220315;全文 *

Also Published As

Publication number Publication date
CN114936471A (zh) 2022-08-23

Similar Documents

Publication Publication Date Title
CN114936471B (zh) 一种基于并行计算的航天器碰撞预警分层快速筛选方法
CN108820260A (zh) 低轨航天器的中期轨道预报方法、装置、存储介质
CN104751012A (zh) 沿飞行弹道的扰动引力快速逼近方法
CN106679675A (zh) 一种基于相对测量信息的火星最终接近段自主导航方法
CN103645489A (zh) 一种航天器gnss单天线定姿方法
Wolf et al. Toward improved landing precision on Mars
Ahn et al. Noniterative instantaneous impact point prediction algorithm for launch operations
McGee et al. Mighty Eagle: the development and flight testing of an autonomous robotic lander test bed
CN111125874A (zh) 一种动平台高精度测轨预报方法
CN110231619B (zh) 基于恩克法的雷达交接时刻预报方法及装置
Prevereaud et al. Debris aerodynamic interactions during uncontrolled atmospheric reentry
Boone et al. Artificial debris collision risk following a catastrophic spacecraft mishap in lunar orbit
Padilla et al. Assessment of rarefied hypersonic aerodynamics modeling and windtunnel data
CN105183948B (zh) 一种基于二次反射的高精度卫星太阳光压摄动力建模方法
CN115392540A (zh) 一种用于月球轨道交会制导的快速预报方法
Xu et al. A method for calculating collision probability between space objects
CN115422699A (zh) 一种交互式地基空间目标监视传感器模拟仿真系统
Gordienko et al. On choosing a rational flight trajectory to the Moon
Ivanov et al. Computational tools for rarefied aerodynamics
Hazra Autonomous Guidance for Asteroid Descent using Successive Convex Optimization
Williams et al. Global Performance Characterization of the three burn trans-Earth injection maneuver sequence over the lunar nodal cycle
Aubuchon et al. Dynamic stability of the boeing cst-100
Starshak et al. Optical free-flight measurements using GPU-accelerated computer graphics
Johnson et al. Multiobjective optimization of earth-entry vehicle heat shields
CN111123980A (zh) 一种卫星飞临时刻与拍摄范围的计算方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant