三維Tesseroid網(wǎng)格模型重力異常正演方法及并行算法

王博，郭良輝,2，崔亞彤，王祥

(1.中國(guó)地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院，北京 100083；2.地質(zhì)過(guò)程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室 中國(guó)地質(zhì)大學(xué)(北京)，北京 100083)

0 引言

三維網(wǎng)格模型的正反演是重力資料立體解釋的重要手段。其中，正演是反演的基礎(chǔ)，高精度、高效率的正演將促進(jìn)高質(zhì)量、高效率的反演解釋。對(duì)于小范圍、小尺度研究區(qū)，通常將地下三維空間剖分為眾多規(guī)則排列的直立長(zhǎng)方體組合，各長(zhǎng)方體的密度各不相同，然后對(duì)直立長(zhǎng)方體組合的三維網(wǎng)格模型進(jìn)行正演和反演。然而，對(duì)于大面積、大尺度研究區(qū)，由于地球曲率的存在，常規(guī)的基于直立長(zhǎng)方體組合的三維網(wǎng)格模型已不能適用，需采用球坐標(biāo)系的三維網(wǎng)格模型方法。

當(dāng)前球坐標(biāo)系重力異常正演方法較多，主要有下面4類：① 將球坐標(biāo)系下的單元體以點(diǎn)元、線元或者面元的方式對(duì)近似的規(guī)則形體進(jìn)行處理[1]，這類方法主要的問(wèn)題是計(jì)算精度非常低，難以滿足日益增加的精度要求；② 將全球地形剖分為眾多扇形環(huán)或球冠單元體等規(guī)則形體的組合，并給出嚴(yán)格的積分解析式，然后計(jì)算球坐標(biāo)系下各規(guī)則模型體的正演重力效應(yīng)，再累加作為全球地形改正值[2]，這類方法一般用于單點(diǎn)觀測(cè)，且隨著觀測(cè)點(diǎn)的變化，相對(duì)應(yīng)的球坐標(biāo)系以及地形模型需要重新構(gòu)建；③ 基于Tesseroid單元體的泰勒級(jí)數(shù)展開法，在正演模型體的幾何中心點(diǎn)對(duì)模型積分核函數(shù)進(jìn)行泰勒展開，計(jì)算球坐標(biāo)系下各Tesseroid單元體的重力效應(yīng)并疊加[1,3-4]，這種方法的計(jì)算精度與泰勒展開級(jí)數(shù)以及剖分網(wǎng)格的大小密切相關(guān)；④ 基于Tesseroid單元體的高斯—勒讓德(GLQ)積分法，給出橢球形地球的精確剖分方案，然后將基于數(shù)值積分方法的GLQ積分法應(yīng)用于橢球體積分計(jì)算中，同樣通過(guò)疊加得到最終重力效應(yīng)[5-7]，這種方法計(jì)算精度一般高于泰勒級(jí)數(shù)展開法。然而，不管是高斯—勒讓德積分法還是泰勒級(jí)數(shù)展開法，它們?cè)谡萦?jì)算地表或近地表觀測(cè)面的測(cè)點(diǎn)異常時(shí)，都會(huì)產(chǎn)生鄰近效應(yīng)，即誤差會(huì)隨著觀測(cè)距離的減小而急劇增大。王祥等和Hao等就此對(duì)高斯—勒讓德積分法作了改進(jìn),給出了新的單元體自適應(yīng)剖分方案，極大提高該算法在地表觀測(cè)面的計(jì)算精度，進(jìn)而分別給出了基于改進(jìn)高斯—勒讓德積分法(GLQ2D)的單一密度層界面深度反演和視密度填圖方法[8-9]。

本文研究將王祥等和Hao等[8-9]的單一密度層模型正演方法擴(kuò)展到多個(gè)密度層的三維網(wǎng)格模型，為大尺度數(shù)值模擬和三維反演奠定基礎(chǔ)。由于GLQ2D算法細(xì)化每個(gè)Tesseroid單元體剖分，導(dǎo)致剖分后的小單元體個(gè)數(shù)增加及正演計(jì)算量指數(shù)級(jí)增大，造成大規(guī)模三維Tesseroid網(wǎng)格模型的正演效率低下。

