CN105445786A - 用于基于gpu获取叠前逆时偏移的方法及装置 - Google Patents

Abstract

本发明公开了一种用于基于GPU来获取叠前逆时偏移的方法，其包括以下步骤：通过震源点正向外推过程得到各个时刻的震源点正向外推波场数据；通过接收点逆向外推过程得到各个时刻的接收点逆向外推波场数据；将同一时刻下的所述震源点正向外推波场数据与所述接收点逆向外推波场数据进行相关求和得到深度域成像值。在接收点逆向外推过程和震源点正向外推过程中仅采用两个存储空间来存储当前时刻的波场数据和当前计算出的波场数据。本发明所设计的基于GPU的叠前逆时偏移方法，在不改变计算结果及不降低计算效率的前提下，将GPU数据空间需求减少近1/3，提高了GPU空间的利用率，使得该方法更加适应于大规模数据处理。

Description

用于基于GPU获取叠前逆时偏移的方法及装置

技术领域

本发明涉及地质勘探领域，具体说涉及一种用于获取叠前逆时偏移的方法及装置。

背景技术

近年来，随着计算机处理能力的发展以及一些新的并行设备的出现，特别是GPU(图形处理器)的通用计算能力和可编程性的提高，为加速科学计算带来了新的选择和希望。相比CPU(中央处理器)，GPU拥有天生的并行结构，而且功耗低、成本小、管理简单。目前，由NVIDIA公司推出的基于自身GPU的运算平台CUDA(ComputeUnifiedDeviceArchitecture，统一计算设备架构)最引人瞩目。在现有的技术中，已实现了基于CUDA的图像与视频处理、计算生物学和化学、流体力学模拟、CT图像再现、以及光线追踪等。

在勘探领域中，叠前逆时偏移RTM(ReverseTimeMigration)是对单炮记录数据进行逆时偏移，然后将各炮成像结果叠加，得到最终的成像剖面。对单炮记录，首先对震源波场进行正向外推，并存储每个时间的波场边界；然后将炮记录的最后一个采样时刻的波场(x,z,T)作为起始平面，按时间反推，并以存储的当前时间的波场边界及边界条件得到此时的震源波场，再应用互相关成像条件得到偏移结果。

每一步波场外推都是利用三维有限差分算法计算的，计算量巨大；同时需要大量的存储空间记录每一个时间的震源外推波场。所以，RTM是一个计算密集型和数据密集型的地震偏移成像方法。

由于逆时偏移的实现需要庞大的计算量，而且计算过程中会产生低频噪音，这些问题阻碍了逆时偏移在实际生产中的应用。即使引入近年来快速发展的GPU及其编程架构CUDA，但由于GPU的显存空间较小，一般为4-8GB。因此，对于大数据、大运算量的RTM计算是其仍然存在制约。

因此，在实际的GPU应用中，迫切需要一种能够节约存储空间又不降低计算效率的叠前逆时偏移计算的实现方法。

发明内容

为解决上述问题，本发明提供了一种基于GPU获取叠前逆时偏移的方法，其包括以下步骤：

通过震源点正向外推过程得到各个时刻的震源点正向外推波场数据；

通过接收点逆向外推过程得到各个时刻的接收点逆向外推波场数据；

将同一时刻下的所述震源点正向外推波场数据与所述接收点逆向外推波场数据进行相关求和得到深度域成像值；

其中，所述震源点正向外推过程包括以下步骤：

S101、模拟震源点的前一时刻和当前时刻的波场数据；

S102、分配两个存储空间分别存储所述震源点的前一时刻和当前时刻的波场数据；

S103、基于两个存储空间中所存储的所述震源点的前一时刻和当前时刻的波场数据根据全波波动方程计算出下一时刻的波场数据，并将所述下一时刻的波场数据保存到用于存储所述前一时刻的波场数据的存储空间中；

S104、判断当前时刻是否为最终时刻，如果是，则输出用于存储所述前一时刻的波场数据的存储空间中的数据，如果否，继续回到步骤S103。

根据本发明的一个实施例，所述接收点逆向外推过程包括以下步骤：

S201、获取接收点的下一时刻和当前时刻的波场数据；

S202、分配两个存储空间分别存储所述接收点的下一时刻和当前时刻的波场数据；

S203、基于两个存储空间中所存储的所述接收点的下一时刻和当前时刻的波场数据根据全波波动方程计算出前一时刻的波场数据，并将所述前一时刻的波场数据保存到用于存储所述下一时刻的波场数据的存储空间中；

S204、判断当前时刻是否为起始时刻，如果是，则输出用于存储所述下一时刻的波场数据的存储空间中的数据，如果否，继续回到步骤S203。

根据本发明的一个实施例，所述全波波动方程为：

$\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2}+\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{v^2(x,y,z)}\frac{{\∂}^{2}u}{{\∂t}^2}$

其中，u(x,y,z,t)为地表记录的压力波场；v(x,y,z)为纵横向可变的介质速度。

根据本发明的一个实施例，该方法还包括对压力波场进行变形处理，并代入所述全波波动方程中，得到截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维正演模拟高阶差分方程为：

$\begin{array}{c}{u}_{i,j,k}^{n+1}={2u}_{i,j,k}^n-u_{i,j,k}^{n-1}+\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δx}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i+m,j,k}^n+u_{i-m,j,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δy}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j+m,k}^n+u_{i,j-m,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δz}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j,k+m}^n+u_{i,j,k-m}^n)]\end{array}$

在步骤S103中，基于所述截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维正演模拟高阶差分方程来计算出震源点的下一时刻的波场数据。

根据本发明的一个实施例，对压力波场进行变形处理并代入所述全波波动方程中，得到截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维逆时深度偏移高阶差分程为：

$\begin{array}{c}{u}_{i,j,k}^{n-1}={2u}_{i,j,k}^n-u_{i,j,k}^{n+1}+\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δx}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i+m,j,k}^n+u_{i-m,j,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δy}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j+m,k}^n+u_{i,j-m,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δz}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j,k+m}^n+u_{i,j,k-m}^n)]\end{array}$

在步骤S203中，基于所述截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维逆时偏移高阶差分方程来计算出接收点的前一时刻的波场数据。M为有限差分的阶数。

在本发明的另一个方面中，其还提供了一种用于基于GPU来获取叠前逆时偏移的装置，该装置包括：

震源点正向外推模块，其通过震源点正向外推过程得到各个时刻的震源点正向外推波场数据；

接收点逆向外推模块，其通过接收点逆向外推过程得到各个时刻的接收点逆向外推波场数据；

深度域成像值计算模块，其用于将同一时刻下的所述震源点正向外推波场数据与所述接收点逆向外推波场数据进行相关求和得到深度域成像值；

其中，所述震源点正向外推模块包括以下单元：

模拟单元，其用于模拟震源点的前一时刻和当前时刻的波场数据；

第一存储空间分配单元，其用于分配两个存储空间以分别存储所述震源点的前一时刻和当前时刻的波场数据；

第一波场数据计算单元，其用于基于两个存储空间中所存储的所述震源点的前一时刻和当前时刻的波场数据根据全波波动方程计算出下一时刻的波场数据，并将所述下一时刻的波场数据保存到用于存储所述前一时刻的波场数据的存储空间中；

第一判断单元，其用于判断当前时刻是否为最终时刻，如果是，则输出用于存储所述前一时刻的波场数据的存储空间中的数据，如果否，继续到所述第一波场数据计算单元中计算波场数据。

根据本发明的一个实施例，所述接收点逆向外推模块包括以下单元：

获取单元，其用于获取接收点的下一时刻和当前时刻的波场数据；

第二存储空间分配单元，其用于分配两个存储空间分别存储所述接收点的下一时刻和当前时刻的波场数据；

