基于EEMD-Treelet和高斯过程的起重机齿轮故障诊断

李 勇，钱尼君，陈星斌，黄荔生，彭献永

(1.江西省特种设备检验检测研究院，南昌 330029；2.腾羿数据(上海)科技有限公司，上海 200240；3.上海交通大学 电子信息与电气工程学院，上海 200240)

0 引言

齿轮传动系统的主要机械振动通常归因于齿轮故障，包括磨损/折断的齿或不平衡的齿轮，并可能导致系统故障。众所周知，对于不同的故障类型，其影响将以不同的频率出现。准确而自动地评估其性能下降有两个重要方面：特征提取和故障分类器。Randall[1]提供了有关齿轮诊断技术的综述。目前，循环平稳性分析方法来诊断旋转机械中的齿轮已经被广泛采用。振动信号的频谱分析也已经用以检测齿轮中的剥落。

通常，在机器故障的情况下，测得的振动信号可能会出现不平稳的现象。因此，时频分析方法，如小波变换方法和希尔伯特-黄氏变换(HHT)被用来观察信号时变特性。Saravanan等[3]使用Morlet小波提取齿轮箱的故障特征作为SVM模型输入。Wu and Chan[4]通过小波变换和人工神经网络开发了齿轮故障识别模型。Peter等[5]利用小波分析和包络线检测进行轴承故障诊断。Peng等[6]比较了改进的希尔伯特-黄变换和小波变换在滚动轴承故障诊断中的性能差异。Lei和Zuo[7]用EEMD产生敏感的本征模态函数(intrinsic mode function,IMF)，用于旋转机械的故障诊断。Fan和Zuo[8]利用本征模态函数提取机器故障特征。Yu等[9]用希尔伯特-黄变换的时频熵方法进行齿轮故障诊断。Yang和Wu[10]用EEMD方法和PCA来诊断齿轮性能退化。

机器的故障信号通常显示出很高的复杂性，许多研究人员利用机器学习来识别机械运行状态故障。Samanta和Al-Balushi[11]提出了一种人工神经网络，使用时域特征对滚动轴承进行故障诊断。Wang等[12]使用时间序列分析和神经网络进行旋转机械故障诊断。Shafei等[13]将神经网络和模糊逻辑应用于齿轮的故障诊断。除了人工神经网络，支持向量机(SVM，support vector machine)[2]和决策树方法[14]也被用于不同的故障类型。Bordoloi和Tiwari[17]用SVM对齿轮频域数据进行多故障分类。他们通过在4个不同的故障条件下检查SVM模型的多类能力的性能。Rajeswari等[18]提出了一种基于遗传算法和反向传播神经网络的齿轮不同状态诊断混合框架。Sugumaran和Ramachandran[19]使用齿轮振动信号的时域统计和直方图特征，作为SVM分类器输入，找到了最佳故障特征。Zhang Ying等[2]提出了基于多种特征的故障诊断方法，他们用核主成分分析和粒子群优化支持向量机(来定位轴承故障。

高斯过程(GP)是一种基于核的非线性学习方法，可以处理带有噪声观测值和系统不确定性的建模问题[20]。它已经被广泛的应用于动作识别等。它具有严格的统计理论基础和很好的适应复杂问题的能力。与ANN和SVM相比，它既保持良好性能又具有灵活的非参数推理能力的条件。在GP模型中，关键问题是估计一组适当的超参数，该参数会使边际似然函数最小化。为了提高GP模型参数寻优的准确性[15]，用一种细菌觅食优化(BFO)算法[16]来优化GP模型的超参数。

为了提高起重机齿轮传动系统故障识别性能，提出一种基于EEMD-Treelet的自动特征选择方法，用于诊断齿轮故障。首先，采用EEMD方法来提取通常隐藏在振动信号中的与故障相关的特征。然后，通过Treelet变换对原始特征进行降维，降维后的特征可以表示齿轮的主要故障特征。所得的特征向量作为BFO-GP模型输入，自动识别齿轮故障。

1 特征提取与降维

1.1 边际谱

EEMD方法[3]的本征模式函数是单振荡成分。因此，IMF瞬时频率计算是可行的。这是建立时频分布的必要和关键步骤。为了建立信号的时频分布，第k个IMF可以用其包络和载波部分表示为:

(1)

其中：包络线Ak(t)表示随时间变化的幅度项，可以通过复杂的三次样条拟合计算得出，而被测信号x(t)的时频分布H(ω,t)可以表示为：

(2)

此外，边缘频谱S(ω)可以通过积分沿时间轴的时频分布表示如下：

(3)

其中：Hk(ω,t)表示第k个IMF的时频分布，T表示信号长度，Ω表示所选IMF的索引集。

1.2 Treelets降维

Treelets降维算法是一个自适应的降维算法,它在数据每层结构上将两列最相似的数据用一组逼近信号与一组细节信号来代替，采用PCA方法进行对这两列数据进行分析[21]。

Treelets算法的结构可以表示为一个树结构。图1即为一个5层的树结构信号：部分(主成分结构)代表的为逼近信号，(次成分信号)代表与主成分正交的细节信号。

图1 五层Treelets结构

2 BFO-GP分类器

2.1 GP分类器

一个高斯过程可以通过其均值函数m(x)和它的协方差k(x,x')确定。它的超参数最佳估计可以在对数边际似然函数最大情况下获取。也就是说，可以通过对目标函数进行边际化，获得高斯过程超参数的概率，使用共轭梯度方法优化超参数。高斯过程的对数边际似然函数定义为：

logp(y|x,Θ)=

(4)

超参数的选择对于高斯过程分类器的整体分类性能至关重要。它的超参数{μ,σ,l1,l2,…,ld}可以在很大的值范围内变化。通常的共轭梯度优化方法的性能或准确性取决于初始数值的选择，并且相对于小扰动和局部最优而言可能不稳定。由于选择超参数会极大地影响高斯过程模型的泛化性能，因此进化算法可能是一个可好的选择。在这个工作中，细菌觅食优化算法被用作高斯过程超参数优化算法。细菌觅食优化算法的适应度函数采用式(7)中定义的对数边际似然函数。

2.2 细菌觅食算法

BFO是一种自然启发式优化算法[16]，它能够找目标函数的最小值，并且不需要计算目标函数的梯度。BFO模仿了在实际细菌系统中观察到的4个主要机制：趋化性，群聚，繁殖和迁徙，以解决非梯度优化问题。

1)趋向性操作：

大肠杆菌在整个觅食过程中有两个基本运动：旋转和游动。旋转是为找一个新的方向而转动，而游动是指保持方向不变的向前运动。BFO算法的趋向性操作就是对这两种基本动作的模拟。通常，细菌会在食物丰盛或环境的酸碱性适中的区域中较多地游动，而在食物缺乏或环境的酸碱性偏高的区域则会较多地旋转，即原地不动。

假设要找到目标函数J(θ)的最小值，其中θ∈Rp(即θ是实数p维向量)，并且没有测量值或梯度▽J(θ)的数学解析描述。假设θi(j,k,l)代表第j个趋化性的第i个细菌,第k个繁殖和第l个消除-扩散事件。C(i)表示在翻滚指定的随机方向上执行的步长。在计算趋化性中，细菌的运动可以表示为：

(5)

其中：Δ表示随机方向上的向量，其元素位于[-1,1]。

2)群聚性操作：设目标函数的变化量如下：

