CN112683391A - 一种eemd和gzc机器状态监测方法及装置 - Google Patents

一种eemd和gzc机器状态监测方法及装置 Download PDF

Info

Publication number
CN112683391A
CN112683391A CN202011238756.1A CN202011238756A CN112683391A CN 112683391 A CN112683391 A CN 112683391A CN 202011238756 A CN202011238756 A CN 202011238756A CN 112683391 A CN112683391 A CN 112683391A
Authority
CN
China
Prior art keywords
data
signal
eemd
fractal
curve
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Withdrawn
Application number
CN202011238756.1A
Other languages
English (en)
Inventor
豆春玲
寇兴磊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong Kerishen Intelligent Technology Co ltd
Original Assignee
Shandong Kerishen Intelligent Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shandong Kerishen Intelligent Technology Co ltd filed Critical Shandong Kerishen Intelligent Technology Co ltd
Priority to CN202011238756.1A priority Critical patent/CN112683391A/zh
Publication of CN112683391A publication Critical patent/CN112683391A/zh
Withdrawn legal-status Critical Current

Links

Images

Landscapes

  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种EEMD和GZC机器状态监测方法及装置，利用EEMD算法对设备振动信号进行分解，利用非线性判别算法去除分解噪声分量和趋势项，保留分形信号分量，采用Newton插值函数对极值点进行插值，利用最小二乘法拟合包络，分离频率调制部分，利用GZC算法估计瞬时频率并计算相应的瞬时尺度，根据分析尺度确定振动信号去趋势结果，计算去趋势信号的多重分形谱，提取多重分形谱的左端点、右端点和极值点坐标作为设备运行状态的特征参数，识别设备运行状态，将上述算法部署到设备状态监测装置，能够准确区分设备运行状态，设备状态监测装置具有良好的柔性和便携性，便于工程应用。

Description

