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

具体实施方式

本发明的一个实施例例示在图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)=\delta (\overrightarrow{x}-{\overrightarrow{x}}_s){\int }_0^{t}w\left(t^{\prime }\right){\mathrm{dt}}^{\prime }\\ Q(x,y,z=0;t)=\delta (\overrightarrow{x}-{\overrightarrow{x}}_s){\int }_0^{t}w\left(t^{\prime }\right){\mathrm{dt}}^{\prime }\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\eta )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\eta )v_{\mathrm{nmo}}^2{v}_{p_0}^2{k}_h^4+[(2\eta +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}{\omega }_{\mathrm{qP}\pm }=\pm \sqrt{\frac{B+\sqrt{B^2-4c}}{2}}\\ {\omega }_{\mathrm{qS}\pm }=\pm \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{{\partial }^4}{{\partial t}^4}F-\left\{\right[(1+2\eta )v_{\mathrm{nmo}}^2+{\mathrm{av}}_{p_0}^2](\frac{{\partial }^4}{{\partial x}^2{\partial t}^2}F+\frac{{\partial }^4}{{\partial y}^2{\partial t}^2}F)+(1+a)v_{p_0}^2\left(\frac{{\partial }^4}{{\partial z}^2{\partial t}^2}F\right)\}$

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

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

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

设：

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

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

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

导致：

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

设：

$\left(\begin{array}{c}Q(x,y,z,t)=v_{p_0}^2(\frac{{\partial }^2}{{\partial x}^2}F+\frac{{\partial }^2}{{\partial y}^2}F)={\int }_0^t{\int }_0^{t^{\prime }}(\frac{{\partial }^2}{{\partial x}^2}P+\frac{{\partial }^2}{{\partial y}^2}P)(x,y,z,t^{\prime \prime }){\mathrm{dt}}^{\prime \prime }{\mathrm{dt}}^{\prime }\\ R(x,y,z,t)=v_{p_0}^2\left(\frac{{\partial }^2}{{\partial z}^2}F\right)={\int }_0^t{\int }_0^{t^{\prime }}\left(\frac{{\partial }^2}{{\partial z}^2}P\right)(x,y,z,t^{\prime \prime }){\mathrm{dt}}^{\prime \prime }{\mathrm{dt}}^{\prime }\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{\partial }{\partial t}p=(1+2\eta )v_{\mathrm{nmo}}^2{g}_{2}p+v_{p_0}^2{g}_{1}q\\ \frac{\partial }{\partial 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{\partial }{\partial t}p,$$U=\frac{\partial }{\partial x}p,$$V=\frac{\partial }{\partial y}p,$$Q=\frac{\partial }{\partial t}\left(\frac{(1+2\eta )q-p}{2\eta }\right),$$R=(1+2\eta )\frac{\partial }{\partial z}q.---\left[21\right]$

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

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

然后：

$\frac{\partial }{\partial t}\left(\begin{array}{c}\tilde{P}\\ \tilde{U}\\ \tilde{V}\\ \tilde{Q}\\ \tilde{R}\end{array}\right){=M}_x\frac{\partial }{\partial x}\left(\begin{array}{c}\tilde{P}\\ \tilde{U}\\ \tilde{V}\\ \tilde{Q}\\ \tilde{R}\end{array}\right)+M_y\frac{\partial }{\partial y}\left(\begin{array}{c}\tilde{P}\\ \tilde{U}\\ \tilde{V}\\ \tilde{Q}\\ \tilde{R}\end{array}\right)M_z\frac{\partial }{\partial 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\eta }& 0& 0& 0\\ v_{\mathrm{nmo}}\sqrt{1+2\eta }& 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\eta }& 0& 0\\ 0& 0& 0& 0& 0\\ v_{\mathrm{nmo}}\sqrt{1+2\eta }& 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\eta }}& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ \frac{v_{p_0}}{\sqrt{1+2\eta }}& 0& 0& 0& \frac{v_{p_0\sqrt{2\eta }}}{\sqrt{1+2\eta }}\\ 0& 0& 0& \frac{v_{p_0}\sqrt{2\eta }}{\sqrt{1+2\eta }}& 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\bar{\theta }}{v^2}=\frac{{[f/2+ϵ{\sin }^2\bar{\theta }]}^2-\{\frac{v^2}{{v^2}_{p_0}}-{[1+ϵ\sin }^2\bar{\theta }]+f/2{\}}^2}{{2\mathrm{fv}}^2(ϵ-\delta ){\sin }^2\bar{\theta }}\\ f=1-{\left(\frac{v_{s_0}}{v_{p_0}}\right)}^2\end{array}\right)$

其中，相速v具有两个幅度的根：一个用于准P波，另一个用于准SV波，是波前法线和对称轴之间的角度，以及其它参数定义在方程[1]中。

按照显示在图13中的几何体系，波前法线和对称轴以及其间的角度采取下列形式：

$\left(\begin{array}{c}\overrightarrow{n}=\sin \theta \cos \phi \overrightarrow{i}+\sin \theta \sin \phi \overrightarrow{j}+\cos \theta \overrightarrow{k}\\ \overrightarrow{t}={\sin \theta }_0\cos {\phi }_0\overrightarrow{i}+\sin {\theta }_0\sin {\phi }_0\overrightarrow{j}+\cos {\theta }_0\overrightarrow{k}\\ \cos \bar{\theta }=\frac{\overrightarrow{n}·\overrightarrow{t}}{\left|n\right|\left|\right|t|}=\sin \theta \cos \phi \sin {\theta }_0\cos {\phi }_0+\sin \theta \sin \phi \sin {\theta }_0{\sin \phi }_0+\cos \theta {\cos \theta }_0\end{array}\right)$

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

$\left(\begin{array}{c}\sin \theta \cos \phi =k_{x}v/\omega \\ \sin \theta \sin \phi =k_{y}v/\omega \\ \cos \theta =k_{z}v/\omega \end{array}\right)$

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

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

$+v_{p_0}^4[2f(ϵ-\delta )f_1{f}_2-(f-1)((1+2ϵ)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{\theta }_0\frac{{\partial }^2}{{\partial x}^2}+{\cos }^2{\theta }_0\frac{{\partial }^2}{{\partial z}^2}+\sin 2{\theta }_0\frac{\partial }{\partial x}\frac{\partial }{\partial z})\\ f_2=({\cos }^2{\theta }_0\frac{{\partial }^2}{{\partial x}^2}+\mathrm{si}{n}^2{\theta }_0\frac{{\partial }^2}{{\partial z}^2}-\sin 2{\theta }_0\frac{\partial }{\partial x}\frac{\partial }{\partial 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(专用集成电路)、计算机阵列或网络、或其它适用计算设备。作为进一步的例子，还应该明白，本发明设想，可以尽可能地将任何实施例的一个或多个特征与任何其它实施例的一个或多个特征组合。