基于ARIMA-Kalman濾波器數(shù)據(jù)挖掘模型的油井產(chǎn)量預測

谷建偉，隋顧磊，李志濤，劉 巍，王依科，張以根，崔文富

1)中國石油大學(華東)石油工程學院，山東青島 266580；2)中國石化勝利油田分公司勘探開發(fā)研究院，山東東營 257015；3)中國石化勝利油田分公司勝利采油廠，山東東營 257015

油田產(chǎn)量是反映油田開發(fā)效果的重要指標，是掌握油田動態(tài)變化的重要依據(jù)．在不同的開發(fā)階段，影響產(chǎn)量變化的因素不同，但始終遵從地下滲流規(guī)律，且這些影響因素之間存在各種關聯(lián)．目前，描述油田產(chǎn)量變化的方法有2種：一是基于基本的滲流理論油藏工程類方法[1]，如產(chǎn)量遞減分析[2]．該方法是預測和分析油藏動態(tài)的常用數(shù)理統(tǒng)計方法，也是油藏工程類的典型代表之一，適用于產(chǎn)量遞減階段的油田．由于石油生產(chǎn)不斷受到各種技術的干預，如采取壓裂和酸化等增產(chǎn)增注措施的油氣井，產(chǎn)量遞減分析具有一定的局限性．油藏工程類方法考慮了儲層性質(zhì)、井況和生產(chǎn)控制參數(shù)對產(chǎn)油量造成的影響，但由于目前的滲流理論是在理想滲流環(huán)境下得到的，不能完全反應實際油田滲流的現(xiàn)象和規(guī)律[3]．二是基于數(shù)據(jù)挖掘的機器學習類方法[4-6]．EDIGER[7]采用差分自回歸積分移動平均(autoregressive integrated moving average, ARIMA)模型預測土耳其石油產(chǎn)量，證明ARIMA回歸和置信區(qū)間的吻合度越高越能提高產(chǎn)量預測模型的精度和可靠性．王濱等[8]應用時間序列傳遞函數(shù)模型建立了考慮因素動態(tài)關系的多因素油田產(chǎn)油量數(shù)據(jù)擬合模型，但并未用于產(chǎn)油量預測．FRAUSTO-SOLS等[9]通過應用ARIMA、NARX神經(jīng)網(wǎng)絡、單指數(shù)平滑以及雙指數(shù)平滑等模型預測石油產(chǎn)量，證明了ARIMA模型預測精度良好．但是，ARIMA模型雖然具有高效的時序影響分析能力，卻有一定滯后性[9]．

本研究以單口生產(chǎn)井的產(chǎn)油量為時間序列，根據(jù)該油井的歷史產(chǎn)油量數(shù)據(jù)建立時間序列中的產(chǎn)油量ARIMA模型，并結(jié)合卡爾曼濾波器(Kalman filter)[10]，構(gòu)建基于ARIMA-Kalman濾波器[11]的產(chǎn)量預測模型．ARIMA-Kalman濾波器具有高效的時序影響因素的分析能力，能夠排除非同步性以及滯后性的影響，使識別出的產(chǎn)油量時間序列模型具有精準的擬合結(jié)果和預測能力，并縮短了滯后時間．

1 ARIMA-Kalman濾波器模型原理

ARIMA模型是將非平穩(wěn)時間序列轉(zhuǎn)化為平穩(wěn)時間序列，然后將因變量僅對其滯后值及隨機誤差項的現(xiàn)值和滯后值進行回歸所建立的模型[12]．在如式(1)的差分自回歸移動平均模型ARIMA(p,d,q)中，p為自回歸(auto regressive, AR)多項式的階數(shù)，q為移動平均(moving average, MA)多項式的階數(shù)，d為時間序列成為平穩(wěn)時所做的差分次數(shù)．ARIMA模型根據(jù)原序列是否平穩(wěn)以及回歸中所含部分的不同，包括移動平均過程、自回歸過程、自回歸移動平均過程以及ARIMA過程[12]．如果研究的時間序列Qt是非平穩(wěn)的，則可在標準的ARIMA模型中通過適當?shù)牟罘肢@得平穩(wěn)的時間序列．

(1)

其中，φ(B)為AR多項式；θ(B)為MA多項式；t為時間序列下標，t=1, 2, …；d為差分次數(shù)；e(t)為正常的白噪聲序列，均值為0，方差為δ2；φi(t)為AR參數(shù)，i=1, 2, …；θj(t)為MA參數(shù)，j=1, 2, …．

平穩(wěn)的時間序列用于建立ARIMA模型并確定其參數(shù)．參數(shù)(p,q)通過自相關函數(shù)(auto correlation function, ACF)和偏自相關函數(shù)(partial autocorrelation function, PACF)預先確定，根據(jù)AIC信息準則(Akaike information criterion)[13]最終確定．在時間序列分析中，Yule-Walker方程在模式識別和參數(shù)估計中起著重要作用[14]．ARIMA模型的殘差可以作為模型估計準則，模型估計過程采用最大似然估計．通過求解殘差序列的ACF，可以判斷殘差序列是否為白噪聲．如果殘差序列不是白噪聲，則重新定義(p,q)．

因此，t+1時刻預測的產(chǎn)油量可描述為

Q(t+1)=φ1(t)Q(t)+φ2(t)Q(t-1) +…+

φp(t)Q(t-p+1)+e(t+1)-

θ1(t)e(t)-θ2(t)e(t-1)-…-

θq(t)e(t-q+1)

(2)

其中，Q(t+1),Q(t),Q(t-1), …,Q(t-p+1)為產(chǎn)油量時間序列；e(t+1),e(t),e(t-1), …,e(t-q+1)為殘差時間序列．

ARIMA模型將預測對象隨時間推移而形成的數(shù)據(jù)序列視為一個隨機序列，用數(shù)學模型來近似描述該序列．此模型一旦被識別，就可以從時間序列的過去值以及現(xiàn)在值來預測未來值．

Kalman濾波器屬于時變線性系統(tǒng)的遞歸濾波器，通過遞歸算法獲取變量的最佳估計值，是將過去的測量估計誤差合并到新的測量誤差中來估計將來的誤差．Kalman濾波器預測算法由狀態(tài)方程(3)以及觀測方程(4)組成．

Xt+1=AXt+Wt

(3)

Yt+1=BXt+Vt

(4)

其中，Xt+1為狀態(tài)矢量；Yt+1為觀測矢量；A為狀態(tài)矩陣；B為觀測矩陣；Wt和Vt為對應的白噪聲矢量矩陣．

確定狀態(tài)方程和預測方程的過程十分復雜，本研究擬將ARIMA模型的數(shù)學表達式引入到Kalman濾波器的狀態(tài)方程和測量方程中，以期預測產(chǎn)油量狀態(tài)．

令Q1(t)=Q(t),Q2(t)=Q(t-1), …,Qp(t)=Q(t-p+1);e1(t)=e(t),e2(t)=e1(t-1), …,eq(t)=eq-1(t-1)， 并將產(chǎn)油量數(shù)據(jù)ARIMA模型引入到Kalman濾波器預測算法的狀態(tài)方程和測量方程中，則ARIMA模型可描述為

Q1(t+1)=φ1(t)Q1(t)+φ2(t)Q2(t) +…+

φp(t)Qp(t)+e1(t+1)-

θ1(t)e1(t)-θ2(t)e2(t)-…-

θq(t)eq(t)

(5)

由式(5)得出，Q2(t+1)=Q1(t),Q3(t+1)=Q2(t), …,Qp+1(t+1)=Qp(t);e2(t+1)=e1(t),e3(t+1)=e2(t), …,eq+1(t+1)=eq(t)． 因而可得ARIMA-Kalman濾波器預測算法的狀態(tài)方程為

(6)

根據(jù)式(2)至式(6)，可得ARIMA-Kalman濾波器預測算法的觀測方程為

Y(t+1)= [1, 0, …, 0]×

[Q1(t+1),Q2(t+1), …,

Qp(t+1)]T

(7)

由式(7)可見，ARIMA-Kalman濾波器預測算法基于歷史產(chǎn)油量數(shù)據(jù)構(gòu)建ARIMA模型，將ARIMA模型引入到Kalman濾波器中構(gòu)建狀態(tài)和測量方程，完成對噪聲的濾波過程，最后根據(jù)確定的最優(yōu)模型進行向后預測．圖1是該預測算法實現(xiàn)的流程圖．

圖1 ARIMA-Kalman濾波器預測算法流程圖Fig.1 Predictive algorithm flowchart of ARIMA-Kalman filter

2 試驗處理

2.1 平穩(wěn)性處理

載入目標油井1991年1月—2016年9月月產(chǎn)油量時間序列數(shù)據(jù)，并繪制產(chǎn)油量時間序列曲線，結(jié)果如圖2． 通過觀察該數(shù)據(jù)序列圖無法確定產(chǎn)油量時間序列是否穩(wěn)定，需要對產(chǎn)油量時間序列進行平穩(wěn)性檢驗，同時采用確定性時序分析方法提取產(chǎn)油量時間序列中所蘊涵的確定性信息．

圖2 月產(chǎn)油量時間序列圖Fig.2 Oil production time series curve

