CN112965103A - 一种多重孔隙储层叠前地震概率化多道反演方法 - Google Patents

Abstract

本发明公开了一种多重孔隙储层叠前地震概率化多道反演方法。该方法包括：步骤1，推导含多重孔隙空间的岩石弹性模量表达式；步骤2，推导多重孔隙储层物性参数表征的地震反射系数方程；步骤3，验证多重孔隙储层反射系数的精度和反演可行性；步骤4，构建待反演模型参数的后验概率密度分布及目标泛函；步骤5，研发多马尔科夫链随机采样的叠前地震多道逐级反演算法；步骤6，研发基于逐步模拟策略的岩石物理参数反演方法。本发明考虑了储层孔隙结构对地震反射系数的影响，研发了多重孔隙储层的地震反射系数参数化方法及叠前地震概率化多道反演技术，实现多重孔隙体积分数、流体体积模量及孔隙度等参数的稳定反演。

Description

一种多重孔隙储层叠前地震概率化多道反演方法

技术领域

本发明针对油气地震勘探领域，涉及复杂孔隙储层的岩石物理建模和叠前地震反演方法，特别是致密砂岩、碳酸盐岩及页岩等含油气储层的定量表征。

背景技术

以致密砂岩、碳酸盐岩及页岩等为代表的复杂油气储层是近年来油气勘探开发的重要研究领域，此类储层通常发育多种几何形态的孔隙空间，属于典型的多重孔隙储层。多重孔隙储层的物性参数和含油气性，与油气储量及储层产能密切相关。在地震定量解释中，有效预测及评价复杂孔隙储层的物性参数和含油气性是地球物理学者研究的热点话题。岩石物理理论搭建了储层物性参数与岩石弹性参数之间的桥梁。岩石物理指导下的地震AVO反演(amplitude variation with offset)是含油气储层定量表征的重要技术之一，被广泛应用于弹性参数反演、储层物性参数预测及流体识别等领域。在地震弹性参数反演中，估计的纵横波速度、纵横波阻抗、拉梅参数、杨氏模量和密度等弹性模量，将通过岩石物理经验关系被间接转化为储层物性参数，如孔隙度、流体饱和度、岩性等。然而，这种间接预测模式会增加储层定量表征的累计误差。在储层物性参数预测中，岩石物理反演和地震AVO物性参数反演是最主要的两种方法。岩石物理反演的输入是弹性参数数据体，它的正演算子是含油气储层的岩石物理模型或岩石物理解释模版。地震AVO物性参数反演联合岩石物理模型和地震反射系数方程，搭建了地震反射数据与储层物性参数之间的正演算子，可以实现从地震反射数据到储层物性参数的直接反演。

在地震流体识别中，敏感流体因子的构建和AVO正演模型的选择是最关键的两项内容。然而，地下含油气饱和岩石中通常发育多种几何形态的孔隙空间，并对岩石的弹性模量和地震响应影响较大。现有的岩石物理反演和地震流体识别方法主要被用来预测含油气储层的孔隙度和流体指示因子，却没有考虑孔隙结构对岩石弹性模量和地震AVO反射振幅的影响。然而，多重孔隙储层的物性参数与地震AVO反射系数之间存在非线性耦合关系，导致储层物性参数的地震AVO反演及流体识别方法表现出较大的计算复杂度与不适定性。

发明内容

为了解决多重孔隙储层岩石物理参数反演问题,本发明充分考虑岩石孔隙结构对岩石弹性模量和地震AVO反射振幅的影响，提供了一种多重孔隙储层叠前地震概率化多道反演方法及系统，提高了多重孔隙储层岩石物理参数反演稳定性和油气识别的精度。

为实现上述目的，解决现有方法存在的关键问题，本发明采用下述技术方案：一种多重孔隙储层叠前地震概率化多道反演方法，其包括如下步骤：

步骤1，推导含多重孔隙空间的岩石弹性模量表达式；

步骤2，推导多重孔隙储层物性参数表征的地震反射系数方程；

步骤3，验证地震反射系数方程的近似精度和反演可行性；

步骤4，构建待反演模型参数的后验概率密度分布及目标泛函；

步骤5，研发多马尔科夫链随机采样的叠前地震多道逐级反演算法；

步骤6，研发基于逐步模拟策略的岩石物理参数反演方法。

优选的，所述步骤1包括：

根据比奥特-伽斯曼孔隙弹性理论，饱和岩石的体积模量和剪切模量的表达式为：

K_s≈K_d+[1-K_d/K_m]²φ^-1K_f (1)，

μ_s＝μ_d (2)，

式(1)、式(2)中，K_s和μ_s为饱和岩石的体积模量和剪切模量，K_d和μ_d为干燥岩石骨架的体积模量和剪切模量，K_m为岩石基质的体积模量，φ为孔隙度，K_f为孔隙流体的体积模量；

将科斯-徐模型用于计算干燥岩石骨架的体积模量和剪切模量：

K_d/K_m≈(1-φ)^p (3),

μ_d/μ_m≈(1-φ)^q (4),

式(3)、式(4)中，p和q表示与孔隙纵横比相关的等效极化因子；

假设不同纵横比的孔隙空间划分为小纵横比的软孔隙和大纵横比的硬孔隙，等效极化因子p和q写成，

p＝α_sν_s+(1-α_s)ν_n (5),

式(5)、式(6)中，1-α_s和α_s分别是软孔隙体积分数和硬孔隙体积分数,软孔隙体积分数是软孔隙度与总孔隙度的比值，硬孔隙体积分数是硬孔隙度与总孔隙度的比值，系数ν_s和

