CN104481496A - 一种有杆泵抽油井故障诊断方法 - Google Patents

Abstract

本发明涉及一种有杆泵抽油井故障诊断方法，获取已知和待诊断有杆泵抽油井地面示功图，将有杆泵抽油井地面示功图转化为有杆泵抽油井泵示功图，采用16方向链码对归一化处理后的泵示功图进行边界链码重画，将边界链码重画后的泵示功图进行区域划分，并提取各区域的特征向量，利用基于粒子群优化的谱聚类算法对量化后的泵示功图特征向量进行聚类，完成有杆泵抽油井故障诊断，本发明采用16方向链码的数量对各个泵示功图特征向量进行量化能够有效、细致地描述图形特征，有杆泵抽油井故障诊断是利用基于粒子群优化的谱聚类算法对量化后的泵示功图特征向量进行聚类，不依赖训练样本，可以提高诊断的性能。

Description

一种有杆泵抽油井故障诊断方法

技术领域

本发明属于抽油井故障诊断领域，特别涉及一种有杆泵抽油井故障诊断方法。

背景技术

有杆泵抽油井泵示功图在油田生产中具有十分重要的作用，它是判断抽油系统井下工作状态的主要依据。泵示功图可以用来进行井下生产状况的故障诊断，也可以用来计算井口的出液量，还可以用来计算井下动液面的位置。传统生产中，一般依靠技术人员对采集的泵示功图进行人工分析，从而判断抽油系统的实际工作状况。这种人工方式工作周期长，受主观因素影响较大，越来越无法满足企业不断提高生产效率的需求。

随着计算机、通信、电子技术等领域的不断进步，利用计算机取代人力劳动方式已成为国内外油田生产的发展趋势。鉴于泵示功图在生产中的重要性，采用计算机进行计算和分析已成为该领域研究的重点和难点。利用计算机处理泵示功图，关键在于提取其显著的特征，识别越准确，计算越正确。目前，对于泵示功图的特征提取方法，主要有面积法、匹配法、网格法、矢量法等。但是这些方法都不能很好地结合实际生产过程的特点，从而使得到的特征向量不能准确、细致地反映不同过程的显著特征。

有杆泵抽油井的抽油泵一般工作在数千米的地下，生产环境恶劣，故障发生率很高，采用计算机实现诊断是很有意义的。现有的诊断方法主要以有监督学习为主，根据生产经验得到不同故障类型的训练样本集，然后由这些训练样本训练出不同故障类型的训练集模型，从而实现待诊断样本的诊断。这种有监督学习方法主要存在以下两方面问题：一是诊断的准确率很大程度上取决于所构成训练样本集的样本特征向量是否能够准确反映样本的典型特征；二是这种有监督方式主要依靠人为划分训练样本集，容易受到主观因素的影响。

发明内容

针对现有技术的不足，本发明提供一种有杆泵抽油井故障诊断方法。

本发明的技术方案是：

一种有杆泵抽油井故障诊断方法，包括以下步骤：

步骤1：获取M个已知故障类型的有杆泵抽油井地面示功图和1个待诊断故障类型的有杆泵抽油井地面示功图；

步骤2：将M+1个有杆泵抽油井地面示功图转化为有杆泵抽油井泵示功图，并对有杆泵抽油井泵示功图进行归一化处理，得到归一化处理后的泵示功图；

步骤3：采用16方向链码对归一化处理后的泵示功图进行边界链码重画，得到边界链码重画后的泵示功图，其中，边界链码的类型由z_t表示，t＝1、2、…16，NUM_t表示边界链码类型为z_t的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数；

步骤4：将边界链码重画后的泵示功图进行区域划分，并提取各区域的特征向量，各特征向量以16方向链码形式表示；划分的区域包括：左上区、右上区、左下区和右下区；所述的特征向量包括：左上区欠缺程度特征向量、左上区快速加载程度特征向量、右上区突出程度特征向量、右上区欠缺程度特征向量、右上区快速卸载程度特征向量、左下区突出程度特征向量、左下区欠缺程度特征向量、右下区欠缺程度特征向量和右下区快速卸载程度特征向量；

步骤4.1：确定有杆泵抽油井泵示功图的上冲程和下冲程：将有杆泵抽油井泵示功图曲线的横坐标最大的点作为上冲程的结束点，则有杆泵抽油井泵示功图的上冲程为有杆泵抽油井泵示功图曲线第一个数据点到上冲程结束的数据点，有杆泵抽油井泵示功图的下冲程为有杆泵抽油井泵示功图曲线上冲程结束点的下一个数据点到最后一个数据点；

步骤4.2：提取有杆泵抽油井泵示功图各区域的泵示功图特征向量，包括左上区欠缺程度特征向量、左上区快速加载程度特征向量、右上区突出程度特征向量、右上区欠缺程度特征向量、右上区快速卸载程度特征向量、左下区突出程度特征向量、左下区欠缺程度特征向量、右下区欠缺程度特征向量和右下区快速卸载程度特征向量；

左上区欠缺程度特征向量表示为：

