CN105005072B - 利用cuda的pml边界三维地震波传播模拟方法 - Google Patents

Abstract

本发明提供了一种利用CUDA的PML边界三维地震波传播模拟方法，方法包括：读入初始参数和雷克子波；将三维地震波的第一声波波动方程中的应力分解成相互垂直三个方向的应力分量，并加入PML边界条件形成第二声波波动方程；运用有限差分方法将第二声波波动方程离散化形成有限差分声波波动方程；给第一三维速度模型加边形成一第二三维速度模型；根据有限差分声波波动方程在中央处理器CPU中计算第二三维速度模型的棱、角、面的衰减系数；根据所述衰减系数在图像处理器GPU中计算所述第二三维速度模型网格点的应力；根据所述应力输出最后时刻的三维地震波模拟记录。本发明可以获得较高的加速比，缩短模拟时间。

Description

利用CUDA的PML边界三维地震波传播模拟方法

技术领域

本发明涉及地震波模拟技术，尤其涉及带PML(最佳匹配层)边界的三维GPU(图形处理器)地震波传播模拟技术，具体来说就是一种利用CUDA(Compute Unified DeviceArchitecture，计算统一设备架构)的PML边界三维地震波传播模拟方法。

背景技术

地震勘探是石油行业的上游业务，对胜利油田、大港油田以及华北油田的发现具有重大贡献，现今，在中国领土寻找石油储备日益困难的情况下，研发一种高精度的地震勘探方法是十分必要的。

地震勘探包含地震采集、数据处理环节，怎样从采集的数据中提取出有效信息进而得到地下构造的分布是数据处理的任务。现今常规的处理方式有时间偏移、克西霍夫深度偏移、单程波偏移。但是，由于中国的地质构造比较特殊，需要对规模较小的构造体呈像，比如，对于串珠状构造，精确知道其水平位置和垂直位置是指导钻井的关键。随时技术的进步，近几年出现的逆时偏移技术可以更加精确的成像，但其计算量太大，现在计算机无法实现，因而大幅度减少计算时间是能把较好的技术应用到石油勘探，进而寻找到更多油井的关键。

进一步来说，逆时偏移技术的核心是地震波传播模拟，只有具有快速准确的地震波传播模拟，才能实现成像效果好且实用的逆时偏移。

常规实现地震波模拟是使用中央处理器(CPU)计算，但是这种计算方法即便借助OpenMP和MPI(多点接口)等并行技术，计算依然十分耗时，不能达到使用要求，而且逆时偏移技术又需要多次使用地震波模拟过程，因此，耗时的地震波模拟使得计算量巨大的逆时偏移技术根本无法实际应用。

因此，本领域亟需一种将逆时偏移技术实际应用到地震波模拟过程中，并且模拟数据处理时间短的技术。

发明内容

本发明提供了一种利用CUDA的PML边界三维地震波传播模拟方法，采用边界吸收效果最好的PML技术，将三维速度模型加边形成新的三维速度模型，采用中央处理器(CPU)计算新三维速度模型的面、棱、角的衰减系数，采用图片处理器(GPU)计算新三维速度模型的网格点的质点速度值和应力，充分发挥CPU和GPU的各自优势，解决现有技术中地震波模拟数值计算时间过长的问题。

本发明的一种利用CUDA的PML边界三维地震波传播模拟方法，所述利用CUDA的PML边界三维地震波传播模拟方法包括：读入初始参数和雷克子波；将三维地震波的第一声波波动方程中的应力分解成相互垂直三个方向的应力分量，并加入PML边界条件形成第二声波波动方程；运用有限差分方法将所述第二声波波动方程离散化形成有限差分声波波动方程；给所述第一三维速度模型加边形成一第二三维速度模型；根据有限差分声波波动方程在中央处理器(CPU)中计算所述第二三维速度模型的棱、角、面的衰减系数；根据所述衰减系数在图像处理器(GPU)中计算所述第二三维速度模型网格点的应力；以及根据所述应力输出最后时刻的三维地震波模拟记录。

本发明提供了一种利用CUDA的PML边界三维地震波传播模拟方法，采用边界吸收效果最好的PML技术，将三维速度模型加边形成新的三维速度模型，采用中央处理器(CPU)计算新三维速度模型的面、棱、角的衰减系数，采用图片处理器(GPU)计算新三维速度模型的网格点的质点速度值和应力，同时考虑到GPU内存类型不同特点，当需要计算下一个网格点的值时，仅需要更新部分寄存器内的数值，大大提高了计算效率，充分发挥CPU和GPU的各自优势，减少了机时和能耗，达到较高的加速比。

应了解的是，上述一般描述及以下具体实施方式仅为示例性及阐释性的，其并不能限制本发明所欲主张的范围。

附图说明

下面的所附附图是本发明的说明书的一部分，其绘示了本发明的示例实施例，所附附图与说明书的描述一起用来说明本发明的原理。

图1为本发明实施例提供的一种利用CUDA的PML边界三维地震波传播模拟方法的实施方式一的流程图；

图2为本发明实施例提供的一种利用CUDA的PML边界三维地震波传播模拟方法的实施方式二的流程图；

图3为本发明实施例提供的一种在GPU中计算所述第二三维速度模型网格点的质点运动速度和应力的流程图；

图4为本发明实施例提供的一种利用GPU内不同类型内存的性质计算第二三维速度模型的网格点的质点运动速度和应力的示意图；

图5为本发明实施例提供的一种利用CUDA的PML边界三维地震波传播模拟装置的框图；

图6为本发明实施例提供的利用传统CPU计算出的波场快照图；

图7为本发明实施例提供的利用本发明计算出的波场快照图。

附图符号说明：

10 震源产生单元 20 应力分解单元

30 有限差分处理单元 40 加边单元

50 第一计算单元 60 第二计算单元

70 记录单元 80 分析单元

S101～S107 步骤 S4061～S4064 步骤

S401～S407 步骤

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面将以附图及详细叙述清楚说明本发明所揭示内容的精神，任何所属技术领域技术人员在了解本发明内容的实施例后，当可由本发明内容所教示的技术，加以改变及修饰，其并不脱离本发明内容的精神与范围。

本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。另外，在附图及实施方式中所使用相同或类似标号的元件/构件是用来代表相同或类似部分。

关于本文中所使用的“第一”、“第二”、…等，并非特别指称次序或顺位的意思，也非用以限定本发明，其仅为了区别以相同技术用语描述的元件或操作。

关于本文中所使用的方向用语，例如：上、下、左、右、前或后等，仅是参考附图的方向。因此，使用的方向用语是用来说明并非用来限制本创作。

关于本文中所使用的“包含”、“包括”、“具有”、“含有”等等，均为开放性的用语，即意指包含但不限于。

关于本文中所使用的“及/或”，包括所述事物的任一或全部组合。

关于本文中所使用的用语“大致”、“约”等，用以修饰任何可以微变化的数量或误差，但这些微变化或误差并不会改变其本质。一般而言，此类用语所修饰的微变化或误差的范围在部分实施例中可为20％，在部分实施例中可为10％，在部分实施例中可为5％或是其他数值。本领域技术人员应当了解，前述提及的数值可依实际需求而调整，并不以此为限。

关于本文中所使用的用词(terms)，除有特别注明外，通常具有每个用词使用在此领域中、在此申请的内容中与特殊内容中的平常意义。某些用以描述本申请的用词将于下或在此说明书的别处讨论，以提供本领域技术人员在有关本申请的描述上额外的引导。

图1为本发明实施例提供的一种利用CUDA的PML边界三维地震波传播模拟方法的实施方式一的流程图，如图1所示，一种利用计算统一设备架构(CUDA)的最佳匹配层(PML)边界三维地震波传播模拟方法，所述利用CUDA的PML边界三维地震波传播模拟方法包括：

S101：读入初始参数和雷克子波；

S102：将三维地震波的第一声波波动方程中的应力分解成相互垂直三个方向的应力分量，并加入PML边界条件形成第二声波波动方程；

S103：运用有限差分方法将所述第二声波波动方程离散化形成有限差分声波波动方程；

S104：给所述第一三维速度模型加边形成一第二三维速度模型；

S105：根据有限差分声波波动方程在中央处理器CPU中计算所述第二三维速度模型的棱、角、面的衰减系数；

S106：根据所述衰减系数在图像处理器GPU中计算所述第二三维速度模型网格点的应力；以及

S107：根据所述应力输出最后时刻的三维地震波模拟记录。

参照图1，本发明提供的地震波模拟方法可以精确地模拟地震波在介质中的传播，加入的PML边界条件能有效消除人工边界反射，扩展到三维能满足实际勘探的需求，采用优化的策略运用GPU进行计算可以大大提高计算效率，相比同价位的CPU计算，原来需要7天能完成的计算，运用本申请公开的技术方案所述模拟流程只需要1小时，大大减少了机时和能耗。

本发明一具体实施例中，所述初始参数包括：炮点位置、检波点位置、时间步长和地震波模拟时长等。这些可根据实验人员的需求合理设定，能够满足各种模拟实验需求，扩大了本发明的适用范围。

本发明的另一个具体实施例中，所述三维地震波模拟记录可以为波场快照图，所述波场快照图是根据所述应力的大小绘制出来的。通过波场快照图可以直观看到三维地震波模拟技术是否能够很好地消除边界反射影响；根据应力的大小绘制出波场快照图，实现起来简单方便，可操作性强。

图2为本发明实施例提供的一种利用CUDA的PML边界三维地震波传播模拟方法的实施方式二的流程图，如图2所示，其中，步骤S401与步骤S101相同，步骤S402与步骤S102相同，步骤S403与步骤S103相同，步骤S404与步骤S104相同，步骤S405与步骤S105相同，所述利用CUDA的PML边界三维地震波传播模拟方法还包括：

S406：根据所述衰减系数在图像处理器GPU中计算所述第二三维速度模型网格点的质点运动速度；以及

S407：根据所述质点运动速度分析地质构造。

参照图2，第二三维速度模型网格点的质点运动速度能够很好地反映网格点处的地质构造，从而指导石油勘探，进而寻找到更多油井。

本发明的一具体实施例中，所述第一声波波动方程为：

$\frac{\∂\mathrm{\τ}}{\∂t}=-\mathrm{\ρ}{v}^2(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y}+\frac{\∂v_z}{\∂z}),$

$\frac{\∂v_x}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂\mathrm{\τ}}{\∂x},$

$\frac{\∂v_y}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂\mathrm{\τ}}{\∂y},$

$\frac{\∂v_z}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂\mathrm{\τ}}{\∂z}$

其中，v为网格点的声波速度值；ρ代表介质密度，取值为1；v_x代表x方向的质点运动速度，v_y代表y方向的质点运动速度，v_z代表z方向的质点运动速度，τ代表网格点的应力，t表示地震波声波传播时间；

所述第二声波波动方程为：

τ＝τ_x+τ_y+τ_z,

$\frac{\∂{\mathrm{\τ}}_x}{\∂t}+d\left(x\right){\mathrm{\τ}}_x=-\mathrm{\ρ}{v}^2\frac{\∂v_x}{\∂x},$

