基于GPU/CPU和震源隨機編碼技術的混合域全波形反演

馮海新,劉 洪,孫 軍,胡 婷,劉志偉

(1.中國科學院地質(zhì)與地球物理研究所,北京100029;2.中國科學院油氣資源研究重點實驗室,北京100029;3.中國科學院大學,北京100049;4.唐山學院,河北唐山063020;5.中國地質(zhì)科學院,北京100037)

基于GPU/CPU和震源隨機編碼技術的混合域全波形反演

馮海新1,2,3,劉 洪1,2,3,孫 軍4,胡 婷1,2,3,劉志偉5

(1.中國科學院地質(zhì)與地球物理研究所,北京100029;2.中國科學院油氣資源研究重點實驗室,北京100029;3.中國科學院大學,北京100049;4.唐山學院,河北唐山063020;5.中國地質(zhì)科學院,北京100037)

傳統(tǒng)的全波形反演利用普通炮集進行反演,反演計算量過大;且利用傳統(tǒng)的相位編碼技術進行全波形反演,會產(chǎn)生炮間串擾問題,因此,提出了基于GPU/CPU和震源隨機編碼技術的混合域全波形反演。該方法將參與反演的多個炮集隨機組合并分成炮集數(shù)相同的組,各組炮集疊加形成多個組合炮集,然后將組合炮集代替普通炮集進行反演。與傳統(tǒng)的相位編碼反演方法相比,震源隨機編碼技術在反演效率和收斂速度方面均有優(yōu)勢,且減少了炮間串擾噪聲;并且在GPU的加速下,計算效率會再次提升。Marmousi模型數(shù)據(jù)測試結(jié)果表明:組合炮集方法得到了與普通炮集方法相同的反演效果,但計算效率卻比普通炮集方法明顯提高,且相較于傳統(tǒng)的相位編碼技術,組合炮集方法有效抑制了串擾噪聲。

混合域;全波形反演;GPU/CPU;組合炮集;震源隨機編碼

全波形反演(FWI)綜合利用了地震記錄中振幅、走時和相位等信息,通過擬合實際波場和預測波場來定量提取地下介質(zhì)的彈性參數(shù),進而為勘探地震成像及速度建模、深部大尺度構(gòu)造演化分析[1-2]等提供可靠依據(jù)。

TARANTOLA[3]在1984年提出了基于廣義最小二乘反演理論的時間域全波形反演,PRATT等[4]和SHIN等[5]在其基礎上進一步發(fā)展了頻率域和Laplace-Fourier域波形反演方法。頻率域和Laplace-Fourier域波形反演本質(zhì)上都是多尺度的反演方法,兩者都遵循低頻到高頻的反演策略。利用SIRGUE等[6]的頻率域優(yōu)選策略,基于直接法的二維頻率域反演效果相對于時間域反演效果更加有效。對于三維頻率域波形反演問題,矩陣分解法有著極大的內(nèi)存需求和計算量,使得大尺度三維頻率域全波形反演問題變得困難。迭代法可以解決這個問題,但會失去頻率域算法多炮高效正演的優(yōu)勢,并且隨著頻率的升高和模型復雜程度的增加,其收斂性也會變差。Laplace-Fourier域方法是一種使用復數(shù)頻率的頻率域方法(頻率虛部代表地震記錄在時間方向上的衰減系數(shù)),該方法可以恢復大尺度構(gòu)造,并且能降低反演算法對初始模型的依賴;但是,作為一種頻率域方法的擴展,它繼承了頻率域方法的缺點,即基于直接法的三維模型需要很大的內(nèi)存和計算量。

