一种基于改进HHT的滚动轴承故障诊断方法

摘要

本发明公开了一种基于改进HHT的滚动轴承故障诊断方法，通过软筛分停止准则自适应地确定EMD和NHT的筛分迭代次数并进行筛分，得到瞬时幅值AMi[n]和瞬时频率fi[n]，再利用AMi[n]和fi[n]构建时频谱，并通过时频谱生成快速谱峭度图，选择峭度值最高的频带进行平方包络解调，找出故障特征频率点观测其幅值变化，从而确定故障。

说明书

技术领域

本发明属于滚动轴承故障诊断技术领域，更为具体地讲，涉及一种基于改进型HHT的滚动轴承故障诊断方法。

背景技术

滚动轴承是旋转机械中应用最广泛的一种部件，而且作为旋转机械关键部件之一，其工作状态的好坏会影响整个设备的运行状态。滚动轴承一旦发生故障，必将导致旋转机械结构失效，从而带来经济损失，严重时还会引发安全事故。因此，对滚动轴承进行故障诊断具有重要的工程意义。

当滚动轴承发生故障时，其振动信号为多分量的调幅-调频信号，在进行解调之前，需要将其分解为若干个单分量的调幅-调频信号。利用HHT变换对故障信号的分析，一般先采用经验模态分解(EMD)对故障振动信号进行模态分解，然后对分解所得的固有模态函数(IMF)进行希尔伯特变换(HT)得到时频谱图。但是，HT解调存在负频率问题，因此，采用归一化希尔伯特变换(NHT) 解调IMF信号。然而，EMD和NHT的筛分停止准则几乎都采用硬筛分停止准则，即需要有先验知识的经验丰富的专家预先设定阈值。这种硬筛分方法对于实际的振动信号不具有自适应性，不利于抑制EMD的模态混叠问题，也不能确保NHT是否解调完全。因此，本发明提供一种自适应地确定筛分次数的软筛分停止准则，首先用均方根(RMS))和超峭度(EK)两个特征确定目标函数，然后提出一种启发式机制，自适应地确定最佳的筛分迭代次数。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于改进型HHT的滚动轴承故障诊断方法，先通过软筛分停止准则进行筛分，然后提取谱峭度，选择峭度值最高的频带进行平方包络解调，找出故障特征频率点观测其幅值变化，从而确定故障。

为实现上述发明目的，本发明为一种基于改进型HHT的滚动轴承故障诊断方法，其特征在于，包括以下步骤：

(1)、采集滚动轴承振动信号

使用振动数据采集仪以采样频率F_s采集待检测滚动轴承在运行状态下垂直方向的振动信号，记为x[n]，n＝1,2,…,N_s，N_s为总采样点数；

(2)对采集的振动信号x[n]进行改进型EMD

(2.1)、令每个IMF的初始信号为r_i[n]，每次筛分过程的初始信号为h_ik[n]，>max；初始化i＝1，r_i[n]＝x[n]；

(2.2)、令k＝0，h_ik[n]＝r_i[n]；

(2.3)、令k＝k+1，确定h_ik-1[n]中所有极值点，然后分别采用三次样条曲线连接所有极大值点和极小值点，从而依次形成上包络线Emax_ik[n]和下包络线>ik[n]；

(2.4)、计算h_i(k-1)[n]的包络均值信号m_ik[n]；

(2.5)、计算估计的零均值信号h_ik[n]；

h_ik[n]＝h_i(k-1)[n]-m_ik[n]

(2.6)、计算h_ik[n]的包络均值信号m_i(k+1)[n]，再根据m_i(k+1)[n]计算f_ik；

f_ik＝RMS(m_i(k+1)[n])+|EK(m_i(k+1)[n])|

(2.7)、判断当前筛分次数k是否满足k≥3，如果满足，则执行步骤(2.8)，否则返回步骤(2.3)；

(2.8)、判断当前筛分次数k是否达到最大筛分次数I_max，如果达到，则IMF_i筛分过程终止，且IMF_i[n]＝h_ik[n]；否则执行步骤(2.9)；

