基于GNSS TEC 的電離層自相關(guān)預(yù)報(bào)誤差分析及參數(shù)優(yōu)化

熊雯，王博文，劉裔文,3，朱慶林

(1.中國(guó)電波傳播研究所,山東 青島 266107；2.南昌大學(xué) 信息工程學(xué)院,南昌 330031；3.上饒師范學(xué)院 物理與電子信息學(xué)院,江西 上饒 334001)

0 引言

電離層電子密度的變化是影響無(wú)線電系統(tǒng)工作性能的重要誤差源之一，對(duì)電離層總電子含量(TEC)進(jìn)行預(yù)報(bào)是減緩電離層誤差的主要途徑.尤其對(duì)全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS) 應(yīng)用具有重要使用價(jià)值.例如，對(duì)GNSS 單頻用戶(hù)，TEC 預(yù)報(bào)是有效減緩電離層延遲誤差影響的重要途徑.另外，電離層預(yù)報(bào)也是提高基于GNSS 信標(biāo)觀測(cè)的全球電離層地圖(GIM)精度的重要手段.

從預(yù)報(bào)時(shí)間提前量的角度看，電離層預(yù)報(bào)模型可分為長(zhǎng)期預(yù)報(bào)模型和短期預(yù)報(bào)模型.長(zhǎng)期預(yù)報(bào)模型一般是太陽(yáng)活動(dòng)指數(shù)或地磁活動(dòng)指數(shù)驅(qū)動(dòng)，能反應(yīng)電離層周日、季節(jié)、11 年周期變化等大尺度變化的全球性模型，例如國(guó)際參考電離層模型、NeQuick 模型等.短期預(yù)報(bào)模型通常基于臨近的歷史觀測(cè)數(shù)據(jù)的建模來(lái)預(yù)測(cè)當(dāng)?shù)匚磥?lái)1～24 h 的結(jié)果.因短期預(yù)報(bào)具有較高的區(qū)域預(yù)測(cè)精度，且效率高，故常被用于一些對(duì)時(shí)效性和精度要求較高的導(dǎo)航、通信、雷達(dá)等系統(tǒng)的電離層延遲修正.

早在2000 年左右，就有學(xué)者提出了一些有效的短期預(yù)報(bào)方法，例如自相關(guān)方法、自回歸方法和人工神經(jīng)網(wǎng)絡(luò)等[1-3].目前，工程上用的較多的是自相關(guān)或自回歸方法[4-8].制約電離層短期預(yù)報(bào)精度的主要因素主要是磁暴、亞暴等引起的電離層擾動(dòng)和平靜期電離層存在的一些局地小尺度擾動(dòng).這些擾動(dòng)的變化很復(fù)雜，至今依然是研究的熱點(diǎn)，目前還很難對(duì)其進(jìn)行準(zhǔn)確預(yù)報(bào).總的來(lái)說(shuō)，對(duì)于電離層短期預(yù)報(bào)，至今并沒(méi)有一種絕對(duì)優(yōu)勢(shì)的方法，尋求對(duì)已有傳統(tǒng)方法進(jìn)行優(yōu)化仍然是提高預(yù)測(cè)精度的重要途徑之一[9-11].

自相關(guān)預(yù)報(bào)是一種模型簡(jiǎn)單且精度較高的預(yù)報(bào)方法，本文針對(duì)自相關(guān)預(yù)報(bào)方法中的參數(shù)設(shè)置對(duì)誤差的影響進(jìn)行深入分析，進(jìn)而提出優(yōu)化方案.首先介紹所使用的數(shù)據(jù)和處理方法，其次介紹自相關(guān)預(yù)報(bào)原理和方法，在再次比較自相關(guān)與自回歸方法的基礎(chǔ)上，開(kāi)展誤差分析和參數(shù)優(yōu)化分析；最后進(jìn)行討論和總結(jié).

1 數(shù) 據(jù)

