任意起伏地形下重力異常三維正演及并行計算

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

1 中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室,長沙 410083 2 中南大學有色資源與地質災害探查湖南省重點實驗室,長沙 410083 3 中南大學地球科學與信息物理學院,長沙 410083 4 西南石油大學地球科學與技術學院,成都 610500

0 引言

重力勘探廣泛用于油氣勘探、金屬礦勘查、地球內部構造、區域地質構造以及環境地球物理等研究領域.研究針對大規模復雜條件下重力異常場高效、高精度三維數值模擬對重力異常反演以及人機交互解釋有著重要作用.重力數值模擬方法按大類主要可分為空間域和頻率域方法.

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

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

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

1 理論方法

1.1 空間-波數域重力微分法正演

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

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

(1)

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

(2)

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

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

1.2 任意傅里葉變換

1.2.1 任意傅里葉變換理論

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

(3)

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

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

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

(5)

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

(6)

(7)

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

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

(8)

式中ej表示第j單元,令第j個單元的三個節點坐標為(xj1,xj2,xj3),xj2為單元ej中間節點,滿足關系式xj1+xj3=2xj2.單元內采用二次插值的形函數擬合函數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個單元上三個節點的傅里葉變換積分系數,該積分系數可求得解析解.式(10)可寫為

(12)

式中,

(13)

同理,式(7)可寫為

(14)

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

(15)

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

1.2.2 采樣規則

任意采樣點數FFT算法提出了離散頻率取值方式(陳龍偉等,2016),標準FFT存在截斷效應.擴邊FFT能減小這種誤差,其本質是擴邊之后波數域基頻減小,提高了計算精度;Gauss-FFT作為一種結合高斯點數的偏移采樣傅里葉變換算法,截斷效應小,計算精度高(吳樂園,2014),其原因是Gauss-FFT基于FFT采樣定理,在基頻單元內將采樣點加密并偏移,加密點數為高斯點個數,偏移距離為高斯偏移距離,能有效提高積分計算精度.這兩種FFT算法在計算離散頻譜時均利用了采樣定理,要求空間域和頻率域均勻采樣,無法實現非均勻采樣.此外,為滿足精度要求,擴邊FFT和4個點的Gauss-FFT需要更多的波數,降低了計算效率.

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

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

(16)

基頻Δk作為基本單位:

(17)

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

在求得正波數采樣范圍后,可根據頻譜能量分布總結出波數選取的規律.本文總結波數選取規律的具體準則為:頻譜能量大且變化劇烈的區間加密采樣,頻譜能量小且變化緩慢的區間稀疏采樣.

1.3 CPU-GPU并行理論方法

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

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

圖1 OpenMP并行機制

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

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

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

圖2 x方向一維積分計算過程

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

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

圖3 二維線程層次結構

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

2 數值算例

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

2.1 正確性驗證

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

(18)

設計棱柱體模型,計算區域范圍20 km(x)×20 km(y)×7.5 km(z),異常體水平方向上位于模擬區域中心,頂面埋深為0.5 km,大小為8 km(x)×8 km(y)×2 km(z),異常體密度為2000 kg·m-3.剖分節點個數為251(x)×251(y)×151(z),三個方向上均勻剖分,水平方向網格間距均為80 m,垂直方向為50 m.

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

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

圖4 重力異常場數值解與解析解等值線圖

2.2 計算精度與計算效率

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

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

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

(19)

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

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

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

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

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

3 CPU-GPU并行加速

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

3.1 GPU線程層次加速對比

一維層次、二維層次以及三維層次GPU結構加速任意傅里葉變換,同時使用culbas庫直接加速任意傅里葉變換.保持z方向節點數為151,改變水平方向計算節點個數,波數個數與空間域節點數保持一致,對比四種GPU方法加速任意傅里葉變換的計算時間以及GPU顯存-CPU內存傳輸時間,結果如表2所示.

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

從表中可以看出,不同計算節點個數下,一維、二維、三維層次的計算時間相近,但是二維層次的顯存-內存傳輸耗時最少,二維層次下的計算總時間少于一維和三維層次的計算總時間,同時一維、二維、三維的計算總時間優于CUBLAS庫.因此,任意傅里葉變換GPU加速方案采用二維層次.

3.2 CPU、GPU加速方案對比

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

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

從表中可以看出,五對角方程組求解過程中,CPU加速效果好于GPU加速,這是因為每個波數下的五對角方程組求解過程復雜,并不是簡單的點對點簡單運算,隨著計算節點數增加,加速比最高可達13.6,因此,采用CPU加速五對角方程組求解.

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

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

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

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

4 實際地形應用

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

圖6 長白山火山區域地形圖

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

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

圖7 實際地形重力異常場數值解以及解析解結果

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

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

5 結論

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

(1) 任意傅里葉變換將二維傅里葉變換轉化為兩個一維傅里葉變換,一維傅里葉變換離散為多個單元積分之和,單元內原函數采用二次形函數擬合,得到單元積分的解析解,新方法采樣靈活、積分精度高、計算速度快和傅里葉變換的截斷效應小.

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

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

(4) 將新方法和加速方案應用于長白山火山區域地形重力異常正演,實驗結果表明該算法計算精度高、效率高,適用于任意復雜地形、復雜地質構造,實用性強.

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