四川盆地页岩气资源分布广、资源总量高，埋深3 500～4 500 m的深层页岩气资源占85%以上^[1]。为了实现页岩气自主商业化开发，通过调研国外压裂技术，突破了常规单缝增产改造经典模型，形成了以“提高储层改造体积（SRV）、提升压裂缝网复杂度”为目标的第一代大规模体积压裂改造技术，成功实现了中深层（埋深小于3 500 m）五峰组—龙马溪组页岩气的规模效益开发，其中储层改造体积（SRV）是评价压裂效果的关键指标^[2–3]。随着页岩气开发转向深层和超深层，大规模体积压裂技术不适用于高闭合压力、高应力差和复杂纹层储层。为了突破深层和超深层油气资源的压裂开发，国内自主攻关形成了以“压裂裂缝表面积最大化”为目标的第二代页岩密切割高强度压裂技术。该技术在第一代压裂技术基础上强化了施工参数，通过采用段簇加密、暂堵转向和提高加砂强度等措施大幅提升了深层页岩压裂效果^[4–5]。由于第二代压裂技术和第一代体积压裂技术在理念上明显不同，传统的SRV评价方法不再适用，需要建立储层改造面积（SRA）估算方法，为深层超深层密切割压裂效果评价和施工参数优化提供支撑和依据。

目前的三维缝网构建方法主要分为裂缝扩展方法和缝网重构方法。裂缝扩展方法的本质是将缝网内的基本关系简化成裂缝与裂缝的连通关系，在得到天然裂缝分布之后，根据水力裂缝−天然裂缝扩展准则的假设，基于一定的数值计算方法得到水力裂缝在天然缝网中的扩展路径，从而得到最终的三维复杂缝网。为了更细致地模拟水力裂缝−天然裂缝的交叉行为，要用到数值方法，而目前常用的数值方法有有限元方法^[6]、扩展有限元法^[7]、位移不连续法^[8]、离散元方法^[9]和近年来发展起来的相场模拟方法^[10]等。缝网重构是将缝网内的基本关系简化为事件点和裂缝面的拓扑关系，在得到多级压裂过程中的微地震监测事件点后，通过假设合理的缝网−事件点分布关系，采用具有高鲁棒性的回归方法从事件点中直接识别三维复杂缝网。常用的回归方法有时序拟合法^[11]、迭代连接法^[12]、条件随机拟合法^[13]、分形拟合法^[14]、震源机制拟合法^[15]等。

传统的缝网面积计算方法主要基于裂缝扩展方法和缝网重构方法重建页岩压后的三维缝网，即先利用模拟结果和监测数据建立缝网模型，然后通过缝网去估算裂缝和与储层基质的接触面积。对于裂缝扩展方法，水力裂缝−天然裂缝的扩展准则一般是基于室内的物理模拟试验建立的。裂缝扩展方法得到缝网的过程和实际裂缝的扩展机理比较一致，其缺点是很难精确描述天然裂缝的三维分布，目前基本都是采用条件随机模拟的方式来弥补该方法的这个缺点^[16–17]。由于页岩中含有胶结程度不同的层理弱面，水力裂缝和天然裂缝的作用极为复杂，而水力裂缝和天然裂缝的作用又直接决定了缝网的复杂度，这些因素使扩展模型模拟得到压裂缝网的不确定性大大增加。对于缝网重构方法，缝网重构将缝网内的基本关系假设为事件点和裂缝面的拓扑关系，得到的缝网和原始监测事件点之间具有很高的贴合度，但是只考虑了裂缝和事件点几何上的相关性，因此扩展机理难以与实际情况相符，准确度有限。

裂缝扩展方法和缝网重构方法都存在由于假设带来的缝网误差问题，这些假设降低了缝网面积计算的准确度。为了尽可能降低缝网重构和缝网模拟过程中假设带来的计算误差，笔者将利用微地震事件直接提取SRV的思路再次应用到SRA的估算上，尝试在识别微地震点云空间几何结构的基础上，构建一种高精度的基于微地震点云的SRA估算方法，在不重构缝网的情况下尽可能提高SRA的准确度。

1.   SRA计算模型的建立

1.1   微地震八叉树分解

八叉树分解是一种计算机图形学的点集体素化分解方法^[18–21]，最重要的特点就是能够充分考虑点集在空间密度上的非均匀分布，并根据密度的分布生成非均匀的包围网格。八叉树网格是一种多层次树状网格结构，每一层的节点根据密度的阈值由上一层高密度的节点迭代生成。

