異構并行算法快速構建全球擾動重力梯度全張量圖

譚勖立，王慶賓，馮進凱，黃 炎，黃子炎

信息工程大學地理空間信息學院，鄭州 450001

0 引言

擾動重力位是地球重力位扣除正常橢球影響后的剩余部分，擾動重力梯度是其二階導數，能反映變化的不規則地球產生的高頻信息，可應用于大地測量[1-2]、地球物理學研究[3]、礦產資源探測[4-5]、海洋科學[6-7]和國防建設[8-9]等方面。獲取擾動重力梯度的方法主要有兩種[10]：一是利用高精度的重力梯度儀直接測量獲得；二是根據其與其他重力場元素間的函數關系計算獲得。由于西方更早開展儀器的相關研究，且針對該技術長期對國內進行封鎖，國內的研制進度相對落后，暫時難以開展實測[11]，故需依靠其他重力數據資料計算獲取重力梯度。可利用的資料有重力場模型[12-13]、重力異常數據[14-15]和地形數據[16]等。

在通過高階次位系數模型計算獲取大范圍高分辨率的擾動重力梯度數據時，存在重復運算多、耗時長的問題。曲政豪等[12]通過交換積分次序來提高計算速度；吳星等[13]引入非正常化球諧函數，推導了廣義球諧函數定積分，實現了擾動重力梯度的快速穩定計算。以上研究都從原理公式出發進行改進，但隨著并行計算技術及其平臺的發展，特別是以GPU+CPU為代表的異構并行計算平臺的出現，使得在程序算法上的改進具有更大的可能性與潛力。例如，黃炎等[17]基于CUDA(compute unified device architecture)設計了擾動引力矢量并行計算算法算，達到了擾動引力快速計算的目的。

本文推導了擾動重力梯度張量的簡化計算公式，將部分中間變量提取出來作為全局參數和局部參數，全局參數在整個計算過程中僅計算一次，局部參數也僅需在每個計算點計算一次，有效減少重復運算。在此基礎上提出擾動重力梯度快速異構并行算法，利用CUDA實現梯度張量在GPU端的并行計算。根據Txx、Tyy、Tzz三個分量滿足Laplace條件檢驗了算法可靠性，并與傳統串行算法進行計算效率對比。最后，利用該算法，基于2 190階地球重力場位系數模型EIGEN-6C4(截斷至2 160階)快速構建全球5′×5′分辨率的擾動重力梯度全張量圖，并對擾動重力梯度同部分地球物理屬性間的聯系及其在全球范圍內的數值特征進行分析。

1 基本原理

地球重力位可視為正常重力位與擾動重力位之和：

W=U+T。

(1)

式中：W為地球重力位；U為正常重力位，可由選定的正常橢球參數經過少量運算求得；T為擾動重力位，是位置的函數。在局部指北坐標系下求導可得

(2)

式中，函數下標x、y、z表示函數對對應坐標軸方向的偏導。矩陣中各元素為擾動重力梯度在各個方向上的分量。對于位函數，交換求導次序導數值不變，故共有6個不同分量，這6個分量構成擾動重力梯度全張量。同時，由于式(2)中矩陣對角線元素滿足Laplace條件，3個元素之和為0，因此在全張量中僅有5個獨立的量。

1.1 擾動重力梯度計算

擾動重力位的球諧展開式如下：

(3)

(4)

其中：

(5)

式中,各函數的具體形式可參考文獻[18]中的式(21)—式(25)，由于篇幅有限本文不予展示。

1.2 計算公式簡化

利用1.1節中方法大范圍計算擾動重力梯度全張量時存在大量的重復計算，需對其計算方法進行必要的改進，故本文推導了計算擾動重力梯度全張量的簡化計算公式。

文獻[18]中的式(21)—式(25)中含有大量僅與n、m有關的項，將這些項提取出來，構成全局系數anm、bnm、cn、dnm、enm：

(6)

在計算一次后存儲在內存中，之后通過n、m索引對應的數值，借此大幅減少與坐標無關項的重復計算。

在計算梯度全張量時，可將式(5)中與坐標相關但重復計算了多次的項作為局部變量，設

(7)

局部變量在每個待估點的計算過程中僅計算一次。盡可能地提取出不同項之間的公共被乘數，并將各自不相同的乘數部分相加，以進一步減少乘法運算的次數。最終，當m=0時:

(8)

當m=1時:

(9)

當m=2時:

(10)

當m≥3時:

(11)

