一种LMD多尺度波动分析状态监测方法及装置与流程

一种lmd多尺度波动分析状态监测方法及装置
技术领域
1.本发明涉及设备状态监测与故障诊断领域，具体涉及一种lmd多尺度波动分析状态监测方法及装置。


背景技术：

2.设备振动信号包含丰富的分形特征，这些分形特征能够描述设备的运行状态。盒维数、功率谱分析和重标极差方法可以估计平稳信号的单重分形参数，去趋势波动分析（dfa）能够估计非平稳信号的单重分形维数。然而，设备出现故障时，其振动信号通常是非平稳的，且具有多重分形特征，这时传统的分形维数估计方法会产生比较大的误差。多重分形去趋势波动分析（mfdfa）能够估计非平稳信号的多重分形参数，但是mfdfa方法存在着分析尺度需要人工确定、拟合多项式趋势阶数难以确定和数据段之间不连续的问题。目前，已经有文献提出了基于emd的mfdfa版本（mfdfaemd），用来解决mfdfa存在的问题。然而，mfdfaemd采用的线性滤波方法容易破坏原始信号的分形结构，且存在着负频率现象，这些缺陷严重影响了mfdfaemd的应用效果。综上所述，现有技术难以准确提取设备振动信号的多重分形特征，难以准确检测设备运行状态。


技术实现要素：

