基于在線動態高斯過程回歸抽油井動液面軟測量建模

李翔宇，高憲文，侯延彬

（東北大學信息科學與工程學院，遼寧 沈陽 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.