基于因果检验的非线性系统的预测试验*

张绿夷 王革丽 谭桂容 吴越

1) (中国科学院大气物理研究所,中层大气和全球环境探测重点实验室,北京 100029)

2) (中国气象局广州热带海洋气象研究所,广州 510641)

3) (南京信息工程大学,气象灾害教育部重点实验室,气象灾害预报预警与评估协同创新中心,气候与环境变化国际合作联合实验室,南京 210044)

4) (四川省气候中心,成都 610072)

非线性、非平稳系统的预测是一个具有重要科学意义的研究课题.最近一些工作已将收敛交叉映射算法(convergent cross mapping,CCM)用于检验变量之间的因果关系,由于在CCM 算法中,相空间中相互靠近的点在时间上具有相似的发展趋势和运动轨迹,因此该方法可以尝试应用于非线性、非平稳系统的预测试验研究中.鉴于此,本文将CCM 算法分别应用于Lorenz 系统和实际气候时间序列的预测中,并检测不同相空间重构方法对预测效果的影响.主要结果如下:1)不论是理想Lorenz 模型还是实际气候序列,对于单变量、多变量和多视角嵌入法3 种重构相空间方法而言,多视角嵌入法对变量的预测效果最好,表明对于给定长度的时间序列,重构相空间中包含的信息越多,其预测能力越强;2)将NAM (northern hemisphere annular mode)加入SAT (surface air temperature)的重构相空间中可以改善SAT 的预测效果.在使用单变量、多变量和多视角嵌入法进行预测时,利用复杂系统中变量中共有信息的特性,在时间序列长度一定的情况下,可以利用动力系统的复杂性来增加系统内的信息.基于因果检验的预测建模方法,通过挖掘数据中定量信息的提取,对非线性、非平稳系统预测技巧的改进提供了一个新颖的思路.

1 引言

对于气候预测而言,目前的预测建模理论大多建立在传统的统计预报、数值模式以及机器学习等方法上[1−5].气候系统作为一个非线性、非平稳的系统,其过程比混沌运动更为复杂[6],气候系统多因子间相互作用的复杂性是其预测理论和技术的瓶颈,有关非线性大气动力学的研究及预测依然是一个具有重大科学意义的前瞻性课题.

1980 年,Packard 等[7]提出了时间序列的相空间重构理论,通过一维时间序列的时间延滞来恢复原系统的动力学;1981 年,Takens[8]提出的嵌入定理,表明可以从一维时间序列中重构一个与原动力系统在拓扑意义下等价的相空间,奠定了非线性系统的理论基础.随后,一系列基于这些理论的非线性时间序列分析和预测方法便应运而生.在非线性时间序列预测中,Farmer 等[9]给出了混沌时间序列对单变量时间序列进行预测的方法;Casdagli[10]则比较了局域近似、全局近似和辐射近似3 种非线性预测方法的优劣.Yang 等[11]和Wang 等[12]则参考场时间序列的思想,分别利用神经网络建模方法和“场时间序列”的局域近似预测模型,对臭氧浓度和北半球500 hPa 高度场进行了预测试验,结果表明利用场时间序列信息可有效提高预测技巧.

然而,在一些非线性时间序列的分析和预测中,均存在一个满足遍历性定理的假设[13],即该时间序列所在系统的驱动力不随时间变化.气候系统作为一个复杂的非线性系统,控制真实气候系统的外部条件并非一成不变[14],此时,假定的遍历性定理也就不再成立[15].驱动力的改变不仅破坏了系统的平稳性,而且对时间序列分析和预测理论带来新的障碍[15].特别是,如何从实际的气候信号中识别并提取外部驱动力因子,分析导致气候变化的可能动力机制,并将其引入气候预测中是亟待解决的重要科学问题之一.

