雷雨,蔡宏兵,高玉平

(1.中國科學(xué)院 國家授時中心,西安 710600；2.中國科學(xué)院 時間頻率基準(zhǔn)重點實驗室,西安 710600)

0 引言

地球自轉(zhuǎn)運動可以用地球定向參數(shù)(Earth orientation parameters,EOP)來表征。EOP包括歲差、章動、UT1-UTC和極移的兩個分量xp和yp。ERP是實現(xiàn)地球參考坐標(biāo)系與天球參考坐標(biāo)系相互轉(zhuǎn)換的必需參數(shù),在天文地球動力學(xué)研究、深空探測及衛(wèi)星導(dǎo)航等領(lǐng)域有著重要應(yīng)用[1-2]。甚長干涉基線測量(very long baseline interferometry,VLBI)、多普勒衛(wèi)星定軌和無線電定位(Doppler orbitography and radiopositioning integrated by satellite,DORIS)及全球衛(wèi)星導(dǎo)航系統(tǒng)(Global Navigation Satellite System,GNSS)等空間測地技術(shù)是測量EOP的主要手段,其中UT1-UTC和極移的測量精度可以分別達(dá)到10 μs 和100 μas[3],但復(fù)雜的資料處理過程致使EOP的獲取存在一定的延遲。由于深空探測等領(lǐng)域?qū)OP實時測量值有重要需求,而衛(wèi)星自主導(dǎo)航又對ERP中長期預(yù)報值有重要需求,因此對ERP進(jìn)行高精度的短期與中長期預(yù)報是非常必要的。在EOP的幾個分量中,UT1-UTC是變化最快、最難預(yù)報的一個分量[4-5]。

目前有多種UT1-UTC預(yù)報方法,其中多數(shù)方法結(jié)合最小二乘(least squares,LS)外推和其他模型進(jìn)行UT1-UTC預(yù)報[6-13],這些方法的思路為,首先采用最小二乘外推模型來提取UT1-UTC序列中的周期性成分,并對其進(jìn)行外推,然后利用神經(jīng)網(wǎng)絡(luò)(neural network,NN)或者自回歸(autoregressive,AR)模型等對最小二乘擬合殘差進(jìn)行建模、預(yù)報,最后再將周期項和殘差項的外推值相加以得到UT1-UTC預(yù)報值。在實際應(yīng)用中發(fā)現(xiàn),在利用最小二乘外推模型對UT1-UTC觀測資料進(jìn)行擬合時,在擬合序列的兩端存在發(fā)散畸變現(xiàn)象,這種現(xiàn)象在數(shù)據(jù)處理中稱為端部效應(yīng)[14-15]。端部效應(yīng)使殘差項與周期項的預(yù)報值出現(xiàn)偏差,最終導(dǎo)致UT1-UTC預(yù)報值不準(zhǔn)確。本文針對UT1-UTC預(yù)報中最小二乘擬合出現(xiàn)的端部畸變現(xiàn)象,在采用最小二乘外推模型對UT1-UTC序列進(jìn)行擬合之前,先利用端部延拓方法對UT1-UTC觀測資料進(jìn)行數(shù)據(jù)延拓,即在UT1-UTC序列的兩端增加應(yīng)用時序分析方法延拓出的若干數(shù)據(jù)點,形成一個新序列,然后用新序列求解最小二乘外推模型系數(shù),最后再基于最小二乘外推模型對UT1-UTC序列中的周期項進(jìn)行外推,這樣就可以將最小二乘擬合出現(xiàn)的端部畸變現(xiàn)象移至模擬序列的兩端。數(shù)值分析表明,通過在UT1-UTC觀測序列端部增加統(tǒng)計延拓數(shù)據(jù),可以有效地抑制端部效應(yīng),從而改進(jìn)UT1-UTC預(yù)報效果。

1 預(yù)報方法

本節(jié)首先介紹最小二乘和AR組合模型預(yù)報UT1-UTC的原理,在此基礎(chǔ)之上,建立考慮最小二乘擬合端部效應(yīng)的UT1-UTC預(yù)報算法。

1.1 數(shù)據(jù)預(yù)處理

本文UT1-UTC觀測資料來源于國際地球自轉(zhuǎn)與參考系服務(wù)(International Earth Rotation and Refe-rence Systems Service,IERS)發(fā)布的EOP-08-C04序列,其中UT1-UTC數(shù)據(jù)的采樣間隔為1 d[16]。UT1-UTC觀測資料中含有閏秒及多種周期項、準(zhǔn)周期項,其中,對于周期為5 d～18.6 a的62個固體地球帶諧潮汐項應(yīng)用IERS協(xié)議(IERS Conventions)給出的經(jīng)驗公式予以扣除[17],扣除62個固體地球帶諧潮汐項后的UT1-UTC稱作UT1R-UTC,然后再去掉UT1R-UTC序列中閏秒,獲得UT1R-TAI序列。本文對UT1-UTC的預(yù)測實質(zhì)上是針對UT1R-TAI的預(yù)測。

1.2 最小二乘外推模型

UT1-UTC序列扣除掉62個固體地球帶諧潮汐項后,還含有長期趨勢項、周年項、半周年項等周期性變化成分,長期趨勢項與周年項、半周年項利用如下模型進(jìn)行擬合、外推：

fUT1R-TAI=a+bt+c1cos(2πt/T1)+d1sin(2πt/T1)+c2cos(2πt/T2)+d2sin(2πt/T2),

(1)

式(1)中,T1、T2分別表示半周年項和周年項的振蕩周期,取T1=182.62 d、T2=365.24 d[4],a、b代表趨勢項參數(shù),c1、d1代表半周年項參數(shù),c2、d2代表周年項參數(shù),這些未知參數(shù)可以根據(jù)最小二乘法求解。

1.3 自回歸模型

AR模型是對平穩(wěn)時間序列{zt,t=1,2,…,n}建立的一個概率統(tǒng)計模型,它根據(jù)變量自身的歷史變化規(guī)律來建立統(tǒng)計模型,其數(shù)學(xué)公式為

(2)

式(2)中,p為模型階數(shù),εt為白噪聲,φ1,φ2,…,φp為模型參數(shù),可以通過求解Yule-Walker方程來獲得[18-19]。

運用AR模型的關(guān)鍵之處在于選取模型階數(shù)p。已有的定階準(zhǔn)則包括傳遞函數(shù)準(zhǔn)則、預(yù)測誤差準(zhǔn)則及信息論準(zhǔn)則[11]。本文選用赤池信息量準(zhǔn)則(Akaike information criterion,AIC)來確定模型階數(shù)p,AIC準(zhǔn)則的目標(biāo)函數(shù)可以表示為

(3)

在模型階數(shù)p與模型參數(shù)φ1,φ2,…,φp確定以后,可以通過如下方式對時間序列作多步外推：

(4)

1.4 預(yù)報過程

考慮最小二乘擬合端部效應(yīng)的UT1-UTC預(yù)測方法和常規(guī)方法的差別之處在于,該方法在對UT1-UTC觀測序列建立趨勢項及周期項最小二乘外推模型之前,先應(yīng)用統(tǒng)計學(xué)方法在UT1R-TAI序列首部和尾部進(jìn)行端部數(shù)據(jù)延拓,以抑制最小二乘擬合端部畸變,預(yù)報過程如下：

①首先通過式(1)對UT1R-TAI序列作最小二乘擬合,建立周期項及趨勢項外推模型,然后利用AR模型對最小二乘擬合殘差序列進(jìn)行建模、預(yù)報,最后聯(lián)合最小二乘外推模型和AR模型(LS+AR)在UT1R-TAI序列首部和尾部分別外推適當(dāng)數(shù)量的數(shù)據(jù)點,這樣UT1R-TAI序列加上首部和尾部外推的數(shù)據(jù)點就構(gòu)成了一個新序列。

②利用新序列求解最小二乘外推模型系數(shù),即用新序列重新建立趨勢項及周期項最小二乘外推模型,然后再結(jié)合最小二乘外推模型和AR模型對UT1R-TAI序列作外推預(yù)報。

③將周期為5 d～18.6 a的62個固體地球帶諧潮汐項及閏秒恢復(fù)到UT1R-TAI序列的預(yù)測值中即可得到UT1-UTC預(yù)測值。

圖1為考慮最小二乘擬合端部效應(yīng)的UT1-UTC預(yù)測方法的流程圖。

圖1 考慮最小二乘擬合端部效應(yīng)的UT1-UTC預(yù)測流程

2 數(shù)值分析

選取2001-01-01至2016-06-01期間的UT1-UTC觀測序列進(jìn)行數(shù)值分析,其中2010-01-01至2016-06-01為預(yù)測期,建模數(shù)據(jù)長度為10 a,每隔7 d預(yù)報1次,總共進(jìn)行了286次預(yù)測。圖2(a)和(b)分別給出了2001-01-01至2016-06-01時段的UT1-UTC測量序列及其扣除閏秒和62個固體地球帶諧潮汐項的UT1R-TAI序列。

圖2 UT1-UTC觀測序列及其UT1R-TAI序列

為了驗證端部延拓方法對最小二乘擬合端部效應(yīng)的改善效果,首先比較端部延拓前后最小二乘擬合效果。圖3繪出了端部數(shù)據(jù)延拓前后2001-01-01至2010-12-31時段UT1R-TAI序列的最小二乘擬合序列,其中,端部延拓數(shù)據(jù)點數(shù)為200,即在首尾兩端各延拓100個數(shù)據(jù)點,ECLS(edge-effect corrected least squares)表示端部效應(yīng)修正的最小二乘擬合序列。為了更加直觀地、清楚地展示端部數(shù)據(jù)延拓方法對LS擬合端部畸變的改善效果,對圖3作局部放大處理,圖4(a)繪出了端部數(shù)據(jù)延拓前后最小二乘擬合序列首部前100個歷元的最小二乘擬合殘差數(shù)據(jù)點,圖4(b)繪出了端部數(shù)據(jù)延拓前后最小二乘擬合序列尾部最后100個歷元的最小二乘擬合殘差數(shù)據(jù)點。

圖3 UT1R-TAI的LS與ECLS擬合序列

從圖4可以發(fā)現(xiàn),與直接對原始UT1R-TAI序列進(jìn)行擬合相比,端部數(shù)據(jù)延拓后最小二乘擬合殘差序列在首部和尾部更加接近于零,換言之,端點延拓后最小二乘擬合的UT1R-TAI序列在首部、尾部和原始序列吻合得更好,這說明端點延拓方法能夠有效地抑制最小二乘擬合出現(xiàn)的端部畸變現(xiàn)象。

圖4 UT1R-TAI的LS與ECLS擬合的首尾兩端殘差序列

為了進(jìn)一步檢驗端部數(shù)據(jù)延拓方法對最小二乘擬合端部效應(yīng)的改善效果,分別利用LS+AR模型和ECLS+AR模型對UT1-UTC作1～360 d時長預(yù)報,采用平均絕對誤差(mean absolute error,MAE)作為預(yù)測性能評價指標(biāo),其計算公式為[20]

(5)

式(5)中,Pj和Oj分別為j點的UT1-UTC預(yù)測值及觀測值,i表示預(yù)報長度,σi表示預(yù)報長度為i時的UT1-UTC預(yù)測平均絕對誤差,M表示預(yù)報期數(shù),本文總共作了286期的預(yù)測,即M=286。

圖5給出了LS+AR方法和ECLS+AR方法的UT1-UTC預(yù)測平均絕對誤差對比圖,表1給出了LS+AR方法和ECLS+AR方法在不同預(yù)報跨度下的平均絕對誤差統(tǒng)計結(jié)果,其中將1～15 d跨度的預(yù)報稱為短期預(yù)報,將大于15 d跨度的預(yù)報稱為中長期預(yù)報。從圖5及表1可以發(fā)現(xiàn),對于1～15 d的短期預(yù)報,ECLS+AR模型的預(yù)報精度相對于常規(guī)LS+AR模型的預(yù)報精度并無改善,但從第15 d開始,ECLS+AR模型的預(yù)報精度明顯優(yōu)于常規(guī)LS+AR模型,精度最大提高了34%,且一直保持在15%以上,這說明與常規(guī)LS+AR模型比較而言,ECLS+AR模型對于UT1-UTC中長期預(yù)報具有更明顯的優(yōu)勢,同時也從側(cè)面反映出最小二乘擬合的端部畸變現(xiàn)象對UT1-UTC中長期預(yù)報的影響更大。

圖5 LS+AR模型與ECLS+AR模型的UT1-UTC預(yù)報MAE對比

表1 LS+AR模型與ECLS+AR模型的UT1-UTC預(yù)報MAE統(tǒng)計

3 結(jié)語

本文提出了一種顧及最小二乘擬合端部效應(yīng)的UT1-UTC預(yù)報方法,這種方法與基于常規(guī)最小二乘外推模型的UT1-UTC預(yù)報方法區(qū)別在于,該方法在進(jìn)行最小二乘擬合之前,首先在原始序列兩端增加統(tǒng)計延拓數(shù)據(jù),然后再對數(shù)據(jù)延拓后的新序列進(jìn)行最小二乘擬合,其目的為將最小二乘擬合存在的端部畸變搬移到新序列的兩端,從而抑制原始觀測序列的端部畸變。數(shù)值分析表明,通過在觀測資料的兩端增加用統(tǒng)計學(xué)方法延拓出的外推數(shù)據(jù)點,然后再進(jìn)行最小二乘擬合,能夠有效地抑制端部效應(yīng)的影響；與常規(guī)LS+AR預(yù)報模型相比,基于端部效應(yīng)改善的ECLS+AR模型的UT1-UTC短期預(yù)報精度沒有改善,但對于UT1-UTC中長期預(yù)報而言,ECLS+AR模型精度提升尤為明顯,因此ECLS+AR模型可用于UT1-UTC中長期預(yù)報。