描述了硬孔隙对干燥岩石骨架的影响，系数ν_n和

描述了软孔隙对干燥岩石骨架的影响；

令，G(φ)＝(1-φ)^p，L(φ)＝(1-φ)^q，将式(3)和式(4)代入式(1)、式(2)中，得到饱和岩石的体积模量和剪切模量的表达式，

K_s≈K_mG(φ)+[1-G(φ)]²φ^-1K_f (7)，

μ_s≈μ_mL(φ) (8)，

f≈[1-G(φ)]²φ^-1K_f (9)，

式(9)中，f为比奥特-伽斯曼理论中的伽斯曼流体项。

优选的，所述步骤2包括：

设线性化的平面纵波地震AVO反射系数方程R_pp(θ)为：

式(10)中，θ为平面纵波的入射角，

和

分别为干燥岩石骨架的纵横波速度比的平方和饱和岩石的纵横波速度比的平方，

分别为顶底地层伽斯曼流体项、饱和岩石剪切模量及密度的平均值，Δf、Δμ_s、Δρ_s分别为顶底地层伽斯曼流体项、饱和岩石剪切模量及密度的绝对变化量；

推导式(8)、式(9)的全微分计算形式，

式(11)、式(12)中，

结合式(8)、式(9)、式(11)、式(12)，求得伽斯曼流体项f的地层变化率

和剪切模量的地层变化率

式(13)、式(14)中，

为顶底地层流体体积模量、孔隙度、硬孔隙体积分数、岩石基质剪切模量的平均值，ΔK_f、Δφ、Δα_s、Δμ_m为顶底地层流体体积模量、孔隙度、硬孔隙体积分数、岩石基质剪切模量的绝对变化量；

将式(13)和式(14)代入式(10)中，推导得到体现硬孔隙体积分数α_s的地震AVO反射系数近似方程：

重新组合式(15)中的各项近似公式，将式(15)写成待反演岩石物理参数地层反射率的线性叠加形式：

式(16)中，a(θ)，b(θ)，c(θ)，d(θ)，e(θ)为与入射角θ和岩石物理参数m相关的系数，

优选的，所述步骤3包括：

设置岩石基质的体积模量和剪切模量、孔隙度、不同孔隙的体积分数、流体体积模量、硬孔隙纵横比、软孔隙纵横比；

依据精确Zoeppritz方程、Aki-Richards近似公式、Russell近似公式及式(16)，计算四类地震AVO反射系数及近似误差，确认式(16)满足了地震AVO反演的精度需求；

分析四类地震AVO反射系数的精度，阐明岩石基质体积模量和剪切模量、孔隙度、不同孔隙的体积分数、流体体积模量、硬孔隙纵横比、软孔隙纵横比对地震AVO反射系数的贡献度及反演可行性，确认孔隙度和硬孔隙体积分数的反演不确定性最小，流体体积模量和基质剪切模量的反演不确定性增大，密度的不确定性在所有待反演的岩石物理参数中是最大的。

优选的，所述步骤4包括：

结合式(16)和褶积模型，第x道地震数据S(θ)可以通过式(17)计算得到：

S(θ,x,t)＝W(θ)R_pp(θ,m,x,t)+N(θ,x,t) (17)，

式(17)中，m＝[φ α_s μ_m K_f ρ]^T为待反演岩石物理参数，W(θ)为入射角为θ的子波矩阵，R_pp(θ,m,x,t)为与待反演岩石物理参数、入射角θ相关的反射系数序列，N(θ,x,t)为随机噪声；

贝叶斯公式通过观测数据和待反演模型参数的先验信息，计算得到未知模型参数的后验概率密度分布p(m|S(θ,x,t))，

p(m|S(θ,x,t))∝p(S(θ,x,t)|m)·p(m) (18)，

式(18)中，p(S(θ,x,t)|m)为观测数据的似然概率密度分布，p(m)为未知模型参数的先验知识。令，同时反演的地震道数量为Z，每道N个采样点，待反演模型参数的数量为T，入射角度θ的数量为M；

设未知的模型参数m_i之间不具有相关性，m_i均服从参数为μ_mi和σ_mi的Laplace概率密度分布L(m_i；μ_mi,σ_mi)；设相邻道岩石物理参数之差服从均值为0、协方差为Σ_m的高斯概率密度函数N(Dm；0,Σ_m)；待反演模型参数的先验概率密度分布p(m)被写成，

式(19)中，D为相邻道模型参数在相同时间点位置的横向差分算子，Σ_m为相邻道模型参数之差Dm的协方差矩阵，μ_mi和σ_mi为Laplace概率密度分布中的常系数，且μ_mi＞0；

相邻2道地震数据，i道和i+1道，同时反演，待反演模型参数m和一阶差分平滑算子D被写成：

m＝[m_(i) m_(i+1)]^T,m_(i)＝[φ_(i) α_s(i) K_m(i) μ_m(i) K_f(i) ρ_s(i)]^T (20)，

D＝[E_TN -E_TN] (21)，

式(20)、式(21)中，E_TN为TN阶的单位对角线矩阵；

相邻3道地震数据，i-1道、i道及i+1道，同时反演，二阶差分平滑算子D被改写成：

m＝[m_(i) m_(i+1) m_(i+2)]^T, D_i＝[E_TN -2E_TN E_TN] (22)，

依次类推，n阶差分算子利用n-1阶差分算子的差分计算得到；

将观测噪声的高斯概率分布赋给似然函数p(S(θ,x,t)|m)，

p(S(θ,x,t)|m)＝N(S-WR_pp；0,Σ_S) (23)，

