硬币型裂缝介质的频散与衰减

张繁昌 桑凯恒 路亚威

(中国石油大学(华东)地球科学与技术学院， 山东青岛 266580)

0 引言

天然裂缝型油气藏近年来在勘探和生产中引起了较大的关注。天然裂缝控制储层的渗透率，决定油气的分布以及采收率，因此描述波在裂缝介质中的传播特征对于识别裂缝以及提高采收率具有重要意义。

叠前方位反演是目前最常用的识别裂缝方向、估计裂缝密度以及预测流体分布的方法[1-3]，其中的反射系数方程由准静态等效裂缝模型所得，以HTI介质反射系数公式[4]、傅里叶级数展开[5]为代表。现阶段最常用、最直观的裂缝介质等效理论模型分别为线性滑移模型[6]、Hudson[7]刚度扰动硬币型裂缝模型和Thomson[8]各向异性参数裂缝模型，通过忽略裂缝微观结构以及孔缝间的流体交换模拟岩石弹性模量。上述方法在甚高频率下模拟的饱和岩石弹性模量适用于超声频带。对于低频，一般先求得干燥裂缝介质的等效模量，然后利用各向异性Gassmann方程进行流体替换[9]。虽然人们探讨了单组、两组及多组裂缝等效模型及裂缝间的相互影响[10-13]，但仍不能准确地描述波在裂缝介质中的传播特征。

波在介质中传播时发生的频散和衰减是表征裂缝介质弹性特征的重要属性。除介质本身固有的衰减机制外，波在介质中传播时引起的孔隙间、缝隙间以及孔缝间的流体流动也是一种重要的衰减模式。常规准静态裂缝模型只适用于弛豫和非弛豫状态，因而裂缝介质的多尺度、全频带衰减特征研究成为近年的研究热点。Chapman等[14-15]假定每个硬币型裂缝都与固定数量的球形孔隙相连，提出微观到介观尺度频变裂缝裂隙模型。基于White周期层状等效模型[16]，Brajanovski等[17-18]将裂缝等效为无限长平面薄层，借助等效传播矩阵的求解方法，得到了裂缝的法向频变弹性模量，进一步结合层状介质的Backus平均理论，分别得到高、低频极限裂缝介质的频变弹性模量。Galvin等[19-21]将裂缝的影响视为各向同性背景介质中的扰动，将裂缝视为薄硬币形状，利用多重散射定理，给出了随机裂缝介质的法向频变弹性模量。Brajanovski模型和Galvin模型的衰减在低频范围与频率成正比，在高频范围与频率的平方根成正比，与Johnson[22]的认识一致。Gurevich等[23]将Brajanovski模型和Galvin模型统一，给出了全频带裂缝介质弹性模量的统一表达式；Guo等[24]进一步将其拓展到有限厚度平面裂缝模型。另外，借助Krzikalla 等[25]提出的频变张量近似式，Galvin等[26]进一步分析了硬币型裂缝VTI介质的相速度频散和衰减，Yang等[27]和Kong等[28]研究了平面裂缝介质模型的相速度频散、衰减以及黏弹性AVO响应。

本文利用离散积分方法给出了硬币型裂缝法向弹性模量的数值模拟过程，并结合线性滑移理论，分别得到了高、低频弹性模量，进一步分析了油、气、水饱和情况下的频散、衰减以及黏弹性反射响应。

1 硬币型裂缝孔隙弹性理论

对于硬币型裂缝介质的频变性质，Galvin等[21]提出了数值模拟的Fredholm方程

(1)

式中：x、y为积分变量；S(x)=(2/π)[sin(ax)-axcos(ax)]/x2;R(x,y)=sin[a(x-y)]/(x-y)-sin[a(x+y)]/(x+y),a为硬币型裂缝长度；B为待求函数；p0为流体压力；T(y)为中间函数

(2)

硬币型裂缝介质在垂直方向的频变弹性模量p33与B的关系为[21]

(3)

式中：ρ为介质密度；ε为裂缝密度；k1为快纵波的波数。

式(1)、式(3)即为Galvin等[21]的硬币型裂缝数值模拟表达式，利用常规积分方法求解难度较大。为此，下面深入分析数值求解方法。

2 数值求解方法

由文献[20]可知，除了直接由式(1)求解B外，还可以使用方程

(4)

求取。式中：θ(ζ)为待求函数；Γ(ζ，ξ)为核函数，满足

(5)

式中J为Bessel函数。θ(ζ)与B(x)的关系为

(6)

利用式(4)～式(6)也可以计算B，但相比式(1)，式(4)～式(6)在形式上更复杂。在实际求解中，发现式(6)可以避免奇异点的影响，而且求取极限的步骤可以合并，进而基于Gauss-Lobatto离散积分，通过数值求解得到频变弹性模量。

首先，取Bessel函数的高阶项

对核函数Г化简

(7)

对于区间[0,a]的积分函数f(x)可以离散为

(8)

式中：xj∈[0,a]为求积节点；Aj为求积系数。式(1)可以离散为

(9)

令x=x1,…,xn, 并记θi=θ(xi),Γij=Γ(xi,xj),Fi=F(xi), 则式(9)可写为

(10)

将上式表示为矩阵形式

(11)

可以看出，求解式(11)可以得到每一个节点的θi值。需要注意，Г函数是T(y)sin(ζy)sin(ξy)在区间(0，∞)积分的震荡函数，当积分变量y的值较大后积分函数快速衰减。图1为当频率为10、100、10000Hz时积分函数随lgy的变化曲线。由图可见，当积分变量y>104后，积分函数趋于0，因此为加快计算速度，积分区间一般取(0，104)。

进一步

(12)

上式将求极限问题转变为求积分。利用Gauss-Lobatto方法对上式离散化

(13)

只要将式(11)得到的每一个节点处的θ(xj)值乘以对应的节点值和积分系数，最后求和就可以得到极限值。

图1 当频率为10(a)、100(b)、10000Hz(c)时积分函数随lgy的变化曲线

3 硬币型裂缝介质频变弹性张量

Galvin等[21]只给出裂缝法向的弹性模量，其他方向的弹性模量可以通过Krzikalla 等[25]推导的频变弹性模量关系得到

(14)

式中Cij为刚度系数，下标i,j=1,…,6，上标“high”代表高频，“low”代表低频。可见，除了已知法向弹性模量p33，还需知道高、低频极限条件的其他各个弹性参数。

3.1 低频极限

研究低频条件的流体充填孔隙介质的等效弹性模量时，一般利用Gassmann方程[29]对干岩石骨架进行流体置换分析。

对于干岩石的描述，Galvin和Hudson裂缝模型给出的结果相同。干硬币型裂缝的法向弱度ΔN和切向弱度ΔT分别为[30]

(15)

(16)

结合线性滑移模型，将干岩石和裂缝参数代入Gassmann方程，可以得到其他弹性参数

(17)

式中λ、μ为干岩石的拉梅参数。

3.2 高频极限

当波在介质中传播时，会引起应力变化。在高频条件下惯性占据主导，因此流体没有足够的时间在孔隙和缝隙间流动使介质中的应力达到平衡状态，此时流体流动产生的影响可以忽略，可将裂缝和孔隙视为独立存在、互不影响。岩石的高频弹性模量可模拟为饱和背景岩石与饱和裂缝弹性模量的叠加，流体饱和裂缝法向弱度ΔN由Hudson模型模拟[30]

(18)

式中：μf为流体的剪切模量(近似为零)；c为裂缝厚度。

饱和裂缝的切向弱度等于干岩石裂缝的切向弱度ΔT(式16)。

将饱和背景岩石与饱和裂缝的弹性模量叠加，得到裂缝介质的高频模量

(19)

4 黏弹性VTI介质的相速度和反射透射系数

VTI介质中复相速度VP的频散表达式为[31]

(20)

其中

