基于條件熵擴維的多變量混沌時間序列相空間重構*

張春濤 馬千里 彭 宏 姜友誼

1)(重慶三峽學院數(shù)學與計算機科學學院，萬州 404000)

2)(華南理工大學計算機科學與工程學院，廣州 510006)

基于條件熵擴維的多變量混沌時間序列相空間重構*

張春濤馬千里2)彭 宏2)姜友誼1)

1)(重慶三峽學院數(shù)學與計算機科學學院，萬州 404000)

2)(華南理工大學計算機科學與工程學院，廣州 510006)

(2010年4月12日收到;2010年5月22日收到修改稿)

提出一種多變量混沌時間序列相空間重構的條件熵擴維方法.首先使用互信息法求解每個變量的時間延遲，其次按條件熵最大原則逐步擴展相空間的嵌入維數(shù)，使得重構坐標從低維到高維的轉換保持較強的獨立性，最終的重構相空間具有較低的冗余度，為多變量時間序列的預測構造了有效的模型輸入向量.通過對幾個經(jīng)典多變量混沌時間序列進行數(shù)值實驗，結果表明該方法比單變量預測和已有多變量預測方法具有更好的預測效果，說明了該重構方法的有效性.

多變量混沌時間序列，相空間重構，條件熵，神經(jīng)網(wǎng)絡預測

PACS:05.45.Pq，05.45.Tp

1.引 言

復雜系統(tǒng)中任一分量的演化都可由系統(tǒng)的其他分量所決定，因此任一分量的發(fā)展過程都包含其他分量的信息，即通過該分量發(fā)展的時間序列重構相空間能充分反映原系統(tǒng)的動力學特征，這就是單變量混沌時間序列相空間重構［1—8］.在現(xiàn)實狀況下由于無法從單變量的時間序列獲全系統(tǒng)信息，以及在分析生物學系統(tǒng)時，由于其結構和動力系統(tǒng)的復雜性，很難由一維時間序列通過延遲嵌入法重構出相空間［9］，因此不能保證實際系統(tǒng)中任何給定的單變量時間序列足以重構原系統(tǒng)［10，11］.而多變量時間序列由于包含了更豐富完整的系統(tǒng)信息，故能重構出更為準確的相空間，因此多變量混沌時間序列的研究得到越來越多的關注［12—18］.

目前，對多變量時間序列已取得了許多研究成果.王海燕等［14］提出由單變量推廣的虛假鄰近點法計算多變量時間序列的相空間嵌入維數(shù)，但其缺少每次擴維選擇標準，使重構相空間具有不定性. Broomhead等［15］提出奇異值分解法求取多變量時間序列的嵌入維數(shù).Cao等［10］對多變量的相空間重構問題進行理論和應用研究，并論證了其優(yōu)越性;樊重?。?6］討論了多觀測變量狀態(tài)空間重構技術，為研究多變量時間序列的相空間重構提供了理論基礎.

信息熵是系統(tǒng)整體不確定性的一種量度，在混沌時間序列相空間重構中有重要應用.田玉楚等［19］系統(tǒng)地論述了在不變測度意義下，將混沌系統(tǒng)視為一類特殊的隨機系統(tǒng)加以處理;肖方紅等［6］對一維時間序列使用條件熵求取相空間的嵌入維數(shù);我們也研究了兩個重構參數(shù)的信息熵模型［8］.這些研究都集中在單變量混沌時間序列上，對多變量混沌系統(tǒng)，主要用互信息法求解單個變量的延遲時間.

本文沿用單變量時間序列建模的思想，結合信息熵和神經(jīng)網(wǎng)絡對多變量時間序列進行相空間重構和預測，提出了多變量混沌時間序列的條件熵擴維方法(CEED).首先使用互信息法求解單個變量的時間延遲;然后使用條件熵擴維方法，按條件熵最大原則選擇重構分量，把低維相空間逐漸擴展到高維相空間，使得最終的相空間具有較少的冗余度.為了檢驗該重構方法的有效性，最后使用神經(jīng)網(wǎng)絡對2個典型的多變量混沌系統(tǒng)進行仿真實驗，數(shù)值實驗表明在該重構方法上進行的多變量混沌時間序列預測能取得比單變量和已有的多變量預測更好的效果，說明CEED方法重構了高質量的相空間.

2.多變量混沌時間序列相空間重構

在一定程度上單變量時間序列是多變量時間序列的特殊情況，因此可將單變量的延遲重構應用到多變量中.對單變量時間序列x(n)(n=1，2，…，N)，在確定嵌入維數(shù) m和時間延遲 τ后，重構相空間

設復雜系統(tǒng)的N維多觀測變量序列為 X1，X2，…，XK，其中 K表示觀測變量個數(shù)，Xi=［xi(1)，xi(2)，…，xi(N)］表示第i個變量的時間序列.如果每個變量都選擇恰當?shù)臅r間延遲τi和嵌入維數(shù)mi(i=1，2，…，K)，則多變量時間序列重構的相空間

其中τi和mi是重構相空間成敗的關鍵，目前還沒有一個統(tǒng)一的選擇標準，下面使用信息論知識來選擇這K對重構參數(shù).

2.1.多變量時間序列延遲時間的選擇

設兩個離散型隨機變量X∈{x1，x2，…，xn}和Y∈{y1，y2，…，yk}，其 先 驗 概 率 分 別 為 {p (xi)}i=1，2，…n和{p(yj)}j=1，2，…，k.信息熵可定義為

它描述了隨機變量X的不定性.類似地可以定義聯(lián)合熵

其中p(xi，yj)為聯(lián)合概率.

變量X和Y之間的互信息定義為

根據(jù)(7)式可推出時間序列和其延遲時間序列的互信息

根據(jù)互信息法［5］，計算出互信息函數(shù) I(τi)在不同延遲τi下的值，取I(τi)的第一個局部極小值點對應的時間τi為第 i個變量時間序列 Xi的延遲時間.

2.2.多變量時間序列嵌入維數(shù)的確定

一般認為，如果重構相空間的維數(shù)足夠大，相空間就能夠刻畫出原系統(tǒng)的吸引子.對多變量時間序列，一個很自然的想法就是對每個變量的時間序列都分別求一個嵌入維數(shù)mi，然后將所有變量的重構相空間進行簡單迭加形成多變量時間序列的相空間如(2)式.由于混沌系統(tǒng)任一變量的時間序列包含其他變量的信息，即使單個變量的重構相空間無冗余信息，經(jīng)過簡單迭加而成的多變量時間序列的重構相空間也包含大量的冗余信息，不利于預測.如何減少這些冗余信息，是確定多變量時間序列嵌入維數(shù)的關鍵.主成分分析(PCA)方法［17，18］為我們提供了一種解決途徑，但是PCA方法只能消除各重構分量間的線性關系，對非線性的時間序列作用有限.下面根據(jù)條件熵能反映變量之間的獨立性的特征，通過條件熵消除各重構分量間的非線性關系，使重構相空間包含較少冗余信息，以達到降低相空間維數(shù)的目的.

記多變量重構相空間(2)式的延遲坐標向量的各分量為

相空間重構也就是用重構分量 X11，X12，…，XK1，XK2，…，XKmK近似表征原系統(tǒng)的各變量，同時X11，X12，…，XK1，XK2，…，XKmK也可以看成是重構相空間的坐標，因此原系統(tǒng)表現(xiàn)出來的混沌吸引子也將在重構分量X11，X12，…，XK1，XK2，…，XKmK的重構相空間中表現(xiàn)出來.這樣隨著重構分量的逐漸增多，重構向量就能越來越完整地表征原系統(tǒng)的動力學特征.當重構的分量個數(shù)增加到一定值時，這時重構的分量包含了原系統(tǒng)充分多的動力學信息，進一步增加重構分量不會增加原系統(tǒng)的動力學信息，只會攜帶冗余信息.因此可以設想分量X11，X21，…，XK1的熵H(X11，X21，…，XK1)和各分量之間的條件熵的值隨著重構分量個數(shù)的增加由大到小逐漸減小，直至重構分量包含充分多的原系統(tǒng)的動力學信息，進一步增加重構分量個數(shù)，這時條件熵將出現(xiàn)一個較明顯的轉折點，隨后的條件熵將趨于平穩(wěn).如果取條件熵出現(xiàn)較為明顯的轉折且趨于平穩(wěn)處的重構分量個數(shù)作為多變量時間序列重構相空間的嵌入維數(shù)，可以看出此時重構向量已包含原系統(tǒng)充分多的動力學信息，應該是一個合適的嵌入維數(shù)［6］.

多變量時間序列的重構向量中增加一個重構分量有諸多選擇，即在已知重構分量為 X11，X12，X21的條件下再增加一個分量，有 K種不同選擇:X13，X22，X31，X41，…，XK1.對這 K個條件熵做怎樣的取舍?由上所述，多變量時間序列重構相空間應具有較少的冗余信息，即重構分量之間應最大限度地相互獨立.因此增加一個重構分量的最佳選擇應是K個條件熵的最大值所對應的分量.因為條件熵最大，表明增加的這個分量和原有分量的獨立性最大，這樣不斷增加重構分量到整個相空間就能夠保證有最大的獨立性.下面給出從低維到高維的擴維(增加相空間重構分量個數(shù))算法步驟:

步驟1 對每個變量的時間序列使用互信息法求時間延遲τi，i=1，2，…，K;

步驟2 初始化:令Y=［X11，X21，…，XK1］為初始的重構向量，每個變量的初始嵌入維數(shù)mi=1，i= 1，2，…，K;

步驟3 計算條件熵并擴維:

3)擴維mk=mk+1，此時的重構向量 Y=［Y，Xkmk］;

4)重復1)—3)直至最大條件熵趨于平穩(wěn).

步驟4 循環(huán)結束后的重構向量空間就是多變量時間序列的重構相空間，此時的mi即為第i觀測變量的嵌入維數(shù)，重構相空間的嵌入維數(shù)m=m1+ m2+…+mK.

對條件熵按(6)式轉化為聯(lián)合熵來計算，本文采用矩陣標識法［6］計算聯(lián)合熵.

3.多變量混沌時間序列仿真實驗

為了驗證 CEED的有效性，我們做了 Henon，Lorenz系統(tǒng)的單步和多步預測實驗，并與現(xiàn)有的單變量和多變量預測方法進行比較.

3.1.預測模型的建立

假設重構的狀態(tài)空間如(2)式，相空間中預測的基本思想是:根據(jù)歷史狀態(tài)數(shù)據(jù)集建立第 i時刻狀態(tài)Y(i)到第i+h時刻狀態(tài)Y(i+h)的映射Fh，然后由映射Fh和當前n時刻狀態(tài)Y(n)去推斷未來n+h時刻狀態(tài)Y(n+h)的估計值Z(n+h)，其中h為預測步長.根據(jù)嵌入定理，狀態(tài)映射Fh是存在的，即在m維的重構相空間存在光滑函數(shù)fh:Rm→R.對多變量混沌時間序列，因此有預測模型

其中，zi(n+h)為模型對未來(n+h)時刻第i個變量狀態(tài)xi(n+h)的預測值(i=1，2，…，K)，h為預測步長.由于混沌系統(tǒng)的高度非線性，要確定函數(shù)fih的解析形式是困難的，我們采用神經(jīng)網(wǎng)絡的學習原理來逼近這些映射函數(shù).結合相空間重構和神經(jīng)網(wǎng)絡建立預測模型的結構如圖1所示.

圖1 多變量時間序列預測模型結構

本文采用的神經(jīng)網(wǎng)絡是 Elman遞歸神經(jīng)網(wǎng)絡［20］，它能夠通過內部的遞歸連接反映非線性動力系統(tǒng)的特征，在混沌時間序列的預測中有重要應用［2，3］.

3.2.仿真實驗

為了評價預測效果采用均方誤差(MSE)、均方根誤差(RMSE)、標準化均方誤差(NMSE)以及絕對誤差e(t)作為誤差評價標準:

3.2.1.Henon混沌時間序列

法國天文學家Henon于1976年提出一個簡單的二維平方映像如下:

通過計算得出變量x，y的3000個點的數(shù)據(jù)集，使用互信息法分別求解x和y的延遲時間如圖2和圖3所示，可知τ1=15，τ2=14，然后通過條件熵擴維如圖4所示.

由圖4可知，在嵌入維數(shù)為5時條件熵出現(xiàn)一個明顯的轉折，當嵌入維數(shù)增加到6后條件熵幾乎不發(fā)生變化，所以我們取Henon重構相空間的嵌入維數(shù)為6，再根據(jù)每次的重構分量選擇m=m1+m2=4+2=6.

通過Elman神經(jīng)網(wǎng)絡按(9)式對變量x進行預測.前2000個數(shù)據(jù)點進行訓練，后1000個數(shù)據(jù)進行預測，其中單步預測結果與單變量預測方法:在線最小二乘支持向量機(LS-SVM)［21］，最小二乘支持向量(LS-SVD)［22］，模糊樹Ⅰ(FTⅠ)［23］，F(xiàn)TⅡ［23］，模糊聚類［24］，基于自適應網(wǎng)絡的模糊推理系統(tǒng)(ANFIS)［25］，局域離散余弦變換領域Vloterra(LDCTDV)［26］，模糊邊界模塊化神經(jīng)網(wǎng)絡 FBMNN［2］以及多變量預測方法的比較如表1所示，其中非線性獨立分量分析(NICA)為文獻［27］所取的多變量相空間結構:m1=2，τ1=1;m2=2，τ2=1，然后使用Elman神經(jīng)網(wǎng)絡的預測.

圖2 x的互信息

圖3 y的互信息

圖4 條件熵與嵌入維數(shù)的關系

從表1可以看出CEED取得比單變量預測和多變量預測更好的預測效果.單步預測效果和預測誤差分別如圖5和圖6所示.

從圖6可以看出，對Henon混沌時間序列的單步預測誤差達到10-5數(shù)量級，比文獻［28］的10-2高3個數(shù)量級，說明使用多變量時間序列預測 Henon是高效.

表1 單步預測

圖5 單步預測效果

圖6 單步預測誤差

在進行的多步預測中，預測結果與正交最小二乘徑向基函數(shù)(RBF(OLS))［29］，局域線性模型樹(LoliMoT)［29］，F(xiàn)BMNN和NICA的比較如表2所示，其中節(jié)點數(shù)為隱含層節(jié)點數(shù).

表2 Henon預測結果的NMSE比較

由表2可以看出CEED使用較少的隱含層節(jié)點，從而降低了計算復雜性，卻能取得更好的預測效果.

從表1和表2可以看出，Henon映像是一個最簡單的二維混沌系統(tǒng)，使用單變量預測方法和已有的多變量相空間重構方法進行單步預測大多能取得較好的預測效果，但多步預測效果一般較差，而使用CEED的多變量預測方法在單步和多步預測上都取得了很好的預測效果.FBMNN，NICA和CEED都使用Elman進行預測，CEED上進行的預測取得了最好的預測效果，說明了 CEED重構方法對Henon映像的有效性.

3.2.2.洛倫茲方程

洛倫茲方程如下:

用四階Runge-Kutta方法求解方程(15)，取積分步長h=0.02，得到變量x，y的兩個7500個點的數(shù)值解，去除前5000個點的暫態(tài)過程，最后考慮兩個觀察變量x，y的2500個點的數(shù)據(jù)集，將其歸一化到［0，1］區(qū)間.先通過互信息法分別求出x和y的延遲時間如圖7和圖8所示，由圖7，圖8可知，τ1=8，τ2=8.

然后通過條件熵擴維如圖9所示，由圖9可知，當嵌入維數(shù)增加到6后條件熵不再發(fā)生變化，再根據(jù)擴維分量選擇有m=m1+m2=2+4=6.

圖7 x的互信息

圖8 y的互信息

圖9 條件熵與嵌入維數(shù)的關系

通過Elman神經(jīng)網(wǎng)絡按(9)式對變量x進行預測.前2000個數(shù)據(jù)點進行訓練，后500個數(shù)據(jù)進行預測，其中單步預測結果與單變量時間序列預測和多變量預測的比較如表3所示，其中IFNNMU為單變量推廣的虛假鄰近點法確定的相空間結構:m1= 2，τ1=8;m2=1，τ2=8，然后使用Elman神經(jīng)網(wǎng)絡的預測.

從表3可以看出，對Lorenz系統(tǒng)CEED能取得比已有的單變量時間序列預測和多變量預測更好的單步預測效果.

為了與多變量預測方法［30］進行比較，將歸一化數(shù)據(jù)還原為原始數(shù)據(jù)，根據(jù)eRMSE的定義有，對Lorenz系統(tǒng)歸一化數(shù)據(jù)后的 eRMSE是原始數(shù)據(jù)的 1/40，即CEED原始數(shù)據(jù)的 eRMSE≤2×10-5，比文獻［30］的0.063高3個數(shù)量級.CEED的嵌入維數(shù)是6，而文獻［30］的嵌入維數(shù)為11，這樣CEED能簡化神經(jīng)網(wǎng)絡結構，更利于實現(xiàn).

單步預測效果如圖10所示，預測誤差如圖11所示.

圖10 單步預測效果

圖11 單步預測誤差

同時我們也進行了多步預測實驗，預測步長取為6和8，與共進化遞歸神經(jīng)網(wǎng)絡(CERNN)［3］等相關預測方法的比較如表4和表5所示.

表4 6步預測

表5 8步預測

從表4和表5可以看出，對Lorenz系統(tǒng)，在CEED上進行的多步預測精度也高于單變量預測和多變量預測方法.從表 3—表 5可看出，F(xiàn)BMNN，CERNN，IFNNMU和CEED都使用了Elman網(wǎng)絡進行預測，在CEED重構相空間上進行的預測取得了更好的預測效果，說明CEED重構方法對Lorenz系統(tǒng)的有效性.

從表1—表5可以看出，對多維混沌系統(tǒng)，不管是與單變量時間序列預測方法，還是與已有的多變量預測方法相比，CEED重構相空間上進行單步和多步預測都能取得更好的預測效果，說明CEED重構的多維混沌系統(tǒng)的相空間是有效性.

4.結 論

本文提出一種多變量混沌時間序列相空間重構的新方法.該方法通過條件熵擴張相空間維數(shù)，按條件熵最大選擇重構分量，使得相空間從低維到高維的轉變過程中既保證重構坐標的獨立性又保持系統(tǒng)的動力學特征.因此通過CEED確定的重構相空間既能充分包含系統(tǒng)信息，又能有效避免冗余信息的產(chǎn)生.最后在CEED上進行多變量混沌時間序列的單步和多步預測，都取得了比使用單變量預測方法和已有的多變量預測方法更好的預測效果，表明CEED重構方法的有效性.

［1］Peng S G，Yu S M 2009 Chin.Phys.B 18 3758

［2］Ma Q L，Zheng Q L，Peng H，Qin J W 2009 Acta Phys.Sin.58 1410(in Chinese)［馬千里、鄭啟倫、彭 宏、覃姜維2009物理學報58 1410］

［3］Ma Q L，Zheng Q L，Peng H，Zhong T W，Qin J W 2008 Chin. Phys.B 17 536

［4］Zhang S Q，Jia J，Gao M，Han X 2010 Acta Phys.Sin.59 1576 (in Chinese)［張淑清、賈 健、高 敏、韓 敘 2010物理學報59 1576］

［5］Fraser A M，Swinney H L 1986 Phys.Rev.A 33 11340

［6］Xiao F H，Yan G R，Han Y H 2005 Acta Phys.Sin.54 550(in Chinese)［肖方紅、閻桂榮、韓宇航2005物理學報54 550］

［7］Yan H，Wei P，Xiao X C 2009 Chin.Phys.B 18 3287

［8］Zhang C T，Ma Q L，Peng H 2010 Acta Phys.Sin.7623(in Chinese)［張春濤、馬千里、彭 宏2010物理學報7623］

［9］Kantz H，Schreiber T 1997 Nonlinear Time Series Analysis (Cambridge:Cambridge University Press)p7

［10］Cao L Y，Mees A，Judd K 1998 Physica D 121 75

［11］Rombowts S A R B，Keunen R WM，Stom C J 1995 Phys.Lett. A 202 352

［12］Lu S，Wang H Y 2006 Acta Phys.Sin.55 572(in Chinese)［盧 山、王海燕2006物理學報55 572］

［13］Zhang Y，Guan W 2009 Acta Phys.Sin.58 756(in Chinese)［張 勇、關 偉2009物理學報58 756］

［14］Wang H Y，Sheng Z H，Zhang J 2003 J.Southeast Univ. (Natural Science Edition)33 115(in Chinese)［王海燕、盛昭瀚、張 進2003東南大學學報(自然科學版)33 115］

［15］Broomhead D S，King G P 1986 Physica D 20 217

［16］Fan C J 2009 J. UniversityofShanghaiforScienceand Technology 31 273(in Chinese)［樊重俊2009上海理工大學學報31 273］

［17］Han M，F(xiàn)an M M，Xi J H 2005 Lecture Notes in Computer Science 3497 618

［18］Han M，Wang Y 2009 Expert Systems with Applications 36 1280

［19］Tian Y C，F(xiàn)u X T，Lu Y Z 1993 Control and Decision 8 345(in Chinese)［田玉楚、符雪桐、呂勇哉1993控制與決策8 345］

［20］Elman J L 1990 Cognitive Science 14 179

［21］Ye M Y，Wang X D，Zhang H R 2005 Acta Phys.Sin.54 2568 (in Chinese)［葉美盈、汪曉東、張浩然 2005物理學報 54 2568］

［22］Ren R，Xu J，Zhu S H 2006 Acta Phys.Sin.55 555(in Chinese)［任 韌、徐 進、朱世華2006物理學報55 555］

［23］Mao J Q，Yao J，Ding H S 2009 Acta Phys.Sin.58 2220(in Chinese)［毛劍琴、姚 健、丁海山2009物理學報58 2220］

［24］Chiu S 1994 J.Intell.Fuzzy Syst.2 267

［25］Jang J S R 1993 IEEE Trans.Syst.Man 23 665

［26］Zhang J S，Li H C，Xiao X C 2005 Chin.Phys.14 49

［27］Jan S，Guatavo D 1997 Physics D 108 335

［28］Zhang J F，Hu S S 2008 Acta Phys.Sin.57 2708(in Chinese)［張軍峰、胡壽松2008物理學報57 2708］

［29］Gholpour A，Araabi B N，Lucas C 2006 Neural Proc.Lett.24 217

［30］Han M，Wei R 2008 14th International Conference of Neural Information Processing Kitakyushu，Japan，November 13—16，2007 p415

PACS:05.45.Pq，05.45.Tp

Multivariate chaotic time series phase space reconstruction based on extending dimension by conditional entropy*

Zhang Chun-TaoMa Qian-Li2)Peng Hong2)Jiang You-Yi1)
1)(College of Mathematic and Computer Science，Chongqing Three Gorges University，Wanzhou 404000，China)
2)(College of Computer and Engineering，South China University of Technology，Guangzhou 510006，China)

12 April 2010;revised manuscript

22 May 2010)

For multivariate chaotic time series，a method of conditional entropy extending dimension(CEED)in the reconstructed phase space is proposed.First，the delay time of any variable time series is selected by mutual information method，and then the embedding dimension of phase space is extended by the conditional entropy.This method can ensure the independence of reconstructed coordinates from low space to high space and eliminate the redundancy of phase space，because the largest condition entropy is choosen.The effective input vector for the prediction of multivariate time series is given.Simulations of the Lorenz system and Henon system show that the neural network predictions of multivariate time series are much better than the prediction of univariate and existing multivariate.Therefore，CEED is effective for multivariate chaotic systems.

multivariate chaotic time series，phase space reconstruction，conditional entropy，neural network prediction

*廣東省自然科學基金(批準號:9451064101003233)、華南理工大學中央高?；究蒲袠I(yè)務費專項資金(批準號:2009ZM0125，2009ZM0189，2009ZM0255)和重慶三峽學院重點項目(批準號:10ZD-16)資助的課題.

*Project supported by the Natural Science Foundation of Guangdong Province，China(Grant No.9451064101003233)，the Fundamental Research Funds for the Central Universities，South China University of Technology(Grant Nos.2009ZM0125，2009ZM0189，2009ZM0255)and the Key Program of the Chongqing Three Gorges University，China(Grant No.10ZD-16).