本发明涉及旋转机械监测领域，具体设计一种基于振动信号的旋转机械预警方法。
背景技术：
国内炼油与化工生产装置规模大型化发展趋势明显，与其配套的旋转机械设备也向大型化、高速化、自动化和智能化方向发展，设备故障发生导致的非计划停机不仅会造成巨大的经济损失，而且可能会带来灾难性的火灾、爆炸等安全事故，实现预测性维修对于确保设备运行安全、可靠具有重要作用。故障类型按故障发生、发展的过程分为突发性故障和渐变性故障，一般渐变性故障具有可检测性。研究设备早期故障检测预警技术，提前检测、告警设备即将发生的轻微或不正常故障征兆，使运行维护人员来预防故障或为故障的发生做好充足准备，并最大限度地减少计划外维修带来的损失具有重要的工程应用价值和实践意义。
目前，我国工业企业在役在线监测故障诊断系统设备故障告警，采用当振动达到某一规定的振动幅值或振动发生显著变化时进行报警的方法，一般无法提前发现早期故障征兆，难以及时判断设备早期故障，存在较多的错误报警、漏报警，给设备操作维护人员造成了“报警疲劳”；在固定阈值报警线以下运行的设备，往往缺乏有效的状态劣化趋势告警，从设备报警到联锁停机，有时p-f间隔期很短，往往来不及采取预防性维修措施，非计划停机屡次发生造成巨大经济财产和安全损失。在工业企业，实现设备预测性维修还存在一定的技术挑战。
技术实现要素：
针对现有技术中的不足，本发明提供了一种基于振动信号的旋转机械设备预警方法，可以准确可靠的探测设备的早期故障并报警。结合上述背景技术的不足，本发明致力于实现的技术目的，主要采用了哪些技术点，可以归纳整理下，以对应下述技术方案实施描述；使审查员更充分了解此技术。
为了达到上述发明的目的，本发明采用的技术方案为：
提供一种基于振动信号的旋转机械设备预警方法，其包括以下步骤：
s1、获取需监测的设备历史运行的且判定为“运行正常”振动信号数据。
s2、采用小波包分解将振动信号进行分解，得到某一分解层下各频带的相对能量值构成特征矩阵。
s3、采用动态核主成分分析方法，将s2得到的特征矩阵分解为特征子空间与残差子空间。
s4、采用t2统计分析的方法处理特征子空间，求得一种表征设备健康状况的指标，t2统计量。
s5、采用基于beta分布自学习控制限构建的方法，自学习正常历史数据的控制限。
s6、将需要监测的振动信号数据采用s2步骤处理。
s7、将s6处理完成的监测振动信号数据的特征矩阵采用动态核主成分分析的方法，求得其t2统计量。
s8、若s7得到的统计量超出s5构建的自学习正常数据控制限则报警。
进一步的，步骤s2中根据小波包能量值分解得到特征矩阵的具体方法为：
将振动信号采用某一小波进行分解，最终在某一分解层数j上划分为不同频带的小波包系数
式中，j为小波包的分解层数，i∈(0,1,…,2j-1)，d(j,i)为第j层第i 1个子频带的小波包系数。
振动信号的能量被分解在各个子频带中，不同的故障特征在各个频带上的能量占比也不同，因此定义小波包相对能量为：
式中xj,i为相对能量值，反映了不同子频带的能量占比，选取某层分解后的小波包各子频带相对能量作为该信号的特征矩阵。
步骤s2中根据小波包能量值分解得到特征矩阵小波包的选择和分解层数选择的具体方法为：
在小波包分解过程中，小波形状需要根据所分析信号的特征与设备类型进行选择，对于机械设备而言，daubechies系列小波是工程上应用最广泛、最成熟的紧支集正交实小波族，简称dbn小波系(n为小波序号)。分解层数的选择与振动信号采样频率以及故障特征频率被调制到高频区间段的位置均有关系，一般工程应用分解层数不宜超过8层，一般选择3-6层。
步骤s3中采用动态核主成分分析方法，将s2得到的特征矩阵分解为特征子空间与残差子空间的具体方法为：
动态核主成分分析(dkpca)相对于传统的核主成分分析算法做出改进。dkpca适用于非线性动态过程监控方法，为了考虑时间相关性，在应用kpca之前执行数据矩阵的时滞扩展。假设某时刻下信号经特征提取后求得的特征矩阵为xt，则用前l个时刻的观测数据扩展该时刻下的样本数据扩展当前的样本数据，扩展后的动态样本振动信号的特征数据矩阵为
x＝[xtxt-1…xt-l]t(3)
式中，x为t时刻下振动信号的动态化特征矩阵，xt-1为t-1时刻下振动信号的特征矩阵。
dkpca的基本思想是动态化处理数据后，采用非线性映射的方法把输入信号映射到特征空间f中，然后在特征空间f内采用pca技术。假设某振动信号经过小波包分解得到动态能量特征矩阵xn×m，存在某变换φ，使得矩阵内某向量xi→φ(xi)，计算特征空间内n个φ(x)的样本协方差阵：
式中，φ(xi)为振动信号的特征矩阵变换，可使对振动信号特征的样本协方差内的c进行特征值分解，得到的特征值λ和特征向量v满足
λv＝cv(5)
上式两边同乘φ(xi)，得
λ(φ(xi)·v)＝(φ(xi)·cv)(6)
可以求解振动信号特征的协方差矩阵c的特征值所对应的特征向量v
式中，αi为相关系数，结合上面三个方程并构造一个n×n矩阵，kj,i＝k<φ(xj),φ(xi)>，并中心化。则：
λnα＝kα(8)
上式中的振动信号特征矩阵特征值λi(i＝1,2,…,n)及其对应的特征向量αi应满足下面约束条件：
λi(αi·αi)＝1(9)
因此，该振动信号特征矩阵的核主元的求取变为：
选取振动信号特征矩阵核主元所携带的原始特征信息量的大小是由其对特征矩阵贡献r的大小来决定的。
式中，λi为特征矩阵k的特征值，p为振动信号特征矩阵核主元数量。因此，确定某一特征矩阵贡献r后，振动信号的特征矩阵就被分解为特征子空间和残差子空间。
步骤s3中采用动态核主成分分析方法，将s2得到的特征矩阵分解为特征子空间与残差子空间时所选择的时滞参数、核函数、贡献值的具体方法为：
在动态核主成分分析中，时滞参数l的选择应根据数据采集器采样间隔和故障检测类型需求来确定，在保证考虑时序相关性的同时也不能污染本时刻数据、降低本时刻数据的特征信息含量；采用径向基核函数，且采用径向基核函数中最常用的高斯核函数做特征映射，选取高斯核函数宽度为70时，在选取特征子空间时，核主元所携带信息的大小由是由其对特征矩阵贡献r的大小来决定的，r取85％即可。
步骤s4采用t2统计分析的方法处理特征子空间，求得一种可以表征设备健康状况的指标的具体方法为：
在特征子空间中利用t2统计量来衡量核主元方法内部波动，它描述了每个采样数据在变化趋势和幅值上与给定方法的偏离程度，t2统计量的定义如下：
t2＝[t1,t2,…tp]λ-1[t1,t2,…tp]t(12)
式中，tk(k＝1,2,…,p)由式(13)确定，λ-1为与得分向量所对应的特征值构成的对角阵的逆矩阵。
步骤s5采用基于beta分布自学习控制限构建的方法，自学习正常历史数据的控制限的具体方法为：
随机变量x服从参数为α,β的beta分布写做：
x～be(α,β)(13)
形状参数α,β是决定beta分布性质的重要参数。建立自学习控制限，根据先验知识或专家选定“正常运行”工况，估计出形状参数，然后通过计算求得控制限。自学习的过程如下：
首先将“正常运行”数据做归一化处理，其次采用最大似然估计计算正常状态下的统计量数据的beta分布形状参数，然后通过确定其双侧分位数对应的阈值来确定归一化的控制限，最终反归一化求得自学习控制限。
步骤s5采用基于beta分布自学习控制限构建的方法，自学习正常历史数据的控制限的参数选择为
在自学习控制限构建过程中，外部影响在采集过程中产生的尖峰误差一般情况下为百分之五，所以双侧分位数取0.05。
步骤s7采用的方法与步骤s3、s4提到的方法相同.
本发明的有益效果为：本发明采用小波包分析、动态核主成分分析分析、t2统计分析、自学习控制限构建等技术综合形成一种旋转机械设备故障预警方法，通过大量的实验室数据以及工程案例验证，本发明提出的设备故障预警方法可以更灵敏的探测到设备的早期故障，能够降低错误报警率和漏报警率，并且又很好的泛化性。
附图说明
图1为本发明的流程示意图
图2t2统计量与自学习控制限监控图
图3频谱分析验证图，其中(a)为样本编号532点频谱图，(b)为样本编号533点频谱图，(c)为样本编号534点频谱图。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
采用美国辛辛那提大学nsfi/ucr中心轴承试验台实验数据进行举例，该实验台第二组实验轴承运行时间为2004年2月12日10:32:39-2004年2月19日06:22:39，在失效实验结束时，1号轴承外圈发现裂纹，振动信号发生剧烈变化，变化反映了1号轴承发生了快速磨损。因此提取1号轴承中的数据进行分析。
s1、选择该实验数据的前400组作为“运行正常”的信号输入方法。
s2、采用db4小波分解至第三层，求得该分解层下8个子频带的相对能量值，构成特征矩阵。
s3、采用动态核主成分分析方法，将s2得到的特征矩阵分解为特征子空间与残差子空间。
s4、采用t2统计分析的方法处理特征子空间，求得一种可以表征设备健康状况的指标，t2统计量。
s5、采用基于beta分布自学习控制限构建的方法，自学习正常历史数据的控制限。
s6、将所有数据按照采用s2步骤处理。
s7、将s6处理完成的所有的特征矩阵采用动态核主成分分析的方法，求得其t2统计量。
s8、若s7得到的统计量超出s5构建的基于beta分布自学习正常数据控制限则报警。
如图2所示，正常数据与故障数据已经发生了分离，前533组数据除刚开始运行的不平稳数据外，均处于自学习控制限(值为12.09)以下，判断534点为早期故障点。分析各编号点的频谱信息，设备在532点以前处于正常运行阶段，如图3所示，分析样本编号为532、533、534点的轴承频谱图可以发现，534、533点的频谱图中有明显的轴承故障特征频率，532点未出现故障特征频率，因此可以确定本文提出的方法算法有效检测出了轴承早期故障并实现了告警。
1.一种基于振动信号的旋转机械设备预警方法，其特征在于：包括以下步骤，
s1、获取需监测的设备历史运行的且判定为“运行正常”振动信号数据；
s2、采用小波包分解将振动信号进行分解，得到某一分解层下各频带的相对能量值构成特征矩阵；
s3、采用动态核主成分分析方法，将s2得到的特征矩阵分解为特征子空间与残差子空间；
s4、采用t2统计分析的方法处理特征子空间，求得一种表征设备健康状况的指标，t2统计量；
s5、采用基于beta分布自学习控制限构建的方法，自学习正常历史数据的控制限；
s6、将需要监测的振动信号数据采用s2步骤处理；
s7、将s6处理完成的监测振动信号数据的特征矩阵采用动态核主成分分析的方法，求得其t2统计量；
s8、若s7得到的统计量超出s5构建的自学习正常数据控制限则报警。
2.根据权利要求1所述的一种基于振动信号的旋转机械设备预警方法，其特征在于：步骤s2中根据小波包能量值分解得到特征矩阵的具体方法为，
将振动信号采用某一小波进行分解，最终在某一分解层数j上划分为不同频带的小波包系数
式中，j为小波包的分解层数，i∈(0,1,…,2j-1)，d(j,i)为第j层第i 1个子频带的小波包系数；
振动信号的能量被分解在各个子频带中，不同的故障特征在各个频带上的能量占比也不同，因此定义小波包相对能量为：
式中xj,i为相对能量值，反映了不同子频带的能量占比，选取某层分解后的小波包各子频带相对能量作为该信号的特征矩阵。
3.根据权利要求1所述的一种基于振动信号的旋转机械设备预警方法，其特征在于：步骤s2中根据小波包能量值分解得到特征矩阵小波包的选择和分解层数选择的具体方法为，
在小波包分解过程中，小波形状需要根据所分析信号的特征与设备类型进行选择，分解层数选择3-6层。
4.根据权利要求1所述的一种基于振动信号的旋转机械设备预警方法，其特征在于：步骤s3中采用动态核主成分分析方法，将s2得到的特征矩阵分解为特征子空间与残差子空间的具体方法为，
在应用kpca之前执行数据矩阵的时滞扩展；假设某时刻下信号经特征提取后求得的特征矩阵为xt，则用前l个时刻的观测数据扩展该时刻下的样本数据扩展当前的样本数据，扩展后的动态样本振动信号的特征数据矩阵为
x＝[xtxt-1…xt-l]t(3)
式中，x为t时刻下振动信号的动态化特征矩阵，xt-1为t-1时刻下振动信号的特征矩阵；
假设某振动信号经过小波包分解得到动态能量特征矩阵xn×m，存在某变换φ，使得矩阵内某向量xi→φ(xi)，计算特征空间内n个φ(x)的样本协方差阵：
式中，φ(xi)为振动信号的特征矩阵变换，使对振动信号特征的样本协方差内的c进行特征值分解，得到的特征值λ和特征向量v满足
λv＝cv(5)
上式两边同乘φ(xi)，得
λ(φ(xi)·v)＝(φ(xi)·cv)(6)
求解振动信号特征的协方差矩阵c的特征值所对应的特征向量v
式中，αi为相关系数，结合上面三个方程并构造一个n×n矩阵，kj,i＝k<φ(xj),φ(xi)>，并中心化；则：
λnα＝kα(8)
上式中的振动信号特征矩阵特征值λi(i＝1,2,...,n)及其对应的特征向量αi应满足下面约束条件：
λi(αi·αi)＝1(9)
因此，振动信号特征矩阵的核主元的求取变为：
选取振动信号特征矩阵核主元所携带的原始特征信息量的大小是由其对特征矩阵贡献r的大小来决定的；
式中，λi为特征矩阵k的特征值，p为振动信号特征矩阵核主元数量；因此，确定某一特征矩阵贡献r后，振动信号的特征矩阵就被分解为特征子空间和残差子空间。
5.根据权利要求1所述的一种基于振动信号的旋转机械设备预警方法，其特征在于：步骤s3中采用动态核主成分分析方法，将s2得到的特征矩阵分解为特征子空间与残差子空间时所选择的时滞参数、核函数、贡献值的具体方法为：在动态核主成分分析中，时滞参数l的选择应根据数据采集器采样间隔和故障检测类型需求来确定，在保证考虑时序相关性的同时也不能污染本时刻数据、降低本时刻数据的特征信息含量；采用径向基核函数，且采用径向基核函数中的高斯核函数做特征映射，选取高斯核函数宽度为70时，在选取特征子空间时，核主元所携带信息的大小由是由其对特征矩阵贡献r决定，r取85％。
6.根据权利要求1所述的一种基于振动信号的旋转机械设备预警方法，其特征在于：步骤s4采用t2统计分析的方法处理特征子空间，求得一种表征设备健康状况的指标的具体方法为：
在特征子空间中利用t2统计量来衡量核主元方法内部波动，描述每个采样数据在变化趋势和幅值上与给定方法的偏离程度，t2统计量的定义如下：
t2＝[t1,t2,...tp]λ-1[t1,t2,...tp]t(12)
式中，tk(k＝1,2,…,p)由式(13)确定，λ-1为与得分向量所对应的特征值构成的对角阵的逆矩阵。
7.根据权利要求1所述的一种基于振动信号的旋转机械设备预警方法，其特征在于：步骤s5采用基于beta分布自学习控制限构建的方法，自学习正常历史数据的控制限的具体方法为：
随机变量x服从参数为α,β的beta分布写做：
x～be(α,β)(13)
形状参数α,β是决定beta分布性质的重要参数；建立自学习控制限，根据先验知识或专家选定“正常运行”工况，估计出形状参数，然后通过计算求得控制限。
8.根据权利要求7所述的一种基于振动信号的旋转机械设备预警方法，其特征在于：
自学习的过程如下：
将“正常运行”数据做归一化处理，采用最大似然估计计算正常状态下的统计量数据的beta分布形状参数，通过确定其双侧分位数对应的阈值来确定归一化的控制限，反归一化求得自学习控制限。
9.根据权利要求7所述的一种基于振动信号的旋转机械设备预警方法，其特征在于：
步骤s5采用基于beta分布自学习控制限构建的方法，自学习正常历史数据的控制限的参数选择：
在自学习控制限构建过程中，外部影响在采集过程中产生的尖峰误差为百分之五，所以双侧分位数取0.05。
技术总结