数据驱动的旋转设备性能退化趋势预测方法

王庆锋,刘家赫,刘晓金,许述剑

(1.北京化工大学 高端机械设备健康监控及自愈化北京市重点实验室,北京 100029；2.中国石油化工股份有限公司青岛安全工程研究院，山东 青岛 266000)

1 问题的提出

中国流程工业生产装置正向大型化、高速化方向发展，旋转设备事故造成的非计划停工生产损失巨大，有时还会造成灾难性的安全、环保及人身伤害事故发生，以早期故障预测、健康状态评估、故障模式自动辨识和剩余工作寿命或性能退化趋势预测为支撑技术的预测性维修受到越来越多企业家和研究者的关注。旋转设备性能退化分为突发性(如图1a)和渐发性(如图1b～图1d)两类，前者性能退化发生、发展过程无明显的征兆;后者性能退化发生、发展过程具有一定的规律性。旋转设备突发性故障(如图1a)难以预防，围绕固定阈值报警线反复穿越报警(如图1b),给操作人员带来“报警疲劳”，固定报警阈值线下设备性能退化状态信息难以追踪(如图1c)，LH报警线和HH联锁报警线时间间隔过短(如图1d),来不及维护引发联锁停机等问题，是工程上预测性维修迫切需要解决的难题。面向工程实际需求，本文旨在研究提出一种数据驱动的旋转设备性能退化趋势预测方法，为开展旋转设备预测性维修提供技术支撑[1]。

当前，国内学者早期故障预测研究方法主要包括机理模型法、数据驱动模型法、机理和数据驱动混合模型法。基于机理模型的早期故障识别或诊断成果有：最优小波尺度循环谱滚动轴承早期故障识别方法[2]；基于最小熵解卷积与稀疏分解的滚动轴承微弱故障识别方法[3]；基于维度排列熵与支持向量机的轴承早期故障诊断方法[4]；变分模态分解识别滚动轴承早期故障方法[5]；基于连续小波变换的维护特征提取与故障识别方法[6]；无量纲指标分析法[7]和多元角域指标离群检测法[8]等风电机组齿轮箱早期故障诊断方法；采用符号时间序列分析的轴承早期微弱故障预警方法[9]；基于Gabor重排对数时频脊流行早期故障预测方法[10]；基于支持向量机(Support Vector Machines, SVM)和双时频微弱故障特征增强的轴承早期故障识别方法[11]；基于混合特征提取的双维度增强型K近邻(Enhanced K-Nearest Neighbor, EKNN)的滚动轴承早期故障分类方法[12]等，这些早期故障预测、识别或诊断方法建立在旋转设备复杂故障机理分析、故障诊断、信号分析处理和外部专家先验知识基础之上。数据驱动模型法早期故障预测相关研究成果有：基于变分模态分解(Variational Mode Decomposition, VMD)和深度置信网络的早期故障预测方法[13]；基于深度学习和大数据分析的旋转设备早期故障检测方法[14]；基于规范变量分析监测差异性的早期故障检测方法[15]；动态自学习阈值算法和l1趋势滤波算法相结合的智能预警方法[16]；基于正交邻域保持嵌入(Orthogonal Neighborhood Preserving Embedding, ONPE)与多核相关向量机(multicore Relevance Vector Machines, RVM )的滚动轴承早期故障诊断方法[17]等，这些早期故障预测方法依靠信号处理和特征工程技术，无需知道旋转设备零部件的劣化机制和故障机理模型，但模型训练仍然需要大量的故障案例数据。目前，机理模型法和数据驱动模型法研究成果仍然未有工程应用成果，国内流程工业服役状态监测系统在早期故障检测、故障诊断等方面智能化程度低，还需要依赖故障机理、故障特征信号等先验知识和外部专家经验。

