CN105467443B - 一种三维各向异性弹性波数值模拟方法及系统 - Google Patents

Abstract

本发明涉及一种三维各向异性弹性波数值模拟方法及系统，其中方法包括：步骤1：建立介质模型，对介质模型进行网格离散得到多个网格点；步骤2：计算震源函数，根据震源函数计算每个网格点上的压力值；步骤3：将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；步骤4：根据波场值确定每个网格点的计算区域，进行分区并对分区边界数据进行数据交换，完成弹性波数值模拟。本发明实现了利用GPU加速复杂介质下的三维弹性波数值模拟，提出了利用GPU Direct技术加速数据传输的实现方案，避免了大量的CPU到GPU，GPU到CPU的数据拷贝，实现了优化通信瓶颈的问题。

Description

一种三维各向异性弹性波数值模拟方法及系统

技术领域

本发明涉及一种三维各向异性弹性波数值模拟方法及系统，属于地球物理勘探领域。

背景技术

地震波数值模拟以地震波在地下介质中的传播理论为基础，在勘探地震学与天然地震学中得到了广泛的应用。目前常规的三维声波方程、弹性波各向同性数值模拟已经广泛的应用于数值模拟、成像、以及反演等各个地球物理领域中。但是对于如何有效的实现大尺度采集方式(例如宽方位采集)和复杂各向异性介质(例如水平横向各向同性HTI或者是正交各向异性)的三维弹性波波场数值模拟研究仍然存在很大的挑战，在实际应用中未得到广泛的使用。此外常规的基于CPU的大规模三维声波、弹性波数值模拟，通常需要大量专用集群计算资源。无论从硬件成本还是计算能耗上来说其成本都十分昂贵。因此快速的提高计算性能、显著的降低计算成本对于实现三维弹性波各向异性数值模拟的实际应用具有重要意义。

近十年中，利用图形处理单元(GPU)进行计算密集型应用的加速实现已经得到了突飞猛进式的发展。图形处理单元(GPU)由于其具有高速的内存带宽，相较于CPU至少高出两个数量级的计算处理核心，更适合并行计算的单指令多数据(SIMD)计算模式，以及更低的能耗成本，正广泛的应用到计算科学的相关领域。对于勘探地球物理领域，对于使用图形处理单元(GPU) 的兴趣也在显著增强，越来越多的研究已经将GPU用于加速地震处理中的核心算法，例如地震数值模拟、地震成像、地震高精度反演等。

发明内容

经过研究发现，现有技术及其存在以下问题：

将时间域有限差分方法应用到GPU设备上实现复杂介质波传播算法上的研究很少，尤其对于如何处理存在大量多节点数据交换的大规模三维问题研究就更少(Heinrichet al，2014)。

目前存在一些GPU地震波传播模拟的实现方案，对于三维GPU集群的实现方案我们以(龙桂华等)的方案为例进行描述。该方案利用区域分解技术将单个GPU上不能计算的地质体沿Z轴方向进行粗粒度分解，采用消息传递接口(MPI)交换边界数据，从而运用MPI+CUDA的方式实现了大尺度三维地震波场的数值模拟。但是该方法存在很大的一个问题是，采用GPU集群计算的加速比相较于单个GPU与CPU来说显著下降，造成这一结果的原因是GPU计算的大部分耗时消耗在GPU集群节点间的GPU到CPU以及CPU到GPU设备之间的数据拷贝上。

本发明所要解决的技术问题是，针对现有技术的不足，提供一种将GPU Direct技术与有限差分数值模拟方法结合得到的基于GPU Direct优化通信的三维各向异性弹性波数值模拟方法及系统，为大规模数值模拟提供有力保证。

本发明解决上述技术问题的技术方案如下：一种三维各向异性弹性波数值模拟方法，具体包括以下步骤：

步骤1：建立介质模型，对介质模型进行网格离散得到多个网格点；

步骤2：计算震源函数，根据震源函数计算每个网格点上的压力值；

步骤3：将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

步骤4：根据波场值确定每个网格点的计算区域，进行分区并对分区边界数据进行数据交换，完成弹性波数值模拟。

