羅文相，萬 麗，2，賴斯敏

(1.廣州大學數學與信息科學學院，廣東 廣州 510006; 2.廣州大學數學與交叉科學廣東普通高校重點實驗室，廣東 廣州 510006)

隨著非線性科學的發(fā)展，時間序列突變檢測越來越得到人們的關注。傳統(tǒng)的序列突變檢測方法，檢測結果嚴重依賴于時間序列的長度[1]。因此，許多學者將熵的理論引入到時間序列突變檢測領域。熵是系統(tǒng)復雜度的一種度量，可用于時間序列動力學結構突變的檢測。在信息論建立后，關于熵的理論和應用得到了迅速發(fā)展[2]。Pincus等[3-4]基于邊緣概率分布統(tǒng)計衡量時間序列復雜度提出的近似熵,被廣泛地應用于生理序列分析、礦化識別及徑流突變分析等領域[5-7]; 針對近似熵存在自身匹配特性的缺點，Richman等[8]對近似熵作了改進，提出了樣本熵。進一步，Bandt等[9]提出了排列熵(permutation entropy，PE)，與近似熵、樣本熵及分形維數、Lyapunov指數等復雜度參數相比，具有算法簡單，抗噪能力強等特點，能有效地放大序列的微小變化，廣泛應用于信號突變檢測，并在機械故障診斷、醫(yī)學和氣候等領域取得良好的應用效果[10-15]。

目前，對排列熵的應用大多只是單純地計算出時間序列的單個熵值，根據熵值的大小來判斷序列的復雜度。本文將排列熵與滑動窗口和滑動移除數據技術相結合，比較了滑動排列熵(moving permutation entropy，M-PE) 和滑動移除排列熵 (moving cut data-permutation entropy，MC-PE ) 方法在非線性時間序列動力學結構突變檢測的性能，為排列熵的實際應用提供更好的參考方法[16-17]。

1 時間序列突變檢測方法——排列熵

1.1 PE方法

PE算法的基本原理如下:

設一長度為N的時間序列X= {x(n)，n=1，2,…，N}，對時間序列X進行相空間重構，得到N- (m- 1)τ行，m列的矩陣

Xk=[x(i)，x(i+τ),…，x(i+ (m-1)τ)] ，i=1，2,…，N-(m-1)τ。

(1)

(1)式中：m是嵌入維，τ是延遲時間。Xk中的每一行為一個重構分量，共有N-(m-1)τ個。將第i個重構分量 (x(i)，x(i+τ),…，x(i+ (m- 1)τ)) 按照升序重新排列，得到

x(i+(j1-1)τ) ≤x(i+ (j2-1)τ) ≤…≤x(i+ (jm- 1)τ)。

(2)

(2)式中：j1、j2、…、jm表示元素所在重構分量中列的索引，若元素相等，則按照元素索引的大小進行排列。因此，矩陣Xk中的每個重構分量重新按照升序排列，可得到一個N-(m-1)τ行m列的符號序列矩陣A。記A(i)=(j1,j2,…，jm)，為(1，2,…，m)的某一種排列，顯然m個元素最多有m!種不同排列。以時間序列{5，3，8，6，4}為例，當m=3，τ= 1時，得到的第一個重構分量為(5，3，8)，由于xt+1

對所有A(i)的排列情況進行統(tǒng)計，記錄A(i)中某種排列A(r)的個數為N(r)，其中r≤m!，0

(3)

則排列熵(E)的定義為:

E(m)=-∑P(r)log(P(r))。

(4)

(4)式中：log表示為以2為底的對數。當P(r)=1/(m!) 時，E(m)取得最大值log(m!)。在實際應用中，通常會對E(m)進行歸一化處理，即:

0 ≤E(m)=E(m)/log(m!) ≤ 1。

(5)

E(m)值的大小與時間序列X的隨機程度表現(xiàn)一致，E(m)值越大，說明時間序列X的隨機程度越高，反之則X越規(guī)則。

1.2 M-PE方法

M-PE方法是一種將排列熵與滑動窗口技術相結合的一種突變檢測方法，步驟是對一時間序列，以長度為W的滑動窗口選取時間序列的數據構成子序列，采用PE方法計算該子序列的排列熵值，保持滑動窗口W不變，以步長為S逐步滑動選取新的子序列，同樣采用PE方法計算新子序列的排列熵值，直到時間序列結束，得到一個隨步長逐步變化的排列熵值序列。根據排列熵值序列的波動變化，判斷序列的突變點或突變區(qū)間。

1.3 MC-PE方法

MC-PE方法是一種將排列熵與滑動移除數據技術相結合的一種突變檢測方法,具體的步驟如下:

1) 選擇滑動移除數據的窗口長度M；

2) 從時間序列{x(n)}的第t(t=1，2,…，N-M+ 1，N為時間序列{x(n)}的長度)個數據開始連續(xù)移除M個數據，再將剩余的N-M個數據連在一起，得到一個長度為N-M的新子序列；

3) 利用排列熵方法計算新子序列的排列熵值；

4) 保持數據移除的窗口長度不變，取滑動移除步長為M，逐步移動窗口，重復2～3步操作，可得到一個長度為int(N/M)(int表示取整)的排列熵序列；

5) 基于相同動力學性質的數據復雜度差異不大，而不同動力學性質的數據復雜度不相同這一特點，結合步驟4得到的排列熵序列確定突變點或突變區(qū)間。

對于同一個具有相關性時間序列，任意移除動力學性質相同的數據，對其E(m)值的影響幾乎可以忽略不計。因此，可以通過步驟1～5觀察E(m)值的變化來檢測時間序列的動力學結構突變。

2 M-PE和MC-PE方法的比較

2.1 非線性理想時間序列的構造

構造由Logistic映射產生的非線性理想時間序列IS，序列總長度為2 000，其方程如下:

xn+1=λxn(1-xn)，xn∈[0,1] 。

(6)

圖1 非線性理想時間序列ISFig. 1 Nonlinear ideal time series IS

(6)式中：初始值為x0= 0.7，λ∈[0,4] ，表示蟲口模型的生長率，IS的前1 000個數據的參數為λ= 3.6，后1 000個數據的參數變?yōu)棣? 3.7，兩者均處于混沌狀態(tài)，其演化曲線由圖1給出。從IS的構造可知，在n=1 001時，系統(tǒng)的演化由較低的生長率3.6變?yōu)檩^高的生長率3.7。排列熵算法的參數均取m=3,τ=2,得出序列IS前1 000個數據的排列熵值為0.575 2；后1 000個數據的排列熵值為0.661 2。IS后1 000個數據的排列熵值比前1 000個數據大，由排列熵的物理意義可知,在n=1 001后系統(tǒng)的復雜度增大，可預測性變小。

2.2 M-PE的非線性時間序列狀態(tài)識別

M-PE方法中滑動窗口W需要選取適當的大小，如果滑動窗口過小，則因所含數據量小，不能很好地反應原始序列的特征；如果滑動窗口過大，則得到的滑動排列熵值個數少，忽略了原始序列的部分細致特點。經過多組數值試驗，發(fā)現(xiàn)選取大小為原始序列長度的5%～10%的滑動窗口能取得比較好的效果。選取滑動窗口W=150,對理想時間序列IS進行滑動排列熵的計算，滑動步長分別為S=10，20，30和50的檢測結果由圖2給出。

圖2 非線性時間序列IS的M-PE檢測結果Fig. 2 M-PE detection results of nonlinear time series IS

由圖2可見，對不同的滑動步長，M-PE方法得到的E(m)值變化趨勢基本一致，在n=1 000前后，E(m)值呈現(xiàn)出兩種不同的穩(wěn)定狀態(tài)。前半部分的E(m)值處于較低水平，而后半部分的E(m)值處于較高水平，由排列熵物理意義，可知前半部分數據的復雜度比后半部分的復雜度低，與理想時間序列IS的構造基本一致，說明M-PE方法具有識別序列動力學結構突變的能力。從圖2可以清楚地看到，在前后不同水平的E(m)值中間存在一個長度與滑動窗口相近的上升區(qū)域，這是因為隨著滑動窗口的移動，滑動窗口中所包含的數據由完全來自IS前1 000的數據，逐漸變?yōu)榍? 000數據和后1 000數據的混合，最后只含有來自IS后1 000的數據。區(qū)間(1 000，1 150)是兩種數據混合的階段，在這個階段，由于后1 000數據的個數逐漸增加，滑動窗口數據的復雜度也隨之增加，直到滑動窗口數據的復雜度接近于后1 000數據的復雜度，因此形成了與滑動窗口大小相近的上升突變區(qū)間。從上升區(qū)間的分析中，可以判斷突變點的位置大致在上升區(qū)間的開始處，即約在n=1 001處，說明M-PE方法能檢測時間序列的動力學結構突變。

2.3 MC-PE的非線性時間序列狀態(tài)識別

作為對比，MC-PE方法中滑動移除步長(即滑動移除窗口)同樣分別選擇M= 10，20，30和50，其檢測結果如圖3所示。

圖3 非線性時間序列IS的MC-PE檢測結果Fig.3 MC-PE detection results of nonlinear time series IS

4種滑動移除步長得到的MC-PE檢測結果幾乎是一致的，以n=1 000為界，移除相應數據后得到的排列熵值分為兩種不同的演化階段，在n>1 000的E(m)值明顯小于n≤1 000時的情形，而前后兩種階段的E(m)值各自處于相同的狀態(tài)。根據排列熵的物理意義，排列熵值越小表示時間序列的復雜度越小，而移除數據后得到較小的E(m)值意味著所移除的數據的復雜度較大，即移除復雜度較大的數據使得剩余數據的復雜度變小，說明序列在n>1 000處的復雜度高于前半部分的復雜度，這與理想時間序列IS的構造完全一致。圖3(a)～3(d)均顯示出E(m)值在n=1 000前后存在明顯差異，而且隨著滑動移除窗口M的增大，移除前后不同動力學結構的數據后所得的E(m)值之間的差異逐漸增大，因此可以判斷IS的突變點大約在n=1 001處，說明MC-PE方法能有效地識別時間序列的動力學結構突變。

2.4 M-PE和MC-PE結果對比分析

M-PE和MC-PE方法均能識別時間序列的動力學結構突變，但M-PE方法的結果存在一個大小與滑動窗口相近的上升突變區(qū)間，當滑動窗口增大時，上升區(qū)間也隨之增大，因此比較難精確地判定上升突變區(qū)間的開始位置，即M-PE方法得出的突變點位置誤差比較大。另外，滑動窗口W的選取與時間序列的長度有關,選取滑動窗口W=30和50，滑動步長S=1對IS進行M-PE檢測，結果顯示，在n=1 000附近的突變區(qū)域分界不明確，而且前半部分和后半部分各自的波動比較大(圖略)，不能有效地判斷突變點的位置。因此，不能為了提高M-PE突變檢測結果的精確度而刻意減少滑動窗口W的大小。

相對于M-PE方法依賴于滑動窗口的大小，導致突變點位置誤差較大這一缺點，MC-PE方法的突變檢測結果就更為精確。從圖3的結果可以看出，前半部分和后半部分的E(m)值之間的突變區(qū)域與滑動移除窗口M(即滑動移除步長)相近，換言之，MC-PE方法得出的突變點位置可以精確到滑動移除窗口M范圍內，這大大提高了突變點的精確度。另外，MC-PE方法幾乎不依賴于滑動移除窗口的大小，對于選取較小的滑動移除窗口，MC-PE方法依然能識別出突變點的位置。圖4給出了滑動移除窗口M=1和3時，MC-PE方法對時間序列IS的突變檢測結果。

圖4 IS的MC-PE檢測結果Fig.4 MC-PE detection results of IS

如圖4(a)所示,可以十分清晰地看到E(m)值以n=1 000為界，前后分別處于兩種不同的狀態(tài)，由于M=1，因此可以精確地判斷IS的突變點為n=1 001。在實際中，一般不會選取M=1，而是選取相對大一點的滑動移除窗口，因為每次只移除一個數據，計算成本太高，而且只移除一個數據，對實測數據的影響不是十分明顯。圖4(b)的結果表明，滑動移除窗口M=3比M=1更能反映序列IS后1 000數據的復雜度高于前1 000數據的復雜度，因為圖4(b)后半部分的E(m)值絕大部分都處于較低水平。綜上所述，M-PE方法對時間序列突變具有一定的識別能力，但該方法依賴于滑動窗口的選取，阻礙了其在實際中的應用。 MC-PE方法的檢測結果幾乎不依賴于滑動移除窗口的大小，能夠更為精準地得出時間序列的突變位置，明顯優(yōu)于M-PE方法。

3 結語

基于衡量時間序列復雜度的參數——排列熵，結合滑動窗口和滑動移除數據技術，比較了M-PE和MC-PE方法在非線性時間序列動力學結構突變檢測的性能。分析表明：M-PE方法對時間序列突變具有一定的識別能力，但檢測結果存在一個大小與滑動窗口相近的上升突變區(qū)間，得出的突變點位置誤差比較大，且該方法依賴于滑動窗口的選取，阻礙了在實際中的應用。與M-PE相比，MC-PE方法的檢測結果幾乎不依賴于滑動移除窗口的大小，能夠更為精準地得出時間序列的突變位置，即便是選擇非常小的滑動移除窗口 (M=1)，MC-PE方法也能有效地識別出序列IS的突變點，明顯優(yōu)于M-PE方法，具有良好的應用前景。