簡化的混合域全波形反演方法及GPU加速

桂 生,劉 洪,李 飛

(1.中國科學(xué)院地質(zhì)與地球物理研究所,北京100029;2.中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室,北京100029)

簡化的混合域全波形反演方法及GPU加速

桂 生1,2,劉 洪1,2,李 飛1,2

(1.中國科學(xué)院地質(zhì)與地球物理研究所,北京100029;2.中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室,北京100029)

全波形反演(FWI)方法綜合利用疊前地震波場的動(dòng)力學(xué)和運(yùn)動(dòng)學(xué)信息,能夠高精度地重建地下介質(zhì)模型參數(shù)場,但巨大的計(jì)算量一直是制約其發(fā)展的一個(gè)重要因素。GPU組成的高性能計(jì)算集群為提高全波形反演計(jì)算效率提供了重要的硬件基礎(chǔ)?；贕PU平臺(tái),采用簡化的混合域全波形反演算法實(shí)現(xiàn)了更快速的三維全波形反演計(jì)算。首先簡單介紹了GPU加速技術(shù)應(yīng)用于簡化的混合域全波形反演時(shí)的一些優(yōu)化技巧,包括線程調(diào)度、GPU之間數(shù)據(jù)傳輸以及共享內(nèi)存的使用等,然后通過多GPU全波形反演測試了簡化的混合域全波形反演的效果,證明了GPU加速技術(shù)能夠有效地提高全波形反演的計(jì)算效率,相比CPU具有十幾倍的加速比。

GPU技術(shù);高性能計(jì)算;GPU技術(shù)優(yōu)化;全波形反演

全波形反演最初由LAILLY[1]和TARANTOLA[2]等人提出,他們以波動(dòng)方程為基礎(chǔ),建立模擬數(shù)據(jù)和觀測數(shù)據(jù)差值的目標(biāo)函數(shù),通過使其最小化來獲得模型的更新量,其模型更新量由梯度和搜索步長決定,梯度可由入射波場和殘差反傳的波場計(jì)算得到。目前已經(jīng)發(fā)展出了多種全波形反演(FWI)方法,按照不同求解波動(dòng)方程的計(jì)算域,大致可分為時(shí)間域反演[3-8]、頻率域反演[9-11]和Laplace域反演[12-13]。按照求解方式不同又可分為局部梯度類迭代反演[2-3,14-15]和全局優(yōu)化搜索反演[16-19]等。

無論是時(shí)間域、頻率域還是Laplace域,巨大的計(jì)算量是FWI瓶頸之一。近年來高性能計(jì)算技術(shù)的發(fā)展使這一問題得以緩解。高性能計(jì)算的主要發(fā)展思路是并行計(jì)算。相對(duì)于CPU而言,GPU的多線程、多處理器的并行結(jié)構(gòu)更能夠滿足高性能計(jì)算的發(fā)展需要[20],而且在相同計(jì)算能力的情況下,GPU具有功耗小、價(jià)格低等特點(diǎn)。在GPU中,用來進(jìn)行數(shù)據(jù)處理的晶體管占據(jù)大部分空間,只有少量晶體管用來進(jìn)行指令流控制和數(shù)據(jù)緩存[21]。利用GPU進(jìn)行密集型科學(xué)計(jì)算是一種潛力巨大且性價(jià)比高的解決方案[22]。

