應用HP濾波的衛星遙測數據預測方法

李志強 張香燕 田華東

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

衛星遙測數據包含了星上各部件在軌運行的豐富信息,地面通常根據下傳的遙測數據對其運行情況進行監視和評估。研究分析在軌衛星的遙測數據變化規律,在此基礎上對遙測數據進行長期和短期預測,不僅有助于發現各部件和性能的老化衰退程度,還能根據預測數據進行提前預警和故障預示,以便提前發現非預期的異常變化并及時有效地處置,對降低衛星在軌運行風險、確保長期安全可靠運行具有重要意義[1-2]。

由于衛星面臨復雜的空間環境,影響遙測參數變化的空間因素較多,如輻射粒子、中性大氣、地磁等空間要素,還有運行軌道、姿態等影響,這些不確定或者周期性因素與衛星各部件相互作用產生效應,造成遙測數據呈現不同規律的波動,給遙測數據預測帶來了困難。遙測數據的變化規律是時間序列的各種統計特征的疊加,既有趨勢性,又有季節性,還有周期性的和不規則的變動,這些都反映了遙測數據變動規律的復雜性。研究時間序列數據的統計規律,并假定這種規律在未來時刻仍將發揮作用,就可以建立合適的時間序列數學模型來描述這種規律,并利用它作出合理的預測。

傳統上,對遙測數據進行預測一般采用幾種常見的預測方法,如曲線擬合、返向傳播(BP)神經網絡、小波分析、灰色系統預測等。曲線擬合方法一般適用于單一或者固定變化的簡單預測,且在長期預測時容易偏離原有的趨勢而發散[3];BP神經網絡預測依賴于大量樣本訓練和算法的有效性,且一般僅用于短期預測,對長期預測誤差較大[4];小波分析方法在進行剔除噪聲和提取周期處理時,容易遺漏數據突變和引起信息失真[5];灰色系統預測一般僅對單調變化具有較好的預測精度,對疊加了周期性和隨機性的數據預測精度不夠,且預測有效期較短[6-7]。HP(Hodrick-Prescott)濾波是一種常用的時間序列分析工具[8-9]?；跁r間序列分析原理,采用HP濾波將數據分解成相應時間序列的趨勢成分、波動成分,可弱化序列的趨勢性、波動性等因素之間的相互影響,用來分析不同序列的規律,同時保留原序列的完整數據特性,為進一步的數據預測提供了基礎。HP 濾波已在工業能源分析、經濟預測等領域應用,將HP濾波應用在衛星遙測數據分析中尚未看到相關報道。

本文將HP 濾波引入到遙測數據分析與預測中,提出應用HP濾波進行時間序列分解的遙測參數預測方法,即應用HP濾波對數據序列進行趨勢成分和波動成分的分解,分開建模、組合預測,最終通過疊加得到參數的預測值,具有較高的預測精度,可為遙測參數的分析和預測、故障診斷和預警提供一種新的技術方法。

1 預測方法

衛星遙測數據的變動規律一般既包含有趨勢性的長期規律,又包含有季節性、周期性和不規則變動的短期波動規律?；诖?本文提出的預測方法是:首先,應用HP濾波將衛星遙測數據Y 分解成趨勢成分數據G 和波動成分數據C;其次,根據趨勢成分和波動成分的不同特點,對兩者分別建模進行預測;最后,將2個預測結果進行疊加,得到原始遙測數據的組合預測結果。預測方法過程見圖1。對遙測數據應用HP濾波,并沒有消除原始遙測數據本身的趨勢性、周期性和季節性,這些性質被保留在分解出的趨勢成分G 與波動成分C 中。觀察兩類數據成分的特征,發現趨勢成分G 帶有顯著的線性趨勢性,因此考慮建立多元自回歸模型進行預測;對波動成分C進行檢驗,發現C 明顯帶有周期性和季節性,考慮建立季節型單整自回歸移動平均(ARIMA)模型進行預測。

圖1 預測方法過程Fig.1 Process of prediction method

1.1 方法步驟

預測方法的具體步驟如下。

(1)根據HP 濾波,將原始遙測數據Y＝{y1,y2,…,yn}進行分解,得到帶長期趨勢性的趨勢成分G＝{g1,g2,…,gn}與帶短期波動性的波動成分C＝{c1,c2,…,cn}[8-9]。它們與原始遙測數據的關系為

式中:t＝1,2,…,n,其中,n 為采樣的數據個數。

HP濾波目的在于將數據的趨勢成分從不平滑的原始遙測數據中分解出來,分離過程必須滿足損失函數最小原則,即

式中:平滑參數λ＝σ12/σ22,σ12和σ22分別為趨勢成分G 和波動成分C 的方差。

(2)將得到的趨勢成分G 視作t的函數,G 具有明顯的線性趨勢,定義各個數據點到自回歸曲線的誤差項為ε1,ε2,…,εn。根據min{ε12＋ε22＋…＋εn2}確定擬合多元自回歸多項式。

式中:G(t－i)表示滯后項的函數,即gt同時與gt－i的線性組合有關。

(3)根據擬合多項式,對趨勢成分G 進行預測,記預測值為。

(4)檢驗波動成分C 的平穩性:若C 平穩,則進行相關性分析;否則,對C 進行差分運算直至平穩,并記為差分數據{wC}。

(5)對差分數據{wC}進行相關性分析,檢驗數據的季節性,確定模型的階數,建立季節型ARIMA模型。ARIMA 模型是一種精度較高的時間序列預測方法,其應用于遙測數據yt的表達式為[10]

式中:εt為白噪聲;L 為滯后算子;通過差分得到的平穩時間序列Δdyt＝(1－L)dyt;Φ(L)＝1－φ1Lφ2L2－…－φpLp和Θ(L)＝1＋θ1L＋θ2L2＋…＋θqLq分別為自回歸算子和移動平均算子多項式,其中,φ 為自回歸系數,θ 為移動平均系數,p 和q 分別為非季節自回歸和移動平均算子的最大滯后階數。

針對遙測數據的季節性變動規律,考慮應用季節型ARIMA 模型[11-12],即

式中:φp(L)為非季節自回歸算子多項式;ΦP(Ls)為季節自回歸算子多項式;θq(L)為非季節移動平均算子多項式;ΘQ(Ls)為季節移動平均算子多項式;D 為季節差分次數;P,Q 分別為季節自回歸和移動平均算子的最大滯后階數;s 為周期。

式(5)即季節型ARIMA(p,d,q)(P,D,Q)s模型,依據這一模型可得到波動成分C 的預測公式為

(6)檢驗殘差數據{εt}是否為白噪聲,通過對殘差進行單位根擴展的迪克富勒(ADF)檢驗,當相應的t統計量小于一定置信水平的臨界值時,表明殘差序列不存在單位根,序列是平穩的,即為白噪聲過程;判斷季節型ARIMA 模型的合理性。

(7)根據季節型ARIMA 模型對波動成分C 進行預測,預測值記為。

(8)根據HP濾波原理,得到原始遙測數據的組合預測結果為

(9)對模型預測結果進行誤差分析,檢驗預測數據的相對誤差為

1.2 關鍵問題

本文預測方法需要解決的關鍵問題如下。

(1)HP濾波平滑參數λ 取值。λ 用來調節趨勢成分對原始遙測數據的跟蹤程度與趨勢成分的光滑程度。根據不同周期的數據和計算需求,λ 取值不同,隨著λ 值的增加,估計的趨勢更加平滑;當λ 無窮大時,估計的趨勢將接近線性函數。一般根據經驗和試算結果,綜合考慮趨勢的平滑程度并有利于波動成分ARIM A 建模來確定λ。

