孟 臻 胡 宇 唐恒光高级工程师

(1.北京信息科技大学 自动化学院，北京 100101；2.中国福利彩票发行管理中心，北京 100101)

0 引言

随着现代工业技术和信息技术的发展，新能源、通信、计算机和工业等各个领域的生产制造系统日趋复杂。由于结构的复杂性和各种潜在因素的影响，这些系统不可避免的存在故障隐患，隐患一旦诱发事故会造成不同程度的人员伤亡及经济损失。因此，系统的故障诊断方法成为研究者关注的焦点。目前常见的故障诊断方法主要包括以下3类：基于控制模型的故障诊断[1]、基于统计方法的故障模型诊断[2]、基于人工智能的故障诊断[3]。目前，大量研究着眼于优化故障诊断模型的准确率，通常由诊断有度衡量其改进结果。然而，实际应用中模型稳定性越高，成本代价付出越大。因此，有必要分析模型稳定性对故障诊断有效性的影响。目前，对模型准确率的衡量仍缺乏统一体系，主要方法有相对偏差、残差平方和法[4]等误差分析法，以及灰色相关理论[5]和基于统计数据的ED检验、置信区间[6]等方法。灰色相关理论可以实现多数据输入诊断，但只能应用于同一特征参数的范围内，不能分析不同特征范围的参数关联性，不适用不同范围特征参数故障诊断模型间的相对比较；残差平方和法是适用于回归模型的评估方法，受因变量和自变量绝对值大小的影响，不利于在不同故障诊断模型间进行相对比较；置信区间指标是基于测试数据的诊断推理结果，同时符合正态分布的假设，但是在数据规模较小的情况下将产生较大误差，且置信区间上限不随准确度的增加而收敛于1，所以不适用于故障诊断准确度高的模型。在上述研究的基础上，考虑故障诊断输出结果的真实分布情况，提出基于椭圆傅利叶描述子的故障诊断模型(Clustering Stability Index Method Based on Elliptic Fourier Descriptor，EFDCSI)稳定性评估方法，能够在不了解故障诊断模型的参数等情况下提供客观的评价，并采用基于特征聚类故障诊断模型的电动机轴承诊断应用实例进行对比验证，证明所提方法有效。

1 算法描述

1.1 椭圆傅利叶描述子

形状是目标物重要的视觉特征之一。现有的形状表示方法可以分为2大类：基于区域特征的形状表示方法和基于轮廓特征的形状表示方法。基于轮廓的方法主要利用目标覆盖区域边界的像素信息来描述形状[7-8]。

傅利叶描述子是一种经典的基于轮廓的形状表示方法，最早由COSGRIFF(1960)提出。其主要思想是利用表示形状整体频率的一组数据来描述轮廓特征，对旋转、平移等操作具有不变性，是形状研究最广泛的描述子。在算法研究方面，研究人员在改进基于傅利叶算子的形状表示方面做了很多工作，以提高对形状的描述能力。如ZHANG等[9]提出一种增强的通用傅利叶描述子抽取图像内容描述关键部分的方法，用以解决大数形状描述不适合泛型形状描述的缺点；LI等[10]提出一个基于区域的仿射不变量环傅利叶描述子进行仿射不变特征提取的方法，可用在具有多个组件的对象轮廓特征提取；KASAUDHAN等[11]提出一个增强版本的网格距离傅利叶描述子计算图像相似度的方法，提高图像匹配率；BELKHAOUI等[12]将傅利叶描述子与分水岭算法(Watershed Algorithm，WS)相结合，提出一种利用逆合成孔径雷达图像进行目标自动识别的过程和方法，解决雷达图像的目标识别问题。

1.2 EFDCSI算法

1.2.1 边缘检测

Canny边缘检测算子相对于Roberts Cross算子、Prewitt算子、Sobel算子、Kirsch算子具有明显优势。因此运用Canny算子进行边缘检测，以识别其轮廓边界。

首先为了平滑图像，以减少边缘检测器上明显的噪声影响，对图像采用式(1)所示大小为的高斯滤波器核对图像进行高斯滤波。

(1)

式中：

σ—标准差；

K0—核矩阵的维数，是正整数；

d1、d2—像素点的位置。

若图像中一个3×3的窗口为W，要滤波的像素点为e，则经过高斯滤波后e，像素点的亮度值为：

(2)

其中，*为卷积体符号，sum表示矩阵中所有元素相加求和。Canny算法采用4个算子来检测图像中的水平、垂直和对角边缘，边缘检测的算子返回水平方向梯度Gx和垂直方向梯度Gy的一阶导数，由此便可以确定像素点的梯度强度G和方向梯度θ，如式(3)。

θ=arctan(Gy/Gx)

