一种地震叠前时间偏移走时处理方法及装置

摘要

本发明公开了一种地震叠前时间偏移走时处理方法及装置，涉及计算机高性能计算、并行计算领域、石油地震勘探领域。本发明公开的地震叠前时间偏移走时处理方法包括：对成像空间建立三维X-Y-Z坐标系，每个点按照(x，y，z)坐标唯一确定，对X、Y坐标相同的点所构成与Z轴平行的竖线之间并行进行走时计算，对每条竖线内的点并行进行走时计算；建立计算统一设备架构(CUDA)线程模型和CUDA内存模型以进行上述并行走时计算。本发明还公开了一种地震叠前时间偏移走时处理装置。本申请技术方案解决了叠前时间偏移走时计算的瓶颈，令叠前时间偏移处理性能获得大大提升。

说明书

技术领域

本发明涉及计算机高性能计算、并行计算领域、石油地震勘探领域，具 体涉及一种针对NVIDIA的Fermi GPU设计的地震叠前时间偏移走时计算并 行算法。

背景技术

地震资料处理是石油和天然气勘探开发领域中非常重要的环节。传统的 方法是通过人工放炮方式产生地震波，通过地面检波器将地下不同地质层反 射回来的地震波信号收集后，利用大型计算机通过多套专业处理软件和一套 完整的叠前时间偏移、叠前深度偏移软件系统进行资料处理，从而得到地下 的构造以及成像。为石油钻井提供更加可靠的勘探数据，用于勘探专家进行 下一步的分析和解释，掌握地下的油气构造。

叠前时间偏移是构造复杂成像最有效的方法之一，它能适应纵向速度变 化较大的情况，适用于大倾角的偏移成像。而其中走时计算是计算地震波从 震源传到观测点所经过的时间，即计算地震波从震源传到成像点所经过的时 间与地震波从成像点传到观测点所走过的时间之和。因此，走时计算是叠前 时间偏移中最关键、最费时的一个环节，其采用的算法大致分为三类，即直 射线、弯曲射线以及非对称走时计算。由于叠前时间偏移动辄需要处理数以 TB的海量数据，以串行方式实现的走时计算算法，处理性能十分低，已经 严重制约工业生产。

发明内容

本发明所要解决的技术问题是，提供一种地震叠前时间偏移走时处理方 法及装置，以提高叠前时间偏移处理性能。

为了解决上述技术问题，本发明公开了一种地震叠前时间偏移走时处理 方法，包括：

对成像空间建立三维X-Y-Z坐标系，每个点按照(x，y，z)坐标唯一确 定，对X、Y坐标相同的点所构成与Z轴平行的竖线之间并行进行走时计算， 对每条竖线内的点并行进行走时计算；

建立计算统一设备架构(CUDA)线程模型和CUDA内存模型以进行上 述并行走时计算。

较佳地，上述方法中，建立CUDA线程模型指：

将整个成像空间XY平面划分为NX*NY/4块，每一个Block的线程计 算4条垂直于地面的竖线的点的走时计算，每一个Block包括128线程，其 中每32个线程计算一条竖线。

较佳地，上述方法中，所建立的CUDA内存模型包括全局内存、纹理内 存、共享内存和常量内存。

较佳地，上述方法中，对每条竖线内的点并行进行走时计算指：

对于每条竖线内的点，先并行计算出索引值能被8整除的点的走时，对 不能被8整除的点则依赖于与其临近的已精确计算的两个点的值作线性插值 得到该点的走时。

本发明还公开了一种地震叠前时间偏移走时处理装置，包括：

第一单元，对成像空间建立三维X-Y-Z坐标系，每个点按照(x，y，z) 坐标唯一确定；

第二单元，对三维X-Y-Z坐标系中X、Y坐标相同的点所构成与Z轴平 行的竖线之间并行进行走时计算，对每条竖线内的点并行进行走时计算，以 及建立计算统一设备架构(CUDA)线程模型和CUDA内存模型以进行上述 并行走时计算。

较佳地，上述装置中，第二单元，建立CUDA线程模型指，将整个成像 空间XY平面划分为NX*NY/4块，每一个Block的线程计算4条垂直于地 面的竖线的点的走时计算，每一个Block包括128线程，其中每32个线程计 算一条竖线。

较佳地，上述装置中，第二单元，所建立的CUDA内存模型包括全局内 存、纹理内存、共享内存和常量内存。

较佳地，上述装置中，所述第二单元，对每条竖线内的点并行进行走时 计算指：对于每条竖线内的点，先并行计算出索引值能被8整除的点的走时， 对不能被8整除的点则依赖于与其临近的已精确计算的两个点的值作线性插 值得到该点的走时。

