CN111175824B - 岩相驱动下的时频联合域地震反演方法 - Google Patents

Abstract

本发明提供一种岩相驱动下的时频联合域地震反演方法，包括：步骤1，根据联合域卷积模型和部分角度叠加地震道集的时‑频域响应，构建时频联合域地震时频联合域地震正演模型；步骤2，利用后验概率分布显式形式量化模型参数的最优解，推导出混合后验概率模型的显式解；步骤3，根据混合后验概率密度分布的均值和协方差的显式表达式，计算获得三类已知数据集

约束下的地震反演最优解；步骤4，利用储层测井、地质信息中模型参数中的低界和高界信息，进行基于边界约束的时频联合域地震反演，协同估算储层离散岩相。该岩相驱动下的时频联合域地震反演方法提高了模型参数反演的稳定性，时频联合域反演、序贯模拟采样策略改善了预测结果的分辨能力。

Description

岩相驱动下的时频联合域地震反演方法

技术领域

本发明涉及地震反演技术领域，特别是涉及到一种岩相驱动下的时频联合域地震反演方法。

背景技术

地震反演是获取地下复杂介质模型参数、地层岩性信息最主要的定量评价手段。地层弹性参数与岩性之间往往存在着确定性的联系，后验概率密度分布的随机解是叠前地震概率化反演中不确定性分析、可靠性评价的重要依据。基于贝叶斯推断的最大后验概率解是概率化反演最常见的方法，其引入了待估计模型参数的先验统计分布信息(即，对未知模型的经验认知)，对观测地震数据的不足起到了补充作用，或可以理解为对最大似然估计解的修正，对缓解地震反演的不适定性有重要意义(Bayes，1973)。由于稀疏概率密度分布(如，Cauchy、Laplace、指数分布)具有“长尾巴-集中于均值”的特征，即，假定地层反射系数集中分布于均值附近，距离均值越远、概率分布随模型参数的变化梯度逐渐平缓，从正则化的角度考虑，该分布相当于“变权重正则化”方法，其具有压制弱小反射(弱小噪声)、突出强反射率的稀疏反演特性(Downton等,2005；Buland等,2003；杨培杰等；Yuan等,2013；Yin等，2008&2016)。Downton等(2005)提出了利用Cauchy概率密度分布最大后验概率解的叠前地震AVO(振幅随偏移距或入射角变化)线性反演方法。Theune等(2010)实现了模型参数Cauchy、Laplace先验分布情况下，地层反射系数的反演，通过道积分获取绝对模型参数，表明柯西分布对提高稀疏反射率恢复、提高反演分辨率有一定的帮助。Alemie和Sacchi等[80](2011)、张世鑫等(2011)基于三变量Cauchy概率分布开展了的叠前地震三参数同步反演方法，都采用反复重加权最小二乘迭代方法求解目标泛函。Yang等(2008)、Zong等(2012)、Yin等(2014&2016)、刘晓晶等(2016)、Li K等(2017&2018)对柯西概率模型分布、修正柯西概率分布、双边指数概率分布或Laplace概率分布的先验假设情况下的最大后验概率解，作出详细的阐述和总结，主要采迭代复重加权最小二乘(IRLS)、和二次线性规划求解地震反问题(如，GPSR梯度投影算法)等。

上面所述常规地震反演方法往往将“地层岩性”、“模型参数”两者独立预测或间接预测，忽视了地层岩性对模型参数的影响，由此引入的先验知识误差会严重影响地震反演的精度，增加“连续弹性参数”预测、“离散岩相”判识的累计误差，导致岩性预测的精度降低。为此我们发明了一种新的岩相驱动下的时频联合域地震反演方法，解决了以上技术问题。

发明内容

本发明的目的是提供一种面向储层离散岩性、模型参数的，改善了预测结果的分辨能力的岩相驱动下的时频联合域地震反演方法。

本发明的目的可通过如下技术措施来实现：岩相驱动下的时频联合域地震反演方法，该岩相驱动下的时频联合域地震反演方法包括：步骤1，根据联合域卷积模型和部分角度叠加地震道集的时-频域响应，构建时频联合域地震时频联合域地震正演模型；步骤2，利用后验概率分布显式形式量化模型参数的最优解，推导出混合后验概率模型的显式解；步骤3，根据混合后验概率密度分布的均值和协方差的显式表达式，计算获得三类已知数据集H约束下的地震反演最优解；步骤4，利用储层测井、地质信息中模型参数中的低界和高界信息，进行基于边界约束的时频联合域地震反演，进而协同估算储层离散岩相。

本发明的目的还可通过如下技术措施来实现：

在步骤1中，首先将联合域卷积模型被视为正演算子，将部分角度叠加地震道集的时-频域响应相结合，构建了时频联合域地震时频联合域地震正演模型，将已知的纵波阻抗的低频先验信息引入正演方程中，以补偿观测地震数据中缺失的低频分量。

在步骤1中，将联合域卷积模型视为正演算子，

式中，S(f)、W(f)分别为入射角为θ_i时的地震频谱和子波频谱，r_pp(τ_z)为Russell线性近似AVO方程，f为频率，τ_z为时间序列，N(f)为高斯随机噪声的频谱向量；为了更清楚的阐述地震反问题，方程(1)重组为，

S^f＝W^f·E·C·m+N_f (2)

式中，W^f＝diag[W(f)]为角度地震子波的对角线频谱矩阵，C为一阶差分矩阵，E代表Fourier变换算子或相移算子，

为纵波阻抗的自然对数形式；考虑到时间域AVO反演的强稳定性和频率域AVO反演的高分辨率，将部分角度叠加地震道集的时-频域响应相结合，构造了时频联合域正演算子，旨在稳定性和分辨率达到平衡，

式中，S^t、W^t和N_t分别是时间域部分角度叠加地震数据、时间域子波矩阵和高斯随机噪声，λ_t和λ_f表示时频域地震反演的加权系数；在不同信噪比SNRs的情况下，通过调整加权系数(λ_t，λ_f)、选择高信噪比的频率分量，进而控制地震反演的稳定性和分辨率；当SNRs高的情况下，提高频率域反演的权重λ_f，提高反演的分辨率；在SNRs较低的情况下，增加时域权重λ_t，选择信噪比较高的频率分量，进而提高预测的稳定性；将已知的纵波阻抗的低频先验信息引入正演方程(3)中，以补偿观测地震数据中缺失的低频分量，其罚函数J_lo_w(m)即通过L2范数

引入目标函数中；η为通过实际测井数据的局部克里格插值构建的纵波阻抗的低频先验，λ_L为低频背景的正则化因子，D为正则化对角线矩阵；联合方程(3)和低频先验背景得到方程(4)，

式中，式中，H为输入的三种类型已知数据集(时间域地震、频率域地震及低频先验数据)，P为正演核矩阵，C_ff、C_tt、C_ηη、C_ft、C_fη及C_tη是3种已知数据集(频率域地震数据S^f、时间域地震数据S^t和纵波阻抗的低频背景η)的先验协方差矩阵，C_H为三种数据集H的综合先验协方差矩阵，λ_t和λ_f表示时间域、频率域地震数据的加权系数，W^f为角度地震子波的对角线频谱矩阵，C为一阶差分矩阵，E代表Fourier变换算子或相移算子，W^t是时间域子波矩阵，D为低频正则化的对角线矩阵。

步骤2中，考虑到不同岩相下纵波阻抗的先验概率密度函数存在差异，将多个高斯分量的混合概率密度分布作为待反演模型的先验概率密度分布；假设目标工区中存在的储层岩性类型数量为M，针对线性地震反问题的贝叶斯推理，则模型参数的最优解利用后验概率分布显式形式进行量化，推导得到三类已知数据集(S^f、S^t和η)约束下的第k_th高斯分量p_k(m|H)的显式表达式

为，

另外，协方差矩阵

为：

式中，

式中，

和k分别为待反演模型参数m的先验均值、先验协方差及第k_th个高斯分量的上标；S_ff、S_tt、S_ηη、S_ft、S_fη及S_tη表示已知模型参数m情况下正演观测数据(频率域地震数据S^f、时间域地震数据S^t和低频背景η)的协方差矩阵。H为输入的三种类型已知数据集(时间域地震、频率域地震及低频先验数据)，P为正演核矩阵，C_ff、C_tt、C_ηη、C_ft、C_fη及C_tη是3种已知数据集(频率域地震数据S^f、时间域地震数据S^t和纵波阻抗的低频背景η)的先验协方差矩阵，C_H为三种数据集H的综合先验协方差矩阵，λ_t和λ_f表示时间域、频率域地震数据的加权系数，W^f为角度地震子波的对角线频谱矩阵，C为一阶差分矩阵，E代表Fourier变换算子或相移算子，W^t是时间域子波矩阵，D为低频正则化的对角线矩阵。此外，待反演模型参数的混合后验高斯概率密度分布p(m|H)重新表征为一个混合概率模型

其中不同岩相对应的高斯分量后验权重

等于

λ_k为第k_th个高斯分量的先验权重，

分别为已知数据集H的先验均值和先验协方差，M为储层岩性的分类数量，N_l(·)表示第l个高斯概率模型,N_k(·)表示第k个高斯概率模型

在步骤3中，通过混合后验概率密度分布的均值和协方差的显式表达式，计算获得三类已知数据集H约束下的地震反演最优解；将实际测井数据和先前的模拟点作为地震反演的约束条件，并基于单点实现的序列模拟算法构建离散相序贯采样算法流程。

在步骤3中，针对第k_th个高斯分量p_k(m|H)，在待模拟点m_i位置的后验均值和协方差改写为：

式中，

分别为待模拟点m_i位置的后验均值和先验均值，

为第k个高斯概率分量情况下4种类型(m_s、S^f、S^t及η)条件数据集的协方差矩阵，m_s代表已知数据点和先前模拟点；

为第k个高斯概率分量的初始模型，U_i＝[0 L 0 1 L 0]为待模拟点m_i位置的采样矩阵即，第ith个元素为1，其余元素均为0，U为已知或先前模拟数据m_s的提取矩阵，S^f为频率域地震数据，S^t为时间域地震数据，η为低频背景。λ_t和λ_f表示时间域、频率域地震数据的加权系数，λ_L为低频背景的正则化系数，

为先前模拟数据m_s的正则化系数，W^f为角度地震子波的对角线频谱矩阵，C为一阶差分矩阵，E代表Fourier变换算子或相移算子，W^t是时间域子波矩阵，D为低频正则化的对角线矩阵，

为待反演模型参数m的先验协方差。当λ_f＝0、λ_t＝0和λ_L＝0的情况下，方程(8)等价于kriging插值和地质统计学中的序贯高斯模拟，

表示第四类已知数据集m_s的权重系数；此外，待模拟点m_i的协方差

改写为：

式中，U_i＝[0 L 0 1 L 0]为待模拟点m_i位置的采样矩阵即，第ith个元素为1，其余元素均为0，已知或先前模拟数据m_s的提取矩阵为U＝[U_i；U_j；U_l；L；U_e]^t，

为四种类型(m_s、S^f、S^t及η)条件数据集的协方差矩阵，m_s代表已知数据点和先前模拟点；

为第k个高斯概率分量的初始模型，U_i＝[0 L 0 1 L 0]为待模拟点m_i位置的采样矩阵即，第ith个元素为1，其余元素均为0，U为已知或先前模拟数据m_s的提取矩阵，S^f为频率域地震数据，S^t为时间域地震数据，η为低频背景。λ_t和λ_f表示时间域、频率域地震数据的加权系数，λ_L为低频背景的正则化系数，

为先前模拟数据m_s的正则化系数，W^f为角度地震子波的对角线频谱矩阵，C为一阶差分矩阵，E代表Fourier变换算子或相移算子，W^t是时间域子波矩阵，D为低频正则化的对角线矩阵，

为待反演模型参数m的先验协方差；令F＝W^f·E·C，G＝W^t·C，

的显式表达式简化为，

根据基于贝叶斯分类，推导得到四类已知数据作为约束条件下第k_th个高斯分量的后验权重

更新为：

式中，

为序贯模拟过程中，四类条件数据的先验均值，

为四种类型(m_s、S^f、S^t及η)条件数据集的协方差矩阵，m_s代表已知数据点和先前模拟点；

为第k个高斯概率分量的初始模型，S^f为频率域数据，S^t为时间域数据，η为低频背景，N_l(·)表示第l个高斯概率模型，N_k(·)表示第k个高斯概率模型。

在步骤3中，方程(11)中，待模拟点m_i的后验权重

用来判识离散岩相；当离散相态分类即优选混合概率模型中的高斯分量之后，纵波阻抗即连续型模型参数从被选择的高斯分量中随机采样获得；当访问下一个模拟点时，先前的模拟点将被视为已知数据，并且已知数据的邻域选择对于提高地震序贯模拟的计算效率具有重要作用，选用矩形邻域或圆形邻域，邻域半径通常设置为一个波长至两个波长之间(λ＜r＜2λ)。

在步骤4中，由于方程(11)中并未考虑待反演弹性参数的下限和上限，由于实际地震数据受噪声干扰，反演结果往往出现一些超边界的解，非线性边界转化算法为

式中，x、m^min、m^max及N分别是待反演模型参数的非线性中间变量、下界先验、上界先验和质控常数；令更新后的弹性参数Δm＝Δm^k＝m^k-m^k-1，则第kth次迭代后的相对弹性参数更新为：

