CN102062875A - 起伏地表弹性波波动方程正演方法 - Google Patents

Abstract

本发明公开了一种起伏地表弹性波波动方程正演方法，涉及石油地震勘探的地震勘探资料处理与解释领域，利用纵向坐标变换将不规则计算区域转变为新坐标系下的规则计算区域，在零牵引力自由边界的基础上，导出新坐标系下的自由边界条件，并采用隐式差分格式计算自由边界上的地震波场。本方法能处理复杂地形，适应地表起伏的同时提高了计算效率，且能更好地实施自由边界条件。

Description

起伏地表弹性波波动方程正演方法

技术领域

本发明涉及石油地震勘探的地震勘探资料处理与解释领域，确切地说涉及一种用于复杂地表的石油地震勘探中的地下介质正演模拟方法。

背景技术

地震波传播数值模拟是在一定的地震波传播理论基础上，通过数值计算来模拟地震波的传播。尽管数值模拟技术所依据的地震波传播理论是实际地震波传播规律的不同程度的近似，但该技术灵活方便、模拟快速等优点使其在理论研究和实际中得到广泛应用，近年来计算机技术的日新月异进一步推动了数值模拟技术的发展。所以，地震波传播数值模拟是地震数据处理的基石，也是研究地震波传播理论、指导地震数据采集、进行地震资料解释的有力工具。

在地形起伏情况下的地震波传播数值模拟方面，第一个关键问题是如何选择适合不规则地表的数值计算方法，另外一个关键问题是有限模型在地下半无限空间的边界处理问题，如何更好的处理这两个问题是正演模拟方法的关键。

在模拟方法问题上，有限元（FE）方法计算效率低，但处理不规则地表比较方便，成为模拟起伏地表情况下地震波传播的有效方法。为此，发展了一些有限元和其它方法相结合的混合方法，例如，Moczo等（1997）用离散波数方法模拟震源激发和下部介质中地震波的传播，而通过有限元方法来模拟沿起伏地表波的传播（DW-FE）；黄自萍等（2003）用有限元和有限差分（FE-FD）相结合的方法来模拟起伏地表地震波传播。这类方法的主要问题是FE和其它方法计算区域交界处易产生人为反射，吸收边界也不易处理。

边界积分或边界元（BE）方法将散射波场通过地表的一个半解析的积分来表示，其中积分项中Green函数一般在频率波数域中计算。这种方法在研究起伏地表地震波传播时使用比较多，例如，Bouchon(1996)、Pedersen(1994)、Durand（1999）等就是用这种方法研究了地形的区域变化对天然地震波场的影响，分析了地表散射以及地表振动随地形的变化，而Ru-Shan Wu等（1998）将边界元和屏传播算子相结合提出了一种混合模拟方法，在处理起伏地表的同时提高了模拟效率。BE方法的半解析性质决定了该方法不能适用于地表速度变化较大情况，而实际情况是，由于后期地质作用造成浅部地层速度变化更为剧烈，因此限制了BE方法的实际应用。

有限差分（FD）方法计算效率高，在模拟复杂模型中地震波传播时应用最为广泛，但该方法的一个重要缺陷是处理复杂地形比较困难。为此，Tessmer,Kosloff and Behle (1992)提出了一种新的思路，即通过坐标变换将具有起伏地表的模型及弹性波方程变换到新的具有水平地表的坐标系中，在新坐标系中求解弹性波方程，时间上用差分法，空间上横向用Fourier法，纵向采用Chebyshev方法计算波场对空间的导数，使用的吸收边界和自由边界条件均由Gottlieb（1982）所提出。Hestholm and Ruud (1994，1998，1999) 借鉴了Tessmer 等的思路，通过坐标变换将起伏的地表转换成水平地表后，完全通过差分方法求解弹性波方程，在适应地表起伏的同时提高了计算效率。甘文权（2001）利用类似方法对位移波动方程进行了数值求解。

关于地形起伏情况下弹性波模拟的另一个关键问题是自由边界条件的实施方法，这个问题实际上是和模拟方法有机结合在一起的。即使是水平地表，为了研究面波等地震波现象，对自由边界条件实施方式也进行了大量的研究（Jin.,1988;Robertsson,1996;Ohminato,1997;Gottschammer and Olsen, 2001），但每一种方法都存在一定的缺陷，不是理论上不完备，就是和具体的数值方法结合时实施困难。

垂直地表的应力为零的自由边界条件理论上非常简单，但离散波动方程的数值方法对此困难很大（Zahradnik and Moczo,1993）。由于一般的数值计算方法以波场连续性为前提，因此，仅仅在模型之上增加一空气层的办法是不正确的。差分算子越复杂，实现难度越大（Kosloff,1990），而伪谱法几乎难以处理，仅仅令地表应力为零也容易引起算法不稳定（Byliss,1986）。另外，这个问题还和地表的起伏联系在一起，更增加了这个问题的难度。因此，除了数值频散和吸收边界外，自由边界条件的实施成为地震波模拟中又一个顽固的问题。

    公开号为CN101369024，公开日为2009年2月18日的中国专利文献公开了一种地震波动方程生成方法及系统，该地震波动方程生成方法包括：

 获取声波方程数据；

获取地质参数信息；

根据所述的声波方程数据和所述的地质参数信息进行任意差分精细积分，得出地震传播方程；

求出稳定性条件、初始输入条件、边界处理条件；

由所述的地震传播方程以及所述的稳定性条件、初始输入条件、边界处理条件生成地震波动方程。本发明有较高的精度和较好的数值稳定性，可以在几乎不增加计算量的情况下较大地提高了计算精度和稳定性。

所述获取地质参数信息的步骤包括：  

采集地震数据；  

利用promax软件处理所述地震数据得出地震剖面；  

对所述地震剖面进行分析得出对应的地质参数信息。

所述进行任意差分精细积分的步骤包括：    

进行时间有限差分和空间精细积分的交替来进行所述的任意差分精细积分。

所述边界条件是指完全匹配层吸收边界条件，以在完全匹配边界条件下进行任意差分精细积分。

所述地震波动方程为一半解析解。

用任意差分精细积分法导出的完全匹配层边界条件对边界反射有很好的衰减作用。本发明可以为复杂区地震波传播规律研究提供了一个高精度、稳定性能好的数据体。但上述技术方案仍然没有解决自由边界条件的实施问题。

发明内容

为解决上述技术问题，本发明提出了一种用于复杂地表的石油地震勘探中起伏地表弹性波波动方程的正演方法。本方法能处理复杂地形，适应地表起伏的同时提高了计算效率，且能更好地实施自由边界条件。

本发明是通过下述技术方案实现的：

一种起伏地表弹性波波动方程正演方法，其特征在于：利用纵向坐标变换将不规则计算区域转变为新坐标系下的规则计算区域，在零牵引力自由边界的基础上，导出新坐标系下的自由边界条件，并采用隐式差分格式计算自由边界上的地震波场。

所述的利用纵向坐标变换将不规则计算区域转变为新坐标系下的规则计算区域，是通过对二维弹性波方程进行纵向坐标变换实现的，具体步骤如下：

a、对空间垂直方向和水平方向进行不同的剖分，水平方向上将观测区域等分为nx－1段；垂直方向从模型底部到地表进行剖分，共nz个采样点，即从模型底部到地表平分为nz－1等份，这样得到一个nx×nz个采样点的非矩形数据体；

b、通过坐标变换将曲面上拉平到同一水平面上，并令模型底部保持水平，虽然高程值在每一个x下是不同的，但是由于每一个地表下都沿深度进行等分，所以原数据体上的采样点将落在新的“均匀剖分”的矩形网格点上；

在“均匀剖分”的矩形网格点形成的新坐标系下，进行纵向变换：变换保持x轴方向不变，在z轴上根据高程进行变化，变换时引入与高程有关的项，所述高程为由模型底部到地表的深度。

所述在零牵引力自由边界的基础上，导出新坐标系下的自由边界条件是指：通过坐标旋转得到自由边界条件。

与现有技术相比，本发明所能达到的技术效果如下：

本发明采用“利用纵向坐标变换将不规则计算区域转变为新坐标系下的规则计算区域，在零牵引力自由边界的基础上，导出新坐标系下的自由边界条件，并采用隐式差分格式计算自由边界上的地震波场”这样的技术方案，能处理复杂地形，针对山地等复杂区域的波场信息进行模拟，在适应地表起伏的同时提高了计算效率，且能更好地实施自由边界条件。

附图说明

下面将结合说明书附图和具体实施方式对本发明作进一步的详细说明，其中：

图1为纵向坐标变换示意图

附图说明：

图1中（a）为 曲网格坐标系                                                

下的网格剖分

图1中（b）为 进行纵向坐标变换之后的

坐标系下的网格剖分。

具体实施方式

利用纵向坐标变换思路，将不规则计算区域转变为新坐标系下的规则计算区域，在零牵引力自由边界概念基础上，导出了新坐标系下的自由边界条件，并采用隐式差分格式计算自由边界上的地震波场。

（1）所述的利用纵向坐标变换将不规则计算区域转变为新坐标系下的规则计算区域，是通过对二维弹性波方程进行纵向坐标变换实现的，具体步骤如下：

a、对空间垂直方向和水平方向进行不同的剖分，水平方向上将观测区域等分为nx－1段；垂直方向从模型底部到地表进行剖分，共nz个采样点，即从模型底部到地表平分为nz－1等份，这样得到一个nx×nz个采样点的非矩形数据体；

b、通过坐标变换将曲面上拉平到同一水平面上，并令模型底部保持水平，虽然高程值在每一个x下是不同的，但是由于每一个地表下都沿深度进行等分，所以原数据体上的采样点将落在新的“均匀剖分”的矩形网格点上；

我们的数值模拟算法将在这样一个“规则网格”上进行，当然在其中所应用的地震波传播方程、自由边界条件、吸收边界条件也需要与坐标系一同进行相应的纵向坐标变换。这样做最大的优点是，在整个模型区域内部各网格点可以通过相同的计算方式，特别是起伏地表的自由边界条件具有统一形式。

在新坐标系下，这种方法应用的地震波传播方程、自由边界条件、以及吸收边界条件都需要进行纵向变换。变换保持x轴方向不变，在z轴上根据高程进行变化，这样的变换将引入与高程有关的项，其中高程被定义为由模型底部到地表的深度。因此为了方便表示，方法中向上为z轴正方向，以模型底部为z=0。

纵向坐标变换如图1所示。

二维情况下，一阶的速度-应力方程为：

                                 (1-1)

                                  (1-2)

                             (1-3)

                            (1-4)

                                     (1-5)

采用我们的坐标变换方法后，新坐标的方程为：

  

                (1-6)

                   (1-7)

                (1-8)

                (1-9)

                 (1-10)

（2）起伏地表自由边界条件的纵向变换

    自由边界条件可以表示为下式，

                                         (2-1)

                                   (2-2)

通过坐标旋转得到自由边界条件。在

坐标系下的向量

旋转到

坐标下得到新向量

，两个向量之间有如下的关系，

为曲面的外法线方向，

                                            (2-3)

                                           (2-4)

                                (2-5)

                               (2-6)

                                       (2-7)

为

轴与的夹角，新的坐标系可以通过坐标轴旋转得到。于是在坐标系

下，自由边界条件可以表示为，

                                         (2-8)

                                  (2-8)

通过推导整理可以得到在坐标系下的自由边界条件，

                      (2-10)

     (2-11)

（3）起伏地表自由边界条件隐格式实现

时间方向上每一步都需要式

，这样一个矩阵方程。矩阵方程是由对z=nz上的每个x，用二阶差商代替微商的办法构造的，其中

         （3-1）

有当前时刻需要求取的自由表面波场值每一点的两个分量依次排列构成。

          

   (3-2)                                                                

中需要用到z=nz-1层上的当前时刻的波场记录，所以虽然

的形式不变，但是每个时

的值是变化的。

    系数矩阵

只与

，及

、

、

、

在地表上的值有关，而当模型给定后，这些值也都是一定的，所以矩阵

也是不变的。

的表达式如(3-3)。

（3-3）

Claims (3)

1.一种起伏地表弹性波波动方程正演方法，其特征在于：利用纵向坐标变换将不规则计算区域转变为新坐标系下的规则计算区域，在零牵引力自由边界的基础上，导出新坐标系下的自由边界条件，并采用隐式差分格式计算自由边界上的地震波场。

2.根据权利要求1所述的起伏地表弹性波波动方程正演方法，其特征在于：所述的利用纵向坐标变换将不规则计算区域转变为新坐标系下的规则计算区域，是通过对二维弹性波方程进行纵向坐标变换实现的，具体步骤如下：

a、对空间垂直方向和水平方向进行不同的剖分，水平方向上将观测区域等分为nx－1段；垂直方向从模型底部到地表进行剖分，共nz个采样点，即从模型底部到地表平分为nz－1等份，这样得到一个nx×nz个采样点的非矩形数据体；

b、通过坐标变换将曲面上拉平到同一水平面上，并令模型底部保持水平，虽然高程值在每一个x下是不同的，但是由于每一个地表下都沿深度进行等分，所以原数据体上的采样点将落在新的“均匀剖分”的矩形网格点上；

在“均匀剖分”的矩形网格点形成的新坐标系下，进行纵向变换：变换保持x轴方向不变，在z轴上根据高程进行变化，变换时引入与高程有关的项，所述高程为由模型底部到地表的深度。

3.根据权利要求1或2所述的起伏地表弹性波波动方程正演方法，其特征在于：所述在零牵引力自由边界的基础上，导出新坐标系下的自由边界条件是指：通过坐标旋转得到自由边界条件。

Images