基于辛几何模式分解与图结构增强动态时间规整的轴承故障诊断方法

本发明属于机械设备故障诊断领域，特别是指一种基于辛几何模式分解与图结构增强动态时间规整的轴承故障诊断方法。

背景技术：

近年来国内制造业发展势头迅猛。技术不断地更新迭代，机械设备不断向着连续化、高速化、重载化以及智能化的方向发展，这导致整个机械系统越来越复杂，各个部件之间的关联性更加紧密。当机械系统中某一个零部件发生故障时，往往会影响整个系统的正常运行。而旋转机械常作为整个机械系统的动力源，保证旋转机械的正常运作就显得尤为重要。其中，轴承作为旋转机械的核心零部件，也是易损零部件，对轴承的工作状态进行实时监测，以便早期发现异常故障是十分有价值和现实意义的。

多基于傅里叶变换来实现，傅里叶变换对于平稳信号具有很好的稳定性，而对于经过强噪声干扰导致信号产生扭曲的情况处理效果不佳。而现代信号处理方法在机械故障诊断领域迫切需要解决的问题之一就是有效提取非平稳信号特征的轴承故障特征提取方法。

辛几何模式分解(symplecticgeometrymodedecomposition,sgmd)，本质上是一种信号分解方法，它通过信号构建轨迹矩阵，然后用轨迹矩阵构建哈密顿矩阵，最后对其分解并对角平均，依据成分的相似性合并最后构成各个独立分量，从而实现主要成分与噪声分离。图结构方法(基于图形变化点检测方法)是构造正常信号的邻接矩阵作为图结构，对图结构分解获得标准特征向量。未知信号用同样步骤构造图结构，并用标准特征向量求解未知信号的图结构，获得矩阵。求矩阵的非对角部分的f范数作为指标，当超过设定的阈值则判断为有故障。然而，图结构方法易受输入信号的影响，当输入信号非平稳时，会影响最后f范数指标的表现，直接影响到故障诊断的成功和准确与否。而且，图结构方法无法区分故障类型，不能判断轴承是哪个部分发生了故障。动态时间规整(dynamictimewarping,dtw)，本质上就是用来计算两个时间序列之间的欧氏距离来衡量两者之间相似度的方法。它通过找寻两个时间序列之间的最优匹配后再计算它们之间的欧式距离，是一个计算不等长时间序列、找寻时间序列之间相似性的方法。将该方法用于计算不同故障类型信号之间的相似度，根据相似度来判断故障类型。但是，该方法同样容易收到强噪声的影响，动态时间规整方法很难区分非平稳信号的类型。

为此，本发明人为了克服动态时间规整易受强噪声影响的不足，对信号采用了两步特征增强，提出了一种基于辛几何模式分解与图结构增强动态时间规整的轴承故障诊断方法，有关这方面研究，目前尚无报道。

技术实现要素：

本发明的目的是为了克服现有技术存在的缺点和不足，而提供一种基于辛几何模式分解与图结构增强动态时间规整的轴承故障诊断方法。

为实现上述目的，本发明的技术方案是包括以下步骤：

①对n个正常轴承诊断样本信号做辛几何模式分解，依据均方根选取出能量最大的分量，构建其图结构xadji(i＝1，2，...，n)；

②根据已建立的图结构计算平均图结构，以克服采样过程中相位变化造成的缺陷，分解图结构，获得特征值和特征向量；

③当采集到一个新的时域信号时，提取一个周期的段，用辛几何模式分解提取能量最大的分量，就可以建立相应的图结构xadj，new；

④用所述的特征向量γ来分解xadj，new，获得新的矩阵g1，最后再提取矩阵g1中的非对角矩阵，求其f范数f1；

⑤从已知标签的各类数据中分别随机选取30个周期段求f范数作为各类模板样本；

⑥当采集到一个新的时域信号时，随机提取其中的30个周期段，计算30个周期段的f范数作为测试样本；

⑦最后，用动态时间规整来计算测试样本和各类模板样本之间的相似度距离，最终通过测试样本与各类模板样本之间的相似度距离来判断故障类型。

进一步设置是所述的步骤①中，辛几何模式分解的步骤：

(1)相空间重构

设振动信号s＝(x1，x2，...，xn)，然后重构轨迹矩阵x：

设m＝n-(d-1)τ，其中d为嵌入维数，τ为延迟时间，选择适当的d和τ值以获得相应的轨迹矩阵x，计算振动信号s的功率谱，然后在psd中获得最大峰值的频率fmax，如果归一化频率小于给定的阈值10^-3，d设为n/3，其中n为数据长度，否则嵌入维数设为d＝1.2×(fs/fmax)；

(2)辛正交矩阵的qr分解

为了构造哈密顿矩阵m，应用轨迹矩阵的自相关来获得协方差对称矩阵a＝x^tx，因此m定义为:

设n＝m²，因此n和m都是哈密顿矩阵，辛正交矩阵构造为q，即:

其中b为上三角矩阵，bij＝0(i＞j+1)；b可以通过对矩阵n应用施密特正交化得到，它们的特征值分别是λ1，λ2，...，λd，设σi是矩阵a的特征值，qi是对应σi的特征向量，设其中z是重构轨迹矩阵，则初始单分量轨迹矩阵z表示为：

z＝z1+z2+…+zd(4)

(3)对角平均变换

上一步骤得到的初始单分量矩阵z是大小为m×d的矩阵，所以需要对初始的单个分量重新排序，重构的矩阵z可以通过对角平均变换成一组新的长度为n的时间序列，就可以得到长度为n的d个新的时间序列，d个新的时间序列之和就是原来的时间序列，

对于任一初始单分量矩阵zi，将矩阵的元素定义为zij，其中1≤i≤d，1≤j≤m，且d^*＝min(m，d)，m^*＝max(m，d)，n＝m+(d-1)τ，如果m＜d，设否则，因此，对角平均传递矩阵如下所示:

根据等式(5)，重构矩阵z转化为一系列新的长度为d×n的矩阵yi(y1，y2，...，yn)，将原始时间序列分解成趋势和频带不同的d个独立叠加分量，d个单分量信号通过对角平均获得:

y＝y1+y2+…+yd(6)

然后依据均方根来选取能量最大的分量，均方根rms(rootmeansquare)表达式为：

设最终选出的分量为yrms，构造其图结构，即邻接矩阵，计算时间序列中的每两个点之间的欧氏距离作为两个点之间的权重，最终得到图结构xadj。

进一步设置是所述的步骤②具体为：

根据步骤①构造n个正常信号的图结构xadji(i＝1，2，...，n)，求其平均获得平均图结构分解平均图结构获得特征值y^t和特征向量γ，

其数学表达式为：

进一步设置是所述的步骤③具体为：

采集信号时，随机设置起点，提取一个周期的段nnew，采样长度l可由l＝fs/r确定，其中fs为采样频率，r为转速，随后，xnew经过步骤①得到图结构xadj，new。

进一步设置是所述的步骤④具体为：

用特征向量г来分解xadj，new，得到新的矩阵g1，其数学表达式为：

xadj，new＝гg1г(9)

计算矩阵g1中的非对角矩阵的f范数，其数学表达式为：

xadj，new＝γg1γ＝γ(diag[g1])γ+γ(non-diag[g1])γ(10)

f＝||(non-diag[g1])||f(11)。

进一步设置是所述的步骤⑤具体为：

从已知标签的各类数据中分别随机选取30个周期段，经过步骤③和步骤④分别获得30个f范数，将这30个f范数组成一个特征向量f＝[f1，f2，...，f30]，轴承的典型故障有：内圈故障、外圈故障以及滚动体故障，加上正常信号，这样就可以获得4类信号的模板样本，分别为fnormal、fin、fout、fball。

进一步设置是所述的步骤⑥具体为：

当采集到一个新的时域信号时，随机提取其中的30个周期段，计算30个周期段的f范数作为测试样本，采集一段未知的信号作为测试样本，经过步骤③和步骤④，得到ftest特征向量作为测试样本，

进一步设置是所述的步骤⑦具体为：

动态时间规整方法的原理描述如下：对于两个序列c＝c1，c2，...，cm和q＝q1，q2，...，qm，相应元素之间的距离可以通过距离函数计算为d(ci，qi)，从而形成一个n×m的距离矩阵，在传统的dtw方法中，距离函数是平方欧氏距离，平方欧氏距离的数学表达为：

d(ci，qi)＝(ci-qi)²(12)

然后，通过使累计距离最小化，得到最优匹配路径u＝(u1，u2，...，cl)，其中max(m，n)≤l≤m+n-1，路径应该满足如下所示的局部约束：

a.端点约束：弯曲路径的起点和终点应该恰好是距离矩阵的第一个和最后一个点，也就是u1＝(c1，q1)，ul＝(cm，qn)，

b.连续性约束：一个时刻，路径可以前进一步，也就是当uk＝(ci，qj)，uk+1＝(ci+1，qj+1)，然后ci+1-ci≤1，qi+1-qi≤1，

c.单调性：当uk＝(ci，qj)，uk+1＝(ci+1，qj+1)，则ci≤ci+1，qi≤qi+1，最后，动态时间规整定义为：

