CN105182414A - 一种基于波动方程正演去除直达波的方法 - Google Patents

Abstract

本发明公开了一种基于波动方程正演去除直达波的方法用波动方程正演去除直达波，去除地震反射记录中的直达波，消除直达波对有效反射信号干扰。本发明克服了现有技术需预先知道直达波的特征的缺陷；采用基于波动方程正演去除直达波的方法，不需过多的人工干预，几乎不会损失有效信号，不会产生额外的假象，参数容易确定。另外，本发明采用基于波动方程正演去除直达波的方法，过程简单，易实现，去除直达波彻底，提高了计算的效率和精度，能够有效、干净地去除地震反射记录中的直达波，消除直达波对有效反射信号干扰，并且不损失有效反射信号；可以应用在逆时偏移、全波形反演去除直达波中，也可应用到野外地震数据的去除直达波中。

Description

一种基于波动方程正演去除直达波的方法

技术领域

本发明属于地震资料数据处理技术领域，尤其涉及一种基于波动方程正演去除直达波的方法。

背景技术

在原始地震资料中，直达波与折射波能量较强，对于反射信号来讲是一种干扰波，如不加以处理会严重干扰浅层的反射信号；常规的处理方法是采用切除法来消除直达波与折射波的影响；但是在切除直达波和折射波的同时也切除掉了远偏移距的有效信号，往往会丢失大量有用的浅层地质构造反射信息，随着远偏移距的信息越来越被人们所重视,近年来越来越多地采用滤波法来衰减直达波和折射波，这样可以在衰减直达波和折射波的同时保护远偏移距的有效反射信号，无疑是一种较好的处理方法。根据地震波传播理论,直达波是沿地表传播的,它的线性时距函数与反射波的呈双曲线关系的时距函数存在着较大的差异。SVD理论在地震勘探领域应用广泛。李文杰等提出利用SVD滤波法衰减地震记录中的直达波和折射波而突出浅层反射波场。为利用直达波与折射波的线性特征，张兴岩等提出用τ-p变换衰减直达波。首先在τ-p域中将有效信号滤除，仅保留直达波和折射波信息，然后采用自适应减法从原始数据中减去直达波和折射波信息，这种方法能更有效地衰减掉直达波和折射波，保护远偏移距信号。地震学家们最早用有限差分法来解决波的传播问题要追溯到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)。在原始地震资料中，直达波与折射波传播距离短、能量较强，与反射波之间干涉，对于反射信号来讲是一种干扰波，如不加以处理会严重干扰浅层的反射信号的叠加效果。在进行叠前逆时偏移时，直达波在逆时偏移中必须切除,因为直达波与反射信息无关。从成像条件观点来看,直达波上的每点都满足成像时间,会给成像剖面上造成严重的干扰。同样直达波与反射信息无关，在进行全波形反演时，直达波能量较强，会对反演结果造成干扰，所以必须切除。

目前去除直达波的技术方案：

(1)切除法

拾取直达波，在直达波及以上的记录赋值给成0。(说明：方法简单，但是需人工交互拾取直达波，人工干预过多，效率很低，并且切除了大量远偏移距的有效信号，往往会丢失大量有用的浅层地质构造反射信息。)

(2)滤波法(李文杰，魏修成，刘洋，吴长江，张朝峰，SVD滤波法在直达波和折射波衰减处理中的应用。2004年10月，石油勘探与开发。)

技术方案：

①处理的对象：针对二维地震记录；

②计算流程：定义参数，参数初始化；读取数据，确定直达波滤波带；计算时窗长度；计算奇异值谱；对记录进行高通滤波。(说明：能在一定程度上改善地震记录质量，但是滤波参数较难确定。)

(3)数学变换去除法(张兴岩，方中于，史文英，但志伟，孙雷鸣，基于τ-p变换的直达波与折射波衰减方法研究及应用。西部探矿工程，2015年第1期。)

技术方案：

①处理的对象：针对二维时间域地震记录；

②计算流程：定义参数，参数初始化；读取数据；对时间域地震记录进行内切除，动校正；对内切除后的地震记录进行τ-p变换到τ-p域；对τ-p域数据进行直达波衰减；用τ-p反变换转到时间域；反动校正；自适应相减。(说明：能较好的衰减直达波，但是仍然使用了切除手段，在地震记录质量较差的情况下，效果不理想，并且参数选择不当会产生假象。)

现有的期初直达波方法存在需预先知道直达波的特征，过多的人工干预，容易损失有效信号，产生额外的假象，参数不容易确定，过程繁琐，去除直达波不彻底，计算的效率和精度较低。

发明内容