考慮到這些問題,SIRGUE等[7]通過離散傅里葉變換,將時間域FWI和頻率域FWI結(jié)合起來,提出了混合域FWI方法。該方法將波傳播在時間域進行,即在時間域正演,利用離散傅里葉正變換(DFT)提取響應的頻率成分,進而在頻率域構(gòu)造梯度場。這種方法可用最小的耗時來獲得任意頻率成分,尤其當GPU參與到程序并行計算時,其計算效率會更高。此外,相對頻率域窗函數(shù)在混合域更易添加,從而得到某些波至的特征,如早至波、反射波等,這些波經(jīng)常被用來降低反演的非線性。KIM等[8]將這種混合域的思想應用到了Laplace-Fourier域,并將其應用到了三維寬方位角的實際資料FWI中。由于時間域算法相對頻率域和Laplace-Fourier域有著更少的內(nèi)存占用和計算量,所以混合域的方法使得三維多尺度波形反演具備實現(xiàn)的可能性。JUN等[9]進一步將該思想應用到了二維時間-Laplace-Fourier混合域彈性波全波形反演。

無論是時間域、頻率域、Laplace域,還是混合域反演,FWI都面臨著海量計算的問題。FWI技術已逐漸由理論模型研究發(fā)展到實際數(shù)據(jù)應用研究,然而計算量超大制約了該技術在實際生產(chǎn)中的廣泛應用。針對FWI的大計算量問題,DIAZ等[10]提出了一種改進的策略,利用隨機抽取的炮集在時間域進行全波形反演,相對于未改進的波形反演,在每次迭代時,并不是讓全部炮集參與反演,而是隨機抽取炮集參與反演。在此基礎上,HA等[11]采取對普通炮集進行循環(huán)采樣,在每次迭代時僅僅對循環(huán)采樣選取的炮集進行反演計算。這兩種方法都能在一定程度上提高反演的計算效率,但是難以從根本上解決全波形反演的大計算量問題。KREBS等[12]基于多炮同時反演的思想提出了一種基于相位編碼技術的多炮同時反演方法,利用隨機產(chǎn)生的相位編碼函數(shù)分別與震源實際地震記錄褶積,形成一個超級震源和超級炮記錄,這樣對于所有炮的反演計算即可轉(zhuǎn)換成對一個超級炮集記錄進行反演,大幅度提高了全波形反演的計算效率(計算效率約可提高炮數(shù)比的倍數(shù))。由于該方法是對一個超級炮記錄進行反演,即多炮同時反演,因此會產(chǎn)生大量的炮間串擾噪聲(cross-talk),從而影響整個全波形反演的計算精度和收斂效率。

GPU俗稱“顯卡”,是一款用于計算機顯示的硬件設備,是一種高并行度、多線程、計算能力巨大和存儲器帶寬極高的多核處理器。CUDA(Compute Unified Device Architecture)語言是一種基于指令集架構(gòu)和新的并行編程模型的通用計算架構(gòu),能夠利用GPU的并行計算引擎,比CPU更高效地解決許多細粒度的計算任務。在石油工業(yè)中,尤其是地震數(shù)據(jù)處理領域,GPU受到越來越多的重視。LI等[13]率先將GPU并行計算應用于非對稱走時Kirchhoff疊前時間偏移,大幅度提升了成像計算效率。LIU等[14]在GPU上實現(xiàn)了三維逆時偏移,并給出了程序優(yōu)化策略。MICIKEVICIUS[15]用GPU實現(xiàn)了三維有限差分程序,并提出了單通法和雙通法兩種優(yōu)化策略,也是目前常用的程序優(yōu)化方法。GPU在石油工業(yè)普及的趨勢,已經(jīng)逐漸將GPU推舉為下一代地震資料處理的主要計算平臺。

針對傳統(tǒng)混合域全波形反演策略,本文提出了一種基于GPU/CPU和震源隨機編碼技術的混合域全波形反演方法,在頻率域迭代時,利用震源隨機編碼技術,將實際炮記錄隨機組合,分成若干組,每組炮集數(shù)相同,然后將每組炮記錄疊加,形成一個組合炮集或超級炮集(Super-Shot)。最后,利用Marmousi模型數(shù)據(jù)驗證了本文方法的有效性。

1 方法原理

1.1 混合域反演原理

全波形反演按照波場求取方式可分為:時間域全波形反演、頻率域全波形反演和Laplace-Fourier域全波形反演?；旌嫌蛉ㄐ畏囱菔窃跁r間域進行正演、在頻率域進行反演的一種反演方法。頻率域聲波方程可以寫成:

(1)

式中:f(ω)為震源項;S(v,ω)為系數(shù)矩陣,與模型參數(shù)、差分離散格式、頻率和邊界條件有關;v是速度模型參量,反演即是對v不斷更新;u(ω)為頻率域的聲波波場。

梯度場為:

(2)

式中:F表示所有虛震源的集合;δd*表示殘差波場的復共軛;r表示殘差反傳波場。

(3)

如果S(v,ω)矩陣對稱,則有:

(4)

如果S(v,ω)矩陣不對稱,則需要用到LU分解的轉(zhuǎn)置。

第i個模型參量對應的梯度場可以表示為:

(5)

由公式(1)得出,頻率域聲波方程系數(shù)矩陣的微擾可以表示為:

將公式(7)代入公式(5)得到頻率域梯度求取公式:

(8)

1.2 震源隨機編碼技術

普通震源編碼技術,就是對單炮集進行編碼,選取固定數(shù)目的單炮,形成一個組合炮,然后以組合炮為單位來計算公式(8)的梯度值。該方法的缺點是,當多炮疊加到一起時,炮與炮之間會相互干擾。利用隨機編碼技術可以消除這種不良影響。隨機編碼可以在兩個方向上隨機:①在橫向方向上使炮點位置隨機;②沿著縱向方向給各個炮點加上時延,并使時延長度隨機。本文以炮點位置隨機為例說明隨機編碼的實現(xiàn)方式和有效性。

組合炮集對應的目標函數(shù)可表示為:

式中:Ucal,Uobs分別為模擬數(shù)據(jù)組合炮集、實際數(shù)據(jù)組合炮集;shot為單炮記錄;Ms是組合炮集的個數(shù);xi是組合炮集的位置,由M個普通炮集疊加而成。i為Rand函數(shù)產(chǎn)生1到Ns(單炮個數(shù))的一組隨機數(shù),每一個隨機數(shù)對應一個炮點位置(即震源位置),即對震源進行編碼;利用該隨機數(shù)可以有效壓制炮間串擾現(xiàn)象。模擬組合炮集是由組合震源激發(fā)產(chǎn)生,Ucal(xj,ω;xi)滿足公式:

(11)

其中,DFT表示傅里葉正變換。

(12)

組合炮集的個數(shù)和普通炮集個數(shù)滿足關系式:

(13)

式中:floor代表向下取整運算。

理論上,采用組合炮集方法可使計算量降低為普通炮集方法的1/M倍。但是炮點的組合會引入串擾噪聲,即波場之間的相互重疊造成所在位置對應的模型參量更新有誤。我們對普通炮集進行隨機編碼,形成Ms個組合炮和Ms個組合震源,將組合炮代替普通炮集參與反演計算,使串擾噪聲產(chǎn)生的可能性降低。且每次迭代時隨機編碼改變,這樣每次對應的地下重疊位置會隨機變化,避免相同的串擾噪聲作用于每次迭代,隨著疊加和迭代的進行,串擾噪聲便會被壓制。因此，在反演精度和收斂速度方面,基于震源隨機編碼技術的全波形反演要優(yōu)于傳統(tǒng)的相位編碼技術。當M=1時,便是普通炮集反演技術。當組合炮集個數(shù)較少時,其反演效果要比普通炮集反演效果差,這時,我們可以通過增加組合炮集個數(shù)來改善反演效果。

1.3 GPU/CPU協(xié)同并行計算

GPU將問題分解成線程塊的網(wǎng)格,每塊包含多個線程,如圖1所示,其中黃色部分表示線程塊的網(wǎng)格,綠色部分表示每一個網(wǎng)格中所包含的線程塊,藍色部分表示每一個線程塊所包含的線程束。塊可以按照任意順序執(zhí)行。CUDA內(nèi)部有自建的變量ThreadIdx,BlockIdx來標明線程號和塊號。一個線程可以通過線程號、塊號以及線程維度信息來獲取其在所有線程中的一個標號。在一維情況下標號等于ThreadIdx.x,在Block大小為(Bx,By)的二維時,標號為ThreadIdx.x+ThreadIdx.y×Bx。而在(Bx,By,Bz)的三維Block情況下,線程標號為ThreadIdx.x+ThreadIdx.y×Bx+ThreadIdz.z×By×Bz。在Grid下的不同Block中線程之間的執(zhí)行完全是并行的,沒有同步機制,也不能夠串行執(zhí)行。

圖1 GPU線程組成

基于CUDA進行編程時,主機作為Host,主要有兩個功能:①利用主機函數(shù)(Host Function)負責程序的串行計算;②負責CPU和GPU之間的數(shù)據(jù)傳輸,并行計算時,CPU讀取輸入數(shù)據(jù)并復制到GPU卡中,計算結(jié)束將結(jié)果復制到CPU中,為其串行計算提供數(shù)據(jù)。GPU作為設備(Device),其內(nèi)運行的函數(shù)稱之為核函數(shù)(Kernel Function),主要負責程序的并行部分。

我們所采用的硬件設備是單GPU卡節(jié)點。該節(jié)點由一個CPU和一個NVIDIA GeForce GTX 580顯卡組成。每炮數(shù)據(jù)的梯度求取涉及多個GPU模塊,包括:執(zhí)行空間高階有限差分GPU模塊,PML吸收邊界GPU模塊,地表接收波場提取的GPU模塊,波場正傳時的離散傅里葉變換GPU模塊,時間域殘差求取的GPU模塊,目標函數(shù)值求取的GPU模塊,殘差波場加載的GPU模塊,殘差波場反傳時的離散傅里葉變換GPU模塊,梯度場求取的GPU模塊,求取梯度場最大值的GPU模塊,優(yōu)化步長求取的GPU模塊,更新速度模型的GPU模塊。這些模塊都是在GPU上完成?；贕PU/CPU協(xié)同計算的偽代碼如圖2所示。

圖2 基于GPU/CPU和震源隨機編碼技術的全波形反演偽代碼

圖2中的主要函數(shù)和參數(shù)說明見表1,其中后綴“cuda”表示在GPU上運行的程序。

1.4 實現(xiàn)流程

本文根據(jù)GPU/CPU的加速方法和CUDA代碼的實現(xiàn)方式,結(jié)合震源隨機編碼技術,實現(xiàn)了Host函數(shù)和Kernel函數(shù)。本文方法實現(xiàn)流程如圖3所示,其中綠色部分表示在GPU上實現(xiàn)的模塊,紅色方塊部分表示本文方法相對于傳統(tǒng)混合域方法[16]的改進之處。

表1 圖2中主要函數(shù)及參數(shù)說明

圖3 本文方法實現(xiàn)流程

2 模型應用及效果分析

采用Marmousi模型,對本文方法進行測試。該模型網(wǎng)格大小為400×120,網(wǎng)格間距均為5m?；谠撃Ｐ徒⒌拦潭?、炮移動的觀測系統(tǒng)。起始炮位于網(wǎng)格點(3,2)處,炮點個數(shù)為65,炮間距為30m,每一炮均采用雙邊接收,檢波器個數(shù)固定為400個,檢波器間距為5m,采樣時間長度為2s,時間采樣間隔為0.001s。我們采用空間優(yōu)化12階,時間2階差分格式進行正演,采用PML邊界條件,其中,PML邊界大小為20層。震源采用主頻為18Hz的Richer子波。初始速度模型采用平滑速度模型,反演算法采用擬Hessian矩陣法,步長求取采用步長衰減法,雖然該步長求取方法精度較低,但計算量小,更多的計算資源主要用于梯度場求取。圖4和圖5分別為真實速度模型和初始速度模型。我們分別采用組合炮集和普通炮集進行反演。

