宋　鵬,王修田

(1.中國海洋大學(xué)海洋地球科學(xué)學(xué)院,山東青島266100;2.青島國家海洋科學(xué)與技術(shù)實(shí)驗(yàn)室,山東青島266100;3.海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東青島266100)

?

基于多卡GPU的隨機(jī)炮分配相位編碼全波形反演

宋鵬1,2,3,王修田1,2,3

(1.中國海洋大學(xué)海洋地球科學(xué)學(xué)院,山東青島266100;2.青島國家海洋科學(xué)與技術(shù)實(shí)驗(yàn)室,山東青島266100;3.海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東青島266100)

摘要:針對傳統(tǒng)時間域相位編碼全波形反演的炮間串?dāng)_問題,提出了一種基于隨機(jī)炮分配的相位編碼時間域全波形反演策略,并實(shí)現(xiàn)了基于多卡GPU的加速?；谒目℅PU加速的模型反演實(shí)驗(yàn)表明：在GPU和相位編碼技術(shù)的共同加速下,時間域全波形反演的計算效率可得到顯著提升,同時由于各卡間的隨機(jī)炮分配,基于隨機(jī)炮分配的相位編碼時間域全波形反演策略能更有效地壓制串?dāng)_噪聲,因此其比傳統(tǒng)的相位編碼全波形反演具有更高的反演精度和收斂效率。

關(guān)鍵詞：全波形反演;相位編碼;串?dāng)_;多卡GPU

時間域全波形反演方法由TARANTOLA首先提出[1-2],其基于波動方程理論,充分利用實(shí)測地震記錄的整體信息(走時、振幅、相位等)來重建地下介質(zhì)速度結(jié)構(gòu)。對于理想觀測系統(tǒng),理論上全波形反演的分辨率可達(dá)到波長數(shù)量級[3-5]。

近20年來,時間域全波形反演方法獲得了長足的發(fā)展。BUNKS等[6]提出了時間域的多尺度全波形反演策略,有效避免了局部極值,提高了全波形反演的精度和穩(wěn)定性；MAO等[7]和ZHANG等[8]采用邊界存儲策略顯著降低了全波形反演的內(nèi)存消耗；BROSSIER等[9]、張生強(qiáng)等[10]以及苗永康[11]實(shí)現(xiàn)了基于L-BFGS算法的全波形反演；WANG等[12]提出并實(shí)現(xiàn)了基于混合迭代優(yōu)化方法的全波形反演,大大提高了反演的收斂效率和深部地層的反演精度；ZHANG等[13-14]提出了一種基于地震波能量的梯度預(yù)處理方法,其應(yīng)用聲波方程模擬的波場信息對梯度進(jìn)行預(yù)處理,可在有效提高深層反演精度的同時避開Hessian矩陣或其逆矩陣的計算和存儲。針對由于低頻缺失而導(dǎo)致的反演初始模型難以獲得的問題,WU等[15]提出了一種基于地震包絡(luò)的反演方法,其應(yīng)用地震數(shù)據(jù)的包絡(luò)替代低頻信息來反演速度模型的長波長結(jié)構(gòu),為全波形反演提供了一個較為精確的初始模型。

時間域全波形反演技術(shù)已逐漸由理論模型研究向?qū)嶋H數(shù)據(jù)應(yīng)用方向發(fā)展,然而超大計算量制約了該技術(shù)在實(shí)際生產(chǎn)中的廣泛應(yīng)用。針對全波形反演的大計算量問題,DIAZ等[16]提出了一種基于隨機(jī)炮抽取的時間域全波形反演策略,在每次迭代時僅隨機(jī)抽取部分炮參與反演,而不是讓全部炮集參與計算;HA等[17]提出了一種基于循環(huán)炮采樣的時間域全波形反演,同樣在每次迭代時不對所有炮進(jìn)行反演計算,而是通過循環(huán)采樣模式選取部分炮參與反演。以上兩種方法均減少了每次迭代時參與計算的炮數(shù),能夠在一定程度上提高反演的計算效率,但其顯然僅能適應(yīng)于覆蓋次數(shù)較高的觀測系統(tǒng),且難以從根本上解決全波形反演的大計算量問題。2009年,KREBS等[18]基于多炮同時反演的思想,提出了一種基于相位編碼技術(shù)的多炮同時反演方法(以下簡稱為相位編碼全波形反演),將所有炮的震源褶積上一個隨機(jī)產(chǎn)生的相位編碼函數(shù),生成一個超級震源,并應(yīng)用相同的編碼函數(shù)對實(shí)際地震記錄進(jìn)行褶積,生成超級炮集記錄,這樣對于所有炮的反演計算即可轉(zhuǎn)換成對一個超級炮集記錄進(jìn)行反演,大幅度提高了全波形反演的計算效率(計算效率約可提高炮數(shù)倍)。但相位編碼全波形反演在提高計算效率的同時,也由于多炮的同時反演(即僅對一個超級炮集記錄進(jìn)行反演)而產(chǎn)生大量的炮間串?dāng)_噪聲,從而影響整個全波形反演的計算精度和收斂效率。