在重磁領(lǐng)域，國(guó)內(nèi)外學(xué)者給出了重力異常正反演的一些快速算法，主要有頻率域算法[10]，平移對(duì)稱算法[11]和并行計(jì)算[12-16]等。其中，并行計(jì)算主要分為兩大方向，主機(jī)端的CPU并行和設(shè)備端的GPU并行。CPU端的并行主要有OpenMP(open multi-processing)和MPI(message passing interface)等方式。OpenMP是一種用于共享內(nèi)存并行系統(tǒng)的多線程程序設(shè)計(jì)的一套指導(dǎo)性注釋，它可以使程序員把更多的精力投入到并行算法本身，而非其具體實(shí)現(xiàn)細(xì)節(jié)。而MPI則是信息傳遞接口，是獨(dú)立于語(yǔ)言的通信協(xié)議，它需要程序員手動(dòng)管理數(shù)據(jù)分配，實(shí)現(xiàn)進(jìn)程通信以及維持同步。GPU端的并行方式主流的有NIVIDIA CUDA架構(gòu)。按并行策略的不同，并行算法也分為兩種：一種是任務(wù)并行，另一種則是數(shù)據(jù)并行。任務(wù)并行的并行思想是把同一批數(shù)據(jù)分給for循環(huán)中不同的循環(huán)體，進(jìn)行處理。數(shù)據(jù)并行的思想則是不同的數(shù)據(jù)，用同一個(gè)程序處理。當(dāng)前，國(guó)內(nèi)外一些學(xué)者應(yīng)用并行算法一定程度上提高了重磁正反演效率。比如，陳召曦等[12]給出了基于NIVIDIA CUDA并行計(jì)算的重力及重力梯度數(shù)據(jù)正演算法，Hou等[13-15]給出了基于多級(jí)混合并行的重力梯度三維密度反演算法，周雪等[16]給出了基于MPI和OpenMP并行計(jì)算的重力及重力梯度正演算法。但這些算法都是基于直角坐標(biāo)系進(jìn)行的。

本文針對(duì)大尺度、地表觀測(cè)面研究區(qū)的高精度、高效率重力正演問(wèn)題，研究球坐標(biāo)系三維Tesseroid網(wǎng)格模型重力異常正演算法，采用GLQ2D算法提高大尺度、地表觀測(cè)面的重力異常計(jì)算精度，采用基于OpenMP的MATLAB任務(wù)并行算法提高正演效率。最后通過(guò)理論模型數(shù)據(jù)和華東巖石圈三維密度網(wǎng)格模型作正演試驗(yàn)，評(píng)價(jià)本文方法的有效性。

1 方法原理

基于球坐標(biāo)系的三維網(wǎng)格模型重力異常正演方法是將地下三維空間剖分為規(guī)則排列的Tesseroid單元體組合，各個(gè)單元體具有不同的密度值，再通過(guò)GLQ2D算法正演計(jì)算每個(gè)單元體在觀測(cè)面每個(gè)測(cè)點(diǎn)的重力異常，觀測(cè)面上任意測(cè)點(diǎn)的重力異常值即是所有Tesseroid單元體重力異常值在該測(cè)點(diǎn)處的累加。

1.1 Tesseroid單元體重力正演公式

在球坐標(biāo)系中，對(duì)于任意Tesseroid單元體(圖1)，其重力異常積分公式可由下式表示[1,6]:

圖1 Tesseroid單元體示意[8]

(1)

其中：

cosψ=sinφosinφs+cosφocosφscos(φs-φo)

式中：G表示萬(wàn)有引力常量，G=6.67×10-11N·m2/kg2，ρ表示模型體密度，單位為g/cm3，觀測(cè)點(diǎn)坐標(biāo)Q0(r0,φ0,λ0)，Tesseroid單元體中心點(diǎn)坐標(biāo)Q(rs,φs,λs)，l為觀測(cè)點(diǎn)Q0與單元體中心點(diǎn)Q之間的距離，ψ為Q0點(diǎn)與Q點(diǎn)相對(duì)于球心的夾角。

1.2 改進(jìn)的高斯—勒讓德積分法

將模型積分核函數(shù)積分區(qū)間化為勒讓德多項(xiàng)式，隨后使用近似積分核函數(shù)的勒讓德多項(xiàng)式計(jì)算累加得到模型異常值。對(duì)球坐標(biāo)系下重力異常積分表達(dá)式(1)進(jìn)行GLQ積分分解，得到Tesseroid單元體的勒讓德多項(xiàng)式表達(dá)式如下[5]：

(2)

式中：

Heck和Seitz[1]以及王祥等[8]針對(duì)這一問(wèn)題對(duì)GLQ積分法進(jìn)行了改進(jìn):

(3)

與傳統(tǒng)的GLQ積分法對(duì)比，Heck和Seitz[1]改進(jìn)后的GLQ積分法先對(duì)球坐標(biāo)系下重力異常積分表達(dá)式r方向積分，王祥等[8]進(jìn)一步作GLQ積分分解，實(shí)現(xiàn)了用較少的權(quán)重控制點(diǎn)來(lái)計(jì)算模型重力異常(即改進(jìn)后不在r方向展開)，既提升了模型重力計(jì)算的精度，又減少了模型計(jì)算的時(shí)間。這種改進(jìn)后的GLQ積分法稱之為改進(jìn)的高斯—勒讓德積分，即GLQ2D算法。

為了進(jìn)一步提高地表觀測(cè)面的正演精度，Hao等[9]在王祥等[8]的基礎(chǔ)上對(duì)前人的自適應(yīng)剖分方案進(jìn)行簡(jiǎn)化改進(jìn),引入根據(jù)不同計(jì)算精度要求而給定的距離判定值W。對(duì)一個(gè)徑向、緯向和經(jīng)向的計(jì)算范圍分別為[r1,r2]、[φ1,φ2]、[λ1,λ2]的Tesseroid單元體是否進(jìn)行進(jìn)一步的剖分，可以由判別式(4)來(lái)進(jìn)行判定。

(4)

其中：

L=max(Lφ,Lλ)

式中：Lφ，Lλ根據(jù)Tesseroid單元體的弧長(zhǎng)與周長(zhǎng)的比例得來(lái)。根據(jù)改進(jìn)后的自適應(yīng)剖分方案，有效的解決了近地表觀測(cè)時(shí)鄰近效應(yīng)帶來(lái)的精度不足的問(wèn)題。

1.3 三維Tesseroid網(wǎng)格模型重力正演流程及并行算法

本文球坐標(biāo)系下三維網(wǎng)格模型正演流程與直角坐標(biāo)系下三維網(wǎng)格模型正演的類似。圖2顯示了球坐標(biāo)系三維網(wǎng)格模型重力正演串行算法流程圖。

圖2 三維Tesseroid網(wǎng)格模型重力正演串行流程

本文正演串行算法步驟如下：① 輸入三維Tesseroid網(wǎng)格模型相關(guān)數(shù)據(jù)以及觀測(cè)面相關(guān)數(shù)據(jù)，模型數(shù)據(jù)包括模型經(jīng)度、緯度、深度方向的坐標(biāo)最小值、縱橫向步長(zhǎng)、Tesseroid網(wǎng)格數(shù)量(nx,ny,nz)以及各單元體密度值，觀測(cè)面數(shù)據(jù)包括測(cè)點(diǎn)總數(shù)(ns)及每個(gè)測(cè)點(diǎn)的經(jīng)度、緯度、高度。② 根據(jù)模型深度軸層數(shù)nz逐層遍歷；根據(jù)觀測(cè)面測(cè)點(diǎn)數(shù)ns逐點(diǎn)遍歷；進(jìn)一步，每一層根據(jù)模型規(guī)則網(wǎng)nx×ny個(gè)單元體按照先y軸方向、再x軸方向依次遍歷計(jì)算。按照上述遍歷順序，利用GLQ2D積分法計(jì)算單個(gè)單元體對(duì)單個(gè)觀測(cè)點(diǎn)的重力效應(yīng)，然后疊加單層所有單元體對(duì)單點(diǎn)重力效應(yīng)，最后疊加所有層單元體對(duì)單點(diǎn)重力效應(yīng)，也即得到模型在觀測(cè)面各測(cè)點(diǎn)的重力異常。③ 輸出該模型對(duì)觀測(cè)面重力正演結(jié)果。圖3是球坐標(biāo)系的多層Tesseroid單元體模型剖分示意。

