孔繁澤，葉 東，柳子然，孫兆偉

(哈爾濱工業(yè)大學(xué) 航天學(xué)院,哈爾濱 150001)

據(jù)美國(guó)空間監(jiān)測(cè)網(wǎng)站給出的數(shù)據(jù)，截至2013年，地球軌道上尺寸大于10 cm的可編目空間物體數(shù)量已超過16 000多個(gè).針對(duì)日益增多的空間物體，在軌衛(wèi)星在進(jìn)行變軌機(jī)動(dòng)時(shí)，需對(duì)大量空間物體進(jìn)行高精度的軌道預(yù)報(bào)計(jì)算，以預(yù)防在軌航天器與其他空間天體的碰撞.由于衛(wèi)星高精度的軌道模型具有較高的計(jì)算復(fù)雜度，而目前計(jì)算機(jī)上常用的CPU(Central Processing Unit)處理器僅能進(jìn)行串行計(jì)算，進(jìn)行軌道預(yù)報(bào)需消耗大量的時(shí)間，難以滿足軍事上及民用上衛(wèi)星的快速響應(yīng)需求，因此有必要對(duì)軌道模型計(jì)算加速方法展開研究.如表1所示為各類星載計(jì)算機(jī)的CPU型號(hào)與性能參數(shù).

表1 On board computer(OBC)及CPU參數(shù)

對(duì)此，國(guó)內(nèi)外都進(jìn)行了相關(guān)的研究.計(jì)天陽(yáng)[1]采用模塊化加速方法及英偉達(dá)公司提供的CUDA庫(kù)，配套GPU(Graphics Processing Unit)對(duì)SGP4/SDP4軌道模型進(jìn)行了加速試驗(yàn)，填補(bǔ)了國(guó)內(nèi)在軌道模型GPU并行加速領(lǐng)域研究的空白. Liu等[2]利用GPU并行計(jì)算加速方法進(jìn)行大量空間碎片的軌道預(yù)報(bào)，通過PMC簡(jiǎn)化模型(The global point mascon model)，可進(jìn)行15 000顆空間碎片軌道的并行計(jì)算，其精度較SGP4/SDP4低，運(yùn)算速度快，加速比可達(dá)10～15倍. Lin等[3]基于CUDA的并行架構(gòu)，采用塊分解算法，對(duì)空間碎片進(jìn)行分批次處理，相鄰批次的計(jì)算和數(shù)據(jù)傳輸操作互相交疊.基于此算法，150個(gè)飛行器3 d內(nèi)的碰撞預(yù)報(bào)不超過3 h，在大型高性能GPU顯卡上計(jì)算速度可達(dá)CPU加速的30倍.郭松等[4]在FPGA平臺(tái)上設(shè)計(jì)并實(shí)現(xiàn)了基于VLIW的SGP4/SDP4軌道預(yù)報(bào)加速器，在單片F(xiàn)GPA上用8個(gè)處理單元組成計(jì)算陣列，采用指令級(jí)并行和任務(wù)級(jí)并行相結(jié)合的方式，將軌道預(yù)報(bào)速度提升到了CPU的8倍以上.馮森[5]利用CUDA并行架構(gòu)和GPU硬件作為CPU的協(xié)處理器，實(shí)現(xiàn)了SGP4軌道預(yù)報(bào)模型計(jì)算時(shí)間的大幅縮小，實(shí)驗(yàn)結(jié)果證明其加速比達(dá)到15倍.綜上所述，衛(wèi)星軌道預(yù)報(bào)并行遞推面臨解算速度以及解算精度這兩個(gè)相互沖突的問題，也即如何在保證解算速度的同時(shí)盡可能的提高軌道預(yù)測(cè)精度.目前精度較高的解算方法都為外推法，如四階龍格庫(kù)塔法以及HPOP高精度外推模型[6].但由于外推法計(jì)算量極大，且其外推的特點(diǎn)在保證精度的同時(shí)卻無(wú)法進(jìn)行并行計(jì)算，外推法難以應(yīng)用于快速解算之中. 因此，本文基于SGP4/SDP4模型，考慮攝動(dòng)模型的高精度性質(zhì)及解析特性，即同一衛(wèi)星在不同時(shí)間點(diǎn)的軌道模型解可獨(dú)立解得，符合并行計(jì)算的特點(diǎn)，提出GPU并行計(jì)算方法，使不同衛(wèi)星的星歷、同一顆星不同時(shí)間點(diǎn)的星歷在不同的線程中同步計(jì)算，以達(dá)到數(shù)倍于CPU處理器速度的加速效果.

Lin等[3]和馮森[5]都在高性能工作站上進(jìn)行并行解算實(shí)驗(yàn)，加速性能達(dá)到較高程度.本文擬在小型板卡NVIDIA TX2上完成并行軌道模型加速解算的部署與優(yōu)化，大大有利于推進(jìn)基于GPU的衛(wèi)星軌道并行計(jì)算的學(xué)術(shù)研究與工程應(yīng)用.

1 軌道預(yù)報(bào)模型

軌道預(yù)報(bào)是指已知衛(wèi)星某一時(shí)刻位置、速度狀態(tài)，根據(jù)衛(wèi)星的動(dòng)力學(xué)微分方程來(lái)預(yù)測(cè)衛(wèi)星未來(lái)一段時(shí)間內(nèi)的位置和速度.為了表示衛(wèi)星的位置、速度信息，北美防空司令部(The north amerian aerospace defense command，NORAD)開發(fā)了雙行元(two-line element，TLE)[7].其綜合考慮了大氣阻力、日月引力的長(zhǎng)期和周期攝動(dòng)、地球扁率等影響.TLE作為平均根數(shù)，使用特定的方法消除了周期擾動(dòng)項(xiàng)，解算時(shí)必須使用同樣的方法重構(gòu)，因此TLE只適用于特定的軌道預(yù)報(bào)模型.

針對(duì)不同軌道高度的空間目標(biāo)，NORAD開發(fā)了與TLE相應(yīng)的近地或深空模型，目前定期更新的TLE針對(duì)的模型是SGP4/SDP4模型. SGP4模型(Simplified general perturbation version 4)能夠?qū)壍栏叨容^低或繞地飛行周期小于225 min的航天器進(jìn)行軌道預(yù)報(bào)，SDP4模型(Simplified deep space perturbation version 4)能夠?qū)壍栏叨容^高或繞地飛行周期大于225 min的航天器進(jìn)行軌道預(yù)報(bào). SGP4/SDP4模型預(yù)報(bào)精度較高，模型為解析模型，計(jì)算復(fù)雜度較遞推模型低，綜合性能高，因此其在衛(wèi)星軌道預(yù)報(bào)領(lǐng)域被人們廣泛使用[8-9].

SGP4預(yù)報(bào)模型解算采用解析法，需要將輸入數(shù)據(jù)轉(zhuǎn)換為初始軌道位置.由于TLE中采用平均運(yùn)動(dòng)，解算時(shí)需要將其重新轉(zhuǎn)換為半長(zhǎng)軸，轉(zhuǎn)換方法與正常軌道根數(shù)轉(zhuǎn)換方法有差異.

首先，求解平均運(yùn)動(dòng)對(duì)應(yīng)半長(zhǎng)軸,即

(1)

式中：μ為地球引力常數(shù);n0為平均運(yùn)動(dòng).

J2攝動(dòng)下的軌道傾角和偏心率分別為：

(2)

(3)

由此，對(duì)半長(zhǎng)軸進(jìn)行修正：

(4)

(5)

(6)

式中：a′0、n′0分別為改進(jìn)的半長(zhǎng)軸和周期參數(shù)[10-11].

對(duì)于低軌道衛(wèi)星，即地心距小于98 km，SGP4功率函數(shù)中常變量s的表達(dá)式為

(7)

計(jì)算出如下參數(shù)：

θ=cosi0,

(8)

(9)

η=a″0e0ζ.

(10)

利用初始參數(shù)，進(jìn)行SGP4軌道預(yù)報(bào)模型的遞推求解.SGP4模型包含大氣阻力攝動(dòng)、非球形軌道攝動(dòng)、長(zhǎng)周期攝動(dòng)項(xiàng)和短周期攝動(dòng)項(xiàng)的計(jì)算[12].對(duì)于低軌道目標(biāo)，主要攝動(dòng)項(xiàng)為大氣阻力攝動(dòng)和非球形攝動(dòng)，表達(dá)式如下：

(11)

(12)

(13)

式中：M0為某段時(shí)間內(nèi)的平均平近角；Ω0為某段時(shí)間內(nèi)的平均升交點(diǎn)赤經(jīng).

大氣阻力攝動(dòng)和非球形攝動(dòng)造成軌道根數(shù)發(fā)生變化，有：

δω=B*C3(cosω0)(t-t0),

(14)

(1+ηcosM0)3],

(15)

MP=MDF+δω+δM,

(16)

ω=ωDF-δω-δM,

(17)

(18)

e=e0-B*C4(t-t0)-B*C5(sinMP-sinM0),

(19)

D3(t-t0)3-D4(t-t0)4]2,

(20)

(21)

(22)

(23)

式中，B*為SGP4典型阻力系數(shù).在近心點(diǎn)距小于220 km時(shí)，a和IL的表達(dá)式中C1后面的項(xiàng)可舍去，含有C5、δω、δM的項(xiàng)可忽略[13-15].

2 基于GPU的并行計(jì)算流程設(shè)計(jì)

2.1 GPU并行計(jì)算

GPU是一種在個(gè)人電腦、工作站、游戲機(jī)和一些移動(dòng)設(shè)備(如平板電腦、智能手機(jī)等)上專門進(jìn)行圖像運(yùn)算工作的微處理器.如圖1所示，GPU與CPU組成相同，都由運(yùn)算基本單元、控制器、寄存器組成. CPU內(nèi)大部分空間被控制器和寄存器占據(jù)，其可利用的臨時(shí)緩存較大，擅長(zhǎng)邏輯控制和串行計(jì)算，多用于處理包含不同數(shù)據(jù)類型的復(fù)雜操作.相比于CPU處理器，GPU擁有更多的ALU(算術(shù)邏輯單元，是能實(shí)現(xiàn)多組算術(shù)運(yùn)算和邏輯運(yùn)算的組合邏輯電路)用于數(shù)據(jù)處理，而非數(shù)據(jù)高速緩存和流控制，此類結(jié)構(gòu)適合于對(duì)密集型數(shù)據(jù)進(jìn)行并行處理，例如本文中所處理的多衛(wèi)星軌道數(shù)據(jù).

圖1 CPU與GPU結(jié)構(gòu)對(duì)比示意

因此，GPU可同時(shí)調(diào)用大量相同類型的數(shù)據(jù)，并將它們分配入不同線程中進(jìn)行相同類型的計(jì)算操作，其耗費(fèi)的時(shí)間和單獨(dú)計(jì)算一個(gè)數(shù)據(jù)的時(shí)間一致，這就是GPU并行計(jì)算的概念，即“以空間換時(shí)間”的處理方式.

2.2 CUDA計(jì)算平臺(tái)

CUDA(Compute unified device architecture)是由NVIDIA公司推出的一種通用并行計(jì)算框架，僅適用于NVIDIA公司的產(chǎn)品(以顯卡產(chǎn)品為主).其架構(gòu)包含軟硬件體系，硬件為CUDA支持的GPU處理器，軟件則包括顯卡驅(qū)動(dòng)程序、nvcc編譯器、調(diào)試模塊等組成的開發(fā)環(huán)境，支持C語(yǔ)言及C++語(yǔ)言編寫，通用性較強(qiáng)，易于上手.可將CUDA框架簡(jiǎn)單理解為一個(gè)可通過特定編程模型實(shí)現(xiàn)GPU并行加速的C語(yǔ)言函數(shù)庫(kù).本文中測(cè)試程序均于Visual Studio 2013與CUDA 9.0開發(fā)平臺(tái)完成.

