向量自回归和T2控制图在变工况齿轮箱故障检测中的应用

薛 成,李 鑫,2,焦 静,毛雨琳

(1.南京工程学院 汽车与轨道交通学院, 南京 211167;2.南京航空航天大学 民航学院, 南京 211106)

0 引言

齿轮箱以其结构紧凑、传动比高、承载能力大的独特特点被广泛地集成到各种旋转机械中。由于工作环境恶劣,齿轮箱容易损坏,导致整个机械系统性能下降甚至停机。为了有效地提高机器运行的可靠性和降低维护成本,对齿轮箱进行在线监测和故障诊断是十分必要的[1-2]。

基于振动信号的分析是最常见的齿轮箱故障诊断方法,但在面对多工况交替等情况时,其信号波动幅度大,各个工况的振幅不同,导致监测诊断能力有所下降。然而,振动信号的特征阶次在非平稳工况下,不会随着转速的变化而发生改变。因此,应用阶次分析的方法可以有效进行非平稳工况下齿轮箱的故障诊断。赵登利等[3]推导并总结了齿轮箱特征阶次,并通过真实风场数据进行了故障诊断分析与实验验证。王博等[4]根据齿轮转速对故障频率的影响,总结了齿轮箱各部件的故障阶次特征表,实现了变转速工况下的齿轮箱故障诊断。唐敏[5]也运用阶次跟踪技术将变速齿轮箱的时域非平稳特征转化为平稳的瞬态特征,优化了变速齿轮箱的状态评估效率。杜丽君等[6]提出了一种时域稀疏调频字典的齿轮复合故障诊断方法,从变速齿轮箱复杂的信号中准确提取出故障信号分量,确保了后续信号的幅值重构精度,完成了复合故障信号的分离诊断。张小刚等[7]提出了一种深度残差网络的方法,能够更加有效地挖掘信号特征之间的信息,提高诊断精度。

程龙欢等[8]则认为,单源的振动信号不足以体现设备的运行状态,应该将多个振动传感器进行组合使用,融合多个传感器的振动信号,才能够正确表明设备的状态。为了解决单一故障特征诊断准确性低的问题,张智禹等[9]提出了融合注意力机制的改进深度置信网络的变工况齿轮箱故障诊断方法,以融合时域、频域、小波包时频域特征作为故障诊断参数,有效提高了故障诊断精确度和泛化能力。Kumar等[10]将多传感器信息融合技术应用在齿轮箱故障诊断中,使用振动加速度计、麦克风和声发射传感器从齿轮箱中获取不同加载工况和转速变化下的机器状态数据,对采集到的振动、声音和声发射信号进行小波特征提取,识别出突出特征,并采用机器学习算法以实现故障诊断过程的自动化。

统计过程控制可以有效监测数据的随机波动和异常波动状况。马晓原[11]运用统计过程控制方法对数控机床加工单零件的能耗进行控制,实现了单零件能耗的提前预警。梁坤等[12]基于统计过程方法,采用指数加权移动平均(exponential weighted moving average,EWMA)控制图,提高了故障预警的准确度。

本研究将时间序列模型与统计过程控制方法相结合,进行了变工况下齿轮箱的异常检测。在时域同步平均算法(time synchronous averaging,TSA)预处理振动信号的基础上,选取健康数据建立向量自回归模型(vector autoregression,VAR),提取特征参数,从而建立统计过程控制图,监测数据的波动情况,评估提取的故障特征值,确定齿轮箱在正常工作状态下故障特征值的幅值,准确预测、判断齿轮箱异常发生点,及时地进行故障预警。

1 数据获取

1.1 数据来源

试验数据来源于如图1所示的宾夕法尼亚大学的机械故障诊断试验台[13]。通过此试验台读取转速与扭矩数据,可对齿轮箱寿命进行可靠性验证试验。

图1 试验台设备图

试验所用的齿轮箱为14号齿轮箱,该齿轮箱额定输入转速为1 750 r/min,最大额定输出扭矩为555 in-lbs(1 in-lbs=0.113 N·m),最大额定输入功率为4.66 hp(1 hp=745.7 W),包含一个齿数为70的驱动齿轮和一个齿数为21的小齿轮,齿轮啮合频率为613 Hz。齿轮箱振动数据来源于三轴加速度传感器。测试的振动数据以离散的间隔采样,每个样本被保存为一个文件[14]。

