0 引言

连通的裂缝不仅能够为油气的储存提供孔隙空间，而且增加了岩石的渗透率，因此识别裂缝型储层的位置，估测裂缝特征参数，实现裂缝特征的有效描述，对裂缝储层的勘探和开发具有尤为重要的作用(Liu and Martinez, 2012；陈怀震等，2014a).在长波长假设条件下，各向同性背景地层中发育一组平行排列的垂直裂缝可等效为具有水平对称轴的横向各向同性(HTI)介质(Rüger, 1997, 1998；Thomsen，2002).

随着宽方位地震采集与处理技术的成熟，利用宽方位地震信息获取裂缝储层各向异性参数信息受到了广大地球物理学者的青睐.基于方位地震数据振幅随入射角和方位角的变化(AVAZ)的各向异性反演方法能够充分利用方位地震信息，获取裂缝储层的各向异性参数，实现裂缝储层裂缝发育特征的描述(Mallick et al., 1998；Shaw and Sen, 2006；Bachrach et al., 2009；Downton and Roure, 2010；张广智等，2012；陈怀震等，2014b).然而，纵波AVAZ反射系数方程中各向同性背景反射系数项与各向异性扰动系数项存在明显的表达差异，而参数表达差异则严重影响弹性参数与各向异性参数反演的精度.为了消除各向同性背景反射系数项与各向异性扰动项的参数表达差异，我们引入了表征HTI介质各向异性强弱特征的拟裂缝弱度参数，并推导了相应的纵波AVAZ反射系数近似公式，基于AVAZ反演实现储层裂缝岩石弹性参数与各向异性参数的有效估测.

储层物性参数是进行储层预测与评价、估算含油气储量、确定钻井井位的重要依据.目前为止，基于神经网络的反演方法(印兴耀等, 1994, 1998)和基于地质统计学的随机模拟方法(González et al., 2008；Bosch et al., 2010)最为常用，但这两种方法过度依赖于测井资料.基于贝叶斯理论的概率反演是一套结合统计岩石物理模型、贝叶斯反演理论和蒙特卡罗抽样方法的储层物性参数和岩性/流体类型的地震反演技术，能够有效地实现储层物性特征的描述(Mukerji et al., 2001；Avseth et al., 2005；Spiles，2008；Bachrach et al., 2009；Grana and Rossa, 2010；Rimstad et al，2012；胡华锋等，2012；田玉坤等，2013；印兴耀等，2014；袁成等，2016).宽方位地震资料不仅具有丰富的各向异性信息，同时也蕴含着丰富的孔隙度等储层物性信息，其作为衡量储层质量及孔隙承载含流体状况的重要参数，为储层预测与孔隙流体识别提供了强有力的参考依据(Spikes，2008；Tsvankin and Grechka, 2011).因此，本文结合实际地震资料提出了贝叶斯框架下岩石物理驱动的各向异性介质多种储层参数概率地震联合反演方法：首先基于AVAZ反演裂缝岩石的弹性参数与各向异性参数，并在此基础上通过统计岩石物理模型表征孔隙度、裂缝密度等各向异性介质储层参数与裂缝岩石弹性参数及各向异性参数的相互关联，并采用马尔科夫链蒙特卡洛(Markov chain Monte Carlo，MCMC)抽样方法进行大量样本的随机模拟，使用期望最大化(Expectation Maximization，EM)算法估计后验条件概率分布，最终寻找最大后验条件概率对应的孔隙度、裂缝密度等各向异性介质储层参数即为反演结果.

1 方法原理 1.1 裂缝岩石弹性参数与各向异性参数AVAZ反演

裂缝弱度参数是储层裂缝描述的特征参数，有助于指导地下的裂缝识别.Hsu和Schoenberg(1993)引入了一组无量纲的参数，用来描述由于裂缝引起的岩石刚度变化特征：

(1)

式中，M₀和μ₀是背景介质的纵波模量与横波模量；正定裂缝柔度参数Z_N与Z_T分别表征法向裂缝柔度和切向裂缝柔度；Δ_N与Δ_T分别表征法向裂缝弱度和切向裂缝弱度，代表垂直于裂缝面方向和平行于裂缝面方向的由裂缝引起的弹性参数差值.

裂缝弱度参数的有效估算有助于各向异性介质的特征描述，而弹性逆反射理论是介质参数反演的有效途径(印兴耀等，2013).基于弹性逆反射理论和稳相法(Shaw and Sen, 2006)，推导了各向异性HTI介质的纵波AVAZ反射系数方程(详见附录A)：

