陳略,唐歌實(shí),胡松杰,平勁松,許雪晴,夏金超

(1.北京航天飛行控制中心航天飛行動力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室;北京100094; 2.中國科學(xué)院國家天文臺;北京100012;3.中國科學(xué)院上海天文臺;上海200030)

高精度UT1-UTC差分預(yù)報(bào)方法研究

陳略1,唐歌實(shí)1,胡松杰1,平勁松2,許雪晴3,夏金超1

(1.北京航天飛行控制中心航天飛行動力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室;北京100094; 2.中國科學(xué)院國家天文臺;北京100012;3.中國科學(xué)院上海天文臺;上海200030)

針對衛(wèi)星導(dǎo)航所需的高精度地球定向參數(shù)(EOP)中的UT1-UTC預(yù)報(bào)問題,提出了基于雙差分LS+ AR的UT1-UTC參數(shù)預(yù)報(bào)方法。對UT1-UTC觀測數(shù)據(jù)進(jìn)行跳秒檢測、固體地球帶諧潮汐項(xiàng)改正,然后對改正后的UT1-UTC數(shù)據(jù)進(jìn)行雙差分處理,增強(qiáng)數(shù)據(jù)平穩(wěn)性;采用最小二乘擬合(LS)與自回歸(AR)分析方法對差分處理后的數(shù)據(jù)進(jìn)行分析與預(yù)報(bào);對預(yù)報(bào)結(jié)果進(jìn)行逆差分處理與潮汐項(xiàng)改正外推、跳秒恢復(fù),獲取高精度的UT1-UTC預(yù)報(bào)值。通過與國際EOP_PCC預(yù)報(bào)結(jié)果對比表明,UT1-UTC短期預(yù)報(bào)精度與EOP_PCC較優(yōu)的預(yù)報(bào)精度相當(dāng),其中1天UT1-UTC預(yù)報(bào)精度優(yōu)于0.03 ms,優(yōu)于EOP_PCC預(yù)報(bào)結(jié)果。介紹了北京航天飛行控制中心的UT1-UTC每日例行預(yù)報(bào)情況。

地球定向參數(shù);UT1-UTC預(yù)報(bào);差分處理;AR模型;潮汐項(xiàng)改正

0 引 言

高精度地球定向參數(shù)(Earth orientation parameters,EOP)是地球參考架與天球參考架實(shí)時轉(zhuǎn)換時的基礎(chǔ)參數(shù)[1]。VLBI、GNSS、DORIS、LLR、SLR等現(xiàn)代空間大地測量技術(shù)可精確地測量地球定向參數(shù)[2-3]。高精度EOP在載人航天、深空探測等任務(wù)中具有重要價值,尤其是在航天器實(shí)時、高精度導(dǎo)航方面不可或缺。目前,由于EOP解算結(jié)果在時間上有滯后,因此EOP參數(shù)高精度預(yù)報(bào)對于航天器導(dǎo)航顯得十分重要,且EOP預(yù)報(bào)對于與EOP相關(guān)的地球物理動力學(xué)理論分析具有重要價值。

EOP預(yù)報(bào)中的極移預(yù)報(bào)、UT1-UTC預(yù)報(bào)是EOP預(yù)報(bào)中的熱點(diǎn)研究內(nèi)容。目前,有多種預(yù)報(bào)方法應(yīng)用于EOP的預(yù)報(bào)中,如最小二乘(LS)外推法[4]、最小二乘聯(lián)合自回歸分析法(LS+AR)[5]、神經(jīng)網(wǎng)絡(luò)預(yù)報(bào)法[6]、譜分析最小二乘外推法[7]、小波分解與自回歸分析方法[8]等。國際地球自轉(zhuǎn)與參考系服務(wù)(IERS)組織于2005年10月發(fā)起了首次全球地球定向參數(shù)預(yù)報(bào)比賽活動(EOP_PCC),旨在號召全世界EOP預(yù)報(bào)人員運(yùn)用不同方法來預(yù)報(bào)EOP,用以分析現(xiàn)有方法的特點(diǎn)及使用情況。LS+AR預(yù)報(bào)方法是預(yù)報(bào)精度最高的方法之一。但運(yùn)用LS+ AR預(yù)報(bào)方法時,LS模型參數(shù)的選擇,AR模型階數(shù)的確定,非平穩(wěn)序列AR分析,EOP預(yù)報(bào)精度評估等是進(jìn)行EOP預(yù)報(bào)的關(guān)鍵問題?，F(xiàn)有LS+AR預(yù)報(bào)方法在對數(shù)據(jù)平穩(wěn)性要求方面的研究較為欠缺,在一定程度上限制了EOP預(yù)報(bào)精度。

