基于新型信噪比指标的随机共振滚动轴承故障诊断方法与流程

本发明涉及滚动轴承微弱故障特征信息提取领域，尤其涉及一种基于新型信噪比指标的自适应随机共振滚动轴承故障诊断方法。

背景技术

机械设备的振动信号包含故障特征信息，但振动信号通常受到强噪声的干扰，从而导致振动信号中的故障特征信息被强噪声淹没，造成难以提取故障特征信息及判别故障类型。

因此，提取强噪声背景下的故障特征信息是振动故障诊断领域关键问题之一。在故障振动特征信息提取方面，传统的方法是通过抑制噪声来提取故障特征信息及判别故障类型，但是该方法会导致原信号中故障特征信息的破坏，对故障特征信息的提取及故障类型判别造成极大影响。针对该问题，前人提出了随机共振的方法。该方法能将部分噪声能量会转化为信号能量，以便提取信号的故障特征信息及判别故障类型。目前，双稳态自适应随机共振的研究比较普遍，该方法通常以经典信噪比公式作为评价指标，来选取最优系统参数。尽管经典信噪比指标在自适应随机共振中起到了重要作用，但是计算该指标时需要提前知道故障实际特征频率的精确值，而在实际应用中是不可能提前知道故障实际特征频率的精确值，尤其在扰动工况下，故障特征频率是波动的，该指标就更不能发挥作用。

所以，在实际工程应用中要解决传统自适应随机共振方法需要提前知道故障实际特征频率精确值的问题，就必须提出一种新型信噪比的计算方法，使其在不知道故障实际特征频率精确值的前提下实现特征频率提取，并将其应用于滚动轴承故障类型判别及故障特征信息提取技术中。

技术实现要素：

发明目的：为了克服现有技术中存在的不足，本发明提供一种为了克服经典信噪比需要提前知道故障实际特征频率精确值的不足，本发明提供了一种基于新型信噪比指标的自适应随机共振滚动轴承故障诊断方法。该方法不用提前知道故障特征频率，就能实现故障特征的提取及故障类型判别，且能够更好的抑制边频，而且不受外载荷波动、轴承转速以及轴承型号的影响。同时，本发明也为由载荷变化等因素引起的小范围内故障特征频率振荡的情况下提取故障特征信息及判别故障类型提供了一种有效的方法。

技术方案：为实现上述目的，本发明的基于新型信噪比指标的随机共振滚动轴承故障诊断方法，其特征在于：包括如下几个步骤：

步骤1：信号采集；

步骤2：信号预处理；

步骤3：滚动轴承故障理论特征频率计算；

步骤4：基于新型信噪比指标的自适应随机共振方法；

步骤5：故障实际特征频率提取及故障类型判别。

进一步的，所述步骤2包括采用高通滤波滤去低频成分以及采用普通变尺度方法使采集到的振动信号满足经典随机共振的小参数要求。

进一步的，所述步骤3包括根据轴承的结构参数和实际转速计算外圈、内圈、滚动体和保持架四种典型故障的理论特征频率。

进一步的，所述步骤4中的新型信噪比作为评价指标对应的计算公式为：

式中：snri表示该新型信噪比；ft表示故障理论特征频率；kt表示ft对应的数字序列；x(·)表示幅值谱；s(ft)表示以ft为中心，在ft-lδft和ft+lδft之间的信号能量，2l表示计算带宽，δft表示频率分辨率；n(ft)表示噪声能量；m表示计算的数据序列的最大位置。

进一步的，所述步骤4中的自适应随机共振对应的朗之万方程为：

式中：u(x)表示双稳态势函数；s(t)表示输入信号；n(t)表示噪声强度为d的高斯白噪声；δ(t)表示狄拉克函数；a、b表示正系统参数，可以通过调节系统参数a、b来改变势阱势垒的深度和宽度；所述朗之万方程采用四阶龙格库塔算法进行求解，其离散公式如下：

式中：sj、nj以及xj分别表示输入信号、噪声以及输出信号的第j个采样值；h表示迭代步长。

进一步的，系统参数a、b采用增加收缩因子的粒子群优化算法，其具体步骤如下所示：

1)设置初始条件：学习因子c1和c2都设置为2.05，最大迭代次数设置为40，初始化群体个数设置为40。