GPU憑借其優(yōu)異的浮點(diǎn)運(yùn)算能力和更快的訪存速度已成為新一代的高性能計算工具。為方便開發(fā)者編程以實(shí)現(xiàn)算法的GPU加速,NVIDIA公司于2007年推出了統(tǒng)一計算設(shè)備架構(gòu)(CUDA)平臺,大大簡化了基于GPU的算法研發(fā)過程。自CUDA推出后,GPU的高性能計算在地球物理界得到廣泛應(yīng)用[19-21],并且目前已有學(xué)者對于時間域全波形反演的GPU加速進(jìn)行了研究和探討[7,22-23],其研究成果表明,即使是采用單卡GPU,時間域全波形反演也能獲得幾十倍的加速比;而當(dāng)前多卡GPU集群系統(tǒng)已逐漸普及,因此發(fā)展基于多卡GPU的全波形反演算法,以借助于多卡GPU集群的高性能并行計算能力來提高全波形反演的計算效率將是未來的發(fā)展趨勢。

本文提出了一種基于多卡GPU的隨機(jī)炮分配相位編碼時間域全波形反演策略,在每次迭代時，首先根據(jù)GPU的卡數(shù)將所有炮隨機(jī)分成若干組,再將每組炮進(jìn)行相位編碼形成獨(dú)立的超級炮集參與反演計算,最后將各個獨(dú)立超級炮集的計算結(jié)果合成總梯度和總目標(biāo)函數(shù)。實(shí)驗(yàn)表明,這種對隨機(jī)生成的多個獨(dú)立超級炮集的反演技術(shù)可有效壓制炮間串?dāng)_的影響,提高迭代精度和收斂效率;同時在GPU和相位編碼技術(shù)的共同加速下,時間域全波形反演的計算效率也得到大幅度提升。

1基于隨機(jī)炮分配的相位編碼全波形反演

1.1隨機(jī)炮分配相位編碼全波形反演目標(biāo)函數(shù)

設(shè)在一個計算平臺中能夠應(yīng)用的GPU卡數(shù)為M,隨機(jī)分配到第m個卡(組)的炮集數(shù)為Nm,則最小平方意義下的隨機(jī)炮分配相位編碼全波形反演的目標(biāo)函數(shù)可寫為：

(1)

式中：c表示速度模型;sn為第m組中第n炮的震源子波;p為相應(yīng)的合成地震記錄;dn表示第m組中第n炮的觀測地震記錄;en表示第m組中第n炮的編碼函數(shù)序列(其與常規(guī)的相位編碼序列[18]類似,為保證各組中的炮間串?dāng)_得到有效壓制,en需隨機(jī)給出,一般情況下,當(dāng)j≠k時,ej≠ek);“?”表示褶積。

當(dāng)p(c,sn)是震源的線性函數(shù)時,(1)式可寫為：

(2)

由(2)式可知,基于隨機(jī)炮分配的相位編碼全波形反演將所有炮的反演計算轉(zhuǎn)換成由隨機(jī)分配形成的M個超級震源與M個超級記錄的反演,通過隨機(jī)炮的分配減少了參與編碼的實(shí)際炮數(shù),降低了串?dāng)_產(chǎn)生的幾率,如果在每次迭代時均重新隨機(jī)分配各組中的炮號(數(shù))則還能避免相同的串?dāng)_作用于每次迭代中,因此基于隨機(jī)炮分配的相位編碼全波形反演比傳統(tǒng)的相位編碼全波形反演有更高的反演精度和收斂效率。顯然,當(dāng)M=1時(即所有炮僅為1組),基于隨機(jī)炮分配的相位編碼全波形反演即轉(zhuǎn)化為常規(guī)的相位編碼全波形反演。

1.2隨機(jī)炮分配相位編碼全波形反演實(shí)現(xiàn)步驟

基于多卡GPU的隨機(jī)炮分配相位編碼全波形反演的具體實(shí)現(xiàn)步驟如下。

1) 根據(jù)GPU卡數(shù)將所有炮隨機(jī)分為若干組(每次迭代時均重新隨機(jī)分配各卡對應(yīng)的炮號(數(shù)))。

2) 將每組炮分別與一個隨機(jī)的編碼函數(shù)序列褶積,生成若干個超級炮集,并將其對應(yīng)的震源子波分別與相同的編碼函數(shù)序列褶積,生成若干個超級震源。

3) 在每一個GPU卡上,分別以一個超級震源作為正時擾動進(jìn)行正演模擬得到正時波場u(x,t,Si)(這里x表示空間坐標(biāo),t表示時間,i表示卡號)和模擬記錄,并進(jìn)一步得到目標(biāo)函數(shù)和超級記錄殘差,然后以該超級記錄殘差作為逆時擾動計算得到逆時波場ψ(x,t,Si),并按(3)式計算得出各個超級記錄對應(yīng)的梯度值：

(3)

式中：di表示第i個GPU卡計算得到的梯度;c(x)表示速度模型;“·”表示時間的一階微分;T表示最大時間。

4) 將各個GPU卡計算得出的梯度值合成總梯度同時得到總目標(biāo)函數(shù),并進(jìn)一步計算得到迭代步長,然后修改模型,即完成本次迭代。

5) 重復(fù)步驟1)到步驟4),直到目標(biāo)函數(shù)滿足要求或達(dá)到最大迭代次數(shù)。

基于多卡GPU的隨機(jī)炮分配相位編碼全波形反演的具體反演流程見圖1(以四卡為例)。

圖1　基于四卡GPU的隨機(jī)炮分配相位編碼反演流程

2相位編碼全波形反演的GPU加速

本文首先基于消息傳遞接口(MPI),根據(jù)卡數(shù)啟動相應(yīng)數(shù)目的進(jìn)程,并將每個進(jìn)程與對應(yīng)的GPU卡綁定;而在每個GPU卡內(nèi)則基于CUDA,實(shí)現(xiàn)相位編碼全波形反演的GPU加速,這樣基于多卡GPU的隨機(jī)炮分配相位編碼全波形反演的計算效率與基于單卡加速的傳統(tǒng)相位編碼全波形反演相當(dāng)。

基于CUDA進(jìn)行編程時,CPU作為主機(jī),GPU被視作設(shè)備。主機(jī)上運(yùn)行的函數(shù)稱為主機(jī)函數(shù),設(shè)備上運(yùn)行的函數(shù)稱為內(nèi)核函數(shù)。主機(jī)函數(shù)負(fù)責(zé)串行部分,并負(fù)責(zé)將數(shù)據(jù)由內(nèi)存復(fù)制到顯存中以供GPU使用。內(nèi)核函數(shù)負(fù)責(zé)執(zhí)行高度線程化的并行處理任務(wù),并以線程網(wǎng)格(grid)的形式組織。每個線程網(wǎng)格又可分為若干個線程塊(block),每個線程塊包含若干個線程(thread)。通過這樣的組織形式,內(nèi)核函數(shù)可同時啟動成千上萬個線程。內(nèi)核函數(shù)執(zhí)行完畢后需將處理完的數(shù)據(jù)由顯存復(fù)制到內(nèi)存中。

本文根據(jù)相位編碼時間域全波形反演的計算模式和CUDA特點(diǎn),創(chuàng)建了相關(guān)的主機(jī)函數(shù)和核函數(shù)(具體名稱和功能見表1,其中核函數(shù)前綴為“__global__”),實(shí)現(xiàn)了相位編碼全波形反演的GPU加速。

表1　主機(jī)函數(shù)及核函數(shù)

應(yīng)用表1中的函數(shù),基于CUDA的相位編碼時間域全波形反演的一段偽碼如圖2所示。其中,“P”為內(nèi)存中的波場數(shù)組,“P_GPU”為顯存中的波場數(shù)組,“V”為內(nèi)存中的模型數(shù)組,“V_GPU”為顯存中的模型數(shù)組,“=>”表示由內(nèi)存復(fù)制到顯存,“<=”表示由顯存復(fù)制到內(nèi)存。

圖2　基于CUDA的相位編碼全波形反演偽碼

3模型反演效果分析

3.1層狀斷塊速度模型反演效果分析

