近年,声弹性理论[1-2]的发展和应用为地下应力的检测提供了新途径.声弹性理论被广泛应用于井孔问题[3-8],实验室测量[9]以及地下矿井[10].声弹性理论和交叉偶极声测井技术[11]的结合为地应力检测开拓新途径,Sinha等[4-5]建立了以井孔模式波声弹性模型判断和反演地应力的方法和应用研究工作,并给出了从井孔模式波同时反演三阶弹性模量的方法.这一技术将油气勘探开发中的测井资料很好地利用了起来.
然而,利用交叉偶极声测井检测异常地应力依据的现行井孔声弹理论和推断地应力的判据,都是在假定未受到异常应力作用前地层是各向同性介质的前提下导出的.而实际地层却常常是未受到外力作用前就已经是各向异性介质,这是利用交叉偶极声测井方法检测地应力时不能回避的问题,因此有必要在理论上进行研究.由于对于油气藏具有重要意义的周期薄互层与裂隙定向分布的介质模型在长波假设下均等效于横向各向同性介质模型[11-12],使在勘探地理物理领域普遍注重于含有5个弹性模量的横向各向同性的介质模型.自Norris 等[3]和Sinha[4]等提出井孔体系的声弹性理论以来,无论对裸眼井[4-9]还是套管井[13-14],利用正交偶极声测井资料求取地层异常应力时,都假定地层在未受到偏离应力作用之前是各向同性介质.这些研究是否能直接应用到参考状态为横向各向同性的介质呢? 显然,对参考状态为横向各向同性地层的井孔声弹问题,有必要进行理论和数值的研究与分析.李刚等[15]在零级场和三阶弹性模量都取各向同性近似的条件下,对横向各向同性对称轴与井轴垂直(HTI)的井孔声弹问题进行了研究;与HTI井况不同,横向各向同性对称轴与井轴一致(VTI)的情况,在实际的测井数据处理中,不会出现横波分裂,但可以进一步根据斯通利波来判断[16]或者根据弯曲波和斯通利波的频散曲线识别[17].VTI井况的弯曲波的声弹效应如何?在理论上尚未见论及.对于各向同性介质,独立的三阶弹性模量只有3个,而对于横向各向同性介质其独立的三阶弹性模量是9个[18],因此横向各向同性介质的声弹效应比各向同性情况复杂.本文就针对VTI井况的声弹效应对弯曲波的影响进行了研究和分析,为在VTI井况利用交叉偶极声测井检测地应力提供理论基础.
2 横向各向同性(VTI)介质井孔弯曲波声弹性公式声弹性理论的前提是在有限静形变(偏置)上叠加小(线性)扰动.描述介质的位移可以用参考坐标、中间坐标和当前坐标,这三种坐标分别对应于介质的三种不同状态,即未变形状态、静形变偏置状态和在该状态上迭加小扰动的形变状态.与其对应有三种不同形式的描述介质质点的运动方程,针对不同的问题,可以应用其中任意一个和相应的边界条件进行求解.一般井下岩层所受的应力状态比较复杂,可以分解成均匀静水压力与偏离应力的和.均匀静水压力的大小可以等于在某一深度点主应力的平均压力或上覆岩压力,偏离应力可以是产生井周方位各向同性的井孔压力及引起井孔应力集中的单轴应力(或双轴应力).本文将岩层受均匀静水压力时的状态选为参考状态,偏离应力与均匀静水压力相比其量值要小得多[6].本文选参考坐标及相应的运动方程来描述形变物理量.
假设所选的参考状态具有横向各向同性(TI)介质特征.类似于参考态为各向同性介质情况[6],根据沿井轴方向传播的模式波的主要能量束缚在井壁周围的特点,可以利用摄动积分法求解和分析井孔模式波的频散受外加偏置应力的影响,得到由偏置静应力和应变引起的弯曲波相速度变化的近似公式
(1) |
角标fl代表弯曲波,vR fl和vfl分别表示在参考状态和偏置应力诱导的当前状态的弯曲波的相速度;重复的张量下角标代表求和,逗号代表对空间位置的偏导数,上标*号代表取复共轭.V是积分体积,它包括模式波在流体和固体中延伸的体积;ρ0 是未受应力状态的密度,uγfl 表示参考状态下弯曲波与特征频率ω 对应的位移分量.参考态弹性模量的增量为
(2) |
(3) |
(4) |
其中,δγυ 是Kronecker符号,TLM、EAB和Wγ,K分别为偏置应力、应变和静位移梯度.cLγMυ 是二阶弹性模量,cLγMνAB是三阶弹性模量,对于横向各向同性介质相应独立的二阶弹性模量分别为c11,c13,c33,c44和c66,独立的三阶弹性模量是c111,c112,c113,c123,c133,c144,c155,c333和c344[18].
假设井内压力是P0,它是相对于参考状态的井内流体静压力的反作用力.沿x轴和y轴的最大和最小水平主应力为SH和Sh,r0 是井孔半径,θ 是以x方向为极轴确定的空间场点的方位角.根据静应力作用平衡方程及平面应变的假设,可以得到以井轴为z轴的柱坐标系下,受井孔压力和水平面内双轴应力的井周应力分布.由于径向、环向和切向应力与各向同性情况相同[8],与介质的弹性模量无关,限于篇幅,没有给出.轴向应力Tzz有别于各向同性情况:
(5) |
这里ν′ =c13/(c11 +c12).相应的静应变为
(6) |
这里ν=c12/(c11+c12)是泊松比.可以看出,横向各向同性介质相应的静应变Err和Eθθ 不同于各向同性介质[6].在柱坐标系下,对(1)式分子、分母的积分,可得到弯曲波相速度相对改变量随外加应力的变化公式,即地层为VTI情况的井孔弯曲波声弹性公式
(7) |
(8) |
(9) |
这里Cji(i=P,S,D,j=0~8)是应力诱导的弯曲波速度改变量的灵敏系数[8],它们与5个二阶弹性模量有关,是在参考状态弯曲波位移场分布uγfl 已知时由体积积分得到的与频率有关的量.上角标P,S和D 分别代表井内压力、主应力的和与差,下角标R 代表参考状态.φ 是相对于x轴偶极源偏振方位角.由公式(7)可知受水平地应力作用后,井孔弯曲波速度的变化依赖井内压力、主应力的和、差及方位角;式(8)表明由内压引起的井孔弯曲波的声弹性公式与6个独立的三阶弹性模量有关,而且与三阶弹性模量c111和c112、c113和c123、c144和c155有关的灵敏系数分别大小相等,符号相反.式(9)的弯曲波声弹公式与7个独立的三阶弹性模量有关,当不顾及三阶弹性模量的各向异性时,根据各向同性三阶弹性模量的关系[1],上式(8)和(9)可简化为
(10) |
(11) |
这里C7PI =C7P -2C5P,C8PI =4C1P -C7P,C5iI =C1i +C2i +C3i +C5i +C6i,C8iI =4C1i+C8i,C7iI =2C1i+2C2i+2C3i+C6i+C7i.上角标I代表三阶弹性模量为各向同性情况的灵敏系数.
与各向同性情况相似[13],由平面波解和平面应变近似,可以导出VTI井况受双轴应力作用时沿井轴z传播的平面纵、横波声弹公式
(12) |
这里
式中ρ0 是参考态的物质密度,Trr,Tθθ,Trθ 是应力,θ 是与最大水平主应力的方位角,ν′ =c13/(c11 +c12).VIJ 指沿XI方向传播沿XJ 方向偏振的平面波速度.由公式(12)可知沿井轴传播的横波主要与三阶弹性模量c144和c155有关.我们知道在低频弯曲波受远场速度支配,高频弯曲波受近井壁区域支配,因此为保证三阶弹性模量c144和c155在低频的主导地位,弯曲波声弹性公式灵敏系数中的Ci7和C8i(i=S,D)在低频应比其它的灵敏系数Cji(i= S,D,j=1~6)大.
3 数值计算分析与各向同性情况相同[6, 8],横向各向同性介质(VTI井况)井周应力集中的复杂性及介质参数的影响也体现在声弹性公式速度改变量的灵敏系数Cji(i=P,S,D,j=0~8)中.本文数值计算了由应力引起的灵敏系数的频率响应、两种偏振的偶极弯曲波的相速度频散曲线和平面横波速度随半径的变化曲线.
计算采用的地层参数和三阶弹性模量见表 1和表 2.其它参数:井内流体密度为ρf=1000kg/m3,井内流体声速为vf=1500m/s,ρ 为井外固体密度,井孔半径取0.1016m.双轴应力SH =-9 MPa,Sh=-3 MPa,井内流体压力P0= -2 MPa (简称内压),负号表示外加应力是压力.参数A是Berea砂岩[19],参数B是页岩[20],它们的三阶弹性模量分别是Sarkar 等[19] 和Prioul 等[20] 在假定其为各向同性前提下测得的.表 2的参数是根据文献[19-20]提供的误差分析适当调整形成的横向各向同性介质的三阶弹性模量.
图 1和2分别是地层参数A和B的弯曲波速度改变量的灵敏系数随频率变化曲线.其中图a是与内压有关的灵敏系数,图b是与双轴应力和有关的灵敏系数,图c是与双轴应力差有关的灵敏系数;图中黑线是顾及了三阶弹性模量的各向异性,红线(上角标I)代表与各向同性三阶弹性模量(c144,c155和c123)有关的灵敏系数.对井内压力系数(QP)中与TI介质的三阶弹性模量有关的灵敏系数考察发现:6个灵敏系数在低频较小,高频较大.对双轴应力和(QS)与差(QD)中与TI介质的三阶弹性模量有关的灵敏系数考察发现:在≤5kHz的低频区(这正是交叉偶极声测井仪的激发频带),与c144和c155有关的灵敏系数要比其它的灵敏系数大1~3个数量级,这与公式(12)的认识一致;双轴应力差(QD)中与c144、c155有关的灵敏系数(C7D ,C8D )符号相反;而且在高频,特别是双轴应力和(见图b)c144的灵敏系数比其它灵敏系数甚至大1 个数量级.当不顾及三阶弹性模量的各向异性时,考察结果表明C7iI 和C8iI(i=S,D)与TI介质的c144、c155有关的灵敏系数C7iC8i(i= S,D)变化规律相似,特别在2kHz左右的低频几乎重合.与各向同性介质相似[6],三阶弹性模量无关的灵敏系数C0P,C0S 和C0D (纵坐标见图 1,2左侧)和三阶弹性模量有关的灵敏系数在同一数量级,因此对相速度起控制作用的是与三阶弹性模量有关的项.由以上分析可知当c144和c155与其它三阶弹性参数为同一量级甚至小一量级的参数值时,在低频弯曲波速度的变化主要受c144和c155的影响,而高频弯曲波的变化依赖于7个三阶弹性模量的大小.
图 3是由公式(7)计算的受水平双轴应力作用,弯曲波偏振方向与最大主应力一致(φ=0°)和垂直(φ=90°)的弯曲波频散曲线.图a是地层参数A;图b是地层参数B.粗虚线代表参考状态为横向各向同性弯曲波的频散曲线,实线代表顾及三阶弹性模量各向异性的频散曲线,点线是仅顾及c144和c155两个三阶弹性模量的频散曲线.由图 3可知,双轴应力使弯曲波相速度的改变量随频率增加明显减小,对φ=90°双轴应力引起的频散曲线的变化也近似于参考状态频散曲线的平移;由双轴应力诱导的互相垂直偏振的弯曲波频散曲线有交叉.以上认识与各向同性介质受应力作用相似[6],因此交叉现象仍是判断地应力存在的标志.由点线和实线的对比可知忽略五个系数(Cji,i=D,S,j=0,2,3,5,6)和由方程(7)求得的频散曲线在低频(参数A约≤3.5kHz,参数B约≤5kHz)基本一致,这也证实了弯曲波速度的变化在低频区主要受c144和c155的影响.图 3a低频在方位φ=0°横波速度约为1927 m/s,在方位φ=90°的横波速度约为1779 m/s.图 3b 低频在方位φ=0°横波速度约为1565 m/s,在方位φ=90°的横波速度约为1524m/s.
图 4是根据式(12)计算的井孔应力集中情况下速度随半径的变化曲线,其中图a为地层参数A,图b为地层参数B.由图 4 可知远场的速度几乎与角度θ 无关,而在小半径范围内,速度与角度θ 有关,径向偏振的横波速度对应于90°方位的V32(细线)及0°方位的V31(粗线).在近井孔附近快偏振方向垂直于施加最大水平主应力的方向,而在离开井孔并大于某径向位置,快偏振方向与最大水平主应力的方向一致,径向偏振的横波出现交叉.由图 4a在远场方位θ=0°横波速度约为1922 m/s,方位θ=90°的横波速度约为1982m/s,这与图 3a低频的弯曲波速度基本一致.由图 4b在远场方位θ=0°横波速度约为1560m/s,在方位θ=90°的横波速度约为1522m/s,这与图 3b在低频的弯曲波基本一致.数值结果再次验证了在低频弯曲波速度的变化主要受c144和c155影响.
图 5是受水平双轴应力作用,但不顾及三阶弹性模量的各向异性(取表 2中的c144,c155和c123的参数值)的弯曲波频散曲线.其中粗虚线代表参考状态为横向各向同性弯曲波的频散曲线,实线代表顾及三阶弹性模量各向异性,点划线是各向同性的三阶弹性模量.对比可知顾及和不顾及三阶弹性模量各向异性的频散曲线在约小于5kHz的低频一致性非常好,这与各向同性的灵敏系数的考察结果一致(见图 1和2中红线,上角标I).因此,在缺乏足够实验条件情况下,对TI介质(VTI井况)以c144、c155 和c123三个独立的量进行测量,然后可暂不考虑三阶弹性模量的各向异性,建立简化的反演公式,从而实现利用交叉偶极声测井检测异常地应力.
图 6 是井孔压力对偶极弯曲波频散曲线的影响.粗线为参考状态为横向各同性介质的弯曲波相速度,虚线是受井孔压力,双轴应力共同作用的两相互垂直方向偏振弯曲波相速度,细线是仅受双轴应力作用的弯曲波的频散曲线.本文仅就地层参数A的情况进行了对比.由图 6可知,受双轴应力和内压共同作用和仅受双轴应力作用的的频散曲线都有交叉现象,交叉点都约为7kHz,有内压作用的交叉点的相速度值变大.这与图 2a与内压有关的灵敏系数其低频小、高频大一致.由此可以看到,与各向同性介质相似[6],横向各向同性介质(VTI井况),井孔压力的大小会改变交叉点的相速度,而该点对应的频率不变.
对井外具有9个独立三阶弹性模量的横向各向同性介质(VTI井况)受到井孔静压力和井外水平双轴应力的作用,本文给出了偶极弯曲波声弹性公式,即利用摄动积分给出了井孔弯曲波的相速度变化与内压P0、双轴应力和(SH +Sh)/2、双轴应力差(SH -Sh)/2以及与声源取向φ 的关系式.理论推导的结果是弯曲波的声弹公式与5个独立的二阶弹性模量以及7 个独立的三阶弹性模量c111、c112、c113、c123、c133、c144和c155有关;由井内流体压力引起的井孔弯曲波的声弹性公式与c133无关;根据平面应变近似推导了横向各向同性介质(VTI井况)受水平双轴应力作用的平面纵、横波声弹公式,此公式表明沿井轴传播的横波速度变化主要与三阶弹性模量c144和c155有关.
对弯曲波的考察表明:由与双轴应力有关的灵敏系数和弯曲波相速度的变化分析可知弯曲波速度的变化在低频主要与c144和c155有关,而高频弯曲波相速度的变化受7 个独立的三阶弹性模量共同影响;两互相垂直方向偏振的弯曲波频散曲线的交叉仍是判断地应力存在的标志;而且上述认识与径向偏振的平面横波一致;不顾及三阶弹性模量的各向异性,其灵敏系数C7iI 和C8iI (i=P,S,D)与TI介质的c144、c155有关的灵敏系C7i 和C8i(i=P,S,D)变化规律相似,在低频(≤5kHz)弯曲波的频散曲线与顾及三阶弹性模量各向异性基本一致;对仅有井内压力时与三阶弹性模量有关的灵敏系数,在低频较小,高频相对较大;内压对低频影响小,对高频影响大,内压会影响交叉点的速度,但不影响交叉点的频率.显然,各向同性地层井孔声弹效应的诸多性质也适应于横向各向同性介质(VTI井况).而无论是各向同性介质还是各向异性介质在利用交叉偶极声测井检测异常地应力时,三阶弹性模量c144和c155的测量至关重要.在缺乏足够实验条件情况下,对TI介质(VTI井况)以c144、c155和c123三个独立的量进行测量,然后可暂不考虑三阶弹性模量的各向异性,建立简化的应力反演公式.反之,如果已知地层的地应力信息,由简化的声弹公式可以反演三阶弹性模量c144、c155和c123.参考态为横向各向同性(VTI井况)的研究结果对利用测井数据估测地应力或三阶弹性模量提供了有意义的指导.
致谢感谢中国科学院声学研究所张海澜研究员和陈浩博士在本文工作中给予的建议.
[1] | Pao Y H, Sachse W, Fukuoka H. Acoustoelasticity and ultrasonic measurements of residual stresses. //Mason W P, Thurston R N, eds. Physics Acoustics. New York: Academic Press, 1984, 17: 61-143. |
[2] | 田家勇, 胡莲莲. 固体声弹性理论、实验技术及应用研究进展. 力学进展 , 2010, 40(6): 652–662. Tian J Y, Hu L L. Progress in theories and experimental technologies of solid acoustoelasticity and its application. Advances in Mechanics (in Chinese) , 2010, 40(6): 652-662. |
[3] | Norris A N, Sinha B K, Kostek S. Acoustoelasticity of solid/fluid composite systems. Geophysics. J. Int. , 1994, 118(2): 439-446. DOI:10.1111/j.1365-246X.1994.tb03975.x |
[4] | Sinha B K, Kostek S. Stress-induced azimuthal anisotropy in borehole flexural waves. Geophysics , 1996, 61(6): 1899-1907. DOI:10.1190/1.1444105 |
[5] | Sinha B K. Inversion of borehole dispersion for formation stresses, Proc. 1997 IEEE International Untrasonic Symposium , 1997: 781-786. |
[6] | 曹正良, 王克协, 李刚, 等. 井孔压力与地层应力集中诱导的偶极弯曲波的分裂. 地球物理学报 , 2003, 46(5): 712–718. Cao Z L, Wang K X, Li G, et al. Dipole flexural waves splitting induced by borehole pressurization and formation stress concentration. Chinese J. Geophys. (in Chinese) , 2003, 46(5): 712-718. |
[7] | 陈浩. 预应力地层井孔声场研究及其应用. 北京: 中国科学院声学研究所, 2007. Chen H. Borehole Acoustic Wave Fields in Pre-Stressed Formation and Its Application (in Chinese). Beijing: Institute of Acoustics, Chinese Academic of Sciences, 2007. |
[8] | Cao Z L, Wang K X, Ma Z T. Acoustoelastic effects on guided waves in a fluid-filled pressurized borehole in a prestressed formation. J. Acoust. Soc. Am. , 2004, 116(3): 1406-1415. DOI:10.1121/1.1777857 |
[9] | Winkler K W, Liu X Z. Measurements of third-order elastic constants in rocks. J. Acoust. Soc. Am. , 1996, 100(3): 1392-1398. DOI:10.1121/1.415986 |
[10] | Bakulin V N, Bakulin A V. Acoustopolarizational method of measuring stress in rock mass and determination of Murnaghan constants. 69th Annual International Meeting, SEG, Expanded Abstracts , 1999: 1971-1974. |
[11] | 陶果, ChengC H, ToksözM N. 应用正交偶极子测井资料测量EDA地层的横波各向异性. 地球物理学报 , 1999, 42(2): 277–286. Tao G, Cheng C H, Toksöz M N. Measurements of shear-wave azimuthal anisotropy with cross-dipole logs. Chinese J. Geophys. (in Chinese) , 1999, 42(2): 277-286. |
[12] | Zhang B X, Dong H F, Wang K X. Multipole sources in a fluid-filled borehole surrounded by a transversely isotropic elastic solid. J. Acoust. Soc. Am. , 1994, 96(4): 2546-2555. DOI:10.1121/1.411383 |
[13] | 刘金霞, 王克协, 曹正良, 等. 套管井应力集中诱导井周波速各向异性的研究. 地球物理学报 , 2005, 48(3): 717–723. Liu J X, Wang K X, Cao Z L, et al. The effects of stress concentrations of cased hole on induced azimuthal anisotropy velocity distribution. Chinese J. Geophys. (in Chinese) , 2005, 48(3): 717-723. |
[14] | 李刚, 王克协, 刘金霞, 等. 套管井应力集中下正交偶极声测井的声弹理论、数值分析与实例. 地球物理学报 , 2006, 49(1): 717–723. Li G, Wang K X, Liu J X, et al. Acoustoelasticity theory, numerical analysis and practical example of the cased hole with stress concentration by cross-dipole acoustic logging. Chinese J. Geophys. (in Chinese) , 2006, 49(1): 717-723. |
[15] | 李刚. 套管井正交偶极声测井异常地应力正、反演研究. 长春: 吉林大学, 2005. Li G. Acoustoelastic theory of guide waves in a cased hole and inversion of formation stress by cross-dipole acoustic logging (in Chinese). Changchun: Jilin University, 2005. |
[16] | 唐晓明, 郑传汉. 定量测井声学. 北京: 石油工业出版社, 2004 . Tang X M, Zheng C H. Quantitative Borehole Acoustic Methods (in Chinese). Beijing: Petroleum Industry Press, 2004 . |
[17] | Valero H P, IKegami T, Sinha B K, et al. Sonic dispersion curves identify TIV anisotropy in vertical wells. 79th Annual International Meeting, SEG, Expanded Abstracts , 2009: 322-325. |
[18] | Fritz I J. Third-order elastic constants for materials with transversely isotropic symmetry. Journal of Applied Physics , 1977, 48(2): 812-814. DOI:10.1063/1.323647 |
[19] | Sarkar D, Bakulin A, Kranz R. Anisotropic inversion of seismic data for stressed media: Theory and a physical modeling study on Berea Sandstone. Geophysics , 2003, 68(2): 690-704. DOI:10.1190/1.1567240 |
[20] | Prioul R, Bakulin A, Bakulin V. Nonlinear rock model for estimation of 3D subsurface stress in anisotropic formation: Theory and laboratory verification. Geophysics , 2004, 69(2): 415-425. DOI:10.1190/1.1707061 |