小波-EMD分解改進(jìn)電離層總電子短期預(yù)報(bào)ARIMA模型

陳雨田,劉立龍,黎峻宇,田祥雨

(桂林理工大學(xué) a.測(cè)繪地理信息學(xué)院；b.廣西空間信息與測(cè)繪重點(diǎn)實(shí)驗(yàn)室,廣西 桂林 541006)

0 引 言

電離層是地球大氣受太陽(yáng)輻射和宇宙射線影響而形成的一個(gè)電離區(qū)域。 在衛(wèi)星通訊、 導(dǎo)航和定位等應(yīng)用領(lǐng)域, 由于電磁波穿過(guò)電離層時(shí)發(fā)生折射、 散射、 吸收等現(xiàn)象[1], 產(chǎn)生延遲造成誤差, 而誤差大小與傳播路徑上的總電子含量(total electron content,TEC)成正比[2]。 一些學(xué)者通過(guò)研究地震多發(fā)區(qū)上空的電離層電子濃度, 發(fā)現(xiàn)震前震區(qū)上空的電離層會(huì)產(chǎn)生異常擾動(dòng), 其TEC較平時(shí)發(fā)生明顯增加或減少[3-4]。 因此，對(duì)電離層總電子含量的高精度預(yù)報(bào)既是導(dǎo)航定位的迫切需要, 又是對(duì)地震等自然災(zāi)難分析研究的重要手段。

目前, 常用的電離層總電子預(yù)報(bào)模型主要有經(jīng)典電離層模型[5-6](Bent模型和Klobuchar模型等)和線性、 非線性預(yù)測(cè)模型[7-8](時(shí)間序列和神經(jīng)網(wǎng)絡(luò)等)。 經(jīng)典電離層模型在本領(lǐng)域經(jīng)過(guò)較長(zhǎng)時(shí)間的發(fā)展和改進(jìn)具有易于操作、 計(jì)算簡(jiǎn)單等優(yōu)點(diǎn), 但其固有的缺陷并沒(méi)有得到很好的彌補(bǔ), 模型預(yù)報(bào)精度較低, 一般情況下預(yù)報(bào)值精度僅約50%～60%, 且受地域影響較大[9-10]。 自回歸積分滑動(dòng)平均模型(autoregressive integrated moving average model,ARIMA)是時(shí)間序列中預(yù)測(cè)精度較高、 實(shí)際應(yīng)用廣泛的一種模型, 能夠?qū)Ψ瞧椒€(wěn)、 非線性的時(shí)間序列進(jìn)行有效預(yù)報(bào)[11-12], 但電離層TEC的時(shí)空分布差異較大, 單一ARIMA模型的預(yù)報(bào)精度較低。 國(guó)內(nèi)一些學(xué)者通過(guò)利用小波分解或經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)對(duì)電離層的TEC變化進(jìn)行探索，并改進(jìn)了ARIMA預(yù)報(bào)模型[13-14], 獲得了一些有益的結(jié)論, 然而改進(jìn)模型的精度仍有提升的空間。 本文在過(guò)往研究的基礎(chǔ)上進(jìn)一步探索, 將小波分解與EMD方法相結(jié)合, 共同納入ARIMA模型,構(gòu)建WEARIMA預(yù)報(bào)模型, 并采用IGS提供的GIM數(shù)據(jù)驗(yàn)證改進(jìn)模型的預(yù)報(bào)精度。

1 基本原理與模型改進(jìn)

1.1 小波分析基本原理

小波分析是一種對(duì)信號(hào)的局部化分析,通過(guò)小波分析可以將信號(hào)序列分解為一個(gè)低頻序列和多個(gè)高頻序列。低頻序列表現(xiàn)為高頻率分辨率和低時(shí)間分辨率,高頻序列表現(xiàn)為低頻率分辨率和高時(shí)間分辨率。在進(jìn)行小波分析時(shí),定義平方可積函數(shù)ψ(t)為基本小波函數(shù),滿足條件:

(1)

將基本小波函數(shù)ψ(t)伸縮和平移后即可得到一個(gè)小波序列:

(2)

