基于GM(1,1)和M ECM組合模型的GPS衛(wèi)星鐘差中短期預(yù)報(bào)?

于 燁張慧君李孝輝肖 波陳婧亞

(1中國科學(xué)院國家授時(shí)中心西安710600)

(2中國科學(xué)院精密導(dǎo)航定位與定時(shí)技術(shù)重點(diǎn)實(shí)驗(yàn)室西安710600)

(3中國科學(xué)院大學(xué)北京100049)

(4中國科學(xué)院大學(xué)天文與空間科學(xué)學(xué)院北京100049)

1 引言

在全球衛(wèi)星導(dǎo)航定位系統(tǒng)應(yīng)用中,時(shí)間同步是一項(xiàng)關(guān)鍵技術(shù),其精度直接影響導(dǎo)航、定位和授時(shí)的精度[1?2].在衛(wèi)星自主導(dǎo)航中,系統(tǒng)要求地面預(yù)報(bào)一定時(shí)長的衛(wèi)星星歷和鐘差作為先驗(yàn)信息,以便導(dǎo)航衛(wèi)星在沒有地面站支撐的條件下,通過一些手段自主完成軌道的確定和廣播星歷的播發(fā)[3?4].但是,由于星載原子鐘極易受到外界和自身因素的影響而很難掌握其復(fù)雜細(xì)致的變化規(guī)律.因此,建立精確的原子鐘運(yùn)行模型變得非常困難,相應(yīng)地準(zhǔn)確預(yù)報(bào)衛(wèi)星鐘差也非常困難.

目前,國內(nèi)外許多學(xué)者在衛(wèi)星鐘差預(yù)報(bào)方面開展了廣泛深入的研究,提出的衛(wèi)星鐘差預(yù)報(bào)模型主要有:二次多項(xiàng)式模型(Quadratic Polynomial Model,QPM)、灰色預(yù)報(bào)模型(Grey Model,GM(1,1))、最小二乘支持向量機(jī)(Least Squares Support Vector Machines,LS-SVM)、自回歸滑動(dòng)平均(Auto-Regressive Moving Average,ARMA)模型、人工神經(jīng)網(wǎng)絡(luò)(Artificial Neural Network,ANN)模型、卡爾曼濾波(Kalman Filter,KF)模型等[5?10].這些預(yù)報(bào)方法分別適用于不同條件下導(dǎo)航衛(wèi)星原子鐘鐘差的短期、中期和長期預(yù)報(bào),但也均有各自的適用范圍和局限性.

本文針對(duì)衛(wèi)星鐘差呈現(xiàn)趨勢(shì)項(xiàng)和隨機(jī)項(xiàng)變化的特點(diǎn),提出基于GM(1,1)和MECM的組合預(yù)報(bào)模型.該模型首先采用GM(1,1)預(yù)報(bào)鐘差的趨勢(shì)項(xiàng)部分,然后利用MECM模型對(duì)GM(1,1)殘差序列進(jìn)行建模后預(yù)報(bào),最后將GM(1,1)和MECM模型的預(yù)報(bào)結(jié)果對(duì)應(yīng)相加得到鐘差的最終預(yù)報(bào)值.此外,采用IGS公布的事后精密衛(wèi)星鐘差產(chǎn)品進(jìn)行預(yù)報(bào)試驗(yàn),通過與衛(wèi)星鐘差預(yù)報(bào)中常用的二次多項(xiàng)式模型和使用MECM模型預(yù)報(bào)結(jié)果的對(duì)比分析,驗(yàn)證了該方法的優(yōu)越性和有效性.

2 組合預(yù)報(bào)模型

2.1 GM(1,1)預(yù)報(bào)模型

在灰色系統(tǒng)理論中,GM(1,1)預(yù)報(bào)模型是最常用的一種灰色系統(tǒng)模型,它是由一個(gè)僅包含單變量的1階微分方程所構(gòu)成的預(yù)報(bào)模型,適合對(duì)自身數(shù)據(jù)的預(yù)報(bào),而且具有只需要少量的數(shù)據(jù)建模的優(yōu)點(diǎn)[11].設(shè)有序列為X(0)={x(0)(1),x(0)(2),···,x(0)(n)},通過一次累加生成新的數(shù)據(jù)序列X(1):

對(duì)累加序列X(1)建立1階微分方程:

其中,參數(shù)a為發(fā)展系數(shù),參數(shù)u為灰作用量.對(duì)(2)式進(jìn)行離散化處理,可得矩陣方程:

其中

根據(jù)最小二乘法可得矩陣方程(3)的最小二乘解為:

求解時(shí)間響應(yīng)函數(shù)(5)式的解為:

由于序列X(1)是序列X(0)的累加序列,所以原始序列X(0)的預(yù)報(bào)模型為:

其中,k為參與預(yù)報(bào)的原始序列數(shù)據(jù)的個(gè)數(shù).原始序列的一般形式為:

其中,p為預(yù)報(bào)點(diǎn)數(shù).由以上預(yù)報(bào)模型即可對(duì)未來任意時(shí)刻的數(shù)據(jù)序列進(jìn)行預(yù)報(bào).

2.2 MECM預(yù)報(bào)模型

修正指數(shù)曲線法預(yù)報(bào)模型描述的是數(shù)據(jù)序列在一段時(shí)間內(nèi)按指數(shù)曲線增長,隨著時(shí)間的推移,增長趨勢(shì)會(huì)減緩[12].設(shè)觀測(cè)序列為:y1,y2,y3,···,yn,所對(duì)應(yīng)的時(shí)刻為:t1,t2,t3,···,tn.可以建立修正指數(shù)曲線模型如下:

其中,這3個(gè)參數(shù)K、d、b均需要用歷史數(shù)據(jù)來確定.這里當(dāng)K值可預(yù)先確定時(shí),采用最小二乘法確定模型中的參數(shù).而當(dāng)K值不能預(yù)先確定時(shí),可采用三和法來估計(jì)參數(shù)K、d、b.

把n個(gè)觀測(cè)數(shù)據(jù)三等分,每個(gè)部分有m個(gè),且n=3m.具體劃分如下:

令每部分的趨勢(shì)值之和等于相應(yīng)的觀測(cè)值之和.三和法參數(shù)估計(jì)的基本步驟可歸結(jié)為如下過程:

其中,Si(i=1,2,3)表示數(shù)據(jù)觀測(cè)值的各部分之和,且

因?yàn)?/p>

所以,

由(12)–(14)式可得:

利用該模型即可對(duì)未來任意時(shí)刻的數(shù)據(jù)序列進(jìn)行預(yù)報(bào).

2.3 預(yù)報(bào)過程

衛(wèi)星鐘差數(shù)據(jù)序列可以認(rèn)為是由趨勢(shì)項(xiàng)和隨機(jī)項(xiàng)兩部分疊加而成,針對(duì)它的這一特點(diǎn),提出一種基于GM(1,1)和MECM的組合預(yù)報(bào)模型.運(yùn)用此組合模型去預(yù)報(bào)衛(wèi)星鐘差的預(yù)報(bào)流程如圖1所示,具體步驟如下所述:

(1)衛(wèi)星鐘差觀測(cè)數(shù)據(jù)預(yù)處理.首先剔除衛(wèi)星鐘差觀測(cè)序列中的異常鐘差數(shù)據(jù),然后采用拉格朗日插值法把缺失的鐘差數(shù)據(jù)補(bǔ)齊.

(2)利用GM(1,1)預(yù)報(bào)模型提取衛(wèi)星鐘差序列的趨勢(shì)項(xiàng)部分,并對(duì)其進(jìn)行外推預(yù)報(bào).

(3)利用MECM預(yù)報(bào)模型對(duì)GM(1,1)模型的建模殘差,即衛(wèi)星鐘差中的隨機(jī)項(xiàng)部分進(jìn)行建模并外推預(yù)報(bào).

(4)將利用GM(1,1)模型預(yù)報(bào)的趨勢(shì)項(xiàng)序列和利用MECM模型預(yù)報(bào)的對(duì)應(yīng)隨機(jī)項(xiàng)序列相加,即得到衛(wèi)星鐘差序列最終的預(yù)報(bào)值.

圖1 基于GM(1,1)和MECM組合模型的衛(wèi)星鐘差預(yù)報(bào)流程Fig.1 Satellite clock bias prediction processing based on a combined model of GM(1,1)and MECM combination model

3 算例與分析

3.1 數(shù)據(jù)來源

為了驗(yàn)證該組合模型的有效性,從IGS服務(wù)器(ftp://cddis.gsfc.nasa.gov)上下載了2017年3月5日至6日共兩天的IGS事后精密鐘差產(chǎn)品進(jìn)行預(yù)報(bào)試驗(yàn),其采樣間隔為15min.考慮到目前在軌的GPS衛(wèi)星有30多顆且我國北斗二代系統(tǒng)搭載的是銣原子鐘,所以我們選3顆裝載銣原子鐘的衛(wèi)星進(jìn)行預(yù)報(bào)試驗(yàn).本算例隨機(jī)選取了1顆GPS Block IIF-M型衛(wèi)星PRN17、1顆GPSBlock IIF型衛(wèi)星PRN03和1顆GPSBlock IIR型衛(wèi)星PRN16的鐘差數(shù)據(jù)進(jìn)行預(yù)報(bào)試驗(yàn).

3.2 建模方案

3.2.1 方案一:用12h的鐘差數(shù)據(jù)建模后預(yù)報(bào)

選用二次多項(xiàng)式模型(QPM)、修正指數(shù)曲線法預(yù)報(bào)模型(MECM)、基于GM(1,1)和MECM的組合模型(GM(1,1)+MECM),用12h的鐘差數(shù)據(jù)建模,分別對(duì)未來6h、12h、18h和24h的衛(wèi)星鐘差進(jìn)行中短期預(yù)報(bào).

3.2.2 方案二:用24h的鐘差數(shù)據(jù)建模后預(yù)報(bào)

選用二次多項(xiàng)式模型(QPM)、修正指數(shù)曲線法預(yù)報(bào)模型(MECM)、基于GM(1,1)和MECM的組合模型(GM(1,1)+MECM),用24h的鐘差數(shù)據(jù)建模,分別對(duì)未來6h、12h、18h和24h的衛(wèi)星鐘差進(jìn)行中短期預(yù)報(bào).

3.3 預(yù)報(bào)結(jié)果與分析

使用方案一和方案二對(duì)PRN17號(hào)(IIR-M-Rb)、PRN03號(hào)(IIF-Rb)和PRN16號(hào)(IIRRb)衛(wèi)星鐘差進(jìn)行建模后預(yù)報(bào).將接下來6h、12h、18h和24h的實(shí)際觀測(cè)鐘差數(shù)據(jù)與各模型的預(yù)報(bào)鐘差數(shù)據(jù)相減就得到了預(yù)報(bào)誤差,取各模型預(yù)報(bào)誤差序列的均方根誤差(RMS)(具體計(jì)算公式見(17)式所示)為衡量預(yù)報(bào)精度的標(biāo)準(zhǔn),結(jié)果如圖2–8和表1所示,其中圖2、圖4和圖6分別為用12h鐘差數(shù)據(jù)建立二次多項(xiàng)式模型、修正指數(shù)曲線法預(yù)報(bào)模型和基于GM(1,1)和MECM的組合模型去預(yù)報(bào)未來6h(a)、12h(b)、18h(c)和24h(d)鐘差數(shù)據(jù)的預(yù)報(bào)誤差變化圖,圖3、圖5和圖7分別為用24h鐘差數(shù)據(jù)建立二次多項(xiàng)式模型、修正指數(shù)曲線法預(yù)報(bào)模型和基于GM(1,1)和MECM的組合模型預(yù)報(bào)未來6h(a)、12h(b)、18h(c)和24h(d)鐘差數(shù)據(jù)的預(yù)報(bào)誤差變化圖.圖8為方案一(a)和方案二(b)各模型預(yù)報(bào)衛(wèi)星鐘差的平均預(yù)報(bào)精度的對(duì)比圖.

由于本試驗(yàn)使用的是IGS服務(wù)器上公布的事后精密鐘差產(chǎn)品,其自身的誤差小于0.1ns,故可以作為“真值”,使用RMS作為統(tǒng)計(jì)量,去檢驗(yàn)二次多項(xiàng)式模型、修正指數(shù)曲線法預(yù)報(bào)模型和基于GM(1,1)與MECM的組合模型所預(yù)報(bào)結(jié)果的好壞程度,計(jì)算公式如下:

圖2 PRN17用12 h數(shù)據(jù)建模預(yù)報(bào)6 h(a)、12 h(b)、18 h(c)和24 h(d)的預(yù)報(bào)誤差Fig.2 PRN17 prediction errors of 6 h(a),12 h(b),18 h(c),and 24 h(d)by using 12 h data modeling

圖3 PRN17用24 h數(shù)據(jù)建模預(yù)報(bào)6 h(a)、12 h(b)、18 h(c)和24 h(d)的預(yù)報(bào)誤差Fig.3 PRN17 prediction errors of 6 h(a),12 h(b),18 h(c),and 24 h(d)by using 24 h data modeling

圖4 PRN03用12 h數(shù)據(jù)建模預(yù)報(bào)6 h(a)、12 h(b)、18 h(c)和24 h(d)的預(yù)報(bào)誤差Fig.4 PRN03 prediction errors of 6 h(a),12 h(b),18 h(c),and 24 h(d)by using 12 h data modeling

圖5 PRN03用24 h數(shù)據(jù)建模預(yù)報(bào)6 h(a)、12 h(b)、18 h(c)和24 h(d)的預(yù)報(bào)誤差Fig.5 PRN03 prediction errors of 6 h(a),12 h(b),18 h(c),and 24 h(d)by using 24 h data modeling

圖6 PRN16用12 h數(shù)據(jù)建模預(yù)報(bào)6 h(a)、12 h(b)、18 h(c)和24 h(d)的預(yù)報(bào)誤差Fig.6 PRN16 prediction errors of 6 h(a),12 h(b),18 h(c),and 24 h(d)by using 12 h data modeling

圖7 PRN16用24 h數(shù)據(jù)建模預(yù)報(bào)6 h(a)、12 h(b)、18 h(c)和24 h(d)的預(yù)報(bào)誤差Fig.7 PRN16 prediction errors of 6 h(a),12 h(b),18 h(c),and 24 h(d)by using 24 h data modeling

表1 精密衛(wèi)星鐘差預(yù)報(bào)結(jié)果的RMS(單位:ns)Table 1 RMS of pred iction results for precise satellite clock bias(unit:ns)

圖8 方案一(a)和方案二(b)鐘差預(yù)報(bào)精度比較Fig.8 Accuracy comparisons of clock bias prediction for scheme 1(a)and scheme 2(b)

結(jié)合圖2–8和表1可知:

(1)采用小數(shù)據(jù)量(12h鐘差數(shù)據(jù))建模,對(duì)未來6h、12h、18h和24h的鐘差進(jìn)行預(yù)報(bào)時(shí),二次多項(xiàng)式模型預(yù)報(bào)誤差的平均均方差分別為1.01ns、2.08ns、3.79ns和5.96ns,修正指數(shù)曲線法模型預(yù)報(bào)誤差的平均均方差分別為0.87ns、1.77ns、3.19ns和5.00ns,而基于GM(1,1)和MECM組合模型預(yù)報(bào)誤差的平均均方差分別為0.43ns、0.63ns、0.74ns和0.79ns.

(2)采用中等數(shù)據(jù)量(24h鐘差數(shù)據(jù))建模,對(duì)未來6h、12h、18h和24h的鐘差進(jìn)行預(yù)報(bào)時(shí),二次多項(xiàng)式模型預(yù)報(bào)誤差的平均均方差分別為0.85ns、1.37ns、2.27ns和3.16ns,修正指數(shù)曲線法模型預(yù)報(bào)誤差的平均均方差分別為1.10ns、1.81ns、2.95ns和4.11ns,而基于GM(1,1)和MECM組合模型預(yù)報(bào)誤差的平均均方差分別為0.57ns、0.61ns、1.02ns和1.48ns.

(3)建模所采用的數(shù)據(jù)量的多少對(duì)二次多項(xiàng)式模型的預(yù)報(bào)性能有一定影響,建模所采用的數(shù)據(jù)量增大,二次多項(xiàng)式模型的預(yù)報(bào)性能有一定提升.建模所采用的數(shù)據(jù)量對(duì)修正指數(shù)曲線法模型的預(yù)報(bào)性能也有一定影響,建模所采用的數(shù)據(jù)量越大,修正指數(shù)曲線法的預(yù)報(bào)性能反而下降,這說明對(duì)于修正指數(shù)曲線法模型來說,不是建模所采用的數(shù)據(jù)量越多越好.而建模所采用的數(shù)據(jù)量的多少對(duì)基于GM(1,1)和MECM組合模型的預(yù)報(bào)性能影響較小,當(dāng)增加一定量的建模數(shù)據(jù)時(shí),基于GM(1,1)和MECM組合模型的預(yù)報(bào)性能變化很小,依然能夠進(jìn)行高精度的中短期預(yù)報(bào).

4 結(jié)語

針對(duì)衛(wèi)星鐘差呈現(xiàn)趨勢(shì)項(xiàng)和隨機(jī)項(xiàng)變化的特點(diǎn),本文提出了基于GM(1,1)和MECM的組合預(yù)報(bào)模型.該模型首先采用GM(1,1)預(yù)報(bào)鐘差的趨勢(shì)項(xiàng),然后利用MECM模型對(duì)GM(1,1)模型殘差序列進(jìn)行建模和預(yù)報(bào),最后將GM(1,1)和MECM模型的預(yù)報(bào)結(jié)果對(duì)應(yīng)相加得到衛(wèi)星鐘差的最終預(yù)報(bào)值.這種組合模型可以把兩種模型的預(yù)報(bào)優(yōu)勢(shì)有機(jī)地結(jié)合起來,進(jìn)一步提高了模型的預(yù)報(bào)性能,從而提高了衛(wèi)星鐘差預(yù)報(bào)的精度.經(jīng)理論和預(yù)報(bào)試驗(yàn)分析,結(jié)果也表明了該組合模型的有效性和可行性.