曾祥健, 萬 麗, 劉 慧

(廣州大學 數(shù)學與信息科學學院, 廣東 廣州 510006)

熵是度量系統(tǒng)復雜度的一種指標，一個系統(tǒng)狀態(tài)的改變情況能夠通過熵的變化來描述.自20世紀香農(nóng)(Shannon)從概率統(tǒng)計的角度出發(fā)，提出了信息熵的概念之后，關(guān)于熵的理論和應用得到了迅速發(fā)展[1]；Pincus[2]于1991年提出的近似熵克服了Kolmogorov熵在應用中的困難和不足；Richman等[3]于2000年提出的樣本熵(Sample Entropy, SampEn)解決了近似熵存在的自身數(shù)據(jù)段匹配的問題；而Bandt等[4]于2002年基于編碼思想提出了具有較高魯棒性和快速簡單特征的排列熵(Permutation Entropy, PE).Li等[5]于2015年提出了分布熵(Distribution Entropy, DistEn)，并取得了良好的應用.近似熵、樣本熵和排列熵等方法的計算結(jié)果嚴重依賴于嵌入維數(shù)和相似容限等預設參數(shù)，當嵌入維數(shù)或相似容限的取值超過某一特定范圍會造成計算結(jié)果的惡化[6-7]，而分布熵不僅對預設參數(shù)的依賴很小，還具有這些方法不具有的高穩(wěn)定性及對參數(shù)的低敏感性[8-9].本文將分布熵與滑動窗口技術(shù)相結(jié)合，提出滑動分布熵(Moving Distribution Entropy，M-DistEn)，將其運用在序列的混沌狀態(tài)識別，并與滑動樣本熵[10](Moving Sample Entropy，M-SampEn)和滑動排列熵[11](Moving Permutation Entropy，M-PE)相比較，探討滑動分布熵在混沌狀態(tài)識別的性能，為分布熵的實際應用提供參考.

1 理論與方法

1.1 DistEn算法

設時間序列為{u(i)，i=1,2,…,N}，分布熵的算法步驟如下[6]：

(1)相空間重構(gòu)：選擇嵌入維數(shù)m，構(gòu)造向量X(i)={u(i),u(i+1),…,u(i+m-1)}(1≤i≤N-m+1)；

(2)構(gòu)造向量距離矩陣：定義向量X(i)和X(j)間的切比雪夫距離di,j，di,j=max{|u(i+k)-u(j+k)|, 0≤i,j≤N-m-1, 0≤k≤m-1}，距離矩陣D={di,j}；

(3)建立經(jīng)驗概率密度函數(shù)：利用直方圖法生成距離矩陣D的經(jīng)驗概率密度函數(shù)ePDF，設M為直方圖的分組數(shù)，記落入第h組的概率為Ph(1≤h≤M)；

1.2 M-DistEn算法

M-DistEn方法是將滑動窗口技術(shù)與DistEn方法相結(jié)合的一種新方法，步驟是對原始時間序列從第一個數(shù)據(jù)開始選取長度為w的滑動窗口構(gòu)成子序列，然后固定窗口w，以步長s滑動選取新的子序列，采用DistEn方法計算每個子序列的DistEn值，直到滑動窗口完全覆蓋原序列，得到DistEn值序列，根據(jù)DistEn值序列的變化，判斷原始時間序列的狀態(tài)變化.

2 混沌狀態(tài)識別

2.1 混沌序列的構(gòu)造

Logistic映射是最常見的描述生物種群系統(tǒng)演化的數(shù)學模型，其方程如下：xn+1=xn(1-xn),xn∈[0,1].當3≤≤4時，系統(tǒng)由倍周期通向混沌，當參數(shù)達到極限值∞=3.569 945 673時，Logistic映射迭代時間序列終態(tài)為周期2∞，即系統(tǒng)此時進入混沌狀態(tài)[12]，見圖1(a).

(a)2.8<<4.0 (b)3.57<<3.70

Henon映射是被廣泛應用的二維混沌映射，其方程如下[13]：

當a∈[1.07, 1.4]，b=0.3時，系統(tǒng)處于混沌狀態(tài)中，見圖2(a).

(a)0.0

