林文輝，杜彥煒，趙 鵬

(1.西安工業(yè)大學(xué) 機(jī)電工程學(xué)院 ,西安 710021;2.西安工業(yè)大學(xué) 新生院,西安 710021)

能源是現(xiàn)代經(jīng)濟(jì)社會發(fā)展的基礎(chǔ),過去100年世界能源結(jié)構(gòu)發(fā)生了巨大變化,天然氣在能源結(jié)構(gòu)中所占比例正穩(wěn)步上升[1]。相比于即將到來的原油峰值,距中國天然氣產(chǎn)量(下文簡稱天然氣產(chǎn)量)達(dá)峰仍有數(shù)十年[2],結(jié)合未來天然氣水合物勘探開發(fā)的突破,天然氣產(chǎn)量增長潛力較大[3]。一方面我國提出2030年碳達(dá)峰的目標(biāo),另一方面碳達(dá)峰不能以損害能源安全為代價[4]。天然氣作為化石能源中的清潔能源,在相當(dāng)長一段時間內(nèi)將是我國能源戰(zhàn)略和碳減排實施的重要保障[5]。而研判能源發(fā)展形勢,謀劃天然氣戰(zhàn)略布局又對天然氣的產(chǎn)量預(yù)測提出迫切要求。

天然氣產(chǎn)量預(yù)測方法主要包括生命旋回模型、儲采比控制法、產(chǎn)量構(gòu)成法等。生命旋回模型中,廣義翁氏模型和Hubbert模型最為常見,其將產(chǎn)量的周期性變化視作一個生命旋回,通過構(gòu)建循環(huán)模型進(jìn)行預(yù)測。文獻(xiàn)[6]提出一種指數(shù)修正翁氏模型用于四川盆地油氣預(yù)測,文獻(xiàn)[7]采用多循環(huán)Hubbert模型預(yù)測了中國未來天然氣產(chǎn)量。儲采比即剩余可采儲量與當(dāng)年實際產(chǎn)量的比值,反映了油氣資源的保證程度。文獻(xiàn)[8]運(yùn)用油藏數(shù)值模擬等方法,得到儲采比與累計產(chǎn)量、采出程度等之間的關(guān)系。產(chǎn)量構(gòu)成法依據(jù)不同類型的氣藏在不同階段的生產(chǎn)規(guī)律、剩余可采儲量等指標(biāo),計算每個單位的產(chǎn)量潛力,按照時間累加得到總產(chǎn)量的規(guī)模。文獻(xiàn)[9]使用產(chǎn)量構(gòu)成法預(yù)測分別預(yù)測未來中國常規(guī)天然氣產(chǎn)量、非常規(guī)天然氣產(chǎn)量以及溶解氣產(chǎn)量。

近年來國內(nèi)外學(xué)者開始將自回歸滑動平均求和模型(Autoregressive Integrated Moving Average Model,ARIMA)等計量經(jīng)濟(jì)學(xué)模型及機(jī)器學(xué)習(xí)方法用于天然氣產(chǎn)量預(yù)測。文獻(xiàn)[10]將自回歸滑動平均求和模型與長短期記憶神經(jīng)網(wǎng)絡(luò)(Long Short-Term Memory,LSTM)結(jié)合起來建立天然氣產(chǎn)量的集成預(yù)測模型,文獻(xiàn)[11]提出一種季節(jié)性灰色模型用于對天然氣產(chǎn)量進(jìn)行預(yù)測,使用粒子群算法優(yōu)化模型參數(shù),文獻(xiàn)[12]基于多目標(biāo)隨機(jī)森林模型對頁巖氣產(chǎn)量進(jìn)行預(yù)測,文獻(xiàn)[13]基于Logistic模型,將原模型中的“環(huán)境容量”和“種群內(nèi)稟增長率”兩個常量參數(shù)函數(shù)化以改進(jìn)模型,使用改進(jìn)Logistic模型對中國天然氣產(chǎn)量進(jìn)行了預(yù)測。

以往研究未充分考慮天然氣產(chǎn)量序列的時序趨勢、不同預(yù)測步長下的模型變化以及對過去預(yù)測誤差信息的利用?；诖?本文先選取能夠較好捕捉中短期時序特征的ARIMA模型,對序列異常波動具有強(qiáng)抗擾動能力、長期預(yù)測效果穩(wěn)定的線性回歸(Liner Regression, LR)模型,和針對信息不完全的情形有較好表現(xiàn)的GM(1,1)模型進(jìn)行組合,以進(jìn)行模型之間的互補(bǔ)。再加入自適應(yīng)學(xué)習(xí)因子與組合學(xué)習(xí)因子,分別對單個模型和組合模型加以改進(jìn),實現(xiàn)對過去預(yù)測誤差信息的充分利用,使得最終的集成預(yù)測模型具有更高的預(yù)測精度。

1 相關(guān)方法

1.1 ARIMA模型

ARIMA模型是由Box和Jenkins提出的一類時間序列分析方法[14],其通過差分對序列平穩(wěn)化,可用于非平穩(wěn)時間序列的預(yù)測,在模型建立充分使用觀測到時間序列的過去項、當(dāng)前項和誤差擾動項,從而使模型具有較高的預(yù)測精度。

用{Yt}表示觀測到的時間序列,εt表示未觀測到的白噪聲序列,即一組服從獨(dú)立同分布、均值為零的隨機(jī)變量。φ表示白噪聲變量的權(quán)重,則{Yt}可表示為一系列白噪聲變量的加權(quán)線性組合[15]:

Yt=εt+φ1et-1+φ2et-2,

(1)

當(dāng)φ不為零時,用1,-θ1,-θ2,-θ3,…,-θq作為et,et-1,et-2,…,et-q的權(quán)重加平均得到{Yt},即得到q階滑動平均過程(MA),記為MA(q):

Yt=et-θ1et-1-θ2et-2-…-θqet-1,

(2)

當(dāng)用{Yt}自身做回歸變量時,用ψ1,ψ2,ψ3,…,ψp作為Yt-1,Yt-2,…,Yt-p的權(quán)重,{Yt}對應(yīng)的白噪聲變量et表示無法用自身解釋的新信息項,得到p階滑動平均過程,記為AR(p):

Yt=φ1Yt-1+φ2Yt-2+…+φpYt-p+et,

(3)

若序列中部分是自回歸,部分是滑動平均,可得到一個混合過程,即自回歸滑動平均過程,記為ARMA(p,q):

Yt=et-θ1et-1-θ2et-2-…-θqet-q+φ1Yt-1+
φ2Yt-2+…+φpYt-p+et,

(4)

ARMA要求時間序列具有平穩(wěn)性,在序列非平穩(wěn)的情形下,需通過差分計算使序列平穩(wěn)化,表示向后差分算子,Yt表示Yt的一階差分,即:

dYt=θ0+θq(B)εt,

(5)

如果{Yt}的d次差分Wt服從ARMA(p,q)模型,則稱{Yt}為ARIMA(p,d,q)過程。參數(shù)p,d,q用來分析時間序列,即描述自回歸階數(shù)p、差分次數(shù)d和移動平均階數(shù)q,表示為[15]:

dYt=θ0+θq(B)εt,

(6)

式中:θq為移動平均算子;θ0為常數(shù)項;B為后移算子,BYt=Yt-1,定義θ0=u(1-φ1-φ2-…-φp),其中u為平均數(shù)。

由于傳統(tǒng)的通過觀察模型自相關(guān)函數(shù)的差分?jǐn)?shù)據(jù)圖(auto correlation function,ACF)和偏自相關(guān)函數(shù)的偏自相關(guān)圖(patial auto correlation function,PACF)進(jìn)行估計模型參數(shù)的方法在多次預(yù)測時耗時長、帶有一定的主觀性,本文使用統(tǒng)計學(xué)中的貝葉斯信息準(zhǔn)則(Bayesian Information Criterion,BIC)指標(biāo),通過對指定參數(shù)范圍內(nèi)的模型進(jìn)行比較,實現(xiàn)對模型的自動定階。

