一种三维各向异性弹性波数值模拟方法及系统的制作方法
【技术领域】
[0001] 本发明涉及一种三维各向异性弹性波数值模拟方法及系统，属于地球物理勘探领 域。
【背景技术】
[0002] 地震波数值模拟以地震波在地下介质中的传播理论为基础，在勘探地震学与天然 地震学中得到了广泛的应用。目前常规的三维声波方程、弹性波各向同性数值模拟已经广 泛的应用于数值模拟、成像、以及反演等各个地球物理领域中。但是对于如何有效的实现大 尺度采集方式(例如宽方位采集)和复杂各向异性介质（例如水平横向各向同性HTI或者是 正交各向异性）的三维弹性波波场数值模拟研究仍然存在很大的挑战，在实际应用中未得 到广泛的使用。此外常规的基于CPU的大规模三维声波、弹性波数值模拟，通常需要大量专 用集群计算资源。无论从硬件成本还是计算能耗上来说其成本都十分昂贵。因此快速的提 高计算性能、显著的降低计算成本对于实现三维弹性波各向异性数值模拟的实际应用具有 重要意义。
[0003] 近十年中，利用图形处理单元(GPU)进行计算密集型应用的加速实现已经得到了 突飞猛进式的发展。图形处理单元(GPU)由于其具有高速的内存带宽，相较于CPU至少高出 两个数量级的计算处理核心，更适合并行计算的单指令多数据(SMD)计算模式，以及更低 的能耗成本，正广泛的应用到计算科学的相关领域。对于勘探地球物理领域，对于使用图形 处理单元(GPU)的兴趣也在显著增强，越来越多的研究已经将GPU用于加速地震处理中的核 心算法，例如地震数值模拟、地震成像、地震高精度反演等。
【发明内容】
[0004] 经过研究发现，现有技术及其存在以下问题：
[0005] 将时间域有限差分方法应用到GHJ设备上实现复杂介质波传播算法上的研究很 少，尤其对于如何处理存在大量多节点数据交换的大规模三维问题研究就更少(Heinrich et al，2014)〇
[0006] 目前存在一些GPU地震波传播模拟的实现方案，对于三维GPU集群的实现方案我们 以（龙桂华等）的方案为例进行描述。该方案利用区域分解技术将单个GPU上不能计算的地 质体沿Z轴方向进行粗粒度分解，采用消息传递接口（MPI)交换边界数据，从而运用MPI + CUDA的方式实现了大尺度三维地震波场的数值模拟。但是该方法存在很大的一个问题是， 采用GHJ集群计算的加速比相较于单个GPU与CPU来说显著下降，造成这一结果的原因是GPU 计算的大部分耗时消耗在GPU集群节点间的GPU到CPU以及CPU到GPU设备之间的数据拷贝 上。
[0007] 本发明所要解决的技术问题是，针对现有技术的不足，提供一种将GPUDirect技术 与有限差分数值模拟方法结合得到的基于GPU Direct优化通信的三维各向异性弹性波数 值模拟方法及系统，为大规模数值模拟提供有力保证。
[0008] 本发明解决上述技术问题的技术方案如下:一种三维各向异性弹性波数值模拟方 法，具体包括以下步骤：
[0009] 步骤1:建立介质模型，对介质模型进行网格离散得到多个网格点；
[0010] 步骤2:计算震源函数，根据震源函数计算每个网格点上的压力值；
[0011] 步骤3:将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带 入传播方程进行计算，得到每一刻的波场值；
[0012] 步骤4:根据波场值确定每个网格点的计算区域，进行分区并对分区边界数据进行 数据交换，完成弹性波数值模拟。
[0013] 本发明的有益效果是:从应力应变方程出发，实现了利用图形处理单元(GPU)加速 复杂介质下的三维弹性波数值模拟，并针对三维问题引入GPU设备所面对的由于区域分解 所造成的多节点间、节点内通信瓶颈问题进行了深入的研究和分析，提出了利用GPU Direct技术加速数据传输的实现方案，避免了大量的CPU到GPU，GPU到CPU的数据拷贝，实现 了优化通信瓶颈的问题。通过GPU图形加速设备并引入GPU Direct通信优化策略，可以显著 的提升整体计算性能，可以用更低的硬件成本和更少的时间实现的三维各向异性弹性数值 模拟，为各种依赖于波动方程正演模拟的算法如逆时偏移，全波形反演的应用提供有力的 保证。
[0014] 在上述技术方案的基础上，本发明还可以做如下改进。
[0015] 进一步，所述步骤4具体包括以下步骤：
[0016] 步骤4.1:根据所有波场值确定每个网格点的计算区域，并采用吸收边界的方式确 定计算区域边界，在计算区域内模拟地下介质中波的传播；
[0017] 步骤4.2:对所有计算区域进行分区，对相邻分区的边界数据进行数据交换，得到 模拟弹性波数据，完成弹性波数值模拟。
[0018] 进一步，所述步骤4.2中边界数据的采用GPU-Direct技术进行数据交换。
[0019] 采用上述进一步方案的有益效果是，利用GPU加速三维弹性波数值模拟，显著提高 计算速度，相较于常规CHJ集群实现，可以达到20-30倍的速度提升。
[0020] 进一步，所述步骤4.2中的数据交换具体包括：
[0021] 沿介质模型中计算区域边界数据变化最慢的轴向进行区域分解，得到多个分区， 将所有分区分配到多个GPU中，每个分区独立在一个GPU上执行计算;每两个相邻分区的边 界数据再进行交换。
[0022]采用上述进一步方案的有益效果是，利用GPU Direct技术加速节点间、节点内数 据传输的实现方案，我们的方案因其避免了大量的CPU到GPU，GPU到CPU的数据拷贝，从而实 现了优化通信瓶颈的问题。
[0023] 进一步，所述步骤1具体包括以下步骤：
[0024] 步骤1.1:根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立 介质模型；
[0025]步骤1.2:采用形状规则的三维网格对介质模型进行网格离散，得到网格点。
[0026]采用上述进一步方案的有益效果是，将介质模型离散成网格点，网格点作为波场 值的载体，由三维结构的多个网格点构成波场，作为波场传播的初始条件。
[0027]进一步，所述震源函数在空间上采用高斯函数，在时间上采用Ricker子波，所述震 源函数的具体公式为：
[0028] s(x，y，z，t)=g(x，y，z) · f(t)公式（1)
[0029] 其中，以〇 =(卜2(31侃)2)6叉？（-(31侃）2)公式（2)
[0030]
[0031 ] 式中：t表示时间，fo表示Ricker子波的中心频率，模型计算中fo = 15Ηζ，β为常数； (义0,7〇,2())为震源的空间位置，1、7和2分别为1轴、7轴和2轴方向上的位置。
[0032] 进一步，所述步骤3具体包括以下步骤：
[0033]步骤3.1:将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限 差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长满足该数值格式的稳 定性条件；
[0034]步骤3.2:采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算， 得到每一刻的波场值。
[0035]本发明解决上述技术问题的技术方案如下:一种三维各向异性弹性波数值模拟系 统，包括建模模块、震源模块、传播模块和数据交换模块；
[0036] 所述建模模块用于建立介质模型，对介质模型进行网格离散得到多个网格点；
[0037] 所述震源模块用于计算震源函数，根据震源函数计算每个网格点上的压力值； [0038]所述传播模块用于将三维各向异性弹性波方程转换为传播方程，将每个网格点上 的压力值带入传播方程进行计算，得到每一刻的波场值；
[0039]所述数据交换模块用于根据波场值确定每个网格点的计算区域，进行分区并对分 区边界数