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

潘昕濃 王革麗 楊培才

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