層狀斷塊速度模型見圖3,該模型橫向距離為6000m,深度為1800m,共5層,其速度從上到下依次為2700,3200,3900,5000,5500m/s。基于該模型建立道固定、炮移動的觀測系統(tǒng),炮間距為25m,共241炮,道間距也為25m,共241道。

本文基于BUNKS等[6]提出的多尺度分頻反演策略,共分5個頻段(高截頻分別為5.0,9.0,16.3,29.6,50.0Hz)進(jìn)行反演,每個頻段各25次迭代,初始速度模型如圖4所示。

基于一個含4個卡的GPU工作站測試反演計算效率(GPU型號為NVIDIA Tesla K20)。測試結(jié)果顯示,基于單卡的傳統(tǒng)相位編碼反演(以下簡稱單卡反演)和基于四卡的隨機(jī)炮分配相位編碼反演(以下簡稱四卡反演)的單次迭代時間均僅為6.5s。

圖5為原始炮集記錄和超級炮記錄(包括單卡反演第1次迭代時合成的超級記錄、四卡反演時第1個卡第1次和第2次迭代時所合成的超級記錄);單卡反演和四卡反演每個頻帶各25次迭代的反演結(jié)果分別見圖6和圖7。

圖3　層狀斷塊速度模型

圖4　層狀斷塊初始速度模型

圖5　原始炮集記錄及超級炮記錄a 原始單炮記錄(炮號=120); b 單卡反演第1次迭代時的超級記錄; c 四卡反演中第1個卡第1次迭代時的超級記錄;d 四卡反演中第1個卡第2次迭代時的超級記錄

由圖5可知,四卡反演時每卡參與編碼的炮數(shù)要遠(yuǎn)遠(yuǎn)小于單卡反演,且每次迭代所參與編碼的炮都為隨機(jī)產(chǎn)生,這最大限度地降低了炮間串?dāng)_對于反演的影響。對比圖6和圖7可知,在同樣迭代次數(shù)下,四卡反演的精度明顯高于單卡反演,串?dāng)_噪聲大幅減少。

3.2逆掩斷層二維模型反演效果分析

圖6　層狀斷塊速度模型單卡25次迭代反演結(jié)果

圖7　層狀斷塊速度模型四卡25次迭代反演結(jié)果

圖8　逆掩斷層二維速度模型

本文采用的逆掩斷層二維模型(圖8)長度為10km,深為2.325km,橫、縱向網(wǎng)格步長均為12.5m。設(shè)計接收道固定、炮移動的觀測方式,每炮接收道相同,自左向右等間隔排放,共401道,道間距為25m,檢波點(diǎn)深度為25m;炮間距也為25m,炮點(diǎn)深度為25m,共401炮?；贐UNKS等[6]提出的多尺度分頻反演策略和圖9所示的初始模型分4個頻段(高截頻分別為5.0,11.0,23.7和47.5Hz)進(jìn)行反演,每個頻段各25次迭代。此外,還采用了ZHANG等[14]提出的基于地震波能量的梯度預(yù)處理方法以提高深層的反演精度。

反演同樣基于型號為Tesla K20的四卡GPU工作站,單卡與四卡反演的單次迭代時間均僅為13s。

圖9　逆掩斷層二維初始速度模型

圖10　逆掩斷層二維速度模型單卡25次迭代反演結(jié)果

圖11　逆掩斷層二維速度模型四卡25次迭代反演結(jié)果

圖12　逆掩斷層二維速度模型在x=4000m(a)與x=6525m(b)處25次迭代反演速度曲線(黑線代表真實(shí)速度,藍(lán)線代表初始模型速度,綠線代表單卡反演結(jié)果,紅線代表四卡反演結(jié)果)

圖10和圖11分別為單卡反演和四卡反演每頻帶各25次迭代的反演結(jié)果;圖12為水平方向4000和6525m處單卡和四卡反演的速度曲線。對比圖10和圖11以及圖12中單卡和四卡反演的速度曲線可知,25次迭代后四卡反演的結(jié)果明顯優(yōu)于單卡反演結(jié)果。

圖13　逆掩斷層二維速度模型單卡50次迭代反演結(jié)果

圖14　逆掩斷層二維速度模型四卡50次迭代反演結(jié)果

圖15　逆掩斷層二維速度模型在x=4000m(a)與x=6525m(b)處50次迭代反演速度曲線(黑線代表真實(shí)速度,藍(lán)線代表初始模型速度,綠線代表單卡反演結(jié)果,紅線代表四卡反演結(jié)果)

