基于GPU并行加速的VSP数据逆时偏移

郭雪豹，王建民，王维红，王云专，柯 璇，刘诗竹

（1.东北石油大学 地球科学学院，黑龙江 大庆 163318； 2.大庆油田有限责任公司 勘探开发研究院，黑龙江 大庆163712）

0 引言

随着我国各大油田勘探的不断深入，探区内有利储藏大多已勘探完毕，而剩余的油气藏多具有油水分布复杂、薄储层变化快等特点.对于日趋复杂多样化的油气藏类型，常规的偏移方法已无法对地下深层次的复杂构造准确成像.基于波动理论的逆时偏移与同样基于波动理论的克希霍夫偏移和单程波偏移相比，波动方程不存在任何假设，既不会受到速度横向剧烈变化的影响，也不存在成像倾角的限制，理论上能够对多次波、回转波等准确成像，是目前成像精度最高的方法.

逆时偏移方法最早是在1983年被提出的，Baysal E等［1］提出不同的逆时偏移概念；Whitmore N D、Mc Mechan G A、Loeuenthal D等［2-4］将它应用于叠后偏移.随着油气资源的复杂多样化，常规的地面地震逆时偏移方法已经无法满足精细化勘探的需求，尤其对地下深层一些微幅构造及井旁的小断层等信息.VSP技术是与地面地震观测方法相对应的，它将检波器置于井中垂向布置.因此，比地面检波器更加接近目的层，减少地下浅层对地震信号的干扰，直接接收来自目的层的反射，同样还能够接收到一些地面检波器所接收不到的陡倾角信息，具有比地面地震记录更高的分辨率.郭建［5］分析VSP技术的应用现状.朱金明等［6］通过按时间逆的次序求解双程无反射波动方程，实现VSP逆时偏移.孙文博等［7］采用伪谱法实现逆时偏移方法，分别对分离波场偏移处理，得到较好的效果.Alok K S等［8］采用包含多次波在内的全波场成像，提高照明范围.Sun Wenbo等［9］在角度域实现VSP逆时偏移.Mark E W等［10］利用walkaway VSP资料在盐丘侧翼上取得较好的成像效果.VSP资料在对地下深层构造成像方面具有优势，并且其内部丰富的波场信息也有助于更高精度成像.

采用VSP数据进行逆时偏移，利用VSP资料的高分辨率、高信噪比的特点，实现地下微幅构造的准确成像.根据设计的VSP观测系统，笔者正演VSP地震记录，利用VSP资料的丰富波场信息进行逆时偏移，以弥补常规地面地震在深层构造及井旁小构造方面的成像不足；逆时偏移算法中波场计算部分由GPU执行，其余逻辑判断部分由CPU执行，充分挖掘两者的计算特性来获得最高加速比.算法经GPU加速后，计算效率得到了显著提升.

1 VSP观测系统

VSP观测系统是一种由地表震源激发，井中检波器接收的观测方式（见图1）.首先，与传统地面检波器接收到的反射信号路径相比，VSP检波器接收的反射波路径明显更短，避免了反射波回到地表时所要经过的浅层部分，减少地表低速带对反射波信息中高频成分的吸收，振幅畸变相对更小.其次，由于检波器布置在地下，与地表的检波器相比更加贴近目的层，因此能够直接接收到来自目的层的反射信息，有助于研究来自目的层的单一反射信息；在地下存在陡倾角构造时，地表检波器范围有限，可能损失部分反射信息，而VSP的检波器则能接收到.根据波场信息，VSP波场较地面地震波场更为丰富，不仅存在上行波（图中以上行一次反射波为例），还有下行波（图中以下行直达波为例），将有利于逆时偏移方法对VSP波场信息的充分利用.因此，与地面地震记录相比，VSP数据的分辨率和信噪比更高，将给地下深层成像中带来明显优势.

2 逆时偏移原理

逆时偏移过程可以分为震源波场正传和检波点处波场反传.首先进行震源波场正传，在正传过程中为削弱用计算机有限存储模拟地下无限介质的波场传播所带来的边界反射，需要添加边界条件，削弱边界反射对波场的影响，并且保存每一时刻的波场信息；然后进行检波点处波场反传，与对应时刻的震源波场应用成像条件得到逆时偏移结果，对逆时偏移结果去噪处理得到最后的成像剖面.

二维声波方程：

式中：p为波场值；v为介质速度.

通过有限差分法对波场进行数值模拟，时间二阶、空间2L阶精度的有限差分格式［12］为

式中：Δx、Δz分别为x、z方向的网格大小；Δt为时间步长；i、j分别为x、z方向网格点的位置；n为时间点；al为有限差分系数.

3 GPU加速技术

尽管逆时偏移方法成像精度极高，但需解决计算量大和存储需求大的问题.GPU即图形处理单元，经NVIDIA公司发布CUDA后，使它可以被开发人员编写指令程序，便于让GPU执行CPU所分布的计算任务.GPU具有比CPU更多的计算核心，更易处理大量的并行计算任务，对于逆时偏移过程中的大量波场值计算更加具有优势［11］，因此可以显著的提高算法的计算效率.

在逆时偏移过程中，计算量主要体现在模拟波场传播计算中.对于一个二维的剖面，将它按纵横向进行离散化，得到用网格点表示的剖面.在偏移过程中，不断利用式（2）计算每一个网格点上不同时刻的波场值.对于同一时刻的波场剖面，每个网格点上的波场值计算并无先后顺序，即同一时刻的所有网格点的计算是并行的.首先在GPU上开辟同一时刻所需要的线程数，然后将每一个网格点上的波场值计算放到对应的线程中，GPU上的计算核心多（计算核心即流多处理器，其本身包含标量流处理器，每个线程都是被发射到一个标量流处理器上执行），因此比CPU更加适合并行计算［13-15］.

由式（2）可知，在计算一点的波场值时需要多个点的波场值参与运算，因此需要在GPU上开辟部分存储以存放波场值.在GPU的存储器中，全局存储器是存储空间最大的，但访问速度远不如其内部的共享存储器.由于共享存储器存储有限，首先将计算的数据先读入全局存储器；然后在每个线程运算时，将所需要的数据再从全局存储器复制到共享存储器中，在计算时可以直接访问共享存储器，有效加快计算式的访问速度，从而进一步提高计算效率.GPU加速计算流程（见图2）.

4 模型测试

4.1 盐丘模型

采用盐丘模型测试算法（见图3），网格大小及时间步长根据稳定性条件选取.模型横、纵向网格点数分别为800、600，横、纵向网格大小分别为5、5m.采用地面放炮、井中接收的观测方式，震源从地表左端50m处开始，向右每隔50m放一炮，共80炮；在分别距模型左端0、4 000m处布置2口井，井中每隔5m布置一个检波器，每口井放置600个；震源处雷克子波频率为40Hz，时间步长为0.5ms，时间采样点为7 500.

在数值模拟过程中，采用保存部分波场信息的完全匹配层边界条件［16-19］，吸收由计算机有限存储所带来的边界反射衰减，成像条件为互相关成像条件［20］，对逆时偏移结果采用拉普拉斯去噪［21-23］处理.盐丘模型运算采用GPU加速节点运行，硬件配置：CPU：E5-2630（2.3GHz／15M／6核心）＊2；内存：64 GB；系统硬盘：1TB；GPU：Nvidia Tesla K10＊3（3 072核心／4.57T单精度／0.19T双精度／8G显存）.

VSP逆时偏移的成像剖面见图4.由图3和图4可知，地下构造基本成像清晰，浅层部分成像略显不足，井旁地层成像准确且清晰，盐丘轮廓清楚，证实VSP逆时偏移对井旁高分辨率成像的优势，能够有效利用来自地下深层的反射信息.

4.2 Marmous模型

利用Marmous模型进行测试，网格大小及时间采样由模型发布方参数决定.模型纵、横向网格点数分别为750、993，纵、横向网格大小分别为6.25、4.00m（见图5）；采用地面放炮、井中接收的观测方式，震源从地表左端62.50m处开始，向右每隔62.50m放一炮，共99炮；在分别位于模型的62.50、3 125.00、6 187.50m处3口井，井中每隔4.00m布置一个检波器，每口井放置750个；震源处雷克子波频率为40Hz，时间步长为0.4ms，时间采样点为10 000.

VSP逆时偏移的成像剖面见图6，地面地震的偏移结果见图7.由图5可知，Marmous模型地下地层较多，地势起伏大，给成像带来一定难度.由图7可知，地面地震逆时偏移可以较清晰地对它进行成像，尤其在浅层部分，地层清晰连续，但在2 000m以下的部位开始模糊.在图6中2 000m以下的部位，构造形态基本正确，成像清晰，突出VSP逆时偏移对于深层构造的成像优势.在计算过程中，采用CPU偏移一炮需要2 844s，在相同参数下采用GPU偏移一炮仅需36 s，有效提高偏移的速度70倍以上.

5 结论

（1）逆时偏移方法能够充分利用VSP数据的丰富波场信息及高分辨率、高信噪比的优势，在井旁及深层构造成像方面具有比地面地震逆时偏移更好的成像效果.但由于对浅层反射信号接收不足，在浅层部分成像效果不如地面地震逆时偏移成像清晰.

（2）逆时偏移算法对计算机的计算能力要求较高，充分利用GPU上的共享存储器，使其达到显著的加速比，增强算法的实用性.

（3）二维VSP逆时偏移在一定程度上弥补地面地震逆时偏移对深层成像不足的问题，但随着勘探目标的日趋复杂，二维VSP逆时偏移也难以达到复杂构造精细勘探的需求，开展三维VSP逆时偏移技术是未来的研究方向.

［1］Baysal E，Kosloff D D，Sherwood W C.Reverse time migration［J］.Geophysics，1983，48（11）：1514-1524.

［2］Whitmore N D.Iterative depth migration by backward time propagation［C］∥SEG Technical Program Expanded Abstracts，1983：382-385.

［3］McMechan G A.Migration by extrapolation of time-dependent boundary values［J］.Geophysical Prospecting，1983，31（3）：413-420.

［4］Loewenthal D，Mulfti I R.Reversed time migration in spatial frequency domain［J］.Geophysics，1983，48（5）：627-635.

［5］郭建.VSP技术应用现状及发展趋势［J］.勘探地球物理进展，2004，27（1）：1-8.Guo Jian.The application status and development trends of VSP technology［J］.Progress in Exploration Geophysics，2004，27（1）：1-8.

［6］朱金明，颜俊华.VSP逆时偏移［J］.石油地球物理勘探，1991，26（5）：564-570.Zhu Jinming，Yan Junhua.VSP reverse-time migration［J］.Oil Geophysical Prospecting，1991，26（5）：564-570.

［7］孙文博，孙赞东.基于伪谱法的 VSP逆时偏移及其应用研究［J］.地球物理学报，2010，53（9）：2196-2203.Sun Wenbo，Sun Zandong.VSP reverse time migration based on the pseudo-spectral method and its applications［J］.Chinese J.Geophys，2010，53（9）：2196-2203.

［8］Alok K S，Xander S，Eric V.VSP imaging using all multiples：Full wavefield migration approach［C］∥SEG Technical Program Expanded Abstracts，2012：1-6.

［9］Sun Wenbo，Sun Zandong.Angle-domain VSP reverse-time migration and its application［C］∥SEG Technical Program Expanded Abstracts，2010：4339-4343.

［10］Mark E W，Lu Rongrong，Xander C，et al.A novel application of time-reversed acoustic：Salt-dome flank imaging using walkaway VSP survey［J］.Geophysics，2006，71（2）：A7-A11.

［11］刘红伟，李博，刘洪.地震叠前逆时偏移高阶有限差分算法及 GPU实现［J］.地球物理学报，2010，53（7）：1725-1733.Liu Hongwei，Li Bo，Liu Hong.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation［J］.Chinese J.Geophys，2010，53（7）：1725-1733.

