基于凡尔工作点的有杆泵抽油井井下工况诊断方法与流程

本发明涉及一种油井井下工况诊断技术，具体地说是一种基于凡尔工作点的有杆泵抽油井井下工况诊断方法。

背景技术：

有杆泵是一种常用的采油设备，它具有综合成本低、结构简单和操作方便等优点。至今，国内80％以上的油井依然使用这种方法进行生产。采油过程通常是在地下数千米进行，井下的工作状况无法使用传感器进行检测，并且恶劣条件经常导致各种故障，严重影响了油田的生产效率和经济效益。所以对有杆泵的井下工作状态进行准确地诊断很有必要。

目前通过悬点示功图进行井下工作状态的判断。悬点示功图，也称地面示功图或光杆示功图，是抽油井采油现场采集的第一手资料，它是由载荷和位移变化构成的封闭曲线，可以全面地展示柱塞的井下工作状态。有经验的工程师对其进行分析可为解除油井故障、保证油井正常生产或提高油井产量等提供依据。但是这种方法不能进行实时监控并且成本也很高昂，人工诊断无法满足企业不断提高生产效率的需求。

技术实现要素：

针对现有技术井下工作状态的判断不能进行实时监控、无法满足企业不断提高生产效率的需求等不足，本发明要解决的问题是提供一种能够全面地、准确地体现井下各种工况特点的基于凡尔工作点的有杆泵抽油井井下工况诊断方法。

为解决上述技术问题，本发明采用的技术方案是：

本发明一种基于凡尔工作点的有杆泵抽油井井下工况诊断方法，包括以下步骤：

采集有杆泵抽油机地面示功图共k+1个，前k个示功图对应n种工况，用于建立工况诊断模型，最后一个为待诊断工况的示功图；

将k+1个有杆泵抽油井地面示功图进行归一化处理；

采用重心分割法将归一化后的示功图划分为四部分，左上区域、右上区域、右下区域和左下区域；

提取凡尔的四个工作点以及地面示功图的七个几何特征；

利用有杆泵抽油井地面示功图的七个几何特征特征作为输入，建立连续隐马尔科夫(chmm)的工况诊断模型，该模型用来描述工作点及示功图特征与有杆泵抽油井井下工况之间的映射关系，完成有杆泵抽油井工况诊断。

采用重心分割法将归一化后的示功图划分为四部分，左上区域、右上区域、右下区域和左下区域为：

对归一化后的示功图进行多边形分解，将示功图分为248个三角形，并根据三角形的重心和面积公式计算出每个三角形的重心和面积，然后通过下式求出示功图的重心(x^*，y^*)：

其中xi和yi为示功图上第i个点的横坐标和纵坐标；

并将归一化后的示功图按x＝x^*和y＝y^*划分为左上区域、右上区域、右下区域和左下区域四部分。

提取活塞的四个工作点以及七个示功图的几何特征为：

提取活塞的四个工作点即固定凡尔打开点、固定凡尔关闭点、游动凡尔打开点和游动凡尔关闭点，以及七个几何特征即四个工作点围成的四边形面积、左上区域缺失面积、右上区域缺失面积、右下区域缺失面积、左下区域缺失面积、固定凡尔工作距离以及游动凡尔工作距离，其中，固定凡尔打开点和游动凡尔打开点，分别为左上区域曲率的变化量最大的点和右下区曲率的变化量最大的点，对于由离散点(xi,yi)组成的示功图，采用有限差分法来计算每一点的曲率：

其中，k(t)为第t个点的曲率，x(t)为示功图第t个点的横坐标，y(t)为示功图第t个点的纵坐标；

相邻点曲率之差的绝对值来计算曲率的变化：

δk＝|k(t)-k(t-1)

采用均值滤波法对δk进行处理：

曲率变化最大的点按下式计算：

左上区域曲率的变化量最大的点和右下区曲率的变化量最大的点即为固定凡尔打开点和游动凡尔打开点，提取固定凡尔关闭点和游动凡尔关闭点，分别为右上区最大位移和左下区最小位移。

凡尔工作点以及示功图的七个几何特征，

根据凡尔的四个工作点，提取由四个工作点围成的四边形面积；

根据重心分解结果，提取左上区域缺失面积，左上区域曲线与x＝0和y＝1所围成封闭图形的面积；

根据重心分解结果，提取右上区域缺失面积，右上区域曲线与x＝1和y＝1所围成封闭图形的面积；

根据重心分解结果，提取右下区域缺失面积，右上区域曲线与x＝1和y＝0所围成封闭图形的面积；

