基于机器学习方法的海洋天然气水合物水平井降压开采模拟—优化耦合模型

2020-09-08王海彬罗建男袁益龙夏盈莉朱慧星

天然气工业 2020年8期

辛欣王海彬罗建男于涵袁益龙夏盈莉朱慧星陈强

1.吉林大学地下水资源与环境教育部重点实验室 2.自然资源部天然气水合物重点实验室·中国地质调查局青岛海洋地质研究所 3.青岛海洋科学与技术试点国家实验室海洋矿产资源评价与探测技术功能实验室

0 引言

天然气水合物（以下简称水合物）以填充、胶结等多种形式赋存于沉积物中，通过降压开采打破水合物形成的压力平衡条件，可以使固态水合物分解为气态甲烷和液态水[1]，从而使含水合物沉积层的抗剪强度和承载能力降低[2]，同时，降压开采使含水合物沉积层压力降低，有效应力增加。基于上述原因，降压开采易导致含水合物沉积层变形，从而产生海底面沉降、变形等地质灾害[3-4]。因此，在水合物的开采过程中，如何保证在获得较高产气量的同时，将海底面沉降量控制在一定的安全范围内，成为实现水合物安全、高效规模化开采的前提条件。

目前国内外学者多运用数值模拟的方法来评价不同开采方案下含水合物沉积层的产气能力及海底面沉降量。Rutqvist等[5]在水合物数值模拟程序TOUGH+HYDRATE上链接岩土力学软件FLAC3D，针对加拿大Mallik冻土区和美国阿拉斯加冻土区模拟预测了降压开采水合物过程中的地层力学响应，结果表明影响含水合物沉积层剪切破坏程度的主控因素是降压幅度。Kimoto等[6]利用化学—热学—力学耦合模型，模拟开采水合物引起的沉积层力学变形过程，结果表明地层变形是基于水合物相态变化引起的松散沉积物应力强度降低。Gupta等[7]将含水合物沉积层力学场变化考虑为线弹性变形，基于有限元法建立了模拟水合物开采的流动、传热、力学多场耦合模型，并运用多个解析解和实验数据对该模型进行了验证。Jin等[8-9]利用自主开发的HYDRATEBiot程序模拟预测了南海神狐海域水合物在降压开采过程中的产气能力及海底面沉降量，结果表明在降压开采初期，海底面发生快速沉降，而在开采后期，沉降速率变慢且趋于平缓。

上述研究均采用数值模拟手段预测了既定开采方案下含水合物沉积层的力学响应特征及其对甲烷累计产气量的影响，但无法解决考虑含水合物沉积层力学变形影响下的水合物开采方案优选问题。而将数值模拟技术与运筹学中的优化技术相结合，构建模拟—优化耦合模型，既能预测水合物开采引起的多相流体渗流过程和温度场、压力场、应力场的变化趋势，又能够在满足各种环境、经济、技术等要求的前提下获得最优的开采方案。因此，将数值模拟模型（以下简称模拟模型）与目标优化模型（以下简称优化模型）耦合的技术（Simulation and Optimization Technology）[10]为解决以水合物安全、高效开采过程为目标的开采方案优选问题提供了科学途径。而目前关于应用模拟—优化耦合技术进行水合物安全、高效开采方案优选的报道还鲜见。

为了找到甲烷累计产气量最优值与地层稳定性的关系，基于机器学习方法形成模拟—优化耦合技术，构建起水合物降压开采传热—流动—力学数值模拟模型、可以替代数值模拟模型的机器学习模型和以甲烷累计产气量最优为目标的混合整数非线性规划优化模型；在此基础上，选取南海神狐海域厚层Ⅱ类水合物藏W11站位为研究对象，获得了海底面沉降量约束下的水合物储层甲烷累计产气量以及对应的最优开采方案参数，以期为水合物安全、高效开采方案的制订提供支撑。

1 基于机器学习方法的模拟—优化耦合技术原理

模拟—优化耦合技术包括模拟模型和优化模型两个部分，其中模拟模型是核心，能够描述真实地质体输入、输出响应关系。建立水合物开采模拟模型包含以下3个步骤：①建立描述水合物储层空间结构的概念模型；②建立水合物传热—流动—力学耦合数学模型；③运用数值模拟程序HYDRATEBiot求解。优化模型是基于运筹学理论的最优方案求解模型，能够突破既定开采方案的限制并且求得全局最优解[10-11]。优化模型的构建方法包括线性规划、非线性规划、整数规划等多种方法[10-11]。水合物开采引起的多相流体渗流及地层力学特征的响应过程具有典型的非线性特征，因此运用混合整数非线性规划方法来构建优化模型，然后采用遗传算法进行求解。

将模拟模型作为一个目标函数或约束条件编入优化模型中，在利用优化模型求解全局最优解时需要反复多次调用模拟模型，模拟模型的非线性特征越强，调用计算的负荷越大，严重时甚至出现无法运算的情况[12]，从而限制了模拟—优化耦合技术的应用。近年来，机器学习的理念逐渐渗透进入越来越多的研究领域，通过学习模拟模型计算得到的输入—输出特征响应关系，建立起模拟模型的替代模型，而该替代模型具备与模拟模型相同的计算能力，此即为采用机器学习方法的主要思路。建立了替代模型以后，优化模型对模拟模型的调用即可以更改为对替代模型的调用，由于替代模型结构简单、求解容易，从而解决了计算负荷大、求解难度高等问题。替代模型的本质是一种机器学习方法，常用的建模方法包含多项式回归法、人工神经网络法及克里格法等[13]。

综上，模拟—优化耦合技术的运用不但能够考虑复杂的水合物分解行为（包含多相流体渗流行为、地层变形行为等），对于给定的约束条件（如海底面沉降量、开采时间、甲烷累计产气量等），还能够反演求得既定目标下的最优开采方案。因此，在考虑影响地层力学稳定性的限制条件下，采用模拟—优化耦合技术可以计算得到最优累计产气量，主要步骤分为以下3步：①建立水合物开采模拟模型并求解；②对水合物开采模拟模型进行机器学习，即构建替代模型；③建立优化模型并求解。

2 水合物开采地层力学响应模拟模型

2.1 概念模型与网格剖分

南海北海神狐海域是我国海域水合物勘探开发的重点靶区[14]。该区域在新生代时期地层沉积速率较大，最大沉积厚度达到11 km，油气资源丰富[15]。多期次构造运动使该区域具有泥底辟活动带，在其上覆地层中形成高角度断裂和垂向裂隙，为深部流体提供了流动通道，进而形成了特有的水合物成藏系统[16-18]。由GMGS3钻探航次调查结果可知，W11站位海底面水深为1 312 m，水合物储层最大厚度达80 m，平均水合物饱和度（SH）为30%，甲烷气体体积分数超过99%。水合物储层孔隙度为50%，渗透率为10 mD，地温梯度为54.9 ℃/km。水合物储层下伏地层无明显的游离气层特征，为Ⅱ类水合物藏[19]。

为求得W11站位水合物开采过程中的累计产气量及地层力学响应特征，针对该站位建立水平井降压开采概念模型。含水合物沉积层为水平延伸地层，由于传热—流动—力学耦合数学模型的求解对计算机内存消耗大，且在降压开采过程中沿整个水平井筒的甲烷气体压力损失较小。因此考虑沿水平井段的网格尺寸为1 m。设置水合物储层厚度为80 m，其上覆弱透水地层厚度为200 m，其下伏饱水地层厚度为20 m（图1），模型顶板即为海底面。模拟区域在x方向上长度为500 m，网格尺寸随着与生产井距离的增加而增加，分别为 1 m、2 m、3 m、4 m 和 5 m ；z方向上高度为300 m，针对水合物层划分的网格尺寸为2 m，下伏饱水地层的网格尺寸也为2 m，上覆地层临近水合物层的20 m范围内网格尺寸为2 m，其余区域则为4 m。

图1 W11站位天然气水合物水平井开采传热—流动—力学耦合概念模型示意图

2.2 初始条件与边界条件