(2)

式中，

(3)

θ是入射角；ϕ是方位角；R_pp^iso(θ)与R_pp^ani(θ, ϕ)分别表示各向同性背景反射系数项与各向异性扰动系数项；I_P，I_S，ρ，Δ_N与Δ_T分别表示纵横波阻抗、密度、法向裂缝弱度及切向裂缝弱度参数；ΔI_P=I_P2-I_P1，ΔI_S=I_S2-I_S1，Δρ=ρ₂-ρ₁，δ_ΔN=Δ_N2-Δ_N1与δ_ΔT=Δ_T2-Δ_T1分别表示反射界面两侧纵横波阻抗、密度及裂缝弱度参数的差值；I_P0，I_S0，ρ₀，Δ_N0与Δ_T0分别表示背景纵横波阻抗、密度及裂缝弱度参数的均值属性，由测井信息或岩石物理实验信息获取；g=I_S0²/I_P0²表示背景横、纵波阻抗的平方比；R_IP=ΔI_P/2I_P0，R_IS=ΔI_S/2I_S0与R_ρ=Δρ/2ρ₀分别表示纵横波阻抗与密度的反射系数.

方程(2)建立了地震波反射系数与背景介质和扰动介质纵横波阻抗、密度及裂缝弱度参数间的定量关系，为裂缝弱度参数反演奠定了公式基础.但同时可看出，各向同性背景反射系数项R_pp^iso(θ)与各向异性扰动项R_pp^ani(θ, ϕ)存在明显的表达差异，而参数表达差异则严重影响弹性参数与裂缝弱度参数反演的精度.为了消除各向同性背景反射系数项与各向异性扰动项的参数表达差异，引入了表征HTI介质各向异性强弱特征的拟裂缝弱度参数(具体见附录B)：

(4)

则各向异性HTI介质纵波反射系数可表示为

(5)

式中，与表示拟裂缝弱度参数反射系数.方程(5)消除了各向同性背景反射系数项与各向异性扰动项的参数表达差异.

假设已有方位叠前道集包含m个入射角，n个方位角，以矩阵形式描述地震记录正演过程：

(6)

式中，d表示地震记录，wvlt表示子波矩阵.上式改为矩阵形式为

(7)

本文采用基于初始模型约束的阻尼最小二乘反演方法，估测得到纵波阻抗反射系数、横波阻抗反射系数、密度反射系数以及拟裂缝弱度反射系数，即：

(8)

式中，m_mod表示待反演参数初始模型，G^T表示矩阵G的转置，σ表示阻尼因子，I表示单位矩阵.本文使用的基于初始模型约束的阻尼最小二乘反演方法依赖于构建的储层弹性参数与各向异性参数的初始模型，若初始模型构建不合理，则求解目标函数耗时较长，容易陷入局部极值.本文构建的初始模型使用测井成像资料和裂缝岩石物理模型估测结果，同时根据实际工区方位叠前地震资料的AVAZ特征属性分析，构建能够有效提高反演结果的初始模型(陈怀震等，2014a).

本文通过求解不同方位的地震反射系数，从而实现拟裂缝弱度参数Δ_N^r与Δ_T^r的反演，进而求取裂缝弱度参数Δ_N与Δ_T.

1.2 贝叶斯框架下岩石物理驱动的储层参数概率地震反演

建立方位叠前地震资料与孔隙度、裂缝密度等裂缝介质储层参数的直接关系难度很大，因此本文研究的裂缝介质储层参数反演并非直接从地震资料进行，而是首先基于AVAZ反演裂缝岩石弹性参数与各向异性参数，并在此基础上进行孔隙度、裂缝密度等裂缝介质储层参数的反演，即，已知参数为反演获得的裂缝岩石弹性参数(纵、横波阻抗与密度)与各向异性参数(裂缝弱度)，待反演目标参数为裂缝介质储层参数(孔隙度与裂缝密度).

(1) 贝叶斯框架下裂缝介质储层参数后验表征

基于贝叶斯理论框架，孔隙度、裂缝密度等裂缝介质储层参数后验条件概率密度可表示为

(9)

式中，R=[ϕ, e]分别表示各向异性储层孔隙度与裂缝密度参数；m= [I_P, I_S, ρ, Δ_N, Δ_T]分别表示裂缝岩石弹性参数与各向异性参数；p(·)表示概率密度函数(PDF)；p(R)表示孔隙度、裂缝密度等裂缝介质储层参数的先验PDF；p(m|R)表示已知裂缝介质储层参数情况下，裂缝岩石参数的条件似然PDF；p(R|m)表示已知裂缝岩石弹性参数与各向异性参数情况下，裂缝介质储层参数的后验PDF；α=p(R)p(m|R)表示裂缝岩石参数的全概率，对p(R|m)起到正则化因子的作用.