由于GPU具有對(duì)密集型計(jì)算任務(wù)加速的優(yōu)勢,因此GPU技術(shù)被廣泛應(yīng)用于正演模擬、疊前偏移、逆時(shí)偏移、多次波壓制以及全波形反演等地震數(shù)據(jù)處理。MICIKEVICIUS[23]用GPU實(shí)現(xiàn)了三維有限差分法計(jì)算,并提出了單通法和雙通法兩種優(yōu)化策略,也是目前常用的程序優(yōu)化方法。WANG等[24]研究了基于隨機(jī)邊界條件的GPU加速時(shí)間域全波形反演算法,LUO等[25]實(shí)現(xiàn)了多GPU的兩級(jí)并行的全波形反演,MAO等[26]利用GPU加速實(shí)現(xiàn)了二維多尺度時(shí)間域全波形反演。KIM等[27-28]利用GPU加速實(shí)現(xiàn)了三維聲波混合域全波形反演的快速計(jì)算。SHIN等[29]實(shí)現(xiàn)了單GPU的三維Laplace域的全波形反演,GOKHBERG等[30]用譜元法在CPU/GPU異構(gòu)集群實(shí)現(xiàn)了全波形反演,張猛[31]對(duì)CPU/GPU平臺(tái)上全波形反演的各種實(shí)用問題進(jìn)行了分析。LIU等[32]實(shí)現(xiàn)了三維混合域全波形反演,但其寬度較窄。YANG等[33]在GPU平臺(tái)上實(shí)現(xiàn)了二維時(shí)間域的全波形反演。由于GPU的顯存限制,以往的研究多是針對(duì)二維或者小規(guī)模的三維全波形反演,為克服這一限制,本文采用多GPU動(dòng)態(tài)加載并行方案,針對(duì)不同大小的數(shù)據(jù)采用不同的并行方案,當(dāng)計(jì)算數(shù)據(jù)過大時(shí),采用區(qū)域分解的方法,從而可以實(shí)現(xiàn)大規(guī)模三維全波形反演;采用簡化的混合域全波形反演,能夠更快速的計(jì)算從以往的單GPU加速到現(xiàn)在的多GPU加速;最后利用Overthrust模型測試了本文方法的可行性和有效性。

1 GPU架構(gòu)及程序優(yōu)化

1.1 GPU架構(gòu)

GPU體系架構(gòu)如圖1所示。

圖1 GPU運(yùn)算單元結(jié)構(gòu)[21]

首先從硬件架構(gòu)上對(duì)GPU進(jìn)行并行層次分析。GPU主要由兩部分組成,一部分是流處理器陣列,另一部分是存儲(chǔ)器系統(tǒng),兩部分連接在同一片板卡上。如圖1所示,其中多個(gè)線程處理器群(Thread Processing Cluster,TPC)組成了一個(gè)流處理器陣列,其中每個(gè)TPC中有2～3個(gè)流多處理器(Streaming Multiprocessor,SM),每個(gè)流多處理器中包含8個(gè)流處理器(Streaming Processor,SP)[34]。在GPU中最基本單元是SP,SP是由算術(shù)邏輯運(yùn)算器、加法器、乘法器以及寄存控制器組成。SP可以在一個(gè)時(shí)鐘周期內(nèi)完成一次邏輯運(yùn)算,或一次32位浮點(diǎn)數(shù)加法或乘法運(yùn)算等。從硬件架構(gòu)上看,一個(gè)流多處理器相當(dāng)于一個(gè)8路單指令多數(shù)據(jù)流(Single Instruction Multiple Data)處理器,稱之為單指令流多線程(Single Instruction Multiple Thread,SIMT)[21]。

圖2為統(tǒng)一計(jì)算設(shè)備架構(gòu)(Compute Unified Device Architecture,CUDA)編程模式示意圖,可以看出,核函數(shù)將線程劃分為網(wǎng)格(Grid)形式,每個(gè)網(wǎng)格又被分成無數(shù)個(gè)線程塊(Block),一定數(shù)量的線程(Thread)組成(一般不超過1024個(gè)線程)一個(gè)線程塊。ThreadIdx和BlockIdx是用來標(biāo)明線程號(hào)和線程塊號(hào)的CUDA內(nèi)置變量。在核函數(shù)中每個(gè)線程都有一個(gè)唯一的偏移量,偏移量可以通過線程號(hào)、線程塊號(hào)以及線程維度信息計(jì)算得到。在一維情況下偏移量等于線程號(hào)即為ThreadIdx.x;在二維情況下且Block大小為(Bx,By)的時(shí),計(jì)算線程偏移量公式為ThreadIdx.x+ThreadIdx.y×Bx;而在Block大小為(Bx,By,Bz)的三維情況下,線程偏移量計(jì)算公式為ThreadIdx.x+ThreadIdx.y×Bx+ThreadIdz.z×By×Bz。在同一個(gè)網(wǎng)格下的不同線程塊中線程不能進(jìn)行通訊,意味著不同線程塊中線程無法進(jìn)行數(shù)據(jù)傳輸,它們完全并行且獨(dú)立執(zhí)行。

與CPU相比,GPU在計(jì)算能力和存儲(chǔ)器帶寬上有顯著優(yōu)勢,面對(duì)相同量計(jì)算任務(wù),GPU所需要的硬件成本和計(jì)算功耗都低于CPU,因此,GPU是大規(guī)模高性能計(jì)算的明智選擇。

