CN110162804A - 基于cpu加速的波场正演模拟优化方法 - Google Patents

Abstract

本发明公开了一种基于CPU加速的波场正演模拟优化方法，其包括以下步骤：S1)提出波场正演具体物理模型；S2)进行数值建模；在空间维度上使用高阶差分来模拟二阶微分，从而减小内存的使用；在时间维度上减少延拓步长；S3)进行静态分块及分配：根据线程数量k，将图像沿长边方向平均分成k个子图，然后将其按顺序标记后作为基础子图；S4)并行化及通信时间隐藏；S5)模型拆分和任务分解。本发明能够有效提高波场正演模拟的计算速度，并且能动态分块调配节点，从模型分块和通信时间隐藏两个方面优化，进行并行加速，充分利用多核计算资源，并有效的同步模拟实际物理传播过程。

Description

基于CPU加速的波场正演模拟优化方法

技术领域

本发明涉及计算机数值模拟领域，尤其涉及一种基于CPU加速的波场正演模拟优化方法。

背景技术

通过计算机模拟获得物理过程，已经成为现代物理学和工程学不可或缺的技术手段，由于很多物理现象空间或时间尺度过大或过小，有的条件过于复杂，实验成本过高，直接通过物理实验的方式来获得实验结果往往变得很不实际。先通过数值模拟的方式进行实验可以节约时间和物质成本，并且实验具有灵活性和可重复性。近年来，随着计算机集群计算能力的提高，数值模拟在天气预报、石油勘探、航天器设计等领域取得了丰硕的成果。

依据波动理论进行正演的方法是把所进行正演的地质模型空间首先进行网格化，这样网格上每个点就能获得对应的数值，然后用差分的方法或者有限元法等，把波动方程离散形式表达出来，然后根据时间步长进行地震波场的迭代，获得波场正演的模拟结果。波动理论正演方法因为能够准确描述地震波的动力学特征，即地震波的振幅、频率、相位属性的变化，因而可以作为一种研究波场的有效方式。所以以波动方程为理论基础的正演模拟，对认识地震波传播规律、进行地震属性研究、地震资料地质解释、储层评价等,均具有重要的意义和价值。以波动方程为基础的数值模拟包含了地震波的运动学特征和动力学特征，因此产生的波包含了种类繁多的信息，使我们在研究波的传播规律和复杂地层时有了更多的依据。简言之，地震波场数值模拟就是把所谓的地球物理问题进行分析、简化，进而形成一个能用数学形式表达的模型，并通过数值计算方式得到相应的地震相应。它是研究在地震条件下，地层构造和大量的岩石物性和岩性等信息与地震波包含的波动信息相对应的一门技术。地震波传播理论是地震勘探的基础，地震波在地下地质构造的传播中所包含的波动信息对探明地层内部结构是一种很重要的手段。

有限差分方法是最早被建立起来的数值模拟方法，该方法具有计算成本低，实现便捷的特点，在网格剖分足够密集的情况下，可以清晰的刻画出各种物理场的细节。由于其计算法则，便于并行化和各种处理器加速，至今仍然被广泛应用于波动力学、流体力学以及热力学领域的数值计算模拟工作中。由于在计算机集群实现过程中模型规模太大，超出了单个节点的处理能力，必须对模型按节点进行分割，导致沿着时间步长延拓的过程中，节点之间的边界数据需要通信。特别是采用显格式计算的时候，由于时间步长被控制得很小，加大了通信成本。

现有技术方案中存在波场正演模型在集群并行计算效率低的问题，此外还存在以下不足：

1、只对原始串行算法进行了表层的静态并行化。

波场正演算法比较复杂，涉及到物理数学和数值计算的知识，且存在过分割问题，这需要进行分割后处理。波场正演运算过程中涉及到非常多的递归运算，且物理波场模拟的数据量庞大，因此在模型分割过程中计算量很大，耗时长，所需计算资源较多，故在分割算法计算过程中主要对其进行并行化改造。如何进行算法的并行化改造是需要技巧的，需要进行合理的算法设计。然而现有算法仅仅进行了最简单的静态并行化改造，对算法的速度提升有限，相当于把模型分割进行并行计算，并没有充分考虑到各个部分之间的联系，没有充分利用计算资源，造成资源的闲置、浪费。因此需要根据算法运行中，各任何和计算资源的状态，设计动态并行算法，进一步提高算法并行效率。

2、在对过分割小模型区域合并时没有考虑算法的通信噪声抗干扰能力。

现有算法在进行过分割小模型区域合并时，基本采用对两个模型边缘区域之间数值的直接应用替换，按照时间序列进行计算。这种方式对无前一个模型对后一个模型影响较小的数值计算的处理有较好效果，但是对于波场正演的物理模型，则会受到严重的干扰，鲁棒性较差。主要原因是震源对四周的扩散模拟有较大幅度的影响，使得后面的结果会发生大幅偏离。需要改进过分割小区域合并算法，使之具有一定的抗噪声干扰能力。

3、对算法并行化后产生的额外的缝合步骤没有进行并行化。

分割并行化和其他基于数值计算的处理并行化不同，存在典型的＂缝合线＂问题，这是并行化后产生的额外步骤。然而现有算法只将次步骤串行处理。随着模型规模的提高、并行分块的数量增大，必然导致最后需要缝合子模块的数量增大，计算量也相应提高，成为整个算法效率提升的瓶颈。为此需要设计合适的处理方法，将缝合拼接部分也并行化，隐藏模块之间的通信时间，深层地提高算法效率。

发明内容

针对现有技术之不足，本发明提出了一种基于CPU加速的波场正演模拟优化方法，其包括以下步骤：

S1：提出波场正演具体物理模型；

对Navier方程进行简化，只保留压缩场分量从而得到声波方程，其偏微分方程表达形式为

将求解的区域进行网格剖分，用差分算子代替微分算子以便于计算机处理；(1)式的差分格式可表示为：

S2：进行数值建模；在空间维度上使用高阶差分来模拟二阶微分，从而减小内存的使用；在时间维度上减少延拓步长，根据傅里叶分析的原理得到显格式的稳定性条件为：

S3：进行静态分块及分配：根据线程数量k，将图像沿长边方向平均分成k个子图，然后将其按顺序标记后作为基础子图；

S4：并行化及通信时间隐藏；

在时间维度上，在计算时分配三段内存来分别存储上一时刻、当前时刻和下一时刻的波场，当一次正演完成后交换指针地址，从而通过三段内存交替使用来完成整个波场的不断推演；

在空间维度上，使用计算机集群做并行计算以进行有限差分正演模拟中空间维度上的差分。

此外，本发明还包括：模型拆分和任务分解；

将模型沿着空间的某一个方向进行拆分，把任务均匀的分配到各个节点中，每个节点只处理相应部分区域的数据；每个节点处理的区域包括D1和D2部分，D0区域为差分算法边界需要用到的相邻节点的数据；

计算D1区域，同时采用多线程技术将相邻节点的D0区域的数据传输给相邻节点的D2区域；

在D1区域计算完成并且D0区域的数据通信也完成时，对D2区域做有限差分计算。

本发明具有以下有益效果：

由于现有算法大部分只是对并行实现的简单应用，所以在最后拼接阶段基本采用单线程/进程串行处理的方式，没有充分利用并行优势。本发明将串行拼接改为层叠拼接，将待拼接子模块进行并行分配，提升了算法效率。

同时现有的波场正演并行分割算法均没有考虑对图像进行平均划分和线程分配这一点，每个分块小模型的运算时间也不相同。此问题会导致部分任务较轻的线程结束自己静态分配的任务后处于闲置状态，浪费计算资源。而任务重的线程则成为整个算法的瓶颈，算法的最终耗时依赖于最后一个运算结束的线程。本发明采用了通信隐藏的方法，充分利用每条线程的计算能力，达到负载的动态均衡。

附图说明

图1示例性地示出了进行波场正演模拟时的工作流程图；

图2示例性地示出了进行静态分块的示意图；

图3示例性地示出了正演过程中的内存使用示意图；

图4示例性地示出了拉普拉斯算子作用于差分方程的示意图；

图5示例性地示出了模型拆分和任务分解示意图；

图6示例性地示出了延迟隐藏算法的内存区域划分示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

图1为采用本发明方法优化后，在进行波场正演模拟时的工作流程图。需要说明的是，该流程图是示例性的，其用于理解本发明，并非构成对本发明优化方法的限制。

本发明的基于CPU加速的波场正演模拟优化方法包括以下步骤：

步骤S1：提出波场正演具体物理模型。

步骤1.1：对Navier方程进行简化，只保留压缩场分量从而得到声波方程，其偏微分方程表达形式为：

(1)式中P表示波场，c表示模型空间X、Y、Z处的介质速度，即该点声音在介质中的传播速度，我们通常把(1)式右边关于X、Y、Z的偏微分项称为拉普拉斯算子。

步骤1.2：将求解的区域进行网格剖分，用差分算子代替微分算子以便于计算机处理；(1)式的差分格式可表示为：

