基于层次聚类的状态监测数据衰退模式挖掘*

刘柏兵， 宋 东， 李春晓

(西北工业大学航空工程系 西安，710000)

引 言

随着现代复杂系统的功能要求、性能要求和使用环境要求的不断提升，对系统的故障预测和健康管理的要求也在不断提升。在PHM技术中，一项基本而重要的任务就是了解系统的健康状况及其发生性能衰退/健康状态降级的相关规律，积累系统的性能衰退/健康降级模式，进而建立其全面的系统衰退模型[1]。建立衰退模型的一种方法，就是对状态监测数据进行时间序列数据挖掘[2]。

国外较早开展该研究的机构是IBM公司的Pazzani和Agrawal 研究小组。在时间序列特征表示方面，目前已有不少方法，其中分段线性表示就是一种简单直观的特征表示方法，它使用线性模型对时间序列进行分割， 根据不同的分割方法可以使用不同的分割策略来实现[3]， 如滑动窗口、自底向上和自顶向下。滑动窗口和自底向上方法的时间复杂度为序列长度的平方阶， 而自顶向下的时间复杂度为线性阶。其中针对滑动窗口方法对时间序列拟合效果不佳这一问题，Shatkay作了对应的解释和分析。Park等[4]提出通过采用单调变换来改进该方法，虽然能有效地分割平滑的时间序列，但对时间序列含有噪声则无能为力。李爱国等[5]提出了一种对时间序列进行在线分割的递推算法，能够实时发现和预测时态模式，包括了分段聚合近似方法、符号化表示方法、基于域变换的表示方法、奇异值分解和基于模型的表示方法等。

针对基于状态监测数据的衰退模式挖掘中存在的噪声、序列时域形变或漂移不识别等问题，笔者提出了一种P-D-H聚类方法,其中：P是指通过分段聚合近似(PAA)方法对由状态监测数据形成的退化轨迹时间序列进行模式表示，这样能够有效地对时间序列进行降维，在降低了时间消耗的基础上充分反映了时间序列的信息，提高算法的效率；D是指采用动态时间弯曲距离(DTW)代替欧式距离作为模式序列的相似性度量，克服了欧式距离当序列在时间轴上发生轻微偏移时变化很大的缺点，提高了算法的准确性；H是指采用层次聚类方法实现衰退模式聚类，可以表征聚类的全过程，并且详尽展示了聚类簇间的距离关系。该方法充分利用现有数据资源，对系统衰退模式进行挖掘，为系统的健康预测提供依据，具有一定的理论价值和应用价值。

1 基于P-D-H的系统衰退模式挖掘方法

1.1 退化轨迹的建立

产品在存储或使用中发生的性能衰退表现为产品的某些性能特征参数不断退化并超过其规定范围，这种用于描述产品衰退过程的性能特征参数变化曲线/序列，称为观测产品的“退化轨迹”[6]。基于状态数据的退化轨迹P-D-H聚类方法的衰退模式挖掘，首先就是要从状态监测数据中提取系统的健康特征参量，特征参量形成的时间序列趋势成分，可以作为退化轨迹表征系统的衰退过程[7]。处理过程如图1所示(其中MA表示采用移动平均的方法)：a.对每组信号进行特征提取，得到不同特征的变化序列；b.提取特征序列，采用平滑法进行趋势分析；c.确定健康特征，获得退化轨迹。

图1 退化轨迹建立流程图Fig.1 Deletion trajectory build process

此处给出本研究所采用特征时间序列趋势成分提取方法以及斯皮尔曼等级相关系数衰退特征选取方法。

1.1.1移动平滑法

假设时间序列

F(t)={Fi(t)|i=1,2,…,m;t=1,2,…,n}

是状态监测参数在t个时间点下，m个特征参数的时间序列集，它可以表示为Fi(t)=Xi(t)+μi(t)，其中：Xi(t)为特征序列Fi(t)采用移动平均平滑法进行处理得到的趋势成分序列；μi(t)为去除趋势成分后的Fi(t)残差序列。则平滑后的趋势序列Xi(t)满足

(1)

1.1.2斯皮尔曼等级相关系数衰退特征选取方法