式中，m^k-1表示第(k-1)_th次迭代后的估计结果，式(13)中，N主要用于控制弹性参数的梯度

N的范围为(0.1,2.0]，控制超过边界的预测结果m^min和m^max的更新幅度。

本发明中的岩相驱动下的时频联合域地震反演方法，考虑时间域地震反演和频率域地震反演在分辨率和抗噪性方面独有的特征，构建了由时域地震、频域地震、低频整合先验、及已知模型数据四类观测数据集协同约束的混合后验概率模型(GMM)的显式解，利用非线性转换函数的形式将模型参数的低界和高界约束信息引入到地震纵波阻抗反演中，提出了序贯模拟随机采样的岩相驱动时频联合域地震概率化反演方法，实现了对地层“模型参数”、“离散岩性”的时-频联合域地震同步预测，整合先验和边界约束算法的引入提高了模型参数反演的稳定性，时频联合域反演、序贯模拟采样策略改善了预测结果的分辨能力。

附图说明

图1为本发明的一具体实施例中无噪声干扰情况下理论测试，纵波阻抗、离散岩相的协同估算结果的示意图；

图2为本发明的一具体实施例中信噪比SNR＝5情况下理论测试，纵波阻抗、离散岩相的协同估算结果的示意图；

图3为本发明的一具体实施例中信噪比SNR＝2情况下理论测试，纵波阻抗、离散岩相的协同估算结果的示意图；

图4为本发明的一具体实施例中信噪比SNR＝1情况下理论测试，纵波阻抗、离散岩相的协同估算结果的示意图；

图5为本发明的一具体实施例中中国东部勘探工区的某条地震剖面的示意图；

图6为本发明的一具体实施例中常规时间域地震反演的纵波阻抗剖面的示意图；

图7为本发明的一具体实施例中时频联合域地震反演的纵波阻抗剖面的示意图；

图8为本发明的一具体实施例中岩相驱动下的时频联合域地震反演的离散岩相协同估算结果的示意图；

图9为本发明的一具体实施例中10次模拟的砂岩相预测概率的示意图；

图10为本发明的岩相驱动下的时频联合域地震反演方法的一具体实施例的流程图。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。

如图10所示，图10为本发明的岩相驱动下的时频联合域地震反演方法的流程图。

本发明的岩相驱动下的时频联合域地震反演方法，首先(1)构建了时频联合域地震数据协同约束的最大后验概率解(MAP)解。推导了由时-频域地震、低频整合先验及已知模型数据点四类条件数据集协同约束的混合后验概率分布的显式解。

(2)提出了基于序贯采样的地层“模型参数”、“离散岩性”的时-频联合域高分辨协同判识。通过混合后验概率分布的均值和协方差的显式表达式，即可以利用序贯模拟方法实现混合后验概率密度分布的随机采样，计算获得四类已知数据集协同约束下的地震反演最优解。

(3)探索基于非线性边界约束的岩相驱动时频联合域地震反演，利用非线性转换函数的形式将模型参数的低界和高界约束信息引入到时频联合域地震反演中，边界约束策略可以有效的将待反演参数严格控制在模型低界和高界先验信息之内，将随机采样过程中出现的“小概率解”进行修正，避免了不现实解的产生。

步骤101，首先将联合域卷积模型被视为正演算子，将部分角度叠加地震道集的时-频域响应相结合，构建了时频联合域地震时频联合域地震正演模型，将已知的纵波阻抗的低频先验信息引入正演方程(3)中，以补偿观测地震数据中缺失的低频分量。

不同变换域地震反演方法是考虑不同域地震响应与模型合成地震数据间的匹配程度来搜索最优模型参数的过程。频率域反演是在频率域惩罚函数的约束下，发展的一种仅仅采用部分频谱分量去预测模型参数的方法。频率域反演相比常规时间域反演具有频率分量自动解耦、频率选择自由、多分量逐级寻优、及分辨率高的特点，然而，频率域反演的稳定性比时间域反演低，因此，将联合域卷积模型被视为正演算子，

式中，S(f)、W(f)分别为入射角为θ_i时的地震频谱和子波频谱，r_pp(τ_z)为Russell线性近似AVO方程，f为频率，τ_z为时间序列，N(f)为高斯随机噪声的频谱向量。为了更清楚的阐述地震反问题，方程(1)可以重组为，

s^f＝W·E·C·m+N_f (2)

式中，W^f＝diag[W(f)]为角度地震子波的对角线频谱矩阵，C为一阶差分矩阵，E代表Fourier变换算子或相移算子，

