CN103675895A - 一种利用gpu提高波场延拓计算效率的方法 - Google Patents

一种利用gpu提高波场延拓计算效率的方法 Download PDF

Info

Publication number
CN103675895A
CN103675895A CN201210315284.4A CN201210315284A CN103675895A CN 103675895 A CN103675895 A CN 103675895A CN 201210315284 A CN201210315284 A CN 201210315284A CN 103675895 A CN103675895 A CN 103675895A
Authority
CN
China
Prior art keywords
wave field
field
gpu
carry out
fourier transform
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.)
Pending
Application number
CN201210315284.4A
Other languages
English (en)
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.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
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 China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201210315284.4A priority Critical patent/CN103675895A/zh
Publication of CN103675895A publication Critical patent/CN103675895A/zh
Pending legal-status Critical Current

Links

Images

Landscapes

  • Image Analysis (AREA)

Abstract

本发明提供了一种利用GPU提高波场延拓计算效率的方法，属于地震资料处理领域。所述方法首先扩大震源波场和记录波场的数据网格，然后利用GPU进行多线程并行计算，并应用CUDA提供的二维傅里叶变换函数实现波场延拓。本发明的方法简单实用且效率高，利用本发明得到的成像效果与CPU单程波动方程叠前深度偏移得到的成像效果一致，但是计算效率提高了30倍以上。

Description

