基于波动方程外推算子的VTI介质中叠前时间偏移方法

摘要

本发明公开了一种基于波动方程外推算子的VTI介质中叠前时间偏移方法，包括：(1)输入CMP道集以及时间域速度场；(2)对CMP道集进行局部倾斜叠加，将CMP道集映射到截距时间-射线参数域的平面波数据体；并按照射线参数对平面波数据体分选，形成共射线参数剖面；(3)对射线参数循环；(4)进程判断所在计算机节点是否为主控节点，由此进入主-从并行模式；(6)所有进程同步；(7)主控节点将当前射线参数对应的成像数据体输出；(8)判断射线参数循环是否完成；若否，转向步骤3；(9)输出时间偏移剖面以及成像道集。本发明可以提供保真的成像道集，为后续的AVO/AVA分析提供更好的基础数据，有利于更精确的油气储层岩性预测。

说明书

技术领域

本发明涉及一种叠前时间偏移技术，尤其涉及一种一种VTI介质中叠前时间偏移方法。

背景技术

在地质构造不甚复杂的区域，叠前时间偏移后的成像道集上的AVO(振幅随偏移距的变化)/AVA(振幅随入射角的变化)分析是油气储层预测的关键技术。叠前时间偏移技术在提高地震成像精度，为叠前反演的实现等方面发挥了重要作用。目前的VTI(具有垂直对称轴的横向各向异性)介质时间偏移基本上是使用Kirchhoff积分法VTI叠前时间偏移算法，采用的是用等效各项异性参数表示的四阶方程(Alkhalifah，T，1997，Velocity analysisusing nonhyperbolic moveout in transversely isotropic media，.Geophysics，62，1839-1854)，其本质上还是基于层状介质假设下推导的时距关系。通常所采用的商用软件系统中基本上都是用Kirchhoff积分法叠前时间偏移给出Offset成像道集，并在此基础上进行AVO分析。实际上，Kirchhoff偏移方法使用的是均方根速度场，适应速度横向变化的能力较弱，当速度场横向变化时，基于常速介质中的格林函数计算得到的振幅加权因子通常是不精确的。因此，Kirchhoff积分法给出的成像道集振幅不保真、成像道集中存在假象干扰，这样的成像道集不利于后续的AVO分析和储层预测。吸收衰减和各向异性的对成像道集和AVO分析的影响也要在成像过程中加以消除。还有，随着各油田探区油气勘探程度的不断加深，对偏移成像精度的要求越来越高，各向同性处理技术已经难以满足勘探开发的要求。

迄今为止，还未见有基于波动方程外推算子的叠前时间偏移算法。与Kirchhoff偏移方法相比而言，波动方程外推算子对振幅的处理更加自然，并且对速度横向变化更强的适应能力。因此，有必要发展一套基于波动方程外推算子的VTI介质中叠前时间偏移方法。

发明内容

针对上述现有技术，从适应缓横向变速、提高成像分辨率(吸收衰减补偿)和考虑各项异性等方面综合考虑，本发明提供一种基于波动方程外推算子的VTI介质中叠前时间偏移方法。是将偏移距平面波波动方程叠前时间偏移方法至VTI介质中，从而得到VTI介质中基于波动方程的叠前时间偏移算法，本发明具有更好的保振幅能力，可以提供保真的成像道集，为后续的AVO/AVA分析提供更好的基础数据，有利于更精确的油气储层岩性预测。本发明是一种更精确高效的波动方程叠前时间偏移方法，能补偿大地的吸收衰减和处理各向异性介质中的成像，可以提供很好的叠前时间偏移叠加剖面，输出射线参数或者角度域成像道集用于后续的剩余速度分析和AVO分析。

为了解决上述技术问题，本发明一种基于波动方程外推算子的VTI介质中叠前时间偏移方法予以实现的技术方案包括以下步骤：

步骤1：输入CMP道集以及时间域速度场；

步骤2：对CMP道集进行局部倾斜叠加，将CMP道集映射到截距时间-射线参数域的平面波数据体；并按照射线参数对平面波数据体分选，形成共射线参数剖面；

步骤3：对射线参数循环；

步骤4：进程判断所在计算机节点是否为主控节点，由此进入主-从并行模式；即：若进程在主控节点上运行，则执行下述(5-A)过程，若进程未在主控节点上运行，而是在计算节点上运行，则执行下述(5-B)过程；

(5-A)首先，读取当前射线参数对应的共射线参数剖面，然后沿时间方向作一维快速傅立叶变换，生成频率域的共射线参数平面波数据体；与此同时，设定总任务数为总频率数目与计算节点数目之中的最小值；主控节点控制着并行任务的分发与收集，并执行以下步骤：

(5-A1)对总任务数循环；

(5-A2)将一个频率切片的平面波数据体发给计算节点，发送任务的计数器自增；

(5-A3)判断任务发送是否结束，若否，转向(5-A1)；若是，转向(5-A4)；

(5-A4)对总频率数循环；

