基于小波包和樣本熵的水泵機組振動特征提取

宋 媛，彭利鴻，賴冠文，張嘉勛，肖志懷，宋麗波

(1.武漢大學動力與機械學院，武漢 430072；2.廣東省水科院，廣州 510000；3.黃河水利水電開發總公司，河南 濟源 454681)

隨著我國經濟水平快速穩定發展，南水北調工程的開工，一大批大、中型泵站將投入運行，水泵一旦發生故障將會造成巨大的經濟損失和嚴重的社會影響。泵站運行可靠性受到越來越多的重視，因此，有必要對水泵故障進行及時診斷,以提高水泵機組運行的可靠性[1]。

水泵機組振動信號中承載著豐富的機械設備運行狀態信息，80%的機組故障在振動信號中有所反映[2]。但泵站的工作環境惡劣，由于采集的振動信號受到影響，其中的有用信息和噪聲摻雜在一起，有用信息被抑制。如何從海量狀態信息中提取出有用的機組故障特征是機組故障診斷的一大難點。因此，從振動信號中提取故障特征在機組故障診斷領域獲得了廣泛應用，如：彭文季[3]等利用快速傅里葉變換把振動信號從時域轉換到頻域，然后提取其幅值特征進行機組故障診斷；趙道利[4]等對機組振動信號經小波分解后的單尺度信號進行傅里葉變換，提取信號的故障特征。但是，傅里葉變換實質上是信號的整體變換，其應用是建立在信號平穩性假設上的，而振動信號一般是非平穩信號，其波動性比較大，信號中的干擾比較多，其中的有用信息在信號中的成分一般是低頻部分或平穩部分，傅里葉變換無法實現信號在時間域上的特征，這使得利用傅里葉變換處理非平穩信號受到限制。

為了克服傅里葉變換的使用缺陷，近年來，一種基于時頻分析的方法——小波分析開始在機械設備的狀態監測與故障診斷領域有了較快的發展和進步。小波分析相對于傅里葉分析的優越之處在于它對復雜、非平穩信號的處理時，具有局部化分析能力，可以根據需要對分析對象進行任意細節的分析，是對傳統傅里葉變換的巨大改進。但小波分析同樣存在缺陷，利用小波變換處理信號時，只是進一步分解信號的低頻部分，而不再分解信號的高頻部分，它不能很好地分解和表示包含大量細節信息的中、高頻信號，如機械非平穩振動信號、地震信號等[5]。因此在小波分析的基礎上，后來的學者又提出了小波包分析。小波包分析是小波分析的進一步發展，克服了小波分析在分析信號高頻部分時的不足，可以更加精細分解信號的高頻部分，且自動匹配信號特征，因此小波包分析對非平穩信號局部的細節分析精度更好，能更加精細的展現信號的特征。樣本熵[6]對數據長度的依賴性很小，同時具有很好的抗干擾能力。因此將小波包變換和樣本熵相結合的分析方法引入到水泵機組的振動信號故障特征提取中，為信號故障診斷提供一個更為準確可靠的特征提取算法。

1 小波包分析

1.1 小波包理論

小波包變換是基于小波變換的進一步發展，能夠提供比小波變換更高的分辨率。設φ(t)正交尺度函數，ψ(t)是相應的正交小波，{Vk}是由φ(t)生成的多分辨分析空間序列，Wk-1是Vk中關于Vk-1正交補空間。這時L2(R)能寫成正交和分解，即：

L2(R)=…⊕W-1⊕W0⊕W1⊕…

(1)

而{2k/2ψ(2kt-l):l∈Z}是Wk的一個規范正交基。

為了對分解后的各個頻帶進一步進行細分，即把時頻進一步局部化，設：

ψ0(t)=φ(t)，ψ1(t)=ψ(t)

則φ和ψ的兩尺度關系可寫成：

(2)

式中：g(k)=(-1)kh(1-k)。由上式確定的函數ψn(t)，n=2l或2l+1，l=0，1，…，稱為關于尺度函數φ(t)的小波包。

1.2 小波包算法

設gnj(t)∈ψnj，則gnj可表示為：

(3)

小波包分解算法：由{dj+1,nl}求{dj,2nl}與{dj,2n+1l}。

(4)

小波包重構算法：由{dj+1,nl}以及{dj,2nl}求{dj,2n+1l}。

(5)