3.本发明要解决的问题是针对以上不足，提出一种lmd多尺度波动分析状态监测方法及装置（本发明提出的方法简称为mfdfaolmd）。采用本发明所提出的方法对设备振动信号进行分析，能够有效提取设备振动信号的多重分形特征，克服mfdfa方法存在的分析尺度需要人工确定、拟合多项式趋势阶数难以确定和数据段之间不连续的问题，解决mfdfaemd方法存在的原始信号分形结构破坏和负频率现象，具有分析结果准确度和精确度高，设备运行状态识别结果正确率高等优点。
4.为解决以上技术问题，本发明采取的技术方案如下：一种lmd多尺度波动分析状态监测方法，其特征在于：包括以下步骤：步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2,
ꢀ…
,n，n为采样信号的长度；步骤2：采用局部均值分解（local mean decomposition, lmd）算法将信号x（k）分解成n个分量和一个趋势项之和，即，其中，c
i
（k）代表由lmd算法得到的第i个分量，r
n
（k）代表由lmd算法得到的趋势项，本例中，n=10；步骤3：采用非线性判别算法从lmd分解结果中排除噪声分量和趋势项，保留包含分形特征的分量c
f
（k）, f=1,2,
…
,p，p代表滤波后剩余分量的数量；步骤4：确定c
f
（k）的局部极大值和局部极小值，采用newton插值函数分别对c
f
（k）的局部极大值和局部极小值进行插值，采用最小二乘法分别拟合c
f
（k）的上包络u(k)和下
包络l(k)，则c
f
（k）的包络定义为，符号|x|表示对x取绝对值；步骤5：重复执行公式 m次，j=1,2,
…
,m，直到，得到c
f
（k）的频率调制部分fm
m
(k)，e
j
(k)代表c
j
(k)的包络，c
j
(k)=fm
(j
‑
1)
(k)，c1(k)= c
f
（k）；步骤6：采用teager能量算子（teager energy operator, teo）计算fm
m
(k)的瞬时频率，获得c
f
（k）的瞬时频率instf
f
(k)，得到c
f
（k）的瞬时尺度；步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为；步骤8：将y
s
(k)分成不重叠的n
s
段长度为s的数据，由于数据长度n通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2n
s
段数据；步骤9：计算每段数据的方差：步骤9：计算每段数据的方差：；步骤10：计算q阶函数：；步骤11：改变s的取值，s=s
f
，f=1,2,
…
,p，重复上述步骤3到步骤10，得到关于q和s的方差函数f
q
(s);步骤12：如果x(k)存在分形特征，则f
q
(s)和尺度s之间存在幂律关系：f
q
(s)~s
h(q)
，h(q)代表x(k)的广义hurst指数；当q=0时， h(0)通过下式所定义的对数平均过程来确定：；步骤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：将上述步骤所述方法部署在状态监测装置上，对设备状态进行监测。
5.进一步的，所述步骤2中lmd算法包括以下步骤：1）计算局部均值函数：确定原始信号x(k)所有的局部极值点n
i
，计算相邻两个极值点n
i
和n
i+1
的平均值m
i
，即将所有相邻两个极值点的平均值m
i
用折线连接，然后采用移动平均方法进行平滑处理，得到局部均值函数m
11
(k)；本例中，移动平均方法中的平滑步长设置为5；2）估计信号的包络值：采用局部极值点n
i
计算包络估计值a
i
同样，将所有相邻两个包络估计值a
i
用折线连接，然后采用移动平均方法进行平滑处理，得到包络估计函数a
11
(k)；3）将局部均值函数m
11
(k)从原始信号x（k）中分离出来， 得到4）用h
11
(k) 除以包络估计函数a
11
(k)从而对h
11
(k)进行解调，得到理想地，s
11
(k)是一个纯调频信号，即它的包络估计函数a
12
(k)满足a
12
(k)=1，如果s
11
(k)不满足该条件，则将s
11
(k)作为新数据重复以上迭代过程m次，直到得到一个纯调频信号s
1m
(k)，即s
1m
(k)满足
‑
1≤s
1m
(k) ≤1，它的包络估计函数a
1(m+1)
(k)满足a
1(m+1)
(k)=1，因此，有式中迭代终止的条件为在实际应用中，可以设定一个变动量δ，当满足1
‑
δ≤a
1m
(k) ≤1+δ时，迭代终
止；本例中变动量δ=0.01；5）把迭代过程中产生的所有包络估计函数相乘便可以得到包络信号6）将包络信号a1(k) 和纯调频信号s
1m
(k)相乘便可以得到原始信号的第1个分量；7）将第1个分量c1(k) 从原始信号x(k)中分离出来，得到一个新的信号r1(k)，将r1(k)作为新数据重复以上步骤，循环n次，直到r
n
(k) 为一个单调函数为止至此，将原始信号x(k)分解为n个分量和一个单调函数r
n
(k)之和，即。
6.进一步的，所述步骤3非线性判别算法包括以下步骤：1) 对信号c（k）执行重排操作和替代操作，经重排操作得到的数据用c
shuf
（k）表示，替代操作后得到数据用c
surr
（k）表示；2) 对c（k）、c
shuf
（k）和c
surr
（k）分别执行多重分形去趋势波动分析（multifractal detrended fluctuation analysis, mfdfa），得到广义hurst指数曲线，c（k）的广义hurst指数曲线用h（q）表示；c
shuf
（k）的广义hurst指数曲线用h
shuf
（q）表示；c
surr
（k）的广义hurst指数曲线用h
surr
（q）表示；3) 定义两个参数e1和e2，，，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由lmd算法得到的信号分量。
7.进一步的，所述步骤1)中数据重排操作包括以下步骤：随机打乱分量c（k）的排列顺序。
8.进一步的，所述步骤1)中数据替代操作包括以下步骤：1） 对分量c（k）执行离散傅里叶变换，获得分量c（k）的相位；2） 用一组位于(
‑
π,π)区间内的伪独立同分布数来代替分量c（k）的原始相位；3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据c
ifft
（k），求取数据c
ifft
（k）的实部。
9.进一步的，所述步骤2)中mfdfa方法包括以下步骤：
1）构造x(k)，k=1,2,
…
,n，的轮廓y(i)： ；x(k)代表权利要求3所述步骤2）中的c（k）或c
shuf
（k）或c
surr
（k）；2）将信号轮廓y(i)分成不重叠的n
s
段长度为s的数据，由于数据长度n通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2n
s
段数据；3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：；y
v
(i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（mf
‑
）dfam；本例中，m=1；4） 计算第q阶波动函数的平均值：；5）如果x(k)存在自相似特征，则第q阶波动函数的平均值f
q
(s)和时间尺度s之间存在幂律关系：f
q
(s)~s
h(q)
；当q=0时，步骤4）中的公式发散，这时h(0)通过下式所定义的对数平均过程来确定：；6）对步骤5）中的公式两边取对数可得ln[f
q
(s)]=h(q)ln(s)+c ，c为常数，由此可以获得直线的斜率h(q)。
[0010]
进一步的，所述步骤4中的最小二乘法包括以下步骤：对x（t），t=1, 2,
ꢀ…
,n，x（t）代表步骤4中对c
f
（k）的局部极大值进行插值产生的序列或对c
f
（k）的局部极小值进行插值产生的序列，n代表插值序列的长度；1）事先选定一组函数r
k
(t)，k=1,2,
…
,m，m<n，构造函数f(t)=a1r1(t)+a2r2(t)+
…
+ a
m
r
m
(t)，其中r
k
(t)代表二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线或复合函数曲线；
2）计算最小二乘指标 ；3）令j对的a
k
偏导数 ，k=1,2,
…
,m，这时(a1,a2,
…
,a
m
)
t =(r
t
r)
‑1r
t
x， x=(x(1),x(2),
…
,x(n)) t
， ，其中，r
t
代表r的转置矩阵，r
‑1代表r的逆矩阵；4）r
k
(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标j，选取j最小时所对应的曲线形式为r
k
(t)的形式。
[0011]
进一步的，所述步骤6中teager能量算子方法包括以下步骤：1）对信号c(k)，k=1,2,
…
,n，c(k)= fm
m
(k)，构建函数 ψ(c(k))=c2(k)
‑ꢀ
c(k+1) c(k
‑
1);2）令d(k)=c(k)
‑
c(k
‑
1)，信号c(k)的瞬时频率instf(k)定义为：。
[0012]
基于以上一种lmd多尺度波动分析状态监测方法，实现该方法的装置，所述步骤16中状态监测装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现以上所述算法。
[0013]
各步骤作用：第1)步：利用加速度传感器采集振动信号；第2)步：将原始信号分解成不同分量和的形式，其中有些分量对应噪声和趋势项，有些分量包含分形特征；第3步：利用非线性判别算法去除信号分解结果中的噪声分量和趋势项，只保留包含分形特征的信号分量；第4)~6)步：分离每个分形信号分量的频率调制部分，利用teo估计每个分形信号分量的瞬时频率和瞬时尺度；第7)步：根据分析尺度选择合适的分形信号分量，将选取的分形信号分量求和作为该分析尺度所对应的信号去趋势结果；第8)~14)步：对每个分析尺度所对应的信号去趋势结果执行波动分析，得到原始信号的多重分形谱；第15)步：提取多重分形谱的左端点、右端点和极值点所对应的奇异指数，将这三个参数作为设备运行状态的特征参数；
第16)步：将上述算法部署在设备状态监测装置上，对设备状态进行监测。本发明采用以上技术方案，与现有技术相比，本发明具有以下优点：1) 采用lmd方法分解振动信号，根据非线性滤波方法去除噪声分量和趋势项，能够保护原始信号的分形结构，避免线性滤波方法对原始信号分形结构的破坏；2) 分离信号分量的频率调制部分，利用teo估计信号分量的瞬时频率，能够确保瞬时频率保持正值，避免了负频率现象；3) 根据信号分量的瞬时频率计算相应的瞬时尺度，根据信号分量的瞬时尺度进行波动分析，避免了人工设定尺度的缺陷；4) 利用lmd方法自动确定信号趋势的类型，并保证信号趋势的连续性，有效解决了现有技术的缺陷；5) 分析结果准确度和精确度高，设备运行状态识别结果正确率高。
[0014]
下面结合附图和实施例对本发明做进一步说明。
附图说明
[0015]
附图1为本发明实施例中本发明方法的流程图；附图2为本发明实施例中设备状态监测装置的示意图；附图3为本发明实施例中由多重分形级联模型产生的多重分形仿真信号；附图4为本发明实施例中采用mfdfaemd方法获得的多重分形仿真信号的瞬时频率，信号分量数目为10；附图5为本发明实施例中采用mfdfaolmd方法获得的多重分形仿真信号的瞬时频率，信号分量数目为10；附图6为本发明实施例中分别采用mfdfa、mfdfaemd和mfdfaolmd方法对多重分形仿真信号分析结果的对比图；附图7为本发明实施例中两个非线性判别参数的计算结果，“圆圈”和“方形”符号分别代表e1和e2；附图8为本发明实施例中分别采用mfdfa、mfdfaemd和mfdfaolmd方法对含噪多重分形仿真信号分析结果的对比图；附图9为本发明实施例中由emd得到各个信号分量与原始信号的相关系数；附图10为本发明实施例中分别采用mfdfa、基于相关滤波的mfdfaemd和基于相关滤波的mfdfaolmd方法对含噪多重分形仿真信号分析结果的对比图；附图11为本发明实施例中的四种齿轮箱振动信号，(a)~(d)分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；附图12为本发明实施例中采用mfdfa得到的这四种齿轮箱振动信号的多重分形谱；附图13为本发明实施例中采用mfdfaemd得到的这四种齿轮箱振动信号的多重分形谱；附图14为本发明实施例中采用mfdfaolmd得到的这四种齿轮箱振动信号的多重分形谱；附图15为本发明实施例中由mfdfa得到的多重分形谱的左端点、右端点和极值点
所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；附图16为本发明实施例中由mfdfaemd得到的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；附图17为本发明实施例中由mfdfaolmd得到的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态。
具体实施方式
[0016]
实施例，如图1、图2所示，一种lmd多尺度波动分析状态监测方法，其特征在于：包括以下步骤：步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2,
ꢀ…
,n，n为采样信号的长度；步骤2：采用局部均值分解（local mean decomposition, lmd）算法将信号x（k）分解成n个分量和一个趋势项之和，即，其中，c
i
（k）代表由lmd算法得到的第i个分量，r
n
（k）代表由lmd算法得到的趋势项，本例中，n=10；步骤3：采用非线性判别算法从lmd分解结果中排除噪声分量和趋势项，保留包含分形特征的分量c
f
（k）, f=1,2,
…
,p，p代表滤波后剩余分量的数量；步骤4：确定c
f
（k）的局部极大值和局部极小值，采用newton插值函数分别对c
f
（k）的局部极大值和局部极小值进行插值，采用最小二乘法分别拟合c
f
（k）的上包络u(k)和下包络l(k)，则c
f
（k）的包络定义为，符号|x|表示对x取绝对值；步骤5：重复执行公式 m次，j=1,2,
…
,m，直到，得到c
f
（k）的频率调制部分fm
m
(k)，e
j
(k)代表c
j
(k)的包络，c
j
(k)=fm
(j
‑
1)
(k)，c1(k)= c
f
（k）；步骤6：采用teager能量算子（teager energy operator, teo）计算fm
m
(k)的瞬时频率，获得c
f
（k）的瞬时频率instf
f
(k)，得到c
f
（k）的瞬时尺度；步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为；步骤8：将y
s
(k)分成不重叠的n
s
段长度为s的数据，由于数据长度n通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2n
s
段数据；步骤9：计算每段数据的方差：
；步骤10：计算q阶函数：；步骤11：改变s的取值，s=s
f
，f=1,2,
…
,p，重复上述步骤3到步骤10，得到关于q和s的方差函数f
q
(s);步骤12：如果x(k)存在分形特征，则f
q
(s)和尺度s之间存在幂律关系：f
q
(s)~s
h(q)
，h(q)代表x(k)的广义hurst指数；当q=0时， h(0)通过下式所定义的对数平均过程来确定：；步骤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：将上述步骤所述方法部署在状态监测装置上，对设备状态进行监测。
[0017]
所述步骤2中lmd算法包括以下步骤：1）计算局部均值函数：确定原始信号x(k)所有的局部极值点n
i
，计算相邻两个极值点n
i
和n
i+1
的平均值m
i
，即将所有相邻两个极值点的平均值m
i
用折线连接，然后采用移动平均方法进行平滑处理，得到局部均值函数m
11
(k)；本例中，移动平均方法中的平滑步长设置为5；2）估计信号的包络值：采用局部极值点n
i
计算包络估计值a
i
同样，将所有相邻两个包络估计值a
i
用折线连接，然后采用移动平均方法进行平滑处理，得到包络估计函数a
11
(k)；3）将局部均值函数m
11
(k)从原始信号x（k）中分离出来， 得到
4）用h
11
(k) 除以包络估计函数a
11
(k)从而对h
11
(k)进行解调，得到理想地，s
11
(k)是一个纯调频信号，即它的包络估计函数a
12
(k)满足a
12
(k)=1，如果s
11
(k)不满足该条件，则将s
11
(k)作为新数据重复以上迭代过程m次，直到得到一个纯调频信号s
1m
(k)，即s
1m
(k)满足
‑
1≤s
1m
(k) ≤1，它的包络估计函数a
1(m+1)
(k)满足a
1(m+1)
(k)=1，因此，有式中迭代终止的条件为在实际应用中，可以设定一个变动量δ，当满足1
‑
δ≤a
1m
(k) ≤1+δ时，迭代终止；本例中变动量δ=0.01；5）把迭代过程中产生的所有包络估计函数相乘便可以得到包络信号6）将包络信号a1(k) 和纯调频信号s
1m
(k)相乘便可以得到原始信号的第1个分量；7）将第1个分量c1(k) 从原始信号x(k)中分离出来，得到一个新的信号r1(k)，将r1(k)作为新数据重复以上步骤，循环n次，直到r
n
(k) 为一个单调函数为止
至此，将原始信号x(k)分解为n个分量和一个单调函数r
n
(k)之和，即。
[0018]
所述步骤3非线性判别算法包括以下步骤：1) 对信号c（k）执行重排操作和替代操作，经重排操作得到的数据用c
shuf
（k）表示，替代操作后得到数据用c
surr
（k）表示；2) 对c（k）、c
shuf
（k）和c
surr
（k）分别执行多重分形去趋势波动分析（multifractal detrended fluctuation analysis, mfdfa），得到广义hurst指数曲线，c（k）的广义hurst指数曲线用h（q）表示；c
shuf
（k）的广义hurst指数曲线用h
shuf
（q）表示；c
surr
（k）的广义hurst指数曲线用h
surr
（q）表示；3) 定义两个参数e1和e2，，，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由lmd算法得到的信号分量。
[0019]
所述步骤1)中数据重排操作包括以下步骤：随机打乱分量c（k）的排列顺序。
[0020]
所述步骤1)中数据替代操作包括以下步骤：1） 对分量c（k）执行离散傅里叶变换，获得分量c（k）的相位；2） 用一组位于(
‑
π,π)区间内的伪独立同分布数来代替分量c（k）的原始相位；3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据c
ifft
（k），求取数据c
ifft
（k）的实部。
[0021]
所述步骤2)中mfdfa方法包括以下步骤：1）构造x(k)，k=1,2,
…
,n，的轮廓y(i)：,n，的轮廓y(i)：；x(k)代表权利要求3所述步骤2）中的c（k）或c
shuf
（k）或c
surr
（k）；2）将信号轮廓y(i)分成不重叠的n
s
段长度为s的数据，由于数据长度n通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2n
s
段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：；y
v
(i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（mf
‑
）dfam；本例中，m=1；4） 计算第q阶波动函数的平均值：；5）如果x(k)存在自相似特征，则第q阶波动函数的平均值f
q
(s)和时间尺度s之间存在幂律关系：f
q
(s)~s
h(q)
；当q=0时，步骤4）中的公式发散，这时h(0)通过下式所定义的对数平均过程来确定：；6）对步骤5）中的公式两边取对数可得ln[f
q
(s)]=h(q)ln(s)+c ，c为常数，由此可以获得直线的斜率h(q)。
[0022]
所述步骤4中的最小二乘法包括以下步骤：对x（t），t=1, 2,
ꢀ…
,n，x（t）代表步骤4中对c
f
（k）的局部极大值进行插值产生的序列或对c
f
（k）的局部极小值进行插值产生的序列，n代表插值序列的长度；1）事先选定一组函数r
k
(t)，k=1,2,
…
,m，m<n，构造函数f(t)=a1r1(t)+a2r2(t)+
…
+ a
m
r
m
(t)，其中r
k
(t)代表二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线或复合函数曲线；2）计算最小二乘指标 ；3）令j对的a
k
偏导数 ，k=1,2,
…
,m，这时(a1,a2,
…
,a
m
)
t =(r
t
r)
‑1r
t
x， x=(x(1),x(2),
…
,x(n)) t
， ，其中，r
t
代表r的转置矩阵，r
‑1代表r的逆矩阵；4）r
k
(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标j，选取j最小时所对应的
曲线形式为r
k
(t)的形式。
[0023]
所述步骤6中teager能量算子方法包括以下步骤：1）对信号c(k)，k=1,2,
…
,n，c(k)= fm
m
(k)，构建函数 ψ(c(k))=c2(k)
‑ꢀ
c(k+1) c(k
‑
1);2）令d(k)=c(k)
‑
c(k
‑
1)，信号c(k)的瞬时频率instf(k)定义为：。
[0024]
基于以上一种lmd多尺度波动分析状态监测方法，实现该方法的装置，所述步骤16中状态监测装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现以上所述算法。
[0025]
实验1 采用由多重分形级联模型产生的多重分形仿真信号对本发明所述算法的性能进行验证。
[0026]
首先采用多重分形级联模型产生的多重分形仿真信号对mfdfa、mfdfaemd和mfdfaolmd的性能进行验证。本例中，p=0.375 ，n=14，产生的多重分形仿真信号如图3所示。采用mfdfaemd计算多重分形仿真信号的瞬时频率，结果如图4所示。从图4可以看出，由mfdfaemd计算得到的瞬时频率存在着许多负频率，由于负频率没有物理意义，因此mfdfaemd分析结果存在着较大的误差。采用mfdfaolmd计算多重分形仿真信号的瞬时频率，结果如图5所示。从图5可以看出，由mfdfaolmd计算得到的瞬时频率全部是正频率，因此mfdfaolmd分析结果符合实际情况。接下来，分别采用mfdfa、mfdfaemd和mfdfaolmd计算多重分形仿真信号的多重分形谱，结果如图6所示。根据图6所示结果，经计算可知， 由mfdfa得到的多重分形谱与理论值的绝对误差平均值为0.064，相对误差平均值为12.21%，由mfdfaemd得到的多重分形谱与理论值的绝对误差平均值为0.035，相对误差平均值为6.57%，由mfdfaolmd得到的多重分形谱与理论值的绝对误差平均值为0.023，相对误差平均值为5.13%，因此由mfdfaolmd得到的多重分形谱比由mfdfa得到的多重分形谱绝对误差平均值减小64.06%，相对误差平均值减小57.99%，由mfdfaolmd得到的多重分形谱比由mfdfaemd得到的多重分形谱绝对误差平均值减小34.29%，相对误差平均值减小21.92%。图7所示为本发明实施例中两个非线性判别参数的计算结果，可以看出所有的信号分量都包含分形信号分量。然后，通过向上述多重分形仿真信号添加高斯白噪声的方法构造信噪比为20db的含噪信号。分别采用mfdfa、mfdfaemd和mfdfaolmd计算含噪多重分形仿真信号的多重分形谱，结果如图8所示。根据图8所示结果，由mfdfa获得的多重分形谱完全偏离了理论值，由mfdfaemd得到的多重分形谱与理论值的绝对误差平均值为0.084，相对误差平均值为11.81%，由mfdfaolmd得到的多重分形谱与理论值的绝对误差平均值为0.031，相对误差平均值为5.36%，因此由mfdfaolmd得到的多重分形谱比由mfdfaemd得到的多重分形谱绝对误差平均值减小63.10%，相对误差平均值减小54.61%。由此可见，与mfdfa和mfdfaemd相比，
mfdfaolmd具有更好地抗噪性。图9所示为本发明实施例中由emd得到各个信号分量与原始信号的相关系数，可以看出第7个分量与原信号具有最弱的相关性，应该从原信号中去除。图10所示为本发明实施例中分别采用mfdfa、基于相关滤波的mfdfaemd和基于相关滤波mfdfaolmd方法对含噪多重分形仿真信号分析结果的对比图。从图10可以看出，基于相关滤波的mfdfaemd和基于相关滤波的mfdfaolmd方法对含噪多重分形仿真信号分析结果完全偏离了理论值，因此相关滤波方法容易破坏原始信号的分形结构。
[0027]
实验2 采用齿轮箱实验信号对本发明所述算法的性能进行验证。
[0028]
本发明所用齿轮箱振动数据一个齿轮箱故障模拟实验台。这个实验通过在一个轮齿的齿根上制造不同程度的划痕，直到最后完全破坏掉这个轮齿，来模拟轮齿从正常状态到失效的过程。实验中使用的齿轮箱为二级齿轮传动，从输入端到输出端的齿轮齿数分别为25、40、22和55，故障轮齿位于输入轴齿轮上，驱动电机转速为2000rmp，齿轮箱振动信号由位于输入端外壳上的加速度传感器测取。采集到的振动信号包含着正常、轻度划痕、重度划痕和断齿四种故障状态，这些振动信号在一定程度上代表了轮齿从正常到失效的过程。振动信号采样频率为16384hz，在每种齿轮箱状态下采集长度为10000点的20段数据，，这四种齿轮箱振动信号如图11所示。首先采用mfdfa方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图12所示，可以看出轻度划痕和重度划痕所对应的多重分形谱严重重叠。然后采用mfdfaemd方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图13所示，可以看出这四种齿轮箱状态所对应的多重分形谱严重重叠。最后采用mfdfaolmd方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图14所示，可以看出断齿状态振动信号的多重分形谱与其他三种齿轮箱状态所对应的多重分形谱明显不同，而正常、轻度划痕和重度划痕齿轮状态所对应的多重分形谱能够在α<0.4时清楚分离。分别提取由mfdfa、mfdfaemd和mfdfaolmd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态进行分类，结果分别如图15~17所示。从图15可以看出，利用由mfdfa方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分正常状态和断齿状态，但是不能区分轻度划痕和重度划痕状态，因此齿轮箱状态识别率为50%。从图16可以看出，利用由mfdfaemd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分正常状态和断齿状态，但是不能区分轻度划痕和重度划痕状态，因此齿轮箱状态识别率为50%。从图17可以看出，利用由mfdfaolmd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分这四种齿轮箱状态，因此齿轮箱状态识别率为100%。可以看出，采用mfdfaolmd方法可以将齿轮箱状态识别正确率提高50%。
[0029]
根据试验结果，分析后认为：1) 采用lmd方法自适应分解振动信号，根据非线性滤波方法去除噪声分量和趋势项，能够保护原始信号的分形结构，避免线性滤波方法对原始信号分形结构的破坏；2) 分离信号分量的频率调制部分，利用teo估计信号分量的瞬时频率，能够确保瞬时频率保持正值，避免了负频率现象；3) 根据信号分量的瞬时频率计算相应的瞬时尺度，根据信号分量的瞬时尺度进行波动分析，避免了人工设定尺度的缺陷；
4) 利用lmd方法自动确定信号趋势的类型，并保证信号趋势的连续性，有效解决了现有技术的缺陷；5) 分析结果准确度和精确度高，设备运行状态识别结果正确率高。
[0030]
本领域技术人员应该认识到，上述的具体实施方式只是示例性的，是为了使本领域技术人员能够更好的理解本发明内容，不应理解为是对本发明保护范围的限制，只要是根据本发明技术方案所作的改进，均落入本发明的保护范围。