CN114813107A - 一种降采样多周期微分均值的故障特征增强方法 - Google Patents

Abstract

本发明公开了一种降采样多周期微分均值的故障特征增强方法，属于故障诊断技术与信号处理分析技术领域；本发明方法针对在编码器信号中故障特征微弱的难题，提出本发明方法，该方法基于均值降采样抑制噪声特性、降低计算成本和减小存储空间的优势，结合多周期微分均值的累积特性，提出一种降采样多周期微分均值(Down‑sampling multi‑perioddifferential equalization，DSMPDM)方法对原始IAS信号中故障分量进行增强处理，进而抑制编码器安装误差、瞬时角速度估计误差和测量噪声等分量的干扰。本发明方法提出的DSMPDM方法可有效增强故障分量，进而消除编码器安装误差分量和测量噪声分量对故障特征辨识的干扰。

Description

一种降采样多周期微分均值的故障特征增强方法

技术领域

本发明涉及一种降采样多周期微分均值的故障特征增强方法，属于故障诊断技术与信号处理分析技术领域。

背景技术

齿轮和轴承是旋转机械的关键零部件，在运行中起传递扭矩和支撑作用，其健康程度直接影响旋转机械的运行精度、效率和寿命。因此，齿轮和轴承的状态检测成为故障诊断领域研究热点之一。

当齿轮、轴承等部件发生局部点蚀、局部剥落等故障时，齿轮间或滚动体与滚道间在故障位置处的接触刚度较无故障状态发生突变，即通过单位角度的时间发生规律性变化，对应的瞬时角速度(Instantaneous angular speed,IAS)产生规律性波动。此外，轴承作为旋转机械的支撑部件，其并不传递扭矩，但在径向载荷的作用下，滚动体通过故障位置的刚度变化也将引起IAS信号变化。

另一方面，由于编码器内径与轴装配误差，编码器安装误差在工程应用中无法避免；编码器安装误差的能量幅值与转速成正相关，即编码器安装误差对故障信号的调制作用随着转速上升而显著提升，即对故障分量的干扰程度是不同的。因此，编码器安装误差的消除对于有效揭示故障特征至关重要。

发明内容

编码器信号中早期故障特征较弱以及编码器安装误差的调制作用随着转速的上升而增大，进而导致故障特征可能无法有效辨识，为解决该问题，本发明提供了一种降采样多周期微分均值的故障特征增强方法，该方法基于均值降采样抑制噪声特性、降低计算成本和减小存储空间的优势，结合多周期微分均值的累积特性，提出降采样多周期微分均值(Down-sampling multi-period differential equalization，DSMPDM)方法对瞬时角速度原始IAS_i信号中故障特征进行增强。

本发明基于降采样多周期微分均值的故障特征增强方法如下：

步骤1：获得包含旋转机械故障信息的瞬时角速度信号；

通过PicoScope信号采集系统获取光学编码器中包含旋转机械故障信息的瞬时角位移和对应的时间，采用向前差分法计算获得瞬时角速度IAS_i信号，其计算式如下

式中，IAS_i表示第i时刻的瞬时角速度，i＝1,2,3,…，Δφ＝2π/N，Δt_i＝t_i+1-t_i；N表示编码器的光栅格数，t为时间。

步骤2：对瞬时角速度原始信号IAS_i进行平均降采样处理，获得平均降采样数据IAS_g ^D，进而抑制随机噪声的干扰和降低后续算法的计算成本，其计算式为：

式中，D为降采样倍数，g＝1,D+1,2D+1,...，IAS_g ^D表示在降采样倍数为D时计算的IAS信号。

步骤3：在平均降采样数据IAS_g ^D的基础上，基于故障时会引起IAS信号的规律性变化及其对应IAS分量较无故障状态具有突变趋势，采用具有对缓变信号值较小对突变信号值较大特点的微分算法，并考虑以理论故障角度周期作为加窗的角度间隔，结合多周期累计特性，采用降采样多周期微分均值(Down-sampling multi-period differentialequalization，DSMPDM)方法以增强故障特征，其计算式为：

式中N_s ^m表示相邻故障冲击间实际角度间隔与理论角度间隔的最大差值，N_w表示微分窗长，K为感兴趣故障冲击的周期数，微分算子

