基于移频技术的短时傅里叶变换阶比分析

D0I:10.13374/.issn1001-053x.2012.10.014 第34卷第10期 北京科技大学学报 Vol.34 No.10 2012年10月 Journal of University of Science and Technology Beijing 0ct.2012 基于移频技术的短时傅里叶变换阶比分析 李修文阳建宏黎敏徐金梧四 北京科技大学机械工程学院，北京100083 ☒通信作者，E-mail:jwxu@usth.cu.cn 摘要提出了一种基于移频技术的短时傅里叶变换阶比分析算法.该算法利用傅里叶变换在频域的卷积性质，对原始信号 在时域乘以ε咏使∫的频谱能量搬迁到零频处，按一定的频率间隔改变∫就可以在零频处得到其他频率的频谱能量，以此 来提高短时傅里叶变换在时频分析中的频率分辨率.然后在时频面上进行局部阈值降噪，同时跟踪转速的变化，最终应用到 变速机械的阶比分析中.与短时傅里叶变换分析结果对比表明，本文方法可以更加准确地跟踪到实际的转速.实际降速过程 中轴承信号利用本文方法进行阶比分析，成功提取到轴承的故障特征频率 关键词故障诊断：傅里叶变换：阶比分析：轴承：移频 分类号TH165.3:TN911.6 Order analysis based on FS-STFT and its application to varying speed machinery fault diagnosis LI Xiu-ten,YANG Jian-hong,LI Min,XU Jin-tou School of Mechanical Engineering,University of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:jwxu@ustb.edu.en ABSTRACT An order analysis method was proposed which utilizes short time Fourier transforms (STFT)based on frequency shift (FS-STFT).In the algorithm,by using the convolution properties of Fourier transforms in the frequency domain,the original signal is multiplied by ein the time domain,and the spectrum energy off is shifted to the position of zero frequency.Other spectrum energies can be got at the position of zero frequency when f is changed according to a certain frequency interval.As a result,the frequency resolution of STFT is improved in the time-frequency analysis.Then local threshold de-noising was used in the time- frequency domain,and the change of speed signals was tracked.Eventually,the mothod was applied to the order analysis of varying speed machinery.A comparison of analysis results between STFT and FS-STFT shows that the method can track actual speed signals accurately.The order analysis for faulty bearings using the method in the deceleration process also demonstrates that the fault frequency can be extracted successfully. KEY WORDS fault diagnosis;Fourier transforms:order analysis:bearings:frequency shift 设备在变速情况下的振动信号包含更丰富的故 处理方法.对非平稳信号无外乎有两种处理方法： 障信息，一些在平稳运行时不易反映的故障特征可 一是将非平稳变成平稳信号，二是研究新的非平稳 能会充分地表现出来0.实际生产中也存在一些变 特征计算方法.对于第一种处理方法，目前常用的 速设备，这些设备在正常工作情况下，转速在不同时 是阶比分析.阶比指的是基准频率（转轴转速）的倍 刻是变化的，如卷取机和转炉耳轴.如果对于这类 数，它通过等角度采样将时域的非平稳信号转换成 时变信号直接进行频谱分析，会出现谱峰“模糊 角度域的平稳信号.在故障诊断中，旋转机械的故 化”回.因此，对这类非平稳信号必须研究新的信号 障频率大多与转频相关，利用阶比分析可以较清晰 收稿日期：201108-30 基金项目：国家自然科学基金资助项目(50905013：51004013)：高等学校博士学科点专项科研基金资助项目(20090006120007)：中央高校基 本科研业务费专项(FRF-TPO9O14A:FRF-AS09008B)

第 34 卷 第 10 期 2012 年 10 月 北京科技大学学报 Journal of University of Science and Technology Beijing Vol． 34 No． 10 Oct． 2012 基于移频技术的短时傅里叶变换阶比分析 李修文 阳建宏 黎 敏 徐金梧! 北京科技大学机械工程学院，北京 100083 !通信作者，E-mail: jwxu@ ustb． edu． cn 摘 要 提出了一种基于移频技术的短时傅里叶变换阶比分析算法． 该算法利用傅里叶变换在频域的卷积性质，对原始信号 在时域乘以 e － j2πfit 使 fi 的频谱能量搬迁到零频处，按一定的频率间隔改变 fi 就可以在零频处得到其他频率的频谱能量，以此 来提高短时傅里叶变换在时频分析中的频率分辨率． 然后在时频面上进行局部阈值降噪，同时跟踪转速的变化，最终应用到 变速机械的阶比分析中． 与短时傅里叶变换分析结果对比表明，本文方法可以更加准确地跟踪到实际的转速． 实际降速过程 中轴承信号利用本文方法进行阶比分析，成功提取到轴承的故障特征频率． 关键词 故障诊断; 傅里叶变换; 阶比分析; 轴承; 移频 分类号 TH165 + . 3; TN911. 6 Order analysis based on FS-STFT and its application to varying speed machinery fault diagnosis LI Xiu-wen，YANG Jian-hong，LI Min，XU Jin-wu! School of Mechanical Engineering，University of Science and Technology Beijing，Beijing 100083，China !Corresponding author，E-mail: jwxu@ ustb． edu． cn ABSTRACT An order analysis method was proposed which utilizes short time Fourier transforms ( STFT) based on frequency shift ( FS-STFT) ． In the algorithm，by using the convolution properties of Fourier transforms in the frequency domain，the original signal is multiplied by e － j2πfi t in the time domain，and the spectrum energy of fi is shifted to the position of zero frequency． Other spectrum energies can be got at the position of zero frequency when fi is changed according to a certain frequency interval． As a result，the frequency resolution of STFT is improved in the time-frequency analysis． Then local threshold de-noising was used in the time￾frequency domain，and the change of speed signals was tracked． Eventually，the mothod was applied to the order analysis of varying speed machinery． A comparison of analysis results between STFT and FS-STFT shows that the method can track actual speed signals accurately． The order analysis for faulty bearings using the method in the deceleration process also demonstrates that the fault frequency can be extracted successfully． KEY WORDS fault diagnosis; Fourier transforms; order analysis; bearings; frequency shift 收稿日期: 2011--08--30 基金项目: 国家自然科学基金资助项目( 50905013; 51004013) ; 高等学校博士学科点专项科研基金资助项目( 20090006120007) ; 中央高校基 本科研业务费专项( FRF--TP--09--014A; FRF--AS--09--008B) 设备在变速情况下的振动信号包含更丰富的故 障信息，一些在平稳运行时不易反映的故障特征可 能会充分地表现出来［1］． 实际生产中也存在一些变 速设备，这些设备在正常工作情况下，转速在不同时 刻是变化的，如卷取机和转炉耳轴． 如果对于这类 时变信号直接进行频谱分析，会 出 现 谱 峰“模 糊 化”［2］． 因此，对这类非平稳信号必须研究新的信号 处理方法． 对非平稳信号无外乎有两种处理方法: 一是将非平稳变成平稳信号，二是研究新的非平稳 特征计算方法． 对于第一种处理方法，目前常用的 是阶比分析． 阶比指的是基准频率( 转轴转速) 的倍 数，它通过等角度采样将时域的非平稳信号转换成 角度域的平稳信号． 在故障诊断中，旋转机械的故 障频率大多与转频相关，利用阶比分析可以较清晰 DOI:10.13374/j.issn1001-053x.2012.10.014

