0 引言

在钻井过程中，尤其是在钻大位移井的情况下，钻柱与井壁、钻头与井底的非线性接触会引起钻柱的黏滑振动^[1]。在钻柱黏滞阶段，由于非线性摩擦阻力的存在，钻头会停止转动，随着顶部驱动装置如转盘或顶驱旋转，钻柱中开始积累扭转弹性势能。当钻柱顶部驱动扭矩增大到大于非线性摩擦扭矩时，钻头开始加速旋转，转速迅速增加。当到达某一时刻，在非线性摩擦阻力矩作用下，钻头再次发生黏滞而停止转动。当上述过程周期性发生时，就形成了钻柱的黏滑振动^[2]。钻柱的黏滑振动往往会给钻井过程带来许多不利影响。例如，黏滑振动会缩短钻头的生命周期，同时会在钻柱当中产生周期性的交变应力，进而可能导致钻柱发生疲劳破坏引发钻井事故，危害钻井工作人员的生命安全，极大地增加钻井成本。另外，黏滑振动还会显著降低钻头的破岩效率以及影响井下随钻测量设备的测量准确性，甚至会对昂贵的随钻测量设备造成损害。

为了有效减小黏滑振动产生的危害，国内外学者从理论分析、现场试验以及数值模拟等角度做了大量研究。汤历平等^[3]通过建立钻柱黏滑振动力学模型，研究了黏滑振动的自激振动特性。付蒙等^[4]通过仿真研究了钻柱系统的黏滑振动特性与振动产生的机理。杨全近等提出了一种利用井下MWD磁传感器数据来分析钻具黏滑现象的新方法。黄根炉等^[5]提出了一种通过拉普拉斯变换解非零极点方程判断钻柱扭转振动衰减快慢的新方法, 研究了顶部负扭矩减振。A.A.BESAISOW等^[6]建立了钻柱黏滑振动的力学模型，讨论了阻尼、转盘惯性、钻柱转速控制系统对黏滑振动的影响。Y.Q.LIN等^[7]通过黏滑振动的数值模拟，解释了钻柱的大幅扭转振动。J.F.BRETT^[8]建立了一个模型，证明了金刚石钻头的固有特性决定其能够引发黏滑振动，并将试验结果和仿真结果进行了对比。N.CHALLAMEL等^[9]同时考虑了现场钻井数据和现场岩石力学特性，提出了减轻黏滑振动的措施。T.RICHARD等^[10]利用一种离散模型研究了刮刀钻头的黏滑振动。D.BRESCH-PIETRI和G.W.HALSEY等^[11-12]利用PI型速度控制器有效地抑制了钻柱的黏滑振动。另外，扭矩负反馈的方法被提出抑制钻柱的黏滑振动。

目前还未发现一种简单、有效、成本低和实时性高的黏滑振动等级评估方法。因此，笔者提出一种新颖的基于因子分析(FA)和支持向量机(SVM)的黏滑振动等级评估方法。该方法不依赖井下随钻测量工具，直接通过井口扭矩数据对黏滑振动等级进行评估，避免了由于信号通过钻井液从井底向井口传输造成的时间延迟，提高了信号处理的实时性。虽然井口扭矩数据可通过井口扭矩测量传感器实时获取，但为了获取更多的扭矩样本数据进行机器学习，首先建立了钻柱黏滑振动力学模型，然后采用仿真计算方法求得井口扭矩数据。通过对大量的扭矩数据样本进行时域和频域分析，得到了扭矩信号的一系列特征值。采用因子分析方法可减少扭矩特征的冗余信息，获取了最终的用于SVM监督式机器学习的输入特征向量。最后通过SVM对黏滑振动等级进行评估，并通过混淆矩阵和准确率评价了支持向量机分类器的分类精度。

1 力学模型与求解

