基于混沌时间序列的Volterra自适应滤波极移预报方法*

雷 雨 赵丹宁 乔海花 徐劲松 蔡宏兵

(1 西安邮电大学计算机学院 西安 710121)(2 宝鸡文理学院电子电气工程学院 宝鸡 721016)(3 中国科学院国家授时中心 西安 710600)(4 江苏师范大学江苏圣理工学院 徐州 221116)

1 引言

地球自转变化可以用一组地球定向参数(Earth Orientation Parameters,EOP)来表征,包括极移(polar motion,PM)、日长变化、岁差与章动.其中极移是地球瞬时自转轴相对于地球本体内的运动,是实现地球参考坐标系和天球参考坐标系相互转换的必要参数之一,在卫星导航、深空探测及跟踪等空间科学工程应用领域具有不可或缺的作用.此外,极移与地球上物质分布等地球物理现象存在密切关系,因此,极移可以为全球性的地球物理现象提供重要信息.虽然甚长基线干涉测量(Very Long Baseline Interferometry,VLBI)、全球卫星导航系统(Global Navigation Satellite System,GNSS)以及卫星激光测距(Satellite Laser Ranging,SLR)等空间测地技术能获取高精度的极移数据,但因复杂的资料处理过程,在极移观测时效性上存在一定的滞后性,难以获得实时的极移观测量,所以只能通过预报手段实现极移参数的实时获取.卫星导航和深空探测器跟踪、导航等空间科学工程应用领域对高精度极移预报数据有重大需求.其中深空探测及跟踪通常需要1–7 d的极移预报数据以支持深空探测器跟踪和导航,而北斗卫星导航一般需要1–60 d的极移预报数据作为卫星自主定轨的先验信息.因此,如何提高极移短期和中期预报精度成为天文地球动力学领域的一个研究热点[1–2].

国内外学者提出了多种极移预报方法,例如最小二乘(Least Squares,LS)外推[3–5]、自回归(Autoregressive,AR)[6–8]、神 经 网 络[9–11]、小 波 分析[12]和奇异谱分析模型[13]等.国际地球定向参数预报比较竞赛(Earth Orientation Parameters Prediction Comparison Campaign,EOP PCC)结果表明,最小二乘和自回归模型的组合是预报极移的较优模型,但该模型无法对任意跨度的极移预报均取得最优[14].当前,国际地球自转服务局(International Earth Rotation and Reference Systems Service,IERS)A公报(Bulletin A)中发布的极移预报产品就是利用最小二乘和自回归模型组合获取的[15].

受大气和海洋等多种激发因素的影响,极移序列不仅包含周年极移与钱德勒极移等确定周期性成分,还包含准周期性甚至非线性不规则成分,其中,周期性成分预报精度较高,而准周期性和非线性成分具有非线性、非平稳特征,是影响极移预报精度的主要因素[16].混沌现象是确定性系统中因系统随机性而产生的外在的、复杂及貌似无规律的行为,混沌系统对初始值敏感的特性使混沌系统输入的变化可以快速地反映在系统输出中,因此混沌理论提供一种更符合现实世界的非线性建模方法.混沌预报法作为一种新颖的非线性预报理论,通过对时间序列进行相空间重构,能够挖掘时间序列中隐含的规律,非常适合于总体呈现确定性但又具有随机性的复杂系统,所以混沌时间序列预报方法用于极移预报是合理可行的.混沌时间序列预报方法可以分为局部预报法[17]、全局预报法[18]和自适应预报法[19].其中,局部预报法只利用部分样本数据拟合非线性函数,计算速度快但预报精度低;全局预报法利用全部样本数据拟合函数,当样本数据量大时计算速度慢且易过拟合;自适应预报法能够自适应地跟踪混沌运动轨迹,可以根据混沌时间序列的变化自动调整模型参数,具有较高的预报精度.而针对混沌时间序列的Volterra自适应预报法能够有效地预测低维混沌时间序列,已在电力、交通和天气预报等领域取得成功应用[20–22].