(2.9)、如果h_ik[n]的极值点的个数和零交叉点的个数相等或相差1个，则执行步骤(2.10)，否则返回步骤(2.3)；

(2.10)、如果f_i(k-2)＜f_i(k-1)且f_i(k-1)＜f_ik，则IMF_i筛分过程终止，此时>i[n]＝h_i(k-2)[n]；否则返回步骤(2.3)；

(2.11)、令i＝i+1，计算IMF_i的初始信号r_i[n]，计算公式如下式：

r_i[n]＝r_i-1[n]-IMF_i-1[n]

(2.12)、判断r_i[n]是否为单调函数，如果是单调函数，则改进型EMD结束；否则重复步骤(2.2)—(2.11)，直到r_i[n]变成一个单调函数，则改进型EMD>

(3)、对每个IMF进行改进型NHT

(3.1)、令每个IMF的初始信号为IMF_i[n]，每次筛分过程的初始信号为>ik[n]，i表示第i个IMF，k表示第k次筛分；设置每个IMF的最大筛分次数为I_max；初始化i＝1；

(3.2)、令k＝0，s_ik[n]＝IMF_i[n]；

(3.3)、令k＝k+1，确定s_i(k-1)[n]中所有极大值点，然后采用滑动平均法连接所有的极大值点，形成包络信号a_ik[n]；

(3.4)、计算s_ik[n]；

s_ik[n]＝s_i(k-1)[n]/a_ik[n]

(3.5)、确定s_ik[n]的包络信号a_i(k+1)[n]，再通过下式计算零均值信号q_i(k+1)[n]；

q_i(k+1)[n]＝a_i(k+1)[n]-1

(3.6)、根据q_i(k+1)[n]计算f_ik；

f_ik＝RMS(q_i(k+1)[n])+|EK(q_i(k+1)[n])|

(3.7)、判断当前筛分次数k是否满足k≥3，如果满足，则执行步骤(3.8)，否则返回步骤(3.3)；

(3.8)、判断当前筛分次数k是否达到最大筛分次数I_max，如果达到，则IMF_i解调过程终止，FM_i[n]＝s_ik[n]，AM_i[n]＝IMF_i[n]/s_ik[n]；否则执行步骤(3.9)；

(3.9)、如果f_i(k-2)＜f_i(k-1)且f_i(k-1)＜f_ik，则IMF_i解调过程终止，此时>i[n]＝s_i(k-2)[n]，AM_i[n]＝IMF_i[n]/s_i(k-2)[n]；否则返回步骤(3.3)；

(3.10)、计算FM_i[n]信号的希尔伯特变换对则解析信号表示为：

(3.11)、计算FM_i[n]的包络信号AM_HT[n]；

(3.12)、更新FM_i[n]和AM_i[n]；

FM_i[n]＝FM_i[n]/AM_HT_i[n]

AM_i[n]＝AM_i[n]×AM_HT_i[n]

(3.13)、计算FM_i[n]的瞬时相位φ_i[n]；

(3.14)、计算FM_i[n]的瞬时频率f_i[n]；

(3.15)、令i＝i+1，再返回步骤(3.2)，直到所有的IMF解调结束；

(4)、计算频率段特征

(4.1)、通过步骤(3)计算出的所有瞬时幅值AM_i[n]和瞬时频率f_i[n]得到改进型HHT的时频谱；

(4.2)、设置分割深度j，j＝1,2,…,h；h为分割层数，将时频谱在频域上按照分割深度j进行分割，得到一系列谱峭度其中，τ表示在同一分割深度上的频带分割位置；

(5)、获得快速谱峭度图

(5.1)、将步骤(4)得到的一系列谱峭度构成快速谱峭度矩阵，再通过 MATLAB中的imagesc函数将快速谱峭度矩阵生成快速谱峭度图；

(5.2)、记录下快速谱峭度图中峭度最大的频带所在的中心频率和带宽；

(6)、处理目标频带

