基于GM-ARIMA模型的遙測數(shù)據(jù)中長期預(yù)測

李志強(qiáng),張香燕,李鴻飛,田華東

(北京空間飛行器總體設(shè)計(jì)部,北京 100094)

1 引言

衛(wèi)星遙測數(shù)據(jù)是地面監(jiān)視和評估衛(wèi)星在軌運(yùn)行狀態(tài)的主要依據(jù),它包含了各器件在軌運(yùn)行的豐富信息。如果衛(wèi)星器件在軌運(yùn)行期間發(fā)生故障,相應(yīng)的遙測數(shù)據(jù)的趨勢也會(huì)發(fā)生改變。開展遙測數(shù)據(jù)長期和短期預(yù)測研究,對遙測數(shù)據(jù)的變化規(guī)律進(jìn)行建模,預(yù)測未來時(shí)段的遙測數(shù)據(jù)變化趨勢,就可以提前預(yù)警衛(wèi)星潛在的故障,增強(qiáng)診斷預(yù)警系統(tǒng)的故障早期發(fā)現(xiàn)能力,為地面人員進(jìn)行提前干預(yù)和處置故障提供決策依據(jù),對保障衛(wèi)星長期安全穩(wěn)定運(yùn)行、開展衛(wèi)星長期性能研究等方面具有重要意義[1,2]。

衛(wèi)星面臨復(fù)雜的空間環(huán)境,影響遙測數(shù)據(jù)變化的因素具有不確定性,造成遙測數(shù)據(jù)表現(xiàn)出復(fù)雜的波動(dòng)規(guī)律,對遙測數(shù)據(jù)的預(yù)測面臨種種困難。一般地,常用的遙測數(shù)據(jù)預(yù)測方法有如曲線擬合法、小波分析法、灰色系統(tǒng)預(yù)測法、神經(jīng)網(wǎng)絡(luò)等方法[3-5],這些方法多適用于單調(diào)趨勢的短期預(yù)報(bào),未能充分考慮數(shù)據(jù)的季節(jié)性和隨機(jī)性變化特點(diǎn),對于中長期預(yù)報(bào)往往因?yàn)榫炔粔蚧蛘哳A(yù)測發(fā)散而存在缺陷。借助于時(shí)間序列分析原理,遙測數(shù)據(jù)可以分解成相應(yīng)時(shí)間序列的趨勢成分、波動(dòng)成分,對不同成分的分別預(yù)測可能成為提高預(yù)測精度的關(guān)鍵。HP濾波(Hodrick-Prescott Filter)作為一種常用的時(shí)間序列數(shù)據(jù)分析工具[6,7],可以將序列分解成趨勢項(xiàng)和波動(dòng)項(xiàng),已在能源、經(jīng)濟(jì)周期預(yù)測等多個(gè)領(lǐng)域廣泛應(yīng)用,取得了較好的應(yīng)用效果。基于此,本文嘗試將HP濾波應(yīng)用在衛(wèi)星遙測數(shù)據(jù)中長期預(yù)測中。

本文提出了一種基于HP濾波進(jìn)行時(shí)間序列分解的遙測數(shù)據(jù)預(yù)測方法,即應(yīng)用HP濾波法將數(shù)據(jù)分解成趨勢項(xiàng)和波動(dòng)項(xiàng),分開建模、組合預(yù)測,通過疊加二者得到數(shù)據(jù)的最終預(yù)測值。該方法有別于以往把數(shù)據(jù)作為整體進(jìn)行一次預(yù)測的方法,具有如下技術(shù)改進(jìn)優(yōu)勢:1)創(chuàng)新性地引入HP濾波分解出趨勢項(xiàng)和波動(dòng)項(xiàng),弱化了趨勢、波動(dòng)間的相互影響,分別進(jìn)行趨勢和波動(dòng)的預(yù)測;2)趨勢項(xiàng)預(yù)測采用灰色GM(1,1)模型,對單調(diào)趨勢性預(yù)測效果良好,規(guī)避了灰色預(yù)測對波動(dòng)性預(yù)測的不足;3)波動(dòng)項(xiàng)采用季節(jié)型ARIMA模型,充分考慮了季節(jié)性和隨機(jī)性特點(diǎn),對波動(dòng)性進(jìn)行較好的刻畫和預(yù)測。因此,該方法綜合了各建模方式的優(yōu)點(diǎn),具有較高的預(yù)測精度,可用于數(shù)據(jù)中長期預(yù)測,為遙測參數(shù)的分析和預(yù)測、故障診斷和預(yù)警提供了一種全新的技術(shù)方法。

