基于二维各向同性介质的弹性波波场数值的模拟方法与流程

本发明属于地球物理勘探中弹性波波场数值模拟技术领域，具体地涉及基于二维各向同性介质的弹性波广角吸收边界条件的确定方法。

背景技术：

弹性波波场数值模拟技术是地震勘探中的重要技术，利用弹性波波场数值模拟技术中的正演技术，可以检验经过地震资料处理、地质解释最终得到的地下物性参数以及反演的速度模型是否能真实地反映地下的储油构造，如果离开了弹性波波场数值模拟技术，我们目前的许多地震资料的处理和解释工作将寸步难行。目前广泛应用于实际工作中的弹性波波场数值模拟技术主要是弹性波波场有限差分波场数值模拟技术，在弹性波有限差分波场数值模拟中，由于计算模型的限制，人们所碰到的一个不可避免的问题是如何解决由计算网格边界所引起的人为边界反射能量，这些来自边界的人为反射扭曲了波在无限介质中传播的真实性，因此有必要提出一种波在边界上传播的算法能使往外传播的波的能量透过边界，否则，我们将不得不扩大模型的边界，大大增加计算的工作量。

为了解决边界反射问题，几十年来，人们进行了大量的研究，提出了许多不同的有关边界波场的计算方法，但至今仍没有一种十分完美的解决方法。而与声波波场数值模拟相比，弹性波有限差分波场数值模拟中的边界问题更复杂、也更难解决。目前，弹性波有限差分波场数值模拟中的边界问题研究仍是该领域里的一个重要课题，同时也相继诞生了一些针对不同边界问题的算法，其中通过波场吸收消除来自边界反射的方法比较通用且效果也较明显，国内国外都有不少这方面的研究成果出现和文献发表，目前勘探中广泛使用两类吸收边界条件：海绵吸收边界条件和傍轴近似吸收边界条件。海绵吸收边界条件利用黏滞边界或靠近边界的条带范围内对入射波进行衰减；傍轴近似吸收边界条件基于单程波传播的方法原理，将不同近似条件下的单程波方程作为不同边界处的吸收边界条件。在国内，董良国在文献中(董良国.弹性波数值模拟中的吸收边界条件.石油地球物理勘探，1999，34(1):45～56)利用特征分析方法将前人提出的一维情况下无边界反射概念推广至三维各向异性介质弹性波的数值模拟中，得到ti介质中的吸收边界条件。在国外，clayton和engquist等人在文献中(robertclayton，bjornengquist.absorbingboundaryconditionsforacousticandelasticwaveequations.bull.seis.soc.am.，1977，67(6):1529～1540)利用弹性波传播方程傍轴近似理论得到了各向同性介质中弹性波数值模拟的吸收边界条件。至今，他们所提出的吸收边界条件仍是在各向同性介质中进行弹性波有限差分数值模拟用来解决边界问题最常用的方法。

此后，clayton和engquist在文献中(robertclayton，bjornengquist，absorbingboundaryconditionsforwaveequationsmigration.geophysics，1980，45(3):895～904)又对弹性波有限差分波场数值模拟中的吸收边界条件进行了更深入的研究。在国内，也有不少研究人员对利用有限差分法进行弹性波场数值模拟情况下的吸收边界条件做过较详细的研究，得到了许多有价值的研究成果，他们的算法对某些模型得到了一些较好的结果，本发明第一发明人在文献中(李文杰等.一种新的弹性波数值模拟吸收边界条件.石油地球物理勘探，2009，44(4):501～507)利用二维各向同性介质的弹性波传播方程傍轴近似理论得到了另一种的吸收边界条件，该方法与clayton和engquist等人提出的吸收边界方法不一样，经过数值模拟结果的证明，该吸收边界条件具有较好的吸收效果。

但遍观已有的弹性波波场数值模拟算法，虽然较好地解决了模型四边的反射能量问题，却没有解决好模型角点的反射能量问题，从模拟效果来看，这些边界条件不能有效消除来自任意角度角点的反射能量，所以，基于各向同性介质的弹性波傍轴近似吸收边界条件有待进一步完善。

技术实现要素：

本发明的目的是为了改进现有的基于二维各向同性介质的弹性波有限差分波场数值模拟技术，进而达到提高地震勘探技术精度的目的。目前基于二维各向同性介质的弹性波有限差分波场数值模拟技术所用的计算速度模型角点波场的算法存在缺陷，这些算法不能有效消除来自任意角度角点的反射能量，而只能消除来自45度角(即角点的入射波场方向与模型的顶边之间的夹角为45度角)的角点反射能量，从而导致计算得到的模型波场值不能正确反映波在无限介质中传播的真实性。