$T_1=\frac{{\mathrm{NUM}}_1+3*{\mathrm{NUM}}_2+2*{\mathrm{NUM}}_3+{\mathrm{NUM}}_4-2*{\mathrm{NUM}}_5}{\mathrm{NUM}};$

左上区快速加载程度特征向量表示为：

$T_2=\frac{{\mathrm{NUM}}_3+2*{\mathrm{NUM}}_4+{3*\mathrm{NUM}}_5-{\mathrm{NUM}}_1-2*{\mathrm{NUM}}_2}{\mathrm{NUM}};$

右上区突出程度特征向量表示为：

其中，为有杆泵抽油井泵示功图的上冲程数据点纵坐标的平均值，为边界链码的类型为z₃、z₄和z₅的边界链码重画后的泵示功图的数据点纵坐标的平均值；

右上区欠缺程度特征向量表示为：

$T_4=\frac{3*{\mathrm{NUM}}_{15}+2*{\mathrm{NUM}}_{16}+{\mathrm{NUM}}_{14}-3*{\mathrm{NUM}}_1-{\mathrm{NUM}}_2}{\mathrm{NUM}};$

右上区快速卸载程度特征向量表示为：

$T_5=\frac{3*{\mathrm{NUM}}_{13}+{\mathrm{NUM}}_{12}-2*{\mathrm{NUM}}_{11}}{\mathrm{NUM}};$

左下区突出程度特征向量表示为：

其中，为有杆泵抽油井泵示功图的下冲程数据点纵坐标的平均值，为边界链码的类型为z₁₁、z₁₂、z₁₃、z₁₄、z₁₅和z₁₆的边界链码重画后的泵示功图的数据点纵坐标的平均值；

左下区欠缺程度特征向量表示为：

$T_7=\frac{3*{\mathrm{NUM}}_6+2*{\mathrm{NUM}}_7-{\mathrm{NUM}}_8-2*{\mathrm{NUM}}_9-{\mathrm{NUM}}_{10}}{\mathrm{NUM}};$

右下区欠缺程度特征向量表示为：

$T_8=\frac{2*{\mathrm{NUM}}_9+3*{\mathrm{NUM}}_{10}+{\mathrm{NUM}}_{11}-2*{\mathrm{NUM}}_{12}-3*{\mathrm{NUM}}_{13}}{\mathrm{NUM}};$

右下区快速卸载程度特征向量表示为：

$T_9=\frac{3*{\mathrm{NUM}}_9+{\mathrm{NUM}}_{12}+2*{\mathrm{NUM}}_{13}-2*{\mathrm{NUM}}_{11}}{\mathrm{NUM}};$

步骤4.3：将各区域的泵示功图特征向量进行量化，得到量化后的泵示功图特征向量；

步骤5：利用基于粒子群优化的谱聚类算法对量化后的泵示功图特征向量进行聚类，完成有杆泵抽油井故障诊断；

步骤5.1：设将M+1个泵示功图特征向量聚类为k类，σ为尺度参数，其中，σ∈[0，1]；

步骤5.2：将聚类数k和尺度参数σ作为粒子进行聚类优化，得到最优的聚类数k和最优的尺度参数σ；

步骤5.2.1：将聚类数k和尺度参数σ作为粒子，随机产生若干个粒子，初始化粒子的速度和位置；

步骤5.2.2：M+1个量化后的泵示功图特征向量由d₁、d₂…d_M+1表示，由高斯函数s(d_p，d_q)＝exp(-||d_p-d_q||²/(2σ²))计算任意两个量化后的泵示功图特征向量的相似度，各相似度值构成相似度矩阵S∈R^(M+1)×(M+1)，其中，p、q＝1、2…M+1；

步骤5.2.3：将每个量化后的泵示功图特征向量看作顶点V，根据量化后的泵示功图特征向量间的相似度对各顶点间的边E赋权重值W，得到一个基于量化后的泵示功图特征向量相似度的无向加权图G＝(V，E)；

步骤5.2.4：将相似度矩阵S的每行元素相加，即得到对应顶点的度，以所有顶点的度作为对角元素构成的对角矩阵，即度矩阵D；

步骤5.2.5：根据度矩阵、各顶点间的边的权重及单位矩阵，建立拉普拉斯矩阵L_rw＝I-D^-1W，其中，I为单位矩阵；

步骤5.2.6：计算拉普拉斯矩阵L_rw前k个最小的特征值对应的特征向量，k个特征向量组成特征向量矩阵U；

步骤5.2.7：将特征向量矩阵U的第g行向量(e_g)_{g＝1，2，…，M+1}作为聚类数据集；

步骤5.2.8：用CritC指标函数的倒数作为各个粒子的适应度值；

步骤5.2.9：若当前粒子的适应度值小于其历史适应度值，则将当前粒子的适应度值作为其历史适应度值；

步骤5.2.10：设定粒子的适应度值终止条件，若当前粒子的适应度值小于其终止条件，则当前的聚类数k和尺度参数σ为最优的聚类数k和最优的尺度参数σ，当前的聚类数据集(e_g)_{g＝1，2，…，M+1}为优化后的聚类数据集(e_g)_{g＝1，2，…，M+1}，否则，更新粒子的速度和位置，返回步骤5.2.3：

