CN102279415B

CN102279415B - 基于图形处理器计算Fourier积分法单程波深度偏移的方法

Info

Publication number: CN102279415B
Application number: CN 201110190524
Authority: CN
Inventors: 刘红伟; 刘洪�; 李博; 佟小龙; 刘钦
Original assignee: Institute of Geology and Geophysics of CAS
Current assignee: Institute of Geology and Geophysics of CAS
Priority date: 2011-07-08
Filing date: 2011-07-08
Publication date: 2013-07-24
Anticipated expiration: 2031-07-08
Also published as: CN102279415A

Abstract

本发明公开了基于图形处理器计算Fourier积分法单程波深度偏移的方法，包括下述步骤：（1）设计震源子波；（2）读入单炮数据，依据观测系统确定偏移孔径；（3）将偏移孔径内的速度模型、单炮数据、震源子波传递给图形处理器；（4）将数据变换到频率波数域；（5）进行相移和Fourier积分运算；（6）进行成像运算；（7）重复步骤（4）-（6），至延拓到最深层，保存成像结果，并进行下一炮偏移；（8）重复步骤（2）-（7），至最后一炮，输出成像结果；上述步骤（4）-(7)在图形处理器中进行。本发明的方法充分利用了图形处理器在矩阵运算中的极大优势，改善了算法的处理效率。

Description

基于图形处理器计算Fourier积分法单程波深度偏移的方法

技术领域

本发明涉及计算Fourier积分法单程波深度偏移，尤其涉及基于图形处理器计算Fourier积分法单程波深度偏移的方法。

背景技术

叠前单程波波动方程偏移技术是解决地震成像问题的有力工具，自其诞生之日起，究竟选择何种近似方式以提高波场延拓算子的精度，一直是学术界研究的热点问题。单程波方法是Claerbout教授上世纪70年代引入地震偏移领域的，因其能够处理多波至问题以及内存要求小等优点在偏移领域倍受青睐。对于单程波偏移方法而言，一个关键的问题是寻找一种单程波延拓算子，使其在速度横向变化剧烈的介质中，能够精确地模拟波场的大角度传播。过去30多年间涌现出很多种单程波延拓算子，其中包括频率空间域算子，比如单程波有限差分法(Finite Difference, Claerbout, 1985)；频率波数域算子，比如相移法 (Phase-Shift, Gazdag, 1978)和傅里叶变换法 (Fourier Transform, Stolt, 1978)；频率波数空间混合域算子；比如相移加插值法(Phase Shift Plus Interpolation, Gazdag and Sguazzero, 1984), 分裂步相移法 (Split-Step Fourier, Stoffa etc., 1990), 傅里叶有限差分 (Fourier Finite Difference, Ristow and Rühl, 1994), 广角相位屏法(Wide-angle Screen, Xie and Wu, 1998)，广义屏方法(Generalized Screen, Le Rousseau and de Hoop, 2001)，最优分裂步方法(Optimum Split-step, Liu, L. and J. Zhang, 2006)，最优可分近似方法 (Optimization Approximation, Chen, J.-B., and H. Liu, 2004, 2006)。