对于特定的事件点集，首先根据第N层节点集分布的最大边界识别出点集对应的最小外接长方体，也就是AABB包围盒。其次，在第N层的划分中将根节点分成8个大小相同的子节点，获取第N层划分的节点信息，并存储每个节点内部的点集信息。然后，根据活动节点阈值$t$判定第N层网格中的活动节点和非活动节点，如果第N层中的任意节点满足活动节点判定条件，即节点内部包含的事件点数量大于或等于活动节点阈值$t$，则认定该节点为第N+1层活动节点，其余不满足条件的认定为第N+1层非活动节点。在得到第N+1层的非活动节点之后，继续对节点进行非空的判定，当节点满足非空节点判定条件时，得到第N+1层非活动非空节点；当节点不满足非空阈值时，得到第N+1层非活动空节点。这样，第N层节点通过八叉树分解划分为第N+1层活动节点、第N+1层非活动非空节点和第N+1层非活动空节点（见图1）。第N+1层非活动节点保持不变，同样将每个节点均匀划分成8个相同的子节点，在对第N+1层中的活动节点进行划分后得到了包围网格的第N+2层的节点，迭代划分，直至所有活动网格都被划分完毕，停止网格的划分，得到如图2中右图的多层次网格结构和图3的节点层次结构。图2中，左侧的立方体网格和右侧的树状节点分布完全对应，高密度事件点集在划分中被进一步细分成高层次的网格，低密度事件点集主要存在于较低层次的网格中，网格中的空白区域在划分过程中被过滤清除，最终得到的网格和点集具有良好的紧密度，能够用近似点集的结构提取计算点集的空间参数。

图  1  八叉树节点划分过程示意

Figure  1.  Octree node decomposition process

下载: 全尺寸图片 幻灯片

图  2  微地震事件点集的八叉树网格划分示意

Figure  2.  Octree grid decomposition of microseismic event point set

下载: 全尺寸图片 幻灯片

图  3  微地震事件点集对应的多层次网格结构

Figure  3.  Multi-level grid structure corresponding to microseismic event points

下载: 全尺寸图片 幻灯片

图3更直观地给出了每一层节点的分类和划分过程，设活动网格阈值为$t$，以节点内部包围的点集数量为依据，每一层的节点根据活动网格阈值t可以划分为活动节点和非活动节点，而非活动节点又可以进一步划分为非活动非空节点和非活动空节点。树状结构采取分层编码的方式，令$i$为节点层次，$j$为节点在单层内部的编码，$g(i,j)$为第$i$层第$j$个节点内部包围的事件点数量，则4种节点的布尔判断标准为：

式中：$f(n(i,j))$为每一层节点的非活动节点和活动节点判断指标，为1时代表活动节点，为0时代表非活动节点；$g(n(i,j))$为非活动节点中非空节点和空节点的判断指标，为1时代表非活动非空节点，为0时代表非活动空节点。

令事件点集为$P$，对应第$i$层划分得到的节点集合为${P_i}$，随着节点划分层次的增加，可以得到：

式中：${A_i}$为第$i$层节点集合${P_i}$中的活动节点集合；${I_i}$为第$i$层节点集合${P_i}$中非活动节点集合；${D_i}$为非活动节点集合${I_i}$中的非空节点集合；${E_i}$为非活动节点集合${I_i}$中的空节点集合。

由于设置空节点的判断标准为内部至少包含一个事件点，显然当${E_i}{\text{ = }} \varnothing$时，${I_i} = {D_i}$。如果空节点的判断标准为至少包含${t_1}$个事件点，且满足$1 \leqslant {t_1} < t$，则${E_i}$不为空。

虽然对于树状网格结构来说，下一层活动节点加上非活动节点等于上一层的活动节点。第0层的节点是根节点，采用了AABB包围盒包围所有的事件点，因为全部参与了第二层的划分，第0层的节点都是活动节点；最后一层的节点是节点划分的终点，全部节点都不满足活动节点的判断标准，因此最后一层全部都是非活动节点，即所有的节点最终都会被划分成非活动节点，而点集最终都会存储在非活动节点当中。从图2可以看出，随着节点划分层次的增多，所有的活动节点最终都会被划分成非活动节点。因此，可以把每一层的非活动节点的累加当作对点集的一种分解，可以得到：

去除非活动节点内部的空节点，进一步可以得到：

至此，完成了点集空间参数提取的最关键一步，建立了和事件点集等效的体素化网格，即用树状网格结构去近似表征点集不规则的复杂空间结构。

1.2   基于八叉树分解的SRA估算模型

对于单个节点，八叉树网格结构去除了空节点，网格结构和原始点集具有很高的紧密度，节点内部点集所对应的特征尺度和节点的尺度相同，即点集占有的特征区域完全充满了整个节点。因此，可以节点为基本单位进行节点内部点集空间参数的估算，如图4所示。从图4可以看出，对于单个节点，上文提出的模型采用该节点的特征面积去估算该处点集所对应的裂缝面积，微地震点集所对应的三维复杂缝网面积SRA则可以等效为点集八叉树网格中所有非活动非空节点的特征面积之和，这样就可通过八叉树体素化分解实现无需重构缝网情况下对SRA的有效估算。

图  4  单节点裂缝面积和节点底面积示意

Figure  4.  Fracture area and base area of single node

下载: 全尺寸图片 幻灯片

采用节点的3个侧面面积的几何平均值作为节点的特征面积，由于根节点对应的底面积为点集AABB包围盒的特征面积，其节点层次为第0层，可以得到根节点的特征面积为：

下一层节点都来源于上一层的八等分，因此对应第$i$层的节点的特征面积${S_i}$可以根据根节点得到：

因此，对于点集$P$，采用布尔判断方法，可以将对应的三维复杂缝网面积SRA等效为所有非活动非空节点的特征面积之和，得到$S_P$为：

式中：S₀为根节点的特征面积；$\alpha (i)$为第$i$层的节点总数量；$S_P$为点集$P$对应的等效缝网面积。

2.   SRA模型模拟验证

2.1   模型验证方法

验证SRA模型时，首先需要确定的是微地震事件所对应的真实缝网，只有将模型得到的SRA估算结果与真实缝网面积进行对比，才能定量评价SRA模型的计算误差。然而，受当前技术的限制，无论是测斜仪、示踪剂还是微地震监测技术，都不可能测得实际多级水力压裂改造中真实的三维复杂缝网，只能得到局部的缝网参数。由于无法获取真实的缝网，也就无法利用真实的微地震监测数据来定量评价模型的误差。为此，为了定量评价SRA模型的准确性和适应性，笔者开发了一套基于模拟微地震事件点的SRA模型模拟验证和定量分析流程，如图5所示。

图  5  八叉树缝网面积计算模型模拟验证流程

Figure  5.  Simulation and verification process of fracture network area calculation model based on octree grid

下载: 全尺寸图片 幻灯片

首先，采用蒙特卡罗方法生成三维复杂的离散裂缝网格，利用高斯混合模型对缝网进行离散，得到和缝网分布一致的模拟事件点集，同时可以得到三维离散裂缝网格的真实缝网面积。其次，把事件点集代入SRA计算模型。利用模拟事件点集生成对应的多层次树状八叉树网格，根据活动节点判断阈值和空节点判断阈值提取每一层节点中的非活动非空节点，这个过程去除了网格中不含事件点的空网格，得到和事件点集分布具有较高贴合度的三维等效网格。在获取原始缝网对应的等效网格之后，提取等效网格的特征面积就可以得到事件点集对应的等效面积。最后，对比离散前的事件点集对应的真实缝网面积和采用八叉树网格计算的等效缝网面积，就可以计算出模型的相对误差。

上文模型基于开源的ADFNE裂缝模型^[22]生成缝网，单裂缝模型采用上述的多边形裂缝，裂缝的方位参数采用von-Mises分布进行模拟。在计算事件点集等效缝网面积过程中，最关键的一步是提取事件点集对应的等效网格结构。图6所示为不同数量裂缝缝网事件点集对应等效网格的提取过程。

图  6  不同数量裂缝缝网事件点集对应等效网格的提取

Figure  6.  Extraction of equivalent grids corresponding to fracture network event point sets with different numbers of fractures

下载: 全尺寸图片 幻灯片

从图6可以看出：相比全部的八叉树网格，去除了空网格的八叉树网格和原始事件点集的分布高度一致；随着模拟缝网中裂缝数量的增多，网格节点的划分层也随之增加，网格划分也更加精细，得到的非空八叉树网格和事件点集的紧密度也会更好。由于八叉树网格划分过程中充分考虑了事件点密度的分布，对高密度事件点区域网格加密，使整体面积的估算更加准确。

2.2   模型参数设置

工区体积设置为边长为1 km的立方体，裂缝数量设置为1～30条。根据活动节点阈值和事件面密度的模拟结果，将活动节点阈值设置为1，以确保网格划分的层次最大且和事件点具有最高的贴合度；事件点面密度为2 500个/km²，这接近缝网层次稳定对应的最小事件点密度。通过改变对应缝网的裂缝数量，分别模拟计算不同缝网条件下的等效缝网面积S_e和真实缝网面积S_t。

2.3   模型模拟验证结果

三维缝网的复杂度和裂缝数量有着直接的关系，裂缝数量越多，缝网之间的交叉关系就越复杂，缝网的复杂度就越高，计算过程中裂缝之间的干扰就越大。相对于工区尺度，这里采用中等尺度的缝网进行模拟事件点的生成，裂缝数量变化范围为1～30，分别模拟计算不同裂缝数量下等效缝网面积的计算误差，结果如图7所示。

图  7  不同裂缝数量缝网等效缝网面积和真实缝网面积的模拟计算结果

Figure  7.  Simulated calculation results of equivalent fracture network area and actual fracture network area with different numbers of fractures

下载: 全尺寸图片 幻灯片

从图7可以看出：随着裂缝数量增多，等效缝网面积和真实缝网面积高度一致；30条裂缝缝网等效面积与真实面积的相对误差为12.31%，具有较高的计算精度。图8给出了裂缝数量小于15条和裂缝数量大于15条缝网等效面积计算误差的分布，其中1～15条裂缝缝网等效面积与真实面积的平均相对误差为20.01%，16～30条裂缝缝网的平均相对误差为4.62%，1~30条裂缝缝网的平均相对误差为12.32%。对比图7也可以明显发现，误差主要集中分布在缝网规模比较小的区间，当裂缝数量比较多时，计算误差反而保持在非常小的范围内，平均误差仅为4.62%。相比较少的裂缝数量，当离散事件点面密度为2 500个/km²时，15～30条裂缝缝网三维网格结构和事件点集的等效程度更好，计算精度也更高，为高密度事件点集的划分提供了依据。

图  8  不同裂缝数量缝网等效面积的平均相对误差

Figure  8.  Average relative errors of equivalent fracture network area with different numbers of fractures

下载: 全尺寸图片 幻灯片

从模拟计算结果看，在不同裂缝数量和交叉程度条件下，上文提出的SRA计算模型具有良好的计算精度，同时也验证了通过等效的八叉树网格结构来近似提取事件点对应真实缝网参数的可行性。

3.   应用实例

3.1   地质和工程背景

水平井平台B位于四川盆地，该平台有水平井B1、水平井B2和水平井B3等3口水平井，平均完钻井深4 650.00 m，水平段长度平均1 500.00 m，目的层为龙马溪组黑灰色页岩。平台采用工厂化拉链式压裂，3口井的段间距为60.00～100.00 m，以低黏滑溜水为压裂液，以70/140目粉砂+40/70目陶粒为支撑剂，按照12～14 m³/min的排量施工，单段液量控制在1 800 m³左右，以提高缝内静压力，促进水力裂缝和天然裂缝的沟通。为了获取B平台水力压裂过程中裂缝的几何形态、展布方位、裂缝复杂度和气藏增产体积等关键信息，布设三分量检波器阵列监测压裂过程中的微地震信号，将采集到的微地震信号进行去噪处理和定位解释，发现3口水平井63个压裂段共有16293个有效事件点。

图9所示为B平台3口水平井压裂过程中监测到的微地震事件点。从图9（a）和图9（b）可以看出，水力缝网沟通了储层中的天然缝网，形成了复杂度较高的裂缝网络。从空间上，裂缝在横向上主要沿最大主应力方向延伸，和井筒存在夹角，在纵向上延伸较窄。从B平台3口井的事件点平面分布看，B1井的事件点和B3井的事件点分布最广，受周围天然裂缝影响，压裂过程中压力的传递激发了一些远井地区的裂缝，事件点中包含了大量的远井区域的低密度事件点。B1井和B3井的高密度事件点主要集中在前半段和后半段，中间段的事件点相对稀疏，压裂效果较差，两端的压裂效果较好。相对B1井和B3井，B2井的高密度事件点主要集中于中间段，两端的压裂效果差一些，事件点的空间延伸范围要比B1井和B3井都要小，事件点覆盖范围最小。

图  9  B平台水平井压裂过程中监测得到的微地震事件点

Figure  9.  Microseismic event points monitored during the fracturing process of horizontal wells on platform B

下载: 全尺寸图片 幻灯片

为了定量对比3口井的事件点，分别统计了3口井的事件点数量，并采用核密度统计分析方法对3口井微地震事件点密度进行了空间估计，结果如图10所示。从图10可以看出，B2井事件点的数量和密度都显然大于B1和B3井。核密度分析图显示，高密度事件点主要围绕在B2井的周围分布，以B1和B3井为边界，B1井和B3井激发的远井事件点的密度都比较低，虽然B2井事件点的空间分布范围要小于B1井和B3井，但B2井的缝网体密度要大于B1井和B3井。

图  10  B平台3口水平井事件点的数量和密度分布

Figure  10.  Number and density distribution of event points for three horizontal wells on platform B

下载: 全尺寸图片 幻灯片

图11所示为B平台3口水平井分段压裂后150 d的累计产量。从图11可以看出，随着生产时间增长，3口水平井的累计产量趋于线性增加，表明压裂后3口水平井都形成了稳定的日产气量。对比发现B2井的产量要明显大于B1和B3井的产量，说明在同等压裂施工条件下，虽然B2井的缝网空间规模要小于B1井和B3井，但是B2井形成的缝网更为复杂，缝网体密度更大，总体缝网的渗流能力要优于B1井和B3井。

图  11  B平台3口水平井150 d的累计产量

Figure  11.  Cumulative gas production of three horizontal wells on platform B over 150 days

下载: 全尺寸图片 幻灯片

3.2   SRA模型与SRA模型的对比

通过微地震监测资料得到事件点集，利用SRV模型和上文SRA模型定量计算B平台3口井复杂缝网的SRV和SRA。为了提高SRV计算的准确性，分别采用常用的包络体模型、均匀网格模型以及新提出的非均匀网格模型计算B平台3口水平井的SRV^[23]，结果如图12所示。从图12可以看出，3种SRV模型计算出的SRV均满足B1井>B3井>B2井，说明B1井具有最大的改造体积，B3井次之，而B2井的改造体积最小。然而，从图11可以看出，B2井的产量最高，B3井的产量次之，B1井的产量最低，SRV的评价结果显然与实际产量结果不一致，甚至完全相反。从SRV的计算过程可知，SRV仅仅基于微地震事件点集的分布计算了压后缝网体的体积，虽然B1井和B3井相对B2井具有大的缝网体积，但是改造区域内事件点密度偏小，SRV中很大一部分的体积由远井区域的低密度事件点带来的，这些事件点很可能并不是由于水力缝网的扩展延伸导致的，而是压力波动对远井天然裂缝的激发导致的局部错动，并未形成具有渗流能力的缝网。虽然B2井的缝网体积最小，但是其具有最高的事件点密度和最大的缝网体密度，SRV只考虑了缝网体积，忽略了缝网体密度对产量的贡献。

图  12  B平台3口水平井SRV模型的计算结果

Figure  12.  SRV model calculation results of three horizontal wells on platform B

下载: 全尺寸图片 幻灯片

在对B平台3口水平井进行SRV评价之后，采用上文提出的基于八叉树网格缝网体面积计算方法计算SRA，其中活动节点阈值设置为1，计算结果如图13所示。从图13可以看出，B2井的SRA最大，B3井次之，B1井最小，与真实产量结果完全一致。相对于传统的SRV只能评价缝网体体积和规模对产量的贡献，SRA更能评价缝网体面积和质量对产量的贡献，得到的评价参数也更为可靠，与真实的产量结果也更加一致。

图  13  B平台3口水平井SRA模型的计算结果

Figure  13.  SRA model calculation results of three horizontal wells on platform B

下载: 全尺寸图片 幻灯片

4.   结　论

1）基于八叉树分解的SRA模型可以适用于不同数目和复杂度缝网的面积估算，模拟发现等效缝网面积曲线和真实缝网面积曲线高度一致，模型对于30条裂缝缝网的平均相对误差为12.31%，具有较高的计算精度。

2）为了分析模型的计算误差，基于蒙特卡罗方法和高斯混合模型开发了一套基于模拟事件点的验证流程，实现了缝网面积模型计算精度的定量评价，可进一步用于其他微地震解释模型的验证。

3）B平台3口井应用结果表明，与传统的SRV评价模型相比，SRA模型的计算结果更能够表征缝网体质量对产量的贡献，对于深层页岩密切割强加砂压裂效果评价和施工参数优化具有重要的价值。