CN103234748B - 一种基于敏感imf的克林根贝尔格锥齿轮故障诊断方法 - Google Patents

Abstract

本发明公开了一种基于敏感IMF的克林根贝尔格锥齿轮故障诊断方法，该方法包括：1、利用加速度传感器对克林根贝尔格锥齿轮箱进行测量，采集加速度振动信号；2、将采集的信号导入Matlab中，得到原始信号，对其进行EMD分解，得到一系列IMF分量；3、根据敏感度评估算法计算各IMF分量的敏感度，选择出敏感IMF分量；4、计算敏感IMF分量的瞬时能量谱并绘制出瞬时能量谱图，根据谱图中幅值的分布，准确地提取故障特征。本发明方法是一种有效的故障特征提取方法，能够运用于克林根贝尔格锥齿轮的故障诊断中，并能快速准确地提取出故障信息；同时为克林根贝尔格锥齿轮的故障诊断和特征提取提供重要的理论依据。

Description

一种基于敏感IMF的克林根贝尔格锥齿轮故障诊断方法

技术领域

本发明属于故障诊断技术领域，涉及一种克林根贝尔格锥齿轮故障诊断方法，更具体涉及一种基于敏感IMF的克林根贝尔格锥齿轮故障诊断方法。

背景技术

克林根贝尔格锥齿轮（Klingelnberg Spiral Bevel Gear）作为螺旋锥齿轮的两大齿制之一，具有传动平稳、承载能力高、硬齿面刮削技术等特点，从而特别适用于大功率和大扭矩重载传动领域，是重型高档数控机床、汽车传动系统、航空航天装备等重要领域中的核心传动部件。该种齿轮通常在重载、冲击、变载荷等复杂工况下工作，常产生点蚀、胶合、崩裂甚至断齿等现象，导致齿轮发生故障。若未能及时发现和排除故障，可能导致整台机械设备发生故障或中断，也可能带来巨大的安全隐患和巨额经济损失。因此，研究该齿轮的故障特征提取方法，尽早发现并及时消除故障对于保障机械设备安全、避免事故和巨额经济损失及提高设备使用性能，具有十分重要的意义。

国内外许多学者对齿轮故障诊断技术进行了较广泛而深入的研究，但多数都是考虑直齿圆柱齿轮的故障诊断技术，而克林根贝尔格锥齿轮故障特征提取的研究和方法很少。针对克林根贝尔格锥齿轮的故障诊断，可利用的诊断信息很多，包括温度，振动，噪声等，但振动信号最能直接反应其状态，且实现方法也简单快捷，因此，将振动信号的特征作为该种齿轮的状态反应。但如何从克林根贝尔格锥齿轮的振动信号中有效的提取出故障信息，目前的机械设备故障诊断技术仍未能很好的满足实际需求，仍然具有很大的研究潜力，探索和研究新的诊断方法、技术，仍是该领域的重要研究内容之一。实现该种齿轮故障的诊断，不仅能达到克林根贝尔格锥齿轮状态监测和故障诊断的目的，还为该种齿轮的故障诊断提供了理论依据。

发明内容

本发明的目的是提供一种基于敏感IMF的克林根贝尔格锥齿轮故障诊断方法，通过对振动信号进行分析和处理，准确快速地提取出齿轮故障特征，为克林根贝尔格锥齿轮的故障诊断和特征提取提供重要的理论依据，为齿轮的故障诊断技术提供有效借鉴。

本发明是采用以下技术手段实现的：

1、利用加速度传感器对克林根贝尔格锥齿轮箱进行测量，采集齿轮加速度振动信号作为待分析信号；

2、将采集的待分析信号导入Matlab软件中，得到原始信号x(t)，对原始信号进行EMD（Empirical Mode Decomposition）分解，得到一系列IMF（Intrinsic Mode Function）分量；

各IMF分量应满足两个条件：(1)在整个时间序列信号中，极值点的个数和过零点的个数相等，或至多相差一个；(2)在任一时间点上，由信号的局部极大值点形成的包络线和由局部极小值点形成的包络线的平均值为零，即信号关于时间轴局部对称。

