CN105158797A - 一种基于实际地震资料的交错网格波动方程正演的方法 - Google Patents

Abstract

本发明公开了一种基于实际地震资料的交错网格波动方程正演的方法，采用实际地震数据的子波进行正演模拟，克服现有技术对于地震资料匹配性较差的缺陷；应用高阶差分格式以及利用多源信息对正演模型参数进行约束，提高正演精度，克服频散现象；紧密结合已知多源资料(地震、测井、地质等资料)，优化正演技术流程与地质建模流程，提前设置好观测系统，在正演过程中就不需要手动去修改观测系统参数，提高计算效率与计算精度。本发明采用基于实际地震、测井、地质资料，利用实际地震资料的子波，先叠后逐步迭代地质模型再叠前正演的方法对实际地质模型进行正演模拟，可以灵活地匹配各个地区的地震资料，具有较高的计算精度和计算效率。

Description

一种基于实际地震资料的交错网格波动方程正演的方法

技术领域

本发明属于勘探地震正演模拟技术领域，尤其涉及一种基于实际地震资料的交错网格波动方程正演的方法。

背景技术

地震学家们最早用有限差分法来解决波的传播问题最早要追溯到40多年前，大部分早期的研究地震波传播的有限差分方法基于二阶位移方程。Alterman和Karal(1968)最早将有限差分方法应用于均匀介质的地震波数值模拟。Boore(1972)将有限差分法用于模拟非均匀介质中地震波的传播。Alford(1974)研究了声波方程有限差分模拟的精确性。Kelly等人(1976)研究了用有限差分制作合成地震记录的方法。Mdariga(1976)提出了一种较为先进的交错网格有限差分方法,并首先将其应用于模拟弹性介质内圆形扩展破裂产生的波动。virieux(1984,1986)提出了在交错网格中用一阶速度一应力方程来代替二阶位移方程并分别对SH波和P-SV波进行了数值模拟。计算中采用计算节点位置周围半个网格长度的值来进行差分计算,而传统的规则网格二阶方程需要在计算节点周围一个网格长度的值来差分计算,因此在模型参数相同的情况下模拟精度有了很大提高。并且一阶速度—应力方程不需要对弹性参数进行求导,因此也更加适用于复杂介质的模拟。由于交错网格的这些优点,目前时间域有限差分地震波模拟中大多数都采用了交错网格的方法。针对Virieux所用的二阶差分精度,Levander(1988)采用了四阶空间算子模拟P-SV波地震记录,使得模拟精度得到进一步提高。Dablain(1986)在求解二阶声波方程时,采用Taylor级数展开并将差分节点对称的表达式求和发现可以消去除了二阶导数外直到任意阶差分误差内的所有阶的导数项,并且差分格式也非常简单。这一发现使得有限差分的精度得到极大提高,并从模拟结果中指出时间四阶,空间十阶的模拟精度与时间四阶频率域模拟方法相当。高阶差分模拟的不足是增加了一定的计算时间,稳定性也比低阶要差些,不过从整体计算精度和计算时间权衡来考虑,高阶差分的效率相对于低阶差分来说要高得多。Crase(1990)将高阶差分方法运用到求解二阶弹性波方程。Graves(1990)给出了三维速度一应力方程交错网格有限差分法弹性波传播的模拟方法。董良国(2000)用类似Dablain的方法得到了一阶导数的高阶差分近似,并将其应用到一阶速度—应力弹性波方程的模拟中。由于交错网格的弹性参数不在一个网格上定义,因此在计算时实际上是用不同位置的弹性参数,弹性参数之间位置不匹配的在不均匀介质中实际上是波场的一种各向异性的行为,当波场遇到强烈的密度差时,就会产生不稳定的问题。为避免这个问题,Saenger(2000)提出了一种新型的交错网格来模拟弹性介质具有强烈差异的弹性波传播。与Virieux提出的交错网格不同,Virieux提出的交错网格将一个位置的弹性参数,应力和位移速度分布在四个位置,可以看作在一个单元矩形的四个顶点上,Saenger提出的交错网格将单元矩形旋转了45度,并且将位移速度和密度放置于一个位置,拉梅系数和应力放置于一个位置,这样在分别计算位移速度及应力时都只用到一个位置上弹性参数。这样就可以稳定地计算具有强弹性性质差异的介质。现在通常把Virieux提出的交错网格形式称为标准交错网格(SSG),把Saenger提出的交错网格形式称为旋转交错网格(RSG)。为解决计算机模拟区域的有限截断问题，Engquist、Majda、Clayton(1977)提出了一种基于Pade序列的单程波近似的边界条件，其对于垂直入射的波可以很好地吸收，但对于大角度入射的波而言，仍会产生强反射；为了解决垂直入射角度限制的问题，Higdon(1986、1987、1991)提出一种多角度吸收的边界条件，其可以通过设置角度参数来吸收某几个特定方向的入射波。宋鹏、王修田(2008)将Higdon吸收边界提高到四阶，取得了较理想的吸收效果。LiuYang、SenM.K.(2010)提出混合吸收边界，通过设置过渡层并在其中计算单程波和双程波然后通过线性加权得到最终波场，LiuYang、SenM.K.(2011)将这种吸收边界应用到三维波场的模拟中，数值上取得了较好的效果。针对人工边界处理问题，Berenger(1994)针对计算电磁提出的完全匹配层(PML)边界处理方法，其是通过修改吸收层中的方程来使得其振动幅值可以迅速衰减达到吸收目的，Berenger(1996)将其应用于三维计算电磁中，Collino(2001)将PML引入弹性波的模拟中，取得了较好的吸收效果。地震子波是激发点激发地震波以后，传播到一定距离以后，波形逐渐稳定下来，具有多个相位、延续60-100ms的地震波，我们称之为地震子波。雷克子波是零相位子波，具有一个主瓣两个对称的旁瓣。雷克子波在地震正演、反演、处理和解释等方面广泛应用。由于受原始数据的频带和信噪比的限制，从预处理到偏移成像之间的每个处理步骤所得的结果都是具有带限子波的记录，以至最终的处理结果仍然是反射系数与带限子波的褶积。这种子波在一定程度上是可选的，处理中常用的子波是带通子波和Ricker子波，前者延续时间长，旁瓣波形复杂，后者旁瓣幅度比较大，都不是冲击函数的最好近似。为此，俞寿朋教授提出了一种新型子波，它是由不同宽度的Ricker子波台成的，它的主瓣窄、旁瓣幅度小、波形简单，在分辨率、保真度和信噪比方面都优于常用的两种子波。当采用理论子波进行正演时就会造成正演出的结果与实际地震剖面不匹配，因此影响我们对地震资料的认识。我们对采用雷克子波等理论子波正演所得到的正演记录进行处理时，由于所采用子波与实际子波相差较大，会造成一些错误的处理认识。进行叠前、叠后逆时偏移时，正演子波选取不当会造成偏移后的剖面与实际位置有较大的误差。正演子波越接近实际地震资料子波，相对应的偏移成像位置越准确。当我们进行全波形反演时，正演子波选取不当会造成反演出的参数同样有较大的误差。

