CN110174698B

CN110174698B - 一种基于方位傅里叶系数的弹性阻抗反演方法及系统

Info

Publication number: CN110174698B
Application number: CN201910569030.7A
Authority: CN
Inventors: 潘新朋; 张广智; 李林; 王东琰; 周顺心
Original assignee: China University of Petroleum East China; Central South University
Current assignee: China University of Petroleum East China; Central South University
Priority date: 2019-06-27
Filing date: 2019-06-27
Publication date: 2020-06-30
Anticipated expiration: 2039-06-27
Also published as: CN110174698A

Abstract

本公开提出了一种基于方位傅里叶系数的弹性阻抗反演方法及系统，包括：利用方位部分叠加地震数据，通过约束稀疏脉冲反演得到方位弹性阻抗；利用反演得到的所有方位弹性阻抗提取不同入射角下的傅里叶系数，利用提取的二阶项及四阶项傅里叶系数进一步反演得到裂缝法向弱度和裂缝切向弱度，用于表征裂缝发育强度。本公开所提出的方法具有一定的抗噪性。模型测试和实际应用验证了所提出的方法的有效性。

Description

一种基于方位傅里叶系数的弹性阻抗反演方法及系统

技术领域

本公开涉及地震监测技术领域，特别是涉及一种基于方位傅里叶系数的弹性阻抗反演方法及系统。

背景技术

利用地震反射振幅随偏移距和方位角变化(AVOAz，Amplitude variation withoffset and azimuth)开展裂缝弱度参数反演对于裂缝储层表征具有重要意义。然而当地震数据具有较低信噪比时，利用传统的AVOAz反演方法提取裂缝弱度参数具有巨大的挑战性。

多参数AVOAZ反演问题是病态的，传统的AVOAZ反演是直接利用方位地震振幅进行各向同性参数和裂缝弱度参数的同步反演，由于在反射系数表达式中，裂缝弱度参数项的系数相比于各向同性参数项的系数较小，裂缝弱度对反射系数的贡献远小于各向同性部分对反射系数的贡献，且各向同性参数与裂缝弱度参数之间存在耦合作用，使得各向同性参数和裂缝弱度参数的同步反演变得不稳定。

发明内容

本说明书实施方式的目的是提供一种基于方位傅里叶系数的弹性阻抗反演方法，该方法用于实现裂缝弱度参数预测，能够有效的表征裂缝发育强度。

本说明书实施方式提供一种基于方位傅里叶系数的弹性阻抗反演方法，通过以下技术方案实现：

包括：

利用方位部分叠加地震数据，通过约束稀疏脉冲反演得到方位弹性阻抗；

利用反演得到的所有方位弹性阻抗提取不同入射角下的傅里叶系数，利用提取的二阶项及四阶项傅里叶系数进一步反演得到裂缝法向弱度和裂缝切向弱度，用于表征裂缝发育强度。

进一步的技术方案，对于方位部分叠加地震数据具体为：每个入射角对应的方位地震子波以及通过裂缝岩石物理建模估测的法向裂缝弱度和切向裂缝弱度。

进一步的技术方案，反演得到裂缝法向弱度和裂缝切向弱度时，在贝叶斯反演框架下，利用提取的傅里叶系数，通过柯西约束正则化和低频模型约束正则化的迭代重加权最小二乘法反演得到裂缝弱度参数。

进一步的技术方案，将方位弹性阻抗方程两边同时取对数，进行傅里叶级数展开，对于X个规则采样的方位地震数据，傅里叶系数可根据离散傅里叶变换计算得到。

进一步的技术方案，基于方位弹性阻抗方程获得各阶傅里叶系数与弹性参数及裂缝弱度参数之间的关系。

进一步的技术方案，考虑时间采样点为M个，入射角为N个，基于各阶傅里叶系数与弹性参数及裂缝弱度参数之间的关系获得矩阵表达式，基于该矩阵表达式获得待反演模型参数的后验概率密度函数；

假设似然函数p(d|m)服从于高斯分布；

假设模型参数先验概率分布服从柯西分布；

求解得到最大后验概率初始目标函数；

结合从裂缝岩石物理建模估测得到裂缝弱度初始低频信息，得到最终的反演目标函数，基于目标函数反演得到裂缝弱度参数。

本说明书实施方式提供一种基于方位傅里叶系数的弹性阻抗反演系统，通过以下技术方案实现：

包括：

方位弹性阻抗获得模块，被配置为：利用方位部分叠加地震数据，通过约束稀疏脉冲反演得到方位弹性阻抗；

裂缝发育强度表征模块，被配置为：利用反演得到的所有方位弹性阻抗提取不同入射角下的傅里叶系数，利用提取的二阶项及四阶项傅里叶系数进一步反演得到裂缝法向弱度和裂缝切向弱度，用于表征裂缝发育强度。

与现有技术相比，本公开的有益效果是：

本公开首先基于傅里叶级数展开方法，将HTI介质的方位弹性阻抗表征为傅里叶系数表达式，然后利用方位地震数据反演得到所有的方位弹性阻抗，进一步提取傅里叶系数，最后在贝叶斯反演框架下，采用柯西约束正则化和低频模型约束正则化的迭代重加权最小二乘法反演得到裂缝弱度参数。所提出的方法具有一定的抗噪性。模型测试和实际应用验证了所提出的方法的有效性。

本公开的输入为利用方位地震数据反演得到的稳定可靠的方位弹性阻抗数据，具有较高信噪比的优势；其次我们利用傅里叶级数展开方法将方位弹性阻抗重写为傅里叶系数的表达式，表达式中的直流分量包含了各向同性参数和裂缝弱度参数的共同影响，而二阶项和四阶项傅里叶系数仅包含了裂缝弱度参数的影响，该方法消除了各向同性参数与裂缝弱度参数之间的耦合作用，有效地减少了待反演参数个数。此外，相比于反射系数，二阶项和四阶项傅里叶系数对裂缝弱度的变化更加敏感，利用二阶项和四阶项傅里叶系数，在贝叶斯反演框架下，采用柯西约束正则化和低频模型约束正则化的迭代重加权最小二乘法反演得到稳定可靠的裂缝弱度参数，提高了反演的稳定性和准确性。最后，该方法将方位弹性阻抗通过傅里叶变换转化为傅里叶系数作为后续反演的输入，可以有效地节约存储空间，提高反演运算效率。

附图说明

构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。

图1为本公开实施例子的反演流程示意图；

图2为本公开实施例子的测井估算的纵、横波阻抗及密度曲线图；

图3为本公开实施例子的裂缝岩石物理建模估算的法向裂缝弱度和切向裂缝弱度图；

图4(a)-图4(c)为本公开实施例子的(入射角为10°、入射角为20°、入射角为30°)合成的无噪音的对数域归一化方位弹性阻抗图；

图5(a)-图5(c)为本公开实施例子的(入射角为10°、入射角为20°、入射角为30°)合成的含有30％高斯噪音的对数域归一化方位弹性阻抗图；

图6为本公开实施例子的无噪音情况下反演的法向和切向裂缝弱度图；

图7为本公开实施例子的含有30％高斯噪音情况下反演的法向和切向裂缝弱度图；

图8(a)-图8(c)为本公开实施例子的不同方位部分角度(平均入射角为10°、平均入射角为20°及平均入射角为30°)叠加地震剖面图；

图9(a)-图9(c)为本公开实施例子的不同平均入射角(平均入射角为10°、平均入射角为20°、平均入射角为30°)反演的方位弹性阻抗剖面图；

图10为本公开实施例子的反演的法向裂缝弱度和切向裂缝弱度。

具体实施方式

应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

实施例子一

该实施例公开了一种基于方位傅里叶系数的弹性阻抗反演方法，首先基于傅里叶级数展开方法，将HTI介质的方位弹性阻抗表征为傅里叶系数表达式，然后利用方位地震数据反演得到所有的方位弹性阻抗，进一步提取傅里叶系数，最后在贝叶斯反演框架下，采用柯西约束正则化和低频模型约束正则化的迭代重加权最小二乘法反演得到裂缝弱度参数。

首先关于HTI介质中的方位弹性阻抗方程如下：

其中，

a(θ)＝sec²θ，b(θ)＝-8g sin²θ，c(θ)＝4gsin²θ-tan²θ，

e(θ，φ)＝2g cos²(φ-φ_sym)sin²θ[1-sin²(φ-φs_ym)tan²θ]。

方程(1)中，g＝V_S ²/V_P ²为上下地层横波速度均值平方与纵波速度均值平方之比，I_P，I_S和r为地层的纵、横波阻抗及密度。下标0表示弹性参数的均值。D_N和D_T表示法向裂缝弱度和切向裂缝弱度。q为入射角，f为地震测线方位角，f_sym为裂缝对称轴方位角。

EI表示方位弹性阻抗，与各向同性参数和裂缝弱度有关，其表达式如下：

井中的EI曲线可以通过井中的各向同性参数和由裂缝各向异性岩石物理建模估测的裂缝弱度参数计算得到。

在利用方位地震数据反演方位弹性阻抗时，可以利用井中计算的方位弹性阻抗曲线构建待反演工区的方位弹性阻抗低频模型。

将方程(1)两边同时取对数，可得到如下方程：

L_EI(φ，θ)＝a(θ)L_P+b(θ)L_S+c(θ)L_ρ+d(φ，θ)Δ_N+e(φ，θ)Δ_T (2)

