基于多卡GPU集群的多次波逆時(shí)偏移成像技術(shù)

朱　博，宋　鵬，2，李金山，2*，譚　軍，2

（1.中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院，山東青島266100；2.中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島266100）

基于多卡GPU集群的多次波逆時(shí)偏移成像技術(shù)

朱博1，宋鵬1，2，李金山1，2*，譚軍1，2

（1.中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院，山東青島266100；2.中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島266100）

相對(duì)于常規(guī)逆時(shí)偏移，多次波逆時(shí)偏移成像技術(shù)具有更高的成像精度。在分析多次波逆時(shí)偏移成像原理和條件的基礎(chǔ)上，對(duì)多次波逆時(shí)偏移成像中的多次波預(yù)測(cè)、波場(chǎng)延拓等步驟的實(shí)現(xiàn)過(guò)程進(jìn)行了研究。在此基礎(chǔ)上，利用Sigsbee2B地質(zhì)模型，應(yīng)用基于反饋環(huán)理論的自由界面多次波衰減方法預(yù)測(cè)得到純多次波炮集，實(shí)現(xiàn)了基于多卡GPU集群的多次波逆時(shí)偏移成像處理。處理結(jié)果顯示，多次波逆時(shí)偏移的中淺層成像清晰，構(gòu)造的成像精度明顯高于常規(guī)逆時(shí)偏移。多卡GPU集群的應(yīng)用可顯著提高多次波逆時(shí)偏移的計(jì)算速度和效率，使得基于地震大數(shù)據(jù)體進(jìn)行多次波逆時(shí)偏移成像處理成為可能。

多次波波場(chǎng)延拓逆時(shí)偏移成像多卡GPU成像精度

逆時(shí)偏移方法于20世紀(jì)80年代由Whitmore等［1-3］首先提出。其核心為基于雙程波動(dòng)方程，應(yīng)用有限差分等數(shù)值模擬手段來(lái)實(shí)現(xiàn)波場(chǎng)的正時(shí)與逆時(shí)延拓，并通過(guò)正時(shí)波場(chǎng)與逆時(shí)波場(chǎng)的互相關(guān)來(lái)實(shí)現(xiàn)成像。經(jīng)過(guò)幾十年的發(fā)展，逆時(shí)偏移技術(shù)無(wú)論是在理論上還是在實(shí)用性方面都取得了長(zhǎng)足的進(jìn)步。但是，由于常規(guī)逆時(shí)偏移多是基于一次反射波來(lái)實(shí)現(xiàn)成像處理的，其間會(huì)將多次波視為干擾進(jìn)行剔除。而常規(guī)逆時(shí)偏移在一次波無(wú)法到達(dá)的陰影區(qū)或一次波照明度較低的區(qū)域會(huì)遇到成像困難的問(wèn)題。事實(shí)上多次波也是地震波在地下介質(zhì)中傳播形成的反射波。與一次反射波相比，多次波在地下傳播的射線路徑更長(zhǎng)，覆蓋的區(qū)域更廣，地下照明度更加均衡，甚至其能夠照射到一次反射波無(wú)法到達(dá)的陰影區(qū)，從而克服成像困難的難題。近年來(lái)研究人員已經(jīng)對(duì)多次波逆時(shí)偏移成像技術(shù)開(kāi)展了諸多卓有成效的研究工作。當(dāng)前應(yīng)用多次波進(jìn)行逆時(shí)偏移成像的方法主要有2類：第1類方法需先將多次波轉(zhuǎn)化為準(zhǔn)一次波，然后應(yīng)用常規(guī)逆時(shí)偏移算法對(duì)準(zhǔn)一次波進(jìn)行成像［4-6］。此類方法在準(zhǔn)一次波數(shù)據(jù)構(gòu)建過(guò)程中，伴隨著準(zhǔn)一次波有效同相軸的形成，不可避免地會(huì)出現(xiàn)由于不成對(duì)的一次波和多次波互相關(guān)而形成干涉噪音，進(jìn)而降低了準(zhǔn)一次波數(shù)據(jù)和后續(xù)成像結(jié)果的信噪比；第2類方法首先基于原始地震記錄應(yīng)用基于反饋環(huán)理論的自由界面多次波衰減方法（SRME）預(yù)測(cè)得到純多次波記錄，然后以原始地震記錄作為正時(shí)擾動(dòng)，將預(yù)測(cè)得到的多次波記錄作為逆時(shí)擾動(dòng)進(jìn)行偏移成像［7-16］，該類方法不需要求取準(zhǔn)一次波，且其能夠較好的保持偏移剖面的動(dòng)力學(xué)特征。

通常情況下多次波逆時(shí)偏移一次成像計(jì)算相當(dāng)于2次正演模擬的計(jì)算量，對(duì)于大模型數(shù)據(jù)或?qū)嶋H數(shù)據(jù)意味著巨大的計(jì)算消耗，早期的單卡CPU已經(jīng)無(wú)法滿足處理這類海量數(shù)據(jù)的要求。而隨著計(jì)算機(jī)技術(shù)的發(fā)展，多卡GPU集群在處理逆時(shí)偏移時(shí)可獲得幾十倍的加速比。筆者基于多卡GPU集群對(duì)第2類多次波逆時(shí)偏移方法進(jìn)行了研究，以期使得地震大數(shù)據(jù)體的多次波逆時(shí)偏移成像處理成為可能。

1　多次波逆時(shí)偏移成像原理

常規(guī)逆時(shí)偏移成像的原理是基于雙程波動(dòng)方程以震源點(diǎn)S的地震子波作為正時(shí)擾動(dòng)進(jìn)行正演模擬獲得正時(shí)波場(chǎng)，以接收點(diǎn)P接受的一次波地震記錄作為逆時(shí)擾動(dòng)，應(yīng)用有限差分等數(shù)值模擬手段進(jìn)行波場(chǎng)延拓，最后依據(jù)成像條件對(duì)反射點(diǎn)rP實(shí)現(xiàn)成像（圖1a）。多次波逆時(shí)偏移成像則是先以P點(diǎn)處的一次波地震記錄作為正向擾動(dòng)進(jìn)行正演模擬獲得正時(shí)波場(chǎng)，同時(shí)以一階多次波記錄點(diǎn)M1處的記錄作為逆向擾動(dòng)進(jìn)行波場(chǎng)逆時(shí)延拓得到逆時(shí)波場(chǎng)，然后根據(jù)成像條件實(shí)現(xiàn)一階多次波記錄對(duì)反射點(diǎn)的成像；再以M1點(diǎn)處的一階多次波記錄作為正向擾動(dòng)進(jìn)行正演模擬獲得正時(shí)波場(chǎng)，并以M2點(diǎn)處的二階多次波記錄作為逆向擾動(dòng)進(jìn)行波場(chǎng)逆時(shí)延拓得到逆時(shí)波場(chǎng)，實(shí)現(xiàn)二階多次波記錄對(duì)反射點(diǎn)rM2的成像（圖1b）。

圖1　不同類型逆時(shí)偏移成像原理示意Fig.1　Different typesofprinciplesof reverse timemigration imaging

多次波逆時(shí)偏移原理與常規(guī)逆時(shí)偏移相似，也是先進(jìn)行正時(shí)、逆時(shí)波場(chǎng)的延拓，再基于成像條件提取深度成像值。但二者不同的是，常規(guī)逆時(shí)偏移的正時(shí)波場(chǎng)和逆時(shí)波場(chǎng)擾動(dòng)使用的是震源子波和原始地震記錄，而多次波逆時(shí)偏移的正時(shí)波場(chǎng)和逆時(shí)波場(chǎng)擾動(dòng)分別為原始地震記錄和基于原始地震記錄預(yù)測(cè)得到的多次波記錄。

2　多次波逆時(shí)偏移成像的實(shí)現(xiàn)

2.1基于反饋環(huán)理論的多次波預(yù)測(cè)

由多次波逆時(shí)偏移成像的原理可知，要想獲得精確的成像結(jié)果需要先準(zhǔn)確地預(yù)測(cè)出多次波記錄。由于SRME可有效地預(yù)測(cè)出地震記錄中的自由界面多次波記錄，且對(duì)復(fù)雜構(gòu)造區(qū)域的多次波和衰減都可達(dá)到較高的預(yù)測(cè)精度，因此，應(yīng)用該方法進(jìn)行多次波的預(yù)測(cè)。