CUDA框架結(jié)構(gòu)中，存在兩個(gè)必要的程序組成部分：主機(jī)程序以及內(nèi)核函數(shù)kernal.主機(jī)程序運(yùn)行在CPU處理器中，實(shí)現(xiàn)了CPU部分的計(jì)算、CPU與內(nèi)存間數(shù)據(jù)傳輸、CPU與GPU間內(nèi)存交互、設(shè)置GPU并行計(jì)算參數(shù)的功能.內(nèi)核函數(shù)僅運(yùn)行于GPU各線程中，為實(shí)際并行計(jì)算的運(yùn)算部分，各GPU線程內(nèi)的相同內(nèi)核函數(shù)帶入不同參數(shù)同步進(jìn)行計(jì)算.如圖2所示，CUDA框架內(nèi)GPU線程結(jié)構(gòu)分為3個(gè)層次：網(wǎng)格(Grid)、線程塊(Block)和線程(Thread).每當(dāng)一個(gè)內(nèi)核函數(shù)被發(fā)送到GPU時(shí)，即生成一個(gè)自定義網(wǎng)格，網(wǎng)格內(nèi)存在若干自定義維度的線程塊，每個(gè)線程塊內(nèi)存在同等數(shù)量的線程.每個(gè)線程都分配有相應(yīng)的寄存器和本地內(nèi)存，在其所在的線程塊內(nèi)有一共享內(nèi)存空間可供該線程塊內(nèi)所有線程讀取，同時(shí)對(duì)于整個(gè)網(wǎng)格有可供所有線程訪問讀取的全局常量及內(nèi)存，GPU并行計(jì)算可根據(jù)此原理進(jìn)行參數(shù)設(shè)置以達(dá)到最優(yōu)加速效果.

圖2 CUDA線程結(jié)構(gòu)示意

2.3 集成式加速流程設(shè)計(jì)

本文集成式框架的輸入項(xiàng)為一TXT格式的雙行元文件，其包含任意數(shù)目衛(wèi)星的某一時(shí)刻雙行元數(shù)據(jù)，表示衛(wèi)星初始的軌道數(shù)據(jù)，多顆衛(wèi)星的雙行元數(shù)據(jù)被統(tǒng)一讀取并存入內(nèi)存.

衛(wèi)星初始數(shù)據(jù)讀入后，從內(nèi)存中依次調(diào)出每一顆星的雙行元數(shù)據(jù)，并根據(jù)雙行元數(shù)據(jù)進(jìn)行軌道初始參數(shù)(軌道六根數(shù))的計(jì)算并存入內(nèi)存，由于此過程單顆衛(wèi)星參數(shù)計(jì)算時(shí)間小于SGP4模型解算時(shí)間的1/10，因此本文中并未對(duì)此計(jì)算流程進(jìn)行并行加速.此過程為在CPU運(yùn)算單元中按順序進(jìn)行每一顆星的軌道參數(shù)初始化操作.

如圖3所示，完成所有雙行元數(shù)據(jù)的軌道初始參數(shù)計(jì)算后，程序進(jìn)入第1次數(shù)據(jù)傳輸階段，即CPU至GPU的數(shù)據(jù)傳輸.根據(jù)雙行元數(shù)據(jù)的個(gè)數(shù)開辟相應(yīng)大小的GPU內(nèi)存空間，每個(gè)雙行元數(shù)據(jù)對(duì)應(yīng)的GPU內(nèi)存空間內(nèi)分別儲(chǔ)存以下4種數(shù)據(jù)：起始時(shí)間、結(jié)束時(shí)間、解算步長(zhǎng)以及衛(wèi)星星歷.其中衛(wèi)星星歷在程序中是以一個(gè)結(jié)構(gòu)體變量形式存儲(chǔ)的，它包含了每一顆星位置、速度與軌道參數(shù)信息.CPU到GPU的數(shù)據(jù)傳輸過程是將CPU內(nèi)存中每顆星起始時(shí)間、結(jié)束時(shí)間、解算步長(zhǎng)及起始時(shí)刻表示衛(wèi)星狀態(tài)的結(jié)構(gòu)體變量數(shù)據(jù)打包傳輸?shù)紾PU內(nèi)存空間中，等待GPU調(diào)用.