(2)季節型ARIMA 模型的模型識別及確定。通過分析波動成分的自相關函數和偏自相關函數的截尾性和拖尾性來識別模型的階數,并檢驗模型參數估計后的殘差,通過比較不同模型的調整擬合優度R2、赤池信息準則(AIC)值和施瓦茨準則(SC)值,篩選確定最優預測模型。其中:調整的擬合優度也稱為修正樣本決定系數,它表示總離差平方和中由回歸方程可以解釋的部分所占的比例,這一比例越大,表明回歸方程可以解釋的部分越多,模型越精確,回歸的效果越顯著;它一般介于0～1,越接近1說明回歸擬合效果越好。AIC 和SC 均是基于對數似然函數的統計量,它們的值越小越好,越小意味著模型的簡潔性和精確性都較好。選擇調整擬合優度R2較大、同時其AIC 值和SC 值較小的模型,即確定了最優預測模型。

本文預測方法根據遙測數據的趨勢成分G 和波動成分C 的不同特點分別建立模型,既有效地擬合了原數據的趨勢性,也降低了季節性對趨勢性的影響,從而建立了較為精確的組合預測模型。

2 預測方法實現及實例驗證

依據本文方法,設計遙測數據預測和預警系統方案,并進行實例驗證。

2.1 預測和預警系統方案

系統方案基于預測方法,采用分布式模塊化設計,由數據接收分發子系統、數據預測子系統、預警管理子系統、任務配置和中心控制子系統、故障診斷子系統及數據庫管理子系統組成,如圖2所示。用戶終端可接收任務配置和中心控制子系統發送的預測結果信息,用于在軌管理日常監視和預警。

圖2 系統組成Fig.2 System composition

數據預測子系統中是系統方案中的核心子系統,構成了預測方法的實現過程,主要包括時間序列分析模塊、預測算法配置模塊、實時預測模塊和歷史數據診斷模塊。其中:時間序列分析模塊配置HP濾波方法,對遙測數據進行趨勢成分和波動成分的分解;預測算法模塊接收來自時間序列分析模塊的分解數據,配置多項式自回歸算法模型和ARIMA模型,分別進行建模預測,并疊加二者預測值,最終得到組合預測結果;實時在軌遙測數據和歷史數據分別調用時間序列分析模塊和預測算法模塊,分別進行實時預測和歷史數據診斷,輸出預測和診斷信息。各模塊的結果發送給任務配置和中心控制子系統。

預警管理子系統獲取預測數據和提前預警閾值配置文件,并進行預測數據和預警閾值的比對,當預測結果超出預警閾值時,輸出相應的預警信息,并將信息發送給任務配置和中心控制子系統。故障診斷管理子系統獲取預測數據和歷史數據,配置診斷知識,并完成預測數據和歷史數據的比對,當預測結果與歷史數據規律不一致時,輸出相應的診斷信息,并根據診斷提示故障信息,同時將信息發送給任務配置和中心控制子系統。

任務配置和中心控制子系統完成整個系統的控制和運行管理,主要完成預測參數的選擇、預測任務的配置、平滑參數設定、ARIMA 模型識別、輸出結果管理等,尤其要解決預測方法實現中的兩個難點:一是時間序列分析模塊中的HP 濾波參數λ 的選擇,需要根據遙測數據分解結果和數據特點,反復迭代、綜合考慮來確定;二是ARIMA 模型的識別,在預測算法模塊計算數據自相關函數和偏自相關函數,通過分析其截尾性和拖尾性識別模型的階數,完成模型的確定和篩選。

2.2 實例驗證

為驗證預測方法的有效性,本文選取某衛星行波管陽壓(陽極電壓)2012年1月－2021年4月的在軌遙測數據用來驗證。其中:2012年1月－2020年10月的數據作為模型識別和建模的參考數據;再用此模型預測2020年11月－2021年4月的數據,并與同期實際數據比較并計算預測誤差。本文利用Matlab和Eviews軟件對數據進行建模分析。

2012－2021年陽壓隨時間變化趨勢如圖3 所示。陽壓隨時間一直有緩慢上升趨勢,且每年陽壓有季節性變動。這是由于隨著行波管工作時間增加,其內部電極有老化衰退現象,表現為陽壓有逐漸緩慢上升的趨勢;同時又受到運行軌道和光照影響,陽壓又與行波管的溫度環境有關,每年表現出一定的波峰波動性,即季節性變動。針對此遙測數據,首先采用HP濾波分解出趨勢成分和波動成分。濾波處理的關鍵是參數λ 的選擇,經過反復試算擬合,綜合考慮趨勢的光滑程度并有利于波動成分建模,最終確定λ 取值為14 400。

