0 引言

近年来，利用地震方法进行岩石各向异性的研究受到广泛关注。裂缝的存在导致地下介质呈现各向异性，在适当的地质条件下，这些裂缝会成为潜在的油气聚集场所。四川盆地下三叠统飞仙关组鲕滩气藏是重要的产层之一，也是近几年勘探开发的热点层位^[1-2]。飞仙关组储集岩主要为鲕粒灰岩、鲕粒云岩及粉晶云岩，储层类型以裂缝-孔隙型为主，根据岩心观察及成像测井资料可知，研究区主要分布高角度或者垂直裂缝^[3]。具有水平同相轴的横向各向同性(HTI)介质可以用来等效表征研究区的裂缝孔隙型储层，根据HTI介质理论，我们可以依靠不同方位的叠前地震资料来检测裂缝。Ruger和Tsvankin等^[4-5]通过研究认为在HTI介质中可用反射振幅随偏移距变化(AVO)梯度来描述裂缝的走向信息，并提出结合方位AVO数据、正常时差校正(NMO)所用的速度和横波分裂分析来计算各向异性参数。随后围绕着各向异性介质中的NMO速度、AVO梯度和层间旅行时等属性在裂缝预测中的应用日趋火热。Gray等^[6]利用方位数据实现裂缝密度和走向的直接反演；Bachrach等^[7]对裂缝储层的弹性参数和各向异性参数进行重构后利用叠前反演预测裂缝参数。前人为了模拟裂缝提出了多种裂缝等效介质理论。其中，镇晶晶等^[8]通过对比不同的裂缝等效介质理论，认为线性滑动等效介质理论用法向和切向弱度来表征裂缝具有较好的应用效果。Shaw等^[9]根据Born一阶近似并在反射界面之上为各向同性介质假设的前提下，提出了利用最小二乘反演方法来计算法向和切向裂隙弱度的方法。Alhussain等^[10]在上覆介质是各向同性岩石假设的裂缝储层中反演了裂隙弱度和裂缝密度，并与基于方位属性的裂缝预测结果进行了对比，结果表明利用振幅随偏移距和方位角变化(AVOA)的裂缝密度反演效果较好。Narhari等^[11]利用Ruger方程实现了分方位叠前地震反射系数随入射角和方位角的变化(AVAZ)反演在深层碳酸盐岩储层裂缝预测中的应用。卢明辉等^[12]基于Ruger提出的HTI介质模型P波反射系数公式，引入小生境遗传算法，得到了较高精度的Thomsen各向异性参数。马中高等^[13]利用多方位速度和AVOA同步分析，通过不断调整速度求取了相应的AVO参数，并将AVO参数用于速度分析中的相似系数计算，确定了最佳的速度与AVO参数，该方法使得各向异性速度分析更精确。薛姣等^[14]根据裂隙弱度差等于一阶差分矩阵与裂隙弱度的乘积，推导出了反射振幅与裂隙弱度的直接关系式，并将裂隙弱度非负值约束条件加入裂隙弱度反演方程中，改进了AVOA裂缝密度反演方法。赵才顺等^[15]基于叠前纵波地震资料的方位各向异性特性进行裂缝预测，并通过岩石物理正演模拟实验以及成像测井资料确定裂缝方位和纵波反射系数的关系，提高了裂缝预测的精度。吴珊珊等^[16]以Ruger近似式为基础建立模型，将模型中各向异性参数和岩性参数设置为概率密度分布函数，通过Monte Carlo随机方法进行叠前AVAZ正演模拟拟合椭圆，进而寻找最佳的裂缝参数。李博南等^[17]通过频变AVO理论建立目标函数并使用全局最优化算法反演裂缝参数，考虑了反射系数频变响应对裂缝密度和时间尺度因子的敏感性，有效描述了裂缝型储层。罗腾等^[18]针对HTI介质模型研究了裂缝型储层的各向异性参数地震振幅随方位角变化的反演方法，并在四川盆地页岩储层中取得了较好的应用效果。

纵波各向异性裂缝预测方法作为目前应用比较广泛、效果较好的一种裂缝预测方法，大量国内外学者都进行了深入的探索，其中最早出现的较为经典的是传统的椭圆拟合法^[19-21]。该方法原理简单、易于实现，被广泛应用于各种商业化软件中，并成功用于油田实际资料的解释中。但是，传统椭圆拟合方法的准确性受到一系列因素的影响，当入射角较大、信噪比较低、方位角分布不均匀、盖层物性和各向异性横向变化、上覆介质透射各向异性较强时，该方法的裂缝预测结果误差较大。在实际的裂缝预测中，由于各种因素的影响，裂缝的密度及发育方向很难准确估计，尤其是裂缝密度的定量预测。传统椭圆拟合方法由于受到各种因素的制约，不能准确得到裂缝密度和方向，而后发展的基于三项方程的反演方法虽然弥补了椭圆拟合法的多种不足，但其仍不能实现裂缝密度的定量预测。因为Ruger方程能够反映各向异性参数对地震响应特征的影响，后续理论和技术主要是对该方程针对研究区实际情况的修编。因此，本次研究推导了一种根据Ruger方程精确式计算储层裂缝密度的方法，以提高裂缝定量预测的精度。

1 方法原理 1.1 裂缝线性滑动等效介质理论

为了建立裂缝型储层的岩石物理模型，准确识别裂缝中流体的性质，国内外很多学者开展了大量的研究。其中Schoenberg^[22]构建了大尺度裂缝线性滑动模型。

Hsu等^[23]提出用法向和切向各向异性系数来表征裂缝，主要是为了便于研究裂缝线性滑动等效介质条件下储层的弹性性质。裂缝法向和切向各向异性系数的值域范围在0~1之间，裂缝线性滑动模型的刚度矩阵为

其中：

式中：λ和μ为不含裂隙各向同性岩石的拉梅参数；r=λ/(λ+2μ)；Δ_N和Δ_T分别是法向和切向各向异性系数。

可以利用切向各向异性系数Δ_T计算裂缝密度e：

在两个各向异性介质的分界面处，地震波会产生反射纵横波和透射纵横波，且反射和透射系数也会随着炮检距(入射角)和方位角的变化而变化。由于反射和透射系数的精确解析表达式非常复杂，一般很难直接用于AVO计算，为此Ruger^[24]通过各种假设在对称轴方向一致的HTI各向异性介质分界面储层中用摄动法推导出了随入射角和方位角变化的纵波反射系数：

式中：R_P^HTI(i, φ)为HTI介质中的纵波反射系数；i为入射角；φ为方位角；ΔZ为界面上下的波阻抗之差；Z为平均波阻抗；α为纵波速度；β为横波速度；G=ρβ²为横波切向模量(ρ为介质密度)；Thomsen各向异性参数ε为纵波各向异性指数、δ为纵波各向异性变异系数、γ为横波各向异性指数；φ_sym为HTI介质对称轴方向的方位角；Δ[·]表示上、下界面物理量之差；[·]表示上、下界面物理量的均值。

式(3)即Ruger方程精确式。由式(3)可以得到不同参数条件下HTI介质的纵波反射系数随入射角和方位角的变化情况，可进一步探讨和研究HTI介质的各向异性。

1.2 基于Ruger方程精确式的裂缝预测方法

椭圆拟合法是指当入射角较小时，通过忽略Ruger方程精确式中的高阶项部分，进一步简化Ruger方程，以方便得到随方位角变化的AVO梯度项，它们的具体表达式分别为：

式中：φ_k为第k次观测到的地震方位角；B(φ_k)为随φ_k变化的AVO梯度项；B^iso、B^ani分别为B(φ_k)的各向同性部分和各向异性部分。事实上，式(4)可以近似为一个椭圆，并且随着地下介质中裂缝强度的增加，拟合出的方位椭圆的扁率也会增大，可以用长轴或短轴方向来代表裂缝的走向。因此，理论上只需知道3个或3个以上方位的反射地震数据，就可以实现对目的地层中任意一点裂缝发育密度和方位的预测。

根据HTI介质中各向异性参数与裂缝岩石物理参数之间的关系(式(2))，对Ruger方程精确式(式(3))进行变换，得到纵波反射系数简化方程：

其中：

式中，g=μ/(λ+2μ)，为横纵波速度比的平方。

将Ruger方程精确式所述的方位角道集及其对应的纵波反射系数代入式(7), 得到矩阵方程

其中：

式中：R_k=R_Pk^HTI(i_k, φ_k)为第k次观测到的方位角道集对应的纵波反射系数；k∈[1，n]；n为方位角道集的数量。

根据矩阵方程(8)可建立用于反演x的最优化目标函数：

对式(10)进行Tikhonov正则化处理后得

式中，τ为正则化参数。

计算经过Tikhonov正则化处理后最优化目标函数的最小二乘解x =(A^TA+τI)^-1A^Tb，即得Δ_T。

