CN104570090B - 全波形反演噪音滤波算子的提取及使用其噪音滤波的方法 - Google Patents

Abstract

本发明公开了全波形反演噪音滤波算子的提取及使用其噪音滤波的方法。本发明在全波形反演中获取单炮的全波形反演梯度之后，对二维梯度进行二维傅里叶变换将数据从空间域转化到波数域；针对全波形反演速度噪音在波数域中的可存在区间，以及有效全波形反演速度在波数域中的分布，确立波数域中垂直波数与水平波数的关系；采用全波形反演噪音滤波算子避免数据在空间域与波数域之间转换产生的吉普斯现象；最后利用二维傅里叶反变换得到适应于地质特征的具有高信噪比、高保真度的全波形反演结果。本发明能够结合地质构造特征确定全波形反演噪声存在区间，保持有效全波形反演速度的同时使反演噪音得以消除，效果稳定可靠。

Description

全波形反演噪音滤波算子的提取及使用其噪音滤波的方法

技术领域

本发明属于地震反演领域，具体涉及全波形反演噪音滤波算子的提取及使用其噪音滤波的方法。

背景技术

速度在地震资料处理及解释中起着关键性的作用。地震波的走时信息仅考虑了地震波传播中短波长信息，走时反演仅能得到具有简单几何构造地层的模糊宏观速度模型或背景速度模型。走时反演不能反映有效中短波长地震波信息，无法提供高分辨率准确速度模型信息。借助于测井资料高纵向分辨率速度参数反演已发展成熟。具备高纵向和横向分辨率信息的速度构建则需求助于全波形反演。全波形反演是一种基于全波波动方程，在数据域开展的基于模拟数据与观测数据误差泛函最小的非线性反演方法。全波形反演利用的是整个或部分地震记录波形，此波形不仅包含地震波记录震相的运动学特征(如震相的到时、波速等)，还包含地震波动力学特征(如震相的振动方向、振幅、周期、相位等)，通过迭代反演实现速度参数构建。

上世纪80年代Tarantola等人提出了基于广义最小二乘反演理论的时间域全波形反演方法，对近20多年多维地震反演理论的发展产生了深远的影响。该方法使用理论波场与实际波场误差的L2泛函作为反演的目标函数，对模型进行网格离散化并利用双程波动方程数值模拟方法(如有限差分法、伪谱法等)模拟波场传播过程，通过梯度寻优实现模型参数迭代更新。该方法可适用于井间地震、VSP、广角地震、反射地震等多种观测系统，反演时能够充分利用多种波形信息，而不局限于只利用反射波信息。为提高计算效率，80年代末90年代初Pratt等人将全波形反演理论推广到频率域，形成了频率域全波形反演方法，也称波形层析成像方法。为了避免反演陷入局部极小，全波形反演方法对初始模型精度要求严格。针对地震数据频带宽度有限、反演初始速度模型获取困难等问题，Shin提出利用阻尼波场零频分量反演低频模型作为频率域波形反演的初始模型，即Laplace域全波形反演方法，为全波形反演理论与应用研究注入了新的活力。近几年来国内外全波形反演方法应用研究发展迅速，例如Bunks et al.(1995)研究了时间域内的多尺度波形反演。Pratt et al.(1990)发展了全波形反演成像的并进行油气开发CO₂流体注入前面破裂的监测应用，Smithyman et al.(2009)使用波形层析成像进行了近地表探测，Operto et al.(2006)进行全波形反演处理多次覆盖的海洋数据得到地壳尺度的结构。这些研究为区域深部构造及演化分析、浅表层环境调查、宏观速度场建模与成像、岩性参数反演提供了新的有力手段。目前全波形反演在理论研究方面已经成熟，取得理想反演效果。理论模型研究方面，由于正演地震数据基本上是完备的:观测系统规则，覆盖次数一致，地震波连续性强，反演结果稳定可靠。

