CN102636808B - 一种地震叠前时间偏移走时处理方法及装置 - Google Patents

一种地震叠前时间偏移走时处理方法及装置 Download PDF

Info

Publication number
CN102636808B
CN102636808B CN201210026341.7A CN201210026341A CN102636808B CN 102636808 B CN102636808 B CN 102636808B CN 201210026341 A CN201210026341 A CN 201210026341A CN 102636808 B CN102636808 B CN 102636808B
Authority
CN
China
Prior art keywords
walking
point
points
parallel
thread
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.)
Active
Application number
CN201210026341.7A
Other languages
English (en)
Other versions
CN102636808A (zh
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.)
Inspur Beijing Electronic Information Industry Co Ltd
Original Assignee
Inspur Beijing Electronic Information Industry Co Ltd
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 Inspur Beijing Electronic Information Industry Co Ltd filed Critical Inspur Beijing Electronic Information Industry Co Ltd
Priority to CN201210026341.7A priority Critical patent/CN102636808B/zh
Publication of CN102636808A publication Critical patent/CN102636808A/zh
Application granted granted Critical
Publication of CN102636808B publication Critical patent/CN102636808B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Image Processing (AREA)

Abstract

本发明公开了一种地震叠前时间偏移走时处理方法及装置,涉及计算机高性能计算、并行计算领域、石油地震勘探领域。本发明公开的地震叠前时间偏移走时处理方法包括:对成像空间建立三维X-Y-Z坐标系,每个点按照(x,y,z)坐标唯一确定,对X、Y坐标相同的点所构成与Z轴平行的竖线之间并行进行走时计算,对每条竖线内的点并行进行走时计算;建立计算统一设备架构(CUDA)线程模型和CUDA内存模型以进行上述并行走时计算。本发明还公开了一种地震叠前时间偏移走时处理装置。本申请技术方案解决了叠前时间偏移走时计算的瓶颈,令叠前时间偏移处理性能获得大大提升。

Description