圖13和圖14分別為單卡反演和四卡反演每頻段各50次迭代后的反演結(jié)果;圖15為4000m和6525m處的反演速度曲線。對比圖13和圖14以及圖15中的單卡和四卡反演速度曲線可知,50次迭代后四卡反演的精度也明顯高于單卡反演結(jié)果。

本文應(yīng)用BEN-HADJ-ALI等[24]定義的一個反演精度評估因子來評估反演結(jié)果和真實(shí)模型的誤差,反演精度評估因子的表達(dá)式為：

(4)

式中：mresult表示最終反演結(jié)果;mreal表示真實(shí)速度。根據(jù)(4)式計算得到單卡反演25次迭代、四卡反演25次迭代、單卡反演50次迭代以及四卡反演50次迭代的Eres值見表2。

表2　逆掩斷層二維模型反演結(jié)果的Eres值

由表2可知,對于相同的迭代次數(shù),四卡反演的精度總是高于單卡反演,甚至四卡25次迭代的反演精度都要高于單卡50次迭代的反演結(jié)果,因此逆掩斷層二維模型基于四卡GPU的隨機(jī)炮分配相位編碼反演的收斂效率和反演精度都遠(yuǎn)遠(yuǎn)超過單卡相位編碼反演。

4結(jié)論與展望

本文提出并實(shí)現(xiàn)了基于多卡GPU的隨機(jī)炮分配相位編碼全波形反演,理論分析和模型實(shí)驗(yàn)得出如下結(jié)論：

1) 基于多卡GPU的隨機(jī)炮分配相位編碼全波形反演可通過卡間隨機(jī)炮的選取進(jìn)一步壓制炮間的串?dāng)_,從而獲得更高的反演精度和收斂效率。理論上講,所用的GPU卡數(shù)目越多,對炮間串?dāng)_的壓制效果越好,相應(yīng)的反演精度和收斂效率越高。

2) 基于多卡GPU的隨機(jī)炮分配相位編碼全波形反演在相位編碼技術(shù)的基礎(chǔ)上,借助于多卡GPU的加速能力,大幅度提高時間域全波形反演的計算效率。

顯然,當(dāng)GPU卡數(shù)達(dá)到實(shí)際炮數(shù)時,隨機(jī)炮分配相位編碼全波形反演即轉(zhuǎn)化為通常意義下的全波形反演,此時各卡僅處理其中一個炮集數(shù)據(jù),也就不存在炮間串?dāng)_問題。但在當(dāng)前(甚至今后相當(dāng)長的一段時間內(nèi))的計算環(huán)境下,海量的地震數(shù)據(jù)使得“GPU卡數(shù)等于實(shí)際炮數(shù)”的條件在實(shí)際生產(chǎn)應(yīng)用中很難滿足。因此,基于多卡GPU的隨機(jī)炮分配相位編碼全波形反演不失為一種較為理想的能適合于大規(guī)模數(shù)據(jù)的反演策略。

參考文獻(xiàn)

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

[2]TARANTOLA A.A strategy for nonlinear elastic inversion of seismic reflection data[J].Geophysics,1986,51(10):1893-1903

[3]WU R S,TOKSOZ M N.Diffraction tomography and multisource holography applied to seismic imaging[J].Geophysics,1987,52(1):11-25

[4]PRATT R G,SHIPP R M.Seismic waveform inversion in the frequency domain.part2:fault delineation in sediments using crosshole data[J].Geophysics,1999,64(3):902-914

[5]卞愛飛,於文輝,周華偉.頻率域全波形反演方法研究進(jìn)展[J].地球物理學(xué)進(jìn)展,2010,25(3):982-993

BIAN A F,YU W H,ZHOU H W.Progress in the frequency-domain full waveform inversion method[J].Progress in Geophysics,2010,25(3):982-993

[6]BUNKS C,SALECK F M,ZALESKI S,et al.Multiscale seismic waveform inversion[J].Geophysics,1995,60(5):1457-1473

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

