边界显存不足

第53卷第12期2010年12月地球物理学报CHINESEJOURNALOFGEOPHYSICSVol.
53,No.
12Dec.
,2010李博,刘红伟,刘国峰等.
地震叠前逆时偏移算法的CPU/GPU实施对策.
地球物理学报,2010,53(12):2938~2943,DOI:10.
3969/j.
issn.
00015733.
2010.
12.
017LiB,LiuHW,LiuGF,etal.
ComputationalstrategyofseismicprestackreversetimemigrationonCPU/GPU.
犆犺犻狀犲狊犲犑.
犌犲狅狆犺狔狊.
(inChinese),2010,53(12):2938~2943,DOI:10.
3969/j.
issn.
00015733.
2010.
12.
017地震叠前逆时偏移算法的犆犘犝/犌犘犝实施对策李博1,刘红伟1,刘国峰1,佟小龙2,刘洪1,郭建1,裴江云31中国科学院地质与地球物理研究所,北京1000292北京吉星吉达科技有限公司,北京1000293中石油大庆油田有限责任公司勘探开发研究院,大庆163712摘要相较于单程波偏移算法而言,逆时偏移成像方法以其物理基础为依托优势,几十年来一直备受国内外地球物理学家的青睐.
目前的逆时偏移(RTM)若直接采用双程波动方程进行延拓,尽管可以回避上下行波的分离处理,然就已有算法而言,其计算量和I/O(输入/输出)量却是最大的.
针对此问题,本文在分析现行逆时偏移的多种算法基础上,提出利用CPU/GPU(中央处理器/图形处理器)作为数值计算核心,建立随机边界模型,从而克服存储I/O难题和提高计算效率.
在实际的数据测试中,本文的方法可以大幅度的提高计算效率和减少存储单元,从而促使其高效地应用于生产实际.
关键词逆时偏移,波动方程,随机边界,中央处理器,图形处理器DOI:10.
3969/j.
issn.
00015733.
2010.
12.
017中图分类号P631收稿日期20100208,20101125收修定稿基金项目国家重大油气专项(2008ZX05008006031),国家重大专项(2008ZX05023005009),国家重大科研装备研制项目(O823011119)资助.
作者简介李博,男,1982年生,中国科学院地质与地球物理研究所在读博士,从事地震数据叠前成像相关研究工作.
Email:libo@mail.
iggcas.
ac.
cn通信作者刘洪,研究员.
Email:liuhong@mail.
igcas.
ac.
cn犆狅犿狆狌狋犪狋犻狅狀犪犾狊狋狉犪狋犲犵狔狅犳狊犲犻狊犿犻犮狆狉犲狊狋犪犮犽狉犲狏犲狉狊犲狋犻犿犲犿犻犵狉犪狋犻狅狀狅狀犆犘犝/犌犘犝LIBo1,LIUHongWei1,LIUGuoFeng1,TONGXiaoLong2,LIUHong1,GUOJian1,PEIJiangYun31犐狀狊狋犻狋狌狋犲狅犳犌犲狅犾狅犵狔犪狀犱犌犲狅狆犺狔狊犻犮狊,犆犺犻狀犲狊犲犃犮犪犱犲犿狔狅犳犛犮犻犲狀犮犲狊,犅犲犻犼犻狀犵100029,犆犺犻狀犪2犅犲犻犼犻狀犵犌犲狅狊狋犪狉犛犮犻犲狀犮犲犪狀犱犜犲犮犺狀狅犾狅犵狔犔狋犱.
犅犲犻犼犻狀犵100029,犆犺犻狀犪3犈狓狆犾狅狉犪狋犻狅狀犪狀犱犇犲狏犲犾狅狆犿犲狀狋犐狀狊狋犻狋狌狋犲,犇犪狇犻狀犵犗犻犾犉犻犲犾犱犔狋犱,犘犲狋狉狅犆犺犻狀犪,犇犪狇犻狀犵163712,犆犺犻狀犪犃犫狊狋狉犪犮狋Comparingwithonewaywavemigrationalgorithm,reversetimemigration(RTM)ismoreattractivebecauseofthetheoryadvantages.
TwowaywaveequationhasbeenusedtoextrapolatewavefieldinRTM,insteadofseparatingtheupgoingwaveanddowngoingwave.
However,duetothelargeamountofcomputationandI/O,RTMismosttimeconsuminginindustrialapplications.
Inthisarticle,weanalyzeseveralcomputationalstrategiesandproposeourmethod,whichusesCPU/GPUascomputationalcoreandbuildsrandomvelocityboundary,forsolvingI/Oproblemandcomputationalefficiencyproblem.
Intheactualtest,ithasbeenprovedthatthismethodcanlargelydecreasestoragememoryunitsandimprovecomputationalefficiency.
犓犲狔狑狅狉犱狊Reversetimemigration(RTM),Waveequation,Randomboundarycondition,Centralprocessingunit(CPU),Graphicprocessingunit(GPU)12期李博等:地震叠前逆时偏移算法的CPU/GPU实施对策1引言逆时偏移的思想早已由Hemon在1978年提出[1],继而,由Baysal、Whitemore、McMechan等[2~4]将之应用于纵波资料的叠后偏移,他们在获得好效果的同时,明确告诫尚存在艰难的数值计算效率和存储等关键问题待解决.
当时的逆时偏移方法虽然只能限于在2D或者叠后偏移中应用,却引发了业界的巨大兴趣.
逆时偏移方法是通过双程波波动方程在时间域上对人工给予的震源子波正向传播和接收到的地震资料进行反向传播,并结合成像条件实现偏移[5].
一方面逆时偏移需同时计算震源在正向传播和地震资料的反向传播过程中的各个时刻的波场,因其各自时间的延拓方向不同,致使实际应用中必须存储其中一个方向的传播过程,这就是何以需要巨大额外存储空间.
如今的地震资料采集已经发展到了三维甚至四维,其数据量更是十分庞大,逆时偏移在推向实际应用时尚需面对此类困难.
另一方面,相较于单程波方程的波场延拓而言,逆时偏移运用双程波波动方程进行波场延拓,避免了上下行波的分离处理,因而成为最准确的成像算法,且不受倾角的限制,并能实现回转波和多次波的成像.
正基于此,研究的焦点也集中于提高逆时偏移的计算效率和减少偏移中的存储I/O问题[6].
近年来,随着计算机技术的发展和计算能力的大幅度提高,业界在逆时偏移研究工作中不断取得突破性的成果.
例如各向异性的逆时偏移、弹性波逆时偏移以及逆时偏移的噪音衰减技术等[7,8].
2008年以来,快速发展的GPU(GraphicProcessingUnit)以及其编程架构CUDA(ComputingUnifiedDeviceArchitecture)已为深度偏移的快速计算提供更加有力的工具.
本文着重由声波方程的高阶有限差分近似为切入点,分析现行多种逆时偏移的实现策略,提出一种利用随机边界的CPU/GPU为编程环境的地震逆时偏移算法实施的对策,并经数值模型和实际资料验证其有效性.
2波动方程的高阶有限差分近似三维声波方程的表达如下:1犞22犘狋2=2犘狓2+2犘狔2+2犘狕2,(1)其中狋代表时间坐标,狓,狔,狕代表空间坐标,犘表示位移函数,犞表示速度函数.
利用时间二阶中心有限差分近似2犘狋2,利用空间高阶中心有限差分近似2犘狓2+2犘狔2+2犘狕2可得:犘狋+Δ狋狓,狔,狕-2犘狋狓,狔,狕+犘狋-Δ狋狓,狔,狕=Δ狋2犞2[犔2(狓)+犔2(狔)+犔2(狕)]犘狋狓,狔,狕,(2)其中,犔(狓)=狓,犔(狔)=狔,犔(狕)=狕.
利用犖阶中心差分近似得:犔2(狓)犘狋狓,狔,狕=∑犖犾=-犖犪犾Δ狓2犘狋狓+犾,狔,狕,犔2(狔)犘狋狓,狔,狕=∑犖犾=-犖犪犾Δ狔2犘狋狓,狔+犾,狕,犔2(狕)犘狋狓,狔,狕=∑犖犾=-犖犪犾Δ狕2犘狋狓,狔,狕+犾.
(3)其中,犪犾=(-1)犾+1犾2(-1)犾+1∏犖犻=1,犻≠犾犻2∏犾-1犻=1,犻≠犾(犾2-犻2)∏犖犻=犾+1(犻2-犾2),犾=1,2,…,犖犪犾=犪-犾,犪0=-2∑犖犾=1犪犾.
分别用Δ狓,Δ狔,Δ狕表示差分网格的间距,犺表示最大网格间距,狀表示一个波形中的采样点数,Δ狋表示延拓步长,犳max表示子波频率的最大值,犞min表示速度模型的最小值,那么差分格式的频散条件可表示为犺≤犞min狀犳max,犺=max(Δ狓,Δ狔,Δ狕),(4)稳定性条件:Δ狋犞1Δ狓2+1Δ狔2+1Δ狕槡2≤1∑犖犾=1犪2犾-烄烆烌烎11/2.
(5)3成像条件成像条件的制定是地震偏移成像算法的关键之一,它直接影响成像的效果和计算成本,本文应用Claerbout的互相关成像原理.
互相关成像条件的成像公式为9392地球物理学报(ChineseJ.
Geophys.
)53卷犐(狓,狕)=∫狋max0犝(狓,狕,狋)犇(狓,狕,狋)d狋,(6)其中,犝(狓,狕,狋)为上行波波场,犇(狓,狕,狋)为下行波波场,d狋为延拓时间步长.
双程波逆时偏移中以外推观测波场犘R(狓,狕,狋)取代上行波,以震源外推波场犘S(狓,狕,狋)取代下行波,(6)式变为犐(狓,狕)=∫狋max0犘R(狓,狕,狋)犘S(狓,狕,狋)d狋.
(7)(7)式中的被积函数犘R(狓,狕,狋)犘S(狓,狕,狋)表示狋时刻对波场做一次成像运算,积分说明如空间犐(狓,狕)中的像是多个时间步成像的叠加.
所以,互相关成像条件充分利用了成像信息,在增强成像信号的同时也有效压制了成像噪音.
逆时偏移互相关成像条件的成像步骤是:(1)首先要完成一次波动方程正演计算,并保存每一时间步的波场信息;(2)逆时外推记录波场,存储每一时间步的波场信息;(3)分别读取保存的同一时刻的震源波场和记录波场做成像运算,即犘R(狓,狕,狋)犘S(狓,狕,狋),再累加入成像空间.
4算法实现策略分析逆时偏移的成像条件需要使用在同一时刻的震源波场(经过激发和正演模拟波的传播)和记录波场(逆时反传播回地下的波场).
因为一个是正传波场,另一个是反传波场,如果想同时得到相同时刻的两个波场,则必须存储其中一个波场的传播过程,即每一时刻的波场分布,需要消耗甚大的存储资源,这个要求在实际操作中是难以满足的.
为此,许多学者找到一些解决办法.
其中包括:策略一震源和记录同时从犜狀(表示狀d狋时刻)延拓到犜0,这样就不需要额外的存储空间,并且可以达到最小计算量犗(2犖)的量级,然而犜狀时刻的震源波场并不是已知数据,所以不能直接实现.
策略二存储犜0到犜狀的波场,然后反传记录波场,当传播的某个时刻时读取那个时刻的激发点波场.
此方法也可以达到犗(2犖)计算量的要求,但是需要巨大的磁盘空间来存储波场,即使磁盘空间可以满足需求,但是因此产生的I/O时间消耗十分庞大.
策略三只存储少数几个时刻的震源波场,在记录波场反向延拓的过程中,利用这些存储的震源波场插值来近似当前时刻的波场然后成像.
这种方法对波场的插值首先是不准确的,并且计算量也随之增加.
策略四首先把震源波场正向延拓到犜狀和犜狀-1两个时刻,把这两个波场作为初始条件,与记录波场同时反向延拓,并且随延拓过程中应用成像条件,不需要额外的存储空间.
这种策略只可以应用在密度为常数的介质中,并且激发点边界条件是Dirichlet边界条件.
但是这种波场的延拓假设是在无限大空间内传播的,而通常的计算方法是加入了人工边界的,包括自由边界条件、吸收边界条件、指数衰减边界条件和完美匹配层边界条件(PML)等[9~11].
其中除自由边界条件外,其他三种边界条件都会破坏波场的完整性,使得波场的延拓不可逆.
为此只能在适当的时刻在适当的空间位置上插入波场能量以补充波场的完整性[12,13],这样就需要存储每个时刻的位于人工边界处的波场值,当计算规模增大后,存储量也会快速上升,尽管小于策略二中需要的存储量,但是也是通常不能满足实际的要求.
策略五首先正向延拓波场,并且记录多组检查点(犘狀,犘狀-1)作为中间时刻的波场初始条件,然后反向延拓检波点波场到某时刻,用最近的检查点来正向传播得到该时刻的震源波场,然后应用成像条件,当应用PML边界条件的时候波场还需要额外的存储空间,并且需要大量的重复计算过程,最大的重复计算可能会达到犗(犿犖)的量级,其中犿为检查点的个数.
Griewank利用这个思路提出优化的检查点的选取方法,可以使得重复计算的次数减少到犗(犖log犿)的量级[14].
尽管如此,逆时偏移的消耗时间和存储的矛盾也并未得到完全解决,但是该方法是目前最经济的计算策略.
在上述5种逆时偏移的计算策略中可以看出计算量级与存储规模是一对矛盾,若要得到最小的计算量就会达到最大的存储量,反之亦然.
然而在策略四中不需要额外的存储空间,而其他策略都要涉及到存储规模与计算量级的平衡问题.
在策略四中,人工边界是主要的影响因素,然而自由边界条件又会引入很多的成像假象,其他的边界条件又不符合该策略的要求只能舍弃.
针对策略四,本文依据2009年Robert提出的随机速度边界模型提出一种利用GPU优化的逆时偏移计算方法.
5随机边界的建立2009年Robert[15]提出随机边界模型,其思想049212期李博等:地震叠前逆时偏移算法的CPU/GPU实施对策是消除人工边界自由边界条件反射波的相干性,使边界反射不能成像.
具体实现过程为将有限空间外扩一定的距离,然后在外扩的空间内填充随机速度,从而形成随机边界速度模型(如图1b),当波传到随机速度区域时波前面将被随机化,使波场变成随机噪音反传回真实速度区域,破坏了边界反射的相干性使得边界反射不能成像.
本文构造随机边界函数如下:珟犞(狓,狕)=犞(狓,狕)-狉犱,(8)其中珟犞(狓,狕)为边界点的随机速度函数,犞(狓,狕)为边界点的原始速度函数,狉是随机数,犱为速度点与内层边界的空间距离.
运用CPU/GPU作为偏移计算核心,采用随机速度边界模型,将激发点波场传播至犜狀时刻,然后利用犜狀时刻的波场作为初始条件,同时反传激发点波场和检波点波场,并且同时利用成像条件进行成像.
进而避免了额外的存储空间,虽然激发点的波场重复计算一次,但利用GPU作为波场延拓的数值计算核心其耗时与大规模的磁盘I/O相比还是非常经济的.
另外,从文后的实际测试中也可表明,利用了GPU作为波场延拓函数的计算核心,延拓过程的计算已经不是主要瓶颈.
存储检查点的计算策略五是目前认为最经济的计算策略,所以利用该策略与本文的计算策略做对比.
6数据试算为了验证上述方法的有效性,我们设计了一个空间采样为256*200的网格(如图1所示),取边界宽度犖犔=30,网格间距为10m,在网格的三个边界生成随机速度模型,在模型的中间放置震源,然后记录不同时刻的波前快照.
图1原始速度模型(a)及随机速度边界模型(b)Fig.
1Originalvelocitymodel(a)andrandomvelocityboundarymodel(b)计算结果如图2所示,波前面正传到120*4ms时的波场尚未到达模型的边界(图2a),正传到180*4ms时的波场已经到达随机边界并且有部分反射能量(图2b),但是与自由边界条件(图2e)相比在同一时刻的反射能量的相干性已经削弱很多.
当正传到400*4ms时的波场已经完全成为随机噪音.
我们将400*4ms时的波场作为初始值再反传到120*4ms时刻,其波场示于图2d,利用正传波场与反传波场在120*4ms时刻的差示于图2f,其误差值量级在10-7以内,与单精度浮点运算的精度范围相同,所以我们认为这样的误差是可以接受的.
7实际数据计算利用常用的BP2004模型作为实例,运用本文提供的方法做数值试验,分别对该模型利用自由边界条件(FBC)、吸收边界条件(ABC)和本文的随机边界条件(RBC)做偏移测试,其中自由边界采用成像策略四,吸收边界采用成像策略五,随机边界采用本文的成像策略.
使用吸收边界条件时的时间延拓步长与其他两种方法相同,但成像时间步长数据的采样率为4ms,这样使得内存的使用量可以达到最1492地球物理学报(ChineseJ.
Geophys.
)53卷图2随机边界模型的波场正演快照(a)正传到120*4ms时的波场;(b)正传到180*4ms时的波场;(c)正传到400*4ms时的波场;(d)反传到120*4ms时的波场;(e)应用自由边界条件时正传到180*4ms时的波场;(f)正传波场(a)与反传波场(d)在同一时刻的差.
Fig.
2randomboundarymodelforwardpropagationwavefieldsnapshot(a)Forwardwavefieldat120*4ms;(b)Forwardwavefieldat180*4ms;(c)Forwardwavefieldat400*4ms;(d)Backwardwavefieldat120*4ms;(e)Forwardwavefieldat120*4mswithfreereflectionboundarycondition;(f)Differencebetween(a)and(d).
图3三种边界条件的偏移结果(a)BP2004速度模型;(b)自由边界偏移结果;(c)随机边界偏移结果;(d)吸收边界偏移结果.
Fig.
3Threeboundaryconditionmigrationresult(a)BP2004velocitymodel;(b)FBCmigrationresult;(c)RBCmigrationresult;(d)ABCmigrationresult.
249212期李博等:地震叠前逆时偏移算法的CPU/GPU实施对策表1计算用时和存储量犜犪犫犾犲1犆狅犿狆狌狋犻狀犵狋犻犿犲犪狀犱犿犲犿狅狉狔狉犲狇狌犻狉犲犿犲狀狋狊计算核心与边界条件单炮偏移所需时间1348炮偏移所用时间所需内存使用量CPUABC:60min80880min>8GBCPU/GPUABC:100s153069s>8GBCPU/GPURBC:90s121358s383MBCPU/GPUFBC:88s118654s383MB小.
为了不失公平性,偏移其他参数完全相同.
偏移结果示于图3,计算所需时间和存储需求示于表1.
从图3中的计算结果来看,随机边界的成像质量大大优于自由边界条件的成像效果,其原因是显而易见的.
随机边界的成像质量略优于吸收边界,笔者认为,本文采用的吸收边界为1977年Clayton&Engquist提出的单程波方程的二阶旁轴近似,其吸收效果并不是十分理想,仍然有边界反射能量反传回来而造成成像结果较差,若采用完美匹配边界条件(PML)可能效果会有所改善.
随机边界条件使得反射波成为随机噪音反传回来,经过多炮的随机叠加有消除噪音的作用.
从表1中可见,若利用CPU作为计算核心,完成逆时偏移的耗时非常巨大,而利用CPU/GPU联合作为计算核心,计算耗时已经大幅度降低,这样更加方便应用于实际当中.
若利用CPU/GPU联合作为计算核心,不同边界条件的耗时基本相当,但是吸收边界仍然需要大量的存储空间,而随机边界条件和自由边界条件都只需要很少的内存空间即可完成,并且而随机边界用时要少于吸收边界的计算时间,略大于自由边界条件的计算用时.
笔者认为,随机边界尽管要重复计算一次震源波场的逆向传播,但是相比于吸收边界的存储IO时间来说还是值得的.
8结语本文利用CPU/GPU作为计算核心,并将计算量最大的波场延拓放入GPU中计算.
在基本保证成像质量的前提下,凭借GPU的计算优势,大大减少了计算周期,并且利用随机速度边界的思想解决了大规模的存储I/O问题.
不过在模型很大的情况下,还有可能遇到GPU显存不足的问题,还有待继续研究.
致谢感谢北京吉星吉达科技有限公司允许发表本文;感谢BP公司提供模型数据.
参考文献(References)[1]HemonC.
Equationsd′ondeetmodeles.
犌犲狅狆犺狔狊犻犮犪犾犘狉狅狊狆犲犮狋犻狀犵,1978,26:790~821[2]WhitmoreDW.
Iterativedepthmigrationbybackwardtimepropagation.
53rdAnnualInternationalMeeting,SEG,ExpandedAbstracts,1983.
382~385[3]BaysalE.
KosloffDD,SherwoodJWC.
Reversetimemigration.
犌犲狅狆犺狔狊犻犮狊,1983,48:1514~1524[4]McMechanGA.
Migrationbyextrapolationoftimedependentboundaryvalues.
犌犲狅狆犺狔狊犻犮犪犾犘狉狅狊狆犲犮狋犻狀犵,1983,31:413~420[5]ClaerboutJ.
Towardaunifiedtheoryofreflectormapping.
犌犲狅狆犺狔狊犻犮狊,1971,36:467~481[6]YoonK,MarfurtKJ,StarrW.
Challengesinreversetimemigration:74thAnnualInternationalMeeting,SEG,ExpandedAbstracts,2004.
1057~1060[7]YuZhang,JamesSun,SamuelGray.
Reversetimemigration:Amplitudeandimplementationissues.
78thAnnualInternationalMeeting,SEG,ExpandedAbstracts,2009.
2145~2149[8]Houzhu(James)Zhang,YuZhang.
Reversetimemigrationin3DheterogeneousTTImedia.
78thAnnualInternationalMeeting,SEG,ExpandedAbstracts,2009.
2196~2200[9]RobertW.
ClaytonandEngquist,Absorbingboundaryconditionsforwaveequationmigration.
犌犲狅狆犺狔狊犻犮狊,1980,45:95~101[10]VirieuxJ.
Psvwavepropagationinheterogeneousmedia:Velocitystressfinitedifferencemethod.
犌犲狅狆犺狔狊犻犮狊,1986,51:889~901[11]MarfurtKJ.
Accuracyoffinitedifferenceandfiniteelementmodelingofthescalarandelasticwaveequations.
犌犲狅狆犺狔狊犻犮狊,1984,49:533~549[12]SymesWW.
MathematicalFoundationsofReflectionSeismology.
TechnicalReport,RiceUniversity.
1995[13]SymesWW.
Reversetimemigrationwithoptimalcheckpointing.
犌犲狅狆犺狔狊犻犮狊,2007,72:213~221[14]GriewankA.
Achievinglogarithmicgrowthoftemporalandspatialcomplexityinreverseautomaticdifferentiation.
犗狆狋犻犿犻狕犪狋犻狅狀犕犲狋犺狅犱狊犪狀犱犛狅犳狋狑犪狉犲,1992,1:35~54[15]RobertG.
Clapp,Reversetimemigrationwithrandomboundaries.
79thAnnualInternationalMeeting,SEG,ExpandedAbstracts,2009.
2809~2813(本文编辑汪海英)3492