將式(6)—(11)代入式(4)中，即可獲得擾動重力梯度全張量的簡化計算公式。在實際編程計算時，采用交換積分次序的方法，將式(3)中三角函數的值存為中間變量，同一緯圈僅計算一次與緯度相關的量，同一經圈僅計算一次與經度相關的量。

2 異構并行算法設計

異構并行算法基于1.2節中的計算方法進行設計，利用CUDA實現GPU端的并行計算。CUDA是由英偉達(NVIDIA)推出的通用并行計算架構，可將用C、C++、FORTRAN等語言編寫的程序在GPU端高性能并行執行。在設計并行算法時，需要考慮計算任務的劃分粒度和運算量之間的平衡。若粒度過小，則單個任務的運算量較小，不利于發揮GPU核心的計算能力，并且會浪費過多的算力資源在線程的控制與調配上；若粒度過大，則單個任務的運算量較大，GPU核心的計算能力難以達到負載要求，同時會致使計算過程的并行化程度不高，浪費了線程資源。因此，本文算法綜合考慮GPU核心的計算能力和擾動重力梯度計算過程的特點，采用分治的思想對計算任務進行遞歸劃分[19](圖1)。

第一層為待計算的擾動重力梯度格網；第二層為同緯度格網點的擾動重力梯度計算；第三層為單個格網點的擾動重力梯度計算。

在計算單個格網點的擾動重力梯度時，若計算的階次較高，則單個GPU核心難以一次性完成計算。因此，進一步劃分出更小的計算任務單元，也就是對函數fn(ρ,φ,λ)進行計算：

(12)

每個GPU僅計算其對應的fn(ρ,φ,λ)，最后將各GPU的計算結果求和，即可求得單個格網點的擾動重力梯度，即

(13)

式中，N為n的最大值，也稱為最大階數。眾多GPU同時參與單個格網點的擾動重力梯度計算，可大幅減少計算耗時，進而減少整個擾動重力梯度格網的計算耗時。

3 算法可靠性與效率檢驗

本文基于EIGEN-6C4模型進行了擾動重力梯度的計算，計算時將模型截斷至前2 160階，利用擾動引力三分量Txx、Tyy、Tzz滿足如式(14)所示Laplace條件的特性檢核計算結果，并對比本文算法與傳統串行算法的計算效率。

Txx+Tyy+Tzz=0。

(14)

EIGEN-6C4綜合了來自LAGEOS(laser geodynamics satellite)、GRACE(gravity recovery and climate experiment)、GOCE(gravity field and steady-state ocean circulation explorer)衛星和DTU12、EGM2008、EIGEN-6C3stat模型的多源數據，最高階次到達2 190，具有較高的精度[20]。實驗計算平臺為戴爾Precision 7530移動工作站，CPU為頻率2.9 GHz的Intel Xeon，內存容量為32 GB，GPU為NVIDIA Quadro P3200，擁有1 792個CUDA核心。編譯環境為Visual Studio 2015-Visual C++，CUDA版本10.0。

3.1 可靠性檢驗

利用本文算法計算了40°×60°區域內5′×5′分辨率的擾動重力梯度全張量。根據Laplace條件，即Txx、Tyy、Tzz三個分量之和為0，檢核計算結果，通過求和結果相對于0值的偏差情況評判計算結果是否可靠。計算獲得的三分量求和結果見表1。

由表1可以看出，擾動重力梯度三分量之和在0值附近波動，其最大值為7.816×10-13E，最小值為-9.948×10-13E，均方根為3.454×10-14E，大部分數值為10-15～10-14E。總體上三分量之和在10-12E的精度下接近于0，可認為本文算法所計算的擾動重力梯度Txx、Tyy、Tzz三個分量滿足Laplace條件，算法具有較高的可靠性與準確性。

表1 三分量檢核值

3.2 計算效率檢驗

為量化并行加速效果，引入并行加速比Sn作為參考量，其計算方式為[21]

(15)

式中：t1為串行計算所消耗的時間；tn為n個線程參與計算時消耗的時間。

分別利用本文算法與傳統串行算法在Debug模式下完成如表2所示的計算任務。耗時對比結果如表3所示。由表3可見，本文算法完成各個計算任務的耗時都遠低于串行算法，耗時分別縮短了91.425%、98.187%、98.444%和97.325%，并行加速比最高可達64.281。

表2 計算任務

表3 算法完成任務耗時

