CN107462924B - 一种不依赖于测井资料的绝对波阻抗反演方法 - Google Patents

Abstract

本发明公开了一种不依赖于测井资料的绝对波阻抗反演方法。首先，从叠前炮集地震数据出发，利用地震数据全部信息通过全波形反演的手段得到长波长P波速度模型；其次，基于长波长P波速度模型通过深时转换以及Gardner经验公式得到低频波阻抗模型；最后，使用低频波阻抗模型约束基于多组变异差分进化算法的波阻抗反演，得到绝对波阻抗模型。本发明直接从地震数据出发，通过利用全波形反演与全局优化算法，实现了在不依赖于测井资料的前提下，对横向强非均质介质绝对波阻抗参数的反演。因为使用了叠前地震数据全部信息，所以可以更好的构建横向连续性好的低频波阻抗模型；并且因为使用了全局优化算法求解波阻抗反演问题，能够反演得到全局最优波阻抗参数。

Description

一种不依赖于测井资料的绝对波阻抗反演方法

技术领域

本发明属于地球物理勘探领域，涉及一种地球模型物性参数反演方法，特别涉及一种不依赖于测井资料的绝对波阻抗反演方法。

背景技术

石油和天然气是事关国家经济发展及国家安全的战略资源，进一步推进我国油气资源的勘探与开发属于国家重大战略。岩性油气藏是我国今后相当长时间内最主要的勘探开发对象。地震勘探方法是油气勘探的主要方法，绝对波阻抗是油气检测的重要参数。我国岩性油气藏大多在地质上纵向呈薄互层结构、横向上具有强非均质性，如何对岩性油气藏的绝对波阻抗参数进行精细描述是我国油气勘探面临的关键问题。地震反演基于地震反射数据来定量的估计地下介质的岩石物理参数，进而进行储层预测和油藏特征描述。因此，地震反演是刻画岩性油气藏绝对波阻抗的有效方法。

在已有的波阻抗反演方法中，基于模型的反演是最常用的一类波阻抗反演方法。这一类方法的目标是在参数空间中寻找最佳的波阻抗参数使得其对应的正演地震数据与实际观测数据能够最佳的匹配。由于地震数据往往带宽有限并且含噪，这些因素导致波阻抗反演是不适定的，很难唯一恢复出地下介质的波阻抗参数。通常为了克服这些问题，基于模型的波阻抗反演方法依赖于引入适当的约束来使得反演更加稳定并且能够得到绝对波阻抗。最为常用的约束来自于测井资料。基于模型的波阻抗反演方法往往使用通过测井资料沿着层位进行插值得到波阻抗参数低频成分作为约束，进而最终反演得到绝对波阻抗。测井资料具有纵向分辨率高这一优势，但其横向的不连续性制约着构造低频波阻抗模型的准确性。这一方法在井网分布较密的工区或地下介质横向变化不剧烈的工区有着良好的效果，但是在少井甚至无井的工区，特别是地下介质横向非均质性强的工区则会失效。这些不利因素使得完全依赖测井资料来构造低频波阻抗模型在很多情况下是不可行的。

基于模型的波阻抗反演是一类典型的非线性反演问题，当使用线性化方法求解这类反问题时，算法容易陷入局部极值点从而无法得到最优的模型参数。为了解决这一问题，一些学者使用全局优化类方法(例如：模拟退火法、差分进化算法等)来求解这一反问题。这类方法可以有效的避免算法陷入局部极值点，但其固有的缺陷在于计算量巨大。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种不依赖于测井资料的绝对波阻抗反演方法，其直接从地震数据出发，通过利用全波形反演与全局优化算法，实现了在不依赖于测井资料的前提下，对横向强非均质介质绝对波阻抗参数的反演。

本发明的目的是通过以下技术方案来实现的：

这种不依赖于测井资料的绝对波阻抗反演方法为：首先，从叠前炮集地震数据出发，利用地震数据全部信息通过全波形反演的手段得到长波长P波速度模型；其次，基于长波长P波速度模型通过深时转换以及Gardner经验公式得到低频波阻抗模型；最后，使用低频波阻抗模型约束基于多组变异差分进化算法的波阻抗反演，得到绝对波阻抗模型。

