丁建新 李雪松 宋先知 张诚恺 马宝东 刘子豪 祝兆鹏
(1.昆仑数智科技有限责任公司 2.中国石油大学(北京)油气资源与工程全国重点实验室 3.中国石油大学(北京)石油工程学院)
近年来,我国油气对外依存度持续处于高位,2022年石油和天然气对外依存度分别达到71.2%和40.5%,远超国际能源安全预警标准[1]。因此,加大油气勘探开发力度,保障国家能源安全,是我国油气行业未来发展的主要任务。钻井工程是油气资源勘探开发的关键环节,保证安全高效钻进是实现油气勘探开发、降本增效的重要保障。其中钻井过程包含钻进破岩、管柱延伸和井筒流动等多个复杂过程[2-4],国内外专家学者在钻速优化、摩阻计算和井筒流动等多个方面进行了深入研究。
钻井过程中的机械钻速是反映钻井效率的关键指标,而准确预测钻速是破岩提速的基础。目前钻速预测主要有统计回归和智能预测2种方法。国外学者基于室内试验或者现场数据,通过分析钻具、工程和地层等各类因素对机械钻速的影响机理,构建了Maurer、Warren、Bourgoyne-Young、Hareland钻速方程以及修正的杨格钻速方程等多种方程。国内学者同样基于多元回归方法分析求解钻速方程,为钻速预测提供了理论支撑[5]。随着人工智能技术与传统油气行业的深度融合,钻速智能预测方法逐渐成为行业研究热点,支持向量机[6]、随机森林[7]等机器学习方法为钻井工程参数与钻速之间映射关系的构建提供了新的方法和思路,而全连接神经网络(FCNN)[8]、循环神经网络(RNN)[9]等深度学习模型借助其非线性建模能力和高维特征分析能力,在精度、泛化性提升等方面已取得显著效果。
管柱摩阻是评价水平井钻井风险的重要因素,摩阻过大会降低钻压、扭矩传递效率,影响岩石破碎效果,严重时可能导致卡钻等风险。摩阻扭矩的传统计算方法主要有软杆模型[10]和刚杆模型[11],分别适用于直井和水平井。有限元、有限差分、加权残值等方法也逐步应用到管柱摩阻扭矩相关研究中,为管柱下入能力评价及卡钻分析提供理论基础。近年来,有学者采用回归分析、随机森林等机器学习算法[12]对摩阻扭矩进行实时预测,利用改进的BP神经网络[13]建立摩阻系数与影响因素之间的隐含关系,实现了摩阻系数的高效计算。例如,李紫璇等[14]考虑了岩屑床对管柱摩阻扭矩的附加影响,提高了卡钻预测精度。
井眼清洁是影响水平井钻井风险的另一个关键因素。钻井过程中产生的岩屑在重力作用下极易在斜井段、水平段沉积并形成岩屑床,掩埋管柱,导致钻杆扭矩增大[15]。因此学者们开展了相关试验,研究环空流速、钻井液密度和流变性、岩屑粒径、机械钻速等参数对岩屑床高度等井眼清洁关键参数的影响规律,随后提出了经典的分层模型和临界速度模型。例如:马志忠等[16]探索了井眼清洁计算新方法,提高了环空岩屑浓度和井筒ECO的计算精度;张菲菲等[17]结合传统水力学模型和人工智能方法,提出了数据驱动定量评价岩屑动态分布的新方法,为解决大位移井井眼清洁不充分提供了详细井下信息。
综上,前人针对破岩提速、管柱减阻和井筒清洁3个方面开展了大量研究,为钻井效率优化和风险管控提供了坚实的理论基础。但是,钻井过程中这3个子系统相互耦合、相互制约,优化钻井过程需要综合考虑风险和效率,现有单过程、单目标优化方法难以取得全局最优效果[18-19]。例如,提高钻速会增大岩屑量,岩屑床堆积过高可能会增大钻柱摩阻造成卡钻,也会降低钻压导致钻速下降。因此在钻井多过程优化方面,挪威eDrilling公司集成了水力学、管柱力学、机械钻速等多个模块,搭建了钻井数字孪生系统[20-21],实现了钻井全生命周期的动态监测与参数优化。挪威NORCE研究所构建了自动化钻井系统[22-23],能够实现钻井优化、轨迹调控和风险处理等功能。尽管上述工作考虑了部分钻井过程子系统的耦合,但在协同钻速、摩阻和岩屑运移等方面仍然有待深入研究。
为此,本文提出一种考虑破岩提速-管柱减阻-井筒清洁3个子系统协同作用的水平井钻井多目标优化方法。该方法以提高机械钻速、降低机械比能为目标,以钻压、转速和排量为调控参数,以减小管柱摩阻、促进井筒清洁为约束条件,通过耦合钻速预测、摩阻系数反演和岩屑运移模型模拟,实现水平井钻井提速-减阻-清屑多目标协同优化,为现场安全高效钻进提供保障。
本文选择某油田1口水平井进行实例分析。该井井深7 788 m,造斜点位于井深7 105 m,表1列举了部分测井和录井参数。在数据处理方面,通过KNNImputer函数[24]计算各数据的欧几里德距离,寻找最近邻样本,使用最近邻样本非空数值的均值来填补缺失值;基于箱线图法[25]计算中位数和四分位数,获取数据正常值的上、下边界,并以此为依据剔除上、下边界外的异常值。
在钻速智能预测方面,神经网络、随机森林和支持向量机是当前研究最充分、使用最广泛的模型之一[26-27]。其中,随机森林[28]是一个集成了许多决策树的智能学习机器,模型示意图见图1。
决策树是相对简单的决策模型,具有类似流程图的树状结构,每棵树都由不同节点(根节点、内部节点和叶节点)组成,可将输入数据分解为更小的子集。作为一个包含许多决策树的集成学习器,随机森林通过综合所有决策树的预测结果来给出输出。与神经网络和支持向量机相比,随机森林的一个优点是预测结果更稳健,输出会始终处于训练标签的范围内,不会像神经网络一样给出不合理的异常预测结果[29]。此外,C.SOARES等[30]、ZHANG C.K.等[31]分别对比了3种模型在钻速实时预测与优化方面的准确性和稳定性,认为随机森林优于神经网络和支持向量机。故本文采用随机森林预测机械钻速。
模型超参数是指在模型训练开始前应设置的参数。例如,在训练随机森林之前,应人工确定树的数量、每棵树的最大深度等。人工智能模型的性能在很大程度上依赖于模型的超参数设置,因此,需要在模型开始训练前进行超参数调优[32]。网格搜索法和随机搜索法是2种常见的超参数调优方法[33]。其中,网格搜索法遍历比较模型在所有超参数组合下的性能,其结果最为准确,但计算代价高昂。例如,要为具有m个超参数且每个超参数都具有n个候选值的模型选择最佳超参数时,模型必须运行mn次。相反,随机搜索法仅从所有超参数组合中随机选择若干进行比较,与网格搜索法相比,在效果相似的情况下可大大缩短训练时间。因此,本文采用随机搜索法对随机森林的超参数进行自动调优。表2列出了模型超参数及其取值范围[30]。
表2 随机森林模型超参数及其取值范围Table 2 Hyperparameters of random forest model
管柱摩阻扭矩计算机理模型主要有软杆模型和刚杆模型,其中刚杆模型相较于软杆模型考虑了钻柱刚度,更适合大位移井、水平井的摩阻扭矩计算,因此本文采用刚杆模型进行水平井的摩阻系数反演[34]。
钻柱整体受力微分方程如下:
(1)
(2)
式中:F为钻柱上的轴向压力,N;q为单位长度钻柱重力,N/m;s为井深,m;“∓”代表钻柱状态,负号代表提拉,正号代表下入;α为井斜角,rad;EI为钻柱的抗弯刚度,N·m2;nt为钻柱与井壁的接触力,N/m;μ1为钻柱轴向摩阻系数;μ2为周向摩阻系数;M为钻柱所受扭矩,N·m;kb为井眼轴线曲率,m-1;Do为钻柱外径,m。
井眼曲率kb的计算公式为:
(3)
钻柱上的接触力nt计算公式为:
(4)
(5)
(6)
(7)
利用有限差分数值求解方法,可将钻柱整体受力模型写成以下形式(以压力为正):
(8)
(9)
式中:Fi和Fi+1分别为第i段钻柱单元靠近地面和靠近井下两端的轴向力,N;MTi和MTi+1分别为第i段钻柱单元两端的扭矩,N·m;Kbi和Kbi+1分别为第i段钻柱单元两端的井眼曲率,m-1;qi、EIi、nti、Δsi和Dbi分别为第i段钻柱单元的线重力(N/m)、抗弯刚度(N·m2)、与井壁的接触力(N/m)、单元端的长度(m)和外径(m)。
本文采用二分法进行钻柱周向摩阻系数的反演,计算流程如图2所示。根据计算流程,首先在0~1之间选取初始摩阻系数,将其与预测得到的井底钻压、扭矩带入到钻柱整体受力模型中[35];利用传递方程从钻头处向上依次计算每个钻柱单元的轴向力和扭矩;比较计算得到的地面大钩载荷、转盘扭矩与实际测量值的误差,若误差在允许范围内,则说明给定的摩阻系数较为准确,否则利用二分法对摩阻系数重新赋值计算。最终输出钻柱整体轴向力、扭矩及弯矩等参数。
图2 二分法反演摩阻系数计算流程图Fig.2 Flowchart for calculating friction coefficient by inversion using dichotomic method
本文采用岩屑运移2层稳态模型进行岩屑床高度的计算。模型将井筒内的岩屑床截面分为悬浮层和移动岩屑床2层运移,如图3所示。其控制方程包括:岩屑质量守恒方程、液相质量守恒方程、悬浮层动量守恒方程、移动岩屑床动量守恒方程、悬浮层岩屑扩散方程[36]。
图3 悬浮层和岩屑移动床2层运移示意图Fig.3 Schematic diagram of two-layer transport of suspended layer and cuttings bed
岩屑质量守恒方程:
UhChAh+UbCbAb=UsCsA
(10)
液相质量守恒方程:
Uh(1-Ch)Ah+Ub(1-Cb)Ab=Us(1-Cs)A
(11)
悬浮层动量守恒方程:
(12)
移动岩屑床动量守恒方程:
(13)
悬浮层岩屑扩散方程:
(14)
依据岩屑床高度计算环空几何结构参数,采用迭代法和割线法进行求解,求解流程如图4所示。
图4 2层岩屑床运移模型求解流程图Fig.4 Flowchart for solving of two-layer cuttings bed transport model
首先给定岩屑床高度(yb)和悬浮层质量分数(Ch)的初始值,根据质量守恒方程计算岩屑床速度(Ub)和悬浮层速度(Uh);根据扩散方程得到悬浮层岩屑质量分数。与初始值对比,若误差超出允许范围,则重新赋值悬浮层质量分数,迭代计算;根据2个动量守恒方程,计算单位长度压降,判断二者误差是否在允许范围内,若不在允许的范围内,则采用割线法重新假设岩屑床高度进行计算。经反复迭代计算后,最终得到岩屑床高度等参数。
提速降本是油气钻井的优先目标,但优化钻速等钻井效率指标时还需考虑到潜在的钻井风险。本文以提高机械钻速、降低机械比能为目标,以钻压、转速和排量为调控参数,以减小管柱摩阻、促进井筒清洁为约束条件,基于钻速智能预测模型、摩阻系数智能反演模型和水平井岩屑运移模型,构建了钻井提速-减阻-清屑多目标协同优化框架,如图5所示。多目标协同优化具体流程(步骤)如下。
图5 水平井钻井提速-减阻-清屑多目标协同优化框架Fig.5 Multi-objective collaborative optimization framework for ROP improvement, drag reduction and hole cleaning in horizontal well drilling
(1)处理数据,训练钻速智能预测模型。采用1.1节所述方法处理数据缺失值和异常值,基于1.2节中的随机森林模型和超参数调优方法训练、优化机械钻速智能预测模型。
(2)设定工程参数上、下限,得到候选参数组合。人工智能模型是从给定训练数据集上学习目标变量与输入变量的映射关系,所以当训练好的智能模型应用于非训练集数据时,其输入参数一般不能与训练集参数差别过大,否则将难以取得准确的预测结果。钻压、转速和排量作为司钻调控钻井过程的主要变量,若根据钻井设备上标定的范围取值,有可能取到不在训练集内的参数。因此,为保证预测模型的准确性,本文分别从训练数据集中取钻压、转速和排量的最大值和最小值作为其取值上、下限,在上下限之间随机取值,得到候选调控参数(钻压、转速和排量)组合。
(3)基于NSGA-Ⅱ算法进行多目标优化。本文以机械钻速和机械比能作为优化目标,属于多目标优化问题,文中采用改进非支配排序遗传算法(Non-dominated sorting genetic algorithm -II,NSGA-Ⅱ)[31,37-40]求解该问题。NSGA-Ⅱ是目前广泛使用的多目标遗传算法之一,它降低了非劣排序遗传算法的复杂性,具有运行速度快、收敛性好的特点。将步骤(2)中候选调控参数组合输入NSGA-Ⅱ算法中,求解后可获得帕累托解集以及其对应的参数组合,此解集和参数组合将作为后续进行风险约束和多方案决策的基础。
(4)基于摩阻系数反演模型和岩屑运移模型施加风险约束。步骤(3)计算帕累托解集时仅考虑了机械钻速和比能2个目标,其对应的调控参数(钻压、转速和排量)有可能造成管柱屈曲、卡钻以及井眼清洁等问题。因此,需要排除这些有可能造成钻井风险的解。将优化后的参数组合输入摩阻系数反演模型和岩屑运移模型,分别计算该参数组合下的摩阻系数、摩阻力和岩屑床高度,根据摩阻系数判断是否有卡钻风险,根据岩屑床高度判断是否符合井眼清洁条件。本文以周向摩阻系数大于0.5为判断有卡钻风险的依据[34],以岩屑床高度大于井眼直径10%为判断不满足井眼清洁条件的依据[41]。淘汰不满足风险约束的参数组合后,即得到可行的参数组合。
(5)基于TOPSIS算法进行多方案决策。步骤(4)淘汰了帕累托解集中有可能导致钻井风险的解,得到了所有可行的参数调控方案。本文采用逼近理想解排序法(Technique for Order Preference by Similarity to an Ideal Solution,TOPSIS)[31,42],从多个可行方案中决策给出最佳方案。该算法是一种广泛使用的综合评价方案,它能够充分利用原始数据信息,通过计算各方案与最优方案、最劣方案间的距离来进行评分,并以评分优劣为依据进行排序。基于TOPSIS算法对步骤(4)得到的可行参数集合进行决策,可获得最优参数组合方案。此方案能够取得机械钻速和机械比能的相对最优值,且不会引起屈曲、卡钻和岩屑床堆积等钻井风险。
将本文提出的水平井钻井提速-减阻-清屑多目标协同优化框架应用于该井7 600~7 650 m的水平段,分析其在机械钻速预测、钻速-比能协同优化、风险约束等方面的效果。
图6展示了协同优化方法在训练集(7 600 m前数据)上的机械钻速预测效果。该井训练数据集上的平均机械钻速为5.6 m/h,与本文提出的机械钻速智能预测模型平均绝对误差约为0.2 m/h,平均预测误差约为3.9%,且95%情况下其绝对和相对预测误差分别小于0.7 m/h和10%。
图6 部分训练集(7 600 m前数据)机械钻速预测效果Fig.6 Prediction effect of ROP of some training sets (data before 7 600 m)
图7对比了7 600~7 650 m的实际机械钻速、预测机械钻速和优化后机械钻速。钻进7 600~7 650 m时的平均机械钻速为2.7 m/h,该模型预测的平均机械钻速为2.5 m/h;其平均绝对预测误差为0.3 m/h,平均相对预测误差为11.6%。优化后,平均机械钻速提高至3.4 m/h,相较优化前提高了32%,且仅在7 638 m、7 639 m和7 640 m处由于钻速预测模型误差过大(约15%)而低于优化前。
图7 7 600~7 650 m水平段实际钻速、预测钻速和优化后钻速对比Fig.7 Comparison of actual,predicted and optimized ROPs in 7 600 to 7 650 m horizontal section
图8对比了7 600~7 650 m的真实机械比能和优化后机械比能。优化后机械比能相比优化前降低了约450 MPa,降幅达19%,且仅在7 639、7 640 和7 641 m处未达成优化效果,其原因与钻速优化情况类似。因此,本文提出的协同优化方法在钻速预测和钻速-比能多目标优化方面都取得了较好的效果。
图8 7 600~7 650 m水平段实际比能和优化后比能对比Fig.8 Comparison of actual and optimized specific energies in 7 600 to 7 650 m horizontal section
在风险约束方面,图9与图10为优化前后平均摩阻系数和平均岩屑床高度的对比。
图9 7 600~7 650 m水平段优化前、后平均摩阻系数Fig.9 Average friction coefficient before and after optimization of 7 600 to 7 650 m horizontal section
图10 7 600~7 650 m水平段优化前、后平均岩屑床高度Fig.10 Average cuttings bed height before and after optimization of 7 600 to 7 650 m horizontal section
由图9可知,优化前平均摩阻系数约为0.32,优化后为0.24,降低了25%。由图10可知,优化前平均岩屑床高度为约1.47 cm,优化后约1.36 cm,降低了7.5%。表3对比了7 649 m处优化前、后的工程参数、效率指标和风险指标。相比优化前,钻压降低了11 kN,转速基本保持不变,排量增大了14 L/s。可见,优化后机械钻速和机械比能有少许提升,且由于钻压降低和排量提升,优化后相关风险指标有较大改善,表现在:优化后平均摩阻系数降低0.12,平均摩阻力减小38 kN,平均岩屑床高度降低0.15 cm。由此可见,在保证钻进效率的同时减小了卡钻、井眼清洁不足等风险。
表3 某井7 649 m处优化前、后工程参数和相关指标对比Table 3 Engineering parameters and related indicators before and after optimization at 7 649 m of a well
本文综合考虑水平井钻井破岩提速、管柱减阻和井筒清洁3个子系统协同作用,提出了钻井过程协同优化方法。基于NSGA-Ⅱ和TOPSIS算法进行多目标优化决策,给出钻压、转速和排量3个工程参数的推荐值。最后,以某水平井7 600~7 650 m水平段为例分析了其在机械钻速预测、机械钻速比能协同优化、风险约束等方面的效果。结果表明,该协同优化方法能够提高机械钻速、降低机械比能,且保证摩阻系数和岩屑床高度在安全范围内,实现安全高效钻进。
综上所述,本文提出的水平井钻井提速-减阻-清屑多目标协同优化方法能够在优化钻井效率的同时降低卡钻和井眼清洁不足等风险,可为现场安全高效钻进提供保障。建议后续可在以下方面改进:①采用3层模型,细化描述岩屑运移过程中流型变化的影响;②探索更高效的多目标优化与决策策略;③考虑井眼轨迹等其他钻井子系统。