为建立钻柱黏滑振动的力学模型，做如下基本假设：①井眼为理想直井。②钻头由转盘或顶驱驱动。③将钻杆和底部钻具组合(BHA)看作具有集中参数的质量块，并由具有扭转刚度的弹簧和黏性阻尼与转盘或顶驱连接。④忽略钻柱的其他动力学行为对黏滑振动的影响。⑤井口输入转速恒定^[2]。根据以上假设，建立钻柱黏滑振动力学模型，如图 1所示。

对该力学模型进行动力学分析，得到钻柱的动力学方程如下:

(1)

式中：J为钻杆和底部钻具组合相对于钻柱轴线的转动惯量；c为阻尼系数；k为钻杆的扭转刚度； 表示钻头与井底的非线性摩擦阻力矩(在此考虑到钻井液的润滑作用，忽略钻柱与井壁的摩擦阻力矩作用)；θ为钻头旋转角度；w₀为转盘旋转的角速度。

为了对式(1)进行求解，将其转化为动力学方程标准化形式：

(2)

式中：ξ为阻尼比；p为钻柱的固有频率。

其余各变量定义如下：

(3)

将式(3)代入式(2)，采用龙格库塔数值积分方法对式(2)进行求解，求得钻头转速及井口扭矩随时间的变化，并对相应数据进行记录。黏性阻尼的力学模型以及式(2)中的关键参数见参考文献[7]。通过改变式(2)当中的阻尼比、固有频率和转盘转速3个参数，可以产生不同程度的黏滑振动。E.W.ROBNETT等^[13]介绍了一种通过井底随钻测量参数评价钻柱黏滑振动程度的方法，他们认为在至少一个黏滑振动周期内，钻柱转速的最大值、最小值之差与转速平均值的比值表征了黏滑振动严重程度。通过该方法可以求出不同参数导致的不同程度的黏滑振动等级，并将该振动等级下的扭矩样本数据与黏滑振动等级对应，对其进行有监督的机器学习，进而可以通过井口扭矩数据评估井下钻柱的黏滑振动等级。

2 扭矩数据时域与频域特征值提取

通过求解动力学方程，可以获得不同黏滑振动等级下的扭矩样本。仿真结果与文献结果对比如图 2所示。将通过龙格库塔方法求出的扭矩曲线与参考文献[7]中的扭矩曲线进行对比，此时式(2)取与参考文献相同的参数，即阻尼比为0.03，固有频率为0.77 Hz，钻盘转速为4.19 rad/s。对比发现，通过龙格库塔方法求得的扭矩数据可靠，可用于进一步分析。

对求得的扭矩时域信号进行时域和频域分析，可以求得扭矩信号的特征值。通常情况下，信号的时域特征值主要有最大值、最小值、平均值、有效值、方根幅值、均方差、标准差、偏度、偏态因数、峭度、峰态、波形因数、脉冲因数、峰值因数和裕度因数等。其中，前9个为有量纲特征，后5个无量纲特征。其相关定义可参考信号处理相关书籍。频域信号的主要特征有重心频率F_c、均方频率F_MS、均方根频率R_F、频率方差V_F和频率标准差R_V等，其定义如表 1所示。表中T(f)表示信号的功率谱。通常情况下，信号的功率谱求解方法主要有周期图法、间接法、Bartlett法和Welch法等。其中，Welch法在标准周期图法和Bartlett法的基础上做了改进，现已成为求解信号功率谱的主要方法。因此，文中选用Welch法求取信号功率谱，进而编程求解出5种频率特征参数。

3 因子分析 3.1 因子分析方法

因子分析假设观察到的变量是一些潜在的隐含的因子的线性组合。在这些隐含的因子当中，有一些因子被认为是两个或者更多变量共同的因子，另一些因子则被认为是各变量的独有因子。这些独有的因子被进一步假定为(至少在探索性因子分析中)其保持互成正交关系。因此，认为这些独有的因子对变量之间的共变关系没有贡献。换言之，只有公共因子对观察到的变量之间的共变关系有贡献^[14]。因子分析的一般模型如下^[15]：

