CN107808114B - 一种基于信号时频分解的幅值谱峭度图的实现方法 - Google Patents

Abstract

本发明公开一种基于信号时频分解的幅值谱峭度图的实现方法，该方法首先用频率切片小波变换对振动信号时频分解获得其时频分解矩阵；然后，在每一尺度分别抽取该尺度每个给定的子时频空间对应的子矩阵，按时间方向求取各个子矩阵的时域均值，把该时域均值作为该子时频空间信号分量幅值随时间变化的趋势，用该时域均值的幅值平方谱峭度按照尺度—频带排列，得到信号的谱峭度尺度—频带平面图，即时频幅值谱峭度图，用颜色深度表示谱峭度的大小，峭度值越大对应谱平面区域的颜色越深，颜色最深的区域对应的频带为特征频带；可用于提取故障特征。

Description

一种基于信号时频分解的幅值谱峭度图的实现方法

技术领域

本发明属于机械设备信号处理领域，具体涉及一种基于信号时频分解的幅值谱峭度图的实现方法。

背景技术

基于峭度色度图确定特征频带提取损伤特征频率是轴承、齿轮等零部件故障诊断的方法之一。文献[1]提出一种谱峭度图构造方法，采用一对互补的滤波器组，对原始信号采用二进塔式分解，在每一分解尺度得到一系列分量，再求这些分量的峭度，然后把每个尺度的各个峭度按频带映射到尺度—频带平面图上；为了改善上述方法的尺度分辨率，文献[1]还提出一种滤波器组的1/3二进分解树的改进方法，在更多的尺度上分解以得到更加细致的子频带划分，这些子频带的信号分量改善了二进塔式分解的尺度—频带平面图分辨率；文献[2]用互补的滤波器组分解得到分量信号的包络幅值平方的频谱峭度指标代替文献[1]的峭度指标，构造另一种称之为Inforgram的尺度—频带平面图，用来检测信号中的重复冲击信号特征。文献[3]采用小波包分解方法构造谱峭度图，利用各个尺度上小波包分解得到的小波包分量来求峭度，再按小波包频带进行映射，得到尺度—频带平面图。

上述方法改进了文献[4]的计算效率，在应用过程中存在以下问题：

(1)采用滤波器组方法，需要确定互补的高通和低通滤波器参数及滤波器的长度，需要一定的先验知识。

(2)小波包分解方法在分解过程中的下抽样操作，使得(第3尺度分解以后)小波包对应的频带发生倒置，在形成尺度—频带平面图时，需要对小波包分量对应的频带进行重排。另外，小波包采用二进塔式分解方法，尺度分辨率较低。

(3)不管是上述哪一种方法，因为采用互补滤器组分解，互补滤器组的特性确定了频带分割范围，在每一尺度分解信号时，高通滤波器和低通滤波器把信号分解到固定的频带范围内，因此，往往不能根据需求确定所需要精确频带以获取有效的分量。

以下是申请人检索的与本申请相关的参考文献：

[1]J.Antoni.Fast computation of the kurtogram for the detection oftransient faults[J]，Mechanical Systems and Signal Processing,21:108-124，2007。

[2]J Antoni.The infogram:Entropic evidence of the signature ofrepetitive transients[J].Mechanical Systems and Signal Processing，74:73-94，2016。

[3]Yaguo Lei，Jing Lin，Zhengjia He，Yanyang Zi.Application of animproved kurtogram method for fault diagnosis of rolling element bearings[J].Mechanical Systems and Signal Processing，25:1738-1749，2011。

[4]J.Antoni.The spectral kurtosis:a useful tool for characterizingnon-stationary signals[J].Mechanical Systems and Signal Processing，20:282-307，2006。

发明内容

针对目前谱峭度图构造方法的不足，本发明的目的在于提供一种基于信号时频分解的幅值谱峭度图的实现方法。

为了实现上述任务，本发明采用如下的技术解决方案：

一种基于信号时频分解的幅值谱峭度图的实现方法，其特征在于，该方法首先用频率切片小波变换对振动信号时频分解获得其时频分解矩阵；然后，在每一尺度分别抽取该尺度每个给定的子时频空间对应的子矩阵，按时间方向求取各个子矩阵的时域均值，把该时域均值作为该子时频空间信号分量幅值随时间变化的趋势，用该时域均值的幅值平方谱峭度按照尺度—频带排列，得到信号的谱峭度尺度—频带平面图，即时频幅值谱峭度图，用颜色深度表示谱峭度的大小，峭度值越大对应谱平面区域的颜色越深，颜色最深的区域对应的频带为特征频带；其中：

所述的频率切片小波变换是一种时频分解方法，振动信号变换之前，需要选择频率切片函数，变换结果为二维时频矩阵；

所述时域均值是按照矩阵的行或列求所有行元素或列元素的平均值。

具体实现方法如下：

(1)根据故障特征频率和振动信号的采样频率确定峭度谱的尺度以及子时频空间的个数，即最大尺度对应的子时频空间个数；

(2)去除采集的振动信号中的直流分量；

