CN102565854B

CN102565854B - 一种海量数据gpu波动方程逆时偏移成像方法

Info

Publication number: CN102565854B
Application number: CN 201110443697
Authority: CN
Inventors: 李博
Original assignee: Institute of Geology and Geophysics of CAS
Current assignee: Institute of Geology and Geophysics of CAS
Priority date: 2011-12-27
Filing date: 2011-12-27
Publication date: 2013-07-31
Anticipated expiration: 2031-12-27
Also published as: CN102565854A

Abstract

本发明公开了一种海量数据GPU波动方程逆时偏移成像方法，逆时偏移的整体步骤包括三个模块，数据分配模块、逆时偏移模块、实时数据合并模块，所述数据分配模块将单炮数据分配给不同的GPU节点，所述逆时偏移模块将单炮数据计算成像并输出成像文件，所述实时数据合并模块对成像文件实时扫描，并统计分类，计算用户所需要的信息文件个数和大小，以百分比的形式通知用户需要多少时间完成数据合并。本发明的海量数据GPU波动方程逆时偏移成像方法，三个模块之间可以独立运行并行实施，更适合于海量数据的逆时偏移实施，计算效率高，计算速度快，容错能力强，并有效克服了海量成像时GPU内存不足的问题。

Description

一种海量数据GPU波动方程逆时偏移成像方法

技术领域

本发明涉及一种海量数据GPU波动方程逆时偏移成像方法，应用于石油勘探、天然气勘探、煤田勘探、矿产勘探等方面，属于地震勘探与开发领域。

背景技术

随着油气田勘探程度的不断深化，盐下油气藏勘探日益受到关注。采用传统的偏移方法无法圈定盐丘边界并对盐下构造进行成像，逆时偏移(RTM)采用双程波动方程，可以精确地描述波的传播过程，利用多次波、回转波等对陡倾角及反转构造进行准确成像，应用目前有限的偏移方法进行复杂地质构造成像。

逆时偏移技术作为复杂条件地震成像难题的解决方案，可以在遇到复杂险峻倾斜的反射面例如盐丘侧翼构造时，克服目前偏移方法的局限性，能够使大于70度的倾斜构造准确成像。逆时偏移技术采用双程声波波动方程，能够使用偏移校正包括垂直倾斜在内的所有地质倾斜，还可以实现倾斜构造的成像并保持振幅不变。逆时偏移提供的是一种可供选择的偏移方法：可在数据采集之后运行逆时偏移，选用适当的波场通过包括盐下、构造倾角大于70度等复杂的速度区域，以及在边界可能产生复杂的多次反射区域进行准确成像。

逆时偏移方法是通过双程波波动方程在时间域上对人工给予的震源子波正向传播和接受到的地震资料进行反向传播，结合成像条件实现偏移(Claerbout，1971)。由于逆时偏移需在同时计算震源正向传播和地震资料的反向传播过程中各个时刻的波场，因其各自时间的延拓方向不同，而致使实际应用中必须存储其中一个方向的传播过程，这需要甚大额外存储空间。并且如今的地震资料采集已经发展到了三维甚至四维，其数据量更是十分庞大，已经达到TB级甚至10TB级的叠前单炮数据。如此庞大的数据的IO时间将成为逆时偏移在推向海量数据的应用时不可回避的困难。例如，计算周期过长、容错能力弱、不易反复试验等问题。

另一方面，相较于单程波方程的波场延拓而言，逆时偏移运用的是双程波波动方程进行波场延拓，因其避免了上下行波的分离处理，因而成为最准确的成像算法，且不受倾角的限制，并能实现回转波和多次波成像。正基于此，地球物理学家从来没有放弃过对推进逆时偏移的追求。

GPU英文全称Graphic Processing Unit，中文翻译为“图形处理器”，是从CPU(中央处理器)衍生出来的的一个概念，它是显卡的心脏，和显卡板载内存成为一个子系统，共同决定PC系统的图形处理性能。图形处理器(GPU)用于通用计算(GPGPU)及其相关方面的问题目前已成为一个热门话题。GPU专为图像处理设计，存储系统实际上是一个二维的分段存储空间，包括一个区段号(从中读取图像)和二维地址(图像中的X、Y坐标)。此外，没有任何间接写指令。输出写地址由光栅处理器确定，而且不能由程序改变。这对于自然分布在存储器之中的算法而言是极大的挑战，存储器的访问效率是算法效率的关键问题所在。不同碎片的处理过程间不允许通信是一个SIMD数据并行执行单元，在所有碎片中独立执行代码。类似CUDA和OpenCL这样的高级编程语言的问世使编程新手也能够掌握GPU的计算性能优势。