为纵波阻抗的自然对数形式。考虑到时间域AVO反演的强稳定性和频率域AVO反演的高分辨率，将部分角度叠加地震道集的时-频域响应相结合，构造了时频联合域正演算子，旨在稳定性和分辨率达到平衡，

式中，S^t、W^t和N_t分别是时间域部分角度叠加地震数据、时间域子波矩阵和高斯随机噪声，λ_t和λ_f表示时频域地震反演的加权系数。在不同信噪比(SNRs)的情况下，通过调整加权系数(λ_t，λ_f)、选择高信噪比的频率分量，进而控制地震反演的稳定性和分辨率。当SNRs高的情况下，可以提高频率域反演的权重λ_f，提高反演的分辨率。在SNRs较低的情况下，增加时域权重λ_t，选择信噪比较高的频率分量，进而提高预测的稳定性。将已知的纵波阻抗的低频先验信息引入正演方程(3)中，以补偿观测地震数据中缺失的低频分量，其罚函数J_low(m)即通过L2范数

引入目标函数中。η为通过实际测井数据的局部克里格插值构建的纵波阻抗的低频先验(＜12Hz)，λ_L为低频背景的正则化因子，D为正则化对角线矩阵。联合方程(3)和低频先验背景得到方程(4)，

式中，H为输入的三种类型已知数据集，P为正演核矩阵，C_ff、C_tt、及C_ηη是三种已知数据集(S^f、S^t和η)的协方差矩阵，C_H为三种数据集H的联合协方差矩阵。

在步骤102中，利用后验概率分布显式形式量化模型参数的最优解，推导出混合后验概率模型的显式解。

考虑到不同岩相下纵波阻抗的先验概率密度函数存在差异，本发明将多个高斯分量的混合概率密度分布作为待反演模型的先验概率密度分布。假设目标工区中存在的储层岩性类型数量为M，针对线性地震反问题的贝叶斯推理，则模型参数的最优解可以利用后验概率分布显式形式进行量化，推导得到三类已知数据集(S^f、S^t和η)约束下的第k_th高斯分量p_k(m|H)的显式表达式

为，

另外，协方差矩阵

为：

式中，

式中，

和k分别为待反演模型参数m的先验均值、先验协方差及第k_th个高斯分量的上标。此外，待反演模型参数的混合后验高斯概率密度分布p(m|H)可以重新表征为一个混合概率模型

其中不同岩相对应的高斯分量后验权重

等于

λ_k为第k_th个高斯分量的先验权重，

分别为三类已知数据集H的先验均值和先验协方差。

在步骤103中，通过混合后验概率密度分布的均值和协方差的显式表达式，计算获得三类已知数据集H约束下的地震反演最优解。将实际测井数据和先前的模拟点作为地震反演的约束条件，并基于单点实现的序列模拟算法构建离散相序贯采样算法流程。

针对第k_th个高斯分量p_k(m|H)，在待模拟点m_i位置的后验均值和协方差可以改写为：

式中，

分别为待模拟点m_i位置的后验均值和先验均值，m_s代表已知数据点和先前模拟点。当λ_f＝0、λ_t＝0和λ_L＝0的情况下，方程(8)等价于kriging插值和地质统计学中的序贯高斯模拟，

表示第四类已知数据集m_s的权重系数。此外，待模拟点m_i的协方差

可以改写为：

式中，U_i＝[0 L 0 1 L 0]为待模拟点m_i位置的采样矩阵(即，第ith个元素为1，其余元素均为0)，已知(或先前模拟)数据m_s的提取矩阵为U＝[U_i；U_j；U_l；L；U_e]^t，

为四种类型(m_s、S^f、S^t及η)条件数据集的协方差矩阵。令F＝W^f·E·C，G＝W^t·C，

的显式表达式可以简化为，

根据基于贝叶斯分类，推导得到四类已知数据作为约束条件下第k_th个高斯分量的后验权重可以更新为：

式中，

为序贯模拟过程中，四类条件数据的先验均值。方程(11)中，待模拟点m_i的后验权重

可以用来判识离散岩相。当离散相态分类(优选混合概率模型中的高斯分量)之后，纵波阻抗(连续型模型参数)从被选择的高斯分量中随机采样获得。值得注意的是，当我们访问下一个模拟点时，先前的模拟点将被视为已知数据，并且已知数据的邻域选择对于提高地震序贯模拟的计算效率具有重要作用，通常选用矩形邻域或圆形邻域，邻域半径通常设置为一个波长至两个波长之间(λ＜r＜2λ)。

在步骤104中，利用储层测井、地质等信息中模型参数中的低界和高界信息，实现基于边界约束的时频联合域地震反演方法。