本文從分析數(shù)據(jù)的平穩(wěn)特性出發(fā),提出了雙差分LS+AR預(yù)報(bào)方法,用于實(shí)現(xiàn)對EOP中UT1-UTC高精度預(yù)報(bào),并獲取了高精度的UT1-UTC短期預(yù)報(bào)結(jié)果。與此同時,同時介紹了與EOP_PCC進(jìn)行精度比較的結(jié)果,此外,介紹了北京航天飛行控制中心EOP參數(shù)高精度預(yù)報(bào)軟件每日例行情況。

1 預(yù)報(bào)理論與方法

1.1 最小二乘(LS)模型

對UT1-UTC觀測數(shù)據(jù)進(jìn)行最小二乘模型擬合時,模型包括線性項(xiàng)與周期項(xiàng)(周年項(xiàng)、半年項(xiàng)、9.3年項(xiàng)、18.6年項(xiàng)等),其最小二乘模型如式(1)所示。

式中:t為序列的UTC時間,單位為年;A,B,C,D1, D2,E1,E2,…為模型的擬合參數(shù);p1,p2,…為擬合周期,可根據(jù)實(shí)際情況確定模型形狀。

1.2 AR模型

1.2.1 AR模型定義

AR模型表示平穩(wěn)隨機(jī)序列xt(t=1,2,…,N)與t時刻以前的規(guī)律性變化和t時刻白噪聲的關(guān)系,其數(shù)學(xué)模型如式(2)所示

式中:φ1,φ2,…,φp為模型參數(shù);a為白噪聲;p為模型階數(shù),稱上式為p階自回歸模型,簡記為AR(p);為白噪聲的方差。

1.2.2 AR模型階數(shù)確定

運(yùn)用AR模型關(guān)鍵在于確定模型階數(shù)p,可采用最終預(yù)報(bào)誤差(final prediction error criterion, FPE)準(zhǔn)則、AIC(akake information criterion)準(zhǔn)則、BIC準(zhǔn)則、矩陣奇異值分解(singular value decomposition,SVD)準(zhǔn)則來確定模型的階數(shù)。

本文采用解Yule-Walker方程方法來估計(jì)模型階數(shù),其基本原理[9]如下:

根據(jù)自協(xié)方差函數(shù)的定義,一個平穩(wěn)、正態(tài)、零均值的隨機(jī)過程{xt}的自協(xié)方差定義為

其中,E表示集合平均算子,即數(shù)學(xué)期望,{xt}的零均值性導(dǎo)致其自協(xié)方差函數(shù)與自相關(guān)函數(shù)是一致的,而且作為Rk的特例,當(dāng)k=0時,得{xt}的方差函數(shù)

自相關(guān)系數(shù)定義為

顯然,ρ0=1。由許瓦茲不定式,知,的取值范圍為

對于平穩(wěn)時間序列,當(dāng)ρk=0時,表示序列前后無關(guān),即為白噪聲;當(dāng)時,表示序列前后取值的關(guān)系不但完全確定,而且前后取值相等,即序列為常值序列。