在冲击周期k时随机滑动角度范围[h,H]＝[q-ηN_s ^m,q+ηN_s ^m]，η为相邻故障冲击间滑动比例，在数据点为m和冲击周期数为k时的微分位置q＝m+(k-1)N/f_fault，k＝1,2,...,K，f_fault为理论故障特征阶次，Q＝m+(k-1)N/f_fault+N_w，m＝2N_s ^m,2N_s ^m+1,...,M，M＝length(IAS_g ^D)-KN/f_fault-N_w-2N_s ^m为处理后数据长度，N/f_fault表示故障理论角度间隔，length(IAS_g ^D)表示均值降采样信号IAS_g ^D的长度，N_s ^m的表达式为：

公式(3)需要设置6个参数：D、K、R、f_fault、N_w和η，其中，D是降低算法计算成本的主要参数之一，通常取值2～4；K是提升算法鲁棒性的主要参数，通常K≥3，K值越大，算法对感兴趣分量能量幅值的增强效果越好，但需要更高的计算成本；根据传动系统参数、齿轮齿数或者轴承尺寸参数，R和f_fault可以被计算获得；为有效包含故障波形分量，2floor(N_s ^m)<N_w<10floor(N_s ^m)，floor(·)为向下取整操作；相邻故障冲击间滑动比例η通常取值为(0,2]以满足轴承随机滑动、降采样误差和离散采样误差的干扰。

步骤4：对DSMPDM算法获得的滤波信号DS_MPDM(m)进行傅里叶变换，进而揭示故障特征，其计算式为：

F_DSMPDM＝FFT(DS_MPDM(m)) (5)

式中，FFT(·)为傅里叶变换。

本发明的有益效果是：

(1)本发明采用均值降采样操作既可以抑制随机噪声的干扰，又可以降低DSMPDM算法的计算成本，为提高DSMPDM算法的计算效率提供有效途径；

(2)本发明通过DSMPDM算法可有效增强故障特征分量的能量幅值，进而抑制编码器安装误差、测量噪声和背景噪声的干扰；

(3)本发明提出基于DSMPDM算法和阶次谱分析实现了故障特征增强提取。

附图说明

图1为本发明方法流程示意图；

图2为实施例1中仿真轴承故障数据和对应的包络阶次谱，其中图(a)为IAS_i波形，图(b)是图(a)的包络谱；

图3为实施例1中所提方法的分析结果，其中图(a)为DSMPDM算法的滤波信号，图(b)为DSMPDM算法获得的阶次谱；

图4为实施例1中不同理论获得的结果，其中图(a)为FK确定的优化解调频带，图(b)为基于图(a)优化解调频带的阶次谱，图(c)为MOMEDA确定的优化周期，图(d)为图(c)的阶次谱；

图5为实施例2中测试实验台；

图6为实施例2中获得的齿轮故障数据和对应的阶次谱，其中图(a)为IAS波形，图(b)是图(a)的包络谱；

图7为实施例2中所提方法的分析结果，其中图(a)为DSMPDM算法的滤波信号，图(b)为DSMPDM算法获得的阶次谱；

图8为实施例2中不同理论获得的结果，其中图(a)为FK确定的优化解调频带，图(b)为基于图(a)优化解调频带的阶次谱，图(c)为MOMEDA确定的优化周期，图(d)为图(c)的阶次谱；

图9为实施例3中获得的轴承故障数据和对应的阶次谱，其中图(a)为IAS波形，图(b)是图(a)的包络谱；

图10为实施例3中所提方法的分析结果，其中图(a)为DSMPDM算法的滤波信号，图(b)为DSMPDM算法获得的阶次谱；

图11为实施例3中不同理论获得的结果，其中图(a)为FK确定的优化解调频带，图(b)为基于图(a)优化解调频带的阶次谱，图(c)为MOMEDA确定的优化周期，图(d)为图(c)的阶次谱。

具体实施方式

下面结合本发明实施例中的附图，对本文发明实例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员没有做出创造性劳动前提下所获得的所有其他实施例，都属于本文发明保护的范围，本实施例中如无特殊说明的方法，均为常规方法。

