钻井模型与机器学习耦合的实时卡钻预警技术*

李紫璇 张菲菲 祝钰明 吴墨染 刘泳敬 于 琛

(1.长江大学石油工程学院 2.油气钻井技术国家工程实验室防漏堵漏研究室 3.油气钻采工程湖北省重点实验室4.青海油田分公司 5.渤海钻探第一钻井分公司 6.渤海钻探第三钻井分公司 7.渤海钻探工程技术研究院)

0 引言

统计数据表明[1]，钻进过程中大约有70%的时间损失与钻具阻卡有关。如不能及时发现，一旦发生卡钻事故，往往需要停钻处理，将会花费大量的时间及处理成本，并有可能引起更为严重的次生危害[2]。目前解决此类问题最有效的方法是对钻井过程中的卡钻风险进行准确预测分析，并对风险的处理提供合理的决策方案以规避卡钻，缩短非作业时间[3]。

目前，卡钻预警方法多以基于大数据的统计分析预警方法和基于钻井模型的分析模型为主。前者主要分为判别分析法和模式识别法。该方法以收集大量有效数据为前提，通过对卡钻发生原因、卡钻时的作业状况等进行分析，利用神经网络、多元统计分析、模糊逻辑及层次分析法等方法建立模型或专家系统判断是否发生卡钻[4－7]。要想实现实时卡钻预测，还需建立时间序列ARMA 模型[8]。用于统计分析方法的数据多来自于地层信息和邻井的卡钻事故，若缺乏邻井资料，则难以建立有效的卡钻检测系统和方法。同时，收集到的数据往往受地域限制，使得利用统计分析的卡钻预警模型难以进行大范围推广[9]。卡钻预警方法使用的物理模型多为摩阻扭矩模型，使用摩阻扭矩模型对摩阻系数进行反算，将摩阻系数随深度的异常变化作为卡钻风险监测的指标[10]；分析对比扭矩预测值和扭矩实测值，对钻井过程中的异常进行监测[11－12]。卡钻发生时会导致钻井参数出现异常变化，但传统物理模型所需计算量大且难以完全反映卡钻发生时的变化规律[13]，在实际应用中误报率高、准确率低、实时性较差，且难以对发生的卡钻进行分类[14]。

本文将瞬态岩屑运移模型、摩阻扭矩模型和机器学习相结合，研究基于实时录井数据的卡钻风险监测及预警技术，利用瞬态岩屑运移模型对岩屑床分布和井眼清洁状况进行实时监测；将瞬态岩屑运移模型的计算结果作为输入，建立考虑井筒岩屑分布对钻杆受力影响的改进摩阻扭矩模型；利用贝叶斯优化训练钻井模型，使模型能够不断进行自我调整以适应新的工况，即使改变井况或井型，也能不断提高模型精度，且模型训练速度快，仅需当前井的录井数据便能完成模型训练。通过分析对比混合模型实时预测结果和实测数据，并进行时序数据分析，最终得到卡钻风险指数。

1 钻井模型

1.1 实时井眼清洁监测模型

井眼清洁不充分是造成卡钻的常见原因之一，约占比例的50%[3]。在钻进过程中对整个井眼的岩屑分布情况进行准确且实时地监测是有效预测卡钻事故的关键。

采用瞬态岩屑运移模型对岩屑沿井眼轨迹分布进行一维数值模拟[15]，将岩屑运移过程在时间域上进行离散化，通过实时计算每个离散单元的岩屑堆积量和环流压力梯度实现井眼清洁实时监测。模型公式基于质量和动量守恒方程，公式为:

式中:角标α和β分别为流型中的不同部分(也可称为层)；n为给定流型的总层数；A为层的横截面积，mm2；C为岩屑密度，g/cm3；U 为流速，m/s；S为层的周长，mm；p 为流动压力，Pa；ρ为层的平均密度，g/cm3；τ 为切应力，Pa；Ex为不同层之间的岩屑体积交换率；g 为重力加速度，m/s2。

