基于叠前岩石力学参数反演的三维地应力预测研究——以长宁区块页岩气为例

金文辉,井 翠,文 冉,宋修艳,陈珂磷,李雪松,张 奎

(1.四川长宁天然气开发有限责任公司,成都 610051;2.北京普瑞斯安能源科技有限公司,北京 100085)

0 引言

地应力是存在于地层中的天然应力,地表发生的一切板块运动、构造变形都是地应力作用的结果[1]。地应力的大小和方向是随着空间和时间变化的,因而是一个非均质的应力场。在油气田勘探开发过程中,地应力的大小和方向控制着油气的生成、运移和聚集,同时,天然裂缝和压裂缝的展布、地层破裂与坍塌、井网布置与优化、储层压裂设计等都与地应力密切相关[2],其研究工作是油气勘探开发的关键环节之一,同时也是沟通地质与工程的重要桥梁[3],在非常规油气勘探开发中更是如此。以页岩气为例,超低孔、特低渗的储层物性严重制约了天然气的流动性,因此需要通过人工干预手段将页岩储层改造为缝网系统,从而改善页岩的渗流能力[4]。解决上述一系列问题的前提是精确了解地下的地应力状态。目前,深部地层地应力大小和方向预测技术已经成为油气勘探开发领域的一大技术难题。

前人在早期的矿山工程建设中对地应力有了初步认识,同时发展了地应力测量技术。最早的原位地应力测量起始于上世纪30年代[5],至今已有近百年历史。目前,地应力研究方法多达二十余种。常用的方法有水压致裂法、岩石声发射法、测井资料估算法、差应变法、应力解除法等[6-11]。但这些方法多依赖于定点取芯试验或钻井测井资料,均为点数据,只能代表井点或井周有限范围内的地应力状态,不能有效反映井场或区域地应力分布特征,故需要探索新方法开展区域应力场分布特征研究。

近年来,国内外学者们就基于地震方法的地应力预测开展了探索性研究,该方法能够得到目标区连续的地应力剖面,反映三维空间内任一位置的地应力分布特征。Dilien[12]通过研究地震反射系数与地应力大小的关系,为实现地应力三维预测提供了可能;Tigrek[13]在此基础上开展了进一步研究,利用地震反射系数获得岩石力学动态参数,通过动静态转换,建立地质力学模型,从而获得地应力预测结果;Hunt[14]通过地震数据计算曲率和杨氏模量等,对地应力分布进行了预测;何英[15]利用构造应力场模拟,通过计算曲率预测区域构造应力场;宗兆云[16]利用地层裂缝的相关参数预测三维地应力;马妮[17]基于正交各向异性介质理论,提出了利用方位叠前地震数据预测地应力的方法;熊晓军[18]等通过叠前弹性参数反演,获得计算地应力分布的岩石力学参数。

笔者在结合前人研究成果的基础上,充分利用三维地震资料的优势,采用井震结合方法,利用构造解释得到的构造信息(层位与断层)和地震预测得到的地层信息(速度与密度)、地质力学信息(杨氏模量与泊松比),对地应力方向、最大和最小水平主应力大小及应力差异系数及进行了预测。通过开展长宁区块区域应力机制研究,并对3口水平井的微地震监测结果进行分析,从而验证了本文地应力预测结果的合理性,在页岩气水平井开发中应用效果良好,具有一定的推广应用意义。

1 区域地质背景

我国有较为丰富的页岩气资源,近年来页岩气藏勘探取得了重大进展。川南长宁地区五峰-龙马溪组深水陆棚相优质页岩分布稳定、保存条件好,早期开发的宁201井区页岩气勘探取得良好效果,页岩气勘探潜力大[19-22]。

长宁区块构造位于四川盘地与云贵高原结合部,川南古坳中隆低陡构造和滇黔北坳陷娄山褶皱带之间(图1)。长宁区块构造以南由娄山褶皱带控制,以北则主要受川南低褶带构造演化的影响,以东受四川盆地北东向边界延伸的影响,以西主要由华蓥山断裂带控制,整体表现出纵横错落的构造格局。研究区埋深适中,在2 100 m ～3 500 m。主力产气层位于下志留统龙马溪组和五峰组,可划分为五峰组、龙一11、龙一21、龙一31、龙一41共5个小层,主体厚度稳定分布在30 m～50 m,是中国重要的页岩气建产区之一[23-24]。

图1 长宁地区构造区划与研究区位置图Fig.1 Structural zoning and study area location map of Changning area

2 地应力预测