2.1.确定信号x(t)的所有局部极值点，用三次样条曲线将所有的极大值点连接起来形成x(t)的上包络线，同样连接所有的极小值点形成x(t)的下包络线，上下包络线应该包络所有的数据点。记上下包络线的平均值为m₁(t)，并计算x(t)和m₁(t)的差值，得到：

h₁(t)=x(t)-m₁(t)。

2.2.对于不同的x(t)，h₁(t)可能满足IMF的条件，也可能不满足。若不满足IMF条件，此时将h₁(t)作为原始信号，重复得到h₁(t)的分解步骤，得到：

h₁₁(t)=h₁(t)-m₁₁(t)

其中，m₁₁(t)为h₁(t)的上、下包络线均值。

2.3.如果h₁₁(t)不满足IMF的条件，则分解继续，重复上述步骤2.1中的分解k次，得：

h_1k(t)=h_1(k-1)(t)-m_1k(t)

2.4.判定h_1k(t)是否为一个IMF分量，必须要有分解终止准则，它可以定义为连续两个分解结果之间的标准差SD值：

$\mathrm{SD}=\frac{\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{|h_{1(k-1)}\left(t\right)-h_{1k}\left(t\right)|}^2}{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{h}_{1(k-1)}^2\left(t\right)}$

2.5.当h_1k(t)满足IMF的条件或SD值小于某一设定值，即认为h_1k(t)是一个IMF分量，记C₁(t)=h_1k(t)，得到第一个IMF分量。

2.6.令r₁(t)=x(t)-C₁(t)，将r₁(t)作为新的待分析信号重复2.1至2.5的步骤，以得到第二个IMF分量，记为C₂(t)，重复循环n次，至此，得到信号x(t)的n个IMF分量，

$\left\{\begin{array}{c}{r}_1\left(t\right)-C_2\left(t\right)=r_2\left(t\right)\\ r_2\left(t\right)-C_3\left(t\right)=r_3\left(t\right)\\ .\\ .\\ .\\ r_{n-1}\left(t\right)-C_n\left(t\right)=r_n\left(t\right)\end{array}\right.$

当r_n(t)成为一个单调函数或不可再从中提取满足IMF条件的分量时，分解结束。至此，信号x(t)被分解为：

$x\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_i\left(t\right)+r_n\left(t\right)$

3、根据敏感度评估算法计算各IMF分量的敏感度，剔除与故障无关或噪声干扰的IMF分量，选择出敏感IMF分量；

3.1.计算信号x(t)和各IMF分量C_i(t)之间的相似性系数，并记为α_i，其中

${\mathrm{\α}}_i=\frac{{\∫}_{-\∞}^{+\∞}x\left(t\right)C_i\left(t\right)\mathrm{dt}}{\sqrt{{\∫}_{-\∞}^{+\∞}{x}^2\left(t\right)\mathrm{dt}\·{\∫}_{-\∞}^{+\∞}{C}_i^2\left(t\right)\mathrm{dt}}},i=\mathrm{1,2,3},...,n$

3.2.计算正常信号x_nor(t)和各IMF分量C_i(t)之间的相似性系数，并记为β_i，其中

${\mathrm{\β}}_i=\frac{{\∫}_{-\∞}^{+\∞}{x}_{\mathrm{nor}}\left(t\right)C_i\left(t\right)\mathrm{dt}}{\sqrt{{\∫}_{-\∞}^{+\∞}{x}_{\mathrm{nor}}^2\left(t\right)\mathrm{dt}\·{\∫}_{-\∞}^{+\∞}{C}_i^2\left(t\right)\mathrm{dt}}},i=\mathrm{1,2,3},...,n$

3.3.定义敏感度因子ξ，并计算各IMF分量C_i(t)的敏感度；

${\mathrm{\ξ}}_i=\frac{|{\mathrm{\α}}_i-{\mathrm{\β}}_i|}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}|{\mathrm{\α}}_i-{\mathrm{\β}}_i|}\×100\%,i=\mathrm{1,2,3},...,n$

3.4.根据ξ值的大小来选择最能反应故障特征信息的IMF分量，定义为敏感IMF分量。通常情况下选择ξ值较大者为敏感IMF分量。

4、计算敏感IMF分量的瞬时能量谱，采用Matlab软件绘制出瞬时能量谱图，根据瞬时能量谱图中幅值的分布，准确地提取故障特征。

本发明的特点在于基于EMD理论，提出了一种齿轮故障诊断方法，通过此方法可以快速准确地提取出克林根贝尔格锥齿轮故障特征，达到该种齿轮故障诊断的目的。发明内容包括四部分。在第一部分中，主要是利用加速度传感器对克林根贝尔格锥齿轮箱进行测量，采集齿轮加速度振动信号作为待分析信号；在第二部分中，将采集的信号导入Matlab软件中，得到原始信号x(t)，对原始信号进行EMD分解，得到一系列IMF分量；在第三部分中，主要是根据敏感性评估算法计算各IMF分量的敏感度，剔除与故障无关或噪声干扰的IMF分量，选择出敏感IMF分量；在第四部分中，通过计算敏感IMF分量的瞬时能量谱，采用Matlab软件绘制出瞬时能量谱图，根据瞬时能量谱图中幅值的分布，准确地提取故障特征。最后通过实例论证了本发明提出的方法。

通过下面的描述并结合附图说明，本发明会更加清晰，附图说明用于解释本发明方法及实施例。

附图说明

图1基于敏感IMF的克林根贝尔格锥齿轮故障诊断方法流程图

图2本发明信号采集系统结构简图

图3本发明实施例正常信号时域波形

图4本发明实施例故障信号时域波形

图5本发明实施例故障信号经EMD分解后得到的IMF分量（前四个IMF）

图6本发明实施例敏感IMF分量时域波形

图7本发明实施例敏感IMF分量的IE谱

具体实施方式

本发明实施例的一种基于敏感IMF的克林根贝尔格锥齿轮故障诊断方法流程图如图1所示，下面结合流程图对本发明的步骤作详细说明。具体实施步骤如下：

第一步：利用加速度传感器对克林根贝尔格锥齿轮箱进行测量，采集齿轮加速度振动信号作为待分析信号；

步骤（1）：信号采集系统的布置及相关参数

图2为信号采集系统布置结构简图。如图2中所示，采集系统主要包括三部分，第一部分是为系统提供动力的电动机M，且电机的速度由速度控制器控制，以保证系统中齿轮的转速要求；第二部分是动力传动部分，为一对克林根贝尔格锥齿轮Z1、Z2，传动比为1:2.5，齿轮重合度在2和3之间；最后一部分是一个交流电动机G，提供系统中的负载扭矩。齿轮副由电机驱动，齿轮Z1安装在轴I（输入轴）上，齿轮Z2安装在轴II（输出轴）上，加速度传感器布置在输入轴轴承盖处，采集齿轮振动加速度信号。本次采集中控制轴I（输入轴）转速为1440r/min，由齿轮的各种参数及啮合频率计算公式可以计算出各轴的转速、转频及啮合频率如表1所示。

表1齿轮副有关参数

步骤（2）：振动加速度信号的采集

先采集正常齿轮副的振动加速度信号作为参考信号x_nor(t)，然后在齿轮Z2上模拟故障，且模拟故障类型为磨损。此时再通过安装在输入轴处的加速度传感器采集发生故障后的齿轮副振动加速度信号，作为故障特征提取分析的原始信号x(t)。实验采样频率为10240Hz，在相同工况下，采集得到两种情况下的振动加速度信号，部分振动信号幅值如表2所示。图3和图4分别为采集的正常信号和故障信号的时域波形。

表2采集的部分实验数据

第二步：将采集的信号导入Matlab软件中，得到原始信号x(t)，对原始信号进行EMD分解，得到一系列IMF分量；

1）、记故障信号为x(t)，找出该信号的所有极值点，用三次样条曲线将所有的极大值点连接起来形成x(t)的上包络线，同样连接所有的极小值点形成x(t)的下包络线，上下包络线应该包络所有的数据点。记上下包络线的平均值为m₁(t)，计算x(t)和m₁(t)的差值得到：