与有限差分相比，基于傅里叶变换的方法在应用到三维情形时不存在双向分裂误差(Brown, 1983)，并且对于粗网格速度模型和高频率地震数据几乎不存在数值频散。要保证傅里叶变换法的高效性，波场在空间域和波数域变换时必须借助快速傅里叶变换(FFT)来实现，这就要求在这些方法应用时必须对速度模型或者单程波延拓算子做不同程度的近似或假设。Stolt的方法假设速度模型为常速，是最快速的单程波方法；Gazdag的相移法允许速度沿深度方向变化，但不允许横向变化；为了处理速度的横向变化，Gazdag提出相移加插值的方法采用不同的参考速度进行相移运算，在不同相移结果之间插值，可以处理速度的横向变化，但计算成本太高；Staffa的分裂步相移法在相移法的基础上引入了透镜相作为零阶修正项，随后，Le Rousseau的广义屏方法在分裂步相移法的基础上引入高阶修正项，然而，这些修正项的引入只对速度横向变化较弱的模型有效，对于速度横向变化剧烈的模型，偏离垂直方向大角度传播的波场误差仍然很大，原因在于广义屏算子对单程波延拓算子做了双重近似：垂直波数的高阶泰勒级数截断以及高阶修正项指数函数的一阶泰勒级数截断。为了提高近似的精度，Liu L.提出最优分裂步方法，在保持广义屏计算结构的基础上，直接对单程波延拓算子做近似，运用最优近似理论求取高阶修正项的系数，在几乎相同的计算成本的情况下，有效改进了波场延拓的精度；Chen Jingbo的最优可分近似方法采用函数最优可分近似理论将单程波延拓算子近似为空间变量函数与波数变量函数乘积的形式，这种方法收敛速度快并且不受分支点问题的影响，是一种快速有效且高精度的单程波方法。

上面提到的傅里叶方法有一个共同的特点，即为了应用快速傅里叶变换，对单程波延拓算子做了不同程度的近似，使空间域变量和波数域变量分离。Margrave在1999年提出一种非稳态相移法(Non-Stationary Phase Shift)，这种方法不对单程波算子做近似，但是Robert J. Ferguson在2002年指出这种方法仅对速度模型在横向分段为常数时才有实际应用价值。

发明内容

针对上述传统单程波方法的缺点，本发明提供了基于图形处理器计算傅里叶积分法(Fourier Integral)单程波深度偏移。

在本发明中包括下述两个技术概念，申请人在此作出解释：

傅立叶积分法计算单程波深度偏移，这种方法不对单程波算子和速度模型做近似，适用于速度横向任意变化的速度模型，能精确模拟偏离垂直方向接近90度传播的波场。

其计算方法描述如下：以二维情形为例，频率空间域的单程波方程如下：

Figure 2011101905248100002DEST_PATH_IMAGE001

(1)

运用傅里叶变换，可以得到单程波方程频率—空间—波数域(混合域)的积分解形式如下：

(2)

对于积分解公式(2)，由于波场延拓算子

Figure 2011101905248100002DEST_PATH_IMAGE003

中，波数变量

与空间变量是耦合在一起的，不能借助快速傅里叶变换(FFT)求解，每个空间位置x处的波场都需要求一次傅里叶积分，计算量巨大。

由积分解公式(2)可知，傅里叶积分变换实际上是一个矩阵乘法运算，频率空间域的波场

和频率波数域的波场分别是一个二维矩阵，而波场延拓算子是一个三维矩阵，如图1所示。求解

中的第i行第j列元素需要矩阵

中的第i行元素与波场延拓算子

中的第i页第j列元素做乘加运算得到。

和二维问题类似，三维情形下，频率空间域的波场和频率波数域的波场分别是两个三维矩阵；而波场延拓算子是一个五维矩阵，我们将三维矩阵按二维矩阵处理，将五维矩阵按三维矩阵处理，则在形式上，三维问题退化成二维问题，如图2所示。

然而如果仅仅采用传统的傅立叶积分法用单程波深度偏移，即便考虑到高速的中央处理器CPU，其处理方法也太慢，计算成本太高。

在上述的基础上，本发明在傅立叶积分法单程波深度偏移的基础上，通过图形处理器(Graphic Processing Unit，GPU)在矩阵乘法方面的巨大优势，大幅度地提高算法的计算效率，使得傅里叶积分法在单程波深度偏移的实际应用成为一种快速、成本低、可应用的技术方案。

在本发明中，GPU的全称是图形处理器(Graphic Processing Unit)，俗称“显卡”，在市场上目前有多家公司有研发生产此类硬件，包括NVIDIA、AMD（原ATI）、VIA等x86体系显卡生产商、POWERVR、ARM等基于RISC体系的显卡生产商。

