CN109143328A - 一种叠后地震反演方法 - Google Patents

Abstract

本发明提出了一种叠后地震反演方法，包括如下步骤：S1、构建地震数据拟合差项；S2、构建模型趋势约束项；S3、结合模型趋势约束项和地震数据拟合差项，建立基追踪反演的目标函数；S4、通过基追踪分解算法从反演目标函数中获取波阻抗。本发明克服了常用的地震反演方法无法补偿工区的低频成分，导致反演结果的横向连续性较差以及反演结果缺失有效的背景趋势，使其在后续的储层预测中无法有效展开的缺陷。

Description

一种叠后地震反演方法

技术领域

本发明涉及油气勘探中地震资料解释与反演技术领域，尤其是涉及一种叠后地震反演方法。

背景技术

地震反演技术可以将反射界面的信息转换为储层内部的信息，因此可以更加直观的揭示储层的物性和含油气性。然而由于地震数据的带限特性，地震反演是个“病态问题”，通过加入合适的约束项，可以有效改善稳定性。

常规地震反演通常是从地震分界面的反射系数出发，在Lp范数最小化的约束下或者反射系数服从“长尾巴”分布的约束下，对目标函数进行最优化求解。但是常规地震反演无法突破1/4波长的极限分辨率，无法对子波进行充分压缩获取高分辨率的反演结果。Castagna和Chopra通过对楔形模型进行研究，提出了反射系数奇偶分解理论，将对反射系数为研究对象的传统思路转化为对薄层厚度和薄层位置的研究。Partyka、Portniaguine、Puryear、Chopra等借助谱分解和反射系数奇偶分解构建了谱反演的目标函数，通过模拟退火实现了高分辨率反演。

随着信号稀疏表示理论和压缩感知理论在地震勘探领域中应用，Zhang R通过反射系数奇偶分解理论和褶积模型理论构建了冗余的楔形字典，将地震反演问题转化为标准的基追踪问题，通过基追踪分解，可以得到一组最稀疏的解，进而可以获取高分辨率的反演结果。然而该方法无法补偿工区的低频成分，导致反演结果的横向连续性较差以及反演结果缺失有效的背景趋势，使其在后续的储层预测中无法有效展开。

发明内容

针对现有技术中所存在的上述技术问题，本发明提出了一种叠后地震反演方法。本发明通过结合反射系数奇偶分解理论、褶积模型理论以及基追踪分解理论，能够反演出高频的反射系数，最后结合模型趋势约束项，可以从初始波阻抗模型中补偿低频信息，进而可以获取横向上较为连续的全频带波阻抗反演结果。克服了常用的地震反演方法无法补偿工区的低频成分，导致反演结果的横向连续性较差以及反演结果缺失有效的背景趋势，使其在后续的储层预测中无法有效展开的缺陷。

本申请提出了一种叠后地震反演方法，包括如下步骤：

S1、构建地震数据拟合差项；

S2、构建模型趋势约束项；

S3、结合模型趋势约束项和地震数据拟合差项，建立基追踪反演的目标函数；

S4、通过基追踪分解算法从反演目标函数中获取波阻抗。

作为对本方法的进一步改进，所述步骤S1具体包括：

S11、根据地震波速度和地震波主频获取薄层的最大时间厚度，并根据薄层的最大时间厚度获取反射系数分解矩阵。其中，薄层的最大时间厚度对应的采样点个数决定了楔形字典以及反演参数的维数。

具体的，依据反射系数奇偶分解理论，一个薄层可以用一组反射系数对表示，其上、下两个反射系数分别表示薄层的项和底，反射系数之间的时间差则表示薄层的时间厚度，该反射系数对可以分解为一个反射系数偶分量和一个反射系数奇分量。

r_e(t，i，j，Δt)＝δ(t-iΔt)+δ(t-iΔt-jΔt) (5)

r_o(t，i，j，Δt)＝δ(t-iΔt)-δ(t-iΔt-jΔt) (6)