本文采用等差頻率進行反演,從3.5Hz開始,每次增加1.5Hz,直到39.5Hz,總共25個頻率,分成24個頻率組,每個頻率的最大迭代次數(shù)為18。我們選取第6組、第12組、第18組及第24組的反演結(jié)果進行對比。圖6,圖7,圖8和圖9分別為采用普通炮集方法和本文組合炮集方法在第6組、第12組、第18組和第24組的反演結(jié)果。分別抽取真實速度模型、初始速度模型、普通炮集方法和組合炮集方法的反演結(jié)果在水平網(wǎng)格點200處(水平方向1000m處)的速度曲線進行對比,結(jié)果如圖10所示。對比圖6,圖7,圖8和圖9,并分析圖10,可見,組合炮集反演結(jié)果與普通炮集反演結(jié)果相差無幾,兩種方法的反演速度曲線擬合得較好。

我們采用的迭代終止條件為:|fmod-fmod_pre|/fmod<0.001,利用(‖V實際-V反演‖2/‖V實際‖2)×100%來評價反演精度。兩種方法得到的迭代終止條件值和反演精度如表2所示。

分析表2可知,組合炮集方法的第12組和普通炮集方法的第18組、第24組的迭代終止條件值均小于0.001,表示反演沒有達到最大迭代次數(shù),便會終止迭代,進入下一組頻率反演計算。并且,第6組、第12組、第18組以及第24組,每組組合炮集方法的反演精度均低于普通炮集反演方法的反演精度,但隨著反演頻率組數(shù)的增加,兩種方法的反演精度之差在逐漸縮小,到最后一組反演頻率,反演精度之差已經(jīng)降低到0.26%,充分說明本文方法的反演效果與普通炮集方法的反演結(jié)果相當。

圖4 真實速度模型

圖5 初始速度模型

圖6 第6組速度反演結(jié)果a 普通炮集方法; b 組合炮集方法

圖7 第12組速度反演結(jié)果a 普通炮集方法; b 組合炮集方法

圖8 第18組速度反演結(jié)果a 普通炮集方法; b 組合炮集方法

圖9 第24組速度反演結(jié)果a 普通炮集方法;b 組合炮集方法

圖10 第6組(a)、第12組(b)、第18組(c)和第24組(d)反演速度曲線

由于實驗模型較小,本文在單GPU卡上進行計算,組合炮集混合域波形反演所用時間為105s左右,普通炮集混合域波形反演所用時間為530s左右,組合炮集反演所用時間約為普通炮集反演所用時間的1/5。

與GPU代碼對應,我們編寫了相應的CPU代碼,從而測試了兩種代碼的執(zhí)行效率。此次測試我們所采用的硬件設備是一個聯(lián)想ThinkCentre小型工作站,CPU為主頻3.0Hz的8核i7處理器,內(nèi)存8GB,并且配備了Kelper架構(gòu)的GTX680顯卡。該顯卡共有1536個處理器,2GB的共享內(nèi)存,時鐘頻率為1.06Hz。CUDA版本為4.2。針對該模型,GPU代碼需要105s完成一次迭代,而CPU代碼需要2971s,加速比為28.3。

表2 迭代終止條件值及反演精度

3 結(jié)論

本文提出并實現(xiàn)了基于GPU/CPU和震源隨機編碼技術的混合域全波形反演,利用GPU/CPU協(xié)同計算,提高了計算效率,并在此基礎上對混合域全波形反演由利用普通炮集反演變?yōu)槔秒S機組合炮集進行反演,實現(xiàn)了震源隨機編碼技術在混合域全波形反演中的應用。通過理論分析和模型實驗得出如下結(jié)論:

1) 震源隨機編碼技術通過對隨機炮的選取,壓制炮間串擾噪聲,從而獲得比傳統(tǒng)相位編碼技術更高的反演精度和收斂效率。理論上講,在相同的迭代次數(shù)和反演頻率的情況下,可以通過增加隨機組合炮集數(shù)來達到與普通炮集反演相同的精度。

2) 利用GPU卡進行全波形反演大幅度提高了混合域全波形反演計算效率。