依据式(2)，切向各向异性系数与裂缝密度之间的关系可表示为

实际计算中Δ_T与e之间的关系是根据Hudson币状裂隙模型^[25]及Schenberg线性滑动模型^[26]确定的。这样得到的裂缝密度既不受地层各向异性、地应力各向异性等因素的影响，也不会因裂缝充填流体类型的变化而变化，它是地下裂缝真实的定量描述，是常规裂缝预测方法所不能做到的。同时，由于在推导反演方程时保留了Ruger方程精确式中的高阶项，因此本次研究的计算公式(式(12))比传统的裂缝反演方法和椭圆拟合法更加精确。

当忽略反射系数公式——Ruger方程精确式中的高阶项时，可以得出拟合椭圆扁率的近似式：

可以看出，拟合出椭圆的扁率不仅与界面上下介质的各向异性参数有关，还与界面上下介质的纵横波速度、纵波速度差、纵波阻抗差、横波切向模量差等因素有关。

假定Thomsen各向异性三参数可以代表地下的裂缝密度。计算得到单一Thomsen各向异性参数中的ε、δ及γ对拟合椭圆扁率的影响，如图 1所示。从图 1可以得出两点结论。1)介质的各向异性越强(即地下介质的裂缝密度越大)时，拟合出的各向异性椭圆的扁率越大，这与通常的认识相吻合。但是，本实例只研究了单一变量对拟合椭圆扁率的影响，当多个变量同时变化时，情况将变得十分复杂，此时拟合椭圆的扁率不一定随着介质各向异性的增强而变大。2)当介质的各向异性参数为正各向异性和弱各向异性时，曲线更加接近线性化。因此，在正各向异性和弱各向异性时的预测结果要好于负各向异性和强各向异性两种情况，同时这也决定了使用椭圆扁率只能定性而非定量代表地下的裂缝密度。

图 2为盖层弹性参数与拟合椭圆扁率的交会图。从图 2中可以看出，盖层的弹性参数变化会影响拟合椭圆扁率的变化，因此椭圆扁率大并不一定是由各向异性引起的，进行裂缝预测时还应消除盖层弹性参数的影响。

由于各向异性梯度项B^ani是各向异性参数ε、δ、γ的线性函数，这在一定程度上要优于用拟合椭圆扁率来表征裂缝密度的情况。同时从式(6)可知，虽然各向异性梯度项也会受到盖层速度的影响，但利用Thomsen各向异性三参数表征地下裂缝密度消除了盖层密度和纵波阻抗等因素的影响。因此相较于利用拟合椭圆扁率来表征裂缝密度而言，利用各向异性梯度项来表征裂缝密度在很大程度上具有优势。

截至目前，很多学者开展了裂缝岩石物理模型的研究，其中较为常用的包括Hudson^[25]提出的薄币状裂隙模型和Schenberg^[26]提出的线性滑动模型。通过研究两个模型的关系，得到岩石物理参数中法向各向异性系数和切向各向异性系数与裂缝参数之间的关系：

式中：K′和μ′分别为裂隙中充填物的体积模量和剪切模量；χ为裂缝的纵横比。Bakulin等^[27]又推导了HTI介质中各向异性参数与裂缝岩石物理参数之间的关系：

这样就建立了裂缝参数(裂缝密度、裂缝纵横比及裂缝充填流体类型)和Thomsen各向异性参数之间的关系。

图 3给出了裂缝饱含水、含油水混合物和含气水混合物3种情况下，裂缝岩石物理参数随横纵波速度比和裂缝密度变化的特征。从图 3可以看出，在这3种情况下，Δ_N和Δ_T都随e的增大呈上升趋势。分析图 3a、c、e中Δ_N的曲线变化可发现，裂缝含气、含油水和含水时，Δ_N的变化幅度均有差异。分析图 3b、d、f可以发现，不论裂缝含气、含油水还是含水，Δ_T的变化幅度不明显；因此我们分析认为Δ_T的变化主要是受裂缝密度影响。如果能从地震数据中计算得到Δ_T，则可以直接预测岩石的裂缝密度。由于Δ_N不仅受裂缝密度的影响还受充填物类型的影响，因此计算得到Δ_N后结合裂缝密度的预测结果，则可以进一步预测裂缝中充填物的类型，即可以用Δ_N做为裂缝流体指示因子为裂缝型储层的勘探和开发提供有利的技术支持。

2 应用实例 2.1 基于理论模型的裂缝参数反演