利用GPU实现三维逆时偏移目前比较流行的方法是采用随机边界条件以计算换存储的方法来实现，该方法是利用了波场计算的可逆性的特点将边界反射随机化。但是这种方法会引入边界随机反射噪音在成像质量上难以上升。另外，若采用吸收边界条件则会由于存储的不连续性而大大降低GPU的实际计算效率。

发明内容

针对现有技术的缺陷，本发明公开了一种海量数据GPU波动方程逆时偏移成像方法，不仅克服了传统的逆时偏移方法在处理海量数据的应用时不可回避的困难，如计算周期过长、容错能力弱、不易反复试验等问题，并且也克服了近些年新兴起的利用GPU采用随机边界条件以计算换存储的方法来实现三维逆时偏移所导致成像质量差、计算效率低的缺陷。

为实现上述目的，本发明是通过如下技术方案实现的：

一种海量数据GPU波动方程逆时偏移成像方法，逆时偏移的整体步骤包括三个模块，数据分配模块、逆时偏移模块、实时数据合并模块，所述数据分配模块将单炮数据分配给不同的GPU节点，所述逆时偏移模块将单炮数据计算成像并输出成像文件，所述实时数据合并模块对成像文件实时扫描，并统计分类，计算用户所需要的信息文件个数和大小，以百分比的形式通知用户需要多少时间完成数据合并，通过此方法，计算密集的波场传播、边界条件和成像部分实现了借助GPU并行计算。

其中，所述数据分配模块利用轮换计算节点的方式把单炮数据分配给不同的GPU节点，每个节点领取一个单炮的计算任务。例如假设有12个计算节点，48个计算GPU，总共需要偏移的总炮数为480炮，则每个节点平均分配40炮。

进一步的，所述数据分配模块设有任务堆栈，所述任务堆栈存取分配好的任务信息，所述数据分配模块继续分配任务，所述逆时偏移模块从任务堆栈去领取偏移任务，这样就实现了不需要等待后续模块的完成就可以继续分配任务，可以实现数据任务分配模块独立化，有利于充分实现负载均衡；所述任务堆栈保存每个单炮数据被分配到的计算节点、节点内的GPU卡编号以及相同GPU卡的单炮顺序号。

在上述基础上，任务堆栈可以长期保存，只要输入数据和使用的计算资源不发生变化，任务的分配方式可以多次使用。

如本领域技术人员所知道的，逆时偏移部分包括成像区域的划分、炮点波场延拓、吸收边界的独立存储和计算、炮点和检波点的反向传播、相关成像等处理过程，在本发明中，所述逆时偏移模块的成像区域划分时采用先划分深度方向然后划分水平方向的方式。

优选的，为了适应地层速度的变化规律为层状变化，本发明将深度方向作为区域划分的首选方向。所述成像区域深度划分完成满足每一个成像层的的存储空间小于GPU的显存容量，然后在每个成像层中做水平方向的划分；所述成像区域水平方向按照方形网格划分，与GPU内核数目相适应，在每个方向上保留差分阶数所需要的边界宽度。

针对已有的GPU的运算单元特点，所述GPU内核以16个线程为一组，所述水平方向按照16x16的方形网格划分。

在本发明中，所用的划分方式是为了适应GPU内核线程组的特点而设计的：

用公式表示为：W＝16×N+2n，其中，W是成像子区域的宽度，N为任意大于0的整数，n是差分阶数。在满足GPU的运算特点的同时，上述的划分方法也满足空间高阶差分法的局部化的计算要求，例如三维8阶空间差分需要读取周围25个网格点的数据。由于GPU全局内存读取的延迟很高，在实际运行时可以先把需要的数据并行装载到共享存储器中，然后利用GPU对共享存储器的读取速度快的特点实现所述划分方法，大幅度提高计算效率。

如本领域技术人员所了解的，逆时偏移的波场延拓过程会涉及到网格边界的影响，已有的GPU逆时偏移方法成像网格包括边界一起计算，并且存储在一个连续内存空间内。在本发明中，所述逆时偏移模块的波场延拓过程涉及的网格边界的网格存储方式采用边界网格和成像网格分离，把二者分别存储在两块独立而且连续的内存空间中，对边界部分单独处理。

