CN105910805B - 一种用于转子碰摩故障诊断的小波局部均值分解方法 - Google Patents

一种用于转子碰摩故障诊断的小波局部均值分解方法 Download PDF

Info

Publication number
CN105910805B
CN105910805B CN201610257375.5A CN201610257375A CN105910805B CN 105910805 B CN105910805 B CN 105910805B CN 201610257375 A CN201610257375 A CN 201610257375A CN 105910805 B CN105910805 B CN 105910805B
Authority
CN
China
Prior art keywords
signal
wpf
components
component
frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201610257375.5A
Other languages
English (en)
Other versions
CN105910805A (zh
Inventor
陈中柘
朱圆
傅思程
张志威
吴迎杰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201610257375.5A priority Critical patent/CN105910805B/zh
Publication of CN105910805A publication Critical patent/CN105910805A/zh
Application granted granted Critical
Publication of CN105910805B publication Critical patent/CN105910805B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Abstract

本发明公开了一种用于转子碰摩故障诊断的小波局部均值分解方法，对局部均值分解产生的PF分量进行小波再分解，最终得到一系列的WPF分量和一个残余分量，用于缓解局部均值分解过程中的模态混淆现象，并将此方法用于转子的碰摩故障诊断。本发明能有效抑制模态混淆现象，提取转子故障特征，从而能更好地进行转子的碰摩故障诊断。此外，本发明还具有算法简单，程序运行时间短，自适应性较好等优点。

Description

一种用于转子碰摩故障诊断的小波局部均值分解方法
技术领域
本发明属于旋转机械故障诊断技术领域，具体涉及一种用于转子碰摩故障诊断的小波局部均值分解方法的设计。
背景技术
旋转机械在现在工业生产中得到广泛应用，它主要由转子、定子、轴承、齿轮传动件等构成。转子是旋转机械设备的重要部件，一旦发生故障造成停机，就可能会带来巨大的经济损失。因此，对转子进行状态监测与故障诊断就显得非常重要。
碰摩故障是转子运行过程中常见的一种故障形式。在高转速高效率的要求下，旋转机械中转子与定子之间的间隙越来越小，这就导致转子系统的碰摩现象不断发生，由此就会产生转轴断裂、机器死锁等一系列非常严重的后果。转子发生碰摩故障时，它的振动信号常常会具有某种特点。因此，可以通过传感器获取转子振动信号，然后对振动信号进行分析处理，从振动信号中提取有关转子碰摩故障特征信息，从而对转子进行故障诊断。
转子系统发生碰摩故障时的振动信号通常表现出非平稳、非线性的特征，因此需要采用时频分析法这一重要的非平稳信号分析手段。在传统的傅里叶变换的基础上，人们研究了不少处理非平稳信号的时频分析方法。典型的有：短时傅里叶变换、Wigner-Ville分布、小波变换、经验模态分解等。经验模态分解是近些年兴起的比较热门的时频分析方法，它有很多优点，但存在诸多的缺陷，例如：模态混淆现象严重、本征模态函数判据不清、端点效应明显以及欠包络和过包络现象。针对经验模态分解方法中的存在的不足，人们又提出局部均值分解的方法。相比较前面几种时频分析方法，该方法具有许多优势，且发展前景广大。局部均值分解方法是目前比较新的一种时频分析方法，该方法最大的优点是具有自适应性，它在机械故障诊断领域得到了广泛应用。但是，和经验模态分解一样，局部均值分解的模态混淆现象依然很严重。所谓的模态混淆现象是指，局部均值分解结果中同一个PF(Product function，乘积函数)分量中出现了尺度或频率差异较大的信号，或者同一尺度或频率的信号被分解到多个不同的PF分量当中。局部均值分解中存在的严重的模态混淆现象对信号分解的准确性产生严重的影响，从而影响故障诊断的准确率。因此，对局部均值分解方法进行改进，抑制其模态混淆现象，将会明显提高故障诊断的准确率。
目前针对局部均值分解过程的存在的模态混淆问题的而提出的改进方法较少，主要的有基于辅助噪声分析的总体局部均值分解方法(ELMD)。该方法的主要过程是：任意给定一个信号x0(t)，向待分解的原始信号x0(t)多次添加不同的白噪声，每次添加完噪声之后，对信号进行局部均值分解，最终得到多次不同的分解结果，将分解得到的多个对应的PF分量的平均值作为最终分解的结果。基于辅助噪声分析的总体局部均值分解方法虽然在一定程度上能抑制局部分解过程的发生的端点效应，但它存在一些问题，例如，添加白噪声的次数和所加白噪声的幅值的选择具有很强的主观性，方法自适应性较差。当添加的白噪声的次数和幅值选择比较合理，能抑制高频模态混淆，却人为造成低频的模态混淆(或者能抑制低频模态混淆，却人为造成高频的模态混淆)。此外，基于辅助噪声分析的总体局部均值分解方法还存在一个比较严重的缺陷，那就是该方法的算法较为复杂，程序运行时间很长，这必然会限制该方法在实时性要求比较高的信号处理的领域应用。
发明内容
本发明的目的是为了解决现有技术中局部均值分解方法中存在严重的模态混淆现象，严重地影响了对信号分解的准确性的问题，提出了一种用于转子碰摩故障诊断的小波局部均值分解方法。
本发明的技术方案为：一种用于转子碰摩故障诊断的小波局部均值分解方法，包括以下步骤：
S1、搭建模拟转子碰摩故障实验台，利用位移传感器采集转子碰摩故障的振动信号x(t)；
S2、对采集的振动信号x(t)进行局部均值分解，得到k个PF分量和一个残余分量uk(t)；
S3、采用小波分解对所有的PF分量进行再分解与重构，得到k个WPF分量和一个残余分量vk(t)；
S4、计算每一个WPF分量与初始振动信号x(t)的互相关系数，从所有的WPF分量中选三个连续WPF分量作为重要WPF分量，求出这三个重要WPF分量的包络信号；分别画出三个重要WPF分量的幅值谱和它们所对应的包络信号的幅值谱，并从中找到转子碰摩故障的特征频率。
进一步地，步骤S2包括以下分步骤：
S21、找到振动信号x(t)的所有局部极值点，求出任意两个相邻的局部极值点的平均值，将所有相邻局部极值点的平均值用直线连接，得到局部均值线段，然后采用滑动平均法进行平滑处理，得到局部均值函数m11(t)；
S22、计算任意两个相邻局部极值点见得包络估计值，将所有相邻包络估计值用直线连接，得到包络估计线段，然后采用滑动平均法进行平滑处理，得到包络估计函数a11(t)；
S23、将局部均值函数m11(t)从振动信号x(t)中分离出来，得到分离信号h11(t)；
S24、用分离信号h11(t)除以包络估计函数a11(t)，对h11(t)进行解调，得到解调信号s11(t)；
S25、判断s11(t)是否为纯调频信号，若是则进入步骤S26，否则用s11(t)代替x(t)，重复步骤S21-S24，迭代n次，直到s1n(t)为纯调频信号为止，进入步骤S26；
S26、把迭代过程中产生的所有包络估计函数相乘，得到包络信号a1(t)；
S27、将包络信号a1(t)和纯调频信号s1n(t)相乘，得到x(t)的第一个PF分量PF1(t)；
S28、将PF1(t)从x(t)中分离出来，得到一个新的信号u1(t)，将u1(t)作为原始数据代替x(t)，重复步骤S21-S27，循环k次，直到uk(t)为一个单调函数为止，将初始振动信号x(t)分解为k个PF分量和一个残余分量uk(t)之和。
进一步地，步骤S21中对局部均值线段进行平滑处理时，滑动平均跨度为相邻极值点最大距离的三分之一。
进一步地，步骤S22中对包络估计线段进行平滑处理时，滑动平均跨度为相邻极值点最大距离的三分之一。
进一步地，步骤S25中判断s1n(t)是否为纯调频信号的具体方法为：
计算s1n(t)的包络估计函数a1(n+1)(t)，设置一个增减量δ>0，判断是否满足条件1-δ≤a1(n+1)(t)≤1+δ，若是则说明s1n(t)是纯调频信号，否则s1n(t)不是纯调频信号。
进一步地，δ的取值在0.001～0.1之间。
进一步地，δ的取值为0.05。
进一步地，步骤S3包括以下分步骤：
S31、对PF1(t)进行若干层小波分解，得到一个小波低频分量和若干个小波高频分量；
S32、剔除PF1(t)的小波低频分量，并把它给PF2(t)，PF1(t)剩余的信号为第一个WPF分量，记作WPF1(t)，PF2(t)获得PF1(t)的小波低频分量之后更新得到PF2 *(t)；
S33、用PF2 *(t)替换PF1(t)，重复步骤S31-S32，直到获得PFk(t)的更新分量PFk *(t)；
S34、对PFk *(t)进行若干层小波分解，得到一个小波低频分量和若干个小波高频分量，剔除PFk *(t)的小波低频分量，并把它给残余分量uk(t)，PFk *(t)的剩余信号为第k个WPF分量，记作WPFk(t)，残余分量uk(t)获得PFk *(t)的小波低频分量，得到残余分量vk(t)。
进一步地，步骤S4中选择三个连续WPF分量作为重要WPF分量的原则为：
(1)这三个WPF分量与x(t)的互相关系数之和最小；
(2)这三个WPF分量的模态最高。
本发明的有益效果是：
(1)本发明对局部均值分解产生的PF分量进行小波再分解，最终得到一系列的WPF分量和一个残余分量，能较好的抑制局部均值分解过程中的模态混淆现象，使信号分解的准确性得到明显提高，有利于故障特征的提取。
(2)本发明算法较为简单，程序运行时间短，抑制局部均值分解的模态混淆效果较好，自适应性较好。
附图说明
图1为本发明提供的一种用于转子碰摩故障诊断的小波局部均值分解方法流程图。
图2为本发明实施例模拟转子单点局部碰摩故障实验台实物图。
图3为本发明实施例振动信号x(t)及其幅值谱示意图。
图4为本发明步骤S2的分步骤流程图。
图5为本发明实施例振动信号x(t)局部均值分解的结果及幅值谱示意图。
图6为本发明步骤S3的分步骤流程图。
图7为本发明实施例PF2 *(t)幅值谱示意图。
图8为本发明实施例PF3 *(t)幅值谱示意图。
图9为本发明实施例PF4 *(t)幅值谱示意图。
图10为本发明实施例PF5 *(t)幅值谱示意图。
图11为本发明实施例振动信号小波局部均值分解结果及幅值谱示意图。
图12为本发明实施例的三个重要WPF分量的包络幅值谱示意图。
具体实施方式
下面结合附图对本发明的实施例作进一步的说明。
本发明提供了一种用于转子碰摩故障诊断的小波局部均值分解方法，如图1所示，包括以下步骤：
S1、搭建模拟转子碰摩故障实验台，利用位移传感器采集转子碰摩故障的振动信号x(t)。
如图2所示为模拟转子发生单点局部碰摩故障实验装置。转子的转速由输入电机控制，本发明实施例中，设定转子的转速为3000r/min。位移传感器水平安装，用于测取转子发生单点局部碰摩故障的振动信号。设定信号采样频率为8000hz，采样时间长度选取为1s。采集的振动信号为x(t)，如图3(a)所示。对振动信号x(t)进行快速傅里叶变换，求它的幅值谱，如图3(b)所示。从图3(b)可以看出，振动信号的幅值谱中主要频率成分为转频50Hz及其3倍频，而高频碰摩成分并不明显。
S2、对采集的振动信号x(t)进行局部均值分解，得到k个PF分量和一个残余分量uk(t)。
如图4所示，该步骤具体包括以下分步骤：
S21、找到振动信号x(t)的所有局部极值点(局部极大值点和局部极小值点)，求出任意两个相邻的局部极值点的平均值，即：
式中，ni为第i个极值，mi为第i个极值点和第i+1个极值点的平均值。将所有相邻局部极值点的平均值用直线连接，得到局部均值线段，然后采用滑动平均法进行平滑处理，得到局部均值函数m11(t)。本发明实施例中，对局部均值线段进行平滑处理的滑动平均跨度为相邻极值点最大距离的三分之一。
S22、计算任意两个相邻局部极值点见得包络估计值ai：
同样地，将所有相邻包络估计值用直线连接，得到包络估计线段，然后采用滑动平均法进行平滑处理，得到包络估计函数a11(t)。本发明实施例中，对包络估计线段进行平滑处理的滑动平均跨度为相邻极值点最大距离的三分之一。
S23、将局部均值函数m11(t)从振动信号x(t)中分离出来，得到分离信号h11(t)：
h11(t)＝x(t)-m11(t) (3)
S24、用分离信号h11(t)除以包络估计函数a11(t)，对h11(t)进行解调，得到解调信号s11(t)：
s11(t)＝h11(t)/a11(t) (4)
S25、判断s11(t)是否为纯调频信号，若是则进入步骤S26，否则用s11(t)代替x(t)，重复步骤S21-S24，迭代n次，直到s1n(t)为纯调频信号为止，进入步骤S26。
这里判断s1n(t)是否为纯调频信号的具体方法如下：
计算s1n(t)的包络估计函数a1(n+1)(t)，若a1(n+1)(t)＝1则说明s1n(t)是一个纯调频信号，即有-1≤s1n(t)≤1，所以有：
其中：
但是，-1≤s1n(t)≤1且a1(n+1)(t)＝1是s1n(t)为纯调频信号的理想条件，该条件在实际的迭代分解过程中是无法实现的。为了获得较为理想的纯调频信号，需要给出一个合理的迭代终止条件。在实际的计算过程中，可以设置一个增减量δ>0，当满足条件1-δ≤a1(n+1)(t)≤1+δ，则说明s1n(t)是纯调频信号，否则s1n(t)不是纯调频信号。δ的取值范围需要根据不同的信号和不同的精度要求来设定，通常δ取值越小，计算量越大，局部均值分解的精度越高。根据大量试用经验，δ取值在0.001～0.1比较合理。本发明实施例中δ取值为0.05。
S26、把迭代过程中产生的所有包络估计函数相乘，得到包络信号a1(t)：
a1(t)＝a11(t)*a12(t)*a13(t)*…a1(n+1)(t) (7)
S27、将包络信号a1(t)和纯调频信号s1n(t)相乘，得到x(t)的第一个PF分量PF1(t)：
PF1(t)＝a1(t)*s1n(t) (8)
S28、将PF1(t)从x(t)中分离出来，得到一个新的信号u1(t)，将u1(t)作为原始数据代替x(t)，重复步骤S21-S27，循环k次，直到uk(t)为一个单调函数为止，即：
至此，将初始振动信号x(t)分解为k个PF分量和一个残余分量uk(t)之和，即：
本发明的实施例中，k＝5，即转子碰摩故障振动信号x(t)经过局部均值分解，得到5个PF分量和一个残余分量，分别记作PF1(t)，PF2(t)，PF3(t)，PF4(t)，PF5(t)和u5(t)。振动信号x(t)进行局部均值分解的结果(忽略残余分量)如图5(a)所示。分别对PF1(t)，PF2(t)，PF3(t)，PF4(t)，PF5(t)做快速傅里叶变换，得到各自幅值谱，如图5(b)所示。局部均值分解所得的每一个分量(除残差分量外)视为一个模态，振动信号x(t)分解得到5个模态，且各阶模态的平均频率依次减少。
由图5我们可以知道：PF1(t)分量是由大量高频噪声成分和部分中低频噪声成分组成，当然还有部分转子碰摩故障特征频率成分混在PF1(t)分量中。PF1(t)分量的成分非常复杂，从0～4000Hz都存在，PF1(t)分量显然发生模态混淆；PF2(t)分量主要由转频50Hz，转频的三倍频150Hz成分，中低频噪声以及部分高频局部碰摩故障频率成分组成，显然PF2(t)分量也发生模态混淆；PF3(t)分量主要由转频50Hz成分和一些低频噪声构成；PF4(t)由低频噪声成分构成；PF5(t)是伪PF分量(由局部均值分解的端点效应，信号采样频率不足等引起的)。由局部均值分解的结果我们可以看出，模态1，模态2均发生模态混淆现象，因此整个局部分解的过程产生严重的模态混淆现象。
另外，模态1和模态2产生模态混淆现象会使我们无法提取转子碰摩故障的信息。局部碰摩故障频率成分本身比较微弱，又被分解到两个模态，即模态1和模态2中。模态1中有大量噪声信号，而模态2中有强大的背景信号(转频信号和3倍转频信号)，局部碰摩故障信息淹没在强大的背景信号(转频信号和3倍转频信号)和噪声信号中。因此我们无法获得转子高频碰摩故障的信息。
对信号x(t)采用某种方法进行分解，如果各阶模态都没有发生混淆，那么应该得到这样理想结果：模态1由高频噪声构成，模态2由高频碰摩成分构成，模态3由3倍转频成分构成，模态4由转频成分构成，模态5由低频造成构成，模态6(7，8…)是伪分量。由此可见，如果信号分解过程未发生模态混淆，那么高频碰摩成分就会独占一个模态，而不是被分到两个不同的模态中。这样，我们可以对高频碰摩成分所在模态进行分析，提取故障的特征。
为了能达到上述效果，我们采用小波分解对局部均值分解的结果进行再分解与重构，也就是对所得的PF分量进行修正，使最终能得到一系列频率成分相对简单的分量。
S3、对所得的所有PF分量进行修正，即采用小波分解对所有的PF分量进行再分解与重构，得到k个WPF(wavelet product function component，小波乘积函数)分量和一个残余分量vk(t)。
如图6所示，该步骤具体包括以下分步骤：
S31、对PF1(t)进行若干层小波分解，得到一个小波低频分量和若干个小波高频分量；
S32、剔除PF1(t)的小波低频分量，并把它给PF2(t)，PF1(t)剩余的信号为第一个WPF分量，记作WPF1(t)，PF2(t)获得PF1(t)的小波低频分量之后更新得到PF2 *(t)；
S33、用PF2 *(t)替换PF1(t)，重复步骤S31-S32，直到获得PFk(t)的更新分量PFk *(t)；
S34、对PFk *(t)进行若干层小波分解，得到一个小波低频分量和若干个小波高频分量，剔除PFk *(t)的小波低频分量，并把它给残余分量uk(t)，PFk *(t)的剩余信号为第k个WPF分量，记作WPFk(t)，残余分量uk(t)获得PFk *(t)的小波低频分量，得到残余分量vk(t)。
至此，将x(t)分解为了k个WPF分量和一个残余分量vk(t)之和，即：
本发明实施例中，k＝5，即得到5个PF分量和一个残余分量，分别记作PF1(t)，PF2(t)，PF3(t)，PF4(t)，PF5(t)和u5(t)。
将PF1(t)中的频率低于某一值(截断频率)的成分给“抽取”出来给下一个分量PF2(t)，剩余的信号的频率成分就会简单些。为此我们对PF1(t)进行若干层小波分解(小波分解的层数由设置的截断频率决定的)，小波分解之后得到一个小波低频分量和多个小波高频分量。剔除PF1(t)的小波低频分量，并把它给PF2(t)。PF1(t)剩余的信号就是第一个WPF分量，记作WPF1(t)。PF2(t)获得PF1(t)的小波低频分量之后得更新得到PF2 *(t)。
上述过程涉及到两个关键问题：
(1)设置PF1(t)截断频率fb1。
我们要把PF1(t)分量中的频率低于某一值(截断频率)的成分给抽取出来给下一个分量PF2(t)，PF1(t)的截断频率fb1的设置是依据PF1(t)的幅值谱和局部均值分解表现的滤波组结构特征。
对纯粹白噪声信号进行局部均值分解的滤波器组结构特征表现为频域上的一个高通滤波器和一系列连续的带通滤波器。对x(t)进行局部均值分解得到第一个PF分量的滤波器结构特征表现出来是频域上的一个非理想高通滤波器，它的衰减速率慢，过渡带很宽。从PF1(t)的幅值谱估计这个高通滤波器的截止频率大概800hz左右，所以截断频率fb1选择800hz。由于小波分解表现的滤波器组的特点，在实际中并不能完全把PF1(t)中频率小于截断频率fb1的低频部分完全“拉到”PF2(t)中去。
(2)确定PF1(t)小波分解层数pp1。
对信号进行小波分解首先要确定小波分解的层数。对离散白噪声信号进行小波分解的滤波器组结构特征表现为频域上的一个高通滤波器和一系列连续的带通滤波器，滤波器截止频率依次是fs/2，fs/4，fs/8…，fs为待分解信号的采样频率。小波分解的层数pp可按以下公式计算：
式中，fs为待分解信号的采样频率，fb为待分解信号的截断频率。把fb1＝800带入公式中，得pp1向下取整取2。
PF1(t)的截断频率fb1和小波分解层数pp1确定后，对PF1(t)分量进行2层小波分解，小波基函数选用db7。小波分解之后，得到一个小波低频分量cA12和2个小波高频分量cD12，cD11。剔除PF1(t)的小波低频分量cA12，并把它给PF2(t)。PF1(t)剩余的信号就是第一个WPF分量，记作WPF1(t)。PF2(t)获得PF1(t)的小波低频分量cA12之后得更新得到PF2 *(t)，即：
PF1(t)＝cA12+cD12+cD11 (13)
WPF1(t)＝PF1(t)-cA12 (14)
PF2 *(t)＝PF2(t)+cA12 (15)
同理，将PF2 *(t)中的频率低于某一值(截断频率)的成分给抽取出来给下一个分量PF3(t)，剩余的信号的频率成分就会简单些。为此我们对PF2 *(t)进行若干层小波分解(小波分解的层数由设置的截断频率决定的)，小波分解之后得到一个小波低频分量和多个小波高频分量。剔除PF2 *(t)的小波低频分量，并把它给PF3(t)。PF2 *(t)剩余的信号就是第二个WPF分量，记作WPF2(t)。PF3(t)获得PF2 *(t)的小波低频分量之后得更新得到PF3 *(t)。同样，我们需要设置PF2 *(t)截断频率fb2，确定PF2 *(t)小波分解层数pp2。
PF2 *(t)截断频率fb2由其幅值谱确定的。PF2 *(t)的幅值谱如图7所示，从幅值谱图我们可以看出，PF2 *(t)这个分量有三个频率主要频率成分，分别是50hz，150hz，500hz。设置一个截断频率fb2，将PF2 *(t)的频率小于截断频率fb2的分量“拉到”PF3(t)这个分量中去。显然PF2 *(t)截断频率设定范围为150～500hz。令fb2＝300hz，PF2 *(t)小波分解层数仍然采用公式来计算，把fb2＝300hz带入公式中，得到
pp2向下取整取3。
PF2 *(t)的截断频率fb2和小波分解层数pp2确定后，对PF2 *(t)分量进行3层小波分解，小波基函数选用db7。小波分解之后，得到一个小波低频分量cA23和三个小波高频分量cD23，cD22，cD21。剔除PF2 *(t)的小波低频分量cA23，并把它给PF3(t)。PF2 *(t)剩余的信号就是第二个WPF分量，记作WPF2(t)。PF3(t)获得PF2 *(t)的小波低频分量cA23之后得更新得到PF3 *(t)，即
PF2 *(t)＝cA23+cD23+cD22+cD21 (16)
PF3 *(t)＝PF3(t)+cA23 (18)
同理，将PF3 *(t)中的频率低于某一值(截断频率)的成分给抽取出来给下一个分量PF4(t)，剩余的信号的频率成分就会简单些。为此我们对PF3 *(t)进行若干层小波分解(小波分解的层数由设置的截断频率决定的)，小波分解之后得到一个小波低频分量和多个小波高频分量。剔除PF3 *(t)的小波低频分量，并把它给PF4(t)。PF3 *(t)剩余的信号就是第三个WPF分量，记作WPF3(t)。PF4(t)获得PF3 *(t)的小波低频分量之后得更新得到PF4 *(t)。同样，我们需要设置PF3 *(t)截断频率fb3，确定PF3 *(t)小波分解层数pp3。
PF3 *(t)截断频率fb3由其幅值谱确定的。PF3 *(t)的幅值谱如图8所示，从幅值谱图我们可以看出，PF3 *(t)这个分量主要有2个频率成分，分别是50hz和150hz。设置一个截断频率fb3，将PF3 *(t)的频率小于截断频率fb3的分量“拉到”PF4(t)这个分量中去。显然PF3 *(t)截断频率设定范围为50～150hz。令fb3＝100hz，PF3 *(t)小波分解层数仍然采用公式来计算，把fb3＝100hz带入公式中，求得pp3的值后并向下取整得pp3＝5。
PF3 *(t)的截断频率fb3和小波分解层数pp3确定后，对PF3 *(t)分量进行5层小波分解，小波基函数选用db7。小波分解之后，得到一个小波低频分量cA35和五个小波高频分量cD35，cD34，cD33，cD32，cD31。剔除PF3 *(t)的小波低频分量cA35，并把它给PF4(t)。PF3 *(t)剩余的信号就是第三个WPF分量，记作WPF3(t)。PF4(t)获得PF3 *(t)的小波低频分量cA35之后得更新得到PF4 *(t)，即：
同理，将PF4 *(t)中的频率低于某一值(截断频率)的成分给抽取出来给下一个分量PF5(t)，剩余的信号的频率成分就会简单些。为此我们对PF4 *(t)进行若干层小波分解(小波分解的层数由设置的截断频率决定的)，小波分解之后得到一个小波低频分量和多个小波高频分量。剔除PF4 *(t)的小波低频分量，并把它给PF5(t)。PF4 *(t)剩余的信号就是第四个WPF分量，记作WPF4(t)。PF5(t)获得PF4 *(t)的小波低频分量之后得更新得到PF5 *(t)。同样，我们需要设置PF4 *(t)截断频率fb4，确定PF4 *(t)小波分解层数pp4。
PF4 *(t)截断频率fb4由其幅值谱确定的。PF4 *(t)的幅值谱如图9所示，从幅值谱图我们可以看出，PF4 *(t)这个分量主要有2个频率成分，分别是50hz和10hz。设置一个截断频率fb4，将PF4 *(t)的频率小于截断频率fb4的分量“拉到”PF5(t)这个分量中去。显然PF4 *(t)截断频率设定范围为10～50hz。令fb3＝30hz，PF4 *(t)小波分解层数仍然采用公式来计算，把fb4＝30hz带入公式中，求得pp4的值后并向下取整得pp4＝7。
PF4 *(t)的截断频率fb4和小波分解层数pp4确定后，对PF4 *(t)分量进行7层小波分解，小波基函数选用db7。小波分解之后，得到一个小波低频分量cA47和七个小波高频分量cD47，cD46，cD45，cD44，cD43，cD42，cD41。剔除PF4 *(t)的小波低频分量cA47，并把它给PF5(t)。PF4 *(t)剩余的信号就是第四个WPF分量，记作WPF4(t)。PF5(t)获得PF4 *(t)的小波低频分量cA47之后得更新得到PF5 *(t)，即：
用同样的方法处理PF5 *(t)，就可以得到第五个WPF分量WPF5(t)和最终的残差分量v5(t)。
整个过程可以表述为振动信号x(t)经小波局部均值分解之后，得到5个WPF分量。PF5 *(t)的幅值谱如图10所示，小波局部均值分解的结果及幅值谱如图11所示。
下面我们将对比局部均值分解的结果和小波局部均值分解的结果，来说明小波局部均值分解方法的有效性，为了方便，我们对比它们分解结果前三个分量。
对比局部均值分解的结果及幅值谱(图5)与小波局部均值分解的结果及幅值谱(图11)，我们可以看出，局部均值分解得到的第一个分量PF1(t)是由大量高频噪声成分和部分中低频噪声成分组成，当然还有部分转子碰摩故障特征频率成分混在PF1(t)分量中。而小波局部均值分解得到的第一个分量WPF1(t)由高频噪声和部分中频噪声组成，WPF1(t)比PF1(t)的成分简单些。
局部均值分解得到的第二个分量PF2(t)分量主要由转频50Hz，转频的三倍频150Hz成分，中低频噪声以及部分高频局部碰摩故障频率成分组成，而小波局部均值分解得到的第二个分量WPF2(t)是主要由高频碰摩信号构成和少量噪声，WPF2(t)比PF2(t)的成分简单些，更重要的是转子故障振动信号经过小波局部均值后高频碰摩成分被单独分到一个模态中，而转子故障振动信号经过局部均值后高频碰摩成分被分到两个模态中去。因此，可以直接分析WPF2(t)来获取转子碰摩故障信息。
局部均值分解得到的第三个分量PF3(t)分量主要由转频50Hz成分和一些低频噪声构成，小波局部均值分解得到的第三个分量WPF3(t)分量也是由转频50Hz成分和一些低频噪声构成。
由此可以看出，小波局部分解得到前2个分量(WPF1(t)，WPF2(t))要比局部均值分解得到前2个分量(PF1(t)，PF2(t))好些，前者的模态混淆没有后者严重。小波局部分解能把碰摩成分分离出来，而局部均值分解不能。
S4、计算每一个WPF分量与初始振动信号x(t)的互相关系数，从所有的WPF分量中选三个连续WPF分量作为重要WPF分量，求出这三个重要WPF分量的包络信号；分别画出三个重要WPF分量的幅值谱和它们所对应的包络信号的幅值谱，并从中找到转子碰摩故障的特征频率。
其中，互相关系数的计算公式如下：
式中，x，y为两个维数相同的离散信号。和分别表示信号x与y的均值。
本发明实施例中，计算每个WPF分量与初始振动信号x(t)的互相关系数，得γ1＝0.057，γ2＝0.030，γ3＝0.185，γ4＝0.979，γ5＝0.049。小波局部均值分解所得的每一个模态(WPF1(t)，WPF2(t)，WPF3(t)…)的平均频率是依次从高到低。局部碰摩故障信号主要是幅值调制信号，载波频率是比较高，另外局部碰摩故障的振动信号比较微弱，结合这两个特点可知，局部碰摩故障信息集中在前面某1个或者某2个WPF分量，且局部碰摩故障信息所在的WPF分量与x(t)的互相关系数比较小。为此，我们从所有的WPF分量选择三个连续的WPF分量作为三个重要WPF分量。这三个重要WPF分量会有有一个或者两个包含高频碰摩故障信息。
选择三个连续WPF分量作为重要WPF分量的原则为：
(1)这三个WPF分量与x(t)的互相关系数之和最小；
(2)这三个WPF分量的模态最高(假设被选的三个WPF分量是WPFs(t)，WPFs+1(t)，WPFs+2(t)，那么s应该尽量小些)。根据这样的原则，WPF1(t)，WPF2(t)和WPF3(t)是三个重要WPF分量。接着，我们分析这三个重要WPF分量，并从中找到转子碰摩故障的特征频率。
分别求出这三个重要WPF分量的包络信号，画出它们的包络幅值谱，如图12所示。从三个重要WPF分量的幅值谱和包络幅值谱中找到碰摩特征频率。由转子局部碰摩故障的特点，高频的碰摩故障信号主要为调幅信号。调幅信号载波频率较高，在400～600hz，而调制频率较低，为转子的转频50hz。从第二个重要的WPF分量的幅值谱中我们看到一个500hz的频率成分，且它的包络幅值谱中存在一个转频50hz的频率成分，因此第二个重要WPF分量包含有转子高频碰摩故障的特征，因此我们可以判断转子发生碰摩故障。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (7)

1.一种用于转子碰摩故障诊断的小波局部均值分解方法，其特征在于，包括以下步骤：
S1、搭建模拟转子碰摩故障实验台，利用位移传感器采集转子碰摩故障的振动信号x(t)；
S2、对采集的振动信号x(t)进行局部均值分解，得到k个PF分量和一个残余分量uk(t)；步骤S2包括以下分步骤：
S21、找到振动信号x(t)的所有局部极值点，求出任意两个相邻的局部极值点的平均值，将所有相邻局部极值点的平均值用直线连接，得到局部均值线段，然后采用滑动平均法进行平滑处理，得到局部均值函数m11(t)；
S22、计算任意两个相邻局部极值点见得包络估计值，将所有相邻包络估计值用直线连接，得到包络估计线段，然后采用滑动平均法进行平滑处理，得到包络估计函数a11(t)；
S23、将局部均值函数m11(t)从振动信号x(t)中分离出来，得到分离信号h11(t)；
S24、用分离信号h11(t)除以包络估计函数a11(t)，对h11(t)进行解调，得到解调信号s11(t)；
S25、判断s11(t)是否为纯调频信号，若是则进入步骤S26，否则用s11(t)代替x(t)，重复步骤S21-S24，迭代n次，直到s1n(t)为纯调频信号为止，进入步骤S26；
S26、把迭代过程中产生的所有包络估计函数相乘，得到包络信号a1(t)；
S27、将包络信号a1(t)和纯调频信号s1n(t)相乘，得到x(t)的第一个PF分量PF1(t)；
S28、将PF1(t)从x(t)中分离出来，得到一个新的信号u1(t)，将u1(t)作为原始数据代替x(t)，重复步骤S21-S27，循环k次，直到uk(t)为一个单调函数为止，将初始振动信号x(t)分解为k个PF分量和一个残余分量uk(t)之和；
S3、采用小波分解对所有的PF分量进行再分解与重构，得到k个WPF分量和一个残余分量vk(t)；步骤S3包括以下分步骤：
S31、对PF1(t)进行若干层小波分解，得到一个小波低频分量和若干个小波高频分量；
S32、剔除PF1(t)的小波低频分量，并把它给PF2(t)，PF1(t)剩余的信号为第一个WPF分量，记作WPF1(t)，PF2(t)获得PF1(t)的小波低频分量之后更新得到PF2 *(t)；
S33、用PF2 *(t)替换PF1(t)，重复步骤S31-S32，直到获得PFk(t)的更新分量PFk *(t)；
S34、对PFk *(t)进行若干层小波分解，得到一个小波低频分量和若干个小波高频分量，剔除PFk *(t)的小波低频分量，并把它给残余分量uk(t)，PFk *(t)的剩余信号为第k个WPF分量，记作WPFk(t)，残余分量uk(t)获得PFk *(t)的小波低频分量，得到残余分量vk(t)；
S4、计算每一个WPF分量与初始振动信号x(t)的互相关系数，从所有的WPF分量中选三个连续WPF分量作为重要WPF分量，求出这三个重要WPF分量的包络信号；分别画出三个重要WPF分量的幅值谱和它们所对应的包络信号的幅值谱，并从中找到转子碰摩故障的特征频率。
2.根据权利要求1所述的小波局部均值分解方法，其特征在于，所述步骤S21中对局部均值线段进行平滑处理时，滑动平均跨度为相邻极值点最大距离的三分之一。
3.根据权利要求1所述的小波局部均值分解方法，其特征在于，所述步骤S22中对包络估计线段进行平滑处理时，滑动平均跨度为相邻极值点最大距离的三分之一。
4.根据权利要求1所述的小波局部均值分解方法，其特征在于，所述步骤S25中判断s1n(t)是否为纯调频信号的具体方法为：
计算s1n(t)的包络估计函数a1(n+1)(t)，设置一个增减量δ>0，判断是否满足条件1-δ≤a1(n+1)(t)≤1+δ，若是则说明s1n(t)是纯调频信号，否则s1n(t)不是纯调频信号。
5.根据权利要求4所述的小波局部均值分解方法，其特征在于，δ的取值在0.001～0.1之间。
6.根据权利要求5所述的小波局部均值分解方法，其特征在于，δ的取值为0.05。
7.根据权利要求1所述的小波局部均值分解方法，其特征在于，所述步骤S4中选择三个连续WPF分量作为重要WPF分量的原则为：
(1)这三个WPF分量与x(t)的互相关系数之和最小；
(2)这三个WPF分量的模态最高。
CN201610257375.5A 2016-04-25 2016-04-25 一种用于转子碰摩故障诊断的小波局部均值分解方法 Expired - Fee Related CN105910805B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610257375.5A CN105910805B (zh) 2016-04-25 2016-04-25 一种用于转子碰摩故障诊断的小波局部均值分解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610257375.5A CN105910805B (zh) 2016-04-25 2016-04-25 一种用于转子碰摩故障诊断的小波局部均值分解方法

Publications (2)

Publication Number Publication Date
CN105910805A CN105910805A (zh) 2016-08-31
CN105910805B true CN105910805B (zh) 2018-06-01

Family

ID=56751717

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610257375.5A Expired - Fee Related CN105910805B (zh) 2016-04-25 2016-04-25 一种用于转子碰摩故障诊断的小波局部均值分解方法

Country Status (1)

Country Link
CN (1) CN105910805B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110531137A (zh) * 2018-05-25 2019-12-03 许继集团有限公司 一种电能质量调节、电压暂变检测方法及动态电压恢复器
CN112183344B (zh) * 2020-09-28 2021-06-01 广东石油化工学院 基于波形和无量纲学习的大机组摩擦故障分析方法及系统
CN113514144B (zh) * 2021-07-28 2022-07-26 郑州轻工业大学 基于电涡流位移传感器的不平衡-碰摩耦合故障检测方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102866027A (zh) * 2012-08-13 2013-01-09 燕山大学 基于lmd和局域时频熵的旋转机械故障特征提取方法
CN103994062A (zh) * 2014-05-13 2014-08-20 山东理工大学 液压泵故障特征信号提取方法
CN104236905A (zh) * 2014-08-26 2014-12-24 中国直升机设计研究所 一种轴承故障诊断方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102866027A (zh) * 2012-08-13 2013-01-09 燕山大学 基于lmd和局域时频熵的旋转机械故障特征提取方法
CN103994062A (zh) * 2014-05-13 2014-08-20 山东理工大学 液压泵故障特征信号提取方法
CN104236905A (zh) * 2014-08-26 2014-12-24 中国直升机设计研究所 一种轴承故障诊断方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
LMD 和改进小波阈值去噪的轴承声发射信号故障诊断;艾延廷等;《科学技术与工程》;20141130;第14卷(第33期);第88-90页 *
应用改进的LMD和小波降噪于滚动轴承故障诊断;刘涛涛等;《噪声与振动控制》;20140430;第34卷(第2期);第152-157页 *

Also Published As

Publication number Publication date
CN105910805A (zh) 2016-08-31

Similar Documents

Publication Publication Date Title
Li et al. An optimized VMD method and its applications in bearing fault diagnosis
Bozchalooi et al. A joint resonance frequency estimation and in-band noise reduction method for enhancing the detectability of bearing fault signals
Hu et al. A new wind turbine fault diagnosis method based on ensemble intrinsic time-scale decomposition and WPT-fractal dimension
Qin et al. Adaptive bistable stochastic resonance and its application in mechanical fault feature extraction
Qin et al. Weak transient fault feature extraction based on an optimized Morlet wavelet and kurtosis
CN105910805B (zh) 一种用于转子碰摩故障诊断的小波局部均值分解方法
CN102866027B (zh) 基于lmd和局域时频熵的旋转机械故障特征提取方法
Gao et al. Spare optimistic based on improved ADMM and the minimum entropy de-convolution for the early weak fault diagnosis of bearings in marine systems
CN110220708A (zh) 一种基于改进小波算法的轴承信号降噪方法
CN108447503B (zh) 基于Hilbert-Huang变换的电机异音检测方法
CN110806315B (zh) 一种基于倒位编辑的齿轮箱复合故障诊断方法
CN105588717A (zh) 一种齿轮箱故障诊断方法
CN106141815A (zh) 一种基于ar模型的高速铣削颤振在线辨识方法
Zeng et al. Normalized complex Teager energy operator demodulation method and its application to fault diagnosis in a rubbing rotor system
CN105046025B (zh) 一种核磁共振多相流测量中各相分离的方法
CN105928701A (zh) 基于相关分析的emd过程中有效imf的判定方法
CN104330258A (zh) 一种基于特征参量的滚动轴承故障灰色关联度辨识方法
Kong et al. Fault feature extraction of planet gear in wind turbine gearbox based on spectral kurtosis and time wavelet energy spectrum
CN111504640B (zh) 一种加权滑动窗二阶同步压缩s变换轴承故障诊断方法
CN112485028B (zh) 振动信号的特征频谱提取方法及机械故障诊断分析方法
Yu et al. Adaptive multiple second-order synchrosqueezing wavelet transform and its application in wind turbine gearbox fault diagnosis
CN109900447B (zh) 一种基于谐波信号分解的循环冲击振动检测方法及系统
Zhang et al. Wind turbine planetary gearbox fault diagnosis via proportion-extracting synchrosqueezing chirplet transform
CN103234750A (zh) 一种基于改进倒频谱法的等高齿锥齿轮故障诊断方法
CN111307426A (zh) 一种基于FrFT-EWT原理的旋转机械故障特征提取方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20180601

Termination date: 20210425

CF01 Termination of patent right due to non-payment of annual fee