h₁(t)=x(t)-m₁(t)

2）、用h₁(t)代替原来的信号x(t)重复1）中的过程，得到：

h_1k(t)=h_1(k-1)(t)-m_1k(t)

定义标准差SD， $\mathrm{SD}=\frac{\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{|h_{1(k-1)}\left(t\right)-h_{1k}\left(t\right)|}^2}{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{h}_{1(k-1)}^2\left(t\right)};$ 经验表明SD一般取0.2-0.3。当h_1k(t)满足IMF的条件或SD值小于某一设定值，即认为h_1k(t)是一个IMF分量，记C₁(t)=h_1k(t)，得到第一个IMF分量。

3）、令r₁(t)=x(t)-C₁(t)，对r₁(t)重复以上两个步骤，得到r_i(t)=r_i-1(t)-C_i(t)(i=1,2,...,n)直到r_n(t)成为一个单调函数或不可再分解为止。至此，信号x(t)被分解为：

$x\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_i\left(t\right)+r_n\left(t\right)$

本发明实施例故障信号经过EMD分解得到10阶IMF分量，其中前四个IMF分量的时域波形如图5中所示。

第三步：根据敏感度评估算法计算各IMF分量的敏感度，剔除与故障无关或噪声干扰的IMF分量，选择出敏感IMF分量；

计算信号x(t)和C_i(t)之间的相似性系数α_i；并计算正常信号x_nor(t)和C_i(t)之间的相似性系数β_i。由敏感度评估算法计算各阶C_i(t)的敏感度，得到结果如表3中所示。根据表3中ξ的值最终选择IMF1分量和IMF2分量为敏感IMF分量。

表3故障信号各IMF分量的α、β和ξ值

第四步：克林根贝尔格锥齿轮故障特征提取

1）、对选择的敏感IMF分量作Hilbert变换，得到：

$Y_i\left(t\right)=H\left[C_i\left(t\right)\right]=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}\frac{C_i\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ},i=\mathrm{1,2}$

2）、构建解析信号Z_i(t)，如下：

$Z_i\left(t\right)=C_i\left(t\right)+j{Y}_i\left(t\right)=a_i\left(t\right)e^{j{\mathrm{\θ}}_i\left(t\right)},i=\mathrm{1,2}.$

式中，a_i(t)和θ_i(t)分别为幅值函数和相位函数，并通过下面等式定义；

$a_i\left(t\right)=\sqrt{C_i^2\left(t\right)+Y_i^2\left(t\right)},i=\mathrm{1,2}$

${\mathrm{\θ}}_i\left(t\right)=\arctan \left[\frac{Y_i\left(t\right)}{C_i\left(t\right)}\right],i=\mathrm{1,2}$

3）、计算敏感IMF的Hilbert谱，记为H(f,t)，其中f为瞬时频率，则：

$H(f,t)=\mathrm{Re}\left\{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{a}_i\left(t\right)e^{[j\∫f_i\left(t\right)\mathrm{dt}]}\right\}$

4）、将H(f,t)的幅值平方对频率进行积分，得到瞬时能量谱，其反应信号在某个局部各种频率成分的能量总和。定义瞬时能量谱IE(t)，即；

$\mathrm{IE}\left(t\right)={\∫}_0^f{H}^2(f,t)\mathrm{df}$

图7为敏感IMF分量的IE谱图。从图7中明显地看出存在周期性峰值点，且相邻两峰值之间的时间间隔大约为0.104s，接近系统中齿轮轴II的旋转周期（1/9.6=0.10417s）。由此，准确地提取出了齿轮故障特征，可以推断出轴II上的齿轮2发生了故障，准确地诊断出故障发生位置，实现了克林根贝尔格锥齿轮的故障诊断。

