CN105445798A - 一种基于梯度处理的全波形反演方法和系统 - Google Patents

Abstract

本发明公开了一种基于梯度处理的全波形反演方法和系统。该方法包括：残差计算步骤，由初始速度模型正演得到的计算数据和野外观测数据确定数据残差；梯度计算步骤，根据数据残差在初始速度模型中的反传波场以及正演得到的正传波场求取梯度场；扰动模型计算步骤，对梯度场进行滤波得到保留地层边界信息的速度扰动模型；更新步骤，根据速度扰动模型更新初始速度模型。本发明对全波形反演中的梯度进行低通和自适应两次滤波，通过低通滤波压制随机噪声，并通过自适应滤波恢复在低通滤波过程中损失的边界信息，达到边界保真的效果。

Description

一种基于梯度处理的全波形反演方法和系统

技术领域

本发明涉及地质勘探技术领域，具体地说，涉及一种基于梯度处理的全波形反演方法和系统。

背景技术

近年来全波反演理论在地质勘探领域得到广泛关注。全波反演可以提高常规地震分辨率并获得优化数据，提高对地质资源的评价能力，从而确定可采区并提出有利的井位建议。由于地质勘探所针对的目标不仅地表地质条件复杂，并且地下地质构造也复杂，因此，提高反演精度进而获得定量的储层表征参数非常重要。

现有技术中陆地地震资料的品质较差，噪音干扰严重。噪音的存在使得求取过程中的梯度存在较大误差，导致模型更新过程中加入了错误信息，从而使反演陷入局部极值，导致收敛失败。此外，噪音的存在也严重影响反演的精度。

目前对于对于噪音的压制，多采用目标泛函约束方法，加入正则化及先验模型信息等方法实现。但是现有技术中这些方法是在数据域计算，将误差引入目标泛函，导致计算量较大，计算过程繁琐。并且，在实际生产中很难得到准确的先验信息。

因此，亟需一种简单高效的全波形反演梯度处理方法和系统。

发明内容

针对上述技术问题，本发明提供一种基于梯度处理的全波形反演方法，包括以下步骤：

残差计算步骤，由初始速度模型正演得到的计算数据和野外观测数据确定数据残差；

梯度计算步骤，根据数据残差在初始速度模型中的反传波场以及正演得到的正传波场求取梯度场；

扰动模型计算步骤，对梯度场进行滤波得到保留地层边界信息的速度扰动模型；

更新步骤，根据速度扰动模型更新初始速度模型。

根据本发明的一个实施例，还包括

判断步骤，判断更新后的初始速度模型是否满足预设的精度要求，若不满足精度要求，则迭代执行残差计算步骤、梯度计算步骤、扰动模型计算步骤和更新步骤；若满足精度要求，则将更新后的初始速度模型确定为反演模型。

根据本发明的一个实施例，所述扰动模型计算步骤包括：

低通滤波步骤，对梯度场进行低通滤波以消除数据残差中的随机噪音；

自适应滤波步骤，对低通滤波之后的梯度场进行自适应滤波以恢复数据残差中的地层边界信息。

根据本发明的一个实施例，所述自适应滤波步骤包括：

针对梯度场中的每个目标点确定包括多个子区域的自适应滤波区域；

分别计算子区域的方差；

选择方差最小的子区域进行低通滤波以保留目标点处的地层边界信息。

根据本发明的一个实施例，所述自适应滤波步骤中自适应滤波方法为最小均方误差LMS滤波方法或者递归最小平方RLS滤波方法。

根据本发明的一个实施例，所述低通滤波步骤中采用高斯低通滤波器对梯度场进行低通滤波，其中，高斯低通滤波器的高斯算子表示为：

${\mathrm{Wg}}_{i,j}=\frac{1}{{2\mathrm{\π\σ}}^2}{e}^{-\frac{{(i-k-1)}^2+{(j-k-1)}^2}{{2\mathrm{\σ}}^2}},$

其中，Wg_i,j为离散高斯卷积核函数，σ²为方差，k为高斯卷积核矩阵的维数，i为低通滤波区域的纵向坐标，j为低通滤波区域的横向坐标。

根据本发明的另一方面，还提供一种基于梯度处理的全波形反演系统，包括：

残差计算模块，其由初始速度模型正演得到的计算数据和野外观测数据确定数据残差；

梯度计算模块，其根据数据残差在初始速度模型中的反传波场以及正演得到的正传波场求取梯度场；

扰动模型计算模块，其对梯度场进行滤波得到保存地层边界信息的速度扰动模型；

更新模块，其根据速度扰动模型更新初始速度模型。

根据本发明的一个实施例，进一步包括：

判断模块，其判断更新后的初始速度模型是否满足预设的精度要求，若不满足精度要求，则调用残差计算模块、梯度计算模块、扰动模型计算模块和更新模块；若满足精度要求，则将更新后的初始速度模型确定为反演模型。

根据本发明的一个实施例，所述扰动模型计算模块进一步包括：

低通滤波单元，其对梯度场进行低通滤波以消除数据残差中的随机噪音；

自适应滤波单元，其对低通滤波之后的梯度场进行自适应滤波以恢复数据残差中的地层边界信息。

根据本发明的一个实施例，所述自适应滤波单元进一步包括：

区域划分子单元，其针对梯度场中的每个目标点确定包括多个子区域的自适应滤波区域；

方差计算子单元，其分别计算子区域的方差；

低通滤波子单元，选择方差最小的子区域进行低通滤波以保留目标点处的地层边界信息。

本发明对全波形反演中的数据残差模型进行梯度滤波处理，利用低通滤波器抑制随机噪声，利用自适应滤波器恢复低通滤波处理过程中损失的地质边界信息。从而在全波形反演的模型更新过程中，保留了地震资料中原有的地质构造信息，使得反演更加稳定，提高反演精度。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例共同用于解释本发明，并不构成对本发明的限制。

图1是本发明实施例一的基于梯度处理的全波形反演方法的步骤流程图；

图2是本发明实施例一中自适应滤波区域示例图；

图3是原始单炮记录的一个示例；

图4是包含噪声的单炮记录的一个示例；

图5是真实地质模型的一个示例；

图6是初始速度模型的一个示例；

图7是现有技术中的全波形反演结果；

图8是根据本发明实施例一的方法进行全波形反演的结果；

图9是本发明实施例二的基于梯度处理的全波形反演系统的结构示意图；

图10是本发明实施例二的自适应滤波单元的结构示意图。

具体实施方式

现有技术中陆地地震资料的噪音干扰严重，使得在全波反演过程中确定的扰动模型并不准确。现有的噪声处理方法可以抑制随机噪声，但是数据运算量大，计算繁琐，并且会导致原有地震资料中包含的地质构造边界模糊。

本发明的实施例中对全波形反演中的梯度进行低通和自适应两次滤波，通过低通滤波手段压制随机噪声，并通过自适应滤波恢复在低通滤波过程中损失的边界信息，达到边界保真的效果。进行得到高精度的扰动模型，使得反演更加稳定，反演精度更准确。

为使本发明的目的、技术方案和优点更加清楚，以下结合附图对本发明作进一步地详细说明。

实施例一

图1是本实施例的全波形反演梯度处理方法的步骤流程图。

首先，在步骤S101中，由初始速度模型正演得到的计算数据和野外观测数据确定数据残差。该数据残差表示为最小二乘目标函数的形式：

$C\left(m\right)=\frac{1}{2}\left[L^2\left(\mathrm{\Δd}\right)\right]=\frac{1}{2}{\mathrm{\Δd}}^{+}\mathrm{\Δd}---\left(1\right)$

其中，△d为残差矢量，即基于初始速度模型m₀计算出的数据与地震数据之差。上标+表示伴随算子。

随后执行步骤S102，根据数据残差在初始速度模型中的反传波场以及正演得到的正传波场求取梯度场。