步骤5.3：利用最优的聚类数k和最优的尺度参数σ进行聚类优化，得到最优的聚类中心；

步骤5.4：根据最优的聚类中心将M+1个泵示功图特征向量聚类为k类，得到最终的聚类结果：H₁、H₂…H_k，若待诊断故障类型的泵示功图特征向量与某知故障类型的泵示功图特征向量聚为一类，则两者为同类故障，即得到有杆泵抽油井故障诊断结果，若待诊断故障类型的泵示功图特征向量单独聚为一类，则该待诊断故障类型为新故障类型。

本发明的有益效果是：泵示功图曲线能够很好地描述曲线的轮廓特征；所提取的泵示功图曲线不同区域的典型泵示功图特征向量能够精确地描述泵示功图特征，很好地反映实际生产特点，从而提高泵示功图的识别准确率；采用16方向链码的数量对各个泵示功图特征向量进行量化能够有效、细致地描述图形特征；利用基于粒子群优化的谱聚类算法对量化后的泵示功图特征向量进行聚类，完成有杆泵抽油井故障诊断，不依赖训练样本，可以提高诊断的性能。

附图说明

图1为本发明具体实施方式中一种有杆泵抽油井故障诊断方法的流程图；

图2为本发明具体实施方式中的有杆泵抽油井泵示功图；

图3为本发明具体实施方式中的归一化处理后的泵示功图；

图4为本发明具体实施方式中采用的16个方向链码示意图；

图5为本发明具体实施方式中16个方向链码按逆时针顺序构成的区域边界示意图；

图6为本发明具体实施方式中采用16方向链码进行边界链码重画后的泵示功图；

图7为本发明具体实施方式中左上区欠缺程度特征链码示意图；

图8为本发明具体实施方式中左上区快速加载程度特征链码示意图；

图9为本发明具体实施方式中右上区欠缺程度特征链码示意图；

图10为本发明具体实施方式中右上区快速卸载程度特征链码示意图；

图11为本发明具体实施方式中左下区欠缺程度特征链码示意图；

图12为本发明具体实施方式中右下区欠缺程度特征链码示意图；

图13为本发明具体实施方式中右下区快速卸载程度特征链码示意图；

图14为本发明具体实施方式中对量化后的泵示功图特征向量进行聚类为k类后的示意图。

具体实施方式

下面结合附图对本发明具体实施方式加以详细的说明。

一种有杆泵抽油井故障诊断方法，如图1所示，包括以下步骤：

步骤1：获取M个已知故障类型的有杆泵抽油井地面示功图和1个待诊断故障类型的有杆泵抽油井地面示功图。

本实施方式中，选用108个已知故障类型的有杆泵抽油井地面示功图和1个待诊断故障类型的有杆泵抽油井地面示功图，其中，已知故障类型的有杆泵抽油井地面示功图有5种，将“正常”也看作一种故障类型，则故障类型分别为：24个“正常”类型的有杆泵抽油井地面示功图、20个“气体影响”类型的有杆泵抽油井地面示功图、24个“供液不足”类型的有杆泵抽油井地面示功图、20个“泵上碰”类型的有杆泵抽油井地面示功图、20个“泵下碰”类型的有杆泵抽油井地面示功图。

步骤2：将109个有杆泵抽油井地面示功图转化为有杆泵抽油井泵示功图，并对有杆泵抽油井泵示功图进行归一化处理，得到归一化处理后的泵示功图。

建立用来描述有杆泵抽油井抽油杆运动和应力传播的一维粘滞阻尼波动方程，采用傅里叶系数法求解该一维粘滞阻尼波动方程，得到抽油杆任意截面处的位移和载荷，进而得到有杆泵抽油井抽油泵处的位移和载荷，从而确定有杆泵抽油井泵示功图。

如图2所示，有杆泵抽油井泵示功图的每一个数据都由离散数据对(x_i，y_i)构成，其中：i＝1，2，…，N，N为构成有杆泵抽油井泵示功图的离散数据对数；x_i为有杆泵抽油井泵示功图位移数据，y_i为有杆泵抽油井泵示功图载荷数据。

对有杆泵抽油井泵示功图进行归一化处理，得到归一化处理后的泵示功图。

对有杆泵抽油井泵示功图按如下归一化公式进行归一化处理：

$\stackrel{\‾}{x_i}=\frac{x_i-x_{\min }}{x_{\max }-x_{\min }}---\left(1\right)$

$\stackrel{\‾}{y_i}=\frac{y_i-y_{\min }}{y_{\max }-y_{\min }}---\left(2\right)$

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

x_i为有杆泵抽油井泵示功图位移数据；

y_i为有杆泵抽油井泵示功图载荷数据；

为归一化处理后的泵示功图的位移数据；

为归一化处理后的泵示功图的载荷数据。

归一化处理后的泵示功图如图3所示。

步骤3：采用16方向链码对归一化处理后的泵示功图进行边界链码重画，得到边界链码重画后的泵示功图，其中，边界链码的类型由z_t表示，t＝1、2、…16，NUM_t表示边界链码类型为z_t的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