式中:a和b分別為伸縮因子和平移因子,a、b∈R，a≠0， 則任意函數(shù)x(t)的連續(xù)小波變換及其逆變換可表示為

(3)

(4)

1.2 EMD方法簡(jiǎn)介

經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)是一種能夠有效處理非線性、非平穩(wěn)信號(hào)的新型信號(hào)處理方法。根據(jù)信號(hào)的特征,通過(guò)EMD分解可以將復(fù)雜信號(hào)序列x(t)從高頻到低頻進(jìn)行多尺度分解, 從而得到頻率由高到低排列的若干個(gè)IMF分量和殘余函數(shù)R(n), 其實(shí)質(zhì)是將信號(hào)序列由高頻至低頻進(jìn)行平穩(wěn)化處理[15],即

(5)

EMD分解要求分解得到的每個(gè)IMF分量滿足兩個(gè)條件[16]:第一,在整個(gè)信號(hào)序列中極值點(diǎn)與零點(diǎn)數(shù)目相同或最多相差一個(gè);第二,由極大值和極小值構(gòu)成的上下包絡(luò)關(guān)于時(shí)間軸局部對(duì)稱。由此可以看出,對(duì)于幅度與頻率相對(duì)較為對(duì)稱的信號(hào),EMD分解會(huì)取得較好的效果,而對(duì)時(shí)域波形發(fā)生畸變的信號(hào)分解效果則較差。

1.3 小波-EMD改進(jìn)ARIMA模型(WEARIMA)

自回歸積分滑動(dòng)平均模型(ARIMA)是對(duì)自回歸滑動(dòng)平均模型(ARMA)的擴(kuò)展。ARMA(p,q)模型的一般形式為

xt=φ1xt-1+φ2xt-2+…+φpxt-p+εt+θ1εt-1+θ2εt-2+…+θqεt-q，

(6)