式(23)中，Σ_S为地震随机噪声N(θ,x,t)的协方差矩阵；

将式(19)和式(23)代入式(18)中，得到待反演岩石物理参数的后验概率密度分布p(m|S(θ,x,t))，

对式(24)两边取自然对数运算，得到与式(24)等价的线性表达式，

优选的，所述步骤5包括：

设当前模拟地震道集为第i个道集([i-1,i,i+1])，马尔科夫链条数量为P，总迭代次数为U，第t条马尔科夫链在第u+1次迭代时的更新状态为：

式(26)中，

和

分别为第t、t1和t2条马尔科夫链经过u次迭代后的最终状态，

是均值为

协方差为∑_p的高斯概率密度分布，F_t ^u表示与迭代次数相关的自适应尺度因子；

根据式(23)中的等价目标泛函O(m|S(θ,x,t))，Metropolis算法的接受概率的对数形式可以被改写为：

式(27)中，

为第t条马尔科夫链在第u+1次迭代时的接受概率。

优选的，所述步骤6包括：

将总迭代次数为U分为三阶段，第一阶段(U₁)、第二阶段(U₂)、第三阶段(U₃)。在u＜U₁的迭代过程中，固定岩石基质剪切模量μ_m、流体体积模量K_f、密度ρ不变，优化孔隙度φ、硬孔隙体积分数α_s，模拟后验概率分布、数值解；

在U₁＜u＜U₂的迭代过程中，固定孔隙度φ、硬孔隙体积分数α_s及密度ρ不变，优化流体体积模量K_f、岩石基质剪切模量μ_m，模拟后验概率分布、数值解；

在u＞U₃的迭代过程中，固定岩石基质剪切模量μ_m、流体体积模量K_f、孔隙度φ、硬孔隙体积分数α_s不变，优化岩石密度ρ，模拟待反演岩石物理参数的后验概率密度分布、数值解及不确定性。

本发明提出了一种孔隙度、硬孔隙体积分数及流体体积模量的叠前地震概率化多道反演方法。我们联合科斯-徐模型和Russell线性近似公式，推导了含硬孔隙体积分数和流体等效体积模量的地震AVO反射系数方程。待反演岩石物理参数的后验分布及目标泛函，是在模型参数纵向的稀疏约束及横向连续性约束的基础上建立的。孔隙度和硬孔隙体积分数的反演稳定性较高。岩石基质剪切模量、流体体积模量及密度的反演不确定性相对增大。逐级反演算法主要包含了3个反演阶段，在高信噪比的情况下，能够实现5个岩石物理参数的逐步预测，改善了流体体积模量和密度参数反演稳定性差的问题。多道反演策略引入了模型参数的横向约束信息，相比单道反演策略的反演横向连续性更好。实例应用表明，孔隙度、硬孔隙体积分数及孔隙流体模量的同时反演有利于碎屑岩储层物性和含流体性的定量解释。

附图说明

附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。为了更清楚地说明本发明的具体实施过程或技术方案，附图如下：

图1为本发明一种多重孔隙储层叠前地震概率化多道反演方法的流程示意图。

图2a依据精确Zoeppritz方程、Aki-Richards公式及式(16)计算得到的四类地震AVO反射系数示意图。

图2b多重孔隙储层反射系数方程的误差对比示意图。

图3a孔隙度对地震反射系数的影响分析示意图。

图3b硬孔隙体积分数对地震反射系数的影响分析示意图。

图4a流体体积模量对地震反射系数的影响分析示意图。

图4b基质剪切模量对地震反射系数的影响分析示意图。

图5密度对地震反射系数的影响分析示意图。

图6基于随机路径的地震概率化多道随机模拟示意图。

图7叠前地震概率化反演的逐级模拟策略示意图。

图8无噪声情况下，本发明反演预测的岩石物理参数结果。

图9利用常规叠前地震反演方法估计的岩石物理参数结果。

图10信噪比等于4的情况下，本发明反演预测的岩石物理参数结果。

图11根据实际测井数据解释得到的岩石物理参数示意图。

图12井旁道位置，本发明反演预测的岩石物理参数示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明充分考虑岩石孔隙结构对岩石弹性模量和地震AVO反射振幅的影响，提供了一种多重孔隙储层叠前地震概率化多道反演方法(详见图1流程示意图)，提高了多重孔隙储层反演稳定性和油气识别精度。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

本发明一种多重孔隙储层叠前地震概率化多道反演方法，包括如下步骤：

步骤1，推导含多重孔隙空间的岩石弹性模量表达式；

步骤2，推导多重孔隙储层物性参数表征的地震反射系数方程；

步骤3，验证地震反射系数方程的近似精度和反演可行性；

步骤4，构建待反演模型参数的后验概率密度分布及目标泛函；

步骤5，研发多马尔科夫链随机采样的叠前地震多道逐级反演算法；

步骤6，研发基于逐步模拟策略的岩石物理参数反演方法。

步骤1中，联合科斯-徐(Keys-Xu)干燥岩石骨架模型和比奥特-伽斯曼(Biot-Gassmann)孔隙弹性理论，将岩石孔隙空间等效为具有两类孔隙纵横比的椭球形孔隙，推导含多重孔隙空间的岩石弹性模量表达式，分析多重孔隙体积分数、孔隙度及流体体积模量对饱和岩石弹性模量的影响。

在均匀各向同性假设、孔隙连通不排液情况下，比奥特-伽斯曼孔隙弹性理论描述了饱和岩石弹性模量(K_s和μ_s)与岩石基质弹性模量(K_m和μ_m)、干燥岩石骨架弹性模量(K_d和μ_d)、孔隙度φ及孔隙流体模量K_f之间的关系。由于岩石基质模量K_m远大于孔隙流体模量K_f，针对中高孔碎屑岩储层(φ＞0.10)，

