兩種滑動(dòng)排列熵序列突變檢測(cè)的對(duì)比分析

羅文相，萬(wàn) 麗，2，賴斯敏

(1.廣州大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，廣東 廣州 510006; 2.廣州大學(xué)數(shù)學(xué)與交叉科學(xué)廣東普通高校重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510006)

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

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

1 時(shí)間序列突變檢測(cè)方法——排列熵

1.1 PE方法

PE算法的基本原理如下:

設(shè)一長(zhǎng)度為N的時(shí)間序列X= {x(n)，n=1，2,…，N}，對(duì)時(shí)間序列X進(jìn)行相空間重構(gòu)，得到N- (m- 1)τ行，m列的矩陣

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

(1)

(1)式中：m是嵌入維，τ是延遲時(shí)間。Xk中的每一行為一個(gè)重構(gòu)分量，共有N-(m-1)τ個(gè)。將第i個(gè)重構(gòu)分量 (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表示元素所在重構(gòu)分量中列的索引，若元素相等，則按照元素索引的大小進(jìn)行排列。因此，矩陣Xk中的每個(gè)重構(gòu)分量重新按照升序排列，可得到一個(gè)N-(m-1)τ行m列的符號(hào)序列矩陣A。記A(i)=(j1,j2,…，jm)，為(1，2,…，m)的某一種排列，顯然m個(gè)元素最多有m!種不同排列。以時(shí)間序列{5，3，8，6，4}為例，當(dāng)m=3，τ= 1時(shí)，得到的第一個(gè)重構(gòu)分量為(5，3，8)，由于xt+1

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

(3)

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

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

(4)

(4)式中：log表示為以2為底的對(duì)數(shù)。當(dāng)P(r)=1/(m!) 時(shí)，E(m)取得最大值log(m!)。在實(shí)際應(yīng)用中，通常會(huì)對(duì)E(m)進(jìn)行歸一化處理，即:

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

(5)

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

1.2 M-PE方法

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

1.3 MC-PE方法

MC-PE方法是一種將排列熵與滑動(dòng)移除數(shù)據(jù)技術(shù)相結(jié)合的一種突變檢測(cè)方法,具體的步驟如下:

1) 選擇滑動(dòng)移除數(shù)據(jù)的窗口長(zhǎng)度M；

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

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

4) 保持?jǐn)?shù)據(jù)移除的窗口長(zhǎng)度不變，取滑動(dòng)移除步長(zhǎng)為M，逐步移動(dòng)窗口，重復(fù)2～3步操作，可得到一個(gè)長(zhǎng)度為int(N/M)(int表示取整)的排列熵序列；

5) 基于相同動(dòng)力學(xué)性質(zhì)的數(shù)據(jù)復(fù)雜度差異不大，而不同動(dòng)力學(xué)性質(zhì)的數(shù)據(jù)復(fù)雜度不相同這一特點(diǎn)，結(jié)合步驟4得到的排列熵序列確定突變點(diǎn)或突變區(qū)間。

對(duì)于同一個(gè)具有相關(guān)性時(shí)間序列，任意移除動(dòng)力學(xué)性質(zhì)相同的數(shù)據(jù)，對(duì)其E(m)值的影響幾乎可以忽略不計(jì)。因此，可以通過(guò)步驟1～5觀察E(m)值的變化來(lái)檢測(cè)時(shí)間序列的動(dòng)力學(xué)結(jié)構(gòu)突變。

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

2.1 非線性理想時(shí)間序列的構(gòu)造

構(gòu)造由Logistic映射產(chǎn)生的非線性理想時(shí)間序列IS，序列總長(zhǎng)度為2 000，其方程如下:

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

(6)

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

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

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

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

圖2 非線性時(shí)間序列IS的M-PE檢測(cè)結(jié)果Fig. 2 M-PE detection results of nonlinear time series IS

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

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

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

圖3 非線性時(shí)間序列IS的MC-PE檢測(cè)結(jié)果Fig.3 MC-PE detection results of nonlinear time series IS

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

2.4 M-PE和MC-PE結(jié)果對(duì)比分析

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

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

圖4 IS的MC-PE檢測(cè)結(jié)果Fig.4 MC-PE detection results of IS

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

3 結(jié)語(yǔ)

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