基于多小波的水电机组振动信号降噪方法研究

张广涛，宋丽波，张晓静，孙召辉，梁红梅，肖志怀，程远楚

(1.武汉大学动力与机械学院，武汉 430072；2.国网新源山东泰山抽水蓄能电站有限责任公司,山东 泰安 271000；3.新疆昌吉职业技术学院电气工程分院，新疆 昌吉 831100)

0 引 言

据统计，在水电机组故障或事故中，约80%在振动信号中有所反映；同时，异常振动也是设备损坏的一项重要原因。因此，振动信号分析方法在水电机组故障诊断领域获得了广泛应用。但是，水电机组运行工况复杂，运行环境恶劣，现场振源丰富，实际采集的振动信号往往被淹没在大量的背景噪声中，给振动信号故障特征提取造成了不良影响。

为了有效地去除噪声，同时考虑到振动信号的非平稳特征，众多学者将小波分析作为信号降噪的主要方法。Donoho[1]首次提出了硬阈值和软阈值降噪方法。孙涛[2]等采用连续小波变换系数模极大值法，对水轮机低频振动信号中故障点的奇异性程度进行估算并取得了很好的诊断效果。Xiangzhong Meng[3]等将小波多分辨率分析引入振动信号滤波中，运用Mallat算法对振动信号进行降噪处理，转子平台实验表明，该方法适用于振动信号分析，降低了计算量，提高了计算效率。Zhao Ting[4]等应用一种基于小波变换的改进阈值降噪方法，通过MATLAB仿真实验证明了当机械振动信号受到高斯白噪声严重干扰时，此方法比其他降噪方法更有效，可以有效是地消除噪声并保留信号的详细信息。Isayed B[5]等将小波技术应用于振动监测和故障检测中，使用小波技术分析旋转机械振动信号，通过小波包变换揭示信号的瞬态信息，利用细节信号中的能量平均值从振动信号中对故障进行有效检测。Coifman R R[6]等针对小波降噪过程中在不连续点可能出现的Gibbs现象，提出了平移不变小波变换方法，采用不同平移量对原始信号进行的多次“平移-降噪-逆平移”处理，然后对所有降噪结果取平均值，消除小波平移依赖性。

尽管与小波分析相关的信号处理方法获得了极大的发展，但是由于小波分析自身的缺陷 以及信号处理技术飞度发展，多小波分析方法应运而生，多小波同时满足对称性，短(紧)支撑性，二阶消失矩和正交性这些特性，可以匹配信号中不同特征。近些年来，多小波在图像处理、数据压缩和故障诊断等领域都获得了极大的发展，并确立了相对单小波的优势[7]。

Strela[8]等将多小波单变量阈值降噪方法应用于二维图像信号的降噪，取得了良好的效果。谢荣生[9]等对信号和噪声的多小波变换进行了研究，在多小波噪声方差阈值的基础上，提出了一种新的多小波信号滤波方法，仿真分析表明，该方法具有很好的滤波效果。刘志刚[10]等探讨了基于多小波的信号去噪方法，提出了一种基于自适应阈值的多小波去噪方法，将其应用于电力系统输电线路故障暂态信号的去噪，仿真结果表明，该方法可以根据实际信号自适应改变阈值的大小，在去噪效果上优于传统多小波去噪方法。这些方法都具有良好的效果，但是忽视了多小波系数间的相关性。文献[11,12]引入了相邻系数的概念，提出了多小波相邻系数降噪方法，通过与其他方法的比较，显示该方法可以取得更好的效果。

尽管多小波相邻系数降噪方法已证明其优越性，也在一些领域得到了应用[13-15]。但在水电机组振动信号降噪方面的应用还很少见。针对水电机组振动信号的特点，本文将围绕多小波相邻系数降噪方法完成对振动信号的降噪处理，针对硬阈值和软阈值两种情况和不同多小波类型和预处理方法，将降噪效果和当前广泛使用的小波阈值降噪方法进行对比。通过降噪前后的定量分析，表明当使用合适的多小波与处理方法时，多小波相邻系数降噪方法在水电机组振动信号降噪方面具有良好的效果。

1 多小波理论

1.1 多小波多分辨率分析

多小波理论是在小波理论的基础上发展而来的，同小波一样，多小波同样满足多分辨率分析。不同之处在于多小波的多尺度函数是以多个尺度函数{φ，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多小波。该多小波包含两个尺度函数和两个小波函数，多尺度函数的支撑区间分别为[0,1]和[0,2]，多小波函数的支撑区间为[0,2]。

GHM多小波的滤波器系数矩阵为：

1.2.2CL3多小波

CL多小波是由Chui C K和Lian J[18]于1995年提出的。考虑不同长度的支撑区间，又可以分为CL3和CL4两类。其中，CL3多小波具有正交性、对称性、二阶消失矩和紧支撑的特性，其尺度函数和小波函数的支撑区间为[0,2]，滤波器长度为3；而CL4多小波的尺度函数和小波函数的支撑区间为[0,3]，滤波器长度为4。为方便与GHM多小波对比，在此选取CL3多小波。

CL3多小波的滤波器系数矩阵为：

2 多小波相邻系数法研究

多小波分析是在小波分析的基础上发展而来的，是对小波分析的继承和发展，。与小波阈值降噪类似，多小波降噪仍然是基于如下理论展开：假定故障信号经过小波分解，信号特征成分的能量会集中在少数大的小波系数上，而噪声分解后的小波系数一般都非常小。因此，经过小波分解后，信号的小波变换系数要大于噪声的小波变换系数，选择合适的λ作为阈值，当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矩阵。在进行多小波变换时，将形成一个多输入多输出系统。然而，我们采集到的信号通常为一维序列。为适应多小波函数的多维特性，需要对信号进行矢量化处理，使其变换成与滤波器维数相匹配的数据流，这个变换过程就称为多小波的预处理。与此相对的，多小波重构后，需要将得到的数据序列还原为一维序列，为对应预处理的逆过程，通常称这个还原过程为多小波的后处理。

针对r=2的情况，Strela[20-22]等将预处理方法两种：过采样预处理方法和Xia提出的预处理方法。

2.1.1过采样预处理方法

过采样预滤波属于基于信号本身的预处理方法，是最简单的预处理方法之一，该方法只需将原始信号乘以一个系数α，并将其作为多小波变换的第二个输入矢量。

(7)

式中：x(i)位原始信号；s(i)为预处理后的二维信号。

2.1.2Xia法

Xia法是Xia在1996年提出的，由时间连续小波的逼近特性导出的预处理方法[23]，现以GHM为例对其进行说明。

设f∈v0，则有：

(8)

φ1(t)的支撑区间为[0,1]，φ2(t)的支撑区间为[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)计算出数据序列，完成信号还原。

多小波相邻系数降噪法是基于信号变换的连续性原理，即如果当前的系数包含信号的某些成分时，则其相邻系数也可能包含这些成分。其关键在于相邻系数的求解，这将对降噪效果产生直接影响。定义一个新的变量sj,k，将相邻系数之间的关系考虑在内。

多小波变换得到的是多个二维多小波系数序列和一个二维多尺度系数序列，假设wj,k为多小波变换所得到的第j层、第k个二维多小波系数，且满足：

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

(13)

式中：w*j,k表示无噪声干扰的多小波系数；ej,k表示多小波系数中的噪声成分，满足二变量正态分布N(0,vj)，这里，vj表示ej,k的协方差矩阵。只有噪声信号时，θj,k=wTj,kv-1jwj,k服从χ2n分布。

则多小波相邻系数降噪法阈值处理过程如下：

(1)根据鲁棒协方差矩阵估计方法计算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提出的预滤波法。

可见，不同多小波和不同的多小波预处理方法都会对多小波相邻系数法的降噪效果产生一定的影响；当选择合适的预处理方法时，多小波相邻系数法方能显示出其优越性，产生明显优于小波阈值降噪法的降噪效果。

3 实际信号降噪分析

本文使用某电站采集的机组振动信号，采样频率fs=333.3 Hz,机组转速为500 rpm,采样点为机组上机架X向水平振动数据，采集信号如图7所示。

图8为Db4小波进行小波阈值降噪后得到的振动信号，图9为GHM多小波使用过采样预处理方法进行相邻系数降噪后得到的振动信号。通过对比可以看出，小波降噪方法在一定程度上消除了噪声信号，但是存在连续性问题，并且存在有用信息被去除的问题；而多小波相邻系数降噪法不仅滤除了噪声，还较好保留了振动波形特征。

图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高频区能量增加明显，造成信息丢失；信号经多小波相邻系数法降噪后在频带1、2区能量增加，3、4高频区能量减少，既滤除了噪声，又保留了相关信息的振动特征。

同时，可以通过几个特定频段幅值去噪前后的差别来衡量去噪效果[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区域放大后的结果。通过图10我们可以看出，各种降噪信号和采集信号的频谱基本相同，信号的基本特征得到保留。相对于小波阈值降噪，多小波相邻系数降噪法得到的降噪信号的曲线与采集信号的重合度更高，对信号特征保留的更完整。

图10 FFT频谱图Fig.10 FFT spectrogram

经过上述各个方面的对比可以看出，经多小波相邻系数降噪法处理后，机组振动信号中噪声被滤除的更彻底，信号特征保留程度更高，是一种更加有效的降噪方法。

4 结 语

本文对多小波理论进行简要介绍，对多小波相邻系数降噪方法进行简单分析，并将其应用于水电机组振动信号降噪过程，通过实验室仿真信号和现场采集数据进行试验，对比不同降噪方法的效果，发现多小波相邻系数降噪方法可以更好地抑制噪声，保留振动信号中的特征信息，表明多小波相邻系数降噪法是一种更加有效的降噪方法。

□

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

[2] 孙 涛,黄天戌,孙颖慧,等. 连续小波变换识别水轮机故障信号孤立奇异点[J]. 哈尔滨工业大学学报，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].水力发电学报，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] 刘志刚，钱清泉. 自适应阈值多小波故障暂态信号去噪方法[J]. 系统工程与电子技术，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] 袁 静，何正嘉，王晓东，等. 平移不变多小波相邻系数降噪方法及其在监测诊断中的应用[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]. 水力发电学报，2010,29(9)：229-236.