CN106324667A - 基于gpu的三维地震波场模拟计算方法及系统 - Google Patents

Abstract

公开了一种基于GPU的三维地震波场模拟计算方法及系统。该方法可以包括以下步骤：调用GPU核，载入前一时间的X方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的X方向上偏导数；调用GPU核，载入前一时间的Y方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Y方向上偏导数；调用GPU核，载入前一时间的Z方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Z方向上偏导数；以及对X、Y、Z三个方向上计算的偏导数求和得到当前位置上的散度，以完成当前时间波场的计算。调用GPU核计算X、Y、Z方向的偏导数是没有先后之分，而且在GPU资源许可的情况下可以同时进行。

Description

基于GPU的三维地震波场模拟计算方法及系统

技术领域

本公开涉及石油勘探应用领域，更具体地，涉及一种基于GPU的三维地震波场模拟计算方法及系统。

背景技术

在石油勘探应用领域，对地震资料的成像是至关重要的，其中，三维地震波场的模拟计算是现代地震波成像技术的核心。在现有技术中，中央处理器(CPU，Central Processing Unit)以及由成千上万个CPU组成的高性能计算机集群是用于实现地震波场模拟的主要工具。然而，近年来，采用图形处理器(GPU，Graphic Processing Unit)协同CPU来实现地震波场模拟和成像已经成为当前提高三维地震波场模拟和成像计算效率的发展主题方向。

发明人发现，现有的基于GPU的三维地震波场模拟计算方法通常采用三维空间偏导数(Laplacian算子)的整体实现方案，然而该计算过程中对每个前一个时间波场计算数据的三个方向相邻点的内存访问通常是非连续的，这导致GPU无法对三个方向相邻点实现合并访问，不能够充分发挥GPU的计算性能，效率较低。因此，有必要开发一种高效的基于GPU的三维地震波场模拟计算方法和系统。

公开于本公开背景技术部分的信息仅仅旨在加深对本公开的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。

发明内容

本公开提供了一种基于GPU的三维地震波场模拟计算方法及系统，其能够分三次调用GPU核分别对任何一个当前波场计算数据在空间X、Y、Z三个方向进行偏导数计算，使每次GPU调用都实现合并访问，从而提高了计算效率。

根据本公开的一个方面，提出了一种基于GPU的三维地震波场模拟计算方法，该方法可以包括以下步骤：调用GPU核，载入前一时间的X方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的X方向上偏导数；调用GPU核，载入前一时间的Y方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Y方向上偏导数；调用GPU核，载入前一时间的Z方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Z方向上偏导数；以及对X、Y、Z三个方向上计算的偏导数求和得到当前位置上的散度，以完成当前时间波场的计算。

根据本公开的另一方面，提出了一种基于GPU的地震波场模拟计算系统，该系统可以包括：用于调用GPU核，载入前一时间的X方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的X方向上偏导数的单元；用于调用GPU核，载入前一时间的Y方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Y方向上偏导数的单元；用于调用GPU核，载入前一时间的Z方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Z方向上偏导数的单元；以及用于对X、Y、Z三个方向上计算的偏导数求和得到当前位置上的散度，以完成当前时间波场的计算的单元。

本公开的方法和装置具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施方案中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方案中进行详细陈述，这些附图和具体实施方案共同用于解释本公开的特定原理。

附图说明

通过结合附图对本公开的示例性实施例进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开的示例性实施例中，相同的参考标号通常代表相同部件。

图1示出了适于用来实现本公开实施例的示例性计算机系统/服务器12的框图。

图2示出了根据本公开的一个实施例的基于GPU的三维地震波场模拟计算方法的流程图。

图3示出了现有技术的CPU调用GPU对波场计算数据同时在空间的X、Y、Z三个方向上进行偏导数计算的示意图。

图4示出了根据本公开的一个实施例的CPU分别调用三次GPU核对波场计算数据在空间的X、Y、Z三个方向上进行偏导数计算的示意图。

图5示出了根据本公开的一个实施例的计算所有时间的波场计算数据以完成波场计算的示意图。

图6A和图6B示出了采用Micikevcius方法和本公开的方法的结果对比的示意图。

图7示出了本公开相对于Micikevcius方法计算时间比率。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施例。虽然附图中显示了本公开的优选实施例，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

图1示出了适于用来实现本公开实施例的示例性计算机系统/服务器12的框图。图1显示的计算机系统/服务器12仅仅是一个示例，不应对本公开实施例的功能和使用范围带来任何限制。

如图1所示，计算机系统/服务器12以通用计算设备的形式表现。计算机系统/服务器12的组件可以包括但不限于：一个或者多个处理器或者处理单元16，系统存储器28，连接不同系统组件(包括系统存储器28和处理单元16)的总线18。