故可忽略掉弱小项。因此，突出流体体积模量和孔隙度的饱和岩石弹性模量(体积模量和剪切模量)的表达式可以写成，

K_s≈K_d+[1-K_d/K_m]²φ^-1K_f (1)，

μ_s＝μ_d (2)。

式(1)、式(2)中，K_s和μ_s为饱和岩石的体积模量和剪切模量，K_d和μ_d为干燥岩石骨架的体积模量和剪切模量，K_m为岩石基质的体积模量，φ为孔隙度，K_f为孔隙流体的体积模量。为研究多重孔隙空间对岩石弹性模量和地震AVO反射系数的影响，将科斯-徐近似模型用于计算干燥岩石骨架的弹性模量，

K_d/K_m≈(1-φ)^p (3)，

μ_d/μ_m≈(1-φ)^q (4),

式(3)、式(4)中，p和q表示与孔隙纵横比相关的等效极化因子。

针对多重孔隙砂岩储层，假设不同纵横比的孔隙空间，可以划分为小纵横比的软孔隙和大纵横比的硬孔隙。基于此，等效极化因子p和q可以写成，

p＝α_sν_s+(1-α_s)ν_n (5)，

式(5)、式(6)中，1-α_s和α_s分别是软孔隙体积分数和硬孔隙体积分数。这些分数分别是软孔隙度和硬孔隙度与总孔隙度的比值。系数ν_s和

描述硬孔隙对干燥岩石骨架的影响。系数ν_n和

描述了软孔隙对干燥岩石骨架的影响。

令，G(φ)＝(1-φ)^p，L(φ)＝(1-φ)^q，将式(3)和式(4)代入式(1)、(2)中，得到饱和岩石体积模量和剪切模量的表达式，

K_s≈K_mG(φ)+[1-G(φ)]²φ^-1K_f (7)，

μ_s≈μ_mL(φ) (8)，

f≈[1-G(φ)]²φ^-1K_f (9)，

式(9)中，f为比奥特-伽斯曼理论中的伽斯曼(Gassmann)流体项。

步骤2中，基于比奥特-伽斯曼孔隙弹性理论，推导含孔隙度、硬孔隙体积分数、岩石基质剪切模量和流体体积模量的地震反射系数方程。

以多重孔隙储层为研究对象，本发明基于科斯-徐模型和伽斯曼近似方程，推导了考虑多重孔隙空间的地震AVO反射系数方程。设线性化的平面纵波地震AVO(amplitudevariation with offset)反射系数方程R_pp(θ)为，

式(10)中，θ为平面纵波的入射角，

和

分别为干燥岩石骨架和饱和岩石的纵横波速度比的平方，

分别为顶底地层伽斯曼流体项、饱和岩石剪切模量及密度的平均值，Δf，Δμ_s，Δρ_s分别为顶底地层伽斯曼流体项、饱和岩石剪切模量及密度的绝对变化量。

为降低弹性模量的非线性程度，推导式(8)-(9)的全微分计算形式，

式中，

结合式(8)、式(9)、式(11)、式(12)，进一步求得伽斯曼流体项f的地层变化率

和剪切模量的地层变化率

式(13)、式(14)中，

为顶底地层岩石物理参数(流体体积模量、孔隙度、硬孔隙体积分数、岩石基质剪切模量)的平均值。ΔK_f，Δφ，Δα_s和Δμ_m为顶底地层岩石物理参数的绝对变化量。

将式(13)和式(14)代入公式(10)中，推导得到体现硬孔隙体积分数α_s的地震反射系数方程，

重新组合式(15)中的各项近似公式，将式(15)写成待反演岩石物理参数地层反射率的线性叠加形式，如下所示，

式(16)中，a(θ)，b(θ)，c(θ)，d(θ)，e(θ)是与入射角θ和岩石物理参数m相关的系数，

步骤3中，分析多重孔隙储层四类地震AVO(amplitude variation with offset)反射系数，验证地震反射系数方程的近似精度和反演可行性，请参阅图2a、图2b、图3a、图3b、图4a、图4b、图5。

本发明在开展叠前地震反演之前，必须验证多重孔隙储层反射系数的精度。为了阐明式(16)的近似精度，开展四套含气砂岩储层的地震反射特征分析。

预设四套含气砂岩的岩石物理参数，即，设置岩石基质的体积模量和剪切模量、孔隙度、不同孔隙的体积分数、流体体积模量、硬孔隙纵横比、软孔隙纵横比；

依据精确Zoeppritz方程、Aki-Richards近似公式、Russell近似公式及式(16)，计算四类地震AVO反射系数及近似误差，确认式(16)满足了地震AVO反演的精度需求；

分析四类地震AVO反射系数的精度，阐明岩石基质体积模量和剪切模量、孔隙度、不同孔隙的体积分数、流体体积模量、硬孔隙纵横比、软孔隙纵横比对地震AVO反射系数的贡献度及反演可行性，确认孔隙度和硬孔隙体积分数的反演不确定性最小，流体体积模量和基质剪切模量的反演不确定性增大，密度的不确定性在所有待反演的岩石物理参数中是最大的。

本发明一实施例针对目标工区的实际测井数据，预设四套含气砂岩的岩石物理参数开展了四套含气砂岩储层的地震反射特征分析。