一种利用GPU提高波场延拓计算效率的方法
技术领域
本发明属于地震资料处理领域，具体涉及一种利用GPU提高波场延拓计算效率的方法。
背景技术
波动方程叠前深度偏移最核心的工作是地震波场的深度外推。深度外推过程必须依赖于波场外推算子，而波场外推算子实质上是求解单向波动方程。单炮道集波动方程叠前深度偏移成像要分别对炮点波场和检波点波场进行下行波和上行波外推，利用激发时间成像条件提取每一外推层的成像值。
波动方程叠前深度偏移成像计算是以单炮为一个计算单元，单炮的波场延拓，每外推一步，需要两次二维傅里叶正变换，一次相移，两次二维傅里叶反变换，一次时移，其中二维傅里叶变换计算量最大，因此需要提高其计算效率。
发明内容
本发明的目的在于解决上述现有技术中存在的难题，提供一种利用GPU提高波场延拓计算效率的方法，根据GPU(Graphic Processing Unit图形处理器)的计算特点，应用CUDA(Compute Unified Device Architecture，统一计算设备架构)提供的傅里叶变换函数，有效的提高单程波动方程叠前深度偏移成像的计算效率，缩短单程波动方程叠前深度偏移成像的处理周期，快速地进行波场延拓计算，提高波场延拓的计算效率，提高炮域叠前深度偏移的效率。
本发明是通过以下技术方案实现的：
一种利用GPU提高波场延拓计算效率的方法，所述方法首先扩大震源波场和记录波场的数据网格，然后利用GPU进行多线程并行计算，并应用CUDA提供的二维傅里叶变换函数实现波场延拓。
所述方法包括以下步骤：
(1)扩大震源波场和记录波场的数据网格得到处理后的震源波场和记录波场；
(2)应用裂步傅里叶法(SSF)对所述处理后的震源波场和记录波场进行波场延拓；
(3)设延拓深度iz＝1；
(4)利用CUDA提供的二维傅里叶变换函数对震源波场进行二维正傅里叶变换；
(5)利用CUDA提供的二维傅里叶变换函数对记录波场进行二维正傅里叶变换；
(6)参考速度场的相移计算：利用步骤用(4)和步骤(5)的结果进行相移计算，参考速度场是平均速度，即输入速度场的平均值(请参考Stoffa P L.Split-step Fourier migration.Geophysics，1990，55(2)；410～421)；
(7)利用CUDA提供的二维傅里叶变换函数对震源波场进行二维反傅里叶变换；
(8)利用CUDA提供的二维傅里叶变换函数对记录波场进行二维反傅里叶变换；
(9)扰动速度场的波场计算与相关成像：利用步骤(7)和步骤(8)的结果进行相移计算，扰动速度场是输入速度场与参考速度场的差值(请参考StoffaP L.Split-step Fourier migration.Geophysics，1990，55(2)；410～421)；
(10)判断iz+1是否大于NZ，所述NZ＝延拓深度/延拓步长，如果是，则转入步骤(11)，如果否，则返回步骤(4)；
(11)结束。
所述步骤(1)是这样实现的：对震源波场和记录波场的长度进行镶边处理，即当震源波场和记录波场的长度小于1000时，将其长度镶边处理为128的倍数，对于不足的数据处填充零；当震源波场和记录波场的长度大于1000而小于2048时，将其长度镶边处理成256的倍数，对于不足的数据处填充零。
所述步骤(2)至(11)都是在GPU上进行的。
与现有技术相比，本发明的有益效果是：本发明的方法简单实用且效率高，利用本发明得到的成像效果与CPU单程波动方程叠前深度偏移得到的成像效果一致，但是计算效率提高了30倍以上。
附图说明
图1是本发明方法的波场延拓的示意图。
图2是本发明方法的GPU运算的流程图。
图3-1是本发明方法实施例中利用CPU计算得到的盐丘模型偏移剖面。
图3-2是本发明方法实施例中利用GPU计算得到的盐丘模型偏移剖面。
图4-1是本发明方法实施例中利用CPU计算得到的负向构造偏移剖面。
图4-2是本发明方法实施例中利用GPU计算得到的负向构造偏移剖面。
图5是利用CPU和GPU进行波场延拓的计算效率对比图。
图6是本发明方法的波场延拓的程序框图。
具体实施方式
下面结合附图对本发明作进一步详细描述：
如图6所示(图6中的NW表示频率片)，本发明利用GPU进行多线程并行计算，应用CUDA(CUDA是一种将GPU作为数据并行计算设备的软硬件体系，采用类C语言进行软件开发。)提供的二维(2D)傅里叶变换函数，结合波场延拓的计算特点(即二维傅立叶计算密集度大)，采用扩大波场数据网格和多信号同时变换的策略，从而达到提高计算效率的目的，该方法简单实用，效率高。
具体来说，在GPU计算中，内核函数是以线程网格(Grid)的形式组织的，每个线程网格由若干个线程块(Block)组成，而每个线程块又由多个线程(thread)组成。在每次调用GPU计算时，都要指定线程的网格结构<<<dimGrid，dimBlock>>>。在相移和时移的计算中，根据炮点波场和检波点波场数据体分配线程结构：每个线程计算三维体中的一个点，所有的点可以在由一个线程网格计算完成。因为提高了算法的并行度，所以计算效率大幅度提升。同时，将扰动速度场(时移)计算和相关成像计算进行合并计算，可以减少数据一次读写和一次线程调用，进一步提高了计算效率。
本发明对炮域叠前深度偏移算法中广泛使用的二维傅里叶变换CUFFT(由CUDA数学函数库提供)进行了深入的测试分析，发现CUFFT对信号长度比较敏感，当信号的长度为奇数时，运算速度较慢；当信号的长度为偶数时，速度较快；当信号长度是长度为2的幂次时，速度最快。
CUFFT的计算速度非常快，计算时间与信号长度基本上呈正比关系。从这个结果出发，本发明对CUFFT的信号长度做镶边处理，以提高CUFFT的计算效率。但是，如果将信号长度都镶边成2n，那么在提高CUFFT计算效率的同时，也会增加数据的存储量，因此要找一个折中点，也就是在稍微增加数据量的基础上大幅提高CUFFT的计算效率；当信号长度小于1000时，信号长度为128的倍数的有较优的计算效率，当信号长度在2048以内，信号长度为256的倍数的有较优的计算效率。因此本发明提出扩大波场数据网格以提高CUFFT的计算效率。
在CUFFT中，一次计算NW个2DFFT要比NW次计算2DFFT快一倍，因此可以将原来的波场延拓的NW次计算震源波场2D正FFT，记录波场2D正FFT，相移计算，震源波场2D反FFT，记录波场2D反FFT，扰动速度场(时移)计算，相关成像七步计算，更改为一次震源波场2D正FFT(NW)，记录波场2D正FFT(NW)，相移计算，震源波场2D反FFT(NW)，记录波场2D反FFT(NW)，扰动速度场(时移)计算+相关成像六步计算，从而提高波场延拓的计算速度。
用CUDA提供的cufftPlanMany(CUDA傅立叶变换批量处理)进行多信号同时变换(是在(4)(5)(7)(8)这四个步骤中实现的)，具体如下：
cufftPlanMany(&plan，2，dims，NULL，1，0，NULL，1，0，CUFFT C2C，nw)；
cufftExecC2C(plan，d_signal_s，d_signal_s，CUFFT_FORWARD)；
cufftExecC2C(plan，d_signal_r，d_signal_r，CUFFT_FORWARD)；
cufftExecC2C(plan，d_signal_s，d_signal_s，CUFFT_INVERSE)；
cufftExecC2C(plan，d_signal_r，d_signal_r，CUFFT_INVERSE)；所述
图1给出的是波场延拓示意图，单炮道集波动方程叠前深度偏移成像要分别对炮点波场和检波点波场进行下行波和上行波外推，利用激发时间成像条件提取每一外推层的成像值。
图2给出的是GPU运算最终流程，具体如下：在每个波场延拓的过程中，依次进行以下计算：fft正变换，相移计算，fft反变换，时移+相关成像。
图3-1是本发明方法实施例中利用CPU计算得到的盐丘模型偏移剖面，图3-2是本发明方法实施例中利用GPU计算得到的盐丘模型偏移剖面，对比图3-1和图3-2可以看出利用CPU和GPU得到的成像结果是一致的。图4-1是本发明方法实施例中利用CPU计算得到的负向构造偏移剖面，图4-2是本发明方法实施例中利用GPU计算得到的负向构造偏移剖面，对比图4-1和图4-2可以看出利用CPU和GPU得到的成像结果是一致的。
图5表示的是图2的每一步分别在CPU和GPU上的计算时间对比，从图5中可以看出，利用GPU后在计算效率提高了30倍以上。
上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (4)

