三维复杂结构声波正演方法
【技术领域】
[0001] 本发明涉及地震勘探技术领域，特别涉及一种三维复杂结构声波正演方法。
【背景技术】
[0002] 石油和天然气的产量关乎国计民生，对地下构造的精细认识是指导油气勘探的关 键环节，现今主要的方法是地震勘探。通过观测和分析人工地震产生的地震波对地下构造 高精度成像，可有效得到地层、构造的分布，进而估算油气层厚度及分布范围。
[0003] 地震勘探包含数据采集、资料处理和资料解释三个环节，数据处理是将地震数据 变成地震剖面或者构造图，确定地下构造，指出有利的含油气地区。地震数据处理需要进行 大量数据运算，对计算效率要求极高。
[0004] 地震数据处理中包含去噪、反褶积、动校正、静校正、速度分析、叠加、偏移成像。地 震偏移成像是地震数据处理的核心，也是计算量最大的环节，通常也是实际生产的瓶颈。
[0005] 近年来，偏移成像技术不断发展，由人工偏移转变到计算机数字化偏移，由射线偏 移延伸到波动方程偏移，由叠后偏移拓展到叠前偏移，由二维扩展到三维等。
[0006] 逆时偏移能够对任意类型的波进行模拟成像，不存在倾角限制，并适用于任意复 杂介质，成为目前最精确的偏移成像技术之一。其原理是对双程波动方程进行精确求解，相 对于射线方法，不做过分近似。逆时偏移主要的计算量在于炮点的正传和检波点的正传，相 当于两次正演过程，因此三维复杂结构正演的效率决定了逆时偏移计算的效率。目前通过 大规模CHJ集群以及OpenMP(用于共享内存并行系统的多处理器程序设计的一套指导性的 编译处理方案）、MPI(Message Passing Interface，消息传递接口）等并行实现节点并行， 可以很大程度上提高运算效率，但在巨大计算需求的逆时偏移计算面前仍出现较多难题。 因此通过GPU(Graphics Processing Unit，图形处理器)-CPU异构系统进行三维复杂结构 声波正演加速显得尤为重要。
[0007] 近些年，全波形反演作为一种高精度速度模型构建的方法得到人们的普遍关注， 三维全波形反演技术中正演模拟部分计算量是限制其大规模应用的瓶颈之一，因此加速优 化三维复杂结构声波正演对三维全波形反演技术也至关重要。
[0008] 目前主要通过CPU集群以及OpenMP、MPI并行实现节点并行，也有使用GPU-CPU异构 系统加速，但未见使用较好的加速策略来使得GPU-CPU异构系统在三维复杂结构声波正演 中取得较高加速比。
【发明内容】
[0009] 本发明实施例提供了一种三维复杂结构声波正演方法，以解决现有技术中使用 GPU-CPU异构系统加速但加速比低的技术问题。该方法包括:获取改进后的图形处理器GPU-CPU异构系统，其中，该改进后的GPU-CPU异构系统包括以下改进内容:将波动方程离散为6 个差分方程，该6个差分方程包括一个用于计算X轴方向的应力的第一差分方程、一个用于 计算Y轴方向的应力的第二差分方程、一个用于计算Z轴方向的应力的第三差分方程、一个 用于计算X轴方向的质点运动速度的第四差分方程、一个用于计算Y轴方向的质点运动速度 的第五差分方程以及一个用于计算Ζ轴方向的质点运动速度的第六差分方程;在GPU内设计 4个计算程序，其中，一个计算程序用于承载第一差分方程的计算，一个计算程序用于承载 第二差分方程的计算，一个计算程序用于承载第三差分方程的计算，一个计算程序用于承 载第四差分方程、第五差分方程以及第六差分方程的计算;采用所述改进后的GPU-CPU异构 系统对待处理的地震数据进行正演。
[0022]其中，νΑΧ轴方向的质点运动速度;vAY轴方向的质点运动速度;Vz*Z轴方向的 质点运动速度;τχ为X轴方向的应力；Ty为Y轴方向的应力；τζ为Z轴方向的应力;P代表介质密 度，取值为l;a m为常量；At为模拟过程中的时间间隔；Δχ为X轴方向网格间隔；Ay为y轴方 向网格间隔；△ ζ轴为ζ方向网格间隔；i，j，k代表网格点的坐标位置;X，y，z分别代表X轴，Y 轴和Z轴方向;m为计算一个网格点所需要的基准网格点个数;L为计算一个网格点所需要的 基准网格点的最大个数;t为时刻，t+Ι为t+Ι时刻，表示向后推一个时刻，t-i代表t-i时刻，
表不向前推一个时刻：
，表不后推半个时刻： 前推半个时刻;d (X)是X轴方向的衰减系数，d (X) = do (i/pML )η，η为指数阶数，η为自然数，do =(111(1/1〇)(11+1)以(2?此），丨代表从区域边界到三维速度模型最内层的网格点数，？此为匹 配层的网格厚度;d(y)是Y轴方向的衰减系数，(1(7)=(1()(1八^)11;(1(2)是2轴方向的衰减系 数，d (z) = do (i /pml )η; ν为网格点的声波速度值;τ是网格点的应力。
[0023] 在一个实施例中，设计用于承载第一差分方程计算的计算程序，包括:将输入的待 处理的地震数据按照X轴方向连续的方式进行存储，数据块的大小为X轴方向上待处理的网 格数;采用连续编号的GHJ线程计算X轴方向连续的网格的应力，网格与线程一一对应，通过 与当前计算网格在X轴方向上相邻的网格的速度分量来计算当前计算网格的应力。
[0024] 在一个实施例中，设计用于承载第二差分方程计算的计算程序，包括:将输入的待 处理的地震数据按照X轴方向连续的方式进行存储，将输入的待处理的地震数据在X轴和Ζ 轴组成的平面上划分为二维数据块;采用连续编号的GPU线程计算Υ轴方向的网格的应力， 网格与线程-对应，通过与当前计算网格在Y轴方向上相邻的网格的速度分量来计算当 前计算网格的应力。
[0025] 在一个实施例中，设计用于承载第四差分方程、第五差分方程以及第六差分方程 计算的计算程序，包括:将输入的待处理的地震数据按照X轴方向连续的方式进行存储，将 输入的待处理的地震数据在X轴和Y轴组成的平面上划分为二维数据块;采用连续编号的 GPU线程分别计算X轴、Y轴以及Z轴方向的网格的质点运动速度，网格与线程一一对应，通过 与当前计算网格在X轴、Y轴以及Z轴方向上相邻的网格的应力来计算当前计算网格的质点 运动速度。
[0026] 在一个实施例中，在采用所述改进后的GPU-CPU异构系统对待处理的地震数据进 行正演之前，还包括:对所述待处理的地震数据，在X轴方向数据的首地址处填充字节，填充 后的X轴方向数据的字节数是32的倍数;在所述待处理的地震数据的头部填充数据，填充后 的待处理的地震数据的首地址按128字节对齐。
[0027] 在一个实施例中，通过以下公式计算在X轴方向数据的首地址处填充字节的长度： pad_x = 32-(nx+2pML+2r) %32;其中^&(1_1是在父轴方向数据的首地址处填充字节的长度， nx是X轴方向的网格数，Pml是匹配层厚度，r是计算当前计算网格时需要的相邻网格层数，r = d/2，d是空间导数离散精度。
[0028] 在一个实施例中，通过以下公式计算在所述待处理