CN106321072A - 一种基于泵功图的抽油井故障诊断方法 - Google Patents

Abstract

本发明涉及一种基于泵功图的抽油井故障诊断方法，采集地面示功图，并将示功图转化为泵功图；将得到的泵功图进行故障诊断。本发明在深入理解有杆抽油系统的工作原理基础上设计，有很好的全面性与实用性，能有效适应不同井况，使之高效准确的利用WIA-PA无线示功仪采集的地面示功图进行有杆抽油系统的故障诊断。

Description

一种基于泵功 图 的抽油井故障诊断方法

技术领域

本发明涉及石油化工以及无线传感器领域，具体地说是一种将地面示功图转化为泵功图并进行故障诊断的方法。

背景技术

石油作为一种不可再生资源，对于国家的政治、经济、军事有着不可替代的战略意义。随着我国经济持续、稳定、快速的发展，石油消费也持续增长，原油需求持续上升，石油供需缺口有增大的趋势，并逐渐成为制约我国经济发展的重要因素之一。截止2009年底，中国石油探明剩余可采储量为27.9亿吨，剩余可采储量储采比为14.8，然而在这些剩余可采储量中，石油资源储备质量较差，低渗、特低渗或超低渗油、稠油和埋深大于3500m的石油资源超过50％，不仅开发勘探难度在逐渐加大，开发成本增加而且老油区综合含水高，普遍进入产量递减阶段。面对这样的严峻形势，各油田企业都意识到了提高原油生产效率的重要性，并在降低投资、节约成本、提高油田数字化管理水平等方面做了大量的努力。其中，将WIA技术应用于油田生产数字化管理是一个重要体现。

WIA-PA(Wireless Networks for Industrial Automation Process Automation)技术是中国科学院沈阳自动化研究所自主开发的面向工业过程自动化的工业无线网络标准，是中国工业无线联盟针对过程自动化领域制定的WIA子标准，是基于IEEE 802.15.4标准的用于工业过程测量、监视与控制的无线网络系统。WIA-PA网络协议遵循ISO/OSI的七层结构，但只定义了数据链路层、网络层、应用层。WIA-PA网络由主控计算机、网关设备、路由设备、现场设备和手持设备5类物理设备构成，此外，还定义了两类逻辑设备：网络管理器、安全管理器。WIA-PA网络采用星形和网状结合的两层网络拓扑结构，第一层的网络结构由网关和路由设备组成，第二层的星型结构由路由设备及现场设备或手持设备组成。

油田生产由油井、计量站、管汇阀组、转油站、联合站、原油外输系统以及其他分散设施组成。由于老油田产量低，生产成本高，其投入产出比低，目前仍然采用人工每日定时检查设备运行情况并测量、统计采油数据。随着油田压缩编制及夜班无人值守，采油现场发生意外故障不能及时发现和排查，直接影响着油田的生产管理和安全。基于以上问题，将基于WIA-PA技术的远程无线监控及上位机优化软件应用于实际，是很有意义的。WIA技术通信网络覆盖面广(两点间通信无障碍2KM)、可靠性高、耗电量低、基本无运行费用，并且具有远程无线传输和拓展功能。

目前关于有杆泵抽油系统井下工况计算机智能诊断方法的研究已经取得了一些成果，例如:专家系统、人工神经网络、粗糙集理论、支持向量机、频谱分析和滤波技术。但是，由于井下工作状况的不可见性，很难得到足够多的学习样本，很难从已有的信息精确地判断待诊断样本。

蔡文等首创的可拓学，采用形式化设计从定性和定量两个角度同时考虑解决问题的规律和方法，利用信息的可拓展性去解决问题。所建立的物元模型，用一个三元组来描述事物，其中训练样本特征向量的特征值由一个点拓展到一个区间，通过计算待诊断点与特征区间的距离来判断其与每个特征向量的关度。通过关联函数，将点与区间的位置关系用定量的形式刻画出来，用以描述关联度的区别。

发明内容

针对现有技术的不足，本发明提供一种基于WIA-PA无线示功仪采集的地面示功图，利用吉布斯求解波动方程得到泵功图，利用四点法和曲线矩理论进行基于泵功图的有杆抽油系统故障诊断的方法。

本发明为实现上述目的所采用的技术方案是：

一种基于泵功 图 的抽油井故障诊断方法，采集地面示功图，并将示功图转 化为泵功图；将得到的泵功图进行故障诊断。

所述将示功图转化为泵功图包括以下步骤：

步骤1：利用有杆抽油系统动态参数预测模型，预测悬点动载荷的变化：

$\frac{{\∂}^{2}U(x,t)}{\∂t^2}=a^2\frac{{\∂}^{2}U(x,t)}{\∂x^2}-C\frac{\∂U(x,t)}{\∂t}$

其中，a为声波在抽油杆中的速度；C为吉布斯粘滞阻尼系数；

步骤2：采用傅里叶变换的方法，求出抽油杆上任意深度x、任意时间t的载荷与位移，得到泵功图：

$U(x,t)=\frac{{\mathrm{\σ}}_0}{2E{A}_r}x+\frac{v_0}{2}+\underset{n=1}{\overset{n}{\mathrm{\Σ}}}(O_n\left(x\right)\cos \mathrm{n\ωt}+P_n\left(x\right)\sin \mathrm{n\ωt})$

$F(x,t)=\frac{{\mathrm{\δ}}_0}{2}x+{\mathrm{EA}}_r\underset{n=1}{\overset{n}{\mathrm{\Σ}}}(O_n^{\′}\left(x\right)\cos \mathrm{n\ωt}+P_n^{\′}\left(x\right)\sin \mathrm{n\ωt})$

其中，U(x,t)为抽油杆的位移；F(x,t)为抽油杆的载荷；σ₀,δ₀为傅氏系数；E为抽油杆的弹性模量；A_r为抽油杆截面积；O_n,O'_n,P_n,P′_n为特殊函数；n为傅氏级数；ω为曲柄角速度；n为傅氏级数。

所述特殊函数为：

O_n(x)＝(k_nchβ_nx+δ_nshβ_nx)sinα_nx+(μ_nchβ_nx+v_nshβ_nx)cosα_nx

P_n(x)＝(k_nchβ_nx+δ_nchβ_nx)cosα_nx+(μ_nchβ_nx+v_nshβ_nx)sinα_nx

$\begin{array}{c}{O}_n^{\′}\left(x\right)=[\frac{{\mathrm{\τ}}_n}{{\mathrm{EA}}_r}\mathrm{sh}{\mathrm{\β}}_{n}x+({\mathrm{\δ}}_n{\mathrm{\β}}_n-v_n{\mathrm{\α}}_n)\mathrm{ch}{\mathrm{\β}}_{n}x]\sin {\mathrm{\α}}_{n}x\\ +[\frac{{\mathrm{\δ}}_n}{{\mathrm{EA}}_r}\mathrm{ch}{\mathrm{\β}}_{n}x+(v_n{\mathrm{\β}}_n-{\mathrm{\δ}}_n{\mathrm{\α}}_n)\mathrm{sh}{\mathrm{\β}}_{n}x]\cos {\mathrm{\α}}_{n}x\end{array}$

$\begin{array}{c}{P}_n^{\′}\left(x\right)=[\frac{{\mathrm{\τ}}_n}{{\mathrm{EA}}_r}\mathrm{ch}{\mathrm{\β}}_{n}x+({\mathrm{\δ}}_n{\mathrm{\β}}_n-v_n{\mathrm{\α}}_n)\mathrm{sh}{\mathrm{\β}}_{n}x]\cos {\mathrm{\α}}_{n}x\\ +[\frac{{\mathrm{\σ}}_n}{{\mathrm{EA}}_r}\mathrm{sh}{\mathrm{\β}}_{n}x+(v_n{\mathrm{\β}}_n-{\mathrm{\δ}}_n{\mathrm{\α}}_n)\mathrm{ch}{\mathrm{\β}}_{n}x]\sin {\mathrm{\α}}_{n}x\end{array}$

$k_n=\frac{{\mathrm{\σ}}_n{\mathrm{\α}}_n+{\mathrm{\τ}}_n{\mathrm{\β}}_n}{{\mathrm{EA}}_r({\mathrm{\α}}_n^2+{\mathrm{\β}}_n^2)},{\mathrm{\μ}}_n=\frac{{\mathrm{\σ}}_n{\mathrm{\β}}_n-{\mathrm{\τ}}_n{\mathrm{\α}}_n}{{\mathrm{EA}}_r({\mathrm{\α}}_n^2+{\mathrm{\β}}_n^2)}$

${\mathrm{\α}}_n=\frac{\mathrm{n\ω}}{a\sqrt{2}}\sqrt{1+\sqrt{1+{\left(\frac{C}{\mathrm{n\ω}}\right)}^2}},{\mathrm{\β}}_n=\frac{\mathrm{n\ω}}{a\sqrt{2}}\sqrt{-1+\sqrt{1+{\left(\frac{C}{\mathrm{n\ω}}\right)}^2}}$

其中，O_n,O'_n,P_n,P′_n为特殊函数；σ_n,δ_n,τ_n,v_n为傅氏系数；E为抽油杆的弹性模量；A_r为抽油杆截面积；α_n,β_n,k_n,μ_n,a_n,b_n,c_n,d_n均为特殊常数；ω为曲柄角速度；a为声波在抽油杆中的速度；C为吉布斯粘滞阻尼系数；n为傅氏级数；

所述傅氏系数为：

$\begin{array}{cc}{\mathrm{\σ}}_n=\frac{2}{k}\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{D}_p\cos \left(\frac{2\mathrm{n\π}}{k}p\right)& (n=\mathrm{0,1},...,\stackrel{\‾}{n})\end{array}$

$\begin{array}{cc}{\mathrm{\τ}}_n=\frac{2}{k}\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{D}_p\cos \left(\frac{2\mathrm{n\π}}{k}p\right)& (n=\mathrm{0,1},...,\stackrel{\‾}{n})\end{array}$

$\begin{array}{cc}{v}_n=\frac{2}{k}\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{U}_p\cos \left(\frac{2\mathrm{n\π}}{k}p\right)& (n=\mathrm{0,1},...,\stackrel{\‾}{n})\end{array}$

$\begin{array}{cc}{\mathrm{\δ}}_n=\frac{2}{k}\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{U}_p\cos \left(\frac{2\mathrm{n\π}}{k}p\right)& (n=\mathrm{0,1},...,\stackrel{\‾}{n})\end{array}$

其中，D(p)为示功图离散点的载荷；U(p)为示功图离散点的位移；k为示功图离散点的个数；

所述吉布斯粘滞阻尼系数为：

$C=\frac{2\mathrm{\π\μ}}{{\mathrm{\ρ}}_r{A}_r}\{\frac{1}{\ln m}+\frac{2}{B_2}(B_1+1)[B_1+\frac{2}{(\mathrm{\ωL}/a)\sin (\mathrm{\ωL}/a)+\cos (\mathrm{\ωL}/a)}\left]\right\}$

m＝D_t/D_r, $B_1=\frac{m^2-1}{2\ln m}-1,B_2=m^4-1-\frac{(m^2-1)}{\ln m}$

其中：m、B₁、B₂均为中间变量；D_t为油管内径；D_r为抽油杆外径；L为抽油杆的长度；μ油井液体粘度；ρ_r为抽油杆密度；

所述将得到的泵功图进行故障诊断包括以下过程：

步骤1：采用曲线矩的特征提取方法，提取特征向量；

步骤2：给出典型故障集并建立典型故障集的物元模型：

步骤3：建立待诊断泵功图的物元模型；

步骤4：计算待诊断泵功图与标准工况的关联度函数，并计算关联度：

步骤5：对求出的关联度进行归一化处理，确定待诊断样本相对于各故障类型的发生程度；

步骤6：将归一化的关联度与油田实际生产数据进行比较分析，并结合实际的油田生产数据，确定油井的实际工况。

所述曲线矩的特征提取方法包括以下步骤：

步骤1：首先对泵功图进行归一化处理：

x'＝(x_i-x_min)/(x_max-x_min)

y'＝(y_i-y_min)/(y_max-y_min)

其中，X_i为泵功图的位移离散点；Y_i为泵功图的载荷离散点；X_max为泵功图位移的最大离散点；X_min为泵功图的最小离散点；Y_max为泵功图载荷的最大离散点；Y_min为泵功图载荷的最小离散点；

步骤2：将归一化后的泵功图根据直线l和直线m进行分区，

直线l和直线m的方程分别为：

m＝(y_max+y_min)/2

l＝str/2

其中，str为冲程；

步骤3：分区后的每一部分都是一条由离散点(x_i,y_i)，i＝1,2,...n组成的曲线，其p+q阶曲线矩m_pq定义为：

$m_{\mathrm{pq}}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i^p{y}_i^q\mathrm{\Δ}{l}_i$

$\mathrm{\Δ}{l}_i=\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}(p,q=\mathrm{0,1,2},...)$

其中：m_pq为p+q阶曲线矩；x_i为离散点横坐标；y_i为离散点纵坐标；i为曲线边缘上第i个离散点；N为离散点的总数；ΔL_i为相邻两个离散点之间的直线距离；

相应的p+q阶中心矩μ_pq的定义为：

${\mathrm{\μ}}_{\mathrm{pq}}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^p{(y_i-\stackrel{\‾}{y})}^q\mathrm{\Δ}{l}_i$

式中，μ_pq为p+q阶中心矩；点为曲线的重心坐标； $\stackrel{\‾}{y}=m_{01}/m_{00};$

步骤4：计算各阶中心距：

μ₀₀＝m₀₀,μ₁₀＝0,μ₀₁＝0,

${\mathrm{\μ}}_{11}=m_{11}-\stackrel{\‾}{x}{m}_{01}=m_{11}-\stackrel{\‾}{y}{m}_{10},$

${\mathrm{\μ}}_{20}=m_{20}-\stackrel{\‾}{x}{m}_{10},$

${\mathrm{\μ}}_{02}=m_{02}-\stackrel{\‾}{x}{m}_{01},$

${\mathrm{\μ}}_{30}=m_{30}-3\stackrel{\‾}{x}{m}_{20}+2m_{10}{\stackrel{\‾}{x}}^2,$

${\mathrm{\μ}}_{03}=m_{03}-3\stackrel{\‾}{y}{m}_{02}+2m_{01}{\stackrel{\‾}{y}}^2,$

${\mathrm{\μ}}_{12}=m_{12}-2\stackrel{\‾}{y}{m}_{11}-\stackrel{\‾}{x}{m}_{02}+2{\stackrel{\‾}{y}}^2{m}_{10},$

${\mathrm{\μ}}_{21}=m_{21}-2\stackrel{\‾}{x}{m}_{11}-\stackrel{\‾}{y}{m}_{02}+2{\stackrel{\‾}{x}}^2{m}_{01},$

其中，m₀₀为零阶矩；m₁₀,m₀₁为一阶矩；m₂₀,m₀₂,m₁₁为二阶矩；m₀₃,m₃₀,m₁₂,m₂₁为三阶矩；μ₀₀为零阶中心距；μ₀₁,μ₁₀为一阶中心距；μ₀₂,μ₂₀,μ₁₁为二阶中心距；μ₀₃,μ₃₀,μ₁₂,μ₂₁为三阶中心距；

零阶矩和零阶中心矩表示曲线的长度；一阶矩m₁₀和m₀₁可以用来确定曲线的灰度重心；二阶中心距μ₂₀，μ₁₁和μ₀₂称为惯性矩，用来衡量曲线的大小和方向；三阶中心距μ₃₀和μ₀₃表示图形的不对称性，用来衡量关于均值分布的偏差程度，μ₀₃表示曲线关于垂直轴线不对称度量；

步骤5：对中心矩进行规范化，得到如下尺度规范化矩：

η_pq＝μ_pq/(μ₀₀)^p+q+1

其中，p,q＝0,1,2,...

步骤6：在每个分区中构造7个不变矩参数：

其中，为分区的7个特征向量；η₀₂,η₂₀,η₁₁为规范化后的二阶中心距；η₀₃,η₃₀,η₁₂,η₂₁为规范化后的三阶中心距；

步骤7：计算由各个分区的不变曲线矩构成的特征向量，并利用修正公式对不变曲线矩特征向量进行修正：

其中为特征向量；为修正后的特征向量，i＝1,2,...,28。

所述典型故障集的物元模型为：

$R_i=\left[\begin{array}{ccc}{F}_i& C_1& [v_{i1a},v_{i1b}]\\ & C_2& [v_{i2a},v_{i2b}]\\ & \·& \·\\ & \·& \·\\ & \·& \·\\ & C_j& [v_{\mathrm{ija}},v_{\mathrm{ijb}}]\\ & \·& \·\\ & \·& \·\\ & \·& \·\\ & C_{28}& [v_{i28a},v_{i28b}]\end{array}\right]$

其中，R_i为典型故障集的物元模型，i＝1,2,...12，表示第i类标准工况；C_j(j＝1,2,...,28)表示第j个不变曲线矩特征；v_ij＝[v_ija,v_ijb]表示第i类工况的第j个不变曲线矩特征区间值；F_i为标准类型工况类型集合。

所述待诊断泵功图的物元模型为：

$R=\left[\begin{array}{ccc}F& C_1& v_1\\ & C_2& v_2\\ & \·& \·\\ & \·& \·\\ & \·& \·\\ & C_j& v_j\\ & \·& \·\\ & \·& \·\\ & \·& \·\\ & C_{28}& v_{28}\end{array}\right]$

其中，R为待诊断功图的物元模型，C_j(j＝1,2,...,28)表示第j个不变曲线矩特征；v_j表示第j个不变曲线矩特征值；F为待诊断泵功图的工况类型集合。

所述计算待诊断泵功图与标准工况的关联函数，并计算关联度包括以下过程：

$I\left(F_i\right)=\underset{j=1}{\overset{28}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{ij}}R{e}_{\mathrm{ij}}\left(v_j\right)$

$R{e}_{\mathrm{ij}}\left(v_j\right)=\left\{\begin{array}{c}\frac{\mathrm{\ρ}(v_j,v_{\mathrm{ij}})}{\mathrm{\ρ}(v_j,X)-\mathrm{\ρ}(v_j,v_{\mathrm{ij}})}\\ -\frac{\mathrm{\ρ}(v_j,v_{\mathrm{ij}})}{\left|v_{\mathrm{ij}}\right|}\end{array}\right.$

$\mathrm{\ρ}(v_j,v_{\mathrm{ij}})=|v_j-\frac{v_{\mathrm{ija}}+v_{\mathrm{ijb}}}{2}|-\frac{v_{\mathrm{ijb}}-v_{\mathrm{ija}}}{2}$

ω_ij＝(v_j/v_ijb)/∑(v_j/v_ijb)

其中，I(F_i)为关联度；i＝1,2,...7；j＝1,2,...,28；Re_ij(v_j)表示第j个不变曲线矩特征值下待诊断泵功图与第i类标准工况的关联函数；ρ(v_j,v_ij)表示其间距；|v_ij|表示区间之间的距离；ω_ij为各不变曲线矩特征值的权重，v_j表示第j个不变曲线矩特征值，v_ija表示第i类标准工况的第j个不变曲线矩特征区间值的下界值，v_ijb表示第i类标准工况的第j个不变曲线矩特征区间值的上界值。

所述待诊断样本相对于各故障类型的发生程度为：

$I^{\′}\left(F_i\right)=\frac{2I\left(F_i\right)-\max \left(I\left(F_i\right)\right)-\min \left(I\left(F_i\right)\right)}{\max \left(I\left(F_i\right)\right)-\min \left(I\left(F_i\right)\right)}$

其中，i＝1,2,...7；I(F_i)为关联度；I’(F_i)为修正后的关联度。

所述将归一化的关联度与油田实际生产数据进行比较分析包括以下过程：

若I(F_i)＜0，则待诊断泵功图没有发生该类型故障；

若I(F_i)≥0，则待诊断泵功图有可能属于该工况类型，并认为待诊断泵功图属于具有最大关联度的故障工况类型。

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

本发明在深入理解有杆抽油系统的工作原理基础上设计，有很好的全面性与实用性，能有效适应不同井况，使之高效准确的利用WIA-PA无线示功仪采集 的地面示功图进行有杆抽油系统的故障诊断。

附图说明

图 1是本发明的泵功图故障诊断流程图；

图 2是本发明供液不足的地面功图和泵功图；

图 3是本发明正常工况下的地面功图和泵功图；

图 4是本发明泵上碰工况的地面功图和泵功图；

图 5是本发明的泵功图分区图。

具体实施方式

下面结合附图及实施例对本发明做进一步的详细说明。

地面功图转泵功图

1)利用吉布斯提出的有杆抽油系统动态预测模型，采用数值求解的方法预测悬点动载荷的变化：

$\frac{{\∂}^{2}U(x,t)}{\∂t^2}=a^2\frac{{\∂}^{2}U(x,t)}{\∂x^2}-C\frac{\∂U(x,t)}{\∂t}$

a为声波在抽油杆中的速度，取4970m/s；C为吉布斯粘滞阻尼系数；

2)采用傅里叶变换的方法，可以求出抽油杆上任意深度、任意时间的载荷与位移：

$U(x,t)=\frac{{\mathrm{\σ}}_0}{2E{A}_r}x+\frac{v_0}{2}+\underset{n=1}{\overset{n}{\mathrm{\Σ}}}(O_n\left(x\right)\cos \mathrm{n\ωt}+P_n\left(x\right)\sin \mathrm{n\ωt})$

$F(x,t)=\frac{{\mathrm{\δ}}_0}{2}x+{\mathrm{EA}}_r\underset{n=1}{\overset{n}{\mathrm{\Σ}}}(O_n^{\′}\left(x\right)\cos \mathrm{n\ωt}+P_n^{\′}\left(x\right)\sin \mathrm{n\ωt})$

以上二式中的4个傅氏系数由实测得D(t)和U(t)曲线数值积分求得，下面以σ_n为例：首先需将D(t)和U(t)离散化，采样点数为160，等时采样。之后，令θ＝ωt, $T=\frac{2\mathrm{\π}}{\mathrm{\ω}},\mathrm{dt}=\frac{\mathrm{d\θ}}{\mathrm{\ω}},$ 可得到：

${\mathrm{\σ}}_n=\frac{1}{\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}D\left(\mathrm{\θ}\right)\cos \mathrm{n\θdt}$

考虑到θ是离散变量：

$\begin{array}{cc}\mathrm{\θ}=\frac{2\mathrm{\πp}}{k}& p=\mathrm{1,2},...,k\end{array}$

用下面符号表示：

$D\left(\frac{2\mathrm{\πp}}{k}\right)=D_p$

3)采用梯形法则进行数值积分：

$\begin{array}{c}{\mathrm{\σ}}_n=\left[\right\{\frac{D_0\cos \left(\frac{2\mathrm{n\π}\·0}{k}\right)+D_1\cos \left(\frac{2\mathrm{n\π}\·1}{k}\right)}{2}+\frac{D_1\cos \left(\frac{2\mathrm{n\π}\·1}{k}\right)+D_2\cos \left(\frac{2\mathrm{n\π}\·2}{k}\right)}{2}\\ +\frac{D_{k-1}\cos \left(\frac{2\mathrm{n\π}\·(k-1)}{k}\right)+D_k\cos \left(\frac{2\mathrm{n\π}\·k}{k}\right)}{2}\left\}\frac{2\mathrm{\π}}{k}\right]\·\frac{1}{\mathrm{\π}}\end{array}$

合并同类项得到：

${\mathrm{\σ}}_n=\frac{2}{k}[\frac{D_0\cos 0}{2}+D_1\cos \frac{2\mathrm{n\π}\·1}{k}+...+D_{k-1}\cos \left(\frac{2\mathrm{n\π}\·(k-1)}{k}\right)+\frac{D_k\cos 2\mathrm{n\π}}{2}]$

由函数的周期性可得：D₀＝D_k,cos0＝cos2nπ，因此：

$\begin{array}{cc}{\mathrm{\σ}}_n=\frac{2}{k}\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{D}_p\cos \left(\frac{2\mathrm{n\π}}{k}p\right)& (n=\mathrm{0,1},...,\stackrel{\‾}{n})\end{array}$

同理可求得其他三个傅氏级数：

$\begin{array}{cc}{\mathrm{\τ}}_n=\frac{2}{k}\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{D}_p\cos \left(\frac{2\mathrm{n\π}}{k}p\right)& (n=\mathrm{0,1},...,\stackrel{\‾}{n})\end{array}$

$\begin{array}{cc}{v}_n=\frac{2}{k}\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{U}_p\cos \left(\frac{2\mathrm{n\π}}{k}p\right)& (n=\mathrm{0,1},...,\stackrel{\‾}{n})\end{array}$

$\begin{array}{cc}{\mathrm{\δ}}_n=\frac{2}{k}\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{U}_p\cos \left(\frac{2\mathrm{n\π}}{k}p\right)& (n=\mathrm{0,1},...,\stackrel{\‾}{n})\end{array}$

4)用分离变量法求解吉布斯波动方程便可得出抽油杆深度x断面的位移随时间的变化。求解步骤如下：

首先令U(x,t)＝X(x)T(t)，X(x)和T(t)分别为仅含x和t的函数，将其带入波动方程并整理得：

$\frac{T^{\′\′}\left(t\right)}{a^{2}T\left(t\right)}+\frac{C{T}^{\′}\left(t\right)}{a^{2}T\left(t\right)}=\frac{X^{\′\′}\left(x\right)}{X\left(x\right)}$

上式每一侧仅含有一个独立变量，因而他等于一个常数，令其为代入上式并将其分离成两个常微分方程：

$\left\{\begin{array}{c}{T}^{\′\′}\left(t\right)+C{T}^{\′}\left(t\right)+{\mathrm{\λ}}_n^2{a}^{2}T\left(t\right)=0\\ X^{\′\′}\left(x\right)+{\mathrm{\λ}}_n^{2}X\left(x\right)=0\end{array}\right.$

方程组(1.12)中(1)式的周期解为：T(t)＝e^{in ω t}，代入方程组第一个等式可得方程(1.11)的特解λ_n：

λ_n＝-α_n+iβ_n

式中α_n与β_n为实数，其值为：

${\mathrm{\α}}_n=\frac{\mathrm{n\ω}}{a\sqrt{2}}\sqrt{1+\sqrt{1+{\left(\frac{C}{\mathrm{n\ω}}\right)}^2}}$

${\mathrm{\β}}_n=\frac{\mathrm{n\ω}}{a\sqrt{2}}\sqrt{-1+\sqrt{1+{\left(\frac{C}{\mathrm{n\ω}}\right)}^2}}$

当n＝0时，λ₀＝0，方程组(1.12)变为：

$\left\{\begin{array}{c}{T}^{\′\′}\left(t\right)+C{T}^{\′}\left(t\right)=0\\ X^{\′\′}\left(x\right)=0\end{array}\right.$

其解为：

$\left\{\begin{array}{c}T\left(t\right)=\mathrm{\ξ}\\ X\left(x\right)=\mathrm{\ξ}+\mathrm{\ηx}\end{array}\right.$

方程组(1.12)中第二个等式解为谐波方程：

X(x)＝φ_n sinλ_nx+θ_n cosλ_nx

综合以上各式，可求出方程(1.11)的解为：

$U(x,t)=\mathrm{\ξ}(\mathrm{\ζ}+\mathrm{\ηx})+\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}({\mathrm{\φ}}_n\sin {\mathrm{\λ}}_{n}x+{\mathrm{\θ}}_n\cos {\mathrm{\λ}}_{n}x)e^{\mathrm{in\ωt}}$

在经过变换可得到Gibbs方程的解析解为：

$U(x,t)=\frac{{\mathrm{\σ}}_0}{2E{A}_r}+\frac{v_0}{2}+\underset{n=1}{\overset{n}{\mathrm{\Σ}}}(O_n\left(x\right)\cos \mathrm{n\ωt}+P_n\left(x\right)\sin \mathrm{n\ωt})$

5)由胡克定律，抽油杆住任意深度x断面上的动载荷随时间的变化为：

$F(x,t)={\mathrm{EA}}_r\frac{\∂U(x,t)}{\∂x}$

即：

$F(x,t)=\frac{{\mathrm{\δ}}_0}{2}+{\mathrm{EA}}_r\underset{n=1}{\overset{n}{\mathrm{\Σ}}}(O_n^{\′}\left(x\right)\cos \mathrm{n\ωt}+P_n^{\′}\left(x\right)\sin \mathrm{n\ωt})$

式(1.20)和(1.22)中的函数和系数如下：

O_n(x)＝(k_nchβ_nx+δ_nshβ_nx)sinα_nx+(μ_nchβ_nx+v_nshβ_nx)cosα_nx

P_n(x)＝(k_nchβ_nx+δ_nchβ_nx)cosα_nx+(μ_nchβ_nx+v_nshβ_nx)sinα_nx

$\begin{array}{c}{O}_n^{\′}\left(x\right)=[\frac{{\mathrm{\τ}}_n}{{\mathrm{EA}}_r}\mathrm{sh}{\mathrm{\β}}_{n}x+({\mathrm{\δ}}_n{\mathrm{\β}}_n-v_n{\mathrm{\α}}_n)\mathrm{ch}{\mathrm{\β}}_{n}x]\sin {\mathrm{\α}}_{n}x\\ +[\frac{{\mathrm{\δ}}_n}{{\mathrm{EA}}_r}\mathrm{ch}{\mathrm{\β}}_{n}x+(v_n{\mathrm{\β}}_n-{\mathrm{\δ}}_n{\mathrm{\α}}_n)\mathrm{sh}{\mathrm{\β}}_{n}x]\cos {\mathrm{\α}}_{n}x\end{array}$

$\begin{array}{c}{P}_n^{\′}\left(x\right)=[\frac{{\mathrm{\τ}}_n}{{\mathrm{EA}}_r}\mathrm{ch}{\mathrm{\β}}_{n}x+({\mathrm{\δ}}_n{\mathrm{\β}}_n-v_n{\mathrm{\α}}_n)\mathrm{sh}{\mathrm{\β}}_{n}x]\cos {\mathrm{\α}}_{n}x\\ +[\frac{{\mathrm{\σ}}_n}{{\mathrm{EA}}_r}\mathrm{sh}{\mathrm{\β}}_{n}x+(v_n{\mathrm{\β}}_n-{\mathrm{\δ}}_n{\mathrm{\α}}_n)\mathrm{ch}{\mathrm{\β}}_{n}x]\sin {\mathrm{\α}}_{n}x\end{array}$

$k_n=\frac{{\mathrm{\σ}}_n{\mathrm{\α}}_n+{\mathrm{\τ}}_n{\mathrm{\β}}_n}{{\mathrm{EA}}_r({\mathrm{\α}}_n^2+{\mathrm{\β}}_n^2)}$

${\mathrm{\μ}}_n=\frac{{\mathrm{\σ}}_n{\mathrm{\β}}_n-{\mathrm{\τ}}_n{\mathrm{\α}}_n}{{\mathrm{EA}}_r({\mathrm{\α}}_n^2+{\mathrm{\β}}_n^2)}$

6)张琪等人(1984)提出了用摩擦功确立的粘滞阻尼系数公式：

$C=\frac{2\mathrm{\π\μ}}{{\mathrm{\ρ}}_r{A}_r}\{\frac{1}{\ln m}+\frac{2}{B_2}(B_1+1)[B_1+\frac{2}{(\mathrm{\ωL}/a)\sin (\mathrm{\ωL}/a)+\cos (\mathrm{\ωL}/a)}\left]\right\}$

式中：m＝D_t/D_r, $B_1=\frac{m^2-1}{2\ln m}-1,B_2=m^4-1-\frac{(m^2-1)}{\ln m}$

在油田的实际生产过程中，一般采用多级杆，所以对抽油杆柱的动载荷需要分级计算，逐级消除抽油杆柱动载荷。首先计算第一级杆柱末端载荷与位移；然后以第一级杆柱末端为第二级初始及边界条件，划分第二级抽油杆，再计算其末端载荷与位移。以此类推，直至求出最后一级末端的载荷与位移，得出最终的泵功图。

如图 2所示为利用本发明的地下供液不足时的地面示功图与地下泵功图。从图中可以看出，在上冲程部分曲线良好，基本正常，但在下冲程卸载段，曲线有明显缺失，是供液不足的典型趋势。

如图 3所示为利用本发明的泵工况正常时的地面示功图与地下泵功图。地面功图纵坐标表示悬点载荷，横坐标表示位移。曲线形状基本为平行四边形，符合正常工作的判定条件，通过地面到地下的转换，可以看出泵的工作状态良好。

如图 4所示为利用本发明的泵上碰的地面示功图与地下泵功图。从图中可以看出，在地面功图中，右上角有明显的短时间内的载荷的增大，转换成地下泵功图，可以看出，在上死点处，泵明显载荷增大，符合泵上碰的判断条件。

二故障诊断

如图 1为本发明的泵功图故障诊断流程图。包括以下步骤：

步骤1：采用曲线矩的特征提取方法，提取特征向量；

步骤2：给出典型故障集并建立典型故障集的物元模型：

步骤3：建立待诊断泵功图的物元模型；

步骤4：计算待诊断泵功图与标准工况的关联度函数，并计算关联度：

步骤5：对求出的关联度进行归一化处理，确定待诊断样本相对于各故障类型的发生程度；

步骤6：将归一化的关联度与油田实际生产数据进行比较分析，并结合实际的油田生产数据，确定油井的实际工况。

具体为：1)采用曲线矩的特征提取方法，将泵功图分为4部分，每一部分提取7个不变矩参数。

2)由于不同油井泵功图数据的精度和量纲可能不一样，需要对数据进行归一化处理：

$\stackrel{\‾}{x}=(x_i-x_{\min })/(x_{\max }-x_{\min })$

$\stackrel{\‾}{y}=(y_i-y_{\min })/(y_{\max }-y_{\min })$

3)如图 5所示为本发明的泵功图分区图，泵功图被直线l和m划分为四个部分，分别为左下、左上、右下、右上，而反映抽油杆工作状况的A(游动凡尔关 闭)、B(固定凡尔打开)、C(固定凡尔关闭)和D(固定凡尔打开)分别在这四部分中。

分别为左下，左上，右上，右下，而反映抽油泵工作状况的A(游动凡尔关闭点)，B(固定凡尔打开点)，C(固定凡尔关闭点)，D(游动凡尔打开点)分别在这四部分中；

直线l和m的方程分别为：

m＝1/2×(最大载荷+最小载荷)

l＝1/2×(冲程)

4)分别求四个分区的7个p+q阶曲线矩：

$m_{\mathrm{pq}}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i^p{y}_i^q\mathrm{\Δ}{l}_i$

$\mathrm{\Δ}{l}_i=\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}$

5)相应的求出p+q阶中心矩：

${\mathrm{\μ}}_{\mathrm{pq}}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^p{(y_i-\stackrel{\‾}{y})}^q\mathrm{\Δ}{l}_i$

η_pq＝μ_pq/(μ₀₀)^p+q+1

式中： $\stackrel{\‾}{x}=m_{10}/m_{00};\stackrel{\‾}{y}=m_{01}/m_{00};$ 点为曲线的重心坐标。

计算各阶中心距：

μ₀₀＝m₀₀,μ₁₀＝0,μ₀₁＝0,

${\mathrm{\μ}}_{11}=m_{11}-\stackrel{\‾}{x}{m}_{01}=m_{11}-\stackrel{\‾}{y}{m}_{10},$

${\mathrm{\μ}}_{20}=m_{20}-\stackrel{\‾}{x}{m}_{10},$

${\mathrm{\μ}}_{02}=m_{02}-\stackrel{\‾}{x}{m}_{01},$

${\mathrm{\μ}}_{30}=m_{30}-3\stackrel{\‾}{x}{m}_{20}+2m_{10}{\stackrel{\‾}{x}}^2,$

${\mathrm{\μ}}_{03}=m_{03}-3\stackrel{\‾}{y}{m}_{02}+2m_{01}{\stackrel{\‾}{y}}^2,$

${\mathrm{\μ}}_{12}=m_{12}-2\stackrel{\‾}{y}{m}_{11}-\stackrel{\‾}{x}{m}_{02}+2{\stackrel{\‾}{y}}^2{m}_{10},$

${\mathrm{\μ}}_{21}=m_{21}-2\stackrel{\‾}{x}{m}_{11}-\stackrel{\‾}{y}{m}_{02}+2{\stackrel{\‾}{x}}^2{m}_{01}.$

其中，m₀₀为零阶矩；m₁₀,m₀₁为一阶矩；m₂₀,m₀₂,m₁₁为二阶矩；m₀₃,m₃₀,m₁₂,m₂₁为三阶矩；μ₀₀为零阶中心距；μ₀₁,μ₁₀为一阶中心距；μ₀₂,μ₂₀,μ₁₁为二阶中心 距；μ₀₃,μ₃₀,μ₁₂,μ₂₁为三阶中心距。

所提取的泵功图曲线的各阶矩都是具有物理意义的，零阶矩和零阶中心矩表示曲线的长度；一阶矩m₁₀和m₀₁可以用来确定曲线的灰度重心；二阶中心距μ₂₀，μ₁₁和μ₀₂称为惯性矩，用来衡量曲线的大小和方向；三阶中心距μ₃₀和μ₀₃表示图形的不对称性，用来衡量关于均值分布的偏差程度，μ₀₃表示曲线关于垂直轴线不对称度量。

6)构造28个不变矩：

由各个分区的不变曲线矩构成的28个特征向量可以表示为：

利用修正公式对28个不变曲线矩特征进行修正，以调整其取值范围，所用修正公式为：

其中i＝1,2,...,28。

7)给出典型故障集并建立故障集的物元模型：

$R_i=\left[\begin{array}{ccc}{F}_i& C_1& [v_{i1a},v_{i1b}]\\ & C_2& [v_{i2a},v_{i2b}]\\ & \·& \·\\ & \·& \·\\ & \·& \·\\ & C_j& [v_{\mathrm{ija}},v_{\mathrm{ijb}}]\\ & \·& \·\\ & \·& \·\\ & \·& \·\\ & C_{28}& [v_{i28a},v_{i28b}]\end{array}\right]$

其中，i＝1,2,...7，表示第i类标准工况；C_j(j＝1,2,...,28)表示第j个不变曲线矩特征；v_ij＝[v_ija,v_ijb]表示第i类工况的第j个不变曲线矩特征区间值。

8)给出待诊断泵功图物元模型

$R=\left[\begin{array}{ccc}F& C_1& v_1\\ & C_2& v_2\\ & \·& \·\\ & \·& \·\\ & \·& \·\\ & C_j& v_j\\ & \·& \·\\ & \·& \·\\ & \·& \·\\ & C_{28}& v_{28}\end{array}\right]$

其中，C_j(j＝1,2,...,28)表示第j个不变曲线矩特征；v_j表示第j个不变曲线矩特征值。

9)计算待诊断泵功图与各类型标准工况的关联函数：

$R{e}_{\mathrm{ij}}\left(v_j\right)=\left\{\begin{array}{c}\frac{\mathrm{\ρ}(v_j,v_{\mathrm{ij}})}{\mathrm{\ρ}(v_j,X)-\mathrm{\ρ}(v_j,v_{\mathrm{ij}})}\\ -\frac{\mathrm{\ρ}(v_j,v_{\mathrm{ij}})}{\left|v_{\mathrm{ij}}\right|}\end{array}\right.$

$\mathrm{\ρ}(v_j,v_{\mathrm{ij}})=|v_j-\frac{v_{\mathrm{ija}}+v_{\mathrm{ijb}}}{2}|-\frac{v_{\mathrm{ijb}}-v_{\mathrm{ija}}}{2}$

其中，i＝1,2,...7；j＝1,2,...,28；Re_ij(v_j)表示第j个不变曲线矩特征值下待诊断示功图与第i类标准工况的关联函数；ρ(v_j,v_ij)表示其间距；|v_ij|表示区间之间的距离。

10)计算关联度

$I\left(F_i\right)=\underset{j=1}{\overset{28}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{ij}}{K}_{\mathrm{ij}}$

其中，ωⁱ _j为各不变曲线矩特征值的权重。

ω_ij＝(v_j/v_ijb)/∑(v_j/v_ijb)

其中，：i＝1,2,...7；j＝1,2,...,28；v_j表示第j个不变曲线矩特征值，；v_ijb表示第i类标准工况的第j个不变曲线矩特征区间值的上界值。

11)对关联度进行归一化处理，确定待诊断样本相对于各故障类型的发生程度；

$I^{\′}\left(F_i\right)=\frac{2I\left(F_i\right)-\max \left(I\left(F_i\right)\right)-\min \left(I\left(F_i\right)\right)}{\max \left(I\left(F_i\right)\right)-\min \left(I\left(F_i\right)\right)}$

其中，i＝1,2,...7。

12)将计算的关联度与油田实际生产数据相结合

若I(F_i)＜0，则认为待诊断泵功图没有发生该类型故障；若I(F_i)≥0，说明待诊断泵功图有可能属于该工况类型；若所有的I(F_i)＜0，认为待诊断泵功图不在现有的类型之中，将I(F_i)≥0的工况类型与实际的油田生产数据相结合，确定最终的油井实际工况。

实施例：

表 1中给出了各故障子集下每个特征元的取值区间，其中故障类型1-15分别代表“正常工作”，“气体影响”，“供液不足”，“抽油杆断脱”，“油稠”，“游动凡尔漏失”，“泵上碰”，“泵下碰”，“固定凡尔漏失”，“柱塞脱出工作筒”，“气锁”，“双凡尔漏失”，“卡泵”，“油井出砂”，“结蜡”。

表 1各类故障类型的特征元区间

设故障级为F＝{F₁,F₂,…F₁₅}，每个故障子集的特征为：c＝{c₁,c₂,…c₁₅}。构造个故障子集的物元模型如下：

给出待诊断样本，提取待诊断样本的不变曲线矩特征如表 2所示：

表 2待诊断样本的不变曲线矩特征向量

建立待诊断样本的物元模型如下：

计算关联度结果如表 3所示：

表 3待诊断样本关联度计算

由表可知，对于样本一来说，I(F_i)<0表示不属于此类型故障。I(F₄)和I(F₁₃)都大于0，判断待诊断样本有可能属于故障类型4(“供液不足”)或故障类型(“卡泵”)，程度分别为：0.0182和0.0066。结合油田实际生产数据可知，该井的动液面为1304m，并且示功图形状与供液不足时的典型工况一样，因此该口井的工况可以诊断为“供液不足”。对于样本二来说，I(F_i)<0表示不属于此类型故障。I(F₇)和I(F₁₅)都大于0，判断待诊断样本有可能属于故障类型7(“泵上碰”)或故障类型12(“凡尔漏失”)，程度分别为0.0212和0.0102。对示功图图形进行分析可以得到，图形的右上角明显“凸出”，这正是“泵上碰”所具有的特征，结合油田实际生产数据可知，该井的最大理论载荷为49.9KN，最小理论载荷为30.7KN，因此不具有“凡尔漏失”所具有的特征。因此该口井的工况可以诊断为“泵上碰”。

Claims (10)

1.一种基于泵功图的抽油井故障诊断方法，其特征在于：采集地面示功图，并将示功图转化为泵功图；将得到的泵功图进行故障诊断。

2.根据权利要求1所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述将示功图转化为泵功图包括以下步骤：

步骤1：利用有杆抽油系统动态参数预测模型，预测悬点动载荷的变化：

其中，a为声波在抽油杆中的速度；C为吉布斯粘滞阻尼系数；

步骤2：采用傅里叶变换的方法，求出抽油杆上任意深度x、任意时间t的载荷与位移，得到泵功图：

其中，U(x,t)为抽油杆的位移；F(x,t)为抽油杆的载荷；σ₀,δ₀为傅氏系数；E为抽油杆的弹性模量；A_r为抽油杆截面积；O_n,O'_n,P_n,P'_n为特殊函数；n为傅氏级数；ω为曲柄角速度；n为傅氏级数。

3.根据权利要求2所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述特殊函数为：

O_n(x)＝(k_nchβ_nx+δ_nshβ_nx)sinα_nx+(μ_nchβ_nx+v_nshβ_nx)cosα_nx