其中，

对方程(2)进行傅里叶级数展开，可得到如下表达式：

L_EI(φ，θ)＝A₀(θ)+A₂(θ)cos(2φ)+B₂(θ)sin(2φ)

+A₄(θ)cos(4φ)+B₄(θ)sin(4φ) (3)

其中，A_n(θ)和B_n(θ)(n＝0，2，4)称为傅里叶系数。

对于X个规则采样的方位地震数据，傅里叶系数可根据离散傅里叶变换计算得到：

各阶傅里叶系数与弹性参数及裂缝弱度参数之间的关系可通过推导方程(1)得到：

A₀(θ)＝a(θ)L_P+b(θ)L_S+c(θ)L_ρ+f(θ)Δ_N+g(θ)Δ_T (6)

A₂(θ)＝h(θ)cos(2φ_sym)Δ_N+i(θ)cos(2φ_sym)Δ_T (7)

B₂(θ)＝h(θ)sin(2φ_sym)Δ_N+i(θ)sin(2φ_sym)Δ_T (8)

A₄(θ)＝j(θ)cos(4φ_sym)Δ_N+k(θ)cos(4φ_sym)Δ_T (9)

B₄(θ)＝j(θ)sin(4φ_sym)Δ_N+k(θ)sin(4φ_sym)Δ_T (10)

其中，

h(θ)＝g(2g-1)sin²θ+g(g1)sin²θ tan²θ，

i(θ)＝g sin²θ，

考虑时间采样点为M个，入射角为N个，则联立方程(7)，(8)，(9)和(10)可得如下矩阵表达式：

其中，

i(θ_m)＝diag[i¹(θ_m)，...，i^N(θ_m)]^T，

k(θ_m)＝diag[k¹(θ_m)，...，k^N(θ_m)]^T，

上标T表示矩阵的转置，符号diag表示对角阵，下标m和n代表第m个入射角和n阶傅里叶系数。

方程(11)可简化为：

d_NM×1＝G_NM×4N·m_4N×1

(12)

其中，

基于贝叶斯理论，待反演模型参数的后验概率密度函数可表示为：

假设似然函数p(d|m)服从于高斯分布：

柯西分布可以最大限度地提高垂直分辨率和保留弱小反射系数，假设模型参数先验概率分布服从柯西分布：

其中，

和

分别表示噪音方差和模型参数方差。

结合方程(14)和方程(15)，可求解得到最大后验概率初始目标函数如下：

结合从裂缝岩石物理建模估测得到裂缝弱度初始低频信息，得到最终的反演目标函数：

其中，

和

分别为法向裂缝弱度和切向裂缝弱度的约束系数，Δ_N0,Δ_T0分别为法向裂缝弱度和切向裂缝弱度的初始模型，可以从裂缝岩石物理建模估测得到。

具体实施时反演流程：

(1)数据准备：包括方位部分角度叠加地震数据，每个入射角对应的方位地震子波以及通过裂缝各向异性岩石物理建模估测的法向裂缝弱度和切向裂缝弱度。

(2)反演方位弹性阻抗：利用方位部分角度叠加地震数据，通过约束稀疏脉冲反演得到方位弹性阻抗。

(3)估测傅里叶系数：利用反演得到的所有方位弹性阻抗提取不同入射角下的傅里叶系数。

(4)反演裂缝弱度参数：在贝叶斯反演框架下，利用提取的傅里叶系数，通过柯西约束正则化和低频模型约束正则化的迭代重加权最小二乘法反演得到裂缝弱度参数。

本公开的另一实施例子公开了模型测试：使用单口井数据进行合成实验来验证我们提出方法的合理性。如图2和图3所示，测试井中包含有测井估算的纵、横波阻抗，密度以及经过裂缝岩石物理建模估算得到的法向裂缝弱度和切向裂缝弱度。选用35Hz的雷克子波，六个方位(分别为0°，30°，60°，90°，120°及150°)和三个入射角(分别为10°，20°及30°)，利用方程(2)合成方位弹性阻抗，并添加30％的高斯噪音，分别进行反演测试。图4(a)-图4(c)和图5(a)-图5(c)分别为合成的无噪音情况下和含有30％高斯噪音情况下的对数域归一化方位弹性阻抗。图6和图7分别为对数域归一化方位弹性阻抗无噪音情况下和含有30％高斯噪音情况下的裂缝弱度参数反演结果，其中，黑色实线表示真实值，黑色虚线表示初始模型，黑色“+”线表示反演结果。从图6和图7可以看出，在无噪音情况下，裂缝弱度参数反演结果与真实值完全一致；而在含有30％高斯噪音情况下，裂缝弱度参数反演结果与真实值变化趋势基本一致，表明所提出的方法具有一定的抗噪性。

