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

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

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

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

1 ARIMA-Kalman濾波器模型原理

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

(1)

其中，φ(B)為AR多項(xiàng)式；θ(B)為MA多項(xiàng)式；t為時(shí)間序列下標(biāo)，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)的時(shí)間序列用于建立ARIMA模型并確定其參數(shù)．參數(shù)(p,q)通過自相關(guān)函數(shù)(auto correlation function, ACF)和偏自相關(guān)函數(shù)(partial autocorrelation function, PACF)預(yù)先確定，根據(jù)AIC信息準(zhǔn)則(Akaike information criterion)[13]最終確定．在時(shí)間序列分析中，Yule-Walker方程在模式識(shí)別和參數(shù)估計(jì)中起著重要作用[14]．ARIMA模型的殘差可以作為模型估計(jì)準(zhǔn)則，模型估計(jì)過程采用最大似然估計(jì)．通過求解殘差序列的ACF，可以判斷殘差序列是否為白噪聲．如果殘差序列不是白噪聲，則重新定義(p,q)．

因此，t+1時(shí)刻預(yù)測(cè)的產(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)油量時(shí)間序列；e(t+1),e(t),e(t-1), …,e(t-q+1)為殘差時(shí)間序列．

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

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

Xt+1=AXt+Wt

(3)

Yt+1=BXt+Vt

(4)

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

確定狀態(tài)方程和預(yù)測(cè)方程的過程十分復(fù)雜，本研究擬將ARIMA模型的數(shù)學(xué)表達(dá)式引入到Kalman濾波器的狀態(tài)方程和測(cè)量方程中，以期預(yù)測(cè)產(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濾波器預(yù)測(cè)算法的狀態(tài)方程和測(cè)量方程中，則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濾波器預(yù)測(cè)算法的狀態(tài)方程為

(6)

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

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

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

Qp(t+1)]T

(7)

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

圖1 ARIMA-Kalman濾波器預(yù)測(cè)算法流程圖Fig.1 Predictive algorithm flowchart of ARIMA-Kalman filter

2 試驗(yàn)處理

2.1 平穩(wěn)性處理

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

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

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

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

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

表1 產(chǎn)油量時(shí)間序列ADF檢驗(yàn)結(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)．

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

(10)

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

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

(11)

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

(12)

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

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

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

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

2.2 模型識(shí)別

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

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

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

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

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

確定產(chǎn)油量時(shí)間序列ARIMA模型結(jié)構(gòu)識(shí)別為(1, 1, 2)，檢驗(yàn)ARIMA(1, 1, 2)模型獲取自回歸以及移動(dòng)平均參數(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濾波器預(yù)測(cè)算法后，得出產(chǎn)油量時(shí)間序列ARIMA-Kalman濾波器模型為

(14)

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

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

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

2.4 預(yù) 測(cè)

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

3 試驗(yàn)結(jié)果分析

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

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

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

4 結(jié) 論

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

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

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