本发明的有益效果是：从应力应变方程出发，实现了利用图形处理单元 (GPU)加速复杂介质下的三维弹性波数值模拟，并针对三维问题引入GPU设备所面对的由于区域分解所造成的多节点间、节点内通信瓶颈问题进行了深入的研究和分析，提出了利用GPUDirect技术加速数据传输的实现方案，避免了大量的CPU到GPU，GPU到CPU的数据拷贝，实现了优化通信瓶颈的问题。通过GPU图形加速设备并引入GPU Direct通信优化策略，可以显著的提升整体计算性能，可以用更低的硬件成本和更少的时间实现的三维各向异性弹性数值模拟，为各种依赖于波动方程正演模拟的算法如逆时偏移，全波形反演的应用提供有力的保证。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述步骤4具体包括以下步骤：

步骤4.1：根据所有波场值确定每个网格点的计算区域，并采用吸收边界的方式确定计算区域边界，在计算区域内模拟地下介质中波的传播；

步骤4.2：对所有计算区域进行分区，对相邻分区的边界数据进行数据交换，得到模拟弹性波数据，完成弹性波数值模拟。

进一步，所述步骤4.2中边界数据的采用GPU-Direct技术进行数据交换。

采用上述进一步方案的有益效果是，利用GPU加速三维弹性波数值模拟，显著提高计算速度，相较于常规CPU集群实现，可以达到20-30倍的速度提升。

进一步，所述步骤4.2中的数据交换具体包括：

沿介质模型中计算区域边界数据变化最慢的轴向进行区域分解，得到多个分区，将所有分区分配到多个GPU中，每个分区独立在一个GPU上执行计算；每两个相邻分区的边界数据再进行交换。

采用上述进一步方案的有益效果是，利用GPU Direct技术加速节点间、节点内数据传输的实现方案，我们的方案因其避免了大量的CPU到GPU，GPU 到CPU的数据拷贝，从而实现了优化通信瓶颈的问题。

进一步，所述步骤1具体包括以下步骤：

步骤1.1：根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立介质模型；

步骤1.2：采用形状规则的三维网格对介质模型进行网格离散，得到网格点。

采用上述进一步方案的有益效果是，将介质模型离散成网格点，网格点作为波场值的载体，由三维结构的多个网格点构成波场，作为波场传播的初始条件。

进一步，所述震源函数在空间上采用高斯函数，在时间上采用Ricker 子波，所述震源函数的具体公式为：

s(x，y，z，t)＝g(x，y，z)·f(t) 公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²) 公式(2)

公式(3)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀，y₀，z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

进一步，所述步骤3具体包括以下步骤：

步骤3.1：将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长满足该数值格式的稳定性条件；

步骤3.2：采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值。

本发明解决上述技术问题的技术方案如下：一种三维各向异性弹性波数值模拟系统，包括建模模块、震源模块、传播模块和数据交换模块；

所述建模模块用于建立介质模型，对介质模型进行网格离散得到多个网格点；

所述震源模块用于计算震源函数，根据震源函数计算每个网格点上的压力值；

所述传播模块用于将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

所述数据交换模块用于根据波场值确定每个网格点的计算区域，进行分区并对分区边界数据进行数据交换，完成弹性波数值模拟。

本发明的有益效果是：从应力应变方程出发，实现了利用图形处理单元 (GPU)加速复杂介质下的三维弹性波数值模拟，并针对三维问题引入GPU设备所面对的由于区域分解所造成的多节点间、节点内通信瓶颈问题进行了深入的研究和分析，提出了利用GPUDirect技术加速数据传输的实现方案，避免了大量的CPU到GPU，GPU到CPU的数据拷贝，实现了优化通信瓶颈的问题。通过GPU图形加速设备并引入GPU Direct通信优化策略，可以显著的提升整体计算性能，可以用更低的硬件成本和更少的时间实现的三维各向异性弹性数值模拟，为各种依赖于波动方程正演模拟的算法如逆时偏移，全波形反演的应用提供有力的保证。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述数据交换模块包括区域确定模块和分区交换模块；

所述区域确定模块用于根据所有波场值确定每个网格点的计算区域，并采用吸收边界的方式确定计算区域边界，在计算区域内模拟地下介质中波的传播；

所述分区交换模块用于对所有计算区域进行分区，对相邻分区的边界数据进行数据交换，得到模拟弹性波数据，完成弹性波数值模拟。

进一步，所述分区交换模块中边界数据的采用GPU-Direct技术进行数据交换。