1.2 灰色GM(1,1)

灰色指的是信息不完全,由于天然氣產(chǎn)量變化與許多因素有關(guān),而僅使用歷史天然氣產(chǎn)量數(shù)據(jù)構(gòu)建的一個預(yù)測系統(tǒng)可視為一個信息不完全的系統(tǒng),即灰色系統(tǒng)?；疑獹M(1,1)是通過對原始數(shù)據(jù)進(jìn)行累加生成來進(jìn)行建模,再對預(yù)測結(jié)果進(jìn)行累減還原得到預(yù)測結(jié)果[16],其適用于對小樣本、貧信息的不確定系統(tǒng)進(jìn)行預(yù)測[17]。

設(shè)原始序列Y(0)=(y(0)(1),y(0)(2),…,y(0)(n)),其中y(0)(t)≥0,t=1,2,…,n;Y(1)為Y(0)的一階累加序列:

Y(1)=(y(1)(1),y(1)(2),…,y(1)(n),

(7)

y(0)(t)+ax(1)(t)=y,

(8)

為GM(1,1)模型[17]。式中,a稱為發(fā)展灰數(shù),反映Y(1)及Y(0)的發(fā)展趨勢,u稱為內(nèi)生控制灰數(shù),體現(xiàn)數(shù)據(jù)間的變化關(guān)系。利用最小二乘法求解得到a、u,繼續(xù)求解微分方程式

(9)

作累減還原后,便可得到預(yù)測值

(10)

1.3 LR模型

回歸是通過一定的數(shù)學(xué)表達(dá)是將變量之間的關(guān)系描述出來的一種方法,當(dāng)所要描述的因變量y和自變量x之間為線性關(guān)系時,則稱為線性回歸(Liner Regression,LR)。自變量個數(shù)為1時,稱為一元線性回歸。以下是一元線性回歸模型的表示,ε表示誤差項[18]:

y=β0+β1x+ε,

(11)

(12)

高斯提出用最小化離差平方和來估計參數(shù)β0和β1,這一估計方法稱為最小二乘法,根據(jù)最小二乘法,令

(13)

解上式得

(14)

2 GM(1,1)-ARIMA-LR模型的構(gòu)建

2.1 自適應(yīng)學(xué)習(xí)因子

單個模型在進(jìn)行預(yù)測的時候,可能由于時滯以及模型本身特性的影響,導(dǎo)致預(yù)測結(jié)果總是偏大或偏小。此時盡管模型的預(yù)測誤差平均值較大,但預(yù)測誤差的方差較小,即預(yù)測結(jié)果較為穩(wěn)定。若是直接將這樣的幾個模型進(jìn)行集成,由于單個模型的預(yù)測誤差較大,集成模型的效果將受到影響。文中提出一種自適應(yīng)學(xué)習(xí)因子F1,如式(15)所示。先通過模型對歷史數(shù)據(jù)進(jìn)行預(yù)測后得到一個預(yù)測誤差百分比列表E,然后將當(dāng)前時間減去預(yù)測步長得到上個周期對應(yīng)時間t-k。

將上個周期的預(yù)測誤差百分比進(jìn)行變換后與模型本身相乘,從而得到修正后的預(yù)測結(jié)果。

(15)

2.2 組合學(xué)習(xí)因子

通過加權(quán)平均對模型進(jìn)行集成,關(guān)鍵在于模型的選取以及權(quán)重系數(shù)的確定。權(quán)重系數(shù)大多采用兩種方式確定,一是平均策略,對每個模型賦予一樣的權(quán)重,二是通過不斷嘗試選取較佳的權(quán)重。但這樣的做法以下問題:一是難以取到最優(yōu)權(quán)重,二是權(quán)重的選取具有較強(qiáng)的主觀和偶然性,三是權(quán)重選取的方式是靜態(tài)的,難以適應(yīng)預(yù)測步長的變化。在一步預(yù)測誤差較小的權(quán)重系數(shù),可能預(yù)測步長改變時誤差放大。為提升集成模型預(yù)測精度,本文根據(jù)預(yù)測步長的不同設(shè)置可變的權(quán)重系數(shù),并且建立一個目標(biāo)規(guī)劃模型,將求解得到的權(quán)重系數(shù)命名為組合學(xué)習(xí)因子F2。

本文通過改變參與預(yù)測的數(shù)據(jù)集來實現(xiàn)預(yù)測步長對權(quán)重系數(shù)的調(diào)節(jié)。若對天然氣產(chǎn)量進(jìn)行單步預(yù)測,則選取起始年份至預(yù)測年份前一年的數(shù)據(jù)。若對天然氣產(chǎn)量進(jìn)行多步預(yù)測,則先用預(yù)測年份減去預(yù)測步數(shù)得到一個終止年份,再選取起始年份到終止年份的數(shù)據(jù)作為數(shù)據(jù)集。例如,要進(jìn)行一步預(yù)測得到2023年天然氣產(chǎn)量的預(yù)測值,則選取2002～2022年數(shù)據(jù)作為輸入數(shù)據(jù)集。要進(jìn)行二步預(yù)測得到2023年天然氣產(chǎn)量的預(yù)測值,則選取2002～2021年數(shù)據(jù)作為輸入數(shù)據(jù)集。

目標(biāo)規(guī)劃模型以近5年天然氣產(chǎn)量預(yù)測平均誤差百分比的絕對值最小為目標(biāo)為

(16)

式中:Ri為天然氣產(chǎn)量數(shù)據(jù)的實際值;Pi為集成模型的預(yù)測值。

由于目標(biāo)函數(shù)中帶有絕對值,為便于模型求解,引入輔助變量Q,加上第二、三約束條件,可去掉目標(biāo)函數(shù)中的絕對值,進(jìn)而目標(biāo)函數(shù)可表示為

(17)

為了更清晰地表示,可將上式等價為

(18)

式中:y為時間片段的長度,即取y年平均誤差百分比作為目標(biāo)值。y取1、2時容易受異常值的影響,所以分別將y取3、4或5,最后取對應(yīng)目標(biāo)值最小的解作為最優(yōu)解。所以目標(biāo)函數(shù)可進(jìn)一步表示為

(19)

接下來是約束部分。權(quán)重系數(shù)優(yōu)化模型的第一個約束,是Pt的求解:

(20)

式中:a為ARIMA模型的權(quán)重系數(shù);b為灰色GM(1,1)的權(quán)重系數(shù);c為LR模型的權(quán)重系數(shù);A為ARIMA模型的預(yù)測值;G為灰色GM(1,1)模型的預(yù)測值;L為LR模型的預(yù)測值。

第二、 三個約束,是引入輔助變量Q后新增的約束:

(Ri-Pi)/Ri≤Qj,

(21)

(Ri-Pi)/Ri≤Qj,

(22)

第四個約束即權(quán)重系數(shù)的非負(fù)約束為

a,b,c>0,

(23)

綜上所述,權(quán)重系數(shù)的優(yōu)化模型表示為

(24)

組合學(xué)習(xí)因子即上述權(quán)重系數(shù)優(yōu)化模型求解得到的系數(shù)a、b、c,應(yīng)用組合學(xué)習(xí)因子對所要預(yù)測年份的各單模型預(yù)測值進(jìn)行加權(quán)求和,即可得到集成模型的預(yù)測值Pt:

F2=[a,b,c]

Pt=[a,b,c]×[At,Gt,Lt]T,

(25)

2.3 GM(1,1)-ARIMA-LR集成預(yù)測模型

GM(1,1)-ARIMA-LR集成模型預(yù)測框架如圖1所示。首先獲取數(shù)據(jù)集,根據(jù)步長選取數(shù)據(jù),分別使用ARIMA、GM(1,1)、LR模型進(jìn)行預(yù)測得到初步結(jié)果,再加入自適應(yīng)學(xué)習(xí)因子F1對結(jié)果進(jìn)行調(diào)整得到改進(jìn)后的結(jié)果。將得到的單個模型預(yù)測結(jié)果輸入已建立的權(quán)重系數(shù)優(yōu)化模型,求解后得到組合學(xué)習(xí)因子F2。最后對模型進(jìn)行融合,輸出GM(1,1)-ARIMA-LR模型預(yù)測結(jié)果。

圖1 GM(1,1)-ARIMA-LR預(yù)測框架

3 實例預(yù)測結(jié)果與分析

3.1 數(shù)據(jù)獲取及處理

從國家統(tǒng)計局網(wǎng)站獲取1973-2022年中國天然氣產(chǎn)量(億立方米)年度數(shù)據(jù)。首先進(jìn)行ARIMA、灰色GM(1,1)和LR單模型的預(yù)測。天然氣產(chǎn)量數(shù)據(jù)的趨勢圖如圖2所示,整體上天然氣產(chǎn)量變化呈穩(wěn)定上升趨勢,曲線斜率接近。但是在一些年份,如1980年、1993年、2009年、2017年等,天然氣產(chǎn)量的曲線斜率發(fā)生了較大的變化。

圖2 1973～2022年天然氣產(chǎn)量

預(yù)測算法程序采用Anaconda3開發(fā),軟件環(huán)境為Python3.9.7,硬件運(yùn)行環(huán)境為Windows10、64位操作系統(tǒng),處理器為Intel(R)Core(TM)i5-7300HQCPU@2.50 GHz,內(nèi)存8 GB。

3.2 預(yù)測建模

將數(shù)據(jù)集分別代入ARIMA、GM(1,1)、LR模型進(jìn)行預(yù)測,接著根據(jù)自適應(yīng)學(xué)習(xí)因子F1對模型進(jìn)行調(diào)整。圖3為一步至八步預(yù)測的ARIMA模型與改進(jìn)ARIMA模型預(yù)測誤差的比較,觀察圖3可見,在一到四步預(yù)測,使用自適應(yīng)學(xué)習(xí)因子對模型進(jìn)行改進(jìn)顯著降低了模型預(yù)測誤差,而在五到八步預(yù)測,該改進(jìn)反而使預(yù)測誤差增大。因而在集成模型的構(gòu)建中,加入變量SA,決定是否加入自適應(yīng)學(xué)習(xí)因子F1。

圖3 ARIMA模型改進(jìn)前后誤差對比

當(dāng)預(yù)測步長k≤4時,取SA=1,表示加入F1,對ARIMA模型進(jìn)行改進(jìn),并用改進(jìn)ARIMA模型進(jìn)行預(yù)測;當(dāng)預(yù)測步長k>4時,取SA=0,表示不加入F1,直接使用原ARIMA模型得到預(yù)測結(jié)果。

同理,觀察圖4和圖5可見,當(dāng)預(yù)測步長為一到四時,改進(jìn)后的GM(1,1)預(yù)測誤差大幅降低,預(yù)測步長為五到八時,改進(jìn)后GM(1,1)預(yù)測誤差反而增大。在一到八步預(yù)測,改進(jìn)后LR模型預(yù)測誤差均顯著降低,因而可一直應(yīng)用F1對LR模型進(jìn)行改進(jìn),加入變量SG,決定是否加入自適應(yīng)學(xué)習(xí)因子F1對GM(1,1)模型進(jìn)行改進(jìn)。

圖4 GM(1,1)模型改進(jìn)前后誤差對比

圖5 LR模型改進(jìn)前后誤差對比

當(dāng)預(yù)測步長k≤4時,取SG=1,表示加入F1,對GM(1,1)模型進(jìn)行改進(jìn),并用改進(jìn)GM(1,1)模型進(jìn)行預(yù)測;當(dāng)預(yù)測步長k>4時,取SG=0,表示不加入F1,直接使用原GM(1,1)模型得到預(yù)測結(jié)果。

繼續(xù)將上述得到的初步預(yù)測結(jié)果代入權(quán)重系數(shù)優(yōu)化模型,求解得到對應(yīng)的權(quán)重系數(shù),即組合學(xué)習(xí)因子F2。通過對比發(fā)現(xiàn),在一、二步預(yù)測中,改進(jìn)GM(1,1)模型預(yù)測誤差最小,分別為1.187%,2.695%;在三到五步預(yù)測中,改進(jìn)GM(1,1)-ARIMA-LR模型表現(xiàn)較優(yōu),預(yù)測誤差分別為4.243%,2.698%,3.129%;在六到八步預(yù)測中,改進(jìn)LR模型效果較好,預(yù)測誤差分別為6.550%,7.423%,9.880%。

通過設(shè)置權(quán)重系數(shù)的方式,可令一到二步的短期預(yù)測中系數(shù)a,b,c的值為0,1,0,六到八步的長期預(yù)測中,系數(shù)a,b,c的值為0,0,1,三到五步a,b,c的值則繼續(xù)使用權(quán)重系數(shù)優(yōu)化模型得到的結(jié)果,從而進(jìn)一步調(diào)整得到最終的改進(jìn)ARIMA-GM(1,1)-LR模型。以2022年為例,權(quán)重系數(shù)的變化如圖6所示。

圖6 權(quán)重系數(shù)變化

可見,隨著預(yù)測步長的增加,GM(1,1)模型的權(quán)重占比越來越小,LR模型的權(quán)重占比越來越大。說明在進(jìn)行不同步長的預(yù)測過程中,ARIMA,GM(1,1)、LR模型實現(xiàn)了互相補(bǔ)充,從而使得改進(jìn)ARIMA-GM(1,1)-LR模型具有更高的預(yù)測精度。

3.3 預(yù)測效果評價及對比分析

以一步預(yù)測為例,表1為一步預(yù)測的結(jié)果對比,可展示為如圖7所示。觀察圖7可見,在一步預(yù)測中,改進(jìn)GM(1,1)的預(yù)測誤差最小,因而選擇改進(jìn)GM(1,1)模型進(jìn)行一步預(yù)測。

表1 一步預(yù)測結(jié)果對比

圖7 一步預(yù)測結(jié)果對比

圖8為多步預(yù)測的模型誤差對比?？梢园l(fā)現(xiàn),隨著預(yù)測步長的增長,模型的預(yù)測誤差也在不斷增加。集成模型ARIMA-GM(1,1)-LR的預(yù)測效果最好,改進(jìn)LR模型預(yù)測誤差較為穩(wěn)定,隨預(yù)測步長增加變化小。改進(jìn)GM(1,1)模型在一到二步的短期預(yù)測中表現(xiàn)優(yōu)異,之后預(yù)測誤差驟增。

圖8 各步長下預(yù)測誤差對比

3.4 2023-2030年天然氣產(chǎn)量預(yù)測

根據(jù)求解得到的自適應(yīng)學(xué)習(xí)因子F1和組合學(xué)習(xí)因子F2,對ARIMA、GM(1,1)、LR進(jìn)行加權(quán)求和,得到2023-2030年中國天然氣產(chǎn)量預(yù)測結(jié)果,見表2。2023年的天然氣產(chǎn)量為2022年天然氣產(chǎn)量的一步預(yù)測值,2024年的天然氣產(chǎn)量為2023年天然氣產(chǎn)量的二步預(yù)測值,以此類推。

表2 2023-2030年天然氣產(chǎn)量預(yù)測結(jié)果

4 結(jié) 論

1) 文中從對預(yù)測誤差信息的利用角度出發(fā),設(shè)計了自適應(yīng)學(xué)習(xí)因子和組合學(xué)習(xí)因子,對ARIMA、GM(1,1)、LR模型進(jìn)行改進(jìn),使得集成后的模型誤差大幅降低,在天然氣產(chǎn)量的短期、中期與長期預(yù)測均具有優(yōu)異表現(xiàn)。

2) 文中主要聚焦于小樣本情形下提高天然氣產(chǎn)量的預(yù)測精度。在樣本容量大、數(shù)據(jù)維度高的情形下,將影響天然氣產(chǎn)量變化的眾多因素考慮在內(nèi),應(yīng)用深度學(xué)習(xí)等方法構(gòu)建多變量預(yù)測模型或能增加預(yù)測精度,這是本研究之后改進(jìn)的方向。