CN102105900A - 在各向异性介质中传播伪声准p波的方法 - Google Patents

Abstract

一种用于在具有可变倾斜的各向异性介质中保持稳定并且不限于弱各向异性条件的伪声准P波传播的计算机实现方法。该方法包括获取感兴趣的地下区域的地震勘探体积，并且确定所述地震勘探体积的模拟几何体系。该方法进一步包括将所述模拟几何体系用于初始条件来传播至少一个波场通过所述地震勘探体积，防止沿着所述地震勘探体积的对称轴累积能量，以及通过利用有限准S波速度来保证应力-应变关系中的正刚度系数，从而产生稳定波场。该方法包括利用所述稳定波场来生成感兴趣的地下区域的地下图像。

Description

在各向异性介质中传播伪声准P波的方法

技术领域

本发明一般涉及使用地震信号的地球物理勘探，尤其涉及在可变倾斜各向异性介质中传播伪声准P波并将传播的波场用于地下特性表征的方法。

背景技术

由于与应力有关的矿物和缺陷的择优排序，在许多油气勘探区域(例如，墨西哥湾、北海、和西非海岸)中处处可以观察到各向异性。在这些区域中，往往可以将岩石特性表征成对称轴垂直或倾斜的横向各向同性(“TI”)介质。各向异性介质中的波传播呈现与各向同性介质中的波传播不同的运动学和动力学，因此，需要各向异性模拟和偏移方法来适当地成像储层用于油气勘探。

然而，三维(“3D”)各向异性地震模拟和偏移是计算强度大的任务。与完全弹性的现有技术解决方案相比，基于分散关系的模拟和偏移是计算效率高的替代品。在一种现有技术的方法(Alkhalifah(2000))中，引入了垂直横向各向同性(“VTI”)介质的伪声近似。在那个现有技术方法的近似中，沿着垂直对称轴将剪切波的相速设置成零。如Grechka等人(2004)所述，这种简化不会消除其它方向上的剪切波。根据Alkhalifah的近似，已经为VTI介质中的地震模拟和偏移提出了若干时空域伪声偏微分方程(PDE)(Alkhalifah，2000；Zhou等人，2006；以及Du等人，2008)。PDE的这些系统在运动学上非常接近于涉及矢量场的完全弹性解决方案。

作为VTI介质的一种扩展，正如在与背斜结构和/或逆冲岩席相关联的区域中观察到的那样，可以倾斜TI介质的对称轴(“TTI”)。Zhou等人(2006)通过应用绕对称轴的旋转，将他们的VTI伪声方程推广到用于2D TTI介质的系统。因此，准SV波的相速在与倾斜轴平行或垂直的方向上是零。Lesage等人根据相同相速近似，进一步将Zhou的TTI系统从2D推广到3D。然而，由于倾斜和/或某些岩石特性的迅速横向变化(当垂直速度大于水平速度时)，这些现有技术的伪声模拟和偏移方法可能在数值上变得不稳定，并且导致不稳定的波传播。

本领域的普通技术人员应该懂得，各向同性介质中的平面波极化矢量与缓慢矢量平行(P波)或正交(S波)。除了特定传播方向以外，在各向异性介质中没有纯粹的纵波和剪切波。由于这种原因，在各向异性波理论中，快模往往被称为“准P”波，而慢模往往被称为“准S₁”和“准S₂”。

发明内容

本发明为各向异性介质提供了伪声模拟方法和伪声偏移方法。本发明实施例的一些方面包括一种用于在可变倾斜各向异性介质中保持稳定并且不限于弱各向异性条件的伪声准P波传播的计算机实现方法。该方法还包括获取感兴趣的地下区域的地震勘探体积，并且确定所述地震勘探体积的模拟几何体系。该方法进一步包括：将所述模拟几何体系用于初始条件来传播至少一个波场通过所述地震勘探体积，防止沿着对称轴累积能量，以及通过利用小的有限准S波速度来保证应力-应变关系中的正刚度系数，从而产生稳定波场。该方法还包括利用所述稳定波场来生成感兴趣的地下区域的地下图像。

本发明的另一个实施例包括一种地球物理地震偏移方法，包含下列步骤：建立与地震勘探体积相对应的地震数据集和速度/各向异性模型；以及对于每个共炮/接收器记录，将边界条件设置成包括来自源位置的激励。该实施例进一步包括按照伪声波方程或其等效物前向传播波场：

$\left\{\begin{array}{c}\frac{{\∂}^2}{{\∂t}^2}P=v_{p_0}^2[((1+2\mathrm{\ϵ})f_2+f_1)-(f-1)f_3]P+v_{p_0}^4[2f(\mathrm{\δ}-\mathrm{\ϵ})f_1{f}_2+(f-1)((1+2\mathrm{\ϵ})f_2+f_1)f_3]Q\\ \frac{{\∂}^2}{{\∂t}^2}Q=P\end{array}\right.---\left[1\right]$

其中，

$\left\{\begin{array}{c}f=1-{\left(\frac{v_{s_0}}{v_{p_0}}\right)}^2\\ f_1={\sin }^2{\mathrm{\θ}}_0({\cos }^2{\mathrm{\φ}}_0\frac{{\∂}^2}{{\∂x}^2}+{\sin }^2{\mathrm{\φ}}_0\frac{{\∂}^2}{{\∂y}^2}+\sin 2{\mathrm{\φ}}_0\frac{\∂}{\∂x}\frac{\∂}{\∂y})+{\cos }^2{\mathrm{\θ}}_0\frac{{\∂}^2}{{\∂z}^2}+\sin 2{\mathrm{\θ}}_0({\cos \mathrm{\φ}}_0\frac{\∂}{\∂x}\frac{\∂}{\∂z}+\sin {\mathrm{\φ}}_0\frac{\∂}{\∂y}\frac{\∂}{\∂z})\\ f_2=(1-{\sin }^2{\mathrm{\θ}}_0{\cos }^2{\mathrm{\φ}}_0)\frac{{\∂}^2}{{\∂x}^2}+(1-{\sin }^2{\mathrm{\θ}}_0{\sin }^2{\mathrm{\φ}}_0)\frac{{\∂}^2}{{\∂y}^2}+s{\mathrm{in}}^2{\mathrm{\θ}}_0\frac{{\∂}^2}{\∂z^2}-\mathrm{si}{n}^2{\mathrm{\θ}}_0{\sin 2\mathrm{\θ}}_0\frac{\∂}{\∂x}\frac{\∂}{\∂y}-\sin 2{\mathrm{\θ}}_0({\cos \mathrm{\φ}}_0\frac{\∂}{\∂x}\frac{\∂}{\∂z}+\sin {\mathrm{\φ}}_0\frac{\∂}{\∂y}\frac{\∂}{\∂z})\\ f_3=f_1+f_2=\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2}+\frac{{\∂}^2}{{\∂z}^2}\end{array}\right.$

Vs₀是准SV波的垂直速度，Vp₀是准P波的垂直速度，θ_o是TI介质中对称轴相对于垂直方向的倾斜角，φ_o是对称轴的方位角，ε、δ是汤姆森各向异性参数，P是标量波场，以及Q是辅助函数。该实施例还包括，对于每个共炮/接收器记录，将边界条件设置成反向传播记录的炮记录，并且按照上面的伪声波方程来反向传播地震数据。该实施例包括应用诸如计算的前向波场和反向波场或它们的等效格林函数之间的互相关(但不限于此)的成像条件来导出地下图像。

本发明的另外实施例还包括如下步骤：通过逆时偏移(RTM)、高斯束偏移、Kirchhoff偏移或其它基于波方程的偏移来传播波场或计算格林函数。

本发明的另外实施例还包括如下步骤：作为选项，除了互相关之外，还应用涉及照明标准化和/或反射角域道集生成和/或相位-振幅补偿的成像条件。

本发明的另外实施例还包括如下步骤：处理共炮/接收器信号，并且在其它相关域中以及以其它模拟和偏移形式来传播波场，所述其它相关域包括(但不限于)共偏置、共方位角和共反射角，所述其它模拟和偏移形式包括(但不限于)延迟炮、平面波和相位编码。

本发明的另外实施例还包括如下步骤：取代汤姆森参数，使用诸如正常时差速度、水平速度的其它等效术语来传播波场或计算格林函数。

本发明的另外实施例包括一种地球物理地震偏移方法，包含如下步骤：建立与地震勘探体积相对应的地震数据集和速度/各向异性模型；以及对于每个共炮/接收器记录，将边界条件设置成包括来自源位置的激励。该实施例还包括按照下列伪声波方程或其等效物来前向传播波场：

$\left\{\begin{array}{c}\frac{{\∂}^2}{{\∂t}^2}P=v_{p_0}^2[((1+2\mathrm{\ϵ})f_2+f_1)-(f-1)f_3]P+v_{p_0}^4[2f(\mathrm{\δ}-\mathrm{\ϵ})f_1{f}_2+(f-1)((1+2\mathrm{\ϵ})f_2+f_1)f_3]Q\\ \frac{{\∂}^2}{{\∂t}^2}Q=v_{p_0}^{2}P\end{array}\right.---\left[2\right]$

该实施例进一步包括：对于每个共炮/接收器记录，将边界条件设置成反向传播记录的炮记录，并且按照上面的伪声波方程反向传播地震数据。该实施例包括应用诸如(但不限于)计算的前向波场和反向波场或它们的等效格林函数之间的互相关的成像条件来导出地下图像。

在地球物理地震偏移中，本发明的不同实施例可以利用其它伪声波方程来前向传播波场。例如，本发明的一个实施例包括按照下面的伪声波方程或其等效物来前向传播波场：

$\left\{\begin{array}{c}{\mathrm{\ω}}^{2}P=v_{p_0}^2[((1+2\mathrm{\ϵ})f_2+f_1)-(f-1)f_3]P+v_{p_0}^4[2f(\mathrm{\δ}-\mathrm{\ϵ})f_1{f}_2+(f-1)((1+2\mathrm{\ϵ})f_2+f_1)f_3]Q\\ {\mathrm{\ω}}^{2}Q=P\end{array}\right.---\left[3\right]$

其中，ω是角频率。

用于地球物理地震偏移的本发明的进一步实施例包括按照下面的伪声波方程或其等效物来前向传播波场：

$\left\{\begin{array}{c}\frac{{\∂}^2}{{\∂t}^2}P=v_{p_0}^2[((1+2\mathrm{\ϵ})f_2+f_1)-(f-1)f_3]P+v_{p_0}^{4}f[(1+2\mathrm{\δ})-(1-2\mathrm{\ϵ})]f_{1}Q+v_{p_{.0}}^4(f-1)[(1+2\mathrm{\ϵ})f_2+f_1]R\\ \frac{{\∂}^2}{{\∂t}^2}Q=f_{2}P\\ \frac{{\∂}^2}{{\∂t}^2}R=f_{3}P\end{array}\right.---\left[4\right]$

其中，Q和R是辅助函数。

本发明的另一个实施例包括一种地球物理地震偏移方法，包含如下步骤：建立与地震勘探体积相对应的地震数据集和速度/各向异性模型；以及对于每个共炮/接收器记录，将边界条件设置成包括来自源位置的激励。该实施例还包括按照用于倾斜介质的伪声波方程及其等效公式来前向传播波场：

$\left\{\begin{array}{c}\frac{\∂}{\∂t}P=(1+2\mathrm{\η})v_{\mathrm{nmo}}^2(\frac{\∂}{\∂x}U+\frac{\∂}{\∂y}V)+\frac{v_{p_0}^2}{1+2\mathrm{\η}}\frac{\∂}{\∂z}R\\ \frac{\∂}{\∂t}U=\frac{\∂}{\∂x}P\\ \frac{\∂}{\∂t}V=\frac{\∂}{\∂y}P\\ \frac{\∂}{\∂t}Q=\frac{v_{p_0}^2}{1+2\mathrm{\η}}\frac{\∂}{\∂z}R\\ \frac{\∂}{\∂t}R=\frac{\∂}{\∂z}P+2\mathrm{\η}\frac{\∂}{\∂z}Q\end{array}\right.---\left[5\right]$

是准P波的垂直速度，

是准P波的正常时差速度，η＝(ε-δ)/(1+2δ)是Alkhalifah-Tsvankin各向异性参数(按照汤姆森各向异性参数ε和δ表达)，P是标量波场，以及U、V、Q和R是辅助函数。该实施例进一步包括：对于每个共炮/接收器记录，将边界条件设置成反向传播记录的炮记录，并且按照上面的伪声波方程来反向传播地震数据。该实施例包括应用诸如计算的前向波场和反向波场或它们的等效格林函数之间的互相关的成像条件来导出地下图像。

用于地球物理地震偏移的本发明的不同实施例可以利用其它伪声波方程为倾斜介质前向传播波场。例如，本发明的一个实施例包括按照用于倾斜介质的伪声波方程及其等效公式来前向传播波场：

$\left\{\begin{array}{c}\frac{{\∂}^2}{{\∂t}^2}P=[(1+2\mathrm{\η})v_{\mathrm{nmo}}^2+{\mathrm{av}}_{p_0}^2]g_{2}P+(1+a)v_{p_0}^2{g}_{1}P-a(1+2\mathrm{\η})v_{\mathrm{nmo}}^2{g}_{2}Q-[(2\mathrm{\η}+a)v_{\mathrm{nmo}}^2+{\mathrm{av}}_{p_0}^2]g_{1}Q-{\mathrm{av}}_{p_0}^2{g}_{1}R\\ \frac{{\∂}^2}{{\∂t}^2}Q=v_{p_0}^2{g}_{2}P\\ \frac{{\∂}^2}{{\∂t}^2}R=v_{p_0}^2{g}_{1}P\end{array}\right.---\left[6\right]$

其中，

$\left\{\begin{array}{c}a=1-f={(v_{s_0}/v_{p_0})}^2\\ g_1=\frac{{\∂}^2}{{\∂z}^2}\\ g_2=\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2}\end{array}\right.$

是准P波的垂直速度，

是准P波的正常时差速度，a是剪切波与P波速度比的平方，η＝(ε-δ)/(1+2δ)是Alkhalifah-Tsvankin各向异性参数(按照汤姆森各向异性参数ε和δ表达)，以及Q和R是辅助函数。

用于地球物理地震偏移的本发明的进一步实施例包括：按照用于倾斜介质的伪声波方程及其等效公式或它们的导出公式/等效物来前向传播波场：

$\frac{{\∂}^4}{{\∂t}^4}F-\left\{\right[(1+2\mathrm{\η})v_{\mathrm{nmo}}^2+{\mathrm{av}}_{p_0}^2](\frac{{\∂}^4}{{\∂x}^2{\∂\mathrm{rt}}^2}F+\frac{{\∂}^4}{{\∂y}^2{\∂t}^2}F)+(1+a)V_{p_0}^2\left(\frac{{\∂}^4}{{\∂z}^2{\∂t}^2}F\right)\}$

$+a(1+2\mathrm{\η})v_{\mathrm{nmo}}^2+a{v}_{p_0}^2(\frac{{\∂}^4}{{\∂x}^4}F+\frac{{\∂}^4}{{\∂x}^2{\∂y}^2}F+\frac{{\∂}^4}{{\∂y}^4}F)---\left[7\right]$

$+[(2\mathrm{\η}+a)v_{\mathrm{nmo}}^2+v_{p_0}^2+{\mathrm{av}}_{p_0}^2](\frac{{\∂}^4}{{\∂x}^2{\∂z}^2}F+\frac{{\∂}^4}{{\∂y}^2{\∂z}^2}F)+{\mathrm{av}}_{p_0}^4\left(\frac{{\∂}^4}{{\∂z}^2}F\right)=0$

其中，F是标量波场。

本发明的另一个实施例包括一种地球物理地震偏移方法，包含如下步骤：建立与地震勘探体积相对应的地震数据集和速度/各向异性模型；以及对于每个共炮/接收器记录，将边界条件设置成包括来自源位置的激励。该实施例还包括按照伪声波方程或其导出公式/等效物来前向传播波场：

$\frac{{\∂}^{4}P}{{\∂t}^4}-v_{p_0}^2[(1+2\mathrm{\ϵ})f_2+f_1-(f-1)f_3]\frac{{\∂}^{2}P}{{\∂t}^2}$ [8]

$+v_{p_0}^4[2f(\mathrm{\ϵ}-\mathrm{\δ})f_1{f}_2-(f-1)((1+2\mathrm{\ϵ})f_2+f_1)f_3]P=0$

该实施例进一步包括：对于每个共炮/接收器记录，将边界条件设置成反向传播记录的炮记录，并且按照上面的伪声波方程反向传播地震数据。该实施例包括应用诸如(但不限于)计算的前向波场和反向波场或它们的等效格林函数之间的互相关的成像条件来导出地下图像。

本发明的另一个实施例包括一种地球物理地震模拟方法，包含如下步骤：建立与地震勘探体积相对应的速度/各向异性模型；以及对于每个炮，设置波场的初始条件。该实施例还包括按照伪声波方程或其等效物来前向传播波场：

$\left\{\begin{array}{c}\frac{{\∂}^2}{{\∂t}^2}P=v_{p_0}^2[((1+2\mathrm{\ϵ})f_2+f_1)-(f-1)f_3]P+v_{p_0}^4[2f(\mathrm{\δ}-\mathrm{\ϵ})f_1{f}_2+(f-1)((1+2\mathrm{\ϵ})f_2+f_1)f_3]Q+\mathrm{\δ}(\stackrel{\→}{x}-{\stackrel{\→}{x}}_s)w\left(t\right)\\ \frac{{\∂}^2}{{\∂t}^2}Q=P+\mathrm{\δ}(\stackrel{\→}{x}-{\stackrel{\→}{x}}_s)w\left(t\right)\end{array}\right.---\left[9\right]$

其中，w(t)是源函数，以及是源位置的矢量。可以改变源项及其插入形式而不影响支配PDE。

用于地球物理地震模拟的本发明的另一个实施例包括：按照用于倾斜介质的伪声波方程(方程5)及其等效公式来前向传播波场。

用于地球物理地震模拟的本发明的另一个实施例包括：按照用于倾斜介质的伪声波方程(方程6)及其等效公式来前向传播波场。

还应当理解，本发明打算与通常包括电子配置的系统一起使用，所述电子配置包括至少一个处理器、用于存储程序代码或其它数据的至少一个存储设备、可选的视频监视器或其它显示设备(即，液晶显示器)以及至少一个输入设备。所述处理器最好是能够显示图像和处理复杂数学算法的基于微处理器或微控制器的平台。所述存储设备可以包括用于存储在与本发明相关联的特定处理期间生成或使用的事件或其它数据的随机访问存储器(RAM)。所述存储设备还可以包括用于存储用于本发明的控制和处理的程序代码的只读存储器(ROM)。

一个这样的实施例包括一种被配置成进行在可变倾斜各向异性介质中保持稳定并且不限于弱各向异性条件的伪声准P波传播的系统。所述系统包括：数据存储设备，其具有包括感兴趣的地下区域的地震勘探体积的计算机可读数据；以及处理器，其被配置和安排成执行存储在处理器可访问存储器中的机器可执行指令以便执行一种方法。这个特定实施例的方法包括：确定所述地震勘探体积的模拟几何体系；以及将所述模拟几何体系用于初始条件来传播至少一个波场通过所述地震勘探体积，防止沿着所述地震勘探体积内的各向异性区域的对称轴累积能量，以及保证应力-应变关系中的正刚度系数，从而产生稳定波场。该方法进一步包括：利用所述稳定波场来生成感兴趣的地下区域的地下图像。

通过参照形成本说明书的一部分、相同标号在各个图形中表示相应部件的附图对如下描述和所附权利要求书加以研究，本发明的这些和其它目的、特征、和特性，以及操作方法、结构的相关元件的功能、部件的组合、和制造成本将变得更加显而易见。然而，不言而喻，这些附图只用于例示和描述的目的，而无意作为限制本发明的定义。正如用在说明书和权利要求书中的那样，除非上下文另有明确指明，单数形式“一个”、“一种”、和“该”也包括复数指示物。

附图说明

图1是例示依照本发明一个或多个实施例的方法的流程图；

图2是例示依照本发明一个或多个实施例的方法的流程图；

图3是例示依照本发明一个或多个实施例的方法的流程图；

图4例示了按照现有技术的Alkhalifah近似的示范性波传播模拟，其中f＝1；

图5例示了按照本发明一个实施例的示范性波传播模拟；

图6例示了按照本发明一个实施例的示范性波传播模拟，其中V_so/V_po=0.01；

图7例示了按照现有技术的Alkhalifah近似的示范性波传播模拟，其中V_so＝0；

图8例示了按照现有技术的Alkhalifah近似的示范性相速分布；

图9例示了按照现有技术的Alkhalifah近似的示范性群速分布；

图10例示了本发明一个实施例的示范性相速分布，其中V_so/V_po＝0.01；

图11例示了本发明一个实施例的示范性群速分布，其中V_so/V_po＝0.01；

图12例示了按照利用一阶5×5PDF系统的本发明一个实施例，具有可变倾斜对称轴的介质中的示范性波传播模拟；

图13例示了用在本发明一个实施例中的几何体系的示意图；以及

图14是执行依照本发明实施例的方法的系统的实施例的示意性例示图。

具体实施方式

本发明的一个实施例例示在图1中，其中，流程图10描述了传播在具有可变倾斜的各向异性介质中保持稳定的准P波的方法。本发明不限于弱各向异性条件。这个特定实施例包括获取感兴趣的地下区域的地震勘探体积12，并且确定地震勘探体积的模拟几何体系14。该实施例进一步包括将模拟几何体系用于初始条件来传播至少一个波场通过所述地震勘探体积，防止沿着所述地震勘探体积的对称轴累积能量，以及利用有限准S波速度来保证应力-应变关系中的正刚度系数，从而产生稳定波场16。然后可以利用所述稳定波场来生成感兴趣的地下区域的地下图像18。

本领域的普通技术人员应该懂得，本发明的不同实施例可以为各向异性介质提供伪声模拟方法或伪声偏移方法。例如，图2例示了用于具有可变倾斜的各向异性介质中的波传播的伪声模拟方法的一个实施例的流程图，其中该方法不局限于弱各向异性条件。该实施例包括获取感兴趣的地下区域的地震勘探体积22，并且确定所述地震勘探体积的模拟几何体系24。该实施例还包括将模拟几何体系用于初始条件来传播至少一个波场通过所述地震勘探体积，其中，沿着所述地震勘探体积的对称轴的人为准剪切波速度大于或等于零，从而防止沿着对称轴累积能量，从而产生稳定波场26。然后可以利用稳定波场来生成感兴趣的地下区域的地下图像28。

图3例示了可以用于伪声偏移的本发明另一个实施例的流程图。该实施例包括获取感兴趣的地下区域的地震勘探体积32，并且确定地震勘探体积的模拟几何体系34。该实施例还包括将模拟几何体系用于初始条件来传播至少一个波场通过地震勘探体积，其中，准剪切波能量不沿着地震勘探体积的对称轴累积，从而产生稳定波场36。然后可以利用稳定波场来生成感兴趣的地下区域的地下图像38。

本发明相对于传统的声各向异性模拟和偏移具有若干优点。本发明在具有可变倾斜的TI介质中提供了稳定的波传播方式，因此可以获得仿真波场传播和反射率图像。现有技术的伪声模拟和偏移方法基于Alkhalifah近似，在所述Alkhalifah近似中，沿着对称轴将剪切波的相速设置成零。尽管现有技术的方法可以用在恒定倾斜的TI介质中，但在倾斜变化可以在对称轴附近局部集中高能量的区域中，零速度剪切波可以使波传播不稳定(即，振幅变成无界)。图4示出了在可变倾斜介质中(例如，在背斜结构的顶部附近)现有技术的方法(f＝1)是不稳定的40。相反，图5示出了基于本发明的波传播(f＝0.98)在相同介质中保持稳定42。另外，本发明可以提供控制剪切波与P波速度比以优化模拟和偏移的结果的灵活性。例如，可以将剪切波与P波速度比设置成接近实际值，以近似弹性波传播中的运动学。而且，在某些岩石中，相对于对称轴，垂直速度可以大于水平速度。在这样的情况下，基于Alkhalifah近似的波方程将导致负刚度矩阵，从而与数值实现算法无关地产生不稳定波场。本发明可以使用有限剪切波速度来保证应力-应变关系中的正刚度系数，从而生成稳定波传播。

在基于Alkhalifah近似的现有技术方法中，使剪切波相速等于零并不会消除剪切波。而是高能量集中在对称轴附近。唯一例外是剪切波处处都消失的椭圆各向异性(即，ε＝δ)。在本发明中，将垂直剪切波速度从等于零中解脱出来，因此，能量较少集中在对称轴附近。由于即使满足椭圆各向异性条件也存在附加交叉导数，所以剪切波将不会消失。

就计算成本而言，本发明的实施例所使用的PDE与现有技术的方法相比涉及要计算附加空间导数项。在具有可变倾斜的区域中，需要与非零Vs₀相关联的附加工作负荷以实现地震模拟和偏移所需的稳定性和可靠性。在倾斜几乎恒定或非常平缓的区域中，可以跳过所述附加工作负荷。

对于本领域的普通技术人员来说清楚的是，可以不偏离本发明的范围地以许多方式改变上面的实施例。例如，本领域的普通技术人员可以明显看出，在本发明中，可以将不同初始条件或边界条件或PDE的不同线性组合方便地用在模拟和偏移中。

在本发明的一个实施例中，各向异性模拟方法包括：建立与地震勘探体积相对应的速度和各向异性模型；设置诸如源激励的初始条件；按照方程[1]或其等效物，在具有倾斜对称轴或垂直对称轴的横向各向同性介质中传播波。对于正向模拟，在方程[1]或方程[2]中，需要在方程的右侧引入形式的源函数，其中

是源位置，以及w(t)是源小波。

在本发明的上述实施例中，与现有技术方法中将f四舍五入成1的近似相反，方程[3]中的垂直剪切波速度可以非零(因此，f可以不同于1)。因此，在本发明中，剪切波的相速在与对称轴平行和垂直的两个方向上都可以非零。在具有可变倾斜的介质中，准剪切波的有限速度可以避免往往发生在对称轴附近的高能量的局部集中。本发明不需要弱各向异性假设。

利用上面的PDE，可以为各向异性介质导出本发明的其它实施例。如果倾斜角θ_o＝0，那么上面的PDE简化成3D VTI系统，类似地，当φ_o＝0时，简化成3D HTI系统。用于3D VTI介质的二阶3×3系统可以采取方程6中的形式。通过用在方程[1]中给出的f₁和f₂取代在方程[6]中给出的g₁和g₂，可将PDE的这种系统推广到用于倾斜TI介质的其等效公式。

作为f＝1时的方程[1]或[2]或[4]中的PDE的替代物，方程[5]中的PDE的一阶5×5系统是双曲线的，并且在具有可变倾斜的TI介质中是稳定的。本发明的这个实施例是对称双曲线的(适定的，甚至对于可变系数)。这个系统也可推广到可变倾斜TTI。3D中的PDE的上述完整一阶5×5系统在2D中减小成4×4。

如上所述，本发明的另外实施例还提供了伪声偏移方法。一个实施例包括如下步骤：建立与地震勘探体积相对应的地震数据集和速度/各向异性模型；设置波传播的边界条件；按照方程[1]方程[2]、方程[4]或方程[6]或它们的等效物，在各向异性介质中分开传播来自源激励的波和记录的地震数据；以及应用诸如(但不限于)两个传播波场之间的互相关的成像条件来获取地下图像。可以应用不同的初始和/或边界条件，而不影响本发明的范围。用于传播源小波的示范性边界条件(例如，根据方程[1])如下：

$\left\{\begin{array}{c}P(x,y,z=0;t)=\mathrm{\δ}(\stackrel{\→}{x}-{\stackrel{\→}{x}}_s){\∫}_0^{t}w\left(t^{\′}\right){\mathrm{dt}}^{\′}\\ Q(x,y,z=0;t)=\mathrm{\δ}(\stackrel{\→}{x}-{\stackrel{\→}{x}}_s){\∫}_0^{t}w\left(t^{\′}\right){\mathrm{dt}}^{\′}\end{array}\right.---\left[10\right]$

以及用于地震数据的逆时外推的边界条件如下：

$\left\{\begin{array}{c}P(x,y,z=0;t)=D(x,y,x_s,y_s;t)\\ Q(x,y,z=0;t)=D(x,y,x_s,y_s;t)\end{array}\right.---\left[11\right]$

其中，w(t’)是源函数，x_s是源位置，以及D(x，y，x_s，y_s；t)是要偏移的炮记录。

下列例子例示了本发明的进一步实施例：

1.为VTI介质中的准P波建立四阶分散关系

v_so未设置成零的用于VTI介质的Tsvankin相速关系导致分散关系：

ω⁴-Bω²+C＝0                                             [12]

其中：

$\left\{\begin{array}{c}B=[(1+2\mathrm{\η})v_{\mathrm{nmo}}^2+{\mathrm{av}}_{p_0}^2]k_h^2+(1+a)v_{p_0}^2{k}_z^2\\ C=a(1+2\mathrm{\η})v_{\mathrm{nmo}}^2{v}_{p_0}^2{k}_h^4+[(2\mathrm{\η}+a)v_{\mathrm{nmo}}^2{v}_{p_0}^2+{\mathrm{av}}_{p_0}^4]k_h^2{k}_z^2+{\mathrm{av}}_{p_0}^4{k}_z^4\end{array}\right.---\left[13\right]$

ω是角频率，k_z是垂直波数，以及

是水平波数矢量(k_x，k_y)的幅度(的平方)。方程[12]允许两对解：

$\left\{\begin{array}{c}{\mathrm{\ω}}_{\mathrm{qP}\±}=\±\sqrt{\frac{B+\sqrt{B^2-4c}}{2}}\\ {\mathrm{\ω}}_{\mathrm{qS}\±}=\±\sqrt{\frac{B-\sqrt{B^2-4C}}{2}}\end{array}\right.---\left[14\right]$

ω_qP±对应于准P波；ω_qS±对应于准SV波。

2.为VTI介质中的准P波建立四阶PDE

在傅立叶域中将方程[12]应用于波场

并且采取逆傅立叶变换(F(x，y，z，t))提供：

$\frac{{\∂}^4}{{\∂t}^4}F-\left\{\right[(1+2\mathrm{\η})v_{\mathrm{nmo}}^2+{\mathrm{av}}_{p_0}^2](\frac{{\∂}^4}{{\∂x}^2{\∂t}^2}F+\frac{{\∂}^4}{{\∂y}^2{\∂t}^2}F)+(1+a)v_{p_0}^2\left(\frac{{\∂}^4}{{\∂z}^2{\∂t}^2}F\right)\}$

$+a(1+2\mathrm{\η})v_{\mathrm{nmo}}^2{v}_{p_0}^2(\frac{{\∂}^4}{{\∂x}^4}F+\frac{{\∂}^4}{{\∂x}^2{\∂y}^2}F+\frac{{\∂}^4}{{\∂y}^4}F)+[(2\mathrm{\η}+a)v_{\mathrm{nmo}}^2{v}_{p_0}^2---[15]$

$+a{v}_{p_0}^2](\frac{{\∂}^4}{{\∂x}^2{\∂z}^2}F+\frac{{\∂}^4}{{\∂y}^2{\∂z}^2}F)+a{v}_{p_0}^4\left(\frac{{\∂}^4}{{\∂z}^2}F\right)=0$

3.为VTI介质建立PDE的二阶3×3系统

设：

$P(x,y,z,t)=\frac{{\∂}^2}{{\∂t}^2}F(x,y,z,t)---\left[16\right]$

其中，F(x，y，z，t)是满足方程[15]的波场。假设初始条件为：

$F(x,y,z,t=0)\≡\frac{\∂}{\∂t}F(x,y,z,t=0)\≡0---\left[17\right]$

导致：

$F(x,y,z,t)={\∫}_0^t{\∫}_0^{t^{\′}}P(x,y,z,t^{\′\′}){\mathrm{dt}}^{\′\′}{\mathrm{dt}}^{\′}---\left[18\right]$

设：

$\left\{\begin{array}{c}Q(x,y,z,t)=v_{p_0}^2(\frac{{\∂}^2}{{\∂x}^2}F+\frac{{\∂}^2}{{\∂y}^2}F)={\∫}_0^t{\∫}_0^{t^{\′}}(\frac{{\∂}^2}{{\∂x}^2}P+\frac{{\∂}^2}{{\∂y}^2}P)(x,y,z,t^{\′\′}){\mathrm{dt}}^{\′\′}{\mathrm{dt}}^{\′}\\ R(x,y,z,t)=v_{p_0}^2\left(\frac{{\∂}^2}{{\∂z}^2}F\right)={\∫}_0^t{\∫}_0^{t^{\′}}\left(\frac{{\∂}^2}{{\∂z}^2}P\right)(x,y,z,t^{\′\′}){\mathrm{dt}}^{\′\′}{\mathrm{dt}}^{\′}\end{array}\right.---\left[19\right]$

然后，通过方程[4]让方程[15]等效于PDE的二阶3×3系统。

图6示出了将上面的PDE用于本发明的那个特定实施例的VTI介质中的波前传播。与基于现有技术方法(例示在图7中)的波前相比，外部qP波前(图6中的44和图7中的48)保持几乎相同，但内部qSV波前(图6中的46和图7中的50)具有不同于菱形的形式。图8和图9分别示出了按照Alkhalifah(现有技术)近似的相速和群速。相反，图10和图11分别示出了按照本发明实施例的相速和群速。与Alkhalifah近似相比，qSV波的相速从沿着对称轴等于零中解脱出来。因此，群速的最大值或高能量不像现有技术的方法中那样沿着对称轴集中。通过应用绕倾斜轴的旋转，可对恒定倾斜的TTI介质应用相同观察。

4.为VTI介质建立PDE的一阶5×5系统

Du等人(2008)为v_s0＝0的VTI介质提出了下列PDE的一阶2×2系统：

$\left\{\begin{array}{c}\frac{\∂}{\∂t}p=(1+2\mathrm{\η})v_{\mathrm{nmo}}^2{g}_{2}p+v_{p_0}^2{g}_{1}q\\ \frac{\∂}{\∂t}q=v_{\mathrm{nmo}}^2{g}_{2}p+v_{p_0}^2{g}_{1}q\end{array}\right.---\left[20\right]$

其中，g₁和g₂在方程[6]之后给出，p是标量波场，以及q是辅助函数。新波场P以及新辅助函数U、V、Q和R通过下式定义：

$P=\frac{\∂}{\∂t}p,$ $U=\frac{\∂}{\∂x}p,$ $V=\frac{\∂}{\∂y}p,$ $Q=\frac{\∂}{\∂t}\left(\frac{(1+2\mathrm{\η})q-p}{2\mathrm{\η}}\right),$ $R=(1+2\mathrm{\η})\frac{\∂}{\∂z}q.---\left[21\right]$

然后方程5是PDE的完整一阶5×5系统。这个系统可以通过对称化而显示成双曲线的。设：

$\tilde{P}=P,$ $\tilde{U}=v_{\mathrm{nmo}}\sqrt{1+2\mathrm{\η}}U,$ $\tilde{V}=v_{\mathrm{nmo}}\sqrt{1+2\mathrm{\η}}V,$ $\tilde{Q}=\sqrt{2\mathrm{\η}}Q,$ $\tilde{R}=\frac{V_{p_0}}{\sqrt{1+2\mathrm{\η}}}R.---\left[22\right]$

然后：

$\frac{\∂}{\∂t}\left[\begin{array}{c}\tilde{P}\\ \tilde{U}\\ \tilde{V}\\ \tilde{Q}\\ \tilde{R}\end{array}\right]{=M}_x\frac{\∂}{\∂x}\left[\begin{array}{c}\tilde{P}\\ \tilde{U}\\ \tilde{V}\\ \tilde{Q}\\ \tilde{R}\end{array}\right]+M_y\frac{\∂}{\∂y}\left[\begin{array}{c}\tilde{P}\\ \tilde{U}\\ \tilde{V}\\ \tilde{Q}\\ \tilde{R}\end{array}\right]M_z\frac{\∂}{\∂z}\left[\begin{array}{c}\tilde{P}\\ \tilde{U}\\ \tilde{V}\\ \tilde{Q}\\ \tilde{R}\end{array}\right]---\left[23\right]$

其中，

$M_x=\left[\begin{array}{ccccc}0& v_{\mathrm{nmo}}\sqrt{1+2\mathrm{\η}}& 0& 0& 0\\ v_{\mathrm{nmo}}\sqrt{1+2\mathrm{\η}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]---\left[24\right]$

$M_y=\left[\begin{array}{ccccc}0& 0& v_{\mathrm{nmo}}\sqrt{1+2\mathrm{\η}}& 0& 0\\ 0& 0& 0& 0& 0\\ v_{\mathrm{nmo}}\sqrt{1+2\mathrm{\η}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]---\left[25\right]$

$M_y=\left[\begin{array}{ccccc}0& 0& \frac{v_{p_0}}{\sqrt{1+2\mathrm{\η}}}& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ \frac{v_{p_0}}{\sqrt{1+2\mathrm{\η}}}& 0& 0& 0& \frac{v_{p_0\sqrt{2\mathrm{\η}}}}{\sqrt{1+2\mathrm{\η}}}\\ 0& 0& 0& \frac{v_{p_0}\sqrt{2\mathrm{\η}}}{\sqrt{1+2\mathrm{\η}}}& 0\end{array}\right]---\left[26\right]$

在2D中，消除变量V，并删除第三个方程，得出一阶4×4系统。

图12示出了在可变倾斜介质中受这种一阶PDE支配的稳定波前传播。

5.在TTI介质中建立四阶分散关系

通过解脱沿着对称轴V_so＝0(或f＝1)的Alkhalifah近似，可以从Tsvankin相速关系(2001)中导出下列方程：

$\left\{\begin{array}{c}\frac{{\cos }^2\stackrel{\‾}{\mathrm{\θ}}}{v^2}=\frac{{[f/2+\mathrm{\ϵ}{\sin }^2\stackrel{\‾}{\mathrm{\θ}}]}^2-\{\frac{v^2}{{v^2}_{p_0}}-{[1+\mathrm{\ϵ}\sin }^2\stackrel{\‾}{\mathrm{\θ}}]+f/2{\}}^2}{{2\mathrm{fv}}^2(\mathrm{\ϵ}-\mathrm{\δ}){\sin }^2\stackrel{\‾}{\mathrm{\θ}}}\\ f=1-{\left(\frac{v_{s_0}}{v_{p_0}}\right)}^2\end{array}\right.$

其中，相速v具有两个幅度的根：一个用于准P波，另一个用于准SV波，

是波前法线和对称轴之间的角度，以及其它参数定义在方程[1]中。

按照显示在图13中的几何体系，波前法线

和对称轴

以及其间的角度采取下列形式：

$\left\{\begin{array}{c}\stackrel{\→}{n}=\sin \mathrm{\θ}\cos \mathrm{\φ}\stackrel{\→}{i}+\sin \mathrm{\θ}\sin \mathrm{\φ}\stackrel{\→}{j}+\cos \mathrm{\θ}\stackrel{\→}{k}\\ \stackrel{\→}{t}={\sin \mathrm{\θ}}_0\cos {\mathrm{\φ}}_0\stackrel{\→}{i}+\sin {\mathrm{\θ}}_0\sin {\mathrm{\φ}}_0\stackrel{\→}{j}+\cos {\mathrm{\θ}}_0\stackrel{\→}{k}\\ \cos \stackrel{\‾}{\mathrm{\θ}}=\frac{\stackrel{\→}{n}\·\stackrel{\→}{t}}{\left|n\right|\left|\right|t|}=\sin \mathrm{\θ}\cos \mathrm{\φ}\sin {\mathrm{\θ}}_0\cos {\mathrm{\φ}}_0+\sin \mathrm{\θ}\sin \mathrm{\φ}\sin {\mathrm{\θ}}_0{\sin \mathrm{\φ}}_0+\cos \mathrm{\θ}{\cos \mathrm{\θ}}_0\end{array}\right.$

其中，θ₀是TI介质中对称轴相对于垂直方向的倾斜角，以及φ_o是对称轴的方位角。认识到：

$\left\{\begin{array}{c}\sin \mathrm{\θ}\cos \mathrm{\φ}=k_{x}v/\mathrm{\ω}\\ \sin \mathrm{\θ}\sin \mathrm{\φ}=k_{y}v/\mathrm{\ω}\\ \cos \mathrm{\θ}=k_{z}v/\mathrm{\ω}\end{array}\right.$

可以导出下列四阶分散关系：

${\mathrm{\ω}}^4-v_{p_0}^2[(1+2\mathrm{\ϵ})f_2+f_1-(f-1)f_3]{\mathrm{\ω}}^2$

$+v_{p_0}^4[2f(\mathrm{\ϵ}-\mathrm{\δ})f_1{f}_2-(f-1)((1+2\mathrm{\ϵ})f_2+f_1)f_3]=0---\left[27\right]$

6.在TTI介质中建立四阶PDE

将上述四阶分散关系的两侧乘以标量波场P，并且将频率-波数算子转换到时空域，TTI/VTI介质的四阶PDE采取方程[8]的形式。

7.为TTI介质建立PDE的二阶2×2系统

上述TTI介质的四阶伪声PDE可以经由方程[3]，通过2×2时空域PDE系统来求解：

其中，

PDE的2×2系统在水平速度v_h和正常时差(NMO)速度vn方面也可以采取等效形式。

在作为特殊情况的2D介质中，通过下列简化偏导数算子，上面的PDE仍然保持有效：

$\left\{\begin{array}{c}{f}_1=({\sin }^2{\mathrm{\θ}}_0\frac{{\∂}^2}{{\∂x}^2}+{\cos }^2{\mathrm{\θ}}_0\frac{{\∂}^2}{{\∂z}^2}+\sin 2{\mathrm{\θ}}_0\frac{\∂}{\∂x}\frac{\∂}{\∂z})\\ f_2=({\cos }^2{\mathrm{\θ}}_0\frac{{\∂}^2}{{\∂x}^2}+\mathrm{si}{n}^2{\mathrm{\θ}}_0\frac{{\∂}^2}{{\∂z}^2}-\sin 2{\mathrm{\θ}}_0\frac{\∂}{\∂x}\frac{\∂}{\∂z})\end{array}\right.$

8.为TTI介质建立PDE的二阶3×3系统

作为替代物，也可以使用不同的线性组合，通过方程[4]中的3×3时空域PDE或其等效物来求解TTI介质的四阶伪声PDE。

本发明的实施例可以在诸如现场可编程门阵列(FPGA)、图形处理单元(GPU)、存储单元的协处理器加速架构或通用计算机上实现。本发明提供了编有执行本发明的方法的指令的装置和通用计算机和/或协处理器，以及编码执行本发明的方法的指令的计算机可读介质。在图14中示意性地例示了执行本发明实施例的系统。该系统52包括数据存储设备或存储器54。可以使存储的数据用于诸如可编程通用计算机的处理器56。处理器56可以包括诸如显示器58和图形用户界面60的接口组件。图形用户界面(GUI)可以用于显示数据和处理后的数据产物，以及允许用户在实现该方法的各个方面的选项当中作出选择。数据可以经由总线62直接从数据获取设备或从中间存储设备或处理设施(未示出)传送给系统52。

尽管为了例示的目的，根据当前认为最实用的优选实施例对本发明作了详细描述，但应该明白，这样的细节仅仅为了那个目的，本发明不局限于所公开的实施例，而是相反，打算涵盖在所附权利要求书的精神和范围之内的所有修改和等效安排。例如，尽管本文提到一台计算机，但这台计算机可以包括通用计算机、专用计算机、编程成执行这些方法的ASIC(专用集成电路)、计算机阵列或网络、或其它适用计算设备。作为进一步的例子，还应该明白，本发明设想，可以尽可能地将任何实施例的一个或多个特征与任何其它实施例的一个或多个特征组合。

Claims (15)

1.一种用于在可变倾斜各向异性介质中保持稳定并且不限于弱各向异性条件的伪声准P波传播的计算机实现方法，所述方法包括：

获取感兴趣的地下区域的地震勘探体积；

确定所述地震勘探体积的模拟几何体系；

将所述模拟几何体系用于初始条件来传播至少一个波场通过所述地震勘探体积，防止沿着所述地震勘探体积内的各向异性区域的对称轴累积能量，以及利用有限准S波速度来保证应力-应变关系中的正刚度系数，从而产生稳定波场；以及

利用所述稳定波场来生成感兴趣的地下区域的地下图像。

2.如权利要求1所述的方法，其中，传播至少一个波场通过所述地震勘探体积包括：人为准剪切波速度沿着所述地震勘探体积内的各向异性区域的对称轴大于或等于零。

3.如权利要求1所述的方法，其中，传播至少一个波场通过所述地震勘探体积包括：限制准剪切波沿着所述地震勘探体积的所述对称轴累积。

4.如权利要求1所述的方法，其中，多个波场被传播通过所述地震勘探体积。

5.如权利要求1所述的方法，其中，传播至少一个波场通过所述地震勘探体积包括：利用逆时偏移、基于波方程的偏移、高斯束偏移或Kirchhoff偏移之一。

6.如权利要求1所述的方法，包括：前向和反向传播波场通过所述地震勘探体积，以及将成像条件应用于前向和反向波场或等效格林函数以导出地下图像。

7.如权利要求6所述的方法，其中，应用成像条件的步骤包括：应用前向和反向波场或等效格林函数之间的互相关来导出地下图像。

8.如权利要求7所述的方法，其中，应用成像条件的步骤包括：照明标准化、反射角域道集生成和相位-振幅补偿中的至少一个。

9.如权利要求1所述的方法，其中，成像输出几何体系包括共偏置、共方位角和共反射角域。

10.如权利要求1所述的方法，其中，传播至少一个波场包括：利用延迟炮、平面波和相位编码中的至少一个。

11.如权利要求1所述的方法，其中，传播至少一个波场包括：利用正常时差速度、水平速度和汤姆森参数中的至少一个。

12.一种配置成执行在可变倾斜各向异性介质中保持稳定并且不限于弱各向异性条件的伪声准P波传播的系统，所述系统包含：

数据存储设备，具有包括感兴趣的地下区域的地震勘探体积的计算机可读数据；以及

处理器，被配置和安排成执行存储在处理器可访问的存储器中的机器可执行指令以便执行包含如下步骤的方法：

确定所述地震勘探体积的模拟几何体系；

将所述模拟几何体系用于初始条件来传播至少一个波场通过所述地震勘探体积，防止沿着所述地震勘探体积内的各向异性区域的对称轴累积能量，以及利用有限准S波速度来保证应力-应变关系中的正刚度系数，从而产生稳定波场；以及

利用所述稳定波场来生成感兴趣的地下区域的地下图像。

13.如权利要求12所述的系统，其中，传播至少一个波场通过所述地震勘探体积包括：人为准剪切波速度沿着所述地震勘探体积内的各向异性区域的对称轴大于或等于零。

14.如权利要求12所述的系统，其中，传播至少一个波场通过所述地震勘探体积包括：限制准剪切波沿着所述地震勘探体积的所述对称轴累积。

15.如权利要求1所述的系统，其中，多个波场被传播通过所述地震勘探体积。

Images