第二波场数据计算单元，其用于基于两个存储空间中所存储的所述接收点的下一时刻和当前时刻的波场数据根据全波波动方程计算出前一时刻的波场数据，并将所述前一时刻的波场数据保存到用于存储所述下一时刻的波场数据的存储空间中；

第二判断单元，其用于判断当前时刻是否为起始时刻，如果是，则输出用于存储所述下一时刻的波场数据的存储空间中的数据，如果否，继续到所述第二波场数据计算单元中计算波场数据。

本发明带来了以下有益效果：

根据本发明的方法，在不改变计算结果以及不降低计算效率的前提下，将GPU数据空间需求减少近1/3。因此，本发明提高了GPU空间的利用率，使得GPU技术更适应于大规模数据处理。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要的附图做简单的介绍：

图1是根据本发明的方法计算叠前逆时偏移并进行成像的方法总体流程图；

图2是根据本发明的一个实施例进行震源点正向外推波场数据计算的方法流程图；以及

图3是是根据本发明的一个实施例进行接收点逆向外推波场数据计算的方法流程图。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

另外，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

地震波叠前逆时偏移RTM在处理多波多分量地震资料方面比常规偏移方法有很大的优势。叠前逆时偏移是基于双程波动方程的偏移方法，因此在偏移过程中不会受到地层倾角和介质横向速度变化的限制，同时逆时偏移在计算中保留的波场的矢量特征，可以得到更高精度的成像。

逆时偏移方法是通过双程波波动方程在时间上对地震资料进行反向外推并结合成像条件实现偏移的。如图1所示，逆时偏移在算法实现上主要包括三部分内容：震源点波场的正向延拓；接收点波场的反向延拓；恰当的成像条件。震源点波场和接收点波场都根据以下方程进行双程外推：

$\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2}+\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{v^2(x,y,z)}\frac{{\∂}^{2}u}{{\∂t}^2}---\left(1\right)$

其中：u(x,y,z,t)为地表记录的压力波场；v(x,y,z)为纵横向可变的介质速度。

首先，对压力波场进行泰勒Taylor展开：

$u(t+\mathrm{\Δt})=u\left(t\right)+\frac{\∂u}{\∂t}\mathrm{\Δt}+\frac{1}{2!}\frac{{\∂}^{2}u}{{\∂t}^2}{\left(\mathrm{\Δt}\right)}^2+\frac{1}{3!}\frac{{\∂}^{3}u}{{\∂t}^3}{\left(\mathrm{\Δt}\right)}^3+\frac{1}{4!}\frac{{\∂}^{4}u}{{\∂t}^4}{\left(\mathrm{\Δt}\right)}^4+...\left(2\right)$

$u(t-\mathrm{\Δt})=u\left(t\right)-\frac{\∂u}{\∂t}\mathrm{\Δt}+\frac{1}{2!}\frac{{\∂}^{2}u}{{\∂t}^2}{\left(\mathrm{\Δt}\right)}^2-\frac{1}{3!}\frac{{\∂}^{3}u}{{\∂t}^3}{\left(\mathrm{\Δt}\right)}^3+\frac{1}{4!}\frac{{\∂}^{4}u}{{\∂t}^4}{\left(\mathrm{\Δt}\right)}^4-...\left(3\right)$

由式(2)和式(3)相加，可得

$\frac{{\∂}^{2}u}{{\∂t}^2}=\frac{1}{{\mathrm{\Δt}}^2}\left\{\right[u(t+\mathrm{\Δt})-2u\left(t\right)+u(t-\mathrm{\Δt})]-\frac{2}{4!}\frac{{\∂}^{4}u}{{\∂t}^4}{\left(\mathrm{\Δt}\right)}^4+...\}---\left(4\right)$

在利用式(4)进行正演模拟或偏移成像过程中，差分方程所涉及的时间层越多，所需的内存也就越大。

为避开此问题，利用声波方程(4)把对时间的高阶微分转移到空间微分上去。因此，有：

