圆柱滚子轴承振动信号时频特征提取及状态识别

刘湘楠，赵学智，何宽芳

（1.华南理工大学机械与汽车工程学院，广东广州 510641；2.佛山科技学院机电工程与自动化学院，广东佛山 528225）

引 言

随着旋转机械向智能化、精密化方向发展，旋转机械中某一部件产生故障可能导致较大的安全事故和经济损失。滚动轴承具有摩擦系数小、传递效率高和径向承载能力大等优点被广泛应用于旋转机械设备。据相关统计，旋转机械失效的案例中，45%～55%是由于滚动轴承发生故障而导致的［1］。实际工程应用中，滚动轴承通常在变转速、变负载和强冲击等恶劣工况下进行作业，磨损、疲劳、过载等因素都可能导致滚动轴承产生局部损伤，进而影响旋转机械设备运行的安全性和可靠性［2］。滚动轴承故障测试实验中，传感器采集的振动信号中含有大量的轴承运行状态信息，对振动信号进行分析可有效实现滚动轴承故障诊断［3］。传统的基于振动信号分析的滚动轴承故障诊断方法主要针对恒定工况，对变工况下的滚动轴承进行故障诊断的案例较少［4］。变工况下的振动信号更为复杂，包含了更多的运行状态信息［5］。因此，开展变工况下滚动轴承故障诊断方法研究，对旋转机械可靠性及安全性具有十分重要的意义。

由于旋转机械复杂的结构，导致传感器采集的滚动轴承振动信号通常具有非平稳特性［6］。在特征提取方面，傅里叶变换法作为一种全局变换方法，缺乏时间和频率“定位”功能，只适用于分析平稳信号［7］。时频分析法能够提供信号在时域和频域的联合分布信息，适用于分析非平稳信号［8］。常用的时频分析方法主要有：短时傅里叶变换（Short Time Fournier Transform，STFT）、连续小波变换（Con⁃tinue Wavelet Transform，CWT）、S 变换和广义S变换等［9⁃12］。STFT 作为一种加窗傅里叶变换方法，窗函数的长度决定了信号的时间分辨率和频率分辨率［9］。另外，由于受Heisenberg 不确定准则的限制，采用STFT 对信号进行时频分析时，难以同时获得良好的时域分布和频域分布［13］。CWT 作为一种窗函数可变的时频分析方法，克服了STFT 的不足，具有多分辨率分析特性［10］。但是小波基函数的选取缺乏自适应性，采用CWT 对信号进行分析时，其时频谱图存在频率干扰、能量泄漏和边界失真等不足［14］。S 变换是在STFT 和CWT 基础上发展形成的，既含有STFT 的相位信息，同时又保留了CWT 的多尺度分辨能力，并且其逆变换完全无损［15］。S 变换窗函数尺度与频率成反比，即在低频处频率分辨率较高，高频处时间分辨率较高。然而，由于S 变换中采用的小波函数是固定的，采用S 变换对信号进行分析时，存在时频分辨率调节力度不够，能量聚集性不足等缺点［16］。为克服S 变换的不足，Pinnegar 等［12］在S 变换的基础上引入窗函数调节系数，提出了GST 方法。GST 作为一种新的时频分析方法，在地震信号［17］、焊接裂纹声发射信号［18］中得到了广泛的应用，部分学者也将GST 应用于机械故障诊断领域［19⁃20］。但是采用GST 获得的时频矩阵存在维数过高、冗余信息过多等不足［21⁃22］。直接将GST 获得的时频矩阵作为智能诊断算法的输入，多的特征量会训练出更为复杂的模型，从而影响状态识别精度。

奇异值分解（Singular Value Decomposition，SVD）具有良好的稳定性，适用于提取时频矩阵中的特征信息［23］。TIAN 等［24］提出一种基于LMD⁃SVD 的变工况下滚动轴承故障特征信息的方法。郭凤仪等［25］采用S 变换对回路电流进行时频域变换，采用SVD 对时频矩阵进行分解，实现了不同电流条件下的串联型故障电弧特征信息提取。将GST 与SVD 相结合用于提取变工况下滚动轴承振动信号特征信息研究的成果尚未可见。