本文研究所采用的數(shù)據(jù)主要來(lái)自麻省理工大學(xué)CEDAR Madrigal 數(shù)據(jù)庫(kù)，其數(shù)據(jù)來(lái)自于分布全球的地基衛(wèi)星接收機(jī)網(wǎng)絡(luò)(Distributed Ground Based Satellite Receivers)中的GNSS 接收機(jī)網(wǎng)絡(luò)(World-wide GNSS Receiver Network) (http://cedar.openmadrigal.org/list).

Madrigal 收集的TEC 測(cè)量值均為垂直TEC，單位為T(mén)ECU，1 TECU 表示每平方米有1016個(gè)電子.垂直TEC 是由斜TEC 通過(guò)電離層薄層假設(shè)(高度為350 km)轉(zhuǎn)換得到.最后通過(guò)格網(wǎng)化高精度TEC 地圖重構(gòu)算法生成并提供1°×1°×5 min 時(shí)空分辨率的全球垂直TEC 地圖[12].

但是，Madrigal 數(shù)據(jù)庫(kù)并沒(méi)有直接提供單站垂直TEC 數(shù)據(jù).為開(kāi)展單站電離層預(yù)報(bào)方法研究，本文在全球TEC 地圖中從觀測(cè)站網(wǎng)密集區(qū)域任意選取了5 個(gè)格網(wǎng)點(diǎn)的TEC 數(shù)據(jù)，然后通過(guò)1 h 平均得到單站垂直TEC 日變化時(shí)間序列，以此作為本文的實(shí)測(cè)TEC 時(shí)間序列.圖1 給出了所選取的5 個(gè)網(wǎng)格點(diǎn)分布圖.圖1 中5 個(gè)網(wǎng)格點(diǎn)的坐標(biāo)分別為：1(45°N，120°W)；2(30°N，105°W) ；3(40°N，90°W) ；4(30°S，25°E)；5(28°N，85°E).

圖1 選取的5 個(gè)網(wǎng)格點(diǎn)坐標(biāo)示意圖

2 自相關(guān)預(yù)報(bào)原理和方法

在1 個(gè)月左右的時(shí)間尺度上，電離層參量的時(shí)間序列通常是平穩(wěn)時(shí)間序列，在統(tǒng)計(jì)學(xué)上可將其視為具有非零自相關(guān)函數(shù)的隨機(jī)穩(wěn)態(tài)過(guò)程的實(shí)現(xiàn)[1].假設(shè)某電離層參量小時(shí)值的時(shí)間序列為z(t)，其在t時(shí)刻的值可表示為前n個(gè)測(cè)量值的加權(quán)平均，如式(1)所示：

式(1)中的加權(quán)系數(shù) λj通過(guò)以下方程組求解得到：

式中：μ 為拉格朗日系數(shù)；ρ(τ) 為z(t) 的自相關(guān)系數(shù)，其均值為0.對(duì)于1 h 分辨率的電離層時(shí)間序列，引入式(3)所示的擬合函數(shù)來(lái)進(jìn)行加權(quán)系數(shù)的估計(jì)[1].

上述處理是從信號(hào)的角度將電離層參量的時(shí)間序列看成是一個(gè)周期信號(hào)疊加一個(gè)非周期性的擾動(dòng).對(duì)于周期性信號(hào)，我們熟知電離層受日照影響存在明顯的周日變化，根據(jù)周期信號(hào)的傅里葉級(jí)數(shù)理論，將自相關(guān)系數(shù)的周期成分用周期為24 h 的基波和2 次諧波的和來(lái)表示.而對(duì)于非周期的擾動(dòng)信號(hào)，式(3)實(shí)際采用的是一個(gè)幅度受限的雙邊指數(shù)信號(hào)來(lái)模擬其變化，其中 τ0是時(shí)間常數(shù).式(3)中的三個(gè)系數(shù)A24、A12、τ0通過(guò)最小二乘法來(lái)確定，然后，選取最近4 天的12 個(gè)最大的自相關(guān)系數(shù)代入式(2)求解，即可得到對(duì)應(yīng)所選時(shí)間的權(quán)重系數(shù)λj[1].最后，將權(quán)重系數(shù)代入式(1)即可得到t時(shí)刻的預(yù)測(cè)值.

3 結(jié)果

3.1 自相關(guān)與自回歸預(yù)報(bào)方法比較

電離層TEC 觀測(cè)序列是一種典型的周期性時(shí)間序列，根據(jù)時(shí)間序列預(yù)測(cè)理論，對(duì)于數(shù)據(jù)量較少的情況，通常采用自相關(guān)或自回歸模型進(jìn)行預(yù)報(bào).本文首先對(duì)兩類(lèi)模型進(jìn)行了比較，自回歸模型選用了常用的ARIMA 模型.圖2 給出了一個(gè)結(jié)果示例，所用數(shù)據(jù)為網(wǎng)格點(diǎn)5 的數(shù)據(jù)，預(yù)報(bào)時(shí)間為2020 年2 月1 日至2 月10 日.從圖2 可以看出，自相關(guān)和ARIMA 預(yù)報(bào)結(jié)果均與觀測(cè)結(jié)果吻合較好，很好地再現(xiàn)了電離層TEC 的周日變化特征.

圖2 自相關(guān)與自回歸預(yù)報(bào)方法結(jié)果比較示例

進(jìn)一步地，我們采用均方根偏差(RMSE)和歸一化均方根偏差(NRMSE)來(lái)衡量預(yù)報(bào)絕對(duì)誤差和相對(duì)誤差，計(jì)算公式為：

上兩式中：εRMSE和 εNRMSE分別為均方根偏差和歸一化均方根偏差；和zi分別為第i個(gè)點(diǎn)位的預(yù)測(cè)值和觀測(cè)值；N為總的點(diǎn)位數(shù).

經(jīng)估算，圖2 所示結(jié)果中，自相關(guān)預(yù)報(bào)的絕對(duì)誤差約為0.716 5 TECU，相對(duì)誤差約為10.87%；ARIMA預(yù)報(bào)絕對(duì)誤差和相對(duì)誤差分別為0.876 2 TECU 和13.29%.此外，我們還利用網(wǎng)格點(diǎn)1～4 的數(shù)據(jù)對(duì)兩種方法進(jìn)行了比較分析，計(jì)算結(jié)果顯示，對(duì)于網(wǎng)格點(diǎn)1～4，自相關(guān)方法的絕對(duì)誤差分別為0.414 5 TECU、0.654 3 TECU、0.515 2 TECU、1.182 1 TECU，ARIMA方法絕對(duì)誤差分別為0.436 9 TECU、0.607 8 TECU、0.522 6 TECU、1.230 1 TECU；自相關(guān)方法的相對(duì)誤差分別為8.3%、9.52%、8.52%、13.28%，ARIMA 方法相對(duì)誤差分別為8.75%、8.84%、8.64%、13.82%.

從上述結(jié)果對(duì)比可以看出，總體上自相關(guān)方法的絕對(duì)誤差和相對(duì)誤差均略小于ARIMA 方法.此外，我們還對(duì)兩種預(yù)報(bào)方法需要花費(fèi)的時(shí)間進(jìn)行了測(cè)算.結(jié)果顯示，預(yù)報(bào)如圖2 所示的10 天結(jié)果，在相同的計(jì)算環(huán)境下，自相關(guān)預(yù)報(bào)時(shí)間約為10 s，ARIMA 預(yù)報(bào)時(shí)間約為50 s.還需指出的是，本文采用的ARIMA 模型使用之前已經(jīng)定好了階數(shù)，后續(xù)的預(yù)測(cè)均使用固定階數(shù)，若采用動(dòng)態(tài)自動(dòng)定階，則預(yù)報(bào)時(shí)間將大大增加.

綜上分析，自相關(guān)預(yù)報(bào)方法預(yù)報(bào)誤差相對(duì)略?xún)?yōu)于ARIMA 方法，且自相關(guān)預(yù)報(bào)所需時(shí)間比ARIMA少得多.因此，本文后續(xù)將主要分析自相關(guān)方法.如前所述，自相關(guān)方法雖然已經(jīng)取得了一些應(yīng)用，但其依然有改進(jìn)的空間，3.2 小節(jié)將在誤差分析的基礎(chǔ)上嘗試對(duì)自相關(guān)方法進(jìn)行參數(shù)優(yōu)化.