初始温度场由海底面温度（4.82 ℃）按地温梯度（5.49 ℃/100 m）逐渐增加[19]，水合物储层底部初始温度为15.373 ℃、初始压力为16.195 MPa，此时水合物储层中固态水合物和液态水处于相平衡状态，其中水合物饱和度为30%、含水饱和度为70%。上覆、下伏地层初始含水饱和度为100%，盐度为3%。饱和水及不含水的含水合物沉积层导热系数分别为3.1 W/（m·K）、1.0 W/（m·K）。

由于缺乏研究区相关站位的原位地应力数据，模型中初始有效地应力场由储层骨架颗粒密度（2 600 kg/m3）和深度计算获取，同时假设在x、y、z三个方向上初始地应力相等[2]。

将模型顶底板（其中模型顶板即为海底面）设定为可发生流体运移和热量交换的定压定温边界，其温度和压力值均为未开采前的初始值。考虑模型的对称性，侧向边界设定为隔水隔热边界。水平井为模型内边界，内边界压力值即为水平井开采压力值。

含水合物沉积层的相对渗透率采用Stone模型，如式（1）、（2）所示；毛细管压力的计算采用van Genuchten模型，如式（3）所示。地质力学特征参数如表1所示[20-21]。

液相相对渗透率（KrA）计算式为：

气相相对渗透率（KrG）计算式为：

式中SA表示液相饱和度；SirA表示残余水饱和度，取0.60；nA表示液相衰减指数，取4.5；SG表示气相饱和度；SirG表示残余气饱和度，取0.02；nG表示气相衰减指数，取3.5。

毛细管压力（pc）计算式为：

式中pc表示毛细管压力，最大毛细管压力为2.0×106Pa；p0表示毛细管进气压力，取 1.0×104Pa；m表示van Genuchten指数，取值为0.45。

表1 含水合物沉积层地质力学特征参数表

2.3 模拟模型的样本空间

确定模拟模型样本空间的目的是为以建立替代模型为目的的机器学习过程提供训练样本，样本空间的覆盖性越强、样本数量越多，训练效果则越佳，从而使替代模型具有与模拟模型相似的输入输出运算关系。

为确定不同水合物开采方案下甲烷累计产气量与地层力学稳定性的平衡关系，以甲烷累计产气量最佳为优化目标，以海底面沉降量小于某限值为约束条件，以降压幅度、开采时间、布井位置（水平井段在垂向上的位置）、水平井段长度为自变量，构建样本空间。依据水合物试采工程经验来确定自变量取值（表2）[22]，得到样本空间样本总数为750组（自变量样本数量的积）。

2.4 模拟模型求解

水合物传热—流动—力学耦合数学模型采用自主开发的数值模拟程序HYDRATEBiot求解，该程序是基于Biot固结理论建立的TOUGH+HYDRATE软件框架的扩展程序。HYDRATEBiot既可以求解水合物分解引起的多相多组分系统中复杂的相态变化、流体输运、热量传输等过程，还可以求解考虑了孔隙水压力和岩土体变形相互作用的含水合物沉积层力学特征响应过程[23]。此次采用该程序求解W11站位含水合物沉积层水平井降压开采时的甲烷累计产气量及地层力学响应特征。

表2 样本空间自变量参数取值范围统计表

2.5 模拟结果分析

2.5.1 降压幅度和布井位置对单位长度水平段甲烷累计产气量的影响

依据前述样本空间自变量的参数取值范围，设定开采时间为180天，考虑降压幅度和布井位置变化，研究单位长度水平段甲烷累计产气量和地层力学响应特征。由模拟计算结果可知，水平井降压幅度为10 MPa，连续开采180天，不同布井位置下单位长度水平段甲烷累计产气量分别为6.75 m3、7.82 m3和6.28 m3。在设定的开采时间（180天）内单位长度水平段甲烷累计产气量随降压幅度升高而升高，且水平段布置在储层中下部可以获得相对较高的累计产气量（图2），这受控于下伏地层良好的热动力学条件与水合物储层的自封闭效应。一方面，对于厚层Ⅱ类水合物藏而言，水平井段在水合物储层垂向上的位置越靠近水合物层底界（BSR），其热动力学条件更有利于水合物发生分解；但另一方面，对于Ⅱ类水合物藏并非将水平井布置在邻近BSR位置处最佳，在开采过程中还需考虑下伏地层涌水的影响。随着开采进行，水平井下方水合物储层中水合物逐渐分解，水平井与下伏透水地层取得良好的水力连通后将引起下伏地层水大量涌入水平井，从而导致水合物储层内部的压力降受到抑制，使水合物分解速率降低，不利于水合物的高效开采。因此，在以较大降压幅度进行水合物开采时，水平段部署在储层中下部为宜。