本文探索用新型的混沌动力学理论来研究极移变化.首先,把只有一维的极移时间序列通过相空间重构方法重构出高维空间,恢复极移内在的非线性特性;其次,对呈现非线性特性变化的极移序列应用混沌理论分析与预报,通过对极移序列的相空间重构计算Lyapunov指数,确定极移为混沌时间序列的前提下,结合Taken嵌入定理,确定Volterra自适应预测滤波器的长度,利用Volterra自适应滤波对极移时序进行预测.文中对极移进行了未来60 d的预报,并将预报结果与EOP PCC和IERS A公报的预报结果进行了对比,对比结果验证了本文方法的预报效果.

2 时间序列相空间重构及混沌识别

常规的时间序列通常是在时间域中进行时间序列的分析与处理,对于混沌时间序列的研究在相空间中进行.所谓混沌时间序列,可以视作是考察混沌系统所得到的一组随着时间而变化的观察值.根据Taken嵌入定理,可以从一维混沌时间序列中重构一个和原动力系统在拓扑意义下相同的相空间,混沌时间序列的识别、分析以及预测均在重构的相空间中进行,因此,混沌时间序列的相空间重构是混沌时间序列分析与处理的关键.

2.1 相空间重构

Packard等人于1980年提出了时间序列相空间重构的两种方法[23]:导数重构法与坐标延迟重构法,其中,坐标延迟重构是通过一维时间序列的不同延迟时间τ来构建D维的相空间矢量.设有长度为N的一维时间序列x(t),t=1,2,···,N,嵌入维数为m,则时间序列x(t)重构的相空间为

式中,M=N-(m-1)τ,X为M×m维矩阵,任一相点X(t)均包含m个元素,相空间的演化轨迹就是M个相点在m维空间的连线.相空间重构的关键是选取适当的τ与m.τ的选取方法有互信息法和自相关法,m的选取方法有最近邻点法与饱和关联维数法.近年研究表明,τ与m是相关的,联合选取可以提高相空间重构的质量[24].本文利用C-C方法联合求取τ与m,该方法使用关联积分同时估计出时延与嵌入窗,具体过程见文献[24].

2.2 混沌识别

时间序列的混沌识别主要是通过计算混沌特性参数来识别序列的混沌特性,混沌参数有K熵、关联维数与Lyapunov指数等.本文选择Lyapunov指数作为混沌参数来识别极移时间序列的混沌特性,如果最大Lyapunov指数(Largest Lyapunov Exponent,LLE)为正数表明时间序列具有混沌特性,且最大Lyapunov指数越大,混沌特性越强.

基于小数量法计算最大Lyapunov指数的步骤如下[25]:

(1)对时间序列进行快速傅里叶变换估计平均周期T;

(2)基于坐标延迟重构法重构时间序列的相空间,对相空间中的每个点X(t),计算该相点和最邻近点X(ˆt)的l个离散时间步长后的距离dt(l),t为相点X(t)的历元,ˆt为相点X(t)的最邻近点的历元,dt(l)可以表示为

式中,l=1,2,···,min(M-t,M-ˆt),ˆt=1,2,···,M,且t/=ˆt;

(3)取dt(l)的对数lndt(l),对每个l求所有非零lndt(l)的均值再除以时间序列的采样周期Δt,即

式中,q为非零dt(l)的个数,i为非零dt(l)的序号;

3 Volterra自适应预测滤波器

对于极移预报,其本质上是建立样本“输入–输出”之间的映射关系,极移预报的准确与否在一定程度上也是由样本“输入–输出”之间的拟合程度决定的.在极移混沌时间序列预报中,设样本输入为重构相空间中的相点X(n)=[x(n),x(n-τ),···,x(n-(m-1)τ)],输出为x(n+1),n为历元,则理论上存在一个映射模型F(·)使得:

由于Volterra级数能够同时对时间序列中的线性行为和非线性行为进行拟合,本文利用Volterra级数构建样本“输入–输出”映射模型来拟合F(·),Volterra级数的展开表达式为[26]

