买专利卖专利找龙图腾，真高效！ 查专利查商标用IPTOP,全免费！专利年费监控用IP管家,真方便！

申请/专利权人：潍坊学院

摘要：本发明公开了一种滚动轴承的NMD、谱峭度和平滑迭代包络分析方法，该方法首先利用非线性模式分解方法对原始信号进行分解，然后利用数据的重排和替代操作排除分解结果中的噪声分量和趋势项，接着再采用谱峭度方法对第一次滤波后的信号进行分析，得到最优滤波器的中心频率和带宽，然后利用该滤波器对第一次滤波后的信号再进行第二次滤波，然后采用平滑迭代包络分析方法对第二次滤波后的信号进行包络分析，最后根据包络谱确定滚动轴承的故障类型。本发明适合于处理复杂的滚动轴承故障信号，能够准确地判定出滚动轴承的故障类型，具有良好的抗噪性和鲁棒性，便于工程应用。

主权项：一种滚动轴承的NMD、谱峭度和平滑迭代包络分析方法，其特征在于，包括以下步骤：步骤1：利用加速度传感器以采样频率fs测取滚动轴承的振动信号x（k）, （k=1, 2, …,N），N为采样信号的长度；步骤2：采用非线性模式分解Nonlinear Mode Decomposition，NMD算法将信号x（k）分解成n个分量之和，即，其中，ci（k）代表由非线性模式分解算法得到的第i个分量，非线性模式分解已公知，见文献Dmytro Iatsenko, Peter V. E. McClintock, Aneta Stefanovska. Nonlinear mode decomposition: A noise‑robust, adaptive decomposition method[J]. PHYSICAL REVIEW E, 2015, 92: 032916；步骤3：对ci（k）执行重排操作和替代操作，经重排操作得到的数据用cishuffle（k）表示，替代操作后得到数据用ciFTran（k）表示；步骤4：对ci（k）、cishuffle（k）和ciFTran（k）分别执行多重分形去趋势波动分析（Multifractal Detrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，ci（k）的广义Hurst指数曲线用Hi（q）表示；cishuffle（k）的广义Hurst指数曲线用Hishuffle（q）表示；ciFTran（k）的广义Hurst指数曲线用HiFTran（q）表示；步骤5：如果Hi（q） 与Hishuffle（q）或Hi（q） 与HiFTran（q）之间的相对误差小于5%，或者Hi（q） 、Hishuffle（q） 和HiFTran（q）三者都不随q而变化，则抛弃对应的ci（k）分量；步骤6：对剩余的ci（k）分量求和，将该和记为信号经重排和替代滤波后的结果xf1（k）；步骤7：对xf1（k）执行谱峭度分析，求出信号峭度最大处所对应的中心频率f0和带宽B；步骤8： 根据中心频率f0和带宽B对xf1（k）进行带通滤波，得到xf2（k）；步骤9：对信号xf2（k）执行平滑迭代包络分析，得到信号包络eov（k）；步骤10：对得到的信号包络eov（k）执行离散傅里叶变换得到包络谱，根据包络谱特征频率判断机器的故障类型。

全文数据：一种滚动轴承的NMD、谱峭度和平滑迭代包络分析方法技术领域[0001]本发明涉及旋转机械状态监测与故障诊断领域，具体涉及一种滚动轴承的NMD、谱峭度和平滑迭代包络分析方法。背景技术[0002]包络分析技术广泛应用于齿轮和滚动轴承的故障诊断中。现有的包络分析技术有下面三个缺陷：①现有的包络分析技术或者是直接对原始信号进行分析，或者是仅对原始信号进行简单的滤波后再进行分析，因此现有的方法容易受到噪声、趋势及其它成分的干扰，从而导致现有技术的分析精度较低;②现有的包络分析技术是以Hilbert变换为基础，而Hilbert变换要求被分析的信号必须是单分量的窄带信号，否则信号的频率调制部分将要污染信号的幅值包络分析结果，但是目前待分析的信号都不严格满足单分量且窄带的条件，这样就会导致现有技术因精度不高而容易出现误判问题;③由传统方法得到的包络谱存在着端点效应。发明内容[0003]本发明要解决的问题是针对以上不足，提出一种滚动轴承的NMD、谱峭度和平滑迭代包络分析方法，采用本发明的包络分析方法后，具有分析结果准确度和精确度高，并能准确地检测出滚动轴承故障类型的优点。[0004]为解决以上技术问题，本发明采取的技术方案如下：一种滚动轴承的匪D、谱峭度和平滑迭代包络分析方法，其特征在于，包括以下步骤：[0005]步骤1:利用加速度传感器以采样频率fs测取滚动轴承的振动信号xk，（k=l，2,…，吣』为采样信号的长度；[0006]步骤2:采用非线性模式分解NonlinearModeDecomposition，NMD算法将信号Xk分解成η个分量之和，即_，其中，Cik代表由非线性模式分解算法得到的第i个分量，非线性模式分解已公知，见文献[0007]DmytroIatsenko,PeterV.E.McClintock,AnetaStefanovska.Nonlinearmodedecomposition:Anoise-robust,adaptivedecompositionmethod[J].PHYSICALREVIEWE,2015,92:032916；[0008]步骤3:对^k执行重排操作和替代操作，经重排操作得到的数据用Clshufflek表示，替代操作后得到数据用ClFTrank表示；[0009]步骤4:对Cik、Cishufflek和CiFTl:ank分别执行多重分形去趋势波动分析MultifractalDetrendedFluctuationAnalysis,MFDFA，得到广义Hurst指数曲线，Cik的广义Hurst指数曲线用Hiq表示；Cishufflek的广义Hurst指数曲线用Hishuffleq表示;CiFTrank的广义Hurst指数曲线用HiFTranq表示；[0010]步骤5:如果Hiq与Hishuffleq或Hiq与HiFTl:anq之间的相对误差小于5%，或者出⑷、妃_161和H^^q三者都不随q而变化，则抛弃对应的^00分量；[0011]步骤6:对剩余的^10分量求和，将该和记为信号经重排和替代滤波后的结果k;[0012]步骤7:对Xflk执行谱峭度分析，求出信号峭度最大处所对应的中心频率fo和带宽B;[0013]步骤8:根据中心频率f〇和带宽B对Xfik进行带通滤波，得至l」Xf2k;[0014]步骤9:对信号Xf2k执行平滑迭代包络分析，得到信号包络e〇Vk;[0015]步骤10:对得到的信号包络eovk执行离散傅里叶变换得到包络谱，根据包络谱特征频率判断机器的故障类型。[0016]—种优化方案，所述步骤3中数据重排操作包括以下步骤：[0017]随机打乱分量Cik的排列顺序。[0018]进一步地，所述步骤3中数据替代操作包括以下步骤：[0019]1对分量Cik执行离散傅里叶变换，获得分量Cik的相位；[0020]2用一组位于-31,31区间内的伪独立同分布数来代替分量^化的原始相位；[0021]3对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据C1ifftOO，求取数据C1ifftOO的实部。[0022]进一步地，所述步骤4中MFDFA方法包括以下步骤：[0023]1构造xkk=l，2，...，N的轮廓Yi:[0025]Xk代表权利要求1所述步骤4中的Cik或Cishufflek或CiFTrank;[0026]2将信号轮廓Y⑴分成不重叠的Ns段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；[0027]为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2NS段数据；[0028]3利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：[0031]yvi为拟合的第V段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为MF-DFAm;本例中，m=l;[0032]4计算第q阶波动函数的平均值：[0034]5如果Xk存在自相似特征，则第q阶波动函数的平均值Fqs和时间尺度s之间存在舉律关系：[0036]当q=0时，步骤4中的公式发散，这时HO通过下式所定义的对数平均过程来确定：[0038]6对步骤5中的公式两边取对数可得In[FqS]=HqIns+cc为常数），由此可以获得直线的斜率Hq。[0039]进一步地，所述步骤7中的谱峭度方法包括以下步骤：[0040]1构造一个截止频率为fc=0.125+ε的低通滤波器hη;ε〇,本例中fc=〇.3;[0041]2基于hη构造通频带为[0，0.25]的准低通滤波器h〇η和通频带为[0.25，0.5]的准高通滤波器hiη，[0043]3信号C1kη经hoη、Inη滤波并降采样后分解成低频部分C2V1η和高频部分C21+V1η，降采样的因子为2，再经多次迭代滤波后形成滤波器树，第k层有2k个频带，其中Λη表示滤波器树中第k层上的第i个滤波器的输出信号，i=0，…，2k-l，0彡k彡K-I，本例中K=8;Coη代表权利要求1所述步骤7中Xfik;[0044]4分解树中第k层上的第i个滤波器的中心频率fki和带宽Bk分别为[0046]5计算每一个滤波器结果^〇11=0，‘"，21"-1的峭度[0047]6将所有的谱峭度汇总，得到信号总的谱峭度。[0048]进一步地，所述步骤9中的平滑迭代包络分析方法包括以下步骤：[0049]1计算局部均值函数:确定信号Xk所有的局部极值点m，计算相邻两个极值点m和m+i的平均值mi，即[0051]将所有相邻两个极值点的平均值m用折线连接，然后采用移动平均方法进行平滑处理，得到局部均值函数mnk;本例中，移动平均方法中的平滑步长设置为5;在第1次迭代中，Xk代表权利要求1所述步骤9中Xf2k;[0052]2估计信号的包络值:采用局部极值点m计算包络估计值ai[0054]同样，将所有相邻两个包络估计值1用折线连接，然后采用移动平均方法进行平滑处理，得到包络估计函数ank;[0055]3将局部均值函数mnk从原始信号xk中分离出来，得到[0057]4用hnk除以包络估计函数ank从而对hnk进行解调，得到[0059]理想地，S11k是一个纯调频信号，即它的包络估计函数a12k满足a12k=1;如果Snk不满足该条件，则将S11k作为新数据重复以上迭代过程m次，直到得到一个纯调频信号81„11〇，即811111〇满足-1彡811111〇彡1，它的包络估计函数31〇11+：〇1〇满足31〇11+：〇1〇=1，因此有[0061]式中[0063]迭代终止的条件为[0065]在实际应用中，可以设定一个变动量Δ，当满足1-Δ0,本例中fc=0.3;[0120]2基于hη构造通频带为[0，0.25]的准低通滤波器hoη和通频带为[0.25，0.5]的准高通滤波器hiη，[0122]3信号C1kη经h〇η、Iuη滤波并降采样后分解成低频部分C2V1η和高频部分C21+V1η，降采样的因子为2，再经多次迭代滤波后形成滤波器树，第k层有2k个频带，其中Λη表示滤波器树中第k层上的第i个滤波器的输出信号，i=0，…，2k-l，0彡k彡K-I，本例中K=8;Coη代表权利要求1所述步骤7中Xfik;[0123]4分解树中第k层上的第i个滤波器的中心频率fki和带宽Bk分别为[0125]5计算每一个滤波器结果Ληi=0,…，2k_l的峭度[0126]6将所有的谱峭度汇总，得到信号总的谱峭度。[0127]步骤9中的平滑迭代包络分析方法包括以下步骤：[0128]1计算局部均值函数:确定信号xk所有的局部极值点m，计算相邻两个极值点m和rii+i的平均值mi，即[0130]将所有相邻两个极值点的平均值Hi1用折线连接，然后采用移动平均方法进行平滑处理，得到局部均值函数mnk;本例中，移动平均方法中的平滑步长设置为5;在第1次迭代中，Xk代表权利要求1所述步骤9中Xf2k;[0131]2估计信号的包络值:采用局部极值点m计算包络估计值ai[0133]同样，将所有相邻两个包络估计值1用折线连接，然后采用移动平均方法进行平滑处理，得到包络估计函数ank;[0134]3将局部均值函数mnk从原始信号xk中分离出来，得到[0136]4用hnk除以包络估计函数ank从而对hnk进行解调，得到[0138]理想地，S11k是一个纯调频信号，即它的包络估计函数a12k满足a12k=1;如果Snk不满足该条件，则将S11k作为新数据重复以上迭代过程m次，直到得到一个纯调频信号81„11〇，即811111〇满足-1彡811111〇彡1，它的包络估计函数31〇11+：〇1〇满足31〇11+：〇1〇=1，因此有[0142]迭代终止的条件为[0143][0144]在实际应用中，可以设定一个变动量Δ，当满足1-Δaimk1+Δ时，迭代终止;本例中变动量Δ=〇.〇1;[0145]5把迭代过程中产生的所有包络估计函数相乘便可以得到包络信号[0147]试验I，利用具有内圈故障的滚动轴承振动数据对本发明所述算法的性能进行验证。[0148]实验所用轴承为6205-2RSJEMSKF，利用电火花加工方法在轴承内圈上加工深度为0.2794mm、宽度为0.3556mm的凹槽来模拟轴承内圈故障，本实验负载约为0.7457kW，驱动电机转频约为29.5Hz，轴承内圈故障特征频率约为160Hz，采样频率为4.8KHz，信号采样时长为Is。[0149]采集到的内圈故障信号如图4所示。[0150]首先采用传统的包络分析方法对图4所示的信号进行分析，得到的分析结果如图5所示。从图5可以看出，轴承的故障特征完全被掩盖，因此传统的包络分析方法不能有效地提取轴承的故障特征;此外，图5所示包络谱的左端点存在着异常高值，这说明由传统方法得到的包络谱存在着端点效应。[0151]采用本发明所提出的方法对图4所示的信号进行分析，得到的分析结果如图6所示。从图6可以看出，160Hz和320Hz所对应的谱线明显高于其它谱线，这两个频率分别对应轴承内圈故障特征频率的1倍频和2倍频，据此可以判断轴承具有内圈故障；从图6可以看出，由本发明得到的包络谱没有端点效应。[0152]经多次实验表明，在负载和故障尺寸深度不变的情况下，本发明能够可靠识别的最小内圈故障尺寸宽度约为0.26mm，而传统方法能够可靠识别的最小内圈故障尺寸宽度约为0.53mm，精度提高50.9%。[0153]试验2,利用具有外圈故障的滚动轴承振动数据对本发明所述算法的性能进行验证。[0154]实验所用轴承为6205-2RSJEMSKF，利用电火花加工方法在轴承外圈上加工深度为0.2794mm、宽度为0.5334mm的凹槽来模拟轴承外圈故障，本实验负载约为2.237kW，驱动电机转频约为28.7Hz，轴承外圈故障特征频率约为103Hz，采样频率为4.8KHz，信号采样时长为Is。[0155]采集到的外圈故障信号如图7所示。[0156]首先采用传统的包络分析方法对图7所示的信号进行分析，得到的分析结果如图8所示。从图8可以看出，轴承的故障特征完全被掩盖，因此传统的包络分析方法不能有效地提取轴承的故障特征;此外，图8所示包络谱的左端点存在着异常高值，这说明由传统方法得到的包络谱存在着端点效应。[0157]采用本发明所提出的方法对图7所示的信号进行分析，得到的分析结果如图9所示。从图9可以看出，103Hz和206Hz所对应的谱线明显高于其它谱线，这两个频率分别对应轴承外圈故障特征频率的1倍频和2倍频，据此可以判断轴承具有外圈故障；从图9可以看出，由本发明得到的包络谱没有端点效应。[0158]经多次实验表明，在负载和故障尺寸深度不变的情况下，本发明能够可靠识别的最小外圈故障尺寸宽度约为0.33mm，而传统方法能够可靠识别的最小外圈故障尺寸宽度约为0.68mm，精度提高51.5%。[0159]根据试验结果，分析后认为：[0160]1传统的包络分析方法直接对原始信号进行包络分析，或者对仅经过简单处理后的原始信号进行包络分析，与传统的包络分析方法不同，本发明首先利用非线性模式分解对原始信号进行分解，然后利用数据的重排和替代操作排除其中的噪声和趋势分量，仅仅保留信号分量中的有用成分，从而避免了噪声和趋势分量对包络分析结果的影响，提高了准确度和精确度。[0161]2传统的包络分析方法以HiIbert变换为基础，而HiIbert变换要求被分析的信号必须是单分量的窄带信号，否则信号的频率调制部分将要污染信号的包络分析结果，但是目前待分析的信号都不严格满足单分量且窄带的条件，这样就会导致现有技术因精度不高而容易出现误判问题，与传统包络分析方法不同，本发明利用平滑迭代包络分析方法将信号包络与频率调制部分完全分离，能够避免频率调制部分对信号包络分析结果的影响，从而提尚包络分析的精度。[0162]3能够准确地检测出旋转机械的故障类型。[0163]4由传统方法得到的包络谱存在端点效应，而由本发明得到的包络谱能够避免端点效应。[0164]5各步骤作用：[0165]第1步:米集振动ί目号；[0166]第2步:将原始信号分解成不同分量和的形式，其中有些分量对应噪声和趋势项，有些分量对应有用信号；[0167]第3~5步:对上述分解得到的信号执行重排操作和替代操作，剔除其中的噪声分量和趋势项，只保留有用信号；[0168]第6步:将剩余的有用信号求和，将该和作为信号经重排和替代滤波后的结果Xflk；[0169]第7步:对滤波后的信号Xflk执行谱峭度分析，求出信号最大峭度处对应的中心频率fQ和带宽B;[0170]第8步:根据中心频率fo和带宽B对Xflk进行带通滤波，得到信号Xf2k;[0171]第9步:计算信号Xf2k的包络eovk;[0172]第10步:对eovk执行离散傅里叶变换得到包络谱，根据包络谱判断轴承的故障类型。[0173]本领域技术人员应该认识到，上述的具体实施方式只是示例性的，是为了使本领域技术人员能够更好的理解本发明内容，不应理解为是对本发明保护范围的限制，只要是根据本发明技术方案所作的改进，均落入本发明的保护范围。

权利要求：1.一种滚动轴承的NMD、谱峭度和平滑迭代包络分析方法，其特征在于，包括以下步骤：步骤1:利用加速度传感器以采样频率fs测取滚动轴承的振动信号xk，k=l，2,…，N，N为采样信号的长度；步骤2:采用非线性模式分解NonlinearModeDecomposition，NMD算法将信号Xk分解成η个分量之和，即.，其中，Cik代表由非线性模式分解算法得到的第i个分量；步骤3:对^k执行重排操作和替代操作，经重排操作得到的数据用Clshufflek表示，替代操作后得到数据用ClFTrank表示；步骤4:对cik、cishuff13k和ciFTl:ank分别执行多重分形去趋势波动分析MultifractalDetrendedFluctuationAnalysis,MFDFA，得到广义Hurst指数曲线，Cik的广义Hurst指数曲线用Hiq表示；Cishufflek的广义Hurst指数曲线用Hishuffleq表示;CiFTrank的广义Hurst指数曲线用HiFTranq表示；步骤5:如果出如与!1产1^1^或出如与!1严1^如之间的相对误差小于5%，或者出q、Hishuffleq和HiFTranq三者都不随q而变化，则抛弃对应的Cik分量；步骤6:对剩余的^10分量求和，将该和记为信号经重排和替代滤波后的结果Xflk;步骤7:对Xflk执行谱峭度分析，求出信号峭度最大处所对应的中心频率fQ和带宽B;步骤8:根据中心频率f〇和带宽B对Xfik进行带通滤波，得到Xf2k;步骤9:对信号Xf2k执行平滑迭代包络分析，得到信号包络eovk;步骤10:对得到的信号包络eovk执行离散傅里叶变换得到包络谱，根据包络谱特征频率判断机器的故障类型。2.根据权利要求1所述的一种滚动轴承的匪D、谱峭度和平滑迭代包络分析方法，其特征在于，所述步骤3中数据重排操作包括以下步骤：随机打乱分量Clk的排列顺序。3.根据权利要求1所述的一种滚动轴承的匪D、谱峭度和平滑迭代包络分析方法，其特征在于:所述步骤3中数据替代操作包括以下步骤：1对分量Cik执行尚散傅里叶变换，获得分量Cik的相位；2用一组位于-31，31区间内的伪独立同分布数来代替分量^k的原始相位；3对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据C1ifftOO，求取数据Climk的实部。4.根据权利要求1所述的一种滚动轴承的匪D、谱峭度和平滑迭代包络分析方法，其特征在于:所述步骤4中MFDFA方法包括以下步骤：1构造Xyk的轮廓Fi，k=l，2,…，N:Xyk代表权利要求1所述步骤4中的C1k或Clshufflek或ClFTrank;2将信号轮廓ri分成不重叠的队段长度为S的数据，由于数据长度^通常不能整除S，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2队段数据；3利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：hi为拟合的第V段数据的趋势，若拟合的多项式趋势为阶，则记该去趋势过程为MF-DFAffl;本例中，m=l;4计算第Q阶波动函数的平均值：5如果Xyk存在自相似特征，则第Q阶波动函数的平均值FgS和时间尺度3之间存在舉律关系：FqS〜Stm;当Q=O时，步骤4中的公式发散，这时HO通过下式所定义的对数平均过程来确定：6对步骤5中的公式两边取对数可得In[FgS]=丑qInS+C，C为常数，由此获得直线的斜率Q。5.根据权利要求1所述的一种滚动轴承的匪D、谱峭度和平滑迭代包络分析方法，其特征在于:所述步骤7中的谱峭度方法包括以下步骤：1构造一个截止频率为fc=〇.125+ε的低通滤波器hη;ε0，本例中fc=0.3;2基于hη构造通频带为[0，0.25]的准低通滤波器hoη和通频带为[0.25，0.5]的准高通滤波器hη，3信号C1kη经h〇η、Iuη滤波并降采样后分解成低频部分C2V1η和高频部分C21+V1η，降采样的因子为2,再经多次迭代滤波后形成滤波器树，第k层有2k个频带，其中C1kη表示滤波器树中第k层上的第i个滤波器的输出信号，i=0，…，2k-l，0彡k彡K-I，本例中K=8;Coη代表权利要求1所述步骤7中Xfik;4分解树中第k层上的第i个滤波器的中心频率fkl和带宽Bk分别为j5计算每一个滤波器结果Ληi=0,…，2k-l的峭度6将所有的谱峭度汇总，得到信号总的谱峭度。6.根据权利要求1所述的一种滚动轴承的匪D、谱峭度和平滑迭代包络分析方法，其特征在于，所述步骤9中的平滑迭代包络分析方法包括以下步骤：1计算局部均值函数:确定信号zk所有的局部极值点m，计算相邻两个极值点m和m+i的平均值mi，即将所有相邻两个极值点的平均值ΠΗ用折线连接，然后采用移动平均方法进行平滑处理，得到局部均值函数mnk;本例中，移动平均方法中的平滑步长设置为5;在第1次迭代中，zGO代表权利要求1所述步骤9中Xf2k;2估计信号的包络值:采用局部极值点m计算包络估计值1同样，将所有相邻两个包络估计值1用折线连接，然后采用移动平均方法进行平滑处理，得到包络估计函数ank;3将局部均值函数mnk从信号zk中分离出来，得到4用hnk除以包络估计函数ank从而对hnk进行解调，得到理想地，Snk是一个纯调频信号，即它的包络估计函数a12k满足a12k=1;如果snk不满足该条件，则将S11k作为新数据重复以上迭代过程m次，直到得到一个纯调频信号811111〇,8卩811111〇满足-1彡811111〇彡1，它的包络估计函数31111+11〇满足31111+11〇=1，因此有式中迭代终止的条件为在实际应用中，设定一个变动量Δ，当满足I-Agaimk1+Δ时，迭代终止;本例中变动量Δ=〇.〇1;5把迭代过程中产生的所有包络估计函数相乘便得到信号包络

百度查询： 潍坊学院 一种滚动轴承的NMD、谱峭度和平滑迭代包络分析方法