第10期 李修文等：基于移频技术的短时傅里叶变换阶比分析 ·1191· 地表示这些故障成分.此外，等角度采样还可以去 法，通过保留时频域上的最大值进行转速拟合，对信 除与转速无关的噪声成分，有利于特征信号的提取 噪比较大的场合显然不适用，并且对多分量信号容 阶比分析的关键是获取准确的转速信号，即阶 易发生误判，无法跟踪到准确的转速信号 比跟踪.目前阶比跟踪主要有三种方式：一是利用 本文提出了一种基于移频技术的短时傅里叶变 鉴相装置直接进行等角度同步采集；二是鉴相信号 short time Fourier transforms based on frequency 和振动信号分开单独采集，最后再对振动信号进行 shif,FS-STFT),利用它可以更加准确地获得分析窗 重采样，即计算阶比跟踪(computed order tracking, 内的频率变化，并对时频分布利用局部阈值降噪后， COT)):三是在没有鉴相信号的前提下，利用时频 可以拟合成准确的转速信号，有利于最终阶比分析. 分析等其他手段计算并拟合转速信号.第一种方法 1 基于移频技术的短时傅里叶变换阶比 需要专门的设备来完成，硬件不但结构复杂、价格昂 分析 贵，而且在旋转机械转速变化快时，其跟踪精度得不 到保证.目前第二种方法即计算阶比跟踪，测量装 1.1基于移频技术的短时傅里叶变换 置得到较大简化，安装也更加方便，已成为旋转机械 短时傅里叶变换是时频分析中较常用的方法 阶比跟踪应用最广泛的一种方法，B&K、Agilent等 它的基本思路是：为了达到时域上的局部化，在信号 著名测试仪器生产商的相关产品也引进了此功能， 傅里叶变换前乘以一个时间有限的窗函数，通过窗 但它仍需要鉴相装置进行转速跟踪，在不方便安装 在时间轴上移动可以得到信号的一组局部频谱，把 鉴相装置的场合也是无能为力回.第三种方法是近 这些局部频谱组合就得到时间一频率图，如图1 几年才发展起来的-)，实用性强，但对方法的要求 所示 高，本文研究所用的就是这一类方法.常见的时频 分析方法有短时傅里叶变换(short time Fourier x(g(I-T) di)g(i-t) xtr)g(t-T) transforms,STFT),Gabor变换、Wigner--Ville分布和 小波变换(wavelet transformation,WT)等.其中 Wigner-Ville分布会产生交叉项：而小波变换得到的 是时间一尺度平面，需要通过转换，且小波基选择没 有统一原则：Gabor变换本质仍然属于短时傅里叶 变换，由于它计算简单、高效、直观和结果容易解释， 并且没有交叉项的问题已被国内外许多学者作为阶 比跟踪的手段.但是，传统短时傅里叶变换的频率 图1短时傅里叶变换的示意图 分辨率与时间窗宽选择有关，分辨率往往不高.尽 Fig.1 Schematic diagram of short time Fourier transforms (STFT) 管后来有学者图在短时傅里叶变换的基础上提出 了自适应短时傅里叶变换(Adaptive -STFT, 根据傅里叶变换的频移性质： ASTFT),但它也只是在不同的时刻选取不同的窗函 x(t)e-Rx(f+f), (1) 数，仍然不能满足转速拟合的需求。目前也有学 即对时域信号乘以e卫e,相当于对频谱进行了f。 者O提出了基于Viterbi算法的Gabor阶比分析，可 的平移，如图2所示 以达到较好的效果，但其中代价函数在实际中是未 基于移频技术的短时傅里叶变换基本思路是： 知的.此外，现有阶比跟踪多数直接采用峰值跟踪 在传统短时傅里叶变换的基础上，根据需要的分辨 a 4() 0 图2傅里叶变换的频移性质.(a)x(t)的频谱：(b)x(t)e22的频谱 Fig.2 Frequency shift properties of Fourier transforms:(a)frequency spectrum of():(b)frequeney spectrum of ()e