式中,i1,i2,···,ip为滞后历元,hp(i1,i2,···,ip)为p阶Volterra核函数.这种无穷级数展开表达式在实际应用中较难实现,需要运用有限次截断与有限次求和的形式对Volterra级数展开式进行简化,较为常用的是二阶截断和求和的形式,即

式中,N1和N2为滤波器长度.对于混沌时间序列,根据Taken嵌入定理:一个混沌时间序列至少需要满足m≥2D+1才能完全描述原系统的动力学特征.因此,应将N1和N2取为N1=N2=m≥2D+1,则用于极移预报的二阶Volterra滤波预测器为

定义FIR(Finite Impulse Response)滤波器的输入矢量U(n)=[1,x(n),x(n-τ),···,x(n-(m-1)τ),x2(n),x(n)x(n-τ),···,x2(n-(m-1)τ)]T,核函数向量H(n)=[h0,h1(0),h1(1),h1(m-1),h2(0,0),h2(0,1),···,h2(m-1,m-1)T],二 阶Volterra滤波预测器可表示为

上式中的核函数向量可直接采用FIR滤波器的自适应算法来确定,本文选用最小均方自适应算法来确定核函数向量:

式中,e(n)为滤波器输出值和期望值之差,核函数向量根据e(n)自适应调整以获得期望输出;c为收敛控制参数,可控制滤波器收敛于某个期望值.

4 预报方法及精度评定

4.1 预报方法

基于混沌时间序列的极移预报方法的基本流程如下:

(1)极移数据预处理:利用最小二乘拟合扣除极移时间序列中的线性趋势项、周年项以及钱德勒项,获得趋势项、周年项以及钱德勒项的外推值,最小二乘外推模型的数学表达式为

式中,ax和bx表示极移Xp分量的线性趋势项参数,和表示Xp的钱德勒项参数,和表示Xp的周年项参数,p1和p2分别表示钱德勒摆动和周年摆动周期,在拟合中取p1=1.183 yr、p2=1 yr.对于极移Yp分量,各对应参数表示含义和Xp分量的模型相同.极移序列最小二乘拟合后的剩余残差序列包含亚季节性准周期性成分和非周期性成分,这些准周期变化和非线性变化主要由大气和海洋对极移的激发引起;

(2)极移残差序列相空间重构:基于坐标延迟重构法重构最小二乘拟合残差时间序列的相空间,其中,相空间的τ与m通过C-C法确定;

(3)极移残差序列混沌识别:采用小数据量方法计算最小二乘残差时间序列的最大Lyapunov指数,识别极移残差的混沌特征;

(4)极移残差序列预报:利用二阶Volterra自适应滤波器对重构相空间的极移残差时间序列进行预测.在建模阶段,滤波器输入为残差序列重构相空间中的相点X(n)=[x(n),x(n-τ),···,x(n-(m-1)τ)],输出为对应输入向量一步后的下一状态x(n+1),利用最小均方自适应算法建立输入X(n)和输出x(n+1)之间的映射关系;在预报阶段,将滤波器预测值ˆx(n+1)加入到输入向量X(n)中,并删除向量中第1个元素,然后将其更新为X(n+1),以此逐步递推,即可实现极移的多步预测.图1为基于混沌时间序列的极移预报流程.

图1 基于混沌时间序列的极移预报流程Fig.1 Flowchart of PM prediction based on chaotic time series

4.2 精度评定

选取平均绝对误差(mean absolute error,MAE)作为评价本文预报方法精度的指标[14]:

式中,k表示预报跨度,L表示预报期数,和分别表示第j期的第k天的极移预报值与观测值.

5 极移预报分析

5.1 极移数据来源

采用国际地球自转服务局发布的EOP C04序列和A公报资料进行极移分析与预报研究,EOP C04序列包含1962年至今的极移观测数据,该观测数据是通过VLBI、GNSS和SLR等多种空间测地观测资料联合解算得到的,极移观测精度约为50 μas,每周更新两次;A公报包含未来1 yr的极移预报数据,该预报数据是通过最小二乘外推和自回归模型组合外推获得的,由IERS快速服务与产品中心每周更新一次.IERS EOP C04序列和A公报公布的极移数据的采样间隔均为1 d,具有较高的精度和可靠性,是开展极移分析与预报研究的常用资料来源.

