CN105093265A - 一种模拟地震波在ti介质中传播规律的方法 - Google Patents

Abstract

本发明提供了一种模拟地震波在TI介质中传播规律的方法，属于地震勘探领域。本方法包括：步骤1：根据待求取的问题确定六阶矩阵系数；步骤2：求解一阶速度——应力方程；步骤3：求取步骤2中的褶积微分算子；步骤4：求取窗函数系数，对褶积微分算子进行加窗截断；步骤5：根据地质体，给定模型参数、震源参数、震源频率，对模型进行离散化；步骤6：进行时间迭代，得到整个区域内的速度和应力值，若设定时间到，则模拟过程结束，输出波场快照和地震记录。

Description

一种模拟地震波在TI介质中传播规律的方法

技术领域

本发明属于地震勘探领域，具体涉及一种模拟地震波在TI介质中传播规律的方法，可模拟地震波在横向各向同性(TI)介质中传播规律，用于波动方程的正演模拟，来模拟地震波在介质中的传播规律。

背景技术

研究地震波在复杂介质中的传播规律，是地球物理学中的一项重要内容，是帮助人们认识地球介质的有效模拟手段。

地震波数值模拟是根据已知的地下结构和物性参数来模拟研究地震波在地下各种介质中的传播规律，在地震勘探和天然地震领域有广泛的应用，而实际地下介质的各向异性是普遍存在的.描述地下介质波的传播过程可以用波动方程，因此基于波动方程的数值解法如有限差分法、有限元法、伪谱法、谱元法等在各向异性介质模拟中均有广泛应用.这些方法都有各自的优缺点：有限差分法理论基础为基于Taylor展开，计算速度快，效率高，但对于非规则区域问题较难处理；有限元法基于变分理论，模拟精度高，可以处理复杂形状的地质体，但计算过程复杂，计算量大；伪谱法也是解决偏微分方程的数值解法之一，在计算过程中引入了快速Fourier变换，精度可以达到谱的精度，但由于Fourier变换需要全局信息，伪谱法在处理复杂介质中波场传播问题时显得力不从心；将有限元和伪谱法优点结合起来的谱元法也存在着譬如边界和新方法难以匹配等问题。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种可用于TI介质正演模拟的褶积微分方法，更好地描述波在介质中的传播特征，并降低地震波数值模拟中数值频散现象。

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

一种模拟地震波在TI介质中传播规律的方法，包括：

步骤1：根据待求取的问题确定六阶矩阵系数：

步骤2：求解一阶速度——应力方程；

步骤3：求取步骤2中的褶积微分算子；

步骤4：求取窗函数系数，对褶积微分算子进行加窗截断；

步骤5：根据地质体，给定模型参数、震源参数、震源频率，对模型进行离散化；

步骤6：进行时间迭代，得到整个区域内的速度和应力值，若设定时间到，则模拟过程结束，输出波场快照和地震记录。

所述步骤1是这样实现的：

将TI介质各向异性参数表示为六阶矩阵形式：

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{21}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right],$

其中C₁₂＝C₁₁-2C₆₆.式中：C为各向异性介质参数.在各向同性情况下，将该六阶矩阵弱化为用拉梅常数λ和剪切模量μ来描述介质，此时各参数之间满足如下关系式：C₁₁＝C₃₃＝λ+2μ，C₁₃＝λ，C₅₅＝μ。

所述步骤2是这样实现的：

考虑x-z平面的二维问题，速度-应力方程为

$\left(\frac{{\∂\mathrm{\σ}}_{\mathrm{xx}}}{\∂t}\right)=C_{11}{D}_x{v}_x+C_{13}{D}_z{v}_z---(1-a)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{zz}}}{\∂t}=C_{13}{D}_x{v}_x+C_{33}{D}_z{v}_z---(1-b)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{xz}}}{\∂t}=C_{55}(D_z{v}_x+D_x{v}_z)---(1-c)$

$\frac{{\∂v}_x}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xx}}+D_z{\mathrm{\σ}}_{\mathrm{xz}})---(1-d)$

$\frac{{\∂v}_z}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xz}}+D_z{\mathrm{\σ}}_{\mathrm{zz}})---(1-e)$

(6)

式中，σ_xx、σ_xz、σ_zz为应力分量；v_x、v_z分别为速度水平分量和垂直分量，ρ为密度，D_x、D_z分别为空间沿水平和垂直方向的褶积微分算子。

所述步骤3是这样实现的：

将一个函数u(x)表示为奇异核函数与元函数的褶积形式：

$u\left(x\right)={\∫}_{-\∞}^{\∞}\mathrm{\δ}(x-t)u\left(t\right)\mathrm{dt}---\left(2\right)$

用Shannon奇异核函数来逼近(2)式中的δ函数：

${\mathrm{\δ}}_{\mathrm{\Δx}}\left(x\right)=\frac{\sin (\mathrm{\πx}/\mathrm{\Δx})}{\mathrm{\πx}/\mathrm{\Δx}}---\left(7\right)$

其中，Δx为空间步长.当Δx→0时，所得函数即为Delta函数.

对于一个函数的离散序列u(x_m)，通过对奇异核函数求导，得到该函数在全空间的插值及微分形式：

其中表示最接近x的格点，m表示离散点的个数，q表示空间微分阶数，W为算子半宽度.根据公式(2)、(3)、(4)，得到基于Shannon奇异核的一阶褶积微分算子为：

$d_1\left(\mathrm{m\Δx}\right)=\left\{\begin{array}{cc}{\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left(\mathrm{m\Δx}\right),& m=\±1,\±2,...\\ 0,& m=0\end{array}\right.---\left(9\right)$

其交错网格形式为：

${\hat{d}}_s\left(\mathrm{m\Δx}\right)={\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left[(m+1/2)\mathrm{\Δx}\right],m=-W,-W+1,...,W-1---\left(10\right)$

公式(5)和(6)得到的即为公式(1)中的离散形式的褶积微分算子D_x和D_z，其中公式(5)适用于普通网格的地震波数值模拟，公式(6)适用于交错网格情况下的地震波数值模拟。

所述步骤6是这样实现的：按步骤2给出的公式计算应力及速度的离散值：其中公式(1)中(1-a)、(1-b)、(1-c)为由速度分量v_x、v_z求取应力分量σ_xx、σ_xz、σ_zz的公式，(1-d)、(1-e)为由应力分量σ_xx、σ_xz、σ_zz求取速度分量v_x、v_z的公式，这些变量的初始值均赋为0，步骤5中的震源参数与时间有关，随着时间迭代，通过这五个公式的迭代，得到整个区域内的速度和应力值，若设定时间到，则模拟过程结束，输出波场快照和地震记录。

与现有技术相比，本发明的有益效果是：

(1)更有利于描述波在介质中的传播特征；

(2)能够降低地震波数值模拟中数值频散现象。

附图说明

图1是TI介质空间交错网格示意图。

图2-1是优化褶积微分法水平分量。

图2-2是优化褶积微分法垂直分量。

图2-3是有限差分法水平分量。

图2-4是有限差分法垂直分量。

图3是有限差分方法与褶积微分方法在检波器位于(500m，400m)处的地震记录对比图(x为水平分量，z为垂直分量)。

图4是本方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

波动方程的空间微分可以采用多种形式，除了以上所述的有限差分、有限元、伪谱法、谱元法等之外，本发明提出一种褶积微分算子对波动方程的空间进行微分。在20世纪70年代早期，就有学者提出褶积微分算子的概念，但将其应用到地震波场模拟中比较晚。Holberg运用Fourier变换，对褶积微分算子进行了设计，并在三维弹性波数值模拟中取得了较好的结果。Zhou等、张中杰等和戴志阳等也把这一方法进行了深入研究和发展。龙桂华等在前人的工作基础上，根据函数分布理论，引入奇异核褶积微分算子，通过一组经过加窗处理的微分算子和波场变量进行褶积，来求取波场变量关于空间微分。褶积微分算子和伪谱法相比，其优势在于褶积微分算子是一种优化短算子，更能突出空间微分本身的局部属性；而伪谱法在计算微分时须利用全局信息，就复杂非均匀介质中波的传播问题而言，与物理因果律相悖。褶积微分算子和有限差分算子相比，精度更高，能够弥补有限差分算子描述复杂介质的不足。因此褶积微分算子应用于波动方程的空间离散或空间微分近似是非常有效的。

如图4所示，本发明是通过以下技术方案来解决的，以模拟TI介质中波的传播规律为例，具体包含如下步骤：

步骤1：根据待求取的问题确定六阶矩阵系数：

TI介质各向异性参数可以表示为六阶矩阵形式：

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{21}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right],$

其中C₁₂＝C₁₁-2C₆₆.式中：C为各向异性介质参数.在各向同性情况下，该六阶矩阵可以弱化为用拉梅常数λ(弹性模量)和μ(剪切模量)来描述介质，此时各参数之间满足如下关系式：C₁₁＝C₃₃＝λ+2μ，C₁₃＝λ，C₅₅＝μ.

步骤2：求解一阶速度——应力方程：

考虑x-z平面的二维问题，速度-应力方程为

$\left(\frac{{\∂\mathrm{\σ}}_{\mathrm{xx}}}{\∂t}\right)=C_{11}{D}_x{v}_x+C_{13}{D}_z{v}_z---(1-a)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{zz}}}{\∂t}=C_{13}{D}_x{v}_x+C_{33}{D}_z{v}_z---(1-b)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{xz}}}{\∂t}=C_{55}(D_z{v}_x+D_x{v}_z)---(1-c)$

