CN107450102A - 基于分辨率可控包络生成算子的多尺度全波形反演方法 - Google Patents

Abstract

本发明公开了一种基于分辨率可控包络生成算子的多尺度全波形反演方法：首先，建立分辨率可控包络生成算子，通过调整参数控制反演使用地震数据的频带范围；其次，基于分辨率可控包络生成算子构建全波形反演方法，推导目标函数对模型参数的梯度以及迭代步长的计算公式；最后，使用最速下降法对模型参数进行迭代更新，直到满足算法结束条件为止。本发明具有更强的克服跳周问题的能力，更适用于岩性油气藏速度模型建模。特别是当地震数据中低频信息缺乏或主频过高时，优势更为明显。本发明无论地震数据中是否具有充足的低频信息，经过分辨率可控包络生成算子作用后的地震数据具有充足的低频信息；本发明的参数可调，可以实现多尺度反演。

Description

基于分辨率可控包络生成算子的多尺度全波形反演方法

技术领域

本发明属于地球物理勘探技术领域，涉及一种地球模型P波速度参数反演方法，特别涉及一种基于分辨率可控包络生成算子的多尺度全波形反演方法。

背景技术

石油和天然气是事关国家经济发展及国家安全的战略资源，进一步推进我国油气资源的勘探与开发属于国家重大战略。岩性油气藏是我国今后相当长时间内最主要的勘探开发对象。地震勘探方法是油气勘探的主要方法，速度参数是油气检测的重要参数。我国岩性油气藏大多在地质上纵向呈薄互层结构、横向上具有强非均质性，如何对岩性油气藏的速度参数进行精细描述是我国油气勘探面临的关键问题。全波形反演是一类直接从叠前地震数据出发，利用叠前炮集地震数据全部信息来得到高精度地下介质速度模型的方法。因此，全波形反演是刻画岩性油气藏速度参数的有效方法。

传统的全波形反演是一种以得到地下介质高分辨率参数模型为目标的数据匹配方法。全波形反演通过最小化一个衡量观测数据与计算数据匹配程度的目标函数来得到地下介质参数。最早将全波形反演方法推向实用的是Lailly和Tarantola的工作。他们提出可以使用由震源激发的正传波场以及由数据残差激发的反传波场进行互相关来构建目标函数对模型参数的梯度。这一误差反传技术克服了显式求取Frechet导数需要大量波场正演这一缺陷。传统全波形反演使用L²范数来衡量观测数据与计算数据的匹配程度，其应用于反射地震资料构建高分辨率速度模型时面临很多问题，例如：全波形反演要求观测数据中包含充足的低频信息或者初始模型足够准确来保证算法可以克服周波跳跃问题进而避免陷入局部极值。因此，降低全波形反演对初始模型的依赖程度就成为全波形反演领域的一个研究热点。

Bunks在时间域提出了多尺度全波形反演方法，之后Pratt提出了频率域全波形反演方法。这两种方法在反演过程中首先使用低频数据进行反演。由于低频数据构成的目标函数相比于高频数据更加线性，全局极小值的收敛域更大，所以此时算法对初始模型的要求就会更低。基于低频数据的反演结果可以作为基于高频数据的全波形反演的初始模型。经过在反演中逐渐加入高频数据，最终可以得到高精度的速度模型。虽然多尺度方法有一些成功的应用案例，但是这类方法并不能从根本上克服跳周问题。当这类方法应用于实际资料时会遇到如下问题：通常建立一个足够精确的初始模型是很困难的；实际资料的低频信息往往被噪音影响；实际资料的偏移距往往是有限的。

很多学者研究了为全波形反演建立比L²范数更优的度量准则这一问题。他们的目标是建立凸性更好，鲁棒性更强并且能够很好克服跳周问题的目标函数。与传统方法进行波形匹配不同，这一方法强调对波场旅行时信息的匹配。这一思路具有代表性的工作就是Luo和Schuster提出的波动方程旅行时反演。与传统旅行时反演方法依赖于人工拾取波场旅行时信息不同，该方法使用波动方程来得到正演数据，可以被视为一种使用互相关目标函数的全波形反演方法。然而，由于波动方程旅行时反演方法产生的速度模型过于光滑，其不能保证可以为全波形反演建立足够精确的初始模型。Chi等以及Luo等提出了在全波形反演中使用波场的瞬时振幅(包络)来构建反演的目标函数。Wu等证明了地震数据包络的振荡和衰减携带了超低频的信息，这一信息可以被用来估计长波长的速度模型。但是这一目标函数仍然会受到周波跳跃问题的困扰。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种基于分辨率可控包络生成算子的多尺度全波形反演方法。

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

