嚴(yán) 鳳 姚宜斌

(武漢大學(xué)測(cè)繪學(xué)院，武漢 430079)

地球自轉(zhuǎn)參數(shù)短期預(yù)報(bào)方法及其實(shí)現(xiàn)*

嚴(yán) 鳳 姚宜斌

(武漢大學(xué)測(cè)繪學(xué)院，武漢 430079)

根據(jù)極移及日長(zhǎng)變化的特性，采用LS與AR模型相結(jié)合的方法，建立了適合地球自轉(zhuǎn)參數(shù)(ERP)的預(yù)報(bào)模型。為減少相鄰殘差存在的強(qiáng)相關(guān)性，先對(duì)殘差序列進(jìn)行差分預(yù)處理，然后利用AR模型進(jìn)行ERP的預(yù)報(bào)。對(duì)ERP進(jìn)行了1—10天的短期預(yù)報(bào)，其預(yù)報(bào)精度接近國(guó)際最好水平。

地球自轉(zhuǎn)參數(shù);三角多項(xiàng)式模型;自回歸模型;預(yù)報(bào)精度;平穩(wěn)性檢驗(yàn)

1 引言

地球自轉(zhuǎn)參數(shù)是指由極移、日長(zhǎng)變化以及歲差和章動(dòng)構(gòu)成的描述地球自轉(zhuǎn)變化的參數(shù)，簡(jiǎn)稱為ERP。ERP表達(dá)了地固參考系相對(duì)于瞬時(shí)天球參考系的運(yùn)動(dòng)，是實(shí)現(xiàn)天球坐標(biāo)系與地球坐標(biāo)系相互轉(zhuǎn)換的必要參數(shù)。極移分量表示天球歷書軸在地球參考系中的運(yùn)動(dòng)。世界時(shí)與協(xié)調(diào)時(shí)之差或者日長(zhǎng)變化則反映了地球自轉(zhuǎn)的不均勻性。地球定向參數(shù)包含有豐富的地球動(dòng)力學(xué)信息，而且還是導(dǎo)航與衛(wèi)星精密定軌中必需的參數(shù)［1］，對(duì)其進(jìn)行實(shí)時(shí)高精度的預(yù)報(bào)具有重要的科學(xué)意義和應(yīng)用價(jià)值。

近幾十年來，隨著甚長(zhǎng)基線干涉測(cè)量、人衛(wèi)激光測(cè)距和全球定位系統(tǒng)等現(xiàn)代測(cè)量技術(shù)的飛速發(fā)展，極大地推動(dòng)了地球自轉(zhuǎn)變化的研究。然而，由于測(cè)量模型和資料處理的復(fù)雜性，目前VLBI和SLR等技術(shù)獲取的ERP參數(shù)往往需要延遲2～5天，難以滿足導(dǎo)航與衛(wèi)星精密定軌對(duì)ERP參數(shù)的時(shí)效性需求。因此，尋求高精度ERP預(yù)報(bào)方法成為一項(xiàng)值得深入研究的課題。目前國(guó)際上常用預(yù)報(bào)方法有最小二乘外推法、協(xié)方差法、自回歸移動(dòng)平均法、神經(jīng)網(wǎng)絡(luò)法、小波分析、卡爾曼濾波等方法，其預(yù)報(bào)精度已經(jīng)達(dá)到一定水平。國(guó)內(nèi)對(duì)ERP參數(shù)進(jìn)行預(yù)報(bào)起步較晚，但隨著探月工程以及北斗二代衛(wèi)星導(dǎo)航系統(tǒng)等空間計(jì)劃的逐步拓展，提高ERP參數(shù)預(yù)報(bào)的精度和時(shí)效性已變得十分迫切，目前多家科研機(jī)構(gòu)對(duì)此展開了研究。1977年湯家豪教授［2］在AR模型基礎(chǔ)上首次提出門限序列自回歸模型;2007年王琪潔將大氣角動(dòng)量的影響引入地球自轉(zhuǎn)參數(shù)預(yù)報(bào)中，建立了基于神經(jīng)網(wǎng)絡(luò)技術(shù)的地球自轉(zhuǎn)變化預(yù)報(bào);2008年楊元喜等人［3］提出為了帶多適應(yīng)因子的自適應(yīng)抗差濾波的動(dòng)態(tài)數(shù)據(jù)預(yù)報(bào)方法，有效改善了動(dòng)態(tài)預(yù)報(bào)精度;2010年許雪晴［4］等人比較LS+AR模型和人工神經(jīng)網(wǎng)絡(luò)模型的預(yù)報(bào)精度，發(fā)現(xiàn)AR模型在短期的預(yù)報(bào)中精度較高，而ANN模型則在中長(zhǎng)期預(yù)報(bào)中有優(yōu)勢(shì)。為了研究不同技術(shù)和方法對(duì)ERP預(yù)報(bào)精度的影響，2005年10月—2008年2月，國(guó)際上舉行了一次大型的地球定向參數(shù)預(yù)報(bào)競(jìng)賽。這次比賽總共由13組來自不同國(guó)家及不同學(xué)院的專家參與，經(jīng)過兩年多預(yù)報(bào)與比較，發(fā)現(xiàn)采用LS+AR模型對(duì)ERP參數(shù)進(jìn)行短期預(yù)報(bào)的精度很高。因此本文采用LS+AR模型對(duì)地球自轉(zhuǎn)參數(shù)進(jìn)行短期預(yù)報(bào)，并與EOP PCC的結(jié)果進(jìn)行比較，以驗(yàn)證預(yù)報(bào)模型的有效性。

2 LS+AR模型

2.1 LS模型

地球自轉(zhuǎn)參數(shù)刻畫了極移、日長(zhǎng)變化、歲差、章動(dòng)等自然現(xiàn)象的變化，這些變化包括周期性、趨勢(shì)性等復(fù)雜過程，時(shí)間序列中主要的部分是周期項(xiàng)和趨勢(shì)項(xiàng)，本文采用三角多項(xiàng)式模型，并用最小二乘方法來擬合這兩項(xiàng)。極移參數(shù)計(jì)算公式為:

式中a、b為趨勢(shì)項(xiàng)的擬合系數(shù)，c、d為周期項(xiàng)的擬合系數(shù)，P1、P2分別為 Chandler項(xiàng)和周年項(xiàng)的周期［5］，ω1、ω2分別為Chandler項(xiàng)和周年項(xiàng)的初始相位。當(dāng)方程個(gè)數(shù)大于未知數(shù)個(gè)數(shù)時(shí)，可根據(jù)最小二乘原理求解參數(shù)。

UT1-UTC、Δlod計(jì)算公式為:

式中:a、b、c為趨勢(shì)項(xiàng)的擬合系數(shù)，d、e、f為周期項(xiàng)的擬合系數(shù)，P1、P2、P3分別為Chandler項(xiàng)、周年項(xiàng)和半周年項(xiàng)的周期［6］，φ1、φ2、φ3為這3個(gè)周期項(xiàng)的初始相位。

2.2 AR模型

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

剔除周期項(xiàng)和趨勢(shì)項(xiàng)后的殘差序列，沒有明顯的規(guī)律特征。經(jīng)分析，發(fā)現(xiàn)該時(shí)間序列的后一項(xiàng)與前幾項(xiàng)具有相關(guān)性，根據(jù)這一特點(diǎn)采用自回歸(AR)模型來擬合殘差序列。自回歸模型就是根據(jù)變量自身過去的規(guī)律來建立預(yù)測(cè)模型，在動(dòng)態(tài)數(shù)據(jù)處理中存在廣泛的應(yīng)用。本文采用的自回歸模型為:

式中:xt為模型變量;b1，b2，…，bp為模型的回歸系數(shù);εt為模型的隨機(jī)誤差;p為模型階數(shù)［7］。

自回歸模型是時(shí)間序列分析的主要模型之一。在建立自回歸模型時(shí)，難點(diǎn)在于模型階數(shù)的確定，而這個(gè)問題可以利用Matlab高階統(tǒng)計(jì)量工具箱中的函數(shù)來解決。確定階數(shù)后，即可用最小二乘方法估計(jì)出模型參數(shù)。Matlab高階統(tǒng)計(jì)工具箱也提供了相關(guān)函數(shù)求解模型參數(shù)，計(jì)算十分簡(jiǎn)便。

2)外推預(yù)報(bào)

當(dāng)回歸系數(shù)b1，b2，…，bp確定時(shí)，可根據(jù)方程進(jìn)行預(yù)報(bào)。

第一步預(yù)報(bào)值為:

第二步預(yù)報(bào)值為:

一般地，第l步預(yù)報(bào)值為:

每一個(gè)基礎(chǔ)時(shí)間序列確定一個(gè)模型階數(shù)，然后根據(jù)模型階數(shù)求出回歸系數(shù)，再用公式(4)即可進(jìn)行預(yù)報(bào)。

2.3 地球自轉(zhuǎn)參數(shù)的預(yù)報(bào)

極移參數(shù)預(yù)報(bào)流程如圖1，UT1-UTC的預(yù)報(bào)流程如圖2。

從www.iers.org下載eopc04.62-now文件，抽取極移參數(shù)序列、UT1-UTC序列和Δlod序列(圖3)，采用從1995年至今的數(shù)據(jù)進(jìn)行建模分析。

3 差分處理

圖1 極移參數(shù)的預(yù)報(bào)流程Fig.1 Prediction process of polar motion parameters

AR模型要求建模的時(shí)間序列是平穩(wěn)隨機(jī)序列，即滿足平穩(wěn)、正態(tài)、零均值的條件［8］。剔除了周期項(xiàng)和趨勢(shì)項(xiàng)影響的殘差序列，其平穩(wěn)性還不能滿足平穩(wěn)隨機(jī)序列的條件。本文提出先對(duì)殘差序列進(jìn)行一次差分，再對(duì)差分后的殘差序列建立AR模型的處理方法。一次差分后的殘差序列平穩(wěn)性大大增強(qiáng)，以極移數(shù)據(jù)為例，首先利用時(shí)間序列的散點(diǎn)圖來判斷殘差序列以及差分后殘差序列的平穩(wěn)性(圖4)，可以看出，差分前時(shí)間序列都繞其均值上下波動(dòng)，但差分后的時(shí)間序列波動(dòng)更趨于平穩(wěn)，更接近零均值。

計(jì)算以250天作為子時(shí)間序列的極移參數(shù)的殘差和差分殘差的均值統(tǒng)計(jì)量以及方差統(tǒng)計(jì)量，結(jié)果如圖5、6。從圖5、6可以看出差分后時(shí)間序列的均值統(tǒng)計(jì)量及方差統(tǒng)計(jì)量更接近零，并且每段子時(shí)間序列的均值統(tǒng)計(jì)量及方差統(tǒng)計(jì)量基本不變，說明差分后的時(shí)間序列平穩(wěn)性更強(qiáng)。

極移參數(shù)以6年觀測(cè)數(shù)據(jù)為基礎(chǔ)時(shí)間序列，采用對(duì)殘差序列進(jìn)行差分處理和不做差分處理兩種方案分別進(jìn)行解算，圖7為兩種方案解算精度的比較，從圖7可以看出進(jìn)行一次差分處理后，大大提高了預(yù)報(bào)精度。

Δlod可以與UT1-UTC相互轉(zhuǎn)換，UT1-UTC去掉秒跳后，再作差分，結(jié)果的相反數(shù)即為Δlod。因此，對(duì)Δlod的預(yù)報(bào)可以轉(zhuǎn)換為對(duì)UT1-UTC的預(yù)報(bào)，再將UT1-UTC的預(yù)報(bào)值換算為Δlod的值即可。

圖2 UT1-UTC的預(yù)報(bào)流程Fig.2 Flow chart of UT1-UTC prediction

圖3 極移參數(shù)、UT1-UTC、Δlod基礎(chǔ)序列Fig.3 Basic sequences of polar motion parameters、UT1-UTC and Δlod

4 精度檢驗(yàn)及結(jié)果分析

4.1 精度檢驗(yàn)

圖4 極移參數(shù)的殘差序列和差分后殘差序列Fig.4 Residual sequences of the polar motion parameters before and after the difference

圖5 殘差序列和差分后殘差序列的均值統(tǒng)計(jì)量Fig.5 Mean statistics of residual sequences before and after the difference

圖6 殘差序列和差分后殘差序列的方差統(tǒng)計(jì)量Fig.6 Variance statistics of residual sequences before and after the difference

為與國(guó)際結(jié)果進(jìn)行比較，用預(yù)報(bào)誤差絕對(duì)值的平均值(MAE)來評(píng)定預(yù)報(bào)精度，計(jì)算式為:

圖7 做差分與不做差分預(yù)報(bào)精度比較Fig.7 Comparison of prediction accuracy between doing and not doing the difference

式中，εi，n是預(yù)報(bào)值與它在IERS C04中相應(yīng)值的差值，N為預(yù)報(bào)次數(shù)，I為預(yù)報(bào)長(zhǎng)度。

本文采用與EOP PCC相同的預(yù)報(bào)間隔和期數(shù)，并預(yù)報(bào)相同時(shí)間點(diǎn)的ERP參數(shù)。EOP PCC從2005年10月1日起開始預(yù)報(bào)，每隔7天預(yù)報(bào)10天的值，總共預(yù)報(bào)了126期，取預(yù)報(bào)誤差絕對(duì)值的平均值作為預(yù)報(bào)精度(圖8(a)［9］)。本文分別采用2005年10月1日前3年、前6年、前9年的數(shù)據(jù)作為基礎(chǔ)序列，建立預(yù)報(bào)模型，其預(yù)報(bào)精度如圖8(b)。圖8反映了本文方法與EOP PCC所用方法的精度比較，其具體分析見表1。

表1 ERP絕對(duì)誤差(MAE)比較Tab.1 Comparison between ERP mean absolute prediction error

由表1可以看出Xp、Yp方向1～5天的預(yù)報(bào)精度高于EOP PCC絕大部分預(yù)報(bào)方法，6～10天預(yù)報(bào)精度達(dá)到EOP PCC領(lǐng)先水平。UT1-UTC 1～10天預(yù)報(bào)精度均達(dá)到EOP PCC領(lǐng)先水平。對(duì)于Δlod的預(yù)報(bào)，EOP PCC中大部分方法的預(yù)報(bào)精度都不高，精度較高的預(yù)報(bào)方法1～10天的精度為0.13～0.38 ms/day;本文直接由UT1-UTC的預(yù)報(bào)結(jié)果轉(zhuǎn)換而來的Δlod，其預(yù)報(bào)精度也相當(dāng)高，1～10天的預(yù)報(bào)精度為0.17～0.25 ms/day。綜上可見，本文采用三角多項(xiàng)式模型與自回歸模型相結(jié)合的方法，其預(yù)報(bào)精度已經(jīng)達(dá)到國(guó)際領(lǐng)先水平。

圖8 ERP預(yù)報(bào)精度Fig.8 ERP prediction accuracy

4.2 基礎(chǔ)序列長(zhǎng)度對(duì)預(yù)報(bào)精度的影響

從上一節(jié)可以看出，不同長(zhǎng)度的基礎(chǔ)序列其預(yù)報(bào)結(jié)果是不一樣的。目前，國(guó)內(nèi)外對(duì)ERP的預(yù)報(bào)絕大多數(shù)沒有考慮基礎(chǔ)序列長(zhǎng)度對(duì)ERP預(yù)報(bào)精度的影響。本文分別以3年、6年、9年、12年作為基礎(chǔ)序列長(zhǎng)度對(duì)地球自轉(zhuǎn)參數(shù)進(jìn)行預(yù)報(bào)，尋找基礎(chǔ)序列長(zhǎng)度與ERP預(yù)報(bào)精度的聯(lián)系(圖9)。圖9中反映，基礎(chǔ)序列長(zhǎng)度越長(zhǎng)，其預(yù)報(bào)精度反而越差。這主要是因?yàn)橄噜彅?shù)據(jù)之間具有明顯的相關(guān)性，而隨著時(shí)間間隔的增加，相關(guān)性逐漸下降。同時(shí)ERP基礎(chǔ)序列過長(zhǎng)的話，其數(shù)據(jù)的精度不一致，從而影響預(yù)報(bào)精度。因此，在建立預(yù)報(bào)模型時(shí)，基礎(chǔ)序列長(zhǎng)度不宜太長(zhǎng)，這樣一方面可以保證基礎(chǔ)序列精度的一致，提高預(yù)報(bào)精度;另一方面，有效縮短了預(yù)報(bào)時(shí)間。另外，6年接近Chandler周期與周年周期的最小公倍數(shù)［9］，選取6年作為基礎(chǔ)序列長(zhǎng)度，更能體現(xiàn)地球自轉(zhuǎn)參數(shù)的周期特征。

圖9 3、6、9、12年的預(yù)報(bào)精度Fig.9 Prediction accuracy of 3，6，9 and 12 years

5 結(jié)語(yǔ)

本文采用LS+AR模型對(duì)地球自轉(zhuǎn)參數(shù)進(jìn)行預(yù)報(bào)，并與EOP PCC的結(jié)果進(jìn)行比較，以驗(yàn)證預(yù)報(bào)模型的有效性。提出先對(duì)殘差序列進(jìn)行一次差分處理的方法，求一次差分后序列的平穩(wěn)性大大增強(qiáng)。本文通過LS+AR模型預(yù)報(bào)精度達(dá)到國(guó)際領(lǐng)先水平。國(guó)內(nèi)外對(duì)EOP的預(yù)報(bào)絕大多數(shù)沒有考慮基礎(chǔ)序列長(zhǎng)度對(duì)EOP預(yù)報(bào)精度的影響，本文通過試驗(yàn)發(fā)現(xiàn)基礎(chǔ)序列長(zhǎng)度越長(zhǎng)，其預(yù)報(bào)精度反而越差，以6年作為基礎(chǔ)序列長(zhǎng)度比較合適。

1 黨亞民，秘金鐘，成英燕.全球?qū)Ш叫l(wèi)星系統(tǒng)原理與應(yīng)用［M］.北京:測(cè)繪出版社，2007.(Dang Yamin，Mi Jinzhong and Cheng Yingyan.Principles and applications of global navigation satellite system［M］.Beijing:Chinese Surveying and Mapping Press，2007)

2 王琪潔.基于神經(jīng)網(wǎng)絡(luò)技術(shù)的地球自轉(zhuǎn)變化預(yù)報(bào)［D］.中國(guó)科學(xué)院上海天文臺(tái)，2007.(Wang Qijie.The earth rotation changes prediction based on neural network technology［D］.Shanghai Astronomical Observatory，2007)

3 Yang Y X and Cui X Q.Adaptively robust filter with multi adaptive factors［J］.Survey Review，2008，40(309):260 -270.

4 許雪晴，周永宏.地球定向參數(shù)高精度預(yù)報(bào)方法研究［J］.飛行控制學(xué)報(bào)，2010，29(2):70-76(Xu Xueqing and Zhou Yonghong.The research of the earth orientation parameters high accuracy prediction method［J］.Journal of Flight Control，2010，29(2):70-76)

5 錢昌夏.極移速度的周期分析［J］.中國(guó)科學(xué)院上海天文臺(tái)年刊，1995，16:35-40.(Qian Changxia.Periodic analysis of velocity of polar motion［J］.Annals of Shanghai Observatory Academia Sinica，1995，16:35-40)

6 韓永剛，李志安.地球自轉(zhuǎn)速率變化主要中短周期的分析［J］.地球物理學(xué)進(jìn)展，2002，17(2):349-352.(Han Yonggang and Li Zhian.Analysis of main middle and short periods in the variation of the Earth’s rotation［J］.Progress in geophysics，2002，17(2):349-352)

7 張善文，雷英杰，馮有前.MATLAB在時(shí)間序列分析中的應(yīng)用［M］.西安:西安電子科技大學(xué)出版社，2007.(Zhang Shanwen，Lei Yingjie and Feng Youqian.Application of MATLAB in time series analysis［M］.Xi’an:Xi’an University of Electronic Science and Technology Press，2007)

8 王新洲，等.高等測(cè)量平差［M］.北京:測(cè)繪出版社，2006.(Wang Xinzhou，et al.Advanced surveying adjustment［M］.Beijing:Chinese Surveying and Mapping Press，2006)

9 李征航，等.空間大地測(cè)量學(xué)［M］.武漢:武漢大學(xué)出版社，2010.(Li Zhenghang，et al.Space geodesy［M］.Wuhan:Wuhan University Press，2010)

SHORT-TERM PREDICTION METHODS AND REALIZATION OF EARTH ROTATION PARAMETERS

Yan Feng and Yao Yibin
(School of Geodesy and Geomatics，Wuhan University，Wuhan 430079)

According to the respective characteristics of polar motion and length of day，combining LS and AR model，a suitable prediction model for the Earth Rotation Parameter(ERP)was builded.In order to reduce the strong correlation，the difference before establishing AR model was down and ERP 1-10 days were predicted.The results show that the prediction accuracy has reached an international leading level.

the Earth rotation parameters(ERP);trigonometric polynomial model;auto-regressive model;prediction accuracy;stationary test

1671-5942(2012)04-0071-05

2012-01-24

國(guó)家自然科學(xué)基金創(chuàng)新研究群體項(xiàng)目(41021061);國(guó)家自然科學(xué)基金(41174012)

嚴(yán)鳳，女，碩士生，主要研究方向?yàn)榈厍蜃赞D(zhuǎn)參數(shù)短期預(yù)報(bào).E-mail:642304401@qq.com

P227;P207

A