支持向量机在煤层地应力预测中的应用

冯 鹏，李 松，汤达祯，陈 博，钟广浩

(1. 中国地质大学(北京) 能源学院，北京 100083；2. 煤层气开发利用国家工程研究中心中国地质大学(北京) 煤储层物性实验室，北京 100083)

0 引 言

地应力场的研究对煤层气开发意义重大[1-2]。由于煤层本身的低孔低渗特征，决定了其只有经过大规模压裂改造才能获得商业产能[3]，而地应力是影响压裂改造效果的关键因素，煤储层中割理(裂隙)系统的开启程度和增产改造压裂缝的形态扩展均受现今地应力场的控制[4-7]，针对煤层开展地应力研究是进行煤储层开采的必要环节。近年来，我国煤层气勘探开发在地质理论和开发技术方面已取得重大进展，但在煤层地应力方面的研究相对薄弱[8]。前人经过大量研究，提出了数十种地应力测量方法[9]。利用仪器测量得到的地应力比较准确，但受工程、技术及成本等因素的限制，特别是对深部地层而言，不可能进行大规模的地应力测量。为了降低成本，许多学者提出了地应力计算模型[10-14]，如Mattews和Kelly模型、Anderson模型、Mohr模型、Coulomb模型、黄氏经验模型、葛氏经验模型和组合弹簧模型等。还有学者结合技术手段对地应力进行预测：如黄继新等建立了应用成像测井资料进行地应力分析的方法[15]；刘景武、傅海成等基于多极子声波测井探测岩层裂缝、判断地应力方位等进而获得岩石力学性质参数[16-20]；楼一珊、林英松等首先利用声波测井资料计算出岩体完整系数，进一步基于静态力学参数计算了地应力的大小[21-22]；王晓锋、陈书平等利用有限元方法建立了煤储层的地质模型以及模拟了各地质时期的地应力分布特征[23-24]；魏忠文等利用构造应力场数值模拟方法计算了构造应力场的分布特征[25]；赵冰等将有限元线法同岩土工程反分析理论相结合，探讨了用岩土工程反演确定初始地应力和围岩参数方法的可行性[26]；张有天等根据少数测点数据及地表边界条件对地应力场进行趋势分析，提出了二维及三维地应力场的计算方法[27]；刘志强基于三轴压缩和凯撒实验，利用水力压裂和测井数据，建立了一套煤层地应力预测方法[28]。

虽然前人提出的地应力预测模型具有显著成效并得到了较好的应用，但对于煤层地应力的预测仍然存在一些不足。一方面，大多数地应力模型都是经验性的，存在一定的地域局限性；模型中还有大量参数需要利用其他方法或公式进行计算，过程繁琐，且参数的确定进一步增大了误差，实用性受到限制。另一方面，大多数的地应力预测研究都仅针对常规砂岩储层，而煤岩与砂岩存在显著的性质差异，导致上述方法不能直接应用在煤层上[28]。因此，建立一个简便且适用于煤层的地应力预测模型十分必要。由于测井资料具有测量深度较深、信息量大、数据相对连续等特点，基于测井资料计算地应力具有独特的优势[9]。近年来在地球物理反演中应用较广泛的支持向量机算法是一种可用于对数据进行高维特征分析的重要算法[29]。因此，本文基于测井资料，尝试将支持向量机算法应用于煤层地应力的预测中，并进一步借助三维地质建模软件完成地应力分布模型的建立，为解决煤层地应力的预测难题探索新的思路和方法。

1 实测地应力计算方法

1.1 水平主应力计算方法

根据地面垂直钻孔水力压裂测量地应力方法[30]，闭合压力Pc为最小水平主应力σh，即

σh=Pc

(1)

最大水平主应力σH的计算公式为：

σH=3Pc-Pf-Po+T

(2)

式中：T为煤或岩石的抗拉强度，MPa，本次计算中取定值0.8；Pf为煤层破裂压力，MPa；Po为储层压力，MPa。

基于鄂尔多斯盆地典型煤层气区块注入/压降试井数据，利用以上公式求出煤层的最大水平主应力σH和最小水平主应力σh，并对二者进行了相关性分析(图1)。最大与最小水平主应力之间呈正相关关系，且相关性较强(R2=0.93)，故在已知最小水平主应力σh的情况下，可利用如下模型计算最大水平主应力σH：

σH=1.7394σh-3.4279

(3)

图1 最大与最小水平主应力拟合关系图Fig.1 Fitting relation diagram of maximum and minimum horizontal principal stress

1.2 垂向主应力计算方法

垂向主应力σv是由上覆岩石的自重引起的，所以计算公式为：

(4)

其中：ρ为地层密度，kg/m3；g为重力加速度，取9.8 m/s2；h为深度，m。

实际计算时简化积分，可以用如下公式计算[31]:

σv=0.027×h

(5)

2 最小水平主应力预测方法

2.1 基于SVM回归的煤层地应力预测原理

支持向量机(SVM)是Cortes和Vapnik于1995年首先提出的[32]，它在解决小样本、非线性及高维模式识别中有独特优势[33]。研究所用的支持向量回归模型，其基本思想包括升维和线性化。即当数据集Xi在原始特征中不是线性可分时，通过非线性函数(Xi)把原始特征空间映射到一个高维特征空间，SVM算法将在这个高维特征空间中找到一个函数间隔最大的线性超平面来对数据进行分类。最后把高维的分类面映射回低维，就可以在低维空间对数据分类(图2)。

图2 SVM结构示意图Fig.2 Structure schematic diagram of the Support Vector Machine (SVM)

在高维特征空间中构造最优超平面只需要计算支持向量与特征空间中向量的内积，而SVM不能直接计算该内积，因此需要用满足Mercer定理的核函数K(x,y)来代替[34]。本文选取径向基函数作为核函数：

(6)

图3 最小水平主应力预测流程图Fig.3 Prediction flow chart of minimum horizontal principal stress

2.2 煤层地应力主要关联因子的筛选

测井数据类型多样，为了避免引入某些噪声信息，在建立预测模型之前，需要对输入的测井数据类型进行筛选。灰色关联分析方法是一种多因素统计分析方法，以各因素的样本数据为依据，用灰色关联度描述因素间关系的强弱、大小和次序[35]。前人提出了多种关联度模型，本文采用邓氏关联度模型进行分析，Xi与X0的关联度可表示为：

(7)

其中：

(8)

通过对鄂尔多斯盆地东缘典型煤层气区块的测井和注入/压降试井数据的整理，将包含9种不同类型测井数据(密度、自然伽马、中子、自然电位、声波时差、两个方向井径、深侧向和浅侧向电阻率)及对应的最小水平主应力的34个样本数据(表1)导入Matlab软件中，作为34×10的矩阵进行分析。由于各组数据的量纲不同，计算关联系数之前，先将数列进行初值化处理。通过公式(7)计算得到如下结果(图4)：声波时差(DT)和自然电位(SP)测井值与最小水平主应力关联度相对较低，中子测井(CNL)与最小水平主应力关联度最高。以关联度大于0.7为界，可得出影响最小水平主应力的主要因子为两个方向的井径测井(CAL1和CAL2)、深浅侧向电阻率测井(LLD和LLS)、补偿中子测井(CNL)、自然伽马测井(GR)和密度测井(DEN)共7个因素。由于CAL1和CAL2以及LLD和LLS的值基本一致，分别与最小水平主应力的关联度也基本相等，为了进一步减少影响因子的数量，以CAL1和CAL2的平均值(CAL)作为一个影响因子，以LLD和LLS的平均值(R)作为一个影响因子。由此得到最小水平主应力的5个主要影响因子：CAL、R、CNL、GR和DEN。

图4 各类测井数据的关联度柱状图Fig.4 Correlation histogram of various types of logging data

2.3 基于SVM的回归模型的建立

按照预测样本和训练样本1∶3的比例从34个样本中随机选取了9个预测样本，其余25个样本用于训练(表1)。其中井径测井(CAL)、补偿中子测井(CNL)、自然伽马测井(GR)、密度测井(DEN)和深浅侧向电阻率测井值的平均值(R)作为输入数据，最小水平主应力作为输出数据。

在对所有数据利用离差标准化方法进行归一化处理之后，将输入数据和学习目标导入编写好的程序中，对样本进行多次迭代训练。通过不断调整，最终得出参数c和μ的值分别为20.06和1.94时模型最优。对模型进行检验，实测值(由式(1)计算所得)与预测值平方相关系数可达到0.94，平均相对误差为0.07，具有较好的预测能力(图5)。

3 实例应用

韩城区块位于鄂尔多斯盆地东缘南部，渭北煤田东北缘，整体表现为北东走向、北北西倾向的单斜构造。位于韩城区块的H3井组煤层埋深跨度较大，受构造活动影响较大，导致煤层的地应力场在平面和垂向上均存在较大差异，而地应力是煤层气开采的一个重要的影响因素[36]。因此，开展该区H3井组11号煤层地应力的预测，同时借助Petrel 2015软件完成该区地应力场的三维地质建模，厘定研究区煤层地应力在平面和垂向上的展布特征，为该区煤层气勘探开发提供理论依据。

表1 训练样本和测试样本原始数据

图5 最小水平主应力实测值与预测值对比图Fig.5 Comparison between actual and predicted values of the minimum horizontal principal stress

3.1 研究区地应力分布特征

三维地质建模是用模型来表征储层和物性参数的三维空间分布以及变化规律的一种手段[37-39]。本次研究基于SVM回归模型计算韩城区块H3井组11号煤层最小水平主应力，进而利用最小水平主应力和煤层埋深分别计算得到最大水平主应力和垂向应力。以计算得到的地应力数据作为基础数据，运用序贯高斯模拟法构建煤层地应力场的地质模型，厘定煤层地应力场的三维展布特征。计算表明，研究区11号煤层最小水平主应力为4.86～27.12 MPa，最大水平主应力为5.02～43.74 MPa，垂向主应力15.93～24.96 MPa。陈世达等通过232组实测地应力数据统计分析得出[2]，我国煤储层最小水平主应力为2.14～30.26 MPa，最大水平主应力为2.70～45.20 MPa，垂向主应力3.67～37.99 MPa。由此进一步说明，本文通过支持向量机回归建立的最小水平主应力预测模型，反演得到的各应力大小均在合理范围。

地应力场三维地质建模结果显示，从平面展布来看，研究区水平最小主应力的高应力范围较小，且主要集中在研究区东部，西部次之，研究区中部和北部为水平主应力的低应力区；从垂向分布来看，该煤层顶底板附近的最小水平主应力整体偏大，而煤层中部最小水平主应力较小(图6(a))。研究区最大水平主应力平面展布和垂向分布规律与最小水平主应力基本一致(图6(b))。垂向应力在平面展布上分布较为均匀，高应力区以西南部为主，低应力区以东北部为主，这是由煤层东北高、西南低的地势以及区域地质构造所造成；垂向上由于煤层较薄，没有明显变化规律(图6(c))。

图6 韩城矿区H3井组11号煤层地应力分布栅状图(地应力单位：MPa)Fig.6 Grid diagram of in-situ stress distribution in No.11 coal seam of well group H3, Hancheng mining area

3.2 不同深度的煤层地应力分析

为了分析研究区不同深度煤层的应力场分布状态，基于前文建立的地应力预测模型，反演获得研究区295个地应力数据点进行应力与埋深相关性分析(图7)。

图7 韩城矿区H3井组煤层地应力与埋深拟合关系Fig.7 Fitting relations between in-situ coal seam stress and burial depth of H3 well group in Hancheng mining area

图8 韩城矿区H3井组煤层侧压系数与埋深拟合关系Fig.8 Fitting relations between coal seam lateral pressure coefficient and burial depth of H3 well group in Hancheng mining area

由地应力随埋深的变化可以看出，该矿区H3井组煤层埋深为500～950 m，煤层所受三个方向的主应力都随着埋深的增加而增大，而垂向应力随埋深的增加而增加的幅度明显大于水平主应力，三个方向的应力之间的相互关系随着埋深不同也在发生变化，即在垂向上地应力场的类型也存在分异(图8)：

(1)500～600 m，以σH>σv>σh型应力场为主，最小水平主应力介于9.32～19.28 MPa之间，平均为13.10 MPa，最大水平主应力占主导作用，地应力表现为大地动力型。水平方向与垂直方向的应力大小关系可进一步用侧压系数来表示，侧压系数即最大与最小水平主应力的平均值与垂直主应力的比值[4]。该深度范围内侧压系数较大，煤储层主要处于挤压带，部分处于过渡带。

(2)600～750 m，以σH≈σv>σh型应力场为主，为地应力场转化的过渡阶段，各地应力间大小关系变化较为频繁。在此深度范围内，大地动力场型和大地静力场型两种应力场共存。该深度范围煤储层主要处于过渡带，部分处于伸张带。

(3)750～950 m，以σv>σH>σh型应力场为主，最小水平主应力介于7.05～22.76 MPa之间，平均为13.85 MPa，垂向主应力占主导作用，地应力表现为大地静力型。该深度范围煤储层主要处于伸张带，部分处于过渡带。

4 结 论

(1)基于灰色关联法模型分析，井径测井(CAL)、补偿中子测井(CNL)、自然伽马测井(GR)、密度测井(DEN)和深浅侧向电阻率测井值的平均值(R)与煤层最小水平主应力关联度最好。以这5个测井参数作为训练因子，利用基于小样本理论的支持向量机回归方法建立煤层最小水平主应力预测模型，实测值与预测值平方相关系数可达到0.94，因此该模型可以较为高效、精准地对最小水平主应力进行预测。

(2)以鄂尔多斯盆地东缘韩城区块H3井组煤层为研究对象，通过建立煤层地应力场三维地质模型以及地应力与埋深的拟合关系综合分析，研究区煤层埋深范围在500～950 m，地应力随埋深的增加呈现递增的趋势。煤层埋深处于500～600 m时，以σH>σv>σh的大地动力型应力场为主；埋深处于600～750 m时，以σH≈σv>σh型应力场为主，为地应力场转化的过渡阶段；埋深处于750～950 m时，以σv>σH>σh的大地静力型应力场为主。总体而言，随着埋深的增加，煤储层所处的应力环境由挤压带过渡为伸张带。