饒國強， 馮輔周， 司愛威， 謝金良

(裝甲兵工程學(xué)院 機械工程系, 北京 100072)

排列熵是衡量一維時間序列復(fù)雜度的平均熵參數(shù)，其計算簡單、抗噪聲能力強，是一種新的動力學(xué)突變檢測方法，能夠較好地反映時間序列數(shù)據(jù)的微小變化[1]。目前，排列熵算法在醫(yī)學(xué)、環(huán)境等領(lǐng)域應(yīng)用廣泛[2-5]，主要是對復(fù)雜系統(tǒng)進行異常檢測，應(yīng)用效果比較突出。

排列熵算法受自身參數(shù)的影響比較大，參數(shù)的選擇仍然僅憑經(jīng)驗或直覺，從而使計算結(jié)果存在很大的不確定性。本文引入相空間重構(gòu)方法，選擇合適的排列熵算法參數(shù)確定方法，以期提高排列熵計算結(jié)果的準確性和算法效率。

1 排列熵算法

排列熵算法的基本原理如下：

設(shè)一時間序列{X(i),i=1,2,…,n}，對其進行相空間重構(gòu)，得到矩陣：

(1)

式中：m、τ分別為嵌入維數(shù)和延遲時間，K=n-(m-1)τ。矩陣中的每一行可看作一個重構(gòu)分量，共有K個重構(gòu)分量。然后對每個重構(gòu)分量的元素，根據(jù)數(shù)值大小按照升序重新排列，提取各個元素在原重構(gòu)分量中所在列的索引，可以得到不同的符號序列，m維相空間映射不同的符號序列總共有m!。若k種不同符號序列出現(xiàn)的概率為P1,P2,…,Pk，則按照Shannon熵的形式，排列熵可以定義為：

(2)

將HPE(m)進行歸一化處理，即：

0≤HPE=HPE/ln(m!)≤1

(3)

根據(jù)排列熵算法的基本原理，相空間重構(gòu)是模型的第一步，其延遲時間τ和嵌入維數(shù)m是排列熵算法的主要參數(shù)，對于排列熵的計算結(jié)果有較大的影響。

2 相空間重構(gòu)

相空間重構(gòu)方法是由Takens在1981年提出，其基本思想是：系統(tǒng)中任一分量的演化都是由與之相互作用著的其它分量所決定的，因此單一變量的時間序列應(yīng)該隱含整個系統(tǒng)的運動規(guī)律，考察一個分量，測量它在某些固定的時間延遲點上的數(shù)值，然后將其重構(gòu)為多維狀態(tài)空間中的高維向量，并使得重構(gòu)坐標(biāo)從低維到高維轉(zhuǎn)換時保持較強的獨立性，最終的重構(gòu)相空間具有較低的冗余度。相空間重構(gòu)有兩種方法：導(dǎo)數(shù)重構(gòu)法和延遲重構(gòu)法。延遲重構(gòu)法在混沌研究中得到了廣泛的應(yīng)用，而導(dǎo)數(shù)重構(gòu)法卻很少得到實際應(yīng)用。

在延遲重構(gòu)法中，關(guān)于τ和m的選取，主要有兩種觀點[6]:①兩者互不相關(guān)，即τ和m獨立確定；②兩者相關(guān)，即τ和m互相依賴。對于第一種觀點，通常應(yīng)用互信息法和偽鄰近點法相結(jié)合，先利用互信息法確定τ，然后在τ已知的情況下運用偽鄰近點法選取m；第二種觀點的代表性方法是關(guān)聯(lián)積分法(C-C算法)，該方法是通過構(gòu)造統(tǒng)計變量和延遲時間的關(guān)系來確定最佳延遲時間和嵌入窗寬，再根據(jù)嵌入窗寬確定嵌入維數(shù)。

3 排列熵算法參數(shù)優(yōu)化確定

對于設(shè)備不同運行狀態(tài)下的時間序列信號，采用相同參數(shù)對不同狀態(tài)下的時間序列進行相空間重構(gòu)，所得出的時間序列排列熵值將會不同。因此，對設(shè)備在某一運行狀態(tài)下的時間序列進行相空間重構(gòu)所確定的算法參數(shù)，可用于設(shè)備運行全過程時間序列的排列熵特征提取，并以排列熵值的變化來反映設(shè)備的不同運行狀態(tài)，從而達到異常檢測的目的?；谠撍悸?，提出構(gòu)建最佳相空間為目標(biāo)的算法參數(shù)確定方法。

3.1 τ和m的獨立確定方法

(1)互信息法

互信息法的基本原理是：對于系統(tǒng)的某一時間序列{x(t),t=1,2,…,N}，為了建立延遲時間τ的關(guān)系式，選取延遲序列x(t+τ)構(gòu)成新的點列y(t)，通過計算x(t)和y(t)的相關(guān)性來確定延遲時間τ。假設(shè)構(gòu)成點對的總數(shù)為n，n=N-τ。則對于兩個離散時間序列{x(1),x(2),…,x(n)}和{y(1),y(2),…,y(n)}對應(yīng)的系統(tǒng)X和Y，根據(jù)信息論的知識，從兩個系統(tǒng)測量中所獲得的平均信息量，即信息熵分別為：

(4)

(5)

其中，Px(xi)和Py(yi)分別為系統(tǒng)X和Y中事件xi和yi的概率。

在給定X的情況下，可采用式(6)來計算X和Y之間的互信息：

I(X,Y)=H(Y)-H(Y|X)

(6)

因此有：

(7)

其中,Pxy(xi,yj)為事件xi和yi的聯(lián)合分布概率，I(τ)是與時間延遲τ有關(guān)的函數(shù)。I(τ)的大小代表了已知系統(tǒng)X，即x(t)在已知情況下，系統(tǒng)Y也就是y(t)的確定性大小。I(τ)=0，表示y(t)完全不可預(yù)測，即x(t)和y(t)完全不相關(guān)；而I(τ)的極小值則表示x(t)和y(t)最大可能的不相關(guān)，重構(gòu)時使用I(τ)的第一個極小值對應(yīng)的時間τ作為最優(yōu)延遲時間。

在互信息的計算中，F(xiàn)raser等[7]提出的互信息算法的具體過程是在x(橫軸)，y(縱軸)邊緣分布等概率的基礎(chǔ)上劃分網(wǎng)格。分段數(shù)可以通過經(jīng)驗公式確定：

K=1.87(N-1)0.4

(8)

式中,N表示序列長度，K表示分段數(shù)。

(9)

式中，n=1,2,…,Zm，R(n)分別對應(yīng)某個網(wǎng)格(xi,yi)。

(2)偽近鄰法

偽近鄰法(False Nearest Neighbor，F(xiàn)NN)是確定嵌入維數(shù)最常用的方法，其主要思想是:如果m維重構(gòu)相空間中的近鄰點在m+1維空間中仍然是近鄰點,則該點稱作真近鄰點,否則稱為偽近鄰點。當(dāng)重構(gòu)維數(shù)m小于原吸引子的真實維數(shù)時會產(chǎn)生偽近鄰點；反之,偽近鄰點消失,吸引子被打開。在某個維數(shù)m0處,偽近鄰點百分比(False Nearest Neighbor Percent, FNNP)將驟然降至(或接近)0,且不再隨m增大而變化,這個驟變點處的m0值就是最小嵌入維數(shù)。

FNN法確定嵌入維數(shù)的步驟如下:

對于時間序列{x(1),x(2),x(3),…,x(N)},當(dāng)延遲時間為τ,嵌入維數(shù)為m,則相空間中的任一點:

Xm(i)={x(i),x(i+τ),x(i+2τ),…,x(i+(m-1)τ)}

(10)

式中,i=1,2,3,…,N-(m-1)τ。

其最近鄰點:

(11)

式中:j=1,2,3,…,N-(m-1)τ,且j≠i。

(12)

同理,當(dāng)嵌入維數(shù)為m+1時,可得Rm+1(i)?？梢?

現(xiàn)引入偽近鄰點的判據(jù)一:

(14)

閾值Rtol可在[10,50]之間選取。對無限長精確的數(shù)據(jù),用上述標(biāo)準可獲得較好的結(jié)果。對有限長具噪聲的數(shù)據(jù),補充判據(jù)二:

(15)

則對實測時間序列,m從2開始,計算虛假最近鄰點的比例,然后增加m,直到偽近鄰點的比例小于5%或偽近鄰點不再隨著m的增加而減少時, 可以認為完全打開,此時的m即為求得的最小嵌入維數(shù)[8]。

3.2 τ和m的聯(lián)合確定方法

Kim等[9]在實際計算中首次提出了關(guān)聯(lián)積分法(C-C算法)，該方法是一種能夠同時確定最佳時間延遲τ和最佳嵌入維數(shù)m的算法。關(guān)聯(lián)積分法主要是通過嵌入時間序列的關(guān)聯(lián)積分來構(gòu)造統(tǒng)計量，該統(tǒng)計量代表非線性時間序列的相關(guān)性。

對于時間序列{x(t),t=1,2,…,N}，根據(jù)重構(gòu)相空間中的相點：

y(i)={x(i),x(i+τ),…,x(i+(m-1)τ)}

(16)

則嵌入時間序列的關(guān)聯(lián)積分方程為:

(17)

式中：M=N-(m-1)τ，m是嵌入維數(shù),N為原時間序列點數(shù),τ是時間延遲,r是參考半徑,dij=‖y(ti)-y(tj)‖，θ(u)為Heaviside函數(shù)：

(18)

對于給定的任意閾值r，檢查相空間中點對之間的距離小于r的點對數(shù)量，它占所有點對總數(shù)的比例稱為關(guān)聯(lián)積分。對于時間序列{x(t),t=1,2,…,N}，將其分成τ個不相交的時間序列，τ為時間延遲，再由這τ個不相關(guān)的時間序列計算S(m,N,r,τ)。

對于一般的自然數(shù)τ，將原時間序列分解成τ個不相交的時間子序列:{x(1),x(1+τ),x(1+2τ),…}，{x(2),x(2+τ),x(2+2τ),…}，…，{x(τ),x(2τ),x(3τ),…}，然后定義每個子序列的S(m,N,r,τ)為：

(19)

令N→∞有：

(20)

如果相空間中的點是獨立分布且N→∞，S(m,r,τ)對于所有的r將等于0，真實的數(shù)據(jù)序列是有限長度，且序列元素間可能相關(guān)，因此一般S(m,r,τ)≠0，從而局部最佳時間可確定為S(m,r,τ)的首次過零點或表現(xiàn)出對不同r變化最小點，這意味著點的分布最接近均勻分布。對于所有r定義變量為:

ΔS(m,τ)=max{S(m,rj,τ)}-min{S(m,rj,τ)}

(21)

Brock等[10]模擬研究了各種類型的分布后指出，對有限時間序列長度，當(dāng)N≥500時，N,m,r的取值為一般地2≤m≤5，σ/2≤m≤2σ，(σ指時間序列的標(biāo)準差)，可以通過有限序列對漸近分布作很好的近似。具體計算時，取m=2,3,4,5，ri=iσ/2,i=1,2,3,4。計算下列三個統(tǒng)計量:

(22)

(23)

(24)

m=int[(τw/τ)+1]

(25)

4 仿真與實驗驗證

4.1 仿真分析及對比

對于信號f(t)：

f(t)=sin(2*π*t)+sin(2*π*10*t)

(26)

圖1 互信息隨時間延遲的變化曲線

信號在正常狀態(tài)下包含2個頻率成分，對于該狀態(tài)下的時間序列可采用前文所述的兩種不同的觀點對應(yīng)的方法來確定模型參數(shù)。首先采用1 000 Hz的采樣頻率對信號進行離散化，選擇重構(gòu)時間序列長度為1 024點，利用互信息法計算時間序列的時間延遲τ與互信息之間的關(guān)系，如圖1所示。按照互信息法確定延遲時間的原則，可知時間延遲τ=12。

圖2是給出了基于FNN法計算偽近鄰率隨嵌入維數(shù)的變化關(guān)系曲線。判據(jù)一的閾值Rtol選擇為15，判據(jù)二的閾值A(chǔ)tol選擇為2，綜合兩種判據(jù)的效果，可以計算其聯(lián)合判據(jù)曲線，如圖2中紅色圓圈連線所示。

