小波分解與Prophet框架融合的電離層VTEC預(yù)報模型

田 睿,董緒榮

(航天工程大學(xué)航天信息學(xué)院,北京 101407)

0 引 言

在導(dǎo)航定位、無線通信、航空航天、氣象預(yù)測等領(lǐng)域,電離層對電磁波的折射、散射、反射和吸收效應(yīng)影響巨大[1-5]。在全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system,GNSS)領(lǐng)域,電離層垂直總電子含量(vertical total electron content,VTEC)是直接決定電離層延遲誤差的重要參數(shù),尤其在單頻精密定位的實時應(yīng)用中,須構(gòu)建高精度的電離層延遲預(yù)報模型對電離層延遲進(jìn)行實時修正。因此,研究電離層VTEC預(yù)報具有重要的應(yīng)用價值[6]。目前,應(yīng)用比較廣泛的電離層模型大多是傳統(tǒng)經(jīng)驗?zāi)Ｐ?如IRI、Klobuchar、Bent等模型。但傳統(tǒng)經(jīng)驗?zāi)Ｐ偷膽?yīng)用效果并不十分理想[7],如常用的Klobuchar模型的預(yù)報精度僅為50%～60%。因此,國內(nèi)外學(xué)者提出了多種電離層VTEC預(yù)報方法,如時間序列、神經(jīng)網(wǎng)絡(luò)等[8-14]。其中,相比于其他預(yù)報方法,時間序列法具有短期預(yù)報精度高、計算簡單、理論完備、樣本數(shù)據(jù)要求較少等優(yōu)勢[15]。因此,在電離層短期預(yù)報領(lǐng)域,精度較高且相對簡單的時間序列法受到了廣泛關(guān)注。李秀海等學(xué)者最早引入了自回歸(autoregressive,AR)模型來構(gòu)建電離層VTEC時序預(yù)測模型,但AR模型難以準(zhǔn)確地擬合VTEC時序數(shù)據(jù)的周期性變化,預(yù)測精度較低[16]。相比于AR模型,差分AR移動平均(AR integrated moving average,ARIMA)模型對時間序列的周期性變化及趨勢項擬合較好,更適用于包含季節(jié)性變化及短期趨勢項的電離層VTEC時序數(shù)據(jù),其預(yù)測精度明顯提升,因此目前大部分研究均以ARIMA模型為基礎(chǔ)[15,17-19]。然而,由于電離層受地磁擾動、太陽活動、日地相對距離等環(huán)境因素影響較大[20],特別是在磁暴期間,受到強(qiáng)地磁擾動的影響,電離層VTEC的周日變化更為顯著[21]。而在采用ARIMA模型進(jìn)行預(yù)報時,受磁暴等復(fù)雜因素影響,電離層VTEC時序數(shù)據(jù)的非平穩(wěn)性特征與非線性特征顯著增強(qiáng),大幅降低了ARIMA模型的預(yù)報精度。對此,翟篤林等學(xué)者[22]基于2017年Facebook開源的Prophet框架[23]進(jìn)行了電離層VTEC的短期預(yù)報與異常探測,實驗證明基于Prophet框架的預(yù)測精度明顯優(yōu)于ARIMA模型。相比于其他時間序列預(yù)測方法,Prophet框架不僅有較理想的預(yù)測精度,還具有:① 無需復(fù)雜的特征工程;② 可自動處理缺失值;③ 支持大規(guī)模細(xì)粒度數(shù)據(jù);④ 調(diào)參簡單;⑤ 考慮了節(jié)假日效應(yīng)與特殊事件的影響;⑥ 可解釋性強(qiáng)等一系列優(yōu)點(diǎn)。然而,如何將Prophet框架與其他算法進(jìn)一步融合,提高預(yù)測精度,是當(dāng)前Prophet框架研究中亟待解決的問題[23]。

為進(jìn)一步改進(jìn)基于Prophet框架的電離層VTEC短期預(yù)報模型,提高預(yù)測精度,本文將小波分解與Prophet框架相融合,并綜合考慮了地磁擾動對電離層的影響,提出一種小波分解與Prophet框架融合的時間序列預(yù)報模型以進(jìn)行電離層VTEC短期預(yù)報。同時采用國際GNSS服務(wù)(international GNSS service，IGS)發(fā)布的GIM格網(wǎng)數(shù)據(jù)進(jìn)行對比實驗,驗證改進(jìn)模型的預(yù)報精度與適用性。

1 小波分解與Prophet框架融合的時間序列預(yù)報模型

為更詳盡地闡述本文所提改進(jìn)模型,將在第1.1節(jié)～第1.4節(jié)中首先對相關(guān)概念進(jìn)行簡要介紹,并在第1.5節(jié)中詳細(xì)地闡述本文改進(jìn)模型的基本架構(gòu)與具體步驟。