第 10 期 李修文等: 基于移频技术的短时傅里叶变换阶比分析 地表示这些故障成分． 此外，等角度采样还可以去 除与转速无关的噪声成分，有利于特征信号的提取． 阶比分析的关键是获取准确的转速信号，即阶 比跟踪． 目前阶比跟踪主要有三种方式: 一是利用 鉴相装置直接进行等角度同步采集; 二是鉴相信号 和振动信号分开单独采集，最后再对振动信号进行 重采样，即计算阶比跟踪( computed order tracking， COT) ［3］; 三是在没有鉴相信号的前提下，利用时频 分析等其他手段计算并拟合转速信号． 第一种方法 需要专门的设备来完成，硬件不但结构复杂、价格昂 贵，而且在旋转机械转速变化快时，其跟踪精度得不 到保证． 目前第二种方法即计算阶比跟踪，测量装 置得到较大简化，安装也更加方便，已成为旋转机械 阶比跟踪应用最广泛的一种方法，B＆K、Agilent 等 著名测试仪器生产商的相关产品也引进了此功能， 但它仍需要鉴相装置进行转速跟踪，在不方便安装 图 2 傅里叶变换的频移性质． ( a) x( t) 的频谱; ( b) x( t) e － j2πfct 的频谱 Fig． 2 Frequency shift properties of Fourier transforms: ( a) frequency spectrum of x( t) ; ( b) frequency spectrum of x( t) e － j2πfct 鉴相装置的场合也是无能为力［4］． 第三种方法是近 几年才发展起来的［5--7］，实用性强，但对方法的要求 高，本文研究所用的就是这一类方法． 常见的时频 分析 方 法 有 短 时 傅 里 叶 变 换 ( short time Fourier transforms，STFT) 、Gabor 变换、Wigner--Ville 分布和 小 波 变 换 ( wavelet transformation，WT) 等． 其 中 Wigner-Ville 分布会产生交叉项; 而小波变换得到的 是时间--尺度平面，需要通过转换，且小波基选择没 有统一原则; Gabor 变换本质仍然属于短时傅里叶 变换，由于它计算简单、高效、直观和结果容易解释， 并且没有交叉项的问题已被国内外许多学者作为阶 比跟踪的手段． 但是，传统短时傅里叶变换的频率 分辨率与时间窗宽选择有关，分辨率往往不高． 尽 管后来有学者［8］在短时傅里叶变换的基础上提出 了 自 适 应 短 时 傅 里 叶 变 换 ( Adaptive-STFT， ASTFT) ，但它也只是在不同的时刻选取不同的窗函 数，仍然不能满足转速拟合的需求． 目前也有学 者［4］提出了基于 Viterbi 算法的 Gabor 阶比分析，可 以达到较好的效果，但其中代价函数在实际中是未 知的． 此外，现有阶比跟踪多数直接采用峰值跟踪 法，通过保留时频域上的最大值进行转速拟合，对信 噪比较大的场合显然不适用，并且对多分量信号容 易发生误判，无法跟踪到准确的转速信号． 本文提出了一种基于移频技术的短时傅里叶变 换( short time Fourier transforms based on frequency shift，FS-STFT) ，利用它可以更加准确地获得分析窗 内的频率变化，并对时频分布利用局部阈值降噪后， 可以拟合成准确的转速信号，有利于最终阶比分析． 1 基于移频技术的短时傅里叶变换阶比 分析 1. 1 基于移频技术的短时傅里叶变换 短时傅里叶变换是时频分析中较常用的方法． 它的基本思路是: 为了达到时域上的局部化，在信号 傅里叶变换前乘以一个时间有限的窗函数，通过窗 在时间轴上移动可以得到信号的一组局部频谱，把 这些局部频谱组合就得到时间--频 率 图，如 图 1 所示． 图 1 短时傅里叶变换的示意图 Fig． 1 Schematic diagram of short time Fourier transforms ( STFT) 根据傅里叶变换的频移性质: x( t) e － j2πfct X( f + fc ) ， ( 1) 即对时域信号乘以 e － j2πfct ，相当于对频谱进行了 fc 的平移，如图 2 所示． 基于移频技术的短时傅里叶变换基本思路是: 在传统短时傅里叶变换的基础上，根据需要的分辨 ·1191·

·1192· 北京科技大学学报 第34卷 率△，对每一次窗内信号分析时都乘以e卫(i= 原始信号 J: 0,1,…,2公为采样频率)进行移频处理，这样可 时频分布结果 FS-STFT(.) 以分别使频率为i△f的谱线移到零频处，记录每一 1 次零频处幅值大小，得到以△f为分辨率的各频率成 分的能量，直到分析了所有的信号，最终可得到以 x时刻频谱FS-STFT(亿，f)及 新旧闽值允许差值Pmt() △f为频率分辨率的短时傅里叶变换 1.2局部阈值降噪 计算S-STFT(r,f)最大值M、最小值N: 直接得到的时频分布首先要进行降噪处理，这 初始化网值Th=(M+N)/2 有助于转速信号的获取.目前用到的降噪方法主要 FS-STFT任，f)巾所有大于等 是利用单一的门槛值或者带通滤波进行降噪@： 于Th值的平均值记为P 但是对于时变信号，每个时刻的频谱结构是变化的， SSrT亿，f)中所有小于 并且转速不一样，同一频率成分对应的能量幅值也 Th值的平均值记为P, 是不一样的.结合基于移频技术的短时傅里叶变换 新阈值T=(P+PV2: 特点，对原始信号的二维时频分布SF-STFT(r,f), 计算新旧阀值差值d=ahs(Th-: 计算各个r时刻的阈值Th().利用各自的阈值进 更新阀值Th=T 行降噪处理，局部阈值的计算流程如图3所示 1.3基于移频技术的短时傅里叶变换阶比分析 差值d 流程 Y 转速提取的最终目的是要将时间域非平稳的信 时刻局部阀值 Th(r)=Th; 号通过等角度采样转换到角度域，从而可以运用平 t=t+l 稳信号的分析和处理方法.在基于移频技术的短时 B 傅里叶变换基础上，通过局部阈值降噪后，在时频域 D总时长E> 上进行转速的跟踪，确定每个分析窗内的转频值. YV 这里需要借助于先验知识，即已知的某个时刻转速 结宋 估计值.这在实际中还是容易满足的，如果是升速 图3局部阀值计算流程 或者降速过程，也总是可以找到一个相对平稳的转 Fig.3 Calculation flow diagram of the local threshold 速阶段.具体过程如图4所示（如果已知的是最后 时刻或者中间某时刻的转速估计值，计算方法类 时标t:,运用Langrange线性插值的方法得到t:对应 似).假设已知。时刻频率估计为f,估计误差范围 的幅值，这样最终得到的就是等角度采样的波形 为上下f,则取区间听：-f,f+f]的最大值f6作为 整个阶比分析的流程图如图5所示. 此时刻的转频跟踪值，也为下一时刻的频率估计中 t时刻的频估计为+ff们， 心，按照此方法直到搜索完所有的时间点，即得到整 其均值为下一时刻频革估计巾心 个时间的转频值.之后对这些点进行拟合，得到转 速曲线.本文采用最小二乘的多项式拟合方法，这 种拟合已被证明能够满足转速跟踪的要求00 拟合次数是拟合时的主要参数，如果值太大会增 加计算量，太小又难以保证拟合精度.若预先知道 设备是处于升速或者降速过程，n可取1或2.否则 图4转速跟踪示意图 根据拟合点的变化程度适当加大n,在保证计算速 Fig.4 Schematic diagram of speed tracking 度的同时，整个采样时间也不会持续太长，因此以小 于5为宜.如果拟合点波动很大，可采用分段拟合 转速跟踪时上下误差∫，的选取主要取决于两 的方法.最后根据需要的角域采样频率（点·rad-l) 个因素：。时刻转频估计值的精度和转速变化（波 对原始信号进行等角度重采样，搜索对应的等角度 动)程度.估计值越准确，转速波动越小，∫。可选择