(6)

其中:JCC(θ，P(j,k,l))表示目标函数J(θ)的目时变值，S是细菌总数，p是每个细菌中要优化的变量数量.θ=[θ1,θ2,…,θp]T是p维搜索域中的一个点。dattrac,wattrac,hrepel,wrepel是固定学习系数。

3)繁殖性操作：

生物进化过程一直服从达尔文进化准则，即“适者生存、优胜劣汰”。BFO算法执行一段时间的觅食过程后，部分寻找食物源能力弱的细菌会被自然淘汰，而为了维持种群规模不变，剩余的寻找食物能力强(适应度值低)的细菌会进行繁殖。在BFO算法中模拟这种现象称为繁殖性操作。

(7)

4)迁徙性操作：

实际环境中的细菌所生活的局部区域可能会发生逐渐变化(如食物消耗殆尽)或者发生突如其来的变化(如温度突然升高等)。这样可能会导致生活在这个局部区域的细菌种群被迁徙到新的区域中去或者集体被外力杀死。在BFO算法中模拟这种现象称为迁徙性操作。

2.3 BFO-GP 分类器

基于EEMD-Treelet和细菌觅食优化算法-高斯过程的齿轮故障诊断步骤如下：

步骤(1):数据准备，从原始振动信号中提取特征作为特征向量，标签这些特征向量。使用Treelet消除原始特征向量的冗余；

步骤(2):GP分类初始化，使用常数均值函数和平方指数协方差核函数；

步骤(3):使用BFO训练GP分类器：选择不同的超参数选择，使用式(4)计算相应的适应度函数值。设置和使用GP分类器：使用优化的超参数集来设置GP分类器。对于新的分类诊断问题，按照步骤(1)提取特征，并计算分类概率。

3 实验结果与分析

实验设备由起重机齿轮箱组成，模拟不同齿轮故障的工作条件。安装在驱动轴上的1号齿轮与安装在从动轴上的2号齿轮啮合。驱动轴连接到由伺服电动机驱动的挠性联轴器。两个轴由四套轴承支撑。机械振动由加速度计测量。表1显示了此实验中6种齿轮故障情况。驱动马达的旋速为1 500 RPM，驱动轴和从动轴的转动频率分别为25 Hz和11.3 Hz。根据轴的转速和相应的齿轮啮合频率，将数据采集装置的采样率设置为6 400 Hz。此实验中，每个实验案例有50组振动信号，总共记录了300组信号。

表1 实验齿轮故障类型

3.1 特征提取

这里以案例3的信号为例进行说明。将所有振动信号用EEMD分解为本征模态函数 (C1,C2,…,Cn)，如图2所示。然后，计算每个IMF的瞬时频率，表2列出了每个瞬时值的平均值。结果显示本征模态函数C7和C8与主动轴和从动轴的转动频率有关，分别为25 Hz和11.3 Hz。研究结果也显示啮合频率575 Hz位于本征模态函数C2和C3的带宽之内。由于在啮合频率附近的振动信号分量通常由与机器故障有关，因此提取C2和C3相关的特征，更具实际意义。为了提取与故障有关的特征，对C2和C3进行包络分析如下。

图2 案例2中振动信号的EEMD分解产生的IMF

表2 每个IMF的瞬时频率分布的平均值

步骤1：将与啮合频率相关的C2和C3表示为S，提取S信号的包络信号Envs。

步骤2：在Envs上应用EEMD和希尔伯特变换，得到S的希尔伯特时频分布频谱。

步骤3：沿时间轴对步骤2中的频谱进行积分，得到S的希尔伯特边际包络谱。

步骤4：利用三次样条拟合希尔伯特边际包络谱的包络信号。

包络分析结果显示案例3的包络信号在0～50 Hz范围内比案例2的能量集中度更高，这是因为案例3的齿轮故障水平比案例2更为严重。齿轮故障在25 Hz(1号齿轮的旋转频率)处观察到峰值。还要注意的是，驱动齿轮转速峰值的大小取决于齿轮损坏的程度，这意味着严重磨损的牙齿(如案例3)的包络信号具有比轻微磨损的情况更集中的能量。另外，由于断齿的劣化程度比轻度和严重磨损的牙齿更严重，因此在25 Hz的情况下断齿(案例4)的包络谱峰值高于案例2和案例3。2号齿轮的主要特征频率是11.3 Hz。在齿轮的边缘希尔伯特包络谱中观察到严重不平衡情况的能量集中度高于轻微不平衡情况和正常情况(11.3 Hz)的能量集中。这是因为严重不平衡情况下的故障水平比轻微失衡和正常情况下要严重。

根据以上分析，从振动信号中提取的特征如表3所示。

从案例1-5中随机选择30个振动信号，形成样本集，计算表3的特征，获得21维振动信号状态特征。Treelet变换算法层次设置为2，将这21维特征降为为2维。图3显示了经Treelet变换后齿轮特征数据分布。显然，Treelet方法可以将正常数据与故障数据分开。每种案例的数据都在特征空间中较好地聚集。仅有一对不平衡数据特征存在重叠。但是，BFO-GP仍然可以获得95%的诊断精度，这是令人满意的。对于严重磨损的齿轮和断齿的齿轮，也发现了相同的情况。

图3 特征空间中的齿轮数据分布

表3 每个IMFS的瞬时频率分布的平均值

BFO-GP输出特征空间中每个数据点的分类概率。以0.5概率轮廓线作为决策线，将特征空间分为两类。图4显示了齿轮诊断结果。显然，在使用Treelet进行特征降维之后，正常齿轮和退化齿轮数据集在特征空间中得到很好的分类，并且彼此之间有着很好的区分。分类结果是令人满意的。基于BFO-GP的方法是适合齿轮故障诊断。

图4 齿轮诊断结果

进一步，C1-C55案例的故障诊断其结果如表6所示。显然，BFO-GP的诊断精度高于其它模型如GA-GP，PSO-GP和BFO-SVM，但是BFO-GP消耗更多的运行时间成本。测试结果显示对于GP的超参数优化能力，BFO比GA、PSO等其它算法表现得更好。

表6 不同案例的诊断结果

4 结束语

基于EEMD-Treelet-GP的齿轮故障诊断结果发现EEMD-Treelet算法能够学习到更多振动信号精细特征。由EEMD-HT获得齿轮振动信号的希尔伯特频谱，更多齿轮故障特征被提取。通过Treelet变换对这些特征进行筛选，获得线性不相关特征。这些特征输入进BFO-GP模型，实现齿轮故障诊断与退化检测。