地应力预测是长期以来困扰页岩气勘探工作的难题之一。由于页岩储层是低孔低渗岩石, 因此需要经过大规模的人工压裂改造才能够获得工业气流, 而地应力是决定压裂改造效果的关键因素之一。在页岩气层的改造过程中, 地层的地应力状态直接决定了人工压裂造缝的效果, 选择地应力适当的区域进行水平井轨迹设计和压裂施工方案设计, 能够有效提高页岩气水平井单井产能。

在地应力预测方面, 研究者通常采用成像测井统计分析或数值模拟等方法。成像测井统计分析法是收集钻井井眼变化^[1-3]、区域应力场^[4]和成像测井显示发育情况^[5-6]等信息开展统计分析, 对区域的地应力场进行间接分析^[7], 但是将这类方法应用于勘探程度较低、构造应力场复杂的区域时, 存在稳定性差、预测精度低的问题; 数值模拟法则是采用有限元法数值模拟地层在背景应力边界条件下的各地质单元的应变和应力分量, 求取单个单元的三向应力分布和应力方向。BEAUMONT等^[8]和KCKERT^[9]将数值模拟法应用于构建俯冲—碰撞构造的力学模型以及模拟地壳的剪切应力状态; PARSONS^[10]使用GPS观测加利福尼亚地区地形变化, 以有限元法模拟分析了该地区构造应力分布; 刘雄^[11]则在东营原油开采区开展了有限元法地应力预测。这些研究结果表明, 有限元法预测准确度相对较高, 但十分依赖数据样本的数量和基础地质资料, 在勘探程度较低和构造变化复杂的地区, 其预测准确度严重下降。同时, 有限元法在计算时由于区域环境的复杂以及断层位置的形变引起计算网格形变, 致使界面追踪处理困难, 降低了计算效率甚至导致计算失败。此外, 还有研究者提出通过叠前反演求取地层参数, 进而利用水平应力差公式^[12]或支持向量机^[13]等方法求取地应力分布的方法, 但这类方法需要通过经验系数对计算公式进行校正或通过机器学习进行拟合, 对工区研究程度要求较高, 并不适用于勘探程度较低的地区。

物质点法(material point method, MPM)是SULSKY等^[14]提出的一种数值模拟方法, 该方法将携带所有物质信息的质点和背景网格连接, 通过插值函数完成物质点与背景网格之间的信息交互表征区域的运动和变形状态。物质点法计算时采用规则的背景网格计算动量方程, 从而实现了质点间的相互作用与联系, 避免了网格畸变导致计算结果误差极大的问题, 并提高了计算效率^[15]。物质点法最初主要用于模拟材料的大型变形和断裂破碎的问题, 而在页岩气研究中, 物质点法能够结合外部区域构造背景应力对内部断裂发育导致的应力释放所引起的地应力异常区域进行预测。国内外研究者大量的地应力研究结果表明, 断层不论其构造规模大小, 都能对其附近的应力状态产生一定的影响^[16], 因此可以认为, 断裂的发育状态能够对地应力进行表征, 同时这些结果也可以被用来识别通常地震资料无法观测的小型断裂, 也可以被用来模拟压裂后的裂缝发育情况。

本文针对丁山地区构造环境复杂、勘探程度较低的现实问题, 利用物质点法开展本地区地应力差的预测, 对提高丁山地区页岩气井单井产量具有重要的现实意义。

1 方法技术 1.1 计算步骤

地应力差为地层实时最大主应力与最小主应力之差, 是控制裂缝延伸变化和地层破裂压力的关键参数。地应力差越大, 裂缝延伸的方向性越强, 地应力差越小, 裂缝延伸越趋于复杂化, 因此可以认为相对较小的地应力差更为有利于水平井压裂造缝。

利用物质点法进行地应力差模拟的技术路线如图 1所示。其具体计算过程为: ①定义背景网格, 将所有质点岩石物理参数置入背景网格上(图 2a); ②采用背景应力场条件计算质点应变并更新其应力(图 2b); ③利用物质点积分获得背景网格节点的节点力, 在背景网格上积分动量方程(图 2c); ④分别采用背景网格节点的应力场和应变场更新质点的位置和应变量(图 2d); ⑤舍去变形的背景网格并构造新网格, 继续迭代计算直到计算结果收敛, 最终实现应力差预测(图 2e)。