这种基于分辨率可控包络生成算子的多尺度全波形反演方法，包括以下步骤：

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

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

3)分析地震数据，给出反演使用的地震数据频带范围；

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

5)求取目标函数J(v(x))对模型参数v(x)的梯度，使用误差反传技术进行高效求解；

6)求取最速下降法迭代时的迭代步长，在频率域使用线性搜索方法确定；

7)通过最速下降法迭代得到目标函数的全局极小值，进而得到在第 3)步给定的频带范围内，P波速度的反演结果；

8)调整反演时使用的地震数据的频带范围，并使用步骤7)的结果作为初始模型，重复步骤3)到7)，循环完所有的有效频带范围，输出最优P波速度模型。

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

进一步，上述步骤3)中，首先分析接收地震数据的有效频带范围，在反演初期先使用相对低频的地震数据，随着反演的进行，逐渐扩大反演使用的频带范围。

进一步，上述步骤4)中，地质模型P波速度v(x)的初始模型是对真实P波速度模型的一个粗略估计，初始模型设定为横向均匀、纵向线性增加的速度模型；所述目标函数是刻画最优解的标准，以计算地震数据与观测地震数据之间的拟合程度或误差大小为标准。

进一步的，设定衡量观测炮集地震数据与基于模型的计算炮集地震数据之间的分辨率可控包络数据残差的2范数为目标函数，即目标函数为：

式中，T为正演的时间长度；G(.)为分辨率可控包络生成算子；S_s,r为一约束算子来保证只有震源和检波器分别位于x_s和x_r的波场数据被用来构建目标函数；u_cal(t，x；x_s)为震源位于x_s的计算波场；目标函数值J(v(x)) 越小说明速度模型v(x)越准确。

进一步，上述步骤5)中，目标函数J(v(x))对模型参数v(x)的梯度是使用误差反传技术进行高效求解。

进一步，上述步骤6)中，迭代步长是在频率域使用线性搜索方法进行确定。

进一步，上述步骤8)中，当使用某一频带范围内的地震数据进行迭代反演后，下一轮反演使用的地震数据频带范围将比上一轮迭代时宽，使更多高频的信息被用于进行反演。

本发明具有以下有益效果：

本发明所提出方法的步骤可以概括为以下三步：首先，本发明建立了一种分辨率可控包络生成算子，该算子可以通过调整参数控制反演使用地震数据的频带范围；其次，本发明基于分辨率可控包络生成算子，构建了全波形反演方法，推导了其目标函数对模型参数的梯度以及迭代步长的计算公式；最后，使用最速下降法对模型参数进行迭代更新，直到满足算法结束条件为止。

与常用全波形反演方法相比，该方法具有更强的克服跳周问题的能力，更适用于岩性油气藏速度模型建模。特别是当地震数据中低频信息缺乏或主频过高时，该方法的优势更为明显。本发明所提出方法具有上述优势的原因在于：首先，无论地震数据中是否具有充足的低频信息，经过分辨率可控包络生成算子作用后的地震数据具有充足的低频信息；其次，分辨率可控包络生成算子的参数可调，可以实现多尺度反演。

附图说明

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

图2给定三组不同参数时移线性算子的振幅谱；

图3是Marmousi模型算例的结果：(a)真实速度模型；(b)初始速度模型；(c)0-8Hz数据反演速度模型；(d)0-12Hz数据反演速度模型； (e)0-16Hz数据反演速度模型；(f)0-20Hz数据繁衍速度模型。

具体实施方式

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

高精度的速度模型是一种可以被应用于储层刻画以及烃类检测的与岩石性质相关的属性，全波形反演是一种基于叠前炮集地震数据来得到地下介质高精度速度模型的反演技术。本发明提出了一种基于分辨率可控包络生成算子的多尺度全波形反演方法。该方法直接从地震数据出发，通过利用地震数据不同频带的信息以及非线性算子构建了非线性多尺度全波形反演方法，大大提高了克服周波跳跃问题并精细速度模型的能力。

本发明的物质基础是通过野外高分辨率地震采集设备采集到的大量地震数据。本发明所提出的基于分辨率可控包络生成算子的多尺度全波形反演方法框架如图1所示，具体步骤分别如下：

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

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

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

