基于GPU加速策略的稀疏反演表面多次波壓制方法

張文武， 劉玉敏， 楊建華， 張晟瑞， 王 程， 白玉勝， 石 穎

( 1. 東北石油大學(xué) 地球科學(xué)學(xué)院，黑龍江 大慶 163318； 2. 東北石油大學(xué) 電氣信息工程學(xué)院，黑龍江 大慶 163318； 3. 大慶鉆探工程公司 地球物理勘探一公司，黑龍江 大慶 163357； 4. 中國石油集團(tuán)測(cè)井有限公司 生產(chǎn)測(cè)井中心，陜西 西安 710201 )

0 引言

常規(guī)的表面多次波壓制方法有兩種[1]：濾波法和基于波動(dòng)方程的預(yù)測(cè)相減法。濾波法利用一次波與多次波的可分離性壓制多次波[2-6]，對(duì)于復(fù)雜構(gòu)造的地震數(shù)據(jù)，一次波和多次波的可分離性低，不能有效保護(hù)有效信號(hào)?；诓▌?dòng)方程的預(yù)測(cè)相減法以地震物理學(xué)為基礎(chǔ)，包含多次波預(yù)測(cè)和自適應(yīng)相減兩部分，能夠適應(yīng)復(fù)雜地質(zhì)構(gòu)造[7-9]，但是自適應(yīng)相減部分存在傷害有效信號(hào)的風(fēng)險(xiǎn)[10-13]。

稀疏反演法建立在與反饋迭代法相同的理論基礎(chǔ)上，通過迭代反演方式估計(jì)震源子波和地下一次波脈沖響應(yīng)，進(jìn)而重構(gòu)一次波和多次波[14]，省去自適應(yīng)相減部分，可以更好地保護(hù)有效信號(hào)。Vanborselen R G等提出參數(shù)反演方法，在震源子波已知和能量最小的約束條件下估計(jì)一次波[15]；Biersteker J假設(shè)一次波能量最小，利用二維反褶積法估計(jì)缺失的淺層有效信號(hào)[16]；Amundsen L利用上、下行波場的多維除法估計(jì)一次波，雖然不需要已知震源子波信息，但需要直達(dá)波信息，特別是近偏移距的直達(dá)波信息[17]；Wang Y利用基于顯式矩陣的反演方法估計(jì)有效波，可以消除由近偏移距缺失導(dǎo)致的壓制效果不理想，但包括多次波預(yù)測(cè)和自適應(yīng)相減兩部分[18]；Berkhout A J在逆數(shù)據(jù)域估計(jì)一次波和多次波，省去自適應(yīng)相減部分[19]，但大量的矩陣求逆運(yùn)算存在反演風(fēng)險(xiǎn)[20]；Herrmann F J等在曲波域估計(jì)一次波，需要預(yù)測(cè)的多次波指導(dǎo)分離過程[21]。在SRME(Surface-Related Multiple Elimination)理論基礎(chǔ)上，Groenestijn G J A V等提出基于稀疏反演估計(jì)一次波的方法(EPSI, Estimation of Primaries by Sparse Inversion)，可以避免對(duì)有效波的損傷[14,22]，但高計(jì)算成本是制約EPSI方法發(fā)展的主要因素。

基于CUDA架構(gòu)的GPU加速技術(shù)在處理矩陣運(yùn)算方面優(yōu)勢(shì)明顯，廣泛應(yīng)用于地震數(shù)據(jù)處理方面[23-25]。為解決EPSI方法大量的矩陣運(yùn)算，筆者在傳統(tǒng)的GPU加速方法的基礎(chǔ)上進(jìn)行改進(jìn)，實(shí)現(xiàn)所有頻率成分的同步運(yùn)算，并減少主機(jī)端和設(shè)備端的數(shù)據(jù)傳輸，相比傳統(tǒng)GPU加速的EPSI方法和常規(guī)EPSI方法，計(jì)算效率更高。

1 EPSI方法原理

根據(jù)Rayleigh Ⅱ積分和數(shù)據(jù)矩陣的概念[9]，接收到的上行總波場的單頻表達(dá)形式為

P=X0S+X0R-P,

(1)

