一种基于比例频带选择准则的轴承故障诊断方法
技术领域
本发明涉及一种基于比例频带选择准则的轴承故障诊断方法，属于故障诊断技术与信号处理分析技术领域。
背景技术
滚动轴承作为旋转机械的支撑部件，其健康程度直接影响旋转机械的精度和寿命。因此，轴承故障诊断和状态监测在实际应用中尤为重要。
随着故障诊断技术的不断应用于工程实际，近年来，许多轴承故障诊断的快速算法逐渐被提出，比如IESFOgram算法，该算法使用短时傅里叶变换(STFT)作为快速有效的估计器。然而，为有效降低算法的计算成本和包含更多的感兴趣分量，常常设置较小的窗宽Nw和较高的采样频率fs，造成谱频率f分辨率往往较差。另外，IESFOgram算法采用1/3-二叉树结构划分奈奎斯特频带，研究发现，在频率分量f分辨率较差的条件下无法实现奈奎斯特频带的严格、合理划分，即频带分量f数量无法整除频带数量。此外，当共振频率fn在fs/2、fs/4、fs/8、fs/16等为中心的频带范围内时，因其所确定的解调频带属于单边频带，导致IESFOgram算法无法确定干扰较小的解调频带。
另一方面，轴承早期故障冲击相对较弱，常常淹没于背景噪声和齿轮啮合等振动分量中，其故障特征提取较为困难，较强的噪声将对基于边带能量比SER指标确定优化解调频带造成干扰。此外，轴承故障理论特征频率与实际特征频率往往不尽相同，主要由以下方面组成：①由于轴承保持架限制了滚动体的运动，滚道和滚动体之间存在1％～2％的随机滑动，导致轴承故障振动冲击不具有严格的周期特性，即为二阶循环平稳信号；②电机实际运行中存在微小的转速波动，直接导致了轴承各部件的实际特征频率与理论特征频率存在差异。以上原因导致IESFOgram无法有效确定包含轴承故障信息最丰富，且干扰更小的解调频带，即所确定的优化解调频带仍然无法有效揭示轴承故障特征。
发明内容
为解决上述问题，本发明提供了一种基于比例频带选择准则的轴承故障诊断方法，解决了传统IESFOgram算法在较强背景噪声、轴承故障特征频率实际值与理值论存在差异和1/3-二叉树结构划分解调频带的条件下，无法确定包含丰富轴承故障信息且干扰较小解调频带的问题。
本发明基于比例频带选择准则的轴承故障诊断方法在采用短时傅里叶变换(STFT)作为估计器的基础上，基于比例频带选择准则，结合MMAD阈值消噪技术和设置特征频率ft，有效提升IESFOgram算法确定包含轴承故障信息最丰富的解调频带的鲁棒性，实现轴承故障特征的有效提取。
本发明基于比例频带选择准则的轴承故障诊断方法如下：
步骤1：通过加速度传感器获得包含轴承故障信息的原始振动信号；
步骤1中通过加速度传感器获取包含轴承故障信息的原始振动信号，具体是将加速度传感器固定安装在测试对象外壳上，通过信号采集系统获取包含轴承故障信息的振动信号；
步骤2：使用短时傅里叶变换(STFT)将所采集的原始振动信号划分为相同长度的数据块，计算各子数据块的频谱；其计算式如下
式中XSTFT(i，fζ)表示离散信号x(tn)在窗长为Nw范围内的傅里叶变换，时间间隔tn＝n/fs，fs表示采样频率，n＝0,1,…,L-1，L表示Nw的长度，w(m)表示窗函数，离散频率fζ＝ζ△f，其中频率分辨率△f＝fs/Nw，i＝1,…,L/R+1，ζ＝0,1,…,Nw-1，移动点数R＝0.25·Nw；
步骤3：计算相邻数据块频谱的谱相关性，提取出与轴承故障相关的分量，获得谱频率f和循环频率α的双变量谱图；
3-1、在区间[(p-1)△f；(p+1)△f]内计算各数据块谱相关性，其计算式为：
式中DFT表示短时傅里叶变换，K＝(L-Nw+R)/R，p＝1,2,…,Nw/2，α表示循环频率，时间指标N0＝(Nw+1)/2，fs表示采样频率，i＝1,2,…,K-1，
3-2、计算整个奈奎斯特频带数据块的循环谱相关性，获得谱频率f和循环频率α的双变量谱图，其计算式如下：
式中α表示循环频率，△α＝fs/L，f表示谱频率。
步骤4：基于移动中位数绝对偏差(MMAD)的阈值消噪技术，对谱频率f中各频率分量所对应的一维谱函数IES的各阶次故障理论特征频率αfault为中心的边频带消噪；
4-1、获得频率分辨率为△f＝fs/Nw的频率分量f＝{f1,f2,…,fζ}，沿谱频率f轴积分获得一维谱函数(即各频率分量均对应一个一维谱函数(IES))，其计算式为：
式中f1和f2表示各子频带的上下截止频率，|Sx(α,f)|表示绝对值操作；
4-2、计算谱频率f中各频率分量所对应的一维谱函数的各阶次故障理论特征频率αfault为中心边频带的阈值，其计算式为
MMAD(IES)＝med{tdem[IES(η)]γ-med(tdem[IES(η)])}
式中med表示取中位数操作，tdem表示串接操作，信号谐波阶次边频带η＝[q·αfault-fb,q·αfault+fb]，其中q＝1,2…,N，αfault表示理论特征频率，fb表示积分频带宽度，fb取值为1～2倍fr，fr表示轴承转频，γ取值范围为[q·αfault-fb,q·αfault+fb]；
4-3、为有效消除不同子频带所对应噪声强度不同对SER指标的影响，使用软阈值消噪技术进行消噪，改进提高包络谱(IIES)的计算式如下
式中σ表示一个极小数；
步骤5：设置频率容差ft，计算消噪后的改进一维谱函数(IIES)的各阶次边带能量比(SER)，其计算式为
式中fc表示中心频率，bw表示频带宽度，ft表示频率容差，谐波阶次q＝1,2,…,N，N表示正整数，SER表示信号在特定频带范围内的能量比，可用于描述信号的信噪比，即SER值越大所对应的解调频带包含轴承故障信息越丰富，ft的取值依据如下:
(1)当电机微小转速波动可忽略时，为包含轴承的最大滑移，ft取值为0.02·αfault。
(2)当电机微小转速波动不可忽略时，即轴承特征频率实际值与理论值差异由转速波动和轴承随机滑动共同作用结果，ft取值为0.02·αfault+Φ，其中Φ表示电机微小转速波动导致的轴承故障部件实际特征频率与理论特征频率的差值。
步骤6：基于比例频带选择准则，划分奈奎斯特频带，取值范围为[0,fs/2]，选择SER值最大时所对应的子频带参数
6-1、基于步骤5获得的各频率分量f＝{f1,f2,…,fζ}所对应的SER值，将SER值降序排列，从大到小选择比例为ε的SER值所对应的频率分量fj，j∈{1,2,…,ζ}；ε比例的取值依据为：
(1)若信号振动信号的信噪比SNR≥0，需要快速收敛，ε≤30％；
(2)若信号振动信号的信噪比SNR<0，需要确定更细化的频带，30％<ε≤80％；
6-2、对频率分量fj分别向左和向右拓展k个频率分量，得到左频带Hfleft[fj-k,fj]和右频带Hfright[fj,fj+k]，对合并后频带中各子频带沿着谱频率f积分获得IIES，并计算各子频带的SER，至此完成第一次频带拓展；每个等级包括子频带个数为δl·Nw/2，其中等级l＝0,1,2,3,…,N-1，δ<1，各子频带的频带宽度bw为(l+1)k·fs/Nw；其中k取值可根据分析信号的信噪比(SNR)来确定，若SNR>0时，k>3·△f，若SNR<0时，△f≤k≤3·△f；
6-3、将步骤6-2的SER值降序排列，并从大到小选择比例为ε的SER值所对应的子频带，重复步骤6-2的操作进行下一次频带划分；若选择出的子频带个数小于2，即j<2，或者拓展后子频带SER值小于拓展前SER值，则停止解调频带划分；
步骤7、获得SER值最大时所对应的子频带参数，即子频带的上下截止频率f1和f2，并沿着谱频率f积分，获得一维谱函数IES，实现滚动故障冲击特征的有效辨识。
7-1、获得SER值最大时子频带所对应的中心频率fc和频带宽度bw，其计算式为
OB＝arg max(SER(fc,bw))
式中arg max表示取最大时操作，f1＝fc-bw/2，f2＝fc+bw/2。
7-2、获得子频带所对应的上下截止频率f1和f2后，沿着谱频率轴f积分，获得一维谱函数IES，实现滚动故障冲击特征的有效提取，一维谱函数IES表达式为
式中f1和f2分别表示各子频带的上下截止频率，|CSCoh(α,f)|表示绝对值操作。
本发明的有益效果是：
(1)本发明提出一种比例频带选择准则，在短时傅里叶变换为估计器的条件下，可实现在低采样频率和较小窗宽的条件下更加合理、严格划分奈奎斯特频带，可选择出干扰较小和包含故障信息最丰富的解调频带；
(2)本发明提出MMAD的阈值消噪技术，对轴承故障特征频率及各谐波阶次的边频带进行消噪操作，有效抑制背景噪声对基于SER指标确定包含轴承故障信息最丰富解调频带的影响；
(3)本发明通过设置频率容差ft，可有效提升SER指标评价所划分解调频带的轴承故障信息的丰富程度的鲁棒性；
(4)本文发明通过对包含轴承故障特征最丰富的解调频带进行包络分析，实现滚动轴承故障特征提取。
附图说明
图1为本发明的频带划分策略图；
图2为本发明中轴承故障诊断流程图；
图3为仿真信号的时域波形图；
图4为对仿真信号分析得到的双变量谱图；
图5为本发明所提方法所优选的解调频带图；
图6为本发明获得的解调频带参数所对应的频谱图；
图7为传统IESFOgram方法得到的解调频带图；
图8为传统IESFOgram获得的解调频带参数所对应的频谱图；
图9为西储大学数据的时域波形图；
图10对西储大学信号分析得到的双变量谱图；
图11本发明所提方法所优选的解调频带图；
图12为本发明获得的解调频带参数所对应的频谱图；
图13传统IESFOgram方法得到的解调频带图；
图14为传统IESFOgram获得的解调频带参数所对应的频谱图。
具体实施方式
下面结合本发明实施例中的附图，对本文发明实例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员没有做出创造性劳动前提下所获得的所有其他实施例，都属于本文发明保护的范围。
在下面的描述中阐述了很多具体细节以便充分理解本发明，但是本发明还可以采用其他不同于在此描述的其他方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。
实施例1：本实例基于比例频带选择准则的轴承故障诊断方法，用于仿真轴承滚动轴承内圈故障信号的诊断，其技术路线如图2所示，具体过程包括：
仿真分析采用的轴承内圈故障信号模型计算式为：
s(t)＝eBt cos(2πfnt)
A＝cos(2πfAt)+randn(t)
式中T表示故障冲击发生的周期，T取值为0.0125，s(t)表示点蚀故障产生的冲击振荡，Ai为第i次冲击幅值，由于滚动轴承故障存在1％～2％的随机滑动，Ti表示第i次冲击相对于周期T的微小波动，为有效模拟滚动轴承的工作环境，模型中考虑噪声n(t)的干扰，randn(t)表示为零均值的随机信号，fA表示调制频率；仿真实验设定采样频率fs＝51.2kHz，轴承故障理论特征频率αfault＝80Hz，固有频率fn＝12.8kHz，调制频率fA＝3Hz，振动信号的信噪比SNR＝-14，电机微小转速波动导致的轴承故障部件实际特征频率与理论特征频率的差值Φ为0.01·αfault，滚动轴承最大随机滑动为0.02·αfault，ft＝0.03·αfault，即ft＝0.03×80＝2.4Hz。
步骤1：仿真滚动轴承内圈故障信号x(t)，t＝1/51200,2/51200,…,5，其时域波形如图3所示。
步骤2：采样频率fs＝51.2kHz，设置循环频率α＝240Hz，窗宽Nw＝512，移动点数R＝0.25·Nw＝128，信号长度L＝2.56×105，K＝1997，△f＝100，fζ＝ζ·△f，其中ζ＝0,1,…,Nw-1；i＝1,…,1984，m＝0,1,…,511，短时傅里叶变换系数计算式为：
步骤3：基于参数K＝1997，i＝1,2…,1996，p＝1,2,…,256，fζ＝0,100,…,51100，N0＝256，fs＝51.2kHz，计算相邻数据块的谱相关性，提取出与轴承故障分量，获得谱频率f和循环频率α的双变量谱图，其计算如下：
式中DFT表示短时傅里叶变换。
进一步地，计算整个奈奎斯特频带内数据块的循环谱相关性，获得谱频率f和循环频率α的双变量谱图如图4所示，其计算如下
步骤4：基于△f＝100，即计算频率分量f＝0,100,…,25600所对应的一维谱函数，其计算式如下
进一步地，基于αfault＝80Hz，q＝1,2…,4，fb＝fA＝3，η＝γ＝[q·αfault-fb,q·αfault+fb]，对一维谱函数IES中各阶次故障理论特征频率αfault为中心的边频带的阈值，其计算如下
MMAD(IES)＝med{tdem[IES(η)]γ-med(tdem[IES(η)])}
进一步地，为有效消除不同子频带所对应噪声强度不同对SER指标的影响，基于σ＝0.001，使用软阈值消噪技术进行消噪，改进提高包络谱(IIES)的计算式如下
步骤5：结合频率容差ft＝2.4Hz，计算消噪后的一维谱函数各阶次边带能量比(SER值)，如图5所示，计算式为
步骤6：基于k＝2·△f＝200，选择出ε＝20％比例的SER值较大时所对应的频率分量，基于如图1所示的频带划分结构，选择SER值最大时所对应的子频带，即中心频率fc＝12.5kHz和频带宽度bw＝800Hz。
步骤7：基于f1＝fc-bw/2＝12.1kHz，f1＝fc+bw/2＝12.9kHz，计算IES，其结果如图6所示
从图6可以发现，轴承内圈故障实际特征频率79.9Hz，二倍频和三倍频谱线分别为159.9Hz和240Hz可清晰辨识。轴承特征频率实际值(79.9Hz)与理论值(80Hz)的误差为0.13％。通过上述分析，本发明所提技术在轴承故障特征频率理论值与实际值存在差异和较强背景噪声的工况下，可有效确定包含轴承故障分量的解调频带，实现轴承内圈故障的准确诊断。
步骤8：为进一步验证本发明所提方法的有效性，采用传统IESFOgram算法对上述数据进行分析，其中采样频率fs＝51.2kHz，循环频率α＝600Hz，窗宽bw＝512，特征频率αfault为80Hz，fb＝fA＝3Hz。所确定的子频带参数如图7所示，其中，中心频率fc＝13.25kHz，频带宽度bw为1300Hz，下截止频率f1＝12.6kHz，上截止频率f2＝13.9kHz。基于截止频率f1和f2，计算IES的结果如图8所示。可以发现，与轴承故障相关的特征谱线无法有效辨识，轴承故障特征提取失败。
对比本发明所提方法和传统IESFOgram算法结果可以发现，传统IESFOgram算法在较强背景噪声、轴承故障特征频率实际值与理值论存在差异和1/3-二叉树结构划分解调频带的条件下，无法有效确定包含轴承故障信息最丰富和干扰更少的解调频带。
实施例2：本实例描述的是本发明方法用于公开轴承外圈故障特征提取
本实施例中对美国西储大学公开数据进行验证，电机输入转速为1750r/min，转频fr＝29.17Hz，信号的采样频率fs为24.8kHz，振动信号的信噪比为SNR＝-6dB，采用SKF6205型号的深沟球轴承作为测试对象，其中，滚子直径d为7.94mm，节圆直径D为39mm，滚子数目n为9，接触角β为0，为模拟轴承外圈故障，在滚动轴承外圈上用线切割方法加工一宽度约为0.53mm，深度约为0.28mm的小槽。
式中fr表示转频，由上式计算可得轴承外圈故障特征频率αfault为104.3Hz。
步骤1：仿真滚动轴承内圈故障信号x(t)，t＝1/24800,2/24800,…,5，其时域波形如图9所示；
步骤2：采样频率fs＝24.8kHz，设置循环频率α＝600Hz，窗宽Nw＝512，移动点数R＝0.25·Nw＝128，信号长度L＝2.56×105，i＝1,…,1984，m＝0,1,…,511，K＝1997，△f＝48.4，fζ＝ζ△f，其中ζ＝0,1,…,511。获得短时傅里叶变换系数，计算如下
步骤3：基于参数K＝1997，i＝1,2,…,1996，p＝1,2,…,256，fζ＝0,48.4,…,247324，N0＝256，计算相邻数据块的谱相关性，提取出与轴承故障分量，获得谱频率f和循环频率α的双变量谱图，其计算如下
式中DFT表示短时傅里叶变换；
进一步地，计算整个奈奎斯特频带内数据块的循环谱相关性，获得谱频率f和循环频率α的双变量谱图如图10所示，其计算式为
步骤4：基于△f＝48.4，即计算频率分量fj＝0,48.4,…,12800所对应的一维谱函数，其计算式如下
进一步地，基于αfault＝104.3Hz，q＝1,2…,4，fb＝fA＝29.17，η＝γ＝[q·αfault-fb,q·αfault+fb]，对一维谱函数IES中各阶次故障理论特征频率αfault为中心的边频带的阈值，其计算如下
MMAD(IES)＝med{tdem[IES(η)]γ-med(tdem[IES(η)])}
进一步地，为有效消除不同子频带所对应噪声强度不同对SER指标的影响，基于σ＝0.001，使用软阈值消噪技术进行消噪，改进提高包络谱(IIES)的计算式如下
步骤5：结合频率容差ft＝αfault·0.02+Φ＝104.3·0.02+0.3＝2.4Hz，计算消噪后各频率分量f的SER值，如图11所示，计算式为
步骤6：基于k＝1·△f＝48.4，选择出ε＝20％比例的SER值较大时所对应的频率分量，基于如图1所示的频带划分结构，选择SER值最大时所对应的子频带，即中心频率fc＝1968.5Hz和频带宽度bw＝375Hz。
步骤7：基于f1＝fc-bw/2＝1781Hz，f1＝fc+bw/2＝2156Hz，计算IES，其结果如图12所示
从图12可以发现，实际轴承内圈故障特征频率104.6Hz，二倍频和三倍频谱线分别为209Hz和313.7Hz可清晰辨识。轴承特征频率实际值(104.6Hz)与理论值(104.3Hz)的误差为0.03％。通过上述分析，本发明所提技术在轴承故障特征频率理论值与实际值存在差异、较强背景噪声的工况下，基于比例选择频带，可有效确定包含轴承故障分量，且干扰较小的解调频带，实现轴承外圈故障的准确诊断。
步骤8：为进一步验证本发明所提方法的有效性，采用传统IESFOgram算法对上述数据进行分析，其中采样频率fs＝24.8kHz，循环频率α＝600Hz，窗宽bw＝512，特征频率αfault为104.3Hz，fb＝fA＝29.17Hz。所确定的子频带参数如图13所示，其中，中心频率fc＝1640Hz，频带宽度bw为280Hz，下截止频率f1＝1500Hz，上截止频率f2＝1781Hz。基于截止频率f1和f2，计算IES的结果如图14所示。可以发现，在频谱图中无法有效辨识与轴承内圈相关的特征谱线无法有效辨识，即传统IESFOgram算法无法有效确定包含轴承故障信息最丰富的解调频带。
本文中应用了具体的实例对本发明的原理及实施方式进行了参数，以上实例的说明只是用于帮助理解本发明及核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变指出，综上，本说明书内容不应理解为本发明的限制。
