基于高斯过程的日长变化预报∗

雷 雨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)变化具有复杂的时变特性,传统线性模型如最小二乘外推模型、时间序列分析模型等的预报效果往往不甚理想,所以将一种新型的机器学习算法—高斯过程(Gaussian processes,GP)方法用于LOD变化预报,并将预报结果同利用反向传播神经网络(back propagation neural networks,BPNN)和广义回归神经网络(general regression neural networks,GRNN)的预报结果以及地球定向参数预报比较竞赛(Earth Orientation Parameters Prediction Comparison Campaign,EOP PCC)的预报结果进行对比.结果表明,GP用于LOD变化预报是高效可行的.

天体测量,时间,方法:数据分析

1 引言

LOD变化是表征地球自转变化的一个重要参数,它是指天文意义上的1 d和标准日长86 400 s之间的差异,反映了地球自转速率的变化.LOD和极移(polar motion, PM)统称为地球自转参数(Earth rotation parameters,ERP).ERP是实现天球参考系和地球参考系之间相互转换的必需参数,在深空探测、卫星精密定轨和天文地球动力研究等领域都有重要应用[1].现代测地技术(甚长干涉基线(Very Long Baseline Interferometry)、全球卫星导航系统(Global Navigation Satellite Systems,GNSS)和卫星激光测距(Satellite Laser Ranging,SLR)等)被广泛应用于地球自转变化的常规监测中,提供了高时空分辨率和高精度的观测资料.然而,由于复杂的数据处理过程,由现代测地技术获取的ERP往往需要延迟几天甚至2个星期,所以对ERP进行实时快速的预报成为一项值得深入研究的课题.

LOD变化的精确预报是ERP预报中的难点之一,特别是在厄尔尼诺(El Ni˜no)事件发生期间,热带季风的变化导致LOD变化出现大幅振荡.LOD变化的高精度实时快速预报引起了越来越多学者的关注.LOD变化主要由潮汐项和非潮汐项2部分组成,潮汐项可以由国际地球自转与参考系服务(International Earth Rotation and Reference Systems Service,IERS)协议给出的模型精确确定[2],而非潮汐项中的半年项和周年项等季节性变化主要是由固体地球和全球大气、海洋以及地下水之间的角动量交换引起的[3].

学者们在ERP预报方面已经做了许多研究,提出了各种预报模型,包括最小二乘(least squares,LS)外推模型[4]、LS外推模型和自回归(autoregressive,AR)模型的组合(LS+AR)[4−5]、卡尔曼滤波联合大气角动量(Kalman Filter+OAM)[6−7]、LS外推模型和人工神经网络(arti fi cial neural networks,ANN)模型的组合(LS+ANN)[3,8−11]、模糊推理系统(fuzzy-inference systems,FIS)[12]以及离散小波变换(discrete wavelet transform,DWT)和自协方差(autocovariance,AC)模型的组合(DWT+AC)[13]等.为了对比不同模型的预报效果,维也纳理工大学大地测量与地球物理研究所从2005年10月1日至2008年2月28日组织了全球性的地球定向参数预报比较竞赛,2 yr多的预报结果表明,没有一种模型既适合于ERP所有分量的预报又适合于所有跨度的预报[14].

受多种激发因素的影响,地球自转变化呈现复杂的非线性不规则变化特性,因此采用非线性的预报方法对其进行预报在理论上更为合理[8−11].ANN是逼近复杂非线性函数的一种有效工具,所以有许多学者将其应用于ERP预报中,并取得了显著的预报效果[3,5,8−11].但是ANN存在一些缺点,例如网络拓扑结构难以确定、训练过程存在过学习现象、迭代过程易陷入局部最优、收敛速度较慢.此外,ANN的优化目标是基于经验的风险最小化,无法保证网络的泛化能力[15].针对上述缺点很多学者提出了改进措施,例如将ANN和其他人工智能(arti fi cial intelligence,AI)算法结合进行网络优化,如遗传算法与ANN的组合、粒子群算法与ANN的组合等,然而这些算法均存在一定的不足,仍处于不断的尝试研究阶段.

GP是近年来发展起来的一种新型机器学习算法,它有着严格的统计学习理论基础,对处理高维数、小样本、非线性等复杂问题具有很好的适应性,且泛化能力强[16−17].与ANN和支持向量机相比,GP具有容易实现、超参数自适应获取、非参数推断灵活以及输出具有概率意义等优点[16−17].现已成为机器学习领域的研究热点,并在许多领域得到了成功的应用[18−19].本文将GP算法用于LOD变化预报,研究表明:GP算法用于LOD变化预报是可行的,且预报效率和精度较高.

2 GP的基本原理及其用于LOD变化预报的方法

2.1 GP回归的基本原理

高斯过程又称正态随机过程,其任意有限变量集合都有着联合高斯分布的特性,即对于任意的变量x1,x2,···,xn与其对应的函数f(x1),f(x2),···,f(xn)的联合概率分布服从n维高斯分布.高斯分布的全部统计特征完全由它的均值函数m(x)和协方差函数C(x,x′)来确定,一般记为f(x)～GP(m(x),C(x,x′)).

若给定训练样本集D={(xi,yi)|xi∈Rd,yi∈R,i=1,2,···,n},其中d为向量xi的维数,则对于测试样本输入x∗,GP模型的预测值为

其中k(x∗)=[C(x∗,x1),C(x∗,x2),···,C(x∗,xn)]为测试样本输入和训练样本输入值之间的1×n阶协方差矩阵,K是训练样本输入值之间的n×n阶协方差矩阵, Kij=C(xi,xj).(1)式表明,GP模型可以根据协方差函数和测试样本进行预测.

