利用各向异性模型计算水平地应力

张兆谦

(大庆油田有限责任公司勘探开发研究院,黑龙江大庆163712)

0 引 言

目前,测量地应力的方法有许多,如岩石水压裂法、岩石声发射法等,前者代价昂贵、耗时长,后者岩心代表性差,无法得到连续的剖面[1]。利用测井方法计算地应力,具有成本低、准确性高的特点,同时可以得到连续的井剖面,因此在地应力的评价方面有着广泛应用。地层中绝大部分岩石具有不同程度的各向异性[2],而以往测井计算水平地应力大多采用各向同性模型,没有考虑地层各向异性。由于致密油气储层岩石弹性参数及地应力各向异性,基于各向同性假设的储层完井工程品质评价与实际存在着差异,为了接近地层的真实情况,必须考虑地层各向异性。

本文利用阵列声波测井中斯通利波包含的水平横波的信息,根据频谱加权平均慢度定理反演水平横波的慢度(或计算反演水平剪切模量C66),利用这一重要参数得到水平方向的弹性模量和泊松比,进而实现了利用各向异性模型计算水平地应力。

1 运用斯通利波反演水平剪切模量

各向异性地层中(水平各向同性TIV),根据广义胡克定律,完全独立的弹性参数有5个,分别是C66、C11、C13、C33和C44。根据弹性波动力学理论,地层垂直方向的纵、横波速度分别为

(1)

地层水平方向纵、横波速度分别为

(2)

式中,C66、C11、C13、C33和C44分别为地层弹性刚度系数,GPa;vcv、vsv分别为地层垂直方向的纵、横波速度,km/s;vch、vsh分别为地层水平方向的纵、横波速度,km/s;ρ为地层密度,g/cm3。

直井中,阵列声波测井只能测量得到垂向的声波速度,无法得到水平方向的声波速度。井中传播的斯通利波对水平方向传播的横波敏感。在低频情况下,斯通利波的传播引起了井周的扭曲,斯通利波可以看作是井眼上的压力克服水平剪切应力而发生的形变,因此,与水平剪切模量C66有关。

1.1 斯通利波频散方程

在多极子阵列声波测井仪器中,为了隔离偶极子沿着仪器传播的挠曲波,仪器安装有隔声体,这使仪器具有相当大的弹性。仪器的弹性对单极子声波测井中的斯通利波的传播有很大影响。因此,分析声波测井中的斯通利波一定要考虑到仪器的弹性效应。根据唐晓明等[3-7]的理论研究,有在TIV介质中充满流体、仪器居中情况下的单极子声波测井中斯通利波的频散方程

D(k,w,C66,C44,MT,R,a,ρ,vf,ρf,C11,C13,C33)=0

(3)

在井参数、流体参数和纵横波各向异性为已知的情况下,可以简写为

(4)

式中,k为轴向波数;w为角频率,Hz;vf为井中流体声波速度,km/s;ρf为流体密度,g/cm3;MT为仪器模量,Pa;R为井径,cm;a为仪器半径,cm;ρ为地层密度,g/cm3。频散方程(4)中各参数的具体形式及推导详见参考文献[4]。

频散方程中,给定角频率w,该方程是关于k的非线性方程,采用牛顿拉弗森数值方法求解。在求解过程中,把频率范围限制在斯通利波的频带范围内,解上面的频散方程,求出波数k对应于斯通利波的根kst,得到了斯通利波相速度或者相慢度的频散曲线

Sst(w)=1/vst(w)=kst/w

(5)

1.2 频谱加权平均慢度定理

声波的频散意味着对应不同的频率,频散波的速度不同,斯通利波是典型的频散波,但是运用波形相似相关方法处理得到斯通利波慢度曲线,是假设所用的声波不存在频散效应。频谱加权平均慢度定理表明[8-9],从时间域的非频散阵列处理方法得到的慢度实际上是对频率域慢度频散曲线的加权平均,其中权重系数为w2A2(w)。

(6)

从式(6)可见,利用阵列声波测井数据可以有效地计算出水平方向的横波速度vsh或水平剪切模量C66。计算中需要的其他参数可以从测井数据中得到。C11、C13这2个弹性系数可以通过岩石力学实验刻度,如果没有实验数据支持,那么根据参考文献[4],有C11=C33(C66/C44),同时C13=C33-2C44。虽然这些假设都是近似,但是它们不太影响C66的计算,Ellefsen[10]给出了这些参数灵敏度的详细分析结果。

利用公式(6)构造反演的目标函数

(7)

对目标函数求其极小值对应的C66,极小值求解过程采用单变量函数的黄金分割寻优方法。

1.3 仪器模量刻度

测井仪器的存在对斯通利波的传播有着很大影响[4],必须对仪器模量MT进行刻度。对于一定井段,各向异性情况是已知的(利用偶极声波测井各向异性分析得到),特别是找到各向同性的地层段,地层垂向上的纵横波慢度可以用软件的速度分析模块得到。利用式(1)计算得到C33、C44,各向同性处C66=C44、C11=C33,井中流体参数可以通过钻井液测井得到,式(3)中唯一未知的就是仪器模量MT。刻度MT的过程就是找到一个合适的仪器模量值使得其满足式(8)