在实际资料应用中，海上地震资料连续性相对保持较好。对于陆上地震资料采集观测系统由于受到建筑、河流以及复杂地形的影响，地震数据往往存在强烈的非规则性。全波形反演是以正演波场与观测波场的误差泛函作为目标函数，全波形反演的结果严重依赖于实际资料的品质，特别是数据的非完备性在反演结果中引入反演噪音，特别是对于反射波全波形反演严重影响了实际资料的应用。目前在反射波全波形反演噪音的引入还缺乏响应的理论研究和解决方法。

发明内容

为了解决现有的全波形反演算法中噪音的问题，本发明提供了全波形反演噪音滤波方法，并提供能够适应于地质构造特征的全波形反演噪音滤波算子的提取方法，解决由于数据采集不规则、地震数据覆盖次数低，地震数据缺少道集，以及地震数据缺乏连续性等相关的数据不完备性问题在全波形反演过程中引入的全波形反演结果中存在噪音问题。该发明不仅能够消除噪音的影响，又能够不破坏有效信息，算法稳定。

全波形反演噪音滤波算子的提取方法，包括如下步骤：

1)在全波形反演中获取单炮的全波形反演梯度之后，通过二维傅里叶变换得到多尺度二维波数下全波形反演参数的波数域分布特征，从而将数据从空间域转化到波数域；

2)针对全波形反演速度噪音在波数域中的可存在区间，以及有效全波形反演速度在波数域中的分布，确立波数域中垂直波数与水平波数的关系；

3)将地质构造的角度引入了虚假全波形反演信息，当地质构造倾角为θ时，在波数域中满足：tanθ＝k_x/k_z，因此当角度大于地质构造倾角θ时，存在关系式k_x＞k_z tanθ，在该区间内为噪音，将该区间定义为噪音区；在k_x＜k_z tanθ时，有效全波形反演梯度占据主导地位，将该区间定义为有效区；在有效区与噪音区之间建立过渡区，区角度范围为

其中，θ为地层倾角，k_x为水平波数，k_z为垂直波数；为微小扰动角度；

4)建立全波形反演噪音滤波算子

在理想情况下期望大于地质地层角度的噪音完全消除，同时小于地层倾角的有效信息完全保留，该种情况下波数域滤波算子为：

式中：h(k_x,k_z)为波数域滤波算子，θ为地层倾角，k_x为水平波数，k_z为垂直波数；

h(k_x,k_z)一阶导数为非连续函数，傅里叶变换存在严重的吉普斯现象，引入从噪音区到有效区的过渡区域，过渡区角度范围为在该角度区间的断点处波数域滤波算子为:

式中：θ为地层倾角，为微小扰动角度，k_x为水平波数，k_z为垂直波数；

根据线性拟合函数：

式中：G为滤波算子值随水平波数变化的斜率。P为水平波数为零时滤波算子值。

则全波形反演噪音滤波算子可表示为:

引入扰动量ε增强除法运算的稳定性，综上，可得稳定的全波形反演噪音滤波算子：

其中

进一步，步骤1)中单炮的全波形反演梯度采用单炮域全波形反演波场误差泛函梯度，对其进行空间域到波数域转换，具体方法为：

式中：g_s(x,z)是空间域单炮全波形反演波场误差泛函梯度，g_s(k_x,k_z)是空间域单炮全波形反演波场误差泛函梯度，k_x是水平波数，k_z是垂直波数；通过二维傅里叶变换将梯度从空间域转换到波数域。

本发明全波形反演噪音滤波算子对全波形反演中由于地震数据的不完备性以及反演本身存在的多解性的问题而引入的噪音具有良好的压制功能。

全波形反演噪音滤波方法，该方法使用上述滤波算子，包括下列步骤：

1)将滤波算子作用于全波形反演梯度数据，消除数据在空间域与波数域之间转换产生的吉普斯现象并获取有效全波形反演结果在波数域中表征；

2)通过傅立叶反变换将全波形反演梯度从波数域转换到空间域，获取速度更新量并进行速度更新后，得到适应于地质特征的具有高信噪比、高保真度的全波形反演结果。

进一步，将滤波算子作用于全波形反演梯度数据采用下式来完成:

