基于数学形态学分形理论与PNN的轴承状态监测与故障诊断

林 懿， 龙 伟， 李炎炎， 张庆华

(四川大学机械工程学院， 成都 610065)

1 引 言

滚动轴承作为旋转机械中应用最广泛的关键部件，因其恶劣的工作环境，导致故障频发，给工业安全生产带来了严峻挑战.《中国制造2025》已明确将提高大型设备的可靠性与安全性列为重要发展方向.因此在工业4.0背景下，对滚动轴承进行基于数据驱动的实时监测与故障诊断具有重要的理论和现实意义.

机电设备运行时产生的振动信号中包含着丰富的轴承状态信息，且抗干扰能力强，常用于设备故障诊断[1].但由于滚动轴承工况复杂、激励源多、耦合性强，振动特性表征模糊，诊断效果往往不佳，因此故障特征的精确提取是实现高精度故障诊断的重要前提[2].经典的时频域处理法存在特征提取混肴、表征性差、自适应性缺乏等缺点，不适用于智能化监测与诊断[3].近年兴起的基于深度学习的轴承故障诊断方法，虽然取得了不错的诊断效果，但因其硬件成本高、运算量大、实时性差等问题，在工业生产方面同样缺乏使用环境[4].分形几何学是一门以不规则几何形态为研究对象的几何学,通过提取复杂信号的分形维数，定量表征信号的复杂性和非线性.李洋等[5]将分形几何理论引入轴承故障诊断中，成功提取了不同状态轴承的特征信息.在多种分形理论中，基于数学形态学的分形维数计算方法以其鲁棒性强、计算量少等优点被广泛应用于信号处理方向[6].李兵等[7]首次使用数学形态学理论计算一维振动信号的分形维数，成功识别轴承四种基本状态，但由于振动信号中噪声的影响，加之表征维度低，难以区分不同损伤程度的轴承故障.

在实测振动信号中，大量的背景噪声因其随机性，不具备明显的分形特性，必须对原始信号降噪提纯，才能保证分形维数的准确性[8].经验模态分解(Empirical Mode Decomposition，EMD)通过迭代分解将原始信号重构为一系列本征模态函数(Intrinsic Mode Function， IMF)达到降噪目的，因其自适应性、正交性及非线性被广泛应用于信号处理领域，但“端点效应”和“模态混叠”的存在影响了其降噪性能.尽管国内外学者对其进行大量的改进与优化，迭代次数多、运算量大及实时性差等问题的存在使其难以适用于轴承状态的在线监测[9].局部均值分解(Local Mean Decomposition，LMD)于2005年被Jonathan S. Smith提出并首次运用于EEG信号中，它通过自适应的将复杂原始信号分解成一系列纯调频信号与包络信号的乘积函数(PF分量)降噪提纯，不仅延缓了EMD“模态混叠”与“端点效应”的缺陷，并且分解速度更快.张亢等[10]通过LMD分解轴承振动信号降低噪声影响，并提取相关性系数最大分量的分形维数作为分类依据，基本实现了轴承工作状态的识别，但轴承工况复杂，使用单一特征表征轴承状态缺乏空间性与全面性.金成功[11]则通过EMD的改进算法CEEMDAN降低振动信号的噪声，并提升特征维度，但结果同样只能识别4种基本状态，且算法运行时间较长.直到齐咏生等[12]使用LMD算法分解振动信号，通过提取前三个PF分量，将其重构为特征向量进行轴承故障诊断，才对同种故障不同损伤程度状态的识别取得了较好的效果，不过，因为LMD算法过量的筛选迭代，导致PF虚假分量的存在，使得第三个PF分量与原始信号的相关性较差，难以精确表征，因此，在LMD算法的改进与轴承故障特征的选择上，仍有进一步提升的可能.

将提取的特征向量与模式识别理论结合，可以实现轴承的在线智能监测与故障诊断[13].其中，聚类分析难以对高维数据分类识别[14]，SVM[15]算法冗余，缺乏效率，且运算时间长.而概率神经网络(Probabilistic Neural Network，PNN)是一种依据贝叶斯决策理论构建的智能分类器.因其训练简洁、收敛快以及容错高的优点，适用于数据的实时处理[16].

就上述不足，本文通过构建“三次样条插值”与“迭代终止条件”完善LMD算法的不足，优化算法实时性与PF分量的精确性；同时，结合分形盒维数的基本原理，首次引入物理量“形态学覆盖面积”作为特征向量，代替相关性系数较差的PF分量，仍旧从多维度表征轴承状态，提升表征精度.

2 背 景

2.1 数学形态学分形维数

分形维数的计算是用不同尺度表征分形集，以定值其非线性与复杂性，同时量化其结构特征.在一维离散数据的背景下，数学形态学维数的计算是基于盒维数与多尺度数学形态学的覆盖思想，对一维振动信号在不同尺度下进行腐蚀(Θ)和膨胀(⊕)运算，具体原理如下.

设f(n)与g(m)分别为定义在F={0，1，2，…，N-1}和G={0，1，2，…，M-1}上的离散函数，且N≥M，则数学形态学中腐蚀滤波器(Θ)和膨胀滤波器(⊕)的定义如下.

(fΘg)=min[f(n+m)-g(m)]

(1)

(f⊕g)=max[f(n-m)+g(m)]

(2)

其中，f为输入信号；g为结构元素.

在尺度ε(整数，且1≤ε≤εmax)下的结构元素定义为ε-1次腐蚀或者膨胀运算，即

gΘε=gΘgΘ…Θg

(3)

g⊕ε=g⊕g⊕…⊕g

(4)

则尺度ε下，对信号的覆盖面积为

(5)

Maragos于1993年提出，当ε→0时有

(6)

其中，ε=1, 2 …,εmax；DM为Minkowski-Bouligand维数；c为常数；同时，对{ln(Ag(ε)/ε2)，ln(1/ε)}使用最小二乘线性拟合，所得直线斜率则为分形维数DM的极限近似值.

2.2 局部均值分解算法及其改进

2.2.1 基本算法 局部均值分解的本质就是从信号中逐步分离出调频信号和调幅包络信号的过程.通过多次筛选，自适应的将PF分量由高频到低频逐渐分离，直到剩余信号为常数或不再包含振荡为止.具体过程的步骤如下.

(1) 计算原始信号x(t)相邻局部极值点的均值，如式(7).

(7)

将各mi线性连接，再用滑动平滑法处理，得到局部均值函数m11(t).

(2) 局部包络函数，计算方法如式(8).

(8)

用局部均值函数相同的方法处理，得到包络估计函数a11(t).

(3) 将m11(t)从x(t)中分离，得到残差信号h11(t).

h11(t)=x(t)-m11(t)

(9)

(4) 振幅解调.

s11(t)=h11(t)/aii(t)

(10)

我们将s11(t)视为原始信号，重复步骤(1)和步骤(2)得到包络估计函数a12(t)，若a12=1，则判定s11(t)为纯调频函数.否则，以s11为原始信号，重复上述步骤，直到满足条件得到纯调频信号s1n(t).不过在实际应用中，通常设定阈值α，使|a1n-1|<α时，终止迭代.

(5) 所有包络估计函数的乘积为调幅包络信号a1(t)，即

(11)

(6) 求出乘积函数PF1(t).

PF1=a1t(t)s1n(t)

(12)

(7) 用x(t)减去PF1(t)得到u1(t)函数，重复上述步骤，直到uk(t)变成常数或不再振荡.

因此，通过局部均值算法的分解，原始信号x(t)如式(13)构成.

(13)

2.2.2 弊端及改进 尽管局部均值在一定程度上改善了经验模态分解理论支撑缺乏、运算量大以及模式混淆等缺点，但其本身仍存在大量迭代“筛选”，实时性不足.比如，在计算局部均值函数与包络估计函数的过程中，连续的滑动平移动作会增加运算复杂度，损耗大量的计算内存；其次，过多的筛选次数不仅浪费时间，也容易产生相关性系数极差的虚假PF分量.

针对上述问题，我们提出了以下两个改善方案.

(1) 采用“三次样条插值法”代替滑动平移法，只需一次插值，即可获得平滑的理想结果，减少了算法的复杂度，改善了实时性.

(2) 设定合理的阈值γ，作为PF分量的选择标准，当其与原始信号的相关性系数符合ρ<γ，则判定为虚假分量，及时终止迭代，不仅降低了算法运行时间，同时保证了分量的合理性，提高了特征提取的准确率.相关性系数ρ计算公式如下.

(14)

其中，x(t)为原始信号；p(t)为PF分量时序信号，并且经过多次调试，选择0.2作为相关系数阈值，既保证了PF分量的个数，又极大的避免了虚假分量的出现.

2.3 特征提取

对冗余过多的特征向量智能识别，不仅增加硬件负担，不利于算法实时性，还会影响算法识别的正确率.因此，只提取相关性系数最大的两个分量的分形维数作为特征向量是合理的.其次，形态学覆盖面积，作为描述分形集整体空间大小的物理量，同样可定量且直观地描述信号的特性.其理论依据基于分形盒维数的覆盖方法，即假设有一个面积为S的未知平面，用边长为ε的正方形填充该区域，需要正方形的个数为N(ε)，则有下列推论：

(15)

两边取对数：

lnN(ε)=lnS+2 ln (1/ε)

(16)

得到未知平面的维数为：

(17)

因此，形态学覆盖面积S，相较于分形维数D，同样可以作为轴承分类依据，此外，在基于数学形态滤波的分形维数提取中，因腐蚀滤波器与膨胀滤波器的特性，形态学覆盖面积与形态学分形维数同样具有平移不变性，因此，将式(5)中的形态学覆盖面积Ag作为特征向量引入轴承故障诊断也具有科学的理论依据.

综上所述，本文提取相关性系数最大的两个PF分量的数学形态学分形维数，结合信号的数学形态学覆盖面积共同构建特征向量，输入概率神经网络中以达到状态识别的目的.

2.4 算法步骤

(1) 信号去噪.使用改进局部均值算法，对采集到的原始信号x(t)去噪，保留相关性系数最大PF1，PF2分量.

(2) 特征提取.提取PF1，PF2的数学形态学分形维数T1，T2并对信号x(t)的形态学覆盖面积Ag照盒维数覆盖思想取对数，得到第三个特征向量T3.

(3) 状态识别.将特征向量T1，T2，T3输入已经训练好的概率神经网络模型中，完成对轴承工作状态的识别.

3 实 验

为了验证该方法在改善实时性与准确率方面的有效性，采用CWRU(美国凯斯西储大学)轴承数据中心[17]的实测数据.本文选取滚动体、外圈和内圈三种故障进行分析，其中，为了更直观体现所提方法的实用性，滚动轴承、内、外圈分别采用了损伤直径为0.007 in和0.021 in(1 in=2.54 cm)的故障信号.因此，包含正常状态在内，共采用了七种类型共819个样本的振动信号，每种类型的样本数量为117，单个样本中含有1 024个样本点.数据其他信息为：转速为1 750 r/min，采样频率为12 kHz，故障深度为0.011 in.

3.1 信号分解

滚动轴承工作状态发生改变时，其振动信号也会产生相应变化.如图1所示，正常工作状态的滚动轴承，振动信号随机性强，呈现出非线性、非平稳的特点；滚动体故障时，有一定频率的冲击，但整体随机性较强；内、外圈故障状态下，时域信号波形呈现周期性冲击，特征明显.因此，轴承四种状态中，正常状态与滚动体故障状态相似，冲击特性模糊、随机性强，且振幅较小，特征混杂，两者不易区分.而内、外圈故障状态规律强、振幅较大，具有明显的冲击特性，但两者特性相似，仍然不易区分.