ADF檢驗(augmented Dickey-fuller test)[15]假設序列至少存在1個單位根，即非平穩(wěn)，對于一個平穩(wěn)的時序數(shù)據(jù)，需在給定的置信水平上(通常置信區(qū)間取95%)顯著，拒絕原假設．本研究為直觀地分析產(chǎn)油量時間序列，采用對數(shù)變換(以自然對數(shù)為底)減小數(shù)據(jù)的振動幅度，消除異方差問題．經(jīng)過對數(shù)變換的產(chǎn)油量時間序列曲線為圖3中原始曲線，ADF檢驗結(jié)果如表1．由表1可知，ADF=-0.799 304，高于1%顯著水平，P=0.819 326， 遠大于0.05，接受原假設，說明產(chǎn)油量序列是一個非平穩(wěn)序列．

圖3 產(chǎn)油量時間序列平滑曲線Fig.3 Oil production time series smooth curve

運用ARIMA模型進行時間序列分析需要將非平穩(wěn)時間序列轉(zhuǎn)化為平穩(wěn)時間序列．為獲得平穩(wěn)的產(chǎn)油量時間序列，本研究采用式(8)移動平均和式(9)指數(shù)平滑兩種確定性信息提取方法，其中n=12． 對對數(shù)變換的產(chǎn)油量序列進行移動平均和指數(shù)平滑結(jié)果如圖3．

表1 產(chǎn)油量時間序列ADF檢驗結(jié)果Table 1 Oil production time series ADF test results

(8)

Qt=αQt-1+α(1-α)Qt-2+

α(1-α)2Qt-3+…

(9)

其中，α為平滑系數(shù)，本研究取α=2/(n+1)．

然而，移動平均和指數(shù)平滑對確定性信息提取都不夠充分，BOX和JENKINS使用大量案例分析證明了差分處理是一種有效的確定性信息處理方法[16]，Cramer分解定理，如式(10)，則在理論上保證了適當階數(shù)的拆分可以充分提取確定性信息．

(10)

其中，β1,β2, …,βd為常數(shù)系數(shù)；at為一個零均值白噪聲序列；Ψ(B)為隨機性參數(shù)，B為延遲算子．

在Cramer分解定理保證下，式(11)的d階差分就可以將{Qt}中蘊含的確定性信息提取出來．

(11)

其中， c為常數(shù)． d=1時展開得出1階差分公式為

(12)

一階差分后通過ADF檢驗確定序列的平穩(wěn)性，結(jié)果如表2．由表2可見，ADF=-11.883 13，顯著低于1%顯著水平，P值遠小于0.05，拒絕原假設，差分處理后的序列是一個平穩(wěn)序列，即d=1．

表2 差分后產(chǎn)油量時間序列ADF檢驗結(jié)果Table 2 Differential oil production time series ADF test results

由圖3可見，產(chǎn)油量時間序列具有明顯的長期趨勢，同時產(chǎn)油量時間序列也存在周期性的可能．依據(jù)相加分解模型和相乘分解模型理論分解產(chǎn)油量時間序列，可獲得其長期趨勢、周期性趨勢以及殘差曲線(圖4)．由圖4可見，產(chǎn)油量周期成分區(qū)間波動幅度小，可忽視周期性因素的影響．ARIMA差分模型白噪聲檢測P=4.844 24×10-13， 顯著小于0.05，屬于一個白噪聲序列．因此，建模過程不需要考慮周期性組件以及重新定義(p,q)．

圖4 分解后曲線成分Fig.4 Decomposed curve components

2.2 模型識別

產(chǎn)油量時間序列進行平穩(wěn)性處理確定了差分項d， 隨后分析ACF函數(shù)以及PACF函數(shù)，并確定參數(shù)(p,q)[17]． ACF函數(shù)是指有序的隨機變量序列與其自身相比較反映了同一序列在不同時序的取值的相關性．PACF函數(shù)是指在剔除中間k-1個隨機變量的干擾后，Qt-k對Qt影響的相關度量．統(tǒng)計分析對數(shù)變換后的產(chǎn)油量時間序列，獲得ACF函數(shù)(圖5)以及PACF函數(shù)(圖6)．根據(jù)ACF圖中曲線第1次穿過置信區(qū)間時獲取q=1， 又根據(jù)PACF圖中曲線第1次穿過置信區(qū)間時獲取p=1.

圖5 ACF函數(shù)圖Fig.5 ACF function diagram

圖6 PACF函數(shù)圖Fig.6 PACF function diagram