原始遙測數據經HP 濾波分解后如圖3所示,趨勢成分G 曲線代表長期趨勢,且增長率在壽命初期較為明顯,中后期較為平滑。波動成分C 曲線代表包含季節性、周期性的波動成分,每年波動情況不同。

圖3 遙測數據HP濾波分解結果Fig.3 Telemetry data separation results using HP filter

2.2.1 趨勢成分預測

趨勢成分G 近似一條平滑的曲線,記t(t＝1,2,3,…)為累積時間,依據式(3),將G 對t進行回歸擬合,根據多項式曲線擬合的結果調整次數和滯后項。選擇2012年1月－2020年10月數據擬合模型,用來預測2020年11月起未來6個月數據。依據點到擬合曲線的誤差平方和最小原則,得到模型估計系數計入表1,即建立了趨勢成分的預測模型。由表1可知:多項式自回歸曲線對G 的擬合效果非常好,R2和調整后的R2達到了1.000;AIC值和SC值分別為－9.714和－9.561,模型的精確性較好。各參數均通過了顯著性檢驗?？梢?趨勢gt不僅與時序t有關,還與2階滯后項有關。

表1 趨勢模型估計結果Table 1 Estimated results by trend model

根據模型估計結果和式(3),對2020年11月－2021年4月的趨勢成分G 進行預測,得到預測曲線如圖4所示。從圖4中可見:置信區間很窄,采用多項式自回歸曲線對趨勢成分進行預測時,對整體平滑走勢預測效果較好。

圖4 趨勢成分預測值及置信區間Fig.4 Prediction values and confidence interval of trend term

2.2.2 波動成分預測

通過建立季節性ARIMA 模型,對波動成分C進行分析與預測。經分析,HP 濾波后的波動成分C 自相關曲線沒有呈指數快速衰減,表明數據的趨勢性沒有完全消除。為進一步消除趨勢性,對C 進行1階差分,差分數據是平穩的,故d＝1。分析1階差分數據的自相關函數和偏自相關函數,發現C 也具有季節變化規律,周期為12;為了消除季節性,進行滯后12期的季節差分處理,故D＝1,s＝12。分析季節差分數據的偏自相關函數和自相關函數[10],可知p＝3,q＝3或1,P＝1,Q＝1,考慮建立模型ARIMA(3,1,1)(1,1,1)12或ARIMA(3,1,3)(1,1,1)12。在確定模型階數后,要篩選出最優模型。一般,調整的R2越大、AIC值和SC值越小,相應的模型越好。調整的R2較大的模型是ARIMA(3,1,1)(1,1,1)12,同時其AIC值和SC值也相對較小,因此較合適的模型是ARIMA(3,1,1)(1,1,1)12。模型估計結果見表2,其R2和調整后的R2分別為0.595和0.560,AIC值和SC值分別為2.960和3.171,相對其他模型較好,對波動能進行合適的刻畫。

表2 波動模型估計結果Table 2 Estimated results by fluctuation model

在模型估計之后,需要對模型的殘差是否為白噪聲進行檢驗,通過白噪聲檢驗才能進行預測。對殘差進行單位根ADF 檢驗后,可確定殘差為白噪聲。因此,ARIMA(3,1,1)(1,1,1)12模型為理想預測模型,可以用于預測。依據這一模型和式(6),對波動成分2020年11月－2021年4月的數據進行預測,得到預測曲線見圖5。

圖5 波動成分預測值及置信區間Fig.5 Prediction values and confidence interval of fluctuation term

圖5顯示模型可以準確地預測隨時間變化的波動性,并隨著預測期的延長,預測置信區間也變大,表明越往后模型的預測精度越差。這是由于動態模型中的AR 項和MA 項依賴于新的信息,而新信息的不確定性逐漸積累,從而導致新的預測精度下降。

2.2.3 組合預測及驗證

