戴世坤,朱德祥*,張瑩,李昆,陳輕蕊,凌嘉宣,田紅軍

1 中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室,長沙 410083 2 中南大學(xué)有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,長沙 410083 3 中南大學(xué)地球科學(xué)與信息物理學(xué)院,長沙 410083 4 西南石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,成都 610500

0 引言

重力勘探廣泛用于油氣勘探、金屬礦勘查、地球內(nèi)部構(gòu)造、區(qū)域地質(zhì)構(gòu)造以及環(huán)境地球物理等研究領(lǐng)域.研究針對大規(guī)模復(fù)雜條件下重力異常場高效、高精度三維數(shù)值模擬對重力異常反演以及人機(jī)交互解釋有著重要作用.重力數(shù)值模擬方法按大類主要可分為空間域和頻率域方法.

空間域方法中,Nagy(1966)、Li和Chouteau(1998)及Talwani和Ewing(1960)分別推導(dǎo)了長方體、直立棱柱體以及多面體的解析解,徐世浙(1984)、Paul(1974)、Barnett(1976)、Kwok(1991)以及 García-Abdeslem(1992)利用高斯公式將體積分轉(zhuǎn)化為面積分;Ren等(2017)針對大規(guī)模三維重力異常正演開發(fā)了一種多極子算法,對復(fù)雜地形適應(yīng)性好,計(jì)算精度高,速度快;Ren等(2018)通過散度定理,將體積分轉(zhuǎn)化為線積分,導(dǎo)出線積分的精準(zhǔn)解,進(jìn)而得到重力場以及梯度張量精確解;Zhong等(2019)通過體積分與面積分的轉(zhuǎn)化,利用高斯-勒讓德積分公式計(jì)算多項(xiàng)式密度分布的幾何體重力場以及重力梯度張量,能很好地解決全球尺度的重力建模問題;也有學(xué)者利用數(shù)值模擬算法(Farquharson and Mosher,2009; Zhang et al.,2004; Cai and Wang,2005; Jahandari and Farquharson,2013)進(jìn)行三維重力異常數(shù)值模擬計(jì)算,取得了較好的效果.空間域方法對于異常體簡單,存在解析解的模型效果好,對復(fù)雜條件下重力場三維數(shù)值模擬頻率域方法具有更大的優(yōu)勢(Pedersen,1978,1985; 熊光楚,1984; 馮銳,1986;Chai and Hinze,1988; Parker,1973; 徐世浙,2007;王萬銀等,2009; 王萬銀,2012;Wu and Tian,2014; Dai et al.,2018; 戴世坤等,2020).相比于空間域,頻率域表達(dá)式簡潔,計(jì)算效率高效,但其精度受傅里葉變換的影響.快速傅里葉變換算法(FFT)的出現(xiàn)提升了離散傅里葉變換的計(jì)算效率,使得頻率域重力異常正演計(jì)算效率高.但標(biāo)準(zhǔn)傅里葉變換存在截?cái)嘈?yīng),通過擴(kuò)邊能降低這種影響(Caratori Tontini et al.,2009);針對精度低問題,柴玉璞(1996)證明了傅里葉變換離散化定理,該定理揭示了函數(shù)離散值與譜離散值之間的關(guān)系,并據(jù)此給出了傅里葉變換數(shù)值計(jì)算理論-移樣離散傅里葉變換理論(柴玉璞,1988; Chai,2009),并將移樣離散傅里葉變換理論(SFT)應(yīng)用到重磁波數(shù)域正演中,極大提升了波數(shù)域正演精度(柴玉璞和萬海珍,2020);吳樂園(2014)、Wu和Tian(2014)、Wu(2016)將SFT和高斯積分節(jié)點(diǎn)及權(quán)系數(shù)相結(jié)合,利用高斯節(jié)點(diǎn)積分的高精度特點(diǎn),提出了Gauss-FFT方法,避免了零波數(shù)計(jì)算,降低了截?cái)嘈?yīng)的影響,在波數(shù)域正演中取得了很好的效果(Wu et al.,2019; Dai et al.,2022);Zhou等(2020)利用Gauss-FFT和NUFFT相結(jié)合進(jìn)一步提高了頻率域重力異常正演的適應(yīng)性和計(jì)算效率;周印明等 (2020)結(jié)合傅里葉變換積分與三次樣條插值函數(shù),在此基礎(chǔ)上進(jìn)行重力異常正演,降低了連續(xù)介質(zhì)重磁場數(shù)值模擬中截?cái)嘈?yīng)的影響,但文中并未提及波數(shù)選取規(guī)律.

Dai等(2018)提出基于微分法的空間-波數(shù)域重力異常三維數(shù)值模擬方法,采用高斯FFT/擴(kuò)邊標(biāo)準(zhǔn)FFT實(shí)現(xiàn)了三維重力異常場快速精準(zhǔn)數(shù)值模擬.該方法具有很好的并行性,但文獻(xiàn)并未實(shí)現(xiàn)方法的并行加速方案.目前,并行計(jì)算分為CPU并行和GPU并行.CPU并行受限于計(jì)算機(jī)線程規(guī)模以及線程調(diào)度/線程同步耗時(巴振寧等,2022);相比于CPU,GPU具有更多的計(jì)算單元,可以將大規(guī)模節(jié)點(diǎn)計(jì)算分配給GPU線程并行執(zhí)行,但是GPU僅適用于簡單的計(jì)算,不適用于復(fù)雜循環(huán)運(yùn)算(馬琳等,2022).因此,綜合多核多線程CPU以及眾核GPU線程的異構(gòu)體系結(jié)構(gòu)是高性能計(jì)算的發(fā)展趨勢(盧風(fēng)順等,2011;陳召曦等,2012;劉守偉等,2013;Qin et al.,2014; 徐云耘等,2022).

綜上,本文基于微分方程的空間-波數(shù)域重力異常數(shù)值模擬算法,進(jìn)行兩個方面的工作:其一是實(shí)現(xiàn)一種基于二次插值形函數(shù)的任意傅里葉變換算法,可以根據(jù)頻譜變化趨勢,合理采樣,提高傅里葉變換剖分的靈活性,單元內(nèi)原函數(shù)采用二次插值形函數(shù)進(jìn)行擬合,求得單元積分的解析表達(dá)式,積分精度高,能減小傅里葉變換帶來的截?cái)嘈?yīng).其二是利用空間-波數(shù)域算法的高度并行性,解剖基于任意傅里葉變換算法的空間-波數(shù)域重力異常正演算法,實(shí)現(xiàn)CPU-GPU的異構(gòu)體系結(jié)構(gòu)并行算法,在模型正演時由CPU來控制主程序,加速了五對角方程組求解(CPU并行)和任意傅里葉變換(GPU并行),提升了算法效率.設(shè)計(jì)模型驗(yàn)證了算法的正確性,給出的采樣規(guī)律適用于任意復(fù)雜條件下三維重力異常正演;測試了算法的并行加速比(串行計(jì)算耗時/并行計(jì)算耗時),相比CPU串行算法,CPU-GPU并行算法的速度提高了40倍以上.最后將新算法應(yīng)用到實(shí)際地形重力異常正演,計(jì)算規(guī)模為721×721×421個節(jié)點(diǎn)的模型,耗時僅需17 s,計(jì)算效率高,適合大規(guī)模復(fù)雜條件下的重力異常數(shù)值模擬.

1 理論方法

1.1 空間-波數(shù)域重力微分法正演

基于微分法的空間-波數(shù)域重力異常正演方法從泊松方程出發(fā),沿水平方向進(jìn)行二維傅里葉變換,得到不同波數(shù)下的常微分方程,結(jié)合準(zhǔn)確的邊界條件,采用有限單元法進(jìn)行求解,將三維問題轉(zhuǎn)化為多個一維問題,減少了計(jì)算量和存儲量需求,有效地提高了三維重力異常正演的計(jì)算效率.原理詳見文獻(xiàn)(Dai et al.,2018),本文僅簡要介紹.

重力位U與空間域密度滿足泊松方程(Blakely,1995):

(1)

其中,U為重力位,(x,y,z)為空間觀測點(diǎn)坐標(biāo),G為萬有引力常量6.67×10-11N·m3·kg-2,ρ為空間域密度.對式(1)沿水平方向做二維傅里葉正變換,得到空間-波數(shù)域常微分方程:

(2)

式(2)結(jié)合上下邊界條件得到空間-波數(shù)域重力異常正演的邊值問題(Dai et al.,2018),采用基于二次插值的一維有限元方法求解(徐世浙,1994),采用追趕法(續(xù)小磊和馬丁,2013)求解最終構(gòu)成的五對角線性方程組,得到空間-波數(shù)域重力位.根據(jù)重力位和場之間的關(guān)系可得到空間-波數(shù)域重力場,最后利用水平方向二維傅里葉反變換得到空間域重力異常場.

在進(jìn)行空間-波數(shù)域重力異常正演計(jì)算時,正反傅里葉變換和求解方程組兩部分具有高度并行性,且二維傅里葉正反變換耗時占整個算法的90%.Dai等(2018)采用標(biāo)準(zhǔn)FFT和高斯FFT實(shí)現(xiàn)傅里葉變換,這兩種方法都需均勻剖分,采樣不靈活;且文獻(xiàn)(Dai et al.,2018)未實(shí)現(xiàn)算法并行.為了實(shí)現(xiàn)任意采樣及進(jìn)一步提升計(jì)算效率,本文引入一種基于二次插值形函數(shù)的傅里葉變換算法,該算法采樣靈活,通過分析重力異常離散譜確定譜能量集中的范圍,能精準(zhǔn)、靈活選取波數(shù),且基于二次插值形函數(shù)求得傅里葉變換積分的解析解,計(jì)算精度高.將正反傅里葉變換和求解方程組兩部分基于CPU-GPU的異構(gòu)體系結(jié)構(gòu)實(shí)現(xiàn)了并行加速,提高了整個算法的計(jì)算效率.