由于方程(11)中并未考虑待反演弹性参数的下限和上限，由于实际地震数据受噪声干扰，反演结果往往出现一些超边界的解，非线性边界转化算法为

式中，x、m^min、m^max及N分别是待反演模型参数的非线性中间变量、下界先验、上界先验和质控常数；令，更新后的弹性参数Δm＝Δm^k＝m^k-m^k-1，则第kth次迭代后的相对弹性参数可更新为：

式中，m^k-1表示第(k-1)_th次迭代后的估计结果，式(13)中，N主要用于控制弹性参数的梯度

通常来讲，N的范围为(0.1,2.0]，可以控制超过边界的预测结果m^min和m^max的更新幅度。

为了进一步了解该算法的细节，验证该岩相驱动下的时频联合域地震反演方法的可行性，进行理论模型测试，图1—图4为通过30Hz零相位Ricker小波、2ms采样间隔和精确的Zoeppritz方程，合成不同信噪比(SNRs)的角度依赖地震道集，从图中可以看出，该方法反演的岩相最大条件概率解是合理的，与理论岩相具有高度一致性，从而表明岩相驱动下的时频联合域地震反演方法可以有效地实现弹性阻抗和离散岩相的稳定预测。

为验证本发明反演方法的实用性，本发明针对一具体实施例进行研究，图5-图9分别为部分角度叠加实际地震剖面、常规时间域、时频联合域地震反演的纵波阻抗剖面、离散岩相协同估算结果及10次实现的砂岩概率，图6和图7中曲线为实际测井曲线解释结果和地震反演估计结果，对比可以看出时频域联合反演结果与测井数据吻合性更好；图8和图9中黑线表示油砂层的解释结果，黑色椭圆指示油层发育位置。图中可见，砂泥岩模型的岩性判识结果与砂岩解释结果(黑线)高度一致。

本发明的组稀疏叠前四参数反演技术，综合考虑了待反演参数反射系数的向量稀疏性及结构稀疏性，并增加了与吸收衰减相关的地层品质因子参数的反演。其特征是：地层品质因子的倒数能够较真实的反应地震波能量衰减特征，对待反演参数进行分组，构建更符合实际情况的先验信息，同时进行多道反演，便于得到横向连续性更强、更稳定的反演结果。

考虑反射系数的内部结构特征，对待反演参数进行分组，在贝叶斯框架下，将分组结果代入修正柯西分布作为先验信息，采用多道同时反演的方法求取最大后验概率。

用组稀疏反演的方法得到地层品质因子，将该参数的倒数作为地震波能量衰减的定量表征，最终利用衰减识别油气。

本发明中的岩相驱动下的时频联合域地震反演方法，在利用反射系数向量稀疏性的同时，考虑了不同参数反射系数之间的结构稀疏特征，在贝叶斯框架下，实现组稀疏叠前四参数多道反演。不仅反演了弹性参数纵横波速度及密度，而且增加了与地层衰减相关的品质因子的反演结果。由于地层品质因子能够定量表述地震波的衰减，同时组稀疏方法也能得到横向连续性更好、稳定性更强的反演结果，因此该方法在油气识别，尤其是非常规油气藏、隐藏油气藏等中具有更大的优势。

以上具体实施方式仅用于说明本发明，而非用于限定本发明。

Claims (4)

1.岩相驱动下的时频联合域地震反演方法，其特征在于，该岩相驱动下的时频联合域地震反演方法包括：

步骤1，根据联合域卷积模型和部分角度叠加地震道集的时-频域响应，构建时频联合域地震正演模型；

步骤2，利用后验概率分布显式形式量化模型参数的最优解，推导出混合后验概率模型的显式解；

步骤3，根据混合后验概率密度分布的均值和协方差的显式表达式，计算获得三类已知数据集H约束下的地震反演最优解；

步骤4，利用储层测井、地质信息中模型参数中的低界和高界信息，进行基于边界约束的时频联合域地震反演，进而协同估算储层离散岩相；

在步骤2中，考虑到不同岩相下纵波阻抗的先验概率密度函数存在差异，将多个高斯分量的混合概率密度分布作为待反演模型的先验概率密度分布；假设目标工区中存在的储层岩性类型数量为M，针对线性地震反问题的贝叶斯推理，则模型参数m的最优解利用后验概率分布显式形式进行量化，推导得到三类已知数据集，S^f、S^t和η约束下的第k_th高斯分量O_k(m|H)的显式表达式

为，

另外，协方差矩阵

为：

式中，

式中，