在数值模拟中使用基于有限体积的半隐式差分方法来求解方程。在钻井过程中，钻井参数不断变化，岩屑的流动形态也会随之改变，为使模型能同时处理时间域和空间域中岩屑流型的变化，使用基于瞬态力学模型的广义流型[16]，在每次计算前对离散单元的岩屑流型进行判断，根据岩屑流型建立相应的模型方程组。所有流型的每个层都有自己的质量守恒和动量守恒方程，给定流型的方程总数将取决于其总层数。模型计算过程中，首先根据初始状态通过动量守恒方程求解下一个时间步的速度和压力，将收敛的速度和压力带入质量守恒方程求得新的悬浮岩屑密度和沉积岩屑体积，最后把结果带入下一个时间步，依次进行计算。模型建立和求解的详细过程见文献[17]。

1.2 实时摩阻扭矩计算模型

随着井筒中岩屑床的堆积，井筒内井眼清洁状况变差，钻杆扭矩和起下钻时的大钩载荷增加明显，卡钻风险也随之增加[18－19]。若要模拟岩屑床对钻杆运动的精确影响，需要对整个井眼中钻杆在环空中的位置进行复杂的有限元模拟。为便于实时计算，在实际应用中，通过使用一个附加力来表示岩屑床和钻杆之间的复杂相互作用:起下钻时，附加力施加在钻杆轴向上，表示岩屑床对钻杆轴向作用力的影响；钻杆旋转时，附加力与钻杆壁曲面相切，方向与旋转方向相反。

为便于描述空间井眼轨迹，将空间直角坐标系与弗莱纳坐标系相结合，如图1 所示。通过对钻杆单元体进行受力分析，建立旋转钻杆的力和力矩平衡方程，方程如下[20]:

图1 钻杆单元体在3D 井眼中的受力图Fig.1 Force diagram of drill pipe unit body in 3D wellbore

式中:Ft为轴向力，N；wp为单位钻杆在钻井液中重力，N；wc为接触力，N；tz、nz和bz分别为z方向上单位切向、法向和副法向向量；θ为沿平面法线与接触力间的夹角，(°) ； Fc为岩屑床施加的附加力，N；φ为井斜角，(°) ；μ为摩阻系数，无量纲；rp为钻杆半径，m；Mt为钻杆旋转所需的轴向扭矩，N·m；ds为钻杆单元的长度，m；κ为曲率，m－1。

通过上述方程，可得Fc等于钻杆上的附加力。

除了每个单元体的附加摩擦力需要基于局部岩屑床的高度来计算之外，用于解改进摩阻扭矩模型的算法与解标准摩阻扭矩模型的算法相似。

实时井眼清洁监测模型与摩阻扭矩模型结合应用时，以实时操作参数作为实时井眼清洁监测模型的输入，预测岩屑床在井眼中的实时动态分布。摩阻扭矩模型将岩屑床的分布预测结果作为输入，计算出每个离散点的附加力，最后得出钻杆上的摩阻扭矩值。为进一步研究岩屑床对钻杆的影响，该研究利用试验设备[21]，通过一系列试验数据，探究附加力与岩屑床高的关系，试验参数如表1 所示。

表1 试验测试参数Table 1 Experimental test matrix

为了便于研究相关性，分别对岩屑床高和岩屑床高对扭矩的影响进行标准化处理，如下所示:

式中:hc为岩屑床高度，mm；D为井眼直径，mm；hs为标准化岩屑床高度，无量纲；Ts为标准化岩屑床高hs下的扭矩测量值，kN·m；Tc为相同井眼结构下清洁井眼的扭矩计算值，kN·m；TR为由岩屑引起的扭矩变化率，无量纲。

经标准化处理后，测试数据无量纲岩屑床高和扭矩变化率的关系如图2 所示。

由图2 可知，标准化岩屑床高和扭矩变化率的关系大致呈指数形式。因此，在数据集的基础上发展指数函数相关性，公式如下:

式中:a和b均为无量纲参数，图2 所示曲线中a和b分别为10.2 和2.025。

图2 标准化岩屑床高和扭矩变化率的关系图Fig.2 Relationship between standardized cuttings bed height and torque change rate

由于试验是在模拟直井的设施上进行的，扭矩与井眼和钻杆间的接触力、井眼和钻杆间的摩阻系数成正比。因接触力在试验过程中保持不变，故扭矩的变化只与摩阻系数的变化有关。又因参数是无量纲的，故公式(10) 可转换为表示标准岩屑床高与摩阻系数变化率之间的关系，如公式(11)所示。

