由于抽油机井具有油田开采普及率高、工况复杂多变、故障发生率高等特点, 及时精准地识别油井工况对油井极大实现降本增效具有重大意义。目前抽油机井工况识别方法主要有以下几类:①基于示功图识别方法; ②基于电参数识别方法; ③基于多源数据识别方法。多数工况识别方法基于示功图识别技术[1-6], 主要是利用泵功图(实测地面示功图或实测电参数通过相应模型计算得到)或实测地面示功图结合人工智能方法(矩特征、小波包变换、神经网络、支持向量机等)进行工况识别。基于电参数识别方法主要是利用电功图[7](由实测电参数计算得到)或实测电参数[8]进行工况识别。目前对基于多源数据识别方法的研究很少, 主要是利用泵功图结合油井生产信息(产量、抽汲参数、井况数据等)进行工况识别[9-10]。上述研究取得了较好成效, 但仍存在4个方面的问题:一是受阻尼系数、“除零”问题影响, 通过模型计算得到的泵功图和电功图存在精度误差, 影响了对特征参数值的精确计算; 二是对于机电液耦合复杂系统, 同一现象有时可能由不同原因造成, 用单一信息源判断油井工况易产生误报现象; 三是受先期油井海量实时数据采集和存储技术限制、井况复杂多变、生产统计数据不可靠等因素影响, 原有的基于多源信息工况识别的模型鲁棒性较差。另外, 传统的神经网络、支持向量机等方法主要采用特征连接方式进行识别, 对海量多源数据的智能信息处理已难以再取得理想效果[11]; 四是多数工况识别方法需要大量标记工况训练样本, 而实际工程中标记工况样本获取很难。无标记样本训练的工况识别方法又极大浪费专家资源。此外, 多数工况特征提取方法计算量大且对特征数据之间的非相关性要求较高。这些都影响了已有工况识别方法的实际应用。大数据和油气生产物联网环境下, 抽油机井采油生产系统采集和存储了海量多源实时信息, 如实测的地面示功图、电参数、井口温度等, 也获取了海量未知工况样本。多视角学习方法能够有效地解决海量而又复杂目标问题, 突破传统单一特征学习精确度低的局限性, 充分利用目标本身的多特征性质, 各特征之间相互补充完善, 提高学习算法性能。研究表明, 通过学习多视角特征和流形信息可以进一步提高识别结果[12-13]。与传统信息处理方法不同, 半监督学习(semi-supervised learning, SSL)方法大多基于流形假设, 其中Hessian局部线性特征映射算法具有一个明显的优势:当标记样本较少时, 其正则化可以极大改善半监督分类结果[14]。基于此, 针对大数据下抽油机井采油生产系统特点及其工况识别研究目前所存在的问题, 笔者有效利用大数据下抽油机井采油生产系统的海量多源实时信息, 突破单一信息源识别的局限性和传统特征连接多源识别方法的技术瓶颈, 提出一种基于Hessian正则化多视角学习的抽油机井工况识别新方法, 并建立抽油机井工况识别模型, 以进一步提高抽油机井工况识别精准率。
1 Hessian正则化多视角学习方法 1.1 Hessian正则化思想在半监督学习SSL中,Hessian正则化[15-17]反映了样本数据分布流形的高阶信息,能够更准确地描述数据的内在局部几何特征。Hessian正则化算法既可以很好地匹配训练样本区域内数据,也能较好地预测区域外数据。利用Hessian正则化可以在标记样本较少下极大改善半监督分类结果。
在SSL中,假设有l个标记样本L={(xi1, xi2, …, xiVyi)}i=1l, u个未标记样本U={(xi1, xi2, …, xiV)}i=l+1l+u,V代表视角数,xi代表第i个样本,xk,k∈{1, 2, …, V}代表第k个视角的特征,xik代表第i个样本的第k个视角的特征矢量,yi∈{-1, 1}代表样本xi的标记。标记样本概率记为p, 无标记样本边缘概率记为px, px构成一个紧流形M,在M上条件分布概率p(y|x)光滑变化。由于相邻两个样本xi和xj之间要有相似的条件分布p(yi|xi)和p(yj|xj),因此对M的局部几何特性的准确研究非常重要。Hessian不仅有丰富的零空间,还可以驱动分类函数按流形线性变化,因此Hessian正则化能够更好地表征数据分布的内在几何特性。
为便于描述,将本文中重要的变量符号列表说明,如表 1所示。
对于半监督分类问题,其目标函数可表示为
$ \mathop {\min }\limits_{f \in {C^\infty }\left( M \right)} \frac{1}{l}\sum\limits_{i = 1}^l {{L_{{\rm{loss}}}}\left( {{x_i},{y_i},f} \right)} + G\left( {\left\| \mathit{\boldsymbol{f}} \right\|} \right). $ | (1) |
式中,C∞(M)为流形M上的平滑函数集;G(‖f‖)为函数惩罚项与流形惩罚项之和;Lloss(xi, yi, f)为损失函数,是区分各种分类方法的关键因素,包含误差项和正则化项两部分,表示如下:
$ J\left( \omega \right) = \sum\limits_i {{L_{{\rm{loss}}}}\left( {{m_i}\left( \omega \right)} \right)} + \gamma R\left( \omega \right). $ |
式中,Lloss(mi(ω))为误差项;R(ω)为正则化项。
添加Hessian流形正则化项后,式(1)可优化如下:
$ \mathop {\min }\limits_{f \in {H_k}} \frac{1}{l}\sum\limits_{i = 1}^l {\psi \left( {f,{x_i},{y_i}} \right)} + {\gamma _k}\left\| \mathit{\boldsymbol{f}} \right\|_k^2 + {\gamma _I}\left\| \mathit{\boldsymbol{f}} \right\|_I^2. $ | (2) |
式中,ψ为损失函数;在一个恰当的再生核希尔伯特空间Hk中,‖f‖k2为分类复杂度惩罚项,‖f‖I2为流形惩罚项,γk为‖f‖k2的调节参数,γI为‖f‖I2的调节参数。
添加多视角核K和多视角Hessian矩阵H后,目标函数式(2)可实现多视角特征表示。
多视角核K:假设在第k视角Kk是有效的(对称、正定)内积核,k=1, 2, …, V,则
$ \mathit{\boldsymbol{K}} = \sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} , $ | (3) |
$ {\rm{s}}.\;{\rm{t}}.\;\sum\nolimits_{k = 1}^V {{\theta ^k}} = 1,{\theta ^k} \ge 0,k = 1,2, \cdots ,V. $ |
多视角Hessian矩阵H:假设Hj是第j视角的Hessian矩阵,且Hj是半正定的,j=1, 2, …, V,则
$ \mathit{\boldsymbol{K}} = \sum\nolimits_{j = 1}^V {{\beta ^j}{\mathit{\boldsymbol{H}}^j}} , $ | (4) |
$ {\rm{s}}.\;{\rm{t}}.\;\sum\nolimits_{j = 1}^V {{\beta ^j}} = 1,{\beta ^j} \ge 0,j = 1,2, \cdots ,V. $ |
添加式(3)和式(4)到式(2)中,式(2)可表示如下:
$ \begin{array}{l} \mathop {\min }\limits_{f \in {H_k},\theta \in {R^v},\beta \in {R^v}} \frac{1}{l}\sum\limits_{i = 1}^l {\psi \left( {f,{x_i},{y_i}} \right)} + {\gamma _k}\sum\limits_{k = 1}^V {{\theta ^k}\left\| \mathit{\boldsymbol{f}} \right\|_{K\left( k \right)}^2} + \\ {\gamma _I}\sum\limits_{j = 1}^V {{\beta ^j}\left\| \mathit{\boldsymbol{f}} \right\|_{I\left( j \right)}^2} + {\gamma _\theta }\left\| \mathit{\boldsymbol{\theta }} \right\|_2^2 + {\gamma _\beta }\left\| \mathit{\boldsymbol{\beta }} \right\|_2^2. \end{array} $ | (5) |
$ \begin{array}{l} {\rm{s}}.\;{\rm{t}}.\;\sum\nolimits_{k = 1}^V {{\theta ^k}} = 1,{\theta ^k} \ge 0,k = 1,2, \cdots ,V,\\ \sum\nolimits_{j = 1}^V {{\beta ^j}} = 1,{\beta ^j} \ge 0,j = 1,2, \cdots ,V. \end{array} $ |
式中,‖θ‖22和‖β‖22为正则项;γθ, γβ∈R+。
根据表示定理[18],给定一个凸的损失函数,式(5)的最优化解可表示为
$ {f^ * } = \sum\nolimits_{i = 1}^{l + u} {{\alpha _i}\mathit{\boldsymbol{K}}\left( {{x_i},x} \right)} . $ | (6) |
对于多视角核K正则化项和多视角Hessian矩阵H正则化项有:
$ \left\{ \begin{array}{l} \left\| \mathit{\boldsymbol{f}} \right\|_K^2 = {\mathit{\boldsymbol{f}}^{\rm{T}}}\mathit{\boldsymbol{Kf}} = \sum\nolimits_{k = 1}^V {{\theta ^k}\left\| \mathit{\boldsymbol{f}} \right\|_{K\left( k \right)}^2} ,\\ \left\| \mathit{\boldsymbol{f}} \right\|_I^2 = {\mathit{\boldsymbol{f}}^{\rm{T}}}\mathit{\boldsymbol{Hf = }}\sum\nolimits_{j = 1}^V {{\beta ^j}\left\| \mathit{\boldsymbol{f}} \right\|_{I\left( j \right)}^2} . \end{array} \right. $ | (7) |
将式(6)和式(7)代入式(5)中,则目标函数的最优化可表示如下:
$ \begin{array}{l} \mathop {\min }\limits_{f \in {H_k},\theta \in {R^v},\beta \in {R^v}} \frac{1}{l}\sum\limits_{i = 1}^l {\psi \left( {f,{x_i},{y_i}} \right)} + {\gamma _k}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\mathit{\boldsymbol{K\alpha }} + \\ {\gamma _I}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\mathit{\boldsymbol{KHK\alpha }} + {\gamma _\theta }\left\| \mathit{\boldsymbol{\theta }} \right\|_2^2 + {\gamma _\beta }\left\| \mathit{\boldsymbol{\beta }} \right\|_2^2. \end{array} $ | (8) |
$ \begin{array}{l} {\rm{s}}.\;{\rm{t}}.\;\sum\nolimits_{k = 1}^V {{\theta ^k}} = 1,{\theta ^k} \ge 0,k = 1,2, \cdots ,V,\\ \sum\nolimits_{j = 1}^V {{\beta ^j}} = 1,{\beta ^j} \ge 0,j = 1,2, \cdots ,V. \end{array} $ |
其中
$ \mathit{\boldsymbol{K}} = \sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} ,\mathit{\boldsymbol{H}} = \sum\nolimits_{j = 1}^V {{\beta ^j}{\mathit{\boldsymbol{K}}^j}} . $ |
选择合适的损失函数并通过式(6)优化后,对目标函数式(8)的求解就变为对系数求解,采用交替最优化方法[19]即可求解。
2 Hessian正则化多视角logistic回归抽油机井工况识别模型 2.1 视角选取目前抽油机井采油生产管理系统采集和存储了海量的多源实时信息,如实测的地面示功图、电参数、井口温度、井口压力等。这些多源实时信息可以从地面、井筒、地层全方位地及时反映抽油机井工况。实测地面示功图主要反映抽油机井泵及地下情况[20],实测电功率信号可以反映抽油机井地面及地下情况[7-8],二者基本上可以全面反映抽油机井地面、井筒及地层状况,但还有少量抽油机井工况不但示功图形状相似且电功图特点类似[7, 20],如连抽带喷和油管底部严重漏失,可以借助井口温度等实时信息进行准确识别。连抽带喷工况下,井温上升,井压下降,所有漏失工况下,井温和井压都下降,通过井温可以准确识别出这些工况。
综上所述,选择实测地面示功图信号、实测电功率信号和实测井口温度3个信息源作为特征视角。
2.2 视角特征提取对目前抽油机井工况识别研究中实测地面示功图和实测电参数进行特征提取时所存在的问题进行分析,为更好地提高工况识别精准度和适合工程实际应用,充分利用先验知识和专家经验,结合机制分析进行特征提取。
理论示功图如图 1所示,其中横坐标代表位移,用S表示,纵坐标代表载荷,用P表示。Sp为活塞冲程(有效冲程),Pl为活塞上液柱重量,A(E)点为游动凡尔关闭点,也是下死点,B点为固定凡尔打开点,C(F)点为固定凡尔关闭点,也是上死点,D点为游动凡尔打开点。
结合理论示功图对实测地面示功图、实测电功率、实测井口温度的特征提取分析与计算如下:
实测地面示功图依据抽油泵一个工作周期内冲程、冲次、功图面积、活塞上液柱重量、最大和最小载荷、有效冲程及冲程损失的变化进行特征提取,其特征参数值精准计算可以通过对实测示功图上泵凡尔开闭点位置的精确提取实现。
实测地面示功图的特征数据有冲程、冲次、示功图实际面积、最大和最小载荷、最大最小载荷比、活塞上液注重量、有效冲程、加载和卸载冲程损失、提前加载和卸载位置12个特征参数,具体计算如下:
冲程、冲次、最大和最小载荷可以从实时采集的示功图数据中直接获得;功图实际面积等于示功图采集点所围成的封闭曲线面积;最大最小载荷比等于最大载荷与最小载荷的比值;活塞上液柱重量等于最大载荷与最小载荷之差;有效冲程等于游动凡尔打开点与游动凡尔关闭点之间的位移差;加载冲程损失等于固定凡尔打开点与游动凡尔关闭点之间的位移差;卸载冲程损失等于固定凡尔关闭点与游动凡尔打开点之间的位移差;提前加载位置等于从游动凡尔关闭点开始到游动凡尔打开点之间斜率正负反向的第一个点的位移;提前卸载位置等于从固定凡尔关闭点开始到固定凡尔打开点之间斜率正负反向的第一个点的位移。
实测电功率信号依据“功特征”和曲线下面积“AUC(area under curve)特征”[7-8]进行特征提取,其特征参数值精准计算可以通过在实测示功图上精确提取上、下死点位置实现。实测电功率信号的特征数据有上行功、下行功、周期功、上行面积、下行面积、周期面积、平衡率(上、下行程做功比)7个特征参数,具体计算如下:
从实时采集的示功图数据中可以直接得到上死点和下死点,上死点为位移最大点,下死点为位移最小点,即功图数据的起始点。上行功等于上冲程(从下死点到上死点)时间段内所作的功;下行功等于下冲程(从上死点到下死点)时间段内所作的功;周期功等于上行功与下行功之和;上行面积等于上冲程时间段内电功率信号曲线与时间水平轴所围曲线面积;下行面积等于下冲程时间段内电功率信号曲线与时间水平轴所围曲线面积;周期面积等于上行面积与下行面积之和;平衡率等于上行功与下行功的比值。
目前对实测井口温度的特征提取没有任何参考文献。受地面、地层等环境变化因素及数据采集精度影响,现场采集到的井口温度数据大多并不严格遵循对应工况特征状态规律,但从本质上可以反映出对应工况每个冲程的热能耗损。实测井口温度信号特征参数值精准计算可以通过对实测示功图精确提取上、下死点位置实现。实测井口温度信号的特征数据有上行热能(温度)耗损、下行热能(温度)耗损、周期热能(温度)耗损3个特征参数,具体计算如下:
一个冲程内井口温度实时采集点个数通常少于示功图实时采集点个数,可以采用插值拟合方法将二者实时采集点同步,同时从实时采集的示功图数据中得到上死点和下死点。上行热能(温度)耗损等于上冲程时间段内所耗损的热能(温度);下行热能(温度)耗损等于下冲程时间段内所耗损的热能(温度);周期热能(温度)耗损等于上行热能(温度)耗损与下行热能(温度)耗损之和。
2.3 工况识别模型建立基于Hessian正则化多视角logistic回归抽油机井工况识别模型建立流程如图 2所示,主要包括损失函数选择和算法求解两大部分。
损失函数是区别各个分类方法的核心因素,损失函数越小,模型鲁棒性越好,而且损失函数尽量是一个凸函数,便于收敛计算。常用的损失函数有Gold Standard loss(理想状态)、Hinge loss(用于SVM)、Log loss(用于logistic回归)、Squared loss(用于线性回归)、Exponential loss(用于推进法)。其中,Log损失函数采用对数损失函数,具有平滑性、对异常点敏感、可预测概率和适用于大数据实验等优势[21]。
本文中采用logistic损失log(1+e-f)作为损失函数,等同于交叉熵损失函数,用于logistic回归。
2.3.2 算法求解结合logistic的损失函数和上述所提方法,本文中采用Hessian正则化多视角logistic回归算法建模。
将式(6)代入logistic的损失函数中,损失函数的最优化可表示为
$ \log \left( {1 + {{\rm{e}}^{ - {y_i}\mathit{\boldsymbol{K}}\left( {{x_i},x} \right)\mathit{\boldsymbol{\alpha }}}}} \right). $ | (9) |
将式(9)代入式(8)中,则目标函数的最优化可表示为
$ \begin{array}{l} \mathop {\min }\limits_{f \in {H_k},\theta \in {{\bf{R}}^v},\beta \in {{\bf{R}}^v}} \frac{1}{l}\sum\limits_{i = 1}^l {\left( {\log \left( {1 + {{\rm{e}}^{ - {y_i}\mathit{\boldsymbol{K}}\left( {{x_i},x} \right)\mathit{\boldsymbol{\alpha }}}}} \right)} \right)} + {\gamma _k}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\mathit{\boldsymbol{K\alpha }} + \\ {\gamma _I}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\mathit{\boldsymbol{KHK\alpha }} + {\gamma _\theta }\left\| \mathit{\boldsymbol{\theta }} \right\|_2^2 + {\gamma _\beta }\left\| \mathit{\boldsymbol{\beta }} \right\|_2^2. \end{array} $ | (10) |
$ \begin{array}{l} {\rm{s}}.\;{\rm{t}}.\;\sum\nolimits_{k = 1}^V {{\theta ^k}} = 1,{\theta ^k} \ge 0,k = 1,2, \cdots ,V,\\ \sum\nolimits_{j = 1}^V {{\beta ^j}} = 1,{\beta ^j} \ge 0,j = 1,2, \cdots ,V, \end{array} $ |
其中
Hessian正则化多视角logistic回归算法(式(10))求解步骤如下:
(1) 固定θ和β,求解α。固定θ和β后,式(10)可表示为
$ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{\alpha }} \in {{\bf{R}}^{l + u}}} \frac{1}{l}\sum\limits_{i = 1}^l {\left( {\log \left( {1 + {{\rm{e}}^{ - {y_i}\mathit{\boldsymbol{K}}\left( {{x_i},x} \right)\mathit{\boldsymbol{\alpha }}}}} \right)} \right)} + {\gamma _k}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\mathit{\boldsymbol{K\alpha }} + \\ {\gamma _I}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\mathit{\boldsymbol{KHK\alpha }}. \end{array} $ | (11) |
其中
$ \mathit{\boldsymbol{K}} = \sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} ,\mathit{\boldsymbol{H}} = \sum\nolimits_{j = 1}^V {{\beta ^j}{\mathit{\boldsymbol{H}}^j}} . $ |
logistic的损失函数可微,本文中采用梯度下降法求解式(11)。式(11)的一阶导数可表示为
$ \begin{array}{l} \nabla f\left( \mathit{\boldsymbol{\alpha }} \right) = - \frac{{{\rm{loge}}}}{l}\sum\limits_{i = 1}^l {\left( {\frac{{{y_i}}}{{1 + {{\rm{e}}^{{y_i}\mathit{\boldsymbol{K}}\left( {{x_i},x} \right)\mathit{\boldsymbol{\alpha }}}}}}{\mathit{\boldsymbol{K}}^{\rm{T}}}\left( {{x_i},x} \right)} \right)} + {\gamma _k}\left( {\mathit{\boldsymbol{K}} + } \right.\\ \left. {{\mathit{\boldsymbol{K}}^{\rm{T}}}} \right)\mathit{\boldsymbol{\alpha }} + {\gamma _I}\left( {\mathit{\boldsymbol{KHK}} + \mathit{\boldsymbol{KH}}{\mathit{\boldsymbol{K}}^{\rm{T}}}} \right)\mathit{\boldsymbol{\alpha }}. \end{array} $ |
梯度下降算法求解过程如下:
① 初始化α0∈Rl+u,δ,d0=-∇f(α0),0<ε≪1,m=0;
② 当|f(αm+1)-f(αm)|>ε时,
$ {\mathit{\boldsymbol{\alpha }}^{m + 1}} = {\mathit{\boldsymbol{\alpha }}^m} + \delta {d^m}, $ |
$ {d^{m + 1}} = - \nabla f\left( {{\mathit{\boldsymbol{\alpha }}^{m + 1}}} \right) + \frac{{{{\left\| {\nabla f\left( {{\mathit{\boldsymbol{\alpha }}^{m + 1}}} \right)} \right\|}^2}}}{{{{\left\| {\nabla f\left( {{\mathit{\boldsymbol{\alpha }}^m}} \right)} \right\|}^2}}}{d^m}, $ |
$ m = m + 1; $ |
③ α*=αm+1。
(2) 固定α和β,求解θ。固定α和β后,式(10)可表示为
$ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{\alpha }} \in {{\bf{R}}^{l + u}}} \frac{1}{l}\sum\limits_{i = 1}^l {\left( {\log \left( {1 + {{\rm{e}}^{ - {y_i}\left( {\sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} \left( {{x_i},x} \right)} \right)\mathit{\boldsymbol{\alpha }}}}} \right)} \right)} + \\ {\gamma _k}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\left( {\sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} \left( {{x_i},x} \right)} \right)\mathit{\boldsymbol{\alpha }} + \\ {\gamma _I}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\left( {\sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} \left( {{x_i},x} \right)} \right)\mathit{\boldsymbol{H}}\left( {\sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} \left( {{x_i},x} \right)} \right)\mathit{\boldsymbol{\alpha + }}\\ {\gamma _\theta }\left\| \mathit{\boldsymbol{\theta }} \right\|_2^2. \end{array} $ | (12) |
$ {\rm{s}}.\;{\rm{t}}.\;\sum\nolimits_{k = 1}^V {{\theta ^k}} = 1,{\theta ^k} \ge 0,k = 1,2, \cdots ,V. $ |
其中
$ \mathit{\boldsymbol{H}} = \sum\nolimits_{j = 1}^V {{\beta ^j}{\mathit{\boldsymbol{H}}^j}} . $ |
式(12)的一阶导数可表示为
$ \begin{array}{l} \nabla f\left( {{\theta ^k}} \right) = - \frac{{{\rm{loge}}}}{l}\sum\limits_{i = 1}^l {\left( {{\mathit{\boldsymbol{K}}^k}\left( {{x_i},x} \right)\mathit{\boldsymbol{\alpha }}\frac{{{y_i}}}{{1 + {{\rm{e}}^{{y_i}\left( {\sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} \left( {{x_i},x} \right)} \right)\mathit{\boldsymbol{\alpha }}}}}}} \right)} + \\ {\gamma _k}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}{\mathit{\boldsymbol{K}}^k}\mathit{\boldsymbol{\alpha }} + {\gamma _I}{\left( {2\mathit{\boldsymbol{H}}\left( {\sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} \left( {{x_i},x} \right)} \right)\mathit{\boldsymbol{\alpha }}} \right)^{\rm{T}}}{\mathit{\boldsymbol{K}}^k}\mathit{\boldsymbol{\alpha }} + \\ {\gamma _\theta }{\theta ^k}. \end{array} $ |
梯度下降算法求解过程如下:
① 初始化θ0∈Rl+u,δ,0<ε≪1,m=0,dk(0)=-∇f(dk(0)),k=1, 2, …, V;
② 当|f(θm+1)-f(θm)|>ε时,
$ {\theta ^{{k^{m + 1}}}} = {\theta ^{{k^m}}} + \delta {d^{{k^m}}}, $ |
$ {d^{{k^{m + 1}}}} = - \nabla f\left( {{\theta ^{{k^{m + 1}}}}} \right) + \frac{{{{\left\| {\nabla \mathit{\boldsymbol{f}}\left( {{\theta ^{{k^{m + 1}}}}} \right)} \right\|}^2}}}{{{{\left\| {\nabla \mathit{\boldsymbol{f}}\left( {{\theta ^{{k^m}}}} \right)} \right\|}^2}}}{d^{{k^m}}}, $ |
$ m = m + 1; $ |
③ θ*=θm+1。
(3) 固定α和θ,求解β。固定α和θ后,式(10)可表示为
$ \mathop {\min }\limits_{\beta \in {{\bf{R}}^{l + u}}} {\gamma _I}{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}\mathit{\boldsymbol{K}}\left( {\sum\nolimits_{j = 1}^V {{\beta ^j}{\mathit{\boldsymbol{H}}^j}} } \right)\mathit{\boldsymbol{K\alpha }} + {\gamma _\beta }\left\| \mathit{\boldsymbol{\beta }} \right\|_2^2, $ | (13) |
$ {\rm{s}}.\;{\rm{t}}.\;\sum\nolimits_{j = 1}^V {{\beta ^j}} = 1,{\beta ^j} \ge 0,j = 1,2, \cdots ,V, $ |
其中
$ \mathit{\boldsymbol{K}} = \sum\nolimits_{k = 1}^V {{\theta ^k}{\mathit{\boldsymbol{K}}^k}} . $ |
式(13)可看作是多重Hessian最优线性组合的学习过程。
通过上述交替优化算法就可求得式(10)的局部最优解。
具体工况识别过程:对实测地面示功图、实测电功率和实测井口温度信号3个视角原始工况数据进行特征提取得到3个特征集后建立工况样本库并生成训练样本集和测试样本集(训练和测试集中每个工况样本都包含3个特征集)。训练样本分别建立各自的核矩阵和Hessian矩阵,并通过各自的权重系数( θ1+θ2+θ3=1,β1+β2+β3=1)生成各自的多视角核K矩阵和多视角H矩阵,选择合适的模参(γk、γI、γθ、γβ)后通过上述方法训练样本,经过多次迭代收敛后得到最优解集,将这个解集和测试集中测试样本各自的多视角核K矩阵结合就得到测试样本的预测值集,对预测值和实际值进行比较来分类识别抽油机井工况。
3 实验结果分析为验证本文中所提方法及模型的有效性和实用性, 利用胜利油田某区块中近60口抽油机井进行工况识别实验, 区块油藏类型是典型的高压低渗稀油油藏。实验所用抽油机井工况样本库根据油井作业记录建立, 包含11种工况, 每种工况有150个样本, 共1 650个样本。11种典型工况有正常、供液不足、抽油杆断脱、连抽带喷、泵卡、泵游动凡尔失灵、结蜡、油管漏失、泵漏、游动凡尔漏失、固定凡尔漏失。
3.1 全标记训练样本下工况识别结果对比实验时将样本库中的1 650个工况样本平均分成两部分:训练集和测试集各包含825个样本, 其中每种工况各包含75个样本。
采用logistic回归一对多两分类器, 迭代次数为1 200次。测试集中样本总数为825, 其中正例样本数为75, 反例样本数为750。识别准确率的计算方法为测试集中分类正确的样本数占样本总数的比例, 测试样本是按其产生的预测值排序后正例样本排在最前面为截断进行分类[22]。重复做5次实验, 每次实验中训练集和测试集互斥且不一样。
最后, 将本文中所提识别方法(Hessian正则化多视角logistic回归, mHLR)分别与单视角识别方法(logistic回归, LR)、Hessian正则化单视角识别方法(Hessian正则化单视角logistic回归, HLR)、传统特征连接多源识别方法(多特征连接logistic回归, mCLR)对抽油机井工况的识别结果进行对比。其中, 取5次实验结果的平均值作为最终实验结果, 具体如表 2所示。
实测地面示功图、实测电功率和实测井口温度3个视角都采用本文中提取的特征数据。
实验结果分析如下:
(1) 与实测地面示功图识别相比, 本文中所提方法(mHLR)与单视角识别方法(LR、HLR)对抽油机井工况的识别率分别提高约2.4%、2.5%;与实测电功率识别相比, 两种方法的识别率均提高约11%。
(2) 与实测地面示功图和实测电功率两个视角识别相比, 本文中所提方法(mHLR)与传统特征连接多源识别方法(mCLR)对抽油机井工况的识别率提高约11.8%;与实测地面示功图、实测电功率和实测井口温度3个视角识别相比, 两种方法的识别率提高约13.8%。本文中所提方法的识别结果明显优于传统特征连接多源方法的识别结果, 并且在选择合适的特征数据下, 本文中所提方法会随着特征数据的增加识别精准率提高, 而传统特征连接多源识别方法会随着特征数据的增加识别精准率下降。
(3) 分别利用实测地面示功图信号、实测地面示功图和实测电功率信号、实测地面示功图、实测电功率和实测井口温度信号, 采用不同方法对抽油机井工况进行识别, 都可以取得较好的识别结果。
(4) 本文实验中Hessian正则化和单视角结合对识别率影响不大, 但和多视角结合可以提高识别率。
3.2 不同标记训练样本下工况识别结果对比在第一部分实验内容基础上, 分别从单视角和多视角中选取识别率高的识别方法进行实验。将本文中所提识别方法(mHLR)和单视角识别方法(HLR和LR)按5组不同比例标记训练样本对抽油机井工况的识别结果进行对比。每组各重复做5次实验, 分别取其5次实验结果的平均值作为最终实验结果, 具体如表 3所示。
实验结果分析如下:
(1) 受地质、设备等因素影响, 抽油机井工况复杂多变, 工况样本不易获取, 但在油气生产物联网环境下抽油机井可以获取大量未知工况样本。与其他方法相比, 本文中所提方法(mHLR)可以实现在少量已知工况样本下, 充分利用大量未知工况样本来实现对抽油机井工况更精准的识别, 该方法更符合实际工程应用。
(2) 在不同比例标记训练样本下, 本文中所提识别方法(mHLR)对抽油机井工况的识别结果都比单视角识别方法(HLR和LR)对抽油机井工况的识别结果准确率高, 尤其是在少量(低于30%)标记训练样本情况下效果更明显, 识别率平均提高约3%。
(3) 在低于30%少量标记训练样本下, 本文中所提识别方法(mHLR)在识别率的提高幅度上高于单视角识别方法(HLR和LR)。
4 结论(1) 提出一种基于Hessian正则化多视角学习的抽油机井工况识别新方法, 该方法通过集成多视角多核学习算法和多视角图集成Hessian学习算法, 大大提高了算法的泛化能力。实验结果表明该方法比单一信息源识别方法、传统特征连接多源识别方法识别精度高。
(2) 利用实测地面示功图、实测电功率和实测井口温度信号, 结合先验知识、专家经验和机制分析进行特征提取, 并将log损失函数引入Hessian正则化多视角学习方法中, 建立Hessian正则化多视角logistic回归抽油机井工况识别模型, 提高了模型的鲁棒性, 且在全标记和少量标记训练样本下都取得了较好的识别效果, 工程实用性强。
(3) Hessian正则化多视角学习方法还可以推广到螺杆泵、无杆泵采油井以及其他领域的故障诊断与识别方面, 具有较好的推广价值。
[1] |
韩国庆, 吴晓东, 毛凤英, 等. 示功图识别技术在有杆泵工况诊断中的应用[J]. 石油钻采工艺, 2003, 25(5): 70-74. HAN Guoqing, WU Xiaodong, MAO Fengying, et al. Application of dynamometer card identification in diagnosis of working condition for suck rod pump[J]. Oil Drilling & Production Technology, 2003, 25(5): 70-74. |
[2] |
陈培毅. 基于抽油机实测电功率的悬点示功图仿真与工况诊断[D]. 秦皇岛: 燕山大学, 2013. CHEN Peiyi. Simulation of dynamometer card and working condition diagnosis based on pumping unit's measured electric power[D]. Qinhuangdao: Yanshan University, 2013. http://cdmd.cnki.com.cn/Article/CDMD-10216-1013029145.htm |
[3] |
LI K, GAO X W, TIAN Z D, et al. Using the curve moment and the PSO-SVM method to diagnose downhole conditions of a sucker rod pumping unit[J]. Petroleum Science, 2013, 10(1): 73-80. DOI:10.1007/s12182-013-0252-y |
[4] |
WU W, SUN W L, WEI H X. A fault diagnosis of suck rod pumping system based on wavelet packet and RBF network[J]. Advanced Materials Research, 2011(189/190/191/192/193): 2665-2669. |
[5] |
XU P, XU S J, YIN H W. Application of self-organizing competitive neural network in fault diagnosis of suck rod pumping system[J]. Journal of Petroleum Science & Engineering, 2007, 58(1/2): 43-48. |
[6] |
梁华. 有杆抽油系统故障递阶诊断的故障识别研究[J]. 西南石油大学学报(自然科学版), 2015, 37(1): 165-171. LIANG Hua. Hierarchical fault diagnosis of rod pumping system based on fault distinguishing[J]. Journal of Southwest Petroleum University(Science & Technology Edition), 2015, 37(1): 165-171. DOI:10.11885/j.issn.1674-5086.2012.11.19.04 |
[7] |
孙振华. 游梁式抽油机采油系统实时评价方法研究[D]. 青岛: 中国石油大学(华东), 2011. SUN Zhenhua. The real-time evaluation method of the beam-pumping recovery system[D]. Qingdao: China University of Petroleum(East China), 2011. http://d.wanfangdata.com.cn/Thesis/Y1877054 |
[8] |
陈黎兴. 基于电功率数据的抽油机故障诊断方法研究[D]. 青岛: 中国石油大学(华东), 2016. CHEN Lixing. Electrical power data based fault diagnosis of oil pumping units[D]. Qingdao: China University of Petroleum(East China), 2016. |
[9] |
王凯. 基于产生式规则系统的抽油泵故障诊断[J]. 石油勘探与开发, 2010, 37(1): 116-120. WANG Kai. Fault diagnosis of rod-pumping unit based on production rules system[J]. Petroleum Exploration and Development, 2010, 37(1): 116-120. |
[10] |
LIU S, RAGHAVENDR C S, LIU Y, et al. Automatic early fault detection for rod pump systems[R]. SPE 146038, 2011.
|
[11] |
王凯, 刘宏昭, 熊俊, 等. 基于改进的超球支持向量机的有杆抽油泵故障诊断研究[J]. 机械科学与技术, 2011, 30(1): 133-141. WANG Kai, LIU Hongzhao, XIONG Jun, et al. Fault diagnosis of a rod-pumping unit based on improved hyper sphere vector machines[J]. Mechanical Science and Technology for Aerospace Engineering, 2011, 30(1): 133-141. |
[12] |
VEERARAGHAVAN A, CHELLAPPA R, ROY-CHO-WDHURY A K. The function space of an activity[C/OL]//2006 IEEE Computer Society Conference on. IEEE, 2006: computer vision and pattern recognition[2017-03-20]. http://ieeexplore.ieee.org/abstract/document/1640855.
|
[13] |
IOSIFIDIS A, TEFAS A, PITAS I. Multi-view action recognition based on action volumes, fuzzy distances and cluster discriminant analysis[J]. Signal Processing, 2013, 93(6): 1445-1457. DOI:10.1016/j.sigpro.2012.08.015 |
[14] |
LIU W F, LIU H L, TAO D C, et al. Multiview Hessian regularized logistic regression for action recognition[J]. Signal Processing, 2015, 110: 101-107. DOI:10.1016/j.sigpro.2014.08.002 |
[15] |
LIU W F, TAO D C. Multiview Hessian regularization for image annotation[J]. IEEE Transactions on Image Processing, 2013, 22(7): 2676-2687. DOI:10.1109/TIP.2013.2255302 |
[16] |
DONOHO D L, GRIMES C. Hessian eigenmaps:locally linear embedding techniques for high-dimensional data[J]. Proceedings of the National Academy of Sciences, 2003, 100(10): 5591-5596. DOI:10.1073/pnas.1031596100 |
[17] |
TAO D C, JIN L, LIU W F, et al. Hessian regularized support vector machines for mobile image annotation on the cloud[J]. IEEE Transactions on Multimedia, 2013, 15(4): 833-844. DOI:10.1109/TMM.2013.2238909 |
[18] |
SCHÖLKOPF B, HERBRICH R, SMOLA A J. A generalized representer theorem[C/OL]//Computational Learning Theory. Berlin/Heidelberg: Springer, 2001: 416-426[2017-03-20]. http://www.umiacs.umd.edu/~hal/tmp/P139.pdf.
|
[19] |
BEZDEK J C, HATHAWAY R J. Convergence of alternating optimization[J]. Neural Parallel & Scientific Computations, 2003, 11(4): 351-368. |
[20] |
胡广杰, 易斌, 田宝库, 等. 抽油机井实测示功图泵况诊断分析[M]. 北京: 石油工业出版社, 2008, 1-45.
|
[21] |
刘红丽. 多视角学习在智能信息处理中的若干应用研究[D]. 青岛: 中国石油大学(华东), 2015. LIU Hongli. Study on multi-view learning in intelligent information processing[D]. Qingdao: China University of Petroleum(East China), 2015. |
[22] |
周志华. 机器学习[M]. 北京: 清华大学出版社, 2016, 29-35.
|