1.1 小波分解基本原理

(1)

則Vj中任意函數(shù)fj均存在如下的多分辨表示:

(2)

可通過Mallat塔式算法對一維離散時間信號(序列)進(jìn)行小波分解,分解過程的表達(dá)式為

(3)

X=Aj+Dj+Dj-1+…+D1

(4)

1.2 Prophet框架

Prophet框架的基礎(chǔ)模型是一個由季節(jié)項、趨勢項、節(jié)假日或特征事件影響(節(jié)假日效應(yīng))、誤差項4部分組成的時間序列廣義加法模型(可通過用戶設(shè)置調(diào)整為乘積季節(jié)性模型),廣義加法模型的表達(dá)式為

y(t)=g(t)+s(t)+h(t)+εt

(5)

式中,y(t)為時間序列在時刻t的取值;g(t)為趨勢項,用于擬合時間序列的非周期變化;s(t)為季節(jié)項(或稱周期項),用于擬合時間序列各種周期性變化(應(yīng)注意,當(dāng)設(shè)置為乘積季節(jié)性模型時,s(t)為對數(shù)形式);h(t)表示節(jié)假日效應(yīng)特征事件的影響,通過該項,節(jié)假日或特征事件的影響可作為先驗信息融入到模型中;εt為誤差項(或稱噪聲因子),假設(shè)其服從正態(tài)分布。Prophet框架僅以時間作為自變量,且無需復(fù)雜的特征工程即可得到趨勢項、季節(jié)項及節(jié)假日效應(yīng)等組分?？山忉屝詮?qiáng),明確解釋了目標(biāo)序列的時間依賴結(jié)構(gòu)[25]。

1.2.1 趨勢項模型

Prophet框架中有兩種趨勢項預(yù)測模型:飽和增長模型與分段線性模型。飽和增長模型基于邏輯回歸函數(shù)進(jìn)行預(yù)測,其預(yù)測結(jié)果的增長趨勢與人口增長的趨勢類似,適用于預(yù)測一定承載能力下的非線性飽和增長。所謂承載能力即增長可到達(dá)的最大極限值。飽和增長模型的表達(dá)式如下：

(6)

式中,C(t)為承載能力;k為增長率;m為偏移量。應(yīng)注意該模型與普通的邏輯回歸函數(shù)有兩處不同:① Prophet框架中的承載能力C(t)不是一個常數(shù),而是隨時間遷移而變化,需要人為給定該參數(shù)的值;② 增長率k也并非常數(shù),Prophet框架通過給定變異點(diǎn)在模型中引入增長趨勢的變化。

Prophet框架中的變異點(diǎn)位置可由用戶根據(jù)先驗信息人為指定,也可由Prophet框架自動選取。默認(rèn)設(shè)置下,Prophet首先確定大量潛在的變異點(diǎn),然后對潛在變異點(diǎn)上趨勢變化的幅度做稀疏先驗(等同于L1正則化)來選取變異點(diǎn)。實際上Prophet在建模時會識別出很多增長率k發(fā)生突變的潛在變異點(diǎn),但會盡可能少地使用。按照默認(rèn)設(shè)置,Prophet會在前80%的時間序列數(shù)據(jù)中識別出25個變異點(diǎn)。

假設(shè)Prophet在序列中確定了S個變異點(diǎn),變異點(diǎn)位置在時間戳sj(j=1,2,…,S)上,也即在時間戳sj上增長率k發(fā)生改變。定義增長率變化向量δ∈RS,δ中的元素δj表示時間戳sj處增長率的變化。設(shè)序列初始增長率為k,某時間戳上sj的增長率可表示為

(7)

式中,αj(t)為向量α(t)中的元素,則在任意時刻t,k(t)可表示為k+αT(t)δ,進(jìn)而可定義向量α(t)中的元素為

(8)

由于變異點(diǎn)導(dǎo)致了趨勢項的非連續(xù)性,須對偏移量m進(jìn)行自適應(yīng)調(diào)整。則在變異點(diǎn)對應(yīng)的時間戳sj處,通過下式對偏移量進(jìn)行調(diào)整:

(9)

此時偏移量調(diào)整為m+αT(t)γ,并最終得到飽和增長模型:

(10)

用于表示時間序列的趨勢項。而許多時間序列的趨勢項并不符合飽和增長趨勢,對此,Prophet框架采用更簡約的分段線性模型對其趨勢項進(jìn)行擬合:

g(t)=[k+αT(t)δ]t+[(m+αT(t)γ)]

(11)

應(yīng)注意,增長率變化向量δ滿足δ～Laplace(0,τ)。其中,參數(shù)τ是本文主要調(diào)整的超參數(shù),因為其直接控制趨勢項增長率變化的靈活性。增大該參數(shù)會使趨勢擬合更加靈活,但存在過擬合風(fēng)險。

1.2.2 季節(jié)周期性

為對時間序列的季節(jié)性(或稱周期性)進(jìn)行精確建模,Prophet框架采用離散傅里葉級數(shù)對時間序列的季節(jié)項進(jìn)行建模:

(12)

定義向量X(t)為

X(t)=

(13)

則可將季節(jié)項s(t)表示為X(t)與一個參數(shù)向量β的點(diǎn)乘形式為

s(t)=X(t)β

(14)

式中,β用于對模型季節(jié)性進(jìn)行平滑,起到類似L1正則化的作用,其服從正態(tài)分布,即β～Normal(0,σ2)?？赏ㄟ^增大參數(shù)σ的值以適應(yīng)更大的季節(jié)性波動,較小的σ則會抑制季節(jié)性效應(yīng),其默認(rèn)值設(shè)置為10,該值在一般問題中通常是適用的[23]。在默認(rèn)設(shè)置下,Prophet框架的參數(shù)估計方法為最大后驗概率估計(簡稱為MAP),僅能得出趨勢不確定性及觀測噪聲的影響,用戶可通過設(shè)置mcmc.samples參數(shù),采用馬爾可夫鏈蒙特卡羅取樣得到季節(jié)的不確定性。

1.2.3 節(jié)假日或特殊事件影響

考慮到節(jié)假日或特殊事件對時序數(shù)據(jù)的影響,Prophet框架在模型中將節(jié)假日效應(yīng)h(t)作為先驗知識融入到模型中,并認(rèn)為節(jié)假日效應(yīng)對時間序列的影響是獨(dú)立的。應(yīng)注意,節(jié)假日或特殊事件需要作為先驗信息由用戶給定。設(shè)節(jié)假日或特殊事件i對應(yīng)的日期列表為Di,如十一黃金周對應(yīng)的日期列表Di中包含10月1日到10月7日7個日期。Prophet框架通過一個指示函數(shù)表示某時刻t是否處于節(jié)假日或特殊事件i期間,同時需要一個參數(shù)κi來表示節(jié)假日及特殊事件的影響,設(shè)節(jié)假日或特殊事件i共包含L天,則節(jié)假日效應(yīng)可表示為

(15)

式中,矩陣Z(t)表示為[1{t∈D1},1{t∈D2},…,1{t∈DL}];矩陣κ表示為[κ1,κ2,…,κL]T,κ～Normal(0,υ2),其中υ默認(rèn)值設(shè)定為10,該值取值越大表示節(jié)假日對模型的影響越大;該值越小表示節(jié)假日對模型影響越小,用戶可根據(jù)先驗知識調(diào)整該值。

1.3 分解-集成模型

在時間序列預(yù)測中,常采用分解時間序列的方法對現(xiàn)有模型進(jìn)行改進(jìn),即構(gòu)建分解-集成模型[26],以提升預(yù)報精度。常采用的方法有STL分解、EMD分解等[27]。其中,小波分解作為分析非線性及非平穩(wěn)信號的重要數(shù)學(xué)工具[28-29],適合處理具有非線性、非平穩(wěn)特征的電離層VTEC時間序列,在對ARIMA模型的改進(jìn)中取得了良好的效果,如劉立龍等學(xué)者[30]及鮑亞東等學(xué)者[31]均引入了小波分解以改進(jìn)ARIMA模型,提高了電離層VTEC短期預(yù)報精度。

基于小波分解構(gòu)建的分解-集成模型能提高預(yù)報精度的原因在于:小波分解實現(xiàn)了時間序列的時頻局部化,充分挖掘了時間序列中包含的信息[32]。因此,可通過小波分解有效地分離和提取時序數(shù)據(jù)的周期性、非線性及變化趨勢[33],使預(yù)測模型能夠更好地對時序數(shù)據(jù)的周期性變化、變點(diǎn)信息及趨勢項進(jìn)行擬合與建模,從而得到更為精確的預(yù)測結(jié)果。

而電離層VTEC時序數(shù)據(jù)可視作非平穩(wěn)非線性的離散時間信號(序列),通過小波分解對其進(jìn)行多分辨分析。分解所得的低頻分量包含了信號的主體信息,而不同分辨率的高頻分量則描繪了信號的細(xì)節(jié)紋理。因此,可利用小波分解快速高效地對電離層VTEC時間序列的各分量進(jìn)行分離,使樣本序列的周期性變化、變點(diǎn)信息和短期趨勢更加顯著,預(yù)測模型能夠更好地對其擬合與建模,從而提高預(yù)測精度[30]。

1.4 小波基的對比及選取

小波基的選取問題一直是小波技術(shù)應(yīng)用中的難題,眾多學(xué)者針對各自的應(yīng)用需求提出了不同的選取方法[34-36],本文將首先對各種小波基的特性進(jìn)行分析,并結(jié)合電離層VTEC時序數(shù)據(jù)的特點(diǎn),通過理論分析與實驗相結(jié)合的方式進(jìn)行小波基的選取。

一般而言,小波基具有5個重要特性[37]:正交性(或雙正交性)、對稱性、正則性、消失距和緊支撐性。正交性反映了小波基的完善程度,規(guī)范的正交性有利于信號的精確重構(gòu);較好的對稱性可避免信號分解與重構(gòu)時的相位失真;正則性表征小波基的可微性,較好的正則性有利于捕獲信號的奇異點(diǎn),大部分正交小波基正則性越高則消失距越高[38];消失距反映了小波變換后的能量集中程度,支撐寬度反映了小波的局部化能力,支撐越小,小波基局部化能力越強(qiáng)[39],而支撐越大,正則性越好。常用小波基的各項特性如表1所示。

表1 常用小波基的各項特性Table 1 Characteristics of commonly used wavelet bases

結(jié)合前文所述的電離層VTEC時間序列的特點(diǎn)及小波分解的作用,適合本文應(yīng)用場景的小波基應(yīng)滿足如下要求:① 規(guī)范的正交性,應(yīng)選擇具有規(guī)范正交性的小波基,有利于小波分解后的精確重構(gòu)。② 較高的消失距,如前文所述,較高的消失距即意味著較好的正則性。這有利于捕獲序列中的奇異點(diǎn)(有利于Prophet框架中的變異點(diǎn)建模),充分挖掘電離層VTEC時間序列包含的信息,便于擬合與建模。③ 適中的支撐長度,如前文所述,支撐越小,小波基局部化能力越強(qiáng),而支撐越大,正則性越好。

基于上述要求,本文排除了Haar、mever、morlet等小波基,并將在Daubechies、Symelets、Coiffets 3種小波基中選取合適的小波基。

如前文所述,小波分解提高預(yù)報精度的原因可簡要概括為:小波分解可充分挖掘序列包含的信息,使樣本序列的周期性變化、變點(diǎn)信息或短期趨勢更加顯著,利于建模和預(yù)測。因此,為對上述小波基的分解效果進(jìn)行對比,本文基于最大信息熵原理對分解效果進(jìn)行評估,該原理廣泛應(yīng)用于時間序列的分析與預(yù)測中[40-42]。最大信息熵原理指出:在滿足所有已知條件的情況下進(jìn)行預(yù)測,應(yīng)選擇信息熵最大的可行解,這樣所得的解最為客觀、超然且偏差最小[43]。楊薛明等學(xué)者即應(yīng)用此原理,引入信息熵作為優(yōu)選樣本序列的依據(jù)[42]。本文計算了上述小波基分解所得子序列的信息熵,并基于最大信息熵原理進(jìn)行小波基的選取。

實驗過程、結(jié)果及分析詳見第2.2節(jié),本文最終選用db4小波基,取得了良好的改進(jìn)效果,這也與文獻(xiàn)[30-31]的結(jié)論相印證。

1.5 改進(jìn)模型的架構(gòu)與步驟

本文提出了一種小波分解與Prophet框架融合的時間序列預(yù)報模型以進(jìn)行電離層VTEC短期預(yù)報,其總體架構(gòu)如圖1所示。

圖1 改進(jìn)模型的總體架構(gòu)Fig.1 Overall architecture of the improved model

如圖1所示,本文所提改進(jìn)模型的基本思路是:首先按第1.1節(jié)中的方法對電離層VTEC時間序列進(jìn)行小波分解,得到近似分量Aj和一組細(xì)節(jié)分量D1,D2,…,Dj,并分別基于Prophet框架對其進(jìn)行預(yù)測。最后,對預(yù)報序列進(jìn)行重構(gòu)得到最終預(yù)測結(jié)果。該模型是一個典型的分解-集成模型[26]。

為進(jìn)行電離層VTEC的期預(yù)報,除預(yù)測模型外,還必須考慮樣本數(shù)據(jù)的采集與預(yù)處理、超參數(shù)的調(diào)節(jié)、模型的性能度量等方面的內(nèi)容。因此,對基于本文改進(jìn)模型進(jìn)行電離層VTEC短期預(yù)報的具體步驟如下。

步驟 1分別選取平靜期和活躍期時段,下載IGS發(fā)布的GIM格網(wǎng)數(shù)據(jù)文件(IONEX格式文件)。同時,從國家空間科學(xué)數(shù)據(jù)中心上下載對應(yīng)時段的地磁指數(shù)數(shù)據(jù),并從空間環(huán)境預(yù)報中心上下載相應(yīng)的歷史預(yù)報產(chǎn)品。

步驟 2選擇格網(wǎng)點(diǎn),參考日本東京海洋大學(xué)開發(fā)的開源軟件RTKLIB中的相關(guān)源碼,按照一定的時間粒度,編程提取目標(biāo)格網(wǎng)點(diǎn)的VTEC值,獲得原始時間序列數(shù)據(jù)。

步驟 3由于GIM格網(wǎng)數(shù)據(jù)具有可靠的精度和穩(wěn)定性,常作為基準(zhǔn)數(shù)據(jù)使用[44]。因此,無需處理缺失值與異常值。應(yīng)注意,其中明顯區(qū)別于一般數(shù)據(jù)的VTEC值并非錯誤值,不應(yīng)作為“異常值”處理,最終得到實驗數(shù)據(jù)集。

步驟 4將實驗數(shù)據(jù)集輸入模型,同時根據(jù)下載的地磁指數(shù)數(shù)據(jù)以及歷史預(yù)報產(chǎn)品,構(gòu)建特殊事件(如磁暴等)列表作為先驗信息輸入模型,從而在模型中引入地磁擾動的影響。以均方根誤差(root mean square error,RMSE)、平均絕對誤差(mean absolute error,MAE)和平均百分比誤差(mean absolute percentage error,MAPE) 3項指標(biāo)評估預(yù)測結(jié)果，并綜合采用網(wǎng)格搜索法自動調(diào)參、可視化技術(shù)實時交互式調(diào)參等方法調(diào)整模型參數(shù)。

步驟 5最后經(jīng)過多次迭代調(diào)參,可獲得較優(yōu)的最終模型,得到預(yù)測結(jié)果,并對最終模型的預(yù)測結(jié)果進(jìn)行評估。

在特殊事件列表的構(gòu)建過程中,已知時段的特殊事件通過國家空間科學(xué)數(shù)據(jù)中心上下載的地磁指數(shù)數(shù)據(jù)來確定,而預(yù)報時段的特殊事件則通過空間環(huán)境預(yù)報中心上下載的歷史預(yù)報產(chǎn)品確定。這與實際應(yīng)用場景相符合,如已知11天的電離層VTEC數(shù)據(jù),預(yù)測步長為3天,即以11天的已知數(shù)據(jù)作為實驗數(shù)據(jù)集進(jìn)行預(yù)測。在這11天的已知時段上,可通過已知的地磁指數(shù)數(shù)據(jù)來確定特殊事件列表,而在未來3天的預(yù)報時段上,則僅能通過預(yù)報產(chǎn)品確定特殊事件列表。

2 實驗及分析

2.1 原始數(shù)據(jù)及先驗信息的獲取與分析

GIM格網(wǎng)數(shù)據(jù)的時間分辨率為2 h,一天的數(shù)據(jù)包含13張全球電離層格網(wǎng)地圖,空間分辨率為2.5°(緯度)×5°(經(jīng)度)。參照文獻(xiàn)[30]中的實驗方法,本文的實驗區(qū)域集中于中國以及周邊國家和地區(qū),如表2所示。

表2 實驗區(qū)域Table 2 Experimental area

本文按照經(jīng)度間隔10°、緯度間隔7.5°選取格網(wǎng)點(diǎn),共計54個格網(wǎng)點(diǎn),并在單個格網(wǎng)點(diǎn)上按每天時段為00:00～22:00,以2 h為間隔提取VTEC數(shù)據(jù)。實驗數(shù)據(jù)的時段選擇上,根據(jù)文獻(xiàn)[15],本文選取2010年2月16日至3月1日(年積日47～60日)進(jìn)行活躍期電離層預(yù)報實驗;選取2008年2月1日至2月14日(年積日32～45日)進(jìn)行平靜期電離層預(yù)報實驗。參照文獻(xiàn)[15]及文獻(xiàn)[30]的實驗方法,以每個時段前11天數(shù)據(jù)作為實驗用數(shù)據(jù)集,預(yù)測步長為3天,并以GIM格網(wǎng)數(shù)據(jù)為基準(zhǔn)對預(yù)報結(jié)果進(jìn)行評估。

如前文所述,特殊事件列表將作為先驗信息輸入模型,并選取地磁指數(shù)Ap作為衡量地磁活動指標(biāo)。其中地磁指數(shù)Ap表示行星等效日幅度,可作為全天地磁活動水平的度量[45]。為說明磁暴等特殊事件的確定方法,獲取建模所需的先驗信息,對活躍期時段的地磁擾動強(qiáng)度進(jìn)行分析。活躍期時段擾動暴實時(disturbance storm time,DST)與地磁指數(shù)Ap的變化,如圖2所示。

圖2 年積日47～60天的地磁指數(shù)Ap與DST指數(shù)Fig.2 Geomagnetic Ap and DST Index (DOY 47～60)