為了驗證所得的算法參數(shù)對異常檢測的效果，在原始信號f(t)中增加一個頻率為20 Hz的正弦分量，以此模擬信號發(fā)生了異常，發(fā)生時刻為60 s，從而構(gòu)成帶有異常變化的信號f′(t)，如式(27)所示。

(27)

圖2 偽近鄰率隨嵌入維數(shù)的變化曲線

圖和Scor(τ)的變化曲線

圖4 模擬異常信號f′(t)

選擇采樣頻率1 000 Hz，對信號進行離散化，異常信號f′(t)的時域波形如圖4所示，將前文中基于f(t)信號所確定的模型參數(shù)延遲時間和嵌入維數(shù)作為排列熵的模型參數(shù)，相鄰兩個重構(gòu)時間序列之間沒有數(shù)據(jù)重疊，依次計算該異常信號的排列熵，結(jié)果如圖5所示。點的連線代表采用參數(shù)獨立確定方法對應(yīng)的排列熵值隨時間的變化曲線，圓圈連線代表采用參數(shù)聯(lián)合確定方法對應(yīng)的排列熵值隨時間的變化曲線。對比兩條曲線的變化特點可知，前者要優(yōu)于后者，一是前者曲線跳變更加明顯，即異常敏感性更高；二是前者曲線在非異常區(qū)域內(nèi)波動更小，即穩(wěn)態(tài)平穩(wěn)性更好。由此可見，排列熵模型中參數(shù)的獨立確定法能夠獲得更好的異常檢測效果。

圖5 排列熵異常檢測結(jié)果

4.2 軸承異常檢測

為了進一步檢驗方法的有效性，提出對軸承試驗臺全壽命振動信號數(shù)據(jù)[11]進行異常檢測，數(shù)據(jù)來自美國NSFI/UCR的智能維護系統(tǒng)中心(IMS)，本文采用了該試驗第2輪次的測試數(shù)據(jù)。

由于EMD方法具有自適應(yīng)分解特性，比較適用于非線性、非平穩(wěn)的信號分析，在此采用EMD方法對振動信號進行預(yù)處理。該實驗是軸承外圈發(fā)生故障，異常信息主要集中在高頻部分，因此對分解后的IMF1分量利用參數(shù)優(yōu)化方法確定排列熵模型參數(shù)，最后利用排列熵算法對軸承振動信號EMD分解后的IMF1分量進行異常檢測。

(1)模型參數(shù)的確定

選擇軸承在正常狀態(tài)下的振動信號數(shù)據(jù)，長度為1024，經(jīng)EMD分解后選擇IMF1分量，計算其延遲時間與互信息的關(guān)系。當(dāng)τ=1時互信息第一次取得極小值，故取延遲時間τ=1。然后基于FNN方法計算偽近鄰率和嵌入維數(shù)之間的關(guān)系，因此可確定m=4。

(2)排列熵異常檢測

對于實驗中每組數(shù)據(jù)，將其分為20段，每段長度為1024，各段之間沒有數(shù)據(jù)重疊，對每段進行EMD分解，計算IMF1分量的排列熵值，然后取其平均值作為一組數(shù)據(jù)的特征值，依次計算各組數(shù)據(jù)的排列熵值。兩種不同方法確定模型參數(shù)后得出的結(jié)果如圖6所示，橫坐標(biāo)是每組數(shù)據(jù)對應(yīng)的測試時間，縱坐標(biāo)表示排列熵結(jié)果，其中圖6中曲線(a)表示聯(lián)合確定算法參數(shù)后排列熵的計算結(jié)果變化曲線，圖6中曲線(b)表示獨立確定算法參數(shù)后排列熵變化曲線。圖7是任意選定算法參數(shù)后得出的排列熵計算結(jié)果，

顯然，經(jīng)過參數(shù)優(yōu)化后的模型輸出結(jié)果更好，主要表現(xiàn)在：一是異常表現(xiàn)效果更明顯；二是檢測異常的發(fā)生時間要更早。對于不同方法優(yōu)化參數(shù)的異常檢測結(jié)果可知，兩者預(yù)測趨勢基本相同。在異常起始點處，圖6中曲線(b)異常跳變更明顯；在異常區(qū)域內(nèi)，排列熵均經(jīng)歷了先減小而后逐漸穩(wěn)定的過程，但曲線(b)排列熵值變化要快，排列熵最小值與正常狀態(tài)之間的差值要大于0.04，而圖6中曲線(a)不到0.03，因此說明模型參數(shù)τ=1m=4對于異常檢測更敏感，檢測效果更明顯；圖中故障起始點均發(fā)生在同一時刻，排列熵值下降到最低點時刻也相同，但曲線(a)故障跳變要更明顯，并且此后，其排列熵值變化要快，說明參數(shù)τ=4m=4對于故障發(fā)生后檢測更敏感，檢測效果更明顯。因此，關(guān)于τ和m的選取，兩種確定參數(shù)的方法都具有可行性，但對于異常檢測來說，算法參數(shù)的獨立確定方法更有效。

圖6 排列熵計算結(jié)果變化曲線

圖7 τ=22，m=9時排列熵變化曲線

5 結(jié) 論

本文重點研究排列熵算法參數(shù)的優(yōu)化確定方法，根據(jù)模型的基本原理，提出以重構(gòu)時間序列最佳相空間的方法確定模型參數(shù)，克服了以往對排列熵參數(shù)人為主觀或憑經(jīng)驗選擇的不確定性。根據(jù)當(dāng)前相空間重構(gòu)的兩種不同觀點，介紹了不同觀點中比較典型算法的基本原理，利用仿真信號和滾動軸承全壽命周期數(shù)據(jù)直觀地說明了不同確定方法對異常檢測效果的影響，結(jié)果表明模型參數(shù)的獨立確定要優(yōu)于聯(lián)合確定方法，為排列熵算法在不同領(lǐng)域推廣應(yīng)用時，如何確定算法的最佳參數(shù)奠定了基礎(chǔ)。

參 考 文 獻

[1]Cao Y H, Tung W W, Gao J B,et al. Detecting dynamical changes in time series using the permutation entropy[J].Phys RevE,2004,70(4)：1-7.

[2]Nicolaou N, Georgiou J. Detection of epileptic electroencephalogram based on permutation entropy and support vector machines[J]. Expert Systems with Applications,2012,39(1)：202-209.

[3]Nair U, Krishna B M, Namboothiriet V N M. Permutation entropy based real-time chatter detection using audio signal in turning process[J]. International Journal of Advanced Manufacturing Technology, 201046:161-68.

[4]侯 威,封國林,董文杰.利用排列熵檢測近40年華北地區(qū)氣溫突變的研究[J].物理學(xué)報,2006,55(5):2663-2668.

HOU Wei, FENG Guo-ling, DONG Wen-jie. A technique for distinguishing dynamical species in the temperature time series of north china[J]. Acta Physica Sinica,2006,55(5):2663-2668.

[5]劉永斌, 龍 潛, 馮志華, 等. 一種非平穩(wěn)、非線性振動信號檢測方法的研究[J]. 振動與沖擊, 2007, 26(12)：131-134.

LIU Yong-bin, LONG Qian, FENG Zhi-hua, et al. Detection method for nonlinear and non-stationary signals[J]. Journal of Vibration and Shock, 2007, 26(12)：131-134.

[6]謝忠玉,張 立.相空間重構(gòu)參數(shù)選擇方法的研究[J].中國科技信息2009(16):42-43.

XIE Zhong-yu, ZHANG Li. Selection of embedding parameters in phase space reconstruction[J]. China Science and Technology Information,2009(16):42-43.

[7]Fraser A M, Swinney H L. Independent coordinates for strange attractors from mutual information [J]. Phys Rev A, 1986, 33(2):1134-1140．

[8]Kennel M B, Brown R, Abarbanel H D I. Determining embedding dimension for phase-space reconstruction using a geometrical construction[J]. Physical Review A, 1992, 45(6)：3403-3411.

[9]Kim H S,Eykholt R, Salas J D. Nonlinear dynamics, delaytimes and embedding windows [J]. Physica D, 1999, 127(1):48-60.

[10]Brock W A, Hsieh D D A,Lebaron B D al. Statistical theory and economic evidence[M]. M IT Press, Cambridge, MA, 1991.

[11]Feng F Z, Zhu D D, Jiang P C, et al. GA-SVR based bearing condition degradation prediction[J].Key Engineering Material2009413-414431.