式中：g_s(x,z)是空间域单炮全波形反演波场误差泛函梯度，g_s′(k_x,k_z)是实施滤波算子后波数域单炮全波形反演波场误差泛函梯度，k_x是水平波数，k_z是垂直波数，x为水平方向坐标，z为垂直方向坐标，j表示为虚部项。

进一步，将全波形反演梯度从波数域转换到空间域采用下式来完成：

式中：g_s(x,z)是空间域单炮全波形反演波场误差泛函梯度，g_s′(x,z)是实施滤波算子后空间域单炮全波形反演波场误差泛函梯度，k_x是水平波数，k_z是垂直波数，x为水平方向坐标，z为垂直方向坐标，j表示为虚部项。

进一步，所述全波形反演采用的资料为叠前炮域地震资料。

进一步，所述噪音为由于实际资料数据不完备所引起的噪音；所述的资料数据不完备是指地震数据采集不规则，地震数据覆盖次数低，地震数据缺少道集，缺乏连续性。

本发明在有效信息与无效信息之间建立过渡区，并采用平滑滤波因子，避免数据在空间域与波数域之间转换产生的吉普斯现象，同时获取有效全波形反演结果在波数域中表征；利用二维傅里叶反变换得到适应于地质特征的具有高信噪比、高保真度的全波形反演结果。本发明能够结合地质构造特征确定全波形反演噪声存在区间，保持有效全波形反演速度的同时使反演噪音得以消除，效果稳定可靠。

全波形反演噪音滤波算子的功能是消除由于数据不完备性(包含地震数据采集不规则，地震数据覆盖次数低，地震数据缺少道集，缺乏连续性)所引起的在全波形反演过程中产生的噪音。其优势是在不破坏有效反演参数信息的前提下消除全波形反演噪音的影响，提高反演结果的可靠性，建立具有高信噪比、高保真度的全波形反演结果。本发明能够适应于实际资料，推动全波形反演在实际生产的中应用。

附图说明

图1是滤波算子示意图。

图2是波数为5的波数域角度分布图。

图3是全波形反演噪音滤波算子计算及实施流程图。

图4是用于全波形反演的原始观测数据图。

图5原始观测数据得到的全波形反演梯度图。

图6是用全波形反演原始不完备的观测数据图。

图7是不完备的观测数据计算得到的全波形反演梯度图。

图8是全波形反演噪音滤波算子在单炮模型数据中的应用效果图(消除45度倾角以上噪音)。

图9是全波形反演噪音滤波算子在单炮模型数据中的应用效果图(消除15度以上倾角噪音)。

图10是全波形反演在实际资料中的应用效果图(含有噪音的全波形反演速度)。

图11是全波形反演噪音滤波算子在实际资料中的应用效果图(实施全波形反演噪音滤波算子的反演速度)。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合所附图，作详细说明如下。

1)对单炮域全波形反演波场误差泛函梯度，进行空间域到波数域转换：

式中：g_s(x,z)是空间域单炮全波形反演波场误差泛函梯度，g_s(k_x,k_z)是空间域单炮全波形反演波场误差泛函梯度，k_x是水平波数，k_z是垂直波数。通过二维傅里叶变换将梯度从空间域转换到波数域

2)根据研究工区已知的地质构造特征，判定全波形反演噪音滤波算子实施的最大地层倾角。

3)设定过渡区在波数域的区间

4)计算过渡区滤波算子：

二维声波波动方程为：

转换到频率域中：

该式即为波动方程精确频散关系方程。

在均匀介质中，当频率一定时，该方程为以水平波数k_x和垂直波数k_z为正交坐标系，以波数为半径的圆。根据空间域与波数域中角度对应关系可知，该式中包含了角度范围为0～2π的所有信息。对于地面地震角度范围缩小为0～π，该角度范围与地质构造特征无关。因此在大于地质构造的角度引入了虚假全波形反演信息。地质构造倾角为θ时，在波数域中满足：tanθ＝k_x/k_z，因此当角度大于地质构造倾角θ时，存在关系式k_x＞k_ztanθ，在该区间内为噪音。在k_x＜k_ztanθ时，有效全波形反演梯度占据主导地位。在理想状态下的波数域滤波算子为：