斯皮尔曼等级相关系数适合于定序变量或不满足正态分布假设的等间隔数据，它利用两变量的秩次大小作线性相关分析，对原始变量的分布不作要求，属于非参数统计方法，适用范围比较广。 采用斯皮尔曼等级相关系数ρ作为衰退特征选取指标的思路是，计算状态监测参数X(t)={Xi(t),i=1,2,…,m}中每个Xi(t)与时间序列t=1,2,…,n之间的斯皮尔曼等级相关系数ρi，选取ρi值最大的若干Xi(t)作为系统的退化轨迹，记为Di(t)。实际应用中，变量接的连结是无关紧要的，于是可以通过简单的步骤计算ρ。对于样本容量为n的样本，n个原始数据Xi和Yi被转换成等级数据xi和yi[8],得到被观测的两个等级差值di=xi-yi，则ρ为

(2)

1.2 退化轨迹集聚类

通过上述退化轨迹建立方法对状态监测数据建立退化轨迹后，得到样本个数为N的监测数据退化轨迹时间序列集D(t)={Di(t),i=1,2,…,N}，采用PAA的方法对退化轨迹D(t)进行近似表征，然后通过DTW度量近似序列之间的距离，以此作为退化轨迹序列之间的相似度度量准则，最后采用自底向上的层次聚类算法对退化轨迹集合中的样本进行聚类。

1.2.1退化轨迹的分段聚合近似

PAA是一种时间序列的基于划分的近似表示，其意义简单明确，较好地保持了时间序列在时域上的变化趋势特征，降低了时间序列中包含的噪声信息。由于本研究中的衰退模式挖掘目的在于发现时间序列变化趋势的相似性，采用分段线性表示的方法对时间序列降维和特征表示是一种可行的方法[9]。将PAA近似表征方法与DTW距离度量相结合，可以克服范数距离度量对序列时域形变或漂移不识别的问题，其方法原理如下。

(3)

为了使数据的维数由n维向量降为N维向量,首先需要将数据分割成等长的N小段,之后需要计算出每一段的均值,再用这些数值重新组成向量作为时间序列降维后的表示[10]。

1.2.2DTW距离

在进行时间序列样本相似性度量的时候，最基本距离函数主要为普通范数距离和动态时间弯曲距离(DTW)。普通范数距离Lp在计算时要求两条时间序列等长，且计算方法严格按照点与点的顺序一一对应，对不等长的时间序列无法计算其普通范数距离。与普通的范数距离不同，DTW在聚类时不要求时间序列之间的点与点匹配，当两条序列长度不同时，允许时间序列点自我复制后再进行对齐匹配。DTW对时间序列的突变或异常点不敏感，更加适用于此类时间序列的距离度量，且可以实现异步相似性比较，普通范数距离只能实现同步比较。因此，笔者采用DTW距离作为相似性度量的准则,其原理如下。

对于两条退化轨迹序列Q=q1,q2,…,qn和C=c1,c2,…,cm，定义一个n行m列矩阵

(4)

将两个时间序列分别置于二维坐标的两轴，如图2所示，进而可以定义弯曲路径：在两个不同时间序列间的距离矩阵中，定义时间序列间相异性关系的K个连续的矩阵元素的集合，称为弯曲路径。

W=w1,w2,…,wi,…,wk

(5)

图2 弯曲路径示意图Fig.2 Schematic diagram of curved path

将退化轨迹序列Q，C之间的动态时间弯曲距离DTW定义为两序列距离矩阵中弯曲路径最小的路径长度，称为最优弯曲路径，即

(6)

其中：分母K是为了在比较不同长度的路径时有统一的标准。

DTW距离的求解采用累积距离矩阵的动态规划方法来计算。累积距离矩阵实际上是一种递推关系，对于上文中定义的两个时间序列Q和C，累积距离定义为

(i=2,3,…,n;j=2,3,…,m)

(7)

DTW距离为两条时间序列对齐匹配关系的最佳值给出了定义，同时为时间序列的时间弯曲和不同长度的时间序列之间的相似性度量提供了支持[11]。

1.2.3P-D-H层次聚类算法流程

