CN112101245A - 基于频域窗函数的短时傅里叶变换机械冲击特征提取方法 - Google Patents

Abstract

本发明属于机械零件结构健康监测技术领域，公开了一种基于频域窗函数短时傅里叶变换的机械冲击特征提取方法，利用最大相关峭度反卷积的方法对振动信号进行滤波，使信号的质量得到提高；在此基础上，提出基于频域窗函数的短时傅里叶变换的相关理论；通过频域窗函数，实现了二维时频平面中时间的精准定位和冲击特征的准确识别。本发明的可行性通过对多组分仿真模拟信号进行分析得到了验证，并将其应用于旋转机械综合实验台的实测滚动轴承故障数据分析，结果证明该方法在轴承冲击特征的定位和识别中有着较好效果。本发明能够有效解决噪声干扰的问题，提高时频平面的分辨率，进一步凸显故障特征，实现复杂环境下故障的准确识别。

Description

基于频域窗函数的短时傅里叶变换机械冲击特征提取方法

技术领域

本发明属于机械零件结构健康监测技术领域，尤其涉及一种基于频域窗函数的短时傅里叶变换机械冲击特征提取方法。

背景技术

目前，轴承作为旋转机械的重要零部件，起着支撑机械旋转体和降低运动摩擦的作用，同时也是各类旋转机械中应用最广泛的零部件之一。轴承在工作中如果发生局部故障，会对整个设备的性能和精度产生重要的影响，其相应的故障特征如冲击也会表现在振动信号中。因此，如果能成功提取该冲击特征，即可有效地对滚动轴承进行故障诊断。在数字信号处理领域，短时傅里叶变换(Short-Time Fourier Transform，STFT)是目前常用的信号处理工具之一，在时频分析领域具有重要作用。STFT是对信号加上一个沿着时间轴移动的短时窗函数，由该短时窗对信号在各个时刻附近截取非平稳信号，此时可将短时窗内的信号看成是平稳信号，并分别对截取结果进行傅里叶变换(Fourier Transform，FT)，得到各时刻附近的频谱，即时频谱。经STFT处理后的信号具有时域和频域的局部化特性，可以使用其来分析信号的时频特性。STFT通过将非平稳信号分割成许多包含准平稳部分的帧来增加时间维度，并使用窗函数来减少频谱中的旁瓣。但是，STFT是由一系列的平行于时间轴的窗来截取信号序列，然后分别进行傅里叶变换的，显然其不适用于强时变和强调频信号。小波变换(Wavelet Transform，WT)使STFT的一些局部化思想得到了发展，解决了窗口不随频率变化的缺陷。但目前应用小波变换的算法中仍然存在着频率混叠、小波基函数和分解层数选择困难等缺点。同步压缩小波变换(Synchrosqueezing Wavelet Transform，SWT)是基于小波的一种改进时频分析的方法，通过压缩算子对时频域进行重排，将信号时频平面内任意点的时频分布移至能量的重心，使瞬时频率的能量集中得到进一步的增强，可以很好地解决传统时频分析方法的时频模糊问题。但是传统的SWT方法对于快变的频率信号的匹配效果较差，在处理多组分非平稳信号时，不可避免地会产生模糊，对时频表达的可读性造成严重干扰。在通过传感器获得的信号中，常常混淆有许多其它频率的噪声干扰。由于这些干扰的存在，不管是STFT还是SWT，其时频分析效果都容易被影响到，有时正确的测量值无法获得，有时甚至有用的信号在这些干扰噪声中被淹没

通过上述分析，现有技术存在的问题及缺陷为：1)现有方法均存在理论上的不足，例如传统短时傅里叶变换STFT中时域窗函数的宽度往往是固定的，无法实现窗宽根据信号的特征进行自适应变化；同步压缩变换作为一种信号后处理方法，其无法有效处理强调制或者强时变的信号，并且其也是以STFT作为理论基础，因此改进经典时域加窗STFT具有重要的理论价值。2)在工程实际中，噪声等无关因素的干扰是普遍存在的现象，其严重影响了信号的质量，淹没了有用信号的特征，降低了信号处理方法的实用性。。