$\frac{{\∂v}_x}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xx}}+D_z{\mathrm{\σ}}_{\mathrm{xz}})---(1-d)$

$\frac{{\∂v}_z}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xz}}+D_z{\mathrm{\σ}}_{\mathrm{zz}})---(1-e)$

(11)

式中，σ_xx、σ_xz、σ_zz为应力分量；v_x、v_z分别为速度水平分量和垂直分量，ρ为密度，D_x、D_z分别为空间沿水平和垂直方向的褶积微分算子.本发明提出的速度-应力方程与常规速度-应力方程不同之处在于：利用褶积运算代替了常规速度-应力方程中的求导运算。以公式(1-a)为例(以下公式中类似)，D_xv_x表示对v_x在空间上沿水平方向进行褶积运算，在常规速度-应力方程中，需要对v_x在空间上沿水平方向求导，用数学表达式则记为本发明提出的方法避免了空问上的求导运算。

公式(1)为应力与速度之间的关系，若使用程序来具体模拟计算时，上述公式中所有变量均取离散格式。左边均为分量对时间求导，离散格式下可用普通时间差分进行计算。

数值模拟可在普通网格或者交错网格上进行。若采用普通网格形式，这些变量均取网格点上的数值；若采用交错网格形式，这些变量有些需采用网格点上的数值，有些需采用半网格点上的数值(即取周围网格点的平均值)，本例中TI介质的各参数分布情况如图1所示，其中实线与实线的交点为网格点，实线与虚线的交点为半网格点(取值为该点前后两个点的平均值)，虚线与虚线的交点也为半网格点(取值为该点周围的四个点的平均值)。

步骤3：求取步骤2中的褶积微分算子。

具体计算方法为：

将一个函数u(x)根据函数分布理论，可表示为奇异核函数(可取Delta(δ)函数)与元函数(可取u(x)本身)的褶积形式：