假設(shè)地震波在地層中傳播過(guò)程的表達(dá)式［17］為

其中

其中

在式（3）中引入不含表層反射的地下響應(yīng)矩陣，表達(dá)式為

為描述檢波器組合特性及虛反射的影響引入算子，則檢波器接收到的地震波場(chǎng)可表示為

其中

由式（7）和式（8）可以分別求得自由界面上檢波器接收到的一次波和多次波。

引入一個(gè)自由界面算子，其表達(dá)式為

則式（6）可表示為

根據(jù)Neumann級(jí)數(shù)性質(zhì)，將式（10）展開(kāi)

對(duì)于多次波逆時(shí)偏移成像而言，通常是以包含一次波和各階多次波的實(shí)際記錄為輸入數(shù)據(jù)，應(yīng)用式（12）和式（13）預(yù)測(cè)得到純多次波記錄，并以此純多次波記錄作為逆時(shí)擾動(dòng)進(jìn)行波場(chǎng)逆時(shí)延拓以得到逆時(shí)波場(chǎng)。

2.2多次波逆時(shí)偏移的波場(chǎng)延拓

多次波逆時(shí)偏移成像主要通過(guò)應(yīng)用聲波方程有限差分?jǐn)?shù)值模擬方法進(jìn)行波場(chǎng)的正時(shí)和逆時(shí)波場(chǎng)延拓。其中，地震波的傳播一般可用二維聲波方程來(lái)表示

式（14）通常應(yīng)用有限差分?jǐn)?shù)值模擬方法來(lái)求解，對(duì)于基于聲波方程的正演模擬來(lái)講，其時(shí)間二階中心差分、空間任意2N階精度差分格式［18］為

基于聲波方程的波場(chǎng)逆時(shí)延拓的時(shí)間二階精度、空間2N階精度差分格式可表達(dá)為

通過(guò)式（16）中的s（i，j，k）可預(yù)測(cè)得到純多次波記錄。

2.3多次波逆時(shí)偏移成像條件分析

多次波逆時(shí)偏移的最終成像需要借助成像條件來(lái)計(jì)算各個(gè)深度點(diǎn)上的成像值。筆者所采用互相關(guān)成像條件［8］的成像值表達(dá)式為要注意的是，和中分別包含了不同階數(shù)的多次波正時(shí)和逆時(shí)波場(chǎng)，即

將式（18）中的表示正時(shí)波場(chǎng)和逆時(shí)波場(chǎng)中的各階多次波代入式（17）中進(jìn)一步展開(kāi)可得

式（19）中第1項(xiàng)各行分別為N（N≥0）階多次波正時(shí)波場(chǎng)（將一次波視為0階多次波）與N+1階多次波逆時(shí)波場(chǎng)的互相關(guān)成像值；第2項(xiàng)分別為N（N≥0）階多次波正時(shí)波場(chǎng)與N+i（i＞1）階多次波逆時(shí)波場(chǎng)的互相關(guān)成像值；第3項(xiàng)分別為N+j（N≥1，j≥0）階多次波正時(shí)波場(chǎng)與N階多次波逆時(shí)波場(chǎng)的互相關(guān)成像值。Liu等［8］認(rèn)為只有第1項(xiàng)能夠?qū)Φ叵碌刭|(zhì)構(gòu)造正確成像，而第2項(xiàng)會(huì)形成串?dāng)_噪音，嚴(yán)重影響多次波逆時(shí)偏移的深層成像精度，第3項(xiàng)則不能成像。

2.4基于多卡CPU集群的多次波逆時(shí)偏移成像