本发明的目的在于提供一种基于波动方程正演去除直达波的方法，旨在解决现有的期初直达波方法存在需预先知道直达波的特征，过多的人工干预，容易损失有效信号，产生额外的假象，参数不容易确定，过程繁琐，去除直达波不彻底，计算的效率和精度较低的问题。

本发明是这样实现的，一种基于波动方程正演去除直达波的方法，所述基于波动方程正演去除直达波的方法

进一步，所述基于波动方程正演去除直达波的方法具体包括以下步骤：

步骤一，处理、分析、整理得到相关地质模型的纵波速度V_p、横波速度V_s、密度参数DEN；

步骤二，利用已得到的纵波速度V_p、横波速度V_s、密度DEN建立地质模型；

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

步骤四，生成直达波记录；

步骤五，去除直达波。

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

第一步，根据稳定条件下式、频散条件结合原始地震资料设置网格空间步长、时间步长，频散条件最短波长达到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；

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

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

$\begin{array}{c}\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}\\ \frac{\∂{\mathrm{\σ}}_{xz}}{\∂t}=\mathrm{\μ}(\frac{\∂v_x}{\∂z}+\frac{\∂v_z}{\∂x})\end{array};$

第七步，根据下面三个公式计算边界区域应力分量；其中弹性参数λ、μ根据纵波速度、横波速度、密度算出，d(x)、d(z)由公式算出：

$\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};$

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

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

第一步，读取地质模型参数，将地质模型第一层网格处的纵波速度、横波速度、密度设置为整个地质模型纵波速度、横波速度、密度，计算边界区域纵波速度、横波速度、密度，衰减系数，计算方法为读取地质模型参数和给定子波，设置的炮点位置，计算边界区域纵波速度、横波速度、密度，及衰减系数；边界区域的纵波速度、横波速度、密度与离它最近的内部区域网格纵波速度、横波速度、密度相等。内部区域的衰减系数为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；

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

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

$\begin{array}{c}\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_z}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{zz}}{\∂t}=\mathrm{\μ}\left(\frac{\∂v_x}{\∂z}+\mathrm{\λ}\frac{\∂v_z}{\∂x}\right)\end{array};$

第七步，计算边界区域应力分量，根据下面三个公式计算边界区域应力分量；其中弹性参数λ、μ根据纵波速度、横波速度、密度算出，d(x)、d(z)由公式算出：

$\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};$

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

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

第一步，分别读取所生成的叠前正演记录；

第二步，用生成的叠前正演记录减去所生成的叠前正演记录即为去除直达波以后的记录。

本发明的另一目的在于提供一种所述的基于波动方程正演去除直达波的方法的系统，所述系统包括：

直达波计算模块，用于以二维地震数据和地质模型为对象，输入二维地震数据，采用交错网格有限差分的计算方法，计算并去除二维地震数据的直达波；

二维正演记录计算模块，用于根据地质模型采用交错网格有限差分的计算方法计算二维正演记录；

二维直达波记录模块，用于读取地质模型，将地质模型的第一层模型参数做为整个模型的参数，采用交错网格有限差分的计算方法计算二维直达波记录，用二维正演记录减去二维直达波记录得到最终的去除直达波的记录。

进一步，所述二维直达波记录模块包括：

读取单元，用于分别读取所生成的叠前正演记录；

计算单元，用生成的叠前正演记录减去所生成的叠前正演记录即为去除直达波以后的记录。

本发明的另一目的在于提供一种使用所述基于波动方程正演去除直达波方法的逆时偏移系统。

本发明的另一目的在于提供一种使用所述基于波动方程正演去除直达波方法的全波形反演系统。

本发明的另一目的在于提供一种使用所述基于波动方程正演去除直达波方法的野外地震数据系统。

本发明提供的基于波动方程正演去除直达波的方法，采用交错网格有限差分正演方法去除直达波，包含了交错网格有限差分正演方法方法，也包含了地质模型的建立，在计算过程中两次采用高阶交错网格有限差分正演方法分别得到了正演记录和直达波记录，具有较高的精度；本发明的核心计算公式是一阶速度-应力弹性波方程，计算效率高；采用高阶交错网格有限差分求解一阶速度-应力弹性波方程，计算精度高；采用两次采用交错网格有限差分正演方法分别得到了正演记录和直达波记录，计算精度、效率高。本发明针对二维正演记录或野外地震资料进行处理采用波动方程正演去除直达波，能够有效、干净地去除地震反射记录中的直达波，消除直达波对有效反射信号干扰，并且不损失有效反射信号；可以应用在逆时偏移、全波形反演去除直达波中，也可应用到野外地震数据的去除直达波中。进行计算流程的优化组合，同时考计算效率和计算精度。与现有技术相比具有以下优势：