和k分别为先验均值、先验协方差及第k_th个高斯分量的上标；S_ff、S_tt、S_ηη、S_ft、S_fη及S_tη表示已知模型参数m情况下正演观测数据，频率域地震数据S^f、时间域地震数据S^t和低频背景η的协方差矩阵；H为输入的三种类型已知数据集，时间域地震、频率域地震及低频先验数据，P为正演核矩阵，C_ff、C_tt、C_ηη、C_ft、C_fη及C_tη是3种已知数据集，频率域地震数据S^f、时间域地震数据S^t和纵波阻抗的低频背景η的先验协方差矩阵，C_H为三种数据集H的综合先验协方差矩阵，W^f为角度地震子波的对角线频谱矩阵，C为一阶差分矩阵，E代表Fourier变换算子或相移算子，D为低频正则化的对角线矩阵，W^t是时间域子波矩阵；

是矩阵S_ft的转置矩阵，

是矩阵S_fη的转置矩阵，

是矩阵S_tη的转置矩阵；

是矩阵C_ft的转置矩阵，

是矩阵C_fη的转置矩阵，

是矩阵C_tη的转置矩阵；λ_L为低频背景的正则化因子；

此外，待反演模型参数m的混合后验高斯概率密度分布O(m|H)重新表征为一个混合概率模型

其中不同岩相对应的高斯分量后验权重

等于

λ_k为第k_th个高斯分量的先验权重，

分别为已知数据集H的先验均值和先验协方差，M为储层岩性的分类数量，N_l(·)表示第l个高斯概率模型,N_k(·)表示第k个高斯概率模型；

在步骤3中，通过混合后验概率密度分布的均值和协方差的显式表达式，计算获得三类已知数据集H约束下的地震反演最优解；将实际测井数据和先前的模拟点作为地震反演的约束条件，并基于单点实现的序列模拟算法构建离散相序贯采样算法流程；

针对第k_th个高斯分量p_k(m|H)，在待模拟点m_i位置的后验均值和协方差改写为：

式中，

分别为待模拟点m_i位置的后验均值和先验均值，

为第k个高斯概率分量情况下4种类型，m_s、S^f、S^t及η条件数据集的协方差矩阵，m_s代表已知数据点和先前模拟点；

为第k个高斯概率分量的初始模型，U_i＝[0 L 0 1 L 0]为待模拟点m_i位置的采样矩阵即，第ith个元素为1，其余元素均为0，U为已知或先前模拟数据m_s的提取矩阵，S^f为频率域地震数据，S^t为时间域地震数据，η为低频背景，λ_L为低频背景的正则化系数，

为先前模拟数据m_s的正则化系数，W^f为角度地震子波的对角线频谱矩阵，E代表Fourier变换算子或相移算子，

为待反演模型参数m的先验协方差；当λ_f＝0、λ_t＝0和λ_L＝0的情况下，方程(8)等价于kriging插值和地质统计学中的序贯高斯模拟，

表示第四类已知数据集m_s的权重系数；此外，待模拟点m_i的协方差

改写为：

式中，

是第k个高斯分量中m_i的标准差，U_i＝[0 L 0 1 L 0]为待模拟点m_i位置的采样矩阵即，第ith个元素为1，其余元素均为0，已知或先前模拟数据m_s的提取矩阵为U＝[U_i；U_j；U_l；L；U_e]^t，

为四种类型，m_s、S^f、S^t及η条件数据集的协方差矩阵，m_s代表已知数据点和先前模拟点；

为第k个高斯概率分量的初始模型，U_i＝[0 L 0 1 L 0]为待模拟点m_i位置的采样矩阵，第ith个元素为1，其余元素均为0，U为已知或先前模拟数据m_s的提取矩阵，S^f为频率域地震数据，S^t为时间域地震数据，η为低频背景；λ_t和λ_f分别表示时间域地震数据与频率域地震数据的加权系数，λ_L为低频背景的正则化系数，

为先前模拟数据m_s的正则化系数，W^f为角度地震子波的对角线频谱矩阵，E代表Fourier变换算子或相移算子，

为待反演模型参数m的先验协方差；令F＝W^f·E·C，G＝W^t·C，

的显式表达式简化为，

根据基于贝叶斯分类，推导得到四类已知数据作为约束条件下第k_th个高斯分量的后验权重

更新为：

式中，

为序贯模拟过程中，四类条件数据的先验均值，

表示模型参数m在第k个高斯分量的先验均值，

为四种类型，m_s、S^f、S^t及η条件数据集的协方差矩阵，m_s代表已知数据点和先前模拟点；