D={[(p11-p44)sin2θ-(p33-p44)cos2θ]2+

黏弹性VTI介质纵、横波反射、透射系数满足下式[32]

(21)

式中：R为反射系数；T为透射系数；下标P代表纵波，S代表横波，1代表界面上方，2代表界面下方。其中

(22)

(23)

W=p55(γsx+βsz)

(24)

Z=βp13sx+γp33sz

(25)

(26)

式中sz为垂直方向的慢度，和水平方向慢度sx及刚度的关系为

(27)

(28)

(29)

5 数值模拟结果分析

首先利用离散积分法模拟硬币型裂缝介质饱含水的法向频变弹性模量，VTI裂缝介质参数如表1所示，流体参数如表2所示，其频散与衰减如图2所示。

表1 VTI裂缝介质参数

由图可见：①在弛豫状态，由于波长远大于裂缝尺度，由波动引起的应力变化在很短时间达到平衡，速度几乎不变；在非弛豫状态，应力变化加快，流体的惯性逐渐成为主导，速度趋于稳定；在过渡状态，因为岩石中的流体流动减弱，整体柔度减小，对应刚度增加，因而随频率增加，速度不断变大(图2左)。②当裂缝密度增加时，岩石柔度增加，速度变小，对应的衰减整体变大，但衰减峰值频率基本不变(图2右)。③如果裂缝密度不变，裂缝长度由小变大，速度与衰减曲线的趋势基本不变，只是整体左移。

图2 不同裂缝密度(上)、裂缝长度(下)的硬币型裂缝的频散(左)与衰减(右)

表2 流体参数

图3为不同频率的含水介质速度和衰减随入射角的变化，图4为不同入射角的含水介质频散和衰减。由图可见：在垂直于裂缝方向(θ=0°)速度和衰减变化最大，在平行裂缝方向(θ=90°)趋于定值；随入射角变大，各向异性程度减小，速度的频散趋势减缓，衰减峰值也逐渐减小。 图5为含气介质的速度与能量随入射角和频率的变化。由图可见：由于含气介质的密度较含水介质小，速度整体变小(图5左)；由于气体黏度较低，当频率较高时，在短时间内通过流体交换达到平衡，因而弛豫和非弛豫之间的过渡区向高频移动，振幅峰值也向高频移动(图5右)。

图6为含油介质的速度与能量随入射角和频率的变化。由图可见：速度变化与饱水介质相似(图6左)；由于油的黏度较高，在较低频带就达到非弛豫状态，因而峰值频率向低频移动(图6右)。

图7为含水、含油、含气介质频变反射系数与流度。由图可见：①由于上层波阻抗大于裂缝介质，垂直反射系数为负，因此在弛豫与非弛豫过渡区，随频率增加，裂缝刚度增加、岩石的波阻抗增加，导致反射系数幅值减小(图7左)。②当裂缝介质饱含气时，在入射角较小范围，反射系数幅值明显大于含水、含油(图7下左)。③流度的峰值频率与能量的峰值频率相对应，即含油时峰值流度最大(图7中右)，含气时峰值流度最小(图7下右)；随入射角增大，各向异性程度减弱，整体幅值减小。

图3 不同频率的含水介质速度(左)和衰减(右)随入射角的变化

图4 不同入射角的含水介质频散(左)和衰减(右)

图5 含气介质的速度(左)与能量(右)随入射角和频率的变化

图6 含油介质的速度(左)与能量(右)随入射角和频率的变化

6 结论

本文通过岩石物理分析，研究了VTI硬币型裂缝介质频散、衰减规律的数值模拟方法，并分析了不同流体饱和介质的频散、衰减现象，得到以下认识：

(1)通过离散积分和高阶近似方法，将第二类Fredholm积分方程与多重散射方程的远场极限值求解相统一，为数值模拟硬币型裂缝模型的法向频变弹性模量提供快速、准确的方法；结合线性滑移理论和Hudson裂缝模型，得到了频变弹性模量的解析式。

(2)分析VTI介质相速度与黏弹性反射、透射系数方程可知，裂缝长度决定弛豫向非弛豫过渡区的频带，裂缝密度影响衰减峰值；当裂缝介质中饱含不同流体时，流体的黏度成为决定弛豫向非弛豫过渡的又一因素，饱气时反射系数幅值大于饱油、饱水时，流度的变化趋势与衰减的变化趋势相似。