在(2)式中，选择二阶差分中的中间差分格式，用一个二阶的差分方程去近似一个二阶的微分方程，由此带来的问题是离散造成的截断误差，以空间项X方向为例，二阶差分导致二阶截断误差

步骤S2：数值建模。

为了使差分方程更好的接近微分方程，需要尽可能的减少截断误差，否则在数值模拟过程中空间网格过大，波长的群速度相速度不一致会导致频闪误差，这种误差会随时间的增加而向外传播，使最终获得的波场变得混乱不堪。而在空间上，如果空间采样率过大则会导致方程不稳定，波场数值随时间增加而变得异常大。

一个比较简单的解决办法就是减小时间步长和空间步长，对于时间来说减小时间步长意味着增加计算量，在计算机技术高速发展的今天，这是一个可以容忍的问题。然而，在保证模拟目标不变的情况下减小空间步长，则意味着增加网格数量，加大内存开销，而且空间有三个维度，这就会造成内存开销随采样间隔减小成三次方增长。在大规模数值模拟实验中，对于计算机集群来说往往变得不可能忍受。因此，使用高阶差分来模拟二阶微分，这样做会增大计算量，好处就是可以大幅减小内存的使用，这是一种时间换空间的做法。

假如使用5个点的差分格式来拟合一个二阶有限差分，则有

对于(4)式中的各项进行泰勒展开可得到关于P和P的偏导之间的关系：

将(5)式代入(4)式可得到线性方程组

等式(6)满足的充要条件是：

通过简单的线性方程组求解算法解(7)式，可求得各项权重系数的值，上式具有5阶的截断误差，即4阶的差分精度，故被称为4阶差分格式。

在时间方向上没有采用高阶差分格式，而是单纯的减小时间步长。

本技术方案的主要思想之一就是利用时间来换空间，如果时间上采用高阶差分则需要存储更多的时间快照，这与我们的初衷背道而驰。当然，也可以使用crank-nicolson格式增加时间步长，但是隐藏格式需要求解方程，同样会增加计算量。

为了保证稳定，在时间维度上直接采用减少延拓步长的办法，根据傅里叶分析的原理可得到显格式的稳定性条件为

通过比较可知本算法有两大优点，一是大部分模型分割确定区域可以同时计算，而不必等到分割后前一块模型计算完后再次计算接下来的分割模型；二是通过隐藏边界条件有一定的抑制噪声的作用，增加算法鲁棒性。

步骤S3：静态分块及分配。

根据线程数量k，将图像沿长边方向平均分成k个子图并且按顺序进行标记，构成子图集合，其称为基础子图，具体如图2所示。

步骤S4：并行化及通信时间隐藏。

有限差分算子是根据上一时刻的波场和当前时刻的波场推演出下一时刻的波场。把(1)式左边关于时间项的差分展开，并把空间差分项写成拉普拉斯算子，对其进行整理，可得：

从公式(9)中可以看出波场延拓的基本思路是根据上一个时刻的波场，和当前时刻的波场以及速度模型来推演下一时刻的波场。在计算的时候分配三段内存来分别存储上一时刻、当前时刻、下一时刻的波场，当一次正演完成后交换指针地址。有限差分正演属于流计算，当完成一次正演后，上一时刻的波场就不会再被使用。因此，把当前时刻波场的地址传递给上一时刻，把下一时刻的地址传递给当前时刻，把上一时刻的地址传递给下一时刻，在进行下一次推演的时候直接把该段内存覆盖掉。这样只需要三段内存交替使用便可以完成整个波场的不断推演。如果用M1、M2、M3分别代表三段内存，三段内存在正演过程中交替使用的情景可由如图3所示。

有限差分正演模拟除了时间维度的差分还有空间维度上的差分，从(9)式中的第一项可以看出，还需要对当前时刻的波场世家拉普拉斯算子这样才能驱动波场的扩散。根据拉普拉斯算子可知，计算空间上每个节点的拉普拉斯算子需要用到上下左右前后6个方向的领域空间的波场值。

如果采用的是二阶有限差分，则如图4所示，要计算当前结点的空间差分需要用到领域内的7个点的值。图4中a)为当前时刻的波场；b)为拉普拉斯算子的计算结果。如果采用16阶有限差分，使用的数据沿着6个方向延伸，总过用到49个点。少量的数据在单机上计算是没有问题的，找到对应的点取值就可以了。然而，在实际应用中经常会做一些很大规模的数据处理，或者有些应用场景下要求非常高的模拟精度，被迫减少网格的间距，这也会导致实际处理数据的规模增大。这个时候，对于单机来讲无论是内存还是计算的效率都是无法忍受的，只能使用计算机集群做并行计算。