可利用叠前方位数据对基于岩石物理的叠前裂缝定量预测方法进行测试。采用HTI介质来模拟实际地层中的裂缝分布情况。在合成方位地震数据时，使用了20°、30°和40° 3个入射角，假设每个入射角生成5个方位角数据，这5个方位角分别为18°、54°、90°、126°和162°。此时，可以得到15个方位角道集，代表每个共深度点道集(CDP)点具有15次覆盖。图 4为工区中部发育一条东西向展布裂缝发育带的模型，图 5为工区内裂缝随机分布的模型。从图 4和图 5可以看出：利用本文提供的预测方法反演得到的结果与真实模型匹配，表明了本文方法的可靠性和准确性；并且，实例用于反演的数据中含有40°入射角的数据，由于在推导反演方程中没有忽略高阶项，因此本文方法在入射角较大时仍然适用，得到的反演结果也依然很精确。

2.2 基于实际地震数据的裂缝参数反演

四川盆地西北部三叠系飞仙关组储层非均质性强，目前发现的储层致密低孔，以裂缝型为主。已钻井证实四川盆地西北部正坝南地区飞仙关组鲕滩储层为裂缝-孔隙型，大部分钻井鲕滩储层裂缝较发育^[28-29]。

为了使叠前反演的结果更加有效和可靠，选择部分角度的叠加来提高地震资料的信噪比。同时裂缝走向的判别主要利用井信息或地震波属性的椭圆拟合来获取。根据与裂缝方向垂直和平行的关系，从实际地震数据中抽取不同入射角对应的角道集。本次研究中选择部分角度叠加数据体的入射角范围是2°~8°、8°~14°、14°~20°、20°~26°和26°~32°，它们分别代表入射角为5°、11°、17°、23°和29°的角道集。反演的方位AVO数据中舍弃了部分近偏移距数据。

研究区的LJ5井在飞仙关组取心70.2 m，岩心观察可见1 259条裂缝，该井从岩心描述来看，裂缝主要以高角度裂缝或垂直缝为主，裂缝较发育。利用过LJ5井的地震纵波AVOA数据反演法向和切向各向异性系数计算裂缝密度，并与测井和钻井信息相对比。

图 6是法向和切向各向异性系数反演结果以及利用切向各向异性系数计算的裂缝密度。由于切向各向异性系数仅与裂缝密度成正比，裂缝中的充填物质不会影响切向各向异性系数，因此切向各向异性系数较大的地区往往对应裂缝密度较大的地区。从切向各项异性系数的反演结果(图 6a)来看，在钻井LJ5位置Δ_T较大，为0.15~0.20。

法向各向异性系数不仅受裂缝密度的影响，还受裂缝中充填物的影响，特别是当裂缝中含气后，法向各向异性系数会增大。从法向各向异性系数反演结果(图 6b)来看，LJ5井所在的位置飞仙关组顶部Δ_N较大，因此判断该井在飞仙关组顶部不仅裂缝较发育，裂缝中含气性也较好。后在距LJ5井飞仙关组顶部100 m的埋深为2 350~2 375 m处测试获得日产气75.65×10⁴m³的高产工业气流。

图 7为利用切向各向异性系数计算的裂缝密度反演结果。从图 7可见，LJ5井所在的位置裂缝密度较大，表明LJ5井裂缝较发育；这与实际取心结果也比较吻合。

研究区的LJ2井在飞仙关组取心111.4 m，岩心观察仅见271条裂缝，裂缝欠发育(图 8)。从图 8可见：LJ5井所在的位置预测的裂缝密度较大，裂缝较发育；LJ2井所在的位置裂缝密度小，裂缝欠发育。反演结果与实际取心结果吻合。

3 结论

1) 本次工作基于岩石物理叠前裂缝预测方法定量预测裂缝密度，消除了其他各向异性对裂缝密度的干扰。在测井和岩石物理分析前提下，能够进一步计算得到裂缝流体指示因子，这是常规裂缝反演所不能做到的，特别是对含气裂缝有较强的敏感性，从而为后期的开发布井给予一定的指导，可降低油气藏预测的风险。

2) 利用四川盆地正坝南地区实际纵波AVOA数据实现了法向和切向各向异性系数反演，进而利用基于Ruger方程精确式的裂缝预测方法反演了裂缝密度，预测结果与钻井取心资料吻合。因此在裂缝储层发育区可以利用基于裂缝岩石物理模型的AVOA反演方法开展优质储层预测。