基于最小耗能原理的库车坳陷超深致密砂岩裂缝定量预测

汪如军，赵力彬,2，张永灵

(1.中国石油塔里木油田公司，新疆库尔勒 841000；2.中国石油大学(北京)地球科学学院，北京 102249)

近年来，全球油气勘探逐渐从常规走向非常规、从浅层走向深层乃至超深层，深层致密储层受到普遍关注并不断取得突破，其中深层致密砂岩储层裂缝研究已经成为目前研究的热点之一[1-3]。构造裂缝是改善储层物性、提高油气产能的关键因素，构造裂缝的准确描述和有效预测，对裂缝性油气藏勘探开发具有重要意义。在多期构造运动作用下，储层裂缝发育不仅受控于区域应力场和构造等因素，还受到砂-泥岩组合、泥岩隔夹层的影响，进一步加剧裂缝空间分布的非均质性[4-6]，给裂缝的识别和预测带来极大困难。目前常用的储层裂缝的预测方法有地质分析估算法、构造曲率法、纵波各向异性及地震属性检测法、构造应力场模拟法和离散裂缝网络建模法(DFN技术)[7-12]。然而深层致密砂岩储层常面临着钻井密度低、取芯资料少、地震资料品质差的实际情况，仅采用单一手段很难对井间裂缝进行有效预测，因此从力学成因角度结合三维应力场模拟是目前最有效的方法[4,11-16]，其关键在于建立深层复杂岩体的破裂机制模型。因此选取库车坳陷超深致密砂岩A气藏为研究区，在经典力学理论的基础上，从新的理论角度引入反映自然界基本发展规律的新方法，完善并构建统一理论框架下的脆性岩石本构关系和破裂准则，推导裂缝定量表征模型，以解决裂缝表征中“质”(裂缝产生、扩展)的难题。依据构建研究区精细地质力学模型和古应力场有限元模拟，对深部储层裂缝的三维空间分布定量预测进行研究，并进一步完善裂缝预测的理论与方法体系，为这类油气田的勘探开发提供重要科学依据。

1 地质概况

库车坳陷位于天山南缘、塔里木盆地北部，是一个在古生代被动大陆边缘基础上发育演化成的中—新生代叠合型陆内前陆盆地[17-18]，K气田位于库车坳陷克拉苏构造带的东部(图1(a))。克拉苏构造带作为一个在中生代燕山运动和新生代喜马拉雅运动背景下形成的典型前展式叠瓦冲断构造[19]，目前已发现A2、A5、A6、A8、A9、A13、A24等大中型深层、超深层天然气藏(图1(b))。其中A气藏位于克拉苏构造带的中东部、K气田的中部，其南北分别受拜城断裂和克拉苏断裂两条逆冲断裂控制。整体上A气藏的构造形态简单，为一东西向延伸的长轴背斜，主要发育东、西2个构造高点。

A气藏巴什基奇克组(K1bs)储层埋藏一般为6 000～8 000 m，以粉砂岩、细砂岩为主[20]，夹有少量的中砂岩和泥岩薄层，平均厚度约为320 m(图1(c))。自上而下可划分为3个岩性段，简称为巴一段(K1bs1)、巴二段(K1bs2)和巴三段(K1bs3)。第1岩性段包含1～2砂层组，第2岩性段包含4～6砂层组，第3岩性段包含7～8砂层组。该区超深层砂岩储层具有低孔低渗透率特征，平均孔隙度仅为3%～5%，基质渗透率仅为(0.035～0.1) ×10-3μm2，而试井解释渗透率介于(1～10) ×10-3μm2，远高于基质渗透率，说明裂缝发育极大地改善储层渗流能力。截至目前，A气藏已累积生产天然气突破160亿m3，但受褶皱、断层、岩性、岩石力学性质等多种因素的影响和控制，井间裂缝表征难度大、预测精度低，致使气藏的整体高效开发难度增大。

2 裂缝发育特征

2.1 单井裂缝特征

岩心观察和成像测井(FMI)解释表明，A气藏巴什基奇克组(K1bs)发育大量构造缝，单井裂缝开度分布在0.1～0.4 mm。裂缝以NNW-SSE和NEE-SWW向为主，NE-SW和NW-SE向次之。主要以高角度缝(45°～75°)与直立缝(75°～90°)为主，分别约占总数量的60%和40%，低角度缝(15°～45°)和水平缝(0～15°)几乎不发育，仅约占0.1%。各井裂缝发育程度差异明显，A8-3、A8-4、A8-6井裂缝最为发育，平均线密度达到0.6～1.22条/m，其次是A8、A807、A8-1、A8003井，平均线密度约为0.5条/m。平面上，东、西部背斜高部位裂缝线密度略高于鞍部以及东、西两翼。垂向上，巴一段和巴二段裂缝密度较高，巴二段略高于巴一段，巴三段裂缝仅局部发育。

2.2 裂缝形成期次

岩石声发射测试、地层沉积埋藏史和构造演化史分析表明，克拉苏构造带自中生代以来主要经历4期构造运动，即晚白垩世—早古新世的燕山运动、晚渐新世—早中新世的早喜马拉雅运动、晚中新世—早上新世的中喜马拉雅运动和晚上新世—第四纪的晚喜马拉雅运动，其中上新世库车期—第四纪西域期来自于天山南缘的构造挤压作用表现最为强烈[19-20]，造就现今的叠瓦冲断构造格局，也正是A气藏高角度—直立有效构造裂缝形成的主要时期。

进一步，通过对大量岩心裂缝的产状、形态、组系、切割限制关系及充填程度等的统计分析，结合地质力学理论，认为A气藏巴什基奇克组内主要发育3期构造裂缝。第1期裂缝形成于白垩纪末—古近纪的近南北向伸展夹短暂NNW-SSE向挤压作用，最大有效构造应力位于35.2～59.9 MPa[19-20]，主要形成部分SWW-NEE向、近EW向的张性缝及少量近SN向和NW-SE向的共轭剪切缝，多数为无效裂缝。第2期裂缝形成于中新世康村组沉积期，受近南北向的天山挤压作用，最大有效构造应力约为74.8 MPa，主要形成一系列NW-SE向、近SN向的共轭剪切缝，伴随少量NE-SW向的剪切缝，仅少量裂缝保留一定的有效性。第3期裂缝形成于上新世库车组沉积期—第四纪西域期沉积期，此时最大有效构造应力达到80.9 MPa，形成大量NW-SE向、NNW-SSE向的共轭剪切缝以及部分近EW向平面共轭缝。同时导致克拉苏褶皱带发生强烈弯曲隆升，沿背斜长轴发育大量近EW向的高渗流能力的共轭纵张缝。第3期裂缝与深部天然气的充注匹配关系良好，促成工业规模性A气藏最终形成，也是本次裂缝定量预测研究的重点。

3 应力-耗能耦合的裂缝参数表征模型

3.1 最小耗能原理与岩石破裂准则

按照最小耗能原理，材料破坏总是在最薄弱地方以最易破坏的方式发生，即在耗能的任意时刻都满足最小耗能率原则[21-22]。因此从自然界的普遍性定律出发，结合经典力学理论，探寻非线性动态热力学过程中非均质岩石破裂与应力场、能量场的内在本质联系，改进和构建统一理论模式下致密砂岩的强度理论新体系，具有广阔的应用前景。

图2 库车坳陷A气藏构造裂缝形成序列及典型特征Fig.2 Formation sequence and typical characteristics of fractures in gas reservoir A in Kuqa Depression

3.1.1 岩石本构模型的建立

(1)

(2)

将式(2)代入式(1)得到耗能率表达式为

(3)

(4)

其中

G=E/[2×(1+μ)]，K=E/[3×(1-2μ)].

