于 燁張慧君李孝輝肖 波陳婧亞
(1中國科學院國家授時中心西安710600)
(2中國科學院精密導航定位與定時技術(shù)重點實驗室西安710600)
(3中國科學院大學北京100049)
(4中國科學院大學天文與空間科學學院北京100049)
在全球衛(wèi)星導航定位系統(tǒng)應用中,時間同步是一項關(guān)鍵技術(shù),其精度直接影響導航、定位和授時的精度[1?2].在衛(wèi)星自主導航中,系統(tǒng)要求地面預報一定時長的衛(wèi)星星歷和鐘差作為先驗信息,以便導航衛(wèi)星在沒有地面站支撐的條件下,通過一些手段自主完成軌道的確定和廣播星歷的播發(fā)[3?4].但是,由于星載原子鐘極易受到外界和自身因素的影響而很難掌握其復雜細致的變化規(guī)律.因此,建立精確的原子鐘運行模型變得非常困難,相應地準確預報衛(wèi)星鐘差也非常困難.
目前,國內(nèi)外許多學者在衛(wèi)星鐘差預報方面開展了廣泛深入的研究,提出的衛(wèi)星鐘差預報模型主要有:二次多項式模型(Quadratic Polynomial Model,QPM)、灰色預報模型(Grey Model,GM(1,1))、最小二乘支持向量機(Least Squares Support Vector Machines,LS-SVM)、自回歸滑動平均(Auto-Regressive Moving Average,ARMA)模型、人工神經(jīng)網(wǎng)絡(Artificial Neural Network,ANN)模型、卡爾曼濾波(Kalman Filter,KF)模型等[5?10].這些預報方法分別適用于不同條件下導航衛(wèi)星原子鐘鐘差的短期、中期和長期預報,但也均有各自的適用范圍和局限性.
本文針對衛(wèi)星鐘差呈現(xiàn)趨勢項和隨機項變化的特點,提出基于GM(1,1)和MECM的組合預報模型.該模型首先采用GM(1,1)預報鐘差的趨勢項部分,然后利用MECM模型對GM(1,1)殘差序列進行建模后預報,最后將GM(1,1)和MECM模型的預報結(jié)果對應相加得到鐘差的最終預報值.此外,采用IGS公布的事后精密衛(wèi)星鐘差產(chǎn)品進行預報試驗,通過與衛(wèi)星鐘差預報中常用的二次多項式模型和使用MECM模型預報結(jié)果的對比分析,驗證了該方法的優(yōu)越性和有效性.
在灰色系統(tǒng)理論中,GM(1,1)預報模型是最常用的一種灰色系統(tǒng)模型,它是由一個僅包含單變量的1階微分方程所構(gòu)成的預報模型,適合對自身數(shù)據(jù)的預報,而且具有只需要少量的數(shù)據(jù)建模的優(yōu)點[11].設有序列為X(0)={x(0)(1),x(0)(2),···,x(0)(n)},通過一次累加生成新的數(shù)據(jù)序列X(1):
對累加序列X(1)建立1階微分方程:
其中,參數(shù)a為發(fā)展系數(shù),參數(shù)u為灰作用量.對(2)式進行離散化處理,可得矩陣方程:
其中
根據(jù)最小二乘法可得矩陣方程(3)的最小二乘解為:
求解時間響應函數(shù)(5)式的解為:
由于序列X(1)是序列X(0)的累加序列,所以原始序列X(0)的預報模型為:
其中,k為參與預報的原始序列數(shù)據(jù)的個數(shù).原始序列的一般形式為:
其中,p為預報點數(shù).由以上預報模型即可對未來任意時刻的數(shù)據(jù)序列進行預報.
修正指數(shù)曲線法預報模型描述的是數(shù)據(jù)序列在一段時間內(nèi)按指數(shù)曲線增長,隨著時間的推移,增長趨勢會減緩[12].設觀測序列為:y1,y2,y3,···,yn,所對應的時刻為:t1,t2,t3,···,tn.可以建立修正指數(shù)曲線模型如下:
其中,這3個參數(shù)K、d、b均需要用歷史數(shù)據(jù)來確定.這里當K值可預先確定時,采用最小二乘法確定模型中的參數(shù).而當K值不能預先確定時,可采用三和法來估計參數(shù)K、d、b.
把n個觀測數(shù)據(jù)三等分,每個部分有m個,且n=3m.具體劃分如下:
令每部分的趨勢值之和等于相應的觀測值之和.三和法參數(shù)估計的基本步驟可歸結(jié)為如下過程:
其中,Si(i=1,2,3)表示數(shù)據(jù)觀測值的各部分之和,且
因為
所以,
由(12)–(14)式可得:
利用該模型即可對未來任意時刻的數(shù)據(jù)序列進行預報.
衛(wèi)星鐘差數(shù)據(jù)序列可以認為是由趨勢項和隨機項兩部分疊加而成,針對它的這一特點,提出一種基于GM(1,1)和MECM的組合預報模型.運用此組合模型去預報衛(wèi)星鐘差的預報流程如圖1所示,具體步驟如下所述:
(1)衛(wèi)星鐘差觀測數(shù)據(jù)預處理.首先剔除衛(wèi)星鐘差觀測序列中的異常鐘差數(shù)據(jù),然后采用拉格朗日插值法把缺失的鐘差數(shù)據(jù)補齊.
(2)利用GM(1,1)預報模型提取衛(wèi)星鐘差序列的趨勢項部分,并對其進行外推預報.
(3)利用MECM預報模型對GM(1,1)模型的建模殘差,即衛(wèi)星鐘差中的隨機項部分進行建模并外推預報.
(4)將利用GM(1,1)模型預報的趨勢項序列和利用MECM模型預報的對應隨機項序列相加,即得到衛(wèi)星鐘差序列最終的預報值.
圖1 基于GM(1,1)和MECM組合模型的衛(wèi)星鐘差預報流程Fig.1 Satellite clock bias prediction processing based on a combined model of GM(1,1)and MECM combination model
為了驗證該組合模型的有效性,從IGS服務器(ftp://cddis.gsfc.nasa.gov)上下載了2017年3月5日至6日共兩天的IGS事后精密鐘差產(chǎn)品進行預報試驗,其采樣間隔為15min.考慮到目前在軌的GPS衛(wèi)星有30多顆且我國北斗二代系統(tǒng)搭載的是銣原子鐘,所以我們選3顆裝載銣原子鐘的衛(wèi)星進行預報試驗.本算例隨機選取了1顆GPS Block IIF-M型衛(wèi)星PRN17、1顆GPSBlock IIF型衛(wèi)星PRN03和1顆GPSBlock IIR型衛(wèi)星PRN16的鐘差數(shù)據(jù)進行預報試驗.
3.2.1 方案一:用12h的鐘差數(shù)據(jù)建模后預報
選用二次多項式模型(QPM)、修正指數(shù)曲線法預報模型(MECM)、基于GM(1,1)和MECM的組合模型(GM(1,1)+MECM),用12h的鐘差數(shù)據(jù)建模,分別對未來6h、12h、18h和24h的衛(wèi)星鐘差進行中短期預報.
3.2.2 方案二:用24h的鐘差數(shù)據(jù)建模后預報
選用二次多項式模型(QPM)、修正指數(shù)曲線法預報模型(MECM)、基于GM(1,1)和MECM的組合模型(GM(1,1)+MECM),用24h的鐘差數(shù)據(jù)建模,分別對未來6h、12h、18h和24h的衛(wèi)星鐘差進行中短期預報.
使用方案一和方案二對PRN17號(IIR-M-Rb)、PRN03號(IIF-Rb)和PRN16號(IIRRb)衛(wèi)星鐘差進行建模后預報.將接下來6h、12h、18h和24h的實際觀測鐘差數(shù)據(jù)與各模型的預報鐘差數(shù)據(jù)相減就得到了預報誤差,取各模型預報誤差序列的均方根誤差(RMS)(具體計算公式見(17)式所示)為衡量預報精度的標準,結(jié)果如圖2–8和表1所示,其中圖2、圖4和圖6分別為用12h鐘差數(shù)據(jù)建立二次多項式模型、修正指數(shù)曲線法預報模型和基于GM(1,1)和MECM的組合模型去預報未來6h(a)、12h(b)、18h(c)和24h(d)鐘差數(shù)據(jù)的預報誤差變化圖,圖3、圖5和圖7分別為用24h鐘差數(shù)據(jù)建立二次多項式模型、修正指數(shù)曲線法預報模型和基于GM(1,1)和MECM的組合模型預報未來6h(a)、12h(b)、18h(c)和24h(d)鐘差數(shù)據(jù)的預報誤差變化圖.圖8為方案一(a)和方案二(b)各模型預報衛(wèi)星鐘差的平均預報精度的對比圖.
由于本試驗使用的是IGS服務器上公布的事后精密鐘差產(chǎn)品,其自身的誤差小于0.1ns,故可以作為“真值”,使用RMS作為統(tǒng)計量,去檢驗二次多項式模型、修正指數(shù)曲線法預報模型和基于GM(1,1)與MECM的組合模型所預報結(jié)果的好壞程度,計算公式如下:
圖2 PRN17用12 h數(shù)據(jù)建模預報6 h(a)、12 h(b)、18 h(c)和24 h(d)的預報誤差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ù)建模預報6 h(a)、12 h(b)、18 h(c)和24 h(d)的預報誤差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ù)建模預報6 h(a)、12 h(b)、18 h(c)和24 h(d)的預報誤差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ù)建模預報6 h(a)、12 h(b)、18 h(c)和24 h(d)的預報誤差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ù)建模預報6 h(a)、12 h(b)、18 h(c)和24 h(d)的預報誤差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ù)建模預報6 h(a)、12 h(b)、18 h(c)和24 h(d)的預報誤差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)星鐘差預報結(jié)果的RMS(單位:ns)Table 1 RMS of pred iction results for precise satellite clock bias(unit:ns)
圖8 方案一(a)和方案二(b)鐘差預報精度比較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ù))建模,對未來6h、12h、18h和24h的鐘差進行預報時,二次多項式模型預報誤差的平均均方差分別為1.01ns、2.08ns、3.79ns和5.96ns,修正指數(shù)曲線法模型預報誤差的平均均方差分別為0.87ns、1.77ns、3.19ns和5.00ns,而基于GM(1,1)和MECM組合模型預報誤差的平均均方差分別為0.43ns、0.63ns、0.74ns和0.79ns.
(2)采用中等數(shù)據(jù)量(24h鐘差數(shù)據(jù))建模,對未來6h、12h、18h和24h的鐘差進行預報時,二次多項式模型預報誤差的平均均方差分別為0.85ns、1.37ns、2.27ns和3.16ns,修正指數(shù)曲線法模型預報誤差的平均均方差分別為1.10ns、1.81ns、2.95ns和4.11ns,而基于GM(1,1)和MECM組合模型預報誤差的平均均方差分別為0.57ns、0.61ns、1.02ns和1.48ns.
(3)建模所采用的數(shù)據(jù)量的多少對二次多項式模型的預報性能有一定影響,建模所采用的數(shù)據(jù)量增大,二次多項式模型的預報性能有一定提升.建模所采用的數(shù)據(jù)量對修正指數(shù)曲線法模型的預報性能也有一定影響,建模所采用的數(shù)據(jù)量越大,修正指數(shù)曲線法的預報性能反而下降,這說明對于修正指數(shù)曲線法模型來說,不是建模所采用的數(shù)據(jù)量越多越好.而建模所采用的數(shù)據(jù)量的多少對基于GM(1,1)和MECM組合模型的預報性能影響較小,當增加一定量的建模數(shù)據(jù)時,基于GM(1,1)和MECM組合模型的預報性能變化很小,依然能夠進行高精度的中短期預報.
針對衛(wèi)星鐘差呈現(xiàn)趨勢項和隨機項變化的特點,本文提出了基于GM(1,1)和MECM的組合預報模型.該模型首先采用GM(1,1)預報鐘差的趨勢項,然后利用MECM模型對GM(1,1)模型殘差序列進行建模和預報,最后將GM(1,1)和MECM模型的預報結(jié)果對應相加得到衛(wèi)星鐘差的最終預報值.這種組合模型可以把兩種模型的預報優(yōu)勢有機地結(jié)合起來,進一步提高了模型的預報性能,從而提高了衛(wèi)星鐘差預報的精度.經(jīng)理論和預報試驗分析,結(jié)果也表明了該組合模型的有效性和可行性.