渝东南彭水地区页岩储层构造应力场模拟解析

葛 勋，汤济广 ，赵培荣，唐 永，许启鲁

1.页岩油气富集机理与高效开发全国重点实验室，北京 昌平 102206；2.中国石化页岩油气勘探开发重点实验室，北京 昌平 102206；3.中国石化石油勘探开发研究院，北京 昌平 102206；4.长江大学地球科学学院，湖北 武汉 430100；5.中国石油化工集团，北京 朝阳 100728；6.长江大学非常规油气湖北省协同创新中心，湖北 武汉 430100；7.中国石油大学（华东）地球科学与技术学院，山东 青岛 266580

引言

渝东南彭水地区在奥陶统五峰组—下志留统龙马溪组发育一套薄层状灰黑色—黑色炭质页岩夹硅质页岩，页岩气成藏条件优越[1-3]。彭水区块桑柘坪向斜中彭页1 井、彭页3 井和彭页4 井已完成工程压裂，日产量具有工业产能[4-5]。

中燕山期—喜马拉雅期，研究区地层发生大幅度的挤压变形并伴随隆升剥蚀[6-10]，构造改造强度与持续时间对页岩储层含气性丰度具有明显的调整作用，其主要原因是古应力场决定了天然缝网的发育程度，增强了渗流能力，可以有效地改善页岩的储集性能[11-13]，从页岩气井高产区分布来看，构造抬升时间越晚，对页岩气成藏越有利[14-15]。通过对渝东地区节理、断裂及滑动特征进行详细观测，可以厘定构造活动期次、古应力场演化等构造特征[16-17]，采用磷灰石裂变径迹空间分布特征定性分析的方法，确定鄂西渝东地区从97 Ma 开始持续抬升剥蚀[18-19]。通过研究直立褶皱、共轭节理和断裂面擦痕等各种构造行迹，可以重建一个研究区的构造应力场的演化特征，分析其构造活动规律[20-21]，也可以通过钻井诱导缝走向与井壁坍塌方位分别确定古、今构造应力场的方向，采用测井资料计算古、今构造应力场的大小[22]。在确定了古、今应力场的方向和大小后，主要应用Abaqus 和Ansys 有限元分析软件建立三维地质模型，确定边界条件，进行古、今构造应力场数值模拟，进而深入分析古、今应力场的分布规律以及构造裂缝的发育程度[23-26]。

前人对于构造应力场的研究大部分通过节理、褶皱等野外构造行迹对某一地区进行构造演化重建，或通过磷灰石裂变径迹热史模拟升降过程，而针对某一特定研究区页岩储层应力场二维数值模拟的研究较少。因此，本次笔者利用野外构造行迹对彭水地区进行构造应力场演化分析，再应用Abaqus进行构造应力场模拟，较为详尽地梳理了构造应力场模拟的应用流程，目的是揭示研究区内应力展布规律及其区域构造发育的控制效应，进而深入叙述区域构造应力场的性质、大小。构造应力场数值模拟不仅能够定量表征彭水桑柘坪向斜的地应力，也对整个渝东南地区页岩气品质和勘探开采难度的评价具有重要的意义。

1 区域地质概况

彭水地区桑柘坪向斜位于重庆东南部（简称渝东南，以下均同），川东南—湘鄂西褶皱带的过渡地区，背斜向斜交替发育，褶皱整体呈NE 向展布，大地构造上属于中扬子台内拗陷构造单元[27-28]，如图1所示。

图1 渝东南彭水地区位置和构造Fig.1 Tectonic location and geological structural of Pengshui Area in the southeastern Chongqing

渝东南主要经历了加里东期、印支期、燕山期和喜马拉雅期等多期运动，加里东期和印支期的构造运动以隆升为主，燕山期和喜马拉雅期对彭水地区的影响最为显著，褶皱和断裂主要形成在中燕山期，褶皱和断裂走向多为NE 向[29-31]，喜马拉雅期，由于青藏高原的隆升引起近EW 向的挤压，研究区内发育NW 向断裂，见图1c。