依據(jù)AIC準則確定模型(1, 1)是否為最優(yōu)的滯后因子階數(shù)，鼓勵數(shù)據(jù)擬合的優(yōu)良性，但盡量避免出現(xiàn)過度擬合的情況．所以，優(yōu)先考慮的模型應是AIC值最小的那一個．AIC準則的方法是尋找可以最好地解釋數(shù)據(jù)，但包含最少自由參數(shù)的模型．表3給出不同(p,d,q)情況下的AIC統(tǒng)計值，從表3可見，(1, 1, 2)為最優(yōu)模型．

表3 AIC統(tǒng)計結(jié)果Table 3 AIC statistics results

確定產(chǎn)油量時間序列ARIMA模型結(jié)構(gòu)識別為(1, 1, 2)，檢驗ARIMA(1, 1, 2)模型獲取自回歸以及移動平均參數(shù)，則該模型描述為

Q(t+1)= -0.009 6Q(t)+0.776 7Q(t-1)+

e(t+1)-1.352 5e(t)+

0.384 4e(t-1)

(13)

式(13)的ARIMA模型引入Kalman濾波器預測算法后，得出產(chǎn)油量時間序列ARIMA-Kalman濾波器模型為

(14)

2.3 數(shù)據(jù)擬合

根據(jù)產(chǎn)油量時間序列ARIMA-Kalman濾波器模型還原數(shù)據(jù)擬合效果，依據(jù)標準化的過程，將預測的產(chǎn)油量逆向還原，并根據(jù)歷史上的實際月產(chǎn)油量，與實際值進行對比檢驗．從圖7模型擬合結(jié)果可見，實際值與擬合值的吻合程度較高．

圖7 產(chǎn)油量時間序列擬合結(jié)果Fig.7 Oil production time series fitting results

2.4 預 測

采用確定的ARIMA-Kalman濾波器模型預測2016年10月—2017年9月的單井實驗區(qū)的產(chǎn)油量，并與該時間段的實際產(chǎn)油量進行對比(構(gòu)造ARIMA-Kalman濾波器模型未載入2016年10月—2017年9月實際產(chǎn)油量數(shù)據(jù))，結(jié)果如表4．由表4可見，此時間段內(nèi)單井實驗區(qū)的實際產(chǎn)油量為650.00 t，采用ARIMA-Kalman濾波器的預測產(chǎn)量為658.83 t，殘差為8.83 t，相對誤差為1.36%．

3 試驗結(jié)果分析

產(chǎn)油量時間序列的影響因素十分復雜，一個微小的相關因素也可能引起產(chǎn)油量曲線發(fā)生急劇變化，同時各個因素之間存在一定的時序影響關系．同時，各種原因?qū)е碌慕Y(jié)果需要一定的時間過程方能顯現(xiàn)，因此預測的產(chǎn)油量時間序列具有非同步性和滯后性．

表4 2016年10月—2017年9月產(chǎn)油量實際值與預測值Table 4 Actual and forecasted oil production from October 2016 to September 2017

觀察產(chǎn)油量時間序列，產(chǎn)油量曲線整體呈現(xiàn)長期遞減趨勢，主要是因為經(jīng)過長期開采該區(qū)塊儲層剩余油儲量逐漸減少．同時，為了減少層間干擾，工程人員采取分層注水技術提高注水波及系數(shù)，隨后逐步采取了水力壓裂、酸化等增產(chǎn)措施，因而2016年11月—2017年2月的產(chǎn)油量又有了顯著回升．對比預測產(chǎn)量與實際產(chǎn)量，實際產(chǎn)量在2017年1月到達峰值，預測產(chǎn)量在2017年2月達到峰值，具有微小的滯后性，然而， EDIGER的研究模型[7]非同步性與滯后性顯著．總體來說，產(chǎn)油量整體呈現(xiàn)長期遞減趨勢，局部呈現(xiàn)時升時降趨勢，同時不能忽略個別月份油井產(chǎn)油時間較少等因素的影響．

4 結(jié) 論

1)建立了具有高效的時序影響因素的分析能力ARIMA-Kalman濾波器模型，該模型能夠排除非同步性以及滯后性的影響，使識別出的產(chǎn)油量時間序列模型具有精準的擬合結(jié)果和預測能力．

2)運用ARIMA-Kalman濾波器模型進行區(qū)塊產(chǎn)油量效果預測，2016年10月—2017年9月實際產(chǎn)量650.00 t，預測產(chǎn)量658.83 t，殘差8.83 t，相對誤差為1.36%，預測結(jié)果精度較高．

3)ARIMA-Kalman濾波器模型對油田產(chǎn)量預測是一種有效的嘗試，該方法可為國內(nèi)外油田產(chǎn)量及油田開發(fā)過程中含水率、產(chǎn)液量等宏觀預測提供一種思路，為油田開發(fā)提供決策與理論依據(jù)．