基于改進的G-P算法的相空間嵌入維數(shù)選擇

高俊杰,王豪

GAO Junjie,WANG Hao

上海交通大學 電子信息與電氣工程學院自動化系，系統(tǒng)控制與信息處理教育部重點實驗室,上海 200240

Department of Automation, Shanghai Jiao Tong University, and Key Laboratory of System Control and Information Processing, Ministry of Education of China, 200240, China

1 引言

混沌時間序列是由混沌動力系統(tǒng)產生的單變量或多變量的時間序列。相空間重構是混沌時間序列分析的基礎，重構的質量直接影響到后期分析的效果。1981年Takens提出了坐標延遲重構法[1]，認為只要選擇合適的延遲時間和嵌入維數(shù)，就能重構出一個與原系統(tǒng)具有相同拓撲性質的動力學系統(tǒng)。因此，延遲時間τ和嵌入維數(shù)m的選取很重要。

如果m選取過小，吸引子會折疊，甚至自相交，導致鄰域內包含不同軌道的相空間點，將使得重構后的吸引子形狀和原始吸引子不一致，影響后續(xù)研究的準確度；但如果m選取過大，吸引子的幾何結構完全打開，但卻將噪聲的影響放大，也導致計算量急劇增加。綜上，通常選取嵌入維數(shù)其中D是動力系統(tǒng)的關聯(lián)維數(shù)[2]。

確定關聯(lián)維數(shù)最常用的是 Grassberger和Procaccia提出的飽和關聯(lián)維數(shù)法（簡稱G-P算法）[3]。目前對該算法的研究主要集中在無標度區(qū)的合理選擇及自動識別上，比如：汪富泉等利用三折線段逼近法來確定無標度區(qū)[4]；費斌等通過遺傳算法選取標度區(qū)的起止點[5]；黨建武等根據(jù)置信度和相關性指標來選擇無標度區(qū)[6]。還有一部分研究是將G-P算法與其他算法結合來求取嵌入維數(shù)，比如呂威等將自相關法與G-P法構成循環(huán)迭代系統(tǒng)，通過不斷優(yōu)化得到參數(shù)[7]。但對于G-P算法其他方面的不足，則較少有人研究。所以本文針對原G-P算法的其他不足之處，提出了一種改進的算法，通過對鄰域半徑區(qū)間、步長和線性部分進行合理的選擇，更準確客觀地求取嵌入維數(shù)，而且通過在程序上的改進，能夠有效地減少運算量，縮短程序運行時間。

2 G-P算法

2.1 G-P算法原理

對混沌時間序列 x(n), n = 1,2,…,N，以延遲時間τ和嵌入維數(shù)m進行相空間重構，相空間中的相點為。給定一個鄰域半徑r，定義關聯(lián)積分如下：

式中：r為鄰域半徑，M=N-(m-1)τ為相空間中相點數(shù)目（N是時間序列長度），θ(·)為Heaviside單位函數(shù)（定義為θ(x)=0,ifx＜0;θ(x)=1,ifx≥0），為無窮范數(shù)。

關聯(lián)積分是一個累積分布函數(shù)，表征了相空間中任意兩個相點間距離小于鄰域半徑的概率，用來刻畫相點聚集程度。

當r→0時：如果奇異吸引子是一維結構，則C(m, r)正比于r；如果吸引子是二維結構，則C(m, r)正比于 r2；同理類推，當吸引子是D維結構時，則有C(m, r)正比于 rD，即：

即：

稱D為時間序列的關聯(lián)維數(shù)。

所以為了求取關聯(lián)維數(shù)，可以在坐標系中畫出一系列不同嵌入維數(shù)m時，ln C(m, r)～lnr的關系曲線，曲線中線性部分的斜率就是當系統(tǒng)嵌入維數(shù)為m時吸引子的關聯(lián)維數(shù)值 D(m)。一般來說，D (m)會隨m的增加而增大。當D(m)不再增加，而趨于穩(wěn)定時，即為飽和關聯(lián)維，對應的m即為奇異吸引子的最佳嵌入維數(shù)。這也是混沌序列與隨機序列的區(qū)別特征之一，隨機序列的關聯(lián)維數(shù)是不會趨于一個穩(wěn)定值的。

2.2 G-P算法的不足

通過深入分析和仿真試驗，發(fā)現(xiàn)G-P算法有以下幾個不足：

(1) 鄰域半徑r的選取主觀隨意性太強。r的選取直接影響關聯(lián)積分C(m, r)，若r取得太小，則鄰域內沒有包含任何相點，C(m, r)趨于0；若r取得太大，超過了相空間中相點的最大距離，將包含所有相點，C(m, r)趨于1。這兩種情況都將給關聯(lián)維數(shù)的判斷造成困難。

目前對于r值區(qū)間范圍的選取仍沒有統(tǒng)一的方法，大多是通過試湊法、效果觀察法來決定，個人主觀性太強，缺乏客觀性。而ln C(m, r)～ln r關系曲線在不同的r值區(qū)間段上表現(xiàn)出的斜率變化情況是差異較大的，這使得判斷線性部分斜率飽和與否更加困難，甚至誤判關聯(lián)維數(shù)值。如圖1所示，是對同一個樣本數(shù)據(jù)，采用不同的r值區(qū)間范圍，通過G-P算法得出的關系曲線。可以看出曲線的變化之大，從而得出結論：r值區(qū)間范圍選取的規(guī)范性對結果有重大影響。

圖1 同一樣本數(shù)據(jù)，不同r值區(qū)間范圍的關系曲線

(2) 鄰域半徑r變化步長的選取。步長太小，將極大地增加運算量，而且噪聲影響也會被放大，影響正確的結果；步長太大，樣本點太少，精度難以保證。有些文獻為了達到計算點在對數(shù)坐標中等間距的效果，采用按2的次冪來取點，但這樣做樣本點太少，精確度難以保證。

(3) 無標度區(qū)間的判定不客觀，斜率的飽和與否靠主觀判斷。系統(tǒng)關聯(lián)維數(shù)值就是ln C(m, r)～lnr關系曲線上線性部分（稱為無標度區(qū)間）的斜率。而對無標度區(qū)間的判定目前主要是靠主觀觀察，選定主觀所認為的線性區(qū)間，然后主觀判斷該區(qū)間線段的斜率是否達到飽和。但從圖 1可以看出，對斜率飽和與否進行主觀判斷容易有出入。

(4) 求取關聯(lián)積分時運算量巨大，如果不加以優(yōu)化，算法無法滿足實用性要求。按照原算法，如果m取 Km個不同值，r取 Kr個不同值：就要進行 Km次相空間重構以得到相點；每次重構后針對不同的r值要進行 Kr次個相點的距離計算；對于長度為m維的兩個相點，要得到其距離需要做m組相點元素間的減法運算及值比較；并將這些點距離與r比較得關聯(lián)積分。但本文分析發(fā)現(xiàn)這樣的算法原理和多層循環(huán)嵌套的計算結構存在大量重復計算。

3 改進的G-P算法

3.1 改進的G-P算法原理

針對G-P算法的上述不足，本文提出以下改進方法：

(1) 對鄰域半徑r值區(qū)間的選取，本文以相點間距離的最小值為r值下限，以相點間距離的最大值為r值上限，即。如此的區(qū)間選取方法具有兩個好處：第一，算法本身可根據(jù)序列的數(shù)值分布情況，自適應調整區(qū)間范圍，既排除了人為主觀性因素，也克服了區(qū)間固定導致的無法適應不同序列的難題；第二，可以綜合考慮時間序列的全局性質，具有較好的穩(wěn)定性。

(2) 步長采用均勻遞增的方法，保證 C(m, r)的平緩變化及足夠的樣本點，使無標度區(qū)間在雙對數(shù)坐標系中容易判斷，并滿足精度要求。

(3) 針對無標度區(qū)間的確定，本文采用基于BDS統(tǒng)計限定范圍的快速自動判定法。根據(jù) BDS統(tǒng)計結論[8]：如果時間序列是獨立同分布，則當N→∞時，關聯(lián)積分有。此時重構相空間中的點幾乎是均勻分布，重構吸引子軌道在相空間完全展開。 并且r=k·σ/2,k=1,2,3,4（σ是時間序列的均方差）；N≥3000;m=2,3,4,5時，符合這個結論的合理估計。因此，本文選定[ln(σ/2),ln(2σ) ]作為無標度區(qū)間的界限是有數(shù)理統(tǒng)計依據(jù)的。下面的仿真試驗也證明了這種選取的合理性。

再對選定的該范圍內的曲線用最小二乘法分別擬合，得到各條線段的斜率（稱為一次斜率，也即D(m)值），據(jù)此畫出D(m)～m的關系曲線；最后通過對D(m)～m曲線的斜率（稱為二次斜率）計算，當二次斜率小于所設定的閾值、充分接近零時，即說明 D(m)趨于飽和，此時即為系統(tǒng)關聯(lián)維，所對應的m即為最佳嵌入維數(shù)，如圖3、圖4所示。至此可以實現(xiàn)維數(shù)的自動計算。

此處改進利用 BDS統(tǒng)計對無標度區(qū)間的范圍進行了有效限定，可以大大減少后續(xù)擬合求一次和二次斜率的運算量，從而快速得到結果。而且進一步保證了無標度區(qū)間的正確鎖定。

(4) 本文主要從三個方面對算法進行優(yōu)化以達到提高速率的目的。

a) 從算法原理上去除重復計算。

首先，相空間重構是為了形成相點以計算點間距離，完全可以在計算距離時按重構規(guī)則直接從原始數(shù)據(jù)中讀取得到相點的各維元素值，沒必要進行Km次重構，既減少運算量，又解決 Km次重構帶來的內存占用問題，這對于數(shù)據(jù)量大的序列效果尤其顯著。

所以前m次減法運算是重復的，只要用以下遞推公式（5）就能進一步減少運算量：

b) 從程序結構上去除重復計算

在相同的m值條件下，相點間的距離并不會隨r變化，因此由于使用循環(huán)嵌套的計算結構而導致的重復計算rK次相點距離是沒必要的。新算法對程序結構進行調整：先用上一步的方法只計算一次相點間距離；再將點距離與r比較；而且，根據(jù)傳遞原理“”，又能進一步減少大量的重復比較運算。

c) 變“數(shù)據(jù)計算”為“數(shù)據(jù)讀取”。

對于計算機，從內存中讀取數(shù)據(jù)遠比計算數(shù)據(jù)要高效得多。為此本算法在內存中建立 Km個表，按如下規(guī)則存儲相點間的距離：第k個表格的第i行第j列元素表示嵌入維為mk時重構相空間中第i個相點和第j個相點的距離（存儲原理如圖2所示）。而且利用上一步的遞推方法，制作這 Km個表格的運算量并不大。計算關聯(lián)積分時，只要讀取對應的上三角矩陣的元素，并與r比較一次即可，去除了原算法中重復計算、重復比較的冗余運算。

圖2 相點距離存儲的原理示意圖

綜上，雖然由于多個參數(shù)的同時變化，無法精確量化本文的改進算法相對于原G-P算法在運算量上的改進程度，但從分析中是可以明顯看出新算法在減少運算量方面的效果的。

3.2 改進的G-P算法仿真試驗

根據(jù)改進的 G-P算法在 Matlab 7.6.0(R2008a)平臺上編成軟件，限于篇幅，僅以Lorenz、Rossler兩個典型的混沌系統(tǒng)為對象進行仿真試驗[9,10]。結果如表1和圖3、圖4所示。

結果分析：

1，改進算法計算出的關聯(lián)維非常接近理論值，說明是準確有效的；

2，適用于不同的混沌系統(tǒng)，說明改進算法的穩(wěn)定性和適用性都較好；

3，本算法的改進在程序運行速度上提高的程度與參數(shù) Km、 Kr、原始數(shù)據(jù)的長度有關，無法統(tǒng)一測定。但用上述例子進行多次測算，分別取和16，，統(tǒng)計平均運算時間僅為原算法的1/26和1/19。而且原時間序列的數(shù)據(jù)長度越長，運行時間的提高效果越明顯。說明本文在算法效率上的改進是有顯著效果的。

表1 改進的G-P算法仿真試驗

圖3 改進的G-P算法求解Lorenz系統(tǒng)x分量嵌入維數(shù)

圖4 改進的G-P算法求解Rossler系統(tǒng)x分量嵌入維數(shù)

4 結束語

本文針對G-P算法求取混沌時間序列相空間維數(shù)存在的不足，提出了改進。主要對鄰域半徑的區(qū)間選取、步長、無標度區(qū)間的判定三個方面進行了深入分析，并提出了改進方法，實現(xiàn)自動準確計算系統(tǒng)維數(shù)；并通過算法原理和程序結構的改良，去除重復、繁雜計算，提高程序效率。仿真試驗證明了本文提出的改進的G-P算法具有很好的準確性、穩(wěn)定性及實效性，可以為混沌時間序列的進一步分析提供良好的基礎。

[1]Takens F. Detecting strange attractors in turbulence[J].Lecture notes in Math,1981,898:361-381.

[2]Martin Casdagli. Chaos and Deterministic versus Stochastic Non-Linear Modelling[J]. Journal of the Royal Statistical Society:Series B (Methodological),1992,54(2):303-328.

[3]Grassberger P, Procaccia I. Characterization of strange attractors[J]. Physical Review Letters,1983,50(05):346-349.

[4]汪富泉，羅朝盛.G-P算法的改進及其應用[J].計算物理，1993,10(3):345-351.

[5]費斌，蔣莊德，王海容.基于遺傳算法求解分形無標度區(qū)的方法[J].西安交通大學學報，1998,32(7):72-75.

[6]黨建武，施怡，黃建國等.分形研究中無標度區(qū)的計算機識別[J].計算機工程與應用，2003,39(12):25-27.

[7]呂威，王和勇，姚正安等.改進嵌入維數(shù)和時間延遲計算的GP預測算法[J].計算機科學，2009,36(5):187-190.

[8]W A Broock, J A Scheinkman, W D Dechert, et al. A test for independence based on the correlation dimension[J].Econometric Reviews,1996,15(3):197-235.

[9]呂金虎，陸君安，陳士華.混沌時間序列分析及其應用[M].武漢：武漢大學出版社，2002:14.

[10]Argyris J, Andreadis I, Pavlos G, et al. The influence of noise on the correlation dimension of chaotic attractors[J].Chaos, Solitons & Fractals, 1998, 9(3): 343-361.