$\frac{\∂{\mathrm{\τ}}_y}{\∂t}+d\left(y\right){\mathrm{\τ}}_y=-\mathrm{\ρ}{v}^2\frac{\∂v_y}{\∂y},$

$\frac{\∂{\mathrm{\τ}}_z}{\∂t}+d\left(z\right){\mathrm{\τ}}_z=-\mathrm{\ρ}{v}^2\frac{\∂v_z}{\∂z},$

$\mathrm{\ρ}\frac{\∂v_x}{\∂t}+d\left(x\right)v_x=-\frac{1}{\mathrm{\ρ}}(\frac{\∂{\mathrm{\τ}}_x}{\∂x}+\frac{\∂{\mathrm{\τ}}_y}{\∂x}+\frac{\∂{\mathrm{\τ}}_z}{\∂x}),$

$\mathrm{\ρ}\frac{\∂v_y}{\∂t}+d\left(y\right)v_y=-\frac{1}{\mathrm{\ρ}}(\frac{\∂{\mathrm{\τ}}_x}{\∂y}+\frac{\∂{\mathrm{\τ}}_y}{\∂y}+\frac{\∂{\mathrm{\τ}}_z}{\∂y}),$

$\mathrm{\ρ}\frac{\∂v_z}{\∂t}+d\left(z\right)v_z=-\frac{1}{\mathrm{\ρ}}(\frac{\∂{\mathrm{\τ}}_x}{\∂z}+\frac{\∂{\mathrm{\τ}}_y}{\∂z}+\frac{\∂{\mathrm{\τ}}_z}{\∂z})$

其中，d(x)是X方向的衰减系数，能够表示为d(x)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；d(y)是Y方向的衰减系数，能够表示为d(y)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；d(z)是Z方向的衰减系数，能够表示为d(z)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；n为指数阶数；v为网格点的声波速度值；ρ代表介质密度，取值为1；t表示地震波声波传播时间；R为反射系数，一般取值为0.0001；v_x为X方向的质点运动速度；v_y为Y方向的质点运动速度；v_z为Z方向的质点运动速度；τ_x为X方向的应力；τ_y为Y方向的应力；τ_z为Z方向的应力。

所述有限差分声波波动方程为：

$\frac{{\mathrm{\τ}}_{x(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i,j,k)}^{t-\frac{1}{2}}}{\mathrm{\Δt}}+d_i^x\frac{{\mathrm{\τ}}_{x(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i,j,k)}^{t-\frac{1}{2}}}{2}=-\mathrm{\ρ}{v}^2\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{v_{x(i+\frac{1}{2}+m-1,j,k)}^t-v_{x(i-\frac{1}{2}-m+1,j,k)}^t}{\mathrm{\Δx}},$

