基于在线动态高斯过程回归抽油井动液面软测量建模

李翔宇，高宪文，侯延彬

（东北大学信息科学与工程学院，辽宁 沈阳 110819）

引 言

抽油井动液面深度是反映油层供液能力和井下供排关系的重要工况参数，是确定油井工作制度、预测产量和故障诊断的重要依据。目前，油田现场主要使用回声仪探测抽油井动液面深度，其主要缺陷是：（1）由人工定期操作，无法实时在线检测；（2）由于油套环空气体主要成分为甲烷，井口弹药发声装置存在一定安全隐患[1]。软测量方法[2-9]为抽油井动液面实时在线检测提供了一种新思路。

文献[1]以有杆抽油系统地面（光杆或称悬点）示功图和井口套压等数据作为辅助变量，采用机理方法建立动液面软测量模型。但有杆抽油系统具有强非线性、机电液耦合和井下工况不确定等综合复杂特性，难以精确建立动液面机理模型，动液面计算误差大。文献[10]针对潜油柱塞泵抽油井提出了一种基于支持向量机的沉没度（泵挂深度与动液面深度的差值）预测方法。但随时间推移，油井工况会发生不同程度的变化，使用原始数据建立的软测量模型将不再适应新的工况而出现模型老化现象，无法准确预测当前沉没度信息。

高斯过程回归（Gaussian process regression，GPR）是基于贝叶斯统计学习理论发展起来的一种全新机器学习方法，对处理高维数、小样本、非线性等复杂问题具有很好的适应性，与神经网络、支持向量机等方法相比，GPR 具有模型参数明显减少、参数优化相对容易、输出具有概率意义等优 点[11-15]。基于此，并考虑到采油过程的慢时变动态特性，本文采用一种增量学习动态高斯过程回归（incremental dynamic Gaussian process regression）软测量建模方法建立抽油井动液面软测量模型，通过一种增量学习算法对模型进行在线更新，以使其适应采油过程中油井工况变化，从而提高软测量模型的预测精度和泛化能力。以期为抽油井动液面实时在线测量提供一种有效手段。

1 有杆泵采油过程描述

典型有杆抽油系统如图1所示，皮带轮、减速箱和曲柄连杆机构将电动机的高速转动变为游梁的上下摆动，挂在驴头上的悬绳器通过抽油杆柱带动抽油泵柱塞作上下往复运动，将井内液体抽汲至地面[16]。抽油泵是有杆抽油系统的井下设备，主要由工作筒、柱塞、游动阀和固定阀组成。上冲程时，抽油杆柱带着柱塞向上运动，柱塞上的游动阀受油管内液柱压力而关闭，泵内压力降低，固定阀在油套环空内气、液柱压力（泵口压力或称沉没压力）与泵内压力之差的作用下被打开，泵内吸入液体。此时，如果油管内已充满液体，在井口将排出相当于柱塞冲程长度的一段液体。下冲程时，抽油杆柱带着柱塞向下运动，固定阀关闭，泵内压力增高到大于柱塞上液柱压力时，游动阀被顶开，泵向油管内排液。由于有相当于冲程长度的一段光杆从井外进入油管，在井口将排挤出相当于这段光杆体积的液体[17]。

抽油井动液面是采油过程中油气从地层向井筒渗流与深井泵向地面排液相互作用的宏观表征。当地层压力一定时，动液面深度越大，油井生产压差越大，油井产量越高。当泵的抽汲参数（冲程、冲次、泵径）一定时，泵的理论排量一定，此时，动液面深度越小则表明地层压力越大，油井供液越充足，油藏潜力并未得到充分发挥。在微观层面上，动液面位置直接影响采油过程中流体进泵运动规律，从而影响泵的充满程度和泵效，并在地面示功图形状特征上有所反映。上冲程中，在沉没压力的作用下，井内液体克服固定阀阻力进入泵内，作用在柱塞底部而产生向上的载荷，其大小与动液面位置有关，是上冲程悬点静载荷的重要组成部分，直接影响上下冲程悬点静载荷之差。此外，由于动液面位置与采油过程中井内液体的有效举升高度关系密切，因此，动液面位置对反映抽油井做功情况的地面示功图面积有直接影响。而当动液面深度过大时，泵的沉没度过小，沉没压力过小，加之由于井底流动压力降低造成原油中所溶解的天然气大量析出，井内液体中游离气含量增大，易导致泵的充满程度不足，泵效偏低和“液击”现象，此时，地面示功图通常呈“刀把”形状。

