張廣濤，宋麗波，張曉靜，孫召輝，梁紅梅，肖志懷，程遠楚

(1.武漢大學動力與機械學院，武漢 430072；2.國網新源山東泰山抽水蓄能電站有限責任公司,山東 泰安 271000；3.新疆昌吉職業(yè)技術學院電氣工程分院，新疆 昌吉 831100)

0 引 言

據統(tǒng)計，在水電機組故障或事故中，約80%在振動信號中有所反映；同時，異常振動也是設備損壞的一項重要原因。因此，振動信號分析方法在水電機組故障診斷領域獲得了廣泛應用。但是，水電機組運行工況復雜，運行環(huán)境惡劣，現(xiàn)場振源豐富，實際采集的振動信號往往被淹沒在大量的背景噪聲中，給振動信號故障特征提取造成了不良影響。

為了有效地去除噪聲，同時考慮到振動信號的非平穩(wěn)特征，眾多學者將小波分析作為信號降噪的主要方法。Donoho[1]首次提出了硬閾值和軟閾值降噪方法。孫濤[2]等采用連續(xù)小波變換系數模極大值法，對水輪機低頻振動信號中故障點的奇異性程度進行估算并取得了很好的診斷效果。Xiangzhong Meng[3]等將小波多分辨率分析引入振動信號濾波中，運用Mallat算法對振動信號進行降噪處理，轉子平臺實驗表明，該方法適用于振動信號分析，降低了計算量，提高了計算效率。Zhao Ting[4]等應用一種基于小波變換的改進閾值降噪方法，通過MATLAB仿真實驗證明了當機械振動信號受到高斯白噪聲嚴重干擾時，此方法比其他降噪方法更有效，可以有效是地消除噪聲并保留信號的詳細信息。Isayed B[5]等將小波技術應用于振動監(jiān)測和故障檢測中，使用小波技術分析旋轉機械振動信號，通過小波包變換揭示信號的瞬態(tài)信息，利用細節(jié)信號中的能量平均值從振動信號中對故障進行有效檢測。Coifman R R[6]等針對小波降噪過程中在不連續(xù)點可能出現(xiàn)的Gibbs現(xiàn)象，提出了平移不變小波變換方法，采用不同平移量對原始信號進行的多次“平移-降噪-逆平移”處理，然后對所有降噪結果取平均值，消除小波平移依賴性。

盡管與小波分析相關的信號處理方法獲得了極大的發(fā)展，但是由于小波分析自身的缺陷 以及信號處理技術飛度發(fā)展，多小波分析方法應運而生，多小波同時滿足對稱性，短(緊)支撐性，二階消失矩和正交性這些特性，可以匹配信號中不同特征。近些年來，多小波在圖像處理、數據壓縮和故障診斷等領域都獲得了極大的發(fā)展，并確立了相對單小波的優(yōu)勢[7]。

Strela[8]等將多小波單變量閾值降噪方法應用于二維圖像信號的降噪，取得了良好的效果。謝榮生[9]等對信號和噪聲的多小波變換進行了研究，在多小波噪聲方差閾值的基礎上，提出了一種新的多小波信號濾波方法，仿真分析表明，該方法具有很好的濾波效果。劉志剛[10]等探討了基于多小波的信號去噪方法，提出了一種基于自適應閾值的多小波去噪方法，將其應用于電力系統(tǒng)輸電線路故障暫態(tài)信號的去噪，仿真結果表明，該方法可以根據實際信號自適應改變閾值的大小，在去噪效果上優(yōu)于傳統(tǒng)多小波去噪方法。這些方法都具有良好的效果，但是忽視了多小波系數間的相關性。文獻[11,12]引入了相鄰系數的概念，提出了多小波相鄰系數降噪方法，通過與其他方法的比較，顯示該方法可以取得更好的效果。

盡管多小波相鄰系數降噪方法已證明其優(yōu)越性，也在一些領域得到了應用[13-15]。但在水電機組振動信號降噪方面的應用還很少見。針對水電機組振動信號的特點，本文將圍繞多小波相鄰系數降噪方法完成對振動信號的降噪處理，針對硬閾值和軟閾值兩種情況和不同多小波類型和預處理方法，將降噪效果和當前廣泛使用的小波閾值降噪方法進行對比。通過降噪前后的定量分析，表明當使用合適的多小波與處理方法時，多小波相鄰系數降噪方法在水電機組振動信號降噪方面具有良好的效果。

1 多小波理論

1.1 多小波多分辨率分析

多小波理論是在小波理論的基礎上發(fā)展而來的，同小波一樣，多小波同樣滿足多分辨率分析。不同之處在于多小波的多尺度函數是以多個尺度函數{φ，r∈N}為基礎構成的向量函數。

多小波有多個尺度函數和小波函數，即v0由r個尺度函數的平移φ1(t-k)，φ2(t-k)，…，φr(t-k)生成。記φ1=[φ1(t),φ2(t),…,φr(t)]T，為多分辨率分析空間的多尺度函數，與其對應的正交多小波函數ψ(t)=[ψ1(t),ψ2(t),…,ψr(t)]T，滿足其平移和伸縮ψj,k={ψ1(2-jt-k),…,ψr(2-jt-k)}T(j,k∈Z)形成正交補子空間的正交基。類似于小波分析，多小波分析中的多尺度函數和多小波函數滿足兩尺度矩陣方程：

(1)

(2)

式中：hk、gk分別為φ(t)和ψ(t)對應的r×r的濾波器系數矩陣。

設f∈vJ，根據多小波的多分辨率分析理論可以得到：

(3)

將Mallat算法擴展至正交多小波多尺度系數和多小波系數的求解，可得：

多小波分解公式：

(5)

式中：vj,k為多小波分解后第j層第k個尺度系數；wj,k為多小波分解后第j層第k個小波系數。

多小波重構公式：

(6)

式中：vj,k=(v1j,k,v2j,k,…,vrj,k)T,wj,k=(w1j,k,w2j,k,…,wrj,k)T,h和g為r×r單位矩陣，分別代表多小波的低通和高通濾波器矩陣。

在實際應用中通常r=2，多小波分解與重構的Mallat算法流程圖如圖1，圖2所示。

圖1 多小波分解過程Mallat算法流程圖Fig.1 Multiwavelets decomposition flowchart using Mallat algorithm

圖2 多小波重構過程Mallat算法流程圖Fig.2 Multiwavelets reconstruction flowchart using Mallat algorithm

在眾多多小波類型中，GHM多小波[17]和CL3多小波[18]為最常用的多小波。

1.2 常用多小波

1.2.1GHM多小波

多小波領域最具代表性的成果是Geronimo，Hardin和Massopust[17]在1994年利用分形差值方法，成功構造出的正交、緊支撐、實對稱和二階逼近的GHM多小波。該多小波包含兩個尺度函數和兩個小波函數，多尺度函數的支撐區(qū)間分別為[0,1]和[0,2]，多小波函數的支撐區(qū)間為[0,2]。

GHM多小波的濾波器系數矩陣為：

1.2.2CL3多小波

CL多小波是由Chui C K和Lian J[18]于1995年提出的?？紤]不同長度的支撐區(qū)間，又可以分為CL3和CL4兩類。其中，CL3多小波具有正交性、對稱性、二階消失矩和緊支撐的特性，其尺度函數和小波函數的支撐區(qū)間為[0,2]，濾波器長度為3；而CL4多小波的尺度函數和小波函數的支撐區(qū)間為[0,3]，濾波器長度為4。為方便與GHM多小波對比，在此選取CL3多小波。

CL3多小波的濾波器系數矩陣為：

2 多小波相鄰系數法研究

多小波分析是在小波分析的基礎上發(fā)展而來的，是對小波分析的繼承和發(fā)展，。與小波閾值降噪類似，多小波降噪仍然是基于如下理論展開：假定故障信號經過小波分解，信號特征成分的能量會集中在少數大的小波系數上，而噪聲分解后的小波系數一般都非常小。因此，經過小波分解后，信號的小波變換系數要大于噪聲的小波變換系數，選擇合適的λ作為閾值，當wj,k小于該閾值時，認為此時wj,k主要由噪聲引起；而當wj,k大于該閾值時，認為該系數主要是信號引起的[19]。考慮到多小波分解系數中，當前系數與其相鄰系數之間存在一定的相關性，因此采用相鄰系數降噪法進行降噪。

