譚勖立 王慶賓 范 雕 馮進(jìn)凱 黃 炎 黃子炎

1 信息工程大學(xué)地理空間信息學(xué)院，鄭州市科學(xué)大道62號(hào)，450001

隨著衛(wèi)星重力測(cè)量技術(shù)的發(fā)展，利用重力衛(wèi)星觀(guān)測(cè)數(shù)據(jù)反演地球重力場(chǎng)已成為普遍的重力場(chǎng)研究手段[1]。針對(duì)重力衛(wèi)星觀(guān)測(cè)特點(diǎn)，眾多學(xué)者提出不同的重力場(chǎng)反演方法，總體上分為空域法和時(shí)域法[2]，時(shí)域法又可分為能量守恒方法[3-5]、短弧法[1,6-8]、動(dòng)力學(xué)法[9]等。時(shí)域法反演地球重力場(chǎng)時(shí)，在處理海量觀(guān)測(cè)數(shù)據(jù)、獲取觀(guān)測(cè)方程以及構(gòu)建并解算法方程的過(guò)程中存在計(jì)算量巨大、計(jì)算耗時(shí)長(zhǎng)的問(wèn)題[10]。針對(duì)以上問(wèn)題，許多學(xué)者基于并行計(jì)算平臺(tái)，分析反演方法的并行潛力，設(shè)計(jì)了并行算法，以提高計(jì)算效率。文獻(xiàn)[10]總結(jié)并行技術(shù)提高反演效率的關(guān)鍵任務(wù)，提出了OpenMP并行算法；文獻(xiàn)[11]實(shí)現(xiàn)了OpenMP并行快速解算重力場(chǎng)模型，反演得到60階模型，解算耗時(shí)減少75%；文獻(xiàn)[12-13]在曙光集群高性能計(jì)算系統(tǒng)上實(shí)現(xiàn)了MPI并行算法，解算120、240階模型時(shí)，法矩陣求逆耗時(shí)僅為229 s、7 359 s，解算總耗時(shí)分別為40 min、7 h；文獻(xiàn)[14]綜合OpenMP和MKL(Math Kernel Library)庫(kù)二者優(yōu)勢(shì)，有效提高了反演效率；文獻(xiàn)[15]分析基于超級(jí)計(jì)算機(jī)平臺(tái)的并行技術(shù)應(yīng)用于衛(wèi)星重力測(cè)量的相關(guān)問(wèn)題。

現(xiàn)今，反演重力場(chǎng)的并行算法研究主要針對(duì)OpenMP、MPI等CPU同構(gòu)并行展開(kāi)，存在并行算法硬件需求高、未充分挖掘反演方法并行潛力的問(wèn)題。隨著GPU(graphic processing unit)的發(fā)展，GPU單個(gè)核心運(yùn)算能力和芯片核心集成數(shù)量都得到大幅提升，CPU+GPU組成的異構(gòu)并行結(jié)構(gòu)成為常見(jiàn)的異構(gòu)并行運(yùn)算平臺(tái)，單臺(tái)個(gè)人電腦即可構(gòu)成其主體部分。相較于CPU，GPU具有多核心、高集成的特點(diǎn)，能夠同時(shí)開(kāi)辟上千條線(xiàn)程進(jìn)行運(yùn)算，從而實(shí)現(xiàn)密集運(yùn)算的高度并行化。