DST指數(shù)表征環(huán)電流強(qiáng)度,當(dāng)DST指數(shù)小于-30 nT時即可能發(fā)生小磁暴,小于-50 nT時即可能發(fā)生中等磁暴[30]。如圖2所示,年積日47～48日多個歷元DST指數(shù)較低(可能發(fā)生中小磁暴),與之對應(yīng)的地磁指數(shù)Ap均達(dá)到或超過5。即可根據(jù)地磁指數(shù)Ap判斷某日的地磁擾動是否劇烈,并構(gòu)建特殊事件列表,作為先驗信息輸入模型。根據(jù)圖2,預(yù)測時段內(nèi)地磁擾動并不劇烈,預(yù)測中無需考慮特殊事件的影響。

圖3 活躍期格網(wǎng)點(diǎn)(32.5°N,75°E)的原始樣本序列小波分解結(jié)果Fig.3 Wavelet decomposition results of the original sample sequence of lattice nodes (32.5°N,75°E) in the active period

如前文所述,本文改進(jìn)模型還需要對輸入的樣本序列進(jìn)行小波分解,以往研究常采用1～3級分解[12,28-31],分解層數(shù)過多會造成累積誤差過大。本文經(jīng)實驗證明,一級分解的預(yù)報結(jié)果精度最高,且編程實現(xiàn)簡單,更重要的是避免了多層分解中誤差的累積,這也與文獻(xiàn)[30]中的結(jié)論相印證[30]。限于篇幅,圖3僅給出了活躍期時段格網(wǎng)點(diǎn)(32.5°N,75°E)的原始樣本序列小波分解的結(jié)果。

2.2 小波基選取實驗

根據(jù)第1.4節(jié)的有關(guān)討論,本文將基于GIM格網(wǎng)數(shù)據(jù),對sym4、db4、sym8、db8、sym12、db12、coif2和coif3等小波基進(jìn)行實驗對比,并從中優(yōu)選適用的小波基。

本實驗基于上述小波基,依次對不同時期不同格網(wǎng)點(diǎn)對應(yīng)的原始樣本序列進(jìn)行小波分解。然后,計算子序列的信息熵,并基于最大信息熵原理進(jìn)行小波基選取。由于對電離層活躍期數(shù)據(jù)進(jìn)行預(yù)報的難度較大,選取時將優(yōu)先考慮活躍期的性能。不同時期不同實驗區(qū)域的分解序列信息熵均值,如圖4所示。分解所得的序列信息熵差異主要表現(xiàn)在高頻分量上,這是因為信號的細(xì)節(jié)紋理主要包含在高頻分量中[30]。由圖4可見,db4小波基分解所得的子序列信息熵均值在不同時期、不同實驗區(qū)域內(nèi)均較高,特別是在中緯度區(qū)域顯著高于其他小波基。則根據(jù)第1.4節(jié)的有關(guān)討論,基于最大信息熵原理,本文擬選用db4小波基。

圖4 不同時期不同實驗區(qū)域的分解序列信息熵均值Fig.4 Mean information entropy of the decomposition sequences in different experimental regions and different periods

2.3 模型精度及適用性分析

2.3.1 模型評估指標(biāo)

如前文所述,采用RMSE、MAE和MAPE 3個性能度量指標(biāo)對模型結(jié)果進(jìn)行評估。

RMSE計算公式為

(16)

MAE計算公式為

(17)

式中,MAE表示實際輸出值與預(yù)測值絕對差值的平均值。MAE取值越小,說明模型精度越高。

MAPE計算公式為

(18)

式中,MAPE取值越小,說明模型精度越高。

2.3.2 模型參數(shù)設(shè)置

根據(jù)電離層VTEC時間序列的性質(zhì)[15,30]可知,在短期預(yù)報中,序列的周期性變化以比較規(guī)律的周日變化為主,序列的周期無明顯改變,周期的不確定性并不顯著。因此,與周期有關(guān)的模型參數(shù)均設(shè)置為Prophet框架的推薦值,上述推薦值在一般問題中通常是適用的[23]。同時,在變異點(diǎn)篩選中,采用Prophet框架的自動篩選[23],這一方面提高了建模的客觀性,另一方面降低了建模難度。然而,在電離層VTEC時間序列的短期預(yù)報中,序列的趨勢項變化較為劇烈,特別是在電離層活躍期的中低緯度地區(qū)。如前文所述,超參數(shù)τ直接控制趨勢項增長率變化的靈活性。增大該參數(shù)會使趨勢擬合更加靈活,但也存在過擬合風(fēng)險。因此,本文采用網(wǎng)格搜索法,基于外推結(jié)果的RMSE值對參數(shù)τ進(jìn)行調(diào)優(yōu)。如前文所述,在本文提出的改進(jìn)模型中,需利用Prophet框架對高頻分量與低頻分量分別進(jìn)行預(yù)測。對兩種分量的預(yù)測需要設(shè)置不同的參數(shù)τ。因此,采用網(wǎng)格搜索法對改進(jìn)模型進(jìn)行調(diào)參,參數(shù)取值范圍設(shè)置為0.05～0.95,搜索步長設(shè)置為0.05,結(jié)果如圖5所示?；谝陨暇W(wǎng)格搜索結(jié)果并進(jìn)行人工微調(diào)后,實驗中改進(jìn)模型的參數(shù)設(shè)置如表3所示。