基于以上分析，本文将GST 和SVD 方法相结合，提出一种适用于提取变工况下滚动轴承故障振动信号特征信息的方法。搭建了某型特种车辆变速箱圆柱滚子轴承实验台架，并在实验台架上采集不同输入转速作用下圆柱滚子轴承典型状态件振动信号。在此基础上，利用GST⁃SVD 方法提取表征圆柱滚子轴承典型状态件特征信息的奇异值向量组，然后采用支持向量机对圆柱滚子轴承典型状态件进行分类识别。实验分析结果表明：该方法可有效实现变工况下圆柱滚子轴承振动信号特征提取及状态识别，为旋转机械在线监测提供一种有效手段。

1 理论基础

1.1 广义S 变换

对于一维连续信号x（t），其S 变换定义为［26］：

式中f表示频率，τ表示时移因子，g（t）表示高斯窗函数，σ表示尺度因子，且

由式（2）可知，高斯窗既是时间也是频率的函数，这就使得窗函数在低频处能够提供较高的频率分辨率，在高频处，提供较高的时间分辨率。

由于S 变换中高斯窗的尺度因子被限定为频率的导数，导致在实际工程应用中，S 变换存在时频分辨率调节力度不够，能量聚集性不足等缺点［16］。因此，为了克服S 变换的不足，Pinnegar 等［12］在高斯窗的尺度因子中引入调节参数p，改造S 变换的高斯窗函数，进一步加快或减慢时窗宽度随信号频率变化的速度，更好地适应具体信号的分析。

令尺度因子为：

结合方程（1）～（3）可得GST 的数学表达式为：

对于离散信号x（n），其GST 可表示为：

式中T表示采样时间间隔，N表示采样总点数，X（n/NT）表示x（n）的离散时间序列。

结合式（5）和（6）可知：原始信号经GST 时频域变换后，获得的时频矩阵是一个二维复数矩阵，行向量表示不同的频率值，列向量对应不同的时间点，矩阵元素表示信号的幅值和相位角。信号x（t）的GST 变换结果可表示为［25］：

式中Am×n表示幅值矩阵，θm×n表示相位矩阵。

另外，由方程（3）可知：GST 高斯窗的宽度可以通过选取不同的参数p值进行调整，从而改善时频分辨率，提高时频聚焦性能。当p=1 时，GST 等价于S 变换；当p>1时，高斯窗的宽度随着频率的增加而减小，不利于描述低频信号特征；当p<0 时，高斯窗的宽度随着频率的增加而增加，不利于描述高频信号特征。

参数p的取值决定了GST 的时频集聚性能的好坏，为实现参数p的选取，文献［27］提出一种时频分布聚焦性度量准则，根据该准则可有效选取高斯窗口函数最优调节参数p值，定义时频分布聚焦性度量准则为：

式（8）的离散形式可表示为：

评价的标准为：M（p）取最小值表示时频分布聚焦性最好，此时所对应的p值即为高斯窗口函数的最优调节参数。

本文选取参数p值的步骤如下：

（1）对于任意的p∈(0，1 ]，根据方程（4）对信号x（t）进行GST；

（2）对GST 系数进行能量归一化：

（3）根据方程（8）计算GST 的时频聚焦性，本文取q=2；

（4）通过求取M（p）的最小值选取窗口函数最优调节参数popt：

1.2 SVD 基本理论

对于矩阵H∈Rm×n，其SVD 定义为［23］：

式中U和VT分别为m×m和n×n阶矩阵，S为m×n阶对角矩阵，即：

式中 0 表示零矩阵，k=min（m，n），σ1≥σ2≥…≥σk>0，σ1≥σ2≥…≥σk>0 称为矩阵Hm×n的奇异值。

基于以上分析，式（12）可表示为：

式中ui和vi分别表示矩阵U和V中的列向量，且ui∈Rm×1，vi∈Rn×1。

定义一个包含ui和vi的子矩阵Hi=uivTi，因此，式（14）可表示为：

根据式（15）可知，采用SVD 提取矩阵中的有效信息的关键在于选取有效奇异值。本文采用文献［28］所提出的方法选取有效奇异值个数r，定义有效奇异值个数r为：奇异值由迅速衰减向趋于平缓的转折点所对应的序数。

1.3 SVM 模型