(6.1)、将步骤(5)得到的中心频率和带宽作为滤波器的参数，然后用滤波器从原始信号中滤出目标信号；

(6.2)、先对目标信号的包络进行平方，再做傅里叶变换，得到其平方包络谱；

(6.3)、在平方包络谱上先找出轴承的故障特征频率点，然后观测故障特征频率处的幅值变化，如果故障特征频率处的幅值基本保持不变，则判断滚动轴承正常；如果故障特征频率处的幅值变化显著，则判断滚动轴承故障。

本发明的发明目的是这样实现的：

本发明为一种基于改进型HHT的滚动轴承故障诊断方法，通过软筛分停止准则自适应地确定EMD和NHT的筛分迭代次数并进行筛分，得到瞬时幅值 AM_i[n]和瞬时频率f_i[n]，再利用AM_i[n]和f_i[n]构建时频谱，并通过时频谱生成快速谱峭度图，选择峭度值最高的频带进行平方包络解调，找出故障特征频率点观测其幅值变化，从而确定故障。

同时，本发明为一种基于改进型HHT的滚动轴承故障诊断方法还具有以下

有益效果：

(1)、本发明提供一种软筛分停止准则自适应地确定筛分过程的最佳迭代次数；

(2)、提供的软筛分停止准则用于EMD和NHT的筛分过程中，从而改善 EMD的模态混叠问题，以及提高NHT解调的精度；

(3)、由改进型EMD和改进型NHT组合得到的改进型HHT方法作为快速谱峭度图的时频表示，提高了快速谱峭度图的时频估计精度，而且被成功用于高速列车的轴承故障诊断任务中。

附图说明

图1是本发明基于改进型HHT的滚动轴承故障诊断方法流程图；

图2是改进型EMD流程图；

图3是改进型NHT流程图；

图4是快速谱峭度矩阵；

图5是高速列车轮对轴承测试装置；

图6是外圈剥落故障图；

图7是外圈剥落故障下的时域波形；

图8是Rilling’s HHT,Wu’s HHT,Wang’s HHT和本发明的改进型HHT滤波后的平方包络幅值谱。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

为了方便描述，先对具体实施方式中出现的相关专业术语进行说明：

EMD(Empirical Mode Decomposition)：经验模态分解；

NHT(Normalized Hilbert Transform)：归一化希尔伯特变换；

HHT(Hilbert Huang Transform):希尔伯特黄变换；

IMF(intrinsic mode functions):固有模态函数；

AM(Amplitude Modulation):幅值调制；

FM(Frequency Modulation):频率调制；

RMS(Root Mean Square):均方根；

EK(Excess Kurtosis):超峭度。

图1是本发明基于改进型HHT的滚动轴承故障诊断方法流程图。

在本实施例中，如图1所示，本发明一种基于改进型HHT的滚动轴承故障诊断方法，包括以下步骤：

S1、采集滚动轴承振动信号

使用振动数据采集仪以采样频率F_s采集待检测滚动轴承在运行状态下垂直方向的振动信号，记为x[n]，n＝1,2,…,N_s，N_s为总采样点数；

S2、对采集的振动信号x[n]进行改进型EMD

S2.1、如图2所示，令每个IMF的初始信号为r_i[n]，每次筛分过程的初始信号为h_ik[n]，i表示第i个IMF，k表示第k次筛分；设置每个IMF的最大筛分次数为I_max＝30；初始化i＝1，r_i[n]＝x[n]；

S2.2、令k＝0，h_ik[n]＝r_i[n]；

S2.3、令k＝k+1，确定h_ik-1[n]中所有极值点，包括极大值点和极小值点，然后分别采用三次样条曲线连接所有极大值点和极小值点，从而分别依次形成上包络线Emax_ik[n]和下包络线Emin_ik[n]；

S2.4、计算h_i(k-1)[n]的包络均值信号m_ik[n]；

S2.5、计算估计的零均值信号h_ik[n]；

h_ik[n]＝h_i(k-1)[n]-m_ik[n]

S2.6、计算h_ik[n]的包络均值信号m_i(k+1)[n]，再根据m_i(k+1)[n]计算f_ik；

f_ik＝RMS(m_i(k+1)[n])+|EK(m_i(k+1)[n])|

其中，RMS(m_i(k+1)[n])和EK(m_i(k+1)[n])计算方法为：

其中，是m_i(k+1)[n]的算术平均值；

S2.7、判断当前筛分次数k是否满足k≥3，如果满足，则执行步骤S2.8，否则返回步骤S2.3；

S2.8、判断当前筛分次数k是否达到最大筛分次数I_max，如果达到，则IMF_i筛分过程终止，且IMF_i[n]＝h_ik[n]；否则执行步骤S2.9；

S2.9、如果h_ik[n]的极值点的个数和零交叉点的个数相等或相差1个，则执行步骤S2.10，否则返回步骤S2.3；

S2.10、如果f_i(k-2)＜f_i(k-1)且f_i(k-1)＜f_ik，则IMF_i筛分过程终止，此时>i[n]＝h_i(k-2)[n]；否则返回步骤S2.3；

S2.11、令i＝i+1，计算IMF_i的初始信号r_i[n]，计算公式如下式：

r_i[n]＝r_i-1[n]-IMF_i-1[n]

S2.12、判断r_i[n]是否为单调函数，如果是单调函数，则改进型EMD结束；否则重复步骤S2.2—S2.11，直到r_i[n]变成一个单调函数，则改进型EMD结束；

S3、对每个IMF进行改进型NHT

S3.1、如图3所示，令每个IMF的初始信号为IMF_i[n]，每次筛分过程的初始信号为s_ik[n]，i表示第i个IMF，k表示第k次筛分；设置每个IMF的最大筛分次数为I_max＝30；初始化i＝1；

S3.2、令k＝0，s_ik[n]＝IMF_i[n]；

S3.3、令k＝k+1，确定s_i(k-1)[n]中所有极大值点，然后采用滑动平均法连接所有的极大值点，形成包络信号a_ik[n]；

S3.4、计算s_ik[n]；

s_ik[n]＝s_i(k-1)[n]/a_ik[n]

S3.5、确定s_ik[n]的包络信号a_i(k+1)[n]，再通过下式计算零均值信号q_i(k+1)[n]；

q_i(k+1)[n]＝a_i(k+1)[n]-1

S3.6、根据q_i(k+1)[n]计算f_ik；

f_ik＝RMS(q_i(k+1)[n])+|EK(q_i(k+1)[n])|

其中，RMS(q_i(k+1)[n])和EK(q_i(k+1)[n])计算方法为：

其中，是q_i(k+1)[n]的算术平均值；

S3.7、判断当前筛分次数k是否满足k≥3，如果满足，则执行步骤S3.8，否则返回步骤S3.3；

S3.8、判断当前筛分次数k是否达到最大筛分次数I_max，如果达到，则IMF_i解调过程终止，FM_i[n]＝s_ik[n]，AM_i[n]＝IMF_i[n]/s_ik[n]；否则执行步骤S3.9；

S3.9、如果f_i(k-2)＜f_i(k-1)且f_i(k-1)＜f_ik，则IMF_i解调过程终止，此时>i[n]＝s_i(k-2)[n]，AM_i[n]＝IMF_i[n]/s_i(k-2)[n]；否则返回步骤S3.3；

S3.10、计算FM_i[n]信号的希尔伯特变换对则解析信号表示为：

S3.11、计算FM_i[n]的包络信号AM_HT[n]；

S3.12、更新FM_i[n]和AM_i[n]；

FM_i[n]＝FM_i[n]/AM_HT_i[n]

AM_i[n]＝AM_i[n]×AM_HT_i[n]

S3.13、计算FM_i[n]的瞬时相位φ_i[n]；

S3.14、计算FM_i[n]的瞬时频率f_i[n]；

S3.15、令i＝i+1，再返回步骤S3.2，直到所有的IMF解调结束；

S4、计算频率段特征