$\begin{array}{c}\frac{{\∂}^{4}u}{{\∂t}^4}=\frac{\∂}{{\∂t}^2}\left(\frac{{\∂}^{2}u}{{\∂t}^2}\right)=\frac{\∂}{{\∂t}^2}\left[v^2(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2}+\frac{{\∂}^{2}u}{{\∂z}^2})\right]\\ =v^2[\frac{\∂}{{\∂x}^2}\left(\frac{{\∂}^{2}u}{{\∂t}^2}\right)+\frac{\∂}{{\∂y}^2}\left(\frac{{\∂}^{2}u}{{\∂t}^2}\right)+\frac{\∂}{\∂z^2}\left(\frac{{\∂}^{2}u}{{\∂t}^2}\right)]\\ =v^4(\frac{{\∂}^{4}u}{{\∂x}^4}+\frac{{\∂}^{4}u}{{\∂y}^4}+\frac{{\∂}^{4}u}{{\∂z}^4})+2v^4(\frac{{\∂}^{4}u}{{\∂x}^2{\∂y}^2}+\frac{{\∂}^{4}u}{{\∂y}^2{\∂z}^2}+\frac{{\∂}^{4}u}{{\∂z}^2{\∂x}^2})\end{array}---\left(5\right)$

把式(4)和式(5)式代入式(1)得：

$\begin{array}{c}u(t+\mathrm{\Δt})=2u\left(t\right)-u(t-\mathrm{\Δt})+{\left(\mathrm{v\Δt}\right)}^2(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2}+\frac{{\∂}^{2}u}{{\∂z}^2})+\frac{2}{4!}{\left(\mathrm{v\Δt}\right)}^4(\frac{{\∂}^{4}u}{{\∂x}^4}+\frac{{\∂}^{4}u}{\∂y^4}+\frac{{\∂}^{4}u}{{\∂z}^4})\\ +\frac{4}{4!}{\left(\mathrm{v\Δt}\right)}^4(\frac{{\∂}^{4}u}{{\∂x}^2{\∂y}^2}+\frac{{\∂}^{4}u}{{\∂y}^2{\∂z}^2}+\frac{{\∂}^{4}u}{{\∂z}^2{\∂x}^2})+O\left({\mathrm{\Δt}}^4\right)\end{array}---\left(6\right)$

$\begin{array}{c}u(t-\mathrm{\Δt})=2u\left(t\right)-u(t+\mathrm{\Δt})+{\left(\mathrm{v\Δt}\right)}^2(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2}+\frac{{\∂}^{2}u}{{\∂z}^2})+\frac{2}{4!}{\left(\mathrm{v\Δt}\right)}^4(\frac{{\∂}^{4}u}{{\∂x}^4}+\frac{{\∂}^{4}u}{\∂y^4}+\frac{{\∂}^{4}u}{{\∂z}^4})\\ +\frac{4}{4!}{\left(\mathrm{v\Δt}\right)}^4(\frac{{\∂}^{4}u}{{\∂x}^2{\∂y}^2}+\frac{{\∂}^{4}u}{{\∂y}^2{\∂z}^2}+\frac{{\∂}^{4}u}{{\∂z}^2{\∂x}^2})+O\left({\mathrm{\Δt}}^4\right)\end{array}---\left(7\right)$

(6)式和(7)式分别是推导用于正演模拟和逆时偏移RTM的高阶差分方程的起始方程。它在时间方向上截断误差为O(Δt⁴)，而空间微商的差商的截断误差根据需要而定。

波场对空间的高阶近似为：

$2\frac{{\∂}^{2}u}{{\∂x}^2}{\mathrm{\Δx}}^2={\mathrm{\ω}}_{0}u\left(x\right)+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m[u(x+\mathrm{m\Δx})+u(x-\mathrm{m\Δx})]+O\left({\mathrm{\Δx}}^M\right)---\left(8\right)$