2 动态GPR 模型

2.1 高斯过程回归（GPR）

图1 有杆抽油系统示意图Fig.1 Sucker-rod pumping systems

协方差函数的确定是建立高斯过程模型的关键步骤，它包含了对期望函数特性的假设。对任意一组输入，协方差函数应满足其产生的协方差矩阵为对称半正定。考虑到系统的平稳性，下列径向基函数是最常用的一类协方差函数

式中，υ0表示服从高斯分布的噪声方差；υ1表示局部相关性的程度；wd为模型的测度参数；δij是Kronecker 算子。

当上述的协方差函数类型确定后，通常使用极大似然、交叉验证和马尔可夫链蒙特卡罗3 种方法对其超参数进行调整，本文通过极大化对数似然函数的方法求得

优化过程中，要计算对数似然函数对各参数的导数

2.2 动态系统建模

高斯过程回归除了可用于建立静态非线性映射，基于具有外部输入自回归模型结构的高斯过程回归还可用于动态系统建模[14]，系统输出 y(k) 与过去l个输出和过去l个输入服从如下非线性关系

式中，f(⋅)为非线性函数，本文利用高斯过程回归方法对其进行拟合；ε(k)为白噪声。

3 在线学习算法

增量学习是在线学习的一个重要研究方向，是在原有模型基础上，充分利用历史训练结果，通过增加新样本对模型进行再学习的一种方法。一种增量式在线学习算法[18]是基于具有相关子样本迭代结构的贝叶斯在线方法[19]的线性组合。这种近似是通过参数化和投影技术获得的，均值与协方差的近似值由式（7）迭代计算

式中，比例系数qk+1和rk+1分别是对数似然的一、二阶导数

将式(9)和式(10)展开可得高斯过程后验概率的近似值为

式中，kx=[C(x1,x),…,C(x k,x)]T，系数α和Γ由式(13)～式(15)计算

采用这种递归形式，系数的规模会随着新数据的到来不断增加。为了保证系数的规模固定，采用两种方式实现GP 模型的更新。基础更新被执行，当由式(16)定义的近似误差

小于定义的阈值时，否则，执行完全更新。基础更新仅更新系数而不改变其规模，通过在式(15)中使用k维单位向量。完全更新则通过在式(15)中使用，并把现在的数据添加到最相关数据的子集中。如果此操作导致子集规模超出最大限度，则相关性最小的数据被移除，其使得式(17)值最小

式中，Q是K的逆矩阵，可由式(18)迭代计算

假设相关性最小的数据被移除，系数更新过程如下

4 抽油井动液面软测量建模

4.1 辅助变量选取

辅助变量的选择一般从被测过程的工艺机理入手，定性地找出对预测模型输出贡献较大的因素作为辅助变量。通过对有杆泵采油过程机理的深入分析，选取上下冲程悬点静载荷差wl′、地面示功图面积ad、抽油泵实际排量q、理论排量qt、生产气液比rgo、含水率nw、油压po、套压pc的历史时刻值以及动液面深度ld的历史时刻值作为辅助变量。其中，上下冲程静载荷之差和地面示功图面积是依据物理意义对示功图进行特征参数提取后计算得到的。

4.2 样本数据采集与处理