基于以上分析，本文針對(duì)衛(wèi)星跟蹤衛(wèi)星觀(guān)測(cè)模式，在能量守恒方法的基礎(chǔ)上，提出一種重力場(chǎng)反演快速異構(gòu)并行算法(后文簡(jiǎn)稱(chēng)快速反演算法)，利用CUDA在GPU端實(shí)現(xiàn)了并行計(jì)算觀(guān)測(cè)方程設(shè)計(jì)矩陣與自由項(xiàng)向量，結(jié)合MKL庫(kù)與分區(qū)平差法[16]、預(yù)處理共軛梯度法[17]在CPU端完成了低內(nèi)存消耗下的法方程的快速構(gòu)建和解算。用該算法處理GRACE-FO衛(wèi)星2020-01-01～06-30期間共13.77 GB觀(guān)測(cè)數(shù)據(jù)，反演獲得120階重力場(chǎng)模型GM-GraceFO2020h(GM為gravity model縮寫(xiě)，后綴GraceFO表示解算采用GRACE-FO數(shù)據(jù)，2020代表所用數(shù)據(jù)年份，h為heterogeneous parallelism的首字母)，并與現(xiàn)有GRACE模型和算法相對(duì)比。結(jié)果表明，快速反演算法所得模型與現(xiàn)有模型精度相當(dāng)，且相較于傳統(tǒng)串行算法，反演耗時(shí)減少98.479%，內(nèi)存消耗減小1個(gè)量級(jí)。算法在保證反演精度的同時(shí)，能有效提高計(jì)算效率，縮短反演時(shí)間，降低對(duì)計(jì)算平臺(tái)的硬件需求，可以為大規(guī)模處理衛(wèi)星重力數(shù)據(jù)提供經(jīng)濟(jì)、高效的實(shí)施方案。

1 能量守恒方法原理

能量守恒方法反演地球重力場(chǎng)，在慣性系下可得單星能量守恒積分式[5]：

(1)

Vconv=Vsolar-tide+Vlunar-tide+Vsolid-tide+

Vocean-tide+Vatmosphere-tide+Vpole-tide

(2)

雙星積分式由單星積分式作差獲得：

(3)

式中，下標(biāo)1、2為單星的編號(hào)，下標(biāo)12表示1、2星對(duì)應(yīng)項(xiàng)相減。將擾動(dòng)位T球諧展開(kāi)：

(4)

綜合式(1)～(4)，顧及球坐標(biāo)與直角坐標(biāo)的轉(zhuǎn)換關(guān)系，基于單星和雙星能量守恒方程可構(gòu)建統(tǒng)一的觀(guān)測(cè)方程：

(5)

2 算法設(shè)計(jì)

并行算法的設(shè)計(jì)主要圍繞兩個(gè)方面，分別為觀(guān)測(cè)方程的并行計(jì)算和法方程的并行計(jì)算，下面圍繞這兩點(diǎn)進(jìn)行討論。

2.1 觀(guān)測(cè)方程的并行計(jì)算

觀(guān)測(cè)方程的并行計(jì)算由CUDA實(shí)現(xiàn)。CUDA將計(jì)算任務(wù)劃分為如圖1所示格網(wǎng)，每個(gè)GPU線(xiàn)程完成對(duì)應(yīng)格網(wǎng)的任務(wù)，線(xiàn)程通過(guò)索引值index獲取對(duì)應(yīng)歷元的觀(guān)測(cè)數(shù)據(jù)，索引值index的計(jì)算方式如下:

圖1 任務(wù)格網(wǎng)Fig.1 Task-grid

index=blockIdx.x*blockDim.x+threadIdx.x

(6)

式中，blockDim為線(xiàn)程塊格網(wǎng)的維度，blockIdx為線(xiàn)程塊在格網(wǎng)中的編號(hào)，threadIdx為線(xiàn)程在線(xiàn)程塊中的編號(hào)。

GPU可開(kāi)辟與其核心數(shù)同樣多的線(xiàn)程。以GTX1080為例，其擁有2 560個(gè)核心，可開(kāi)辟2 560個(gè)線(xiàn)程同時(shí)進(jìn)行運(yùn)算，故GPU可充分利用能量守恒方法設(shè)計(jì)矩陣各行向量的獨(dú)立性，實(shí)現(xiàn)觀(guān)測(cè)方程構(gòu)建過(guò)程的高度并行化，在單位時(shí)間內(nèi)完成更多的計(jì)算。同時(shí)，為減少計(jì)算過(guò)程中的循環(huán)次數(shù)，參考文獻(xiàn)[19]，將中間變量向量化：

(7)