在下面的描述中阐述了很多具体细节以便充分理解本发明，但是本发明还可以采用其他不同于在此描述的其他方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。

实施例1：本实例描述的是本发明方法用于轴承故障特征的提取，如图1所示，具体过程包括：

仿真分析采用的轴承内圈故障信号模型计算式为：

W_bpfi＝w(θ)+w_bpfi(θ)+w_o(θ)+n(θ)

A_m＝A₀cos(2πf_mθ+φ_A)+C_A

式中，W_bpfi表示具有噪声干扰的轴承内圈故障仿真信号，w_bpfi(θ)表示轴承内圈冲击衰减波形，A₀为冲击幅值，f_m表示调制频率，φ_A为初始角度，C_A为常量；w_o(θ)为编码器安装误差，θ_e和θ_t分别表示初始角度，n(θ)表示为IAS估计误差、高速计数器采样误差、编码盘刻蚀误差的综合，w(θ)表示平均角速度，ξ表示阻尼系数，f_n表示固频，ψ＝θ-jΘ-τ_j，角度序列θ＝2π/N,4π/N,6π/N...，n(θ)表示编码器测量噪声，ρ＝Δr/r是几何中心和旋转中心的偏心距离的比值，r为编码器孔直径，Δr是几何中心和旋转中心的偏心距离，β表示编码器旋轴与转轴之间的倾斜角度，初始角度θ_e∈[φ_e,2kπ+φ_e]，θ_t∈[φ_t,2kπ+φ_t]，C表示轴承故障冲击幅值。在本实施例仿真信号中，N＝10000，在转速为5rad/s的自然阶次频率f_n＝50×，Θ＝floor(N/f_reb)＝1941，ρ＝0.001，β＝0.03，A＝0.0001，ξ＝0.03，<τ_j＞＝0，max(|τ_j|)＝39，f_bpfi＝5.15×，R＝1，转频分量f_r1＝1×，信噪比为-25dB。

步骤1：仿真信号IAS_i波形和对应的包络阶次谱如图2(a)和图2(b)，仿真信号IAS_i的数据长度为10⁵，可见转频特征阶次谱线占主导地位，而故障特征谱线无法有效辨识。因此，需要对轴承故障分量进行增强处理。

步骤2：对原始信号IAS_i进行降采样处理，降采样倍数D＝2，g＝1,3,5,...；计算式为：

步骤3：对平均降采样信号IAS_g ^D中故障特征分量进行增强

3-1、设置参数：K＝7，f_reb＝5.15×，R＝1，N_w＝58，N_s ^m＝NR/49/f_reb＝39；

3-2、采用DSMPDM算法对故障分量进行增强处理，其计算式为：

式中，M＝length(IAS_i)-KN/f_reb-N_w-2N_s ^m＝10⁵-7×5000/5.15-58-78＝93067，m＝2N_s ^m,2N_s ^m+1,...,M＝78,79,...,96951，q＝m+(k-1)N/f_reb，k＝1,2,...,7，[h,H]＝[q-2N_s ^m,q+2N_s ^m]，Q＝m+(k-1)N/f_reb+N_w；DSMPDM算法获得的滤波信号如图3(a)所示；

步骤4：对DSMPDM算法获得的滤波信号DS_MPDM(m)进行傅里叶变换，进而揭示故障特征，其计算式为：

F_DSMPDM＝FFT(DS_MPDM(m))

获得的阶次谱如图3(b)所示。可见，轴承故障特征谱线可被有效辨识，编码器安装误差和测量噪声得到有效抑制。

为进一步显示本实施例方法的有效性，采用FK、MOMEDA算法分别对瞬时角速度IAS_i信号进行分析作为对比，其中FK中分解等级k＝7，MOMEDA中滤波长度为1000，周期点数为1941，搜索范围为[2,3000]。上述两种方法获得的结果如图4所示；可见转频阶次谱线的3倍频占主导地位，而轴承内圈故障特征阶次谱线及其二倍频谱线并不显著。

实施例2：本实施例描述的是本发明方法用于齿轮局部故障特征提取

