本发明属于地震数据数值模拟和偏移成像领域，具体涉及一种基于GPU的波动方程反偏移方法，可用于地震数据正演模拟和最小二乘偏移数据重构。
背景技术：
：：地震数据数值模拟是研究地震波在地下介质中传播特征的重要手段，也是地震成像(包括逆时偏移，最小二乘偏移、全波形反演等)实现过程中的一个重要步骤。基于偏移结果的反偏移是一种地震数据数值模拟方法，它是偏移成像的反过程，由偏移结果通过数学运算得到叠前道集数据。目前常用的反偏移是基于kirchhoff积分算子的射线类反偏移方法(可参考：孙建国，Kirchhoff型真振幅偏移与反偏移[J]，勘探地球物理进展，2002，25(6)：1-5；查树贵，陈洪堤，龚小金，反偏移技术对叠前资料的噪声消除[J]，中外能源，2006，11(3)：16-20；孙建国，Kirchhoff型高保真反偏移理论[J]，地球科学与环境学报，2011，33(2)：207-212，216；杨庆道，蒋兵，岳玉波，郭朝斌，基于Kirchhoff偏移/反偏移的随机噪声压制方法[J]，石油物探，2013，52(5)：530-536)，可以较好地模拟一次波，但对复杂波场无法精确模拟。基于双程波方程的逆时反偏移方法(可参考：YuZhang,LianDuanandYiXie.Astableandpracticalimplementationofleast-squaresreversetimemigration[J].ExpandedAbstractsof83thAnnualInternationalMeeting,SEG,2013,3715-3720；YuZhang,LianDuan.Predictingmultiplesusingareversetimedemigration[J].ExpandedAbstractsof82ndAnnualInternationalMeeting,SEG,2012,1-5)可以预测一次波、多次波，是一种精确的反偏移技术，然而这种方法计算量大，目前基于CPU(CentralProcessingUnit，中央处理器)的实现算法效率较低。近几年，出现一种基于扩展成像结果的波 动方程反偏移方法(可参考：HervéChauris1andMondherBenjemaa.Seismicwave-equationdemigration/migration[J].Geophysics,2010，75(3):S111-S119；WiktorWaldemarWeibullandArntsen.Reverse-timedemigrationusingtheextended-imagingcondition[J].Geophysics,2014,79(3):WA97-WA105)，依据爆炸反射面的思想，利用波动方程算子对每一个偏移距的成像剖面进行反偏移。这种方法效率较快，但无法得到叠前炮集数据体，因而算不上严格意义的反偏移。上述现有反偏移方法技术在效果或效率方面无法满足实际需求，较难得到应用和推广。技术实现要素：本发明针对现有反偏移方法存在的效果不佳或效率较低的问题，提供一种基于GPU(GraphicProcessingUnit，图形处理器)的波动方程反偏移方法，输入光滑速度模型和已有的偏移成像，逐炮按照炮点坐标和接收点坐标，利用时间空间域波动方程高阶差分算子在GPU计算平台对背景波场和散射波场进行波场传播外推计算，每个时刻在地表记录散射波场，得到反偏移的叠前炮集数据体。本发明提供一种基于GPU的波动方程反偏移方法，其特征在于，包括以下步骤：输入速度模型、偏移结果和控制参数；CPU获得炮点位置和震源子波；CPU获得接收点位置，确定炮数据范围及波动方程计算网格，并从输入的速度模型和偏移结果中获取计算网格内的速度c0和偏移值m；将CPU上的炮点位置、震源子波、接收点位置、计算网格内的速度c0和偏移值m传送到GPU；在GPU上进行散射波场延拓并记录反偏移波场；将GPU上记录的反偏移波场传送到CPU；以及CPU将记录的反偏移波场写入存储器。进一步地，在GPU上进行散射波场延拓并记录反偏移波场的步骤包括：在炮点处加载震源子波；利用计算网格内的速度c0进行背景波场u0的传播；利用计算网格内的速度c0和偏移值m进行散射波场u1的传播；在接收点处记录该时刻的 散射波场u1的值，即为反偏移波场。进一步地，利用计算网格内的速度c0进行背景波场u0的传播，是基于输入的速度模型和已知的震源位置，采用时间空间域双程波方程高阶有限差分算子对地震波场进行正向延拓。进一步地，利用计算网格内的速度c0和偏移值m进行散射波场u1的传播，是基于输入的速度模型和偏移结果，采用时间空间域双程波方程高阶有限差分算子对地震波场进行正向延拓。进一步地，时间空间域双程波方程为：其中，t表示时间，c0表示计算网格内的速度，m表示偏移值，u0表示背景波场，u1表示散射波场。优选地，所述利用计算网格内的速度c0进行背景波场u0的传播以及所述利用计算网格内的速度c0和偏移值m进行散射波场u1的传播，是通过利用GPU并行技术进行计算的。在时间空间域双程波方程高阶有限差分计算中，x、y，z三个方向上的网格数为Nx、Ny、Nz，则共需要对Nx×Ny×Nz个网格点进行差分计算，每个网格点的计算是相互独立的，由GPU独立地同时进行并行计算。GPU并行计算中采用每个线程进行一个网格点的差分计算，采用多线程的并行计算技术。进一步地，本发明的基于GPU的波动方程反偏移方法可以逐个炮点计算反偏移。进一步地，本发明的基于GPU的波动方程反偏移方法可以逐个时刻进行散射波场延拓并记录反偏移波场。本发明针对现有反偏移技术效果不佳或效率较低的问题，采用双程波方程进行散射波场传播获得反偏移数据，可由偏移结果模拟一次波和多次波，是一种精确的反偏移方法。此外，为应对双程波方程高阶差分计算效率低的缺点，采取基于GPU多线程的并行计算技术，大幅提高波动方程反偏移的计算速度。 本发明形成了高精度、高效率的反偏移方法，可用于地震数据正演模拟和最小二乘偏移数据重构。本发明的其它特征和优点将在随后的具体实施方式部分予以详细说明。附图说明附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：图1是根据本发明一种实施方式的基于GPU的波动方程反偏移方法的步骤流程图；图2是根据本发明一种实施方式的基于GPU进行散射波场延拓并记录反偏移波场的步骤流程图；图3是在一个2D理论模型数据的反偏移计算中输入的偏移结果；图4是根据一个实施例的反偏移得到的一炮地震数据；图5是根据一个实施例的原始炮集数据；图6是根据一个实施例的伊拉克某工区一条偏移剖面；以及图7是根据一个实施例的伊拉克某工区反偏移得到的一炮数据。具体实施方式以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。根据本发明的一种实施方式，提供一种基于GPU的波动方程反偏移方法，其特征在于，包括以下步骤：输入速度模型、偏移结果和控制参数；CPU获得炮点位置和震源子波；CPU获得接收点位置，确定炮数据范围及波动方程计算网格，并从输入的速度模型和偏移结果中获取计算网格内的速度c0和偏移值m；将CPU上的炮点位置、震源子波、接收点位置、计算网格内的速度c0和偏移值m 传送到GPU；在GPU上进行散射波场延拓并记录反偏移波场；将GPU上记录的反偏移波场传送到CPU；以及CPU将记录的反偏移波场写入存储器。本实施方式的基于GPU的波动方程反偏移方法可以逐个炮点计算反偏移，具体地，可以在输入速度模型、偏移结果和控制参数之前或者之后，判断是否针对所有炮点计算完毕，如果计算完毕则方法流程结束，否则开始进行下一个炮点的计算，计算完毕后CPU将记录的反偏移波场写入存储器，再次进行判断，直至针对所有的炮点计算完毕。进一步地，在GPU上进行散射波场延拓并记录反偏移波场的步骤包括：在炮点处加载震源子波；利用计算网格内的速度c0进行背景波场u0的传播；利用计算网格内的速度c0和偏移值m进行散射波场u1的传播；在接收点处记录该时刻的散射波场u1的值，即为反偏移波场。利用计算网格内的速度c0进行背景波场u0的传播，是基于输入的速度模型和已知的震源位置，采用时间空间域双程波方程高阶有限差分算子对地震波场进行正向延拓。利用计算网格内的速度c0和偏移值m进行散射波场u1的传播，是基于输入的速度模型和偏移结果，采用时间空间域双程波方程高阶有限差分算子对地震波场进行正向延拓。本实施方式的基于GPU的波动方程反偏移方法可以逐个时刻进行散射波场延拓并记录反偏移波场。具体地，在开始进行散射波场延拓时，首先判断时间循环是否结束，如果时间循环结束，则散射波场延拓流程结束；如果时间循环没有结束，则开始进行下一时刻的散射波场延拓。针对一个时刻的散射波场延拓结束后，再次进行时间循环是否结束的判断，直至针对所有时刻的散射波场延拓结束。优选地，时间空间域双程波方程为：其中，t表示时间，c0表示计算网格内的速度，m表示偏移值，u0表示背景波场，u1表示散射波场。优选地，所述利用计算网格内的速度c0进行背景波场u0的传播以及所述利用计算网格内的速度c0和偏移值m进行散射波场u1的传播，是通过利用GPU并行 技术进行计算的。在时间空间域双程波方程高阶有限差分计算中，x、y，z三个方向上的网格数为Nx、Ny、Nz，则共需要对Nx×Ny×Nz个网格点进行差分计算，每个网格点的计算是相互独立的，由GPU独立地同时进行并行计算。GPU并行计算中采用每个线程进行一个网格点的差分计算，采用多线程的并行计算技术。下面参照附图描述本发明的具体实施方式。图1是根据本发明一种实施方式的基于GPU的波动方程反偏移方法的步骤流程图。首先输入速度模型、偏移结果及关键控制参数。反偏移是偏移的逆过程，偏移时输入实际地震采集数据和由该数据分析得到的速度模型，通过运算输出偏移结果。反偏移是把速度模型、偏移结果，以及实际地震采集数据(利用炮点位置和接收点位置的信息)作为输入，输出模拟的地震采集数据。控制参数包括反偏移的炮数、反偏移子波频率、反偏移计算步数、反偏移计算步长、速度模型x、y、z三个方向的网格数、速度模型x、y、z三个方向的网格间隔等。然后判断是否针对所有炮点计算完毕，如果计算完毕则方法流程结束，否则开始进行下一个炮点的计算。本实施例中逐炮进行反偏移计算并保存反偏移结果，所有炮计算完毕则反偏移结束。接下来的流程是针对一个炮点的计算。1)CPU上获得炮点位置和震源子波。炮点位置从实际地震采集数据中获得，震源子波可以是人为给定的。2)CPU上获得接收点位置，确定炮数据范围及波动方程计算网格，并从输入的速度模型和偏移结果中获取计算网格内的速度c0和偏移值m。接收点位置可从实际地震采集数据中获得。炮数据范围及波动方程计算网格是根据接收点位置确定的。通常输入的偏移结果是整体的，一炮范围对应的计算网格是其中的一部分，由炮范围的坐标值可以从整体偏移结果中抽取该炮范围内的偏移结果。而输入的速度模型也是整体的，一个炮范围对应的计算网格是其中的一部分，因此由炮范围的坐标值可以从整体速度模型中抽取该炮范围内的速度模型。3)将CPU上的炮点位置、震源子波、接收点位置和计算网格内的速度c0和 偏移值m等数据传送到GPU。4)在GPU上进行散射波场延拓并记录反偏移波场。5)将GPU上记录的反偏移波场发送到CPU。6)CPU将反偏移波场写入存储器。优选地，在GPU上进行散射波场延拓并记录反偏移波场的步骤4)，可参见附图2。图2是根据本发明一种实施方式的基于GPU进行散射波场延拓并记录反偏移波场的步骤流程图。本实施方式中，逐个时刻进行散射波场延拓和记录偏移波场，每个时刻的实现方法是：41)在炮点处加载震源子波；42)利用速度c0进行背景波场u0的传播；43)利用速度c0和偏移值m进行散射波场u1的传播；44)在接收点处记录该时刻的散射波场u1的值，即为反偏移波场。在利用速度c0进行背景波场u0的传播的步骤42)中，背景波场u0的正向传播，是基于输入的光滑速度模型和已知的震源位置，采用时间空间域双程波方程高阶有限差分算子对地震波场进行正向延拓。可选地，震源传播双程波方程和高阶差分算子可参考“刘守伟、王华忠、陈生昌等，三维逆时偏移GPU/CPU机群实现方案研究[J]，地球物理学报，2013，56(10)：3487-3496”。其中震源位置可从输入的地震采集数据中获得。在利用速度c0和偏移值m进行散射波场u1的传播的步骤43)中，散射波场u1的正向传播，是基于输入的光滑速度模型和偏移结果，采用时间空间域双程波方程高阶有限差分算子对地震波场进行正向延拓。上面所述的时间空间域双程波方程为：1c02∂2u1∂2t-▿2u1=-m∂2u0∂2t]]>其中，t表示时间，c0表示速度，m表示偏移值，相当于速度扰动，u0表示传播的背景波场，u1表示散射波场。在接收点位置记录每一个时刻的偏移波场， 得到的地震记录即为一炮的波动方程反偏移结果。优选地，上述步骤42)、43)是利用GPU并行技术进行计算。在时间空间域双程波方程高阶有限差分计算中，假设x、y，z三个方向上的网格数为Nx、Ny、Nz，则共需要对Nx×Ny×Nz个点进行差分计算，而每个点的计算是相互独立的，由GPU独立地同时进行并行计算。GPU计算中，内核函数是以线程网格(Grid)的形式组织的，每个线程网格由若干个线程块(Block)组成，而每个线程块又由多个线程(thread)组成。因此，采用每个线程进行一个网格点的差分计算，可以实现并行快速运算，提高了计算效率。下面参照图3-图7描述本发明的两个具体实施例。实施例一：一个2D理论模型数据的反偏移计算图3是输入的偏移结果，图4是通过反偏移得到的某炮集记录，图5是切除了直达波后的原始炮集记录。对比图4和图5，两者比较接近，由此说明根据本发明方法的反偏移可以很好地模拟一次波和多次波。实施例二：伊拉克某3D工区地震数据反偏移计算图6是输入的偏移结果中的一条剖面，由此利用本发明方法进行反偏移计算，图7是反偏移得到的一炮地震数据。反偏移计算网格为356×191×801，波场传播时间步数为5000步，时间间隔为1ms，采用NVIDIA公司产品TeslaM2090型号的GPU卡进行计算，一炮反偏移用时23分钟。对于相同的计算量，CPU计算用时105分钟。可见，本发明通过GPU多线程并行运算，相比于现有方法中通过CPU运算，具有很高的计算效率。本发明针对现有反偏移技术效果不佳或效率较低的问题，提出一种基于GPU的波动方程反偏移方法。通过逐炮进行反偏移计算，对每一炮而言，CPU先将炮点、接收点信息和震源子波、速度、偏移值发送到GPU，然后在GPU上逐个时刻进行散射波场延拓并进行采样记录，记录的散射波场即为一炮的反偏移结果，再将其发送到CPU上并写入存储器。通过本发明的方法反偏移精度好，计算效率高，可用于地震数据正演模拟和最小二乘偏移数据重构。以上结合附图详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。当前第1页1 2 3 当前第1页1 2 3