圖3 三維Tesseroid網(wǎng)格模型

本文的并行思路采用基于OpenMP的MATLAB任務(wù)并行算法，即利用parfor對(duì)for循環(huán)進(jìn)行并行(圖4)。MATLAB進(jìn)行parfor循環(huán)時(shí)采用client和worker模式。其中client為編寫和啟動(dòng)該代碼的MATLAB端，而worker指運(yùn)行該代碼的MATLAB端。電腦中的MATLAB軟件可視為一個(gè)進(jìn)程，同一臺(tái)電腦可以同時(shí)運(yùn)行多個(gè)MATLAB進(jìn)程，每個(gè)worker對(duì)應(yīng)的物理單元為處理器或處理器核。每個(gè)MATLAB進(jìn)程之間可以通過(guò)一定方式開展數(shù)據(jù)傳輸。用戶首先在client端編寫所要運(yùn)行的代碼，client端在運(yùn)行該代碼的過(guò)程中，將需要并行的代碼段分配到其他MATLAB進(jìn)程運(yùn)行。

圖4 MATLAB parfor任務(wù)并行算法示意

對(duì)圖2所示的球坐標(biāo)系三維網(wǎng)格模型重力正演計(jì)算的并行策略可以有很多種，但效果基本相當(dāng)。本文根據(jù)最外層(模型深度層)循環(huán)并行加速最優(yōu)的原則，提出并行策略如下：首先明確需要并行的代碼段，即模型分層的循環(huán)結(jié)構(gòu)，然后通過(guò)任務(wù)并行加速在client端實(shí)現(xiàn)代碼設(shè)計(jì)，MATLAB client端自動(dòng)將模型不同層使用算法分配給不同worker處理單元，各處理單元分別計(jì)算單層重力效應(yīng)，最終疊加得到模型正演結(jié)果。

2 SLAB模型試驗(yàn)

2.1 模型參數(shù)

該系列模型為俯沖帶簡(jiǎn)單模型模擬，經(jīng)度、緯度范圍分別為100°E～110.2°E，20°N～30.2°N，深度方向范圍為0～200 km，俯沖帶用以模擬板塊匯聚邊緣，其頂部距觀測(cè)面的垂直距離為30 km，并以東傾60°的方向一直延伸至200 km，剩余密度為0.08 g/cm3。

本文設(shè)計(jì)了3個(gè)不同數(shù)量級(jí)的模型：① 模型大小51°×51°×20 km，經(jīng)度、緯度方向網(wǎng)格間距為0.2°×0.2°，深度方向步長(zhǎng)10 km；② 模型大小51°×51°×40 km，經(jīng)度、緯度方向網(wǎng)格間距為0.2°×0.2°，深度方向步長(zhǎng)5 km；③ 模型大小101°×101°×40 km(圖5)，經(jīng)度、緯度方向網(wǎng)格間距為0.1°×0.1°，深度方向步長(zhǎng)5 km。3個(gè)模型觀測(cè)面高度統(tǒng)一設(shè)置為0 m。

圖5 俯沖帶三維模型示意

2.2 模型正演與效果評(píng)價(jià)

本文在集群機(jī)上實(shí)現(xiàn)基于CPU的并行計(jì)算，集群CPU信息：Intel(R)Xeon(R)CPU E5-2680 v2 @ 2.80GHz，并行使用軟件為：MATLAB R2018a，模型并行計(jì)算單元數(shù)即并行使用CPU核數(shù)：12。圖6為GLQ2D串行正演結(jié)果，并行正演結(jié)果，以及模型24°N剖面結(jié)果對(duì)比。通過(guò)對(duì)比可見并行算法與串行算法正演結(jié)果完全一致。