渝东南地区经历了多期的褶皱隆升作用，广泛缺失泥盆系和石炭系，局部缺失侏罗系。五峰组—龙马溪组在背斜核部广泛出露，页岩气保存条件遭受破坏，不利于气藏保存。

该地区向斜褶皱则较为宽缓，地层剥蚀厚度较小，页岩埋藏适中，天然裂缝呈网状发育，脆性指数高，成为页岩气井的高产区[14]。根据彭页1 井综合柱状图特征，可知在五峰组—龙马溪组底部地层页岩具有很好的含气性显示，TOC 在3%∼4%，Ro普遍大于2%，有效泥页岩层段厚度约90 m（图2）。

图2 彭页1 井综合柱状图Fig.2 Comprehensive column chart of Well PY1

2 古构造应力方向解析

2.1 节理对古应力方向的指示

为探讨彭水地区的构造演化史，需要对节理进行分期与配套，从时间、空间和形成力学上研究一个地区节理的分布产出规律，可以为该地区的构造应力场模拟提供可靠依据。若两组节理互切，表明两组节理同时形成一对共轭剪节理，利用共轭剪节理系统地解析印支期以来古构造应力场的方向和期次，重现渝东南彭水地区及周缘构造演化特征与动力学机制[32]，对该区页岩气保存条件差异性的研究具有重要的指导意义。

通过一定数量的共轭剪节理统计与分析，可以确定区域应力场的古构造方向[33-34]。对研究区进行野外地质详查，共测量了510 个节理的产状数据，这些节理的观测点广泛分布于渝东南彭水地区及周缘，出露的五峰组—龙马溪组有节理发育。为研究方便，对地层进行了拉平处理（在岩层近水平地区进行节理分析时，岩层水平校正可以省略；在岩层倾斜地区进行分析时，要进行倾斜校正）。通过地层拉平处理、统计分析发现，整个研究区发育的节理，主要分为以下5 组（图3）。

图3 彭水地区及周缘节理走向分析Fig.3 Analysis of the strikes of fractures in Pengshui Area and its periphery

1）NNE 向节理（10◦±10◦）。根据统计的玫瑰花图来看，在研究区整个范围内均有发育。砂岩中节理未被填充，灰岩及灰质云岩中有区域可见方解石填充，整体上几何形态较为平直，发育较为稳定，近与层面垂直。

2）NE 向节理（50◦∼70◦）。此方向节理统计点在全区的二叠系—三叠系中均有发育，节理面相对粗糙，延伸距离较短，节理在砂岩中较为发育，泥岩中趋于消失，分析其性质为张节理。

3）近EW 向节理（100◦±10◦）。主要发育在桑柘坪、道真，在节理面上偶见擦痕，说明节理具有左旋走滑的特点。

4）NW 向节理（135◦±10◦）。在全区的侏罗系—三叠系中均有发育。

5）NNW 向节理（340◦±10◦）。在全区所有地层中均有发育，说明该组节理的发育时间相对较晚。

通过野外露头观察剪节理的菱形结环、折尾、分叉以及交切关系[34]，可以发现NNE、EW、NNW和NW 向节理面平直，节理面还有较为明显的擦痕，结合节理之间的交切关系，可以判断NNE 向与EW 向组成一对平面X 型共轭剪节理（图4a，图4b）；NNW 向与NW 向组成一对平面X 型共轭剪节理（图4c，图4d）。NE 向节理延伸范围有限，且节理面相对粗糙，为张节理（图4e）。

图4 彭水地区及周缘五峰组--龙马溪组野外露头照片Fig.4 Field outcrop of Wufeng-Longmaxi Formations in Pengshui Area and its periphery

对节理形态、力学性质和发育程度进行细致的观测，并对地层复平转化，可以求取研究区共轭剪节理的优势产状。同期配套的节理与其所受的主应力方向存在一定的几何关系：一对共轭剪节理的交线平行于中等主应力方向；它们的夹角平分线分别为最大主应力轴和最小主应力轴。一般情况下，共轭剪节理的锐角等分线平行于最大主应力方向；共轭剪节理的钝角等分线平行于最小主应力方向，因此，可以对本区主要的共轭剪节理优势产状及其应力状态进行分析（表1）。

