基于动态自学习阈值和趋势滤波的机械故障智能预警方法
2014-05-17江志农
张 明,冯 坤,江志农
(北京化工大学诊断与自愈工程研究中心,北京100029)
基于动态自学习阈值和趋势滤波的机械故障智能预警方法
张 明,冯 坤,江志农
(北京化工大学诊断与自愈工程研究中心,北京100029)
针对当前机械在线监测系统报警难以实现机械故障早期预警问题,提出一种智能预警方法。基于在线监测系统大量监测数据统计分析,采用动态的自学习阈值算法计算预警阈值,并应用l1趋势滤波技术消除随机误差获取滤波后的趋势。应用动态自学习阈值替代监测系统中的常规报警阈值,比较自学习预警阈值与滤波后的趋势,实现了机械故障早期预警。工程实例表明,该方法能够对机械故障实现早期预警,对预防机械事故的发生有重要的作用。
自学习阈值;故障预警;非参数检验;beta分布;l1趋势滤波
汽轮机、发电机、压缩机等是石化、电力等流程工业中广泛使用的机械,该类设备安全、稳定的运行,会产生良好的经济效益和社会效益。目前,企业中的大型关键设备多数已安装机械在线监测系统,但是,当前的机械在线监测系统并不能实现机械故障的早期预警,主要因为存在以下问题:①报警阈值由主机厂提供并预先设定到监测系统中,当机组报警时故障已经恶化到一定程度,并不能在故障发生的早期实现预警;②若为了使监测系统实现早期预警而将报警阈值调低,则可能会因噪声及采集误差的影响使实时采集的监测数据反复穿越报警线,导致大量错误预警;③若采用传统平滑滤波技术消除噪声及采集误差的影响,则有可能丢失关键故障信息,从而产生严重事故,给企业带来巨大损失。
研究能够充分利用机械在线监测系统的监测参数数据解决上述问题的机械智能预警技术[1]是现实的迫切需要,具有重要的工程意义,同时对推进机械在线监测系统走向智能化具有理论指导意义。
1 智能预警方法的提出
机械在线监测系统当前常规报警方式为:预先设定测点的报警值与危险值(固定值,参考API、ISO等标准),将实时采集的特征值数据与预先设定的报警值与危险值进行比较,若实时数据小于报警值,则判定机组运行正常;若实时数据大于报警值且小于危险值,则判定机组发生故障,此时需要排查故障原因并判定是否需要停车维修;当实时数据大于危险值,则需要立即停车维修。
如图1所示,机械在线监测系统常规报警方式十分简单,未能利用监测参数趋势数据,当超过报警值时机械故障往往已经恶化,研究能够充分利用监测参数趋势数据的智能预警技术是非常必要的。国内外诸多学者对预测故障的智能报警技术进行了研究,Dokas等[2]研究了基于故障树分析和模糊专家系统的预警诊断技术;Widarsson等[3]研究了回收锅炉泄漏故障的贝叶斯神经网络预警诊断方法;Zhao等[4]研究了市政固体废物焚化炉在线早期故障监测与诊断方法;马晋等[5]研究了基于瞬时转速波动率实现了内燃机点火状态的早期预警方法;魏中青等[6]研究了快变报警技术实现了大型旋转机械突发性故障的早期预警诊断;张栋梁等[7]研究了状态监测系统中的智能预警管理系统。
图1 常规报警判断流程图Fig.1 General alarm flow chart
本文基于监测参数趋势数据统计分析提出了自学习预警阈值与趋势滤波相结合的机械故障智能预警方法。研究了动态自学习预警阈值算法,实现监测参数预警阈值的自动学习;使用Kim等提出的l1趋势滤波技术[8]对监测参数趋势数据进行滤波,消除了趋势数据中的采集误差。
2 监测参数趋势数据的统计分析
选取了卡方检验(Chi-square test)、Kolmogorov-Smirnov检验及秩和检验,对机械在线监测系统不同测点的监测参数趋势数据进行了非参数检验[9]。结果表明:平稳状态下监测参数趋势趋数据不服从高斯分布,并且同一机组的不同测点监测数据也不能认为来自同一总体。
2.1 卡方检验
假设一个定性变量Y具有k个可能取值或有k种分类(标为1,2,…,k),Y的概率分布自然地由概率函数P(Y=i)(i=1,2,…,k)所确定。现在要考查已观察到的一组样本(容量为n)与某确定的分布G拟合的程度,相当于研究P(Y=i)(i=1,2,…,k)与G之间的差异,看这个差异是否属于偶然变异,这就是卡方检验的基本原理。卡方检验的步骤如下:
其中:Oi表示落入第i组的样本观测值的实际频数,表示理论频数。
(3)作出判断
2.2 Kolmogorov-Sm irnov检验
Kolmogorov-Smirnov检验用于检验一组样本观测结果的经验分布同某一指定的理论分布(如高斯分布)之间是否一致。K-S检验的基本思路为,将顺序分类数据的理论累计频率分布同观测的经验累计频率分布加以比较,求出他们最大的偏离值,然后在给定的显著性水平上检验这种偏离值是否是偶然出现的。
设理论累计频数分布为F(x),n次观测的随机样本的经验分布函数,K-S检验的步骤如下:
(1)零假设H0:经验分布与理论分布没有显著差别。
(2)把样本观测值从小到大排列为:
x(1),x(2),…,x(n),计算经验累积分布函数Fn(x)和理论累积分布函数F(x)。
双侧显著性水平(P值)根据Smirnov提出的公式计算。
(3)作出判断
若P<α,则拒绝零假设。
2.3 秩和检验
对于分别抽自总体X和Y的两个独立样本,把两者混合在一起,并按其量值由小到大顺序排列,得:
如果xk=zj,则记Rk(x)≜j,则称为xk在混合样本中的秩。
设Y的样本容量n2小于X的样本容量n1,威尔柯克逊提出用Y在混合样本中的秩和总和:
来检验X和Y的分布是否一致。秩和检验的步骤如下:
(1)零假设H0:两个独立样本来自相同的总体。
(2)计算相应检验统计量值或P值。
(3)作出判断
若P<α,则拒绝零假设。
选取同一机组不同测点的监测参数趋势数据1和监测参数趋势数据2,如图2所示,进行卡方检验、Kolmogorov-Smirnov检验及秩和检验,数据单位为μm。
对监测参数趋势数据1进行卡方检验和Kolmogorov-Smirnov检验:
(1)零假设:监测参数趋势数据1服从高斯分布;
(2)卡方检验计算P1=2.595×10-4,Kolmogorov-Smirnov检验计算P2=7.166 4×10-307;
(3)在显著性水平α=0.05时,P1<α,P2<α,则两种检验均拒绝零假设,因此监测参数趋势数据1不服从高斯分布。
同理,对监测参数趋势数据2进行卡方检验和Kolmogorov-Smirnov检验:
(1)零假设:监测参数趋势数据2服从高斯分布;
(2)卡方检验计算P1=6.095 7×10-8,Kolmogorov-Smirnov检验计算p2=7.166 4×10307;
图2 (a)监测参数趋势数据1(b)监测参数趋势数据2Fig.2(a)Monitoring parameters trend data 1(b)Monitoring parameters trend data 2
(3)在显著性水平α=0.05时,P1<α,P2<α,则两种检验均拒绝零假设,因此监测参数趋势数据2也不服从高斯分布。
对监测数据参数数据1与数据2进行秩和检验:
(1)零假设:监测参数趋势数据1与数据2来自相同的总体;
(2)秩和检验计算P=1.042 3×10-115;
(3)在显著性水平α=0.05时,P<α,则拒绝零假设,监测参数趋势数据1和数据2来自不同的总体。
3 基于贝塔分布的预警阈值智能学习方法
近年来国内外对故障诊断自学习报警阈值方法做了诸多研究,Alkaya等[10]研究了基于PCA模型的方差敏感自学习阈值的方法及应用;Raka等[11]研究了基于鲁棒自学习阈值的故障监测方法;Shi等[12]做了基于模型的自学习阈值故障监测方法的进一步研究;魏中青等[13]针对往复压缩机气阀早期故障预警问题研究了MLE(最大似然估计)阈值规则。
由于监测参数趋势数据分布不服从高斯分布,并且不同测点的趋势数据来自不同的总体,说明了不同测点的趋势数据分布形式不同。基于统计学分析对不同分布特性研究发现:Beta分布是一种多参数统计分布,可以通过调整其多参数逼近任意分布形式,使用贝塔分布分别逼近正常状态下各测点监测参数趋势数据的概率密度分布,可得到各测点趋势数据的贝塔概率密度函数,计算出各测点趋势数据的预警阈值区间。
3.1 贝塔分布
式中:μx为x的均值;σ2x为x的方差。
当贝塔分布的上限a=0,下限b=1时,即成为标准贝塔分布:
贝塔分布密度函数的形状可由γ,η来控制,其形状从均匀分布到近似高斯分布,从对称到不对称,而且分布在一个有限的范围内。
3.2 贝塔分布的参数拟合算法
机械在线监测系统的监测参数趋势数据的概率密度可以通过调整贝塔分布的形状参数γ,η进行拟合,本文采用最小二乘法进行拟合[14]。
(1)估计概率密度函数
其中:Xi是将样本数据分区间后,每个区间的中点值;mi是样本数据落在第i个区间内的个数;n是样本数据的个数;hi是每个区间的半宽度;N是所有区间的个数。
(2)利用最小二乘法求解形状参数
其中:Wi为权系数,可根据测量数据情况进行选取;f(Xi,γ,η)为标准贝塔分布的概率密度函数。
3.3 监测参数趋势数据的预警阈值智能学习算法
在线监测参数数据一般情况下不在0到1的范围内,所以使用标准贝塔分布拟合监测参数数据的概率密度分布时需要对振动趋势数据进行归一化处理。
根据机组平稳运行一段时间后得到的监测参数趋势数据,计算得到监测参数趋势数据的概率密度分布,通过最小二乘法拟合求出最逼近实际趋势数据的贝塔分布形状参数γ,η,得到趋势数据的贝塔分布概率密度函数β(γ,η)。预警阈值自学习过程如下:
(1)将一段时间内的监测参数趋势数据Xi归一化处理:
Xi为监测参数趋势数据;Xmin为监测参数趋势数据中的最小值;Xmax为监测参数趋势数据中的最大值;X^i为归一化后的监测参数趋势数据;
(2)依据式(9)计算监测参数趋势数据的概率密度分布;
(3)使用上文提到的最小二乘法拟合监测参数趋势数据的概率密度分布得到标准贝塔分布的两个形状参数γ,η;
α在自学习预警阈值的应用中定义为尖峰噪声引起的采集误差。外部影响导致在采集过程中产生的尖峰噪声的误差一般情况下为5%[18],即α=0.05。在实际应用过程中,可根据采集系统的性能以及传感器特性调整α值。本文在计算自学习预警阈值时,α取0.05。
选取一组旋转机械振动趋势数据,进行预警阈值智能自学习,数据单位为μm(如图3所示)。按照自学习预警阈值计算方法,取α=0.05时,使用最小二乘法拟合监测参数趋势数据的概率密度分布(如图4(a)所示),得到拟合后的贝塔分布(如图4(b)所示),贝塔分布两个形状参数分别为γ=34.062 8、η=46.565 0,计算趋势数据的阈值下限Threshol d1=23.009 5、上限Threshol d2=24.853 8,自学习报警阈值空间[23.009 5,24.853 8]。
图3 监测参数趋势数据Fig.3 Monitoring parameters trend data
4 基于l1趋势滤波技术的机械故障智能预警方法
4.1 趋势滤波
一个标准的时间序列y1,t=1,…,n,假设是由一个基本趋势xt和一个随机变量zt组成。趋势滤波就是从标准的时间序列中估计出基本趋势xt,或者估计出随
图4 (a)趋势数据概率密度分布(b)贝塔概率密度分布Fig.4(a)Trend data probability density distribution(b)beta probability density distribution
机变量zt=yt-xt。趋势滤波实际上就是估计一个标准时间序列的基本趋势xt的过程。
趋势滤波方法包括:滑动平均滤波[15],指数平滑滤波[16],和Hodrick-Prescott(H-P)滤波[17]及l1趋势滤波[8]等。H-P滤波作为趋势滤波方法广泛地应用在许多领域,l1趋势滤波是在H-P滤波基础上发展来的,H-P滤波和l1趋势滤波各有优势,l1趋势滤波更适用于分段线性的趋势的估计。
4.2 l1趋势滤波
l1趋势滤波[8]是Kim等在近几年提出的,是一种H-P滤波的变异,这种改变是使用绝对值(l1范数)和替代了应用在H-P滤波估计趋势中惩罚变量的平方和。l1趋势滤波方法和H-P滤波趋势滤波方法有很多共同的性质,同时还有共同的计算复杂性。两者根本性不同在于,l1趋势滤波得到的估计趋势在分段的线性范围内更平滑。因此,l1趋势滤波特别适合分析滤波含有分段线性时间序列的趋势。实际机械设备故障特征往往表现为趋势的奇点、斜率的改变等变化,因此l1趋势滤波恰恰可以检测和估计这些潜在的变化。
机械在线监测系统中监测参数趋势数据总会存在一定采集误差,它会影响机械故障早期预警的准确性。在线监测系统监测参数趋势数据主要有两种采集误差:一种是连续的高斯白噪声,这种误差无法绝对避免只能尽量减小;另一种是尖峰噪声,这种误差可能是由于系统及传感器问题造成[18]。使用l1趋势滤波对监测参数趋势数据进行滤波处理,能够很好地得到机组真实的振动趋势,同时消除了采集误差对机械故障早期预警的影响。
4.2.1 l1趋势滤波算法
l1趋势滤波算法是在H-P滤波基础上进行改进的新型滤波算法。这种趋势估计是通过最小化加权目标函数实现的,加权目标函数见式(15):
λ是一个非负参数用来控制x的平滑性和平衡余项的大小。加权目标函数对x来说是一个严格的凸函数,所以只有一个最小值,用xlt表示,因此xlt就是滤波后的趋势。
文献[8]指出,l1趋势滤波求解问题可以等同于对正则化l1最小二乘求解问题,
通过最小二乘法求得此问题结果θlt,则l1趋势滤波结果xlt=Aθlt。
4.2.2 l1趋势滤波仿真研究
模拟叶片脱落故障原始趋势如图5(a)所示,在原始趋势的基础上添加噪声如图5(b)所示,分别使用l1趋势滤波和H-P趋势滤波方法对叶片脱落故障的模拟信号进行滤波。调整l1趋势滤波和H-P趋势滤波控制系数λ,保证拟合误差(即:时间序列余项zt=ytxt二范数的均值)相等的情况下,得到图5所示滤波结果,拟合误差为0.141 9。对比图5中(c)、(d)两图,看出l1趋势滤波能够更好地估计出实际机械设备突发型故障特征趋势。
5 本文方法在线监测系统中的工程应用研究
机械智能预警方法将动态自学习预警阈值算法与l1趋势滤波技术相结合,具体步骤如下:
(1)依据动态自学习预警阈值算法计算平稳运行阶段监测参数趋势数据的动态自学习预警阈值。
(2)使用l1趋势滤波技术对监测参数趋势数据进行滤波处理。
图5 (a)叶片脱落原始趋势(b)叶片脱落原始趋势+噪声(c)l1趋势滤波后的趋势(d)H-P趋势滤波后的趋势Fig.5(a)Original trend(b)Original trend+noise(c)l1 filtering trend(d)H-P filtering trend
(3)将实时滤波后的监测参数趋势数据与计算得到的动态自学习预警阈值进行比较,如果趋势数据超出阈值范围则预警。
下面选取国内某石化企业透平压缩机在线监测系统中缓变故障案例和突变故障案例验证本文方法。
5.1 工程案例一
某石化企业烟气轮机的常规报警值为80μm、危险值为100μm。该机组某测点监测趋势从2012年6月初开始缓慢爬升,在此之前该测点趋势一直平稳在73 μm左右,到2012年7月初爬升到95μm左右,一个月内趋势异常爬升了22μm。排查故障原因发现烟气中含有了催化剂,使得烟机叶片上逐渐粘附催化剂导致不平衡量逐渐增大,导致趋势发生缓慢爬升。
根据动态自学习预警阈值算法计算正常运行阶段的动态自学习预警阈值,取α=0.05时,得到拟合后的贝塔分布的两个形状参数分别为γ=2.725 6、η=2.419 1,趋势数据的预警阈值下限Threshol d1=70.790 5、上限Threshol d2=73.698 5,自学习预警阈值空间[70.790 5,73.698 5]。使用l1趋势滤波技术对趋势数据进行滤波处理。当该机组实时滤波后的趋势数据超出动态自学习预警阈值空间后发出预警。
如图6所示,使用常规报警方式发出报警时,故障已经恶化到一定程度。与常规报警方式相比,机械智能预警方法在缓变故障发生的早期就可以发出预警,实现了缓慢变化故障早期预警,弥补常规报警方式在缓变故障预警上的不足。
5.2 工程案例二
某石化企业二氧化碳压缩机的常规报警值为60 μm,危险值为90μm。该机组在正常运行过程中突然发生异常,在5 s的时间里振动趋势从40μm左右跳变到50μm左右。事后排查故障原因发现压缩介质不干净,在长期运行过程中污垢在叶片上积累,污垢脱落瞬间导致转子不平衡量突然变化,从而导致趋势突变。
根据动态自学习预警阈值算法计算正常运行阶段的动态自学习预警阈值,取α=0.05时,得到拟合后的贝塔分布的两个形状参数分别为γ=1.918 1、η=2.642 5,趋势数据的预警阈值下限Threshol d1=34.953 2、上限Threshol d2=41.717 7,自学习预警阈值空间[34.953 2,41.717 7]。使用l1趋势滤波技术对趋势数据进行滤波处理。当该机组实时滤波后的趋势数据超出动态自学习预警阈值空间后发出预警。
图6 (a)缓变故障原始趋势(b)l1滤波后缓变故障趋势Fig.6(a)Original trend of slowly changed fault(b)l1 filtering trend of slowly changed fault
图7 (a)突变故障原始趋势(b)l1趋势滤波后突变故障趋势Fig.7(a)Original trend of fast changed fault(b)l1 filtering trend of fast changed fault
如图7所示,本例突变故障并没有触发常规报警。与常规报警相比,机械智能预警方法在突变故障发生的早期就可以发出预警,实现了突然变化故障早期预警,弥补常规报警方式在突变故障预警上的不足。
6 结 论
本文提出基于自学习预警阈值和趋势滤波的机械故障智能预警方法,得到如下结论:
(1)对监测参数趋势数据进行非参数检验,验证了趋势数据并不服从高斯分布并且不同测点的趋势数据来自不同总体。
(2)使用贝塔分布拟合不同测点监测参数趋势数据的概率密度分布,并提出了动态自学习预警阈值算法。利用l1趋势滤波技术消除了监测参数趋势数据中的采集误差。
(3)提出机械智能预警方法为:将实时l1趋势滤波后的监测参数趋势数据与动态自学习预警阈值进行比较,如果趋势数据超出阈值范围则预警。
工程案例表明:机械智能预警方法能够弥补常规报警方式在早期故障预警方面的不足,在故障发生的早期尚未触发常规报警时就可以监测到趋势异常从而实现早期预警。
[1]高金吉.机械故障诊治与自愈化[M].北京:高等教育出版社,2012.
[2]Dokas IM,Karras D A,Panagiotakopoulos D C.Fault tree analysis and fuzzy expert systems:Early warning and emergency response of landfill operations[J].Environmental Modelling&Software,2009,24(1):8-25.
[3]Widarsson B,Dotzauer E.Bayesian network-based earlywarning for leakage in recovery boilers[J].Applied Thermal Engineering,2008,28(7):754-760.
[4]Zhao JS,Huang JC,Sun W.On-line early fault detection and diagnosis of municipal solid waste incinerators[J].Waste Management,2008,28(11):2406-2414.
[5]马晋,江志农,高金吉.基于瞬时转速波动率的内燃机故障诊断方法研究[J].振动与冲击,2012,31(13):119-124.
MA Jin,JIANG Zhi-nong,GAO Jin-ji.Diesel engine fault diagnosis method based on instantaneous angular speed fluctuation ratio[J].Journalof Vibration and Shock,2012,31(13):119-124.
[6]魏中青,马波,江志农,等.旋转机械振动监测快变报警系统研究及工程应用[J].机电工程技术,2010,39(12):82-84,114.
WEIZhong-qing,MA Bo,JIANG Zhi-nong,et al.Rotating machinery vibration monitoring fast changing alarm systems research and engineering applications[J].Mechanical&Electrical Engineering Technology,2010,39(12):82- 84,114.
[7]张栋梁,魏中青.基于状态监测的预警管理系统研究与应用[J].电子设计工程,2011,12:8-10,14.
ZHANG Dong-liang,WEIZhong-qing.Study and application of early warning management system based on condition monitoring[J].Electronic Design Engineering,2011,12:8-10,14.
[8]Kim S J,Koh K,Boyd B,et al.l1 Trend filtering[J].SIAMReview,Problems and Techniques Section,2009,51(2):339-360.
[9]宋昕,蔡泳,徐刚,等.非参数检验方法概述[J].上海口腔医学,2004,13(6):561-563.
SONG Xin,CAI Yong,XU Gang,et al.A review of nonparametric test[J].Shanghai Journal of Stomatology,2004,13(6):561-563.
[10]Alkaya A,Eker I.Variance sensitive adaptive thresholdbased PCA method for fault detection with experimental application[J].ISA Transactions,2011,50(2):287-302.
[11]Raka S A,Combastel C.Fault detection based on robust adaptive thresholds:A dynamic interval approach[J].Annual Reviews in Control,2013,37(1):119-128.
[12]Shi Z,Gu F,Lennox B,et al.The development of an adaptive threshold for model-based fault detection of a nonlinear electro-hydraulic system[J].Control Engineering Practice,2005,13(11):1357-1367.
[13]魏中青,马波,窦远,等.基于MLE阈值规则的小波特征提取技术在气阀故障诊断中的应用[J].振动与冲击,2011,30(1):237-241.
WEIZhong-qing,MABo,DOU Yuan,etal.Wavelet feature extracting technique based on maximum likelihood estimation threshold rule and its application in fault diagnosis of a gas valve[J].Journal of Vibration and Shock,2011,30(1):237-241.
[14]曹天捷,周则恭.最小二乘法在估计概率分布参数中的应用[J].机械强度,1997(4):34-37.
CAO Tian-jie,ZHOU Ze-gong.Application of the least square method in the estimation of the parameters of the probability distribution[J].Journal of mechanical strength,1997(4):34-37.
[15]Jose A R,Eduardo R,Juan C E.Detrending fluctuation analysis based on moving average filtering[J].Physica A:Statistical Mechanics and its Applications,2005,354:199-219.
[16]Taylor J W.Exponential smoothing with a damped multiplicative trend[J].International Journal of Forecasting,2003,19(4):715-725.
[17]Hodrick JR,Prescott C E,Postwar U S.Business cycles:An empirical investigation[J].Journal of Money,Credit and Banding,1997,29(1):1-29.
[18]Chen Z,Ivanov P C,Hu K,et al.Effect of nonstationarities on detrended fluctuation analysis[J].Physical Review E,2002,65(4):041107.
A mechanical fault early warning methodology based on dynamic self-learning threshold and trend filtering techniques
ZHANGMing,FENG Kun,JIANG Zhi-nong
(Diganosis&Self-Recovery Engineering Research Center,Beijing University of Chemical Technology,Beijing 100029,China)
Because the alarm mode of current online monitoring systems is hard to realize early-warning a mechanical fault,an early warning methodology was proposed here.Based on statistical analysis ofmass data in online monitoring systems,advantages of the dynamic self learning threshold algorithm were taken to calculate a warning threshold,and the l1 trend filtering technology was used to gain the filtered trend eliminating random errors.With dynamic self-learning threshold instead of general alarm threshold in monitoring systems,early warning of a fault could be acquired by comparing the self-learning warning threshold and the filtered trend.Itwas shown from engineering examples that early warning of amechanical fault can be achieved with the proposed method,it plays an important role in preventing occurrence ofmechanical fault.
self-learning threshold;fault early warning;nonparametric test;beta distribution;l1 trend filtering
TH165+.3;TP391
A
10.13465/j.cnki.jvs.2014.24.002
国家自然科学重点基金项目(51135001);国家青年科学基金项目(51305020);国家重点基础研究发展计划(“973”计划)项目(2012CB026000)
2013-10-21 修改稿收到日期:2014-01-02
张明男,博士生,1988年7月生
江志农男,教授,博士生导师,1967年10月生