2019,
Vol. 43
2. 西安理工大学经济与管理学院, 陕西西安 710054
2. School of Economics and Management, Xi'an University of Technology, Xi'an 710054, China
泵功图包含了抽油泵工况的丰富信息。不同井下问题都会通过泵功图形状特征体现, 因此可以根据泵功图来诊断井下设备的故障。诊断过程包含特征提取和特征匹配两部分; 目前常见的泵功图特征提取方法有两种:灰度特征和几何特征[1-3]。常见的泵功图匹配技术有神经网络法、支持向量机法等。通过采用泵功图特征提取方法与匹配技术的结合可诊断出抽油泵的故障情况[4-6]。起初在泵功图的故障诊断技术中多数学者都利用灰度理论对泵功图进行特征提取[7-9]。随后学者又对泵功图的几何特征进行提取。王科科[10]提取不变矩和泵功图形状特征, 使用BP神经网络诊断抽油泵的工况。姜民政等[11]提取泵功图的不变矩特征, 应用SOM神经网络识别故障类型。Li等[12-13]提取泵功图 28个不变曲线矩特征向量, 先后采用粒子群算法和改进的果蝇优化算法诊断抽油系统的故障。形状上下文(shape context)算法[14-16]能够描述轮廓上每个点与轮廓上其余点的空间分布关系, 有很强的形状描述能力, 被广泛应用于各个领域, 例如遥感图像的匹配[17], 道路标识牌的检测[18]以及图像拷贝的检测[19]等。因此, 笔者选择形状上下文算法描述泵功图的形状, 计算待诊断泵功图和故障泵功图之间任意两点的匹配代价值, 获得代价矩阵; 利用回溯法对代价矩阵进行搜索, 计算匹配代价值, 诊断抽油泵的故障类型。
1 泵功图获取鉴于实测泵功图的获取难度, 通常对有杆抽油系统建立诊断模型。通过对模型的求解, 将抽油机悬点示功图转化为泵功图。现有的诊断模型大多是根据吉布斯方程[20]演变而来, 本文中选取多级杆柱运动的波动方程[21], 表示为
| $ \frac{{{\partial ^2}u}}{{\partial {t^2}}} = a_k^2\frac{{{\partial ^2}u}}{{\partial {x^2}}} - {c_k}\frac{{\partial u}}{{\partial t}} - \delta {h_k}{m_{(x)}} + g_k^\prime \cos \alpha ,k = 1,2, \cdots ,K. $ |
其中
| $ \left\{ {\begin{array}{*{20}{l}} {{a_k} = \sqrt {{E_{{\rm{r}}k}}/{\rho _{{\rm{r}}k}},} }\\ {{c_k} = {v_{\rm{e}}}/{\rho _{{\rm{r}}k}}{A_{{\rm{r}}k}},}\\ {{h_k} = {\mu _k}/{\rho _{{\rm{r}}k}}{A_{{\rm{r}}k}},}\\ {g_k^\prime = \rho _{{\rm{r}}k}^\prime g/{\rho _{rk}}.} \end{array}} \right. $ | (1) |
式中, u为抽油杆柱上任意位置x处的截面在任意t时刻沿井深方向的位移, m; Erk为抽油杆柱的弹性模量, Pa; ρrk为抽油杆柱在空气中的密度, kg/m3; m为单位长度上杆柱的黏滞阻尼系数, kg/(m·s), Ark为抽油杆柱横截面积, m2; δ为符号函数, 上、下冲程时分别取-1和+1;μk为第k级抽油杆柱与油管间的摩擦系数; N(x)为油管对微元体的支反力, N; ρ′rk为微元体在油液中的相对密度, kg/m3; K为多级杆级数。
将抽油机的悬点位移和悬点载荷作为波动方程的一个边界条件[21], 利用有限差分法对方程(1)求解, 可获得抽油杆柱第i个节点处在时刻j的位移uij。
抽油泵处柱塞载荷对应于抽油杆柱的末端。根据牛顿后插公式, 可获得抽油泵载荷计算公式[2]为
| $ {F_{N,j}} = {E_{\rm{r}}}{A_{\rm{r}}}{\left( {\frac{{\partial u}}{{\partial x}}} \right)_{N,j}} = \frac{{{E_{\rm{r}}}{A_{\rm{r}}}}}{{\Delta x}}\left( {\frac{3}{2}{u_{N,j}} - 2{u_{N - 1,j}} + \frac{1}{2}{u_{N - 2,j}}} \right). $ | (2) |
式中, Er为抽油杆柱末端的弹性模量, MPa; Ar为抽油杆柱末端的横截面积, m2; Δx为抽油杆柱的离散步长, m。
将获得的抽油杆柱第N-2、N-1和N节点处在时刻j的位移uN-2, j、uN-1, j和uN, j带入式(2)中, 即可获得抽油泵在时刻j受到的载荷FN, j, 进而绘制泵功图。
2 基于形状上下文匹配方法为了计算两个待匹配泵功图的相似程度, 采用基于形状上下文的匹配方法。首先选取轮廓上有规律的不连续的样本点表示轮廓。其次利用形状上下文描述符对轮廓上的点进行描述, 将二维图像转化为矢量集。然后, 计算两个待匹配的形状轮廓点之间的匹配代价值, 获得一个代价矩阵。最后, 采用回溯法在代价矩阵中搜索到最优匹配点集, 计算两个待匹配泵功图的匹配代价值。
2.1 形状上下文算法形状上下文是一种基于轮廓点集空间分布的描述符。利用形状上下文描述符进行形状描述, 可以获得轮廓上点与点之间的距离以及点在形状中的相对位置关系。基本思想是对于N个点的形状轮廓, 任意一点与其余N-1个点的矢量集表示该点相对于整个形状的信息。由此, 一个二维图像可以表示成N个矢量集。考虑到位置分布的二维不变性, 将笛卡尔坐标系上的点进行对数极坐标变换, 以对数极坐标的直方图表示点的位置分布。
形状I(x, y)到I(r, θ)的极坐标变换定义为
| $ r = \sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {y - {y_0}} \right)}^2}} , $ | (3) |
| $ \theta = \arctan \left( {\frac{{y - {y_0}}}{{x - {x_0}}}} \right). $ | (4) |
一般情况下(x0, y0)为坐标原点。
用复数表示为Z=x+iy=r(cos θ+isin θ)=reiθ。令ω=lnZ=lnr+iθ, 则笛卡尔坐标转换为对数极坐标的映射方程为p(r, θ)=lnr, q(r, θ)=θ。因此, 当笛卡尔坐标中的图像相对于坐标原点发生缩放和旋转变化时, 即图像放大r0倍, 旋转θ0度, 相应的极坐标为(r+r0, θ+θ0), 对数极坐标空间ω(μ, v)的图像只产生平移, 即μ=lnr0+lnr, v=θ0+θ。这就是所谓的对数极坐标映射的二维不变性。
给定一个轮廓上的点集P={p1, p2, …, pN}, 计算任意点pi(i=1, 2, …, N)的对数极坐标直方图。以pi为坐标原点进行对数极坐标变换。以点pi为极坐标原点, 将整个二维平面分为4×n个方向(n为正整数), 按照lnr划分为3份, 整个平面被分为12n份, 每一份定义为一“格”; 统计点集P落在12n个格中的个数。以12点方向为基准, 顺时针定义4n个方向为直方图的列, 每个方向沿原点向外的每个格定义为直方图的行。根据统计的个数获得点pi的直方图。这个直方图为pi点的形状上下文。
图 1给出了形状上下文的描述过程。其中图(a)为形状一的轮廓取样点, (b)为形状二的轮廓取样点, (c)、(d)和(e)分别是以(a)中三角点、正方形点和(b)中圆形点为对数极坐标原点的变换图, (f)、(g)和(h)为n取2时这3个点对应的形状直方图。
|
图 1 形状上下文描述过程 Fig.1 Description process of shape contexts |
为了描述两个形状轮廓中任意两点的匹配情况, 定义了两点的匹配代价。设pi为第一幅图轮廓上的一点, qj为第二幅图轮廓上的一点, pi和qj的匹配代价为Cij=C(pi, qj), 运用χ2分布, 有
| $ C\left( {{p_i},{q_j}} \right) = \frac{1}{2}\sum\limits_{k = 1}^K {\frac{{{{\left[ {{h_i}(k) - {h_j}(k)} \right]}^2}}}{{{h_i}(k) + {h_j}(k)}}} . $ | (5) |
式中, hi(k)为pi直方图中第k格的数值; hj(k)为qj直方图中第k格的数值。
根据方程(5)便可计算两幅泵功图轮廓点集中每两个点之间的匹配代价值。
2.2 形状匹配对于两个待匹配的泵功图, 设实况泵功图提取的轮廓点数为N, 故障泵功图提取的轮廓点数为M, 根据这两个形状点集之间的匹配代价函数生成一个代价矩阵(以轮廓点数多的作为横坐标, 轮廓点数少的为纵坐标)。采用搜索法对两幅泵功图的轮廓点进行全局匹配, 使得代价矩阵中第一行第一列元素和最后一行最后一列元素被选择, 获得两幅泵功图的匹配代价值。
考虑到回溯法是以深度优先方式搜索解空间, 并且在搜索过程中能够避免无效搜索, 因此采用回溯法搜索代价矩阵。具体步骤如下:
(1) 从开始结点(根结点)出发, 向纵深方向形成搜索子空间并移至一个新结点(活结点);
(2) 从新结点(活结点)出发, 如果新结点可以在搜索子空间中移动, 那么新结点为活结点, 找到下一个新结点; 否则新结点为死结点, 此时应往回移动(回溯)至最近一个活结点重新搜索;
(3) 不断重复步骤(2)过程, 直至搜索到代价矩阵的边界, 此时搜索结束。
通过回溯法搜索代价矩阵, 最终实现两点集终点到终点、起点到起点的最优对齐, 满足匹配点数最多、空格和不匹配点数最少的要求。定义搜索目标为总代价值最小, 表示为
| $ H(\pi ) = \min \sum\limits_{i \in N} C \left( {{p_i},{q_{\pi (i)}}} \right). $ | (6) |
式中, pi为实况泵功图轮廓点; qπ(i)为与其对齐的故障泵功图轮廓点。
泵功图中的轮廓点是从图中左下角抽油机的下死点处(图 2中的*点)按顺时针顺序形成的, 因此回溯开始在代价矩阵的最后一行最后一列处, 采用当前位置的相邻空间作为搜索子空间。当代价矩阵的行列数相差不大时, 搜索子空间就相对较小, 可以从开始点处以左、上、对角位置为回溯点, 即2×2的子空间; 当代价矩阵的行列数相差较大时, 搜索子空间也相对增大, 从2×2的搜索空间开始逐渐增加空间尺寸, 来尽可能的实现全局最优对齐。然后, 利用回溯法获得代价矩阵的回溯路径, 即两个形状点集的对应关系。最后, 定义对齐的两形状轮廓点的总代价值除以实况泵功图轮廓点数为两形状的匹配代价值, 用F(N, M)表示, 且该值越小表示两个形状越相似, 其中
|
图 2 轮廓点生成顺序 Fig.2 Formation order of contour points |
| $ F(N,M) = H({\rm{ \mathsf{ π} }})/N. $ | (7) |
例如, 两个形状轮廓上的点集分别是P={p1, p2, …, p250}, Q={q1, q2, …, q145}, 利用形状上下文算法对两个形状轮廓进行描述, 并根据任意两个点的匹配代价值生成一个250×145的代价矩阵, 图 3为代价矩阵的一部分。由于代价矩阵的行数和列数相差较多, 在回溯法中采用4×3矩阵的搜索子空间能够实现两个形状的全局最优对齐。图 3中的虚线框为搜索子空间, 实线框为匹配点对, 同时实线框构成了回溯路径。当回溯路径中某一点对应多个点时, 选择其中匹配代价值最小的为匹配点。当代价矩阵中某行或某列中未出现回溯点, 则说明该点不匹配。如果待诊断的泵功图轮廓上出现不匹配的点, 则给该点设置一个较大的匹配代价值(代价矩阵中最大的代价值)。
|
图 3 4×3搜索空间下的回溯路径 Fig.3 Trace-back path in 4×3 search space |
为了验证所提出的基于形状上下文匹配方法的有效性, 以文献[7]中给出的油井1#和7#的实测悬点数据以及辽河油田两口油井作为验证算例。图 4为根据实测悬点载荷绘制的待诊断4口油井悬点示功图。
|
图 4 待诊断油井悬点示功图 Fig.4 Polished rod dynamometer cards of oil wells to be diagnosed |
抽油泵常见的故障类型有气体影响、供液不足、游动阀漏失、固定阀漏失等。以这4种故障泵功图和正常泵功图作为样板图, 分别与4口油井的实况泵功图进行匹配, 诊断这4口井的故障。将已知油井的实测悬点位移和悬点载荷带入方程(1)和(2), 将悬点示功图转化成泵功图; 提取两幅泵功图轮廓上的离散点集, 利用形状上下文算法获得两点集中每个点的上下文信息; 计算两点集中任意两个点之间的匹配代价值, 得到代价矩阵; 使用回溯法搜索代价矩阵, 找到最优的点集对应关系; 计算两幅泵功图的匹配代价值, 判断故障类型。具体流程见图 5。
|
图 5 基于形状上下文匹配方法流程 Fig.5 Flow chart based on shape context matching method |
考虑到匹配过程中产生的数据较多, 只给出油井1#的匹配数据。由实测悬点的位移和载荷, 获得油井1#的泵功图, 如图 6所示。提取泵功图上的离散点, 描述每个点的形状上下文信息。任取泵功图上*和◇两个点为例, 绘制形状上下文直方图:根据形状上下文算法, 分别以这两点作为极坐标原点, 将整个二维平面分为12个方向, 获得这两点的形状直方图, 如图 7所示。根据油井1#泵功图和正常泵功图的上下文信息, 计算这两幅图中任意两个点之间的匹配代价值, 获得这两幅图的代价矩阵。采用回溯法找到代价矩阵中的回溯路径, 见图 8。结合方程(7), 获得两幅泵功图的匹配代价值。
|
图 6 油井1#泵功图 Fig.6 Pumping dynamometer card of well 1# |
|
图 7 两点的形状直方图 Fig.7 Shape histogram of two points |
|
图 8 油井1#泵功图和正常泵功图匹配代价矩阵中的回溯路径 Fig.8 Trace-back path in cost matrix of pumping dynamometer cards between well 1# and normal well |
为了验证本文方法的有效性, 表 1给出了基于灰度理论的诊断结果以及提出的基于形状上下文匹配方法的诊断结果。便于和灰度关联值进行比较, 计算的匹配代价值均已归一化处理。归一化过程为:选取同一油井实测泵功图与样板图的匹配代价最大值并向上取整, 然后用每一个匹配代价值与该整数相除。基于灰度理论方法是以灰度关联值最大为诊断结果, 基于形状上下文方法是以匹配代价值最小为诊断结果。由表 1可知, 本文方法的诊断结果和实际情况相同, 这说明提出的匹配方法是有效的。
|
表 1 基于灰度理论和形状上下文的两种匹配方法的诊断结果 Table 1 Diagnostic results of two matching methods based on gray theory and shape context |
为了说明匹配方法诊断结果的可靠性, 定义了最小对比度:计算结果中最优和次优两个数值差值的绝对值除以二者中较大的数值。最小对比度越小, 说明最优值和次优值越接近, 诊断结果越易受到计算误差等因素的影响。观察表 1可知, 在灰度理论方法中油井1#的泵功图与气体影响的灰度关联值为0.993 9, 与泵正常的灰度关联值为0.990 3、与游动凡尔漏失的灰度关联值为0.989 7, 虽然诊断结果为气体影响, 但是油井1#诊断结果的灰度关联值与其他故障的最小对比度仅为3.6×10-3。而提出的基于形状上下文匹配方法获得的匹配代价值与其他故障的最小对比度为0.547 7。与此同时, 对于油井7#、辽河4#和5#的泵功图, 基于灰度理论诊断结果的最小对比度分别为5×10-4、7.9×10-3和3.7×10-3, 而本文方法诊断结果的最小对比度分别为0.049、0.215 6和0.022 2, 明显大于灰度理论的最小对比度。因此本文提出的匹配方法提高了抽油泵故障诊断的可靠性。
4 结论(1) 针对基于灰度理论的图像匹配技术存在的匹配可靠性问题, 根据泵功图的几何特征, 提出的基于形状上下文的匹配方法能准确地描述细节较多的泵功图, 从形状描述符的选取上保证了匹配结果的准确度。
(2) 新匹配方法以两幅图的匹配代价值表示它们的相似程度, 且匹配代价值越小两幅图越相似, 能够准确诊断出实测泵功图的故障类型。给出的4个算例的诊断结果与实际情况相同验证了本文提出的方法的有效性。
(3) 基于形状上下文的匹配方法的最小对比度明显高于基于灰度理论的匹配方法, 提高了有杆抽油系统故障诊断结果的可靠性, 对细节较多的泵功图诊断问题提供了更准确的匹配方法, 可以广泛应用于油井的故障诊断作业中。
| [1] |
刘炜, 刘宏昭. 三种泵功图特征提取方法及比较[J]. 煤炭技术, 2010, 29(9): 202-204. LIU Wei, LIU Hongzhao. Three methods of feature extraction of pumping cards and their comparison[J]. Coal Technology, 2010, 29(9): 202-204. |
| [2] |
李兆敏, 林日亿, 付路长, 等. 有杆抽油系统效率分析及抽汲参数优化设计[J]. 石油学报, 2005, 26(5): 106-110. LI Zhaomin, LIN Riyi, FU Luchang, et al. Efficiency analysis of sucker-rod pumping system and optimal design of swabbing parameter[J]. Acta Petrolei Sinica, 2005, 26(5): 106-110. |
| [3] |
ZHENG B Y, GAO X W. Sucker rod pumping diagnosis using valve working position and parameter optimal continuous hidden Markov model[J]. Journal of Process Control, 2017, 59: 1-12. DOI:10.1016/j.jprocont.2017.09.007 |
| [4] |
LI K, GAO X W, ZHOU H B, et al. Fault diagnosis for down-hole conditions of sucker rod pumping systems based on the FBH-SC method[J]. Petroleum Science, 2015, 12(1): 135-147. DOI:10.1007/s12182-014-0006-5 |
| [5] |
GAO Q, SUN S B, LIU J C. Fault diagnosis of suck rod pumping system via extreme learning machines: IEEE International Conference on Cyber Technology in Automation, Control, and Intelligent Systems (CYBER), 2015[C]. Shenyang: IEEE, 2015: 503-507.
|
| [6] |
LI K, GAO X W, YANG W B, et al. Multiple fault diagnosis of down-hole conditions of sucker-rod pumping wells based on Freeman chain code and DCA[J]. Petroleum Science, 2013, 10(3): 347-360. DOI:10.1007/s12182-013-0283-4 |
| [7] |
许广繁, 杨善林. 灰色理论在有杆抽油井故障诊断中的应用[J]. 合肥工业大学学报(自然科学版), 2013, 36(10): 1265-1268. XU Guangfan, YANG Shanlin. Application of grey theory in fault diagnosis of rod-pumped well[J]. Journal of Hefei University of Technology(Natural Science), 2013, 36(10): 1265-1268. DOI:10.3969/j.issn.1003-5060.2013.10.023 |
| [8] |
管闯.基于灰度轮廓势能与支持向量机的抽油井故障诊断[D].大庆: 东北石油大学, 2014. GUAN Chuang. The fault diagnosis of oil pumping based on grey profile potential energy and SVM[D]. Daqing: Northeast Petroleum University, 2014. |
| [9] |
LI K, HAN Y, WANG T. A novel prediction method for down-hole working conditions of the beam pumping unit based on 8-directions chain codes and online sequential extreme learning machine[J]. Journal of Petroleum Science and Engineering, 2018, 160: 285-301. DOI:10.1016/j.petrol.2017.10.052 |
| [10] |
王科科.远程监测抽油机井工况智能诊断技术[D].青岛: 中国石油大学(华东), 2009. WANG Keke, Intelligent diagnosis technology for remote monitoring pumping well conditions[J]. Qingdao: China University of Petroleum(East China), 2009. |
| [11] |
姜民政, 段天玉, 张迪, 等. 基于RS-LVQ的同井注采系统故障诊断研究[J]. 石油机械, 2018, 46(3): 95-99. JIANG Minzheng, DUAN Tianyu, ZHANG Di, et al. RS-LVQ-based fault diagnosis of the injection-production system[J]. China Petroleum Machinery, 2018, 46(3): 95-99. |
| [12] |
LI K, GAO X W, TIAN Z D, et al. Using the curve moment and the PSO-SVM method to diagnose downhole conditions of a sucker rod pumping unit[J]. Petroleum Science, 2013, 10(1): 73-80. DOI:10.1007/s12182-013-0252-y |
| [13] |
李琨, 韩莹, 佘东生, 等. 基于IFOA-KELM-MEA模型的游梁式抽油机采油系统井下工况的短期预测[J]. 化工学报, 2017, 68(1): 188-198. LI Kun, HAN Ying, SHE Dongsheng, et al. IFOA-KELM-MEA model based transient prediction on down-hole working conditions of beam pumping units[J]. CIESC Journal, 2017, 68(1): 188-198. |
| [14] |
刘屿, 孙坤, 谢宏威, 等. 基于形状上下文特征和ICP的高精度轮廓视觉检测算法[J]. 华南理工大学学报(自然科学版), 2018, 46(1): 131-138, 144. LIU Yu, SUN Kun, XIE Hongwei, et al. High accurate contour vision algorithm based on shape context and ICP[J]. Journal of South China University of Technology (Natural Science Edition), 2018, 46(1): 131-138, 144. DOI:10.3969/j.issn.1000-565X.2018.01.017 |
| [15] |
LIM M J, HAN C H, LEE S W, et al. Fast shape matching using statistical features of shape contexts[J]. IEICE TRANSACTIONS on Information and Systems, 2011, 94(10): 2056-2058. |
| [16] |
QIAN D H, CHEN T S, QIAO H. Background of sha-pe contexts for point matching[J]. Pattern Re-cognition Letters, 2016, 84: 114-119. DOI:10.1016/j.patrec.2016.07.016 |
| [17] |
JIANG J, ZHANG S M, CAO S X. Rotation and scale invariant shape context registration for remote sensing images with background variations[J]. Journal of Applied Remote Sensing, 2015, 9(1): 095092. DOI:10.1117/1.JRS.9.095092 |
| [18] |
LIU C S, CHANG F L, LIU C Y. Multiview road sign detection via self-adaptive color model and shape context matching[J]. Journal of Electronic Imaging, 2016, 25(5): 051202. DOI:10.1117/1.JEI.25.5.051202 |
| [19] |
ZHOU Z L, WANG Y L, WU Q M J, et al. Effective and efficient global context verification for image copy detection[J]. IEEE Transactions on Information Forensics and Security, 2017, 12(1): 48-63. DOI:10.1109/TIFS.2016.2601065 |
| [20] |
GIBBS S G. Predicting the behavior of sucker-rod pumping systems[J]. Journal of Petroleum Technology, 1963, 15(7): 769-778. DOI:10.2118/588-PA |
| [21] |
王凯. 定向井有杆抽油泵系统诊断模型的建立[J]. 中国石油大学学报(自然科学版), 2010, 34(2): 130-135. WANG Kai. Modeling method for fault diagnosis to sucker rod pumping system in directional well[J]. Journal of China University of Petroleum(Edition of Natural Science), 2010, 34(2): 130-135. DOI:10.3969/j.issn.1673-5005.2010.02.026 |