基于tvd和med的滚动轴承故障诊断方法
技术领域
1.本发明属于机械故障诊断领域，具体涉及一种基于tvd和med的滚动轴承故障诊断方法。
背景技术：2.滚动轴承是机械设备的重要零部件，由于其长期处于受载状态和恶劣环境，无法保证运行工作时的稳定性。通常采集到的滚动轴承故障信号中包含噪声信号，采用单一分析方法和手段无法准确快速地判断故障类型。因此，如何在噪声信号中提取轴承故障特征成为主要研究内容。
3.滚动轴承发生故障所产生的振动信号一般为非线性、非循环平稳的高频信号，使得在提取故障频率时存在一定困难。常见的降噪方法有小波阈值降噪、自相关降噪等。小波阈值降噪需要确定基函数和阈值，其自适应性较差。自相关降噪虽然在一定程度上保证信号的相似性，但只适用于高信噪比信号。总变差降噪(total variation denoising，tvd)是利用正则化技术进行降噪，有效减少噪声的影响。隋文涛等将tvd应用于轴承故障领域，并提出将峭度和相关系数结合作为参数选择的指标。figueiredomat等采用mm算法解决tvd目标函数优化问题，取得较好的降噪效果，然而参数选择对去噪效果有较大影响。但总变差降噪作为一种降噪方法，并不能提取故障特征。endo等结合最小熵解卷积(minimum entropy deco nvolution，med)来放大故障特征，便于更加容易提取故障频率，但滤波器阶数大小会影响故障频率的提取。王晗等针对此问题结合轮对轴承进行深入研究设计最优滤波器来提取故障频率。由于轴承故障特征微弱，在强噪声背景下的故障信号会被噪声所淹没，严重影响故障特征提取。
技术实现要素：4.当信号经过tvd降噪时，如何在提高原始信号降噪效果的同时，保证其相关性和冲击特性，以及如何选取正则化参数a是机械故障信号分析的难点。
5.本发明的技术问题是如何选取正则化参数a提高信号降噪效果以及单一方法难以有效提取故障特征频率，影响轴承故障诊断效果。
6.本发明的目的是针对上述问题，提供一种基于tvd和med的滚动轴承故障诊断方法，该方法有效地剔除滚动轴承原始振动信号中含有的干扰噪声，得到清晰地反映实际故障信息的信号，降噪信号进行最小熵解卷积增强故障冲击特性，减少脉冲响应函数的传播，对故障信号的包络谱进行故障诊断。
7.本发明的技术方案是基于tvd和med的滚动轴承故障诊断方法，包括以下步骤：
8.步骤1：通过加速度传感器获取滚动轴承振动信号；
9.步骤2：计算滚动轴承振动信号的峰值因子初步判断轴承是否产生故障；
10.步骤3：对原始滚动轴承故障振动信号计算一阶总变差；
11.步骤4：利用相关系数和teager能量算子的加权指标选取最优的正则化系数a，并
结合一阶总变差得到降噪信号；
12.步骤5：采用最小熵解卷积增强降噪信号中的故障信号，减少脉冲响应函数的传播；
13.步骤6：最后对其包络解调得到包络谱图，提取清晰的故障特征频率及其倍频频率，与计算所得的滚动轴承故障特征频率理论值进行对比，诊断得到滚动轴承故障类型。
14.进一步地，步骤2中，滚动轴承振动信号的峰值因子的计算式如下：
15.x
p
＝max{|xi|}
ꢀꢀꢀꢀ
(1)
[0016][0017][0018]
式中x
p
表示信号x(t)的峰值，x
rms
表示信号x(t)的均方根值，c表示信号x(t)的峰值因子，xi表示第i个采样点，n表示采样点的个数。
[0019]
步骤3包括以下子步骤：
[0020]
步骤3.1：定义原始信号的一阶总变差：
[0021]
tv＝||dx||1ꢀꢀꢀꢀ
(4)
[0022]
式中tv表示信号一阶总变差；x表示原始振动信号；d表示n-1行n列的矩阵，n表示x的长度；||
·
||1表示l1范数；
[0023]
其中矩阵d定义如下：
[0024][0025]
步骤3.2：对步骤3.1得到的一阶总变差进行降噪，并得到目标优化问题如下：
[0026][0027]
式中tvd(y，a)表示总变差降噪函数；a为正则化参数；y为含有噪声的信号，x为原始振动信号；||
·
||1表示l1范数，||
·
||2表示l2范数；
[0028]
步骤3.3：利用mm算法(majorization-minimization)求解式(6)，得到如下结果：
[0029][0030]
式中xi、x
i+1
分别表示原始信号、降噪信号；diag()表示对角矩阵。
[0031]
步骤4包括以下子步骤：
[0032]
步骤4.1：计算原始信号的相关系数，其公式如下：
[0033][0034]
式中c表示互相关系数；和分别表示x和y的均值；σ
x
和σy分别表示x和y的标准差，e表示期望。
[0035]
步骤4.2：计算teager能量算子，其公式如下：
[0036]
t[x(n)]＝x(n)
2-x(n-1)x(n+1)
ꢀꢀ
(9)
[0037]
式中t()表示teager能量算子；x(n-1)、x(n)、x(n+1)分别表示连续的三个采样点；
[0038]
步骤4.3：为了使得降噪信号同时具备相关性和冲击特性，将二者指标结合得到加权指标如下：
[0039]
tc＝t
×cr
ꢀꢀ
(10)
[0040]
式中tc表示加权指标；t表示teager能量算子；r为调节参数，r在本发明中取1.2。
[0041]
步骤5中，最小熵解卷积的计算过程如下：
[0042]
1)假设振动信号y(n)可以表示如下：
[0043]
y(n)＝h(n)
·
x(n)+e(n)
ꢀꢀ
(11)
[0044]
式中h(n)为传递函数；x(n)为故障信号；e(n)为噪声信号
[0045]
2)通过寻找最优逆滤波器w(n)来消除卷积影响，得到与故障信号x(n)相近的逆卷积信号i(n)，计算公式如下：
[0046][0047][0048]
式中w(i)表示第i个逆波器；y(n-i)表示振动信号；l表示逆波器个数；
[0049]
3)通过计算逆卷积信号i(n)的范数衡量其熵的大小，并将其作为目标函数获得最优输出：
[0050][0051]
式中o2()表示逆卷积信号的熵值；i2(i)表示逆卷积信号幅值的平方；
[0052]
4)为使得式(14)范数最大，令其求导：
[0053][0054]
5)结合式(12)得到如下结果，并采用矩阵形式表示为：
[0055][0056]
c＝wa
ꢀꢀꢀꢀ
(17)
[0057]
w＝ca-1
ꢀꢀ
(18)
[0058]
式中w为逆滤波器的系数矩阵；c为输出信号与输入信号的互相关矩阵；a为输入信号自相关矩阵。
[0059]
步骤6中，滚动轴承内圈故障特征频率的计算式如下：
[0060][0061]
式中f0表示内圈特征频率的理论值，z代表滚动体个数，d表示滚动体直径，d表示节圆直径，α表示轴承接触角，fr表示转频；
[0062]
滚动轴承外圈故障特征频率的计算式如下：
[0063][0064]
式中f1表示外圈特征频率的理论值；
[0065]
转频fr的计算式如下：
[0066][0067]
式中ni表示转速。
[0068]
相比现有技术，本发明的有益效果包括：
[0069]
(1)本发明通过总变差降噪对原始信号降噪，得到清晰的反映实际故障信息的信号后，对降噪信号进行最小熵解卷积，降低脉冲响应函数的传播，利用重构信号的包络谱进行故障诊断，有效提高了滚动轴承故障诊断的准确性；
[0070]
(2)本发明采用相关系数和teager能量算子相结合作为选取正则化参数a的加权指标，在保证降噪信号与原始信号的相关性的同时，也能突出降噪信号的故障冲击特性，加权指标中teager能量算子相比于峭度具有更高运算效率，提高了信号的降噪效率；
[0071]
(3)本发明根据mm算法针对总变差目标优化问题进行求解，得到的结果通过迭代计算不断降低噪声影响，提高降噪效果；
[0072]
(4)本发明相较于tvd单一方法，tvd和med复合方法不仅能提取故障频率，对于高频范围的倍频故障频率也具有较好的提取效果。
附图说明
[0073]
下面结合附图和实施例对本发明作进一步说明。
[0074]
图1为本发明实施例的滚动轴承故障诊断方法的流程示意图。
[0075]
图2a为本发明实施例的仿真的原始信号的时频分布图。
[0076]
图2b为本发明实施例的加噪信号的时频分布图。
[0077]
图3为本发明实施例的仿真信号参数a对加权指标tc的影响。
[0078]
图4a为本发明实施例的tvd信号的时域波形图。
[0079]
图4b为本发明实施例的tvd-med信号的时域波形图。
[0080]
图5a为本发明实施例的tvd信号包络谱图。
[0081]
图5b为本发明实施例的tvd-med信号的包络谱图。
[0082]
图6a为本发明实施例的滚动轴承内圈故障振动信号的时域波形图。
[0083]
图6b为本发明实施例的滚动轴承内圈故障振动信号的频域波形图。
[0084]
图7为本发明实施例的滚动轴承内圈小波阈值去噪信号的时域波形图。
[0085]
图8为本发明实施例的滚动轴承内圈参数a对加权指标tc的影响。
[0086]
图9a为本发明实施例的滚动轴承内圈tvd信号的时域波形图。
[0087]
图9b为本发明实施例的滚动轴承内圈tvd-med信号的时域波形图。
[0088]
图10a为本发明实施例的滚动轴承内圈tvd包络谱图。
[0089]
图10b为本发明实施例的滚动轴承内圈tvd-med包络谱图。
[0090]
图11a为本发明实施例的滚动轴承外圈故障振动信号时域波形图。
[0091]
图11b为本发明实施例的滚动轴承外圈故障振动信号频域波形图。
[0092]
图12为本发明实施例的滚动轴承外圈小波阈值去噪信号时域波形图。
[0093]
图13为本发明实施例的滚动轴承外圈参数a对加权指标tc的影响。
[0094]
图14a为本发明实施例的滚动轴承外圈tvd信号时域波形图。
[0095]
图14b为本发明实施例的滚动轴承外圈tvd-med信号时域波形图。
[0096]
图15a为本发明实施例的滚动轴承外圈tvd包络谱图。
[0097]
图15b为本发明实施例的滚动轴承外圈tvd-med包络谱图。
具体实施方式
[0098]
实施例的轴承为6205-2rs jem skf深沟球轴承，轴承参数如表1所示。
[0099]
表1深沟球轴承参数表
[0100][0101]
实验采样频率fs＝12khz，电机转速为1750r/min，对应的转频fr＝29.17hz。根据式(19)计算处轴承内圈故障频率fi＝157.76hz。本发明的滤波器阶数为56。
[0102]
如图1所示，基于tvd和med的滚动轴承故障诊断方法，包括以下步骤：
[0103]
步骤1：通过加速度传感器获取滚动轴承振动信号；
[0104]
步骤2：计算滚动轴承振动信号的峰值因子初步判断滚动轴承是否产生故障；
[0105]
滚动轴承振动信号的峰值因子的计算式如下：
[0106]
x
p
＝max{|xi|}
ꢀꢀꢀꢀ
(1)
[0107][0108][0109]
式中x
p
表示信号x(t)的峰值，x
rms
表示信号x(t)的均方根值，c表示信号x(t)的峰值因子。
[0110]
步骤3：对原始滚动轴承故障振动信号计算一阶总变差；
[0111]
步骤3.1：定义原始信号的一阶总变差：
[0112]
tv＝||dx||1ꢀꢀꢀꢀ
(4)
[0113]
其中矩阵d定义如下：
[0114][0115]
步骤3.2：对步骤3.1得到的一阶总变差进行降噪，并得到目标优化问题如下：
[0116][0117]
式中a(a＞0)为正则化参数，其目的是控制信号的平滑程度，
[0118]
步骤3.3：利用mm(majorization-minimization)算法求解式(6)，得到如下结果：
[0119][0120]
步骤4：利用相关系数和teager能量算子的加权指标选取最优的正则化系数a，并得到最佳降噪信号；
[0121]
步骤4.1：计算原始信号的相关系数，其公式如下：
[0122][0123]
式中和分别表示x和y的均值；σ
x
和σy分别表示x和y的标准差，e表示期望。
[0124]
步骤4.2：计算teager能量算子，其公式如下：
[0125]
t[x(n)]＝x(n)
2-x(n-1)x(n+1)
ꢀꢀ
(9)
[0126]
步骤4.3：为了使得降噪信号同时具备相关性和冲击特性，将二者指标结合得到加权指标如下：
[0127]
tc＝t
×cr
ꢀꢀ
(10)
[0128]
式中r为调节参数，r在本发明中取1.2。
[0129]
步骤5：采用最小熵解卷积增强降噪信号中的故障信号，减少脉冲响应函数的传播；
[0130]
1)假设振动信号y(n)可以表示如下：
[0131]
y(n)＝h(n)
·
x(n)+e(n)
ꢀꢀ
(11)
[0132]
式中h(n)为传递函数；x(n)为故障信号；e(n)为噪声信号
[0133]
2)通过寻找一个最优逆滤波器w(n)来消除卷积影响，得到与故障信号x(n)相近的逆卷积信号i(n)，计算公式如下：
[0134][0135][0136]
3)通过计算信号i(n)的范数衡量其熵的大小，并将其作为目标函数获得最优输出：
[0137][0138]
4)为使得式(14)范数最大，令其求导：
[0139][0140]
5)结合式(12)得到如下结果，并采用矩阵形式表示为：
[0141][0142]
c＝wa
ꢀꢀꢀꢀ
(17)
[0143]
w＝ca-1
ꢀꢀ
(18)
[0144]
式中w为逆滤波器的系数矩阵；c为输出信号与输入信号的互相关矩阵；a为输入信号自相关矩阵。
[0145]
步骤6：最后对其包络解调得到包络谱图，提取清晰的故障特征频率及其倍频频率，与计算所得的滚动轴承故障特征频率理论值进行对比，诊断得到滚动轴承故障类型。
[0146]
滚动轴承内圈故障特征频率的计算式如下：
[0147][0148]
式中f0表示内圈特征频率的理论值，z代表滚动体个数，d表示滚动体直径，d表示节圆直径，α表示轴承接触角，fr表示转频；
[0149]
滚动轴承外圈故障特征频率的计算式如下：
[0150][0151]
式中f1表示外圈特征频率的理论值；
[0152]
转频fr的计算式如下：
[0153][0154]
式中ni表示转速。
[0155]
通过模拟滚动轴承出现故障后所产生的周期性冲击信号和高斯白噪声合成仿真信号表达式如下：
[0156]
x(t)＝x0exp(-2πfnξt)sin(2πfnt)+n(t)
ꢀꢀ
(22)
[0157]
式中位移常数x0＝5；载波频率fn＝3000hz；阻尼系数ξ＝0.1；采样频率fs＝20khz；采样点数n＝4096；冲击周期t＝0.01s；t为采样时间；n(t)为白噪声。
[0158]
原始信号和加噪信号的时域波形图如图2a、2b所示，加噪信号中周期性冲击特性并不明显，大部分有效信号都被噪声所覆盖。为此，对加噪信号进行tvd降噪，参数a对加权指标tc的影响如图3所示，由图3可以看出，当a＝3.1时，加权指标tc的值达到最大值。选取最佳参数a进行tvd降噪，对去噪后信号采用med方法增强信号的冲击特性，如图4a、4b所示。从图4a、4b中可看出tvd去噪消除了大量噪声的影响，信号出现周期冲击的特性。经med处理后，本发明的方法减少了脉冲信号的分散传播，使得信号的冲击特性更加明显。
[0159]
如图5a、5b所示，经tvd降噪后的频谱图虽然在低频范围能清晰分析出轴承故障频率，但在高频范围并不能提取处故障特征频率，而经过med后，在高频处同样能较为准确分析处故障频率。仿真结果表明本发明所提方法可用于在噪声背景下滚动轴承的故障特征提取。
[0160]
图6a、6b分别为滚动轴承内圈故障振动信号时域波形、频谱图，首先采用小波阈值降噪方法作为对比，得到小波阈值降噪信号的包络谱如图7所示。然后对内圈振动信号进行tvd降噪处理，参数a对加权指标的影响如图8所示。当a＝1.2时，加权指标tc值达到最大值，并选取参数a最大值进行tvd降噪，对降噪后信号采用med增强冲击特性，如图9a、9b所示。最后通过包络解调得到包络谱图，如图10a、10b所示。由图10a、10b可见，tvd和tvd-med方法，虽然二者都能分析出内圈故障频率，但tvd-med方法不仅能提取转频29.80hz、二倍转频
58.50hz和故障频率158.20hz，而且在高频范围也能清晰提取多倍故障频率，据此可判断轴承的内圈故障类型。小波阈值降噪方法所得到的包络谱图相较于本发明所提方法并不能快速准确地识别故障特征频率。同时，小波阈值降噪方法在高频范围提取的内圈故障特征频率也不明显，分析效果欠佳。通过对比，说明本发明方法能够准确诊断出滚动轴承外圈故障。
[0161]
另一个实施例对同类型轴承的外圈故障进行诊断，轴承参数如表1所示。经式(20)计算得轴承外圈故障特征频率为fo＝104.44hz，转频为fr＝29.17hz。
[0162]
同理，依据内圈故障特征提取相同的方式，对轴承外圈故障信号进行分析，轴承外圈故障信号的时域波形、频谱如图11a、11b所示。同样采用小波阈值方法作为对比，如图12所示。参数a对加权指标的影响如图13所示，当a＝3.1时，加权指标tc值到达最大值。首先对外圈故障信号通过tvd降噪，然后采用med增强故障信号周期性的冲击特征，结果如图14a、14b所示，经过med处理过的信号冲击特性更加清晰明显。最后通过包络解调得到包络谱，如图15a、15b所示，从图15a、15b中可以看出采用med后故障信号的转频29.29hz、二倍转频59.59hz和故障特征频率105.46hz与理论值较为接近，另外在高频范围中故障特征频率也较为明显，说明本发明方法能够准确诊断出滚动轴承外圈故障。
技术特征：1.基于tvd和med的滚动轴承故障诊断方法，其特征在于，包括以下步骤：步骤1：获取滚动轴承振动信号；步骤2：计算滚动轴承振动信号的峰值因子，初步判断轴承是否产生故障；步骤3：对原始滚动轴承故障振动信号计算一阶总变差；步骤4：利用相关系数和teager能量算子的加权指标选取最优的正则化系数a，并结合一阶总变差得到降噪信号；步骤5：采用最小熵解卷积增强降噪信号中的故障信号，减少脉冲响应函数的传播；步骤6：最后对其包络解调得到包络谱图，提取清晰的故障特征频率及其倍频频率，与计算所得的滚动轴承故障特征频率理论值进行对比，诊断得到滚动轴承故障类型。2.根据权利要求1所述的滚动轴承故障诊断方法，其特征在于，步骤2中，滚动轴承振动信号的峰值因子的计算式如下：x
p
＝max{|x
i
|}
ꢀꢀꢀꢀ
(1)(1)式中x
p
表示信号x(t)的峰值，x
rms
表示信号x(t)的均方根值，c表示信号x(t)的峰值因子，x
i
表示第i个采样点，n表示x的长度。3.根据权利要求2所述的滚动轴承故障诊断方法，其特征在于，步骤3包括以下子步骤：步骤3.1：定义原始信号的一阶总变差，tv＝||dx||1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)式中tv表示一阶总变差；x表示原始振动信号；d表示n-1行n列的矩阵，n表示x的长度；||
·
||1表示l1范数；其中，矩阵d定义如下：步骤3.2：将步骤3.1得到的一阶总变差归结为目标优化问题式中tvd(y,a)表示总变差降噪函数；a为正则化参数；y为含有噪声的信号，x为原始振动信号，||
·
||1表示l1范数，||
·
||2表示l2范数；步骤3.3：利用mm算法(majorization-minimization)算法求解式(6)得到如下结果：式中x
i
、x
i+1
分别表示原始信号、降噪信号；diag( )表示对角矩阵。
4.根据权利要求3所述的滚动轴承故障诊断方法，其特征在于，步骤4包括以下子步骤：步骤4.1：计算原始信号的相关系数，其公式如下：式中c表示互相关系数；分别表示x、y的均值；σ
x
、σ
y
分别表示x、y的标准差，e表示期望；步骤4.2：计算teager能量算子，其公式如下：t(x(n))＝x(n)
2-x(n-1)x(n+1)
ꢀꢀꢀꢀ
(9)t()表示teager能量算子；x(n-1)、x(n)、x(n+1)分别表示连续的三个采样点；步骤4.3：为使降噪信号同时具备相关性和冲击特性，将相关性指标和冲击特性指标结合得到加权指标，如下：t
c
＝t
×
c
r
ꢀꢀꢀꢀ
(10)式中t
c
表示加权指标；t表示teager能量算子；r为调节参数。5.根据权利要求4所述的滚动轴承故障诊断方法，其特征在于，步骤5中，最小熵解卷积的计算过程如下：1)将振动信号y(n)表示如下：y(n)＝h(n)
·
x(n)+e(n)
ꢀꢀꢀꢀ
(11)式中h(n)为传递函数；x(n)为故障信号；e(n)为噪声信号；2)通过寻找最优逆滤波器w(n)来消除卷积影响，得到与故障信号x(n)相近的逆卷积信号i(n)，计算式如下：号i(n)，计算式如下：式中w(i)表示第i个逆滤波器；y(n-i)表示振动信号；l表示滤波器的个数；3)通过计算逆卷积信号i(n)的范数衡量其熵的大小，并将其作为目标函数获得最优输出：式中o2( )表示逆卷积信号的熵值；i2(i)表示逆卷积信号幅值的平方；4)为使得式(14)范数最大，令其导数为0，5)结合式(12)得到如下结果，并采用矩阵形式表示为：c＝wa
ꢀꢀꢀꢀ
(17)
w＝ca-1
ꢀꢀꢀꢀ
(18)式中w为逆滤波器的系数矩阵；c为输出信号与输入信号的互相关矩阵；a为输入信号自相关矩阵。6.根据权利要求1-5任意一项所述的滚动轴承故障诊断方法，其特征在于，步骤6中，滚动轴承内圈故障特征频率的计算式如下：式中f0表示内圈特征频率的理论值，z代表滚动体个数，d表示滚动体直径，d表示节圆直径，α表示轴承接触角，f
r
表示转频。7.根据权利要求6所述的滚动轴承故障诊断方法，其特征在于，滚动轴承外圈故障特征频率的计算式如下：式中f1表示外圈特征频率的理论值。8.根据权利要求6所述的滚动轴承故障诊断方法，其特征在于，转频fr的计算式如下：式中n
i
表示转速。
技术总结本发明涉及基于TVD和MED的滚动轴承故障诊断方法，包括：获取滚动轴承振动信号；采用峰值系数法初步判断滚动轴承是否产生故障；对原始滚动轴承故障信号进行总变差降噪，利用相关系数和Teager能量算子的加权指标确定正则化参数a的最优值；对降噪信号采用最小熵解卷积滤波，减少脉冲响应函数的传播；最后对增强信号进行包络解调，提取清晰的故障特征频率及其倍频，诊断轴承故障类型。本发明解决了包络解调不能快速识别高频范围内故障特征频率的问题，提出的TVD-MED复合方法，对高频范围倍频故障频率具有更好的提取效果；本发明在对原始轴承振动信号的降噪过程中，同时考虑相关性和冲击特性，利用加权指标选取最优的正则化参数，使得降噪效果更好。使得降噪效果更好。使得降噪效果更好。
技术研发人员：王林军 邹腾枭 徐洲常 蔡康林 刘洋
受保护的技术使用者：三峡大学
技术研发日：2022.07.03
技术公布日：2022/11/1