由于多次波逆時(shí)偏移成像的正時(shí)和逆時(shí)2次波場(chǎng)延拓計(jì)算量較大。相比以往的單卡CPU，GPU具有更多的計(jì)算核心和更高的帶寬和訪存速度［19-21］，而且借助高度并行的架構(gòu)可同時(shí)執(zhí)行成千上萬(wàn)個(gè)線程，因此在處理運(yùn)算密集、高度并行、控制簡(jiǎn)單的計(jì)算問(wèn)題時(shí)有著單卡CPU無(wú)可比擬的優(yōu)勢(shì)，二者計(jì)算效率比約為50∶1。因此，基于多卡GPU集群可進(jìn)一步提高多次波逆時(shí)偏移成像的計(jì)算效率和速度。根據(jù)多次波逆時(shí)偏移按炮進(jìn)行偏移計(jì)算的特點(diǎn)，基于處理系統(tǒng)多點(diǎn)接口（MPI）先實(shí)現(xiàn)按炮進(jìn)行任務(wù)劃分，即根據(jù)參與運(yùn)算的GPU卡的數(shù)目啟動(dòng)同樣數(shù)目的處理進(jìn)程，并將所有炮記錄均衡地分配到各個(gè)進(jìn)程，然后將每個(gè)進(jìn)程與各個(gè)卡一對(duì)一進(jìn)行綁定，保證每個(gè)進(jìn)程對(duì)應(yīng)一個(gè)卡，最后將GPU集群中各卡的計(jì)算結(jié)果歸約到主進(jìn)程即可實(shí)現(xiàn)多次波逆時(shí)偏移成像成果的輸出（圖2）。

圖2　基于多卡GPU集群的多次波逆時(shí)偏移成像流程Fig.2　Flow chartof reverse timemigration ofmultiples based on the acceleration ofmulti-card GPU

3　多次波逆時(shí)偏移成像效果

Sigsbee2B地質(zhì)模型是國(guó)際上用于地震多次波數(shù)據(jù)處理的標(biāo)準(zhǔn)模型，依據(jù)該模型建立一個(gè)長(zhǎng)度為24 386m、深度為9 144m，橫向和縱向網(wǎng)格步長(zhǎng)均為7.62m，含有起伏海底界面和高陡鹽丘的構(gòu)造模型（圖3）進(jìn)行逆時(shí)偏移成像處理。

圖3　基于Sigsbee2B地質(zhì)模型的地震多次波數(shù)據(jù)體Fig.3　Seismic data processing ofmultiplesbased on the Sigsbee2Bgeologicalmodel

應(yīng)用SRME方法預(yù)測(cè)得到該地質(zhì)模型原始記錄的純多次波炮集，然后通過(guò)1個(gè)4卡的GPU工作站實(shí)現(xiàn)了數(shù)據(jù)體的多次波逆時(shí)偏移成像處理。為便于分析比較，將地質(zhì)模型的局部構(gòu)造1和2通過(guò)常規(guī)與多次波逆時(shí)偏移對(duì)成像處理的結(jié)果進(jìn)行放大（圖4）對(duì)比可以看出，多次波逆時(shí)偏移處理的中淺層成像清晰，構(gòu)造細(xì)節(jié)分辨率較高，小高速體成像精度明顯高于常規(guī)一次波逆時(shí)偏移。主要原因是因?yàn)槎啻尾▊鞑ヂ窂礁L(zhǎng)，覆蓋范圍更廣，從而使得多次波逆時(shí)偏移在中淺層復(fù)雜構(gòu)造的成像精度優(yōu)于常規(guī)逆時(shí)偏移。

但深層界面的成像結(jié)果也顯示，常規(guī)一次波逆時(shí)偏移在深層的界面清晰，構(gòu)造連續(xù)性好，與原始模型基本吻合；而多次波逆時(shí)偏移的深層構(gòu)造模糊，成像精度相對(duì)較低。造成二者成像精度差異的主要原因是串?dāng)_噪音的影響。因此，要使得多次波逆時(shí)偏移成像效果兼顧各層段，必須設(shè)法壓制多次波逆時(shí)偏移的串?dāng)_噪音，提高其深層成像精度。

圖4　Sigsbee2B地質(zhì)模型不同位置逆時(shí)偏移成像結(jié)果對(duì)比Fig.4　Comparison of imaging resultsof the reverse time migration atdifferent locations in the Sigsbee2Bgeologicalmodel

4　結(jié)論

