滚动轴承作为主要旋转元件，广泛存在于核电厂的各类泵机、电机、压缩机以及汽轮机等大型机械设备中，其运行状态直接影响设备性能，当复杂的运行环境导致轴承出现故障时，可能引起整个系统设备的失效。据统计，旋转设备的故障中30%以上源于轴承故障^[1]，因此轴承部件故障诊断的方法研究，对提高核电厂设备状态监测精确度，增强核电厂运行安全性具有重要意义。

滚动轴承振动信号一般具有非平稳特性，直接对原始信号进行故障诊断是比较困难的。文献[2]采用包络谱分析技术对CPR1000核电厂泵组轴承故障进行了分析，但是文中轴承振动烈度阈值不具有普适性。文献[3]基于能量解调算法，研究了一种应用于核电厂故障诊断的方法。国外也有学者将递归神经网络应用于核电厂旋转机械的状态监测^[4]，但是文献[5]指出神经网络结构的选择过分依赖于先验知识或经验。相比之下，将振动监测信号处理与机器学习算法相结合的旋转机械故障诊断方法，已成为近年来设备故障诊断的研究热点。其中，基于EMD、ITD、VMD的信号处理方法同样得到了广泛关注，众多学者也根据其局限性作出相应的改进^[6-9]。

本文基于EWT的信号分解方法，将轴承振动信号分解为各阶AM-FM分量，通过K-L散度值筛选出包含原始信号故障特征的主分量，并结合样本熵及LZ复杂度的特征提取技术，采用GG聚类算法对旋转机械轴承进行故障诊断。实验测试分别与基于FCM、GK聚类以及基于EMD的信号分解方法作比较，从而验证了该方法能对旋转设备轴承的不同故障进行精确有效识别的优势。

1 信号预处理及主分量提取方法 1.1 经验小波变换原理

基于EMD自适应分解思想，Gilles^[10]对小波理论做出改进，提出了EWT的多分量信号分解方法。EWT信号分解的效果关键在于对Fourier频谱的分割方式，一般采用基于频谱极大值的边界检测方法，即若将信号频谱划分为N+1个连续区间，则先找到信号频谱的前N-2个最大的极大值点，再将边界集定义为相邻2个极大值点的中点或者最小值点。

然而旋转设备的轴承振动信号频谱较为杂乱，往往前几个最大的极大值集中在一个频率峰群，导致根据频谱极大值的方法分解效果较差。因此本文应用一种基于给定初始边界的自适应边界检测方法，该方法首先根据信号频谱给出N个初始边界集，然后检测出每个边界各邻域内的最小值作为新的边界集。其中检测的邻域范围是初始相邻边界距离的一半，从而避免了新边界出现混叠现象。该方法不仅能根据信号频谱特性进行分解，而且克服了检测边界过于集中的问题。因此针对旋转设备的振动信号，EWT具体实现步骤如下：

1) 应用傅里叶变换得到信号的频谱F(ω)，根据频谱特性给出N个初始边界集{ω_n}，根据本文的自适应边界检测方法，搜索得出最终分割边界{ω_n}，以及N+1个连续分割区间Λ_i=[ω_i-1, ω_i]，其中ω₀及ω_N+1分别为频谱2个端点；

2) 根据所确定的分割边界，定义每个分割区间Λ_i上的带通滤波器，即经验小波。根据Meyer小波构造经验小波，计算经验小波函数${\widehat \psi _n}(\omega)$及经验尺度函数${\widehat \varphi _n}(\omega)$，其表达式为:

$ {\hat \psi _n}(\omega ) = \left\{ \begin{array}{l} 1,\quad (1 + \gamma ){\omega _n} \le |\omega | \le (1 - \gamma ){\omega _{n + 1}};\\ {\rm{cos}}\left( {\frac{\pi }{2}\beta \left( {\frac{1}{{2\gamma {\omega _n}}}(|\omega | - (1 - \gamma ){\omega _{n + 1}})} \right)} \right),\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (1 - \gamma ){\omega _{n + 1}} \le |\omega | \le (1 + \gamma ){\omega _{n + 1}};\\ {\rm{sin}}\left( {\frac{\pi }{2}\beta \left( {\frac{1}{{2\gamma {\omega _n}}}(|\omega | - (1 - \gamma ){\omega _n})} \right)} \right),\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (1 - \gamma ){\omega _n} \le |\omega | \le (1 + \gamma ){\omega _n};\\ 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{其他}} \end{array} \right. $

(1)

$ {\hat \varphi _n}(\omega ) = \left\{ {\begin{array}{*{20}{l}} {1,\quad |\omega | \le (1 - \gamma ){\omega _n};}\\ {{\rm{cos}}\left( {\frac{\pi }{2}\beta \left( {\frac{1}{{2\gamma {\omega _n}}}(|\omega | - (1 - \gamma ){\omega _n})} \right)} \right),}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (1 - \gamma ){\omega _n} \le |\omega | \le (1 + \gamma ){\omega _n};}\\ {0\quad {\rm{ 其他 }}} \end{array}} \right. $

(2)

式中β(x)和γ表达式为：

$ \left\{ {\begin{array}{*{20}{l}} {\beta (x) = {x^4}(35 - 84x + 70{x^2} - 20{x^3});}\\ {\gamma < \mathop {{\rm{min}}}\limits_n \left( {\frac{{{\omega _{n + 1}} - {\omega _n}}}{{{\omega _{n + 1}} + {\omega _n}}}} \right){\rm{ 且 }}0 < \gamma < 1} \end{array}} \right. $

(3)

3) 根据传统小波变换定义方法，将信号与经验小波函数内积产生细节系数，同时信号与经验尺度函数内积产生近似系数；

4) 再由细节系数、近似系数与经验小波函数、经验尺度函数的卷积即可得到各AM-FM分量。

自适应边界检测方法的EWT能够较好分解出多个分辨率的信号分量，但是需要指出，初始边界集可以根据信号各频率峰群给定，同时也可以类似文献[11]将相邻极大值间的极小值作为初始边界集。

1.2 基于K-L散度的主AM-FM分量提取

信号经过EWT分解处理后，得到包含不同频段的AM-FM分量，有些分量包含信号的固有信息，能够较大程度上代表原始信号，但是其余分量则不能反映原始信号特征或者直接是噪声干扰，因此有必要从各AM-FM分量中选出原始信号较多特征的主AM-FM分量，从而提高旋转设备的故障诊断精度。本文应用K-L散度来筛选EWT分解的主AM-FM分量，K-L散度也被称为相对熵，是信息论中较为重要的概念，用来度量2个概率分布间的差异程度^[12]。在本文的主AM-FM分量筛选中，与原信号K-L散度越小，则分量对原信号的近似度越强，反之，则与原信号相似度越小。K-L散度理论计算步骤如下：

1) 首先采取非参数估计法计算各分量概率分布：

$ \left\{ {\begin{array}{*{20}{l}} {p(x) = \frac{1}{{nh}}\sum\limits_{i = 1}^n k \left[ {\frac{{{x_i} - x}}{h}} \right],x \in {\bf{R}}}\\ {k(u) = \frac{1}{{\sqrt {2\pi } }}{{\rm{e}}^{ - \frac{{{u^2}}}{2}}}} \end{array}} \right. $

(4)

得到信号x=(x₁, x₂, …, x_n)的概率分布p(x)，以及信号y=(y₁, y₂, …, y_n)的概率分布q(x)，其中平滑参数h是给定的正数，k(u)为高斯核函数。

2) 根据2个信号的概率分布分别计算其K-L散度值，其中p(x)与q(x)的K-L距离为：

$ \delta (p,q) = \sum\limits_{x \in n} p (x){\rm{log}}\frac{{p(x)}}{{q(x)}} $

(5)

由于δ(p, q)≠δ(q, p)，无法满足距离的概念，故应用定义的p(x)及q(x)的K-L散度值来定量衡量其差异：

$ D(p,q) = \delta (p,q) + \delta (q,p) $

(6)

2 基于EWT-GG聚类的故障诊断方法 2.1 特征提取方法

对于由K-L散度筛选得到的主AM-FM分量，为较完备地获取其特征，选用基于样本熵及LZ复杂度的特征提取方法，组成二维矩阵进行聚类分析，从而提高聚类性能。

2.1.1 样本熵

样本熵是由Richman提出的一种度量时间序列复杂性的方法^[13]。对于n个数据点组成时间序列x={x(1), x(2), …, x(n)}，其样本熵计算步骤如下：

1) 构造一个m维向量X_m={X_m(1), X_m(2), …, X_m(n-m+1)}，其中X_m(i)={x(i), x(i+1), …, x(i+m-1)}。

2)定义X_m(i)与X_m(j)的距离D[X_m(i), X_m(j)]为2个向量中元素最大差值的绝对值，即：

$ D[{X_m}(i),{X_m}(j)] = {\rm{max}}[|x(i + k) - x(j + k)| $

(7)

3) 给定阈值r，统计每个向量距离值小于r的数目，计作B_i，并计算B^m(r)，其表达式为：

$ {B^m}(r) = \frac{1}{{n - m}}\sum\limits_{i = 1}^{n - m} {\frac{1}{{n - m - 1}}} {B_i} $

(8)

4) 将维数增加至m+1，统计每个向量距离值小于r的数目，计作A_i，并计算相应的A^m(r)。从而得到了2个序列在相似容限r下匹配m个点的概率B^m(r)，以及匹配m+1个点的概率A^m(r)。则样本熵的定义为：

$ {\rm{SampEn}}(m,r) = - {\rm{ln}}[{A^m}(r)/{B^m}(r)] $

(9)

可以看出，样本熵与m、r取值有关，本文依据Pincus研究结果^[14]，取m=2，r=0.2E_Std，其中E_Std为原始数据x(i)的标准差。

2.1.2 Lempel-Ziv复杂度

Lempel-Ziv复杂度算法由Lempel及Ziv给出，是一种时间序列复杂度分析方法^[15]。在计算LZ复杂度之前，首先将信号转化为二值符号序列，将序列中大于平均值的点赋值为1，而小于平均值的点赋值为0。

对于二值符号序列S(s₁, s₂, s_n)，在循环之前初始化字符变量P₀、Q₀为空字符，复杂度C(i)=0；随后，分别将P_t-1、Q_t-1与s_i连接成为新的字符串P_i、Q_i；判断若Q_i为P_i-1的子字符串，则表示暂时未出现新模式，复杂度C(i)保持不变，反之复杂度值加1；最后循环遍历n次后，将C(n)归一化得到C，即为该时间序列的LZ复杂度，其归一化表达式为^[16]：

$ 0 \le C = \frac{{C(n)}}{{B(n)}} \le 1 $

(10)

其中，B(n)为C(n)的上限，其表达式为：

$ B(n) = \mathop {{\rm{lim}}}\limits_{n \to \infty } C(n) = \frac{n}{{lb(n)}} $

(11)

2.2 GG模糊聚类

选用GG模糊分类算法对旋转机械振动信号的特征向量进行分析，该聚类算法引入基于模糊最大似然估计的距离测度，适用于非规则分布的数据分析，从而提高聚类性能^[17]，算法具体步骤如下：

1) 假定聚类样本为X={x₁, x₂, …, x_n}，其中每个样本x_j(1≤j≤n)包含d个特征属性，即x_i=(x_i1, x_i2, …, x_id)，对样本X初始划分为c类。设每个聚类中心的向量为V={v₁, v₂, …, v_c}，且隶属度矩阵为U=[μ_ij]_c×n，其中μ_ij∈[0, 1]表示第j个样本对第i类的隶属度；

2) 设定迭代的终止容限ε，有ε>0，并随机初始化隶属度矩阵U；

3) 计算聚类中心：

$ v_i^{(l)} = \sum\limits_{j = 1}^n {{{(\mu _{ij}^{(l - 1)})}^m}} {x_j}/\sum\limits_{j = 1}^n {{{(\mu _{ij}^{(l - 1)})}^m}} $

(12)

4) 计算各样本与类的模糊最大似然估计距离：

$ \left\{ {\begin{array}{*{20}{l}} {D_{ij}^2({\mathit{\boldsymbol{x}}_j},v_i^{(l)}) = \frac{{{{( {\rm{det}} (\mathit{\boldsymbol{A}}_i^{(l)}))}^{1/2}}}}{{{p_i}}}{{\rm{e}}^{{B_{ij}}}}}\\ {{\mathit{\boldsymbol{B}}_{ij}} = \frac{1}{2}{{({\mathit{\boldsymbol{x}}_j} - v_i^{(l)})}^{\rm{T}}}\mathit{\boldsymbol{A}}_i^{ - 1}({\mathit{\boldsymbol{x}}_j} - v_i^{(l)})} \end{array}} \right. $

(13)

式中：p_i为第i个聚类被选中的先验概率，有p_i=$\frac{1}{n}\sum\limits_{j = 1}^n {{\mu _{ij}}} $。m为加权指数，一般取2。A_i为第i个聚类的协方差矩阵，其表达式为：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{A}}_i^{(l)} = \sum\limits_{j = 1}^n {{{(\mu _{ij}^{(l - 1)})}^m}} {\mathit{\boldsymbol{K}}_{ij}}/\sum\limits_{j = 1}^n {{{(\mu _{ij}^{(l - 1)})}^m}} }\\ {{\mathit{\boldsymbol{K}}_{ij}} = ({\mathit{\boldsymbol{x}}_j} - v_i^{(l)}){{({\mathit{\boldsymbol{x}}_j} - v_i^{(l)})}^{\rm{T}}}} \end{array}} \right. $

(14)

5) 更新隶属度矩阵U：

$ \mu _{ij}^{(l)} = \frac{1}{{\sum\limits_{k = 1}^c {{{\left( {\frac{{{D_{ij}}({x_j},{v_i})}}{{{D_{kj}}({x_j},{v_k})}}} \right)}^{2/(m - 1)}}} }} $

(15)

判断条件‖U^(l)-U^(l-1)‖ < ε是否满足，如果满足，则终止迭代计算，否则令l=l+1。重复执行迭代V及U，直到满足条件，使目标函数J取得最小值：

$ J(\mathit{\boldsymbol{V}},\mathit{\boldsymbol{U}}) = \sum\limits_{i = 1}^c {\sum\limits_{j = 1}^n {\mu _{ij}^2} } D_{ij}^2 $

(16)

2.3 基于EWT-GG聚类的多分量故障诊断方法

采用Gath-Geva聚类的方法建立基于信息论及复杂度特征的轴承故障诊断模型，以便更加准确高效地识别故障类型，其故障诊断流程图如图 1所示。

本文故障诊断方法主要有以下几个步骤：

1) 首先采集轴承正常及一些常见故障下的振动信号，将其分为训练及测试样本集。对各样本集序列分别应用EWT进行分解，得到各频段的AM-FM分量;

2) 对每个分量进行基于K-L散度的主AM-FM分量筛选，得到的主分量与原始信号相似度最高。对主分量采用基于信息论及复杂度的特征提取方法，得到二维特征向量空间，从而完备地提取主分量的特征;

3) 采用GG聚类算法对训练样本的特征向量进行聚类，以此识别不同的轴承运行工况。同时，采用分类系数PC以及平均模糊熵CE对GG聚类的分类性能进行对比评价：

$ {{\rm{PC}} = \frac{1}{n}\sum\limits_{i = 1}^c {\sum\limits_{j = 1}^n {\mu _{ij}^2} } } $

(17)

$ {{\rm{CE}} = - \frac{1}{n}\sum\limits_{i = 1}^c {\sum\limits_{j = 1}^n {{\mu _{ij}}} } {\rm{ln}}{\mu _{ij}}} $

(18)

4) 计算测试样本特征向量与步骤3)GG聚类中心的欧式贴近度E(C_k, T)：

$ E({C_k},\mathit{\boldsymbol{T}}) = 1 - \frac{1}{{\sqrt n }}\sqrt {\mathop \sum \limits_{i = 1}^n {{[{C_k}({x_i}) - T({x_i})]}^2}} $

(19)

其中C_k为第k个状态，T为测试样本特征向量，元素x={x₁, x₂, …, x_n}。

5) 依据择近原则进行判断，贴近度最大则表明为同一类运行状态，从而实现了对旋转设备轴承振动信号的故障诊断。

3 轴承故障诊断实例分析

本文针对核电厂旋转机械轴承进行常见故障诊断方法研究，提高核电厂设备监测及诊断精确度。采用美国Case Western Reserve University电气工程实验室的滚动轴承实验数据，该轴承为SKF轴承，实验采用电火花技术加工轴承内圈、外圈以及滚动体单点故障，实验装置如图 2所示。采用的轴承数据样本集分为轴承正常、内圈故障、外圈故障以及滚动体故障4种类型，如图 3所示，损伤直径为0.177 8 mm，损伤深度为0.279 4 mm，轴承转速为1 797 r/min，采样频率为12 kHz。

3.1 主分量的选取

将正常、内圈、外圈以及滚动体故障样本数据分别分割为40组子样本集，样本长度为3 000，即0.25 s内轴承振动信号，并且交叉选取20组作为训练样本、剩余20组作为测试样本。首先对各训练样本进行EWT处理，得到各阶AM-FM分量，并计算其与原始信号K-L散度值以进行主分量筛选。为了对比筛选标准，本文计算分量与原始信号的相关系数。以内圈故障信号为例，将第1组训练样本经EWT分解，得到5个分量时频域结果如图 4所示。

由图 4可以看出，内圈故障的第1组训练数据经EWT分解后，得到的5个分量在时域上出现规律的周期成分，频域没有出现模式混叠现象，从而极大提高了后续故障诊断的精确度。接着把得到的5个分量选取K-L散度及相关性系数评价每种分量代替原始信号的程度，其归一化计算结果如表 1。

由表 1可以看出，分量4与原始信号的K-L散度较小，说明其与原始信号的关联性较大，能较大程度上代表原始信号，从而将分量1、分量2、分量3、分量5视为虚假分量滤除。同时，虽然分量4的相关性系数最大，但与分量2、分量3同属一个数量级，在后续筛选中易与虚假分量混淆。因此，本文采用基于K-L散度的主AM-FM分量筛选标准，可以更加准确识别虚假分量，提高后续GG聚类识别率。

3.2 模糊聚类算法对比测试

对于由K-L散度筛选的80组主分量训练样本数据，分别选用样本熵以及LZ复杂度进行特征提取，以便于完备地提取分量特征，同时将组成的二维特征向量作为输入，进行聚类分析。在GG聚类中，设定聚类中心数目为4，加权指数m为2，迭代终止容限设定为，得到的GG聚类结果如图 5(a)所示。

由图 5可以看出，GG聚类算法将特征向量分为4类，其中黑色的“O”代表聚类中心，其具体坐标如表 2所示。各聚类簇较集中于聚类中心，簇之间没有重叠部分，且等高线边界划分具有较大分类间隔。因此应用K-L散度提取信号主分量，并利用样本熵以及LZ复杂度作为特征向量进行GG聚类分析，极大地提高了聚类效果以及故障诊断精确度。

同时为了验证GG聚类算法中，采用基于模糊最大似然估计距离的优越性，本文应用FCM聚类以及GK聚类算法对同样的特征向量进行分析，得到的聚类效果如图 5(b)、(c)所示。

由图 5(b)、(c)可以看出，各聚类算法的聚类中心相较于GG聚类差别不大，但是GG聚类算法具有最大化分类间隔，其聚类等高线为任意形状，但是FCM聚类以及GK聚类等高线形状分别近似圆形以及椭圆形，表明GG聚类对数据分别最为自适应，能够区分任意分布的特征向量，而FCM聚类以及GK聚类算法仅对于圆形及椭圆形分布的特征向量具有较好的分析效果。

同时计算3种聚类算法的分类系数PC以及平均模糊熵CE以对其分类性能进行定量评价，其中分类系数越接近1，聚类效果越好，平均模糊熵越接近0，聚类效果越好。计算结果如表 3所示，可以看出，作为FCM算法的改进，GK聚类算法的分类性能优于FCM聚类，但是GG聚类算法的分类系数最高，达到了最好的分类效果，并且其平均模糊熵最低，说明该算法对4类故障信号的聚类效果最优。

3.3 信号分解方法对比测试

为进一步分析基于EWT信号分解方法的优越性，本文将信号分解方法替换为经验模态分析(EMD)进行对比测试。同样以内圈故障信号为例，将第1组训练样本经EMD分解，得到10个分量，其中前5阶分量时频域结果如图 6所示。

由图 6可以看出，EMD处理得到的前5阶内圈信号分量频率由高到低依次排列，其中第1阶分量频率最高。相比于图 4的EWT内圈信号分解图，第1阶分量所包含的频域较广，出现了较严重的模态混叠现象，不利于故障特征提取。后续故障诊断的其他流程不变，将训练样本应用GG聚类进行分析，得到的EMD-GG聚类结果如图 7所示。

由图 7可以看出，采用GG聚类方法具有较大的分类边界距离，但是EMD方法容易产生模态混叠现象，且受到虚假成分干扰较大，不能很好地分解出原始信号的主分量，导致聚类簇较为分散，不便于聚类边界的划分。再利用训练完毕的2种聚类器，计算EWT-GG聚类以及EMD-GG聚类的测试样本欧式贴近度，分别得到4类各20组测试样本的欧式贴近度，如图 8，每类样本4种图柱代表相对于4个聚类中心的平均欧式贴近度。

由图 8(a)可知，基于EWT-GG方法的测试样本与聚类中心的最大欧式贴近度更接近于1，且相对于其他聚类中心具有明显的区分度，便于测试样本的故障诊断；而在图 8(b)中，基于EMD-GG方法的最大欧式贴近度较小一些，同时正常与外圈故障样本的贴近度较为接近，不利于故障的识别。因此，本文方法在更加准确地提取出信号特征的同时，还能提高故障诊断精度。

4 结论

1) 应用EWT的信号分解方法能极大抑制模态混叠现象，相比于EMD更适合非平稳信号分析。

2) 采用K-L散度筛选出信号的主AM-FM分量，能较好地获取信号主要特征，便于提高后续故障诊断效率。

3) GG聚类对轴承常见故障诊断的分类性能优于其他模糊聚类算法，最后应用欧式贴近度进行故障识别，诊断过程简单且高效。