为第k个高斯概率分量的初始模型，S^f为频率域数据，S^t为时间域数据，η为低频背景，N_l(·)表示第l个高斯概率模型，N_k(·)表示第k个高斯概率模型；

方程(11)中，待模拟点m_i的后验权重

用来判识离散岩相；当离散相态分类即优选混合概率模型中的高斯分量之后，纵波阻抗即连续型参数从被选择的高斯分量中随机采样获得；当访问下一个模拟点时，先前的模拟点将被视为已知数据，并且已知数据的邻域选择对于提高地震序贯模拟的计算效率具有重要作用，选用矩形邻域或圆形邻域，邻域半径通常设置为一个波长至两个波长之间。

2.根据权利要求1所述的岩相驱动下的时频联合域地震反演方法，其特征在于，在步骤1中，首先将联合域卷积模型视为正演算子，将部分角度叠加地震道集的时-频响应相结合，构建了时频联合域地震正演模型，将已知的纵波阻抗的低频先验信息引入正演方程中，以补偿观测地震数据中缺失的低频分量。

3.根据权利要求2所述的岩相驱动下的时频联合域地震反演方法，其特征在于，在步骤1中，将联合域卷积模型视为正演算子，

式中，S(f)、W(f)分别为入射角为θ_i时的地震频谱和子波频谱，f_v表示频率分量，r_pp(τ_z)为Russell线性近似AVO振幅随偏移距或入射角变化方程，f为频率，τ_z为时间序列，N(f)为高斯随机噪声的频谱向量；为了更清楚的阐述地震反问题，方程(1)重组为，

S^f＝W^f·E·C·m+N_f (2)

式中，W^f＝diag[W(f)]为角度地震子波的对角线频谱矩阵，E代表Fourier变换算子或相移算子，

为纵波阻抗的自然对数形式；Ip指的纵波阻抗，T表示转置运算，τ₁Lτ_n表示时间采样序列；n表示采样点数量，考虑到时间域AVO反演的强稳定性和频率域AVO反演的高分辨率，将部分角度叠加地震道集的时-频域响应相结合，构造了时频联合域正演算子，旨在稳定性和分辨率达到平衡，

式中，S^t、W^t和N_t分别是时间域部分角度叠加地震数据、时间域子波矩阵和高斯随机噪声，在不同信噪比SNRs的情况下，通过调整加权系数λ_t，λ_f、选择高信噪比的频率分量，进而控制地震反演的稳定性和分辨率；当SNRs高的情况下，提高频率域反演的权重λ_f，提高反演的分辨率；在SNRs较低的情况下，增加时域权重λ_t，选择信噪比较高的频率分量，进而提高预测的稳定性；将已知的纵波阻抗的低频先验信息引入正演方程(3)中，以补偿观测地震数据中缺失的低频分量，其罚函数J_low(m)即通过L2范数

引入目标函数中；η为通过实际测井数据的局部克里格插值构建的纵波阻抗的低频先验，D为正则化对角线矩阵；联合方程(3)和低频先验背景得到方程(4)，

式中，H为输入的三种类型已知数据集，时间域地震、频率域地震及低频先验数据，P为正演核矩阵，C_ff、C_tt、C_ηη、C_ft、C_fη及C_tη是3种已知数据集，频率域地震数据S^f、时间域地震数据S^t和纵波阻抗的低频背景η的先验协方差矩阵，C_H为三种数据集H的综合先验协方差矩阵，W^f为角度地震子波的对角线频谱矩阵，E代表Fourier变换算子或相移算子，

表示频率域地震数据与时间域地震数据的协方差矩阵C_ft的转置矩阵，

表示频率域地震数据与纵波阻抗低频背景的协方差矩阵C_fη的转置矩阵，

表示时间域地震数据与纵波阻抗低频背景的协方差矩阵C_tη的转置矩阵。

4.根据权利要求1所述的岩相驱动下的时频联合域地震反演方法，其特征在于，在步骤4中，由于方程(11)中并未考虑待反演弹性参数的下限和上限，由于实际地震数据受噪声干扰，反演结果往往出现一些超边界的解，非线性边界转化算法为

式中，x、m^min、m^max及N分别是待反演模型参数m的非线性中间变量、下界先验、上界先验和质控常数；令更新后的弹性参数Δm＝Δm^k＝m^k-m^k-1，则第kth次迭代后的相对弹性参数更新为：

式中，m^k-1表示第(k-1)_th次迭代后的估计结果，式(13)中，N用于控制弹性参数的梯度

N的范围为(0.1,2.0]，控制超过边界的预测结果m^min和m^max的更新幅度。

Images