本发明公开了一种基于图形处理器计算Fourier积分法单程波深度偏移的方法，包括下述步骤：

（1）设计震源子波；

（2）读入单炮数据，依据观测系统确定偏移孔径；

（3）将偏移孔径内的速度模型、单炮数据、震源子波传递给图形处理器；

（4）将数据变换到频率波数域；

（5）进行相移和Fourier积分运算；

（6）进行成像运算；

（7）重复步骤（4）-（6），至延拓到最深层，保存成像结果，并进行下一炮偏移；

（8）重复步骤（2）-（7），至最后一炮，输出成像结果；

上述步骤（4）-(7)在图形处理器中进行。

在本发明中，步骤（3）所述的偏移孔径内的速度模型、单炮数据、震源子波传递给图形处理器，是指将偏移孔径内的速度模型、单炮数据、震源子波的数据传递到GPU的内存中（通常称之为显存），所述的GPU的内存是直接位于GPU芯片上的内存。

在本发明中，所述图形处理器尤其可以是近些年广泛兴起的支持通用计算功能的图形处理器，此类处理器通常包括基于统一计算设备架构平台CUDA的图形处理器，此图形处理器由显卡厂商NVIDIA推出；基于开放通用计算模型OpenCL的图形处理器，此类图形处理器包括NVIDIA、AMD、Intel、VIA等均有推出；基于微软DirectCompute标准的图形处理器，此类图形处理器包括NVIDIA、AMD、Intel、VIA等均有推出；优选的，基于成熟度和开发快捷的需要，本发明中所用的图形处理器优选为基于统一计算设备架构平台CUDA的图形处理器。本领域技术人员根据本发明的实质精神采用其他类型的支持通用计算功能的图形处理器依旧属于本发明保护范围。

基于上述基础，在本发明中，在步骤（3）中，偏移孔径内的速度模型、单炮数据、震源子波数据由中央处理器通过图形处理器的CUDA编程接口，传送到图形处理器中的渲染管道，到达可编程片段处理器，以多个线程的方式并行计算，通过图形处理器的共享存储器内存为线程提供数据，从而克服了图形处理器与中央处理器内存通讯导致的时间浪费；所述步骤（4）-(7)在图形处理器的可编程片段处理器中进行，其具体的计算过程如下：

步骤（4）调用统一计算设备架构平台CUDA的函数库CUFFT将数据变换到频率波数域；

步骤（5）调用统一计算设备架构平台CUDA的函数kintegralmig进行相移和Fourier积分运算；

步骤（6）调用统一计算设备架构平台CUDA的函数imag_mul进行成像运算。

在完成上述运算后，步骤（7）用于判断是否延拓到最深层，如否，重复4-6步，若已经到最深层，将成像结果保存，并进行下一炮偏移；

步骤（8）用于判断是否是最后一炮，若否，重复2-7步，若是，输出成像结果，偏移完成。

尽管上述给出的是采用基于统一计算设备架构平台CUDA的图形处理器实现本发明目的的实现方法，然而本领域技术人员根据上述实现方法的实质精神采用其它支持通用计算功能的图形处理器并调用其中对应的函数库或相应函数实现相应的正变换，相移加反变换以及成像等功能也是可以为本领域技术人员所能知悉的，依旧属于本发明的保护范围。

在上述所述中，所用的词汇“渲染管道、可编程片段处理器”对于本领域技术人员是广为理解的，其他类似的称呼如管线、ROP、处理单元、流处理器等亦可为本领域技术人员所理解为与本发明所用词汇含义一致，改用不同的名称，只要其针对的是图新处理器上相同的硬件单元结构，依旧属于本发明的保护范围。

在上述基础上，本发明的基于图形处理器计算Fourier积分法单程波深度偏移的方法进一步的还包括将成像结果从图形处理器传输回中央处理器和内存的步骤。