式中：P為接收的上行總波場；X0為一次波脈沖響應(yīng)；S為震源特性；R-為地表反射系數(shù)。對(duì)于海洋數(shù)據(jù)，R-=-1，X0和S的乘積為有效波P0(不在自由表面發(fā)生反射的波，包括一次波和層間多次波)。X0、R-和P三者的乘積為自由表面多次波M。不同于SRME方法，EPSI方法不引入自適應(yīng)算子A，直接把X0和S作為未知量，通過不斷迭代反演估計(jì)兩者的值。盡管假設(shè)反射系數(shù)R-=-1，但式(1)的未知量(X0和S)多于已知量(P)，為欠定方程，有多解，需增加約束條件：(1)一次波脈沖響應(yīng)X0在時(shí)間域是稀疏的；(2)忽略震源的方向特性。

式(1)可看作一個(gè)反演過程，通過依次求取X0和S，使接收的地震記錄與估計(jì)的地震記錄的殘差能量最小。目標(biāo)函數(shù)[26]可表示為

(2)

利用最速下降法，經(jīng)過多次迭代得到X0和S。由于X0和S相互獨(dú)立，整個(gè)過程可分為兩部分：(1)固定S，更新X0；(2)固定X0，更新S。EPSI方法的全部過程為

(1)設(shè)X0和S的初值為0。

(3)對(duì)ΔX0進(jìn)行稀疏約束：

2)變換到時(shí)間域，采用合適的時(shí)窗(隨迭代次數(shù)的增加而增大，可以加快收斂速度)，選取每道振幅的最大值。

2 基于GPU加速策略的EPSI方法

EPSI方法存在大量矩陣運(yùn)算，耗費(fèi)大量時(shí)間。GPU加速技術(shù)在處理矩陣這類并行度高、密集型強(qiáng)的浮點(diǎn)方面運(yùn)算優(yōu)勢(shì)明顯，因此被廣泛應(yīng)用于地震數(shù)據(jù)處理[23-25]。將GPU加速策略應(yīng)用于EPSI方法，可

圖1 傳統(tǒng)GPU加速方法的計(jì)算流程Fig.1 Flow chart of traditional acceleration strategy

以提高計(jì)算效率，降低計(jì)算成本。

CUDA是NVIDIA推出的通用并行計(jì)算架構(gòu)，作為一種運(yùn)算平臺(tái)，可以用C語言寫出在顯示芯片上執(zhí)行的程序，不需要學(xué)習(xí)特定的顯示芯片指令和特殊結(jié)構(gòu)。GPU加速方法的過程為主機(jī)端(CPU及系統(tǒng)的內(nèi)存)的數(shù)據(jù)導(dǎo)入設(shè)備端(GPU及其顯示內(nèi)存)，在設(shè)備端運(yùn)行核函數(shù)[27]，將運(yùn)行的結(jié)果導(dǎo)回主機(jī)端。

在傳統(tǒng)的GPU加速方法基礎(chǔ)上，提出一種新的GPU加速策略。

傳統(tǒng)GPU加速方法的計(jì)算流程為：(1)將x-t域地震數(shù)據(jù)在時(shí)間方向進(jìn)行傅里葉變換，得到x-ω域數(shù)據(jù)。(2)抽取單頻成分?jǐn)?shù)據(jù)導(dǎo)入設(shè)備端。(3)在設(shè)備端完成矩陣運(yùn)算。(4)將運(yùn)算結(jié)果導(dǎo)回到主機(jī)端。(5)進(jìn)行下一頻率成分?jǐn)?shù)據(jù)的運(yùn)算，重復(fù)步驟(2-4)，直至所有的頻率成分計(jì)算結(jié)束。以求取一次波脈沖響應(yīng)X0的搜索步長為例，傳統(tǒng)GPU加速方法的計(jì)算流程見圖1。其中，有_SF的表示抽取單頻后的數(shù)據(jù)，有_d的表示設(shè)備端的數(shù)據(jù)，ω0表示選取的初始頻率，ωmax表示截止頻率。