3)首先分析接收地震数据的有效频带范围，在反演初期先使用相对低频的地震数据，随着反演的进行，逐渐扩大反演使用的频带范围。

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

式中，T为正演的时间长度；S_s,r为一约束算子来保证只有震源和检波器分别位于x_s和x_r的波场数据被用来构建目标函数；u_cal(t,x；x_s)为震源位于x_s的计算波场；G(.)为分辨率可控包络生成算子。目标函数值 J(v(x))越小说明速度模型v(x)越准确。上式中计算波场u_cal(t，x；x_s)满足如下常密度声波方程：

式中，v(x)为空间位置x处的P波速度值；Δ为拉普拉斯算符；f(t)为震源函数。式(3)中使用的分辨率可控包络生成算子G(.)定义如下：

式中，“*”代表时间域卷积运算；H(.)代表希尔伯特变换；D(·)代表时移线性变换；d(t)为某一线性函数，a_n,N_ξ,ξ_n为d(t)的三个参数。式(5) 中d(t)的振幅谱与其参数的关系如下：

式中，ω_s为采样频率。如图2所示为给定三组不同的参数时，d(t)振幅谱的形态。图2验证了d(t)的振幅谱是可以通过调整其参数进行调整的。

5)求取目标函数对模型参数的梯度，可以使用误差反传技术进行高效求解。目标函数J(v(x))对模型参数v(x)的梯度是使用误差反传技术进行高效求解；该求解步骤的关键在于确定误差反传波场的震源项。

使用u_s(t)表示u_cal(t,x；x_s)，d_s,r(t)表示d_obs(t,x_r；x_s)以及f_s(t)表示 f(t)δ(x-x_s)。考虑到初始条件，由于震源f_s(t)激发得到的压力场u_s(t)可以表示为：

式中，m为模型参数，即v(x)，使用伴随状态法，可以得到如下的伴随状态方程：

式中，y(t)即为误差反传的震源项。最终梯度的表达式为：

式中，m为模型参数；q_s满足如式(8)所定义的波动方程。

6)求取最速下降法迭代时的步长，在频率域使用线性搜索方法确定。迭代步长是在频率域使用线性搜索方法进行确定；在频率域求取迭代步长的原因在于反演中只使用了地震数据部分频带的信息，只有用这部分信息来确定迭代时使用的步长才是合理的。为了表述简洁，做如下标记约定，表示模型参数为m(x)，震源f_s(t)激发的波场的频率域表达式；为d_obs(t,x_r；x_s)的频率域表达式。那么第n步迭代的步长α_n可以按照如下公式进行求取：

式中，mⁿ(x)为第n步迭代后的模型参数；dⁿ(x)为第n+1步的负梯度方向；ε为一个较小的数，必须保证：

式中，max(·)为取最大值函数。

7)通过最速下降法迭代得到目标函数的全局极小值，进而得到在步骤3)步给定的频带范围内，P波速度的反演结果。其中最速下降法的迭代公式如下：

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

式中，mⁿ(x),mⁿ⁺¹(x)分别为第n步以及第n+1步迭代的模型参数(此处为速度模型v(x))；α_n为迭代步长；dⁿ(x)为迭代方向，在最速下降法中，迭代方向为负的梯度方向。

通过上述最速下降法迭代，可以得到在当前数据频带范围内最优的P 波速度模型。

8)调整反演时使用的地震数据的频带范围，并使用步骤7)的结果作为初始模型，重复步骤3)到7)，循环完所有的有效频带范围，输出最优P波速度模型。当使用某一频带范围内的地震数据进行迭代反演后，下一轮反演使用的地震数据频带范围将比上一轮迭代时宽，使更多高频的信息被用于进行反演。

一般调整反演使用的地震数据频带范围的方式是随着反演的进行逐渐扩大反演使用的数据频带，逐渐将高频的地震数据加入到反演中来。

模型算例

本发明的具体实施过程被应用于Marmousi速度模型。该模型横向共有5090米，纵向共有1330米，该模型的真实速度模型如图3(a)所示。给定如图3(b)所示的一维速度模型(由浅到深从1500m/s线性增加到 3500m/s)作为初始速度模型，该模型是对地下介质真实速度模型的粗略估计。

试验中共在地表均匀设置48炮，炮间距为100米，每一炮均由510 个均匀布设在地表的检波器进行接收，检波器间距为10米。实验中使用主频为20Hz的Ricker子波作为震源，实验中设定横向与纵向的采样间隔为10米，时间采样间隔为0.001s，正演共进行3001个时刻点，炮集记录时间为3s。实验中，对比了本发明提出的方法与一种公认的具有较强克服周波跳跃能力的全波形反演方法：包络反演。

