高新成 (东北石油大学现代教育技术中心，黑龙江 大庆 163318)石颖 (东北石油大学地球科学学院，黑龙江 大庆 163318)



基于CPU/GPU异构平台的叠前逆时偏移成像系统

高新成 (东北石油大学现代教育技术中心，黑龙江 大庆 163318)石颖 (东北石油大学地球科学学院，黑龙江 大庆 163318)

针对目前地震资料偏移成像过程中代码复杂分散、可移植性差和可视化操作弱等问题，结合叠前逆时偏移业务需求，研究了CPU/GPU协同并行计算技术优化叠前逆时偏移成像算法，设计了逆时偏移系统流程与功能，开发了一套基于CPU/GPU异构平台的叠前逆时偏移成像系统。利用自设模型和传统模型数据对系统进行测试，测试结果表明系统操作方便可行，运行效果良好，能够有效的对复杂构造地质体精确成像。

逆时偏移；CPU/GPU并行加速；PML边界条件；QT

随着地震资料偏移成像应用推广，在逆时偏移成像处理中积累了很多的应用程序和方法，如何将这些程序和方法集成整合起来，实现偏移成像的可视化操作和高效管理成为当前研究的重要课题[1]。目前叠前逆时偏移是获得复杂构造精确成像最有效的方法[2]。为此，笔者结合叠前逆时偏移业务需求研究地震波场正演模拟、成像条件、边界吸收条件和噪音压制算法等，形成了一套叠前逆时偏移成像方法；并研究了CPU/GPU异构平台并行加速技术，将其应用于叠前逆时偏移成像方法中地震波场正演计算、成像条件和波场反传计算等计算环节，开发了一套基于CPU/GPU异构平台的叠前逆时偏移成像系统，为油田开发资料的解释提供技术支持。

1 CPU/GPU异构平台的逆时偏移流程

1.1 叠前逆时偏移成像算法

笔者采用逆时偏移算法，在波场正演过程中采用随机边界条件和PML完全吸收边界条件相结合的2种边界处理方式，成像过程中采用互相关成像条件，成像后采用相应滤波法进行去噪[3～7]。逆时偏移算法实现流程如图1所示。算法实现步骤如下：

步1 使用波动方程有限差分法沿时间正向延拓波场计算到Tmax时刻，保存最后2个时刻波场值。

步2 采用PML边界条件，使用波动方程有限差分法沿时间正向延拓波场计算到Tmax时刻，保存每个时刻的检波点数据。

步3 读取最后2个时刻波场值，使用波动方程有限差分法沿时间反向延拓计算到0时刻，计算每个时刻波场值。

步4 读取检波点数据，使用波动方程有限差分法沿时间反向延拓计算到0时刻，计算每个时刻波场值。

步5 读取步3和步4同一时刻2个波场值，利用互相关成像条件进行成像，从Tmax到0时刻。

步6 循环计算多炮数据到Nshot，完成多炮成像结果叠加。

步7 采用滤波法对成像结果进行压制低频噪声。

步8 输出成像结果。

1.2 CPU/GPU协同加速逆时偏移流程

将CUDA编程思想引入到逆时偏移波场延拓和相关成像算法中，采用基于CUDA架构的CPU/GPU协同工作模式来处理逆时偏移计算任务。其主要思想是GPU端负责完成波场延拓和相关成像等并行计算任务，CPU端负责数据输入和输出等相关操作控制工作。

图2 CPU/GPU协同加速逆时偏移流程

首先，将激发点初始波场值和检波点最大时刻波场值由内部存储器(内存)传至设备存储器(显存)中[8]。其次，设置GPU计算维度，将GPU多核处理器划分为相应个数的计算块Block；同时将每个Block划分为若干个计算线程Thread；开始执行炮点波场正传和检波点波场反传计算任务。根据波场计算相关成像条件，需要同一时刻正向传播的震源波场和反向传播的检波点波场，因此互相关成像计算也是由GPU来完成[9，10]，具体实现流程如图2所示。

2 系统流程与功能设计

2.1 系统设计流程

根据逆时偏移处理流程，将逆时偏移系统数据处理分为工区创建、数据加载、参数设置、作业发送和任务执行等多个步骤，并结合CPU/GPU协同并行计算模式[11]设计了逆时偏移数据处理系统流程(见图3)。

图3 逆时偏移系统处理流程

2.2 系统功能设计

根据逆时偏移的实际业务需求，设计并开发了叠前逆时偏移成像系统，主要完成二维和三维叠前逆时深度偏移成像功能。系统主要包括资源管理、参数设置、数据处理、任务管理、偏移成像和用户管理等6个功能模块(见图4)。

1)资源管理。对所有资源节点进行实时有效地管理，能够实时监控各个计算节点的运行状态，获取节点的性能参数，完成计算节点的获取、添加和删除等操作。

2)参数设置。完成逆时偏移成像前所需要的各类参数设置，包括工区创建、偏移成像前基本参数设置和高级参数设置。

图4 叠前逆时偏移系统功能模块图

3)数据处理。完成偏移成像前数据的预处理，包括导入原始数据和导入速度模型等操作。

4)偏移成像。包括地震波场正演计算、检波点波场反传计算和互相关成像等操作。包括二维数据叠前逆时深度偏移成像和三维数据叠前逆时深度偏移成像功能，其中包括2D CPU RTM、2D GPU RTM、3D CPU RTM和3D GPU RTM。

5)任务管理。对逆时偏移成像任务的执行情况进行实时监测和管理，能够监管整体任务和各个计算节点任务的执行情况。包括任务添加、任务删除、任务重启、任务中断以及查看任务执行时间等操作，同时，也能实时监测出每个节点运行性能参数。

6)用户管理。主要包括用户管理、权限设置和日志查询等功能。用户管理包括用户的添加、删除和修改等操作，能够根据用户使用系统功能分配相应权限，根据权限不同，用户登录后显示的操作界面也不同，并能够记录用户操作日志。

3 系统主要功能实现

根据叠前逆时偏移数据处理流程与实现算法，通过研究可视化平台开发技术和叠前逆时偏移方法，在Linux环境下应用QT图形界面开发平台[12]，设计并实现了一套叠前逆时偏移成像系统。针对逆时偏移实现过程中二维CPU、二维GPU、三维CPU和三维GPU偏移算法和处理方法不相同，分别设计了二维 CPU逆时偏移、二维GPU逆时偏移、三维 CPU逆时偏移和三维 GPU逆时偏移等功能。

3.1 叠前逆时偏移成像实现

叠前逆时偏移成像主要包括2D CPU逆时偏移、2D GPU逆时偏移、3D CPU逆时偏移和3D GPU逆时偏移4个模块(见图5)。笔者以2D GPU逆时偏移模块为例描述系统的实现过程，其他模块操作类似：首先创建工区，导入炮集数据和速度模型文件，设置偏移成像的相关参数；选择2D GPU RTM处理方式和CPU计算节点，完成整体偏移成像。2D GPU逆时偏移具体操作如图6所示。

图5 2D CPU逆时偏移模块

图6 2D CPU逆时偏移成像

偏移成像界面中可以通过Edit按钮修改参数，Inport按钮载入参数卡文件，参数设定后，通过Export输出参数卡文件。所有参数设置完成后，点击Start按钮发送作业开始执行逆时偏移计算，作业执行过程中点击Refresh按钮可以更新进度条，查看作业执行进度。

3.2 多节点任务管理的实现

图7 多节点作业流程图

系统设计了多节点间任务监控功能，支持异构节点间作业的发送与终止。通过对节点文件信息的读取，识别节点设备名称；根据节点名与登录用户名信息生成唯一的shell文件，作业发送或终止命令发出时，各节点执行自身匹配的shell文件完成作业操作。多节点作业发送流程如图7所示。

在完成多节点作业发送功能基础上，添加了多节点作业进度监控功能。在多节点任务执行过程中，作业流会定期更新进度文件并保存在工区指定目录中。通过定期发送更新指令，系统从进度文件中读取作业进度信息，完成进度显示模块的更新。

4 系统应用与效果分析

4.1 系统应用可行性分析

为验证系统的可行性及算法集成后的有效性，试验中选取一个断层模型的200炮数据完成正演计算和偏移成像试验。试验节点为1个Intel E5-2630六核CPU处理器和Nvidia Tesla K10 GPU卡，试验具体参数如表1所示。

表1 系统软件测试

逆时偏移过程中采用时间二阶、空间十阶的有限差分数值算法完成地震波场正传和检波点波场反传计算，采用震源归一化互相关成像条件提取成像点，并采用拉普拉斯算子滤波法对图像进行去噪处理，得到最终成像效果。试验中的速度模型、逆时偏移成像后剖面和采用拉普拉斯算法去噪后剖面分别见图8(a)、(b)和(c)所示。