采用上述进一步方案的有益效果是，利用GPU加速三维弹性波数值模拟，显著提高计算速度，相较于常规CPU集群实现，可以达到20～30倍的速度提升，利用GPU Direct技术加速节点间、节点内数据传输的实现方案，我们的方案因其避免了大量的CPU到GPU，GPU到CPU的数据拷贝，从而实现了优化通信瓶颈的问题。

附图说明

图1为本发明实施例1所述的一种三维各向异性弹性波数值模拟方法流程图；

图2为本发明实施例中利用八阶精度有限差分近似中心点的一阶导数需要25个网格点图；

图3为本发明实施例中三维各向异性弹性波数值模拟在单GPU设备上的实施流程图；

图4为本发明实施例中数据体区域分解及通信情况概述图；

图5a为现有技术中节点内采用MPI数据传输方式示意图；

图5b为本发明实施例中采用GPU Direct P2P数据传输方式示意图；

图6a为现有技术中节点间采用MPI数据传输示意图；

图6b为本发明实施例中采用GPU Direct RDAM方式进行数据传输示意图；

图7为本发明实施例1所述的一种三维各向异性弹性波数值模拟系统结构框图；

图8a-8c为本发明实施例中用于测试的三种具有不同TI对称性的弹性介质刚度系数矩阵示意图；

图9a-9c为本发明实施例中三种具有不同TI对称性的弹性介质的三维脉冲响应示意图；

图10为本发明实施例中节点内常规MPI通信方式与GPU Direct通信方式的测试分析结果对比图；

图11为本发明实施例中节点间采用常规MPI通信方式与采用GPU Direct RDMA通信方式的测试结果对比图。

附图中，各标号所代表的部件列表如下：

1、建模模块，2、震源模块，3、传播模块，4、数据交换模块。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

如图1所示，本发明实施例1所述的一种三维各向异性弹性波数值模拟方法，具体包括以下步骤：

步骤1：建立介质模型，对介质模型进行网格离散得到多个网格点；

步骤2：计算震源函数，根据震源函数计算每个网格点上的压力值；

步骤3：将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

步骤4：根据波场值确定每个网格点的计算区域，进行分区并对分区边界数据进行数据交换，完成弹性波数值模拟。

实施例2中，在实施例1的基础上，所述步骤4具体包括以下步骤：

步骤4.1：根据所有波场值确定每个网格点的计算区域，并采用吸收边界的方式确定计算区域边界，在计算区域内模拟地下介质中波的传播；

步骤4.2：对所有计算区域进行分区，对相邻分区的边界数据进行数据交换，得到模拟弹性波数据，完成弹性波数值模拟。

实施例3中，在实施例1或2的基础上，所述步骤4.2中边界数据的采用GPU-Direct技术进行数据交换。

实施例4中，在实施例1-3任一实施例的基础上，所述步骤4.2中的数据交换具体包括：

沿介质模型中计算区域边界数据变化最慢的轴向进行区域分解，得到多个分区，将所有分区分配到多个GPU中，每个分区独立在一个GPU上执行计算；每两个相邻分区的边界数据再进行交换。

此实施例中，进一步，对每一时刻的各个计算节点的计算边界数据采用 GPU-Direct技术进行数据交换具体实现如下：

当处理求解大规模三维模型时，有限的全局内存使得单个GPU设备无法存储整个模型网格。为此，我们的算法必须将CUDA代码扩展到运行在多GPU 设备，多节点的异构体系结构中，以充分利用GPU计算能力进行高性能计算。基于这个原因，我们采用一个区域分解方案，我们的区域分解方案如图3所示，沿着计算网格变化最慢的Y轴方向进行区域分解，从而将整个计算区域分解为N/M份，N为Y方向上的网格数目，M为当前计算集群的GPU数目。每个分区独立在一个GPU设备执行计算步骤1-3。