3.2 自相關(guān)誤差分析和參數(shù)優(yōu)化

對(duì)于自相關(guān)預(yù)報(bào)方法，在預(yù)報(bào)參數(shù)選取方面，目前普遍的做法是，選取最近4 天相關(guān)系數(shù)最高的12 個(gè)時(shí)間點(diǎn)，然后通過(guò)式(2)解出相應(yīng)的權(quán)重系數(shù).由于這兩個(gè)參數(shù)是早期經(jīng)驗(yàn)值，其對(duì)預(yù)報(bào)誤差的影響目前還沒(méi)有文獻(xiàn)公開(kāi)報(bào)道.為弄清楚這一影響，我們對(duì)具體預(yù)報(bào)過(guò)程和誤差產(chǎn)生原因進(jìn)行了深入分析.

圖3 給出了參數(shù)選取4 天和12 點(diǎn)時(shí)的一個(gè)預(yù)報(bào)結(jié)果示例.圖3 中的TEC 數(shù)據(jù)來(lái)自網(wǎng)格點(diǎn)4 (30°S，25°E)，預(yù)測(cè)日期為2020 年2 月8 日，所用的歷史數(shù)據(jù)為2020 年1 月9 日至2 月7 日共30 天數(shù)據(jù)，分辨率為1 h.圖3 中黑色圓點(diǎn)表示實(shí)測(cè)值，2 月8 日的紅色圓點(diǎn)表示預(yù)測(cè)值.我們以第8 個(gè)預(yù)測(cè)點(diǎn)(如圖中黑色箭頭所指)為例，給出了用來(lái)計(jì)算該點(diǎn)預(yù)測(cè)值的歷史數(shù)據(jù)分布(如圖中藍(lán)色圓圈所示)，同時(shí)在每個(gè)歷史點(diǎn)位左側(cè)標(biāo)出了權(quán)重系數(shù).根據(jù)自相關(guān)預(yù)測(cè)理論，所選取的12 個(gè)采樣點(diǎn)是前4 天相關(guān)系數(shù)最大的采樣點(diǎn)，這些點(diǎn)主要是前4 天相同時(shí)刻和相鄰時(shí)刻的采樣點(diǎn)，且相同時(shí)刻的權(quán)重最大(為0.093)，而相鄰時(shí)刻采樣點(diǎn)的權(quán)重均為0.079.

圖3 2020 年2 月8 日“4+12”參數(shù)方案時(shí)的預(yù)報(bào)結(jié)果(紅色圓點(diǎn))以及2 月4 至8 日實(shí)測(cè)結(jié)果(黑色圓點(diǎn))

對(duì)比圖3 中2 月8 日的預(yù)報(bào)結(jié)果和實(shí)測(cè)結(jié)果可以看出，雖然兩條曲線的變化趨勢(shì)基本一致，但幾乎所有點(diǎn)位的預(yù)測(cè)值均要大于觀測(cè)值.對(duì)黑色箭頭所示第8 個(gè)點(diǎn)位，其觀測(cè)值為9.264 TECU，預(yù)測(cè)值為10.292 TECU，誤差達(dá)到1.029 TECU.從第8 個(gè)預(yù)測(cè)點(diǎn)所用歷史點(diǎn)位信息來(lái)看，距離預(yù)測(cè)點(diǎn)最遠(yuǎn)的2 月4 日的兩個(gè)被選用的點(diǎn)位對(duì)誤差影響較大，這兩個(gè)采樣點(diǎn)數(shù)值明顯高于2 月5 日至2 月8 日其他點(diǎn)位的值，很可能是造成該點(diǎn)預(yù)測(cè)值明顯偏離實(shí)測(cè)值的主要原因.

圖3 結(jié)果反映出距離預(yù)測(cè)點(diǎn)較遠(yuǎn)的歷史數(shù)據(jù)可能與當(dāng)前時(shí)刻觀測(cè)值的偏差較大，從而影響預(yù)報(bào)精度.據(jù)此，我們推測(cè)，縮短參與預(yù)報(bào)的歷史數(shù)據(jù)天數(shù)，有可能能夠減少預(yù)報(bào)誤差.因此，我們將參與預(yù)測(cè)的天數(shù)減少為3 天，相應(yīng)的點(diǎn)數(shù)減為9 個(gè)，對(duì)上述選例重新開(kāi)展自相關(guān)預(yù)測(cè)，結(jié)果如圖4 所示.由圖4 可知，整體上，修改參數(shù)后的預(yù)測(cè)值與觀測(cè)值更接近.對(duì)黑色箭頭所示第8 個(gè)點(diǎn)位，其觀測(cè)值為9.264 TECU，預(yù)測(cè)值為9.863 TECU，誤差減小為0.599 TECU.