(5-A5)接收计算节点返回的任务标识符；

(5-A6)检查已经发送给计算节点的任务序号是否小于总频率数目；若是，转向(5-A7)；若否，转向(5-A8)；

(5-A7)将当前频率切片的平面波数据体发给计算节点，发送任务的计数器自增，转向(5-A4)；

(5-A8)向所有计算节点发送任务结束标识符；

(5-B)首先，判断该计算节点号是否大于总的频率数目；若是，转向步骤6；若否，转向(5-B1)；

(5-B1)进程进入死循环；

(5-B2)进程等待接收从主控节点发送来的数据；

(5-B3)进程判断主控节点是否发送结束标识符；若是，转向(5-B9)；若否，转向(5-B4)；

(5-B4)沿时间方向进行波场外推；

(5-B5)计算偏移距平面波传播的相速度求解如下隐式方程：

$>f\left(\theta \right)=2\frac{\sin \theta }{v_p(\theta ;v_{p0},ϵ,\delta )}-p_h=0---\left(1\right)$>

公式(1)中，平面波相速度v_p(θ；v_p0，ε，δ)可以表示为：

$>\frac{{V_{\mathrm{qP}}}^2\left(\theta \right)}{{v_{p0}}^2}=\frac{1}{2}+ϵ{\sin }^2\theta +\frac{1}{2}\sqrt{{(1+2ϵ{\sin }^2\theta )}^2-2(ϵ-\delta ){\sin }^{2}2\theta }---\left(2\right)$>

公式(2)中，ε、δ为VTI介质各向异性参数，θ为相角，v_p0为VTI介质中沿对称轴(Z轴)方向地震波的相速度；通过公式(1)～(2)，可以得出由射线参数p_h表征的偏移距平面波传播的相速度v_p(p_h)；

(5-B6)对VTI介质中的偏移距平面波方程采用隐式有限差分格式进行波场外推，实现偏移距平面波叠前时间偏移；所采用的偏移方程为：

$>(\frac{4}{v_p^2\left(p_h\right)}\frac{{\partial }^2}{{\partial \tau }^2}+\frac{{\partial }^2}{{\partial x}^2})\tilde{U}(x,\tau ;\omega )+4\frac{{\omega }^2}{v_p^2\left(p_h\right)}(1-\frac{v_p^2\left(p_h\right){p_h}^2}{4})\tilde{U}(x,\tau ;\omega )=0---\left(3\right)$>

公式(3)中，为频率域平面波数据，ω为角频率，v_p(p_h)为(2)中计算出的平面波传播的相速度；

采用零时间成像条件，得到共射线参数域成像道集；

$>I\left({x,\tau ;p}_h\right)=\underset{\omega }{\int }\tilde{U}(x,\tau ;\omega )·{\tilde{U}}^{*}(x,\tau ;\omega )\mathrm{d\omega }---\left(4\right)$>

(5-B7)判断波场外推是否结束；若否，转向(5-B4)；若是，转向(5-B8)；

(5-B8)向主控节点发送任务结束标识符，转向(5-B1)；

(5-B9)所有计算节点上的进程向主控节点规约当前射线参数对应的成像数据体，规约操作符号为求和；

步骤6：所有进程同步；

步骤7：主控节点将当前射线参数对应的成像数据体输出；

步骤8：判断射线参数循环是否完成；若否，转向步骤3；

步骤9：输出时间偏移剖面以及成像道集。

与现有技术相比，本发明的有益效果是：

现在的商用软件系统中的叠前时间方法技术基本上都是Kirchhoff积分方法。主要是提供一种快速的、适应不规则系统的成像方法。它不适于精细的保真成像和后续AVO/AVA分析及储层预测。其所提供的成像道集中存在假象干扰，振幅保真性差，不适于后续的AVO分析。Kirchhoff积分法叠前时间偏移很难与吸收衰减补偿相结合。其他一些基于单平方根或双平方根方程的波动方程偏移方法都具有巨大的计算量，不能大规模用于生产实际，这也是波动方程叠前时间偏移不能在实际应用中普及的原因之一。

本发明VTI介质中叠前时间偏移方法中的基于波动方程的波场外推算子可以提供更加保振幅的成像剖面；在平面波域实现偏移具有很高的计算效率；利用本发明在波场外推过程中利用本发明可以自适应的计算平面波的相速度而无需射线追踪；同时可以输出射线参数或者角度域成像道集用来做AVO分析或者剩余速度分析。

综上，与现有技术相比，本发明一种VTI介质中叠前时间偏移方法的所具有的优点是：

(1)是Offset平面波数据的成像，可以具有很高的计算效率；

(2)是有限差分法成像，使用的是时间域层速度不是均方根速度，可以适应缓横向变速；

(3)是波动方程成像，可以得到较好的振幅保真成像效果；

(4)可以方便地产生用于AVO/AVA分析及偏移速度分析的角度道集；

(5)可以很方便地处理吸收衰减补偿；