一种地震叠前时间偏移走时处理方法及装置
技术领域
本发明涉及计算机高性能计算、并行计算领域、石油地震勘探领域,具体涉及一种针对NVIDIA的Fermi GPU设计的地震叠前时间偏移走时计算并行算法。
背景技术
地震资料处理是石油和天然气勘探开发领域中非常重要的环节。传统的方法是通过人工放炮方式产生地震波,通过地面检波器将地下不同地质层反射回来的地震波信号收集后,利用大型计算机通过多套专业处理软件和一套完整的叠前时间偏移、叠前深度偏移软件系统进行资料处理,从而得到地下的构造以及成像。为石油钻井提供更加可靠的勘探数据,用于勘探专家进行下一步的分析和解释,掌握地下的油气构造。
叠前时间偏移是构造复杂成像最有效的方法之一,它能适应纵向速度变化较大的情况,适用于大倾角的偏移成像。而其中走时计算是计算地震波从震源传到观测点所经过的时间,即计算地震波从震源传到成像点所经过的时间与地震波从成像点传到观测点所走过的时间之和。因此,走时计算是叠前时间偏移中最关键、最费时的一个环节,其采用的算法大致分为三类,即直射线、弯曲射线以及非对称走时计算。由于叠前时间偏移动辄需要处理数以TB的海量数据,以串行方式实现的走时计算算法,处理性能十分低,已经严重制约工业生产。
发明内容
本发明所要解决的技术问题是,提供一种地震叠前时间偏移走时处理方法及装置,以提高叠前时间偏移处理性能。
为了解决上述技术问题,本发明公开了一种地震叠前时间偏移走时处理方法,包括:
对成像空间建立三维X-Y-Z坐标系,每个点按照(x,y,z)坐标唯一确定,对X、Y坐标相同的点所构成与Z轴平行的竖线之间并行进行走时计算,对每条竖线内的点并行进行走时计算;
建立计算统一设备架构(CUDA)线程模型和CUDA内存模型以进行上述并行走时计算。
较佳地,上述方法中,建立CUDA线程模型指:
将整个成像空间XY平面划分为NX*NY/4块,每一个Block的线程计算4条垂直于地面的竖线的点的走时计算,每一个Block包括128线程,其中每32个线程计算一条竖线。
较佳地,上述方法中,所建立的CUDA内存模型包括全局内存、纹理内存、共享内存和常量内存。
较佳地,上述方法中,对每条竖线内的点并行进行走时计算指:
对于每条竖线内的点,先并行计算出索引值能被8整除的点的走时,对不能被8整除的点则依赖于与其临近的已精确计算的两个点的值作线性插值得到该点的走时。
本发明还公开了一种地震叠前时间偏移走时处理装置,包括:
第一单元,对成像空间建立三维X-Y-Z坐标系,每个点按照(x,y,z)坐标唯一确定;
第二单元,对三维X-Y-Z坐标系中X、Y坐标相同的点所构成与Z轴平行的竖线之间并行进行走时计算,对每条竖线内的点并行进行走时计算,以及建立计算统一设备架构(CUDA)线程模型和CUDA内存模型以进行上述并行走时计算。
较佳地,上述装置中,第二单元,建立CUDA线程模型指,将整个成像空间XY平面划分为NX*NY/4块,每一个Block的线程计算4条垂直于地面的竖线的点的走时计算,每一个Block包括128线程,其中每32个线程计算一条竖线。
较佳地,上述装置中,第二单元,所建立的CUDA内存模型包括全局内存、纹理内存、共享内存和常量内存。
较佳地,上述装置中,所述第二单元,对每条竖线内的点并行进行走时计算指:对于每条竖线内的点,先并行计算出索引值能被8整除的点的走时,对不能被8整除的点则依赖于与其临近的已精确计算的两个点的值作线性插值得到该点的走时。
本申请技术方案充分利用了GPU的多核处理能力,解决了叠前时间偏移走时计算的瓶颈,令叠前时间偏移处理性能获得大大提升。即本申请技术方案满足了石油地震勘探资料处理的需求,并且降低了机房构建成本和管理、运行、维护费用。
具体实施方式
图1为本实施例中地震叠前时间偏移走时处理流程图;
图2为本实施例中8层内32个点并行走时计算线程图;
图3为本实施例中1层内32个点并行走时计算线程图;
图4为本实施例中叠前时间偏移串行走时计算成像效果图;
图5为本实施例中叠前时间偏移并行走时计算成像效果图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白,下文将结合附图对本发明技术方案作进一步详细说明。需要说明的是,在不冲突的情况下,本申请的实施例和实施例中的特征可以任意相互组合。
实施例1
本实施例,提供一种地震叠前时间偏移走时处理方法,该方法以Kirchhoff积分法叠前时间偏移走时计算串行算法为基础,并基于NVIDIA的Fermi GPU硬件体系结构,利用CUDA(Computing Unified DeviceArchitecture,计算统一设备架构)技术对原有串行走时计算算法进行并行化设计的。具体地,该方法的过程如图1所示,包括如下步骤:
步骤100,确定并行走时计算算法;
该步骤中,对成像空间建立三维X-Y-Z坐标系,假设第一维大小为NX,第二维大小为NY,第三维大小为NZ,那么需要进行走时计算的总点数为NX*NY*NZ,每个点按照(x,y,z)坐标唯一确定。X、Y坐标相同的点构成一条与Z轴平行的线,物理上是一条垂直于地面的竖线,那么这样的竖线有NX*NY条,并行走时计算采用两级并行,竖线与竖线之间并行进行走时计算;每条竖线内的点并行进行走时计算。
其中,对于同一条竖线上每一个点的走时T可采用多种方式计算。而本实施例中对同一条竖线上每一个点的走时T采用估算方式计算,即首先并行计算出索引值能被8整除的点的走时T,其余点的走时T依赖于与它临近的已精确计算的两个点的值作线性插值得。
步骤200,确定CUDA线程模型;
该步骤中,将整个成像空间XY平面看成一个网格(Grid),将Grid定义为(NX,NY/4),即把整个Grid分为NX*NY/4块(Block),每一个Block的线程计算4条垂直于地面的竖线的点的走时计算;将Block定义为(32,4),即每一个Block包括128线程(Thread),其中每32个线程计算一条竖线。
步骤300,确定CUDA内存模型,主要是根据数据访问特点及NVIDIAFermi GPU内存资源特性,选择不同的内存存放不同的数据,以达到性能最优。
该步骤中,CUDA内存模型至少包括Global memory(全局内存)、Texturememory(纹理内存)、Shared memory(共享内存)和Constant memory(常量内存)。
其中,Global memory指:成像空间数据存储方式是先按照Z方向连续,同时Z方向维度大小总是32的整数倍。Z方向的计算总是以32为单位进行计算。因此在同一时刻,32个线程可以同时访问成像空间的32个点,而且这32个点在内存中是连续的,并且可以做到线程与访问点的一一对应,从而形成对Global Memory的合并访问,提高访存性能;
Texture memory指:由于GPU中的Texture memory有cache,把只读并且频繁访问的大输入道数据存放其中,将提高访存性能;
Shared memory指:由于Shared memory为GPU的片上内存,访问速度快,对于一个Block块中公共的数据,如能被8整除的点的走时、速度场等小数据,可以放入共享内存中,将提高访存性能;
Constant memory指:对于PSTM计算中只读的、频繁被访问的小数据块可以放入Constant memory中,将大大提高访存性能。
下面结合具体应用场景,详细介绍上述方法的实现过程。
首先进行走时计算部分串行算法的并行性分析,将整个走时计算的代码从叠前时间偏移程序中抽取出来,研究走时计算的串行算法,分析其是否具有并行性,研究发现每一点的走时计算是并行的,数据不存在依赖性,可以设计走时计算并行算法,利用CUDA实现,移植到NVIDIA GPU中运行。
其次,需要确定并行走时计算算法,选定一个体偏成像空间,假设其第一维大小为177,第二维大小为1097,第三维大小为1504,那么需要进行走时计算的总点数为177*1097*1504,即垂直于地面的竖线有177*1097=194169条,并行走时计算采用两级并行,上述194169竖线之间互不依赖,可以并行进行走时计算;每条竖线内的1504个点可以并行进行走时计算。并且对于同一条竖线上的1504个点的走时T并不需要精确计算出,而是首先并行计算出索引值能被8整除的点的走时T,其余点的走时T依赖于与它临近的已精确计算的两个点的值作线性插值得。其具体过程如下:
A、如图2所示,以8层为一个计算单元,1个warp内32个线程精确并行计算出索引值能被8整除的33个点的走时T,其中线程0多计算一个点;
B、判断一层内所有点的索引号落在哪两个最近的索引号能被8整除的点之间,然后计算其插值因子;
C、并行进行插值计算,如图3所示,以一层为一个计算单元,1个warp内32个线程并行插值计算一层内所有点(包括索引号能被8整除的点和索引号不能被8整除的点)的走时T,循环8次后,整个8层的所有点的走时T将被计算完成。
接下来,确定CUDA线程模型,将整个上述成像空间XY平面看成一个网格(Grid),其Grid定义为dim3 Grid(177,(1097+3)/4),即把整个Grid分为177*(1097+3)/4块(Block),每一个Block的线程计算4条垂直于地面的竖线的点的走时计算;把Block定义为dim3 Block(32,4),即每一个Block包括128线程(Thread),其中每32个线程计算一条竖线。
最后确定CUDA内存,即根据数据访问特点及NVIDIA Fermi GPU内存资源特性,选择不同的内存存放不同的数据,以达到性能最优。其中,CUDA内存至少包括Global memory、Texture memory、Shared memory和Constantmemory。
下面针对上述具体计算,进行正确性测试,主要通过运行上述体偏作业,获得图4、图5成像效果图。其中,图4为叠前时间偏移走时计算采用串行算法,运行于CPU上所获得的成像效果。图5为叠前时间偏移走时计算采用并行算法,运行于GPU上所获得的成像效果。从结果来看,两幅图像不存在明显差异,证明基于GPU的地震叠前时间偏移走时处理方法是正确的。
再来了解性能测试中的测试环境及测试数据。测试环境包括硬件环境、软件环境、叠前时间偏移运行内核,其中叠前时间偏移运行内核为叠前时间偏移的核心计算部分,即包括走时计算部分;测试数据为输入的测试地震道数据集。对于成像空间而言,其中第一维为X轴方向的大小;第二维为Y轴方向的大小;第三维为Z轴方向的大小,具体各项参数如表1所示。
表1为测试环境和测试数据参数表
最后看测试结果,为了保证测试性能结果的稳定性,对上述体偏作业进行了10次测试,运行于CPU上的串行走时计算10次的平均时间是54320秒,而对于上述同样体偏作业而言,运行GPU上并行走时计算10次的平均时间是1646秒,并行走时计算算法的性能是串行走时计算算法的54320/1646=33倍。
由上述实施例可以看出,基于GPU的地震叠前时间偏移走时处理方案是正确的,成像效果跟CPU串行走时计算算法一样,并且其性能是串行走时计算算法的33倍,这样不但满足了石油地震勘探处理的需求,而且大大降低了功耗,减少了机房构建成本和管理、运行、维护费用,而且这种方法实现简单,需要的开发成本低。
实施例2
本实施例介绍一种地震叠前时间偏移走时处理装置,可实现上述实施例1的方法。该装置包括如下各单元。
第一单元,对成像空间建立三维X-Y-Z坐标系,每个点按照(x,y,z)坐标唯一确定;
第二单元,对三维X-Y-Z坐标系中X、Y坐标相同的点所构成与Z轴平行的竖线之间并行进行走时计算,对每条竖线内的点并行进行走时计算,以及建立计算统一设备架构(CUDA)的线程模型和CUDA内存模型以进行上述并行走时计算。需要说明的是第二单元所建立的CUDA内存模型包括全局内存、纹理内存、共享内存和常量内存。
上述第二单元,建立CUDA线程模型指,将整个成像空间XY平面划分为NX*NY/4块,每一个Block的线程计算4条垂直于地面的竖线的点的走时计算,每一个Block包括128线程,其中每32个线程计算一条竖线。
而上述第二单元,对每条竖线内的点并行进行走时计算指:对于每条竖线内的点,先并行计算出索引值能被8整除的点的走时,对不能被8整除的点则依赖于与其临近的已精确计算的两个点的值作线性插值得到该点的走时。
其他细节可参见实施例1,在此不再赘述。
本领域普通技术人员可以理解上述方法中的全部或部分步骤可通过程序来指令相关硬件完成,所述程序可以存储于计算机可读存储介质中,如只读存储器、磁盘或光盘等。可选地,上述实施例的全部或部分步骤也可以使用一个或多个集成电路来实现。相应地,上述实施例中的各模块/单元可以采用硬件的形式实现,也可以采用软件功能模块的形式实现。本申请不限制于任何特定形式的硬件和软件的结合。
以上所述,仅为本发明的较佳实例而已,并非用于限定本发明的保护范围。凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (6)

