CN114510880A - 一种基于傅里叶变换和几何特征的有杆泵工况诊断方法 - Google Patents

一种基于傅里叶变换和几何特征的有杆泵工况诊断方法 Download PDF

Info

Publication number
CN114510880A
CN114510880A CN202210407067.1A CN202210407067A CN114510880A CN 114510880 A CN114510880 A CN 114510880A CN 202210407067 A CN202210407067 A CN 202210407067A CN 114510880 A CN114510880 A CN 114510880A
Authority
CN
China
Prior art keywords
sucker
rod pump
indicator diagram
load
follows
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210407067.1A
Other languages
English (en)
Other versions
CN114510880B (zh
Inventor
张凯
李敏
张黎明
张华清
严侠
刘丕养
杨勇飞
孙海
姚军
樊灵
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN202210407067.1A priority Critical patent/CN114510880B/zh
Publication of CN114510880A publication Critical patent/CN114510880A/zh
Application granted granted Critical
Publication of CN114510880B publication Critical patent/CN114510880B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • E—FIXED CONSTRUCTIONS
    • E21—EARTH OR ROCK DRILLING; MINING
    • E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B47/00—Survey of boreholes or wells
    • E21B47/008—Monitoring of down-hole pump systems, e.g. for the detection of "pumped-off" conditions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/243—Classification techniques relating to the number of classes
    • G06F18/24323—Tree-organised classifiers
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G06F18/253—Fusion techniques of extracted features

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Mining & Mineral Resources (AREA)
  • Geology (AREA)
  • Fluid Mechanics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geophysics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geochemistry & Mineralogy (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Complex Calculations (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明公开了一种基于傅里叶变换和几何特征的有杆泵工况诊断方法，属于有杆泵工况诊断技术领域，包括以下步骤：对所获有杆泵生产数据进行选择；针对示功图进行基于波动方程的傅里叶系数提取；获取示功图曲线数据，进行示功图简单几何特征提取；采用DCA将傅里叶系数与简单几何特征进行融合；使用XGBoost算法建立有杆泵工况诊断模型，并进行模型训练；进行有杆泵工况诊断模型的参数优化；对优化后的有杆泵工况诊断模型，进行模型性能评价；将训练完成的有杆泵工况诊断模型应用到油田现场。本发明能够提高油田开发现场对有杆泵工况诊断的效率，同时提高对油田现有数据的有效利用，实现高效诊断有杆泵工况。

Description

一种基于傅里叶变换和几何特征的有杆泵工况诊断方法
技术领域
本发明属于有杆泵工况诊断技术领域，具体涉及一种基于傅里叶变换和几何特征的有杆泵工况诊断方法。
背景技术
在石油行业，有杆泵是最常用的人工举升方式，而有杆抽油机大多是在野外工作，地理位置偏僻，自然环境恶劣，而且有杆抽油系统的井下工况复杂，如果有杆泵的故障诊断及处理措施不当，将会直接影响油田的产量，严重时甚至会发生危及人身安全的生产事故。为减少油田开发风险，降低油田开发成本，使用机器学习的方法，实现对有杆泵工况的智能诊断，该技术同时也是实现油井智能管理的有效技术。
通过机器学习的方法基于示功图进行有杆泵工况诊断的方法已经被油田使用，近年来，机器学习在有杆泵工况诊断中的应用越来越多，例如基于改进的超球支持向量机故障诊断模型，基于极限学习机的有杆泵工况诊断。
但是在示功图特征提取方面，常用的方法有：示功图Hu矩、傅里叶描述子、示功图曲率特征向量，以及灰度矩阵，还不能够完全满足有杆泵机器学习诊断的需要，导致有杆泵工况诊断的效率低，准确性不高。
发明内容
为了解决上述问题，本发明提出了一种基于傅里叶变换和几何特征的有杆泵工况诊断方法，提高油田开发现场对有杆泵工况诊断的效率，同时提高对油田现有数据的有效利用，实现高效诊断有杆泵工况。
本发明的技术方案如下：
一种基于傅里叶变换和几何特征的有杆泵工况诊断方法，包括以下步骤：
步骤1、对所获有杆泵生产数据进行选择；
步骤2、针对示功图进行基于波动方程的傅里叶系数提取；
步骤3、获取示功图曲线数据，进行示功图简单几何特征提取；
步骤4、采用DCA将傅里叶系数与简单几何特征进行融合；
步骤5、使用XGBoost算法建立有杆泵工况诊断模型，并进行模型训练；
步骤6、进行有杆泵工况诊断模型的参数优化；
步骤7、对优化后的有杆泵工况诊断模型，进行模型性能评价；
步骤8、将训练完成的有杆泵工况诊断模型应用到油田现场，实时采集油田现场数据，进行现场油井有杆泵工作状况的诊断。
进一步地，步骤1的具体内容为：
步骤1.1、选择出有杆泵工作过程中各周期中的悬点冲程值、载荷值、理论上载荷值、理论下载荷值；
步骤1.2、通过人工方法清洗异常的示功图；若示功图内包括数据越界、离散点、曲线不闭合、曲线反向或空数据的情况则被判定为异常的示功图。
进一步地，步骤2的具体内容为：
步骤2.1、将有杆泵冲程数据与载荷数据提取出来，分别记为：U={S|S h ∈well h }、D={W|W h ∈well h }；式中，S为所有井有杆泵示功图所有点的冲程取值集合，h为第h口井；S h 为第h口井有杆泵示功图所有点的冲程取值集合；W为所有井有杆泵示功图所有点的载荷取值集合；W h 为第h口井有杆泵示功图所有点的载荷取值集合；well h 为第h口井的所有数据信息；
步骤2.2、以第h口井的计算为例，采用吉布斯求解抽油杆运动方程，求解时提出的悬点冲程、载荷的方程见式（1）至（2）：
Figure 100002_DEST_PATH_IMAGE001
(1)
Figure 100002_DEST_PATH_IMAGE002
(2)
其中，t为时间，n取值为1到
Figure 100002_DEST_PATH_IMAGE003
，是累加计算中的遍历变量，
Figure 100002_DEST_PATH_IMAGE004
为系数
Figure 100002_DEST_PATH_IMAGE005
在级数 为0时的取值，
Figure 100002_DEST_PATH_IMAGE006
为系数
Figure 100002_DEST_PATH_IMAGE007
在级数为0时的取值，
Figure 100002_DEST_PATH_IMAGE008
为所确定的傅里叶级数，
Figure 100002_DEST_PATH_IMAGE009
为运动 角速度；
其中，傅里叶系数
Figure 100002_DEST_PATH_IMAGE010
、
Figure 100002_DEST_PATH_IMAGE011
、
Figure 100002_DEST_PATH_IMAGE012
、
Figure 100002_DEST_PATH_IMAGE013
的计算方法见式（3）至（6）；
Figure 100002_DEST_PATH_IMAGE014
(3)
Figure 100002_DEST_PATH_IMAGE015
(4)
Figure 100002_DEST_PATH_IMAGE016
(5)
Figure 100002_DEST_PATH_IMAGE017
(6)
式中，M为示功图所具数据点数量，i为累加和中的遍历变量，取值为0到M-1。
进一步地，步骤3中，简单几何特征的提取包括提取示功图对角线AC、BD斜率，增载 线AB边斜率，卸载线CD边斜率，对角线AC、BD的长度，上静载线BC边长度，下静载线DA边长 度，上静载平均载荷与理论上载荷的差值C BC ，下静载平均载荷与理论下载荷的差值C DA ，以 及示功图最大载荷与最小载荷之间的差值
Figure 100002_DEST_PATH_IMAGE018
；提取的具体过程为：
步骤3.1、提取简单几何特征进行工况的判断，每种简单几何特征都对应一定的工况判别情况；
步骤3.2、示功图A、B、C、D四点冲程载荷值根据有杆泵理论示功图的曲线特点进行计算；当A点冲程最小时，其坐标U a 为最小冲程，D a 为最小冲程对应的载荷；当B点载荷最大时，其坐标D b 为最大载荷，U b 为最大载荷对应的冲程；当C点冲成最大时，其坐标U c 为最大冲程，D c 为最大冲程对应的载荷；当D点载荷最小时，其坐标D d 为最小载荷，U d 为最小载荷对应的冲程；
步骤3.3、以第h口井的计算为例，示功图简单几何特征提取的计算方法见式（7）至（17）：
Figure 100002_DEST_PATH_IMAGE019
(7)
式中，K hAC 为第h口井有杆泵示功图对角线AC的斜率；D ha 为第h口井示功图A点的载荷值；U ha 为第h口井示功图A点的冲程值；D hc 为第h口井示功图C点的载荷值；U hc 为第h口井示功图C点的冲程值；
Figure 100002_DEST_PATH_IMAGE020
(8)
式中，K hBD 为第h口井有杆泵示功图对角线BD的斜率；D hb 为第h口井示功图B点的载荷值；U hb 为第h口井示功图B点的冲程值；D hd 为第h口井示功图D点的载荷值；U hd 为第h口井示功图D点的冲程值；
Figure 100002_DEST_PATH_IMAGE021
(9)
式中，K hAB 为第h口井有杆泵示功图加载线AB线段的斜率；
Figure 100002_DEST_PATH_IMAGE022
(10)
式中，K hCD 为第h口井有杆泵示功图卸载线CD线段的斜率；
Figure 100002_DEST_PATH_IMAGE023
(11)
式中，L hAC 为第h口井有杆泵示功图对角线AC的长度；
Figure 100002_DEST_PATH_IMAGE024
(12)
式中，L hBD 为第h口井有杆泵示功图对角线BD的长度；
Figure 100002_DEST_PATH_IMAGE025
(13)
式中，L hBC 为第h口井有杆泵示功图上冲程BC线段的长度；
Figure 100002_DEST_PATH_IMAGE026
(14)
式中，L hDA 为第h口井有杆泵示功图下冲程DA线段的长度；
Figure 100002_DEST_PATH_IMAGE027
(15)
式中，C hBC 为第h口井示功图B、C两点的平均载荷与理论上载荷的差值，D h理论上为第h口井示功图理论上载荷值；
Figure 100002_DEST_PATH_IMAGE028
(16)
式中，C hDA 为第h口井示功图D、A两点的平均载荷与理论下载荷的差值，D h理论下为第h口井示功图理论下载荷值；
Figure 100002_DEST_PATH_IMAGE029
(17)
式中，
Figure 100002_DEST_PATH_IMAGE030
为第h口井示功图最大载荷与最小载荷的差值；D hmax 为第h口井示功图 最大载荷值；D hmin 为第h口井示功图最小载荷值。
进一步地，步骤4中，DAC特征融合的具体步骤为：
步骤4.1、用F 1表示步骤2中计算出的所有井的傅里叶系数，F 1={F 11，F 21，...，F h1，...}
其中
Figure 100002_DEST_PATH_IMAGE031
，用F 2表 示步骤3中计算出的所有井的简单几何特征，F 2={F 12，F 22，...，F h2，...}，其中
Figure 100002_DEST_PATH_IMAGE032
；
步骤4.2、假设数据矩阵中的样本是从c个单独的类中收集的；相应地，数据矩阵的N列被分成c个单独的组，其中N v 列属于第v类；分别求出F 1，F 2的协方差，计算方法见式（18）与式（19）：
Figure 100002_DEST_PATH_IMAGE033
(18)
Figure 100002_DEST_PATH_IMAGE034
(19)
式中，
Figure 100002_DEST_PATH_IMAGE035
为F 1的协方差，
Figure 100002_DEST_PATH_IMAGE036
为F 2的协方差，p为F 1的维数，q为F 2的维数；
Figure 100002_DEST_PATH_IMAGE037
表示第v类工况的傅里叶系数所构成向量的均值；
Figure 100002_DEST_PATH_IMAGE038
表示所有井所有类工况傅里叶系数 构成向量的均值；
Figure 100002_DEST_PATH_IMAGE039
表示第v类工况的简单几何特征所构成向量的均值；
Figure 100002_DEST_PATH_IMAGE040
表示所有井所 有类工况简单几何特征构成向量的均值；
Figure 100002_DEST_PATH_IMAGE041
、
Figure 100002_DEST_PATH_IMAGE042
分别通过
Figure 100002_DEST_PATH_IMAGE043
、
Figure 100002_DEST_PATH_IMAGE044
计算获得；
步骤4.3、利用协方差分别计算出F 1，F 2变换后的矩阵
Figure 100002_DEST_PATH_IMAGE045
，
Figure 100002_DEST_PATH_IMAGE046
，以
Figure 976831DEST_PATH_IMAGE045
计算为例，计算 方法见式（20）至式（24）：
Figure 100002_DEST_PATH_IMAGE047
(20)
式中，P是正交矩阵；
Figure 100002_DEST_PATH_IMAGE048
是按降序排列的实非负特征值的对角矩阵；
Figure 100002_DEST_PATH_IMAGE049
(21)
式中，Q由P中前r个特征向量组成，对应于矩阵P中最大的r个非零特征值；
Figure 100002_DEST_PATH_IMAGE050
(22)
记：
Figure 100002_DEST_PATH_IMAGE051
；
Figure 100002_DEST_PATH_IMAGE052
(23)
Figure 100002_DEST_PATH_IMAGE053
(24)
其中，
Figure 891829DEST_PATH_IMAGE045
为F 1变换后的矩阵；I为单位阵；
同理，得到F 2变换后的矩阵
Figure 805558DEST_PATH_IMAGE046
；
得到新傅里叶系数
Figure 611840DEST_PATH_IMAGE045
和新简单几何特征
Figure 661836DEST_PATH_IMAGE046
后，得到
Figure 100002_DEST_PATH_IMAGE054
， 对
Figure 100002_DEST_PATH_IMAGE055
进行SVD分解，得到式（25），并根据式（26）计算出融合空间内特征
Figure 100002_DEST_PATH_IMAGE056
，
Figure 100002_DEST_PATH_IMAGE057
；下面以
Figure 403527DEST_PATH_IMAGE056
计算为例；
Figure 100002_DEST_PATH_IMAGE058
(25)
其中，H为r×r的左酉矩阵，V为r×r的右酉矩阵；J是一个对角元素非零的对角矩阵；
令
Figure 100002_DEST_PATH_IMAGE059
，
Figure 100002_DEST_PATH_IMAGE060
，则
Figure 100002_DEST_PATH_IMAGE061
(26)
同理，得到融合空间内的特征
Figure 100002_DEST_PATH_IMAGE062
；
将得到的
Figure 945498DEST_PATH_IMAGE056
，
Figure 430837DEST_PATH_IMAGE062
进行连接，得到融合后的特征X，计算方法见式（27）；
Figure 100002_DEST_PATH_IMAGE063
(27)
其中，
Figure 866497DEST_PATH_IMAGE056
为融合空间中的傅里叶级数；
Figure 903724DEST_PATH_IMAGE062
为融合空间中的简单几何特征。
进一步地，步骤5的具体内容为：
步骤5.1、将融合后的特征，以及对应的有杆泵工况类型进行组合，记为：{(X，Y) |x m = (特征向量)，y m = (工况类型)}；m表示第m个样本，x m 为第m个样本的特征向量，y m 为第m个样本的工况类型；再将训练集与测试集按照8:2的比例划分；
步骤5.2、针对两个数据集，对X分别进行标准化，标准化方法见式（28）：
Figure 100002_DEST_PATH_IMAGE064
(28)
式中，
Figure 100002_DEST_PATH_IMAGE065
为各数据集中特征值的平均值；
Figure 100002_DEST_PATH_IMAGE066
为各数据集中特征值的标准差；
步骤5.3、以标准化处理后的训练集X_std与有杆泵工况类型标签Y为输入，输入到XGBoost算法中，再次将训练集与测试集按照8:2的比例划分，并进行有杆泵工况诊断模型的训练。
进一步地，步骤6的具体内容为：
优化时，以模型分类准确率为模型评价指标，最高准确率对应参数的取值即所需要的优化后的参数值；模型分类准确率的计算见式（29）：
Figure 100002_DEST_PATH_IMAGE067
(29)
式中，TP为实例是正类且被预测成正类的样本个数；FP为实例是负类且被预测成正类的样本个数；TN为实例是负类且被预测成负类的样本个数；FN为实例是正类且被预测成负类的样本个数；
使用网格搜索法进行有杆泵工况诊断模型参数的优化，主要对其中的lambda、 max_depth，以及learning_rate进行优化，优化范围分别设置为：
Figure 100002_DEST_PATH_IMAGE068
、
Figure 100002_DEST_PATH_IMAGE069
、
Figure 100002_DEST_PATH_IMAGE070
；其中，
Figure 100002_DEST_PATH_IMAGE071
、
Figure 100002_DEST_PATH_IMAGE072
、
Figure 100002_DEST_PATH_IMAGE073
、
Figure 100002_DEST_PATH_IMAGE074
、
Figure 100002_DEST_PATH_IMAGE075
、
Figure 100002_DEST_PATH_IMAGE076
分别为待优化参数取 值区间的上下界。
进一步地，步骤7的具体内容为：
步骤7.1、按照步骤2-步骤6的步骤方法对测试集数据进行特征提取、特征融合、构建出{(X_test，Y_test) | x m_test = (特征向量)，y m_test = (工况类型)}的数据集合、对集合进行标准化，将X_test 作为输入，经过有杆泵工况诊断模型诊断后得到Y_pred；
步骤7.2、针对有杆泵工况诊断模型输出的Y_pred，进行模型诊断结果输出；
以模型分类的准确率Accuracy和召回率Recall作为模型评价指标，模型分类召回率的计算见式（30）：
Figure 100002_DEST_PATH_IMAGE077
(30)
通过绘制出混淆矩阵，计算出该模型对有杆泵各工况诊断的准确率与召回率；Recall值和Accuracy值越接近1，表示分类器特异识别能力和整体分类性能越好。
进一步地，XGBoost算法的目标函数由两部分构成，损失函数和正则项，正则项是用来刻画树复杂度的，损失函数是迭代次数下误差的叠加；在树模型结构中，第e棵树针对样本数据特征向量x m 的预测结果表达，见式（31）：
Figure 100002_DEST_PATH_IMAGE078
(31)
式中，
Figure 100002_DEST_PATH_IMAGE079
表示经过e-1棵树后对样本数据特征向量x m 的预测结果；f e (x m )表示 第e棵树的模型预测结果；k表示从1到e的遍历变量，f k 表示第k棵树的模型预测结果；
XGBoost的目标函数，见式（32）：
Figure 100002_DEST_PATH_IMAGE080
(32)
式中，m表示第m个样本，z为总的样本数量；j表示建立的树模型；T为e迭代次数下 树的数量，l为损失函数，
Figure 100002_DEST_PATH_IMAGE081
为树的复杂度；
损失函数对每一个样本都进行一次损失计算，这里的损失为第e棵树的预测值与真实值之差，复杂度计算是对每棵树的复杂度进行累加；树的复杂度越小，模型的泛化能力越强，其中复杂度函数的表达见式（33）；
Figure 100002_DEST_PATH_IMAGE082
(33)
式中，
Figure 100002_DEST_PATH_IMAGE083
为叶子权重；
Figure 100002_DEST_PATH_IMAGE084
、
Figure 100002_DEST_PATH_IMAGE085
为超参数；
XGBoost的目标函数是关于
Figure 100002_DEST_PATH_IMAGE086
的二次方程，所以损失关于
Figure 710487DEST_PATH_IMAGE086
的导数是线性的，通过 导数等于零即可求得最优解
Figure 100002_DEST_PATH_IMAGE087
；通过训练模型找到一组能使目标函数最小化的系数
Figure 530675DEST_PATH_IMAGE087
， 由此构建出有杆泵工况诊断的XGBoost网络模型。
本发明所带来的有益技术效果：
本发明基于示功图进行有杆泵工况诊断，同时使用了傅里叶系数和简单几何特征两种有杆泵示功图特征提取的方法，保证将示功图的曲线信息全面提取；通过DCA方法进行特征融合，考虑了数据集中类结构间的关系，通过最大化特征集之间的相关性，达到将每组特征中的类分开的目的；基于XGBoost算法进行模型训练，提高了计算精度。本发明既能有效利用油田现有数据，又能提高油田开发现场对有杆泵工况诊断的效率，实现高效诊断有杆泵工况。
附图说明
图1为本发明有杆泵工况诊断方法的流程示意图；
图2为本发明示功图傅里叶系数提取流程示意图；
图3为本发明示功图简单几何特征提取流程示意图；
图4为本发明DCA特征融合的流程示意图；
图5为本发明有杆泵理论示功图；
图6为本发明实施例中有杆泵工况诊断模型训练过程中损失的变化曲线；
图7为本发明实施例中根据混淆矩阵计算出的有杆泵各工况诊断的查准率；
图8为本发明实施例中根据混淆矩阵计算出的有杆泵各工况诊断的召回率。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
本发明提出了一种基于傅里叶变换和几何特征的有杆泵工况诊断方法，主要通过提取示功图的傅里叶系数和简单几何特征作为特征向量。其中，简单几何特征是通过比较油田现场示功图与有杆泵典型故障示功图图像特点总结出来的，能够更有效地区分不同工况下的示功图；傅里叶系数是通过波动方程推导出的符合有杆泵工作特点的特征，能够更加全面的反映示功图曲线的细节特点。
如图1所示，一种基于傅里叶变换和几何特征的有杆泵工况诊断方法，包括以下步骤：
步骤1、对所获有杆泵生产数据进行选择；具体步骤为：
步骤1.1、选择出有杆泵工作过程中各周期中的悬点冲程值、载荷值、理论上载荷值、理论下载荷值；
步骤1.2、通过人工方法清洗异常的示功图；若示功图内包括数据越界、离散点、曲线不闭合、曲线反向或空数据情况即判定为异常的示功图；其中，
数据越界，表现为示功图载荷或冲程数据的数值明显不合逻辑，分为上限超出和下限超出，分别表示为s hg ＜0和s hg ＞L 光杆；其中，L 光杆为光杆长度；s hg 为第h口井第g个冲程值的大小。
离散点，示功图数据发生漂移，破坏示功图曲线的光滑性，曲线中的离散点有三种存在形式，分别为：
①D(t)在t0处无定义；
②D(t)在t0处有定义，但极限
Figure DEST_PATH_IMAGE088
不存在；
③D(t)在t0处有定义，极限
Figure DEST_PATH_IMAGE089
也存在，但
Figure DEST_PATH_IMAGE090
；
式中，D(t)为悬点载荷方程；t0表示离散点存在的位置。
曲线不闭合是示功图数据采集周期出现错误，导致所采数据样本不足一周期，曲 线不闭合可以表示为：
Figure DEST_PATH_IMAGE091
；式中，S h 为第h口井有杆 泵示功图所有点的冲程取值集合；W h 为第h口井有杆泵示功图所有点的载荷取值集合；s hg 为 第h口井有杆泵示功图第g个点的冲程值；w hg 为第h口井有杆泵示功图第g个点的载荷值。
曲线反向是示功图数据传输过程出现错误，导致示功图的面积为负，表示为：
Figure DEST_PATH_IMAGE092
。
空数据是数据未采集或数据传输失败，导致无法形成示功图曲线，表示为：S h =[Null]、W h = [Null]，式中，S h 为第h口井有杆泵示功图所有点的冲程取值集合；W h 为第h口井有杆泵示功图所有点的载荷取值集合；
步骤2、针对示功图进行基于波动方程的傅里叶系数提取；如图2所示，具体步骤为：
步骤2.1、将有杆泵冲程数据与载荷数据提取出来，分别记为：U={S|S h ∈well h }、D={W|W h ∈well h }。
步骤2.2、以第h口井的计算为例，采用吉布斯求解抽油杆运动方程，求解时提出的悬点冲程、载荷的方程见式（1）至（2）：
Figure 696209DEST_PATH_IMAGE001
(1)
Figure 169915DEST_PATH_IMAGE002
(2)
其中，t为时间，n取值为1到
Figure 8515DEST_PATH_IMAGE003
，是累加计算中的遍历变量，
Figure 694711DEST_PATH_IMAGE004
为系数
Figure 839385DEST_PATH_IMAGE005
在级数 为0时的取值，
Figure 483993DEST_PATH_IMAGE006
为系数
Figure 815748DEST_PATH_IMAGE007
在级数为0时的取值，
Figure 571215DEST_PATH_IMAGE008
为所确定的傅里叶级数，
Figure DEST_PATH_IMAGE093
为运动 角速度；
其中，傅里叶系数
Figure DEST_PATH_IMAGE094
、
Figure DEST_PATH_IMAGE095
、
Figure DEST_PATH_IMAGE096
、
Figure DEST_PATH_IMAGE097
的计算方法见式（3）至（6）；
Figure DEST_PATH_IMAGE098
(3)
Figure DEST_PATH_IMAGE099
(4)
Figure DEST_PATH_IMAGE100
(5)
Figure DEST_PATH_IMAGE101
(6)
式中，M为示功图所具数据点数量；i为累加和中的遍历变量，取值为0到M-1。
步骤3、获取示功图曲线数据，进行示功图简单几何特征提取，包括提取示功图对 角线AC、BD斜率，增载线AB边斜率，卸载线CD边斜率，对角线AC、BD的长度，上静载线BC边长 度，下静载线DA边长度，上静载平均载荷与理论上载荷的差值C BC ，下静载平均载荷与理论下 载荷的差值C DA ，以及示功图最大载荷与最小载荷之间的差值
Figure DEST_PATH_IMAGE102
；如图3所示，具体步骤为：
步骤3.1、通过比较有杆泵典型与工况示功图之间的差异，发现每种简单几何特征相对应一定的工况发生情况，故可以提取简单几何特征进行工况的判断；每种简单几何特征与工况判别的对应关系如表1所示，
表1有杆泵示功图与工况判别的比较结果
Figure DEST_PATH_IMAGE103
表中，K AC 为有杆泵示功图对角线AC的斜率；K BD 为有杆泵示功图对角线BD的斜率；K AB 为有杆泵示功图加载线AB线段的斜率；K CD 为有杆泵示功图卸载线CD线段的斜率；L AC 为有 杆泵示功图对角线AC的长度；L BD 为有杆泵示功图对角线BD的长度；L BC 为有杆泵示功图上冲 程BC线段的长度；L DA 为有杆泵示功图下冲程DA线段的长度；C BC 为示功图B、C两点的平均载 荷与理论上载荷的差值；C DA 为示功图D、A两点的平均载荷与理论下载荷的差值；
Figure 993231DEST_PATH_IMAGE030
为示 功图最大载荷与最小载荷的差值；
步骤3.2、示功图A、B、C、D四点冲程载荷值根据如图5所示的有杆泵理论示功图的曲线特点进行计算。当A点冲程最小时，其坐标U a 为最小冲程，D a 为最小冲程对应的载荷；当B点载荷最大时，其坐标D b 为最大载荷，U b 为最大载荷对应的冲程；当C点冲成最大时，其坐标U c 为最大冲程，D c 为最大冲程对应的载荷；当D点载荷最小时，其坐标D d 为最小载荷，U d 为最小载荷对应的冲程。
步骤3.3、以第h口井的计算为例，示功图简单几何特征提取的计算方法见式（7）至（17）：
Figure DEST_PATH_IMAGE104
(7)
式中，K hAC 为第h口井有杆泵示功图对角线AC的斜率；D ha 为第h口井示功图A点的载荷值；U ha 为第h口井示功图A点的冲程值；D hc 为第h口井示功图C点的载荷值；U hc 为第h口井示功图C点的冲程值；
Figure DEST_PATH_IMAGE105
(8)
式中，K hBD 为第h口井有杆泵示功图对角线BD的斜率；D hb 为第h口井示功图B点的载荷值；U hb 为第h口井示功图B点的冲程值；D hd 为第h口井示功图D点的载荷值；U hd 为第h口井示功图D点的冲程值；
Figure DEST_PATH_IMAGE106
(9)
式中，K hAB 为第h口井有杆泵示功图加载线AB线段的斜率；
Figure DEST_PATH_IMAGE107
(10)
式中，K hCD 为第h口井有杆泵示功图卸载线CD线段的斜率；
Figure DEST_PATH_IMAGE108
(11)
式中，L hAC 为第h口井有杆泵示功图对角线AC的长度；
Figure DEST_PATH_IMAGE109
(12)
式中，L hBD 为第h口井有杆泵示功图对角线BD的长度；
Figure DEST_PATH_IMAGE110
(13)
式中，L hBC 为第h口井有杆泵示功图上冲程BC线段的长度；
Figure DEST_PATH_IMAGE111
(14)
式中，L hDA 为第h口井有杆泵示功图下冲程DA线段的长度；
Figure DEST_PATH_IMAGE112
(15)
式中，C hBC 为第h口井示功图B、C两点的平均载荷与理论上载荷的差值，D h理论上为第h口井示功图理论上载荷值；
Figure DEST_PATH_IMAGE113
(16)
式中，C hDA 为第h口井示功图D、A两点的平均载荷与理论下载荷的差值，D h理论下为第h口井示功图理论下载荷值；
Figure DEST_PATH_IMAGE114
(17)
式中，
Figure DEST_PATH_IMAGE115
为第h口井示功图最大载荷与最小载荷的差值；D hmax 为第h口井示功图最 大载荷值；D hmin 为第h口井示功图最小载荷值。
步骤4、进行傅里叶系数与简单几何特征的融合，即采用DiscriminantCorrelation Analysis(DCA)进行特征融合，如图4所示；DAC特征融合的具体步骤为：
步骤4.1、用F1表示步骤2中计算出的所有井的傅里叶系数，F 1={F 11，F 21，...，F h1，...}，其中
Figure DEST_PATH_IMAGE116
，用F2 表示步骤3中计算出的所有井的简单几何特征，F 2={F 12，F 22，...，F h2，...}，其中
Figure DEST_PATH_IMAGE117
。
步骤4.2、假设数据矩阵中的样本是从c个单独的类中收集的。相应地，数据矩阵的N列被分成c个单独的组，其中N v 列属于第v类。分别求出样本F 1，F 2的协方差，计算方法见式（18）与式（19）：
Figure DEST_PATH_IMAGE118
(18)
Figure DEST_PATH_IMAGE119
(19)
式中，
Figure DEST_PATH_IMAGE120
为F 1的维数；
Figure DEST_PATH_IMAGE121
为F 2的维数；p为F 1的维数，q为F 2的维数；
Figure DEST_PATH_IMAGE122
表示 第v类工况的傅里叶系数所构成向量的均值；
Figure DEST_PATH_IMAGE123
表示所有井所有类工况傅里叶系数构成向 量的均值；
Figure DEST_PATH_IMAGE124
表示第v类工况的简单几何特征所构成向量的均值；
Figure DEST_PATH_IMAGE125
表示所有井所有类工 况简单几何特征构成向量的均值；
Figure DEST_PATH_IMAGE126
、
Figure DEST_PATH_IMAGE127
分别通过
Figure DEST_PATH_IMAGE128
、
Figure DEST_PATH_IMAGE129
计算获得；
步骤4.3、利用协方差分别计算出F 1，F 2变换后的矩阵
Figure DEST_PATH_IMAGE130
，
Figure DEST_PATH_IMAGE131
，以
Figure 71390DEST_PATH_IMAGE130
计算为例，计算 方法见式（20）至式（24）：
Figure DEST_PATH_IMAGE132
(20)
式中，P是正交矩阵；
Figure DEST_PATH_IMAGE133
是按降序排列的实非负特征值的对角矩阵。
Figure DEST_PATH_IMAGE134
(21)
式中，Q由P中前r个特征向量组成，对应于矩阵P中最大的r个非零特征值。
Figure DEST_PATH_IMAGE135
(22)
记：
Figure DEST_PATH_IMAGE136
；
Figure DEST_PATH_IMAGE137
(23)
Figure DEST_PATH_IMAGE138
(24)
其中，
Figure 969070DEST_PATH_IMAGE130
为F 1变换后的矩阵；I为单位阵。
同理，也可以得到F 2变换后的矩阵
Figure 528227DEST_PATH_IMAGE131
。
得到新傅里叶系数
Figure 116334DEST_PATH_IMAGE130
和新简单几何特征
Figure 837166DEST_PATH_IMAGE131
后，可以得到
Figure DEST_PATH_IMAGE139
， 对
Figure DEST_PATH_IMAGE140
进行SVD分解，得到式（25），并根据式（26）计算出融合空间内特征
Figure DEST_PATH_IMAGE141
，
Figure DEST_PATH_IMAGE142
；下面以
Figure 284459DEST_PATH_IMAGE141
计算为例。
Figure DEST_PATH_IMAGE143
(25)
其中，H为r×r的左酉矩阵，V为r×r的右酉矩阵；J是一个对角元素非零的对角矩阵；
令
Figure DEST_PATH_IMAGE144
，
Figure DEST_PATH_IMAGE145
，则
Figure DEST_PATH_IMAGE146
(26)
同理可以得到融合空间内的特征
Figure DEST_PATH_IMAGE147
。
将得到的
Figure 681197DEST_PATH_IMAGE141
，
Figure 514023DEST_PATH_IMAGE147
进行连接，得到融合后的特征X，计算方法见式（27）。
Figure DEST_PATH_IMAGE148
(27)
其中，
Figure 546702DEST_PATH_IMAGE141
为融合空间中的傅里叶级数；
Figure 605924DEST_PATH_IMAGE147
为融合空间中的简单几何特征。
步骤5、使用XGBoost算法建立有杆泵工况诊断模型，并进行模型训练；具体步骤为：
步骤5.1、将融合后的特征，以及对应的有杆泵工况类型进行组合，记为：{(X，Y) |x m = (特征向量)，y m = (工况类型)}；m表示第m个样本，x m 为第m个样本的特征向量，y m 为第m个样本的工况类型；再按照一定比例（训练集：测试集=8：2）进行训练集与测试集的划分；
步骤5.2、针对两个数据集，对X分别进行标准化，标准化方法见式（28）：
Figure DEST_PATH_IMAGE149
(28)
式中，
Figure DEST_PATH_IMAGE150
为各数据集中特征值的平均值；
Figure DEST_PATH_IMAGE151
为各数据集中特征值的标准差。
步骤5.3、以标准化处理后的训练集X_std与有杆泵工况类型标签Y为输入，输入到XGBoost算法中，再次按照一定比例（训练集：验证集=8：2）划分训练集与验证集，并进行有杆泵工况诊断模型的训练。
步骤6、进行有杆泵工况诊断模型的参数优化；
优化时，以模型分类准确率为模型评价指标，最高准确率对应参数的取值即所需要的优化后的参数值；模型分类准确率的计算见式（29）：
Figure DEST_PATH_IMAGE152
(29)
式中，TP为实例是正类且被预测成正类的样本个数；FP为实例是负类且被预测成正类的样本个数；TN为实例是负类且被预测成负类的样本个数；FN为实例是正类且被预测成负类的样本个数。
使用网格搜索法进行有杆泵工况诊断模型参数的优化，主要对其中的lambda、 max_depth，以及learning_rate进行优化，优化范围分别设置为：
Figure DEST_PATH_IMAGE153
、
Figure DEST_PATH_IMAGE154
、
Figure DEST_PATH_IMAGE155
；其中，
Figure DEST_PATH_IMAGE156
、
Figure DEST_PATH_IMAGE157
、
Figure DEST_PATH_IMAGE158
、
Figure DEST_PATH_IMAGE159
、
Figure DEST_PATH_IMAGE160
、
Figure DEST_PATH_IMAGE161
分别为待优化参数取值区 间的上下界。
步骤7、对优化后的有杆泵工况诊断模型，进行模型性能评价；具体步骤为：
步骤7.1、按照步骤2-步骤6的步骤方法对测试集数据进行特征提取、特征融合、构建出{(X_test，Y_test) | x m_test = (特征向量)，y m_test = (工况类型)}的数据集合、对集合进行标准化，将X_test 作为输入，经过有杆泵工况诊断模型诊断后得到Y_pred；
步骤7.2、针对有杆泵工况诊断模型输出的Y_pred，进行模型诊断结果输出。
以模型分类的准确率（Accuracy）和召回率(Recall)作为模型评价指标，模型分类召回率的计算见式（30）：
Figure DEST_PATH_IMAGE162
(30)
通过绘制出混淆矩阵，计算出该模型对有杆泵各工况诊断的准确率与召回率。Recall值和Accuracy值越接近1，表示分类器特异识别能力和整体分类性能越好。
步骤8、将训练完成的有杆泵工况诊断模型应用到油田现场，实时采集油田现场数据，进行现场油井有杆泵工作状况的诊断。
另，上述XGBoost算法的目标函数由两部分构成，损失函数和正则项，正则项是用来刻画树复杂度的，损失函数是迭代次数下误差的叠加。在树模型结构中，第e棵树针对样本数据x m 的预测结果可以表达，见式（31）：
Figure DEST_PATH_IMAGE163
(31)
式中，
Figure DEST_PATH_IMAGE164
表示经过e-1棵树后对样本e的预测结果；f e (x m )表示第e棵树的模型 预测结果；k表示从1到e的遍历变量，f k 表示第k棵树的模型预测结果。
XGBoost的目标函数，见式（32）：
Figure DEST_PATH_IMAGE165
(32)
式中，m表示第m个样本，S为总的样本数量；j表示建立的树模型；T为e迭代次数下 树的数量，l为损失函数，
Figure DEST_PATH_IMAGE166
为树的复杂度。
损失函数对每一个样本都进行一次损失计算，这里的损失为第e棵树的预测值与真实值之差，复杂度计算是对每棵树的复杂度进行累加。树的复杂度越小，模型的泛化能力越强，其中复杂度函数的表达见式（33）。
Figure DEST_PATH_IMAGE167
(33)
式中，
Figure DEST_PATH_IMAGE168
为叶子权重；
Figure DEST_PATH_IMAGE169
、
Figure DEST_PATH_IMAGE170
为超参数。
XGBoost的目标函数是关于
Figure DEST_PATH_IMAGE171
的二次方程，所以损失关于
Figure 753222DEST_PATH_IMAGE171
的导数是线性的，通过 导数等于零即可求得最优解（
Figure DEST_PATH_IMAGE172
）。通过训练模型找到一组能使目标函数最小化的系数
Figure 581501DEST_PATH_IMAGE172
， 由此构建出有杆泵工况诊断的XGBoost网络模型。
实施例
下面结合具体油田数据对本发明方法进行描述，同时验证本发明方法的可行性和优越性。本实施例的数据来自某油田的某区块，该区块有杆泵的生产数据共有7542条，按照一定比例（训练集：测试集=8：2）将数据集随机划分为训练集与测试集，其中训练集数据6042条，测试集数据1500条。
该区块所包含的有杆泵的工况类型有：泵工作正常、供液不足、连抽带喷、抽油杆断、气影响、泵漏失、油管漏、活塞脱出工作筒。
本实施例中，使用python编程软件，进行有杆泵工况诊断模型的程序编写。
使用训练集进行有杆泵工况诊断模型的初步建立时，XGBoost的参数取值如表2所示。
表2 XGBoost参数
Figure DEST_PATH_IMAGE173
按照本发明的工况诊断方法对有杆泵现场数据进行建模训练，并预测有杆泵工况情况；具体过程如下：
步骤1、从7542条生产数据中选择出有杆泵工作过程中各周期中悬点冲程值、载荷值、理论上载荷值、理论下载荷值，为有杆泵工况诊断模型的训练做准备；
通过人工方法进行明显异常示功图的清洗，包括数据越界、离散点、曲线不闭合、曲线反向以及空数据；
对于异常示功图剔除后的有杆泵示功图样本曲线数据，对其进行统一，保证每个示功图曲线具有相同数量的数据点；
步骤2、基于所获得的有杆泵工作资料数据，提取计算傅里叶系数特征值；
步骤3、同样，基于所获得的有杆泵工作资料数据，提取简单几何特征；
步骤4、根据计算出的傅里叶系数以及示功图简单几何特征，采用DCA方法进行特征融合；
步骤5、基于XBGoost算法建立有杆泵工况诊断模型，并进行模型训练；
其中，按照公式(28)对X进行标准化时，
Figure DEST_PATH_IMAGE174
和
Figure DEST_PATH_IMAGE175
的具体计算公式为：
Figure DEST_PATH_IMAGE176
，
Figure DEST_PATH_IMAGE177
；
以标准化处理后的6042条训练集的X_std与有杆泵工况类型标签Y为输入，输入到XBGoost算法中，并将该训练集再次按照一定比例（训练集：验证集=8：2）划分训练集与验证集，使用表2所述XBGoost参数值进行有杆泵工况诊断模型的初步建立；
同时，训练时采用对数损失函数进行模型稳定性的验证，模型训练过程的损失曲线如图6所示，横坐标是样本数量，纵坐标是负对数损失，训练时和验证时的损失函数值分别稳定在0.00和0.25附近；
步骤6、设定lambda，max_pepth,learning_rate的取值范围，分别为lambda∈[0,50]、max_depth∈[4,10]、learning_rate∈[0.03,0.3]，并使用网格搜索法进行参数值的优化；
网格搜索法进行参数优选时，每次仅针对一个参数类型进行优选，设置参数变化的步长，进行搜索，最终得到的最优参数为：lambda=50；max_depth=6；learning_rate=0.1；
步骤7、按照前述同样的方法对1500条测试集数据进行特征提取、特征融合、构建出{(X_test，Y_test) | x m_test = (特征向量)，y m_test = (工况类型)}的数据集合、对集合进行标准化，将X_test作为输入，经过优化后的有杆泵工况诊断模型诊断后得到Y_pred；
针对有杆泵工况诊断模型输出的Y_pred，进行模型诊断结果输出；
根据输出结果绘制出混淆矩阵，如表3所示，因此混淆矩阵中对角线上的数字即为测试样本中被准确分类的样本的数量，其对应的横纵坐标为该工况的类型；
表3 输出结果的混淆矩阵
Figure DEST_PATH_IMAGE178
并根据混淆矩阵计算出该模型对有杆泵各工况诊断的结果，如图7、8所示；
图7横坐标是工况名称，纵坐标是准确率；从图7可以看出，使用我们提出的方法建立的有杆泵工况诊断模型对于各种工况预测的准确率都带到了98%及以上；
图8横坐标是工况名称，纵坐标是召回率；从图8可以看出，该模型对于样本数量较多的供液不足、正常工况预测的召回率达到98%及以上，对于样本量少的工况预测的召回率能够达到80%左右。
基于上述性能评价，证明该模型诊断效果极好。
步骤8、输出该有杆泵工况诊断模型，并利用该模型，实时采集油田现场数据，进行现场油井有杆泵工作状况的诊断。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (9)

1.一种基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，包括以下步骤：
步骤1、对所获有杆泵生产数据进行选择；
步骤2、针对示功图进行基于波动方程的傅里叶系数提取；
步骤3、获取示功图曲线数据，进行示功图简单几何特征提取；
步骤4、采用DCA将傅里叶系数与简单几何特征进行融合；
步骤5、使用XGBoost算法建立有杆泵工况诊断模型，并进行模型训练；
步骤6、进行有杆泵工况诊断模型的参数优化；
步骤7、对优化后的有杆泵工况诊断模型，进行模型性能评价；
步骤8、将训练完成的有杆泵工况诊断模型应用到油田现场，实时采集油田现场数据，进行现场油井有杆泵工作状况的诊断。
2.根据权利要求1所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤1的具体内容为：
步骤1.1、选择出有杆泵工作过程中各周期中的悬点冲程值、载荷值、理论上载荷值、理论下载荷值；
步骤1.2、通过人工方法清洗异常的示功图；若示功图内包括数据越界、离散点、曲线不闭合、曲线反向或空数据的情况则被判定为异常的示功图。
3.根据权利要求2所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤2的具体内容为：
步骤2.1、将有杆泵冲程数据与载荷数据提取出来，分别记为：U={S|S h ∈well h }、D={W|W h ∈well h }；式中，S为所有井有杆泵示功图所有点的冲程取值集合，h为第h口井；S h 为第h口井有杆泵示功图所有点的冲程取值集合；W为所有井有杆泵示功图所有点的载荷取值集合；W h 为第h口井有杆泵示功图所有点的载荷取值集合；well h 为第h口井的所有数据信息；
步骤2.2、以第h口井的计算为例，采用吉布斯求解抽油杆运动方程，求解时提出的悬点冲程、载荷的方程见式（1）至（2）：
Figure DEST_PATH_IMAGE001
(1)
Figure DEST_PATH_IMAGE002
(2)
其中，t为时间，n取值为1到
Figure DEST_PATH_IMAGE003
，是累加计算中的遍历变量，
Figure DEST_PATH_IMAGE004
为系数
Figure DEST_PATH_IMAGE005
在级数为0时 的取值，
Figure DEST_PATH_IMAGE006
为系数
Figure DEST_PATH_IMAGE007
在级数为0时的取值，
Figure DEST_PATH_IMAGE008
为所确定的傅里叶级数，
Figure DEST_PATH_IMAGE009
为运动角速度；
其中，傅里叶系数
Figure DEST_PATH_IMAGE010
、
Figure DEST_PATH_IMAGE011
、
Figure DEST_PATH_IMAGE012
、
Figure DEST_PATH_IMAGE013
的计算方法见式（3）至（6）；
Figure DEST_PATH_IMAGE014
(3)
Figure DEST_PATH_IMAGE015
(4)
Figure DEST_PATH_IMAGE016
(5)
Figure DEST_PATH_IMAGE017
(6)
式中，M为示功图所具数据点数量，i为累加和中的遍历变量，取值为0到M-1。
4.根据权利要求3所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在 于，所述步骤3中，简单几何特征的提取包括提取示功图对角线AC、BD斜率，增载线AB边斜 率，卸载线CD边斜率，对角线AC、BD的长度，上静载线BC边长度，下静载线DA边长度，上静载 平均载荷与理论上载荷的差值C BC ，下静载平均载荷与理论下载荷的差值C DA ，以及示功图最 大载荷与最小载荷之间的差值
Figure DEST_PATH_IMAGE018
；提取的具体过程为：
步骤3.1、提取简单几何特征进行工况的判断，每种简单几何特征都对应一定的工况判别情况；
步骤3.2、示功图A、B、C、D四点冲程载荷值根据有杆泵理论示功图的曲线特点进行计算；当A点冲程最小时，其坐标U a 为最小冲程，D a 为最小冲程对应的载荷；当B点载荷最大时，其坐标D b 为最大载荷，U b 为最大载荷对应的冲程；当C点冲成最大时，其坐标U c 为最大冲程，D c 为最大冲程对应的载荷；当D点载荷最小时，其坐标D d 为最小载荷，U d 为最小载荷对应的冲程；
步骤3.3、以第h口井的计算为例，示功图简单几何特征提取的计算方法见式（7）至（17）：
Figure DEST_PATH_IMAGE019
(7)
式中，K hAC 为第h口井有杆泵示功图对角线AC的斜率；D ha 为第h口井示功图A点的载荷值；U ha 为第h口井示功图A点的冲程值；D hc 为第h口井示功图C点的载荷值；U hc 为第h口井示功图C点的冲程值；
Figure DEST_PATH_IMAGE020
(8)
式中，K hBD 为第h口井有杆泵示功图对角线BD的斜率；D hb 为第h口井示功图B点的载荷值；U hb 为第h口井示功图B点的冲程值；D hd 为第h口井示功图D点的载荷值；U hd 为第h口井示功图D点的冲程值；
Figure DEST_PATH_IMAGE021
(9)
式中，K hAB 为第h口井有杆泵示功图加载线AB线段的斜率；
Figure DEST_PATH_IMAGE022
(10)
式中，K hCD 为第h口井有杆泵示功图卸载线CD线段的斜率；
Figure DEST_PATH_IMAGE023
(11)
式中，L hAC 为第h口井有杆泵示功图对角线AC的长度；
Figure DEST_PATH_IMAGE024
(12)
式中，L hBD 为第h口井有杆泵示功图对角线BD的长度；
Figure DEST_PATH_IMAGE025
(13)
式中，L hBC 为第h口井有杆泵示功图上冲程BC线段的长度；
Figure DEST_PATH_IMAGE026
(14)
式中，L hDA 为第h口井有杆泵示功图下冲程DA线段的长度；
Figure DEST_PATH_IMAGE027
(15)
式中，C hBC 为第h口井示功图B、C两点的平均载荷与理论上载荷的差值，D h理论上为第h口井示功图理论上载荷值；
Figure DEST_PATH_IMAGE028
(16)
式中，C hDA 为第h口井示功图D、A两点的平均载荷与理论下载荷的差值，D h理论下为第h口井示功图理论下载荷值；
Figure DEST_PATH_IMAGE029
(17)
式中，
Figure DEST_PATH_IMAGE030
为第h口井示功图最大载荷与最小载荷的差值；D hmax 为第h口井示功图最大载 荷值；D hmin 为第h口井示功图最小载荷值。
5.根据权利要求4所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤4中，DAC特征融合的具体步骤为：
步骤4.1、用F 1表示步骤2中计算出的所有井的傅里叶系数，F 1={F 11，F 21，...，F h1，...}
其中
Figure DEST_PATH_IMAGE031
，用F 2表示步骤3 中计算出的所有井的简单几何特征，F 2={F 12，F 22，...，F h2，...}，其中
Figure DEST_PATH_IMAGE032
；
步骤4.2、假设数据矩阵中的样本是从c个单独的类中收集的；相应地，数据矩阵的N列被分成c个单独的组，其中N v 列属于第v类；分别求出F 1，F 2的协方差，计算方法见式（18）与式（19）：
Figure DEST_PATH_IMAGE033
(18)
Figure DEST_PATH_IMAGE034
(19)
式中，
Figure DEST_PATH_IMAGE035
为F 1的协方差，
Figure DEST_PATH_IMAGE036
为F 2的协方差，p为F 1的维数，q为F 2的维数；
Figure DEST_PATH_IMAGE037
表示 第v类工况的傅里叶系数所构成向量的均值；
Figure DEST_PATH_IMAGE038
表示所有井所有类工况傅里叶系数构成向 量的均值；
Figure DEST_PATH_IMAGE039
表示第v类工况的简单几何特征所构成向量的均值；
Figure DEST_PATH_IMAGE040
表示所有井所有类工 况简单几何特征构成向量的均值；
Figure DEST_PATH_IMAGE041
、
Figure DEST_PATH_IMAGE042
分别通过
Figure DEST_PATH_IMAGE043
、
Figure DEST_PATH_IMAGE044
计算获得；
步骤4.3、利用协方差分别计算出F 1，F 2变换后的矩阵
Figure DEST_PATH_IMAGE045
，
Figure DEST_PATH_IMAGE046
，以
Figure 731518DEST_PATH_IMAGE045
计算为例，计算方 法见式（20）至式（24）：
Figure DEST_PATH_IMAGE047
(20)
式中，P是正交矩阵；
Figure DEST_PATH_IMAGE048
是按降序排列的实非负特征值的对角矩阵；
Figure DEST_PATH_IMAGE049
(21)
式中，Q由P中前r个特征向量组成，对应于矩阵P中最大的r个非零特征值；
Figure DEST_PATH_IMAGE050
(22)
记：
Figure DEST_PATH_IMAGE051
；
Figure DEST_PATH_IMAGE052
(23)
Figure DEST_PATH_IMAGE053
(24)
其中，
Figure 71495DEST_PATH_IMAGE045
为F 1变换后的矩阵；I为单位阵；
同理，得到F 2变换后的矩阵
Figure 382390DEST_PATH_IMAGE046
；
得到新傅里叶系数
Figure 184124DEST_PATH_IMAGE045
和新简单几何特征
Figure 264076DEST_PATH_IMAGE046
后，得到
Figure DEST_PATH_IMAGE054
， 对
Figure DEST_PATH_IMAGE055
进 行SVD分解，得到式（25），并根据式（26）计算出融合空间内特征
Figure DEST_PATH_IMAGE056
，
Figure DEST_PATH_IMAGE057
；下面以
Figure 159307DEST_PATH_IMAGE056
计算为 例；
Figure DEST_PATH_IMAGE058
(25)
其中，H为r×r的左酉矩阵，V为r×r的右酉矩阵；J是一个对角元素非零的对角矩阵；
令
Figure DEST_PATH_IMAGE059
，
Figure DEST_PATH_IMAGE060
，则
Figure DEST_PATH_IMAGE061
(26)
同理，得到融合空间内的特征
Figure DEST_PATH_IMAGE062
；
将得到的
Figure 442652DEST_PATH_IMAGE056
，
Figure 172711DEST_PATH_IMAGE062
进行连接，得到融合后的特征X，计算方法见式（27）；
Figure DEST_PATH_IMAGE063
(27)
其中，
Figure 451377DEST_PATH_IMAGE056
为融合空间中的傅里叶级数；
Figure 975899DEST_PATH_IMAGE062
为融合空间中的简单几何特征。
6.根据权利要求5所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤5的具体内容为：
步骤5.1、将融合后的特征，以及对应的有杆泵工况类型进行组合，记为：{(X，Y) | x m =(特征向量)，y m = (工况类型)}；m表示第m个样本，x m 为第m个样本的特征向量，y m 为第m个样本的工况类型；再将训练集与测试集按照8:2的比例划分；
步骤5.2、针对两个数据集，对X分别进行标准化，标准化方法见式（28）：
Figure DEST_PATH_IMAGE064
(28)
式中，
Figure DEST_PATH_IMAGE065
为各数据集中特征值的平均值；
Figure DEST_PATH_IMAGE066
为各数据集中特征值的标准差；
步骤5.3、以标准化处理后的训练集X_std与有杆泵工况类型标签Y为输入，输入到XGBoost算法中，再次将训练集与测试集按照8:2的比例划分，并进行有杆泵工况诊断模型的训练。
7.根据权利要求6所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤6的具体内容为：
优化时，以模型分类准确率为模型评价指标，最高准确率对应参数的取值即所需要的优化后的参数值；模型分类准确率的计算见式（29）：
Figure DEST_PATH_IMAGE067
(29)
式中，TP为实例是正类且被预测成正类的样本个数；FP为实例是负类且被预测成正类的样本个数；TN为实例是负类且被预测成负类的样本个数；FN为实例是正类且被预测成负类的样本个数；
使用网格搜索法进行有杆泵工况诊断模型参数的优化，主要对其中的lambda、max_ depth，以及learning_rate进行优化，优化范围分别设置为：
Figure DEST_PATH_IMAGE068
、
Figure DEST_PATH_IMAGE069
、
Figure DEST_PATH_IMAGE070
；其中，
Figure DEST_PATH_IMAGE071
、
Figure DEST_PATH_IMAGE072
、
Figure DEST_PATH_IMAGE073
、
Figure DEST_PATH_IMAGE074
、
Figure DEST_PATH_IMAGE075
、
Figure DEST_PATH_IMAGE076
分别为待优化参数 取值区间的上下界。
8.根据权利要求7所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤7的具体内容为：
步骤7.1、按照步骤2-步骤6的步骤方法对测试集数据进行特征提取、特征融合、构建出{(X_test，Y_test) | x m_test = (特征向量)，y m_test = (工况类型)}的数据集合、对集合进行标准化，将X_test 作为输入，经过有杆泵工况诊断模型诊断后得到Y_pred；
步骤7.2、针对有杆泵工况诊断模型输出的Y_pred，进行模型诊断结果输出；
以模型分类的准确率Accuracy和召回率Recall作为模型评价指标，模型分类召回率的计算见式（30）：
Figure DEST_PATH_IMAGE077
(30)
通过绘制出混淆矩阵，计算出该模型对有杆泵各工况诊断的准确率与召回率；Recall值和Accuracy值越接近1，表示分类器特异识别能力和整体分类性能越好。
9.根据权利要求8所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，XGBoost算法的目标函数由两部分构成，损失函数和正则项，正则项是用来刻画树复杂度的，损失函数是迭代次数下误差的叠加；在树模型结构中，第e棵树针对样本数据特征向量x m 的预测结果表达，见式（31）：
Figure DEST_PATH_IMAGE078
(31)
式中，
Figure DEST_PATH_IMAGE079
表示经过e-1棵树后对样本数据特征向量x m 的预测结果；f e (x m )表示第e棵 树的模型预测结果；k表示从1到e的遍历变量，f k 表示第k棵树的模型预测结果；
XGBoost的目标函数，见式（32）：
Figure DEST_PATH_IMAGE080
(32)
式中，m表示第m个样本，z为总的样本数量；j表示建立的树模型；T为e迭代次数下树的 数量，l为损失函数，
Figure DEST_PATH_IMAGE081
为树的复杂度；
损失函数对每一个样本都进行一次损失计算，这里的损失为第e棵树的预测值与真实值之差，复杂度计算是对每棵树的复杂度进行累加；树的复杂度越小，模型的泛化能力越强，其中复杂度函数的表达见式（33）；
Figure DEST_PATH_IMAGE082
(33)
式中，
Figure DEST_PATH_IMAGE083
为叶子权重；
Figure DEST_PATH_IMAGE084
、
Figure DEST_PATH_IMAGE085
为超参数；
XGBoost的目标函数是关于
Figure DEST_PATH_IMAGE086
的二次方程，所以损失关于
Figure 383091DEST_PATH_IMAGE086
的导数是线性的，通过导数 等于零即可求得最优解
Figure DEST_PATH_IMAGE087
；通过训练模型找到一组能使目标函数最小化的系数
Figure 57786DEST_PATH_IMAGE087
，由此 构建出有杆泵工况诊断的XGBoost网络模型。
CN202210407067.1A 2022-04-19 2022-04-19 一种基于傅里叶变换和几何特征的有杆泵工况诊断方法 Active CN114510880B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210407067.1A CN114510880B (zh) 2022-04-19 2022-04-19 一种基于傅里叶变换和几何特征的有杆泵工况诊断方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210407067.1A CN114510880B (zh) 2022-04-19 2022-04-19 一种基于傅里叶变换和几何特征的有杆泵工况诊断方法

Publications (2)

Publication Number Publication Date
CN114510880A true CN114510880A (zh) 2022-05-17
CN114510880B CN114510880B (zh) 2022-07-12

Family

ID=81554925

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210407067.1A Active CN114510880B (zh) 2022-04-19 2022-04-19 一种基于傅里叶变换和几何特征的有杆泵工况诊断方法

Country Status (1)

Country Link
CN (1) CN114510880B (zh)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106837301A (zh) * 2017-03-03 2017-06-13 中国石油化工股份有限公司胜利油田分公司胜利采油厂 一种电功图与示功图整合的抽油机井工况诊断方法
CN106951662A (zh) * 2017-04-12 2017-07-14 东北大学 基于凡尔工作点的有杆泵抽油井井下工况诊断方法
CN108764361A (zh) * 2018-06-01 2018-11-06 北京中油瑞飞信息技术有限责任公司 基于集成学习的游梁式抽油机示功图的工况识别方法
CN111144548A (zh) * 2019-12-23 2020-05-12 北京寄云鼎城科技有限公司 抽油机井工况的识别方法及装置
CN112949196A (zh) * 2021-03-11 2021-06-11 中国石油大学(北京) 一种基于残差神经网络的抽油机井故障诊断方法及系统
WO2022011754A1 (zh) * 2020-07-16 2022-01-20 苏州大学 一种基于自适应流形嵌入动态分布对齐的故障诊断方法
WO2022037012A1 (zh) * 2020-08-19 2022-02-24 江苏大学 一种适用于大规模数据的降维、关联分析方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106837301A (zh) * 2017-03-03 2017-06-13 中国石油化工股份有限公司胜利油田分公司胜利采油厂 一种电功图与示功图整合的抽油机井工况诊断方法
CN106951662A (zh) * 2017-04-12 2017-07-14 东北大学 基于凡尔工作点的有杆泵抽油井井下工况诊断方法
CN108764361A (zh) * 2018-06-01 2018-11-06 北京中油瑞飞信息技术有限责任公司 基于集成学习的游梁式抽油机示功图的工况识别方法
CN111144548A (zh) * 2019-12-23 2020-05-12 北京寄云鼎城科技有限公司 抽油机井工况的识别方法及装置
WO2022011754A1 (zh) * 2020-07-16 2022-01-20 苏州大学 一种基于自适应流形嵌入动态分布对齐的故障诊断方法
WO2022037012A1 (zh) * 2020-08-19 2022-02-24 江苏大学 一种适用于大规模数据的降维、关联分析方法
CN112949196A (zh) * 2021-03-11 2021-06-11 中国石油大学(北京) 一种基于残差神经网络的抽油机井故障诊断方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
DU QIAOLIAN等: "Condition monitoring and fault diagnosis of hydraulic pump based on inherent vibration signals", 《TRANSACTIONS OF THE CHINESE SOCIETY OF AGRICULTURAL ENGINEERING》 *
孙婷婷等: "基于LIBSVM的融合傅里叶幅值与相位的示功图识别方法", 《计算机测量与控制》 *
罗仁泽等: "基于功图分析的有杆抽油泵故障诊断技术研究", 《计算机测量与控制》 *

Also Published As

Publication number Publication date
CN114510880B (zh) 2022-07-12

Similar Documents

Publication Publication Date Title
CN107288617B (zh) 一种提高抽油机井示功图量油精度的方法及系统
CN109460574A (zh) 一种航空发动机剩余寿命的预测方法
CN106321072B (zh) 一种基于泵功图的抽油井故障诊断方法
CN110751076B (zh) 车辆检测方法
CN114429009B (zh) 一种基于元迁移学习的小样本有杆泵井工况诊断方法
CN114241522B (zh) 现场作业安全穿戴识别方法、系统、设备及存储介质
Abdalla et al. Identification of downhole conditions in sucker rod pumped wells using deep neural networks and genetic algorithms (includes associated discussion)
CN106022352A (zh) 基于支持向量机的潜油柱塞泵故障诊断方法
CN114092697A (zh) 注意力融合全局和局部深度特征的建筑立面语义分割方法
CN114879628B (zh) 基于对抗局部最大均值差异的多模态工业过程故障诊断方法
CN111461067A (zh) 基于先验知识映射及修正的零样本遥感影像场景识别方法
Zhou et al. Identification of working condition from sucker-rod pumping wells based on multi-view co-training and hessian regularization of SVM
CN112664185A (zh) 一种基于示功图的抽油机井工况预测方法
CN114510880B (zh) 一种基于傅里叶变换和几何特征的有杆泵工况诊断方法
Yin et al. Imbalanced working states recognition of sucker rod well dynamometer cards based on data generation and diversity augmentation
CN117909881A (zh) 多源数据融合的抽油机的故障诊断方法及装置
CN113495800A (zh) 基于扩展多属性决策的诊断预测数据和特征再认知方法
CN112418223A (zh) 一种基于改进优化的野生动物图像显著性目标检测方法
CN114718861A (zh) 基于深度学习的螺杆泵井工况智能诊断方法
CN117271996A (zh) 一种基于半监督学习的油井堵水效果预测方法
CN114596445B (zh) 一种提升采油机故障诊断精度方法
CN114169594A (zh) 一种基于LSTM-LightGBM变权组合模型的瓦斯浓度预测方法
CN109236277A (zh) 一种基于产生式规则的抽油井故障诊断专家系统
CN118366104B (zh) 一种基于姿态估计的视频监控识别抽油机工况分析方法
Riyadi et al. Comparison of ResNet50V2 and MobileNetV2 Models in Building Architectural Style Classification

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant