本发明属于信号处理技术领域，具体涉及一种适用于提取轴承故障特征的自适应滤波方法。
背景技术：
滚动轴承是旋转机械的传动关键部件，也是极易发生故障的部件，其常见的故障类型包括：内环故障，外环故障，滚动体故障和保持架故障等，每种故障类型都对应一种理论特征频率。振动分析法是轴承故障诊断的常用方法，其诊断原理为：首先采集轴承工作时的振动信号，再利用信号处理方法提取振动信号中反映故障特征的成分，最后对提取的故障特征成分进行频谱分析，通过频谱峰值处频率与理论故障特征频率对比，最终诊断出轴承的故障类型。
如何提取故障特征是轴承故障诊断的关键。现有方法中，时域分析利用振动信号的时域参数(如峰值，均值，方差，峭度等)组成特征向量；频域分析利用傅里叶变换进行频谱分析以观察频域特征；小波分析利用多尺度小波基对振动信号进行小波分解提取特征成分；稀疏表示方法利用冗余的原子字典对振动信号进行稀疏分解提取特征成分等等。这些方法在轴承故障特征提取具体应用中取得良好效果，但受工况复杂性影响，这些方法也存在一些缺陷，如：轴承早期故障阶段，振动信号中反映故障特征的成分很微弱，往往被转频和其它谐波以及噪声所覆盖，难以从时域或频域提取故障特征，小波分析的恒品质因子性限制了提取故障特征的灵活性和自适应性，稀疏分解方法的选择对提取故障特征的准确性有一定的影响。
技术实现要素：
针对上述现有技术中描述的不足，本发明提供一种适用于提取轴承故障特征的自适应滤波方法，本发明具有自动确定滤波分解总次数和自适应筛选滤波频带分量的特点，在提取轴承早期故障特征的自适应性方面和准确性方面具有一定的优势。
为解决上述技术问题，本发明所采用的技术方案如下：一种适用于提取轴承故障特征的自适应滤波方法，步骤如下：
s1，采集时域轴承振动信号x(n)；
s2，对轴承振动信号x(n)进行离散傅里叶变换，得到频域轴承振动信号x(k)；
式中：n为轴承振动信号x(n)的长度；i为虚单位，满足i2＝-1。
s3，计算滤波分解的总次数j；
式中：round(·)表示取不大于●的正整数。
s4，进行j次滤波分解，并得到时域滤波频带分量集v(n)。
具体步骤为：s4.1，给循环次数j赋初值为1，并令初始输入信号v0,0(k)＝x(k)。
s4.2，计算滤波分解时的低频分量长度n0,j和高频分量长度n1,j：
s4.3，计算滤波分解时的低频分量v0,j(k)和高频分量v1,j(k)：
低频分量v0,j(k)的计算公式为：
高频分量v1,j(k)的计算公式为：
其中，θ(·)函数为：
v0,0(k)＝x(k)(7)。
s4.4，输出高频分量v1,j(k)，并将低频分量v0,j(k)作为下一次滤波分解的输入信号。
s4.5，循环步骤s4.2-s4.4，直至滤波分解次数j等于总次数j时结束，且使输出的各高频分量和最后一次滤波分解得到的低频分量构成频域频带分量集v(k)，
v(k)＝(v1,1(k),v1,2(k),…,v1,j(k),…,v1,j(k),v0,j(k))。
s4.6，对步骤s4.5中得到的频域频带分量集v(k)的各频带分量元素分别进行离散傅里叶逆变换，得到时域滤波频带分量集v(n)，且
v(n)＝(v1,1(n),v1,2(n),…,v1,j(n),…v1,j(n),v0,j(n))；
离散傅里叶逆变换的通式为：
式中：vq(k)为频域频带分量集v(k)中第q个频带分量，nvq为频带分量vq(k)的长度，vq(n)为时域滤波频带分量集v(n)中第q个滤波频带分量；且q＝1,2,...,j,j+1。
s5，对步骤s4得到的时域滤波频带分量集v(n)中的各时域滤波频带分量vq(n)，分别求取峭度值kq，通式为：
式中：vq为待求峭度值的滤波频带分量，μ为滤波频带分量vq的平均值，σ为滤波频带分量vq的标准差。
s6，从时域滤波频带分量集v(n)中利用峭度值自适应筛选出m个反映轴承故障特征的滤波频带分量vq。
具体步骤为：判断每个峭度值kq是否显著大于3，若成立，则筛选出峭度值kq显著大于3所对应的滤波频带分量vq；若不成立，则舍弃；
而显著性由莱以特准则判别，判别式为：
|kq-μk|＞3σk(10)；
kq为第q个峭度值，μk为j+1个峭度值的均值；σk为j+1峭度值的标准差。
s7，对步骤s6筛选出的各滤波频带分量分别进行阈值处理。
阈值处理的步骤为：
s7.1，计算
s7.2，判断
s7.3，重复步骤s7.1-s7.2，直至滤波频带分量vq中的每个分量值都处理一次。
本发明首先利用离散傅里叶变换将轴承振动信号转化至频域，再在频域进行多次滤波分解，利用振动信号本身长度自动确定滤波分解的总次数，并应用莱以特显著性判别准则筛选出峭度值显著大于3的滤波频带分量，作为提取的轴承故障特征成分，最后通过阈值处理对筛选出的故障特征成分进行降噪。本发明避免了现有方法中人为选择参数的影响，根据操作步骤对时域轴承振动信号进行自适应滤波处理，提取的轴承故障特征成分时域波形的周期性明显，包络谱图上的峰值特征频率与理论特征频率匹配准确，凭借此发明可以准确诊断出滚动轴承的常见故障类型。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一部分实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明的系统流程图。
图2是本发明采集的故障轴承时域振动信号波形图。
图3是本发明各滤波频带分量的峭度值柱形图。
图4是本发明筛选出的峭度值显著大于3的滤波频带分量波形图。
图5是本发明阈值处理后的滤波频带分量波形图。
图6是图5中滤波频带分量波形的包络谱分布图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
一种适用于提取轴承故障特征的自适应滤波方法，如图1所示，步骤如下：
s1，采集时域轴承振动信号x(n)。
s2，对轴承振动信号x(n)进行离散傅里叶变换，得到频域轴承振动信号x(k)；
式中：n为轴承振动信号x(n)的长度；i为虚单位，满足i2＝-1。
s3，计算滤波分解的总次数j；
式中：round(·)表示取不大于·的正整数。
s4，进行j次滤波分解，并得到时域滤波频带分量集v(n)。
具体步骤为：s4.1，给循环次数j赋初值并令初始输入信号v0,0(k)＝x(k)。s4.2，计算滤波分解时的低频分量长度n0,j和高频分量长度n1,j：
s4.3，计算滤波分解时的低频分量v0,j(k)和高频分量v1,j(k)：
低频分量v0,j(k)的计算公式为：
高频分量v1,j(k)的计算公式为：
其中，θ(·)函数为：
v0,0(k)＝x(k)(7)。
s4.4，输出高频分量v1,j(k)，并将低频分量v0,j(k)作为下一次滤波分解的输入信号。
s4.5，循环步骤s4.2-s4.4，直至滤波分解次数j等于总次数j时结束，且使输出的各高频分量和最后一次滤波分解得到的低频分量构成频域频带分量集v(k)，且v(k)＝(v1,1(k)，v1,2(k),…,v1,j(k),…,v1,j(k),v0,j(k))。
s4.6，对步骤s4.5中得到的频域频带分量集v(k)的各频带分量元素分别进行离散傅里叶逆变换，得到时域滤波频带分量集v(n)，且
v(n)＝(v1,1(n)，v1,2(n),…,v1,j(n),…v1,j(n),v0,j(n))；
离散傅里叶逆变换的通式为：
式中：vq(k)为频域频带分量集v(k)中第q个频带分量，nvq为频带分量vq(k)的长度，vq(n)为时域滤波频带分量集v(n)中第q个滤波频带分量；且q＝1,2,...,j,j+1。
s5，对步骤s4得到的时域滤波频带分量集v(n)中的各时域滤波频带分量vq(n)，分别求取峭度值kq，通式为：
式中：vq为待求峭度值的滤波频带分量，μ为滤波频带分量vq的平均值，σ为滤波频带分量vq的标准差。
s6，从时域滤波频带分量集v(n)中利用峭度值自适应筛选出m个反映轴承故障特征的滤波频带分量vq。
具体步骤为：判断每个峭度值kq是否显著大于3，若成立，则筛选出峭度值kq显著大于3所对应的滤波频带分量vq；若不成立，则舍弃；
而显著性由莱以特准则判别，判别式为：
|kq-μk|＞3σk(10)；
kq为第q个峭度值，μk为j+1个峭度值的均值；σk为j+1峭度值的标准差。
s7，对步骤s6筛选出的各滤波频带分量分别进行阈值处理。
阈值处理的步骤为：
s7.1，计算
s7.2，判断
s7.3，重复步骤s7.1-s7.2，直至滤波频带分量vq中的每个分量值都处理一次。
下面以具体实例对本发明进一步说明。
采集型号为skf6205-2rs轴承发生早期外环故障时(故障点为直径为0.007英尺，深度为0.011英尺的小凹点)的振动数据，采样频率为12khz，图2为采样时间为1s的振动数据波形。
步骤1：对图2所示的振动数据按照公式1进行离散傅里叶变换，其中，振动数据的总点数为n＝12000，得到x(k)。
步骤2：对x(k)进行滤波分解，其中滤波总次数j由公式2自适应确定，本实例中得到滤波总次数j＝18。
具体滤波分解过程为：在第1次滤波分解时，即，j＝1，输入x(k)，先分别确定低频分量v0,1(k)的长度和高频分量v1,1(k)的长度，再分别计算出低频分量v0,1(k)和高频分量v1,1(k)，计算完毕后，输出高频分量v1,1(k)，低频分量v0,1(k)作为下一次滤波分解的输入；
在第2次滤波分解时，即，j＝2，输入v0,1(k)，先分别确定低频分量v0,2(k)和高频分量v1,2(k)的长度，再分别计算出低频分量v0,2(k)和高频分量v1,2(k)，计算完毕后，输出高频分量v1,1(k)，低频分量v0,2(k)作为下一次滤波分解的输入。
以此类推，直到分解到j次，即，18次结束。
分解结束，依次得到高频分量v1,1(k),v1,2(k),…,v1,18(k)和最后一次分解得到的低频分量v0,18(k)。
并对高频分量v1,1(k),v1,2(k),…,v1,18(k)和最后一次分解得到的低频分量v0,18(k)进行离散傅里叶逆变换，得到v1,1(n)，v1,2(n),…,v1,18(n)和v0,18(n)。
步骤3：求v1,1(n)，v1,2(n),…,v1,18(n)和v0,18(n)各自的峭度值，得到k1，k1,…,k18和k19，其中：k19为v0,18(n)的峭度值，各峭度值柱形图如图3所示。
步骤4：利用峭度值自适应筛选出反映轴承故障特征的滤波频带分量。
具体筛选过程为：首先计算k1，k1,…,k18和k19的平均值μk和标准差σk；再逐一计算|k1-μk|,|k2-μk|,…,|k18-μk|和|k19-μk|；最后利用莱以特显著性判别准则，筛选出峭度值显著大于3所对应的滤波频带分量v10和v11。
v10和v11的波形如图4所示。
步骤5：对筛选出的峭度值显著大于3所对应的滤波频带分量v10和v11进行阈值处理，具体操作步骤为:
步骤5.1：计算
特殊情况：当m＝1时(即起点情况)，令dq,m-1＝dq,m；当m＝nq时，即末端点情况，令dq,m+1＝dq,m；且nq为滤筛选出的滤波频带分量vq的长度。
步骤5.2：判断
步骤5.3：重复步骤s:5.1-s5.2，直至滤波频带分量vq中的每个分量值都处理一次。
阈值处理后的滤波频带分量v10和v11，就是所提取的轴承故障特征，波形如图5所示。从图5可以看出，经过本发明提取的两个滤波频带分量所含冲击成分的周期性明显体现了出来。进一步分析所提取的这两个滤波频带分量的包络谱，如图6所示。从图6中可看出，两者谱图中最大峰值处对应的频率均为105.6hz，此值与轴承理论外环故障特征频率105.3hz非常接近，由此可诊断该轴承存在外环故障点，这一诊断结论与事实相符。
上面所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。