石 穎，陸加敏，柯 璇，田東升，王 菲

（1.東北石油大學(xué) 地球科學(xué)學(xué)院，黑龍江 大慶 163318； 2.大慶油田有限責(zé)任公司 勘探開(kāi)發(fā)研究院，黑龍江 大慶163712； 3.中國(guó)石油大學(xué)（北京）地球科學(xué)學(xué)院，北京 102200）

基于GPU并行加速的疊前逆時(shí)偏移方法

石 穎1，陸加敏2，柯 璇1，田東升1，王 菲3

（1.東北石油大學(xué) 地球科學(xué)學(xué)院，黑龍江 大慶 163318； 2.大慶油田有限責(zé)任公司 勘探開(kāi)發(fā)研究院，黑龍江 大慶163712； 3.中國(guó)石油大學(xué)（北京）地球科學(xué)學(xué)院，北京 102200）

為了提高復(fù)雜地下介質(zhì)的成像精度和偏移算法的計(jì)算效率，提出可高效對(duì)地下復(fù)雜構(gòu)造進(jìn)行準(zhǔn)確成像的GPU加速疊前逆時(shí)偏移方法.該方法采用雙程聲波方程進(jìn)行波場(chǎng)延拓，突破傾角限制，借助于高階有限差分方法實(shí)現(xiàn)疊前逆時(shí)偏移成像；利用GPU（Graphic Processing Unit）并行加速技術(shù)對(duì)波場(chǎng)延拓和成像進(jìn)行計(jì)算，相比于傳統(tǒng)算法，其計(jì)算效率有較大提高，可以解決疊前逆時(shí)偏移算法計(jì)算量過(guò)大問(wèn)題；在獲取波場(chǎng)信息過(guò)程中，也采用隨機(jī)邊界條件，實(shí)施以計(jì)算換存儲(chǔ)策略，解決逆時(shí)偏移計(jì)算中的海量存儲(chǔ)問(wèn)題.模型測(cè)試結(jié)果表明，該方法能夠高效和高精度地對(duì)地下復(fù)雜地質(zhì)體成像.

逆時(shí)偏移；GPU；加速；高階有限差分；隨機(jī)邊界條件；復(fù)雜構(gòu)造

0 引言

隨著勘探目標(biāo)的日趨復(fù)雜化，疊前時(shí)間偏移[1]無(wú)法完全滿足地震勘探對(duì)成像技術(shù)的需求.在疊前深度偏移方法中，疊前逆時(shí)偏移作為一種高精度地震成像方法，可以對(duì)回轉(zhuǎn)波、棱柱波和多次波等成像，突破傾角的限制，能夠?qū)?fù)雜地質(zhì)構(gòu)造進(jìn)行精確成像.逆時(shí)偏移應(yīng)用于地震勘探領(lǐng)域始于20世紀(jì)80年代[2-5]，受到當(dāng)時(shí)硬件計(jì)算和存儲(chǔ)條件約束影響，并未得到充分發(fā)展.近年來(lái)，隨著低成本的并行計(jì)算設(shè)備和有效存儲(chǔ)硬件的出現(xiàn)，逆時(shí)偏移在地震成像領(lǐng)域獲得新的發(fā)展，逐漸成為主要的成像復(fù)雜構(gòu)造的疊前深度偏移方法.然而，計(jì)算成本高、存儲(chǔ)量大以及低頻噪音[6]問(wèn)題是困擾逆時(shí)偏移方法發(fā)展的瓶頸.劉紅偉[7]和李博[8]等通過(guò)有限差分算法及GPU[9]實(shí)現(xiàn)逆時(shí)偏移，提高逆時(shí)偏移算法的計(jì)算效率.Robert G C[10]提出隨機(jī)邊界的思想，使得邊界反射波以隨機(jī)噪音的形式成像，進(jìn)而可以通過(guò)波場(chǎng)反向傳播計(jì)算模擬各個(gè)時(shí)刻的波場(chǎng)，該方法在很大程度上節(jié)約了存儲(chǔ)量.Zhang Y[11]和劉紅偉[12]利用拉普拉斯算子濾波方法消除低頻成像噪音.