其中，反传波场为检波点处数据残差在初始速度模型中的波场数据，正传波场表示震源在初始速度模型中的波场数据。梯度场表示速度模型更新量的方向。

假设公式(1)表示的目标函数是局部二次型，使用局部最优化方法在初始模型m₀附近处求得极值，其第一步是线性化问题。在m₀处对目标函数进行二阶泰勒公式展开舍去高次项可得：

$\frac{\∂C\left(m\right)}{\∂m}\≈\frac{\∂C\left(m_0\right)}{\∂m}+\frac{{\∂}^{2}C\left(m_0\right)}{{\∂m}^2}\mathrm{\Δm}---\left(2\right)$

其中，m₀处的梯度定义为表示在m₀处速度更新的方向，目标函数的二阶偏导数为海信矩阵。

在本步骤中，m₀处的梯度表示速度的误差方向，沿梯度方向更新速度模型可更快地接近真实的速度分布情况。

随后，在步骤S103中对梯度场进行滤波得到保留地层边界信息的速度扰动模型。

由于在目标函数的极值处为零，因此根据公式(2)可得速度扰动模型矢量△m如下

$\mathrm{\Δm}=-{\left[\frac{{\∂}^{2}C\left(m\right)}{{\∂m}^2}{|}_{m=m_0}\right]}^{-2}\frac{\∂C\left(m\right)}{\∂m}{|}_{m=m_0}---\left(3\right)$

由于地震数据采集过程中存在各种噪声污染，例如测量环境周围的震动、检波仪器本身的噪声以及数据处理过程中产生的噪声，求取的梯度存在严重的干扰。本实施例中对梯度采取低通滤波以消除随机噪声干扰，继而通过自适应滤波对低通滤波损害的构造边界予以恢复。

具体的，首先对梯度场进行低通滤波以消除数据残差中的随机噪音。

采用低通滤波器W_g对梯度进行滤波以抑制随机噪声，得到：

$\mathrm{\Δ}{m}_1=-{\left[\frac{{\∂}^{2}C\left(m\right)}{{\∂m}^2}{|}_{m=m_0}\right]}^{-2}{W}_g\frac{\∂C\left(m\right)}{\∂m}{|}_{m=m_0}---\left(4\right)$

优选的，采用高斯低通滤波器，其高斯算子表示为：

${\mathrm{Wg}}_{(i,j)}=\frac{1}{{2\mathrm{\π\σ}}^2}{e}^{-\frac{{(i-k-1)}^2+{(j-k-1)}^2}{{2\mathrm{\σ}}^2}}---\left(5\right)$

其中，Wg_i,j为离散高斯卷积核函数，σ²为方差，用于调节滤波的平滑程度，k为高斯卷积核矩阵的维数，i为低通滤波区域的纵向坐标，j为低通滤波区域的横向坐标。

如此以来可以抑制梯度场中的随机噪声，但同时也会使地震剖面图中的边缘和细节处产生严重模糊，损失地质边界构造信息。这些地质边界构造信息可在下述步骤中得到恢复。

然后，对低通滤波之后的梯度进行自适应滤波，以恢复数据残差中的地层边界信息，这样可得到最终的速度扰动模型矢量：

$\mathrm{\Δ}{m}_2=-{\left[\frac{{\∂}^{2}C\left(m\right)}{{\∂m}^2}{|}_{m=m_0}\right]}^{-2}{W}_a{W}_g\frac{\∂C\left(m\right)}{\∂m}{|}_{m=m_0}---\left(6\right)$

其中，W_a为自适应滤波器。

优选的，本实施例中采用的自适应滤波方法为最小均方误差LMS滤波方法。不限于此，还可以采用其他保护边界的自适应滤波器方法，比如递归最小平方RLS方法等。

以下结合图2对本实施例的自适应滤波方法进行说明。

