桂 生,劉 洪,李 飛

(1.中國科學院地質與地球物理研究所,北京100029;2.中國科學院油氣資源研究重點實驗室,北京100029)

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

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

(1.中國科學院地質與地球物理研究所,北京100029;2.中國科學院油氣資源研究重點實驗室,北京100029)

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

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

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

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

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

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

1.1 GPU架構

GPU體系架構如圖1所示。

圖1 GPU運算單元結構[21]

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

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

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

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

圖2 CUDA編程模式[21]

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

1.2 線程調度

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

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

圖3 CPU的微觀架構(a)和GPU的微觀架構(b)對比

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

圖4 不同線程塊計算時間對比

1.3 多GPU優(yōu)化

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

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

圖5 多GPU計算性能對比

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

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

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

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

圖6 P2P技術原理

圖7 共享內存示意

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

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

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

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

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

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

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

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

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

讀取各個節(jié)點的硬件信息;

if(小于單個GPU內存)

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

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

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

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

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

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

}

else(大于單個GPU內存)

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

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

并將每個節(jié)點上單炮數(shù)據(jù)進行區(qū)域分解,載入對應的GPU

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

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

}}

3 模型測試

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

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

圖11 Overthrust真實模型

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

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

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

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

4 結論

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

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

[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] 趙改善.地球物理高性能計算的新選擇:GPU 計算技術[J].勘探地球物理進展,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] 李博,劉國峰,劉洪.地震疊前時間偏移的一種圖形處理器提速實現(xiàn)方法[J].地球物理學報,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異構平臺的全波形反演及其實用化分析[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].科研信息化技術與應用,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—),男,博士,主要從事速度建模和計算地球物理研究。

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

國家高技術研究發(fā)展計劃(863計劃)(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).