為了便于檢測M-DistEn方法的參數(shù)影響及其在混沌狀態(tài)識別中的應用，構(gòu)造由Logistic映射產(chǎn)生的理想序列IS1.取初始值x0=0.4，理想序列IS1的前2 000個數(shù)據(jù)的參數(shù)1=3.574，2 001至4 000數(shù)據(jù)的參數(shù)2=3.58，4 001至6 000數(shù)據(jù)的參數(shù)3=3.7，總長度為6 000.從IS1的構(gòu)造可知，系統(tǒng)的演化由生長率較低的1變成較高的2，再到更高的3，見圖1(b)，系統(tǒng)的復雜度也越來越高，其中2-1明顯小于3-2，即第一段與第二段處于相差程度較小的混沌狀態(tài)，第二段與第三段處于相差程度較大的混沌狀態(tài).

同樣構(gòu)造由Henon映射產(chǎn)生的理想序列IS2.取初始值x0=0.0，y0=0.0，b=0.3，理想序列IS2的前2 000個數(shù)據(jù)的參數(shù)a1=1.072，2 001至4 000數(shù)據(jù)的參數(shù)a2=1.078，4 001至6 000數(shù)據(jù)的參數(shù)a3=1.2，總長度為6 000.理想序列IS2中，a2-a1明顯小于a3-a2，見圖2(b)，與IS1類似，IS2中第一段與第二段處于相差程度較小的混沌狀態(tài)，而第二段與第三段則處于相差程度較大的混沌狀態(tài).

2.2 算法參數(shù)的選取

首先，分析不同步長對M-DistEn計算結(jié)果的影響.選取滑動窗口為w=400，步長分別為10、20、50，并計算理想序列IS1、IS2的DistEn值，其中IS1的計算結(jié)果見圖3(a)～圖3(b)，IS2的計算結(jié)果見圖3(d)～圖3(f).圖3(a)、圖3(d)的步長為10，可以更為細微地觀察到熵值的變化，有利于發(fā)現(xiàn)原序列中細微的狀態(tài)變化，但計算時間相對較長，圖3(b)、圖3(e)對應步長為20，圖3(c)、圖3(f)對應步長為50.從圖中可見，隨著步長的增大，DistEn值的曲線變化更為平滑，但總體趨勢是不變的，均能準確識別出IS1和IS2的三種混沌狀態(tài)變化位置分別在2 001和4 001處，這說明M-DistEn方法對序列混沌狀態(tài)的識別能力受步長影響較小，在步長較大的情況下也能準確識別序列的狀態(tài)變化，在實際應用中可以靈活的根據(jù)數(shù)據(jù)長度選取合適的步長，提高計算速度.

圖3 理想序列IS1及IS2的M-DistEn檢測結(jié)果(基于不同步長)

其次，分析不同滑動窗口大小對混沌狀態(tài)的識別能力，選取滑動步長s=20，滑動窗口大小分別為150、300、400，計算理想序列IS1、IS2的DistEn值，IS1的計算結(jié)果見圖4(a)～圖4(c)，IS2的計算結(jié)果見圖4(d)～圖4(f).可見，滑動窗口大小為150時，其熵值雖有波動，但也能準確識別出三種混沌狀態(tài)，隨著窗口大小的增加，其熵值也越來越穩(wěn)定，對混沌狀態(tài)的識別效果也更好，這說明雖然M-DistEn方法在滑動窗口較小的情況下仍然具有較為出色的識別能力，但為了保證DistEn值在不同混沌狀態(tài)下的準確性和穩(wěn)定性，滑動窗口不宜過少，理想滑動窗口大小應盡量在400以上.

圖4 理想序列IS1及IS2的M-DistEn檢測結(jié)果(基于不同窗口)

2.3 識別結(jié)果與討論

為了說明M-DistEn方法在不同混沌狀態(tài)識別上的優(yōu)勢，將M-DistEn、M-SampEn和M-PE這三種滑動熵方法對理想序列IS1和IS2的計算結(jié)果進行比較.

圖5 IS1的三種滑動熵方法檢測結(jié)果

圖5(a)、圖5(d)分別是M-DistEn方法在嵌入維數(shù)m=2和m=3時的計算結(jié)果，從圖中可見，M-DistEn在不同嵌入維數(shù)下均能準確識別出序列IS1中相差較小的混沌狀態(tài)和相差較大的混沌狀態(tài)，且狀態(tài)發(fā)生變化位置分別在n=2 001和n=4 001處，而隨著參數(shù)的增大，序列的復雜度也增大，其DistEn值也增大，計算結(jié)果與實際意義一致.

圖5(b)、圖5(e)分別是M-SampEn方法在嵌入維數(shù)m=2和m=3時的計算結(jié)果，從圖中可見，在不同嵌入維數(shù)下，對于相差較小的混沌狀態(tài)，雖然樣本熵值有增大的趨勢，但在n=2 001處無法準確識別出序列的混沌狀態(tài)發(fā)生了變化，而在n=4 001處熵值發(fā)生了明顯變化，后2 000個數(shù)據(jù)的熵值明顯增大，說明序列的狀態(tài)發(fā)生了明顯變化，序列復雜度增大，這說明M-SampEn方法可以準確識別IS1中相差較大的混沌狀態(tài)[14]，但無法準確識別相差較小的混沌狀態(tài)，且樣本熵值在不同混沌狀態(tài)下的波動較大，不如分布熵值穩(wěn)定.

圖5(c)、圖5(f)分別是M-PE方法在嵌入維數(shù)m=2和m=3時的計算結(jié)果，從圖5(c)可知，M-PE在嵌入維數(shù)m=2時，前4 000個數(shù)據(jù)的排列熵值保持不變，在n=4 001處排列熵值開始下降，錯誤地認為序列的復雜度降低，計算結(jié)果與實際意義不符，這是因為在嵌入維數(shù)m=2時重構(gòu)向量中包含的狀態(tài)太少，算法失去有效性，不能準確檢測序列的狀態(tài)變化[15-16]，因此，M-PE方法應避免在嵌入維數(shù)m=2時使用；從圖5(f)可知，在嵌入維數(shù)m=3時，前4 000個數(shù)據(jù)的排列熵值保持不變，錯誤地認為前4 000個數(shù)據(jù)來源于同一參數(shù)的混沌狀態(tài)，無法識別出前面兩種相差較小的混沌狀態(tài)，在n=4 001處排列熵值開始增大，說明序列的復雜度變高，與實際意義一致.因此，M-DistEn方法對差別較大的混沌狀態(tài)和差別較小的混沌狀態(tài)的識別能力均優(yōu)于M-SampEn方法和M-PE方法.

三種滑動熵方法對理想序列IS2進行計算，結(jié)果見圖6.從圖6(a)、圖6(d)可知M-DistEn方法不僅能準確識別出IS2中第二段與第三段相差較大的混沌狀態(tài)，第一段與第二段相差較小的混沌狀態(tài)也能準確識別；而圖6(b)、圖6(e)以及圖6(c)、圖6(f)的結(jié)果則表明，M-SampEn方法和M-PE方法只識別出了第二段與第三段相差較大的混沌狀態(tài)，未能準確識別出IS2中第一段與第二段相差較小的混沌狀態(tài)，即三種滑動熵方法對IS2的計算結(jié)果基本與IS1的計算結(jié)果一致.

3 結(jié) 論

將DistEn方法與滑動窗口結(jié)合提出M-DistEn，基于Logistic映射和Henon映射構(gòu)建了不同參數(shù)的混合序列(IS1，IS2)，以此檢驗M-DistEn方法對不同混沌狀態(tài)的識別效果，并與M-SampEn方法、M-PE方法的結(jié)果進行對比.研究發(fā)現(xiàn)：M-DistEn在不同步長下和不同滑動窗口下均保持良好的混沌識別能力，表明該方法受步長的影響很小，而且DistEn值在滑動窗口較大時則更為穩(wěn)定；對于理想序列IS1和IS2中差別較大的混沌狀態(tài)，三種滑動熵方法在不同嵌入維數(shù)下均能準確識別出兩種相差較大的混沌狀態(tài)，但DistEn值的穩(wěn)定性要優(yōu)于SampEn值與PE值；而對于IS1和IS2中差別較小的混沌狀態(tài)，只有M-DistEn方法能準確識別出混沌狀態(tài)的微小差異，M-SampEn方法與M-PE方法則無法識別，且M-DistEn值序列波動小于M-SampEn值序列與M-PE值序列，指示M-DistEn方法的混沌狀態(tài)識別最為有效，其識別結(jié)果的穩(wěn)定性、準確性都優(yōu)于M-SampEn方法和M-PE方法，為進一步應用提供參考.