1.2 任意傅里葉變換

1.2.1 任意傅里葉變換理論

本文空間-波數(shù)域重力異常場數(shù)值模擬算法中采用的二維傅里葉變換對(Blakely,1995)如下,

(3)

式中,f(x,y)表示空間域場,F(kx,ky)表示對應(yīng)的二維頻譜.

本文以二維傅里葉正變換式(3)為例,說明本文基于二次插值形函數(shù)的任意傅里葉變換算法的計(jì)算過程.

在實(shí)際應(yīng)用中,模擬區(qū)域有限,設(shè)x方向范圍:xmin～xmax,y方向范圍:ymin～ymax,此時式(3)可以寫為

(5)

本文提出的基于二次插值形函數(shù)的傅里葉變換算法首先將式(5)的二重積分轉(zhuǎn)化為兩個一重積分,可寫為

(6)

(7)

式中,Fx(kx,y)為對x積分得到的譜函數(shù).計(jì)算式(6)和式(7)兩個一重積分,可得二維傅里葉變換的二維頻譜.

將計(jì)算區(qū)域采用均勻或非均勻的剖分方式沿著x和y方向離散為Mx×My個單元,每個單元內(nèi)采用二次插值的形函數(shù)擬合空間域場,則節(jié)點(diǎn)總數(shù)為Nx×Ny.以式(6)為例,可寫為

(8)

式中ej表示第j單元,令第j個單元的三個節(jié)點(diǎn)坐標(biāo)為(xj1,xj2,xj3),xj2為單元ej中間節(jié)點(diǎn),滿足關(guān)系式xj1+xj3=2xj2.單元內(nèi)采用二次插值的形函數(shù)擬合函數(shù)f(x,y)(徐世浙,1994),則有

f(x,y)=Nj1f(xj1,y)+Nj2f(xj2,y)+Nj3f(xj3,y),

(9)

(10)

令

(11)

式(11)中I(xj1,kx),I(xj2,kx),I(xj3,kx)為沿x方向第j個單元上三個節(jié)點(diǎn)的傅里葉變換積分系數(shù),該積分系數(shù)可求得解析解.式(10)可寫為

(12)

式中,

(13)

同理,式(7)可寫為

(14)

式中,J(yj1,ky),J(yj2,ky),J(yj3,ky)分別為沿y方向上第j個單元的傅里葉變換節(jié)點(diǎn)積分系數(shù),式中(yj1,yj2,yj3)分別為y方向第j個單元的三個節(jié)點(diǎn)坐標(biāo).節(jié)點(diǎn)積分系數(shù)推導(dǎo)計(jì)算方法與I(xj1,kx),I(xj2,kx),I(xj3,kx)相同,此處不再贅述,直接寫出表達(dá)式:

(15)

從式(13)和式(15)可以看出,當(dāng)x和y方向剖分固定時,單元形函數(shù)不變,兩個方向的節(jié)點(diǎn)積分系數(shù)可提前計(jì)算出并存儲,重復(fù)利用,減少冗余計(jì)算,提高計(jì)算效率,并且不同波數(shù)下積分系數(shù)相互獨(dú)立,可并行性高.對于二維傅里葉反變換式(4),推導(dǎo)過程與正變換類似,不再贅述.

1.2.2 采樣規(guī)則

任意采樣點(diǎn)數(shù)FFT算法提出了離散頻率取值方式(陳龍偉等,2016),標(biāo)準(zhǔn)FFT存在截?cái)嘈?yīng).擴(kuò)邊FFT能減小這種誤差,其本質(zhì)是擴(kuò)邊之后波數(shù)域基頻減小,提高了計(jì)算精度;Gauss-FFT作為一種結(jié)合高斯點(diǎn)數(shù)的偏移采樣傅里葉變換算法,截?cái)嘈?yīng)小,計(jì)算精度高(吳樂園,2014),其原因是Gauss-FFT基于FFT采樣定理,在基頻單元內(nèi)將采樣點(diǎn)加密并偏移,加密點(diǎn)數(shù)為高斯點(diǎn)個數(shù),偏移距離為高斯偏移距離,能有效提高積分計(jì)算精度.這兩種FFT算法在計(jì)算離散頻譜時均利用了采樣定理,要求空間域和頻率域均勻采樣,無法實(shí)現(xiàn)非均勻采樣.此外,為滿足精度要求,擴(kuò)邊FFT和4個點(diǎn)的Gauss-FFT需要更多的波數(shù),降低了計(jì)算效率.