图1 滚动轴承四种状态时域波形图Fig.1 Time domain waveform diagram of rolling bearing in four states

振动特征容易淹没在摩擦、刚度、载荷等复杂工况环境带来的噪声中，因此，需要对原始振动信号分解提纯.以0.07 in内圈故障为例，分别使用改进局部均值分解(ILMD)与标准局部均值分解(LMD)处理原始信号，分解结果如图2和图3，且两种算法各PF分量的相关系数ρ如表1所示.可以看到，因为阈值γ的缘故，当ILMD分解出的PF分量与原始信号的相关系数ρ小于0.2时，迭代终止，不仅减少了运行时间，也有效的避免了虚假分量的产生；此外，三次样条插值法的采用，一定程度的改善了“端点效应”的影响，并且显著提升了前两阶PF分量的相关系数，最大程度保留了原始信号的振动特征，提高了特征提取的精度.同时，为更进一步验证所提算法的优势，将其与经验模态分解(EMD)、自适应噪声的完全集合经验模态分解(CEEMDAN)比较.

图2 ILMD时域图Fig.2 Temporal waveform of ILMD

图3 LMD时域图Fig.3 Temporal waveform of LMD

表1 PF分量相关系数对比

通过对0.007 in轴承内圈故障信号的处理，分别对比前两阶分量相关性系数ρ、程序运行时间t，为避免偶然因素，结果为重复20次实验求取平均值.具体对比结果如表2所示，显然，ILMD前两阶分量的相关系数比EMD、CEEMDAN大，对轴承状态信息表征更好，故障特征提取更精确.此外，CEEMDAN算法复杂，2.94 s的运行时间，远超ILMD算法与EMD算法.综合比较，ILMD算法在表征精确性与运行实时性两方面，都优于EMD及CEEMDAN算法，因此，选择ILMD对原始振动信号分解是合理且具有优势的.

表2 各算法比较

3.2 特征提取

在信号多尺度形态学分析的过程中，结构元素G和最大形态学尺度εmax的选取非常关键.结构元素G的类型和长度是影响形态学覆盖面积与分形维数计算的重要参数，将“零高度”的扁平型G作为单位结构元素，一方面保持了信号的原始特性，即其幅值免受结构元素高度的影响；另一方面，削减了部分计算量，增加算法实时性[18].因此，选取的单位结构元素为[0 0 0].而最大形态学尺度εmax的选取也会对结果产生重大影响，一般而言，εmax的标准为[1,N/2]之间的整数，通过文献参考与多次调试，再权衡算法实时性与准确率之后，选择εmax=8最为合适.

使用前文提到的ILMD算法分解七种工作状态下的轴承振动信号，计算第一个PF分量的数学形态学的分形维数，并在一维空间中进行轴承状态识别.当识别类型为正常状态及0.007 in损伤直径的三种故障时，比较结果如图4所示.可见不同状态的分形维数呈现出明显的分形区间，故障诊断效果较好.当识别类型增加了0.021 in损伤直径的三种故障后，比较结果如图5所示，可以看到，正常状态与滚动体故障有一定的分形区间，分形特性明显，区分容易，但不同损伤程度滚动体故障的分形维数交叉重叠，难以区分；此外，对于0.007 in与0.021 in 的内、外圈故障，冲击特性模糊，辨识不易，需要进一步的特征提取.综上可知，利用分形特性对轴承工作状态进行识别的方法是可行的，并且对不同轴承故障的识别效果良好，但对同种故障但不同损伤程度的状态分辨效果差.因此，使用单一维度特征表征轴承状态缺乏空间性与全面性，可以增加特征维数以更好的表征轴承状态.

图4 不同故障一维PF分量分形维数Fig.4 Fractal dimension of one-dimensional PF component of different faults

图5 七种状态下一维PF分量分形维数Fig.5 Fractal dimension of one-dimensional PF component in seven states

图6为使用两个PF分量分形维数的二维散点图，由于向量维度的增加，尽管“特征混叠”现象有所改善，但0.021 in损伤直径的内、外圈故障及不同损伤直径的滚动体故障的特征仍然有部分混叠交联，不易区分.图7为三个PF分量的三维散点图，相较于二维散点图，“特征重叠”现象几乎没有改善.一方面是因为单一尺度特征向量的表征缺少完备性，另一方面则是因为第三阶PF分量与原始信号的相关性系数较小，难以精准表征原始信号的分形特性.因此，需要采用其他特征替代第三阶PF分量以改善识别效果.

图6 七种状态下二维PF分量分形维数Fig.6 Fractal dimension of two-dimensional PF component in seven states

图7 七种状态下三维PF分量分形维数Fig.7 Fractal dimension ofthree-dimensional PF component in seven states

作为描述分形集整体空间大小的物理量，“形态学覆盖面积”也可以定值定量的表征轴承工作状态，其理论依据已在前文给出.因此，将形态学覆盖面积的对数作为第三维特征向量引入聚类识别中，其结果如图8所示.相较于前文三种方法，该方法类聚合度高，“特征混叠”现象改善明显，不仅在识别轴承不同工作状态方面表现优异，在区分同种故障的不同损伤程度方面也表现良好.

图8 引入形态学覆盖面积的三维特征Fig.8 Three-dimensional features of the morphological coverage area introduced

上述实验证明了本文提出的特征提取方法具有显著优势，并为轴承状态监测与故障诊断提供了一种有效可行的新指标.此外，对不同损伤程度状态的精确识别，也有助于轴承寿命的预测.

3.3 智能识别与分类

为减少人工干预，实现在线监测，引入概率神经网络(PNN)对轴承状态智能分类.实验时，将不同状态的样本随机划分为两份，第一份约占样本数3/4，作为训练样本，另一份为测试样本，因此在7种不同状态下，训练集共609组，测试集共210组.此外，在PNN网络训练中，较小的spread系数会优化分类结果[19]，经过多次调试，该参数设定为0.01.

确定参数后，用本文提出的三维特征向量构建输入矩阵，并训练生成神经网络模型，以实现对轴承状态的识别分类.为保证分类结果的准确性与科学性，随机抽取样本并重复实验20次，统计其平均准确率.图9为20次重复实验中分类效果最优的一次，7种状态下的识别准确率达到100%，不仅能够识别不同状态的轴承，并且成功对同种故障不同损伤程度的轴承状态进行了分类.为进一步体现本文算法在实时性与准确率方面的提升，将其与不同算法比较，结果如表3所示.

图9 七种状态下的识别结果Fig.9 Recognition results in seven states

表3 不同算法各状态的识别率与识别时间

综上所述，针对滚动轴承状态的在线监测与故障的智能诊断问题，本文提出算法在监测实时性与诊断准确度上面都有大幅改善，为滚动轴承的智能管理提供了一种有效的新方法.

4 结 论

本文提出了一种基于ILMD与形态学分形理论的特征提取方法，并结合PNN成功实现对轴承状态的监测与故障诊断.该方法首先通过ILMD分解原始振动信号，抑制噪声，提高特征维度；其次，使用数学形态学覆盖求取前两阶PF分量的分形维数，并结合盒维数分形理论，引入形态学覆盖面积构成三维特征矩阵；最后，通过CWRU中心轴承真实数据构建PNN神经网络模型完成识别，并与其余算法比较，验证了该方法的可行性与优异性.结果表明，本文所提方法不仅能识别不同状态下的轴承，也可精确表征不同损伤程度的轴承.其平均识别时间只要0.21 s，平均识别准确率高达99.6 %，远超其余算法，为滚动轴承状态在线监测与智能故障诊断提供了一种高效、科学的方法.