CN106133556A - 确定波场的分量 - Google Patents

Abstract

描述了涉及一种确定地球的各向异性地下地层中的波场的方法的实施例。所述方法包括基于空间变化的各向异性参数对去耦的准声学单波模式波动方程进行数值求解，以确定各向异性地下地层中的波场。

Description

确定波场的分量

技术领域

本发明涉及确定地球地下地层、尤其是各向异性地下地层的波场或波场分量。

背景技术

可以通过“波动方程”描述波通过介质的传播。通常可以求解波动方程来确定波的作为时间和空间的函数的幅度或其它性质。

可以通过执行地球物理勘测来对地球的地下进行勘探。例如，可以执行地震或电磁勘测，其中，从源向地下地层发射声波或电磁(EM)波，并使用探测器来检测响应。可以基于所检测到的响应来确定地球的性质。例如，所述响应可以包括具有与地下的界面相关联的特征的时间序列数据。

在地震反射勘测中，地震数据中的高振幅事件可以与岩石性质发生突变的界面有关，并且可以将显著的量的能量反射回探测器。因此，数据中的高振幅可以指示地下地层界面。

在地震数据的典型处理中，往往将地震波的发射与高振幅事件的检测之间的时间看作是地震波抵达处于源与接收器之间的中点处的界面并因反射而返回的行进时间。因此，通常将行进时间看作是深度的代表，因此将地震数据标绘为行进时间的函数可以指示振幅特征的深度关系，继而指示地下地层界面的深度关系。

然而实际上，波穿过地下层的地震速度不是恒定的(通常随深度而增大)，因此为了获得对界面位置的更加实际的表达，将基于预先确定的速度分布图或速度模型来对数据进行深度转换，即将其标绘为深度的函数。

在源与接收器之间的中点处对源-接收器对的数据进行标绘。因此，来自所有这样的对的数据可以呈现为该中点处的时间或深度记录，其中，振幅事件所呈现出的样子就像是从正入射源-接收器布置导出的一样。尽管这能够给出对各层为水平且平面的地球结构的适当表示，但是应当认识到地震波可能沿其传播路径受到折射，受到倾斜界面的反射，以及与入射和反射方向内的不同层发生相互作用。因此，行进时间数据和深度转化数据记录可能表明错误的几何结构。例如，可能在特定深度位置处的数据中出现高振幅事件，而实际上在该深度位置处根本没有反射体。

这是在地震数据处理和对地震反射数据的解释中，尤其是在地质条件复杂的结构中已知的问题。可以通过偏移过程来尝试对数据进行校正。地震偏移用来将地震数据重新定位到其在数据集中的代表性的、几何结构正确的位置。已经开发出了用于对地震数据进行偏移的偏移算法。理解地震波通过地下地层的波传播以执行偏移是有用的。简言之，通过跟踪波从源到反射体再返回到地表所采取的路径，有可能识别出与该反射体相关联的能量所抵达的接收器以及该能量的抵达时间。为了达到该目的，需要地下地层的速度模型。该模型可以由模型建立程序包获得，并且可以基于井数据或者可以通过其它方式获得，例如，或许可以通过基于区域地质知识来进行估算。

另一个复杂化的因素是地下地层的各向异性。例如，特定的岩石类型可能具有很强的传播方向性。不同的(例如，正交的)方向内的地震速度可能存在显著的差异，而且可能在各层之间存在差异。

已经开发出了全波场偏移技术，以及时地并在3D空间中的跨越地下地层区域的密集位置上对波的发展进行建模。这样的技术可能涉及求解上文提及的种类的波动方程，从而关于与不同位置上的源信号发射有关的时间来获得不同位置上的波场。

这些偏移算法以标准地震数据作为输入。可以将地震数据作为共炮点道集加以组织，或者在诸如平面波或延迟炮点域的其它域内对地震数据加以组织。可以在经变换的域(与最初采集的地震数据相比)内，例如在针对延迟炮点偏移的炮点位置的线性相位的组合域内提供地震数据，在所述域内通过对炮点位置进行tau-p变换来执行数据的组合。也可以对接收器位置执行所述tau-p变换，以执行平面波偏移。波动方程的波场解在每一目标位置/点处提供了一组复值，从而提供(例如)振幅、相位、波数等。这些值是用于偏移算法的关键输入分量。

已经开发出了很多不同的波动方程技术，而且感兴趣的是开发出在各向异性的地下地层的情况下在地震偏移中使用的波动方程技术。现在更详细地描述一些已知的偏移技术。

当今，准确的地震各向异性模型被认为是对复杂结构的深度挑战区域进行成像的关键点，因为这些模型描述地震波的方向性传播速度的差异，这对于浮盖层地质结构(例如具有断裂的层里沉积物)而言要实际得多。倾斜的横向各向同性(TTI)或倾斜的正交各向异性(TOA)通常是对这种复杂的各向异性结构进行成像所必需的(Zhang and Zhang,2008；Zhang and Zhang,2011；Fowler and King,2011)。假定准确的各向异性模型，偏移算法注定提供无伪影的地下图像。这可能需要应用最小二乘法偏移算法(Lambaré等,1992；Nemeth等,1999)，并且它在从事于对弹性模型进行成像时甚至可能变得更加复杂(Jin等,1992)。已知当前的基于积分的偏移算法(例如基尔霍夫偏移)是非常有效的，因为射线跟踪可以容易地适应于各向异性的程函方程(Han and Xu,2012)。这些算法被开发为可应用于单模式波，例如准P波，这允许使偏移算法线性化，以便降低计算成本。然而，对于基于波动方程的算法，对诸如准P波或准S波等单模式波进行有效率且准确的模拟是开发成功的偏移算法(例如，逆时偏移(RTM))所必不可少的(Baysal等,1983；McMechan,1983；Whitmore,1983)。

历史上，存在不同的方式对单模式波进行数值模拟。这些方式中的第一种方式是求解全弹性波动方程，之后对波场进行拆分，以分离出准P波来进行进一步分析。波场分离对于各向同性情况可能是有效的(Sun等,2004)，但是对于异质各向异性情况则并不是容易的任务，其可能导致计算效率极为低下(Yan and Sava,2009；Cheng and Fomel,2013)。

另一种在仍然保持横向各向同性(TI)的各向异性波传播的同时降低计算成本的方式是通过将沿对称轴的切变波速度设为零来对TI方程应用声学波近似(Alkhalifah,2000)。这得到标量四阶微分方程。然而，Alkhalifah(2000)并未提出用于由该四阶微分方程获得波场的任何可工作的数值技术。Zhou等(2006)将该四阶微分方程分解成2×2二阶微分方程的耦合方程组。这产生了比原始弹性方程更具计算效率的方案。然而，对具有可变TTI对称轴的TTI介质的应用(尤其是在TTI对称轴上存在突变的情况下)证实该方程组不是数值稳定的：出现了弱不稳定性，并且噪声随时间线性增长(Liu等,2009,Zhang等,2011)。这些不稳定性已经得到了很好的分析，并且也提出了解决方案(Bakker and Duveneck,2011；Zhang等,2011；Bube等,2012)。尽管能够解决不稳定性，但是将沿对称轴的切变波速度设为零不排除伪切变波的存在，伪切变波是该耦合方程组的固有解。人们已经尝试了各种方式来衰减伪切变波(Zhang等,2009；Guan等,2011)；或者简单地忽视这些切变波伪影的存在，就像TTI RTM的当前行业实践所做的那样，并且预期通过成像条件和偏移叠加来抵消这些伪影。

另一种方式是直接计算解耦的准声波P波方程。假定恒定的各向异性参数ε和δ，可以利用伪谱算法求解解耦的方程(Etgen and Brandsberg-Dahl,2009)。因此，其允许我们采用具有不同的恒定ε和δ的一组解，从而继续进行内插方案，以在数值上求解出具有在空间域内和缓地变化的各向异性参数的准P波(Chu等,2013)。对于VTI(具有垂直对称轴的横向各向同性)介质，内插方案关于ε和δ是二维的。但是对于TTI介质，如果对称轴还发生空间变化，那么内插需要另外的两个维度。该伪谱与内插结合的方法可能需要计算很大的一组各向异性参数的组合，并且在应用于非常复杂的模型时效率可能非常低。在各向异性模型的复杂度不高时，可以通过近似优化的可分离方法(Liu等,2009)来求解所述去耦方程，该方法可以被看作低秩近似的特例(Fomel等,2012)。计算成本可能比伪谱与内插结合的方法更低，但是效率仍然很低，而且不准确。

发明内容

提供了如在所附权利要求中阐述的本发明的各个方面。

任何方面都可以具有联系任何其它方面所描述的进一步的特征，并且可以使各个方面之间的特征组合以及互换。

有利地，本发明的实施例提供用于对在各向异性介质中传播的波的波场进行数值计算的技术。具体而言，就降低的数值或计算成本以及实施的简便性而言所述技术是有利的。

现有的技术通常使用描述波传播的耦合方程组，并为几个波模式一起提供波场解。通常使用针对准P波和准S模式的耦合方程组。这样的方程组固有地具有分别与准P波和准S波模式相关联的两个特征值。在现有的耦合方程技术中，发现由于包含准S波模式的原因而在波场解中存在伪影和噪声。耦合方程解总体上具有高数值成本。

一直缺乏用于对在地球的各向异性区域中传播的地震波的波场进行数值计算的技术。

附图说明

现在将参考附图仅通过举例的方式描述本发明的实施例，在附图中：

图1是根据本发明的实施例的确定波场的方法的流程图表示；

图2是根据本发明的实施例的用于确定波场中的计算机装置的示意性表示；

图3是用于采集数据的设备的示意性表示；

图4是VTI模型波场解的波场幅度的等值线图表；

图5是用于计算图4的模型波场的算子的幅度的等值线图表；

图6是图4和图5的波场中的复算子的幅度的等值线图表；

图7是TTI模型波场解的波场幅度的等值线图表；

图8是用于计算图7的模型波场的算子的幅度的等值线图表；

图9是图7和图8的波场中的复算子的幅度的等值线图表；

图10是SEAM TTI盐丘模型的速度模型的等值线图表；

图11是在时刻t＝1.6s的图10的模型的波场幅度的等值线图表，

图12是在时刻t＝2.4s的图10的模型的波场幅度的等值线图表，

图13示出了三条频散曲线，

图14示出了渐近校正项的振幅，

图15示出了两个波场，以及

图16示出了利用椭圆模型执行的计算。

具体实施方式

波动方程模型

我们首先针提出用于复杂的各向异性模型的纯粹的影片上准确的(cinematically-accurate)准P波方程。通过将Alkhalifah(2000)的原始伪微分算子分解成两个数值可解项：微分算子和标量算子来获得该方程。所得到的新的准P波方程将各向异性模型参数与空间伪微分算子分离，因而能够借助于传统的数值方案来求解所述方程而不牺牲其准确度，即，并未应用任何近似，因此其保持了与初始伪微分波动方程相同的频散关系。方程中的微分算子被设计为自共轭，因此其能够在波传播时节约能量。之后，我们讨论所提出的算法的数值实现。所提出的准P波方程具有与声波方程类似的形式并且能够被类似地实施。在我们的当前实施方式中，我们将伪谱算法用于空间导数，并且将二阶准确度有限差分方案用于时间导数。新算法向TTI和TOA介质的扩展简单明了，并且在文中已经简要地提到过。

我们借助于VTI/TTI模型的脉冲响应的示例对我们的方案进行了验证，并且示出了新的纯粹的准P波方程中的项的功能性。最终，我们借助于SEAM TTI盐丘模型中的准P波传播的示例来证实了我们的准P方程的有效性和鲁棒性。

伪微分算子的分解

在下文中，为简单起见我们将首先研究VTI情况。这要求我们求解下述标量伪微分方程(Alkhalifah,2000)：

$\frac{{\∂}^{2}u}{\∂t^2}-\frac{v_0^2}{2}\left(\frac{{\∂}^2}{\∂h^2}\right(1+2\mathrm{\ϵ})+\frac{{\∂}^2}{\∂z^2}+\sqrt{{\left(\frac{{\∂}^2}{\∂h^2}\right(1+2\mathrm{\ϵ})+\frac{{\∂}^2}{\∂z^2})}^2-8(\mathrm{\ϵ}-\mathrm{\δ})\frac{{\∂}^2}{\∂h^2}\frac{{\∂}^2}{\∂z^2}})u=0---\left(1\right)$

其中，是水平拉普拉斯算子；v₀是沿对称轴的速度；ε和δ是Thomsen(1986)各向异性参数，如在Thomsen(1986)中所定义的。通过引用将Thomsen(1986；参见下文的参考文献部分)内提供的对这些参数的定义并入本公开内容。注意，所有的这些参数都是空间变化的。方程1控制准P波传播，但是其解的振幅可能与实际的弹性准P波相当不同。这是因为全弹性方程的解耦，在解耦之后省略了不同波模式的所有相互作用，包括与其它波模式之间的相互转化。然而，方程1共享与弹性准P波的频散相同的频散，这指示其相位与弹性准P波相比是准确的。尽管方程1对于单模式波传播而言是非常有前景的，但是令人遗憾的是，其代表不能用传统的有限差分数值方案来求解的伪微分算子方程。Alkhalifah(2000)未在复杂介质中以数值方式求解方程1。

该伪微分算子的公式表明如下的数值解，其中，首先确定当前波场的空间梯度分量，形成了平方根内的梯度分量的组合，之后取得该组合的平方根。这在传统的方案中行不通。我们改写该方程，并用两个不再是伪微分算子的算子来代替伪微分算子。

为了描述我们的方案，我们从它的方程1的对应频散关系开始，该频散关系被如下表示：

$({\mathrm{\ω}}^2-\frac{v_0^2}{2}(k_h^2\left(1+2\mathrm{\ϵ}\right)+k_z^2+\sqrt{{\left(k_h^2\right(1+2\mathrm{\ϵ})+k_z^2)}^2-8(\mathrm{\ϵ}-\mathrm{\δ})k_h^2{k}_z^2}\left)\right)u=0.---\left(2\right)$

这里，ω是角频率。空间波数向量就像通常情况那样被定义为因此，k_h是水平波数，其中，k_z是垂直波数。还将方程2改写为：

$({\mathrm{\ω}}^2-\frac{v_0^2{\stackrel{\→}{k}}^2}{2}(n_h^2\left(1+2\mathrm{\ϵ}\right)+n_z^2+\sqrt{{\left(n_h^2\right(1+2\mathrm{\ϵ})+n_z^2)}^2-8(\mathrm{\ϵ}-\mathrm{\δ})n_h^2{n}_z^2}\left)\right)u=0,---\left(3\right)$

其中，是相位方向的单位向量，并且被定义为

让我们将辅助标量算子S定义为

$S=\frac{1}{2}\left(n_h^2\right(1+2\mathrm{\ϵ})+n_z^2+\sqrt{{\left(n_h^2\right(1+2\mathrm{\ϵ})+n_z^2)}^2-8(\mathrm{\ϵ}-\mathrm{\δ})n_h^2{n}_z^2}).---\left(5\right)$

现在，方程3变为

$({\mathrm{\ω}}^2-v_0^2{\stackrel{\→}{k}}^{2}S)u=0---\left(6\right)$

在空间域中，将算子表达为其中，符号表示散度，表示梯度。将方程6从频率—波数域转换回到时间—空间域得到了下述方程：

$\frac{{\∂}^{2}u}{\∂t^2}-\▿\·(v_0^{2}S\▿u)=0.---\left(7\right)$

该方程表示方程6在时间—空间域中的对应偏微分方程。其定义了在各向异性介质(例如地球的地下地层)中的准P波传播的模型，可以用易于实施并且有效率的数值方法来求解该方程，以在空间和时间中提供波场。其不再是伪微分方程，因为我们已经从导数中分离出了各向异性项。但是该方程是非线性波方程，因为标量算子S还取决于波场的解。可以在空间域中计算该标量算子，因为向量具有其物理含义：波前的相位方向。由于方程5对的分量的符号不敏感(因为所有的分量都具有奇数次幂)，我们可能仅使用的分量来代替方程5中的的分量。标量算子S在我们的方案中相当关键；其实际上起着整个波传播的频散关系的作用。其控制准P波的传播速度(各向异性传播算子)，其不仅取决于所述传播在每个空间点处的相位方向，还取决于每个个体空间位置处的各向异性参数。并且该算子的计算仅需要存在空间变化的各向异性参数、以及当前波场的梯度。

可能要注意，以时间步长来计算波传播。为了计算出下一时间步长的波场，需要当前波场。定义起始波场，为计算提供了边界条件或起始点。可以定义在震源位置处发射的预定义震源子波，并使用该预定义震源子波来允许确定该起始波场。

速度模型是预先确定的，例如是由其它算法估算的和/或是在模型建立程序包中提供的。模型建立程序包可以包括很多算法，例如层析反演、全波反演、盐丘解释等。速度是所需的输入，并且其对于偏移算法而言被认为是已知的。各向异性参数是波传播模型(输入)的部分，并且还包含在速度模型中。可以通过模型建立程序包来估算这些参数。

方程7中的微分算子是自共轭算子，并且因此节约能量。其应该确保波传播的稳定性，即使对于模型和各向异性参数存在突然的空间变化的情况而言仍然如此。

在我们的方案中，从VTI向TTI或TOA的推广是非常简单的。对于TTI而言，我们只需将波场的梯度向量投射至局部坐标，其中，局部z轴是各向异性的对称轴；之后向方程5应用完全相同的过程。可以将类似的分析推广至TOA介质。

来自方程6的替代的解可以用其如下所述的方程来实施：

$\frac{{\∂}^{2}u}{\∂t^2}-v_0^{2}S{\▿}^{2}u=0---\left(8\right)$

这里，将方程1中的伪微分算子分解成两个算子：标量算子和拉普拉斯算子。该方程直接使用拉普拉斯算子，并且具有与声波方程类似的形式。

分解成球面项

在波数域中对标量算子S的数值计算要求方程7中的各向异性参数是恒定的，并且因此可以在对相位方向进行渐进近似的情况下在空间域中计算所述的数值计算。

为了对渐进近似项的作用进行理论分析，我们将方程7中的标量算子改写为：

$S=1+\mathrm{\Δ}S,\mathrm{\Δ}S=\frac{1}{2}(\sqrt{{\left(n_h^2\right(1+2\mathrm{\ϵ})+n_z^2)}^2-8(\mathrm{\ϵ}-\mathrm{\δ})n_h^2{n}_z^2}+2{\mathrm{\ϵn}}_h^2-1).---\left(9\right)$

使用标量算子的这一表达式，波动方程变为

$\frac{{\∂}^{2}u}{\∂t^2}=v_0^2{\▿}^{2}u+v_0^2\▿\·(\mathrm{\Δ}S\▿u).---\left(10\right)$

方程10中的右侧第一项为背景波动方程。可以将其看作微分算子——拉普拉斯算子，其不包含近似。可以将方程10的右侧第二项看作校正项。ΔS的计算取决于波传播方向，该计算为渐进近似。可以将方程10看作初始波动方程1的球面分解。

在图13中，频散曲线131对应于方程1的解，频散曲线132是对应于方程10的右侧第一项的背景频散曲线。下文将联系椭圆模型讨论频散曲线133。针对恒定的各向异性参数ε＝0.25和δ＝0.1示出了背景频散132。渐近校正旨在将相位从曲线132校正至曲线131。

在图14中，曲线141示出了0～360°的角度之间的渐进校正ΔS的幅度。注意，与背景频散相比，校正相对较大：对于该示例，最大校正约为背景频散的50％。因此，该方案要求非常准确地估算波场的方向。可以使用比奈奎斯特抽样大的对计算的横向空间抽样。因此，波场的高波数部分被混叠，并且可能在波场的较低波数范围内的方向向量的计算中引入误差，这使得该方案易受方向误差的影响。

分解成椭圆项

为了提高算法中的方向向量的数值误差的容限，我们提出了在下面的方程11中定义的椭圆分解，以代替球面分解。在该方程中，我们仍然将伪微分算子分解成两个算子：微分算子和标量算子。然而，用椭圆微分算子代替原始分解中的拉普拉斯算子，同时还对标量算子进行相应地修改，以确保波传播的准确相位。新的分解的目的在于减小渐近项的幅度。该分解为：

$\frac{{\∂}^{2}u}{\∂t^2}=(v_0^2\frac{{\∂}^{2}u}{\∂z^2}+(1+2\mathrm{\ϵ}\left)v_0^2\frac{{\∂}^{2}u}{\∂x^2}\right)S_e,S_e=\frac{1}{2}(1+\sqrt{1-\frac{8(\mathrm{\ϵ}-\mathrm{\δ})n_x^2{n}_z^2}{{\left(\right(1+2\mathrm{\ϵ})n_x^2+n_z^2)}^2}})---\left(11\right)$

其中，S_e为椭圆标量。项可以被写作并被解释为微分算子乘以u，并且该微分算子(大括号内的项)是椭圆微分算子。为了进一步分析方程5中的渐近项，我们将其明确地改写为(方程12)：

$\frac{{\∂}^{2}u}{\∂t^2}=v_0^2\frac{{\∂}^{2}u}{\∂z^2}+\left(1+2\mathrm{\ϵ}\right)v_0^2\frac{{\∂}^{2}u}{\∂x^2}+\left(v_0^2\frac{{\∂}^{2}u}{\∂z^2}+\left(1+2\mathrm{\ϵ}\right)v_0^2\frac{{\∂}^{2}u}{\∂x^2}\right){\mathrm{\ΔS}}_e,{\mathrm{\ΔS}}_e=\frac{1}{2}(\sqrt{1-\frac{8(\mathrm{\ϵ}-\mathrm{\δ})n_x^2{n}_z^2}{{\left(\right(1+2\mathrm{\ϵ})n_x^2+n_z^2)}^2}}-1)$

图13示出了球面方案与椭圆方案之间的频散曲线的比较。对于恒定的各向异性参数ε＝0.25和δ＝0.1，曲线133示出了方程12的右侧第一项的频散曲线。将该曲线与球面分解方案进行比较，我们注意到，椭圆分解背景更加接近被示为曲线131的期望的精确解。在图14中，线142示出了0～360°内的渐进校正ΔS_e的幅度。ΔS_e的最大幅度为0.068，其比球面分解小7倍。因此，椭圆分解对方向误差的容限要高得多。

首先，我们用图15中的简单脉冲的示例来证实该算法的效果。该示例是简单的TTI(倾斜的横向各向同性)模型，该模型是同质的，其中，在倾斜角为30°并且方位角为135°的情况下将垂直速度定义为2000m/s，并且将Thomsen(1986)各向异性参数定义为ε＝0.24和δ＝0.1。该示例的震源子波是具有24Hz的最高频率的雷克子波。计算网格是具有6.0km的边长的3D立方体，并且空间抽样在所有的3个方向内均为15m。将点源置于网格的中央。图15a)标绘了根据球面分解方案的在时刻t＝1.0s处的3D波场快照的位于Y方向的中间的2D地震片层；图15b)标绘了利用椭圆分解方案获得的在时刻t＝1.0s处的3D波场快照的位于Y方向的中间的2D地震片层。两种方案仅生成了伪P波，而不存在切变波。两波场取得了相同的传播相位，但是椭圆分解方案给出了更加均衡的振幅。我们还注意到，球面分解和椭圆分解的数值成本几乎等同。

第二示例是利用SEG SEAM模型的偏移测试。我们为该测试选择了炮线。图16在左侧示出了该炮线的位置，该炮线包含342炮。我们利用所提出的纯粹的准P波方程和发生偏移的这些炮来构建TTI RTM。图16的右侧示出了通过使用文中公开的椭圆算法所得到的图像结果，该结果与密度模型叠加。显然，用新的准P波方程生成的图像与密度模型匹配得非常好，并且呈现出了清晰且准确的结果。

数值法

对方程7进行数值求解，以在不同的波传播时间处获得地下地层内的每个预定义位置的波场分量。这是通过数值估算过程完成的。求解方程7的数值过程的实施起来相对简单直接。方程10和11的数值解相似，但是标量的值是不同的，因为在与方程7相比时微分算子是不同的。所述过程具有用于确定波场的步骤S1到S3，如下文所阐述以及图1所图示的：

S1.计算当前波场的梯度。例如，当前波场是针对前一时间步长所确定的波场。其最初为边界时间零点处的波场。

S2.计算标量算子S。这是按照方程5中指示的完成的。通过使用来自步骤S1的梯度场来在每个空间位置处计算算子S，以获得相位方向分量(在投射至局部坐标之后忽略符号)。将标量算子和速度的平方与波场的梯度相乘。这些全部是数值。

S3.计算来自步骤2的结果的散度，以确定在给定时间步长处的波场。获得来自当前波场的时间变化率的值，该值又用于典型地通过使用积分法来给出新时刻处的场。

步骤S3中确定的波场用作接下来的传播时间步长的步骤S1中的当前波场，如图1中的环所指示的。在相继的传播时间步长处重复步骤S1到S3，使得能够确定关于时间和空间的准确的波场。

从图1中可以看出，在该过程中可以存在用于提供初始波场的初始步骤S0。在图1中还示出了另一步骤S4，在该步骤中使用所确定的波场使地震数据偏移。一旦针对所有时间步长完成了每个空间位置处的波场计算就进行该步骤。将标准偏移算法布置为使用波场或其分量。可以由上述求解方法获得针对预期行进时间和感兴趣的地下地层位置的适当波场。

类似地，可以通过具有下述步骤T1和T2的数值估算过程求解方程8。

T1.计算当前波场的梯度。例如，当前波场是针对前一时间步长所确定的波场。其最初为在时间零点处的边界处的波场。

T2.通过使用来自步骤T1的梯度场来针对每个空间位置计算在方程5中看到的标量算子S，以获得相位方向分量(在投射至局部坐标之后忽略符号)。将标量算子和速度的平方与波场的拉普拉斯梯度相乘。获得了时间变化率的值，该值又用于典型地通过使用积分法来给出新时刻处的场。

可以使用标准有限差分算法，或者替代地使用标准快速傅里叶变换(FFT)法来执行上述方法的在相关步骤S1到S3或者T1到T2内的梯度计算。

在使用有限差分数值技术时确定波场的成本最多是使用标准各向同性声波方程解的二倍。这种方法比各向异性波传播中的其它方案快得多。

为了获得上面的步骤S1到S3中的梯度和散度，必须计算波场的一阶导数(del和散度)，并且可以使用优化的数值方案来同时计算这些一阶导数，这样可以是有效率的。类似地，为了获得方程8中的拉普拉斯(del的平方)，需要一阶导数和二阶导数，并且可以使用优化的数值方案，其中，可以同时且有效率地计算波场的一阶导数和二阶导数。例如，可以使用快速傅里叶变换(FFT)来计算这些空间导数。在该情况下，其只需要一次正向FFT和两次反向FFT就能同时获得一阶导数和二阶导数两者。而且，标准FFT算法是使用的。

在使用FFT数值技术时，以仅比针对各向同性介质的标准声波方程的解多50％的额外数值成本执行对波场的确定。我们注意到，在从VTI移至TTI或TOA介质时数值成本最低限度地增大。

方程8直接使用拉普拉斯算子。因此，有效率的数值方案可以更易于实施。与方程7相比，方程8具有相同的运动行为，但是具有不同的振幅效应。

用于获得步骤S3或T2中的波场的积分法可以是标准的时间积分数值法，例如瞬时积分的二阶准确度的有限差分方案或快速展开法(REM)(Kosloff等，1989)。

总之，波动方程(方程8)具有的形式。在数值上，波场取决于空间点X和时间t，可以将其表达为U(x,t)。求解波场的任务涉及使用当前时间样本处的波场来计算下一时间样本的波场。因而，对于方程的左侧(关于时间的二阶导数)，我们使用的方法可以是被我们称作积分法的标准方案。对于该积分的计算，我们必须知道方程的右侧的值，其包含波场的空间导数(方程7中的一阶导数和方程8中的二阶导数)。可以使用为有限差分(FD)或FFT的(FD更有效率，而FFT则更加准确)常规数值方案来用于该目的，同时注意，总体上FD和FFT两者都是标准的，并且凭其自身的能力是很好理解的算法。

现在参考图2，描绘了用于如上文所述执行确定波场或其分量的方法的计算机装置。计算机装置10具有输入/输出装置11、微处理器12和存储器13。计算机程序14a和14b存储在存储器13中。波场计算机程序14a具有用于执行求解方程7和方程8的数值法的指令，以获得不同传播时间和空间位置处的波场。微处理器12被布置为读取并执行波场计算机程序中所包含的指令，以确定波场。所计算出的波场还优选地存储在存储器中，并作为输入被传送至也存储在存储器中的偏移程序，所述偏移程序也可由处理器12运行，以便使用所计算出的波场来执行地震数据的偏移。所述装置还可以具有用于查看存储在处理器中和/或经由程序计算出的数据的显示器。可以如上文(包括背景技术部分)所述地执行偏移。

输入/输出装置11用于向读取计算机装置中读入数据或者从计算机装置输出数据。具体而言，可以通过输入/输出装置接收可能在地震勘测中获得的地震数据，并且将这样的数据存储在存储器13中。

在图3中，描绘了包括地震勘测设备1的设备。所述设备包括通过水体牵引地震源4和地震探测器5的地震勘测船。使用地震源发射通过地下地层2的地震波。所述波与界面3相互作用，并将部分能量朝探测器反射回来。探测器被布置为检测在探测器处接收到的能量。通常可以使用探测器获得包含相对于与生成地震波的震源事件有关的行进时间的振幅记录的数据。高振幅事件则可以与来自地下地层界面的反射相关联。可以通过计算机装置读取来自勘测的数据，并对其进行处理，以提供地下地层图像，从而有助于揭示地质构造。例如，所述设备可以包括图2中描述的计算机装置，并且可以由该计算机装置读取来自勘测的数据并使用处理器来处理该数据。然后使数据发生偏移，以提供偏移的地震勘测数据报告或图像，例如偏移的地震剖面。可以应用这样的系统来提供在上文中的别处描述的数据。

计算机装置可以是分布式装置，因为输入/输出装置11、微处理器12、存储器13和显示器15中的任何一个或多个可以分布在不同位置。其间的通信可以是如所指示的通过网络(例如，无线网络)发生的。在某些实施例中，程序、波场数据和/或偏移的数据可以存储在可移除存储介质(例如存储棒或光盘)中，可以在与计算机装置和/或处理器连接时由它们来运行。可以提供经由网络来传递的信号，该信号包含上文所述产生的程序、其机器可读指令、波场数据和/或偏移的数据。

可以在需要使用所预测的单模式波的波场的其它应用(例如，地震建模、全波形反演等)中使用当前提出的波动方程7和波动方程8的波场解。

示例和结果

已经对上文描述的方程7和方程8的波动方程模型执行了各种测试，从图4到12中可以看到其结果。在所有的下述示例中，我们对空间导数应用伪谱技术，对时间导数应用具有二阶准确度的有限差分方案。也就是说，将FFT(由于其具有更高的准确度)用于空间导数，而将具有二阶有限差分方案的积分法用于时间导数。

选择数值网格尺寸，以避免空间频散，即，其最高可达奈奎斯特波数，并且选择传播时间步长以满足稳定性条件。

第一示例是简单的VTI模型，其为同质的，并且垂直速度被定义为2000m/s，Thomsen(1986)各向异性参数被定义为ε＝0.2和δ＝0.1。该示例的震源子波为具有24Hz的最高频率的雷克子波。计算网格是具有6.0km的边长的3D立方体，其中，在所有的3个方向内空间抽样为15m。将点源置于网格的中间。图4标绘了在时刻t＝0.8s处在3D波场快照的位于Y方向的中间处的2D地震片层。显然，只出现了伪P波并且不存在切变波。图5示出了在同一快照时刻处的算子S的对应图像。该算子在我们的算法中起着关键的作用。我们观察到该图并不像波场那样平滑、清晰。这是因为由表示的传播方向在梯度函数趋于零的位置周围丧失了其准确性。幸运的是，可以通过复合算子找回该丧失的准确性，该算子是在应用方程7中定义的速度模型之前在散度算子内的项的幅度。在图6中示出了两个量的乘积的组合效果。平滑、清晰的图像证实了，方向的不准确性几乎未向波场计算引入可觉察的误差。

为了针对TTI情况验证我们的方案，我们利用与第一示例中的VTI模型相同的计算参数，并将其扩展到简单的TTI模型中。因此，除了相同的VTI各向异性参数之外，我们还引入了45°的倾角和0°的方位角，即，使各向异性对称轴倾斜，倾斜是由倾角和方位角定义的(就VTI而言，对称轴是垂直的(非倾斜的))。图7示出了在与第一示例相同的位置和时间片层处的快照波场的2D片层。注意波前的旋转是TTI参数的结果。类似地，图8是算子S的对应输出。图9标绘了两个算子的乘积的组合效果。

我们针对SEAM TTI盐丘模型测试了我们的方案。已知模型尺寸为：nx＝ny＝864，nz＝768。网格抽样率在所有的三个方向内为10m。我们将震源位置放在(x,y,z)＝(17,23.0,0.0)km的位置处。我们再次使用雷克子波作为震源子波，但是这次具有75Hz的最高频率。传播时间步长为0.5ms。图10描绘了3D速度模型的中线；并且图11示出了t＝1.6s处的3D时间快照的中央线内片层，并且图12示出了t＝2.4s处的时间快照。注意，传播器(propagator)将复杂的波场处理得非常好。其生成发射波、反射波和首波，但不生成任何切变波。

优点

与传统的波动方程相比，我们的方案具有一些显著的优点。首先是方程7和方程8简单。它们对VTI、TTI和TOA保持相同的形式。其次是数值效率。与TTI介质中的波传播相比，常规方案所需的计算资源一般比各向同性波的模拟所需的计算资源多3-5倍，而我们提出的方案只引入50％的额外成本。这是比现有方案有效率得多的方案。此外，其数值性能对于横向各向同性、倾斜的横向各向同性、正交各向异性或倾斜的正交各向异性几乎是相同的。再次是方程的稳定性。与声波情况相似，在我们的新方程中没有出现常规的2×2二阶微分方程组中的TTI的弱不稳定性。我们的解决方案对于非常复杂的模型，例如具有复杂盐丘结构和各向异性对称轴发生突变的上冲结构的模型是数值稳定的。由于仅使用一个微分方程，新提出的方案比常规算法的方案更有效率。

总之，用于确定各向异性地下地层的地震波场的当前方法：

·更有效率

·更加准确

·无切变噪声

·实施简单

·数值稳定

尽管我们只讨论了用于准P波的算法，但是可以容易地推广所提出的方案，以解决准SV波传播问题，乃至弹性衰减问题等。其将来可以在S波成像或转化波成像方面存在价值。

该当前解决方案提供了求解伪微分方程的方式。衰减方程是另一个这种方程的一个示例。

方程的缩写用于表示“方程”，并且术语在文中可互换使用。

可以做出各种修改和改进而不背离文中描述的本发明的范围。

参考文献

Alkhalifah T.1998.Acoustic approximations for processing intransversely isotropic media.Geophysics 63,623–631.

Alkhalifah T.2000.An acoustic wave equation foranisotropicmedia.Geophysics 65,1239–1250.

Bakker,P.M.,and E.Duveneck,2011,Stability analysis for acoustic wavepropagation in tilted TI media by finite differences:Geophysical JournalInternational,185,no.2,911–921.

Bube,K.P.,T.Nemeth,J.P.Stefani,R.Ergas,W.Liu,K.T.Nihei,and L.Zhang,2012,On the instability in second-order systems for acoustic VTI and TTImedia:Geophysics,77,no.5P.T171-T186.

Cheng J.and S.Fomel,2013,Fast algorithms of elastic wave modeseparation and vector decomposition using low-rank approximation fortransverse isotropic media:Expanded Abstracts of 83rd Ann.Internat.Mtg.,Soc.Expl.Geophys.

Chu C.,B.K.Macy and P.D.Anno,2013,Pure acoustic wave propagation intransversely isotropic media by the pseudospectral method:GeophysicalProspecting:61(3),pages 556–567.

Duveneck E.and Bakker P.M.2011.Stable P-wave modeling for reverse-time migration in tilted TI media.Geophysics 76,S65–S75.

Etgen J.T.and Brandsberg-Dahl S.2009.The pseudo-analytical method:Application of pseudo-Laplacians to acoustic and acoustic anisotropic wavepropagation.79th Annual International Meeting,SEG,Expanded Abstracts,2552–2556.

Fomel S.,Ying L.and Song X.,2012,Seismic wave extrapolation usinglowrank symbol approximation:Geophysical Prospecting,Volume 61,Issue 3,pages526–536.

Fowler P.J.,Du X.and Fletcher R.P.2010a.Coupled equations for reversetime migration in transversely isotropic media.Geophysics,75,S11–S22.

Fowler,P.J.and R.King,2011,Modeling and reverse time migration oforthorhombic pseudo‐acoustic P‐waves:SEG,Expanded Abstracts,30,no.1,190-195.

Guan,H.,E.Dussaud,B.Denel,P.Williamson,2011,Techniques for anefficient implementation of RTM in TTI media:81st Annual InternationalMeeting,SEG Expanded Abstracts,3393-3397.

Jin,S.,R.Madariaga,J.Virieux and G.Lambaré,1992,Two-dimensionalnonlinear inversion of seismic waveforms,numerical results:Geophysics,51,1387-1403.

Kosloff,D.,Filho A.Q.,Tessmer E.,Behle A.,1989,Numerical solution ofthe acoustic and elastic wave equations by a new rapid expansion method:Geophysical prospecting；37(4):383-394.

Lambaré,G.,J.Virieux,R.Mandariaga,and S.Jin,1992,Iterative asymptoticinversion in the acoustic approximation:Geophysics,57,1138-1154.

Liu F.,Morton S.A.,Jiang S.,Ni L.and Leveille J.P.2009.Decoupled waveequations for P and SV waves in an acoustic VTI media.79th AnnualInternationalMeeting,SEG,Expanded Abstracts,2844–2848.

Nemeth,T.,Wu,C.,and Schuster,G.T.,1999,Least-squares migration ofincomplete reflection data:Geophysics,64,208-221.

Sun R.,G.A.McMechan,H.Hsiao and J.Chow,2004.Separating P-and S-wavesin prestack 3D elastic seismograms using divergence and curl:Geophysics,69,286–297.

Thomsen,L.,1986,Weak elastic anisotropy:Geophysics,51,1954–1966.

Yan J.and P.Sava,2009.Elasticwave-mode separation for VTImedia.Geophysics 74,WB19–WB32.

Zhan G.,R.Pestana and P.L.Stoffa,2013,An efficient hybridpseudospectral/finite-difference scheme for solving the TTI pure P-waveequation:Geophys.Eng.10 025004

Zhang,H.and Y.Zhang,2008,Reverse time migration in 3D heterogeneousTTI media:78th Annual International Meeting,SEG,Extended Abstracts,2196-2200.

Zhang H.,G.Zhang and Y.Zhang,2009.Removing S-wave noise in TTIreverse time migration.79th Annual International Meeting,SEG,ExpandedAbstracts,2849–2853.

Zhang Y.,H.Zhang and G.Zhang,2011.A stable TTI reverse time migrationand its implementation.Geophysics 76,WA3–WA11.

Zhang,H.and Y.Zhang,2011,Reverse time migration in vertical andtilted orthorhombic media:SEG,Expanded Abstracts,30,no.1,185-189.

Zhou H.,G.Zhang and R.Bloor,2006a.An anisotropic acoustic waveequation for VTI media.68th Annual Conference and Exhibition,EAGE,ExpandedAbstracts,H033.

Zhou,H.,G.Zhang and R.Bloor,2006b,An anisotropic acoustic waveequation for modeling and migration in 2D TTI media:76th Annual InternationalMeeting,SEG,Expanded Abstracts,194–198.

Claims (40)

1.一种确定地球的各向异性地下地层中的波场的方法，包括：基于空间变化的各向异性参数对去耦的准声学单波模式波动方程进行数值求解，以确定所述各向异性地下地层中的所述波场。

2.一种确定地球的各向异性地下地层中的波场或波场的分量的方法，所述方法包括：

提供具有基于地下地层各向异性和波场相位方向的分量的算子S，所述相位方向是从预定波场的估算的梯度获得的；以及

使用所述算子S来确定所述波场或者所述波场的所述分量。

3.根据权利要求2所述的方法，其中，所述分量包括以下任何一个或多个：

-振幅；

-相位；

-压力幅度或位移；以及

-电磁场幅度。

4.根据权利要求2或3所述的方法，其还包括以下任何操作：

-估算波场；

-确定所述波场的梯度；

-提供预定义的各向异性参数；

-使用所述梯度和所述各向异性参数来计算所述算子S；

-将所述算子与速度的平方和所述波场的所述梯度组合，以获得组合的结果；以及

-使用所述组合的结果来确定所述波场的所述分量。

5.根据任何前述权利要求所述的方法，其中，所述算子S是标量。

6.根据任何前述权利要求所述的方法，其中，所述算子控制各向异性波传播的不同方向内的传播速度。

7.根据任何前述权利要求所述的方法，其中，基于各向异性的所述分量包括至少一个Thomsen各向异性参数。

8.根据任何前述权利要求所述的方法，其中，基于各向异性的所述分量是基于所述地下地层的预期的各向异性区域而估算出的或者以其它方式预先确定的。

9.根据任何前述权利要求所述的方法，其中，所述波场是压力波场。

10.根据任何前述权利要求所述的方法，其中，所述波场的所述分量是在所述地下地层中的预定位置处确定的。