由于克林根贝尔格锥齿轮具有重合度高、传动平稳等特点，给该种齿轮的故障特征提取及故障诊断带来一定的难度。且故障信号经EMD分解后得到的一系列IMF分量，通常只有少数分量与故障特征紧密相关（称与故障紧密相关的IMF分量为敏感IMF分量），而其他则是与故障无关或者是噪声干扰成分，这同样给特征提取带来一定难度。而本发明方法能够准确选择出敏感IMF分量，减弱与故障无关或噪声干扰成份的干扰，加强故障特征信息的提取，结果直接反应故障特征。

通过以上实例分析总结出：本发明方法能够运用于克林根贝尔格锥齿轮的故障诊断中，并能快速准确地提取出故障特征信息，最终实现该种齿轮的故障诊断。本发明方法不仅为克林根贝尔格锥齿轮的故障诊断和特征提取提供重要的理论依据，且为齿轮的故障诊断技术提供了有效借鉴。

Claims (1)

1.一种基于敏感IMF的克林根贝尔格锥齿轮故障诊断方法，其特征在于，该方法包括以下步骤：

1)利用加速度传感器对克林根贝尔格锥齿轮箱进行测量，采集齿轮加速度振动信号作为待分析信号；其中分别采集正常克林根贝尔格锥齿轮齿轮副的振动加速度信号作为参考信号x_nor(t)，发生故障后的齿轮副振动加速度信号作为故障特征提取的待分析信号x(t)；

2)将采集的待分析信号x(t)导入Matlab软件中，得到原始信号，对原始信号进行EMD(Empirical ModeDecomposition)分解，得到一系列IMF(Intrinsic Mode Function)分量，其中，各IMF分量应当满足两个条件：

(a)在整个时间序列信号中，极值点的个数和过零点的个数相等，或至多相差一个，

(b)在任一时间点上，由信号的局部极大值点形成的包络线和由局部极小值点形成的包络线的平均值为零，即信号关于时间轴局部对称；

2.1)确定待分析信号x(t)的所有局部极值点，用三次样条曲线将所有的极大值点连接起来形成x(t)的上包络线，同样连接所有的极小值点形成待分析信号x(t)的下包络线，上下包络线应该包络所有的数据点，记上下包络线的平均值为m₁(t)，并计算待分析信号x(t)和m₁(t)的差值，得到：

h₁(t)＝x(t)-m₁(t)

2.2)对于不同的待分析信号x(t)，h₁(t)可能满足IMF的条件，若不满足IMF分量的条件，此时将h₁(t)作为原始信号，计算h₁(t)和m₁₁(t)的差值，得到：

h₁₁(t)＝h₁(t)-m₁₁(t)

其中，m₁₁(t)为h₁(t)的上、下包络线均值，

若此时h₁₁(t)不满足IMF分量的条件，则重复上述过程k次，得h_1k(t)＝h_1(k-1)(t)-m_1k(t)，直到h_1k(t)满足IMF分量的条件，或者连续两个分解结果之间的标准差SD值小于某一设定值，得到一个IMF分量，记为：

C₁(t)＝h_1k(t)

其中，标准差SD为：

2.3)令r₁(t)＝x(t)-C₁(t)，将r₁(t)作为新的待分析信号重复2.1)、2.2)的步骤，得到信号x(t)的n个IMF分量，直到r_n(t)成为一个单调函数或不可再分解为止，至此，信号x(t)被分解为：

3)根据敏感性评估算法计算各IMF分量的敏感度，剔除与故障无关或噪声干扰的IMF分量，选择出敏感IMF分量；

3.1)计算信号x(t)和各IMF分量C_i(t)之间的相似性系数，并记为α_i，其中，

3.2)计算正常信号x_nor(t)和各IMF分量C_i(t)之间的相似性系数，并记为β_i，其中，

3.3)定义敏感度因子ξ，并计算各IMF分量C_i(t)的敏感度；

3.4)根据ξ值的大小来选择最能反应故障特征信息的IMF分量，定义为敏感IMF分量，其中，选择ξ值较大者为敏感IMF分量；

4)计算敏感IMF分量的瞬时能量谱，采用Matlab软件绘制出瞬时能量谱图，根据瞬时能量谱图中幅值的分布，准确地提取故障特征。