利用慢特征分析法提取層次結構系統中的外強迫?

潘昕濃 王革麗 楊培才

1)(中國科學院大氣物理研究所,中層大氣和全球環(huán)境探測重點實驗室,北京100029)2)(中國科學院大學,北京100049)

利用慢特征分析法提取層次結構系統中的外強迫?

潘昕濃1)2)王革麗1)?楊培才1)

1)(中國科學院大氣物理研究所,中層大氣和全球環(huán)境探測重點實驗室,北京100029)2)(中國科學院大學,北京100049)

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

在大量真實的動力系統中,外部驅動力總是隨時間發(fā)生變化,正是這種變化導致了非平穩(wěn)行為的產生.因此,從此類系統的觀測數據中提取和分析外強迫(也稱驅動力)信號引起了人們越來越多的關注.慢特征分析法(slow feature analysis,SFA)是從非平穩(wěn)時間序列中提取外強迫信息的一種有效算法.在其基礎上利用變參數的Logistic映射產生的非平穩(wěn)時間序列,通過數值試驗進一步討論了該方法的應用前景,并發(fā)展了一些相應的分析技術.試驗結果表明,對于模型中包含兩個時變驅動力參數的系統,經過一次SFA處理之后,可以進一步利用子波分析技術檢索出外強迫信號中的兩個參數;對于模型中有兩個疊加驅動力層次的三層動力系統,可先通過一次SFA處理,提取出次慢層外強迫信號,對該信號進行二次SFA處理,可提取出最慢層外強迫信號.

非平穩(wěn)系統,非線性系統,驅動力,慢特征分析

1 引言

長期以來,在絕大多數非線性時間序列分析或預報中,幾乎都暗含著一個缺省的假定,即被分析的時間序列滿足遍歷性定理.在這樣的假定下,人們可以重建系統的動力學[1,2],并在此基礎上討論系統吸引子的幾何結構和可預報性等[3].遍歷性定理的成立,意味著時間序列來自于一個驅動力不隨時間變化的系統,也就是說,被分析或預報的系統是平穩(wěn)的.

然而,對于大多數實際的動力系統來說,其驅動力不可能一成不變.實際上,人們在實際觀測中已經發(fā)現了此類系統的存在[4].驅動力的變化不僅破壞了系統的平穩(wěn)性,而且破壞了當前時間序列分析和預報的理論基礎.因此,有關平穩(wěn)性問題的研究已經引起了人們的廣泛關注[5?12].

在此背景下,一些氣候學、生物學和經濟學領域中的科學家開始轉向非平穩(wěn)問題的研究,并且取得了許多重要的進展[13?16].特別值得提出的是,在21世紀初Verdes等[17]和W iskott[18?20]分別提出了不同的從非平穩(wěn)時間序列中提取外強迫因子的方法.Verdes等的方法建立在交叉預報誤差的基礎上,利用系統在兩個不同時段動力學規(guī)律的變化,反演由此引起的外強迫因子的變化.W iskott則是將觀測信號在一個給定的函數空間展開,并借助奇異譜分析(singular spectruManalysis,SSA)技術,找出信號變化最慢的分量.由他們的方法重建的外強迫與真實驅動力的值相比,都只差一個振幅因子和一個平移因子.從理論上講,他們的方法都只能給出一個單一的外強迫因子,但是這個因子究竟包含哪些信息,這些方法是否還有提供更多外強迫信息的潛在能力,都是值得進一步討論的問題.

Konen和Koch[21]曾經在一些數值試驗中發(fā)現,在不同的情況下,W iskott[18?20]的方法可能只給出變化最慢的一個驅動力因子分量,也可能給出所有分量的一個組合.這些問題的解決在實際應用中是十分重要的,如果被提取的外強迫來自一個擁有更多層次或更多參數的系統,那么慢特征分析(slow feature analysis,SFA)算法是否可以被推廣到這些更為復雜的情況,都值得深入研究[22].

本文將借助一些常見的非線性理論模式[23],利用SFA方法和子波分析方法,從這些模式所產生的更復雜的非平穩(wěn)信號中重建其包含的外強迫,破解和討論SFA方法在推廣和應用中的一些問題.

2 慢特征分析算法

慢特征分析方法的主要目的,是從一個給定的非平穩(wěn)信號中提取它所包含的變化最慢的分量.快變和慢變意味著在信號分量之間存在尺度和能量上的差異.

K?rner[13]曾這樣描述非平穩(wěn)信號:在這樣的系統中,小尺度分量是非平穩(wěn)的,而大尺度分量是平穩(wěn)的.楊培才和周秀驥[16]則用層次概念描述非平穩(wěn)系統的動力學結構,他們認為非平穩(wěn)系統可以看作由多個處在不同層次上的子系統串級而成.處在較高層次上的子系統是一個慢變過程,而處在低層次上的子系統是一個快變過程.由于不同的子系統在尺度上和能量上存在巨大差異(前者遠大于后者),所以它們之間的相互作用是不對稱的,前者控制著后者.上述系統的定義可以在數學上描述為

式中t表示時間,α(t)和xi(t)(i=1,2,···,m)分別代表高層次分量和低層次分量,并且滿足

式中符號?·?表示對時間的平均.(2)式表示高層次分量α(t)是一個變化最慢的信號,(3)式則表示α(t)單向地控制或影響著所有低層次信號分量xi(t).

SFA算法的思想建立在(1)式的基礎上,其目標是利用(1)式給出的一個時間序列,反演出它的高層次分量α(t).在W iskott提出的算法中,(1)式中每一個分量都假定為多元二次多項式函數,則提取變化最慢的信號分量的問題,可轉化為尋找滿足(2)式條件的多項式函數的變分優(yōu)化問題.

假定已知的非平穩(wěn)時間序列為{x(t)}t=t1,t2,···,tn(n表示序列的長度),則SFA的計算步驟可簡要描述如下.

將{x(t)}t=t1,t2,···,tn嵌入一個嵌入維數為m,時滯參數τ=1的狀態(tài)空間:

式中N=n?m+1.

利用(4)式所產生的全部一階和二階單項式,構建一個k維函數空間:

為了方便起見,將(5)式寫為

式中k=m+m(m+1)/2,繼續(xù)對(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映射產生的非平穩(wěn)時間序列{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映射所產生的非平穩(wěn)時間序列,來討論SFA方法的某些推廣應用問題.

3 數值試驗

考慮雙參數Logistic映射

式中a和b為兩個驅動力參數.我們可以將其改造成經典的Logistic映射的形式.令yt=xt/b,代入(12)式,可以得到

借助于經典的Logistic映射與控制參數之間的依賴關系,可以得到(13)式的解在參數空間中的分布:當ab∈(0,1)時,系統有一個穩(wěn)定的零解;當ab∈(1,3)時,系統為非零的穩(wěn)定平衡態(tài);繼續(xù)下去可有:ab∈(3,3.5699),為穩(wěn)定周期態(tài),系統產生倍周期分叉結構;ab∈(3.5699,4),系統進入混沌區(qū).圖2為上述結果的直觀表示,也代表雙參數Logistic映射(12)式的解在參數空間中的分布.當參數a和b依賴于時間時,我們還可以從中看到變參數Logistic映射的狀態(tài)隨時間的大致變化.

圖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)非平穩(wěn)時間序列{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}的變化規(guī)律既不與{at}相同,也和{bt}不一致,應該是{at}和{bt}的組合.利用子波分析,我們得到了{yt}的時間平均功率譜,并發(fā)現兩個峰值功率(分別記為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)區(qū)間內.這種情況下,(17)式處于混沌體制內.將(17)式迭代2100次,取后2000個數據作為試驗信號,各參數變化如圖5所示,最慢層信號{rt}比次慢層信號{bt}慢4倍,{xt}呈現非平穩(wěn)特征.

圖5 (a)真實驅動力{r t};(b)真實驅動力{b t};(c)非平穩(wěn)時間序列{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},發(fā)現其與次慢層驅動力信號{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方法結合子波分析技術,可以從更復雜的非平穩(wěn)時間序列中提取系統的外強迫信號,并探測它們的結構.需要注意的是,本文構建了非平穩(wěn)(非自治)系統,所以傳統的基于遍歷性定理的求取嵌入維數m與時滯參數τ的方法不再適用,需要依靠經驗與試探來選取m與τ,試驗者可以通過功率譜分析對提取的結果進行物理性質的判斷,進而選取合理的參數.本文的主要結論如下:1)對于雙變參數的Logistic映射產生的非平穩(wěn)時間序列,可以利用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氣候與環(huán)境研究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中國環(huán)境科學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氣象與環(huán)境科學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