本文提出的基于二次插值形函數(shù)的任意傅里葉變換算法,在空間域和波數(shù)域任意剖分,且能根據(jù)波譜能量分布靈活采樣.其中,波數(shù)選取范圍可延用采樣定理,由于波數(shù)關(guān)于零值點(diǎn)的對稱性,所以在波數(shù)選取時令正波數(shù)與負(fù)波數(shù)互為相反數(shù),這里只給出正波數(shù)選取規(guī)則.

由采樣定理可知,波數(shù)最大值由最小網(wǎng)格間距決定.設(shè)空間域最小網(wǎng)格間距為Δxmin,Δymin,可知kx波數(shù)的正波數(shù)采樣范圍為[Δk,kxmax],ky波數(shù)的正波數(shù)采樣范圍為[Δk,kymax],式中:

(16)

基頻Δk作為基本單位:

(17)

式(17)中tmax為空間域x或y方向計(jì)算范圍的最大值,tmin為空間域x或y方向計(jì)算范圍的最小值.

在求得正波數(shù)采樣范圍后,可根據(jù)頻譜能量分布總結(jié)出波數(shù)選取的規(guī)律.本文總結(jié)波數(shù)選取規(guī)律的具體準(zhǔn)則為:頻譜能量大且變化劇烈的區(qū)間加密采樣,頻譜能量小且變化緩慢的區(qū)間稀疏采樣.

1.3 CPU-GPU并行理論方法

本文提出的基于任意傅里葉變換的空間-波數(shù)域三維重力異常正演算法中,采用任意傅里葉變換計(jì)算時不同波數(shù)下離散單元積分相互獨(dú)立,不同波數(shù)下求解五對角方程組相互獨(dú)立,這兩個關(guān)鍵過程均具有良好的并行性.本文采用OpenMP并行實(shí)現(xiàn)追趕法求解五對角方程組,采用GPU并行實(shí)現(xiàn)正、反二維傅里葉變換過程,實(shí)現(xiàn)整體算法的CPU-GPU并行加速.下面具體介紹這兩個并行理論和結(jié)構(gòu).

每個波數(shù)需進(jìn)行方程組求解,在求解五對角方程組時,追趕法算法雖然計(jì)算量小,但存在大量循環(huán)和迭代運(yùn)算,計(jì)算復(fù)雜,GPU相對于CPU控制單元數(shù)量少,指令控制能力差,且GPU適用于大量的輕量級運(yùn)算(巴振寧等,2022),追趕法算法存在復(fù)雜的指令流控制和復(fù)雜的公式計(jì)算,因此將五對角方程追趕法求解計(jì)算在CPU上采用OpenMP并行計(jì)算,并行過程如圖1所示:

圖1 OpenMP并行機(jī)制

如圖所示,主線程開始運(yùn)行,分配多個子線程執(zhí)行并行計(jì)算任務(wù),主線程等待子線程完成所有的并行計(jì)算任務(wù),然后繼續(xù)后續(xù)工作(蔡佳佳等,2007).OpenMP并行有豐富的指令,易于實(shí)現(xiàn),算法改造簡單高效(劉凱和寇正,2003).

相對于CPU,GPU可以同時開辟成千上萬個線程,且存儲器讀寫效率更高,在NIVIDIA CUDA編程框架中,將大規(guī)模的簡單數(shù)據(jù)計(jì)算分配到GPU線程上并行執(zhí)行,能有效地提升計(jì)算速度.本文中的基于二次插值形函數(shù)的任意傅里葉變換算法計(jì)算過程主要以乘法和加法為主,且每個節(jié)點(diǎn)的傅里葉變換之間相互獨(dú)立,適用于GPU并行計(jì)算.下面具體介紹GPU并行方案.

以二維形函數(shù)傅里葉變換正變換為例,將空間域x和y方向離散為Nx×Ny個節(jié)點(diǎn),節(jié)點(diǎn)的積分系數(shù)提前算出并存儲,剩下的計(jì)算分為兩步:①如圖2所示,先對每個波數(shù)下f(x,y)沿x方向做Nx次乘法和加法,得到Fx(kx,y);②再對Fx(kx,y)沿y方向做Ny次乘法和加法得到F(kx,ky),得到該波數(shù)下的頻譜值.若波數(shù)選取總數(shù)為Nkx×Nky,模型垂向剖分為Nz層,則算法中傅里葉正變換所需要加法和乘法的計(jì)算次數(shù)均為:Nz×(Ny×Nkx×Nx+Nkx×Nky×Ny).由于不同波數(shù)之間的傅里葉變換積分相互獨(dú)立,可以將外面三層循環(huán)(Nz×Ny×Nkx或Nz×Nkx×Nky)分配給GPU線程并行計(jì)算,每個GPU線程僅負(fù)責(zé)一個波數(shù)的積分運(yùn)算:Nx或Ny次乘法和加法運(yùn)算.

圖2 x方向一維積分計(jì)算過程

GPU并行部分稱為核函數(shù),CPU上啟動GPU上執(zhí)行,本文采用的GPU通用技術(shù)通過線程網(wǎng)格grid以及線程塊block來組織核函數(shù)并行計(jì)算部分,grid以及block有三個方向參數(shù):x,y,z.常用的線程層次有一維層次(一維線程塊與一維線程網(wǎng)格)、二維層次(二維線程塊與二維線程網(wǎng)格)和三維層次(三維線程塊與三維線程網(wǎng)格)(董廷星等,2009;賈格等,2017).

以二維線程層次為例,如圖3所示,每個線程有各自的核函數(shù)計(jì)算任務(wù),二維線程層次中線程通過(threadIdx%x,threadIdx%y)來定位自己在線程塊中的索引,線程塊通過(blockIdx%x,blockIdx%y)來定位自己在線程網(wǎng)格中的索引,兩個索引準(zhǔn)確指向線程數(shù)據(jù)對應(yīng)的邏輯位置,實(shí)現(xiàn)并行部分和主程序的對接.

圖3 二維線程層次結(jié)構(gòu)

圖3表示f(x,y)沿x方向做傅里葉變換積分得到Fx(kx,y)的過程.從圖3中可以看出,每一個線程塊中有a×b個線程,線程網(wǎng)格中有m×n個線程塊,線程塊中的一個GPU線程負(fù)責(zé)一個波數(shù)的積分計(jì)算,因此有m=(a+Nkx-1)/a,n=(b+Ny-1)/b,Nkx、Ny分別為沿x方向積分中kx和y方向網(wǎng)格剖分節(jié)點(diǎn)數(shù).一維層次中,每一個線程塊中有a個線程,線程網(wǎng)格中有m個線程塊,有m=(a+Nkx×Ny-1)/a.三維層次增加了z方向的維度,用來計(jì)算不同z平面的任意傅里葉變換積分,每一個線程塊中有a×b×c個線程,線程網(wǎng)格中有m×n×k個線程塊,有m=(a+Nkx-1)/a,n=(b+Ny-1)/b,k=(z+Nz-1)/c,Nz為z方向節(jié)點(diǎn)個數(shù).

2 數(shù)值算例

本節(jié)主要介紹任意傅里葉變換在重力正演時的波數(shù)采樣規(guī)則,驗(yàn)證算法的正確性與效率.測試所用PC參數(shù)為CPU Intel(R) Core(TM) i9-9980XE,主頻3.0 GHZ(18核),128 GB RAM.GPU顯卡型號NVIDIA TITAN RTX,顯存24G,計(jì)算能力7.5,CUDA版本號11.6,基于Linux編譯環(huán)境.

2.1 正確性驗(yàn)證

引入相對均方根誤差(RRMS)來研究基于任意傅里葉變換的重力異常數(shù)值模擬精度,當(dāng)RRMS<1%時,即認(rèn)為精度滿足要求.計(jì)算公式(Wu,2016)如下:

(18)

設(shè)計(jì)棱柱體模型,計(jì)算區(qū)域范圍20 km(x)×20 km(y)×7.5 km(z),異常體水平方向上位于模擬區(qū)域中心,頂面埋深為0.5 km,大小為8 km(x)×8 km(y)×2 km(z),異常體密度為2000 kg·m-3.剖分節(jié)點(diǎn)個數(shù)為251(x)×251(y)×151(z),三個方向上均勻剖分,水平方向網(wǎng)格間距均為80 m,垂直方向?yàn)?0 m.

根據(jù)1.2節(jié)總結(jié)的采樣規(guī)則,本節(jié)模型選取kx和ky波數(shù)相同,基頻大小均為Δk=3.14×10-4,波數(shù)最大值為3.9×10-2,在頻譜能量變化劇烈的區(qū)間進(jìn)行加密,加密范圍為前三個基頻區(qū)間:1×10-5～9×10-4,選取25個波數(shù),其他區(qū)間均勻分布,正負(fù)波數(shù)選取181×181個.將三維重力異常全息傅數(shù)值模擬異常體內(nèi)部(z=0 m)、異常體外部(z=500 m)結(jié)果與解析解對比.

圖4為重力異常場在地面z=0 m和異常體內(nèi)部z=500 m數(shù)值解與解析解的平面等值線圖.由圖可知,兩個平面上數(shù)值解和解析解等值線圖重合,場值吻合得很好,且在z=0 m平面gx、gy以及gz三個分量的相對均方根誤差分別為0.48%,0.48%,0.52%,z=500 m平面gx、gy以及gz三個分量的相對均方根誤差分別為0.61%,0.61%,0.65%,異常體內(nèi)部和外部的重力異常場的相對均方根誤差均小于1%,表明基于任意傅里葉變換的空間-波數(shù)域三維重力異常算法正確,波數(shù)選取規(guī)則對異常體內(nèi)外部的場均適用.

圖4 重力異常場數(shù)值解與解析解等值線圖

2.2 計(jì)算精度與計(jì)算效率

設(shè)計(jì)連續(xù)密度模型,研究基于任意傅里葉變換空間-波數(shù)域三維重力異常算法的計(jì)算效率與精度.文獻(xiàn)表明(吳樂園,2014;Dai et al.,2018)Gauss-FFT法對頻率域重力異常數(shù)值模擬具有非常高的精度,本節(jié)以高斯點(diǎn)數(shù)NG=4的高斯快速傅里葉變換的空間-波數(shù)域重力異常場數(shù)值模擬結(jié)果(Dai et al.,2018)作為連續(xù)介質(zhì)的參考解.研究本文算法基于不同波數(shù)個數(shù)的計(jì)算精度和計(jì)算效率,算法在CPU上串行計(jì)算.

連續(xù)密度模型垂直方向密度不變,在水平方向上密度變化可表示為

ρ=2000×e-6×10-8(x2+y2).

(19)

模擬區(qū)域范圍為20 km(x)×20 km(y)×7.5 km(z).異常體區(qū)域大小為20 km(x)×20 km(y)×2 km(z),頂面埋深為0.5 km.區(qū)域剖分節(jié)點(diǎn)個數(shù)為:201(x)×201(y)×151(z),水平方向網(wǎng)格間距均為100 m,垂直方向50 m.kx和ky的波數(shù)范圍為-3.14×10-2～3.14×10-2,基頻大小均為3.14×10-4,波數(shù)個數(shù)變化如表1所示.

表1 任意傅里葉變換與高斯快速傅里葉變換計(jì)算精度與計(jì)算效率對比

采用任意傅里葉變換算法進(jìn)行數(shù)值模擬,表1為任意傅里葉變換和Gauss-FFT空間-波數(shù)域重力異常正演計(jì)算精度和計(jì)算效率對比表.表中統(tǒng)計(jì)了不同波數(shù)個數(shù)下:地面z=0 m重力場的相對均方根誤差、整個區(qū)域計(jì)算耗時和內(nèi)存占用.從表中可以看出,①隨著波數(shù)增多,重力異常場三個分量的RRMS逐漸減少,基于任意傅里葉變換的空間-波數(shù)域算法正演耗時增多,內(nèi)存占用也增多,當(dāng)波數(shù)大于99×99時,相對均方根誤差基本不變;②波數(shù)選取個數(shù)為99×99時,以Gauss-FFT的數(shù)值解為參照,任意傅里葉變換正演得到gx、gy、gz三分量的RRMS分別為0.0298%,0.0298%和0.02%,與Gauss-FFT數(shù)值解結(jié)果非常接近;③模型剖分節(jié)點(diǎn)為201×201×151、計(jì)算節(jié)點(diǎn)數(shù)為201×201×151時,任意傅里葉變換計(jì)算耗時約3 s,占用內(nèi)存約469 MB,Gauss-FFT法耗時23 s,占用內(nèi)存約6144 MB,相比于Gauss-FFT,任意傅里葉變換算法效率提高了約8倍.

圖5是波數(shù)為99×99時任意傅里葉變換與Gauss-FFT兩種數(shù)值解的重力異常場z=0 m、z=500 m平面等值線圖.從圖中可以看出,異常體內(nèi)外部兩者結(jié)果基本吻合,表明本文算法對連續(xù)模型適應(yīng)強(qiáng)且計(jì)算效率高.

圖5 重力異常場任意傅里葉變換數(shù)值解與高斯數(shù)值解等值線圖

3 CPU-GPU并行加速

通過1.3節(jié)分析,并行加速有兩部分:第一部分:五對角方程組求解;第二部分:二維任意傅里葉正、反變換.采用2.2節(jié)中的模型,該模型已驗(yàn)證精度,改變計(jì)算節(jié)點(diǎn)數(shù),對比不同GPU線程層次加速任意傅里葉變換效率;使用多核CPU、GPU分別加速五對角方程組求解與任意傅里葉變換,改變計(jì)算節(jié)點(diǎn)數(shù),對比CPU與GPU的加速比.

3.1 GPU線程層次加速對比

一維層次、二維層次以及三維層次GPU結(jié)構(gòu)加速任意傅里葉變換,同時使用culbas庫直接加速任意傅里葉變換.保持z方向節(jié)點(diǎn)數(shù)為151,改變水平方向計(jì)算節(jié)點(diǎn)個數(shù),波數(shù)個數(shù)與空間域節(jié)點(diǎn)數(shù)保持一致,對比四種GPU方法加速任意傅里葉變換的計(jì)算時間以及GPU顯存-CPU內(nèi)存?zhèn)鬏敃r間,結(jié)果如表2所示.