2)初始化种群个体：随机初始化粒子的位置和速度。

3)计算每个粒子的适应度并找出局部最优和全局最优。

4)进入主循环：首先，更新粒子的速度和位置。接着，重新计算各粒子的适应度，并且更新局部最优值和全局最优值：然后判断最优值是否在0到2之间。最后判断是否达到最大迭代次数，若达到则输出最优解，否则继续循环。

进一步的，所述更新粒子的速度是指根据该粒子群优化算法的计算公式来更新，计算公式如下所示：

式中，为收缩因子，vi.j(t)和xi.j(t)分别表示速度和位置，r1和r2表示随机权重，pi.j表示局部最优，pg.j表示全局最优。

进一步的，所述适应度对应的是所述新型信噪比指标。

有益效果：相比于经典信噪比指标需要提前知道故障实际特征频率，该方法不需要提前知道故障实际特征频率的精确值，可以通过故障理论特征频率在其周围搜索到故障实际特征频率。

能较好的抑制边频，提高输出的信噪比，进而较好的提取故障特征信息及判别故障类型。

该方法不受外载荷、轴承转速以及轴承型号的影响，能普遍适用于机械设备中的滚动轴承故障诊断。

该方法可以适用于波动工况下引起的特征频率小范围振荡的故障特征频率及故障类型的提取。

附图说明

图1为本发明专利中基于新型信噪比指标的自适应随机共振滚动轴承故障诊断方法的总体流程图；

图2a为本实验室的实验信号的时域图和频谱图；

图2b为本实验室的实验信号通过高通滤波后的时域图和频谱图；

图3a为本实验室的实验信号用该方法在外圈故障理论特征频率处计算新型信噪比输出的频谱图；

图3b为本实验室的实验信号用该方法在内圈故障理论特征频率处计算新型信噪比输出的频谱图；

图3c为本实验室的实验信号用该方法在滚动体故障理论特征频率处计算新型信噪比输出的频谱图；

图3d为本实验室的实验信号用该方法在保持架故障理论特征频率处计算新型信噪比输出的频谱图；

图4a为本实验室的实验信号在工况1下的频谱图；

图4b为本实验室的实验信号在工况1下用该方法处理后输出的频谱图；

图4c为本实验室的实验信号在工况2下的频谱图；

图4d为本实验室的实验信号在工况2下用该方法处理后输出的频谱图；

图4e为本实验室的实验信号在工况3下的频谱图；

图4f为本实验室的实验信号在工况3下用该方法处理后输出的频谱图；

图4g为本实验室的实验信号在工况4下的频谱图；

图4h为本实验室的实验信号在工况4下用该方法处理后输出的频谱图；

图5a为本实验室的实验信号在转速为900rpm下的频谱图；

图5b为本实验室的实验信号在转速为900rpm下用该方法处理后输出的频谱图；

图5c为本实验室的实验信号在转速为1200rpm下的频谱图；

图5d为本实验室的实验信号在转速为1200rpm下用该方法处理后输出的频谱图；

图5e为本实验室的实验信号在转速为1500rpm下的频谱图；

图5f为本实验室的实验信号在转速为1500rpm下用该方法处理后输出的频谱图；

图6a为凯斯西储大学实验室的实验信号的时域图和频谱图；

图6b为凯斯西储大学实验室的实验信号通过高通滤波后的时域图和频谱图；

图7a为凯斯西储大学实验室的实验信号用该方法在外圈故障理论特征频率处计算新型信噪比输出的频谱图；

图7b为凯斯西储大学实验室的实验信号用该方法在内圈故障理论特征频率处计算新型信噪比输出的频谱图；

图7c为凯斯西储大学实验室的实验信号用该方法在滚动体故障理论特征频率处计算新型信噪比输出的频谱图；

图7d为凯斯西储大学实验室的实验信号用该方法在保持架故障理论特征频率处计算新型信噪比输出的频谱图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如附图1至7所示，下面结合附图对本发明作更进一步的说明。应当注意的是，以下具体实施方式中描述的技术特征或者技术特征的组合应当被相互结合从而达到更好的技术效果，而不应当被孤立。

如附图1所示，本发明的一种基于新型信噪比指标的自适应随机共振滚动轴承故障诊断方法，包括如下步骤：

步骤1)信号采集。

步骤2)信号预处理：

