地层裂缝参数研究在油气田勘探开发、水利工程、土木工程等领域具有重要意义。随着全球油气资源勘探程度的不断提高以及需求量的不断增大, 油气勘探开发从常规油气储层转向致密砂岩、复杂碳酸盐岩、火山岩、页岩等复杂油气藏, 这些油气藏都涉及到裂缝参数预测。在裂缝型油气藏中, 由于定向构造应力作用, 自然裂缝通常呈现一种定向排列现象, 裂缝体的走向、发育密度、充填物等裂缝参数与油气运移、聚集、成藏、储量有着密切的联系, 另外裂缝参数也是储层工程“甜点”的主要评价参数之一。所以储层裂缝参数的定量预测和评价是油气勘探开发的重要环节之一。

目前, 裂缝型储层特征分析和预测技术可以分为三大类, 分别是地质岩心综合分析预测法、测井识别技术以及地震预测技术。地质岩心综合分析法主要是通过地表露头、岩心直接观察裂缝的发育状态, 统计裂缝参数, 进而预测、标定与之相对应地层的裂缝分布^[1-5]。这种技术简单、直接、高效, 但露头容易受风化等外动力地质作用的影响而改变其原始地层赋存状态, 且露出地表, 原始地层应力已经释放, 裂缝状态也相应发生了变化。岩心裂缝状态会受到取心机械外力的影响, 且随赋存应力变化而发生变化, 不能直接反映地层裂缝分布, 而且取心率低, 难以覆盖全井段。测井技术识别裂缝分布, 最初是利用常规测井曲线在裂缝层的异常响应或者基于非线性理论的多曲线联合组构新特征参数曲线来识别井周裂缝参数^{[3, 6]}, 后来随着成像测井技术的发展, 可以获取能够直观、真实地刻画井壁周围地层裂缝产状、开度、密度等参数的图像^[7-9]。近些年发展的偶极声波远探测技术在裂缝、断层的探测中也有着良好的应用效果, 使测井探测实现了从近井筒到远井筒的突破, 探测直径提高到近百米^[10-13]。测井技术具有较高的分辨率, 是原始地层状态的真实反映, 但测井数据易受到钻井液、井壁坍塌等因素的干扰, 预测精度受到了一定的影响, 成像测井成本高、解释结果的主观性大, 且测井识别技术难以准确地预测整个工区的裂缝分布。地震预测虽然精度不如测井技术, 但在探测整个工区或者深部地层裂缝赋存状态方面, 有着先天的优势^[14-18]。

裂缝参数地震预测技术主要包括三维地震属性预测技术^[19-25]、多波多分量预测技术^[26-27]、纵波方位各向异性预测技术^[28-35]。利用三维地震属性预测地层裂缝参数, 即利用地震数据体的一阶导数——倾角^[19]、二阶导数——曲率^{[21, 25, 36-37]}以及地震波相邻道相似性——相干体^[38-40]等几何属性实现地层裂缝刻画。反射面曲率代表对应地层的变形程度, 某一层位曲率值越大, 代表该层位褶皱程度越大, 通常裂缝越发育, 即曲率预测裂缝的基本假设是认为地层构造程度与裂缝发育程度之间存在单调映射关系, 因此对于裂缝发育的构造空白区, 曲率属性则无法应用。相干属性发展主要经历了三代, 基于三道互相关的C1算法^[38], 基于多道相似性的C2算法^[38], 基于本征结构的相似算法C3^[40]。C3算法较前两种算法无论在分辨率上, 还是抗噪性上都有着明显的优势, 但其计算速度较慢。后来许多学者提出了一系列的改进算法, 如基于梯度结构张量倾角估计算法^[41]、基于方位叠前地震数据的相干算法^[42]等。多波多分量预测技术, 主要利用裂缝储层中的横波分裂现象^[27]实现裂缝参数的预测, 但高昂的采集处理成本, 纯横波勘探激发难等问题极大地限制了其推广应用。转换PS波(纵-横转换波)勘探使得多分量勘探有了进一步的发展, 但相较于纵波勘探, 依旧存在高成本、对地震数据质量要求高、对裂缝充填物不敏感等问题。P波(纵波)方位各向异性裂缝预测技术主要利用裂缝诱导的纵波速度、振幅、频率等属性在观测方位角域的变化规律实现裂缝密度、方位、填充物等参数预测, 包括基于纵波运动学特性的预测技术和基于纵波动力学特性的预测技术。前者有NMO(normal moveout)速度方位变化裂缝预测技术^[43-46]、正交地震测线纵波时差预测技术^[47-49]和基于VSP数据的各向异性参数预测^[50]等, 这3项技术均是基于速度在方位角域和偏移距(入射角)域的规律变化实现裂缝参数预测, 但分辨率不足是其最难解决的问题。后者主要利用地震P波振幅^[51]、阻抗^[52-54]、衰减^[55-56]等属性随方位角的变化实现储层预测, 具有更高的分辨率。整体来讲, 地震预测技术的探测深度大、范围广等特征是地质岩心技术和测井技术所不具备的, 而且基于P波属性随入射角和方位角变化(五维地震数据)的技术较其它裂缝地震预测技术具有更高的分辨率和准确度, 可获得更多的储层信息且纵波数据更易获取。与此同时, 裂缝储层的等效地震岩石物理理论是连接储层微观参数与宏观地震响应之间的桥梁, 也是反射系数参数化和裂缝参数稳定预测的理论基础。所以我们将从裂缝型储层等效岩石物理理论和五维地震数据裂缝预测技术两个方面介绍岩石物理驱动的裂缝预测技术研究现状与进展。本文主要综述裂缝型储层等效各向异性岩石物理理论的研究现状与进展。五维地震数据裂缝预测技术研究现状与进展将在第二部分(Ⅱ)中详细介绍。

1 裂缝储层等效各向异性岩石物理理论

地震岩石物理模型建立了地层岩石物性参数和弹性参数之间的关系, 是地震反演和裂缝参数预测的桥梁, 所以了解岩石物理裂缝模型, 厘清各模型的理论前提和物理机理是利用地震资料预测裂缝参数的基础^[57]。所以, 本节将从3个方面综述裂缝型储层地震岩石物理理论: ①考虑孤立孔缝的地震岩石物理理论, 包括Hudson理论^[58-62]、Schoenberg线性滑动理论^[63-68]、Eshelby-Cheng模型^[69]、Aniso-DEM模型(anisotropic differential effective model)^[70]、Aniso-SCA模型(anisotropic self-consistent approximation model)^[70]等; ②衰减各向异性地震岩石物理理论, 包括修正Shoenberg线性滑动模型^[71-75]、Chapman多尺度裂隙模型^[76-77]等; ③极限频率下考虑孔缝连通性的地震岩石物理理论, 包括Thomsen含等径孔隙裂隙模型^[78]、基于广义Gassmann方程的饱和裂缝岩石模型^[79-80]等。

1.1 考虑孤立孔缝的地震岩石物理理论

以静力学为基础导出的裂缝型储层地震岩石物理理论, 主要研究不考虑孔隙-裂隙或裂隙-裂隙之间流体流动的含孤立孔缝岩石。该理论将裂缝作为岩石的包含物, 采用等效平均思想, 推导参考元的等效岩石模量。下面重点介绍岩石弹性模量解析表达的Hudson理论和Schoenberg线性滑动理论。

1.1.1 Hudson理论

如图 1所示, HUDSON^[58-59]基于平均波场散射理论研究了各向同性背景中包含一组水平硬币型裂缝介质的等效弹性刚度矩阵C, 提出了Hudson理论(本文坐标系统均采用Cartesian坐标系, 并规定z或x3轴正方向垂直向下, n为裂隙的法向单位向量):

$ \boldsymbol{C}=\boldsymbol{C}_{0}+\boldsymbol{C}_{1}+\boldsymbol{C}_{2}+o\left(e^{2}\right) $

(1)

式中: C为裂缝介质等效刚度张量; C₀为各向同性背景刚度张量; C₁, C₂分别是由裂缝引起的刚度张量的一阶校正项和二阶校正项; o(e²)为裂隙密度e²的高阶无穷小量。刚度矩阵C的非零元素分别为:

$ \begin{gathered} C_{11}=\lambda_{b}+2 \mu_{b}-\frac{e}{\mu_{b}}\left(\lambda_{b}+2 \mu_{b}\right)^{2} U_{33}+ \\ \frac{e^{2}}{15}\left(\lambda_{b}+2 \mu_{b}\right) q U_{33}^{2} \end{gathered} $

(2a)