通过上述方法，解决了GPU处理非连续边界存储时的内存不连续所导致的计算效率降低的问题。尽管边界条件的中间结果需要借助通讯保存在CPU内存或硬盘上，但存储量已经大幅度的降低，易于实现；所用的方法即解决海量数据IO问题，又能保证成像质量，避免随机噪音的产生。

其中，在逆时偏移过程中以日期、时间、CPU节点号、GPU节点号、单炮的炮号、运行的时长、开始运行的时间来命名输出的单个成像文件，这样有利于数据合并时的实时统计功能的实现，并且适合于作业的断点恢复功能。

例如命名为tmp20110910-1020-CPU0-GPU0-shot123-2h-0910-0820。

在本发明中，实时数据合并模块采用指定目录的实时扫描方式，扫描文件名获取文件信息，通过对目录内的文件做统计分类，计算用户所需要的信息文件个数和大小，以百分比的形式通知用户需要多少时间完成数据合并。通过这种方法，带给用户很方便的使用体验。

例如，可以随时统计某一时间段内的偏移结果，可以按照计算节点统计偏移结果，可以利用炮号统计偏移结果等等。无需等待所有偏移完成后再做数据合并，这样可以方便用户在偏移过程中做质量监控。另外，所述的数据合并方案可以方便用户在偏移完成后随意的从偏移结果中抽取所需的成像道集或剖面，无需再次偏移，具体实现方式如下：

文件目录统计的目的是通过用户给定的通配文件名称，例如“tmp20110910*”的方式统计当前目录中的文件能够匹配该通配文件名的个数，将所有的单炮偏移结果做全部扫描，给出扫描结果。包括需要合并的单炮数据总数，在实际合并时用当前读取的文件数除以单炮总数作为数据合并的进度百分比显示给用户。这样用文件通配名称的方式和偏移输出文件名的命名方式可以实现用户不同需求的文件合并作业提交。例如统计2011年9月10日完成的偏移单炮数据的成像结果。

在本发明中，实时数据合并模块中用户可以自由定义输出成像空间的范围，例如主测线(inline)范围、联络测线(crossline)范围、偏移距(offset)范围。通过成像空间的并行分配给不同的节点可以实现数据合并的并行化处理。

通过本发明所述的方法，用于波动方程逆时偏移成像，不仅可以实现对海量数据的逆时偏移的计算处理，并且计算效率高，计算速度快，容错能力强，并有效克服了海量成像时GPU内存不足的问题，成像质量好。

附图说明

图1为本发明的用GPU实现逆时偏移流程图；

图2为本发明的方法数据分配模块利用任务堆栈实现负载均衡的炮数据分配示意图，其中不同颜色代表不同节点分配到的单炮；

图3为本发明的方法逆时偏移模块差分成像区域划分图，其中，GPU16个线程为一组；

图4为本发明的方法逆时偏移模块边界独立存储示意图；

图5为本发明的方法观测系统建立示意图；

图6为常规逆时偏移算法实现流程图；

图7为本专利的方法逆时偏移算法实现流程图。

具体实施方式

结合附图和具体实施可以更好的理解本发明的发明精神，本领域技术人员可以理解，所提供的具体实施仅为优选，而非仅限于采用此实现。本领域技术人员在理解本发明技术方案的基础上进行的变更，例如采用不同的GPU、CPU内核线程组数目改变等，依旧属于本发明的保护范围，本发明的保护涵盖于其权利要求及权利要求的等同变换。

在以下实施例中，采用的GPU型号为NVIDIA公司生产的Tesla C1060，并且下述方法适用于Tesla C1060之后发布的所有支持CUDA的GPU。本领域技术人员可以预期采用其他的通用计算方案平台在理解本发明原理的基础上也可以实现本发明所述的功能，例如MS的DirectCompute通用计算平台等，这依旧属于本发明的保护范围。

参考图1，公开了本发明的方法的流程框图，其中，“地震数据输入分配任务节点，并写入堆栈”、“任务数据堆栈”构成了数据分配模块的工作内容；“逆时偏移”、“单炮成像数据输出”构成了逆时偏移模块的工作内容；“合并数据抽取偏移距(offset)共像点道集”构成了实时数据合并模块的工作内容。

