饒國(guó)強(qiáng)， 馮輔周， 司愛(ài)威， 謝金良

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

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

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

1 排列熵算法

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

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

(1)

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

(2)

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

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

(3)

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

2 相空間重構(gòu)

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

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

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

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

3.1 τ和m的獨(dú)立確定方法

(1)互信息法

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

(4)

(5)

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

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

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

(6)

因此有：

(7)

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

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

K=1.87(N-1)0.4

(8)

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

(9)

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

(2)偽近鄰法

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

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

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

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

(10)

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

其最近鄰點(diǎn):

(11)

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

(12)

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

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

(14)

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

(15)

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

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

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

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

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

(16)

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

(17)

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

(18)

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

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

(19)

令N→∞有：

(20)

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

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

(21)

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

(22)

(23)

(24)

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

(25)

4 仿真與實(shí)驗(yàn)驗(yàn)證

4.1 仿真分析及對(duì)比

對(duì)于信號(hào)f(t)：

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

(26)

圖1 互信息隨時(shí)間延遲的變化曲線

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

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

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

(27)

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

圖和Scor(τ)的變化曲線

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

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

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

4.2 軸承異常檢測(cè)

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

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

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

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

(2)排列熵異常檢測(cè)

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

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

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

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

5 結(jié) 論

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

參 考 文 獻(xiàn)

[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]侯 威,封國(guó)林,董文杰.利用排列熵檢測(cè)近40年華北地區(qū)氣溫突變的研究[J].物理學(xué)報(bào),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)、非線性振動(dòng)信號(hào)檢測(cè)方法的研究[J]. 振動(dòng)與沖擊, 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].中國(guó)科技信息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.