對式(2)中左右兩邊乘以xt-k,再取數(shù)學(xué)期望并除以R0,得到

分別令k=1,2,…,p,并注意自協(xié)方差函數(shù)是偶函數(shù)的性質(zhì),R-k=Rk,因而ρ-k=ρk,可得線性方程組如下

分別令k=1,2,…,n,則可得方程的個數(shù)等于未知數(shù)φ1,φ2,…,φn的個數(shù),即上式等效為式(9)。

此式成為Yule-Walker方程,其中ρ為n×1的列矩陣,T為n階仿陣,稱為Toeplitz矩陣,如式(10)所示。

T中各元素可按定義估計(jì)出。由于T是滿秩的,因此可容易求出模型參數(shù)如式(11)所示。

本文采用FPE準(zhǔn)則來確定AR模型階數(shù)。FPE準(zhǔn)則函數(shù)如式(12)所示。

1.3 精度評估

為了評估預(yù)報(bào)精度,采用EOP預(yù)報(bào)中的通用精度評估準(zhǔn)則-平均絕對誤差(mean absulute eroor, MAE)作為預(yù)報(bào)結(jié)果的精度指標(biāo),其計(jì)算公式如下:

式中:o為實(shí)際觀測值;p為預(yù)報(bào)值;i為預(yù)報(bào)跨度;n為預(yù)報(bào)期數(shù)。

2 雙差分LS+AR預(yù)報(bào)流程

雙差分LS+AR的UT1-UTC預(yù)報(bào)方法的實(shí)現(xiàn)流程如圖1所示。

圖1 雙差分LS+AR UT1-UTC預(yù)報(bào)流程圖Fig.1 UT1-UTC prediction process of differential LS+AR method

3 計(jì)算分析

利用IERS公布的EOP數(shù)據(jù)進(jìn)行預(yù)報(bào)與精度驗(yàn)證。選取EOP 05C04的極移數(shù)據(jù)進(jìn)行分析,便于與EOP_PCC結(jié)果進(jìn)行比較。時間跨度為1989年1月1日～2009年12月31日,數(shù)據(jù)采樣頻率為1點(diǎn)/天。預(yù)報(bào)的基礎(chǔ)時間序列為1989年1月1日～2009年10月2日,對2009年10月3日至2009年12月31日的EOP極移進(jìn)行預(yù)報(bào),并將預(yù)報(bào)值與真實(shí)觀測值進(jìn)行比較,評估預(yù)報(bào)精度。預(yù)報(bào)跨度分別從1～30天,預(yù)報(bào)期數(shù)為90期。具體預(yù)報(bào)結(jié)果如下所述。

圖2為1989年—2009年的UT1-UTC的時間變化圖,其中圖2(a)為原始的UT1-UTC波形圖;圖2(b)為剔除跳秒后的UT1-UTC圖;圖2(c)為固體地球帶諧潮汐項(xiàng)改正圖,依據(jù)IERS2010規(guī)范進(jìn)行改正;圖2(d)是經(jīng)過跳秒剔除與潮汐項(xiàng)改正后的UT1-UTC時間序列。

圖3為1989年—2009年UT1-UTC雙差分后的LS擬合結(jié)果。從圖中可以看出,經(jīng)過雙差分后的UT1-UTC序列具有明顯的平穩(wěn)性。利用LS所得參數(shù)進(jìn)行外推,利用AR模型對LS擬合后的殘差進(jìn)行預(yù)報(bào),將LS外推值與AR模型外推值進(jìn)行相加,即得跨度為1天的預(yù)報(bào)值。其中,在進(jìn)行AR預(yù)報(bào)過程中,AR模型階數(shù)的確定是AR預(yù)報(bào)最關(guān)鍵步驟,本文采用FPE準(zhǔn)則確定模型階數(shù),在階數(shù)確定過程中采用迭代搜索方法確定最優(yōu)模型階數(shù),階數(shù)搜索范圍為1～100,采用FPE準(zhǔn)則計(jì)算出每一項(xiàng)FPE值,其中最小FPE值對應(yīng)階數(shù)即為最優(yōu)的AR模型階數(shù)。