式中：εt為白噪聲序列, 并且與t時(shí)刻前的原始序列x(k)(k

ARIMA(p,d,q)×(P,D,Q)考慮了季節(jié)因素影響模型[11], 具有建模簡(jiǎn)單、 操作方便等特點(diǎn),在對(duì)具有季節(jié)性的非平穩(wěn)時(shí)間序列的預(yù)報(bào)中具有廣泛的應(yīng)用。 模型中的p、d、q分別為自回歸階數(shù)、 差分階數(shù)和移動(dòng)平均數(shù); 而參數(shù)P為季節(jié)性自回歸階數(shù),D為常規(guī)差分階數(shù),Q為季節(jié)性移動(dòng)平均數(shù)。

WEARIMA模型的建模流程如圖1所示,建模步驟如下:

① 提取IGS中心發(fā)布的30天電離層TEC格網(wǎng)數(shù)據(jù),將前25天作為樣本序列,利用db4正交小波對(duì)樣本序列進(jìn)行一級(jí)分解(用以提高序列的平穩(wěn)性及相對(duì)于時(shí)間軸的對(duì)稱性),對(duì)分解出的參數(shù)選取合適的閾值進(jìn)行分析,通過(guò)逆小波變換對(duì)分解結(jié)果重構(gòu),得到高頻序列d1和低頻序列a1;

② 將小波分解得到的低頻序列a1進(jìn)行EMD分解,得到若干個(gè)IMF分量和單調(diào)殘余函數(shù)R(n)(EMD分解后的趨勢(shì)余量);

③ 將高頻序列d1和EMD分解后的各個(gè)分量分別代入ARIMA模型中,根據(jù)不同分量表現(xiàn)出的自相關(guān)與偏相關(guān)特性,確定各模型的參數(shù);

④ 對(duì)各分量的預(yù)報(bào)值進(jìn)行求和,得到最終的預(yù)報(bào)結(jié)果。

圖1 小波-EMD改進(jìn)ARIMA模型(WEARIMA)建模流程Fig.1 Modeling process of improved ARIMA model based on wavelet-EMD

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

2.1 實(shí)驗(yàn)數(shù)據(jù)準(zhǔn)備

由EMD分解后IMF分量的特性可知，當(dāng)原始信號(hào)序列的幅度與頻率較為對(duì)稱時(shí)，分解效果較好,利用ARIMA模型預(yù)報(bào)的精度也較高。太陽(yáng)輻射對(duì)大氣分子的作用是形成電離層的主要因素[17],在太陽(yáng)活動(dòng)的活躍期往往會(huì)伴隨磁暴等現(xiàn)象,造成電離層總電子含量(TEC)異常擾動(dòng),不利于EMD分解后的預(yù)報(bào),且使得預(yù)報(bào)結(jié)果具有偶然性,故選取太陽(yáng)活動(dòng)較為平靜的2009年年積日122～151共30天的GIM電離層格網(wǎng)數(shù)據(jù)進(jìn)行實(shí)驗(yàn),利用前25天的TEC觀測(cè)值作為樣本序列，預(yù)報(bào)后5天的TEC值。

圖2給出了預(yù)報(bào)期內(nèi)表征地磁活動(dòng)強(qiáng)度的Kp指數(shù)與表征全球環(huán)電流強(qiáng)度的Dst指數(shù)。Kp指數(shù)可以反映3 h的地磁活動(dòng)強(qiáng)度,一般認(rèn)為，Kp值小于5時(shí)不會(huì)發(fā)生不同等級(jí)的地磁風(fēng)暴；每日Kp指數(shù)之和小于30表示地磁活動(dòng)平靜。而Dst指數(shù)的時(shí)間分辨率為1 h,當(dāng)Dst指數(shù)不低于-30 nT時(shí)不會(huì)發(fā)生任何等級(jí)的磁暴。預(yù)報(bào)期的地磁指數(shù)均滿足上述條件,故可排除磁暴與地磁活動(dòng)對(duì)預(yù)報(bào)結(jié)果的影響。

2.1 實(shí)驗(yàn)與精度分析

圖3是對(duì)年積日122～146的數(shù)據(jù)在格網(wǎng)點(diǎn)(45°N, 110°E)處的小波-EMD分解, 采用db4正交小波對(duì)原始序列進(jìn)行一級(jí)分解, 得到低頻序列a1與高頻序列d1, 高頻序列d1的極值介于-2～2 TECu, 使用ARIMA模型直接進(jìn)行預(yù)報(bào)。 低頻序列a1極值介于5～15 TECu, 利用EMD方法對(duì)其進(jìn)一步分解得到5個(gè)IMF分量和1個(gè)趨勢(shì)余量R(6), 將得到的各分量分別使用ARIMA模型對(duì)后5天的TEC進(jìn)行預(yù)報(bào),結(jié)果如圖4所示。除IMF1與IMF2以外,其余各分量的預(yù)報(bào)值與真實(shí)值之間相差很小,其中IMF5與趨勢(shì)余量R(6)和真實(shí)值基本一致。

實(shí)驗(yàn)中為了方便表示,將僅用小波分解改進(jìn)的ARIMA模型記為WARIMA模型,小波-EMD改進(jìn)的ARIMA模型記為WEARIMA模型,通過(guò)ARIMA、WARIMA和WEARIMA分別對(duì)上述格網(wǎng)點(diǎn)(45°N, 110°E)進(jìn)行5天預(yù)報(bào),得到各模型預(yù)報(bào)結(jié)果與IGS觀測(cè)值的對(duì)比情況,如圖5所示?？梢钥闯?僅采用ARIMA模型得到的預(yù)報(bào)結(jié)果在每日TEC的極值附近與觀測(cè)值有較大偏差,而WARIMA與WEARIMA相對(duì)于原模型能夠有效削弱極值點(diǎn)附近預(yù)報(bào)值的誤差,在預(yù)報(bào)性能上有較大的提升,且WEARIMA模型的改進(jìn)效果更好。對(duì)預(yù)報(bào)值殘差的絕對(duì)值進(jìn)行統(tǒng)計(jì)分析,得到表1。

從后5天的總體預(yù)報(bào)來(lái)看，ARIMA模型在改進(jìn)前僅有60%左右的殘差絕對(duì)值小于1 TECu,通過(guò)兩種方法改進(jìn)后誤差絕對(duì)值小于1 TECu的可增大約10%;改進(jìn)后誤差絕對(duì)值大于3 TECu的較原模型明顯減少,且WEARIMA模型的預(yù)報(bào)值在3天中沒(méi)有出現(xiàn)大于3 TECu的誤差。

圖2 預(yù)報(bào)期Dst與Kp指數(shù)(年積日122～151)Fig.2 Dst and Kp index during forecast period(DOY 122-151)

圖3 樣本序列小波-EMD分解圖(年積日122～146)Fig.3 Wavelet-EMD decomposition of sample sequence(DOY 122-146)

僅通過(guò)單個(gè)格網(wǎng)點(diǎn)數(shù)據(jù)與預(yù)報(bào)值殘差的絕對(duì)值比例并不能很好地體現(xiàn)出改進(jìn)模型的優(yōu)勢(shì),為進(jìn)一步說(shuō)明小波-EMD改進(jìn)模型在電離層TEC預(yù)報(bào)上較原模型與WARIMA模型有更高的預(yù)報(bào)精度,本文將IGS中心發(fā)布的觀測(cè)值視為真實(shí)值,采用平均相對(duì)精度(P)和均方根誤差(RMSE)進(jìn)一步評(píng)定3種模型預(yù)報(bào)結(jié)果的優(yōu)劣。P、RMSE指標(biāo)為

(7)

(8)

式中:IPRE和IIGS分別為模型的預(yù)報(bào)值與IGS觀測(cè)值;m和n分別為當(dāng)天觀測(cè)的起始?xì)v元和結(jié)束歷元。

