CN110542546A - 基于谐波分形分解快速算法的机械故障特征识别方法 - Google Patents

Abstract

本发明涉及一种基于谐波分形分解快速算法的机械故障特征识别方法，对信号进行深度为的分形分解提出了基于频域的快速算法，实现了分形“频率‑尺度”网格上各尺度各子空间的时域信号获取，在算法的执行中只需要计算网格最后一层的小波包子空间分解，其它尺度的子空间都可以直接由这些子空间进行直接计算，避免了对每个尺度每个子空间进行的信号重构，减少计算强度，提交分析效率，还提出了一种强暂态冲击抑制指标，用于抑制信号中的偶发性冲击特征对分析结果的干扰，对分形分解结果进行最优特征评价，提高故障识别的准确性。

Description

基于谐波分形分解快速算法的机械故障特征识别方法

技术领域

本发明涉及机械故障诊断领域，具体涉及一种基于谐波分形分解快速算法的机械故障特征识别方法。

背景技术

旋转机械及以转子系统为主要组成部分的机械设备在长期的运行中容易发生部件的局部疲劳损伤，进而诱发机械故障。机械零部件的局部疲劳损伤在振动测试信号中往往表现为周期性的冲击特征。然而，由于机械设备其它部件振动成分的影响、其它设备振动经由地基传达到被测试设备，以及振动测试中可能出现的其它随机噪声的影响，无法从振动信号中直接识别故障特征成分。

针对此类周期性冲击特征成分，出现了许多基于先进信号处理方法的识别和提取技术，主要依赖于信号分解及特征优选。丁康等的“一种谱峭度和Morlet小波的滚动轴承微弱故障诊断方法[J]，振动工程学报，2014，27(1):128-135”，采用Morlet小波对振动信号进行分解，结合峭度统计指标进行特征优选，但这种方法容易误将信号中偶发性的强冲击作为故障特征进行错误的识别；王晓龙等的“自适应可调品质因子小波变换在轴承早期故障诊断中的应用[J]，航空动力学报，2017，32(10):2467-2475”采用参数化的小波变换对振动信号进行精细的分解，但这种自适应方法在频域上进行，参数调整的过程需要进行庞大的计算；陈彬强等“机械故障诊断的衍生增强离散解析小波分析框架[J]，机械工程学报.2014，50(17):77-86”提出的衍生增强离散解析小波分析框架，提出了信号分形的“频率-尺度”构型，但在每个尺度的每个子空间都必须进行单枝重构，分析效率较低。

发明内容

针对现有技术的缺陷，本发明的目的在于提供一种基于谐波分形分解快速算法的机械故障特征识别方法，能够避免对每个尺度每个子空间进行的信号重构，提出的强暂态冲击抑制指标能够抑制信号中偶发性的强冲击对分析结果的干扰，具有分析效率高，计算强度低，识别准确的特点，具有较好的实用性及工程应用推广价值。

为了达到上述目的，本发明所采用的技术方案是：

基于谐波分形分解快速算法的机械故障特征识别方法，包括以下步骤：

步骤1)、在装备机械传动链的轴承座上安装振动传感器以采集动态信号，信号的采样频率为f_s，采样长度为N，N必须为偶数，对信号进行去均值处理，得到：

x:＝{x(n)|n＝1,2,…,N}

根据监测设备的机械结构计算系统可能出现的若干故障特征频率且i≤M},其中表示正整数，M表示故障特征频率的总个数；

步骤2)、对动态信号x进行J层谐波分形分解，J≥2，具体步骤为：

步骤2.1).对x进行快速傅里叶变换，得到频谱函数

步骤2.2).将原信号分解为2^J个小波包子空间wp_J:＝{wp_J,k(n)|k＝1,2,…,2^J且k＝0,1,2,…,N-1}，wp_J中的各元素按照频带的能量重心从小到大进行排列，第k个小波包子空间wp_J,k(n)所占据的频带为

步骤2.3).从原信号的频谱构建小波包子空间wp_J,k(n)的频谱函数定义为

步骤2.4).对进行快速傅里叶逆变换得到对应的时域信号{wp_J,k(n)}表示为：

其中IFFT{·}表示傅里叶逆变换算子；

步骤3)、采用谐波分形快速分析算法对输入信号进行多尺度分解；当分析深度为J时，最终生成的二维“频率-尺度”划分平面由2J个尺度构成，分别标记为{'1e','1p','2p','2e',…,'Je','Jp'}，其中标记为‘ip’(i＝1,2,…J)的尺度为常规二进尺度；标记为‘ie’(i＝1,2,…J)的尺度为间接联合尺度，每个尺度由若干个子空间构成；各尺度各子空间的计算步骤为：

步骤3.1).按照步骤2)对输入信号进行层数为J的谐波分形分解，得到标记为‘Jp’的常规二进尺度子空间信号wp_J＝{wp_J,k|k＝1,2,…,2^J}

步骤3.2).计算其它常规二进尺度的子空间信号：

标记为'jp'的第j层常规二进尺度拥有2^j个子空间{wp_j,k(n)|k＝1,2,…,2^j且n＝0,1,2,…,N-1},处于该尺度的子空间wp_j,k(n)可由相邻较高层的常规二进尺度子空间wp_j+1＝{wp_j+1,k|k＝1,2,…,2^j+1}直接构建，表示为

wp_j,k(n)＝wp_j+1,2k-1(n)+wp_j+1,2k(n)；

因此可以由最高层(第J层)的常规二进尺度wp_J中的子空间元素构建wp_J-1中的各元素，进而不断迭代间接计算其它尺度(j＝J-2,J-3,…,1)的子空间信号；

步骤3.3).计算间接联合尺度的子空间信号：

标记为'je'的第j层间接联合尺度拥有2^j-1个子空间ewp_j:＝{ewp_j,k|k＝1,…,2^{j -1}}。其中的元素ewp_j,k占据的频带为：

其时域信号由由同一层的常规二进尺度子空间ewp_j,k(n)计算：

步骤4)、对各尺度的子空间信号定义一种“强暂态冲击抑制指标”；

步骤5)、采用所提出的“强暂态冲击抑制指标”对分析深度为J的分形分解结果进行最优特征评价，选取该指标值最大的子空间，通过观察该最优子空间的时域、频域及包络解调域特征确认其中是否存在频率与f_c(i)相接近的故障特征。

进一步地，所述“强暂态冲击抑制指标”的计算步骤为：

步骤4.1).输入为子空间时域信号w(n)，依次选取故障特征频率集合中的{f_c(i),1≤i≤M}一个元素f_c(i)；

步骤4.2)、将输入信号w(n)从第一个取样数据点开始按照间隔L＝int{f_s/f_c(i)}划分为首尾相接且不重叠的int{N/L}个时间分区{seg_i|i＝1,2,…,m}，算子int{·}为取整函数，输出不大于输入实数的最大整数，必须保证L是偶数，即当L是奇数时，令L＝L-1。舍弃长度超出int{N/L}·L的取样点，将w(n)截断后的信号定义为

步骤4.3)、依次从{seg_l|l＝1,2,…,m}中抽取第k个元素，将剩余(m-1)个时域分块的采样作为一个信号，再计算其方差kurt_k；

步骤4.4)、计算方差集合{kurt_k|k＝1,2,…,m}的平均值mean，如果满足判别式(1)

令强暂态冲击抑制指标imp(i)＝0，结束计算并进入步骤4.5)；如果不满足上述判别式(1)，令强暂态冲击抑制指标imp(i)＝kurtosis(w)并进入步骤4.5)；

步骤4.5).针对步骤4.2)中得到的信号删去最左端及最右端各int{L/2}个取样数据，剩余的信号称之为将从第一个取样数据点开始按照间隔L将划分为首尾相接且不重叠的(m-1)个时间分区

步骤4.6).依次从中抽取第k个元素，将剩余(m-2)个时域分块的采样作为一个信号，再计算其方差

步骤4.7).计算方差集合的平均值如果满足判别式(2)：

令强暂态冲击抑制指标结束计算并进入步骤4.8)；如果不满足上述判别式(2)，令强暂态冲击抑制指标结束计算并进入步骤4.8)；

