基于STFT时频谱系数收缩的信号降噪方法*

2015-06-13郭远晶魏燕定周晓军

振动、测试与诊断 2015年6期

关键词：时域步长频谱

郭远晶，魏燕定，周晓军

(浙江大学浙江省先进制造技术重点研究实验室杭州，310027)

基于STFT时频谱系数收缩的信号降噪方法*

郭远晶，魏燕定，周晓军

(浙江大学浙江省先进制造技术重点研究实验室杭州，310027)

针对旋转机械故障振动信号的降噪问题，提出一种基于短时Fourier变换(short time Fourier transform,简称STFT)时频谱系数收缩的信号降噪方法。先将信号进行STFT，得到其时频谱。由于谱系数为复数，故根据模值大小进行谱系数收缩，并利用步长迭代算法在0到谱系数最大模值的区间内估计最优阈值。迭代运算过程中，首先，分别采用基本的硬阈值函数和软阈值函数进行系数收缩；然后，以改进风险函数为阈值评价标准，估计最优阈值；最后，利用最优阈值重新进行谱系数收缩，对得到的新谱进行STFT逆变换，重构降噪后的时域信号。仿真信号与试验数据的处理结果表明，利用所估计的最优阈值，STFT时频谱系数硬、软阈值函数收缩方法均能够实现噪声混合信号的降噪。

故障诊断；信号降噪；短时Fourier变换；步长迭代算法；改进风险函数；最优阈值估计；谱系数收缩

引言

对于旋转机械故障振动信号的降噪，较为实用有效的方法是阈值去噪，其中小波阈值去噪的研究与应用最为广泛。Donoho等[1-2]提出了基于小波变换系数的“VisuShrink”阈值去噪方法，引入了硬阈值函数和软阈值函数这两种基本的系数收缩函数(系数收缩函数又可称为阈值函数)，并建立了通用的阈值估计公式。随后，又提出了小波阈值去噪的“SUREShrink”方法[3]以及Minimax阈值估计方法[4]。国内外学者提出了各种形式的阈值函数以及阈值估计方法[5-10]，改进了小波阈值去噪的性能，扩展了其应用领域，并研究了信号阈值去噪另一种实现方式的联合时频域去噪。文献[11-14]借鉴了小波阈值去噪的思想，研究了基于Gabor变换时频谱系数的阈值去噪方法。Parolai[15]针对地震图的消噪处理，提出了基于S变换谱系数的自适应阈值函数的去噪方法，但其主要是现有的小波阈值函数和通用阈值估计公式的应用。Hon等[16]针对超声弹性成像，提出了基于短时Fourier变换谱系数阈值滤波的图像去噪算法，该算法本质上是时频滤波。联合时频域阈值去噪研究的核心和难点依然是设计合适的阈值函数以及估计最优阈值，以期获得最大限度接近理想无噪信号的估计信号。

笔者针对一维时域信号的降噪，提出基于STFT时频谱系数收缩的信号降噪方法。该方法先利用STFT获取一维时域信号的时频谱系数，然后将谱系数按模值大小、根据阈值函数进行收缩，最后利用STFT逆变换，重构得到降噪后的时域信号。系数收缩过程中，阈值函数可以为硬阈值函数、软阈值函数或者其他的阈值函数。对于最优阈值估计，以所提出的改进风险函数为评价指标，采用步长迭代算法获取。

1 STFT

给定一个时间宽度很短的窗函数γ(t)，令窗滑动，则连续时间信号x(t)的STFT定义[17]为

(1)

其中：“*”表示复数共轭。

STFT既是时间的函数，又是频率的函数，它可以看作是信号x(t)在分析时刻τ附近的局部频谱。

STFT是一种可逆变换。若窗函数γ(t)满足完全重构条件

(2)

即能量归一化条件，则STFT逆变换定义为

(3)

对于任何实际应用而言，都需将STFT(τ,f)离散化。考虑STFT在等间隔时频网格点(mT,nF)处采样，其中T>0和F>0分别是时间变量和频率变量的采样周期，m和n为整数。为了简便，记STFT(mT,nF) = STFT (m,n)，则很容易得到离散时间信号x(k)的STFT离散化形式及STFT逆变换的离散化形式

(4)

(5)

2 STFT时频谱系数收缩的信号降噪方法

在STFT时频谱系数收缩的信号降噪方法实现过程中，首先需要将采样得到的长度为N的离散时间信号x(k)进行STFT变换，得到其时频谱STFT(m,n)。考虑到STFT(m,n)为复数矩阵，如果对谱系数的实部和虚部分别进行收缩，则必然改变系数的相位，从而使STFT逆变换重构得到的时域信号中产生附加噪声，因此，笔者根据模值大小进行系数收缩。

1) 将目标模值区间或者模值系数目标区间分成M个步长，则每个系数步长Δd=(β-α)/M，对应的阈值步长ΔD=Δdmax{|STFT(m,n)|}。阈值变量δth初始化为αmax{|STFT(m,n)|} 。

2) 使阈值变量δth按阈值步长ΔD增大，每增大一个步长，获得一个新阈值δnew。

3) 以δnew为阈值(初始值为αmax{|STFT(m,n)|})，根据模值大小对谱系数进行收缩。其中，阈值函数可以为硬阈值函数、软阈值函数或者其他的阈值函数。最后对系数收缩后得到的新谱进行STFT逆变换，重构时域信号。

假设信号x(k)中原始无噪的信号为s(k)，并混有加法性噪声u(k)，则x(k)表示为

x(k)=s(k)+u(k)

(6)

令STFT时频谱经过系数收缩后得到的新谱为STFT′[m,n]，那么新谱经STFT逆变换重构得到的时域信号为

(7)

(8)

但是，在处理实际信号时，原始无噪信号s(k)是不可知的，L2-RISK无法计算，因此，笔者提出一种改进风险函数RM(r)，其定义为

(9)

这里，xr(k)为阈值变量δth增加了r个阈值步长ΔD，即δth=αmax{|STFT(m,n)|}+rΔD时，经过STFT时频谱系数收缩后重构得到的时域信号。这样，RM(r)就定义为第r步迭代计算得到的时域信号xr(k)与第r-1步迭代计算得到的时域信号xr-1(k)的均方误差。因此，重构得到时域信号后，进一步计算改进风险函数RM(r)。

4) 重复步骤2和3，直到阈值变量δth增加到βmax{|STFT(m,n)|}，此时可获取改进风险函数RM(r)关于步数M或阈值变量δth的关系曲线。

3 仿真研究

为验证STFT时频谱系数收缩方法对于信号降噪的有效性，构造仿真信号x(t)进行试验。仿真信号x(t)由频率为余弦信号分量s(t)叠加噪声分量u(t)得到，其信噪比为0，即

(10)

仿真振动信号x(t)与其中的余弦信号分量s(t)分别如图1、图2所示。

图1 仿真振动信号x(t)

图2 余弦信号分量s(t)

下面对信号x(t)进行基于STFT时频系数收缩的降噪，以获取余弦信号分量s(t)的最佳估计信号s*(t)。

(11)

软阈值函数的表达式为

(12)

对于硬阈值函数的STFT时频谱系数收缩，RM(r)关于步长数M的关系曲线如图3所示。

图3 仿真振动信号硬阈值函数系数收缩得到的改进风险函数RM(r)曲线

对于软阈值函数的STFT时频谱系数收缩，RM(r)关于步长数M的关系曲线如图5所示。

仿真试验中，对比图2、图4和图7可知，基于STFT时频谱系数硬、软阈值函数收缩的信号降噪方法均能够从噪声混合信号x(t)中获取余弦信号分量s(t)的最佳估计信号s*(t)。虽然这两种方法得到的s*(t)幅值存在一定程度上的失真，且后者幅值失真更为严重，但这主要是由阈值函数所决定

图4 硬阈值函数系数收缩得到余弦信号分量s(t)的最佳估计信号s*(t)

Fig.4 The optimal estimated signals*(t)of cosine signal components(t) obtained by hard threshold coefficients shrinkage

图5 仿真振动信号软阈值函数系数收缩得到的改进风险函数RM(r)曲线

Fig.5 Emulational vibration signal modified risk functionRM(r) curve obtained by soft threshold function coefficients shrinkage

图6 由仿真振动信号改进风险函数RM(r)获得的差分谱dRM(r)曲线的，而作为最重要的信号频率信息以及信号形态，两种方法均能理想地提取。因此，仿真试验结果验证了基于STFT时频谱系数收缩方法对于信号降噪的有效性。