步骤S5：模型拆分和任务分解。

波场是沿着时间步长进行正演，不同的时刻数据之间存在因果关系，无法做并行处理。因此，只能将模型沿着空间的某一个方向进行拆分，把任务均匀的分配到各个节点中，每个节点只处理部分区域的数据，如图5所示。

在模拟的过程中，每一个时刻在计算这个位置的拉普拉斯算子的时候都要从相邻的其他节点上获取数据，一个比较简单的做法就是，先让各个节点做计算。当计算到边界领域的时候就和相邻的节点通信，把数据拷贝过来再继续计算。如果差分使用的阶数越高，从相邻的节点上获取的数据就会越多。这样在通信的过程中，计算机处于等待状态，就会降低计算的效率。

从图6中可以看出，每个节点处理的区域为D1和D2部分；D0区域是差分算法边界需要用到的相邻节点的数据，它对应着相邻节点的D2部分，也就是需要计算的边界领域。这一部分的通信就产生通信延迟，而且每一次时间步长的延拓走需要做通信。为了保证显格式的稳定性，时间步长往往取得很小，导致正演次数变得很大，同时每一次延拓都伴随着通信，那么要完成整个计算，通信时间就会被累加得特别大，以至于影响计算效率。为了解决这个问题，我们采取两步来隐藏通信延迟：

(1)计算D1区域的部分，因为D1部分区域最多延拓到D2区域，不需要通信直接就可以计算。同时，使用多线程技术，将相邻节点的D0区域的数据传输给相邻节点的D2部分。由于是多线程，所以计算和传输是可以同时进行的。

(2)当D1部分计算完，D0区域通信也已经完成，这个时候在对D2区域做有限差分计算。这种方法的好处在于通信的同时也在做计算，通过D1区域的计算把通信时间隐藏掉这样就达到提高计算效率的效果。特别是在GPU或者FPGA平台上由于节点的计算效率太高，通信延迟问题就显得尤为突出。

下表为采用本发明方法后的通信隐藏优化结果表：

本发明能够有效提高波场正演算法的计算速度，并且能动态分块调配节点，从模型分块和通信时间隐藏两个方面优化，进行并行加速，充分利用多核计算资源，并有效的同步模拟实际物理传播过程。此外，本发明的方法可扩展到数值模拟的很多相关领域。

需要注意的是，上述具体实施例是示例性的，本领域技术人员可以在本发明公开内容的启发下想出各种解决方案，而这些解决方案也都属于本发明的公开范围并落入本发明的保护范围之内。本领域技术人员应该明白，本发明说明书及其附图均为说明性而并非构成对权利要求的限制。本发明的保护范围由权利要求及其等同物限定。

Claims (2)

1.一种基于CPU加速的波场正演模拟优化方法，其特征在于，包括以下步骤：

S1：提出波场正演具体物理模型；

对Navier方程进行简化，只保留压缩场分量从而得到声波方程，其偏微分方程表达形式为：

将求解的区域进行网格剖分，用差分算子代替微分算子以便于计算机处理；(1)式的差分格式可表示为：

S2：进行数值建模；在空间维度上使用高阶差分来模拟二阶微分，从而减小内存的使用；在时间维度上减少延拓步长，根据傅里叶分析的原理得到显格式的稳定性条件为：

S3：进行静态分块及分配：根据线程数量k，将图像沿长边方向平均分成k个子图，然后将其按顺序标记后作为基础子图；

S4：并行化及通信时间隐藏；

在时间维度上，在计算时分配三段内存来分别存储上一时刻、当前时刻和下一时刻的波场，当一次正演完成后交换指针地址，从而通过三段内存交替使用来完成整个波场的不断推演；

在空间维度上，使用计算机集群做并行计算以进行有限差分正演模拟中空间维度上的差分。

2.如权利要求1所述的方法，其特征在于，还包括：S5：模型拆分和任务分解；

将模型沿着空间的某一个方向进行拆分，把任务均匀的分配到各个节点中，每个节点只处理相应部分区域的数据；每个节点处理的区域包括D1和D2部分，D0区域为差分算法边界需要用到的相邻节点的数据；

计算D1区域，同时采用多线程技术将相邻节点的D0区域的数据传输给相邻节点的D2区域；

在D1区域计算完成并且D0区域的数据通信也完成时，对D2区域做有限差分计算。