表1 节理产状统计及应力场分析（水平恢复后）Tab.1 Joints occurrence statistics and analysis of stress field（after level recovery）

基于表1 中主应力的分析可知，主要具有两个方向的古构造应力。

1）NW 向褶皱带的主应力方向为32◦∼81◦，倾伏角2◦∼29◦。

2）NE 向褶皱带中等主应力方向为120◦∼151◦，倾伏角2◦∼11◦。武隆县城东、南川和道真旧城最大主应力方向为213◦∼319◦，最大主应力逐渐向NNW方向转变，说明研究区受SE—NW 向挤压扩展的过程中受到一定的阻挡呈现出较为明显的旋转特征。整体可以反映出NE—SW 向和NW—SE 向两期挤压作用。

2.2 断层擦痕对古构造应力方向的指示

断层擦痕是断层两盘相对错动在断层面上因摩擦作用留下的痕迹，根据断面法线和擦痕确定的平面可以确定主应力方向。研究区内断裂较为发育，断面多表现为挤压特征，并在结构面以及滑抹晶体上表现出多期擦痕，这些擦痕是断层多期活动的直接证据（图4f）。通过对研究区内大量断层擦痕的测量，可以恢复不同时期构造应力的主要方向[32]。

武隆、桑柘坪、道真向斜的擦痕大圆圈主要反映出两个方向的古构造应力见图5（其中，σ1—最大主应力，MPa；σ2—中等主应力，MPa；σ3—最小主应力，MPa）和表2。

表2 断层擦痕统计及应力场分析Tab.2 Fault striations statistics and analysis of stress field

图5 彭水地区及周缘断层擦痕构造应力解析Fig.5 Analysis of the fault striations tectonic stress in Pengshui Area and its periphery

1）SE—NW 向的挤压，最大主应力的方向为125◦∼160◦，倾伏角30◦∼50◦。

2）NE—SW 向的挤压，最大主应力方向为58◦∼72◦，倾伏角33◦∼36◦。

通过断层擦痕恢复出的古应力方向与共轭剪节理恢复出的方向一致，证明此结果具有较好的可靠性。

3 古应力场构造演化解析

3.1 古应力场构造演化时期

3.1.1 流体包裹体均一温度分析

流体包裹体是矿物结晶过程中捕获的流体，它具有固流体当时的温度、压力及成分等方面的信息。如果包裹体的宿主矿物是构造作用中构造带的充填矿物，据此可以基本判断构造作用的时间期次[35-36]。

通过对彭水地区及周缘断裂（裂缝、断裂岩等）填充矿物的流体包裹体进行测试（图6），获得了包裹体形成的均一温度。对研究区测得的盐水包裹体均一温度数据分析显示，研究区包裹体的均一温度分布在149.6∼119.9◦C和96.0∼82.5◦C两个区间，见表3（其中，包裹体宿主矿物为方解石，原生成因，主要成分为盐水），说明在渝东南彭水地区曾经历了两期构造演化阶段。

表3 彭水地区及周缘包裹体均一温度Tab.3 Homogenization temperature of inclusions in Pengshui Area and its periphery

图6 彭水地区纯气相包裹体与纯液相包裹体Fig.6 Pure gas inclusions and liquid inclusions in Pengshui Area

3.1.2 磷灰石裂变径迹证据

研究区的主体构造方向以NE 和NW 向为主，其形成与演化经历了两次构造运动。根据渝东南地区褶皱带的地层接触关系，中下侏罗统与下伏地层（中下三叠统间）呈角度不整合，上白垩统与下伏地层间呈高角度不整合。基于前人磷灰石裂变径迹分析结果的热史模拟（图7），4 个样品分别取自方斗山（SJ-2 井）、万县（WE-8 井）、利川（X5 井）和石柱（JL-1井），样品均在中燕山期达到最大埋深，之后迅速发生隆升剥蚀，晚燕山期—喜马拉雅期再次抬升至地表接受隆升剥蚀。热史模拟结果显示，渝东南地区在燕山期—喜马拉雅期至少经历两次强烈的构造运动[37-38]。