目前有限差分波动方程正演的技术方案：

(1)基于二阶位移方程有限差分法(刘洋，李承楚，牟永光，任意偶数阶精度有限差分法数值模拟，1998年2月，石油地球物理勘探。)

技术方案：

①处理的对象：针对二维均匀复杂介质的地质模型；

②计算流程：定义参数，参数初始化；波场初始化，并在波场上加载震源(采用点震源)；计算模拟区域内的波场；计算边界区域的波场；波场更新；进行下一个时间循环直到设定的最大时间，输出数据。(说明：使用理论子波即点震源，进行的是二维均匀复杂介质的地质模型正演。时间空间参数、炮点位置需要手动调节。)

(2)交错网格有限差分法(周学明，交错网格高阶差分数值模拟及叠前逆时偏移，2010年6月，长安大学硕士论文。)

技术方案：

①处理的对象：针对二维非均匀复杂介质的地质模型；

②计算流程：定义参数，参数初始化；应力分量初始化，并在应力或速度分量加载震源(采用雷克子波)；计算模拟区域内的应力分量；计算边界区域的应力分量；计算模拟区域的速度分量；计算边界区域的速度分量；应力分量、速度分量更新；进行下一个时间循环直到设定的最大时间，输出数据。(说明：使用理论子波(雷克子波)，进行的是二维非均匀复杂介质的地质模型正演。时间空间参数、炮点位置需要手动调节。)

