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

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

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

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

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

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

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

逆時偏移；GPU；加速；高階有限差分；隨機邊界條件；復雜構造

0 引言

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

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

1 逆時偏移原理

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

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

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

利用規則網格的高階有限差分求解聲波方程，其中U（x，z，t）對x，z的二階偏導數的2 N階精度中心差分格式分別近似為

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

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

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

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

完成震源波場的正向延拓和檢波點波場的反向延拓后，下步計算的關鍵是選擇合適的成像條件，互相關成像條件可為成像提供準確的動力學特性，且實現方法簡單.對于高陡界面，逆時偏移結果產生大量強振幅的低頻噪音，從而降低地震數據體的成像精度，需采用適當的方法壓制成像數據中產生的低頻噪音.

2 協同并行加速逆時偏移算法

2.1 GPU加速技術

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

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

GPU的編程平臺是CUDA（統一計算設備架構），它包括一個硬件驅動程序和一個應用程序接口（API）[13]，降低編程的難度，能夠通過函數調用訪問顯存和在GPU上執行指令，進行大規模的并行計算，這些改進使CUDA架構更加適合進行GPU通用計算.

2.2 實現方法

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

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

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

3 模型試算

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

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

圖2 Marmousi速度模型

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

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

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

圖4 Marmousi模型疊前逆時偏移結果

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

4 結論

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

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

（3）在提高算法計算精度的同時，根據GPU硬件的結構特點改進并行算法，優化存儲器訪問，使數據傳輸的同時，GPU也能夠正常進行運算，隱藏數據傳輸及訪問延遲，是未來的研究方向.

[1]袁剛，蔣波，曾華會.各向異性疊前時間偏移在塔里木碳酸鹽巖資料處理中的應用[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，Faria 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實現[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]李博，劉國峰，劉洪.地震疊前時間偏移的一種圖形處理器提速實現方法[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”高技術研究發展計劃項目（2012AA061202）；國家自然科學基金青年基金項目（41104088，41004057）；黑龍江省教育廳科學技術研究項目（12511025）；黑龍江省博士后科學基金項目（LBH-Z11272）

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