如图4所示，16方向链码的角度分别为：0°、22.5°、45°、67.5°、90°、112.5°、135°、157.5°、180°、202.5°、225°、247.5°、270°、292.5°、315°、337.5°。

边界链码的类型由z_t表示，t＝1、2、…16。

采用的16方向链码分别为：0°(1)、22.5°(2)、45°(3)、67.5°(4)、90°(5)、112.5°(6)、135°(7)、157.5°(8)、180°(9)、202.5°(10)、225°(11)、247.5°(12)、270°(13)、292.5°(14)、315°(15)、337.5°(16)。

16方向链码在0°、90°、180°、270°的长度为1，在45°、135°、225°、315°的长度为在22.5°、67.5°、12.5°、157.5°、202.5°、247.5°、292.5°、337.5°的长度为

如图5所示，为根据16个方向链码重画泵示功图的曲线轮廓，其按照逆时针方向顺序构成的区域边界示意图。

由ω_i表示第i个离散数据对到第i+1个离散数据对的角度，其计算公式如下：

${\mathrm{\ω}}_i=\arctan \frac{\stackrel{\‾}{y_{i+1}}-\stackrel{\‾}{y_i}}{\stackrel{\‾}{x_{i+1}}-\stackrel{\‾}{x_i}}---\left(3\right)$

那么，边界链码的类型z_t可以由下式表示：

$z_t=\mathrm{Arg}\underset{t}{\min }|({\mathrm{\ω}}_i+{\mathrm{\ω}}_{i-1}-z_{t-1})-\frac{\mathrm{\π}}{8}t|---\left(4\right)$

其中，表示取计算的最小值的类型，t＝1、2、…16。

将归一化处理后的泵示功图上的两个相邻的数据点转换成边界链码，每一个方向的链码表示前一个数据点指向后一个数据点的角度。将归一化处理后的泵示功图上所有方向的边界链码连接成曲线，得到边界链码重画后的泵示功图。

采用16方向链码进行边界链码重画后的泵示功图如图6所示。

NUM_t表示边界链码类型为z_t的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

步骤4：将边界链码重画后的泵示功图进行区域划分，并提取各区域的特征向量，各特征向量以16方向链码形式表示；划分的区域包括：左上区、右上区、左下区和右下区；所述的特征向量包括：左上区欠缺程度特征向量、左上区快速加载程度特征向量、右上区突出程度特征向量、右上区欠缺程度特征向量、右上区快速卸载程度特征向量、左下区突出程度特征向量、左下区欠缺程度特征向量、右下区欠缺程度特征向量和右下区快速卸载程度特征向量。

步骤4.1：确定有杆泵抽油井泵示功图的上冲程和下冲程：将有杆泵抽油井泵示功图曲线的横坐标最大的点作为上冲程的结束点，则有杆泵抽油井泵示功图的上冲程为有杆泵抽油井泵示功图曲线第一个数据点到上冲程结束的数据点，有杆泵抽油井泵示功图的下冲程为有杆泵抽油井泵示功图曲线上冲程结束点的下一个数据点到最后一个数据点。

为有杆泵抽油井泵示功图的上冲程数据点纵坐标的平均值；

为有杆泵抽油井泵示功图的下冲程数据点纵坐标的平均值。

步骤4.2：提取有杆泵抽油井泵示功图各区域的泵示功图特征向量，包括左上区欠缺程度特征向量、左上区快速加载程度特征向量、右上区突出程度特征向量、右上区欠缺程度特征向量、右上区快速卸载程度特征向量、左下区突出程度特征向量、左下区欠缺程度特征向量、右下区欠缺程度特征向量和右下区快速卸载程度特征向量。

由于系统加载时变得缓慢，使有杆泵抽油井泵示功图曲线的增载线后半部分曲率变化较小，造成有杆泵抽油井泵示功图的图形左上区欠缺，因此左上区欠缺程度特征向量表示为：

$T_1=\frac{{\mathrm{NUM}}_1+3*{\mathrm{NUM}}_2+2*{\mathrm{NUM}}_3+{\mathrm{NUM}}_4-2*{\mathrm{NUM}}_5}{\mathrm{NUM}}---\left(5\right)$

其中，NUM₁表示边界链码类型为z₁的边界链码个数，NUM₂表示边界链码类型为z₂的边界链码个数，NUM₃表示边界链码类型为z₃的边界链码个数，NUM₄表示边界链码类型为z₄的边界链码个数，NUM₅表示边界链码类型为z₅的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

左上区欠缺程度特征链码示意图如图7所示，其中，实线代表系统正常加载，虚线代表系统加载变缓。

由于系统突然急剧加载，使有杆泵抽油井泵示功图曲线的增载线前半部分曲率变化较大，因此，左上区快速加载程度特征向量表示为：

$T_2=\frac{{\mathrm{NUM}}_3+2*{\mathrm{NUM}}_4+{3*\mathrm{NUM}}_5-{\mathrm{NUM}}_1-2*{\mathrm{NUM}}_2}{\mathrm{NUM}}---\left(6\right)$