进一步，以上方法具体包括以下步骤：

1)采集原始地震数据，然后对采集到的地震数据进行预处理，处理后得到叠前炮集地震数据，记为d_obs(t,x_r；x_s)，其中t表示时间变量；x_r,x_s分别为检波器和震源的位置；处理得到叠后反射地震数据，记为其中t表示时间变量；n为共深度点道集序号；

2)构建长方形网格地质模型，设定离散网格的空间离散间隔、网格规模以及正演模型时间采样间隔；

3)给定地质模型P波速度v(x)的初始模型，并指定待优化的目标函数J(v(x))；

4)通过最速下降法迭代得到目标函数的全局极小值，进而得到地质模型P波速度的长波长分量v_l(x)；

5)通过“深时转换”得到低频P波速度模型v_lf(t,n)，进而通过Gardner公式得到低频波阻抗模型I_lf(t,n)；

6)构建水平层状地质模型，给定地质模型的层数和层的厚度，每层的介质模型参数为波阻抗I；

7)通过使用步骤5)的结果确定参数I的搜索空间，并给定波阻抗反演待优化的目标函数；

8)使用全局优化算法对步骤7)的目标函数进行优化，求取其全局极小值点；

9)选出步骤8)迭代结束后全局优化算法群体中目标函数值最小的个体，将所述数值最小的个体解码后即得到最终搜索到的最优绝对波阻抗模型。

进一步，以上步骤2)中，离散网格的空间离散间隔、网格规模以及正演模型时间采样间隔是根据实际地震数据的炮检距、有效频带范围、采样时间决定的；选定网格参数的标准是使得基于该网格进行有限差分正演模拟时，不仅满足稳定性条件而且有效压制数值频散。

进一步，以上步骤3)中，地质模型P波速度v(x)的初始模型是对真实P波速度模型的估计，初始模型设定为横向均匀、纵向线性增加的速度模型；所述目标函数是刻画最优解的标准，以计算地震数据与观测地震数据之间的拟合程度或误差大小为标准。设定衡量观测炮集地震数据与基于模型的计算炮集地震数据之间的时移非线性数据残差的2范数为目标函数，即目标函数为：

式中，T为正演的时间长度；P⁽²⁾(.)为二阶时移非线性算子；S_s,r为一约束算子来保证只有震源和检波器分别位于x_s和x_r的波场数据被用来构建目标函数；u_cal(t,x；x_s)为震源位于x_s的计算波场；目标函数值J(v(x))越小说明速度模型v(x)越准确。

进一步，以上步骤4)中，P波速度的长波长分量是通过将一个各向同性高斯滤波器作用于反演得到的P波速度模型，滤掉其短波长分量得到。

进一步，以上步骤7)中，波阻抗参数的搜索空间是通过将步骤5)得到的低频波阻抗模型作为背景值，在背景值的基础上增加一设定值后作为搜索空间上界，在背景值的基础上减去一设定值后作为搜索空间下界；波阻抗反演的目标函数是以观测叠后地震数据与由模型合成的计算叠后地震数据之间误差绝对之和为目标函数，对于一道叠后地震数据，目标函数表示为：

式中，d_cal(t)为给定波阻抗参数时得到的正演单道叠后反射地震数据。

进一步，以上步骤8)中，使用的全局优化算法为多组变异差分进化算法，所述多组变异差分进化算法包括：初始化、变异、交叉和选择。所述多组变异差分进化算法基于一个规模为NP的群体来求解优化问题，群体中的每一个个体均为与模型参数维度相同的矢量；所述多组变异差分进化算法的流程如下：

8.1)在搜索空间范围内随机生成NP个初始个体构成群体；

8.2)对群体中的每一个个体进行变异操作；

8.3)对群体中的每一个个体进行交叉操作；

8.4)对群体中的每一个个体进行选择操作；

8.5)判断是否满足算法结束条件，若满足输出群体中的最优个体为优化问题的解；否则，跳回8.2)，进行下一次迭代。

与现有技术相比，本发明具有以下有益效果：