式中，α、β为材料参数，由试验数据拟合得到；φ12、φ23、φ13为内摩擦角，(°)；G为岩石的体积模量，MPa；K为岩石的剪切模量，MPa。

按照最小耗能理论，岩石在损伤演化过程中，耗能率都能在对应条件下取驻值，引入Lagrange乘子λ，由最小耗能原理得：

(5)

联立式(1)～(5)，得到三轴压缩情况下的岩石损伤方程为

(6)

且

(7)

式中，A为定义的公式，仅作公式缩写用。

当应变达到阈值后岩石损伤开始发展，在损伤阈值前处于弹性变形阶段。因而引入损伤阈值应变ε0后的岩石损伤方程为

(8)

其中

式中，λ、C为与岩性有关的常数。由此结合包含损伤变量D的广义胡克定律，得到岩石的本构模型为

(9)

3.1.2 岩石破裂准则建立

前人大量试验结果表明，在轴向应力达到峰值强度(用σc表示，MPa)的0.75、0.85和0.92倍时，脆性岩石中现存的大量微裂缝开始扩展联结并最终于峰值强度处形成宏观裂缝，之后则是沿着破裂面滑移形成断层[25]。因此针对具有强烈脆性特征的深层致密砂岩，采用应力-应变曲线的强度峰值点来辨识模型参数可使参数具有明确的物理含义。假定围压下岩石的应力-应变曲线峰值点处应力、应变分别为σ1s、σ2s、σ3s、ε1s、ε2s、ε3s，在峰值点处，

(10)

整理公式(6)并代入公式(10)得到常数λ、C0的值：

(11)

以及损伤临界值：

(12)

其中

故

(13)

从而得到最小耗能原理下的耗能率表达式为

(14)

根据最小耗能或最小耗能率原理，当φ>φ0时，岩石失去稳态，发生破裂。

3.2 裂缝参数表征模型建立

3.2.1 裂缝破裂角参数模型

图3 极坐标系下的裂缝破裂角计算示意图Fig.3 Diagram for calculating fracture angle in polar coordinate system

如图3所示，在主应力坐标系下，以坐标原点为极点，最大主应力轴σ1正半轴为极轴，建立合适的极坐标系，其中θ表示裂缝面与最大主应力之间的夹角，即破裂角。首先确定裂缝表面(r，θ)点的应力函数表达式F(r，θ)，设应力函数为

F=r2(acos2θ+bsin2θ+cθ+d).

(15)

故应力分量为

(16)

式中，σr和σθ分别为极坐标系下沿径向和环向两个方向的应力分量，MPa；τrθ为极坐标系下的剪应力分量，MPa。

将式(16)代入边界条件：

(17)

式中,σmax为最大主应力，MPa。

得到常数A、B、C、D的值，即

(18)

故应力函数表达式为

(19)

结合式(2)得到极坐标系下最小耗能原理约束的耗能率表达式为

(20)

通过最小耗能原理求解得到的裂缝破裂角表达式为

(21)

其中

3.2.2 裂缝走向、倾角参数模型

(22)

式中，αij(i,j=1,2,3)分别表示基向量eσ与x,y,z轴之间的夹角，(°)。

如图4(b)所示，裂缝面法向量n在空间直角坐标系中可表示为

(23)

式中，l、m、n表示法向量n的空间坐标；βij(i,j=1,2,3)分别表示法向量n与x,y,z轴之间的夹角，(°);θ为裂缝面与最大主应力之间的夹角，(°)。

在三维空间直角坐标系中，裂缝的走向和倾角一般用投影法来确定。空间直角坐标系的x轴与大地坐标的x轴(东)重合，z轴正方向与y轴负方向(南)重合，y轴和z轴重合。因此若将裂缝面法向量n投影到xOz平面，如图5(a)所示，则通过最小耗能原理预测的裂缝走向为

(24)

再从地质角度分析可知，裂缝倾角即为裂缝面与xOz平面的夹角，如图5(b)所示，故通过最小耗能原理预测的裂缝倾角为

(25)