$ C_{33}=\lambda_{b}+2 \mu_{b}-\frac{e}{\mu_{b}} \lambda_{b}^{2} U_{33}+\frac{e^{2}}{15} \cdot \frac{\lambda_{b}^{2} q}{\lambda_{b}+2 \mu_{b}} U_{33}^{2} $

(2b)

$ C_{44}=\mu_{b}, C_{66}=\mu_{b}-\frac{e}{\mu_{b}} \mu_{b}^{2} U_{11}+\frac{e^{2}}{15} N_{b} U_{11}^{2} $

(2c)

$ C_{13}=\lambda_{b}-\frac{e}{\mu_{b}} \lambda_{b}\left(\lambda_{b}+2 \mu_{b}\right) U_{33}+\frac{e^{2}}{15} \lambda_{b} q U_{33}^{2} $

(2d)

$ q=15\left(\frac{\lambda_{b}}{\mu_{b}}\right)^{2}+28 \frac{\lambda_{b}}{\mu_{b}}+28, N_{b}=\frac{2 \mu_{b}\left(3 \lambda_{b}+8 \mu_{b}\right)}{\lambda_{b}+\mu_{b}} $

(2e)

其中, C_ij代表等效介质刚度系数; λ_b, μ_b为各向同性背景拉梅参数; e=υ_ca³为裂缝密度且υ_c为单位体积岩石中裂缝的条数, a为硬币型裂缝的半长轴长度。U_ii的物理意义是硬币型水平裂缝表面受到i方向的单位应力时, 在i方向产生的位移不连续性, 其值由裂缝自身的参数以及内部充填物的性质决定。由公式(2)可知, U_ii真正独立的分量为U₁₁和U₂₂, 可表达为^[59]U₁₁=(16/3)[(λ_b+2μ_b)/(3λ_b+4μ_b)](1+M)^-1和U₃₃=(4/3)[(λ_b+2μ_b)/(λ_b+μ_b)](1+K)^-1。当裂缝中充填空气时, M=0, K=0;当裂缝中充填无粘性液体时, M=0, K→∝⇒U₃₃=0;当裂缝中充填弱介质或粘性流体时, M=[4μ′(λ_b+2μ_b)]/[πα_ratioμ_b(λ_b+μ_b)], K={[k′+(4/3)μ′)](λ_b+2μ_b)}/[πα_ratioμ_b(λ_b+μ_b)], 其中α_ratio代表裂隙横纵比, k′和μ′分别为充填物的体积模量和剪切模量。

Hudson理论有着较为苛刻的假设条件: ①裂缝为硬币型, 纵横比较小; ②裂缝单组、定向排列; ③裂缝稀疏分布于各向同性背景中且裂缝之间孤立, 即忽略裂缝之间的弹性互动及流体流动; ④裂缝尺度远小于地震波长。在此基础上, HUDSON等^{[60-61, 82]}继续发展该模型, 不断突破假设条件的限制: 针对假设②, 提出了各向同性背景中发育多组定向裂缝系统的等效模型^[60]; 针对假设③, 推导了饱和流体的连通裂隙-裂隙和孔隙-裂隙介质模型^[82]和各向异性背景中包含定向裂缝模型^[61], 后者利用各向同性介质包含虚拟定向裂缝体来模拟各向异性背景, 这是基于弱各向异性背景假设, 本质上并没有考虑真实裂缝系统与各向异性背景之间的作用, 所以对于裂缝型页岩, 由于强各向异性页岩背景, 该模型很难适用。另外, HUDSON^[83]还研究了不饱和流体填充时, 孤立裂缝介质的弹性性质; LIU等^[84]进一步考虑裂缝表面形态, 提出了孤立滑移面裂隙模型、非连续接触裂隙模型和液体薄层裂隙模型; TOD等^[85]针对裂隙面与裂缝面不共面但裂隙中心位于裂缝表面的情况, 提出了裂缝法向与裂隙法向不平行的裂缝模型。

1.1.2 Schoenberg线性滑动理论

SCHOENBERG忽略了裂缝的形状和微观结构, 将其看作一个非常薄的滑动面且表面位移不连续, 这种位移不连续性与穿过裂缝的应力呈线性关系。滑动面与背景介质按照体积比例构成层状介质, 并假设滑动面的刚度模量与其厚度为等价无穷小, 然后利用平滑方法实现裂缝介质等效模量的求取^{[63, 65-66, 86-87]}, 形成Schoenberg线性滑动理论。由该理论发展的模型是一种适用于具有线性连续边界、充满弱强度(充填物模量小)介质的平行层模型。按照Schoenberg线性滑动理论, 含单组平行定向裂缝岩石(如图 2a所示)的有效柔度张量可写为^[67]:

$ \begin{gathered} S_{i j k l}=S_{b_{i j k l}}+\frac{1}{4}\left(Z_{i k} n_{l} n_{j}+Z_{j k} n_{l} n_{i}+\right. \\ \left.Z_{i k} n_{k} n_{j}+Z_{i l} n_{k} n_{i}\right) \end{gathered} $

(3)

式中: 下标b代表背景介质; S_ijkl表示岩石等效柔度张量的第ijkl个元素; Z_ik是裂缝柔度矩阵的第ik个元素; n_i是裂缝法向矢量的第i个分量。介质的柔度张量即为其刚度张量的逆S=C^-1。对于包含多组定向裂缝系统的岩石, 如图 2b所示, 其等效柔度张量可进一步表示为^[67]:

$ S_{i j k l}=S_{b_{i j k l}}+\sum\limits_{m} S_{f_{i j k l}}^{(m)} $

(4)

式中: S_{f ijkl}^(m)=(1/4)(Z_ik^(m)n_l^(m)n_j^(m)+Z_jk^(m)n_l^(m)n_i^(m)+Z_ik^(m)n_k^(m)n_j^(m)+Z_il^(m)n_k^(m)n_i^(m))为第m组定向裂缝柔度张量。

假设一组稀疏、旋转不变的垂直裂缝镶嵌在各向同性背景岩石中, 可用法向、切向柔度Z_N, Z_T对裂缝进行完全表征, 当裂缝法向与x轴一致时, 裂缝岩石等效刚度矩阵C_HTI可用Voigt缩写形式表达为^[67]:

$ \boldsymbol{C}_{\mathrm{HTI}}=\left[\begin{array}{cccccc} M_{b}\left(1-\delta_{\mathrm{N}}\right) & \lambda_{b}\left(1-\delta_{\mathrm{N}}\right) & \lambda_{b}\left(1-\delta_{\mathrm{N}}\right) & 0 & 0 \\ \lambda_{b}\left(1-\delta_{\mathrm{N}}\right) & M_{b}\left(1-\chi_{b}^{2} \delta_{\mathrm{N}}\right) & \lambda_{b}\left(1-\chi_{b} \delta_{\mathrm{N}}\right) & 0 & 0 & 0 \\ \lambda_{b}\left(1-\delta_{\mathrm{N}}\right) & \lambda_{b}\left(1-\chi_{b} \delta_{\mathrm{N}}\right) & M_{b}\left(1-\chi_{b}^{2} \delta_{\mathrm{N}}\right) & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu_{b} & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu_{b}\left(1-\delta_{\mathrm{T}}\right) & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu_{b}\left(1-\delta_{\mathrm{T}}\right) \end{array}\right] $

(5)

式中: M_b, λ_b和μ_b分别为各向同性背景岩石的纵波模量、第一拉梅参数和剪切模量; χ_b=1-2γ_b, γ_b=μ_b/M_b代表背景介质横纵波速度比的平方; δ_N=Z_NM_b/(1+Z_NM_b)和δ_T=Z_Tμ_b/(1+Z_Tμ_b)分别为裂缝法向和切向弱度参数。

综上, Schoenberg线性滑动理论的基本假设可以总结为: ①长波长限制——长波长限制是指裂缝的大小远小于地震波的波长, 相比裂缝间距, 裂缝分布细节及张开度可以忽略, 也就是说, Schoenberg模型将裂缝看成满足线性光滑边界条件的无限薄且非常松软的地层或平面, 忽略裂缝的形状和微结构; ②对于有多套裂缝存在的情况, 忽略裂缝间的相互作用; ③Schoenberg模型中假设裂缝面为位移间断, 但是界面两侧应力保持连续。

实验证明Schoenberg理论中的裂缝和Hudson理论中的裂隙在地震响应上是分不开的^[59], (1)式和(5)式是等效的, 可得^[65]:

$ \delta_{\mathrm{T}}=\frac{16}{3\left[3-2 \gamma_{b}+\frac{4 \mu^{\prime}}{{\rm{ \mathsf{ π} }} \alpha_{\text {ratio }} \mu_{b}}\right]} e^{4} $

(6a)