對(duì)比圖6a、b，從數(shù)值來(lái)看，其異常結(jié)果均呈現(xiàn)中間異常高，兩側(cè)異常低，且俯沖帶右側(cè)異常降低較左側(cè)慢的特征。然后進(jìn)一步對(duì)比圖中24°N剖面正演曲線，如圖6c所示，重力異常中間高兩側(cè)低，俯沖帶左側(cè)異常降低較右側(cè)快，重力異常串行結(jié)果和并行結(jié)果完全吻合，驗(yàn)證了對(duì)于該數(shù)據(jù)模型使用并行計(jì)算不會(huì)改變它的計(jì)算結(jié)果和計(jì)算精度。

注：白色虛線方框?yàn)楦_帶沿傾角方向投影到觀測(cè)面上位置；黑色虛線為對(duì)比剖面位置

表1與表2分別對(duì)比了不同核數(shù)并行俯沖帶模型重力正演計(jì)算效率以及不同數(shù)量級(jí)俯沖帶模型的重力正演計(jì)算效率對(duì)比。在這里為了解釋方便，引入加速比的定義，即串行運(yùn)行用時(shí)比并行運(yùn)行時(shí)間[17]。

表1 不同并行核數(shù)重力正演計(jì)算效率對(duì)比

表2 不同數(shù)據(jù)量模型重力正演計(jì)算效率對(duì)比

從表2可以看出并行計(jì)算相比串行計(jì)算有明顯且穩(wěn)定的速度優(yōu)勢(shì)，且隨并行核數(shù)的增加，加速比不斷增加，在12核并行計(jì)算時(shí)得到最大平均加速比9.580。因此本文中其他正演算法均采用12核并行計(jì)算。

對(duì)一個(gè)并行系統(tǒng)(包含但不限于算法，程序)來(lái)說(shuō)，如果其性能(如加速比)可以隨處理單元數(shù)目的增加而按比例增加，我們稱該并行系統(tǒng)具有可擴(kuò)放性[17]。對(duì)于球坐標(biāo)系MATLAB重力異常正演并行算法，保持程序的可擴(kuò)放性是不可缺少的。對(duì)比表1中不同并行核數(shù)模型計(jì)算運(yùn)行時(shí)間及其加速比以及表2中不同數(shù)據(jù)量模型計(jì)算運(yùn)行時(shí)間及其加速比，可以得知該并行策略穩(wěn)定可行，即具有良好的可擴(kuò)放性，且隨計(jì)算數(shù)據(jù)量的增加，計(jì)算加速比不斷增加并趨于處理單元數(shù)。

3 華東三維巖石圈模型試驗(yàn)

本次正演試驗(yàn)所選取的研究區(qū)為中國(guó)東部，經(jīng)度范圍99.75°E～122.25°E,緯度范圍20.75°N～45.25°N。華東是我國(guó)地質(zhì)構(gòu)造和動(dòng)力學(xué)研究及資源能源勘查的重要區(qū)域，其區(qū)域地球物理場(chǎng)數(shù)值模擬和殼幔結(jié)構(gòu)成像具有重要的科學(xué)意義。本文從Shen等[18]中提取出研究區(qū)巖石圈三維橫波速度結(jié)構(gòu)模型，根據(jù)Brocher[19]的速度Vp、Vs和密度經(jīng)驗(yàn)式換算得到研究區(qū)巖石圈三維密度模型(圖7)。該華東巖石圈三維密度模型網(wǎng)格大小111°×121°×95 km，經(jīng)度、緯度方向網(wǎng)格間距為0.2°×0.2°，深度方向步長(zhǎng) 2 km。最小密度值為2.26 g/cm3,最大密度值為 3.45 g/cm3，平均密度為3.144 g/cm3。針對(duì)該模型，傳統(tǒng)的直角坐標(biāo)系計(jì)算受地球曲率的影響會(huì)產(chǎn)生明顯的誤差，因此需要球坐標(biāo)系計(jì)算。而該模型的球坐標(biāo)系正演串行算法耗時(shí)較長(zhǎng)(約15 h)，因此，本文采用并行算法來(lái)提高正演效率。