2.1 多小波相鄰系數降噪算法

多小波相鄰系數降噪法的步驟如圖3所示，原始信號經過多小波分解后，按照設計的閾值處理方案對多小波系數進行處理，之后進行信號重構。

圖3 多小波相鄰系數法降噪流程Fig.3 Flowchart of multiwavelets neighboring coefficients de-noising

多小波具有r個尺度函數和r個小波函數(r≥2)，多小波變換的濾波器系數為r×r矩陣。在進行多小波變換時，將形成一個多輸入多輸出系統(tǒng)。然而，我們采集到的信號通常為一維序列。為適應多小波函數的多維特性，需要對信號進行矢量化處理，使其變換成與濾波器維數相匹配的數據流，這個變換過程就稱為多小波的預處理。與此相對的，多小波重構后，需要將得到的數據序列還原為一維序列，為對應預處理的逆過程，通常稱這個還原過程為多小波的后處理。

針對r=2的情況，Strela[20-22]等將預處理方法兩種：過采樣預處理方法和Xia提出的預處理方法。

2.1.1過采樣預處理方法

過采樣預濾波屬于基于信號本身的預處理方法，是最簡單的預處理方法之一，該方法只需將原始信號乘以一個系數α，并將其作為多小波變換的第二個輸入矢量。

(7)

式中：x(i)位原始信號；s(i)為預處理后的二維信號。

2.1.2Xia法

Xia法是Xia在1996年提出的，由時間連續(xù)小波的逼近特性導出的預處理方法[23]，現(xiàn)以GHM為例對其進行說明。

設f∈v0，則有：

(8)

φ1(t)的支撐區(qū)間為[0,1]，φ2(t)的支撐區(qū)間為[0,2]。因此，可以從式(8)得到式(9)、(10)。

f(k)=φ2(1)c1,0,k-1

(9)

(10)

(12)

多小波初始系數c0,k=[c1,0,k,c2,0,k]T可根據式(11)、(12)確定。相應地，進行多小波重構后可以利用式(9)、(10)計算出數據序列，完成信號還原。

多小波相鄰系數降噪法是基于信號變換的連續(xù)性原理，即如果當前的系數包含信號的某些成分時，則其相鄰系數也可能包含這些成分。其關鍵在于相鄰系數的求解，這將對降噪效果產生直接影響。定義一個新的變量sj,k，將相鄰系數之間的關系考慮在內。

多小波變換得到的是多個二維多小波系數序列和一個二維多尺度系數序列，假設wj,k為多小波變換所得到的第j層、第k個二維多小波系數，且滿足：

wj,k=w*j,k+ej,k

(13)

式中：w*j,k表示無噪聲干擾的多小波系數；ej,k表示多小波系數中的噪聲成分，滿足二變量正態(tài)分布N(0,vj)，這里，vj表示ej,k的協(xié)方差矩陣。只有噪聲信號時，θj,k=wTj,kv-1jwj,k服從χ2n分布。

則多小波相鄰系數降噪法閾值處理過程如下：

(1)根據魯棒協(xié)方差矩陣估計方法計算vj[12]，具體方法如下：

定義：mad(y)=1.482 6*median{abs[y-median(y)]}，設變量a1、a2、b1、b2為實數，vj為2×2實數矩陣，row1與row2分別為多小波系數wj,k的第一行和第二行數據序列;

求解a1、a2、b1、b2：

a1=1.0/mad(row1)

a2=1.0/mad(row2)

b1=mad(a1*row1+a2*row2)

b2=mad(a1*row1-a2*row2)

求解j：

Vj[1][1]=1/(a1*a1)

Vj[2][2]=1/(a2*a2)

Vj[1][2]=(b1-b2)/[(b1+b2)*a1*a2]

Vj[2][1]=Vj[1][2]

(2)利用θj,k=wTj,kV-1jwj,k，求解θj,k；

(3)利用Smj,k=θmj,k-1+θmj,k+θmj,k+1將θmj,k與其相鄰的系數結合，得到包含相鄰系數信息的Smj,k，其中m為非負整數時降噪效果較好，研究表明當m=2時效果較好[12]。

(4)利用閾值函數對多小波系數進行處理，得到去噪后的多小波系數，閾值函數主要包括硬閾值函數和軟閾值函數。

硬閾值函數：

(14)

軟閾值函數：

(15)

式中：λ=2logN。

2.2 基于仿真信號的降噪算法驗證

水電機組振動頻率與機組參數存在一定關系，根據這些關系構建水導軸承處主軸徑向振動仿真信號[24]。

s(t)=f(t)+z(t)=0.39sin(πt)+0.21sin(4.2πt)+

0.1sin(33.4πt)+0.11sin(58.4πt)+0.13sin(100πt)+

0.06sin(200πt)+z(t)

(16)

構造的標準仿真信號f(t)如圖4(a)所示,仿真振動信號s(t)信噪比SNR=12，其時域波形圖如圖4(b)所示。

圖4 故障仿真信號Fig.4 Fault simulation signal

為衡量多小波相鄰系數降噪法的降噪效果，針對式(16)的仿真信號，分別使用小波閾值降噪法和多小波相鄰系數降噪法對其進行降噪處理，同時考慮軟硬閾值函數和不同多小波預處理方法對降噪效果產生的影響。

此外，為了定量評價多小波的降噪效果，選用信噪比(SNR)和均方根誤差(RMSE)兩個評價指標對比降噪效果。

信噪比的定義為：

(17)

均方根誤差的定義為：

(18)

降噪信號信噪比SNR越大，降噪信號的均方根誤差RMSE越小，則降噪信號就越接近原始信號，降噪效果越好。

分別選取Haar，DB4兩種典型的小波和GHM，CL3兩種多小波對式(16)的仿真信號進行降噪處理，采樣頻率fs=1 024 Hz,分解層數L=2。經2類方法降噪處理后，部分結果如圖5，圖6所示，其SNR和RMSE兩項指標結果如表1所示。

圖5 Db4小波降噪Fig.5 Db4 wavelets de-noising

圖6 GHM多小波過采樣預處理相鄰系數法降噪Fig.6 GHM-rr neighboring coefficients de-noising

降噪方法SNR硬閾值/h軟閾值/sRMSE硬閾值/h軟閾值/shaar12．272412．27240．08570．0857Db413．358813．35880．07560．0756Ghm＿rr15．619915．46810．05830．0593Ghm＿X13．682013．52180．07280．0742Cl＿rr15．873615．67290．05660．0579Cl＿X13．430413．33260．07280．0742

注：rr代表過采樣預濾波，X表示Xia提出的預濾波法。

可見，不同多小波和不同的多小波預處理方法都會對多小波相鄰系數法的降噪效果產生一定的影響；當選擇合適的預處理方法時，多小波相鄰系數法方能顯示出其優(yōu)越性，產生明顯優(yōu)于小波閾值降噪法的降噪效果。

3 實際信號降噪分析

本文使用某電站采集的機組振動信號，采樣頻率fs=333.3 Hz,機組轉速為500 rpm,采樣點為機組上機架X向水平振動數據，采集信號如圖7所示。

圖8為Db4小波進行小波閾值降噪后得到的振動信號，圖9為GHM多小波使用過采樣預處理方法進行相鄰系數降噪后得到的振動信號。通過對比可以看出，小波降噪方法在一定程度上消除了噪聲信號，但是存在連續(xù)性問題，并且存在有用信息被去除的問題；而多小波相鄰系數降噪法不僅濾除了噪聲，還較好保留了振動波形特征。

圖7 采集振動信號Fig.7 Factual signal

圖8 Db4小波降噪Fig.8 Db4 wavelets de-nosing

圖9 GHM多小波過采樣預處理相鄰系數法降噪Fig.9 GHM-rr neighboring coefficients de-noising

實際采集信號無法得到原始理想無噪信號，因此無法通過式(17)和式(18)進行衡量，需要其他評價指標對降噪效果進行評價。將采集信號進行小波分解后，各頻段的能量特征包含著信號時頻域的信息，基于此可以將保留信號特征的能量作為衡量去噪效果的依據[25]。先計算各頻段內的信號能量，再將各頻段的能量進行歸一化處理。將利用不同降噪方法得到的降噪信號在不同頻段的信號能量與采集信號的信號能量進行對比，以衡量其降噪效果。本文對信號進行3層小波分解，由奈奎斯特定理可知數據的頻率上限為fs/2=166.66 Hz,其頻帶范圍見表2，各頻帶能量分布見表3。

表2 三層小波分解后各頻段的頻率范圍Tab.2 Frequency range after three layer wavelet decomposition

表3 不同降噪方法得到的降噪信號能量對比Tab.3 Energy results contrast on different threshold functions

注：h表示硬閾值降噪；s表示軟閾值降噪。

從表3也可以看出，信號經小波降噪后在頻帶1中能量減少，在頻帶2中能量有所增加，在頻帶3、4高頻區(qū)能量增加明顯，造成信息丟失；信號經多小波相鄰系數法降噪后在頻帶1、2區(qū)能量增加，3、4高頻區(qū)能量減少，既濾除了噪聲，又保留了相關信息的振動特征。

同時，可以通過幾個特定頻段幅值去噪前后的差別來衡量去噪效果[26]。對原始信號和降噪信號進行傅里葉變換，運用FFT頻譜對其進行對比，比較降噪信號幅值與理想幅值的接近程度。本文選擇1倍頻，2倍頻和3倍頻3個特征量對降噪效果對比分析，機組轉速為500 rpm,則其轉頻fn=8.33 Hz。分析結果由表4可以看出，通過多小波相鄰系數降噪法去噪后，振動信號各個倍頻幅值與小波降噪相比要更接近采集的振動信號，較好保留了振動信號的特征信息，與能量分析結果一致。

表4 不同降噪方法的降噪效果對比Tab.4 De-noising results contrast on different threshold functions

對采集信號和各降噪信號的FFT頻譜分析如圖10，圖10(a)為采集信號和各將噪信號的頻譜分析圖，圖10(b)為將圖10(a)中X區(qū)域放大后的結果。通過圖10我們可以看出，各種降噪信號和采集信號的頻譜基本相同，信號的基本特征得到保留。相對于小波閾值降噪，多小波相鄰系數降噪法得到的降噪信號的曲線與采集信號的重合度更高，對信號特征保留的更完整。

圖10 FFT頻譜圖Fig.10 FFT spectrogram

經過上述各個方面的對比可以看出，經多小波相鄰系數降噪法處理后，機組振動信號中噪聲被濾除的更徹底，信號特征保留程度更高，是一種更加有效的降噪方法。

4 結 語

本文對多小波理論進行簡要介紹，對多小波相鄰系數降噪方法進行簡單分析，并將其應用于水電機組振動信號降噪過程，通過實驗室仿真信號和現(xiàn)場采集數據進行試驗，對比不同降噪方法的效果，發(fā)現(xiàn)多小波相鄰系數降噪方法可以更好地抑制噪聲，保留振動信號中的特征信息，表明多小波相鄰系數降噪法是一種更加有效的降噪方法。

□

[1] Donoho D L. De-noising by soft-thresholding [J]. Information Theory IEEE Transactions on, 1995,41(3):613-627.

[2] 孫 濤,黃天戌,孫穎慧,等. 連續(xù)小波變換識別水輪機故障信號孤立奇異點[J]. 哈爾濱工業(yè)大學學報，2003,(1).

[3] Xiangzhong Meng,Jinping Ni,Yanbo Zhu. Research on Vibration Signal Filtering Based on Wavelet Multi-resolution Analysis[C]∥ Artificial Intelligence and Computational Intelligence (AICI),2010:115-118.

[4] Zhao Ting,Yang Ping. Threshold Denoising Analysis of Machinery Vibrating Signal Based on Wavelet Transform[C]∥ Electronic Measurement and Instruments, ICEMI '07 8thInternational Conference on, 2007,(3):1 041-1 044.