现有的有限差分波动方程正演存在对于地震资料匹配性较差，计算时间较长，计算效率较低。

发明内容

本发明的目的在于提供一种基于实际地震资料的交错网格波动方程正演的方法，旨在解决现有的有限差分波动方程正演存在对于地震资料匹配性较差，计算时间较长，计算效率较低的问题。

本发明是这样实现的，一种基于实际地震资料的交错网格波动方程正演的方法，其特征在于，所述基于实际地震资料的交错网格波动方程正演的方法

进一步，所述基于实际地震资料的交错网格波动方程正演的方法具体包括以下步骤：

步骤一，综合研究区已有地震资料获取时间采样率、空间采样间隔、覆盖次数为后面的程序参数设置提供参考；

步骤二，根据实际地震资料提取子波并保存，有井的地方可以联合地震和测井提取子波，没井的地方可以用地震数据做相关来提取子波；

步骤三，对工区已有声波、密度测井资料进行处理、整理得到工区相关地质模型的纵波速度、横波速度、密度参数；

步骤四，利用已得到的纵波速度、横波速度、密度和地质模型轮廓进行地质建模并保存；

步骤五，对地质模型进行叠后交错网格正演；

步骤六，对地质模型进行叠前交错网格正演。

进一步，所述步骤五具体包括：

第一步，根据稳定条件、频散条件结合原始地震资料设置网格大小、采样时间；

第二步，根据实际地震资料定义二维模型的大小，定义参数并初始化；

第三步，读取纵波速度、密度、子波，并计算边界区域纵波速度、密度及衰减系数；

第四步，计算模型的反射系数；

第五步，利用爆炸反射面原理，设置初始条件；

第六步，计算内部区域正应力分量；

第七步，计算边界区域正应力分量；

第八步，计算内部区域速度分量；

第九步，计算边界区域速度分量；

第十步，更新应力分量和速度分量，进行下一次时间循环重复第五步-第十步直到最大时间为止，输出叠后正演记录；

第十一步，对输出的叠后正演记录进行叠后逆时偏移，输出偏移记录；

第十二步，偏移记录和实际的叠后地震数据对比，调节地质模型重复第二步-第十步，直到达到期望的输出为止；结束叠后正演，并保存最后的地质模型参数和子波。

进一步，所述步骤六具体包括：

第一步，根据稳定条件、频散条件结合原始地震资料自动设置网格大小、采样时间；

第二步，根据实际地震资料定义二维模型的大小，定义参数并初始化；

第三步，读取保存地质模型参数和子波，并计算边界区域纵波速度、横波速度、密度，及衰减系数；

第四步，设置炮点的起始位置、根据跑间距，设置炮点循环变量；

第五步，设置应力变量速度变量的初始值为0；

第六步，根据设置的炮点位置，对相应网格点处对应力或速度给相应时间点的子波值；

第七步，根据设置的炮点位置，对相应网格点处对应力或速度给相应时间点的子波值；

第八步，计算内部区域正应力分量；

第九步，计算边界区域正应力分量；

第十步，计算内部区域速度分量；

第十一步，计算边界区域速度分量；

