基于阶频谱相干的变转速滚动轴承内外圈复合故障特征分离提取

杨新敏，郭 瑜，华健翔

(昆明理工大学 机电工程学院,昆明 650500)

滚动轴承是旋转机械的核心部件之一，轴承故障诊断在机械设备的监测与维护中具有重要意义。相比单一故障，复合故障相互耦合，增加了诊断难度，复合故障的诊断近来备受关注。研究表明，内、外圈故障占轴承故障的90%[1]，因此，本文的研究基于内、外圈复合故障展开。目前，复合故障特征提取方法主要有EMD(empirical mode decomposition)[2]、VMD(variational mode decomposition)[3]、CYCBD(cyclostationarity blind deconvolution)[4]和MCKD(maximum correlated kurtosis deconvolution)[5]等，但上述方法多用于转速较稳定的工况，而在实际中变转速工况也较为常见，比如：设备的启动和停止、机器人臂的重复性动作等。在变转速工况下，故障引起的冲击在时域上不再是周期的，常规的诊断方法不再适用。变转速工况下滚动轴承复合故障诊断方法多基于深度学习[6]，传统特征提取方法的研究较少。因此，滚动轴承变转速工况下的复合故障诊断研究具有普遍的研究价值。

在变转速旋转机械设备的故障诊断中，计算阶比跟踪(computed order tracking,COT)为经典方法[7]，其本质是将时域信号等角度重采样为角域信号，从而消除转速波动的影响。但是COT需要获得准确的参考轴键相信号对信号进行等角度重采样，其限制了该方法的应用。与COT原理不同，近年来，Abboud等[8]在谱相关(spectral correlation,SC)基础上提出循环非平稳信号的谱分析方法，应用阶频谱相干(order-frequency spectral coherence,OFSCoh)算法实现变转速工况下轴承信号的故障特征提取且该方法无需对原信号进行等角度重采样，有利于避免重采样带来的误差。国内，夏均忠等[9]提出基于集成阶频谱相关的方法应用于变转速工况下轴承的故障特征提取，在特定频带内对OFSCoh函数积分，将OFSCoh的表达由三维变为二维，有利于突显故障阶次，更利于故障识别。但上述研究在降维过程中积分频带通过观察法确定，不具有自适应性且当轴承内外圈复合故障信息相互耦合时，积分频带难以确定。

另一方面，共振解调是轴承故障特征提取的有效方法，其核心在于准确选择合适的解调频带，滚动轴承发生复合故障时，不同故障对应冲击往往在不同共振频带激起共振，且具有不同的能量[10]。分别选择包含不同故障信息丰富的共振区作为解调频带对信号进行分离，是复合故障特征提取的有效方法[11]。OFSCoh函数是以双变量的形式表示载波频率f和循环阶次α(或故障特征阶次)的阶频谱相关性，其值越大，表明故障信息越集中在以f为中心的一定带宽的频带内，因此，其可作为确定不同故障解调频带中心频率的指标。本文结合OFSCoh的优势，提出一种基于OFSCoh解调频带确定的复合故障特征分离提取方法。首先，对信号进行OFSCoh计算；然后，以轴承内、外圈对应的故障阶次区间分别对OFSCoh函数进行积分获得特征频带谱，将该谱线最大值对应的频率确定为解调频带中心频率，以最大转频对应的3倍故障频率作为解调带宽；最后，分别对信号进行带通滤波，并计算其改进包络谱(improved envelope spectrum,IES)，从而实现轴承复合故障特征分离提取。

1 理论基础

变转速工况下，故障冲击在时域上随着速度的变化而增加或衰减，此类信号称为角域/时间循环平稳 (angle/time cyclostationary signals,AT-CS)信号，假设X(t)是AT-CS信号，可通过傅里叶级数的形式表示该信号，即

(1)

(2)

式中：ω(t)为瞬时角速度，rad/s。将时变自相关函数的时间变量替换成角度/时间域，同时考虑角度和时间变量，角度/时间自相关函数(angle/time correlation function,ATCF)表示为

R2X(τ,θ)=E{X(t(θ))X(t(θ)-τ)*}

(3)

式中：E为总体平均；τ为延时变量。对ATCF进行两次傅里叶变换(第一次将时域映射到频域，第二次实现角域到阶次域的映射)即可得到阶频谱相关函数(OFSC)

S2X(α,f)=Fτ→fθ→α{R2X(τ,θ)}

(4)

式中：α为循环阶次(在轴承振动中，对应安装轴每转一圈的冲击次数)；f为谱频率(对应载波频率)。由式(4)可知，OFSC是关于f和α的双变量函数。为抵消信号功率在频域的不均匀分布，将OFSC进行归一化处理，定义OFSCoh函数为

(5)