图7 磷灰石裂变径迹的热史模拟[38]Fig.7 Thermal history simulation of apatite fission track

3.2 古应力场构造演化序列

古应力场构造演化序列是根据各构造活动阶段所形成的构造行迹的交切关系来判定不同构造演化阶段的先后顺序[6]。研究区主要发育NE 向、近NS向和NEE 向褶皱体系，基于褶皱的叠加特征来判断渝东南彭水地区及周缘的构造演化序列（图8）。

图8 彭水地区及周缘褶皱叠加特征Fig.8 Fold superposition characteristics in Pengshui Area and its periphery

马武背斜和桑柘坪向斜（图8a，图8b）两期水平褶皱的褶皱轴呈大角度相交或直交叠加。NE 向褶皱长于NW 向褶皱且被NW 向褶皱截切，说明NW向构造晚于NE 向构造形成；洛龙背斜（图8c）由两期直立水平褶皱叠加改造形成，相互连接呈“L”型。说明NW 向构造限制了近NNE 向构造的发展，NW向构造晚于近NNE 向构造形成；利川向斜、南川向斜和焦石坝背斜（图8d，图8e，图8f）由两期直立水平褶皱叠加改造形成，相互连接呈“T”型，早期形成的NE 向被晚期形成的NW 向构造改造，说明NW向构造晚于近NE 向构造形成。

综合以上对研究区内褶皱叠加特征的分析，可以确定彭水地区及周缘残留向斜的构造演化序列，NE 向构造早于NW 向构造形成。其整体的构造演化特征为中燕山期，江南隆起由SE--NW 的递进扩展作用，致使渝东南地区经历了强烈的构造变形作用，发育了大量NE 向和NNE 向褶皱与断裂[39]；晚燕山期—喜马拉雅期，印—亚板块的碰撞引起青藏高原的隆升，使研究区处于SW—NE 的挤压应力场而形成NW 向褶皱[40]。

4 多期古构造应力场模拟方法

采用有限元的基本思想，在明确构造应力期次、大小、方向及演化时序的前提下，对构造应力场进行数值模拟，能够量化分析构造应力场的大小及展布规律[41]。本文借助Abaqus 软件精准建立了地质模型和力学模型，确定了边界条件，最终对渝东南彭水地区中燕山期—喜马拉雅期构造应力大小进行了定量数值模拟[42-43]。

4.1 地质模型

彭水地区页岩气勘探目的层为五峰组—龙马溪组，所以只需要对这个层位进行二维构造应力场模拟即可。将二维地质模型划分为两个力学区域：页岩区（主要为黑色页岩、硅质页岩和粉砂质泥岩）和断裂带。根据彭水地区断裂的走向以及构造演化特征将构造应力场分为两个阶段进行模拟，中燕山期发育NE 向断裂，晚燕山期—喜马拉雅期发育NW 向断裂。

4.2 力学模型

地质模型仅展示了彭水地区地层以及断层的空间形态，需要对地质模型进行材料力学参数的填充变为力学模型后，才能够参与计算。在渝东南地区选取了8 个岩石样品进行力学参数分析，三轴实验样品被加工成ϕ25 mm×50 mm 的圆柱形试样，使岩样的长径比≥1.5，将圆柱试样两端磨平，两端面的不平行度小于0.015 mm。三轴岩石力学实验的围压设定为30 MPa，通过三轴压缩实验，得到了岩石样品的杨氏弹性模量和泊松比，再通过杨氏弹性模量与泊松比计算抗压强度、内聚力及内摩擦角，最终采用加权平均法和类比法得到实验结果。

