利用慢特征分析法提取层次结构系统中的外强迫∗

潘昕浓 王革丽 杨培才

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

利用慢特征分析法提取层次结构系统中的外强迫∗

潘昕浓1)2)王革丽1)†杨培才1)

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

(2016年10月8日收到;2017年1月16日收到修改稿)

在大量真实的动力系统中,外部驱动力总是随时间发生变化,正是这种变化导致了非平稳行为的产生.因此,从此类系统的观测数据中提取和分析外强迫(也称驱动力)信号引起了人们越来越多的关注.慢特征分析法(slow feature analysis,SFA)是从非平稳时间序列中提取外强迫信息的一种有效算法.在其基础上利用变参数的Logistic映射产生的非平稳时间序列,通过数值试验进一步讨论了该方法的应用前景,并发展了一些相应的分析技术.试验结果表明,对于模型中包含两个时变驱动力参数的系统,经过一次SFA处理之后,可以进一步利用子波分析技术检索出外强迫信号中的两个参数;对于模型中有两个叠加驱动力层次的三层动力系统,可先通过一次SFA处理,提取出次慢层外强迫信号,对该信号进行二次SFA处理,可提取出最慢层外强迫信号.

非平稳系统,非线性系统,驱动力,慢特征分析

1 引言

长期以来,在绝大多数非线性时间序列分析或预报中,几乎都暗含着一个缺省的假定,即被分析的时间序列满足遍历性定理.在这样的假定下,人们可以重建系统的动力学[1,2],并在此基础上讨论系统吸引子的几何结构和可预报性等[3].遍历性定理的成立,意味着时间序列来自于一个驱动力不随时间变化的系统,也就是说,被分析或预报的系统是平稳的.

然而,对于大多数实际的动力系统来说,其驱动力不可能一成不变.实际上,人们在实际观测中已经发现了此类系统的存在[4].驱动力的变化不仅破坏了系统的平稳性,而且破坏了当前时间序列分析和预报的理论基础.因此,有关平稳性问题的研究已经引起了人们的广泛关注[5−12].

在此背景下,一些气候学、生物学和经济学领域中的科学家开始转向非平稳问题的研究,并且取得了许多重要的进展[13−16].特别值得提出的是,在21世纪初Verdes等[17]和W iskott[18−20]分别提出了不同的从非平稳时间序列中提取外强迫因子的方法.Verdes等的方法建立在交叉预报误差的基础上,利用系统在两个不同时段动力学规律的变化,反演由此引起的外强迫因子的变化.W iskott则是将观测信号在一个给定的函数空间展开,并借助奇异谱分析(singular spectruManalysis,SSA)技术,找出信号变化最慢的分量.由他们的方法重建的外强迫与真实驱动力的值相比,都只差一个振幅因子和一个平移因子.从理论上讲,他们的方法都只能给出一个单一的外强迫因子,但是这个因子究竟包含哪些信息,这些方法是否还有提供更多外强迫信息的潜在能力,都是值得进一步讨论的问题.

Konen和Koch[21]曾经在一些数值试验中发现,在不同的情况下,W iskott[18−20]的方法可能只给出变化最慢的一个驱动力因子分量,也可能给出所有分量的一个组合.这些问题的解决在实际应用中是十分重要的,如果被提取的外强迫来自一个拥有更多层次或更多参数的系统,那么慢特征分析(slow feature analysis,SFA)算法是否可以被推广到这些更为复杂的情况,都值得深入研究[22].

本文将借助一些常见的非线性理论模式[23],利用SFA方法和子波分析方法,从这些模式所产生的更复杂的非平稳信号中重建其包含的外强迫,破解和讨论SFA方法在推广和应用中的一些问题.

2 慢特征分析算法

慢特征分析方法的主要目的,是从一个给定的非平稳信号中提取它所包含的变化最慢的分量.快变和慢变意味着在信号分量之间存在尺度和能量上的差异.

Kärner[13]曾这样描述非平稳信号:在这样的系统中,小尺度分量是非平稳的,而大尺度分量是平稳的.杨培才和周秀骥[16]则用层次概念描述非平稳系统的动力学结构,他们认为非平稳系统可以看作由多个处在不同层次上的子系统串级而成.处在较高层次上的子系统是一个慢变过程,而处在低层次上的子系统是一个快变过程.由于不同的子系统在尺度上和能量上存在巨大差异(前者远大于后者),所以它们之间的相互作用是不对称的,前者控制着后者.上述系统的定义可以在数学上描述为

式中t表示时间,α(t)和xi(t)(i=1,2,···,m)分别代表高层次分量和低层次分量,并且满足

式中符号⟨·⟩表示对时间的平均.(2)式表示高层次分量α(t)是一个变化最慢的信号,(3)式则表示α(t)单向地控制或影响着所有低层次信号分量xi(t).

SFA算法的思想建立在(1)式的基础上,其目标是利用(1)式给出的一个时间序列,反演出它的高层次分量α(t).在W iskott提出的算法中,(1)式中每一个分量都假定为多元二次多项式函数,则提取变化最慢的信号分量的问题,可转化为寻找满足(2)式条件的多项式函数的变分优化问题.

假定已知的非平稳时间序列为{x(t)}t=t1,t2,···,tn(n表示序列的长度),则SFA的计算步骤可简要描述如下.

将{x(t)}t=t1,t2,···,tn嵌入一个嵌入维数为m,时滞参数τ=1的状态空间:

式中N=n−m+1.

利用(4)式所产生的全部一阶和二阶单项式,构建一个k维函数空间:

为了方便起见,将(5)式写为

式中k=m+m(m+1)/2,继续对(6)式实施归一化和正交化,并将得到的结果记为

且Z(t)的每个分量都满足¯zj=0,zj=1,j=1,2,···,k.至此,(7)式中每一个输出信号都可以表示为zj的线性组合:

式中(a1,a2,···,ak)是一组待定的非零权重函数.

式中r和c是任意常数,分别为求解特征向量W1和积分导函数(9)式所产生.

至此已经简要地说明了高层次外强迫的提取步骤.我们将给出一个由上述方法得到的算例.

考虑一个时变的Logistic映射

式中控制参数at=C cos(2πt/T1)exp(−t/T2)+B是一个缓慢变化的衰减振荡(图1).出现于其中的常数分别为C=−0.2,B=3.75,T1=160,T2=200.这组常数所决定的at使得(11)式处于混沌体制之下.令(11)式迭代1000次,截取最后500次的数据段,将其作为试验时间序列(图1).将SFA算法应用于此试验序列,并令嵌入维数m=6,时滞参数τ=1,可以获得隐藏于其中的外强迫信号.当此外强迫与真实驱动力参数at都进行标准化后,它们之间的相关系数可以达到0.998(图1).

图1 (a)Logistic映射产生的非平稳时间序列{x t};(b)SFA方法提取得到的外强迫(红线)及真实驱动力{a t}(黑线)Fig.1.(a)The non-stationary tiMe series{x t}generated by logistic Mapping;(b)the d riving force extracted by SFA(red line)and the true d riving force{a t}(b lack line).

我们将利用几个双参数Logistic映射所产生的非平稳时间序列,来讨论SFA方法的某些推广应用问题.

3 数值试验

考虑双参数Logistic映射

式中a和b为两个驱动力参数.我们可以将其改造成经典的Logistic映射的形式.令yt=xt/b,代入(12)式,可以得到

借助于经典的Logistic映射与控制参数之间的依赖关系,可以得到(13)式的解在参数空间中的分布:当ab∈(0,1)时,系统有一个稳定的零解;当ab∈(1,3)时,系统为非零的稳定平衡态;继续下去可有:ab∈(3,3.5699),为稳定周期态,系统产生倍周期分叉结构;ab∈(3.5699,4),系统进入混沌区.图2为上述结果的直观表示,也代表双参数Logistic映射(12)式的解在参数空间中的分布.当参数a和b依赖于时间时,我们还可以从中看到变参数Logistic映射的状态随时间的大致变化.

图2 双参数Logistic映射的解在参数空间中的分布Fig.2.The solu tion distribution of logistic Mapp ing w ith two paraMeters in paraMetric space.

3.1 双驱动力参数Logistic映射

考虑带有两个时变参数的Logistic映射:

设定参数遵从

它们分别为一个慢变信号和一个快变信号.可以得到,它们满足条件3.02

令初值x1=0.6,并对(14)式迭代15000次,取后2000个数据作为试验信号{xt}.令嵌入维数m=6,时滞参数τ=1,对{xt}进行SFA处理,得到外强迫信号{yt},各信号的变化如图3所示.

图3 (a)真实驱动力{a t};(b)真实驱动力{b t};(c)非平稳时间序列{x t};(d)当嵌入维数m=6时SFA给出的外强迫信号{y t}Fig.3.(a)The true d riving force{a t};(b)the true d riving force{b t};(c)the non-stationary tiMe series{x t};(d)the derived d riving force for{x t}by using SFA(m=6).

图4 (a){y t}的时间平均功率谱(黑色实线)及其95%的信度检验(红色实点),功率谱上的两个峰值W ave1和W ave2(蓝色实点);(b)标准化的真实驱动力信号{a t}(蓝线)和W ave2的滤波信号{w2t}(绿线);(c)标准化的真实信号{b t}(蓝线)和W ave1的滤波信号{w 1t}(绿线)Fig.4.(a)The tiMe-averaged power spectruMof{y t}(b lack line)and 95%con fidence level(red dots),two peak powers W ave1 and W ave2(b lue points);(b)the norMalized true signal{a t}(b lue line)and the band-pass fi ltered signal{w2t}(green line);(c)the norMalized true signal{b t}(b lue line)and the band-pass fi ltered signal{w 1t}(green line).

由图3可见,{yt}的变化规律既不与{at}相同,也和{bt}不一致,应该是{at}和{bt}的组合.利用子波分析,我们得到了{yt}的时间平均功率谱,并发现两个峰值功率(分别记为Wave1和Wave2).利用带通滤波技术,我们又进一步得到了与峰值功率对应的滤波信号(分别记为{w1t}和{w2t}),如图4所示.

进一步分析表明,滤波信号{w1t}和{w2t}与真实驱动力信号{bt}和{at}高度相似.它们之间的相关系数皆达到0.99,相关系数记为|C(at,w2t)|=0.99,|C(bt,w1t)|=0.99.这说明子波分析可以成功地分离外强迫信号.

3.2 具有多层驱动力结构的Logistic映射

考虑层次结构,设置一个时变Logistic映射的三层模型:

式中{bt}为次慢层变化,{rt}为最慢层变化.设a=15,t=1,2,···,N,两个初值x1=0.5,b1=0.5.在此模型中,3.5 6 rt6 4.1,在rt的控制下,bt的变化在(3.6,3.97)区间内.这种情况下,(17)式处于混沌体制内.将(17)式迭代2100次,取后2000个数据作为试验信号,各参数变化如图5所示,最慢层信号{rt}比次慢层信号{bt}慢4倍,{xt}呈现非平稳特征.

图5 (a)真实驱动力{r t};(b)真实驱动力{b t};(c)非平稳时间序列{x t}Fig.5.(a)The true d riving force{r t};(b)the true d riving force{b t};(c)the non-stationary tiMe series{x t}.

取嵌入维数m=6,时滞参数τ=1,对{xt}进行SFA处理,得到信号{y1t},发现其与次慢层驱动力信号{bt}相关度很高,相关系数|C(y1t,bt)|=0.99;取m=36,再对{y1t}进行SFA处理,得到的信号{y2t}与最慢层信号{rt}相关度很高,相关系数|C(y2t,rt)|=0.96,即可以通过两次应用SFA的方法分层次提取出{xt}中的外强迫信号.标准化后的各变量如图6所示.结果表明,通过两次应用SFA,可以分别提取出次慢层外强迫及最慢层外强迫信号.

图6 (a)嵌入维数m=6时{x t}的标准化外强迫信号{y1t}(蓝线)以及标准化后的真实信号{b t}(绿线);(b)m=36时{y1t}的标准化外强迫信号{y2t}(蓝线)以及标准化后的真实信号{r t}Fig.6.(a)The derived d riving force{y1t}for{x t}(b lue line)and the true d riving force{b t}(green line)as m=6;(b)the ex tracted signal{y2t}for{y1t}(b lue line)and the true d riving force{r t}(green line)as m=36.

4 结论

本文通过两组试验说明,SFA方法结合子波分析技术,可以从更复杂的非平稳时间序列中提取系统的外强迫信号,并探测它们的结构.需要注意的是,本文构建了非平稳(非自治)系统,所以传统的基于遍历性定理的求取嵌入维数m与时滞参数τ的方法不再适用,需要依靠经验与试探来选取m与τ,试验者可以通过功率谱分析对提取的结果进行物理性质的判断,进而选取合理的参数.本文的主要结论如下:1)对于双变参数的Logistic映射产生的非平稳时间序列,可以利用SFA方法提取模型中的外强迫信息,再结合功率谱分析的方法,重建真实驱动力信号;2)对于隐含最慢层外强迫的三层结构Logistic映射,可以通过逐次应用SFA的方法提取系统中的外强迫信号,第一次应用SFA可以提取次慢层外强迫,第二次应用SFA可以提取最慢层外强迫.

[1]Packard N H,Cru tch field J P,FarMer J D,Shaw R S 1980 Phys.Rev.Lett.45 712

[2]David R,Takens F 1971 ComMun.Math.Phys.20 167

[3]Bian J C,Yang P C 2003 P lateau Meteor.22 315(in Chinese)[卞建春,杨培才2003高原气象22 315]

[4]Tsonis A A 1996 Nature 382 700

[5]W ang SW,Zhu J H 2000 Quart.J.Appl.Meteor.11 1(in Chinese)[王绍武,朱锦红2000应用气象学报11 1]

[6]W ang H J,Zhou G Q,Lin Z H 2002 CliMatic Environ.Res.7 220(in Chinese)[王会军,周广庆,林朝晖2002气候与环境研究7 220]

[7]D ing Y H,Ren G Y,Zhao Z C,Xu Y,Luo Y,Li Q P,Zhang J 2007 Desert Oasis Meteor.1 1(in Chinese)[丁一汇,任国玉,赵宗慈,徐影,罗勇,李巧萍,张锦2007沙漠与绿洲气象1 1]

[8]W ang SW 1998 Adv.Earth Sci.13 8(in Chinese)[王绍武1998地球科学进展13 8]

[9]Chen B M,Ji L R,Yang P C,Zhang D M,W ang G L 2003 Chin.Sci.Bu ll.48 513(in Chinese)[陈伯民,纪立人,杨培才,张道民,王革丽2003科学通报48 513]

[10]W an SQ,Feng G L,Dong W J,Li J P 2005 Acta.Phys.Sin.54 5487(in Chinese)[万仕全,封国林,董文杰,李建平2005物理学报54 5487]

[11]W ang G L,Yang P C,Zhou X J 2016 Theor.Appl.C liMatol.124 985

[12]Chen X X,W ang G L,Jin L J 2015 China Environ.Sci.35 694(in Chinese)[陈潇潇,王革丽,金莲姬2015中国环境科学35 694]

[13]Kärner O 2002 J.Geophs.Res.107 ACL1-1

[14]Davis A,Marshak A,W iscombeW,Cahalan R 1996 J.A tMos.Sci.53 1538

[15]Yang P C,Bian J C,W ang G L,Zhou X J 2003 Chin.Sci.Bu ll.48 1470(in Chinese)[杨培才,卞建春,王革丽,周秀骥2003科学通报48 1470]

[16]Yang P C,Zhou X J 2005 Acta.Meteor.Sin ica 63 556(in Chinese)[杨培才,周秀骥2005气象学报63 556]

[17]Verdes P F,G ranitto P M,Navone H D,Ceccatto H A 2001 Phys.Rev.Lett.87 124101

[18]W iskott L 2003 Neural CoMput.15 2147

[19]W iskott L,Sejnowski T J 2002 Neural CoMput.14 715

[20]W iskott L http://www.arxiv.org/abs/cond-Mat/0312317[2016-9-21]

[21]Konen W,Koch P http://arxiv.org/abs/0911.4397v1[2016-9-21]

[22]Pan X N,W ang G L,W ang P F,Zhu K Y 2016 Meteor.Environ.Sci.39 96(in Chinese)[潘昕浓,王革丽,王鹏飞,朱克云2016气象与环境科学39 96]

[23]Robert M1976 Nature 261 459

(Received 8 October 2016;revised Manuscrip t received 16 January 2017)

PACS:05.45.–aDOI:10.7498/aps.66.080501

*Pro ject supported by the National Natural Science Foundation of China(G rant No.41575058).

†Corresponding author.E-Mail:wgl@Mail.iap.ac.cn

Ex tracting the d riving force signal froMh ierarchy systeMbased on slow featu re analysis∗

Pan Xin-Nong1)2)Wang Ge-Li1)†Yang Pei-Cai1)

1)(K ey Laboratory ofMidd le A tMosphere and G lobal EnvironMent Observation,Institute of A tMospheric Physics,Chinese AcadeMy of Sciences,Beijing 100029,China)2)(University of Chinese AcadeMy of Sciences,Beijing 100049,China)

Extracting the signals froMnon-stationary time series isa diffi cu lt task inmany fieldssuch as physics,econoMics,and atMospheric sciences.The theory of hierarchy suggests that varying d riving force leads to the non-stationary behavior,so extracting and analyzing the slow ly varying features can help to study non-stationary dynaMical system,which has become a coMpelling question recently.Slow feature analysis(SFA)is an eff ective technique for extracting slow ly varying d riving forces froMquick ly varying non-stationary tiMe series.The basic idea of SFA is to nonlinearly extend the reconstructive signal into a combination forMw ith one or higher order polynoMials,and to app ly the p rincipal coMponent analysis to this extended signal and its time derivatives.The algorithMis guaranteed to seek an op timal solution froMa group of functions directly and can extract a lot of uncorrelated features that are ordered by slowness.A series of studies has shown its superiority in extracting the driving force of non-stationary tiMe series.The extracted signal is found to be highly correlated w ith the real driving force.Results based on idealmodels show that either the slow driving force itself or a slower subcoMponent can be detected by SFA.Yet despite all that,the further investigating of SFA is still needed to reduce its uncertainty.In this study,we create two types of non-stationary models by the logistic map w ith tiMe-varying paraMeters:one includes two varying driving forces w ith diff erent tiMe periods constraining the evolution of tiMe series in a non-stationary way;and the other is a three-layer structure encoMpassing two superiMposed signals in which the slower signalof d riving force ismodulated by the lowest one.According to the idealmodeland SFA,we conduct the nuMerical experiMents to develop corresponding analysisMethod and discuss its app lication prospect in extracting driving force signals.W e find that for the systeMof fi rst kind,either the slowest signal or the combination of two d riving forces constructed by SFA contains some uncertain information.However,we can detect the two independent d riving forces froMthe constructed signal by wavelet analysis.For the three-hierarchy systeMthat includes two superiMposed signals of driving force,successive app lications through SFA on the original tiMe series and the constructed SFA signal w ill in turn detect the slower varying driving force signal and the slowest varying driving forces signal.The successful app lication of SFA show s its p roMising prospect in analyzing the external driving forces in non-stationary systeMand understanding relevant dynaMic Mechanism.

non-stationary system,nonlinear system,driving force,slow feature analysis

10.7498/aps.66.080501

∗国家自然科学基金(批准号:41575058)资助的课题.

†通信作者.E-Mail:w gl@Mail.iap.ac.cn

©2017中国物理学会C h inese P hysica l Society

http://w u lixb.iphy.ac.cn