刘慕臣，宋先知，李大钰，朱 硕，付 利，祝兆鹏，张诚恺，潘 涛

(1.中国石油大学 (北京)人工智能学院，北京 102249；2.中国石油大学 (北京)石油工程学院，北京 102249；3.油气资源与工程全国重点实验室，北京 102249；4.中石油江汉机械研究所有限公司，湖北 武汉 430024；5.中国石油集团工程技术研究院有限公司，北京 102206；6.中国石油大学 (北京)机械与储运工程学院，北京 102249)

摩阻扭矩是钻井管柱受力的重要分析参数，可用于大钩载荷预测、管柱可下入性分析、井眼轨道优化设计与卡钻风险预测等。在钻井过程中钻柱与井壁直接接触产生巨大的接触力和摩擦力，使钻柱受到巨大磨损，为钻进过程埋下潜在安全隐患[1]。此外，钻柱摩阻过大将导致钻压和扭矩不足，降低钻进效率，钻柱与井壁间的巨大摩擦力会严重破坏井壁稳定性，诱发井塌等井下事故，对钻井效率与安全影响巨大[2]。因此，准确预测钻柱摩阻扭矩对油气钻井的提质增效具有重要意义。

对于摩阻扭矩计算方法，目前国内外已有许多学者进行了研究，包括机理模型方法和智能化方法两个方面。在机理模型方面，最经典且应用最广泛的是C.A.Johancsik 等[3]提出的软杆模型，该模型忽略了管柱的抗弯刚度，为摩阻扭矩分析奠定了理论基础。M.Lesage 等[4]在Johancsik 模型的基础上分析了钻井工况的影响，通过软杆模型反演起下钻、旋转摩擦因数。J.F.Brett 等[5]将软杆模型进行实际应用，通过反演摩阻系数分析井眼状况。黄文君等[6]考虑了钻柱接头的影响，建立了修正的摩阻扭矩模型。钻柱摩阻扭矩机理模型发展愈发成熟，摩阻系数作为其中的关键参数，反映了阻碍钻柱与井壁之间相互运动的量化程度[7]，是对摩阻扭矩预测的关键依据，但受到岩屑床高度、扶正器等诸多因素影响，尚未形成直接计算钻柱摩阻系数的有效方法，目前工程应用以经验估算或钻后反演为主，准确性和时效性有待提升[8-9]。在智能化方法方面，国内外相关研究刚刚起步。朱硕[2,10]等通过实时智能计算井底钻压扭矩反演摩阻系数，进行摩阻扭矩实时分析而非超前预测，且钻压预测准确性直接影响摩阻系数反演准确性。C.Hegde 等[11]通过地面实测工程数据及测井资料，结合机器学习方法(MLR、SVM、RF)预测摩阻扭矩。丁培宇[12]、申静波[13]等基于实钻数据使用BP 网络预测摩阻系数，但是未考虑摩阻系数随井深分布的序列特征、未分析黑箱模型可解释性。祝兆鹏等[14]通过提取特征及时序的权重系数表征参数重要性，实现了井底压力预测结果的可解释。

人工智能方法在复杂的非线性映射问题求解方面具备显著优势，引领了新一代变革性钻井技术并形成了油气钻完井人工智能应用场景体系[15-16]。笔者提出了一种机理-数据融合的钻柱摩阻扭矩预测方法，基于已钻井数据使用钻柱摩阻扭矩模型反演摩阻系数，为智能预测模型提供数据基础，进而使用考虑数据序列特征的LSTM 网络建立摩阻系数智能预测模型，结合实时数据和黑箱模型解释方法对预测模型进行合理性验证，形成机理数据融合的钻柱摩阻扭矩智能预测方法，为钻井管柱力学的准确表征与动态分析提供理论与技术支撑。

1 智能预测方法流程

摩阻扭矩智能预测整体工作流程如图1 所示：综合录井数据处理，包括时间域、深度域数据；针对起下钻工况使用软杆模型反演摩阻系数，进行特征工程确定摩阻系数数据集；基于LSTM 的钻柱摩阻系数智能预测；基于SHAP 方法的摩阻系数智能预测模型解释方法；钻柱摩阻扭矩预测及验证。

