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

高俊杰,王豪

GAO Junjie,WANG Hao

上海交通大學 電子信息與電氣工程學院自動化系，系統控制與信息處理教育部重點實驗室,上海 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 引言

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

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

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

2 G-P算法

2.1 G-P算法原理

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

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

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

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

即：

稱D為時間序列的關聯維數。

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

2.2 G-P算法的不足

通過深入分析和仿真試驗，發現G-P算法有以下幾個不足：

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

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

圖1 同一樣本數據，不同r值區間范圍的關系曲線

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

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

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

3 改進的G-P算法

3.1 改進的G-P算法原理

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

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

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

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

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

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

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

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

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

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

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

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

c) 變“數據計算”為“數據讀取”。

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

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

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

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

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

結果分析：

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

2，適用于不同的混沌系統，說明改進算法的穩定性和適用性都較好；

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

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

圖3 改進的G-P算法求解Lorenz系統x分量嵌入維數

圖4 改進的G-P算法求解Rossler系統x分量嵌入維數

4 結束語

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

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

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

[7]呂威，王和勇，姚正安等.改進嵌入維數和時間延遲計算的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.