第十二步，更新应力分量和速度分量，进行下一次时间循环重复第六步-第十二步直到最大时间为止，输出相应炮点位置的叠前正演记录；

第十三步，进行下一炮循环重复第五步-第十二步进行下一炮叠前正演，直到最大炮点位置为止；

第十四步，对所有炮点正演记录进行常规处理、偏移并叠加形成最终剖面。

本发明的另一目的在于提供一种使用所述基于实际地震资料的交错网格波动方程正演的方法的油气地震勘探正演数值模拟系统。

本发明的另一目的在于提供一种使用所述基于实际地震资料的交错网格波动方程正演的方法的野外地震数据采集系统。

本发明提供的基于实际地震资料的交错网格波动方程正演的方法，是在MicrosoftVisualC++编译环境下用C语言实现，采用高阶交错网格有限差分法进行叠后、叠前正演并采用基于实际地震、测井、地质等资料，应用实际地震子波而不是理论子波，提出了一套先叠后逐步迭代地质模型再叠前正演的优化正演模型的方法对实际地质模型进行正演模拟，其可以灵活地匹配各个地区的地震资料，具有较高的计算精度和计算效率。本发明的核心是应用实际地震子波而不是理论子波，利用多源信息(测井、地震、地质等资料)对正演模型参数进行约束，提高正演精度，克服频散现象，提出了先叠后正演再叠前正演优化正演模型的方法。与现有技术相比，具有以下优势：

(1)采用实际地震数据的子波而不是理论子波进行正演模拟，克服现有技术对于实际地震资料匹配性较差的缺陷。

(2)应用高阶差分格式以及利用多源信息(测井、地震、地质等资料)对正演模型参数进行约束，提高正演精度，克服频散现象。

(3)紧密结合已知多源资料(地震、测井、地质等资料)，提出了一套先叠后逐步迭代地质模型再叠前正演的优化正演模型的方法，优化正演技术流程与地质建模流程，提前设置好观测系统，在正演过程中就不需要手动去修改观测系统参数，提高计算效率与计算精度。

附图说明

图1是本发明实施例提供的基于实际地震资料的交错网格波动方程正演的方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明提供一种基于实际地震资料的交错网格波动方程有限差分正演方法，它针对二维地质模型进行波动方程正演，可用于逆时偏移、全波形反演等方面，也可有效地指导油气地球物理勘探的采集、处理、解释。

下面结合附图对本发明的应用原理作详细的描述。

如图1所示，本发明实施例的基于实际地震资料的交错网格波动方程正演的方法包括以下步骤：

S101：综合研究区已有地震资料获取时间采样率、空间采样间隔、覆盖次数为后面的程序参数设置提供参考；

S102：根据实际地震资料提取子波并保存，有井的地方可以联合地震和测井提取子波，没井的地方可以用地震数据做相关来提取子波(杨培杰，地震子波提取方法综述，石油地球物理勘，2008)；

S103：对工区已有声波、密度测井资料进行处理、整理得到工区相关地质模型的纵波速度、横波速度、密度参数；

S104：利用已得到的纵波速度、横波速度、密度和地质模型轮廓进行地质建模并保存；

S105：对地质模型进行叠后交错网格正演(周学明，交错网格高阶差分数值模拟及叠前逆时偏移，2010年6月，长安大学硕士论文。)；

S106：对地质模型进行叠前交错网格正演(周学明，交错网格高阶差分数值模拟及叠前逆时偏移，2010年6月，长安大学硕士论文。)。

本发明的具体步骤包括：

(1)综合研究区已有地震资料获取时间采样率、空间采样间隔、覆盖次数为后面的程序参数设置提供参考。

(2)根据实际地震资料提取子波并保存，有井的地方可以联合地震和测井提取子波，没井的地方可以用地震数据做相关来提取子波。

(3)对工区已有井资料进行处理、整理得到工区相关地质模型的纵波速度、横波速度、密度参数；

