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

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

Info

Publication number
CN112697479A
CN112697479A CN202011240502.3A CN202011240502A CN112697479A CN 112697479 A CN112697479 A CN 112697479A CN 202011240502 A CN202011240502 A CN 202011240502A CN 112697479 A CN112697479 A CN 112697479A
Authority
CN
China
Prior art keywords
signal
data
itd
fractal
scale
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
CN202011240502.3A
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 CN202011240502.3A priority Critical patent/CN112697479A/zh
Publication of CN112697479A publication Critical patent/CN112697479A/zh
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M99/00—Subject matter not provided for in other groups of this subclass
    • G01M99/004—Testing the effects of speed or acceleration
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/02—Gearings; Transmission mechanisms
    • G01M13/021—Gearings
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/02—Gearings; Transmission mechanisms
    • G01M13/028—Acoustic or vibration analysis
    • 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
    • G06F17/15—Correlation function computation including computation of convolution operations

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Acoustics & Sound (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

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

Description

一种ITD多尺度波动分析状态监测方法及装置
技术领域
本发明涉及设备状态监测与故障诊断领域，具体涉及一种ITD多尺度波动分析状态监测方法及装置。
背景技术
设备振动信号包含丰富的分形特征，这些分形特征能够描述设备的运行状态。盒维数、功率谱分析和重标极差方法可以估计平稳信号的单重分形参数，去趋势波动分析（DFA）能够估计非平稳信号的单重分形维数。然而，设备出现故障时，其振动信号通常是非平稳的，且具有多重分形特征，这时传统的分形维数估计方法会产生比较大的误差。多重分形去趋势波动分析（MFDFA）能够估计非平稳信号的多重分形参数，但是MFDFA方法存在着分析尺度需要人工确定、拟合多项式趋势阶数难以确定和数据段之间不连续的问题。目前，已经有文献提出了基于EMD的MFDFA版本（MFDFAemd），用来解决MFDFA存在的问题。然而，MFDFAemd采用的线性滤波方法容易破坏原始信号的分形结构，且存在着负频率现象，这些缺陷严重影响了MFDFAemd的应用效果。综上所述，现有技术难以准确提取设备振动信号的多重分形特征，难以准确检测设备运行状态。
发明内容
本发明要解决的问题是针对以上不足，提出一种ITD多尺度波动分析（本发明提出的方法简称为MFDFAoitd）状态监测方法及装置。采用本发明所提出的方法对设备振动信号进行分析，能够有效提取设备振动信号的多重分形特征，克服MFDFA方法存在的分析尺度需要人工确定、拟合多项式趋势阶数难以确定和数据段之间不连续的问题，解决MFDFAemd方法存在的原始信号分形结构破坏和负频率现象，具有分析结果准确度和精确度高，设备运行状态识别结果正确率高等优点。
为解决以上技术问题，本发明采取的技术方案如下：一种ITD多尺度波动分析状态监测方法，其特征在于;包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用本征时间尺度分解（Intrinsic Time-scale Decomposition, ITD）算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure DEST_PATH_IMAGE002
，其中，ci（k）代表由ITD算法得到的第i个分量，rn（k）代表由ITD算法得到的趋势项，本例中，n=10；
步骤3：采用非线性判别算法从ITD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量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，本例中，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中ITD算法包括以下步骤：
1) 对于信号xt， t=1, 2, …,N，xt代表权利要求1所述步骤2中的信号x(k)，定义一个算子
Figure DEST_PATH_IMAGE022
用于抽取低频基线信号，即：
Figure DEST_PATH_IMAGE024
其中
Figure DEST_PATH_IMAGE026
是基线信号，
Figure DEST_PATH_IMAGE028
是一个固有旋转分量, 假设
Figure DEST_PATH_IMAGE030
是一个实值信号，
Figure DEST_PATH_IMAGE032
代表xt的局部极值所对应的时刻，为方便起见定义
Figure DEST_PATH_IMAGE034
，如果xt在某个区间上具有恒定值，考虑到邻近的信号存在着波动，我们仍然认为xt在这个区间上包含着极值，这时
Figure DEST_PATH_IMAGE036
是该区间的右端点，为方便起见，定义
Figure DEST_PATH_IMAGE038
，
Figure DEST_PATH_IMAGE040
; 假设
Figure DEST_PATH_IMAGE042
和
Figure DEST_PATH_IMAGE044
在
Figure DEST_PATH_IMAGE046
上有定义，xk在
Figure DEST_PATH_IMAGE048
上有定义，在区间
Figure DEST_PATH_IMAGE050
上的连续极值点之间定义一个分段线性的基线信号抽取算子
Figure 781460DEST_PATH_IMAGE022
，即：
Figure DEST_PATH_IMAGE052
其中
Figure DEST_PATH_IMAGE054
这里参数
Figure DEST_PATH_IMAGE056
是一个线性增益，
Figure DEST_PATH_IMAGE058
，本例中
Figure DEST_PATH_IMAGE060
；
2) 定义一个固有旋转分量抽取算子
Figure DEST_PATH_IMAGE062
，即：
Figure DEST_PATH_IMAGE064
其中，
Figure DEST_PATH_IMAGE066
为第1次迭代分解得到的分量；
Figure DEST_PATH_IMAGE068
为第1次迭代分解得到的基线信号；在第1次迭代中，xt代表权利要求1所述步骤2中x(k)；
3）再把
Figure 885551DEST_PATH_IMAGE068
作为新数据，重复上述步骤，可以分离出频率依次降低的固有旋转分量
Figure DEST_PATH_IMAGE070
，直到基线信号变得单调为止，这样xt的整个分解过程可以写为：
Figure DEST_PATH_IMAGE072
其中
Figure DEST_PATH_IMAGE074
代表第i次迭代分解得到的分量，
Figure DEST_PATH_IMAGE076
代表第i次迭代分解得到的基线信号。
进一步的，所述步骤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_IMAGE078
，
Figure DEST_PATH_IMAGE080
，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由ITD算法得到的信号分量。
进一步的，所述步骤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 DEST_PATH_IMAGE082
Figure DEST_PATH_IMAGE084
；
x(k)代表权利要求3所述步骤1）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；
为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure DEST_PATH_IMAGE086
Figure DEST_PATH_IMAGE088
；
y v (i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；本例中，m=1；
4） 计算第q阶波动函数的平均值：
Figure DEST_PATH_IMAGE090
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值F q (s)和时间尺度s之间存在幂律关系：
F q (s)~s H(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 374212DEST_PATH_IMAGE020
；
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 DEST_PATH_IMAGE092
；
3）令J对的ak偏导数
Figure DEST_PATH_IMAGE094
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure DEST_PATH_IMAGE096
，其中，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 DEST_PATH_IMAGE098
。
基于以上所述的一种ITD多尺度波动分析状态监测方法，实现该方法的装置，所述装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现以上所述算法。
各步骤作用：
第1)步：采集振动信号；
第2)步：将原始信号分解成不同分量和的形式，其中有些分量对应噪声和趋势项，有些分量包含分形特征；
第3步：利用非线性判别算法去除信号分解结果中的噪声分量和趋势项，只保留包含分形特征的信号分量；
第4)~6)步：分离每个分形信号分量的频率调制部分，利用TEO估计每个分形信号分量的瞬时频率和瞬时尺度；
第7)步：根据分析尺度选择合适的分形信号分量，将选取的分形信号分量求和作为该分析尺度所对应的信号去趋势结果；
第8)~14)步：对每个分析尺度所对应的信号去趋势结果执行波动分析，得到原始信号的多重分形谱；
第15)步：提取多重分形谱的左端点、右端点和极值点所对应的奇异指数，将这三个参数作为设备运行状态的特征参数；
第16)步：将上述算法部署在设备状态监测装置上，对设备状态进行监测。
本发明采用以上技术方案，与现有技术相比，本发明具有以下优点：
1) 采用ITD方法分解振动信号，根据非线性滤波方法去除噪声分量和趋势项，能够保护原始信号的分形结构，避免线性滤波方法对原始信号分形结构的破坏；
2) 分离信号分量的频率调制部分，利用TEO估计信号分量的瞬时频率，能够确保瞬时频率保持正值，避免了负频率现象；
3) 根据信号分量的瞬时频率计算相应的瞬时尺度，根据信号分量的瞬时尺度进行波动分析，避免了人工设定尺度的缺陷；
4) 利用ITD方法自动确定信号趋势的类型，并保证信号趋势的连续性，有效解决了现有技术的缺陷；
5) 分析结果准确度和精确度高，设备运行状态识别结果正确率高。
下面结合附图和实施例对本发明做进一步说明。
附图说明
附图1为本发明实施例中本发明方法的流程图；
附图2为本发明实施例中设备状态监测装置的示意图；
附图3为本发明实施例中由多重分形级联模型产生的多重分形仿真信号；
附图4为本发明实施例中采用MFDFAemd方法获得的多重分形仿真信号的瞬时频率，
信号分量数目为10；
附图5为本发明实施例中采用MFDFAoitd方法获得的多重分形仿真信号的瞬时频率，
信号分量数目为10；
附图6为本发明实施例中分别采用MFDFA、MFDFAemd和MFDFAoitd方法对多重分形仿真信号分析结果的对比图；
附图7为本发明实施例中两个非线性判别参数的计算结果，“圆圈”和“方形”符号分别代表e1和e2；
附图8为本发明实施例中分别采用MFDFA、MFDFAemd和MFDFAoitd方法对含噪多重分形仿真信号分析结果的对比图；
附图9为本发明实施例中由EMD得到各个信号分量与原始信号的相关系数；
附图10为本发明实施例中分别采用MFDFA、基于相关滤波的MFDFAemd和基于相关滤波MFDFAoitd的方法对含噪多重分形仿真信号分析结果的对比图；
附图11为本发明实施例中的四种齿轮箱振动信号，(a)~(d)分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；
附图12为本发明实施例中采用MFDFA得到的这四种齿轮箱振动信号的多重分形谱；
附图13为本发明实施例中采用MFDFAemd得到的这四种齿轮箱振动信号的多重分形谱；
附图14为本发明实施例中采用MFDFAoitd得到的这四种齿轮箱振动信号的多重分形谱；
附图15为本发明实施例中由MFDFA得到的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；
附图16为本发明实施例中由MFDFAemd得到的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态；
附图17为本发明实施例中由MFDFAoitd得到的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态的分类结果，“圆圈”、“方形”、“加号”和“菱形”符号分别代表正常、轻度划痕、重度划痕和断齿齿轮状态。
具体实施方式
实施例，如图1、图2所示，一种ITD多尺度波动分析状态监测方法，其特征在于；包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用本征时间尺度分解（Intrinsic Time-scale Decomposition, ITD）算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure 346585DEST_PATH_IMAGE002
，其中，ci（k）代表由ITD算法得到的第i个分量，rn（k）代表由ITD算法得到的趋势项，本例中，n=10；
步骤3：采用非线性判别算法从ITD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量cf（k）, f=1,2,…,p，p代表滤波后剩余分量的数量；
步骤4：确定cf（k）的局部极大值和局部极小值，采用有理样条插值函数分别对cf（k）的局部极大值和局部极小值进行插值，采用最小二乘法分别拟合cf（k）的上包络u(k)和下包络l(k)，则cf（k）的包络定义为
Figure DEST_PATH_IMAGE099
，符号|x|表示对x取绝对值；
步骤5：重复执行公式
Figure 410224DEST_PATH_IMAGE006
m次，j=1,2,…,m，直到
Figure 295004DEST_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_IMAGE100
；
步骤7：当尺度为s时，则振动信号x(k)的去趋势结果为
Figure 375086DEST_PATH_IMAGE012
；
步骤8：将Y s (k)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
步骤9：计算每段数据的方差：
Figure 573987DEST_PATH_IMAGE014
Figure DEST_PATH_IMAGE101
；
步骤10：计算q阶函数：
Figure 305182DEST_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 610130DEST_PATH_IMAGE020
；
步骤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中ITD算法包括以下步骤：
1) 对于信号xt， t=1, 2, …,N，xt代表权利要求1所述步骤2中的信号x(k)，定义一个算子
Figure 364460DEST_PATH_IMAGE022
用于抽取低频基线信号，即：
Figure DEST_PATH_IMAGE102
其中
Figure 82230DEST_PATH_IMAGE026
是基线信号，
Figure 917200DEST_PATH_IMAGE028
是一个固有旋转分量, 假设
Figure 550307DEST_PATH_IMAGE030
是一个实值信号，
Figure 385408DEST_PATH_IMAGE032
代表xt的局部极值所对应的时刻，为方便起见定义
Figure DEST_PATH_IMAGE103
，如果xt在某个区间上具有恒定值，考虑到邻近的信号存在着波动，我们仍然认为xt在这个区间上包含着极值，这时
Figure 270318DEST_PATH_IMAGE036
是该区间的右端点，为方便起见，定义
Figure 490953DEST_PATH_IMAGE038
，
Figure 154015DEST_PATH_IMAGE040
; 假设
Figure 86199DEST_PATH_IMAGE042
和
Figure 509221DEST_PATH_IMAGE044
在
Figure 476040DEST_PATH_IMAGE046
上有定义，xk在
Figure DEST_PATH_IMAGE104
上有定义，在区间
Figure 310004DEST_PATH_IMAGE050
上的连续极值点之间定义一个分段线性的基线信号抽取算子
Figure 837806DEST_PATH_IMAGE022
，即：
Figure 517049DEST_PATH_IMAGE052
其中
Figure 479320DEST_PATH_IMAGE054
这里参数
Figure 890710DEST_PATH_IMAGE056
是一个线性增益，
Figure 656541DEST_PATH_IMAGE058
，本例中
Figure DEST_PATH_IMAGE105
；
2) 定义一个固有旋转分量抽取算子
Figure 654321DEST_PATH_IMAGE062
，即：
Figure DEST_PATH_IMAGE106
其中，
Figure 189208DEST_PATH_IMAGE066
为第1次迭代分解得到的分量；
Figure 505920DEST_PATH_IMAGE068
为第1次迭代分解得到的基线信号；在第1次迭代中，xt代表权利要求1所述步骤2中x(k)；
3）再把
Figure 775358DEST_PATH_IMAGE068
作为新数据，重复上述步骤，可以分离出频率依次降低的固有旋转分量
Figure 468508DEST_PATH_IMAGE070
，直到基线信号变得单调为止，这样xt的整个分解过程可以写为：
Figure DEST_PATH_IMAGE107
其中
Figure DEST_PATH_IMAGE108
代表第i次迭代分解得到的分量，
Figure DEST_PATH_IMAGE109
代表第i次迭代分解得到的基线信号。
进一步的，所述步骤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_IMAGE110
，
Figure DEST_PATH_IMAGE111
，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由ITD算法得到的信号分量。
进一步的，所述步骤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 966223DEST_PATH_IMAGE082
Figure 594782DEST_PATH_IMAGE084
；
x(k)代表权利要求3所述步骤1）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；
为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure 476150DEST_PATH_IMAGE086
Figure 97624DEST_PATH_IMAGE088
；
y v (i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；本例中，m=1；
4） 计算第q阶波动函数的平均值：
Figure 748048DEST_PATH_IMAGE090
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值F q (s)和时间尺度s之间存在幂律关系：
F q (s)~s H(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 780464DEST_PATH_IMAGE020
；
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 149128DEST_PATH_IMAGE092
；
3）令J对的ak偏导数
Figure 308714DEST_PATH_IMAGE094
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure 79224DEST_PATH_IMAGE096
，其中，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 DEST_PATH_IMAGE112
。
基于以上所述的一种ITD多尺度波动分析状态监测方法，实现该方法的装置，所述装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现以上所述算法。
实验1 采用由多重分形级联模型产生的多重分形仿真信号对本发明所述算法的性能进行验证。
首先采用多重分形级联模型
Figure DEST_PATH_IMAGE114
产生的多重分形仿真信号对MFDFA、MFDFAemd和MFDFAoitd的性能进行验证。本例中，p=0.375 ，n=14，产生的多重分形仿真信号如图3所示。采用MFDFAemd计算多重分形仿真信号的瞬时频率，结果如图4所示。从图4可以看出，由MFDFAemd计算得到的瞬时频率存在着许多负频率，由于负频率没有物理意义，因此MFDFAemd分析结果存在着较大的误差。采用MFDFAoitd计算多重分形仿真信号的瞬时频率，结果如图5所示。从图5可以看出，由MFDFAoitd计算得到的瞬时频率全部是正频率，因此MFDFAoitd分析结果符合实际情况。接下来，分别采用MFDFA、MFDFAemd和MFDFAoitd计算多重分形仿真信号的多重分形谱，结果如图6所示。根据图6所示结果，经计算可知， 由MFDFA得到的多重分形谱与理论值的绝对误差平均值为0.064，相对误差平均值为12.21%，由MFDFAemd得到的多重分形谱与理论值的绝对误差平均值为0.035，相对误差平均值为6.57%，由MFDFAoitd得到的多重分形谱与理论值的绝对误差平均值为0.024，相对误差平均值为5.14%，因此由MFDFAoitd得到的多重分形谱比由MFDFA得到的多重分形谱绝对误差平均值减小62.50%，相对误差平均值减小57.90%，由MFDFAoitd得到的多重分形谱比由MFDFAemd得到的多重分形谱绝对误差平均值减小31.43%，相对误差平均值减小21.77%。图7所示为本发明实施例中两个非线性判别参数的计算结果，可以看出所有的信号分量都包含分形信号分量。然后，通过向上述多重分形仿真信号添加高斯白噪声的方法构造信噪比为20dB的含噪信号。分别采用MFDFA、MFDFAemd和MFDFAoitd计算含噪多重分形仿真信号的多重分形谱，结果如图8所示。根据图8所示结果，由MFDFA获得的多重分形谱完全偏离了理论值，由MFDFAemd得到的多重分形谱与理论值的绝对误差平均值为0.084，相对误差平均值为11.81%，由MFDFAoitd得到的多重分形谱与理论值的绝对误差平均值为0.032，相对误差平均值为5.37%，因此由MFDFAoitd得到的多重分形谱比由MFDFAemd得到的多重分形谱绝对误差平均值减小61.90%，相对误差平均值减小54.53%。由此可见，与MFDFA和MFDFAemd相比，MFDFAoitd具有更好地抗噪性。图9所示为本发明实施例中由EMD得到各个信号分量与原始信号的相关系数，可以看出第7个分量与原信号具有最弱的相关性，应该从原信号中去除。图10所示为本发明实施例中分别采用MFDFA、基于相关滤波的MFDFAemd和基于相关滤波MFDFAoitd的方法对含噪多重分形仿真信号分析结果的对比图。从图10可以看出，基于相关滤波的MFDFAemd和基于相关滤波的MFDFAoitd方法对含噪多重分形仿真信号分析结果完全偏离了理论值，因此相关滤波方法容易破坏原始信号的分形结构。
实验2 采用齿轮箱实验信号对本发明所述算法的性能进行验证。
本发明所用齿轮箱振动数据一个齿轮箱故障模拟实验台。这个实验通过在一个轮齿的齿根上制造不同程度的划痕，直到最后完全破坏掉这个轮齿，来模拟轮齿从正常状态到失效的过程。实验中使用的齿轮箱为二级齿轮传动，从输入端到输出端的齿轮齿数分别为25、40、22和55，故障轮齿位于输入轴齿轮上，驱动电机转速为2000RMP，齿轮箱振动信号由位于输入端外壳上的加速度传感器测取。采集到的振动信号包含着正常、轻度划痕、重度划痕和断齿四种故障状态，这些振动信号在一定程度上代表了轮齿从正常到失效的过程。振动信号采样频率为16384Hz，在每种齿轮箱状态下采集长度为10000点的20段数据，，这四种齿轮箱振动信号如图11所示。首先采用MFDFA方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图12所示，可以看出轻度划痕和重度划痕所对应的多重分形谱严重重叠。然后采用MFDFAemd方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图13所示，可以看出这四种齿轮箱状态所对应的多重分形谱严重重叠。最后采用MFDFAoitd方法对这四种齿轮箱振动信号进行分析，得到这四种齿轮箱振动信号所对应的多重分形谱如图14所示，可以看出断齿状态振动信号的多重分形谱与其他三种齿轮箱状态所对应的多重分形谱明显不同，而正常、轻度划痕和重度划痕齿轮状态所对应的多重分形谱能够在α<0.4时清楚分离。分别提取由MFDFA、MFDFAemd和MFDFAoitd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数对这四种齿轮箱状态进行分类，结果分别如图15~17所示。从图15可以看出，利用由MFDFA方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分正常状态和断齿状态，但是不能区分轻度划痕和重度划痕状态，因此齿轮箱状态识别率为50%。从图16可以看出，利用由MFDFAemd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分正常状态和断齿状态，但是不能区分轻度划痕和重度划痕状态，因此齿轮箱状态识别率为50%。从图17可以看出，利用由MFDFAoitd方法获得的多重分形谱的左端点、右端点和极值点所对应的奇异指数，可以正确区分这四种齿轮箱状态，因此齿轮箱状态识别率为100%。可以看出，采用MFDFAoitd方法可以将齿轮箱状态识别正确率提高50%。
根据试验结果，分析后认为：
1) 采用ITD方法自适应分解振动信号，根据非线性滤波方法去除噪声分量和趋势项，能够保护原始信号的分形结构，避免线性滤波方法对原始信号分形结构的破坏；
2) 分离信号分量的频率调制部分，利用TEO估计信号分量的瞬时频率，能够确保瞬时频率保持正值，避免了负频率现象；
3) 根据信号分量的瞬时频率计算相应的瞬时尺度，根据信号分量的瞬时尺度进行波动分析，避免了人工设定尺度的缺陷；
4) 利用ITD方法自动确定信号趋势的类型，并保证信号趋势的连续性，有效解决了现有技术的缺陷；
5) 分析结果准确度和精确度高，设备运行状态识别结果正确率高。
本领域技术人员应该认识到，上述的具体实施方式只是示例性的，是为了使本领域技术人员能够更好的理解本发明内容，不应理解为是对本发明保护范围的限制，只要是根据本发明技术方案所作的改进，均落入本发明的保护范围。

Claims (9)

1.一种ITD多尺度波动分析状态监测方法，其特征在于；包括以下步骤：
步骤1：利用加速度传感器以采样频率fs测取设备振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；
步骤2：采用本征时间尺度分解（Intrinsic Time-scale Decomposition, ITD）算法将信号x（k）分解成n个分量和一个趋势项之和，即
Figure DEST_PATH_IMAGE001
，其中，ci（k）代表由ITD算法得到的第i个分量，rn（k）代表由ITD算法得到的趋势项；
步骤3：采用非线性判别算法从ITD分解结果中排除噪声分量和趋势项，保留包含分形特征的分量cf（k）, f=1,2,…,p，p代表滤波后剩余分量的数量；
步骤4：确定cf（k）的局部极大值和局部极小值，采用有理样条插值函数分别对cf（k）的局部极大值和局部极小值进行插值，采用最小二乘法分别拟合cf（k）的上包络u(k)和下包络l(k)，则cf（k）的包络定义为
Figure 900906DEST_PATH_IMAGE002
，符号|x|表示对x取绝对值；
步骤5：重复执行公式
Figure DEST_PATH_IMAGE003
m次，j=1,2,…,m，直到
Figure 886222DEST_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 728277DEST_PATH_IMAGE006
；
步骤8：将Y s (k)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
步骤9：计算每段数据的方差：
Figure DEST_PATH_IMAGE007
Figure 185803DEST_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)存在分形特征，则F q (s)和尺度s之间存在幂律关系：F q (s)~s H(q)，H(q)代表x(k)的广义Hurst指数；
当q=0时， H(0)通过下式所定义的对数平均过程来确定：
Figure 62492DEST_PATH_IMAGE010
；
步骤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所述的一种ITD多尺度波动分析状态监测方法，其特征在于；所述步骤2中ITD算法包括以下步骤：
1) 对于信号xt， t=1, 2, …,N，xt代表权利要求1所述步骤2中的信号x(k)，定义一个算子
Figure DEST_PATH_IMAGE011
用于抽取低频基线信号，即：
Figure 947271DEST_PATH_IMAGE012
其中
Figure DEST_PATH_IMAGE013
是基线信号，
Figure 276621DEST_PATH_IMAGE014
是一个固有旋转分量, 假设
Figure DEST_PATH_IMAGE015
是一个实值信号，
Figure 98691DEST_PATH_IMAGE016
代表xt的局部极值所对应的时刻，为方便起见定义
Figure DEST_PATH_IMAGE017
，如果xt在某个区间上具有恒定值，考虑到邻近的信号存在着波动，我们仍然认为xt在这个区间上包含着极值，这时
Figure 847464DEST_PATH_IMAGE018
是该区间的右端点，为方便起见，定义
Figure DEST_PATH_IMAGE019
，
Figure 903145DEST_PATH_IMAGE020
; 假设
Figure DEST_PATH_IMAGE021
和
Figure 719791DEST_PATH_IMAGE022
在
Figure DEST_PATH_IMAGE023
上有定义，xk在
Figure 581437DEST_PATH_IMAGE024
上有定义，在区间
Figure DEST_PATH_IMAGE025
上的连续极值点之间定义一个分段线性的基线信号抽取算子
Figure 167139DEST_PATH_IMAGE011
，即：
Figure 659300DEST_PATH_IMAGE026
其中
Figure DEST_PATH_IMAGE027
这里参数
Figure 963243DEST_PATH_IMAGE028
是一个线性增益，
Figure DEST_PATH_IMAGE029
；
2) 定义一个固有旋转分量抽取算子
Figure 425317DEST_PATH_IMAGE030
，即：
Figure DEST_PATH_IMAGE031
其中，
Figure 131105DEST_PATH_IMAGE032
为第1次迭代分解得到的分量；
Figure DEST_PATH_IMAGE033
为第1次迭代分解得到的基线信号；在第1次迭代中，xt代表权利要求1所述步骤2中x(k)；
3）再把
Figure 528588DEST_PATH_IMAGE033
作为新数据，重复上述步骤，可以分离出频率依次降低的固有旋转分量
Figure 319826DEST_PATH_IMAGE034
，直到基线信号变得单调为止，这样xt的整个分解过程可以写为：
Figure DEST_PATH_IMAGE035
其中
Figure 460958DEST_PATH_IMAGE036
代表第i次迭代分解得到的分量，
Figure DEST_PATH_IMAGE037
代表第i次迭代分解得到的基线信号。
3.根据权利要求1所述的一种ITD多尺度波动分析状态监测方法，其特征在于：所述步骤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 21252DEST_PATH_IMAGE038
，
Figure DEST_PATH_IMAGE039
，如果e1和e2都小于10%，则信号c（k）被判别为噪声分量或趋势项，c（k）代表由ITD算法得到的信号分量。
4.根据权利要求3所述的一种ITD多尺度波动分析状态监测方法及装置，其特征在于，所述步骤1)中数据重排操作包括以下步骤：随机打乱分量c（k）的排列顺序。
5.根据权利要求3所述的一种ITD多尺度波动分析状态监测方法及装置，其特征在于：所述步骤1)中数据替代操作包括以下步骤：
1） 对分量c（k）执行离散傅里叶变换，获得分量c（k）的相位；
2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量c（k）的原始相位；
3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据cIFFT（k），求取数据cIFFT（k）的实部。
6.根据权利要求3所述的一种ITD多尺度波动分析状态监测方法，其特征在于：所述步骤2)中MFDFA方法包括以下步骤：
1）构造x(k)，k=1,2,…,N，的轮廓Y(i)：
Figure 861075DEST_PATH_IMAGE040
Figure DEST_PATH_IMAGE041
；
x(k)代表权利要求3所述步骤1）中的c（k）或cshuf（k）或csurr（k）；
2）将信号轮廓Y(i)分成不重叠的N s 段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；
为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N s 段数据；
3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：
Figure 936348DEST_PATH_IMAGE042
Figure DEST_PATH_IMAGE043
；
y v (i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；
4） 计算第q阶波动函数的平均值：
Figure 943487DEST_PATH_IMAGE044
；
5）如果x(k)存在自相似特征，则第q阶波动函数的平均值F q (s)和时间尺度s之间存在幂律关系：
F q (s)~s H(q)；
当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：
Figure 295971DEST_PATH_IMAGE010
；
6）对步骤5）中的公式两边取对数可得ln[F q (s)]=H(q)ln(s)+c，c为常数，由此可以获得直线的斜率H(q)。
7.根据权利要求1所述的一种ITD多尺度波动分析状态监测方法，其特征在于：所述步骤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_IMAGE045
；
3）令J对的ak偏导数
Figure 363153DEST_PATH_IMAGE046
，k=1,2,…,m，这时(a1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure DEST_PATH_IMAGE047
，其中，RT代表R的转置矩阵，R-1代表R的逆矩阵；
4）rk(t)分别选取二阶多项式、三阶多项式、双曲线、指数曲线、对数曲线和复合函数曲线进行计算，然后对比各种曲线形式所产生的最小二乘指标J，选取J最小时所对应的曲线形式为rk(t)的形式。
8.根据权利要求1所述的一种ITD多尺度波动分析状态监测方法，其特征在于：所述步骤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 863404DEST_PATH_IMAGE048
。
9.实施如权利要求1-8之一所述的一种ITD多尺度波动分析状态监测方法的装置，其特征在于：所述装置包括以下部分：数据线，加速度传感器，数据采集卡，机箱、笔记本电脑和信号分析软件，加速度传感器通过数据线与数据采集卡连接，数据采集卡安装在机箱内，机箱通过数据线与笔记本电脑连接，信号分析软件安装在笔记本电脑上，信号分析软件用来实现所述算法。
CN202011240502.3A 2020-11-09 2020-11-09 一种itd多尺度波动分析状态监测方法及装置 Withdrawn CN112697479A (zh)

Priority Applications (1)

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

Applications Claiming Priority (1)

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

Publications (1)

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

Family

ID=75506429

Family Applications (1)

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

Country Status (1)

Country Link
CN (1) CN112697479A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117541020A (zh) * 2024-01-04 2024-02-09 山东合能科技有限责任公司 一种城市排水泵站调度管理方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117541020A (zh) * 2024-01-04 2024-02-09 山东合能科技有限责任公司 一种城市排水泵站调度管理方法及系统
CN117541020B (zh) * 2024-01-04 2024-03-22 山东合能科技有限责任公司 一种城市排水泵站调度管理方法及系统

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
CN112697484A (zh) 一种ssd多尺度波动分析状态监测方法及装置
CN112697483A (zh) 一种lmd多尺度波动分析状态监测方法及装置
CN112697266A (zh) 一种itd和gzc机器状态监测方法及装置
CN108801630B (zh) 单通道盲源分离的齿轮故障诊断方法
CN110991055B (zh) 一种旋转类机械产品剩余寿命预测系统
CN112683393A (zh) 一种lmd设备故障诊断方法及系统
CN113375939A (zh) 基于svd和vmd的机械件故障诊断方法
CN111912521B (zh) 一种非平稳信号的频率检测方法和存储介质
CN112697479A (zh) 一种itd多尺度波动分析状态监测方法及装置
Zhang et al. Improved local cepstrum and its applications for gearbox and rolling bearing fault detection
CN112697472A (zh) 一种wd多尺度波动分析状态监测方法及装置
CN112697482A (zh) 一种vmd多尺度波动分析状态监测方法及装置
CN112697473A (zh) 一种emd和gzc机器状态监测方法及装置
CN112697263A (zh) 一种eemd多尺度波动分析状态监测方法及装置
CN112697477A (zh) 一种emd设备故障诊断方法及系统
CN112697470A (zh) 一种ssd设备故障诊断方法及系统
CN112881046A (zh) 一种vmd和gzc机器状态监测方法及装置
CN112683395A (zh) 一种elmd和gzc机器状态监测方法及装置
CN112697474A (zh) 一种nmd多尺度波动分析状态监测方法及装置
CN113435281A (zh) 基于改进hht变换的波纹补偿器故障诊断方法
CN112683392A (zh) 一种elmd多尺度波动分析状态监测方法及装置
CN112697265A (zh) 一种自适应多重分形方法及设备工况监测装置
CN112629849A (zh) 一种fdm多尺度波动分析状态监测方法及装置
CN112697476A (zh) 一种wd和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: 20210423

WW01 Invention patent application withdrawn after publication