(3)选择频率切片函数，对去除直流分量的信号做频率切片小波变换得到其时频分解矩阵；

(4)由时频分解矩阵求出信号的时频幅值矩阵；

(5)在每一个尺度，求出每个子时频空间对应的时频幅值子矩阵的时域幅值平均向量，再求该向量平方的频谱，然后求出该频谱的幅值峭度，每一尺度的所有子时频空间的幅值峭度构成谱峭度向量；

(6)由每个尺度的谱峭度向量构造谱峭度矩阵，再将谱峭度矩阵映射到尺度—频率平面图，以色度深浅来表示谱峭度的大小，得到信号的谱峭度图；

(7)信号分析时，在谱峭度图上颜色最深的区域为特征频带，用频率切片小波逆变换分离出该特征频带的信号分量，以备进一步分析。

本发明的基于信号时频分解的幅值谱峭度图的实现方法，带来的有益效果是在于，采用基于信号时频分解构造谱峭度图，通过提取子时频空间上的信号幅值随时间变化的趋势来计算等效的谱峭度，具有较高的计算效率，可以用于机械设备、电气系统的信号分析，尤其是对具有调频调幅特征的信号，该方法可以直观地图示出调频调幅分量的最佳解调频带，有益于信号的特征提取。

附图说明

图1是信号s(t)的L尺度频带分割示意图；

图2是峭度色谱排列形式；

图3为一个轴承内圈损伤的振动加速度信号及其频谱图，其中图(1)为振动信号的时域波形图，图(2)为振动信号的频谱图；

图4是幅值峭度谱图；

图5是特征子频带重构信号时域波形图；

图6是特征子频带重构信号的包络谱图。

以下结合附图和实施例对本发明作进一步的详细说明。

具体实施方式

本实施例给出一种基于信号时频分解的幅值谱峭度图的实现方法，利用频率小波变换的全频带时频连续分解得到时频分解矩阵，在给定尺度把信号的时频分解空间以等宽频带的方式进行分割，得到一系列子时频空间的时域特征分量，以该特征分量为基础构造谱峭度图。

首先用频率切片小波变换对振动信号时频分解获得其时频分解矩阵；然后，在每一尺度分别抽取该尺度每个给定的子时频空间对应的子矩阵，按时间方向求取各个子矩阵的时域均值，把该时域均值作为该子时频空间信号分量幅值随时间变化的趋势，用该时域均值的幅值平方谱峭度按照尺度—频带排列，得到信号的谱峭度尺度—频带平面图，即时频幅值谱峭度图，用颜色深度表示谱峭度的大小，峭度值越大对应谱平面区域的颜色越深，颜色最深的区域对应的频带为特征频带；其中：

所述的频率切片小波变换是一种时频分解方法，振动信号变换之前，需要选择频率切片函数，变换结果为二维时频矩阵；

所述时域均值是按照矩阵的行或列求所有行元素或列元素的平均值。

本实施例中，上述幅值平方谱峭度是描述信号的一种无量纲指标，谱峭度尺度—频带平面图是描述峭度在尺度—频率平面上的色度图。

具体实现方法如下：

(1)根据故障特征频率和振动信号的采样频率确定峭度谱的尺度以及子时频空间的个数，即最大尺度对应的子时频空间个数；

(2)去除采集的振动信号中的直流分量；

(3)选择频率切片函数，对去除直流分量的信号做频率切片小波变换得到其时频分解矩阵；

(4)由时频分解矩阵求出信号的时频幅值矩阵；

(5)在每一个尺度，求出每个子时频空间对应的时频幅值子矩阵的时域幅值平均向量，再求该向量平方的频谱，然后求出该频谱的幅值峭度，每一尺度的所有子时频空间的幅值峭度构成谱峭度向量；

(6)由每个尺度的谱峭度向量构造谱峭度矩阵，再将谱峭度矩阵映射到尺度—频率平面图，以色度深浅来表示谱峭度的大小，得到信号的谱峭度图；

(7)信号分析时，在谱峭度图上颜色最深的区域为特征频带，用频率切片小波逆变换分离出该特征频带的信号分量，以备进一步分析。

以下给出具体的实现过程。

设p(t)为频率切片函数，时域信号s(t)在频带[f_bgn,f_end]的频率切片小波变换时频分解为

式中κ＞0，

为

的共轭函数,

是s(t)傅里叶变换。

设κ为任意给定值时，由于ω＝2πf，把信号s(t)的频率切片小波变换时频分解W(t,ω,κ)写为W_κ(t,f)。

设s(t)为振动信号，去除振动信号中的直流分量：

其中

为s(t)在时间区间[0,t_end]的均值。

设采样频率为f_s，采用公式(1)对信号x(t)进行频率切片小波变换，得到x(t)在时频区间[0,t_end,0,f_s/2]的时频分解矩阵为W_κ(t,f)＝{W(t_k,f_k),t_k＝0～t_end,f_k＝0～f_s/2}。

设故障特征频率为f_k，定义最小分解带宽为：Δ_BW＝(3～5)f_k,则定义谱的最大尺度为：

其中，f_N＝f_s/2，f_N为奈奎斯特频率。