解决以上问题及缺陷的难度为：1)如何提出一种有效的方法，可以降低噪声的干扰，提高信号的质量；2)如何改进传统经典时域加窗STFT的理论模型，通过频域加窗的方法，提高复杂信号中代表机械故障的冲击成分在时频域中的分辨率。

解决以上问题及缺陷的意义为：1)理论上的意义：提出一种新的基于频域加窗的短时傅里叶变换方法，可以提高其对复杂多组分信号的分辨率，可以更好地刻画信号的细节特征，该方法可以进一步拓展到其他领域；2)强噪声干扰下的冲击特征识别一直是机械故障诊断的核心问题，该方法的提出可以有效地提高对冲击特征的识别精度，对于设备的健康管理和结构健康监测领域，具有突出的应用价值。

发明内容

针对现有技术存在的问题，本发明提供了一种基于频域窗函数的短时傅里叶变换机械冲击特征提取方法。

本发明是这样实现的，一种基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，所述基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，包括：

步骤一，利用最大相关峭度反卷积(MCKD)的方法对振动信号进行滤波，使信号的质量得到提高；

步骤二，提出基于频域窗函数的短时傅里叶变换(STFT-FD)的相关理论；

步骤三，通过频域窗函数，实现了二维时频平面中时间的精准定位和冲击特征的准确识别。

进一步，所述步骤一中，最大相关峭度反卷积(MCKD)的方法，包括：

y(t)是属于实数域平方可积函数空间的信号，其傅里叶变换FT定义为：

短时傅里叶变换STFT定义为：

上式中，g(t)为短时窗函数，当窗函数g(t)＝1时，STFT则简化成了传统的傅里叶变换；

在采集信号的过程中，传感器收集的振动信号里会包含很多噪声；对特征信号进行降噪处理，将相关峭度作为测量尺度，通过降低有用信号中噪声所占比例来提高轴承故障信号的峭度，从而凸显有用信号中代表故障的脉冲成分；

w(n)是输入信号y(n)的周期性冲击分量，MCKD在不考虑噪声的背景下，选择一个合适的滤波器ψ(k)，使冲击分量w(n)相关峭度达到最大，同时起到降噪的效果：

上式中，ψ＝[ψ₁,ψ₂,…,ψ_L]^T，L是滤波器的长度；

信号的相关峭度(CK)需要由MCKD来进行最大化，使信号处理的结果达到最优；相关峭度的表达式为：

上式中，T为信号周期，M是位移阶次。

本发明实施例中MCKD算法的目标函数为：

即求解方程：

式(6)的解即为滤波器的系数，求解方程后得到的滤波器系数组合可表示为：

上式中，p＝0,T,2T,…,mT，将滤波器系数代入式(3)中，获得振动信号的冲击分量w(n)。

进一步，所述步骤二中，基于频率窗函数的短时傅里叶变换的方法包括：

对于一个冲击信号w(n)，MS为信号样本数量，以下窗函数：

φ(l),l∈{1,…,MW(f)} (11)

上式中，利用时间索引l定义窗口，时间索引l从1开始，MW(f)表示窗口的样本数；采样间隔为Ts，则与每个周期的给定采样数q相对应的频率f表示为：

将窗口大小设置为频率相关，定义MC为窗口函数内的周期数或循环数；每个频率的窗口大小由以下公式确定：

基于频域窗函数的短时傅里叶变换的相关理论为：

在t附近选择大小为MW(f)的时间间隔，使用索引l∈{1,…,MW(f)}覆盖此窗口函数的时间间隔{t-MW(f)/2+1,…,t+MW(f)/2}；

下式(14)表示应用窗函数之后的加窗信号

其中，窗口信号仅在时间间隔{MW(f)/2,…,MS-MW(f)/2}中定义；

上式中，l∈{1,…,MW(f)}；由于离散傅里叶变换(DFT)的第一项是常量，加窗信号的采样数等于MW(f)，因此对于时频表示W(t,f)，仅需计算DFT的(1+MC)项；

在标准STFT中，窗口大小以秒为单位固定，可以在给定的时刻使用相同的DFT计算STFT的所有项；但是，窗口大小与频率相关，必须根据不同窗口大小DFT计算的项；

如前所述，将STFT-FD加窗信号DFT计算的(1+MC)项作为参考；

上式为解释所提出方法的关键，由于窗口的大小取决于变量q，因此用于获取结果的FFT随频率分量而有所不同；则式(16)可以进一步表示为：

上式中，

若将MC设置为偶数，则MW(f)/2始终是整数；

通过在式(12)中代入式(17)可知，研究的STFT-FD的计算公式，如式(18)所示；

上式中，

进一步，所述步骤三中，先利用MCKD对振动信号进行预处理，然后基于STFT-FD对降噪后的信号进行冲击特征提取。

本发明另一目的在于提供一种基于频域窗函数的短时傅里叶变换机械冲击特征提取系统包括：

滤波模块，利用最大相关峭度反卷积的方法对振动信号进行滤波；

短时傅里叶变换模块，对滤波后数据进行基于频域窗函数的短时傅里叶变换；

识别模块，通过频域窗函数，进行二维时频平面中时间的精准定位和冲击特征识别。

本发明另一目的在于提供一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行如下步骤：

利用最大相关峭度反卷积的方法对振动信号进行滤波；

对滤波后数据进行基于频域窗函数的短时傅里叶变换；

通过频域窗函数，进行二维时频平面中时间的精准定位和冲击特征识别。

本发明另一目的在于提供一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如下步骤：

利用最大相关峭度反卷积的方法对振动信号进行滤波；

对滤波后数据进行基于频域窗函数的短时傅里叶变换；

通过频域窗函数，进行二维时频平面中时间的精准定位和冲击特征识别。

本发明另一目的在于提供一种实施所述基于频域窗函数的短时傅里叶变换机械冲击特征提取方法机械轴承运行信息监测终端。

结合上述的所有技术方案，本发明所具备的优点及积极效果为：

(1)本发明区别于经典的短时傅里叶变换方法和小波分析、同步压缩变换等方法，提出了联合MCKD滤波与基于频域窗函数的短时傅里叶变换的机械故障特征提取方法，首先利用MCKD实现信号噪声的抑制，然后将待分析信号在频域内利用变化窗宽的窗函数进行时频分析。由于在频率中调整了窗口大小，因此能更好地适应信号的变化。本发明中基于频域窗函数的短时傅里叶变换机械冲击特征提取方法中的可行性通过对多组分仿真模拟信号进行分析得到了验证，并将其应用于旋转机械综合实验台的实测滚动轴承故障数据分析，结果证明该方法在轴承冲击特征的定位和识别中有着较好效果。

为了突出有用信号，抑制噪声干扰，需要先进行降噪处理。最大相关峭度反卷积(MCKD)具有优秀的消噪能力和特征提取能力，已在多个领域得到广泛应用。与常用的降噪技术如最小熵反褶积相比，MCKD能得到更好的实际效果。针对多组分的复杂振动信号处理，首先利用MCKD的手段对采集到的振动信号进行降噪，然后提出基于频域固定窗函数的短时傅里叶变换(STFT-FD)，在利用STFT基本概念和理论框架的同时，将窗口大小在频域中固定。不同的频率使用不同的窗口大小，高频时使用小窗口，低频时使用大窗口。利用STFT-FD，可以获得更好的频率分辨率，并且不需要多分辨率技术的带通滤波器组，也不需要评估自适应技术的局部信号特性。本发明中对多组分信号进行的仿真模拟分析，并应用于实测滚动轴承故障，使该方法在定位和识别轴承冲击特征的有效性和优越性方面得到了验证。

对比的技术效果或者实验效果包括：

将本发明提供了一种基于频域窗函数的短时傅里叶变换机械冲击特征提取方法与传统STFT及同步压缩变换进行比较。可以看出，传统STFT对时变信号的处理效果较差，时频模糊的现象比较严重，冲击信号的时频特征无法准确地识别。同步压缩变换地计算结果尽管实现了频率方向上的时频系数重排，但是没有考虑时间的特性和定位，其对数值仿真冲击信号的分析效果较差。本发明提供了一种基于频域窗函数的短时傅里叶变换机械冲击特征提取方法计算所得的时频平面，可以清楚地识别出冲击特征，能够更加清晰的刻画非平稳信号的时间特性，对于冲击特征具有更好的时频分辨率。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的基于频域窗函数的短时傅里叶变换机械冲击特征提取方法流程图。

