基于GPU集群的大规模三维有限差分正演模拟并行策略

引用本文

廉西猛, 张睿璇. 基于GPU集群的大规模三维有限差分正演模拟并行策略[J]. 物探与化探, 2015,39(3): 615-620.
LIAN Xi-Meng, ZHANG Rui-Xuan. Parallel strategy of large-scale 3D seismic forward by finite difference method on GPU cluster[J]. Geophysical and Geochemical Exploration, 2015,39(3): 615-620.
Doi:10.11720/wtyht.2015.3.31 复制到剪切板

Permissions

《物探与化探》编辑部

基于GPU集群的大规模三维有限差分正演模拟并行策略

廉西猛, 张睿璇

中国石化胜利油田分公司物探研究院, 山东东营 257022

作者简介: 廉西猛(1983-),男,山东济宁人,博士,工程师,主要从事地震勘探软件研发和地球物理方法研究。

收稿日期: 2014-10-09

基金: 中石化集团公司项目(P13074、P14083)

摘要

三维弹性波动方程有限差分正演模拟的低效率问题是导致该算法无法在大规模实际生产中应用的最重要的原因,使用GPU平台进行加速又面临GPU显存的限制。针对这一问题,提出了一种基于CPU/GPU的异构集群环境的并行加速解决方案。通过使用CPU/GPU协同并行模式和区域分解方法,实现了该算法的多级并行策略,并提出了非阻塞通讯、存储优化和基于MPI-IO的并行读写等方法,对程序的执行效率和存储使用进行了优化,为大规模实际生产应用提供了有效支撑。

关键词: 大规模三维正演; CPU/GPU集群; 区域分解; 存储优化; 并行读写

中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2015)03-0615-06

Parallel strategy of large-scale 3D seismic forward by finite difference method on GPU cluster

LIAN Xi-Meng, ZHANG Rui-Xuan

Shengli Geophysical Research Institute of SINOPEC,Dongying 257022,China

Abstract

Because of low efficiency,seismic forward simulation of 3D wave equation by finite difference method cannot be applied to real large-scale product.To deal with this problem,the authors present a parallel accelerating solution based on multi-GPUs heterogeneous cluster.By using CPU/GPU cooperation parallel mode and domain decomposition method,the authors carry out a multi-level parallel strategy of this algorithm.Furthermore,non-blocking MPI communications,storage optimization and parallel I/O mechanism using MPI-IO API are presented to optimize computation efficiency and memory usage.This parallel 3D forward algorithm can effectively support large-scale practical production.

Keyword: large-scale 3D seismic forward; CPU/GPU cluster; domain decomposition; storage optimization; parallel I/O

Show Figures

众所周知, 正演模拟在地震资料采集、处理和解释、储层建模、油气藏开发等环节, 都扮演着重要的角色, 其应用贯穿了石油勘探开发的整个流程。特别是大规模三维正演模拟, 随着勘探开发的深入, 其应用需求越来越巨大。

有限差分方法是波动方程地震波正演模拟的重要方法之一, 其实现简单, 精度较高, 因此成为目前使用最广泛的方法。三维波动方程有限差分正演模拟方法的理论基础已经十分成熟, 包括了声波, 弹性波、粘弹性以及各向异性、非均匀等各种波动方程, 普通网格, 交错网格以及旋转交错网格等多种网格类型, 都有大量的文献进行了深入的研究^{[1, 2, 3]}。然而理论的成熟与生产实用化之间还存在巨大的鸿沟— — 计算效率。在实际生产中, 需要进行模拟的区块往往十分庞大, 由此定义的网格剖分点数量也特别巨大。由于有限差分正演算法需要在空间每个网格点上都应用差分格式计算, 因此计算量十分庞大。经测试, 计算一个面积为1 010 m× 1 010 m, 采样点数为2 000的区块, 剖分网格大小为10 m× 10 m× 10 m, 使用CPU串行计算一个单炮的时间就超过24 h, 这样的效率在实际生产中是不能容忍的。

