基于图滤波理论的齿轮箱故障特征提取应用研究

Application Research of Gearbox Fault Feature Extraction Based on Graph Filtering Theory

摘要: 风力发电机、汽车、轮船等机械设备都离不开齿轮箱。传统的故障特征提取方法有时难以正确检测故障。本文采用图信号处理方法，提出了基于图滤波理论的齿轮箱故障特征提取应用研究。首先，将齿轮箱的振动信号转化为路图，建立路图拉普拉斯矩阵(Laplacian Matrix)，得到图信号的特征值和特征向量；其次，进行图傅里叶变换(Graph Fourier Transform, GFT)，选择图常数滤波器和图理想带通滤波器对图傅里叶变换后的信号进行重构；最后，将重构信号再次进行图傅里叶变换，从而提取齿轮箱故障特征。研究结果表明：经过图信号滤波的信号相比没有经过图信号滤波的信号故障特征提取更明显，具有一定的故障诊断效果。

Abstract: Wind turbines, cars, ships and other mechanical equipment are inseparable from the gear box. Traditional fault feature extraction methods are difficult to detect faults correctly sometimes. Using graph signal processing method, this paper puts forward the application research of gearbox fault feature extraction based on graph filtering theory. Firstly, the vibration signal of gearbox is trans-formed into road diagram signal, where its Laplacian Matrix is established, then, extracting its ei-genvalues and eigenvectors; Secondly, the Graph Fourier transform (GFT) is established, the graph constant filter and the graph ideal band-pass filter are selected to filter the graph signal after Fou-rier transform, and then, the graph signal on convolution is reconstructed; Finally, the reconstruct-ed signal is transformed into graph Fourier transform, which extracts the fault features of gearbox from different angles. The results show that the signal filtered by graph signal is better than that without graph signal filtering, and the fault feature extraction is more obvious.

1. 引言

齿轮箱有许多优点，如传动平稳、传动比大、便于控制等，广泛应用于机械设备上。但其承受载荷大，振动幅度大，很容易发生故障。传统的齿轮箱检测方法主要靠听、看和通过大量的实验，如拆解查看、空运转试验、传动噪声检测、密封性和渗漏检查等 [1] 。齿轮箱故障特征提取的方法也有很多。如：温度法，通过检测齿轮箱轴承座的温度来判断齿轮箱是否正常；油样分析法，通过分析齿轮箱内润滑油的金属颗粒大小和形状来分析齿轮的状态等。目前，齿轮箱故障特征提取的方法主要是振动法，即采集并分析齿轮箱的振动信号来监测故障状态。振动信号分析方法最常用的是时域分析和频域分析。时域分析法对冲击信号、周期信号等有明显效果，但对于齿轮箱复杂信号效果并不理想；频域分析法则是基于傅里叶变换原理，此方法只适用于线性平稳的振动信号，然而，齿轮箱振动信号具有非平稳特性。其他常用的齿轮故障特征提取方法有：Wigner-Ville分布、主成分分析、奇异值分解循环平稳性分析和经验模式分解等。常用的统计指标有均方根、峰度、峰因子等。本文将运用图信号滤波理论，将齿轮箱振动的时域信号转化为图信号，进而进行滤波分析，提出基于图滤波理论的齿轮箱故障特征提取应用研究。

近年来，基于时域和频域滤波理论广泛运用在图论中，使得图信号处理的应用越来越广泛。2012年，Arun Venkitaraman, Saikat Chatterjee，Peter Händel提出了图信号的Hilbert变换 [2] ；同年，K. Narang提出了双通道小波滤波器组，用于分析定义在任意有限加权无向图顶点上的函数 [3] ；2013年，Kenji Nakahira和Atsushi Miyamoto提出了多链路小波，可以精确地模拟数据 [4] ；同年Aliaksei Sandryhaila和José M. F. Moura提出了图信号滤波器，还将傅里叶变换运用到了图信号处理技术中，证明了其与图邻接矩阵的广义特征向量的关系，并研究了它们的性质 [5] ；2015年，David I Shuman和Christoph Wiesmeyr等人提出频谱自适应图小波和倒频谱框架，对于顶点具有更好的辨别能力 [6] 。2017年韩墨研究了图信号的采样和重构，并且将采样、重构后的图信号进行不同种类的滤波 [7] 。目前关于图滤波理论在齿轮箱故障诊断中的应用相对较少。本文运用图信号滤波和重构理论，基于图变换的方法，提取齿轮箱的故障特征，此方法比一般的时域分析和频域分析范围更广，尤其是对于非平稳的振动信号有较好的效果。

2. 图信号及图信号处理理论

图是由几个顶点以及连接顶点的边组合而成，图的顶点由V表示，边由E表示，图的集合记作 $G=\left\{V,{\rm E}\right\}$ ，顶点的集合记作 $V=\left\{v_1,v_2,v_3,\cdots ,v_N\right\}$ ，边的集合则表示成 $E=\left\{e_1,e_2,e_3,\cdots ,e_N\right\}$ ，其中， $v_i$ 表示图的第i个顶点， $e_i$ 表示图的第i条边 [2] 。在图中，每个顶点 $v_i$ 都对应着一个信号值 $f\left(i\right)$ ，图顶点集合 $V=\left\{v_1,v_2,v_3,\cdots ,v_N\right\}$ 对应的图信号值： $f=\left\{f\left(1\right),f\left(2\right),f\left(3\right),\cdots ,f\left(N\right)\right\}$ ，该图信号值可以表达图的顶点信号 [8] 。对于连接顶点的边有无权重，可以将图分为有权图和无权图，所谓无权图，就是所有边的权重都为1的图，即

$W_{ij}=1$ (1)

对于有权图，有三种常用的计算权重的方式，分别是：

$W_{ij}={\text{e}}^{-\frac{\Vert f\left(i\right)-f\left(j\right)\Vert }{2t^2}}$ (2)

$W_{ij}=f\left(i\right)-f\left(j\right)$ (3)

$W_{ij}=1-\frac{\langle f\left(i\right),f\left(j\right)\rangle }{\Vert f\left(i\right)\Vert \Vert f\left(j\right)\Vert }$ (4)

式中热核函数加权如式(2)，其中t表示热核的宽度， $f\left(i\right)$ 和 $f\left(j\right)$ 表示图顶点的值。欧式距离加权如式(3)，这种计算权重的方式能在一定程度上体现出顶点与顶点之间的差别。余弦距离加权舅式(4)。

路图为相邻顶点由一条无向边连接，不相邻顶点之间没有边组成，如图1所示。

Figure 1. Path graph signal

图1. 路图信号

离散的振动信号可以由路图表示，其中路图的顶点顺序即为离散信号点的时间采样点的顺序，顶点值即为离散信号的幅值。此路图顶点的集合为 $V=\left\{v_1,v_2,v_3,v_4,v_5,v_6\right\}$ ，边的集合为 $E=\left\{\left(v_1,v_2\right),\left(v_2,v_3\right),\left(v_3,v_4\right),\left(v_4,v_5\right),\left(v_5,v_6\right)\right\}$ 。图的属性可以用矩阵表达，对于图中顶点的度，用度对角矩阵来表示，记作D。度对角矩阵的对角线记作 $d_{ii}$ ，表示顶点的度，即 $\text{diag}\left(v_i\right)$ ，度对角矩阵见式(5)：

$D_{n\times n}=\{\begin{array}{l}\text{diag}\left(v_i\right),i=j\\ 0,i\ne j\end{array}$ (5)

度对角矩阵可表示图顶点之间的连接关系， $d_{ii}$ 也可以表示这个点 $v_i$ 在图中的重要程度，对于路图的度对角矩阵，见式(6)：

$D_{n\times n}=\{\begin{array}{l}1,i=j=1,i=j=n\\ 2,i=j\\ 0,i\ne j\end{array}$ (6)

邻接矩阵可以表示图中各顶点之间的距离，记作W，在矩阵中，第i行，第j列的元素记为 $w_{ij}$ 。当顶点 $v_i$ 和 $v_j$ 有边连接时，在本文中用欧式距离加权，若两个顶点之间没有直接连接，则记 $w_{ij}=0$ 。

图的拉普拉斯矩阵记作L。L包含图中顶点和边的信息，既有顶点连接边的数目，又有连接两顶点见的权重，因此可以很好的阐述图的性质，此矩阵的定义为度对角矩阵D与邻接矩阵W之差，L矩阵定义见式(7)：

$L=D-W$ (7)

由公式(7)可知，L是半正定矩阵，L的特征向量记为U， $U=\left[u_1,u_2,\cdots ,u_N\right]$ ，且

$U^{\text{T}}U=I$ (8)

特征值记为 $\Lambda$ ， $\Lambda =\text{diag}\left[{\lambda }_1,{\lambda }_2,{\lambda }_3,\cdots ,{\lambda }_N\right]$ ，且 $\Lambda$ 中所有的元素 ${\Lambda }_{ij}\ge 0$ 。

另有对L进行标准化，记作 $\tilde{L}$ ，将L标准化后，其特征值 $\lambda$ 的取值范围为[0,1]， $\tilde{L}$ 的定义为：

$\tilde{L}=I-D^{-\frac{1}{2}}W{D}^{\frac{1}{2}}$ (9)

图傅里叶变换(Graph Fourier transform, GFT)可以把顶点域的图信号转化到图谱域上，所谓顶点域；而图谱域则表示图的特征值与图幅值的关系 [9] 。由公式(7)可知，图的拉普拉斯矩阵L的特征向量为正交矩阵，而L的每一个特征值 ${\lambda }_i$ 都对应着一个特征向量 $u_1$ ，所以GFT可以由特征向量推导，且特征向量 $u_1$ 为GFT的基函数，与傅里叶变换的基函数不同的是GFT的基函数是离散且变化的，因此GFT的定义为：

$F\left({\lambda }_i\right)={\displaystyle \underset{k=0}{\overset{N-1}{\sum }}f\left(k\right)u_i^{\ast }\left(k\right)}$ (10)

其中，i为图中顶点的顺序，k为离散信号和特征向量中元素的顺序。同理，由傅里叶逆变换或GFT可推得图傅里叶逆变换(Graph inverse Fourier transform, GIFT)：

$f\left(k\right)={\displaystyle \underset{i=0}{\overset{N-1}{\sum }}{u}_i\left(k\right)F\left({\lambda }_i\right)}$ (11)

3. 图滤波理论及信号重构方法

3.1. 图滤波理论

类似与经典信号滤波器，如低通、带通、高通滤波器和高斯滤波器。在图滤波理论中，对应相应的滤波器，有：图谱域滤波和顶点域滤波。

图谱域滤波：在振动信号的滤波中，针对频域的滤波公式为：

$F_{out}\left(\xi \right)=F_{in}\left(\xi \right)H\left(\xi \right)$ (12)

其中 $H\left(\xi \right)$ 为滤波器，根据频域的滤波公式，即可推得图谱域的滤波公式：

$F_{out}\left({\lambda }_i\right)=F_{in}\left({\lambda }_i\right)H\left({\lambda }_i\right)$ (13)

其中 。 $H\left({\lambda }_i\right)=\text{diag}\left(H\left({\lambda }_0\right),H\left({\lambda }_1\right),\cdots ,H\left({\lambda }_{N-1}\right)\right)$ 。

顶点域滤波：将经典信号处理的频域滤波进行GFT，即可得到时域滤波。对于式(13)进行GIFT，得：

$f_{out}\left(i\right)={\displaystyle \underset{l=0}{\overset{N-1}{\sum }}{F}_{in}\left({\lambda }_l\right)H\left({\lambda }_l\right)u_l\left(i\right)}$ (14)

上式即为图顶点域的滤波公式，但是由于符号不同，公式如下：

$f_{out}\left(k\right)={\displaystyle \underset{i=0}{\overset{N-1}{\sum }}{F}_{in}\left({\lambda }_i\right)H\left({\lambda }_i\right)u_{{\lambda }_i}\left(k\right)}$ (15)

在本文中，采用顶点域滤波为图信号滤波公式。

由公式(15)可知，需要建立一个滤波器 $H\left({\lambda }_i\right)$ ，在图信号处理中，也有许多滤波器的选择：

常量图滤波器

$H=CI$ (16)

其中，C为常数，可根据GFT图像进行调整，I为单位矩阵。此滤波器相当于是一个常数。

理想图带通滤波器

当L进行标准化后 ${\lambda }_i$ 大大减小，而在图信号处理中， $\lambda$ 相当于振动信号处理的 $\omega$ ，路图顶点的k相当于振动信号处理的时间t，因此在图信号处理中的低通滤波器为

$H\left(\lambda \right)=\{\begin{array}{l}0,{\lambda }_a\ge 0\\ {\lambda }_n,{\lambda }_a\le {\lambda }_n\le {\lambda }_b\\ 0,{\lambda }_b\le 2\end{array}$ (17)

其中， ${\lambda }_a$ 和 ${\lambda }_b$ 为截止频率，理想图带通滤波器可以在各个频段任意设置，过滤任意的频段，与之类似的，还有图低通滤波器和图高通滤波器。

除了以上几个图滤波器，还有图热核滤波器、图高斯滤波器等，滤波器的选择需要根据信号的类型来判断，有时为了较好的分析图信号，需要用到不同的滤波器进行滤波。

3.2. 图信号的重构

图信号的重构有许多方法，如卷积重构、插值重构等。下面介绍几个简单的重构方法：

1) 最邻近插值：

最邻近插值法是最简单的插值方法，即对于时间t，在相邻的离散信号中取点，设相邻离散信号的点为 $x_1$ 和 $x_2$ ，则对应的信号值为 $f\left(x_1\right)$ 和 $f\left(x_2\right)$ ，在 $x_1$ 和 $x_2$ 中取点 $x_{1.5}$ ，则

$f\left(x_{1.5}\right)=f\left(x_1\right)$ (18)

这两种此方法要求 $x_1$ 和 $x_2$ 距离较近，若是采样周期过长，则不适合用这两种方法。

2) 卷积重构

卷积的物理意义为一组数据乘以其相应权重的和，卷积的定义公式为

$\left(f\ast g\right)\left(x\right)={\displaystyle {\int }_{-\infty }^{+\infty }f\left(t\right)g\left(x-t\right)\text{d}t}$ (19)

可以把插值离散信号看成是离散信号与某一函数的卷积，这个函数称为核函数。对于线性插值，其核函数可表达为；

$w\left(x\right)=\{\begin{array}{l}1-x,\left|x\right|<1\\ 0,\left|x\right|\ge 1\end{array}$ (20)

将其看成是离散信号，且取值频率是1时，只有在 $x=0$ 处 $w\left(x\right)=1$ ，其他部分 $w\left(x\right)=0$ 。

4. 齿轮箱故障特征提取方法

图滤波理论的应用非常广泛，常应用于二维图像中，如盒式滤波器、梯度滤波器和二项滤波器等能一定程度地去除图像噪点，使图像像素之间的过度更加地平滑自然。本文研究针对于单个齿轮的故障特征提取。在齿轮箱运行过程中，采集离散的振动信号，将其转化为图信号，经过图信号滤波，然后重构图信号，最后分别对其进行FT、GFT，以提取并齿轮箱故障特征。基于齿轮箱的图滤波故障诊断方法如下：

1) 振动信号采样：为了避免信号混叠，必须用大于原信号两倍的频率对原信号进行采样，本课题的采样频率远大于原信号的频率；

2) 图信号建立：由于经典的信号处理对于复杂、非平稳的振动信号效果不佳，因此采用图信号处理，将振动信号转化为路图信号，并建立相应的度矩阵D和邻接矩阵W，由此得出拉普拉斯矩阵L，求其特征值和特征向量；

3) 图信号滤波：本文采用常数图滤波器和图带通滤波器，对路图信号进行GFT，以便将图信号和滤波器进行顶点域滤波，能对于转换后的图信号降噪，或是将图信号重新排布，减小噪音的干扰；

4) 信号重构：进行滤波后的图信号为了可以更精准地还原出原始信号，需要进行信号重构，根据式(19)的核函数与滤波后的信号进行信号卷积，以实现信号的插值重构；

5) 提取故障特征：对重构后的信号进行GFT、FT和特征值提取，并建立其相应的图谱，以便观察；

6) 分析故障：将无故障的齿轮和轴承与有故障的齿轮和轴承的图谱进行相应的对比，分析器故障特征。

提取齿轮箱故障特征流程图如图2所示：

Figure 2. Flow diagram of Gearbox fault feature extraction

图2. 齿轮箱故障特征提取流程图

5. 应用实例

本文研究的实验数据来自湖南大学实验数据中心，选取1200转的齿轮，采样周期 $T=0.00012\text{s}$ ，因

此采样频率 $f=\frac{1}{T}=8333.3\text{Hz}$ ，取采样个数为1024个点，实验数据有正常的齿轮和人为制造的断齿，其

振动信号时域波形如图3所示：

Figure 3. Time domain waveform of gear vibration signal

图3. 齿轮振动信号时域波形

由图3可以看出，在1200转的断齿和正常的时域波形有细微差别，分别用图滤波理论对齿轮箱中的齿轮和轴承的故障进行故障特征提取，用GFT图、FT图和振动信号的时域图像来对滤波后的效果进行描述。

5.1. 基于常数滤波器的齿轮箱故障特征提取

常数滤波器提取故障特征提取的流程图如图4所示：

Figure 4. Flow chart of fault feature extraction using constant filter

图4. 使用常数滤波器提取故障特征的流程图

仍取上述的实验数据进行对比分析，由式(16)的常数滤波器，选取常数 $C=-0.8$ ，建立矩阵L后，进行常数滤波器的顶点域滤波，并根据式(19)对其进行信号重构，由于选取点数为1024，为了使得重构后的信号更加集中，因此取 $x\in \left[-512,511\right]$ ，最后将重构后的信号再次作GFT，以方便观察故障的位置。1200转正常齿轮和断齿的图信号处理结果如图5所示。

由图5可知，正常齿轮的振动频率呈均匀分布，且幅值较小，而断齿的振动频率分布并不均匀，且在特征值 $\lambda =2$ 左边的幅值有一个突变，且相对于特征值为2的右边来说，左边的幅值较大，达到了0.6，而在正常齿轮的幅值中，幅值基本上为0.16左右，相较于图3，可以非常明显地看出两者的区别所在，便于后续的故障诊断。

图常数滤波器可以更改GFT图特征值的范围以及幅值的大小，当常数 $\Vert C\Vert$ 越大，其图像的能量越分

Figure 5. 1200 rpm gear diagram signal comparison

图5. 1200转齿轮图信号对比

散，幅值越大；而常数 $\Vert C\Vert$ 越小则其能量更加集中，但是幅值减小。采用有波动断齿的实验数据，采样频率 $f=8192\text{Hz}$ ，选取其中任意1500个点作为采样长度，分别取常数为±1，±50和±0.1，对比图如图6所示。

由图6可知，不同的常数对GFT图的相对形状没什么太大的影响，但是对总体形状有较大改变，如图(b)所示，当 $\Vert C\Vert$ 变大时，其GFT图像的特征值范围变大，且幅值对比于图(a)来说有巨大变化，但是故障的特征值更为分散；(c)图的情况与(b)图正好相反，其特征值范围减小，能量更加聚集，更容易看出故障的所在，但其幅值随之减小。当 $\Vert C\Vert =1$ 时，其与GFT图几乎相同。

用常数滤波器后的GFT图对比FT图也有较大优势，依然用上述1200转齿轮数据，采样频率，任意采取其中1500个采样点作为采样长度，齿轮断齿和正常的FT图如图7所示。对比图6和图7可知，在常数滤波器滤波后的情况下，同样是断齿信号，常数为50，其效果较FT更好，因为其幅值较FT图更大，若是比较相同的幅值，当常数滤波器的常数为0.1时，其能量更加集中，相较于FT更加明显。

(a) (b) (c)

Figure 6. Comparison of constant values of constant filters

图6. 常数滤波器的常数值比较

Figure 7. Comparison of FT between normal gear and broken tooth

图7. 正常齿轮和断齿的FT比较图

5.2. 基于图带通滤波器的齿轮箱故障特征提取

由前述可知，故障大多发生在中高频部分，而常数滤波器并不能将一些频段的信号过滤，只能起到调整特征值和幅值的作用，且由图3可知，当矩阵L标准化后，其能量较为分散，方便用带通滤波器进行滤波，而低通滤波器和高通滤波器可能会把故障滤去，或是效果不明显，采用的方法为用带通滤波器对图信号进行滤波后并重构，将图信号转化为振动信号，并将其与没有滤波的振动信号作对比。

理想带通滤波器提取故障特征的操作软件为Matlab，其流程图如图8所示。

选取另一组1200转齿轮断齿的故障，采样频率 $f=8192\text{Hz}$ ，采样个数为1024个，将其故障的GFT图与经过高通滤波器的GFT图进行对比，先出GFT图，后根据其故障的所在区域，以确定带通滤波器的特征值选取范围，高通滤波器的特征值选取范围为 $0.9\le \lambda \le 0.99$ ，实验结果如图9所示。

由图9可知，带通滤波器将图像的低频和高频部分过滤，而留下了中频部分，在 $\lambda <0.8$ 的部分幅值减弱了许多，并且在 $\lambda >2.2$ 部分减少了数量，但是其整体幅值同样下降。

Figure 8. Flow chart of fault feature extraction for idealized bandpass filter

图8. 理想图带通滤波器的故障特征提取流程图

Figure 9. Comparison diagram of gear broken tooth with or without bandpass filter

图9. 齿轮断齿有无带通滤波器对比图

接下来将其转化到振动信号与原始信号作比较，如图10：

Figure 10. Comparison diagram of gear broken tooth with or without bandpass filter

图10. 齿轮断齿带通滤波与原始振动信号的比较

由图10可知，经过图带通滤波后的振动信号的故障特征要比原始信号明显许多，且在故障处能量依然较为密集，幅值也没什么变化，具有较高的优势。同时，想要信号更加分散，也可以利用图常数滤波器将其能量分散开来，可以更好地判断其故障特征值所在，图像如图11所示：

Figure 11. Vibration diagram of bandpass filter with constant filter

图11. 加了常数滤波器的带通滤波器振动图

由图11可知，在时间为0.1s右边的振动信号相比于单纯的带通滤波器后的振动信号，其幅值大幅度上升，且振动的频率不变，可以更好地判断出故障所在位置。

经过这两个图滤波器后的FT图同样有明显效果，选取另一组1200转断齿，将滤波前的FT图与滤波后的FT图作对比，如图12所示。

由图12可明显看出再频率为2500 Hz至3500 Hz的波形被过滤，且故障特征的峰值不变，滤波后的图像比滤波前更加明显。

Figure 12. Comparison of FT images before and after filtering

图12. 滤波前后FT图像的对比

6. 结论

1) 常数滤波器可以将特征值进行缩放，并改变其幅值，可以较为方便地提取故障特征；

2) 图理想带通滤波器可以对特征值进行筛选，从而较为精确地提取出故障特征值，但是特征值需要人为调节；

3) 图常数滤波器可以和图理想带通滤波器一起使用，以得到更好的效果；

4) 当故障频率非常小时，用GFT可能不如FT更加直观，因此在滤波后，可以两种方法同时运用，以便互相对比。

基金项目

湖南省教育厅一般项目(19C1473)，谱图理论在机械故障诊断中的应用研究；怀化学院一般项目(HHUY2019-26)，图滤波理论在机械关键零部件故障特征提取中的应用。