EPSI方法的不同頻率成分之間的矩陣運(yùn)算相互獨(dú)立、互不相關(guān)，因此在理論上可以同步進(jìn)行?；趥鹘y(tǒng)GPU加速策略的EPSI方法，沒有實(shí)現(xiàn)不同頻率成分矩陣運(yùn)算的同步進(jìn)行，除初始頻率成分的計(jì)算外，其他頻率成分的計(jì)算需要等待之前的所有頻率成分的計(jì)算結(jié)束后才能進(jìn)行，在很大程度上影響算法的計(jì)算效率。另外，得到的所有的中間結(jié)果需要導(dǎo)回主機(jī)端，若它們?cè)俅螀⑴c計(jì)算，仍然需要從主機(jī)端導(dǎo)入設(shè)備端，中間結(jié)果在主機(jī)端和設(shè)備端之間傳導(dǎo)耗費(fèi)大量的時(shí)間。

圖2 新的GPU加速策略的計(jì)算流程Fig.2 Flow chart of new acceleration strategy

在對(duì)傳統(tǒng)的GPU加速方法改進(jìn)后，提出一種新的GPU加速策略，以求取一次波脈沖響應(yīng)X0的搜索步長為例，計(jì)算流程見圖2。該策略相比于傳統(tǒng)GPU加速方法的優(yōu)點(diǎn)：(1)所有頻率成分整體導(dǎo)入的耗時(shí)較抽取單頻導(dǎo)入耗時(shí)少。(2)可實(shí)現(xiàn)多個(gè)頻率成分運(yùn)算的同步進(jìn)行。(3)得到的中間結(jié)果除部分必須導(dǎo)回主機(jī)端完成運(yùn)算外，其他的中間結(jié)果直接存儲(chǔ)于設(shè)備端，減少數(shù)據(jù)在主機(jī)端和設(shè)備端之間的傳導(dǎo)。由于設(shè)備端的存儲(chǔ)空間有限，當(dāng)數(shù)據(jù)體較大、設(shè)備端的存儲(chǔ)空間不能保存所有中間變量時(shí)，需引入多卡GPU的編程策略，處理大的數(shù)據(jù)體分塊，每個(gè)卡完成部分?jǐn)?shù)據(jù)的運(yùn)算，減輕單卡的存儲(chǔ)負(fù)擔(dān)，提高新的GPU加速策略的適用性[26]。

3 數(shù)據(jù)測(cè)試

3.1 水平層狀模型測(cè)試

實(shí)驗(yàn)?zāi)Ｐ蜑閱螌铀綄訝罱橘|(zhì)，正演得到的地震記錄有200炮，每炮200道接收，空間方向采樣間隔為10.0 m，時(shí)間方向1 000個(gè)采樣點(diǎn)，采樣間隔為4 ms，包含兩階多次波。采用EPSI方法，得到不同迭代次數(shù)估計(jì)的一次波單炮記錄(見圖3)。由圖3可知，隨迭代次數(shù)增加，一次波同相軸由淺到深被估計(jì)出來。抽取該單炮記錄的某一道進(jìn)行對(duì)比分析(見圖4)，EPSI方法估計(jì)的一次波和多次波的曲線與真實(shí)的曲線幾乎重合。

圖3 水平層狀模型不同迭代次數(shù)估計(jì)的一次波的單炮記錄Fig.3 Single-shot record of estimated primary of horizontal layered by different iterations

3.2 SMAART模型測(cè)試

對(duì)構(gòu)造相對(duì)復(fù)雜的SMAART模型進(jìn)行測(cè)試。地震記錄有200炮，每炮200道接收，空間方向采樣間隔為22.9 m，時(shí)間方向626個(gè)采樣點(diǎn)，采樣間隔為8 ms。采用EPSI方法，得到不同迭代次數(shù)的單炮記錄(見圖5)。由圖5可知，隨迭代次數(shù)增加，一次波的同相軸由淺到深被估計(jì)出來，擬合殘差逐漸減小。抽取該單炮記錄的某一道進(jìn)行分析(見圖6)，EPSI方法估計(jì)的一次波和多次波的曲線與真實(shí)的曲線幾乎重合。

圖4 水平層狀模型單道記錄Fig.4 Single-trace record of horizontal layer model

圖5 SMAART模型不同迭代次數(shù)估計(jì)的一次波的單炮記錄Fig.5 Single-shot record of SMAART model in different iterations

