CN107884188B - 一种基于模式转换的变转速/变负载滚动轴承故障诊断方法 - Google Patents

Abstract

本发明公开了一种基于模式转换的变转速/变负载滚动轴承故障诊断方法，其包括以下步骤：通过假设故障滚动轴承振动信号至少由两个具有不同分布特性的信号源组成，建立轴承振动信号模型，进而基于阈值技术及尾估计方法对信号进行初步分离分析以获得具有不同分布特性的信号模式，然后结合期望‑最大化理论建立信号的隐马尔可夫模型以获得噪声脉冲，并将原始信号与之相减从而进一步提高信号的信噪比，最后利用传统包络谱分析信号，评估滚动轴承是否有故障以及该故障的类型。本发明一方面克服了变负载/变转速条件下振动信号难以诊断的问题；另一方面，解决了当非高斯噪声存在的情况下，高能脉冲改变信号的包络性时导致传统包络法无法应用的问题。

Description

一种基于模式转换的变转速/变负载滚动轴承故障诊断方法

技术领域

本发明属于机械设备故障诊断技术领域，尤其是一种于模式转换的变转速/变负载滚动轴承故障诊断方法。

背景技术

现代社会对机械设备的需求日益增大，作为机械设备关键核心部件的滚动轴承，在电力、石化、冶金、机械、航空航天以及一些军事工业部门中被广泛使用, 同时也是最易损伤的部件之一。据有关资料统计,机械故障的70％是振动故障,而振动故障中有30％是由滚动轴承引起的。这是因为滚动轴承在机械设备中起着承受载荷和传递载荷的作用,而且工作条件比较恶劣。其运行状态对整台机器的安全运行影响极大，一旦滚动轴承发生故障或者失效，将对机器运行带来极大的安全隐患。因此，研究滚动轴承的故障诊断方法，对提高结构系统整体运行安全性和可靠性，防止由于损伤积累而引起重大事故的发生有着重要的现实意义。

对滚动轴承进行故障诊断的关键就是提取故障所产生的周期性冲击特征及频率信息。通常，除固定的外圈之外，滚动轴承中的运动部件(包括内圈、保持架和滚动体)局部损伤造成的振动激励源与传感器之间的位置相对变化，而且在轴和轴上多种零部件振动的影响作用下，信号中的干扰激励多，滚动轴承振动信号的成分比较复杂，因此，相对于外圈故障来说，内圈、保持架和滚动体故障特征微弱，不易提取。所以提高信号的信噪比对于获得特征信号尤为重要。

近几十年来，随着信号处理理论的不断发展，相关研究人员提出了多种降噪方法，然而针对变负载/变转速条件下的机械系统使得振动信号不稳定难以诊断的问题，一直缺乏适用的方法。

发明内容

为了克服现有技术的不足，本发明提供了一种基于模式转换的变转速/变负载滚动轴承故障诊断方法，本发明在处理变转速/变负载滚动轴承故障信号中具有良好的效果，适用于在机械系统中去噪及故障检测。

为了实现上述目的，本发明采用的技术方案是：一种基于模式转换的变速/变负载滚动轴承故障诊断方法，其特征在于，包括以下内容：首先假设提取的振动信号至少包含两个具有不同分布特性的信号源并据此建立振动信号模型；其次基于尾估计对信号进行初步分离分析，获得具有不同分布特性的信号模式；然后利用模式转换迭代参数以建立隐马尔可夫模型获得噪声脉冲，并将原始信号与之相减，进一步提高信号的信噪比；最后利用传统包络谱分析信号，评估滚动轴承是否有损伤以及该损伤的类型；

具体步骤为：

(1)信号初步分离

设定硬阈值，当固定邻域内信号X₁，X₂，...，X_n中的正负峰值超过给定阈值则被识别，并定义模式R₁为正峰，R₃为负峰，R₂为移除正峰和负峰后的信号模式，

计算R₂的高斯累计分布与经验累计分布曲线之间的距离并选取最小的距离作为新的硬阈值，并重新带入原信号再次分类得到新的R₁、R₂、R₃；距离采用 Kolmogorov-Smirnov统计公式计算：

上式中Φ(·)为标准高斯分布的累计分布函数，F_n(·)为经验分布函数，其定义如下：

上式中1_A为A集的指示器；

模式识别后使用不同的统计方法分别分析R₁、R₂、R₃并获得其分布状态；对 R₂使用Kolmogorov-Smirnov测试验证其服从高斯分布，如果P值超过置信程度 0.05，则样本的高斯性需验证；

R₁、R₃的分布特性分析基于其经验右尾的行为状态，对于样本 X₁，X₂，...，X_n，用公式1-F_n(t)计算其经验右尾，

上式中F_n(t)为经验累积分布函数，并将其与幂函数bt^-γ进行比较，