圖5 改進(jìn)模型超參數(shù)τ的網(wǎng)格搜索結(jié)果Fig.5 Grid search results for the hyperparameter τ of the improved model

表3 參數(shù)設(shè)置Table 3 Parameter settings

2.3.3 實驗結(jié)果分析

按照前文所述實驗方法,分別采用本文改進(jìn)模型、未改進(jìn)的Prophet框架、ARIMA模型進(jìn)行對比實驗,并對預(yù)測結(jié)果進(jìn)行統(tǒng)計,即可計算RMSE、MAE和MAPE 3項指標(biāo)。在電離層平靜期與活躍期,不同實驗區(qū)域內(nèi)各項指標(biāo)的均值分別如表4和表5所示。

表4 電離層平靜期不同實驗區(qū)域內(nèi)各項指標(biāo)的均值Table 4 Mean values of each index in different experimental areas during the ionosphere quiet period

表5 電離層活躍期不同實驗區(qū)域內(nèi)各項指標(biāo)的均值Table 5 Mean values of each index in different experimental areas in the ionosphere active period

如表4和表5所示,無論在電離層平靜期還是活躍期,本文改進(jìn)模型在低、中、高緯度3個實驗區(qū)域的預(yù)測誤差均比較小,且各項評估指標(biāo)的均值優(yōu)于未改進(jìn)的Prophet框架,這驗證了本文改進(jìn)的有效性。同時,本文改進(jìn)模型與未改進(jìn)的Prophet框架明顯優(yōu)于傳統(tǒng)的ARIMA模型,這也與文獻(xiàn)[22]中的結(jié)論相印證[22]。

表4與表5進(jìn)行對比可知,在電離層活躍期,3種模型的預(yù)測精度均有所下降,各項指標(biāo)大多明顯增加。這是因為在電離層活躍期,受到地磁活動以及日地相對距離、太陽活動等復(fù)雜的不確定因素影響,電離層VTEC周日變化更加劇烈,導(dǎo)致時間序列的非平穩(wěn)性與非線性特征顯著增強(qiáng),造成建模難度增大,預(yù)報精度受到嚴(yán)重影響。類似的,由于VTEC的周日變化幅度隨緯度降低而總體增大,其非平穩(wěn)性非線性特征更加顯著,由表4和表5可知,低緯度地區(qū)預(yù)報精度普遍低于中、高緯度地區(qū),且改進(jìn)的效果相對較差,這是因為序列的非平穩(wěn)性受各項復(fù)雜因素影響顯著增大,增加了分解難度與累積誤差。總體而言,本文改進(jìn)模型在中、高緯度地區(qū)適用性更好。

由于表4和表5給出的是實驗區(qū)域內(nèi)所有格網(wǎng)點(diǎn)的均值,不能反映各格網(wǎng)點(diǎn)上的預(yù)測效果。RMSE能夠更好的反映預(yù)報值的精度及可靠性,常作為預(yù)測結(jié)果評估的主要指標(biāo)[21,30-31]。因此,基于蘭勃特投影,繪制了實驗區(qū)域內(nèi)各格網(wǎng)點(diǎn)上每天的RMSE分布圖。本文僅給出中緯度區(qū)域的分布圖,電離層平靜期與活躍期的分布圖分別如圖6與圖7所示。

圖6 電離層平靜期中緯度區(qū)域內(nèi)格網(wǎng)點(diǎn)上的RMSE分布Fig.6 RMSE distribution at grid nodes of the mid-latitude region during the ionospheric quiet period

圖7 電離層活躍期中緯度區(qū)域內(nèi)格網(wǎng)點(diǎn)上的RMSE分布Fig.7 RMSE distribution at grid nodes of the mid-latitude region during the ionospheric active period

結(jié)合圖6和圖7可直觀地看出,在中緯度地區(qū),本文改進(jìn)模型的RMSE峰值低于未改進(jìn)的Prophet框架,且兩者均顯著優(yōu)于ARIMA模型。特別是在活躍期第3日的預(yù)測中,相比于ARIMA模型,本文改進(jìn)模型能將RMSE峰值降低約4 TECu。此外,平靜期以及活躍期前2日的預(yù)報中,本文改進(jìn)模型的RMSE峰值低于2 TECu,精度較為理想。而在活躍期第三日,RMSE峰值有所增加,約為3 TECu?？偟膩碚f,相比其他模型,本文改進(jìn)模型預(yù)報結(jié)果的RMSE總體較低,精度較為理想,在實驗區(qū)域內(nèi)有良好的適用性,進(jìn)一步驗證了本文改進(jìn)模型的有效性。