圖6 SMAART模型單道記錄Fig.6 Single-trace record of SMAART model

EPSI方法壓制多次波前后的零偏移剖面的結(jié)果見圖7，其中箭頭標(biāo)記處表明EPSI方法能夠有效的壓制多次波。EPSI方法與SRME方法壓制多次波后的結(jié)果的局部放大(2.2～2.6 s)見圖8。由圖8可知，采用EPSI方法壓制多次波后同相軸更加連續(xù)，可以更好地保護(hù)有效信號(hào)。

圖7 EPSI方法壓制多次波前后零偏移距剖面Fig.7 Zero-offset profiles before and after multiple suppression by the EPSI method

圖8 EPSI和SRME方法壓制多次波后的零偏移距剖面的局部放大Fig.8 Zoomed-in sections of zero Offset profile after multiple supression by EPSI and SRME methods

3.3 加速比測(cè)試

分別用常規(guī)EPSI方法、基于傳統(tǒng)GPU加速的EPSI方法和基于新GPU加速策略的EPSI方法，處理SMAART模型的地震記錄，迭代30次計(jì)算耗時(shí)見圖9。其中：CPU表示常規(guī)EPSI方法；GPU_Global表示存儲(chǔ)器為全局存儲(chǔ)器的、基于傳統(tǒng)GPU加速的EPSI方法；GPU_Share表示存儲(chǔ)器為共享存儲(chǔ)器的、基于傳統(tǒng)GPU加速的EPSI方法；GPU_New表示基于新GPU加速策略的EPSI方法。由圖9可知，基于新GPU加速策略的EPSI方法的計(jì)算效率較常規(guī)EPSI方法的計(jì)算效率提升11倍，較基于傳統(tǒng)GPU加速的EPSI方法提升4倍；不同存儲(chǔ)器類型的、基于傳統(tǒng)GPU加速的EPSI方法的運(yùn)行時(shí)間幾乎相等。通常情況下，存儲(chǔ)器為共享存儲(chǔ)器的運(yùn)算速度是存儲(chǔ)器為全局存儲(chǔ)器的運(yùn)算速度的2倍，由于數(shù)據(jù)傳輸?shù)臅r(shí)間遠(yuǎn)遠(yuǎn)大于在設(shè)備端計(jì)算的時(shí)間，因此兩種方法總計(jì)算耗時(shí)相似。

圖9 不同方法的計(jì)算耗時(shí)Fig.9 Time cost of different methods

以SMAART模型為例，根據(jù)基于新GPU加速策略的EPSI方法和常規(guī)EPSI方法的運(yùn)行結(jié)果，驗(yàn)證前者正確性。兩種方法估計(jì)的一次波單炮記錄見圖10。由圖10可知，兩種方法估計(jì)一次波的同相軸形態(tài)和位置幾乎一致。兩種方法估計(jì)的一次波頻譜見圖11。由圖11可知，兩者頻譜曲線幾乎重合。兩種方法估計(jì)的一次波單道見圖12。由圖12可知，兩種方法的單道曲線形態(tài)幾乎一致，表明新GPU加速策略的EPSI方法正確有效。

圖10 兩種方法估計(jì)的一次波單炮記錄Fig.10 Single-shot record of estimated primary by two methods

圖11 兩種方法估計(jì)的一次波頻譜Fig.11 Spectrum comparsion of the estimated primary by two methods

圖12 兩種方法估計(jì)的一次波單道記錄Fig.12 Single-trace record of the estimated primary by two methods

4 結(jié)論

(1)基于新GPU加速策略的EPSI方法將所有頻率成分整體導(dǎo)入設(shè)備端，并將中間結(jié)果存儲(chǔ)于設(shè)備端，實(shí)現(xiàn)所有頻率成分的同步運(yùn)算，減少主機(jī)端和設(shè)備端的數(shù)據(jù)傳輸時(shí)間，相比于傳統(tǒng)GPU加速的EPSI方法，計(jì)算效率更高。

(2)SMAART模型測(cè)試結(jié)果表明，基于新GPU加速策略的EPSI方法的計(jì)算結(jié)果正確有效，可以應(yīng)用于地震資料的表面多次波的壓制。