基于高斯過程的日長變化預報?

雷 雨123 趙丹寧13 高玉平12 蔡宏兵12

(1中國科學院國家授時中心西安710600)

(2中國科學院時間頻率基準重點實驗室西安710600)

(3中國科學院大學北京100049)

基于高斯過程的日長變化預報?

雷 雨1,2,3?趙丹寧1,3 高玉平1,2 蔡宏兵1,2

(1中國科學院國家授時中心西安710600)

(2中國科學院時間頻率基準重點實驗室西安710600)

(3中國科學院大學北京100049)

由于日長(length-of-day,LOD)變化具有復雜的時變特性,傳統(tǒng)線性模型如最小二乘外推模型、時間序列分析模型等的預報效果往往不甚理想,所以將一種新型的機器學習算法—高斯過程(Gaussian processes,GP)方法用于LOD變化預報,并將預報結(jié)果同利用反向傳播神經(jīng)網(wǎng)絡(back propagation neural networks,BPNN)和廣義回歸神經(jīng)網(wǎng)絡(general regression neural networks,GRNN)的預報結(jié)果以及地球定向參數(shù)預報比較競賽(Earth Orientation Parameters Prediction Comparison Campaign,EOP PCC)的預報結(jié)果進行對比.結(jié)果表明,GP用于LOD變化預報是高效可行的.

天體測量,時間,方法:數(shù)據(jù)分析

1 引言

LOD變化是表征地球自轉(zhuǎn)變化的一個重要參數(shù),它是指天文意義上的1 d和標準日長86 400 s之間的差異,反映了地球自轉(zhuǎn)速率的變化.LOD和極移(polar motion, PM)統(tǒng)稱為地球自轉(zhuǎn)參數(shù)(Earth rotation parameters,ERP).ERP是實現(xiàn)天球參考系和地球參考系之間相互轉(zhuǎn)換的必需參數(shù),在深空探測、衛(wèi)星精密定軌和天文地球動力研究等領域都有重要應用[1].現(xiàn)代測地技術(shù)(甚長干涉基線(Very Long Baseline Interferometry)、全球衛(wèi)星導航系統(tǒng)(Global Navigation Satellite Systems,GNSS)和衛(wèi)星激光測距(Satellite Laser Ranging,SLR)等)被廣泛應用于地球自轉(zhuǎn)變化的常規(guī)監(jiān)測中,提供了高時空分辨率和高精度的觀測資料.然而,由于復雜的數(shù)據(jù)處理過程,由現(xiàn)代測地技術(shù)獲取的ERP往往需要延遲幾天甚至2個星期,所以對ERP進行實時快速的預報成為一項值得深入研究的課題.

LOD變化的精確預報是ERP預報中的難點之一,特別是在厄爾尼諾(El Ni?no)事件發(fā)生期間,熱帶季風的變化導致LOD變化出現(xiàn)大幅振蕩.LOD變化的高精度實時快速預報引起了越來越多學者的關(guān)注.LOD變化主要由潮汐項和非潮汐項2部分組成,潮汐項可以由國際地球自轉(zhuǎn)與參考系服務(International Earth Rotation and Reference Systems Service,IERS)協(xié)議給出的模型精確確定[2],而非潮汐項中的半年項和周年項等季節(jié)性變化主要是由固體地球和全球大氣、海洋以及地下水之間的角動量交換引起的[3].

學者們在ERP預報方面已經(jīng)做了許多研究,提出了各種預報模型,包括最小二乘(least squares,LS)外推模型[4]、LS外推模型和自回歸(autoregressive,AR)模型的組合(LS+AR)[4?5]、卡爾曼濾波聯(lián)合大氣角動量(Kalman Filter+OAM)[6?7]、LS外推模型和人工神經(jīng)網(wǎng)絡(arti fi cial neural networks,ANN)模型的組合(LS+ANN)[3,8?11]、模糊推理系統(tǒng)(fuzzy-inference systems,FIS)[12]以及離散小波變換(discrete wavelet transform,DWT)和自協(xié)方差(autocovariance,AC)模型的組合(DWT+AC)[13]等.為了對比不同模型的預報效果,維也納理工大學大地測量與地球物理研究所從2005年10月1日至2008年2月28日組織了全球性的地球定向參數(shù)預報比較競賽,2 yr多的預報結(jié)果表明,沒有一種模型既適合于ERP所有分量的預報又適合于所有跨度的預報[14].