筆者利用GPU/CPU協(xié)同并行加速實(shí)現(xiàn)逆時(shí)偏移算法中的波場(chǎng)延拓和相關(guān)成像計(jì)算，與傳統(tǒng)方法相比，能夠大幅度提高逆時(shí)偏移的計(jì)算效率，縮短計(jì)算時(shí)間，節(jié)省地震偏移成像的處理周期，可以解決疊前逆時(shí)偏移計(jì)算量巨大的問(wèn)題；利用隨機(jī)邊界解決存儲(chǔ)問(wèn)題，利用拉普拉斯算子濾波方法去除低頻噪音，使得疊前逆時(shí)偏移技術(shù)在工業(yè)領(lǐng)域得以廣泛應(yīng)用成為可能.

1 逆時(shí)偏移原理

疊前逆時(shí)偏移算法的實(shí)現(xiàn)步驟包括震源波場(chǎng)沿時(shí)間方向的正向延拓、檢波點(diǎn)波場(chǎng)沿時(shí)間方向的反向延拓和正確應(yīng)用成像條件.在反向延拓中，需要計(jì)算從最大時(shí)刻到零時(shí)刻的所有時(shí)間點(diǎn)上地下介質(zhì)的波場(chǎng).逆時(shí)偏移的成像條件通常包括激發(fā)時(shí)刻成像條件、互相關(guān)成像條件及振幅比值成像條件.較為常用的互相關(guān)成像條件需要使用在同一時(shí)刻的震源波場(chǎng)和檢波點(diǎn)波場(chǎng).由于前者是正傳波場(chǎng)，后者是反傳波場(chǎng)；若想同時(shí)得到相同時(shí)刻的2個(gè)波場(chǎng)，則必須存儲(chǔ)其中一個(gè)波場(chǎng)的整個(gè)傳播過(guò)程的波場(chǎng)信息，即每一時(shí)刻的波場(chǎng)分布，需要消耗甚大的存儲(chǔ)資源，這在實(shí)際操作中是難以滿足的.為此，在疊前逆時(shí)深度偏移計(jì)算中采用隨機(jī)邊界條件，只需存儲(chǔ)最大2個(gè)時(shí)刻的震源波場(chǎng)，通過(guò)計(jì)算獲得其他時(shí)刻的波場(chǎng)信息，然后將震源波場(chǎng)和檢波點(diǎn)波場(chǎng)同時(shí)反傳，再進(jìn)行相關(guān)成像，無(wú)需存儲(chǔ)震源波場(chǎng)正向傳播的中間結(jié)果，避免對(duì)巨大存儲(chǔ)量的需求，是一種以計(jì)算換存儲(chǔ)、以時(shí)間換空間的策略.因此，該方法在利用隨機(jī)邊界解決存儲(chǔ)問(wèn)題的同時(shí)，對(duì)計(jì)算能力也提出更高的要求.

疊前逆時(shí)偏移算法通過(guò)直接求解雙程聲波方程獲取波場(chǎng)傳播信息，從聲波方程出發(fā)，利用高階有限差分算法，計(jì)算震源波場(chǎng)正傳和檢波點(diǎn)波場(chǎng)反傳的地震波場(chǎng)信息.二維均勻橫向各向同性介質(zhì)聲波方程為

式中：V為速度，是空間變量的函數(shù)；U是位移函數(shù)；t為時(shí)間坐標(biāo)；x，z為空間坐標(biāo).

利用規(guī)則網(wǎng)格的高階有限差分求解聲波方程，其中U（x，z，t）對(duì)x，z的二階偏導(dǎo)數(shù)的2 N階精度中心差分格式分別近似為

U（x，z，t）對(duì)時(shí)間t的二階精度中心差分格式為

式（2-4）中：Uik，j＝U（iΔx，jΔz，nΔt），Δx，Δz分別為沿x，z方向的空間采樣間隔，Δt為時(shí)間步長(zhǎng)；Cn為差分系數(shù)，且不同階數(shù)的Cn的求解公式為

綜合式（1-5），波場(chǎng)正向和反向延拓的公式分別表示為

式中：v為當(dāng)前時(shí)刻點(diǎn)所對(duì)應(yīng)的速度.利用式（6-7），結(jié)合隨機(jī)邊界條件，進(jìn)行波場(chǎng)的正向與逆向延拓計(jì)算，以避免因存儲(chǔ)每一時(shí)刻的波場(chǎng)信息而占用巨大存儲(chǔ)量.

完成震源波場(chǎng)的正向延拓和檢波點(diǎn)波場(chǎng)的反向延拓后，下步計(jì)算的關(guān)鍵是選擇合適的成像條件，互相關(guān)成像條件可為成像提供準(zhǔn)確的動(dòng)力學(xué)特性，且實(shí)現(xiàn)方法簡(jiǎn)單.對(duì)于高陡界面，逆時(shí)偏移結(jié)果產(chǎn)生大量強(qiáng)振幅的低頻噪音，從而降低地震數(shù)據(jù)體的成像精度，需采用適當(dāng)?shù)姆椒▔褐瞥上駭?shù)據(jù)中產(chǎn)生的低頻噪音.