步骤4.8).将原始信号w(n)对应于故障特征频率f_c(i)的“强暂态冲击抑制指标”定义为

步骤4.9).根据步骤4.1)～步骤4.8)分别计算各故障特征频率的“强暂态冲击抑制指标”，并将子空间w(n)的“强暂态冲击抑制指标”定义为

有益效果：

1、一方面，本发明公开的一种基于谐波分形分解快速算法的机械故障特征识别方法，对信号进行深度为J的分形分解提出了基于频域的快速算法，实现了分形“频率-尺度”网格上各尺度各子空间的时域信号获取，在算法的执行中只需要计算网格最后一层的小波包子空间分解，其它尺度的子空间都可以直接由这些子空间进行直接计算，避免了对每个尺度每个子空间进行的信号重构，减少计算强度，提交分析效率。

2、另一方面，本发明针对振动测试中的偶发性强冲击在时域上的稀疏分布提出了一种强暂态冲击抑制指标，用于抑制信号中的偶发性冲击特征对分析结果的干扰，该指标考虑了偶发性冲击在时域上的稀疏分布特性，通过时域分块及选择性遮掩计算信号的峭度值，根据峭度值的波动性判断该尺度是干扰成分还是故障特征成分，最后将具有强暂态冲击抑制指标的子空间作为最优分析特征，从其衍生出的时域、频域、包络域进行故障确认，相比现有技术，本发明能够有效抑制信号中的偶发性冲击特征对分析结果的干扰，提高故障识别的准确性。

附图说明

图1为本发明采集振动信号的时域波形图；

图2为图1中所采集振动信号的傅里叶频谱图；

图3为分析深度为J＝3常规二进尺度分解所产生8个子空间的频带划分方式；

图4对图1中振动信号谐波分形分解结果的强暂态冲击抑制指标分布图；

图5为所提取的最优特征结果：其中(a)时域波形图；(b)包络解调谱。

具体实施方式

下面结合附图和具体实施例对本发明做进一步详细的说明。

本发明采用机械设备振动信号作为诊断分析媒介，通过信号分析方法挖掘其中的故障特征，需要机械系统潜在故障特征频率作为先验知识，本发明的基于谐波分形分解快速算法的机械故障特征识别方法具体包括以下步骤。

步骤1)、在装备机械传动链的轴承座上安装振动传感器以采集动态信号，如图1所示，信号的采样频率为f_s,采样长度为N(N必须为偶数),对信号进行去均值处理，得到

x:＝{x(n)|n＝1,2,…,N}；

根据监测设备的机械结构计算系统可能出现的若干故障特征频率且i≤M},其中表示正整数，M表示故障特征频率的总个数；

步骤2)、对动态信号x进行J层谐波分形分解(J≥2)，其具体步骤为：

2.1)对x进行快速傅里叶变换，得到频谱频谱图如图2所示。

2.2)将原信号分解为2^J个小波包子空间wp_J:＝{wp_J,k(n)|k＝1,2,…,2^J且k＝0,1,2,…,N-1}。wp_J中的各元素按照频带的能量重心从小到大进行排列，第k个子空间wp_J,k(n)所占据的频带为

常规二进尺度的子空间频带划分方式如图3所示。

2.3)从原信号的频谱构建小波包子空间wp_J,k(n)的频谱函数定义为

2.4)对进行快速傅里叶逆变换得到对应的时域信号{wp_J,k(n)}，表示为：

其中IFFT{·}表示傅里叶逆变换算子。

步骤3)、采用谐波分形快速分析算法对输入信号进行多尺度分解。当分析深度为J时，最终生成的二维“频率-尺度”划分平面由2J个尺度构成，分别标记为{'1e','1p','2p','2e',…,'Je','Jp'}。其中标记为‘ip’(i＝1,2,…J)的尺度为常规二进尺度；标记为‘ie’(i＝1,2,…J)的尺度为间接联合尺度。每个尺度由若干个子空间构成。各尺度各子空间的计算步骤为：

3.1)按照步骤2)对输入信号进行层数为J的谐波分形分解，得到标记为‘Jp’尺度的常规二进子空间wp_J＝{wp_J,k|k＝1,2,…,2^J}；