协方差函数在GP回归模型中起到关键作用,它表达了一种样本间的相似性,对所要学习的函数提供了假设信息,协方差函数必须是半正函数.常用的协方差函数为平方指数函数[17],即

2.2 基于GP回归模型的LOD变化预报

2.2.1 数据预处理

本文所用的LOD变化数据来自IERS发布的EOP C04序列,采样间隔为1 d.LOD变化序列中周期为5 d～18.6 yr的固体地球带谐潮汐项可以通过IERS协议给出的经验模型精确确定[2],近周日和半周日海洋潮汐项不作修正,LOD变化的长期趋势项、季节性变化的周年和半年项等根据下述线性模型确定[4]:

其中LODR表示经过固体地球带谐潮修正后的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为协调世界时(Coordinated Universal Time,UTC),在拟合时单位转换为d.a、b表示长期趋势项的参数,c1,1、c1,2表示半年项的参数, c2,1、c2,2表示周年项的参数,c3,1、c3,2表示9.3 yr项的参数,c4,1、c4,2表示18.6 yr项的参数,这10个未知参数通过最小二乘法求得.

经过上述数据预处理后的剩余部分为含有非线性成分的残差序列,主要包括海洋近周日、半日潮项以及不规则的短周期成分.图1从上至下依次绘出了1990—2010年期间LOD变化的原始序列、带谐潮序列、线性模型拟合序列以及残差序列.本文采用GP回归模型对残差序列进行建模和预测,将线性模型的预报结果和残差序列的预报结果相加即可获得最终的LOD变化预报值.

2.2.2 建模和预报

GP模型的建模过程就是通过对样本数据D的训练,确定协方差函数的超参数.超参数的选取方法主要有交叉检验法、贝叶斯推理法和最大似然法[13].本文采用最大似然法选取超参数,即任意给定超参数的初值,采用共轭梯度优化算法求取训练样本对数似然函数的最大值,从而得到似然函数最大值所对应的超参数即为最优超参数.对数似然函数的形式为

图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

除了协方差函数及其超参数,样本的输入和输出方式也非常重要.按以下方式构建样本的输入和输出:

在训练阶段,样本的输入和输出方式为

其中{ξ(i),i=1,2,···,n},表示LOD变化序列经数据预处理后的残差序列,根据经验确定,本文取d=4.

在预报阶段,预报跨度为k=1,2,···,d,d+1,···时样本的输入和输出方式分别为

2.2.3 精度评定指标

采用均方根误差(root mean square error,RMSE)和平均绝对误差(mean absolute error,MAE)作为预报结果的精度评定指标,其计算公式分别为

其中i为预报跨度,N为预报期数,分别表示第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变化残差的预报结果(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变化预报结果同Schuh等[3]使用的反向传播神经网络(back propagation neural networks,BPNN)和张晓红等[8−9]使用的改进的BPNN以及广义回归神经网络(general regression neural networks,GRNN)的预报结果进行对比,预报RMSE值见表1.

从图2(a)、(b)和表1可以看出,GP用于LOD变化预报是可行和有效的,随着预报跨度的增大,预报精度有所降低.

表1 GP预报结果与BPNN预报结果(Schuh等[3]的预报结果)、改进的BPNN及GRNN预报结果(张晓红等[8−9]的预报结果)的比较(单位: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的预报精度仍优于BPNN的预报精度,但低于改进BPNN和GRNN的预报精度.在预报效率上,由于GP模型参数(协方差参数)可以自适应获取,而不必像ANN技术需要对训练样本反复训练才能得到最优网络参数,故训练速度较ANN技术要快.此外,因为本文采用的GP预报模式只需建模一次便可实现LOD变化的多天连续预报,因此预报所用时间较少,对于1～360 d的连续预报,训练时间和预报时间之和一般在20 min以内,保证了算法的实时性.

图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的预报结果进行比较,选取了1990年1月1日至2005年9月30日的LOD变化数据作为基础序列,预报2005年10月1日到2008年2月28日(与EOP PCC预报时间段相同)1～360 d跨度的LOD变化值,统计了预报结果的MAE,并与EOP PCC的结果进行了对比,对比结果见图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的预报精度优于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 讨论与总结

本文根据实验验证了GP用于LOD变化的可行性和有效性.与ANN方法相比,GP方法较容易实现,并且它不需要太多的先验信息,只需事先选择适当的协方差函数,其超参数在训练过程中便可以自适应地确定,从而可以避免预报的人为主观性,提高预报结果的可信度.通过实例发现,GP方法用于LOD变化预报可以取得较好的预报效果.通过与ANN预报结果以及与EOP PCC预报结果的比较发现,GP方法的短期(1～30 d)预报精度较高,但中期(1～360 d)预报精度则不如ANN预报精度和EOP PCC预报精度,这可能是由本文所使用的递推预报模式的误差累积效应引起的,对此可以尝试以下两种方法对GP中期预报精度进行改进,一种是改进样本输入方式,如采用连续输入方式或者跨度输入方式[20−21],另外一种方法是对GP模型进行在线训练,本文暂不对此进行讨论.在预报效率方面,由于本文采用递推预报模式进行预报,故只需一次模型计算便可以实现多步预报,极大提高了预报效率.对于跨度为1～360 d的预报,应用GP方法预报LOD变化只需要20 min左右的时间,而应用ANN则需要数小时,预报效率大大提高,这对于ERP的实时快速预报具有重要的现实意义.此外,基于GP方法预报LOD变化所需训练样本数量远远小于ANN所需样本量,这在历史数据较少的情况下进行LOD变化的预报具有更高的现实意义.

由于协方差函数类型、超参数选取方法以及样本输入方式等对GP方法的预报效果都有一定的影响,因此,如何从上述角度来优化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