图2a与图2b分别为依据精确Zoeppritz方程、Aki-Richards近似公式、Russell近似公式及式(16)计算得到的四类地震AVO反射系数R_pp(θ)和各类反射系数的近似误差。图中可见，式(16)与精确Zoeppritz方程、Aki-Richards近似及Russell近似的计算结果保持了较高的一致性。在入射角θ≤30°的情况下，近似误差得到了较好的控制。含气砂岩的四类振幅随入射角变化特征能够被有效的体现出来。因此，式(16)满足了地震AVO反演的精度需求。

图3a、图3b、图4a、图4b、图5中，本发明测试了式(16)中岩石物理参数对地震AVO反射系数的贡献度。本发明以第三类含气砂岩为例，图3a、图3b、图4a、图4b和图5分别展示了孔隙度、硬孔隙体积分数、流体体积模量、基质剪切模量及岩石密度对地震AVO反射系数的影响结果。图中可见，孔隙度和硬孔隙体积分数对地震反射系数的影响最大。即，在总孔隙度不变的情况下，硬孔隙体积分数不同，地震AVO反射系数同样会发生剧烈变化。此外，流体体积模量和密度对地震AVO反射系数的影响相对较小，基本被淹没在硬孔隙体积分数的影响中。因此，从AVO反射系数理论分析的角度，孔隙度和硬孔隙体积分数的反演不确定性最小，流体体积模量和基质剪切模量的反演不确定性增大，密度的不确定性在所有待反演的岩石物理参数中是最大的。

步骤4中，利用岩石物理参数具有垂直相关性的先验Laplace分布和横向相关性的先验高斯分布，构建待反演模型参数的后验概率密度分布及目标泛函。

针对式(16)所示非线性正演算子R_pp(θ,m)的地震反演和优化问题，本发明在贝叶斯推断框架下发展了蒙特卡洛-马尔科夫链概率化叠前地震反演算法。结合式(16)和褶积模型，第x道地震数据S(θ)可以通过下式计算得到：

S(θ,x,t)＝W(θ)R_pp(θ,m,x,t)+N(θ,x,t) (17)

式(17)中，m＝[φ α_s μ_m K_fρ]^T为待反演岩石物理参数，W(θ)为入射角为θ的子波矩阵，R_pp(θ,m,x,t)为与待反演岩石物理参数、入射角θ相关的反射系数序列，N(θ,x,t)为随机噪声。贝叶斯公式可以通过观测数据和待反演模型参数的先验信息，计算得到未知模型参数的后验概率密度分布p(m|S(θ,x,t))，

p(m|S(θ,x,t))∝p(S(θ,x,t)|m)·p(m) (18)

式(18)中，p(S(θ,x,t)|m)为观测数据的似然概率密度分布，p(m)为未知模型参数的先验知识。令，同时反演的地震道数量为Z，每道N个采样点，待反演模型参数的数量为T，入射角度θ的数量为M。本发明假设未知的模型参数m_i之间不具有相关性，m_i均服从参数为μ_mi和σ_mi的Laplace概率密度分布L(m_i；μ_mi,σ_mi)。为了提高随机模拟结果的横向连续性，假设相邻道岩石物理参数之差服从均值为0、协方差为Σ_m的高斯概率密度函数N(Dm；0,Σ_m)。因此，待反演模型参数的先验概率密度分布p(m)可以被写成，

式(19)中，D为相邻道模型参数在相同时间点位置的横向差分算子，Σ_m为相邻道模型参数之差Dm的协方差矩阵，μ_mi和σ_mi为Laplace概率密度分布中的常系数，且μ_mi＞0。以相邻2道地震数据(i道和i+1道)同时反演为例，待反演模型参数m和一阶差分平滑算子D可以被写成：

m＝[m_(i) m_(i+1)]^T,m_(i)＝[φ_(i) α_s(i) K_m(i) μ_m(i) K_f(i) ρ_s(i)]^T (20)

D＝[E_TN -E_TN] (21)

式(20)、式(21)中，E_TN为TN阶的单位对角线矩阵。以相邻3道地震数据(i-1道、i道及i+1道)同时反演为例，二阶差分平滑算子D可以被改写成：

m＝[m_(i) m_(i+1) m_(i+2)]^T, D_i＝[E_TN -2E_TN E_TN] (22)

依次类推，n阶差分算子可以利用n-1阶差分算子的差分计算得到。将观测噪声的高斯概率分布赋给似然函数p(S(θ,x,t)|m)，

p(S(θ,x,t)|m)＝N(S-WR_pp；0,Σ_S) (23)

式(23)中，Σ_S为地震随机噪声N(θ,x,t)的协方差矩阵。将式(19)和式(23)代入式(18)中，得到待反演岩石物理参数的后验概率密度分布p(m|S(θ,x,t))，

对式(24)两边取自然对数运算，得到与式(24)等价的线性表达式，

步骤5中，研发了多马尔科夫链随机采样的叠前地震多道逐级反演算法,请参阅图6(即，基于随机路径的叠前地震多道随机模拟的示意图)。

图6中，本发明以相邻3道地震数据同时反演为例，将3道地震数据编为一组，基于随机路径发展概率化叠前地震多道反演算法。此外，随机游走的地震道集将存在交叉的地震道，本方法将先前模拟的地震道作为当前模拟道集的约束条件。本发明将蒙特卡洛马尔科夫链方法应用到模型空间的随机抽样，进而获取模型参数的后验统计特征。该算法在提议分布

产生新状态

的过程中，融合了差分进化算法的自适应进化性质。假设当前模拟地震道集为第i个道集([i-1,i,i+1])，马尔科夫链条数量为P，总迭代次数为U。第t条马尔科夫链在第u+1次迭代时的更新状态为：

式(26)中，