3.2)计算其它常规二进尺度‘jp’(j＝1,2,…J-1)的子空间信号；

标记为'jp'的第j层常规二进尺度拥有2^j个子空间wp_j＝{wp_j,k(n)|k＝1,2,…,2^j},处于该尺度的子空间wp_j,k(n)可由相邻较高层尺度的子空间wp_j+1＝{wp_j+1,k|k＝1,2,…,2^j+1}直接构建,表示为：

wp_j,k(n)＝wp_j+1,2k-1(n)+wp_j+1,2k(n)。

因此可以由最高层(第J层)的常规二进尺度wp_J中的子空间元素构建wp_J-1中的各元素，进而不断迭代间接计算其它尺度(j＝J-2,J-3,…,1)的子空间信号；

3.3)计算间接联合小波包尺度的子空间信号

标记为'je'的第j层间接联合尺度拥有2^j-1个子空间ewp_j:＝{ewp_j,k|k＝1,…,2^{j -1}}。其中的元素ewp_j,k占据的频带为：

其时域信号可由同一层的常规二进尺度子空间ewp_j,k(n)计算：

步骤4)对各尺度的子空间信号定义一种“强暂态冲击抑制指标”，采用“强暂态冲击抑制指标”对分析深度为J的分形分解结果进行最优特征评价，“强暂态冲击抑制指标”的计算步骤为：

4.1).输入为子空间时域信号w(n)，依次选取故障特征频率集合中的{f_c(i),1≤i≤M}一个元素f_c(i)。

4.2).将输入信号w(n)从第一个取样数据点开始按照间隔L＝int{f_s/f_c(i)}划分为首尾相接且不重叠的int{N/L}个时间分区{seg_i|i＝1,2,…,m}。算子int{·}为取整函数，输出不大于输入实数的最大整数。必须保证L是偶数，即当L是奇数时，令L＝L-1。舍弃长度超出int{N/L}·L的取样点。将w(n)截断后的信号定义为

4.3).依次从{seg_i|i＝1,2,…,m}中抽取第k个元素，将剩余(m-1)个时域分块的采样作为一个信号，再计算其方差kurt_k；

4.4).计算方差集合{kurt_k|k＝1,2,…,m}的平均值mean,如果满足判别式(1)：

令强暂态冲击抑制指标imp(i)＝0，结束计算并进入步骤4.5)；如果不满足上述判别式(1)，令强暂态冲击抑制指标imp(i)＝kurtosis(w)并进入步骤4.5)。

4.5).针对步骤4.2)中得到的信号删去最左端及最右端各int{L/2}个取样数据，剩余的信号称之为将从第一个取样数据点开始按照间隔L将划分为首尾相接且不重叠的(m-1)个时间分区

4.6).依次从中抽取第k个元素，将剩余(m-2)个时域分块的采样作为一个信号，再计算其方差

4.7).计算方差集合的平均值如果满足判别式(2)：

令强暂态冲击抑制指标结束计算并进入步骤4.8)。如果不满足上述判别式(2)，令强暂态冲击抑制指标结束计算并进入步骤4.8)。

4.8).将w(n)对应于故障特征频率f_c(i)的“强暂态冲击抑制指标”定义为

4.9).根据步骤4.1)～步骤4.8)分别计算各故障特征频率的“强暂态冲击抑制指标”，并将子空间w(n)的“强暂态冲击抑制指标”定义为

步骤5)、采用所提出的“强暂态冲击抑制指标”对分析深度为J的分形分解结果进行最优特征评价分析。根据步骤4)计算各尺度上子空间的强暂态冲击抑制指标，常规二进尺度的表示为间接联合尺度的表示为将这些计算结果绘制在“频率-尺度”划分网格上，如图4所示。将该指标值最大的子空间定义为最优分析子空间，通过观察该最优分析子空间的时域、频域及包络解调域特征确认其中是否存在频率与f_c(i)相接近的故障特征，如图5所示。

下面结合具体的工程例对图1～5的内容进行进一步说明。