(4)利用已得到的纵波速度、横波速度、密度和地质模型轮廓进行地质建模并保存。

(5)对地质模型进行叠后交错网格正演。

(5-1)根据稳定条件、频散条件结合原始地震资料设置网格大小、采样时间。

(5-2)根据实际地震资料定义二维模型的大小，定义参数并初始化；

(5-3)读取纵波速度、密度、子波，并计算边界区域纵波速度、密度及衰减系数；

(5-4)计算模型的反射系数；

(5-5)利用爆炸反射面原理，设置初始条件；

(5-6)计算内部区域正应力分量；

(5-7)计算边界区域正应力分量；

(5-8)计算内部区域速度分量；

(5-9)计算边界区域速度分量；

(5-10)更新应力分量和速度分量，进行下一次时间循环重复(5-5)-(5-10)直到最大时间为止，输出叠后正演记录；

(5-11)对输出的叠后正演记录进行叠后逆时偏移，输出偏移记录；

(5-12)偏移记录和实际的叠后地震数据对比，调节地质模型重复(5-2)-(5-10)，直到达到期望的输出为止；结束叠后正演，并保存最后的地质模型参数和子波。

(6)对地质模型进行叠前交错网格正演。

(6-1)根据稳定条件、频散条件结合原始地震资料自动设置网格大小、采样时间。

(6-2)根据实际地震资料定义二维模型的大小，定义参数并初始化；

(6-3)读取(5-12)保存地质模型参数和子波，并计算边界区域纵波速度、横波速度、密度，及衰减系数；

(6-4)设置炮点的起始位置、根据跑间距，设置炮点循环变量；

(6-5)设置应力变量速度变量的初始值为0；

(6-6)根据设置的炮点位置，对相应网格点处对应力或速度给相应时间点的子波值；

(6-7)根据设置的炮点位置，对相应网格点处对应力或速度给相应时间点的子波值；

(6-8)计算内部区域正应力分量；

(6-9)计算边界区域正应力分量；

(6-10)计算内部区域速度分量；

(6-11)计算边界区域速度分量；

(6-12)更新应力分量和速度分量，进行下一次时间循环重复(6-6)-(6-12)直到最大时间为止，输出相应炮点位置的叠前正演记录；

(6-13)进行下一炮循环重复(6-5)-(6-12)进行下一炮叠前正演，直到最大炮点位置为止；

(6-14)对所有炮点正演记录进行常规处理、偏移并叠加形成最终剖面。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于实际地震资料的交错网格波动方程正演的方法，其特征在于,所述基于实际地震资料的交错网格波动方程正演的方法采用基于实际地震、测井、地质资料，利用实际地震资料的子波，先叠后逐步迭代地质模型再叠前正演的方法对实际地质模型进行正演模拟。

2.如权利要求1所述的基于实际地震资料的交错网格波动方程正演的方法，其特征在于，所述基于实际地震资料的交错网格波动方程正演的方法具体包括以下步骤：

步骤一，综合研究区已有地震资料获取时间采样率、空间采样间隔、覆盖次数为后面的程序参数设置提供参考，程序参数通过读取地震数据获得；

步骤二，根据实际地震资料提取子波并保存，有井的地方联合地震和测井提取子波，没井的地方用地震数据做相关来提取子波

步骤三，对工区已有声波和密度测井数据进行统计分析，剔除异常数据得到工区相关地质模型的纵波速度、横波速度、密度参数，对于没有横波测井的地方根据测井数据进行横波反演得到；

步骤四，利用已得到的纵波速度、横波速度、密度和地质模型轮廓进行地质建模并保存；

步骤五，读取提取的子波，对地质模型进行叠后交错网格正演；

步骤六，对地质模型进行叠前交错网格正演。

3.如权利要求2所述的基于实际地震资料的交错网格波动方程正演的方法，其特征在于，所述步骤五具体包括：