图2是本发明实施例提供的基于MCKD-STFT-FD的轴承故障诊断流程图。

图3是本发明实施例提供的数值仿真信号时域图。

图中：a、原始冲击信号；b、含噪声的冲击信号。

图4是本发明实施例提供的小波时频分析的结果示意图。

图5是本发明实施例提供的不同时频分析方法对仿真信号处理的结果示意图。

图中：a、STFT时频表达结果；b、SST时频表达结果。

图6是本发明实施例提供的基于频域窗函数的短时傅里叶变换示意图。

图7是本发明实施例提供的旋转机械振动故障诊断试验台示意图。

图中：a、实验装置实物图；b、实验装置结构简图。

图8是本发明实施例提供的实测轴承振动信号的时域波形和频谱分析结果

图中：1、振动信号的时域图；2、振动信号的频谱图。

图9是本发明实施例提供的轴承故障信号包络谱分析结果示意图。

图10是本发明实施例提供的比较不同方法的时频分析结果示意图。

图中：(a)STFT时频表达的结果；(b)SST时频表达的结果示意图。

图11是本发明实施例提供的对轴承故障信号分析示意图。

图中：a、时间重排同步压缩变换；b、时域冲击波形。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

针对现有技术存在的问题，本发明提供了一种基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，下面结合附图对本发明作详细的描述。

如图1所示，本发明实施例提供的基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，包括：

S101：利用最大相关峭度反卷积(MCKD)的方法对振动信号进行滤波，使信号的质量得到提高；

S102：提出基于频域窗函数的短时傅里叶变换(STFT-FD)的相关理论；

S103：通过频域窗函数，实现了二维时频平面中时间的精准定位和冲击特征的准确识别。

本发明实施例提供的S101中，最大相关峭度反卷积(MCKD)的方法，包括：

y(t)是属于实数域平方可积函数空间的信号，其傅里叶变换FT定义为：

短时傅里叶变换STFT定义为：

上式中，g(t)为短时窗函数，当窗函数g(t)＝1时，STFT则简化成了传统的傅里叶变换；

在采集信号的过程中，传感器收集的振动信号里会包含很多噪声；对特征信号进行降噪处理，将相关峭度作为测量尺度，通过降低有用信号中噪声所占比例来提高轴承故障信号的峭度，从而凸显有用信号中代表故障的脉冲成分；

w(n)是输入信号y(n)的周期性冲击分量，MCKD在不考虑噪声的背景下，选择一个合适的滤波器ψ(k)，使冲击分量w(n)相关峭度达到最大，同时起到降噪的效果：

上式中，ψ＝[ψ₁,ψ₂,…,ψ_L]^T，L是滤波器的长度；

信号的相关峭度(CK)需要由MCKD来进行最大化，使信号处理的结果达到最优；相关峭度的表达式为：

上式中，T为信号周期，M是位移阶次。

本发明实施例中MCKD算法的目标函数为：

即求解方程：

式(6)的解即为滤波器的系数，求解方程后得到的滤波器系数组合可表示为：

上式中，p＝0,T,2T,…,mT，将滤波器系数代入式(3)中，获得振动信号的冲击分量w(n)。

本发明实施例提供的S102中，基于频域窗函数的短时傅里叶变换(STFT-FD)的相关理论为：

对于一个冲击信号w(n)，定义MS为信号样本数量，考虑以下窗函数：

φ(l),l∈{1,…,MW(f)} (11)

上式中，利用时间索引l定义了窗口，时间索引l从1开始，MW(f)表示窗口的样本数。如果采样间隔为Ts，则与每个周期的给定采样数q相对应的频率f可以表示为：

将窗口大小设置为频率相关，定义MC为窗口函数内的周期数或循环数。每个频率的窗口大小(即窗口的样本数)由以下公式确定：

