专利名称：基于图形处理器计算起伏地表直接叠前逆时偏移的方法
技术领域：
本发明涉及计算计算起伏地表直接叠前逆时偏移的方法，尤其涉及基于图形处理器计算计算起伏地表直接叠前逆时偏移的方法。
背景技术：
随着地震勘探技术的不断发展，油气地震勘探的重点正转向起伏地表条件和复杂地质条件的区域，如山地、滩�：驼釉蟮厍龋舛缘卣鹂碧焦ぷ骷白柿洗硖岢隽诵碌奶粽�。传统地震成像方法在复杂地表、复杂构造地区主要存在两方面的问题一是地表起伏大，表层速度结构十分复杂；二是地下构造复杂，如褶皱强烈、断层发育、构造陡峻、地层产状变化大等。这种双复杂结构使得传统地震成像方法在这类地区难以准确成像。起伏地表地区地震资料成像主要采用两种方法一种是先进行表层波场校正，再偏移成像，该法在实际生产中占主导地位；另一种是直接从起伏地表进行深度域的偏移成像。表层波场校正一直是实际地震资料处理中一项关键处理技术，通常采用高程基准面静校正将地震数据校正到固定基准面或者浮动基准面上的办法解决地形起伏的影响。这种方法隐含着一个明显的基本假设就是地表一致性假设，即在地表起伏不大，低速带横向速度变化缓慢的地区，地下浅、中、深层的反射经过低速带时，几乎遵循同一路径近乎垂直入射至地表，这时它们的静校正量基本相等。在地形平缓的地区，近地表速度比地下速度慢很多，且射线路径出射角较�。捎蒙鲜鲂Ｕ椒ū冉虾鲜剩诟丛拥乇砗透丛拥叵碌厍� 近地表速度和地下速度相差不大，且射线路径出射角较大，因此，简单的校正可能会扭曲波�。档偷卣鸪上竦闹柿浚耸背９娲聿荒懿返某上�。波动方程基准面延拓方法为表层波场校正提供了一种更精确的技术手段。只要能够得到比较准确的近地表速度模型， 以“波动方程基准面延拓”替换常规“时移静校正”，对提高山前带复杂构造成像精度非常重要。直接从起伏地表进行深度域偏移成像的方法，将野外静校正隐含地包括在其中，这种校正不但包含了旅行时的纵向分量，同时还包含了旅行时的横向分量，因此，基于起伏地表的叠前深度偏移可以同时解决上文中提到的复杂地表和复杂地下构造这种地质上的双重复杂性。起伏地表直接叠前偏移对山前带复杂构造成像的价值在于一方面，在近地表反射信息较好的地区，可基于起伏地表叠前深度偏移进行速度分析，自地表逐层向下优化深度一速度模型；另一方面，如果前期工作能够提供满足构造成像要求的宏观速度模型，起伏地表叠前深度一步法偏移则可提供地表以下全面、精细的构造图像。接下来的问题是选择何种偏移算法进行起伏地表直接偏移。KirchhofT积分法简便、高效，适应任意观测方式，且易于对局部目标成像。但Kirchhoff积分法是一种高频近似，不能对复杂波场中的焦散、干涉现象进行很好的处理，在强横向非均勻介质中成像精度不高；与Kirchhoff偏移法相比， 起伏地表单程波波动方程叠前深度偏移法对复杂地区构造成像精度高，且能处理多波至问题，对陡倾界面、断层和断面的刻画更清晰。然而，单程波方法同样是一种近似方法，成像角度限制在90度以内，直接求解双程波波动方程的逆时偏移方法突破了成像倾角的限制，是最精确地成像方法。因此，在复杂地表和复杂地下构造这种双重复杂性的地区应该开展起伏地表直接逆时偏移工作。对于起伏地表双程波波动方程的求解，有限元法处理不规则地表比较方便，但计算效率低，计算精度难以提高，使用也不方便。有限元和其他方法相结合的离散波数和区域分裂等方法，尽管在一定程度上提高了计算效率，但计算精度低、使用不方便等问题并没有根本解决。边界积分或边界元法的半解析性质决定了该方法不能适用于地表速度变化较大情况，限制了其在勘探实际中的应用。
发明内容
针对上述传统方法的缺点，本发明提供了基于图形处理器计算起伏地表直接叠前逆时偏移的方法。在本发明中包括下述两个技术概念，基于图形处理器的计算，所采用的计算方法是有限差分方法。申请人在此作出解释
采用有限差分(FD)方法计算起伏地表直接叠前逆时偏移，此方法计算效率高，使用方便。然而，受算法稳定性和频散关系的影响，有限差分法波场模拟计算量巨大。本发明中波场传播采用的是时间域二阶差分、空间域高阶差分格式。高阶差分格式表明，波场传播和应用成像条件时空间中每个网格点都是解耦的，即独立的。此种方法下所有网格点可以并行计算，并行粒度很小。针对该问题，由于有限差分每个网格点并行粒度很�。夹未砥鞫嘤纱罅康丶扑愫诵淖槌桑虼似浯碓怂闶泵扛鱿叱讨患扑阋桓鐾竦慊蛘呒父鐾竦愕闹担哉庵中×６炔⑿械奈侍庀啾却车闹醒氪砥鞲杏攀�。如附图1所示，显示了本发明的方法算实现逆时偏移的原理流程图，计算密集的波场传播与成像部分借助图形处理器并行计算，中间结果需要借助通讯保存在中央处理器内存或硬盘上。由于运用空间高阶差分法需要大量的内存读写，以三维8阶差分格式为例， 每计算一个网格点的值需要读取网格点周围25个网格点的数据，内存读取冗余度很高。针对这一问题，与传统的计算方法直接使用图形处理器上的内存相比，本发明的方法进一步的通过图形处理器共享存储器(share memory)机制重用内存数据，可以大幅度降低内存读取冗余度，提高算法效率。在上述基础上，有限差分法波场模拟的一个缺陷是处理起伏地表比较困难自由边界处的网格点需要特殊处理，即需要对起伏地表网格点进行细致的分类，在程序实现上需要大量的逻辑判断，这将大大影响图形处理器的运行效率。本发明的方法通过起伏地表边界条件核心函数克服了此一问题。在本发明中，GPU的全称是图形处理器(Graphic Processing Unit)，俗称“显卡”， 在市场上目前有多家公司有研发生产此类硬件，包括NVIDIA、AMD (原ATI)、VIA等x86体系显卡生产商、POWERVR、ARM等基于RISC体系的显卡生产商。本发明公开了基于图形处理器计算起伏地表直接叠前逆时偏移的方法，包括下述步骤
(1)、读取偏移参数及起伏地表边界参数；
(2)、读入单炮数据，并根据观测系统确定偏移孔径；
(3)、确定本炮数据的时间延拓步长；
(4 )、将偏移孔径内的速度模型以及震源子波上传到图形处理器；(5)、利用有限差分的方法正演震源波场到最大接收时间，记录不同时刻的正演波�。�
(6)、将单炮数据读入图形处理器并进行数据规则化；
(7)、利用有限差分的方法沿时间方向反延拓一步，并运用起伏地表边界条件；
(8)、读取步骤(5)中存储的相应时刻的波场并运用互相关成像条件进行成像运算；
(9)、重复步骤(7)、(8)，至反延拓到t=0时刻；
(10)、重复步骤(2)-(9)，至最后一炮数据;
(11)、对偏移结果去除低频成像噪音，输出偏移结果，将偏移结果传输回CPU内存； 其中，上述步骤(5)、(7)、(8)、(9)、(10)是在图形处理器中进行的。 在本发明中，上述步骤(9 )、( 10 )如果在步骤(8 )完成后一次反延拓到t=0或达到最后一炮，则步骤(9)、(10)可省略，对于本领域技术人员而言，这种情况是显而易见可以理解的。在本发明中，步骤(4)、(6)所述的偏移孔径内的速度模型、单炮数据、震源子波传递给图形处理器，是指将偏移孔径内的速度模型、单炮数据、震源子波的数据传递到GPU的内存中(通常称之为显存)。所述的GPU的内存是直接位于GPU芯片上的内存。在本发明中，所述图形处理器尤其可以是近些年广泛兴起的支持通用计算功能的图形处理器，词类处理器通常包括基于统一计算设备架构平台CUDA的图形处理器，此图形处理器由显卡厂商NVIDIA推出；基于开放通用计算模型OpenCL的图形处理器，此类图形处理器包括NVIDIA、AMD、Intel, VIA等均有推出；基于微软DirectCompute标准的图形处理器，此类图形处理器包括NVIDIA、AMD、htel、VIA等均有推出；优选的，基于成熟度和开发快捷的需要，本发明中所用的图形处理器优选为基于统一计算设备架构平台CUDA的图形处理器。本领域技术人员根据本发明的实质精神采用其他类型的支持通用计算功能的图形处理器依旧属于本发明保护范围。基于上述基�。诒痉⒚髦校诓街�(4)、(6)的偏移孔径内的速度模型、震源子波数据、单炮数据由中央处理器通过图形处理器的CUDA编程接口，传送到图形处理器中的渲染管道，到达可编程片段处理器；所述步骤(5)、(7)、(8)、(9)、(10)在图形处理器的多个可编程片段处理器中的多个线程并行进行，通过图形处理器的共享存储器为线程提供、存储、 传输数据，其具体的计算过程如下
步骤(5)调用统一计算设备架构平台CUDA的有限差分核心函数； 步骤(7)调用统一计算设备架构平台CUDA的起伏地表边界条件核心函数； 步骤(8)调用统一计算设备架构平台CUDA的互相关成像条件函数； 在上述基础上，步骤(9)用于判断是否反延拓到t=0时刻，若否，则重复(7)，(8)两步； 若是，则到步骤(10);步骤(10)用于判断是否是最后一炮数据，若否，则重复(2)- (9)步骤；若是，则将偏移结果传输回CPU内存。尽管上述给出的是采用基于统一计算设备架构平台CUDA的图形处理器实现本发明目的的实现方法，然而本领域技术人员根据上述实现方法的实质精神采用其它支持通用计算功能的图形处理器并调用其中对应的函数库或相应函数实现相应功能也是可以为本领域技术人员所能知悉的，依旧属于本发明的保护范围。在上述所述中，所用的词汇“渲染管道、可编程片段处理器”对于本领域技术人员是广为理解的，其他类似的称呼如管线、R0P、处理单元、流处理器等亦可为本领域技术人员所理解为与本发明所用词汇含义一致，改用不同的名称，只要其针对的是图新处理器上相同的硬件单元结构，依旧属于本发明的保护范围。其它的技术用词如共享存储器用于指位于图形处理器上的share memory，在硬件架构中与图形处理器片上缓存Ll Cache类似，所有的共享存储器与Cache相通，用于为Cache提供有效和高速的数据支持；由于有限差分的方法需要大量的内存读写，内存冗余度高，尤其需要高速的数据传输支持，本发明的方法通过使用图形处理器共享存储器重用内存数据，从而大幅度降低内存读取冗余度。与传统方法计算起伏地表逆时偏移相比，本发明的方法具有如下技术优势
1、运用逆时偏移方法直接从起伏地表偏移，与传统的起伏地表处理方法相比较，成像效果优势明显，并且近地表的层位都能够准确的成像，因此在复杂地表和复杂地下构造这种地质上的双重复杂性的区域具有广泛的应用前景。2、针对传统的双程波波动方程求解时计算效率低，计算精度难以提高，使用不方便等缺点，本发明的方法在具体求解上采用有限差分法求解。针对有限差分法波场模拟的处理起伏地表困难，自由边界处的网格点需要特殊处理这一缺陷，本发明在处理自由边界网格点时，采用统一计算设备架构平台CUDA的起伏地表边界条件核心函数将网格点统一处理，既解决了自由边界反射问题，又避免了大量的逻辑判断，从而使得图形处理器得以快速应用。3、本发明的方法利用图形处理器实现计算方法的加速，解决了逆时偏移计算量甚巨的问题。
附图1为本发明基于图形处理器计算起伏地表直接叠前逆时偏移的方法原理图； 附图2为本发明基于图形处理器计算起伏地表直接叠前逆时偏移的方法流程图； 附图3为SEG起伏地表速度模型；
附图4为采用本发明的方法计算基于起伏地表直接叠前逆时偏移结果。
具体实施例方式下面结合附图结合实施例对发明进行进一步阐述，本领域技术人员应知，本发明也可以有多种不同形式或使用不同的图形处理器或者不同的通用计算函数实施，因此不应认为它局限于说明书列出的实施例。本发明的实际应用效果与具体的图形处理器硬件相关，在一定范围内，图形处理器所拥有的流处理器越多、速度越快，速度提升的效果将越好。结合附图1、2，可以理解本发明的基于图形处理器计算起伏地表直接叠前逆时偏移的方法的原理和流程，即利用图形处理器可以同时处理众多并发线程的原理实现高速计算，具体的包含下述步骤
(1)、读取偏移参数及起伏地表边界参数；
(2)、读入单炮数据，并根据观测系统确定偏移孔径；
(3)、确定本炮数据的时间延拓步长；
(4 )、将偏移孔径内的速度模型以及震源子波上传到图形处理器；
(5)、利用有限差分的方法正演震源波场到最大接收时间，记录不同时刻的正演波�。�
(6)、将单炮数据读入图形处理器并进行数据规则化；(7)、利用有限差分的方法沿时间方向反延拓一步，并运用起伏地表边界条件；
(8)、读取步骤(5)中存储的相应时刻的波场并运用互相关成像条件进行成像运算；
(9)、重复步骤(7)、(8)，至反延拓到t=0时刻；
(10)、重复步骤(2)-(9)，至最后一炮数据;
(11)、对偏移结果去除低频成像噪音，输出偏移结果，将偏移结果传输回CPU内存； 其中，上述步骤(5)、(7)、(8)、(9)、(10)是在图形处理器中进行的。附图3为国际上标准的SEG起伏地表模型(Amoco和BP公司设计的加拿大起伏地表逆掩断层模型)数据，其地质原型是加拿大英属哥伦比亚(British Columbia)东北部的逆掩断层构造。该数据共277炮，每炮480道，道距15 m，中间激发，记录长度8s，采样率 4ms,偏移距-3600nT3600m，地形最大高差1527m，速度360(T6000m/s。速度模型剖面长度 25 km，深度 10 km。附图4为采用本发明的方法对原始数据进行起伏地表直接逆时偏移的偏移结果， 如图4所示，逆掩断层成像效果优势明显，并且近地表的层位都能够准确的成像。并且与单纯采用中央处理器的方法相比，速度提升了 100倍以上，具有显著的优势。
权利要求
1.基于图形处理器计算起伏地表直接叠前逆时偏移的方法，包括下述步骤(1)、读取偏移参数及起伏地表边界参数；(2)、读入单炮数据，并根据观测系统确定偏移孔径；(3)、确定本炮数据的时间延拓步长；(4 )、将偏移孔径内的速度模型以及震源子波上传到图形处理器；(5)、利用有限差分的方法正演震源波场到最大接收时间，记录不同时刻的正演波�。�(6)、将单炮数据读入图形处理器并进行数据规则化；(7)、利用有限差分的方法沿时间方向反延拓一步，并运用起伏地表边界条件；(8)、读取步骤(5)中存储的相应时刻的波场并运用互相关成像条件进行成像运算；(9)、重复步骤(7)、(8)，至反延拓到t=0时刻；(10)、重复步骤(2)-(9)，至最后一炮数据；(11)、对偏移结果去除低频成像噪音，输出偏移结果，将偏移结果传输回中央处理器内存；上述步骤(5)、(7)、(8)、(9)、(10)在图形处理器中进行。
2.根据权利要求1所述的基于图形处理器计算起伏地表直接叠前逆时偏移的方法，其特征在于步骤(4)的偏移孔径内的速度模型以及震源子波上传于图形处理器的内存中、步骤(6)的单炮数据读入图形处理器的内存中。
3.根据权利要求1所述的基于图形处理器计算起伏地表直接叠前逆时偏移的方法，其特征在于所述图形处理器是基于统一计算设备架构平台CUDA的图形处理器。
4.根据权利要求3所述的基于图形处理器计算起伏地表直接叠前逆时偏移的方法，其特征在于步骤(4)、(6)的偏移孔径内的速度模型、震源子波数据、单炮数据由中央处理器通过图形处理器的CUDA编程接口，传送到图形处理器中的渲染管道，到达可编程片段处理器；所述步骤(5)、(7)、(8)、(9)、(10)在图形处理器的多个可编程片段处理器的多个线程并行运算，通过图形处理器共享内存为线程提供数据。
5.根据权利要求3所述的基于图形处理器计算起伏地表直接叠前逆时偏移的方法，其特征在于所述步骤(5)调用统一计算设备架构平台CUDA的有限差分核心函数。
6.根据权利要求3所述的基于图形处理器计算起伏地表直接叠前逆时偏移的方法，其特征在于所述步骤(7)调用统一计算设备架构平台CUDA的起伏地表边界条件核心函数。
7.根据权利要求3所述的基于图形处理器计算起伏地表直接叠前逆时偏移的方法，其特征在于所述步骤(8)调用统一计算设备架构平台CUDA的互相关成像条件函数。
8.根据权利要求1所述的基于图形处理器计算起伏地表直接叠前逆时偏移的方法，其特征在于采用有限差分法计算起伏地表直接叠前逆时偏移，其中波场传播采用的是时间域二阶差分、空间域高阶差分格式，波场传播和应用成像条件时空间中每个网格点都是解耦的，独立的。
全文摘要
本发明公开了基于图形处理器计算起伏地表直接叠前逆时偏移的方法，通过使用运用逆时偏移方法直接从起伏地表偏移，并采用有限差分法求解，利用图形处理器的并行计算特性显著提高了计算速度和效率。
文档编号G01V1/28GK102353988SQ201110190520
公开日2012年2月15日 申请日期2011年7月8日 优先权日2011年7月8日
发明者佟小龙, 刘洪 , 刘红伟, 刘钦, 李博 申请人:中国科学院地质与地球物理研究所, 北京吉星吉达科技有限公司
专业数控铣床产品供应商
20年生产经验，实力工厂，精工品质，质量有保障