考虑到实际采油过程的慢时变动态特性，日常生产时作业区每日对油井的生产数据统计一次，因此，在本文中以日为周期对样本数据进行采集。对于软测量模型的主导变量动液面数据由作业区利用回声仪每天测量一次，采样时做多次测量取平均值处理，对于明显异常的数据，由相邻数据取平均值代替[10]。前述软测量模型的辅助变量则是由井口配备的示功图采集器、翻斗式计量器、压力传感器等测井仪器和计量装置实现自动在线检测，并通过嵌入式工控机无线传输至上位机形成数据库。

值得注意的是，上述数据间数量级相差较大，需要对输入数据进行归一化处理，避免训练时算法侧重修改数值大的变量所对应的权值，从而主导小数值变量的情况，同时可以避免计算过程中的“数值灾难”[10,21]。因此，在模型训练之前需将样本集映射到[-1,1]内。

4.3 示功图特征参数提取

由于实测地面示功图中除静载荷外，还包含惯性载荷、振动载荷和摩擦载荷，因此，实测示功图常奇形怪状，各不相同。而通过对示功图形成机理的定性分析与定量计算，可以从中提取出与动液面信息关系密切的上下冲程悬点静载荷之差wl′和地面示功图面积ad这两个特征参数。以架5-3 井为例，图2为该油井2012年9月15日实测地面示功图，示功图由250 组悬点载荷-位移数据对构成，传感器采样周期为抽油机运动周期的1/250。图3为与上述示功图对应的，悬点、抽油泵柱塞和泵筒内液面位置变化曲线。

图2示出，架5-3 井存在井下供液不足，泵的充满程度较低的现象，其示功图呈“刀把”状，开采中后期的油井普遍长期处于这种工况下运行。针对此类示功图，可以按图2所示“两点八段”原则进行悬点静载荷分析与特征参数提取。

图2 架5-3 井实测地面示功图Fig.2 Measured surface dynamometer card of J5-3

（1）冲程下死点区：如图2所示，此区域为上、下冲程换向区，摩擦力方向发生变化，悬点静载荷PRL 为

式中，Wr′为抽油杆在井液中重力；F为摩擦载荷。

（2）悬点加载区：在此区域内，固定阀仍处于关闭状态，悬点载荷随抽油杆柱的拉伸而增大。

（3）托浮力下降区：如图2所示，示功图上特征点A为固定阀开启点。固定阀开启后，井内液体克服固定阀阻力进入泵内，随着流量增大，固定阀阻力上升，泵腔内压力降低，作用于柱塞底部向上的托浮力减小，悬点静载荷为

式中，Wr抽油杆柱在空气中的重力；Wl为作用于柱塞上的液柱载荷；Pi为泵腔内压力作用于柱塞底部产生向上的载荷。

（4）上冲程柱塞离油区：如图3所示，在沉没度过低和气体影响等因素作用下，进泵液体运动滞后于柱塞运动，即柱塞脱离于泵筒内液面，泵腔内压力作用于柱塞底部向上的托浮力可忽略，悬点静载荷可按下式计算

图3 架5-3 井悬点、柱塞和泵筒内液面位置变化Fig.3 Position of polished rod、plunger and liquid level in chamber of J5-3

（5）冲程上死点区：与过程（1）类似，该区域亦为上、下冲程换向区，摩擦力方向发生变化，悬点静载荷为

（6）下冲程悬点未卸载区：在此区域，由于游动阀尚未开启，悬点不能卸载，相比过程(4)，仅摩擦力换向，悬点静载荷为

（7）下冲程悬点卸载区(柱塞入油)：如图2所示，示功图上特征点B为游动阀开启点。当柱塞下行遇到液面，游动阀开启，抽油杆柱收缩，悬点载荷迅速减小，泵筒开始向油管排液。

（8）游动阀阻力下降区：排油过程即将结束，经过游动阀液体流量减少，游动阀阻力下降，悬点载荷增大。

根据上述分析，架5-3 井上下冲程中悬点静载荷差wl′ (k)可由式(27)计算

式中，PRLi为上死点区内n1个悬点载荷，PRLj为下死点区内n2个悬点载荷。示功图面积ad(k)则可通过悬点载荷对位移的积分得到。

4.4 软测量模型的建立