圖3 程序流程圖

第1次數(shù)據(jù)傳輸完成后，GPU將每一顆星對(duì)應(yīng)的結(jié)構(gòu)體變量帶入多個(gè)獨(dú)立線程中的核函數(shù)進(jìn)行并行計(jì)算.GPU內(nèi)每一個(gè)運(yùn)算單元為一獨(dú)立線程，每一獨(dú)立線程中計(jì)算一顆星在某一時(shí)間點(diǎn)的星歷，則每一顆星的結(jié)構(gòu)體變量需帶入計(jì)算的線程數(shù)量為

(24)

式中：n為單顆星所需線程數(shù)量；ts為結(jié)束時(shí)間；t0為起始時(shí)間；tp為步長(zhǎng)

在并行計(jì)算過程中，多顆星的計(jì)算步驟不是順序進(jìn)行的，而是同步進(jìn)行的，即多顆星在獨(dú)立線程中的解算是同時(shí)進(jìn)行的，以此達(dá)到加速效果.所需的并行線程數(shù)量為

N=ns×n.

(25)

式中:N為總線程數(shù)量；ns為衛(wèi)星數(shù)；n為單星線程數(shù)量.

以計(jì)算一顆衛(wèi)星24 h的星歷為例，以1 s為解算步長(zhǎng)，所需線程數(shù)為86 500個(gè)，計(jì)算200顆星即為17 300 000個(gè)線程，所需的運(yùn)算單元數(shù)超過CPU串行運(yùn)算2～3個(gè)數(shù)量級(jí)，符合本文 “以空間換時(shí)間”的理念.但因此，使用本方法進(jìn)行并行解算的星數(shù)會(huì)受GPU內(nèi)存空間大小的限制.例如，在Jetson TX2設(shè)備上可同步解算500顆星，但在PC電腦端本方法只可同步遞推50顆星的軌道.

GPU并行計(jì)算完成后得到衛(wèi)星星歷數(shù)據(jù)，進(jìn)行第2次數(shù)據(jù)傳輸(GPU至CPU),將星歷數(shù)據(jù)傳輸回CPU內(nèi)存，再輸出到文件中結(jié)束解算.

3 實(shí)驗(yàn)加速效果與分析

為了驗(yàn)證所提集成式加速方法的有效性和適應(yīng)性，即本文提出方法是否達(dá)到加速效果，若達(dá)到加速效果，在何種計(jì)算條件及硬件設(shè)備在能達(dá)到最好的使用效果，本文分別在兩種仿真環(huán)境下進(jìn)行GPU并行加速運(yùn)行試驗(yàn)：筆記本電腦PC端環(huán)境以及JETSON TX2開發(fā)板環(huán)境.

在筆記本電腦PC端進(jìn)行的實(shí)驗(yàn)對(duì)比了CPU解算速度與GPU集成式加速方法解算速度，驗(yàn)證了集成式加速方法的有效性.又將兩種環(huán)境下的集成式加速結(jié)果與模塊化加速結(jié)果進(jìn)行比較，得出兩種加速方法各自的適用場(chǎng)景.最后通過對(duì)比兩種環(huán)境下集成式加速方法的加速效果，討論了GPU并行加速方法的設(shè)備適用性，提供了衛(wèi)星在軌解算平臺(tái)的硬件選擇建議.

3.1 集成式加速方法有效性分析

