微熵率算法分析及實(shí)證研究

黃奕，謝維波

（1.華僑大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，福建 廈門361021；2.華僑大學(xué) 廈門軟件園嵌入式技術(shù)開(kāi)放實(shí)驗(yàn)室，福建 廈門361008）

在時(shí)間序列的分析中，決定序列的可觀測(cè)因素很多，且相互作用的動(dòng)力學(xué)方程往往是非線性的.20世紀(jì)80年代以來(lái)，由于Takens對(duì)Whitney早期在拓?fù)鋵W(xué)方面工作的發(fā)展，使得深入分析時(shí)間序列的背景和動(dòng)力學(xué)機(jī)制成為可能.在確定性的基礎(chǔ)上對(duì)序列動(dòng)力學(xué)因素的分析，目前廣泛采用的是相空間重構(gòu)法.微熵率法是Gautam 等提出的一個(gè)基于樣本時(shí)間序列及其替代數(shù)據(jù)的相空間重構(gòu)方法［1］.熵率是指隨機(jī)源（1個(gè)隨機(jī)過(guò)程）隨時(shí)間變化的平均不確定性；1個(gè)隨機(jī)過(guò)程的熵率是該過(guò)程平均每產(chǎn)生1個(gè)隨機(jī)變量其不確定度大小的度量.微熵率法中的替代數(shù)據(jù)方法為iAAFT.iAAFT 是一種性能穩(wěn)定的替代數(shù)據(jù)產(chǎn)生方法，能很好匹配原始數(shù)據(jù)的傅里葉幅度譜和概率密度分布，在數(shù)據(jù)的非線性檢驗(yàn)中被廣泛采用.本文以Henon map混沌系統(tǒng)為數(shù)據(jù)源，以Henon map混沌特性的理論為依據(jù)，實(shí)證研究微熵率算法的各個(gè)環(huán)節(jié).

1 Henon map混沌系統(tǒng)

1976年，Michel Henon給出了Henon map混沌系統(tǒng)［2－3］.此后，許多學(xué)者就Henon map混沌系統(tǒng)的全局結(jié)構(gòu)進(jìn)行數(shù)值與仿真研究，并將Henon map系統(tǒng)擴(kuò)展到更高維［4］.

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

1.1 Henon混沌序列

Henon混沌序列為

式（1）中：a，b為系統(tǒng)參數(shù)；d為系統(tǒng)延遲.當(dāng)a＝1.4，b＝0.3時(shí)，系統(tǒng)具有混沌特性.經(jīng)過(guò)足夠多次的迭代，Henon序列呈現(xiàn)“混沌”現(xiàn)象.

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

1.2 Henon map奇異吸引子

奇異吸引子是混沌運(yùn)動(dòng)的主要特征.令

由式（1），（2）可得

式（2），（3）構(gòu)成了Henon map混沌系統(tǒng)［5］.點(diǎn)集（xk，yk）組成一條不封閉的曲線，即Henon map混沌系統(tǒng)的奇異吸引子相圖.

圖1兩種序列的奇異吸引子相圖，如圖2所示.由圖2可知：兩種序列幾乎“完全重合”.盡管圖1兩種序列存在急劇的差異，點(diǎn)集（xk，yk）卻呈現(xiàn)“相同的軌跡”.顯然，由圖1可知：對(duì)應(yīng)點(diǎn)出現(xiàn)的“順序”是不同的，但是軌跡卻幾乎“完全重合”.奇異吸引子給出的“確定性軌跡”，體現(xiàn)了Henon map混沌系統(tǒng)的穩(wěn)定性，奠定了混沌系統(tǒng)發(fā)展規(guī)律的研究基礎(chǔ).

圖2 Henon map混沌系統(tǒng)奇異吸引子相圖Fig.2 Henon map chaotic system attractor

2 替代數(shù)據(jù)

為實(shí)現(xiàn)時(shí)間序列分析的統(tǒng)計(jì)學(xué)檢驗(yàn)，時(shí)間序列需要足夠多樣本.樣本的獲取，可通過(guò)構(gòu)造產(chǎn)生時(shí)間序列的系統(tǒng)，或者直接構(gòu)造時(shí)間序列本身.后者產(chǎn)生時(shí)間序列樣本的方法，稱為替代數(shù)據(jù)法［6］.替代數(shù)據(jù)能夠盡可能精確地復(fù)制原始數(shù)據(jù)的性質(zhì)（包括時(shí)間概率分布和自相關(guān)函數(shù)），但同時(shí)又是盡量隨機(jī)的.

2.1 AAFT生成算法

振幅調(diào)節(jié)傅里葉變換（AAFT）有以下4個(gè)步驟［7］.

步驟1原始數(shù)據(jù)｛x0n｝，排列序號(hào)｛rank0n｝，當(dāng)x0n是｛x0n｝中第k小時(shí)，rank0n.

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

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

步驟4求｛s′n｝的排列序號(hào)｛ranksn｝，原始數(shù)據(jù)｛x0n｝按｛ranksn｝重排得替代數(shù)據(jù)｛sn｝.

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

然而，原始數(shù)據(jù)功率譜密度（自相關(guān)函數(shù)）的性質(zhì)被改變了，因?yàn)榘l(fā)生在步驟2和步驟4的兩個(gè)重排在嚴(yán)格意義上不是彼此的逆操作.這種差異導(dǎo)致替代數(shù)據(jù)的功率譜密度在原始數(shù)據(jù)的基礎(chǔ)上被白化.二者的差異程度取決于原始數(shù)據(jù)的時(shí)間概率分布與高斯分布的相似程度，簡(jiǎn)而言之，AAFT 算法適用于類高斯分布的時(shí)間序列.

2.2 iAAFT生成算法

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

步驟1原始數(shù)據(jù)｛x0n｝，排列序號(hào)｛rank0n｝，傅里葉變換的幅度值|Xk|，計(jì)算AAFT 替代數(shù)據(jù)｛sn｝.

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

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

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

確保｛s″n｝和原始數(shù)據(jù)的統(tǒng)計(jì)性質(zhì)（時(shí)間概率分布、功率譜密度）趨于一致，該算法有兩個(gè)基本假設(shè).

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

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

Schreiber等［7］證明：對(duì)于非線性自相關(guān)過(guò)程，替代數(shù)據(jù)的功率譜將會(huì)逐漸趨近原數(shù)據(jù)的功率譜.

3 微熵率法確定時(shí)間序列最佳嵌入?yún)?shù)

3.1 微熵率法

時(shí)間序列可觀測(cè)性的決定因素很多，其相互作用的動(dòng)力學(xué)方程往往是非線性的，甚至是混沌的.同時(shí)，計(jì)算的復(fù)雜性、有限的測(cè)量精度，以及可能存在的本質(zhì)上的非確定性等多方面困難，嚴(yán)重制約著人們對(duì)時(shí)間序列內(nèi)在機(jī)制的理解.20世紀(jì)80年代，Takens［9］對(duì)Whitney早期在拓?fù)鋵W(xué)方面工作的發(fā)展，為深入分析時(shí)間序列的背景和動(dòng)力學(xué)機(jī)制奠定基礎(chǔ).在確定性的基礎(chǔ)上，對(duì)序列動(dòng)力學(xué)因素的分析，目前廣泛采用的是延遲坐標(biāo)狀態(tài)空間重構(gòu)法（相空間重構(gòu)法）［10］.一般來(lái)說(shuō)，非線性系統(tǒng)的相空間可能維數(shù)很高，甚至無(wú)窮，在大多數(shù)情況下維數(shù)并不知道.

對(duì)于給定的時(shí)間序列，相空間重構(gòu)法表明存在一個(gè)最優(yōu)的嵌入維數(shù)m和時(shí)延τ.如果τ太小，為了使m·τ覆蓋（大于）“捕捉信號(hào)的動(dòng)力學(xué)性質(zhì)所需的最小時(shí)間距”，m將變得相當(dāng)大；相反，如果τ大于最佳值，模型的性質(zhì)將變得太離散，導(dǎo)致捕捉不到信號(hào)的動(dòng)力學(xué)性質(zhì).Gautama等提出基于樣本時(shí)間序列及其替代數(shù)據(jù)的微熵率方法，用于確定相空間的最佳嵌入維數(shù)m和時(shí)延τ.

3.2 微熵率法詳細(xì)步驟

微熵率法有以下4點(diǎn)詳細(xì)步驟［1］.

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

步驟2嵌入維數(shù)m和時(shí)延τ的相空間.原始數(shù)據(jù)的相空間為

相空間的狀態(tài)向量記為

替代數(shù)據(jù)的相空間（共Ns個(gè)），即

相空間的狀態(tài)向量記為

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

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

其中：I（m，τ）＝H（x，m，τ）／（H（xs，i，m，τ））i.由原始數(shù)據(jù)構(gòu)成的相空間，其熵為

其中：歐拉常數(shù)CE＝0.577 2；ρj是第j個(gè)延遲矢量與其最近鄰點(diǎn)的歐氏距離，即

相應(yīng)有

由Ns組替代數(shù)據(jù)構(gòu)成的Ns個(gè)相空間，其平均熵為可見(jiàn)，熵率體現(xiàn)為原始數(shù)據(jù)相空間的熵與Ns個(gè)替代數(shù)據(jù)相空間平均熵之比.

步驟4計(jì)算Rent（m，τ）的最小值，相應(yīng)的m和τ即最佳嵌入維數(shù)mopt和時(shí)延τopt.

4 數(shù)值實(shí)驗(yàn)

4.1 實(shí)驗(yàn)數(shù)據(jù)的選取

Henon序列混沌特性的理論結(jié)果：當(dāng)a＝1.4，b＝0.3時(shí)，式（1）具有混沌特性，其參數(shù)d為時(shí)間延遲，對(duì)應(yīng)最佳嵌入時(shí)延τopt.

不當(dāng)?shù)某踔颠x擇會(huì)造成式（1）的發(fā)散，分別選取兩組初值以驗(yàn)證Henon序列混沌特征的穩(wěn)定性.經(jīng)過(guò)足夠多次的迭代，Henon序列才會(huì)進(jìn)入“穩(wěn)定的”混沌狀態(tài)［11］.取Henon序列10 000個(gè)數(shù)據(jù)后的500個(gè)作為數(shù)據(jù)源，才得出“穩(wěn)定的”結(jié)果.

4.2 實(shí)驗(yàn)結(jié)果

Henon序列的時(shí)延d＝4，初值取（x0＝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），生成兩組上述的數(shù)據(jù)源（N＝500）應(yīng)用微熵率法分別求最佳嵌入維數(shù)mopt和時(shí)延τopt.其中：Ns的選取，以（H（xs，i，m，τ））i趨于穩(wěn)定為準(zhǔn)，文中取Ns＝10.

圖3，表1，2對(duì)應(yīng)d＝4的Henon序列.圖3和表1對(duì)應(yīng)的初值?。▁0＝0.2，x1＝0.6，x2＝0.4，x3＝0.3，x4＝0.1，x5＝0.5，x6＝0.8，x7＝0.7）的結(jié)果.表2對(duì)應(yīng)的初值?。▁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）的結(jié)果.從表1，2可以看出：m＝3，τ＝4時(shí)，熵率Rent（m，τ）均取得最小值，即最佳嵌入維數(shù)mopt＝3和時(shí)延τopt＝4（對(duì)應(yīng)d＝4），實(shí)證了不同的初值下Henon序列混沌特征的穩(wěn)定性.

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

表1 第一組初值對(duì)應(yīng)的RentTab.1 Rentof the first set initial value

表2 第二組初值對(duì)應(yīng)的RentTab.2 Rentof the second set initial value

對(duì)應(yīng)Henon序列的不同時(shí)延d，表3給出基于微熵率的辨識(shí)結(jié)果（初值?。?，1］的隨機(jī)數(shù)），包括mopt，τopt，Rent（mopt，τopt）.由表3可知：mopt的值恒為3，τopt的值始終與d保持一致；有效地驗(yàn)證了Henon序列混沌特征的穩(wěn)定性，這種穩(wěn)定性奠定了混沌系統(tǒng)的研究基礎(chǔ).

表3 不同時(shí)延的微熵率法辨識(shí)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］王光義，鄭艷，劉敬彪.一個(gè)超混沌Lorenz吸引子及其電路實(shí)現(xiàn)［J］.物理學(xué)報(bào)，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］王海燕，盛昭瀚.混沌時(shí)間序列相空間重構(gòu)參數(shù)的選取方法［J］.東南大學(xué)學(xué)報(bào)：自然科學(xué)版，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.