北 京 科 技 大 学 学 报 第 34 卷 率 Δf，对每一次窗内信号分析时都乘以 e － j2πiΔ ( ft i = 0，1，…，fs 2Δf ，fs 为采样频率 进行移频处理 ) ，这样可 以分别使频率为 iΔf 的谱线移到零频处，记录每一 次零频处幅值大小，得到以 Δf 为分辨率的各频率成 分的能量，直到分析了所有的信号，最终可得到以 Δf 为频率分辨率的短时傅里叶变换． 1. 2 局部阈值降噪 直接得到的时频分布首先要进行降噪处理，这 有助于转速信号的获取． 目前用到的降噪方法主要 是利用单一的门槛值或者带通滤波进行降噪［9--10］; 但是对于时变信号，每个时刻的频谱结构是变化的， 并且转速不一样，同一频率成分对应的能量幅值也 是不一样的． 结合基于移频技术的短时傅里叶变换 特点，对原始信号的二维时频分布 SF-STFT( τ，f) ， 计算各个 τ 时刻的阈值 Th( τ) ． 利用各自的阈值进 行降噪处理，局部阈值的计算流程如图 3 所示． 1. 3 基于移频技术的短时傅里叶变换阶比分析 流程 转速提取的最终目的是要将时间域非平稳的信 号通过等角度采样转换到角度域，从而可以运用平 稳信号的分析和处理方法． 在基于移频技术的短时 傅里叶变换基础上，通过局部阈值降噪后，在时频域 上进行转速的跟踪，确定每个分析窗内的转频值． 这里需要借助于先验知识，即已知的某个时刻转速 估计值． 这在实际中还是容易满足的，如果是升速 或者降速过程，也总是可以找到一个相对平稳的转 速阶段． 具体过程如图 4 所示( 如果已知的是最后 时刻或者中间某时刻的转速估计值，计算方法类 似) ． 假设已知 t0 时刻频率估计为 fi，估计误差范围 为上下 fa，则取区间［fi － fa，fi + fa］的最大值 f0 作为 此时刻的转频跟踪值，也为下一时刻的频率估计中 心，按照此方法直到搜索完所有的时间点，即得到整 个时间的转频值． 之后对这些点进行拟合，得到转 速曲线． 本文采用最小二乘的多项式拟合方法，这 种拟合已被证明能够满足转速跟踪的要求［10--11］． 拟合次数 n 是拟合时的主要参数，如果值太大会增 加计算量，太小又难以保证拟合精度． 若预先知道 设备是处于升速或者降速过程，n 可取 1 或 2． 否则 根据拟合点的变化程度适当加大 n，在保证计算速 度的同时，整个采样时间也不会持续太长，因此以小 于 5 为宜． 如果拟合点波动很大，可采用分段拟合 的方法． 最后根据需要的角域采样频率( 点·rad － 1 ) 对原始信号进行等角度重采样，搜索对应的等角度 图 3 局部阈值计算流程 Fig． 3 Calculation flow diagram of the local threshold 时标 ti，运用 Langrange 线性插值的方法得到 ti 对应 的幅值，这样最终得到的就是等角度采样的波形． 整个阶比分析的流程图如图 5 所示． 图 4 转速跟踪示意图 Fig． 4 Schematic diagram of speed tracking 转速跟踪时上下误差 fa 的选取主要取决于两 个因素: t0 时刻转频估计值的精度和转速变化( 波 动) 程度． 估计值越准确，转速波动越小，fa 可选择 ·1192·