首先，从图3(b)所示的初始模型出发，实验中使用本发明所提出的非线性多尺度全波形反演方法来对地下介质速度模型进行估计。实验中共使用了地震数据4个不同的频带范围：0-8Hz、0-12Hz、0-16Hz以及 0-20Hz，进行反演，得到的反演结果分别如图3(c)-(f)所示。其次，从图 3(b)所示的初始模型出发，实验中使用包络反演方法来对地下介质速度模型进行估计，图3(g)所示为使用包络反演方法反演得到的速度模型。

从图3中所示的实验结果可以看出，包络反演方法无法有效的解决这一反演问题，其得到的反演速度模型与真实速度模型无论在结构以及数值上有明显的差异。与包络反演方法形成鲜明的对比，本发明所提出的非线性多尺度能够得到对速度模型的良好刻画，并且其刻画精度随着使用地震数据频带的扩展而逐步增高。这一实验结果证实了本发明所提出方法具有比包络反演方法更强的克服周波跳跃的能力，进一步证实了本发明所提出的方法的优越性。

Claims (8)

1.一种基于分辨率可控包络生成算子的多尺度全波形反演方法，其特征在于，包括以下步骤：

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

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

3)分析地震数据，给出反演使用的地震数据频带范围；

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

5)求取目标函数J(v(x))对模型参数v(x)的梯度，使用误差反传技术进行高效求解；

6)求取最速下降法迭代时的迭代步长，在频率域使用线性搜索方法确定；

7)通过最速下降法迭代得到目标函数的全局极小值，进而得到在第3)步给定的频带范围内，P波速度的反演结果；

8)调整反演时使用的地震数据的频带范围，并使用步骤7)的结果作为初始模型，重复步骤3)到7)，循环完所有的有效频带范围，输出最优P波速度模型。

2.根据权利要求2所述的基于分辨率可控包络生成算子的多尺度全波形反演方法，其特征在于，所述步骤2)中，离散网格的空间离散间隔、网格规模以及正演模型时间采样间隔是根据实际地震数据的炮检距、有效频带范围和采样时间决定的；选定网格参数的标准是使得基于该网格进行有限差分正演模拟时，不仅满足稳定性条件而且有效压制数值频散。

3.根据权利要求2所述的基于分辨率可控包络生成算子的多尺度全波形反演方法，其特征在于，所述步骤3)中，首先分析接收地震数据的有效频带范围，在反演初期先使用相对低频的地震数据，随着反演的进行，逐渐扩大反演使用的频带范围。

4.根据权利要求2所述的基于分辨率可控包络生成算子的多尺度全波形反演方法，其特征在于，所述步骤4)中，地质模型P波速度v(x)的初始模型是对真实P波速度模型的一个粗略估计，初始模型设定为横向均匀、纵向线性增加的速度模型；所述目标函数是刻画最优解的标准，以计算地震数据与观测地震数据之间的拟合程度或误差大小为标准。

5.根据权利要求5所述的基于分辨率可控包络生成算子的多尺度全波形反演方法，其特征在于，设定衡量观测炮集地震数据与基于模型的计算炮集地震数据之间的分辨率可控包络数据残差的2范数为目标函数，即目标函数为：

式中，T为正演的时间长度；G(.)为分辨率可控包络生成算子；S_s,r为一约束算子来保证只有震源和检波器分别位于x_s和x_r的波场数据被用来构建目标函数；u_cal(t,x；x_s)为震源位于x_s的计算波场；目标函数值J(v(x))越小说明速度模型v(x)越准确。

6.根据权利要求2所述的基于分辨率可控包络生成算子的多尺度全波形反演方法，其特征在于，所述步骤5)中，目标函数J(v(x))对模型参数v(x)的梯度是使用误差反传技术进行高效求解。

7.根据权利要求2所述的基于分辨率可控包络生成算子的多尺度全波形反演方法，其特征在于，所述步骤6)中，迭代步长是在频率域使用线性搜索方法进行确定。

8.根据权利要求2所述的基于分辨率可控包络生成算子的多尺度全波形反演方法，其特征在于，所述步骤8)中，当使用某一频带范围内的地震数据进行迭代反演后，下一轮反演使用的地震数据频带范围将比上一轮迭代时宽，使更多高频的信息被用于进行反演。