2 協(xié)同并行加速逆時(shí)偏移算法

2.1 GPU加速技術(shù)

GPU即圖形處理單元，其功能已不局限于圖形渲染，隨著GPU的可編程性不斷提高，應(yīng)用范圍愈加廣泛.由于GPU擁有比CPU更加強(qiáng)大的并行計(jì)算能力，為許多領(lǐng)域的科學(xué)計(jì)算提供了新的選擇.與CPU相比，GPU設(shè)計(jì)更多的晶體管用于數(shù)據(jù)處理或執(zhí)行單元[9]，也引入片內(nèi)共享存儲(chǔ)器，極大地提高計(jì)算效率.

運(yùn)行在GPU上的程序稱為kernel（內(nèi)核函數(shù)）.一個(gè)kernel函數(shù)中存在2個(gè)層次的并行，即Grid中block間并行和block中的thread間并行[9].

GPU的編程平臺(tái)是CUDA（統(tǒng)一計(jì)算設(shè)備架構(gòu)），它包括一個(gè)硬件驅(qū)動(dòng)程序和一個(gè)應(yīng)用程序接口（API）[13]，降低編程的難度，能夠通過(guò)函數(shù)調(diào)用訪問(wèn)顯存和在GPU上執(zhí)行指令，進(jìn)行大規(guī)模的并行計(jì)算，這些改進(jìn)使CUDA架構(gòu)更加適合進(jìn)行GPU通用計(jì)算.

2.2 實(shí)現(xiàn)方法

疊前逆時(shí)偏移算法的核心由震源波場(chǎng)正傳、檢波點(diǎn)波場(chǎng)反傳和應(yīng)用相關(guān)成像條件組成，均是CPU串行計(jì)算最為耗時(shí)的部分.采用基于CUDA架構(gòu)的GPU和CPU聯(lián)合并行計(jì)算方法，對(duì)波場(chǎng)傳播和相關(guān)成像進(jìn)行加速，先將激發(fā)點(diǎn)的初始波場(chǎng)值和檢波點(diǎn)的最大時(shí)間波場(chǎng)值由內(nèi)部存儲(chǔ)器（內(nèi)存）傳至設(shè)備存儲(chǔ)器（顯存）中，以避免由多次重復(fù)對(duì)內(nèi)部存儲(chǔ)器的數(shù)據(jù)讀寫(xiě)所引起的時(shí)間延遲；然后把GPU的多核處理器劃分為相應(yīng)個(gè)數(shù)的計(jì)算塊（block），同時(shí)每個(gè)計(jì)算塊又可以劃分為若干個(gè)線程（thread），為了使GPU更高效運(yùn)行，定義每個(gè)計(jì)算塊中分配256個(gè)線程，利用GPU的多線程實(shí)現(xiàn)大規(guī)模的并行計(jì)算，在計(jì)算過(guò)程中涉及到的數(shù)據(jù)讀寫(xiě)和運(yùn)算均在設(shè)備存儲(chǔ)器和圖形處理器（GPU）中進(jìn)行，將波場(chǎng)延拓及相關(guān)成像的計(jì)算并行化，提高計(jì)算效率；最后將數(shù)據(jù)傳回內(nèi)部存儲(chǔ)器，通過(guò)CPU進(jìn)行結(jié)果的I/O操作.疊前逆時(shí)偏移主要流程見(jiàn)圖1.

圖1 GPU加速疊前逆時(shí)偏移流程

利用相關(guān)成像條件計(jì)算疊前逆時(shí)偏移算法時(shí)，需要同一時(shí)刻正向傳播的震源波場(chǎng)和反向傳播的檢波點(diǎn)波場(chǎng)，因此在利用相關(guān)成像條件之前，需要事先存儲(chǔ)震源波場(chǎng)或檢波點(diǎn)波場(chǎng).為了節(jié)約存儲(chǔ)空間，文中采用隨機(jī)邊界條件的方法，首先利用GPU加速震源波場(chǎng)正傳，只保留最大時(shí)刻的波場(chǎng)，而不保留中間時(shí)刻的波場(chǎng)；然后利用隨機(jī)邊界條件，由最大時(shí)刻的波場(chǎng)反向傳播震源波場(chǎng)，在傳播的每一時(shí)刻與反傳的檢波點(diǎn)波場(chǎng)進(jìn)行相關(guān)成像.該方法在增加少量計(jì)算量的前提下，較好地解決逆時(shí)偏移成像中的存儲(chǔ)問(wèn)題，在GPU加速技術(shù)大幅度降低計(jì)算成本的情況下，增加少量的計(jì)算量對(duì)算法的計(jì)算效率沒(méi)有太大影響.