因为有限差分方法(如图2所示，以8阶精度为例)需要y方向上的向前和向后四个点，因此在每个子区域中进行有限差分操作的GPU设备需要获取其他设备上与其计算区域相邻点的值，因此相邻子区域的边界数据必须在每个时间步上进行交换，其数据传输路径如图3所示，即分为节点内和节点间两种类型，如图4所示，数据体区域分解及通信情况概述，图中术语说明， PCI-E，是PCI-Express的简称，指计算机总线和接口标准。InfiniBand，是一种专门用于服务器端网络连接的多并发网络架构，GPU Direct P2P，GPU Direct的节点内实现方式(称为GPU直连点对点技术)，GPU Direct RDMA， RDMA技术全称为远程直接数据存取，GPU Direct RDMA指GPU Direct技术的节点间实现方式(成为GPU远程直连技术)。节点内利用GPU Direct P2P 通信，节点间利用GPU Direct RDAM通信。通常，这种不同GPU设备间的通信由于数据需要进行GPU到CPU以及CPU到GPU的数据拷贝，整个通信过程十分耗时，这是常规GPU集群计算的一个主要瓶颈。

对此，本发明利用GPU Direct优化这些数据通信，其具体实施如下：

a)保证CUDA版本在6.0之上，并且显卡计算能力在2.5之上，以保证对GPU Direct的支持。

b)采用CUDA-Aware MPI，CUDA-Aware MPI为支持GPU Direct特性的针对GPU计算优化的MPI实现，现阶段支持CUDA-Aware MPI的MPI实现有， MAVPCH21.9以上版本或者OpenMPI 1.7以上版本，我们以MAVPICH2为例。

c)cudaSetDevice函数需要在MPI_Init之前调用，并利用MPI环境变量实现对GPU的指定，MVAPICH2：MV2_COMM_WORLD_LOCAL_RANK。

d)MV2_USE_CUDA需要在执行mpirun时添加。

e)对于节点内的数据交换采用GPU Direct P2P方式(GPU直连点对点)，在完成上述配置后可以直接利用cudaMemcoycopy进行GPU到GPU的数据传输，其具体描述参见图5。

f)对于跨节点的数据交换采用GPU Direct RDMA方式(RDMA技术全称远程直接数据存取，就是为了解决网络传输中服务器端数据处理的延迟而产生的，利用Infiniband网络直接进行GPU直连称之为GPU Direct RDMA)，在完成上述配置后可以直接采用MPI_Isend与MPI_Irecv发送和接受设备内存，而非像以往一样需要先将设备内存做显式的cudaMemcpy到主机内存再通过MPI_Isend与MPI_Irecv发送和接收，再通过显式的cudaMemcpy到设备端，进行通信，两种通信区别可见图6。

实施例5中，在实施例1-4任一实施例的基础上，所述步骤1具体包括以下步骤：

步骤1.1：根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立介质模型；

步骤1.2：采用形状规则的三维网格对介质模型进行网格离散，得到网格点。

实施例6中，在实施例1-5任一实施例的基础上，所述震源函数在空间上采用高斯函数，在时间上采用Ricker子波，所述震源函数的具体公式为：

s(x，y，z，t)＝g(x，y，z)·f(t) 公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²) 公式(2)

公式(3)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀，y₀，z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

实施例6中，在实施例1-5任一实施例的基础上，所述步骤3具体包括以下步骤：

步骤3.1：将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长满足该数值格式的稳定性条件；

步骤3.2：采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值。

此实施例中，进一步，所述步骤3.2中采用区域分解的方式利用GPU在各个计算节点对三维各向异性弹性波方程进有限差分计算，得到每一刻的波场值的具体方法如下：

首先，给出三维各向异性介质下的弹性波波动方程的数学描述，三维各向异性弹性波动方程其形式为：

公式(4)

式中，u为波场函数；σ为柯西应力张量；F_i是单位体积的体力项，其可以被等效为应力源，ρ是介质物质密度，表示二阶时间偏导。对于应力张量利用线性弹性理论可以得到：

σ_ij＝c_ijklε_kl， 公式(5)

式中ε_kl为线性应变张量，c_ijkl为四阶张量的弹性刚度矩阵，该参数需要根据介质的实际情况给出；而对于ε_kl我们有

公式(6)

其中ε_kl代表线性应变张量，代表对kth维求空间导数，u_l表示l th波场位移。根据此式我们就可以得行弹性波动方程的数值模拟，利用GPU进行的弹性波数值模拟的具体计算步骤如下：