开展设备早期故障预测一般基于原始数据抽取的故障敏感特征参数或构建的健康指数模型而不是原始数据[18]。现代信号处理和特征工程技术为旋转设备故障特征信号提取提供了丰富的手段，包括时域特征提取法、频域特征提取法和时—频域特征提取法3种。短时傅里叶变换[19]、小波包变换[20]、经验模态分解[21]等时频域特征提取方法能够通过时间和频率的联合函数来表示信号，虽然该类方法对信号进行诊断分析有良好的效果，但是时频域信号处理技术需要对原始信号进行预处理且计算一般相当复杂。无需对原始信号进行预处理、无需知道旋转设备故障机理及损伤机制、无需各种故障案例样本数据，直接利用旋转设备运行正常状态原始信号数据和实时监测信号数据，通过简单时域、频域或时频域特征参数计算实现旋转设备性能退化早期检测是本文研究探索的方法。利用运行可靠性函数[22]，即时域互相关函数、频域凝聚函数和谱距离指标函数，分别提取互相关系数、相干系数和谱距离指标，用来表征旋转设备当前运行状态监测数据与参考数据的差异性，并采用归一化度量区间[0,1]的数字表征旋转设备性能退化状态。本文拟提出一种基于谱距离指标运行可靠性曲线l1趋势滤波[23]后曲线斜率变化动态追踪的旋转设备性能退化趋势预测方法，构建一种数据驱动的旋转设备性能退化趋势预测模型，来解决工程上设备性能退化趋势预测及跟踪的难题。

2 设备性能退化状态表征

何正嘉等[22]提出基于旋转设备正常状态振动监测数据和实时振动监测数据，利用相关系数、相干系数、谱距离指标等方法计算运行可靠性，建立旋转设备运行状态与可靠性关联映射模型，实现运行可靠性评估的方法。旋转设备从“运行正常”到“运转到坏”，理论上其可靠性会历经从“1”到“0”的变化过程，同样也预示旋转设备的性能退化演化过程。

(1)相关系数运行可靠性模型

对于同一个设备，以正常运行状态振动监测原始波形数据为参考，采用互相关函数表征其运行可靠度，互相关系数从“大”到“小”表示旋转设备运行状态从“正常”到“故障”的性能退化过程。

(1)

式中：x表示运行正常状态；y表示待评估状态；Cxy(τ)为运行正常状态信号x(t)和待评估状态信号y(t)之间的正协方差函数；δx和δy分别为信号x(t)和y(t)的标准差；R(Cxy(τ))表示互相关系数，其取值范围在(0,1)之间。

(2)相干系数运行可靠性模型

对于同一个设备，以正常运行状态振动监测原始波形数据为参考，采用相干函数表征其运行可靠度，相干系数从“大”到“小”表示旋转设备运行状态从“正常”到“故障”的性能退化过程。

(2)

(3)谱距离指标运行可靠性模型

对于同一个设备，以正常运行状态振动监测原始波形数据为参考，采用谱距离指标函数表征其运行可靠度，谱距离指标从“大”到“小”表示旋转设备运行状态从“正常”到“故障”的性能退化过程。

(3)

式中：α为灵敏度系数，由设备性能退化趋势决定；R(Jxy)表示谱距离指标，其取值范围在(0,1)之间；Jxy表示运行正常状态信号x(t)和待评估状态信号y(t)之间的J散度，

(4)

式中：Sx(k)和Sy(k)分别为信号x(t)和y(t)的自功率谱；N为功率谱线个数。

本文选择谱距离指标运行可靠性模型表征设备性能退化状态，提取谱距离指标作为设备性能退化表征参量，进而构建数据驱动的旋转设备性能退化预测模型。

3 旋转设备性能退化预测模型构建

3.1 l1趋势滤波

最小绝对收缩和选择算子(Least Absolute Shrinkage and Selection Operator, LASSO)回归由Tibshirani于1996年提出，这是一种线性模型的压缩估计方法，有时也叫作线性回归的l1正则化，KIM等[23]借鉴趋势滤波的思想，对内在趋势的二次差分项进行l1范数处理，提出了l1趋势滤波，类似于LASSO回归原理，l1范数使得很多二次差分项缩减到0，由此得到的趋势项是分段线性的，即有折点的趋势。因此，l1趋势滤波非常适用于分段线性时间序列的趋势分析，用以估计曲线趋势、拐点或斜率的变化。