幸运的是, 有限差分方法具有很高的可并行性, 在进行时间步更新时, 各个网格点的计算是相互独立的, 都只依赖于已经计算出的上一次的若干个网格点的值, 因此各网格点的计算可以并行执行。因此, 许多文献都提出了基于CPU+GPU异构平台的优化加速技术, 如王延光^[4], 刘红伟^[5], 李博^[6]等对叠前逆时偏移有限差分算法在GPU设备上的并行实现策略进行了研究, 可是只采用了单个GPU模拟一个单炮的方式。GPU具有极强的计算能力, 适合于密集型的计算, 缺点是GPU的内存相比CPU较小, 无法满足大规模三维正演模拟的巨大内存需求。为了解决这一两难局面, 区域分解技术被广泛的应用。该技术的核心思想是将计算区域分割并分配到多个计算机节点上计算, 以达到分摊单个GPU内存压力的目的。Miché a^[7]针对三维弹性波常规网格高阶有限差分正演算法, 提出了多GPU并行策略, 使用了MPI非阻塞通讯实现了消息传递时间的隐藏。Miché a同时研究了该正演算法的CPML边界条件, 但并未讨论CPML边界条件在GPU上的实现方法。

Komatitsch^[8]应用区域分解方法研究了有限元方法在多GPU环境上的并行实现策略。龙桂华^[9]参考Miché a的区域剖分方法, 将研究对象扩展到三维交错网格有限差分正演模拟算法尺度三维地震波场模拟问题。刘守伟^[10]将区域分解方法应用到三维声波方程逆时偏移算法中, 根据GPU设备是否具备卡间直连(GPU Direct)通讯功能, 分别讨论了其区域分解方法的实现方式。针对吸收边界条件在GPU上的实现, 刘守伟选用了NPML(Nearly PML)边界条件。在该条件下, 吸收层的控制方程的二阶微分部分与内部区域一致, 从而避免了一部分分支判断, 一定程度上提高了GPU的实现效率。

笔者根据有限差分方法的特点, 在以上研究的基础上, 对多GPU集群并行策略进行了进一步的优化。使用区域分解方法和CPU+GPU的异构并行模式, 通过MPI+CUDA的多核异构并行编程方式, 实现了基于GPU集群的多级并行策略, 并提出了MPI非阻塞通讯、存储优化和基于MPI-IO的并行读写等方法对算法进行了优化。在下一节, 将详细地论述这一策略, 第3节中将通过应用实例验证策略的并行加速效果。

1 基于GPU集群的多级并行策略

基于弹性波方程高阶有限差分正演数值模拟算法, 使用的高阶有限差分格式为

$\begin{matrix} \begin{matrix} \frac{\partial f}{\partial x} = \frac{1}{Δx} \overset{N}{\sum_{n = 1}} C_{n}^{(N)} \{f [x + \frac{Δx}{2} (2 n - 1)] - \\ f [x - \frac{Δx}{2} (2 n - 1)]\} + O (Δ x^{2 N}) \end{matrix} \end{matrix}$

其中:差分系数 $\begin{matrix} C_{n}^{(N)} \end{matrix}$ 可以通过求解一个N阶方程组确定。

三维弹性波正演模块的主要流程如图1所示, 其核心算法主要包括单炮循环、时间步循环以及对空间网格结点的循环等多个级别的循环。各单炮的计算是各自独立的, 因此这一级循环可以进行并行, 目前的正演算法大都采用MPI实现了这一级别的并行, 多个CPU进程同时运算, 每个进程完成一个单炮的计算。时间步的循环是不能进行并行化的, 因为每个时间层的计算都依赖于其上一时间层的结果。依据有限差分正演算法的特点, 空间网格结点的循环包含了程序95%以上的计算量, 并且各节点的计算都是独立的, 因此适合使用GPU进行并行加速。但是GPU的显存较小, 而生产中的大规模正演的单炮模拟需要占用大量的内存, 单个GPU无法满足需求, 因此我们采用了区域分解方法, 即将要模拟的区域进行分割, 分成若干个子区域, 这样每个子区域上的模拟计算的内存需求都不会太大, 单个GPU可以满足。这种方案既分摊了内存压力, 又可以并发执行各个子区域上的计算。

	Figure Option View Download New Window
	图1 限差分正演算法CPU串行流程(N_x, N_y, N_z分别表示在空间三个方向上的网格点个数, N_t表示时间步个数)