本次研究采用的三维地应力预测技术是在地层压力预测基础上,借助井筒岩石力学实验和水力压裂等数据进行质控分析,并优化了应力计算模型,具体研究路线见图2。

图2 地应力预测流程图Fig.2 Flow chart of ground stress prediction

1)三维地层压力预测:通过井筒实测的DIFT地层压力数据和测井资料,建立多种方法的地层压力预测模型,在此基础上优选出最佳的地层压力预测模型。关键步骤包括:①利用测井资料和速度场建立3D地层速度模型,通过时深转换,将叠前反演得到的时间域纵波速度、密度、杨氏模量、泊松比等成果数据转换至深度域;②利用测井资料拟合Amoco系数,补充浅层密度缺失段,拼接高精度密度反演数据体,求取密度积分,得到上覆地层压力;③基于有效应力理论,获取垂直有效应力与地层速度(高精度速度反演体)之间的关系,通过对比井上实测压力资料,对地层压力模型参数进行优化,预测孔隙流体压力和地层压力系数。这一环节的目标是获得上覆地层压力Sv、孔隙压力Pp,从而为计算由上覆地层泊松效应产生的地应力提供基础数据。

2)构造应力场模拟:利用构造解释层位获取地层的几何信息,利用储层反演获取地层岩性信息(包括厚度、速度和密度),基于广义胡克定律和薄板理论,求得层位面的地应力参数,包括地层面的曲率张量、变形张量和应力场张量,最终得到主曲率、主应变和主应力,进而获取目的层应力场展布特征[25-26]。这一环节的目标是预测地应力方向,同时获取沿着最大、最小水平主应变,从而计算构造应力及残余应力。地应力大小最终结果来自泊松效应、构造应力和残余应力的总和。

3)地应力计算模型方法优选与参数拟合:结合区域构造特征与演化,结合实测地应力数据(包含岩心综合实验、水压致裂法等地应力测量结果),分析研究区的区域地应力机制,选取适合本区的地应力大小预测模型,并利用实测地应力数据刻度模型中的未知参数,比如奥特系数、构造应力系数等。利用地质统计学反演得到的高精度岩石力学参数泊松比、杨氏模量数据,结合三维地层压力预测得到的上覆地层压力Sv、孔隙压力Pp和构造应力场模拟得到的最大、最小水平主应变成果,计算本区三维地应力大小。

2.1 地应力方向预测

应力场研究,即研究地应力分布的规律性。地应力方向预测的理论基础是构造力学[27-30],应力场模拟需要较为理想的模型:假设地层介质是均匀连续、各向同性且完全弹性的,地层的构造变形完全由地应力所致。通常发生弯曲变形地层的宽度和长度都远远大于其厚度,因此用薄板弯曲理论可以模拟构造面附近的应力状态[31-33],前人基于薄板理论预测应力方向的研究与应用较多,笔者不再赘述其原理方法。

本次研究共收集到长宁区块6口探井的FMI成像测井数据,统计研究区地应力方向数据发现:研究区地应力方向在不同区域内变化较大,NW向、NWW向、近EW向均有分布,反映出本区地应力的复杂多变性(表1)。

表1 长宁区块最大水平主应力方向预测符合情况统计Tab.1 Statistical table for prediction of maximum horizontal principal stress direction in Changning Block

基于构造力学的应力场数值模拟理论中指出,按构造位置的不同,构造缝发育程度与方向存在较大的差异性,具体表现为:褶皱顶部受拉张应力,底部受挤压应力,顶部(浅层)拉张缝较为发育,底部(深层)构造缝不发育。在外部挤压应力作用下,背斜顶部拉张缝发育,方向垂直于最大主应力;向斜构造缝相对不发育,方向平行于最大主应力[34-37]。

应用应力场数值模拟方法,得到了目的层地应力方向分布特征和构造缝分布特征,见图3。构造应力场模拟结果显示:长宁区块目的层五峰组～龙马溪组在外部挤压应力作用下,最大水平主应力方向与井筒所在构造位置(图4)及是否位于断层附近有关:以D井区为代表的背斜(隆起区)最大主应力方向与构造缝方向垂直。以A井、B井、F井区为代表的向斜(低洼区)最大主应力多与构造缝方向平行,靠近断裂区域最大水平主应力多与构造缝方向垂直。

图3 五峰组之上30 m构造缝及最大水平主应力方向叠合图Fig.3 Overlapping map of 30 m structural joints and maximum horizontal principal stress direction above Wufeng Formation

图4 五峰组底面构造特征Fig.4 Structural characteristics of the bottom of Wufeng Formation