假设一个标准时间序列yi(i=1,…,n)由一个基本趋势xi和一个随机变量zi组成。趋势滤波就是从标准的时间序列中估计出基本趋势xi，或者估计出随机变量zi=yi-xi，l1趋势滤波是Kim等提出的一种H-P趋势滤波的变异，这种趋势估计是通过最小加权目标函数实现的，加权目标函数式为：

(5)

矩阵形式为：

(6)

(7)

λ是一个非负参数，用来控制趋势线的平滑性和平衡余项的大小，控制估计的趋势和信号冗余之间的平衡[24]。加权目标函数对x来说是一个严格的凸函数，因此只存在一个最小值xl1，xl1就是滤波后的趋势。使用l1趋势滤波对监测参数信号进行滤波处理，能够很好地得到机组真实的振动趋势；同时，可消除采集误差对设备早期故障预测准确性的影响。

3.2 数据驱动的旋转设备性能退化预测模型构建

图2为构建的数据驱动的旋转设备性能退化预测模型。

(1)筛选参考数据集

获取振动监测原始波形数据集X(tl)，l为收集到原始波形数据集的数组的数量。

1)对原始波形数据进行包络谱分析，分析是否有故障特征频率，若有故障特征频率，该原始波形数据不保存；若无故障特征频率，则该原始波形数据被保存到新构建的正常原始波形数据集X0(ts)中。s为正常原始波形数据集X0(ts)的长度，l≥s。

2)将长度为s的正常原始波形数据集X0(ts)的谱距离指标按照时间序列分解成长度为m的若干数组[xj,xj+1,…,xj+m-1](j=1,2,…,s-m+1)。

3)令j=1，将长度为m的数组[x1,x2,…,xm]放到一个固定窗w1内，求取窗口内谱距离指标数据l1趋势滤波后曲线上导数的平均值km。

4)j分别取2,3,…,s-m+1，采用类似办法组建w2,w3,…,ws固定窗，分别求取固定窗口内谱距离指标数据l1趋势滤波后曲线上导数的平均值km+1,km+2,…,ks。

5)ki(i=m,m+1,…,s)代表窗口内数据的波动程度，找出ki中绝对值最小的值并记为k0，则k0所对应的原始数组即为筛选的正常运行状态振动监测原始波形数据，记为参考数据集XR(tm)。

(2)计算实时监测数据的谱距离指标运行可靠性

基于模型参考数据和实时监测原始数据，应用谱距离指标运行可靠性模型计算实时可靠性。

(3)谱距离指标运行可靠性曲线趋势滤波并计算滤波后曲线上点的导数值

谱距离指标运行可靠性曲线趋势滤波后曲线上点的导数值的计算过程表述如下。

1)构建一个长度为m+n的滑移窗。其中：m代表固定数组的数据个数，由筛选出的正常运行状态振动监测原始波形数据的谱距离指标数据组成；n代表活动数组的数据个数，由参考数据和实时监测原始数据计算得到的实时运行数据的谱距离指标数据组成。

2)对长度为m+n的滑移窗内谱距离指标运行可靠性曲线进行l1趋势滤波。

3)求取长度为m+n的滑移窗内谱距离指标运行可靠性曲线l1趋势滤波后曲线上最后一点的导数值，滑移窗口活动数组最后一点的导数值记录为ki。

以N个谱距离指标运行可靠性数据为例，取滑移窗口长度l=m+n=50(m=45,n=5)，n=5为N个实时谱距离指标运行可靠性数据的前5个。滑移窗口自左(第1个实时运行可靠性数据)向右(第N个实时运行可靠性数据)滑移，逐点求取运行可靠性曲线l1趋势滤波后曲线上点的导数值，其过程如图3所示。

