CN114166507A - 一种基于快速谱峭度的谐波识别方法 - Google Patents

Abstract

本发明提出的一种基于快速谱峭度的谐波识别方法，属于设备状态监测和故障诊断领域。首先对采集的振动信号进行快速谱峭度计算，得到对应谱峭度较大的p个滤波器，并对振动信号进行滤波和降噪处理，得到对应包络谱序列；再对每个包络谱序列进行区间划分，确定各子区间的基础频率；然后根据谐波搜索区间计算各子区间的有效谐波个数和谐波强度；最后计算各子区间的谐波强度指数，输出谐波强度指数最大的前n个子区间所对应的频率。该方法通过自动识别设备振动信号中的谐波成分，可以准确识别出调制信号的有效频率成分，为智能诊断轴承或齿轮等设备提供重要的故障征兆参数，相比人工识别调制信号的有效频率成分，明显简化工作复杂度且能提高诊断效率。

Description

一种基于快速谱峭度的谐波识别方法

技术领域

本发明涉及一种基于快速谱峭度的谐波识别方法，属于设备状态监测和故障诊断领域。

背景技术

在设备状态监测和故障诊断领域，轴承、齿轮等设备发生故障时振动信号伴随调制现象，频谱图出现多个边带驼峰。这些设备在出现故障时都与谐波存在对应关系，不同设备故障会对不同(频率)谐波造成影响，因此在对这类故障时振动信号信号进行调制解调时，往往展现出某种故障特征频率占主导的谐波现象。

通常，人工通过包络解调等方法逐一查看频谱图边带驼峰，再根据设备关键部件的特征频率逐个匹配包络谱图上的谐波成分，判断设备是否存在故障。但在监控设备数据较多时，这种人工识别谐波成分的方法将暴露出工作繁杂且诊断效率低等问题。

发明内容

本发明的目的是提供一种基于快速谱峭度的谐波识别方法，以解决目前人工识别过程中效率较低的问题。

本发明提出一种基于快速谱峭度的谐波识别方法，该方法包括以下步骤：

1)采集设备原始振动信号，通过短时傅里叶和滤波器组对设备原始振动信号进行快速谱峭度计算，选出前p个较大谱峭度对应的滤波器，其中各个滤波器包含对应的中心频率和带宽；

2)使用得到的p个滤波器对设备原始振动信号进行滤波和降噪处理，得到对应的p个包络谱序列；

3)对每个包络谱序列进行区间划分，得到若干个子区间，并确定每个子区间通频幅值最大的频率，将该频率作为该子区间的基础频率；

4)根据每个子区间的基础频率设定谐波搜索区间，对得到的每一个包络谱序列的各子区间按照设定的谐波搜索区间进行查找，确定每个子区间是否包含有有效谐波以及有效谐波个数和对应谐波的强度；

5)根据得到的每一个包络谱序列的各子区间的有效谐波个数和对应谐波强度，计算每一个包络谱序列中各子区间的谐波强度指数，输出谐波强度指数最大的前n个子区间的基础频率。

本发明将快速谱峭度应用到谐波识别中，通过筛选较大的快速谱峭度确定得到若干个滤波器，根据得到的滤波器对设备原始振动信号进行滤波生成若干个包络谱序列，对每个包络谱序列进行子区间划分，按照设定谐波搜索区间对个子区间进行搜索，实现谐波的自动查找。该方法通过自动识别设备振动信号中的谐波成分，可以准确识别出调制信号的有效频率成分，为智能诊断轴承或齿轮等设备提供重要的故障征兆参数；相比人工识别调制信号的有效频率成分，明显简化了工作复杂度且能提高诊断效率。

进一步地，所述步骤1)中快速谱峭度计算过程为：设计一个高通滤波器和一个低通滤波器，采用所述的高通滤波器和低通滤波器分别对设备原始振动信号进行滤波，对每个滤波结果做2倍降采样，不断反复迭代将滤波结果分成K层，并对K层的滤波结果进行谱峭度计算，计算公式为：

其中，

表示第K层第i个滤波器所得到的滤波结果，E(·)表示求取期望。