3.2.1我们定义一个N_x×N_y×N_z的计算网格，并且指定N_t个时间步长，这样在空间和时间的每一个网格点就可以利用一个四元组进行表示：即， [x，y，z|t]＝[pΔx，qΔy，rΔz|nΔt]。因此在指定点的连续的波场位移记录就可以利用离散的网格点进行表示：

u_i|_x，y，z|t≈u_i ^p，r，q|n 公式(7)

我们近似一阶偏导通过一个如下的M阶精度的的中心差商算子D_x[·]，可以得到：

公式(8)

式中W_α是一个权重系数多项式，即有限差分系数，M为有限差分计算阶数；

空间差分算子D_y[·]和D_z[·]与D_x[·]相似分别表示y和z方向的偏导数。将所有这三者差分算子代入方程(6)，用以求解本构方程中的偏导数。

3.2.2对于时间离散，我们使用一个标准的空间二阶精度近似来计算方程(6)中的

公式(9)

将上述的差分算子代入方程(4)并且通过重新排列这些离散项，我们就可以通过时间步增的方式进行未知波场求解。例如对于我们需要给定其当前和之前的波场值 以及根据有限差分计算空间计算模版所需要x-，y-，z-方向相邻点。图3表述了空间有限差分计算模版的在当前时间步长关于点所需要的计算网格点。

如图7所示，为本发明实施例1所述的一种三维各向异性弹性波数值模拟系统，包括建模模块1、震源模块2、传播模块3和数据交换模块4；

所述建模模块1用于建立介质模型，对介质模型进行网格离散得到多个网格点；

所述震源模块2用于计算震源函数，根据震源函数计算每个网格点上的压力值；

所述传播模块3用于将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

所述数据交换模块4用于根据波场值确定每个网格点的计算区域，进行分区并对分区边界数据进行数据交换，完成弹性波数值模拟。

实施例2中，在实施例1的基础上，所述数据交换模块包括区域确定模块和分区交换模块；

所述区域确定模块用于根据所有波场值确定每个网格点的计算区域，并采用吸收边界的方式确定计算区域边界，在计算区域内模拟地下介质中波的传播；

所述分区交换模块用于对所有计算区域进行分区，对相邻分区的边界数据进行数据交换，得到模拟弹性波数据，完成弹性波数值模拟。

实施例3中，在实施例2的基础上，所述分区交换模块中边界数据的采用GPU-Direct技术进行数据交换。

本发明的基础是地震波的传播理论及其高性能GPU数值模拟方法，本发明基于GPUDirect通信技术实现了高性能的三维各向异性介质弹性波数值模拟方法，在具体示例中实施步骤分别为：

1.介质模型建立：

以不同TI对称性(各向同性、VTI、HTI)的三维弹性介质作为测试模型，其弹性参数如下：相同的各向同性参数，v_p＝2.0km/s，以及ρ＝2000kg/m³，但使用不同的Thomsen各向异性参数。令VTI介质， [ε₁，δ₁，γ₁]＝[0.2，-0.1，0.2]，HTI介质，[ε₂，δ₂，γ₂]＝[0.2，-0.1，0.2]。这些参数都通过相应的变换公式(Thomsen，1986)转换进刚度矩阵中。图8a-c分别表示用于脉冲响应的ISO、VTI、HTI介质6×6的刚度矩阵C_ij，用不同的颜色块来表示不同的值。

2.模型离散化：

对设计的介质模型进行网格离散，测试地震数据维度大小为 N_x×N_y×N_z＝200³，网格间距为Δx＝Δy＝Δz＝0.005km。

3.震源函数：

震源函数在空间上采用高斯函数，时间上采用Ricker子波，形式如下：

s(x，y，z，t)＝g(x，y，z)·f(t) 公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²) 公式(2)

公式(3)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀，y₀，z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

4.弹性波动方程的数值格式：

采用8阶精度进行有限差分数值模拟，则对应的我们利用公式(4)-(6)进行每一步的波场计算。

5.区域分解及GPU Direct通信：

采用前述的区域分解方案，在Y方向上进行区域分解，节点内采用GPU Direct P2P进行通信，节点间采用GPU Direct RDMA进行通信。测试环境为使用本机是三台双路GPU服务节点作为测试平台。每个节点配置有双路10 核Xeon E5-2680V2 2.80GHz处理器，128GBDDR3内存。GPU设备为NVIDIA Tesla K40c(Kepler)，Infiniband适配卡为MellanoxConnectX-3 IB QDR MT26428，所有节点运行在RHEL 6.3，Infiniband的驱动版本为 OFED-2.3-1.0.1。所有测试基于MVAPICH2-2.1a发行版本的MPI。