公式(5)和公式(6)中δ(*)表示狄拉克函数，t-iΔt对应于薄层顶部的时间，t-iΔt-jΔt对应于薄层底部的时间，薄层的时间厚度为jΔt；r_e和r_o分别表示反射系数偶分量和奇分量。可以看出，偶分量是由符号相同，幅值为1的两组脉冲构成；而奇分量是由符号相反，幅值为1的两组脉冲构成。单个薄层情况下，反射系数奇偶分解的数学表达为：

r(t)＝ar_e+br_o (7)

公式(7)中r(t)表示反射系数，a和b分别表示反射系数偶分量和奇分量的系数，假设反射系数的采样个数为N，薄层可能的最大时间厚度为MΔt，则多个薄层情况下，反射系数奇偶分解可以表示为：

公式(8)可以表达为矩阵形式：

r＝Dx (9)

公式(9)中，x表示反射系数奇偶分量系数：

x＝[a b]^T (10)

D表示反射系数分解矩阵，其数学表达式为：

D＝[D_e D₀] (13)

利用公式(9)～(15)可以建立起反射系数与薄层位置和薄层厚度的数学关系。

S12、根据地震子波构建地震子波褶积矩阵，将所述反射系数分解矩阵结合地震子波褶积矩阵构建楔形字典。

具体的，结合测井数据和井旁地震道，提取混合相位的地震子波w(t)。结合地震褶积模型理论以及S21步骤中获取的反射系数奇偶分解矩阵D，建立地震楔形字典G。

根据地震褶积模型理论，地震数据、地震子波以及反射系数满足如下数学关系：

s(t)＝w(t)*r(t)+n(t) (16)

公式(16)中*表示褶积运算，s(t)表示地震数据，w(t)表示地震子波，n(t)表示地震噪声。公式(16)同样可以用矩阵形式表示：

s＝wr+n (17)

结合公式(9)和公式(17)，地震数据的正演过程可以表示为：

s＝wDx+n (18)

其中，G＝wD，G表示楔形字典，该楔形字典建立了地震反射振幅与薄层厚度和薄层位置的关系。

S13、结合所述楔形字典和地震数据构建地震数据拟合差项。

公式(18)可以表示为：

s＝Gx+n (19)

假设随机噪声为高斯分布，地震数据的拟合差项可以表示为：||s-Gx||。

作为对本方法的进一步改进，所述步骤S2具体包括：

S21、测井数据中的波阻抗曲线进行环境校正和一致性校正后，在地震层位控制下进行横向插值，建立波阻抗模型。

S22、对地震数据进行频谱分析，获取地震数据中缺失的低频成分和截止频率f_low；根据所述截止频率f_low对所述波阻抗模型的自然对数进行低通滤波，获取波阻抗的自然对数的低频趋势I_Trend。

S23、根据所述截止频率f_low构建低通滤波矩阵L。

L＝F^-1ΛF (20)

公式(20)中F和F^-1分别为离散傅立叶变换矩阵，其表达式为：

其中，人为由汉宁窗函数组成的对角阵：

对角线元素Λ_fi表示汉宁窗函数表示的滤波器：

S24、根据所述低通滤波矩阵和波阻抗自然对数的低频趋势构建模型趋势约束项：

Constr_{mod el_trend}＝||LCDx-I_Trend|| (25)

该模型趋势约束项可以为基追踪反演提供量化的低频趋势约束和低频成分的精准补偿，有效提高反演结果的横向连续性。

作为对该方法的进一步改进，步骤S3中，再次参考图2，结合地震数据拟合差项以及模型趋势约束项，在最小化反演参数L1范数的约束下，建立反演目标函数为：

arg min{||s-Gx||+λ|x|+α||LCDx-I_Tremd||}，s.t.x_i≥0 (26)