进一步地，为了快速实现对包络谱序列的降噪处理，所述步骤2)中降噪处理方法为：根据每一个滤波后的包络谱序列的最大幅值和中位数设定对应的噪声阈值，将每一个滤波后的包络谱序列中的频率幅值小于对应噪声阈值的频率幅值设为0。

进一步地，为了将时域信号转换成频域信号，所述降噪处理前需对滤波处理结果进行快速傅里叶变换和希尔伯特黄变换。

进一步地，所述步骤4)中若包络谱序列中满足有效谐波个数阈值的子区间个数小于n，则输出包络谱序列中通频幅值最大的前n个频率。

进一步地，为了保证查找谐波的准确性，所述步骤3)中的谐波搜索区间为[Cf-1.5，Cf+1.5]，其中f为基础频率，C为大于等于2的整数。

进一步地，所述步骤4)中谐波强度指数的计算公式为：

f(m,s)＝β₁s+β₂m

其中，f(m,s)为谐波强度指数，m为有效谐波个数，s_i(i＝1,2,…,m)为i个有效谐波对应的谐波强度，s为标准化处理后的谐波强度，β₁和β₂为指数系数，且β₁>β₂，H_zi(max)为第i个包络谱新序列中的最大通频幅值。

附图说明

图1为本发明基于快速谱峭度的谐波识别方法的整体流程图；

图2为本发明实施例中设备故障时的波形频谱图；

图3为本发明实施例的快速谱峭度图；

图4为本发明实施例中使用中心频率为640Hz和带宽为1280的滤波器得到的处理结果图；

图5为本发明实施例中使用中心频率为1600Hz和带宽为640的滤波器下得到的处理结果图；

图6为本发明实施例中使用中心频率为1493Hz和带宽为426的滤波器下得到的处理结果图；

图7为本发明实施例中使用中心频率为1440Hz和带宽为320的滤波器下得到的处理结果图；

图8为本发明实施例中使用中心频率为960Hz和带宽为213的滤波器滤下得到的处理结果图；

图9为本发明实施例中使用中心频率为1173Hz和带宽为213的滤波器下得到的处理结果图。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步地说明。

本发明提出的基于快速谱峭度的谐波识别方法，具体流程图如图1所示。首先采集设备原始振动信号，通过短时傅里叶和滤波器组对设备原始振动信号进行快速谱峭度计算，选出前p个较大谱峭度对应的滤波器；再使用得到的p个滤波器对设备原始振动信号进行滤波处理，得到对应的p个包络谱序列；然后将每个包络谱序列划分成若干个子区间，并确定每个子区间通频幅值最大的频率，将该频率作为该子区间的基础频率；之后根据每个子区间的基础频率设定谐波搜索区间，对得到的每一个包络谱序列的各子区间按照设定的谐波搜索区间进行查找，确定每个子区间是否包含有有效谐波以及有效谐波个数和对应谐波的强度；最后根据得到的每一个包络谱序列的各子区间的有效谐波个数和对应谐波强度，计算每一个包络谱序列中各子区间的谐波强度指数，输出谐波强度指数最大的前n个子区间的基础频率。该方法通过自动识别设备振动信号中的谐波成分，可以准确识别出调制信号的有效频率成分，为智能诊断轴承或齿轮等设备提供重要的故障征兆参数；相比人工识别调制信号的有效频率成分，明显简化了工作复杂度且能提高诊断效率。

步骤1.谱峭度计算及筛选最优滤波器

首先采集设备的原始振动信号，通过快速傅里叶和滤波器组对采集的设备振动信号进行快速谱峭度计算，选出前p个较大谱峭度对应的滤波器，各个滤波器包含对应的中心频率和带宽。在本实施例中采集的数据为某工厂扇风机电机轴承发生内圈故障时的振动信号，采样频率为f_s＝5120Hz，采样点数为N＝4096，转速为3000r/min，该轴承发生故障时的波形频谱图如图2所示，图2中上图为波形图，图2中下图为频谱图，从波形图上看到波形有明显的冲击成分，从频谱图上可以看出有多个驼峰带。

