阅读说明：本技术 弹性波场数值模拟方法及系统 (Elastic wave field numerical simulation method and system ) 是由 李文杰 刘来详 陈世军 于 2018-06-15 设计创作，主要内容包括：公开了一种弹性波场数值模拟方法及系统。该方法可以包括：步骤1：建立速度模型,确定空间和时间的采样间隔,进行有限差分网格离散化；步骤2：根据弹性波动方程计算内部波场数据；步骤3：根据角点波场数值计算公式计算角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算边界波场数据；步骤5：获得速度模型的波场数据,判断波场边界的反射情况是否满足预定标准,若是,执行步骤6,若否,重新确定波场吸收系数,重复步骤4-5；步骤6：输出最终波场数据。本发明的模型波场值能准确反映出弹性波在二维各向同性介质中传播的真实性,达到提高弹性波数值模拟结果精度的目的。(An elastic wave field numerical simulation method and system are disclosed. The method can comprise the following steps: step 1: establishing a speed model, determining sampling intervals of space and time, and performing finite difference grid discretization; step 2: calculating internal wave field data according to an elastic wave equation; and step 3: calculating angular point wave field data according to an angular point wave field numerical calculation formula; and 4, step 4: calculating boundary wave field data according to a wave field absorption boundary equation and a wave field absorption coefficient; and 5: acquiring wave field data of the velocity model, judging whether the reflection condition of the wave field boundary meets a preset standard, if so, executing the step 6, otherwise, re-determining the wave field absorption coefficient, and repeating the steps 4-5; step 6: and outputting final wave field data. The model wave field value of the invention can accurately reflect the reality of the elastic wave propagating in the two-dimensional isotropic medium, thereby achieving the purpose of improving the accuracy of the numerical simulation result of the elastic wave.)

弹性波场数值模拟方法及系统

技术领域

本发明涉及弹性波波场数值模拟技术领域，更具体地，涉及一种弹性波场数值模拟方法及系统。

背景技术

弹性波波场数值模拟技术是地震勘探中的重要技术，利用弹性波波场数值模拟技术中的正演技术，可以检验经过地震资料处理、地质解释最终得到的地下物性参数以及反演的速度模型是否能真实地反映地下的储油构造，如果离开了弹性波波场数值模拟技术，我们目前的许多地震资料的处理和解释工作将寸步难行，目前广泛应用于实际工作中的弹性波波场数值模拟技术主要是弹性波波场有限差分波场数值模拟技术，在弹性波有限差分波场数值模拟中，由于计算模型的限制，人们所碰到的一个不可避免的问题是如何解决由计算网格边界所引起的人为边界反射能量，这些来自边界的人为反射扭曲了波在无限介质中传播的真实性，因此有必要提出一种波在边界上传播的算法能使往外传播的波的能量透过边界，否则，我们将不得不扩大模型的边界，大大增加计算的工作量。为了解决边界反射问题，几十年来，人们进行了大量的研究，提出了许多不同有关边界波场的计算方法，但至今，仍没有一种十分完美的解决方法，而与声波波场数值模拟相比，弹性波有限差分波场数值模拟中的边界问题更复杂、也更难解决，目前，弹性波有限差分波场数值模拟中的边界问题研究仍是该领域里的一个重要课题，同时也相继诞生了一些针对不同边界问题的算法，其中通过波场吸收消除来自边界反射的方法比较通用且效果也较明显，国内国外都有不少这方面的研究成果出现和文献发表，目前勘探中广泛使用两类吸收边界条件：海绵吸收边界条件和傍轴近似吸收边界条件。海绵吸收边界条件利用黏滞边界或靠近边界的条带范围内对入射波进行衰减；傍轴近似吸收边界条件基于单程波传播的方法原理，将不同近似条件下的单程波方程作为不同边界处的吸收边界条件。在国内，董良国在文献中(董良国.弹性波数值模拟中的吸收边界条件.石油地球物理勘探,1999,34(1):45～56)利用特征分析方法将前人提出的一维情况下无边界反射概念推广至三维各向异性介质弹性波的数值模拟中，得到TI介质中的吸收边界条件，在国外，Clayton和Engquist等人在文献中(RobertClayton,Bjorn Engquist.absorbing boundary conditions for acoustic and elasticwave equations.Bull.Seis.Soc.Am.,1977,67(6):1529～1540)利用弹性波传播方程傍轴近似理论得到了各向同性介质中弹性波数值模拟的吸收边界条件，至今，他们所提出的吸收边界条件仍是在各向同性介质中进行弹性波有限差分数值模拟用来解决边界问题最常用的方法，此后，Clayton和Engquist在文献中(Robert Clayton,Bjorn Engquist,absorbing boundary conditions for wave equations migration.Geophysics,1980,45(3):895～904)又对弹性波有限差分波场数值模拟中的吸收边界条件进行了更深入的研究。在国内，也有不少研究人员对利用有限差分法进行弹性波场数值模拟情况下的吸收边界条件做过较详细的研究，得到了许多有价值的研究成果，他们的算法对某些模型得到了一些较好的结果，本发明第一发明人在文献中(李文杰等.一种新的弹性波数值模拟吸收边界条件.石油地球物理勘探,2009,44(4):501～507)利用二维各向同性介质的弹性波传播方程傍轴近似理论得到了另一种的吸收边界条件,该方法与Clayton和Engquist等人提出的吸收边界方法不一样，经过数值模拟结果的证明，该吸收边界条件具有较好的吸收效果。遍观已有的弹性波波场数值模拟算法，每种方法都有自己各自的特点，但所有这些吸收边界条件都不能完全吸收来自边界的反射波场能量，只是部分吸收来自边界的反射能量。因此，对于弹性波数值模拟来说，确实需要发明一种可以完全吸收来自边界反射能量的吸收边界条件，这样才能提高弹性波数值模拟结果的精度，使数值模拟结果更加符合实际波场的传播规律。因此，有必要开发一种弹性波场数值模拟方法及系统。