CUDA的出現(xiàn)使得基于GPU的通用計(jì)算(GeneralPurpose GPU,GPGPU)得到了高速發(fā)展?，F(xiàn)階段GPU的浮點(diǎn)運(yùn)算能力已經(jīng)達(dá)到每秒數(shù)十萬億次級(jí)別,是X86處理器計(jì)算能力的數(shù)十倍。但是要發(fā)揮GPU巨大的計(jì)算優(yōu)勢,就需要對(duì)程序進(jìn)行一定的優(yōu)化,未經(jīng)優(yōu)化的程序反而會(huì)降低程序的計(jì)算效率。

圖2 CUDA編程模式[21]

圖3對(duì)比了CPU和GPU的微觀硬件架構(gòu),可以看出CPU和GPU計(jì)算能力差異巨大的原因。圖中ALU(Arithmetic Logic Unit,算術(shù)邏輯單元)是計(jì)算單元,CPU只具有很少的ALU,而GPU則具有大量的ALU;CPU具有大量的數(shù)據(jù)緩存和邏輯控制單元,而GPU的數(shù)據(jù)緩存和邏輯控制單元較少,這就導(dǎo)致GPU的計(jì)算能力要遠(yuǎn)遠(yuǎn)高于CPU的計(jì)算能力,但數(shù)據(jù)讀取以及邏輯運(yùn)算能力較差,因此GPU非常適合處理計(jì)算密集型和并行度較高的計(jì)算問題。在FWI中,波場的正傳、反傳以及梯度計(jì)算的計(jì)算量所占比例最大。本文的波場正傳和反傳均采用時(shí)間域有限差分算法,有限差分算法中,每個(gè)點(diǎn)都是離散的,其計(jì)算是獨(dú)立的,可以同時(shí)進(jìn)行,因此GPU非常適合這種細(xì)粒度并行計(jì)算。根據(jù)計(jì)算階數(shù)的不同需要不同個(gè)數(shù)的鄰近網(wǎng)格點(diǎn)上一個(gè)時(shí)刻的值,因此在計(jì)算時(shí)需要采用共享內(nèi)存方式存儲(chǔ)鄰近網(wǎng)格點(diǎn)的數(shù)值,隱藏訪存延遲。GPU高性能計(jì)算用于地震波有限差分?jǐn)?shù)值模擬能大大提高計(jì)算速度,并且具有更小的能耗。梯度的計(jì)算也每個(gè)點(diǎn)獨(dú)立計(jì)算,不涉及和其它網(wǎng)格點(diǎn)的數(shù)據(jù)交換,屬于細(xì)粒度的并行,因此采用GPU同樣能提高梯度的計(jì)算效率。如何根據(jù)FWI的求解區(qū)域分配計(jì)算任務(wù)、如何傳輸數(shù)據(jù)以及如何發(fā)揮GPU的最優(yōu)計(jì)算性能都需要根據(jù)硬件信息以及算法進(jìn)行優(yōu)化,其中包括線程塊的大小、多GPU之間的數(shù)據(jù)傳輸、共享內(nèi)存的使用等。

1.2 線程調(diào)度

在GPU中最基本的執(zhí)行單元是線程,其次是線程塊,每個(gè)線程塊的大小影響著計(jì)算性能,線程塊中的每一個(gè)線程被發(fā)送到一個(gè)SP上。每個(gè)SM上有8個(gè)SP,每個(gè)SM最多可以同時(shí)執(zhí)行8個(gè)線程塊。SM最多能執(zhí)行的線程數(shù)量和線程能使用的資源相互矛盾,一個(gè)SM里同時(shí)執(zhí)行的線程越多,就越能隱藏線程訪存延遲,這樣每個(gè)線程能使用的資源相對(duì)就越少,因此需要選取適當(dāng)?shù)木€程塊大小,使兩者平衡,從而得到最優(yōu)的執(zhí)行效率。不同的硬件擁有的計(jì)算資源不一樣,為充分利用硬件資源需要進(jìn)行測試分析。

本文采用12階規(guī)則網(wǎng)格有限差分法進(jìn)行波場正傳和反傳計(jì)算。為此,我們使用NVIDIA Tesla K20m GPU對(duì)此算法進(jìn)行了測試分析。三維模型大小256×256×256,網(wǎng)格間距10m。

圖3 CPU的微觀架構(gòu)(a)和GPU的微觀架構(gòu)(b)對(duì)比