$2\frac{{\∂}^{2}u}{{\∂y}^2}{\mathrm{\Δy}}^2={\mathrm{\ω}}_{0}u\left(y\right)+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m[u(y+\mathrm{m\Δy})+u(y-\mathrm{m\Δy})]+O\left({\mathrm{\Δy}}^M\right)---\left(9\right)$

$2\frac{{\∂}^{2}u}{{\∂z}^2}{\mathrm{\Δz}}^2={\mathrm{\ω}}_{0}u\left(z\right)+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m[u(z+\mathrm{m\Δz})+u(z-\mathrm{m\Δz})]+O\left({\mathrm{\Δz}}^M\right)---\left(10\right)$

把式(8)、(9)和(10)代入式(6)或式(7)，可以写出具有任意截断误差的高阶差分方程。

截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维正演模拟高阶差分方程为：

$\begin{array}{c}{u}_{i,j,k}^{n+1}={2u}_{i,j,k}^n-u_{i,j,k}^{n-1}+\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δx}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i+m,j,k}^n+u_{i-m,j,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δy}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j+m,k}^n+u_{i,j-m,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δz}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j,k+m}^n+u_{i,j,k-m}^n)]\end{array}---\left(11\right)$

截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维逆时深度偏移高阶差分程为：

$\begin{array}{c}{u}_{i,j,k}^{n-1}={2u}_{i,j,k}^n-u_{i,j,k}^{n+1}+\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δx}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i+m,j,k}^n+u_{i-m,j,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δy}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j+m,k}^n+u_{i,j-m,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δz}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j,k+m}^n+u_{i,j,k-m}^n)]\end{array}---\left(12\right)$

本发明的方法对有限差分外推算子进行优化，在不降低计算效率及不改变计算结果的条件下减少存储。其具体实现如下。

一、现有技术的波场外推实现过程

震源点波场外推是根据前一时刻(it-1)的波场和当前时刻(it)的波场，由有限差分外推算子(公式(11))计算出下一时刻(it+1)的波场，以此循环往复(从时刻1到IT)，直到将最后时刻(IT)波场算出(图2)。

it＝1时，dw1+dw2→dw3；

it＝2时，dw2+dw3→dw1；

it＝3时，dw3+dw1→dw2；

it＝4时，dw1+dw2→dw3；

it＝5时，dw2+dw3→dw1；

…………

根据这一实现过程，原方法是开设三个数组空间dw1，dw2，dw3，它们分别存放前一时刻(it-1)的波场数据、当前时刻(it)的波场数据以及下一时刻(it+1)的波场数据。然后，根据当前时间it对3整除取余数来进行判断：

若it％3＝1，即余数为1，则dw1+dw2→dw3；

若it％3＝2，即余数为2，则dw2+dw3→dw1；

若it％3＝0，即余数为0，则dw3+dw1→dw2；

然而，接收点波场逆向外推与震源点波场正向外推正相反，它是根据后一时刻(it+1)与当前时刻(it)的波场数据由有限差分外推算子(公式(12))计算出前一时刻(it-1)的波场数据，直到将最初时刻(it＝0)的波场数据计算出。原方法是开设三个数组空间dr1，dr2，dr3，它们分别存放后一时刻(it+1)的波场数据、当前时刻(it)的波场数据以及前一时刻(it-1)的波场数据。然后，根据当前时间it对3整除取余数来进行判断：

若it％3＝1，则dr1+dr2→dr3；

若it％3＝2，则dr2+dr3→dr1；

若it％3＝0，则dr3+dr1→dr2；

二、根据本发明的实施例的波场外推实现过程

上述实现过程中完全可以通过两个差分函数解决数据互换更替。以it＝1时刻为例，此时需要dw1和dw2的波场数据来求出dw3波场数据，即dw1+dw2→dw3。而此时dw1数组原数据已经没有用了。因此，可以将dw3波场的数据仍然放在dw1数组中。这样，就可以节省1/3的数组空间。因此，在震源点波场外推计算时开设两个数据空间dw1，dw2，分别存放it-1和it时刻的波场数据，而下一时刻的波场数据存放在it-1波场数组中。

