基于样本熵的改进小波包阈值去噪算法*

向北平， 周 建， 倪 磊， 艾攀华

(西南科技大学制造过程测试技术教育部重点实验室 绵阳，621000)

引 言

长期以来，小波包滤波是对信号进行分析预处理的主要工具，该方法对信号进行小波包变换，然后利用小波包阈值函数对小波包系数进行阈值收缩处理以达到去噪目的。然而由于含噪信号中的噪声分布往往不均匀，传统的软、硬阈值方法对信号小波包系数进行固定格式的阈值处理，无法很好地满足信号去噪要求[1]。另外，在阈值的选择方面，常用的Heursure阈值、Donoho阈值等能够对信噪比较高的信号实现噪声与信号的最优分离，而对于强噪声信号，去噪效果并不理想[2]。近年来，针对这两个问题，许多学者进行了研究。

文献[3]提出了一种自适应对数小波阈值函数去噪算法，结合自适应对数阈值函数对每一层小波系数设置最优阈值，且应用于动压信号，增加了信号信噪比且减少了计算时间，但其去噪算子形式依然固定不变。文献[4]分析了传统软、硬阈值方法的局限，对阈值函数进行改进，使其具有能量分布自适应性，但信号小波包系数的能量分布并不能准确表达小波包系数的噪声分布，因此该方法仍然存在一些不足。文献[5]对Donoho阈值方法中的噪声标准差估计方法进行改进，进行仿真且取得了良好的去噪效果。Richman等[6]提出一种新的时间序列复杂度表征参数即样本熵(sample entropy,简称SE)算法，被广大学者所关注且近年来被常用于机械信号分析与故障诊断领域[7]。

基于以上分析，笔者提出将样本熵与小波包阈值去噪算法相结合，且将其应用于高速深沟球轴承振动信号去噪分析，通过分析信号的小波包系数噪声分布情况及其对应的样本熵值，使阈值去噪算子具有噪声分布自适应性以达到最优去噪效果；同时以噪声估计信号样本熵值为基准，提出了一种最优阈值估计方法。仿真分析以及实验结果皆对所提方法进行了验证。

1 样本熵算法及其噪声表征

1.1 样本熵的计算

设有长度为N的时间序列Xi={x1,x2,…,xN}，其样本熵的计算方法如下：

1) 确定嵌入维数为m，对Xi的元素按顺序进行排列，即可得到一组维数为m的向量{xm(1),…,xm(Ν-m+1)}，且

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

(1≤i≤N-m+1)

(1)

2) 定义向量Xm(i)与Xm(j)之间的间隔d[Xm(i),Xm(j)]为两向量之间对应元素求差的绝对值的最大值，即

d[Xm(i),Xm(j)]=

maxk=0,…,m-1(|x(i+k)-x(j+k)|)

(2)

3) 对于固定的Xm(i)，统计Xm(i)与Xm(j)之间距离小于等于参数r的j(1≤j≤N-m,j≠i)的个数，且记为Bi，则当1≤i≤N-m时定义

(3)

4) 定义B(m)(r)为

(4)

(5)

6) 定义Am(r)为

(6)

据上述分析可知，B(m)(r)是两个序列在相似容限r下匹配m个点的概率，而Am(r)是两个序列匹配m+1个点的概率。则该时间序列样本熵定义为

(7)

实际信号中N无法趋近于无穷，因此可将样本熵设为

(8)

上述算法中的嵌入维数m和相似容限r通常取为m=1或2，r=0.1，Std～0.25Std(Std为原始数据Xi={x1,x2,…,xN}的标准差)。文中取m=2，r=0.2Std。

1.2 样本熵与噪声

当信号受到噪声干扰后，其状态取值不确定性增加，即信号无序程度与复杂程度增加。而样本熵作为信号复杂度表征参数，当信号中噪声增加，样本熵值也应增加。为验证以上分析，设定1 kHz采样率与1 000总采样数的仿真信号f(t)=y(t)+n(t)，式中n(t)为标准差σ∈[0,3]的带限500 Hz高斯白噪声信号。

y(t)=5sin(10πt)sin(2πt) (0≤t≤1)

(9)

分析信号f(t)的样本熵值在不同采样数下与噪声标准差的关系，其结果如图1(a)所示(图中100采样点为信号前100采样点，下同)。由于实际工程信号中的噪声很少为单纯的白噪声，因此，利用反幂律滤波器对功率谱密度分布均匀的白噪声n(t)上色，选定反频谱指数为1，此时理想数字滤波器的幅度平方响应为1/f1。由于该滤波器处理后导致噪声幅值衰减，因此，对加色后的噪声信号进行10倍处理，得到有色噪声n1(t)及含噪仿真信号f1(t)=y(t)+n1(t)。设n2(t)为1 kHz采样率，1 000采样数的幅值区间为[0,3]的带限500Hz周期性随机噪声，得到含噪仿真信号f2(t)=y(t)+n2(t)。同样分析信号f1(t)与f2(t)的样本熵值随噪声大小的变化关系，其结果分别如图1(b)，(c)所示。

图1 样本熵与几种常见噪声大小关系Fig.1 Correlation between sample entropy and several common noise

由图1(a)看出，信号的样本熵值与噪声标准差成正相关，且在采样点相差巨大的情况下其样本熵值变化依然相近。图1(b)，(c)较图1(a)而言，样本熵随噪声标准差变化有少许波动，熵值大小区间有变化，但总体趋势仍呈正相关，且数据长度对其趋势影响不大。这说明，虽然对仿真信号施加不同种类的噪声，其样本熵变化曲线有所不同，但由于对信号增加噪声的结果导致了信号的随机度与复杂度增加，其结果仍然是样本熵增加。图1所示说明了样本熵算法可用来表征时间序列的含(多种)噪声情况且受数据长度影响较小。

2 改进小波包阈值去噪算法

2.1 传统小波包去噪及其不足

传统的阈值函数主要有软阈值函数和硬阈值函数两种。软阈值函数

(10)

硬阈值函数

(11)

工程实践中常用的阈值函数还有形如文献[8]中提出的一种介于软、硬阈值函数之间的改进阈值函数

(12)

由以上函数式可知：硬阈值由于其函数不连续性可能产生伪Gibbs现象且导致有用信号缺失；软阈值函数虽然连续性好，但存在较大偏差[9-10]。文献[8]阈值函数虽进行了改进，但只是软硬阈值函数的折中算法，通过下文分析可知，该阈值函数仅为本研究提出的自适应阈值函数的一种取值情况(即当调节参数s=0.5时)。在实际应用中，阈值函数的选择并无固定标准，对于不同的信号选用不同的去噪算子达到的去噪效果也不同[11]，因此，研究一种能依据信号小波包系数的含噪情况自适应调整的新阈值函数具有实际的工程意义。

2.2 自适应改进阈值函数

根据以上分析，笔者对传统阈值函数作改进，使其不受限于固定去噪形式，得到的改进阈值函数为

(13)

其中：0

由于阈值函数偏硬时对大于阈值的小波包系数保留较好，而偏软时对小波包系数压缩更大。因此，对于含噪较多的小波包系数，阈值处理方式应偏软即s较大；而受噪声影响小的小波包系数阈值函数应偏硬，即s应较小。

由1.2节分析可知，样本熵能较好地反应时间序列的噪声变化情况，且适用于分析短时间序列，因此设定调节参数s的确定方法如下：

1) 对信号的小波包系数序列(w1,w2,…,wn)按顺序分割为n-l个相互之间最大重叠且长度均为l的子序列(k1,k2,…,kn-l)，即ki向后移动一位数据得到ki+1；

3) 将该序列进行极值归一化后带入改进阈值函数中，可使其具有小波包系数噪声分布自适应性。

2.3 阈值估计

(14)

图2 阈值与的关系Fig.2 Correlation between threshold and

如图2所示，按本方法计算得到的阈值为λ=0.6，为进行对比，设定其他阈值分别为λ=0.5,0.7，利用文中自适应改进阈值函数(sym8小波，分解层数为3层)分别选定上述3个阈值进行去噪分析，得到结果如图3所示。由图3可知，当阈值λ=0.5时，阈值过小，信号毛刺较多，噪声去除不完全；当阈值取0.7时，原始信号被过度压缩；而通过本方法确定的阈值取得了更好的去噪效果，证明了该阈值估计方法的有效性。

图3 不同阈值去噪效果对比Fig.3 Denoising effect comparison of different threshold

为了验证该去噪方法受噪声标准差的影响，同样分析上述信号，分别加入不同标准差的混合噪声，对其进行去噪，计算去噪前后的信噪比(signal to noise ratio,简称SNR)与均方根误差(root mean square error,简称RMSE)，且将结果列于表1。由表1可知，对不同噪声标准差情况下的含噪信号进行去噪，去噪后信号的SNR与RMSE变化幅度不大，证明该方法受噪声标准差影响较小，适用于对不同噪声尺度情况下的信号进行去噪。

表1 不同噪声尺度下的去噪结果

3 实 验

3.1 振动实验装置与信号获取

为了对本算法进一步验证，搭建了轴承振动实验平台如图4所示，该实验平台主要包括中间的小型高速直流电机，两端测试轴承为分子泵中使用的微小型高速深沟球轴承QC0011286(轴承参数为：内径为4 mm，外径为13 mm，轴承节圆直径D为8.5 mm，滚珠直径d为2.3 mm，滚珠个数N为7个，接触角β为0°)以及固定于轴承外圈的加速度传感器。测试时电机转速设定为60 kr/min，采样率为20 kHz，采样时间为0.1 s，采样时前置加上10 kHz的低通抗混滤波。为模拟轴承外圈故障，在轴承外圈内壁加工宽为0.15 mm、深为0.2 mm的横向沟槽，此时轴承基频为fr=1 kHz，轴承外圈故障频率根据公式(15)可得

(15)