圖4中橫坐標(biāo)表示線程塊大小,縱坐標(biāo)表示某一個(gè)時(shí)刻的計(jì)算時(shí)間(ms),從圖中可以看出,剛開始隨著線程塊大小的增加,計(jì)算時(shí)間逐漸減少,這是因?yàn)槊總€(gè)SM上執(zhí)行的線程數(shù)量隨著線程塊大小的增加而增加;當(dāng)線程塊中的線程數(shù)量增加到32×32時(shí),計(jì)算時(shí)間反而急劇增加,這是因?yàn)檫^大的線程塊增加了寄存器的使用數(shù)量,使得能同時(shí)計(jì)算的線程塊數(shù)量減少,同時(shí)執(zhí)行的線程數(shù)量也減少所致。綜合考慮線程塊大小及其使用資源數(shù)量,我們進(jìn)行12階交錯(cuò)網(wǎng)格有限差分計(jì)算時(shí)采用的線程塊大小為32×16。

圖4 不同線程塊計(jì)算時(shí)間對(duì)比

1.3 多GPU優(yōu)化

當(dāng)計(jì)算的數(shù)據(jù)體十分巨大時(shí),單個(gè)GPU顯存無法滿足計(jì)算需求,需采用多GPU并行方案將數(shù)據(jù)進(jìn)行區(qū)域分解,每個(gè)GPU運(yùn)算一個(gè)區(qū)域。這樣不僅能克服單個(gè)GPU顯存較小的問題,同時(shí)由于多個(gè)GPU并行計(jì)算,提高了計(jì)算效率,減少了計(jì)算耗時(shí)。我們以規(guī)則網(wǎng)格12階有限差分進(jìn)行測試,網(wǎng)格大小256×256×256,網(wǎng)格間距10m。

圖5中橫坐標(biāo)代表GPU數(shù)量,縱坐標(biāo)代表計(jì)算一個(gè)時(shí)刻波場的計(jì)算耗時(shí),從圖中可以看出,隨著GPU的增多,計(jì)算時(shí)間成線性減少,理想情況下應(yīng)該成倍數(shù)減少,但是由于區(qū)域分解后需要對(duì)分解后的數(shù)據(jù)進(jìn)行擴(kuò)邊用來交換數(shù)據(jù),保證數(shù)據(jù)的一致性;同時(shí)由于各個(gè)GPU之間也需要進(jìn)行數(shù)據(jù)交換,增加了計(jì)算耗時(shí),因此無法達(dá)到理想的情況。

圖5 多GPU計(jì)算性能對(duì)比

利用多GPU進(jìn)行有限差分計(jì)算時(shí),計(jì)算邊界處需要進(jìn)行數(shù)據(jù)交換,數(shù)據(jù)交換分為兩種：一種是數(shù)據(jù)傳輸?shù)紺PU端,再從CPU端傳輸?shù)搅硪粔KGPU端;第二種是2塊GPU在一條PCIE線上時(shí),可以直接利用對(duì)等網(wǎng)絡(luò)(peer-to-peer,P2P)技術(shù)進(jìn)行數(shù)據(jù)交換(圖6),不需要經(jīng)過CPU,大大減少了傳輸時(shí)間。利用P2P技術(shù)傳輸4GB大小數(shù)據(jù),所耗時(shí)間為0.11s,不使用P2P技術(shù)時(shí)的傳輸時(shí)間為0.28s,可以看出P2P技術(shù)大大提高了傳輸速度,特別是針對(duì)多GPU全波形反演時(shí)頻繁的數(shù)據(jù)交換情況。

1.4 共享內(nèi)存優(yōu)化

在使用GPU時(shí),需要將數(shù)據(jù)傳輸?shù)紾PU上,有兩種方法提高傳輸效率,一是使用鎖頁內(nèi)存(其存放的數(shù)據(jù)傳輸?shù)紾PU不需要CPU干涉)提高傳輸速度;第二種方法是改進(jìn)硬件條件,如最新Nvlink技術(shù)的傳輸速度達(dá)到80GB/s,是現(xiàn)有PCIE3.0技術(shù)傳輸速度的4倍以上。本文主要是在程序上進(jìn)行方法的改進(jìn),因此采用鎖頁內(nèi)存方式提高數(shù)據(jù)的傳輸效率。經(jīng)過測試,使用鎖頁內(nèi)存?zhèn)鬏?GB大小的數(shù)據(jù)用時(shí)0.13s,而不使用鎖頁內(nèi)存?zhèn)鬏斖瑯哟笮〉臄?shù)據(jù),則需要0.17s。

前文談到,在進(jìn)行12階規(guī)則網(wǎng)格的有限差分計(jì)算時(shí),需要讀取網(wǎng)格點(diǎn)鄰近點(diǎn)上一時(shí)刻的波場,而且需要重復(fù)讀取。為了避免每次計(jì)算都需要單獨(dú)讀取波場值,我們采用共享內(nèi)存的方式加以解決。計(jì)算時(shí),直接在開辟的共享內(nèi)存區(qū)域重復(fù)進(jìn)行高速訪存,這樣很大程度上提高了有限差分法的計(jì)算效率。圖7 為共享內(nèi)存示意圖,其中藍(lán)色部分是需要計(jì)算的點(diǎn),黃色部分是邊界上需要多余讀取的數(shù)據(jù)點(diǎn)。由圖可知共享內(nèi)存實(shí)際存儲(chǔ)的數(shù)據(jù)要比計(jì)算數(shù)據(jù)大,以12階規(guī)則網(wǎng)格為例,實(shí)際存儲(chǔ)的數(shù)據(jù)大小為12nx+12ny+ny×nx,其中nx為計(jì)算區(qū)域x方向的網(wǎng)格點(diǎn)數(shù),ny為y方向的網(wǎng)格點(diǎn)數(shù)。測試發(fā)現(xiàn),使用共享內(nèi)存比未使用共享內(nèi)存的程序加速比高2.6倍。

圖6 P2P技術(shù)原理

圖7 共享內(nèi)存示意

2 簡化的混合域FWI方法GPU實(shí)現(xiàn)

對(duì)于三維全波形反演,頻率域FWI的正傳和反傳均需要存儲(chǔ)巨大的阻抗矩陣,而時(shí)間域FWI梯度計(jì)算太復(fù)雜[35]。因此,本文基于混合域思想對(duì)其進(jìn)行簡化。簡化的方法是將波場的正傳和反傳都放到時(shí)間域求取,利用離散傅里葉變換(DFT)抽取相應(yīng)波場的頻率成分,在頻率域求取殘差,再用離散傅里葉反變換(IDFT)將頻率域殘差返回時(shí)間域進(jìn)行反傳;然后將反傳時(shí)間域波場再次用DFT轉(zhuǎn)換到頻率域波場,在頻率域求取梯度場。為了方便GPU實(shí)現(xiàn),簡化的頻率域方法直接在時(shí)間域進(jìn)行殘差計(jì)算,然后將殘差波場反傳并用DFT轉(zhuǎn)換到頻率域,在頻率域求取梯度。簡化的混合域全波形反演流程如圖8所示。

圖8 簡化的混合域全波形反演流程

圖8中藍(lán)色部分是利用GPU實(shí)現(xiàn)的,需要注意的是對(duì)于單個(gè)節(jié)點(diǎn)來說數(shù)據(jù)加載是動(dòng)態(tài)的,根據(jù)計(jì)算數(shù)據(jù)的大小和現(xiàn)有計(jì)算資源動(dòng)態(tài)選取并行方案。在炮循環(huán)之前,程序會(huì)讀取各個(gè)節(jié)點(diǎn)的硬件信息,主要是GPU的個(gè)數(shù)和單個(gè)GPU的顯存大小,當(dāng)需要計(jì)算的數(shù)據(jù)小于GPU的顯存時(shí),采用一個(gè)GPU計(jì)算一炮地震數(shù)據(jù),如圖9所示。

當(dāng)需要計(jì)算的數(shù)據(jù)大于單個(gè)GPU的顯存時(shí),系統(tǒng)會(huì)根據(jù)實(shí)際需要的顯存數(shù)量以及GPU個(gè)數(shù)動(dòng)態(tài)地將數(shù)據(jù)進(jìn)行區(qū)域分解,以獲得最佳的計(jì)算性能,這時(shí)每個(gè)GPU將計(jì)算一小塊反演區(qū)域的數(shù)據(jù),如圖10a 所示(本文指定z方向劃分)。

同時(shí)由于整個(gè)反演區(qū)域分解成數(shù)個(gè)小的反演區(qū)域,那么區(qū)域邊界處就需要進(jìn)行數(shù)據(jù)交換,以保證數(shù)據(jù)的一致性,這樣在小數(shù)據(jù)的邊界處增加了交換區(qū)域,見圖10b。