本实施例采用如图5所示实验台进行验证，通过采样率为10⁷Hz的自制高速计数器采集系统获取光学编码器的瞬时角位移IAD信息和对应的时间信息，采用向前差分法计算瞬时角速度IAS_i，编码器型号为ETF100-H851007B，编码器每圈光栅线数N为5000，编码器安装在大齿轮转轴上。齿轮箱中小齿轮齿数为32，大齿轮齿数为48，为模拟齿轮局部故障，小齿轮表面加工局部点蚀故障，齿轮啮合频率为48×，小齿轮故障阶次为1.5×。

步骤1：获得的齿轮局部故障IAS_i波形如图6(a)所示，IAS_i数据长度为238853，进而直接对原始IAS_i信号进行包络阶次谱分析，如图6(b)所示。可以发现，啮合频率f_mesh＝48×位置处的转频调制谱线1×及其倍频谱线占主导地位。因此，需要对齿轮故障数据进行处理以抑制转频调制谱线1×及其倍频谱线的干扰。

步骤2：对原始信号IAS_i进行降采样处理，降采样倍数D＝4，g＝1,5,9,...；计算式为：

步骤3：对平均降采样信号IAS_g ^D中故障分量进行增强

3-1、设置参数：K＝5，f_reb＝1.5×，R＝1，N_w＝10，N_s ^m＝2；

3-2、采用DSMPDM算法对故障分量进行增强处理，其计算式为：

式中，M＝length(IAS_i)-KN/f_reb-N_w-2N_s ^m＝238853-5×5000/5.15-10-4＝233984，m＝2N_s ^m,2N_s ^m+1,...,M＝4,5,...,233984，q＝m+(k-1)N/f_reb，k＝1,2,...,5，[h,H]＝[q-2N_s ^m,q+2N_s ^m]，Q＝m+(k-1)N/f_reb+N_w。DSMPDM算法获得的滤波信号如图7(a)所示。

步骤4：对DSMPDM算法获得的滤波信号DS_MPDM(m)进行傅里叶变换，进而揭示故障特征，其计算式为：

F_DSMPDM＝FFT(DS_MPDM(m))

获得的阶次谱如图7(b)所示。可见小齿轮故障阶次谱线1.5×及其谐波谱线可有效辨识，编码器安装误差和测量噪声得到有效抑制。

为进一步显示本实施例方法的有效性，采用FK算法确定包含丰富齿轮故障信息的解调频带，如图8(a)所示，确定的优化中心频率为2187.5×，频带宽度为625×，优化解调频带对应的阶次谱如图8(b)所示。可以发现，转频调制阶次谱线1×占优，而与小齿轮故障相关的阶次谱线1.5×无法有效辨识。另外，采用MOMEDA算法对齿轮故障数据的分析结果如图8(c)和8(d)所示，其中滤波长度为1000，周期点数为3333，搜索范围为[5,3500]，大齿轮转频阶次谱线1×及其倍频谱线占优。

实施例3：本实施例描述的是本发明方法用于轴承故障特征提取

本实施例采用如图5所示的轴承实验台进行验证，通过采样率为10⁷Hz的自制高速计数器采集系统获取光学编码器的瞬时角位移IAD信息和对应的时间信息，并采用向前差分法估计IAS_i信号，编码器型号为ETF100-H851007B，编码器每圈光栅线数N为5000。使用的故障轴承型号为6202-2RZ，为模拟滚动轴承外圈故障，在轴承外圈上用线切割方法加工一宽度约为0.3mm，深度约为0.28mm的小槽，滚子直径d为9.52mm，节圆直径D为46mm，滚子数目n为9，接触角α为0。轴承故障特征阶次f_fault＝f_bpfo＝5.15×。

步骤1：获得的齿轮局部故障IAS_i波形如图9(a)所示，IAS_i数据长度为79617，进而直接对原始IAS_i信号进行包络阶次谱分析，如图9(b)所示。可以发现，啮合频率f_mesh＝48×位置处的转频调制谱线1×及其倍频谱线占主导地位。因此，需要对齿轮故障数据进行处理以抑制转频调制谱线1×及其倍频谱线的干扰。

步骤2：对原始信号IAS_i进行降采样处理，降采样倍数D＝4，g＝1,5,9,...；计算式为：