2 理論基礎(chǔ)

2.1 HP濾波

HP濾波是由Hodrick和Prescot提出的一種將時(shí)間序列進(jìn)行分解的方法[6]。它的基本原理是將原序列Y={y1,y2,…,yn}分為兩組:趨勢序列X={x1,x2,…,xn}和波動(dòng)序列C={c1,c2,…,cn},其與原序列的關(guān)系如式(1)[6,7]

yt=xt+ct,t=1,2,…n

(1)

HP濾波的目的是將平滑的時(shí)間序列即趨勢成分從不平滑的時(shí)間序列中分解出來,分離過程必須滿足損失函數(shù)最小原則

(2)

2.2 灰色預(yù)測GM(1,1)模型

GM(1,1)模型是鄧聚龍教授在1982年提出的一種應(yīng)用最廣泛的灰色預(yù)測模型[8]。它是一種研究少數(shù)據(jù)、貧信息、不確定性問題的方法,基本思路是首先利用累加的方法使數(shù)據(jù)具備指數(shù)規(guī)律,然后建立一階微分方程并對其求解,并將所求結(jié)果再累減還原,即為灰色預(yù)測值,從而對未來進(jìn)行預(yù)測。GM(1,1)模型建模步驟如下[9,10]

設(shè)原始數(shù)據(jù)序列為X(0),X(0)=(x(0)(1),x(0)(2),…,x(0)(n)),對X(0)作一次累加操作(1-AGO),得到序列X(1)=(x(1)(1),x(1)(2),…,x(1)(n)),其中

(3)

對X(1)作緊鄰均值生成,得到均值生成序列Z(1)={z(1)(2),z(1)(3),…,z(1)(n)},其中

(4)

則GM(1,1)灰色微分方程為

x(0)(k)+az(1)(k)=b

(5)

其中a為發(fā)展系數(shù),b為灰色作用量。

將式(3)帶入式(4),得到白化的微分方程

(6)

(7)

運(yùn)用最小二乘法對參數(shù)a和b進(jìn)行估計(jì),得到

(8)

通過式(8)計(jì)算出參數(shù)a和b的值,將其代入式(7),從而可得預(yù)測模型。

最后,對預(yù)測數(shù)據(jù)經(jīng)累減處理可得原始數(shù)據(jù)的預(yù)測公式

k=2,3,…,n

(9)

GM(1,1)模型常用的精度檢驗(yàn)方法有:殘差檢驗(yàn)法和后驗(yàn)差檢驗(yàn)法[10]。

1)殘差檢驗(yàn)法

殘差檢驗(yàn)是把預(yù)測數(shù)據(jù)和實(shí)際數(shù)據(jù)進(jìn)行比較,檢驗(yàn)每個(gè)點(diǎn)誤差的大小,一般用平均相對誤差作為檢驗(yàn)標(biāo)準(zhǔn)。

原始數(shù)據(jù)序列為X(0), 根據(jù)GM(1,1)灰色模型進(jìn)行預(yù)測,得到的擬合數(shù)據(jù)為(0)=((0)(1),(0)(2),…,(0)(n))。計(jì)算預(yù)測序列與原始序列的殘差序列為E=X(0)-(0)=(e(1),e(2),…,e(n)),其中計(jì)算相對誤差,即

(10)

GM(1,1)模型的平均相對誤差即為

(11)

2)后驗(yàn)差檢驗(yàn)法

(12)

(13)

計(jì)算的后驗(yàn)差比值如下

c=S2/S1

(14)

計(jì)算的小誤差概率如下

(15)

后驗(yàn)差檢驗(yàn)兩個(gè)比較重要的指標(biāo)分別是c和p。指標(biāo)c越小越好,c小就表明盡管原始數(shù)據(jù)很離散,但是經(jīng)過模型后數(shù)據(jù)序列的模擬值與實(shí)際值的差并不是太離散。指標(biāo)p越大越好,p越大,表明殘差與殘差平均值之差小于給定值0.6745S1的點(diǎn)較多,即擬合值(或預(yù)測值)分布比較均勻。

表1 模型精度等級

2.3 季節(jié)型ARIMA模型

單整自回歸移動(dòng)平均模型即ARIMA模型(AutoregressiveIntegratedMovingAverageModel),是由Box和Jenkins提出的一種時(shí)間序列預(yù)測模型[11],是一種精度較高的時(shí)間序列預(yù)測方法。

ARIMA模型的一般表達(dá)式為

Φ(L)Δdyt=Θ(L)εt

(16)

其中:εt為白噪聲序列,L是滯后算子,Δdyt=(1-L)dyt是通過差分得到的平穩(wěn)時(shí)間序列;Φ(L)=1-φ1L-φ2L2-…-φpLp和Θ(L)=1+θ1L+θ2L2+…+θqLq分別為自回歸算子和移動(dòng)平均算子多項(xiàng)式。

季節(jié)型單整自回歸移動(dòng)平均模型是以傳統(tǒng)ARIMA模型為基礎(chǔ),考慮進(jìn)去季節(jié)性因素調(diào)整的模型[12],即:

φp(L)ΦP(Ls)(1-L)d(1-Ls)Dyt=θq(L)ΘQ(Ls)εt

(17)

其中,φp(L)為非季節(jié)自回歸算子(AR)多項(xiàng)式;ΦP(Ls)為季節(jié)自回歸算子(SAR)多項(xiàng)式;θq(L)為非季節(jié)移動(dòng)平均算子(MA)多項(xiàng)式;ΘQ(Ls)為季節(jié)移動(dòng)平均算子(SMA)多項(xiàng)式;D、d分別為季節(jié)與非季節(jié)差分次數(shù);P、Q、p、q分別為季節(jié)自回歸、季節(jié)移動(dòng)平均、非季節(jié)自回歸、非季節(jié)移動(dòng)平均算子的最大滯后階數(shù)。周期為S的非平穩(wěn)季節(jié)性時(shí)間序列{yt}經(jīng)D次季節(jié)差分之后,即可建立季節(jié)型ARIMA(p,d,q)(P,D,Q)s模型。

建立季節(jié)型ARIMA模型有如下步驟:

1)對原序列進(jìn)行平穩(wěn)性檢驗(yàn),如果非平穩(wěn),則對序列進(jìn)行差分變換,直至其平穩(wěn)為止;

2)通過分析原序列或者差分序列的自相關(guān)函數(shù)(ACF)或者偏自相關(guān)函數(shù)(PACF),根據(jù)其截?cái)嗪屯衔蔡卣?確定模型的階數(shù),從而對ARIMA模型進(jìn)行識別;

3)對模型參數(shù)進(jìn)行估計(jì),檢驗(yàn)?zāi)Ｐ凸烙?jì)的殘差序列是否滿足隨機(jī)性要求;

4)優(yōu)化模型,綜合比較各種模型,根據(jù)不同篩選準(zhǔn)則,選擇最優(yōu)模型;

5)對未來數(shù)據(jù)進(jìn)行預(yù)測,根據(jù)最小均方誤差原則得到預(yù)測值。

