兩種線(xiàn)性地球自轉(zhuǎn)參數(shù)短期預(yù)報(bào)方法

郭忠臣, 徐 波

(1. 宿州學(xué)院 環(huán)境與測(cè)繪工程學(xué)院, 安徽 宿州 234000; 2. 河南省測(cè)繪工程院, 河南 鄭州 450003)

兩種線(xiàn)性地球自轉(zhuǎn)參數(shù)短期預(yù)報(bào)方法

郭忠臣1, 徐 波2

(1. 宿州學(xué)院 環(huán)境與測(cè)繪工程學(xué)院, 安徽 宿州 234000; 2. 河南省測(cè)繪工程院, 河南 鄭州 450003)

介紹了兩種常用線(xiàn)性地球自轉(zhuǎn)參數(shù)短期預(yù)報(bào)的原理和方法,處理了地球自轉(zhuǎn)參數(shù)相關(guān)數(shù)據(jù),并進(jìn)行了分析.

地球自轉(zhuǎn)參數(shù); 最小二乘; ARMA; 定階

地球自轉(zhuǎn)參數(shù)(ERP)在高精度衛(wèi)星導(dǎo)航、深空探測(cè)等領(lǐng)域具有不可替代的作用,對(duì)其實(shí)時(shí)性的要求也越來(lái)越高[1-2].當(dāng)前高精度ERP的觀測(cè)手段主要有VLBI、SLR、GNSS等,但這些方法并不能實(shí)時(shí)地給出ERP,根據(jù)測(cè)定技術(shù)的不同,往往都需要延遲幾個(gè)小時(shí)至數(shù)天,而且觀測(cè)的精度也不甚相同,如GNSS技術(shù)雖然解算速度較SLR、VLBI等技術(shù)快,但是其解算結(jié)果有時(shí)還需VLBI和SLR技術(shù)的解算值來(lái)進(jìn)行修正[2-4],因此想要實(shí)時(shí)測(cè)定高精度的ERP目前并不可行.基于此,高精度的ERP預(yù)測(cè)技術(shù)就顯得尤為重要,但因ERP受到眾多可測(cè)及不可測(cè)的激發(fā)源的影響,導(dǎo)致其預(yù)報(bào)工作較為困難.因此,構(gòu)建一種較為合理、有效的高精度預(yù)報(bào)模型也成為了當(dāng)下大地測(cè)量工作者的主要任務(wù)之一[5].本文對(duì)當(dāng)前主要使用的預(yù)報(bào)方法原理進(jìn)行分析,并依據(jù)ERP相關(guān)性質(zhì),構(gòu)造出相應(yīng)的預(yù)報(bào)模型,以供后續(xù)研究使用.

1 ERP短期預(yù)報(bào)方法

為推進(jìn)地球自轉(zhuǎn)參數(shù)預(yù)測(cè)的研究,IERS于2005年10月1日—2008年2月28日組織了地球自轉(zhuǎn)參數(shù)預(yù)測(cè)方案比較大會(huì)戰(zhàn)(earth orientation parameters prediction comparison campaign, EOP PCC)來(lái)確定EOP預(yù)測(cè)的現(xiàn)狀,活動(dòng)在兩年半的時(shí)間內(nèi)收到了11位參與人員的近6500份預(yù)報(bào)值,采用的預(yù)報(bào)方法主要有:最小二乘外推法和自回歸模型(lst squares extrapolation of the harmonic model and autoregressive prediction)、譜分析和最小二乘外推法(sectral analysis and least squares extrapolation)、卡爾曼濾波法(Kalman filter)、神經(jīng)網(wǎng)絡(luò)模型(nural networks)、聯(lián)合大氣角動(dòng)量預(yù)報(bào)(ading AAM in LOD prediction)、小波分解和自協(xié)方差預(yù)測(cè)(wvelet decomposition an auto-ovariance prediction)等方法[6-7].此次活動(dòng)得出:不同方法組合的預(yù)報(bào)精度優(yōu)于單一方法,沒(méi)有一種預(yù)報(bào)方法既適合所有跨度又適合EOP所有參數(shù).國(guó)內(nèi)外學(xué)者一致認(rèn)為,在極移預(yù)報(bào)方面:LS+AR模型、頻譜分析與最小二乘方法及神經(jīng)網(wǎng)絡(luò)較優(yōu)于其他方法;在UT1-TC/LOD方面,小波分析和自協(xié)方差、卡爾曼濾波、聯(lián)合AAM預(yù)報(bào)優(yōu)于其他預(yù)報(bào)方法.