一种EEMD和GZC机器状态监测方法及装置
技术领域
本发明涉及设备状态监测与故障诊断领域，具体涉及一种EEMD和GZC机器状态监测方法及装置。
背景技术
设备振动信号包含丰富的分形特征，这些分形特征能够描述设备的运行状态。盒维数、功率谱分析和重标极差方法可以估计平稳信号的单重分形参数，去趋势波动分析（DFA）能够估计非平稳信号的单重分形维数。然而，设备出现故障时，其振动信号通常是非平稳的，且具有多重分形特征，这时传统的分形维数估计方法会产生比较大的误差。多重分形去趋势波动分析（MFDFA）能够估计非平稳信号的多重分形参数，但是MFDFA方法存在着分析尺度需要人工确定、拟合多项式趋势阶数难以确定和数据段之间不连续的问题。目前，已经有文献提出了基于EMD的MFDFA版本（MFDFAemd），用来解决MFDFA存在的问题。然而，MFDFAemd采用的线性滤波方法容易破坏原始信号的分形结构，且存在着负频率现象，这些缺陷严重影响了MFDFAemd的应用效果。综上所述，现有技术难以准确提取设备振动信号的多重分形特征，难以准确检测设备运行状态。
发明内容
本发明要解决的问题是针对以上不足，提出一种EEMD和GZC机器状态监测方法及装置（本发明提出的方法简称为MFDFAoeemd）。采用本发明所提出的方法对设备振动信号进行分析，能够有效提取设备振动信号的多重分形特征，克服MFDFA方法存在的分析尺度需要人工确定、拟合多项式趋势阶数难以确定和数据段之间不连续的问题，解决MFDFAemd方法存在的原始信号分形结构破坏和负频率现象，具有分析结果准确度和精确度高，设备运行状态识别结果正确率高等优点。
为解决以上技术问题，本发明采取的技术方案如下： 一种EEMD和GZC机器状态监测方法，其特征在于：包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用集合经验模式分解（Ensemble Empirical Mode Decomposition,EEMD）算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure 217877DEST_PATH_IMAGE001
，其中，ci（k）代表由EEMD算法得到的第i个分量，rn（k）代表由EEMD算法得到的趋势项，本例中，n=10；
步骤3：采用非线性判别算法从EEMD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量cf（k）, f=1,2,…,p，p代表滤波后剩余分量的数量；
步骤4：确定cf（k）的局部极大值和局部极小值，采用Newton插值函数分别对cf（k）的局部极大值和局部极小值进行插值，利用最小二乘法分别拟合cf（k）的上包络u(k)和下包络l(k)，则cf（k）的包络定义为
Figure 603859DEST_PATH_IMAGE002
，符号|x|表示对x取绝对值；
步骤5：重复执行公式
Figure 852437DEST_PATH_IMAGE003
m次，j=1,2,…,m，直到
Figure 982067DEST_PATH_IMAGE004
，
得到cf（k）的频率调制部分FMm(k)，ej(k)代表cj(k)的包络，cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
步骤6：采用广义零交叉法（Generalized zero-crossing,GZC）计算FMm(k)的平均周期，得到cf（k）的瞬时尺度sf；
步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为
Figure 100002_DEST_PATH_IMAGE005
；
步骤8：将Y s (k)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
步骤9：计算每段数据的方差：
Figure 999702DEST_PATH_IMAGE007
Figure 789541DEST_PATH_IMAGE009
；
步骤10：计算q阶函数：
Figure 790995DEST_PATH_IMAGE011
；
步骤11：改变s的取值，s=sf，f=1,2,…,p，重复上述步骤3到步骤10，得到关于q和s的方差函数Fq(s);
步骤12：如果x(k)存在分形特征，则F q (s)和尺度s之间存在幂律关系：F q (s)~s H(q)，H(q)代表x(k)的广义Hurst指数；
当q=0时， H(0)通过下式所定义的对数平均过程来确定：
Figure 724316DEST_PATH_IMAGE012
；
步骤13：计算信号x(k)的标准标度指数τ(q)=qH(q)-1，本例中，q在(-5, 5)范围内取值;
步骤14：计算信号x(k)的奇异指数α和多重分形谱f(α):
α=H(q)+q H’(q),
f(α)=q(α-H(q))+1，其中H’(q)代表H(q)的一阶导数；
步骤15：提取多重分形谱f(α)的左端点、右端点和极值点所对应的奇异指数，利用这3个参数来描述设备的运行状态；
步骤16：将上述步骤所述方法部署在状态监测装置上，对设备状态进行监测。
进一步的，所述步骤2中EEMD算法包括以下步骤：
1）向数据x0(k)添加白噪声序列产生一个新数据xj(k) ：
Figure 862036DEST_PATH_IMAGE014
std[x0(k)]代表数据x0（k）的标准差，wnj(k)代表wnj中的第k个数据，wnj代表第j个随机产生的白噪声序列，wnj幅值为1，1≤j≤K；x0(k)代表权利要求1所述步骤2中x(k)；本例中，K=100；
2）对xj(k)执行经验模式分解，得到n个分量和一个趋势项
Figure 324242DEST_PATH_IMAGE015
cij(k)代表对xj(k)执行经验模式分解得到的第i个分量，rnj(k)代表对xj(k)执行经验模式分解得到的趋势项；
3）计算K次分解结果的平均值
Figure 812992DEST_PATH_IMAGE016
ci(k)表示对x0(k)进行集合经验模式分解得到的第i个分量，rn(k)表示对x0(k)进行集合经验模式分解得到的趋势项。
进一步的，所述步骤3非线性判别算法包括以下步骤：
1) 对信号c（k）执行重排操作和替代操作，经重排操作得到的数据用cshuf（k）表示，替代操作后得到数据用csurr（k）表示；
2) 对c（k）、cshuf（k）和csurr（k）分别执行多重分形去趋势波动分析（MultifractalDetrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，c（k）的广义Hurst指数曲线用H（q）表示；cshuf（k）的广义Hurst指数曲线用Hshuf（q）表示；csurr（k）的广义Hurst指数曲线用Hsurr（q）表示；
3) 定义两个参数e1和e2，
Figure 284425DEST_PATH_IMAGE017
，
Figure 276651DEST_PATH_IMAGE018
，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由EEMD算法得到的信号分量。
进一步的，所述步骤1)中数据重排操作包括以下步骤：随机打乱分量c（k）的排列顺序。
进一步的，所述步骤1)中数据替代操作包括以下步骤：
1） 对分量c（k）执行离散傅里叶变换，获得分量c（k）的相位；
2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量c（k）的原始相位；
3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据cIFFT（k），求取数据cIFFT（k）的实部。
进一步的，所述步骤2)中MFDFA方法包括以下步骤：
1）构造x(k)，k=1,2,…,N，的轮廓Y(i)：
Figure 440916DEST_PATH_IMAGE019
Figure 151383DEST_PATH_IMAGE020
；
x(k)代表权利要求3所述步骤2）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure 160928DEST_PATH_IMAGE021
Figure 771775DEST_PATH_IMAGE022
；
y v (i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；本例中，m=1；
4） 计算第q阶波动函数的平均值：
Figure 841363DEST_PATH_IMAGE023
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值F q (s)和时间尺度s之间存在幂律关系：
F q (s)~s H(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 304705DEST_PATH_IMAGE012
；
6）对步骤5）中的公式两边取对数可得ln[F q (s)]=H(q)ln(s)+c ，c为常数，由此可以获得直线的斜率H(q)。
进一步的，所述步骤4中的最小二乘法包括以下步骤：对x（t），t=1, 2, …,n，x（t）代表步骤4中对cf（k）的局部极大值进行插值产生的序列或对cf（k）的局部极小值进行插值产生的序列，n代表插值序列的长度，
1）事先选定一组函数rk(t)，k=1,2,…,m，m<n，构造函数
f(t)=a1r1(t)+a2r2(t)+…+ amrm(t)，其中rk(t)代表二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线或复合函数曲线；
2）计算最小二乘指标
Figure 852361DEST_PATH_IMAGE024
；
3）令J对的ak偏导数
Figure 84759DEST_PATH_IMAGE025
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure 590827DEST_PATH_IMAGE026
，其中，RT代表R的转置矩阵，R-1代表R的逆矩阵；
4）rk(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标J，选取J最小时所对应的曲线形式为rk(t)的形式。
进一步的，所述步骤6中GZC方法包括以下步骤：
1）确定信号c(k)，k=1,2,…,N，的局部极大值点、局部极小值点和零点，c(k)= FMm(k)；
2）按照时间顺序统计连续两个局部极大值点、连续两个局部极小值点、连续两个上升零点和连续两个下降零点之间的时间间隔T1j，j=1,2,…,N1，N1代表时间间隔的数量，按照时间顺序统计相邻局部极大值点和局部极小值点、相邻局部极小值点和局部极大值点、相邻上升零点和下降零点、相邻下降零点和上升零点之间的时间间隔T2j，j=1,2,…,N2，N2代表时间间隔的数量，按照时间顺序统计相邻局部极大值点和下降零点、相邻下降零点和局部极小值点、相邻局部极小值点和上升零点、相邻上升零点和局部极大值点之间的时间间隔T3j，j=1,2,…,N3，N3代表时间间隔的数量，计算信号c(k)的平均周期T
Figure 10307DEST_PATH_IMAGE027
，则信号c(k)的瞬时尺度为sf=T。
基于以上所述的一种EEMD和GZC机器状态监测方法的装置，其特征在于：所述步骤16中状态监测装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现以上所述算法。
各步骤作用：
第1)步：采集振动信号；
第2)步：将原始信号分解成不同分量和的形式，其中有些分量对应噪声和趋势项，有些分量包含分形特征；
第3步：利用非线性判别算法去除信号分解结果中的噪声分量和趋势项，只保留包含分形特征的信号分量；
第4)~6)步：分离每个分形信号分量的频率调制部分，利用GZC估计每个分形信号分量的瞬时频率和瞬时尺度；
第7)步：根据分析尺度选择合适的分形信号分量，将选取的分形信号分量求和作为该分析尺度所对应的信号去趋势结果；
第8)~14)步：对每个分析尺度所对应的信号去趋势结果执行波动分析，得到原始信号的多重分形谱；
第15)步：提取多重分形谱的左端点、右端点和极值点的纵坐标，将这三个参数作为设备运行状态的特征参数；
第16)步：将上述算法部署在设备状态监测装置上，对设备状态进行监测；
本发明采用以上技术方案，与现有技术相比，本发明具有以下优点：
1) 采用EEMD方法自适应分解振动信号，根据非线性滤波方法去除噪声分量和趋势项，能够保护原始信号的分形结构，避免线性滤波方法对原始信号分形结构的破坏；
2) 分离信号分量的频率调制部分，利用GZC估计信号分量的瞬时频率，能够确保瞬时频率保持正值，避免了负频率现象；
3) 根据信号分量的瞬时频率计算相应的瞬时尺度，根据信号分量的瞬时尺度进行波动分析，避免了人工设定尺度的缺陷；
4) 利用EEMD方法自动确定信号趋势的类型，并保证信号趋势的连续性，有效解决了现有技术的缺陷；
5) 分析结果准确度和精确度高，设备运行状态识别结果正确率高。
下面结合附图和实施例对本发明做进一步说明。
附图说明
附图1为本发明实施例中本发明方法的流程图；
附图2为本发明实施例中设备状态监测装置的示意图；
附图3为本发明实施例中由多重分形级联模型产生的多重分形仿真信号；
附图4为本发明实施例中采用MFDFAemd方法获得的多重分形仿真信号的瞬时频率，信号分量数目为10；
附图5为本发明实施例中采用MFDFAoeemd方法获得的多重分形仿真信号的瞬时频率，信号分量数目为10；
附图6为本发明实施例中分别采用MFDFA、MFDFAemd和MFDFAoeemd方法对多重分形仿真信号分析结果的对比图；
附图7为本发明实施例中两个非线性判别参数的计算结果，“圆圈”和“方形”符号分别代表e1和e2；
附图8为本发明实施例中分别采用MFDFA、MFDFAemd和MFDFAoeemd方法对含噪多重分形仿真信号分析结果的对比图；
附图9为本发明实施例中由EEMD得到各个信号分量与原始信号的相关系数；
附图10为本发明实施例中分别采用MFDFA、基于相关滤波的MFDFAemd和基于相关滤波的MFDFAoeemd方法对含噪多重分形仿真信号分析结果的对比图；
附图11为本发明实施例中的五种齿轮箱振动信号，(a)~(e)分别代表正常、磨损、点蚀、断齿和磨损+点蚀齿轮状态；
附图12为本发明实施例中采用MFDFA得到的这五种齿轮箱振动信号的多重分形谱；
附图13为本发明实施例中采用MFDFAemd得到的这五种齿轮箱振动信号的多重分形谱；
附图14为本发明实施例中采用MFDFAoeemd得到的这五种齿轮箱振动信号的多重分形谱；
附图15为本发明实施例中由MFDFA得到的多重分形谱的左端点、右端点和极值点坐标对这五种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”、“菱形”和“左三角”符号分别代表正常、磨损、点蚀、断齿和磨损+点蚀齿轮状态；
附图16为本发明实施例中由MFDFAemd得到的多重分形谱的左端点、右端点和极值点坐标对这五种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”、“菱形”和“左三角”符号分别代表正常、磨损、点蚀、断齿和磨损+点蚀齿轮状态；
附图17为本发明实施例中由MFDFAoeemd得到的多重分形谱的左端点、右端点和极值点坐标对这五种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”、“菱形”和“左三角”符号分别代表正常、磨损、点蚀、断齿和磨损+点蚀齿轮状态。
具体实施方式
实施例，如图1、图2所示，一种EEMD和GZC机器状态监测方法，其特征在于：包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用集合经验模式分解（Ensemble Empirical Mode Decomposition,EEMD）算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure 627233DEST_PATH_IMAGE001
，其中，ci（k）代表由EEMD算法得到的第i个分量，rn（k）代表由EEMD算法得到的趋势项，本例中，n=10；
步骤3：采用非线性判别算法从EEMD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量cf（k）, f=1,2,…,p，p代表滤波后剩余分量的数量；
步骤4：确定cf（k）的局部极大值和局部极小值，采用Newton插值函数分别对cf（k）的局部极大值和局部极小值进行插值，利用最小二乘法分别拟合cf（k）的上包络u(k)和下包络l(k)，则cf（k）的包络定义为
Figure 448559DEST_PATH_IMAGE002
，符号|x|表示对x取绝对值；
步骤5：重复执行公式
Figure 125528DEST_PATH_IMAGE003
m次，j=1,2,…,m，直到
Figure 297883DEST_PATH_IMAGE004
，
得到cf（k）的频率调制部分FMm(k)，ej(k)代表cj(k)的包络，cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
步骤6：采用广义零交叉法（Generalized zero-crossing,GZC）计算FMm(k)的平均周期，得到cf（k）的瞬时尺度sf；
步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为
Figure 452921DEST_PATH_IMAGE028
；
步骤8：将Y s (k)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
步骤9：计算每段数据的方差：
Figure 158446DEST_PATH_IMAGE007
Figure 6317DEST_PATH_IMAGE030
；
步骤10：计算q阶函数：
Figure 400389DEST_PATH_IMAGE031
；
步骤11：改变s的取值，s=sf，f=1,2,…,p，重复上述步骤3到步骤10，得到关于q和s的方差函数Fq(s);
步骤12：如果x(k)存在分形特征，则F q (s)和尺度s之间存在幂律关系：F q (s)~s H(q)，H(q)代表x(k)的广义Hurst指数；
当q=0时， H(0)通过下式所定义的对数平均过程来确定：
Figure 359118DEST_PATH_IMAGE012
；
步骤13：计算信号x(k)的标准标度指数τ(q)=qH(q)-1，本例中，q在(-5, 5)范围内取值;
步骤14：计算信号x(k)的奇异指数α和多重分形谱f(α):
α=H(q)+q H’(q),
f(α)=q(α-H(q))+1，其中H’(q)代表H(q)的一阶导数；
步骤15：提取多重分形谱f(α)的左端点、右端点和极值点所对应的奇异指数，利用这3个参数来描述设备的运行状态；
步骤16：将上述步骤所述方法部署在状态监测装置上，对设备状态进行监测。
所述步骤2中EEMD算法包括以下步骤：
1）向数据x0(k)添加白噪声序列产生一个新数据xj(k) ：
Figure 155035DEST_PATH_IMAGE032
std[x0(k)]代表数据x0（k）的标准差，wnj(k)代表wnj中的第k个数据，wnj代表第j个随机产生的白噪声序列，wnj幅值为1，1≤j≤K；x0(k)代表权利要求1所述步骤2中x(k)；本例中，K=100；
2）对xj(k)执行经验模式分解，得到n个分量和一个趋势项
Figure 173807DEST_PATH_IMAGE033
cij(k)代表对xj(k)执行经验模式分解得到的第i个分量，rnj(k)代表对xj(k)执行经验模式分解得到的趋势项；
3）计算K次分解结果的平均值
Figure 55175DEST_PATH_IMAGE016
ci(k)表示对x0(k)进行集合经验模式分解得到的第i个分量，rn(k)表示对x0(k)进行集合经验模式分解得到的趋势项。
所述步骤3非线性判别算法包括以下步骤：
1) 对信号c（k）执行重排操作和替代操作，经重排操作得到的数据用cshuf（k）表示，替代操作后得到数据用csurr（k）表示；
2) 对c（k）、cshuf（k）和csurr（k）分别执行多重分形去趋势波动分析（MultifractalDetrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，c（k）的广义Hurst指数曲线用H（q）表示；cshuf（k）的广义Hurst指数曲线用Hshuf（q）表示；csurr（k）的广义Hurst指数曲线用Hsurr（q）表示；
3) 定义两个参数e1和e2，
Figure 817595DEST_PATH_IMAGE034
，
Figure 733598DEST_PATH_IMAGE035
，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由EEMD算法得到的信号分量。
所述步骤1)中数据重排操作包括以下步骤：随机打乱分量c（k）的排列顺序。
所述步骤1)中数据替代操作包括以下步骤：
1） 对分量c（k）执行离散傅里叶变换，获得分量c（k）的相位；
2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量c（k）的原始相位；
3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据cIFFT（k），求取数据cIFFT（k）的实部。
所述步骤2)中MFDFA方法包括以下步骤：
1）构造x(k)，k=1,2,…,N，的轮廓Y(i)：
Figure 657692DEST_PATH_IMAGE019
Figure 291936DEST_PATH_IMAGE020
；
x(k)代表权利要求3所述步骤2）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure 592467DEST_PATH_IMAGE021
Figure 362977DEST_PATH_IMAGE022
；
y v (i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；本例中，m=1；
4） 计算第q阶波动函数的平均值：
Figure 962366DEST_PATH_IMAGE023
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值F q (s)和时间尺度s之间存在幂律关系：
F q (s)~s H(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 83906DEST_PATH_IMAGE012
；
6）对步骤5）中的公式两边取对数可得ln[F q (s)]=H(q)ln(s)+c ，c为常数，由此可以获得直线的斜率H(q)。
所述步骤4中的最小二乘法包括以下步骤：对x（t），t=1, 2, …,n，x（t）代表步骤4中对cf（k）的局部极大值进行插值产生的序列或对cf（k）的局部极小值进行插值产生的序列，n代表插值序列的长度，
1）事先选定一组函数rk(t)，k=1,2,…,m，m<n，构造函数
f(t)=a1r1(t)+a2r2(t)+…+ amrm(t)，其中rk(t)代表二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线或复合函数曲线；
2）计算最小二乘指标
Figure 188128DEST_PATH_IMAGE036
；
3）令J对的ak偏导数
Figure DEST_PATH_IMAGE037
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure 281986DEST_PATH_IMAGE038
，其中，RT代表R的转置矩阵，R-1代表R的逆矩阵；
4）rk(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标J，选取J最小时所对应的曲线形式为rk(t)的形式。
所述步骤6中GZC方法包括以下步骤：
1）确定信号c(k)，k=1,2,…,N，的局部极大值点、局部极小值点和零点，c(k)= FMm(k)；
2）按照时间顺序统计连续两个局部极大值点、连续两个局部极小值点、连续两个上升零点和连续两个下降零点之间的时间间隔T1j，j=1,2,…,N1，N1代表时间间隔的数量，按照时间顺序统计相邻局部极大值点和局部极小值点、相邻局部极小值点和局部极大值点、相邻上升零点和下降零点、相邻下降零点和上升零点之间的时间间隔T2j，j=1,2,…,N2，N2代表时间间隔的数量，按照时间顺序统计相邻局部极大值点和下降零点、相邻下降零点和局部极小值点、相邻局部极小值点和上升零点、相邻上升零点和局部极大值点之间的时间间隔T3j，j=1,2,…,N3，N3代表时间间隔的数量，计算信号c(k)的平均周期T
Figure 813462DEST_PATH_IMAGE027
，则信号c(k)的瞬时尺度为sf=T。
基于以上所述的一种EEMD和GZC机器状态监测方法的装置，所述步骤16中状态监测装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现以上所述算法。
实验1 采用由多重分形级联模型产生的多重分形仿真信号对本发明所述算法的性能进行验证。
首先采用多重分形级联模型
Figure DEST_PATH_IMAGE039
产生的多重分形仿真信号对MFDFA、MFDFAemd和MFDFAoeemd的性能进行验证。本例中，p=0.375 ，n=14，产生的多重分形仿真信号如图3所示。采用MFDFAemd计算多重分形仿真信号的瞬时频率，结果如图4所示。从图4可以看出，由MFDFAemd计算得到的瞬时频率存在着许多负频率，由于负频率没有物理意义，因此MFDFAemd分析结果存在着较大的误差。采用MFDFAoeemd计算多重分形仿真信号的瞬时频率，结果如图5所示。从图5可以看出，由MFDFAoeemd计算得到的瞬时频率全部是正频率，因此MFDFAoeemd分析结果符合实际情况。接下来，分别采用MFDFA、MFDFAemd和MFDFAoeemd计算多重分形仿真信号的多重分形谱，结果如图6所示。根据图6所示结果，经计算可知， 由MFDFA得到的多重分形谱与理论值的绝对误差平均值为0.064，相对误差平均值为12.21%，由MFDFAemd得到的多重分形谱与理论值的绝对误差平均值为0.035，相对误差平均值为6.57%，由MFDFAoeemd得到的多重分形谱与理论值的绝对误差平均值为0.024，相对误差平均值为5.28%，因此由MFDFAoeemd得到的多重分形谱比由MFDFA得到的多重分形谱绝对误差平均值减小62.50%，相对误差平均值减小56.76%，由MFDFAoeemd得到的多重分形谱比由MFDFAemd得到的多重分形谱绝对误差平均值减小31.43%，相对误差平均值减小19.63%。图7所示为本发明实施例中两个非线性判别参数的计算结果，可以看出所有的信号分量都包含分形信号分量。然后，通过向上述多重分形仿真信号添加高斯白噪声的方法构造信噪比为20dB的含噪信号。分别采用MFDFA、MFDFAemd和MFDFAoeemd计算含噪多重分形仿真信号的多重分形谱，结果如图8所示。根据图8所示结果，由MFDFA获得的多重分形谱完全偏离了理论值，由MFDFAemd得到的多重分形谱与理论值的绝对误差平均值为0.084，相对误差平均值为11.81%，由MFDFAoeemd得到的多重分形谱与理论值的绝对误差平均值为0.035，相对误差平均值为5.51%，因此由MFDFAoeemd得到的多重分形谱比由MFDFAemd得到的多重分形谱绝对误差平均值减小58.33%，相对误差平均值减小53.34%。由此可见，与MFDFA和MFDFAemd相比，MFDFAoeemd具有更好地抗噪性。图9所示为本发明实施例中由EMD得到各个信号分量与原始信号的相关系数，可以看出第7个分量与原信号具有最弱的相关性，应该从原信号中去除。图10所示为本发明实施例中分别采用MFDFA、基于相关滤波的MFDFAemd和基于相关滤波MFDFAoeemd方法对含噪多重分形仿真信号分析结果的对比图。从图10可以看出，基于相关滤波的MFDFAemd和基于相关滤波MFDFAoeemd方法对含噪多重分形仿真信号分析结果完全偏离了理论值，因此相关滤波方法容易破坏原始信号的分形结构。
实验2 采用齿轮箱实验信号对本发明所述算法的性能进行验证。
本发明所用齿轮箱振动数据来自一个齿轮箱故障模拟实验台。这个实验模拟了三种单点齿轮故障：磨损、点蚀和断齿，和一种复合齿轮故障：磨损+点蚀。采集到的振动信号包含着正常、磨损、点蚀、断齿和磨损+点蚀五种齿轮运行状态。马达转速为1500RPM，振动信号采样频率为5120Hz，在每种齿轮状态下采集长度为10000点的5段数据，这五种齿轮箱振动信号如图11所示。首先采用MFDFA方法对这五种齿轮箱振动信号进行分析，得到这五种齿轮箱振动信号所对应的多重分形谱如图12所示，可以看出正常、磨损、断齿和复合故障状态所对应的多重分形谱严重重叠，而点蚀状态所对应的多重分形谱形状异常。然后采用MFDFAemd方法对这五种齿轮箱振动信号进行分析，得到这五种齿轮箱振动信号所对应的多重分形谱如图13所示，可以看出正常、磨损、断齿和复合故障状态所对应的多重分形谱严重重叠。最后采用MFDFAoeemd方法对这五种齿轮箱振动信号进行分析，得到这五种齿轮箱振动信号所对应的多重分形谱如图14所示，可以看出这五种齿轮状态所对应的多重分形谱可以分离。分别提取由MFDFA、MFDFAemd和MFDFAoeemd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这五种齿轮箱状态进行分类，结果分别如图15~17所示。从图15可以看出，利用由MFDFA方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分点蚀齿轮状态，但是不能区分剩余的四种齿轮状态，因此齿轮箱状态识别率为20%。从图16可以看出，利用由MFDFAemd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分点蚀齿轮状态，但是不能区分剩余的四种齿轮状态，因此齿轮箱状态识别率为20%。从图17可以看出，利用由MFDFAoeemd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分这五种齿轮状态，因此齿轮箱状态识别率为100%。可以看出，采用MFDFAoeemd方法可以将齿轮箱状态识别正确率提高80%。
根据试验结果，分析后认为：
1) 采用EEMD方法自适应分解振动信号，根据非线性滤波方法去除噪声分量和趋势项，能够保护原始信号的分形结构，避免线性滤波方法对原始信号分形结构的破坏；
2) 分离信号分量的频率调制部分，利用GZC估计信号分量的瞬时尺度，能够确保瞬时频率保持正值，避免了负频率现象；
3) 根据信号分量的瞬时尺度进行波动分析，避免了人工设定尺度的缺陷；
4) 利用EEMD方法自动确定信号趋势的类型，并保证信号趋势的连续性，有效解决了现有技术的缺陷；
5) 分析结果准确度和精确度高，设备运行状态识别结果正确率高。
本领域技术人员应该认识到，上述的具体实施方式只是示例性的，是为了使本领域技术人员能够更好的理解本发明内容，不应理解为是对本发明保护范围的限制，只要是根据本发明技术方案所作的改进，均落入本发明的保护范围。