受多種激發(fā)因素的影響,地球自轉(zhuǎn)變化呈現(xiàn)復雜的非線性不規(guī)則變化特性,因此采用非線性的預報方法對其進行預報在理論上更為合理[8?11].ANN是逼近復雜非線性函數(shù)的一種有效工具,所以有許多學者將其應用于ERP預報中,并取得了顯著的預報效果[3,5,8?11].但是ANN存在一些缺點,例如網(wǎng)絡拓撲結(jié)構(gòu)難以確定、訓練過程存在過學習現(xiàn)象、迭代過程易陷入局部最優(yōu)、收斂速度較慢.此外,ANN的優(yōu)化目標是基于經(jīng)驗的風險最小化,無法保證網(wǎng)絡的泛化能力[15].針對上述缺點很多學者提出了改進措施,例如將ANN和其他人工智能(arti fi cial intelligence,AI)算法結(jié)合進行網(wǎng)絡優(yōu)化,如遺傳算法與ANN的組合、粒子群算法與ANN的組合等,然而這些算法均存在一定的不足,仍處于不斷的嘗試研究階段.

GP是近年來發(fā)展起來的一種新型機器學習算法,它有著嚴格的統(tǒng)計學習理論基礎,對處理高維數(shù)、小樣本、非線性等復雜問題具有很好的適應性,且泛化能力強[16?17].與ANN和支持向量機相比,GP具有容易實現(xiàn)、超參數(shù)自適應獲取、非參數(shù)推斷靈活以及輸出具有概率意義等優(yōu)點[16?17].現(xiàn)已成為機器學習領域的研究熱點,并在許多領域得到了成功的應用[18?19].本文將GP算法用于LOD變化預報,研究表明:GP算法用于LOD變化預報是可行的,且預報效率和精度較高.

2 GP的基本原理及其用于LOD變化預報的方法

2.1 GP回歸的基本原理

高斯過程又稱正態(tài)隨機過程,其任意有限變量集合都有著聯(lián)合高斯分布的特性,即對于任意的變量x1,x2,···,xn與其對應的函數(shù)f(x1),f(x2),···,f(xn)的聯(lián)合概率分布服從n維高斯分布.高斯分布的全部統(tǒng)計特征完全由它的均值函數(shù)m(x)和協(xié)方差函數(shù)C(x,x′)來確定,一般記為f(x)～GP(m(x),C(x,x′)).

若給定訓練樣本集D={(xi,yi)|xi∈Rd,yi∈R,i=1,2,···,n},其中d為向量xi的維數(shù),則對于測試樣本輸入x?,GP模型的預測值為

其中k(x?)=[C(x?,x1),C(x?,x2),···,C(x?,xn)]為測試樣本輸入和訓練樣本輸入值之間的1×n階協(xié)方差矩陣,K是訓練樣本輸入值之間的n×n階協(xié)方差矩陣, Kij=C(xi,xj).(1)式表明,GP模型可以根據(jù)協(xié)方差函數(shù)和測試樣本進行預測.

協(xié)方差函數(shù)在GP回歸模型中起到關(guān)鍵作用,它表達了一種樣本間的相似性,對所要學習的函數(shù)提供了假設信息,協(xié)方差函數(shù)必須是半正函數(shù).常用的協(xié)方差函數(shù)為平方指數(shù)函數(shù)[17],即

2.2 基于GP回歸模型的LOD變化預報

2.2.1 數(shù)據(jù)預處理

本文所用的LOD變化數(shù)據(jù)來自IERS發(fā)布的EOP C04序列,采樣間隔為1 d.LOD變化序列中周期為5 d～18.6 yr的固體地球帶諧潮汐項可以通過IERS協(xié)議給出的經(jīng)驗模型精確確定[2],近周日和半周日海洋潮汐項不作修正,LOD變化的長期趨勢項、季節(jié)性變化的周年和半年項等根據(jù)下述線性模型確定[4]:

其中LODR表示經(jīng)過固體地球帶諧潮修正后的LOD變化序列,p1、p2、p3和p4分別表示半年項、周年項、9.3 yr項和18.6 yr項的周期,取p1=182.62 d、p2=365.24 d、p3=3396.732d、p4=6793.464d,t為協(xié)調(diào)世界時(Coordinated Universal Time,UTC),在擬合時單位轉(zhuǎn)換為d.a、b表示長期趨勢項的參數(shù),c1,1、c1,2表示半年項的參數(shù), c2,1、c2,2表示周年項的參數(shù),c3,1、c3,2表示9.3 yr項的參數(shù),c4,1、c4,2表示18.6 yr項的參數(shù),這10個未知參數(shù)通過最小二乘法求得.