与传统的单程波偏移方法相比，本发明公开的基于图形处理器计算Fourier积分法单程波深度偏移的方法具有如下优势：

1、与国内外的现有单程波偏移方法相比，本方法对波场延拓算子不做任何近似，直接延用其原始形式，因此，在存有横向变速的介质中，可精确模拟偏离垂直方向接近90度传播的波场。数值试验表明，傅里叶积分法单程波偏移不存在数值频散，算法稳定，且将之扩展到三维问题时不存在方向分裂误差。

2、与传统的完全基于CPU的算法相比，本发明的方法整个偏移过程包括正变换，相移加反变换以及成像条件都是在图形处理器上实现的，最大限度的减少中央处理器内存和图形处理器内存之间的通讯，大大提升程序算法运行的效率。数据传输只发生在将速度文件以及地震数据传输到图形处理器内存和将成像结果传输回中央处理器内存过程中，借助图形处理器在矩阵乘法方面的巨大优势，大幅度地提高了算法的计算效率。并且本发明的方法利用图形处理器的共享存储器为计算线程提供数据，降低了中央处理器和图形处理器的通信次数和通信代价，省掉了图形处理器与中央处理器及其内存频繁通讯导致的速度延迟，显著提高了方法的执行效率和速度。

附图说明

附图1为波场延拓傅里叶积分与矩阵乘积运算示意图，其中频率空间域的波场和频率波数域的波场分别是一个二维矩阵，而波场延拓算子是一个三维矩阵；

附图2为三维波场延拓傅里叶积分与矩阵乘积运算示意图，其中频率空间域的波场和频率波数域的波场分别是两个三维矩阵，按二维矩阵处理；而波场延拓算子是一个五维矩阵，按三维矩阵处理；

附图3为本发明基于图形处理器计算Fourier积分法单程波深度偏移的方法流程图；

附图4为Marmousi模型；

附图5为采用本发明基于图形处理器计算Fourier积分法单程波深度偏移的方法进行Marmousi模型单炮偏移测试偏移结果。

具体实施方式

下面结合附图结合实施例对发明进行进一步阐述，本领域技术人员应知，本发明也可以有多种不同形式或使用不同的图形处理器或者不同的通用计算函数实施，因此不应认为它局限于说明书列出的实施例。本发明的实际应用效果与具体的图形处理器硬件相关，在一定范围内，图形处理器所拥有的流处理器越多、速度越快，速度提升的效果将越好。

参考附图1、2，公开了本发明的基于图形处理器计算Fourier积分法单程波深度偏移的方法对于二维模型、三维模型的傅里叶积分与矩阵乘积运算的数学方法，显示了本发明方法的数学原理和计算方式。

参考附图3，公开了本发明的基于图形处理器计算Fourier积分法单程波深度偏移的方法流程图，图3中用设备一词代指图形处理器。包括下述步骤：（1）设计震源子波；（2）读入单炮数据，依据观测系统确定偏移孔径；（3）将偏移孔径内的速度模型、单炮数据、震源子波传递给图形处理器；（4）将数据变换到频率波数域；（5）进行相移和Fourier积分运算；（6）进行成像运算；（7）重复步骤（4）-（6），至延拓到最深层，保存成像结果，并进行下一炮偏移；（8）重复步骤（2）-（7），至最后一炮，输出成像结果；上述步骤（4）-(7)在图形处理器中进行。

参考图4所示，Marmousi模型(Bourgeois, 1991)是法国石油研究所(Institut Français du Pétrole，简称IFP)1988年制造的，地质原型是穿过Cuanza盆地Quenguela地区北部的一条地质剖面。Marmousi模型数据在1990年第52届EAGE年会地震数据反演讨论会上被使用，Marmousi模型已经成为工业界标准的经典模型。