3 模型試算

對(duì)Marmousi模型進(jìn)行疊前逆時(shí)深度偏移測(cè)試計(jì)算，所采用的平臺(tái)為GPU Nvidia Geforce GTX560，顯存為1 024 M，顯存位寬為256 bit，核心頻率為850 MHz，顯存頻率為4 500 MHz，流處理器為336個(gè).另一計(jì)算環(huán)境為Intel i3的CPU，內(nèi)存為8 G.

Marmousi模型主要由傾斜地層、高角度逆沖斷層、角度不整合地層及地層隆起構(gòu)成，速度模型見(jiàn)圖2.模型大小為737×750，其中采樣點(diǎn)為750，網(wǎng)格大小為12.5 m×4.0 m，介質(zhì)速度為1 500～5 500 m/s，震源采用Ricker子波，共240炮，每炮96道.

圖2 Marmousi速度模型

利用文中疊前逆時(shí)偏移并行加速算法計(jì)算Marmousi模型數(shù)據(jù)，隨機(jī)抽取第70，100，150，200炮的GPU單炮逆時(shí)偏移結(jié)果（見(jiàn)圖3）.將單炮的偏移結(jié)果進(jìn)行疊加，即可得到Marmousi模型成像剖面.

圖3 隨機(jī)抽取單炮疊前逆時(shí)偏移結(jié)果

與單程波偏移結(jié)果不同，由于存在強(qiáng)反射界面內(nèi)的反射波，逆時(shí)偏移結(jié)果產(chǎn)生低頻噪音.常見(jiàn)的壓制逆時(shí)偏移低頻噪音的方法有高通濾波方法、坡印廷矢量成像條件和拉普拉斯算子濾波法等，文中采用拉普拉斯算子濾波法對(duì)逆時(shí)偏移結(jié)果去噪，去噪后的偏移結(jié)果見(jiàn)圖4.由圖4可見(jiàn)，同相軸清晰，高陡構(gòu)造成像精度較高.

圖4 Marmousi模型疊前逆時(shí)偏移結(jié)果

計(jì)算成本高也是阻礙逆時(shí)偏移方法發(fā)展的主要瓶頸之一，GPU加速技術(shù)可以大幅度地提高逆時(shí)偏移方法的計(jì)算效率，利用CPU進(jìn)行Marmousi模型單炮偏移耗時(shí)53 min，而GPU耗時(shí)37 s；完成Marmousi模型240炮數(shù)據(jù)偏移，CPU耗時(shí)12 720 min，GPU耗時(shí)8 800 s，計(jì)算效率提高約86倍.可見(jiàn)，在逆時(shí)偏移算法中引入GPU并行加速技術(shù)，是提高逆時(shí)偏移計(jì)算效率的有效途徑.

4 結(jié)論

（1）針對(duì)疊前逆時(shí)偏移算法計(jì)算成本高的問(wèn)題，利用CPU/GPU協(xié)同并行計(jì)算，由GPU負(fù)責(zé)計(jì)算量較大的波場(chǎng)延拓和成像運(yùn)算，能夠提高逆時(shí)偏移的運(yùn)算效率.相比于傳統(tǒng)的CPU串行計(jì)算方法，GPU并行加速計(jì)算使計(jì)算效率提高86倍左右.

（2）疊前逆時(shí)深度偏移算法計(jì)算精度較高，可以對(duì)復(fù)雜構(gòu)造地震勘探數(shù)據(jù)準(zhǔn)確成像；借助于GPU加速計(jì)算的優(yōu)勢(shì)，隨機(jī)邊界方法可有效解決數(shù)據(jù)存儲(chǔ)問(wèn)題.

（3）在提高算法計(jì)算精度的同時(shí)，根據(jù)GPU硬件的結(jié)構(gòu)特點(diǎn)改進(jìn)并行算法，優(yōu)化存儲(chǔ)器訪問(wèn)，使數(shù)據(jù)傳輸?shù)耐瑫r(shí)，GPU也能夠正常進(jìn)行運(yùn)算，隱藏?cái)?shù)據(jù)傳輸及訪問(wèn)延遲，是未來(lái)的研究方向.