表2 不同GPU方法加速任意傅里葉變換計(jì)算時間

從表中可以看出,不同計(jì)算節(jié)點(diǎn)個數(shù)下,一維、二維、三維層次的計(jì)算時間相近,但是二維層次的顯存-內(nèi)存?zhèn)鬏敽臅r最少,二維層次下的計(jì)算總時間少于一維和三維層次的計(jì)算總時間,同時一維、二維、三維的計(jì)算總時間優(yōu)于CUBLAS庫.因此,任意傅里葉變換GPU加速方案采用二維層次.

3.2 CPU、GPU加速方案對比

針對并行加速兩部分,分別采用CPU單核、18核CPU以及GPU進(jìn)行數(shù)值模擬,保持z方向節(jié)點(diǎn)數(shù)為151,改變水平方向計(jì)算節(jié)點(diǎn)個數(shù),波數(shù)個數(shù)與空間域節(jié)點(diǎn)數(shù)保持一致,對比不同加速方法的加速比,CPU、GPU并行加速五對角方程組求解時間如表3所示.

表3 CPU、GPU加速五對角方程組計(jì)算時間

從表中可以看出,五對角方程組求解過程中,CPU加速效果好于GPU加速,這是因?yàn)槊總€波數(shù)下的五對角方程組求解過程復(fù)雜,并不是簡單的點(diǎn)對點(diǎn)簡單運(yùn)算,隨著計(jì)算節(jié)點(diǎn)數(shù)增加,加速比最高可達(dá)13.6,因此,采用CPU加速五對角方程組求解.

CPU、GPU并行加速任意傅里葉變換時間如表4所示.

表4 CPU、GPU加速任意傅里葉變換計(jì)算時間

從表中可以看出,任意傅里葉變換時間占整個計(jì)算過程的90%,計(jì)算節(jié)點(diǎn)較少時,CPU加速效果優(yōu)于GPU加速,這是因?yàn)橛?jì)算量小,GPU處理性能并沒有體現(xiàn)出來(陳召曦等,2012),隨著計(jì)算節(jié)點(diǎn)數(shù)的增加,18核的CPU加速比最高達(dá)到12.59,GPU的加速比越來越大,最高達(dá)到48.77,CPU線程核數(shù)與加速比并不是線性關(guān)系,這是因?yàn)镺penMP并行時不同線程間數(shù)據(jù)存在依賴關(guān)系,線程同步、線程調(diào)度消耗時間,而GPU擁有眾多計(jì)算單元,每個GPU線程負(fù)責(zé)一個波數(shù)的積分運(yùn)算,隨著計(jì)算節(jié)點(diǎn)的增加,加速比也會越來越大,GPU加速效果優(yōu)于CPU加速.

因此本文采用OpenMP并行實(shí)現(xiàn)追趕法求解五對角方程組,采用二維層次GPU并行實(shí)現(xiàn)二維任意傅里葉正、反變換,實(shí)現(xiàn)整體算法的CPU-GPU并行加速方案.

4 實(shí)際地形應(yīng)用

長白山火山地區(qū)具有豐富的地?zé)豳Y源.基于長白山火山區(qū)域的DEM數(shù)字高程數(shù)據(jù)(經(jīng)度為127°E至129°E,緯度為41°N至43°N)設(shè)計(jì)實(shí)際地形應(yīng)用模型,采用基于任意傅里葉變換的空間-波數(shù)域重力異常正演CPU-GPU并行算法進(jìn)行數(shù)值模擬.地形數(shù)據(jù)如圖6所示,圖中黑框內(nèi)的范圍為目標(biāo)區(qū)域,范圍為100×100 km,其他區(qū)域?yàn)閿U(kuò)邊區(qū)域,總體范圍為200×200 km.為了將地形描述清楚,在地形起伏大的目標(biāo)區(qū)域(x方向范圍:-20～20 km,y方向范圍:-20～20 km)采用100 m網(wǎng)格間距,其他區(qū)域采用500 m網(wǎng)格間距.海拔最高為2727 m,最低為257 m.z方向向下為正,z方向范圍-4～4 km,-4～0 km垂直方向網(wǎng)格大小為10 m,0～4 km網(wǎng)格大小為200 m,模擬區(qū)域節(jié)點(diǎn)個數(shù)為721(x)×721(y) ×421(z),密度設(shè)為常密度,數(shù)值為2580 kg·m-3.

圖6 長白山火山區(qū)域地形圖

根據(jù)1.2節(jié)中的波數(shù)選取規(guī)律,本節(jié)模型選取的kx和ky波數(shù)相同,兩個方向上基頻大小均為Δk=3.14×10-5,波數(shù)最大值為3.14×10-2,在前10個基頻1×10-6～3×10-4進(jìn)行波數(shù)加密,選取201×201、401×401兩種波數(shù)進(jìn)行數(shù)值模擬.將起伏地形看成721×721個密度為2580 kg·m-3的常密度棱柱體模型組成,可近似求得起伏地形重力異常的解析解,z=-4 km為觀測面.