图1 智能预测方法工作流程Fig.1 Workflow of intelligent prediction method

2 综合录井数据处理

首先对时间域和深度域数据进行收集与处理，然后反演钻柱摩阻系数，将反演值作为样本标签结合相关参数建立摩阻系数深度域数据集。

2.1 时间域数据处理

数据来源为国内某油田，共计74 口井，时间域数据库数据为每2 s 采集1 个点，贯穿钻井全流程，包含起钻、下钻、钻进、短起下、下套管、固井、电测以及复杂事故处理。对于钻柱摩阻扭矩模型，管柱井口与井底的边界条件是摩阻系数反演的必要条件，井口钩载有实测值，为了避免钻压预测不准导致摩阻系数反演不准，本研究主要针对钻压为0 的起下钻工况进行摩阻系数反演,用以预测分析钻进工况。首先需要从时间域数据库提取起下钻工况数据。根据钻头深度与井深大小关系及运动方向初步识别管柱运动工况，某井识别效果如图2 所示。为了摩阻系数反演的连续性和完整性，选取钻头从井底上提至井口或从井口下放至井底的完整过程，图3 显示了该井2019 年2 月15日和16 日2 d 数据。提取该时间段内下放工况数据，将钻头深度和大钩载荷的时间域数据进行深度域转化，具体方法包括对齐、删除、拼接等。转换前后结果如图4 所示，最终得到不同井深处摩阻系数反演的边界条件，为摩阻系数反演提供数据基础。

图2 工况识别效果Fig.2 Condition identification result

图3 起下钻工况数据Fig.3 Tripping condition data

图4 大钩载荷深度域转换Fig.4 Depth-domain coversion of hook load

2.2 深度域数据处理

深度域数据库中包括5 种数据表，各数据表的采样频率与数据类型见表1。

表1 各数据表采样频率与数据类型Table 1 Sampling frequency and data type of various data sheets

首先需要将5 种数据的采样频率统一为1 次/ m。对于井斜数据，通过计算各点井眼曲率将轨迹划分为直线和圆弧部分。直线部分每一米处的井斜角、方位角均一致。而圆弧段利用三次样条曲线进行插值，最终得到采样频率1 次/ m 的轨迹数据。泥浆性能数据主要是由现场工作人员定期记录，采样频率均为几十米，由于钻井液数据不会存在较大的变化，认为在某段内每米处的钻井液数据均与该段内测量点相同。类似地，钻头记录数据和岩屑描述记录数据也作如此处理。

统一数据采样频率后，需要对已有数据进行分析与处理，为下一步特征分析与智能建模提供基础。井斜数据与工程参数数据类型均为数值型数据。泥浆性能参数主要包括：漏斗黏度、塑性黏度、屈服度、钾离子含量、3 转读数、6 转读数、100 转读数、钻井液体系等。钻井液体系为字符型，其余均为数值型。通过数据分析发现钻井液体系细分种类较多，因此，粗分类为水基、油基、气基3 类大体系，其中水基包含：不分散、分散、聚合物、低固相、钙处理、饱和盐水钻井液等，气基包括空气与泡沫钻井液，见表2。

表2 钻井液体系数据处理Table 2 Data processing of drilling fluid system

钻头记录数据：主要包括钻头直径、钻头型号、入井新度、出井新度、钻压、转速、钻头水眼、钻具组合信息等。岩屑描述记录数据：原始岩屑描述记录中包含层位、井段、岩性定名、岩性及含油气水描述等信息，均是文本信息，需要利用python 中的正则表达式提取信息。

2.3 数据融合与清洗

2.3.1 数据融合与编码

经过统一采样频率、数据预处理后，以井深为索引将以上5 种数据进行数据拼接处理，形成综合录井数据集。由于神经网络不能直接识别字符型数据，而岩性等对摩阻系数的影响不能忽略，因此需要对字符进行编码。独热编码采用稀疏矩阵方式存储数据，需要占用更多神经网络参数空间，不利于模型稳定性，因此采用索引编码方法。以某一口井为例，编码结果见表3。

2.3.2 数据归一化