經(jīng)過上述數(shù)據(jù)預處理后的剩余部分為含有非線性成分的殘差序列,主要包括海洋近周日、半日潮項以及不規(guī)則的短周期成分.圖1從上至下依次繪出了1990—2010年期間LOD變化的原始序列、帶諧潮序列、線性模型擬合序列以及殘差序列.本文采用GP回歸模型對殘差序列進行建模和預測,將線性模型的預報結(jié)果和殘差序列的預報結(jié)果相加即可獲得最終的LOD變化預報值.

2.2.2 建模和預報

GP模型的建模過程就是通過對樣本數(shù)據(jù)D的訓練,確定協(xié)方差函數(shù)的超參數(shù).超參數(shù)的選取方法主要有交叉檢驗法、貝葉斯推理法和最大似然法[13].本文采用最大似然法選取超參數(shù),即任意給定超參數(shù)的初值,采用共軛梯度優(yōu)化算法求取訓練樣本對數(shù)似然函數(shù)的最大值,從而得到似然函數(shù)最大值所對應的超參數(shù)即為最優(yōu)超參數(shù).對數(shù)似然函數(shù)的形式為

圖1 LOD變化的原始序列(a);帶諧潮項(b);線性模型擬合項(c);殘差項(d)Fig.1 The raw series(a);tidal terms(b); fi tting terms of linear model(c);and residual terms(d)of the LOD variations

除了協(xié)方差函數(shù)及其超參數(shù),樣本的輸入和輸出方式也非常重要.按以下方式構(gòu)建樣本的輸入和輸出:

在訓練階段,樣本的輸入和輸出方式為

其中{ξ(i),i=1,2,···,n},表示LOD變化序列經(jīng)數(shù)據(jù)預處理后的殘差序列,根據(jù)經(jīng)驗確定,本文取d=4.

在預報階段,預報跨度為k=1,2,···,d,d+1,···時樣本的輸入和輸出方式分別為

2.2.3 精度評定指標

采用均方根誤差(root mean square error,RMSE)和平均絕對誤差(mean absolute error,MAE)作為預報結(jié)果的精度評定指標,其計算公式分別為

其中i為預報跨度,N為預報期數(shù),分別表示第j期的第id LOD變化的預測值和實際值.

3 實驗分析

首先將1990年1月1日至1999年12月31日的LOD變化殘差序列用于GP模型的訓練,然后用訓練好的GP回歸模型對2000—2001年的LOD變化殘差序列進行1～10 d、15 d、20 d、25 d、30 d、60 d、90 d、120 d、150 d、···、360 d跨度的預報(與Schuh等[3]和張曉紅等[8?9]的預報時間段相同).圖2給出了基于GP方法的預報跨度為1 d的LOD變化殘差的預報曲線(a)和預報誤差(b)圖,圖2(a)中虛線和實線分別代表殘差預報值和觀測值.

圖2 跨度為1 d的LOD變化殘差的預報結(jié)果(a)和預報誤差(b)Fig.2 The prediction results of the residual of LOD variations(a)and the predicted errors(b)at the prediction horizon of 1 d

同時本文將基于GP方法的LOD變化預報結(jié)果同Schuh等[3]使用的反向傳播神經(jīng)網(wǎng)絡(back propagation neural networks,BPNN)和張曉紅等[8?9]使用的改進的BPNN以及廣義回歸神經(jīng)網(wǎng)絡(general regression neural networks,GRNN)的預報結(jié)果進行對比,預報RMSE值見表1.

從圖2(a)、(b)和表1可以看出,GP用于LOD變化預報是可行和有效的,隨著預報跨度的增大,預報精度有所降低.

表1 GP預報結(jié)果與BPNN預報結(jié)果(Schuh等[3]的預報結(jié)果)、改進的BPNN及GRNN預報結(jié)果(張曉紅等[8?9]的預報結(jié)果)的比較(單位:ms)Table 1 The comparison of the prediction results of GP with those of the BPNN (Schuh et al.[3]),the modi fi ed BPNN,and the GRNN(Zhang et al.[8?9])(unit:ms)

為了更加直觀地比較4種方法的預報精度,圖3繪出了不同跨度的預報精度.從圖中可以看出,在短期(1～30 d)預報中,除了當跨度為1～3 d時GP的預報精度略低于BPNN預報精度外,其它跨度的預報精度均高于另外3種方法的預報精度.對于中期(1～360 d)預報,GP的預報精度仍優(yōu)于BPNN的預報精度,但低于改進BPNN和GRNN的預報精度.在預報效率上,由于GP模型參數(shù)(協(xié)方差參數(shù))可以自適應獲取,而不必像ANN技術(shù)需要對訓練樣本反復訓練才能得到最優(yōu)網(wǎng)絡參數(shù),故訓練速度較ANN技術(shù)要快.此外,因為本文采用的GP預報模式只需建模一次便可實現(xiàn)LOD變化的多天連續(xù)預報,因此預報所用時間較少,對于1～360 d的連續(xù)預報,訓練時間和預報時間之和一般在20 min以內(nèi),保證了算法的實時性.

圖3 GP預報精度與BPNN、改進BPNN和GRNN預報精度對比.(a)短期(1～30 d)預報,(b)中期(1～360 d)預報Fig.3The comparison of the prediction accuracies of GP with those of the BPNN,the modi fi ed BPNN, and the GRNN.(a)The short-term(1～30 d)prediction,and(b)the medium-term(1～360 d)prediction

為了與EOP PCC的預報結(jié)果進行比較,選取了1990年1月1日至2005年9月30日的LOD變化數(shù)據(jù)作為基礎序列,預報2005年10月1日到2008年2月28日(與EOP PCC預報時間段相同)1～360 d跨度的LOD變化值,統(tǒng)計了預報結(jié)果的MAE,并與EOP PCC的結(jié)果進行了對比,對比結(jié)果見圖4～6.

圖4～6中不同顏色和不同形狀線條分別代表參與EOP PCC的不同團隊所得的預報誤差,參與此項競賽的團隊詳細情況參見文獻[14].在圖4～5中紅色實線代表Gross團隊的預報誤差,粉色實線代表Kalarus團隊的預報誤差,藍色虛線代表Akyilmaz團隊的預報誤差,藍色點劃線代表Kosek團隊的預報誤差,綠色實線、綠色虛線和綠色點劃線代表Zotov團隊的預報誤差,黃色實線代表Pasynok團隊的預報誤差,藍色實線代表Mendes Cerveira團隊的預報誤差,黑色實線代表本文預報誤差;在圖6中黑色方形線條代表Mendes Cerveira團隊的預報誤差,黑色三角形線條代表Kosek團隊的預報誤差,黑色五角星線條代表Gross等的預報誤差,黑色圓形線條代表本文預報誤差,其中超短期(1～10 d)預報精度較高的團隊是Gross和Kalarus團隊,短期(1～30 d)預報精度較高的團隊是Gross、Kalarus和Kosek團隊,而參與中期(1～360 d)預報競賽的只有Gross、Kosek以及Mendes Cerveira 3個團隊.

從圖4～6的比較中可以看出,對于1～4 d的預報,GP方法的預報精度低于排在第1位的Gross等和排在第2位的Kalarus等的預報精度,從第5 d開始,GP的預報精度優(yōu)于Kalarus等的預報精度,但仍低于排在第1位的Gross等的預報精度;對于短期(1～30 d)預報,GP的預報精度僅次于排在第1位的Gross等的預報精度,與并列排在第2位的Kalarus等和Kosek等的預報精度大致相當;對于中期預報,GP的預報效果則不如EOP PCC.

圖4 超短期(1～10 d)MAE對比Fig.4 The comparison of the predicted MAE for the ultra short-term(1～10 d)

圖5 短期(1～30 d)MAE對比Fig.5 The comparison of the predicted MAE for the short-term(1～30 d)

圖6 中期(1～360 d)MAE對比Fig.6 The comparison of the predicted MAE for the medium-term(1～360 d)

4 討論與總結(jié)

本文根據(jù)實驗驗證了GP用于LOD變化的可行性和有效性.與ANN方法相比,GP方法較容易實現(xiàn),并且它不需要太多的先驗信息,只需事先選擇適當?shù)膮f(xié)方差函數(shù),其超參數(shù)在訓練過程中便可以自適應地確定,從而可以避免預報的人為主觀性,提高預報結(jié)果的可信度.通過實例發(fā)現(xiàn),GP方法用于LOD變化預報可以取得較好的預報效果.通過與ANN預報結(jié)果以及與EOP PCC預報結(jié)果的比較發(fā)現(xiàn),GP方法的短期(1～30 d)預報精度較高,但中期(1～360 d)預報精度則不如ANN預報精度和EOP PCC預報精度,這可能是由本文所使用的遞推預報模式的誤差累積效應引起的,對此可以嘗試以下兩種方法對GP中期預報精度進行改進,一種是改進樣本輸入方式,如采用連續(xù)輸入方式或者跨度輸入方式[20?21],另外一種方法是對GP模型進行在線訓練,本文暫不對此進行討論.在預報效率方面,由于本文采用遞推預報模式進行預報,故只需一次模型計算便可以實現(xiàn)多步預報,極大提高了預報效率.對于跨度為1～360 d的預報,應用GP方法預報LOD變化只需要20 min左右的時間,而應用ANN則需要數(shù)小時,預報效率大大提高,這對于ERP的實時快速預報具有重要的現(xiàn)實意義.此外,基于GP方法預報LOD變化所需訓練樣本數(shù)量遠遠小于ANN所需樣本量,這在歷史數(shù)據(jù)較少的情況下進行LOD變化的預報具有更高的現(xiàn)實意義.

由于協(xié)方差函數(shù)類型、超參數(shù)選取方法以及樣本輸入方式等對GP方法的預報效果都有一定的影響,因此,如何從上述角度來優(yōu)化GP模型以進一步提高LOD變化的預報精度是我們下一步的研究重點.對此我們將另行文討論.

致謝 感謝IERS提供的LOD變化資料,對中國科學院上海天文臺鄭大偉研究員提供的幫助表示由衷的感謝!

[1]Gambis D,Luzum B.Metro,2011,48:165

[2]McCarthy D D,Petit G.ITN,2003,13:605

[3]Schuh H,Ulrich M,Egger D,et al.JGeod,2002,76:247

[4]Tomasz N,Kose W.JGeod,2008,82:83

[5]許雪晴,周永宏.飛行器測控學報,2010,29:70

[6]Freedman A P,Steppe J A,Dickey J O,et al.JGR,1994,99:6981

[7]Gross R S,Eubanks T M,Steppe J A,et al.JGeod,1998,72:215

[8]張曉紅,王琪潔,朱建軍,等.天文學報,2011,52:322

[9]Zhang X H,Wang Q J,Zhu J J,et al.ChA&A,2012,36:86

[10]王琪潔,廖德春,周永宏.科學通報,2007,52:1728

[11]Wang Q J,Liao D C,Zhou Y H.ChSBu,2008,53:969

[12]Akyilmaz O,Kutterer H.JGeod,2004,78:82

[13]Kosek W,Kalarus M,Johnson T J.ArtSa,2005,40:119

[14]Kalarus M,Schuh H,Kosek W,et al.JGeod,2010,84:587

[15]Samanta B,Al-Balushi K R,Al-Araimi S A.JAGI,2003,16:657

[16]Seeger M.IJNS,2004,14:69

[17]何志昆,劉光斌,趙曦晶,等.控制與決策,2013,28:1121

[18]Brahim-Belhouari S,Bermak A.Computational Statistics&Data Analysis,2004,47:705

[19]劉冬,張清華.測繪學報,2011,40:59

[20]張曉紅,王琪潔,朱建軍,等.中國科學院上海天文臺年刊,2011,32:147

[21]Akyilmaz O,Kutterer H,Shum C K,et al.Applied Soft Computing,2011,11:837

The Prediction of Length-of-day Variations Based on Gaussian Processes

LEI Yu1,2,3ZHAO Dan-ning1,3GAO Yu-ping1,2CAI Hong-bing1,2
(1 National Time Service Center,Chinese Academy of Sciences,Xi’an 710600)
(2 Key Laboratory of Time and Frequency Primary Standards,Chinese Academy of Sciences,Xi’an
710600)
(3 University of Chinese Academy of Sciences,Beijing 100049)

Due to the complicated time-varying characteristics of the length-of-day (LOD)variations,the accuracies of traditional strategies for the prediction of the LOD variations such as the least squares extrapolation model,the time-series analysis model, and so on,have not met the requirements for real-time and high-precision applications. In this paper,a new machine learning algorithm—the Gaussian process(GP)model is employed to forecast the LOD variations.Its prediction precisions are analyzed and compared with those of the back propagation neural networks(BPNN),general regression neural networks(GRNN)models,and the Earth Orientation Parameters Prediction Comparison Campaign(EOP PCC).The results demonstrate that the application of the GP model to the prediction of the LOD variations is efficient and feasible.

astrometry,time,methods:data analysis

P127;

A

10.15940/j.cnki.0001-5245.2015.01.007

2014-07-02收到原稿,2014-08-04收到修改稿?國家自然科學基金項目(10573019)資助

?leiyu@ntsc.ac.cn