(2) 裂缝介质储层参数先验分布的构建与估计

孔隙度、裂缝密度等裂缝介质储层参数先验PDF的构建需对测井资料进行统计分析获知其大体分布，然后基于先验认知建立与测井储层参数较为相近的先验分布.由于混合高斯分布的可解析性，假设反演目标储层参数服从混合高斯(Gaussian Mixture，GM)分布，即不同分量的高斯分布的线性组合：

(10)

式中，N表示先验分布中高斯分量的个数；R=[ϕ, e]表示各向异性储层孔隙度与裂缝密度参数，可视为随机变量；β_j表示先验分布中第j个高斯分量所占权重，满足；μ_R^j与Σ_R^j分别表示第j个高斯分量的均值和方差；N(·; μ, Σ)表示均值为μ、方差为Σ的高斯分布函数.

采用期望最大化(EM)算法估计每一项高斯分量的均值μ_R^j、方差Σ_R^j及所占权重β_j，从而构建孔隙度、裂缝密度等裂缝介质储层参数的先验分布.

(3) 统计性岩石物理模型的构建

地下储层条件是复杂多变的，确定性的岩石物理模型难以准确描述裂缝岩石弹性参数及各向异性参数与孔隙度、裂缝密度等各向异性介质储层参数之间的关系.确定性岩石物理模型添加随机误差构成了统计性岩石物理模型，可用于弥补确定性岩石物理模型本身精确度的缺陷，描述测井上没有体现的所有可能的储层条件实现(Mukerji et al., 2001；Bachrach et al., 2006).若将所有变量看为随机变量，统计岩石物理模型可表示为

(11)

式中，f_RPM(·)表示裂缝岩石参数与各向异性介质储层参数之间确定性的关系函数，即确定性岩石物理模型；ε是用来描述模型精确度的随机误差，可通过确定性岩石物理模型与实际测井资料之间的相对差异估计获得，通常定义其服从截断高斯分布.为了寻求一种与实际裂缝型各向异性储层资料匹配最佳的岩石物理模型，本文通过构建裂缝型各向异性岩石物理等效模型，用于表征裂缝岩石参数与裂缝介质储层参数之间的关系，具体构建步骤可参考相关文献(张广智等，2013).基于估测的Thomsen各向异性参数，可计算相应的法向及切向裂缝弱度参数，两者之间近似的参数关系式为

(12)

式中，ε^(V)，δ^(V)与γ^(V)或γ表示HTI介质的Thomsen各向异性参数.

典型的岩石物理模型都是基于测井数据的，是深度域模型，但地震反演得到的裂缝岩石参数是时间域的，其分辨率低于测井数据.在实际应用中，需考虑由资料信息来源不同造成的有效数据尺度和数值范围的影响(袁成等，2016).使用Backus平均，实现测井资料与地震资料尺度的一致性处理.

(4) 联合分布的构建与估计

MCMC方法是一种启发式的全局寻优算法，基于MCMC抽样可以得到大量来自于后验概率分布的样本，不仅可以得到每个未知参数的估计值，而且可以得到与之相关的各种不确定性信息(Tarantola，2005；张广智等，2011).通过MCMC抽样，可从裂缝介质储层参数先验PDF中生成多个储层参数的采样样本，构成裂缝介质储层参数的先验分布随机样本空间.对裂缝介质储层参数的先验样本空间采用确定性岩石物理模型，即构建的裂缝型各向异性岩石物理等效模型，从而计算每个储层参数样点所对应的裂缝岩石弹性参数与各向异性参数，并添加相应的随机误差，进行横向延拓，得出裂缝岩石参数的随机分布样本空间.已知储层参数的情况下，裂缝岩石弹性参数与各向异性参数的条件PDF p(m|R)可表示为

(13)

式中，μ(R)=f_RPM(R)表示裂缝岩石参数的均值；Σ_ε表示随机误差项的方差，并假设其只与随机误差ε有关，而与储层参数R无关.

在实际应用中，通过MCMC抽样方法对裂缝介质储层参数先验PDFp(R)抽取一系列样点，得到储层参数随机样本{ R₁, R₂, …, R_Ns}，其中N_s表示抽样样本的数量.利用确定性岩石物理模型f_RPM(R)求取储层参数R_k(k=1, 2, …, N_s)对应的裂缝岩石参数均值μ(R_k)：