用动态时间规整来分别计算测试样本ftest和各类模板样本fnormal、fin、fout、fball之间的相似度距离,最终，在故障距离图中可以直观地找到与测试样本ftest相似度距离最近的模板样本，则可以直接判断测试样本是什么故障类型。

本发明的优点是：

本发明方法采取了两步特征增强，加强动态时间规整对故障类型的区分能力。一方面利用辛几何模式分解剔除噪声，有效地提高信噪比，使图结构的输入更加稳定，更有效的提取特征，另一方面计算多个f范数组成特征向量，增强稳定性，提高最终动态时间规整的鲁棒性，从而可以可靠地检测出轴承故障类型。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，根据这些附图获得其他的附图仍属于本发明的范畴。

图1本发明的轴承故障诊断方法流程图；

图2图结构xadj示意图；

图3正常(a)、内圈故障(b)、外圈故障(c)、滚动体故障(d)的原始信号时域波形图；

图4实施例中，未知健康状况轴承的原始信号时域波形图；

图5本发明实施例的故障距离图；

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述。

如图1所示，本发明采用技术方案包括：

①对n个正常样本做辛几何模式分解，依据均方根选取出能量最大的分量，构建其图结构xadji(i＝1，2，...，n)。

②根据已建立的图结构计算平均图结构，以克服采样过程中相位变化造成的缺陷。分解图结构，获得特征值和特征向量。

③当采集到一个新的时域信号时，提取一个周期的段，用辛几何模式分解提取能量最大的分量，就可以建立相应的图结构xadj，new。

④用前面获得的特征向量γ来分解xadj，new，获得新的矩阵g1，最后再提取矩阵g1中的非对角矩阵，求其f范数f1。

⑤从已知标签的各类数据中分别随机选取30个周期段求f范数作为各类模板样本。

⑥当采集到一个新的时域信号时，随机提取其中的30个周期段，计算30个周期段的f范数作为测试样本。

⑦最后，用动态时间规整来计算测试样本和各类模板样本之间的相似度距离。最终通过测试样本与各类模板样本之间的相似度距离来判断故障类型。

上述各步骤具体地包括以下步骤：

①对n个正常样本做辛几何模式分解，依据均方根选取出能量最大的分量，构建其图结构xadji(i＝1，2，...，n)。

辛几何模式分解的步骤：

(1)相空间重构

设振动信号s＝(x1，x2，...，xn)，然后重构轨迹矩阵x：

设m＝n-(d-1)τ，其中d为嵌入维数，τ为延迟时间。选择适当的d和τ值以获得相应的轨迹矩阵x。计算振动信号s的功率谱(powerspectraldensity,psd)，然后在psd中获得最大峰值的频率fmax。如果归一化频率小于给定的阈值10^-3，d设为n/3，其中n为数据长度。否则嵌入维数设为d＝1.2×(fs/fmax)。

(2)辛正交矩阵的qr分解

为了构造哈密顿矩阵m，应用轨迹矩阵的自相关来获得协方差对称矩阵a＝x^tx，因此m可以定义为:

设n＝m²，因此n和m都是哈密顿矩阵。辛正交矩阵构造为q，即:

其中b为上三角矩阵，bij＝0(i＞j+1)。b可以通过对矩阵n应用施密特正交化得到，它们的特征值分别是λ1，λ2，...，λd。设σi是矩阵a的特征值，qi是对应σi的特征向量。设zi＝qisi，其中z是重构轨迹矩阵，则初始单分量轨迹矩阵z表示为：

z＝z1+z2+…+zd(4)(3)对角平均变换

上一步骤得到的初始单分量矩阵z是大小为m×d的矩阵，所以需要对初始的单个分量重新排序，重构的矩阵z可以通过对角平均变换成一组新的长度为n的时间序列，就可以得到长度为n的d个新的时间序列，d个新的时间序列之和就是原来的时间序列。

对于任一初始单分量矩阵zi，将矩阵的元素定义为zij，其中1≤i≤d，1≤j≤m，且d^*＝min(m，d)，m^*＝max(m，d)，n＝m+(d-1)τ。如果m＜d，设否则，因此，对角平均传递矩阵如下所示:

根据等式(5)，重构矩阵z转化为一系列新的长度为d×n的矩阵yi(y1，y2，...，yn)。将原始时间序列分解成趋势和频带不同的d个独立叠加分量。d个单分量信号通过对角平均获得:

y＝y1+y2+…+yd(6)

然后依据均方根来选取能量最大的分量。均方根rms(rootmeansquare)表达式为：

设最终选出的分量为yrms，构造其图结构(邻接矩阵)，计算时间序列中的每两个点之间的欧氏距离作为两个点之间的权重，最终得到如图2所示的图结构xadj。