進(jìn)而設(shè)計(jì)矩陣行向量Ai及自由項(xiàng)元素li，皆可由向量運(yùn)算獲得(其中符號(hào)“·”表示向量?jī)?nèi)積，符號(hào)“°”表示同維向量對(duì)應(yīng)元素相乘)。

(8)

2.2 法方程的并行計(jì)算

MKL(math kernel library)庫(kù)是由Intel提供的函數(shù)庫(kù)，經(jīng)過(guò)了高度優(yōu)化，并實(shí)現(xiàn)了CPU上的并行計(jì)算，能提供包括線(xiàn)性代數(shù)在內(nèi)的各類(lèi)數(shù)學(xué)運(yùn)算，具有較高的性能和穩(wěn)定性。利用MKL庫(kù)中的cblas_dgemm函數(shù)可實(shí)現(xiàn)大型矩陣的乘法快速并行運(yùn)算。

在構(gòu)建法方程時(shí)，傳統(tǒng)算法需要將設(shè)計(jì)矩陣整體記錄在內(nèi)存，當(dāng)面對(duì)海量觀(guān)測(cè)數(shù)據(jù)時(shí)，設(shè)計(jì)矩陣將占用巨大的內(nèi)存空間。以衛(wèi)星182 d數(shù)據(jù)、60 s采樣間隔反演120階模型為例，此時(shí)將獲得262 080條數(shù)據(jù)記錄，設(shè)計(jì)矩陣需28.12 G內(nèi)存空間進(jìn)行存儲(chǔ)，個(gè)人電腦難以提供足夠大小的內(nèi)存。針對(duì)大型設(shè)計(jì)矩陣在存儲(chǔ)時(shí)存在內(nèi)存需求大的問(wèn)題，可采用分區(qū)平差的方法，對(duì)觀(guān)測(cè)數(shù)據(jù)按時(shí)序進(jìn)行分區(qū)處理。將觀(guān)測(cè)數(shù)據(jù)分區(qū)后，可得各分區(qū)的誤差方程：

(9)

式中，Ai、Bi和Li分別為i分區(qū)的設(shè)計(jì)矩陣和自由項(xiàng)向量。各分區(qū)分別構(gòu)建法方程，消去區(qū)域性未知數(shù)Yi，得到僅含聯(lián)系未知數(shù)的約化法方程：

(10)

將各聯(lián)系未知數(shù)的約化法方程系數(shù)和自由項(xiàng)相加，獲得總體約化法方程的法矩陣與自由項(xiàng)：

(11)

解算總體約化法方程便能求得聯(lián)系未知數(shù)X，若需進(jìn)一步求取區(qū)域性未知數(shù)，則將X回代至各分區(qū)法方程中即可。利用分區(qū)平差法，每次計(jì)算僅需原有內(nèi)存需求的ni/n(n為總觀(guān)測(cè)數(shù)，ni為各分區(qū)觀(guān)測(cè)數(shù))。與MKL庫(kù)結(jié)合后，在減小對(duì)硬件資源消耗的同時(shí)能夠快速完成法矩陣及自由項(xiàng)計(jì)算。

在解算法方程時(shí)，采用預(yù)處理共軛梯度法[17]，并引入MKL庫(kù)并行加速相關(guān)線(xiàn)性代數(shù)運(yùn)算，可以兼顧解算精度和計(jì)算效率。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)平臺(tái)與數(shù)據(jù)準(zhǔn)備

實(shí)驗(yàn)計(jì)算平臺(tái)為戴爾Precision 7530移動(dòng)工作站，CPU為Intel Xeon(頻率為2.9 GHz，6核12線(xiàn)程)，內(nèi)存容量32 GB，GPU為NVIDIA Quadro P3200(6 G顯存，1 792個(gè)CUDA核心)。程序編譯環(huán)境為Visual Studio 2015-Visual C++，CUDA10.0，OpenMP2.0并行庫(kù)與MKL19.0科學(xué)計(jì)算庫(kù)。實(shí)驗(yàn)所涉及算法都在Debug模式下進(jìn)行測(cè)試。

