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

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

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

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

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

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

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

0 引言

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

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

1 逆時偏移原理

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

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

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

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

U（x，z，t）對時間t的二階精度中心差分格式為

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

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

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

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

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

2.1 GPU加速技術(shù)

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

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

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

2.2 實現(xiàn)方法

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

圖1 GPU加速疊前逆時偏移流程

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

3 模型試算

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

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

圖2 Marmousi速度模型

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

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

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

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

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

4 結(jié)論

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

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

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

[1]袁剛，蔣波，曾華會.各向異性疊前時間偏移在塔里木碳酸鹽巖資料處理中的應(yīng)用[J].大慶石油學院學報，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]陳康，吳國忱.逆時偏移拉普拉斯算子濾波改進算法[J].石油地球物理勘探，2012，47（2）：249-255.

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

[8]李博，劉紅偉，劉國峰，等.地震疊前逆時偏移算法的CPU/GPU 實施對策[J].地球物理學報，2010，53（12）：2938-2943.

[9]張舒，褚艷利，趙開勇，等.GPU高性能運算之CUDA[M].北京：中國水利水電出版社，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]劉紅偉，劉洪，鄒振.地震疊前逆時偏移中的去噪與存儲[J].地球物理學報，2010，53（9）：2171-2180.

[13]李博，劉國峰，劉洪.地震疊前時間偏移的一種圖形處理器提速實現(xiàn)方法[J].地球物理學報，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；編輯：任志平

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

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