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

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

0 引言

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

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

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

1 基本原理

地球重力位可視為正常重力位與擾動(dòng)重力位之和：

W=U+T。

(1)

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

(2)

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

1.1 擾動(dòng)重力梯度計(jì)算

擾動(dòng)重力位的球諧展開(kāi)式如下：

(3)

(4)

其中：

(5)

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

1.2 計(jì)算公式簡(jiǎn)化

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

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

(6)

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

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

(7)

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

(8)

當(dāng)m=1時(shí):

(9)

當(dāng)m=2時(shí):

(10)

當(dāng)m≥3時(shí):

(11)

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

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

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

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

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

(12)

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

(13)

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

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

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

Txx+Tyy+Tzz=0。

(14)

EIGEN-6C4綜合了來(lái)自L(fǎng)AGEOS(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)計(jì)算平臺(tái)為戴爾Precision 7530移動(dòng)工作站，CPU為頻率2.9 GHz的Intel Xeon，內(nèi)存容量為32 GB，GPU為NVIDIA Quadro P3200，擁有1 792個(gè)CUDA核心。編譯環(huán)境為Visual Studio 2015-Visual C++，CUDA版本10.0。

3.1 可靠性檢驗(yàn)

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

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

表1 三分量檢核值

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

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

(15)

式中：t1為串行計(jì)算所消耗的時(shí)間；tn為n個(gè)線(xiàn)程參與計(jì)算時(shí)消耗的時(shí)間。

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

表2 計(jì)算任務(wù)

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

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

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

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

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

圖2 全球擾動(dòng)重力梯度全張量

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

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

5 結(jié)論

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

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

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

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

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