图2是从A10、A11以及A12三轴加速度计中采集的齿轮箱原始的振动数据曲线,其中A10设置为轴向,A11与地面方向平行,而A12则与地面方向垂直[15]。观察原始振动数据,发现数据点的数量十分巨大,前96 h,输出扭矩恒定在100%时,振幅稳定,但工况开始变化以后,振幅变化不明显,无法提取有用的信号来进行异常检测。

图2 整个试验过程的原始振动数据曲线

1.2 采样过程

齿轮箱全寿命周期试验在不同工况时分别采样一次,采样频率为20 kHz,每次采样时间为10 s。试验在2种工况下进行:转速恒定不变,齿轮箱在100%额定扭矩下运行96 h;然后以50%、100%、150%、200%、250%、300%扭矩和功率运行了19.3 h,直到出现故障,失效停机。本试验在100%额定扭矩工况下,齿轮箱正常工作。在变工况状态下,齿轮箱出现异常直至故障停机,所以只需选取变工况下齿轮箱的振动信号进行分析。TR#14齿轮箱在不同负载扭矩下共采样152个文件。将故障停机的齿轮箱拆开进行检查,检查发现小齿轮正常,驱动齿轮有5个完全损坏的齿和2个部分损坏的齿。

2 数据预处理

2.1 TSA

TSA算法选取信号时应按照齿轮的旋转周期进行,提取齿轮旋转的旋转频率与啮合频率,然后将同步相加、平均[16]。随着同步平均次数增加,噪声干扰信号逐渐减少,直到不足以对故障诊断产生影响。假设每个振动数据文件中有n个数据点。研究齿轮旋转一周的完整的采样点数由以下公式求得:

(1)

式中:fs为取样频率,fm为研究齿轮的基本啮合频率,Nt为齿轮的齿数,「·⎤表示对计算结果向上取整。

需要平均的循环次数可以由下公式求得:

(2)

式中:N为每个数据文件中采样点的总数,⎣·」为对计算结果向下取整。

时域信号中,振动数据中的TSA信号可由下式计算:

(3)

式中:K为齿轮每运行一个周期的采样点数,M0为单个采样文件所包含的齿轮的转数。

将原始振动信号导入Matlab软件中,利用TSA算法得到TSA处理曲线,如图3所示。通过TSA信号的健康数据可以建立VAR模型来拟合数据,再将原始数据与拟合数据相减得到残差。

图3 变工况下A10、A11、A12 TSA信号

2.2 分工况处理

在对原始振动信号进行TSA处理之后,由于不同工况的扭矩负载不同,因此需要对数据进一步按不同的工况进行分类。文件序号对应的输出扭矩值如图4所示。

图4 齿轮箱输出扭矩曲线

将相同工况的文件中的振动数据组合,即可得到各个工况下的TSA振动数据,其中比较典型的A12分工况处理后的TSA振动数据如图5所示。

图5 分工况TSA处理后A12的振动数据曲线

2.3 VAR模型及残差计算

VAR模型是一种常用的计量经济模型,常使用模型中所有当期变量对所有变量的若干滞后变量进行回归。VAR模型可以在没有预先限制的情况下,对联合内生变量进行动态预测,实现对2个或多个相关联的时间序列的预测。VAR模型描述在同一样本期间内的n个变量(内生变量)可以作为它们过去值的线性函数[17]。一个VAR(p)模型可以写成:

yt=c+A1yt-1+A2yt-2+…+Apyt-p+εt

(4)

式中:c为常数向量,yt为n维内生变量,Ai(i=1,2,…,p)为n×n矩阵,εt为误差向量。

选取经TSA处理后健康状态下的数据作为训练数据,使用Eviews软件建立VAR模型。因为齿轮箱大部分时间处于正常情况下运行,为减小计算量,只需选取一部分平稳数据进行建模。选取前5组数据TSA信号建立VAR模型。

数据序列经过检验后证明平稳。有代表性的300%工况下建立的VAR模型的阶数确定如表1所示,根据各个信息准则的参数,综合考虑确定该VAR模型的滞后阶数为47。

表1 300%工况建立的VAR模型的阶数

通过单位圆检验对该模型进行稳定性判断,如图6所示,所有的点均处于单位圆内,因此判断该VAR模型稳定。

图6 300%工况建立的VAR模型稳定性判断

模型的表达式如下:

A10=1.334 969 002 31×A10(-1)-