在每一谱尺度，对信号的奈奎斯特频带做等宽频带分割，尺度带宽为：

其中，L_i＝1,2,...,L。

信号s(t)的L尺度频带分割示意图如图1所示。

则L_i尺度的第j个时频区间的频带为：(j-1)Δ_i～jΔ_i，j＝1,2,...,L_i。

时频区间[0,t_end,(j-1)Δ_i,jΔ_i]的时域幅值平均值为：

令

对式(7)做傅里叶变换

式中ω＝2πf。

设

为|FA_i,j(ω)|的均值，σ_A为|FA_i,j(ω)|的均方差。|FA_i,j(ω)|的峭度为

幅值峭度计算的方法用伪代码描述如下：

for i＝1to L

for j＝i to i

(1)计算尺度子频带的起始BW_bgn和终止频率BW_end：

起始频率：BW_bgn＝[(j-1)/i]×f_N

终止频率：BW_end＝[(j)/i]×f_N

(2)计算时频区间[0,tend,BW_bgn,BW_end]的时域幅值均值：

(3)计算幅值平方：

(4)求幅值平方的傅里叶变换：

(5)求幅值平方峭度

End

End

峭度矩阵KA＝{KA_i,j,i＝1～L,j＝1～L}的形式如下：

设L尺度每个尺度子时频空间的频带

对应的像素宽度为Δu，则峭度色谱像素矩阵Ax＝{Ax_i,j,i＝1～L,j＝1～L×Δ_L},峭度色谱像素矩阵构造方法用伪代码表示如下：

for i＝1to L

for j＝1to j

填充峭度色谱像素矩阵

end

end

峭度色谱排列形式如图2所示。

得到信号的尺度平方幅值峭度色谱后，选择色谱图上颜色最深的区域对应的尺度频带，采用下式(10)重构该区域的时域信号。

式中，f₁、f₂分别为所选尺度频带的起始频率和终止频率，t₁、t₂分别为原始信号的起始时刻和结束时刻。

具体应用实例：

测试轴承为SKF6203－2RS型深槽滚动轴承，滚动体个数为8，内圈直径为17mm，外圈直径为40mm，滚动体直径为6.75mm，滚道节径为28.5mm。设轴承的回转频率为f_r，内圈、外圈、滚动体及支架特征频率为4.9469f_r、3.0530f_r、3.9874f_r和0.3817f_r。

图3为采用加速度传感器从轴承支座采集的轴承内圈损伤的振动信号及其频谱，采样频率为12000Hz，信号长度为8192，此时转轴回转频率f_r为29.95Hz的一组振动信号。此时，内圈、外圈、滚动体及支架特征频率分别为148.16Hz、91.44Hz、119.42Hz和11.43Hz。

取轴承内圈故障特征频率的3倍来计算最大尺度，得到最大分解尺度为14，图4为幅值谱峭度图，图中用圆圈标出了颜色最深的子频带。图5为该频带的重构信号，该重构信号的包络谱见图6，图6中谱峰对应的频率为148Hz(图中标记为f_in)，该频率为轴承内圈的特征频率。

Claims (1)

1.一种基于信号时频分解的幅值谱峭度图的实现方法，其特征在于，该方法首先用频率切片小波变换对振动信号时频分解获得其时频分解矩阵；然后，在每一尺度分别抽取该尺度每个给定的子时频空间对应的子矩阵，按时间方向求取各个子矩阵的时域均值，把该时域均值作为该子时频空间信号分量幅值随时间变化的趋势，用该时域均值的幅值平方谱峭度按照尺度—频带排列，得到信号的谱峭度尺度—频带平面图，即时频幅值谱峭度图，用颜色深度表示谱峭度的大小，峭度值越大对应谱平面区域的颜色越深，颜色最深的区域对应的频带为特征频带；其中：

所述的频率切片小波变换是一种时频分解方法，振动信号变换之前，需要选择频率切片函数，变换结果为二维时频矩阵；

所述时域均值是按照矩阵的行或列求所有行元素或列元素的平均值；

具体实现方法如下：

(1)根据故障特征频率和振动信号的采样频率确定峭度谱的尺度以及子时频空间的个数，即最大尺度对应的子时频空间个数；

(2)去除采集的振动信号中的直流分量；

(3)选择频率切片函数，对去除直流分量的信号做频率切片小波变换得到其时频分解矩阵；

(4)由时频分解矩阵求出信号的时频幅值矩阵；

(5)在每一个尺度，求出每个子时频空间对应的时频幅值子矩阵的时域幅值平均向量，再求该向量平方的频谱，然后求出该频谱的幅值峭度，每一尺度的所有子时频空间的幅值峭度构成谱峭度向量；

(6)由每个尺度的谱峭度向量构造谱峭度矩阵，再将谱峭度矩阵映射到尺度—频带平面图，以色度深浅来表示谱峭度的大小，得到信号的谱峭度图；

(7)信号分析时，在谱峭度图上颜色最深的区域为特征频带，用频率切片小波逆变换分离出该特征频带的信号分量，以备进一步分析。

Images