SuperGraph: Spatial-temporal graph-based feature extraction for rotating machinery diagnosis

引用和代码链接

C. Yang, K. Zhou and J. Liu, “SuperGraph: Spatial-Temporal Graph-Based Feature Extraction for Rotating Machinery Diagnosis,” in IEEE Transactions on Industrial Electronics, vol. 69, no. 4, pp. 4167-4176, April 2022, doi: 10.1109/TIE.2021.3075871.
代码链接

摘要

振动信号本身是不规则的且往往是含有噪声的。这使得频谱分析难以提取高阶特征。随着图论被应用于频谱分析中来提高特征提取的性能。本文提出一种基于时空图的旋转机械故障特征提取方法超图法。具体说就是基于图论的频谱分析被用于构建时空图。然后，从构建的时空图中提取基于拉普拉斯矩阵的特征向量。通过这种方式，将时空图转换为一维向量以进一步构造超图，其中超图的每个节点表示时空图，并且超图由许多局部图组成。只有相同类型的节点才能完全连接形成一张局部图。因此，图分类的任务可以转化为对超图中的节点进行分类。在建立图卷积网络学习和获取深度特征后，根据SoftMax模型识别节点标签。通过将原始数据转换为图形得到隐藏的结构和拓扑信息的方法不仅滤除区分了振动信号中的噪声，也提高了特征提取的质量。在两个基准数据集和一个实际的实验平台上进行实验，来验证所提出故障诊断方法的有效性。

一、Introduction

1. 机械故障大背景介绍
2. 常用谱分析介绍。说明基于傅里叶变换的频谱分析的短板，引出短时傅里叶变换（STFT）。STFT可以通过采用窗口功能来处理震动信号中的噪声，对非平稳信号进行操作。通常，采集的振动信号通常带有噪声和不规则性，这给使用STFT提取高级特征带来了困难。为了提高基于频谱分析的特征提取质量，本文将图论引入频谱分析。
3. 介绍现有图模型，并引出本文中提出的一种结合谱分析和图论的时空图模型。
4. 将图卷积神经网络（GCN）引入机器故障检测。

贡献：

1） 原始信号被转换成时空图来构造拉普拉斯矩阵。
2） 采用基于拉普拉斯矩阵的特征提取方法，将图形转化为特征向量，将图形分类转化为节点分类。
3） 为了提高图数据较大时的诊断性能，开发了一种超图。
4） 将GCN模型应用于轴承故障诊断，在小样本复杂信号故障识别中取得了较好的效果。

二、时空图模型

A． 短时周期图提取

这里是对原始数据进行STFT，再计算信号的周期图得到信号的功率谱估值
给定的振动信号表示为 { x ( n ) } , n ∈ [ 0 , N − 1 ] \{x(n)\},n∈[0,N-1] {x(n)},n∈[0,N−1]，具有N个观察值。在离散时间域和频域中的STFT计算是：

$△M$表示时间跳变，$△f$表示频率的跳变，$Y(k,m)$是频率$k△f$和时间$m△M$下的输出，$\ast$号表示共轭，合适的窗口$\omega$，也称为加权函数，可以根据处理信号的特性选择合适的窗口。本文选择的窗函数为$Hanning$，它是最常用的窗函数之一。
然后通过以下公式计算短时周期图：

其中T是窗口的长度。类似的，$p(k,m)$是频率指数$1\le k\le K$和时间指数$1\le m\le M$时的输出。

B． 时空图构造

短时周期图是一种方便的指示器，包括时域和频域信息，用于描述机器的运行状态。需要用合适的模型评估具有高动态特性的条件。因此，时空图模型被用来解决上述问题。与传统的特征提取方法相比，时空图不仅包含了时频域信息，还包含了图的结构信息，进一步提高了提取特征的判别能力。执行以下步骤来构造时空图。
1）图形表示
图表示为 G = { V , E , W } G=\{V, E, W\} G={V,E,W}。$V$表示图的节点，$E$表示这些节点之间的边缘连接，W是图的权重矩阵。
2）权重矩阵构造
在本文中$W_{ij}$由以下公式计算得到

$Dis()$表示某一时间点t上，点$F_i$和点$F_j$之间的欧式距离。欧式距离是一个常用的距离函数。
3）图形构造
如第二、A节所述，STFT根据振动信号运行，其频率分为$K$个部分。以短时周期图中的每个频率$F_1$，$F_2$，…，$F_k$为节点构建的图，以及构建时空图的过程如图1所示。每对没有属性的节点相互连接，构成时空图。在本文中，时空域中的节点数为33。

章节小结和代码

作者对一个故障样本信号划分了33个窗口，并对每个窗口进行傅里叶变换（短时傅里叶变换的本质就是找到最合适的窗口大小，对窗口内的信号进行傅里叶变换）。然后通过周期图计算33个窗口各自的功率谱估值。并将33个窗口视为33个节点，节点之间相互连接，节点属性为各自窗口的功率谱估值。最后计算节点与节点之间的欧式距离作为权重矩阵。

def calpkm(data):
    data = np.array(data)  # 将数据转化为数组格式
    freqs, times, Sxx = signal.stft(data, fs=25000, window='hanning', nperseg=64, noverlap=0)  
    # 短时傅里叶变换（数据，采样频率，窗口函数，窗口长度，窗口重叠数
    # f频率，t时间，Sxx：STFT时频数据
    pkm = np.zeros((len(Sxx),len(Sxx[0]),len(Sxx[0][0])))   # pkm,短时周期图计算
    pkmn = np.zeros((len(Sxx),len(Sxx[0]),len(Sxx[0][0])))
    for k in range(len(Sxx)):       # k为频率指数
        for i in range(len(Sxx[0])):   # i为时间指数
            for j in range(len(Sxx[0][0])):   # j为幅值
                pkm[k][i][j] = abs(Sxx[k][i][j])**2/64   # 计算公式，周期T为64，abs（）为绝对值
    return pkm

def weight(matrix):
    W = np.zeros((len(matrix),len(matrix)))   # 创建一个n*n的矩阵
    for i in range(len(matrix)):   #
        for j in range(i,len(matrix)):
            for k in range(len(matrix[0])):
                W[i][j] += np.sqrt(np.square(matrix[i][k] - matrix[j][k]))   # 求两点间的欧式距离
                # np.square为求括号中的平方，np.sqrt为求括号中的平方根
            W[j][i] = W[i][j]   # 对称阵（构建的时空图为无向图）
    return W

三、基于图论的特征提取

A． 图论分析

利用图中隐藏的结构和拓扑信息的优点，提出了一种基于图论的特征提取方法。图论主要涉及两类矩阵，即邻接矩阵和拉普拉斯矩阵，本文只考虑拉普拉斯矩阵。为了得到拉普拉斯矩阵，在第二节中首先定义了权重矩阵。然后，将拉普拉斯矩阵的正交分解得到的特征值作为模式识别模型的输入。
在谱图分析中，图拉普拉斯算子是一个关键运算，其定义如式（4）所示。为了获得拉普拉斯矩阵，度矩阵的定义如式（5）所示。

其中$D$是对角矩阵，$L$是拉普拉斯矩阵，$d_{i,i}$是对角矩阵对应位置的元素。

B． 特征向量提取

拉普拉斯矩阵是一个实对称矩阵，可以通过如下正交分解操作
$L=U\Lambda U^T$ (6)
其中$\Lambda =diag({\lambda }_0,{\lambda }_1,\dots ,{\lambda }_{n-1})$是特征值矩阵，$U=[u_0,u_1,\dots ,u_{n-1}]$是特征向量。
每个时空图都可以用一维特征向量$F$表示，它由特征值组成，表示为$F=[{\lambda }_0,{\lambda }_1,\dots ,{\lambda }_{n-1}]$。

章节小结和代码

主要就通过拉变换将上一章得到的$W$矩阵转变为了一维特征向量$F$。文章最精髓的地方。（没有这一步直接进入图卷积的话，会走一个图分类，而做了这一步之后直接转变成了节点分类。）

def calDifLaplacian(W1):
    n = len(W1)
    D = np.zeros((n,n))    # D为对角矩阵，大小和W1相同
    I = np.identity(n)     # 大小为n的单位阵
    for i in range(n):
        D[i][i] = sum(W1[i])  # D矩阵的对角元素为W当行元素之和
    L = D - W1            # 获得拉普拉斯矩阵
    return L

四、图卷积网络

构造的图卷积网络有两个主要步骤：输入和卷积。将由时空图组成的超图视为输入。为了获得超图，最重要的操作是确定每个图之间的相关性。本节介绍如何获取超图并构造GCN。

A． 超图构造

如第二节所述，超图也可以表示为 G S = { V S , E S , A } G_S = \{V_S,E_S,A\} GS​={VS​,ES​,A}。$V_S$表示超图的节点集，其数目为$M$，其中$M$为总样本数；$E_S$表示超图的边级，其中表示节点的相同类型信号以完全连接的无相图形式连接。$A\in R^{m\times m}$表示任意两个节点之间连接的相邻矩阵，相邻矩阵中的元素的值为0或1。相邻矩阵的定义如式（7）所示。

图2显示了三种类型的节点和三种类型的边连接的示例。可以看出，超图由三个局部图组成，这三个局部图在空间上是独立的。对于标记的节点，颜色表示故障类型，相同颜色的节点相互连接，并且不同颜色的节点没有边缘连接。对于两个未标记的样本，它们的数据点是在同一时间段内收集的，因此它们属于相同的故障类型，但我们不知道具体的故障。然后将相同故障类型的未标记节点互连以形成局部图。当GCN聚集相邻特征时，信息通过边缘传播。在超级图中，边上的两个节点属于相同的故障类型，因此聚合的特征可以包含有关相应故障的更多信息。

在图谱分析中，对称归一化拉普拉斯矩阵被广泛用于图卷积运算，其定义式如式（8）所示。图拉普拉斯矩阵$L_S$也可以在$L_s=U_s{\Lambda }_s{U}_s^T$中操作。

其中$I_n$是单位矩阵。

B． 图卷积网络

图卷积层对图像信号的处理可分为两部分，第一部分是图卷积，第二部分是特征变换。
1）图卷积
由于传统的卷积运算不能直接应用于图，因此[25]中开发了一种基于图拉普拉斯算子的泛化卷积方法，如式（9）所示。根据式（6）和式（8），图卷积可以表示为式（10）。

其中$X\in R^{m\times s}$表示输入矩阵，$g_{\theta }()=diag()$是傅里叶域中的$\theta \in R^n$参数化的滤波器，$g_{\theta }({\Lambda }_S)$可以看作是特征值的函数。$U_s\in R^{m\times m}$表示特征向量矩阵，$Y\in R^{m\times s}$是图卷积的输出。
然而，在普通图卷积中计算$g_{\theta }({\Lambda }_S)$是一个困难且计算代价昂贵的过程[26]。为了避免这个问题，在[27]中引入了切比雪夫多项式展开式。它表明，$g_{\theta }({\Lambda }_S)$可以很好地用切比雪夫多项式$T_k(x)$到K阶的展开式来逼近。

其中${\theta }_k$是切比雪夫系数的向量，$T_K(x)$是一个递归计算式，他的计算式如式（12）。${\tilde{\Lambda }}_S$ 是${\Lambda }_S$的归一化版本，可根据式（13）获得。

其中${\lambda }_{max}$是${\Lambda }_s$的最大元素，$\tilde{\Lambda }$中的元素在-1到1的范围内。
基于信号$x$与滤波器$g_{\theta }^{\prime }$的卷积定义，图卷积可通过公式（14）计算。

2）特征变换
图卷积运算可以看作是在不改变特征维数的情况下对信号进行滤波。将图形信号乘以可训练权重矩阵后，可改变特征尺寸，如等式（15）所示。

其中$W^{(i)}\in R^{m\times j}$是第i张图卷积层的可训练权矩阵，$Cheb()$表示切比雪夫卷积函数，$X^{\prime }\in R^{m\times j}$是第$i$个图形卷积的最终输出。

章节小结和代码

上一章最后，本文得到了每个故障样本的节点属性（对W矩阵做了拉变换之后得到的特征向量）。这一章节看论文的代码会有一个问题，就是本章节一开始说的节点连接部分（邻接矩阵的构建）同类样本之间相互连接。这一块光看论文的代码会感觉有问题。其实，他选取了每种故障类型的30个样本，然而30个样本是分为3段取出的，这三段时间取出的10个样本理论上来说他们属于一个类别（同一时间一个机械一般不会出现两种故障吧，我觉得。）所以，论文代码将这10个连续采集的样本之间产生一条无向边。接下来的图神经网络就是对谱卷积的一个优化的网络（切比雪夫图卷积）。对于图卷积这部分，涉及内容还是挺多的，光看这篇论文是不太够的，还是需要系统的学习一下，所以论文中所有的图卷积部分没有着重介绍，请谅解。

edge_index = [[], []]  # edge connection
for i in range(0, len(a), 10):
    for j in range(i, i+10):
        for k in range(i, i+10):
            if j != k:
                edge_index[0].append(j) 
                edge_index[1].append(k)
# edge_index[0]中有一条边的起始节点，edge_index[1]中存放的是一条边的结束节点，
# edge_index[0]的大小代表边数。

五、基于图卷积网络的故障诊断

A． GCN模型构建

构建的GCN模型，如图3所示，可在等式（16）中简单描述。根据Zhang等人[22]讨论的结果，本文选择了三层切比雪夫卷积，GCN模型还包括两层激活函数校正线性单元$（ReLU）$和一层$SoftMax$。

其中$\sigma$是$ReLU$函数； $W^{(0)}\in R^{C\times H}$是具有$C$个特征的输入层到具有$H$个特征的隐藏层1的权重矩阵；$W^{(1)}\in R^{H\times I}$是具有$H$个特征的隐藏层1到具有$I$个特征的隐藏层2的权重矩阵；$W^{(2)}\in R^{I\times J}$是具有I个特征的隐藏层2到具有$J$个特征的隐藏层3的权重矩阵；最后通过隐藏层2输出$F$个特征。

B． GCN模型的输入

通过构造时空图和提取特征向量，将振动信号转换为以特征向量为携带信息的节点。然后，根据每对节点之间的相关性构造超图。GCN模型的输入是一个表示超图的矩阵。图4示出了如何获得输入矩阵。

C． GCN模型的输入

在最后一层，将最后一层卷积输出的学习特征输入到$SoftMax$层进行故障诊断。从一般数据域到图形域的整个故障诊断过程如图5所示。

章节小结和代码

对模型的一个总写。

# 定义网络
class Net(torch.nn.Module):   # GCN model
    def __init__(self):  # 定义隐藏层和卷积核k
        super(Net, self).__init__()
        self.conv1 = ChebConv(33, 30, K=3)
        self.conv2 = ChebConv(30, 25, K=3) 
        self.conv3 = ChebConv(25, 11, K=3)

    def forward(self):    # 传播网络
        x, edge_index, edge_weight = data.x, data.edge_index, data.edge_attr
        x = F.relu(self.conv1(x, edge_index, edge_weight))
        x = F.relu(self.conv2(x, edge_index, edge_weight))
        x = self.conv3(x, edge_index, edge_weight)
        return F.log_softmax(x, dim=1)