(8)

一旦仪器模量参数得到,可以利用式(7)反演整个井段的水平剪切模量C66。

2 方法实现

地层的纵横波、斯通利波慢度及斯通利波的波形曲线利用商业软件包可以得到,该研究主要工作集中在求解频散方程及利用频散曲线反演C66。图1反演C66的流程图,

在A井中取一深度点, 根据上述已知参数利用牛顿拉弗森数值方法求解频散方程,可以得到斯通利波的频散曲线。图2为A井2 170.8 m处的斯通利波波形频谱及频散曲线。

图1 反演C66流程图

图2 斯通利波波形频谱及频散曲线

利用斯通利波频散曲线及频谱数据,结合构造的目标反演函数,可以反演出水平剪切模量C66。垂直方向纵横波慢度由MAAS软件的速度分析处理模块得到,水平方向参数由反演方法得到,其中C11=C33(C66/C44)。利用式(2),可以得到水平方向纵横波慢度。C66代表了水平剪切模量,就是代表了水平横波的速度,准确得到这个参数,为下一步计算水平方向的弹性模量和泊松比,进而利用各向异性模型计算水平地应力打下基础。

3 各向异性模型计算水平地应力

3.1 各向异性地应力计算模型

各向同性地应力计算模型认为垂直和水平方向弹性模量是一致的,而各向异性考虑了垂直和水平方向的纵横波差异,即弹性模量和泊松比不同。利用按照上述方法反演的水平剪切模量C66及其他弹性模量,分别是C12、C11、C13、C33和C44,其中C12=C11-2C66,利用这些参数,计算垂直及水平方向弹性模量、泊松比,计算公式为

(9)

(10)

(11)

(12)

目前常用的计算地应力的模型是孔隙介质模型,该模型垂向上考虑了上覆岩石压力及孔隙压力,水平上考虑了构造残余应力的作用,其水平地应力计算公式为[11]

(13)

(14)

式中,Ev、Eh分别为垂直及水平方向弹性模量,MPa;μv、μh分别为垂直及水平方向泊松比,无量纲;σh、σH分别为水平最小、最大地应力,MPa;pp为地层孔隙压力,MPa;α为Biot系数,无量纲;εH和εh分别为水平方向最大、最小构造应力系数。

3.2 计算实例

AD为大庆深层致密气区块,储层岩性种类较多,非均质性强,岩石造岩矿物成分变化大,主要岩石类型包括砾岩、砂砾岩、砂质砾岩等, 有效孔隙度主要集中在4%～9%,平均5.6%;空气渗透率主要集中在0.03～0.5 mD(2)非法定计量单位,1 mD=9.87×10-4 μm2,下同,平均0.147 mD,为致密储层,单井必须压裂才有产能。因次,准确计算该区块地应力大小将为压裂施工起到很大的作用。该区块有12块水平地应力实验测量数值(见表1),另外有56块样品动静态模量参数,38块样品的Biot系数。由式(13)和式(14)结合实验室测量的最大、最小水平地应力值及其他参数可以得到εH和εh这2个构造应力系数,该区块最大、最小构造应力系数分别为0.000 16和0.000 78。

图3为致密气区块SS4井测井计算水平地应力深度剖面。第1、第2道为常规测井剖面,第4道为岩石体积剖面。第5道、第6道分别为各向异性模型计算与岩心分析最大、最小水平地应力,可以看出数值一致,证明测井计算的合理、准确。第7道为利用式(9)和式(10)计算的垂直和水平方向的弹性模量,可以看出,2 769～2 783 m储层段内弹性模量的值比较接近。第8道、第9道为地层弹性刚度参数,第10道中,ΔtCh和ΔtCv分别代表横向和垂向的纵波慢度,第11道中ΔtSh和ΔtSv分别代表横向和垂向的横波慢度,可以很明显的看出,水平和垂直方向的纵横波慢度不一致,这代表了地层真实的情况。2 769～2 783 m段为该井气层段,测井计算最小水平地应力为62 MPa,上下围岩为泥岩,其最小水平地应力为85 MPa,纵向应力差为23 MPa,围岩起到很好的应力遮挡作用,压裂后日产气4.3×104m3,为工业气层。

表1 AD区块实验测量水平最大(最小)地应力及地应力梯度

图3 各向异性模型计算与岩心分析水平地应力精度对比

4 结 论

(1)直井中,利用阵列声波测井资料反演水平剪切模量C66是可行的,利用反演得到的C66和测井得到的地层弹性刚度系数可以计算垂直和水平方向的弹性模量、泊松比,进而利用各向异性模型计算地层水平地应力。

(2)反演水平剪切模量C66之前,必须掌握待反演C66的井段各向异性情况,找出各向同性的区间,刻度出仪器模量MT。

(3)对大庆AD区块的5口井计算了最大、最小水平地应力,计算结果与岩心分析的结果进行了比较,数据显示非常好的一致性,验证了该方法计算结果的可靠性。