一种基于迭代扩展本征模态分解的旋转机械故障诊断方法
技术领域
本发明涉及旋转机械故障诊断领域，尤其是涉及一种基于迭代扩展本征模态分解的旋转机械故障诊断方法。
背景技术
变转速下机械设备的动态信号十分复杂，信号会出现调频、调幅、调相等非平稳特征，导致严重的“频谱混叠”而使故障诊断变得十分困难。然而，实际工程中变转速工况几乎无处不在，尤其是在启停车等大转速波动工况下更是难以分析和辨识。
阶次跟踪(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)，该定位表达式为：
式中，TFR(t,f)为时频表达式，为关键时频脊线初始估计值的起点；
沿时间方向双向搜索Δf范围内频率维度的能量极大值点：该搜索表达式为:
fR＝fem,fL＝fem
式中，为关键时频脊线的初始估计值，N为采样点个数；
步骤303：获得关键时频脊线的初始估计值
进一步地，所述步骤4具体包括以下步骤：
步骤401：初始化参数：正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L，并获取振动加速度信号x(t)和关键时频脊线初始估计值其中傅里叶展开阶数L的计算表达式为：
式中，BW为带宽，N为采样点数，fs为采样频率；
步骤402：对所述振动加速度信号x(t)进行希尔伯特变换，构造关键模态mkey(t)的解析式
其中，k＝0,...,end为迭代次数，为关键瞬时频率的k次迭代值，η(t)为干扰噪声，t为离散采样时间，为解析幅值：
其中，为初始相位，为关键瞬时频率的真实值，akey(t)为关键模态的瞬时幅值；
然后建立解析幅值的冗余傅里叶展开式：
其中， 中的l为傅里叶级数的阶数；
步骤403：将所述展开模型带入中构造矩阵表达式：
Kk＝Ek·T为N×(2L+1)的核矩阵，
式中，d为矩阵索引位置；
然后求解最优化问题：
式中，为信号解析式；
获得傅里叶参数矩阵并重构得
式中，I为单位矩阵，H为矩阵转置；
步骤404：基于欧拉公式展开所述重构取实部和虚部，并通过反正切解调法估计得到误差对进行滤波光滑并更新瞬时频率
步骤405：将步骤404中更新后的所述代入步骤403更新核矩阵Kk+1、参数矩阵和解析幅值矩阵重复计算预设的end次后停止，然后求解重构关键模态mkey(t)及其幅值akey(t)。
进一步地，所述的欧拉公式为：
所述误差的估计表达式为：
式中，Rk为复包络实部，Ik为复包络虚部。
进一步地，所述瞬时频率的更新表达式为：
其中，为改进二阶差分算子。
进一步地，所述关键模态mkey(t)的计算表达式为：
mkey＝u·A+v·B
所述幅值akey(t)的计算表达式为：
进一步地，所述步骤5具体包括以下步骤：
步骤501：对所述关键模态mkey(t)进行希尔伯特变换，反正切其虚部与实部的比值并解卷绕处理，获得瞬时相位曲线
步骤502：通过最小二乘拟合法，建立所述瞬时相位曲线的时间-弧度映射关系ψkey(t)；
步骤503：确定采样阶次Om和等弧度间隔Δθ，反求所述映射关系得到等弧度时刻序列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)：
式中，TFR(t,f)为时频表达式，为关键时频脊线初始估计值的起点；
并沿时间方向双向搜索Δf范围内频率维度的能量极大值点：
fR＝fem,fL＝fem
式中，为关键时频脊线的初始估计值N为采样点个数。
步骤303：获得关键时频脊线的初始估计值
步骤4：通过迭代扩展本征模态分解，根据所述振动加速度信号和关键时频脊线初始值，分离重构关键模态及其瞬时频率和瞬时幅值，具体包括以下步骤：
步骤401：初始化正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L，输入振动加速度信号x(t)和关键时频脊线初始估计值
式中，BW为带宽，N为采样点数，fs为采样频率；
步骤402：对所述x(t)进行希尔伯特变换，构造关键模态mkey(t)的解析式
其中，k＝0,...,end为迭代次数，为关键瞬时频率的k次迭代值，η(t)为干扰噪声，t为离散采样时间，为解析幅值：
其中，为初始相位，为关键瞬时频率的真实值，akey(t)为关键模态的瞬时幅值；然后建立解析幅值的冗余傅里叶展开式：
其中， 中的l为傅里叶级数的阶数；
步骤403：将所述展开模型带入中构造矩阵表达式：
Kk＝Ek·T为N×(2L+1)的核矩阵：
式中，d为矩阵索引位置；
然后求解最优化问题：
式中，为信号解析式；
获得傅里叶参数矩阵并重构得
其中，I为单位矩阵，H为矩阵转置；
步骤404：基于欧拉公式展开所述重构为：
取实部和虚部，并通过反正切解调法估计得到误差
式中，Rk为复包络实部，Ik为复包络虚部。
对进行滤波光滑并更新瞬时频率
其中为改进二阶差分算子；
步骤405：将所述跳转至步骤403更新核矩阵Kk+1、参数矩阵和解析幅值矩阵重复计算end次后停止，然后求解重构关键模态mkey(t)及其幅值akey(t)：
mkey＝u·A+v·B
其中，
步骤5：对所述关键模态进行等角度重采样，获得稳态角域信号，具体包括以下步骤：
步骤501：对所述关键模态mkey(t)进行希尔伯特变换，反正切其虚部与实部的比值并解卷绕处理，获得瞬时相位曲线
步骤502：通过最小二乘拟合法，建立所述瞬时相位曲线的时间-弧度映射关系ψkey(t)；
步骤503：确定采样阶次Om和等弧度间隔Δθ，反求所述映射关系得到等弧度时刻序列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提取的关键时频脊线与真实脊线匹配度极高；如图3(d)所示，计算误差曲线整体平稳且均在0.1％以内，Vold-Kalman滤波误差曲线的两端偏差大。均方误差(MSE)和相对误差(RE)如表1所示，与同类方法相比，瞬时频率提取精度提高21％，尤其与Vold-Kalma滤波相比提升近3倍。
表1 MSE和RE指标
如图4(a)所示，IEIMD重构所得关键模态mkey(t)及其瞬时幅值akey(t)波形质量良好，物理特征明显。重构信噪比(SNR)如表2所示，从原-5dB提升至6.35dB，有效抑制了随机噪声的干扰，与同类方法相比，模态重构精度提高14％，尤其与Vold-Kalma滤波相比提高43％。
表2 SNR指标
如图5(a)所示，定位小于一阶的能量显著最低阶次，求其倒数为：
可知关键瞬时频率即为第二啮合瞬时频率，将阶次序列O(k)放大35.91倍后得到标准的阶次谱图，如图5(b)所示，啮合特征阶次周围有明显的边频带，有效信息丰富。特征阶次计算误差如表3所示，本方法计算误差在0.5％以内，与同类方法相比，计算精度提高19％，尤其与Vold-Kalma滤波相比提高40％。
表3特征阶次计算误差
测试过程中软硬件环境，硬件环境：CPU为Intel(R)Core(TM)i5-6300HQ，主频为2.30GHz，内存为8.00GB，64位操作系统。软件环境：Windows 10企业版，MATLAB R2018a。计算消耗如表4所示，与同类方法相比，计算速度提升2倍左右。
表4计算消耗
综上所述，与同类方法相比，本方法的各项性能指标均有可观的提升，尤其与Vold-Kalman阶次跟踪相比，本方法更适用于大转速波动的工况。
实施例2
如图6所示，为MSF-PK5M机械故障模拟试验台，ICP加速度传感器型号为623C01，滚动轴承型号均为ER16K，具体参数如表5所示，本次试验选用内圈故障轴承。
表5轴承参数
如图7(a)所示，由于交叉项限制和噪声干扰，时频谱图分辨率低。如图8所示，IEIMD重构所得关键模态mkey(t)及其幅值akey(t)波形清晰，物理特性好，有利于深入分析故障。如图9(a)所示，定位小于一阶的能量显著最低阶次，求其倒数为：
可知关键瞬时频率即为内圈故障瞬时频率，将阶次序列O(k)放大5.42倍后得到标准的阶次谱图，如图9(b)所示，故障特征阶次呈簇状，其2倍、3倍阶次均能量显著，故障信息丰富。
特征阶次计算误差如表6所示，本方法计算误差在0.2％以内，与同类方法相比，计算精度提高38％，尤其与COT相比提高近82％。这是因为COT计算精度受转速信号质量的影响，所以计算误差最大。
表6特征阶次计算误差
测试过程中软硬件环境，硬件环境：CPU为Intel(R)Core(TM)i5-6300HQ，主频为2.30GHz，内存为8.00GB，64位操作系统。软件环境：Windows 10企业版，MATLAB R2018a。本方法的计算消耗如表7所示，相比于同类方法，计算速度提升3倍左右。
表7计算消耗
综上所述，一种基于迭代扩展本征模态分解的无键相阶次跟踪技术的旋转机械故障诊断方法具有显著优越性，可以对变转速尤其是大转速波动工况的机械设备实现高效的信号处理与故障诊断，并有利于嵌入式开发和实际工程应用。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。