首先，对采集到的振动信号采用高通滤波技术滤去低频成分；然后，采用普通变尺度方法使其满足经典随机共振的小参数要求。

步骤3)滚动轴承故障理论特征频率计算：

根据轴承的设计参数和实际转速计算外圈、内圈、滚动体和保持架四种典型故障的理论特征频率。

步骤4)基于新型信噪比的自适应随机共振：

以新型信噪比作为随机共振响应的评价指标，其对应的计算公式为：

式中：snri表示该新型信噪比；ft表示故障理论特征频率；kt表示ft对应的数字序列；x(·)表示幅值谱；s(ft)表示以ft为中心，在ft-lδft和ft+lδft之间的信号能量，2l表示计算带宽，δft表示频率分辨率；n(ft)表示噪声能量；m表示计算的数据序列的最大位置。

随机共振技术对应的朗之万方程为：

式中，u(x)表示双稳态势函数；s(t)表示输入信号；n(t)表示噪声强度为d的高斯白噪声；δ(t)表示狄拉克函数；a、b表示正系统参数，可以通过调节系统参数a、b来改变势阱势垒的深度和宽度，进而实现随机共振。

上述朗之万方程可以采用四阶龙格库塔算法进行求解，其离散公式如下：

式中，sj、nj以及xj分别表示输入信号、噪声以及输出信号的第j个采样值；h表示迭代步长。

在基于新型信噪比指标的自适应随机共振中，采用增加收缩因子的粒子群优化算法对系统参数a、b进行同步优化，进而实现自适应；以新型信噪比指标为适应度函数的增加收缩因子的粒子群优化算法的具体步骤如下所示：

(1)设置初始条件：学习因子c1和c2都设置为2.05，最大迭代次数4设置为40，初始化群体个数设置为40；

(2)初始化种群个体：随机初始化粒子的位置和速度；

(3)计算每个粒子的适应度并找出局部最优和全局最优；

(4)进入主循环：首先，根据该粒子群优化算法的寻优速度更新公式来更新粒子的位置和速度，计算公式如下所示；接着，重新计算各粒子的适应度值，并且更新局部最优值和全局最优值；然后，判断最优值是否在0到2之间；最后，判断是否达到最大迭代次数，若达到则输出最优解，否则继续循环。

其中，为收缩因子，vi.j(t)和xi.j(t)分别表示速度和位置，r1和r2表示随机权重，pi.j表示局部最优，pg.j表示全局最优。

步骤5)依据自适应随机共振最优输出实现故障实际特征频率提取及故障类型判别。

下面采用本实验室的实验信号对本发明方法进行测试：

实验参数：采样点数10240，采样频率2048hz，设置电动机转速为1500rpm，制动扭矩30nm，径向力150n。为了模拟强噪声背景下的随机共振，在实验信号中加入了强度为1的高斯白噪声。用红外线测速仪测得的实际转速为1466rpm，根据轴承外圈、内圈、滚动体以及保持架的故障理论特征频率计算公式，计算得外圈故障理论特征频率为108.54hz，内圈故障理论特征频率为176.78hz，滚动体故障理论特征频率为122.35hz，保持架故障理论特征频率为14.57hz。

图2a、2b分别为含噪实验信号在高通滤波前后的时域图和频谱图，可以看出高通滤波前后故障特征频率都被完全淹没，难以识别。