近年来,Verdes 等[16]和Wiskott[17]分别提出了从非平稳时间序列中提取外强迫因子的理论方法.Verdes 等[16]提出的“交叉预测法”是根据局部线性映射来反演由此引起的外强迫因子;Wiskott[17]提出的“慢特征分析法”则是通过提取快速变化的信号中的慢变特征,从而评估一个单一的外强迫因子[13,14].这两种方法都被应用于一些非平稳时间序列的分析中并得到了较好的结果[18,19].

此外,建立在因果关系上的驱动力分析,近年来也得到了长足的发展.Wiener[20]提出了一种因果关系的哲学概念,即因必须有助于改善果的预测.在此概念基础上,Granger[21]提出了著名的格兰杰因果关系(Granger causality),然而此种方法并不适用于复杂的非线性系统[22,23].2012 年,生物学家Sugihara 等[24]提出了基于相空间重构和Takens定理的收敛交叉映射算法(convergent cross mapping,CCM),该方法可以检验自然界中非线性动力系统中的因果关系,并已得到广泛的应用[25−27].例如,Zhang 等[28]利用 CCM 算法,探讨了北半球环状模(northern hemisphere annular mode,NAM)与东北亚地区冬季地面气温(surface air temperature,SAT)的信息传递,结果表明,二者存在单向因果关系,NAM 作为驱动力因子影响东北亚地区冬季SAT.

同时,由于在CCM 算法中相空间中相互靠近的点在时间上具有相似的发展趋势和运动轨迹,还可以尝试利用此方法对变量进行预测.因此,本文运用CCM 方法建立预测模型,并以Lorenz 系统以及东北亚地区冬季地面温度时间序列为例,将NAM 信号加入SAT 的重构相空间中,检验对SAT的预测效果.借助因果检验的手段识别影响气候变化要素的外强迫因子,并将其应用在实际的气候预测中,检验预测建模效果.本文首先简要介绍收敛交叉映射算法以及预测建模的思路,给出理想序列的预测检验和包含NAM 信息的东北亚地区冬季SAT 时间序列的预测分析,对非线性、非平稳系统预测技巧的改进提供了一个新颖的思路.

2 收敛交叉映射算法在预测上的应用

2.1 CCM 算法

CCM 算法基于相空间重构理论和Takens 嵌入定理,其基本思想是通过两变量间的重构相空间的相互预测来判定二者之间的因果关系[24].若变量Y为变量X的驱动力因子(Y影响X),利用变量X的重构相空间MX预测Y时,随着时间序列长度L的增大,邻近点间的距离不断减小,逐渐收敛于Y(t),并且与观测值Y(t) 的相关系数会逐渐增大,且二者的相关系数不断增大且达到收敛.

假设变量X,Y的时间序列长度为L,重构相空间的嵌入维度为E,采样间隔为τ,在t时刻X和Y的重构相空间MX,MY坐标为

2.2 收敛交叉映射算法中三种重构相空间的方法

利用CCM 算法可以使用历史轨迹上的信息来预测未来值,而不是依赖于求解固定的方程式或方程组.相空间重构的方法是利用CCM 算法预测的关键.到目前为止,在动态经验模型的框架中,已提出三种不同的方法分别为:单变量嵌入法[29,30]、多变量嵌入法[31,32]多视角嵌入法[33].

单变量嵌入法使用单个变量的时间滞后值来重构相空间,首先将预测变量的时间序列分为X,Y两部分,然后用X的重构相空间MX来预测Y.t时刻重构相空间MX的坐标为

其中,E为重构相空间的嵌入维数,τ为采样间隔.由(7)式得到重构流形MX={x(t)}.定义(t+1)|MX为利用重构相空间MX预测得到的Y(t+1).

确定t时刻MX上的向量x(t),并在MX上找到距离其最近的E+1 个邻近点,离其最近的点记为x(t1),第二近的点为x(t2),依次类推.根据相空间的运动轨迹,找到E+1 个邻近点在t+1 时刻的状态点x(t1+1),x(t2+1)···x(t(E+1)+1).则Y(t+1)为

其中wi是MX上x(t)和其第i个邻近点的距离的权重.

