微熵率算法分析及實證研究

黃奕，謝維波

（1.華僑大學 計算機科學與技術學院，福建 廈門361021；2.華僑大學 廈門軟件園嵌入式技術開放實驗室，福建 廈門361008）

在時間序列的分析中，決定序列的可觀測因素很多，且相互作用的動力學方程往往是非線性的.20世紀80年代以來，由于Takens對Whitney早期在拓撲學方面工作的發展，使得深入分析時間序列的背景和動力學機制成為可能.在確定性的基礎上對序列動力學因素的分析，目前廣泛采用的是相空間重構法.微熵率法是Gautam 等提出的一個基于樣本時間序列及其替代數據的相空間重構方法［1］.熵率是指隨機源（1個隨機過程）隨時間變化的平均不確定性；1個隨機過程的熵率是該過程平均每產生1個隨機變量其不確定度大小的度量.微熵率法中的替代數據方法為iAAFT.iAAFT 是一種性能穩定的替代數據產生方法，能很好匹配原始數據的傅里葉幅度譜和概率密度分布，在數據的非線性檢驗中被廣泛采用.本文以Henon map混沌系統為數據源，以Henon map混沌特性的理論為依據，實證研究微熵率算法的各個環節.

1 Henon map混沌系統

1976年，Michel Henon給出了Henon map混沌系統［2－3］.此后，許多學者就Henon map混沌系統的全局結構進行數值與仿真研究，并將Henon map系統擴展到更高維［4］.

圖1 Henon混沌序列Fig.1 Henon chaotic sequence

1.1 Henon混沌序列

Henon混沌序列為

式（1）中：a，b為系統參數；d為系統延遲.當a＝1.4，b＝0.3時，系統具有混沌特性.經過足夠多次的迭代，Henon序列呈現“混沌”現象.

Henon混沌序列（d＝1），如圖1所示.圖1中：對應初值為（x0＝0.4，x1＝0.6）和（x0＝0.4＋10－8，x1＝0.6）的Henon混沌序列利用Matlab仿真的情況（實驗均采用Matlab完成）.大約在0＜k＜40的范圍內，xk序列基本一致；當k＞40之后，初始值微小擾動（10－8）的兩組xk序列呈現急劇的差異，體現了混沌序列“初值敏感性”的特征.

1.2 Henon map奇異吸引子

奇異吸引子是混沌運動的主要特征.令

由式（1），（2）可得

式（2），（3）構成了Henon map混沌系統［5］.點集（xk，yk）組成一條不封閉的曲線，即Henon map混沌系統的奇異吸引子相圖.

圖1兩種序列的奇異吸引子相圖，如圖2所示.由圖2可知：兩種序列幾乎“完全重合”.盡管圖1兩種序列存在急劇的差異，點集（xk，yk）卻呈現“相同的軌跡”.顯然，由圖1可知：對應點出現的“順序”是不同的，但是軌跡卻幾乎“完全重合”.奇異吸引子給出的“確定性軌跡”，體現了Henon map混沌系統的穩定性，奠定了混沌系統發展規律的研究基礎.

圖2 Henon map混沌系統奇異吸引子相圖Fig.2 Henon map chaotic system attractor

2 替代數據

為實現時間序列分析的統計學檢驗，時間序列需要足夠多樣本.樣本的獲取，可通過構造產生時間序列的系統，或者直接構造時間序列本身.后者產生時間序列樣本的方法，稱為替代數據法［6］.替代數據能夠盡可能精確地復制原始數據的性質（包括時間概率分布和自相關函數），但同時又是盡量隨機的.

2.1 AAFT生成算法

振幅調節傅里葉變換（AAFT）有以下4個步驟［7］.

步驟1原始數據｛x0n｝，排列序號｛rank0n｝，當x0n是｛x0n｝中第k小時，rank0n.

步驟2高斯白噪聲｛g0｝，按｛rank0n｝重排得｛g′n｝.

步驟3對｛g′n｝傅里葉變換和相位隨機化處理，得序列｛s′n｝.

步驟4求｛s′n｝的排列序號｛ranksn｝，原始數據｛x0n｝按｛ranksn｝重排得替代數據｛sn｝.

替代數據｛sn｝是原始數據｛x0n｝的重排，保證了替代數據與原始數據有相同的時間概率分布.因此，也有一樣的均值方差等一、二階統計量.替代數據｛sn｝的隨機性體現在高斯白噪聲｛gn｝，及其傅里葉變換的相位隨機化處理.

然而，原始數據功率譜密度（自相關函數）的性質被改變了，因為發生在步驟2和步驟4的兩個重排在嚴格意義上不是彼此的逆操作.這種差異導致替代數據的功率譜密度在原始數據的基礎上被白化.二者的差異程度取決于原始數據的時間概率分布與高斯分布的相似程度，簡而言之，AAFT 算法適用于類高斯分布的時間序列.

2.2 iAAFT生成算法

為了解決AAFT 替代數據的功率譜白化問題，1996年Schreiber提出了AAFT 迭代生成算法（iAAFT）.iAAFT 是一種性能穩定的替代數據產生方法，能很好匹配原始數據的傅里葉頻譜和概率密度分布，在數據的統計學檢驗中被廣泛采用.具體有以下4個步驟［8］.

步驟1原始數據｛x0n｝，排列序號｛rank0n｝，傅里葉變換的幅度值|Xk|，計算AAFT 替代數據｛sn｝.

步驟2記｛sn｝傅里葉變換Sk＝|Sk|exp（jφ（k）），保持相位不變，幅度值用|Xk|代替，得S′k＝|Xk|exp（jφ（k））.

步驟3對S′k進行傅里葉反變換得｛s′n｝，再按｛rank0n｝重排得｛s″n｝.

步驟4重復步驟2，3，直至所得數據｛s″n｝和原始數據有相近的功率譜密度.

確保｛s″n｝和原始數據的統計性質（時間概率分布、功率譜密度）趨于一致，該算法有兩個基本假設.

1）步驟2對傅里葉變換幅度值的矯正，造成時間概率分布的扭曲都比上一次迭代的小.

2）步驟3的重排，造成功率譜密度的扭曲都比上一次迭代的小.

Schreiber等［7］證明：對于非線性自相關過程，替代數據的功率譜將會逐漸趨近原數據的功率譜.

3 微熵率法確定時間序列最佳嵌入參數

3.1 微熵率法

時間序列可觀測性的決定因素很多，其相互作用的動力學方程往往是非線性的，甚至是混沌的.同時，計算的復雜性、有限的測量精度，以及可能存在的本質上的非確定性等多方面困難，嚴重制約著人們對時間序列內在機制的理解.20世紀80年代，Takens［9］對Whitney早期在拓撲學方面工作的發展，為深入分析時間序列的背景和動力學機制奠定基礎.在確定性的基礎上，對序列動力學因素的分析，目前廣泛采用的是延遲坐標狀態空間重構法（相空間重構法）［10］.一般來說，非線性系統的相空間可能維數很高，甚至無窮，在大多數情況下維數并不知道.

對于給定的時間序列，相空間重構法表明存在一個最優的嵌入維數m和時延τ.如果τ太小，為了使m·τ覆蓋（大于）“捕捉信號的動力學性質所需的最小時間距”，m將變得相當大；相反，如果τ大于最佳值，模型的性質將變得太離散，導致捕捉不到信號的動力學性質.Gautama等提出基于樣本時間序列及其替代數據的微熵率方法，用于確定相空間的最佳嵌入維數m和時延τ.

3.2 微熵率法詳細步驟

微熵率法有以下4點詳細步驟［1］.

步驟1原始數據｛x（k）∶k＝1，2，…，N｝，計算其Ns組替代數據，記為

步驟2嵌入維數m和時延τ的相空間.原始數據的相空間為

相空間的狀態向量記為

替代數據的相空間（共Ns個），即

相空間的狀態向量記為

其中：X（k）和Xs，i（k）又稱為延遲矢量，延遲矢量的個數為M＝N－（m－1）τ.

步驟3基于原始數據及其替代數據，確定原始數據｛x（k）∶k＝1，2，…，N｝的熵率為

其中：I（m，τ）＝H（x，m，τ）／（H（xs，i，m，τ））i.由原始數據構成的相空間，其熵為

其中：歐拉常數CE＝0.577 2；ρj是第j個延遲矢量與其最近鄰點的歐氏距離，即

相應有

由Ns組替代數據構成的Ns個相空間，其平均熵為可見，熵率體現為原始數據相空間的熵與Ns個替代數據相空間平均熵之比.

步驟4計算Rent（m，τ）的最小值，相應的m和τ即最佳嵌入維數mopt和時延τopt.

4 數值實驗

4.1 實驗數據的選取

Henon序列混沌特性的理論結果：當a＝1.4，b＝0.3時，式（1）具有混沌特性，其參數d為時間延遲，對應最佳嵌入時延τopt.

不當的初值選擇會造成式（1）的發散，分別選取兩組初值以驗證Henon序列混沌特征的穩定性.經過足夠多次的迭代，Henon序列才會進入“穩定的”混沌狀態［11］.取Henon序列10 000個數據后的500個作為數據源，才得出“穩定的”結果.

4.2 實驗結果

Henon序列的時延d＝4，初值?。▁0＝0.2，x1＝0.6x2＝0.4，x3＝0.3，x4＝0.1，x5＝0.5，x6＝0.8，x7＝0.7）和（x0＝0.2＋10－8，x1＝0.6，x2＝0.4，x3＝0.3，x4＝0.1，x5＝0.5，x6＝0.8，x7＝0.7），生成兩組上述的數據源（N＝500）應用微熵率法分別求最佳嵌入維數mopt和時延τopt.其中：Ns的選取，以（H（xs，i，m，τ））i趨于穩定為準，文中取Ns＝10.

圖3，表1，2對應d＝4的Henon序列.圖3和表1對應的初值?。▁0＝0.2，x1＝0.6，x2＝0.4，x3＝0.3，x4＝0.1，x5＝0.5，x6＝0.8，x7＝0.7）的結果.表2對應的初值?。▁0＝0.2＋10－8，x1＝0.6，x2＝0.4，x3＝0.3，x4＝0.1，x5＝0.5，x6＝0.8，x7＝0.7）的結果.從表1，2可以看出：m＝3，τ＝4時，熵率Rent（m，τ）均取得最小值，即最佳嵌入維數mopt＝3和時延τopt＝4（對應d＝4），實證了不同的初值下Henon序列混沌特征的穩定性.

圖3 Rent（m，τ）三維示意圖Fig.3 3D－diagram of Rent（m，τ）

表1 第一組初值對應的RentTab.1 Rentof the first set initial value

表2 第二組初值對應的RentTab.2 Rentof the second set initial value

對應Henon序列的不同時延d，表3給出基于微熵率的辨識結果（初值?。?，1］的隨機數），包括mopt，τopt，Rent（mopt，τopt）.由表3可知：mopt的值恒為3，τopt的值始終與d保持一致；有效地驗證了Henon序列混沌特征的穩定性，這種穩定性奠定了混沌系統的研究基礎.

表3 不同時延的微熵率法辨識Tab.3 Differential entropy of different time delay

［1］GAUTAMA T，MANDIC D P，VANHULLE M M.A differential entropy based method for determing the optimal embedding parameters of signal［C］∥Processing of the Int Conf on a Coustics，Speech and Signal Processing.Hong Kong：［s.n.］，2003：29－32.

［2］王光義，鄭艷，劉敬彪.一個超混沌Lorenz吸引子及其電路實現［J］.物理學報，2007，56（6）：3113－3120.

［3］GALLAS J A C.Structure of the parameter space of the Hénon map［J］.Physical Review Letters，1993，70（18）：2714.

［4］HENON M.A two dimensional mapping with stranger attractor［J］.Commun Math Phys，1976，50（1）：69－77.

［5］BENEDICKS M，CARLESON L.The dynamics of the Henon map［J］.Annals of Mathematics.1991，163（3）：749－841.

［6］THEILER J，EUBANK S，LONGTIN A，et al.Testing for nonlinearity in time series：The method of surrogate data［J］.Physica D：Nonlinear Phenomena，1992，58（1）：77－94.

［7］SCHREIBER T，SCHMITZ A.Improved surrogate data for nonlinearity tests［J］.Physical Review Letters，1996，77（4）：635.

［8］KEYLOCK C J.Constrained surrogate time series with preservation of the mean and variance structure［J］.Physical Review E，2006，73（3）：036707.

［9］王海燕，盛昭瀚.混沌時間序列相空間重構參數的選取方法［J］.東南大學學報：自然科學版，2000，30（5）：113－117.

［10］TAKENS F.Detecting strange atractors in turbulence［J］.Lecture Notes in Mathematics，1982，898（12）：198－366.

［11］GRASSBERGER P，KANTZ H，MOENIG U.On the symbolic dynamics of the Hénon map［J］.Journal of Physics A：Mathematical and General，1989，22（24）：5217.