滚动轴承故障诊断的最优加权路图GFT冲击提取方法
王好将 于德介 高艺源
摘要：滚动轴承发生局部故障时将产生由共振频率调制的周期瞬态冲击，有效提取冲击特征是诊断滚动轴承故障的关键。图信号处理方法（Graph Signal Processing，GSP）是基于圖谱理论发展起来的新研究领域，将振动信号转换为图信号进行分析，能有效揭示振动信号特性。对高斯函数加权下的路图拉普拉斯矩阵进行特征分解，发现代数连通度（Algebraic Connectivity）以内的特征向量存在明显的冲击，因此提出利用代数连通度以内的特征向量结合逆图傅里叶变换（Graph Fourier Transform，GFT）重构故障信号中冲击分量的方法。高斯加权函数中的热核宽度决定冲击特征向量的分布，直接影响重构结果，为解决热核宽度的选择问题，提出结合粒子群算法（ParticleSwarm Optimization，PSO）确定最优热核宽度;然后利用最优特征向量组重构冲击信号，并进行包络解调;最后实现滚动轴承故障的有效诊断。算法仿真和应用实例表明，基于最优加权的路图GFT方法能有效地重构滚动轴承故障冲击特征，诊断故障类型。
关键词：故障诊断;滚动轴承;冲击提取;路图;高斯加权函数
中图分类号：TH165+。3;TH133.33;7N911.7文献标志码：A 文章编号：1004-4523（2020）03-0604-10
DOI：10.16385/j.cnki.issn.1004-4523.2020.03.020
引言
滚动轴承作为旋转机械重要的承载零部件，其发生故障时直接影响整台机械的性能，为避免重大故障的发生，需要监测其运行状态。轴承内圈、外圈和滚动体存在局部故障时，伴随着轴承的周期性旋转，会以一定故障通过频率产生由共振频率调制的周期瞬态冲击。因此，滚动轴承的故障信息存在于周期性冲击成分中，从轴承故障振动信号中提取周期性冲击特征是有效诊断轴承故障的关键。但在实际工程中，受到工作环境、振动传输路径和多振动源相互耦合等因素影响，导致轴承振动信号中往往存在严重的背景噪声干扰及工频振动、谐波振动等强故障特征信号，此时直接对故障信号进行解调分析无法发现弱故障特征，容易造成故障的漏诊。因此，寻找行之有效的冲击特征提取方法一直是滚动轴承故障诊断研究的重点。
共振解调方法是提取滚动轴承故障冲击特征的有效方法。陈祥龙等提出了基于平方包络谱相关峭度指标的最优共振解调方法，稳健地确定最优频带中心频率和带宽，准确提取循环瞬态冲击信号。Antoni提出了基于FIR带通滤波器的快速谱峭度算法（SK），通过设计二叉树结构的带通滤波器对信号进行滤波，对各层滤波结果计算峭度值得到快速谱峭度图，根据峭度最大原则选取最优解调分析频带，提取冲击特征信号。目前谱峭度方法已广泛应用于滚动轴承故障诊断。
近年来，以图信号为研究对象的图信号处理方法开拓了信号处理方法新领域，逐渐引起国内外学者的关注。图谱理论（Graph Spectrum Theo-ry）是图信号处理方法的基础，其主要内容是结合代数方法分析图信号。目前，图谱理论已经广泛应用于机器视觉、网络学习、多传感器融合等领域，不但在故障诊断领域中的应用研究则相对较少。图傅里叶变换方法将图信号正交投影于拉普拉斯矩阵的特征向量上，为图结构数据的分析提供了类似“频率”的概念，从而将图信号研究从顶点域拓展到图谱域，并为图信号顶点域中卷积、平移及调制运算的实现提供了可能。为了将图信号处理方法引入到振动信号分析中，需要将振动信号时间序列转换为合适的图信号。路图仅考虑相邻时间点间存在边连接关系，与时间序列存在一一对应关系。研究发现，无权路图的拉普拉斯矩阵特征向量与离散余弦变换的基函数相同，路图傅里叶变换也与经典FT结果有对应关系。
Ou等研究发现，在欧氏距离加权下，路图傅里叶变换将冲击成分投影在高阶次区域，并利用这一特点实现了冲击特征提取。然而，欧氏距离加权的形式与一般加权的定义相违背，即距离越远权值越小，另外冲击成分出现在高阶次的现象也未从理论上解释，限制了该方法的进一步应用。本文研究发现，在高斯函数加权下，与路图拉普拉斯矩阵零特征值对应的特征向量呈现冲击现象，这一特点与图上连通分量相关联，对应的路图傅里叶变换将冲击成分投影在低阶次区域，尤其是代数连通度（Alge-braic Connectivity）以内的低阶次区域。针对上述特点，可利用代数连通度以内的拉普拉斯矩阵特征向量通过路图傅里叶逆变换重构冲击信号。此外，热核宽度是高斯函数中的重要参数，影响零特征值数目及冲击特征在阶次域的能量分布，不合理的热核宽度将导致无法重构冲击特征。因此，选择合适的热核宽度值是有效重构冲击信号分量的关键。
基于高斯函数加权下路图傅里叶变换的特点，本文提出一种结合粒子群优化算法的最优加权路图傅里叶变换冲击提取方法。该方法先将振动信号转换为路图;再以重构冲击信号包络谱故障频率处的幅值为适应度函数，通过粒子群算法寻找最优热核宽度值，在热核宽度最优时，利用代数连通度以内的拉普拉斯矩阵特征向量重构冲击信号分量;最后对重构冲击信号进行包络解调分析，诊断滚动轴承故障类型。滚动轴承数值仿真故障信号和应用实例表明，与文献方法及谱峭度方法比较，该方法能更有效地提取轴承故障信号中的冲击特征，实现背景噪声干扰情况下的轴承故障的有效诊断。
1 GFT与高斯函数加权路图
路图仅考虑时间序列中相邻时间点的连接关系，与时间序列的连续发展特点一致。因此，本文将振动信号转换为路图。不同加权方式下，图拉普拉斯矩阵的特征向量代替了传统三角函数基函数，为基于图傅里叶变换的特殊信号分量提取提供了可能。本节首先介绍图信号处理的基础概念，基于图给出GFT的定义;再结合典型路图拉普拉斯矩阵特征分解的特点，解释可以用于冲击信号提取的原因;最后对高斯加权热核宽度参数对冲击信号重构的影响作了分析。
1.1 图与图傅里叶变换
图由顶点通过顶点间连边构成，一般的，无向、加权图可表示为G=（V，E，A），其中V表示图顶点的有限集合（顶点数量|V|=N），E为连边的集合，A表示邻接矩阵，其中元素aij为顶点vi与vj间连边eij的加权值。特别地，对于无向图，A为对称矩阵。假定x ∈ RN为定义在图G上的信号，Xn为图G中第，n个顶点的信号值。则顶点间权值通过顶点间距离dij=xi-xj进行构造，本文采用高斯函数加权
图傅里叶变换（GFT）定义为图信号关于拉普拉斯矩阵特征向量函数的展开，其主要思想是将图信号正交投影于拉普拉斯矩阵的特征向量上，从而实现图信号的“频率”分析，将图信号研究从顶点域拓展到图谱域。
按式（2）构建拉普拉斯矩阵并进行正交分解
Lul=λlul，l=0，1，2…，N-1（3）
式中 l为阶次;λl为第l阶特征值;ul为第l阶特征向量。记特征向量矩阵U=[u0u1…uN-1;特征值对角矩阵A=diag Lλ0λ1…λN-I]。
结合拉普拉斯矩阵的结构可知其特征分解有以下特点：
（1）L为对称半正定矩阵，因此，特征向量矩阵为单位正交矩阵，UTU=i.
（2）L的最小特征值为0，其数量为k（k≥1），将其按升序排列，0=λ0…=λk-1<λk≤…≤λn，其中λk为第一个大于零的特征值，称其为代数连通度（Algebraic Connectivity）。
（3）零特征值的数量k等于图中连通分量的个数。连通分量定义为图的一个子图，且子图中任意两个顶点连通。
基于拉普拉斯矩阵及其特征向量可定义GFT，与传统傅里叶变换定义类似，将图拉普拉斯矩阵特
1.2 路图拉普拉斯矩阵特征分解与冲击特征
路图是时间序列最为吻合的构图对象，仅相邻时间点间存在边连接。考虑如图1所示包含9个顶点的简单路图，对应时间序列为
x={8，1，1，8，1，1，8，1，1} （7）
该时间序列可看作冲击信号。理想情况下，设定距离较远的顶点间不连接且权值为0，否则为1.则仅顶点2与3，5与6，8与9存在边连接，且权值aij=1.显然，图1中存在6个连通分量，分别为{{v1}，{v2，v3}，{v4}，{v5，v6}，{v7}，{v8，v9}}
。
按式（2）计算拉普拉斯矩阵并进行标准正交分解。特征值对角矩阵如下式所示
A=diag[000000222]（8）
可见零特征值的数量与图1中连通分量数目相同，代数连通度λk=2此外，图1中顶点v1，v4，v7，呈现的冲击特征，分别对应第1，3，5个孤立连通成分，在特征向量矩阵中体现为特征向量u0，u2，u4为单位脉冲向量，且冲击发生时刻与顶点位置一致，特征向量的冲击特性如图2所示。因此，原始时间序列中的冲击特征可以由特征向量u0，u2，u4进行路图傅里叶逆变换重构。对于实际信号，由于存在噪声与其他分量信号干扰，考虑利用零特征值对应的所有特征向量重构信号中的冲击分量，如下式所示
1.3 热核宽度参数对冲击特征重构结果的影响
对实际信号，由于无法直接判定顶点间的连接关系，故采用式（1）所示的加权方式衡量顶点间连接关系与权重，热核宽度值δ将影响冲击特征的重构结果。因此，需要分析热核宽度参数对冲击特征重构结果的影响。以图3所示仿真冲击信号为例，冲击特征频率为50Hz，定义代数连通度能量占比为路图傅里叶变换中代数连通度对应阶次以内能量与所有阶次能量的比值
比较在不同信噪比下（SNR=-5，0，5），零特征值数目尾、代数连通度能量占比r及重构信号包络谱在特征频率处的幅值随热核宽度（0.001≤δ≤0.06）变化的情况，结果如图4所示。
從图4（a）中可看出，随着热核宽度增加，零特征值数目逐渐减少，因此冲击特征也逐渐集中于较低阶次的有限特征向量中;从图4（b）中可看出，随着热核宽度增加，代数连通度能量占比逐渐减小;图4（c）中特征频率处幅值存在最大值，此时热核宽度为最优热核宽度值。综合比较，热核宽度最优时，既要保证零特征值的数目较小又要保证代数连通度能量占比较大;相同冲击信号在不同信噪比下的最优热核宽度值不同，且信噪比越小取值越大，对应幅值也随噪声增加而减小。因此，针对具体信号最优热核宽度值是一个变动的值，不同信号取值不同，本文中以重构信号包络谱在特征频率处的幅值为适应度函数，结合粒子群优化算法对热核宽度参数进行寻优，以达到最优的冲击信号重构效果。
2 基于最优加权路图GFT的滚动轴承故障诊断
滚动轴承存在局部故障时会产生周期性的瞬态冲击特征，通常使用包络解调方法即可诊断出轴承故障。然而，对于滚动轴承的早期故障，振动信号中包含轴承故障信息的冲击成分能量很小，往往淹没在机械系统谐波振动与背景噪声中，此时直接使用包络解调方法难以取得好的效果。因此，需要对包含轴承故障信息的冲击成分进行提取。
本文提出基于最优加权路图GFT冲击提取的滚动轴承故障诊断方法，该方法将振动信号转换为路图，在高斯加权下使用有限拉普拉斯矩阵特征向量重构冲击信号;以重构冲击信号故障特征处的包络谱值为适应度函数，利用粒子群算法寻优最优高斯热核宽度参数;对最优冲击分量进行包络解调诊断轴承故障。故障诊断流程如图5所示，步骤如下：
（1）采集滚动轴承故障振动信号，采样点数为N，并将其转换为N顶点路图。
（2）使用粒子群优化算法对高斯热核宽度参数进行寻优。首先确定热核宽度值的寻优范围。本文中选取0.001≤δ≤δmax，δmax为信号方差;
（3）初始化粒子群。设置粒子群参数维度、粒子群种群大小及最大迭代次数;
（4）计算各个粒子的适应度。首先构建路图对应的图拉普拉斯矩阵;然后对拉普拉斯矩阵进行特征分解，通过代数连通度λk确定用于冲击信号重构的特征向量个数k;最后，根据式（9）进行信号重构，并计算其故障特征频率处对应的包络谱值作为粒子的适应度值;
（5）计算粒子群中的个体最优值和全局最优值;
（6）更新粒子的速度和位置，达到最大迭代次数时停止更新，返回最优热核宽度参数值;
（7）根据最优热核宽度δ，重新按步骤（4）重构冲击分量，对重构信号进行包络解调分析，并诊断滚动轴承故障类型。
3 仿真分析
为了验证本文提出方法的有效性，利用滚动轴承数值仿真模型，模拟滚动轴承局部故障，数值仿真模型表示如下
式中h（t）为冲击信号分量;T为仿真故障冲击的时间间隔，则故障特征频率f=1/T;M为仿真冲击的个数;Bm为第m个冲击的幅值;β为阻尼比，fre表示共振频率;u（t）为单位阶跃函数。s（t）为两种频率叠加的谐波信号，频率分别为f1，f2。x（t）为仿真合成信号，n（t）为添加的随机噪声。本文中，设置故障特征频率f=50Hz，采样频率fs=4096Hz，采样点数取N=4096，添加噪声后合成信号的信噪比为-10dB，其他参数设置如表1所示，仿真冲击信号，谐波信号及合成信号的时域波形如图6所示。
从图6（a）中可以明显看出周期瞬态冲击成分及发生时刻，叠加上周期干扰分量及随机噪声后，图6（c）合成信号中已无法识别出任何瞬态冲击成分。对仿真合成信号直接作包络谱分析，结果如图7所示，受谐波分量的干扰，包络解调分析方法在对信号进行解调时将两种余弦频率之差作为解调结果解出，因此图中出现了频率峰值220Hz，而故障特征频率及其倍频处无明显峰值，说明轴承故障冲击特征已完全淹没。
应用本文提出的最优加权路图GFT方法提取冲击信号并进行解调分析。按照图5中的算法流程，设置热核宽度最大寻优值为信号方差δmax=var（x）=0.75，则寻优范围为[0.01，0.75]。设置粒子种群大小为10，最大迭代次数为50，适应度随进化代数曲线如图8（a）所示，在第19次迭代运算时取得最大值，其后保持稳定。此时计算得到的最优热核宽度值为δ*=0.279，使用最优热核宽度值提取冲击信号并做包络解调分析，包络谱图如图8（b）所示。从图中可以能够清晰地识别仿真信号的故障特征频率f以及其2，3倍频，以上结果说明了该方法提取冲击特征的可行性。
作为对比，分别使用文献提出的欧氏距离加权下的冲击特征信号提取方法和谱峭度（SK）方法提取冲击信号并做包络谱分析，结果分别如图9和10所示。从图9中虽然可以看出1倍及2倍故障特征频率，但均存在54和103Hz的干扰谱线;此外，其故障特征频率处幅值也明显低于本文提出方法。使用SK方法计算快速Kurtogram图，谱峭度指标最大时的中心频率为fc=1184Hz，层阶数为6.5.从图10可以看出，冲击信号包络谱仅有1倍故障频率，且存在其他较大的干扰频率。以上结果表明，最优加权路图GFT方法能有效重构周期瞬态冲击信号。
为了更为客观地对比不同方法提取滚动轴承故障信号中冲击特征的性能，引入包络谱故障特征比率指标Rf对本文提出方法、欧氏距离加权GFT方法和SK方法提取的冲击信号进行评价。该指标表示为
式中 5为信号Hilbert包络谱;S（f），S（2f）和S（3f）分别为故障特征频率及其2，3倍频谱线的幅值。该指标直观地反映了故障特征频率在包络谱中的幅值比率。使用三种方法提取的冲击特征信号的Rf值如表2所示，易知本文提出方法提取的沖击特征信号对应的Rf值取得最大值，因此更为明确地诊断了故障。
4 应用实例
为验证基于最优加权路图GFT方法提取冲击信号的有效性，分别设置滚动轴承内、外圈故障进行实验。实验在轴承振动模拟实验台进行，实验装置如图11所示。内、外圈故障均采用激光在6307E型滚动轴承的内圈和外圈上切割宽为0.15mm、深为0.13mm的槽，以模拟内圈和外圈故障。滚动轴承内圈故障实验时的转轴转速为1200r/min，外圈故障实验时的转轴转速为1500r/min.采用加速度传感器采集轴承振动信号，采样频率设置为4096Hz，采样点数取1000.通过计算得到内圈故障特征频率为98.8Hz，外圈故障特征频率为76.5Hz。
图12（a）为滚动轴承内圈故障振动信号的时域波形图，虽然能够看到冲击特征，但冲击信号的周期性不明显。直接对时域信号作包络解调分析，包络谱如图12（b）所示，可以看出，故障特征频率谱线fi被其他与轴承内圈故障无关的干扰谱线淹没，因此不能识别故障类型。
使用本文提出方法对内圈故障振动信号中的冲击成分进行重构。首先设置热核宽度最大寻优值为信号方差δmax=var（xi）=1107，则寻优范围取[0.01，1107];然后设置粒子群大小为10，最大迭代次数为50，如图13（a）所示适应度值在第6次迭代运算后取得最大值，之后保持稳定，此时计算得到的最优热核宽度值为δ*=519.3;最后使用最优热核宽度值重新提取冲击信号并做包络解调分析，包络谱图如图13（b）所示。从图中可以明顯看出1倍和2倍轴承内圈故障特征频率谱线，准确诊断出了内圈故障。
同样作为对比，使用基于欧氏距离加权的冲击特征提取方法提取冲击信号并作包络解调分析，结果如图14所示。从图14中可知工倍及2倍故障特征频率谱线并不明显，且存在较多的干扰谱线，故障特征频率及其倍频处的幅值也均小于本文提出方法得到的结果。
使用SK方法对内圈故障信号进行滤波后并作包络谱，首先计算快速Kurtogram图，谱峭度指标最大时的中心频率为fc=1792Hz，层阶数为3.5.滤波信号包络谱如图15所示，易知图中存在2倍故障频率谱线，故障频率处无明显峰值。
滚动轴承外圈的故障特征频率为fo=76.5Hz，振动信号的时域波形如图16（a）所示，从图中看不到明显冲击特征。对外圈故障时域信号作解调分析，包络谱如图16（b）所示，图中外圈故障特征频率谱线f。被周围的64Hz等干扰谱线淹没，因此不能够直接诊断出轴承故障类型。
使用本文提出方法对外圈故障振动信号中的冲击成分进行提取。首先设置热核宽度最大寻优值为信号方差δmax=var（xo）=16.62，则范围取[0.01，16.62];然后初始化粒子群寻优算法，设置粒子群大小为10，最大迭代次数为50，如图17（a）所示，适应度值在第35次进化后取得最大值，其后全局最优值保持稳定，此时计算得到的最优热核宽度值δ*=4.8;最后使用最优热核宽度值重新提取冲击信号并做包络解调分析，包络谱图如图17（b）所示。从图中可以明显看出轴承外圈故障特征频率及其倍频谱线，准确诊断出了外圈故障。
作为对比，使用基于欧氏距离的冲击特征信号提取方法和SK方法提取冲击信号并做包络谱分析，结果分别如图18和19所示。从图18中可以看出，故障特征频率及其倍频谱线f。和2f。并不明显。使用SK方法，谱峭度指标最大时的中心频率为fo=341Hz，层阶数为3.5.滤波信号包络谱如图19所示，图中仅存在1倍故障频率谱线。从上述结果可知，针对滚动轴承故障的冲击特征，本文提出方法的提取结果更为有效。
为了进一步对比三种不同方法提取内、外圈故障信号中冲击特征的效果，使用Rf指标进行对比分析，结果如表3所示。从表中可知，对于内、外圈故障，本文提出方法提取的冲击特征信号的Rf均取得最大值，说明了该方法提取冲击特征的有效性。
5 结论
（1）本文提出了基于最优加权路图GFT的冲击分量提取方法，将图信号处理方法运用于冲击信号的提取中，并结合包络解调实现了轴承的故障诊断。该方法首先将振动信号转换为路图，在最优高斯加权下，利用有限图拉普拉斯矩阵特征向量结合GFT方法重构冲击特征信号，实现了滚动轴承故障的有效诊断。
（2）实验结果表明，本文提出方法能够在背景噪声或其他强故障特征分量干扰下，有效重构故障冲击特征，并准确诊断滚动轴承故障类型。与欧氏距离加权下的重构结果和谱峭度方法重构结果对比，验证了最优加权路图GFT冲击提取方法的准确性和优越性。
