多周期微分均值和循环谱相干的轴承故障特征提取方法

1.本发明涉及一种多周期微分均值和循环谱相干的轴承故障特征提取方法，属于故障诊断技术与信号处理分析技术领域。


背景技术：

2.轴承作为旋转机械的支撑部件，其健康程度直接影响旋转机械的精度和寿命。轴承发生故障时，滚动体与滚道在故障位置处的接触刚度发生变化，对应的瞬时角速度(ias)将发生规律性变化。因此，基于ias信号的轴承故障特征提取是故障诊断领域的热点之一。
3.轴承作为旋转机械的支撑部件，其并不传递扭矩，但在径向载荷的作用下，滚动体与滚道在故障位置处的接触刚度将发生规律性变化，通过单位角度的时间发生规律性变化，对应的ias信号产生规律性波动。然而轴承故障在低径向载荷工况引起的ias变化较齿轮故障来说更微弱。此外，由于滚动轴承存在1％～2％的随机滑动，即滚动轴承相邻角度间隔不尽相同。因此，轴承故障对应分量具有二阶循环平稳特性。
4.另一方面，由于编码器内径与轴装配误差，编码器安装误差在工程应用中无法避免；编码器安装误差的能量幅值与转速成正相关，即编码器安装误差对轴承信号的调制作用随着转速上升而显著提升，即对轴承故障分量的干扰程度是不同的。因此，编码器安装误差的消除对于有效揭示轴承故障特征至关重要。


技术实现要素：