1.一种利用GPU提高波场延拓计算效率的方法，其特征在于：所述方法首先扩大震源波场和记录波场的数据网格，然后利用GPU进行多线程并行计算，并应用CUDA提供的二维傅里叶变换函数实现波场延拓。
2.根据权利要求1所述的利用GPU提高波场延拓计算效率的方法，其特征在于：所述方法包括以下步骤：
(1)扩大震源波场和记录波场的数据网格得到处理后的震源波场和记录波场；
(2)应用裂步傅里叶法对所述处理后的震源波场和记录波场进行波场延拓；
(3)设延拓深度iz＝1；
(4)利用CUDA提供的二维傅里叶变换函数对震源波场进行二维正傅里叶变换；
(5)利用CUDA提供的二维傅里叶变换函数对记录波场进行二维正傅里叶变换；
(6)参考速度场的相移计算：利用步骤用(4)和步骤(5)的结果进行相移计算，参考速度场是平均速度，即输入速度场的平均值；
(7)利用CUDA提供的二维傅里叶变换函数对震源波场进行二维反傅里叶变换；
(8)利用CUDA提供的二维傅里叶变换函数对记录波场进行二维反傅里叶变换；
(9)扰动速度场的波场计算与相关成像：利用步骤(7)和步骤(8)的结果进行相移计算，扰动速度场是输入速度场与参考速度场的差值；
(10)判断iz+1是否大于NZ，所述NZ＝延拓深度/延拓步长，如果是，则转入步骤(11)，如果否，则返回步骤(4)；
(11)结束。
3.根据权利要求2所述的利用GPU提高波场延拓计算效率的方法，其特征在于：所述步骤(1)是这样实现的：对震源波场和记录波场的长度进行镶边处理，即当震源波场和记录波场的长度小于1000时，将其长度镶边处理为128的倍数，对于不足的数据处填充零；当震源波场和记录波场的长度大于1000而小于2048时，将其长度镶边处理成256的倍数，对于不足的数据处填充零。
4.根据权利要求2所述的利用GPU提高波场延拓计算效率的方法，其特征在于：所述步骤(2)至(11)都是在GPU上进行的。
CN201210315284.4A 2012-08-30 2012-08-30 一种利用gpu提高波场延拓计算效率的方法 Pending CN103675895A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210315284.4A CN103675895A (zh) 2012-08-30 2012-08-30 一种利用gpu提高波场延拓计算效率的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210315284.4A CN103675895A (zh) 2012-08-30 2012-08-30 一种利用gpu提高波场延拓计算效率的方法

Publications (1)

Publication Number Publication Date
CN103675895A true CN103675895A (zh) 2014-03-26

Family

ID=50314014

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210315284.4A Pending CN103675895A (zh) 2012-08-30 2012-08-30 一种利用gpu提高波场延拓计算效率的方法

Country Status (1)