圖9 單GPU未進(jìn)行區(qū)域分解

圖10 多GPU區(qū)域分解(a)和交換區(qū)域增加(b)

GPU實(shí)現(xiàn)程序偽代碼如下：

讀取各個(gè)節(jié)點(diǎn)的硬件信息;

if(小于單個(gè)GPU內(nèi)存)

{for(炮循環(huán))

根據(jù)節(jié)點(diǎn)數(shù)以及GPU個(gè)數(shù)信息MPI向各節(jié)點(diǎn)發(fā)送多少炮數(shù)據(jù);

for(單GPU需要計(jì)算的炮數(shù)循環(huán))

為每個(gè)GPU動(dòng)態(tài)發(fā)送炮數(shù)據(jù)

for(時(shí)間循環(huán))

每個(gè)GPU同時(shí)計(jì)算一炮數(shù)據(jù)

}

else(大于單個(gè)GPU內(nèi)存)

{for(炮循環(huán)){

根據(jù)節(jié)點(diǎn)數(shù)以及GPU個(gè)數(shù)信息MPI向各節(jié)點(diǎn)發(fā)送多少炮數(shù)據(jù);

并將每個(gè)節(jié)點(diǎn)上單炮數(shù)據(jù)進(jìn)行區(qū)域分解,載入對(duì)應(yīng)的GPU

for(時(shí)間循環(huán))

每個(gè)GPU計(jì)算一小塊炮數(shù)據(jù)

}}

3 模型測試

我們利用SEG的Overthrust模型(圖11)對(duì)本文方法進(jìn)行了測試,網(wǎng)格大小801×801×187,網(wǎng)格間距均為10m。地表放炮,共100炮,炮間距40m,x方向和y方向各10炮,第一炮位置位于x=3800m,y=380m處。使用主頻為35Hz的零相位雷克子波,采集系統(tǒng)是地表全接收。時(shí)間采樣率2ms,采樣點(diǎn)數(shù)為2000。共正演地震數(shù)據(jù)478GB。反演頻率從1.5Hz開始,每次增加1.5Hz,2個(gè)頻率組成一個(gè)頻率組,最大頻率60Hz。由于計(jì)算數(shù)據(jù)太大,這里動(dòng)態(tài)分配為4個(gè)GPU進(jìn)行處理。

根據(jù)節(jié)點(diǎn)信息,程序首先將初始模型數(shù)據(jù)分解成4個(gè)小數(shù)據(jù)體,每個(gè)數(shù)據(jù)體大小如圖12a所示,每個(gè)GPU負(fù)責(zé)一個(gè)小數(shù)據(jù)體的計(jì)算,其反演結(jié)果見圖12b,可以看出反演結(jié)果和真實(shí)模型(圖11)能相互對(duì)應(yīng),在反演結(jié)果中能清楚地識(shí)別斷層之間褶皺成因的背斜。圖13a和圖13b分別為初始模型和反演結(jié)果在x=4000m,y=4000m,z=400m處的切片。由圖13a

圖11 Overthrust真實(shí)模型

可以看出,河道信息比較模糊,無法識(shí)別;由圖13b可看出,反演結(jié)果很好恢復(fù)出了河道信息,能有效識(shí)別出逆斷層,雖然局部信息有差別,但大體都恢復(fù)得很好。

采用4節(jié)點(diǎn)8核CPU與同樣數(shù)量的GPU進(jìn)行對(duì)比。在相同任務(wù)量(一次迭代計(jì)算)的情況下,GPU耗時(shí)996s,而CPU耗時(shí)14458s,同樣數(shù)量的GPU相對(duì)于CPU其程序加速比達(dá)到了14.5倍,可見GPU在FWI的快速計(jì)算中具有明顯優(yōu)勢。

圖12 使用4塊GPU時(shí),初始模型數(shù)據(jù)分解成的小數(shù)據(jù)體(a)及小數(shù)據(jù)體反演結(jié)果(b)(z方向上分解成的4個(gè)小數(shù)據(jù)體之一)

圖13 使用4塊GPU時(shí),初始模型(a)和反演結(jié)果(b)在x,y,z方向上的切片及其連片顯示

4 結(jié)論

簡化混合域方法直接利用時(shí)間域數(shù)據(jù)來實(shí)現(xiàn)多尺度反演,與以前的混合域方法相比,流程每次迭代少了一次DFT和IDFT,并且可以獲得良好的反演效果。