本文數(shù)據(jù)來(lái)源如下：1)由德國(guó)GFZ(German Research Centre for Geosciences) ISDC(Information System and Data Center for Geoscientific Data)數(shù)據(jù)中心提供的GRACE-FO衛(wèi)星1 Hz精密軌道數(shù)據(jù)GNV1B、加速度計(jì)數(shù)據(jù)ACT1B、衛(wèi)星姿態(tài)數(shù)據(jù)SCA1B以及0.5 Hz星間激光干涉測(cè)距數(shù)據(jù)LRI1B、0.2 Hz Ka/K波段測(cè)距數(shù)據(jù)KBR1B(測(cè)距數(shù)據(jù)主要用于雙星動(dòng)能差的計(jì)算[6,17])，各類(lèi)數(shù)據(jù)都以每天一個(gè)文件的形式保存；2)海洋潮汐數(shù)據(jù)來(lái)自于EOT11a模型；3)日月星歷數(shù)據(jù)由DE430模型計(jì)算獲得。

對(duì)于衛(wèi)星數(shù)據(jù)，首先利用觀(guān)測(cè)數(shù)據(jù)標(biāo)準(zhǔn)差依據(jù)3σ準(zhǔn)則進(jìn)行粗差剔除；然后用5階牛頓插值填補(bǔ)跨度小于30 s的缺失歷元數(shù)據(jù)，對(duì)于無(wú)法內(nèi)插的歷元?jiǎng)t用0補(bǔ)足；同時(shí)舍棄缺失歷元數(shù)過(guò)半的數(shù)據(jù)，最終獲得166 d的可用數(shù)據(jù)，每天有86 400個(gè)觀(guān)測(cè)歷元，總數(shù)據(jù)量為13.77 GB。

3.2 反演精度分析

利用本文算法處理GRACE-FO衛(wèi)星2020-01-01～06-30共166 d的可用數(shù)據(jù)，反演獲得120階重力場(chǎng)模型GM-GraceFO2020h。為檢驗(yàn)反演精度，將ITU_GRACE16、HUST-Grace2016s、Tongji-Grace02s、GGM05S等GRACE重力場(chǎng)模型作為對(duì)比組，與GM-GraceFO2020h進(jìn)行對(duì)比分析。同時(shí)選取2 190階模型EIGEN-6C4(該模型綜合了LAGEOS-1/2、GRACE、GOCE以及地形數(shù)據(jù)，具有較高的精度[20])作為基準(zhǔn)，對(duì)比檢驗(yàn)?zāi)Ｐ蛢?nèi)符合精度。

用GM-GraceFO2020h與4個(gè)對(duì)比組模型(截?cái)嘀?20階)分別計(jì)算全球1°×1°高程異常與重力異常格網(wǎng)，并求取各自相較于EIGEN-6C4的差值，統(tǒng)計(jì)結(jié)果如表1所示。

表1 各重力場(chǎng)模型與EIGEN-6C4模型間高程異常、重力異常差值

由表1可知，GM-GraceFO2020h與EIGEN-6C4的模型高程異常、重力異常差值相較于對(duì)比組，在數(shù)值上相近，差值最小值的差異分別在7 cm和1.2 mGal以?xún)?nèi)，最大值的差異分別在15 cm和11 mGal以?xún)?nèi)，平均值和標(biāo)準(zhǔn)差的差異分別在0.1 cm、0.001 mGal和0.2 cm、0.003 mGal以?xún)?nèi)。圖2為GM-GraceFO2020h與對(duì)比組模型相較于EIGEN-6C4的階方差。由圖2可見(jiàn)，GM-GraceFO2020h與GRACE重力場(chǎng)模型在各階位系數(shù)上的差異較小，總體上精度相近。圖3為GM- GraceFO2020h階方差曲線(xiàn)，其與Kaula曲線(xiàn)基本吻合，說(shuō)明該模型與Kaula準(zhǔn)則符合度較高，具有較高的內(nèi)符合精度。圖4為該模型的重力異常、緯向垂線(xiàn)偏差和高程異常。