基于以上分析，可以得到信号R₁、R₃的分布状态服从帕累托分布，帕累托分布可用其概率密度函数定义，公式如下：

上式中γ为尺度参数，x_m为转移参数。对于帕累托分布，理论右尾用如下公式表示：

R₁、R₃可视为具有相似参数的帕累托分布；

由以上分析，可得到两种分布状态，分别为高斯分布和帕累托分布；

(2)模式转换

经过初步分离，确定两种模式分别为帕累托分布与高斯分布，分别对应着 R₁、R₂，用下式表达：

上式中R_i为第i个信号的状态变量，且所观察到信号的统计特性取决于实际状态R_i；sP(γ,x_m)表示帕累托分布，可用概率密度函数定义如下：

通过过渡矩阵P将信号从一个状态转换到另一个状态，其中包含在时刻i时的状态l过渡到时刻i+1时状态j的概率p_1j＝P(R_i+1＝j|R_i＝l)，结果获得分别来自状态R₁和状态R₂的两个概率，分别为P(R_i＝1|x₁,x₂,...,x_n)和P(R_i＝2|x₁,x₂,...,x_n)，分离方法正是基于这些概率。

模型参数的估计是基于期望-最大化算法，其包含二步，分别是期望值化 (E-Step)和最大化(M-Step)；每次迭代中，状态值先根据先前计算的参数进行推论；然后根据这些推论估计新的参数集；重复迭代，直到到达局部最大似然函数，该算法的具体过程如下：

初始化：设置初始参数矢量j的取值为1和 2，

期望值化(E-Step)：假设θ^(k)是先前迭代中最大化迭代法(M-Step)所累积的参数向量；具体包含以下步骤：

1)向前过滤：当i＝1,2，...，n时迭代方程如下：

上式中x_i＝(x₁,x₂,...,x_i)为信号矢量，gj(x_i；θ^(k))为当测量值来自状态j时时刻i下的信号概率密度函数(g₁和g₂分别为对称的帕累托和高斯密度)；当 P(R_n＝j|x_n；θ^(k))累计时有下式：

迭代的起点为

3)向后过滤：当i＝n-1,n-2，...，1时迭代方程如下：

最大化(M-Step)：新参数θ^(k+1)的估计可以通过最大似然函数推导。实际考虑以下：

由于x_m不存在最大似然估计，因此其可任意选择为常数，

最后，得到过渡概率根据以下公式计算：

上式中是从先前迭代中的转移概率。在期望值化的下次迭代中，最大化步骤中得到的所有值构成一个新的参数向量如下式：

上式中j＝1,2；l＝1,2，当|θ^(k+1)-θ^(k)＜δ|时算法终止；

即，当P(R_i＝j|x₁,x₂,...,x_n)＞50％时(其中j＝1,2)，信号x_i认为属于j状态。

采用上述方案，该发明方法结合计算机、加速度传感器及多通道数据采集分析仪组成的滚动轴承故障诊断系统能够对变转速/变负载故障轴承振动信号进行准确分析处理，实现对轴承的故障诊断。该系统能够有效地克服传统消噪方法难以诊断变转速/变负载条件下的振动信号的问题，并解决在非高斯噪声存在的情况下高能脉冲改变信号包络性使得传统包络法无法应用的问题，具有一定的工程应用价值。

下面结合附图对本发明作进一步描述。

附图说明

附图1为本发明的滚动轴承故障诊断系统图；

附图2为本发明的诊断方法流程图；

附图3为本发明噪声脉冲识别效果图；

附图4为本发明轴承故障信号处理结果，(a)为未使用模式转换方法处理结果，(b)为使用模式转换方法处理结果。

具体实施方式

本发明的具体实施例如图1-4所示是基于模式转换的变转速/变负载滚动轴承故障诊断系统，该轴承故障诊断系统可实现轴承故障智能诊断。

滚动轴承故障诊断系统用于检测滚动轴承的故障现象，并将故障振动信号上传到PC机。具体的方案中，通过故障轴承端盖上的加速度传感器采集参数信号，将这些采集的参数信号通过模式转换技术进行处理，从而判别出设备的状况，对这种设备状况的描述即为故障现象，如轴承的内圈、外圈、滚动体等故障现象。具体的计算方法可以下列公式计算：

外圈故障公式：

内圈故障公式：

滚动体故障公式：

式中：f_r为旋转频率，n为轴承滚动体数，φ为径向方向接触角，d为滚动体

平均直径，D为轴承的平均直径。

实例案例：

本案例处理复杂的变转速/变负载故障轴承，通过公式(15)-(17)计算本实验中的轴承故障为复合故障，其中包括内圈、外圈和滚动体故障，根据测试轴承的数据可计算出其内圈、外圈、滚动体的故障频率分别为197.05Hz、121.51Hz和 79.25Hz。