但是，公式(26)中反演参数的维数大于观测数据的维数，属于欠定方程，不是基追踪问题的标准形式，将其转化为基追踪问题的标准形式：

其中，G表示楔形字典，L表示低通滤波矩阵，D表示反射系数奇偶分解矩阵，C表示积分矩阵，α表示稀疏约束项权重，s表示地震数据，I_Trend表示波阻抗自然对数的低频趋势，x表示反射系数奇偶分量系数，λ表示地震波长。

作为对该方法的进一步改进，所述步骤S4具体包括：

S41、获取待反演的地震数据，设置稀疏约束项权重；

S42、对反演目标函数进行基追踪分解获取反射系数奇偶分量系数；

S43、将获取的反射系数奇偶分量系数转换为反射系数并进行道积分获取波阻抗。

作为对该方法的进一步改进，所述步骤S43采用以下公式获取波阻抗：

其中，r＝Dx，r表示反射系数，D表示反射系数奇偶分解矩阵，x表示反射系数奇偶分量系数，I₀表示起始波阻抗。

作为对该方法的进一步改进，步骤S4中所述的基追踪分解算法采用原始-对偶对数障碍法或者梯度投影稀疏重构法。

与现有技术相比，本发明通过结合反射系数奇偶分解理论、褶积模型理论以及基追踪分解理论，能够反演出高频的反射系数，最后结合模型趋势约束项，可以从初始波阻抗模型中补偿低频信息，进而可以获取横向上较为连续的全频带波阻抗反演结果。克服了常用的地震反演方法无法补偿工区的低频成分，导致反演结果的横向连续性较差以及反演结果缺失有效的背景趋势，使其在后续的储层预测中无法有效展开的缺陷。

附图说明

下面将结合附图来对本发明的优选实施例进行详细地描述。在图中：

图1显示了根据本发明的实施例所述的叠后地震反演方法的流程图一。

图2显示了根据本发明的实施例所述的叠后地震反演方法的流程图二。

图3显示了一种理论地震数据。

图4(a)显示了采用本发明的实施例所述的叠后地震反演方法对图3地震理论数据进行反演的结果示意图。

图4(b)显示了采用常规叠后地震反演方法对图3地震理论数据进行反演的结果示意图。

图5(a)显示了采用本发明的实施例所述的叠后地震反演方法对图3理论地震数据进行反演获取的反射系数与实测井曲线结果对比示意图。

图5(b)显示了采用常规叠后地震反演方法对图3理论地震数据进行反演获取的反射系数与实测井曲线结果示意图。

图6(a)显示了一种实际地震数据。

图6(b)显示了采用本发明的实施例所述的叠后地震反演方法对图6(a)实际地震数据进行反演获取的反射系数结果示意图。

图7(a)显示了采用本发明的实施例所述的叠后地震反演方法对图6(a)实际地震数据进行反演获取的波阻抗剖面示意图。

图7(b)显示了采用常规叠后地震反演方法对图6(a)实际地震数据进行反演获取的波阻抗剖面示意图。

图8(a)显示了采用本发明的实施例所述的叠后地震反演方法对6(a)中共深度点道集75处进行反演获取的井旁道反演结果与实际井曲线剖面对比示意图。

图8(b)显示了采用常规叠后地震反演方法对对6(a)中共深度点道集75处进行反演获取的井旁道反演结果与实际井曲线剖面对比示意图。

在附图中，相同的部件使用相同的附图标记。附图并未按照实际的比例绘制。

具体实施方式

下面将结合附图对本发明做进一步说明。

参考图1，本申请提出了一种叠后地震反演方法，包括如下步骤：

S1、构建地震数据拟合差项；

S2、构建模型趋势约束项；

S3、结合模型趋势约束项和地震数据拟合差项，建立基追踪反演的目标函数；

S4、通过基追踪分解算法从反演目标函数中获取波阻抗。

作为对本方法的进一步改进，参考图2，所述步骤S1具体包括：