3 建模預(yù)測方法

通過觀察HP濾波后的兩序列的特征,發(fā)現(xiàn)趨勢序列X帶有顯著的線性趨勢性,考慮建立GM(1,1)模型;對波動(dòng)序列C進(jìn)行相關(guān)性檢驗(yàn),序列C明顯帶有周期性和季節(jié)性,考慮建立季節(jié)型ARIMA模型。

具體步驟如下:

1)根據(jù)HP濾波,將原序列Y進(jìn)行分解,得到趨勢序列X與波動(dòng)序列C;

2)記分解后的趨勢序列X(0)=X,對序列X(0)建立GM(1,1)模型,利用最小二乘法對參數(shù)a和b進(jìn)行估計(jì),得到模型的擬合序列(0),相應(yīng)的擬合公式:

(18)

其中t表示時(shí)間期數(shù);

3)采用殘差檢驗(yàn)法和后驗(yàn)差檢驗(yàn)法對GM(1,1)模型進(jìn)行檢驗(yàn),檢驗(yàn)?zāi)Ｐ偷木?

5)檢驗(yàn)序列C的平穩(wěn)性,若序列平穩(wěn)則進(jìn)行相關(guān)性分析,否則對序列進(jìn)行差分運(yùn)算直至平穩(wěn),并記為序列{wc};

6)對序列{wc}進(jìn)行相關(guān)性分析,檢驗(yàn)序列的季節(jié)性,確定模型的階數(shù),建立季節(jié)型ARIMA模型

(19)

8)根據(jù)HP濾波原理,得到原序列的預(yù)測結(jié)果:

(20)

9)對模型預(yù)測結(jié)果進(jìn)行誤差分析,檢驗(yàn)的誤差類型為相對誤差(RelativeError,RE)。

(21)

組合模型的預(yù)測流程如圖1所示。組合模型預(yù)測的優(yōu)勢是:灰色GM(1,1)模型對近似線性指數(shù)趨勢有良好的預(yù)測效果,時(shí)間序列ARIMA模型對于反映細(xì)節(jié)隨機(jī)影響的波動(dòng)性有良好的預(yù)測作用。故將灰色預(yù)測GM(1,1)模型和季節(jié)型ARIMA預(yù)測模型結(jié)合起來分別對趨勢項(xiàng)和波動(dòng)項(xiàng)進(jìn)行預(yù)測,充分發(fā)揮兩種預(yù)測方法的長處,有利于提高預(yù)測的精度。

圖1 基于HP濾波的組合預(yù)測流程圖

4 實(shí)證分析

為了驗(yàn)證該預(yù)測方法,選取某衛(wèi)星行波管陽壓(陽極電壓)2013年1月～2021年7月的在軌數(shù)據(jù)作為研究對象,其中2013年1月～2020年12月的數(shù)據(jù)作為模型識別和建模的參考數(shù)據(jù);再用此模型預(yù)測2021年1月～2021年7月的數(shù)據(jù),并與同期實(shí)際數(shù)據(jù)比較并計(jì)量預(yù)測誤差。

陽壓數(shù)據(jù)隨時(shí)間變化趨勢如圖2所示。陽壓一直有緩慢上升趨勢,且每年陽壓有季節(jié)性變動(dòng):這是由于隨著行波管工作時(shí)間增加而出現(xiàn)內(nèi)部電極老化衰退,表現(xiàn)為陽壓有緩慢上升的趨勢;同時(shí)又受到軌道和姿態(tài)影響,陽壓還與行波管的溫度環(huán)境有關(guān),表現(xiàn)出一定的波峰波動(dòng)性(季節(jié)性變動(dòng))。

圖2 HP濾波序列分解圖

