一种基于异步流的地震数据多GPU快速正演计算方法

摘要

本发明公开了一种基于异步流的地震数据多GPU快速正演计算方法。所述方法将地震数据按照慢维域进行分割，分割的子域个数与GPU数目相同，将每个子域的数据分配给一个GPU进行计算；实现各GPU两两之间双向实时通信；采用多异步流技术控制计算流程。通过所述方法既可以保证高效率实施大数据体正演计算，又可以使双GPU计算效率与单GPU相比达到2:1。

说明书

技术领域

本发明属于地球物理勘探技术领域，具体涉及一种基于异步流的地震数据多GPU快速正演计算方法。

背景技术

波动方程正演模拟在地震数据的反演解释中发挥着重要意义。对于三维地震资料的处理和解释需要有效的正演模型来验证。同时波动方程正演模拟是逆时偏移成像技术和全波形反演技术的基础。最常用的波动方程正演模拟方法是基于全波方程的，如有限差分、伪谱法和有限元方法，可以模拟地震波在地下的传播特征，适用复杂介质模型，但是计算量和内存开销很大，尤其是三维波动方程模拟。有限差分是通过对时间和空间导数的网格差分离散近似来实现波动方程的模拟。伪谱法是通过正、反傅里叶变换来实现空间导数的精确求解，计算量相比有限差分大。有限差分方法在地震波正演模拟和偏移领域已经得到了长足的发展和广泛的应用，它的优点是差分格式稳定性强，能适应速度任意变化的介质，并且有限差分方法适合GPU并行加速计算。但是有限差分方法存在数值频散问题，这直接影响到有限差分方法在波动方程中的应用。数值频散问题是由于对时间和空间偏导数的网格化离散造成的，导致地震波模拟的数值相速度不等于实际介质的真实相速度，使波场模拟的精度降低。

为了改善或压制网格数值频散问题，一般采用有以下几种方法：一是采用低阶差分格式，给定一个很小的模拟时间步长和空间步长，这极大的加大了计算量，尤其三维情况，所以该方法一般不适用三维数据；二是采用高阶有限差分格式，一般时间导数采用二阶差分格式，空间导数采用高阶有限差分格式；三是采用更高精度的算法，比如以交错网格替代常规的均匀网格算法。Dablain和Chen提出了时间方向导数的四阶有限差分格式，提高了波场模拟的精度。总之，一般高精度的数值模拟对应更大的计算量。

综上所述，随着数据量和计算量的增加，计算效率是制约正演工业界广泛应用的重要因素。基于PC集群的并行计算越来越不能满足计算需求，并且成本高，而基于GPU计算的众核并行计算技术因其强大的计算能力并且成本较低、占用空间少、方便等特点在地震资料正演以及偏移成像中被广泛应用。Micikevicius（2008）提出了GPU平台的三维有限差分正演算法，李博等（2010）对比分析了基于CPU/GPU平台的逆时偏移实现策略，刘红伟等（2010）讨论了地震叠前逆时偏移高阶有限差分算法及GPU实现，刘国峰等（2013）给出了基于CPGPU系统逆时偏移实现的快速计算算法。但是，GPU相对于普通的PC集群计算，内存较小，对于大规模三维数据来说，单个GPU的内存不足以完成计算。因此，迫切需要研究提出利用多个GPU并行实施地震数据正演模拟计算的方法。

发明内容

为了克服在大规模地震数据推演计算需求中存在的PC集群计算能力不足，单个GPU内存不足等问题，本发明提供一种基于异步流的地震数据多GPU快速正演计算方法，所述方法既可以保证高效率实施大数据体正演计算，又可以使双GPU计算效率与单GPU相比达到2:1。

为实现上述目标，本发明采用以下技术方案：

1、计算原理

波动方程正演计算的原理可以描述为：以子波模拟实际震源，沿时间实现波动方程的延拓计算，最后抽取地表的波场值得到地震正演记录。其本质是求解波动方程，进行波场的延拓。声波方程可以写作如下形式：

经2N阶有限差分离散后得到：

其中fdcoe[l]为有限差分计算的差分系数，s(it)为震源子波。

2、方法步骤

基于多GPU计算的快速正演算法核心是基于异步流的慢维域的分解，本发明的慢维对应深度。正演计算时，沿着深度将模型空间等分为n份（n对应GPU的个数），每个GPU计算一个子模型域的网格波场值，并且在计算过程中实现子模型域间的实时通讯。每个GPU是通过将对应模型网格点波场的值分配到GPU的不同线程进行计算。

多卡计算的核心问题之一是两两GPU卡间直接的通讯，这里以双卡并行计算为例，图1展示了两GPU实现正演模拟的模型分解示意图。具体过程如下：整个模型区域沿着深度Z方向一分为二（S0、S1），分别分配给GPU0、GPU1计算；为了保证两个GPU之间的实时通讯，同时设置交换区域（transfer）及对应的接收区域（图中0-1与1-1），各GPU接收区域0-1（1-1）用于接收对方GPU的transfer区域1-0（0-0）计算更新的波场，以实现两卡之间的实时通信，保证波场的连续传播。

