CN112666602B - 地震纵波吸收边界条件处理方法及装置 - Google Patents

地震纵波吸收边界条件处理方法及装置 Download PDF

Info

Publication number
CN112666602B
CN112666602B CN201910982447.6A CN201910982447A CN112666602B CN 112666602 B CN112666602 B CN 112666602B CN 201910982447 A CN201910982447 A CN 201910982447A CN 112666602 B CN112666602 B CN 112666602B
Authority
CN
China
Prior art keywords
wave field
field value
moment
value
calculating
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910982447.6A
Other languages
English (en)
Other versions
CN112666602A (zh
Inventor
侯思安
王春明
胡英
崔栋
张才
李萌
张征
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Petrochina Co Ltd
Original Assignee
Petrochina Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Petrochina Co Ltd filed Critical Petrochina Co Ltd
Priority to CN201910982447.6A priority Critical patent/CN112666602B/zh
Publication of CN112666602A publication Critical patent/CN112666602A/zh
Application granted granted Critical
Publication of CN112666602B publication Critical patent/CN112666602B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Geophysics And Detection Of Objects (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请公开了一种地震纵波吸收边界条件处理方法及装置，该方法包括：确定三维速度模型中每一点的波场值的初始值；根据当前时刻和当前时刻的前一时刻的波场值，计算当前时刻的下一时刻的第一波场值；根据第一波场值，利用GPU的寄存器阻塞算法计算当前时刻的下一时刻的单程波波场值；利用第一波场值和单程波波场值计算GPU内存的最慢速方向和第二快速方向的吸收边界条件；根据第一波场值，利用GPU的共享内存阻塞算法计算当前时刻的下一时刻的单程波波场值；利用第一波场值和单程波波场值计算GPU内存的最快速方向的吸收边界条件。本申请可以减少全局内存的读取次数，提高地震纵波波动方程数值模拟的总体运算效率。

Description

地震纵波吸收边界条件处理方法及装置
技术领域
本申请涉及地球物理勘探地震波数值模拟技术领域，尤其涉及一种地震纵波吸收边界条件处理方法及装置。
背景技术
本部分旨在为权利要求书中陈述的本发明实施例提供背景或上下文。此处的描述不因为包括在本部分中就承认是现有技术。
地震纵波波动方程模拟是一种通过数值方法求解波动方程，进而估计纵波在地球介质中传播特征的方法。研究人员和工程师可以通过该方法指导地震勘探中的施工设计，分析采集信号中的有效成分和识别地震剖面中的构造特征。此外，该方法还是许多新方法、新技术，如全波形反演和最小二乘偏移的核心基础。因此，高效的地震纵波波动方程模拟技术有助于提高石油地球物理勘探的效率。
在进行地震纵波波动方程三维数值模拟时，需要同时计算偏微分算子以及处理吸收边界条件，其中吸收边界是为了衰减由于计算机数值边界引起的虚假反射而人为设置的额外计算区域，在真实地质体中是不存在的。
为了提高地震纵波波动方程模拟的效率，许多学者基于图形处理器(GraphicsProcessing Unit，GPU)开展了波动方程微分算子的并行算法研究，并取得了显著的效果。但是针对GPU的吸收边界条件处理算法的研究是非常少的。这是因为边界条件的处理难度远大于微分算子，每个数值计算点的微分算子是完全相同的，容易实现并行计算，而吸收边界条件需要根据不同模型位置分为边、角、面以及上下、前后、左右等关系，加大了并行算法设计的难度，所以，目前的研究重点集中在计算微分算子上，而针对边界条件并行算法的研究很少。
当前已有的针对吸收边界条件的处理方法，往往会遇到下述问题：由于GPU并行计算的特殊性，即全局内存的读取延迟对总的计算时长起到决定性的影响，如果采用的吸收边界条件算法不合理，会造成全局内存的多次读取，严重情况下边界条件的计算时长会占到总时长的50％左右，极大降低了地震纵波波动方程数值模拟的总体运算效率。
发明内容
本申请实施例提供一种地震纵波吸收边界条件处理方法，用以减少全局内存的读取次数，提高地震纵波波动方程数值模拟的总体运算效率，该方法包括：
确定地震纵波传播的三维速度模型在三个维度的吸收边界条件，三维速度模型的三个维度分别表示GPU内存的最慢速方向、第二快速方向和最快速方向；初始化三维速度模型中每一点在三维地震波数值模拟的初始时刻以及初始时刻的前一时刻的波场值；针对三维数值模拟的初始时刻至终止时刻之间的每个时刻，执行下述方法计算吸收边界条件内每个点的波场值：根据当前时刻和当前时刻的前一时刻的波场值，利用声波动方程的数值算法计算当前时刻的下一时刻的第一波场值；根据第一波场值，利用GPU的寄存器阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第二波场值；利用第一波场值和第二波场值计算GPU内存的最慢速方向和第二快速方向的吸收边界条件内每个点的波场值；根据第一波场值，利用GPU的共享内存阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第三波场值；利用第一波场值和第三波场值计算GPU内存的最快速方向的吸收边界条件内每个点的波场值。
本申请实施例还提供一种地震纵波吸收边界条件处理装置，用以减少全局内存的读取次数，提高地震纵波波动方程数值模拟的总体运算效率，该装置包括：
边界条件确定模块，用于确定地震纵波传播的三维速度模型在三个维度的吸收边界条件，三维速度模型的三个维度分别表示GPU内存的最慢速方向、第二快速方向和最快速方向；初始化模块，用于初始化三维速度模型中每一点在三维地震波数值模拟的初始时刻以及初始时刻的前一时刻的波场值；针对三维数值模拟的初始时刻至终止时刻之间的每个时刻，处理模块，用于：根据当前时刻和当前时刻的前一时刻的波场值，利用声波动方程的数值算法计算当前时刻的下一时刻的第一波场值；根据第一波场值，利用GPU的寄存器阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第二波场值；利用第一波场值和第二波场值计算GPU内存的最慢速方向和第二快速方向的吸收边界条件内每个点的波场值；根据第一波场值，利用GPU的共享内存阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第三波场值；利用第一波场值和第三波场值计算GPU内存的最快速方向的吸收边界条件内每个点的波场值。
本申请实施例中，采用GPU的寄存器阻塞算法计算地震波在三维速度模型中最慢速方向和第二快速方向的吸收边界条件，采用GPU的共享内存阻塞算法计算地震波在三维速度模型中最快速度方向的吸收边界条件，通过上述两种GPU的阻塞算法可以减少计算吸收边界条件时全局内存的读取次数，从而减少全局内存的读取延迟，提高地震纵波波动方程数值模拟的总体运算效率。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为本申请实施例中一种地震纵波吸收边界条件处理方法的流程图；
图2为本申请实施例中一种三维速度模型的示意图；
图3为本申请实施例中一种y＝w-1,w-2,···,0的区域示意图；
图4为本申请实施例中一种z＝w-1,w-2,···,0的区域示意图；
图5为本申请实施例中一种地震纵波吸收边界条件处理装置的结构示意图。
具体实施方式
为使本申请实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本申请实施例做进一步详细说明。在此，本申请的示意性实施例及其说明用于解释本申请，但并不作为对本申请的限定。
本申请实施例提供了一种地震纵波吸收边界条件处理方法，如图1所示，该方法包括步骤101至步骤103：
步骤101、确定地震纵波传播的三维速度模型在三个维度的吸收边界条件。
其中，三维速度模型的三个维度分别表示GPU内存的最慢速方向、第二快速方向和最快速方向。
如图2所示，为三维速度模型的示意图。在图2中，y轴方向代表GPU内存的最慢速方向；x轴方向代表GPU内存的第二快速方向；z轴方向代表GPU内存的最快速方向；ny、nx和nz分别表示三维速度模型在三个维度上的大小。
吸收边界区域的宽度为w，具体可以划分为y轴方向吸收边界(y＝w-1,w-2,···,0区域和y＝ny-w-1,ny-w-2,···,ny-1区域)，x轴方向吸收边界(x＝w-1,w-2,···,0区域和x＝nx-w-1,nx-w-2,···,nx-1区域)和z轴方向吸收边界(z＝w-1,w-2,···,0区域和z＝nz-w-1,nz-w-2,···,nz-1区域)。
步骤102、初始化三维速度模型中每一点在三维地震波数值模拟的初始时刻以及初始时刻的前一时刻的波场值。
在本申请实施例中，将初始时刻以及初始时刻的前一时刻的波场值设置为0。
同时，还可以设置两个相邻计算点在空间上的间隔Δs、两个相邻计算点在时间上的间隔Δt、震源点坐标(y0,x0,z0)等关键参数。
需要说明的是，在任意方向上的两个相邻计算点在空间上的间隔都相同，在时间上的间隔也相同。
步骤103、针对三维数值模拟的初始时刻至终止时刻之间的每个时刻，执行下述步骤1031至步骤1033计算吸收边界条件内每个点的波场值。
在本申请实施例中，针对三维数值模拟的初始时刻至终止时刻之间的每个时刻，均计算吸收边界条件内每个点的波场值，也即，在n时刻，计算吸收边界条件内每个点的波场值，在n+1时刻，再计算吸收边界条件内每个点的波场值。
步骤1031、根据当前时刻和当前时刻的前一时刻的波场值，利用声波动方程的数值算法计算当前时刻的下一时刻的第一波场值。
需要说明的是，当前时刻包括初始时刻，当计算初始时刻的下一时刻的第一波场值时，利用的是步骤102中确定的初始时刻以及初始时刻的前一时刻的初始化波场值。当计算得到初始时刻的下一时刻的波场值后，将初始时刻的下一时刻作为当前时刻，初始时刻作为当前时刻的上一时刻，利用初始时刻的下一时刻的波场值及初始时刻的波场值计算当前时刻的下一时刻的波场值，如此循环，直至计算出从初始时刻开始至终止时刻为止，三维速度模型中每个点的第一波场值。
以二阶时间差分和二阶空间差分为例，在本申请实施例中，可以根据如下公式计算当前时刻(n时刻)的下一时刻(n+1时刻)三维速度模型范围内每一点的第一波场值
其中，＝为赋值号，y、x和z分别表示当前计算点在空间坐标系中的位置，y表示GPU内存的最慢速方向，x表示GPU内存的第二快速方向，z表示GPU内存的最快速方向，y＝0,1,2,···,ny-1，x＝0,1,2,···,nx-1，z＝0,1,2,···,nz-1，ny、nx和nz分别用于表示三维速度模型在三个维度上的大小；n表示当前时刻，n-1表示当前时刻的上一时刻，n+1表示当前时刻的下一时刻；Δs表示两个相邻计算点在空间上的间隔；Δt表示两个相邻计算点在时间上的间隔；c表示地震纵波在介质中的传播速度；f表示震源函数。
震源函数可以为现有技术中公开的任意震源函数，示例性的，震源函数可以为雷克子波函数，其函数表达式为：
其中，A表示雷克子波的振幅；ω表示雷克子波的主频；α表示控件衰减因子；y0、x0、z0表示震源中心点坐标。
步骤1032、根据第一波场值，利用GPU的寄存器阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第二波场值；利用第一波场值和第二波场值计算GPU内存的最慢速方向和第二快速方向的吸收边界条件内每个点的波场值。
三维速度模型和波场值在GPU内存上沿着y轴方向和x轴方向都不是连续存储的，无法通过合并访问进行连续读取，但是通过每个GPU线程的寄存器阻塞算法可以实现数据的一次读取多次利用，减少全局内存的读取次数。
在本申请实施例中，可以定义GPU线程块大小为32×16，即每个线程块有512个独立线程，每个线程计算速度模型上的一个网格点，如图3所示，每个大网格表示一个线程块，每个小方块表示一个线程。
需要说明的是，还可以定义其他大小的GPU线程块，例如32×8或32×4等，GPU线程块的大小可以根据读取数据的实际情况确定，对于其具体值，在此不做限定。
具体的，根据当前时刻与当前时刻的前一时刻吸收边界条件中每个点的第一波场值，定义寄存器的初始值；根据寄存器的初始值，确定当前时刻的下一时刻的单程波波场值，作为第二波场值。
其中，定义的寄存器为和/>
针对GPU内存的最慢速方向，寄存器的初始值为 和/>
针对GPU内存的第二快速方向，寄存器的初始值为 和/>
其中，y1∈y，y1＝w-1,w-2,···,0；x1∈x，x1＝w-1,w-2,···,0，w用于表示吸收边界条件的长度。
针对GPU内存的最慢速方向，根据如下公式计算当前时刻的下一时刻的单程波波场值
根据计算GPU内存的最慢速方向的吸收边界条件内点(y1,x,z)的波场值；其中，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs。(y1,x,z)点的第一波场值的计算方法即为步骤1031中/>的计算方法。
针对GPU内存的第二快速方向，根据如下公式计算当前时刻的下一时刻的单程波波场值：
根据计算GPU内存的第二快速方向的吸收边界条件内点(y,x1,z)的波场值；其中，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs。
步骤1033、根据第一波场值，利用GPU的共享内存阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第三波场值；利用第一波场值和第三波场值计算GPU内存的最快速方向的吸收边界条件内每个点的波场值。
因为三维速度模型和波场值在GPU内存上沿着z轴方向是连续存储的，可以通过合并访问进行连续读取，因此通过每个GPU线程块的共享内存阻塞算法可以实现数据的一次读取多次利用，减少全局内存的读取次数。
其中，所使用的GPU线程块的大小为w×1，即每个线程块有w个线程，每个线程计算三维速度模型上的一个网格点。如图4所示，每个网格表示一个线程块，每个小方块表示一个线程。
具体的，根据当前时刻与当前时刻的前一时刻吸收边界条件中每个点的第一波场值，定义共享内存的初始值；根据共享内存的初始值，确定当前时刻的下一时刻的单程波波场值。
其中，定义的共享内存为SC和Sp，共享内存的初始值为
确定共享内存的值之后，根据如下公式计算当前时刻的下一时刻的单程波波场值
其中，z1∈z，z1＝w-1,w-2,···,0，i＝0,1,2,···,w+1，j＝0,1,2,···,w+3，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs，w用于表示吸收边界条件的长度。
之后，根据计算GPU内存的最快速方向的吸收边界条件内每个点(y,x,z1)的波场值/>
由于在步骤103中需要计算从初始时刻至终止时刻之间每个时刻的吸收边界条件内每个点的波场值，而在计算n+1时刻的波场值时需要使用n时刻和n-1时刻的波场值，在计算n+2时刻的波场值时需要使用n时刻和n+1时刻的波场值，因此，可以使用循环迭代的方法来实现步骤1032和步骤1033的计算。具体方法介绍如下：
1、采用8个寄存器的GPU的寄存器阻塞算法计算地震波在y轴方向和x轴方向的吸收边界条件。
由于x轴方向和y轴方向的计算方法相似，因此，在本申请实施例中仅以y轴方向y＝w-1,w-2,···,0为例来对吸收边界条件的计算进行介绍。
第一步，定义GPU线程块大小为32×16，即每个线程块有512个独立线程。
第二步，从y1＝w-1开始计算，定义8个寄存器 每个寄存器的初始值为
第三步，根据单程波计算公式计算n+1时刻的单程波波场值计算公式如下：
第四步，更新吸收边界条件内的波场值，更新公式为：
第五步，通过已有寄存器值更新6个寄存器值通过读取全局内存更新额外两个寄存器值(/>和/>)，更新公式表示为：
第六步，如图3中指示的虚线所示(虚线围成的矩形区域即为由y1＝w-2，x＝0,1,2,···,nx-1以及z＝0,1,2,···,nz-1形成的平面)，重复执行第三步至第五步，直至y1＝w-1,w-2,···,0的吸收边界条件全部计算完毕。
2、采用GPU的共享内存阻塞算法计算地震波在z轴方向的吸收边界条件。
以z1＝w-1,w-2,···,0区域为例：
第一步，定义GPU线程块的大小为w×1，即每个线程块有w个线程。
第二步，如图4中实线(图4中以“实线”标示)区域所示，为y＝0，z1＝w-1,w-2,···,0以及x＝0,1,2,...,ny-1构成的矩形区域，下面将从y＝0开始计算，定义2个共享内存(SC和Sp)，来计算该计算区域中每个点的波场值。每个共享内存的初始化值为：
第三步，根据单程波计算公式计算n+1时刻的单程波波场值计算公式如下：
第四步，更新吸收边界条件内的波场值，更新公式为：
第五步，更新y＝y+1的共享内存(SC和Sp)：
第六步，如图4中的虚线所示，重复执行第三步至第五步，直至y＝0,1,2,···,ny-1的吸收边界条件全部计算完毕。
本申请实施例中，采用GPU的寄存器阻塞算法计算地震波在三维速度模型中最慢速方向和第二快速方向的吸收边界条件，采用GPU的共享内存阻塞算法计算地震波在三维速度模型中最快速度方向的吸收边界条件，通过上述两种GPU的阻塞算法可以减少计算吸收边界条件时全局内存的读取次数，从而减少全局内存的读取延迟，提高地震纵波波动方程数值模拟的总体运算效率。
本申请实施例提供了一种地震纵波吸收边界条件处理装置，如图5所示，该装置500包括边界条件确定模块501、初始化模块502和处理模块503。
其中，边界条件确定模块501，用于确定地震纵波传播的三维速度模型在三个维度的吸收边界条件，三维速度模型的三个维度分别表示GPU内存的最慢速方向、第二快速方向和最快速方向。
初始化模块502，用于初始化三维速度模型中每一点在三维地震波数值模拟的初始时刻以及初始时刻的前一时刻的波场值。
针对三维数值模拟的初始时刻至终止时刻之间的每个时刻，处理模块503，用于：
根据当前时刻和当前时刻的前一时刻的波场值，利用声波动方程的数值算法计算当前时刻的下一时刻的第一波场值；
根据第一波场值，利用GPU的寄存器阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第二波场值；利用第一波场值和第二波场值计算GPU内存的最慢速方向和第二快速方向的吸收边界条件内每个点的波场值；
根据第一波场值，利用GPU的共享内存阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第三波场值；利用第一波场值和第三波场值计算GPU内存的最快速方向的吸收边界条件内每个点的波场值。
在本申请实施例的一种实现方式中，处理模块503，用于：
将初始时刻作为当前时刻，计算当前时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值；
在初始时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值计算完毕后，将初始时刻的下一时刻作为当前时刻，重新计算当前时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值。
在本申请实施例的一种实现方式中，处理模块503，用于：
根据计算当前时刻的下一时刻的第一波场值/>
其中，＝为赋值号，y、x和z分别表示当前计算点在空间坐标系中的位置，y表示GPU内存的最慢速方向，x表示GPU内存的第二快速方向，z表示GPU内存的最快速方向，y＝0,1,2,···,ny-1，x＝0,1,2,···,nx-1，z＝0,1,2,···,nz-1，ny、nx和nz分别用于表示三维速度模型在三个维度上的大小；n表示当前时刻，n-1表示当前时刻的上一时刻，n+1表示当前时刻的下一时刻；Δs表示两个相邻计算点在空间上的间隔；Δt表示两个相邻计算点在时间上的间隔；c表示地震纵波在介质中的传播速度；f表示震源函数。
在本申请实施例的一种实现方式中，处理模块503，用于：
根据当前时刻与当前时刻的前一时刻吸收边界条件中每个点的第一波场值，定义寄存器的初始值；
根据寄存器的初始值，确定当前时刻的下一时刻的单程波波场值。
在本申请实施例的一种实现方式中，定义的寄存器为和/>针对GPU内存的最慢速方向，寄存器的初始值为/> 和/>针对GPU内存的第二快速方向，寄存器的初始值为/> 和/>其中，y1∈y，y1＝w-1,w-2,···,0；x1∈x，x1＝w-1,w-2,···,0，w用于表示吸收边界条件的长度。
在本申请实施例的一种实现方式中，处理模块503，用于：
根据计算当前时刻的下一时刻的单程波波场值/>
利用第一波场值和第二波场值计算GPU内存的最慢速方向的吸收边界条件内每个点的波场值，包括：
根据计算GPU内存的最慢速方向的吸收边界条件内点(y1,x,z)的波场值；
其中，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs。
在本申请实施例的一种实现方式中，处理模块503，用于：
根据计算当前时刻的下一时刻的单程波波场值/>
利用第一波场值和第二波场值计算GPU内存的最慢速方向的吸收边界条件内每个点的波场值，包括：
根据计算GPU内存的第二快速方向的吸收边界条件内点(y,x1,z)的波场值；
其中，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs。
在本申请实施例的一种实现方式中，处理模块503，用于：
根据当前时刻与当前时刻的前一时刻吸收边界条件中每个点的第一波场值，定义共享内存的初始值；
根据共享内存的初始值，确定当前时刻的下一时刻的单程波波场值。
在本申请实施例的一种实现方式中，定义的共享内存为SC和Sp，共享内存的初始值为
处理模块503，用于：
根据计算当前时刻的下一时刻的单程波波场值/>
其中，z1∈z，z1＝w-1,w-2,···,0，i＝0,1,2,···,w+1，j＝0,1,2,···,w+3，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs，w用于表示吸收边界条件的长度。
在本申请实施例的一种实现方式中，处理模块503，用于：
根据计算GPU内存的最快速方向的吸收边界条件内每个点(y,x,z1)的波场值/>
本申请实施例中，采用GPU的寄存器阻塞算法计算地震波在三维速度模型中最慢速方向和第二快速方向的吸收边界条件，采用GPU的共享内存阻塞算法计算地震波在三维速度模型中最快速度方向的吸收边界条件，通过上述两种GPU的阻塞算法可以减少计算吸收边界条件时全局内存的读取次数，从而减少全局内存的读取延迟，提高地震纵波波动方程数值模拟的总体运算效率。
本申请实施例还提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行计算机程序时实现权利要求步骤101至步骤103任一方法。
本申请实施例还提供了一种计算机可读存储介质，计算机可读存储介质存储有执行步骤101至步骤103任一方法的计算机程序。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例，对本申请的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本申请的具体实施例而已，并不用于限定本申请的保护范围，凡在本申请的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (6)

1.一种地震纵波吸收边界条件处理方法，其特征在于，所述方法包括：
分别确定地震纵波传播的三维速度模型在三个维度的吸收边界条件，三维速度模型的三个维度分别表示GPU内存的最慢速方向、第二快速方向和最快速方向；
初始化三维速度模型中每一点在三维地震波数值模拟的初始时刻以及初始时刻的前一时刻的波场值；
针对三维数值模拟的初始时刻至终止时刻之间的每个时刻，执行下述步骤计算吸收边界条件内每个点的波场值：
根据当前时刻和当前时刻的前一时刻的波场值，利用声波动方程的数值算法计算当前时刻的下一时刻的第一波场值；
根据第一波场值，利用GPU的寄存器阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第二波场值；利用第一波场值和第二波场值计算GPU内存的最慢速方向和第二快速方向的吸收边界条件内每个点的波场值；
根据第一波场值，利用GPU的共享内存阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第三波场值；利用第一波场值和第三波场值计算GPU内存的最快速方向的吸收边界条件内每个点的波场值；
所述利用声波动方程的数值算法计算当前时刻的下一时刻的第一波场值，包括：
根据计算当前时刻的下一时刻的第一波场值/>
其中，＝为赋值号，y、x和z分别表示当前计算点在空间坐标系中的位置，y表示GPU内存的最慢速方向，x表示GPU内存的第二快速方向，z表示GPU内存的最快速方向，y＝0,1,2,···,ny-1，x＝0,1,2,···,nx-1，z＝0,1,2,···,nz-1，ny、nx和nz分别用于表示三维速度模型在三个维度上的大小；n表示当前时刻，n-1表示当前时刻的上一时刻，n+1表示当前时刻的下一时刻；Δs表示两个相邻计算点在空间上的间隔；Δt表示两个相邻计算点在时间上的间隔；c表示地震纵波在介质中的传播速度；f表示震源函数；
所述利用GPU的寄存器阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第二波场值，包括：
根据当前时刻与当前时刻的前一时刻吸收边界条件中每个点的第一波场值，定义寄存器的初始值；
根据寄存器的初始值，确定当前时刻的下一时刻的单程波波场值；
定义的寄存器为和/>针对GPU内存的最慢速方向，所述寄存器的初始值为/> 和/>针对GPU内存的第二快速方向，所述寄存器的初始值为/> 和/>其中，y1∈y，y1＝w-1,w-2,···,0；x1∈x，x1＝w-1,w-2,···,0，w用于表示吸收边界条件的长度；
针对GPU内存的最慢速方向，所述根据寄存器的初始值，确定当前时刻的下一时刻的单程波波场值，包括：
根据计算当前时刻的下一时刻的单程波波场值/>
利用第一波场值和第二波场值计算GPU内存的最慢速方向的吸收边界条件内每个点的波场值，包括：
根据计算GPU内存的最慢速方向的吸收边界条件内点(y1,x,z)的波场值；
其中，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs；
针对GPU内存的第二快速方向，所述根据寄存器的初始值，确定当前时刻的下一时刻的单程波波场值，包括：
根据计算当前时刻的下一时刻的单程波波场值/>
利用第一波场值和第二波场值计算GPU内存的第二快速方向的吸收边界条件内每个点的波场值，包括：
根据计算GPU内存的第二快速方向的吸收边界条件内点(y,x1,z)的波场值；
其中，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs；
根据第一波场值，利用GPU的共享内存阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第三波场值，包括：
根据当前时刻与当前时刻的前一时刻吸收边界条件中每个点的第一波场值，定义共享内存的初始值；
根据共享内存的初始值，确定当前时刻的下一时刻的单程波波场值；
定义的共享内存为和/>共享内存的初始值为/>
根据共享内存的初始值，确定当前时刻的下一时刻的单程波波场值，包括：
根据计算当前时刻的下一时刻的单程波波场值/>
其中，z1∈z，z1＝w-1,w-2,···,0，i＝0,1,2,···,w+1，j＝0,1,2,···,w+3，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs，w用于表示吸收边界条件的长度；
利用第一波场值和第三波场值计算GPU内存的最快速方向的吸收边界条件内每个点的波场值，包括：
根据计算GPU内存的最快速方向的吸收边界条件内每个点(y,x,z1)的波场值/>
2.根据权利要求1所述的方法，其特征在于，在初始化三维速度模型中的每一点在三维地震波数值模拟的初始时刻及初始时刻的前一时刻的波场值之后，所述方法包括：
将初始时刻作为当前时刻，计算当前时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值；
在初始时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值计算完毕后，将初始时刻的下一时刻作为当前时刻，重新计算当前时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值。
3.一种地震纵波吸收边界条件处理装置，其特征在于，所述装置包括：
边界条件确定模块，用于确定地震纵波传播的三维速度模型在三个维度的吸收边界条件，三维速度模型的三个维度分别表示GPU内存的最慢速方向、第二快速方向和最快速方向；
初始化模块，用于初始化三维速度模型中每一点在三维地震波数值模拟的初始时刻以及初始时刻的前一时刻的波场值；
针对三维数值模拟的初始时刻至终止时刻之间的每个时刻，处理模块，用于：
根据当前时刻和当前时刻的前一时刻的波场值，利用声波动方程的数值算法计算当前时刻的下一时刻的第一波场值；
根据第一波场值，利用GPU的寄存器阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第二波场值；利用第一波场值和第二波场值计算GPU内存的最慢速方向和第二快速方向的吸收边界条件内每个点的波场值；
根据第一波场值，利用GPU的共享内存阻塞算法计算当前时刻的下一时刻的单程波波场值，作为第三波场值；利用第一波场值和第三波场值计算GPU内存的最快速方向的吸收边界条件内每个点的波场值；
所述处理模块，用于：
根据计算当前时刻的下一时刻的第一波场值/>
其中，＝为赋值号，y、x和z分别表示当前计算点在空间坐标系中的位置，y表示GPU内存的最慢速方向，x表示GPU内存的第二快速方向，z表示GPU内存的最快速方向，y＝0,1,2,···,ny-1，x＝0,1,2,···,nx-1，z＝0,1,2,···,nz-1，ny、nx和nz分别用于表示三维速度模型在三个维度上的大小；n表示当前时刻，n-1表示当前时刻的上一时刻，n+1表示当前时刻的下一时刻；Δs表示两个相邻计算点在空间上的间隔；Δt表示两个相邻计算点在时间上的间隔；c表示地震纵波在介质中的传播速度；f表示震源函数；
所述处理模块，用于：
根据当前时刻与当前时刻的前一时刻吸收边界条件中每个点的第一波场值，定义寄存器的初始值；
根据寄存器的初始值，确定当前时刻的下一时刻的单程波波场值；
定义的寄存器为和/>针对GPU内存的最慢速方向，所述寄存器的初始值为/> 和/>针对GPU内存的第二快速方向，所述寄存器的初始值为/> 和/>其中，y1∈y，y1＝w-1,w-2,···,0；x1∈x，x1＝w-1,w-2,···,0，w用于表示吸收边界条件的长度；
所述处理模块，用于：
根据计算当前时刻的下一时刻的单程波波场值/>
利用第一波场值和第二波场值计算GPU内存的最慢速方向的吸收边界条件内每个点的波场值，包括：
根据计算GPU内存的最慢速方向的吸收边界条件内点(y1,x,z)的波场值；
其中，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs；
所述处理模块，用于：
根据计算当前时刻的下一时刻的单程波波场值/>
利用第一波场值和第二波场值计算GPU内存的第二快速方向的吸收边界条件内每个点的波场值，包括：
根据计算GPU内存的第二快速方向的吸收边界条件内点(y,x1,z)的波场值；
其中，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs；
所述处理模块，用于：
根据当前时刻与当前时刻的前一时刻吸收边界条件中每个点的第一波场值，定义共享内存的初始值；
根据共享内存的初始值，确定当前时刻的下一时刻的单程波波场值；
定义的共享内存为和/>共享内存的初始值为/>
所述处理模块，用于：
根据计算当前时刻的下一时刻的单程波波场值/>
其中，z1∈z，z1＝w-1,w-2,···,0，i＝0,1,2,···,w+1，j＝0,1,2,···,w+3，T11＝(2-r)(1-r)/2，T12＝r(2-r)，T13＝r(r-1)/2，r＝cΔt/Δs，w用于表示吸收边界条件的长度；
所述处理模块，用于：
根据计算GPU内存的最快速方向的吸收边界条件内每个点(y,x,z1)的波场值/>
4.根据权利要求3所述的装置，其特征在于，在初始化模块初始化三维速度模型中的每一点在三维地震波数值模拟的初始时刻及初始时刻的前一时刻的波场值之后，所述处理模块，用于：
将初始时刻作为当前时刻，计算当前时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值；
在初始时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值计算完毕后，将初始时刻的下一时刻作为当前时刻，重新计算当前时刻的下一时刻三维速度模型的吸收边界条件内每个点的波场值。
5.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至2任一所述方法。
6.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有执行权利要求1至2任一所述方法的计算机程序。
CN201910982447.6A 2019-10-16 2019-10-16 地震纵波吸收边界条件处理方法及装置 Active CN112666602B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910982447.6A CN112666602B (zh) 2019-10-16 2019-10-16 地震纵波吸收边界条件处理方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910982447.6A CN112666602B (zh) 2019-10-16 2019-10-16 地震纵波吸收边界条件处理方法及装置

Publications (2)

Publication Number Publication Date
CN112666602A CN112666602A (zh) 2021-04-16
CN112666602B true CN112666602B (zh) 2023-10-31

Family

ID=75400621

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910982447.6A Active CN112666602B (zh) 2019-10-16 2019-10-16 地震纵波吸收边界条件处理方法及装置

Country Status (1)

Country Link
CN (1) CN112666602B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103616721A (zh) * 2013-11-25 2014-03-05 中国石油天然气股份有限公司 基于二阶偏微分波动方程的pml吸收边界条件的方法
CN109725351A (zh) * 2018-12-18 2019-05-07 中国石油天然气集团有限公司 一种3d弹性波混合吸收边界条件的确定方法、装置及系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103616721A (zh) * 2013-11-25 2014-03-05 中国石油天然气股份有限公司 基于二阶偏微分波动方程的pml吸收边界条件的方法
CN109725351A (zh) * 2018-12-18 2019-05-07 中国石油天然气集团有限公司 一种3d弹性波混合吸收边界条件的确定方法、装置及系统

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
A GPU-accelerated nodal discontinuous Galerkin method with high-order absorbing boundary conditions and corner/edge compatibility;A. Modave et al.;Int J Numer Meth Engng.;第112卷;第1659-1686页 *
Accelerating a three-dimensional finite-difference wave propagation code using GPU graphics cards;David Michea et al.;Geophysical Journal International;第182卷;第389-402页 *
一种新的弹性波数值模拟吸收边界条件;李文杰;魏修成;宁俊瑞;;石油地球物理勘探(第04期);第501-507页 *
基于CPU/GPU协同加速叠前逆时偏移方法研究;高新成;李春生;;陕西理工学院学报(自然科学版)(第01期);第44-49页 *
声波方程逆时偏移中的无分裂PML吸收边界条件;王鹏飞;何兵寿;;工程地球物理学报(第05期);第583-590页 *
组合吸收边界条件下VTI介质地震波场模拟;杜启振;李宾;;石油地球物理勘探(第02期);第173-178页 *

Also Published As

Publication number Publication date
CN112666602A (zh) 2021-04-16

Similar Documents

Publication Publication Date Title
Huthwaite Accelerated finite element elastodynamic simulations using the GPU
Ingram et al. Developments in Cartesian cut cell methods
Marrone et al. Study of ship wave breaking patterns using 3D parallel SPH simulations
Komatitsch et al. High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster
CN107066720B (zh) 一种基于piv技术的可压缩流体压力场的计算方法和装置
CN109709602B (zh) 一种远探测声波偏移成像方法、装置及系统
Giuliani et al. Adaptive mesh refinement on graphics processing units for applications in gas dynamics
CN112666602B (zh) 地震纵波吸收边界条件处理方法及装置
CN109490948B (zh) 地震声学波动方程矢量并行计算方法
CN112859162B (zh) 基于伪谱法的三维声波方程空间域正演方法及装置
Weinbub et al. Shared-memory parallelization of the fast marching method using an overlapping domain-decomposition approach
CN115327626A (zh) 一种用于三维ati介质的弹性波场矢量分解方法及系统
CN113866827A (zh) 一种解释性速度建模地震成像方法、系统、介质和设备
CN110162804B (zh) 基于cpu加速的波场正演模拟优化方法
Tang et al. A fast RTM implementation in TTI media
CN106353801A (zh) 三维Laplace域声波方程数值模拟方法及装置
CN112836327B (zh) 三维波动方程有限差分数值模拟方法及装置
CN114624766B (zh) 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法
Moreira et al. An Architecture Using a Finite Difference Method to Calculate Realistic Sound Equalization in Games
CN105372704B (zh) 一种获取地震波传播方向的方法及装置
Calandra et al. Recent advances in numerical methods for solving the wave equation in the context of seismic depth imaging
CN114924313B (zh) 基于行波分离的声波最小二乘逆时偏移梯度求取方法
Liu et al. Using pseudo-spectral method on curved grids for SH-wave modeling of irregular free-surface
CN113126151B (zh) 一种基于纯波延拓方程的弹性反射波旅行时反演方法
CN114879249B (zh) 基于四面体单元走时扰动插值的地震波前走时计算方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant