CN105257277A - 基于多变量灰色模型的有杆泵抽油井井下故障预测方法 - Google Patents

Abstract

一种基于多变量灰色模型的有杆泵抽油井井下故障预测方法，由示功图远程数据采集***采集地面示功图；将地面示功图转化为井下泵功图；根据曲线矩理论提取能够表征井下泵功图图形特征的曲线矩特征向量；由多变量灰色模型对每一个特征向量进行预测；计算预测样本与标准集中每一种故障类型的灰关联度，最大灰关联度所对应的故障类型即为预测样本的故障类型。可合理调整抽油井的工作制度，从而提高抽油井的生产效率，可实现对有杆泵抽油井井下工作状况的故障预测；计算简单，所需要的样本少，并且对于任意分布的数据集都适用；预测的准确度高，特别对于渐变型的故障类型具有较高的预测精度。

Description

基于多变量灰色模型的有杆泵抽油井井下故障预测方法

技术领域

本发明属于有杆泵抽油井故障预测方法领域，具体涉及一种基于多变量灰色模型的有杆泵抽油井井下故障预测方法。

背景技术

有杆泵抽油井的故障是油田生产所面临的一个主要问题，会影响井下抽油泵的运行状况和油井的产油量。由于抽油泵工作在数千米的井下，工作状况十分复杂，工作环境极其恶劣，故障发生率很高，很大程度会影响油田的正常生产。一旦抽油井发生了故障而不能进行及时诊断，就会造成能源的浪费，并且影响抽油井生产，给企业带来损失。目前已有很多研究关注于有杆泵抽油井井下的故障诊断方法，但是这是建立在故障已经发生的基础上，已经给抽油井的正常生产带来了一定的影响。

发明内容

本发明要解决的技术问题是提供一种基于多变量灰色模型的有杆泵抽油井井下故障预测方法，实现对有杆泵抽油井井下工作状况的故障预测。

本发明的技术解决方案是：

一种基于多变量灰色模型的有杆泵抽油井井下故障预测方法，其特殊之处是：

1)、采用示功图远程数据采集***，首先，将无线示功仪安装在有杆泵抽油井的抽油机的驴头上，采集光杆的载荷和位移数据；无线RTU以无线的方式接收无线示功仪采集的数据，由数据采集模块通过无线AP从站模块以无线网络远程传送到无线AP主站模块上；示功图监测及数据处理服务器接收采集到的光杆的载荷和位移数据，绘制采集到的地面示功图图形；

2)、将地面示功图转化为井下泵功图，以得到能够真实反映抽油泵工作状况的示功图图形；

3)、提取井下泵功图的曲线矩特征作为特征向量，将井下泵功图划分为上冲程部分和下冲程部分，然后根据曲线矩理论分别提取上冲程部分的7个曲线矩特征向量和下冲程部分的7个曲线矩特征向量，将得到的14个曲线矩特征向量作为下一步预测的变量，具体步骤如下：

a)将井下泵功图划分为上冲程部分和下冲程部分，其中上冲程部分曲线反映了有杆泵抽油井的上冲程过程，是抽油杆柱在动力装置的作用下向上拉动柱塞，抽油泵内吸入油液而井口排出油液的过程；下冲程部分曲线反映了有杆泵抽油井的下冲程过程，是抽油杆柱在动力装置的作用下向下推动柱塞，抽油泵向油管内排出油液的过程；假设井下泵功图由N个采样点组成，那么从第1个点到位移最大的点(假设为第l个点)为上冲程部分；从第l+1个点到第N个点为下冲程部分；

b)划分后的上冲程部分曲线由l个采样点组成，每个采样点的坐标为(x_i,y_i)，其s+r阶曲线矩m_sr定义为：

$m_{\mathrm{sr}}=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{x_i}^s{y_i}^r\mathrm{\Δ}{z}_i---\left(1\right)$

其中i＝1,…,l；s,r＝0,1,2；Δz_i为两个连续的采样点间的距离，

μ_sr为s+r阶中心矩，定义为：

${\mathrm{\μ}}_{\mathrm{sr}}=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\left(x_i\text{-}\stackrel{\‾}{x}\right)}^s{(y_i-\stackrel{\‾}{y})}^r\mathrm{\Δ}{z}_i---\left(2\right)$

其中 $\stackrel{\‾}{x}=\frac{m_{10}}{m_{00}},\stackrel{\‾}{y}=\frac{m_{01}}{m_{00}};$ 点为曲线的重心坐标；

计算出如下各阶中心矩：μ₀₀＝m₀₀，μ₁₀＝0，μ₀₁＝0， ${\mathrm{\μ}}_{11}=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};$

对提取的各阶中心矩进行规范化处理，由η_sr表示，定义为：

η_sr＝μ_sr/(μ₀₀)^s+r+1(3)

那么构造的上冲程部分的7个曲线矩特征向量分别定义如下：

采用如下修正公式使构造的上冲程部分的7个曲线矩特征向量的取值范围统一，记为ψ₁—ψ₇，定义为：

其中q＝0,1,…,7；

c)提取下冲程部分的曲线矩特征向量，划分后的下冲程部分曲线由N-l个采样点组成，每个采样点的坐标为(x_j,y_j)，其s+r阶曲线矩m'_sr定义为：

$m_{\mathrm{sr}}^{\′}=\underset{j=l+1}{\overset{N}{\mathrm{\Σ}}}{x_j}^s{y_j}^r\mathrm{\Δ}{z}_j---\left(12\right)$

其中j＝l+1,l+2,…,N；s,r＝0,1,2；Δz_j为两个连续的采样点间的距离，

μ'_sr为s+r阶中心矩，定义为：

${\mathrm{\μ}}_{\mathrm{sr}}^{\′}=\underset{j=l+1}{\overset{N}{\mathrm{\Σ}}}{(x_j-\stackrel{\‾}{x^{\′}})}^s{(y_j-\stackrel{\‾}{y^{\′}})}^r\mathrm{\Δ}{z}_j---\left(13\right)$

其中 $\stackrel{\‾}{x^{\′}}=\frac{m_{10}}{m_{00}},\stackrel{\‾}{y^{\′}}=\frac{m_{01}}{m_{00}};$ 点为曲线的重心坐标；

计算出如下各阶中心矩：μ′₀₀＝m′₀₀，μ′₁₀＝0，μ′₀₁＝0， ${\mathrm{\μ}}_{11}^{\′}=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}^{\′};$

对提取的各阶中心矩进行规范化处理，由η'_sr表示，定义为：

η'_sr＝μ'_sr/(μ'₀₀)^s+r+1(14)

那么构造的下冲程部分的7个曲线矩特征向量分别定义如下：

采用如下修正公式使构造的下冲程部分的7个曲线矩特征向量的取值范围统一，记为ψ₈—ψ₁₄，定义为：

其中q’＝8,9,…,14；

将提取得到的能够表征井下泵功图图形特征的14个特征向量ψ₁—ψ₁₄，作为下一步预测的变量；

4)、根据多变量灰色模型MGM(1,n)的基本理论建立有杆泵抽油井井下故障的多变量灰色预测模型，利用MGM(1,n)模型建立多个变量之间的动态相关性，弥补单变量方法的局限性；将一段连续时间内的泵功图所提取的14个曲线矩特征向量作为多变量灰色模型(MGM(1,14))的输入，建立灰时间序列，每一个曲线矩特征向量作为灰时间序列的一个变量；由所建立的多变量灰色预测模型分别预测14个曲线矩特征向量的下一个时间点的特征向量，所述14个预测的曲线矩特征向量即为能够表征预测样本特征的特征向量；

将所建立的泵功图曲线矩特征向量灰时间序列表示如下：

${\mathrm{\ψ}}^{\left(0\right)}\left(k\right)=\left[\begin{array}{cccccc}{\mathrm{\ψ}}_1^{\left(0\right)}\left(1\right)& {\mathrm{\ψ}}_1^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_1^{\left(0\right)}\left(k\right)& ...& {\mathrm{\ψ}}_1^{\left(0\right)}\left(m\right)\\ {\mathrm{\ψ}}_2^{\left(0\right)}\left(1\right)& {\mathrm{\ψ}}_2^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_2^{\left(0\right)}\left(k\right)& ...& {\mathrm{\ψ}}_2^{\left(0\right)}\left(m\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ {\mathrm{\ψ}}_p^{\left(0\right)}\left(1\right)& {\mathrm{\ψ}}_p^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_p^{\left(0\right)}\left(k\right)& ...& {\mathrm{\ψ}}_p^{\left(0\right)}\left(m\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(1\right)& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(k\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(m\right)\end{array}\right]---\left(23\right)$

其中p＝1,2,…,14；k＝1,2,…,m，m为所采用的泵功图的数量；表示泵功图曲线矩特征向量灰时间序列中的第p个变量中的第k个影响因素；

令为相应的一次累加生成序列，有：

${\mathrm{\ψ}}_p^{\left(1\right)}\left(k\right)=\underset{g=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\ψ}}_p^{\left(0\right)}\left(g\right)---\left(24\right)$

建立14元一阶常微分方程，有：

$\left\{\begin{array}{c}\frac{d{\mathrm{\ψ}}_1^{\left(1\right)}}{\mathrm{dt}}=a_{1-1}{\mathrm{\ψ}}_1^{\left(1\right)}+a_{1-2}{\mathrm{\ψ}}_2^{\left(1\right)}+...+a_{1-14}{\mathrm{\ψ}}_{14}^{\left(1\right)}+b_1\\ \frac{d{\mathrm{\ψ}}_2^{\left(1\right)}}{\mathrm{dt}}=a_{2-1}{\mathrm{\ψ}}_1^{\left(1\right)}+a_{2-2}{\mathrm{\ψ}}_2^{\left(1\right)}+...+a_{2-14}{\mathrm{\ψ}}_{14}^{\left(1\right)}+b_2\\ .\\ .\\ .\\ \frac{d{\mathrm{\ψ}}_{14}^{\left(1\right)}}{\mathrm{dt}}=a_{14-1}{\mathrm{\ψ}}_1^{\left(1\right)}+a_{14-2}{\mathrm{\ψ}}_2^{\left(1\right)}+...+a_{14-14}{\mathrm{\ψ}}_{14}^{\left(1\right)}+b_{14}\end{array}\right.---\left(25\right)$

记

$A=\left[\begin{array}{cccc}{a}_{1-1}& a_{1-2}& ...& a_{1-14}\\ a_{2-1}& a_{2-2}& ...& a_{2-14}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{14-1}& a_{14-2}& ...& a_{14-14}\end{array}\right]---\left(26\right)$

B＝(b₁,b₁,…,b_n)^T(27)

A和B为辨识参数，那么将公式(25)简记为

$\frac{d{\mathrm{\Ψ}}^{\left(1\right)}}{\mathrm{dt}}=A{\mathrm{\Ψ}}^{\left(1\right)}+B---\left(28\right)$

其中 ${\mathrm{\Ψ}}^{\left(1\right)}={({\mathrm{\ψ}}_1^{\left(1\right)},{\mathrm{\ψ}}_2^{\left(1\right)},...,{\mathrm{\ψ}}_{14}^{\left(1\right)})}^T;$

在[0,t]区间上积分，得到连续时间响应为

${\mathrm{\Ψ}}^{\left(1\right)}\left(t\right)=e^{\mathrm{At}}{\mathrm{\Ψ}}^{\left(1\right)}\left(0\right)+A^{-1}(e^{\mathrm{At}}-I)B---\left(29\right)$

其中 $e^{\mathrm{At}}=I+\underset{u=1}{\overset{\∞}{\mathrm{\Σ}}}\frac{A^u}{u!}{t}^u;$

将公式(25)进行离散化，得到：

${\mathrm{\ψ}}_p^{\left(0\right)}\left(k\right)=\underset{w=1}{\overset{14}{\mathrm{\Σ}}}\frac{a_{p-w}}{2}({\mathrm{\ψ}}_w^{\left(1\right)}\left(k\right)+{\mathrm{\ψ}}_w^{\left(1\right)}(k-1))---\left(30\right)$

其中p＝1,2,…,14；k＝2,3,…,m；

记

$E=\left[\begin{array}{cccc}{a}_{1-1}& a_{2-1}& ...& a_{14-1}\\ a_{1-2}& a_{2-2}& ...& a_{14-2}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{1-14}& a_{2-14}& ...& a_{14-14}\\ b_1& b_2& ...& b_{14}\end{array}\right]---\left(31\right)$

由最小二乘法得到E的辨识值E，有：

E＝(L^TL)^-1L^TY(32)

其中，

$Y=\left[\begin{array}{cccc}{\mathrm{\ψ}}_1^{\left(0\right)}\left(2\right)& {\mathrm{\ψ}}_2^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(2\right)\\ {\mathrm{\ψ}}_1^{\left(0\right)}\left(3\right)& {\mathrm{\ψ}}_2^{\left(0\right)}\left(3\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(3\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\ψ}}_1^{\left(0\right)}\left(m\right)& {\mathrm{\ψ}}_2^{\left(0\right)}\left(m\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(m\right)\end{array}\right]---\left(34\right)$

$E={\left[\begin{array}{cc}A& B\end{array}\right]}^T=\left[\begin{array}{cccc}{a}_{1-1}& a_{2-1}& ...& a_{14-1}\\ a_{1-2}& a_{2-2}& ...& a_{14-2}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{1-14}& a_{2-14}& ...& a_{14-14}\\ {\tilde{b}}_1& {\tilde{b}}_2& ...& {\tilde{b}}_{14}\end{array}\right]---\left(35\right)$

由公式(35)得到辨识参数A和B的辨识值A和B，分别为：

$A=\left[\begin{array}{cccc}{a}_{1-1}& a_{1-2}& ...& a_{1-14}\\ a_{2-1}& a_{2-2}& ...& a_{2-14}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{14-1}& a_{14-2}& ...& a_{14-14}\end{array}\right]---\left(36\right)$

$B=\left[\begin{array}{c}{\tilde{b}}_1\\ {\tilde{b}}_2\\ .\\ .\\ .\\ {\tilde{b}}_{14}\end{array}\right]---\left(37\right)$

根据得到的辨识值A和B，将公式(29)写成如下离散形式，

Ψ⁽⁰⁾(1)＝Ψ⁽⁰⁾(1)(38)

Ψ⁽¹⁾(k)＝e^A(k-1)Ψ⁽¹⁾(1)+A^-1(e^A(k)-I)B(39)

其中Ψ⁽¹⁾(k)和Ψ⁽¹⁾(1)可以将公式(36)和公式(38)代入公式(30)求得；

得到预测模型为：

Ψ⁽⁰⁾(k)＝Ψ⁽¹⁾(k)-Ψ⁽¹⁾(k-1)(40)

其中k＝2,3,…,m；

5、最后，判断预测样本属于哪类故障类型通过已有的示功图信息建立每一种故障类型的标准集，即对已确定故障类型的一定数量的示功图通过将地面示功图转化为井下泵功图和分别提取其14个曲线矩特征向量，得到每一种故障类型中每一个曲线矩特征向量的取值区间；所得到的取值区间由表示，p＝1,2,…,14，其中表示第F_M种故障类型的第p个曲线矩特征向量的下限值，表示第F_M种故障类型的第p个曲线矩特征向量的上限值；

计算预测样本的每一个曲线矩特征向量值与每一种故障类型中相应的曲线矩特征向量值的距离，由下式进行计算：

$d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})=|U_{{\mathrm{\ψ}}_p}-\frac{V_{F_M\mathrm{pa}}+V_{F_M\mathrm{pb}}}{2}|+\frac{V_{F_M\mathrm{pb}}-V_{F_M\mathrm{pa}}}{2}---\left(41\right)$

然后由下式计算预测样本与每一种故障类型的灰关联系数，

${\mathrm{\ξ}}_{p{F}_M}=\frac{\underset{M}{\min }\underset{p}{\min }\left\{d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})\right\}+\mathrm{\ρ}\underset{M}{\max }\underset{p}{\max }\left\{d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})\right\}}{d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})+\mathrm{\ρ}\underset{M}{\max }\underset{p}{\max }\left\{d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})\right\}}---\left(42\right)$

其中p＝1,2,…,14；M＝1,2,3,…。ρ∈[0,1]为分辨系数，在这里取0.5；

再由下式计算预测样本与每一种故障类型的灰关联度，

$\mathrm{Gre}{y}_{F_M}=\frac{1}{p}\underset{p=1}{\overset{14}{\mathrm{\Σ}}}{\mathrm{\ξ}}_{p{F}_M}---\left(43\right)$

最大的灰关联度所对应的故障类型即为预测样本的故障类型。

进一步的，步骤3)中井下泵功图的上冲程部分为有杆泵抽油井抽油杆的悬点由最低点运动到最高点的部分；下冲程部分为有杆泵抽油井抽油杆的悬点由最高点运动到最低点的部分。

进一步的，步骤4)对于有杆泵抽油井井下故障的多变量灰色预测模型中所采用的一段连续时间内的样本，连续时间的单位可以为小时或天，所述样本的数量为50－200个，以提高预测模型的计算效率和准确度。

进一步的，步骤5)由已有的信息建立每一种故障类型的标准集，已有的信息由已知道故障类型的示功图得到，用来建立每一种故障类型的标准集的示功图的数量大于2个。

进一步的，所述有杆泵抽油井井下故障类型有以下9种：“正常”、“气体影响”、“供液不足”、“抽油杆断落”、“游动阀漏失”、“固定阀漏失”、“泵上碰”、“泵下碰”、“油井出砂”。

本发明的有益效果是：

1、有杆泵抽油井是一个复杂的机械***，具有很强的非线性，其故障具有随机性、不确定性和相对性。本发明发掘并利用有杆泵抽油井生产中的有用信息，对有杆泵抽油井的故障进行预测，进而合理调整抽油井的工作制度，从而提高抽油井的生产效率，这对保证油田正常生产是具有十分重要的意义的。

2、示功图是油田生产中分析有杆泵抽油井井下工作状况的主要方法，是采集光杆位移变化一个完整周期内载荷和位移数据并绘制的光杆载荷和光杆位移关系的曲线，示功图图形的不同形状反映了有杆泵抽油井井下不同的工作状况。由于有杆泵抽油井的生产过程是一个连续的过程，其井下的工作状况也是一个渐变的过程，因此，可以通过采集一段时间内的有杆泵抽油井的示功图，根据其图形形状的渐变过程实现对有杆泵抽油井井下工作状况的故障预测。

3、本发明所建立的有杆泵抽油井井下故障的多变量灰色预测模型，计算简单，所需要的样本少，并且对于任意分布的数据集都适用。根据所建立的各故障类型的标准集进行预测样本的判断，准确度高，特别对于渐变型的故障类型具有较高的预测精度。

附图说明

图1是本发明的示功图远程数据采集***示意图；

图2是本发明的工作原理图；

图3是采集的地面示功图和转化后的井下泵功图示意图；

图4是井下泵功图的划分示意图；

图5是曲线矩特征向量Ψ1的灰色模型预测结果示意图；

图6是曲线矩特征向量Ψ2的灰色模型预测结果示意图；

图7是曲线矩特征向量Ψ3的灰色模型预测结果示意图；

图8是曲线矩特征向量Ψ4的灰色模型预测结果示意图；

图9是曲线矩特征向量Ψ5的灰色模型预测结果示意图；

图10是曲线矩特征向量Ψ6的灰色模型预测结果示意图；

图11是曲线矩特征向量Ψ7的灰色模型预测结果示意图；

图12是曲线矩特征向量Ψ8的灰色模型预测结果示意图；

图13是曲线矩特征向量Ψ9的灰色模型预测结果示意图；

图14是曲线矩特征向量Ψ10的灰色模型预测结果示意图；

图15是曲线矩特征向量Ψ11的灰色模型预测结果示意图；

图16是曲线矩特征向量Ψ12的灰色模型预测结果示意图；

图17是曲线矩特征向量Ψ13的灰色模型预测结果示意图；

图18是曲线矩特征向量Ψ14的灰色模型预测结果示意图。

具体实施方式

一种基于多变量灰色模型的有杆泵抽油井井下故障预测方法，其步骤如下：

1、如图1和图2所示，采用示功图远程数据采集***，首先将无线示功仪安装在有杆泵抽油井的抽油机的驴头上，采集光杆的载荷和位移数据；无线RTU以无线的方式接收无线示功仪采集的数据，由数据采集模块经交换机2通过无线AP从站模块以无线网络远程传送到无线AP主站模块上，再经交换机1传至示功图监测及数据处理服务器；示功图监测及数据处理服务器接收采集到的光杆的载荷和位移数据，绘制采集到的地面示功图图形。

2、如图3所示，将地面示功图转化为井下泵功图，以得到能够真实反映抽油泵工作状况的示功图图形；由于受到抽油柱变形及振动等因素的影响，地面示功图并不能真实反映有杆泵抽油井井下的实际工作状况。因此，通过将地面示功图转化为井下泵功图，以消除这些影响。

3、由于井下泵功图的图形形状反映了有杆泵抽油井井下的工作状况，因此，若要实现对井下故障的预测，可以先提取能够表征井下泵功图图形特征的若干特征向量，通过一段时间内各特征向量的变化趋势来预测未来的特征向量，从而得到该特征向量所对应的故障类型。本发明提取井下泵功图的曲线矩特征作为特征向量，具体步骤如下：

a)将图3中的井下泵功图划分为上冲程部分和下冲程部分，井下泵功图的上冲程部分为有杆泵抽油井抽油杆的悬点由最低点运动到最高点的部分；下冲程部分为有杆泵抽油井抽油杆的悬点由最高点运动到最低点的部分。如图4所示，其中上冲程部分曲线反映了有杆泵抽油井的上冲程过程，是抽油杆柱在动力装置的作用下向上拉动柱塞，抽油泵内吸入油液而井口排出油液的过程；下冲程部分曲线反映了有杆泵抽油井的下冲程过程，是抽油杆柱在动力装置的作用下向下推动柱塞，抽油泵向油管内排出油液的过程。假设井下泵功图由N个采样点组成，那么从第1个点到位移最大的点(假设为第l个点)为上冲程部分；从第l+1个点到第N个点为下冲程部分。

b)划分后的上冲程部分曲线由l个采样点组成，每个采样点的坐标为(x_i,y_i)，其s+r阶曲线矩m_sr定义为：

$m_{\mathrm{sr}}=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{x_i}^s{y_i}^r\mathrm{\Δ}{z}_i---\left(1\right)$

其中i＝1,…,l；s,r＝0,1,2；Δz_i为两个连续的采样点间的距离，

μ_sr为s+r阶中心矩，定义为：

${\mathrm{\μ}}_{\mathrm{sr}}=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\left(x_i\text{-}\stackrel{\‾}{x}\right)}^s{(y_i-\stackrel{\‾}{y})}^r\mathrm{\Δ}{z}_i---\left(2\right)$

其中 $\stackrel{\‾}{x}=\frac{m_{10}}{m_{00}},\stackrel{\‾}{y}=\frac{m_{01}}{m_{00}};$ 点为曲线的重心坐标；

计算出如下各阶中心矩：μ₀₀＝m₀₀，μ₁₀＝0，μ₀₁＝0， ${\mathrm{\μ}}_{11}=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}.$

对提取的各阶中心矩进行规范化处理，由η_sr表示，定义为：

η_sr＝μ_sr/(μ₀₀)^s+r+1(3)

那么构造的上冲程部分的7个曲线矩特征向量分别定义如下：

采用如下修正公式使构造的上冲程部分的7个曲线矩特征向量的取值范围统一，记为ψ₁—ψ₇，定义为：

其中q＝0,1,…,7。

c)提取下冲程部分的曲线矩特征向量，划分后的下冲程部分曲线由N-l个采样点组成，每个采样点的坐标为(x_j,y_j)，其s+r阶曲线矩m'_sr定义为：

$m_{\mathrm{sr}}^{\′}=\underset{j=l+1}{\overset{N}{\mathrm{\Σ}}}{x_j}^s{y_j}^r\mathrm{\Δ}{z}_j---\left(12\right)$

其中j＝l+1,l+2,…,N；s,r＝0,1,2；Δz_j为两个连续的采样点间的距离，

μ'_sr为s+r阶中心矩，定义为：

${\mathrm{\μ}}_{\mathrm{sr}}^{\′}=\underset{j=l+1}{\overset{N}{\mathrm{\Σ}}}{(x_j-\stackrel{\‾}{x^{\′}})}^s{(y_j-\stackrel{\‾}{y^{\′}})}^r\mathrm{\Δ}{z}_j---\left(13\right)$

其中 $\stackrel{\‾}{x^{\′}}=\frac{m_{10}}{m_{00}},\stackrel{\‾}{y^{\′}}=\frac{m_{01}}{m_{00}};$ 点为曲线的重心坐标；

计算出如下各阶中心矩：μ′₀₀＝m′₀₀，μ′₁₀＝0，μ′₀₁＝0， ${\mathrm{\μ}}_{11}^{\′}=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}^{\′}.$

对提取的各阶中心矩进行规范化处理，由η'_sr表示，定义为：

η'_sr＝μ'_sr/(μ'₀₀)^s+r+1(14)

那么构造的下冲程部分的7个曲线矩特征向量分别定义如下：

采用如下修正公式使构造的下冲程部分的7个曲线矩特征向量的取值范围统一，记为ψ₈—ψ₁₄，定义为：

其中q’＝8,9,…,14。

将提取得到的能够表征井下泵功图图形特征的14个特征向量ψ₁—ψ₁₄，作为下一步预测的变量。对图4所示的划分后的井下泵功图进行曲线矩特征向量的提取，得到的14个特征向量如表1所示。

表1提取的泵功图曲线矩特征向量

上冲程部分	ψ1	ψ2	ψ3	ψ4	ψ5	ψ6	ψ7
								曲线矩特征向量值	1.7060	0.2577	1.6842	0.7475	1.8887	0.1878	1.6955
下冲程部分	ψ8	ψ9	ψ10	ψ11	ψ12	ψ13	ψ14
								曲线矩特征向量值	1.3639	0.4036	2.2777	1.3236	3.1241	1.1877	1.3260

4、根据多变量灰色模型MGM(1,n)的基本理论建立有杆泵抽油井井下故障的多变量灰色预测模型。

所建立的泵功图曲线矩特征向量灰时间序列如下：

${\mathrm{\ψ}}^{\left(0\right)}\left(k\right)=\left[\begin{array}{cccccc}{\mathrm{\ψ}}_1^{\left(0\right)}\left(1\right)& {\mathrm{\ψ}}_1^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_1^{\left(0\right)}\left(k\right)& ...& {\mathrm{\ψ}}_1^{\left(0\right)}\left(m\right)\\ {\mathrm{\ψ}}_2^{\left(0\right)}\left(1\right)& {\mathrm{\ψ}}_2^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_2^{\left(0\right)}\left(k\right)& ...& {\mathrm{\ψ}}_2^{\left(0\right)}\left(m\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ {\mathrm{\ψ}}_p^{\left(0\right)}\left(1\right)& {\mathrm{\ψ}}_p^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_p^{\left(0\right)}\left(k\right)& ...& {\mathrm{\ψ}}_p^{\left(0\right)}\left(m\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(1\right)& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(k\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(m\right)\end{array}\right]---\left(23\right)$

其中p＝1,2,…,14；k＝1,2,…,m，m为所采用的泵功图的数量；表示泵功图曲线矩特征向量灰时间序列中的第p个变量中的第k个影响因素；

令为相应的一次累加生成序列，有：

${\mathrm{\ψ}}_p^{\left(1\right)}\left(k\right)=\underset{g=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\ψ}}_p^{\left(0\right)}\left(g\right)---\left(24\right)$

建立14元一阶常微分方程，有：

$\left\{\begin{array}{c}\frac{d{\mathrm{\ψ}}_1^{\left(1\right)}}{\mathrm{dt}}=a_{1-1}{\mathrm{\ψ}}_1^{\left(1\right)}+a_{1-2}{\mathrm{\ψ}}_2^{\left(1\right)}+...+a_{1-14}{\mathrm{\ψ}}_{14}^{\left(1\right)}+b_1\\ \frac{d{\mathrm{\ψ}}_2^{\left(1\right)}}{\mathrm{dt}}=a_{2-1}{\mathrm{\ψ}}_1^{\left(1\right)}+a_{2-2}{\mathrm{\ψ}}_2^{\left(1\right)}+...+a_{2-14}{\mathrm{\ψ}}_{14}^{\left(1\right)}+b_2\\ .\\ .\\ .\\ \frac{d{\mathrm{\ψ}}_{14}^{\left(1\right)}}{\mathrm{dt}}=a_{14-1}{\mathrm{\ψ}}_1^{\left(1\right)}+a_{14-2}{\mathrm{\ψ}}_2^{\left(1\right)}+...+a_{14-14}{\mathrm{\ψ}}_{14}^{\left(1\right)}+b_{14}\end{array}\right.---\left(25\right)$

记

$A=\left[\begin{array}{cccc}{a}_{1-1}& a_{1-2}& ...& a_{1-14}\\ a_{2-1}& a_{2-2}& ...& a_{2-14}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{14-1}& a_{14-2}& ...& a_{14-14}\end{array}\right]---\left(26\right)$

B＝(b₁,b₁,…,b_n)^T(27)

A和B为辨识参数，那么公式(25)可以简记为

$\frac{d{\mathrm{\Ψ}}^{\left(1\right)}}{\mathrm{dt}}=A{\mathrm{\Ψ}}^{\left(1\right)}+B---\left(28\right)$

其中 ${\mathrm{\Ψ}}^{\left(1\right)}={({\mathrm{\ψ}}_1^{\left(1\right)},{\mathrm{\ψ}}_2^{\left(1\right)},...,{\mathrm{\ψ}}_{14}^{\left(1\right)})}^T;$

在[0,t]区间上积分，得到连续时间响应为

${\mathrm{\Ψ}}^{\left(1\right)}\left(t\right)=e^{\mathrm{At}}{\mathrm{\Ψ}}^{\left(1\right)}\left(0\right)+A^{-1}(e^{\mathrm{At}}-I)B---\left(29\right)$

其中 $e^{\mathrm{At}}=I+\underset{u=1}{\overset{\∞}{\mathrm{\Σ}}}\frac{A^u}{u!}{t}^u;$

将公式(25)进行离散化，得到：

${\mathrm{\ψ}}_p^{\left(0\right)}\left(k\right)=\underset{w=1}{\overset{14}{\mathrm{\Σ}}}\frac{a_{p-w}}{2}({\mathrm{\ψ}}_w^{\left(1\right)}\left(k\right)+{\mathrm{\ψ}}_w^{\left(1\right)}(k-1))---\left(30\right)$

其中p＝1,2,…,14；k＝2,3,…,m；

记

$E=\left[\begin{array}{cccc}{a}_{1-1}& a_{2-1}& ...& a_{14-1}\\ a_{1-2}& a_{2-2}& ...& a_{14-2}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{1-14}& a_{2-14}& ...& a_{14-14}\\ b_1& b_2& ...& b_{14}\end{array}\right]---\left(31\right)$

由最小二乘法得到E的辨识值E，有：

E＝(L^TL)^-1L^TY(32)

其中，

$Y=\left[\begin{array}{cccc}{\mathrm{\ψ}}_1^{\left(0\right)}\left(2\right)& {\mathrm{\ψ}}_2^{\left(0\right)}\left(2\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(2\right)\\ {\mathrm{\ψ}}_1^{\left(0\right)}\left(3\right)& {\mathrm{\ψ}}_2^{\left(0\right)}\left(3\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(3\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\ψ}}_1^{\left(0\right)}\left(m\right)& {\mathrm{\ψ}}_2^{\left(0\right)}\left(m\right)& ...& {\mathrm{\ψ}}_{14}^{\left(0\right)}\left(m\right)\end{array}\right]---\left(34\right)$

$E={\left[\begin{array}{cc}A& B\end{array}\right]}^T=\left[\begin{array}{cccc}{a}_{1-1}& a_{2-1}& ...& a_{14-1}\\ a_{1-2}& a_{2-2}& ...& a_{14-2}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{1-14}& a_{2-14}& ...& a_{14-14}\\ {\tilde{b}}_1& {\tilde{b}}_2& ...& {\tilde{b}}_{14}\end{array}\right]---\left(35\right)$

由公式(35)得到辨识参数A和B的辨识值A和B，分别为：

$A=\left[\begin{array}{cccc}{a}_{1-1}& a_{1-2}& ...& a_{1-14}\\ a_{2-1}& a_{2-2}& ...& a_{2-14}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{14-1}& a_{14-2}& ...& a_{14-14}\end{array}\right]---\left(36\right)$

$B=\left[\begin{array}{c}{\tilde{b}}_1\\ {\tilde{b}}_2\\ .\\ .\\ .\\ {\tilde{b}}_{14}\end{array}\right]---\left(37\right)$

根据得到的辨识值A和B，可以将公式(29)写成如下离散形式，

Ψ⁽⁰⁾(1)＝Ψ⁽⁰⁾(1)(38)

Ψ⁽¹⁾(k)＝e^A(k-1)Ψ⁽¹⁾(1)+A^-1(e^A(k)-I)B(39)

其中Ψ⁽¹⁾(k)和Ψ⁽¹⁾(1)，将公式(36)和公式(38)代入公式(30)求得；

则得到预测模型为：

Ψ⁽⁰⁾(k)＝Ψ⁽¹⁾(k)-Ψ⁽¹⁾(k-1)(40)

其中k＝2,3,…,m。

对于有杆泵抽油井井下故障的多变量灰色预测模型中所采用的一段连续时间内的样本，连续时间的单位可以为小时或天，所述样本的数量为50－200个，以提高预测模型的计算效率和准确度。本实施例采集实际生产中某口有杆泵抽油井在50天内的50幅示功图作为样本，首先将地面示功图转化为井下泵功图；然后分别提取50幅井下泵功图的14个曲线矩特征作为特征向量；再根据公式(23)—公式(40)进行预测模型的计算，得到的结果如图5－图18所示。由图5－图18，采用多变量灰色模型分别对提取的井下泵功图的14个曲线矩特征向量进行预测，预测值都能达到满意的效果。

采用上述50个样本通过多变量灰色模型预测第51个样本，得到其14个曲线矩特征向量值如表2所示。

表2预测样本的曲线矩特征向量

上冲程部分	ψ1	ψ2	ψ3	ψ4	ψ5	ψ6	ψ7
								曲线矩特征向量值	1.1870	0.5246	1.4833	0.9148	1.0945	0.3579	1.0045
下冲程部分	ψ8	ψ9	ψ10	ψ11	ψ12	ψ13	ψ14
								曲线矩特征向量值	0.9699	0.6187	3.0128	1.7286	1.9872	1.5046	0.8967

5、最后，判断所预测的样本属于哪类故障类型。假设有杆泵抽油井井下故障的类型由F_M表示，其中M＝1,2,3,…。可以通过已有的示功图信息建立每一种故障类型的标准集，已有的信息由已知道故障类型的示功图得到，用来建立每一种故障类型的标准集的示功图的数量大于2个。即对已确定故障类型的一定数量的示功图通过将地面示功图转化为井下泵功图和分别提取其14个曲线矩特征向量，得到每一种故障类型中每一个曲线矩特征向量的取值区间。所得到的取值区间由表示，p＝1,2,…,14，其中表示第F_M种故障类型的第p个曲线矩特征向量的下限值，表示第F_M种故障类型的第p个曲线矩特征向量的上限值。由已有的信息建立每一种故障类型的标准集。

本发明中有杆泵抽油井井下故障类型有以下9种：“正常”、“气体影响”、“供液不足”、“抽油杆断落”、“游动阀漏失”、“固定阀漏失”、“泵上碰”、“泵下碰”、“油井出砂”。所建立的各故障类型的标准集如表3所示。

表3各故障类型的标准集

ψ1

ψ2

ψ3

ψ4

ψ5

ψ6

ψ7

正常

[0.05,0.45]

[0.87,2.45]

[0.76,2.49]

[3.14,5.78]

[0.03,0.17]

[0.66,2.06]

[0.66,2.24]

气体影响

[0.06,0.58]

[1.13,1.44]

[1.23,1.94]

[2.42,2.71]

[0.12,1.08]

[1.04,2.19]

[1.12,3.16]

供液不足

[0.22,0.36]

[0.79,2.88]

[0.8,3.16]

[3.57,5.91]

[1.31,1.73]

[2.8,5.58]

[1.47,2.86]

抽油杆断落

[0.27,0.39]

[0.7,0.91]

[0.78,0.94]

[3.1,3.35]

[0.04,0.11]

[1.97,2.65]

[3.75,5.19]

游动阀漏失

[0.05,0.68]

[1.14,3.02]

[1.17,3.26]

[3.55,5.65]

[0.09,0.2]

[0.64,1.1]

[0.67,1.26]

固定阀漏失

[0.06,0.64]

[1.16,3.08]

[1.16,3.45]

[0.04,0.17]

[0.64,1.77]

[1.27,3.33]

[3.96,5.68]

泵上碰

[0.94,3.11]

[0.94,3.98]

[0.94,4.19]

[0.03,0.47]

[0.74,2.3]

[0.76,2.41]

[2.97,5.34]

泵下碰

[0.02,0.06]

[0.74,1.28]

[1.47,2.56]

[0.77,1.28]

[3.65,5.51]

[0.73,1.99]

[1.46,4.07]

油井出砂

[0.08,0.45]

[0.29,1.14]

[0.3,1.16]

[0.02,0.14]

[0.34,2.3]

[0.36,2.34]

[0.53,1.92]

ψ8

ψ9

ψ10

ψ11

ψ12

ψ13

ψ14

正常

[3.91,5.66]

[0.01,1.09]

[0.52,2.62]

[0.53,2.69]

[0.65,1.03]

[1.31,2.11]

[4.5,5.06]

气体影响

[2.1,3.98]

[0.56,0.9]

[2.56,3.99]

[2.41,2.83]

[0.47,0.67]

[1.5,2.42]

[1.99,2.49]

供液不足

[0.01,0.19]

[0.67,2.56]

[0.7,2.74]

[0.66,3.01]

[1.35,5.93]

[0.71,2.97]

[0.02,0.18]

抽油杆断落

[4.16,5.93]

[0.22,0.38]

[0.55,0.89]

[0.58,0.91]

[0.08,0.09]

[1.69,1.99]

[0.66,1.07]

游动阀漏失

[1.5,2.67]

[3.02,8.12]

[1.52,3.98]

[0,0.23]

[0.64,2.46]

[0.64,2.6]

[2.93,4.94]

固定阀漏失

[0.09,0.84]

[0.48,1.55]

[0.49,1.56]

[0.07,0.17]

[1.29,2.64]

[1.3,2.68]

[3.82,6.37]

泵上碰

[1.38,2.77]

[2.77,5.66]

[5.87,7.29]

[0.64,1.16]

[1.3,2.33]

[0.65,1.17]

[4.28,5.65]

泵下碰

[4.58,5.81]

[3.17,4.33]

[1.27,3.89]

[1.27,4.09]

[2.54,8.08]

[1.27,4.21]

[6.8,9.38]

油井出砂

[1.06,3.97]

[0.53,2.08]

[0.4,1.76]

[0.41,1.8]

[0.81,3.57]

[0.41,1.81]

[0.53,2]

多变量灰色模型所预测的第51个样本的曲线矩特征向量值由表示，那么预测样本的每一个曲线矩特征向量值与每一种故障类型中相应的曲线矩特征向量值的距离可以由下式进行计算，

$d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})=|U_{{\mathrm{\ψ}}_p}-\frac{V_{F_M\mathrm{pa}}+V_{F_M\mathrm{pb}}}{2}|+\frac{V_{F_M\mathrm{pb}}-V_{F_M\mathrm{pa}}}{2}---\left(41\right)$

然后由下式计算预测样本与每一种故障类型的灰关联系数，

${\mathrm{\ξ}}_{p{F}_M}=\frac{\underset{M}{\min }\underset{p}{\min }\left\{d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})\right\}+\mathrm{\ρ}\underset{M}{\max }\underset{p}{\max }\left\{d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})\right\}}{d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})+\mathrm{\ρ}\underset{M}{\max }\underset{p}{\max }\left\{d(U_{{\mathrm{\ψ}}_p},V_{F_{M}p})\right\}}---\left(42\right)$

其中p＝1,2,…,14；M＝1,2,3,…。ρ∈[0,1]为分辨系数，在这里取0.5；

再由下式计算预测样本与每一种故障类型的灰关联度，

$\mathrm{Gre}{y}_{F_M}=\frac{1}{p}\underset{p=1}{\overset{14}{\mathrm{\Σ}}}{\mathrm{\ξ}}_{p{F}_M}---\left(43\right)$

根据公式(41)—公式(43)，计算预测样本与每一种故障类型的灰关联度，得到的结果如表4所示。

表4预测样本与各故障类型的灰关联度

由表4，预测样本与“气体影响”故障类型的灰关联度最大，认为属于“气体影响”故障类型。

以上仅为本发明的具体实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于多变量灰色模型的有杆泵抽油井井下故障预测方法，其特征是：

1.1)采用示功图远程数据采集***，首先，将无线示功仪安装在有杆泵抽油井的抽油机的驴头上，采集光杆的载荷和位移数据；无线RTU以无线的方式接收无线示功仪采集的数据，由数据采集模块通过无线AP从站模块以无线网络远程传送到无线AP主站模块上；示功图监测及数据处理服务器接收采集到的光杆的载荷和位移数据，绘制采集到的地面示功图图形；

1.2)将地面示功图转化为井下泵功图，以得到能够真实反映抽油泵工作状况的示功图图形；

1.3)提取井下泵功图的曲线矩特征作为特征向量，将井下泵功图划分为上冲程部分和下冲程部分，然后根据曲线矩理论分别提取上冲程部分的7个曲线矩特征向量和下冲程部分的7个曲线矩特征向量，将得到的14个曲线矩特征向量作为下一步预测的变量，具体步骤如下：

a)将井下泵功图划分为上冲程部分和下冲程部分，其中上冲程部分曲线反映了有杆泵抽油井的上冲程过程，是抽油杆柱在动力装置的作用下向上拉动柱塞，抽油泵内吸入油液而井口排出油液的过程；下冲程部分曲线反映了有杆泵抽油井的下冲程过程，是抽油杆柱在动力装置的作用下向下推动柱塞，抽油泵向油管内排出油液的过程；假设井下泵功图由N个采样点组成，那么从第1个点到位移最大的点(假设为第l个点)为上冲程部分；从第l+1个点到第N个点为下冲程部分；

b)划分后的上冲程部分曲线由l个采样点组成，每个采样点的坐标为(x_i,y_i)，其s+r阶曲线矩m_sr定义为：

其中i＝1,…,l；s,r＝0,1,2；Δz_i为两个连续的采样点间的距离，

μ_sr为s+r阶中心矩，定义为：