在线动态高斯过程回归软测量建模步骤如下。

（1）样本数据归一化后，经过稳健估计去除离群点[2]，随机误差噪声采用均值滤波算法处理[3]。

（2）针对部分样本数据，采用动态GPR 方法建立基本软测量模型，协方差函数如式(3)所示。软测量模型超参数Θ中w1～w9的初始值均为1，υ1和υ0分别取1 和0.001，即假设每个辅助变量对输出预测的贡献相同，使用共轭梯度法，通过式(5)搜索得到与输入对应的超参数最优值。根据油田现场采油工艺专家经验，以油井前一天相关生产参数作为软测量模型的输入，即式(6)中取l=1，则有

（3）随着时间推移，当模型不再适应新的工况而出现预测误差增大现象时，采用增量学习算法对软测量模型进行在线更新。对于所有新数据(xk+1,yk+1)，计算标量qk+1、rk+1，若ξk+1小于定义的阈值，确定矢量kk+1、sk+1，进而更新GP 参数(αk+1,Γk+1)；若ξk+1大于等于定义的阈值，由式(17)确定相关性最小的数据，按式(19)、式(20)更新GP参数。从上述流程可知，IDGPR 针对每个新增样本均能自适应地进行模型更新，并建立相应的软测量模型。

5 工业应用

2011年10月以来，在东北某油田对架5-3 等多口油井进行了旨在实现数字化油田的，包括抽油井关键生产参数实时检测与有杆抽油系统故障检测与诊断的集散监控管理系统建设，通过安装在井口的无线传感器及抽油井智能测控终端设备同步采集了大量的油压、套压、产液量和光杆示功图等反映抽油井工况信息的样本数据。在现场同步测试的基础上，应用上述方法进行了抽油井动液面软测量建模研究，现以2013年架5-3 井为例进行分析。

架5-3 井机械装置及有关参数如下：抽油机型号CYJ10-3-53HB，电机型号Y250M-6，光杆冲程3 m，泵径44 mm，泵深1801.15 m，油管内/外径62/73 mm，抽油杆柱组合19 mm×994.25 m＋22 mm×498.6 m＋25 mm×292.47 m。正常生产时，该油井套管闸门关闭，气体全部通过泵抽出。此外，该井配备有变频调速装置，可以方便地实现抽油机冲次无级调整。

以平均相对误差、均方根误差作为软测量模型性能评价标准

式中，N为实际测试样本数。其中，平均相对误差与均方根误差越小，模型的预测精度越高。

选取架5-3 井全年365 组数据，经过数据预处理后，前100 组数据用于批量学习建立基本动态高斯过程回归(DGPR)软测量模型，从第101 组数据开始对模型进行在线更新。为了考察增量学习动态高斯过程回归(IDGPR)软测量模型对油井生产工况和过程参数变化的在线更新适应能力，在2013年6月8日至6月25日期间对有杆泵抽汲参数进行调整，使有杆抽油系统工况发生变化，抽油机冲次在6月8日由4 次改为3 次，在6月18日又由3 次改为4.5 次。图4给出了GPR 和IDGPR 软测量模型对测试数据进行预测的结果，图5为软测量值与真值的相对误差曲线。从图4、图5可以看出，GPR模型的预测曲线已偏离真值曲线，预测相对误差较大，模型已经开始失效，而IDGPR 模型能够适应工况的变化，在整个测试阶段模型预测输出更接近真实值，具有更好的拟合效果。为进一步验证IDGPR 模型性能，以2013年8月1日至2013年11 月8日100 组数据进行测试，图6、图7分别给出了IDGPR 模型和GPR 模型对验证数据的预测输出和相对误差。其中，图4、图6中的IDGPR±2σ是IDGPR 模型预测值的95%置信概率对应的两倍标准差置信区间曲线。

图4 短期冲次改变动液面预测效果比较Fig.4 Comparison between predictive output and real value during changing pump cycle

图5 短期冲次改变动液面预测相对误差Fig.5 Predictive error during changing pump cycle

