一种基于gpu加速的空间碎片预警方法
【专利摘要】一种基于GPU加速的空间碎片预警方法，步骤如下：(1)在GPU上移植SGP4/SDP4模型，在CPU上初始化SGP4/SDP4模型参数并分配GPU的内存，设置GPU的线程块的数量以及每个线程块的线程数；(2)在CPU上，获取当前预警的空间碎片的TLE数据，并将获取的数据复制到上述分配的GPU内存中；(3)在CPU上初始化零拷贝内存指针，将该指针映射到上述GPU内存上，根据当前预警的空间碎片数量结合GPU的可用最大内存以及线程数，计算空间碎片的批次大小；(4)根据上述批次大小，使用零拷贝技术对空间碎片进行分批，以流水线的方式得到每个碎片在每个时间点的碰撞概率并回传CPU并存储；(5)所有空间碎片处理完成后，释放分配的GPU内存。
【专利说明】
一种基于GPU加速的空间碎片预警方法
技术领域
[0001]本发明涉及一种空间碎片预警方法，更特别地说，是利用英伟达公司的GPU通用计 算平台开发一种并行算法进行空间碎片预警计算。
【背景技术】
[0002] 空间碎片预警，即预分析在轨空间碎片与目标航天器位置关系，判断碎片与航天 器碰撞的危险程度是进行空间活动的先决条件。
[0003] 然而在短时间内完成常规预警任务需要巨大的计算资源。使用高性能计算设备， 开发高性能算法来加速空间碎片预警计算是必然的选择。目前国内外大都选择多台服务器 多个计算核心组成计算机集群来做大规模空间碎片碰撞预警计算。
[0004] 中国科学院的杨大陶在其公开的论文《空间碎片预警并行计算系统设计与实现》 中公开了一种空间碎片预警算法并完成开发了预警系统的设计开发。该系统使用曙光公司 TC2600高性能计算设备，内置八台CB85-F刀片服务器共128个计算核心组成计算机机群对 我国航天器进行72小时空间碎片预警计算。其工作包括对空间物体的接近分析和目标筛选 (使用几何筛选法和时间筛选法），空间物体的轨道预报(使用SGP4/SDP4模型）和误差计算 (全轨道周期误差生产）。该算法较好的解决了短时间内完成空间碎片碰撞预警的问题。但 是计算代价巨大，功耗高，效费比低，无法用于大规模普及应用。
【发明内容】
[0005] 本发明的技术解决问题是:在保证计算性能的同时，大幅度降低计算成本，降低计 算功耗，提出一种新的基于GPU和SGP4/SDP4模型的空间碎片预警算法。
[0006] 本发明的技术解决方案是:一种基于GPU加速的空间碎片预警方法，步骤如下：
[0007] (1)在GPU上移植SGP4/SDP4模型，在CPU上初始化SGP4/SDP4模型参数并分配GPU的 内存，设置GHJ的线程块的数量以及每个线程块的线程数；
[0008] (2)在CPU上，获取当前预警的空间碎片的TLE数据，并将获取的数据复制到上述分 配的GHJ内存中；
[0009] (3)在CPU上初始化零拷贝内存指针，将该指针映射到上述GPU内存上，根据当前预 警的空间碎片数量结合GPU的可用最大内存以及线程数，计算空间碎片的批次大小；
[0010] (4)根据上述批次大小，使用零拷贝技术对空间碎片进行分批，以流水线的方式执 行如下处理：
[0011] (4.1)将每个线程块代表一个碎片，每个时间点对应一个线程；
[0012] (4.2)利用GPU上的SGP4/SDP4模型计算每个碎片每个时间点的轨道参数，进而该 碎片在该时间点的碰撞概率；
[0013 ] (4.3)将每个碎片在每个时间点的碰撞概率回传 CPU并存储；
[0014] (5)所有空间碎片处理完成后，释放分配的GPU内存。
[0015] 所述步骤⑶中批次大小& = Y X 士/ (HW);
[0016] 其中，空间碎片总数为N个，以△ t为时间步长，计算t时间内的所有可能碰撞概率 事件;总共开启m个线程块，每个线程块开启η条线程。
[0017] 线程块m的最优范围为512-2048之间且为16的倍数。
[0018]线程η的最优范围为64-512之间且为64的倍数。
[0019]对于TLE两行轨道要素数据的处理采用字符串处理的方式，同时利用宏定义的方 式减少动态分配内存。
[0020] 本发明与现有技术相比有益效果为：
[0021] (1)本发明将传统的SGP4/SDP4轨道预报模型移植到NVIDA的GPU平台，然后使用块 分解策略，将所有空间碎片分成多个批次作流处理，利用GHJ的零拷贝内存将计算和数据传 输操作重叠，从而隐藏内存读取延迟，获得更快的加速比。
[0022] (2)由于一个GPU的成本远远低于一台服务器的价格，本发明与现有开发的算法相 比成本可以极大的降低。一个较高性能的GPU如Tesla C2075便可以在五六百s内完成我国 所有空间目标(有TLE编目的目标）的空间碎片预警计算，加速比最大可达25倍左右（见图 3)，完全满足常规预警任务的时效性要求。具有很高的效费比。
[0023] (3)通过合理设置线程块和每个线程块线程数的小，既可以充分利用GPU上的所有 寄存器，避免每个流处理器上的计算资源浪费，也能获得最大的计算效率而不需要把计算 空间碎片预警计算中的一部分存到显卡内存而降低执行效率。
【附图说明】
[0024]图1本发明流程图；
[0025] 图2本发明块分解策略示意图（X =最优线程数目，y =可用内存）；
[0026] 图3为CPU和GPU计算时长对比曲线；
[0027]图4为本发明方法处理实施例的加速比曲线。
【具体实施方式】
[0028] 下面将结合附图和实施例对本发明做进一步的详细说明。
[0029] 本发明的空间碎片预警算法主要包括两大部分:1)精确轨道预报和考虑误差的碰 撞概率计算，2)使用块分解策略设计基于GPU的预警并行算法。
[0030] 步骤一，将SGP4/SDP4模型和概率计算模型移植至GPU;
[0031] 由于在GPU上运行的程序是基于英伟达公司的⑶DA平台编写的，现场的算法不能 直接应用于GPU并行计算，首先要做的工作是将SGP4/SDP4模型和碰撞概率计算模型移植到 GPU平台。其碰撞预警计算算法包括基本分两个部分，包括精确轨道预报和考虑误差的碰撞 概率计算。首先通过TLE数据库计算空间碎片与航天器的接近距离、接近时刻和此时各自的 位置速度矢量和相对速度。另外，由于轨道预报和数据库本身数据存在误差，因此需要在精 确数据计算的基础上，考虑轨道预报误差进行碰撞概率计算。
[0032] 步骤1 · I :SGP4/SDP4是由北美防空司令部（North American Aerospace Defense CommancUNORAD)开发用于轨道预报的模型，结合美国全球观测网的观测资料生产了全球最 大的空间目标编目数据库，并以两行轨道要素(Two-line element，TLE)形式发布。配套使 用TLE根数与SGP4/SDP4模型可以获得很好的轨道预报精度。其中SGP4由Ken Cranford于 1970年开发，适用于中低轨目标，在大气摄动修正模型中运用含二次方项的密度函数;SDP4 为SGP4的外延，由Hujsak于1979年开发，加入了日月引力项和地球引力势函数中的某些扇 形和田形谐函数项，专用于深空轨道预报。应用时，首先对卫星轨道根数进行修正：
[0033]
[0034]
[0035]
[0036]
[0037]
[0038] H [0039]式中为修正后的轨道半长轴，为修正后的周期频率。
[0040]然后根据卫星高度对大气阻力进行相应的修正。最后将卫星轨道参数代入Kepler 方程迭代获得卫星运动的轨道方程，带入相应的轨道要素及单位向量便可确定卫星的位置 和速度矢量：
[0041]
[0042] (2)
[0043] 式中r和V为天球坐标系下的位置和速度矢量，U、V为单位向量，由下式确定：
[0044]
[0045]
[0046]
[0047]
[0048]对于SDP4模型，计算方法基本与SGP4-致，只是在常量、大气及引力模型上有所不 同。
[0049] 步骤1.2:由于航天器与空间碎片的轨道预报存在误差，因此存在一定概率的碰撞 可能。碰撞概率的计算基于定轨得出的轨道数据和误差，即采用轨道预报算法和相应的概 率理论。
[0050] 为了计算碰撞概率，首先需要获得以下4组数据:航天器位置矢量η;碎片位置矢 量K;航天器位置误差矩阵A1;碎片位置误差矩阵A 2。其中前两项通过SGP4/SDP4轨道预报算 法获得，后两项通过结合TLE历史真实数据和SGP4/SDP4预报数据获得。
[0051] 这里假设航天器和空间碎片的位置误差满足三维正态分布，有
[0052]
[0053] 式中，a是一个由真实数据确定的三维向量，误差矩阵
?(1)指数部 分为一椭球方程，用这样的概率密度椭球来描述空间物体分布，定义为误差椭球。
[0054]碰撞概率的计算在交会坐标系内进行，交会参考系原心在航天器上，ζ轴为相对速 度矢量方向，与相对速度矢量垂直的平面即为交会平面，交会平面上X轴为复合位置误差椭 球的主轴在交会平面上的投影，y轴由右手规则确定。
[0055] 在交会平面内，航天器和碎片的位置满足正态分布，其概率密度表示为
[0056] P)
[0057] 式中XQ，yQ分别为航天器与碎片的相对位置关系，〇x，〇y为相对位置标准差。假定航 天器的安全半径为R，则碰撞概率就是此半径内对概率密度函数的积分：
[0058]
(3)
[0059] 将计算的碰撞概率和预先设定的碰撞概率值进行比较，当碰撞概率值小于设定的 预警值时(例如10-5)，则说明航天器与此碎片无碰撞风险，反之则说明航天器与碎片有危 险交会，需要进一步精确测轨和计算以确定是否采取规避措施。
[0060] 由于美国国防部公布的SGP4/SDP4模型只有对应的FORTRAN和C++串行版本，其中 用到的许多库都无法直接用于GHJ计算，特别是目前CUDA不支持C++的string类，对TLE两行 轨道要素数据的处理采用能够在GPU上运行的字符串处理算法，同时在GPU上动态分配内存 会极大的降低GHJ执行速度，而TLE数据量巨大，因此在TLE数据处理上尽量少用动态分配内 存，以大量使用宏定义来代替。另外对于程序异常的处理代码也应相应地改为符合CUDA编 程规范。
[0061] 步骤二，设计高效的空间碎片预警并行算法；
[0062]步骤二是本发明的核心部分，步骤一只是将模型移植与GPU，而步骤二是设计高效 低成本算法的关键。GPU以CUDA(compute Unified Device Architecture)作为新的并行计 算计算设备软件体系。CUDA的编程模型将CPU作为终端(Host)，负责逻辑性较强的事务和串 行计算，而GPU作为设备(Device)或协处理器(Coprocessor )，从而让GPU运行一些能被高度 线程化的程序。运行在GPU上的CUDA并行计算函数定义为核函数Kernel〈〈〈DimGrid， DimBlock>>>(Parameters)，其中DimGrid为并行计算的线程块(block)数目，DimBlock为每 个线程块包含的线程（thread)数。同一线程块内的线程可以通过共享存储器（Shared Memory)进行共享内存和同步实现通信。为便于并行控制，CUDA使用了四个内建变量 (blockIdx，gridDim, threadIdx，blockDim)来描述线程块和线程格的维数和索引。
[0063]为充分发挥GPU性能，软件架构设计如下:每个线程块代表一个碎片，每个线程计 算一个碎片一个时间步长的碰撞概率计算。每个线程块内的所有线程共享同一个空间碎片 的初始TLE参数。同时使用零拷贝（Zero-Copy)技术进行CPU与GPU之间的数据传输。因此对 GPU上内存解引用时，如果是基于GPU的，那么可以获得全局内存的高速带宽（144GB/S)。对 于空间碎片碰撞预警计算这种高度计算密集型任务，零拷贝技术可以极大的节省设备显式 传输的时间。事实上，使用零拷贝，将传输和内核操作分成更小的块，然后以流水线方式执 行密集型的计算任务相当于将计算与数据传输操作重叠，隐藏了 PCI-E传输延迟从而获得 计算性能的极大提升。
[0064] 基于空间碎片预警计算所具有的并行性，利用GPU单指令多数据流的并行优势，设 计一种计算碎片预警的并行异构算法。其中SGP4/SDP4模型轨道预报和碎片碰撞的概率计 算在GPU上完成。而CPU负责初始化SGP4/SDP4模型参数，从文本文件获取主航天器和碎片初 始时刻的TLE参数，定义仿真的初始条件并管理调度与GPU之间的数据传输。使用GPU的块分 解技术，从而允许进行任意数量的空间碎片任意长时间的仿真。为了便于应用于其他类似 问题，将算法设计得尽可能通用。
[0065] 图1详细描述了空间碎片轨道通用计算模型。算法首先在CPU上初始化碎片TLE参 数，SGP4/SDP4模型参数，然后使用块分解策略将所有空间碎片数分成多批次。算法每次循 环处理一定数量的空间碎片在一定时间内的碰撞概率计算。块分解策略如图2所示，每个块 的分配大小取决于两个参数，即每个块最优线程数目和可用的最大内存，通常由GPU的性能 参数确定。
[0066] 每个线程块(block)处理一个空间碎片，每个线程(thread)块下的线程处理这个 空间碎片某个时间点与主航天器危险交会的计算。因此若假设共有空间碎片N个，以△ t为 时间步长，计算t时间内的所有可能碰撞概率事件。总共开启m个线程块，每个线程块开启η 条线程，则每个批次能处理mXη条线程，需要分成的处理批次是? = 每个批次 〇 每个空间碎片的TLE参数通过⑶DA的零拷贝内存技术复制到GPU。一个批次的每个核函数 (kernel)被启动后进入循环开始碎片与主航天器的碰撞概率计算。由于块分解策略的使 用，使得的空间碎片数量以及轨道预报时间长度不再被GHJ内存所限制，因此大规模数量长 时间的仿真成为可能。具体步骤概括如下：
[0067] (1)在CPU上初始化SGP4/SDP4模型参数并分配GPU的内存，设置GPU的线程块的数 量m以及每个线程块的线程数η;线程块m和线程η的确定，由于需要计算的可编目空间碎片 的数量高达四五千个，同时计算碰撞概率过程较为复杂，每个线程块里的每个线程占用较 多的寄存器，根据GPU线程分配原则和仿真试验结果，随GPU计算能力不同，最优范围m范围 为512-2048之间且为16的倍数，最优范围η范围为64-512之间且为64的倍数。
[0068] (2)在CPU上，获取当前预警的空间碎片的TLE数据，并将获取的数据复制到上述分 配的GHJ内存中；
[0069] (3)在CPU上初始化零拷贝内存指针，将该指针映射到上述GPU内存上，根据当前预 警的空间碎片数量结合GPU的可用最大内存以及线程数，计算空间碎片的批次大小；
[0070] (4)根据上述批次大小，使用零拷贝技术对空间碎片进行分批，以流水线的方式执 行如下处理：
[0071] (4.1)将每个线程块代表一个碎片，每个时间点对应一个线程；
[0072] (4.2)利用GPU上的SGP4/SDP4模型计算每个碎片每个时间点的轨道参数，进而该 碎片在该时间点的碰撞概率；
[0073 ] (4.3)将每个碎片在每个时间点的碰撞概率回传 CPU并存储；
[0074] (5)所有空间碎片处理完成后，释放分配的GPU内存。
[0075] 实施例1
[0076] 本发明以对单个航天器为例计算与空间碎片群的碰撞概率并进行精度和加速性 能分析。首先采用的GHJ为NVIDIA的Tesla C2075,共有448个加速核心，具体配置如表1所示 [0077]表1数值仿真硬件配置
[0078]表1数值仿真硬件配置
[0080] 选取要进行预警的目标航天器为空间实验室"天宫一号"，计算天宫一号与2009年 美俄卫星相撞形成的碎片群（共1382个碎片）可能发生碰撞的概率。其TLE数据从 CelesTrack网站(http://www. celestrak.com/)下载，预警计算时间从2016年1月15号零点 开始，时间步长为60s，共2.5E5s (大约三天）。
[0081] 如图3、4所示，当数据规模小于10X50S时，采用并行算法比采用CPU传统计算方法 的速度要慢，这是由于数据量较少，GPU设备的启动、与CPU之间的数据传输通信消耗的时间 远远大于数据处理的耗时，无法体现GPU的数据处理能力。当仿真时长大于50s时(数据规模 为10 X 50s)，GPU处理速度与CPU相比明显加快，加速比快速增加，当数据规模进一步变大时 (10 XI. 5E5s)，加速比基本不再变化，说明GPU计算能力已经饱和，由图4可知，此时最大加 速比可达25倍左右。而《空间碎片预警并行计算系统的设计与实现》用曙光公司TC2600高性 能计算设备开发的预警算法也只达到约15倍的加速比。相比较本发明的算法具有更高的效 率而只需要更便宜的计算设备。这样，对天宫一号作碰撞预警计算时长2.5E5秒时间只需要 花费3.787s，理论上完成我国所有编目航天器约150个只需要568.05s。
[0082]本发明未详细说明部分属于本领域技术人员公知常识。
【主权项】
1. 一种基于GPU加速的空间碎片预警方法，其特征在于步骤如下： (1) 在GPU上移植SGP4/SDP4模型，在CHJ上初始化SGP4/SDP4模型参数并分配GPU的内 存，设置GPU的线程块的数量以及每个线程块的线程数； (2) 在CPU上，获取当前预警的空间碎片的TLE数据，并将获取的数据复制到上述分配的 GHJ内存中； (3) 在CPU上初始化零拷贝内存指针，将该指针映射到上述GPU内存上，根据当前预警的 空间碎片数量结合GHJ的可用最大内存以及线程数，计算空间碎片的批次大小； (4) 根据上述批次大小，使用零拷贝技术对空间碎片进行分批，以流水线的方式执行如 下处理： (4.1) 将每个线程块代表一个碎片，每个时间点对应一个线程； (4.2) 利用GPU上的SGP4/SDP4模型计算每个碎片每个时间点的轨道参数，进而该碎片 在该时间点的碰撞概率； (4.3) 将每个碎片在每个时间点的碰撞概率回传 CRJ并存储； (5) 所有空间碎片处理完成后，释放分配的GRJ内存。2. 根据权利要求1所述的一种基于GPU加速的空间碎片预警方法，其特征在于:所述步 骤(3)中批次其中，空间碎片总数为N个，以At为时间步长，计算t时间内的所有可能碰撞概率事件； 总共开启m个线程块，每个线程块开启n条线程。3. 根据权利要求2所述的一种基于GPU加速的空间碎片预警方法，其特征在于:线程块m 的最优范围为512-2048之间且为16的倍数。4. 根据权利要求2所述的一种基于GPU加速的空间碎片预警方法，其特征在于:线程n的 最优范围为64-512之间且为64的倍数。5. 根据权利要求1所述的一种基于GPU加速的空间碎片预警方法，其特征在于:对于TLE 两行轨道要素数据的处理采用字符串处理的方式，同时利用宏定义的方式减少动态分配内 存。
【文档编号】G06F17/50GK106055780SQ201610364134
【公开日】2016年10月26日
【申请日】2016年5月26日
【发明人】徐 明, 林名培, 胡肖肖, 何艳超, 付小宇, 罗通
【申请人】北京航空航天大学