0.563 750 025 639×A10(-2) +…+

0.001 687 210 172 87×

A12(-47)-1.718 841 198 06

A11=0.088 596 424 477×A10(-1)-

0.163 006 450 186×A10(-2) +…-

0.007 271 001 622 06×

A12(-47)+6.807 399 845 83

A12=-0.065 850 411 545×A10(-1) +

0.061 530 633 529 9×A10(-2) +…+

0.079 860 198 146 3×

A12(-47)-1.704 486 071 05

分析残差需要选取全部的文件,通过VAR模型拟合的数据与原始数据进行对比,整体分析判断齿轮箱的健康状态,得到的残差图如图7所示。

图7 残差曲线

2.4 特征值提取

常用的时域信号的有量纲特征值有最大值、最小值、峰值、整流平均值、均值、方差、均方根值;无量纲特征值有值因子、脉冲因子、裕度因子、峭度因子、波形因子和偏度等。

整流平均值可表示为信号绝对值积分的平均值:

(5)

均方根(root mean square,RMS)又称有效值。将所有数据平方求和,求其均值,再开平方,就得到均方根值,其计算公式如下:

(6)

有量纲指标尽管对信号特性非常敏感,但是在工作条件(例如负荷)变化时也会发生改变,而且容易受到外界的干扰,因此其性能不稳定。相对来说,无量纲指标可以消除干扰,因而在特征抽取方面得到了广泛的应用。

峭度因子是四阶中心矩和标准差的四次方的比值,也是对振动信号冲击特性的反应。偏度和峭度有一定的相关性,偏度因子是三阶中心矩和标准差的三次方的比值,描述的是分布。

为了确保提取的特征值能够较为准确地反映出齿轮箱的运行状态,选取了有量纲特征值整流平均值、均方根以及无量纲特征值偏度、峭度因子进行比较。以齿轮箱变工况运行时A10的振动数据为例,提取了其VAR模型与TSA数据残差的特征值,如图8—11所示。

图8 整流平均值曲线

齿轮箱在变工况状态下运行,其振动数据应具有一定的规律性,同一种工况时的振动数据应大致相等或随时间而略有增加,振动数据应随着扭矩增加或减小而相应的增加或减小。根据上述判断,观察选取符合实际情况的特征值图像。观察图8数据,曲线具有规律性,随时间稳定上升,符合推断。图9曲线与图8曲线特征相似,符合推断。

图9 均方根曲线

图10曲线波动幅度基本相同,不符合推断。图11曲线复杂无规律,难以找到特征。

图10 偏度因子

图11 峭度因子

因此,选用整流平均值以及均方根这2个特征指标对数据进行提取。

残差是实际数据与拟合数据的差值,采用时域分析将残差数据中的故障特征提取可以进一步实现对齿轮箱的故障诊断。通过特征值的波动情况可以大概判断齿轮箱的异常情况。利用时域分析计算残差的整流平均值(图12)和均方根(图13)指标来初步判断齿轮的故障发生点。

图12 残差的整流平均值

图13 残差的均方根值

通过对图12、图13的观察,发现齿轮箱在3个方向上的振动信号基本稳定,但在曲线末端都有明显的上升,曲线十分陡峭,可以确定有故障出现,但在其他位置处是否产生故障无法判断,接下来进行统计过程分析,从而有效地找到异常点。

3 统计过程控制

控制图是通过测量、记录和评价工艺质量特性,从而监测工艺是否在控制中的一种图表。它可以区别正常和不正常的波动,是一种重要的统计手段。在这张图表上,有3条与横向轴线平行的线:中心线(central line,CL)、上控制限(upper control limit,UCL)、下控制限(lower control limit,LCL)以及按时间顺序提取的采样统计数字的描点序列。UCL、CL和LCL共同称为控制极限,一般将控制极限设置在±3个标准偏差处。CL是统计数据的平均值,其上、下控制边界和CL之间有几个标准偏差。如果在控制图中的描点处在UCL和LCL以外,或者在UCL和LCL之间没有任意排列,那么说明该过程是不正常的[17]。

EWMA控制图可以检测过程均值中的小偏差,而不会受到低值和高值的影响。由图3中TSA信号与图7中残差值可以明显得出,A10振动信号振幅大于A11与A12,波动强于A11与A12。选取A10的VAR模型与TSA残差的整流平均值与均方根数据,绘制EWMA控制图。在EWMA控制图选项中使用默认设定,即距离CL大于3个标准差为检验不合格。EWMA控制图如图14、图15所示。