圖4 2020 年2 月8 日“3+9”參數(shù)方案時(shí)的預(yù)報(bào)結(jié)果(紅色圓點(diǎn))以及2 月4 至8 日實(shí)測(cè)結(jié)果(黑色圓點(diǎn))

此外，我們還分析了天數(shù)增加或減小時(shí)的結(jié)果，圖5 給出了不同參數(shù)設(shè)置下的所有24 個(gè)點(diǎn)位預(yù)測(cè)誤差曲線比較，圖5 中黑色圓圈代表參數(shù)設(shè)置為“4+12”時(shí)的誤差曲線(對(duì)應(yīng)圖3 所示結(jié)果)，紅色圓圈代表參數(shù)設(shè)置為“3+9”時(shí)的誤差曲線(對(duì)應(yīng)圖4 所示結(jié)果).對(duì)比上述兩條誤差曲線可以看出，除12 UT、13 UT、14 UT 三個(gè)點(diǎn)位外，其余點(diǎn)位的預(yù)測(cè)誤差都有所減小，尤其是10 UT 之前的點(diǎn)位減少更明顯.總體上，原方法“4+12”平均絕對(duì)預(yù)報(bào)誤差為0.8 TECU，參數(shù)調(diào)整為“3+9”后平均絕對(duì)預(yù)報(bào)誤差為0.63 TECU，平均絕對(duì)誤差減小了0.17 TECU.此外，圖5 還給出了參數(shù)設(shè)置分別為“5+12”(紫色圓圈)和“2+6”(綠色圓圈)時(shí)的預(yù)報(bào)誤差曲線，但是從對(duì)比結(jié)果看，兩種方案均存在誤差明顯增加的點(diǎn)位，而誤差減少的點(diǎn)位較少.經(jīng)估計(jì)，圖5 中“5+12”案的平均絕對(duì)誤差為0.885 TECU，“2+6”方案的平均絕對(duì)誤差為0.892 TECU，相對(duì)于原方法“4+12”誤差還分別增加了0.085 TECU和0.092 TECU.

圖5 2020 年2 月8 日不同參數(shù)方案預(yù)報(bào)誤差比較

上述結(jié)果表明，在TEC 自相關(guān)短期預(yù)報(bào)方法中，參與預(yù)測(cè)的天數(shù)設(shè)置為3 天(相應(yīng)的點(diǎn)位設(shè)置為9)可能是更優(yōu)的方案.這從側(cè)面說(shuō)明電離層當(dāng)前的狀態(tài)可能主要取決于前3 天的狀態(tài)，或者說(shuō)前3 天的電離層狀態(tài)對(duì)當(dāng)前狀態(tài)的關(guān)聯(lián)最大.時(shí)間跨度越大，關(guān)聯(lián)性越小，引起誤差的可能性越大.

為進(jìn)一步檢驗(yàn)“3+9”參數(shù)方案的預(yù)報(bào)精度改善效果，我們對(duì)圖1 中所有5 個(gè)觀測(cè)點(diǎn)的數(shù)據(jù)進(jìn)行了分析，且每個(gè)觀測(cè)點(diǎn)從2 月1 日至3 月1 日進(jìn)行30 次1～24 h 短期預(yù)報(bào).對(duì)每個(gè)觀測(cè)站，我們利用30×24 個(gè)時(shí)間點(diǎn)的觀測(cè)值和預(yù)測(cè)值來(lái)計(jì)算該站的均方根誤差(RMSE)和歸一化均方根偏差(NRMSE)，以進(jìn)行兩種參數(shù)設(shè)置方案的誤差比較.兩種方案的綜合誤差比較結(jié)果如表1 所示.從表1 結(jié)果可以看出，對(duì)于本文隨機(jī)選取的5 個(gè)觀測(cè)站，“3+9”綜合預(yù)報(bào)誤差均要優(yōu)于“4+12”方案.將5 個(gè)觀測(cè)站誤差值求平均值可得，所有站RMSE 平均值從0.825 TECU 降低至0.805 TECU，降低了0.02 TECU，NRMSE 平均值從11.756%降低至11.452%，降低了約0.3%.需要指出的是，由于本文所選用TEC 數(shù)據(jù)的取值水平較低，所以改進(jìn)效果從絕對(duì)值上看較小，但本文經(jīng)過(guò)一定規(guī)模的統(tǒng)計(jì)分析顯示所有站點(diǎn)均有一定改進(jìn)，因此我們認(rèn)為這一改進(jìn)方案是依然有一定參考價(jià)值的.