S11、根据地震波速度和地震波主频获取薄层的最大时间厚度，并根据薄层的最大时间厚度获取反射系数分解矩阵。其中，薄层的最大时间厚度对应的采样点个数决定了楔形字典以及反演参数的维数。

具体的，依据反射系数奇偶分解理论，一个薄层可以用一组反射系数对表示，其上、下两个反射系数分别表示薄层的顶和底，反射系数之间的时间差则表示薄层的时间厚度，该反射系数对可以分解为一个反射系数偶分量和一个反射系数奇分量。

r_e(t，i，j，Δt)＝δ(t-iΔt)+δ(t-iΔt-jΔt) (5)

r_o(t，i，j，Δt)＝δ(t-iΔt)-δ(t-iΔt-jΔt) (6)

公式(5)和公式(6)中δ(*)表示狄拉克函数，t-iΔt对应于薄层顶部的时间，t-iΔt-jΔt对应于薄层底部的时间，薄层的时间厚度为jΔt；r_e和r_o分别表示反射系数偶分量和奇分量。可以看出，偶分量是由符号相同，幅值为1的两组脉冲构成；而奇分量是由符号相反，幅值为1的两组脉冲构成。单个薄层情况下，反射系数奇偶分解的数学表达为：

r(t)＝ar_e+br_o (7)

公式(7)中r(t)表示反射系数，a和b分别表示反射系数偶分量和奇分量的系数，假设反射系数的采样个数为N，薄层可能的最大时间厚度为MΔt，则多个薄层情况下，反射系数奇偶分解可以表示为：

公式(8)可以表达为矩阵形式：

r＝Dx (9)

公式(9)中，x表示反射系数奇偶分量系数：

x＝[a b]^T (10)

D表示反射系数分解矩阵，其数学表达式为：

D＝[D_e D₀] (13)

利用公式(9)～(15)可以建立起反射系数与薄层位置和薄层厚度的数学关系。

S12、根据地震子波构建地震子波褶积矩阵，将所述反射系数分解矩阵结合地震子波褶积矩阵构建楔形字典。

具体的，结合测井数据和井旁地震道，提取混合相位的地震子波w(t)。结合地震褶积模型理论以及S21步骤中获取的反射系数奇偶分解矩阵D，建立地震楔形字典G。

根据地震褶积模型理论，地震数据、地震子波以及反射系数满足如下数学关系：

s(t)＝w(t)*r(t)+n(t) (16)

公式(16)中*表示褶积运算，s(t)表示地震数据，w(t)表示地震子波，n(t)表示地震噪声。公式(16)同样可以用矩阵形式表示：

s＝wr+n (17)

结合公式(9)和公式(17)，地震数据的正演过程可以表示为：

s＝wDx+n (18)

其中，G＝wD，G表示楔形字典，该楔形字典建立了地震反射振幅与薄层厚度和薄层位置的关系。

S13、结合所述楔形字典和地震数据构建地震数据拟合差项。

公式(18)可以表示为：

s＝Gx+n (19)

假设随机噪声为高斯分布，地震数据的拟合差项可以表示为：||s-Gx||。

作为对本方法的进一步改进，所述步骤S2具体包括：

S21、测井数据中的波阻抗曲线进行环境校正和一致性校正后，在地震层位控制下进行横向插值，建立波阻抗模型。

S22、对地震数据进行频谱分析，获取地震数据中缺失的低频成分和截止频率f_low；根据所述截止频率f_low对所述波阻抗模型的自然对数进行低通滤波，获取波阻抗的自然对数的低频趋势I_Trend。

S23、根据所述截止频率f_low构建低通滤波矩阵L。

L＝F^-1ΛF (20)

公式(20)中F和F^-1分别为离散傅立叶变换矩阵，其表达式为：

其中，Λ为由汉宁窗函数组成的对角阵：

对角线元素Λ_fi表示汉宁窗函数表示的滤波器：