根據HP濾波原理yt＝gt＋ct,利用前述模型進行組合預測,得到2020年11月－2021年4月的陽壓預測結果,并與實際值比較,進行可靠性驗證?；疑碚擃A測是一種應用廣泛的數據預測方法,一般可應用于對趨勢性數據的預測。為驗證本文預測方法的有效性,利用灰色模型GM(1,1)對數據進行同期預測,結果與本文方法對比,如表3 所示。從表3中看出:利用組合模型進行預測,得到的預測結果非常精確,預測值與實際值差距非常小,絕對誤差在－1.6～－0.1,相對誤差都在0.04%之內,最小相對誤差僅為－0.003%,如圖6所示,在半年時間的預測區間內,達到了很高的精度。作為對比,灰色理論預測值絕對誤差明顯偏大。同時,本文預測方法獲得的預測值也準確模擬跟蹤了實際曲線的波動,跟蹤實時性較強,這主要是由于ARIMA 模型對波動細節和隨機性的刻畫程度高。相比較而言,灰色模型GM(1,1)模型的預測值僅僅反映了趨勢性增長,而對波動性預測不足,導致預測結果不能緊密跟蹤實際值變化,其預測絕對誤差也遠大于組合預測誤差。因此,本文預測方法的組合模型識別及驗證結果理想,用近10年的數據預測半年的數據,計算時間約5 min,在行波管狀態不發生突變的情況下,可利用此模型對行波管在軌陽壓數據進行較為快速準確的預測。

圖6 組合模型預測、灰色模型預測及相對誤差Fig.6 Combined model prediction,GM mode prediction and relative errors

表3 2020年11月-2021年4月陽壓預測結果Table 3 Prediction results of anode voltage from November 2020 to April 2021

事實上,陽壓遙測數據是由遙測原碼通過數據處理方法計算得到的,其原碼變化的最小單位是一個分層值,而一個分層值對應計算得到的陽壓數值約為2.00 V;因此該組合模型的預測精度(約1.65 V)已經達到預測一個分層值的微小變化,甚至先于遙測值更本質地反映其變化趨勢。由于陽壓受空間溫度環境、行波管自身性能衰退等諸多外界和自身因素的影響,趨勢成分和波動成分相互疊加,使得基于物理機制的預測得到的預測精度往往不佳。而本文預測方法的驗證結果表明:該方法不需要考慮多種復雜的物理和外界因素,僅僅依據數據本身的變化規律,在較長時間的跨度內能夠實現優于單個分層值的高精度預測,可以有效提升對遙測參數的預測能力。

本文預測方法還可以應用于在軌故障診斷和壽命預測中。衛星數據長期表現包含有各種因素共同作用的影響,評估性能時分離出各種成分的數據是必要的,可以更精準地發現數據變化的規律。在實際衛星的在軌監視中,利用本文預測方法對衛星歷史遙測數據進行回溯建模,可以發現參數理論預測值與實際值的差別,當實際值偏離預測值較大時,就可能診斷出相應設備的潛在故障狀態;同時,針對老化部件和性能衰退部件,利用此方法可以對設備長期性能進行診斷和預警,設置動態預警監視范圍,提前預警異常問題,進而實現對運行壽命的有效預測。

3 結束語

本文將HP濾波應用到衛星遙測數據領域,提出了一種基于HP濾波分解的遙測參數預測方法。將原始遙測數據應用HP濾波按照趨勢成分、波動成分分解,弱化了數據的趨勢性、波動性之間的相互影響,從而更精確地建立組合預測模型。通過在軌陽壓數據實例分析,驗證了該組合預測方法的有效性。該方法具有擬合精度高、預測精度高的特點,同時也能保證在長時間跨度內多步預測的精度?；谠摲椒ㄌ岢龅膽抿炞C方案,可在衛星故障診斷預警、壽命預測等方面實際應用。同時,將HP濾波思想應用到遙測數據,引入時間序列分析的成熟工具和模型,可為衛星狀態監視、健康評估與故障預警提供新的思路和技術方法,為在軌管理提供高效有力的分析手段,具有重要的工程應用價值。