S4.1、通过步骤S3计算出的所有瞬时幅值AM_i[n]和瞬时频率f_i[n]得到改进型HHT的时频谱；

S4.2、设置分割深度j，j＝1,2,3,4,5，如图4所示，将时频谱在频域上按照分割深度j进行分割，得到一系列谱峭度其中，τ表示在同一分割深度上的频带分割位置；

S5、获得快速谱峭度图

S5.1、如图4所示，将步骤S4得到的一系列谱峭度构成快速谱峭度矩阵，再通过MATLAB中的imagesc函数将快速谱峭度矩阵生成快速谱峭度图；

S5.2、记录下快速谱峭度图中峭度最大的频带所在的中心频率和带宽；

S6、处理目标频带

S6.1、将步骤S5得到的中心频率和带宽作为滤波器的参数，然后用滤波器从原始信号中滤出目标信号；

S6.2、先对目标信号的包络进行平方，再做傅里叶变换，得到其平方包络谱；

S6.3、在平方包络谱上先找出轴承的故障特征频率点，然后后观测故障特征频率处的幅值变化，如果故障特征频率处的幅值基本保持不变，则判断滚动轴承正常；如果故障特征频率处的幅值变化显著，则判断滚动轴承故障；

其中，在平方包络谱上找出轴承的故障特征频率点的方法为：

设外圈固定，内圈随轴转动，D为轴承节径，d为滚动体直径，α为接触角，N为滚动体个数，f_r为转轴转频，轴承各部件的故障特征频率计算如下：

1)、内圈故障特征频率f_a：

2)、外圈故障特征频率f_o：

3)、滚动体故障特征频率f_b：

4)、保持架旋转频率f_c：

实例

在实施例中，依托四方所的轮轴轴承试验室中的轴承试验台，试验台如图5 所示。具体相关信息如下：

该轮对轴承故障诊断试验台由驱动电机、皮带传动系统、垂直加载装置、横向加载装置、两个风扇电机和控制系统组成。垂直和侧向载荷加载装置的设计用于模拟列车实际运行中轮对轴承承载的轴向和侧向的载荷。两个风扇电机可以产生与列车运行方向相反的风。通过两个加速度计确保轮对轴承水平方向和垂直方向的振动都能被检测到。我们选择了一个外圈故障数据，其具有 10×30mm的剥落面积，如图6所示。该数据是在236kN的垂直载荷以及604rpm 的速度下由垂直方向的加速度传感器以12800Hz的采样频率采集得到。计算得知外圈故障特征频率是f_BPFO＝137.7Hz。外圈故障的加速度数据的时域波形如图>

图8是Rilling’s HHT,Wu’s HHT,Wang’s HHT和本发明的改进型HHT滤波后的平方包络幅值谱。

为了证明本发明的有效性和优势，本发明选择了三个传统的HHT方法和本发明对比，分别是Rilling’s HHT,Wu’s HHT,Wang’s HHT，这三个方法的代码来自于网上公开的MATLAB代码。经过快速谱峭度图后，基于Rilling’s HHT的快速谱峭度图的最大峭度所对应的频带中心频率是300Hz,带宽为200Hz；基于 Wu’s HHT的快速谱峭度图的最大峭度所对应的频带中心频率是300Hz，带宽为 200Hz；基于Wang’s HHT的快速谱峭度图的最大峭度所对应的频带中心频率是 4666.67Hz，带宽为266.67Hz；基于本发明改进型HHT的快速谱峭度图的最大峭度所对应的频带中心频率是2300Hz，带宽为200Hz。然后使用每个方法得到的滤波器参数对原始信号滤波得到目标信号。如图8所示，四种方法基于快速谱峭度图得到的目标信号的平方包络幅值谱。显然，本发明改进型HHT方法经过快速谱峭度图得到的目标的平方包络谱在轴承外圈故障特征频率处比其他三种对比方法具有更大的幅值和更突出的冲击成分，证实了轴承的外圈存在故障。因此，本发明改进型HHT以及快速谱峭度图能够被成功地用于轴承的故障诊断中。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。