當(dāng)前對(duì)ERP進(jìn)行預(yù)報(bào)的主要方法可分為線(xiàn)性方法和非線(xiàn)性方法,本文主要對(duì)兩種常用線(xiàn)性方法的預(yù)報(bào)原理進(jìn)行介紹.

2 最小二乘外推法

最小二乘法是通過(guò)最小化誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配,利用最小二乘法可以簡(jiǎn)單地求出未知參數(shù)[8].早期對(duì)ERP進(jìn)行預(yù)測(cè)時(shí),多采用最小二乘法對(duì)ERP的趨勢(shì)項(xiàng)和周期項(xiàng)進(jìn)行擬合,并通過(guò)擬合參數(shù)對(duì)ERP外推預(yù)報(bào).

最小二乘法的一般模型可表示為[9-10]:

式中:L為觀測(cè)向量;B為觀測(cè)矩陣.

將式(1)轉(zhuǎn)換成誤差方程即:

式中:i=1,2,…,n,令

則式(2)可化為

根據(jù)最小二乘原理

式(4)可表達(dá)為

通過(guò)式(5)即可求得:

由文獻(xiàn)[8]可知,極移受長(zhǎng)期趨勢(shì)項(xiàng)和周期項(xiàng)(Chandler項(xiàng)、周年項(xiàng))的影響比較大,因此,根據(jù)最小二乘原理,在只考慮趨勢(shì)項(xiàng)和周期項(xiàng)的條件下可構(gòu)造極移PMX和PMY的最小二乘模型如下:

式中ai、bi為趨勢(shì)項(xiàng)的擬合系數(shù),ci、di為各周期項(xiàng)的振幅,P1、P2分別為周年項(xiàng)和Chandler項(xiàng)的周期,ωi為各周期項(xiàng)的相位,t為時(shí)間.

與極移類(lèi)似,LODR也主要受到長(zhǎng)期趨勢(shì)項(xiàng)和周期項(xiàng)(周年項(xiàng)、半周年項(xiàng))的影響,可構(gòu)造其最小二乘模型如下:

式中a、b為趨勢(shì)項(xiàng)的擬合系數(shù);c、d為各周期項(xiàng)的振幅;P1、P2分別為周年項(xiàng)和半周年項(xiàng)的周期;ωi為各周期項(xiàng)的相位;t為時(shí)間.

圖1為用最小二乘方法對(duì)2005年1月1日—2015年1月1日(MJD:53371-57023)期間的地球自轉(zhuǎn)參數(shù)(PMX、PMY和LODR)的擬合結(jié)果.

圖1 地球自轉(zhuǎn)參數(shù)擬合結(jié)果Fig.1 The fitting result of ERP(a)—PMX; (b)—PMY; (c)—LODR.

最小二乘法只能對(duì)ERP中的周期項(xiàng)和趨勢(shì)項(xiàng)等主項(xiàng)進(jìn)行外推預(yù)報(bào),由于該模型并不能將各影響因素完整的反映出來(lái),所以通過(guò)該式擬合后會(huì)得到一組殘差序列,該殘差序列包含了其他較小因素的影響程度,需通過(guò)其他模型對(duì)該殘差項(xiàng)進(jìn)行處理預(yù)報(bào).

3 ARMA模型

3.1 模型簡(jiǎn)介

對(duì)包含測(cè)量隨機(jī)噪聲以及不甚明顯的周期性變化的時(shí)間序列進(jìn)行分析,并建立其數(shù)學(xué)模型的方法和過(guò)程稱(chēng)為時(shí)間序列分析,簡(jiǎn)稱(chēng)時(shí)序分析[11-12].時(shí)間序列分析的ARMA模型就是對(duì)平穩(wěn)時(shí)間序列Zt(t=1,2,…,n)建立理想的統(tǒng)計(jì)模型,其數(shù)學(xué)模型表示為[13]

式(9)被稱(chēng)為p階自回歸,q階滑動(dòng)平均的時(shí)間序列模型,簡(jiǎn)記為ARMA(p,q),其中:φ1,φ2,…,φp為模型的自回歸系數(shù),θ1,θ2,…,θq為模型的滑動(dòng)平均系數(shù),at為t=t,t-1,…,t-q時(shí)刻的白噪聲序列.

由式(9)可知,zt主要受序列內(nèi)部規(guī)律和白噪聲的影響,即式(9)可表示為

若zt僅與t時(shí)刻的白噪聲以及之前的規(guī)律性有關(guān),則式(10)可表示為

此時(shí),稱(chēng)式(11)為p階自回歸模型,記為ARMA(p,0)或者AR(p).

若zt僅與t時(shí)刻以及之前的白噪聲有關(guān),此時(shí)式(10)可表示為

這時(shí),稱(chēng)式(12)為q階滑動(dòng)平均模型,記為ARMA(0,q)或者M(jìn)A(q).

假設(shè)一時(shí)間序列zt(t=1,2,…,n)為AR(p)模型,即有

令t=n+1,即可根據(jù)式(13)得到zn+1時(shí)刻的狀態(tài):

3.2 模型識(shí)別

使用ARMA模型對(duì)時(shí)間序列進(jìn)行分析時(shí),需選擇與該序列契合度較高的模型,通常根據(jù)時(shí)間序列的自/偏自相關(guān)函數(shù)的性質(zhì)來(lái)判斷.

設(shè)有平穩(wěn)時(shí)間序列Zt(t=1,2,…,n),其自/偏自相關(guān)函數(shù)計(jì)算公式如下:

通過(guò)證明可知,平穩(wěn)時(shí)間序列的相關(guān)函數(shù)存在“截尾”和“拖尾”的特性.“截尾”是指當(dāng)k大于某個(gè)值時(shí),相關(guān)函數(shù)全部為“零”,“拖尾”是指當(dāng)k大于某個(gè)值時(shí),相關(guān)函數(shù)趨近于“零”,但不會(huì)恒等于“零”.時(shí)間序列相關(guān)函數(shù)的“截尾”和“拖尾”現(xiàn)象,為模型選擇提供了依據(jù),表1為模型選擇原則.

表1 平穩(wěn)時(shí)間序列ARMA模型識(shí)別原則

根據(jù)表1中模型識(shí)別原則,分別求出ERP各項(xiàng)殘差序列的自/偏自相關(guān)函數(shù),通過(guò)分析其“截尾”“拖尾”特性即可判斷使用何種模型.因兩個(gè)極移方向的變化相同,LOD與UT1-UTC又可互相轉(zhuǎn)化得到,因此本文只對(duì)PMX和LOD殘差序列的相關(guān)函數(shù)進(jìn)行分析.通過(guò)圖2和圖3可以看出,PMX和LOD的自相關(guān)函數(shù)均存在“拖尾”現(xiàn)象,偏自相關(guān)函數(shù)均存在“截尾”現(xiàn)象,依據(jù)表1中模型判別原則,即可知PMX、LOD的殘差序列均適用于AR模型,又因PMY和UT1-UTC分別與PMX和LOD類(lèi)似,所以PMY、UT1-UTC的殘差序列也均適用于AR模型.

圖2 PMX殘差序列相關(guān)函數(shù)圖Fig.2 The correlation function figure of PMX residual

圖3 LOD殘差序列相關(guān)函數(shù)圖Fig.3 The correlation function figure of LOD residual

3.3 模型定階

確定模型類(lèi)別之后還需確定模型階數(shù),通常情況下,模型階數(shù)越高,擬合精度越好,但隨著階數(shù)的增加,需要的信息也隨之增多,當(dāng)信息量不變的情況下,階數(shù)過(guò)多反而會(huì)導(dǎo)致擬合誤差變大,而階數(shù)過(guò)少,預(yù)測(cè)結(jié)果會(huì)容易受到觀測(cè)粗差的影響,因此需要一種合理的方法來(lái)確定模型階數(shù)[12].常用的AR模型定階準(zhǔn)則有信息論準(zhǔn)則(AIC)、傳遞函數(shù)準(zhǔn)則(CAT)和最終預(yù)測(cè)誤差準(zhǔn)則(FPE),在使用時(shí),這三個(gè)準(zhǔn)則是等效的,本文主要介紹AIC準(zhǔn)則來(lái)確定AR模型階數(shù)p.

AIC準(zhǔn)則由赤池弘次于1973年提出,其一般形式為

對(duì)于AR模型,假設(shè)其階數(shù)為k,AIC準(zhǔn)則可表示為

根據(jù)AIC準(zhǔn)則判定AR模型最優(yōu)階數(shù)時(shí)首先需合理的確定階數(shù)上限p(p?n),理論上當(dāng)AIC(k)取最小值時(shí),階數(shù)k即為最佳階數(shù),但因該準(zhǔn)則只是模型優(yōu)化的一種宏觀度量,并不能單純的以AIC(k)最小來(lái)確定階數(shù),根據(jù)分析,當(dāng)階數(shù)大于一定值后,AIC(p)的變化幅度較小.選取一組數(shù)據(jù)進(jìn)行試驗(yàn),考慮到基礎(chǔ)時(shí)間序列的長(zhǎng)度以及AIC(p)的變化情況,本次試驗(yàn)取階數(shù)上限為30,當(dāng)序列長(zhǎng)度或跨度時(shí)間不同時(shí),模型最優(yōu)階數(shù)的選取可能不同,隨機(jī)選取多種情況中的一種,分別計(jì)算不同階數(shù)對(duì)應(yīng)的AIC值.通過(guò)圖4可知,當(dāng)階數(shù)為28時(shí),此時(shí)AIC值最小,但當(dāng)階數(shù)大于8時(shí),AIC值的變化趨于平緩,故在確定AR模型階數(shù)時(shí),可根據(jù)基礎(chǔ)序列的長(zhǎng)度以及計(jì)算時(shí)間效率等方面綜合考慮,以取得較優(yōu)的階數(shù).

圖4 AIC準(zhǔn)則Fig.4 AIC criterion

假設(shè)有零均值平穩(wěn)時(shí)間序列Zt(t=1,2,…,n),確定其最佳階數(shù)為p,則該時(shí)間序列可表示為

此時(shí)有

即

又可表示為

式(21)用矩陣表示為

4 結(jié) 語(yǔ)

當(dāng)前有多種空間大地測(cè)量技術(shù)可精確測(cè)定ERP,但由于數(shù)據(jù)處理過(guò)程復(fù)雜,并不能實(shí)時(shí)得到高精度的ERP測(cè)定結(jié)果.考慮到ERP在現(xiàn)代空間導(dǎo)航中的作用日漸顯著,高精度ERP預(yù)報(bào)方法也成為了當(dāng)今亟待解決的問(wèn)題之一.本文主要對(duì)兩種常用線(xiàn)性預(yù)報(bào)方法的原理進(jìn)行介紹,并通過(guò)實(shí)驗(yàn)對(duì)其在ERP預(yù)報(bào)中的可行性進(jìn)行分析,也為后續(xù)使用該方法進(jìn)行預(yù)報(bào)提供基礎(chǔ).

[ 1 ] WANG Q X,DANG Y M,XU TI H. The method of earth rotation parameter determination using GNSS observations and precision analysis[C]∥Lecture Notes in Electrical Engineering. Berlin: Springer, 2013:247-256.

[ 2 ] YAO Y B,YUE S Q,PENG C. A new LS+AR model with additional error correction for polar motion forecast[J]. Science China Earth Science, 2013,56(5):818-828.

[ 3 ] 魏二虎,常亮,劉經(jīng)南. 我國(guó)進(jìn)行激光測(cè)月的研究[J]. 測(cè)繪信息與工程, 2006,31(3):1-3.

WEI E H,CHANG L,LIU J N. Research on the lunar laser ranging for China[J]. Journal of Geomatics, 2006,31(3).

[ 4 ] 張志斌,王廣利,劉祥,等. 中國(guó)VLBI網(wǎng)觀測(cè)地球定向參數(shù)能力分析[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2013,38(8):911-915.

ZHANG Z B,WANG G L,LIU X,et al. Analysis of EOPdetermination via Chinese VLBI network[J]. Geomatics amp; Information Science of Wuhan University, 2013,38(8):911-915.

[ 5 ] 王潛心. 基于GNSS觀測(cè)數(shù)據(jù)的高精度地球自轉(zhuǎn)參數(shù)測(cè)定理論與方法研究[D]. 西安:西安測(cè)繪研究所, 2015.

WANG Q X. The theory and method of earth rotation parameters determination using GNSS observation data[D]. Xi’an: Xi’an Institute of Surveying and mapping, 2015.

[ 6 ] KALARUS M,SCHUH H,KOSEK W,et al. Achievements of the earth orientation parameters prediction comparison campaign[J]. Journal of Geodesy, 2010,84(10):587-596.

[ 7 ] XU X Q,ZHOU Y H . EOP prediction using least square fitting and autoregressive filter over optimized data intervals[J]. Advances in Space Research, 2015,56(10):2248-2253.

[ 8 ] 賈小勇,徐傳勝,白欣. 最小二乘法的創(chuàng)立及其思想方法[J]. 西北大學(xué)學(xué)報(bào)(自然科學(xué)版), 2006,36(3):507-511.

JIA X Y,XU C S,BAI X. The invention and way of thinking on least sqeares[J]. Jounnal of Northwest University (Natural Science Edition), 2006,36(3):507-511.

[ 9 ] XU X Q,ZHOU Y H. EOP prediction using least square fitting and autoregressive filter over optimized data intervals[J]. Advances in Space Research, 2015,56(10):2248-2253.

[10] XU X Q,ZHOU Y H,LIAO X H. Short-term earth orientation parameters predictions by combination of the least-squares, AR model and Kalman filter[J]. Journal of Geodynamics, 2012,62(8):83-86.

[11] 武偉,劉希玉,楊怡,等. 時(shí)間序列分析方法及ARMA,GARCH兩種常用模型[J]. 計(jì)算機(jī)技術(shù)與發(fā)展, 2010,20(1):247-249.

WU W,LIU X Y,YANG Y,et al. Analysis method of time array and two mosels of ARMA and GARCH[J]. Computer Technology and Development, 2010,20(1):247-249.

[12] 陳曉鋒. AIC準(zhǔn)則及其在計(jì)量經(jīng)濟(jì)學(xué)中的應(yīng)用研究[D]. 天津:天津財(cái)經(jīng)大學(xué), 2012.

CHEN X F. Akaike information criterion and its application in econometrics[D]. Tianjin: Tianjin University of Finance amp; Economics, 2012.

[13] 劉海洋,郝哲. 基于時(shí)序分析的邊坡變形預(yù)報(bào)與變形行為特征[J]. 沈陽(yáng)大學(xué)學(xué)報(bào)(自然科學(xué)版), 2012,24(2):81-86.

LIU H Y,HAO Z. Prediction and behavior characteristic resarach of slope deformation based on time series analysis[J]. Journal of Shenyang University(Natutral Science), 2012,24(2):81-86.

【責(zé)任編輯:肖景魁】

TwoLinearShort-TermPredictionMethodsofEarthRotationParaments

GuoZhongchen1,XuBo2

(1. Faculty of Environment Science and Surveying Engineering, Suzhou University, Suzhou 234000, China; 2. Henan Surveying and Mapping Engineering Institute, Zhengzhou 450003, China)

The principle of two short-time forecasting method of linear earth rotation parameters was introduced, and the data of the Earth’s rotation parameter was analyzed.

earth rotation parameters; least square; ARMA; order

P 228.6

A

2017-09-25

郭忠臣(1992-),男,安徽阜陽(yáng)人,宿州學(xué)院助教.

2095-5456(2017)06-0505-06