图6 长期动液面预测效果比较Fig.6 Comparison between predictive output and real value for a long time

图7 长期动液面预测相对误差Fig.7 Predictive error for a long time

表1示出了在相同测试数据集下进行动液面深度预测，上述两种软测量模型的性能指标。可以看出，使用增量学习动态高斯过程软测量模型可以获得更小的平均相对误差、均方根误差，进一步说明该软测量建模方法预测精度较高、不确定度较小。架5-3 井日数据报表与动液面预测对比见表2。

表1 架5-3 井动液面预测误差分析Table 1 Predictive error analysis of J5-3 well

6 结 论

针对现有回声仪测试法、光杆示功图法和软测量建模法在确定抽油井动液面深度方面存在的不足，采用一种增量学习动态高斯过程回归软测量建模方法，实现动液面实时在线检测。通过引入具有外部输入的自回归建模策略，构造动态高斯过程软测量模型，并采用一种增量学习算法使原始模型能够根据油井工况变化不断进行在线更新，自适应获得更加准确的软测量模型。现场应用结果表明IDGPR 软测量建模方法能够对抽油井动液面深度给出较准确的实时在线预测，模型泛化能力强，可以满足工程应用要求。

表2 架5-3 井日数据表与动液面预测对比Table 2 Dynamic production data and dynamic liquid level prediction of J5-3 well

符 号 说 明

ad——地面示功图面积，m2

C(xi,xj)——协方差函数

D——训练样本集

ek+1——k+1 阶单位向量

F——摩擦载荷，kN

f(·)——非线性函数

K——训练样本输入间的n×n协方差矩阵

k(x∗)——测试输入和训练样本输入间n×1协方差向量

L(Θ)——超参数对数似然函数

ld——动液面深度，m

nw——原油积含水率，%

PRL——悬点载荷，kN

pc——井口套管压力，MPa

pi——吸入压力作用在柱塞底部产生的载荷，kN

po——井口油管压力，MPa

Q——K的逆矩阵

q——抽油泵实际排量，m3·d-1

qt——抽油泵理论排量，m3·d-1

qk+1——对数似然函数的一阶导数

rgo——生产气液体积比，m3·m-3

rk+1——对数似然函数的二阶导数

T k+1——把k维向量扩展成k+1 维的算子，通过在向量最后一行加零的方式

U k+1——把k维矩阵扩展成k+1 维的算子，通过在矩阵最后一列添加零的方式

Wl——作用在柱塞上的液柱载荷，kN

Wr——抽油杆柱在空气中的重力，kN

lw′——上下冲程悬点静载荷差，kN

wd——模型的测度参数

yi——第i个训练输出

y——训练输出构成n×1 的向量

xi——第i个训练输入

Z=E{P(y k+1|fk+1)}k——完备数据似然

α——预测均值的迭代更新系数

Γ——预测方差的迭代更新系数

δij——Kronecker 算子

ε(k)——白噪声

Θ——协方差函数超参数

μ——预测均值

ξk+1——模型更新近似误差

σ2——预测方差

υ0——白噪声方差

υ1——局部相关性的程度

[1]Zhang Shengli (张胜利),Luo Yi (罗毅),Wu Zanmei (吴赞美),et al.Corrected algorithm for calculating dynamic fluid level with indicator diagram for rob-pumped well [J].Oil Drilling & Production Technology(石油钻采工艺),2011,33 (6):122-124.

[2]Li Haibo (李海波),Chai Tianyou (柴天佑),Yue Heng (岳恒).Soft sensor of technical indices based on KPCA-ELM and application for flotation process [J].CⅠESC Journal(化工学报),2012,63 (9):2892-2898.

[3]Li Lijuan (李丽娟),Pan Lei (潘磊),Zhang Shi (张湜).Soft sensor modeling for mobility of jig bed based on AP-clustering algorithm [J].CⅠESC Journal(化工学报),2012,63 (9):2675-2680.