本研究采用的层次聚类不需要预先设定聚类数目，允许从任意维度对数据集进行划分，它的实现建立于对数据集中的两两样本进行相似性(距离)计算的基础上，通过不断合并相似性最高(距离最小)的样本来构成新的集群，将合并得到的集群看成一个新的样本，再次进行相似性计算和样本的合并，如此反复进行运算，直至所有样本合并为一类。层次聚类方法的优点在于可以表征聚类的全过程，并且详尽地展示了聚类簇间的距离关系[12]。其具体流程如下：

1) 对输入样本集D={x1,x2,…,xm},定义样本间距离函数dist(xi,xj)和类间距离函数d(Ci,Cj);

2) 将样本集中的每一个样本化为一簇，即Cj={xj}；

3) 对于聚类集合Cj={xj}进行两两的距离计算，构建距离矩阵M，M(i,j)=M(j,i)=d(Ci,Cj)；

4) 找出距离最近的两个聚类簇Ci*和Cj*，形成新的聚类簇及编号；

5) 删除距离矩阵M的相应行列，计算新的聚类簇与其他簇的距离，更新距离矩阵；

6) 循环进行第3步～第5步，直到全部样本聚为一类；

7) 确定划分的聚类数目k。

这里样本间距离为dist(xi,xj)，采用DTW距离计算方法，记为distDTW；类间距离d(Ci,Cj)采用关于集合的某种距离即可，可选用的类间距离有最小距离、最大距离和平均距离。

由上述可知，最小距离由两个类别间的最近样本决定，最大距离由两个类别间的最远样本决定，而平均距离则由两个类别中的所有样本共同决定。最小距离和最大距离的计算比较简单，但对簇的整体特性刻画比较片面，平均距离能更好地刻画簇的特征，但计算量稍大一些。

1.2.4数据标准化

1) Z-score标准化

(8)

其中：μ和σ分别为X的均值和标准差。

2) Min-max标准化

(9)

3) 规范化

(10)

2 试验分析

通过对轴承磨损衰退模式挖掘来验证基于P-D-H的系统衰退模式挖掘方法的有效性，本节的数据来源于IEEE可靠性学会与FEMTO-ST Institute共同组织的2012IEEE PHM大会，目的是用于轴承剩余寿命预测比赛。数据采集于FEMTO-ST Institute的PRONOSTIA试验台[13]进行的17个轴承加速退化试验[14]。

试验在数小时内完成了17个轴承的加速退化，失效判据统一为监测量的振动加速度，超过20g时停止试验。工况与试验件的信息如下。

工况1：1800r/min，4000N；试验轴承编号1-1～1-7，共7个；

工况2：1650r/min，4200N；试验轴承编号2-1～2-7，共7个；

工况3：1500r/min，5000N；试验轴承编号3-1～3-3，共3个。

采样方式和参数单位如表1所示，即振动加速度信号每间隔10s进行1次采集，每次记录2560个数据点；温度信号每秒记录1个数据点。

表1 采集参数详情表

2.1 退化轨迹的建立

2.1.1 特征提取

基于对振动加速度信号的观测和一般了解，采取如下的特征提取方法：首先，对每组信号直接提取2个有量纲统计特征和2个为无量纲统计特征，包括均方根(root mean square,简称RMS)、峰值(peak)、峰值因子(crest factor，简称CF)、谱峭度(spectrum kurtosis，简称SK)；其次，对振动信号进行带宽为500 Hz的数字带通滤波，得到频段0～12 kHz的24个滤波信号，对每个滤波信号同样提取4个统计量。峰值特征在计算的时候，将每次观测数据中的3个最大值进行平均来得到峰值特征，以消除异常值的影响[15]。各个统计特征量计算公式如表2所示。

表2 统计特征量计算方法

表2中：xi为观测值；N为观测点数，此处为2 560；i∈[1,N]，为样本均值；Xord(ti)为降序排列的观测序列。

2.1.2 平滑趋势

对特征序列进行平滑处理，以得到其平滑趋势。以轴承1-1为例，观测记录434组数据，对未滤波的单次2 560个采样值提取以上4个统计量，得到的特征序列如图3所示。此处选取窗宽n=20，平滑后序列计算方法参见式(1)，平滑后的趋势序列如图4所示。

图3 轴承1-1特征时间序列Fig.3 Bearing 1-1 characteristic time series

图4 轴承1-1特征平滑时间序列Fig.4 Bearing 1-1 features smooth time series