原始数据各个特征的尺度差异较大，易引起梯度问题，进行归一化或z-score 标准化，可以消除量纲的影响，避免梯度问题。数据集是以井深为索引的序列性数据，并无明显的正态分布性，因此，更适合归一化处理。

2.3.3 数据划分与模型评价

按照8∶2 比例划分数据，选取均方根误差ERMS(Root Mean Squared Error，RMSE)、平均相对误差δ、最大相对误差δmax、训练用时作为模型对大钩载荷预测效果的评价指标，其中训练用时可以间接反映模型的空间复杂度。

式中：yit为第i条数据的目标真实值；yip为第i条数据的目标预测值；N为样本总数量；δi为第i个值的相对误差。

3 钻柱摩阻系数数据集建立

3.1 钻柱摩阻系数反演

摩阻系数是软/刚杆模型计算摩阻扭矩的重要参数，智能模型准确预测摩阻系数的前提是准确的样本标签值。为了避免钻进工况下井底钻压预估不准确导致摩阻系数反演精度低、波动大，针对钻压为0 的起下钻工况反演摩阻系数。

3.1.1 钻柱摩阻扭矩模型

首先建立钻柱摩阻扭矩模型。摩阻扭矩模型已经有了大量的研究并逐渐成熟。目前根据是否考虑钻柱的弯矩，分为软杆模型和刚杆模型。采用经典的软杆模型进行摩阻系数反演。软杆模型差分方程如下：

式中：Fi、Fi+1分别为第i段钻柱微元段上下处的轴向力，N；q为第i段钻柱线重，N/m；为第i段钻柱平均井斜角，(°)；μ1与μ2分别为第i段钻柱轴向和周向摩阻系数；ni为第i段钻柱与井壁的单位长度接触力，N/m；Δsi为第i段钻柱长度，m；Mi、Mi+1分别为第i段钻柱上下处的扭矩，N·m；Db为第i段钻柱外径，m；αi、αi+1为第i段钻柱上下处井斜角，(°)；φi、φi+1为第i段钻柱上下处方位角，(°)。

摩阻扭矩求解流程如图5 所示。

图5 软杆模型求解流程Fig.5 Solution process of soft string model

3.1.2 钻柱摩阻系数反演方法

使用二分法反演全井段摩阻系数。起下钻时，钻头轴向受力为0。对钻柱微元进行受力分析，如图6所示。

图6 钻柱微元受力分析Fig.6 Force analysis of drill string element

可得钻柱微元上端轴向力Fi的计算公式(5)。在起下钻工况下认为钻压为0，因此，将该式从井底累积求和至井口可得井口轴向力计算如下式。

式中：|F|为井口轴向力大小，N；Wi为第i段钻柱微元沿切线方向的重力大小，N；μi为第i段钻柱微元的轴向摩阻系数；nti为第i段钻柱微元与井壁的接触力大小，N；j为单元总数目；“ ∓”中，负号代表上提，正号代表下入；每段的Wi已知，nti可以通过式(4)求得，大钩载荷现场可以实际测量。

式(9)中只有摩阻系数μi是未知的。由于总共有j个μi，无法一次反演求解，需要从井口不断重复反演计算至井底，才能得到j个μi。反演方法采用二分法。以钻头下入工况第i段钻柱微元为例，设定初始摩阻系数取值范围为(0，1)，计算步骤为：第一步，令μ1=0，μ2=1，精度阈值ε=0.01。第二步，令μa=(μ1+μ2)/2。第三步，将μa代入式(9)，计算井口钻柱轴向力Ta。以管柱下入为例，若地面钩载实测值TR＞Ta则表明真实摩阻系数μR＜μa，设μ2=μa；若Ta＞TR则表明真实摩阻系数μR＞μa，则μ1=μa。第四步，若 |1-Ta/TR|＞ε，则返回第二步；否则停止程序，μa为最终摩阻系数反演值。继续令式(9)中μi=μa，然后继续反演第i+1 段钻柱微元的摩阻系数μi+1。摩阻系数反演流程框图如图7 所示，某两口井的摩阻系数反演结果如图8 所示。

图7 钻柱下入摩阻系数反演流程Fig.7 Inversion process of friction coefficient for tripping of drill string

图8 钻柱摩阻系数反演值Fig.8 Inversed value of friction coefficient of drill string