在t附近选择大小为MW(f)的时间间隔，使用索引l∈{1,…,MW(f)}覆盖此窗口函数的时间间隔{t-MW(f)/2+1,…,t+MW(f)/2}；

下式(14)表示应用窗函数之后的加窗信号

其中，窗口信号仅在时间间隔{MW(f)/2,…,MS-MW(f)/2}中定义；

上式中，l∈{1,…,MW(f)}；由于离散傅里叶变换(DFT)的第一项是常量，加窗信号的采样数等于MW(f)，因此对于时频表示W(t,f)，仅需计算DFT的(1+MC)项；

在标准STFT中，窗口大小以秒为单位固定，可以在给定的时刻使用相同的DFT计算STFT的所有项；但是，窗口大小与频率相关，必须根据不同窗口大小DFT计算的项；

如前所述，将STFT-FD加窗信号DFT计算的(1+MC)项作为参考；

上式为解释所提出方法的关键，由于窗口的大小取决于变量q，因此用于获取结果的FFT随频率分量而有所不同；则式(16)可以进一步表示为：

上式中，

若将MC设置为偶数，则MW(f)/2始终是整数；

通过在式(12)中代入式(17)可知，研究的STFT-FD的计算公式，如式(18)所示；

上式中，

本发明实施例提供的S103中，先利用MCKD对振动信号进行预处理，然后基于STFT-FD对降噪后的信号进行冲击特征提取。

下面结合具体实施例对本发明的技术方案作进一步的描述。

1理论分析

1.1最大相关峭度反卷积滤波

y(t)是属于实数域平方可积函数空间的信号，其傅里叶变换FT定义为：

短时傅里叶变换STFT定义为：

上式中，g(t)为短时窗函数，当窗函数g(t)＝1时，STFT则简化成了传统的傅里叶变换。

在采集信号的过程中，传感器收集的振动信号里会包含很多噪声。因此，为了更好地提取冲击特征，使轴承故障被有效识别，首先需要对特征信号进行降噪处理。最大相关峭度反卷积(MCKD)是信号去噪有效的技术之一。其将相关峭度作为测量尺度，通过降低有用信号中噪声所占比例来提高轴承故障信号的峭度，从而凸显有用信号中代表故障的脉冲成分。

w(n)是输入信号y(n)的周期性冲击分量，MCKD在不考虑噪声的背景下，选择一个合适的滤波器ψ(k)，使冲击分量w(n)相关峭度达到最大，同时起到降噪的效果：

上式中，ψ＝[ψ₁,ψ₂,…,ψ_L]^T，L是滤波器的长度。

为了使周期性脉冲更加突出，信号的相关峭度(CK)需要由MCKD来进行最大化。这样信号处理的结果才将达到最优。相关峭度的表达式为：

上式中，T为信号周期，M是位移阶次。

为了突出周期成分的冲击，MCKD算法的目标函数为：

为了找到最合适的滤波器，即求解方程：

式(6)的解即为滤波器的系数，求解方程后得到的滤波器系数组合可表示为：

上式中，p＝0,T,2T,…,mT，将滤波器系数代入式(3)中，可以获得振动信号的冲击分量w(n)。

1.2基于频率窗函数的短时傅里叶变换

对于一个冲击信号w(n)，定义MS为信号样本数量，考虑以下窗函数：

φ(l),l∈{1,…,MW(f)} (11)

上式中，利用时间索引l定义了窗口，时间索引l从1开始，MW(f)表示窗口的样本数。如果采样间隔为Ts，则与每个周期的给定采样数q相对应的频率f可以表示为：

将窗口大小设置为频率相关，定义MC为窗口函数内的周期数或循环数。每个频率的窗口大小(即窗口的样本数)由以下公式确定：

窗口大小将取决于信号的频率。例如将MC设置为4，则对于f＝50Hz的频率和周期T＝20ms，窗口大小将为4个周期(80ms)；对于f＝100Hz，T＝10ms的频率，窗口大小将为4个周期(40ms)。参数MC决定了基于频域窗函数的短时傅里叶变换STFT-FD的计算效率，低MC意味着只计算每个频率分量的几个周期，高MC值意味着需要考虑多个周期来计算变换。

与传统的STFT一样，STFT-FD也可以使用多种类型的窗口。由于应用于信号的窗口并不一定为正方形窗口，因此不使用MC＝1的值。在后面的分析中，使用N点对称的汉明窗口。

在给定的瞬间t时实现计算STFT。在t附近选择大小为MW(f)的时间间隔。使用索引l∈{1,…,MW(f)}覆盖此窗口函数的时间间隔{t-MW(f)/2+1，…，t+MW(f)/2}。下式(14)表示应用窗函数之后的加窗信号

其中，窗口信号仅在时间间隔{MW(f)/2,…,MS-MW(f)/2}中定义。

上式中，l∈{1,…,MW(f)}。由于离散傅里叶变换(DFT)的第一项是常量，加窗信号的采样数等于MW(f)，因此对于时频表示W(t,f)，仅需计算DFT的(1+MC)项。

在标准STFT中，窗口大小以秒为单位固定，因此可以在给定的时刻使用相同的DFT计算STFT的所有项。但是，窗口大小与频率相关，因此必须根据不同窗口大小DFT计算的项。如前所述，将STFT-FD加窗信号DFT计算的(1+MC)项作为参考。

上式为解释所提出方法的关键。由于窗口的大小取决于变量q，因此用于获取结果的FFT随频率分量而有所不同。则式(16)可以进一步表示为：

上式中，

若将MC设置为偶数，则MW(f)/2始终是整数。

通过在式(12)中代入式(17)可知，研究的STFT-FD的计算公式，如式(18)所示。

上式中，

1.3MCKD-STFT-FD故障诊断方法

先利用MCKD对振动信号进行预处理，然后基于STFT-FD对降噪后的信号进行冲击特征提取。具体流程如图2所示。

2数值仿真分析

为了对提出方法的有效性进行说明，首先要对数值仿真信号进行分析。主要研究强噪声环境下的冲击特征的提取，冲击信号的数学表达式如下：

z＝w(t)+b (20)

上式中，固有频率为f_m，其对应取值为f_m＝3000Hz。仿真信号的采样频率为f_s＝20000Hz，阻尼系数为a，采样时刻为t，b代表信噪比为-5dB的随机噪声。图3所示为数值仿真信号的时域图。从图3(a)中可以看出原始不含噪声信号具有明显的冲击特征。但是，添加高斯白噪声后如图3(b)所示，由于噪声的幅值较大，其从客观上造成对冲击特征识别的困难。随后，通过小波分析的手段，对数值仿真信号x进行时频分析。小波分析选用的小波基函数为morse，其计算结果如图4所示。显然，小波分析的结果能够在一定程度上发现冲击特征存在的现象，但是其时频分辨率较低，冲击特征的时间定位不够准确，时频模糊的现象非常明显。

随后，经典的时频分析方法如短时傅里叶变换STFT和近年来热门的时频分析方法同步压缩变换SST，被用于处理数值仿真信号。其时频分析结果如图5所示，由于STFT可以理解为短时窗内的傅里叶变换，其对时变信号的处理效果较差，时频分析结果如图5(a)所示。从图中我们可以发现时频模糊的现象比较严重，冲击信号的时频特征无法准确地识别。随后，同步压缩变换地计算结果如图5(b)所示，尽管其实现了频率方向上的时频系数重排，但是没有考虑时间的特性和定位，其对数值仿真冲击信号的分析效果较差。

最后，对提出的基于MCKD滤波和频域加窗短时傅里叶变换，图6即是其对数值仿真信号分析的结果。通过图中所示的时频平面，可以清楚地识别出冲击特征。究其原因，基于频域窗函数的短时傅里叶变换不同于传统的同步压缩变换方法和经典的短时傅里叶变换方法，它主要考虑了频域内的窗函数，因此能够更加清晰的刻画非平稳信号的时间特性，其对于冲击特征具有更好的时频分辨率。

3实验数据分析