基于模式转换的变速/变负载滚动轴承故障诊断方法，其具体步骤为：

(1)信号初步分离

设定硬阈值，当固定邻域内信号X₁，X₂，...，X_n中的正负峰值超过给定阈值则被识别，并定义模式R₁为正峰，R₃为负峰，R₂为移除正峰和负峰后的信号模式，

计算R₂的高斯累计分布与经验累计分布曲线之间的距离并选取最小的距离作为新的硬阈值，并重新带入原信号再次分类得到新的R₁、R₂、R₃；距离采用 Kolmogorov-Smirnov统计公式计算：

上式中Φ(·)为标准高斯分布的累计分布函数，F_n(·)为经验分布函数，其定义如下：

上式中1_A为A集的指示器；

模式识别后使用不同的统计方法(即本领域常规使用的方法)分别分析 R₁、R₂、R₃并获得其分布状态；对R₂使用Kolmogorov-Smirnov测试验证其服从高斯分布，如果P值超过置信程度0.05，则样本的高斯性需验证；

R₁、R₃的分布特性分析基于其经验右尾的行为状态，对于样本 X₁，X₂，...，X_n，用公式1-F_n(t)计算其经验右尾，

上式中F_n(t)为经验累积分布函数，并将其与幂函数bt^-γ进行比较，

基于以上分析，可以得到信号R₁、R₃的分布状态服从帕累托分布，帕累托分布可用其概率密度函数定义，公式如下：

上式中γ为尺度参数，x_m为转移参数。对于帕累托分布，理论右尾用如下公式表示：

R₁、R₃可视为具有相似参数的帕累托分布；

由以上分析，可得到两种分布状态，分别为高斯分布和帕累托分布；

(2)模式转换

经过初步分离，确定两种模式分别为帕累托分布与高斯分布，分别对应着 R₁、R₂，用下式表达：

上式中R_i为第i个信号的状态变量，且所观察到信号的统计特性取决于实际状态R_i；sP(γ,x_m)表示帕累托分布，可用概率密度函数定义如下：

通过过渡矩阵P将信号从一个状态转换到另一个状态，其中包含在时刻i时的状态l过渡到时刻i+1时状态j的概率p_1j＝P(R_i+1＝j|R_i＝l)，结果获得分别来自状态R₁和状态R₂的两个概率，分别为P(R_i＝1|x₁,x₂,...,x_n)和P(R_i＝2|x₁,x₂,...,x_n)，分离方法正是基于这些概率。

模型参数的估计是基于期望-最大化算法，其包含二步，分别是期望值化 (E-Step)和最大化(M-Step)；每次迭代中，状态值先根据先前计算的参数进行推论；然后根据这些推论估计新的参数集；重复迭代，直到到达局部最大似然函数，该算法的具体过程如下：

初始化：设置初始参数矢量j的取值为1和 2，

期望值化(E-Step)：假设θ^(k)是先前迭代中最大化迭代法(M-Step)所累积的参数向量；具体包含以下步骤：

1)向前过滤：当i＝1,2，...，n时迭代方程如下：

上式中x_i＝(x₁,x₂,...,x_i)为信号矢量，gj(x_i；θ^(k))为当测量值来自状态j时时刻i下的信号概率密度函数(g₁和g₂分别为对称的帕累托和高斯密度)；当 P(R_n＝j|x_n；θ^(k))累计时有下式：

迭代的起点为

4)向后过滤：当i＝n-1,n-2，...，1时迭代方程如下：

最大化(M-Step)：新参数θ^(k+1)的估计可以通过最大似然函数推导。实际考虑以下：

由于x_m不存在最大似然估计，因此其可任意选择为一些常数(该常数的选择为本领域常规技术手段，故不赘述)，

最后，得到过渡概率根据以下公式计算：

上式中是从先前迭代中的转移概率。在期望值化的下次迭代中，最大化步骤中得到的所有值构成一个新的参数向量如下式：

上式中j＝1,2；l＝1,2，当|θ^(k+1)-θ^(k)＜δ|时算法终止；

即，当P(R_i＝j|x₁,x₂,...,x_n)＞50％时(其中j＝1,2)，信号x_i认为属于j状态，噪声脉冲识别结果如图3所示。

将原始信号与识别的噪声脉冲相减并用包络法进行分析处理，识别结果如图4所

示，相比于未利用模式转换的故障信号来说，处理过后的轴承信号信噪比进一步提高，更容易提取周期信号。

从图4中可以看出外圈故障频率匹配为120Hz,滚动体故障的二倍频符合的频率则为160Hz,内圈故障与202.5Hz也基本符合。此结果进一步说明本发明对于具有较好的处理效果，值得推广应用。