在具体计算时, 该方法(公式(1))的核心是在任意物质点处计算背景网格节点的动量方程, 用以求取背景网格的应变量。

$ \dot{P}_{i I}=f_{i I}^{\mathrm{ext}}+f_{i I}^{\mathrm{int}} $

(1)

式中: I代表网格节点; $\dot{P}_{i I} $代表动量; f_iI^int和f_iI^ext分别代表背景网格节点的内力(表示物体断裂后的应力)和外力(表示背景在网格施加的应力)。

第I个网格节点在i方向的动量为:

$ P_{i I}=\sum\limits_{J=1}^{n_g} m_{I J} \dot{u}_{i J} $

(2)

式中: n_g为计算网格的节点总数; m表示质量; J表示区域; u_iJ表征节点在区域J的i方向的位移。

背景网格的质量矩阵为:

$ m_{I J}=\sum\limits_{p=1}^{n_p} m_p N_{I p} N_{J p} $

(3)

式中: p表示物质点; n_p为物质点的总数; m_p为物质点所代表区域的质量; N_IP为节点I的形函数在物质点p上的值; N_JP为区域J的形函数在物质点p上的值。

$ f_{i I}^{\mathrm{int}}=-\sum\limits_{p=1}^{n p} N_{I p, j} \sigma_{i j p} \frac{m_p}{\rho_p} $

(4)

$ f_{i I}^{\mathrm{ext}}=\sum\limits_{p=1}^{n_p} m_p N_{I p} b_i+\sum\limits_{p=1}^{n_p} N_{I p} h^{-1} \bar{t}_i \frac{m_p}{p_p} $

(5)

式中: σ_ijp为物质点p的应力; ρ_p为物质点p的密度; b_i为体力; h为层厚度; t_i为i方向的面力。

1.2 关键模型构建 1.2.1 裂缝模型

利用物质点法求取应变量时, 第一步首先需要构建一个等效天然裂缝模型(equivalent fracture model, EFM), 无天然裂缝时, 将该模型整体视为一个等效均匀介质(图 3a), 而当网格中存在裂缝属性时, 则将裂缝位置的质点设置为空, 用空的位置代表存在天然裂缝(图 3b), 无裂缝发育的网格则离散为64个物质点。

1.2.2 计算模型

第二步则是利用由三维地震资料求取的目的层的弹性参数为每个物质点赋予平均杨氏模量、平均泊松比、平均密度和平均断裂韧性等岩石物理参数信息^[16], 并将最大水平主应力方向、平均最大水平主应力值和水平应力差异系数作为背景应力场信息赋予整个计算模型(图 4)。

2 应用实例及分析 2.1 实验区研究背景

丁山构造位于四川盆地东南部(图 5a), 齐曜山断裂带以东, 主体为推覆断背斜构造, 区内主要勘探目的层为志留系龙马溪—五峰组, 普遍发育黑色硅质页岩(图 5b), 为硅质深水陆棚相沉积^[17], 是良好的页岩气储层。丁山地区处于推覆背斜的北翼斜坡, 推覆断层纵向未延伸至志留系地层内。受川南燕山期以及燕山—喜山两期构造影响, 区域内构造应力十分复杂, 断裂和裂缝广泛发育。前人研究认为, 丁山地区整体最大主应力呈双簇特征^[18], 分别为NEE向和NW向(图 6)。

目前丁山地区龙马溪组开发水平井压裂后产量较低, 前人分析认为造成这种现象的主要原因是由于区域内地应力情况不明导致造缝效果较差^[19-20]。因此寻找适用于丁山地区的地应力预测方法十分重要。

2.2 地应力差场的计算

第一步, 将丁山三维五峰组深度域层位(图 7)抽稀为100 m×100 m的构造网格, 使用丁山三维地震资料提取I3D裂缝属性(图 8, Inflated 3D ConvNet, 一种基于图形深度学习的裂缝预测属性)用以构建研究区地应力背景网格, 将五峰组底界的I3D属性赋值到构造面网格上, 即可得到含有连续裂缝模型(CFM)信息的构造面网格(图 9)。