多变量嵌入法为在单变量嵌入的基础上,使用多个变量的时间滞后来重构相空间.例如,在M系统中包含C1和C2两个变量,利用多变量嵌入法重构M的相空间,t时刻重构相空间的坐标为{C1(t),C1(t–τ),C2(t)},其重构相空间与M具有微分同胚的特点.

多视角嵌入法则通过组合多种与预测变量有关的变量来增加重构相空间中的信息[33],根据嵌入理论[8],即使系统中只存在几个变量,通过多种有效嵌入也是可预测的.在包含n个变量的动态系统中,假设嵌入维数为E,允许最大的滞后时间步长为l,利用多视角嵌入法重构相空间时共有种组合形式.根据重构相空间中嵌入变量的不同,不同重构相空间中包含的信息有所差异.不同的重构相空间是从不同的角度对原始动力系统进行重构,其重构相空间均与原始动力系统的相空间微分同胚.利用CCM 算法进行预测时, 根据m个重构相空间预测能力大小进行排序, 通过对前k=构相个空预间测预值测求能平力均大得小到进最终的预测结果.尽管所有变量组合都是有效的嵌入,但是使用有限的数据可能无法很好地解决系统动力学问题.因此,多样视角嵌入中仅使用了按样本内预测能力排名的前k个进行重构,然后将从top-k重构中预测的值平均,计算单个预测值.

至此,本文已经简要说明了CCM 算法以及该算法中包含的三种嵌入方法.接下来将用一组理想数据实验以及实际的气候时间序列来讨论CCM算法在预测上的应用.

3 理想数据实验

为了讨论CCM 算法在预测上的应用以及上述3 种嵌入方法的异同,本文将Lorenz 系统[34]作为理想模型来进行实验设计.Lorenz 模型是描述大气对流等问题的简化模型,其方程的表达式如下:

当σ=10,β=8/3,ρ=28 时,系统出现混沌现象并且产生奇异吸引子,其中ρ为Rayleigh 数,是系统的控制变量.令Lorenz 系统的初值x(1)=y(1)=z(1)=10−17,步长 ∆t=0.1,t取[0,8000],迭代80000 次,取后1000 个数据作为实验序列,并对系统中的x分量进行预测.

在使用单变量嵌入法进行预测时,首先要确定嵌入维数,比较不同的嵌入维数E下x分量对应的时间序列重构相空间的预测的平均绝对误差MAE (图1).结果表明,当嵌入维数E=3 时,系统的平均绝对误差最小,将x分量嵌入到一个三维相空间中可以很好地恢复其原始系统的动力学特征.因此,当对Lorenz 系统中x分量进行预测时,选择最佳嵌入维数E=3 (见图1),预测步数tp=1.

图1 最佳嵌入维数E 的选取Fig.1.Selection of the best embedding dimension E..

图2(a)—(c)给出了分别使用单变量嵌入法、多变量嵌入法、多视角嵌入法对上述Lorenz 系统中x分量重构相空间的预测结果,预测步数tp=1.实验结果表明,在CCM 算法中三种重构相空间方法的预测能力从大到小依次为:多视角嵌入法、多变量嵌入法、单变量嵌入法,且预测能力越强,对应的平均绝对误差就越小(图2(d)).多视角嵌入法的预测能力优于多变量嵌入法和单变量嵌入法,说明利用复杂系统中变量中共有信息的特性,可以提高对系统的预测能力.

图2 (a)单变量嵌入法预测结果;(b)多变量嵌入法预测结果;(c)多视角嵌入法预测结果;(d)三种嵌入法预测结果的平均绝对误差比较Fig.2.Forecast results:(a) Univariate embedding;(b) multivariate embedding;(c) multiview embedding;(d) average absolute error of the prediction results of the three embedding methods.

图3 给出了利用三种嵌入方法对Lorenz 系统中x分量的多步预测结果,可以看出,随着预测步数tp的增加,三种方法的的预测效果均逐渐减弱.同时,考虑到真实气候系统数据的长度,作为初步探索,下文仅考虑了向后一步的预测.

图3 三种嵌入方法预测能力随步长的变化 (a)单变量嵌入法;(b)多变量嵌入法;(c)多视角嵌入法Fig.3.The prediction ability of three embedding methods varies with step size:(a)Univariate embedding;(b) multivariate embedding;(c) multiview embedding;.

利用Lorenz 系统对比了三种重构相空间方法的预测能力并简单讨论了预测能力随预测步数tp的变化.接下来为了更好地说明CCM 算法在预测上的应用,将其应用于实际大气时间序列中,并对其预测试验结果进行检验比较.

4 东北亚地区冬季SAT 时间序列的 预测分析

4.1 所用资料

1)美国气候预测中心(Climate Prediction Center,CPC)提供的逐月北半球环状模指数(NAM index,NAMI).

2)美国气象环境预报中心/美国国家大气研究中心(NCEP/NCAR)提供的 1000 hPa 逐月地面温度,水平分辨率为2.5° × 2.5°.

本文主要研究的区域为东北亚区域:40°—50°(N),90°—130°(E).时间跨度为1950 年11 月—2019 年2 月.冬季定义为11 月至次年2 月.

4.2 单变量、多变量嵌入法的预测试验研究

Zhang 等[28]的工作已经表明,NAM 与东北亚地区冬季地面气温SAT 存在因果关系,NAM 作为驱动力因子能够影响东北亚地区冬季SAT. 因此,本文在此工作基础上,使用CCM 并结合单变量、多变量嵌入法及多视角嵌入法对东北亚地区冬季SAT 进行预测试验研究.

首先使用单变量嵌入法进行预测(图4(a)),单变量嵌入法为使用单个变量的时滞特性来重建相空间.在预测试验中,选取东北亚地区冬季SAT 系统中的2 月份重构系统的相空间.允许最大的滞后时间步长max_lag=4,采样间隔τ=1,预测步数tp=1,t时刻重构相空间的坐标为{X2(t),X2(t–τ),X2(t–2τ),X2(t–3τ)}.理论上,其重构相空间虽然在形状、大小方面发生了变化,但是两个不动点以及其基本的动力特征没有发生改变[8].然而由于时间序列长度的限制,在较短的时间序列中,SAT 的重构x吸引子非常稀疏,阻碍了对来自临近点的动力状态的准确推断.使用单变量嵌入法对SAT 进行预测的预测能力仅为0.12,加入NAM 信号后预测能力稍有提升,但预测效果仍不佳,其预测能力为0.18.

图4 (a)单变量嵌入法和在目标变量中加入NAM 信号后的预测结果;(b)多变量嵌入法的预测结果(黑色圆圈代表仅利用SAT 序列进行预测,红色三角代表加入NAM 信号后)Fig.4.(a) Univariate embedding method and prediction result after adding NAM signal to the target variable;(b) prediction result of multivariate embedding method (black circle represents prediction using only SAT sequence,red triangle represents after adding NAM signal).

然后运用多变量嵌入法来重构东北亚地区冬季SAT 的相空间(图4(b)).与单变量嵌入法不同的是,多变量嵌入法使用多个变量来重构相空间,而不是运用原变量的滞后值.因此在对2 月份的SAT 的相空间进行重构时,可将其前三个月(即11,12,1 月)的值视作为系统中的另外三个变量,并且由于在冬季系统中,NAM 信号与东北亚冬季SAT 存在单向因果关系,NAM 是东北亚地区冬季SAT 的驱动力因子,二者微分同胚,所以可将变量NAM 中的信息加入SAT 的重构相空间中来预测SAT.分别选取1950—2018 年11 月-次年2 月的SAT 以及NAM 信号,利用多变量嵌入法重构2 月的SAT 相空间,该相空间中吸引子同样具有微分同胚的特点.图4(b)中黑色圆圈代表仅用多个月份的SAT 信号嵌入进行预测的结果,红色三角代表加入NAM 后运用多变量嵌入法对2 月SAT 进行预测的结果.