$u\left(x\right)={\∫}_{-\∞}^{\∞}\mathrm{\δ}(x-t)u\left(t\right)\mathrm{dt}---\left(2\right)$

而(2)式中的δ函数可用Shannon奇异核函数来进行逼近：

${\mathrm{\δ}}_{\mathrm{\Δx}}\left(x\right)=\frac{\sin (\mathrm{\πx}/\mathrm{\Δx})}{\mathrm{\πx}/\mathrm{\Δx}}---\left(12\right)$

其中，Δx为空间步长.当Δx→0时，所得函数即为Delta函数.

对于一个函数的离散序列u(x_m)，通过对奇异核函数求导，可得该函数在全空间的插值及微分形式：

其中表示最接近x的格点，m表示离散点的个数，q表示空间微分阶数，W为算子半宽度.根据公式(2)、(3)、(4)，可得基于Shannon奇异核的一阶褶积微分算子为：

$d_1\left(\mathrm{m\Δx}\right)=\left\{\begin{array}{cc}{\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left(\mathrm{m\Δx}\right),& m=\±1,\±2,...\\ 0,& m=0\end{array}\right.---\left(14\right)$

其交错网格形式为：

${\hat{d}}_s\left(\mathrm{m\Δx}\right)={\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left[(m+1/2)\mathrm{\Δx}\right],m=-W,-W+1,...,W-1---\left(15\right)$

公式(5)和(6)得到的即为公式(1)中的离散形式的褶积微分算子D_x(水平方向)和D_z(垂直方向)，其中公式(5)适用于普通网格的地震波数值模拟，公式(6)适用于交错网格情况下的地震波数值模拟。

步骤4：求取窗函数系数，对微分算子进行加窗截断：

为了消除因截断算子引起的Gibbs现象，可对算子进行加窗处理。以交错网格情况为例，对公式(6)进行加窗：

$\stackrel{\‾}{d_s}\left(\mathrm{m\Δx}\right)={\hat{d}}_s\left(\mathrm{m\Δx}\right)w(m+1/2),m=-W,-W+1,...,W-1---\left(16\right)$

所加窗函数为Hanning窗：

$w\left(m\right)={[2\mathrm{\α}-1+2(1-\mathrm{\α}){\cos }^2\frac{\mathrm{m\π}}{2(W+2)}]}^{\frac{\mathrm{\β}}{2}}---\left(17\right)$

其中，α、β为刻画窗函数性质的参数.当算子半宽度W一定时，通过调节窗函数的系数α和β可进一步调节数值模拟的精度.关于α、β的求取数学推导过程较为复杂，与本发明内容关联度较小，此处不在赘述。本算例中，算子长度W为9，用最优化方法(如模拟退火方法即可)求得的窗函数权系数为α＝0.4015，β＝2.7069。当然也可以根据目标问题的需要选择合适的算子长度、通过优化方法来求取新的窗函数权系数。

步骤5：根据地质体，对模型进行离散化。给定模型参数、震源参数、震源频率。

本例中，设定模型参数为C₁₁＝25.5×10⁹N.m，C₁₃＝1.0×10⁹N.m，C₃₃＝18.4×10⁹N.m，C₅₅＝5.6×10⁹N.m，密度为2440kg/m³.模型网格点数为301×301，网格大小为Δx＝Δz＝10m，时间采样间隔为1ms，震源位于模型正中央，坐标为(1500m，1500m).震源采用垂直集中力源，为主频25Hz的雷克子波，地震记录的接收排列位于地下500m深度水平排列.以上参数只需符合常规地质体的物理性质即可，本发明对以上参数的取值范围不做限制。

此步骤与常规波动方程正演模拟方法类似。

步骤6：按步骤2给出的公式计算应力及速度的离散值。其中公式(1)中(1-a)、(1-b)、(1-c)为由速度分量v_x、v_z求取应力分量σ_xx、σ_xz、σ_zz的公式，(1-d)、(1-e)为由应力分量σ_xx、σ_xz、σ_zz求取速度分量v_x、v_z的公式。这些变量的初始值均赋为0，步骤5中震源参数与时间有关，随着时间迭代，通过这五个公式的迭代，得到整个区域内的速度和应力值。若设定时间到，则模拟过程结束，输出波场快照和地震记录。

图2-1至图2-4为该介质400ms时弹性波传播的波场快照，其中图2-1、图2-2为9点优化褶积微分算子法的模拟结果，前者为水平分量，后者为垂直分量，图2-3、2-4为相同算子长度(8阶)有限差分法的模拟结果，前者为水平分量，后者为垂直分量.9点优化褶积微分算子法的模拟结果与8阶交错网格有限差分的模拟结果基本一致，不同之处在于交错网格有限差分法在qSV波模拟时出现了较为明显的频散现象，而用优化褶积微分算子进行模拟时频散得到了很好的压制.这表明，将这种优化褶积微分算子用于TI介质中弹性波的模拟是可行的，而且精度高，在压制数值频散效应方面要优于相同算子长度的有限差分法.因此在计算过程中可以根据需要加大网格间距从而减少计算量，进而提高计算效率.

为进一步分析优化褶积微分算子法的精度，算例中选取了接收器的位置为(500m，400m)处的单道地震记录，如图3所示为有限差分方法与褶积微分方法在检波器位于(500m，400m)处的地震记录对比图(x为水平分量，z为垂直分量)，图中实线为8阶有限差分法的模拟结果，虚线为9点优化褶积微分算子法的模拟结果，x为水平分量，z为垂直分量.可以看出，除在振幅的波峰和波谷处两种方法略有偏差外，其它地方均吻合得较好.此外，在0.9s后交错网格的有限差分法地震记录图上出现了明显的频散现象，而优化褶积微分算子法的记录图上却并未观察到频散效应.这与前面图2-1至图2-4中的波场快照相吻合，验证了优化褶积微分算子法模拟精度高于同算子长度的有限差分算子方法.而在实际模拟中，还可以根据精度的需要灵活地调节窗函数的系数，从而调节褶积微分算子的精度.这说明，这种褶积微分算子法除了精度高，还具有使用灵活的特点。

本方法通过对窗函数参数进行优化实现了对基于Shannon奇异核理论的交错网格褶积微分算子的优化过程.应用这种优化褶积微分算子方法对各向异性介质进行了数值模拟，分析了弹性波在此类介质中的传播特征，并与高阶交错网格有限差分方法进行了对比.数值实验结果表明，该方法适用于各向异性介质中弹性波场模拟，精度高，稳定性好，是一种研究复杂介质中地震波传播的有效数值方法。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (5)

1.一种模拟地震波在TI介质中传播规律的方法，其特征在于：所述方法包括：

步骤1：根据待求取的问题确定六阶矩阵系数：

步骤2：求解一阶速度——应力方程；

步骤3：求取步骤2中的褶积微分算子；

步骤4：求取窗函数系数，对褶积微分算子进行加窗截断；

步骤5：根据地质体，给定模型参数、震源参数、震源频率，对模型进行离散化；

步骤6：进行时间迭代，得到整个区域内的速度和应力值，若设定时间到，则模拟过程结束，输出波场快照和地震记录。

2.根据权利要求1所述的模拟地震波在TI介质中传播规律的方法，其特征在于：所述步骤1是这样实现的：

将TI介质各向异性参数表示为六阶矩阵形式：

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{21}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right],$

其中C₁₂＝C₁₁-2C₆₆.式中：C为各向异性介质参数.在各向同性情况下，将该六阶矩阵弱化为用拉梅常数λ和剪切模量μ来描述介质，此时各参数之间满足如下关系式：C₁₁＝C₃₃＝λ+2μ，C₁₃＝λ，C₅₅＝μ。

3.根据权利要求2所述的模拟地震波在TI介质中传播规律的方法，其特征在于：所述步骤2是这样实现的：

考虑x-z平面的二维问题，速度-应力方程为

$\left(\frac{{\∂\mathrm{\σ}}_{\mathrm{xx}}}{\∂t}\right)=C_{11}{D}_x{v}_x+C_{13}{D}_z{v}_z---(1-a)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{zz}}}{\∂t}=C_{13}{D}_x{v}_x+C_{33}{D}_z{v}_z---(1-b)$

