专利名称：一种海量数据gpu波动方程逆时偏移成像方法
技术领域：
本发明涉及一种海量数据GPU波动方程逆时偏移成像方法，应用于石油勘探、天然气勘探、煤田勘探、矿产勘探等方面，属于地震勘探与开发领域。
背景技术：
随着油气田勘探程度的不断深化，盐下油气藏勘探日益受到关注。采用传统的偏移方法无法圈定盐丘边界并对盐下构造进行成像，逆时偏移(RTM)采用双程波动方程，可以精确地描述波的传播过程，利用多次波、回转波等对陡倾角及反转构造进行准确成像，应用目前有限的偏移方法进行复杂地质构造成像。逆时偏移技术作为复杂条件地震成像难题的解决方案，可以在遇到复杂险峻倾斜的反射面例如盐丘侧翼构造时，克服目前偏移方法的局限性，能够使大于70度的倾斜构造准确成像。逆时偏移技术采用双程声波波动方程，能够使用偏移校正包括垂直倾斜在内的所有地质倾斜，还可以实现倾斜构造的成像并保持振幅不变。逆时偏移提供的是一种可供选择的偏移方法可在数据采集之后运行逆时偏移，选用适当的波场通过包括盐下、构造倾角大于70度等复杂的速度区域，以及在边界可能产生复杂的多次反射区域进行准确成像。逆时偏移方法是通过双程波波动方程在时间域上对人工给予的震源子波正向传播和接受到的地震资料进行反向传播，结合成像条件实现偏移(ClaerbOUt，1971)。由于逆时偏移需在同时计算震源正向传播和地震资料的反向传播过程中各个时刻的波�。蚱涓髯允奔涞难油胤较虿煌率故导视τ弥斜匦氪娲⑵渲幸桓龇较虻拇ス蹋庑枰醮蠖钔獯娲⒖占�。并且如今的地震资料采集已经发展到了三维甚至四维，其数据量更是十分庞大，已经达到TB级甚至IOTB级的叠前单炮数据。如此庞大的数据的IO时间将成为逆时偏移在推向海量数据的应用时不可回避的困难。例如，计算周期过长、容错能力弱、不易反复试验等问题。另一方面，相较于单程波方程的波场延拓而言，逆时偏移运用的是双程波波动方程进行波场延拓，因其避免了上下行波的分离处理，因而成为最准确的成像算法，且不受倾角的限制，并能实现回转波和多次波成像。正基于此，地球物理学家从来没有放弃过对推进逆时偏移的追求。GPU英文全称Graphic Processing Unit，中文翻译为“图形处理器”，是从CPU(中央处理器)衍生出来的的一个概念，它是显卡的心脏，和显卡板载内存成为一个子系统，共同决定PC系统的图形处理性能。图形处理器(GPU)用于通用计算(GPGPU)及其相关方面的问题目前已成为一个热门话题。GPU专为图像处理设计，存储系统实际上是一个二维的分段存储空间，包括一个区段号(从中读取图像)和二维地址(图像中的X、Y坐标)。此外， 没有任何间接写指令。输出写地址由光栅处理器确定，而且不能由程序改变。这对于自然分布在存储器之中的算法而言是极大的挑战，存储器的访问效率是算法效率的关键问题所在。不同碎片的处理过程间不允许通信是一个SIMD数据并行执行单元，在所有碎片中独立执行代码。类似CUDA和OpenCL这样的高级编程语言的问世使编程新手也能够掌握GPU的
3计算性能优势。利用GPU实现三维逆时偏移目前比较流行的方法是采用随机边界条件以计算换存储的方法来实现，该方法是利用了波场计算的可逆性的特点将边界反射随机化。但是这种方法会引入边界随机反射噪音在成像质量上难以上升。另外，若采用吸收边界条件则会由于存储的不连续性而大大降低GPU的实际计算效率。
发明内容
针对现有技术的缺陷，本发明公开了一种海量数据GPU波动方程逆时偏移成像方法，不仅克服了传统的逆时偏移方法在处理海量数据的应用时不可回避的困难，如计算周期过长、容错能力弱、不易反复试验等问题，并且也克服了近些年新兴起的利用GPU采用随机边界条件以计算换存储的方法来实现三维逆时偏移所导致成像质量差、计算效率低的缺陷。为实现上述目的，本发明是通过如下技术方案实现的一种海量数据GPU波动方程逆时偏移成像方法，逆时偏移的整体步骤包括三个�？椋莘峙淠？�、逆时偏移�？椤⑹凳笔莺喜⒛？椋鍪莘峙淠？榻ヅ谑莘峙涓煌腉PU节点，所述逆时偏移�？榻ヅ谑菁扑愠上癫⑹涑龀上裎募鍪凳笔莺喜⒛？槎猿上裎募凳鄙瑁⑼臣品掷啵扑阌没枰男畔⑽募鍪痛笮。园俜直鹊男问酵ㄖ没枰嗌偈奔渫瓿墒莺喜ⅲü朔椒ǎ扑忝芗牟ǔ〈�、边界条件和成像部分实现了借助GPU并行计算。其中，所述数据分配�？槔寐只患扑憬诘愕姆绞桨训ヅ谑莘峙涓煌腉PU 节点，每个节点领取一个单炮的计算任务。例如假设有12个计算节点，48个计算GPU，总共需要偏移的总炮数为480炮，则每个节点平均分配40炮。进一步的，所述数据分配模块设有任务堆栈，所述任务堆栈存取分配好的任务信息，所述数据分配模块继续分配任务，所述逆时偏移�？榇尤挝穸颜蝗チ烊∑迫挝瘢庋褪迪至瞬恍枰却笮？榈耐瓿删涂梢约绦峙淙挝瘢梢允迪质萑挝穹峙淠？槎懒⒒欣诔浞质迪指涸鼐猓凰鋈挝穸颜槐４婷扛龅ヅ谑荼环峙涞降募扑憬诘�、节点内的GPU卡编号以及相同GPU卡的单炮顺序号。在上述基础上，任务堆栈可以长期保存，只要输入数据和使用的计算资源不发生变化，任务的分配方式可以多次使用。如本领域技术人员所知道的，逆时偏移部分包括成像区域的划分、炮点波场延拓、 吸收边界的独立存储和计算、炮点和检波点的反向传播、相关成像等处理过程，在本发明中，所述逆时偏移�？榈某上袂蚧质辈捎孟然稚疃确较蛉缓蠡炙椒较虻姆绞健Ｓ叛〉模耸视Φ夭闼俣鹊谋浠媛晌阕幢浠痉⒚鹘疃确较蜃魑蚧值氖籽》较�。所述成像区域深度划分完成满足每一个成像层的的存储空间小于GPU的显存容量，然后在每个成像层中做水平方向的划分；所述成像区域水平方向按照方形网格划分，与GPU内核数目相适应，在每个方向上保留差分阶数所需要的边界宽度。针对已有的GPU的运算单元特点，所述GPU内核以16个线程为一组，所述水平方向按照16x16的方形网格划分。在本发明中，所用的划分方式是为了适应GPU内核线程组的特点而设计的
用公式表示为W= 16XN+2n，其中，W是成像子区域的宽度，N为任意大于0的整数，η是差分阶数。在满足GPU的运算特点的同时，上述的划分方法也满足空间高阶差分法的局部化的计算要求，例如三维8阶空间差分需要读取周围25个网格点的数据。由于 GPU全局内存读取的延迟很高，在实际运行时可以先把需要的数据并行装载到共享存储器中，然后利用GPU对共享存储器的读取速度快的特点实现所述划分方法，大幅度提高计算效率。如本领域技术人员所了解的，逆时偏移的波场延拓过程会涉及到网格边界的影响，已有的GPU逆时偏移方法成像网格包括边界一起计算，并且存储在一个连续内存空间内。在本发明中，所述逆时偏移�？榈牟ǔ⊙油毓躺婕暗耐癖呓绲耐翊娲⒎绞讲捎帽呓缤窈统上裢穹掷耄讯叻直鸫娲⒃诹娇槎懒⒍伊哪诖婵占渲校员呓绮糠值ザ来怼Ｍü鲜龇椒ǎ饩隽� GPU处理非连续边界存储时的内存不连续所导致的计算效率降低的问题。尽管边界条件的中间结果需要借助通讯保存在CPU内存或硬盘上，但存储量已经大幅度的降低，易于实现；所用的方法即解决海量数据IO问题，又能保证成像质量， 避免随机噪音的产生。其中，在逆时偏移过程中以日期、时间、CPU节点号、GPU节点号、单炮的炮号、运行的时长、开始运行的时间来命名输出的单个成像文件，这样有利于数据合并时的实时统计功能的实现，并且适合于作业的断点恢复功能。例如命名为tmp20110910-1020-CPU0-GPU0-shotl23-2h-0910-0820。在本发明中，实时数据合并�？椴捎弥付柯嫉氖凳鄙璺绞剑栉募袢∥募畔ⅲü阅柯寄诘奈募鐾臣品掷啵扑阌没枰男畔⑽募鍪痛笮。园俜直鹊男问酵ㄖ没枰嗌偈奔渫瓿墒莺喜ⅰＭü庵址椒ǎ没Ш芊奖愕氖褂锰逖�。例如，可以随时统计某一时间段内的偏移结果，可以按照计算节点统计偏移结果， 可以利用炮号统计偏移结果等等。无需等待所有偏移完成后再做数据合并，这样可以方便用户在偏移过程中做质量监控。另外，所述的数据合并方案可以方便用户在偏移完成后随意的从偏移结果中抽取所需的成像道集或剖面，无需再次偏移，具体实现方式如下文件目录统计的目的是通过用户给定的通配文件名称，例如“tmp20110910*”的方式统计当前目录中的文件能够匹配该通配文件名的个数，将所有的单炮偏移结果做全部扫描，给出扫描结果。包括需要合并的单炮数据总数，在实际合并时用当前读取的文件数除以单炮总数作为数据合并的进度百分比显示给用户。这样用文件通配名称的方式和偏移输出文件名的命名方式可以实现用户不同需求的文件合并作业提交。例如统计2011年9月10 日完成的偏移单炮数据的成像结果。在本发明中，实时数据合并�？橹杏没Э梢宰杂啥ㄒ迨涑龀上窨占涞姆段В缰鞑庀�(inline)范围、联络测线(crossline)范围、偏移距(offset)范围。通过成像空间的并行分配给不同的节点可以实现数据合并的并行化处理。通过本发明所述的方法，用于波动方程逆时偏移成像，不仅可以实现对海量数据的逆时偏移的计算处理，并且计算效率高，计算速度快，容错能力强，并有效克服了海量成像时GPU内存不足的问题，成像质量好。
图1为本发明的用GPU实现逆时偏移流程图；图2为本发明的方法数据分配�？槔萌挝穸颜皇迪指涸鼐獾呐谑莘峙涫疽馔迹渲胁煌丈聿煌诘惴峙涞降牡ヅ冢煌�3为本发明的方法逆时偏移�？椴罘殖上袂蚧滞迹渲校珿PU16个线程为一组；图4为本发明的方法逆时偏移�？楸呓缍懒⒋娲⑹疽馔迹煌�5为本发明的方法观测系统建立示意图；图6为常规逆时偏移算法实现流程图；图7为本专利的方法逆时偏移算法实现流程图。
具体实施例方式结合附图和具体实施可以更好的理解本发明的发明精神，本领域技术人员可以理解，所提供的具体实施仅为优选，而非仅限于采用此实现。本领域技术人员在理解本发明技术方案的基础上进行的变更，例如采用不同的GPU、CPU内核线程组数目改变等，依旧属于本发明的保护范围，本发明的保护涵盖于其权利要求及权利要求的等同变换。在以下实施例中，采用的GPU型号为NVIDIA公司生产的I1esIa C1060，并且下述方法适用于Tesla C1060之后发布的所有支持CUDA的GPU。本领域技术人员可以预期采用其他的通用计算方案平台在理解本发明原理的基础上也可以实现本发明所述的功能，例如MS 的DirectCompute通用计算平台等，这依旧属于本发明的保护范围。参考图1，公开了本发明的方法的流程框图，其中，“地震数据输入分配任务节点， 并写入堆栈”、“任务数据堆栈”构成了数据分配�？榈墓ぷ髂谌荩弧澳媸逼啤�、“单炮成像数据输出”构成了逆时偏移模块的工作内容；“合并数据抽取偏移距(offset)共像点道集”构成了实时数据合并�？榈墓ぷ髂谌�。参考图2，显示了本发明的方法数据分配�？榈娜挝穸颜皇迪侄缘ヅ诩扑闳挝穹峙涞慕峁迪至瞬煌扑憬诘愕呐谑峙涓涸鼐�。参考图3，显示了本发明的方法逆时偏移�？槌上袂蚧值姆椒ǎ质辈捎昧讼然稚疃确较蛉缓蠡炙椒较虻姆绞剑渲兴玫腉PU内核是以16个线程为一组；本领域技术人员可以理解对于其他不同类型的GPU，对划分方式相应调整即可。参考图4，显示了本发明的方法逆时偏移�？槎杂谖毡呓绲耐翊娲⒎绞剑捎帽呓缤窈统上裢穹掷氲姆桨福讯叻直鸫娲⒃诹娇槎懒⒍伊哪诖婵占渲�。参考图6和图7，可以充分理解本发明的方法相对于已有的方法的不同以及所进行的改进和改进的技术效果。为了更好的说明本发明方法的实现方式，申请人在此提供一个计算示例运行环境=Linux操作系统输入数据地震数据的炮记录(按炮点分选，su格式)、地质速度模型(二进制格式)输出数据地震偏移成像剖面数据体(su格式)
其中，观测系统定义利用原点坐标和方位角信息定义矩形处理工区网格，原点用 fx = 0和fy = 0来指定，方位角用参数sita = 30来指定。其中fx，fy分别代表速度模型起始点的地面坐标，sita表示inline方向与X轴的夹角。如附图5所示。数据规则化和任务分配将地震数据的炮记录作为输入，并给出所用的节点个数和观测系统中的原点坐标和方位角参数，运行rtmspd。然后，用desdata. sh将准备好的数据分发到各个节点去。任务的分配，首先手工建立一个标准的参数文件，该参数文件的实例如下
权利要求
1.一种海量数据GPU波动方程逆时偏移成像方法，其特征在于逆时偏移的整体步骤包括三个�？椋莘峙淠？�、逆时偏移�？�、实时数据合并�？椋鍪莘峙淠？榻ヅ谑莘峙涓煌腉PU节点，所述逆时偏移�？榻ヅ谑菁扑愠上癫⑹涑龀上裎募鍪凳笔莺喜⒛？槎猿上裎募凳鄙瑁⑼臣品掷啵扑阌没枰男畔⑽募鍪痛笮。园俜直鹊男问酵ㄖ没枰嗌偈奔渫瓿墒莺喜�。
2.根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述数据分配�？槔寐只患扑憬诘愕姆绞桨训ヅ谑莘峙涓煌腉PU节点，每个节点领取一个单炮的计算任务。
3.根据权利要求2所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述数据分配�？樯栌腥挝穸颜唬鋈挝穸颜淮嫒》峙浜玫娜挝裥畔ⅲ鍪莘峙淠？榧绦峙淙挝瘢鋈挝穸颜槐４婷扛龅ヅ谑荼环峙涞降募扑憬诘�、节点内的GPU卡编号以及相同GPU卡的单炮顺序号。
4.根据权利要求3所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述逆时偏移�？榇尤挝穸颜蝗チ烊∑迫挝�。
5.根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述逆时偏移�？榈某上袂蚧质辈捎孟然稚疃确较蛉缓蠡炙椒较虻姆绞�。
6.根据权利要求5所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述成像区域深度划分完成满足每一个成像层的的存储空间小于GPU的显存容量，然后在每个成像层中做水平方向的划分；所述成像区域水平方向按照方形网格划分，与GPU内核数目相适应，在每个方向上保留差分阶数所需要的边界宽度。
7.根据权利要求6所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述水平方向按照16x16的方形网格划分，所述GPU内核以16个线程为一组。
8.根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述逆时偏移�？榈牟ǔ⊙油毓躺婕暗耐癖呓绲耐翊娲⒎绞讲捎帽呓缤窈统上裢穹掷耄讯叻直鸫娲⒃诹娇槎懒⒍伊哪诖婵占渲校员呓绮糠值ザ来�。
9.根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述逆时偏移�？槟媸逼乒讨幸匀掌�、时间、CPU节点号、GPU节点号、单炮的炮号、运行的时长、开始运行的时间来命名输出的单个成像文件。
10.根据权利要求1所述的海量数据GPU波动方程逆时偏移成像方法，其特征在于所述实时数据合并模块采用指定目录的实时扫描方式，通过文件名对目录内的文件做统计分类，计算用户所需要的信息文件个数和大�。园俜直鹊男问酵ㄖ没枰嗌偈奔渫瓿墒莺喜�。
全文摘要
本发明公开了一种海量数据GPU波动方程逆时偏移成像方法，逆时偏移的整体步骤包括三个模块，数据分配�？�、逆时偏移�？�、实时数据合并�？椋鍪莘峙淠？榻ヅ谑莘峙涓煌腉PU节点，所述逆时偏移模块将单炮数据计算成像并输出成像文件，所述实时数据合并�？槎猿上裎募凳鄙瑁⑼臣品掷啵扑阌没枰男畔⑽募鍪痛笮。园俜直鹊男问酵ㄖ没枰嗌偈奔渫瓿墒莺喜�。本发明的海量数据GPU波动方程逆时偏移成像方法，三个�？橹淇梢远懒⒃诵胁⑿惺凳屎嫌诤Ａ渴莸哪媸逼剖凳扑阈矢撸扑闼俣瓤欤荽砟芰η浚⒂行Э朔撕Ａ砍上袷盙PU内存不足的问题。
文档编号G01V1/28GK102565854SQ20111044369
公开日2012年7月11日 申请日期2011年12月27日 优先权日2011年12月27日
发明者李博 申请人:中国科学院地质与地球物理研究所
专业数控铣床产品供应商
20年生产经验，实力工厂，精工品质，质量有保障
