本发明属于机械设备故障诊断技术领域，具体涉及一种基于信号时频分解的频域幅值谱峭度图的实现方法。
背景技术：
采用振动信号在线监测是机械设备状态监测和故障诊断常用的方法之一，由于工作过程中，机械设备中的机械传动、环境噪声、电磁干扰等因素，从现场设备监测点采集的振动信号含有较强的噪声，有用的信息通常被淹没，因此，从含有噪声的振动信号中分离出含有故障特征的特征分量对正确识别故障具有重要的作用。
在机械设备中，齿轮和轴承是传动系中最常用的零部件，在设备工作过程中，当它们出现损伤时，振动信号会出现与故障的特征频率相关的信号分量。在某一时刻，当转速已知时，不同的故障特征频率可以通过理论计算确定下来，因此，通过检查振动信号中是否含有特征频率分量来确定故障类型及损伤位置。
齿轮和轴承出现故障时，故障的特征频率与传动系的结构及其回转频率有关，故障的发生使振动信号出现幅值调制现象，幅值调制的频率与故障特征频率相关，因此，共振解调方法是分离出这些分量的常用方法。共振解调关键是确定包含幅值调制分量的特征频带及其中心频率[1]。文献[2]提出了一种谱峭度确定特征频带的方法，这种方法对给定带宽的振动信号进行短时傅立叶变换以获取该频带的信号分量，然后对这个信号提取该分量的包络，并计算包络信号的峭度，最后通过最大峭度确定最佳频带。在此基础上，提出了快速谱峭度表示计算方法[3]，采用一对互补的滤波器组，对原始信号进行二进分解，在每一分解尺度得到一系列分量，再求这些分量的峭度，然后把每个尺度的各个峭度按频带映射到尺度—频带平面图上；为了改善上述方法的尺度分辨率，文献[3]还提出一种滤波器组的1/3二进分解树的改进方法，增加了分解次数已获得更多的分解尺度，改善了尺度—频带平面图分辨率。文献[4]采用小波包分解方法把信号分解到不同的频带中，利用各个频带的小波包分量求取峭度，然后把小波包分量的峭度按小波包频带映射，得到尺度—频带平面图。
上述方法在应用过程中均存在以下问题：
(1)采用滤波器组方法，需要确定互补的高通和低通滤波器参数及滤波器的长度，需要一定的先验知识。
(2)采用小波包分解方法，由于在分解过程中的下抽样操作，大于3尺度分解时，在其后的每一尺度小波包对应的频带发生倒置，在形成尺度—频带平面图时，需要对小波包分量对应的频带进行重排。另外，小波包采用二进塔式分解方法，尺度分辨率较低。
(3)不管是上述哪一种方法，因为采用互补滤器分解，信号分解时，高通滤波器和低通滤波器的特性限制了每一分解尺度的频带分割模式，因此，往往不能根据需求确定所需要精确频带以获取有效的分量。
以下是申请人检索的相关参考文献：
[1]、r.b.randall，j.antoni.rollingelementbearingdiagnostics—atutorial[j].mechanicalsystemsandsignalprocessing，2011，25：485–520。
[2]、j.antoni，r.b.randall.thespectralkurtosis:applicationtothevibratorysurveillanceanddiagnosticsofrotatingmachines[j].mechanicalsystemsandsignalprocessing，2006，20：308-331。
[3]、j.antoni.fastcomputationofthekurtogramforthedetectionoftransientfaults[j].mechanicalsystemsandsignalprocessing，2007，21：108–124。
[4]、yaguolei，jinglin，zhengjiahe，yanyangzi.applicationofanimprovedkurtogrammethodforfaultdiagnosisofrollingelementbearings[j]..mechanicalsystemsandsignalprocessing，2011，25：1738–1749。
技术实现要素：
针对目前现有技术中的谱峭度方法存在的缺陷或不足，本发明的目的在于，提供一种基于信号时频分解的频域幅值谱峭度图的实现方法。
为了实现上述任务，本发明采用如下的技术解决方案：
一种基于信号时频分解的频域幅值谱峭度图的实现方法，其特征在于，该方法首先对振动信号进行频率切片小波变换得到时频分解矩阵，以此为基础，按照每一尺度给定的子频带进行信号重构，再求取重构信号的平方幅值，经傅里叶变换得到重构信号的平方谱，对平方谱的幅值序列求频谱峭度，最后将频谱峭度值按照尺度—频带排列，得到频谱峭度图，在该频谱峭度图中用颜色深度表示频谱峭度的大小，频谱峭度值越大对应谱平面区域的颜色越深；在应用时，选取颜色最深的平面区域所对应的频带进行重构，得到时域特征信号用于进一步分析；其中：
所述的频率切片小波变换是一种时频分解方法，信号变换之前，需要选择频率切片函数，变换结果为二维时频矩阵；
所述的信号重构是指对所选取的频带进行频率切片小波逆变换分离出该频带的信号分量的过程；
所述的尺度是指对信号进行频率切片小波变换时对信号奈奎斯特频带分割的子频带数。
具体的实现方法如下：
(1)根据故障特征频率和振动信号的采样频率确定峭度谱的尺度以及最大子频带个数，即最大尺度对应的频带个数；
(2)去除采集的振动信号中的直流分量；
(3)选择频率切片函数；
(4)在每一尺度，针对每个子频带，对去除直流分量的振动信号做频率切片小波变换分离出每个子频带的信号分量，然后，求出每个子频带信号分量的幅值平方的频谱峭度，得到频谱峭度矩阵；
(5)由频谱峭度矩阵构造频域幅值谱峭度尺度—频域平面分布图，即频谱峭度图；
(6)在频谱峭度图寻找颜色最深的区域，确定该区域所对应的频率范围，即特征频带；
(7)针对上述特征频带，用频率切片小波变换分离出该频带的信号分量，即特征信号分量，以备进一步分析。
本发明的基于信号时频分解的频域幅值谱峭度图的实现方法，带来有益效果在于，改善了尺度的分辨率，可以用于机械设备、电气系统的信号分析和故障特征提取，特别是对具有调频调幅特征的信号，该方法可以图示出调频调幅分量的最佳解调频带。
附图说明
图1是频带分割示意图；
图2是幅值峭度计算方法流程图；
图3是频谱峭度色像素矩阵构造方法流程图；
图4是频谱峭度图排列形式图；
图5是振动信号的时域波形图；
图6是振动信号的频谱；
图7是尺度为8的频谱峭度图；
图8是谱值最大频带的重构信号分量的时域波形图；
图9是重构信号分量的包络谱；
图10是尺度为14的频谱峭度图。
以下结合附图和实施例对本发明作进一步的详细说明。
具体实施方式
本实施例给出一种基于信号时频分解的频域幅值谱峭度图的实现方法，采用频率切片小波逆变换可以重构任意给定子频带的信号分量的特点，根据分析需求把信号的频带分成若干个期望的子频带，通过获取子频带信号分量的幅值平方频率谱的幅值峭度，构造尺度—频率谱峭度平面色度图，即谱峭度图，为特征分量的提取提供依据。
该方法首先对振动信号进行频率切片小波变换得到时频分解矩阵，以此为基础，按照每一尺度给定的子频带进行信号重构，再求取重构信号的平方幅值，经傅里叶变换得到重构信号的平方谱，对平方谱的幅值序列求频谱峭度，最后将频谱峭度值按照尺度—频带排列，得到频谱峭度图，在该频谱峭度图中用颜色深度表示频谱峭度的大小，频谱峭度值越大对应谱平面区域的颜色越深；在应用时，选取颜色最深的平面区域所对应的频带进行重构，得到时域特征信号用于进一步分析；其中：
所述的频率切片小波变换是一种时频分解方法，信号变换之前，需要选择频率切片函数，变换结果为二维时频矩阵；
所述的信号重构是指对所选取的频带进行频率切片小波逆变换分离出该频带的信号分量的过程；
所述的尺度是指对信号进行频率切片小波变换时对信号奈奎斯特频带分割的子频带数；
本实施例中，上述频谱峭度图为信号的时频分解分量的包络平方频域幅值谱峭度色度图，它也称为信号的一种信号分量平方的频域幅值的谱峭度尺度—频域平面分布图。
具体实现方法如下：
(1)根据故障特征频率和振动信号的采样频率确定峭度谱的尺度以及最大子频带个数，即最大尺度对应的频带个数；
(2)去除采集的振动信号中的直流分量；
(3)选择频率切片函数；
(4)在每一尺度，针对每个子频带，对去除直流分量的信号做频率切片小波变换分离出每个子频带的信号分量，然后，求出每个子频带信号分量的平方幅值的频谱峭度，得到频谱峭度矩阵；
(5)由频谱峭度矩阵构造频域幅值谱峭度尺度—频域平面分布图，即频谱峭度图；
(6)在频谱峭度图寻找颜色最深的区域，确定该区域所对应的频率范围，即特征频带；
(7)针对上述特征频带，用频率切片小波变换分离出该频带的信号分量，即特征信号分量，以备进一步分析。
以下给出具体的过程。
设p(t)为频率切片函数，时域信号s(t)在频带[fbgn,fend]的频率切片小波变换时频分解为：
式中κ>0，
设κ为任意给定值时，由于ω＝2πf，把信号s(t)的频率切片小波变换时频分解w(t,ω,κ)写为wκ(t,f)，则在时频空间上通过逆变换重构时频区域[t1,t2,f1,f2]的信号分量sx(t)的公式为：
设s(t)振动信号的采样频率为fs，得到长度为n信号序列s＝{si,i＝1～n}，它对应时间区间为[0,...,tj,...,tn-1]，其中，tj＝j/fs,j＝0～n-1；
去除振动信号中的直流分量公式为：
其中
采用公式(1)对信号序列x进行频率切片小波变换，得到x在时频区间[0,tn-1,0,fs/2]的时频分解矩阵为wκ(t,f)＝{w(tk,fk),tk＝0～tn-1,fk＝0～fs/2}。
设故障特征频率为fk，定义最小分解带宽为：δbw＝(3～5)fk,则谱尺度为：
其中，fn＝fs/2，fn为奈奎斯特频率。
在每一谱尺度，对信号的奈奎斯特频带做等宽频带分割，尺度带宽为：
其中，li＝1,2,...,l。
信号s(t)的l尺度频带分割示意图如图1所示。
则li尺度的第j个时频区间的频带为：(j-1)δi～jδi，j＝1,2,...,li。
该频带信号时域分量为：
时域分量的幅值平方为：
ai,j(t)＝|fi,j(t)|2(7)
对式(7)做傅里叶变换
式中ω＝2πf。
设
幅值峭度计算的方法如图2所示，具体包括：
步骤i、设置尺度l；
步骤ii、初始化i＝1，j＝1；
步骤iii、计算尺度子频带的起始频率和终止频率：
其中，起始频率：bw_bgn＝[(j-1)/i]×fn
终止频率：bw_end＝[(j)/i]×fn
步骤iv、利用式(1)计算时频分解矩阵wk(t,f)；
利用式(2)重构出频带[bw_bgn,bw_end]的信号分量fi,j(t)；
步骤v、由式(7)、(8)和(9)计算信号分量fi,j(t)的平方幅值谱的峭度kai,j；
步骤vi、若j>i，则执行步骤vii；若j<＝i，执行步骤iii；
步骤vii，若j>l，则算法结束；若j<＝l，执行步骤iii。
幅值峭度矩阵ka＝{kai,j,i＝1～l,j＝1～l}的形式如下：
设l尺度每个尺度频带
步骤a、设置尺度l，像素宽度△u；
步骤b、初始化i＝1，j＝1；
步骤c、填充峭度色谱像素矩阵；
步骤d、若j>i，则执行步骤e，若j<＝i，则执行步骤c；
步骤e，若j>l，则算法结束，若j<＝l，则执行步骤c。
频谱峭度色像素排列如图4所示。
得到信号的尺度平方幅值频谱峭度图后，选择频谱峭度图上颜色最深的区域对应的尺度频带，采用式(2)重构该区域的时域信号，再通过包络分析，观察综合包络谱曲线峰值对应频率，确定损伤部位。
具体应用实例：
测试轴承为skf6203－2rs型深槽滚动轴承，滚动体个数为8，内圈直径为17mm，外圈直径为40mm，滚动体直径为6.75mm，滚道节径为28.5mm。设轴承的回转频率为fr，内圈、外圈、滚动体及支架特征频率为4.9469fr、3.0530fr、3.9874fr和0.3817fr。
图5为采用加速度传感器从轴承支座采集的轴承内圈存在损伤的振动信号，图6为振动信号的频谱采样频率为12000hz，信号长度为8192，此时转轴回转频率fr为29.95hz的一组振动信号。此时，内圈、外圈、滚动体及支架特征频率分别为148.16hz、91.44hz、119.42hz和11.43hz。
取轴承内圈故障特征频率的5倍来计算最大尺度，得到最大分解尺度为8，图7为频谱峭度图，图中第5尺度第5子频带的颜色最深(图中圆圈所在)，它对应的频带为4800hz～6000hz。图8为该频带的重构信号，该重构信号的包络谱见图9，图9谱峰对应的频率为148hz，该频率为轴承内圈的特征频率。
取轴承内圈故障特征频率的3倍来计算最大尺度，得到最大分解尺度为14，图10为频谱峭度图，图中依然是第5尺度第5子频带的颜色最深(图中圆圈所在)，它对应的频带为4800hz～6000hz。