我们为了进一步提高计算效率，提出多异步流的计算策略，这是一种任务级别的并行，旨在提高并行度，最大化占用带宽，充分发挥并行优势。具体步骤如下：

1.首先图1中所示的0-0和1-0所指示的交换区域进行波场递推，与此同时GPU0、GPU1分别计算模型主体中的S0、S1所示区域，四者并发执行，这四块区域的计算采用四个不同的异步流实现并发计算，如图2所示的流Ct_S0、Tf_00、Tf_10、Ct_S1。

2. 上述步骤计算过程中，发送交换区域拷贝的任务进程，0-1所示接收区域等待1-0区域更新的波场拷入，即将1-0区域计算更新后的波场传递给0-1，类似的，1-1所示接收区域等待0-0区域更新的波场拷入，即将0-0计算更新的波场传递给1-1，这种存在等待，具有先后顺序的计算，可以采用流技术控制，具体指通过同一个流来控制这些顺序任务，使之位于一个队列，如图2中的流Tf_00与Tf_10，这种GPU0与GPU1间的直接通信，就是我们所说的P2P技术，如图2所示的虚线箭头。

3、发明要点

一种基于异步流的地震数据多GPU快速正演计算方法，所述方法将地震数据按照慢维域进行分割，分割的子域个数与GPU数目相同，将每个子域的数据分配给一个GPU进行计算；实现各GPU两两之间双向实时通信；采用多异步流技术控制计算流程。

一种基于异步流的地震数据多GPU快速正演计算方法，包括以下步骤：

1）基于异步流的慢维域的分解；

2）实现两两GPU卡间直接的通讯；

3）多异步流的计算策略。

所述步骤1）包括以下子过程：

a）沿着慢维域将模型空间等分为n份，n等于GPU的个数；

b）将每一个子模型域的网格波场值分配给一个GPU实施计算；

c）将对应子模型域的网格点波场值分配到每个GPU的不同线程进行并行计算。

所述步骤2）的实现过程是：在每个子模型域同时设置交换区域及对应的接收区域，用于保证两个GPU之间的实时通讯；各GPU的接收区域用于接收对方GPU的交换区域计算更新的波场，以实现两卡之间的实时通信，保证波场的连续传播。

所述步骤3）包括以下子过程：

a）N个交换区域进行波场递推，与此同时N个GPU分别计算模型主体中的对应区域，2N个进程并发执行，这2N个区域的计算采用2N个不同的异步流实现并发计算；

b）计算过程中，发送交换区域拷贝的任务进程，即各GPU的接收区域，等待对方GPU的交换区域更新的波场拷入，这种具有先后顺序的计算，采用流技术控制，具体指通过同一个流来控制这些顺序任务，使之位于一个队列，这种GPU之间的直接通信，就是P2P技术。

参见附图2，单炮具体正演模拟的主要步骤如下：

1）设置时间初值it=0；

2）添加震源子波；

3）炮点波场延拓；

4）交换不同时刻波场；

5）增加时间值it=it+1；

6）判别时间是否达到阈值，若达到，则抽取正演记录，否则，返回步骤2）。

本发明的优点和有益效果为：和传统的PC集群计算方法相比，该方法计算效率高，成本较低。与单GPU计算方法相比，多卡计算扩大了显存，以实现精细的大规模数据计算，同时提高了计算的并行度，进一步提高了计算效率。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1为本发明所述的双卡正演计算模型域分解示意图。

图中，(a)是整个模型空间；(b)是GPU0上的计算子模型空间；(c)是GPU1上的计算子模型空间。

图2为本发明所述的多卡正演的流程图。

图3为本发明所述的单卡与双卡正演耗时结果对比图。

图4为本发明所述的双卡正演的加速比示意图。

图5为本发明所述的正演模拟记录。

图中，左图是单卡结果，右图是双卡结果。

图6为本发明所述的随机抽取的单道记录对比图。

图中，左图是单卡结果，中间是双卡结果，右图是二者误差。

具体实施方式

实施例

参见附图3-6，根据上述步骤与方法，我们采用数据进行测试。测试的GPU为Tesla K80，CPU主频为2.6GHz。测试结果如附图3-6所示。其中横坐标表示正演模型尺寸大小，附图3的纵坐标表示计算耗时。粗线表示双卡GPU耗时，细线表示单卡GPU耗时。附图4表示加速比的对比图，以单卡时间为参考，得到单卡和双卡计算时间的比值，从图中可以看出，双卡的加速比超过了2:1，证明了我们异步流方法对效率的提高效果。最后对比单卡与双卡的计算模拟结果。附图5是单卡与双卡的正演模拟结果；附图6是随机抽取的单道结果对比，对结果进行作差处理，各时刻差值均为0，表明双卡模拟结果与单卡结果完全一致，没有误差。

最后应说明的是：显然，上述实施例仅仅是为清楚地说明本发明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明的保护范围之中。