本实验台为镇江千鹏机械设备故障综合实验台，可快速模拟旋转机械多种状态及振动，特别是可以模拟齿轮和轴承的多个故障以及混合故障。图7所示为实验台实体照片，其中，在图7(a)中所示的是装置的实物图，所用到装置的结构简图为图7(b)所示。整个实验装置由变速电机1、传动带2、联轴器3、单级齿轮变速箱4、制动器5、加载装置6、圆盘7和轴承座8等组成。选取滚动轴承(型号为NU205)作为故障轴承，其包含13个圆柱滚子，实测滚动体直径：7.5mm，节径：39mm。设51200Hz为采样频率，12.5Hz为旋转频率，f_o＝65Hz为经理论计算得到的外圈故障频率，f_i＝97.5Hz为内圈故障频率，其关键参数如表1所示。如果存在外圈故障，则在时域图中两个连续冲击之间的时间间隔为15.4ms。

表1实验参数及轴承故障频率

轴承故障实验台测得的振动信号的时域波形和频谱分析结果如图8所示。从振动信号的时域图可以看出，振动信号中包含有冲击特征，但是其幅值被其他成分所掩盖，不易被识别。从频谱分析结果中，依然无法发现外圈故障特征频率或者内圈故障特征频率的峰值或者倍频谐波成分。

包络谱分析是常用的振动信号处理方法，其计算结果如图9所示。从中，故障的特征频率及倍频成分也一样无法确定。

随后，如图10所示结果，将振动信号通过时频分析的手段来进行处理。图10(a)为STFT计算的结果，从中可以看出冲击特征的存在性，但是其时频能量发散的现象比较严重，时频模糊的问题没有得到解决。利用SST对实测振动信号进行时频分析的结果，如图10(b)所示，其主要用于频率方向时频系数的重排，但是对时频平面中冲击特征的识别效果并不理想。

最后，提出的基于MCKD滤波和频域加窗短时傅里叶变换对轴承振动信号的时频分析结果如图11所示。从中可以识别出冲击特征对应的时频脊线，两条邻近脊线的时间间隔是15.3ms，这与振动信号的时域特征基本对应。根据两条连续冲击特征谱线的间隔，故障特征频率65.3HZ可由计算得出，同时，将其表明为外圈故障。上述实现结果充分证明了本文方法的有效性。

4、区别于经典的短时傅里叶变换方法和小波分析、同步压缩变换等方法，本文提出了联合MCKD滤波与基于频域窗函数的短时傅里叶变换的机械故障特征提取方法，首先利用MCKD实现信号噪声的抑制，然后将待分析信号在频域内利用变化窗宽的窗函数进行时频分析。由于在频率中调整了窗口大小，因此能更好地适应信号的变化。该方法的可行性通过对多组分仿真模拟信号进行分析得到了验证，并将其应用于旋转机械综合实验台的实测滚动轴承故障数据分析，结果证明该方法在轴承冲击特征的定位和识别中有着较好效果。

在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上；术语“上”、“下”、“左”、“右”、“内”、“外”、“前端”、“后端”、“头部”、“尾部”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”等仅用于描述目的，而不能理解为指示或暗示相对重要性。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。

Claims (9)

1.一种基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，其特征在于，所述基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，包括：

利用最大相关峭度反卷积的方法对振动信号进行滤波；

对滤波后数据进行基于频域窗函数的短时傅里叶变换；

通过频域窗函数，进行二维时频平面中时间的精准定位和冲击特征识别。

2.如权利要求1所述的基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，其特征在于，最大相关峭度反卷积的方法，包括：

y(t)是属于实数域平方可积函数空间的信号，其傅里叶变换FT定义为：

短时傅里叶变换STFT定义为：

上式中，g(t)为短时窗函数，当窗函数g(t)＝1时，STFT则简化成了传统的傅里叶变换；

在采集信号的过程中，传感器收集的振动信号里会包含很多噪声；对特征信号进行降噪处理，将相关峭度作为测量尺度，通过降低有用信号中噪声所占比例来提高轴承故障信号的峭度，从而凸显有用信号中代表故障的冲击成分；

w(n)是输入信号y(n)中包含的周期性冲击分量，最大相关峭度反卷积(MCKD)在不考虑噪声的背景下，选择一个合适的滤波器ψ(k)，使冲击分量w(n)相关峭度达到最大，同时起到降噪的效果：