根据以上分析, 基于异构平台, 结合正演算法的特点, 文中设计了如图2所示的多级并行策略。按照并行的层次和方式, 将其分成了单炮级、区域级和结点级三个级别, 分别使用了MPI进程组间、MPI进程组内进程和GPU线程进行实现。

	Figure Option View Download New Window
	图2 多级并行策略

1.1 单炮级并行

单炮级并行使用MPI进程组间并行来实现。按照每个CPU进程(对应一个GPU设备, 如果GPU设备具备卡间直连通讯功能, 可以对应多个可直连的GPU设备)负责计算一个子区域的原则, 根据子区域的个数确定每一个单炮计算需要的进程数, 然后将所有进程按此进程数分成若干个进程组, 将需要计算的所有炮平均分配到每个进程组上计算。每个进程组可以视为单个进程, 进程组之间相互独立, 互不依赖; 进程组内各进程的计算并不完全独立, 需要在通讯域内相互传递数据。组内进程协调合作, 共同完成一个单炮的计算。

1.2 区域级并行

区域级并行使用区域分解方法以及MPI进程组内并行来实现。文中以一个方向上的区域分解方法为例, 该方法可以推广到在三个方向上都进行分解。根据计算区域的大小, 确定分解子区域的个数, 保证每个子区域的计算所需内存小于GPU的显存。在对区域进行平均划分之后, 进程组中的每个进程负责一个子区域的运算。

但是子区域的计算并非完全独立并行执行的。由于差分格式的限制, 每个结点的计算都需要用到该结点上下、左右、前后六个方向的若干个点的波场值, 因此对于每个区域的边界附近的结点, 其计算需要用到相邻区域的结点的波场值。因此我们在每个子区域上附加一个辅助区域(如图3中A标记区域), 宽度等于差分格式阶数的1/2个网格间距。该区域用来接收相邻区域中B标记的部分传递来的波场值。图3给出了各个区域之间的数据传递示意。一旦区域间的数据传递完成, 各个子区域上的计算就可以相互独立地并发执行。

	Figure Option View Download New Window
	图3 MPI通讯域组内各进程间的数据交换(图中标记中的数字表示子区域的编号, A、B、C分别代表用于接收数据的辅助分块、用于发送数据的分块和剩余分块, R(L)表示分块位于子区域的右(左)侧)

	Figure Option View Download New Window
	图4 为隐藏消息传递时间设计的实现流程

数据传递的时间将会大大降低计算效率, 因此我们设计了如图4所示的方案, 让数据传递时间与计算时间同时进行, 从而隐藏数据传递的时间。执行流程为:首先在GPU上计算需要交换的部分(B标记部分)的波场值, 传回CPU后由CPU启动MPI非阻塞发送(接收)。启动后, 程序并不等待发送(接收)操作完成, 而是立即返回, 执行下一指令--在GPU进行子区域剩余部分(C标记部分)的计算, 而CPU同时在进行发送(接收)数据的任务, 这样CPU收发数据的时间与GPU计算C标记部分的时间就相互重叠, 从而达到隐藏数据传递时间的目的。

无论GPU设备是否具备卡间直连(GPU Direct)通讯功能, 此种数据传递方式都可以使用, 因此适合配置了多种不同系列GPU设备的集群环境。

区域分解引发的另一个问题是集群多节点对同一文件的读写访问。一般采用的方式(如图5左图所示)是设置一个节点专门负责读写文件。读文件时, 读写节点从文件读取所有数据, 然后根据需求将数据分发给其他计算节点。计算结束后, 读写节点从各个计算进程收集数据, 合并之后写出到文件。由于需要进行文件传递, 这样的读写方式效率较低。根据区域分解方法的特点, 我们使用了MPI-IO机制实现了多节点并行读写(如图5右图所示), 即不再设置专门的读写节点, 而是为每个节点在文件中指定其读写位置, 各个节点并行的从该位置读写数据, 这样可以省去数据传递的过程, 从而提升读写的效率。

	Figure Option View Download New Window
	图5 读写方式示意