(3)

将当前像素的梯度强度与沿正负梯度方向上的2个像素进行比较，如果当前像素的梯度强度最大，则该像素点保留为边缘点，否则被抑制掉，称之为最大值抑制。在施加最大值抑制后，仍然存在由噪声和颜色变化引起的一些边缘像素，为解决这些杂散响应，通过设立选择高低阈值来实现。如果边缘像素的梯度值高于高阈值，则标记为强边缘像素；如果边缘像素的梯度值小于高阈值大于低阀值，则记为弱边缘像素点，但如果弱边缘像素的8个邻域像素只要有一个为强边缘像素，则可以作为边缘点保留下来；如果边缘像素的梯度值小于低阈值，则被抑制。

1.2.2 边缘特征的傅利叶描述子

确定目标边界的一个起始点(x0，y0)沿逆时针方向以一定速度移动后，就可以用边界点的坐标对序列来描述边界。边界曲线被用复数形式定义：

s(t)=x(t)+jy(t)，t=0，1，…，N-1

(4)

式中：

t—沿边界曲线移动的单位弧长；

x(t)—t时刻移动点位置的横坐标；

y(t)—t时刻移动点位置的纵坐标；

j—表示虚数单位；

N—正整数，表示坐标对的个数。

为描述图像轮廓的需要，选定的起始点需要沿边界曲线移动一周。因此，s(t)是一个周期2π的周期函数。为取得边界曲线的椭圆傅利叶描述子，首先进行傅利叶级数展开，对s(t)可将其用一维傅利叶级数展开为：

(5)

式中：

l—第l个分量；

w—角频率；

T—周期。

根据式(4)可以有椭圆系数表达式：

sk=sxk+jsyk

(6)

(7)

根据三角函数与指数函数的关系，有：

(8)

则：

(9)

m是轮廓曲线的采样点数量，通常为轮廓曲线像素点数量的1/2。xi和yi是x(t)和y(t)第个样本点值。根据式(8)和式(9)，sk可以看作复数对之和。即

sk=Ak-jBk，s-k=Ak+Bk

(10)

(11)

式(5)可以被表达为：

(12)

称式(13)为轮廓曲线的椭圆描述子。

(13)

1.2.3 基于傅利叶描述子的故障诊断模型稳定性评价

如果基于特征聚类故障诊断模型的稳定性好，则类中心的凝聚力较强，即向某一类中添加或移除同类数据时，聚类结果的形状变化极小，反之则形状变化较大。假设类内元素数据变化前后映射图像的边界描述子分别为EFD1和EFD2，将2个故障分类结果轮廓形状的相似性定义为故障诊断模型的稳定性评价标准，如式(14)所示：

(14)

式中，

cov()—2个描述子的协方差；

σ—描述子向量的标准差。

sim的值越接近1，说明故障分类模型的稳定性越好。

2 实验与分析

2.1 数据描述

为了电动机轴承故障新技术、新理论的实验与验证，一个由ROCKWELL开发的电机轴承状态评估系统获取一系列可以用来验证或改进电机性能评价技术的电机性能数据库[13]。平台包括一个1.5kW的电动机，一个扭矩传感器/译码器，一个功率测试计，还有电子控制器。待检测的轴承支撑着电动机的转轴，驱动端轴承为SKF6205，风扇端轴承为SKF6203，轴承用电火花加工单点损伤，损伤直径分为0.177 8、0.355 6、0.533 4、1.016、1.016mm。其中，轴承外圈的损伤点在时钟：3、6、12点钟3个不同位置进行设置。电动机风扇端和驱动端的轴承座上方各放置一个加速度传感器，用来采集故障轴承的振动加速度信号。共采集正常轴承、单点驱动端和风扇端故障的数据。对于驱动端轴承，分别使用12kHz和48kHz的频率进行数据采集。对于风扇端轴承，仅使用12kHz的频率进行数据采集。数据文件均为Matlab格式(.mat格式)。每一个文件包含风扇端和驱动端振动数据，也包含电机转速。对于所有文件，变量名字的含义如下：DE：驱动端加速度数据；FE：风扇端加速度数据；BA：基本加速度数据；Time：时间序列数据；RPM： Rpm During Testing。部分数据，见表1。

2.2 基于特征聚类的故障诊断