第10期 李修文等：基于移频技术的短时傅里叶变换阶比分析 ·1193· 小些 般∫选为起始估计值f的10%可满足 4: 要求 原始 FS-STFT 局部國值 转速限踪 转速 0.700.75 0.80 0.85 0.900.95 1.00 信号 降噪 (SD 拟合 % 等角度 0.4 阶比谱 角域平稳 原始信号 信号 插值 采样时标 图5基于移频技术的短时傅里叶变换阶比分析流程图 10 20 30405060 70 Fig.5 Flow diagram of order analysis based on FS-STFT 图6仿真信号时域波形(a)及频谱(b) 1.4基于移频技术的短时傅里叶变换的分辨率确 Fig.6 Waveform of the simulation signal (a)and its frequency 定原则 spectrum (b) 基于移频技术的短时傅里叶变换阶比分析中， 短时傅里叶变换窗宽取1024点，利用本文1.4 △f的选取决定时频分析的分辨率，从而影响转速跟 方法确定最优的△f值.图7是方差σ随M值的增 踪时的精度这个值理论上是越小越好，但是过小 大的变化趋势.可以看出当M等于5时，继续增大 的△f大大地增加计算量.因此，需要在保证一定分 /128001.1 辨率的前提下，确定一个合适的值.本文利用的方 分辨率影响不再明显，所以取△=(1024)×了= 法是：如果继续增大分辨率△f时，得到的转频值与 2.5Hz进行移频处理.最终得到的基于移频技术的 上一分辨率下得到的比较接近，则当前△f为最优 短时傅里叶变换时频分布结果如图8所示. 值.具体做法是在传统的短时傅里叶变换的分辨率 2.07 4的基瑞上依次计算4=普（V为整数1,2，… 1.6 Z)时基于移频技术的短时傅里叶变换结果，得到各 1.2 时刻转频的跟踪值f(M):,并计算与M-1时的转频 跟踪值结果f(M-1):的方差σ.如果σ较小并趋 0.4 于稳定，此时M对应的△f确定为最佳值 5 89 0= f(M),-f(M-1),]2N-1(M≥2). 次数，M 图7方差σ随次数M的变化趋势 (2) Fig.7 Trend of variance o with times M 其中，f(M):和f(M-1),分别对应取M和M-1时 得到的第i时刻转频值，N为转频值总数.当M=1 利用本文1.2的方法进行转速跟踪并拟合转速 时，即为传统的短时傅里叶变换 的变化曲线，设起始时刻的转频估计f。=50Hz,允 许波动范围f。=5Hz.图9是短时傅里叶变换和基 2数值仿真与实际应用 于移频技术的短时傅里叶变换分析的转速跟踪结 2.1数值仿真 果.可以看出本文的方法能较好跟踪到实际的转速 对于旋转机械的振动信号，转频总是会存在 曲线。用跟踪转频值与实际转频值的方差D来衡量 转速跟踪的精度 的：但如果转频的能量较大，往往预示着设备的某 种故障，如不平衡和碰磨2-围.模拟转速呈正弦变 D- f-F)2N-1 (3) 化的振动信号，根据瞬时频率表达地=w=2(b 其中，∫是跟踪转频值，F:是实际转频值，N是数据 dt 和w分别表示相位和角频率)，仿真信号s(t)=sin 长度.基于移频技术的短时傅里叶变换分析的转频 跟踪值与实际转频值的方差D=0.0918,而短时傅 (2πr。50+10sin(2mx0.5w）du）+rands(),信噪 里叶变换的方差D=5.0403. 比为1、采样频率为12800Hz.图6是仿真信号时域 最后，按本文1.3阶比分析得到其阶比幅值谱， 波形及频谱图.由于转速存在波动，直接频谱分析 结果如图10所示，其中角域采样率为100点· 能量泄露严重，出现谱峰的“模糊化”，无法准确分 rad-,多项式拟合次数选择4.阶比幅值谱可以看 析到相应的频率成分 到清晰的1倍阶比，也就是说原始信号中的主要成

第 10 期 李修文等: 基于移频技术的短时傅里叶变换阶比分析 小些． 一 般 fa 选为起始估计值 fi 的 10% 可 满 足 要求． 图 5 基于移频技术的短时傅里叶变换阶比分析流程图 Fig． 5 Flow diagram of order analysis based on FS-STFT 1. 4 基于移频技术的短时傅里叶变换的分辨率确 定原则 基于移频技术的短时傅里叶变换阶比分析中， Δf 的选取决定时频分析的分辨率，从而影响转速跟 踪时的精度． 这个值理论上是越小越好，但是过小 的 Δf 大大地增加计算量． 因此，需要在保证一定分 辨率的前提下，确定一个合适的值． 本文利用的方 法是: 如果继续增大分辨率 Δf 时，得到的转频值与 上一分辨率下得到的比较接近，则当前 Δf 为最优 值． 具体做法是在传统的短时傅里叶变换的分辨率 Δf0 的基础上，依次计算 Δf = Δf0 M ( M 为整数 1，2，… Z) 时基于移频技术的短时傅里叶变换结果，得到各 时刻转频的跟踪值 f( M) i，并计算与 M － 1 时的转频 跟踪值结果 f( M － 1) i 的方差 σ． 如果 σ 较小并趋 于稳定，此时 M 对应的 Δf 确定为最佳值． σ = ∑ N i = 1 ［f( M) i － f( M － 1) i ］2 N － 1 ( M≥2) ． ( 2) 其中，f( M) i 和 f( M － 1) i 分别对应取 M 和 M － 1 时 得到的第 i 时刻转频值，N 为转频值总数． 当 M = 1 时，即为传统的短时傅里叶变换． 2 数值仿真与实际应用 2. 1 数值仿真 对于旋转机械的振动信号，转频总是会存在 的; 但如果转频的能量较大，往往预示着设备的某 种故障，如不平衡和碰磨［12--13］． 模拟转速呈正弦变 化的振动信号，根据瞬时频率表达dΦ dt = ω = 2πf( Φ 和 ω 分别表示相位和角频率) ，仿真信号s( t) ( = sin 2π ∫ t 0 50 + 10sin( 2π × 0. 5u) du ) + rands( t) ，信噪 比为 1、采样频率为 12 800 Hz． 图 6 是仿真信号时域 波形及频谱图． 由于转速存在波动，直接频谱分析 能量泄露严重，出现谱峰的“模糊化”，无法准确分 析到相应的频率成分． 图 6 仿真信号时域波形( a) 及频谱( b) Fig． 6 Waveform of the simulation signal ( a ) and its frequency spectrum ( b) 短时傅里叶变换窗宽取 1 024 点，利用本文 1. 4 方法确定最优的 Δf 值． 图 7 是方差 σ 随 M 值的增 大的变化趋势． 可以看出当 M 等于 5 时，继续增大 分辨率影响不再明显，所以取 Δf = ( 12 800 ) 1 024 × 1 5 = 2. 5 Hz 进行移频处理． 最终得到的基于移频技术的 短时傅里叶变换时频分布结果如图 8 所示． 图 7 方差 σ 随次数 M 的变化趋势 Fig． 7 Trend of variance σ with times M 利用本文 1. 2 的方法进行转速跟踪并拟合转速 的变化曲线，设起始时刻的转频估计 f0 = 50 Hz，允 许波动范围 fa = 5 Hz． 图 9 是短时傅里叶变换和基 于移频技术的短时傅里叶变换分析的转速跟踪结 果． 可以看出本文的方法能较好跟踪到实际的转速 曲线． 用跟踪转频值与实际转频值的方差 D 来衡量 转速跟踪的精度． D = ∑ N i = 1 ( fi － Fi ) 2 N － 1 ． ( 3) 其中，fi 是跟踪转频值，Fi 是实际转频值，N 是数据 长度． 基于移频技术的短时傅里叶变换分析的转频 跟踪值与实际转频值的方差 D = 0. 091 8，而短时傅 里叶变换的方差 D = 5. 040 3． 最后，按本文 1. 3 阶比分析得到其阶比幅值谱， 结果 如 图 10 所 示，其中角域采样率为 100 点· rad － 1 ，多项式拟合次数选择 4． 阶比幅值谱可以看 到清晰的 1 倍阶比，也就是说原始信号中的主要成 ·1193·

·1194· 北京科技大学学报 第34卷 100p (a) 80 20 0.2 0.40.60.81.012141.6 0 020.40.60.81.012141.6 图8短时傅里叶变换()和基于移频技术的短时傅里叶变换(b)时频分析结果 Fig.8 Time-frequency analysis results of STFT (a)and FS-STFT (b) 70r 分析无法分析到外圈故障，也看不到其他相关的频 ·STFT 60 一实际曲线 率成分 FS-STFT 是50 电机 40 30%0.2040.6081i012141.6 t/9 图9短时傅里叶变换和基于移频技术的短时傅里叶变换分析 故障轴承 的转速跟踪结果对比 Fig.9 Comparison of speed tracking between STFT and FS-STFT 图11轴承试验台 分是随转速变化的成分 Fig.11 Bearing test-bed 4r 表1故障轴承参数 MWAMMWW Table 1 Parameters of the fault bearing 里-2 0 25 轴承 滚动体 外圈故障 f对应的 20 30 (角度2 rad 型号 个数，Z 频率，fu 阶比 10 6207E 9 f=0.4乐.(f表示转频) 3.6 g0.5 3 5 6 阶比 人 AlldiWalhwhwshe 图10基于移频技术的短时傅里叶变换阶比分析结果.(a)时 hm m mn mmmr m 域图：(b)阶比谱 0.5 1.01.5 2.0 2.53.0 Fig.10 Order analysis result based on FS-STFT:(a)waveform of 0.3 the time domain:(b)order spectrum 0.2 0.1 2.2实际应用 0 0 50 0 150 试验装置由一台550W(220V,50Hz)交流电 机带动，通过联轴节带动轴系运转.在轴系上装有 图12时域波形(a)及包络谱(b) 两个滚动轴承，如图11所示，其中远离电机的轴承 Fig.12 Waveform of the time domain (a)and its envelope spectrum 为可更换的6207E型故障轴承，相关参数见表1.利 () 用加速度传感器拾取降速（起始转速约为1200r· 利用本文的方法进行阶比分析，同样先确定最 min-)过程中带有外圈故障（电火花加工）的轴承 优的△f值.图13(a)是方差o随M值的增大变化 振动信号，采样频率为5120Hz.图12是原始波形 趋势.当M等于4时，继续增大分辨率影响不再明 和包络谱.由于是减速过程，转速单调减小，信号表 现为明显的非平稳特点.可以看出直接进行包络谱 显，因此取4=（器）×=1.25进行移频处

北 京 科 技 大 学 学 报 第 34 卷 图 8 短时傅里叶变换( a) 和基于移频技术的短时傅里叶变换( b) 时频分析结果 Fig． 8 Time-frequency analysis results of STFT ( a) and FS-STFT ( b) 图 9 短时傅里叶变换和基于移频技术的短时傅里叶变换分析 的转速跟踪结果对比 Fig． 9 Comparison of speed tracking between STFT and FS-STFT 分是随转速变化的成分． 图 10 基于移频技术的短时傅里叶变换阶比分析结果． ( a) 时 域图; ( b) 阶比谱 Fig． 10 Order analysis result based on FS-STFT: ( a) waveform of the time domain; ( b) order spectrum 2. 2 实际应用 试验装置由一台 550 W ( 220 V，50 Hz) 交流电 机带动，通过联轴节带动轴系运转． 在轴系上装有 两个滚动轴承，如图 11 所示，其中远离电机的轴承 为可更换的6207E 型故障轴承，相关参数见表1． 利 用加速度传感器拾取降速( 起始转速约为 1 200 r· min － 1 ) 过程中带有外圈故障( 电火花加工) 的轴承 振动信号，采样频率为 5 120 Hz． 图 12 是原始波形 和包络谱． 由于是减速过程，转速单调减小，信号表 现为明显的非平稳特点． 可以看出直接进行包络谱 分析无法分析到外圈故障，也看不到其他相关的频 率成分． 图 11 轴承试验台 Fig． 11 Bearing test-bed 表 1 故障轴承参数 Table 1 Parameters of the fault bearing 轴承 型号 滚动体 个数，Z 外圈故障 频率，fout fout对应的 阶比 6207E 9 fout = 0. 4Zfr ( fr 表示转频) 3. 6 图 12 时域波形( a) 及包络谱( b) Fig． 12 Waveform of the time domain ( a) and its envelope spectrum ( b) 利用本文的方法进行阶比分析，同样先确定最 优的 Δf 值． 图 13( a) 是方差 σ 随 M 值的增大变化 趋势． 当 M 等于 4 时，继续增大分辨率影响不再明 显，因此取 Δf = ( 5 120 ) 1 024 × 1 4 = 1. 25 Hz 进行移频处 ·1194·

第10期 李修文等：基于移频技术的短时傅里叶变换阶比分析 ·1195· 理.根据起始转速，取f6=20Hz,允许波动范围∫。= 图13(b)所示.由于是减速过程（近似），转速曲线 2Hz,对时频面进行转速跟踪，最终得到的转速值如 拟合次数选择2，拟合结果如图13(b)所示. 20 ◆一转浊跟踪值 (a) 6 一拟合结果（转速曲线） 18 16 4 三4 2 2 4.68.0.1214.16820 0.5 1.01.52.0 2.5 3.0 次数，M 图13转速跟踪结果.(a)方差σ随次数M的变化趋势：(b)转速跟踪值与拟合结果 Fig.13 Results from speed tracking:(a)trend of variance with times M:(b)speed tracking value and its fitting result 根据转速曲线对原始信号进行等角度重采样， (李辉，郑海起，杨绍普.基于幅值和相位解调分析的齿轮箱 阶比分析的结果如图14所示.从图可以看到明显 起动过程故障诊断.振动与冲击，2008,27(2)：8) 的外圈故障频率阶比3.6以及转频成分，并且还可 2]Lopatinskaia E,Zhu J,Mathew J.Monitoring varying speed ma- chinery vibrations:I.The use of non-stationary recursive filters. 以发现外圈故障频率带有转频的边带 Mech Syst Signal Process,1995,9(6):635 0.4 B] Bossley K M,McKendrick R J,Harris C J,et al.Hybrid compu- 0.3 ted order tracking.Mech Syst Signal Process,1999,13(4):627 0.2 4]Zhao X P,Hou R T.Gabor order tracking based on Viterbi algo- 22.6 7) rithm.J Mech Eng,2009.45(11)247 (赵晓平，侯荣涛.基于Viterbi算法的Gabor阶比跟踪技术. 2 3 45 6 阶比 机械工程学报，2009,45(11)：247) 5] Wang K S,Heyns P S.Application of computed order tracking, 图14基于移频技术的短时傅里叶变换阶比分析结果 Vold-Kalman filtering and EMD in rotating machine vibration. Fig.14 Result from order analysis based on FS-STFT Mech Syst Signal Process,2011,25(1):416 [6]Gao Y,Guo Y,Chi Y L,et al.Order tracking based on robust 3结论 peak search instantaneous frequency estimation.J Phys Conf Ser 2006,48:479 (1)提出了一种基于移频技术的短时傅里叶 [7]Liu W B,Guo Y,Li Z X.Removal of crossing-items of WVD 变换分析方法.在传统短时傅里叶变换的基础上， based on gabor expansion.J Vib Shock,2008,27 (10):121 利用傅里叶变换的移频性质，根据一定的频率间 (刘文彬，郭瑜，李之雄.基于Gabor展开的Wigner-Ville分布 隔依次将各频率成分进行搬移，减少能量泄露，算 的交叉项消除.振动与冲击，2008,27(10)：121) 法高效的同时，可以大大提高时频分析的频率分 8] Jones D L,Parks T W.A high resolution data-adaptive time-fre- quency representation.IEEE Trans Acoust Speech Signal Process, 辨率. 1990,38(12):2127 (2)研究了无转速信号下的阶比分析方法.首 ] Kong Q P,Song K C,Chen Y.Study of time-frequency order 先对基于移频技术的短时傅里叶变换分析的结果进 tracking of vibration signals in engine speed changing stage.I Vib 行局部阈值降噪，不同时刻选取不同的阈值，利用峰 Eng,2005,18(4):448 值搜索并跟踪转速值，之后采用最小二乘多项式拟 (孔庆鹏，宋开臣，陈鹰.发动机变速阶段振动信号时频分析 合得到转速曲线，最终用于阶比分析.对仿真和实 阶比跟踪研究.振动工程学报，2005,18(4)：448) [10]Guo Y.Research on Virtual Characteristics Analyzer System of Ro- 际轴承信号的阶比分析表明本文方法可以取得较好 tating Machinery Based on Time-Frequency Analysis [Disserta- 效果 tion].Chongqing:Chongqing University,2003:2 (郭瑜.基于时一频分析的虚拟式旋转机械特征分析仪系统 参考文献 的研究[学位论文].重庆：重庆大学，2003：2) [1]LI H,Zheng H Q,Yang S P.Gear fault diagnosis based on ampli- [11]Xu MQ.Huang W H,Zhang JZ.Application of haar wavelet on tude and phase demodulation during run-up.J Vib Shock,2008, analysis of vibration signal of rotating machinery in fast run up 27(2):8 state.J Vib Eng,2000,13(2):216