其中，NUM₁表示边界链码类型为z₁的边界链码个数，NUM₂表示边界链码类型为z₂的边界链码个数，NUM₃表示边界链码类型为z₃的边界链码个数，NUM₄表示边界链码类型为z₄的边界链码个数，NUM₅表示边界链码类型为z₅的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

左上区快速加载程度特征链码示意图如图8所示，其中，实线代表系统正常加载，虚线代表系统急剧加载。

在上冲程末端由于系统内部相关器件发生碰撞，使载荷突然增大，造成泵示功图曲线的上冲程尾部的曲率变化较大，因此右上区突出程度特征向量表示为：

其中，NUM₃表示边界链码类型为z₃的边界链码个数，NUM₄表示边界链码类型为z₄的边界链码个数，NUM₅表示边界链码类型为z₅的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数，为有杆泵抽油井泵示功图的上冲程数据点纵坐标的平均值，为边界链码的类型为z₃、z₄和z₅的边界链码重画后的泵示功图的数据点纵坐标的平均值。

由于系统卸载提前，使泵示功图曲线的上冲程线后半部分曲率变化较小，因此右上区欠缺程度特征向量表示为：

$T_4=\frac{3*{\mathrm{NUM}}_{15}+2*{\mathrm{NUM}}_{16}+{\mathrm{NUM}}_{14}-3*{\mathrm{NUM}}_1-{\mathrm{NUM}}_2}{\mathrm{NUM}}---\left(8\right)$

其中，NUM₁表示边界链码类型为z₁的边界链码个数，NUM₂表示边界链码类型为z₂的边界链码个数，NUM₁₄表示边界链码类型为z₁₄的边界链码个数，NUM₁₅表示边界链码类型为z₁₅的边界链码个数，NUM₁₆表示边界链码类型为z₁₆的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

右上区欠缺程度特征链码示意图如图9所示，其中，实线代表系统正常加载，虚线代表系统卸载变缓。

系统在上冲程结束时快速卸载，使泵示功图曲线的上死点处曲率变化较大，因此右上区快速卸载程度特征向量表示为：

$T_5=\frac{3*{\mathrm{NUM}}_{13}+{\mathrm{NUM}}_{12}-2*{\mathrm{NUM}}_{11}}{\mathrm{NUM}}---\left(9\right)$

其中，NUM₁₁表示边界链码类型为z₁₁的边界链码个数，NUM₁₂表示边界链码类型为z₁₂的边界链码个数，NUM₁₃表示边界链码类型为z₁₃的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

右上区快速卸载程度特征链码示意图如图10所示，其中，实线代表系统正常加载，虚线代表系统快速卸载。

在下冲程末端，由于系统内部相关器件发生碰撞，使载荷突然减小，造成泵示功图曲线的下冲程尾部的曲率变化较大，因此左下区突出程度特征向量表示为：

其中，NUM₁₁表示边界链码类型为z₁₁的边界链码个数，NUM₁₂表示边界链码类型为z₁₂的边界链码个数，NUM₁₃表示边界链码类型为z₁₃的边界链码个数，NUM₁₄表示边界链码类型为z₁₄的边界链码个数，NUM₁₅表示边界链码类型为z₁₅的边界链码个数，NUM₁₆表示边界链码类型为z₁₆的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数，为有杆泵抽油井泵示功图的下冲程数据点纵坐标的平均值，为边界链码的类型为z₁₁、z₁₂、z₁₃、z₁₄、z₁₅和z₁₆的边界链码重画后的泵示功图的数据点纵坐标的平均值。

由于系统的加载提前，使泵示功图曲线的下冲程线后半部分曲率变化较小，因此左下区欠缺程度特征向量表示为：

$T_7=\frac{3*{\mathrm{NUM}}_6+2*{\mathrm{NUM}}_7-{\mathrm{NUM}}_8-2*{\mathrm{NUM}}_9-{\mathrm{NUM}}_{10}}{\mathrm{NUM}}---\left(11\right)$

其中，NUM₆表示边界链码类型为z₆的边界链码个数，NUM₇表示边界链码类型为z₇的边界链码个数，NUM₈表示边界链码类型为z₈的边界链码个数，NUM₉表示边界链码类型为z₉的边界链码个数，NUM₁₀表示边界链码类型为z₉的边界链码个数，NUM₁₀表示边界链码类型为z₁₀的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

左下区欠缺程度特征链码示意图如图11所示，其中，实线代表系统正常加载，虚线代表系统加载提前。

由于系统在卸载时变缓，使泵示功图曲线的卸载线后半部分曲率变化较小，因此右下区欠缺程度特征向量表示为：

$T_8=\frac{2*{\mathrm{NUM}}_9+3*{\mathrm{NUM}}_{10}+{\mathrm{NUM}}_{11}-2*{\mathrm{NUM}}_{12}-3*{\mathrm{NUM}}_{13}}{\mathrm{NUM}}---\left(12\right)$

其中，NUM₁₀表示边界链码类型为z₉的边界链码个数，NUM₁₀表示边界链码类型为z₁₀的边界链码个数，NUM₁₁表示边界链码类型为z₁₁的边界链码个数，NUM₁₂表示边界链码类型为z₁₂的边界链码个数，NUM₁₃表示边界链码类型为z₁₃的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