总线18表示几类总线结构中的一种或多种，包括存储器总线或者存储器控制器，外围总线，图形加速端口，处理器或者使用多种总线结构中的任意总线结构的局域总线。举例来说，这些体系结构包括但不限于工业标准体系结构(ISA)总线，微通道体系结构(MAC)总线，增强型ISA总线、视频电子标准协会(VESA)局域总线以及外围组件互连(PCI)总线。

计算机系统/服务器12典型地包括多种计算机系统可读介质。这些介质可以是任何能够被计算机系统/服务器12访问的可用介质，包括易失性和非易失性介质，可移动的和不可移动的介质。

系统存储器28可以包括易失性存储器形式的计算机系统可读介质，例如随机存取存储器(RAM)30和/或高速缓存存储器32。计算机系统/服务器12可以进一步包括其它可移动/不可移动的、易失性/非易失性计算机系统存储介质。仅作为举例，存储系统34可以用于读写不可移动的、非易失性磁介质(图1未显示，通常称为“硬盘驱动器”)。尽管图1中未示出，可以提供用于对可移动非易失性磁盘(例如“软盘”)读写的磁盘驱动器，以及对可移动非易失性光盘(例如CD-ROM，DVD-ROM或者其它光介质)读写的光盘驱动器。在这些情况下，每个驱动器可以通过一个或者多个数据介质接口与总线18相连。存储器28可以包括至少一个程序产品，该程序产品具有一组(例如至少一个)程序模块，这些程序模块被配置以执行本公开各实施例的功能。

具有一组(至少一个)程序模块42的程序/实用工具40，可以存储在例如存储器28中，这样的程序模块42包括——但不限于——操作系统、一个或者多个应用程序、其它程序模块以及程序数据，这些示例中的每一个或某种组合中可能包括网络环境的实现。程序模块42通常执行本公开所描述的实施例中的功能和/或方法。

计算机系统/服务器12也可以与一个或多个外部设备14(例如键盘、指向设备、显示器24等)通信，还可与一个或者多个使得用户能与该计算机系统/服务器12交互的设备通信，和/或与使得该计算机系统/服务器12能与一个或多个其它计算设备进行通信的任何设备(例如网卡，调制解调器等等)通信。这种通信可以通过输入/输出(I/O)接口22进行。并且，计算机系统/服务器12还可以通过网络适配器20与一个或者多个网络(例如局域网(LAN)，广域网(WAN)和/或公共网络，例如因特网)通信。如图所示，网络适配器20通过总线18与计算机系统/服务器12的其它模块通信。应当明白，尽管图中未示出，可以结合计算机系统/服务器12使用其它硬件和/或软件模块，包括但不限于：微代码、设备驱动器、冗余处理单元、外部磁盘驱动阵列、RAID系统、磁带驱动器以及数据备份存储系统等。

实施例1

图2示出了根据本公开的实施例的基于GPU的三维地震波场模拟计算方法的流程图。

在该实施例中，该方法可以包括以下步骤：步骤201，调用GPU核，载入前一时间的X方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的X方向上偏导数；步骤202，调用GPU核，载入前一时间的Y方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Y方向上偏导数；步骤203，调用GPU核，载入前一时间的Z方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Z方向上偏导数；以及步骤204，对X、Y、Z三个方向上计算的偏导数求和得到当前位置上的散度，以完成当前时间波场的计算。

该实施例通过分三次调用GPU分别对当前波场计算数据在空间的X、Y、Z三个方向上进行偏导数计算，使每次GPU调用都实现合并访问，从而提高计算效率。调用GPU核计算X、Y、Z方向的偏导数是没有先后之分，而且在GPU资源许可的情况下可以同时进行。

波场计算

在一个示例中，可以通过本领域技术人员已知的任何现有技术手段来获取波场计算数据。例如，可以对地震资料进行预处理，以获得波场计算数据。其中，可以使用当前波场计算数据表示当前正在计算的波场计算数据。

在三维地震波场模拟时，计算当前时间在当前位置的的波场值可以通过计算前一个时间的波场在该位置的散度来得到。在计算波场在某个位置的散度时，可以通过CPU分三次调用GPU核对当前波场计算数据分别在空间的X、Y、Z三个方向上进行偏导数计算。然后，可以对X、Y、Z三个方向上计算的偏导数求和得到该位置上的散度，以完成当前时间波场的计算。其中，调用GPU核计算X、Y、Z方向的偏导数可以没有先后之分，而且在GPU资源许可的情况下可以同时进行。

为了便于理解本实施例的原理和技术效果，图3和图4给出了对波场计算数据在X、Y、Z三个方向上进行偏导数计算的一个示例，其中图3示意了现有技术的方法，图4示意了本公开的实施例的方法。本领域技术人员应理解，该示例仅为了便于理解，而不意在以任何方式限制本发明。例如，本领域技术人员应理解，虽然示例中给出了利用两个相邻数据来求某一方向上的偏导数，但也可以利用其它数量的相邻数据来计算。并且本领域技术人员可以根据需要采用任何已知的具体方式来利用相邻数据求偏导数。此外，本领域技术人员还应理解，虽然该示例中在计算每个方向的偏导数时调用了一个GPU，然而在分布式计算环境下，也可以分别调用一组GPU来计算每个方向的偏导数。

图3示出了现有技术的CPU调用GPU核对波场计算数据同时在空间的X、Y、Z三个方向上进行偏导数计算的示意图。如图3所示，在现有的GPU整体访问内存的方案中，需要同时计算三个方向的偏导数(Laplacian算子)，对点0的单次偏导数计算涉及三个方向的相邻点(1，2，3，4，5，6)，而由于相邻点5和6在内存中与其它点相距甚远，导致无法实现GPU的合并访问，从而导致整体计算效率受到限制。此外，同时计算三个方向的偏导数，很容易使得每一个线程的寄存器使用数量达到极限，进一步限制了此类方法的计算效率和计算能力，特别是在采用高阶方法的时候，例如16阶有限差分方法。

图4示出了根据本公开的一个实施例的CPU分别调用三次GPU核对波场计算数据在空间的X、Y、Z三个方向上进行偏导数计算的示意图。如图4所示，可以分别在X、Y、Z三个方向上计算偏导数。在对X方向计算偏导数时，对点0的偏导数计算可以仅使用前一时间的X方向的相邻点(5，6)的波场数据；在对Y方向计算偏导数时，对点0的偏导数计算可以仅使用前一时间的Y方向的相邻点(3，4)的波场数据；在对Z方向计算偏导数时，对点0的偏导数计算可以仅使用前一时间的Z方向的相邻点(1，2)的波场数据。然后，可以对X、Y、Z三个方向上的计算结果求和，从而完成波场计算。通过上述步骤，在每个方向上都可以实现GPU全局内存的合并访问，从而可以提高计算效率。同时，每次GPU内核计算，仅限于对一个方向的偏导数执行，使得每一个线程的寄存器使用数量大幅度降低，从而消除了因寄存器数目有限带来的计算效率限制。尤其是在计算高阶偏导数的时候，本公开的方法可以显著提高计算效率和计算能力。

如上所述，与现有技术相比，根据本公开的实施例的基于GPU的三维地震波场模拟计算方法可以显著提高波场计算的效率。

在一个示例中，可以通过对所有时间的波场计算数据的计算来实现整个波场计算过程。图5示出了根据本公开的一个实施例的计算所有时间的波场计算数据以完成波场计算的示意图。其中，可以将开始计算时的时间it值设为零，在对当前时间it的波场计算数据进行计算之后，可以使当前时间it增加一个时间增量，即前进到下一时间，例如使it＝it+1(其中1可以代表一个时间增量，例如1个单位采样时间)，然后可以对下一时间的波场计算数据进行计算。当时间it达到总时间nt时，可以停止计算。通过对所有时间的波场计算数据的计算，可以实现整个波场计算过程。

其中，时间增量的值可以依据本领域技术人员的经验设定为任意合理的值，例如1个单位采样时间。本领域技术人员应理解，以上所列出的仅仅是示例，并非是穷举了所有的可以设定的时间增量。

在一个示例中，可以基于波场计算，对每个时间计算得到的地震波场数据进行保存，从而获得整个地震波场。

应用示例

以下给出了本公开的实施例的一个应用示例，来说明本公开的实施例在保证计算结果的同时，提高了计算效率的有益效果。本领域技术人员应理解该应用示例仅为了便于理解本公开，其中的数值及其他细节仅为示例性的，而非意在限制本公开。

图6A和图6B示出了分别采用传统的Micikevcius方法和本公开的方法，在其他计算参数相同的条件下，获得的地震波场结果图。在图6A和图6B所示的示例中，采用256×256×256的三维计算网格，在单个GPU上进行了2000步的地震波场模拟计算。图7是在其他计算参数相同的条件下，对于不同网格数相对于传统的Micikevcius方法的计算时间比率图，横坐标为整体三维网格数，纵坐标为相对于传统的Micikevcius方法的计算时间的比率。如图7所示，在整体三维网格数为1x10⁷的情况下，本公开的方法的计算时间只占传统的Micikevcius方法的67％，在整体三维网格数为9x10⁷的情况下，本公开的方法的计算时间只占传统的Micikevcius方法的58％。