如图3a所示，固定数组内有45个谱距离指标运行可靠性数据；活动数组内有5个谱距离指标运行可靠性数据。对固定数组和活动数组共50个谱距离指标运行可靠性数据进行l1趋势滤波，并求取滤波后曲线上最后一个点的可靠性数据的导数值(对应于第5个实时谱距离指标可靠性数据点)，该点的导数值记为k5。

如图3b所示，滑移窗口向右滑移一个数据，原来滑移窗内活动数组最左边的第1个数据去掉，剩下的第2～5个数据和第6个数据组成新的活动数组；滑移窗口内新的50个运行可靠性数据进行l1趋势滤波，并求取滤波后曲线上最后一个点的导数值(对应于第6个实时谱距离指标可靠性数据点)，该点的导数值记为k6。

如图3c～图3e所示，以此类推，滑移窗口滑移到第N个实时谱距离指标可靠性数据点处，滑移窗口内新的50个谱距离指标运行可靠性数据进行l1趋势滤波，并求取滤波后曲线上最后一个点的导数值(对应于第N个实时谱距离指标可靠性数据点)，该点的导数值记为kN。

(4)旋转设备性能退化判别准则

|ki～i+n-1|>M·CF·|kr|。

(8)

本研究将辛辛那提大学智能维护系统中心(Intelligent Maintenance Systems, IMS)实验平台上收集的滚动轴承数据集作为模型训练数据。在工程上，实时监测振动信号不可避免地出现尖峰噪声，同样实时监测振动原始数据的变化会导致谱距离指标波形曲线异常波动，当滑移窗滑过活动数据异常点时l1趋势滤波后的曲线就会恢复正常。因此，本文设置一条报警规则，当在某一时间点及其之后的时间点连续报警次数超过活动数组数据长度的1～2倍时，即|ki|～|ki+n-1|均满足旋转设备性能退化判别准则时，将报警开始时间点视为有效的早期故障报警点。

(5)旋转设备性能退化趋势加速劣化预警判别方法

若满足式(9)，则系统给出设备性能退化加速劣化的预警，它适合已经检测到设备性能退化早期故障点，工程上需要重点跟踪设备性能退化趋势加速劣化的场合，尤其是能解决设备长期工作在如图1c所示的报警线以下，采用固定阈值报警方法无法追踪旋转设备性能退化趋势加速劣化发展的情况。

|ki+1|>A·|ki|。

(9)

式中：i∈(5,6,…,N-1,N);A为管理系数，A一般取大于1的正整数，由管理者根据旋转设备故障的风险程度主观制定。

(6)旋转设备性能退化趋势预测模型参数选择

滑移窗内固定数组数据长度m过长，将导致检测早期故障点或性能退化趋势加速劣化点不敏感；滑移窗内活动数组数据n过长，将造成较多的性能退化早期误报警或性能退化趋势加速劣化误报警。本文基于大量实验数据验证、优化得出结论：滑移窗口中固定数组长度m与活动数组n长度之比维持在9:1最适宜。滑移窗口长度视模型训练样本数据总量及信号平稳性而定，本研究中窗口长度取50个数据。

采用辛辛那提大学IMS中心轴承实验数据集进行模型训练，分析模型中平滑系数λ与调整参数M的取值范围，研究表明：当λ在10～100之间取值时，M取值范围为[36,70]；M取值一定时，λ的取值对滑移窗口活动数据最后一点的导数值的计算结果影响不大。详细说明见4.1节。

4 实验数据模型训练和验证

4.1 旋转设备性能退化预测模型轴承实验数据训练

选用公开的辛辛那提大学智能维护系统中心实验平台上收集的滚动轴承数据集作为模型训练数据[25]。如图4所示，实验平台主轴上安装了4个ZA-2115双列轴承，并对其进行了强制润滑，同时保持2 000转/分钟的恒定速度。弹簧机构用于向轴承2和轴承3施加26 689 N的径向载荷，同时沿x轴和y轴在每个轴承上安装了两个PCB 353B33高灵敏度加速度传感器。使用采样率20 kHz的NI DAQ卡6062E每10分钟收集一次长度为20 480点的振动信号。