本发明不依赖于测井资料的绝对波阻抗反演方法与常用依赖测井资料约束的绝对波阻抗方法相比，更加适用于岩性油气藏波阻抗参数的刻画；当测井资料缺乏时，本发明的方法依旧能够直接从地震数据出发得到对绝对波阻抗参数的准确刻画。本发明具有上述优势的原因在于：首先，因为使用了叠前地震数据全部信息，所以该方法可以更好的构建横向连续性好的低频波阻抗模型；其次，因为使用了全局优化算法求解波阻抗反演问题，所以该方法能够反演得到全局最优波阻抗参数。

附图说明

图1是本发明流程示意图；

图2是目标区不同速度模型对比图：(a)真实模型；(b)初始模型；(c)反演速度模型的长波长分量；

图3是目标区不同速度模型对应的角道集对比图：(a)真实模型；(b)初始模型；(c)反演速度模型；

图4是目标区不同绝对波阻抗模型对比图：(a)真实波阻抗模型；(b)某商业软件反演波阻抗模型；(c)本发明所提出方法反演波阻抗模型。

具体实施方式

下面结合附图对本发明做进一步详细描述：

波阻抗参数是一种可以被应用于储层刻画以及烃类检测的与岩石性质相关的属性，波阻抗反演是一种基于地震反射数据来得到地下介质波阻抗参数的反演技术。本发明提出了一种不依赖于测井资料的绝对波阻抗反演方法。该方法首先，从叠前炮集地震数据出发，利用地震数据全部信息通过全波形反演的手段得到长波长P波速度模型；其次，基于长波长P波速度模型通过深时转换以及Gardner经验公式得到低频波阻抗模型；最后，使用低频波阻抗模型约束基于多组变异差分进化算法的波阻抗反演，得到绝对波阻抗模型。

本发明的物质基础是通过野外高分辨率地震采集设备采集到的大量地震数据。本发明所提出的不依赖于测井资料的绝对波阻抗反演方法框架如图1所示，具体步骤分别如下：

1)采集原始地震数据，然后对采集到的地震数据进行常规预处理，包括：静校正、动校正、噪声压制、去除面波等处理后得到叠前炮集地震数据(处理海上资料时还需做去除鬼波操作)，记为d_obs(t,x_r；x_s)，其中t表示时间变量；x_r,x_s分别为检波器和震源的位置；除此之外，经过偏移成像得到叠后反射地震数据，记为其中t表示时间变量；n为共深度点道集序号；

2)构建长方形网格地质模型，设定离散网格的空间离散间隔、网格规模以及正演模拟时间采样间隔；以二维反演为例，离散网格需要确定的参数有：横向网格间隔dx，纵向网格间隔dz，横向网格数Nx以及纵向网格数Nz，除此之外还需确定正演模拟时间采样间隔dt。实际中可以根据地震数据有效频带范围、地震记录时间长度、地震数据观测系统以及Courant-Friedrichs-Lewy(CFL)条件来确定相关参数；

3)地质模型P波速度v(x)的初始模型是对真实P波速度模型的一个粗略估计，一般可以通过走时反演得到，或者将初始模型设定为横向均匀、纵向线性增加的速度模型；所述目标函数是刻画最优解的标准，以计算地震数据与观测地震数据之间的拟合程度或误差大小为标准。本发明中设定衡量观测炮集地震数据与基于模型的计算炮集地震数据之间的时移非线性数据残差的2范数为目标函数，即目标函数为：

式中，T为正演的时间长度；S_s,r为一约束算子来保证只有震源和检波器分别位于x_s和x_r的波场数据被用来构建目标函数；u_cal(t,x；x_s)为震源位于x_s的计算波场；P⁽²⁾(.)为二阶时移非线性算子。上式中计算波场u_cal(t,x；x_s)满足如下常密度声波方程：

式中，v(x)为空间位置x处的P波速度值；Δ为拉普拉斯算符；f(t)δ(x-x_s)为位置x_s处的震源函数。式(5)中使用的二阶时移非线性算子P⁽²⁾(.)定义如下：

式中，u(t)为一时间函数；“*”代表时间域卷积运算；D(·)代表时移线性变换；d(t)为一线性函数，a_n,N_ξ,ξ_n为d(t)的三个参数。本发明采取二阶时移非线性算子的原因在于该非线性算子具有将地震数据向低频集中并保持地震数据旅行时信息不被破坏的特性，适用于反演速度模型的长波长分量。

4)通过最速下降法迭代得到如式(5)所示的目标函数的极小值。最速下降法进行迭代时遵循如下的模型更新方式：

mⁿ⁺¹(x)＝mⁿ(x)+α_ndⁿ(x) (8)

式中，mⁿ(x),mⁿ⁺¹(x)分别为第n步以及第n+1步迭代的模型参数(此处为速度模型v(x))；α_n为迭代步长；dⁿ(x)为迭代方向，在最速下降法中，迭代方向为负的梯度方向。步长α_n可以按照Pica等提出的方法进行求取。为了表述简洁，做如下标记约定，u_s(t；m(x))表示模型参数为m(x)，震源f_s(t)激发的波场。则步长α_n的求取公式如下：

其中，ε为一个较小的数，必须保证：

式(8)中的迭代方向dⁿ(x)是负的梯度方向，其中目标函数对模型参数的梯度可以表示为：

式中，q_s(t)被称为残差反传波场，其满足如下的波动方程：

其中u_s(t)＝u_cal(t,x；x_s)，d_s,r(t)＝d_obs(t,x_r；x_s)。

目标函数对模型参数的梯度可以采用误差反传方法进行高效的求解。

通过上述最速下降法迭代，可以得到最优的P波速度模型，将一个各向同性高斯滤波器作用于反演得到的P波速度模型，滤掉其短波长分量，从而得到长波长P波速度模型。

5)基于第4)步得到的长波长P波速度模型，首先使用“深时转换”将深度域的速度模型转化成时间域的低频速度模型，进而使用Gardner公式计算低频波阻抗模型。根据Gardner公式，地层波阻抗参数与速度参数的关系可以表示为：

I＝aV^m+1 (13)

式中，I为波阻抗参数；V为速度参数；a和m为两个参数，当密度单位取做g/cc并且速度单位取做m/s时，一般取a＝0.31,m＝0.25。

6)构建水平层状地质模型，给定地质模型的层数和层的厚度，每层的介质模型参数为波阻抗。其中层的厚度是时间的函数，一般与叠后地震数据的采样间隔有关；层数等于叠后地震记录长度除以层的厚度。

7)模型空间的搜索范围通过将第5)步得到的低频波阻抗模型作为背景值，在背景值的基础上增加某一设定值后作为搜索空间上界，在背景值的基础上减去某一设定值后作为搜索空间下界。波阻抗反演的目标函数是以观测叠后地震数据与由模型合成的计算叠后地震数据之间误差绝对之和为目标函数，对于某一道叠后地震数据，目标函数可以表示为：

式中，d_cal(t)为给定波阻抗参数时得到的正演单道叠后反射地震数据。

8)本发明中使用多组变异差分进化算法对式(14)定义的目标函数进行优化，求解最优波阻抗参数。多组变异差分进化算法主要包含：初始化、变异、交叉和选择操作。该算法基于一个规模为NP的群体来求解优化问题，群体中的每一个个体均为与模型参数维度相同的矢量。算法的基本流程简述如下：

第一步：在搜索空间范围内随机生成NP个初始个体构成群体；

第二步：对群体中的每一个个体进行变异操作；

第三步：对群体中的每一个个体进行交叉操作；

第四步：对群体中的每一个个体进行选择操作；

第五步：判断是否满足算法结束条件，若满足输出群体中的最优个体为优化问题的解；否则，跳回第二步，进行下一次迭代。

具体来说，如若搜索空间的第j维的上下界分别为和那么初始化群体中的第i个个体的第j维可以表示如下：

式中，r_i,j是一个在[0,1]之间均匀分布的随机数。在初始化后，多组变异差分进化算法会按照变异，交叉和选择三个基本操作来迭代的更新和进化群体中的所有个体，直到满足迭代停止条件为止。