$\frac{{\∂\mathrm{\σ}}_{\mathrm{xz}}}{\∂t}=C_{55}(D_z{v}_x+D_x{v}_z)---(1-c)$

$\frac{{\∂v}_x}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xx}}+D_z{\mathrm{\σ}}_{\mathrm{xz}})---(1-d)$

$\frac{{\∂v}_z}{\∂t}=\frac{1}{\mathrm{\ρ}}(D_x{\mathrm{\σ}}_{\mathrm{xz}}+D_z{\mathrm{\σ}}_{\mathrm{zz}})---(1-e)$

(1)

式中，σ_xx、σ_xz、σ_zz为应力分量；v_x、v_z分别为速度水平分量和垂直分量，ρ为密度，D_x、D_z分别为空间沿水平和垂直方向的褶积微分算子。

4.根据权利要求3所述的模拟地震波在TI介质中传播规律的方法，其特征在于：所述步骤3是这样实现的：

将一个函数u(x)表示为奇异核函数与元函数的褶积形式：

$u\left(x\right)={\∫}_{-\∞}^{\∞}\mathrm{\δ}(x-t)u\left(t\right)\mathrm{dt}---\left(2\right)$

用Shannon奇异核函数来逼近(2)式中的δ函数：

${\mathrm{\δ}}_{\mathrm{\Δx}}\left(x\right)=\frac{\sin (\mathrm{\πx}/\mathrm{\Δx})}{\mathrm{\πx}/\mathrm{\Δx}}---\left(2\right)$