和

分别为第t、t1和t2条马尔科夫链经过u次迭代后的最终状态，

是均值为

协方差为∑_p的高斯概率密度分布，F_t ^u表示与迭代次数相关的自适应尺度因子。根据式(23)中的等价目标泛函O(m|S(θ,x,t))，Metropolis算法的接受概率的对数形式可以被改写为：

式(27)中，

为第t条马尔科夫链在第u+1次迭代时的接受概率。

基于多马尔科夫链随机采样算法的叠前地震概率化反演伪代码如下所示。

步骤6中，研发基于逐步模拟策略的岩石物理参数反演方法，逐级预测岩石物理参数的后验概率密度分布、数值解及不确定性,请参阅图7(即，叠前地震概率化反演的逐级模拟策略示意图)。

步骤6包括三个阶段。孔隙度和硬孔隙体积分数是第一阶段需要反演的岩石物理参数，第二阶段估算了流体体积模量和岩石基质剪切模量，由于岩石密度在所有的岩石物理参数中，对地震AVO反射系数的影响最小，它将作为第三阶段待反演的模型参数。

步骤6的具体步骤为：

将总迭代次数为U分为三阶段，第一阶段(U₁)、第二阶段(U₂)、第三阶段(U₃)。在u＜U₁的迭代过程中，固定岩石基质剪切模量μ_m、流体体积模量K_f、密度ρ不变，优化孔隙度φ、硬孔隙体积分数α_s，模拟后验概率分布、数值解；

在U₁＜u＜U₂的迭代过程中，固定孔隙度φ、硬孔隙体积分数α_s及密度ρ不变，优化流体体积模量K_f、岩石基质剪切模量μ_m，模拟后验概率分布、数值解；

在u＞U₃的迭代过程中，固定岩石基质剪切模量μ_m、流体体积模量K_f、孔隙度φ、硬孔隙体积分数α_s不变，优化岩石密度ρ，模拟待反演岩石物理参数的后验概率密度分布、数值解及不确定性。

该策略，在高信噪比的情况下，能够实现5个岩石物理参数的逐步预测，避免了多参数同时模拟不稳定的问题，改善了流体体积模量和密度参数反演稳定性差的问题，提高了模型参数的反演精度。

本发明一种多重孔隙储层叠前地震概率化多道反演方法与国际上现有技术相比，具有下列优点：1)本发明充分考虑岩石孔隙结构、多重孔隙体积分数对弹性模量和地震反射系数的影响，常规方法欠缺考虑岩石多重孔隙空间对地震反射系数的影响；2)本发明中设计的多道反演策略，引入了模型参数的横向约束信息，相比单道反演策略的反演横向连续性更好；3)多马尔科夫链逐级模拟算法避免了多参数同时模拟不稳定的问题，提高了反演稳定性和预测精度。

为进一步说明本发明的可行性和有效性，下面举两个实施例：

实施例1：理论模型测试，详见图8、图9、图10。

为了更清晰地阐述本文AVO反演算法的抗噪性和可靠性，本专利利用井旁道重采样后的测井数据开展了理论模型测试。图8中白色虚线表示根据测井数据计算的岩石物理参数。图8和图9分别为无噪音和信噪比等于4情况下，估计得到的岩石物理参数，其中，白色散点为60次随机模拟的后验均值。

从图8中可以看出，孔隙度和硬孔隙体积分数反演结果的稳定性是最高的，与真实模型具有较强的相似性，从后验概率密度分布中产生的随机样本点的离散程度较低。此外，相比孔隙度和硬孔隙体积分数的反演结果，岩石基质剪切模量和流体体积模量的反演不确定性增大，其随机模拟结果的离散程度增大。由于岩石密度对地震AVO反射系数的贡献度最小，密度参数的反演不确定性在待反演岩石物理参数中是最大的。因此，在实际资料解释过程中，密度参数反演结果的可信程度是最低的。这种现象与理论岩石物理分析的结果是一致的。在信噪比等于4情况下(图9)，所有岩石物理参数均受到随机噪音的影响，反演不确定性也同样有所增大。然而，本方法的反演结果与真实的岩石物理参数仍然保持了较高的吻合度，进一步验证了本文反演算法的稳定性。

在不考虑碎屑岩多重孔隙效应的情况下，本专利开展了常规地震反演方法的测试。图10中黑色实线为真实的岩石物理参数，虚线展示了常规地震反演方法估计的岩石物理参数结果，三角形散点为本发明估计得到的岩石物理参数。图中可见，孔隙度的常规反演结果仍然保持了较高的精度。流体体积模量的预测精度却有较大幅度的下降。即使在地震数据匹配较好的情况下，岩石基质剪切模量和密度参数已经难以被有效地恢复出来。我们得出结论，孔隙结构或硬孔隙体积分数对地震反演的影响是不可以直接被忽略的，否则，其它岩石物理参数在反演过程中将会出现较大的误差。

实施例2：实际资料测试，详见图11、图12。

图11展示了根据实际测井数据解释得到的岩石物理参数，矩形虚线框位置表示四个含油砂岩储层的位置。图12为本发明反演预测的岩石物理参数，灰线表示20条马尔可夫链的随机模拟结果，白色散点表示本发明估计的后验均值解，虚线表示反演结果的95％置信区间，白色实线表示使用常规叠前地震反演方法估计的流体体积模量。

在验证方法可行性和算法稳定性的基础上，本发明将该方法应用于Z油田的地震勘探实例中，以验证方法对陆上资料的实用性。该工区发育着储层物性较好的高孔高渗碎屑岩油藏，埋藏深度约1160m～1320m。图11展示了根据实际测井数据解释得到的储层岩石物理参数。该井目的层段主要发育了4套含油砂岩储层和1套含油水层，如图中的矩形框所示。这口井将作为盲井来验证本发明的可靠性。