$ \delta_{\mathrm{N}}=\frac{4}{3 \gamma_{b}\left\{1-\gamma_{b}+\frac{{\left[\kappa^{\prime}+\left(\frac{4}{3}\right) \mu^{\prime}\right]}}{{\rm{ \mathsf{ π} }} \alpha_{\text {ratio }} \mu_{b}}\right\}} e $

(6b)

式中: α_ratio=c/a表示裂缝横纵比。设置μ′=0, 可研究饱和孤立裂缝岩石弱度参数特征; 设置μ′=κ′=0, 可研究干裂缝岩石弱度参数特征。

1.2 衰减各向异性地震岩石物理理论

对于包含饱和、相互连通的定向裂缝及背景孔隙的地下岩石, 当地震波到达时, 会导致孔缝中的流体产生振荡扩散^[88], 这种波致流将诱导地震波衰减各向异性, 发生频散现象, 表现出动力学特征。所以考虑这种波致流衰减的岩石物理模型在本文中称为衰减各向异性地震岩石物理理论, 该类理论已经被多位学者研究, 并给出了各种情况下的岩石弹性模量复值表达式。本节将主要综述修正的Shoenberg线性滑动模型及相关理论和Chapman多尺度裂隙模型。

1.2.1 修正的Schoenberg线性滑动模型

HUDSON根据地下实际情况, 提出了3种多孔裂缝岩石的波致流机制, 如图 3所示。图 3a表示裂缝与裂缝之间连通性引起的流体振荡扩散原理; 图 3b表示裂缝与背景等径孔隙之间的波致流扩散机制; 图 3c表示波场导致部分饱和裂缝内部产生压力梯度, 促使液体振荡扩散占据原来气体的位置, 从而产生衰减的原理。针对这3种波致流机制, HUDSON等^[82-83]展开了详细讨论, POINTER等^[71]在该研究基础上进行总结和发展, 给出不同机制下、与频率相关的$\widetilde{U}_{i i}$(i=1或3)解析表达, 其中字母上方的波浪线表示复值。在此基础上, CHICHININA等^{[72, 75]}假设: ①岩石背景本征衰减远远小于线性滑动裂缝面处诱导的各向异性衰减; ②背景本征衰减是常值, 与地震波传播方向无关, 对衰减各向异性不产生影响, 从而忽略背景本征衰减, 简化各向异性衰减的研究过程; ③背景介质弹性模量为实数。然后通过引入复弱度参数, 提出考虑波致流衰减的修正Shoenberg理论的表达式如下:

$ \widetilde{\boldsymbol{C}}_{\mathrm{VTI}}=\left[\begin{array}{cccccc} M_{b}\left(1-\chi_{b}^{2} \tilde{\delta}_{\mathrm{N}}\right) & \lambda_{b}\left(1-\chi_{b} \tilde{\delta}_{\mathrm{N}}\right) & \lambda_{b}\left(1-\tilde{\delta}_{\mathrm{N}}\right) & 0 & 0 & 0 \\ \lambda_{b}\left(1-\chi_{b} \tilde{\delta}_{\mathrm{N}}\right) & M_{b}\left(1-\chi_{b}^{2} \tilde{\delta}_{\mathrm{N}}\right) & \lambda_{b}\left(1-\tilde{\delta}_{\mathrm{N}}\right) & 0 & 0 & 0 \\ \lambda_{b}\left(1-\tilde{\delta}_{\mathrm{N}}\right) & \lambda_{b}\left(1-\tilde{\delta}_{\mathrm{N}}\right) & M_{b}\left(1-\tilde{\delta}_{\mathrm{N}}\right) & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu_{b}\left(1-\tilde{\delta}_{\mathrm{T}}\right) & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu_{b}\left(1-\tilde{\delta}_{\mathrm{T}}\right) & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu_{b} \end{array}\right] $

(7a)

$ \tilde{\delta}_{\mathrm{N}} \equiv \delta_{\mathrm{N}}-\mathrm{i} \delta_{\mathrm{N}}^{I}=e \frac{2 \widetilde{U}_{33}}{1-\chi_{b}} $

(7b)

$ \tilde{\delta}_{\mathrm{T}} \equiv \delta_{\mathrm{T}}-\mathrm{i} \delta_{\mathrm{T}}^{I}=e \widetilde{U}_{11} $

(7c)

式中: χ_b=1-2γ_b=λ_b/M_b; $\widetilde{U}_{i i}$(ω)表示频率相关的裂隙间断面的位移不连续性, 对于不同的波致流机制, 其表达式均可表达为:

$ \widetilde{U}_{11}=\frac{16}{3} \cdot \frac{\lambda_{b}+2 \mu_{b}}{3 \lambda_{b}+4 \mu_{b}} /(1+\widetilde{M}) $

(8a)

$ \widetilde{U}_{33}=\frac{4}{3} \cdot \frac{\lambda_{b}+2 \mu_{b}}{\lambda_{b}+\mu_{b}} /(1+\widetilde{K}) $

(8b)

但公式(8)中的$\widetilde{M}$和$\widetilde{K}$依赖于波致流机制。

① 对于机制1, 如图 3a所示:

$ \widetilde{M}=\frac{4 a}{{\rm{ \mathsf{ π} }} c}\left(\frac{{\rm i} \omega \eta_{f}}{\mu_{b}}\right)\left(\frac{\lambda_{b}+2 \mu_{b}}{3 \lambda_{b}+4 \mu_{b}}\right) $

(9a)

$ \widetilde{K}=\left[\frac{a}{{\rm{ \mathsf{ π} }} c} \cdot \frac{\kappa_{f}}{\mu_{b}}\left(\frac{\lambda_{b}+2 \mu_{b}}{\lambda_{b}+\mu_{b}}\right)\right]\left[1-\frac{3 \mathrm{i} \kappa_{f} k^{2} K_{r}}{4 {\rm{ \mathsf{ π} }} \upsilon a^{2} c \omega \eta_{f}}\right]^{-1} $

(9b)

式中: η_f和κ_f分别表示裂缝流体的粘滞系数和体积模量; K_r代表岩石渗透率; ω和k表示角频率和波数。

② 对于机制2, 如图 3b所示:

$ \widetilde{M}=\frac{4 a}{{\rm{ \mathsf{ π} }} c}\left(\frac{\mathrm{i} \omega \eta_{f}}{\mu_{b}}\right)\left(\frac{\lambda_{b}+2 \mu_{b}}{3 \lambda_{b}+4 \mu_{b}}\right) $

(10a)

$ \widetilde{K}=\left[\frac{a}{{\rm{ \mathsf{ π} }} c} \frac{\kappa_{f}}{\mu_{b}}\left(\frac{\lambda_{b}+2 \mu_{b}}{\lambda_{b}+\mu_{b}}\right)\right]\left[1+\frac{3(1-\mathrm{i}) J}{2 c}\right]^{-1} $

(10b)

式中: J²=φ_mκ_fK_m/2ωη_f(J>0);φ_m和K_m分别为无裂缝背景基质的孔隙度和渗透率。

③ 对于机制3, 如图 3c所示:

$ \widetilde{M}=\frac{4 a}{{\rm{ \mathsf{ π} }} c} \frac{\mathrm{i} \omega}{\mu_{b}}\left[q_{l} \eta_{l}+\left(1-q_{l}\right) \eta_{g}\right]\left(\frac{\lambda_{b}+2 \mu_{b}}{3 \lambda_{b}+4 \mu_{b}}\right) $

(11a)

式中: q_m和η_m分别代表部分饱和裂缝m相流体的饱和度和粘滞系数, m=l或g, 表示液态相或气态相, q_l+q_g=1。此时, 公式(8b)可以改写为:

$ \widetilde{U}_{33}=\frac{4}{3}\left(\frac{\lambda_{b}+2 \mu_{b}}{\lambda_{b}+\mu_{b}}\right)\left[1-{\rm i} K_{2} /\left(1+K_{1}\right)\right] /\left(1+K_{1}\right) $

(11b)

其中,

$ K_{1}=\frac{a}{{\rm{ \mathsf{ π} }} c \mu_{b}}\left(\frac{\lambda_{b}+2 \mu_{b}}{\lambda_{b}+\mu_{b}}\right)\left(\frac{q_{l}}{\kappa_{l}}+\frac{1-q_{l}}{\kappa_{g}}\right)^{-1} $

(11c)

