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

摘要

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

说明书

技术领域

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

背景技术

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

上世纪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_ztanθ，在该区间内为噪音，将该区间定义为噪音区；在k_x＜k_ztanθ时，有效 全波形反演梯度占据主导地位，将该区间定义为有效区；在有效区与噪音区之间 建立过渡区，区角度范围为

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

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

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

$h(k_x,k_z)=\left(\begin{array}{cc}1& k_x\le k_z\tan \theta \\ 0& k_x>k_z\tan \theta \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为水平波数为零时滤波算子 值。

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

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

其中。

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

$g_s(k_x,k_z)=\int \int g_s(x,z)e^{-j{k}_{x}x}{e}^{-j{k}_{z}z}\mathrm{dxdz}$

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

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

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

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

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

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

$g_s^{\prime }(k_x,k_z)=h(k_x,k_z){\int \int g}_s(x,z)e^{-j{k}_{x}x}{e}^{-j{k}_{z}z}\mathrm{dxdz}---\left(10\right)$

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

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

$g_s^{\prime }(x,z)=\frac{1}{4{\pi }^2}\int \int h(k_x,k_z){\int \int g}_s(x,z)e^{-j{k}_{x}x}{e}^{-j{k}_{z}z}\mathrm{dxdz}{e}^{j{k}_{x}x}{e}^{j{k}_{z}z}d{k}_{x}d{k}_z---\left(11\right)$

式中：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(k_x,k_z)=\int \int g_s(x,z)e^{-j{k}_{x}x}{e}^{-j{k}_{z}z}\mathrm{dxdz}$

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

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

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

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

二维声波波动方程为：

$\frac{1}{v^2}\frac{{\partial }^{2}u}{\partial t^2}=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial z^2}---\left(1\right)$

转换到频率域中：

$\frac{{\omega }^2}{v^2}u\%=k_x^{2}u\%+k_z^{2}u\%\Rightarrow \frac{{\omega }^2}{v^2}=k_x^2+k_z^2---\left(2\right)$

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

在均匀介质中，当频率一定时，该方程为以水平波数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)=\left(\begin{array}{cc}1& k_x\le k_z\tan \theta \\ 0& k_x>k_z\tan \theta \end{array}\right)---\left(3\right)$

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

根据线性拟合函数：

形成过渡区的滤波算子。

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

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

$g_s^{\prime }(k_x,k_z)=h(k_x,k_z)\int \int g_s(x,z)e^{-j{k}_{x}x}{e}^{-j{k}_{z}z}\mathrm{dxdz}---\left(10\right)$

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

$g_s^{\prime }(x,z)=\frac{1}{4{\pi }^2}\int \int h(k_x,k_z){\int \int g}_s(x,z)e^{-j{k}_{x}x}{e}^{-j{k}_{z}z}\mathrm{dxdz}{e}^{j{k}_{x}x}{e}^{j{k}_{z}z}d{k}_{x}d{k}_z---\left(11\right).$

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

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

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

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

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

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

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

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

步骤八，进行速度更新。

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

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

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

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

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