基于稳相区间理论的逆时偏移成像方法与流程

本发明属于地球物理
技术领域：
，具体涉及一种逆时偏移成像方法，特别涉及一种基于稳相区间理论的逆时偏移成像方法。
背景技术：
：逆时偏移技术克服了单成波成像传播倾角限制和kirchhoff成像中的高频近似的假设，在振幅恢复上优于时间偏移方法，可以更好地对陡倾成像，是现有工业上的标准成像方法。逆时偏移成像主要由三部分组成，包括震源波场的正向延拓、检波点波场的反向延拓以及成像条件的应用，这样可以获得地下每一点每一时刻所有的地震记录。成像条件是逆时偏移的核心算法，成像条件可以直接影响逆时偏移算法的计算效率以及成像质量，这两点也都是工业生产所关心的核心问题。因此，从逆时偏移算法被提出以后，许多地球物理学家都致力于研究成像条件，随着偏移技术的发展多种成像条件被提出。1971年，claerbout提出零延迟互相关成像条件，这种成像条件是逆时偏移中应用最广泛的成像条件。它可以将所有的成像信息包含到成像点中，而且算法实现简单。但是，这种算法需要地下每一点的震源波场和检波点波场在每一时刻都要被同时获得，因此需要大量的存储空间，很不利于工业生产中应用。1986年，changandmcmechan提出激发时间成像条件，这种成像条件计算速度快，并且可以大大减少逆时偏移算法对计算内存的需求，但是在复杂模型中，会产生多路径问题，此时初至波携带的成像信息的波场通常会很弱，从而导致成像模糊。2009年，nguyenandmcmechan提出源波场重建技术，解决了互相关成像条件应用时的存储问题，但是却增加了50％的计算量。2013年，nguyenandmcmechan又提出激发振幅成像条件，这种成像条件在简单模型中很稳定，克服了激发时间成像条件的缺陷，但是仍然无法解决多路径问题，因为整个成像过程，这种成像条件只是对携带最大振幅的信息成像，所以在复杂模型中的成像结果比不上归一化互相关成像条件的成像结果。与此同时，归一化成像条件要求震源波场信息和检波点波场信息完全存储在计算内存中，工业生产中，大模型，大数据偏移是无法实现的。技术实现要素：本发明的目的就在于针对上述现有成像技术的不足，提供一种基于稳相区间理论的逆时偏移成像方法，其计算效率高，计算内存需求低，并且可以克服偏移成像中的多路径的问题。本发明的目的是通过以下技术方案实现的：首先，将稳相点理论与多种成像条件理论结合，提出稳相区间理论。使用程函方程计算初至波走时，与最大振幅走时场和模型复杂程度结合，计算出适合的成像时间范围，称为稳相区间，这是本方法的核心思想，取代传统逆时偏移源波场重建过程，以此节省计算时间。在此稳相区间范围内的信号是模型中每个成像点整条观测记录中成像信息最丰富的一段。我们使用这个范围内的信号成像就实现了一种高效的逆时偏移成像方法。一种基于稳相区间理论的逆时偏移成像方法，包括以下步骤：a、基于程函方程的走时计算初至波的走时场tref(x,z)，这个走时场与稳相区间结合，就可以确定地下每一点的成像范围。程函方程中，tref为时间(单位s)，s(x,z)为二维模型中的慢度分布函数(单位为s/m)，x,z为空间坐标；b、从近偏移距到远偏移距炮点分布中，选取3炮计算最大走时场tmax(x,z)，选取对应炮的初至波走时计算，消除计算系统误差，(两种算法之间的误差)。根据模型复杂程度和走时场分布的相似程度选取每一个成像时间范围nt。n与模型复杂程度和tmax(x,z)-tref(x,z)均匀程度有关,t是子波周期。结合走时场tref(x,z)，确定地下每一点的稳相区间成像时间范围{tref(x,z)-(n-1)t，tref(x,z)+nt}；c、使用高阶有限差分算法和pml边界算法结合波动方程进行正演，计算逆时偏移震源波场p(x,z,t)，模拟过程中直接记录稳相区间内的波场ps(x,z,t)；d、逆时偏移检波点波场延拓过程中，将地震记录u(x,t)，进行插值，作为边值条件，使用步骤c的波动方程对于波场进行反时间方向传播，同样提取{tref(x,z)-(n-1)t，tref(x,z)+nt}稳相区间范围内的检波点波场pr(x,z,t)；e、c和d步骤我们获的ps(x,z,t)和pr(x,z,t)使用互相关算法成像，获得i(x,z)，这个过程就是稳相区间成像条件算法。公式中i(x,z)为偏移成像。f、对于每一炮数据重复进行步骤a-步骤e，即可获得每一炮对应的成像结果，选取适当的有效偏移孔径，将所有的炮集成像结果累加在一起即可得到最终的成像结果。与现有技术相比，本发明的有益效果在于：本发明成功地将稳区间概念与多种逆时偏移成像条件结合，提出稳相区间成像条件，这种成像条件可以在消耗较少的计算机内存条件下，高效的实现逆时偏移成像技术。克服了传统逆时偏移低效率，高内存消耗的缺陷，同时解决了激发振幅成像条件无法解决的多路径问题。此外，由于我们在成像过程只是使用稳相区内的信号，可以压制一些背景低频噪声。本发明主要包含以下方面特点：1、将稳相点概念与偏移成像条件结合，提出稳相区间概念，这是本发明的核心部分。对于模型中每一点的记录，稳相区间内的信号都是以较短的时窗包含整条记录中最有效的成像部分；2、通过将程函方程计算初至走时场，与最大振幅走时场和模型复杂度结合，确定适合的稳相区间。既可以高效完成逆时偏移成像过程，保证成像质量，又可以减少偏移过程对于内存的需求；3、与传统逆时偏移的互相关成像条件方法相比较，适当的增加了较少内存需求，节省了约1/3的计算时间；4、与激发振幅成像条件方法相比较，解决了激发振幅成像条件中由于只使用最多振幅作为成像值无法克服多路径成像问题的缺陷；5、使用稳相区间成像条件，使得归一化成像条件，可以在较小的内存消耗条件下完成，可以满足工业需求，适合在实际生产中推广。附图说明图1为基于稳相区间理论的逆时偏移成像方法流程图；图2a为三层水平层状速度模型示意图；图2b为使用程函方程计算初至波的结果示意图(初至波走时场fldr＝1)；图2c为使用程函方程计算初至波的结果示意图(初至波走时场fldr＝51)；图2d为使用程函方程计算初至波的结果示意图(初至波走时场fldr＝101)；图3a为三层水平模型对应的最大振幅走时场fldr＝1；图3b为三层水平模型对应的最大振幅走时场fldr＝51；图3c为三层水平模型对应的最大振幅走时场fldr＝101；图4a为第51炮稳相区间内的震源波场的成像波场ps(x,z,t)；图4b为图4a对应的局部放大图；图4c为模型中(300,200)成像位置的震源波场记录；图5a为第51炮稳相区间内的检波点波场的成像波场pr(x,z,t)；图5b为图5a对应的局部放大图；图5c为为模型中(300,200)成像位置的检波点波场记录；图6为使用稳相区间成像条件的三层模型的偏移成像结果；图7a为倾斜层状速度模型；图7b为倾斜模型第51炮初至波走时场；图7c为倾斜模型第51炮最大振幅走时场；图8a为倾斜模型的互相关成像条件结果；图8b为倾斜模型的稳相区间成像条件结果；图8c为倾斜模型的激发振幅成像条件结果；图9a为marmousi速度模型；图9b为光滑的marmousi速度模型；图9c为marmousi模型第51炮初至波走时场；图9d为marmousi模型第51炮最大振幅走时场；图10a为marmousi模型的互相关成像条件结果；图10b为marmousi模型的稳相区间成像条件结果；图10c为marmousi模型的激发振幅成像条件结果。具体实施方式下面结合附图和实例对本发明进一步的详细说明，以简单的三层状模型为例。基于稳相区间理论的逆时偏移成像方法，包括以下步骤：a、通过基于程函方程的走时计算程序计算初至波的走时场tref(x,z)程函方程中，tref为时间(单位s),s(x,z)为二维模型中的慢度分布函数(单位为s/m),x,z为空间坐标。图2a为水平层状模型，图2b-图2d为通过程函方程计算出的初至波走时场tref(x,z)，每一炮都有对应的初至波走时场。b、从近偏移距到远偏移距的炮点分布中，选取3炮计算最大走时场场tmax(x,z)，如图3a-图3c。选取对应炮的初至波走时计算，消除计算系统误差。根据模型复杂程度和走时场的相似程度选取每一个成像时间范围nt。n与模型复杂程度和tmax(x,z)-tref(x,z)均匀程度有关,t是子波周期。通常模型很简单，例如层状模型，一个子波周期即可。而类似于marmousi的复杂模型，波场传播过程产生多路径情况，只需要对复杂区域增加n，选3个周期即可。结合走时场tref(x,z)，确定地下每一点的稳相区间成像时间范围{tref(x,z)-(n-1)t，tref(x,z)+nt}；c、将高阶有限差分算法和pml边界算法，与交错网格结合对波动方程进行离散，用于下面逆时偏移计算波场p(x,z,t)。其中，波动方程使用一阶应力速度方程：对一阶应力速度方程进行离散，结合pml边界，我们可以展开为：其中p,px,pz分别质点应力波场，以及波场延x方向和z方向的分量，vx,vz分别为质点震动速度延x方向和z方向的分量。d(x),d(z)分别为pml吸收边界x,z方向的吸收系数，完成模拟过程。在整个模拟过程直接记录每一个成像点的稳相区间内的成像波场ps(x,z,t)。我们选取模型中深度200*5m的位置，以第51炮数据为例。如图4a所示，我们可以看出稳相区间的上下两界将正传波场记录的有效成像信息包含在内，而区间以外的记录几乎为零，图4b是局部放大更加明显。我们选取(300,200)这个成像点的记录进行详细分析(图4c)。带有坐标的黑点即为稳相区间的范围，我们可以清楚看到，这个范围将整条记录中最有效的成像部分包含在内，并且包含最大振幅点(图4c中*所指示的位置),说明了稳相区间内的信息可以充分的成像，几乎不会造成有效信息浪费。d、逆时偏移检波点波场延拓过程中，将地震记录u(x,t)，进行插值，作为边值条件，使用步骤c的波动方程对于波场进行反时间方向传播，同样提取{tref(x,z)-(n-1)t，tref(x,z)+nt}稳相区间范围内的检波点波场pr(x,z,t)；如图5a-图5c所示，稳相区间把成像信息同相轴包含在稳相区间的上下两界内，而非成像信息的另一条同相轴排除在外。e、c和d步骤我们获的ps(x,z,t)和pr(x,z,t)使用互相关算法成像，获得i(x,z)，这个过程就是稳相区间成像条件算法。公式中i(x,z)为偏移成像。f、对于每一炮数据重复进行步骤a-步骤e，即可获得每一炮对应的成像结果，选取适当的有效偏移孔径，将所有的炮集成像结果累加在一起即可得到最终的成像结果。如图6所示，三层层状速度模型的偏移成像结果。本发明基于稳相区间理论的逆时偏移成像方法是在linux系统下机群上使用c语言程序进行数据处理的。实施例1我们使用倾斜模型验证我们发明的可行性。模型参数：nx＝601,nz＝601,dx＝dz＝5m,上层速度2000m/s,下层速度2200m/s。如图7a所示。观测系统：总炮数nsrc＝101，炮间距dsrc＝6,检波器数量nrecv＝601,道间距drecv＝1。使用简单的连续观测系统，检波点分布在模型表面不动，炮点从左至右等间隔排布。使用程函方程计算初至波走时场，有限差分算法计算最大振幅走时场，我们展示第51炮数据，如图7b和图7c所示。从图中我们看出两种走时场分布都很均匀，形态相似，所以稳相区间选取范围可以很小，一个子波周期就可以。我们分别使用互相关成像条件，激发振幅成像条件和我们发明的基于稳相区间理论的成像条件，完成101炮的偏移过程，比较计算时间和偏移成像效果。表13种成像条件计算时间对比结果20炮101炮互相关成像条件1333450000064583290000激发振幅成像条件941710000046072270000我们发明的成像条件951620000046983780000注：表中时间没有单位，只是计数，使用c语言自带计时函数理论上逆时偏移主要耗时过程是波场正、反外推过程。互相关成像条件，由于使用源波场重建技术，耗时应该为普通正演过程3倍，而其他两种算法应该是正演过程的两倍耗时。从表1我们可以看出使用互相关成像条件的逆时偏移计算效率最低；使用激发振幅成像条件计算效率最高，而我们的发明的方法，时间略微增加，符合理论估计。图8a-图8c为三种成像条件的结果。从图中我们可以看出，对于简单模型，三种成像条件都可以准确成像，说明我们发明的成像条件理论上是可行的，具有较高的计算效率和成像质量。实施例2我们使用marmousi模型验证我们发明的优势。模型参数：nx＝1001,nz＝401,dx＝dz＝5m,选取的部分模型，主要为层状结构，中间位置出现低速透镜体，如图9a。用于偏移的光滑模型如图9b所示。由于透镜体的存在，透镜体以下的波场传播过程会产生多路径情况。观测系统：总炮数nsrc＝101，炮间距dsrc＝10,检波器数量nrecv＝501,道间距drecv＝2。使用简单的连续观测系统，检波点分布在模型表面不动，炮点从左至右等间隔排布。使用程函方程计算初至波走时场，有限差分算法计算最大振幅走时场，我们展示第51炮数据，如图9c和图9d所示。从图中我们看出两种走时场分布形态相似，但是模型中心由于有低速透镜体出现，导致下部的走时场差异较大，所以这个区域稳相区间选取范围要适当加，使得多路径的成像信息全部包含在内，其他区域可以选择只一个子波周期就可以。我们使用互相关成像条件，激发振幅成像条件和我们发明的基于稳相区间理论的成像条件成像分别成像，如图10a-图10c所示。互相关成像条件成像过程使用了每一点每一时刻的全部成像信息，可以自动解决多路径问题，所以透镜体周围的模型界面位置图像的同相轴连续。使用我们发明的成像条件在图像质量上比的上互相关成像条件成图的图像质量，并且大大缩减了计算时间。激发振幅成像条件成图在透镜体以下和周围的区域成图质量欠佳，出现绕射偏移画弧和不连续现象。图中红色方框标注区域。本发明将稳相点理论与多种成像条件结合，提出稳相区间理论。分析多路径产生原因，结合稳相区间理论，提出基于稳相区间理论的逆时偏移成像方法，这种成像方法计算效率高，计算内存需求低，并且可以克服偏移成像中的多路径的问题。当前第1页12