公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。

发明内容

本发明提出了一种弹性波场数值模拟方法及系统，其能够通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。

根据本发明的一方面，提出了一种弹性波场数值模拟方法。所述方法可以包括：步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出所述速度模型的最终波场数据。

优选地，所述速度模型的内部各个网格点的各个时间的波场值的表达式为：

其中，

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

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

优选地，所述波场吸收边界方程通过以下方式获得：对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将所述化简方程转换至时间域，获得波场吸收边界方程。

优选地，傅里叶变换后的弹性波动方程为：

优选地，所述化简方程为：

优选地，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，

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

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

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

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

其中，

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

r₁，r₂、r₃，r₄为波场吸收系数，1≥r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄≥0，α、β分别为纵、横波速度。

优选地，通过方程(5a)-(5d)计算所述速度模型四个角点所对应的各个时刻的波场值，

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

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

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

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

其中：

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

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

其中：

θ为波的入射角。

根据本发明的另一方面，提出了一种弹性波场数值模拟系统，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出所述速度模型的最终波场数据。

优选地，所述速度模型的内部各个网格点的各个时间的波场值的表达式为：

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

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

优选地，所述波场吸收边界方程通过以下方式获得：对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将所述化简方程转换至时间域，获得波场吸收边界方程。

优选地，傅里叶变换后的弹性波动方程为：

优选地，所述化简方程为：

优选地，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，

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

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

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

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

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

r₁，r₂、r₃，r₄为波场吸收系数，1≥r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄≥0，α、β分别为纵、横波速度。

优选地，通过方程(5a)-(5d)计算所述速度模型四个角点所对应的各个时刻的波场值，

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

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

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

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

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

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

其中：

θ为波的入射角。

本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的

具体实施方式

中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述，这些附图和具体实施方式共同用于解释本发明的特定原理。

附图说明

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

图1示出了根据本发明的弹性波场数值模拟方法的步骤的流程图。

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