第一步，根据稳定条件下式、频散条件结合原始地震资料设置网格空间步长、时间步长，频散条件最短波长达到6个网格以上：

$\mathrm{\Δ}t\sqrt{\frac{V_p^2}{{\mathrm{\Δx}}^2}+\frac{V_s^2}{{\mathrm{\Δz}}^2}}\≤1;$

其中，Δt，Δx，Δz为时间步长和x,z方向的步长，V_p，V_s为纵、横波速度；

第二步，根据实际地震资料定义二维模型的大小，定义参数并初始化；

第三步，读取地质模型参数和给定子波，计算边界区域纵波速度、横波速度、密度，及衰减系数；边界区域的纵波速度、横波速度、密度与离它最近的内部区域网格纵波速度、横波速度、密度相等，内部区域的衰减系数为0，边界区域的衰减系数根据下式计算：

d(x)＝(3Vmax/2PML)*log₁₀(1/R)*(x/PML)*(x/PML)

d(z)＝(3Vmax/2PML)*log₁₀(1/R)*(z/PML)*(z/PML)；

其中V_max为最大纵波速度，PML为匹配层宽度，x为横向距离，z为纵向距离，R为理想的边界反射系数，R取为0.000001，d(x)、d(z)不等于零时表示衰减，当d(x)、d(z)等于零时表示不衰减；

第四步，计算模型的反射系数，反射系数的计算公式根据下式计算：

R＝(r2*v2-r1*v1)/(r2*v2+r1*v1)；

其中R为反射系数，v1为第一层的速度，r1为第一层的密度，v2为第二层的速度，r2为第二层的密度；

第五步，利用爆炸反射面原理(牟永光，2007)，设置初始条件；

第六步，根据下式计算内部区域应力分量；其中弹性参数λ、μ根据纵波速度、横波速度、密度算出：

$\frac{\∂{\mathrm{\σ}}_{xx}}{\∂t}=\left(\mathrm{\λ}+2\mathrm{\μ}\right)\frac{\∂v_x}{\∂x}+\mathrm{\λ}\frac{\∂v_z}{\∂z};$

$\frac{\∂{\mathrm{\σ}}_{zz}}{\∂t}=\left(\mathrm{\λ}+2\mathrm{\μ}\right)\frac{\∂v_z}{\∂z}+\mathrm{\λ}\frac{\∂v_x}{\∂x}$

第七步，根据下面三个公式计算边界区域应力分量；其中弹性参数λ、μ根据纵波速度、横波速度、密度算出，d(x)、d(z)由公式 $\begin{array}{c}d\left(x\right)=(3Vmax/2PML)*{\log }_{10}(1/R)*(x/PML)*(x/PML)\\ d\left(z\right)=(3Vmax/2PML)*{\log }_{10}(1/R)*(z/PML)*(z/PML)\end{array}$ 算出：

$\begin{array}{c}{\mathrm{\σ}}_{zz}={\mathrm{\σ}}_{zz}^x+{\mathrm{\σ}}_{zz}^z\\ \frac{\∂{\mathrm{\σ}}_{zz}^x}{\∂t}+d\left(x\right){\mathrm{\σ}}_{zz}^x=\mathrm{\λ}\frac{\∂v_x}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{zz}^z}{\∂t}+d\left(z\right){\mathrm{\σ}}_{zz}^z=\left(\mathrm{\λ}+2\mathrm{\μ}\right)\frac{\∂v_z}{\∂z}\end{array};$

$\begin{array}{c}{\mathrm{\σ}}_{xx}={\mathrm{\σ}}_{xx}^x+{\mathrm{\σ}}_{xx}^z\\ \frac{\∂{\mathrm{\σ}}_{xx}^x}{\∂t}+d\left(x\right){\mathrm{\σ}}_{xx}^x=\left(\mathrm{\λ}+2\mathrm{\μ}\right)\frac{\∂v_x}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xx}^z}{\∂t}+d\left(z\right){\mathrm{\σ}}_{xx}^z=\mathrm{\λ}\frac{\∂v_z}{\∂z}\end{array};$