$ \begin{gathered} K_{2}=\frac{\omega}{{\rm{ \mathsf{ π} }} \mu_{b}}\left(\frac{a}{c}\right)^{3}\left(\frac{\lambda_{b}+2 \mu_{b}}{\lambda_{b}+\mu_{b}}\right)\left(\frac{1}{\kappa_{l}}-\frac{1}{\kappa_{g}}\right)^{2}\left(\frac{q_{l}}{\kappa_{l}}+\right. \\ \left.\frac{1-q_{l}}{\kappa_{g}}\right)^{-2} \times\left\{\eta_{l} F_{l}\left(q_{l}\right)+\eta_{g} F_{g}\left(1-q_{l}\right)\right\} \end{gathered} $

(11d)

并且:

$ F_{l}(q)=F_{1}(q) \approx A_{1}(1-q)[1+\cos ({\rm{ \mathsf{ π} }}-{\rm{ \mathsf{ π} }} q)] $

(11e)

$ F_{g}(q)=F_{2}(q) \approx A_{2}(1-q)[1+\cos ({\rm{ \mathsf{ π} }}-{\rm{ \mathsf{ π} }} q)] $

(11f)

当液态相占据裂缝中心位置时, A₁=0.053, A₂=0.058;当液态相分布在裂缝边缘时, F₁(q)=F₂(q), F_g(q)=F_l(q)。

1.2.2 Chapman多尺度裂缝模型

不同尺度裂缝的定向排列均能诱导地层呈现各向异性特征, 但只有足够大尺度的裂缝才对油气的运移、成藏及生产起主要作用^[89-90]。常规等效介质理论采用裂缝密度、纵横比等参数来考虑裂缝对等效弹性模量的影响, 并未研究裂缝尺度与弹性模量的关系。多尺度非均匀性的引入诱导地震反射/透射的频散现象^[91-92], 所以多尺度定向裂缝的引入也会产生各向异性的频变特征, LIU等^[93]和MAULTZSCH等^[94]在实际数据处理中已经发现了各向异性与频率相关的现象。为了分析多尺度裂缝的地震各向异性响应特征, 实现多尺度裂缝的识别, CHAPMAN^[77]提出频变多尺度裂缝岩石物理模型, 接下来对该模型进行简要阐述。

如图 4所示, Chapman多尺度裂缝模型^[76-77]假设岩石中含球状等径孔隙、随机排列的硬币型微裂隙以及定向排列的硬币型裂缝, 其中基质等径孔隙和微裂隙与岩石颗粒尺度相当, 而定向排列的裂缝尺度远大于颗粒尺度, 但小于地震波长尺度。同时, 假设微裂隙与微裂隙间、背景等径孔隙与微裂隙间相互连通构成单元c₁, 定向排列的裂缝与多个等径孔隙或微裂隙连通构成单元c₂, 但每个等径孔隙或微裂隙最多与一条定向排列的裂缝连通, 且定向裂缝之间不连通。当地震波到达含孔隙裂隙岩石时, 认为压力梯度将诱导两种不同尺度流体扩散效应, 即定向排列的裂缝与背景等径孔隙、微裂隙之间的达西流效应, 及微裂隙之间或与背景等径孔隙之间的喷射流效应。因此Chapman模型的等效刚度参数C_chapman可表示为:

$ \boldsymbol{C}_{\text {chapman }}=\boldsymbol{C}_{0}^{\mathrm{ch}}-\varphi_{p} \boldsymbol{C}_{1}^{\mathrm{ch}}-e_{m c} \boldsymbol{C}_{2}^{\mathrm{ch}}-e_{c} \boldsymbol{C}_{3}^{\mathrm{ch}} $

(12)

式中: C₀^ch是背景基质岩石的弹性刚度张量; C₁^ch、C₂^ch及C₃^ch分别为背景等径孔隙、微裂隙及定向裂缝的弹性刚度贡献张量, 是流体体积模量、裂缝参数、频率及喷射流弛豫时间参数的函数; φ_p是基质等径孔隙度; e_c是定向裂缝密度; e_mc是微裂隙密度。在低频假设下, Chapman多尺度裂隙模型与各向异性Gassmann方程一致; 在高频条件下, 与Hudson孤立裂隙模型相一致。需要注意的是, 公式(12)是在交互能量法的基础上推导而来, 所以只有在稀疏孔裂隙的情况下, 公式(12)才成立。

中观尺度定向裂缝和微观尺度孔裂隙的存在引入两个尺度的弛豫时间参数τ_c和τ_m, 且两者存在近似函数关系:

$ \tau_{c}=\frac{a_{c}}{a_{g}} \tau_{m} $

(13)

式中: a_c和a_g分别是定向排列裂缝与颗粒尺寸。公式(13)物理含义为: 随着裂缝半径增大, 裂缝表面积与体积的比值变小, 达到流体压力平衡时其单位表面积将有更多的流体通过, 所需时间也更长; 而与定向裂缝相关的弛豫时间τ_c将导致地震频率范围内出现速度的频散与衰减, 也即定向裂缝诱导的各向异性具有频率相关性。由于实际应用中, 无缝多孔背景的模量较易获得, 为此, CHAPMAN等^[96]将公式(12)进行扩展, 得:

$ \begin{gathered} \boldsymbol{C}_{\text {chapman }}(f)=\boldsymbol{C}_{0}^{\mathrm{ch}}(\varGamma, Z)-\varphi_{p} \boldsymbol{C}_{1}^{\mathrm{ch}}\left(\lambda_{g}, \mu_{g}, f\right) \\ e_{m c} \boldsymbol{C}_{2}^{\mathrm{ch}}\left(\lambda_{g}, \mu_{g}, f\right)-e_{c} \boldsymbol{C}_{3}^{\mathrm{ch}}\left(\lambda_{g}, \mu_{g}, f\right) \end{gathered} $

(14)

式中: f是频率; 固体颗粒拉梅常数分别为λ_g=Γ-Φ_{mc, p}(λ_g, μ_g, f_ω)与μ_g=Z-Φ_{mc, p}(λ_g, μ_g, f_ω), 其中, Γ和Z是频率f_ω下经背景等径孔隙与微裂隙校正后的拉梅常数, Φ_{mc, p}是校正项。对于大多数的裂缝型储层, 微裂隙孔隙度要远远小于等径孔隙度。当勘探频率低于微观喷射流频率时, 可忽略公式(14)中的微裂隙项。

CHAPMAN^[97]在单组定向多尺度裂缝频变弹性特征的研究基础上, 提出了多组定向裂缝频率依赖的岩石物理模型。主要解决思路为: 首先求出由于等径孔隙存在而产生的弹性参数修正量, 然后求出本征坐标系下, 每组裂缝单独存在时产生的弹性修正量, 接着通过坐标旋转将裂缝引起的修正量变换至观测坐标系下, 最后将背景基质张量、等径孔隙修正量、多组裂缝修正量求和, 便可得到最终的弹性张量, 可表达为:

$ \boldsymbol{C}_{i j k l}(f)=\boldsymbol{C}_{i j k l}^{0}-\boldsymbol{C}_{i j k l}^{p}(f)-\sum\limits_{n=1}^{l} \boldsymbol{B}_{i p}^{n} \boldsymbol{B}_{j q}^{n} \boldsymbol{B}_{k r}^{n} \boldsymbol{B}_{l s}^{n} c_{p q r s}^{n}(f) $

(15)

式中: C_ijkl⁰为背景基质弹性模量; C_ijkl^p(f)为与等径孔隙相关的模量修正量; B_ip为坐标变换矩阵; c_pqrsⁿ(f)为第n组定向裂缝在本征坐标系下的弹性模量修正量; C_ijkl(f)为含多组定向裂缝的多孔介质的频率依赖刚度张量。

1.3 极限频率下考虑孔缝连通性的地震岩石物理理论

上节表明: 地球岩石受到外界波场扰动时, 岩石内流体在压力梯度的作用下会发生振荡扩散, 且这种扩散特征与波场频率有关。地震频率相对较低, 一般认为在地震波扰动下, 岩石内流体在连通孔缝间有充足时间流动扩散以达到压力均衡, 表现为极限频率下的动力学特征。所以本节将以极限频率下考虑孔缝连通性的地震岩石物理理论为主题, 重点介绍Thomsen含等径孔隙裂隙模型和基于广义Gassmann方程的饱和裂缝岩石模型。

1.3.1 Thomsen含等径孔隙裂隙模型

Thomsen将裂隙看作一个横纵比较小的椭球体, 然后假设一组水平定向裂隙稀疏地分布在各向同性多孔岩石中, 裂隙与背景等径孔隙之间是连通的, 在压力梯度作用下, 流体发生振荡扩散, 如图 5所示。在先前弱各向异性参数的研究基础上^[98], THOMSEN^[78]提出了考虑流体压力均衡的孔隙裂隙介质模型:

$ \varepsilon_{\mathrm{sat}}^{\mathrm{VTI}}=\frac{8}{3}\left(1-\frac{\kappa_{f}}{K_{g}}\right) D_{c i} \frac{\left[1-\left(\sigma_{b}^{\mathrm{dry}}\right)^{2}\right] E_{b}^{\mathrm{sat}}}{\left[1-\left(\sigma_{b}^{\mathrm{sat}}\right)^{2}\right] E_{b}^{\mathrm{dry}}} \mathrm{e} $

(16a)

$ \gamma_{\mathrm{sat}}^{\mathrm{VTI}}=\frac{8}{3} \cdot\left(\frac{1-\sigma_{b}^{\mathrm{dry}}}{2-\sigma_{b}^{\mathrm{dry}}}\right) e $

(16b)

$ \delta_{\mathrm{sat}}^{\mathrm{VTI}}=2\left(1-\sigma_{b}^{\mathrm{sat}}\right) \varepsilon_{\mathrm{sat}}^{\mathrm{VTI}}-2 \frac{1-2 \sigma_{b}^{\mathrm{sat}}}{1-\sigma_{b}^{\mathrm{sat}}} \gamma_{\mathrm{sat}}^{\mathrm{VTI}} $

(16c)

式中: ε_sat^VTI, γ_sat^VTI, δ_sat^VTI分别为Thomsen各向异性参数; K_g表示固体颗粒(基质)的有效体积模量; E和σ分别表示杨氏模量和泊松比; 下标b、上标dry和sat分别表示各向同性背景及干孔和流体饱和孔条件。D_ci为流体影响因素, 在不同频段, 其形式也不同。对于地震勘探, 我们常用其低频形式:

$ \begin{gathered} D_{c i}\left(f_{\text {low }}\right)=\left\{1-\frac{\kappa_{f}}{K_{g}}+\frac{\kappa_{f}}{K_{b}^{\mathrm{dry}}\left(\varphi_{c}+\varphi_{p}\right)} \cdot\right. \\ \left.\left\{\left(1-\frac{K_{b}^{\mathrm{dry}}}{K_{g}}\right)+\frac{16}{9} \cdot\left[\frac{1-\left(\sigma_{b}^{\mathrm{dry}}\right)^{2}}{1-2 \sigma_{b}^{\mathrm{dry}}}\right] e\right\}\right\}^{-1} \end{gathered} $

(16d)

此时裂缝与背景孔隙之间有充足的时间通过流体流动达到局部压力均衡; 在中高频条件下, 背景孔隙与裂缝之间没有足够的时间达到流体压力均衡, 则有:

$ \begin{gathered} D_{c i}\left(f_{\mathrm{m}-\mathrm{h}}\right)=\left\{1-\frac{\kappa_{f}}{K_{g}}+\frac{\kappa_{f}}{K_{b}^{\mathrm{dry}}} \cdot\left\{\frac{16}{9} \cdot\left[\frac{1-\left(\sigma_{b}^{\mathrm{dry}}\right)^{2}}{1-2 \sigma_{b}^{\mathrm{dry}}}\right] \cdot\right.\right. \\ \left.\left.\frac{e}{\varphi_{c}}\left(\frac{1-\frac{\kappa_{f}}{K_{g}}}{1-\frac{\kappa_{f}}{K_{b}^{\mathrm{sat}}}}\right)\right\}\right\}^{-1} \end{gathered} $

(16e)

式中: f_low和f_m-h分别代表低频和中高频条件; φ_c为裂隙孔隙度; φ_p为背景等径孔隙度。

1.3.2 基于广义Gassmann方程的饱和裂缝岩石模型

GASSMANN^[99]提出了等效的各向异性岩石流体替换方程^[100]:

$ C_{i j}^{\mathrm{sat}}=C_{i j}^{\mathrm{dry}}+\beta_{i} \beta_{j} K_{p}^{\mathrm{ani}} $

(17a)

$ i, j=1,2, \cdots, 6 $

式中: C_ij^sat和C_ij^dry分别为各向异性饱和裂缝岩石和干岩石的刚度系数; β_i表示各向异性Biot系数, 即:

$ \begin{gathered} \beta_{m}=\left(\begin{array}{cc} 1-\frac{\sum\limits_{n=1}^{3} C_{m n}^{\mathrm{dry}}}{3 K_{g}} \end{array}\right) \vartheta_{m} \\ \vartheta_{m}= \begin{cases}1 & m=1,2,3 \\ 0 & m=4,5,6\end{cases} \end{gathered} $

(17b)

K_p^ani表示各向异性Gassmann孔隙空间模量, 可表达为:

$ K_{p}^{\mathrm{ani}}=\frac{K_{g}}{\left(1-\frac{K_{\mathrm{dry}}^{\mathrm{ani}}}{K_{g}}\right)-\varphi\left(1-\frac{K_{g}}{\kappa_{f}}\right)} $

(17c)

式中: φ=φ_c+φ_p为裂隙岩石总孔隙度; K_dry^ani=$\sum\limits_{i=1}^{3} \sum\limits_{j=1}^{3}$C_ij^dry/9被定义为各向异性干岩石抗压缩模量。BROWN等^[101]针对非均质各向异性介质, 以柔度矩阵的形式给出了公式(17)的等效表达式。

对于垂直定向裂缝诱导的水平横向各向同性介质(horizontally transverse isotropy, HTI), 联合Schoenberg线性滑动模型和广义Gassmann方程, 公式(17b)和公式(17c)可表达为^{[79, 102]}:

$ \beta_{1}=\beta_{0}+\left(1-\beta_{0}\right) \delta_{\mathrm{N}}^{\mathrm{dry}} $

(17d)

$ \beta_{2}=\beta_{3}=\beta_{0}+\left(1-\beta_{0}\right) \chi_{b}^{\mathrm{dry}} \delta_{\mathrm{N}}^{\mathrm{dry}} $

(17e)

$ K_{p}^{\text {ani }}=\frac{K_{g}}{1-\varphi-{\rm{ \mathsf{ π} }}\left(1-\beta_{0}\right)\left(1-\frac{K_{b}^{\text {dry }}}{M_{b}^{\text {dry }}} \delta_{\mathrm{N}}^{\text {dry }}\right)+\varphi \frac{K_{g}}{\kappa_{f}}} $

(17f)

式中: β₀=1-$\frac{K_{b}^{\mathrm{dry}}}{K_{g}}$为各向同性Biot系数; K_b^dry代表各向同性背景干岩石有效体积模量; δ_N^dry为各向异性干岩石弱度参数。结合公式(17a)至公式(17f), 含孔隙裂缝饱和岩石的刚度参数可表达为^{[79, 102]}:

$ \begin{aligned} C_{11}^{\mathrm{sat}}=&M_{b}^{\mathrm{dry}}\left(1-\delta_{\mathrm{N}}^{\mathrm{dry}}\right)+\\ & \frac{\left[\beta_{0}+\left(1-\beta_{0}\right) \delta_{\mathrm{N}}^{\mathrm{dry}}\right]^{2}}{\frac{\varphi}{\kappa_{f}}+\frac{1-\varphi}{K_{g}}-\frac{1-\beta_{0}}{K_{g}}+\frac{\left(1-\beta_{0}\right)^{2}}{M_{b}^{\mathrm{dry}}} \delta_{\mathrm{N}}^{\mathrm{dry}}} \end{aligned} $

(18a)

$ \begin{aligned} C_{13}^{\mathrm{sat}}=&\lambda_{b}^{\mathrm{dry}}\left(1-\delta_{\mathrm{N}}^{\mathrm{dry}}\right)+\\ & \frac{\left[\beta_{0}+\left(1-\beta_{0}\right) \delta_{\mathrm{N}}^{\mathrm{dry}}\right]\left[\beta_{0}+\left(1-\beta_{0}\right) \chi_{b}^{\mathrm{dry}} \delta_{\mathrm{N}}^{\mathrm{dry}}\right]}{\frac{\varphi}{\kappa_{f}}+\frac{1-\varphi}{K_{g}}-\frac{1-\beta_{0}}{K_{g}}+\frac{\left(1-\beta_{0}\right)^{2}}{M_{b}^{\mathrm{dry}}} \delta_{\mathrm{N}}^{\mathrm{dry}}} \end{aligned} $

(18b)

$ \begin{aligned} C_{23}^{\mathrm{sat}}=&\lambda_{b}^{\text {dry }}\left(1-\chi_{b}^{\text {dry }} \delta_{\mathrm{N}}^{\mathrm{dry}}\right)+\\ &\frac{\left[\beta_{0}+\left(1-\beta_{0}\right) \chi_{b}^{\mathrm{dry}} \delta_{\mathrm{N}}^{\mathrm{dry}}\right]^{2}}{\frac{\varphi}{\kappa_{f}}+\frac{1-\varphi}{K_{g}}-\frac{1-\beta_{0}}{K_{g}}+\frac{\left(1-\beta_{0}\right)^{2}}{M_{b}^{\mathrm{dry}}} \delta_{\mathrm{N}}^{\mathrm{dry}}} \end{aligned} $