图2 不同布井位置、降压幅度下单位长度水平井段甲烷累计产气量柱状图

2.5.2 海底面沉降量

水平井在Ⅱ类水合物藏中降压开采30天时海底面沉降量与开采180天的总沉降量的占比介于15%～20%；开采60天时，不同开采制度下的海底面沉降量与总沉降量的占比介于40%～50%（图3）。由此可见，水平井降压开采Ⅱ类水合物藏时地层沉降主要发生在开采早期阶段，地层的沉降主要是由于孔隙压力降低导致储层骨架压缩，这是深部地层内部发生的一种体积应变。因此在降压开采早期阶段应采取必要的措施来防止井筒周围地层产生强烈的变形作用进而导致严重的工程问题出现。

海底面沉降量随降压幅度增大而增大，以生产井位于水合物储层中部为例，当降压幅度依次为4 MPa、6 MPa、8 MPa、10 MPa、12 MPa，相应最大海底面沉降量分别为 0.035 m、0.060 m、0.083 m、0.101 m、0.118 m（图3）。假定海底面沉降0.1 m为开采水合物保持地层长期稳定的临界条件，那么针对属于Ⅱ类水合物藏的W11站位，降压幅度超过10 MPa的强降压方案不能保证水合物储层的长期、稳定开采。进行开采方案设计时，应结合考虑采取适当的措施来避免较大的海底面沉降量产生，同时需做好海底变形的多场监测。

图3 不同开采时间下海底面沉降量占比及不同降压幅度下最大海底面沉降量统计图

3 替代模型的建立与验证

3.1 替代模型的建立

采用径向基函数人工神经网络（Radial Basis Function Artificial Neural Network,简写为 RBFANN）方法建立模拟模型的替代模型。径向基函数是对中心点径向对称衰减或增强的非线性高斯函数。RBFANN是一种前馈式神经网络，包括输入层、隐含层和输出层，从输入层到隐含层为非线性变换，从隐含层到输出层为线性变换[24]。总体来说，RBFANN模型是一个黑箱模型，即不需要考虑复杂的物理过程，仅基于输入输出样本数据集来构建内部映射过程，因此其特点是运算速度快，对于高度非线性问题具有较强拟合能力[24]。由此，RBFANN方法可以用于建立描述水合物开采的多相渗流非线性数学模型的替代模型。

利用前述750组样本的输入输出数据作为替代模型的数据库，随机抽取600组为训练样本、150组为检验样本，通过MATLAB基础平台编写径向基函数调用程序，对径向基函数人工神经网络进行训练。在训练过程中，将降压幅度、开采时间、水平井段长度、布井位置作为替代模型的输入变量，再分别将甲烷累计产气量和海底面沉降量作为替代模型的输出变量，进而建立起两个替代模型，其中替代模型2的输入变量不包括水平井段长度（表3），两个模型训练目标绝对误差均为0.001。

表3 替代模型参数统计表单位：个

3.2 替代模型的验证

将基于RBFANN方法建立的替代模型的计算结果与150组检验样本进行对比，如图4、5所示，替代模型的计算精度较高。

运用R2（确定性系数）、MRE（平均相对误差）两个指标对替代模型的计算精度进行检验，其计算式如式（4）、（5）所示。R2数值越接近于1，表示替代模型的计算精度越高；MRE数值越接近于0，表示替代模型的计算精度越高。由检验结果可知，两个替代模型的计算精度均较高（表4）。

图4 替代模型与基于检验样本的模拟模型计算结果对比图（甲烷累计产气量为输出变量）

图5 替代模型与基于检验样本的模拟模型计算结果对比图（海底面沉降量为输出变量）

式中n表示样本个数；i表示样本序号；yi表示模拟模型响应值；表示替代模型响应值；表示模拟模型响应值的平均值。

表4 替代模型计算精度统计表

4 混合整数非线性规划优化模型的建立与求解

4.1 优化模型的建立

4.1.1 目标函数和决策变量

设置甲烷累计产气量为目标函数，降压幅度、开采时间、水平井段长度、布井位置为决策变量，则目标函数式为：

式中Q表示甲烷累计产气量，m3；L表示水平井段长度，m；Z表示与布井位置相关的整数变量（取值0或1）；p表示降压幅度，MPa；t表示开采时间，d。

4.1.2 约束条件

约束条件包括每个决策变量的约束条件和附加约束条件，共考虑5个约束条件。

将海底面沉降量作为约束条件，其条件式为：

式中g表示海底面沉降量，m；M表示最大允许沉降量，m。

将水平井段长度作为约束条件，其条件式为：

式中L0表示水平井段长度最大值，取值为500 m。

将降压幅度作为约束条件，其条件式为：

式中p0表示降压幅度的最大值，取值为12 MPa。

将开采时间作为约束条件，其条件式为：

式中t0表示开采时间最大值，取值为180 d。

将布井位置作为约束条件，共考虑3种布井位置，Z1代表水平井段位于储层中部，Z2代表水平井段位于储层中下部，Z3代表水平井段位于储层下部，上述3个变量取值为1，表示井存在，若取值为0，则表示井不存在，满足的条件式为：

将式（6）～（11）联立，则建立起优化模型，该优化模型具有两个特点：①包含目标函数、决策变量、约束条件的两个替代模型；②约束条件中含有整数变量。

4.2 优化模型的求解

基于MATLAB平台编写遗传算法求解程序，求解前述混合整数变量的非线性规划优化模型。求解步骤为：①选用二进制形式进行编码；②设定遗传算法参数（表5）；③通过选择—交叉—变异—再选择过程的反复迭代，使个体适应度不断提高，最终得到满足所有约束条件的最优解。在编码、选择、交叉、变异各过程中选用了不同的策略，如采用二进制形式进行编码，选择过程采用随机遍历抽样，交叉方法采用单点交叉，采用二进制形式进行变异。基于上述策略，通过提高初始种群数和遗传代数来保证遗传算法的全局寻优能力。通过求解优化模型，即可获得在海底面沉降量约束下的水合物储层甲烷气累计产气量，以及对应的最优开采方案参数（表6）。

表5 遗传算法参数统计表

4.3 优化模型计算结果与分析

如图6所示，随着最大允许沉降量数值增大，甲烷累计产气量也增大，二者满足正相关关系；当最大允许沉降量介于0.04～0.05 m时，是开采方案优化结果的转折点。当最大允许沉降量较小时（小于0.045 m），影响甲烷累计产气量和海底面沉降量的决策变量是开采时间，海底面沉降量随开采时间增长而增大；当最大允许沉降量较大时（大于0.045 m），影响甲烷累计产气量和海底面沉降量的决策变量是降压幅度，海底面沉降量随降压幅度增大而增大。由此可知，因水合物开采引起海底面沉降主要发生在开采初期，如期望获得较高甲烷气产量同时控制最大海底面沉降量在0.05 m以下，在开采初期应该减小降压幅度。

表6 优化模型计算甲烷累计产气量与最优开采方案参数统计表

图6 甲烷累计产气量与海底面沉降量平衡关系曲线图

以甲烷累计产气量最高为目标，并且以最大允许沉降量为约束条件，优选水平井段部署在水合物储层下部，这与模拟模型的计算结果（水平井段部署在水合物储层中下部）不一致。原因在于模拟模型是在既定方案下，以甲烷累计产气量最高为单一目标的求解结果，没有考虑海底面沉降量的约束。另外，与模型设置的开采时间也有关，模拟模型和优化模型的水合物开采时间均为180天，如果生产时间再延长，布井位置的优选结果还会发生变化。

将表6列出的开采方案参数代入到前述模拟模型中进行求解，对比结果可知，优化模型的求解结果与模拟模型的计算结果相对误差均低于±5%（表7），表明优化模型的计算结果可靠。