(6)可以进行VTI介质的叠前时间偏移，实际数据偏移成像效果明显优于类似的Kirchhoff积分偏移方法。

附图说明

附图是本发明VTI介质中叠前时间偏移方法流程图。

具体实施方式

下面结合具体实施方式对本发明作进一步详细地描述。

如附图所示，本发明一种基于波动方程外推算子的VTI介质中叠前时间偏移方法，包括以下步骤：

步骤1：输入CMP道集以及时间域速度场；

步骤2：对CMP道集进行局部倾斜叠加，将CMP道集映射到截距时间-射线参数域的平面波数据体；并按照射线参数对平面波数据体分选，形成共射线参数剖面；

步骤3：对射线参数循环；

步骤4：进程判断所在计算机节点是否为主控节点，由此进入主-从并行模式；

即：若进程在主控节点上运行，则执行下述(5-A)过程，若进程未在主控节点上运行，而是在计算节点上运行，则执行下述(5-B)过程；

(5-A)首先，读取当前射线参数对应的共射线参数剖面，然后沿时间方向作一维快速傅立叶变换，生成频率域的共射线参数平面波数据体；与此同时，设定总任务数为总频率数目与计算节点数目之中的最小值；主控节点控制着并行任务的分发与收集，并执行以下步骤：

(5-A1)对总任务数循环；

(5-A2)将一个频率切片的平面波数据体发给计算节点，发送任务的计数器自增；

(5-A3)判断任务发送是否结束，若否，转向(5-A1)；若是，转向(5-A4)；

(5-A4)对总频率数循环；

(5-A5)接收计算节点返回的任务标识符；

(5-A6)检查已经发送给计算节点的任务序号是否小于总频率数目；若是，转向(5-A7)；若否，转向(5-A8)；

(5-A7)将当前频率切片的平面波数据体发给计算节点，发送任务的计数器

自增，转向(5-A4)；

(5-A8)向所有计算节点发送任务结束标识符；

(5-B)首先，判断该计算节点号是否大于总的频率数目；若是，转向步骤6；若否，转向(5-B1)；

(5-B1)进程进入死循环；

(5-B2)进程等待接收从主控节点发送来的数据；

(5-B3)进程判断主控节点是否发送结束标识符；若是，转向(5-B9)；若否，转向(5-B4)；

(5-B4)沿时间方向进行波场外推；

(5-B5)计算偏移距平面波传播的相速度求解如下隐式方程：

$>f\left(\theta \right)=2\frac{\sin \theta }{v_p(\theta ;v_{p0},ϵ,\delta )}-p_h=0---\left(1\right)$>

公式(1)中，平面波相速度v_p(θ；v_p0，ε，δ)可以表示为：

$>\frac{{V_{\mathrm{qP}}}^2\left(\theta \right)}{{v_{p0}}^2}=\frac{1}{2}+ϵ{\sin }^2\theta +\frac{1}{2}\sqrt{{(1+2ϵ{\sin }^2\theta )}^2-2(ϵ-\delta ){\sin }^{2}2\theta }---\left(2\right)$>

公式(2)中，ε、δ为VTI介质各向异性参数，θ为相角，v_p0为VTI介质中沿对称轴(Z轴)方向地震波的相速度；通过公式(1)～(2)，可以得出由射线参数p_h表征的偏移距平面波传播的相速度v_p(p_h)；

(5-B6)对VTI介质中的偏移距平面波方程采用隐式有限差分格式进行波场外推，实现偏移距平面波叠前时间偏移；所采用的偏移方程为：

$>(\frac{4}{v_p^2\left(p_h\right)}\frac{{\partial }^2}{{\partial \tau }^2}+\frac{{\partial }^2}{{\partial x}^2})\tilde{U}(x,\tau ;\omega )+4\frac{{\omega }^2}{v_p^2\left(p_h\right)}(1-\frac{v_p^2\left(p_h\right){p_h}^2}{4})\tilde{U}(x,\tau ;\omega )=0---\left(3\right)$>

公式(3)中，为频率域平面波数据，ω为角频率，v_p(p_h)为(2)中计算出的平面波传播的相速度；

采用零时间成像条件，得到共射线参数域成像道集；

$>I\left({x,\tau ;p}_h\right)=\underset{\omega }{\int }\tilde{U}(x,\tau ;\omega )·{\tilde{U}}^{*}(x,\tau ;\omega )\mathrm{d\omega }---\left(4\right)$>

(5-B7)判断波场外推是否结束；若否，转向(5-B4)；若是，转向(5-B8)；

(5-B8)向主控节点发送任务结束标识符，转向(5-B1)；

(5-B9)所有计算节点上的进程向主控节点规约当前射线参数对应的成像数据体，规约操作符号为求和；

步骤6：所有进程同步；

步骤7：主控节点将当前射线参数对应的成像数据体输出；

步骤8：判断射线参数循环是否完成；若否，转向步骤3；

步骤9：输出时间偏移剖面以及成像道集。

尽管上面结合图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以作出很多变形，这些均属于本发明的保护之内。