SVM 主要通过寻找一个超平面来对样本进行分割，在基于间隔最大化的分割原则下，其最终问题归类于在约束条件下求解二次规划问题［29］。

设样本集{xi，yi}，i=1，2，…，N，其中xi为第i个样本，yi∈{-1，1}，那么最优超平面方程定义为：

式中ω表示该平面的法向量，b为常数项。

根据SVM 的定义，将上述最优超平面选取问题转换为二次规划问题：

通过Langrange 乘子法解决上述规划问题，那么原问题的对偶问题为：

式中ai表示Langrange 乘子，Q（a）的最大值取决于训练集{xiT，xj}和{yi，yj}。

假定a*i为最优Langrange 乘子，x表示测试数据，那么最优超平面函数（fx）定义为：

1.4 本文方法

本文结合GST，SVD 和SVM，提出一种变工况下圆柱滚子轴承振动信号特征提取及状态识别方法，其基本流程如图1所示。

图1 本文方法诊断流程Fig.1 The proposed method diagnostic process

该方法的具体步骤如下：

（1）搭建某型特种车辆变速箱圆柱滚子轴承实验台架，在实验台架上采集不同转速工况下的滚动轴承不同状态件振动信号，并从中随机选取多个样本构成训练集和测试集；

（2）采用GST 分别对训练集中的单个样本信号进行时频域转换，将得到时频矩阵作为特征矩阵；

（3）对特征矩阵进行SVD 分析，选取有效奇异值构建特征向量集；

（4）将得到的特征集输入SVM 模型中对其进行训练；对测试集中的样本信号重复上述步骤（2）和（3）进行特征参数提取，然后输入到已训练好的SVM 模型中实现故障状态分类识别。

2 时频分析方法仿真分析

构造了非平稳仿真信号模拟滚动轴承局部损伤引起的冲击信号，分别采用STFT，CWT，S 变换和GST 方法对仿真信号进行分析，以验证GST 良好的时频聚集性。仿真信号的数学表达式如下：

式中t1=mod（t，1/33），故障引起的冲击信号特征频率为33 Hz，mod 表示求余函数，f1=3500 Hz，f2=450 Hz，f3=150 Hz，f4=48 Hz，r（t）表示均值为0，方差为0.16 的高斯白噪声。

仿真信号采样频率为8 kHz，采样点数为4000个，图2为仿真信号时域波形及其频谱图。

由图2可知，由于噪声信号的干扰，传统的频谱分析难以获取故障特征频率，缺乏“定位功能”，仅从信号时域波形图中无法获取仿真信号中的冲击特征信息。

图2 仿真信号波形图Fig.2 Waveform of the simulated signal

采用STFT，CWT，S 变换（p=1）和GST（popt=0.76）等时频方法分别对仿真信号进行分析，其中STFT 采用Hamming 窗，CWT 采用尺度为64 的Morlet 小波。图3为仿真信号经不同时频分析方法处理后获得的时频谱图。

图3 不同时频分析方法处理后获得的时频谱图Fig.3 Time-frequency spectrum obtained by different timefrequency analysis methods

如图3（a）所示，采用STFT 对仿真信号进行分析，能够提取低频成分（f2，f3，f4），但由于时间分辨率的限制，难以有效提取高频成分（f1）和冲击信号频率成分（f0）。如图3（b）所示，采用CWT 对仿真信号进行分析，能够提出提取低频成分（f2，f3，f4）、高频成分（f1）和冲击信号频率成分（f0），但高频成分和冲击信号频率成分分辨率较低，存在能量泄露。如图3（c）所示，采用S 变换对仿真信号进行分析，能够提出提取低频成分（f2，f3，f4）、高频成分（f1）和冲击信号频率成分（f0），但由于其窗函数尺度与频率成反比，因此，低频成分具有较高的频率分辨率和较低的时间分辨率，高频成分和冲击信号具有较低的频率分辨率和较高的时间分辨率。如图3（d）所示，采用GST 对仿真信号进行分析，能够有效提取低频成分（f2，f3，f4），且在高频成分（f1）附近检测到明显的冲击特征，冲击时间间隔约为0.03 s，与冲击信号频率成分（f0=33 Hz）相对应。

基于以上分析可知，GST 相较于传统的时频分析方法在冲击信号时频分辨率方面具有显著的优越性，且GST 时频矩阵中包含了冲击信号特征信息。