$\frac{{\mathrm{\τ}}_{y(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{y(i,j,k)}^{t-\frac{1}{2}}}{\mathrm{\Δt}}+d_j^y\frac{{\mathrm{\τ}}_{y(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{y(i,j,k)}^{t-\frac{1}{2}}}{2}=-\mathrm{\ρ}{v}^2\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{v_{y(i,j+\frac{1}{2}+m-1,k)}^t-v_{y(i,j-\frac{1}{2}-m+1,k)}^t}{\mathrm{\Δy}},$

$\frac{{\mathrm{\τ}}_{z(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{z(i,j,k)}^{t-\frac{1}{2}}}{\mathrm{\Δt}}+d_k^z\frac{{\mathrm{\τ}}_{y(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{y(i,j,k)}^{t-\frac{1}{2}}}{2}=-\mathrm{\ρ}{v}^2\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{v_{z(i,j,k+\frac{1}{2}+m-1)}^t-v_{z(i,j,k-\frac{1}{2}-m+1)}^t}{\mathrm{\Δz}},$

$\frac{v_{x(i+\frac{1}{2},j,k)}^{t+1}-v_{x(i+\frac{1}{2},j,k)}^t}{\mathrm{\Δt}}+d_i^x\frac{v_{x(i+\frac{1}{2},j,k)}^{t+1}-v_{x(i+\frac{1}{2},j,k)}^t}{2}=-\frac{1}{\mathrm{\ρ}}\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{{\mathrm{\τ}}_{x(i+1+m-1,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i-m+1,j,k)}^{t+\frac{1}{2}}}{\mathrm{\Δx}},$

$\frac{v_{x(i,j+\frac{1}{2},k)}^{t+1}-v_{x(i,j+\frac{1}{2},k)}^t}{\mathrm{\Δt}}+d_j^y\frac{v_{x(i+\frac{1}{2},j,k)}^{t+1}-v_{x(i+\frac{1}{2},j,k)}^t}{2}=-\frac{1}{\mathrm{\ρ}}\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{{\mathrm{\τ}}_{y(i,j+1+m-1,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i,j-m+1,k)}^{t+\frac{1}{2}}}{\mathrm{\Δy}},$

$\frac{v_{x(i,j,k+\frac{1}{2})}^{t+1}-v_{x(i,j,k+\frac{1}{2})}^t}{\mathrm{\Δt}}+d_k^z\frac{v_{x(i,j,k+\frac{1}{2})}^{t+1}-v_{x(i,j,k+\frac{1}{2})}^t}{2}=-\frac{1}{\mathrm{\ρ}}\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{{\mathrm{\τ}}_{z(i,j,k+1+m-1)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i,j,k-m+1)}^{t+\frac{1}{2}}}{\mathrm{\Δz}}$

其中，v_x为X方向的质点运动速度；v_y为Y方向的质点运动速度；v_z为Z方向的质点运动速度；τ_x为X方向的应力；τ_y为Y方向的应力；τ_z为Z方向的应力；ρ代表介质密度，取值为1；a_m为常量；Δt为模拟过程中的时间间隔；Δx为x方向网格间隔；Δy为y方向网格间隔；Δz为z方向网格间隔；i,j,k代表网格点的坐标位置，例如x(1,2,3)表示X方向的坐标为1，2，3；x,y,z分别代表x轴,y轴和z轴方向；m为计算一个网格点所需要的基准网格点个数；L为计算一个网格点所需要的基准网格点的最大个数，实际操作中最大值为16；t为时刻，t+1为t+1时刻，即向后推一个时刻，t-1代表t-1时刻，即向前推一个时刻，代表时刻，即后推半个时刻，代表时刻，即向前推半个时刻；d(x)是X方向的衰减系数，能够表示为d(x)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度，其中，d(x)在此方程中还可以表示为即X轴方向i点处的衰减系数；d(y)是Y方向的衰减系数，能够表示为d(y)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度，其中，d(y)在此方程中还可以表示为即Y轴方向j点处的衰减系数；d(z)是Z方向的衰减系数，能够表示为d(z)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度，其中，d(z)在此方程中还可以表示为即Z轴方向k点处的衰减系数；n为指数阶数；v为网格点的声波速度值。

本发明的另一具体实施例中，在CPU中计算所述第二三维速度模型的棱、角、面的衰减系数具体包括：

在X轴方面的两个表面上，仅在X方向上衰减，此时，d(x)≠0，d(y)＝0，d(z)＝0；

在Y轴方面的两个表面上，仅在Y方向上衰减，此时，d(x)＝0，d(y)≠0，d(z)＝0；

在Z轴方面的两个表面上，仅在Z方向上衰减，此时，d(x)＝0，d(y)＝0，d(z)≠0；

与X轴平行的四条棱，仅在X方向上没有衰减，此时，d(x)＝0，d(y)≠0，d(z)≠0；

与Y轴平行的四条棱，仅在Y方向上没有衰减，此时，d(x)≠0，d(y)＝0，d(z)≠0；

与Z轴平行的四条棱，仅在Z方向上没有衰减，此时，d(x)≠0，d(y)≠0，d(z)＝0；

角在X、Y、Z方向上均有衰减，此时，d(x)≠0，d(y)≠0，d(z)≠0。

由于衰减系数是由所要计算的网格点的位置(在第二三维速度模型上的坐标)以及该网格点处的速度值确定，因此在第二三维速度模型上，不同的棱、角、边具有不同的衰减表达。举例来说，在上表面d(x)＝0，即在x方向不需要衰减，d(y)＝0，即在y方向也不需要衰减，d(z)≠0，即在z方向衰减，具体的d(z)值由公式d(z)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)计算，同样的道理，进而求出第二三维速度模型的下表面、左边的面、右边的面、前面、后面、上部靠前的棱、上部靠后的棱、上部左边的棱、上部右边的棱、竖着的四条棱、下面的四条棱以及八个顶点分别x，y，z的衰减系数，因此，共有12条棱、8个角和6个面需要分别计算衰减系数。由于衰减系数的计算设计过多的判断循环，不适合在图像处理器(GPU)中计算，因此衰减系数的计算在中央处理器(CPU)中计算，整体计算完成后传入GPU内存，这样处理使得内存占用略微增加但是换来的是计算效率的显著提高。

图3为本发明实施例提供的一种在GPU中计算所述第二三维速度模型网格点的质点运动速度和应力的流程图，如图3所示，所述GPU的内存类型包括本地存储器(LocalMemory)和共享存储器(Shared Memory)，可以根据GPU中不同类型内存的性质合理使用本地存储器和共享存储器，来计算所述第二三维速度模型的网格点的质点运动速度和应力，在GPU中计算所述第二三维速度模型网格点的质点运动速度和应力具体包括：

S4061：将需要计算的当前网格点的属性值读入本地存储器；

S4062：基于所述本地存储器中存储的属性值计算当前网格点的质点运动速度和应力；

S4063：所述本地存储器从共享存储器中读入下一网格点与当前网格点不相同的属性值；以及

S4064：基于所述本地存储器中当前存储的属性值计算下一网格点的质点运动速度和应力。

参照图3，GPU的内存类型包括本地存储器(Local Memory)和共享存储器(SharedMemory)，可以本地存储器和共享存储器的不同性质，充分加速所述第二三维速度模型网格点的质点运动速度和应力的计算，本发明的一具体实施例中，所述属性值可以为网格点的质点运动速度和应力，本发明不以此为限。

图4为本发明实施例提供的一种利用GPU内不同类型内存的性质计算第二三维速度模型的网格点的质点运动速度和应力的示意图，如图4所示，由于图像处理器(GPU)中不同内存的差异，需要合适的策略才能充分使用GPU的计算能力，做差分的过程是逐个网格点的计算，其间要用到周围的网格点值，Local Memory读取速度最快，首先把需要计算的网格点读入Local Memory中，当需计算下一个网格点的值时，只需从Shared Memory中读入一个点即可计算下一个网格点的值，利用该策略可以大大提高计算效率，如图4所示，假设根据8个相邻网格点的质点运动速度值和应力可以获得一个网格点的质点运动速度值和应力，图4中，Reg[4]所指的黑色方块为存储计算出来的当前网格点的质点运动速度值和应力，当前网格点的质点运动速度值和应力的计算需要借助其余8个相邻网格点的质点运动速度值和应力，左侧中Reg[4]所指的黑色方块为寄存器Reg[4]存储的计算出的当前网格点的质点运动速度值和应力，当需要计算下一网格点的质点运动速度值和应力时，寄存器(Reg[0]，Reg[1]，Reg[2]，…Reg[8])所有指针上移一位，只有寄存器Reg[8]从Shared Memory中读入一个点(图4右侧最上方的方块)，即可获得下一个网格点的质点运动速度值和应力(图4右侧Reg[4]所指的黑色方块)，极大节省了读取时间，进而节省了三维速度模型网格点数值的处理时间。

参照图3和图4，根据不同计算类型选择使用不同类型GPU内存，其中，ConstantMemory存在于显存，只读，可被grid中所有thread读取；Shared Memory是GPU中最重要的内存，类似于CPU中的cache，可读写，非常快；Register 3.0硬件有63个寄存器，是thread私有；Global Memory是访问最慢(400-600周期)的GPU内存。计算第二三维速度模型的面、棱、角的衰减系数时，需要计算大量的判断循环，这种计算并不是GPU所擅长的，但这种大量逻辑计算却是CPU所擅长的，因此，首先在CPU中计算衰减系数，衰减系数计算完成后，转入GPU中计算第二三维速度模型的网格点的质点运动速度和应力。因此，本发明根据PML边界三维地震波传播模拟过程中不同阶段数据运算的性质，选择CPU和GPU进行配合运算，极大提高了计算效率。另外，由于GPU的内存类型包括本地存储器(Local Memory)和共享存储器(Shared Memory)，本地存储器(图4中由多个寄存器Reg[0]，Reg[1]，Reg[2]，…Reg[8]组成)读取速度最快，共享存储器类似于CPU中的cache，基本上可以配合本地存储器运算，本地存储器从共享存储器中读取数据(图4中箭头所标识的方块表示共享存储器中存储的数据)，由于计算第二三维速度模型的网格点的质点运动速度和应力是迭代运算，因此每次只需更新一个寄存器(Register)内数据值即可，这样就促使GPU的运算速度快的优势得以最大发挥，进一步提高了计算效率。

本发明的一个具体实施例中，给所述第一三维速度模型加边形成第二三维速度模型具体包括：

所述第二三维速度模型的面的PML区域由所述第一三维速度模型的表面最上层的网格点的声波速度值映射到PML区域，棱的PML区域由所述第一三维速度模型对应棱的声波速度值映射到PML区域，所述第二三维速度模型的角的PML区域由所述第一三维速度模型的对应角的声波速度值映射到PML区域。

由于第一三维速度模型的面、棱、角的速度值都是已知的，由于本发明采用边界吸收效果最好的PML技术，需要完全吸收来自边界的反射，同时需要全面分析第一三维速度模型内部及面、棱、角的网格点的质点运动速度和应力，进而分析第一三维速度模型的内部地质构造，所以需要给第一三维速度模型加边，即在第一三维速度模型外围再形成一个三维速度模型(即第二三维速度模型)，并且两个三维速度模型均为正立方体，中心点重合，而且对应面平行，这样通过第二三维速度模型就可以很好地反映第一三维速度模型面、棱、角的网格点的质点运动速度和应力了，从而达到对第一三维速度模型的全面分析。为了让地震波的质点运动速度和应力在第一三维速度模型和第二三维速度模型之间的区域内刚好衰减到0(这样就不会出现边界反射，而且能够进行地震波模拟)，需要将第一三维速度模型面、棱、角的速度值映射到第二三维速度模型对应的面、棱、角上，从而将PML技术成功扩展至三维。

图5为本发明实施例提供的一种利用CUDA的PML边界三维地震波传播模拟装置的框图，如图5所示，所述利用CUDA的PML边界三维地震波传播模拟装置包括震源产生单元10、应力分解单元20、有限差分处理单元30、加边单元40、第一计算单元50、第二计算单元60、记录单元70和分析单元80，其中，震源产生单元10用于读入初始参数和雷克子；应力分解单元20将三维地震波的第一声波波动方程中的应力分解成相互垂直三个方向的应力分量，并加入PML边界条件形成第二声波波动方程；有限差分处理单元30用于运用有限差分方法将所述第二声波波动方程离散化形成有限差分声波波动方程；加边单元40用于给所述第一三维速度模型加边形成一第二三维速度模型；第一计算单元50用于根据有限差分声波波动方程在中央处理器CPU中计算所述第二三维速度模型的棱、角、面的衰减系数；第二计算单元60用于根据所述衰减系数在图像处理器GPU中计算所述第二三维速度模型网格点的质点运动速度和应力；记录单元70用于根据所述应力输出最后时刻的三维地震波模拟记录；分析单元80用于根据所述质点运动速度分析地质构造。

参照图5，本发明提供的地震波模拟装置可以精确地模拟地震波在介质中的传播，加入的PML边界条件能有效消除人工边界反射，扩展到三维能满足实际勘探的需求，采用优化的策略运用GPU进行计算可以大大提高计算效率，相比同价位的CPU计算，原来需要7天能完成的计算，运用本申请公开的技术方案所述模拟流程只需要1小时，大大减少了机时和能耗；根据GPU的内存类型采取优化手段，进一步加速地震波模拟，使得地震波模拟时间大大缩短。

本发明的一具体实施例中，所述第一声波波动方程为：

$\frac{\∂\mathrm{\τ}}{\∂t}=-\mathrm{\ρ}{v}^2(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y}+\frac{\∂v_z}{\∂z}),$

$\frac{\∂v_x}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂\mathrm{\τ}}{\∂x},$

$\frac{\∂v_y}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂\mathrm{\τ}}{\∂y},$

$\frac{\∂v_z}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂\mathrm{\τ}}{\∂z}$

其中，v为网格点的声波速度值；ρ代表介质密度，取值为1；v_x代表x方向的质点运动速度，v_y代表y方向的质点运动速度，v_z代表z方向的质点运动速度，τ代表网格点的应力，t表示地震波的声波传播时间；

所述第二声波波动方程为：

τ＝τ_x+τ_y+τ_z,

$\frac{\∂{\mathrm{\τ}}_x}{\∂t}+d\left(x\right){\mathrm{\τ}}_x=-\mathrm{\ρ}{v}^2\frac{\∂v_x}{\∂x},$

$\frac{\∂{\mathrm{\τ}}_y}{\∂t}+d\left(y\right){\mathrm{\τ}}_y=-\mathrm{\ρ}{v}^2\frac{\∂v_y}{\∂y},$

$\frac{\∂{\mathrm{\τ}}_z}{\∂t}+d\left(z\right){\mathrm{\τ}}_z=-\mathrm{\ρ}{v}^2\frac{\∂v_z}{\∂z},$

$\mathrm{\ρ}\frac{\∂v_x}{\∂t}+d\left(x\right)v_x=-\frac{1}{\mathrm{\ρ}}(\frac{\∂{\mathrm{\τ}}_x}{\∂x}+\frac{\∂{\mathrm{\τ}}_y}{\∂x}+\frac{\∂{\mathrm{\τ}}_z}{\∂x}),$

$\mathrm{\ρ}\frac{\∂v_y}{\∂t}+d\left(y\right)v_y=-\frac{1}{\mathrm{\ρ}}(\frac{\∂{\mathrm{\τ}}_x}{\∂y}+\frac{\∂{\mathrm{\τ}}_y}{\∂y}+\frac{\∂{\mathrm{\τ}}_z}{\∂y}),$

$\mathrm{\ρ}\frac{\∂v_z}{\∂t}+d\left(z\right)v_z=-\frac{1}{\mathrm{\ρ}}(\frac{\∂{\mathrm{\τ}}_x}{\∂z}+\frac{\∂{\mathrm{\τ}}_y}{\∂z}+\frac{\∂{\mathrm{\τ}}_z}{\∂z})$

其中，d(x)是X方向的衰减系数，能够表示为d(x)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；v为网格点的速度值；d(y)是Y方向的衰减系数，能够表示为d(y)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；d(z)是Z方向的衰减系数，能够表示为d(z)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；n为指数阶数；v为网格点的声波速度值；ρ代表介质密度，取值为1；t表示地震波声波传播时间；R为反射系数，一般取值为0.0001；v_x为X方向的质点运动速度；v_y为Y方向的质点运动速度；v_z为Z方向的质点运动速度；τ_x为X方向的应力；τ_y为Y方向的应力；τ_z为Z方向的应力；

所述有限差分声波波动方程为：

$\frac{{\mathrm{\τ}}_{x(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i,j,k)}^{t-\frac{1}{2}}}{\mathrm{\Δt}}+d_i^x\frac{{\mathrm{\τ}}_{x(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i,j,k)}^{t-\frac{1}{2}}}{2}=-\mathrm{\ρ}{v}^2\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{v_{x(i+\frac{1}{2}+m-1,j,k)}^t-v_{x(i-\frac{1}{2}-m+1,j,k)}^t}{\mathrm{\Δx}},$

$\frac{{\mathrm{\τ}}_{y(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{y(i,j,k)}^{t-\frac{1}{2}}}{\mathrm{\Δt}}+d_j^y\frac{{\mathrm{\τ}}_{y(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{y(i,j,k)}^{t-\frac{1}{2}}}{2}=-\mathrm{\ρ}{v}^2\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{v_{y(i,j+\frac{1}{2}+m-1,k)}^t-v_{y(i,j-\frac{1}{2}-m+1,k)}^t}{\mathrm{\Δy}},$

$\frac{{\mathrm{\τ}}_{z(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{z(i,j,k)}^{t-\frac{1}{2}}}{\mathrm{\Δt}}+d_k^z\frac{{\mathrm{\τ}}_{y(i,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{y(i,j,k)}^{t-\frac{1}{2}}}{2}=-\mathrm{\ρ}{v}^2\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{v_{z(i,j,k+\frac{1}{2}+m-1)}^t-v_{z(i,j,k-\frac{1}{2}-m+1)}^t}{\mathrm{\Δz}},$

$\frac{v_{x(i+\frac{1}{2},j,k)}^{t+1}-v_{x(i+\frac{1}{2},j,k)}^t}{\mathrm{\Δt}}+d_i^x\frac{v_{x(i+\frac{1}{2},j,k)}^{t+1}-v_{x(i+\frac{1}{2},j,k)}^t}{2}=-\frac{1}{\mathrm{\ρ}}\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{{\mathrm{\τ}}_{x(i+1+m-1,j,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i-m+1,j,k)}^{t+\frac{1}{2}}}{\mathrm{\Δx}},$

$\frac{v_{x(i,j+\frac{1}{2},k)}^{t+1}-v_{x(i,j+\frac{1}{2},k)}^t}{\mathrm{\Δt}}+d_j^y\frac{v_{x(i+\frac{1}{2},j,k)}^{t+1}-v_{x(i+\frac{1}{2},j,k)}^t}{2}=-\frac{1}{\mathrm{\ρ}}\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{{\mathrm{\τ}}_{y(i,j+1+m-1,k)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i,j-m+1,k)}^{t+\frac{1}{2}}}{\mathrm{\Δy}},$

$\frac{v_{x(i,j,k+\frac{1}{2})}^{t+1}-v_{x(i,j,k+\frac{1}{2})}^t}{\mathrm{\Δt}}+d_k^z\frac{v_{x(i,j,k+\frac{1}{2})}^{t+1}-v_{x(i,j,k+\frac{1}{2})}^t}{2}=-\frac{1}{\mathrm{\ρ}}\underset{m=1}{\overset{L}{\mathrm{\Σ}}}{a}_m\frac{{\mathrm{\τ}}_{z(i,j,k+1+m-1)}^{t+\frac{1}{2}}-{\mathrm{\τ}}_{x(i,j,k-m+1)}^{t+\frac{1}{2}}}{\mathrm{\Δz}}$

其中，v_x为X方向的质点运动速度；v_y为Y方向的质点运动速度；v_z为Z方向的质点运动速度；τ_x为X方向的应力；τ_y为Y方向的应力；τ_z为Z方向的应力；ρ代表介质密度，取值为1；a_m为常量；Δt为模拟过程中的时间间隔；Δx为x方向网格间隔；Δy为y方向网格间隔；Δz为z方向网格间隔；i,j,k代表网格点的坐标位置，例如x(1,2,3)表示X方向的坐标为1，2，3；x,y,z分别代表x轴,y轴和z轴方向；m为计算一个网格点所需要的基准网格点个数；L为计算一个网格点所需要的基准网格点的最大个数，实际操作中最大值为16；t为时刻，t+1为t+1时刻，即向后推一个时刻，t-1代表t-1时刻，即向前推一个时刻，代表时刻，即后推半个时刻，代表时刻，即向前推半个时刻；d(x)是X方向的衰减系数，能够表示为d(x)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度，其中，d(x)在此方程中还可以表示为即X轴方向i点处的衰减系数；d(y)是Y方向的衰减系数，能够表示为d(y)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度，其中，d(y)在此方程中还可以表示为即Y轴方向j点处的衰减系数；d(z)是Z方向的衰减系数，能够表示为d(z)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度，其中，d(z)在此方程中还可以表示为即Z轴方向k点处的衰减系数；n为指数阶数；v为网格点的声波速度值。

图6为本发明实施例提供的利用传统CPU计算出的波场快照图；图7为本发明实施例提供的利用本发明计算出的波场快照图，参照图6和图7，在256×256×256规模的三维速度模型上进行地震波传播模拟，在只用CPU的情况下，计算机配置为Inter(R)Core(TM)i7-2600CPU 3.40GHz,8G内存，电脑系统为Linux 2.6.32-220.el6，PML为32层，时间步数500，三个方向采样间隔均为3米，子波主频28Hz，需要3个小时，达到的模拟效果如图6所示；如果利用本发明的模拟方法，电脑配置一样，GPU显卡型号为GTX680，仅仅需要一分钟，达到的模拟效果如图7所示，图6和图7是一样的。具体来说，波场快照图为网格点上的应力，因此，本发明的一具体实施例中，所述三维地震波模拟记录可以为波场快照图。

本发明的一具体实施例中，在CPU中利用OpenMP加速计算所述第二三维速度模型的棱、角、面的所述衰减系数。其中，OpenMP为一种用于共享内存并行系统的多线程程序设计的库，特别适合于多核CPU上的并行程序开发设计。利用OpenMP并行运算第二三维速度模型的棱、角、面的所述衰减系数，极大地缩短了处理时间。

本发明提供了一种利用CUDA的PML边界三维地震波传播模拟方法及装置，采用边界吸收效果最好的PML技术，将三维速度模型加边形成新的三维速度模型，采用中央处理器(CPU)计算新三维速度模型的面、棱、角的衰减系数，采用图片处理器(GPU)计算新三维速度模型的网格点的质点运动速度和应力，同时考虑到GPU内存类型不同特点，当需要计算下一个网格点的值时，仅需要更新部分寄存器内的数值，大大提高了计算效率，充分发挥CPU和GPU的各自优势，减少了机时和能耗，达到较高的加速比。

本发明至少还具有以下有益效果：

1.实现三维情况下PML边界地震波GPU模拟是难点；

2.合理利用GPU的内部计算资源，利用GPU内存优化策略，达到了较高的加速比，使得原本需要以天计的机时缩短到1小时，如果扩展到集群，将是巨大的资源节约；

3.实现三维PML边界的地震波GPU模拟过程以及获得较高加速比。

上述的本发明实施例可在各种硬件、软件编码或两者组合中进行实施。例如，本发明的实施例也可为在数据信号处理器(Digital Signal Processor,DSP)中执行的执行上述程序的程序代码。本发明也可涉及计算机处理器、数字信号处理器、微处理器或现场可编程门阵列(Field Programmable Gate Array,FPGA)执行的多种功能。可根据本发明配置上述处理器执行特定任务，其通过执行定义了本发明揭示的特定方法的机器可读软件代码或固件代码来完成。可将软件代码或固件代码发展为不同的程序语言与不同的格式或形式。也可为了不同的目标平台编译软件代码。然而，根据本发明执行任务的软件代码与其他类型配置代码的不同代码样式、类型与语言不脱离本发明的精神与范围。

以上所述仅为本发明示意性的具体实施方式，在不脱离本发明的构思和原则的前提下，任何本领域的技术人员所做出的等同变化与修改，均应属于本发明保护的范围。

Claims (10)

1.一种利用计算统一设备架构CUDA的最佳匹配层PML边界三维地震波传播模拟方法，其特征在于，包括：

读入初始参数和雷克子波；

将三维地震波的第一声波波动方程中的应力分解成相互垂直三个方向的应力分量，并加入PML边界条件形成第二声波波动方程；

运用有限差分方法将所述第二声波波动方程离散化形成有限差分声波波动方程；

给第一三维速度模型加边形成一第二三维速度模型；

根据所述有限差分声波波动方程在中央处理器CPU中计算所述第二三维速度模型的棱、角、面的衰减系数；

根据所述衰减系数在图像处理器GPU中计算所述第二三维速度模型网格点的应力；以及

根据所述应力输出最后时刻的三维地震波模拟记录。

2.如权利要求1所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，还包括：

根据所述衰减系数在图像处理器GPU中计算所述第二三维速度模型网格点的质点运动速度。

3.如权利要求2所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，还包括：

根据所述质点运动速度分析地质构造。

4.如权利要求3所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，所述第一声波波动方程为：

其中，v为网格点的声波速度值；ρ代表介质密度，取值为1；v_x代表x方向的质点运动速度值，v_y代表y方向的质点运动速度值，v_z代表z方向的质点运动速度值，τ代表网格点的应力，t表示地震波声波传播时间；

所述第二声波波动方程为：

τ＝τ_x+τ_y+τ_z,

其中，d(x)是X方向的衰减系数，能够表示为d(x)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；d(y)是Y方向的衰减系数，能够表示为d(y)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；d(z)是Z方向的衰减系数，能够表示为d(z)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；n为指数阶数；v为网格点的声波速度值；ρ代表介质密度，取值为1；t表示地震波声波传播时间；R为反射系数，一般取值为0.0001；v_x为X方向的质点运动速度；v_y为Y方向的质点运动速度；v_z为Z方向的质点运动速度；τ_x为X方向的应力；τ_y为Y方向的应力；τ_z为Z方向的应力；

所述有限差分声波波动方程为：

其中，v_x为X方向的质点运动速度；v_y为Y方向的质点运动速度；v_z为Z方向的质点运动速度；τ_x为X方向的应力；τ_y为Y方向的应力；τ_z为Z方向的应力；ρ代表介质密度，取值为1；a_m为常量；Δt为模拟过程中的时间间隔；Δx为x方向网格间隔；Δy为y方向网格间隔；Δz为z方向网格间隔；i,j,k代表网格点的坐标位置；x,y,z分别代表x轴,y轴和z轴方向；m为计算一个网格点所需要的基准网格点个数；L为计算一个网格点所需要的基准网格点的最大个数；t为时刻，t+1为t+1时刻，即向后推一个时刻，t-1代表t-1时刻，即向前推一个时刻，代表时刻，即后推半个时刻，代表时刻，即向前推半个时刻；d(x)是X方向的衰减系数，能够表示为d(x)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；d(y)是Y方向的衰减系数，能够表示为d(y)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数，P_ML为匹配层的网格厚度；d(z)是Z方向的衰减系数，能够表示为d(z)＝d₀(i/P_ML)ⁿ,i＝0,1,2....，d₀＝(ln(1/R))(n+1)v/(2P_ML)，i代表从区域边界到三维速度模型最内层的网格点数， P_ML为匹配层的网格厚度；n为指数阶数；v为网格点的声波速度值。

5.如权利要求4所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，在CPU中计算所述第二三维速度模型的棱、角、面的衰减系数具体包括：

在X轴方向的两个表面上，仅在X方向上衰减，此时，d(x)≠0，d(y)＝0，d(z)＝0；

在Y轴方向的两个表面上，仅在Y方向上衰减，此时，d(x)＝0，d(y)≠0，d(z)＝0；

在Z轴方向的两个表面上，仅在Z方向上衰减，此时，d(x)＝0，d(y)＝0，d(z)≠0；

与X轴平行的四条棱，仅在X方向上没有衰减，此时，d(x)＝0，d(y)≠0，d(z)≠0；

与Y轴平行的四条棱，仅在Y方向上没有衰减，此时，d(x)≠0，d(y)＝0，d(z)≠0；

与Z轴平行的四条棱，仅在Z方向上没有衰减，此时，d(x)≠0，d(y)≠0，d(z)＝0；

角在X、Y、Z方向上均有衰减，此时，d(x)≠0，d(y)≠0，d(z)≠0。

6.如权利要求3所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，所述GPU的内存类型包括本地存储器和共享存储器，根据所述衰减系数在GPU中计算所述第二三维速度模型网格点的质点运动速度和应力具体包括：

将需要计算的当前网格点的属性值读入本地存储器；

基于所述本地存储器中存储的属性值计算当前网格点的质点运动速度和应力；

所述本地存储器从共享存储器中读入下一网格点与当前网格点不相同的属性值；以及

基于所述本地存储器中当前存储的属性值计算下一网格点的质点运动速度和应力。

7.如权利要求1所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，给所述第一三维速度模型加边形成第二三维速度模型具体包括：

所述第二三维速度模型的面的PML区域由所述第一三维速度模型的表面最上层的网格点的声波速度值映射到PML区域，棱的PML区域由所述第一三维速度模型对应棱的声波速度值映射到PML区域，所述第二三维速度模型的角的PML区域由所述 第一三维速度模型的对应角的声波速度值映射到PML区域。

8.如权利要求1所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，所述初始参数包括炮点位置、检波点位置、时间步长和地震波模拟时长。

9.如权利要求1所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，在CPU中利用OpenMP加速计算所述第二三维速度模型的棱、角、面的所述衰减系数。

10.如权利要求1所述的利用计算统一设备架构CUDA的PML边界三维地震波传播模拟方法，其特征在于，所述三维地震波模拟记录为波场快照图。