第二步, 将杨氏模量和泊松比等参数赋值于CFM构造面网格, 最终建立基于连续裂缝模型提取的包含岩石物理参数的等效天然裂缝模型(图 10), 具体参数见表 1。

第三步, 在使用基于包含背景应力场信息的模型估算赋予岩石物理信息的物质点时, 含有天然裂缝的地层会在背景应力作用情况下发生微弱形变。模拟含有裂缝的地层在远景应力场作用下发生的应变, 根据最大水平应力方向和最小水平应力方向的应变差异(具体参数见表 1), 即可估算得到目的层相对应力差平面展布。提取这种应力梯度微弱变化与初始应力梯度分布之间的差异数据, 即可形成相对地应力差分布图。图 11为丁山构造相对应力差场平面分布。

2.3 预测有效性分析 2.3.1 成像测井

在利用成像测井对地应力进行解释时通常采用椭圆拟合法或图像观察法(诱导缝、井壁崩落)。研究区成像测井显示丁山地区最大水平主应力方向为近东西向, 局部有差异(图 12)。两者结果对比显示, 物质点法所获取的地应力差分布, 与实际钻井显示的最大水平主应力方向结果基本一致, 具体预测结果及有关参数见表 2。

与成像测井单井分析相比, 基于物质点法的地应力差预测经济成本更低, 预测结果能够有效指导勘探早期井轨迹设计和压裂施工方案设计。

2.3.2 Zoback应力估算法

在Zoback法估算中, 摩擦系数通常使用工区经验参数, 因此可以认为, 埋深决定了地应力场的数值大小, 而弹性和地层压力等参数主要用于求取地层孔隙压力, 以调整地应力预测细节。利用叠前反演泊松比、杨氏模量、地层压力预测等成果, 通过在井点确定的应力计算模型进行三维地应力估算(图 13), 认为丁山地区构造高部位应力差异系数均小于0.13, 储层压裂改造时易于形成网状缝, 具体预测结果及有关参数见表 2。

与Zoback应力估算法相比, 基于物质点法的地应力差预测结果细节更丰富, 预测结果能够有效描述地应力差的方向变化细节, 为水平井轨迹设计和压裂施工方案设计提供更准确的依据。

2.4 预测结果讨论

从地应力差模拟结果看, 由于研究区内受多期次构造运动影响, 构造主体表现出主体地应力差相对较低, 周缘地应力差较大的特点。但是由于区内受多期次天然裂缝发育影响, 构造主体内地应力沿NE向和NW向两个方向断层发育的变化趋势, 造成这种情况的原因主要受到断裂发育的影响, 区内地应力释放程度差异不同, 近断层的区域地应力变化相对远断层区域较大。

统计单井平台无阻流量和测试产量(表 2), 可以看出, A1H井、A3H井和A2H井相对靠近断裂部位, 地应力差相对较高, 导致压裂改造中人工缝发育不规则造成单井产量相对较低。而A4H井、A5H井附近基本不存在大型断裂, 地应力差相对适中, 压裂后人工缝情况相对较好, 测试产量明显较高。

分析结果表明, 在大型挤压、拉张等构造活动背景下能够形成一定的地应力发育倾向, 而断裂的发育则决定区内地应力释放的具体情况, 断层部位地应力差相对较大, 不利于人工压裂缝均匀发育。因此可以认为在平行断裂走向且距离断裂相对较远的位置更有利于水平井轨迹设计及后期开展压裂改造。

3 结论

1) 基于物质点法的地应力场差异预测技术是一种考虑宏观构造背景和微观断裂缝发育的数值模拟方法。与常规使用的地应力分析和预测手段相比, 本方法对构造应力场复杂且断裂发育地区页岩气层地应力场差异性评价具有良好的效果, 预测结果对下一步开展水平井轨迹设计具有现实指导意义;

2) 丁山地区受多期次构造运动影响, 构造主体地应力差相对较低, 主体内部NE向和NW向两个方向地应力差异较大。同时由于区内地应力释放程度存在差异, 近断层轴的区域的地应力变化相对远断层区域较大, 平行断裂走向且距离断裂相对较远的位置更有利于水平井轨迹设计及后期开展压裂改造。