對(duì)多次波逆時(shí)偏移成像原理進(jìn)行了系統(tǒng)剖析，結(jié)果表明，相對(duì)于一次波，多次波在地下傳播的射線路徑更長(zhǎng)，覆蓋的區(qū)域更廣，地下照明度更加均衡，甚至能夠照明到一次反射波無(wú)法到達(dá)的陰影區(qū)，可以克服常規(guī)逆時(shí)偏移成像困難的難題。

在剖析多次波逆時(shí)偏移成像原理和條件的基礎(chǔ)上，分析多次波逆時(shí)偏移成像中的多次波預(yù)測(cè)、波場(chǎng)延拓等步驟的實(shí)現(xiàn)過(guò)程，利用Sigsbee2B地質(zhì)模型應(yīng)用SRME方法預(yù)測(cè)得到純多次波炮集，實(shí)現(xiàn)了基于多卡GPU集群的多次波逆時(shí)偏移成像處理。處理結(jié)果顯示，多次波逆時(shí)偏移中淺層成像清晰，構(gòu)造分辨率高，小高速體成像精度明顯高于常規(guī)一次波逆時(shí)偏移處理結(jié)果。由于受串?dāng)_噪音的影響，多次波逆時(shí)偏移結(jié)果的深層成像相對(duì)較模糊，成像精度遜色于常規(guī)一次波逆時(shí)偏移，進(jìn)一步研發(fā)高精度的多次波逆時(shí)偏移方法以提高深層的成像精度將是下一步的工作重點(diǎn)。

［1］ Whitemore N D.Iterative depth imaging by backward time propagation［C］.Proceedings of the 53rd Annual International Meeting，1983：382-385.

［2］ Baysal E，Kosloff D D，Sherwood JW.Reverse timemigration［J］. Geophysics，1983，48（11）：1 514-1 524.

［3］ Levin S A.Principle of reverse-time migration［J］.Geophysics，1984，49（5）：581-583.

［4］ Verschuur D J，Berkhout A J.Transformingmultiples into primaries：Experience with field data［C］.Expanded Abstracts of the 75th AnnualMeeting of the Society of Exploration Geophysicists，2005.

［5］ Berkhout A J，Verschuur D J.Imaging ofmultiple reflections［J］. Geophysics，2006，71（4）：I209-I220.

［6］ 單國(guó)健.地表多次波應(yīng)用研究［J］.石油物探，2007，46（6）：604-610. Shan Guojian.Surface-related multiplemigration［J］.Geophysical Prospecting for Petroleum，2007，46（6）：604-610.

［7］ 郭書(shū)娟，李振春，仝兆岐，等.基于廣義的炮偏移方法實(shí)現(xiàn)地表多次波和一次波聯(lián)合成像［J］.地球物理學(xué)報(bào)，2011，54（4）：1 098-1 105. Guo Shujuan，Li Zhenchun，Tong Zhaoqi，et al.Joint imaging of primaries and surface-related multiples based on generalized shot-profilemigration［J］.Chinese Journal of Geophysics，2011，54（4）：1 098-1 105.

［8］ Liu Y，Chang X，Jin D，et al.Reverse timemigration ofmultiplesfor subsaltimaging［J］.Geophysics，2011，76（5）：B209-B216.

［9］ Wang Y，Chang X，Hu H.Simultaneous reverse timemigration of primaries and free-surface relatedmultipleswithoutmultiple prediction［J］.Geophysics，2014，79（1）：S1-S9.

［10］Anstey N A.The Sectional auto-correlogram and the sectional retro-correlogram［J］.Geophysical Prospecting，1966，14（4）：389-426.

［11］Riley DC，Claerbout JF.2-Dmultiple reflections［J］.Geophysics，1976，41（4）：592-620.

［12］Kennett B.The suppression of surfacemultipleson seismic records ［J］.Geophysical Prospecting，1979，27（3）：584-600.

［13］Berkhout A J.Seismic migration：Imaging of acoustic energy by wave field extrapolation［M］.Amsterdam of Holland：Elsevier，1984.

［14］Verschuur D J，BerkhoutA J，Wapenaar C.Adaptive surface-related multiple elimination［J］.Geophysics，1992，57（9）：1 166-1 177.

［15］Verschuur D J，Berkhout A J.Estimation ofmultiple scattering by iterative inversion，Part II：Practical aspects and examples［J］. Geophysics，1997，62（5）：1 596-1 611.