(14)

将式(14)代入式(13)，可得到裂缝介质储层参数所对应裂缝岩石弹性参数与各向异性参数的条件概率分布，使用MCMC抽样方法对其进行采样，获取储层参数样本空间对应的裂缝岩石参数样本空间分布情况m_k(k=1, 2, …, N_s).裂缝介质储层参数的先验分布空间样本及其对应的裂缝岩石弹性参数与各向异性参数空间样本构成了联合分布空间样本.假设联合分布p(m, R)亦可视为GM分布，同样采用EM算法估计该联合分布的参数，则其可表示为

(15)

式中，

(16)

π_j表示联合分布中第j个高斯分量所占权重；μ_{[m, R]}^j与Σ_{[m, R]}^j分别表示联合分布第j个高斯分量的均值和协方差矩阵.

(5) 裂缝介质储层参数后验分布的解析计算

与联合分布p(m, R)类似，后验PDF p(R|m)同样可视为GM分布，其高斯分量的参数可通过先验分布及联合分布估计的参数解析求解.因此，裂缝介质储层参数后验条件概率p(R|m)可表示为

(17)

式中，

(18)

η_j表示后验PDF中第j个高斯分量所占权重；μ_R|m^j与Σ_R|m^j分别表示后验PDF第j个高斯分量的均值和协方差矩阵.

根据贝叶斯估计方法，寻找后验PDF p(R|m)的最大值所对应的裂缝介质储层参数R，该储层参数便可视为反演结果.

2 模型试算与实际资料处理 2.1 模型试算

如图 1所示，取经Backus平均的实际裂缝型碳酸盐岩储层工区井A上深度范围为5150~5250 m的一段作为地层模型，图中曲线自左至右分别表示孔隙度、裂缝密度、纵横波阻抗、密度及法向、切向裂缝密度参数.图 2是基于构建的裂缝型各向异性岩石物理模型计算的储层弹性参数(纵横波阻抗、密度)与各向异性参数(法向、切向裂缝密度参数)结果，黑色实线表示原始测井值，红色虚线表示计算值，结果表明构建的确定性岩石物理模型能够很好的表征裂缝岩石弹性参数、各向异性参数与裂缝介质储层参数之间的岩石物理关系，两者之间的准确性较高，验证了所构建的确定性岩石物理等效模型的合理性.由于岩石物理建模的不确定性，裂缝岩石弹性参数与各向异性参数计算值与实际值存在相应差异.为了削弱此差异，给确定性岩石物理模型添加一项随机误差，图 3为随机误差的直方分布，随机误差取均值为零的高斯截断分布，方差则由确定性岩石物理模型的计算值与实际测井值的差异大小决定.

为了建立与实际测井资料吻合度较高的储层参数先验分布，假设孔隙度与裂缝密度均服从三项高斯分量的GM分布，采用EM算法分别估计孔隙度与裂缝密度的GM分布参数，并采用MCMC抽样方法进行采样.图 4b与图 5b分别表示孔隙度与裂缝密度参数的GM先验分布样本，与测井先验信息图 4a与图 5a吻合度较高.基于构建的统计性岩石物理模型，求取孔隙度与裂缝密度先验分布样本中每一个样点所对应的裂缝岩石弹性参数(纵横波阻抗及密度)与各向异性参数(裂缝弱度)的随机分布样本空间，并与孔隙度及裂缝密度等裂缝介质储层参数样本相结合，获得联合分布样本空间.图 6与图 7分别是孔隙度及裂缝密度与裂缝岩石参数的联合分布概率图.基于EM算法分别估算孔隙度及裂缝密度与裂缝岩石参数联合分布中高斯分量的均值、协方差及所占权重，根据公式(18)与(19)便可计算相应的孔隙度及裂缝密度的后验PDF.图 8a及图 9a分别表示孔隙度及裂缝密度的后验PDF分布图，图中红色调表示孔隙度及裂缝密度的后验条件概率较大，而蓝色调则表示后验条件概率较小.根据式(9)所示的裂缝介质储层参数反演目标泛函，估计孔隙度及裂缝密度后验条件概率的最大值，并寻找该最大值所对应的孔隙度值及裂缝密度值，即为最终的反演结果.图 8b及图 9b分别表示孔隙度及裂缝密度的反演结果，结果表明反演的孔隙度及裂缝密度曲线与实际测井资料数据吻合度较高，验证了该反演方法的可行性.