1.3 结点级并行

结点级的并行依赖GPU线程实现。子区域上的计算主要是对区域内的每个网格结点应用差分格式, 其中包含大量的数学计算, 属于密集型计算任务, 正符合GPU擅长处理密集计算的优点。

GPU使用线程格(Grid)和并行线程块(Block)来管理线程, 目前CUDA仅支持二维的线程块数组。而有限差分方法对网格结点的循环计算涉及了三个维度, 在使用CUDA将算法迁移到GPU设备上时, 需要将三维空间沿x方向进行切片, 每个yz平面为一个切片。线程块Block中的每个线程完成切片上一个结点的计算。GPU每次并行计算一个切片, 沿x方向循环计算, 直到所有切片计算完成。

在吸收边界条件上, 我们使用了CPML边界条件^[11], 该边界条件使用了复频移技术, 可以有效地瞬逝波和掠射波。在CPML边界条件中, 内部区域和各个PML区域的计算所基于的微分方程形式都不相同, PML区域的微分方程比内部区域的微分方程多了若干了辅助变量。目前的处理方式是在每个时间步将各个PML区域与内部区域分开进行计算和更新, 这种方式下程序中存在许多分支语句, 不利于向GPU平台移植。Toivanen^[12]深入讨论了这种处理方式在GPU中实现策略, 提出了几种优化方案, 但这些方案或者分支过多引起线程发散, 或者各个区域串行实现, 计算效率都不够理想。笔者提出了一种基于存储优化的GPU实现方法。首先通过在内部区域设置零值辅助变量 , 统一内部区域和PML区域的微分方程形式, 然后使用存储优化方法来降低内存需求。这里的存储优化指的是内部计算区域中的引入的零值辅助变量并不全部存储, 而只为一个切片上的零值辅助变量开辟内存。较之前一种方法, 该方法具有显著的优势, 在只增加了少量内存需求的情况下, 达到了避免分支判断, 提高效率的目的。

2 应用实例

我们配置了一个CPU/GPU异构机群进行实例测试。使用了5个计算节点, 每个节点配置了12核CPU, 并载两个GPU卡设备。CPU均为Intel Xeon X5650, 主频为2.67 GHZ, 内存为24 GB; GPU设备均为Tesla M2090, 每个GPU的显存为6 GB。每一个单炮都设置了6条检波线, 检波线间距为100 m。每条线有20个检波点, 检波点间距为50 m, 炮间距为100 m, 采样点数为2 000。计算使用的网格大小为10 m× 10 m× 10 m, 时间步长为0.5 ms。若只使用单个CPU进程进行一个单炮的模拟需要约192.5 小时。

对单炮进行模拟时, 我们测试了将计算区域分解成不同数目的子区域的情况下的计算效果。在上述观测系统条件下, 如果不进行区域分解, 将需要约11.5GB的存储空间, 远远大于GPU设备的显存, 程序会因为内存溢出而终止。因此我们分别测试了将区域分解成2~10个子区域的情形, 使用与子区域个数相同的CPU进程数, 每个进程控制一个GPU。

图6展示了使用CPU串行和使用本文算法正演得到的炮记录, 可以发现二者是一致的。更进一步, 我们使用下面的误差公式计算得到二者的误差满足

$\begin{matrix} ε = \max_{1 < i < n_{t}} \frac{[\overset{n_{s}}{\sum_{j = 1}} (U_{cpu}^{ij} - U_{gpu}^{ij})^{2}]^{\frac{1}{2}}}{[\overset{n_{s}}{\sum_{j = 1}} (U_{cpu}^{ij})^{2}]^{\frac{1}{2}}} < 1.3584 e^{- 7}, \end{matrix}$

式中:n_t表示炮记录中的总道数, n_s表示每一道的采样点数。此误差在允许误差的范围内, 不会对结果产生影响。