3 案例分析

3.1 实验方案及数据获取

如图4所示，搭建了某型特种车辆变速箱圆柱滚子轴承实验台架。在台架上分别采集了某特种车辆变速箱中N218 圆柱滚子轴承外圈磨损、滚动体故障和轴承正常三种典型状态件的振动信号。表1为圆柱滚子轴承N218 基本参数。

图4 台架模拟故障试验台Fig.4 The bench simulation fault test bench

表1 N218 圆柱滚子轴承基本参数Tab.1 Basic parameters of the cylindrical roller bearings

在轴承故障件制作过程中，外圈磨损状态是在轴承外圈外表面一圈磨掉0.15 mm 而成；滚动体故障状态是对轴承的1 个滚动体磨损0.15 mm 而成。图5为圆柱滚子轴承故障件。

图5 圆柱滚子轴承故障件Fig.5 Faulty parts of the cylindrical roller bearing

试验过程中，利用加速度传感器采集圆柱滚子轴承三种典型状态件在不同转速下的振动信号。

本文的研究目的在于分析GST⁃SVD 方法应用于提取变工况下圆柱滚子轴承故障特征信息的可行性，因此，仅对传感器A1 采集的振动信号进行分析。传感器A1 垂直贴在变速箱输出轴端箱盖表面。

在变速箱处于三挡时进行信号采集，三挡齿轮传动比为1.9286。采样频率为20 kHz，采样点数为20000 个，负载为100 HP，输入转速分别设定为500，800，1000 和1200 r/min。

表2为本文选取的实验数据集，A，B，C 和D 分别表示四种不同工况，每种工况随机选取150个样本（包含三种典型状态件各50 个样本），共600 个样本。其中训练集包含360个样本、测试集包含240个样本。

表2 实验数据集Tab.2 Experimental data set

3.2 圆柱滚子轴承典型状态件时频特征分析

以变速箱处于三挡，输入转速为500 r/min，负载为100 HP 时，传感器A1 检测到的振动信号为例，对GST 应用于提取故障特征信息的有效性进行分析。图6为传感器A1 检测到的圆柱滚子轴承三种典型状态件振动信号波形图。分别采用STFT，CWT，S 变换和GST 对传感器A1 检测到的实验数据进行时频分析，分析结果如图7所示。

图6 圆柱滚子轴承典型状态件振动信号波形Fig.6 Vibration signal waveform of cylindrical roller bearing in typical state

根据文献［28］中所提及的滚动轴承故障特征频率计算公式计算，可得外圈故障频率为31.32 Hz，滚动体故障频率为14.34 Hz。如图7（a）所示，采用STFT 对圆柱滚子轴承典型状态件进行时频分析，由于时间分辨率的限制难以有效识别冲击信号特征。如图7（b）所示，采用CWT 对圆柱滚子轴承典型状态件进行时频分析，冲击信号频率成分分辨率较低，存在能量泄露。如图7（c）所示，采用S 变换对圆柱滚子轴承典型状态件进行时频分析，能够有效识别电机输入转频（8.333 Hz）、外圈故障频率（31.0559 Hz）及滚动体故障频率（14.34 Hz），但频率分辨率较低。由图7（d）可知，正常轴承的时频谱图中，在频率轴约2237 Hz 处存在周期性冲击特性，其周期约为0.1056 s，对应的频率约为9.4697 Hz，刚好与电机输入转频（8.333 Hz）相近，正常轴承时频谱中不存在故障特征信息；外圈故障的时频谱图中，在频率轴约2237 Hz 处存在周期性冲击特性，其周期约为0.0322 s，对应的频率约为31.0559 Hz，刚好与外圈故障频率（31.32 Hz）相近；滚动体故障的时频谱图中，在频率轴约2237 Hz 处存在周期性冲击特性，其周期约为0.1053 和0.0657 s，对应的频率约为9.4967 和15.2207 Hz，刚好与电机输入转频（8.333 Hz）和滚动体故障频率（14.34 Hz）相近。

图7 不同时频分析方法处理后获得的时频谱图Fig.7 Time-frequency spectrum obtained by different time-frequency analysis methods

基于以上分析可知，GST 相较于传统的时频分析方法在冲击信号时频分辨率方面具有显著的优越性，且GST 时频矩阵中包含了冲击信号特征信息。