［12］刘洋，李承楚，牟永光.任意偶数阶精度的有限差分数值模拟方法［J］.石油地球物理勘探，1998，33（1）：1-10.Liu Yang，Li Chengchu，Mu Yongguo.Finite-difference numerical modeling of any even-order accuracy［J］.Oil Geophysical Prospecting，1998，33（1）：1-10.

［13］李博，刘红伟，刘国峰，等.地震叠前逆时偏移算法的CPU／GPU实施对策［J］.地球物理学报，2010，53（12）：2938-2943.Li Bo，Liu Hongwei，Liu Guofeng，et al.Computational strategy of seismic pre-stack reverse time migration on CPU／GPU［J］.Chinese J.Geophys，2010，53（12）：2938-2943.

［14］石颖，陆加敏，柯璇，等.基于 GPU 并行加速的叠前逆时偏移方法研究［J］.东北石油大学学报，2012，36（4）：111-115.Shi Ying，Lu Jiamin，Ke Xuan，et al.Prestack reverse time migration based on GPU parallel accelerating algorithm［J］.Journal of Northeast Petroleum University，2012，36（4）：111-115.

［15］柯璇，石颖，刘诗竹.地震资料逆时偏移中的图像处理器加速算法［J］.计算机系统应用，2013，22（11）：115-118.Ke Xuan，Shi Ying，Liu Shizhu.GPU acceleration algorithm of the reverse time migration with seismic data［J］.Computer Systems＆ Applications，2013，22（11）：115-118.

［16］郭念民，吴国忱.基于PML边界的变网格高阶有限差分声波方程逆时偏移［J］.石油地球物理勘探，2012，47（2）：256-265.Guo Nianmin，Wu Guochen.Highorder finite difference method in reverse-time migration with variable grids based on PML boundary condition［J］.Oil Geophysical Prospecting，2012，47（2）：256-265.

［17］王守东.声波方程完全匹配层吸收边界［J］.石油地球物理勘探，2003，38（1）：31-34.Wang Shoudong.Absorbing boundary condition for acoustic wave equation by perfectly matched layer［J］.Oil Geophysical Prospecting，2003，38（1）：31-34.

［18］王永刚，邢文军，谢万学，等.完全匹配层吸收边界条件的研究［J］.中国石油大学学报：自然科学版，2007，31（1）：19-24.Wang Yonggang，Xing Wenjun，Xie Wanxue，et al.Study of absorbing boundary condition by perfectly matched layer［J］.Journal of China University of Petroleum：Natural Science Edition，2007，31（1）：19-24.

［19］胡昊，刘伊克，常旭，等.逆时偏移计算中的边界处理分析及应用［J］.地球物理学报，2013，56（6）：2033-2042.Hu Hao，Liu Yike，Chang Xu，et al.Analysis and application of boundary treatment for the computation of reverse-time migration［J］.Chinese J.Geophys，2013，56（6）：2033-2042.

［20］薛东川，张云鹏，朱振宇，等.几种叠前逆时偏移成像条件的比较［C］.中国地球物理，2012：581.Xue Dongchuan，Zhang Yunpeng，Zhu Zhenyu，et al.Comparison of several imaging condition for prestack reverse-time migration［C］.The Chinese Geophysical，2012：581.

［21］刘红伟，刘洪，邹振，等.地震叠前逆时偏移中的去噪与存储［J］.地球物理学报，2010，53（9）：2171-2180.Liu Hongwei，Liu Hong，Zou Zhen，et al.The problems of denoise storage in seismic reverse time migration［J］.Chinese J.Geophys.，2010，53（9）：2171-2180.

［22］Zhang Y，Sun J.Practical issues in reverse time migration：true amplitude gathers，noise removal and harmonic source encoding［J］.First Break，2009，26（1）：29-35.

［23］陈可洋.基于拉普拉斯算子的叠前逆时噪声压制方法［J］.岩性油气藏，2011，23（5）：87-95.Chen Keyang.Pre-stack reverse-time noise suppressing method based on Laplacian operator［J］.Lithologic Reservoirs，2011，23（5）：87-95.