图3a和图3b分别示出了根据本发明的一个实施例的刚性边界条件下传播到120毫秒时波场X分量与Z分量的示意图。

图4a和图4b分别示出了根据本发明的一个实施例的吸收边界条件下传播到120毫秒时波场X分量与Z分量的示意图。

图5示出了根据本发明的一个实施例的双层介质速度模型的示意图。

图6a和图6b分别示出了根据本发明的一个实施例的弹性波双重吸收边界条件下传播到260毫秒时水平位移与垂直位移的示意图。

图7a和图7b分别示出了根据本发明的一个实施例的Clayton吸收边界条件下传播到260毫秒时波场水平位移与垂直位移的示意图。

具体实施方式

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

图1示出了根据本发明的弹性波场数值模拟方法的步骤的流程图。

在该实施例中，根据本发明的弹性波场数值模拟方法可以包括：步骤1：建立速度模型，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出速度模型的最终波场数据。

在一个示例中，速度模型的内部各个网格点的各个时间的波场值的表达式为：

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

α,β分别为纵、横波速度，t为时间。

在一个示例中，波场吸收边界方程通过以下方式获得：对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将化简方程转换至时间域，获得波场吸收边界方程。

在一个示例中，傅里叶变换后的弹性波动方程为：

在一个示例中，化简方程为：

在一个示例中，通过方程(4a)-(4d)计算速度模型的边界上的网格点所对应的各个时刻的波场值，

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

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

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

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

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

r₁，r₂、r₃，r₄为波场吸收系数，1≥r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄≥0，α、β分别为纵、横波速度。

在一个示例中，通过方程(5a)-(5d)计算速度模型四个角点所对应的各个时刻的波场值，

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

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

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

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

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

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

其中：

θ为波的入射角。

具体地，根据本发明的弹性波场数值模拟方法可以包括：

步骤1：建立速度模型，确定炮点(即激发源)位置，选择地震子波，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化。

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

假定密度ρ为常数，这样就可以把公式(6)看成是随空间位置变化的P波和SV波速度的函数。其中λ和μ与纵、横波速度α(x,z)和β(x,z)的关系为公式(6)：

进而将公式(6)改写为矩阵形式，即为公式(1)，计算速度模型的内部各个网格点的各个时刻的波场值为公式(1)，获得内部波场数据。

步骤3：根据角点波场数值计算公式计算速度模型的四个角点的波场值分别为公式(5a)-(5d)，获得角点波场数据，根据炮点位置坐标以及四个角点所对应的位置坐标，计算出sinθ和cosθ，由此得到J矩阵向量中的各项值。

步骤4：根据弹性波动方程与波场吸收系数获得波场吸收边界方程，根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据，包括：

对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程为公式(2)；为了能对公式(2)的左边项进行分解，考虑将公式(2)的左边项变为公式(8)：

其中：

对(8)式进行整理，得到公式(9)：

由(9)式得到(10)式：

令：则有公式(11)：

进而获得公式(12)：

其中：1≥r₁≥0；1≥r₂≥0。

于是，(10)式可以写成下面的推导式子：

可以将(13)式可改写成下面的形式：

I+C(k_z/ω)-A(k_x/ω)-B(k_z/ω)＝0 (14)

其中：

这样，通过因式分解，公式(2)便可以将改写成下面的形式：

令G＝I+C(k_z/ω)+A(k_x/ω)+B(k_z/ω)，在(15)式的两边乘上G^-1，则获得化简方程为公式(3)；将化简方程转换至时间域，获得波场吸收边界方程；根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值为公式(4a)-(4d)，获得边界波场数据。

步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5，其中，本领域技术人员可以根据速度模型的不同来确定预定标准。

步骤6：输出速度模型的最终波场数据。

本方法通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。

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

应用示例1

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

步骤1：建立二维各向同性介质常速度模型，如图2所示，纵波速度为3000米/秒，横波速度为2000米/秒，激发源为纵波震源，位于速度模型的中央。选取主频为50Hz的雷克子波作为我们的地震子波，空间采样间隔为1米，时间采样间隔为0.25ms，然后根据所选的空间采样对所选的速度模型进行网格离散化。