[1]袁剛，蔣波，曾華會(huì).各向異性疊前時(shí)間偏移在塔里木碳酸鹽巖資料處理中的應(yīng)用[J].大慶石油學(xué)院學(xué)報(bào)，2010，34（3）：23-28.

[2]Whitmore D N.Iterative depth imaging by back time propagation[C].53th Annual International Meeting，SEG Expanded Abstracts，1983：382-385.

[3]Baysal E，Kosloff D D，Sherwood J W C.Reverse time migration[J].Geophysics，1983，48（11）：1514-1524.

[4]Mc Mechan G A.Migration by exploration of time-dependent boundary values[J].Geophysical Prospecting，1983，31：413-420.

[5]Loewenthal D，Stoffa P A，F(xiàn)aria E L.Suppressing the unwanted reflections of the full wave equation[J].Geophysics，1987，52（7）：1007-1012.

[6]陳康，吳國(guó)忱.逆時(shí)偏移拉普拉斯算子濾波改進(jìn)算法[J].石油地球物理勘探，2012，47（2）：249-255.

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

[8]李博，劉紅偉，劉國(guó)峰，等.地震疊前逆時(shí)偏移算法的CPU/GPU 實(shí)施對(duì)策[J].地球物理學(xué)報(bào)，2010，53（12）：2938-2943.

[9]張舒，褚艷利，趙開(kāi)勇，等.GPU高性能運(yùn)算之CUDA[M].北京：中國(guó)水利水電出版社，2009.

[10]Robert G C.Reverse time migration with random boundaries[C].79th Annual International Meeting，SEG Expanded Abstracts，2009：2809-2813.

[11]Zhang Y，Sun James.Practical issues in reverse time migration：true amplitude gathers，noise removal and harmonic source encoding[J].First Break，2009，1：53-59.

[12]劉紅偉，劉洪，鄒振.地震疊前逆時(shí)偏移中的去噪與存儲(chǔ)[J].地球物理學(xué)報(bào)，2010，53（9）：2171-2180.

[13]李博，劉國(guó)峰，劉洪.地震疊前時(shí)間偏移的一種圖形處理器提速實(shí)現(xiàn)方法[J].地球物理學(xué)報(bào)，2009，52（1）：245-252.

Prestack reverse time migration based on GPU parallel accelerating algorithm/2012，36（4）：111-115

SHI Ying1，LU Jia-min2，KE Xuan1，TIAN Dong-sheng1，WANG Fei3
（1.School of Geosciences，Northeast Petroleum University，Daqing，Heilongjiang 163318，China；2.Exploration and Development Research Institute，Daqing Oilfield Co.Ltd.，Daqing，Heilongjiang 163712，China；3.School of Earth Science，China Petroleum University （Beijing），Beijing 102200，China）

In order to improve the complex subsurface imaging accuracy and computational efficiency of the algorithm，this paper presents an algorithm of prestack reverse time migration based on GPU（Graphic Processing Unit）accelerating which can image the underground complex structure effectively and accurately.By two-way wave equation to calculate wave field extrapolation，prestack reverse-time migration can overcome the dip limit，and the imaging algorithm is performed by high order finite difference in the paper.Wave field extrapolation and imaging condition are calculated by GPU parallel accelerating technology，comparing to conventional algorithm，its computation efficiency has been greatly improved，and it meets large amount of computation requirement in prestack reverse-time migration.The random boundary condition approach is adopted to obtain wavefield information，which reduces the memory demand but sacrifices the computation cost，and it solves the massy memory problem in reverse time migration.The tests on model illustrate that this approach can imaging complicated geological body efficiently and precisely.

reverse time migration；GPU；acceleration；high order finite difference；random boundary condition；complicated structure

TE132.1

A

2095-4107（2012）04-0111-05

DOI 10.3969/j.issn.2095-4107.2012.04.020

2012-05-18；編輯：任志平

國(guó)家“863”高技術(shù)研究發(fā)展計(jì)劃項(xiàng)目（2012AA061202）；國(guó)家自然科學(xué)基金青年基金項(xiàng)目（41104088，41004057）；黑龍江省教育廳科學(xué)技術(shù)研究項(xiàng)目（12511025）；黑龍江省博士后科學(xué)基金項(xiàng)目（LBH-Z11272）

石 穎（1976-），女，博士，副教授，主要從事地震資料處理方面的研究.