[8]ZHANG Q C,ZHOU H,WANG J.Efficient boundary storage strategy for three-dimensional elastic full-waveform inversion in time domain[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:1142-1146

[9]BROSSIER R,OPERTO S,VIRIEUX J.Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[J].Geophysics,2009,74(6):WCC105-WCC118

[10]張生強(qiáng),劉春成,韓立國,等.基于L-BFGS算法和同時激發(fā)震源的頻率多尺度全波形反演[J].吉林大學(xué)學(xué)報(地球科學(xué)版),2013,43(3):1004-1012

ZHANG S Q,LIU C C,HAN L G,et al.Frequency multi-scale full waveform inversion based on L-BFGS algorithm simultaneous sources approach[J].Journal of Jilin University(Earth Science Edition),2013,43(3):1004-1012

[11]苗永康.基于L-BFGS算法的時間域全波形反演[J].石油地球物理勘探,2015,50(3):469-474

MIAO Y K.The time-domain full waveform inversion based on L-BFGS algorithm[J].Oil Geophysical Prospecting,2015,50(3):469-474

[12]WANG Y,DONG L G,LIU Y Z.Improved hybrid iterative optimization method for seismic full waveform inversion[J].Applied Geophysics,2013,10(3):265-277

[13]ZHANG Z G,LIN Y Z,HUANG L J.Full-waveform inversion in the time domain with an energy-weighted gradient[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2772-2776

[14]ZHANG Z G,HUANG L J,LIN Y Z.A wave-energy-based precondition approach to full-waveform inversion in the time domain[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-5

[15]WU R S,LUO J R,WU B Y.Seismic envelope inversion and modulation signal model[J].Geophysics,2014,79(3):WA13-WA24

[16]DIAZ E,GUITTON A.Fast full waveform inversion with random shot decimation[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2804-2808

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

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

[19]李博,劉國峰,劉洪.地震疊前時間偏移的一種圖形處理器提速實(shí)現(xiàn)方法[J].地球物理學(xué)報,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

[20]劉紅偉,李博,劉洪,等.地震疊前逆時偏移高階有限差分算法及GPU實(shí)現(xiàn)[J].地球物理學(xué)報,2010,53(7):1725-1733

LIU H W,LI B,LIU H,et al.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J].Chinese Journal of Geophysics,2010,53(7):1725-1733

[21]劉國峰,劉洪,李博,等.山地地震資料疊前時間偏移方法及其GPU實(shí)現(xiàn)[J].地球物理學(xué)報,2009,52(12):3101-3108

LIU G F,LIU H,LI B,et al.Method of prestack time migration of seismic data of mountainous regions and its GPU implementation[J].Chinese Journal of Geophysics,2009,52(12):3101-3108

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

[23]張猛,王華忠,任浩然,等.基于CPU/GPU異構(gòu)平臺的全波形反演及其實(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

[24]BEN-HADJ-ALI H,OPERTO S,VINEUX J.An efficient frequency-domain full waveform inversion method using simultaneous encoded sources[J].Geophysics,2011,76(4):R109-R124

(編輯：陳杰)

The phase-encoding full waveform inversion with random shots allocation based on the acceleration of multi-card GPU

SONG Peng1,2,3,WANG Xiutian1,2,3

(1.CollegeofMarineGeo-science,OceanUniversityofChina,Qingdao266100,China;2.QingdaoNationalLaboratoryforMarineScienceandTechnology,Qingdao266100,China;3.KeyLabofSubmarineGeosciencesandProspectingTechniquesMinistryofEducation,Qingdao266100,China)

Abstract:To deal with the crosstalk problem of the traditional time-domain,a strategy of phase-encoded full waveform inversion (FWI) with random shots allocation is proposed and implemented based on the acceleration of multi-card GPU in the paper.Numerical experiments based on a four-card GPU show that,under the acceleration of phase encoding technique and GPU,the computational efficiency of time-domain FWI can be improved greatly.Meanwhile,by the random allocation of shots in each card,this strategy can suppress the crosstalk more effectively.Thus,the time-domain phase-encoding FWI with random shots allocation has higher inversion precision and convergent rate than the traditional method.

Keywords:full waveform inversion,phase-encoding,crosstalk,multi-card GPU

文章編號：1000-1441(2016)02-0251-10

DOI:10.3969/j.issn.1000-1441.2016.02.011

中圖分類號：P631

文獻(xiàn)標(biāo)識碼：A

基金項(xiàng)目：國家自然科學(xué)基金項(xiàng)目(41574105)資助。

作者簡介：宋鵬(1979—),男,講師,主要從事波動方程成像及反演方法研究工作。通訊作者：王修田(1961—),男,教授,主要從事地震數(shù)據(jù)處理方法研究與軟件系統(tǒng)研發(fā)工作。

收稿日期：2014-12-05;改回日期：2015-08-15。

This research is financially supported by the National Natural Science Foundation of China (Grant No.41574105).