右下区欠缺程度特征链码示意图如图12所示，其中，实线代表系统正常加载，虚线代表系统加载提前。

由于系统发生急剧卸载，使泵示功图曲线的卸载线的曲率变化较大，因此右下区快速卸载程度特征向量表示为：

$T_9=\frac{3*{\mathrm{NUM}}_9+{\mathrm{NUM}}_{12}+2*{\mathrm{NUM}}_{13}-2*{\mathrm{NUM}}_{11}}{\mathrm{NUM}}---\left(13\right)$

其中，NUM₉表示边界链码类型为z₉的边界链码个数，NUM₁₁表示边界链码类型为z₁₁的边界链码个数，NUM₁₂表示边界链码类型为z₁₂的边界链码个数，NUM₁₃表示边界链码类型为z₁₃的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数。

右下区快速卸载程度特征链码示意图如图13所示，其中，实线代表系统正常加载，虚线代表系统加载提前。

步骤4.3：将各区域的泵示功图特征向量进行量化，得到量化后的泵示功图特征向量。

本发明具体实施方式中得到的一个量化后的泵示功图特征向量如表1所示：

表1 量化后的泵示功图特征向量

步骤5：利用基于粒子群优化的谱聚类算法对量化后的泵示功图特征向量进行聚类，完成有杆泵抽油井故障诊断。

步骤5.1：设将109个泵示功图特征向量聚类为k类，σ为尺度参数，其中，σ∈[0，1]。

步骤5.2：将聚类数k和尺度参数σ作为粒子进行聚类优化，得到最优的聚类数k和最优的尺度参数σ。

步骤5.2.1：将聚类数k和尺度参数σ作为粒子，随机产生若干个粒子，初始化粒子的速度和位置。

步骤5.2.2：109个量化后的泵示功图特征向量由d₁、d₂…d₁₀₉表示，由高斯函数s(d_p，d_q)＝exp(-||d_p-d_q||²/(2σ²))计算任意两个量化后的泵示功图特征向量的相似度，构成相似度矩阵S∈R^109×109，其中，p、q＝1、2…109。

步骤5.2.3：将每个量化后的泵示功图特征向量看作顶点V，根据量化后的泵示功图特征向量间的相似度对各顶点间的边E赋权重值W，得到一个基于量化后的泵示功图特征向量相似度的无向加权图G＝(V，E)。

步骤5.2.4：将相似度矩阵S的每行元素相加，即得到对应顶点的度，以所有顶点的度作为对角元素构成的对角矩阵，即度矩阵D。

步骤5.2.5：根据度矩阵、各顶点间的边的权重及单位矩阵，建立拉普拉斯矩阵L_rw＝I-D^-1W，其中，I为单位矩阵。

步骤5.2.6：计算拉普拉斯矩阵L_rw前k个最小的特征值对应的特征向量u₁、u₂…u_k，k个特征向量组成特征向量矩阵U∈R^109×k。

步骤5.2.7：将特征向量矩阵U∈R^109×k的第g行向量(e_g)_{g＝1、2...109}作为聚类数据集。

步骤5.2.8：用CritC指标函数的倒数作为各个粒子的适应度值。

CritC指标函数的定义如下：

$\mathrm{CritC}={(a\·F)}^{({\log }_2(k+1)+1)}---\left(14\right)$

其中： $F=\left(\frac{1}{1+V/B}\right);$

$B={\mathrm{\Σ}}_{f=1}^k\left|C_f\right|\mathrm{\γ}(C_f,O);$

γ(C_f，O)表示聚类数据集(e_g)_{g＝1、2...109}中每个聚类中心C_f和整个聚类数据集(e_g)_{g＝1、2...109}的中心0之间的距离；

C_f由每个聚类中所有数据的平均值求得；

0由聚类数据集(e_g)_{g＝1、2...109}中所有数据的平均值求得；

|C_f|表示γ(C_f，O)所对应的权值，等于每个聚类中数据的数量；

γ(C_f，o)表示聚类数据集(e_g)_{g＝1、2...109}中每个数据点和它所在聚类的中心C_f之间的距离；

a为常数系数，本实施方式中取a＝0.8；

CritC指标函数在[0，1]区间内取值。

步骤5.2.9：若当前粒子的适应度值小于其历史适应度值，则将当前粒子的适应度值作为其历史适应度值。

步骤5.2.10：设定粒子的适应度值终止条件，若当前粒子的适应度值小于其终止条件，则当前的聚类数k和尺度参数σ为最优的聚类数k和最优的尺度参数σ，当前的聚类数据集(e_g)_{g＝1、2...109}为优化后的聚类数据集(e_g)_{g＝1、2...109}，否则，更新粒子的速度和位置，返回步骤5.2.3。

本实施方式中，设定代表聚类数k和尺度参数σ的粒子的适应度值终止条件为0.5％，聚类数k和尺度参数σ的最优值为5和0.1237。

步骤5.3：利用最优的聚类数k和最优的尺度参数σ进行聚类优化，得到最优的聚类中心。