聚类分析是一种无监督学习，无需事先定义类别或给出训练样例等先验来说明数据关系，可以直接将一个数据集划分为若干差异性较大的类，而类内数据具有极高的相似性。这对于无法确定标准信息征兆群的场景是非常适用的，如故障诊断。由于真实环境中系统参数、环境干扰与噪声很难准确确认，造成很难建立起精确的故障诊断模型，基于数据驱动的方法避免了对过程进行数学建模，在不清楚诊断对象机理信息的情况下，可以通过历史数据进行学习建模完成故障诊断。常用聚类算法有K-means[14]、Birch[15]、Em[15]、Dbscan[16]、Clara[17]等。本文主要研究故障诊断模型稳定性评价方法EFDCSI。因此，直接采用Clara、K-means、Dbscan方法进行基于特征聚类的电机驱动端轴承故障诊断。表1中数据DE是具有时序特征的驱动端加速度数据，如图1。

表1 部分驱动端轴承故障数据

图1 DE时间序列图

在电机轴承故障诊断过程中，故障特征提取的效果决定最终诊断率的高低，振动数据的峰均比(Peak-to-AverageRatio，PAR)、峭度(Kurtosis)、偏度(Skewness)涵盖振动的分布特征、统计特征、线性特征，能有效反映振动事件主要特性。因此，将这3个特征作为故障诊断的依据。具体计算方法如式(15)-(17)：

(15)

式中：

E—数学期望；

p—峰值功率，W。

(16)

(17)

式中：

U(t)—瞬时振幅；

q(x)—概率密度；

由于振动的部分特征向量可能具备一定的相关性，会影响聚类故障诊断模型的稳定性。所以，去除相关特征向量是完成故障诊断的第一步。利用主成分分析法提取不相关的特征向量，计算结果显示PAR、Kurtosis可以表征表1中DE的特征向量，见表2。

表2 DE的特征向量表

为有效验证所提的EFDCSI评价方法的可靠性，采用如戴维森堡丁指数(Davies-Bouldin，DB)[18]、邓恩指数(Dunn Validity Index，DVI)[19]、轮廓系数(Silhouette Coefficient，SC)[20]等10余项聚类评价指标对DE的特征向量的最佳聚类个数进行预评估，将电机轴承故障分为2类。即将采用Clara、K-means、Dbscan分别将故障分为2类。将特征集F按照3：1的比例从中随机抽取测试集T，然后分别对测试集T和特征集F进行二分类。分类结果，如图2。

(a)特征集F的Clara聚类

分别将测试集T与特征集F相对应的第一类数据提取出来，如图3。

(a)特征集F的Clara聚类结果的I类

为利用图像边缘检测技术评价基于特征聚类的故障诊断模型的稳定性，需要将二维空间点根据一定规则映射为图像。采用狄利克雷图对空间点坐标进行图形化映射，以任一点作为三角形的顶点，将与顶点欧氏距离最近的2个点与顶点连接，经过若干次迭代得到Delaunay三角网。测试集T与特征集F的第I类数据的泰森图，如图4。

(a)特征集F的Clara聚类泰森图

2.3 故障诊断模型的性能评估

经过Canny算子边缘检测后，图4的边缘轮廓，如图5。

(a)特征集F的Clara聚类泰森图轮廓

按照式(13)可得图5所示的故障分类轮廓曲线的椭圆傅利叶描述子，并将其归一化，部分数据，见表3。

表3 故障分类边缘的椭圆傅利叶描述子

聚类结果的验证通常采用2种技术，一类是簇内距离，如簇内总和的平方(within.cluster.ss)计算的是每个聚内部的距离平方，值越小说明类内数据的特征越相似，聚类效果越好；一类是簇间距离，如平均轮廓宽度(avg.silwidth)计算的是平均轮廓系数，值越大说明不同类的数据特征差距越大，聚类算法的数据区越好。将聚类算法评估指标：within.cluster.ss、avg.silwidth与提出的指标sim进行比较，见表4。由表4可知，sim值表明3种聚类算法的故障诊断效果K-means算法最好，Clara算法次之，Dbscan算法最差，与评估指标avg.silwidth和within.cluster.ss的结论一致，证明提出的EFDCSI方法对于基于特征聚类的故障诊断模型稳定性的评价是有效的。

表4 Clara、K-means与 Dbscan算法性能评价

3 结论

EFDCSI针对基于特征聚类的故障诊断模型稳定性效果评价问题，从聚类问题的核心思想出发，结合图像轮廓检测技术，将故障诊断模型稳定性评估问题转化为图形结构相似度计算问题，通过提取故障分类结果轮廓曲线的特征向量计算轮廓形状相似度来评价故障诊断效果，是一种新型的基于特征聚类故障诊断模型稳定性的评价方法。该方法可以在数据样本及数据方法未知的情况下，仅根据故障分类结果的可视化图形轮廓特征向量来评价模型故障诊断效果的稳定性。与经典聚类评估算法K-means、Clara、Dbscan的评估结果一致，且计算简单。从实验方法与原理来计，该评价适用于基于特征聚类的故障诊断模型稳定性评价。