6.数值模拟结果分析：

图9为三个测试模型得到的波场快照，可以清楚的看到各向异性对波场造成的影响。图10为进行的GPU Direct P2P节点内的通信测试，图中的纵坐标为加速比，是通过GPU的计算时间与单颗CPU计算时间而得到的，可以清楚的看到，采用了GPU Direct P2P的节点内数据传输相较于传统的MPI 方法有着显著的计算效率提升。图11为进行的GPU DirectRDMA节点间的通信测试，其中图例所示的MPI-GDR表示的是GPU Direct RDMA的意思，进行了两个平台的测试，平台A每个节点1块K40C共三个节点，平台B每个节点3块K40C共三个阶段。测试结果表明GPU Direct RDMA对跨节点的通信有显著的提升效果，而且随着通信数量的增加其效果越发明显。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种三维各向异性弹性波数值模拟方法，其特征在于，具体包括以下步骤：

步骤1：建立介质模型，对介质模型进行网格离散得到多个网格点；

步骤2：计算震源函数，根据震源函数计算每个网格点上的压力值；

步骤3：将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

步骤4：根据波场值确定每个网格点的计算区域，进行分区并对分区边界数据进行数据交换，完成弹性波数值模拟；

所述步骤4具体包括以下步骤：

步骤4.1：根据所有波场值确定每个网格点的计算区域，并采用吸收边界的方式确定计算区域边界，在计算区域内模拟地下介质中波的传播；

步骤4.2：对所有计算区域进行分区，对相邻分区的边界数据进行数据交换，得到模拟弹性波数据，完成弹性波数值模拟。

2.根据权利要求1所述的一种三维各向异性弹性波数值模拟方法，其特征在于，所述步骤4.2中边界数据的采用GPU-Direct技术进行数据交换。

3.根据权利要求2所述的一种三维各向异性弹性波数值模拟方法，其特征在于，所述步骤4.2中的数据交换具体包括：

沿介质模型中计算区域边界数据变化最慢的轴向进行区域分解,得到多个分区，将所有分区分配到多个GPU中，每个分区独立在一个GPU上执行计算；每两个相邻分区的边界数据再进行交换。

4.根据权利要求1-3任一项所述的一种三维各向异性弹性波数值模拟方法，其特征在于，所述步骤1具体包括以下步骤：

步骤1.1：根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立介质模型；

步骤1.2：采用形状规则的三维网格对介质模型进行网格离散，得到网格点。

5.根据权利要求4所述的一种三维各向异性弹性波数值模拟方法，其特征在于，所述震源函数在空间上采用高斯函数，在时间上采用Ricker子波，所述震源函数的具体公式为：

s(x,y,z,t)＝g(x,y,z)·f(t) 公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²) 公式(2)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀,y₀,z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

6.根据权利要求4所述的一种三维各向异性弹性波数值模拟方法，其特征在于，所述步骤3具体包括以下步骤：

步骤3.1：将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长满足该数值格式的稳定性条件；

步骤3.2：采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值。

7.一种三维各向异性弹性波数值模拟系统，包括建模模块、震源模块、传播模块和数据交换模块；

所述建模模块用于建立介质模型，对介质模型进行网格离散得到多个网格点；

所述震源模块用于计算震源函数，根据震源函数计算每个网格点上的压力值；

所述传播模块用于将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

所述数据交换模块用于根据波场值确定每个网格点的计算区域，进行分区并对分区边界数据进行数据交换，完成弹性波数值模拟；

所述数据交换模块包括区域确定模块和分区交换模块；

所述区域确定模块用于根据所有波场值确定每个网格点的计算区域，并采用吸收边界的方式确定计算区域边界，在计算区域内模拟地下介质中波的传播；

所述分区交换模块用于对所有计算区域进行分区，对相邻分区的边界数据进行数据交换，得到模拟弹性波数据，完成弹性波数值模拟。

8.根据权利要求7所述的一种三维各向异性弹性波数值模拟系统，其特征在于，所述分区交换模块中边界数据的采用GPU-Direct技术进行数据交换。