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

信息工程大學(xué)地理空間信息學(xué)院，鄭州 450001

0 引言

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

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

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

1 基本原理

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

W=U+T。

(1)

式中：W為地球重力位；U為正常重力位，可由選定的正常橢球參數(shù)經(jīng)過少量運(yùn)算求得；T為擾動重力位，是位置的函數(shù)。在局部指北坐標(biāo)系下求導(dǎo)可得

(2)

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

1.1 擾動重力梯度計算

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

(3)

(4)

其中：

(5)

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

1.2 計算公式簡化

利用1.1節(jié)中方法大范圍計算擾動重力梯度全張量時存在大量的重復(fù)計算，需對其計算方法進(jìn)行必要的改進(jìn)，故本文推導(dǎo)了計算擾動重力梯度全張量的簡化計算公式。

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

(6)

在計算一次后存儲在內(nèi)存中，之后通過n、m索引對應(yīng)的數(shù)值，借此大幅減少與坐標(biāo)無關(guān)項的重復(fù)計算。

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

(7)

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

(8)

當(dāng)m=1時:

(9)

當(dāng)m=2時:

(10)

當(dāng)m≥3時:

(11)

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

2 異構(gòu)并行算法設(shè)計

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

第一層為待計算的擾動重力梯度格網(wǎng)；第二層為同緯度格網(wǎng)點(diǎn)的擾動重力梯度計算；第三層為單個格網(wǎng)點(diǎn)的擾動重力梯度計算。

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

(12)

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

(13)

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

3 算法可靠性與效率檢驗(yàn)

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

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)衛(wèi)星和DTU12、EGM2008、EIGEN-6C3stat模型的多源數(shù)據(jù)，最高階次到達(dá)2 190，具有較高的精度[20]。實(shí)驗(yàn)計算平臺為戴爾Precision 7530移動工作站，CPU為頻率2.9 GHz的Intel Xeon，內(nèi)存容量為32 GB，GPU為NVIDIA Quadro P3200，擁有1 792個CUDA核心。編譯環(huán)境為Visual Studio 2015-Visual C++，CUDA版本10.0。

3.1 可靠性檢驗(yàn)

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

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

表1 三分量檢核值

3.2 計算效率檢驗(yàn)

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

(15)

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

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

表2 計算任務(wù)

表3 算法完成任務(wù)耗時

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

4 全球擾動重力梯度全張量圖構(gòu)建

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

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

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

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

表4 全球擾動重力梯度全張量數(shù)值統(tǒng)計

5 結(jié)論

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

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

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

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

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