第八步，根据下式计算内部区域速度分量；

$\mathrm{\ρ}\frac{\∂v_z}{\∂t}=\frac{\∂{\mathrm{\σ}}_{xz}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}$

第九步，根据下式计算边界区域速度分量；

$\begin{array}{c}{v}_z=v_z^x+v_z^z\\ \frac{\∂v_z^x}{\∂t}+d\left(x\right)v_z^x=\frac{1}{\mathrm{\ρ}}\frac{\∂{\mathrm{\σ}}_{xz}}{\∂x}\\ \frac{\∂v_z^z}{\∂t}+d\left(z\right)v_z^z=\frac{1}{\mathrm{\ρ}}\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}\end{array};$

第十步，更新应力分量和速度分量，进行下一次时间循环重复第五步-第十步直到最大时间为止，输出叠后正演记录；

第十一步，对输出的叠后正演记录进行叠后逆时偏移，输出偏移记录；

第十二步，偏移记录和实际的叠后地震数据对比，调节地质模型重复第二步-第十步，直到达到期望的输出为止；结束叠后正演，并保存最后的地质模型参数和子波。

4.如权利要求2所述的基于实际地震资料的交错网格波动方程正演的方法，其特征在于，所述步骤六具体包括：

第一步，根据稳定条件下式、频散条件结合原始地震资料设置网格空间步长、时间步长，频散条件最短波长达到6个网格以上：

$\mathrm{\Δ}t\sqrt{\frac{V_p^2}{{\mathrm{\Δx}}^2}+\frac{V_s^2}{{\mathrm{\Δz}}^2}}\≤1;$

其中，Δt，Δx，Δz为时间步长和x,z方向的步长，V_p，V_s为纵、横波速度；

第二步，根据实际地震资料定义二维模型的大小，定义参数并初始化；

第三步，读取地质模型参数和给定子波，设置的炮点位置，计算边界区域纵波速度、横波速度、密度，及衰减系数；边界区域的纵波速度、横波速度、密度与离它最近的内部区域网格纵波速度、横波速度、密度相等，内部区域的衰减系数为0，边界区域的衰减系数根据下式计算：

d(x)＝(3Vmax/2PML)*log₁₀(1/R)*(x/PML)*(x/PML)

d(z)＝(3Vmax/2PML)*log₁₀(1/R)*(z/PML)*(z/PML)；

其中V_max为最大纵波速度，PML为匹配层宽度，x为横向距离，z为纵向距离，R为理想的边界反射系数，R取为0.000001，d(x)、d(z)不等于零时表示衰减，当d(x)、d(z)等于零时表示不衰减；

第四步，设置应力变量、速度变量的初始值为0；

第五步，根据设置的炮点位置，对相应网格点处对应力或速度给相应时间点的子波值；

第六步，根据下式计算内部区域应力分量；其中弹性参数λ、μ根据纵波速度、横波速度、密度算出：

$\frac{\∂{\mathrm{\σ}}_{xx}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_x}{\∂x}+\mathrm{\λ}\frac{\∂v_z}{\∂z}$

$\frac{\∂{\mathrm{\σ}}_{zz}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_z}{\∂z}+\mathrm{\λ}\frac{\∂v_x}{\∂x};$

$\frac{\∂{\mathrm{\σ}}_{xz}}{\∂t}=\mathrm{\μ}(\frac{\∂v_x}{\∂z}+\frac{\∂v_z}{\∂x})$

第七步，根据下面三个公式计算边界区域应力分量；其中弹性参数λ、μ根据纵波速度、横波速度、密度算出，d(x)、d(z)由公式 $\begin{array}{c}d\left(x\right)=(3Vmax/2PML)*{\log }_{10}(1/R)*(x/PML)*(x/PML)\\ d\left(z\right)=(3Vmax/2PML)*{\log }_{10}(1/R)*(z/PML)*(z/PML)\end{array}$ 算出：