采用本发明所述的基于新型信噪比指标的自适应随机共振方法对滚动轴承实验信号进行处理。在预处理中，由于故障实际特征频率应在故障理论特征频率的某个邻域范围内，将通带频率和阻带频率分别设为105hz和100hz，普通变尺度的时间尺度为1000。在计算新型信噪比时，设置l＝15，j＝2，m＝2/n。由于外圈的故障理论特征频率大约是108.6hz，所以在105.6hz-111.6hz的范围内搜索故障实际特征频率；由于内圈的故障理论特征频率大约是176.8hz，所以在173.8hz-179.8hz的范围内搜索故障实际特征频率；由于滚动体的故障理论特征频率大约是122.4hz，所以在119.4hz-125.4hz的范围内搜索故障实际特征频率；由于保持架的故障理论特征频率大约是14.6hz，所以在11.6hz-17.6hz的范围内搜索故障实际特征频率。为了判别该信号的故障特征频率及故障类型，对以上四种典型的故障类型分别用本发明所提的方法进行处理，最终输出结果如图3a、3b、3c、3d所示。图3a是用所述方法在外圈故障理论特征频率处计算新型信噪比输出的频谱图，从图中可以发现故障实际特征频率107.8hz比较明显，在故障理论特征频率108.6hz的附近，且在105.6hz-111.6hz的范围内，因此初步判定该故障为外圈故障且故障实际特征频率为107.8hz。图3b、3c是分别用所述方法在内圈及滚动体故障理论特征频率处计算新型信噪比输出的频谱图，从图中依然可以发现频率107.8hz比较明显，但是该频率不在内圈故障理论特征频率176.8hz附近的搜索范围内，也不在滚动体故障理论特征频率122.4hz附近的搜索范围内，所以不是内圈或滚动体故障。图3d是用该方法在保持架故障理论特征频率处计算新型信噪比输出的频谱图，从图中可以看出本发明所提方法失效，并没有搜索到故障实际特征频率，所以不是保持架故障。因此最终判定该故障为外圈故障，且故障实际特征频率为107.8hz。

为了验证该发明在不同载荷下的有效性，采用了四种不同载荷工况来验证。工况1为空载，工况2设置制动扭矩30nm，工况3设置径向力300n，工况4设置制动扭矩30nm和径向力300n。实验参数：采样点数64000，采样频率12800hz，设置电动机转速为1500rpm。为了模拟强噪声背景下的随机共振，在实验信号中加入了强度为1的高斯白噪声。在工况1时，用红外线测速仪测得实际转速为1495rpm，根据轴承外圈故障特征频率计算公式得外圈故障理论特征频率为110.687hz；在工况2时，用红外线测速仪测得实际转速为1421rpm，计算得外圈故障理论特征频率为105.208hz；在工况3时，用红外线测速仪测得实际转速为1494rpm，计算得外圈故障理论特征频率为110.613hz；在工况4时，用红外线测速仪测得实际转速为1423rpm，计算得外圈故障理论特征频率为105.357hz。图4a、4c、4e、4g分别为四种工况下的含噪实验信号，从图中可以看出故障特征频率都被完全淹没，难以识别。

采用本发明所述的基于新型信噪比指标的自适应随机共振方法对四种工况下的信号进行处理。在预处理中，工况1和工况3下将通带频率和阻带频率分别设为106hz和102hz，工况2和工况4下将通带频率和阻带频率分别设为100hz和95hz，普通变尺度的时间尺度为1000。在计算新型信噪比时，设置l＝15，j＝2，m＝2/n。工况1时，由于外圈的故障理论特征频率大约是110.6hz，所以在107.6hz-113.6hz的范围内搜索故障实际特征频率；工况2时，由于外圈的故障理论特征频率大约是105.2hz，所以在102.2hz-108.2hz的范围内搜索故障实际特征频率；工况3时，由于外圈的故障理论特征频率大约是110.6hz，所以在107.6hz-113.6hz的范围内搜索故障实际特征频率；工况4时，由于外圈的故障理论特征频率大约是105.4hz，所以在102.4hz-108.4的范围内搜索故障实际特征频率。最终输出结果如图4b、4d、4f、4h所示，四种工况下的故障实际特征频率均能清晰地突显出来。因此，该发明所述方法不受外载荷的影响，依然可以检测到故障实际特征频率及故障类型。

为了验证该发明诊断效果不受轴承转速的影响，分别设置电动机转速为900rpm、1200rpm和1500rpm来验证。在这三种转速下，均设置为空载。在电动机转速为900rpm和1500rpm时，设置采样点数64000，采样频率12800hz。在电动机转速为1200rpm时，设置采样点数10240，采样频率2048hz。为了模拟强噪声背景下的随机共振，在实验信号中加入了强度为1的高斯白噪声。在设置电动机转速为900rpm时，用红外线测速仪测得实际转速为896.1rpm，根据轴承外圈故障特征频率计算公式得外圈故障理论特征频率为66.346hz；在设置电动机转速为1200rpm时，用红外线测速仪测得实际转速为1195rpm，计算得外圈故障理论特征频率为88.48hz；在设置电动机转速为1500rpm时，用红外线测速仪测得实际转速为1495rpm，计算得外圈故障理论特征频率为110.687hz。图5a、5c、5e分别为三种转速下的含噪实验信号，从图中可以看出故障特征频率都被完全淹没，难以识别。