1.一种地震叠前时间偏移走时处理方法,其特征在于,该方法包括:
对成像空间建立三维X-Y-Z坐标系,每个点按照(x,y,z)坐标唯一确定,对X、Y坐标相同的点所构成与Z轴平行的竖线之间并行进行走时计算,对每条竖线内的点并行进行走时计算;
建立计算统一设备架构CUDA线程模型和CUDA内存模型以进行上述并行走时计算;
其中,对每条竖线内的点并行进行走时计算指:
对于每条竖线内的点,先并行计算出索引值能被8整除的点的走时,对不能被8整除的点则依赖于与其临近的已精确计算的两个点的值作线性插值得到该点的走时,具体包括:
A、以8层为一个计算单元,1个线程组织单元warp内32个线程精确并行计算出索引值能被8整除的33个点的走时T,其中线程0多计算一个点;
B、判断一层内所有点的索引号落在哪两个最近的索引号能被8整除的点之间,然后计算其插值因子;
C、并行进行插值计算,以一层为一个计算单元,1个线程组织单元warp内32个线程并行插值计算一层内所有点的走时T,循环8次后,整个8层的所有点的走时T将被计算完成,其中,所述所有点包括索引号能被8整除的点和索引号不能被8整除的点。
2.如权利要求1所述的方法,其特征在于,建立计算统一设备架构CUDA线程模型指:
将整个成像空间XY平面划分为NX*NY/4块,每一个线程块Block的线程计算4条垂直于地面的竖线的点的走时计算,每一个Block包括128线程,其中每32个线程计算一条竖线;
其中,所述NX*NY是指:对成像空间建立三维X-Y-Z坐标系,第一维大小为NX,第二维大小为NY,第三维大小为NZ,那么需要进行走时计算的总点数为NX*NY*NZ,每个点按照(x,y,z)坐标唯一确定;X、Y坐标相同的点构成一条与Z轴平行的线,物理上是一条垂直于地面的竖线,那么这样的竖线有NX*NY条;
所述NX*NY/4是指:对于实际需要处理的所述NX*NY条测线,CUDA线程模型中需要的线程块数为NX*NY/4。
3.如权利要求1或2所述的方法,其特征在于,
所建立的CUDA内存模型包括全局内存、纹理内存、共享内存和常量内存。
4.一种地震叠前时间偏移走时处理装置,其特征在于,该装置包括:
第一单元,对成像空间建立三维X-Y-Z坐标系,每个点按照(x,y,z)坐标唯一确定;
第二单元,对三维X-Y-Z坐标系中X、Y坐标相同的点所构成与Z轴平行的竖线之间并行进行走时计算,对每条竖线内的点并行进行走时计算,以及建立计算统一设备架构CUDA线程模型和CUDA内存模型以进行上述并行走时计算;
其中,所述第二单元,对每条竖线内的点并行进行走时计算指:对于每条竖线内的点,先并行计算出索引值能被8整除的点的走时,对不能被8整除的点则依赖于与其临近的已精确计算的两个点的值作线性插值得到该点的走时,具体包括:
A、以8层为一个计算单元,1个线程组织单元warp内32个线程精确并行计算出索引值能被8整除的33个点的走时T,其中线程0多计算一个点;
B、判断一层内所有点的索引号落在哪两个最近的索引号能被8整除的点之间,然后计算其插值因子;
C、并行进行插值计算,以一层为一个计算单元,1个warp内32个线程并行插值计算一层内所有点的走时T,循环8次后,整个8层的所有点的走时T将被计算完成,其中,所述所有点包括索引号能被8整除的点和索引号不能被8整除的点。
5.如权利要求4所述的装置,其特征在于,
第二单元,建立计算统一设备架构CUDA线程模型指,将整个成像空间XY平面划分为NX*NY/4块,每一个线程块Block的线程计算4条垂直于地面的竖线的点的走时计算,每一个Block包括128线程,其中每32个线程计算一条竖线;
其中,所述NX*NY是指:对成像空间建立三维X-Y-Z坐标系,第一维大小为NX,第二维大小为NY,第三维大小为NZ,那么需要进行走时计算的总点数为NX*NY*NZ,每个点按照(x,y,z)坐标唯一确定;X、Y坐标相同的点构成一条与Z轴平行的线,物理上是一条垂直于地面的竖线,那么这样的竖线有NX*NY条;
所述NX*NY/4是指:对于实际需要处理的所述NX*NY条测线,CUDA线程模型中需要的线程块数为NX*NY/4。
6.如权利要求4或5所述的装置,其特征在于,
第二单元,所建立的CUDA内存模型包括全局内存、纹理内存、共享内存和常量内存。
CN201210026341.7A 2012-02-07 2012-02-07 一种地震叠前时间偏移走时处理方法及装置 Active CN102636808B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210026341.7A CN102636808B (zh) 2012-02-07 2012-02-07 一种地震叠前时间偏移走时处理方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210026341.7A CN102636808B (zh) 2012-02-07 2012-02-07 一种地震叠前时间偏移走时处理方法及装置