(1)采用基于波动方程正演去除直达波的方法，克服现有技术需预先知道直达波的特征的缺陷。

(2)采用基于波动方程正演去除直达波的方法，不需过多的人工干预，几乎不会损失有效信号，不会产生额外的假象，参数容易确定。

(3)采用基于波动方程正演去除直达波的方法，原理、过程简单，易实现，能够有效、干净地去除地震反射记录中的直达波，消除直达波对有效反射信号干扰，并且几乎不损失有效反射信号，提高计算的效率和精度。

(4)本发明采用波动方程正演去除直达波，能够有效、干净地去除地震反射记录中的直达波，消除直达波对有效反射信号干扰，并且不损失有效反射信号；可以应用在逆时偏移、全波形反演去除直达波中，也可应用到野外地震数据的去除直达波中。

附图说明

图1是本发明实施例提供的基于波动方程正演去除直达波的方法流程图。

具体实施方式

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

本发明是要提供一种基于波动方程有限差分正演去除直达波方法，可用于逆时偏移、全波形反演去除直达波等方面，也可用于实际野外地震资料处理，有效地减少了直达波的影响。

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

如图1所示，本发明实施例的基于波动方程正演去除直达波的方法包括以下步骤：

S101：根据已有资料进行处理、分析、整理得到相关地质模型的纵波速度、横波速度、密度参数；

S102：利用已得到的纵波速度、横波速度、密度参数建立地质模型；

S103：对地质模型进行叠前交错网格正演；

S104：生成直达波记录；

S105：用生成的叠前正演记录减去所生成的叠前正演记录即为去除直达波以后的记录。

本发明的具体步骤包括：

(1)根据已有资料进行处理、分析、整理得到相关地质模型的纵波速度V_p、横波速度V_s、密度参数DEN(理论模型可以直接给出V_p、V_s、DEN)；

(2)利用已得到的纵波速度V_p、横波速度V_s、密度DEN建立地质模型。

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

(3-1)根据稳定条件(公式1)、频散条件(最短波长达到6个网格以上)结合原始地震资料设置网格空间步长、时间步长。

$\mathrm{\Δ}t\sqrt{\frac{V_p^2}{{\mathrm{\Δx}}^2}+\frac{V_s^2}{{\mathrm{\Δz}}^2}}\≤1---\left(1\right)$

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

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

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

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

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

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

(3-4)设置应力变量、速度变量的初始值为0；

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

(3-6)根据公式(3)计算内部区域应力分量；其中弹性参数λ、μ可以根据纵波速度、横波速度、密度算出。

$\begin{array}{c}\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}\\ \frac{\∂{\mathrm{\σ}}_{xz}}{\∂t}=\mathrm{\μ}\left(\frac{\∂v_x}{\∂z}+\mathrm{\λ}\frac{\∂v_z}{\∂x}\right)\end{array}---\left(3\right)$

(3-7)根据公式(4)-(6)计算边界区域应力分量；其中弹性参数λ、μ可以根据纵波速度、横波速度、密度算出。d(x)、d(z)由公式(2)算出。

$\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}---\left(4\right)$

$\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}---\left(5\right)$

$\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}---\left(6\right)$

(3-8)根据公式(7)计算内部区域速度分量；

$\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}---\left(7\right)$

(3-9)根据公式(8)-(9)计算边界区域速度分量；

$\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}---\left(8\right)$

$\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}---\left(9\right)$

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

(4)生成直达波记录。

(4-1)除了纵波速度、横波速度、密度这3个参数外，其他一切参数同(3)中的参数。

(4-2)读取地质模型参数。将地质模型第一层网格处的纵波速度、横波速度、密度设置为整个地质模型纵波速度、横波速度、密度(第一层网格处的参数可以根据测井资料给定，也可利用由静校正以后的直达波记录、速度分析、层析速度反演等方法等给定)，计算边界区域纵波速度、横波速度、密度，衰减系数。计算方法同(3-3)；

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

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

(4-6)计算内部区域应力分量，方法同(3-6)；

(4-7)计算边界区域应力分量，方法同(3-7)；

(4-8)计算内部区域速度分量，方法同(3-8)；

(4-9)计算边界区域速度分量，方法同(3-9)；

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

(5)去除直达波。

(5-1)分别读取(3-10)、(4-10)所生成的叠前正演记录。

(5-2)用(3-10)生成的叠前正演记录减去(4-10)所生成的叠前正演记录即为去除直达波以后的记录。

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

Claims (10)

1.一种基于波动方程正演去除直达波的方法，其特征在于，所述基于波动方程正演去除直达波的方法用波动方程正演去除直达波，去除地震反射记录中的直达波，消除直达波对有效反射信号干扰，具体包括以下步骤：