实际应用：实际数据来自于中国西部某工区，用于验证提出方法的有效性。在进行反演之前需要对方位叠前道集进行保幅处理，然后将方位偏移距道集转化为方位入射角道集，经过处理后每个CDP道集的入射角范围为5°～35°，最后共划分方位角6个，即0°(-15°～15°)，30°(15°～45°)，60°(45°～75°)，90°(75°～105°)，120°(105°～135°)及150°(135°～165°)，入射角3个，即10°(5°～15°)，20°(15°～25°)及30°(25°～35°)，共18个方位部分角度叠加数据体。首先利用18个方位部分角度叠加数据体，通过约束稀疏脉冲反演得到18个方位弹性阻抗数据体，然后提取三个角度的二阶及四阶项傅里叶系数，最后利用提取的傅里叶系数估测得到裂缝弱度参数。图8(a)-图8(c)展示了方位部分角度叠加地震剖面。图9(a)-图9(c)展示了反演得到的方位弹性阻抗。图10展示了法向裂缝弱度和切向裂缝弱度的反演结果。红色表示高值，代表裂缝发育区域。从图10的裂缝弱度反演结果可以看出，在井位置处的法向裂缝弱度和切向裂缝弱度表现出高值，表明该位置处裂缝较为发育，这与井上的裂缝解释结果一致，进一步验证了我们所提出的方法的有效性。

本公开提出了一种新的方位傅里叶系数弹性阻抗反演方法，用于实现裂缝弱度参数预测。首先利用方位部分叠加地震数据，通过约束稀疏脉冲反演得到方位弹性阻抗，然后提取每个入射角对应的傅里叶系数，最后利用提取的二阶项及四阶项傅里叶系数进一步反演得到裂缝法向弱度和裂缝切向弱度，用于表征裂缝发育强度。模型测试和实际应用验证了所提出的方法的有效性。

实施例子二

包括：

还系统中模块的具体实现方式参见实施例子一中的基于方位傅里叶系数的弹性阻抗反演方法，此处不再进行具体说明。

实施例子三

本说明书实施方式提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现实施例子一中的基于方位傅里叶系数的弹性阻抗反演方法的步骤。

实施例子四

本说明书实施方式提供一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现实施例子一中的基于方位傅里叶系数的弹性阻抗反演方法的步骤。

可以理解的是，在本说明书的描述中，参考术语“一实施例”、“另一实施例”、“其他实施例”、或“第一实施例～第N实施例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims

1.一种基于方位傅里叶系数的弹性阻抗反演方法，其特征是，包括：

基于方位弹性阻抗方程获得各阶傅里叶系数与弹性参数及裂缝弱度参数之间的关系；

2.如权利要求1所述的一种基于方位傅里叶系数的弹性阻抗反演方法，其特征是，对于方位部分叠加地震数据具体为：每个入射角对应的方位地震子波以及通过裂缝岩石物理建模估测的法向裂缝弱度和切向裂缝弱度。

3.如权利要求1所述的一种基于方位傅里叶系数的弹性阻抗反演方法，其特征是，反演得到裂缝法向弱度和裂缝切向弱度时，在贝叶斯反演框架下，利用提取的傅里叶系数，通过柯西约束正则化和低频模型约束正则化的迭代重加权最小二乘法反演得到裂缝弱度参数。

4.如权利要求1所述的一种基于方位傅里叶系数的弹性阻抗反演方法，其特征是，将方位弹性阻抗方程两边同时取对数，进行傅里叶级数展开，对于X个规则采样的方位地震数据，傅里叶系数可根据离散傅里叶变换计算得到。

5.如权利要求1所述的一种基于方位傅里叶系数的弹性阻抗反演方法，其特征是，考虑时间采样点为M个，入射角为N个，基于各阶傅里叶系数与弹性参数及裂缝弱度参数之间的关系获得矩阵表达式，基于该矩阵表达式获得待反演模型参数的后验概率密度函数；

假设似然函数p(d|m)服从于高斯分布；

假设模型参数先验概率分布服从柯西分布；

求解得到最大后验概率初始目标函数；

6.一种基于方位傅里叶系数的弹性阻抗反演系统，其特征是，包括：

裂缝发育强度表征模块，被配置为：利用反演得到的所有方位弹性阻抗提取不同入射角下的傅里叶系数，利用提取的二阶项及四阶项傅里叶系数进一步反演得到裂缝法向弱度和裂缝切向弱度，用于表征裂缝发育强度；

基于方位弹性阻抗方程获得各阶傅里叶系数与弹性参数及裂缝弱度参数之间的关系。

7.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现权利要求1-5任一所述的基于方位傅里叶系数的弹性阻抗反演方法的步骤。

8.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现权利要求1-5任一所述的基于方位傅里叶系数的弹性阻抗反演方法的步骤。