3.2 相关性分析与特征优选

预测摩阻系数首先需要确定与摩阻系数有潜在映射关系的参数。由于摩阻系数是钻柱与井壁间的固有属性，在选取参数时，应更多考虑井筒内的固有属性参数。首先是轨迹数据，选取井深、井斜角、方位角、井眼曲率、闭合距、闭合方位共6 个参数。钻井液性能参数对摩阻系数同样有重要影响，尤其是钻井液体系，选取钻井液密度、黏度、屈服度、塑性黏度、600 转读数与钻井液体系共6 个参数。井身结构、钻具参数与地质数据同样有不可忽略的影响，选取钻柱内径、钻柱外径、井径与岩性共4 个参数。基于16 个参数结合摩阻系数建立钻柱摩阻系数数据集，为智能预测提供样本支撑。

3.2.1 特征定量化分析与优选

首先进行降维处理，提取强相关特征，起到简化模型，降低过拟合，增强泛化能力的作用。采用Filter 方法[17]中的距离相关系数(Distance correlation coefficient)，定义如下。

式中：Dcorr(X,Y)为X,Y的距离相关系数；Dcov(X,Y)为X,Y的距离协方差；Dcov(X,X)为X的距离方差；Dcov(Y,Y)为Y的距离方差。

摩阻系数是按照井深索引排列的，且某一井深处的摩阻系数与上部地层的摩阻系数存在某种复杂的影响关系，包括地层数据的连续性等，属于序列数据范畴。因此，需要利用相关性分析探究序列性质以便指导模型选择。图9 和图10 为18 个参数的相关性系数。μ1m和μ10m分别表示上1 m 和上10 m 处的摩阻系数，由图发现，二者与μ相关系数较大，说明数据内部含有较强的序列性，因此，相比于BP 网络，考虑序列性质的LSTM 更适用于本研究。为了减少模型参数量，增强模型鲁棒性，选取与摩阻系数相关性强，而与其他输入参数相关性弱的参数。钻柱内径、钻柱外径与摩阻系数的相关系数均为0.2，但二者彼此相关系数高达1.0，故仅保留钻柱外径。最终选取15 种参数，包括标准井深、井斜角、方位角、井眼曲率、闭合距、闭合方位、钻井液密度、黏度、屈服度、塑性黏度、600 转读数、钻井液体系、井径、钻柱外径、岩性。

图9 各参数与μ 的距离相关系数Fig.9 Distance correlation coefficient between parameters and μ

图10 各参数间的距离相关系数热力图Fig.10 Thermodynamic diagram of the distance correlation coefficient between parameters

4 钻柱摩阻系数智能预测方法

4.1 基于LSTM 的摩阻系数智能预测

LSTM(Long Short-Term Memory)擅长处理和预测序列数据[18]。神经元包含输入门i、遗忘门f、输出门o、记忆单元C。“门”结构使用sigmoid 作为激活函数，输出一个0 到1 之间的数值，描述通过此结构的信息量，进而影响神经网络状态，避免产生梯度爆炸或梯度消失问题。十分适合处理钻井过程中的动态变化问题。本文利用LSTM 的非线性拟合和长期、短期记忆的能力表征15 种参数与摩阻系数的映射关系，同时考虑工程参数自身在深度维度上的变化关系，结构如图11 所示。

图11 LSTM 模型结构Fig.11 Structure of LSTM model

4.2 基于正交实验的超参数优选

正交实验法是研究多影响因素的实验设计方法，它依据伽罗瓦理论选择部分具有代表性的参数组合进行实验，可以有效减少实验次数。对于LSTM 网络，设计3 种不同的神经网络结构，每种网络结构设计4种超参数组合，包括神经元个数、激活函数、学习率、Dropout，共计12 种网络模型。为了测试鲁棒性，训练过程不设置固定的随机种子参数，以保证充分的随机性。每组进行10 次训练与测试，记录实验结果的平均值与标准差。每次训练的Epoch 为3 000。结果见表4。