附图5是利用本发明的基于图形处理器计算Fourier积分法单程波深度偏移的方法进行Marmousi模型单炮偏移测试偏移结果。申请人分别使用了目前主流的CPU和GPU做单炮偏移测试（CPU为Intel (R) E5200 2.5GHz，GPU为NVIDIA Tesla C1060），CPU版本程序的运行时间是7412s(2h3m32s)，而GPU版本的运行时间是3s，加速比是2470，证实本发明的发明显著改善了Fourier积分法单程波深度偏移的计算效率。

Claims

1.基于图形处理器计算Fourier积分法单程波深度偏移的方法，包括下述步骤：

(1)设计震源子波；

(2)读入单炮数据，依据观测系统确定偏移孔径；

(3)将偏移孔径内的速度模型、单炮数据、震源子波传递给图形处理器；

(4)将数据变换到频率波数域；

(5)进行相移和Fourier积分运算；

(6)进行成像运算；

(7)重复步骤(4)-(6)，至延拓到最深层，保存成像结果，并进行下一炮偏移；

(8)重复步骤(2)-(7)，至最后一炮，输出成像结果；

上述步骤(4)-(7)在图形处理器中进行，所述图形处理器是基于统一计算设备架构平台CUDA的图形处理器。

2.根据权利要求1所述的基于图形处理器计算Fourier积分法单程波深度偏移的方法，其特征在于所述步骤(3)偏移孔径内的速度模型、单炮数据、震源子波上传于图形处理器的内存中。

3.根据权利要求1所述的基于图形处理器计算Fourier积分法单程波深度偏移的方法，其特征在于步骤(3)的偏移孔径内的速度模型、单炮数据、震源子波数据由中央处理器通过图形处理器的CUDA编程接口，传送到图形处理器中的渲染管道，到达可编程片段处理器；所述步骤(4)-(7)在图形处理器的可编程片段处理器中进行，以多个线程的方式并行计算，通过图形处理器的共享存储器为线程提供数据。

4.根据权利要求1所述的基于图形处理器计算Fourier积分法单程波深度偏移的方法，其特征在于所述步骤(4)调用统一计算设备架构平台CUDA的函数库CUFFT将数据变换到频率波数域。

5.根据权利要求1所述的基于图形处理器计算Fourier积分法单程波深度偏移的方法，其特征在于所述步骤(5)调用统一计算设备架构平台CUDA的函数kintegralmig进行相移和Fourier积分运算。

6.根据权利要求1所述的基于图形处理器计算Fourier积分法单程波深度偏移的方法，其特征在于所述步骤(6)调用统一计算设备架构平台CUDA的函数imag_mul进行成像运算。

7.根据权利要求1所述的基于图形处理器计算Fourier积分法单程波深度偏移的方法，进一步的包括将成像结果从图形处理器传输回中央处理器和内存的步骤。

8.根据权利要求1所述的基于图形处理器计算Fourier积分法单程波深度偏移的方法，针对二维模型，频率空间域的单程波方程为：

运用傅里叶积分变换，所得单程波方程频率-空间-波数域(混合域)的积分解形式为：

其中，波场延拓算子exp(ik_zΔz)中，波数变量k_x与空间变量x是耦合在一起的，频率空间域的波场D(x，z+Δz，ω)和频率波数域的波场d(k_x，z，ω)分别是一个二维矩阵，波场延拓算子exp(ik_zΔz)是一个三维矩阵，求解D(x，z+Δz，ω)中的第i行第j列元素通过矩阵d(k_x，z，ω)中的第i行元素与波场延拓算子exp(ik_zΔz)中的第i页第j列元素做乘加运算得到。

9.根据权利要求8所述的基于图形处理器计算Fourier积分法单程波深度偏移的方法，针对三维模型，频率空间域的波场和频率波数域的波场分别是两个三维矩阵；而波场延拓算子是一个五维矩阵，将三维矩阵按二维矩阵处理，将五维矩阵按三维矩阵处理。