5.2 极移混沌特性分析

为验证极移的混沌特性,选取EOP 14C04序列中1997年1月1日至2020年12月31日的极移数据进行分析,将该时间段内的极移数据每隔6 yr划分为一个观测时段,共分为4个观测时段,以极移Xp分量为例.图2给出了每个观测时段内的极移观测序列、最小二乘拟合序列及其最小二乘拟合残差序列,其中,实线表示极移观测序列,点线表示极移最小二乘拟合序列,虚线表示极移最小二乘拟合残差序列.从图中可以看到,扣除线性趋势项、周年项以及钱德勒项后的极移最小二乘拟合残差呈现为明显的非线性变化,主要包含亚季节性准周期成分和剩余的不规则成分.

图2 极移Xp分量观测序列、最小二乘拟合序列及其残差序列Fig.2 Observed,fitted and residual time series of the PM Xp component

根据小数据量方法分别计算上述4个时段的极移最小二乘拟合残差序列的最大Lyapunov指数,计算方法为对(3)式的计算结果¯d(t)进行最小二乘拟合,拟合直线的斜率即为最大Lyapunov指数.限于篇幅,本文仅给出1997年1月1日至2002年12月31日的极移Xp和Yp序列的¯d(t)及其拟合结果,如图3所示,图中的阴影部分为所选的线性拟合区间,拟合斜率分别为0.0275和0.0159,其余时段的¯d(t)拟合斜率(最大Lyapunov指数)见表1,表中LLE表示最大Lyapunov指数.由表1可知,扣除线性趋势项、周年项以及钱德勒项后的极移残差序列具有低强度的混沌特性,说明极移Xp、Yp分量时间序列可视为混沌序列.本文将扣除线性趋势项、周年项以及钱德勒项后的极移最小二乘拟合残差作为一组混沌序列,利用二阶Volterra自适应滤波器对混沌序列中的有序行为进行外推预报,然后将最小二乘外推结果和Volterra预报结果合并即为最终极移预测序列.为验证本文方法的极移预报效果,将本文 预报结果与EOP PCC和IERS A公报的预报结果进行比较分析.

图3 ¯d随时间步长的变化及其拟合区间Fig.3variations with the time step and the fitted interval

表1 不同观测时段的极移序列的最大Lyapunov指数Table 1 Largest Lyapunov exponent of the PM in different observed periods

5.3 与EOP PCC比较

为评估各种方法的EOP预报效果,维也纳理工大学于2005年10月1日至2008年2月28日组织了国际性的EOP预报比较竞赛,国际上共有12个小组提交EOP预报数据,涉及20种预报方法.维也纳理工大学对20种EOP预报方法进行分析评估,通过分析评估发现极移预报精度较高的方法是最小二乘外推和自回归组合模型,其中,超短期(1–10 d)预报精度最高的小组是Kosek 051小组和Kalarus 061小组,使用的是最小二乘外推与自回归模型组合方法;短期(1–30 d)预报精度最高的小组是Kosek 051小组、Kalarus 061小组和Zotov 091小组,Zotov 091小组使用的是自回归模型;中长期(30 d以上)预报效果最好的小组是Kosek 051小组和Kalarus 061小组,其他小组编号及其使用的预报方法详见文献[10],此处不再赘述.

首先随机选择两个时段进行极移预报,选择的两个时段分别是2006年11月25日至2007年1月23日和2007年5月19日至2007年7月17日,极移两个分量的预报结果和预报误差如图4–5所示,图中虚线表示本文极移预报结果,实线表示EOP 05C04序列极移观测值,点线表示预报值和观测值之差.从图4–5中可以看出,本文方法1–30 d的极移预报值和观测值符合得很好,绝对误差小于5 mas,当预报跨度超过30 d时,预报值逐渐偏离观测值,预报误差表现出发散趋势,绝对误差可以达到12 mas.

图4 本文方法Xp预报值和EOP 05C04序列极移观测值的比较Fig.4 Comparison between the Xp predictions by the proposed method and the EOP 05C04 observations

为了与EOP PCC预报结果进行客观的比较,本文预报了与EOP预报比较竞赛同时段的极移数据,在进行预报结果比较时,采用EOP 05C04序列中的极移数据作为基础数据,基础数据长度取为周年极移和钱德勒极移的最小公倍数6 yr,选取平均绝对误差作为精度评估指标.混沌时间序列预测方法的1–10 d、1–30 d和1–60 d预报结果和EOP PCC预报结果的平均绝对误差如图6所示,EOP PCC预报结果由维也纳理工大学提供.由于各小组提交的EOP预报数据长度不一致,如Zotov 091小组仅提交超短期和短期预报结果,而未提交中长期预报结果,所以图6中1–10 d、1–30 d和1–60 d预报结果比较中所涉及到平均绝对误差曲线数量并不相同.图6中,Kosek 051小组预报结果的平均绝对误差用红色实线表示,Kalarus 061小组的平均绝对误差用黑色实线表示,Zotov 091小组用蓝色实线表示,本文方法的平均绝对误差在图6中用红色虚线表示,图中所涉及的1–10 d平均绝对误差根据约100期预报结果获得,而1–30 d和1–60 d平均绝对误差根据约50期预报结果获得.表2给出了Kosek 051小组、Kalarus 061小组、Zotov 091小组和本文预报结果的平均绝对误差数值.

从图6中和表2可以发现,本文方法的极移1–10 d超短期预报精度和Kosek 051小组、Kalarus 061小组相当,显著优于其他小组;对于10–30 d短期预报,本文方法的平均绝对预报误差与预报精度最高的Kosek 051小组、Kalarus 061小组和Zotov 091小组基本相当;在30–60 d中期预报中,本文方法的预报效果不如排在第1位Kalarus 061小组,极移Xp分量的平均绝对预报误差与排在第2位Kosek 051小组的预报误差处于同一水平,低于其他小组的平均预报误差;对于极移Yp分量,30–50 d的平均绝对预报误差与Kosek 051小组大体相当,当预报跨度超过50 d时,极移Yp分量的预报误差发散趋势明显,预报精度不及Kosek小组.从表2还可以看出,没有任何一种预报方法对于任意跨度和任何分量的预报精度都能够达到最优,说明各种方法均存在各自的局限性.

表2 本文预报结果和EOP PPC极移预报结果的平均绝对误差Table 2 Mean absolute error of the polar motion predictions obtained by the proposed method and EOP PPC

图5 本文方法Yp预报值和EOP 05C04序列极移观测值的比较Fig.5 Comparison between the Yp predictions by the proposed method and the EOP 05C04 observations

图6 本文方法和EOP PCC极移预报误差的比较Fig.6 Comparison of the PM prediction errors between the proposed method and EOP PCC

5.4 与IERS A公报比较

自2017年1月起IERS利用EOP 14C05序列作为极移基础数据发布常规极移预报值,本文选取EOP 14C05序列中的极移观测资料作为建模数据进行极移未来60 d预报,建模数据长度为6 yr,将预报结果与IERS A公报2021年第20周和25周发布的极移预报数据进行对比,并将EOP 14C04序列极移观测数据作为参考值.图7给出了本文方法与IERS A公报极移预报值以及EOP 14C04序列极移观测值的对比结果,图中点线表示IERS A公报极移预报值,虚线表示本文方法的极移预报值,实线表示EOP 14C04序列极移观测值.从图7可以看到,本文方法在预报前期与EOP 14C04极移序列符合得很好,和IERS A公报预报效果相当甚至更好,而在预报后期本文预测值逐渐偏离实际值,其中,本文方法对于极移Xp分量1–40 d预报精度较高,预报跨度超过40 d时误差较大;本文方法对于极移Yp分量预报的改进没有Xp分量那样明显,1–20 d预报结果与IERS A公报比较接近,预报跨度超过20 d时本文方法的预报误差开始增大.整体而言,本文方法对于极移短期预报效果较好,当预报跨度增加时预报结果不如IERS A公报.