根据重心分解结果，提取左下区域缺失面积，左上区域曲线与x＝0和y＝0所围成封闭图形的面积；

根据固定凡尔关闭点，提取固定凡尔工作距离，固定凡尔打开点到关闭点的直线距离；

根据游动凡尔关闭点，提取游动凡尔工作距离，游动凡尔打开点到关闭点的直线距离。

将有杆泵抽油井地面示功图的七个几何特征作为输入，建立连续隐马尔科夫(chmm)的工况诊断模型为：

确定chmm的模型结构；

对chmm进行训练，将每种工况类型作为隐藏状态的输入，并将不同工况类型示功图的七个几何特征作为观测向量的输入，共训练出n个chmm模型，每一个模型对应一种工况，训练是为了让模型能够根据观测序列o＝{o1,o2,…,om}来确定模型参数(a，c，μ，u，π)，其中，a为状态转移矩阵，c为混合权重，u为协方差矩阵，μ为平均向量，π为初始状态矩阵；

将待诊示功图七个几何特征作为观测向量分别输入到n个chmm中，然后分别计算出待诊断示功图与n个chmm中的相似概率，与待诊断示功图相似概率最大的模型所对应的工况即为待诊断示功图的工况。

确定chmm的模型结构为：

每一个状态对应一个工况种类，由不同状态构成的集合为：

s＝{s1,s2,…,sn}

其中n为隐藏状态的数量；

每一个观测向量oi由示功图的七个几何特征组成，i＝1，2，……m，观测序列为：

o＝{o1,o2,…,om}

其中m是观测向量的数量；

a＝{aij}为状态转移矩阵，其中aij是指在时间t到时间t+1的过程中由第i种工况转移到第j种工况的概率；

aij＝p(qt+1＝sj|qt＝si)1≤i,j≤n

根据功图及几何特征的变化为连续过程这一特点，通过混合高斯概率密度函数对观测值转移概率进行描述：

其中bn(ot)表示在时间t时第n个状态的高斯概率密度函数；c是混合权重；u是协方差矩阵；μ是平均向量；

初始状态矩阵π＝{πi}表示在初始时刻处于各种工况的概率，即；

π＝p(q0＝si)1≤i≤n

其中，p为条件概率，q0为在初始时刻的工况。

本发明具有以下有益效果及优点：

1.本发明方法根据不同工况下示功图图形的特点对其进行合理的区域划分，并对示功图的特征进行定性和定量的分析，提取出了活塞的四个工作点，以及在每个工作点周围曲线的几何特征，利用这些特征对井下的工作过程进行全面、准确地描述，并利用连续隐马尔科夫模型建立工况诊断模型，实现对杆泵抽油井井下工况的实时诊断，很大程度上少了提高了诊断效率，降低了生产成本，可保证杆泵抽油井安全、高效的进行生产。

2.本发明方法能够全面地、准确地体现井下各种工况的特点，并利用这些特征建立工况实时诊断模型，有效提高了有杆泵抽油井生产的安全性和高效性。

附图说明

图1为本发明方法流程图；

图2为本发明方法中待诊断示功图；

图3为本发明方法中归一化后的待诊断示功图；

图4为本发明方法中重心划分结果示功图；

图5为本发明方法中理论示功图；

图6为本发明方法中示功图的载荷-时间关系曲线；

图7为本发明方法中正常工况工作点分布示功图；

图8为本发明方法中气体影响工况工作点分布示功图；

图9为本发明方法中供液不足工况工作点分布示功图；

图10为本发明方法中泵下碰工况工作点分布示功图；

图11为本发明方法中泵上碰况工作点分布示功图；

图12为本发明方法中游动凡尔漏失工况工作点示功图；

图13为本发明方法中固定凡尔漏失工况工作点分布示功图；

图14为本发明方法中出砂工况工作点分布示功图；

图15为本发明方法中油稠工况工作点分布示功图。

具体实施方式

下面结合说明书附图对本发明作进一步阐述。

如图1所示，基于凡尔工作点的有杆泵抽油井井下工况诊断方法，包括以下步骤：

1)采集有杆泵抽油机地面示功图共k+1个，前k个示功图对应n种工况，用于建立工况诊断模型，最后一个为待诊断工况的示功图；

2)将k+1个有杆泵抽油井地面示功图进行归一化处理；

3)采用重心分割法将归一化后的示功图划分为四部分，左上区域、右上区域、右下区域和左下区域；

4)提取凡尔的四个工作点以及地面示功图的七个几何特征；