h(k_x,k_z)函数具有强非连续性，导致数据在空间域与波数域之间转换时产生吉普斯现象，因此需要引入从噪音区到有效区的过渡区域。假设过渡区角度范围为则需要满足以下条件:

根据线性拟合函数：

形成过渡区的滤波算子。

5)形成稳定的全波形反演噪音滤波算子:

6)将滤波算子作用于全波形反演梯度数据:

7)将全波形反演梯度从波数域转换到空间域

如图3所示，全波形反演噪音滤波算子计算及实施流程是：

步骤一：从初始速度场出发，获得地震反演波场；

步骤二：根据观测测场，获得误差传播波场；

步骤三:根据上述数据，获得空间域的全波形反演梯度；

步骤四：通过傅里叶变换，将数据从空间域转化到波数域；

步骤五：建立滤波算子，进行全波形反演，获得有效全波形反演梯度；

步骤六：根据傅里叶反变换的对称性，进行傅里叶反变换，将数据从波波数域转换到空间域；

步骤七：获取速度更新量；

步骤八，进行速度更新。

图1是滤波算子示意图。其中：①有效区，②过渡区，③噪音区。噪音区位于高角度范围内，有效信息则分布于低角度范围内。

图2是波数为5的波数域角度分布图。横坐标为水平波数，纵坐标为垂直波数。颜色深浅程度代表角度值的大小。

试验例1。本发明全波形反演噪音滤波算子及其滤波方法，应用于油田某块区的实际地震数据，取得了理想的成像效果。

对比图5与图7，由于地震观测数据的非连续性导致全波形反演梯度存在严重的噪音。对比图7与图9，全波形反演滤波算子的有效实施使得全波反演梯度在有效地信息不被破坏的同时高角度噪音得以消除，效果明显。

由图10和图11对比可以看出，经过全波形反演噪音滤波算子处理后，全波形反演结果得到明显改善，增强了反演结果的连续性，提高了全波形反演结果的信噪比和保真度，同时浅层的全波形反演结果得以恢复。

Claims (7)

1.全波形反演噪音滤波算子的提取方法，其特征在于包括如下步骤：

1)在全波形反演中获取单炮的全波形反演梯度之后，通过二维傅里叶变换得到多尺度二维波数下全波形反演参数的波数域分布特征，从而将数据从空间域转化到波数域；

2)针对全波形反演速度噪音在波数域中的存在区间，以及有效全波形反演速度在波数域中的分布，确立波数域中垂直波数与水平波数的关系；

3)在大于地层倾角时引入了虚假全波形反演信息，当地层倾角为θ时，在波数域中满足：tanθ＝k_x/k_z，在k_x＞k_z tanθ时，噪音占主导地位，将相应区间定义为噪音区；在k_x＜k_ztanθ时，有效全波形反演梯度占据主导地位，将相应区域定义为有效区；在有效区与噪音区之间建立过渡区，过渡区角度范围为

其中，θ为地层倾角，k_x为水平波数，k_z为垂直波数；为微小扰动角度；

4)建立全波形反演噪音滤波算子

在理想情况下期望大于地层倾角的噪音完全消除，同时小于地层倾角的有效信息完全保留，该种情况下波数域噪音滤波算子为：

$h(k_x,k_z)=\left\{\begin{array}{cc}1& k_x\≤k_{z}tan\mathrm{\θ}\\ 0& k_x>k_{z}tan\mathrm{\θ}\end{array}\right.---\left(1\right)$

式中：h(k_x,k_z)为波数域噪音滤波算子，θ为地层倾角，k_x为水平波数，k_z为垂直波数；

h(k_x,k_z)一阶导数为非连续函数，傅里叶变换存在严重的吉普斯现象，引入从噪音区到有效区的过渡区域，过渡区角度范围为在该过渡区角度范围的断点处波数域噪音滤波算子为:

式中：θ为地层倾角，为微小扰动角度，k_x为水平波数，k_z为垂直波数；

根据线性拟合函数：

式中：G为噪音滤波算子值随水平波数变化的斜率，P为水平波数为零时噪音滤波算子值；

则全波形反演噪音滤波算子可表示为:

其中

引入扰动量ε，综上，可得稳定的全波形反演噪音滤波算子：

其中

2.如权利要求1所述的全波形反演噪音滤波算子的提取方法，其特征在于：

步骤1)中单炮的全波形反演梯度采用空间域全波形反演波场误差泛函梯度，对其进行空间域到波数域转换，具体方法为：