本申请技术方案充分利用了GPU的多核处理能力，解决了叠前时间偏移 走时计算的瓶颈，令叠前时间偏移处理性能获得大大提升。即本申请技术方 案满足了石油地震勘探资料处理的需求，并且降低了机房构建成本和管理、 运行、维护费用。

具体实施方式

图1为本实施例中地震叠前时间偏移走时处理流程图；

图2为本实施例中8层内32个点并行走时计算线程图；

图3为本实施例中1层内32个点并行走时计算线程图；

图4为本实施例中叠前时间偏移串行走时计算成像效果图；

图5为本实施例中叠前时间偏移并行走时计算成像效果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下文将结合附图对 本发明技术方案作进一步详细说明。需要说明的是，在不冲突的情况下，本 申请的实施例和实施例中的特征可以任意相互组合。

实施例1

本实施例，提供一种地震叠前时间偏移走时处理方法，该方法以 Kirchhoff积分法叠前时间偏移走时计算串行算法为基础，并基于NVIDIA的 Fermi GPU硬件体系结构，利用CUDA(Computing Unified Device  Architecture，计算统一设备架构)技术对原有串行走时计算算法进行并行化 设计的。具体地，该方法的过程如图1所示，包括如下步骤：

步骤100，确定并行走时计算算法；

该步骤中，对成像空间建立三维X-Y-Z坐标系，假设第一维大小为NX， 第二维大小为NY，第三维大小为NZ，那么需要进行走时计算的总点数为 NX*NY*NZ，每个点按照(x，y，z)坐标唯一确定。X、Y坐标相同的点 构成一条与Z轴平行的线，物理上是一条垂直于地面的竖线，那么这样的竖 线有NX*NY条，并行走时计算采用两级并行，竖线与竖线之间并行进行走 时计算；每条竖线内的点并行进行走时计算。

其中，对于同一条竖线上每一个点的走时T可采用多种方式计算。而本 实施例中对同一条竖线上每一个点的走时T采用估算方式计算，即首先并 行计算出索引值能被8整除的点的走时T，其余点的走时T依赖于与它临近 的已精确计算的两个点的值作线性插值得。

步骤200，确定CUDA线程模型；

该步骤中，将整个成像空间XY平面看成一个网格(Grid)，将Grid定 义为(NX，NY/4)，即把整个Grid分为NX*NY/4块(Block)，每一个 Block的线程计算4条垂直于地面的竖线的点的走时计算；将Block定义为 (32，4)，即每一个Block包括128线程(Thread)，其中每32个线程计算 一条竖线。

步骤300，确定CUDA内存模型，主要是根据数据访问特点及NVIDIA  Fermi GPU内存资源特性，选择不同的内存存放不同的数据，以达到性能最 优。

该步骤中，CUDA内存模型至少包括Global memory(全局内存)、Texture  memory(纹理内存)、Shared memory(共享内存)和Constant memory(常 量内存)。

其中，Global memory指：成像空间数据存储方式是先按照Z方向连续， 同时Z方向维度大小总是32的整数倍。Z方向的计算总是以32为单位进行 计算。因此在同一时刻，32个线程可以同时访问成像空间的32个点，而且 这32个点在内存中是连续的，并且可以做到线程与访问点的一一对应，从而 形成对Global Memory的合并访问，提高访存性能；

Texture memory指：由于GPU中的Texture memory有cache，把只读并 且频繁访问的大输入道数据存放其中，将提高访存性能；

Shared memory指：由于Shared memory为GPU的片上内存，访问速度 快，对于一个Block块中公共的数据，如能被8整除的点的走时、速度场等 小数据，可以放入共享内存中，将提高访存性能；

Constant memory指：对于PSTM计算中只读的、频繁被访问的小数据块 可以放入Constant memory中，将大大提高访存性能。

下面结合具体应用场景，详细介绍上述方法的实现过程。

首先进行走时计算部分串行算法的并行性分析，将整个走时计算的代码 从叠前时间偏移程序中抽取出来，研究走时计算的串行算法，分析其是否具 有并行性，研究发现每一点的走时计算是并行的，数据不存在依赖性，可以 设计走时计算并行算法，利用CUDA实现，移植到NVIDIA GPU中运行。

其次，需要确定并行走时计算算法，选定一个体偏成像空间，假设其第 一维大小为177，第二维大小为1097，第三维大小为1504，那么需要进行走 时计算的总点数为177*1097*1504，即垂直于地面的竖线有177*1097＝ 194169条，并行走时计算采用两级并行，上述194169竖线之间互不依赖， 可以并行进行走时计算；每条竖线内的1504个点可以并行进行走时计算。并 且对于同一条竖线上的1504个点的走时T并不需要精确计算出，而是首先 并行计算出索引值能被8整除的点的走时T，其余点的走时T依赖于与它临 近的已精确计算的两个点的值作线性插值得。其具体过程如下：

A、如图2所示，以8层为一个计算单元，1个warp内32个线程精确并 行计算出索引值能被8整除的33个点的走时T，其中线程0多计算一个点；

B、判断一层内所有点的索引号落在哪两个最近的索引号能被8整除的 点之间，然后计算其插值因子；

C、并行进行插值计算，如图3所示，以一层为一个计算单元，1个warp 内32个线程并行插值计算一层内所有点(包括索引号能被8整除的点和索引 号不能被8整除的点)的走时T，循环8次后，整个8层的所有点的走时T 将被计算完成。

接下来，确定CUDA线程模型，将整个上述成像空间XY平面看成一个 网格(Grid)，其Grid定义为dim3 Grid(177，(1097+3)/4)，即把整个Grid 分为177*(1097+3)/4块(Block)，每一个Block的线程计算4条垂直于 地面的竖线的点的走时计算；把Block定义为dim3 Block(32，4)，即每一个 Block包括128线程(Thread)，其中每32个线程计算一条竖线。

最后确定CUDA内存，即根据数据访问特点及NVIDIA Fermi GPU内存 资源特性，选择不同的内存存放不同的数据，以达到性能最优。其中，CUDA 内存至少包括Global memory、Texture memory、Shared memory和Constant  memory。

下面针对上述具体计算，进行正确性测试，主要通过运行上述体偏作业， 获得图4、图5成像效果图。其中，图4为叠前时间偏移走时计算采用串行 算法，运行于CPU上所获得的成像效果。图5为叠前时间偏移走时计算采用 并行算法，运行于GPU上所获得的成像效果。从结果来看，两幅图像不存在 明显差异，证明基于GPU的地震叠前时间偏移走时处理方法是正确的。

再来了解性能测试中的测试环境及测试数据。测试环境包括硬件环境、 软件环境、叠前时间偏移运行内核，其中叠前时间偏移运行内核为叠前时间 偏移的核心计算部分，即包括走时计算部分；测试数据为输入的测试地震道 数据集。对于成像空间而言，其中第一维为X轴方向的大小；第二维为Y轴 方向的大小；第三维为Z轴方向的大小，具体各项参数如表1所示。

表1为测试环境和测试数据参数表

最后看测试结果，为了保证测试性能结果的稳定性，对上述体偏作业进 行了10次测试，运行于CPU上的串行走时计算10次的平均时间是54320 秒，而对于上述同样体偏作业而言，运行GPU上并行走时计算10次的平均 时间是1646秒，并行走时计算算法的性能是串行走时计算算法的54320/1646 ＝33倍。

由上述实施例可以看出，基于GPU的地震叠前时间偏移走时处理方案是 正确的，成像效果跟CPU串行走时计算算法一样，并且其性能是串行走时计 算算法的33倍，这样不但满足了石油地震勘探处理的需求，而且大大降低了 功耗，减少了机房构建成本和管理、运行、维护费用，而且这种方法实现简 单，需要的开发成本低。

实施例2

本实施例介绍一种地震叠前时间偏移走时处理装置，可实现上述实施例 1的方法。该装置包括如下各单元。

第一单元，对成像空间建立三维X-Y-Z坐标系，每个点按照(x，y，z) 坐标唯一确定；

第二单元，对三维X-Y-Z坐标系中X、Y坐标相同的点所构成与Z轴平 行的竖线之间并行进行走时计算，对每条竖线内的点并行进行走时计算，以 及建立计算统一设备架构(CUDA)的线程模型和CUDA内存模型以进行上 述并行走时计算。需要说明的是第二单元所建立的CUDA内存模型包括全局 内存、纹理内存、共享内存和常量内存。

上述第二单元，建立CUDA线程模型指，将整个成像空间XY平面划分 为NX*NY/4块，每一个Block的线程计算4条垂直于地面的竖线的点的 走时计算，每一个Block包括128线程，其中每32个线程计算一条竖线。

而上述第二单元，对每条竖线内的点并行进行走时计算指：对于每条竖 线内的点，先并行计算出索引值能被8整除的点的走时，对不能被8整除的 点则依赖于与其临近的已精确计算的两个点的值作线性插值得到该点的走 时。

其他细节可参见实施例1，在此不再赘述。

本领域普通技术人员可以理解上述方法中的全部或部分步骤可通过程序 来指令相关硬件完成，所述程序可以存储于计算机可读存储介质中，如只读 存储器、磁盘或光盘等。可选地，上述实施例的全部或部分步骤也可以使用 一个或多个集成电路来实现。相应地，上述实施例中的各模块/单元可以采用 硬件的形式实现，也可以采用软件功能模块的形式实现。本申请不限制于任 何特定形式的硬件和软件的结合。

以上所述，仅为本发明的较佳实例而已，并非用于限定本发明的保护范 围。凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等， 均应包含在本发明的保护范围之内。