为了解决这个问题，在本发明第一发明人提出的弹性波吸收边界条件的基础上，本发明提出了弹性波广角吸收边界条件的算法，利用本发明的波场算法，无论激发源位于模型内的什么位置，都可以很好地解决速度模型任意边界角点的波场反射问题，使得到的模型波场值准确地反映波在无限介质中传播的真实性，这样便可以提高波场数值模拟的准确性，从而达到提高地震勘探精度的目的。

本发明提供了一种基于二维各向同性介质的弹性波波场数值的模拟方法，该方法包括以下步骤：

根据给定的速度模型，确定波场激发源的位置所对应的坐标值；

选择地震子波，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化；

计算速度模型内部各个网格点的各个时刻的波场值，计算速度模型四条边上的网格点所对应的各个时刻的波场值，计算速度模型四个角点所对应的各个时刻的波场值；

由四个角点的波场值、四条边上的网格点的波场值以及速度模型内部各个网格点的波场值得到整个速度模型的各个网格点在各个时刻的波场值。

进一步地，利用公式(3)计算速度模型内部各个网格点的各个时刻的波场值：

其中：u、w分别为弹性波场所对应的水平位移场和垂直位移场；

α、β分别为纵、横波速度，x和z是二维介质情况下沿x轴和z轴的坐标。

进一步地，利用公式(10a-10d)计算速度模型四个角点所对应的各个时刻的波场值：

右下角角点的波场计算方程为：

左上角角点的波场计算方程为：

右上角角点的波场计算方程为：

左下角角点的波场计算方程为：

其中：u、w分别为弹性波场所对应的水平位移场和垂直位移场；

α、β分别为纵、横波速度，x和z是二维介质情况下沿x轴和z轴的坐标；

其中：

θ为波的入射角。

进一步地，根据激发源的位置坐标以及四个角点所对应的位置坐标，计算出sinθ和cosθ，由此得到c矩阵向量中的各项值。

进一步地，利用公式(10a-10d)计算速度模型四条边上的网格点所对应的各个时刻的波场值：

速度模型底边所对应的吸收边界方程为：

速度模型顶边所对应的吸收边界方程为：

速度模型左边所对应的吸收边界方程为：

速度模型右边所对应的吸收边界方程为：

其中，式(11a)和(11b)所对应的参数项为：

其中，式(11c)和(11d)所对应的参数项为：

其中，α、β分别为纵、横波速度，x和z是二维介质情况下沿x轴和z轴的坐标。

将本发明提出的弹性波广角吸收边界条件方法应用于二维各向同性介质的弹性波波场数值模拟，可以有效地消除来自有限的速度模型边界及任意角度的边界角点所产生的反射能量。本发明所提供的边界波场以及边界节点波场的算法能使往外传播的波的能量透过边界和边界角点，使计算得到的模型波场值能准确地反映出弹性波在二维各向同性无限介质中传播的真实性，从而达到提高地震勘探技术精度的目的。

附图说明

通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施方式中，相同的参考标号通常代表相同部件。

图1显示了本发明基于二维各向同性介质的弹性波波场数值的模拟方法的流程图。

图2显示了根据一个实施例的二维各向同性介质常速度模型。

图3显示了根据一个实施例的刚性边界条件下t＝120ms时的波场传播快照。

图4显示了根据一个实施例的吸收边界条件下t＝120ms时的波场传播快照。

图5显示了根据一个实施例的利用常规吸收边界条件方法得到的t＝120ms时的波场传播快照。

图6显示了根据一个实施例的利用本发明的方法得到的t＝120ms时的波场传播快照。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施方式。虽然附图中显示了本公开的优选实施方式，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

本发明提出的弹性波广角吸收边界条件算法是对现有的基于二维各向同性介质弹性波傍轴近似理论的吸收边界条件的一种改进，它解决了来自速度模型内任意角度的角点的反射能量问题，使通过弹性波数值模拟得到的模型波场值能正确地反映弹性波在二维各向同性无限介质中传播的真实性。

参照图1，描述本发明提供了一种基于二维各向同性介质的弹性波波场数值的模拟方法，该方法包括以下步骤：

根据给定的速度模型，确定波场激发源的位置(即炮点位置)所对应的坐标值。优选地，速度模型可以是各向同性弹性介质的速度模型。

选择地震子波，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化(子波一般选用雷克子波。优选地，时间采样间隔δt与空间采样间隔δx、δz满足如下的公式：其中：h＝min(δx,δz),α、β分别为纵、横波速度)。

计算速度模型内部各个网格点的各个时刻的波场值，计算速度模四条边上的网格点所对应的各个时刻的波场值，计算速度模型四个角点所对应的各个时刻的波场值；

由四个角点的波场值、四条边上的网格点的波场值以及速度模型内部各个网格点的波场值得到整个速度模型的各个网格点在各个时刻的波场值。

根据本发明具体实施例，从二维各向同性介质的弹性波传播方程出发，通过系列推导得到了速度模型任意角度的角点的弹性波波场计算公式，其方法原理如下：

在各向同性介质条件下，假设x和z是二维介质情况下沿x轴和z轴的坐标，x轴的方向向右，z轴的方向朝下，可以利用两个共轭的二阶偏微分方程来描述介质中p波的运动以及垂向偏振的sv波的运动，在这里不考虑沿水平方向偏振的sh波，并假定u、w分别为水平方向与垂直方向的位移，ρ为介质的密度，t为时间，λ和μ为具体介质中的拉梅系数，那么可以得到关于各向同性非均匀介质的弹性波动方程如下式表示：

现在假定密度ρ为常数，这样就可以把上面的方程组看成是随空间位置变化的p波和sv波速度的函数。其中λ和μ与纵、横波速度α(x,z)和β(x,z)的关系如下所示：

这样便可以把(1)式写成下列形式：

其中：分别为弹性波场所对应的水平位移场和垂直位移场；

α、β分别为纵、横波速度。

方程(3)的傅氏变换可以写成下列形式：

为了能利用傍轴近似，将式(4)的左边项写成下面的形式：

d2(kz/ω)²+h(kz/ω)(kx/ω)+d1(kx/ω)²＝i(5)

对(5)式进行整理，得：

令：

其中：

很显然：

令：

其中：

θ为波的入射角；

根据激发源的位置坐标以及四个角点所对应的位置坐标，计算出sinθ和cosθ，由此得到c矩阵向量中的各项值。

令：则有：

于是，(6)式可以写成下面式子：

由(7)式便可以将(4)式改写成下面的式子：

令d＝i+c(kz/ω)(kx/ω)+a(kz/ω)+b(kx/ω)，在(8)式的两边乘上d^-1，则得到下面式子：

将(9)式变换到时间域，可以得到下面的方程：

式(10a)为右下角角点的波场计算式，我们得到其他三个角点的波场计算式如下：

左上角角点的波场计算方程为：

右上角角点的波场计算方程为：

左下角角点的波场计算方程为：

在(10a-10d)各个表达式中，其中：u、w分别为弹性波场所对应的水平位移场和垂直位移场；

α、β分别为纵、横波速度。

本发明计算边界的波场采用本发明第一发明人基于二维各向同性介质弹性波傍轴近似理论提出的吸收边界条件，其对应的四边方程如下所示：

速度模型底边所对应的吸收边界方程为：

速度模型顶边所对应的吸收边界方程为：

其中，式(11a)和(11b)所对应的参数项为：

速度模型左边所对应的吸收边界方程为：

速度模型右边所对应的吸收边界方程为：

其中，式(11c)和(11d)所对应的参数项为：

为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。下面结合附图和实施例对本发明进一步说明。

实现弹性波波场正演，首先需要确定速度模型，图2为本实施例所选的二维各向同性介质常速度模型，纵波速度为3000米/秒，横波速度为2000米/秒，激发源为纵波震源，位于速度模型的中央。

在正演过程中，选取主频为50hz的雷克子波作为本实施例的地震子波，空间采样间隔为1米，时间采样间隔为0.25ms，然后根据所选的空间采样对所选的速度模型进行网格离散化。

图3为选用刚性边界条件(即波场传播到边界时出现全反射的情况)波场传播到120毫秒时整个模型的波场快照；图4为使用本发明的弹性波广角吸收边界条件，波场传播到120毫秒时刻时整个模型的波场快照。对比图3和图4，很容易看出，本发明所使用的吸收边界条件可以对传播到边界的波场起到很好的吸收作用。

为了展示本发明在波场正演中的优越性，我们将模型1中的激发源位置进行一下变动：将激发源的位置移到距模型顶边12.5米、距模型的左边102.5米的位置。图5为使用常规吸收边界条件所得到的波场传播到120毫秒时刻时整个模型的波场快照；图6为使用本发明的广角吸收边界条件所得到的波场传播到120毫秒时刻时整个模型的波场快照。

对比图5和图6，不难发现常规吸收边界条件不能对来自角点的波场反射起到很好的吸收作用，而使用本发明的广角吸收边界条件则可以对来自角点的波场反射起到较好的吸收作用。

本发明提出了弹性波广角吸收边界条件波场算法，利用本发明的波场算法可以有效地消除来自边界及任意角度的边界角点的反射能量，使计算的模型波场值能准确地反映波在无限介质中传播的真实性，从而达到提高勘探精度的目的。

以上已经描述了本公开的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。