Publications (2)

Publication Number Publication Date
CN102636808A CN102636808A (zh) 2012-08-15
CN102636808B true CN102636808B (zh) 2014-11-05

Family

ID=46621270

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210026341.7A Active CN102636808B (zh) 2012-02-07 2012-02-07 一种地震叠前时间偏移走时处理方法及装置

Country Status (1)

Country Link
CN (1) CN102636808B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105572730B (zh) * 2015-12-15 2017-11-14 中国科学院地质与地球物理研究所 三维复杂结构声波正演方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102243321A (zh) * 2011-03-15 2011-11-16 浪潮(北京)电子信息产业有限公司 一种地震叠前时间偏移的处理方法及系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6912468B2 (en) * 2003-08-14 2005-06-28 Westerngeco, L.L.C. Method and apparatus for contemporaneous utilization of a higher order probe in pre-stack and post-stack seismic domains

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102243321A (zh) * 2011-03-15 2011-11-16 浪潮(北京)电子信息产业有限公司 一种地震叠前时间偏移的处理方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于GPU实现叠前时间偏移走时计算的并行算法;张清等;《计算机系统应用》;20111231;第20 卷(第8期);42-46 *
张清等.基于GPU实现叠前时间偏移走时计算的并行算法.《计算机系统应用》.2011,第20 卷(第8期),42-46. *

Also Published As

Publication number Publication date
CN102636808A (zh) 2012-08-15

Similar Documents

Publication Publication Date Title
CN109800863B (zh) 一种基于模糊理论和神经网络的测井相识别方法
CN105334542B (zh) 任意密度分布复杂地质体重力场快速、高精度正演方法
US9291735B2 (en) Probablistic subsurface modeling for improved drill control and real-time correction
CN104635262B (zh) 一种基于增强型矩形网格的正逆断层等值线自动生成方法
CN103713314B (zh) 一种叠前时间偏移并行处理方法
CN105005072B (zh) 利用cuda的pml边界三维地震波传播模拟方法
CN105954802A (zh) 一种岩性数据体的转换方法及装置
CN109100795A (zh) 一种面元的炮检点布设方法、装置及系统
Noble et al. High-performance 3D first-arrival traveltime tomography
US20160377752A1 (en) Method of Digitally Identifying Structural Traps
CN105894439A (zh) 基于CUDA的海洋涡旋及Argo浮标交集数据快速提取算法
CN105116447A (zh) 一种基于曲率异常条带的地质河道方向判别方法
CN102636808B (zh) 一种地震叠前时间偏移走时处理方法及装置
CN111413741A (zh) 一种砂岩型铀矿资源量计算方法和装置
CN103886129B (zh) 将测井数据离散到储层网格模型的方法和装置
CN109782355A (zh) Obs检波点漂移的检测方法及装置
CN107526104A (zh) 基于多机多核的裂缝介质地震波场数值模拟方法
CN104422953A (zh) 一种提高地震叠前时间偏移计算效率的方法
CN107797148B (zh) 一种基于三维地质建模的航磁异常场分离方法及系统
Weinbub et al. Shared-memory parallelization of the fast marching method using an overlapping domain-decomposition approach
US10705235B2 (en) Method of characterising a subsurface volume
CN109752758B (zh) 一种地震数据分解方法、系统及存储介质与终端
CN110967754B (zh) 一种基于偏移速度寻优的缝洞储层充填与流体识别方法
CN104407383B (zh) 三维地震照明分析并行实现方法
CN112363221B (zh) Walkaway VSP测线的布设方法及装置

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant