李心一 谢志江 罗久飞

1.重庆大学机械工程学院，重庆，400044 2.重庆邮电大学先进制造工程学院，重庆，400065

0 引言

滚动轴承是旋转机构的重要部件之一，其运行状态的好坏直接影响旋转机构的工作状态和使用寿命，因此，滚动轴承故障诊断具有重要的工程现实意义。

在滚动轴承故障诊断中，由损伤点产生的低频振动称为轴承的故障特征频率，采用频谱分析法进行轴承诊断时，通过查看轴承振动信号频谱中是否包含故障特征频率，可以判断轴承故障类型。研究表明，在滚动轴承故障诊断中，由故障产生的振动信号处于强噪声背景，且往往出现不同程度的调制现象，因此，准确地提取故障特征频率是判断轴承故障类型的关键。Hilbert变换是一种常用的信号解调方法，可用于轴承信号的解调，但有一定局限［1］，文献［2-3］提出了基于迭代Hilbert变换的滚动轴承故障诊断方法。Teager能量算子是在研究非线性语音建模时用于分析和跟踪窄带信号的能量所引入的算法，具有比Hilbert变换更好的解调效果，非常适合检测信号中的冲击成分，并广泛应用于轴承故障诊断中。文献［4］利用Teager能量算子提取轴承故障特征频率，判断轴承故障类型。文献［5-6］将经验模式分解与Tea⁃ger能量算子解调方法相结合，用Teager能量算子对单分量的固有模态函数进行解调后分析，从而可有效判断轴承和齿轮的故障类型。文献［7］利用Teager能量算子增强倒阶次谱提取轴承微弱故障特征。在强噪声背景下，为了能有效提取滚动轴承故障特征频率，文献［8-9］将最小熵解卷积（minimum entropy deconvolution，MED）应用于滚动轴承和齿轮箱的故障诊断，获得了良好的效果。文献［10］对MED在滚动轴承故障特征提取中的应用进行了详细的分析。文献［11］将Teager能量算子与MED相结合用于滚动轴承的故障诊断中，首先利用MED对振动信号进行降噪预处理，然后通过Teager能量算子求得包络谱，可有效提取轴承故障特征。

加窗插值快速傅里叶变换（fast Fourier trans⁃form，FFT）方法最早用于提高FFT的计算精度［12］，后广泛应用于电力谐波分析中。文献［13-14］对基于最大旁瓣衰减窗（maximum sidelobe decay window，MSDW）插值FFT做了大量的研究，此类方法较为简单，易于实现。

为了提高滚动轴承故障频率的识别精度，本文针对轴承故障信号的特点，并结合MED和Tea⁃ger能量算子解调的优点，采用D.Belega加窗插值方法，利用常用的MSDW对解调信号进行加权处理并进行频谱校正。

1 最小熵解卷积

当滚动轴承发生故障时，采集到的信号可表示为

由式（1）可以看出，故障轴承的振动信号 y(n)是故障轴承脉冲冲击信号u(n)、背景噪声n(n)与传递路径h(n)的卷积结果。当轴承发生故障时，由于受随机噪声、传递路径等的影响，输出信号y(n)往往不能真实地反映输入信号u(n)的特性。由文献［15］可知，熵最小时，信号的峭度最大，所以MED核心思想是利用熵最小原理构造L阶逆滤波器 fi(n)，使输出信号 y(n)经此逆滤波还原为输入信号u(n)，从而突出冲击脉冲，即

MED求解方法主要有目标函数法和特征向量法，具体计算步骤见文献［16］。

由以上分析可知，通过选择合理的滤波器参数，采用MED对信号进行处理，可降低背景噪声的干扰，增强信号冲击成分，非常适合轴承振动信号的滤波预处理。

2 Teager能量算子解调

频率为 f0的单频信号以采样频率 fs进行采样，得到离散时间信号

式中，A为幅值；φ为初相位。

则信号 x(n)的Teager能量算子ψ定义为［17］

由式（2）可以看出，Teager能量算子的输出与信号的频率和幅值的平方成正比，因此Teager能量算子的输出对瞬态冲击成分更敏感。与Hilbert变换相比，Teager能量算子计算简单，只需计算相邻3个采样点，因此具有良好的时间分辨率。