由于OFSCoh结果是以f和α的双变量的三维谱图形式呈现，通过其获取故障信息需要一定的信号处理方面的技能，为了增强识别故障的可读性，将f进行积分，从而将三维图降为二维图，积分后的阶频谱定义为IES，定义式为

(6)

式中，F1和F2分别为f的积分下限和上限，关于OFSCoh和IES的更多细节可参阅Abboud等的研究。

2 基于OFSCoh的故障特征提取

本文基于OFSCoh对变转速工况下轴承复合故障进行分离提取。在故障轴承振动信号的OFSCoh函数中，故障阶次(循环阶次α)集中分布在以某f为中心的区域，函数值越大，表明感兴趣成分(故障阶次)越集中在该值对应的f处。根据这一特性，结合已知的理论故障阶次沿f轴进行切片，可得到f与OFSCoh值的关系谱线，选择值最大时对应的f作为包含故障信息丰富的解调频带中心频率。但由于滚动轴承会产生随机滑移，导致理论故障频率与实际故障频率可能存在1%～2%的误差[12]。为避免理论故障阶次与实际故障阶次存在误差导致感兴趣成分与f映射关系无意义的问题。本文通过对OFSCoh函数积分的形式获得感兴趣的成分与f的关系谱线。根据轴承故障阶次选择特定的区间对OFSCoh函数积分。轴承外圈、内圈理论故障阶次计算公式为

(7)

(8)

式中：Oo和Oi分别为外圈和内圈故障阶次；n为滚动体个数；d为滚动体直径；D为节圆直径；φ为接触角。本文考虑理论故障阶次与实际故障阶次存在 2%的误差，以理论故障阶次为区间中值，区间下限为0.98Ox，上限为1.02Ox(下标x表示o和i，分别代表外圈和内圈对应的故障阶次)。因此，积分区间可由式(9)计算

αx=[0.98Qx,1.02Qx]

(9)

式中，αx为内、外圈故障阶次对应的积分区间。

根据积分区间对OFSCoh函数积分，获得f和OFSCoh值的关系谱线，称其为特征频带谱，选取特征频带谱最大值对应的f作为解调频带中心频率，fnx可由式(10)计算

(10)

本文所提方法实现流程如图1所示。主要实现步骤包括：

步骤1对振动信号进行OFSCoh分析，获得关于f和α的双变量图(OFSCoh分析的体计算细节参阅Abboud等的研究)。

步骤2将轴承参数分别代入式(7)和式(8)计算出外圈和内圈的理论故障阶次，根据式(9)分别计算内、外圈特征频带谱的积分区间，然后对OFSCoh函数积分求出内、外圈对应的特征频带谱，根据式(10)分别求出内、外圈对应的解调频带中心频率fnx，根据转速计算解调带宽Δfx，从而获得滤波参数组{fnx,Δfx}。

步骤3根据步骤2计算出的滤波参数组分别对信号进行带通滤波，然后根据式(6)求出滤波后信号的IES，从而实现变转速工况下轴承内、外圈复合故障特征分离提取。

图1 复合故障特征分离提取流程Fig.1 Compound fault feature separation and extraction process

3 仿真信号分析

仿真轴承内、外圈复合故障在变速工况下振动信号对本文方法进行验证。传感器拾取的信号具有多个能量集中的共振区域的特征，不同区域包含着不同大小的故障能量。为了让仿真信号更接近实际，在仿真中，采用多共振带具有不同能量的内、外圈成分的复合故障仿真模型，仿真模型如式(11)所示。

(11)

式中：x1(t)为同一共振区外圈幅值较大、内圈幅值较小的仿真信号；x2(t) 为同一共振区外圈幅值较小、内圈幅值较大的仿真信号；fr为转频；n(t)为高斯白噪声；fno为外圈故障信息丰富的共振频率；fni为内圈故障信息丰富的共振频率；Si为外圈第i次冲击的幅值，由于是变转速，幅值大小受转频影响，其值与转频满足关系：Si=c×fr；Sj为内圈第j次冲击的幅值，其值与转频满足关系：Sj=d×fr；Ti和Tj分别为第i次和第j次冲击发生的时刻；β为衰减系数；τ为随机滑移的时间。仿真参数如表1所示。

表1 仿真参数Tab.1 Simulation parameters

仿真中设轴承外圈阶次Oo为3.5，即每转发生3.5次冲击；内圈故障阶次Oi为6.0，仿真时长为2 s。设转频方程为fr(t)=8t+5；在设定时间内转过的转数可由式(12)计算

(12)

式中，θi(i=1,2,3,…,I，I为设定时间内所转圈数)根据式可求出每转对应的时间T。因此，外圈、内圈每次冲击对应的时间分别表示为Ti=T/Oo，Tj=T/Oi。综上分析，根据仿真模型式和表1参数得到轴承内、外圈变转速(升速)工况下复合故障仿真信号的时域波形如图2所示。图2(a)、图2(b)分别表示式(11)中的x1(t)和x2(t)；图2(c)表示x(t)，即：加入高斯白噪声后信噪比为-10 dB的轴承复合故障仿真信号的时域波形。

图2 变转速复合故障仿真信号Fig.2 Simulation signal of variable speed compound fault

首先采用Guo等提出的基于快速谱峭度(fast-kurtogram)的包络阶比跟踪方法与本文方法进行对比，验证本文方法的有效性。对仿真信号进行快速谱峭度处理，结果如图3所示。

图3 仿真信号的快速谱峭度 Fig.3 Fast-kurtogram of simulation signal

图3中确定的解调频带为{2 203 Hz,94 Hz}，但谱峭度值次大的解调区难以辨识。根据上述确定的解调频带提取仿真信号滤波后的包络并对其进行等角度重采样(重采样频率为1 024)，然后求其包络阶比谱如图4所示。图4中无法识别内圈和外圈的故障阶次，说明受噪声干扰，快速谱峭度算法解调频带的选取受限。

图4 仿真信号的包络阶比谱Fig.4 Order envelope spectrum of simulation signal

为进一步证明本文算法分离信号的有效性，采用VMD与本文方法对比，本文采用文献[14]所提的观察中心频率法确定模态数为6，分别求出各IMF分量的包络阶次谱，结果如图5所示。从图5中可以看出，IFM2的包络阶次谱中能较易识别出内圈故障阶次，而IMF3的包络阶次谱中能观察到外圈故障阶次，但高次谐波受噪声干扰严重，不易识别。

应用本文的基础方法对仿真信号进行处理。首先，求信号的OFSCoh和IES，结果如图6所示。图6(a)OFSCoh图中白色矩形虚线框理论上应该显示1阶外圈故障阶次(左虚线框)和1阶内圈故障阶次(右虚线框)，但从图中难以识别故障阶次。对信号进一步处理，求仿真信号IES，结果如图6(b)，图中内、外圈故障阶次相互混淆，受噪声干扰严重，不利于识别。

用本文所提方法对信号进行分离，在OFSCoh分析的基础上，根据仿真的理论故障阶次，通过式(9)分别计算各故障对应的积分区间，外、内圈对应的积分区间分别为：[3.43，3.57]，[5.88，6.12](如图6(a)白色虚线框标注区域)。

图5 IMFs的包络阶比谱Fig.5 Envelope spectrum of IMFs

图6 仿真信号的OFSCoh分析结果Fig.6 Order-frequency spectral coherence analysis results of the simulation signal

根据式(10)求得特征频带谱如图7所示。图7中内圈参数(故障阶次区间)对应的特征频带谱的最大值，即内圈故障信息丰富的解调频带中心频率fni为1 500 Hz，仿真的最大转频为26 Hz，可计算得解调带宽为468 Hz，带通滤波参数组为{1 500 Hz,468 Hz}；同理，外圈参数对应的带通滤波参数组即为{4 125 Hz,273 Hz}。需要注意的是，若存在不同故障激起的共振频率重叠的情况，此时，优先将特征频带谱值较大的元件对应的f作为其故障信息丰富的解调频带中心频率，而另一元件对应的解调频带应在其他频段内搜索，以达到信号分离提取的目的，具体实例见下文试验分析部分。

图7 特征频带谱Fig.7 Characteristic frequency band spectrum

根据所得的滤波参数组分别对信号进行带通滤波，然后分别求分离后信号的IES，结果如图8所示。对比图4、图5和图8可知，受噪声影响，快速谱峭度对解调频带的识别受限，无法提取故障。VMD分解的敏感分量IMF2的包络阶比谱中能识别内圈故障阶次及其倍频，但外圈故障阶次难以在其他分量中识别。而图8中能较好识别外圈故障阶次和内圈故障阶次及其倍频。仿真试验说明本文方法有利于识别包含各故障信息丰富的解调频带，对复合故障特征分离提取有效。

图8 本文所提方法对信号分离后的IESFig.8 IES after signal separation by the method of this paper

4 试验信号分析

为了进一步验证所提方法的有效性，在QPZZ-Π试验平台上设置外圈、内圈复合故障对本文算法进行验证，试验台如图9(a)所示。以N205EM型号轴承作为研究对象，轴承参数(节圆直径D/mm：38.5；滚动体直径d/mm：7.94；滚动体个数n：12；接触角φ：0°)。采用线切割技术在轴承外圈切一个宽约0.5 mm,深约0.2 mm的槽模拟轴承外圈故障，在内圈处切一个宽约0.5 mm，深约1.5 mm的槽模拟轴承内圈故障，故障如图9(b)所示。