本文在筆記本電腦PC端環(huán)境下首先進(jìn)行了CPU處理器上的SGP4/SDP4軌道模型解算實(shí)驗(yàn)，又于同一臺(tái)設(shè)備上的GPU處理器進(jìn)行集成式方法的并行加速實(shí)驗(yàn)，兩次實(shí)驗(yàn)結(jié)果如表2和圖4所示.

圖4 CPU和GPU解算速度對(duì)比

表2 仿真環(huán)境1下仿真結(jié)果

仿真過程中每顆星解算步長(zhǎng)為1 s，解算時(shí)長(zhǎng)為24 h，則每顆星計(jì)算86 400步，輸出86 400組解算結(jié)果.

試驗(yàn)驗(yàn)證環(huán)境1(筆記本電腦)：

操作系統(tǒng)：Windows 10中文版.

CPU：英特爾 core i7-6700HQ，2.6 Hz主頻，8 GB內(nèi)存.

GPU：NVIDIA GeForce GTX 950 M，包含640個(gè)CUDA單元，4 096 MB顯存.

集成式方法需較大的數(shù)據(jù)儲(chǔ)存空間進(jìn)行運(yùn)算，鑒于GPU顯卡顯存及運(yùn)算單元數(shù)量限制，筆記本電腦環(huán)境僅能并行解算50顆星星歷數(shù)據(jù)，而TX2開發(fā)板能并行解算500顆星的數(shù)據(jù).顯存空間的大小及GPU運(yùn)行單元數(shù)量決定并行計(jì)算所能計(jì)算的星數(shù)上限，并影響加速效果.

由實(shí)驗(yàn)結(jié)果如表2及圖4可知，在星數(shù)較少時(shí)，總運(yùn)算步數(shù)也較少，CPU解算時(shí)間要低于GPU并行解算時(shí)間.這是由于CPU處理器在進(jìn)行單次解算時(shí)解算速度為GPU處理器速度的數(shù)倍；且進(jìn)行GPU加速時(shí)，CPU調(diào)用GPU核函數(shù)進(jìn)入線程的時(shí)間為主要消耗時(shí)間，真正的GPU線程解算時(shí)間都僅為解算一個(gè)星歷時(shí)間點(diǎn)的解算時(shí)間.兩者的共同因素導(dǎo)致在總運(yùn)算步數(shù)較低時(shí)，采用CPU處理器運(yùn)算效率較高.

當(dāng)星數(shù)增加時(shí)，數(shù)據(jù)量增大，CPU調(diào)用核函數(shù)的時(shí)間隨之增加，但調(diào)用函數(shù)耗費(fèi)時(shí)間增長(zhǎng)率較無(wú)并行計(jì)算的CPU解算時(shí)間增長(zhǎng)率低，此時(shí)GPU并行解算開始呈現(xiàn)速度優(yōu)勢(shì)，加速比(即GPU解算速度除以CPU解算速度)最終逐漸穩(wěn)定于4.6，達(dá)到較高加速比.

3.2 集成式加速方法適應(yīng)性分析

上文分析了集成式方法的加速有效性，以下對(duì)集成式方法和模塊化方法進(jìn)行對(duì)比，提出兩種方法各自的適用范圍，以及加速方法的設(shè)備適用性.

本文于NVIDIA Jetson TX2小型開發(fā)板上進(jìn)行了集成式加速方法的SGP4/SDP4軌道遞推模型并行計(jì)算實(shí)驗(yàn)，將其結(jié)果與計(jì)天陽(yáng)[1]中模塊化加速方法實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，得出兩者方法的適用計(jì)算范圍.又通過集成式方法在試驗(yàn)環(huán)境1、2上加速效果的對(duì)比，分析提出軌道模型GPU并行加速的設(shè)備適用性.試驗(yàn)結(jié)果對(duì)比見圖5和表3.

表3 不同GPU解算方法運(yùn)算時(shí)間表

圖5 不同方法GPU解算速度對(duì)比

試驗(yàn)驗(yàn)證環(huán)境2(NVIDIA TX2小型開發(fā)板)：

操作系統(tǒng)：ARM-LINUX Ubuntu 16.04