首先采用HP濾波法分解出趨勢項(xiàng)和波動(dòng)項(xiàng)。經(jīng)過反復(fù)試算擬合,綜合考慮趨勢的光滑程度并有利于波動(dòng)項(xiàng)建模,最終確定λ取值為14400。如圖2所示,序列X曲線代表長期趨勢,呈上升趨勢;序列C曲線代表了季節(jié)性和周期性波動(dòng)項(xiàng),具有規(guī)律性變化。

4.1 趨勢項(xiàng)預(yù)測

HP濾波分解后的趨勢序列X近似一條平滑的線性曲線,對其建立GM(1,1)模型。

記初始數(shù)據(jù)序列為X(0),X(0)=(x(0)(1),x(0)(2),…,x(0)(t))=(4835.197,4835.603,…,4848.70),其中t(t=1,2,3…)為時(shí)間期數(shù),對數(shù)據(jù)序列X(0)作一次累加運(yùn)算,構(gòu)造數(shù)據(jù)矩陣B以及數(shù)據(jù)向量Y,利用最小二乘法求解得到參數(shù)a=-0.0001,b=4835.3395。根據(jù)白化微分方程可得到GM(1,1)模型X(0)的時(shí)間響應(yīng)式

t=2,3,…,n

(22)

根據(jù)GM(1,1)檢驗(yàn)結(jié)果和式(22),對2021年1月至2021年7月的序列X進(jìn)行預(yù)測,得到模型的擬合和預(yù)測結(jié)果見圖3,預(yù)測結(jié)果計(jì)入表3?？梢?采用GM(1,1)模型對趨勢項(xiàng)進(jìn)行預(yù)測時(shí),對整體平滑走勢預(yù)測效果較好。

圖3 趨勢項(xiàng)擬合及預(yù)測值

4.2 波動(dòng)項(xiàng)預(yù)測

為建立波動(dòng)項(xiàng)ARIMA模型,對HP濾波后的波動(dòng)序列C進(jìn)行分析。為進(jìn)一步消除趨勢性,對序列C進(jìn)行一階差分,差分序列是平穩(wěn)的,故d=1;序列C也具有季節(jié)變化規(guī)律,周期為12,為了消除季節(jié)性,進(jìn)行滯后12期的季節(jié)差分處理,故D=1,S=12。

通過觀察季節(jié)差分序列的PACF和ACF來進(jìn)行模型的識別,可知p=12、q=4、P=0或1,Q=1。從中選擇最優(yōu)模型,調(diào)整后的擬合優(yōu)度R2最大的是ARIMA(12,1,4)(0,1,1)12,同時(shí)其AIC值和SC值也相對較小,因此較合適的模型是ARIMA(12,1,4)(0,1,1)12。此模型估計(jì)結(jié)果見表2。

表2 模型估計(jì)結(jié)果

在模型參數(shù)估計(jì)之后,需要對模型的殘差是否為白噪聲過程進(jìn)行檢驗(yàn)。對模型ARIMA(12,1,4)(0,1,1)12進(jìn)行擬合,圖4為季節(jié)差分序列的實(shí)際值、模擬值和殘差的曲線圖,可見模型的模擬值與實(shí)際值的變動(dòng)具有較好的一致性,殘差值較小,殘差比較平穩(wěn),初步判定其為白噪聲過程。對殘差進(jìn)行ADF單位根檢驗(yàn)后,可確定殘差序列為白噪聲過程。故ARIMA(12,1,4)(0,1,1)12模型為理想預(yù)測模型,可以用于預(yù)測。

圖4 季節(jié)差分序列的實(shí)際值、模擬值、殘差曲線圖

依據(jù)該模型和式(19),對波動(dòng)項(xiàng)2021年1月-7月的數(shù)據(jù)進(jìn)行預(yù)測,得到預(yù)測結(jié)果計(jì)入表3,預(yù)測曲線見圖5。可知,模型可以準(zhǔn)確的預(yù)測隨時(shí)間變化的波動(dòng)性,刻畫出波動(dòng)特性。預(yù)測置信區(qū)間很窄,并隨著預(yù)測期的擴(kuò)展,預(yù)測置信區(qū)間也變大,表明越往后模型的預(yù)測精度越差。