(18c)

$ \begin{aligned} C_{33}^{\mathrm{sat}}=& M_{b}^{\mathrm{dry}}\left[1-\left(\chi_{b}^{\mathrm{dry}}\right)^{2} \delta_{\mathrm{N}}^{\mathrm{dry}}\right]+\\ & \frac{\left[\beta_{0}+\left(1-\beta_{0}\right) \chi_{b}^{\mathrm{dry}} \delta_{\mathrm{N}}^{\mathrm{dry}}\right]^{2}}{\frac{\varphi}{\kappa_{f}}+\frac{1-\varphi}{K_{g}}-\frac{1-\beta_{0}}{K_{g}}+\frac{\left(1-\beta_{0}\right)^{2}}{M_{b}^{\mathrm{dry}}} \delta_{\mathrm{N}}^{\mathrm{dry}}} \end{aligned} $

(18d)

$ C_{44}^{\mathrm{sat}}=\mu_{b} $

(18e)

$ C_{55}^{\mathrm{sat}}=\mu_{b}\left(1-\delta_{\mathrm{T}}\right) $

(18f)

$ C_{22}^{\mathrm{sat}}=C_{33}^{\mathrm{sat}} $

(18g)

$ C_{55}^{\mathrm{sat}}=C_{66}^{\mathrm{sat}} $

(18h)

$ C_{12}^{\mathrm{sat}}=C_{13}^{\mathrm{sat}} =C_{21}^{\mathrm{sat}}=C_{31}^{\mathrm{sat}} $

(18i)

$ C_{23}^{\mathrm{sat}}=C_{32}^{\mathrm{sat}} $

(18j)

其它刚度系数均为零。由于调用了广义Gassmann方程, 因此公式(18)是低频假设下, 含连通孔隙裂缝岩石的等效刚度系数, 即要求在地震波的扰动下, 背景等径孔隙与定向裂缝中的流体有充足的时间达到压力均衡。GUREVICH^[79]研究表明, 由于该模型考虑了背景孔隙和裂缝之间的流体流动, 所以不再满足HTI型岩石总柔度是各向同性背景柔度与裂缝柔度的叠加。CIZ等^[103]通过假设在加载围压和孔隙空间应力时岩石骨架和孔隙充填材料变形较小, 将各向异性流体替换方程推广为各向异性固体替换方程, 以实现孔隙空间充填固体或固-液混合物的各向异性储层刚度参数的估算, 可用柔度张量S_ijkl表达为:

$ S_{i j k l}^{\mathrm{sat}}=S_{i j k l}^{\mathrm{dry}}-\frac{\left(S_{i j k l}^{\mathrm{dry}}-S_{i j k l}^{0}\right)\left(S_{m n p q}^{\mathrm{dry}}-S_{m n p q}^{0}\right)}{\left[\left(S^{\mathrm{dry}}-S^{0}\right)+\varphi\left(S^{i f}-S^{\varphi}\right)\right]_{m n p q}} $

(19)

式中: 上标sat, dry和0分别代表与饱和岩石、无充填多孔骨架和背景基质相关的量。上标if和φ分别代表与孔隙充填物体平均应变有关的量和与骨架孔隙空间有关的量, 且CIZ等^[103]也给出了S^if和S^φ的估计方式。

岩石物理模型是利用五维地震数据预测裂缝参数的理论基础。本文根据储层流体在地震波扰动下的响应特征, 从3个角度分别介绍了当前国内外常用的裂缝型储层地震岩石物理理论。每种模型都有其适用的地质环境: ①对于喉道不发育的储层, 当微裂隙/裂缝饱和均匀流体时, 可以采用孤立孔缝Hudson理论^[58-59]、线性滑动模型^{[65, 67-68]}或高频Thomsen含等径孔隙裂隙模型^[78]来模拟裂缝、流体等对岩石属性的影响; 当裂隙/裂缝充填多相流体时, 可以采用部分饱和Hudson模型^[83]等; ②对于喉道充分发育的裂缝型地层, 在地震数据超低频假设下, 可以采用基于广义Gassmann方程的饱和裂缝岩石模型^{[79, 102]}或低频Thomsen含等径孔隙裂隙模型^[78]计算地层弹性参数, 但研究岩石的等效弹性模量在全频带的特征时, 应该采用连通孔-裂隙的Hudson模型^{[82, 85]}等。此时, 我们并没有强调基质孔隙(或硬孔隙)的发育情况, 但可以通过调节每种模型的相关参数来改变基质孔隙的含量以适应更广泛的地质情况。另外, 当上述两类储层发育多组相同尺度裂缝体时, 多组裂隙Hudson模型^[60]、线性滑动模型^[104]等将是较好的选择。但对于基质孔隙-随机裂隙-定向裂缝发育的储层, 地震波经过时诱导多尺度流体流动, 此时可采用Chapman波致流衰减理论等多尺度裂缝模型^{[77, 105]}来研究储层宏观弹性特征。总之, 地下介质非均质性强, 孔隙、裂缝结构多样, 油气水关系复杂, 一种岩石物理理论难以全面刻画地下岩石的弹性特征, 应该针对不同的地质环境, 研究和构建适用的地震岩石物理模型。

2 裂缝型储层地震响应特征分析

研究不同岩石物理模型之间的联系和差异, 分析裂缝型储层微观因素对宏观地震响应(反射系数)的影响, 将为五维地震裂缝参数预测方法及稳定性研究提供理论基础。所以本节在前文介绍的岩石物理理论基础上, 主要讨论Hudson模型、Thomsen含等径孔隙裂隙模型和基于广义Gassmann方程的饱和裂缝岩石模型(后面将其简称为Gurevich模型)间的联系和差异, 并分析裂缝、流体及背景孔隙度等参数对地震响应特征(反射系数)的影响, 从而为裂缝和流体参数反演策略的制定提供思路。

为了便于分析裂缝诱导各向异性特征, 我们设计一个两层模型, 如图 6所示, 如无特殊说明, 上层介质始终为各向同性致密介质, 下层为发育定向垂直裂缝的储层, 裂缝法向与x轴相同, 上下层介质的矿物组成相同且分布均匀。孔-裂隙空间中饱含均匀流体(水的体积模量2.2GPa, 油的体积模量1.4GPa, 气的体积模量0.133GPa)。表 1展示了两类不同矿物的速度和密度参数, 这两类矿物的区别是横纵波速度比的平方不一样, 矿物A为0.475, 矿物B为0.350。采用RÜGER^[106]给出的反射系数方程研究裂缝密度、充填物和孔隙度等参数对地震响应的影响, 方程表达如下:

$ R=R_{\mathrm{PP}}^{\text {iso }}+R_{\mathrm{PP}}^{\mathrm{ani}} $

(20a)

$ \begin{gathered} R_{\mathrm{PP}}^{\text {iso }}=\frac{1}{2}\left(\frac{\Delta Z}{\bar{Z}}\right) +\frac{1}{2}\left[\frac{\Delta \alpha}{\bar{\alpha}}-\left(\frac{2 \bar{\beta}}{\bar{\alpha}}\right)^{2} \frac{\Delta G}{\bar{G}}\right] \sin ^{2} \theta+\\ \frac{1}{2}\left(\frac{\Delta \alpha}{\bar{\alpha}}\right) \sin ^{2} \theta \tan ^{2} \theta \end{gathered} $

(20b)

$ \begin{gathered} R_{\mathrm{PP}}^{\mathrm{ani}}= \frac{1}{2}\left\{\left[\Delta \delta^{(V)}+2\left(\frac{2 \bar{\beta}}{\bar{\alpha}}\right)^{2} \Delta \gamma\right] \sin ^{2} \theta+\left[\Delta \varepsilon^{(V)}\cdot\right.\right.\\ \left.\left.\cos ^{2} \varphi+\Delta \delta^{(V)} \sin ^{2} \varphi\right] \sin ^{2} \theta \tan ^{2} \theta\right\} \cos ^{2} \varphi \end{gathered} $

(20c)

公式(20)中, G=ρβ²和Z=ρα。R代表反射系数, θ为入射角, φ为方位角。α和β分别为垂直方向的纵、横波速度, ρ为密度, Δδ^(V), Δε^(V)和Δγ为上下两层介质Thomsen各向异性参数差值。接下来, 我们将从背景无等径孔隙和背景包含等径孔隙两种情况展开讨论。

