2017, Vol. 37
地球自转轴相对于地球表面的运动被称为极移。极移是实现高精度天球坐标系与地球坐标系转换的必要参数，也是卫星超快速定轨、卫星实时定轨以及深空大地测量等技术所必需的参数[1]。在18世纪人们就对极移进行研究，天文学家采用传统的光学仪器对极移进行检测，由于仪器容易受到外界环境的影响，导致测定的极移精度不太理想。随着观测技术的发展，多种现代空间大地测量技术被应用到极移的测定当中，极大地提高了极移的测定精度。虽然各种空间大地测量技术能够获取高精度的极移参数[1-2]，但都是事后处理数据，故极移参数在获取上存在一定的时延性，使得极移参数难以实时获取，只能采取参数预报的方法来满足相关技术对极移的需求[2-3]。
目前，国际上有大量的极移预报模型，这些预报模型主要分为两类，即线性预报模型和非线性预报模型。线性预报模型主要包括卡尔曼滤波(Kalman filter)[4-5]、谐波最小二乘外推(LS)[6]、自回归(AR)[7]等，而非线性预报模型主要包括人工神经网络(ANN)组合[8]、机器学习[9]等。为比较各预报模型的预报效果，Schuh等[1]发起了一次地球定向参数预报竞赛(EOP_PCC)，从2005-10-01开始到2008-02对地球定向参数进行了超短期、短期和中长期预报，结果显示，没有一个模型适用于所有间隔的极移预报, 就短期预报而言，LS+AR预报模型最佳。
采用LS+AR模型进行极移预报时，分别对PMX(极移X方向)数据以及PMY(极移Y方向)数据进行最小二乘拟合外推预报，然后再分别采用AR模型对LS拟合残差进行预报，残差的预报与LS外推值之和为最终的预报值。但极移是表示地球自转轴在地球表面的运动，被人为分解为PMX、PMY两个方向，且PMX和PMY受到的激发源的机理是一致的，PMX和PMY之间具有良好的相关性。因此，对PMX和PMY的残差分别进行预报是不恰当的，没有考虑两者之间的相关性信息。本文针对采用LS+AR模型进行极移预报的这一缺陷进行改进，提出了采用多元回归模型与最小二乘组合进行极移预报，即采用LS+MAR模型进行预报。
1 模型 1.1 LS模型极移(PMX, PMY)序列可分为线性和非线性两部分，其中线性模型包含周期项、趋势项等复杂过程，如式(1)所示[10]，可根据式(1)构建模型，进行最小二乘拟合及LS外推预报，而非线性部分指LS拟合残差值。
|
|
(1) |
式中，P(t)表示t(UTC时间)时刻极移变化量, 其中A、B、C1、C2、D1、D2表示待拟合参数，在最小二乘拟合中PSA、PA、PS通常分别取为0.5 a、1 a、1.183 a[8]。
1.2 AR模型自回归模型(auto regressive model, AR)是利用数据本身作为回归变量的过程，即利用前面若干期数据的线性组合来描述以后某时刻变量的线性回归模型[11-12]。假设一个时间序列y1, y2, y3, …, yN，若yt是包含前面p个序列的线性组合以及误差项的函数，则称该模型为p阶自回归模型(简称AR(p)), 其数学模型如下：
|
(2) |
式中，φ1、φ2、…、φn为模型参数，et是均值为0、方差为σ的白噪声。
模型参数φ最常用的估计方法为最小二乘方法[13]，这是因为该方法的参数估计比较简单，具有无偏性且精度高。因此本文选用最小二乘法求解AR模型参数，采用FPE准则来确定AR模型的最佳阶数[14]。
1.3 MAR模型多元线性回归(MAR)模型旨在在多变量的残差随机成分之间找到自动和互相关的参数，并将其应用于随机过程中建模和预测[15-16]。因此，MAR与AR不同，在建模中考虑了多随机变量之间的相关性。假设存在变量xti(i=1, 2, …, n; t=1, 2…, N)为第i个随机变量时间为t时刻的变量值，依据多元自回归模型可知，变量xti不仅与自身存在时间上的相关性，还与其他随机变量相关，受到其他随机变量的影响。MAR模型的数学表达式为：
|
(3) |
式中，anpn为第n个随机变量延迟第p步的自回归系数，etn是均值为零、标准差为σ的白噪声。对于MAR，本文采用最小二乘法求解其模型参数。MAR模型阶数的确定一般采用基于Schwarz贝叶斯准则(SBC)[16]，计算公式如下所示：
|
(4) |
式中，
由上述可知，LS+MAR模型与LS+AR模型在进行极移预报时既有相同部分又有不同部分，相同部分为对极移数据(PMX, PMY)进行最小二乘拟合外推，不同的部分在极移的残差预报方面。LS+AR模型分别对PMX、PMY的残差进行预报，而LS+MAR模型则考虑PMX与PMY拟合残差的相关性信息，结合PMX残差与PMY残差数据共同构建回归模型，实现PMX和PMY残差的同时预报。LS+MAR模型的预报流程如图 1所示。
|
图 1 LS+MAR模型的预报流程 Fig. 1 Forecast flow chart of LS+MAR model |
采用基础序列长度为6 a的极移数据进行LS拟合，获取PMX和PMY的拟合残差序列，计算两个残差序列之间的相关性系数。从2010-01-01起进行LS拟合，滚动拟合至2014-12-31，共进行了1 825期实验，相关系数统计见图 2。由图可见，绝大部分的相关系数在0.5以上，其中相关系数大于0.5的比例达到81.92%，相关系数大于0.8的比例达到55.45%。因此PMX拟合残差与PMY拟合残差具有较强的相关性。
|
图 2 相关系数 Fig. 2 Chart of correlation coefficient |
为检验本文提出的LS+MAR模型是否能够修正LS+AR模型的不足以及提高极移的预报精度，实验设计如下。
当基础序列(参与构建LS模型的数据量)为6 a时，LS+AR模型预报的极移精度最佳，这主要是因为相邻数据之间具有明显的相关性，而随着时间间隔的增加，相关性逐渐下降，但数据量过短则不能体现极移的周期特征[17]。因此，本文选取6 a的基础序列，根据LS+MAR模型进行极移短期预报，极移样本数据为IERS EOP 08C04 (http://hpiers.obspm.fr/iers/eop/eopc04/Ieopc04.62-now)，极移参数的样本间隔为1 d。从2010-01-01起进行预报，每次预报的时间跨度为30 d，每隔10 d预报一次(相邻预报的起始历元相差10 d)，滚动预报至2015-12-31。为验证LS+MAR模型的优越性，将其与LS+AR模型进行对比，预报精度分别采用误差绝对值的平均值(mean absolute error, MAE)作为评定预报精度的指标，其计算公式如下：
|
(5) |
式中，i为预测的天数，M为预测的期数，εi, j为预测值与观测值之差。
LS+MAR模型与LS+AR模型对比结果见表 1(单位mas)和图 3。可以看出，无论是PMX预报还是PMY预报，对超短期(1~10 d)或者短期(1~30 d)的预报，LS+MAR模型的预报结果明显优于LS+AR模型的预报结果，说明LS+MAR模型能够修正LS+AR模型的不足，有效提高预报精度。此外，LS+MAR模型相对于LS+AR模型的预报优势随着时间跨度的增加逐渐变大, 就预报精度提高幅度而言，PMY要大于PMX。
|
表 1 LS+MAR与LS+AR预报结果的精度(MAE)统计 Tab. 1 Accurate statistics of the forecast results (MAE) of LS+MAR and LS+AR |
|
图 3 LS+MAR与LS+AR预报结果对比 Fig. 3 Comparison of LS+MAR and LS+AR prediction results |
为与EOP_PCC的预报结果进行比较，采用与EOP_PCC相同的预报时间段进行预报，即采用LS+MAR模型从2005-10-01开始进行时间跨度为30 d的短期极移预报，滚动预报至2008-03-08。由于EOP_PCC的预报结果是以图的形式给出，故只能以图的形式进行比较，见图 4、图 5所示。
|
图 4 LS+MAR模型与EOP_PCC在超短期预报结果的对比 Fig. 4 Comparison of prediction results between LS+MAR model and EOP_PCC in ultra-short term polar motion prediction |
|
图 5 LS+MAR模型与EOP_PCC在短期预报结果的对比 Fig. 5 Comparison of prediction results between LS+MAR model and EOP_PCC in short term polar motion prediction |
在EOP_PCC的预报结果中，不同曲线分别表示不同参赛小组的预报结果。对极移预报而言，Kalarus、EOP产品中心和联合不同参赛小组预报方法的预报结果的精度最高。通过比较可知，在极移PMX方向上，对于超短期(1~10 d)预报，LS+MAR模型的预报精度要略优于EOP_PCC的超短期预报结果。对于短期(1~30 d)预报，LS+MAR模型的预报精度与在EOP_PCC中排在第4位(EOP产品中心)的预报精度相当。在极移PMY方向上，无论超短期还是短期的预报，LS+MAR模型的预报精度仅次于EOP_PCC中排在第1位(Kalarus)和第2位(联合预报)的预报结果。
由上述分析可知，LS+MAR模型的预报结果能够达到与EOP_PCC预报精度相当的水平，且LS+MAR模型的预报精度优于大部分EOP_PCC参赛小组的预报结果，其短期极移预报精度能够达到国际水平。
4 结语本文对LS+AR模型在预报极移中存在的问题进行了分析，即LS+AR模型在预报极移时是分别对PMX和PMY的残差进行预报，没有考虑两者之间的相关性信息。对此，本文采用MAR模型对PMX和PMY的残差同时进行预报，考虑了PMX残差和PMY残差之间的相关性。对比结果表明，LS+MAR模型的预报结果要优于LS+AR模型的预报结果，体现了LS+MAR模型的优越性。此外，通过与EOP_PCC的预报结果的对比也说明，LS+MAR模型在短期极移预报中能够获取与国际最好预报精度相当的预报结果。
致谢: 感谢国际地球自转服务(IERS)以及国际GNSS检测评估系统(iGMAS)对本文实验提供支持。
| [1] |
Dick W R, Richter B. IERS Annual Report 2004[R]. 2005
(  0) |
| [2] |
魏二虎, 刘文杰, Wei Jianan, 等. VLBI和GPS观测联合解算地球自转参数和日长变化[J]. 武汉大学学报:信息科学版, 2016, 41(1): 66-71 (Wei Erhu, Liu Wenjie, Wei Jianan, et al. Estimation of Earth Rotation Parameters and Delta LOD with Combining VLBI and GPS Observations[J]. Geomatics and Information Science of Wuhan University, 2016, 41(1): 66-71)
( 0) |
| [3] |
Kosek W, Popiski W, Kalarus M, et al. A Comparison of LOD and UT1-UTC Forecasts by Different Combination Prediction Techniques[J]. Artificial Satellites, 2005, 40(2): 119-125
( 0) |
| [4] |
Gross R S, Eubanks T M, Steppe J A, et al. A Kalman-Filter-Based Approach to Combining Independent Earth-Orientation Series[J]. Journal of Geodesy, 1998, 72(4): 215-235 DOI:10.1007/s001900050162
( 0) |
| [5] |
Akulenko L D, Kumakshev S A, Markov Y G. Motion of the Earth's Pole[J]. Doklady Physics MAIK Nauka/Interperiodica, 2002, 47(1): 78-84 DOI:10.1134/1.1450668
( 0) |
| [6] |
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
( 0) |
| [7] |
Schuh H, Ulrich M, Egger D, et al. Prediction of Earth Orientation Parameters by Artificial Neural Networks[J]. Journal of Geodesy, 2002, 76(5): 247-258 DOI:10.1007/s00190-001-0242-5
( 0) |
| [8] |
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 DOI:10.1007/s00190-010-0387-1
( 0) |
| [9] |
王小辉.改进极移预报的研究[D].长沙: 中南大学, 2013 (Wang Xiaohui. Improvement of Pole Motion Prediction[D]. Changsha: Central South University, 2013) http://cdmd.cnki.com.cn/Article/CDMD-10533-1014147423.htm
( 0) |
| [10] |
张昊.地球定向参数极移的预报理论与方法研究[D].长沙: 中南大学, 2012 (Zhang Hao. Researches on the Prediction Theories and Algorithms of Polar Moion of Earth Orientation Paramters[D]. Changsha: Central South University, 2012) http://cdmd.cnki.com.cn/Article/CDMD-10533-1012477699.htm
( 0) |
| [11] |
徐慧娟.自回归AR模型的整体最小二乘分析研究[D].抚州: 东华理工大学, 2012 (Xu Huijuan. The Autoregressive Model of Total Least Square Analysis[D]. Fuzhou: East China University of Technology, 2012) http://cdmd.cnki.com.cn/Article/CDMD-10405-1012030479.htm
( 0) |
| [12] |
蔺海新.自回归模型的平稳性研究[D].成都: 电子科技大学, 2008 (Lin Haixin. Research on the Stability of the Autoregressive Model[D]. Chengdu: University of Electronic Science and Technology of China, 2008) http://cdmd.cnki.com.cn/Article/CDMD-10614-2008122794.htm
( 0) |
| [13] |
陈国强, 赵俊伟, 黄俊杰, 等. 基于Matlab的AR模型参数估计[J]. 工具技术, 2005, 4(4): 39-40 (Chen Guoqiang, Zhao Junwei, Huang Junjie, et al. Matlab-Based Parameter Estimation of AR Model[J]. Tool Engineering, 2005, 4(4): 39-40 DOI:10.3969/j.issn.1000-7008.2005.04.011)
( 0) |
| [14] |
李平. 应用FPEC、AIC、FAIC准则确定AR模型阶的体会[J]. 辽宁石油化工大学学报, 1988(1): 39-43 (Li Ping. Application of FPEC, AIC and FAIC criteria to determine the order of AR model[J]. Journal of Liaoning University of Petroleum, 1988(1): 39-43)
( 0) |
| [15] |
Reinsel G C. Elements of Multi-Variate Time-Series Analysis[M]. Berlin: Springer, 1993
( 0) |
| [16] |
Schneider T, Neumaier A. Estimation of Parameters and Eigen Modes of Multivariate Autoregressive Models[J]. Acm Transactions on Mathematical Software, 2001, 27(1): 27-57 DOI:10.1145/382043.382304
( 0) |
| [17] |
严凤, 姚宜斌. 地球自转参数短期预报方法及其实现[J]. 大地测量与地球动力学, 2012, 32(4): 71-75 (Yan Feng, Yao Yibin. Short-Term Prediction Methods and Realization of Earth Rotation Parameters[J]. Journal of Geodesy and Geodynamics, 2012, 32(4): 71-75)
( 0) |