(4)

式中：f₁、f₂、……、f_m为m个公共因子，ε_i为x_i特有的特殊因子，它们均为不可观测隐含变量。

将式(4)写成矩阵的形式，则有：

(5)

式中：A为因子载荷矩阵，f为公共因子向量，α为样本均值向量，ε为特殊因子向量。

考虑到公共因子对变量之间的关系起主要作用，因此只保留公共因子，这样可以在保存原有数据样本大部分信息的情况下对数据样本进行降维处理，能够有效地消除信息冗余。因子分析模型的求解主要有主成分法、主因子法以及最大似然法。采用MATLAB软件自带的因子分析工具函数，该函数采用最大似然法，通过调用该因子分析函数，可得到因子模型的载荷矩阵、方差的最大似然估计、因子旋转矩阵、包含模型信息的检验结构体变量以及因子得分矩阵。

3.2 扭矩特征信号降维

通过对扭矩样本的时频与频域分析，对每组扭矩样本信号都可以提取19个时域与频域的特征值。对动力学方程(2)求解，并改变其中的3个关键参数阻尼比ξ、固有频率p以及转盘转速w₀，可以得到不同黏滑振动等级的扭矩样本。阻尼比的改变范围为0.01~0.30，改变步长为0.01;固有频率的改变范围为0.5~2.5 Hz，改变步长为0.1 Hz；转盘转速的改变范围为1~50 rad/s，改变步长为0.5 rad/s。编写龙格库塔数值积分程序进行分析计算，累计得到62 370组不同黏滑振动等级的扭矩样本。一方面，考虑到每组扭矩样本的所有19个特征值当中，其内部存在一定的相关性，即存在信息冗余；另一方面，如果将所有的19个时域与频域特征全部送入SVM网络进行监督式机器学习，计算量大，同时计算效率较低。因此，需采用因子分析对62 370×19扭矩特征值矩阵进行特征降维处理。因子分析之前，需要对扭矩矩阵每维特征值进行标准化处理。标准化处理的一般方法如下：

(6)

式中：x_i表示每维特征值中的第i个数据，μ为每维特征值的均值，σ为每维特征值的方差。

标准化处理之后，对扭矩矩阵进行降维处理。由于因子分析之前并不知道取几个公共因子才能保留大部分的样本信息，所以先取较多的公共因子进行数据降维处理，然后根据分析结果调整公共因子数量，最终确定最优的公共因子数目。分析后先取9个公共因子进行计算，也就是取大概1/2特征值数目的公共因子。编程计算此情况下的因子载荷矩阵，根据因子载荷矩阵，计算每个公共因子的贡献率以及累计贡献率。贡献率计算结果如图 3所示。

从图 3可以看出，前面5个公共因子对原始方差的贡献率分别为28.52%、26.38%、16.05%、9.23%和8.85%，该5个公共因子的累计方差贡献率达到了89.03%，并且随着公共因子数目的进一步增加，因子累计贡献率变化十分缓慢。这说明该因子分析模型中所取的公共因子数目可以进一步减少，只考虑前5个公共因子应该最合适。对该扭矩数据样本重新进行因子分析，最终确定公共因子数量为5，经过计算，此时的因子载荷矩阵如表 2所示。通过因子分析的降维处理，最终将62 370×19的扭矩特征值矩阵降维为62 370×5的扭矩特征值矩阵，有效消除了扭矩特征值当中的信息冗余，同时可显著提高机器学习效率，节省计算机硬件资源，为下一步进行SVM机器监督学习做好数据准备工作。

4 支持向量机黏滑振动等级评估 4.1 支持向量机

