CN112697263A - 一种eemd多尺度波动分析状态监测方法及装置 - Google Patents

一种eemd多尺度波动分析状态监测方法及装置 Download PDF

Info

Publication number
CN112697263A
CN112697263A CN202011238759.5A CN202011238759A CN112697263A CN 112697263 A CN112697263 A CN 112697263A CN 202011238759 A CN202011238759 A CN 202011238759A CN 112697263 A CN112697263 A CN 112697263A
Authority
CN
China
Prior art keywords
data
signal
eemd
scale
fractal
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
CN202011238759.5A
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 CN202011238759.5A priority Critical patent/CN112697263A/zh
Publication of CN112697263A publication Critical patent/CN112697263A/zh
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H17/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the preceding groups
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations

Abstract

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

Description

一种EEMD多尺度波动分析状态监测方法及装置
技术领域
本发明涉及设备状态监测与故障诊断领域，具体涉及一种EEMD多尺度波动分析状态监测方法及装置。
背景技术
设备振动信号包含丰富的分形特征，这些分形特征能够描述设备的运行状态。盒维数、功率谱分析和重标极差方法可以估计平稳信号的单重分形参数，去趋势波动分析（DFA）能够估计非平稳信号的单重分形维数。然而，设备出现故障时，其振动信号通常是非平稳的，且具有多重分形特征，这时传统的分形维数估计方法会产生比较大的误差。多重分形去趋势波动分析（MFDFA）能够估计非平稳信号的多重分形参数，但是MFDFA方法存在着分析尺度需要人工确定、拟合多项式趋势阶数难以确定和数据段之间不连续的问题。目前，已经有文献提出了基于EMD的MFDFA版本（MFDFAemd），用来解决MFDFA存在的问题。然而，MFDFAemd采用的线性滤波方法容易破坏原始信号的分形结构，且存在着负频率现象，这些缺陷严重影响了MFDFAemd的应用效果。综上所述，现有技术难以准确提取设备振动信号的多重分形特征，难以准确检测设备运行状态。
发明内容
本发明要解决的问题是针对以上不足，提出一种EEMD多尺度波动分析状态监测方法及装置（本发明提出的方法简称为MFDFAoeemd）。采用本发明所提出的方法对设备振动信号进行分析，能够有效提取设备振动信号的多重分形特征，克服MFDFA方法存在的分析尺度需要人工确定、拟合多项式趋势阶数难以确定和数据段之间不连续的问题，解决MFDFAemd方法存在的原始信号分形结构破坏和负频率现象，具有分析结果准确度和精确度高，设备运行状态识别结果正确率高等优点。
为解决以上技术问题，本发明采取的技术方案如下：一种EEMD多尺度波动分析状态监测方法，其特征在于：包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用集合经验模式分解(Ensemble Empirical Mode Decomposition,EEMD)算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure DEST_PATH_IMAGE001
，其中，ci（k）代表由EEMD算法得到的第i个分量，rn（k）代表由EEMD算法得到的趋势项，本例中，n=10；
步骤3：采用非线性判别算法从EEMD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量cf（k）, f=1,2,…,p，p代表滤波后剩余分量的数量；
步骤4：确定cf（k）的局部极大值和局部极小值，采用三次样条插值函数分别对cf（k）的局部极大值和局部极小值进行插值，采用最小二乘法分别拟合cf（k）的上包络u(k)和下包络l(k)，则cf（k）的包络定义为
Figure 863068DEST_PATH_IMAGE002
，符号|x|表示对x取绝对值；
步骤5：重复执行公式
Figure DEST_PATH_IMAGE003
m次，j=1,2,…,m，直到
Figure 221369DEST_PATH_IMAGE004
，得到cf（k）的频率调制部分FMm(k)，ej(k)代表cj(k)的包络，cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
步骤6：采用Teager能量算子（Teager Energy Operator, TEO）计算FMm(k)的瞬时频率，获得cf（k）的瞬时频率instff(k)，得到cf（k）的瞬时尺度
Figure DEST_PATH_IMAGE005
；
步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为
Figure 506988DEST_PATH_IMAGE006
；
步骤8：将Ys(k)分成不重叠的Ns段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2Ns段数据；
步骤9：计算每段数据的方差：
Figure DEST_PATH_IMAGE007
Figure 682754DEST_PATH_IMAGE008
；
步骤10：计算q阶函数：
Figure DEST_PATH_IMAGE009
；
步骤11：改变s的取值，s=sf，f=1,2,…,p，重复上述步骤3到步骤10，得到关于q和s的方差函数Fq(s);
步骤12：如果x(k)存在分形特征，则Fq(s)和尺度s之间存在幂律关系：Fq(s)~sH(q)，H(q)代表x(k)的广义Hurst指数；
当q=0时， H(0)通过下式所定义的对数平均过程来确定：
Figure 418104DEST_PATH_IMAGE010
；
步骤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 DEST_PATH_IMAGE011
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 704729DEST_PATH_IMAGE012
cij(k)代表对xj(k)执行经验模式分解得到的第i个分量，rnj(k)代表对xj(k)执行经验模式分解得到的趋势项；
3）计算K次分解结果的平均值
Figure DEST_PATH_IMAGE013
ci(k)表示对x0(k)进行集合经验模式分解得到的第i个分量，rn(k)表示对x0(k)进行集合经验模式分解得到的趋势项。
进一步的，所述步骤3非线性判别算法包括以下步骤：
1) 对信号c（k）执行重排操作和替代操作，经重排操作得到的数据用cshuf（k）表示，替代操作后得到数据用csurr（k）表示；
2) 对c（k）、cshuf（k）和csurr（k）分别执行多重分形去趋势波动分析（Multifractal Detrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，c（k）的广义Hurst指数曲线用H（q）表示；cshuf（k）的广义Hurst指数曲线用Hshuf（q）表示；csurr（k）的广义Hurst指数曲线用Hsurr（q）表示；
3) 定义两个参数e1和e2，
Figure 579275DEST_PATH_IMAGE014
，
Figure DEST_PATH_IMAGE015
，如果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 925943DEST_PATH_IMAGE016
Figure DEST_PATH_IMAGE017
；
x(k)代表权利要求3所述步骤2）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的Ns段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2Ns段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure 151519DEST_PATH_IMAGE018
Figure DEST_PATH_IMAGE019
；
yv(i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；本例中，m=1；
4） 计算第q阶波动函数的平均值：
Figure 54884DEST_PATH_IMAGE020
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值Fq(s)和时间尺度s之间存在幂律关系：
Fq(s)~sH(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 298784DEST_PATH_IMAGE010
；
6）对步骤5）中的公式两边取对数可得ln[Fq(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 DEST_PATH_IMAGE021
；
3）令J对的ak偏导数
Figure 322014DEST_PATH_IMAGE022
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure DEST_PATH_IMAGE023
，其中，RT代表R的转置矩阵，R-1代表R的逆矩阵；
4）rk(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标J，选取J最小时所对应的曲线形式为rk(t)的形式。
进一步的，所述步骤6中Teager能量算子方法包括以下步骤：
1）对信号c(k)，k=1,2,…,N，c(k)= FMm(k)，构建函数
ψ(c(k))=c2(k)- c(k+1) c(k-1);
2）令d(k)=c(k)-c(k-1)，信号c(k)的瞬时频率instf(k)定义为：
Figure 97203DEST_PATH_IMAGE024
。
基于以上所述的一种EEMD多尺度波动分析状态监测方法，实现该方法的装置，所述步骤16中状态监测装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现以上所述算法。
各步骤作用：
第1)步：采集振动信号；
第2)步：将原始信号分解成不同分量和的形式，其中有些分量对应噪声和趋势项，有些分量包含分形特征；
第3步：利用非线性判别算法去除信号分解结果中的噪声分量和趋势项，只保留包含分形特征的信号分量；
第4)~6)步：分离每个分形信号分量的频率调制部分，利用TEO估计每个分形信号分量的瞬时频率和瞬时尺度；
第7)步：根据分析尺度选择合适的分形信号分量，将选取的分形信号分量求和作为该分析尺度所对应的信号去趋势结果；
第8)~14)步：对每个分析尺度所对应的信号去趋势结果执行波动分析，得到原始信号的多重分形谱；
第15)步：提取多重分形谱的左端点、右端点和极值点所对应的奇异指数，将这三个参数作为设备运行状态的特征参数；
第16)步：将上述算法部署在设备状态监测装置上，对设备状态进行监测；
本发明采用以上技术方案，与现有技术相比，本发明具有以下优点：
1) 采用EEMD方法分解振动信号，根据非线性滤波方法去除噪声分量和趋势项，能够保护原始信号的分形结构，避免线性滤波方法对原始信号分形结构的破坏；
2) 分离信号分量的频率调制部分，利用TEO估计信号分量的瞬时频率，能够确保瞬时频率保持正值，避免了负频率现象；
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为本发明实施例中由EMD得到各个信号分量与原始信号的相关系数；
附图10为本发明实施例中分别采用MFDFA、基于相关滤波的MFDFAemd和基于相关滤波的MFDFAoeemd方法对含噪多重分形仿真信号分析结果的对比图；
附图11为本发明实施例中的四种齿轮箱振动信号，(a)~(d)分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；
附图12为本发明实施例中采用MFDFA得到的这四种齿轮箱振动信号的多重分形谱；
附图13为本发明实施例中采用MFDFAemd得到的这四种齿轮箱振动信号的多重分形谱；
附图14为本发明实施例中采用MFDFAoeemd得到的这四种齿轮箱振动信号的多重分形谱；
附图15为本发明实施例中由MFDFA得到的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；
附图16为本发明实施例中由MFDFAemd得到的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；
附图17为本发明实施例中由MFDFAoeemd得到的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态。
具体实施方式
实施例，如图1、图2所示，一种EEMD多尺度波动分析状态监测方法，包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用集合经验模式分解(Ensemble Empirical Mode Decomposition,EEMD)算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure 663313DEST_PATH_IMAGE001
，其中，ci（k）代表由EEMD算法得到的第i个分量，rn（k）代表由EEMD算法得到的趋势项，本例中，n=10；
步骤3：采用非线性判别算法从EEMD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量cf（k）, f=1,2,…,p，p代表滤波后剩余分量的数量；
步骤4：确定cf（k）的局部极大值和局部极小值，采用三次样条插值函数分别对cf（k）的局部极大值和局部极小值进行插值，采用最小二乘法分别拟合cf（k）的上包络u(k)和下包络l(k)，则cf（k）的包络定义为
Figure 27299DEST_PATH_IMAGE002
，符号|x|表示对x取绝对值；
步骤5：重复执行公式
Figure 387873DEST_PATH_IMAGE003
m次，j=1,2,…,m，直到
Figure 588041DEST_PATH_IMAGE004
，得到cf（k）的频率调制部分FMm(k)，ej(k)代表cj(k)的包络，cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
步骤6：采用Teager能量算子（Teager Energy Operator, TEO）计算FMm(k)的瞬时频率，获得cf（k）的瞬时频率instff(k)，得到cf（k）的瞬时尺度
Figure 957843DEST_PATH_IMAGE005
；
步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为
Figure 910755DEST_PATH_IMAGE006
；
步骤8：将Ys(k)分成不重叠的Ns段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2Ns段数据；
步骤9：计算每段数据的方差：
Figure 520859DEST_PATH_IMAGE007
Figure 129695DEST_PATH_IMAGE008
；
步骤10：计算q阶函数：
Figure 365504DEST_PATH_IMAGE009
；
步骤11：改变s的取值，s=sf，f=1,2,…,p，重复上述步骤3到步骤10，得到关于q和s的方差函数Fq(s);
步骤12：如果x(k)存在分形特征，则Fq(s)和尺度s之间存在幂律关系：Fq(s)~sH(q)，H(q)代表x(k)的广义Hurst指数；
当q=0时， H(0)通过下式所定义的对数平均过程来确定：
Figure 376186DEST_PATH_IMAGE010
；
步骤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 154261DEST_PATH_IMAGE011
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 250393DEST_PATH_IMAGE012
cij(k)代表对xj(k)执行经验模式分解得到的第i个分量，rnj(k)代表对xj(k)执行经验模式分解得到的趋势项；
3）计算K次分解结果的平均值
Figure 961997DEST_PATH_IMAGE013
ci(k)表示对x0(k)进行集合经验模式分解得到的第i个分量，rn(k)表示对x0(k)进行集合经验模式分解得到的趋势项。
所述步骤3非线性判别算法包括以下步骤：
1) 对信号c（k）执行重排操作和替代操作，经重排操作得到的数据用cshuf（k）表示，替代操作后得到数据用csurr（k）表示；
2) 对c（k）、cshuf（k）和csurr（k）分别执行多重分形去趋势波动分析（Multifractal Detrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，c（k）的广义Hurst指数曲线用H（q）表示；cshuf（k）的广义Hurst指数曲线用Hshuf（q）表示；csurr（k）的广义Hurst指数曲线用Hsurr（q）表示；
3) 定义两个参数e1和e2，
Figure 889502DEST_PATH_IMAGE014
，
Figure 841409DEST_PATH_IMAGE015
，如果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 424837DEST_PATH_IMAGE016
Figure 736869DEST_PATH_IMAGE017
；
x(k)代表权利要求3所述步骤2）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的Ns段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2Ns段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure 722143DEST_PATH_IMAGE018
Figure 844951DEST_PATH_IMAGE019
；
yv(i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；本例中，m=1；
4） 计算第q阶波动函数的平均值：
Figure 977992DEST_PATH_IMAGE020
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值Fq(s)和时间尺度s之间存在幂律关系：
Fq(s)~sH(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 31399DEST_PATH_IMAGE010
；
6）对步骤5）中的公式两边取对数可得ln[Fq(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 871179DEST_PATH_IMAGE021
；
3）令J对的ak偏导数
Figure 164888DEST_PATH_IMAGE022
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure 457329DEST_PATH_IMAGE023
，其中，RT代表R的转置矩阵，R-1代表R的逆矩阵；
4）rk(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标J，选取J最小时所对应的曲线形式为rk(t)的形式。
所述步骤6中Teager能量算子方法包括以下步骤：
1）对信号c(k)，k=1,2,…,N，c(k)= FMm(k)，构建函数
ψ(c(k))=c2(k)- c(k+1) c(k-1);
2）令d(k)=c(k)-c(k-1)，信号c(k)的瞬时频率instf(k)定义为：
Figure 376743DEST_PATH_IMAGE024
。
基于以上所述的一种EEMD多尺度波动分析状态监测方法，实现该方法的装置，所述步骤16中状态监测装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现以上所述算法。
利用传动齿轮箱振动数据对本发明所述算法的性能进行验证。
实验1 采用由多重分形级联模型产生的多重分形仿真信号对本发明所述算法的性能进行验证。
首先采用多重分形级联模型
Figure DEST_PATH_IMAGE025
产生的多重分形仿真信号对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.022，相对误差平均值为5.19%，因此由MFDFAoeemd得到的多重分形谱比由MFDFA得到的多重分形谱绝对误差平均值减小65.63%，相对误差平均值减小57.49%，由MFDFAoeemd得到的多重分形谱比由MFDFAemd得到的多重分形谱绝对误差平均值减小37.14%，相对误差平均值减小21.00%。图7所示为本发明实施例中两个非线性判别参数的计算结果，可以看出所有的信号分量都包含分形信号分量。然后，通过向上述多重分形仿真信号添加高斯白噪声的方法构造信噪比为20dB的含噪信号。分别采用MFDFA、MFDFAemd和MFDFAoeemd计算含噪多重分形仿真信号的多重分形谱，结果如图8所示。根据图8所示结果，由MFDFA获得的多重分形谱完全偏离了理论值，由MFDFAemd得到的多重分形谱与理论值的绝对误差平均值为0.084，相对误差平均值为11.81%，由MFDFAoeemd得到的多重分形谱与理论值的绝对误差平均值为0.037，相对误差平均值为5.42%，因此由MFDFAoeemd得到的多重分形谱比由MFDFAemd得到的多重分形谱绝对误差平均值减小55.95%，相对误差平均值减小54.11%。由此可见，与MFDFA和MFDFAemd相比，MFDFAoeemd具有更好地抗噪性。图9所示为本发明实施例中由EMD得到各个信号分量与原始信号的相关系数，可以看出第7个分量与原信号具有最弱的相关性，应该从原信号中去除。图10所示为本发明实施例中分别采用MFDFA、基于相关滤波的MFDFAemd和基于相关滤波的MFDFAoeemd方法对含噪多重分形仿真信号分析结果的对比图。从图10可以看出，基于相关滤波的MFDFAemd和基于相关滤波的MFDFAoeemd方法对含噪多重分形仿真信号分析结果完全偏离了理论值，因此相关滤波方法容易破坏原始信号的分形结构。
实验2 采用齿轮箱实验信号对本发明所述算法的性能进行验证。
本发明所用齿轮箱振动数据一个齿轮箱故障模拟实验台。这个实验通过在一个轮齿的齿根上制造不同程度的划痕，直到最后完全破坏掉这个轮齿，来模拟轮齿从正常状态到失效的过程。实验中使用的齿轮箱为二级齿轮传动，从输入端到输出端的齿轮齿数分别为25、40、22和55，故障轮齿位于输入轴齿轮上，驱动电机转速为2000RMP，齿轮箱振动信号由位于输入端外壳上的加速度传感器测取。采集到的振动信号包含着正常、轻度划痕、重度划痕和断齿四种故障状态，这些振动信号在一定程度上代表了轮齿从正常到失效的过程。振动信号采样频率为16384Hz，在每种齿轮箱状态下采集长度为10000点的20段数据，，这四种齿轮箱振动信号如图11所示。首先采用MFDFA方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图12所示，可以看出轻度划痕和重度划痕所对应的多重分形谱严重重叠。然后采用MFDFAemd方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图13所示，可以看出这四种齿轮箱状态所对应的多重分形谱严重重叠。最后采用MFDFAoeemd方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图14所示，可以看出断齿状态振动信号的多重分形谱与其他三种齿轮箱状态所对应的多重分形谱明显不同，而正常、轻度划痕和重度划痕齿轮状态所对应的多重分形谱能够在α<0.4时清楚分离。分别提取由MFDFA、MFDFAemd和MFDFAoeemd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态进行分类，结果分别如图15~17所示。从图15可以看出，利用由MFDFA方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分正常状态和断齿状态，但是不能区分轻度划痕和重度划痕状态，因此齿轮箱状态识别率为50%。从图16可以看出，利用由MFDFAemd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分正常状态和断齿状态，但是不能区分轻度划痕和重度划痕状态，因此齿轮箱状态识别率为50%。从图17可以看出，利用由MFDFAoeemd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分这四种齿轮箱状态，因此齿轮箱状态识别率为100%。可以看出，采用MFDFAoeemd方法可以将齿轮箱状态识别正确率提高50%。
根据实验结果，分析后认为：
1) 采用EEMD方法自适应分解振动信号，根据非线性滤波方法去除噪声分量和趋势项，能够保护原始信号的分形结构，避免线性滤波方法对原始信号分形结构的破坏；
2) 分离信号分量的频率调制部分，利用TEO估计信号分量的瞬时频率，能够确保瞬时频率保持正值，避免了负频率现象；
3) 根据信号分量的瞬时频率计算相应的瞬时尺度，根据信号分量的瞬时尺度进行波动分析，避免了人工设定尺度的缺陷；
4) 利用EEMD方法自动确定信号趋势的类型，并保证信号趋势的连续性，有效解决了现有技术的缺陷；
5) 分析结果准确度和精确度高，设备运行状态识别结果正确率高。
本领域技术人员应该认识到，上述的具体实施方式只是示例性的，是为了使本领域技术人员能够更好的理解本发明内容，不应理解为是对本发明保护范围的限制，只要是根据本发明技术方案所作的改进，均落入本发明的保护范围。

