本发明涉及自动化检测技术，特别设计一种振动信号的特征频谱提取方法及机械故障诊断分析方法。
背景技术：
工业生产和日常生活领域中，旋转机械的应用极为广泛，因此旋转机械的故障是关系到生命和生产安全的重要问题，对旋转机械的故障诊断具有突出的实际价值。利用振动信号可以对不同的机械故障进行分析，根据特征频谱的不同判断机械故障类型。
但是采集的振动信号中不可避免地存在噪声，而噪声的频谱存在于整个频域范围内，在故障分析中，旋转类机械结构的故障信号能量往往较低，噪声信号的存在降低了故障信号的信噪比，增加了故障诊断的难度。
现有滤波方法中，滤波的中心频率一般是先验性的或通过搜索方式确定，前者主观性影响较为严重，而搜索方式的计算量相当大。
技术实现要素：
本发明要解决的技术问题是提供一种振动信号的特征频谱提取方法，可以提高振动信号故障诊断的准确性，有助于旋转机械故障诊断分析。
为解决上述技术问题，本发明提供的振动信号的特征频谱提取方法，其包括以下步骤：
一.计算振动信号x(t)的时频分布v(f,t)；
二.在时频域中，根据交错定理获得切比雪夫范数意义下的最佳逼近函数m(t)；
三.利用最佳逼近函数，确定分数阶傅里叶变换的最佳阶次p；
四.计算最佳阶次下的振动信号x(t)的分数阶傅里叶变换结果xp(u)，其中u为分数阶傅里叶域的横轴坐标，确定振动信号在分数阶域内的聚集中心位置，以聚集中心位置为中心在分数阶域内对振动信号分数阶傅里叶变换结果xp(u)进行滤波；
五.对滤波后的振动信号分数阶傅里叶变换结果xp(u)进行分数阶傅里叶逆变换，获得时域下的振动信号滤波结果，对其进行包络解调，利用离散傅氏变换的快速算法变换得到振动信号特征频谱。
较佳的，步骤一中采用一维小波变换计算振动信号的时频分布，以morlet小波作为母小波，进行计算得到振动信号x(t)的时频分布结果。
较佳的，步骤二中，所述切比雪夫范数意义下的最佳逼近函数m(t)的具体计算过程为：
1)选取振动信号x(t)的时频分布v(f,t)幅值最高的三个点作为交错点，分别为v1(f1,t1)、v2(f2,t2)、v3(f3,t3)，其中t为时间变量，ti表示第i个交错点横坐标位置，f为频率变量，fi表示第i个交错点纵坐标位置，i＝1,2,3；
2)最佳逼近函数m(t)取一次多项式，m(t)＝c0 c1*t，t为时间变量，c0为常数项系数，c1为一次项系数；
3)设切比雪夫范数‖r(t)‖c＝α，以f1-m(t1)＝-α，f2-m(t2)＝α，f3-m(t3)＝-α三个方程组成方程组，计算得到最佳逼近函数m(t)的常数项系数c0及一次项系数c1和α的值，即获得切比雪夫范数意义下的最佳逼近函数m(t)。
较佳的，步骤三中，利用最佳逼近函数m(t)的一次项系数c1，确定分数阶傅里叶变换的最佳阶次p，
较佳的，步骤四中，将计算得到的最佳阶次p的数值增加2后除以4得到余数re，根据余数数值确定对应的分数阶傅里叶变换结果。
较佳的，当re＝2时，分数阶傅里叶变换结果为xp(u)＝x(t)；
当re＝0时，分数阶傅里叶变换结果为xp(u)＝x(-t)；
当re＝3时，分数阶傅里叶变换为普通傅里叶变换；
re＝1时，分数阶傅里叶变换为普通傅里叶逆变换；
re为其余值时，分数阶傅里叶变换结果xp(u)为
其中，β为旋转角度，j为虚数单位。
较佳的，在对振动信号进行分数阶傅里叶变换后，取振动信号分数阶傅里叶变换结果xp(u)幅值的绝对值，在分数阶域内确定幅值最高点的位置，以该位置为振动信号在分数阶域内的聚集中心位置。
较佳的，步骤五中，对滤波后的振动信号分数阶傅里叶变换结果xp(u)进行分数阶傅里叶逆变换，获得时域下的振动信号滤波结果，随后对时域下的信号滤波结果进行经验模态分解，对分解得到的固有模态分量进行包络解调，然后利用离散傅氏变换的快速算法变换得到振动信号频谱。
为解决上述技术问题，本发明提供的一种根据上述振动信号的特征频谱提取方法得到的振动信号频谱进行机械故障诊断分析方法，所述机械为旋转机械。
本发明的振动信号的特征频谱提取方法，基于分数阶傅里叶变换自适应地对振动信号进行滤波，然后通过包络分析得到振动信号特征频谱，可以提高振动信号故障诊断的准确性，改善包络分析方法在旋转结构故障诊断的应用效果，有助于旋转机械故障诊断分析。
附图说明
为了更清楚地说明本发明的技术方案，下面对本发明所需要使用的附图作简单的介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明的时域振动信号的特征频谱提取方法一实施例的流程图；
图2是仿真信号x(t)的时域波形图；
图3是对仿真信号未滤波直接进行经验模态分解的固有模态分量频谱图；
图4是图3的固有模态分量包络解调后的频谱图；
图5是对仿真信号经过滤波后处理后进行经验模态分解的固有模态分量频谱图；
图6是图5的固有模态分量包络解调后的频谱图。
具体实施方式
下面将结合附图，对本发明中的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其它实施例，都属于本发明保护的范围。
实施例一
如图1所示，振动信号的特征频谱提取方法包括以下步骤：
一.计算振动信号x(t)的时频分布v(f,t)；
二.在时频域中，根据交错定理获得切比雪夫范数意义下的最佳逼近函数m(t)；
三.利用最佳逼近函数，确定分数阶傅里叶变换(frft)的最佳阶次p；
四.计算最佳阶次下的振动信号x(t)的分数阶傅里叶变换结果xp(u)，其中u为分数阶傅里叶域的横轴坐标，确定振动信号在分数阶域内的聚集中心位置，以聚集中心位置为中心在分数阶域内对振动信号分数阶傅里叶变换结果xp(u)进行滤波；
五.对滤波后的振动信号分数阶傅里叶变换结果xp(u)进行分数阶傅里叶逆变换，获得时域下的振动信号滤波结果，对其进行包络解调，利用fft(离散傅氏变换的快速算法)变换得到振动信号特征频谱。
实施例一的振动信号的特征频谱提取方法，基于分数阶傅里叶变换自适应地对振动信号进行滤波，然后通过包络分析得到振动信号特征频谱，可以提高振动信号故障诊断的准确性，改善包络分析方法在旋转结构故障诊断的应用效果，有助于旋转机械故障诊断分析。
实施例二
基于实施例一的振动信号的特征频谱提取方法，步骤一中采用一维小波变换计算振动信号的时频分布，以morlet小波作为母小波，进行计算得到振动信号x(t)的时频分布结果。
实施例三
基于实施例一的振动信号的特征频谱提取方法，步骤二中，所述切比雪夫范数意义下的最佳逼近函数的具体计算过程为：
1)选取振动信号x(t)的时频分布v(f,t)幅值最高的三个点作为交错点，分别为v1(f1,t1)、v2(f2,t2)、v3(f3,t3)，其中t为时间变量，ti表示第i个交错点横坐标位置，f为频率变量，fi表示第i个交错点纵坐标位置，i＝1,2,3；
2)最佳逼近函数m(t)取一次多项式，m(t)＝c0 c1*t，t为时间变量，c0为常数项系数，c1为一次项系数；
3)在交错点位置，最佳逼近函数m(ti)与振动信号频率值fi间的偏差值r(ti)可以表示为:
r(ti)＝fi-m(ti)(1)；
根据交错定理有：r(ti)＝-r(ti-1)＝±‖r(t)‖c(2)
其中‖r(t)‖c为切比雪夫范数，表达式为：‖r(t)‖c＝max|fi-m(ti)|，tmin≤t≤tmax(3)；
设α＝max|fi-m(ti)|，根据三个交错点可以获得三个方程：
f1-m(t1)＝-α，f2-m(t2)＝α，f3-m(t3)＝-α(4)
以(4)中的三个方程组成方程组，以此计算得到最佳逼近函数m(t)的常数项系数c0及一次项系数c1和α的值，即获得切比雪夫范数意义下的最佳逼近函数m(t)。
例如，采用仿真信号x(t)＝x1(t) x2(t) x3(t)，其中x1(t)为谐波信号，具体表达式为cos(2*π*20*t) 2*cos(2*π*55*t)，对应的频率值为20hz和55hz，x2(t)为衰减冲击信号，以此模拟旋转机械故障信号，具体表达式为e(-2*t)*sin(2*π*12*t)，对应的频率值为12hz，x3(t)为随机白噪声，仿真信号的信噪比为-9，信号采样频率为1024hz，数据长度为1024，可以计算得到c1＝4.71。图2为仿真信号x(t)的时域波形图。
实施例三的振动信号的特征频谱提取方法，选取时频图中幅值最高的三个点作为交错点，根据交错定理确定切比雪夫范数意义下的最佳逼近函数。
实施例四
基于实施例一的振动信号的特征频谱提取方法，步骤三中，利用最佳逼近函数m(t)的一次项系数c1，确定分数阶傅里叶变换(frft)的最佳阶次p，
实施例五
基于实施例一的振动信号的特征频谱提取方法，步骤四中，将计算得到的最佳阶次p的数值增加2后除以4得到余数re，根据余数数值确定对应的分数阶傅里叶变换结果。
较佳的，设振动信号原时域上的信号为x(t)，分数阶傅里叶变换后的结果为xp(u)，其中u为分数阶傅里叶域的横轴坐标；
当re＝2时，分数阶傅里叶变换结果为xp(u)＝x(t)；
当re＝0时，分数阶傅里叶变换结果为xp(u)＝x(-t)；
当re＝3时，分数阶傅里叶变换为普通傅里叶变换；
re＝1时，分数阶傅里叶变换为普通傅里叶逆变换；
re为其余值时，分数阶傅里叶变换结果xp(u)为
其中，β为旋转角度，j为虚数单位。
例如，最佳逼近函数m(t)＝c0 c1*t，当c1＝4.71，
较佳的，在对振动信号进行分数阶傅里叶变换后，取振动信号分数阶傅里叶变换结果xp(u)幅值的绝对值，在分数阶域内确定幅值最高点的位置，以该位置为振动信号在分数阶域内的聚集中心位置。
实施例六
基于实施例一的振动信号的特征频谱提取方法，步骤五中，对滤波后的振动信号分数阶傅里叶变换结果xp(u)进行分数阶傅里叶逆变换，获得时域下的振动信号滤波结果，随后对时域下的信号滤波结果进行经验模态分解，对分解得到的固有模态分量进行包络解调，然后利用fft(离散傅氏变换的快速算法)变换得到振动信号频谱。
如果对仿真振动信号未滤波直接进行经验模态分解，取前三个固有模态分量进行fft，得到的频谱结果如图3所示，从图中可以看出，固有模态分量1(imf1)中噪声干扰严重，无法获得有效的信号频率，固有模态分量2(imf2)和固有模态分量3(imf3)能够获得谐波信号的频率——55hz和20hz，但是未能得到关心的故障频率。对上述得到的固有模态分量采用包络的方法计算各固有模态分量的频谱，结果如图4所示，其中固有模态分量1(imf1)的噪声干扰严重，无法分辨出明显的特征频谱，固有模态分量2(imf2)在故障信号频率附近存在峰值，但是噪声干扰仍明显存在，固有模态分量3(imf3)在故障信号频率的1/2倍频附近存在峰值。可见由于噪声干扰严重，采用包络得到的结果仍难以有效分析故障特征频率。
如果对仿真振动信号经过滤波后处理后再进行经验模态分解，取前三个固有模态分量进行fft，得到的频谱结果如图5所示，从图中可以看出，噪声干扰已经基本被滤除，固有模态分量1(imf1)和固有模态分量2(imf2)清晰地获得了谐波信号的频率值，同时固有模态分量3(imf3)清晰地获得了模拟故障信号的频率——12hz，可见基于分数阶傅里叶变换的自适应滤波的良好效果。对上述得到的固有模态分量采用包络的方法计算各固有模态分量的频谱，结果如图6所示，固有模态分量1(imf1)在故障信号频率处存在峰值，通过滤波后，噪声干扰被较好地滤除掉，相较于未滤波的包络解调结果，滤波后的结果特征频率清晰可见，获得了较好的分析结果。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明保护的范围之内。
1.一种振动信号的特征频谱提取方法，其特征在于，包括以下步骤：
一.计算振动信号x(t)的时频分布v(f，t)；
二.在时频域中，根据交错定理获得切比雪夫范数意义下的最佳逼近函数m(t)；
三.利用最佳逼近函数，确定分数阶傅里叶变换的最佳阶次p；
四.计算最佳阶次下的振动信号x(t)的分数阶傅里叶变换结果xp(u)，其中u为分数阶傅里叶域的横轴坐标，确定振动信号在分数阶域内的聚集中心位置，以聚集中心位置为中心在分数阶域内对振动信号分数阶傅里叶变换结果xp(u)进行滤波；
五.对滤波后的振动信号分数阶傅里叶变换结果xp(u)进行分数阶傅里叶逆变换，获得时域下的振动信号滤波结果，对其进行包络解调，利用离散傅氏变换的快速算法变换得到振动信号特征频谱。
2.根据权利要求1所述的振动信号的特征频谱提取方法，其特征在于，
步骤一中采用一维小波变换计算振动信号的时频分布，以morlet小波作为母小波，进行计算得到振动信号x(t)的时频分布结果。
3.根据权利要求1所述的振动信号的特征频谱提取方法，其特征在于，
步骤二中，所述切比雪夫范数意义下的最佳逼近函数m(t)的具体计算过程为：
1)选取振动信号x(t)的时频分布v(f，t)幅值最高的三个点作为交错点，分别为v1(f1，t1)、v2(f2，t2)、v3(f3，t3)，其中t为时间变量，ti表示第i个交错点横坐标位置，f为频率变量，fi表示第i个交错点纵坐标位置，i＝1，2，3；
2)最佳逼近函数m(t)取一次多项式，m(t)＝c0 c1*t，t为时间变量，c0为常数项系数，c1为一次项系数；
3)设切比雪夫范数||r(t)||c＝α，以f1-m(t1)＝-α，f2-m(t2)＝α，f3-m(t3)＝-α三个方程组成方程组，计算得到最佳逼近函数m(t)的常数项系数c0及一次项系数c1和α的值，即获得切比雪夫范数意义下的最佳逼近函数m(t)。
4.根据权利要求1所述的振动信号的特征频谱提取方法，其特征在于，
步骤三中，利用最佳逼近函数m(t)的一次项系数c1，确定分数阶傅里叶变换的最佳阶次p，
5.根据权利要求1所述的振动信号的特征频谱提取方法，其特征在于，
步骤四中，将计算得到的最佳阶次p的数值增加2后除以4得到余数re，根据余数数值确定对应的分数阶傅里叶变换结果。
6.根据权利要求5所述的振动信号的特征频谱提取方法，其特征在于，
当re＝2时，分数阶傅里叶变换结果为xp(u)＝x(t)；
当re＝0时，分数阶傅里叶变换结果为xp(u)＝x(-t)；
当re＝3时，分数阶傅里叶变换为普通傅里叶变换；
re＝1时，分数阶傅里叶变换为普通傅里叶逆变换；
re为其余值时，分数阶傅里叶变换结果xp(u)为
其中，β为旋转角度，j为虚数单位。
7.根据权利要求1所述的振动信号的特征频谱提取方法，其特征在于，
在对振动信号进行分数阶傅里叶变换后，取振动信号分数阶傅里叶变换结果xp(u)幅值的绝对值，在分数阶域内确定幅值最高点的位置，以该位置为振动信号在分数阶域内的聚集中心位置。
8.根据权利要求1所述的振动信号的特征频谱提取方法，其特征在于，
步骤五中，对滤波后的振动信号分数阶傅里叶变换结果xp(u)进行分数阶傅里叶逆变换，获得时域下的振动信号滤波结果，随后对时域下的信号滤波结果进行经验模态分解，对分解得到的固有模态分量进行包络解调，然后利用离散傅氏变换的快速算法变换得到振动信号频谱。
9.一种根据权利要求1到8所述的振动信号的特征频谱提取方法得到的振动信号频谱进行机械故障诊断分析方法，其特征在于，所述机械为旋转机械。
技术总结