以二维地震数据和地质模型为对象，输入二维地震数据，采用交错网格有限差分的计算方法，计算并去除二维地震数据的直达波；

根据地质模型采用交错网格有限差分的计算方法计算二维正演记录；

读取地质模型，将地质模型的第一层模型参数做为整个模型的参数，采用交错网格有限差分的计算方法计算二维直达波记录；

用二维正演记录减去二维直达波记录得到最终的去除直达波的记录。

2.如权利要求1所述的基于波动方程正演去除直达波的方法，其特征在于，所述基于波动方程正演去除直达波的方法具体包括以下步骤：

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

步骤二，利用已得到的纵波速度V_p、横波速度V_s、密度DEN建立地质模型，并保存地质模型数据；

步骤三，读取地质模型数据，对地质模型进行叠前交错网格正演并保存正演记录；

步骤四，根据上述建立的地质模型，把第一层网格点处的纵波速度V_p、横波速度V_s、密度DEN作为整个模型的V_p、V_s、DEN，然后进行叠前交错网格正演生成直达波记录并保存正演记录；

步骤五，读取步骤三和步骤四的记录，并用步骤三的记录数据减去步骤四的数据，得到去除直达波的记录。

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)等于零时表示不衰减；

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

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

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

$\begin{array}{c}\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})\end{array};$

第七步，根据下面三个公式计算边界区域应力分量；其中弹性参数λ、μ根据纵波速度、横波速度、密度算出，d(x)、d(z)由公式

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

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

$\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=(\mathrm{\λ}+2\mathrm{\μ})\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=(\mathrm{\λ}+2\mathrm{\μ})\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(z\right)v_x^x=\frac{1}{\mathrm{\ρ}}\frac{\∂{\mathrm{\σ}}_{xx}}{\∂z}\\ \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};$

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

4.如权利要求2所述的基于波动方程正演去除直达波的方法，其特征在于，所述步骤四具体包括：

第一步，读取地质模型参数，将地质模型第一层网格处的纵波速度、横波速度、密度设置为整个地质模型纵波速度、横波速度、密度，计算边界区域纵波速度、横波速度、密度，衰减系数，计算方法为读取地质模型参数和给定子波，设置的炮点位置，计算边界区域纵波速度、横波速度、密度，及衰减系数；边界区域的纵波速度、横波速度、密度与离它最近的内部区域网格纵波速度、横波速度、密度相等，内部区域的衰减系数为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；

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

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

$\begin{array}{c}\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})\end{array};$

第七步，计算边界区域应力分量，根据下面三个公式计算边界区域应力分量；其中弹性参数λ、μ根据纵波速度、横波速度、密度算出，d(x)、

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

d(z)由公式d(z)＝(3Vmax/2PML)*log₁₀(1/R)*(z/PML)*(z/PML)算出：

$\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=(\mathrm{\λ}+2\mathrm{\μ})\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=(\mathrm{\λ}+2\mathrm{\μ})\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(z\right)v_x^x=\frac{1}{\mathrm{\ρ}}\frac{\∂{\mathrm{\σ}}_{xx}}{\∂z}\\ \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.如权利要求2所述的基于波动方程正演去除直达波的方法，其特征在于，所述步骤五具体包括：

第一步，分别读取所生成的叠前正演记录；

第二步，用生成的叠前正演记录减去所生成的叠前正演记录即为去除直达波以后的记录。

6.一种如权利要求1-5任意一项所述的基于波动方程正演去除直达波的方法的系统，其特征在于，所述系统包括：

直达波计算模块，用于以二维地震数据和地质模型为对象，输入二维地震数据，采用交错网格有限差分的计算方法，计算并去除二维地震数据的直达波；

二维正演记录计算模块，用于根据地质模型采用交错网格有限差分的计算方法计算二维正演记录；

二维直达波记录模块，用于读取地质模型，将地质模型的第一层模型参数做为整个模型的参数，采用交错网格有限差分的计算方法计算二维直达波记录，用二维正演记录减去二维直达波记录得到最终的去除直达波的记录。

7.如权利要求6所述的系统，其特征在于，所述二维直达波记录模块包括：

读取单元，用于分别读取所生成的叠前正演记录；

计算单元，用生成的叠前正演记录减去所生成的叠前正演记录即为去除直达波以后的记录。

8.一种使用权利要求1-5任意一项所述基于波动方程正演去除直达波方法的逆时偏移系统。

9.一种使用权利要求1-5任意一项所述基于波动方程正演去除直达波方法的全波形反演系统。

10.一种使用权利要求1-5任意一项所述基于波动方程正演去除直达波方法的野外地震数据系统。