GPU：NVIDIA Pascal，256 CUDA 核心

在進(jìn)行GPU并行解算過程中，較GPU的解算速度，CPU調(diào)用核函數(shù)的速度更為緩慢，如圖6所示為一顆星86 400步數(shù)下使用Visual Profiler分析程序得出的各進(jìn)程占用時(shí)間比圖，可見建立加速進(jìn)程的過程較核函數(shù)運(yùn)行時(shí)間長(zhǎng).為提高解算速度，本文提出減少調(diào)用GPU并行計(jì)算核函數(shù)次數(shù)以達(dá)到更高速度的集成式加速方法.

圖6 GPU加速后進(jìn)程時(shí)間線(左邊黃色為調(diào)用創(chuàng)建進(jìn)程的時(shí)間)

集成式解算方法將SGP4解算模型(包括初始常量定義)作為整體加入核函數(shù)計(jì)算過程中.計(jì)算機(jī)內(nèi)存僅需與GPU內(nèi)存進(jìn)行一次交互數(shù)據(jù)傳輸和一次核函數(shù)調(diào)用及顯存空間分配，相比模塊化方法省去多次參數(shù)交互和函數(shù)調(diào)用的時(shí)間.但此加速方法單次計(jì)算操作復(fù)雜，占用顯存空間大，當(dāng)并行計(jì)算的星數(shù)達(dá)到一定數(shù)量時(shí)，會(huì)由于緩存不足出現(xiàn)溢出錯(cuò)誤，造成解算失敗.

模塊化解算方法則將解算過程分為4個(gè)模塊：重力攝動(dòng)常數(shù)初始化模塊、雙行元與定軌參數(shù)轉(zhuǎn)換模塊、SGP4軌道模型初始化模塊、SGP4軌道預(yù)報(bào)模塊.除重力攝動(dòng)常數(shù)初始化模塊外，其余3個(gè)模塊在GPU中進(jìn)行并行計(jì)算.分模塊進(jìn)行解算對(duì)顯存使用需求小，可以實(shí)現(xiàn)大規(guī)模衛(wèi)星軌道并行計(jì)算.但此方法在CPU與GPU間參數(shù)交互的過程中耗費(fèi)大量時(shí)間，加速效果不夠理想.

由圖5及表3，綜合分析集成式和模塊化兩種并行計(jì)算方法的試驗(yàn)數(shù)據(jù)可發(fā)現(xiàn)，在運(yùn)算數(shù)較小(本文試驗(yàn)中運(yùn)算組數(shù)*步數(shù)<400萬(wàn)步)時(shí)，宜采用集成式解算方法，在運(yùn)算數(shù)較大(本文試驗(yàn)中運(yùn)算組數(shù)*步數(shù)>400萬(wàn)步)時(shí)，宜采用模塊化解算方法.兩種加速方法對(duì)比及適用范圍可見表4.

在得出加速方法最佳適用范圍后，還需討論加速方法所適用的計(jì)算設(shè)備，為搭建衛(wèi)星在軌解算平臺(tái)提供參考.

由實(shí)驗(yàn)結(jié)果圖5和表4可知，低計(jì)算量時(shí)相同計(jì)算量下集成式加速方法解算時(shí)間低于模塊化加速方法，集成式加速方法在TX2設(shè)備上實(shí)現(xiàn)解算速度最高.TX2開發(fā)板GPU計(jì)算性能優(yōu)越，CPU性能較普通筆記本計(jì)算機(jī)低，其搭載的專用ARM-LINUX操作系統(tǒng)內(nèi)核使得其上內(nèi)存與GPU內(nèi)存數(shù)據(jù)傳輸速率較普通計(jì)算高，因此在GPU并行運(yùn)算時(shí)間差距不大的情況下，TX2開發(fā)板能夠?qū)崿F(xiàn)最高的解算速度.

表4 兩種解算方法對(duì)比及適用范圍

由圖5可見，兩種解算方法解算時(shí)間與計(jì)算量都呈現(xiàn)出線性關(guān)系，解算時(shí)間隨計(jì)算量線性增長(zhǎng).集成式加速在筆記本上解算時(shí)間隨時(shí)間增長(zhǎng)速率最高，經(jīng)分析認(rèn)為是WINDOWS操作系統(tǒng)下Visual Studio開發(fā)環(huán)境下，CUDA程序中由CPU內(nèi)存至GPU內(nèi)存的數(shù)據(jù)傳輸過程需經(jīng)由數(shù)次內(nèi)存轉(zhuǎn)移到顯卡設(shè)備，導(dǎo)致解算速度較低，低于LINUX系統(tǒng)環(huán)境效率.此現(xiàn)象也說(shuō)明，在進(jìn)行在軌衛(wèi)星軌道遞推解算時(shí)，平臺(tái)應(yīng)盡量選用搭載結(jié)構(gòu)簡(jiǎn)約的LINUX系統(tǒng)的計(jì)算機(jī)，以降低數(shù)據(jù)傳輸帶來(lái)的時(shí)間耗費(fèi).

4 GPU運(yùn)算數(shù)據(jù)精度分析

獲取在試驗(yàn)環(huán)境1下CPU解算和并行解算得出的兩份星歷遞推數(shù)據(jù)，遞推時(shí)間為24 h，即86 400步.將數(shù)據(jù)導(dǎo)入MATLAB中對(duì)比，進(jìn)行二者衛(wèi)星星歷xyz軸上位置誤差及直線距離誤差計(jì)算，因本文實(shí)驗(yàn)中筆記本計(jì)算機(jī)上顯卡設(shè)備與TX2設(shè)備的GPU雙精度運(yùn)算位數(shù)相同，試驗(yàn)環(huán)境1與試驗(yàn)環(huán)境2并行計(jì)算解算得到的星歷數(shù)據(jù)一致.用于仿真的雙行元數(shù)據(jù)見表5，精度仿真結(jié)果如圖8所示.

表5 仿真實(shí)驗(yàn)中使用的雙行元數(shù)據(jù)

圖7 GPU運(yùn)算位置與CPU運(yùn)算對(duì)比

圖8 GPU計(jì)算速度與CPU計(jì)算對(duì)比

與CPU運(yùn)算相比，GPU運(yùn)算的位置誤差在10-6km量級(jí)，速度誤差在10-9km量級(jí)，誤差較小，精度較高，但誤差并不收斂，隨著運(yùn)行時(shí)間的增長(zhǎng)逐漸發(fā)散.誤差發(fā)散的主要原因在于CPU和GPU中的變量精度不同.本文對(duì)CPU變量及GPU核函數(shù)中變量都以雙精度形式保存.CPU中雙精度運(yùn)算位數(shù)為80位，GPU中雙精度運(yùn)算位數(shù)為64位，隨著解算過程的進(jìn)行，代入解析解的數(shù)值增大，GPU雙精度計(jì)算產(chǎn)生的誤差隨之增大，結(jié)果誤差也越大，但總體在可接受范圍.

5 結(jié) 論

1)該方法將完整的SGP4/SDP4軌道預(yù)報(bào)模型解算過程寫入GPU核函數(shù)中進(jìn)行運(yùn)算.在中小規(guī)模計(jì)算量下的加速效果優(yōu)于模塊化加速方法.

2)通過實(shí)驗(yàn)仿真對(duì)比，集成式加速方法較模塊化加速方法在中小規(guī)模計(jì)算量下加速比更高，且精度損失小.

3)由于集成式方法占用較大顯存容量，因此難以用于超大規(guī)模計(jì)算量的加速解算，此外，由于GPU內(nèi)部計(jì)算條件限制，導(dǎo)致誤差會(huì)隨著運(yùn)行時(shí)間增長(zhǎng)而逐漸發(fā)散.針對(duì)這兩個(gè)問題，在后續(xù)工作中將采用新的算法嘗試進(jìn)行優(yōu)化.