变异操作：多组变异差分进化算法第G代群体中第i个个体的变异向量可以表示为：

式中，为“DE/rand/1”方法产生的临时变异向量；为协同变异差分进化算法的变异方法产生的临时变异向量；a∈(0,1)为常数。

交叉操作：交叉操作的目的是在当前个体和其对应的变异个体之间进行“基因”交换从而得到一个试探个体常用的交叉操作定义如下：

式中r_i,j∈U[0,1]，CR∈[0,1]为交叉率，其控制着群体的多样性并且保证算法不陷入局部极值，q∈[1,D]是一个随机选取的整数来保证即使CR＝0，试探个体中也有一维的信息来自于变异个体。

选择操作：选择操作根据当前个体和试探个体的目标函数值来更新当前个体，其具体形式如下所示：

9)选出第8)步迭代结束后全局优化算法群体中目标函数值最小的个体，将该个体解码后即得到最终搜索到的最优绝对波阻抗模型。

模型算例

本发明的具体实施过程被应用于某孔隙介质地质模型。该模型横向共有14500米，纵向共有3200米。从地表自上而下分别为水层、疏松砂岩、泥岩、致密砂岩、泥岩、致密砂岩和泥岩。在第三层致密砂岩中含有一个倒梯形目标层，内含三层，分别为：含气、含油以及含水的孔隙砂岩。目标区的速度模型如图2(a)所示。

给定如图2(b)所示的一维速度模型(从浅到深从1500m/s线性增加到3600m/s)作为初始速度模型，经过本发明所提出的方法的步骤1到步骤4进行反演得到的目标区反演速度模型的长波长分量如图2(c)所示。

图3所示的分别为真实速度模型、初始速度模型以及反演速度模型对应的角度域共成像点道集(角道集)。众所周知，角道集是否拉平是判断偏移速度长波长分量是否准确的准则。如图3所示，真实速度模型对应的角道集处于拉平状态；初始速度模型对应的角道集存在随着角度增大道集向浅层倾斜的现象，这说明初始速度模型是不准确的；反演速度模型对应的角道集与真实速度模型对应的角道集有较好的一致性，这说明反演速度模型的长波长分量是准确的。

图4所示的分别为真实绝对波阻抗模型、某商业软件反演绝对波阻抗模型以及本发明所提出的方法反演得到的绝对波阻抗模型。为了检验反演绝对波阻抗模型的正确性，图4(b)和(c)中投影了5口伪井处的真实绝对波阻抗曲线。如图4所示，商用软件的反演结果严重的依赖于测井资料构建的初始波阻抗模型，当反演中使用的测井数据不足时，商用软件构造的初始波阻抗模型不能准确的刻画出真实波阻抗模型的横向变化，进而得到错误的波阻抗反演结果。本发明所提出的绝对波阻抗反演方法在不依赖于测井数据的前提下，可以得到横向连续性好，与真实波阻抗模型一致性高的波阻抗反演结果。虽然部分区域的反演波阻抗参数值与真实波阻抗值有差异，但是本章所提出的方法能够很好的刻画不同层间的波阻抗差异。

Claims (8)

1.一种不依赖于测井资料的绝对波阻抗反演方法，其特征在于，首先，从叠前炮集地震数据出发，利用地震数据全部信息通过全波形反演的手段得到长波长P波速度模型；其次，基于长波长P波速度模型通过深时转换以及Gardner经验公式得到低频波阻抗模型；最后，使用低频波阻抗模型约束基于多组变异差分进化算法的波阻抗反演，得到绝对波阻抗模型，具体包括以下步骤：

1)采集原始地震数据，然后对采集到的地震数据进行预处理，处理后得到叠前炮集地震数据，记为d_obs(t,x_r；x_s)，其中t表示时间变量；x_r,x_s分别为检波器和震源的位置；处理得到叠后反射地震数据，记为，其中t表示时间变量；n为共深度点道集序号；

2)构建长方形网格地质模型，设定离散网格的空间离散间隔、网格规模以及正演模型时间采样间隔；