Claims (9)

1.一种EEMD和GZC机器状态监测方法，其特征在于：包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用集合经验模式分解（Ensemble Empirical Mode Decomposition, EEMD）算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure 2605DEST_PATH_IMAGE001
，其中，ci（k）代表由EEMD算法得到的第i个分量，rn（k）代表由EEMD算法得到的趋势项；
步骤3：采用非线性判别算法从EEMD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量cf（k）, f=1,2,…,p，p代表滤波后剩余分量的数量；
步骤4：确定cf（k）的局部极大值和局部极小值，采用Newton插值函数分别对cf（k）的局部极大值和局部极小值进行插值，利用最小二乘法分别拟合cf（k）的上包络u(k)和下包络l(k)，则cf（k）的包络定义为
Figure 910518DEST_PATH_IMAGE002
，符号|x|表示对x取绝对值；
步骤5：重复执行公式
Figure 390041DEST_PATH_IMAGE003
m次，j=1,2,…,m，直到
Figure 92417DEST_PATH_IMAGE004
，
得到cf（k）的频率调制部分FMm(k)，ej(k)代表cj(k)的包络，cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
步骤6：采用广义零交叉法（Generalized zero-crossing,GZC）计算FMm(k)的平均周期，得到cf（k）的瞬时尺度sf；
步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为
Figure DEST_PATH_IMAGE005
；
步骤8：将Y s (k)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
步骤9：计算每段数据的方差：
Figure 890347DEST_PATH_IMAGE007
Figure 70792DEST_PATH_IMAGE009
；
步骤10：计算q阶函数：
Figure 670401DEST_PATH_IMAGE011
；
步骤11：改变s的取值，s=sf，f=1,2,…,p，重复上述步骤3到步骤10，得到关于q和s的方差函数Fq(s);
步骤12：如果x(k)存在分形特征，则F q (s)和尺度s之间存在幂律关系：F q (s)~s H(q)，H(q)代表x(k)的广义Hurst指数；
当q=0时， H(0)通过下式所定义的对数平均过程来确定：
Figure 543679DEST_PATH_IMAGE012
；
步骤13：计算信号x(k)的标准标度指数τ(q)=qH(q)-1；
步骤14：计算信号x(k)的奇异指数α和多重分形谱f(α):
α=H(q)+q H’(q),
f(α)=q(α-H(q))+1，其中H’(q)代表H(q)的一阶导数；
步骤15：提取多重分形谱f(α)的左端点、右端点和极值点所对应的奇异指数，利用这3个参数来描述设备的运行状态。
2.根据权利要求1所述的一种EEMD和GZC机器状态监测方法，其特征在于：所述步骤2中EEMD算法包括以下步骤：
1）向数据x0(k)添加白噪声序列产生一个新数据xj(k) ：
Figure 861528DEST_PATH_IMAGE014
std[x0(k)]代表数据x0（k）的标准差，wnj(k)代表wnj中的第k个数据，wnj代表第j个随机产生的白噪声序列，wnj幅值为1，1≤j≤K；x0(k)代表权利要求1所述步骤2中x(k)；
2）对xj(k)执行经验模式分解，得到n个分量和一个趋势项
Figure 845665DEST_PATH_IMAGE016
cij(k)代表对xj(k)执行经验模式分解得到的第i个分量，rnj(k)代表对xj(k)执行经验模式分解得到的趋势项；
3）计算K次分解结果的平均值
Figure 299780DEST_PATH_IMAGE017
ci(k)表示对x0(k)进行集合经验模式分解得到的第i个分量，rn(k)表示对x0(k)进行集合经验模式分解得到的趋势项。
3.根据权利要求1所述的一种EEMD和GZC机器状态监测方法，其特征在于：所述步骤3非线性判别算法包括以下步骤：
1) 对信号c（k）执行重排操作和替代操作，经重排操作得到的数据用cshuf（k）表示，替代操作后得到数据用csurr（k）表示；
2) 对c（k）、cshuf（k）和csurr（k）分别执行多重分形去趋势波动分析（MultifractalDetrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，c（k）的广义Hurst指数曲线用H（q）表示；cshuf（k）的广义Hurst指数曲线用Hshuf（q）表示；csurr（k）的广义Hurst指数曲线用Hsurr（q）表示；
3) 定义两个参数e1和e2，
Figure 78380DEST_PATH_IMAGE018
，
Figure 883525DEST_PATH_IMAGE019
，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由EEMD算法得到的信号分量。
4.根据权利要求3所述的一种EEMD和GZC机器状态监测方法，其特征在于：所述步骤1)中数据重排操作包括以下步骤：随机打乱分量c（k）的排列顺序。
5.根据权利要求3所述的一种EEMD和GZC机器状态监测方法，其特征在于：所述步骤1)中数据替代操作包括以下步骤：
1） 对分量c（k）执行离散傅里叶变换，获得分量c（k）的相位；
2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量c（k）的原始相位；
3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据cIFFT（k），求取数据cIFFT（k）的实部。
6.根据权利要求3所述的一种EEMD和GZC机器状态监测方法，其特征在于：所述步骤2)中MFDFA方法包括以下步骤：
1）构造x(k)，k=1,2,…,N，的轮廓Y(i)：
Figure 671352DEST_PATH_IMAGE021
Figure 478509DEST_PATH_IMAGE022
；
x(k)代表权利要求3所述步骤2）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure 693590DEST_PATH_IMAGE023
Figure 720451DEST_PATH_IMAGE024
；
y v (i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；
4） 计算第q阶波动函数的平均值：
Figure 46391DEST_PATH_IMAGE025
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值F q (s)和时间尺度s之间存在幂律关系：
F q (s)~s H(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 475098DEST_PATH_IMAGE012
；
6）对步骤5）中的公式两边取对数可得ln[F q (s)]=H(q)ln(s)+c ，c为常数，由此可以获得直线的斜率H(q)。
7.根据权利要求1所述的一种EEMD和GZC机器状态监测方法，其特征在于：所述步骤4中的最小二乘法包括以下步骤：对x（t），t=1, 2, …,n，x（t）代表步骤4中对cf（k）的局部极大值进行插值产生的序列或对cf（k）的局部极小值进行插值产生的序列，n代表插值序列的长度，
1）事先选定一组函数rk(t)，k=1,2,…,m，m<n，构造函数
f(t)=a1r1(t)+a2r2(t)+…+ amrm(t)，其中rk(t)代表二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线或复合函数曲线；
2）计算最小二乘指标
Figure 595501DEST_PATH_IMAGE026
；
3）令J对的ak偏导数
Figure DEST_PATH_IMAGE027
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure 109659DEST_PATH_IMAGE028
，其中，RT代表R的转置矩阵，R-1代表R的逆矩阵；
4）rk(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标J，选取J最小时所对应的曲线形式为rk(t)的形式。
8.根据权利要求1所述的一种EEMD和GZC机器状态监测方法，其特征在于：所述步骤6中GZC方法包括以下步骤：
1）确定信号c(k)，k=1,2,…,N，的局部极大值点、局部极小值点和零点，c(k)= FMm(k)；
2）按照时间顺序统计连续两个局部极大值点、连续两个局部极小值点、连续两个上升零点和连续两个下降零点之间的时间间隔T1j，j=1,2,…,N1，N1代表时间间隔的数量，按照时间顺序统计相邻局部极大值点和局部极小值点、相邻局部极小值点和局部极大值点、相邻上升零点和下降零点、相邻下降零点和上升零点之间的时间间隔T2j，j=1,2,…,N2，N2代表时间间隔的数量，按照时间顺序统计相邻局部极大值点和下降零点、相邻下降零点和局部极小值点、相邻局部极小值点和上升零点、相邻上升零点和局部极大值点之间的时间间隔T3j，j=1,2,…,N3，N3代表时间间隔的数量，计算信号c(k)的平均周期T
Figure DEST_PATH_IMAGE029
，则信号c(k)的瞬时尺度为sf=T。
9.实施如权利要求1-8之一所述的一种EEMD和GZC机器状态监测方法的装置，其特征在于：所述装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现所述算法。
CN202011238756.1A 2020-11-09 2020-11-09 一种eemd和gzc机器状态监测方法及装置 Withdrawn CN112683391A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011238756.1A CN112683391A (zh) 2020-11-09 2020-11-09 一种eemd和gzc机器状态监测方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011238756.1A CN112683391A (zh) 2020-11-09 2020-11-09 一种eemd和gzc机器状态监测方法及装置

Publications (1)

Publication Number Publication Date
CN112683391A true CN112683391A (zh) 2021-04-20

Family

ID=75446001

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011238756.1A Withdrawn CN112683391A (zh) 2020-11-09 2020-11-09 一种eemd和gzc机器状态监测方法及装置

Country Status (1)

Country Link
CN (1) CN112683391A (zh)

Similar Documents

Publication Publication Date Title
Yan et al. Improved Hilbert–Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis
CN108802525A (zh) 基于小样本的设备故障智能预测方法
CN110046476B (zh) 滚动轴承故障的三元二进分形小波稀疏诊断方法
CN106168538B (zh) 一种滚动轴承的itd、谱峭度和平滑迭代包络分析方法
CN110717472B (zh) 基于改进的小波阈值去噪的故障诊断方法及系统
CN112683393A (zh) 一种lmd设备故障诊断方法及系统
CN112697484A (zh) 一种ssd多尺度波动分析状态监测方法及装置
CN112697266A (zh) 一种itd和gzc机器状态监测方法及装置
CN113375939A (zh) 基于svd和vmd的机械件故障诊断方法
CN112697483A (zh) 一种lmd多尺度波动分析状态监测方法及装置
CN112697482A (zh) 一种vmd多尺度波动分析状态监测方法及装置
CN112697479A (zh) 一种itd多尺度波动分析状态监测方法及装置
CN112683395A (zh) 一种elmd和gzc机器状态监测方法及装置
CN116773961A (zh) 基于振动信号高频特征分析的输电线路腐蚀检测方法
CN112881046A (zh) 一种vmd和gzc机器状态监测方法及装置
CN112697473A (zh) 一种emd和gzc机器状态监测方法及装置
CN106096199B (zh) 一种滚动轴承的wt、谱峭度和平滑迭代包络分析方法
CN112697477A (zh) 一种emd设备故障诊断方法及系统
CN112697470A (zh) 一种ssd设备故障诊断方法及系统
CN118375603A (zh) 一种真空泵的故障监测方法及系统
CN112697471A (zh) 一种ssd和gzc机器状态监测方法及装置
CN112697265A (zh) 一种自适应多重分形方法及设备工况监测装置
CN112697472A (zh) 一种wd多尺度波动分析状态监测方法及装置
CN112697475A (zh) 一种itd设备故障诊断方法及系统
CN112683391A (zh) 一种eemd和gzc机器状态监测方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
WW01 Invention patent application withdrawn after publication

Application publication date: 20210420

WW01 Invention patent application withdrawn after publication