表3 陽壓預(yù)測結(jié)果

圖5 波動(dòng)項(xiàng)預(yù)測值及置信區(qū)間

4.3 組合預(yù)測及驗(yàn)證

依據(jù)HP濾波原理yt=xt+ct,利用上述模型進(jìn)行預(yù)測,得到 2021年1月-7月的陽壓預(yù)測結(jié)果并計(jì)入表3,與實(shí)際值比較,并進(jìn)行可靠性驗(yàn)證。

由表3中可知,預(yù)測值與實(shí)際值差值非常小,相對誤差都在0.04%之內(nèi),最小相對誤差僅為0.0017%,如圖6所示。預(yù)測結(jié)果非常準(zhǔn)確,達(dá)到了很高的精度;在半年預(yù)測期末,預(yù)測值依然很精確。因此,模型識別及驗(yàn)證效果理想,在行波管狀態(tài)不發(fā)生突變的情況下,可利用此模型對陽壓進(jìn)行較為準(zhǔn)確的中長期預(yù)測。

圖6 組合模型預(yù)測及相對誤差圖

陽壓遙測數(shù)據(jù)的原碼變化的最小單位是1個(gè)分層值,而1個(gè)分層值對應(yīng)計(jì)算得到的陽壓數(shù)值約為1.96,故該方法的預(yù)測精度(最大絕對誤差約1.85)已經(jīng)達(dá)到預(yù)測最小1個(gè)分層值的水平,精度非常高。由于陽壓受空間環(huán)境、自身老化衰退等影響,基于物理規(guī)律的預(yù)測精度往往不佳。而本方法的驗(yàn)證結(jié)果表明,不需要考慮復(fù)雜的物理和空間因素,僅僅利用數(shù)據(jù)自身的變化規(guī)律,能夠在中長期內(nèi)實(shí)現(xiàn)優(yōu)于單個(gè)分層值的高精度預(yù)測。

該組合預(yù)測方法可以推廣應(yīng)用于在軌故障診斷和預(yù)警中。在實(shí)際衛(wèi)星的在軌管理中,利用本方法對衛(wèi)星歷史和實(shí)時(shí)遙測數(shù)據(jù)進(jìn)行建模,可以發(fā)現(xiàn)理論預(yù)測值與實(shí)際值的差別,當(dāng)實(shí)際值偏離預(yù)測值較大時(shí),就可能診斷出相應(yīng)設(shè)備的潛在故障狀態(tài),同時(shí)可以用來對設(shè)備遠(yuǎn)期狀態(tài)進(jìn)行預(yù)警。

5 結(jié)束語

本文提出了一種基于時(shí)間序列分解的參數(shù)預(yù)測方法,通過對陽壓在軌數(shù)據(jù)的預(yù)測,證實(shí)了該預(yù)測方法的有效性,主要結(jié)論如下:

1)對HP濾波后的趨勢項(xiàng)采用灰色GM(1,1)模型預(yù)測,表明模型對線性趨勢預(yù)測效果良好;

2)對HP濾波后的波動(dòng)項(xiàng)采用季節(jié)型ARIMA模型預(yù)測,驗(yàn)證了模型能夠準(zhǔn)確模擬實(shí)際曲線的波動(dòng)性,跟蹤實(shí)時(shí)性較強(qiáng),預(yù)測效果良好;

3)在半年預(yù)測期內(nèi),預(yù)測值與實(shí)際值的相對誤差都在0.04%之內(nèi),表明該方法擬合精度高,預(yù)測精度高,適用于中長期預(yù)測。

應(yīng)用于在軌管理實(shí)踐中,該方法能夠?yàn)樾l(wèi)星狀態(tài)監(jiān)控、健康評估與故障預(yù)警提供新的思路和技術(shù)支撐,具有重要的工程應(yīng)用價(jià)值。