图4 裂缝面法向量表征与主应力坐标系与空间直角坐标系转换示意图Fig.4 Characterization diagram of normal vector of fracture surface and transformation between principal stress coordinate system and space rectangular coordinate system

图5 裂缝走向预测与裂缝倾角计算示意图Fig.5 Schematic diagram of fracture strike prediction and fracture dip angle calculation

3.2.3 裂缝密度参数模型

根据损伤力学理论，连续性损伤变量D(0≤D≤1)可用来直接度量岩体的损伤程度，因此引入损伤变量D，结合最小耗能理论，建立最小耗能原理下的裂缝密度参数模型。对于脆性深层致密砂岩材料，当应力状态达到其破裂强度时，岩石即发生破裂，并释放出能量，结合式(14)计算可得产生裂缝过程中的耗散能(假设t=0时刻耗能为0)为

(26)

式中，ω为新增裂缝释放的的应变能，J ；t为应力作用的时间，s。

岩石破裂释放的能量，一部分用来新增裂缝总表面积，其余则以弹性波形式释放出去。对于裂缝来说，放出的弹性波能量很小，可忽略不计，释放的能量主要用来生成裂缝，即

ω=ωsV=SJ.

(27)

式中，ωs为新增裂缝表面积的应变能密度，J/m3；V为表征单元体的体积，m3；S为新增裂缝表面积，m2；J为产生单位面积裂缝所需要的能量，即裂缝表面能，J/m2。

将式(27)变化得到最小耗能理论下的裂缝体密度表达式为

(28)

其中Dv为单元体内裂缝密度,m2/m3。

进一步结合裂缝密度转化公式以及最小耗能理论，得到最小耗能原理下的裂缝线密度参数模型为

(29)

式中，Dl为裂缝线密度，条/m；L1、L3为在σ1、σ3方向上的长度，m；Dv为裂缝体积密度，m2/m3；θ为裂缝面与最大主应力之间的夹角，(°)。

4 三维裂缝定量预测

4.1 地质力学模型建立

基于巴什基奇克组顶面构造图，构建研究区的三维力学地质模型(图6)。该地区的大量有效构造裂缝主要形成于上新世库车期—第四纪西域期的NW向强烈挤压，此时A背斜构造形态已基本定型，与现今构造非常接近，可近似以现今地质模型为载体进行古构造应力场模拟。根据现场收集到的地震反演力学参数体和岩石力学试验结果，通过动-静态校正获得A气藏的静态非均质力学参数体，包括岩石密度、弹性模量、泊松比和内聚力、内摩擦角等(图6)。将巴一段、巴二段和巴三段分别设置为32.45°、36.21 °和30.26 °，断层带设置为27.19 °。然后将岩石力学参数按照“网格形心法”赋值到有限元地质模型中，构建地质力学模型，形成34 577个节点129 576个单元，垂向局部网格间距达到10 m。

图6 A气藏地质力学模型构建方法Fig.6 Methodology of building geomechanical model of gas reservoir A

以声发射测试应力值为基础，反复尝试、拟合确定A气藏地质力学模型的边界条件，在模型北部施加100 MPa的正向挤压应力，南部边界固定，以表征盆地的阻挡作用。第3造缝期天山南缘NW向的强烈斜向挤压在A地区产生相应的右旋剪应力，反复尝试确定剪应力值为20 MPa。在模型东西两侧设为自由边界，并施加大小为80 MPa的正向挤压应力，设置重力加速度为9.8 m/s2，同时在模型底部施加垂向位移约束。

4.2 古应力场模拟及裂缝预测

在古构造应力场数值模拟的基础上，将建立的耗能率、裂缝密度及倾角表达模型编成的计算程序导入有限元平台中，完成构造裂缝的三维定量预测。模拟结果显示，最大主应力的分布与构造形态关系密切，在背斜高部位为应力低值区，向翼部逐渐增大转变为挤压应力，拉张应力在背斜顶部零星分布(挤压应力为负值，拉张应力为正值)；由于A背斜受控于前缘主断层，天山南缘的挤压作用下，背斜南翼变形更为强烈，地层更为高陡，因此应力也更为集中，沿主断层出现明显的拉张应力区，在背斜顶部的次级断层下盘出现串珠状分布的应力低值区。总体上最大主应力方向较为稳定，以NNE-SSW向近SN向为主，反应的仍是区域挤压应力方向，仅在较大规模的断裂处有偏转(图7(a))。

图7 A气藏巴什基奇克组三维古应力场模拟及裂缝预测结果Fig.7 3D paleostress field simulation and fracture prediction results of Bashijiqike Formation in gas reservoir A

采用绝对线密度的概念，以1 m为间隔进行裂缝发育程度的统计，将裂缝预测结果与单井统计结果进行有效对比。平面上(图7(b))，预测结果在背斜的高部位呈现高值区，鞍部稍低，在两翼呈明显的低值区。背斜高部位的绝对裂缝线密度局部达到6条/m，一般为3.15～3.55条/m，如A803井、A8003井、A8-6井区，两翼为1.85～3.15条/m，如A8-5井、A8004井、A806井和A8-10井区。围绕断层常发育条带状或串珠状裂缝高值区，且断层规模越大特征越明显。A气藏隆升幅度较高，背斜形态受控于南部逆冲主断层和后缘(北部)反冲主断层，裂缝的发育与后缘反冲断层的强烈活动关系密切，东部背斜高点处的连片裂缝密度高值区也因此向北部迁移(图7(b))。垂向上，巴一段和巴二段裂缝密度明显高于巴三段，前两者的绝对裂缝线密度大部分为1.85～3.55条/m，局部高达6条/m，而后者的绝对线密度为0.15～3.15条/m。由于受钻井完钻深度的影响，单井统计得到的巴三段裂缝密度明显较低，仅在东部背斜顶部的东段，如A8003井、A8-3井和A8-6井区裂缝相对发育，绝对裂缝线密度达到3.55条/m。倾角作为直接反映裂缝有效性和储层连通性的重要参数，由图7(c)可以看出巴一段裂缝的倾角总体上要高于巴二段和巴三段，高角度缝在巴一段内大面积发育，直立缝主要沿着背斜高部位的长轴方向呈串珠状分布，这与巴二段的分布趋势非常相似。

4.3 裂缝预测结果的可靠性验证

A气藏的第3期裂缝以共轭网状缝为主，充填程度低，是关键的油气运移通道和储集空间。尽管第3期裂缝的线密度明显大，但不能完全代表裂缝发育的总密度，因此裂缝预测的可靠性应以井点实测数据矫正，同时考虑实际生产情况。从表1中可以看出，绝大多数井的分层平均裂缝线密度与成像测井解释结果相吻合，总体平均相对误差为23.94%，其中巴一段的平均相对误差为17.45%，巴二段和巴三段的平均相对误差分别为25.47%和23.41%，这很可能是K气田巴什基奇克组储层内存在大量泥岩夹层导致的。另外A8004井的巴二段裂缝预测相对误差达到66.67%，A8井巴一段和克A8003井巴二段的预测相对误差也分别达到46.67%和39.62%，与单井实测结果相差大。但从沿背斜长轴方向的裂缝绝对密度剖面上看，以A8井巴一段为例，其成像测井得到的裂缝平均线密度为0.3条/m，预测结果为0.44条/m，两者存在一定的误差，但从岩心观察统计结果来看，其裂缝平均线密度为0.52条/m，与成像测井解释结果的误差仅为18.18%。原因在于当井壁被泥浆污染或该段岩性变化大时，会造成成像测井显示效果差，对于一些细微的裂缝很难识别。

从沿背斜长轴方向的裂缝密度预测剖面上看，裂缝绝对密度的预测结果与单井成像实测结果吻合度高。尤其是在背斜高部位的A8-3井、A8井、A8-2井、A8003井和A8-6井等，其巴一段和巴二段裂缝发育高值区呈现连片分布趋势，显示出较好的裂缝垂向和平面上连通性。垂向上，在背斜高部位及背斜转折端处，裂缝发育程度有着从浅部向深部降低的趋势，即巴一段、巴二段、巴三段裂缝绝对线密度依次降低，但在背斜翼部具有相反的趋势，即巴二段底部和巴三段顶部的裂缝更为发育。

表1 A气藏裂缝平均线密度预测结果与单井实测结果对比Table 1 Comparison between prediction results of fracture density and measured results of single wells in gas reservoir A

为了进一步有效验证裂缝预测参数的准确性，以气藏整体生产情况为基础，考虑裂缝走向和充填程度等因素，对裂缝密度和倾角预测结果进行分析。裂缝的走向在东西向上变化较大，A气藏西部次级背斜的裂缝走向以NNW-SSE、NW-SE和NE-SW向为主，发育少量近EW走向裂缝。而气藏东部次级背斜的裂缝走向由西往东逐渐变为NWW-SEE向和NEE-SWW向为主，NE-SW向和近SN向次之，其中NEE-SWW向和NWW-SEE向裂缝形成一组共轭裂缝系，成因上属于高角度晚期平面“X”型节理，走向与背斜枢纽近于平行。近EW向裂缝与背斜枢纽近于平行，成因上属于晚期纵张节理。两种成因的裂缝主要形成于喜山运动晚期，倾角大，平面连通性强，充填程度低，有效性好。A气藏无阻流量高于400 m3/d的单井主要是位于西部的A807井、A802井和中部的A8-1井、A8井、A8-3井、A8-2井、A8003井、KeS8-6井和A8-8井，而低产井主要位于背斜的东、西两端。研究区高产井基本上都位于直立裂缝发育区，相应的裂缝发育程度为中等和高。尤其是A8井、A8-3井均位于裂缝密度高值区，裂缝倾角近似直立，矿物充填率分别为0.26和0.21，其无阻流量达到523和508 m3/d。相比之下，位于背斜倾没部位的A8-11井、A801井和A8-10井产能均小于300 m3/d，尽管矿物充填程度较低，但裂缝发育程度为中等—较低，主要为高角度和低角度斜交缝，整体连通性和有效性一般，这些都进一步验证裂缝预测结果的可靠性。综合以上分析表明，裂缝密度、倾角和充填程度是关系有效储层发育和分布的关键因素，A802井西北地区、A8井北部和A8003井西部地区是下一步部署高效开发井的有利区带，当然发育密度高的直立裂缝也是造成气藏见水的重要因素(图8)。

图8 A气藏裂缝预测结果与实际生产动态数据对比Fig.8 Comparison of fracture prediction results with actual production dynamic data in gas reservoir A

5 结 论

(1)将最小耗能原理与经典力学理论相融合，探寻非线性动态热力学过程中岩石破裂与应力场、能量场的内在本质联系，改进和完善统一理论模式下致密砂岩的强度破坏准则和裂缝参数表征模型，适用于陆相超深层致密砂岩储层的裂缝研究。

(2)裂缝在A背斜高部位最为发育，在鞍部稍低，在两翼呈低值；东部背斜高点处裂缝密度高值区连片发育，明显高于西部背斜高点；褶皱是A气藏裂缝发育的首要主控因素，断层次之，岩性仅在局部影响着裂缝发育；垂向上，在背斜高部位以及背斜转折端位置，裂缝发育程度有着从浅部向深部降低的趋势，但在背斜翼部具有相反的趋势，即在巴二段底部和巴三段顶部裂缝更为发育。

(3)裂缝密度、倾角和充填程度是控制有效储层发育和分布的关键因素，高产井基本上都位于直立裂缝发育区，相应的裂缝发育程度为中等和高，A802井西北地区、A8井北部和A8003井西部地区是下一步部署高效开发井的有利区带，高密度直立裂缝也是造成气藏见水的重要因素。