2.1 背景无等径孔隙的岩石(背景孔隙度假设为0)

假设背景介质仅由均匀的矿物A组成, 则上下层介质的差异仅由垂直裂缝和充填的流体而产生。采用Hudson孤立裂隙模型、Thomsen含等径孔隙裂隙模型和Gurevich模型分别计算干岩石和饱含水岩石的刚度模量, 并代入公式(20)可计算反射系数。图 7展示了由Gurevich模型和公式(20)计算的无背景孔隙的两层模型反射系数随入射角和方位角变化的三维图。图 7a和图 7b中裂缝分别无充填(干裂缝)和饱含水, 沿红色箭头方向, 下层裂缝密度依次为0.001, 0.010, 0.025和0.050。从图 7中可以发现, 裂缝引起的方位振幅变化较弱, 且随着充填物模量的增大, 方位振幅逐渐降低。图 8和图 9展示了无背景孔隙的两层模型反射系数随方位角和入射角变化的二维曲线。深粉红色、青蓝色、绿色和深蓝色曲线分别代表利用低频Thomsen理论、中高频Thomsen理论、Gurevich模型和Hudson模型计算得到的反射系数, 图 8和图 9对应的下层裂缝密度分别为0.001和0.050, 裂缝纵横比均为0.005。

对于干裂缝情况, 从图 8a、图 8b和图 9a、图 9b中我们可以发现, Gurevich模型与Hudson模型的结果一致, 但与高、低频Thomsen模型的结果存在差异, 同时, 高频与低频Thomsen模型的结果相同。这是由于Gurevich模型是在Hudson干裂缝模型基础上采用各向异性Gassmann流体替换而推导的, 所以对于干裂缝, 二者是一致的。但此时, Thomsen理论给出的各向异性参数表达是Hudson理论给出的结果关于裂缝密度的线性近似^[107]。因此小裂缝密度情况下, Hudson理论、Gurevich模型和Thomsen理论比较接近, 正如图 8a和图 8b所示; 大裂缝密度情况下, Hudson理论、Gurevich模型和Thomsen理论的差异较大, 这也被图 9a和图 9b证实。另外, 由于干燥岩石不存在流体流动现象, 所以高频与低频Thomsen模型的结果相同。

对于饱含水的岩石, 分析图 8d和图 9d可以发现, 当小裂缝密度时, 基于Hudson理论、Gurevich模型和Thomsen理论的反射系数差异较小; 对于大裂缝密度, Hudson理论和Gurevich模型所得结果的差异较小, 但与Thomsen理论给出的反射系数有较大差异。这是由于在小裂缝密度情况下, 流体的作用大于裂缝的影响, 尽管Gurevich模型与低频Thomsen理论假设激励频率足够低, 但定向裂隙的形状相同, 所以在地震波长范围内孔隙压力均衡, 不会发生裂隙之间的流体流动, 其表现与孤立裂缝情况相似^[71]。因此Gurevich模型、Hudson模型和Thomsen理论在小裂缝密度情况下表现相近。在大裂缝密度时, 裂缝影响大于充填水的作用, 所以Thomsen理论的线性近似导致其与Gurevich模型和Hudson模型间的差异变大。对比图 8c和图 9c可以发现, 无论裂缝密度大小, 在裂缝走向上(观测方位90°), Hudosn模型与Gurevich模型存在差异, 且Gurevich模型与Thomsen理论更接近。这是由于裂缝走向的反射系数由上下层垂直速度和密度的差异决定, 与各向异性参数无关(根据公式(20)可知), 流体充填弱化裂缝的作用, Gurevich模型与Thomsen理论考虑流体作用的方式相同, 所以导致图 8c和图 9c所示的现象。这种现象是理论模型差异而导致, 并非地层属性所造成。

综上, 对于背景孔隙度极其小的地层(例如页岩地层), 地震波的传播不再引起孔-裂隙流体的流动, 和高频极限情况一致, 所以无论喉道结构发育情况, 我们均可采用孤立裂隙模型(Hudson模型和高频Thomsen理论)研究地层的孔弹性特征。且当裂缝密度较小时, 利用Hudson理论、Gurevich模型和Thomsen理论均可获得较好的理论模拟, 但当裂缝密度较大时, Thomsen理论的模拟结果精度较Hudson理论和Gurevich模型的精度低。

在上述分析基础上, 我们采用Hudson模型研究的裂缝密度、充填物类型对地震反射特征的影响。图 10和图 11展示了矿物A和矿物B分别作为基质矿物时的反射系数R(θ, φ)随方位角φ和入射角θ的变化, 其中每种线型代表不同的裂缝密度, 依次为0.01, 0.02, 0.03, 0.05, 裂缝横纵比为0.005, 图 10a、图 10c和图 10e以及图 11a、图 11c和图 11e代表裂缝走向反射系数相对垂直入射反射系数的改变量(即R(θ, 90°)-R(0, 90°))随入射角变化, 图 10b、图 10d和图 10f以及图 11b、图 11d和图 11f代表入射角为30°反射系数相对裂缝法向反射系数的改变量(即R(30°, φ)-R(30°, 0))随方位角变化。图 10a、图 10b、图 11a和图 11b为气充填裂缝, 图 10c、图 10d、图 11c和图 11d为油充填, 图 10e、图 10f、图 11e和图 11f为水充填。图 10和图 11显示地震反射振幅随观测方位角呈现余弦变化规律, 且随着裂缝密度的增加, 由定向裂缝引起的AVO梯度及振幅方位变化幅度均会增加。这是由于裂缝的存在会降低岩石中波的传播速度, 进而增大界面上下岩石的弹性差异, 从公式(20b)可得, AVO梯度将会增加; 公式(20c)表明反射振幅方位变化幅度是各向异性参数的线性函数, 各向异性参数是裂缝密度的线性函数^[108], 则反射振幅方位变化幅度是裂缝密度的近似线性函数, 所以裂缝密度的增加会提高振幅方位变化幅度。很多油气勘探工作者已经利用这种特性开展裂缝密度预测方法研究。

裂缝充填物的模量越小, ε^(V)和δ^(V)的绝对值将会变大, γ与充填物性质无关^[108-109]。比较图 10b和图 10f可以发现, 气充填的裂缝储层反射振幅方位变化幅度(|ΔR|)小于水充填的裂缝储层, 所以各向异性越强并不代表反射振幅方位变化幅度越大。这是由于反射振幅方位变化幅度由ΔR=[Δδ^(V)+8(β/α)²Δγ]sin²θ+Δε^(V)sin²θtan²θ决定(即R(θ, 0)-R(θ, 90°)), 且对于裂缝诱导的各向异性介质, ε^(V)和δ^(V)一般为负值, γ为正值。当水充填裂缝时, ε^(V)和δ^(V)的绝对值非常小, γ起主导作用, 所以ΔR>0, R(θ, 0)>R(θ, 90°), 图 10f证明了这个结论; 随着充填物的体积模量减小, ε^(V)和δ^(V)的绝对值将变大, 则ΔR逐渐变小, 所以图 10d中油充填的裂缝储层反射振幅方位变化幅度小于图 10f中水充填的裂缝储层; 当气充填裂缝时, ε^(V)和δ^(V)起主导作用, 则ΔR<0, R(θ, 0)<R(θ, 90°), 图 10b体现了这一点, 但由于β/α=0.475比较大, 所以|ΔR|比较小, 即矿物A作为基质成分时, 气充填的裂缝储层反射振幅方位变化幅度小于水充填的裂缝储层。图 11与图 10不同的是, 图 11b中气充填的裂缝储层反射振幅方位变化幅度大于图 11f中水充填的情况, 这是由于矿物B的β/α=0.35比较小, 所以气充填时, 反射振幅方位变化幅度|ΔR|将更大。

综上, 对于背景孔隙度非常小或孔喉结构发育较差的地层, 尽管裂缝密度的增大将会导致地震反射振幅方位变化幅度的增大, 但裂缝充填物类型也会影响反射振幅方位变化幅度, 所以直接利用反射振幅方位变化幅度难以准确预测裂缝密度与流体类型。为了准确预测裂缝参数、实现油气识别, 需要研究如何从反射振幅方位变化信息中解耦背景基质、裂缝及油气等信息, 以提高裂缝型储层预测精度。除此, 虽然反射振幅方位变化幅度与充填物类型之间的关系依赖横纵波速度比等参数, 但从图 10和图 11中可以发现充填物的模量越小——从水充填到气充填, 裂缝走向(90°观测方位)反射振幅与裂缝倾向(0观测方位)反射振幅的差(即ΔR)将会越大。这或许为流体识别指明一条好的方向——在获得裂缝方位的前提下, 可利用五维地震数据直接检测裂缝流体类型。