②、根据已建立的图结构计算平均图结构，以克服采样过程中相位变化造成的缺陷。分解图结构，获得特征值和特征向量。

根据步骤1构造n个正常信号的图结构xadji(i＝1，2，...，n)，求其平均获得平均图结构

分解平均图结构获得特征值y^t和特征向量г。其数学表达式为：

③、当采集到一个新的时域信号时，提取一个周期的段，用辛几何模式分解提取能量最大的分量，就可以建立相应的图结构xadj，new。

采集信号时，随机设置起点，提取一个周期的段xnew(采样长度l可由l＝fs/r确定，其中fs为采样频率，r为转速)。随后，xnew经过步骤1可以得到图结构xadj，new

④、用前面获得的特征向量γ来分解xadj，new，获得新的矩阵g1，最后再提取矩阵g1中的非对角矩阵，求其f范数f1。

用特征向量γ来分解xadj，new，得到新的矩阵g1，其数学表达式为：

xadk，new＝γg1γ(9)

计算矩阵g1中的非对角矩阵的f范数，其数学表达式为：

xadj，new＝γg1γ＝γ(diag[g1])γ+γ(non-diag[g1])γ(10)

f＝||(non-diag[g1])||f(11)

⑤、从已知标签的各类数据中分别随机选取30个周期段求f范数作为各类模板样本。

从已知标签的各类数据中分别随机选取30个周期段，经过步骤3和步骤4可以分别获得30个f范数，将这30个f范数组成一个特征向量f＝[f1，f2，...，f30]。轴承的典型故障有：内圈故障、外圈故障以及滚动体故障。加上正常信号，这样就可以获得4类信号的模板样本，分别为fnormal、fin、fout、fball。

⑥、当采集到一个新的时域信号时，随机提取其中的30个周期段，计算30个周期段的f范数作为测试样本。

采集一段未知的信号作为测试样本，经过步骤3和步骤4，得到ftest特征向量作为测试样本。

⑦、最后，用动态时间规整来计算测试样本和各类模板样本之间的相似度距离。最终通过测试样本与各类模板样本之间的相似度距离来判断故障类型。

动态时间规整方法的原理描述如下：对于两个序列c＝c1，c2，...，cm和q＝q1，q2，...，qm，相应元素之间的距离可以通过距离函数计算为d(ci，qi)，从而形成一个n×m的距离矩阵。在传统的dtw方法中，距离函数是平方欧氏距离。平方欧氏距离的数学表达为：

d(ci，qi)＝(ci-qi)²(12)

然后，通过使累计距离最小化，得到最优匹配路径u＝(u1，u2，...，cl)，其中max(m，n)≤l≤m+n-1。路径应该满足如下所示的局部约束：

a.端点约束：弯曲路径的起点和终点应该恰好是距离矩阵的第一个和最后一个点，也就是u1＝(c1，q1)，ul＝(cm，qn)。

b.连续性约束：一个时刻，路径可以前进一步。也就是当uk＝(ci，qj)，uk+1＝(ci+1，qj+1)，然后ci+1-ci≤1，qi+1-qi≤1。

c.单调性：当uk＝(ci，qj)，uk+1＝(ci+1，qj+1)，则ci≤ci+1，qi≤qi+1。

最后，动态时间规整定义为：

用动态时间规整来分别计算测试样本ftest和各类模板样本fnormal、fin、fout、fball之间的相似度距离。最终，在故障距离图中可以直观地找到与测试样本ftest相似度距离最近的模板样本，则可以直接判断测试样本是什么故障类型。

实施案例1：轴承内圈故障诊断

取某真实传动系统的信号，实验用的轴承，型号为er-12k,其中滚动体个数为8个，滚动体直径为7.9375mm，轴承节圆直径d＝33.4772mm,接触角α＝0,实验轴承空载运行，转速1800r/min,fs＝25.6khz。正常轴承采集一部分数据，带有内圈故障的轴承采集一部分数据，带有外圈故障的轴承采集一部分数据，带有滚动体故障的轴承采集一部分数据，将这些数据作为按照提出的方法处理获得正常、内圈故障、外圈故障、滚动体故障的特征向量模板。

正常、内圈故障、外圈故障、滚动体故障的原始信号波形如图3所示。采集一段未知健康状况轴承的信号，其时域波形如图4所示。将其按照所提方法的步骤处理，得到待测特征向量。通过动态时间规整求待测特征向量与各个特征向量模板的相似度距离，并可在故障距离图中直观地判断故障类型。故障距离图如图5所示，最终判断故障为内圈，将轴承拆下，经过对比发现，确实是内圈故障。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，所述的存储介质，如rom/ram、磁盘、光盘等。

以上所揭露的仅为本发明较佳实施例而已，当然不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。