圖7 華東巖石圈三維密度模型

應(yīng)用本文方法對(duì)華東巖石圈三維模型進(jìn)行了正演計(jì)算，并截取模型30°N剖面進(jìn)行效果對(duì)比。圖8顯示了華東模型GLQ2D串行正演結(jié)果和12核并行正演結(jié)果，以及其30°N剖面的正演效果對(duì)比。通過(guò)對(duì)比可見模型并行正演與串行正演結(jié)果完全一致。研究區(qū)巖石圈密度模型的理論重力異常變化整體較為平緩，幅值范圍為-450～650 mGal。其中，內(nèi)蒙古中東部和河北北部及青藏高原一帶呈現(xiàn)大范圍的低值分布，揚(yáng)子克拉通呈現(xiàn)大范圍的高值分布，鄂爾多斯地塊、松遼盆地呈現(xiàn)中等范圍的高值分布，東南沿海一帶異常幅值整體在-100～200 mGal之間。

本文收集研究區(qū)EGM2008[20]的實(shí)際布格重力異常數(shù)據(jù)(圖8d)用于異常比較，數(shù)據(jù)網(wǎng)格間距為0.2°×0.2°。研究區(qū)EGM2008實(shí)際布格重力異?？傮w呈東高西低的趨勢(shì)，幅值范圍為-350～200 mGal。本文正演結(jié)果的整體異常特征與實(shí)際布格重力異常差別較大，究其原因主要有兩點(diǎn)：① 本文正演所依據(jù)的模型為0～200 km深度范圍內(nèi)的華東巖石圈三維密度模型，不包括周邊區(qū)域及200 km深度以下的介質(zhì)模型；② 本文所用的巖石圈三維密度模型是由背景噪聲成像獲得的巖石圈三維橫波速度模型[18]通過(guò)速度—密度經(jīng)驗(yàn)式[19]換算而來(lái)的，該模型沒有得到重力數(shù)據(jù)約束，而且所用的經(jīng)驗(yàn)式并不能準(zhǔn)確反映華東的實(shí)際速度和密度關(guān)系，不能反映研究區(qū)真實(shí)殼幔密度分布，從而導(dǎo)致上述的異常差別。

表3為華東模型運(yùn)行時(shí)間表，從表3可以看出，并行計(jì)算相比串行計(jì)算有明顯且穩(wěn)定的速度優(yōu)勢(shì)，且對(duì)比理論模型，該模型串行運(yùn)行時(shí)間均值比上12核并行計(jì)算運(yùn)行時(shí)間均值得到平均加速比為 11.206，遠(yuǎn)大于前者的9.580。這進(jìn)一步驗(yàn)證了該并行系統(tǒng)的可擴(kuò)放性，即隨計(jì)算數(shù)據(jù)量的增大，并行加速比隨之增大且趨于處理單元數(shù)。

表3 重力正演計(jì)算效率對(duì)比

4 結(jié)論

本文針對(duì)大尺度、地表觀測(cè)面的高精度、高效率重力正演問(wèn)題，給出了基于改進(jìn)的高斯—勒讓德積分法的并行計(jì)算方案，并通過(guò)理論模型和華東巖石圈三維模型數(shù)據(jù)試驗(yàn)驗(yàn)證了本文方法的有效性。數(shù)據(jù)試驗(yàn)表明本文方法可以實(shí)現(xiàn)大尺度、地表觀測(cè)面重力異常的高精度、高效率正演，基于OpenMP的MATLAB任務(wù)并行算法具有可擴(kuò)放性及穩(wěn)定性，即隨數(shù)據(jù)量的增大，加速比會(huì)趨于并行單元數(shù)；數(shù)據(jù)量不變，并行單元增大，加速比也會(huì)增大。本文方法為高效的大尺度重力場(chǎng)模擬和三維反演奠定技術(shù)基礎(chǔ)。