CN113405823A - 一种基于迭代扩展本征模态分解的旋转机械故障诊断方法 - Google Patents

一种基于迭代扩展本征模态分解的旋转机械故障诊断方法 Download PDF

Info

Publication number
CN113405823A
CN113405823A CN202110536535.0A CN202110536535A CN113405823A CN 113405823 A CN113405823 A CN 113405823A CN 202110536535 A CN202110536535 A CN 202110536535A CN 113405823 A CN113405823 A CN 113405823A
Authority
CN
China
Prior art keywords
order
key
frequency
time
expression
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.)
Granted
Application number
CN202110536535.0A
Other languages
English (en)
Other versions
CN113405823B (zh
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.)
Tongji University
Original Assignee
Tongji University
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 Tongji University filed Critical Tongji University
Priority to CN202110536535.0A priority Critical patent/CN113405823B/zh
Publication of CN113405823A publication Critical patent/CN113405823A/zh
Application granted granted Critical
Publication of CN113405823B publication Critical patent/CN113405823B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • 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/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • 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
    • 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Discrete Mathematics (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明涉及一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，包括：获取旋转机械在变转速工况运行时的振动加速度信号；进行离散短时傅里叶变换，获得时频表达式；通过双向快速单线提取法，初步估计能量最显著的关键时频脊线；通过迭代扩展本征模态分解，根据振动加速度信号和关键时频脊线初始估计值，分离重构关键模态；对关键模态进行等角度重采样，获得稳态角域信号；通过离散傅里叶变换将角域转换至阶次域，分析各特征阶次，得到故障诊断结果。与现有技术相比，本发明在不借助控制电路、转速计等硬件设备下，通过准确提取和分析振动信号中关键本征成分，高效地对变转速尤其是大转速波动工况的机械设备进行信号处理和故障诊断。

Description

一种基于迭代扩展本征模态分解的旋转机械故障诊断方法
技术领域
本发明涉及旋转机械故障诊断领域，尤其是涉及一种基于迭代扩展本征模态分解的旋转机械故障诊断方法。
背景技术
变转速下机械设备的动态信号十分复杂，信号会出现调频、调幅、调相等非平稳特征，导致严重的“频谱混叠”而使故障诊断变得十分困难。然而，实际工程中变转速工况几乎无处不在，尤其是在启停车等大转速波动工况下更是难以分析和辨识。
阶次跟踪(Order tracking)是解决转速波动问题最直接和有效的方法，可以分为硬件阶次跟踪(HOT)，计算阶次跟踪(COT)，无键相阶次跟踪(TLOT)三类。其中，硬件阶次跟踪主要存在需要硬件设备支持、难以布置实施、成本高昂等问题；计算阶次跟踪主要存在需要安装键相触发器、计算精度依赖于转速信号质量等问题；无键相阶次跟踪主要存在算法计算消耗大、计算精度低、仅适用于微弱转速波动、人为干预多等问题。
国内有关阶次跟踪技术的发明专利申请有：
授权日为2016年1月20日，授权号为CN103499443B的中国专利中，公开了一种名称为“一种齿轮故障无键相角域平均计算阶次分析方法”的发明专利。其利用平滑伪Wigner-Ville时频分析和Viterbi最优路径搜索算法估计瞬时转频，然后基于估计的键相信号进行等角度重采样，有效地对变转速工况下齿轮箱进行无键相阶次跟踪，摆脱对硬件的依赖，但是存在仅适用于微弱转速波动、易受噪声干扰、无法获得阶次时域信息的缺陷。
授权日为2018年5月1日，授权号为CN105512369B的中国专利中，公开了一种名称为“基于阶次谱的Vold-Kalman滤波带宽优选方法”的发明专利。其利用计算阶次跟踪和Vold-Kalman阶次跟踪分别获得置信区间的阶次谱，通过计算剩余信号的标准差来选择带宽库中最佳的滤波带宽，有效地解决了Vold-Kalman阶次跟踪带宽选择的问题，提高了阶次分析的精度，但是存在仅适用于微弱转速波动、计算消耗较大、需要预知转速信息的缺陷。
公开日为2019年8月2日，公开号为CN110084208A的中国专利中，公开了一种名称为“一种自适应降噪并避免阶次混叠的计算阶次跟踪方法”的发明专利。其利用转速信息定义裕量频率，基于排列熵PE和差分进化算法优化的VMD对信号分解重构，然后进行计算阶次跟踪，有效地自适应降低噪声干扰，克服变工况下阶次跟踪的混叠问题，但是存在着计算消耗大、计算精度依赖于转速信号质量的缺陷；
公开日为2020年6月19日，公开号为CN111307460A的中国专利中，公开了一种名称为“基于计算阶次跟踪与谱峭度的滚动轴承故障诊断方法”的发明专利。其利用谱峭度和包络解调预处理振动信号，然后设定转速脉冲阈值进行计算阶次跟踪，有效地在变工况下提取轴承故障特征阶次，但是存在着需要键相触发器、需要设定脉冲阈值、人为干预较多的缺陷。
发明内容
本发明的目的就是为了克服上述现有技术存在需要预知转速、计算精度受限、算法计算消耗大、仅适用于微弱转速波动、人为干预多等缺陷而提供一种通过高效分析和提取振动信号中关键本征成分，实现变转速工况下机械设备故障诊断的基于迭代扩展本征模态分解的旋转机械故障诊断方法。
本发明的目的可以通过以下技术方案来实现：
一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，包括以下步骤：
步骤1：获取旋转机械在变转速工况运行时的振动加速度信号；
步骤2：对所述振动加速度信号进行离散短时傅里叶变换，获得时频表达式；
步骤3：通过双向快速单线提取法，获取所述时频表达式中能量最显著的关键时频脊线的初始估计值；
步骤4：通过迭代扩展本征模态分解，根据所述振动加速度信号和关键时频脊线初始估计值，分离重构关键模态及其瞬时频率和瞬时幅值；
步骤5：对所述关键模态进行等角度重采样，获得稳态角域信号；
步骤6：通过离散傅里叶变换将所述角域转换至阶次域，分析各特征阶次，得到旋转机械的故障诊断结果。
进一步地，所述步骤2包括以下步骤：
步骤201：设置窗长度Nwindow和傅里叶点数Lfourier，构造相应的窗函数Wf(t)；
步骤202：通过所述窗函数Wf(t)对振动加速度信号x(t)进行离散短时傅里叶变换，获得时频表达式TFR(t,f)。
进一步地，步骤3中，所述双向快速单线提取法具体为：定位所述时频表达式中的能量最大值点，并沿时间方向双向搜索Δf范围内频率维度的能量极大值点，构建所述关键时频脊线的初始估计值，所述Δf由预先设定。
进一步地，所述步骤3具体包括以下步骤：
步骤301：初始化最大频率差值Δf；
步骤302：定位所述时频表达式中能量最大值点(tem,fem)，该定位表达式为：
Figure BDA0003070065190000031
式中，TFR(t,f)为时频表达式，
Figure BDA0003070065190000032
为关键时频脊线初始估计值的起点；
沿时间方向双向搜索Δf范围内频率维度的能量极大值点：该搜索表达式为:
fR＝fem,fL＝fem
Figure BDA0003070065190000033
式中，
Figure BDA0003070065190000034
为关键时频脊线的初始估计值，N为采样点个数；
步骤303：获得关键时频脊线的初始估计值
Figure BDA0003070065190000035
进一步地，所述步骤4具体包括以下步骤：
步骤401：初始化参数：正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L，并获取振动加速度信号x(t)和关键时频脊线初始估计值
Figure BDA0003070065190000036
其中傅里叶展开阶数L的计算表达式为：
Figure BDA0003070065190000037
式中，BW为带宽，N为采样点数，fs为采样频率；
步骤402：对所述振动加速度信号x(t)进行希尔伯特变换，构造关键模态mkey(t)的解析式
Figure BDA0003070065190000038
Figure BDA0003070065190000039
Figure BDA00030700651900000310
其中，k＝0,...,end为迭代次数，
Figure BDA0003070065190000041
为关键瞬时频率的k次迭代值，η(t)为干扰噪声，t为离散采样时间，
Figure BDA0003070065190000042
为解析幅值：
Figure BDA0003070065190000043
其中，
Figure BDA0003070065190000044
为初始相位，
Figure BDA0003070065190000045
为关键瞬时频率的真实值，akey(t)为关键模态的瞬时幅值；
然后建立解析幅值
Figure BDA0003070065190000046
的冗余傅里叶展开式：
Figure BDA0003070065190000047
其中，
Figure BDA0003070065190000048
Figure BDA0003070065190000049
中的l为傅里叶级数的阶数；
步骤403：将所述
Figure BDA00030700651900000410
展开模型带入
Figure BDA00030700651900000411
中构造矩阵表达式：
Figure BDA00030700651900000412
Figure BDA00030700651900000413
Kk＝Ek·T为N×(2L+1)的核矩阵，
Figure BDA00030700651900000414
Figure BDA00030700651900000415
式中，d为矩阵索引位置；
然后求解最优化问题：
Figure BDA00030700651900000416
式中，
Figure BDA00030700651900000417
为信号解析式；
获得傅里叶参数矩阵
Figure BDA00030700651900000418
并重构得
Figure BDA00030700651900000419
Figure BDA00030700651900000420
Figure BDA00030700651900000421
式中，I为单位矩阵，H为矩阵转置；
步骤404：基于欧拉公式展开所述重构
Figure BDA00030700651900000422
取
Figure BDA00030700651900000423
实部和虚部，并通过反正切解调法估计得到误差
Figure BDA0003070065190000051
对
Figure BDA0003070065190000052
进行滤波光滑并更新瞬时频率
Figure BDA00030700651900000520
步骤405：将步骤404中更新后的所述
Figure BDA0003070065190000053
代入步骤403更新核矩阵Kk+1、参数矩阵
Figure BDA0003070065190000054
和解析幅值矩阵
Figure BDA0003070065190000055
重复计算预设的end次后停止，然后求解重构关键模态mkey(t)及其幅值akey(t)。
进一步地，所述
Figure BDA0003070065190000056
的欧拉公式为：
Figure BDA0003070065190000057
所述误差
Figure BDA0003070065190000058
的估计表达式为：
Figure BDA0003070065190000059
Figure BDA00030700651900000510
式中，Rk为复包络实部，Ik为复包络虚部。
进一步地，所述瞬时频率
Figure BDA00030700651900000511
的更新表达式为：
Figure BDA00030700651900000512
Figure BDA00030700651900000513
其中，
Figure BDA00030700651900000514
为改进二阶差分算子。
进一步地，所述关键模态mkey(t)的计算表达式为：
mkey＝u·A+v·B
Figure BDA00030700651900000515
Figure BDA00030700651900000516
所述幅值akey(t)的计算表达式为：
Figure BDA00030700651900000517
进一步地，所述步骤5具体包括以下步骤：
步骤501：对所述关键模态mkey(t)进行希尔伯特变换，反正切其虚部与实部的比值并解卷绕处理，获得瞬时相位曲线
Figure BDA00030700651900000518
步骤502：通过最小二乘拟合法，建立所述瞬时相位曲线
Figure BDA00030700651900000519
的时间-弧度映射关系ψkey(t)；
步骤503：确定采样阶次Om和等弧度间隔Δθ，反求所述映射关系
Figure BDA0003070065190000061
得到等弧度时刻序列tΔθ(k)和弧度序列radΔθ(k)；
步骤504：对所述时刻序列tΔθ(k)和振动加速度信号x(t)进行三次样条插值，获得等弧度幅值序列aΔθ(k)。
进一步地，所述步骤6具体包括以下步骤：
步骤601：对所述稳态角域信号加汉宁窗whanning(k)，根据所述采样阶次Om计算阶次序列O(k)，通过离散傅里叶变换计算幅值序列AΔθ(k)，转换至阶次域；
步骤602：在所述阶次域中，若一阶为能量显著且最低阶次，则直接分析与其成比例的各特征阶次，包括：一阶及其倍阶、故障特征阶次及其倍阶；
在所述阶次域中，若一阶非能量显著且最低阶次，则定位小于一阶的能量显著且最低的阶次，并求其倒数得放大倍数m，将所述阶次序列O(k)等比例放大m倍，然后重复执行步骤602。
与现有技术相比，本发明具有以下优点：
本发明新提出了一种双向快速单线提取法(BFSE)和迭代扩展本征模态分解法(IEIMD)，并与无键相阶次跟踪技术(TLOT)有效结合，实现旋转机械故障诊断，通过准确提取和分析振动信号中的关键本征成分，对变转速尤其是大转速波动工况的机械设备进行信号处理和故障诊断，具有下列区别于常用技术的显著优势：
(1)本发明无需控制电路、键相触发器等硬件设备，显著降低运维成本和应用难度，同时不依赖于转速信息，计算精度不受转速信号质量的影响；
(2)本发明计算速度相比同类方法提升2倍以上，瞬时频率提取精度提高20％以上，模态重构精度提高15％以上，特征阶次误差小于1％并且计算精度提高20％以上，无需较多人为干预，有利于嵌入式开发和实际工程应用；
(3)本发明适用于大转速波动(瞬时转频变化率大于30％)的情况，可以用于启停车、快速升速和降速等复杂工况；
(4)本发明精确重构关键阶次的时域波形、瞬时频率和瞬时幅值，有效地抑制了原信号中的噪声干扰，信息丰富完整，有利于深入分析故障。
附图说明
图1为本发明基于迭代扩展本征模态分解的旋转机械故障诊断方法的流程示意图；
图2(a)为本发明实施例中仿真信号时域图；
图2(b)为本发明实施例中仿真信号频域图；
图3(a)为本发明实施例中仿真信号时频谱及局部放大图；
图3(b)为本发明实施例中仿真信号BFSE关键时频脊线估计图；
图3(c)为本发明实施例中仿真信号IEIMD关键时频脊线提取图；
图3(d)为本发明实施例中仿真信号IEIMD关键时频脊线计算误差曲线图；
图4(a)为本发明实施例中仿真信号IEIMD关键模态波形及局部放大图；
图4(b)为本发明实施例中仿真信号加窗稳态角域图；
图5(a)为本发明实施例中仿真信号更新前阶次谱图；
图5(b)为本发明实施例中仿真信号更新后阶次谱图及局部放大图；
图6为本发明实施例中试验台结构示意图；
图7(a)为本发明实施例中试验台信号时频谱及局部放大图；
图7(b)为本发明实施例中试验台信号IEIMD关键时频脊线提取图；
图8为本发明实施例中试验台信号IEIMD关键模态波形及局部放大图；
图9(a)为本发明实施例中试验台信号更新前阶次谱图；
图9(b)为本发明实施例中试验台信号更新后阶次谱图及局部放大图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
如图1所示，本发明提供基于迭代扩展本征模态分解的旋转机械故障诊断方法，假设变转速振动信号为x(t)，该方法具体包含以下步骤：
步骤1：在旋转机械附近安装振动加速度传感器，采集变转速工况运行时的振动加速度信号x(t)；
步骤2：对所采集的振动加速度信号x(t)进行离散短时傅里叶变换，获得时频表达式，具体包括以下步骤：
步骤201：设置窗长度Nwindow和傅里叶点数Lfourier，构造相应窗函数Wf(t)；
步骤202：通过所述窗函数Wf(t)对振动加速度信号x(t)进行离散短时傅里叶变换(STFT)，获得时频表达式TFR(t,f)。
步骤3：通过双向快速单线提取法，初步估计所述时频表达式中能量最显著的关键时频脊线，具体包括以下步骤：
步骤301：初始化参数最大频率差值Δf；
步骤302：定位所述TFR(t,f)中能量最大值点(tem,fem)：
Figure BDA0003070065190000081
式中，TFR(t,f)为时频表达式，
Figure BDA0003070065190000082
为关键时频脊线初始估计值的起点；
并沿时间方向双向搜索Δf范围内频率维度的能量极大值点：
fR＝fem,fL＝fem
Figure BDA0003070065190000083
式中，
Figure BDA0003070065190000084
为关键时频脊线的初始估计值N为采样点个数。
步骤303：获得关键时频脊线的初始估计值
Figure BDA0003070065190000085
步骤4：通过迭代扩展本征模态分解，根据所述振动加速度信号和关键时频脊线初始值，分离重构关键模态及其瞬时频率和瞬时幅值，具体包括以下步骤：
步骤401：初始化正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L，输入振动加速度信号x(t)和关键时频脊线初始估计值
Figure BDA0003070065190000086
Figure BDA0003070065190000087
式中，BW为带宽，N为采样点数，fs为采样频率；
步骤402：对所述x(t)进行希尔伯特变换，构造关键模态mkey(t)的解析式
Figure BDA0003070065190000088
Figure BDA0003070065190000089
Figure BDA00030700651900000810
其中，k＝0,...,end为迭代次数，
Figure BDA0003070065190000091
为关键瞬时频率的k次迭代值，η(t)为干扰噪声，t为离散采样时间，
Figure BDA0003070065190000092
为解析幅值：
Figure BDA0003070065190000093
其中，
Figure BDA0003070065190000094
为初始相位，
Figure BDA0003070065190000095
为关键瞬时频率的真实值，akey(t)为关键模态的瞬时幅值；然后建立解析幅值
Figure BDA0003070065190000096
的冗余傅里叶展开式：
Figure BDA0003070065190000097
其中，
Figure BDA0003070065190000098
Figure BDA0003070065190000099
中的l为傅里叶级数的阶数；
步骤403：将所述
Figure BDA00030700651900000910
展开模型带入
Figure BDA00030700651900000911
中构造矩阵表达式：
Figure BDA00030700651900000912
Figure BDA00030700651900000913
Kk＝Ek·T为N×(2L+1)的核矩阵：
Figure BDA00030700651900000914
Figure BDA00030700651900000915
式中，d为矩阵索引位置；
然后求解最优化问题：
Figure BDA00030700651900000916
式中，
Figure BDA00030700651900000917
为信号解析式；
获得傅里叶参数矩阵
Figure BDA00030700651900000918
并重构得
Figure BDA00030700651900000919
Figure BDA00030700651900000920
Figure BDA00030700651900000921
其中，I为单位矩阵，H为矩阵转置；
步骤404：基于欧拉公式展开所述重构
Figure BDA00030700651900000922
为：
Figure BDA00030700651900000923
取
Figure BDA00030700651900000924
实部和虚部，并通过反正切解调法估计得到误差
Figure BDA00030700651900000925
Figure BDA0003070065190000101
Figure BDA0003070065190000102
式中，Rk为复包络实部，Ik为复包络虚部。
对
Figure BDA0003070065190000103
进行滤波光滑并更新瞬时频率
Figure BDA0003070065190000104
Figure BDA0003070065190000105
Figure BDA0003070065190000106
其中
Figure BDA0003070065190000107
为改进二阶差分算子；
步骤405：将所述
Figure BDA0003070065190000108
跳转至步骤403更新核矩阵Kk+1、参数矩阵
Figure BDA0003070065190000109
和解析幅值矩阵
Figure BDA00030700651900001010
重复计算end次后停止，然后求解重构关键模态mkey(t)及其幅值akey(t)：
mkey＝u·A+v·B
Figure BDA00030700651900001011
其中，
Figure BDA00030700651900001012
Figure BDA00030700651900001013
步骤5：对所述关键模态进行等角度重采样，获得稳态角域信号，具体包括以下步骤：
步骤501：对所述关键模态mkey(t)进行希尔伯特变换，反正切其虚部与实部的比值并解卷绕处理，获得瞬时相位曲线
Figure BDA00030700651900001014
步骤502：通过最小二乘拟合法，建立所述瞬时相位曲线
Figure BDA00030700651900001015
的时间-弧度映射关系ψkey(t)；
步骤503：确定采样阶次Om和等弧度间隔Δθ，反求所述映射关系
Figure BDA00030700651900001016
得到等弧度时刻序列tΔθ(k)和弧度序列radΔθ(k)；
步骤504：对所述时刻序列tΔθ(k)和振动加速度信号x(t)进行三次样条插值，获得等弧度幅值序列aΔθ(k)；
步骤6：通过离散傅里叶变换将所述角域转换至阶次域，分析各特征阶次，具体包括以下步骤：
步骤601：对所述等弧度幅值序列aΔθ(k)加汉宁窗whanning(k)，根据所述采样阶次Om计算阶次序列O(k)，通过离散傅里叶变换计算幅值序列AΔθ(k)，转换至阶次域；
步骤602：在所述阶次域中，若一阶为能量显著且最低阶次，则直接分析与其成比例的各特征阶次，包括：一阶及其倍阶、故障特征阶次及其倍阶；
步骤603：在所述阶次域中，若一阶非能量显著且最低阶次，则定位小于一阶的能量显著且最低的阶次，并求其倒数得放大倍数m，将所述阶次序列O(k)等比例放大m倍，然后按步骤602所述分析各特征阶次，得到故障诊断结果；
实施例1
如图2(a)所示，齿轮箱振动仿真信号x(t)的信噪比为-5dB，采样频率4000Hz，采样时间5s；如图2(b)所示，频谱X(t)中发生严重“频谱混叠”，无法分析故障信息。如图3(a)所示，2.5秒内转频从400Hz提升至1100Hz左右，瞬时频率变化率大于30％，具有明显大转速波动特性，由于交叉项限制和噪声干扰，时频谱图中边频带和两端时频脊线模糊；如图3(c)所示，IEIMD提取的关键时频脊线
Figure BDA0003070065190000111
与真实脊线匹配度极高；如图3(d)所示，
Figure BDA0003070065190000112
计算误差曲线整体平稳且均在0.1％以内，Vold-Kalman滤波误差曲线的两端偏差大。均方误差(MSE)和相对误差(RE)如表1所示，与同类方法相比，瞬时频率提取精度提高21％，尤其与Vold-Kalma滤波相比提升近3倍。
表1 MSE和RE指标
Figure BDA0003070065190000113
如图4(a)所示，IEIMD重构所得关键模态mkey(t)及其瞬时幅值akey(t)波形质量良好，物理特征明显。重构信噪比(SNR)如表2所示，从原-5dB提升至6.35dB，有效抑制了随机噪声的干扰，与同类方法相比，模态重构精度提高14％，尤其与Vold-Kalma滤波相比提高43％。
表2 SNR指标
Figure BDA0003070065190000114
Figure BDA0003070065190000121
如图5(a)所示，定位小于一阶的能量显著最低阶次，求其倒数为：
Figure BDA0003070065190000122
可知关键瞬时频率
Figure BDA0003070065190000123
即为第二啮合瞬时频率，将阶次序列O(k)放大35.91倍后得到标准的阶次谱图，如图5(b)所示，啮合特征阶次周围有明显的边频带，有效信息丰富。特征阶次计算误差如表3所示，本方法计算误差在0.5％以内，与同类方法相比，计算精度提高19％，尤其与Vold-Kalma滤波相比提高40％。
表3特征阶次计算误差
Figure BDA0003070065190000124
测试过程中软硬件环境，硬件环境：CPU为Intel(R)Core(TM)i5-6300HQ，主频为2.30GHz，内存为8.00GB，64位操作系统。软件环境：Windows 10企业版，MATLAB R2018a。计算消耗如表4所示，与同类方法相比，计算速度提升2倍左右。
表4计算消耗
Figure BDA0003070065190000125
综上所述，与同类方法相比，本方法的各项性能指标均有可观的提升，尤其与Vold-Kalman阶次跟踪相比，本方法更适用于大转速波动的工况。
实施例2
如图6所示，为MSF-PK5M机械故障模拟试验台，ICP加速度传感器型号为623C01，滚动轴承型号均为ER16K，具体参数如表5所示，本次试验选用内圈故障轴承。
表5轴承参数
Figure BDA0003070065190000126
如图7(a)所示，由于交叉项限制和噪声干扰，时频谱图分辨率低。如图8所示，IEIMD重构所得关键模态mkey(t)及其幅值akey(t)波形清晰，物理特性好，有利于深入分析故障。如图9(a)所示，定位小于一阶的能量显著最低阶次，求其倒数为：
Figure BDA0003070065190000131
可知关键瞬时频率
Figure BDA0003070065190000132
即为内圈故障瞬时频率，将阶次序列O(k)放大5.42倍后得到标准的阶次谱图，如图9(b)所示，故障特征阶次呈簇状，其2倍、3倍阶次均能量显著，故障信息丰富。
特征阶次计算误差如表6所示，本方法计算误差在0.2％以内，与同类方法相比，计算精度提高38％，尤其与COT相比提高近82％。这是因为COT计算精度受转速信号质量的影响，所以计算误差最大。
表6特征阶次计算误差
Figure BDA0003070065190000133
测试过程中软硬件环境，硬件环境：CPU为Intel(R)Core(TM)i5-6300HQ，主频为2.30GHz，内存为8.00GB，64位操作系统。软件环境：Windows 10企业版，MATLAB R2018a。本方法的计算消耗如表7所示，相比于同类方法，计算速度提升3倍左右。
表7计算消耗
Figure BDA0003070065190000134
综上所述，一种基于迭代扩展本征模态分解的无键相阶次跟踪技术的旋转机械故障诊断方法具有显著优越性，可以对变转速尤其是大转速波动工况的机械设备实现高效的信号处理与故障诊断，并有利于嵌入式开发和实际工程应用。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，包括以下步骤：
步骤1：获取旋转机械在变转速工况运行时的振动加速度信号；
步骤2：对所述振动加速度信号进行离散短时傅里叶变换，获得时频表达式；
步骤3：通过双向快速单线提取法，获取所述时频表达式中能量最显著的关键时频脊线的初始估计值；
步骤4：通过迭代扩展本征模态分解，根据所述振动加速度信号和关键时频脊线初始估计值，分离重构关键模态及其瞬时频率和瞬时幅值；
步骤5：对所述关键模态进行等角度重采样，获得稳态角域信号；
步骤6：通过离散傅里叶变换将所述角域转换至阶次域，分析各特征阶次，得到旋转机械的故障诊断结果。
2.根据权利要求1所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，所述步骤2包括以下步骤：
步骤201：设置窗长度Nwindow和傅里叶点数Lfourier，构造相应的窗函数Wf(t)；
步骤202：通过所述窗函数Wf(t)对振动加速度信号x(t)进行离散短时傅里叶变换，获得时频表达式TFR(t,f)。
3.根据权利要求1所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，步骤3中，所述双向快速单线提取法具体为：定位所述时频表达式中的能量最大值点，并沿时间方向双向搜索Δf范围内频率维度的能量极大值点，构建所述关键时频脊线的初始估计值，所述Δf由预先设定。
4.根据权利要求1所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，所述步骤3具体包括以下步骤：
步骤301：初始化最大频率差值Δf；
步骤302：定位所述时频表达式中能量最大值点(tem,fem)，该定位表达式为：
Figure FDA0003070065180000011
式中，TFR(t,f)为时频表达式，
Figure FDA0003070065180000012
为关键时频脊线初始估计值的起点；
沿时间方向双向搜索Δf范围内频率维度的能量极大值点，该搜索表达式为：
fR＝fem,fL＝fem
Figure FDA0003070065180000021
式中，
Figure FDA0003070065180000022
为关键时频脊线的初始估计值，N为采样点个数；
步骤303：获得关键时频脊线的初始估计值
Figure FDA0003070065180000023
5.根据权利要求1所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，所述步骤4具体包括以下步骤：
步骤401：初始化参数：正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L，并获取振动加速度信号x(t)和关键时频脊线初始估计值
Figure FDA0003070065180000024
其中傅里叶展开阶数L的计算表达式为：
Figure FDA0003070065180000025
式中，BW为带宽，N为采样点数，fs为采样频率；
步骤402：对所述振动加速度信号x(t)进行希尔伯特变换，构造关键模态mkey(t)的解析式
Figure FDA0003070065180000026
Figure FDA0003070065180000027
Figure FDA0003070065180000028
其中，k＝0,...,end为迭代次数，
Figure FDA0003070065180000029
为关键瞬时频率的k次迭代值，η(t)为干扰噪声，t为离散采样时间，
Figure FDA00030700651800000210
为解析幅值：
Figure FDA00030700651800000211
其中，
Figure FDA00030700651800000212
为初始相位，
Figure FDA00030700651800000213
为关键瞬时频率的真实值，akey(t)为关键模态的瞬时幅值；
然后建立解析幅值
Figure FDA00030700651800000214
的冗余傅里叶展开式：
Figure FDA00030700651800000215
其中，
Figure FDA00030700651800000216
Figure FDA00030700651800000217
中的l为傅里叶级数的阶数；
步骤403：将所述
Figure FDA00030700651800000218
展开模型带入
Figure FDA00030700651800000219
中构造矩阵表达式：
Figure FDA0003070065180000031
Figure FDA0003070065180000032
Kk＝Ek·T为N×(2L+1)的核矩阵，
Figure FDA0003070065180000033
Figure FDA0003070065180000034
式中，d为矩阵索引位置；
然后求解最优化问题：
Figure FDA0003070065180000035
式中，
Figure FDA0003070065180000036
为信号解析式；
获得傅里叶参数矩阵
Figure FDA0003070065180000037
并重构得
Figure FDA0003070065180000038
Figure FDA0003070065180000039
Figure FDA00030700651800000310
式中，I为单位矩阵，H为矩阵转置；
步骤404：基于欧拉公式展开所述重构
Figure FDA00030700651800000311
取
Figure FDA00030700651800000312
实部和虚部，并通过反正切解调法估计得到误差
Figure FDA00030700651800000313
对
Figure FDA00030700651800000314
进行滤波光滑并更新瞬时频率
Figure FDA00030700651800000315
步骤405：将步骤404中更新后的所述
Figure FDA00030700651800000316
代入步骤403更新核矩阵Kk+1、参数矩阵
Figure FDA00030700651800000317
和解析幅值矩阵
Figure FDA00030700651800000318
重复计算预设的end次后停止，然后求解重构关键模态mkey(t)及其幅值akey(t)。
6.根据权利要求5所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，所述
Figure FDA00030700651800000319
的欧拉公式为：
Figure FDA00030700651800000320
所述误差
Figure FDA00030700651800000321
的估计表达式为：
Figure FDA00030700651800000322
Figure FDA00030700651800000323
式中，Rk为复包络实部，Ik为复包络虚部。
7.根据权利要求6所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，所述瞬时频率
Figure FDA0003070065180000041
的更新表达式为：
Figure FDA0003070065180000042
Figure FDA0003070065180000043
其中，
Figure FDA0003070065180000044
为改进二阶差分算子。
8.根据权利要求5所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，所述关键模态mkey(t)的计算表达式为：
mkey＝u·A+v·B
Figure FDA0003070065180000045
Figure FDA0003070065180000046
所述幅值akey(t)的计算表达式为：
Figure FDA0003070065180000047
9.根据权利要求1所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，所述步骤5具体包括以下步骤：
步骤501：对所述关键模态mkey(t)进行希尔伯特变换，反正切其虚部与实部的比值并解卷绕处理，获得瞬时相位曲线
Figure FDA0003070065180000048
步骤502：通过最小二乘拟合法，建立所述瞬时相位曲线
Figure FDA0003070065180000049
的时间-弧度映射关系ψkey(t)；
步骤503：确定采样阶次Om和等弧度间隔Δθ，反求所述映射关系
Figure FDA00030700651800000410
得到等弧度时刻序列tΔθ(k)和弧度序列radΔθ(k)；
步骤504：对所述时刻序列tΔθ(k)和振动加速度信号x(t)进行三次样条插值，获得等弧度幅值序列aΔθ(k)。
10.根据权利要求1所述的一种基于迭代扩展本征模态分解的旋转机械故障诊断方法，其特征在于，所述步骤6具体包括以下步骤：
步骤601：对所述稳态角域信号加汉宁窗whanning(k)，根据所述采样阶次Om计算阶次序列O(k)，通过离散傅里叶变换计算幅值序列AΔθ(k)，转换至阶次域；
步骤602：在所述阶次域中，若一阶为能量显著且最低阶次，则直接分析与其成比例的各特征阶次，包括：一阶及其倍阶、故障特征阶次及其倍阶；
在所述阶次域中，若一阶非能量显著且最低阶次，则定位小于一阶的能量显著且最低的阶次，并求其倒数得放大倍数m，将所述阶次序列O(k)等比例放大m倍，然后重复执行步骤602。
CN202110536535.0A 2021-05-17 2021-05-17 一种基于迭代扩展本征模态分解的旋转机械故障诊断方法 Active CN113405823B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110536535.0A CN113405823B (zh) 2021-05-17 2021-05-17 一种基于迭代扩展本征模态分解的旋转机械故障诊断方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110536535.0A CN113405823B (zh) 2021-05-17 2021-05-17 一种基于迭代扩展本征模态分解的旋转机械故障诊断方法

Publications (2)

Publication Number Publication Date
CN113405823A true CN113405823A (zh) 2021-09-17
CN113405823B CN113405823B (zh) 2022-05-20

Family

ID=77678738

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110536535.0A Active CN113405823B (zh) 2021-05-17 2021-05-17 一种基于迭代扩展本征模态分解的旋转机械故障诊断方法

Country Status (1)

Country Link
CN (1) CN113405823B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113887360A (zh) * 2021-09-23 2022-01-04 同济大学 一种基于迭代扩展频散模态分解的频散波提取方法
CN114216676A (zh) * 2021-11-30 2022-03-22 上海海事大学 一种时变工况下无转速计的行星齿轮箱复合故障诊断方法
CN114942338A (zh) * 2022-05-09 2022-08-26 重庆大学 基于嵌入式重力加速度感知的转子或旋转件的转动参数预估方法及系统
CN113887360B (zh) * 2021-09-23 2024-05-31 同济大学 一种基于迭代扩展频散模态分解的频散波提取方法

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050102116A1 (en) * 2003-11-07 2005-05-12 National Chiao Tung University High-resolution intelligent rotor machine diagnostic system and method
EP2053375A1 (en) * 2007-10-24 2009-04-29 Abb Research Ltd. A method for detection and automatic identification of damage to rolling bearings
CN102542296A (zh) * 2012-01-10 2012-07-04 哈尔滨工业大学 一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法
CN103353396A (zh) * 2013-06-24 2013-10-16 西安交通大学 基于无时标短时相位解调的齿轮箱故障诊断方法
US20150100609A1 (en) * 2013-10-07 2015-04-09 Brigham Young University Compression of time-varying simulation data
WO2015196735A1 (zh) * 2014-06-23 2015-12-30 华南理工大学 基于啮合频率和频谱校正技术的风电齿轮箱阶次跟踪方法
CN107341504A (zh) * 2017-06-07 2017-11-10 同济大学 一种基于时序数据流行学习的机械设备故障诊断方法
CN109084826A (zh) * 2018-07-06 2018-12-25 同济大学 一种用于故障预测与健康管理的智能传感系统
CN109359633A (zh) * 2018-12-10 2019-02-19 西北工业大学 基于希尔伯特-黄变换和小波脊线的信号联合分类方法
CN109682600A (zh) * 2018-09-14 2019-04-26 温州大学 一种用于发动机主轴轴承故障诊断的改进变分模态分解诊断方法
CN111665051A (zh) * 2020-07-01 2020-09-15 天津大学 基于能量权重法的强噪声变转速条件下轴承故障诊断方法
CN112446329A (zh) * 2020-11-30 2021-03-05 广州大学 一种时变结构瞬时频率确定方法、系统、装置及存储介质

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050102116A1 (en) * 2003-11-07 2005-05-12 National Chiao Tung University High-resolution intelligent rotor machine diagnostic system and method
EP2053375A1 (en) * 2007-10-24 2009-04-29 Abb Research Ltd. A method for detection and automatic identification of damage to rolling bearings
CN102542296A (zh) * 2012-01-10 2012-07-04 哈尔滨工业大学 一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法
CN103353396A (zh) * 2013-06-24 2013-10-16 西安交通大学 基于无时标短时相位解调的齿轮箱故障诊断方法
US20150100609A1 (en) * 2013-10-07 2015-04-09 Brigham Young University Compression of time-varying simulation data
WO2015196735A1 (zh) * 2014-06-23 2015-12-30 华南理工大学 基于啮合频率和频谱校正技术的风电齿轮箱阶次跟踪方法
CN107341504A (zh) * 2017-06-07 2017-11-10 同济大学 一种基于时序数据流行学习的机械设备故障诊断方法
CN109084826A (zh) * 2018-07-06 2018-12-25 同济大学 一种用于故障预测与健康管理的智能传感系统
CN109682600A (zh) * 2018-09-14 2019-04-26 温州大学 一种用于发动机主轴轴承故障诊断的改进变分模态分解诊断方法
CN109359633A (zh) * 2018-12-10 2019-02-19 西北工业大学 基于希尔伯特-黄变换和小波脊线的信号联合分类方法
CN111665051A (zh) * 2020-07-01 2020-09-15 天津大学 基于能量权重法的强噪声变转速条件下轴承故障诊断方法
CN112446329A (zh) * 2020-11-30 2021-03-05 广州大学 一种时变结构瞬时频率确定方法、系统、装置及存储介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Y SAKUMA: ""Proper Orthogonal Decomposition of Fluctuating Pressure on Train Side and Train Lateral Vibrations"", 《年会一般演讲》 *
吕中亮: ""基于变分模态分解与优化多核支持向量机的旋转机械早期故障诊断方法研究"", 《中国博士学位论文全文数据库》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113887360A (zh) * 2021-09-23 2022-01-04 同济大学 一种基于迭代扩展频散模态分解的频散波提取方法
CN113887360B (zh) * 2021-09-23 2024-05-31 同济大学 一种基于迭代扩展频散模态分解的频散波提取方法
CN114216676A (zh) * 2021-11-30 2022-03-22 上海海事大学 一种时变工况下无转速计的行星齿轮箱复合故障诊断方法
CN114942338A (zh) * 2022-05-09 2022-08-26 重庆大学 基于嵌入式重力加速度感知的转子或旋转件的转动参数预估方法及系统
CN114942338B (zh) * 2022-05-09 2023-10-20 重庆大学 基于嵌入式重力加速度感知的转子或旋转件的转动参数预估方法及系统

Also Published As

Publication number Publication date
CN113405823B (zh) 2022-05-20

Similar Documents

Publication Publication Date Title
Sun et al. Bearing fault diagnosis based on EMD and improved Chebyshev distance in SDP image
Chen et al. High-accuracy fault feature extraction for rolling bearings under time-varying speed conditions using an iterative envelope-tracking filter
Feng et al. Joint amplitude and frequency demodulation analysis based on intrinsic time-scale decomposition for planetary gearbox fault diagnosis
CN113405823B (zh) 一种基于迭代扩展本征模态分解的旋转机械故障诊断方法
CN109520738B (zh) 基于阶次谱和包络谱的旋转机械滚动轴承故障诊断方法
Miao et al. Application of sparsity-oriented VMD for gearbox fault diagnosis based on built-in encoder information
Jiang et al. Sparse dictionary design based on edited cepstrum and its application in rolling bearing fault diagnosis
Wang et al. Maximum envelope-based autogram and symplectic geometry mode decomposition based gear fault diagnosis method
CN109765055B (zh) 基于ewt、谱有效值和knn的滚动轴承故障检测方法及系统
CN109827777A (zh) 基于偏最小二乘法极限学习机的滚动轴承故障预测方法
Rodopoulos et al. A parametric approach for the estimation of the instantaneous speed of rotating machinery
CN112101245A (zh) 基于频域窗函数的短时傅里叶变换机械冲击特征提取方法
CN110940524B (zh) 一种基于稀疏理论的轴承故障诊断方法
Faysal et al. Noise eliminated ensemble empirical mode decomposition for bearing fault diagnosis
Liu et al. Asymmetric penalty sparse model based cepstrum analysis for bearing fault detections
Chen et al. A visualized classification method via t-distributed stochastic neighbor embedding and various diagnostic parameters for planetary gearbox fault identification from raw mechanical data
CN107831013A (zh) 一种利用概率主分量分析增强循环双谱的轴承故障诊断方法
CN110987431B (zh) 一种基于tqwt辅助spc的轴承状态监测及故障诊断的方法
Lv et al. Longitudinal synchroextracting transform: A useful tool for characterizing signals with strong frequency modulation and application to machine fault diagnosis
Wei et al. Two-level variational chirp component decomposition for capturing intrinsic frequency modulation modes of planetary gearboxes
Abboud et al. Synchronous analysis of cyclo-non-stationary signals: A comprehensive study with aeronautic applications
CN112648220A (zh) 一种基于小波-近似熵的风机故障诊断方法
CN113281047A (zh) 一种基于变尺度Lempel-Ziv的轴承内外圈故障定量趋势诊断方法
Tang et al. Rolling bearing diagnosis based on an unbiased-autocorrelation morphological filter method
Xu et al. Rolling bearing fault feature extraction via improved SSD and a singular-value energy autocorrelation coefficient spectrum

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
GR01 Patent grant
GR01 Patent grant