首先，针对梯度场的每个目标点确定包括多个子区域的自适应滤波区域。图2显示了自适应滤波区域的示例图。滤波区域200为窗口大小为k，包含子区域201、202、203和204的正方形区域，中心点(x,y)为目标点。

对于滤波区域200的中心点(x，y)，分别在子区域201、202、203和204中计算方差。四个子区域都包含中心点(x，y)，方差最小的子区域说明该区域内速度信息变化趋势较小，较为平缓，相对其他几个子区域而言存在地层边界信息的可能性更小。因此，选择方差最小的子区域，即没有地层信息的子区域进行平滑滤波，而方差较大的子区域，即可能存在地层边界的区域中的梯度场保持不变，不做任何处理。从而可以有效保留目标点的自适应滤波区域中的地层边界信息。

对梯度场中的每一个目标点都做同样处理，那么在去除噪音的同时地层边界得到了有效保护。这样，最终的速度扰动模型矢量△m₂的精度更高，包含的有效信息更多。

在步骤S104中，根据速度扰动模型矢量更新初始模型。这样得到的初始模型中，能够保留原地震资料中的地质边界信息，在此基础上得到反演模型精度更高。

优选的，还包括判断步骤S105，判断更新后的初始速度模型是否满足预设的精度要求，若不满足精度要求，则迭代执行残差计算步骤S101、梯度计算步骤S102、扰动模型计算步骤S103和更新步骤S104；若满足精度要求，则将更新后的初始速度模型确定为反演模型。

这样在每次迭代过程中均对梯度场做自适应滤波处理，能够抑制随机噪声并保留原地层构造信息。进而可增强反演稳定性，提高反演精度，得到地层边界清晰的反演模型。

以下提供采用本发明的梯度处理方法对地震资料进行全波形反演的结果的示例。

图3是原始单炮记录，图4是包含噪声的单炮记录。由于噪声干扰，图4中地质边界很模糊。图5是真实地质模型，图7是现有技术中的全波形反演结果。可以看到，图7中的反演结果与真实地质模型相差很远，并不能反应真实的地层构造。

图6是根据本实施例的方法得到的初始模型，图8是根据本实施例的方法得到的全波形反演结果。可以看到噪声信息得到有效压制，构造部分更加清晰，对比度明显提高。

实施例二

本实施例提供一种基于基于梯度处理的全波形反演系统。如图9所示，该系统包括残差计算模块91，梯度计算模块92，扰动模型计算模块93和更新模块94。

其中，残差计算模块91由初始速度模型正演得到的计算数据和野外观测数据确定数据残差。梯度计算模块92与残差计算模块91连接，用于根据数据残差在初始速度模型中的反传波场以及正演得到的正传波场求取梯度场。扰动模型计算模块93对梯度场进行滤波得到保存地层边界信息的速度扰动模型。更新模块94根据速度扰动模型更新初始速度模型。

优选的，该系统还包括判断模块95，其判断更新后的初始速度模型是否满足预设的精度要求，若不满足精度要求，则调用残差计算模块91、梯度计算模块92、扰动模型计算模块93和更新模块94，以再次更新初始速度模型；若满足精度要求，则将更新后的初始速度模型确定为反演模型。

需要说明的是，扰动模型计算模块93进一步包括低通滤波单元931和自适应滤波单元932。其中，低通滤波单元931对梯度场进行低通滤波以消除数据残差中的随机噪音；自适应滤波单元932对低通滤波之后的梯度场进行自适应滤波以恢复数据残差中的地层边界信息。例如，自适应滤波单元932采用最小均方误差LMS滤波方法或者递归最小平方RLS滤波方法进行滤波。

如图10所示，自适应滤波单元932包括区域划分子单元933、方差计算子单元934和低通滤波子单元935。

其中，区域划分子单元933针对梯度场中的每个目标点确定包括多个子区域的自适应滤波区域。方差计算子单元934分别计算子区域的方差，低通滤波子单元935用于选择方差最小的子区域进行低通滤波以保留目标点处的地层边界信息。