[1] BLEIBINHAUS F,HOLE J A,RYBERG T,et al.Structure of the California Coast Ranges and San Andreas Fault at SAFOD from seismic waveform inversion and reflection imaging[J].Journal of Geophysical Research:Solid Earth,2007,112(B6):3015-3018

[2] OPERTO S,VIRIEUX J,DESSA J X,et al.Crustal seismic imaging from multifold ocean bottom seismometer data by frequency domain full waveform tomography:application to the eastern Nankai trough[J].Journal of Geophysical Research:Solid Earth,2006,111:B09306

[3] TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

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

[5] SHIN C,CHA Y H.Waveform inversion in the Laplace—Fourier domain[J].Geophysical Journal International,2009,177(3):1067-1079

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

[7] SIRGUE L,ETGEN J,ALBERTIN U.3D frequency domain waveform inversion using time domain finite difference methods[J].Expanded Abstracts of 70thEAGE Annual Conference,2008:F022

[8] 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

[9] JUN H,KIM Y,SHIN C,et al.2D elastic time-Laplace-Fourier-domain hybrid full waveform inversion[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:1008-1013

[11] HA W,SHIN C.Efficient full waveform inversion using a cyclic shot sampling method[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-5

[12] KREBS J R,ANDERSON J E,HINKLEY D,et al.Fast full-wavefield seismic inversion using encoded sources[J].Geophysics,2009,74(6):WCC105

[13] 李博,劉國峰,劉洪.地震疊前時間偏移的一種圖形處理器提速實現(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

[14] LIU G F,LIU Y N,REN L,et al.3D seismic reverse time migration on GPU[J].Computers & Geosciences,2013,59(3):17-23

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

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

(編輯：顧石慶)

Hybrid domain full waveform inversion based on GPU/CPU and source random coding technique

FENG Haixin1,2,3,LIU Hong1,2,3,SUN Jun4,HU Ting1,2,3,LIU Zhiwei5

(1.InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China;2.KeyLaboratoryofPetroleumResourcesResearch,Beijing100029,China;3.UniversityofChineseAcademyofSciences,Beijing100049,China;4.TangshanCollege,Tangshan063020,China;5.ChineseAcademyofGeologicalSciences,Beijing100037,China)

For the conventional waveform inversion we use the common shot gathers,which consumes too much computation capacity.The conventional phase coding technique to carry out the full waveform inversion will induce the problem of cross talk between the shots.In order to solve these problems,we propose a hybrid domain full waveform inversion based on GPU/CPU and source random coding technique.This method randomly combines the shot gathers and divides the shot gathers into groups with the same number of shot gather.Then,the common shot gathers are replaced by the combined shot gathers to carry out inversion.Compared with the traditional phase coding technique,source random coding technique can effectively increase the efficiency and the convergence rate and reduce the cross talk noise.With the acceleration of the GPU computation,the full waveform inversion efficiency of the hybrid domain is effectively improved further.The testing result with Marmousi model shows that the common shot gathers and the combined shot gathers have the same inversion result,but the computational efficiency of the combined shot gathers is significantly higher than that of the common shot gathers.Compared with the traditional phase encoding technique,crosstalk is effectively suppressed by the combined shot gathers.

hybrid domain,full waveform inversion,GPU/CPU,combined shot gathers,source random coding

2016-07-20;改回日期：2016-11-25。

馮海新(1987-),男,博士在讀,主要從事油儲地球物理方面的研究。

劉洪(1959-),男,研究員,主要從事油儲地球物理方面的研究。

國家重點研發(fā)計劃(2016YFC0601101)和中國石油集團“彈性波地震成像技術合作研發(fā)課題”(2015A-3613)聯(lián)合資助。

P631

A

1000-1441(2017)01-0107-09

10.3969/j.issn.1000-1441.2017.01.013

This research is financially supported by the National Key Basic Research Program of China (Grant No.2016YFC0601101) and the Collaboration Project of Seismic Imaging in Elastic Medium of China National Petroleum Corporation (Grant No.2015A-3613).