S24、根据所述低通滤波矩阵和波阻抗自然对数的低频趋势构建模型趋势约束项：

Constr_{mod el_trend}＝||LCDx-I_Trend|| (25)

该模型趋势约束项可以为基追踪反演提供量化的低频趋势约束和低频成分的精准补偿，有效提高反演结果的横向连续性。

作为对该方法的进一步改进，步骤S3中，再次参见图2，结合地震数据拟合差项以及模型趋势约束项，在最小化反演参数L1范数的约束下，建立反演目标函数为：

arg min{||s-Gx||+λ|x|+α||LCDx-I_Tremd||}，s.t.x_i≥0 (26)

但是，公式(26)中反演参数的维数大于观测数据的维数，属于欠定方程，不是基追踪问题的标准形式，将其转化为基追踪问题的标准形式：

其中，G表示楔形字典，L表示低通滤波矩阵，D表示反射系数奇偶分解矩阵，C表示积分矩阵，α表示稀疏约束项权重，s表示地震数据，I_Trend表示波阻抗自然对数的低频趋势，x表示反射系数奇偶分量系数，λ表示地震波长。

作为对该方法的进一步改进，所述步骤S4具体包括：

S41、获取待反演的地震数据，设置稀疏约束项权重；

S42、对反演目标函数进行基追踪分解获取反射系数奇偶分量系数；

S43、将获取的反射系数奇偶分量系数转换为反射系数并进行道积分获取波阻抗。

作为对该方法的进一步改讲.所述步骤S43采用以下公式获取波阻抗：

其中，r＝Dx，r表示反射系数，D表示反射系数奇偶分解矩阵，x表示反射系数奇偶分量系数，I₀表示起始波阻抗。

作为对该方法的进一步改进，步骤S4中所述的基追踪分解算法采用原始-对偶对数障碍法或者梯度投影稀疏重构法。

图3给出了一种理论地震数据。图3(a)为深度域的声波时差曲线和密度曲线，考虑到测井数据和地震数据尺度差异，借助于Backus滤波在深度域进行500Hz滤波，消除尺度差异，在此基础上，将其重采样至时间域后(如图3(b)所示)，计算反射系数，并于40Hz的雷克子波进行褶积得到合成地震记录(如图3(c)所示)。

图4(a)为基于图2中的理论数据，利用本发明获取的波阻抗结果与实测井曲线对比，其中实线为实测井曲线，点虚线为反演结果；图4(b)为基于图2中的理论数据，利用常规叠后反演获取的波阻抗与实测井曲线对比，其中实线为实测井曲线，点虚线为反演结果。通过对比可以看出，利用本发明的获取的波阻抗反演结果具有较高的垂向分辨率。

图5(a)为基于图2中的理论数据，利用本发明获取的反射系数结果与实测井曲线对比，其中实线为实测井曲线，点虚线为反演结果；图5(b)为基于图2中的理论数据，利用常规叠后反演获取的反射系数与实测井曲线对比，其中实线为实测井曲线，点虚线为反演结果。通过对比可以看出，利用本发明获取的波阻抗反演结果具有较高的垂向分辨率。

图6(a)为验证本发明所选用的实际地震数据，图6(b)是根据本发明获取的反射系数反演结果，通过对比可以看出，反射系数剖面垂向分辨率较高，子波旁瓣得到了有效压制，地层分界面较为清晰。

图7(a)为本发明获取的波阻抗剖面，图7(b)为常规叠后地震反演获取的波阻抗剖面。通过对比可以看出，相对常规反演方法，本发明所获取的波阻抗剖面垂向分辨率更高，展示了更加丰富的细节信息。

图8给出了井A位于CDP 75(共深度点道集75)处的井旁道反演结果与实际井曲线的对比，其中图8(a)为本发明获取的反演结果，图8(b)为常规叠后剖面获取的反演结果，图中实线为实测井曲线，点虚线为反演结果，通过对比可以看出，相对常规叠后反演，本发明获取的反演结果与实测井曲线吻合度更高，垂向分辨率也更高。