采集到的轴承振动时域波形如图5所示。由图5可知，从该时域信号无法直观地得出任何轴承信号特征。

图4 轴承振动实验平台Fig.4 Bearing vibration experiment platform

图5 轴承振动时域信号Fig.5 Bearing vibration time domain signal

3.2 去噪分析与对比

对图5所示的信号选用sym8母小波进行四层小波包分解，获得最大尺度上的小波包系数序列w(i)，为进行样本熵计算，对其分割为若干个子序列，为在不失统计分布性的前提下尽可能准确地表达小波包系数噪声分布情况，进行多次实验，最终取子序列长度为l=127，求得小波包系数序列所对应的样本熵序列，且对其进行归一化即可得到阈值函数调节参数序列s如图6所示。由图6可知，该小波包系数序列中的噪声分布不均匀，因此在对低频与中频小波包系数去噪时，调节参数较小，去噪形式偏硬，而对其他含噪较多(调节参数较大)的系数则去噪形式偏软。

图6 阈值函数调节参数Fig.6 Adjust parameters of threshold function

图7 不同去噪算法结果对比Fig.7 Denoising results comparison of different algorithm

根据以上分析利用本方法对实验信号进行去噪，为进行直观对比，将去噪前后信号进行功率谱分析且与其他阈值去噪方法相对比得到如图7所示的结果。由于轴承外圈故障而产生周期性脉冲信号，其表现形式为调制信号，即以轴承的转动频率为调制频率，外圈故障频率为载波频率，形成边频带。通常由于故障特征微弱且存在噪声干扰，轴承故障特征及其调制特征无法清晰显露。通过信号去噪处理后，利用功率谱分析得到轴承振动信号频谱，且由于不同的去噪方法得到的功率谱分析效果不同，因此为了直观地分析轴承振动信号去噪效果，以去噪后信号特征频率成分的还原情况作为去噪效果评判标准，具体分析如下。

由图7(a)可知，原始含噪轴承振动信号功率谱中仅能分辨出轴承基频1 kHz及其倍频2 kHz的频率成分，其他有用的频率成分被噪声所淹没；通过自相关去噪法去噪后得到的结果如图7(b)所示，由图可知，自相关去噪法去除了部分噪声，还原了信号的3 kHz频率成分，但故障频率依然难以分辨，这表明该信号所含噪声并非单纯的高斯白噪声，想要得到更好的去噪结果，有必要考虑其他方法；通过软、硬阈值函数且选用基于Stein无偏似然估计的阈值确定规则(文献[8]同样采取该规则)，得到去噪结果如图7(c)与(d)所示，硬阈值函数去噪后，能够识别轴承故障特征频率2 552.9 Hz，但在故障频率周围存在着无效频率成分，无法判断该故障特征的有效性，而软阈值函数虽然滤除了大部分噪声成分，但原始频率成分如基频1 kHz，故障特征频率2 552.9 Hz也被扼杀严重；文献[8]阈值函数去噪结果如图7(e)所示，该结果较软硬阈值函数更好，信号的基频及其倍频1,2,3 kHz与轴承外圈故障频率foc=2 552.9 Hz及其调制成分3 552.9 Hz被很好的还原，然而在3 500 Hz与高频区域仍然存在一些无关频率成分，去噪效果有待提升；图7(f)显示自适应阈值函数与阈值估计方法去除噪声明显且有效地还原了原始信号的频率特征(基频及其谐波1,2，3 kHz；故障频率及其倍频与调制成分2 552.9,3 552.9,5 150.8 Hz)，信号失真较少，去噪效果较好，提升了轴承故障诊断的准确性。

事实上，通过对比信号去噪前后0.01 s的细致波形(如图8所示)也可以直观地发现，信号去噪前其波形毛刺较多，即存在较多高频噪声干扰，而去噪后信号细致波形毛刺较少，信号波形较为规律。

通过实验发现，本算法虽然去噪效果优于传统去噪算法，但用时较长，计算速度较传统算法更慢，不适用于在线实时去噪分析与故障诊断。因此，在实际应用中，可利用传统的软、硬阈值去噪算法作为在线实时初步分析工具，而本算法用于离线的进一步精细分析。

图8 去噪前后信号细致波形对比Fig.8 Detail waveform comparison of noisy and denoised signal

4 结 论

1) 样本熵可以表征信号中不同种类的噪声含量的大小，样本熵越大，则信号含噪越多。将样本熵应用于小波包系数序列中，能得到其噪声分布情况，据此对阈值函数进行自适应调整，使得其能够对含噪较多的小波包系数进行大尺度压缩，而使含信号较多的小波包系数得到尽可能的保护。

3) 利用本方法对滚动轴承振动实验信号进行分析，能够获得较其他方法更好的去噪效果，且有效地还原了信号的转动特征频率与故障特征频率，提高了轴承状态监测的准确度。

4) 本方法去噪效果较传统算法更好，但计算速度较慢，实际应用中可将传统算法作为在线初步去噪算法，本方法作为离线精细算法进行配合分析。