5)利用有杆泵抽油井地面示功图的七个几何特征特征作为输入，建立连续隐马尔科夫即chmm的工况诊断模型，该模型用来描述工作点及示功图特征与有杆泵抽油井井下工况之间的映射关系，完成有杆泵抽油井工况诊断。

本实施方式中，步骤1)从现场收集数据，共有121个示功图用于建模和诊断。前120个示功图为已知工况类型，对应9种不同工况，最后1个为待诊断示功图，其中已知工况类型的示功图有9种：15个“正常”类型的示功图、20个“气体影响”类型的示功图、25个“供液不足”类型的示功图、10个泵下碰”类型的示功图、10个“泵上碰”类型的示功图、10个“游动凡尔漏失”类型的示功图、10个“固定凡尔漏失”类型的示功图、10个“出砂”类型的示功图、10个“油稠”类型的示功图；

步骤2)将120个已知工况的示功图和1个待诊断示功图进行归一化处理；

其中i＝1，2，…，250；

yi为示功图的载荷数据；xmin和xmax为示功图的最小和最大位移；ymin和ymax为示功图的最小和最大载荷；为归一化后的横坐标；为归一化后的纵坐标；

以待诊断示功图为例，待诊断归一化前的示功图和归一化后的示功图如图2和图3所示；

步骤3)采用重心分割法对归一化后示功图进行划分；

步骤3.1)当三角形△a1a2a3的顶点为ai(xi,yi)(i＝1，2，3)，三角形的重心公式如下：

其中xg为三角形重心横坐标；yg为三角形重心纵坐标；

步骤3.2)当三角形△a1a2a3的顶点为ai(xi,yi)(i＝1，2，3)，三角形的面积公式如下：

步骤3.3)一个n条边的多边形通常可被分为n-2个三角形，由于示功图是由250个离散点组成，所以将归一化后的示功图分解为248个三角形，并根据式(3)到式(5)计算出每个三角形的重心及面积；

步骤3.4)通常上，在计算多边形重心时，多边形可以被视为一个均匀密度的平面薄片，若要求得由250个离散点构成的示功图的重心，因此式(6)和式(7)中的积分项可被转化为累加和；

其中σi为第i个三角形的面积；xgi是第i个三角形重心的横坐标；ygi是第i个三角形重心的纵坐标；s为多边形面积；

步骤3.5)综上所述，多边形的重心g(x*,y*)按下式计算：

步骤3.6)并将去噪后的示功图按x＝x^*和y＝y^*划分为左上区域、右上区域、右下区域和左下区域四部分；

以待诊断工况的地面示功图为例，划分结果如图4所示。

步骤4)中，提取固定凡尔打开点、固定凡尔关闭点、游动凡尔打开点和游动凡尔关闭点四个活塞的工作点，并提取由工作点围成的四边形面积、左上区域缺失面积、右上区域缺失面积、右下区域缺失面积、左下区域缺失面积、固定凡尔工作距离、游动凡尔工作距离；

步骤4.1)根据示功图的物理意义，如图5所示，根据b、d点的物理意义，柱塞上行到b点时抽油杆弹性变形结束，光杆下行到d点时油管、抽油杆弹性变形结束，即b、d点分别是上冲程和下冲程中载荷变化的最大点，根据上述分析，可从示功图的载荷-时间关系曲线(如图6所示)上寻找b、d点，根据上述对b、d点物理意义的分析和重心分解方法的划分结果，应从左上区域最大载荷附近和右下区的最小载荷附近寻找出载荷变化(载荷随时间的曲率变化)的最大点，计算载荷-时间关系曲率的变化量得最大值，可归为如下非线性无约束最优：

对于连续曲线的曲率可由下式求解：

x′(t)为横坐标曲线的一阶导数，x″(t)为横坐标曲线的二阶导数，y′(t)为纵坐标曲线的一阶导数，y″(t)为纵坐标曲线的二阶导数，t代表某一时刻，t＝1，2，3，…，250。

由离散点(xi,yi)组成的示功图，可采用有限差分法来计算上式中的变量；

相邻点曲率之差的绝对值来计算曲率的变化：

δk＝|k(t)-k(t-1)

为了减少示功图中噪声对曲率变化率计算造成的影响，采用均值滤波法对δk进行处理：

曲率变化最大的点可按下式计算：

根据上式分别求出左上区域曲率的变化量最大的点和右下区曲率的变化量最大的点，即为固定凡尔打开点和游动凡尔打开点；