Claims (9)

1.一种EEMD多尺度波动分析状态监测方法，其特征在于：包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用集合经验模式分解(Ensemble Empirical Mode Decomposition, EEMD)算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure DEST_PATH_IMAGE002
，其中，ci（k）代表由EEMD算法得到的第i个分量，rn（k）代表由EEMD算法得到的趋势项；
步骤3：采用非线性判别算法从EEMD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量cf（k）, f=1,2,…,p，p代表滤波后剩余分量的数量；
步骤4：确定cf（k）的局部极大值和局部极小值，采用三次样条插值函数分别对cf（k）的局部极大值和局部极小值进行插值，采用最小二乘法分别拟合cf（k）的上包络u(k)和下包络l(k)，则cf（k）的包络定义为
Figure DEST_PATH_IMAGE004
，符号|x|表示对x取绝对值；
步骤5：重复执行公式
Figure DEST_PATH_IMAGE006
m次，j=1,2,…,m，直到
Figure DEST_PATH_IMAGE008
，得到cf（k）的频率调制部分FMm(k)，ej(k)代表cj(k)的包络，cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
步骤6：采用Teager能量算子（Teager Energy Operator, TEO）计算FMm(k)的瞬时频率，获得cf（k）的瞬时频率instff(k)，得到cf（k）的瞬时尺度
Figure DEST_PATH_IMAGE010
；
步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为
Figure DEST_PATH_IMAGE012
；
步骤8：将Y s (k)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
步骤9：计算每段数据的方差：
Figure DEST_PATH_IMAGE014
Figure DEST_PATH_IMAGE016
；
步骤10：计算q阶函数：
Figure DEST_PATH_IMAGE018
；
步骤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 DEST_PATH_IMAGE020
；
步骤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多尺度波动分析状态监测方法，其特征在于：所述步骤2的EEMD算法包括以下步骤：
1）向数据x0(k)添加白噪声序列产生一个新数据xj(k) ：
Figure DEST_PATH_IMAGE022
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 DEST_PATH_IMAGE024
cij(k)代表对xj(k)执行经验模式分解得到的第i个分量，rnj(k)代表对xj(k)执行经验模式分解得到的趋势项；
3）计算K次分解结果的平均值
Figure DEST_PATH_IMAGE026
ci(k)表示对x0(k)进行集合经验模式分解得到的第i个分量，rn(k)表示对x0(k)进行集合经验模式分解得到的趋势项。
3.根据权利要求1所述的一种EEMD多尺度波动分析状态监测方法，其特征在于:所述步骤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 DEST_PATH_IMAGE028
，
Figure DEST_PATH_IMAGE030
，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由EEMD算法得到的信号分量。
4.根据权利要求3所述的一种EEMD多尺度波动分析状态监测方法，其特征在于：所述步骤1)中数据重排操作包括以下步骤：随机打乱分量c（k）的排列顺序。
5.根据权利要求3所述的一种EEMD多尺度波动分析状态监测方法，其特征在于：所述步骤1)中数据替代操作包括以下步骤：
1） 对分量c（k）执行离散傅里叶变换，获得分量c（k）的相位；
2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量c（k）的原始相位；
3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据cIFFT（k），求取数据cIFFT（k）的实部。
6.根据权利要求3所述的一种EEMD多尺度波动分析状态监测方法，其特征在于：所述步骤2)中MFDFA方法包括以下步骤：
1）构造x(k)，k=1,2,…,N，的轮廓Y(i)：
Figure DEST_PATH_IMAGE032
Figure DEST_PATH_IMAGE034
；
x(k)代表权利要求3所述步骤2）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure DEST_PATH_IMAGE036
Figure DEST_PATH_IMAGE038
；
y v (i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm； 4） 计算第q阶波动函数的平均值：
Figure DEST_PATH_IMAGE040
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值F q (s)和时间尺度s之间存在幂律关系：
F q (s)~s H(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 350678DEST_PATH_IMAGE020
；
6）对步骤5）中的公式两边取对数可得ln[F q (s)]=H(q)ln(s)+c ，c为常数，由此可以获得直线的斜率H(q)。
7. 根据权利要求1所述的一种EEMD多尺度波动分析状态监测方法，其特征在于:所述步骤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 DEST_PATH_IMAGE042
；
3）令J对的ak偏导数
Figure DEST_PATH_IMAGE044
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure DEST_PATH_IMAGE046
，其中，RT代表R的转置矩阵，R-1代表R的逆矩阵；
4）rk(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标J，选取J最小时所对应的曲线形式为rk(t)的形式。
8. 根据权利要求1所述的一种EEMD多尺度波动分析状态监测方法，其特征在于：所述步骤6中Teager能量算子方法包括以下步骤：
1）对信号c(k)，k=1,2,…,N，c(k)= FMm(k)，构建函数
ψ(c(k))=c2(k)- c(k+1) c(k-1);
2）令d(k)=c(k)-c(k-1)，信号c(k)的瞬时频率instf(k)定义为：
Figure DEST_PATH_IMAGE048
。
9.实施如权利要求1-8之一所述的一种EEMD多尺度波动分析状态监测方法的装置，其特征在于：所述装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现所述算法。
CN202011238759.5A 2020-11-09 2020-11-09 一种eemd多尺度波动分析状态监测方法及装置 Withdrawn CN112697263A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011238759.5A CN112697263A (zh) 2020-11-09 2020-11-09 一种eemd多尺度波动分析状态监测方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011238759.5A CN112697263A (zh) 2020-11-09 2020-11-09 一种eemd多尺度波动分析状态监测方法及装置

Publications (1)

Publication Number Publication Date
CN112697263A true CN112697263A (zh) 2021-04-23

Family

ID=75507000

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011238759.5A Withdrawn CN112697263A (zh) 2020-11-09 2020-11-09 一种eemd多尺度波动分析状态监测方法及装置

Country Status (1)

Country Link
CN (1) CN112697263A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113324646A (zh) * 2021-05-26 2021-08-31 兰州交通大学 一种大风区电气化铁路接触网正馈线舞动定位算法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113324646A (zh) * 2021-05-26 2021-08-31 兰州交通大学 一种大风区电气化铁路接触网正馈线舞动定位算法

Similar Documents

Publication Publication Date Title
CN107505135B (zh) 一种滚动轴承复合故障提取方法及系统
Yan et al. Improved Hilbert–Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis
CN106017926A (zh) 基于变模态分解的滚动轴承故障诊断方法
CN113375939B (zh) 基于svd和vmd的机械件故障诊断方法
CN101587017A (zh) 一种基于局部均值分解循环频率谱的齿轮故障诊断方法
CN112697484A (zh) 一种ssd多尺度波动分析状态监测方法及装置
CN112697482A (zh) 一种vmd多尺度波动分析状态监测方法及装置
CN112697483A (zh) 一种lmd多尺度波动分析状态监测方法及装置
CN110084208B (zh) 一种自适应降噪并避免阶次混叠的计算阶次跟踪方法
CN111175045B (zh) 一种机车牵引电机轴承的振动加速度数据的清洗方法
CN112683393A (zh) 一种lmd设备故障诊断方法及系统
CN112697263A (zh) 一种eemd多尺度波动分析状态监测方法及装置
CN106096199B (zh) 一种滚动轴承的wt、谱峭度和平滑迭代包络分析方法
CN112697472A (zh) 一种wd多尺度波动分析状态监测方法及装置
CN112697479A (zh) 一种itd多尺度波动分析状态监测方法及装置
CN106198009B (zh) 一种滚动轴承的emd、谱峭度和平滑迭代包络分析方法
CN112697473A (zh) 一种emd和gzc机器状态监测方法及装置
CN112697471A (zh) 一种ssd和gzc机器状态监测方法及装置
CN112697470A (zh) 一种ssd设备故障诊断方法及系统
CN112697266A (zh) 一种itd和gzc机器状态监测方法及装置
CN112697477A (zh) 一种emd设备故障诊断方法及系统
CN112881046A (zh) 一种vmd和gzc机器状态监测方法及装置
CN113435281A (zh) 基于改进hht变换的波纹补偿器故障诊断方法
CN112683395A (zh) 一种elmd和gzc机器状态监测方法及装置
CN112697475A (zh) 一种itd设备故障诊断方法及系统

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
WW01 Invention patent application withdrawn after publication

Application publication date: 20210423