因此，现有的波场外推过程可以修改为：

it＝1时，dw1+dw2→dw1；

it＝2时，dw2+dw1→dw2；

it＝3时，dw1+dw2→dw1；

it＝4时，dw2+dw1→dw2；

it＝5时，dw1+dw2→dw1；

……。

在这种情况下，只需两个判断条件即可：

若it％2＝1，则dw1+dw2→dw1；

否则，dw2+dw1→dw2。

同样，在计算接收点波场时开设两个数组空间dr1，dr2，它们分别用于存放it+1时刻和it时刻的波场数据，it-1时刻的波场数据存放在it+1波场数组中，即：

若it％2＝1，则dr1+dr2→dr1；

否则，dr2+dr1→dr2。

综上可知，本发明的总体方法通过计算机来实现时，具有如下步骤，如图1所示。

首先，通过震源点正向外推过程得到各个时刻的震源点正向外推波场数据；然后，通过接收点逆向外推过程得到各个时刻的接收点逆向外推波场数据；最后，将同一时刻下的震源点正向外推波场数据与接收点逆向外推波场数据进行相关求和得到深度域成像值。接下来，根据合适、恰当的成像条件进行叠加成像。

在本发明中，震源点正向外推过程可以这样实现，参见图2。在图2中，过程开始于步骤S101。在步骤S101中，模拟震源点的前一时刻和当前时刻的波场数据。接下来在步骤S102中，分配两个存储空间，其用于分别存储震源点的前一时刻和当前时刻的波场数据。

在步骤S103中，基于两个存储空间中所存储的震源点的前一时刻和当前时刻的波场数据根据全波波动方程计算出下一时刻的波场数据，并将下一时刻的波场数据保存到用于存储前一时刻的波场数据的存储空间中。在本实施例中，该存储空间对应于上面所述的计算机语言中的数组。

在步骤S104中，判断迭代计算结束条件是否满足，如果当前时刻为最终时刻，则输出用于存储前一时刻的波场数据的存储空间中的数据。该数据如前所述即为当前最终时刻的波场数据。如果未到达最终时刻，则继续回到步骤S103进行下一时刻的波场数据的迭代计算。在迭代计算过程中，为了达到一定的精度，采用如上公式(11)中截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维正演模拟高阶差分方程作为外推算子来计算下一时刻的波场数据。

类似地，接收点逆向外推过程包括以下步骤，参见图3。在图3中，过程开始于步骤S201。在步骤S201中，获取接收点的下一时刻和当前时刻的波场数据。接下来，在步骤S202中，与震源点相同，也仅分配两个存储空间，其用于分别存储接收点的下一时刻和当前时刻的波场数据；

在步骤S203中，基于两个存储空间中所存储的接收点的下一时刻和当前时刻的波场数据根据全波波动方程计算出前一时刻的波场数据，并将前一时刻的波场数据保存到用于存储下一时刻的波场数据的存储空间中。同样，在本实施例中，该存储空间对应于上面所述的计算机语言中的数组。

在步骤S204中，判断迭代计算结束条件是否满足，如果当前时刻为起始时刻，则输出用于存储下一时刻的波场数据的存储空间中的数据。该数据为当前起始时刻的波场数据。如果未到起始时刻，则继续回到步骤S203进行前一时刻的波场数据的迭代计算。类似的是，为了达到一定的精度，采用如上公式(12)中截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维逆时深度偏移高阶差分方程作为外推算子来计算前一时刻的波场数据。

本发明所设计的基于GPU的叠前逆时偏移方法，在不改变计算结果及不降低计算效率的前提下，将GPU数据空间需求减少近1/3，提高了GPU空间的利用率，使得该方法更加适应于大规模数据处理。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (10)