轴承工作时由于局部损伤（剥落、腐蚀、裂纹、擦伤等），振动信号具有周期性冲击的特点，而Teager能量算子具有增强信号瞬态特征的优点，因此非常适合检测信号瞬态成分。

在实际工程应用中Hilbert变换和Teager能量算子是2种常用的解调方法，文献［18］对这2种方法进行了详细的对比。与Hilbert变换相比，Tea⁃ger能量算子具有运算量小和解调精度高的优点，因此本文选用Teager能量算子用于轴承振动信号的解调分析。

3 解调谱频谱校正

最大旁瓣衰减窗是具有良好旁瓣性能的余弦组合窗，其表达式如下：

其中，I为窗函数的项数；ai为窗函数系数并满足约束条件：

可计算得到2～5项最大旁瓣衰减窗函数系数见表1。由表1可知，Hanning窗是2项最大旁瓣衰减窗，适用于最大旁瓣衰减窗插值FFT方法。图1所示为I取2、3、4、5，窗长和采样频率均为64时，窗函数的幅频响应曲线，其中横坐标以频率点B表示，B等于采样频率除以采样点。由图1可知，随着项数的增加旁瓣电平变小，旁瓣衰落速率增大。但是，MSDW的主瓣宽度随着项数的增加而变宽，从而降低了频率分辨率。

表1 最大旁瓣衰减窗函数系数Tab.1 Coefficients of MSDW

图1 2～5项MSDW窗幅频特性Fig.1 Amplitude-frequency properties of 2～5 term MSDW

文献［13-14］对基于最大旁瓣衰减窗插值FFT法的频率估计进行了详细阐述，与其他经典窗FFT插值法相比，此类算法具有解析解，并且无需计算一元超越方程，插值公式较为简单也更易实现。不失一般性，考虑式（3）所示的单频信号，用长度为N的Hanning窗w(n)对该信号进行加权处理，得到加窗信号xw(n)=x(n)w(n)。根据卷积定理，其加窗离散傅里叶变换为

其中，W(·)为Hanning窗的离散时间傅里叶变换，当 N＞＞1时，其表达式为

若只考虑正频部分，则Xw(k)可表示为

然而，因非同步采样等原因，在实际情况中k0往往不是整数，因此，信号的精确频率 f0可表示为

式中，l为峰值频点k0附近抽样幅值最大所对应的谱线号，为正整数；q为频率偏差，其取值范围为-0.5≤0＜0.5。

可见，得到精确频率 f0的关键是求出频率偏差值q。有多种方法可以求解q值，文献［13］介绍了一种十分简单的三谱线插值方法，该方法利用峰值频点k0附近的3根离散频谱的幅值求解q。令l-1为l左边谱线，l+1为l右边谱线，记 y1=|X(l-1)|，y2=|X(l)|,y3=|X(l+1)|，另记：

定义比值

根据式（8）和式（9），b可改写为

把式（7）代入式（13），通过计算得

根据式（14），可求得

求得q后，根据式（9）可求出信号的精确频率 f0。

4 轴承故障诊断方法

根据前文分析，基于汉宁窗插值FFT法的轴承故障诊断方法包括以下几个主要步骤：①利用MED方法对振动信号进行降噪预处理；②计算轴承振动信号Teager能量算子输出；③对Teager能量算子的输出进行快速傅里叶变换，得到Teager解调谱；④对解调谱加Hanning窗，然后确定l、l+1、l-1；⑤根据三谱线幅值插值计算求解频率偏值q；⑥得到校正后的轴承故障特征频率值，判断故障类型。详细流程见图2。

图2 滚动轴承诊断流程图Fig.2 Flow chart of rolling bearing diagnosis

5 轴承故障诊断实例

实例分析数据来自美国凯斯西储大学轴承研究中心轴承故障模拟试验台［19］。内圈和外圈试验采用的轴承为安装于电机驱动端的6205-2RS JEM SKF深沟球滚动轴承，其轴承节径dp=39.039 8 mm，滚动体直径dB=7.940 04 mm，滚动体个数Z=9，接触角 β=0°。电机负载为0～2.2 kW，转速约为1720～1797r/min，采样频率为12kHz。

5.1 轴承内外圈故障诊断

轴承内圈故障点直径0.177 8 mm，故障深度0.279 4 mm，电机转速1 797 r/min。计算轴承内圈故障的通过频率（the ball passing frequency inner race，BPFI）为162.185 973 Hz。选取数据段 2 049～3 072共1 024个采样点进行分析，图3所示为内圈故障信号。从图3中可以明显看出，轴承内圈故障信号具有周期性冲击特性。

图3 内圈故障信号图Fig.3 Time domain waveform of inner race fault signal

图4 所示为图3信号的Teager解调谱。由图4可知，Teager能量算子仅能提取出内圈故障特征频率fBPFI=164.062 5 Hz（并且与理论计算值有较大误差），其倍频并不明显。对图3信号进行MED降噪预处理，参考文献［8-10，20］选择合适的MED参数（滤波器阶数L=20，最大循环迭代次数M=30，迭代误差e=0.01）。图5所示为经MED滤波后的内圈信号解调谱。与图3相比，内圈信号经MED预处理后，增强了信号中的冲击成分，突出了轴承内圈故障特征频率及其倍频，但是其精度与理论计算值仍有一定误差。

图4 内圈信号解调谱Fig.4 Demodulation spectrum of inner race fault signal

图5 MED滤波后内圈信号解调谱Fig.5 Demodulation spectrum of inner race fault signal filtered with MED

为了提高轴承故障频特征频率识别精度，计算经MED处理后的轴承内圈信号的Teager能量算子的输出，并对该瞬时Teager能量序列进行傅里叶变换，得到解调谱。然后利用加Hanning窗插值FFT方法对解调频进行频谱校正，得到校正后的轴承故障特征频率见图6。由图6可知，经过MED和频谱校正后，轴承内圈故障特征频率及其倍频的精度都得到了有效提高，与理论计算值相吻合。

图6 内圈故障特征频率校正值Fig.6 Corrected BPFI of inner race

为了验证本文所提方法的普遍适用性，排除人为选择因素和数据样本等因素的影响，分别选取轴承内圈故障点直径为0.177 8 mm的随机数据样本100个和50个，并分别对应每个样本数据长度为1 024和2 048。引入均方根误差（root mean square error，RMSE）指标评估算法有效性，反映校正结果的精密度。RMSE定义如下：

式（16）中，N代表样本数，fci是每个样本的校正故障特征频率值，fj表示计算的理论频率值。表2、表3所示分别为100个样本和50个样本在故障点直径为0.177 8 mm、电机转速为1 721～1 797 r/min的RMSE仿真结果。

表2 100个样本的RMSE（内圈 0.177 8 mm）Tab.2 RMSE of 100 samples（IR 0.177 8 mm）

表3 50个样本的RMSE（内圈0.177 8 mm）Tab.3 RMSE of 50 samples（IR 0.177 8 mm）

对比表2、表3结果可以看出，未校正的故障频率RMSE最大值为5.419 789（100个样本）和2.982 938（50个样本），最小值也接近2，存在很大误差。总体上看，在不同电机运行状态下，经频谱校正后的故障特征频率RMSE明显优于直接由Teager能量算子解调的结果。

轴承外圈选择故障点直径为0.177 8 mm，故障深度为0.279 4 mm，故障点位于6点钟方向。同外圈分析方法一样，表4、表5所示分别为100个样本和50个样本在故障点直径为0.177 8 mm、电机转速为1 725～1 796 r/min的RMSE仿真结果。

表4 100个样本的RMSE（外圈0.177 8mm）Tab.4 RMSE of 100 samples（OR 0.177 8 mm）

表5 50个样本的RMSE（外圈0.177 8 mm）Tab.5 RMSE of 50 samples（OR 0.177 8 mm）

由表4、表5可知，与其他电机转速相比，在电机转速为1 773 r/min时，未校正值和校正值比较接近，这是因为此状态下未校正特征频率值为105.468 75 Hz，与理论计算特征频率值105.930 11 Hz已十分接近。在此状态下，校正后的RMSE值未有明显减小。对于大多数情况，经频谱校正后的RMSE值效果都较好。

5.2 窗函数影响

窗函数的选择对加窗插值FFT的精度有较大影响。分析采用加窗插值FFT进行轴承故障诊断时窗函数对轴承特征频率精度的影响，图7、图8所示为分别选用2～5项MSDW对轴承内外圈进行加窗插值FFT方法的RMSE值曲线（数据长度均为1 024），可以看出，RMSE值随着窗函数项数的增加而增大。由图1可知，2-MSDW具有最窄的主瓣宽度，主瓣宽度决定了加窗序列的频率分辨率，因此较窄的主瓣宽度能带来较高的频率分辨率，从而提高轴承特征频率的识别精度。

图7 2～5项MSDW的RMSE（内圈 100样本）Fig.7 RMSE of 2～5 term MSDW（IR 100 samples）

图8 2～5项MSDW的RMSE（外圈 100样本）Fig.8 RMSE of 2～5 term MSDW（OR 100 samples）

6 结论

（1）根据Teager能量算子和故障滚动轴承振动信号的特点，计算振动信号的Teager能量算子输出，并对其进行傅里叶变换后得到解调谱，可以有效识别轴承的故障特征频率，但在大多数情况下识别精度不高，与理论计算值存在偏差。

（2）本文介绍的三谱线插值FFT方法计算精度高且易实现，所用的最大旁瓣衰减窗为信号处理最常用的一类窗函数，算法易于编程，非常适合嵌入式系统的应用。

（3）主瓣较窄的窗函数具有较高的频率分辨率，对故障特征频率校正有一定影响。

（4）经轴承内外圈诊断实例分析，将MED、Teager能量算子解调和加窗插值FFT法相结合，在大多数情况下选取较少的分析点也能显著提高轴承特征频率的识别精度，性能上优于传统的Teager能量算子解调方法。

［1］ 郑晓溪，丁康，江利旗.希尔伯特变换解调分析在故障诊断中应用的局限性研究［J］.汕头大学学报（自然科学版），1999，14（2）：40-46.ZHENG Xiaoxi，DING Kang，JIANG Liqi.Limitaton Studies on Application of Hilbert Transform Demodula⁃tion Analysis in Fault Diagnosis［J］.Journal of Shan⁃tou University（Natural Science），1999，14（2）：40-46.

［2］ 秦毅，任兵，李宁，等.改进的迭代希尔伯特变换及滚动轴承故障诊断应用［J］.振动与冲击，2013，32（11）：83-88.QIN Yi，REN Bing，LI Ning，et al.An Improved Itera⁃tive Hilbert Transformation and Its Application in Roll⁃ing Bearing Fault Diagnosis［J］.Journal of Vibration and Shock，2013，32（11）：83-88.

［3］ 邓四二，王勇，王恒迪.基于IHT的共振解调技术的滚动轴承故障诊断方法［J］.航空动力学报，2012，27（1）：69-74.DENG Sier，WANG Yong，WANG Hengdi.Resonance Demodulation Technique Based on IHT for Rolling Bearings Fault Diagnosis［J］.Journal of Aerospace Power，2012，27（1）：69-74.

［4］ 王天金，冯志鹏，郝如江，等.基于Teager能量算子的滚动轴承故障诊断研究［J］.振动与冲击，2012，31（2）：1-5.WANG Tianjin，FENG Zhipeng，HAO Rujiang，et al.Fault Diagnosis of Rolling Element Bearings Based on Teager Energy Operator［J］.Journal of Vibration and Shock，2012，31（2）：1-5.

［5］ 李辉，郑海起，杨绍普.基于EMD和Teager能量算子的轴承故障诊断研究［J］.振动与冲击，2008，27（10）：15-17.LI Hui，ZHENG Haiqi，YANG Shaopu.Bearing Fault Diagnosis Based on EMD and Teager Kaiser Energy Operator［J］.Journal of Vibration and Shock，2008，27（10）：15-17.

［6］ 赵军，崔颖，赖欣欢，等.随机共振降噪下的齿轮微弱故障特征提取［J］.中国机械工程，2014，25（4）：539-546.ZHAO Jun，CUI Ying，LAI Xinhuan，et al.Weak Fea⁃ture Extraction of Gear Faults Based on Stochastic Res⁃onance Denoising［J］.China Mechanical Engineering，2014，25（4）：539-546.

［7］ 杨青乐，梅检民，肖静，等.Teager能量算子增强倒阶次谱提取轴承微弱故障特征［J］.振动与冲击，2015，34（6）：1-5.YANG Qingle，MEI Jianmin，XIAO Jing，et al.Weak Fault Feature Extraction for Bearings Based on an Or⁃der Cepstrum Enhanced with Teager Energy Operator［J］.Journal of Vibration and Shock，2015，34（6）：1-5.

［8］ ENDO H，RANDALL R B.Enhancement of Autore⁃gressive Model Based Gear Tooth Fault Detection Technique by the Use of Minimum Entropy Deconvolu⁃tion Filter［J］.Mechanical Systems&Signal Process⁃ing，2007，21：906-919.

［9］ SAWALHI N，RANDALL R B H，ENDO H.The En⁃hancement of Fault Detection and Diagnosis in Rolling Element Bearings Using Minimum Entropy Deconvolu⁃tion Combined with Spectral Kurtosis［J］.Mechanical Systems&Signal Processing，2007，21：2616-2633.

［10］ 龚廷恺，袁晓辉，王细洋.最小熵反褶积的数学形态法在滚动轴承故障特征提取中的应用［J］.中国机械工程，2016，27（18）：2467-2471.GONG Tingkai，YUAN Xiaohui，WANG Xiyang.Ap⁃plications of Mathematical Morphology Method to Fault Feature Extraction of Rolling Bearings Based on Minimum Entropy Deconvolution［J］.China Me⁃chanical Engineering，2016，27（18）：2467-2471.

［11］ 陈海周，王家序，汤宝平，等.基于最小熵解卷积和Teager能量算子直升机滚动轴承复合故障诊断研究［J］.振动与冲击，2017，36（9）：45-50.CHEN Haizhou，WANG Jiaxu，TANG Baoping，et al.Helicopter Rolling Bearing Hybrid Faults Diagnosis Using Minimum Entropy Deconvolution and Teager Energy Operator［J］.Journal of Vibration and Shock，2017，36（9）：45-50.

［12］ JAIN V K，COLLINS W L，DAVIS D C.High Accura⁃cy Analog Measurements via Interpolated FFT［J］.IEEE Transactions on Instrumentation&Measure⁃ment，2007，28（2）：113-122.

［13］ BELEGA D，DALLET D.Frequency Estimation via Weighted Multipoint Interpolated DFT［J］.IET Sci.Meas.Technol.，2008，2（1）：1-8.

［14］ BELEGA D，DALLET D.Multifrequency Signal Anal⁃ysis by Interpolated DFT Method with Maximum Side⁃lobe Decay Windows［J］.Measurement，2009，42（3）：420-426.

［15］ WIGGINS R A.Minimum Entropy Deconvolution［J］.GeophysicalProspectingforPetrole，1980，16（1）：21-35.

［16］ 王宏超，陈进，董广明.基于最小熵解卷积与稀疏分解的滚动轴承微弱故障特征提取［J］.机械工程学报，2013，49（1）：88-94.WANG Hongchao，CHEN Jin，DONG Guangming.Fault Diagnosis Method for Rolling Bearing’s Weak Fault Based on Minimum Entropy Deconvolution and Sparse Decomposition［J］.Journal of Mechanical En⁃gineering，2013，49（1）：88-94.

［17］ KAISER J F.On a Simple Algorithm to Calculate the‘Energy’of a Signal［C］//Proceedings of IEEE Inter⁃national Conference on Acoustics，Speech，Signal Pro⁃cessing.Albuquerque，1990：381-384.

［18］ POTAMIANOS A，MARAGOS P.A Comparison of the Energy Operator and Hilbert Transform Approach⁃es for Signal and Speech Demodulation［J］.Signal Processing，1994，37（1）：95-120.

［19］ Case Western Reserve University Bearing Data Cen⁃ter.Bearing Data Center Fault Test Data［EB/OL］.［2012-11-10］.http：//www.Eecs.Case.edu/laborato⁃ry/bearing.

［20］ JIANG R L，CHEN J，DONG G M，et al.The Weak Fault Diagnosis and Condition Monitoring of Rolling Element Bearing Using Minimum Entropy Deconvolu⁃tion and Envelop Spectrum［J］.Engineering Science Engineers，Part C：Journal of Mechanical Engineer⁃ing Science，2013，227（5）：1116-1129.