$g_s(k_x,k_z)=\∫\∫g_s(x,z)e^{-{\mathrm{jk}}_{x}x}{e}^{-{\mathrm{jk}}_{z}z}dxdz---\left(8\right)$

式中：g_s(x,z)是空间域单炮全波形反演波场误差泛函梯度，x为水平方向坐标，z为垂直方向坐标，j表示为虚部项，g_s(k_x,k_z)是波数域单炮全波形反演波场误差泛函梯度，k_x是水平波数，k_z是垂直波数；通过二维傅里叶变换将单炮全波形反演波场误差泛函梯度从空间域转换到波数域。

3.全波形反演噪音滤波方法，使用权利要求1或2任意一权利要求所述的噪音滤波算子，其特征在于：所述方法包括下列步骤：

1)将噪音滤波算子作用于波数域单炮全波形反演波场误差泛函梯度,消除数据在空间域与波数域之间转换产生的吉普斯现象并获取有效全波形反演结果在波数域中表征；

2)通过傅立叶反变换将单炮全波形反演波场误差泛函梯度从波数域转换到空间域，获取速度更新量并进行速度更新后，得到适应于地质特征的具有高信噪比、高保真度的全波形反演结果。

4.如权利要求3所述的全波形反演噪音滤波方法，其特征在于：将噪音滤波算子作用于波数域单炮全波形反演波场误差泛函梯度,采用下式来完成:

${g_s}^{\′}\left(k_x,k_z\right)=h\left(k_x,k_z\right)\∫\∫g_s\left(x,z\right)e^{-{\mathrm{jk}}_{x}x}{e}^{-{\mathrm{jk}}_{z}z}dxdz---\left(9\right)$

式中：g_s(x，z)是空间域单炮全波形反演波场误差泛函梯度，g_s′(k_x,k_z)是实施噪音滤波算子后波数域单炮全波形反演波场误差泛函梯度，k_x是水平波数，k_s是垂直波数，x为水平方向坐标，z为垂直方向坐标，j表示为虚部项。

5.如权利要求3所述的全波形反演噪音滤波方法，其特征在于：将全波形反演波数域单炮全波形反演波场误差泛函梯度从波数域转换到空间域采用下式：

$g_s^{\′}(x,z)=\frac{1}{4{\mathrm{\π}}^2}\∫\∫h(k_x,k_z)\∫\∫g_s(x,z)e^{-{\mathrm{jk}}_{x}x}{e}^{-{\mathrm{jk}}_{z}z}{\mathrm{dxdze}}^{{\mathrm{jk}}_{x}x}{e}^{{\mathrm{jk}}_{z}z}{\mathrm{dk}}_x{\mathrm{dk}}_z---\left(10\right)$

式中：g_s(x，z)是空间域单炮全波形反演波场误差泛函梯度，g′_s(x,z)是实施噪音滤波算子后波数域单炮全波形反演波场误差泛函梯度，k_x是水平波数，k_s是垂直波数，x为水平方向坐标，z为垂直方向坐标，j表示为虚部项。

6.如权利要求3到5任意一权利要求所述的全波形反演噪音滤波方法，其特征在于：所述全波形反演采用的资料为叠前炮域地震资料。

7.如权利要求3到5任意一权利要求所述的全波形反演噪音滤波方法，其特征在于：所述噪音为由于资料数据不完备所引起的噪音；所述的资料数据不完备是指地震数据采集不规则，地震数据覆盖次数低，地震数据缺少道集，缺乏连续性。