圖2 各模型與EIGEN-6C4位系數(shù)差的階方差Fig.2 Difference degree with EIGEN-6C4

圖3 GM-GraceFO2020h模型階方差Fig.3 Degree variance of GM-GraceFO2020h

圖4 GM-GraceFO2020h模型重力異常、 緯向垂線(xiàn)偏差、高程異常Fig.4 Gravity anomaly, deflection of the vertical and height anomaly derived from GM-GraceFO2020h

3.3 計(jì)算效率分析

從觀(guān)測(cè)方程構(gòu)建、法矩陣計(jì)算、法方程解算3個(gè)方面對(duì)快速反演算法的計(jì)算效率進(jìn)行分析。同時(shí)，引入并行加速比Sn[21]作為量化加速效果的參數(shù)：

(12)

式中，n表示參與運(yùn)算的線(xiàn)程數(shù)，tn表示有n個(gè)線(xiàn)程參與時(shí)完成運(yùn)算所需的時(shí)間。

3.3.1 觀(guān)測(cè)方程構(gòu)建

在對(duì)數(shù)據(jù)進(jìn)行60 s間隔采樣條件下，分別利用串行算法與快速反演算法完成如下計(jì)算任務(wù)：1)處理衛(wèi)星23 d數(shù)據(jù)(共29 011條記錄)，構(gòu)建60階模型3 716個(gè)待求參數(shù)的觀(guān)測(cè)方程；2)處理75 d數(shù)據(jù)(共100 960條記錄)，構(gòu)建96階模型9 404個(gè)待求參數(shù)的觀(guān)測(cè)方程；3)處理166 d數(shù)據(jù)(共225 341條記錄)，構(gòu)建120階模型14 636個(gè)待求參數(shù)的觀(guān)測(cè)方程。表2為各任務(wù)耗時(shí)情況。

表2 不同數(shù)據(jù)規(guī)模的計(jì)算任務(wù)耗時(shí)對(duì)比

對(duì)比結(jié)果顯示，傳統(tǒng)串行算法耗時(shí)隨著計(jì)算量的增大呈線(xiàn)性遞增趨勢(shì)；而本文提出的算法在構(gòu)建觀(guān)測(cè)方法時(shí)，相較于串行算法，能有效提高計(jì)算效率，其并行加速比達(dá)到20以上，但隨著計(jì)算量的增加，加速比會(huì)逐漸減小到10左右。

3.3.2 法矩陣計(jì)算

分別用串行算法、OpenMP算法[11]、快速反演算法計(jì)算956階(設(shè)計(jì)矩陣維度為8 862×956)、3 716階(設(shè)計(jì)矩陣維度為29 011×3 716)、9 404階(設(shè)計(jì)矩陣維度為100 960×9 404)、14 636階(設(shè)計(jì)矩陣維度為225 341×14 636)法矩陣，時(shí)間消耗情況如表3，內(nèi)存空間消耗情況如圖5。

表3 不同階法矩陣計(jì)算耗時(shí)情況

圖5 內(nèi)存空間消耗情況Fig.5 Memory consumption

結(jié)合表3和圖5可見(jiàn)，快速反演算法完成各計(jì)算任務(wù)所需時(shí)間分別是串行算法的7.660%、2.753%、1.763%和1.537%，并行加速比在13～66之間，是OpenMP算法的52.078%、17.756%、11.869%和10.186%，其計(jì)算效率優(yōu)于二者。本文算法的內(nèi)存消耗是串行算法的1/20，較之小一個(gè)量級(jí)，也是OpenMP算法的1/2；同時(shí)，串行算法的內(nèi)存消耗會(huì)隨著法矩陣和設(shè)計(jì)矩陣維度的增加而快速增大，而本文算法的增速較小。通過(guò)以上分析可得，快速反演算法綜合MKL庫(kù)與分區(qū)平差法二者優(yōu)勢(shì)，相較于現(xiàn)有算法，在快速計(jì)算法矩陣的同時(shí)有效限制了對(duì)內(nèi)存資源的消耗，大幅降低了模型反演對(duì)計(jì)算平臺(tái)硬件水平的要求，使得在內(nèi)存資源有限的普通PC上處理大規(guī)模衛(wèi)星重力數(shù)據(jù)成為可能。