与现有技术相比，本发明通过结合反射系数奇偶分解理论、褶积模型理论以及基追踪分解理论，能够反演出高频的反射系数，最后结合模型趋势约束项，可以从初始波阻抗模型中补偿低频信息，进而可以获取横向上较为连续的全频带波阻抗反演结果。克服了常用的地震反演方法无法补偿工区的低频成分，导致反演结果的横向连续性较差以及反演结果缺失有效的背景趋势，使其在后续的储层预测中无法有效展开的缺陷。

以上所述仅为本发明的优选实施方式，但本发明保护范围并不局限于此，任何本领域的技术人员在本发明公开的技术范围内，可容易地进行改变或变化，而这种改变或变化都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求书的保护范围为准。

Claims (9)

1.一种叠后地震反演方法，其特征在于，包括如下步骤：

S1、构建地震数据拟合差项；

S2、构建模型趋势约束项；

S3、结合所述模型趋势约束项和所述地震数据拟合差项，建立基追踪反演的目标函数；

S4、通过基追踪分解算法从所述目标函数中获取波阻抗。

2.根据权利要求1所述的叠后地震反演方法，其特征在于，所述步骤S1具体包括：

S11、根据地震波速度和地震波主频获取薄层的最大时间厚度，并根据薄层的最大时间厚度获取反射系数分解矩阵；

S12、根据地震子波构建地震子波褶积矩阵，将所述反射系数分解矩阵结合地震子波褶积矩阵构建楔形字典；

S13、结合所述楔形字典和地震数据构建地震数据拟合差项。

3.根据权利要求1所述的叠后地震反演方法，其特征在于，所述步骤S2具体包括：

S21、建立波阻抗模型；

S22、根据所述波阻抗模型获取波阻抗自然对数的低频趋势；

S23、构建低通滤波矩阵；

S24、根据所述低通滤波矩阵和波阻抗自然对数的低频趋势构建模型趋势约束项。

4.根据权利要求3所述的叠后地震反演方法，其特征在于，所述步骤S21具体包括：将测井数据中的波阻抗曲线进行环境校正和一致性校正后，在地震层位控制下进行横向插值，建立波阻抗模型。

5.根据权利要求4所述的叠后地震反演方法，其特征在于，所述步骤S22具体包括：

S221、对地震数据进行频谱分析，获取地震数据中缺失的低频成分和截止频率；

S222、根据所述截止频率对所述波阻抗模型的自然对数进行低通滤波，获取波阻抗的自然对数的低频趋势。

6.根据权利要求1所述的叠后地震反演方法，其特征在于，所述步骤S3中所述反演目标函数为：

其中，G表示楔形字典，L表示低通滤波矩阵，D表示反射系数奇偶分解矩阵，C表示积分矩阵，α表示稀疏约束项权重，s表示地震数据，I_trend表示波阻抗自然对数的低频趋势，x表示反射系数奇偶分量系数，λ表示地震波长。

7.根据权利要求6所述的叠后地震反演方法，其特征在于，所述步骤S4具体包括：

S41、获取待反演的地震数据，设置稀疏约束项权重；

S42、对反演目标函数进行基追踪分解获取反射系数奇偶分量系数；

S43、将获取的反射系数奇偶分量系数转换为反射系数并进行道积分获取波阻抗。

8.根据权利要求7所述的叠后地震反演方法，其特征在于，所述步骤S43采用以下

公式获取波阻抗：

其中，r＝Dx，r表示反射系数，D表示反射系数奇偶分解矩阵，x表示反射系数奇偶分量系数，I₀表示波阻抗起始值。

9.根据权利要求1所述的叠后地震反演方法，其特征在于，步骤S4中所述的基追踪分解算法采用原始-对偶对数障碍法或者梯度投影稀疏重构法。