步骤2：根据公式(1)计算速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据。

步骤3：根据角点波场数值计算公式计算速度模型的角点的波场值为公式(5a)-(5d)，获得角点波场数据。

步骤4：根据弹性波动方程与波场吸收系数获得波场吸收边界方程，根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值为公式(4a)-(4d)，获得边界波场数据。

步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5。

步骤6：输出速度模型的最终波场数据。

图3a和图3b分别示出了根据本发明的一个实施例的刚性边界条件(即波场传播到边界时出现全反射的情况)下传播到120毫秒时波场X分量与Z分量的示意图。图4a和图4b分别示出了根据本发明的一个实施例的吸收边界条件下传播到120毫秒时波场X分量与Z分量的示意图。对比图3和图4，很容易看出，本发明所使用的吸收边界条件可以对传播到边界的波场起到很好的吸收作用。

应用示例2

图5示出了根据本发明的一个实施例的双层介质速度模型的示意图。

步骤1：建立二维各向同性层状介质速度模型，如图5所示，第一层介质的纵波速度为2000米/秒，横波速度为1000米/秒，第二层介质的纵波速度为4000米/秒，横波速度为2000米/秒，激发源位于模型顶部以下300米，距左右两边各300米。选取主频为50Hz的雷克子波作为我们的地震子波，空间采样间隔为1米，时间采样间隔为0.25ms，然后根据所选的空间采样对所选的速度模型进行网格离散化。

步骤2：根据公式(1)计算速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据。

步骤3：根据角点波场数值计算公式计算速度模型的角点的波场值分别为公式(5a)-(5d)，获得角点波场数据。

步骤4：根据弹性波动方程与波场吸收系数获得波场吸收边界方程，根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值为公式(4a)-(4d)，获得边界波场数据。

步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5。

步骤6：输出速度模型的最终波场数据。

图6a和图6b分别示出了根据本发明的一个实施例的弹性波双重吸收边界条件下传播到260毫秒时水平位移与垂直位移的示意图。图7a和图7b分别示出了根据本发明的一个实施例的Clayton吸收边界条件下传播到260毫秒时波场水平位移与垂直位移的示意图。在上层介质，这时直达纵波已经传出了模型外，来自两层介质分界面上的反射波也快到达顶部边界；在下层介质，透射波也已传到底部，此时模型的波场分布图上只有透射波和反射波，从图上对比Clayton吸收边界条件和弹性波双重吸收边界条件两种情况下的波场分布状况，可以看到：利用Clayton吸收边界条件得到的波场分布图上出现了来自边界的反射，而利用本发明提出的弹性波双重吸收边界条件却能收到很好的吸收效果。

综上所述，本发明通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。

本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。

根据本发明的另一方面，提出了一种弹性波场数值模拟系统，其上存储有计算机程序，其中，程序被处理器执行时实现以下步骤：步骤1：建立速度模型，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出速度模型的最终波场数据。

在一个示例中，速度模型的内部各个网格点的各个时间的波场值的表达式为：

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

α,β分别为纵、横波速度，t为时间。

在一个示例中，波场吸收边界方程通过以下方式获得：对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将化简方程转换至时间域，获得波场吸收边界方程。

在一个示例中，傅里叶变换后的弹性波动方程为：

在一个示例中，化简方程为：

在一个示例中，通过方程(4a)-(4d)计算速度模型的边界上的网格点所对应的各个时刻的波场值，

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

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

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

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

其中，u、w分别为弹性波场所对应的水平位移和垂直位移； r₁，r₂、r₃，r₄为波场吸收系数，1≥r₁≥0；1≥r₂≥0，1≥r₃≥0；1≥r₄≥0，α、β分别为纵、横波速度。

在一个示例中，通过方程(5a)-(5d)计算速度模型四个角点所对应的各个时刻的波场值，

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

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

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

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

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

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

其中：

θ为波的入射角。

本系统通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。