进行模型训练，共收集了实验数据984组，每组数据有20 480个点，轴承运行前期为运行正常阶段，然后轴承状态逐渐劣化，最终由于外圈故障导致停机。图5为第532组数据的包络谱图，未出现故障特征频率。采用类似的方法对第1～200组数据进行包络谱分析，也未发现轴承故障特征频率，由此可以判定1～200组数据为运行正常数据。首先，使用包络谱分析筛选后的数据构建正常原始波形数据集X0(ts)，其中s≤200；然后，再筛选构建最优的正常运行状态数据集XR(tm)，其中m=45；如图6和图7所示，采用包络谱逐点分析，一直到第533组数据才出现滚动轴承外圈故障特征频率f0，这与WANG等[26]基于Gabor小波系数平方包络频谱分析得到在第533组数据时刻轴承发生早期故障的结论一致。

对筛选的正常运行状态数据集XR(tm)分别计算出谱距离指标数据。取45个谱距离指标数据组成滑移窗口固定数组，并计算出固定数组可靠性数据l1趋势滤波后曲线上点的导数的均值k0=1.39×10-6；取第1～5个实时监测数据，计算谱距离指标数据组成滑移窗口活动数组；构建如图2所示的数据驱动的旋转设备性能退化预测模型；将533点作为轴承性能退化早期检测点，滑移窗移动到第533点窗口内最后一点导数的绝对值k533与k0的比值设定为调整参数M的最大值；若将530点作为轴承性能退化早期检测点，530点导数的绝对值k532与k0的比值设定为调整参数M的最小值。

如图8所示，当λ∈[10,100]取值范围、M为最大值时，研究得出M的取值范围为[69.2,70.4]；当M为最小值时，研究得出M的取值范围为[35.8,36.2]。基于IMS中心轴承实验数据构建并训练模型，该模型平滑系数λ=50、调整参数M=53。

4.2 旋转设备性能退化预测模型轴承实验数据验证

取平滑系数λ=50、调整参数M=53、修正系数CF=1，应用构建的数据驱动的旋转设备性能退化预测模型，得出所有可能的报警点。由表1可以看出：自第15组数据时刻连续2次报警、自第163组数据时刻连续5次报警、在第259组数据时刻1次报警、自第291组数据时刻连续5次报警、自第320组数据时刻连续5次报警、在第345组数据时刻1次报警、自第426组数据时刻连续5次报警、自第532组数据时刻连续报警，自第15组数据时刻、第21组数据时刻、第163组数据时刻、第259组数据时刻、第291组数据时刻、第320组数据时刻、第345组数据时刻、第426组数据时刻的连续报警次数都没有超过模型活动数组长度(n=5)。按照图2构建的数据驱动的旋转设备性能退化预测模型早期检测点的判别规则，第541组数据时刻为模型寻找的设备性能退化趋势早期故障检测点。如图9所示，采用图2构建的模型计算出的设备性能退化预警点为第532组数据时刻，比常规的基于有效值(RMS)判断的第703组数据时刻性能退化预警点早了(703-532)×10=1 710 min。本模型具有较强的工程实践意义。

表1 IMS中心轴承实验台NO.2轴承实验数据计算出报警点对应的滑移窗口活动数组数据序号

5 工程案例数据模型验证

以离心式压缩机转子不平衡故障案例数据对构建的数据驱动的旋转设备性能退化预测模型进行工程验证。中国某石化公司BCL409/B型循环氢离心式压缩机额定转速12 100 转/min，轴振动固定报警限值为63.5 μm(高限报警)/88.9 μm(联锁报警)。在第1 218组数据时刻，振动峰—峰值由37.17 μm上涨到71 μm；在第1 706组数据时刻，振动峰—峰值上涨到92 μm后联锁停机，检修发现转子、流道、过滤网结垢严重，更换新转子后开机运行正常。

5.1 离心式压缩机性能退化早期检测