参考图2，显示了本发明的方法数据分配模块的任务堆栈实现对单炮计算任务分配的结果，实现了不同计算节点的炮数分配负载均衡。

参考图3，显示了本发明的方法逆时偏移模块成像区域划分的方法，划分时采用了先划分深度方向然后划分水平方向的方式，其中所用的GPU内核是以16个线程为一组；本领域技术人员可以理解对于其他不同类型的GPU，对划分方式相应调整即可。

参考图4，显示了本发明的方法逆时偏移模块对于吸收边界的网格存储方式，采用边界网格和成像网格分离的方案，把二者分别存储在两块独立而且连续的内存空间中。

参考图6和图7，可以充分理解本发明的方法相对于已有的方法的不同以及所进行的改进和改进的技术效果。

为了更好的说明本发明方法的实现方式，申请人在此提供一个计算示例：

运行环境：Linux操作系统

输入数据：地震数据的炮记录(按炮点分选，su格式)、地质速度模型(二进制格式)

输出数据：地震偏移成像剖面数据体(su格式)

其中，观测系统定义：利用原点坐标和方位角信息定义矩形处理工区网格，原点用fx＝0和fy＝0来指定，方位角用参数sita＝30来指定。其中fx，fy分别代表速度模型起始点的地面坐标，sita表示inline方向与X轴的夹角。如附图5所示。

数据规则化和任务分配：将地震数据的炮记录作为输入，并给出所用的节点个数和观测系统中的原点坐标和方位角参数，运行rtmspd。然后，用desdata.sh将准备好的数据分发到各个节点去。

任务的分配，首先手工建立一个标准的参数文件，该参数文件的实例如下：

建立好标准作业参数文件后，利用“creatjob标准参数文件名”生成批量作业的参数文件用以不同节点的偏移作业。

偏移成像：单独运行偏移作业可以使用命令“cudartm参数文件名”或者批量运行偏移程序可以使用命令“runrtm.sh参数文件名”。

数据合并输出：将偏移作业运行完毕后的整合叠加，即将所有的数据叠加生成偏移剖面。该指令的运行是：“megdata参数文件名”。

与已有技术(CPU)，对于同样的数据量，其计算效果如下：

Claims

1.一种海量数据GPU波动方程逆时偏移成像方法，其特征在于逆时偏移的整体步骤包括三个模块，数据分配模块、逆时偏移模块、实时数据合并模块，所述数据分配模块将单炮数据分配给不同的GPU节点，所述逆时偏移模块将单炮数据计算成像并输出成像文件，所述实时数据合并模块对成像文件实时扫描，并统计分类，计算用户所需要的信息文件个数和大小，以百分比的形式通知用户需要多少时间完成数据合并；所述逆时偏移模块的成像区域划分时采用先划分深度方向然后划分水平方向的方式，所述成像区域深度划分满足每一个成像层的的存储空间小于GPU的显存容量，然后在每个成像层中做水平方向的划分；所述成像区域水平方向按照方形网格划分，与GPU内核数目相适应，在每个方向上保留差分阶数所需要的边界宽度。

2.根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述数据分配模块利用轮换计算节点的方式把单炮数据分配给不同的GPU节点，每个节点领取一个单炮的计算任务。

3. 根据权利要求2所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述数据分配模块设有任务堆栈，所述任务堆栈存取分配好的任务信息，所述数据分配模块继续分配任务，所述任务堆栈保存每个单炮数据被分配到的计算节点、节点内的GPU卡编号以及相同GPU卡的单炮顺序号。

4. 根据权利要求3所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述逆时偏移模块从任务堆栈去领取偏移任务。

5. 根据权利要求4所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述水平方向按照16x16的方形网格划分，所述GPU内核以16个线程为一组。

6. 根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述逆时偏移模块的波场延拓过程涉及的网格边界的网格存储方式采用边界网格和成像网格分离，把二者分别存储在两块独立而且连续的内存空间中，对边界部分单独处理。

7．根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述逆时偏移模块逆时偏移过程中以日期、时间、CPU节点号、GPU节点号、单炮的炮号、运行的时长、开始运行的时间来命名输出的单个成像文件。

8. 根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述实时数据合并模块采用指定目录的实时扫描方式，通过文件名对目录内的文件做统计分类，计算用户所需要的信息文件个数和大小，以百分比的形式通知用户需要多少时间完成数据合并。