经过20条马尔科夫链的随机模拟，图12展示了井旁道的20条随机反演得到的孔隙度、硬孔隙体积分数及流体体积模量结果。白色散点表示本文方法反演得到的岩石物理参数，虚线为20次随机采样的95％置信区间，白色实线为不考虑多重孔隙效应的反演结果。图中可见，反演结果与实际测井解释结果在储层发育位置保持了较高的一致性。孔隙度和硬孔隙体积分数在砂岩储层位置发育明显高值，能够较好地指示高孔隙度的砂岩储层。流体体积模量反演结果在含油砂岩储层位置呈现明显的低值异常，较好地区分了油层及顶部水层。

通过实例1和实例2的研究表明，多重孔隙的干燥岩石骨架效应对流体识别结果干扰较大，在地震含油气性的定量解释中是不可以被忽略的，进一步验证了本发明方法的有效性和实际应用前景。

以上描述和实施例仅用于说明本发明，凡是在本发明技术方案的基础上进行的等同变化和改进，均应包含在本发明的保护范围之内。

Claims (7)

1.一种多重孔隙储层叠前地震概率化多道反演方法，其特征在于，其包括如下步骤：

步骤1，推导含多重孔隙空间的岩石弹性模量表达式；

步骤2，推导多重孔隙储层物性参数表征的地震反射系数方程；

步骤3，验证地震反射系数方程的近似精度和反演可行性；

步骤4，构建待反演模型参数的后验概率密度分布及目标泛函；

步骤5，研发多马尔科夫链随机采样的叠前地震多道逐级反演算法；

步骤6，研发基于逐步模拟策略的岩石物理参数反演方法。

2.如权利要求1所述多重孔隙储层叠前地震概率化多道反演方法，其特征在于，

所述步骤1包括：

根据比奥特-伽斯曼孔隙弹性理论，饱和岩石的体积模量和剪切模量的表达式为：

K_s≈K_d+[1-K_d/K_m]²φ^-1K_f (1)，

μ_s＝μ_d (2)，

式(1)、式(2)中，K_s和μ_s为饱和岩石的体积模量和剪切模量，K_d和μ_d为干燥岩石骨架的体积模量和剪切模量，K_m为岩石基质的体积模量，φ为孔隙度，K_f为孔隙流体的体积模量；

将科斯-徐模型用于计算干燥岩石骨架的体积模量和剪切模量：

K_d/K_m≈(1-φ)^p (3),

μ_d/μ_m≈(1-φ)^q (4),

式(3)、式(4)中，p和q表示与孔隙纵横比相关的等效极化因子；

假设不同纵横比的孔隙空间划分为小纵横比的软孔隙和大纵横比的硬孔隙，等效极化因子p和q写成，

p＝α_sν_s+(1-α_s)ν_n (5),

式(5)、式(6)中，1-α_s和α_s分别是软孔隙体积分数和硬孔隙体积分数,软孔隙体积分数是软孔隙度与总孔隙度的比值，硬孔隙体积分数是硬孔隙度与总孔隙度的比值，系数ν_s和

描述了硬孔隙对干燥岩石骨架的影响，系数ν_n和

描述了软孔隙对干燥岩石骨架的影响；

令，G(φ)＝(1-φ)^p，L(φ)＝(1-φ)^q，将式(3)和式(4)代入式(1)、式(2)中，得到饱和岩石的体积模量和剪切模量的表达式，

K_s≈K_mG(φ)+[1-G(φ)]²φ^-1K_f (7)，

μ_s≈μ_mL(φ) (8)，

f≈[1-G(φ)]²φ^-1K_f (9)，

式(9)中，f为比奥特-伽斯曼理论中的伽斯曼流体项。

3.如权利要求2所述多重孔隙储层叠前地震概率化多道反演方法，其特征在于，

所述步骤2包括：

设线性化的平面纵波地震AVO反射系数方程R_pp(θ)为：

式(10)中，θ为平面纵波的入射角，

和

分别为干燥岩石骨架的纵横波速度比的平方和饱和岩石的纵横波速度比的平方，

分别为顶底地层伽斯曼流体项、饱和岩石剪切模量及密度的平均值，Δf、Δμ_s、Δρ_s分别为顶底地层伽斯曼流体项、饱和岩石剪切模量及密度的绝对变化量；

推导式(8)、式(9)的全微分计算形式，

式(11)、式(12)中，

结合式(8)、式(9)、式(11)、式(12)，求得伽斯曼流体项f的地层变化率

和剪切模量的地层变化率

式(13)、式(14)中，

为顶底地层流体体积模量、孔隙度、硬孔隙体积分数、岩石基质剪切模量的平均值，ΔK_f、Δφ、Δα_s、Δμ_m为顶底地层流体体积模量、孔隙度、硬孔隙体积分数、岩石基质剪切模量的绝对变化量；

将式(13)和式(14)代入式(10)中，推导得到体现硬孔隙体积分数α_s的地震AVO反射系数近似方程：

重新组合式(15)中的各项近似公式，将式(15)写成待反演岩石物理参数地层反射率的线性叠加形式：

式(16)中，a(θ)，b(θ)，c(θ)，d(θ)，e(θ)为与入射角θ和岩石物理参数m相关的系数，

4.如权利要求3所述多重孔隙储层叠前地震概率化多道反演方法，其特征在于，

所述步骤3包括：

设置岩石基质的体积模量和剪切模量、孔隙度、不同孔隙的体积分数、流体体积模量、硬孔隙纵横比、软孔隙纵横比；

依据精确Zoeppritz方程、Aki-Richards近似公式、Russell近似公式及式(16)，计算四类地震AVO反射系数及近似误差，确认式(16)满足了地震AVO反演的精度需求；

分析四类地震AVO反射系数的精度，阐明岩石基质体积模量和剪切模量、孔隙度、不同孔隙的体积分数、流体体积模量、硬孔隙纵横比、软孔隙纵横比对地震AVO反射系数的贡献度及反演可行性，确认孔隙度和硬孔隙体积分数的反演不确定性最小，流体体积模量和基质剪切模量的反演不确定性增大，密度的不确定性在所有待反演的岩石物理参数中是最大的。

5.如权利要求4所述多重孔隙储层叠前地震概率化多道反演方法，其特征在于，

所述步骤4包括：

结合式(16)和褶积模型，第x道地震数据S(θ)可以通过式(17)计算得到：

S(θ,x,t)＝W(θ)R_pp(θ,m,x,t)+N(θ,x,t) (17)，

式(17)中，m＝[φ α_s μ_m K_f ρ]^T为待反演岩石物理参数，W(θ)为入射角为θ的子波矩阵，R_pp(θ,m,x,t)为与待反演岩石物理参数、入射角θ相关的反射系数序列，N(θ,x,t)为随机噪声；

贝叶斯公式通过观测数据和待反演模型参数的先验信息，计算得到未知模型参数的后验概率密度分布p(m|S(θ,x,t))，

p(m|S(θ,x,t))∝p(S(θ,x,t)|m)·p(m) (18)，

式(18)中，p(S(θ,x,t)|m)为观测数据的似然概率密度分布，p(m)为未知模型参数的先验知识。令，同时反演的地震道数量为Z，每道N个采样点，待反演模型参数的数量为T，入射角度θ的数量为M；

设未知的模型参数m_i之间不具有相关性，m_i均服从参数为μ_mi和σ_mi的Laplace概率密度分布L(m_i；μ_mi,σ_mi)；设相邻道岩石物理参数之差服从均值为0、协方差为Σ_m的高斯概率密度函数N(Dm；0,Σ_m)；待反演模型参数的先验概率密度分布p(m)被写成，

式(19)中，D为相邻道模型参数在相同时间点位置的横向差分算子，Σ_m为相邻道模型参数之差Dm的协方差矩阵，μ_mi和σ_mi为Laplace概率密度分布中的常系数，且μ_mi＞0；

相邻2道地震数据，i道和i+1道，同时反演，待反演模型参数m和一阶差分平滑算子D被写成：

m＝[m_(i) m_(i+1)]^T,m_(i)＝[φ_(i) α_s(i) K_m(i) μ_m(i) K_f(i) ρ_s(i)]^T (20)，

D＝[E_TN -E_TN] (21)，

式(20)、式(21)中，E_TN为TN阶的单位对角线矩阵；

相邻3道地震数据，i-1道、i道及i+1道，同时反演，二阶差分平滑算子D被改写成：

m＝[m_(i) m_(i+1) m_(i+2)]^T, D_i＝[E_TN -2E_TN E_TN] (22)，

依次类推，n阶差分算子利用n-1阶差分算子的差分计算得到；

将观测噪声的高斯概率分布赋给似然函数p(S(θ,x,t)|m)，

p(S(θ,x,t)|m)＝N(S-WR_pp；0,Σ_S) (23)，

式(23)中，Σ_S为地震随机噪声N(θ,x,t)的协方差矩阵；

将式(19)和式(23)代入式(18)中，得到待反演岩石物理参数的后验概率密度分布p(m|S(θ,x,t))，

对式(24)两边取自然对数运算，得到与式(24)等价的线性表达式，

6.如权利要求5所述多重孔隙储层叠前地震概率化多道反演方法，其特征在于，

所述步骤5包括：

设当前模拟地震道集为第i个道集([i-1,i,i+1])，马尔科夫链条数量为P，总迭代次数为U，第t条马尔科夫链在第u+1次迭代时的更新状态为：

式(26)中，

和

分别为第t、t1和t2条马尔科夫链经过u次迭代后的最终状态，

是均值为

协方差为∑_p的高斯概率密度分布，F_t ^u表示与迭代次数相关的自适应尺度因子；

根据式(23)中的等价目标泛函O(m|S(θ,x,t))，Metropolis算法的接受概率的对数形式可以被改写为：

式(27)中，

为第t条马尔科夫链在第u+1次迭代时的接受概率。

7.如权利要求6所述多重孔隙储层叠前地震概率化多道反演方法，其特征在于，

所述步骤6包括：

将总迭代次数为U分为三阶段，第一阶段(U₁)、第二阶段(U₂)、第三阶段(U₃)。在u＜U₁的迭代过程中，固定岩石基质剪切模量μ_m、流体体积模量K_f、密度ρ不变，优化孔隙度φ、硬孔隙体积分数α_s，模拟后验概率分布、数值解；

在U₁＜u＜U₂的迭代过程中，固定孔隙度φ、硬孔隙体积分数α_s及密度ρ不变，优化流体体积模量K_f、岩石基质剪切模量μ_m，模拟后验概率分布、数值解；

在u＞U₃的迭代过程中，固定岩石基质剪切模量μ_m、流体体积模量K_f、孔隙度φ、硬孔隙体积分数α_s不变，优化岩石密度ρ，模拟待反演岩石物理参数的后验概率密度分布、数值解及不确定性。

Images