一种利用gpu提高波场延拓计算效率的方法
【专利摘要】本发明提供了一种利用GPU提高波场延拓计算效率的方法，属于地震资料处理领域。所述方法首先扩大震源波场和记录波场的数据网格，然后利用GPU进行多线程并行计算，并应用CUDA提供的二维傅里叶变换函数实现波场延拓。本发明的方法简单实用且效率高，利用本发明得到的成像效果与CPU单程波动方程叠前深度偏移得到的成像效果一致，但是计算效率提高了30倍以上。
【专利说明】一种利用GPU提高波场延拓计算效率的方法
【技术领域】
[0001]本发明属于地震资料处理领域，具体涉及一种利用GPU提高波场延拓计算效率的方法。
【背景技术】
[0002]波动方程叠前深度偏移最核心的工作是地震波场的深度外推。深度外推过程必须依赖于波场外推算子，而波场外推算子实质上是求解单向波动方程。单炮道集波动方程叠前深度偏移成像要分别对炮点波场和检波点波场进行下行波和上行波外推，利用激发时间成像条件提取每一外推层的成像值。
[0003]波动方程叠前深度偏移成像计算是以单炮为一个计算单元，单炮的波场延拓，每外推一步，需要两次二维傅里叶正变换，一次相移，两次二维傅里叶反变换，一次时移，其中二维傅里叶变换计算量最大，因此需要提高其计算效率。
【发明内容】
[0004]本发明的目的在于解决上述现有技术中存在的难题，提供一种利用GPU提高波场延拓计算效率的方法，根据GPU (Graphic Processing Unit图形处理器)的计算特点,应用CUDA (Compute Unified Device Architecture,统一计算设备架构)提供的傅里叶变换函数，有效的提高单程波动方程叠前深度偏移成像的计算效率，缩短单程波动方程叠前深度偏移成像的处理周期，快速地进行波场延拓计算，提高波场延拓的计算效率，提高炮域叠前深度偏移的效率。
[0005]本发明是通过以下技术方案实现的:
[0006]一种利用GPU提高波场延拓计算效率的方法，所述方法首先扩大震源波场和记录波场的数据网格，然后利用GPU进行多线程并行计算，并应用CUDA提供的二维傅里叶变换函数实现波场延拓。
[0007]所述方法包括以下步骤:
[0008](I)扩大震源波场和记录波场的数据网格得到处理后的震源波场和记录波场；
[0009](2)应用裂步傅里叶法(SSF)对所述处理后的震源波场和记录波场进行波场延拓；
[0010](3)设延拓深度iz = I ;
[0011](4)利用CUDA提供的二维傅里叶变换函数对震源波场进行二维正傅里叶变换；
[0012](5)利用CUDA提供的二维傅里叶变换函数对记录波场进行二维正傅里叶变换；
[0013](6)参考速度场的相移计算:利用步骤用(4)和步骤(5)的结果进行相移计算，参考速度场是平均速度，即输入速度场的平均值(请参考Stoffa P L.Split-step Fouriermigration.Geophysics,1990,55(2) ;410 ?421)；
[0014](7)利用CUDA提供的二维傅里叶变换函数对震源波场进行二维反傅里叶变换；
[0015](8)利用CUDA提供的二维傅里叶变换函数对记录波场进行二维反傅里叶变换；[0016](9)扰动速度场的波场计算与相关成像:利用步骤(7)和步骤(8)的结果进行相移计算，扰动速度场是输入速度场与参考速度场的差值(请参考StofTa P L.Split-stepFourier migration.Geophysics,1990,55(2) ;410 ?421)；
[0017](10)判断iz+1是否大于NZ，所述NZ =延拓深度/延拓步长，如果是，则转入步骤
(11)，如果否，则返回步骤⑷；
[0018](11)结束。
[0019]所述步骤(I)是这样实现的:对震源波场和记录波场的长度进行镶边处理，即当震源波场和记录波场的长度小于1000时，将其长度镶边处理为128的倍数，对于不足的数据处填充零；当震源波场和记录波场的长度大于1000而小于2048时，将其长度镶边处理成256的倍数，对于不足的数据处填充零。
[0020]所述步骤⑵至(11)都是在GPU上进行的。
[0021]与现有技术相比，本发明的有益效果是:本发明的方法简单实用且效率高，利用本发明得到的成像效果与CPU单程波动方程叠前深度偏移得到的成像效果一致，但是计算效率提高了 30倍以上。
【专利附图】
【附图说明】
[0022]图1是本发明方法的波场延拓的示意图。
[0023]图2是本发明方法的GPU运算的流程图。
[0024]图3-1是本发明方法实施例中利用CPU计算得到的盐丘模型偏移剖面。
[0025]图3-2是本发明方法实施例中利用GPU计算得到的盐丘模型偏移剖面。
[0026]图4-1是本发明方法实施例中利用CPU计算得到的负向构造偏移剖面。
[0027]图4-2是本发明方法实施例中利用GPU计算得到的负向构造偏移剖面。
[0028]图5是利用CPU和GPU进行波场延拓的计算效率对比图。
[0029]图6是本发明方法的波场延拓的程序框图。
【具体实施方式】
[0030]下面结合附图对本发明作进一步详细描述:
[0031]如图6所示(图6中的NW表示频率片)，本发明利用GPU进行多线程并行计算，应用CUDA (CUDA是一种将GPU作为数据并行计算设备的软硬件体系，采用类C语言进行软件开发。)提供的二维(2D)傅里叶变换函数，结合波场延拓的计算特点(即二维傅立叶计算密集度大)，采用扩大波场数据网格和多信号同时变换的策略，从而达到提高计算效率的目的，该方法简单实用，效率高。
[0032]具体来说，在GPU计算中，内核函数是以线程网格(Grid)的形式组织的，每个线程网格由若干个线程块(Block)组成，而每个线程块又由多个线程(thread)组成。在每次调用GPU计算时，都要指定线程的网格结构〈〈〈dimGrid，dimBloCk>>>。在相移和时移的计算中，根据炮点波场和检波点波场数据体分配线程结构:每个线程计算三维体中的一个点，所有的点可以在由一个线程网格计算完成。因为提高了算法的并行度，所以计算效率大幅度提升。同时，将扰动速度场(时移)计算和相关成像计算进行合并计算，可以减少数据一次读写和一次线程调用，进一步提高了计算效率。[0033]本发明对炮域叠前深度偏移算法中广泛使用的二维傅里叶变换CUFFT (由CUDA数学函数库提供)进行了深入的测试分析，发现CUFFT对信号长度比较敏感，当信号的长度为奇数时，运算速度较慢；当信号的长度为偶数时，速度较快；当信号长度是长度为2的幂次时，速度最快。
[0034]CUFFT的计算速度非常快，计算时间与信号长度基本上呈正比关系。从这个结果出发，本发明对CUFFT的信号长度做镶边处理，以提高CUFFT的计算效率。但是，如果将信号长度都镶边成2n，那么在提高CUFFT计算效率的同时，也会增加数据的存储量，因此要找一个折中点，也就是在稍微增加数据量的基础上大幅提高CUFFT的计算效率；当信号长度小于1000时，信号长度为128的倍数的有较优的计算效率，当信号长度在2048以内，信号长度为256的倍数的有较优的计算效率。因此本发明提出扩大波场数据网格以提高CUFFT的计算效率。
[0035]在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)，扰动速度场(时移)计算+相关成像六步计算，从而提高波场延拓的计算速度。
[0036]用CUDA提供的cufftPlanMany (CUDA傅立叶变换批量处理)进行多信号同时变换(是在(4) (5) (7) (8)这四个步骤中实现的)，具体如下:
[0037]cufftPlanMany (&plan，2，dims, NULL, 1,0, NULL, 1,0, CUFFT C2C, nw)；
[0038]cufftExecC2C(plan, d_signal_s, d_signal_s,CUFFT_F0RWARD)；
[0039]cufftExecC2C(plan, d_signal_r, d_signal_r, CUFFT_F0RWARD)；
[0040]cufftExecC2C(plan, d_signal_s, d_signal_s, CUFFT_INVERSE)；
[0041]cufftExecC2C (plan, d_signal_r, d_signal_r, CUFFT_INVERSE);所述
[0042]图1给出的是波场延拓示意图，单炮道集波动方程叠前深度偏移成像要分别对炮点波场和检波点波场进行下行波和上行波外推，利用激发时间成像条件提取每一外推层的成像值。
[0043]图2给出的是GPU运算最终流程，具体如下:在每个波场延拓的过程中，依次进行以下计算正变换，相移计算，fft反变换，时移+相关成像。
[0044]图3-1是本发明方法实施例中利用CPU计算得到的盐丘模型偏移剖面，图3-2是本发明方法实施例中利用GPU计算得到的盐丘模型偏移剖面，对比图3-1和图3-2可以看出利用CPU和GPU得到的成像结果是一致的。图4-1是本发明方法实施例中利用CPU计算得到的负向构造偏移剖面，图4-2是本发明方法实施例中利用GPU计算得到的负向构造偏移剖面，对比图4-1和图4-2可以看出利用CPU和GPU得到的成像结果是一致的。
[0045]图5表示的是图2的每一步分别在CPU和GPU上的计算时间对比，从图5中可以看出，利用GPU后在计算效率提高了 30倍以上。
[0046]上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述【具体实施方式】所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。
【权利要求】
1.一种利用GPU提高波场延拓计算效率的方法，其特征在于:所述方法首先扩大震源波场和记录波场的数据网格，然后利用GPU进行多线程并行计算，并应用CUDA提供的二维傅里叶变换函数实现波场延拓。
2.根据权利要求1所述的利用GPU提高波场延拓计算效率的方法，其特征在于:所述方法包括以下步骤: (1)扩大震源波场和记录波场的数据网格得到处理后的震源波场和记录波场； (2)应用裂步傅里叶法对所述处理后的震源波场和记录波场进行波场延拓； (3)设延拓深度iz= I ; (4)利用CUDA提供的二维傅里叶变换函数对震源波场进行二维正傅里叶变换； (5)利用CUDA提供的二维傅里叶变换函数对记录波场进行二维正傅里叶变换； (6)参考速度场的相移计算:利用步骤用⑷和步骤(5)的结果进行相移计算，参考速度场是平均速度，即输入速度场的平均值； (7)利用CUDA提供的二维傅里叶变换函数对震源波场进行二维反傅里叶变换； (8)利用CUDA提供的二维傅里叶变换函数对记录波场进行二维反傅里叶变换； (9)扰动速度场的波场计算与相关成像:利用步骤(7)和步骤(8)的结果进行相移计算，扰动速度场是输入速度场与参考速度场的差值； (10)判断iz+1是否大于NZ,所述NZ=延拓深度/延拓步长，如果是，则转入步骤(11)，如果否，则返回步骤(4)； (11)结束。
3.根据权利要求2所述的利用GPU提高波场延拓计算效率的方法，其特征在于:所述步骤(I)是这样实现的:对震源波场和记录波场的长度进行镶边处理，即当震源波场和记录波场的长度小于1000时，将其长度镶边处理为128的倍数，对于不足的数据处填充零；当震源波场和记录波场的长度大于1000而小于2048时，将其长度镶边处理成256的倍数，对于不足的数据处填充零。
4.根据权利要求2所述的利用GPU提高波场延拓计算效率的方法，其特征在于:所述步骤⑵至(11)都是在GPU上进行的。
【文档编号】G01V1/28GK103675895SQ201210315284
【公开日】2014年3月26日 申请日期:2012年8月30日 优先权日:2012年8月30日
【发明者】孔祥宁, 张慧宇, 段心标, 徐兆涛, 张兵, 孙武亮 申请人:中国石油化工股份有限公司, 中国石油化工股份有限公司石油物探技术研究院