图9 QPZZ-Π试验平台与轴承复合故障Fig.9 QPZZ-Π test platform and compound fault of bearing

试验中以电机自由升速和停机模拟变转速工况，采样频率fs为51.2 kHz。试验中转频曲线和振动数据时域波形如图10所示。对试验信号进行快速谱峭度处理，结果如图11(a)所示。根据图11(a)可知，确定的解调频带为{1 867 Hz,533 Hz}，但谱峭度值次大的解调频带难以辨识。根据上述确定的解调频带提取信号滤波后的包络并对其进行等角度重采样(重采样频率为1 024)，最后求包络谱，其包络阶比谱如图11(b)所示。

图10 试验信号转频曲线与时域波形Fig.10 Experimental signal frequency curve and time domain waveform

用VMD对信号进行处理，选择模态数为4，分别求出各IMF分量的包络阶次谱如图12所示。从图12中可看出，IMF2的包络阶次谱可以较好识别出内圈故障阶次，在IMF4的包络阶次谱中能观察到外圈故障阶次，但受噪声干扰严重，不利于识别。

应用本文基础方法对试验信号进行处理，首先，求信号的OFSCoh和IES，结果如图13所示。

图11 快速谱峭度方法处理结果Fig.11 The results of fast-kurtogram

图12 IMFs的包络阶比谱Fig.12 Envelope spectrum of IMFs

图13 试验信号的OFSCoh分析结果Fig.13 OFSCoh analysis results of the experimental signal

图13(a)中白色矩形虚线框理论上应该显示1阶外圈故障阶次(左虚线框)和1阶内圈故障阶次(右虚线框)，但从图中难以识别出故障阶次；图13 (b) 为信号的IES，图中内圈故障阶次相对明显，基本可以识别，但外圈故障阶次受噪声干扰严重，难以识别。

用本文所提方法对信号进行分离，将轴承参数代入式(7)和式(8)计算轴承外圈、内圈的故障阶次，分别为Oo=4.76和Oi=7.24。通过式(9)分别计算各故障对应的积分区间，外、内圈对应的积分区间分别为：[4.66，4.86]，[7.10，7.38](如图13(a)白色矩形虚线框标注区域)，根据式(10)求得特征频带谱，如图14所示。

图14 特征频带谱Fig.14 Characteristic frequency band spectrum

图14内圈参数的特征频带谱最大值对应的f为2 000 Hz，即:内圈故障信息丰富的解调频带中心频率fni为2 000 Hz，试验中，最大转频为15 Hz，计算出带宽为326 Hz，因此，带通滤波参数组为{2 000 Hz,326 Hz}。需要注意的是，外圈参数的特征频带谱最大值对应的频率也为2 000 Hz，与内圈对应的解调频带中心频率重叠，但在此频率下外圈对应幅值比内圈小，说明该f下，内圈故障信息丰富，应作为内圈的解调频带。因此，外圈对应的解调频带应该在其他频带内搜索，取外圈对应的特征频带谱值的第二峰值点作为其中心频率，即：fno为19 200 Hz，带宽为214 Hz，带通滤波参数组即为{19 200 Hz,214 Hz}。根据上述所得外圈、内圈对应的解调频带，分别对信号进行带通滤波，然后求滤波后信号的IES，结果如图15所示。

图15 本文所提方法对信号分离后的IESFig.15 IES after signal separation by the method of this paper

对比图11、图12和图15，三种方法都能较易识别内圈故障阶次，但较弱的外圈故障受干扰严重，快速谱峭度无法对包含其故障信息丰富的解调频带进行识别，无法提取外圈故障。VMD分解的敏感分量IMF4能观察到外圈故障阶次，但受干扰严重难以识别。从图15(a)中可看出，外圈故障阶次及其倍频明显，容易识别。结果表明，本文所提方法能有效识别外圈、内圈故障信息丰富的解调频带，实现复合故障特征分离提取。

5 结 论

本文根据滚动轴承内、外圈对应的理论故障阶次区间对OFSCoh积分，获得特征频带谱，根据其最大值对应的载波频率，确定其为故障信息丰富的解调频带中心频率，根据转频计算带宽，从而求出内、外圈故障对应的解调频带，分别对信号带通滤波，求滤波后信号的IES，从而实现复合故障特征分离提取。本文方法无需对信号进行等角度重采样，可避免等角度重采样带来的计算误差。采用积分的形式获得特征频带谱，有效避免滑移造成理论故障阶次与实际故障阶次之间的差异导致感兴趣成分与载波频率映射关系无意义的问题。本文方法解调频带选择准确，算法鲁棒性好，通过故障分离的方式可有效增强弱故障信息。仿真和试验表明，与快速谱峭度和VMD方法对比，本文所提方法能更好实现变转速工况下轴承内外圈复合故障特征分离提取，具有一定的应用价值。