图7中比较了CPU内存、GPU内存、检波线方向网格点数、计算时间和加速比等参数, 为了研究其变化趋势, 我们将GPU内存的数值除以10, 加速比(此处等于T_i/T_i-₁, T_i为子区域个数为i时的计算时间)的数值乘上了1 000, 将所有参数绘在同一张图中。从图中可以看到, 随着分解数目的增多, 每个进程计算区域的大小在不断减小, 因此所需的CPU内存和GPU显存也持续递减, 达到了分摊内存的目标。这同时使得每个GPU内核的计算时间也相应的减少, 因此总的计算时间也减少, 这一点从加速比大于1也可知。但是加速比并不与子区域个数的增加成比例, 这是因为区域分解导致了数据传递时间的增加, 抵消了一部分GPU内核计算时间的减少。随着子区域个数的增加(如图6中子区域数目增加到8个), 这部分数据传递时间也显著增加, 增加量逐渐超过GPU内核计算时间的减少量, 进而导致总的计算时间不降反升, 加速比也开始小于1, 并逐渐降低。因此子区域数目过多, 会导致GPU利用率降低, 不能发挥GPU的优势, 从而影响计算效率。

	Figure Option View Download New Window
	图6 正演模拟结果对比

	Figure Option View Download New Window
	图7 使用不同子区域数目的测试实例对比

3 结论

为了解决三维正演模拟算法在大规模实际生产的应用问题, 笔者提出了一套解决方案。首先使用CPU/GPU异构并行方法提高计算效率。但GPU显存较小, 无法计算大尺度问题, 于是使用区域分解方法, 将巨大的内存需求分摊到多个GPU设备上。区域分解方法又增加了数据传递时间, 并导致了多节点竞争读写同一文件的问题, 我们采用MPI非阻塞消息收发和MPI-IO并行读写方法将它们解决。最后通过存储优化方法, 很好的解决了GPU上CPML边界条件的处理问题。这样形成的具有三级并行策略的三维正演算法能够为大规模实际生产应用提供有力支持。