3.3 基于GST‑SVD 的圆柱滚子轴承振动信号特征提取

采用GST 分别对实验数据集600 个样本进行时频域转换，将获得的时频矩阵作为特征矩阵，对特征矩阵进行SVD 分解，获得特征矩阵的奇异值向量。图8为实验数据集600 个样本的奇异值曲线图。

由图8可知，变工况下圆柱滚子轴承振动信号经GST⁃SVD 处理后，获得的奇异值向量具有较好的一致性和稳定性，且每个特征矩阵所包含的特征信息可以由28 个奇异值向量表征。当轴承发生故障时，特征矩阵的奇异值数值大小发生了明显的变化，因此，由多组奇异值向量组成的矩阵包含了原始振动信号中的故障特征信息，可作为变工况下圆柱滚子轴承不同状态件状态类型识别的特征参数。

图8 变工况下圆柱滚子轴承振动信号经GST-SVD 后获得的奇异值曲线Fig.8 Singular value curve obtained by GST-SVD for signal of cylindrical roller bearing under variable operating conditions

3.4 圆柱滚子轴承振动信号状态识别

将GST⁃SVD 所提炼出的特征参数输入SVM。本文采用网格寻优法对SVM 中的核函数参数c和惩罚函数g进行筛选，得到的最优的h和p分别为1.4 和2.2。图9为SVM 模型对变工况下圆柱滚子轴承不同状态件的分类识别结果。标签1～3 分别表示轴承正常、外圈故障和滚动体故障等轴承三种典型状态。

由图9可知，采用本文方法对不同转速工况下圆柱滚子轴承进行诊断时，训练集识别误差为0，测试集识别误差为2.083%。因此，采用本文方法可以有效实现变工况下圆柱滚子轴承故障诊断。

图9 本文方法识别结果Fig.9 Identification results by the proposed method

3.5 方法对比

为了验证本文提出的变工况下圆柱滚子轴承振动信号特征提取方法的优越性和有效性，分别采用以下方法对变工况下圆柱滚子轴承进行故障特征提取。

（1）S 变换和SVD：首先采用S 变换将一维时域信号转变成二维时频矩阵，选取时频矩阵作为特征矩阵；采用SVD 对特征矩阵进行分解，获取表征不同故障类型的特征向量。

（2）局部均值分解（Local Mode Decomposi⁃tion，LMD）和SVD：采用LMD 分解原始信号，获得若干生产函数（Product Functions，PFs）分量；采用SVD 对PFs 分量构成的矩阵进行分解，获取表征不同故障类型的特征向量［24］。

（3）小波包近似熵：选用db5 小波基函数对原始信号进行四层小波包分解，通过计算各层频带的近似熵得到特征向量［30］。

将上述特征提取方法提取的特征向量分别作为SVM 输入，采用SVM 对变工况下圆柱滚子轴承三种典型状态件进行分类识别，表3为不同特征提取方法的SVM 状态识别结果。

表3 不同转速工况下圆柱滚子轴承故障诊断模型的分类结果Tab.3 Classification results of fault diagnosis model for cylindrical roller bearings under different working conditions

由表3可知，本文所提出的方法能有效地提高圆柱滚子轴承故障状态识别的效率和精度。

4 结 论

（1）构造了非平稳仿真信号模拟滚动轴承局部损伤引起的冲击信号，采用STFT，CWT，S 变化和GST 分别对仿真信号进行时频域转换，验证了GST相较于其他时频分析方法具有更高的时频分辨率。

（2）结合GST 和SVD 提出一种变工况下滚动轴承振动信号特征提取方法。搭建了某型特种车辆变速箱圆柱滚子轴承实验台架，在实验台架上采集了不同转速下圆柱滚子轴承三种典型状态件振动信号。采用本文方法对实验获得的振动信号进行分析，结果表明：该方法获得的奇异值向量具有较好的一致性和稳定性，由多组奇异值向量组成的矩阵包含了原始振动信号中的故障特征信息，可作为变工况下圆柱滚子轴承不同状态件模式识别的特征参数。

（3）采用SVM 对不同特征提取方法提取的特征信息进行分类识别，结果表明：本文方法具有较高的识别精度和效率。