Country Link
CN (1) CN103675895A (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106257308A (zh) * 2016-08-22 2016-12-28 中国石油天然气股份有限公司 地震数据处理方法及系统
CN106324667A (zh) * 2015-07-08 2017-01-11 中国石油化工股份有限公司 基于gpu的三维地震波场模拟计算方法及系统
CN107728199A (zh) * 2017-09-22 2018-02-23 中国地质大学（北京） 基于多gpu并行的多分量各向异性叠前时间偏移加速方法
CN107783184A (zh) * 2016-08-31 2018-03-09 中国科学院地质与地球物理研究所 一种基于多流优化的gpu逆时偏移方法及系统
CN108845355A (zh) * 2018-09-26 2018-11-20 中国矿业大学(北京) 地震偏移成像方法及装置
CN110895350A (zh) * 2018-09-13 2020-03-20 中国石油化工股份有限公司 基于gpu的双程波傅里叶有限差分波场传播方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6021094A (en) * 1998-12-03 2000-02-01 Sandia Corporation Method of migrating seismic records
US20050090989A1 (en) * 2003-10-23 2005-04-28 Kelly Steve M. Method for downward extrapolation of pre-stack seismic data

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6021094A (en) * 1998-12-03 2000-02-01 Sandia Corporation Method of migrating seismic records
US20050090989A1 (en) * 2003-10-23 2005-04-28 Kelly Steve M. Method for downward extrapolation of pre-stack seismic data

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
STOFFA L.P. 等: "Split-step Fourier Migration", 《GEOPHYSICS》 *
张兵 等: "地震叠前深度偏移在CUDA平台上的实现", 《勘探地球物理进展》 *
李振春 等: "《地震叠前成像理论与方法》", 30 September 2011, 中国石油大学出版社 *
李斐 等: "裂步傅里叶真振幅叠前深度偏移", 《石油地球物理勘探》 *
陈辉 等: "基于MPI+OpenMP的多层次并行偏移算法研究", 《成都理工大学学报（自然科学版）》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106324667A (zh) * 2015-07-08 2017-01-11 中国石油化工股份有限公司 基于gpu的三维地震波场模拟计算方法及系统
CN106257308A (zh) * 2016-08-22 2016-12-28 中国石油天然气股份有限公司 地震数据处理方法及系统
CN106257308B (zh) * 2016-08-22 2018-09-04 中国石油天然气股份有限公司 地震数据处理方法及系统
CN107783184A (zh) * 2016-08-31 2018-03-09 中国科学院地质与地球物理研究所 一种基于多流优化的gpu逆时偏移方法及系统
CN107783184B (zh) * 2016-08-31 2020-01-21 中国科学院地质与地球物理研究所 一种基于多流优化的gpu逆时偏移方法及系统
CN107728199A (zh) * 2017-09-22 2018-02-23 中国地质大学（北京） 基于多gpu并行的多分量各向异性叠前时间偏移加速方法
CN107728199B (zh) * 2017-09-22 2019-05-31 中国地质大学（北京） 基于多gpu并行的多分量各向异性叠前时间偏移加速方法
CN110895350A (zh) * 2018-09-13 2020-03-20 中国石油化工股份有限公司 基于gpu的双程波傅里叶有限差分波场传播方法
CN108845355A (zh) * 2018-09-26 2018-11-20 中国矿业大学(北京) 地震偏移成像方法及装置

Similar Documents

Publication Publication Date Title
CN103675895A (zh) 一种利用gpu提高波场延拓计算效率的方法
CN107341544A (zh) 一种基于可分割阵列的可重构加速器及其实现方法
CN103713314B (zh) 一种叠前时间偏移并行处理方法
CN104635258B (zh) 一种应用cpu‑gpu平台进行地震波逆时偏移成像方法
CN103135132A (zh) Cpu/gpu协同并行计算的混合域全波形反演方法
Xu et al. An efficient implementation of 3D high-resolution imaging for large-scale seismic data with GPU/CPU heterogeneous parallel computing
CN103076627A (zh) 一种速度模型平滑优化方法
CN111638551A (zh) 地震初至波走时层析方法及装置
CN104281490A (zh) 一种基于多gpu的高速计算全息图的方法
CN106338766A (zh) 基于分步傅里叶算法的叠前时间偏移方法
Hou et al. Fast inversion of probability tomography with gravity gradiometry data based on hybrid parallel programming
CN104318619B (zh) 面向无损检测的自适应压缩感知的重建方法
US11686870B2 (en) Interpretive-guided velocity modeling seismic imaging method and system, medium and device
WO2020063131A1 (zh) 拟微分算子储存方法及装置
CN116520418A (zh) 一种弹性波角度域共成像点道集高效提取方法
CN102778690B (zh) 一种基于混合基离散傅氏变换的波动方程叠前偏移性能优化方法
Han et al. GPU acceleration of amplitude-preserved Q compensation prestack time migration
CN102279415B (zh) 基于图形处理器计算Fourier积分法单程波深度偏移的方法
Tomas et al. Acceleration of the anisotropic PSPI imaging algorithm with Dataflow Engines
Liu et al. Sparse decomposition of seismic data and migration using Gaussian beams with nonzero initial curvature
Gu et al. Towards Efficient Reverse-time Migration Imaging Computation by Pipeline and Fine-grained Execution Parallelization
Gan et al. Million-core-scalable simulation of the elastic migration algorithm on Sunway TaihuLight supercomputer
Liu et al. A fast and accurate elastic reverse-time migration method based on decoupled elastic wave equations
ZHAO et al. Double‐Grid Algorithm for Calculating Travel Times of Wide‐Angle Reflection Waves
CN107783184B (zh) 一种基于多流优化的gpu逆时偏移方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20140326

RJ01 Rejection of invention patent application after publication