［16］Dedem E V，Verschuur D J.3D surface-related multiple prediction：A sparse inversion approach［J］.Geophysics，2005，70（3）：V31-V43.

［17］譚軍.自由界面多次波的預(yù)測(cè)與衰減［D］.青島：中國(guó)海洋大學(xué)，2011. Tan Jun.The prediction and attenuation of surface-related multiple［D］.Qingdao：Ocean University ofChina，2011.

［18］劉洋，李承楚，牟永光.任意偶數(shù)階精度有限差分法數(shù)值模擬［J］.石油地球物理勘探，1998，33（1）：1-10. Liu Yang，LiChengchu，Mu Yongguang.Finite-difference numerical modeling of any even-order accuracy［J］.Oil Geophysical Prospecting，1998，33（1）：1-10.

［19］李博，劉國(guó)峰，劉洪.地震疊前時(shí)間偏移的一種圖形處理器提速實(shí)現(xiàn)方法［J］.地球物理學(xué)報(bào)，2009，52（1）：245-252. Li Bo，Liu Guofeng，Liu Hong.Amethod of using GPU to accelerate seismic pre-stack timemigration［J］.Chinese Journal of Geophysics，2009，52（1）：245-252.

［20］劉紅偉，李博，劉洪，等.地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn)［J］.地球物理學(xué)報(bào)，2010，53（7）：1 725-1 733. Liu Hongwei，LiBo，Liu Hong，etal.Thealgorithm ofhigh order finite difference pre-stack reverse time migrationm and GPU implementation［J］.Chinese Journal of Geophysics，2010，53（7）：1 725-1 733.

［21］劉文卿，王宇超，雍學(xué)善，等.基于GPU/CPU疊前逆時(shí)偏移研究及應(yīng)用［J］.石油地球物理勘探，2012，47（5）：712-716. Liu Wenqing，Wang Yuchao，Yong Xueshan，et al.Prestack reverse timemigration on GPU/CPU co-parallel computation［J］. OilGeophysical Prospecting，2012，47（5）：712-716.

編輯裴磊

Reverse timem igration ofmultip lesbased on theacceleration ofmulti-card GPU

Zhu Bo1，Song Peng1，2，Li Jinshan1，2，Tan Jun1，2

（1.CollegeofMarineGeosciences，Ocean University ofChina，Qingdao City，Shandong Province，266100，China；2.Key Lab ofSubmarineGeosciencesand Prospecting Techniques，Ministry ofEducation，Ocean University ofChina，Qingdao City，Shandong Province，266100，China）

Compared to the conventional reverse timemigration（RTM），RTM ofmultiples has a higher imaging precision. Based on the analysis of the imaging principle and conditions of RTM ofmultiples，the paper studied the steps ofmultiples prediction，wave field extrapolation and so on in the imaging process of RTM ofmultiples.Then，the shotgatherswith only multipleswere obtained by applying the surface-relatedmultipleseliminationmethod based on feedback loop theory to the Sigsbee2B geologicalmodel，and the imaging processing ofRTM ofmultipleswas implemented based on the acceleration of multi-card GPU.The processing resultsshow that the imaging resultof themiddle-shallow layerobtained by RTM ofmultiples is clear，and the imaging precision of shallow structures is obviously higher than thatobtained by conventional RTM. Besides，the computationalspeed and efficiency can be improved significantly by the application ofmulti-card GPU，which makes itpossible for the imaging processingofRTM ofmultiplesbased on big data.

multiples；wave field extrapolation；reverse timemigration imaging；multi-card GPU；imaging precision

P631.4

A

1009-9603（2015）02-0060-06

2015-01-29。

朱博（1988—），男，湖北隨州人，在讀碩士研究生，從事地震資料數(shù)據(jù)處理研究。聯(lián)系電話：13366640301，E-mail：depacino@sina.com。

李金山（1963—），男，山東昌邑人，教授級(jí)高級(jí)工程師，碩士，從事地震資料數(shù)據(jù)處理研究。聯(lián)系電話：13012556790，E-mail：ljs @ouc.edu.cn。

中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)“多次波分階逆時(shí)偏移成像及其GPU集群加速”（201513005）。