[4]Wang Bo (王博),Sun Yukun (孙玉坤),Ji Xiaofu (嵇小辅),et al.Soft-sensor modeling for lysine fermentation processes based on PSO-SVM inversion [J].CⅠESC Journal(化工学报),2012,63 (9):3000-3007.

[5]Zhong Weimin (钟伟民),Li Jie (李杰),Cheng Hui (程辉),et al.A soft sensor multi-modeling for furnace temperature of gasifier based FCM clustering [J].CⅠESC Journal(化工学报),2012,63 (12):3951-3955.

[6]Yang Xiaomei (杨小梅),Liu Wenqi (刘文琦),Yang Jun (杨俊).LSSVM modeling for fermentation process based on dividing stages [J].CⅠESC Journal(化工学报),2013,64 (9):3262-3269.

[7]Wei Yujie (魏宇杰),Shang Chao (尚超),Gao Xinqing (高莘青),et al.Dynamic process based soft sensing of melt flow rate [J].CⅠESC Journal(化工学报),2014,65 (8):3062-3070.

[8]Li Jun (李军),Yue Wenqi (岳文琦).Dynamic soft sensor modeling and its application using leaky-integrator ESN [J].CⅠESC Journal(化工学报),2014,65 (10):4004-4014.

[9]Wang Tong (王通),Gao Xianwen (高宪文),Liu Wenfang (刘文芳).Adaptive soft sensor method and application in determination of dynamic fluid levels [J].CⅠESC Journal(化工学报),2014,65 (12):4898-4904.

[10]Yu Deliang (于德亮),Qi Weigui (齐维贵),Deng Shengchuan (邓盛川),et al.Submergence forecasting of a submersible plunger pump based on the support vector machine [J].Acta Petrolei Sinica(石油学报),2011,32 (3):534-538.

[11]Wang Huazhong (王华忠).Gaussian process and its application to soft-sensor modeling [J].Journalof Chemical Ⅰndustry and Engineering(China) (化工学报),2007,58 (11):2839-2845.

[12]Cao Pengfei (曹鹏飞),Luo Xionglin (罗雄麟).Modeling of soft sensor for chemical process [J].CⅠESC Journal(化工学报),2013,64 (3):788-800.

[13]Lei Yu (雷瑜),Yang Huizhong (杨慧中).Combination model soft sensor based on Gaussian process and Bayesian committee machine [J].CⅠESC Journal(化工学报),2013,64 (12):4434-4438.

[14]He Zhikun (何志昆),Liu Guangbin (刘光斌),Zhao Xijing (赵曦晶),et al.Overview of Gaussion process regression [J].Control andDecision(控制与决策),2013,28 (8):1121-1129.

[15]Yu Tao (于涛),Wang Jianlin (王建林),He Kun (何坤),et al.Staged soft-sensor modeling method for fermentation process based on MPCA-GP [J].Chinese Journal of Scientific Ⅰnstrument(仪器仪表学报),2013,34 (12):2703-2708.

[16]Wan Renpu (万仁溥).Production Technology Manual (采油技术手册:第四分册) [M].Beijing:Petroleum Industry Press,1993:1.

[17]Zhang Qi (张琪).Production Engineering and Design (采油工程原理与设计) [M].Dongying:Press of China University of Petroleum,2000:98-99.

[18]Petelin D,Kocijan J.Control system with evolving Gaussian process models//2011 IEEE Workshop on Evolving and Adaptive Intelligent Systems (EAISs) [C].2011:178-184.

[19]Csato L.Gaussian process-iterative sparse approximations [D].Birmingham:Aston University,2002.

[20]Opper M.A Bayesian Approach to Online Learning [M].Cambridge:Cambridge University Press,1998:363-378.

[21]Zhou Dawei (周大为),Gao Xiang (高翔),Xia Changgao (夏长高).Soft sensor method for estimating engine minimum fuel consumption of based on chaos optimization and SVM [J].Chinese Journal of Scientific Ⅰnstrument(仪器仪表学报),2011,32 (2):463-468.