2.2 实际资料处理

实际方位地震道集来自于中国陆上某勘探工区，进行叠前反演前，需对地震数据进行去噪、保幅等预处理.利用本文提出的方法对方位叠前道集进行了实际资料的试处理.图 10表示四个方位、三个角度的角度叠加道集.在实际工区裂缝介质储层参数的反演应用中，首先基于推导的HTI介质纵波反射系数方程及相应的AVAZ反演方法，使用上述的方位叠前道集反演获取实际工区裂缝岩石的弹性参数及裂缝弱度参数，图 11表示工区弹性参数及裂缝弱度反演结果.然后基于反演获得的裂缝岩石弹性参数及各向异性参数数据体，计算裂缝岩石参数已知条件下孔隙度及裂缝密度等裂缝介质储层参数的后验条件概率.根据贝叶斯估计方法中的最大概率估计方法，从后验条件概率体中搜索出概率最大值，并寻找最大后验条件概率所对应的裂缝介质储层参数值，最终反演出裂缝介质储层参数体，即孔隙度与裂缝密度.图 12及图 13分别表示过井A的孔隙度及裂缝密度反演结果，红色虚线椭圆为裂缝型气藏发育区，反演的孔隙度及裂缝密度均呈现明显的高值，与测井解释结果吻合得较好，并与实际钻遇结果一致，验证了方法的有效性.

3 结论

在裂缝岩石弹性参数与各向异性参数AVAZ反演的基础上，提出基于裂缝型各向异性统计岩石物理模型的贝叶斯框架下岩石物理驱动的裂缝介质储层参数概率地震联合反演方法，综合应用MCMC抽样方法和EM算法，实现孔隙度、裂缝密度等各向异性介质储层参数后验条件概率的表征.该方法的优势主要体现在以下方面：

(1) 引入了表征HTI介质各向异性强弱特征的拟裂缝弱度参数，推导了相应的纵波AVAZ反射系数近似公式，消除了各向同性背景反射系数项与各向异性扰动项的参数表达差异，基于AVAZ反演实现了裂缝岩石弹性参数与各向异性参数的有效估测；

(2) 联合岩石物理、测井及方位地震资料实现了孔隙度与裂缝密度等多参数概率地震联合反演，获得了裂缝介质储层参数的后验概率表征，能够定量描述反演结果的误差；

(3) 在测井资料有限的情况下也能获得很好的反演结果，克服了神经网络反演方法及地质统计学随机模拟方法对测井资料过度依赖的问题.

测井资料及实际方位地震处理结果表明，该方法能够稳定合理地从方位地震资料中获取裂缝岩石弹性参数与各向异性参数，是一种较为可靠的孔隙度、裂缝密度等各向异性介质储层参数概率地震反演方法.

附录A 各向异性HTI介质的纵波反射系数方程推导

基于Born近似，各向异性HTI介质纵波反射系数可表示为(Shaw and Sen, 2006)

(A1)

式中，r₀是水平界面一点，而S(r₀)表征反射函数，其表示为

(A2)

式中，

(A3)

p和t分别表示慢度矢量和极化矢量.ξ与η_mn的具体表达式见Shaw和Sen(2006).弱反射的各向异性HTI介质在均匀各向同性上的扰动刚度矩阵ΔC可表示为

(A4)

式中，λ与μ表示拉梅常数.将式(A3)与(A4)代入式(A2)中，则HTI介质纵波反射函数可表示为

(A5)

式中，α₀²表示背景纵波速度的平方.将式(A5)代入式(A1)中，则HTI介质纵波反射系数方程为

(A6)

当介质为均质各向同性介质时，方程(A6)退化为Fatti反射系数公式(Fatti et al., 1994).

附录B 拟裂缝弱度参数化方法

引入拟法向裂缝弱度参数Δ_N^r，并假定R_ΔN^r=δ_ΔN，即，也即

(B1)

式中，下标0表征拟法向裂缝弱度参数的均值属性.对方程(B1)左右两侧进行积分，得

(B2)

也即，

(B3)

式中，C是一常数.拟法向裂缝弱度参数Δ_N^r可写为法向裂缝弱度参数Δ_N的一阶扰动形式，即

(B4)

式中，m与n均为两个未知量.将方程(B4)代入方程(B3)可得：

(B5)

同样，将方程(B5)代入方程(B4)可得：

(B6)

取C=1，则

(B7)

同样地，对切向裂缝弱度参数Δ_T作类似的变形处理，引入拟切向裂缝弱度参数Δ_T^r，则

(B8)