本实施例中对梯度采取低通滤波以消除随机噪声干扰，继而通过自适应滤波对低通滤波损害的构造边界予以恢复，进而可增强反演稳定性，提高反演精度，得到地层边界清晰的反演模型。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (10)

1.一种基于梯度处理的全波形反演方法，其特征在于，包括以下步骤：

残差计算步骤，由初始速度模型正演得到的计算数据和野外观测数据确定数据残差；

梯度计算步骤，根据数据残差在初始速度模型中的反传波场以及正演得到的正传波场求取梯度场；

扰动模型计算步骤，对梯度场进行滤波得到保留地层边界信息的速度扰动模型；

更新步骤，根据速度扰动模型更新初始速度模型。

2.根据权利要求1所述的方法，其特征在于，还包括：

判断步骤，判断更新后的初始速度模型是否满足预设的精度要求，若不满足精度要求，则迭代执行残差计算步骤、梯度计算步骤、扰动模型计算步骤和更新步骤；若满足精度要求，则将更新后的初始速度模型确定为反演模型。

3.根据权利要求1或2所述的方法，其特征在于，所述扰动模型计算步骤包括：

低通滤波步骤，对梯度场进行低通滤波以消除数据残差中的随机噪音；

自适应滤波步骤，对低通滤波之后的梯度场进行自适应滤波以恢复数据残差中的地层边界信息。

4.根据权利要求3所述的方法，其特征在于，所述自适应滤波步骤包括：

针对梯度场中的每个目标点确定包括多个子区域的自适应滤波区域；

分别计算子区域的方差；

选择方差最小的子区域进行低通滤波以保留目标点处的地层边界信息。

5.根据权利要求3所述的方法，其特征在于，所述自适应滤波步骤中自适应滤波方法为最小均方误差LMS滤波方法或者递归最小平方RLS滤波方法。

6.根据权利要求4所述的方法，其特征在于，所述低通滤波步骤中采用高斯低通滤波器对梯度场进行低通滤波，其中，高斯低通滤波器的高斯算子表示为：

${\mathrm{Wg}}_{i,j}=\frac{1}{{2\mathrm{\π\σ}}^2}{e}^{-\frac{{(i-k-1)}^2+{(j-k-1)}^2}{{2\mathrm{\σ}}^2}},$

其中，Wg_i,j为离散高斯卷积核函数，σ²为方差，k为高斯卷积核矩阵的维数，i为低通滤波区域的纵向坐标，j为低通滤波区域的横向坐标。

7.一种基于梯度处理的全波形反演系统，其特征在于，包括：

残差计算模块，其由初始速度模型正演得到的计算数据和野外观测数据确定数据残差；

梯度计算模块，其根据数据残差在初始速度模型中的反传波场以及正演得到的正传波场求取梯度场；

扰动模型计算模块，其对梯度场进行滤波得到保存地层边界信息的速度扰动模型；

更新模块，其根据速度扰动模型更新初始速度模型。

8.根据权利要求7所述的系统，其特征在于，进一步包括：

判断模块，其判断更新后的初始速度模型是否满足预设的精度要求，若不满足精度要求，则调用残差计算模块、梯度计算模块、扰动模型计算模块和更新模块；若满足精度要求，则将更新后的初始速度模型确定为反演模型。

9.根据权利要求8所述的系统，其特征在于，所述扰动模型计算模块进一步包括：

低通滤波单元，其对梯度场进行低通滤波以消除数据残差中的随机噪音；

自适应滤波单元，其对低通滤波之后的梯度场进行自适应滤波以恢复数据残差中的地层边界信息。

10.根据权利要求9所述的系统，其特征在于，所述自适应滤波单元进一步包括：

区域划分子单元，其针对梯度场中的每个目标点确定包括多个子区域的自适应滤波区域；

方差计算子单元，其分别计算子区域的方差；

低通滤波子单元，选择方差最小的子区域进行低通滤波以保留目标点处的地层边界信息。