基于CEEMDAN与GWO-NLM的滚动轴承故障诊断方法

基于ceemdan与gwo
‑
nlm的滚动轴承故障诊断方法
技术领域
1.本发明涉及轴承故障诊断技术领域，尤其涉及一种基于ceemdan与gwo
‑
nlm的滚动轴承故障诊断方法。


背景技术：

2.滚动轴承被广泛地应用到各种旋转式的机械系统中，其运动状态会对整体机械系统的精确性、可靠度和使用寿命产生巨大的影响。由于轴承工作于高速，高载荷的情况下，使其较为容易发生故障，其工作状态可能会影响整个机器运转的性能和安全性。因此，滚动轴承的故障诊断对于旋转机械系统的健康监测是必不可少的。这将有助于减少紧急维护和生产有关的成本。在轴承工作过程中，由于安装不当、过载及润滑不良等原因都容易使轴承局部缺陷造成损伤，而在轴承早期故障时，及时诊断故障可以有效避免损伤的加剧。信号分析是机械故障诊断研究与应用中的一个关键课题。它是一种提取故障特征，然后识别故障模式的工具。在这一领域，目前被广泛应用的轴承故障检测方法主要有两种，即声信号分析和振动信号分析。其中，基于振动信号的分析因其易于测量而成为最流行的监测技术。然而，从振动数据分析轴承故障的主要障碍是信号易于被背景噪声所淹没，并与轴、齿轮等机械的振动信号混合，因此，设计一种有效的信号增强方法来减轻噪声，突出故障特征就显得十分重要。
3.针对滚动轴承信号具有的脉冲性、非线性的特点，经验模态分解(empirical mode decomposition，emd)可以有效针对轴承信号进行自适应处理提高信噪比，但存在模态混叠的问题。为了解决emd存在的模态混叠问题，集总经验模态分解(ensemble empirical mode decomposition，简称eemd)，与emd不同的是eemd在分解过程中通过多次添加高斯白噪声来解决模态混叠问题。但是由于高斯白噪声的引入，使得最终重构信号受到了残留白噪声的影响，无法准确的重构原始信号，不仅如此，在算法的迭代运行过程中，计算量大，分解时间长。目前通过小波变换和盲源分离已被证实对于分离噪声和故障信号有较好的效果，但存在着局部失真和部分有效信息丢失等问题。因此，在滚动轴承故障信号的提取方面还存在很多的不足，亟需一种更高效准确的轴承故障诊断技术。


技术实现要素：

4.针对上述现有技术的不足，本发明提供一种基于ceemdan与gwo
‑
nlm的滚动轴承故障诊断方法。其中，自适应的噪声完备集成经验模态分解(complete ensemble empirical mode decomposition with adaptive noise,ceemdan)，该方法通过在每次分解过程中加入自适应的高斯白噪声，计算唯一的信号余量。非局部均值滤波(non
‑
local means,nlm)是对传统邻域滤波方法的一种改进滤波，使用了图像的自相似性质，其滤波效果取决于三个参数：m(搜索框半径)、λ(带宽参数)、p(相似框半径)，参数的选取十分重要。
5.为解决上述技术问题，本发明所采取的技术方案是：一种基于ceemdan与gwo
‑
nlm的滚动轴承故障诊断方法，包括如下步骤：
6.步骤1：采用振动传感器采集滚动轴承的振动信号，同时测量轴承的相关参数；
7.所述轴承的相关参数包括：接触角α、滚珠数z、滚珠直径d、节径d。
8.步骤2：采用ceemdan方法对采集得到的滚动轴承的振动信号进行分解，分解出不同的imf分量；
9.步骤3：结合峭度值、相关系数和能量比，对分解后的imf分量进行筛选，并对筛选后的分量进行线性重构；
10.所述结合峭度值、相关系数和能量比，对分解后的imf分量进行筛选的过程如下：
11.步骤3.1：计算各imf分量的峭度值k，如下所示：
[0012][0013]
其中，x为imf分量的振幅，u为imf分量振幅的平均值；σ为imf分量振幅的标准差；
[0014]
步骤3.2：计算各imf分量与原始信号的皮尔逊相关系数r；
[0015]
步骤3.3：计算各imf分量的能量比系数，如下所示：
[0016][0017]
其中，e
x
为imf分量的总能量；e
imf
(i)为第i个imf分量所具有的能量；ε为能量比系数；
[0018]
步骤3.4：对计算后的峭度值k、相关系数r和能量比系数ε进行加权求和，得到综合筛选指标krε值，如下所示：
[0019]
krε＝a1k+a2r+a3ε
[0020]
其中，a1、a2、a3分别为imf分量的峭度值、相关系数和能量比的权重值；
[0021]
步骤3.5：对各imf分量的综合筛选指标krε值进行从大到小排序，筛选出最大的前n个imf分量。
[0022]
步骤4：采用灰狼算法gwo优化nlm的带宽参数和相似框半径，过程如下：
[0023]
步骤4.1：初始化nlm算法，设定带宽参数λ、相似框半径p和搜索框半径m；为了运算速度，选择固定搜索框半径m；
[0024]
步骤4.2：初始化灰狼算法中的灰狼种群、收敛因子a和系数向量a、c；
[0025][0026]
a＝2a
×
r2‑
α
[0027]
c＝2r1[0028]
其中，t是迭代次数，tmax是最大迭代次数，r1和r2的模取[0,1]区间内的随机数；
[0029]
步骤4.3：把峭度值计算函数作为灰狼算法的适应度函数，计算每个灰狼的适应度值，同时保存适应度较大的前三个灰狼的参数；
[0030]
步骤4.4：更新各灰狼的位置；
[0031]
步骤4.5：更新a、a、c；
[0032]
步骤4.6：计算全部灰狼的适应度值，同时更新适应度最好的前三匹灰狼的适应度和位置；
[0033]
步骤4.7：判断是否达到最大的迭代次数，如果达到则结束寻优迭代过程；如果未
达到则返回执行步骤4.4至步骤4.7继续寻优迭代过程，直到迭代结束输出最优的带宽参数λ、相似框半径p。
[0034]
步骤5：采用优化后的gwo
‑
nlm方法对步骤3重构后的信号进行降噪处理，输出最优去噪信号；
[0035]
步骤6：将输出的最优去噪信号通过sg平滑滤波进行二次降噪，对二次降噪后的信号进行包络谱分析，提取出故障特征频率；
[0036]
步骤7：结合轴承的相关参数计算出轴承的理论故障特征频率，将步骤6提取出的故障特征频率与理论故障特征频率进行对比，判断滚动轴承的故障类型。
[0037]
采用上述技术方案所产生的有益效果在于：
[0038]
1、本发明提供的方法所选用的自适应的噪声完备集成经验模态分解，通过在每次分解过程中加入自适应的高斯白噪声，计算唯一的信号余量，不仅使重构误差几乎为0，而且运算效率较高。
[0039]
2、本发明提供的方法所选用的峭度
‑
相关系数
‑
能量比准则筛选出分解的模态分量，可以有效去除无关分量，避免了通过单一指标选择的imf分量的片面性，达到抑制背景噪声的目的。
[0040]
3、本发明提供的方法所选用的运用峭度指标作为灰狼优化算法的目标函数，由于峭度对周期性冲击较为敏感，可以解决由于nlm算法参数的寻优问题，使nlm算法的去噪性能发挥最优。
[0041]
4、本发明提供的方法应用于仿真实验信号和实际滚动轴承振动信号，结果表明，经过该方法去噪得到的信号有效的抑制了噪声的干扰，信噪比得到提高，而且包含十分丰富且明显的故障信息，具有一定的工程应用价值。
附图说明
[0042]
图1为本发明实施例中基于ceemdan与gwo
‑
nlm的滚动轴承故障诊断方法的流程图；
[0043]
图2为本发明实施例中用于仿真轴承信号的去噪包络谱图；
[0044]
图3为本发明实施例中实验轴承内圈故障的去噪包络谱图的整体图；
[0045]
图4为本发明实施例中实验轴承内圈故障的去噪包络谱图的局部细节图；
[0046]
图5为本发明实施例中实验轴承外圈故障的去噪包络谱图的整体图；
[0047]
图6为本发明实施例中实验轴承外圈故障的去噪包络谱图的局部细节图。
具体实施方式
[0048]
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
[0049]
本实施例采用了真实的实验数据进行分析，取自西储大学轴承数据中心。选择用于分析的故障轴承为6205
‑
2rjem skf型深沟球轴承，利用电火花技术对轴承内圈和外圈的故障进行单部位损伤加工。振动数据的采样频率为12000hz。
[0050]
如图1所示，本实施例中基于ceemdan与gwo
‑
nlm的滚动轴承故障诊断方法如下所述。
[0051]
步骤1：采用振动传感器采集滚动轴承的振动信号，同时测量轴承的相关参数；
[0052]
所述轴承的相关参数包括：接触角α、滚珠数z、滚珠直径d、节径d。
[0053]
步骤2：采用ceemdan方法对采集得到的滚动轴承的振动信号进行分解，分解出不同的imf分量；
[0054]
本实施例中在原信号中加入白噪声，对加入噪声的原始信号进行emd分解，将分解得到的imf
1i
分量进行平均得到imf1：
[0055]
x(t)＝y(t)+z0n
i
(t)
[0056][0057]
式中：y为原始信号；z0为所添加白噪声的幅值；i为构造信号次数，i∈z
+
；n为集合平均总次数，n∈z
+
；n
i
为满足n(0,1)分布的白噪声。
[0058]
第一阶段的残余分量可以得出：
[0059]
x1(t)＝x(t)
‑
imf1[0060]
式中：x1(t)为第一阶段的残余分量。
[0061]
定义一种新的运算emd
k
(
·
)意为由emd分解后所得到的第k个imf分量。给剩余信号x1(t)加上白噪声x1(t)+z1emd1(n
i
(t))，然后通过emd分解，对第一个imf取平均值可以得到imf2:
[0062][0063]
式中：imf2为第二阶段模态分量。
[0064]
对于k＝1,2,...,k，第k个残余分量可以得出：
[0065]
x
k
(t)＝x
k
‑1(t)
‑
imf
k
[0066]
反复以上步骤，直到剩余的信号不能被emd分解，分解完成后就能得到许多模态分量和剩余信号，其结果是：
[0067][0068]
式中：k为模态分解过程中模态的总数。
[0069]
步骤3：结合峭度值、相关系数和能量比，对分解后的imf分量进行筛选，并对筛选后的分量进行线性重构；
[0070]
所述结合峭度值、相关系数和能量比，对分解后的imf分量进行筛选的过程如下：
[0071]
步骤3.1：计算各imf分量的峭度值k，如下所示：
[0072][0073]
其中，x为imf分量的振幅，u为imf分量振幅的平均值；σ为imf分量振幅的标准差；
[0074]
步骤3.2：计算各imf分量与原始信号的皮尔逊相关系数r；
[0075]
步骤3.3：计算各imf分量的能量比系数，如下所示：
[0076]
[0077]
其中，e
x
为imf分量的总能量；e
imf
(i)为第i个imf分量所具有的能量；ε为能量比系数；
[0078]
步骤3.4：对计算后的峭度值k、相关系数r和能量比系数ε进行加权求和，得到综合筛选指标krε值，如下所示：
[0079]
krε＝a1k+a2r+a3ε
[0080]
其中，a1、a2、a3分别为imf分量的峭度值、相关系数和能量比的权重值；
[0081]
步骤3.5：对各imf分量的综合筛选指标krε值进行从大到小排序，筛选出最大的前n个imf分量。
[0082]
本实施例中，相关系数的求解采用皮尔逊相关系数。相关系数是[
‑
1,1]之间的数值，量化地描述x和y的相关程度，即相关系数越大，相关度越大，在相关系数为0时对应度最低。
[0083]
设两样本分别为x和y，相关系数为：
[0084][0085]
式中：r为两样本的相关系数；cov(x,y)为两样本的协方差；为x的方差；为y的方差。
[0086]
本实施例中，计算综合筛选指标krε值时，对三种不同的指标取相同的权重进行了分析和计算，并对其进行归一化处理，计算出各imf分量在ceemdan分解后的krε值，并对其进行排序，取最大值的前六个imf分量进行线性重构，得到重构信号y。
[0087]
步骤4：采用灰狼算法gwo优化nlm的带宽参数和相似框半径，gwo算法其灵感来自于狼的分工和协同觅食。它是一种新的群体智能算法，模拟狼的层次结构和狼的捕食行为。排名最高的狼是a种狼，位于食物链的顶端，负责领导、决策和其他行为，其次是b种狼、c种狼和e种狼。虽然b种狼和c种狼不是排名最高的狼种，但当a种狼失去领导权时，它们可以接替a狼成为新的领导者。e种狼是狼群中等级最低的一种，负责平衡种群内部的关系。
[0088]
gwo算法将每只狼视为一个潜在解，其中a种狼为第一最优解，狼种b和c分别为第二最优解和次最优解。gwo算法是一个不断更新狼种a、b、c位置的迭代优化过程。狼群通过下述公式更新距离和位置，完成对猎物的搜索。
[0089]
d＝|c
×
x
p
(t)
‑
x(t)|
[0090]
x(t+1)＝x
p
(t)+d
[0091]
其中，d是灰狼与猎物的距离，t是迭代次数；x
p
是猎物的位置，x是灰狼的位置，其初始位置坐标定义为(c,g)。系数向量a和c的表达式是a＝2a
×
r2‑
α，c＝2r1。当|a|>1表示全局搜索时，即灰狼群体扩大搜索范围以更方便寻找猎物。相比之下，|a|≤1，代表局部搜索，灰狼会收缩包围圈，就近寻找猎物。随迭代数的增加，收敛因子a由2线性减至0，而tmax是最大迭代次数，r1和r2的模分别在[0,1]区间内随机取值。
[0092]
在狼群判断猎物的位置时，头狼a领导狼b和狼c包围着猎物，因为狼a、狼b、狼c是最接近猎物，所以三只狼的位置是逐渐接近猎物，他们描述如下：
[0093]
d
a
＝|c1×
x
a
(t)
‑
x(t)|
[0094]
d
b
＝|c2×
x
b
(t)
‑
x(t)|
[0095]
d
c
＝|c3×
x
c
(t)
‑
x(t)|
[0096]
其中，x
a
代表a狼的当前位置,x
b
代表b狼的当前位置,x
c
代表c狼的当前位置。c1,c2,c3是随机变量。x(t)是狼群当前的位置。狼e到狼a、b、c的步长分别记为x1、x2和x3，如下所示：
[0097]
x1＝|c3×
x
a
‑
a1d
a
|
[0098]
x2＝|c2×
x
b
‑
a2d
b
|
[0099]
x3＝|c3×
x
c
‑
a2d
c
|
[0100]
狼e的最终位置定义如下：
[0101][0102]
在狼群捕猎时，狼a、狼b和狼c的对猎物的适应程度不同。对不同的适应度值进行计算，得出第一最优解、第二最优解和次优解，并保留其当前的位置信息。与此同时，狼群根据这三组位置信息来判断移动方向，以接近猎物完成捕猎。然后，再一次更新灰狼的位置，直到得到最优解。定义最优解对应的位置坐标值为(bestc,bestg)。
[0103]
步骤4的具体过程如下：
[0104]
步骤4.1：初始化nlm算法，设定带宽参数λ、相似框半径p和搜索框半径m；为了运算速度，选择固定搜索框半径m；
[0105]
步骤4.2：初始化灰狼算法中的灰狼种群、收敛因子a和系数向量a、c；
[0106][0107]
a＝2a
×
r2‑
α
[0108]
c＝2r1[0109]
其中，t是迭代次数，tmax是最大迭代次数，r1和r2的模取[0,1]区间内的随机数；
[0110]
步骤4.3：把峭度值计算函数作为灰狼算法的适应度函数，计算每个灰狼的适应度值，同时保存适应度较大的前三个灰狼的参数；
[0111]
步骤4.4：更新各灰狼的位置；
[0112]
步骤4.5：更新a、a、c；
[0113]
步骤4.6：计算全部灰狼的适应度值，同时更新适应度最好的前三匹灰狼的适应度和位置；
[0114]
步骤4.7：判断是否达到最大的迭代次数，如果达到则结束寻优迭代过程；如果未达到则返回执行步骤4.4至步骤4.7继续寻优迭代过程，直到迭代结束输出最优的带宽参数λ、相似框半径p。
[0115]
步骤5：采用优化后的gwo
‑
nlm方法对步骤3重构后的信号进行降噪处理，输出最优去噪信号；
[0116]
假设带噪声信号y由真实信号u加噪声信号n组成：y＝u+n；
[0117]
在样本i处的nlm去噪值是由其邻域ω的所有样本的加权平均得到的。
[0118]
[0119]
其中，是一个归一化常数。权重ω(i,j)，满足0<ω(i,j)<1和σ
j
ω(i,j)＝1，它比较了样本i和j周围的邻域。如果样本i的邻域与样本j的邻域相似，则取大值，否则取小。权重ω(i,j)的计算方法是：
[0120][0121]
其中，δ是一个带宽参数，δ代表一个局部的样本补丁。表示以j为中心的斑块中的样本的点对点差异的平方和。
[0122]
在上述公式中，不考虑i＝j这种特殊情况，因为权重ω(i,j)可以比其他每个样本的权重大得多。因为每个邻域都与自己相似(ω(i,j)＝1)，为了防止样本i对自己的权重过高，让ω(i,j)等于其他样本的最大权重。
[0123][0124]
其中，搜索框的半径越大去噪效果越好，但是过大会导致算法运行时间过长，因此权衡算法运行时间和优化算法时间，选择固定搜索框半径，通过优化带宽参数、相似框半径提高滤波器的降噪效果，其中带宽参数对去噪效果起着决定性作用。
[0125]
步骤6：将输出的最优去噪信号通过sg平滑滤波进行二次降噪，对二次降噪后的信号进行包络谱分析，提取出故障特征频率；
[0126]
步骤7：结合轴承的相关参数计算出轴承的理论故障特征频率，将步骤6提取出的故障特征频率与理论故障特征频率进行对比，判断滚动轴承的故障类型。
[0127]
本实施例中，用于仿真轴承信号的去噪包络谱如图2所示，可以看出降噪后信号的包络谱图突出了故障特征频率和倍频成分。
[0128]
实验轴承内圈故障的去噪包络谱图的整体如图3所示，可以看出经过gwo
‑
nlm
‑
sg降噪后信号包络谱可以提取到十分丰富的信息，故障频率f
i
为161.5hz、二倍频2f
i
＝323.4hz、三倍频3f
i
＝484.9hz、四倍频4f
i
＝646.7hz、五倍频5f
i
＝808.2hz、六倍频6f
i
＝970.1hz、七倍频7f
i
＝1132.0hz。转频f
r
为30.0hz、二倍频2f
r
＝59.7hz、三倍频3f
r
＝89.7hz、四倍频4f
r
＝119.8hz、五倍频5f
r
＝149.4hz、六倍频6f
r
＝179.4hz。在0
‑
1000hz区间，可以明显看出故障频率的峰值最大，其对应倍频的幅值先下降，四倍频达到最低，之后开始上升，在五倍频后开始下降。
[0129]
实验轴承内圈故障的去噪包络谱图的局部细节如图4所示，不仅能够提取到转频、故障频率及它们对应的倍频，还可以提取出以倍频和故障频率为中心丰富的调制频率。其中出现了明显的转频及其倍频nf
r
(其中n＝1，2，3，4，5，6，7)、故障频率及其倍频mf
i
(其中m＝1，2)、及以故障频率和倍频为中心的调制频率f
i
‑
af
r
(其中a＝1，2，3，4，5)、f
i
+bf
r
(其中b＝1，2，3，4，5，6)、2f
i
‑
cf
r
(其中c＝1，2，3，4，5，6，7，8)、2f
i
+df
r
(其中d＝1，2)、3f
i
‑
ef
r
(其中e＝3，4，5，6，7，8，9)。通过综合分析，可以清楚地分析故障为滚动轴承内圈的故障。
[0130]
实验轴承外圈故障的去噪包络谱图的整体如图5所示，可以看出该方法有效的提取出转频及2倍频、3倍频、4倍频、5倍频、内圈故障频率及2倍频到11倍频都被清晰的提取了出来，外圈故障频率的幅值的趋势为先下降，在故障频率的三倍频处达到局部最低值，而后对应其倍频幅值逐渐升高，在故障频率的五倍频处达到局部最高值，随着逐渐下降。综上，
可以判断出是滚动轴承的外圈发生故障。
[0131]
实验轴承外圈故障的去噪包络谱图的局部细节如图6所示，可以明显看出降噪后的信号包含着十分丰富的故障信息都被提取出来。