对比图4(a)和图4(b)可以看出,多变量嵌入法的预测能力大于单变量嵌入法,将NAM 加入冬季SAT 系统中后,对SAT 的预测能力会显着提高,并且预测值和观测值之间的相关系数更高(将红色虚线与黑色虚线进行比较).然而,尽管多变量嵌入法的预测能力大于单变量嵌入法,SAT 的实际观测值和预测值仍存在很大误差,表明系统中仍然存在很多无法解释的变化.此类方法可能会受到序列长度的限制,并且可能会受到噪音的影响.此外,观测误差将导致精度降低,即使时间序列足够长以密集地吸引吸引子,临近点也可能无法形成平滑的曲线,从而影响预测结果.

上述结果表明,单变量嵌入和多变量嵌入的方法十分依赖数据长度,时间序列长度越长,说明序列中可能包含的信息越多,用来重构相空间的点就越多;若数据长度过短,则不能很好地恢复系统中时间序列的动力学特性.因此在时间序列长度较短时,两种方法的预测能力均不佳.但将NAM 信号加入SAT 的重构相空间中,增加了SAT 相空间中的信息,有助于改善其预测效果.

4.3 多视角嵌入法的预测试验

与上述两种嵌入方式不同,多视角嵌入通过组合多种与预测变量有关的变量来增加重构相空间中的信息[31].选用1950—2018 年冬季SAT 和NAM,运用多视角嵌入法对冬东北亚地区冬季SAT 进行预测,并对比其预测能力.

首先仅用东北亚地区冬季SAT 进行预测,选择与上述单变量和多变量嵌入法相同的嵌入维数E=4,允许最大滞后时间步长l=4,预测步数tp=1,使用多视角嵌入法进行相空间重构共有即1325 种重构相空间.根据其预测值和观测值的相关系数从大到小进行排列,挑选出其中预测能力最强的前4 个重构相空间,并将预测值和观测值的相关系数进行比较,如图5 所示,其中k为嵌入的次数.可以看出当k取40—43 时,SAT 的预测效果最好,其预测能力最大可达0.6638.

图5 仅用东北亚地区冬季SAT 重构相空间,(a)—(d)分别表示k 取40,41,42,43 的预测结果Fig.5.Only using Northeast Asia winter surface air temperature reconstruct the phase space:(a)–(d) represent the prediction results of k taking 40,41,42,43 respectively.

同样地,将1950—2018 年11 月至次年2 月的NAM 信号加入SAT 的重构相空间中,选择嵌入维数E=4,允许最大滞后时间步长l=4,预测步数tp=1,使用多视角嵌入法整合所有与预测变量有关的增量进行相空间重构.图6 中仅选择了其中1 个月的NAM 加入SAT 的重构相空间中,共有3480 种重构相空间,挑选出其中预测能力最强的前4 个重构相空间的预测值,并将预测值和观测值的相关系数进行比较.可以看出,当k为30,36,31,34 时预测效果最好,预测值和真实值的相关系数最高可达0.687.

图6 仅将2 月的NAM 信号加入SAT 的重构相空间中,(a)—(d)分别表示k 取30,36,31,34 的预测结果Fig.6.Only adding the NAM signal in February to the reconstructed phase space of the SAT:(a)–(d) represent the prediction results of k taking 30,36,31,and 34,respectively.

图7 所示为将11 月到次年2 月的全部NAM信号加入SAT 的重构相空间中对SAT 的预测结果.在东北亚地区冬季温度的重构相空间中继续加入NAM 信号后,其可重构的相空间数量又有所增多,共有25334 种重构相空间.挑选出其中预测能力最强的前4 个重构相空间的预测值,并将预测值和观测值的相关系数进行比较,当k为38,39,36,37 时对SAT 的预测效果最好,预测值和真实值的相关系数最高可达0.7946.

图7 将12 月至次年2 月的NAM 信号加入SAT 的重构相空间中,(a)—(d)分别表示当k 取38,39,36,37 时的预测结果Fig.7.Adding the NAM signals from December to next February to the reconstructed phase space of the SAT:(a)–(d) show the prediction results when k takes 38,39,36,and 37,respectively.

可以看出,使用多视角嵌入法进行预测时,利用复杂系统中变量中共有信息的特性,在时间序列长度一定的情况下,利用冬季动力系统的复杂性来增加系统内的信息,将NAM 加入SAT 的重构相空间中可以改善对SAT 的预测效果.

4.4 三种重构相空间方法的比较

表1 对比了单变量嵌入、多变量嵌入、多视角嵌入三种重构相空间方法对1950—2018 年2 月份SAT 的预测能力(ρ)、平均绝对误差(MAE)、均方根误差(RMSE),得到三种嵌入方法的预测能力如下:多视角嵌入法 >多变量嵌入法 >单变量嵌入法;使用三种嵌入法预测的平均绝对误差和均方根误差如下:单变量嵌入法>多变量嵌入法>多视角嵌入法.经上述对比可知,使用多视角嵌入法对变量的预测能力最强,其对应的均方根误差和平均绝对误差最小,表明对于已知长度的时间序列,重构相空间中包含的有用信息越多,预测效果越好.此外,对比仅用SAT 重构相空间和将NAM 加入SAT 重构相空间的预测能力可知,将NAM 加入SAT 的重构相空间中可以改善SAT 的预测效果.

表1 多种嵌入方法的预测结果Table 1.Predicted results of multiple embedding methods.

此外,从理想数据结果和实际应用的结果可以看出,就单变量嵌入法而言,系统内部的所有信息最终都反映在单一变量时间序列的全部演化过程中,对于变量较少且系统中能量较大的分量所产生的作用很快就会体现在时间序列的短期变化上.而在实际的非平稳系统中,系统的高阶项或较弱分量产生的影响在短时间无法反映出来,几乎不可能从长度有限的时间序列中提取系统的全部特征.多视角嵌入法通过整合多种与预测变量有关的变量来增加重构相空间中的信息,从而恢复系统的动力学特征,该种方法包含的信息最多,因此在实际应用中,多视角嵌入法的预测结果要明显优于其他两种嵌入方法.

5 结论

将CCM 算法应用于非线性系统的预测,是因果分析及其预测应用上研究进展之一.同时,它也加深了我们对非线性的气候系统内部规律的认识,使进一步探究系统的驱动机理成为可能.本文通过两组预测试验,讨论了结合 CCM 方法,探讨加入驱动力因素对预测结果的影响并比较3 种不同预测建模方法,并得到如下结果.

1)对于理想Lorenz 模型和实际气候时间序列,单变量、多变量和多视角嵌入法三重构相空间方法对变量的预测能力均为:多视角嵌入法 >多变量嵌入法 >单变量嵌入法.利用多视角嵌入法对变量的预测效果最好,表明对于给定长度的时间序列,重构相空间中包含的信息越多,其预测能力越强.

2)将NAM 加入SAT 的重构相空间中有助于改善SAT 的预测效果.在气候系统预测中,可以考虑利用复杂系统中变量共有信息的特性,构建包含多种影响因子的东亚气温预测模型.

由于在CCM 算法中相空间中相互靠近的点在时间上具有相似的发展趋势和运动轨迹,CCM算法在理想模型和东北亚地区冬季温度预测上的应用,进一步展示了基于因果检验在非线性系统的预测能力.然而,文中的试验仍是初步的.本文并没有考虑变量之间的相互作用及其物理过程,此类问题的深入和扩展将是未来工作的方向.同时,实际的时间序列中含有噪声,解决预测模型对噪声的抗干扰能力也是一个重要课题.此外,预测中对于挖掘数据中定量信息的提取,可能对观测数据和质量提出新的要求.