步骤4.2)根据示功图的物理意义，如图5所示，下死点为活塞上行时的起始点a，即位移最小点，上死点是上冲程结束，下冲程开始的点c，即位移最大点，据此，从位移-时间曲线上(如图6所示)可直接确定固定凡尔关闭点和游动凡尔关闭点，游动凡尔关闭点为最小位移点，固定凡尔关闭点为最大位移点；

各种工况的工作点分布如图7到图15所示；

步骤4.3)提取由工作点围成的四边形面积；

步骤4.4)提取左上区域缺失面积，左上区域曲线与x＝0和y＝1所围成封闭图形的面积；

步骤4.5)提取右上区域缺失面积，右上区域曲线与x＝1和y＝1所围成封闭图形的面积；

步骤4.6)提取右下区域缺失面积，右上区域曲线与x＝1和y＝0所围成封闭图形的面积；

步骤4.7)提取左下区域缺失面积，左上区域曲线与x＝0和y＝0所围成封闭图形的面积；

步骤4.8)提取固定凡尔工作距离，固定凡尔打开点到关闭点的直线距离；

步骤4.9)提取游动凡尔工作距离，游动凡尔打开点到关闭点的直线距离；

表1为各种工况下的特征量化平均值

在步骤5)中，将已知工况的示功图按照工况类型分为m类，并将工况的类型作为连续隐马尔科夫模型的隐藏状态输入，提取出步骤4中所描述示功图的几何特征，作为连续应马尔科夫模型的观测序列输入，对模型进行训练，然后用训练完毕的模型对待诊断样本的工况进行诊断，具体为

步骤5.1)确定chmm的模型结构；

步骤5.1.1)每一个状态都对应一个特定的故障种类(例如气体影响、供液不足等)。由状态构成的集合为：

s＝{s1,s2,…,sn}

其中n是隐藏状态的数量；

步骤5.1.2)每一个观测向量oi是由示功图的七个几何特征组成。观测序列为：

o＝{o1,o2,…,om}

其中m是观测向量的数量；

步骤5.1.3)状态转移矩阵a＝{aij}是指各工况之间的转移概率，其中aij是指在时间t到时间t+1的过程中由第i种工况转移到第j种工况的概率；

aij＝p(qt+1＝sj|qt＝si)1≤i，j≤n

步骤5.1.4)由于石油开采是一个连续的过程，且示功图几何特征的变化也为连续过程，所以应该混合高斯概率密度函数对观测值转移概率进行描述：

其中bn(ot)表示在时间t时第n个状态的高斯概率密度函数，a为状态转移矩阵，cnm为第n种工况对应第m个观测向量的混合权重；unm为第n种工况对应第m个观测向量的协方差矩阵；μnm为第n种工况对应第m个观测向量的平均向量；ot为在时刻t的观测向量。

步骤5.1.5)最后一部分是初始状态π＝{πi}，它表示在初始时刻各个状态的概率；

π＝p(q0＝si)1≤i≤n

其中，p为条件概率，q0为在初始时刻的工况。

步骤5.2)对chmm进行训练，将每种工况类型作为隐藏状态的输入，并将不同工况类型示功图的七个几何特征作为观测向量的输入，共训练出n个chmm模型，每一个模型对应一种工况，训练是为了让模型能够根据观测序列o＝{o1,o2,…,om}的特征来确定模型参数(a，c，μ，u，π)；

步骤5.2.1)chmm参数更新的重估公式是根据baum-welch进行的，定义下列变量：

其中αt为前向变量；βt为后向变量；γt,n表示在时刻t时状态qt为si的概率；ξt,ij表示在时刻t时状态qt为si并在时刻t+1时状态sj的概率；γt,nm表示在时刻t且状态为qt时si混合数为xnm的概率；o为观测序列；λ为模型参数的集合，即λ＝(a，c，μ，u，π)。

步骤5.2.2)根据baum-welch和式(11)，模型的参数重估公式如下所示：

πi＝γ1(i)

步骤5.3)将待诊示功图的七个几何特征作为观测向量分别输入到n个chmm中，然后分别计算出待诊断示功图与n个chmm中的相似概率，与待诊断示功图相似概率最大的模型所对应的工况即为待诊断示功图的工况，为了弥补相似概率在计算过程易出现下溢的不足，因此采用数据压缩(取对数)来处理这个不足；

p＝lg|p(o|λ)|

此时p的值越接近0代表相似概率越大。

本实施方式中，待诊断故障类型经过重心分割后提取出了四个工作点及其他几何特征，并将这些提取后的结果输入到连续隐马尔科夫模型工况诊断模型中，相似概率分布如表2所示；

表2待诊断示功图相似概率分布

如表2结果所示，诊断结果为“供液不足”。