由于電離層TEC值受太陽(yáng)直射影響較大,其周日變化幅度隨地理經(jīng)度的變化較小，而隨緯度的升高呈現(xiàn)明顯的降低趨勢(shì),表2給出了北半球低緯(15°N, 110°E)、 中緯(45°N, 110°E)和高緯(75°N, 110°E)3個(gè)緯度地區(qū)ARIMA、WARIMA和WEARIMA模型預(yù)報(bào)值的日平均相對(duì)精度。從緯度變化上看,3種模型預(yù)報(bào)精度隨緯度的升高而增加。 在同一格網(wǎng)點(diǎn)上, 除低緯度WEARIMA模型的相對(duì)精度低于WARIMA以外,其余緯度的平均相對(duì)精度均表現(xiàn)為WEARIMA模型高于WARIMA模型,高于ARIMA模型。各模型5 d預(yù)報(bào)值均方根誤差的計(jì)算結(jié)果見(jiàn)表3,均方根誤差隨緯度的升高而減小,在同一格網(wǎng)點(diǎn)上WEARIMA模型預(yù)報(bào)值的均方根誤差要明顯低于WARIMA和ARIMA模型,具有最為可靠的預(yù)報(bào)性能。

圖4 分解后各分量的5 d預(yù)報(bào)值(點(diǎn)線)與觀測(cè)值(實(shí)線)對(duì)比Fig.4 Comparison between the 5 d prediction value and the observed value of each component obtained after decomposition

圖5 各模型預(yù)報(bào)結(jié)果對(duì)比Fig.5 Comparison of prediction results of each model

圖6為中國(guó)及周邊地區(qū)上空TEC預(yù)報(bào)值的均方根誤差(RMSE)分布情況。 從空間上看,隨緯度的降低，RMSE呈現(xiàn)明顯的增大趨勢(shì), 且在中國(guó)西南部地區(qū)達(dá)到最大值, 即表現(xiàn)為“駝峰”; 從時(shí)間上分析, 除第一天的預(yù)報(bào)結(jié)果在中國(guó)南部地區(qū)有較大誤差, 其余預(yù)報(bào)值隨時(shí)間推移誤差逐漸增大, 在最后一天達(dá)到極大值。 通過(guò)對(duì)比每天的RMSE值以及5天的RMSE平均值, 可以發(fā)現(xiàn)WEARIMA模型在該區(qū)域的預(yù)報(bào)精度相較于ARIMA模型有明顯提升, 在預(yù)報(bào)誤差較大的中國(guó)西南部地區(qū), WEARIMA模型能夠有效削弱RMSE峰值,提高預(yù)報(bào)精度。

表1 5 d預(yù)報(bào)值殘差絕對(duì)值(Δ)統(tǒng)計(jì)

Table 1 Absolute value of the residuals(Δ)of the predicted values over the next 5 days %

預(yù)報(bào)時(shí)間/dARIMA模型/WARIMA模型/WEARIMA模型Δ<1TECu1TECu≤Δ<2TECu2TECu≤Δ<3TECuΔ≥3TECu163.89/66.67/66.6719.44/13.89/19.448.33/11.11/5.568.33/8.33/8.33275.00/86.11/69.4419.44/8.33/19.440.00/0.00/11.115.56/5.56/0.00358.33/72.22/75.0027.78/13.89/13.895.56/8.33/8.338.33/5.56/2.78461.11/72.22/69.4422.22/13.89/16.6755.56/5.56/13.8911.11/8.33/0.00552.78/66.67/77.7819.44/16.67/16.6711.11/8.33/5.5616.67/8.33/0.00

表2 不同模型5 d預(yù)報(bào)值相對(duì)精度對(duì)比

Table 2 Comparison of relative accuracy of predicted values of different models over the next 5 days %

經(jīng)緯度模 型預(yù)報(bào)天數(shù)/d12345平均相對(duì)精度ARIMA82.9791.0084.9485.8078.46846315°N,110°EWARIMA87.9693.4990.3789.6185.2189.33WEARIMA80.3083.5989.9379.4789.9984.66ARIMA91.3992.1590.7789.3689.2690.5945°N,110°EWARIMA88.8492.3391.8793.0094.0292.01WEARIMA95.0394.9193.3195.1195.4094.75ARIMA95.0394.1391.3393.5295.2193.8475°N,110°EWARIMA95.5594.8692.1393.0995.1994.16WEARIMA95.5793.8493.3094.6996.3194.74

表3 不同模型5 d預(yù)報(bào)值均方根誤差(RMSE)

圖6 RMSE空間分布情況對(duì)比Fig.6 Comparison of RMSE spatial distribution

3 結(jié) 論

采用小波與EMD分解相結(jié)合的方法改進(jìn)ARIMA預(yù)報(bào)模型,根據(jù)IGS提供的GIM電離層格網(wǎng)數(shù)據(jù),分析改進(jìn)模型在中國(guó)及周邊地區(qū)的預(yù)報(bào)精度,并與單一小波分解改進(jìn)的模型及原模型進(jìn)行對(duì)比分析,得到如下結(jié)論:

(1)對(duì)小波分解后的低頻序列進(jìn)行EMD分解可得到若干個(gè)IMF分量與一個(gè)趨勢(shì)余量,對(duì)其分別預(yù)報(bào)得到的預(yù)報(bào)結(jié)果表現(xiàn)為隨頻率降低與真值的擬合度升高,并且較單一小波分解具有一定的提升,實(shí)驗(yàn)表明該方法在提高電離層TEC預(yù)報(bào)精度上有一定作用。

(2)分析各模型預(yù)報(bào)值的相對(duì)精度與均方根誤差發(fā)現(xiàn),總體而言，3種模型的預(yù)報(bào)精度隨緯度的升高而增加,改進(jìn)模型的預(yù)報(bào)精度較原模型有較大提高,3種模型預(yù)報(bào)精度表現(xiàn)為WEARIMA模型>WARIMA模型>ARIMA模型。

(3)對(duì)中國(guó)及周邊地區(qū)上空進(jìn)行預(yù)報(bào)并分析區(qū)域均方根誤差分布情況,得到WEARIMA模型能有效削弱中國(guó)西南部地區(qū)預(yù)報(bào)值均方根誤差的峰值,與原模型相比精度有較大的提升。

本文基于小波-EMD分解改進(jìn)ARIMA模型,為電離層總電子預(yù)報(bào)提供了新的思路。但是,電離層各項(xiàng)參數(shù)受太陽(yáng)活動(dòng)及其他因素的影響較大,如何顧及這些因素進(jìn)行預(yù)報(bào)以及如何選取小波基、確定ARIMA模型階數(shù)等問(wèn)題仍有待進(jìn)一步研究。