對整個區(qū)域進(jìn)行數(shù)值模擬,區(qū)域節(jié)點(diǎn)個數(shù)為721×721×421,計(jì)算節(jié)點(diǎn)個數(shù)為521×521×421.當(dāng)波數(shù)個數(shù)為201×201時,CPU串行計(jì)算時間為171.62 s,占用內(nèi)存為12.3 G,CPU-GPU并行計(jì)算時間為7.81 s,CPU+GPU占用內(nèi)存為12.3 G+215 MB,加速比為22,z=-4 km重力場三分量的RRMS為0.34%,0.33%,0.11%;波數(shù)個數(shù)為401×401時,CPU串行計(jì)算時間為470 s,占用內(nèi)存為17.3 G,CPU-GPU并行計(jì)算時間為17.68 s,CPU+GPU占用內(nèi)存為17.3G+276 MB,加速比為26.5,z=-4 km重力場三分量的RRMS為0.0063%,0.0063%,0.0035%.圖7為觀測面z=-4 km上波數(shù)選取401×401的重力異常場gx、gy以及gz三分量數(shù)值解和解析解對比圖,從圖中可以看出,數(shù)值模擬與解析解計(jì)算結(jié)果形態(tài)相似.

圖7 實(shí)際地形重力異常場數(shù)值解以及解析解結(jié)果

為了對比本文算法計(jì)算精度,采用NG=4的Gauss-FFT方法(Dai et al.,2018)進(jìn)行數(shù)值模擬,CPU串行計(jì)算時間為6234.23 s,占用內(nèi)存為113.75 G,與解析解對比,z=-4 km重力場三分量的RRMS為0.035%,0.035%,0.017%,當(dāng)任意傅里葉變換的波數(shù)選取個數(shù)為401×401時,本文算法的計(jì)算精度高于NG=4的Gauss-FFT計(jì)算精度,計(jì)算時間僅需17.68 s,適用于復(fù)雜條件下重力模型計(jì)算.

綜上,計(jì)算結(jié)果證明本文算法及加速方案適用于大規(guī)模起伏地形模型的計(jì)算,證明1.2節(jié)中總結(jié)的波數(shù)選取規(guī)律也適用于空間域非均勻采樣,且基于CPU-GPU加速方案的計(jì)算時間比單核CPU計(jì)算時間快26倍.

5 結(jié)論

本文實(shí)現(xiàn)了基于二次插值形函數(shù)任意傅里葉變換的空間-波數(shù)域重力異常數(shù)值模擬算法,并采用CPU-GPU進(jìn)行加速,進(jìn)一步提高了現(xiàn)有空間-波數(shù)域算法的優(yōu)勢.得到以下結(jié)論:

(1) 任意傅里葉變換將二維傅里葉變換轉(zhuǎn)化為兩個一維傅里葉變換,一維傅里葉變換離散為多個單元積分之和,單元內(nèi)原函數(shù)采用二次形函數(shù)擬合,得到單元積分的解析解,新方法采樣靈活、積分精度高、計(jì)算速度快和傅里葉變換的截?cái)嘈?yīng)小.

(2) 利用棱柱體模型,將數(shù)值解與解析解對比,異常體內(nèi)部和外部場的誤差均滿足精度要求,驗(yàn)證了波數(shù)選取規(guī)律的有效性和方法的正確性;利用連續(xù)介質(zhì)模型,與高斯傅里葉變換NG=4的結(jié)果對比,表明在計(jì)算精度相近的情況下,本文算法所需波數(shù)少、耗時短、占用內(nèi)存少.

(3) 基于NVIDIA TITAN顯卡,對比不同GPU線程層次加速任意傅里葉變換的計(jì)算效率,結(jié)果表明二維層次加速效果最好;對比多核CPU、GPU分別加速五對角方程組求解和任意傅里葉變換的計(jì)算效率,隨著空間域和波數(shù)域計(jì)算節(jié)點(diǎn)數(shù)增加,CPU-GPU的并行效果越好,最高加速比可達(dá)48.

(4) 將新方法和加速方案應(yīng)用于長白山火山區(qū)域地形重力異常正演,實(shí)驗(yàn)結(jié)果表明該算法計(jì)算精度高、效率高,適用于任意復(fù)雜地形、復(fù)雜地質(zhì)構(gòu)造,實(shí)用性強(qiáng).

基于任意傅里葉變換的空間-波數(shù)域算法及其CPU-GPU的并行加速方案同樣適用于弱磁、強(qiáng)磁、直流電和電磁場的數(shù)值模擬,下一步將研究新方法應(yīng)用于電磁法正演計(jì)算中.