式中:fR为摩阻系数的变化率，可以用公式(12)定义。

式中:fs为标准化岩屑床高下的摩阻系数；fc为没有岩屑时的摩阻系数。

将公式(11) 应用于前文所述的摩阻扭矩计算模型中，用于计算有岩屑井的扭矩值。与试验测试不同，在实际钻井应用中需要对参数进行改变。虽然由试验数据得出的相关性可能不适用于所有情况，但是摩阻系数和无量纲岩屑床高的变化趋势应是相同的(即a和b的值可能在不同的情况下有所不同，但公式(11) 仍然有效)。

2 机器学习和模型训练

在摩阻扭矩模型计算过程中，钻杆被分成若干离散单元，由于岩屑床高沿钻杆发生变化，所以每个单元体都需要一个特定的摩阻系数。前文所述摩阻系数与岩屑床高呈指数关系，但a和b的值可能在不同的工况下有所不同，为使模型应用范围更广，更好地适用于现场实际，需要以历史钻井数据为输入，使用参数估值的方法对模型进行训练。

通过现场数据训练改进摩阻扭矩模型的过程，本质上是对考虑井眼清洁影响的摩阻扭矩模型进行参数估值，可以视为寻找最优参数集(即摩阻系数沿井眼的分布情况)，使得井口扭矩与转盘扭矩差值最小，是一个非线性优化问题。贝叶斯优化算法可用较少的评估次数求得复杂目标函数的最优解，是较为常用的优化算法，其关键在于概率代理模型和采集函数的建立。

训练模型优化过程需要最小化扭矩差值，选用纳什效率系数(NSE) 作为模型的评价指标，表示模型拟合方差占总方差的百分比[22]。因极值对NSE影响较大，本文选用lgNSE作为模型评价指标，可有效降低极值对模型总体拟合效果的影响，如公式(13) 所示[23]。式中:yt和分别代表时间步为t时的扭矩实测值和计算值，代表实测值的平均值。

计算结果越接近于1，参数优化效果越好。贝叶斯优化的目的是使lgNSE最大化，以找出最有可能的摩阻系数分布情况。

使用概率代理模型替代评估代价高昂的复杂目标函数，根据贝叶斯定理不断地进行迭代以增加信息量，对先验知识进行修正[24]。高斯过程函数具有高度的灵活性、可扩展性和可分析性，对线性和非线性关系均具有良好的替代性[25]，因而成为贝叶斯优化中应用最广泛的概率代理模型。高斯过程是多元高斯概率分布的范化，由均值函数m(x) 和半正定的协方差函数k(x，x′) 构成，如公式(14)所示。

式中:g(x) 为多元高斯概率分布，GP为高斯过程。

在使用过程中，为方便计算，假定先验均值函数为0，该设定对后验分布的准确性几乎没有影响[26]，协方差函数采用的是一类高度灵活的Matérn 协方差函数簇，可产生二阶可微的样本函数，如公式(15) 所示。

式中:v为平滑参数，l为尺度参数，Kv为第二类变形贝塞尔函数。

采集函数根据概率代理模型的后验结果构造主动选择策略，被用于确定贝叶斯优化的方向，选择采集函数最大值点作为下一个评估点。采集函数有多种类型，如提升概率(PI)、期望提升量(EI)、置信上边界(UCB) 等，其中EI能将提升概率和提升量整合起来，所需参数少且能兼顾全局搜索和局部变化，如公式(16) 所示。

式中:y为当前最优函数值，φ为标准正态分布累计密度函数，μt(x) 为均值，σt(x) 为方差，ξ为平衡参数，无因次。

本文使用贝叶斯优化对改进摩阻扭矩模型进行训练的流程图如图3 所示。图3 中录井数据为实时录井数据，井眼数据为井身结构、井眼信息、钻具组合和钻井液参数等静态数据。

图3 使用贝叶斯优化训练摩阻扭矩模型流程图Fig.3 Flow chart of torque and drag model training by Bayesian optimization

3 实时卡钻风险系数分析

本文所用实时卡钻风险预测方法通过分析3 个不同来源的扭矩数据:考虑岩屑影响的扭矩预测值、不考虑岩屑影响的扭矩预测值以及离底自由旋转扭矩的测量值；生成两个独立的卡钻风险指数，即与井眼清洁相关的卡钻风险指数和与其他因素相关的卡钻风险指数。由于摩阻扭矩模型已通过数据驱动模型以历史数据为基础做好了调试，所以实测值与预测值出现偏差表示井下异常情况的出现。本文选用时序数据分析的方法对实时数据的异常趋势进行检测。首先计算相对偏差值，计算公式如下:

式中:Tc为考虑岩屑影响的扭矩预测值，kN·m；Tn为不考虑岩屑影响的扭矩预测值，kN·m；Tr为离底自由旋转扭矩的测量值，kN·m；Rc为考虑岩屑影响与不考虑岩屑影响的扭矩预测值之间的相对偏差，无量纲；Rm为实测值与考虑岩屑影响的扭矩预测值之间的相对偏差，无量纲。

在趋势分析中，移动平均值偏差方法被用于追踪两个相对偏差值随时间的变化趋势，移动平均偏差值为正表示数据处于增长趋势，为负表示数据处于降低趋势[15]，如公式(19)所示。

最后，计算出表示卡钻事故发生概率的风险指数，为实时监测提供更直观的监测结果和警报。最终卡钻风险指数将被转换成介于0 到1 之间的值，0 表示没有卡钻风险，1 表示钻杆有极大可能会卡死，公式为:

4 案例分析

选用钻井实时数据对模型进行验证。图4 为一段时长为290 min 的卡钻发生前的实际井录井数据，该数据每10 s 记录一次。卡钻发生在钻进过程中，井斜角约为56°。如图4 所示，在早上5:20 左右，当其他参数基本保持不变时，扭矩从15 kN·m 增加到18 kN·m，且增长速度较快。

图4 卡钻发生之前的实时钻井数据Fig.4 Real-time drilling data before pipe sticking

为更好地解释扭矩的变化，将扭矩与井深相对应，如图5 所示。从2 255 m 钻到2 450 m 期间扭矩测量值持续增加。考虑岩屑影响的扭矩预测值与不考虑岩屑影响的扭矩预测值之间的偏差持续增加，表明井眼中的岩屑堆积在阻碍钻进。在钻到大约2 410 m (早上5:20) 时，扭矩从15 kN·m 迅速增加到18 kN·m，井下出现异常状况。由于钻井队没有注意到扭矩的迅速增加，继续向下钻进，在钻杆卡死前又持续钻了40 m。在后续钻的这40 m 中，扭矩依然在增加。

图5 扭矩随深度变化图Fig.5 Change of torque with depth

将实时卡钻风险分析应用于钻后分析和算法验证。卡钻事故发生前5 h 的风险指数如图6 所示。

图6 实时卡钻风险分析结果Fig.6 Analysis results of real-time pipe sticking risks

由图6 可知，与井眼清洁相关的卡钻风险指数在2:05 左右已经达到0.6，并在事发前上升到0.7 以上。5:50 前，与其他因素相关的卡钻风险指数处于极低水平(0.2～0.4 之间)，表明目前引起卡钻的主要因素与井眼清洁相关。当扭矩突然增大，与其他因素相关的卡钻风险系数突然增加到0.8 以上，表明发生卡钻事故的可能性很大。通过实例分析，该模型能够捕捉到卡钻是否发生，同时能够对要发生的卡钻进行分类，并在事故发生前给出预警信号。因此，如果将该工具应用于钻井监测系统中，将有助于避免卡钻事故的发生。

5 结论

(1) 研究了一种实时卡钻预警技术，该技术将钻井模型和机器学习相结合，通过分析实时钻井数据对卡钻进行实时监测，并区分即将发生的卡钻是否与井眼清洁有关。

(2) 将摩阻扭矩模型和瞬态岩屑运移模型相耦合，充分考虑了岩屑床对钻杆的影响。

(3) 使用贝叶斯优化训练钻井模型，使钻井模型能够适应操作参数的变化不断进行自我调整的同时，也能不断提高模型预测精度。

(4) 利用时序数据分析监测实时数据的异常趋势，通过计算卡钻风险指数，为卡钻事件提供直接的实时监控和警报，可帮助工程技术人员及时采取有效措施避免卡钻事故发生。