钱 炜 岳建平 单丽杰 韩宸宇
1 河海大学地球科学与工程学院,南京市佛城西路8号,211100
日长(length-of-day,LOD)是地球自转参数之一,获取LOD时存在着时延问题,因此,对日长变化进行高精度预报具有重要的科学意义和应用价值[1]。传统的日长变化预报采用最小二乘拟合外推与其他残差修正的组合模型[2-4]方法。然而,最小二乘法是线性拟合方法,不能很好地拟合日长变化序列中实际存在的非线性特征,从而影响最终的预报精度。研究表明,大气活动与日长变化密切相关,甚至部分日长变化年际尺度信号可以完全由大气激发来解释[5-6]。文献[7-8]在预报日长变化时,在误差修正模型中引入轴向大气角动量(atmospheric angular momentum,AAM)序列,结果表明,顾及大气角动量能够有效提高日长变化的预报精度。
本文提出应用于日长变化预测的Prophet与向量自回归(VAR)组合模型。Prophet模型适用于拟合信号的非线性特征,VAR模型本质上是自回归模型(AR)在多元变量时间序列领域的推广。首先,利用Prophet模型代替最小二乘法,分别对日长变化参数序列及大气角动量序列进行趋势项及周期项拟合外推;然后,针对所得的两种时序残差,利用VAR模型再次进行拟合;最后,将两种模型预报结果相加,即可得到日长变化的最终预报结果。
Prophet模型是由Facebook在2017年构建的一种用于时序分析的新模型,该模型算法将时间序列y(t)分解为长期趋势项、季节周期项以及节假日效应项等特征,并通过贝叶斯曲线拟合分别对各特征序列进行建模及预测,最后根据如下广义加法模型对特征序列进行累加,以此得到完整的拟合及预测结果[9]:
y(t)=g(t)+s(t)+h(t)+ε(t)
(1)
式中,y(t)为完整时间序列,g(t)、s(t)、h(t)分别为趋势项、周期项及节日项,ε(t)为误差项。趋势项g(t)的逻辑回归函数为[9]:
(2)
式中,C为趋势上限,k为增长率,m为偏移量。随着时间t的延续,曲线以增长率k非线性增加,最后趋于上限值C[10]。周期项s(t)的模型为[11]:
(3)
式中,T为周期,N为频谱个数。
p阶的VAR(p)模型可以表达为[12]:
yt=c+A1yt-1+…+Apyt-p+εt
(4)
式中,c为n×1阶待估常数向量;Ai为n×n阶待估系数矩阵;yi为n维内生变量组成的向量;εt为随机扰动项,与自身滞后项及等式右边各变量严格无关[12]。构建VAR模型最关键的步骤是通过不同的筛选准则确定最优滞后阶数p。以最终预测误差(FPE)为例,其表达式为:
FPE(p)=
(5)
预报流程如下:1)在原始LOD时序中去除固体潮影响,得到LODR序列;2)利用Prophet算法将LODR序列进行特征分解,得到基于Prophet的拟合序列;3)利用LODR拟合的残差序列,以及经过Prophet拟合的AAM残差序列,共同建立VAR模型,从而得到顾及大气角动量的残差修正模型预测结果;4)将Prophet拟合序列和VAR拟合序列的预测结果相加,得到最终的LODR预测值。采用平均绝对误差(MAE)作为预测精度的评价指标,其计算公式为[15]:
(6)
日长变化数据为IERS发布的EOP14 C04地球自转参数序列(http:∥hpiers.obspm.fr/iers/eop),大气角动量序列来源于由欧洲中期气候预报中心(ECMWF)发布的轴向分量χ3序列(http:∥rz-vm115.gfz-potsdam.de:8080/repository)。序列区间为2008-01-01~2020-12-31,区间内序列保持连续,其中2008-01-01~2019-12-31为初始建模序列,预报跨度为150 d。每进行一次预报后将序列整体进行隔天前推,从而在维持预报序列总长度为12 a的情况下进行第二次预报。以此类推,总预报期数为217组,统计所有组1~150 d预报期间的平均绝对误差进行分析。图1为去除潮汐项后的LODR序列及原始AAM序列,可以看出,两个序列之间存在明显的同向运动规律,对其进行皮尔逊相关性分析,得到二者的相关系数为0.751,呈强相关,因此完全可以引入大气角动量序列来预报日长变化。
图1 去除潮汐项后的LODR序列及原始AAM序列Fig.1 LODR sequence after deducting tidal term and original AAM sequence
采用Prophet算法分别对LODR和AAM序列进行拟合,并预测2020年数值,结果如图2和图3所示。由图可见,经过Prophet算法拟合后,拟合曲线可以较好地反映相应原始序列的发展趋势。为更好地显示拟合效果,对LODR序列与其拟合序列作差,并统计残差序列的分布情况,将统计结果与经过最小二乘拟合去除趋势项及周期项的残差统计信息进行对比。表1(单位ms)给出了LODR序列分别经过Prophet、LS拟合后残差序列的均方根误差(RMSE)、平均值、最大值及最小值信息。统计发现,采用Prophet算法来拟合外推LODR序列更为有效。
图2 Prophet模型拟合预报LODR效果Fig.2 The results of fitting and predicting LODR by Prophet model
图3 Prophet模型拟合预报AAM效果Fig.3 The results of fitting and predicting AAM by Prophet model
表1 LODR残差序列的统计信息
图4为拟合后的LODR及AAM残差序列,二者之间的相关系数为0.815。这说明残差序列之间依然存在着较强的相关性,并且两个残差序列均为平稳序列。因此,可以考虑建立联合LODR残差序列及AAM残差序列的VAR误差补偿模型,再次对Prophet模型的LODR预测误差进行拟合修正。
图4 LODR及AAM残差序列Fig.4 Residual sequences of LODR and AAM
建立VAR模型的关键是确定模型的阶数q。利用初始残差序列,分别计算其1~100阶内的FPE、AIC及SC指数,结果如图5所示。图中虚线代表序列最小值处的阶数,分别为63、63、26。根据多数原则,本文选取p值为63作为VAR模型的阶数,从而构建VAR(63)模型用于拟合预测。值得注意的是,这仅仅是初始残差建模序列所构造的预报模型,当完成一期预报后,第二期的建模序列会整体向前推移1 d,此时需要重新构造VAR模型。
图5 3种定阶准则指数信息Fig.5 Index information of three order determination criteria
构建好VAR模型后,需要检验模型是否稳定,判断的准则是模型特征根倒数的模均小于1。检验结果如图6所示,图中点表示特征根的倒数。如果没有点落在单位圆外,就认为通过了平稳检验,VAR误差补偿模型有效。
图6 特征根检验Fig.6 Eigenvalues test
利用Prophet和VAR模型,可得到每期最终的LODR预测序列。将EOP14 C04序列作为真值序列,在217组预报全部完成以后,计算1~150 d跨度的MAE值,并将其作为精度评价指标。为了更好地分析预报效果,本文设计另外两种方案作为对比:
1)Prophet-AR模型。不考虑AAM序列,仅依据LODR序列进行日长预测。此时,VAR模型变为单变量的AR模型。将此方法与本文方法进行对比,可分析顾及大气角动量对预报效果的影响。
2)LS-AR模型。提取原始序列趋势项,选用传统的LS算法来代替Prophet算法。将此方法与本文方法及方案1进行对比,可分析Prophet算法与LS算法的优劣。
3种方案采用同样的预报流程,最终的MAE序列如图7所示。由图可见,在1~150 d的预报跨度内,本文给出的Prophet-VAR预报模型最优。为定量分析预报效果,表2(单位ms)列出3种方案在1~150 d预报跨度的MAE值,并分别计算Prophet-AR模型及Prophet-VAR模型相对于LS-AR模型的精度改善百分比。
图7 3种预报模型的MAE曲线Fig.7 The MAE curves of three prediction models
表2 3种预报模型的MAE值及对比LS-AR的改善程度
由表2可见,Prophet-AR与LS-AR相比,精度提升8%~40%,尤其是在超短期预报(1~10 d)中精度改善幅度更大。Prophet-VAR与LS-AR相比,精度提升16%~57%,在任意预报跨度内的精度提升都非常明显,且显著高于同时段Prophet-AR方法的改善精度,表明在日长变化参数预报中加入大气角动量序列,对预报精度有明显的提升作用。综上所述,顾及大气角动量的Prophet-VAR模型是一种非常有效的日长变化预报方法。
本文提出一种顾及大气角动量的Prophet-VAR日长变化参数预报模型,利用2008-01-01~2020-12-31期间的日长变化序列进行实验,以1~150 d的预报跨度来分析该模型的预报效果。为便于比较,同时设计了传统的LS-AR及不顾及大气角动量的Prophet-AR两种预报方案与本文方法进行对比分析,评价指标采用平均绝对误差。结果表明,Prophet算法代替LS算法来结合AR模型在预报日长变化参数时更有优势,这种优势在超短期预报时最为明显,且加入大气角动量序列对提高预报精度有显著的提升效果。由此说明,在预报日长变化参数时顾及大气角动量的Prophet-VAR模型是一种更有优势的高精度预报模型。
本文在分析过程中发现,尽管Prophet算法相较于LS算法能够更加有效地拟合LODR原始序列中的趋势项,从而提取出更多真实的非线性长期趋势项信号,但在提取周期项时,其与LS算法一样,仍然是提取振幅固定的周期信号,这并不符合日长变化序列的周期振幅时变特征。因此下一步的工作可以从该角度出发,探索如何通过提取更准确的周期项信号来提升日长变化参数预报的效果。