[5] Isayed B, Cheded L, AI-Badour F. Vibration monitoring and faults detection using wavelet techniques[J]. Signal Processing and Its Applications, 2007. ISSPA 2009. 9th International Symposium on, 1-4.

[6] 盧 娜，肖志懷，符向前.基于蟻群初始化小波網絡的水電機組振動故障診斷[J].水力發(fā)電學報，2014，33(2)：251-258.

[7] Coifman R R, Donoho D L. Translation-invariant de-noising [J], Wavelets and Statistics, Lecture Notes in Statistics. 1995:1-26.

[8] 段 汕，何 娟，劉少英. 多小波變換在信號去噪中的應用[J]. 中南民族大學學報；自然科學版，2009,28(2)：99-103.

[9] V Strela, P N Heller, G Strang, et al. The application of multiwavelet filterbanks to signal and image processing[J]. IEEE Transactions on Image Processing, 1998,8(4):548-563.

[10] 謝榮生，李漢杰，孫 楓，等. 基于多小波噪聲方差閾值的信號濾波方法[J]. 哈爾濱工程大學學報，2002,(2).

[11] 劉志剛，錢清泉. 自適應閾值多小波故障暫態(tài)信號去噪方法[J]. 系統(tǒng)工程與電子技術，2004,(7).

[12] Chen G Y, Bui T D. Multiwavelets denoising using neighboring coefficients[J]. IEEE Signal Processing Letters, 2003,10(7)：211-214.

[13] Cai T T, Silverman B W. Incorporating information on neighboring coefficients into wavelet estimation[J]. Indian Journal of Statistics B, 2001,63(2):127-148.

[14] 徐冰雁，黃成軍，錢 勇，等. 多小波相鄰系數法在局部放電去噪中的應用[J]. 電網技術，2005,15:61-64,70.

[15] 袁 靜，何正嘉，王曉東，等. 平移不變多小波相鄰系數降噪方法及其在監(jiān)測診斷中的應用[J]. 機械工程學報，2009,(4):155-160.

[16] Xiaodong Wang, Yanyang Zi, Zhengjia He. Multiwavelet denoising with improved neighboring coefficients for application on rolling bearing fault diagnosis[J]. Machanical Systems and Signal Processing, 2011,25:285-304.

[17] Jeffery S Geronimo, Douglas P Hardin, Peter R. Massopust, Fractal function and wavelet expansions based on several scaling functions[J].Journal of Approximation Theory, 1994,78:373-401.

[18] Charles K Chui, Jian-ao Lian.A study of Orthonormal multiwavelets [J].Applied Numerial Mathematics, 1996,20:273-298.

[19] 夏國榮，徐志勝． 多小波閾值降噪法在鋼絲繩缺陷檢測中的應用[J]． 測試技術學報，2007，21(4): 311-323．

[20] V Strela, A T Walden. Orthogonal and Biorthogonal Multiwavelets for Signal Denoising and Image Compression[C]∥ Proc.SPIE, 1998,3 391:96-107.

[21] V Strela, P N Heller, G Strang P, et al. The application of multiwavelet filter banks to image processing [J]. IEEE Trans.Image Processing, 1999,(8):548-562.

[22] V Strela, A T Walden. Signal and image denoising via wavelet thresholding: Orthogonal and biorthogonal, scalar and multiple wavelet transforms[R]. Stat. Sect. Tech. Rep. TR-98-01, Dept. Math., Imperial Coll., London, U.K., 1998.

[23] Xiang-Gen Xia, Jeffrey S Geronimo, Douglas P Hardin, et al. Design of Prefilter for Discrete Multiwavelet Transforms, IEEE Trans[J]. Singnal Processing, 1996,44(1)：25-35.

[24] 趙道利，馬 薇，梁武科，等. 水電機組振動故障的信息融合診斷與仿真研究[J]. 中國電機工程學報，2004,25(20):137-142.

[25] 唐友懷，張海濤，羅 珊，等．基于小波包能量譜分析的電機故障診斷[J]．應用天地，2008，27(2)：54-57，63．

[26] 白 亮，王 瀚，李 輝，等. 基于時間序列相似性挖掘的水電機組振動故障診斷研究[J]. 水力發(fā)電學報，2010,29(9)：229-236.