本发明属于旋转机械的健康评估领域，尤其是涉及一种基于快速谱峭度分析的风机健康评价方法。
背景技术：
风机是用于输送气体的机械,从能量观点看,它是把原动机的机械能转换为气体能量的一种机械。简单而言,风机是依靠输入的机械能,提高气体压力并排送气体的机械，广泛应用于冶金、石化、电力、城市交通、船舶等领域
风机在工业系统中占有十分重要的作用，以城市交通为例，随着城市化建设的高速发展,地铁成为了缓解城市交通堵塞、建设用地紧张等问题的重要出行方式,是现代化城市的重要交通工具。地铁系统建立在地面以下，地下空间内除活塞风以外的全部通风都完全依靠机械方式进行，风机有效保障了地铁地下内部空间空气与外部空间自然空气的气体交换，尤其当城市地铁系统处于事故工况中时，风机高效地为地下环境排烟、补充新风更成为减少事故损失的重要条件，因此地铁风机在地铁系统中的作用是不可替代的。
风机种类众多，其中包括轴流风机、离心风机、射流风机等，它们是保障工业系统安全、有效运行的关键。为了预防风机发生故障，造成重大的经济损失和安全事故，寻找一种风机的健康评判方法，有效实时地反应风机的健康状态，对预防故障的发生具有重要作用。
所谓的风机健康状况评估，就是运用合理的方法，在风机不停机的情况下，实现智能的风机运行状况估计，合理分析风机运行参数，及时发现风机的异常状态，降低风机发生事故的几率，以减少风机停机时所带来的损失。选择合适的健康评估指标和建立健康评估模型对风机健康评价尤为重要。
一般情况下，风机异常状况可以通过风机的各种运行参数的变化趋势反映出来。一个故障从发生到它造成严重的现场事故，是一个渐变的发展过程，该过程根据一定的趋势轨迹缓慢地变为新的稳态或者逐渐恶化直至风机损毁。传统的健康评估方法是风机各项运行参数出现异常值后才发现，但此时故障已经产生了一定的影响。风机健康状况评估方法可以监测风机故障的早期迹象，使工人有足够的时间进行分析和诊断，从而采取措施防止故障的发生。
目前主流的基于数据驱动的设备健康状况评估方法中，时间序列预测法和神经网络法最为广泛。时间序列预测法因突出时间序列暂不考虑外界因素影响，因而存在着预测误差的缺陷，当遇到外界发生较大变化时，往往会有较大偏差。神经网络对特定问题层数和节点数的未知性，通常导致模型采用了过多的参数，这也就意味着神经网络模型自由度极高，从而使其在训练网络的过程中，会受到噪声影响，导致神经网络模型过拟合，无法准确区分设备的各种运行状态。
技术实现要素：
本发明提供了一种基于快速谱峭度分析的风机健康评价方法，可以从原始振动信号中提取有效的统计学信息，提高对冲击信号的敏感性，进而准确的反应风机的健康程度。
一种基于快速谱峭度分析的风机健康评价方法，其特征在于，包括以下步骤：
(1)采集风机正常工况与不同程度故障工况下的原始振动信号，获得振动数据原始数据库；
(2)对原始数据库中的振动信号x(t)进行短时傅里叶变换，得到时间-频率-振幅的三维图像x(t,f)；
(3)对上述三维图像x(t,f)进行预白化处理，得到预白化之后的时间-频率-振幅三维图像x′(t,f)；
(4)使用预白化后的包络信号作为快速谱峭度的输入，计算最大峭度所在位置及拟合程度最高的快速谱峭度算法分解阶数；
(5)根据最大峭度所在位置及拟合程度最高的快速谱峭度算法分解阶数，确定最优的解调频带宽度和其中心频率，并根据此最优解调频带宽度及其中心频率对白化后的信号进行带通滤波；
(6)对带通滤波后的信号进行傅里叶变换，获得其频谱图，并对所得频谱图进行峰值提取，依次计算所提取的峰值序列中最大值与该序列所有值差值的均值与最大值的比例αi；
(7)对所提取的峰值序列按照幅值降序重新排列获得不同工况峰值提取后的频谱幅值变化趋势数据库；
(8)采集待检测风机振动信号，使用步骤(2)至(7)的方法获得待检测信号所提取的峰值序列中最大值与其余值差值的均值与最大值的比例α0，以及幅值变化趋势序列，使用待检测风机振动信号的幅值变化趋势序列依次与故障工况数据库内每一组幅值变化趋势序列进行相关性分析，获得幅值变化趋势相关系数βi；
(9)构造健康因子γi，
本发明的风机健康评价方法基于快速谱峭度分析，快速谱峭度分析能有效的从含有强噪声信号中发现暂态成分及其在频域中的位置。谱峭度的物理意义就是能够描述信号在频率f处偏离高斯的程度，偏离高斯过程越大，谱峭度的值就越大。瞬态冲击信号占主导的频率段的谱峭度较大，而平稳高斯背景噪声信号占主导频率段的谱峭度很小。因此可以通过计算不同频率位置下信号的峭度值，将冲击成分从原信号中分离出来，从而得到瞬态冲击故障信号，将待测信号与风机振动数据库中不同工况下的数据进行对比分析，得到待测风机的健康因子，反映该风机的健康程度。
步骤(2)中，短时傅里叶变换的输出为：
其中，x(t)为源信号，g(t)为窗函数，g*(t)为g(t)的复共轭函数，t为时间，f为频率，e为自然常数；
为了方便计算机计算，将信号进行离散化处理，具体公式为：
其中，t＞0，f＞0，分别是时间变量和频率变量的采样周期，m、n分别是时间和频率序列。
步骤(3)中，预白化处理在matlab软件中的实现语句为：
defzca_whitening(inputs):
sigma＝np.dot(inputs,inputs.t)/inputs.shape
u,s,v＝np.linalg.svd(sigma)
epsilon＝0.1
zcamatrix＝np.dot(np.dot(u,np.diag(1.0/np.sqrt(np.diag(s) epsilon))),u.t)
returnnp.dot(zcamatrix,inputs)
所述的预白化处理具体步骤为：
(3-1)主分量分析pca预处理：通过协方差矩阵求得原始数据x(t,f)的特征向量u1、u2，然后将每个数据点投影到这两个特征向量上；假设样本集合为{x·j＝[x1j,x1j,…xnj]t|1≤j≤m}，所有样本表示成一个n×m的矩阵，其协方差矩阵为：
其中，m为样本数量；
(3-2)pca白化：对新的每一维坐标做一个标准差归一化处理：
其中x″(t,f)为步骤(3-1)处理的新pca坐标空间，λi为第i维特征对应的特征值，ε是为了避免除数为0的极小值；
(3-3)零相位分量分析zca白化：将步骤(3-2)处理的结果变换到原来的坐标系下。
步骤(4)中，根据试验待处理振动信号的数据量大小，确定拟合程度最高的快速谱峭度算法分解阶数当作信号处理的分解阶数，具体实施方式为：
(4-1)在matlab软件中，根据实际数据量，设置一个初步分解阶数；
(4-2)在该阶数下，使用快速谱峭度算法得到的载波频率及带宽，获得此载波频率及带宽下的信号傅里叶变换后的频谱包络图；
(4-3)根据频谱包络图峰特征，确定分解阶数。
分解阶数确立的原则是根据处理结果的频谱包络图峰的密度来调整的，如果峰密度太低，则降低分解阶数；反之则升高。
步骤(5)中，所述的快速谱峭度算法在matlab软件中为一个以原信号、分解阶数(nlevel)以及采样频率为自变量的函数。
步骤(6)中，傅里叶变换的数学表达式为：
其中，t为时间，f为频率，x(t)为原始信号，x(f)为原始信号的幅值谱。
对傅里叶变换后的频谱图进行峰值提取的限制条件为所提取峰值的值必须大于频谱幅值的均值。
依次计算所提取的峰值序列中最大值与该序列所有值差值的均值与最大值的比例αi，其公式为：
式中，标号i表示数据库中第i组序列，x_max表示所提取的峰值序列中的最大值，n为峰值序列的长度，xj为峰值序列值。
步骤(7)中，对所提取的峰值序列进行降序排列在mantlab中的实现语句为：
x_peaks_new＝sort(x_peaks,'descend')
其中x_peaks_new表示降序后的峰值序列，sort为排序函数，x_peaks表示所提取的峰值序列，'descend'表示降序排列。
步骤(8)中，幅值变化趋势相关系数βi的公式为：
其中，标号i表示第i组序列，βi表示相关系数，x_peaks_new为待检测振动信号峰值的降序序列，x_peaks_new_i为故障数据库中第i个故障信号的峰值降序序列，cov(x_peaks_new,x_peaks_new_i)为二者的协方差，σ(x_peaks_new)，σ(x_peaks_new_i)分别为二者的标准差。
步骤(9)的具体过程为：
按照数据库数据，依次计算健康因子γi，
与现有技术相比，本发明具有以下有益效果：
1.本发明提出的基于谱峭度分析的风机健康评价方法，峭度对冲击信号特别敏感，特别适用于表面损伤类故障，尤其是早期故障的诊断。
2.本发明实现了在实际应用中可以准确、快速地获取风机振动信号的最优解调频带的中心频率与频带宽度进行信号滤波，通过对该段时域信息进行傅里叶变换得到频域信息，将待测风机的振动信息与数据库中不同工况下的振动数据进行对比分析，获得待测风机的健康因子，进而为风机的健康评估提供了一种新的思路与方法。
附图说明
图1为本发明一种基于快速谱峭度分析的风机健康评价方法的流程示意图；
图2为本发明实施例采用快速谱峭度对正常工况下风机分析结果示意图；
图3为本发明实施例采用快速谱峭度对轻微故障工况下风机分析结果示意图；
图4为本发明实施例采用快速谱峭度对中等故障工况下风机分析结果示意图；
图5为本发明实施例采用快速谱峭度对严重故障工况下风机分析结果示意图；
图6为本发明实施例正常工况下基于快速谱峭度滤波处理后的信号的频谱图；
图7为本发明实施例轻微故障工况下基于快速谱峭度滤波处理后的信号频谱图；
图8为本发明实施例中等故障工况下基于快速谱峭度滤波处理后的信号频谱图；
图9为本发明实施例严重故障工况下基于快速谱峭度滤波处理后的信号频谱图；
图10为本发明实施例正常工况下基于快速谱峭度滤波处理后的信号的频谱图峰值提取后的幅值图；
图11为本发明实施例轻微故障工况下基于快速谱峭度滤波处理后的信号的频谱图峰值提取后的幅值图；
图12为本发明实施例中等故障工况下基于快速谱峭度滤波处理后的信号的频谱图峰值提取后的幅值图；
图13为本发明实施例严重故障工况下基于快速谱峭度滤波处理后的信号频谱图峰值提取后的是幅值图；
图14为本发明实施例中4组不同工况组成的频谱幅值变化趋势数据库示意图；
图15为本发明实施例待检测信号的幅值变化趋势图。
具体实施方式
下面结合附图和实施例对本发明做进一步详细描述，需要指出的是，以下所述实施例旨在便于对本发明的理解，而对其不起任何限定作用。
如图1所示，一种基于快速谱峭度分析的风机健康评价方法，包括以下步骤：
s01，采集不同型号风机正常工况、不同故障工况下的原始振动信号，获得振动数据原始数据库。
s02，对s01采集的振动信号x(t)进行短时傅里叶变换，得到时间-频率-振幅的三维图像x(t,f)：
x(t)为源信号，g(t)为窗函数，g*(t)为g(t)的复共轭函数，t为时间，f为频率，e为自然常数。
为了方便计算机计算，通常把信号进行离散化处理，具体公式为：
其中，t＞0，f＞0分别是时间变量和频率变量的采样周期，m、n分别是时间和频率序列。
s03，对s02三维图像进行预白化操作，得到预白化之后的时间-频率-振幅三维图像：
其中，
s04，根据数据量的大小，选取合适的计算分解阶数。在这里，分解阶数确立的原则是根据处理结果的频谱包络图峰的密度来确定的，如果峰密度太低，则降低分解阶数；反之则升高。本次试验数据中，实验采集频率为12800hz，因此采用五阶的分析。
s05，在快速谱峭度的分析结果图中，寻找最大谱峭度值所在位置的中心频率及相应的频率带宽，此数据即为最优解调频带。例如，风机运转频率为30hz，正常工况下的快速快速谱峭度分析结果如图2所示，其最优的解调频带的中心频率为6168.75hz，其频率带宽为262.5hz。轻微故障工况时对应的快速谱峭度分析结果如图3，其最优的解调频带的中心频率为5020.325hz，其频率带宽为196.875hz。中度故障工况时对应的快速谱峭度分析结果如图4，其最优的解调频带的中心频率为4725hz，其频率带宽为3150hz。严重故障时对应的快速谱峭度分析结果如图5，其最优的解调频带的中心频率为1050hz，其频率带宽为2100hz。
s06，依据快速谱峭度方法确定的最优解调频带对原振动数据进行带通滤波，对滤波后的数据进行傅里叶变换得到频率的包络图，例如正常工况下的30hz，正常工况运行的风机对应的滤波后的频谱如图6，轻微故障时对应的滤波后的频谱如图7，中等故障实对应的滤波后的频谱如图8，严重故障时对应的滤波后的频谱如图9。
s07，对频谱幅值进行峰值提取，提取过程的限制条件为峰值的值不得低于频谱幅值的均值，图10即为图6所对应的的频谱峰值提取后的幅值谱，图11即为图7所对应的的频谱峰值提取后的幅值谱，图12即为图8所对应的的频谱峰值提取后的幅值谱，图13即为图9所对应的的频谱峰值提取后的幅值谱。此过程在matlab中的实现语句为：
x_peaks＝findpeaks(x,'minpeakheight',mean(x))；
其中x_peaks为所提取的峰值序列，findpeaks为峰值查找函数，x为原频谱图幅值序列，'minpeakheight'，mean(x)表示峰值高度限制为不小于序列x的均值。
s08，求所提取的峰值序列中，最大值与其余值差值的均值与最大值的比例αi，此过程的计算公式为：
式中，标号i表示数据库中第i组序列，x_max表示所提取的峰值序列中的最大值，n为峰值序列的长度，xj为峰值序列值。
此过程在matlab上的实现语句为：
α＝(sum(x_peaks_max-x_peaks)/(length(x_peaks))*x_peaks_max)
其中，x_peaks_max为幅值序列中最大值，x_peaks为幅值序列，length为数组长度获取函数，sum为求和函数。
s09，对提取到的峰值序列进行降序排列，此过程在mantlab中的实现语句为：
x_peaks_new＝sort(x_peaks,'descend')
其中，x_peaks_new表示降序后的峰值序列，sort为排序函数，x_peaks表示所提取的峰值序列，'descend'表示降序排列。对原始数据库内所有数据均进行以上操作，获得不同工况下滤波后幅值变化趋势数据库，图14即为4组数据组成的幅值变化趋势数据库。
s10，对待检测数据同样进行s02至s08的操作，获得待检测数据所提取的峰值序列中最大值与其余值差值的均值与最大值的比例α0,以及峰值提取后的幅值变化趋势序列。图15即为待检测数据对应的幅值变化趋势图。
s11，使用幅值变化趋势序列与幅值变化趋势数据库内序列依次求相关系数βi，此过程计算公式为：
其中，βi表示相关系数，x_peaks_new为待检测振动信号峰值的降序序列，x_peaks_new_i为故障数据库中第i个故障信号的峰值降序序列，cov(x_peaks_new,x_peaks_new_i)为二者的协方差，σ(x_peaks_new)，σ(x_peaks_new_i)分别为二者的标准差。
s12，计算健康因子γi，
s13，寻找最大γi值所对应的数据库中数据所处工况，即为待检测数据的健康评估。结果如下表1所示，待检测数据的最大γi值为10.3129，对应的数据为数据库中第4组数据，此时显示待检测风机处于严重故障状态。
表1
以上所述的实施例对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的具体实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换，均应包含在本发明的保护范围之内。
1.一种基于快速谱峭度分析的风机健康评价方法，其特征在于，包括以下步骤：
(1)采集风机正常工况与不同程度故障工况下的原始振动信号，获得振动数据原始数据库；
(2)对原始数据库中的振动信号x(t)进行短时傅里叶变换，得到时间-频率-振幅的三维图像x(t,f)；
(3)对上述三维图像x(t,f)进行预白化处理，得到预白化之后的时间-频率-振幅三维图像x′(t,f)；
(4)使用预白化后的包络信号作为快速谱峭度的输入，计算最大峭度所在位置及拟合程度最高的快速谱峭度算法分解阶数；
(5)根据最大峭度所在位置及拟合程度最高的快速谱峭度算法分解阶数，确定最优的解调频带宽度和其中心频率，并根据此最优解调频带宽度及其中心频率对白化后的信号进行带通滤波；
(6)对带通滤波后的信号进行傅里叶变换，获得其频谱图，并对所得频谱图进行峰值提取，依次计算所提取的峰值序列中最大值与该序列所有值差值的均值与最大值的比例αi；
(7)对所提取的峰值序列按照幅值降序重新排列获得不同工况峰值提取后的频谱幅值变化趋势数据库；
(8)采集待检测风机振动信号，使用步骤(2)至(7)的方法获得待检测信号所提取的峰值序列中最大值与其余值差值的均值与最大值的比例α0，以及幅值变化趋势序列，使用待检测风机振动信号的幅值变化趋势序列依次与故障工况数据库内每一组幅值变化趋势序列进行相关性分析，获得幅值变化趋势相关系数βi；
(9)构造健康因子γi，
2.根据权利要求1所述的基于快速谱峭度分析的风机健康评价方法，其特征在于，步骤(2)中，为了方便计算机计算，将短时傅里叶变换过程进行信号的离散化处理，具体公式为：
其中，t＞0，f＞0，分别是时间变量和频率变量的采样周期，m、n分别是时间和频率序列。
3.根据权利要求1所述的基于快速谱峭度分析的风机健康评价方法，其特征在于，步骤(3)中，所述的预白化处理具体步骤为：
(3-1)主分量分析pca预处理：通过协方差矩阵求得原始数据x(t,f)的特征向量u1、u2，然后将每个数据点投影到这两个特征向量上；假设样本集合为{x·j＝[x1j,x1j,…xnj]t|1≤j≤m}，所有样本表示成一个n×m的矩阵，其协方差矩阵为：
其中，m为样本数量；
(3-2)pca白化：对新的每一维坐标做一个标准差归一化处理：
其中x″(t,f)为步骤(3-1)处理的新pca坐标空间，λi为第i维特征对应的特征值，ε是为了避免除数为0的极小值；
(3-3)零相位分量分析zca白化：将步骤(3-2)处理的结果变换到原来的坐标系下。
4.根据权利要求1所述的基于快速谱峭度分析的风机健康评价方法，其特征在于，步骤(4)中，根据试验待处理振动信号的数据量大小，确定拟合程度最高的快速谱峭度算法分解阶数当作信号处理的分解阶数，具体实施方式为：
(4-1)在matlab软件中，根据实际数据量，设置一个初步分解阶数；
(4-2)在该阶数下，使用快速谱峭度算法得到的载波频率及带宽，获得此载波频率及带宽下的信号傅里叶变换后的频谱包络图；
(4-3)根据频谱包络图峰特征，确定分解阶数。
5.根据权利要求1所述的基于快速谱峭度分析的风机健康评价方法，其特征在于，步骤(6)中，对傅里叶变换后的频谱图进行峰值提取的限制条件为所提取峰值的值必须大于频谱幅值的均值。
6.根据权利要求1所述的基于快速谱峭度分析的风机健康评价方法，其特征在于，步骤(6)中，依次计算所提取的峰值序列中最大值与该序列所有值差值的均值与最大值的比例αi，其公式为：
式中，标号i表示数据库中第i组序列，x_max表示所提取的峰值序列中的最大值，n为峰值序列的长度，xj为峰值序列值。
7.根据权利要求1所述的基于快速谱峭度分析的风机健康评价方法，其特征在于，步骤(8)中，幅值变化趋势相关系数βi的公式为：
其中，标号i表示第i组序列，βi表示相关系数，x_peaks_new为待检测振动信号峰值的降序序列，x_peaks_new_i为故障数据库中第i个故障信号的峰值降序序列，cov(x_peaks_new,x_peaks_new_i)为二者的协方差，σ(x_peaks_new)，σ(x_peaks_new_i)分别为二者的标准差。
8.根据权利要求1所述的基于快速谱峭度分析的风机健康评价方法，其特征在于，步骤(9)的具体过程为：
按照数据库数据，依次计算健康因子γi，