图7 本文方法和IERS A公报极移预报值的比较Fig.7 Comparison of the predicted PM values by the proposed method and IERS Bulletin A

为进一步对比分析本文方法与IERS A公报的极移预测结果,将本文预报结果与IERS A公报2021年第1周至第35周发布的极移预报结果进行比较,将IERS EOP 14C04序列作为参考值,分别统计极移Xp分量和Yp分量的平均绝对预报误差,统计结果见表3.从表3可以看到,本文方法和IERS A公报对于不同跨度的极移预报结果存在一定程度的差异,相对于IERS A公报发布的预报值,本文方法的预报结果在预报前期效果有所提高,其中,对于极移Xp分量,未来1 d预报结果提高最为明显,预报精度比A公报提高15%.除第2 d、第3 d和第4 d外,本文方法对于未来30 d以内其他跨度的预报有1%–12%程度的提高.对于极移Yp分量也是如此,只是本文方法对于Yp分量的预报改善程度没有Xp分量明显,改善程度在6%以内.无论是对于Xp分量还是Yp分量,在预报后期,即当预报跨度超过一定的天数后,本文方法的极移预报结果相对于IERS A公报非但没有提高,反而有所下降,其中,当预报跨度超过30 d时,本文方法对极移Xp分量的预报结果开始恶化;当预报跨度超过20 d时,本文方法对极移Yp分量的预报效果逐渐降低,且Yp分量的预报效果降低得更为明显,这说明随着预报跨度的增加,本文方法的极移预报误差比IERS A公报增大得更快.根据第4.1节所述的极移预报过程,结合理论分析发现,出现这种现象的可能原因在于极移预报值是由Volterra自适应滤波器逐步递推得到的,这说明上一步滤波器输出预报值的预报误差会累积到下一步的预报值中,预报误差会出现逐步传递的情况,换言之,预报结果存在误差累积效应,且随着预报跨度的增加,这种误差累积效应会愈发明显.从表3还可以发现,本文方法对于极移Xp分量预报的改善程度比Yp分量更显著,这是由于Xp分量的混沌强度比Yp分量要强(见表1),因此基于混沌时间序列的预测方法对Xp分量具有更好的适用性.

表3 本文预报结果和IERS A公报极移预报结果的平均绝对误差Table 3 Mean absolute error of the polar motion predictions obtained by the proposed method and IERS Bulletin A

6 总结与展望

极移由趋势成分、周期成分与随机成分组成,是多种激发源的物理机制共同作用的结果,考虑到极移这种复杂的非线性时变特征,提出基于混沌时间序列的极移预报方法.首先利用最小二乘拟合算法分离极移的确定项和残差项,随后根据最大Lyapunov指数对极移残差项的混沌特性进行识别,计算结果表明:极移残差项的最大Lyapunov指数为正数,且数值较小,表明极移残差项具有低强度的混沌特性.在此基础上,构建二阶Volterra自适应滤波器对极移残差项进行预测,并叠加极移确定项的外推值得到极移预报值.将本文方法的极移预报结果和EOP PCC预报结果以及IERS A公报发布的极移预报产品分别进行了对比,分析发现:

(1)和EOP PPC极移预报结果相比,本文方法对于1–30 d的极移预报结果与EOP PPC最优预报方法相当,当预报跨度超过30 d时,本文方法的预报精度低于EOP PCC最优预报方法,但优于参与EOP PCC的其他方法;

(2)和IERS A公报发布的极移预报产品相比,本文方法的前期预报结果相对于A公报有不同程度提高,当达到一定的预报跨度以后,本文方法的预报精度出现明显的降低,总体而言,本文方法更适合于极移的短期预报.

通过分析发现,本文方法在后期预报中误差增大的主要原因可能是由Volterra自适应滤波器的误差累积效应导致的,因此,如何设计更为优良的Volterra自适应滤波器,从而提高极移的预报精度是下一步研究的重点.