3.3.3 法方程解算

分別利用串行算法(Gauss-Jordan法)、MKL算法[14](LU分解求逆法)、快速反演算法，求解不同維數(shù)的法方程，計(jì)算耗時(shí)情況如表4。

表4 不同維數(shù)法方程解算耗時(shí)對(duì)比

由表4可見(jiàn)，快速反演算法在法矩陣階數(shù)為116時(shí)的解算耗時(shí)略多于MKL算法0.004 s，其主要原因在于，計(jì)算量較小時(shí)內(nèi)存讀寫(xiě)耗時(shí)具有較大占比；而在其余情況下其耗時(shí)皆低于前兩種算法，說(shuō)明該算法將MKL庫(kù)與預(yù)處理共軛梯度法相結(jié)合，能有效加快法方程的解算速度。但需要注意，共軛梯度法并不對(duì)法矩陣求逆，故需估計(jì)所求位系數(shù)的精度時(shí)，仍推薦使用MKL算法。

3.3.4 反演計(jì)算總體情況

分別利用串行算法和快速反演算法處理衛(wèi)星166 d數(shù)據(jù)，反演120階重力場(chǎng)模型，計(jì)算耗時(shí)總體情況如表5所示。

表5 反演計(jì)算總體耗時(shí)

由表5可見(jiàn)，串行算法反演耗時(shí)超過(guò)19 h，而快速反演算法僅需約17 min，反演耗時(shí)縮減98.479%，并行加速比超過(guò)65。算法耗時(shí)最多的計(jì)算過(guò)程為法方程的構(gòu)建，可采用僅計(jì)算法矩陣下三角部分、引入斯特拉森(Strassen)算法等方式來(lái)進(jìn)一步優(yōu)化算法，提高計(jì)算效率。

4 結(jié) 語(yǔ)

本文基于能量守恒方法，提出重力場(chǎng)反演快速異構(gòu)并行算法。該算法實(shí)現(xiàn)了GPU并行計(jì)算設(shè)計(jì)矩陣與自由項(xiàng)向量，通過(guò)結(jié)合MKL庫(kù)與分區(qū)平差、預(yù)處理共軛梯度法來(lái)快速構(gòu)建并解算法方程。利用本文算法處理GRACE-FO衛(wèi)星2020-01-01～06-30期間數(shù)據(jù)，反演獲得120階重力場(chǎng)模型GM-GraceFO2020h。通過(guò)對(duì)比實(shí)驗(yàn)，分析驗(yàn)證了算法的反演精度與計(jì)算效率。

在反演精度方面，快速反演算法所得模型與現(xiàn)有GRACE重力場(chǎng)模型在前120階精度相當(dāng)，具有較高的內(nèi)符合精度。在計(jì)算效率方面，該算法在重力場(chǎng)反演各主要計(jì)算過(guò)程的計(jì)算效率都優(yōu)于現(xiàn)有算法，相較于傳統(tǒng)串行算法，反演耗時(shí)縮減98.479%，內(nèi)存消耗減小1個(gè)量級(jí)。綜上所述，本文提出的重力場(chǎng)反演快速異構(gòu)并行算法，能在保證反演精度的同時(shí)，有效提高計(jì)算效率，減少反演耗時(shí)，降低硬件需求，可以為大規(guī)模處理衛(wèi)星重力數(shù)據(jù)提供經(jīng)濟(jì)高效的計(jì)算實(shí)施方案。

致謝：感謝德國(guó)GFZ ISDC數(shù)據(jù)中心提供GRACE-FO衛(wèi)星軌道、姿態(tài)、激光干涉測(cè)距及加速度計(jì)數(shù)據(jù)。