本发明不局限于上述具体实施方式，本领域一般技术人员根据本发明公开的内容，可以采用其他多种具体实施方式实施本发明的，或者凡是采用本发明的设计结构和思路，做简单变化或更改的，都落入本发明的保护范围。

Claims (1)

1.一种基于模式转换的变速/变负载滚动轴承故障诊断方法，其特征在于，包括以下内容：首先假设提取的振动信号至少包含两个具有不同分布特性的信号源并据此建立振动信号模型；其次基于尾估计对信号进行初步分离分析，获得具有不同分布特性的信号模式；然后利用模式转换迭代参数以建立隐马尔可夫模型获得噪声脉冲，并将原始信号与之相减，进一步提高信号的信噪比；最后利用传统包络谱分析信号，评估滚动轴承是否有损伤以及该损伤的类型；

具体步骤为：

(1)信号初步分离

设定硬阈值，当固定邻域内信号X₁，X₂，...，X_n中的正负峰值超过给定阈值则被识别，并定义模式R₁为正峰，R₃为负峰，R₂为移除正峰和负峰后的信号模式，

计算R₂的高斯累计分布与经验累计分布曲线之间的距离并选取最小的距离作为新的硬阈值，并重新带入原信号再次分类得到新的R₁、R₂、R₃；距离采用Kolmogorov-Smirnov统计公式计算：

上式中Φ(·)为标准高斯分布的累计分布函数，F_n(·)为经验分布函数，其定义如下：

上式中1_A为A集的指示器；

模式识别后使用不同的统计方法分别分析R₁、R₂、R₃并获得其分布状态；对R₂使用Kolmogorov-Smirnov测试验证其服从高斯分布，如果P值超过置信程度0.05，则样本的高斯性需验证；

R₁、R₃的分布特性分析基于其经验右尾的行为状态，对于样本X₁，X₂，...，X_n，用公式1-F_n(t)计算其经验右尾，

上式中F_n(t)为经验累积分布函数，并将其与幂函数bt^-γ进行比较，

基于以上分析，可以得到信号R₁、R₃的分布状态服从帕累托分布，帕累托分布可用其概率密度函数定义，公式如下：

上式中γ为尺度参数，x_m为转移参数；

对于帕累托分布，理论右尾用如下公式表示：

R₁、R₃可视为具有相似参数的帕累托分布；

由以上分析，可得到两种分布状态，分别为高斯分布和帕累托分布；

(2)模式转换

经过初步分离，确定两种模式分别为帕累托分布与高斯分布，分别对应着R₁、R₂，用下式表达：

上式中R_i为第i个信号的状态变量，且所观察到信号的统计特性取决于实际状态R_i；sP(γ,x_m)表示帕累托分布，可用概率密度函数定义如下：

通过过渡矩阵P将信号从一个状态转换到另一个状态，其中包含在时刻i时的状态l过渡到时刻i+1时状态j的概率p_1j＝P(R_i+1＝j|R_i＝l)，结果获得分别来自状态R₁和状态R₂的两个概率，分别为P(R_i＝1|x₁,x₂,...,x_n)和P(R_i＝2|x₁,x₂,...,x_n)，

模型参数的估计是基于期望-最大化算法，其包含二步，分别是期望值化E-Step和最大化M-Step；每次迭代中，状态值先根据先前计算的参数进行推论；然后根据这些推论估计新的参数集；重复迭代，直到到达局部最大似然函数，该算法的具体过程如下：

初始化：设置初始参数矢量j的取值为1和2，

期望值化E-Step：假设θ^(k)是先前迭代中最大化M-Step所累积的参数向量；具体包含以下步骤：

1)向前过滤：当i＝1,2，...，n时迭代方程如下：

上式中x_i＝(x₁,x₂,...,x_i)为信号矢量，gj(x_i；θ^(k))为当测量值来自状态j时时刻i下的信号概率密度函数，g₁和g₂分别为对称的帕累托和高斯密度；当P(R_n＝j|x_n；θ^(k))累计时有下式：

迭代的起点为

2)向后过滤：当i＝n-1,n-2，...，1时迭代方程如下：

最大化M-Step：新参数θ^(k+1)的估计可以通过最大似然函数推导；实际考虑以下：

由于x_m不存在最大似然估计，因此其可任意选择为常数，

最后，得到过渡概率根据以下公式计算：

上式中是从先前迭代中的转移概率；在期望值化的下次迭代中，最大化步骤中得到的所有值构成一个新的参数向量如下式：

上式中j＝1,2；l＝1,2，当|θ^(k+1)-θ^(k)＜δ|时算法终止；

即，当P(R_i＝j|x₁,x₂,...,x_n)＞50％时(其中j＝1,2)，信号x_i认为属于j状态。