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

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

Info

Publication number
CN114510880B
CN114510880B CN202210407067.1A CN202210407067A CN114510880B CN 114510880 B CN114510880 B CN 114510880B CN 202210407067 A CN202210407067 A CN 202210407067A CN 114510880 B CN114510880 B CN 114510880B
Authority
CN
China
Prior art keywords
sucker
rod pump
indicator diagram
load
formula
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.)
Active
Application number
CN202210407067.1A
Other languages
English (en)
Other versions
CN114510880A (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 362301DEST_PATH_IMAGE001
(1)
Figure DEST_PATH_IMAGE002
(2)
其中，t为时间，n取值为1到
Figure 482704DEST_PATH_IMAGE003
，是累加计算中的遍历变量，
Figure DEST_PATH_IMAGE004
为系数
Figure 714971DEST_PATH_IMAGE005
在级数为0时的取值，
Figure DEST_PATH_IMAGE006
为系数
Figure 47863DEST_PATH_IMAGE007
在级数为0时的取值，
Figure DEST_PATH_IMAGE008
为所确定的傅里叶级数，
Figure 550651DEST_PATH_IMAGE009
为运动角速度；
其中，傅里叶系数
Figure DEST_PATH_IMAGE010
、
Figure 576376DEST_PATH_IMAGE011
、
Figure DEST_PATH_IMAGE012
、
Figure 46671DEST_PATH_IMAGE013
的计算方法见式（3）至（6）；
Figure DEST_PATH_IMAGE014
(3)
Figure 432522DEST_PATH_IMAGE015
(4)
Figure DEST_PATH_IMAGE016
(5)
Figure 570243DEST_PATH_IMAGE017
(6)
式中，M为示功图所具数据点数量，i为累加和中的遍历变量，取值为0到M-1。
进一步地，步骤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 766869DEST_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 478123DEST_PATH_IMAGE021
(9)
式中，K hAB 为第h口井有杆泵示功图加载线AB线段的斜率；
Figure DEST_PATH_IMAGE022
(10)
式中，K hCD 为第h口井有杆泵示功图卸载线CD线段的斜率；
Figure 418397DEST_PATH_IMAGE023
(11)
式中，L hAC 为第h口井有杆泵示功图对角线AC的长度；
Figure DEST_PATH_IMAGE024
(12)
式中，L hBD 为第h口井有杆泵示功图对角线BD的长度；
Figure 394312DEST_PATH_IMAGE025
(13)
式中，L hBC 为第h口井有杆泵示功图上冲程BC线段的长度；
Figure DEST_PATH_IMAGE026
(14)
式中，L hDA 为第h口井有杆泵示功图下冲程DA线段的长度；
Figure 761839DEST_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 206727DEST_PATH_IMAGE029
(17)
式中，
Figure 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 435845DEST_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 751420DEST_PATH_IMAGE033
(18)
Figure DEST_PATH_IMAGE034
(19)
式中，
Figure 539117DEST_PATH_IMAGE035
为F 1的协方差，
Figure DEST_PATH_IMAGE036
为F 2的协方差，p为F 1的维数，q为F 2的维数；
Figure 471300DEST_PATH_IMAGE037
表示第v类工况的傅里叶系数所构成向量的均值；
Figure DEST_PATH_IMAGE038
表示所有井所有类工况傅里叶系数构成向量的均值；
Figure 501180DEST_PATH_IMAGE039
表示第v类工况的简单几何特征所构成向量的均值；
Figure DEST_PATH_IMAGE040
表示所有井所有类工况简单几何特征构成向量的均值；
Figure 202420DEST_PATH_IMAGE041
、
Figure DEST_PATH_IMAGE042
分别通过
Figure 177329DEST_PATH_IMAGE043
、
Figure DEST_PATH_IMAGE044
计算获得；
步骤4.3、利用协方差分别计算出F 1，F 2变换后的矩阵
Figure 580497DEST_PATH_IMAGE045
，
Figure DEST_PATH_IMAGE046
，以
Figure 400686DEST_PATH_IMAGE045
计算为例，计算方法见式（20）至式（24）：
Figure 222011DEST_PATH_IMAGE047
(20)
式中，P是正交矩阵；
Figure DEST_PATH_IMAGE048
是按降序排列的实非负特征值的对角矩阵；
Figure 852975DEST_PATH_IMAGE049
(21)
式中，Q由P中前r个特征向量组成，对应于矩阵P中最大的r个非零特征值；
Figure DEST_PATH_IMAGE050
(22)
记：
Figure 494172DEST_PATH_IMAGE051
；
Figure DEST_PATH_IMAGE052
(23)
Figure 367319DEST_PATH_IMAGE053
(24)
其中，
Figure 574309DEST_PATH_IMAGE045
为F 1变换后的矩阵；I为单位阵；
同理，得到F 2变换后的矩阵
Figure 422180DEST_PATH_IMAGE046
；
得到新傅里叶系数
Figure 816252DEST_PATH_IMAGE045
和新简单几何特征
Figure 774981DEST_PATH_IMAGE046
后，得到
Figure DEST_PATH_IMAGE054
， 对
Figure 352041DEST_PATH_IMAGE055
进行SVD分解，得到式（25），并根据式（26）计算出融合空间内特征
Figure DEST_PATH_IMAGE056
，
Figure 839655DEST_PATH_IMAGE057
；下面以
Figure 721023DEST_PATH_IMAGE056
计算为例；
Figure DEST_PATH_IMAGE058
(25)
其中，H为r×r的左酉矩阵，V为r×r的右酉矩阵；J是一个对角元素非零的对角矩阵；
令
Figure 201552DEST_PATH_IMAGE059
，
Figure DEST_PATH_IMAGE060
，则
Figure 586397DEST_PATH_IMAGE061
(26)
同理，得到融合空间内的特征
Figure DEST_PATH_IMAGE062
；
将得到的
Figure 464485DEST_PATH_IMAGE056
，
Figure 833149DEST_PATH_IMAGE062
进行连接，得到融合后的特征X，计算方法见式（27）；
Figure 602522DEST_PATH_IMAGE063
(27)
其中，
Figure 638611DEST_PATH_IMAGE056
为融合空间中的傅里叶级数；
Figure 717295DEST_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 DEST_PATH_IMAGE064
(28)
式中，
Figure 307676DEST_PATH_IMAGE065
为各数据集中特征值的平均值；
Figure DEST_PATH_IMAGE066
为各数据集中特征值的标准差；
步骤5.3、以标准化处理后的训练集X_std与有杆泵工况类型标签Y为输入，输入到XGBoost算法中，再次将训练集与测试集按照8:2的比例划分，并进行有杆泵工况诊断模型的训练。
进一步地，步骤6的具体内容为：
优化时，以模型分类准确率为模型评价指标，最高准确率对应参数的取值即所需要的优化后的参数值；模型分类准确率的计算见式（29）：
Figure 880740DEST_PATH_IMAGE067
(29)
式中，TP为实例是正类且被预测成正类的样本个数；FP为实例是负类且被预测成正类的样本个数；TN为实例是负类且被预测成负类的样本个数；FN为实例是正类且被预测成负类的样本个数；
使用网格搜索法进行有杆泵工况诊断模型参数的优化，主要对其中的lambda、max_depth，以及learning_rate进行优化，优化范围分别设置为：
Figure DEST_PATH_IMAGE068
、
Figure 722401DEST_PATH_IMAGE069
、
Figure DEST_PATH_IMAGE070
；其中，
Figure 722718DEST_PATH_IMAGE071
、
Figure DEST_PATH_IMAGE072
、
Figure 784083DEST_PATH_IMAGE073
、
Figure DEST_PATH_IMAGE074
、
Figure 160838DEST_PATH_IMAGE075
、
Figure 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 374782DEST_PATH_IMAGE077
(30)
通过绘制出混淆矩阵，计算出该模型对有杆泵各工况诊断的准确率与召回率；Recall值和Accuracy值越接近1，表示分类器特异识别能力和整体分类性能越好。
进一步地，XGBoost算法的目标函数由两部分构成，损失函数和正则项，正则项是用来刻画树复杂度的，损失函数是迭代次数下误差的叠加；在树模型结构中，第e棵树针对样本数据特征向量x m 的预测结果表达，见式（31）：
Figure DEST_PATH_IMAGE078
(31)
式中，
Figure 31153DEST_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 64968DEST_PATH_IMAGE081
为树的复杂度；
损失函数对每一个样本都进行一次损失计算，这里的损失为第e棵树的预测值与真实值之差，复杂度计算是对每棵树的复杂度进行累加；树的复杂度越小，模型的泛化能力越强，其中复杂度函数的表达见式（33）；
Figure DEST_PATH_IMAGE082
(33)
式中，
Figure 229102DEST_PATH_IMAGE083
为叶子权重；
Figure DEST_PATH_IMAGE084
、
Figure 297553DEST_PATH_IMAGE085
为超参数；
XGBoost的目标函数是关于
Figure DEST_PATH_IMAGE086
的二次方程，所以损失关于
Figure 639672DEST_PATH_IMAGE086
的导数是线性的，通过导数等于零即可求得最优解
Figure 445604DEST_PATH_IMAGE087
；通过训练模型找到一组能使目标函数最小化的系数
Figure 164161DEST_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 352697DEST_PATH_IMAGE089
也存在，但
Figure DEST_PATH_IMAGE090
；
式中，D(t)为悬点载荷方程；t0表示离散点存在的位置。
曲线不闭合是示功图数据采集周期出现错误，导致所采数据样本不足一周期，曲线不闭合可以表示为：
Figure 114986DEST_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 388972DEST_PATH_IMAGE001
(1)
Figure 176800DEST_PATH_IMAGE002
(2)
其中，t为时间，n取值为1到
Figure 751000DEST_PATH_IMAGE003
，是累加计算中的遍历变量，
Figure 451234DEST_PATH_IMAGE004
为系数
Figure 478096DEST_PATH_IMAGE005
在级数为0时的取值，
Figure 538456DEST_PATH_IMAGE006
为系数
Figure 701584DEST_PATH_IMAGE007
在级数为0时的取值，
Figure 87566DEST_PATH_IMAGE008
为所确定的傅里叶级数，
Figure 850992DEST_PATH_IMAGE093
为运动角速度；
其中，傅里叶系数
Figure 980622DEST_PATH_IMAGE010
、
Figure 263835DEST_PATH_IMAGE011
、
Figure 820719DEST_PATH_IMAGE012
、
Figure 822173DEST_PATH_IMAGE013
的计算方法见式（3）至（6）；
Figure 237717DEST_PATH_IMAGE014
(3)
Figure 109858DEST_PATH_IMAGE015
(4)
Figure 103222DEST_PATH_IMAGE016
(5)
Figure 591972DEST_PATH_IMAGE017
(6)
式中，M为示功图所具数据点数量；i为累加和中的遍历变量，取值为0到M-1。
步骤3、获取示功图曲线数据，进行示功图简单几何特征提取，包括提取示功图对角线AC、BD斜率，增载线AB边斜率，卸载线CD边斜率，对角线AC、BD的长度，上静载线BC边长度，下静载线DA边长度，上静载平均载荷与理论上载荷的差值C BC ，下静载平均载荷与理论下载荷的差值C DA ，以及示功图最大载荷与最小载荷之间的差值
Figure 797826DEST_PATH_IMAGE018
；如图3所示，具体步骤为：
步骤3.1、通过比较有杆泵典型与工况示功图之间的差异，发现每种简单几何特征相对应一定的工况发生情况，故可以提取简单几何特征进行工况的判断；每种简单几何特征与工况判别的对应关系如表1所示，
表1有杆泵示功图与工况判别的比较结果
Figure DEST_PATH_IMAGE094
表中，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 773741DEST_PATH_IMAGE095
为示功图最大载荷与最小载荷的差值；
步骤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_IMAGE096
(7)
式中，K hAC 为第h口井有杆泵示功图对角线AC的斜率；D ha 为第h口井示功图A点的载荷值；U ha 为第h口井示功图A点的冲程值；D hc 为第h口井示功图C点的载荷值；U hc 为第h口井示功图C点的冲程值；
Figure 406847DEST_PATH_IMAGE097
(8)
式中，K hBD 为第h口井有杆泵示功图对角线BD的斜率；D hb 为第h口井示功图B点的载荷值；U hb 为第h口井示功图B点的冲程值；D hd 为第h口井示功图D点的载荷值；U hd 为第h口井示功图D点的冲程值；
Figure DEST_PATH_IMAGE098
(9)
式中，K hAB 为第h口井有杆泵示功图加载线AB线段的斜率；
Figure 586156DEST_PATH_IMAGE099
(10)
式中，K hCD 为第h口井有杆泵示功图卸载线CD线段的斜率；
Figure DEST_PATH_IMAGE100
(11)
式中，L hAC 为第h口井有杆泵示功图对角线AC的长度；
Figure 80853DEST_PATH_IMAGE101
(12)
式中，L hBD 为第h口井有杆泵示功图对角线BD的长度；
Figure DEST_PATH_IMAGE102
(13)
式中，L hBC 为第h口井有杆泵示功图上冲程BC线段的长度；
Figure 662008DEST_PATH_IMAGE103
(14)
式中，L hDA 为第h口井有杆泵示功图下冲程DA线段的长度；
Figure DEST_PATH_IMAGE104
(15)
式中，C hBC 为第h口井示功图B、C两点的平均载荷与理论上载荷的差值，D h理论上为第h口井示功图理论上载荷值；
Figure 715283DEST_PATH_IMAGE105
(16)
式中，C hDA 为第h口井示功图D、A两点的平均载荷与理论下载荷的差值，D h理论下为第h口井示功图理论下载荷值；
Figure DEST_PATH_IMAGE106
(17)
式中，
Figure 913046DEST_PATH_IMAGE095
为第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 460702DEST_PATH_IMAGE107
，用F2表示步骤3中计算出的所有井的简单几何特征，F 2={F 12，F 22，...，F h2，...}，其中
Figure DEST_PATH_IMAGE108
。
步骤4.2、假设数据矩阵中的样本是从c个单独的类中收集的。相应地，数据矩阵的N列被分成c个单独的组，其中N v 列属于第v类。分别求出样本F 1，F 2的协方差，计算方法见式（18）与式（19）：
Figure 427521DEST_PATH_IMAGE109
(18)
Figure DEST_PATH_IMAGE110
(19)
式中，
Figure 878795DEST_PATH_IMAGE111
为F 1的维数；
Figure DEST_PATH_IMAGE112
为F 2的维数；p为F 1的维数，q为F 2的维数；
Figure 563854DEST_PATH_IMAGE113
表示第v类工况的傅里叶系数所构成向量的均值；
Figure DEST_PATH_IMAGE114
表示所有井所有类工况傅里叶系数构成向量的均值；
Figure 649621DEST_PATH_IMAGE115
表示第v类工况的简单几何特征所构成向量的均值；
Figure DEST_PATH_IMAGE116
表示所有井所有类工况简单几何特征构成向量的均值；
Figure 454635DEST_PATH_IMAGE117
、
Figure DEST_PATH_IMAGE118
分别通过
Figure 866025DEST_PATH_IMAGE119
、
Figure DEST_PATH_IMAGE120
计算获得；
步骤4.3、利用协方差分别计算出F 1，F 2变换后的矩阵
Figure 507222DEST_PATH_IMAGE121
，
Figure DEST_PATH_IMAGE122
，以
Figure 147413DEST_PATH_IMAGE121
计算为例，计算方法见式（20）至式（24）：
Figure 354403DEST_PATH_IMAGE123
(20)
式中，P是正交矩阵；
Figure DEST_PATH_IMAGE124
是按降序排列的实非负特征值的对角矩阵。
Figure 405536DEST_PATH_IMAGE125
(21)
式中，Q由P中前r个特征向量组成，对应于矩阵P中最大的r个非零特征值。
Figure DEST_PATH_IMAGE126
(22)
记：
Figure 799608DEST_PATH_IMAGE127
；
Figure DEST_PATH_IMAGE128
(23)
Figure 476446DEST_PATH_IMAGE129
(24)
其中，
Figure 537943DEST_PATH_IMAGE121
为F 1变换后的矩阵；I为单位阵。
同理，也可以得到F 2变换后的矩阵
Figure 556715DEST_PATH_IMAGE122
。
得到新傅里叶系数
Figure 703662DEST_PATH_IMAGE121
和新简单几何特征
Figure 200503DEST_PATH_IMAGE122
后，可以得到
Figure DEST_PATH_IMAGE130
， 对
Figure 850927DEST_PATH_IMAGE131
进行SVD分解，得到式（25），并根据式（26）计算出融合空间内特征
Figure DEST_PATH_IMAGE132
，
Figure 257244DEST_PATH_IMAGE133
；下面以
Figure 891488DEST_PATH_IMAGE132
计算为例。
Figure DEST_PATH_IMAGE134
(25)
其中，H为r×r的左酉矩阵，V为r×r的右酉矩阵；J是一个对角元素非零的对角矩阵；
令
Figure 660861DEST_PATH_IMAGE135
，
Figure DEST_PATH_IMAGE136
，则
Figure 415059DEST_PATH_IMAGE137
(26)
同理可以得到融合空间内的特征
Figure DEST_PATH_IMAGE138
。
将得到的
Figure 667311DEST_PATH_IMAGE132
，
Figure 523271DEST_PATH_IMAGE138
进行连接，得到融合后的特征X，计算方法见式（27）。
Figure 361914DEST_PATH_IMAGE139
(27)
其中，
Figure 986931DEST_PATH_IMAGE132
为融合空间中的傅里叶级数；
Figure 518406DEST_PATH_IMAGE138
为融合空间中的简单几何特征。
步骤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_IMAGE140
(28)
式中，
Figure 579772DEST_PATH_IMAGE141
为各数据集中特征值的平均值；
Figure DEST_PATH_IMAGE142
为各数据集中特征值的标准差。
步骤5.3、以标准化处理后的训练集X_std与有杆泵工况类型标签Y为输入，输入到XGBoost算法中，再次按照一定比例（训练集：验证集=8：2）划分训练集与验证集，并进行有杆泵工况诊断模型的训练。
步骤6、进行有杆泵工况诊断模型的参数优化；
优化时，以模型分类准确率为模型评价指标，最高准确率对应参数的取值即所需要的优化后的参数值；模型分类准确率的计算见式（29）：
Figure 956527DEST_PATH_IMAGE143
(29)
式中，TP为实例是正类且被预测成正类的样本个数；FP为实例是负类且被预测成正类的样本个数；TN为实例是负类且被预测成负类的样本个数；FN为实例是正类且被预测成负类的样本个数。
使用网格搜索法进行有杆泵工况诊断模型参数的优化，主要对其中的lambda、max_depth，以及learning_rate进行优化，优化范围分别设置为：
Figure DEST_PATH_IMAGE144
、
Figure 436050DEST_PATH_IMAGE145
、
Figure DEST_PATH_IMAGE146
；其中，
Figure 360930DEST_PATH_IMAGE147
、
Figure DEST_PATH_IMAGE148
、
Figure 925904DEST_PATH_IMAGE149
、
Figure DEST_PATH_IMAGE150
、
Figure 575191DEST_PATH_IMAGE151
、
Figure DEST_PATH_IMAGE152
分别为待优化参数取值区间的上下界。
步骤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 158488DEST_PATH_IMAGE153
(30)
通过绘制出混淆矩阵，计算出该模型对有杆泵各工况诊断的准确率与召回率。Recall值和Accuracy值越接近1，表示分类器特异识别能力和整体分类性能越好。
步骤8、将训练完成的有杆泵工况诊断模型应用到油田现场，实时采集油田现场数据，进行现场油井有杆泵工作状况的诊断。
另，上述XGBoost算法的目标函数由两部分构成，损失函数和正则项，正则项是用来刻画树复杂度的，损失函数是迭代次数下误差的叠加。在树模型结构中，第e棵树针对样本数据x m 的预测结果可以表达，见式（31）：
Figure DEST_PATH_IMAGE154
(31)
式中，
Figure 500607DEST_PATH_IMAGE155
表示经过e-1棵树后对样本e的预测结果；f e (x m )表示第e棵树的模型预测结果；k表示从1到e的遍历变量，f k 表示第k棵树的模型预测结果。
XGBoost的目标函数，见式（32）：
Figure DEST_PATH_IMAGE156
(32)
式中，m表示第m个样本，S为总的样本数量；j表示建立的树模型；T为e迭代次数下树的数量，l为损失函数，
Figure 552877DEST_PATH_IMAGE157
为树的复杂度。
损失函数对每一个样本都进行一次损失计算，这里的损失为第e棵树的预测值与真实值之差，复杂度计算是对每棵树的复杂度进行累加。树的复杂度越小，模型的泛化能力越强，其中复杂度函数的表达见式（33）。
Figure DEST_PATH_IMAGE158
(33)
式中，
Figure 756588DEST_PATH_IMAGE159
为叶子权重；
Figure DEST_PATH_IMAGE160
、
Figure 945123DEST_PATH_IMAGE161
为超参数。
XGBoost的目标函数是关于
Figure DEST_PATH_IMAGE162
的二次方程，所以损失关于
Figure 458144DEST_PATH_IMAGE162
的导数是线性的，通过导数等于零即可求得最优解（
Figure 512557DEST_PATH_IMAGE163
）。通过训练模型找到一组能使目标函数最小化的系数
Figure 300384DEST_PATH_IMAGE163
，由此构建出有杆泵工况诊断的XGBoost网络模型。
实施例
下面结合具体油田数据对本发明方法进行描述，同时验证本发明方法的可行性和优越性。本实施例的数据来自某油田的某区块，该区块有杆泵的生产数据共有7542条，按照一定比例（训练集：测试集=8：2）将数据集随机划分为训练集与测试集，其中训练集数据6042条，测试集数据1500条。
该区块所包含的有杆泵的工况类型有：泵工作正常、供液不足、连抽带喷、抽油杆断、气影响、泵漏失、油管漏、活塞脱出工作筒。
本实施例中，使用python编程软件，进行有杆泵工况诊断模型的程序编写。
使用训练集进行有杆泵工况诊断模型的初步建立时，XGBoost的参数取值如表2所示。
表2 XGBoost参数
Figure DEST_PATH_IMAGE164
按照本发明的工况诊断方法对有杆泵现场数据进行建模训练，并预测有杆泵工况情况；具体过程如下：
步骤1、从7542条生产数据中选择出有杆泵工作过程中各周期中悬点冲程值、载荷值、理论上载荷值、理论下载荷值，为有杆泵工况诊断模型的训练做准备；
通过人工方法进行明显异常示功图的清洗，包括数据越界、离散点、曲线不闭合、曲线反向以及空数据；
对于异常示功图剔除后的有杆泵示功图样本曲线数据，对其进行统一，保证每个示功图曲线具有相同数量的数据点；
步骤2、基于所获得的有杆泵工作资料数据，提取计算傅里叶系数特征值；
步骤3、同样，基于所获得的有杆泵工作资料数据，提取简单几何特征；
步骤4、根据计算出的傅里叶系数以及示功图简单几何特征，采用DCA方法进行特征融合；
步骤5、基于XBGoost算法建立有杆泵工况诊断模型，并进行模型训练；
其中，按照公式(28)对X进行标准化时，
Figure 343427DEST_PATH_IMAGE165
和
Figure DEST_PATH_IMAGE166
的具体计算公式为：
Figure 292928DEST_PATH_IMAGE167
，
Figure DEST_PATH_IMAGE168
；
以标准化处理后的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 270855DEST_PATH_IMAGE169
并根据混淆矩阵计算出该模型对有杆泵各工况诊断的结果，如图7、8所示；
图7横坐标是工况名称，纵坐标是准确率；从图7可以看出，使用我们提出的方法建立的有杆泵工况诊断模型对于各种工况预测的准确率都带到了98%及以上；
图8横坐标是工况名称，纵坐标是召回率；从图8可以看出，该模型对于样本数量较多的供液不足、正常工况预测的召回率达到98%及以上，对于样本量少的工况预测的召回率能够达到80%左右。
基于上述性能评价，证明该模型诊断效果极好。
步骤8、输出该有杆泵工况诊断模型，并利用该模型，实时采集油田现场数据，进行现场油井有杆泵工作状况的诊断。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (6)

1.一种基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，包括以下步骤：
步骤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、采用吉布斯求解抽油杆运动方程，求解时提出的悬点冲程、载荷的方程见式（1）至（2）：
Figure DEST_PATH_IMAGE001
(1)
Figure 839016DEST_PATH_IMAGE002
(2)
其中，t为时间，n取值为1到
Figure DEST_PATH_IMAGE003
，是累加计算中的遍历变量，
Figure 35118DEST_PATH_IMAGE004
为系数
Figure DEST_PATH_IMAGE005
在级数为0时的取值，
Figure 611593DEST_PATH_IMAGE006
为系数
Figure DEST_PATH_IMAGE007
在级数为0时的取值，
Figure 537960DEST_PATH_IMAGE008
为所确定的傅里叶级数，
Figure DEST_PATH_IMAGE009
为运动角速度；
其中，傅里叶系数
Figure 899803DEST_PATH_IMAGE010
、
Figure DEST_PATH_IMAGE011
、
Figure 253424DEST_PATH_IMAGE012
、
Figure DEST_PATH_IMAGE013
的计算方法见式（3）至（6）；
Figure 317195DEST_PATH_IMAGE014
(3)
Figure DEST_PATH_IMAGE015
(4)
Figure 797986DEST_PATH_IMAGE016
(5)
Figure DEST_PATH_IMAGE017
(6)
式中，M为示功图所具数据点数量，i为累加和中的遍历变量，取值为0到M-1；
步骤3、获取示功图曲线数据，进行示功图简单几何特征提取；
简单几何特征的提取包括提取示功图对角线AC、BD斜率，增载线AB边斜率，卸载线CD边斜率，对角线AC、BD的长度，上静载线BC边长度，下静载线DA边长度，上静载平均载荷与理论上载荷的差值C BC ，下静载平均载荷与理论下载荷的差值C DA ，以及示功图最大载荷与最小载荷之间的差值
Figure 263602DEST_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、示功图简单几何特征提取的计算方法见式（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 788124DEST_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 89924DEST_PATH_IMAGE022
(10)
式中，K hCD 为第h口井有杆泵示功图卸载线CD线段的斜率；
Figure DEST_PATH_IMAGE023
(11)
式中，L hAC 为第h口井有杆泵示功图对角线AC的长度；
Figure 358094DEST_PATH_IMAGE024
(12)
式中，L hBD 为第h口井有杆泵示功图对角线BD的长度；
Figure DEST_PATH_IMAGE025
(13)
式中，L hBC 为第h口井有杆泵示功图上冲程BC线段的长度；
Figure 412638DEST_PATH_IMAGE026
(14)
式中，L hDA 为第h口井有杆泵示功图下冲程DA线段的长度；
Figure DEST_PATH_IMAGE027
(15)
式中，C hBC 为第h口井示功图B、C两点的平均载荷与理论上载荷的差值，D h理论上为第h口井示功图理论上载荷值；
Figure 639220DEST_PATH_IMAGE028
(16)
式中，C hDA 为第h口井示功图D、A两点的平均载荷与理论下载荷的差值，D h理论下为第h口井示功图理论下载荷值；
Figure DEST_PATH_IMAGE029
(17)
式中，
Figure 165666DEST_PATH_IMAGE030
为第h口井示功图最大载荷与最小载荷的差值；D hmax 为第h口井示功图最大载荷值；D hmin 为第h口井示功图最小载荷值；
步骤4、采用DCA将傅里叶系数与简单几何特征进行融合；
步骤5、使用XGBoost算法建立有杆泵工况诊断模型，并进行模型训练；
步骤6、进行有杆泵工况诊断模型的参数优化；
步骤7、对优化后的有杆泵工况诊断模型，进行模型性能评价；
步骤8、将训练完成的有杆泵工况诊断模型应用到油田现场，实时采集油田现场数据，进行现场油井有杆泵工作状况的诊断。
2.根据权利要求1所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤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 237527DEST_PATH_IMAGE032
；
步骤4.2、假设数据矩阵中的样本是从c个单独的类中收集的；相应地，数据矩阵的N列被分成c个单独的组，其中N v 列属于第v类；分别求出F 1，F 2的协方差，计算方法见式（18）与式（19）：
Figure DEST_PATH_IMAGE033
(18)
Figure 412156DEST_PATH_IMAGE034
(19)
式中，
Figure DEST_PATH_IMAGE035
为F 1的协方差，
Figure 294793DEST_PATH_IMAGE036
为F 2的协方差，p为F 1的维数，q为F 2的维数；
Figure DEST_PATH_IMAGE037
表示第v类工况的傅里叶系数所构成向量的均值；
Figure 554873DEST_PATH_IMAGE038
表示所有井所有类工况傅里叶系数构成向量的均值；
Figure DEST_PATH_IMAGE039
表示第v类工况的简单几何特征所构成向量的均值；
Figure 430425DEST_PATH_IMAGE040
表示所有井所有类工况简单几何特征构成向量的均值；
Figure DEST_PATH_IMAGE041
、
Figure 210293DEST_PATH_IMAGE042
分别通过
Figure DEST_PATH_IMAGE043
、
Figure 778678DEST_PATH_IMAGE044
计算获得；
步骤4.3、利用协方差分别计算出F 1，F 2变换后的矩阵
Figure DEST_PATH_IMAGE045
，
Figure 260475DEST_PATH_IMAGE046
，
Figure 674139DEST_PATH_IMAGE045
计算方法见式（20）至式（24）：
Figure DEST_PATH_IMAGE047
(20)
式中，P是正交矩阵；
Figure 308513DEST_PATH_IMAGE048
是按降序排列的实非负特征值的对角矩阵；
Figure DEST_PATH_IMAGE049
(21)
式中，Q由P中前r个特征向量组成，对应于矩阵P中最大的r个非零特征值；
Figure 47799DEST_PATH_IMAGE050
(22)
记：
Figure DEST_PATH_IMAGE051
；
Figure 16892DEST_PATH_IMAGE052
(23)
Figure DEST_PATH_IMAGE053
(24)
其中，
Figure 982050DEST_PATH_IMAGE045
为F 1变换后的矩阵；I为单位阵；
同理，得到F 2变换后的矩阵
Figure 720199DEST_PATH_IMAGE046
；
得到新傅里叶系数
Figure 99227DEST_PATH_IMAGE045
和新简单几何特征
Figure 555616DEST_PATH_IMAGE046
后，得到
Figure 45504DEST_PATH_IMAGE054
， 对
Figure DEST_PATH_IMAGE055
进行SVD分解，得到式（25），并根据式（26）计算出融合空间内特征
Figure 388891DEST_PATH_IMAGE056
，
Figure DEST_PATH_IMAGE057
；
Figure 469980DEST_PATH_IMAGE056
计算过程为：
Figure 148086DEST_PATH_IMAGE058
(25)
其中，H为r×r的左酉矩阵，V为r×r的右酉矩阵；J是一个对角元素非零的对角矩阵；
令
Figure DEST_PATH_IMAGE059
，
Figure 707243DEST_PATH_IMAGE060
，则
Figure DEST_PATH_IMAGE061
(26)
同理，得到融合空间内的特征
Figure 436296DEST_PATH_IMAGE062
；
将得到的
Figure 891548DEST_PATH_IMAGE056
，
Figure 56950DEST_PATH_IMAGE062
进行连接，得到融合后的特征X，计算方法见式（27）；
Figure DEST_PATH_IMAGE063
(27)
其中，
Figure 419798DEST_PATH_IMAGE056
为融合空间中的傅里叶级数；
Figure 987046DEST_PATH_IMAGE062
为融合空间中的简单几何特征。
3.根据权利要求2所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤5的具体内容为：
步骤5.1、将融合后的特征，以及对应的有杆泵工况类型进行组合，记为：{(X，Y) | x m =(特征向量)，y m = (工况类型)}；m表示第m个样本，x m 为第m个样本的特征向量，y m 为第m个样本的工况类型；再将训练集与测试集按照8:2的比例划分；
步骤5.2、针对两个数据集，对X分别进行标准化，标准化方法见式（28）：
Figure 629511DEST_PATH_IMAGE064
(28)
式中，
Figure DEST_PATH_IMAGE065
为各数据集中特征值的平均值；
Figure 547788DEST_PATH_IMAGE066
为各数据集中特征值的标准差；
步骤5.3、以标准化处理后的训练集X_std与有杆泵工况类型标签Y为输入，输入到XGBoost算法中，再次将训练集与测试集按照8:2的比例划分，并进行有杆泵工况诊断模型的训练。
4.根据权利要求3所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤6的具体内容为：
优化时，以模型分类准确率为模型评价指标，最高准确率对应参数的取值即所需要的优化后的参数值；模型分类准确率的计算见式（29）：
Figure DEST_PATH_IMAGE067
(29)
式中，TP为实例是正类且被预测成正类的样本个数；FP为实例是负类且被预测成正类的样本个数；TN为实例是负类且被预测成负类的样本个数；FN为实例是正类且被预测成负类的样本个数；
使用网格搜索法进行有杆泵工况诊断模型参数的优化，主要对其中的lambda、max_depth，以及learning_rate进行优化，优化范围分别设置为：
Figure 448748DEST_PATH_IMAGE068
、
Figure DEST_PATH_IMAGE069
、
Figure 878025DEST_PATH_IMAGE070
；其中，
Figure DEST_PATH_IMAGE071
、
Figure 206238DEST_PATH_IMAGE072
、
Figure DEST_PATH_IMAGE073
、
Figure 611812DEST_PATH_IMAGE074
、
Figure DEST_PATH_IMAGE075
、
Figure 316463DEST_PATH_IMAGE076
分别为待优化参数取值区间的上下界。
5.根据权利要求4所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，所述步骤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，表示分类器特异识别能力和整体分类性能越好。
6.根据权利要求5所述基于傅里叶变换和几何特征的有杆泵工况诊断方法，其特征在于，XGBoost算法的目标函数由两部分构成，损失函数和正则项，正则项是用来刻画树复杂度的，损失函数是迭代次数下误差的叠加；在树模型结构中，第e棵树针对样本数据特征向量x m 的预测结果表达，见式（31）：
Figure 874614DEST_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 373728DEST_PATH_IMAGE080
(32)
式中，m表示第m个样本，z为总的样本数量；j表示建立的树模型；T为e迭代次数下树的数量，l为损失函数，
Figure DEST_PATH_IMAGE081
为树的复杂度；
损失函数对每一个样本都进行一次损失计算，这里的损失为第e棵树的预测值与真实值之差，复杂度计算是对每棵树的复杂度进行累加；树的复杂度越小，模型的泛化能力越强，其中复杂度函数的表达见式（33）；
Figure 266598DEST_PATH_IMAGE082
(33)
式中，
Figure DEST_PATH_IMAGE083
为叶子权重；
Figure 260093DEST_PATH_IMAGE084
、
Figure DEST_PATH_IMAGE085
为超参数；
XGBoost的目标函数是关于
Figure 922018DEST_PATH_IMAGE086
的二次方程，所以损失关于
Figure 326455DEST_PATH_IMAGE086
的导数是线性的，通过导数等于零即可求得最优解
Figure DEST_PATH_IMAGE087
；通过训练模型找到一组能使目标函数最小化的系数
Figure 457353DEST_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 CN114510880A (zh) 2022-05-17
CN114510880B true 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 (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
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 江苏大学 一种适用于大规模数据的降维、关联分析方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106837301A (zh) * 2017-03-03 2017-06-13 中国石油化工股份有限公司胜利油田分公司胜利采油厂 一种电功图与示功图整合的抽油机井工况诊断方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
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
Condition monitoring and fault diagnosis of hydraulic pump based on inherent vibration signals;Du QiaoLian等;《Transactions of the Chinese Society of Agricultural Engineering》;20070430;第23卷(第4期);全文 *
基于LIBSVM的融合傅里叶幅值与相位的示功图识别方法;孙婷婷等;《计算机测量与控制》;20181023(第10期);全文 *
基于功图分析的有杆抽油泵故障诊断技术研究;罗仁泽等;《计算机测量与控制》;20160125(第01期);全文 *

Also Published As

Publication number Publication date
CN114510880A (zh) 2022-05-17

Similar Documents

Publication Publication Date Title
CN107288617B (zh) 一种提高抽油机井示功图量油精度的方法及系统
CN108985380B (zh) 一种基于聚类集成的转辙机故障识别方法
CN106321072B (zh) 一种基于泵功图的抽油井故障诊断方法
CN114429009B (zh) 一种基于元迁移学习的小样本有杆泵井工况诊断方法
CN109033178B (zh) 一种挖掘能见度多维时空数据之间格兰杰因果关系的方法
CN110751076B (zh) 车辆检测方法
Abdalla et al. Identification of downhole conditions in sucker rod pumped wells using deep neural networks and genetic algorithms (includes associated discussion)
CN109630095A (zh) 一种基于多视角学习的抽油机井工况识别方法及系统
CN106022352A (zh) 基于支持向量机的潜油柱塞泵故障诊断方法
CN106951662A (zh) 基于凡尔工作点的有杆泵抽油井井下工况诊断方法
CN114239404A (zh) 一种基于多目标优化的材料智能优化设计方法
CN110276493A (zh) 一种油井检泵周期预测方法、装置及存储介质
Zhou et al. Identification of working condition from sucker-rod pumping wells based on multi-view co-training and hessian regularization of SVM
CN110288257A (zh) 一种深度超限示功图学习方法
CN112664185A (zh) 一种基于示功图的抽油机井工况预测方法
CN114510880B (zh) 一种基于傅里叶变换和几何特征的有杆泵工况诊断方法
CN114718861A (zh) 基于深度学习的螺杆泵井工况智能诊断方法
CN116011351B (zh) 一种基于聚类算法和WideDeep网络的油井合理沉没度确定方法
CN117172360A (zh) 一种基于mlp和高效pso的钻井机械钻速优化方法、系统、设备及介质
CN117078956A (zh) 一种基于点云多尺度并行特征提取和注意力机制的点云分类分割网络
Sharaf Beam pump dynamometer card prediction using artificial neural networks
CN117171713A (zh) 一种基于轴承寿命的交叉自适应深度迁移学习方法及系统
Xu et al. A progressive fault diagnosis method for rolling bearings based on VMD energy entropy and a deep adversarial transfer network
CN113137211B (zh) 一种基于模糊综合决策的油井生产参数自适应调控方法
CN113495800A (zh) 基于扩展多属性决策的诊断预测数据和特征再认知方法

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