表1 不同觀測(cè)站兩種方案預(yù)測(cè)誤差綜合比較

4 結(jié)果與討論

本文首先比較了自相關(guān)和自回歸滑動(dòng)平均(ARIMA)方法，結(jié)果顯示兩種方法均能較好地預(yù)測(cè)TEC周日變化.對(duì)本文預(yù)報(bào)結(jié)果進(jìn)行了RMSE 和NRMSE估算，結(jié)果表明總體上自相關(guān)預(yù)報(bào)誤差略低于ARIMA預(yù)報(bào)誤差，且在相同計(jì)算環(huán)境下，自相關(guān)預(yù)報(bào)所花費(fèi)的時(shí)間明顯少于ARIMA 模型.

其次，本文對(duì)自相關(guān)預(yù)報(bào)方法的誤差進(jìn)行了深入分析，重點(diǎn)分析了實(shí)際參與加權(quán)的觀測(cè)值覆蓋天數(shù)和觀測(cè)值數(shù)量這兩個(gè)參數(shù)的選取對(duì)預(yù)報(bào)誤差的影響.從Madrigal 數(shù)據(jù)庫(kù)中任意選取了5 個(gè)網(wǎng)格點(diǎn)的觀測(cè)數(shù)據(jù)開(kāi)展了預(yù)報(bào)試驗(yàn)，對(duì)不同參數(shù)設(shè)置方案的預(yù)報(bào)誤差進(jìn)行了對(duì)比分析.2020 年2 月8 日網(wǎng)格點(diǎn)4 (30°S，25°E)的個(gè)例分析結(jié)果表明，“3+9”方案的預(yù)報(bào)誤差相對(duì)最小，天數(shù)增加或減少(例如“2+6”、“4+12”和“5+12”等方案) 都會(huì)使誤差變大.進(jìn)一步地，我們對(duì)5 個(gè)觀測(cè)站分別進(jìn)行了30 天的1～24 h 短期預(yù)報(bào)試驗(yàn)，綜合分析了3 600 個(gè)預(yù)測(cè)點(diǎn)的RMSE 和NRMSE，結(jié)果顯示：相比于傳統(tǒng)的“4+12”方案，“3+9”方案在5 個(gè)觀測(cè)站中均取得了更優(yōu)的預(yù)報(bào)性能.具體地，所有站RMSE 平均值從0.825 TECU 降低至0.805 TECU，降低了0.02 TECU，NRMSE 平均值從11.756%降低至11.452%，降低了約0.3%.

因此，我們推測(cè)，TEC 時(shí)間序列的當(dāng)前狀態(tài)可能主要與前3 天的狀態(tài)有關(guān).當(dāng)參與加權(quán)的觀測(cè)天數(shù)增加(大于3 天)時(shí)，內(nèi)在關(guān)聯(lián)性小的觀測(cè)值參與加權(quán)運(yùn)算，可能會(huì)出現(xiàn)較大的狀態(tài)差異，從而造成誤差增大.當(dāng)參與加權(quán)的觀測(cè)天數(shù)減少時(shí)，可能會(huì)丟失主要的內(nèi)在關(guān)聯(lián)信息，從而也引起誤差增大.但需要指出的是，本文提出的參數(shù)優(yōu)化方案只是對(duì)傳統(tǒng)參數(shù)方案的一個(gè)微小的改進(jìn)，還無(wú)法明顯改善電離層預(yù)報(bào)性能，但可作為工程實(shí)現(xiàn)的參考方案.

電離層隨時(shí)間的變化非常復(fù)雜，除了熟知的周日變化、季節(jié)變化和太陽(yáng)活動(dòng)性11 年周期變化外，還存在很多小尺度的時(shí)間變化，目前我們對(duì)這些小尺度時(shí)間變化規(guī)律和控制因素的認(rèn)識(shí)還不足.若要顯著提高電離層短期預(yù)報(bào)性能，必須深入研究電離層參量隨時(shí)間短期變化的規(guī)律和內(nèi)在控制機(jī)理，這將是我們下一步研究工作的重點(diǎn).