对比图10和图11，采用本文构建的旋转设备性能退化早期预测模型可以判断出第787组数据为性能退化早期故障检测预警点(如表2)，而采用传统的峰—峰值固定阈值报警方法在第1 218组数据时刻才发生报警，早期故障预警方法前者比后者提前了(1218-787)×10=4 310 min，这对于避免非计划停机、制定预测性维修时间具有非常重要的工程价值。

5.2 离心式压缩机性能退化早期检测点验证

如图12所示，在转子不平衡故障发生时，振动信号基频上的振幅会发生变化，压缩机运行正常状态下基频上的幅值稳定在4 μm。按照本文构建的旋转设备性能退化预测模型，以基频上的幅值数据代替谱距离指数运行可靠性数据，计算出在第790组数据时刻基频上幅值开始增加(图12)，比谱距离指数法延迟了(790-787)×10=30 min的时间间隔(谱距离指数表征设备性能退化状态比有效值、峰—峰值灵敏)。采用l1趋势滤波和动态自学习阈值控制线方法[16]，近似确定基频上相位在第825组数据时刻发生了变化(如图13)。依据基频上幅值增大、相位也发生变化可以判定为转子不平衡故障发生的先验知识，证明在第787～825组数据之间某时刻有转子不平衡故障发生。保守测算，采用第825组数据时刻作为离心压缩机性能退化早期故障检测预警点，相比采用传统的峰—峰值固定阈值报警方法在第1 218组数据时刻才发生报警，早期故障预警方法前者比后者提前了(1 218-825)×10=3 930 min，这对于避免非计划停机、实施预测性维修具有非常重要的工程价值。

6 结束语

针对在役状态监测系统采用常规固定阈值报警难以追踪到旋转设备性能退化的发生及发展趋势的工程实际问题，构建了一种数据驱动的旋转设备性能退化预测模型，提出了一种基于运行可靠性曲线l1趋势滤波斜率变化动态监测的旋转设备性能退化趋势预测方法。构建的模型只需要旋转设备运行正常状态监测原始数据，它不需要依赖故障机理、故障特征信号等先验知识以及外部专家经验。

利用辛辛那提大学智能维修系统(IMS)中心轴承实验台No.2轴承运转到坏的加速度数据对构建的旋转设备性能退化预测模型进行了训练，得出模型中l1趋势滤波平滑系数λ∈[10,100]、M∈[36,70]时，对于滚动轴承性能退化早期故障检测点结果影响不大。取平滑系数λ=50、调整参数M=53、修正系数CF=1，应用构建的数据驱动的旋转设备性能退化预测模型可以确定滚动轴承性能退化早期检测点为第532组数据时刻，这种方法比基于有效值(RMS)判断的第703组数据时刻早了1 710 min。

应用中国某石化公司离心压缩机转子不平衡故障案例的振动位移峰—峰值数据，对构建的旋转设备性能退化预测模型进行了验证。应用构建的旋转设备性能退化预测模型，取平滑系数λ=50、调整参数M=53，计算修正系数CF=4.61，确定该离心式压缩机性能退化发生在第787组数据时刻；取平滑系数λ=50、调整参数M=53，计算修正系数CF=5.14，转子基频上幅值增加发生在第790组数据时刻。采用l1趋势滤波和动态自学习阈值控制线方法，近似确定基频上相位在第825组数据时刻发生了变化。基于先验知识判断：离心压缩机在第787～825组数据之间某时刻发生了转子不平衡早期故障，比传统的峰—峰值固定阈值报警方法至少提前3 930 min预测到设备性能退化早期故障的发生。

应用本文构建的数据驱动的旋转设备性能退化预测模型，只需要运行正常状态下的振动原始数据，无需依赖外部专家先验知识，能够准确预测和追踪旋转设备性能退化趋势的发生和发展，这为工业4.0环境下实现基于工业大数据分析的预测性维修提供了技术支撑。

未来将把本文提出的旋转设备性能退化趋势预测方法，应用于炼油厂催化裂化装置烟气轮机、主风机叶片断裂突发性故障的早期预警，并根据模型测试和验证情况不断优化模型参数的设置。