P_n(x)＝(k_nchβ_nx+δ_nchβ_nx)cosα_nx+(μ_nchβ_nx+v_nshβ_nx)sinα_nx

其中，O_n,O'_n,P_n,P'_n为特殊函数；σ_n,δ_n,τ_n,v_n为傅氏系数；E为抽油杆的弹性模量；A_r为抽油杆截面积；α_n,β_n,k_n,μ_n,a_n,b_n,c_n,d_n均为特殊常数；ω为曲柄角速度；a为声波在抽油杆中的速度；C为吉布斯粘滞阻尼系数；n为傅氏级数；

所述傅氏系数为：

其中，D(p)为示功图离散点的载荷；U(p)为示功图离散点的位移；k为示功图离散点的个数；

所述吉布斯粘滞阻尼系数为：

其中：m、B₁、B₂均为中间变量；D_t为油管内径；D_r为抽油杆外径；L为抽油杆的长度；μ油井液体粘度；ρ_r为抽油杆密度。

4.根据权利要求1所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述将得到的泵功图进行故障诊断包括以下过程：

步骤1：采用曲线矩的特征提取方法，提取特征向量；

步骤2：给出典型故障集并建立典型故障集的物元模型：

步骤3：建立待诊断泵功图的物元模型；

步骤4：计算待诊断泵功图与标准工况的关联函数，并计算关联度：

步骤5：对求出的关联度进行归一化处理，确定待诊断样本相对于各故障类型的发生程度；

步骤6：将归一化的关联度与油田实际生产数据进行比较分析，并结合实际的油田生产数据，确定油井的实际工况。

5.根据权利要求4所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述曲线矩的特征提取方法包括以下步骤：

步骤1：首先对泵功图进行归一化处理：

x'＝(x_i-x_min)/(x_max-x_min)

y'＝(y_i-y_min)/(y_max-y_min)

其中，X_i为泵功图的位移离散点；Y_i为泵功图的载荷离散点；X_max为泵功图位移的最大离散点；X_min为泵功图的最小离散点；Y_max为泵功图载荷的最大离散点；Y_min为泵功图载荷的最小离散点；

步骤2：将归一化后的泵功图根据直线l和直线m进行分区，

直线l和直线m的方程分别为：

m＝(y_max+y_min)/2

l＝str/2

其中，str为冲程；

步骤3：分区后的每一部分都是一条由离散点(x_i,y_i)，i＝1,2,...n组成的曲线，其p+q阶曲线矩m_pq定义为：

其中：m_pq为p+q阶曲线矩；x_i为离散点横坐标；y_i为离散点纵坐标；i为曲线边缘上第i个离散点；N为离散点的总数；ΔL_i为相邻两个离散点之间的直线距离；

相应的p+q阶中心矩μ_pq的定义为：

式中，μ_pq为p+q阶中心矩；点为曲线的重心坐标；

步骤4：计算各阶中心距：

μ₀₀＝m₀₀,μ₁₀＝0,μ₀₁＝0,

其中，m₀₀为零阶矩；m₁₀,m₀₁为一阶矩；m₂₀,m₀₂,m₁₁为二阶矩；m₀₃,m₃₀,m₁₂,m₂₁为三阶矩；μ₀₀为零阶中心距；μ₀₁,μ₁₀为一阶中心距；μ₀₂,μ₂₀,μ₁₁为二阶中心距；μ₀₃,μ₃₀,μ₁₂,μ₂₁为三阶中心距；

零阶矩和零阶中心矩表示曲线的长度；一阶矩m₁₀和m₀₁可以用来确定曲线的灰度重心；二阶中心距μ₂₀，μ₁₁和μ₀₂称为惯性矩，用来衡量曲线的大小和方向；三阶中心距μ₃₀和μ₀₃表示图形的不对称性，用来衡量关于均值分布的偏差程度，μ₀₃表示曲线关于垂直轴线不对称度量；

步骤5：对中心矩进行规范化，得到如下尺度规范化矩：

η_pq＝μ_pq/(μ₀₀)^p+q+1

其中，p,q＝0,1,2,...

步骤6：在每个分区中构造7个不变矩参数：

其中，为分区的7个特征向量；η₀₂,η₂₀,η₁₁为规范化后的二阶中心距；η₀₃,η₃₀,η₁₂,η₂₁为规范化后的三阶中心距；

步骤7：计算由各个分区的不变曲线矩构成的特征向量，并利用修正公式对不变曲线矩特征向量进行修正：

其中为特征向量；为修正后的特征向量，i＝1,2,...,28。

6.根据权利要求4所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述典型故障集的物元模型为：

其中，R_i为典型故障集的物元模型，i＝1,2,...12，表示第i类标准工况；C_j(j＝1,2,...,28)表示第j个不变曲线矩特征；v_ij＝[v_ija,v_ijb]表示第i类工况的第j个不变曲线矩特征区间值；F_i为标准类型工况类型集合。

7.根据权利要求4所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述待诊断泵功图的物元模型为：

其中，R为待诊断功图的物元模型，C_j(j＝1,2,...,28)表示第j个不变曲线矩特征；v_j表示第j个不变曲线矩特征值；F为待诊断泵功图的工况类型集合。

8.根据权利要求4所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述计算待诊断泵功图与标准工况的关联函数，并计算关联度包括以下过程：

其中，I(F_i)为关联度；i＝1,2,...7；j＝1,2,...,28；Re_ij(v_j)表示第j个不变 曲线矩特征值下待诊断泵功图与第i类标准工况的关联函数；ρ(v_j,v_ij)表示其间距；|v_ij|表示区间之间的距离；ω_ij为各不变曲线矩特征值的权重，v_j表示第j个不变曲线矩特征值，v_ija表示第i类标准工况的第j个不变曲线矩特征区间值的下界值，v_ijb表示第i类标准工况的第j个不变曲线矩特征区间值的上界值。

9.根据权利要求4所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述待诊断样本相对于各故障类型的发生程度为：

其中，i＝1,2,...7；I(F_i)为关联度；I’(F_i)为修正后的关联度。

10.根据权利要求4所述的基于泵功图的抽油井故障诊断方法，其特征在于：所述将归一化的关联度与油田实际生产数据进行比较分析包括以下过程：

若I(F_i)＜0，则待诊断泵功图没有发生该类型故障；

若I(F_i)≥0，则待诊断泵功图有可能属于该工况类型，并认为待诊断泵功图属于具有最大关联度的故障工况类型。