步骤3：对平均降采样信号IAS_g ^D中故障分量进行增强

3-1、设置参数：K＝13，f_reb＝5.15，R＝1，N_w＝58，N_s ^m＝NR/49/f_reb＝20；

3-2、采用DSMPDM算法对故障分量进行增强处理，其计算式为：

式中，M＝length(IAS_i)-KN/f_reb-N_w-2N_s ^m＝79617-13×5000/5.15-58-40＝66897，m＝2N_s ^m,2N_s ^m+1,...,M＝4,5,...,66897，q＝m+(k-1)N/f_reb，k＝1,2,...,5，[h,H]＝[q-2N_s ^m,q+2N_s ^m]，Q＝m+(k-1)N/f_reb+N_w。DSMPDM算法获得的滤波信号如图10(a)所示。

步骤4：对DSMPDM算法获得的滤波信号DS_MPDM(m)进行傅里叶变换，进而揭示故障特征，其计算式为：

F_DSMPDM＝FFT(DS_MPDM(m))

获得的阶次谱如图10(b)所示。可见，轴承外圈故障特征阶次谱线占优，编码器安装误差分量和齿轮啮合分量得到有效抑制。

为进一步显示本实施例方法的有效性，采用FK算法对原始轴承外圈故障数据进行分析作为对比，最大分解等级k为7，确定的优化解调频带(中心频率为1393×，频带宽度为26×)和对应的阶次谱如图11(a)和11(b)所示，可见轴承故障特征阶次谱线无法有效辨识；此外，采用MOMEDA算法获得结果如图11(c)和图11(d)，其中滤波长度为1000，周期点数为970，搜索范围为[5,3000]，可见在图11(d)中无法有效辨识轴承外圈故障特征阶次谱线。因此，FK算法和MOMEDA算法在本实验中无法有效揭示轴承故障特征。

本文中应用了具体的实例对本发明的原理及实施方式进行了参数，以上实例的说明只是用于帮助理解本发明及核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变指出，综上，本说明书内容不应理解为本发明的限制。

Claims (3)

1.一种降采样多周期微分均值的故障特征增强方法，其特征在于，步骤如下：

(1)采集光学编码器的瞬时角位移和对应的时间，采用向前差分法计算瞬时角速度IAS_i；

(2)对瞬时角速度原始信号IAS_i进行平均降采样处理，抑制随机噪声的干扰、降低存储空间和计算成本；基于均值降采样抑制噪声特性、降低计算成本和减小存储空间的优势，结合多周期微分均值的累积特性，采用降采样多周期微分均值方法对平均降采样处理后数据中故障分量进行增强处理，获得滤波信号；

(3)滤波信号通过阶次谱分析揭示故障特征。

2.根据权利要求1所述的降采样多周期微分均值的故障特征增强方法，其特征在于，对瞬时角速度原始信号IAS_i进行平均降采样处理是采用如下计算公式获得平均降采样数据IAS_g ^D：

式中，D为降采样倍数；g＝1,D+1,2D+1,...；IAS_g ^D表示在降采样倍数为D时计算的IAS信号。

3.根据权利要求1所述的降采样多周期微分均值的故障特征增强方法，其特征在于：采用降采样多周期微分均值方法对平均降采样处理后的数据中故障分量进行增强处理是采用如下公式进行：

式中，N_s ^m表示相邻故障冲击间实际角度间隔与理论角度间隔的最大差值，N_w表示微分窗长，K为感兴趣故障冲击的周期数，微分算子

在冲击周期k时随机滑动角度范围[h,H]＝[q-ηN_s ^m,q+ηN_s ^m]，η为相邻故障冲击间滑动比例，在数据点为m和冲击周期数为k时的微分位置q＝m+(k-1)N/f_fault，k＝1,2,...,K，f_fault为理论故障特征阶次，Q＝m+(k-1)N/f_fault+N_w，m＝2N_s ^m,2N_s ^m+1,...,M，M＝length(IAS_g ^D)-KN/f_fault-N_w-2N_s ^m为处理后数据长度，N/f_fault表示故障理论角度间隔，length(IAS_g ^D)表示均值降采样信号IAS_g ^D的长度；N_s ^m的表达式为：

Images