采用本发明所述的基于新型信噪比指标的自适应随机共振方法对四种转速下的信号进行处理。在预处理中，转速为900rpm时将通带频率和阻带频率分别设为61hz和58hz；转速为1200rpm时将通带频率和阻带频率分别设为83hz和80hz；转速为1500rpm时将通带频率和阻带频率分别设为106hz和102hz，普通变尺度的时间尺度为1000。在计算新型信噪比时，设置l＝15，j＝2，m＝2/n。在电动机转速为900rpm的情况下，由于外圈的故障理论特征频率大约是66.4hz，所以在63.4hz-69.4hz的范围内搜索故障实际特征频率；在电动机转速为1200rpm的情况下，由于外圈的故障理论特征频率大约是88.4hz，所以在85.4hz-91.4hz的范围内搜索故障实际特征频率；在电动机转速为1500rpm的情况下，由于外圈的故障理论特征频率大约是110.6hz，所以在107.6hz-113.6hz的范围内搜索故障实际特征频率。最终输出结果如图5b、5d、5f所示，三种转速下的故障实际特征频率均能清晰地突显出来。因此，该发明所提方法不受轴承转速的限制，依然可以检测到故障实际特征频率及故障类型。

由于以上本实验室的实验数据都是通过型号为n306e的轴承测得的，为了验证该发明所述方法不受轴承型号的影响，采用了凯斯西储大学实验室的轴承故障数据进行验证，该轴承型号为6205-2rsjemskf。实验参数：采样点数60000，采样频率12000hz，转速为1730rpm，为了模拟强噪声背景下的随机共振，在实验信号中加入了强度为1的高斯白噪声。根据轴承外圈、内圈、滚动体以及保持架的故障理论特征频率计算公式，计算得外圈故障理论特征频率为103.36hz，内圈故障理论特征频率为156.14hz，滚动体故障理论特征频率为135.91hz，保持架故障理论特征频率为17.35hz。图6a、6b分别为含噪信号在高通滤波前后的时域图和频谱图，可以看出高通滤波前后故障特征频率都被完全淹没，难以识别。

采用本发明所述的基于新型信噪比指标的自适应随机共振方法对该实验信号进行处理。在预处理中，将通带频率和阻带频率分别设为100hz和98hz，普通变尺度的时间尺度为1000。在计算新型信噪比时，设置l＝15，j＝2，m＝2/n。由于外圈的故障理论特征频率大约是103.4hz，所以在100.4hz-106.4hz的范围内搜索故障实际特征频率；由于内圈的故障理论特征频率大约是156hz，所以在153hz-159hz的范围内搜索故障实际特征频率；由于滚动体的故障理论特征频率大约是136hz，所以在133hz-139hz的范围内搜索故障实际特征频率；由于保持架的故障理论特征频率大约是17.4hz，所以在14.4hz-20.4hz的范围内搜索故障实际特征频率。为了判别该信号的故障特征频率及故障类型，对以上四种典型的故障类型分别用本发明所述的方法进行处理，最终输出结果如图7a、7b、7c、7d所示。图7a是用该方法在外圈故障理论特征频率处计算新型信噪比输出的频谱图，从图中可以发现故障实际特征频率103.2hz比较明显，在故障理论特征频率103.4hz的附近，且在100.4hz-106.4hz的范围内，因此初步判定该故障为外圈故障，且故障实际特征频率为103.2hz。图7b、7c是分别用该方法在内圈及滚动体的故障理论特征频率处计算新型信噪比输出的频谱图，从图中依然可以发现频率103.2hz比较明显，但是该频率不在内圈故障理论特征频率156hz附近的搜索范围内，也不在滚动体故障理论特征频率136hz附近的搜索范围内，所以不是内圈或滚动体故障。图7d是用该方法在保持架故障理论特征频率处来计算新型信噪比输出的频谱图，从图中可以看出本发明所提方法失效，并没有搜索到故障实际特征频率，所以不是保持架故障。因此最终判定该故障为外圈故障，且故障实际特征频率为103.2hz。因此，该发明所述方法不受轴承型号的限制，依然可以检测到故障实际特征频率及判别故障类型。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。