2.2 背景基质发育等径孔隙的岩石

假设上下层基质均由矿物A组成, 下层背景介质包含10%的等径孔隙, 此时如果孔喉发育良好, 在低频激励下, 下层介质流体将在定向裂缝和背景孔隙之间流动以平衡局部压力。我们依旧采用Hudson孤立裂隙模型、Thomsen含等径孔隙裂隙模型和Gurevich模型分别计算饱含气和饱含水岩石的刚度模量, 并代入公式(20)计算反射系数。图 12和图 13展示了反射系数随方位角和入射角变化的曲线, 图中不同线型与图 8中线型代表的意义一致, 且图 12和图 13对应的下层裂缝密度分别为0.001和0.050, 裂缝纵横比均为0.005。

Hudson孤立裂隙模型和高频Thomsen模型均假设地层孔裂隙之间流体无法流动, 低频Thomsen模型和Gurevich模型均考虑了因局部压力差而导致的孔-裂隙间流体流动。比较图 9a和图 13a, 可以发现孤立孔裂隙和连通孔裂隙介质的地震响应存在差异, 即流体流动对地层地震响应的影响不可忽视, 这也证明了地震方法预测储层流体流动性是可能的。图 13a中Gurevich模型和低频Thomsen模型之间存在差异的原因之一是Thomsen模型对裂缝密度的线性近似, 另一部分原因是Thomsen假设裂缝与流体引起的柔度扰动矩阵为对角矩阵, 而Gurevich模型中裂缝与流体引起的柔度扰动矩阵不为对角矩阵。所以, 从图 13a到图 13c, 水充填孔裂隙替换气充填孔裂隙, 减弱了裂缝的影响, 则图 13c中Gurevich模型和低频Thomsen模型之间的差异较图 13a更小。图 12对应的裂缝密度非常小, 流体流动效应也非常微弱, 所以图 12a和图 12c中, 基于这4种模型的地震响应差异非常小。

为了分析流体流动等因素对地震方位响应的影响, 我们绘制地震振幅方位差异(任意观测方位地震振幅与0观测方位的地震振幅的差异)随方位角变化的曲线, 如图 12b, 图 12d和图 13b, 图 13d所示。我们可以发现孤立孔-裂隙介质的方位振幅变化规律(Hudson模型和高频Thomsen理论)与背景无等径孔隙情况(图 10和图 11)相似, 这是由于此时含等径孔隙的背景可作为一个整体, 作用等效于无等径孔隙的背景。除此, 流体流动效应对地震振幅方位响应有着较大的影响。不同于孤立裂缝的情况, 饱含流体连通裂缝情况的裂缝走向(90°观测方位)反射振幅一直大于裂缝倾向(0观测方位)反射振幅。这是因为考虑孔-缝流体流动时, γ不被影响, 但ε^(V)模值增大且达到主导地位, 所以ΔR<0, 即R(θ, 0)<R(θ, 90°)。因此, 此时可以利用反射振幅的极大值点对应的观测方位估计裂缝走向, 预测结果不存在不确定性。而且流体模量越小, 孔-缝连通介质的地震振幅方位变化越剧烈, 该现象可指导喉道比较发达的含孔裂隙地层的流体识别。比较图 12b和图 13b(或比较图 12d和图 13d)也可知, 无论孔-裂隙地层是连通的还是孤立的地层, 裂缝密度越大, 地震振幅方位变化均越明显。

仍然假设上下层基质均由矿物A组成, 下层介质裂缝密度为0.050。首先假设裂缝横纵比为0.005, 采用Gurevich模型和公式(20)绘制不同等径孔隙度下反射振幅方位差异随观测方位变化的曲线, 如图 14所示, 其中图 14a代表孔裂隙充填气, 图 14b代表孔裂隙充填水。然后假设背景孔隙度固定, 采用Gurevich模型和公式(20)绘制不同裂缝纵横比下反射振幅方位差异随观测方位变化的曲线, 如图 15所示, 其中图 15a和图 15b代表背景孔隙度为0.005, 图 15c和图 15d代表背景孔隙度为0.050, 图 15a和图 15c代表孔裂隙充填气, 图 15b和图 15d代表充填水。根据图 14, 我们可以发现, 对于含气储层(图 14a), 当背景孔隙度处于1%~10%时, 反射振幅方位变化几乎独立于背景孔隙度, 对于含水储层(图 14b), 这种独立性成立的孔隙度范围大约变为5%~15%。这指导我们在发育垂直裂缝体的含气页岩或含油水砂岩中, 可利用测井等探测方法估计背景孔隙度, 然后将估计值作为先验常数直接代入反演流程中来预测裂缝密度等参数, 以避免背景孔隙度地震难预测^[110]的问题。分析图 15, 可以总结出含气裂缝储层反射振幅方位变化对裂缝横纵比的依赖度比含水裂缝储层低(即流体模量的降低将减弱反射振幅方位变化对裂缝纵横比的依赖度); 背景孔隙度较大的裂缝储层反射振幅方位变化对裂缝横纵比的依赖度比背景孔隙度较小的裂缝储层低(即背景孔隙度的增大将减弱反射振幅方位变化对裂缝纵横比的依赖度)。所以对于孔隙、喉道发育良好的含气储层, 我们可以将裂缝横纵比设置为常数, 直接利用反射振幅方位变化预测裂缝密度等参数。

综合上述分析, 我们可以总结: 考虑孔-裂隙间的流体流动时, 裂缝走向始终对应方位反射振幅的极大值点, 这为裂缝方位预测提供了理论基础; 但孔、缝孤立时, 裂缝走向与方位反射振幅极值点(极大值点或极小值点, 这两个点的观测方位相差90°)的对应关系受制于充填流体类型等因素, 直接利用反射振幅方位变化规律预测裂缝方位存在不确定性, 此时需要其它预测技术为地震预测裂缝方位提供先验信息。另外, 裂缝密度越大, 地震振幅方位变化越剧烈; 饱含的流体模量越小, 裂缝走向反射振幅与裂缝法向反射振幅的差值越大。这些规律是地球物理学家利用地震方法预测裂缝密度和流体类型的理论基础之一。但地震反射振幅的方位变化规律是裂缝密度、孔裂隙中充填的流体类型、喉道发育情况, 甚至岩石的矿物成分等因素的综合响应, 难以直接利用反射振幅方位变化实现裂缝和流体的精准识别, 所以在分析清楚喉道和等径孔隙发育情况的前提下, 需要研究如何从反射振幅方位变化信息中解耦基质矿物、裂缝及油气等信息, 提出一套“固-液-裂缝”解耦的裂缝储层地震描述及油气检测的方法。除此, 流体流动对地震响应影响较大, 所以我们应该在分析储层孔喉结构和渗透率的基础上, 准确选择岩石物理模型, 防止岩石物理模型的误用带来解释误差。最后, 给出在地震勘探实践中或许有用的两条结论: 对于垂直裂缝发育的含气页岩或含油水砂岩地层, 可将非地震方法的背景孔隙度估计值(甚至常数值)直接代入反演流程中来预测裂缝密度等参数; 对于孔喉发育良好的含气储层, 可将裂缝纵横比设置为常数, 直接利用反射振幅方位变化预测储层参数。这样避免了地震预测背景孔隙度和裂缝纵横比难度大的问题！

3 挑战和机遇

等效各向异性岩石物理理论的快速发展有效推动了裂缝型油气储层的勘探发现, 但随着勘探目标复杂程度的不断提高, 以及油气储层的开发生产活动对地震预测精度的要求不断提高, 裂缝储层等效各向异性岩石物理理论的发展也迎来了新的挑战。结合上文的分析和讨论, 在未来的研究中, 我们应聚焦于以下几方面:

1) 裂缝储层等效介质岩石物理理论目前主要考虑了裂缝夹杂物和各向同性背景之间的弹性互动作用且研究成果得到了应用。但对于页岩等裂缝型储层, 强各向异性背景与裂缝之间的弹性互动鲜有研究;

2) 实际地层中的裂缝形状、大小、表面粗糙度等特征具有较高的复杂性, 应该探索如何克服当前裂缝储层等效介质岩石物理理论的理想假设带来的误差, 建立一套能够更精确地描述实际裂缝特征的模型;

3) 多尺度或跨频带的裂缝介质地震岩石物理理论逐步发展起来, 例如Chapman岩石物理理论, 已经引起地球物理工作者的研究兴趣, 但地震波衰减特征或与地震波频率有关的地震属性和复杂储层裂缝参数的量化关系有待深入研究和应用。