步骤5.3.1：将优化后的聚类数据集(e_g)_{g＝1、2...109}中数据作为粒子，随机选取k＝5个粒子作为初始聚类中心，初始化粒子的速度和位置。

步骤5.3.2：用Jc指标函数作为粒子的适应度值。

Jc指标函数的定义如下：

$\mathrm{Jc}={\mathrm{\Σ}}_{g=1}^{109}{\mathrm{\Σ}}_{\mathrm{Nc}=1}^5{\left|\right|e_g-{\mathrm{\φ}}_{\mathrm{Nc}}\left|\right|}^2---\left(15\right)$

其中：φ_Nc表示优化后的聚类数据集(e_g)_{g＝1、2...109}的聚类中心，Nc＝1、2、...5。

步骤5.3.3：若当前粒子的适应度值小于其历史适应度值时，则把当前粒子的适应度值作为其历史适应度值，当前粒子作为最优粒子，同时将当前粒子位置记为其历史最优位置。

步骤5.3.4：设定粒子的适应度值终止条件，若当前粒子的适应度值小于其终止条件，则当前最优粒子的聚类中心为最优的聚类中心，否则，更新粒子的速度和位置，返回步骤5.3.2。

步骤5.4：根据最优的聚类中心将109个泵示功图特征向量聚类为5类，得到最终的聚类结果：H₁、H₂…H₅，若待诊断故障类型的泵示功图特征向量与某知故障类型的泵示功图特征向量聚为一类，则两者为同类故障，即得到有杆泵抽油井故障诊断结果，若待诊断故障类型的泵示功图特征向量单独聚为一类，则该待诊断故障类型为新故障类型。

如图14所示，为本发明具体实施方式中对量化后的泵示功图特征向量进行聚类为5类后的示意图。

本实施方式中，待诊断故障类型的量化后的泵示功图特征向量(由“☆”标记)和“正常”类型的量化后的泵示功图特征向量聚类为同一类型，因此该待诊断故障类型为“正常”。

Claims (2)

1.一种有杆泵抽油井故障诊断方法，其特征在于，包括以下步骤：

步骤1：获取M个已知故障类型的有杆泵抽油井地面示功图和1个待诊断故障类型的有杆泵抽油井地面示功图；

步骤2：将M+1个有杆泵抽油井地面示功图转化为有杆泵抽油井泵示功图，并对有杆泵抽油井泵示功图进行归一化处理，得到归一化处理后的泵示功图；

步骤3：采用16方向链码对归一化处理后的泵示功图进行边界链码重画，得到边界链码重画后的泵示功图，其中，边界链码的类型由z_t表示，t＝1、2、…16，NUM_t表示边界链码类型为z_t的边界链码个数，NUM表示边界链码重画后的泵示功图的边界链码总数；

步骤4：将边界链码重画后的泵示功图进行区域划分，并提取各区域的特征向量，各特征向量以16方向链码形式表示；划分的区域包括：左上区、右上区、左下区和右下区；所述的特征向量包括：左上区欠缺程度特征向量、左上区快速加载程度特征向量、右上区突出程度特征向量、右上区欠缺程度特征向量、右上区快速卸载程度特征向量、左下区突出程度特征向量、左下区欠缺程度特征向量、右下区欠缺程度特征向量和右下区快速卸载程度特征向量；

步骤4.1：确定有杆泵抽油井泵示功图的上冲程和下冲程：将有杆泵抽油井泵示功图曲线的横坐标最大的点作为上冲程的结束点，则有杆泵抽油井泵示功图的上冲程为有杆泵抽油井泵示功图曲线第一个数据点到上冲程结束的数据点，有杆泵抽油井泵示功图的下冲程为有杆泵抽油井泵示功图曲线上冲程结束点的下一个数据点到最后一个数据点；

步骤4.2：提取有杆泵抽油井泵示功图各区域的泵示功图特征向量，包括左上区欠缺程度特征向量、左上区快速加载程度特征向量、右上区突出程度特征向量、右上区欠缺程度特征向量、右上区快速卸载程度特征向量、左下区突出程度特征向量、左下区欠缺程度特征向量、右下区欠缺程度特征向量和右下区快速卸载程度特征向量；

左上区欠缺程度特征向量表示为：

$T_1=\frac{{\mathrm{NUM}}_1+3*{\mathrm{NUM}}_2+2*{\mathrm{NUM}}_3+{\mathrm{NUM}}_4-2*{\mathrm{NUM}}_5}{\mathrm{NUM}};$

左上区快速加载程度特征向量表示为：

$T_2=\frac{{\mathrm{NUM}}_3+2*{\mathrm{NUM}}_4+{3*\mathrm{NUM}}_5-{\mathrm{NUM}}_1-2*{\mathrm{NUM}}_2}{\mathrm{NUM}};$

右上区突出程度特征向量表示为：

其中，为有杆泵抽油井泵示功图的上冲程数据点纵坐标的平均值，为边界链码的类型为z₃、z₄和z₅的边界链码重画后的泵示功图的数据点纵坐标的平均值；

右上区欠缺程度特征向量表示为：

$T_4=\frac{3*{\mathrm{NUM}}_{15}+2*{\mathrm{NUM}}_{16}+{\mathrm{NUM}}_{14}-3*{\mathrm{NUM}}_1-{\mathrm{NUM}}_2}{\mathrm{NUM}};$

右上区快速卸载程度特征向量表示为：

$T_5=\frac{3*{\mathrm{NUM}}_{13}+{\mathrm{NUM}}_{12}-2*{\mathrm{NUM}}_{11}}{\mathrm{NUM}};$

左下区突出程度特征向量表示为：

其中，为有杆泵抽油井泵示功图的下冲程数据点纵坐标的平均值，为边界链码的类型为z₁₁、z₁₂、z₁₃、z₁₄、z₁₅和z₁₆的边界链码重画后的泵示功图的数据点纵坐标的平均值；

左下区欠缺程度特征向量表示为：

$T_7=\frac{3*{\mathrm{NUM}}_6+2*{\mathrm{NUM}}_7-{\mathrm{NUM}}_8-2*{\mathrm{NUM}}_9-{\mathrm{NUM}}_{10}}{\mathrm{NUM}};$

右下区欠缺程度特征向量表示为：

$T_8=\frac{2*{\mathrm{NUM}}_9+3*{\mathrm{NUM}}_{10}+{\mathrm{NUM}}_{11}-2*{\mathrm{NUM}}_{12}-3*{\mathrm{NUM}}_{13}}{\mathrm{NUM}};$

右下区快速卸载程度特征向量表示为：

$T_9=\frac{3*{\mathrm{NUM}}_9+{\mathrm{NUM}}_{12}+2*{\mathrm{NUM}}_{13}-2*{\mathrm{NUM}}_{11}}{\mathrm{NUM}};$

步骤4.3：将各区域的泵示功图特征向量进行量化，得到量化后的泵示功图特征向量；

步骤5：利用基于粒子群优化的谱聚类算法对量化后的泵示功图特征向量进行聚类，完成有杆泵抽油井故障诊断；

步骤5.1：设将M+1个泵示功图特征向量聚类为k类，σ为尺度参数，其中，σ∈[0，1]；

步骤5.2：将聚类数k和尺度参数σ作为粒子进行聚类优化，得到最优的聚类数k和最优的尺度参数σ；

步骤5.3：利用最优的聚类数k和最优的尺度参数σ进行聚类优化，得到最优的聚类中心；

步骤5.4：根据最优的聚类中心将M+1个泵示功图特征向量聚类为k类，得到最终的聚类结果：H₁、H₂…H_k，若待诊断故障类型的泵示功图特征向量与某知故障类型的泵示功图特征向量聚为一类，则两者为同类故障，即得到有杆泵抽油井故障诊断结果，若待诊断故障类型的泵示功图特征向量单独聚为一类，则该待诊断故障类型为新故障类型，对该故障类型的泵示功图特征向量进行保存。

2.根据权利要求1所述的一种有杆泵抽油井故障诊断方法，其特征在于，所述的步骤5.2按以下步骤执行：

步骤5.2.1：将聚类数k和尺度参数σ作为粒子，随机产生若干个粒子，初始化粒子的速度和位置；

步骤5.2.2：M+1个量化后的泵示功图特征向量由d₁、d₂…d_M+1表示，由高斯函数s(d_p，d_q)＝exp(-||d_p-d_q||²/(2σ²))计算任意两个量化后的泵示功图特征向量的相似度，各相似度值构成相似度矩阵S∈R^(M+1)×(M+1)，其中，p、q＝1、2…M+1；

步骤5.2.3：将每个量化后的泵示功图特征向量看作顶点V，根据量化后的泵示功图特征向量间的相似度对各顶点间的边E赋权重值W，得到一个基于量化后的泵示功图特征向量相似度的无向加权图G＝(V，E)；

步骤5.2.4：将相似度矩阵S的每行元素相加，即得到对应顶点的度，以所有顶点的度作为对角元素构成的对角矩阵，即度矩阵D；

步骤5.2.5：根据度矩阵、各顶点间的边的权重及单位矩阵，建立拉普拉斯矩阵L_rw＝I-D^-1W，其中，I为单位矩阵；

步骤5.2.6：计算拉普拉斯矩阵L_rw前k个最小的特征值对应的特征向量，k个特征向量组成特征向量矩阵U；

步骤5.2.7：将特征向量矩阵U的第g行向量(e_g)_{g＝1，2，…，M+1}作为聚类数据集；

步骤5.2.8：用CritC指标函数的倒数作为各个粒子的适应度值；

步骤5.2.9：若当前粒子的适应度值小于其历史适应度值，则将当前粒子的适应度值作为其历史适应度值；

步骤5.2.10：设定粒子的适应度值终止条件，若当前粒子的适应度值小于其终止条件，则当前的聚类数k和尺度参数σ为最优的聚类数k和最优的尺度参数σ，当前的聚类数据集(e_g)_{g＝1，2，…，M+1}为优化后的聚类数据集(e_g)_{g＝1，2，…，M+1}，否则，更新粒子的速度和位置，返回步骤5.2.3。