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

田 睿,董緒榮

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

0 引 言

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

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

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

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

1.1 小波分解基本原理

(1)

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

(2)

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

(3)

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

(4)

1.2 Prophet框架

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

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

(5)

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

1.2.1 趨勢(shì)項(xiàng)模型

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

(6)

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

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

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

(7)

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

(8)

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

(9)

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

(10)

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

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

(11)

應(yīng)注意,增長(zhǎng)率變化向量δ滿足δ～Laplace(0,τ)。其中,參數(shù)τ是本文主要調(diào)整的超參數(shù),因?yàn)槠渲苯涌刂期厔?shì)項(xiàng)增長(zhǎng)率變化的靈活性。增大該參數(shù)會(huì)使趨勢(shì)擬合更加靈活,但存在過擬合風(fēng)險(xiǎn)。

1.2.2 季節(jié)周期性

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

(12)

定義向量X(t)為

X(t)=

(13)

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

s(t)=X(t)β

(14)

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

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

考慮到節(jié)假日或特殊事件對(duì)時(shí)序數(shù)據(jù)的影響,Prophet框架在模型中將節(jié)假日效應(yīng)h(t)作為先驗(yàn)知識(shí)融入到模型中,并認(rèn)為節(jié)假日效應(yīng)對(duì)時(shí)間序列的影響是獨(dú)立的。應(yīng)注意,節(jié)假日或特殊事件需要作為先驗(yàn)信息由用戶給定。設(shè)節(jié)假日或特殊事件i對(duì)應(yīng)的日期列表為Di,如十一黃金周對(duì)應(yīng)的日期列表Di中包含10月1日到10月7日7個(gè)日期。Prophet框架通過一個(gè)指示函數(shù)表示某時(shí)刻t是否處于節(jié)假日或特殊事件i期間,同時(shí)需要一個(gè)參數(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é)假日對(duì)模型的影響越大;該值越小表示節(jié)假日對(duì)模型影響越小,用戶可根據(jù)先驗(yàn)知識(shí)調(diào)整該值。

1.3 分解-集成模型

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

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

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

1.4 小波基的對(duì)比及選取

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

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

表1 常用小波基的各項(xiàng)特性Table 1 Characteristics of commonly used wavelet bases

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)及分析

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

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

表2 實(shí)驗(yàn)區(qū)域Table 2 Experimental area

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

如前文所述,特殊事件列表將作為先驗(yàn)信息輸入模型,并選取地磁指數(shù)Ap作為衡量地磁活動(dòng)指標(biāo)。其中地磁指數(shù)Ap表示行星等效日幅度,可作為全天地磁活動(dòng)水平的度量[45]。為說明磁暴等特殊事件的確定方法,獲取建模所需的先驗(yàn)信息,對(duì)活躍期時(shí)段的地磁擾動(dòng)強(qiáng)度進(jìn)行分析?；钴S期時(shí)段擾動(dòng)暴實(shí)時(shí)(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時(shí)即可能發(fā)生小磁暴,小于-50 nT時(shí)即可能發(fā)生中等磁暴[30]。如圖2所示,年積日47～48日多個(gè)歷元DST指數(shù)較低(可能發(fā)生中小磁暴),與之對(duì)應(yīng)的地磁指數(shù)Ap均達(dá)到或超過5。即可根據(jù)地磁指數(shù)Ap判斷某日的地磁擾動(dòng)是否劇烈,并構(gòu)建特殊事件列表,作為先驗(yàn)信息輸入模型。根據(jù)圖2,預(yù)測(cè)時(shí)段內(nèi)地磁擾動(dòng)并不劇烈,預(yù)測(cè)中無需考慮特殊事件的影響。

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

2.2 小波基選取實(shí)驗(yàn)

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

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

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

2.3 模型精度及適用性分析

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

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

RMSE計(jì)算公式為

(16)

MAE計(jì)算公式為

(17)

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

MAPE計(jì)算公式為

(18)

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

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

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

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

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

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

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

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

由于表4和表5給出的是實(shí)驗(yàn)區(qū)域內(nèi)所有格網(wǎng)點(diǎn)的均值,不能反映各格網(wǎng)點(diǎn)上的預(yù)測(cè)效果。RMSE能夠更好的反映預(yù)報(bào)值的精度及可靠性,常作為預(yù)測(cè)結(jié)果評(píng)估的主要指標(biāo)[21,30-31]。因此,基于蘭勃特投影,繪制了實(shí)驗(yàn)區(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ù)測(cè)中,相比于ARIMA模型,本文改進(jìn)模型能將RMSE峰值降低約4 TECu。此外,平靜期以及活躍期前2日的預(yù)報(bào)中,本文改進(jìn)模型的RMSE峰值低于2 TECu,精度較為理想。而在活躍期第三日,RMSE峰值有所增加,約為3 TECu?？偟膩碚f,相比其他模型,本文改進(jìn)模型預(yù)報(bào)結(jié)果的RMSE總體較低,精度較為理想,在實(shí)驗(yàn)區(qū)域內(nèi)有良好的適用性,進(jìn)一步驗(yàn)證了本文改進(jìn)模型的有效性。

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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