5.在编码器信号中，轴承早期故障特征较弱以及编码器安装误差的调制作用随着转速的上升而增大，进而导致轴承特征可能无法有效辨识。为解决该问题，本发明提供了一种多周期微分均值和循环谱相干的轴承故障特征增强方法，该方法基于微分技术不受幅值干扰的优势和多周期的累积特性，采用多周期微分均值技术对瞬时角速度iasi信号中轴承故障特征分量进行增强，进而抑制编码器安装误差、估计误差和测量噪声分量的干扰；其次，针对多周期微分均值技术引起的频率模糊问题，结合滚动轴承故障分量具有的二阶循环平稳特性，采用循环谱相干技术消除多周期微分均值技术带来的频率模糊现象并进一步提取轴承故障特征分量。
6.本发明基于多周期微分均值和循环谱相干的轴承故障特征提取方法具体如下：
7.步骤1：获得包含轴承故障信息的瞬时角速度信号；
8.通过picoscope信号采集系统获取光学编码器中包含轴承故障信息的瞬时角位移和对应的时间，采用向前差分法计算获得瞬时角速度iasi信号，其计算式如下
[0009][0010]
式中，iasi表示第i时刻的瞬时角速度，i＝1,2,3,
…
，δφ＝2π/n，δti＝t
i+1-ti；n表示编码器每圈的光栅数。
[0011]
步骤2：采用多周期微分均值(multi-period differential means，mpdm)算法增
强瞬时角速度iasi信号中轴承故障特征分量；
[0012]
基于轴承故障时会引起iasi信号的规律性变化及其对应iasi分量较无故障状态具有突变趋势，采用具有对缓变信号值较小对突变信号值较大特点(例如编码器安装误差，测量误差等)的微分算法，并考虑以轴承理论故障角度周期作为加窗的角度间隔，结合多周期累计特性，提出采用多周期微分均值(multi-period differential means，mpdm)技术以增强滚动轴承故障特征，其计算式为：
[0013][0014]
式中，轴承最大随机滑动角度n
sm
＝nr/49/f
reb
，r为编码器安装轴和故障轴承转轴间的传动比，nw表示微分窗长，k为感兴趣轴承故障冲击的周期数，在数据点为m和冲击周期数为k时的微分位置q＝m+(k-1)n/f
reb
，在冲击周期k时轴承随机滑动角度范围[h,h]＝[q-2n
sm
,q+2n
sm
]，在数据点为m和冲击周期数为k时的微分范围[q,q]＝[m+(k-1)n/f
reb
,m+(k-1)n/f
reb
+nw]；k＝1,2,...,k；f
reb
为轴承故障理论特征阶次，q＝m+(k-1)n/f
reb
+nw；m＝2n
sm
,2n
sm
+1,...,m；处理后数据长度m＝length(iasi)-kn/f
reb-n
w-2n
sm
，n/f
reb
表示轴承故障理论角度间隔，length(iasi)表示原始iasi信号的长度。
[0015]
公式(2)中需要设置四个参数：k、r、f
reb
和nw。k是提升mpdm算法鲁棒性的主要参数，通常k≥3，k值越大，mpdm对感兴趣分量能量幅值的增强效果越好，但需要更高的计算代价；根据传动系统参数和轴承尺寸参数，r和f
reb
可以被计算获得；为有效包含轴承故障微分分量，根据经验2floor(n
sm
)《nw《10floor(n
sm
)可满足基于iasi信号的轴承故障特征提取，floor(
·
)为向下取整操作。
[0016]
步骤3：为解决mpdm算法引入的“频率模糊”问题，基于轴承故障分量的二阶循环平稳特性，采用循环谱相干技术消除多周期微分均值技术带来的频率模糊现象并进一步提取轴承故障分量，对于被mpdm算法增强后信号mpdm(m)，其循环谱相关计算式为：
[0017][0018]
式中，fft[
·
]表示信号mpdm(m)在循环周期t内的傅里叶变换，[
·
]
*
表示函数共轭，e{
·
}表示数学期望，α表示循环频率，f表示谱频率，为抑制背景噪声分布不均匀的影响，使用循环谱相干(cyclic spectral coherence，cscoh)进行归一化处理，其表达式为：
[0019][0020]
基于公式(4)，沿着谱频率f轴进行积分，获得以一维谱函数以实现轴承故障冲击特征的有效辨识，其表达式为：
[0021][0022]
式中f1和f2分别表示各子频带的上下截止频率，|cscoh(α,f)|表示绝对值操作。
[0023]
本发明的有益效果是：
[0024]
(1)本发明基于微分技术不受较大幅值干扰的优势，结合多周期累计特性，提出一种基于iasi信号的mpdm轴承故障特征增强方法，其可有效增强轴承故障冲击特征；
[0025]
(2)本发明采用cscoh技术不仅可有效解决mpdm算法导致“频率模糊”的缺陷，而且基于mpdm算法增强的轴承故障分量二阶循环平稳特性，cscoh算法可进一步有效提取滚动轴承故障分量；
[0026]
(3)本发明提出一种基于iasi信号mpdm和cscoh相结合的轴承故障特征增强方法。
附图说明
[0027]
图1为本发明方法流程示意图；
[0028]
图2为实施例1中仿真信号数据和对应的包络阶次谱，其中图(a)为ias波形，图(b)是图(a)的包络谱；
[0029]
图3为实施例1中cscoh算法和所提方法的分析结果，其中图(a)为直接使用cscoh算法对原始信号分析获得的结果，图(b)为mpdm算法的滤波信号，图(c)为图(b)的包络阶次谱，图(d)为采用cscoh算法在图(b)上进一步处理的结果；
[0030]
图4为实施例1中不同理论获得的结果，其中图(a)为drs的滤波信号，图(b)为图(a)的包络阶次谱，图(c)为cpw的滤波信号，图(d)为图(c)的包络阶次谱；
[0031]
图5为实施例2中轴承测试实验台；
[0032]
图6为实施例2中获得的轴承外圈故障数据和对应的包络阶次谱，其中图(a)为ias波形，图(b)是图(a)的包络谱；
[0033]
图7为实施例2中cscoh算法和所提方法的分析结果，其中图(a)为直接使用cscoh算法获得的结果，图(b)为mpdm算法的滤波信号，图(c)为图(b)的包络阶次谱，图(d)为采用cscoh算法在图(b)上进一步处理的全频带包络阶次谱；
[0034]
图8为实施例2中不同理论获得的结果，其中图(a)为drs的滤波信号，图(b)为图(a)的包络阶次谱，图(c)为cpw的滤波信号，图(d)为图(c)的包络阶次谱。
具体实施方式
[0035]
下面结合本发明实施例中的附图，对本文发明实例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员没有做出创造性劳动前提下所获得的所有其他实施例，都属于本文发明保护的范围，本实施例中如无特殊说明的方法，均为常规方法。
[0036]
实施例1：如图1所示，本实施例是本发明方法用于仿真轴承故障特征的提取，具体过程包括：
[0037]
仿真分析采用的轴承外圈故障信号模型计算式为：
[0038][0039][0040]
式中，w(θ)表示平均角速度，wo(θ)表示编码器安装误差，ξ表示阻尼系数，fn表示固频，ψ＝θ-jθ-τj，角度序列θ＝2π/n,4π/n,6π/n...，n(θ)表示编码器测量噪声，ρ＝δr/r是几何中心和旋转中心的偏心距离的比值，r为编码器孔直径，δr是几何中心和旋转中心的偏心距离，β表示编码器旋轴与转轴之间的倾斜角度，初始角度θe∈[φe,2kπ+φe]，θ
t
∈
[φ
t
,2kπ+φ
t
]，c表示轴承故障冲击幅值。在本仿真信号中，n＝5000，c＝0.001，f
reb
＝3.56
×
，θ＝n/f
reb
＝1404，《τj》＝0，max{τj}＝28，fn＝50
×
，
×
表示特征阶次，ξ＝0.03，w(θ)＝5rad/s，ρ＝0.001，β＝0.03，信噪比snr＝-15，转频分量fr＝1
×
。
[0041]
步骤1：仿真信号iasi波形和对应的包络阶次谱如图2(a)和图2(b)，信号iasi的长度length(iasi)为105；可见轴承故障特征谱线无法有效辨识。因此，需要对轴承故障分量进行增强处理。
[0042]
步骤2：对原始iasi中轴承故障分量进行增强
[0043]
2-1、设置参数：k＝3，f
reb
＝3.56
×
，r＝1，nw＝56，n
sm
＝nr/49/f
reb
＝28；
[0044]
2-2、采用mpdm算法对轴承故障分量进行增强处理，其计算式为：
[0045][0046]
式中，m＝length(iasi)-kn/f
reb-n
w-2n
sm
＝10
5-3
×
5000/3.56-56-56＝95674，m＝2n
sm
,2n
sm
+1,
…
；m＝56,57,
…
,95674；q＝m+(k-1)n/f
reb
，k＝1,2,3，[h,h]＝[q-2n
sm
,q+2n
sm
]，q＝m+(k-1)n/f
reb
+nw。mpdm算法增强处理的结果如图3(b)和图3(c)所示。
[0047]
步骤3：根据采样率n＝5000，循环频率α＝600，窗宽bw＝1024，应用循环谱相关技术获得谱频率f和循环频率α的双变量谱图，其计算式如下：
[0048][0049]
为抑制噪声分布的不均性，对cscoh(α,f)进行归一化操作，获得谱频率f和循环频率α的双变量谱，循环谱相干计算式为：
[0050][0051]
步骤4：沿着谱频率f，范围为[0,n/2]，获得的全频带包络谱如图3(d)所示。可见，轴承故障特征谱线可被有效辨识，编码器安装误差和测量噪声得到有效抑制。
[0052]
为进一步显示本文所提方法的有效性，直接采用cscoh对原始iasi进行分析，其结果如图3(a)所示，可见轴承故障特征谱线无法有效辨识；进一步地，采用drs和cpw算法分别对原始信号进行处理，获得的滤波信号如图4(a)和4(c)所示，对应的包络阶次谱如图4(b)和4(d)，其中drs算法中延迟为100，窗长为53。可见，转频分量被有效抑制，但轴承故障特征谱线依然无法有效辨识。其原因在于轴承故障分量的能量幅值微弱，导致其谱线无法有效辨识。
[0053]
实施例2：本实施例描述的是本发明方法用于实际滚动轴承外圈故障特征的提取
[0054]
本实施例中采用轴承测试试验台，如图5所示，在实验台上安装一个etf100-h851007b光学编码器，编码器线数n＝5000，并采用106采样率的picoscope高速采集装置获取对应的角度信息和时间信息。本试验台的轴承类型是nu206e(nb＝13,eb＝9.525,e
p
＝46,α＝0)，为模拟轴承外圈故障，采用线切割方式在外圈上加工宽度大约0.3mm，深度约为0.28mm的槽；由以下计算式得到轴承外圈故障特征频率f
reb
为5.15
×
。
[0055]
[0056]
步骤1：获得的iasi轴承外圈故障波形如图6(a)所示，信号iasi的长度length(iasi)为477700，首先直接对原始iasi信号进行包络阶次谱分析，如图6(b)所示。可见转频阶次谱线占主导地位，而轴承故障阶次谱线基本淹没于干扰分量中。因此，编码器安装误差的干扰需要被抑制以增强轴承故障特征。
[0057]
步骤2：对原始iasi中轴承故障特征进行增强处理
[0058]
2-1、设置参数：k＝3，f
reb
＝5.15
×
，r＝1，nw＝50，n
sm
＝nr/49/f
reb
＝20。
[0059]
2-2、采用mpdm算法对轴承故障分量进行增强处理，其计算式为：
[0060][0061]
式中，m＝length(iasi)-kn/f
reb-n
w-2n
sm
＝477700-3
×
5000/5.15-50-40＝474697；m＝2n
sm
,2n
sm
+1,
…
；m＝40,41,
…
,474697；q＝m+(k-1)n/f
reb
；k＝1,2,3；[h,h]＝[q-2n
sm
,q+2n
sm
]；q＝m+(k-1)n/f
reb
+nw；mpdm算法增强处理的结果如图7(b)和图7(c)所示。
[0062]
步骤3：根据采样率n＝5000，循环频率α＝256，窗宽bw＝30，应用谱相关技术获得谱频率f和循环频率α的双变量谱图，其计算式如下：
[0063][0064]
为抑制噪声分布的不均性，对cscoh(α,f)进行归一化操作，获得谱频率f和循环频率α的双变量谱，计算式为：
[0065][0066]
步骤4：沿着谱频率f，范围为[0,n/2]，获得的全频带包络谱如图7(d)所示。可见，轴承故障特征谱线可被有效辨识，编码器安装误差和测量噪声得到有效抑制。
[0067]
为进一步显示本文所提方法的有效性，直接采用cscoh对原始iasi进行分析，其结果如图7(a)所示，可见轴承故障特征谱线依然无法有效辨识；进一步地，采用drs技术对原始信号iasi进行分析，其中drs算法中延迟为200，窗长为78，如图8(a)和8(b)所示。可见转频特征谱线占优，然而轴承故障谱线淹没于背景噪声中。此外，如图8(c)和8(d)，采用cpw用于抑制具有一阶循环平稳特性的分量，但轴承故障特征谱线无法有效辨识。因此，在本实验数据下，cscoh、drs和cpw均未有效增强轴承故障分量。
[0068]
本文中应用了具体的实例对本发明的原理及实施方式进行了参数，以上实例的说明只是用于帮助理解本发明及核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变指出，综上，本说明书内容不应理解为本发明的限制。