GPU顯存有限,所以針對(duì)大規(guī)模三維混合域波形反演,特別是疊前大數(shù)據(jù)的計(jì)算,多GPU成為必然選擇。在使用多GPU計(jì)算時(shí)需要對(duì)程序進(jìn)一步優(yōu)化,才能發(fā)揮GPU的最優(yōu)性能,本文提出的多GPU動(dòng)態(tài)并行方案具有良好的全波形反演結(jié)果并具有較好的加速性能,為實(shí)現(xiàn)大規(guī)模FWI的實(shí)際應(yīng)用提供一種可選擇的計(jì)算方法和計(jì)算框架。

[1] LAILLY P.The seismic inverse problem as a sequence of before stack migrations[C]∥BEDNAR J B.Conference on inverse scattering:theory and application.Philadelphia:SIAM,1983:206-220

[2] TARANTOLA A,NERCESSIAN A,GAUTHIER O.Nonlinear inversion of seismic reflection data[J].Geophysics,1984,49(2):645-649

[3] BLEISTEIN N,COHEN J K,HAGIN F G.Two and one-half dimensional Born inversion with an arbitrary reference[J].Geophysics,1987,52(1):26-36

[4] CRASE E,PICA A,NOBLE M,et al.Robust elastic nonlinear waveform inversion:application to real data[J].Geophysics,1990,55(5):527-538

[5] KOLB P,COLLINO F,LAILLY P.Pre-stack inversion of a 1-D medium[J].Proceedings of the IEEE,1986,74(3):498-508

[6] MARFURT K J.Accuracy of finite difference and finite element modeling of the scalar and elastic wave equation[J].Geophysics,1984,49(5):533-549

[7] PICA A,DIET J P,TARANTOLA A.Nonlinear inversion of seismic reflection data in a laterally invariant medium[J].International Journal of Science Education,1990,55(3):284-292