$\begin{array}{c}{\mathrm{\σ}}_{zz}={\mathrm{\σ}}_{zz}^x+{\mathrm{\σ}}_{zz}^z\\ \frac{\∂{\mathrm{\σ}}_{zz}^x}{\∂t}+d\left(x\right){\mathrm{\σ}}_{zz}^x=\mathrm{\λ}\frac{\∂v_x}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{zz}^z}{\∂t}+d\left(z\right){\mathrm{\σ}}_{zz}^z=\left(\mathrm{\λ}+2\mathrm{\μ}\right)\frac{\∂v_z}{\∂z}\end{array};$

$\begin{array}{c}{\mathrm{\σ}}_{xx}={\mathrm{\σ}}_{xx}^x+{\mathrm{\σ}}_{xx}^z\\ \frac{\∂{\mathrm{\σ}}_{xx}^x}{\∂t}+d\left(x\right){\mathrm{\σ}}_{xx}^x=\left(\mathrm{\λ}+2\mathrm{\μ}\right)\frac{\∂v_x}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xx}^z}{\∂t}+d\left(z\right){\mathrm{\σ}}_{xx}^z=\mathrm{\λ}\frac{\∂v_z}{\∂z}\end{array};$

$\begin{array}{c}{\mathrm{\σ}}_{xz}={\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z\\ \frac{\∂{\mathrm{\σ}}_{xz}^x}{\∂t}+d\left(x\right){\mathrm{\σ}}_{xz}^x=\mathrm{\μ}\frac{\∂v_z}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xz}^z}{\∂t}+d\left(z\right){\mathrm{\σ}}_{xz}^z=\mathrm{\μ}\frac{\∂v_x}{\∂z}\end{array};$

第八步，根据下式计算内部区域速度分量；

$\begin{array}{c}\mathrm{\ρ}\frac{\∂v_x}{\∂t}=\frac{\∂{\mathrm{\σ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{xz}}{\∂z}\\ \mathrm{\ρ}\frac{\∂v_z}{\∂t}=\frac{\∂{\mathrm{\σ}}_{xz}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}\end{array};$

第九步，根据下式计算边界区域速度分量；

$\begin{array}{c}{v}_x=v_x^x+v_x^z\\ \frac{\∂v_x^x}{\∂t}+d\left(x\right)v_x^x=\frac{1}{\mathrm{\ρ}}\frac{\∂{\mathrm{\σ}}_{xx}}{\∂x}\\ \frac{\∂v_x^z}{\∂t}+d\left(z\right)v_x^z=\frac{1}{\mathrm{\ρ}}\frac{\∂{\mathrm{\σ}}_{xz}}{\∂z}\end{array};$

$\begin{array}{c}{v}_z=v_z^x+v_z^z\\ \frac{\∂v_z^x}{\∂t}+d\left(x\right)v_z^x=\frac{1}{\mathrm{\ρ}}\frac{\∂{\mathrm{\σ}}_{xz}}{\∂x}\\ \frac{\∂v_z^z}{\∂t}+d\left(z\right)v_z^z=\frac{1}{\mathrm{\ρ}}\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}\end{array};$

第十步，更新应力分量和速度分量，进行下一次时间循环重复第五步-第十步直到最大时间为止，输出相应炮点位置的叠前正演记录；

第十一步，进行下一炮循环重复第五步-第十步进行下一炮叠前正演，直到最大炮点位置为止；

第十二步，对所有炮点正演记录进行常规处理、偏移并叠加形成最终剖面。

5.一种使用权利要求1-4任意一项所述基于实际地震资料的交错网格波动方程正演的方法的油气地震勘探正演数值模拟系统。

6.一种使用权利要求1-4任意一项所述基于实际地震资料的交错网格波动方程正演的方法的野外地震数据采集系统。