综合考虑均方根误差、平均相对误差、最大相对误差以及一口井数据的训练用时共4 个指标对MLP模型进行优选。由表4 和图12 可知，均方根误差和平均相对误差最小为0.000 12 和3.1%，均对应第4 个模型，且最大相对误差为26.4%也处于相对低水平。鲁棒性方面，3 个评价指标的标准差区间均最小，意味着模型在初始参数完全随机的情况下，性能表现较为稳定。但训练用时35 s 相比较长，表明模型复杂度更高。相比于准确度要求，本预测任务对模型反应速度无较高要求，只需在同一数量级内且无指数级增长趋势即可，因此，优选第4 个模型。

图12 12 个模型预测评价指标对比Fig.12 Comparison of evaluation indexes for prediction of 12 models

5 基于SHAP 的摩阻系数智能预测模型解释方法

LSTM 作为一种神经网络是典型的不可解释黑箱模型，其性能能够达到较高水平但无法得知内部决策机理以及输入对输出的影响等，导致模型的可靠性不足。因此，有必要对建立的摩阻系数智能预测模型使用黑箱模型解释方法，进一步分析模型的可靠性。机器学习可解释性方法主要分为全局解释方法和局部解释方法[19]。全局解释方法指的是全面理解黑箱模型内部复杂的运算逻辑与成因，受限于应用场景、具体算法等，当前阶段很难实现全局解释。局部解释方法更专注于分析模型输入与输出之间的逻辑关系，解释输入是如何影响输出变化的，而不是解释决策函数本身。由于局部解释方法发展较成熟且应用更广泛，本文使用该方法进行可解释性分析[20]。

SHAP(SHapley Additive explanation)是一种受博弈论启发的方法，通过计算每个特征对单一预测结果的重要性值分析可解释性。它是一种模型无关的方法，适用于任何单一输出的模型。在解释任意黑箱模型上，SHAP 方法是迄今为止在可视化交互和特征重要性方面最全面和最主导的方法[21]。其本质上是对输入特征进行排列组合，将特征输入依次设为0，根据模型输出的变化，量化评价每个输入特征对输出结果的影响程度，方法简示图如图13 所示。输入对输出的量化影响程度用Shapley 值表示[22]，如下式。

图13 SHAP 方法简示图Fig.13 Simple diagram of SHAP method

式中：φij为基于第j个数据样本分析的第i个特征Shapley 值；xij为第j个数据样本分析的第i个特征；N为所有输入特征的集合；S为包含xij的特征集合；|S|为集合S包含的输入特征个数；n为所有输入特征的个数；f(S)为特征集合S作为模型输入时的模型输出值；为特征集合S排除特征xij后模型输出值。

经过式(11)的计算，得到了第j个数据样本中第i个特征的SHapley 值(影响程度)，从宏观上(考虑所有数据)证明局部可解释性，基于Shapley 值引入SHAP系数，定义第i个特征的SHAP 系数为：

式中：SSHAPi为第i个特征的SHAP 系数，取值[0,1]，即基于第j个数据样本分析的第i个特征SHapley 值；n为所有输入特征的个数；M为所有数据样本的个数。

SHAP 系数基于现有数据集从宏观上表征影响程度，根据Y.Nohara 等的推导[23]，各特征的Shapley 值之和等于模型当前输出和输出期望的差，表示为：

式中：EX(f(X))为模型输出的期望；该式能够在微观角度以线性模型形式进一步理解黑箱模型，辅助判断主控因素。

6 结果分析

6.1 摩阻系数预测结果

基于优选的第4 个模型，图14 展示了某两口井的预测效果。其中蓝色是真实值，黄色是部分训练集预测结果，红色是测试集预测结果。图14a 中，在训练集基本拟合较好的基础上，除了3 000 m 和4 000 m 处的预测存在极端值，测试集能够基本正确预测趋势，测试集平均误差为5.49%，最大误差为25.07%；图14b 中，数据集整体更为平滑，拟合效果较好，测试集平均误差为5.89%，最大误差为26.53%。

图14 两口井的预测效果Fig.14 Prediction results of two wells

6.2 摩阻扭矩预测与准确性验证