支持向量机最先由Vapnik提出，它是一种通用的前馈型神经网络，可用于解决线性与非线性的数据分类和拟合问题^[16]。线性可分超平面如图 4所示。支持向量机的主要思想是构建一个决策超平面(图 4中的黑实线)，使距离该平面最近的不同类型的样本距离该平面的距离达到最大。图 4中点画线是将最大距离标准化之后的结果，位于其上面的点成为支持向量，支持向量在SVM分类中起主要作用。进行支持向量机监督学习，首先需要把训练数据集从一个低维空间映射到一个高维特征空间，其目的是把原来线性不可分的关键问题转化为线性可分的问题^[17]。支持向量机在求解过程中采用了基于KKT条件的拉格朗日乘子法，进而引入了核函数，大大简化了计算过程，常用的核函数主要有高斯核函数、多项式核函数、线性核函数以及双层感知器核函数等。

4.2 黏滑振动等级评估

通过因子分析降维后的扭矩数据样本被分为两部分：一部分作为训练集，另一部分作为测试集。70%的样本数据用于对SVM分类器进行训练，剩下的30%用于评估分类器的性能表现。本文的黏滑振动等级分为4类，故该问题是一个多分类的监督式机器学习问题。多分类的SVM及机器学习问题主要有两种解决办法：一种是直接法，该方法计算复杂度高，实施较困难；第二种是间接法，它主要通过构造多个二分类SVM进行实现，根据其具体实现方法又分为一对多法、一对一法、层次SVM以及无环图SVM等方法^[18]。本文的SVM机器学习采用了台湾大学林智仁教授等人开发的LIBSVM机器学习工具箱，其内部分类过程采用间接法中的一对一法。

核函数的选择对机器学习预测准确率有很大的影响，因此采用不同的核函数依次对SVM进行训练及测试，通过对比可以确定最优的核函数选择结果。不同核函数下SVM测试结果如表 3所示。

从表 3可以看到，对于该黏滑振动等级的分类，基于径向基核函数的支持向量机分类效果最好，测试集预测准确率达到了80%以上。图 5显示了采用径向基核函数时支持向量机钻柱黏滑振动等级分类结果。为了更加直观地显示类别错误分配的结果，采用混淆矩阵的方法。表 4显示了黏滑振动等级预测结果与实际振动等级的对比情况。从表 4中可以看到，4个黏滑振动等级的SVM预测结果准确率分别达到了81.5%、81.2%、70.1%和89.5%, 整体预测正确率超过80%。因此，可以认为该方法可以较为有效地对黏滑振动等级进行判别，可以将其作为调整钻井参数的有力参考。可以将该支持向量机的等级分类模块作为钻井专家系统的一个子模块，当黏滑振动等级超过给定值，就可以提醒司钻及时调整钻井参数，避免长时间的强烈黏滑振动造成钻井工具损坏以及钻头寿命缩短，从而有效提高机械钻速，节省钻井成本，将由黏滑振动造成的断钻具风险降到最低，保证钻井作业现场工作人员的生命财产安全。

5 结论

(1) 建立了钻柱黏滑振动力学模型，采用仿真计算方法求得井口扭矩数据。通过对大量的扭矩数据样本进行时域和频域分析，得到了扭矩信号的一系列特征值。采用因子分析方法可减少扭矩特征的冗余信息，获取最终的用于支持向量机监督式机器学习的输入特征向量。

(2) 基于井口扭矩信号的SVM黏滑振动等级预测方法的整体预测精度超过80%，能够较准确地对黏滑振动强度等级进行预测。因此该方法是一个有效的黏滑振动等级分类方法，运用该方法可以减小对井下随钻测量工具的依赖，消除以钻井液为载体的信号传输延迟，节省黏滑振动等级诊断的成本，同时提高实时性。

(3) 将该方法与钻井专家系统相结合，并作为其子模块，就可以协助井口司钻根据钻头黏滑振动等级及时调整钻井参数，如转盘或顶驱对钻柱的输入转速、钻头钻压等，减轻甚至消除黏滑振动，提高钻井效率与安全性。