其中，快速谱峭度计算过程为：先设计一个高通滤波器和一个低通滤波器，采用设计的高通滤波器和低通滤波器分别对设备原始振动信号进行滤波，得到两个滤波结果，对每个滤波结果做2倍降采样，不断反复迭代将滤波结果分成K层，每个滤波结果包含中心频率和带宽，对K层滤波结果进行谱峭度计算。

其中高通滤波器和低通滤波器的计算公式为：

式中，h₁(t)为高通滤波器，h₀(t)为低通滤波器，h(t)为标准FIR低通滤波器，截止频率为1/8+ε(ε>0)。在本实施例中，截止频率设置为0.78，作为其他实施方式，ε的数值可以根据具体设备震动信号的频率确定。

其中谱峭度计算公式为：

式中，

表示第K层第i个滤波器所得到的滤波结果，E(·)表示求取期望。在本实施例中，K设置为4。作为其他实施方式，K的取值可以根据具体设备震动信号确定。

最后，整合所有计算出的谱峭度，构成设备原始振动信号的快速谱峭度图；将计算出的谱峭度从大到小排序，筛选出前p个较大的谱峭度对应的滤波器，其中各个滤波器包含对应的中心频率和带宽。在本实施例中，p设置为6。作为其他实施方式，p的取值可以根据谱峭度的具体数值确定。在本实施例中，整合所有计算出的谱峭度构成的设备原始振动信号快速谱峭度图如图3所示，从图3可以看出滤波结果最终被分为4层，图中的不同颜色深度分别代表着不同的峭度值，颜色越深的地方峭度值越大，可以看出最大谱峭度在层级为3.5，中心频率为1173Hz，带宽为213的地方，而频率1173Hz正好对应图2频谱图上最大谱峰的位置；将计算出的谱峭度从大到小排序，筛选出前6个谱峭度对应的滤波器，每个滤波器包含对应的中心频率和带宽，具体如表1所示：

表1：

序号	层级Level	中心频率fc/Hz	带宽Bw
				1	1	640	1280
2	2	1600	640
				3	2.5	1493	426
4	3	1440	320
				5	3.5	960	213
6	3.5	1173	213

步骤2.滤波和降噪处理

采用步骤1中筛选出的p个滤波器对设备原始振动信号进行滤波处理，为了将振动信号从时域转换到频域，本发明经过一系列快速傅里叶变换和希尔伯特黄变换等信号处理过程，最终生成对应p个包络谱序列。同时为了提高后续寻找谐波的效率，本发明对得到的p个包络谱序列分别进行降噪处理，具体过程为：首先，计算每个包络谱序列的最大幅值H_i(max)和中位数M_i，将包络谱序列每个频率的幅值表示为H_iv(f)，若H_iv(f)<αH_i(max)+M_i，则令H_iv(f)＝0；最终，得到降噪后的p个包络谱序列H_zi(f)(i＝1,2,…,p)。在本实施例中，α设置为0.2。作为其他实施方式，α的取值可根据包络谱序列最大幅值的数值确定。

在本实施例中，采用步骤1筛选出的6个滤波器对原始信号进行上述滤波和降噪处理，得到的每个滤波器处理结果如图4-图9所示。其中，图4为使用中心频率为640Hz和带宽为1280的滤波器得到的处理结果图，图5为使用中心频率为1600Hz和带宽为640的滤波器得到的处理结果图，图6为使用中心频率为1493Hz和带宽为426的滤波器得到的处理结果图，图7为使用中心频率为1440Hz和带宽为320的滤波器得到的处理结果图，图8为使用中心频率为960Hz和带宽为213的滤波器得到的处理结果图，图9为使用中心频率为1173Hz和带宽为213的滤波器得到的处理结果图；图4从上至下分别为原始波形信号图、包络信号图、包络谱图、降噪后的包络谱图；同理，图5-图9中每幅图从上至下均为对应滤波器下的原始波形信号图、包络信号图、包络谱图、降噪后的包络谱图；从图4-图9中的包络谱图可以看出已经将时域信号转变成频域信号，从图4-图9中降噪后的包络谱图可以看出在已经频率幅值小于设定阈值的频率幅值设为了0，可以减少计算量。

步骤3.包络谱序列子区间划分及基础频率确定

对步骤2得到的每个降噪后的包络谱序列进行区间划分，得到若干个子区间，并确定每个子区间通频幅值最大的频率，将该频率作为该子区间的基础频率f。在本实施例中，取步骤2降噪后的包络谱新序列H_zi(f)对应采样频率的前四分之一为搜索区间，将其等分k个子区间，查找各子区间内通频幅值最大的频率，将该频率作为基础频率。作为其他实施方式，搜索区间的大小和子区间个数可根据具体采样频率大小设定。

步骤4.查找有效谐波个数和谐波强度

在步骤3中每个降噪后的包络谱序列H_zi(f)划分的各个子区间内确定是否包含有有效谐波以及有效谐波个数和对应有效谐波的强度。首先，设定个子区间的谐波搜索区间，根据步骤3得到的各个子区间的基础频率f，对应的谐波频率为Cf(C≥2)，设定对应的谐波搜索区间为[Cf-1.5，Cf+1.5]，在本实施例中，C的取值区间为[2，20]；然后在各个子区间对应的谐波搜索区间内，寻找最大通频幅值A_max，若A_max大于设定阈值，则认为在该子区间内存在有效谐波，计算有效谐波个数及对应谐波强度，最后得到每个降噪后的包络谱序列H_zi(f)中各子区间的有效谐波个数m，本实施例中A_max设置为0.001。作为其他实施方式，C的取值范围和A_max取值大小可以根据具体情况确定。

以第i个降噪后的包络谱序列为例H_zi(f)为例，根据步骤2得到H_zi(f)划分后的各子区间的基础频率，设定谐波频率为2f时的各子区间对应的谐波搜索区间为[2f-1.5，2f+1.5]，在该区间下寻找最大通频幅值A_max，当A_max大于设定阈值，认为在该子区间内存在有效谐波，同理谐波频率为3f、4f、…等时的各子区间下的有效谐波判断方法相同，最后计算第i个降噪后的包络谱序列中各子区间下的有效谐波个数m及对应谐波强度s。在本实施例中，对图4-图9中的每一幅降噪后的包络谱图进行有效谐波识别，确定是否存在有效谐波及有效谐波个数和对应谐波强度。

步骤5.根据谐波强度指数筛选设备可能存在故障频率

根据步骤4计算得到的每个降噪后的包络谱序列H_zi(f)中各子区间的有效谐波个数m和谐波强度s，计算每一个降噪后包络谱序列中各子区间的谐波强度指数，若包络谱序列H_zi(f)中满足有效谐波个数阈值q的子区间个数大于等于n，则输出谐波强度指数最大的前n个子区间的基础频率；若包络谱序列H_zi(f)中满足有效谐波个数阈值q的子区间个数小于n，则输出包络谱序列H_zi(f)中通频幅值最大的前n个频率。本实施例中，q设置为2，n设置为5。作为其他实施方式，q和n的取值可以根据具体实际情况确定。以第i个降噪后的包络谱序列为例H_zi(f)为例，在步骤3中将H_zi(f)分了10个子区间，根据H_zi(f)中10个子区间的对应的有效谐波个数m和对应的谐波强度s，计算出H_zi(f)中10个子区间的谐波强度指数，若有效谐波个数大于等于2的子区间个数大于等于5个，则输出第i个包络谱序列H_zi(f)中谐波强度指数最大的前5个子区间的基础频率，若有效谐波个数大于等于2的子区间个数小于5个，则输出第i个包络谱序列H_zi(f)中通频幅值最大的前n个频率。

其中，谐波强度指数计算公式为：

f(m,s)＝β₁s+β₂m (3)

式中，f(m,s)为谐波强度指数，m为有效谐波个数，s_i(i＝1,2,…,m)为i个有效谐波对应的谐波强度，s为标准化处理后的谐波强度，H_zi(max)为降噪后第i个包络谱序列的最大通频幅值，β₁和β₂为指数系数，且β₁>β₂。本实施例中，β₁设置为1，β₂设置为0.05，β₁远大于β₂。作为其他实施方式，β₁和β₂的取值可根据实际情况确定。

在本实施例中，通过图4-图9中每一幅降噪后的包络谱图识别出的有效谐波个数和对应的谐波强度进行谐波强度指数计算，得到对应6个滤波器去噪后的包络谱图中识别出的谐波信息如表2所示。从图4-图9可以看出降噪后的包络谱图在50Hz处的幅值最大，且存在谐波，因为电机转速为3000r/min，转频正好为50Hz，此征兆正好是轴承内圈故障的重要特征参数；同时从表2可以得到频率为50Hz的谐波强度指数在每个滤波器中皆排在前三，与上述降噪后的包络谱图中信息对应。

表2：

Claims (7)

1.一种基于快速谱峭度的谐波识别方法，其特征在于，该方法包括以下步骤：

1)采集设备原始振动信号，通过短时傅里叶和滤波器组对设备原始振动信号进行快速谱峭度计算，选出前p个较大谱峭度对应的滤波器，其中各个滤波器包含对应的中心频率和带宽；

2)使用得到的p个滤波器对设备原始振动信号进行滤波和降噪处理，得到对应的p个包络谱序列；

3)对每个包络谱序列进行区间划分，得到若干个子区间，并确定每个子区间通频幅值最大的频率，将该频率作为该子区间的基础频率；

4)根据每个子区间的基础频率设定谐波搜索区间，对得到的每一个包络谱序列的各子区间按照设定的谐波搜索区间进行查找，确定每个子区间是否包含有有效谐波以及有效谐波个数和对应谐波的强度；

5)根据得到的每一个包络谱序列的各子区间的有效谐波个数和对应谐波强度，计算每一个包络谱序列中各子区间的谐波强度指数，输出谐波强度指数最大的前n个子区间的基础频率。

2.根据权利要求1所述的基于快速谱峭度的谐波识别方法，其特征在于，所述步骤1)中快速谱峭度计算过程为：

设计一个高通滤波器和一个低通滤波器，采用所述的高通滤波器和低通滤波器分别对设备原始振动信号进行滤波，对每个滤波结果做2倍降采样，不断反复迭代将滤波结果分成K层，并对K层的滤波结果进行谱峭度计算，计算公式为：

其中，

表示第K层第i个滤波器所得到的滤波结果，E(·)表示求取期望。

3.根据权利要求1所述的基于快速谱峭度的谐波识别方法，其特征在于，所述步骤2)中降噪处理方法为：根据每一个滤波后的包络谱序列的最大幅值和中位数设定对应的噪声阈值，将每一个滤波后的包络谱序列中的频率幅值小于对应噪声阈值的频率幅值设为0。

4.根据权利要求3所述的基于快速谱峭度的谐波识别方法，其特征在于，所述降噪处理前需对滤波处理结果进行快速傅里叶变换和希尔伯特黄变换。

5.根据权利要求1所述的基于快速谱峭度的谐波识别方法，其特征在于，所述步骤5)中若包络谱序列中满足有效谐波个数阈值的子区间个数小于n，则输出包络谱序列中通频幅值最大的前n个频率。

6.根据权利要求1所述的基于快速谱峭度的谐波识别方法，其特征在于，所述步骤3)中的谐波搜索区间为[Cf-1.5，Cf+1.5]，其中f为基础频率，C为大于等于2的整数。

7.根据权利要求1所述的基于快速谱峭度的谐波识别方法，其特征在于，所述步骤5)中谐波强度指数的计算公式为：

f(m,s)＝β₁s+β₂m

其中，f(m,s)为谐波强度指数，m为有效谐波个数，s_i(i＝1,2,…,m)为i个有效谐波对应的谐波强度，s为标准化处理后的谐波强度，β₁和β₂为指数系数，且β₁>β₂，H_zi(max)为第i个包络谱新序列中的最大通频幅值。

Images