从图6和图7可以看到，采用CPU调用三次GPU核分别对当前波场计算数据在空间X、Y、Z三个方向进行偏导数计算，其计算结果与传统的CPU调用GPU对当前波场计算数据在X、Y、Z三个方向整体进行偏导数计算相同，但是计算效率得到了很大的提高。

本领域技术人员应理解，上面对本公开的实施例的描述的目的仅为了示例性地说明本公开的实施例的有益效果，并不意在将本公开的实施例限制于所给出的任何示例。

实施例2

根据本公开的实施例，提供了一种基于GPU的地震波场模拟计算系统。该系统可以包括：用于调用GPU核，载入前一时间的X方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的X方向上偏导数的单元；用于调用GPU核，载入前一时间的Y方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Y方向上偏导数的单元；用于调用GPU核，载入前一时间的Z方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的Z方向上偏导数的单元；以及用于对X、Y、Z三个方向上计算的偏导数求和得到当前位置上的散度，以完成当前时间波场的计算的单元。

该实施例通过分三次调用GPU分别对当前波场计算数据在空间的X、Y、Z三个方向上进行偏导数计算，使每次GPU调用都实现合并访问，从而提高计算效率。

在一个示例中，可以通过对所有时间的波场计算数据的计算来实现整个波场计算过程。

在一个示例中，可以对每个时间计算得到的地震波场数据进行保存，从而获得整个地震波场。

本领域技术人员应理解，上面对本公开的实施方案的描述的目的仅为了示例性地说明本公开的实施方案的有益效果，并不意在将本公开的实施方案限制于所给出的任何示例。

本公开可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本公开的各个方面的计算机可读程序指令。

计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是――但不限于――电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如，通过光纤电缆的光脉冲)、或者通过电线传输的电信号。

这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令，并转发该计算机可读程序指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。

用于执行本公开操作的计算机程序指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，所述编程语言包括面向对象的编程语言—诸如Smalltalk、C++等，以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络—包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。在一些实施例中，通过利用计算机可读程序指令的状态信息来个性化定制电子电路，例如可编程逻辑电路、现场可编程门阵列(FPGA)或可编程逻辑阵列(PLA)，该电子电路可以执行计算机可读程序指令，从而实现本公开的各个方面。

这里参照根据本公开实施例的方法、装置(系统)和计算机程序产品的流程图和/或框图描述了本公开的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。

这些计算机可读程序指令可以提供给通用计算机、专用计算机或其它可编程数据处理装置的处理器，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理器执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。

也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。

附图中的流程图和框图显示了根据本公开的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，所述模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

以上已经描述了本公开的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。

Claims (6)

1.一种基于GPU的三维地震波场模拟计算方法，所述方法包括以下步骤：

调用GPU核，载入前一时间的X方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的X方向上偏导数；

调用GPU核，载入前一时间的Y方向上与所述当前位置相邻的波场数据，来计算当前时间在当前位置的Y方向上偏导数；

调用GPU核，载入前一时间的Z方向上与所述当前位置相邻的波场数据，来计算当前时间在当前位置的Z方向上偏导数；以及

对X、Y、Z三个方向上计算的偏导数求和得到所述当前位置上的散度，以完成当前时间波场的计算。

2.根据权利要求1所述的一种基于GPU的三维地震波场模拟计算方法，其中，通过对所有时间的波场计算数据的计算来实现整个波场计算过程。

3.根据权利要求2所述的一种基于GPU的三维地震波场模拟计算方法，其中，对每个时间计算得到的地震波场数据进行保存，从而获得整个地震波场。

4.一种基于GPU的三维地震波场模拟计算系统，所述系统包括以下单元：

用于调用GPU核，载入前一时间的X方向上与当前位置相邻的波场数据，来计算当前时间在当前位置的X方向上偏导数的单元；

用于调用GPU核，载入前一时间的Y方向上与所述当前位置相邻的波场数据，来计算当前时间在当前位置的Y方向上偏导数的单元；

用于调用GPU核，载入前一时间的Z方向上与所述当前位置相邻的波场数据，来计算当前时间在当前位置的Z方向上偏导数的单元；以及

用于对X、Y、Z三个方向上计算的偏导数求和得到所述当前位置上的散度，以完成当前时间波场的计算的单元。

5.根据权利要求4所述的一种基于GPU的三维地震波场模拟计算系统，其中，通过对所有时间的波场计算数据的计算来实现整个波场计算过程。

6.根据权利要求4所述的一种基于GPU的三维地震波场模拟计算系统，其中，对每个时间计算得到的地震波场数据进行保存，从而获得整个地震波场。