對比本文算法完成各任務時的并行加速比可知：1)任務1由于計算量較小，此時線程分配與調度的時間占比較大，故加速效果不明顯；2)由于本文采用的簡化算法旨在減少大量計算梯度全張量時的重復運算，而任務4僅計算了垂向分量，因此加速比低于任務2、3；3)本文算法在采用高階模型大量計算擾動重力梯度全張量時的加速效果最佳，且計算量越大加速效果越明顯。綜合來看，本文算法計算效率可達到傳統串行算法的60倍以上,并且計算量越大越能凸顯本文算法的優勢。

4 全球擾動重力梯度全張量圖構建

利用本文算法基于EIGEN-6C4模型構建了全球擾動重力梯度5′×5′分辨率的全張量圖(圖2)，計算了9 331 200個格網中心點，僅耗時3 h 43 min 56.326 s；而由表3推算傳統串行算法需耗費時間近2月。

為便于展示和分析梯度全張量數值空間分布特征，在圖2中將色標標值上下限固定為±100 E。由圖2可以看出：1)梯度全張量與地形具有較高的相關性，較為清晰地反映出千島海溝、馬納利亞海溝、四川盆地、青藏高原、柴達木盆地、埃塞俄比亞高原、東非大裂谷、大西洋中脊和安第斯山脈(圖2中藍線框住的區域)等構造特征，特別是Tzz分量，除陸地地形以外，其反映出的海底地形與文獻[22]中的模型具有較高的符合度；2)環太平洋區域、喜馬拉雅山脈以及東非大裂谷(圖2中紅線框住的區域)的梯度全張量都表現出較大的數值變化，表明這些區域質量分布變化較大，地質活動較為活躍，與全球主要的地震帶與火山帶的分布情況切合；3)在兩極部分地區，如格林蘭島東南部海岸(圖2中黑線框住的區域)，垂向分量明顯大于其他分量，可認為其變化可能主要由垂向的質量損失，即冰川消融引起的，在一定程度上佐證了文獻[23-25]中的發現。綜合以上分析，高精度高分辨率的全球擾動重力梯度張量圖在海底地形、地質變化和兩極冰川等方面的研究具有一定參考意義，有助于在大尺度下把握變化規律和信息。

圖2 全球擾動重力梯度全張量

表4反映了全球擾動重力梯度全張量的數值特征：Tzz分量數值變化范圍最大，達到-557.843～785.223 E；Txy分量數值變化范圍最小，僅為-138.612～149.166 E；各分量的平均值都接近于0 E，特別是Txy和Tyz分量，其均值在10-16量級上接近于0，表明擾動重力梯度的數值在全球范圍內正負分布較為均勻；Tyy、Tyz、Tzz分量數值變化較為劇烈，能夠較為清晰地反映擾動重力梯度的局部特征，而Txy分量數值變化較為平緩，難以反映出局部特征。

表4 全球擾動重力梯度全張量數值統計

5 結論

本文推導了擾動重力梯度張量的簡化計算公式，提出了擾動重力梯度快速異構并行算法，實現了梯度全張量在GPU端的并行計算。根據Txx、Tyy、Tzz三個分量滿足Laplace條件檢驗了算法可靠性，并與傳統串行算法對比了計算效率。最后利用本文算法，基于EIGEN-6C4模型(截至2 160階)快速構建了5′×5′分辨率的全球擾動重力梯度全張量圖。得出以下結論：

1)利用本文算法計算的擾動引力三分量之和在0值附近的波動極小，振幅小于10-12，擾動引力三分量在10-12量級上滿足Laplace條件，可認為計算結果準確可靠。

2)與串行算法對比，本文算法完成實驗中各計算任務的耗時分別縮短了91.425%、98.187%、98.444%和97.325%，并行加速比最高達到64.281，總體上可將計算效率提高60倍以上，算法快速高效。

3)全球擾動重力梯度全張量的計算結果顯示梯度全張量與地形具有較高的相關性，且能在一定程度上反映出局部區域質量分布變化情況。在數值分布方面，不同分量表現出不同的特征，Tzz分量數值變化范圍最大，Txy分量則最小；Tyy、Tyz、Tzz分量數值變化較為劇烈，Txy分量數值變化則較為平緩；各分量的平均值都接近于0 E。

綜合以上分析，本文算法能有效改善高階模型計算擾動重力梯度張量時效率低下的情況，可為獲取全球高分辨率擾動重力梯度數據提供高效可靠的實施方案；利用本文算法快速構建的全球擾動重力梯度全張量圖可以了解掌握各分量在全球范圍內的數值分布情況，對海洋地形、地質變化和兩極冰川相關研究也具有一定的參考意義。