進(jìn)一步分析3種模型預(yù)報殘差,并統(tǒng)計了以1 TECu、2 TECu、3 TECu為節(jié)點(diǎn)的不同區(qū)間內(nèi)的殘差比例,在電離層平靜期和活躍期的統(tǒng)計結(jié)果分別如表6和表7所示。

表6 電離層平靜期不同實驗區(qū)域內(nèi)預(yù)報殘差Δ統(tǒng)計表Table 6 Statistical table of prediction residuals Δ in different experimental areas during the ionosphere quiet period

表7 電離層活躍期不同實驗區(qū)域內(nèi)預(yù)報殘差Δ統(tǒng)計表Table 7 Statistical table of prediction residuals Δ in different experimental areas during the ionosphere active period

為更直觀地分析表6與表7的統(tǒng)計結(jié)果,對表6與表7中本文改進(jìn)模型在3天預(yù)測時段內(nèi)的預(yù)報殘差統(tǒng)計結(jié)果進(jìn)行了整合,繪制了百分比環(huán)形圖,并在環(huán)形圖內(nèi)部重點(diǎn)給出了本文改進(jìn)模型在3天預(yù)測時段內(nèi)小于3 TECu的殘差所占的百分比,如圖8所示。

結(jié)合表6、表7及圖8以看出,電離層活躍期預(yù)測時段內(nèi),在低緯度地區(qū),本文改進(jìn)模型3 TECu以內(nèi)的預(yù)報殘差百分比優(yōu)于75%;電離層平靜期預(yù)測時段內(nèi),在低緯度地區(qū)3 TECu以內(nèi)的預(yù)報殘差百分比優(yōu)于85%;而無論是在電離層平靜期還是活躍期,中、高緯度地區(qū)預(yù)測時段內(nèi)的預(yù)報殘差總體上均在3 TECu以內(nèi)。從預(yù)報殘差的角度來看,該結(jié)果與IGS本身提供的TEC值精度相當(dāng)[15]。

圖8 本文改進(jìn)模型3天的總體預(yù)報殘差Δ統(tǒng)計圖Fig.8 Statistical chart of the three-day general prediction residuals Δ calculated by the improved model of this paper

從表6和表7可以看出,3天的預(yù)測時段內(nèi),本文改進(jìn)模型在3 TECu以內(nèi)的預(yù)報殘差百分比略優(yōu)于未改進(jìn)的Prophet框架,而在其他分類區(qū)間內(nèi)互有優(yōu)劣,而上述兩模型的預(yù)報殘差百分比顯著優(yōu)于ARIMA模型。總體來說,本文改進(jìn)模型具有較高的預(yù)報精度,是一種相對簡單且比較理想的預(yù)測方法。

對比表6和表7可以看出,活躍期的預(yù)報殘差大于平靜期,且預(yù)報殘差隨緯度降低而增大。其原因在前文中已有詳細(xì)的分析,簡單來說即電離層受復(fù)雜因素影響,非平穩(wěn)性非線性特征增強(qiáng),導(dǎo)致建模難度增加,偏差增大,嚴(yán)重影響了預(yù)報精度。

3 結(jié) 論

本文提出了一種小波分解與Prophet框架融合的時間序列預(yù)測模型,并在建模中考慮了地磁擾動對電離層的影響,應(yīng)用于電離層VTEC短期預(yù)報。利用IGS發(fā)布的GIM格網(wǎng)數(shù)據(jù),分別基于電離層平靜期數(shù)據(jù)與活躍期數(shù)據(jù)進(jìn)行了實驗,通過RMSE、MAE、MAPE 3個性能評估指標(biāo)對預(yù)報結(jié)果進(jìn)行了評估,并從預(yù)報殘差的角度對預(yù)測精度進(jìn)行了進(jìn)一步的分析。實驗結(jié)果表明該模型的預(yù)報精度優(yōu)于未改進(jìn)的Prophet框架,顯著優(yōu)于以往文獻(xiàn)中常采用的ARIMA模型,且在中、高緯度地區(qū)表現(xiàn)出良好的適用性,驗證了改進(jìn)模型的預(yù)報精度與適用性。

為獲取建模所需的先驗信息,討論了磁暴等特殊事件的確定方法,并對電離層活躍期的地磁擾動強(qiáng)度進(jìn)行了分析。通過研究活躍期時段地磁指數(shù)Ap與DST指數(shù)的變化可知,磁暴發(fā)生時DST指數(shù)達(dá)到極小值,而Ap指數(shù)明顯增大。則可根據(jù)Ap指數(shù)判斷某日的地磁擾動是否劇烈,即可獲取建模所需的先驗信息,并構(gòu)建特殊事件列表。