第 10 期 李修文等: 基于移频技术的短时傅里叶变换阶比分析 理． 根据起始转速，取 f0 = 20 Hz，允许波动范围 fa = 2 Hz，对时频面进行转速跟踪，最终得到的转速值如 图 13( b) 所示． 由于是减速过程( 近似) ，转速曲线 拟合次数选择 2，拟合结果如图 13( b) 所示． 图 13 转速跟踪结果． ( a) 方差 σ 随次数 M 的变化趋势; ( b) 转速跟踪值与拟合结果 Fig． 13 Results from speed tracking: ( a) trend of variance σ with times M; ( b) speed tracking value and its fitting result 根据转速曲线对原始信号进行等角度重采样， 阶比分析的结果如图 14 所示． 从图可以看到明显 的外圈故障频率阶比 3. 6 以及转频成分，并且还可 以发现外圈故障频率带有转频的边带． 图 14 基于移频技术的短时傅里叶变换阶比分析结果 Fig． 14 Result from order analysis based on FS-STFT 3 结论 ( 1) 提出了一种基于移频技术的短时傅里叶 变换分析方法． 在传统短时傅里叶变换的基础上， 利用傅里叶变换的移频性质，根据一定的频率间 隔依次将各频率成分进行搬移，减少能量泄露，算 法高效的同时，可以大大提高时频分析的频率分 辨率． ( 2) 研究了无转速信号下的阶比分析方法． 首 先对基于移频技术的短时傅里叶变换分析的结果进 行局部阈值降噪，不同时刻选取不同的阈值，利用峰 值搜索并跟踪转速值，之后采用最小二乘多项式拟 合得到转速曲线，最终用于阶比分析． 对仿真和实 际轴承信号的阶比分析表明本文方法可以取得较好 效果． 参 考 文 献 ［1］ LI H，Zheng H Q，Yang S P． Gear fault diagnosis based on ampli￾tude and phase demodulation during run-up． J Vib Shock，2008， 27( 2) : 8 ( 李辉，郑海起，杨绍普． 基于幅值和相位解调分析的齿轮箱 起动过程故障诊断． 振动与冲击，2008，27( 2) : 8) ［2］ Lopatinskaia E，Zhu J，Mathew J． Monitoring varying speed ma￾chinery vibrations: Ⅰ． The use of non-stationary recursive filters． Mech Syst Signal Process，1995，9( 6) : 635 ［3］ Bossley K M，McKendrick R J，Harris C J，et al． Hybrid compu￾ted order tracking． Mech Syst Signal Process，1999，13( 4) : 627 ［4］ Zhao X P，Hou R T． Gabor order tracking based on Viterbi algo￾rithm． J Mech Eng，2009，45( 11) : 247 ( 赵晓平，侯荣涛． 基于 Viterbi 算法的 Gabor 阶比跟踪技术． 机械工程学报，2009，45( 11) : 247) ［5］ Wang K S，Heyns P S． Application of computed order tracking， Vold-Kalman filtering and EMD in rotating machine vibration． Mech Syst Signal Process，2011，25( 1) : 416 ［6］ Gao Y，Guo Y，Chi Y L，et al． Order tracking based on robust peak search instantaneous frequency estimation． J Phys Conf Ser， 2006，48: 479 ［7］ Liu W B，Guo Y，Li Z X． Removal of crossing-items of WVD based on gabor expansion． J Vib Shock，2008，27( 10) : 121 ( 刘文彬，郭瑜，李之雄． 基于 Gabor 展开的 Wigner-Ville 分布 的交叉项消除． 振动与冲击，2008，27( 10) : 121) ［8］ Jones D L，Parks T W． A high resolution data-adaptive time-fre￾quency representation． IEEE Trans Acoust Speech Signal Process， 1990，38( 12) : 2127 ［9］ Kong Q P，Song K C，Chen Y． Study of time-frequency order tracking of vibration signals in engine speed changing stage． J Vib Eng，2005，18( 4) : 448 ( 孔庆鹏，宋开臣，陈鹰． 发动机变速阶段振动信号时频分析 阶比跟踪研究． 振动工程学报，2005，18( 4) : 448) ［10］ Guo Y． Research on Virtual Characteristics Analyzer System of Ro￾tating Machinery Based on Time-Frequency Analysis ［Disserta￾tion］． Chongqing: Chongqing University，2003: 2 ( 郭瑜． 基于时--频分析的虚拟式旋转机械特征分析仪系统 的研究［学位论文］． 重庆: 重庆大学，2003: 2) ［11］ Xu M Q，Huang W H，Zhang J Z． Application of haar wavelet on analysis of vibration signal of rotating machinery in fast run up state． J Vib Eng，2000，13( 2) : 216 ·1195·