摩阻系数是一种间接存在的物理量，并无实测值，仅依靠现有数据集的验证不足以说明预测结果的准确性与模型的合理性。因此基于摩阻系数智能预测结果，进行摩阻扭矩实例预测分析，作进一步的准确性验证。针对同一区块的某3 口井，收集与处理数据，使用建立的LSTM 模型预测钻柱摩阻系数并代入软杆模型进行摩阻扭矩预测分析，对比大钩载荷预测值与录井实测值，验证摩阻系数预测的准确性。表5 列举了该3口井的大钩载荷智能预测值与基于不同摩阻系数经验值的大钩载荷计算值。结果表明，在不同经验值中，摩阻系数取0.25 误差最小，误差率为8.67%。与不同经验值相比，摩阻系数智能预测值最接近实测值，误差率为4.26%，显著降低了4.41%，具有较好的预测效果。因此，相比于钻井现场摩阻系数取经验值的方法，本方法预测摩阻扭矩更加准确。

表5 大钩载荷预测结果对比Table 5 Comparison of hook load prediction results

6.3 基于SHAP 的可解释性分析

SHAP 系数从宏观上反映了各输入参数对模型输出值的贡献程度。针对现有模型，计算每个特征的SHAP 系数，结果如图15 所示。模型是否具备一定程度的局部可解释性，取决于各参数的SHAP 系数大小关系是否符合真实的物理规律。根据前人的总结与分析[24]，地层数据、钻井液性能数据与部分工程数据是影响钻柱摩阻系数的主要因素。例如，在钻井过程，地层数据中不同的岩性会导致形成不同厚度的泥饼，进而影响摩阻系数；工程数据中的井眼曲率或钻柱外径越大，都会导致钻柱与井壁的接触面积越大，导致摩阻系数增大。而井深与方位角影响程度较小，只有钻进方位角与设计方位角出现偏差时，才会轻微影响摩阻力；钻井液性能同样会很大程度上影响摩阻系数，例如钻井液黏度、钻井液体系的不同会产生不同的润滑效果。图15 中，地层、钻井液性能与部分工程3 类数据基本占据了主要位置，表明这些参数对摩阻系数预测产生了主要影响作用，其中岩性影响程度最大。而井深与方位角整体排在最后，对摩阻系数预测影响最小。这一结果比较符合上文所述的摩阻系数影响因素与规律，说明模型输入输出的响应机制基本符合工程实际，能够表征管柱力学机理，具备一定程度的局部可解释性。此结果进一步提高了模型的稳定性和可靠性。

图15 各参数的SHAP 系数Fig.15 SHAP coefficient of parameters

图16 展示了该井3 882 m 处的微观解释过程，线性分析形式表征了模型在不同参数影响下由均值输出0.175 变为当前输出0.215。图中Shapley 值正负表示了该特征参数对模型输出作正或负贡献，可以辅助判断主控因素。图中该井深处摩阻系数的增大主要由钻井液塑性黏度、井眼曲率贡献，而岩性对摩阻系数产生了降低作用。因此，对于该井需要注意井眼曲率较大的井段，若监测到摩阻异常增大趋势，可及时降低钻井液塑性黏度以避免发生阻卡风险。模型仅达到局部可解释而未达到全局可解释，虽不能完全代替现场专家决策，但可以在异常发生前及时提供参考信息辅助决策，提高作业效率。

图16 3 882 m 处数据的微观解释Fig.16 Microscopic interpretation of data at 3 882 m

7 结论

a.本研究建立了机理-数据融合的钻柱摩阻扭矩预测方法，具有较好的准确性和可靠性，解决了传统方法中摩阻系数选取的盲目性和滞后性问题。

b.选取15 个特征建立了基于LSTM 的摩阻系数智能预测模型，平均误差为5.89%。同时结合软杆模型实现机理-数据融合的摩阻扭矩预测方法，钩载预测值相比经验计算值误差降低4.41%，验证了预测结果的准确性。使用SHAP 方法进行可解释性分析，结果表明，模型表征的输入特征与摩阻系数的映射关系符合管柱力学机理，具备一定可解释性，进而验证了方法的可靠性，同时能够辅助判断主控因素。

c.本研究与先前工作相比，能够实现摩阻扭矩的可靠智能预测。为摩阻扭矩分析提供了新的思路，同时也是智能钻完井理论与方法的进一步探索。但智能模型未达到全局可解释，建议进行更深一步研究，以指导参数优化进而避免风险。