致谢:感谢单联瑜教授和隋志强教授对本文研究成果的指导和帮助。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Xia F, Dong L G, Mz Z T. The numerical modeling of 3-D elastic wave equation using a high-order, staggered-grid, finite difference scheme[J]. Applied Geophysics, 2004, 1(1): 38-41. [本文引用:1]
[2]	张文生, 宋海斌. 三维正交各向异性介质三分量高精度有限差分正演模拟[J]. 石油地球物理勘探, 2001, 36(4): 422-432. [本文引用:1]
[3]	杨仁虎, 常旭, 刘伊克. 基于非均匀各向同性介质的黏弹性波正演数值模拟[J]. 地球物理学报, 2009, 52(9): 2321-2327. [本文引用:1]
[4]	王延光, 匡斌. 起伏地表叠前逆时深度偏移与并行实现[J]. 石油地球物理勘探, 2012, 47(2): 266-273. [本文引用:1]
[5]	刘红伟, 李博, 刘洪, 等. 地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报, 2010, 53(7): 1725-1733. [本文引用:1]
[6]	李博, 刘红伟, 刘国峰, 等. 地震叠前逆时偏移算法的CPU/GPU实施对策[J]. 地球物理学报, 2010, 53(12): 2938-2943. [本文引用:1]
[7]	Michéa D, Komatitsch D. Accelerating a 3D finite-difference wave propagation code using GPU graphics cards[J]. Geophys J Int, 2010, 182(1): 389-402. [本文引用:1]
[8]	Komatitsch D, Erlebacher G, Göddeke D, et al. High-order finite-element seismic wave propagation modeling with MPI on a large GPUcluster[J]. J Comput Phys, 2010, 229(20): 7692-7714. [本文引用:1]
[9]	龙桂华, 赵宇波, 李小凡, 等. 三维交错网格有限差分地震波模拟的GPU集群实现[J]. 地球物理学进展, 2011, 26(6): 1938-1949. [本文引用:1]
[10]	刘守伟, 王华忠, 陈生昌, 等. 三维逆时偏移GPU/CPU机群实现方案研究[J]. 地球物理学报, 2013, 56(10): 3487-3496. [本文引用:1]
[11]	Roden J A, Gedney S D. Convolution PML ( CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media[J]. Microwave and Optical Technology Letters, 2000, 27: 334-339. [本文引用:1]
[12]	Toivanen J I, Stefanski T P, Kuster N, et al. Comparison of CPML implementations for the GPU-accelerated FDTD solver[J]. Progress In Electromagnetics Research M, 2011, 19: 61-75. [本文引用:1]

2004

0.0

... 三维波动方程有限差分正演模拟方法的理论基础已经十分成熟,包括了声波,弹性波、粘弹性以及各向异性、非均匀等各种波动方程,普通网格,交错网格以及旋转交错网格等多种网格类型,都有大量的文献进行了深入的研究^[1,2,3] ...

2001

0.0

张文生, 宋海斌. 三维正交各向异性介质三分量高精度有限差分正演模拟[J]. 石油地球物理勘探, 2001, 36(4): 422-432.

正用有限差分法进行了三维正交各向异性介质中的三分量波场模拟,包括VSP记录和共炮点记录,并合成了三维各向异性模型的零偏移距地震记录,压制了频散,对三维吸收边界条件和稳定性条件作了推导,对各向异性波场分离作了初步探讨,从而使得该方法在三维各向异性介质中的波场分析和解释中成为一种有效方法.

2009

0.0

杨仁虎, 常旭, 刘伊克. 基于非均匀各向同性介质的黏弹性波正演数值模拟[J]. 地球物理学报, 2009, 52(9): 2321-2327.

In this paper, the authors present simplified viscoelastic wave equations in heterogeneous and isotropic media based on generalized linear viscoelasticity. Viscoelastic Lame constants are introduced and Lame difference matrix is put forward. Lame difference matrix that has the similar form with physical property matrix holds a certain quantitative relation with elastic and viscoelastic physical property matrices. It makes the viscoelastic equations more concise. The viscoelastic modeling in 2-D SEG/EAGE salt dome model is carried out by using the simplified viscoelastic equations. A comparison of the viscoelastic single-shot records derived from acoustic equations, Carcione's equations and the simplified viscoelastic equations is made. The results show that the simplified viscoelastic equations not only have the better effects of single-shot records, but also have a better computational efficiency than Carcione's viscoelastic equations and it agrees well with the theoretical analysis.

本文基于广义线性黏弹性体，通过引入黏弹性拉梅常量，提出拉梅差异矩阵，给出了非均匀各向同性介质的黏弹性波简化方程.拉梅差异矩阵具有和物性矩阵相似的形式，与弹性、黏弹性物性矩阵之间具有特定的数量关系，体现了黏弹性与弹性物性参数之间的差异程度，可以大大简化方程.并利用该简化方程对SEG/EAGE二维盐丘模型进行了正演数值模拟，将模拟得到的单炮记录和声波、Carcione黏弹性波方程模拟得到的单炮记录进行了对比分析.结果表明该简化方程模拟得到的单炮记录不仅具有较好的效果，而且比Carcione黏弹性波方程计算效率高，与理论分析基本相符.

2012

0.0

王延光, 匡斌. 起伏地表叠前逆时深度偏移与并行实现[J]. 石油地球物理勘探, 2012, 47(2): 266-273.

本文系统讨论了逆时偏移中常用的波动方程求解方法，并简述了其各自的特点和优势。文中提出用逆时偏移实现基于起伏地表的波动方程叠前深度偏移，理论模型数据的试算结果与单程波方程和Kirchhoff 积分方法结果进行对比，展示了逆时偏移的优越性。自行研发的基于CPU/GPU平台的并行软件，推动了逆时偏移算法实用化的运行，在应用中取得了良好效果。

... 因此,许多文献都提出了基于CPU+GPU异构平台的优化加速技术,如王延光^[4],刘红伟^[5],李博^[6]等对叠前逆时偏移有限差分算法在GPU设备上的并行实现策略进行了研究,可是只采用了单个GPU模拟一个单炮的方式 ...

2010

0.0

刘红伟, 李博, 刘洪, 等. 地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报, 2010, 53(7): 1725-1733.

Pre-stack reverse time migration (RTM) is a very useful tool for seismic imaging. But it has not been widely used because of the highly intensive computation cost, imaging noise and massy memory demand. In this paper, we illustrate the implementation process of RTM and analyze the stability condition and dispersion relation of finite difference (FD) method. For the problem of intensive computation cost, we use the Graphic Processing Unit (GPU) architecture to realize RTM and get an order of magnitude higher speedup ratio compared to the traditional CPU architecture. Compared to the one way wave equation migration methods, RTM does not have the imaging dip limit and the imaging effect is significantly improved. The test on vertical fault and BP models proves this conclusion. The problems of imaging noise removing and massy memory demand will be discussed in other papers.

叠前逆时偏移技术是解决地震成像问题的有力工具,但由于计算量大、成像噪音以及存储量大等原因没有得到广泛的应用.本文给出了逆时偏移的实现过程,分析了高阶有限差分格式的稳定性与频散关系.针对叠前逆时偏移计算量大的问题,使用图形处理器(Graphic Processing Unit,简称GPU)实现算法加速,比传统的CPU计算速度提高了一个数量级.文中对理论模型进行了计算,并与单程波偏移方法做比较,结果表明:叠前逆时偏移有效突破了成像倾角限制,对垂直断层、盐丘空腔内幕等特殊构造成像效果均有显著提高.本文尚未涉及成像噪音去除以及存储量等问题,笔者将另文阐述.

2010

0.0

李博, 刘红伟, 刘国峰, 等. 地震叠前逆时偏移算法的CPU/GPU实施对策[J]. 地球物理学报, 2010, 53(12): 2938-2943.

Comparing with one-way wave migration algorithm, reverse time migration (RTM) is more attractive because of the theory advantages. Two-way wave equation has been used to extrapolate wave field in RTM, instead of separating the up going wave and down going wave. However, due to the large amount of computation and I/O, RTM is most time-consuming in industrial applications. In this article, we analyze several computational strategies and propose our method, which uses CPU/GPU as computational core and builds random velocity boundary,for solving I/O problem and computational efficiency problem. In the actual test, it has been proved that this method can largely decrease storage memory units and improve computational efficiency.

相较于单程波偏移算法而言,逆时偏移成像方法以其物理基础为依托优势,几十年来一直备受国内外地球物理学家的青睐.目前的逆时偏移(RTM)若直接采用双程波动方程进行延拓,尽管可以回避上下行波的分离处理,然就已有算法而言,其计算量和I/O(输入/输出)量却是最大的.针对此问题,本文在分析现行逆时偏移的多种算法基础上,提出利用CPU/GPU(中央处理器/图形处理器)作为数值计算核心,建立随机边界模型,从而克服存储I/O难题和提高计算效率.在实际的数据测试中,本文的方法可以大幅度的提高计算效率和减少存储单元,从而促使其高效地应用于生产实际.

2010

0.0

... a^[7]针对三维弹性波常规网格高阶有限差分正演算法,提出了多GPU并行策略,使用了MPI非阻塞通讯实现了消息传递时间的隐藏 ...

2010

0.0

... Komatitsch^[8]应用区域分解方法研究了有限元方法在多GPU环境上的并行实现策略 ...

2011

0.0

龙桂华, 赵宇波, 李小凡, 等. 三维交错网格有限差分地震波模拟的GPU集群实现[J]. 地球物理学进展, 2011, 26(6): 1938-1949.

As an efficient numerical method for seismic modeling, finite difference has been widely used in computation-intensive waveform inversion and reverse-time migration. Although simple and fast, it is still hard to alleviate giant computation cost in three-dimensional seismic modeling that greatly restricts the industrial application of pre-stack reverse-time migration and inversion. The GPU general-purpose computation technology that is well known for its inherent data parallelism is expected to break through this limitation and makes large-scale three-dimensional reverse-time migration and inversion possible in application. In this paper, we present a method to configure the discrete grids on GPU thread blocks that optimally maps the global memory address on global video memory into inner shared memory to achieve maximum efficiency. The tests applied in different scale models show that we can achieve 79x~108x speedup ratio when compared to traditional single-core CPU. By virtue of domain decomposition, we use coarse-grained algorithm to partition the large-scale geological model along the Z-axis and compute them using hybrid MPI and CUDA codes to realize the large-scale seismic modeling which can hardly be implemented on single GPU. In the experiments, we also analyze some key factors that affect the parallel efficiency of GPU cluster. The methodology we present and the problems we discuss can be used as a basis to promote the industrial application of pre-stack reverse-time migration and inversion that depend heavily on large-scale three-dimensional seismic modeling.

有限差分实现简单、速度快,作为地震波场模拟一种有效数值方法,被广泛用于正演计算密集的波形反演和逆时偏移中.三维地震波正演模拟计算量大,一直以来制约着三维叠前逆时偏移和反演的工业化应用,GPU通用计算技术的产生及其内在的数据并行性有望改变这一现状.本文通过分析三维交错网格有限差分方法在GPU上的实施,利用片内共享存储器实现了三维地震波数值模拟的高效算法,取得了较单核CPU快79x~108x的加速比;通过区域分解技术将单GPU上不能计算的地质体模型沿Z轴方向进行粗粒度分解,采用消息传递接口交换边界数据,运用MPI+CUDA的方式实现了大尺度三维地震波场模拟,并着重分析了影响GPU并行计算效率的一些关键因素.大尺度三维地震波场模拟的加速实现,为促进叠前逆时偏移和波形反演技术的工业化转化提供了可能,因此具有重要的研究意义.

... 龙桂华^[9]参考Mich#cod#x000e9 ...

2013

0.0

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

Prestack reverse time migration(RTM) is currently the most precise seismic imaging method. It needs appropriate implementation strategies and efficient computation platforms, because of the needs for large amount of calculation and large storage. Based on high order finite difference RTM, this paper focuses on solving the questions of the lack of video memory and artificial boundary, which must be solved for Graphic Processing Unit (GPU) RTM. Using the domain decomposition technique, we can efficiently realize three dimensional (3D) RTM of any production scale. It eliminates the limitation of GPU memory size. The control equation of the traditional perfectly matched layer (PML) boundary condition is different to that of the inner area. It is not suitable for GPU calculation. This paper realizes the Nearly PML (NPML) boundary condition on GPU. It makes high order finite difference calculation eliminate the needs of conditional branches. And the amount of auxiliary wavefields in boundary area is small. Therefore, wavefields propagation using NPML is efficient on GPU. 3D theoretical data and real field data RTM results prove the correctness of this paper's methods.

叠前逆时偏移是当前最为准确的地震成像方法,由于计算量大、存储量大等原因需要合适的实现策略和高效的计算平台.本文以高阶有限差分逆时偏移为基础,重点讨论了在GPU上实现需要解决的显存不足问题和人工边界问题.利用区域分解技术可以在当前GPU上高效地实现任意生产规模的三维逆时偏移成像,不会受到GPU显存规模的制约.常规最佳匹配层边界条件边界区域控制方程与内部区域差异较大,不适于GPU高速运算.本文在GPU上实现近似最佳匹配层(NPML)边界条件,使得高阶有限差分计算不需要分支判断,边界区域辅助波场的存储量也较低,保证了在GPU上进行波场传播的高效性.三维理论数据和实际资料成像结果表明了本文方法的正确性.

... 刘守伟^[10]将区域分解方法应用到三维声波方程逆时偏移算法中,根据GPU设备是否具备卡间直连(GPU Direct)通讯功能,分别讨论了其区域分解方法的实现方式 ...

2000

0.0

... 在吸收边界条件上,我们使用了CPML边界条件^[11],该边界条件使用了复频移技术,可以有效地瞬逝波和掠射波 ...

2011

0.0

... Toivanen^[12]深入讨论了这种处理方式在GPU中实现策略,提出了几种优化方案,但这些方案或者分支过多引起线程发散,或者各个区域串行实现,计算效率都不够理想 ...