·1196· 北京科技大学学报 第34卷 (徐敏强，黄文虎，张嘉钟.旋转机械高速启动过程振动信号 为的识别.北京科技大学学报，2007,29（增刊2）：144) 分析方法的研究.振动工程学报，2000,13(2)：216) [13]Wang G D,Li M,Yang J H.et al.Fault pattern recognition of [12]Zhang W J,Xu J W,Li Z M,et al.Bifurcation and chaos iden- rolling bearings based on characteristic waveform sparse matc- tification on the rubbing rotor.J Univ Sci Technol Beijing,2007. hing.J Univ Sci Technol Beijing,2010.32(3):390 29(Suppl2):144 (王国栋，黎敏，阳建宏，等.基于特征波形稀疏匹配的滚动 (张武军，徐金梧，吕志民，等.碰摩转子系统分叉与混沌行 轴承故障模式识别.北京科技大学学报.2010,32(3)：390)

北 京 科 技 大 学 学 报 第 34 卷 ( 徐敏强，黄文虎，张嘉钟． 旋转机械高速启动过程振动信号 分析方法的研究． 振动工程学报，2000，13( 2) : 216) ［12］ Zhang W J，Xu J W，Lü Z M，et al． Bifurcation and chaos iden￾tification on the rubbing rotor． J Univ Sci Technol Beijing，2007， 29( Suppl 2) : 144 ( 张武军，徐金梧，吕志民，等． 碰摩转子系统分叉与混沌行 为的识别． 北京科技大学学报，2007，29( 增刊 2) : 144) ［13］ Wang G D，Li M，Yang J H． et al． Fault pattern recognition of rolling bearings based on characteristic waveform sparse matc￾hing． J Univ Sci Technol Beijing，2010，32( 3) : 390 ( 王国栋，黎敏，阳建宏，等． 基于特征波形稀疏匹配的滚动 轴承故障模式识别． 北京科技大学学报． 2010，32( 3) : 390) ·1196·