1.一种用于基于GPU来获取叠前逆时偏移的方法，其特征在于，所述方法包括：

通过震源点正向外推过程得到各个时刻的震源点正向外推波场数据；

通过接收点逆向外推过程得到各个时刻的接收点逆向外推波场数据；

将同一时刻下的所述震源点正向外推波场数据与所述接收点逆向外推波场数据进行相关求和得到深度域成像值；

其中，所述震源点正向外推过程包括以下步骤：

S101、模拟震源点的前一时刻和当前时刻的波场数据；

S102、分配两个存储空间分别存储所述震源点的前一时刻和当前时刻的波场数据；

S103、基于两个存储空间中所存储的所述震源点的前一时刻和当前时刻的波场数据根据全波波动方程计算出下一时刻的波场数据，并将所述下一时刻的波场数据保存到用于存储所述前一时刻的波场数据的存储空间中；

S104、判断当前时刻是否为最终时刻，如果是，则输出用于存储所述前一时刻的波场数据的存储空间中的数据，如果否，继续回到步骤S103。

2.如权利要求1所述的用于基于GPU来获取叠前逆时偏移的方法，其特征在于，所述接收点逆向外推过程包括以下步骤：

S201、获取接收点的下一时刻和当前时刻的波场数据；

S202、分配两个存储空间分别存储所述接收点的下一时刻和当前时刻的波场数据；

S203、基于两个存储空间中所存储的所述接收点的下一时刻和当前时刻的波场数据根据全波波动方程计算出前一时刻的波场数据，并将所述前一时刻的波场数据保存到用于存储所述下一时刻的波场数据的存储空间中；

S204、判断当前时刻是否为起始时刻，如果是，则输出用于存储所述下一时刻的波场数据的存储空间中的数据，如果否，继续回到步骤S203。

3.如权利要求2所述的用于基于GPU来获取叠前逆时偏移的方法，其特征在于，所述全波波动方程为：

$\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2}+\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{v^2(x,y,z)}\frac{{\∂}^{2}u}{{\∂t}^2}$

其中，u(x,y,z,t)为地表记录的压力波场；v(x,y,z)为纵横向可变的介质速度。

4.如权利要求3所述的用于基于GPU来获取叠前逆时偏移的方法，其特征在于，对压力波场进行变形处理并代入所述全波波动方程中，得到截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维正演模拟高阶差分方程为：

$\begin{array}{c}{u}_{i,j,k}^{n+1}={2u}_{i,j,k}^n-u_{i,j,k}^{n-1}+\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δx}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i+m,j,k}^n+u_{i-m,j,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δy}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j+m,k}^n+u_{i,j-m,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δz}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j,k+m}^n+u_{i,j,k-m}^n)]\end{array}$

在步骤S103中，基于所述截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维正演模拟高阶差分方程来计算出震源点的下一时刻的波场数据，上标n表示当前时刻的数据，n-1表示前一时刻的数据，n+1表示下一时刻的数据。

5.如权利要求3所述的用于基于GPU来获取叠前逆时偏移的方法，其特征在于，对压力波场进行变形处理并代入所述全波波动方程中，得到截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维逆时深度偏移高阶差分程为：

$\begin{array}{c}{u}_{i,j,k}^{n-1}={2u}_{i,j,k}^n-u_{i,j,k}^{n+1}+\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δx}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i+m,j,k}^n+u_{i-m,j,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δy}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j+m,k}^n+u_{i,j-m,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δz}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j,k+m}^n+u_{i,j,k-m}^n)]\end{array}$

在步骤S203中，基于所述截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维逆时偏移高阶差分方程来计算出接收点的前一时刻的波场数据。

6.一种用于基于GPU来获取叠前逆时偏移的装置，其特征在于，所述装置包括：

震源点正向外推模块，其通过震源点正向外推过程得到各个时刻的震源点正向外推波场数据；

接收点逆向外推模块，其通过接收点逆向外推过程得到各个时刻的接收点逆向外推波场数据；