其中点为曲线的重心坐标；

计算出如下各阶中心矩：μ₀₀＝m₀₀，μ₁₀＝0，μ₀₁＝0，

对提取的各阶中心矩进行规范化处理，由η_sr表示，定义为：

η_sr＝μ_sr/(μ₀₀)^s+r+1(3)

那么构造的上冲程部分的7个曲线矩特征向量分别定义如下：

采用如下修正公式使构造的上冲程部分的7个曲线矩特征向量的取值范围统一，记为ψ₁—ψ₇，定义为：

其中q＝0,1,…,7；

c)提取下冲程部分的曲线矩特征向量，划分后的下冲程部分曲线由N-l个采样点组成，每个采样点的坐标为(x_j,y_j)，其s+r阶曲线矩m'_sr定义为：

其中j＝l+1,l+2,…,N；s,r＝0,1,2；Δz_j为两个连续的采样点间的距离，

μ'_sr为s+r阶中心矩，定义为：

其中点为曲线的重心坐标；

计算出如下各阶中心矩：μ'₀₀＝m'₀₀，μ′₁₀＝0，μ'₀₁＝0，

对提取的各阶中心矩进行规范化处理，由η'_sr表示，定义为：

η'_sr＝μ'_sr/(μ'₀₀)^s+r+1(14)

那么构造的下冲程部分的7个曲线矩特征向量分别定义如下：

采用如下修正公式使构造的下冲程部分的7个曲线矩特征向量的取值范围统一，记为ψ₈—ψ₁₄，定义为：

其中q’＝8,9,…,14；

将提取得到的能够表征井下泵功图图形特征的14个特征向量ψ₁—ψ₁₄，作为下一步预测的变量；

1.4)根据多变量灰色模型MGM(1,n)的基本理论建立有杆泵抽油井井下故障的多变量灰色预测模型，利用MGM(1,n)模型建立多个变量之间的动态相关性，弥补单变量方法的局限性；将一段连续时间内的泵功图所提取的14个曲线矩特征向量作为多变量灰色模型(MGM(1,14))的输入，建立灰时间序列，每一个曲线矩特征向量作为灰时间序列的一个变量；由所建立的多变量灰色预测模型分别预测14个曲线矩特征向量的下一个时间点的特征向量，所述14个预测的曲线矩特征向量即为能够表征预测样本特征的特征向量；

将所建立的泵功图曲线矩特征向量灰时间序列表示如下：

其中p＝1,2,…,14；k＝1,2,…,m，m为所采用的泵功图的数量；表示泵功图曲线矩特征向量灰时间序列中的第p个变量中的第k个影响因素；

令为相应的一次累加生成序列，有：

建立14元一阶常微分方程，有：

记

B＝(b₁,b₁,…,b_n)^T(27)

A和B为辨识参数，那么将公式(25)简记为

其中

在[0,t]区间上积分，得到连续时间响应为

Ψ⁽¹⁾(t)＝e^AtΨ⁽¹⁾(0)+A^-1(e^At-I)B(29)

其中

将公式(25)进行离散化，得到：

其中p＝1,2,…,14；k＝2,3,…,m；

记

由最小二乘法得到E的辨识值E，有：

E＝(L^TL)^-1L^TY(32)

其中，

由公式(35)得到辨识参数A和B的辨识值A和B，分别为：

根据得到的辨识值A和B，将公式(29)写成如下离散形式，

Ψ⁽⁰⁾(1)＝Ψ⁽⁰⁾(1)(38)

Ψ⁽¹⁾(k)＝e^A(k-1)Ψ⁽¹⁾(1)+A^-1(e^A(k)-I)B(39)

其中Ψ⁽¹⁾(k)和Ψ⁽¹⁾(1)可以将公式(36)和公式(38)代入公式(30)求得；

得到预测模型为：

Ψ⁽⁰⁾(k)＝Ψ⁽¹⁾(k)-Ψ⁽¹⁾(k-1)(40)

其中k＝2,3,…,m；

1.5)最后，判断预测样本属于哪类故障类型通过已有的示功图信息建立每一种故障类型的标准集，即对已确定故障类型的一定数量的示功图通过将地面示功图转化为井下泵功图和分别提取其14个曲线矩特征向量，得到每一种故障类型中每一个曲线矩特征向量的取值区间；所得到的取值区间由表示，p＝1,2,…,14，其中表示第F_M种故障类型的第p个曲线矩特征向量的下限值，表示第F_M种故障类型的第p个曲线矩特征向量的上限值；

计算预测样本的每一个曲线矩特征向量值与每一种故障类型中相应的曲线矩特征向量值的距离，由下式进行计算：

然后由下式计算预测样本与每一种故障类型的灰关联系数，

其中p＝1,2,…,14；M＝1,2,3,…。ρ∈[0,1]为分辨系数，在这里取0.5；

再由下式计算预测样本与每一种故障类型的灰关联度，

最大的灰关联度所对应的故障类型即为预测样本的故障类型。

2.根据权利要求1所述的基于多变量灰色模型的有杆泵抽油井井下故障预测方法，其特征是：步骤3)中井下泵功图的上冲程部分为有杆泵抽油井抽油杆的悬点由最低点运动到最高点的部分；下冲程部分为有杆泵抽油井抽油杆的悬点由最高点运动到最低点的部分。

3.根据权利要求1所述的基于多变量灰色模型的有杆泵抽油井井下故障预测方法，其特征是：步骤4)对于有杆泵抽油井井下故障的多变量灰色预测模型中所采用的一段连续时间内的样本，连续时间的单位可以为小时或天，所述样本的数量为50－200个，以提高预测模型的计算效率和准确度。

4.根据权利要求1所述的基于多变量灰色模型的有杆泵抽油井井下故障预测方法，其特征是：步骤5)由已有的信息建立每一种故障类型的标准集，已有的信息由已知道故障类型的示功图得到，用来建立每一种故障类型的标准集的示功图的数量大于2个。

5.根据权利要求1所述的基于多变量灰色模型的有杆泵抽油井井下故障预测方法，其特征是：所述有杆泵抽油井井下故障类型有以下9种：“正常”、“气体影响”、“供液不足”、“抽油杆断落”、“游动阀漏失”、“固定阀漏失”、“泵上碰”、“泵下碰”、“油井出砂”。