本工程例用于验证本发明内容的有效性和正确性。振动位移信号采集自本特利转子实验台，信号采样频率2000Hz，采样长度为1024。该试验台由直流电机驱动转轴运动(转速为2200r/min)，转轴上安装有两个质量块，采样摩擦螺栓在转轴的末端模拟动静碰摩故障，记录下的信号时域波形及其傅里叶频谱如图1及图2所示。

为了对信号进行深入的分析，首先对原始信号进行深度J的常规二进尺度分解，在wp_J上得到2^J个子空间，当J＝3时，8个子空间的频带划分方式如图3所示。从wp_J出发，由相邻常规二进尺度的迭代递推公式可以间接计算wp_J-1；再由同一层常规二进尺度及间接联合尺度的公式可以计算ewp_J；由wp_J-1可以间接计算wp_J-2及ewp_J-1。如此往复迭代，最后由wp₁计算ewp₁，结束谐波分形分解。

然后对各尺度上的子空间计算“强暂态冲击抑制指标”。该转子系统的故障特征频率只有一个，为f_c(1)＝36.67Hz。则在计算imp(1)时，对各子空间按照L＝54截取前18个时间分块并舍弃最后52个数据采样得到在计算时舍弃首位各27个数据采样并按照L＝54截取前17个时间分块；最终动态信号谐波分形分解各尺度的强暂态冲击抑制指标分布图如图4所示，其中wp_3,2(n)的强暂态冲击抑制指标数值最大，选为最优子空间。

最优子空间的时域波形及希尔伯特包络解调谱如图5所示，两个波形都出现了频率为36.67Hz的周期性周期，说明本发明提出的方法能够有效地提取机械故障特征。

综上，本发明的一种基于谐波分形分解快速算法的机械故障特征识别方法，对信号进行深度为J的分形分解提出了基于频域的快速算法，首先将信号的频谱分为2^J个频带首尾相接且互不重叠的小波包，通过傅里叶逆变换获得各子空间的时域信号；其它小波包子空间和联合小波包子空间可由这2^J个子空间组合得到，避免了对每个尺度每个子空间进行的信号重构。

此外，为了抑制信号中的偶发性冲击特征对分析结果的干扰，提出一种强暂态冲击抑制指标，该指标考虑了偶发性冲击在时域上的稀疏分布特性，通过时域分块及选择性遮掩计算信号的峭度值，根据峭度值的波动性判断该尺度是干扰成分还是故障特征成分，最后将具有强暂态冲击抑制指标的子空间作为最优分析特征，从起衍生出的时域、频域、包络域进行故障确认。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (2)

1.基于谐波分形分解快速算法的机械故障特征识别方法，其特征在于，包括以下步骤：

步骤1)、在装备机械传动链的轴承座上安装振动传感器以采集动态信号，信号的采样频率为f_s，采样长度为N，N必须为偶数，对信号进行去均值处理，得到：

x:＝{x(n)|n＝1,2,…,N}

根据监测设备的机械结构计算系统可能出现的若干故障特征频率其中表示正整数，M表示故障特征频率的总个数；

步骤2)、对动态信号x进行J层谐波分形分解，J≥2，具体步骤为：

步骤2.1).对x进行快速傅里叶变换，得到频谱函数

步骤2.2).将原信号分解为2^J个小波包子空间：wp_J:＝{wp_J,k(n)|k＝1,2,…,2^J且k＝0,1,2,…,N-1}，wp_J中的各元素按照频带的能量重心从小到大进行排列，第k个小波包子空间wp_J,k(n)所占据的频带为：

步骤2.3).从原信号的频谱构建小波包子空间wp_J,k(n)的频谱函数定义为

步骤2.4).对进行快速傅里叶逆变换得到对应的时域信号{wp_J,k(n)}表示为：

其中IFFT{·}表示傅里叶逆变换算子；

步骤3)、采用谐波分形快速分析算法对输入信号进行多尺度分解；当分析深度为J时，最终生成的二维“频率-尺度”划分平面由2J个尺度构成，分别标记为{'1e','1p','2p','2e',…,'Je','Jp'}，其中标记为‘ip’(i＝1,2,…J)的尺度为常规二进尺度；标记为‘ie’(i＝1,2,…J)的尺度为间接联合尺度，每个尺度由若干个子空间构成，各尺度各子空间的计算步骤为：