2.1.3 退化轨迹的建立

通过式(2)及相关定义，对平滑趋势后的特征序列与其对应的记录时间进行斯皮尔曼等级相关性系数计算。综合所有数字特征和相关系数计算结果，发现均方根和峰值因子的与时间序列相关度系数均在0.5以下，因此不作为特征参数参与建立退化轨迹。在所有轴承信号分析中，峰值信号与时间的相关度系数值大多保持在0.8以上，表明了其良好的单调性，在轴承1～3中高达0.893 2，且各个轴承峰值时间序列数字范围均在0～20之内，具备较好的统一性。谱峭度与时间的相关性低于峰值信号，保持在0.6以上，数值范围起于3止于4。因此选取未滤波信号峰值(以下简称峰值P)和5.5～6 kHz频段下滤波信号的谱峭度(以下简称谱峭度K)作为两个不同的健康特征量。这也符合工程实际中的一般认识：轴承在健康状态下的振动信号是服从高斯分布的，其谱峭度值为3；随着试验的进行，轴承的磨损逐渐加重，其振动信号的谱峭度值越来越偏离3，表现出渐变的退化趋势。图5和图6分别为轴承1-1水平振动加速度信号峰值特征序列和谱峭度特征序列。

图5 峰值特征平滑序列Fig.5 Peak feature smoothing sequence

图6 谱峭度特征平滑序列Fig.6 Spectrum kurtosis feature smoothing sequence

2.2 聚类过程

选取了两种特征量作为健康表征参数，故对两种特征参数下的样本退化轨迹分别进行聚类，聚类的总样本个数为34，样本时间序列长度不完全相等。分段聚合数目N=100聚类过程按照2.2节所述的聚类过程，这里直接给出得到两种特征下的样本聚类谱系分别如图7和图8所示。

图7 基于谱峭度的聚类谱系图Fig.7 Based on the clustering spectrum of spectral kurtosis

图8 基于峰值的聚类谱系图Fig.8 Based on the clustering spectrum of peak

由图7、图8可观察到，在谱峭度和峰值特征表征下，大多数的样本退化轨迹比较相似，基于谱峭度特征的样本在更小的距离上聚合，离群样本与其他样本差异性较大；而基于峰值特征的聚类对样本区分度更加明确。综合以上两图，将最终聚类数目确定为4，得到聚类划分结果如表3和表4所示。

表3 基于谱峭度的模式聚类结果表

表4 基于峰值的模式聚类结果表

由于两种不同的特征参数与不同故障模式相关系数有所不同，且试验中可能存在噪声的干扰，因此，两种不同特征参数下的模式3和模式4聚类的结果可能存在差异。

图9选用样本1,6,20,23来分别表示4种衰退模式之间的差异，其中4条不同颜色的曲线分别代表了不同的样本，图中横轴为观测时间点，纵轴为SK值。同理，图10选用样本1,6,19,20来分别表示4种衰退模式之间的差异，图中横轴为观测时间点，纵轴为峰值。

图9 基于谱峭度的衰退模式样本退化轨迹Fig.9 Sample degradation trajectory of degradation pattern based on spectrum kurtosis

图10 基于峰值的不同模式样本退化轨迹Fig.10 Sample degradation trajectory of different patterns based on peaks

由聚类结果可以看到，样本6和样本20在两种特征下的退化轨迹都表现出与大多数样本不一致的特性。在两种健康特征下，聚类获得的4种退化模式下的退化曲线均有明显的形状差别。样本6基于谱峭度的退化轨迹与基于峰值的退化轨迹区别明显，得到的退化轨迹和聚类结果可用于建立系统的衰退模型，具有较好的效果。

3 结束语

笔者 针对PHM技术中系统衰退模型的建立进行了深入的研究，对比了多种现有的基于状态监测数据的退化轨迹建立和时间序列聚类的模式挖掘方法，总结了其中的缺点和不足。在此基础上提出了一种P-D-H聚类方法，利用此方法对IEEE可靠性学会与FEMTO-ST Institute共同组织的轴承剩余寿命预测比赛数据进行了衰退模式挖掘，获得了4种退化模式下的退化曲线，验证了方法的有效性，为状态监测数据的系统健康衰退模式挖掘提供了一种有效的方法。