3)给定地质模型P波速度v(x)的初始模型，并指定待优化的目标函数J(v(x))；

4)通过最速下降法迭代得到目标函数的全局极小值，进而得到地质模型P波速度的长波长分量v_l(x)；

5)通过“深时转换”得到低频P波速度模型v_lf(t,n)，进而通过Gardner公式得到低频波阻抗模型I_lf(t,n)；

6)构建水平层状地质模型，给定地质模型的层数和层的厚度，每层的介质模型参数为波阻抗I；

7)通过使用步骤5)的结果确定参数I的搜索空间，并给定波阻抗反演待优化的目标函数F(I(t))；

8)使用全局优化算法对步骤7)的目标函数进行优化，求取其全局极小值点；

9)选出步骤8)迭代结束后全局优化算法群体中目标函数值最小的个体，将目标函数值最小的个体解码后即得到最终搜索到的最优绝对波阻抗模型。

2.根据权利要求1所述的不依赖于测井资料的绝对波阻抗反演方法，其特征在于，步骤2)中，离散网格的空间离散间隔、网格规模以及正演模型时间采样间隔是根据实际地震数据的炮检距、有效频带范围、采样时间决定的；选定网格参数的标准是使得基于该网格进行有限差分正演模拟时，不仅满足稳定性条件而且有效压制数值频散。

3.根据权利要求1所述的不依赖于测井资料的绝对波阻抗反演方法，其特征在于，步骤3)中，地质模型P波速度v(x)的初始模型是对真实P波速度模型的估计，初始模型设定为横向均匀、纵向线性增加的速度模型；所述目标函数是刻画最优解的标准，以计算地震数据与观测地震数据之间的拟合程度或误差大小为标准。

4.根据权利要求3所述的不依赖于测井资料的绝对波阻抗反演方法，其特征在于，设定衡量观测炮集地震数据与基于模型的计算炮集地震数据之间的时移非线性数据残差的2范数为目标函数，即目标函数为：

式中，T为正演的时间长度；P⁽²⁾(.)为二阶时移非线性算子；S_s,r为一约束算子来保证只有震源和检波器分别位于x_s和x_r的波场数据被用来构建目标函数；u_cal(t,x；x_s)为震源位于x_s的计算波场；目标函数值J(v(x))越小说明地质模型P波速度v(x)越准确。

5.根据权利要求1所述的不依赖于测井资料的绝对波阻抗反演方法，其特征在于，步骤4)中，P波速度的长波长分量是通过将一个各向同性高斯滤波器作用于反演得到的P波速度模型，滤掉其短波长分量得到。

6.根据权利要求1所述的不依赖于测井资料的绝对波阻抗反演方法，其特征在于，步骤7)中，波阻抗参数的搜索空间是通过将步骤5)得到的低频波阻抗模型作为背景值，在背景值的基础上增加一设定值后作为搜索空间上界，在背景值的基础上减去一设定值后作为搜索空间下界；波阻抗反演的目标函数是以观测叠后地震数据与由模型合成的计算叠后地震数据之间误差绝对之和为目标函数，对于一道叠后地震数据，目标函数表示为：

式中，d_cal(t)为给定波阻抗参数时得到的正演单道叠后反射地震数据。

7.根据权利要求1所述的不依赖于测井资料的绝对波阻抗反演方法，其特征在于，步骤8)中，使用的全局优化算法为多组变异差分进化算法，所述多组变异差分进化算法包括：初始化、变异、交叉和选择。

8.根据权利要求7所述的不依赖于测井资料的绝对波阻抗反演方法，其特征在于，所述多组变异差分进化算法基于一个规模为NP的群体来求解优化问题，群体中的每一个个体均为与模型参数维度相同的矢量；所述多组变异差分进化算法的流程如下：

8.1)在搜索空间范围内随机生成NP个初始个体构成群体；

8.2)对群体中的每一个个体进行变异操作；

8.3)对群体中的每一个个体进行交叉操作；

8.4)对群体中的每一个个体进行选择操作；

8.5)判断是否满足算法结束条件，若满足输出群体中的最优个体为优化问题的解；否则，跳回8.2)，进行下一次迭代。