其中，Δx为空间步长.当Δx→0时，所得函数即为Delta函数.

对于一个函数的离散序列u(x_m)，通过对奇异核函数求导，得到该函数在全空间的插值及微分形式：

其中表示最接近x的格点，m表示离散点的个数，q表示空间微分阶数，W为算子半宽度.根据公式(2)、(3)、(4)，得到基于Shannon奇异核的一阶褶积微分算子为：

$d_1\left(\mathrm{m\Δx}\right)=\left\{\begin{array}{cc}{\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left(\mathrm{m\Δx}\right),& m=\±1,\±2,...\\ 0,& m=0\end{array}\right.---\left(4\right)$

其交错网格形式为：

${\hat{d}}_s\left(\mathrm{m\Δx}\right)={\mathrm{\δ}}_{\mathrm{\Δx}}^{\′}\left[(m+1/2)\mathrm{\Δx}\right],m=-W,-W+1,...,W-1---\left(5\right)$

公式(5)和(6)得到的即为公式(1)中的离散形式的褶积微分算子D_x和D_z，其中公式(5)适用于普通网格的地震波数值模拟，公式(6)适用于交错网格情况下的地震波数值模拟。

5.根据权利要求4所述的模拟地震波在TI介质中传播规律的方法，其特征在于：所述步骤6是这样实现的：按步骤2给出的公式计算应力及速度的离散值：其中公式(1)中(1-a)、(1-b)、(1-c)为由速度分量v_x、v_z求取应力分量σ_xx、σ_xz、σ_zz的公式，(1-d)、(1-e)为由应力分量σ_xx、σ_xz、σ_zz求取速度分量v_x、v_z的公式，这些变量的初始值均赋为0，步骤5中的震源参数与时间有关，随着时间迭代，通过这五个公式的迭代，得到整个区域内的速度和应力值，若设定时间到，则模拟过程结束，输出波场快照和地震记录。