图14 整流平均值EWMA控制图

图15 均方根EWMA控制图

图14与图12中的曲线形状趋势类似,在图14中检测到多个点不合格,A10振动信号TSA数据的整流平均值稳定性较差,可能导致故障检测结果准确率严重下降。图15与图13中的曲线类似,在图15中检测出5个点不合格,且曲线末端的2个点确定发生故障。初步推断均方根数据用于故障检测效果更好,更准确的异常检测结果还需进一步观察三维数据融合后的多变量控制图。

MEWMA控制图是EWMA控制图的多变量形式,可以在EWMA控制图中同时监视2个或多个相关过程特征。MEWMA控制图中,每个标绘点都包括来自所有以前数据的加权,可以检测较小的过程偏移。

MEWMA统计量的计算方法如下:

(7)

式中:Zi表示先前子组的权重:

Zi=rXi+(1-r)Zi

(8)

加权的协方差矩阵为:

(9)

式中:

(10)

选取A10、A11和A12的VAR模型与TSA残差的整流平均值与均方根,绘制三维变量的MEWMA控制图,如图16、图17所示。

图16 整流平均值MEWMA控制图

图17 均方根MEWMA控制图

图16中在曲线的首尾两端检测出多个异常点,且超过控制限距离较大,与实际不符。图17中虽然同样检测出多个异常点,但大部分异常点都接近控制限,这是由于MEWMA控制图常用于检测数据中较小的过程偏移,导致其敏感度较高,最终造成虚警率过高。

广义方差控制图可以同时监视2个或多个相关变量的过程变异性是否受控制,此控制图是极差控制图、标准差控制图和移动极差控制图的多变量形式。选取A10、A11和A12的VAR模型与TSA残差的整流平均值与均方根,绘制三维变量的广义方差控制图,如图18、图19所示。

图18 整流平均值广义方差控制图

图19 均方根广义方差控制图

图18中在曲线首端检测出2个异常点,然而齿轮箱在试验开始阶段正常运行,因此不符合实际。图19中在曲线前段检测出49、60 2个异常点,在曲线后段检测出102、103、114、125 4个异常点,最后在曲线末端检测出151、152 2个异常点。然而,齿轮箱在试验中一直运行直至故障停机,控制图中不应出现较多的异常点。因此,均方根广义方差控制图的虚警率略高。

使用T2控制图可以同时监视2个或多个相关变量的过程位置是否受控制。此控制图是均值控制图和单值控制图的多变量形式。

当数据位于子组中时,T2的计算方式如下:

(11)

式中:

(12)

(13)

T2控制图的控制上限为:

(14)

图20是选取了A10、A11和A12的VAR模型与TSA残差的整流平均值,绘制了三维变量的T2控制图。该控制图在齿轮箱运行早期检测出多个异常点,不符合实际。

图20 整流平均值T2控制图

图21是选取了A10、A11和A12的VAR模型与TSA残差的均方根数据,绘制了三维变量的T2控制图。

图21 均方根T2控制图

均方根T2控制图检验结果如表2所示。

表2 均方根异常结果

表2表明均方根T2控制图检测出了3个异常点,文件编号分别为60、151、152,结合图21进行分析,得出齿轮箱故障诊断结论:齿轮箱在文件编号为60时,出现微小故障,产生了早期微弱的故障信号。在文件编号为151、152时,齿轮箱振动幅度巨大,出现故障,失效停机。

4 结论

1) 以14号齿轮箱为研究对象,采用TSA算法去除振动信号中的背景噪声。对振动文件进行分工况处理后,选取部分健康数据建立VAR模型,得到整体数据的残差,并提取残差的整流平均值、均方根指标。

2) 以残差均方根和整流平均值指标建立统计过程控制曲线,比较MEWMA、广义方差以及T2控制曲线,最终确定均方根T2控制图检测效果最佳,判断齿轮箱异常位置为第60、151、152个文件。以多维数据进行分析,可以更为准确地判断出异常点,不仅能找出故障出现的异常点,也能捕捉到早期微弱的故障信号。

3) 相比于传统的针对单一工况下的齿轮箱故障诊断方法,本研究对齿轮箱振动数据进行了分工况处理,可以判断齿轮箱全运行状态下的健康状态,提升了应用性。