式中：dj,2nl和dj,2n+1l為下一層小波包分解結果；dj+1,nl為上層分解的結果；j為尺度系數；l為位置系數；n為頻率；k為變量；h和g分別為正交共軛低通和高通濾波器。

小波包對高頻小波系數做了進一步的分解，將隨分解尺度j增大而變寬的頻譜窗口進一步分割細化，克服了正交小波基隨著分解尺度j的增大，空間分辨率增大而頻率分辨率降低的缺陷，從而提高了頻率分辨率，更能適應信號的變化，對信號局部特征的刻畫更加清楚，時頻特性更好。

2 樣本熵

樣本熵(SampEn，Sample Entropy)是由Richman和Moornan[6]提出的一種新的時間序列復雜性測度方法，可以用SampEn(m，r，N)來表示。其中，N為長度，r為相似容限，維數為m及m+1。樣本熵旨在降低近似熵的誤差，與已知隨機部分有更加緊密的一致性。樣本熵是一種與近似熵類似但精度更高的方法[6]。樣本熵算法[7]如下：

Step1：設原始數據為u(1)，u(2)，… ，u(N)共N個點。

Step2：按序號連續順序組成一組m維矢量：從Xm(1),到Xm(N-m),其中:Xm(i)={u(i),u(i+1),…,u(i+m-1) },1≤i≤N-m+1。這些矢量代表著從第i個點開始連續的m個u的值。

Step3：定義矢量Xm(i)和Xm(j)間的距離d[Xm(i),Xm(j)]為兩者對應元素中差值最大的絕對值, 即:

d[Xm(i),Xm(j)]=max(|u(i+k)-u(j+k)|)

(6)

其中k=0～m-1;i,j=1～N-m,j≠i。

Step4：給定閾值r，對每個i≤N-m的值，統計d[Xm(i),Xm(j)]小于r的數目(稱為模板匹配數，記為Nmi)及此數目與距離總數N-m-1 的比值，記作：

(7)

對所有i對應的Bmi(r)求平均值，記為：

(8)

Step5：增加維數為m+1，按序號連續順序組成一組m+1維矢量：

從Xm+1(1)到Xm+1(N-m),其中:Xm+1(i)=[u(i),u(i+1),…,u(i+m)]。這些矢量代表著從第i個點開始連續的m+1個u的值。

Step6：定義矢量Xm+1(i)和Xm+1(j)間的距離d[Xm+1(i),Xm+1(j)]為兩者對應元素中差值最大的絕對值，即：

d[Xm+1(i),Xm+1(j)]=max(|u(i+k)-u(j+k)|)

(9)

其中k=0～m;i,j=1～N-m,j≠i。

Step7：給定閾值r，對每個i≤N-m的值，統計d[Xm+1(i),Xm+1(j)]小于r的數目(稱為模板匹配數，記為Nm+1i及此數目與距離總數N-m-1的比值，記作：

(10)

對所有i對應的Bm+1i(r)求平均值，記為：

(11)

理論上，此序列的樣本熵為：

(12)

樣本熵的計算過程很明顯地反映出m，r取值對于樣本熵值的巨大影響。根據Pincus[8]的研究結果，當m=1或2，r=0.1Std～0.25Std(Std是原始數據的標準差)時計算得到的樣本熵值具有較為合理的統計特性。

樣本熵具有以下優點：

(1)樣本熵不包含自身數據段的比較，因此它是條件概率的負平均自然對數的精確值，因此樣本熵的計算不依賴數據長度；

(2)樣本熵具有更好的一致性。即如一時間序列比另一時間序列有較高的值的話，那對于其他m和r值，也具有較高的值；

(3)樣本熵對于丟失數據不敏感。即使數據丟失多達1/3，對計算值影響依然很小。

綜合上述優點可知，樣本熵算法可以應用于隨機過程的研究。

3 基于小波包變換和樣本熵的特征提取

時頻域特征量可以同時反映分析對象的時域和頻域特性。和小波分析相比，小波包分解可以對信號的低頻部分進行處理，同時對信號的高頻部分進行處理，具有更好的時頻分辨率。當機械設備處于異常運行狀態或故障狀態時，對應信號的不規則性和不平穩性可以反映設備的狀態，對故障的發生和發展也會有所體現。首先使用小波包對信號進行分解，得到n個頻帶系數。對小波包頻帶系數進行重構，計算各頻帶的樣本熵值。以各頻帶信號的樣本熵值為元素，構造特征向量V=[SamEn1,SamEn2,SamEn3,…,SamEnn]。

3.1 試驗數據來源及處理

圖1為本文采用的信號采集系統。本采集系統主要由轉子振動試驗臺、轉子臺控制器、前置器以及計算機軟件系統構成。轉子振動試驗臺配備有一臺直流電機，配有DH5600轉軸控制器，通過改變轉子轉速、剛度、質量不平衡、軸的摩擦或沖擊條件以及聯軸器的型式來模擬旋轉機械振動情況。該設備可以模擬水泵機組的常見故障(如不平衡、不對中、碰磨等)。最終的振動信號通過傳感器(測量振動和轉速)傳遞到前置器，進行相關處理后，輸送到計算機軟件系統供存儲、分析使用。

圖1 轉子振動試驗臺

利用轉子試驗臺模擬了水泵機組運行過程中轉子的4種典型工況：①正常工況；②不平衡工況；③不對中工況；④碰磨工況。在轉子試驗臺上，通過對轉子的主軸或者轉盤進行設置來模擬不同工況。讓轉子正常運轉得到的就是正常工況下的振動信號，要模擬不平衡工況，就需要在轉盤上的螺紋孔里嵌入一個兩克重的質量塊；模擬不對中工況需要將法蘭處兩個軸的位置錯開放置；模擬碰磨工況則需要將碰磨螺栓旋入碰磨螺紋箱，采集到的信號即為所需的故障信號。信號轉速和采樣頻率都設置為相同大小，分別為1 200 r/min和2 048 Hz，共采集100組數據，每組數據都包含2 048個采樣點。通過小波改進閾值降噪[9]對采集信號進行3層分解，進行降噪處理。圖2為不同運行狀態下降噪后的信號時域圖。

圖2 不同運行狀態下降噪后的信號

從圖2可以看出，在經過了降噪處理后，4種運行狀態下的振動信號中的大量噪聲成分被去除，留下的是真實有用的適合后續故障診斷特征提取及模式識別的信號。并且在四種運行狀態下，去噪后可以直觀地看到不同運行狀態采集到的信號形狀各不相同，在分辨率允許的情況下，肉眼可以根據信號波形對水泵機組運行狀態進行粗略的分類。但是，為了能精確區分出水泵機組運行狀態，需要將這些信號波形轉化為可以量化的特征量。將這些特征量輸入到模式識別模型中可以進行水泵機組的智能診斷。

3.2 振動信號特征提取

本文使用小波包分析和樣本熵相結合的方法對降噪信號進行特征提取，借此區分水泵機組幾種不同的運行狀態。

首先使用DB8小波[10]對降噪后的信號進行小波包分解，分3層，其小波包分解樹結構如圖3所示。經小波包3層分解后，包含8個節點，每個節點對應于不同的頻帶。信號的采樣頻率為2 048 Hz，由香農采樣定理可知，信號的奈奎斯特頻率為1 024 Hz。可知，經小波包分解后對應節點對應的頻率范圍如表1所示。

圖3 小波包3層分解樹形結構

表1 小波包節點對應的頻率范圍

圖4～圖7是不同運行狀態下的水泵機組振動信號經小波包分解后得到的不同節點的小波包系數。從圖中可以看出，不同的運行狀態下，信號的能量集中在不同的頻帶范圍內，這為信號的特征提取提供了可能。計算各頻帶的樣本熵，作為特征向量，記作V=[SamEn1,SamEn2,SamEn3,SamEn4,SamEn5,SamEn6,SamEn7,SamEn8]。

從表2中的樣本熵值可以看出，不同運行狀態下，相同頻帶的樣本熵值有很大不同，可以很明顯地區分不同運行狀態，對水泵機組的不同運行狀態具有較高的區分度。

4 基于小波包變換和樣本熵的特征提取方法驗證

為了定性定量地驗證小波包變換與樣本熵結合的振動特征提取方法的實際效果。本文分別采集了水泵機組四種不同運行狀態下的振動信號若干組，使用學習向量量化(learning vector quantization，LVQ)神經網絡[11]對樣本進行分類。

LVQ算法的基本思想是：計算距離輸入向量最近的競爭層神經元，從而找到與之相連接的線性輸出層神經元，若輸入向量的類別與線性輸出層神經元所對應的類別一致，則對應的競爭層神經元權值沿著輸入向量的方向移動；反之，則對應的競爭層神經元權值沿著輸入向量的反方向移動。

LVQ神經網絡的結構如圖8所示，它由三層組成：即輸入層、隱含層和輸出層。該網絡在輸入層與隱含層間為完全連接，而在隱含層與輸出層間為部分連接，每個輸出神經元與隱含神經元的不同組相連接。

圖4 正常運行狀態下的振動信號小波包分解

圖5 不平衡運行狀態下的振動信號的小波包分解

圖6 不對中運行狀態下的振動信號的小波包分解

圖7 碰磨運行狀態下的振動信號的小波包分解

狀態SamEn1SamEn2SamEn3SamEn4SamEn5SamEn6SamEn7SamEn8正常不平衡不對中碰磨0.41060.35120.52080.45330.47560.37300.63240.58980.49000.38520.66110.60290.48610.38130.62820.57790.73420.74600.76090.73751.03531.01230.85850.83901.08120.99960.67520.81011.46081.42211.13731.2664

LVQ神經網絡有以下優點：

(1)結構簡單，只需要三層網絡就可以實現模式識別；

(2)不存在BP網絡有可能陷入局部最小問題；

(3)收斂速度比BP網絡更快；

(4)不需要將輸入向量進行歸一化、正交化，只需要直接計算輸入向量與競爭層之間的距離，從而實現模式識別,因此簡單易行，識別效率更高。

圖8 LVQ神經網絡結構圖

使用LVQ神經網絡驗證基于小波包分解與樣本熵的水泵機組振動信號特征提取方法的有效性，隨機選取樣本中60組用于訓練樣本，40組作為測試樣本。利用神經網絡工具箱[12]中的newlvq( )函數創建LVQ神經網絡，其中隱含層神經元個數為20。LVQ神經網絡測試結果如圖9所示。

圖9 測試集實際分類和預測分類圖

從圖9可以看出，40組測試樣本中僅有一組數據分類出現偏差，由此得出測試樣本的分類效果識別率達到97.5%，分類效果較好。由此說明基于小波包變換和樣本熵相結合的特征提取方法對水泵機組不同運行狀態具有較好的識別率。

5 結 語

本文在前人研究的基礎上，在信號特征提取方面提出了新的觀點和方法，該方法將小波包變換和樣本熵兩種工具結合起來，既考慮了信號的細節成分，同時樣本熵算法不依賴數據長度，因而使故障特征提取更精確。通過對轉子試驗臺采集的振動信號的處理與分類驗證表明，所提出的方法具有一定的優越性，可以作為信號處理方法的借鑒。

□

[1] 趙康德.基于小波變換和Lipschitz指數的水泵故障診斷研究[D]. 江蘇揚州：揚州大學，2010.

[2] 梁武科,張彥寧,羅興锜.水電機組故障診斷系統信號特征的提取[J].大電機技術,2003,(4):53-56.

[3] 彭文季,羅興锜,趙道利.基于頻譜法于徑向基函數網絡的水電機組振動故障診斷[J].中國電機工程學報，2006,(9):155-158.

[4] 趙道利,梁武科,羅興锜,等.水電機組振動信號的子帶能量特征提取方法研究[J].水力發電學報，2004,(6):116-119,115.

[5] 單立場.小波變換在吸毒者心率變異性信號分析中的應用[D].重慶:重慶大學,2007.

[6] Richman JS, Moorman JR. Physiological time-series analysis using approximate entropy and sample entropy[J]. American Journal of Physiology-Heart and Circulatory Physiology, 2002,278(6):2 039-2 049.

[7] R Alcaraz, J J Rieta. A review on sample entropy applications for the non-invasive analysis of atrial fibrillation electrocardiograms[J]. Biomedical Signal Processing and Control, 2010

[8] Pincus S M. Assessing serial irregularity and its implications for health [J]. Ann. N. Y. Acad. Sci, 2002,954:245-267.

[9] 安周鵬,肖志懷,孫召輝,等.改進小波閾值降噪算法在水電機組信號處理中的應用[J].中國農村水利水電，2014,(12):165-168,172.

[10] 費佩燕,劉曙光.幾種常見小波應用性能分析[C]∥中國電子學會第七屆學術年會論文集，2001.

[11] 劉文軒.基于LVQ神經網絡的配電網故障定位方法[J].電力系統保護與控制,2012,40(5):90-95.

[12] 史 峰.MATLAB神經網絡30個案例分析[M]. 北京:北京航空航天大學出版社,2010.