为了使模拟结果更好地反映出真实的地应力状态，将发生强烈形变的断裂带视为软弱带。跟页岩区相比，断裂带岩体的密度没有任何改变，杨氏弹性模量相对较小，约为页岩区的50%，泊松比略高于页岩区，两者之间差0.02∼0.10。依据CPS8R 单元，将地质模型划分为多节点组成的多类型单元，包括节点14 019 个，网格单元14 162 个，见图9，模型力学参数见表4。

表4 五峰组—龙马溪组物理模型力学参数Tab.4 Mechanical parameters of physical model of Wufeng-Longmaxi Formations

图9 网格划分Fig.9 Grid mesh

4.3 边界条件

有限元方法中的边界条件由应力方向、应力大小以及边界约束3 个方面组成[23]。利用岩石声发射Kaiser 试验测量了渝东南地区多期构造应力场的强度值。根据区域构造演化分析，渝东南地区主要发生了两次大规模的构造运动，分别发生在中燕山期和晚燕山期—喜马拉雅期。因此，将其分为两个时期进行模拟，中燕山期构造变形作用最强烈，所受最大主应力为79.62 MPa，其次为晚燕山期—喜马拉雅期，所受最大主应力为51.85 MPa（表5）。

表5 利用岩石声发射测量最大主应力值Tab.5 The result of maximum principle stress measured value from the acoustic emissions

在中燕山期，整个彭水地区处于SE—NW 向的挤压应力环境中，对模型的NE、NW 和SW 等3 个方向的边界进行约束，对其SE 方向施加79.62 MPa的梯度作用力；在晚燕山期—喜马拉雅期，对模型的NW、NE 和SE 等3 个方向的边界进行约束，对其SW 方向施加51.85 MPa 的梯度作用力，见图10。

图10 应力模拟边界条件Fig.10 Boundary condition for the stress simulation

4.4 计算输出模型

建立了力学模型之后，应用Abaqus 对彭水地区的五峰组—龙马溪组进行二维平面古应力场模拟。基于有限元的基本思想，将弹性连续体离散成为有限个单元；针对每个单元，通过固体力学变分原理，将各个单元刚度矩阵集合起来，计算求解后便可得到岩石储集层各个单元节点的变形位移。

计算方程主要集中在3 个方程：静力学平衡方程、几何学几何方程和材料力学本构方程。计算节点位移的方程是由总体刚度矩阵所表示的线性代数方程组，由位移场和几何方程计算获得形变场，最后由形变场和本构方程计算应力场。

4.4.1 平衡方程

在外力作用下处于平衡状态时，其内部各单元节点处的应力状态不同。根据理论力学，在外力作用下平衡的物体，其表面各处的应力分量应当与作用在该处的外力分量相平衡，平衡方程将边界静力条件与施加在目标体单元上的应力建立起联系

式中：σx—x向主应力，MPa；

τxy—xy面上的剪切应力，MPa；

τxz—xz面上的剪切应力，MPa；

Fx—沿x向的外力分量，N；

σy—y 向主应力，MPa；

τyz—yz面上的剪切应力，MPa；

Fy—沿y 向的外力分量，N；

σz z—向主应力，MPa；

Fz—沿z向的外力分量，N。

本次研究的是二维平面，有关z方向上的应力和应变不予考虑。

4.4.2 几何方程

当目标岩体受到外界作用力的影响，构造岩体内部各单元点的位置会发生变化。因此，研究岩层在外力作用下的变形规律，只需要研究物体内各点的相对位置变动情况。由于刚体位移不产生变形，故仅讨论变形位移，若岩体在应力作用下变形，其变形位移可以转化为x、y 和z方向的3 个分量u、v和w，则岩体的线应变和剪应变为

式中：εx—x向线应变，无因次；

u—x向位移量，m；

εy—y 向线应变，无因次；

v—y 向位移量，m；

εz—z向线应变，无因次；

w—z向位移量，m；

γxy—xy 面上的剪应变，无因次；

γyz—yz面上的剪应变，无因次；

γxz—xz面上的剪应变，无因次。

4.4.3 本构方程

对某种具体的岩石储集体，在一定围压作用下，应力和应变之间具有一定的量化关系，本次研究过程中将目标岩层在一定范围内视为各向同性体，在此前提下，应力-应变的关系为

式中：E—杨氏弹性模量，MPa；

µ—泊松比，无因次；

G—剪切模量，MPa。

模拟结果显示，中燕山期，彭水地区五峰组—龙马溪组最大主应力集中在22∼129 MPa，最小主应力集中在9∼99 MPa，见图11。

图11 中燕山期五峰组—龙马溪组主应力分布图Fig.11 Distribution of principle stress of Wufeng-Longmaxi Formations in the Middle Yanshanian period

由图11 可见，最大主应力的高值区主要分布在向斜的核部，最大主应力从核部到两翼逐渐降低，断裂带末端和断层的交汇处也处于高应力状态，应力集中在117∼129 MPa。此模拟结果与实际地质状况相符，向斜的转折端裂缝和节理发育程度高，岩石变形严重，断裂带末端岩石破碎程度严重，处于沿断层走线逐渐伸展的状态。

晚燕山期—喜马拉雅期，最大主应力普遍集中在13∼78 MPa，最小主应力集中在4∼45 MPa（图12）。此时，其最大主应力与最小主应力的值均低于中燕山期，与实际地质状况相符，中燕山期构造变形作用最强烈。

图12 晚燕山期—喜马拉雅期五峰组—龙马溪组主应力分布图Fig.12 Distribution of principle stress of Wufeng-Longmaxi Formations in the Late Yanshanian-Himalayan period

桑木场向斜中目前完成5 口页岩井钻探，彭页1 井、彭页3 井和彭页4 井已完成工程压裂。彭页1井测试产量2.5×104m3/d，地层压力系数1.0；彭页3 井测试产量3.8×104m3/d，地层压力系数1.1；彭页4井测试产量1.7×104m3/d，地层压力系数0.9。整体表现为靠近五峰组—龙马溪组露头区，地层压力系数和测试日产量逐渐降低。这与构造应力场模拟的实际情况相符，中燕山期，彭页1 井、彭页3 井、彭页4 井所受的最大主应力相差不大，均处于向斜核部，而在晚燕山期—喜马拉雅期，彭页4 井所受的主应力明显大于彭页1 井与彭页3 井，因此，日产气量与压力系数也明显降低[12]。

5 结论

1）通过野外节理和断层擦痕的地质行迹特征与交差叠加关系，对渝东南彭水地区及周缘的古应力场构造演化进行解析，包括古应力方向指示、构造演化时期以及构造演化时序。研究表明，渝东南彭水地区主要受到两期不同方向的构造应力作用，中燕山期，江南隆起由SE—NW 的递进扩展作用，致使渝东南地区经历了地质历史时期最强烈的一次构造作用，发育了大量NE 向、NNE 向地质行迹；晚燕山期—喜马拉雅期，印—亚板块的碰撞引起青藏高原的隆升，使得研究区处于SW—NE 的挤压应力场中，致使研究区多发育略近于菱形的叠加褶皱。

2）采用有限元的基本思想，借助Abaqus 软件精准建立地质模型和力学模型，确定了边界条件，对渝东南彭水地区五峰组—龙马溪组页岩储层的两期构造应力场进行数值模拟。将地质模型划分为页岩区和断裂带，通过三轴压缩实验得到页岩区与断裂带的力学参数，把力学参数赋值到地质模型中即成为力学模型，最后对力学模型的SE 和SW 方向分别施加79.62 和51.85 MPa 的梯度作用力，计算输出模拟结果。模拟结果表明，无论是在中燕山期还是在晚燕山期—喜马拉雅期，最大主应力的高值区主要分布在向斜核部，主应力值从核部到两翼逐渐降低，断裂带末端和断层的交汇处也处于高应力状态。根据数值模拟结果，结合岩石强度极限或破裂准则来预测裂缝的发育程度，以此进行压裂工程设计，可达到最佳的压裂效果。