按照圖1的處理流程,對預(yù)報(bào)初值進(jìn)行逆雙差分處理、固體地球帶諧潮汐項(xiàng)改正外推、跳秒恢復(fù)后,即可獲得第一天的UT1-UTC預(yù)報(bào)值。在得到第1天的預(yù)報(bào)值之后,將第1天的預(yù)報(bào)值加入到基礎(chǔ)序列中,重新按照圖1的處理流程進(jìn)行第2天的UT1-UTC預(yù)報(bào)。以此類推,可獲取不同預(yù)報(bào)跨度的UT1-UTC預(yù)報(bào)結(jié)果。

圖2 1989-2009年UT1-UTC的時間序列分布圖Fig.2 Time sequence of UT1-UTC from 1989 to 2009

圖3 1989-2009年UT1-UTC雙差分LS擬合結(jié)果Fig.3 Differential LS fitting results of UT1-UTC from 1989 to 2009

基于以上方法進(jìn)行90期預(yù)報(bào),并利用平均絕對誤差(MAE)準(zhǔn)則評估UT1-UTC預(yù)報(bào)精度,其結(jié)果如圖4所示,具體數(shù)值如表1所示。

圖4 BACC的UT1-UTC短期預(yù)報(bào)結(jié)果Fig.4 Short term prediction results of UT1-UTC in Beijing Aerospace Control Center

表1 BACC的UT1-UTC短期預(yù)報(bào)MAE結(jié)果Table 1 MAE of UT1-UTC short term prediction in Beijing Aerospace Control Center

全球性的EOP預(yù)報(bào)比對運(yùn)動(earth orientation parameters prediction comparison campaign,EOP_PCC),共有8個國家12名時間序列分析領(lǐng)域的頂尖級專家學(xué)者參與,涉及20余中預(yù)報(bào)方法。EOP_PCC運(yùn)動分為10天超短期預(yù)報(bào), 30天短期預(yù)報(bào),500天長期預(yù)報(bào)。

本文將預(yù)報(bào)結(jié)果與EOP_PCC的結(jié)果進(jìn)行比對分析,分析進(jìn)行超短期預(yù)報(bào)比對,與短期預(yù)報(bào)比對。圖5為EOP_PCC的預(yù)報(bào)結(jié)果,圖6為本文的預(yù)報(bào)結(jié)果。

從結(jié)果可以看出,本文的UT1-UTC短期預(yù)報(bào)精度與EOP_PCC的預(yù)報(bào)精度相當(dāng)。對于1天預(yù)報(bào)精度,本文UT1-UTC預(yù)報(bào)精度優(yōu)于EOP_PCC預(yù)報(bào)精度。本文獲取的1天UT1-UTC預(yù)報(bào)精度優(yōu)于0.03 ms。EOP_PCC的UT1-UTC一天預(yù)報(bào)誤差[1]在0.08 ms水平。

圖5 EOP_PCC的UT1-UTC短期預(yù)報(bào)結(jié)果Fig.5 Short term prediction results of UT1-UTC in EOP_PCC

圖6 BACC的UT1-UTC短期預(yù)報(bào)結(jié)果Fig.6 Short term prediction results of UT1-UTC in Beijing Aerospace Control Center

4 UT1-UTC每日例行預(yù)報(bào)