步骤3.1).按照步骤2)对输入信号进行层数为J的谐波分形分解，得到标记为‘Jp’的常规二进尺度子空间信号wp_J＝{wp_J,k|k＝1,2,…,2^J}；

步骤3.2).计算其它常规二进尺度的子空间信号；

标记为'jp'的第j层常规二进尺度拥有2^j个子空间

{wp_j,k(n)|k＝1,2,…,2^j且n＝0,1,2,…,N-1}，处于该尺度的子空间wp_j,k(n)可由相邻较高层的常规二进尺度子空间wp_j+1＝{wp_j+1,k|k＝1,2,…,2^j+1}直接构建，表示为：

wp_j,k(n)＝wp_j+1,2k-1(n)+wp_j+1,2k(n)；

因此可以由最高层即第J层的常规二进尺度wp_J中的子空间元素构建wp_J-1中的各元素，进而不断迭代间接计算其它尺度j＝J-2,J-3,…,1的子空间信号；

步骤3.3).计算间接联合尺度的子空间信号；

标记为'je'的第j层间接联合尺度拥有2^j-1个子空间ewp_j:＝{ewp_j,k|k＝1,…,2^j-1}，其中的元素ewp_j,k占据的频带为：

其时域信号由由同一层的常规二进尺度子空间ewp_j,k(n)计算：

步骤4)、对各尺度的子空间信号定义一种“强暂态冲击抑制指标”；

步骤5)、采用所提出的“强暂态冲击抑制指标”对分析深度为J的分形分解结果进行最优特征评价，选取该指标值最大的子空间，通过观察该最优子空间的时域、频域及包络解调域特征确认其中是否存在频率与f_c(i)相接近的故障特征。

2.如权利要求1所述的机械故障特征识别方法，其特征在于，所述“强暂态冲击抑制指标”的计算步骤为：

步骤4.1).输入为子空间时域信号w(n)，依次选取故障特征频率集合中的{f_c(i),1≤i≤M}一个元素f_c(i)；

步骤4.2)、将输入信号w(n)从第一个取样数据点开始按照间隔L＝int{f_s/f_c(i)}划分为首尾相接且不重叠的int{N/L}个时间分区{seg_i|i＝1,2,…,m}，算子int{·}为取整函数，输出不大于输入实数的最大整数；必须保证L是偶数，即当L是奇数时，令L＝L-1，舍弃长度超出int{N/L}·L的取样点，将w(n)截断后的信号定义为

步骤4.3)、依次从{seg_l|l＝1,2,…,m}中抽取第k个元素，将剩余(m-1)个时域分块的采样作为一个信号，再计算其方差kurt_k；

步骤4.4)、计算方差集合{kurt_k|k＝1,2,…,m}的平均值mean，如果满足判别式(1)

令强暂态冲击抑制指标imp(i)＝0，结束计算并进入步骤4.5)；如果不满足上述判别式(1)的条件，令强暂态冲击抑制指标imp(i)＝kurtosis(w)并进入步骤4.5)；

步骤4.5)、针对步骤4.2)中得到的信号删去最左端及最右端各int{L/2}个取样数据，剩余的信号称之为将从第一个取样数据点开始按照间隔L将划分为首尾相接且不重叠的(m-1)个时间分区

步骤4.6)、依次从中抽取第k个元素，将剩余(m-2)个时域分块的采样作为一个信号，再计算其方差

步骤4.7)、计算方差集合的平均值如果满足判别式(2)：

令强暂态冲击抑制指标结束计算并进入步骤4.8)；如果不满足上述判别式(2)，令强暂态冲击抑制指标结束计算并进入步骤4.8)；

步骤4.8)、将原始信号w(n)对应于故障特征频率f_c(i)的“强暂态冲击抑制指标”定义为

步骤4.9)、根据步骤4.1)～步骤4.8)分别计算各故障特征频率的“强暂态冲击抑制指标”，并将子空间w(n)的“强暂态冲击抑制指标”定义为