同时,应力场数值模拟还可反映应力环境的相对强弱,图3中蓝-深蓝色区域代表拉张应力环境,应力强度较小,地应力非均质性较弱,对井筒压裂的影响相对较小;绿-黄-红-粉红区域代表挤压应力环境,应力强度随着色标由冷变暖而逐渐变大,地应力对井筒压裂的影响也逐渐增强,在水平井轨迹设计时,应重点考虑粉红色强应力环境对井筒的影响,尽量避开强应力环境区域。

研究区地应力预测方向与FMI成像测井识别方向吻合情况见表1,研究区6口探井的地应力预测方向基本与井筒FMI识别的最大水平主应力方向一致。研究区各井地应力方向多以NWW、NW向为主,个别井位受局部多组应力共同作用影响,地应力方向呈现多簇特征。地应力方向预测符合率整体较高,证明该方法预测地应力方向具较高的可靠性。

2.2 地应力大小预测

2.2.1 地应力大小预测原理及方法

地应力预测是在一定的假设条件下,利用地应力实测数据建立相对简单的地应力计算模型,并结合地球物理和测井数据进行地应力计算分析的一种方法,所建立的计算模型在一定程度上直接影响计算结果。目前业内计算地应力的模型主要有:①基于地层处于剪切破坏临界状态这一假设,结合最大主应力、最小主应力之间的关系提出的Mohr-Columb模式。②单轴应变模式,具有代表性的计算模型有:Matthews &Kelly模型(1967)、Anderson模型(1973)等。③黄荣樽模式(1984):该模式可以解释水平应力大于垂向应力的现象,同时也考虑了构造应力的影响;④组合弹簧模式(1988年),又称斯伦贝谢模式,在其开发的地应力预测商业软件中广泛应用;⑤多孔弹性水平应变模型[28-30]。

近年来业内应用较多的是组合弹簧模型和多孔弹性水平应变模型。二者综合考虑了岩石力学、孔隙流体压力及构造作用等因素对地应力的影响。模型首先假设岩石为各向同性的均质体,且是满足线弹性的连续地层,并假设在沉积及后期地质构造运动过程中地层和地层之间无相对位移,地层两水平方向的应变为常数。本次研究采用的地应力大小预测方法主要依据这两种理论进行融合与优化,共考虑了3种效应对现今地应力大小的影响,分别是:泊松效应、水平构造应力和残余应力。其中泊松效应是指上覆地层受地心引力作用对下伏地层产生的应力,又称重力应力;水平构造应力是指原岩中存在的促使水平方向发生构造运动的应力;残余应力是古构造应力的残余场,其作用释放缓慢,对现今地层仍会产生一定应变影响[31]。与前人采用的地应力模型相比[33-38],笔者增加了残余应变对现今地应力的影响,在计算公式中增设了构造形变和残余应力影响因子,可提升受多期构造运动影响较大的研究区块地应力预测结果的精确度。本次研究采用的预测模型其计算模型如下:

图5公式中σh为最小水平应力,σH为最大水平应力,σv为上覆地层压力,α为Biot系数,Pp为孔隙压力,v为泊松比,E为杨氏模量,εh为水平最小主应变,εH为水平最大主应变,A和B分别为最小、最大构造应力系数,a和b分别为最小、最大构造应力系数,D为深度值。

图5 采用的地应力大小计算公式Fig.5 The formula for calculating the magnitude of in-situ stress is adopted in this paper

式中α为多孔弹性系数,又称比奥特(Boit)系数,一般高渗岩心为0.9～1;中渗岩芯为0.8～0.9;低渗为0.6～0.75。一般满足:孔隙度<α<1,同样孔隙度的灰岩,其比奥特系数明显大于同孔隙度的砂岩。骨架结构越稳定,比奥特系数越小。

式中A、B分别为最小和最大主应力方向的应力系数,a、b分别为残余最小和最大应变系数,均需要钻井岩心实验得到的水平应力数据刻度,而水平应力测量数据的常见获取方式分为两种:一种来自水力压裂试验,该方法可获得最小水平主应力,再根据胡克定律和Zoback应力多边形,可推导出最大水平主应力计算模型,也需要工程或实验室数据进行刻度;另一种为声发射实验,利用Kaiser方法可测定最小和最大水平主应力。

图6中展示了利用前人地应力模型[33-38],即只考虑泊松效应和水平方向构造应力效应的预测结果,和笔者采用的考虑残余应变对现今地应力的影响后地应力预测结果。二者对比可知,后者在全区横向地应力变化特征的刻画上精度更高,更能精确地反映地应力大小的局部变化特征,同时可进一步提高预测结果的符合率,笔者采用的地应力大小预测模型符合率较前人使用模型符合率提高约3%(表2)。