深度域成像值计算模块，其用于将同一时刻下的所述震源点正向外推波场数据与所述接收点逆向外推波场数据进行相关求和得到深度域成像值；

其中，所述震源点正向外推模块包括以下单元：

模拟单元，其用于模拟震源点的前一时刻和当前时刻的波场数据；

第一存储空间分配单元，其用于分配两个存储空间以分别存储所述震源点的前一时刻和当前时刻的波场数据；

第一波场数据计算单元，其用于基于两个存储空间中所存储的所述震源点的前一时刻和当前时刻的波场数据根据全波波动方程计算出下一时刻的波场数据，并将所述下一时刻的波场数据保存到用于存储所述前一时刻的波场数据的存储空间中；

第一判断单元，其用于判断当前时刻是否为最终时刻，如果是，则输出用于存储所述前一时刻的波场数据的存储空间中的数据，如果否，继续到所述第一波场数据计算单元中计算波场数据。

7.如权利要求6所述的用于基于GPU来获取叠前逆时偏移的装置，其特征在于，所述接收点逆向外推模块包括以下单元：

获取单元，其用于获取接收点的下一时刻和当前时刻的波场数据；

第二存储空间分配单元，其用于分配两个存储空间分别存储所述接收点的下一时刻和当前时刻的波场数据；

第二波场数据计算单元，其用于基于两个存储空间中所存储的所述接收点的下一时刻和当前时刻的波场数据根据全波波动方程计算出前一时刻的波场数据，并将所述前一时刻的波场数据保存到用于存储所述下一时刻的波场数据的存储空间中；

第二判断单元，其用于判断当前时刻是否为起始时刻，如果是，则输出用于存储所述下一时刻的波场数据的存储空间中的数据，如果否，继续到所述第二波场数据计算单元中计算波场数据。

8.如权利要求7所述的用于基于GPU来获取叠前逆时偏移的装置，其特征在于，所述全波波动方程为：

$\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2}+\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{v^2(x,y,z)}\frac{{\∂}^{2}u}{{\∂t}^2}$

其中，u(x,y,z,t)为地表记录的压力波场；v(x,y,z)为纵横向可变的介质速度。

9.如权利要求8所述的用于基于GPU来获取叠前逆时偏移的装置，其特征在于，对压力波场进行变形处理并代入所述全波波动方程中，得到截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维正演模拟高阶差分方程为：

$\begin{array}{c}{u}_{i,j,k}^{n+1}={2u}_{i,j,k}^n-u_{i,j,k}^{n-1}+\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δx}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i+m,j,k}^n+u_{i-m,j,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δy}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j+m,k}^n+u_{i,j-m,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δz}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j,k+m}^n+u_{i,j,k-m}^n)]\end{array}$

在第一波场数据计算单元中，基于所述截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维正演模拟高阶差分方程来计算出震源点的下一时刻的波场数据。

10.如权利要求8所述的用于基于GPU来获取叠前逆时偏移的装置，其特征在于，对压力波场进行变形处理并代入所述全波波动方程中，得到截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维逆时深度偏移高阶差分程为：

$\begin{array}{c}{u}_{i,j,k}^{n-1}={2u}_{i,j,k}^n-u_{i,j,k}^{n+1}+\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δx}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i+m,j,k}^n+u_{i-m,j,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δy}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j+m,k}^n+u_{i,j-m,k}^n)]\\ +\frac{1}{2}{\left(\frac{\mathrm{v\Δt}}{\mathrm{\Δz}}\right)}^2[{\mathrm{\ω}}_0{u}_{i,j,k}^n+\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{\mathrm{\ω}}_m(u_{i,j,k+m}^n+u_{i,j,k-m}^n)]\end{array}$

在第二波场数据计算单元中，基于所述截断误差为O(Δx^M,Δy^M,Δz^M,Δt²)的统一的三维逆时偏移高阶差分方程来计算出接收点的前一时刻的波场数据。