Fig.6 Difference spectrum dRM(r) curve obtained from emulational vibration signal modified risk function RM(r)

图7 软阈值函数系数收缩得到的余弦信号分量s(t)最佳估计信号s*(t)

4 试验数据降噪处理

为进一步验证基于STFT时频谱系数收缩方法对于实际故障振动信号降噪的有效性与实用性，采用美国PHM Society 2009年故障预测与健康管理挑战赛的齿轮箱相关故障数据进行试验[19]。以斜齿轮为试验齿轮，选择齿轮箱的两种故障状况：输入轴弯曲和24T齿轮的一个齿上有破损。试验时的齿轮箱输入转速为3 kr/min，重载。振动信号原始采样频率为66.666 7 kHz，再以6.666 7 kHz的频率对采集到的信号数据进行重采样。

对于第1种故障状况，即输入轴弯曲，所采集到的待分析振动信号x1(t)如图8所示。

图8 输入轴弯曲故障振动信号x1(t)

信号x1(t)的STFT时频谱系数最大模值max{|STFT(m,n)|}=0.060 3。在运用步长迭代算法的过程中，取模值系数区间[α,β]=[0.2,0.7]，并将该区间分成M=100个步长。去噪过程中同样分别采用硬阈值函数和软阈值函数。

对于硬阈值函数去噪，RM(r)关于步长数M的关系曲线如图9所示。

图9 故障振动信号x1(t)硬阈值函数系数收缩得到的改进风险函数RM(r)曲线

对于软阈值函数去噪，RM(r)关于步长数M的关系曲线如图11所示。

对于第2种故障状况，即24T齿轮的一个齿上有破损，所采集到的待分析振动信号x2(t)如图14所示。信号x2(t)的STFT时频谱系数最大模值max{|STFT(m,n)|}=0.028 4。在运用步长迭代算法的过程中，取模值系数区间[α,β]=[0.3,0.8]，并将该区间分成M=100个步长。去噪过程中同样分别采用硬阈值函数和软阈值函数。

对于硬阈值函数去噪，RM(r)关于步长数M的关系曲线如图15所示。

图11 故障振动信号x1(t)软阈值函数系数收缩得到的改进风险函数RM(r)曲线

Fig.11 Fault vibration signalx1(t) modified risk functionRM(r)curve obtained by soft threshold function coefficients shrinkage

图12 由故障振动信号x1(t)改进风险函数RM(r)获得的差分谱dRM(r)曲线

Fig.12 Difference spectrum dRM(r) curve obtained from fault vibration signalx1(t) modified risk functionRM(r)

图13 故障振动信号x1(t)经由软阈值函数系数收缩得到的时域降噪信号

图14 轮齿破损故障振动信号x2(t)

图15 故障振动信号x2(t)硬阈值函数系数收缩得到的改进风险函数RM(r)曲线

对于软阈值函数去噪，RM(r)关于步长数M的关系曲线如图17所示。

图17 故障振动信号x2(t)软阈值系数收缩得到的改进风险函数RM(r)曲线

Fig.17 Fault vibration signalx2(t) modified risk functionRM(r) curve obtained by soft threshold function coefficients shrinkage

图18 故障振动信号x2(t)改进风险函数RM(r)获得的差分谱dRM(r)曲线

Fig.18 Difference spectrum dRM(r) curve obtained from bearing fault vibration signalx2(t) modified risk functionRM(r)

图19 故障振动信号x2(t)经由硬阈值函数系数收缩得到的时域降噪信号

5 结论

2) STFT时频谱系数收缩方法能够从噪声混合信号中恢复时域降噪信号。将此方法应用于齿轮箱相关故障振动信号的处理，可以有效提取出与故障特征信号一致的时域降噪信号，从而指导相关故障的诊断。

[1] Donoho D L， Johnstone I M． Ideal spatial adaptation by wavelet shrinkage[J]．Biometrika，1994，81(3)：425-455．

[2] Donoho D L． De-noising by soft-thresholding [J]．IEEE Transactions on Information Theory，1995，41(3)：613-627．

[3] Donoho D L，Johnstone I M．Adapting to unknown smoothness via wavelet shrinkage[J]． Journal of the American Statistical Association，1995，90(432)：1200-1224．

[4] Donoho D L，Johnstone I M．Minimax estimation via wavelet shrinkage[J]． The Annals of Statistics，1998，26(3)：879-921．

[5] 曲天书，戴逸松，王树勋．基于SURE无偏估计的自适应小波阈值去噪[J]．电子学报，2002，30(2)：266-268．

Qu Tianshu，Dai Yisong，Wang Shuxun．Adaptive wavelet thresholding denoising method based on SURE estimation[J]．Acta Electronica Sinica，2002，30(2)：266-268．(in Chinese)

[6] Nezamoddin N K，Paul F，Edward J．BayesShrink ridgelets for image denoising [C]∥Proceedings of International Conference on Image Analysis and Recognition (ICIAR) . Verlag Berlin Heidelberg, Germany:Springer, 2004：163-170．

[7] 李方，李友荣，王志刚．基于Morlet小波与最大似然估计方法的降噪技术[J]．振动、测试与诊断，2005，25(1)：40-42．

Li Fang， Li Yourong，Wang Zhigang．De-noising technology based on Morlet wavelet transform and maximum likelihood estimation[J]．Journal of Vibration，Measurement & Diagnosis，2005，25(1)：40-42．(in Chinese)

[8] 侯新国，刘开培，魏建华．最佳小波包基改进软阈值的消噪方法及应用[J]．振动、测试与诊断，2008，28(4)：366-368.

Hou Xinguo，Liu Kaipei，Wei Jianhua．Application of improved soft threshold noise eliminating method based on optimal wavelet packet[J]．Journal of Vibration，Measurement & Diagnosis，2008，28(4)：366-368．(in Chinere)

[9] 张弦，王宏力．进化小波消噪方法及其在滚动轴承故障诊断中的应用[J]．机械工程学报，2010，46(15)：76-81．

Zhang Xian，Wang Hongli．Evolutionary wavelet denoising and its application to ball bearing fault diagnosis[J]．Journal of Mechanical Engineering，2010，46(15)：76-81．(in Chinese)

[10]刘文艺，汤宝平，蒋永华．一种自适应小波消噪方法[J]．振动、测试与诊断，2011，31 (1)：74-77．

Liu Wenyi，Tang Baoping，Jiang Yonghua．Research on an adaptive wavelet denoising method[J]．Journal of

Vibration，Measurement & Diagnosis，2011，31(1)：74-77．(in Chinese)

[11]Xia Xianggen．System identification using chirp signals and time-variant filters in the joint time-frequency domain [J]． IEEE Transactions on Signal Processing，1997，45(8)：2072-2084．

[12]Christiansen O．Time-frequency analysis and its applications in denoising [D]．Norway：University of Bergen，2002．

[13]Nezamoddin N K，Paul F．A gabor based technique for image denoising[C]∥Canadian Conference on Electrical and Computer Engineering．Waterloo，Canada：IEEE，2005：980-983．

[14]Walker J S，Chen Y J．Denoising gabor transforms [EB/OL]．[2013-12-05] http:∥www. uwec. edu./walkerjs．

[15]Parolai S．De-noising of seismograms using the S-transform [J]．Bulletin of the Seismological Society of America，2009，99(1)：226-234．

[16]Hon T K，Subramaniam S R，Georgakis A，et al． STFT-based Denoising of Elastogram[C]∥Proceedings of IEEE International Conference on Acoustics，Speech and Signal Processing(ICASSP)． Prague，Czech Republic：IEEE，2011：677-680．

[17]张贤达，保铮．非平稳信号分析与处理[M]．北京：国防工业出版社，1998：20-25．

[18]郭远晶，魏燕定，周晓军，等.基于S变换谱阈值去噪的冲击特征提取方法[J].振动与冲击，2014，33(21)：44-50.

Guo Yuanjing, Wei Yanding, Zhou Xiaojun, et al. An impact feature extracting method based on S-transformation spectrum threshold denoising[S].Journal of Vibration and Shock,2014,33(21):44-50.(in Chinese)

[19]IEEE PHM Society.2009 PHM challenge competition data set[EB/OL]．http:∥[2013-12-05].www.phmsociety.org/references/datasets．