[8] WANG B L,GAO J H.Fast full waveform inversion of multi-shot seismic data[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010：4453-4459

[9] PRATT G,SHIN C,HICKS.Gauss-Newton and full Newton methods in frequency space seismic waveform inversion[J].Geophysical Journal International,1998,133(2):341-362

[10] SIRGUE L,PRATT R G.Efficient waveform inversion and imaging:a strategy for selecting temporal frequencies[J].Geophysics,2004,69(1):231-248

[11] SONG Z M,WILLIAMSON P R,PRATT R G.Frequency-domain acoustic-wave modeling and inversion of crosshole data:part Ⅱ:inversion method,synthetic experiments and real-data results[J].Geophysics,1995,60(S1):796-809

[12] CHANGSOO S,HO C Y.Waveform inversion in the Laplace domain[J].Geophysical Journal International,2008,173(3):922-931

[13] HA W,SHIN C.Laplace-domain full-waveform inversion of seismic data lacking low-frequency information[J].Geophysics,2012,77(5):R199-R206

[14] GAUTHIER O,VIRIEUX J,TARANTOLA A.Two-dimensional nonlinear inversion of seismic waveforms:numerical results[J].Geophysics,1986,51(7):1387-1403

[15] MORA P.Nonlinear two-dimensional elastic inversion of multioffset seismic data[J].Geophysics,1987,52(9):1211-1228

[16] SEN M K,STOFFA P L.Global optimization methods in geophysical inversion[J].Advances in Exploration Geophysics,1995,4：269-277

[17] SAMBRIDGE M,DRIJKONINGEN G.Genetic algorithms in seismic waveform inversion[J].Geophysical Journal International,1992,109(2):323-342

[18] JIN S,MADARIAGA R.Nonlinear velocity inversion by a two-step Monte Carlo method[J].Geophysics,1994,59(4):577-590

[19] SEN M K,STOFFA P L.Nonlinear one-dimensional seismic waveform inversion using simulated annealing[J].Geophysics,1991,56(10):1624-1638

[20] 趙改善.地球物理高性能計(jì)算的新選擇:GPU 計(jì)算技術(shù)[J].勘探地球物理進(jìn)展,2007,30(5):399-404 ZHAO G S.New alternative to geophysical high performance computing:GPU computing [J].Progress in Exploration Geophysics,2007,30(5):399-404

[21] KIRK D.NVIDIA CUDA software and GPU parallel computing architecture[OB/EL].http:∥kr.nvidia.com/content/cudazone/download/showcase/kr/Tutorial-DKIRK.pdf,2007:103-104

[22] 李博,劉國峰,劉洪.地震疊前時(shí)間偏移的一種圖形處理器提速實(shí)現(xiàn)方法[J].地球物理學(xué)報(bào),2009,52(1):245-252 LI B,LIU G F,LIU H.A method of using GPU to accelerate seismic pre-stack time migration[J].Chinese Journal of Geophysics,2009,52(1):245-252

[23] MICIKEVICIUS P.3D finite difference computation on GPUs using CUDA[J].The Workshop on General Purpose Processing on Graphics Processing Units,2009:79-84

[24] WANG B L,GAO J H,ZHANG H L,et al.CUDA-based acceleration of full waveform inversion on GPU[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2528-2533

[25] LUO J R,GAO J H,WANG B L.Multi-GPU based two-level acceleration of full waveform inversion[J].Journal of Seismic Exploration,2012,21(4):377-394

[26] MAO J,WU R S,WANG B L.Multiscale full waveform inversion using GPU[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:7-13

[27] KIM Y,SHIN C,CALANDRA H.3D hybrid waveform inversion with GPU devices[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-5

[28] KIM Y,SHIN C,CALANDRA H,et al.An algorithm for 3D acoustic time-Laplace-Fourier-domain hybrid full waveform inversion[J].Geophysics,2013,78(4):R151-R166

[29] SHIN J,HA W,JUN H,et al.3D Laplace-domain full waveform inversion using a single GPU card[J].Computers & Geosciences,2014,67(1):1-13

[30] GOKHBERG A,FICHTNER A.Full-waveform inversion on heterogeneous HPC systems [J].Computers & Geosciences,2016,89(3):260-268

[31] 張猛,王華忠,任浩然,等.基于CPU/GPU異構(gòu)平臺(tái)的全波形反演及其實(shí)用化分析[J].石油物探,2014,53(4):461-467 ZHANG M,WANG H Z,REN H R,et al,Full waveform inversion on the CPU/GPU heterogeneous platform and its application analysis[J].Geophysical Prospecting for Petroleum,2014,53(4):461-467

[32] LIU L,DING R,LIU H,et al.3D hybrid-domain full waveform inversion on GPU[J].Computers & Geosciences,2015,83(1):27-36

[33] YANG P L,GAO J H,WANG B L.A graphics processing unit implementation of time-domain full-waveform inversion[J].Geophysics,2015,80(3):F31-F39

[34] 王澤寰,王鵬.GPU并行計(jì)算編程技術(shù)介紹[J].科研信息化技術(shù)與應(yīng)用,2013,4(1):81-87 WANG Z H,WANG P.Introduction to GPU parallel programming technology [J].E-Scinence Technology & Application,2013,4(1):81-87

[35] VIRIEUX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

(編輯：朱文杰)

Simplified hybrid domain FWI method and GPU acceleration

GUI Sheng1,2,LIU Hong1,2,,LI Fei1,2

(1.InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China;2.KeyLaboratoryofPetroleumResourcesResearch,ChineseAcademyofSciences,Beijing100029,China)

Full waveform inversion comprehensively utilizes the dynamics and kinematics information of the prestack seismic wave field to rebuild accurately the parameter field of underground model.It is the important research orientation in geophysical exploration field of domestic and aboard,but the huge amount of computation has been an important factor which restricts its development.Computing capacity of high performance computing cluster currently consisting of the development of GPU provides an important basis of hardware for the improvement of the full waveform inversion (FWI) to make a more rapid 3D FWI,using the simplified hybrid domain FWI algorithm.We briefly introduce some optimization techniques of GPU acceleration applied to the simplified hybrid domain FWI,including the thread scheduling,GPU and GPU data transfer and the use of shared memory.The test shows that hybrid domain FWI can get a good inversion results.GPU technology can effectively improve the computation efficiency of the FWI and its speedup ratio is ten times more than CPU.

GPU technology,high performance computing (HPC),GPU technology optimization,full waveform inversion (FWI)

2016-09-30;改回日期：2016-12-08。

桂生(1988—),男,博士,主要從事速度建模和計(jì)算地球物理研究。

劉洪(1959—),男,研究員,主要從事地震波成像和油儲(chǔ)地球物理研究。

國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)(2012AA061202)資助。

P631

A

1000-1441(2017)01-0099-08

10.3969/j.issn.1000-1441.2017.01.012

This research is financially supported by the National High-tech R & D Program (863 Program) (Grant No.2012AA061202).