基于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-），女，博士，副教授，主要从事地震资料处理方面的研究.