图8 逆时偏移系统测试

通过试验成像结果发现，系统成像结果与源代码方式成像结果基本一样，没有明显区别。同时，试验中源代码计算与系统计算的执行时间分别为134.28s和133.46s，系统执行时间略少。可见，系统完成逆时偏移从成像效果和执行时间等方面具有可行性，成像效果良好。

4.2 系统2D RTM成像效果分析

表2 试验参数

为验证系统的成像效果，选取二维盐丘(Salt-2D)模型进行叠前逆时偏移试验。试验环境为4个计算节点，每个节点配备1个Intel E5-2630六核CPU处理器和Nvidia Tesla K10 GPU卡，试验参数如表2所示。

在Salt-2D模型试验中，采用时间二阶、空间十阶的有限差分数值算法完成波场数值计算，采用互相关成像条件提取成像点，成像后采用拉普拉斯算子滤波法去噪。速度模型如图9所示，逆时偏移成像结果如图10所示。

从试验结果可以看出，笔者设计的逆时偏移系统能够有效的对地下结构精确成像，深层构造成像清楚，整体趋势与速度模型数据基本一致，地层分界明显，微幅构造成像清晰。

5 结语

笔者主要围绕逆时偏移系统设计与应用展开研究，通过研究逆时偏移成像算法和CPU/GPU并行加速技术，设计逆时偏移系统数据处理流程，并根据CPU和GPU程序代码、处理方法不同设计了2D CPU逆时偏移、2D GPU逆时偏移、3D CPU逆时偏移和3D GPU逆时偏移等功能模块，利用QT环境开发了基于CPU/GPU异构平台的叠前逆时偏移成像系统。试验结果表明系统具有可行性，与源代码方式相比成像结果一致，能够对复杂地下构造获得良好的成像效果。同时，利用GPU多线程实现了大规模并行计算，大大提高了偏移计算效率。

图9 Salt-2D速度模型

图10 Salt-2D逆时偏移结果

[1]陈可洋. 地震波逆时偏移方法研究综述[J]. 勘探地球物理进展，2010，33(3)：153～159.

[2] Liu Y，Chang X，Jin D，et al. Reverse time migration of multiples for subsalt imaging[J]. Geophysics，2011，76(5)：209～216.

[3] Costa J C，Silva Neto F A.Obliquity-correction imaging condition for reverse time migration[J]. Geophysics，2009，74(3)：57～66.

[4] 高新成，李春生.基于CPU/GPU协同加速叠前逆时偏移方法研究[J].陕西理工学院学报(自然科学版)，2016(2)：44～49.

[5] 石颖，柯璇，张莹莹. 逆时偏移边界条件与存储策略分析[J].地球物理学进展，2015(2)：581～585.

[6] Clapp R G. Reverse time migration with random boundaries[J]. 79th Annual International Meeting，SEG Expanded Abstracts，2009，28：2809～2813.

[7]Jiaze He，Fuh-Gwo Yuan.Damage identification for composite structures using a cross-correlation reverse-time migration technique[J]. Structural Health Monitoring，2015， 146：224～226.

[8] Ying Shi，Yanghua Wang. Reverse time migration of 3D vertical seismic profile data[J]. Geophysics，2016，81(1)：31～38.

[9] 石颖，陆加敏，柯璇，等. 基于GPU并行加速的叠前逆时偏移方法[J].东北石油大学学报，2012，36(4)：111～115.

[10] 刘守伟，王华忠，陈生昌，等. 三维逆时偏移GPU/CPU机群实现方案研究[J]. 地球物理学报，2013，56(10)：3487～3496.

[11] 王维红，柯璇，郭雪豹，等. 基于QT的Walkaway VSP逆时偏移软件开发研究[J]. 石油物探，2015，54(4)：452～458.

[12] 陈楠，祝媛媛，张光德，等. 基于QT的地震勘探可扩展平台研发与应用[J]. 华北地震科学，2013，32(4)：31～36.

[编辑] 辛长静

2016-09-18

国家自然科学基金项目(41574117)；黑龙江省教育科学规划重点课题(GJB1215013)

高新成(1979-)，男，博士，副教授，现主要从事地震数据处理与智能计算方面的教学与研究工作；E-mail：gxc@nepu.edu.cn。

TP399

A

1673-1409(2016)34-0015-06

[引著格式]高新成，石颖.基于CPU/GPU异构平台的叠前逆时偏移成像系统[J].长江大学学报(自科版)，2016,13(34)：15～20.