北京航天飛行控制中心(BACC)飛行動力學(xué)重點(diǎn)實(shí)驗(yàn)室開發(fā)了EOP高精度預(yù)報(bào)軟件,可實(shí)現(xiàn)對UT1-UTC每日例行預(yù)報(bào)與發(fā)布。該軟件每日從互聯(lián)網(wǎng)上自動化下載UT1-UTC解算數(shù)據(jù),而后進(jìn)行UT1-UTC預(yù)報(bào)并生成預(yù)報(bào)結(jié)果,同時輸出BACC/ IERS/美國海軍天文臺(USNO)的UT1-UTC預(yù)報(bào)的比較結(jié)果,用以每日評估BACC的UT1-UTC預(yù)報(bào)。BACC的EOP預(yù)報(bào)軟件輸入的UT1-UTC解算數(shù)據(jù)由兩個部分組成:IERS網(wǎng)站公布的EOP 08 C04序列+USNO最近一個月的解算序列。

圖7 2014-11-03 UT1-UTC 30天預(yù)報(bào)結(jié)果比較Fig.7 UT1-UTC 30 days prediction comparing results on 2014-11-03

圖8 2014-11-04 UT1-UTC 30天預(yù)報(bào)結(jié)果比較Fig.8 UT1-UTC 30 days prediction comparing results on 2014-11-04

圖9 2014-11-05 UT1-UTC 30天預(yù)報(bào)結(jié)果比較Fig.9 UT1-UTC 30 days prediction comparing results on 2014-11-05

圖7～圖9顯示了2014年11月3—5日, BACC的UT1-UTC預(yù)報(bào)結(jié)果與IERS/USNO結(jié)果對比圖。在圖7～圖9中,圓圈表示BACC的預(yù)報(bào)結(jié)果,正方形表示USNO的預(yù)報(bào)結(jié)果,梅花形表示IERS的預(yù)報(bào)結(jié)果。從以上圖中可以看出,BACC的UT1-UTC預(yù)報(bào)結(jié)果趨勢與IERS、USNO基本一致,且BACC 30天的UT1-UTC預(yù)報(bào)結(jié)果基本位于IERS與USNO之間,這說明了北京中心UT1-UTC預(yù)報(bào)結(jié)果有效性。

5 結(jié)束語

本文提出了基于雙差分LS+AR的UT1-UTC預(yù)報(bào)方法,詳細(xì)介紹了理論方法與處理流程,并介紹了BACC的UT1-UTC每日例行預(yù)報(bào)情況。通過對IERS實(shí)際的EOP中UT1-UTC進(jìn)行預(yù)報(bào)與比對分析,本文的UT1-UTC預(yù)報(bào)精度與EOP_PCC較優(yōu)的短期預(yù)報(bào)精度相當(dāng)。本文獲取的1天UT1-UTC預(yù)報(bào)精度優(yōu)于0.03 ms,優(yōu)于EOP_PCC每日預(yù)報(bào)結(jié)果。

相比于傳統(tǒng)LS+AR的地球定向參數(shù)預(yù)報(bào)方法,本文重點(diǎn)從AR模型輸入數(shù)據(jù)的平穩(wěn)性要求方面進(jìn)行考慮,采用了雙差分預(yù)處理方法,用以提高UT1-UTC數(shù)據(jù)的平穩(wěn)性,使其更加適應(yīng)利用AR模型進(jìn)行預(yù)報(bào)。此外,AR模型階數(shù)確定采用了迭代搜索方式選擇最優(yōu)模型階數(shù)。通過以上兩項(xiàng)步驟,可有效提高對UT1-UTC的高精度預(yù)報(bào)。

后續(xù)將進(jìn)一步研究本文提出的雙差分LS+AR預(yù)報(bào)方法在極移、LOD預(yù)報(bào)上的應(yīng)用,同時將進(jìn)一步評估該方法在對EOP參數(shù)長期預(yù)報(bào)中的精度與適用性。

[1]Kalarus M,Schuh H,Kosek W,et al.Achievements of the Earth orientation parameters prediction comparison campaign [J].Journal of Geodynamics,2010(62):587-596.

[2]Tobias Nilsson,Robert Heinkelmann,Maria Karbon,et al. Earth orientation parameters estimated from VLBI during the CONT11 campaign[J].Journal of Geodesy,2014(88): 491-502.

[3]許雪晴,周永宏.地球定向參數(shù)高精度預(yù)報(bào)方法研究[J].飛行器測控學(xué)報(bào),2010,29(2):70-76.[Xu X Q,Zhou Y H. High precision prediction method of Earth orientation parameters[J].Journal of Spacecraft TT&C Technology, 2010,29(2):70-76.]

[4]Xu X,Zhou Y,Liao X.Short-term earth orientation parameters predictions by combination of the least-squares, AR model and Kalman filter[J].Journal of Geodynamics, 2010(62):83-86.

[5]Guo J Y,Li Y B,Dai C L,et al.A technique to improve the accuracy of Earth orientation prediction algorithms based on least squares extrapolation[J].Journal of Geodynamics, 2013(70):36-48.

[6]Wang Q J,Liao D C,Zhou Y H.Real-time rapid prediction of variations of Earth's rotational rate[J].Chin.Sci.Bull, 2008(53):969-973.

[7]Akulenko L D,Kumakshev S A,Rykhlova L V.A model for the polar motion of the deformable Earth adequate for astrometric data[J].Astron.Rep.,2002(46):74-82.

[8]Kosek W,Kalarus M,Johnson T J,et al.A comparison of LOD and UT1-UTC forecasts by different combination prediction techniques[J].Artif.Satell,2005(40):119-125.

[9]楊叔子.時間序列分析的工程應(yīng)用[M].武漢:華中科技大學(xué)出版社,2007:184-189.[Yang S Z.Time series analysis in engineering application[M].Wuhan:Huazhong University of Science&Technology Press,2007:184-189.]

通信地址:北京海淀區(qū)5130信箱120號(100094)

電話:(010)66362814

E-mail:luechen0912@163.com

[責(zé)任編輯:宋宏]

High Accuracy Differential Prediction of UT1-UTC

CHEN Lue1,TANG Geshi1,HU Songjie1,PING Jinsong2,XU Xueqing3,XIA Jinchao1
(1.National Key Laboratory of Science and Technology on Aerospace Flight Dynamic,Beijing Aerospace Control Center, 100094 Beijing,China;2.National Astronomical Observatory,Chinese Academy of Sciences,100012 Beijing,China; 3.Shanghai Astronomical Observatory,Chinese Academy of Sciences,200030 Shanghai,China)

This paper proposes a prediction method of UT1-UTC in Earth orientation parameters(EOP)by dual differential least-squares(LS)and autoregressive(AR)model.Firstly,leap seconds are removed in UT1-UTC observations,and Earth zonal harmonic tidal are corrected.Then,the corrected UT1-UTC are processed by dual differential method,the stationarity of polar motion parameters is improved.LS+AR method is utilized to analyze the dual differential UT1-UTC to obtain the preliminary prediction results.Finally,the preliminary prediction results are processed by inverse dual differential method,and tidal correction are extrapolated and leap seconds are recovered to obtain high accuracy UT1-UTC prediction results.The prediction results are compared with EOP prediction comparison campaign(EOP_PCC)results.It shows that the short-term UT1-UTC parameters prediction error is at the same level of EOP_PCC.The one day prediction accuracy of UT1-UTC is less than 0.03ms,which is better than EOP_PCC one day UT1-UTC prediction accuracy.And,the daily routine UT1-UTC prediction in Beijing Aerospace Control Center is introduced.

Earth orientation parameters;UT1-UTC prediction;dual differential;AR model;earth zonal harmonic tidal correction

P228

:A

:2095-7777(2014)03-0230-06

10.15982/j.issn.2095-7777.2014.03.012

陳略(1983—),男,工程師,主要研究方向:地球自轉(zhuǎn)參數(shù)高精度解算與預(yù)報(bào),無線電測量與科學(xué)應(yīng)用。

2014-07-28

2014-08-19

國家自然科學(xué)基金(41304026);國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)項(xiàng)目(2015CB857101)