11.根据任何前述权利要求所述的方法，其还包括使用所确定的所述波场的分量来使地球物理数据偏移。

12.根据权利要求8所述的方法，其中，所述地球物理数据包括地震数据和电磁数据中的任一个或两者。

13.根据任何前述权利要求所述的方法，其中，使用步骤包括求解波动方程。

14.一种处理数据的方法，包括：

(a)提供与地球的各向异性地下地层相关联的第一数据。

(b)对波动方程进行数值求解，以确定所述各向异性地下地层中的波场或波场的分量；以及

(c)使用所确定的波场或波场的分量对所述第一数据进行处理，以产生与所述地下地层相关联的第二数据。

15.根据权利要求14所述的方法，其中，基于在前文定义的方程7或者方程8来执行步骤(b)。

16.根据权利要求14或15所述的方法，其中，使用各向异性传播算子执行步骤(b)。

17.根据权利要求16所述的方法，其还包括：通过将各向异性项与波场相位方向项组合来确定所述各向异性传播算子。

18.根据权利要求17所述的方法，其中，确定步骤还包括估计所述波场相位方向项。

19.根据权利要求18所述的方法，其中，估计步骤包括使用波场梯度来确定所述波场相位方向。

20.根据权利要求14或15所述的方法，其中，使用具有基于所述地下地层的各向异性和波场相位方向的分量的算子来执行步骤(b)，所述相位方向是从波场梯度获得的。

21.根据任何前述权利要求所述的方法，其中，所述波动方程是非伪微分方程。

22.根据权利要求1所述的方法，其中，对所述准声学单波模式波动方程进行球面分解。

23.根据权利要求22所述的方法，其中，所分解的方程包括球面微分算子和球面标量算子。

24.根据权利要求1所述的方法，其中，对所述准声学单波模式波动方程进行椭圆分解。

25.根据权利要求24所述的方法，其中，所分解的方程包括椭圆微分算子和椭圆标量算子。

26.根据任何前述权利要求所述的方法，其中，所述波动方程具有如下形式：

或

27.根据任何前述权利要求所述的方法，其中，所述波动方程是

$\frac{{\∂}^{2}u}{\∂t^2}-\▿\·(v_0^{2}S\▿u)=0$

或者是

$\frac{{\∂}^{2}u}{\∂t^2}-v_0^{2}S{\▿}^{2}u=0,$

其中，u是波场，t是时间，S是标量算子，并且v₀是所述地下地层中的沿各向异性的对称轴的速度。

28.根据权利要求1所述的方法，其中，所述波动方程为

$\frac{{\∂}^{2}u}{\∂t^2}=v_0^2{\▿}^{2}u+v_0^2\▿\·(\mathrm{\Δ}S\▿u),$

其中，S＝1+ΔS,其中，u是波场，t是时间，S是标量算子，v₀是所述地下地层中沿各向异性的对称轴的速度，n_h是水平相位方向，n_z是垂直相位方向，并且其中，ε和δ是各向异性参数。

29.根据权利要求1所述的方法，其中，所述波动方程为

$\frac{{\∂}^{2}u}{\∂t^2}=v_0^2\frac{{\∂}^{2}u}{\∂z^2}+\left(1+2\mathrm{\ϵ}\right)v_0^2\frac{{\∂}^{2}u}{\∂x^2}+\left(v_0^2\frac{{\∂}^{2}u}{\∂z^2}+\left(1+2\mathrm{\ϵ}\right)v_0^2\frac{{\∂}^{2}u}{\∂x^2}\right){\mathrm{\ΔS}}_e,{\mathrm{\ΔS}}_e=\frac{1}{2}\left(\sqrt{1-\frac{8\left(\mathrm{\ϵ}-\mathrm{\δ}\right)n_x^2{n}_z^2}{{\left(\left(1+2\mathrm{\ϵ}\right)n_x^2+n_z^2\right)}^2}}-1\right)$

其中，u是波场，t是时间，S是标量算子，v₀是所述地下地层中沿各向异性的对称轴的速度，n_x和n_y是水平相位方向，n_z是垂直相位方向，并且其中，ε和δ是各向异性参数。

30.根据任何权利要求14所述的方法，其中，步骤b包括执行FFT或有限差分数值法，以确定所述波动方程中的del项以及del的平方项。

31.一种处理数据的方法，包括：

(a)提供与地球的各向异性地下地层相关联的第一数据。

(b)提供所述各向异性地下地层的波场或波场的分量，所述波场或所述波场的分量是通过求解波动方程而数值确定的；以及

(c)使用所确定的波场对所述第一数据进行处理，以产生与所述地下地层相关联的第二数据。

32.一种方法，包括基于如前文在方程7或方程8中所定义的波动方程来确定地球的各向异性地下地层中的波场。

33.一种方法，包括基于伪微分波方程来确定地球的各向异性地下地层中的波场，但是所述伪微分波方程被布置成数值可解的形式。

34.一种使用波场的分量来对地球物理数据进行处理的方法，所述波场的分量是由估计根据任何前述权利要求所述的波动方程而获得的。

35.根据权利要求14到34中任一项所述的方法，其还包括：通过执行有限差分法或FFT数值法来计算所述波动方程的空间导数。

36.一种确定地球的各向异性地下地层中的波场的方法，其包括以下步骤：

计算第一空间导数分量，以获得第一传播时间处的第一波场的梯度；

计算具有基于地下地层各向异性模型和波场相位方向的分量的算子S，使用所述第一波场的所述梯度来提供所述相位方向；

将所述算子S与所述第一波场的所述梯度以及第二空间导数分量组合，以提供组合的结果；

使用所述组合的结果来计算第二传播时间处的第二波场。

37.一种通过对各向异性波动方程进行数值求解、使用FFT或FD法确定空间导数、并且使用积分法确定时间导数来确定地球的各向异性地下地层中的波场的方法。

38.一种用于执行根据任何前述权利要求所述的方法的设备。

39.一种计算机程序，包括用于执行根据权利要求1到36中任一项所述的方法的机器可读指令。

40.一种计算机装置，包括被布置为运行根据权利要求39所述的计算机程序以执行所述方法的处理器。