表2 长宁区块考虑残余应力前、后最小水平主应力大小预测符合率统计表Tab.2 Statistical table of prediction coincidence rate of minimum horizontal principal stress before and after considering the residual stress in Changning Block

图6 增设残余应力前后最小水平主应力预测结果对比Fig.6 The comparison of prediction results of minimum horizontal principal stress model before and after adding residual stress

2.2.2 地应力大小预测结果

针对区域应力机制研究,笔者采用了基于Zoback应力多边形理论的Anderson模型(如图5(a)):当垂向应力最大,最小水平主应力足够小的时候(NF:Shmin

本次研究对长宁区块岩心综合实验获取的10口井Kaiser实验数据(包含邻区4口井,共计25个样点)进行对比分析,结果表明:长宁区块五峰组～龙马溪组最大水平主应力、最小水平主应力、垂向应力的大小排序关系为:SHmax>SV>Shmin,即上覆地层压力大小处于最大、最小水平主应力之间,为典型的走滑断层模式(图7(b))。

图7 Anderson三种应力环境断层模式图及长宁区块应用结果Fig.7 Anderson three kinds of stress environment fault pattern diagram and Changning Block application results

地应力模型以研究区内6口探井实测地应力数据(Kaiser实验)为依据,精细刻度了龙马溪组～五峰组各个地应力计算模型参数,即公式(1)、(2)中的比奥特系数α、构造应力系数a和b、局部应变系数A和B,最终得到最小、最大水平主应力预测结果(图8、图9)。

图8 长宁区块五峰组之上30 m最小水平主应力预测图Fig.8 Prediction of minimum horizontal principal stress at 30 m above Wufeng Formation in Changning Block

图9 长宁区块五峰组之上30 m最大水平主应力预测图

五峰组之上20 m基本包含了长宁区块页岩气主力产层:五峰组、龙一11和龙一21小层。由预测结果可知:平面上最小水平主应力系数分布范围为57 MPa～83 MPa,最大水平主应力系数分布范围为67 MPa～99 MPa,总体表现为西北和中东部高、中部偏西较低的特征。整体受埋深控制非常明显(图10):两个西北(E井区)和中东部(B、F、C井区以南)向斜深埋区水平主应力较高,两个高值区中间所夹的中西部条带(D、A井区之间)浅埋区水平主应力较低。

图10 长宁区块五峰组底面埋深图

研究区地应力大小与岩心实测地应力数据(Kaiser实验)吻合情况见表3,最大、最小水平主应力的预测符合率均大于90%,证明该方法预测地应力大小具较高的可靠性。

表3 长宁区块地应力大小预测符合率统计表Tab.3 Statistical table of prediction coincidence rate of in-situ stress in Changning Block

为了更好地表征页岩地层压裂改造效果,引入应力差异系数的概念:所谓应力差异系数是指最大与最小水平主应力差值与最小水平主应力的比值。该参数越小,表明地层压裂后形成网状缝的能力越强;该参数越大,表明地层压裂后形成的压裂缝方向较为单一,储层改造效果差。利用预测得到的三维最大、最小水平主应力数据体计算应力差异系数,由图11可知:应力差异系数整体较低,多在0.15～0.2之间,大部分地区的压裂造缝能力较好。

图11 长宁区块五峰组之上30 m应力差异系数预测图Fig.11 Prediction of stress difference coefficient at 30 m above Wufeng Formation in Changning Block

此外,结合同期裂缝预测研究成果对优质页岩的可压性做进一步分析。图12给出了多尺度裂缝预测图,图中红色代表大尺度裂缝,断距>20 m,为构造解释获得的断层结果;蓝色代表中尺度裂缝,8 m≤断距≤20 m,为利用地震资料进行机器学习得到的裂缝分布概率体;绿色代表小尺度裂缝,断距<8 m,为利用地震资料计算的地层曲率体,往往与大中尺度裂缝相伴生,为大中尺度裂缝分支或分布于大中尺度裂缝周缘地区。结合图11可知:应力差异系数的分布受断裂及裂缝的控制较为明显,表现为:断裂(大、中尺度裂缝)发育区域内,应力差异系数为高值(图中蓝色区域),天然裂缝发育,方向较为单一,不易产生复杂压裂网,压裂改造效果受限;小尺度裂缝发育区的应力差异系数为中高值(图中黄、绿色区域),压裂可激活天然裂缝,且形成复杂缝网的能力适中,压裂改造效果中等;多尺度裂缝欠发育区,水平应力差和应力差异系数低值(图中红色区域),弱应力区域有利于复杂压裂缝网的形成,压裂改造效果较好。

图12 长宁区块多尺度断裂预测平面图Fig.12 Multi-scale fracture prediction plane map of Changning Block

3 地应力预测应用

随着石油工业的发展,非常规油气资源逐渐成为各国油气田的增产主体。页岩气藏在全球范围有着广泛的分布,与常规油气藏相比,页岩气藏有着自生自储、低孔低渗的特点,其实现大规模开采的关键在于水平井体积压裂改造技术的应用。页岩气水平井开发,需要对目的层进行压裂、造缝,能否将页岩储层压裂成网状缝是导致水平井能否增产的重要因素。

选取研究区内已压裂投产,且具有微地震资料的3口水平井:X1井、X2井、X3井,结合预测的地应力结果,对主要产气层的可压性做定性分析。图13中分别展示了X3井井周边天然裂缝发育特征、地应力方向分布特征、应力差异系数分布特征和微地震监测结果:该X3井周边发育多条过井天然裂缝,且天然裂缝发育分布较为均衡,井筒东西两侧、北段和南段不存在明显的差异性。

图13 X3井压裂效果分析Fig.13 X3 well fracturing effect analysis

地应力环境整体为拉张应力,地应力非均质性较弱,弱应力环境对压裂造缝较为有利。地应力方向总体与井轨迹方向垂直或高角度斜交,井轨迹设计较为合理,垂直于井筒方向地层易于压开。综上所述,天然裂缝发育程度、地应力环境均不是影响最终压裂效果的主要因素。由微地震检测结果可知:该井微地震事件在靠近北段数目更多,能量更大,而在靠近水平井南段和中段的微地震事件数目更少,能量较小,反映出水平段压裂效果的差异性。由应力差异系数预测结果可知,该井在跟端及水平井中段的应力差异系数较大(图13c中蓝色区域),高于0.2,不易产生网状压裂缝,这与微地震监测结果吻合,进一步印证了地应力预测结果的可靠性。

应用同样的研究方法,对X2、X1井压裂效果进行评价。X2井周边发育多条过井天然裂缝,水平段北侧和南侧天然裂缝发育程度略高,中段天然裂缝欠发育;地应力环境整体为拉张应力,地应力非均质性较弱,弱应力环境对压裂造缝较为有利;且地应力方向总体与井轨迹方向垂直或高角度斜交,井轨迹设计较为合理,垂直于井筒方向地层易于压开;应力差异系数在中段较高(图14c中蓝色区域),高于0.2,不易产生网状压裂缝,而在北侧和南侧应力差异系数较低,易于产生网状缝;该井微地震识别的SRV在北侧和南侧较高,而在水平井中段较低,反映出水平段压裂效果的差异性。综上所述,X2井和X3井的压裂效果受天然裂缝和地应力双因素控制。

图14 X2井压裂效果分析Fig.14 X2 well fracturing effect analysis

图15分别展示了X1井周边天然裂缝发育特征、地应力方向分布特征、应力差异系数分布特征和微地震监测结果:该井周边天然裂缝发育较少,水平段北侧天然裂缝发育程度略高于南侧,中段天然裂缝欠发育;地应力环境整体为拉张应力,地应力非均质性较弱,弱应力环境对压裂造缝较为有利;地应力方向总体与井轨迹方向中等角度斜交,考虑到水平段整体应力环境较弱,压裂缝仍然大概率沿垂直于井筒方向扩展;应力差异系数在靠近北侧和中段较高(图15c中蓝色区域),高于0.2,不易产生网状压裂缝。而在水平段南侧,应力差异系数较低,压易于产生网状缝;该井微地震监测结果显示:在水平段南侧微地震事件更多,局部成团状;而在中段和北段微地震事件相对南侧减少,且呈现明显的单一方向展布,与应力差异系数预测结果吻合。综上所述,该井压裂效果受天然裂缝和地应力双因素控制。

图15 X1井压裂效果分析Fig.15 X1 well fracturing effect analysis

4 结论

笔者采用三维地应力预测技术成功在长宁区块实际资料中进行了应用。长宁区块最大/最小水平主应力均受埋深影响最大,呈现深埋区高应力,浅埋区低应力的特征。全区应力差异系数均受埋深(构造)和断裂的双重影响:深埋区、向斜中心部位及靠近断裂区域内应力差异系数较大,地应力非均质性较高;浅埋区、向斜斜坡与断隆部位及远离断裂部位应力差异系数较小,地应力非均质性较低,易于压裂改造,且应力差异系数越小,越有利于形成压裂缝网,有利于水平井压裂改造。