上式中，ψ＝[ψ₁,ψ₂,…,ψ_L]^T，L是滤波器的长度；

信号的相关峭度(CK)需要由MCKD来进行最大化，使信号处理的结果达到最优，相关峭度的表达式为：

上式中，T为信号周期，M是位移阶次。

本发明实施例中MCKD算法的目标函数为：

即求解方程：

式(6)的解即为滤波器的系数，求解方程后得到的滤波器系数组合可表示为：

上式中，p＝0,T,2T,…,mT，将滤波器系数代入式(3)中，获得振动信号的冲击分量w(n)。

3.如权利要求2所述的基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，其特征在于，基于频率窗函数的短时傅里叶变换的方法包括：

对于一个冲击信号w(n)，MS为信号样本数量，以下窗函数：

φ(l),l∈{1,...,MW(f)} (11)

上式中，利用时间索引l定义窗口，时间索引l从1开始，MW(f)表示窗口的样本数；采样间隔为Ts，则与每个周期的给定采样数q相对应的频率f表示为：

将窗口大小设置为频率相关，定义MC为窗口函数内的周期数或循环数；每个频率的窗口大小由以下公式确定：

4.如权利要求2所述的基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，其特征在于，基于频域窗函数的短时傅里叶变换的方法包括：

在t附近选择大小为MW(f)的时间间隔，使用索引l∈{1,...,MW(f)}覆盖此窗口函数的时间间隔{t-MW(f)/2+1,...,t+MW(f)/2}；

下式(14)表示应用窗函数之后的加窗信号

其中，窗口信号仅在时间间隔{MW(f)/2,...,MS-MW(f)/2}中定义；

上式中，l∈{1,...,MW(f)}；由于离散傅里叶变换(DFT)的第一项是常量，加窗信号的采样数等于MW(f)，因此对于时频表示W(t,f)，仅需计算DFT的(1+MC)项；

在标准STFT中，窗口大小以秒为单位固定，可以在给定的时刻使用相同的DFT计算STFT的所有项；但是，窗口大小与频率相关，必须根据不同窗口大小DFT计算的项；

将STFT-FD加窗信号DFT计算的(1+MC)项作为参考；

上式为解释所提出方法的关键，由于窗口的大小取决于变量q，因此用于获取结果的FFT随频率分量而有所不同；则式(16)可以进一步表示为：

上式中，

若将MC设置为偶数，则MW(f)/2始终是整数；

通过在式(12)中代入式(17)可知，研究的STFT-FD的计算公式，如式(18)所示；

上式中，

5.如权利要求1所述的基于频域窗函数的短时傅里叶变换机械冲击特征提取方法，其特征在于，通过频域窗函数，进行二维时频平面中时间的精准定位和冲击特征识别具体包括：

利用MCKD对振动信号进行预处理；

然后基于STFT-FD对降噪后的信号进行冲击特征提取。

6.一种实施权利要求1～4任意一项所述提取方法的基于频域窗函数的短时傅里叶变换机械冲击特征提取系统，其特征在于，所述基于频域窗函数的短时傅里叶变换机械冲击特征提取系统包括：

滤波模块，利用最大相关峭度反卷积的方法对振动信号进行滤波；

短时傅里叶变换模块，对滤波后数据进行基于频域窗函数的短时傅里叶变换；

识别模块，通过频域窗函数，进行二维时频平面中时间的精准定位和冲击特征识别。

7.一种计算机设备，其特征在于，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行如下步骤：

利用最大相关峭度反卷积的方法对振动信号进行滤波；

对滤波后数据进行基于频域窗函数的短时傅里叶变换；

通过频域窗函数，进行二维时频平面中时间的精准定位和冲击特征识别。

8.一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如下步骤：

利用最大相关峭度反卷积的方法对振动信号进行滤波；

对滤波后数据进行基于频域窗函数的短时傅里叶变换；

通过频域窗函数，进行二维时频平面中时间的精准定位和冲击特征识别。

9.一种实施权利要求1～5任意一项所述基于频域窗函数的短时傅里叶变换机械冲击特征提取方法机械轴承运行信息监测终端。

Images