无转速计的旋转机械Vold-Kalman阶比跟踪研究

邓 蕾，傅炜娜，董绍江，汤宝平

(重庆大学 机械传动国家重点实验室，重庆 400030)

阶比跟踪(OT)技术是旋转机械升降速过程非平稳振动信号分析的有效方法，在旋转机械的非平稳振动信号分析和故障诊断中占有重要地位。

1993年Vold和Leuridan［1］首次提出了基于卡尔曼滤波的阶比跟踪的方法，开创了通过时域滤波方法进行阶比跟踪的新思维。目前，丹麦的B＆K公司将VKF-OT方法作为其阶比分析商业化软件的主要技术，并得到了广泛的应用。文献［2］对VKF-OT理论作了进一步的详细研究，提出了角速度和角位移两种VKFOT方法，对两者之间的区别进行了详细地描述。在此基础上文献［3］提出和实现了自适应VKF-OT方法。另外，国内对此也作了一些研究［4，5］。

然而，VKF-OT方法仍需要通过安装转速计来获取转速信号，在不方便安装鉴相装置的场合限制了应用。针对简化阶比跟踪的实现途径和对硬件安装要求的问题，文献［6］提出了基于短时傅里叶变换的峰值估计瞬时频率的阶比跟踪方法，结合Gabor时频滤波实现了阶比分量的提取。文献［7］提出了基于STFT-VA瞬时频率估计的自适应VKF-OT方法。但这些方法主要是在旋转机械振动信号的时频谱图中通过峰值搜索获得某阶比分量的瞬时频率，进而得到所需要的参考轴转速，所以分析精度受到频率分辨率的限制，只有在特定条件下才能应用于多分量信号。

本文从瞬时频率估计理论出发，结合离散频谱校正技术中的能量重心法，通过研究两者之间的联系，改进了瞬时频率估计方法，再将其应用于VKF-OT方法中，提出了一种无转速计的基于能量重心法瞬时频率估计的VKF-OT方法。

1 Vold-Kalman阶比跟踪原理

阶比定义为参考轴每转内发生的旋转振动次数，可以设定为一个幅值和频率随时间变化的正弦函数，阶比的频率和参考轴的旋转频率相关联。阶比信号可以表示为幅值和载波的乘积［4］:

式中，k为基本频率或参考轴转速的倍数，表示为被跟踪的阶比;ak(t)表示复包络，代表了阶比幅值的变化，a－k(t)是 ak(t)的复共轭 ak(t)=a－k(t)*。Θk(t)为载波，表示为:

式中，ω(τ)为参考轴的角频率，∫t0ω(τ)dτ 为角位移。式(2)的离散形式表示为:

VKF-OT的主要目标是获得复包络ak(t)或者离散形式的ak(n)

由式(1)知，复包络ak(t)是载波Θk(t)的低频幅值调制，由于系统的固有惯量，使得阶比幅值平滑变化，可以用低阶多项式表示，则系统的状态方程为:

式中，▽表示不同的算子，s为给定阶数，ε(n)为非一致项。

实际测得的振动信号y(n)由各阶比分量xk(n)的总和加上测量误差和噪声ξ(n)组成，则系统的观测方程为:

为了跟踪某一阶比分量xk(n)，用VKF-OT方法来计算复包络ak(n)，则状态方程和观测方程给出以下形式:

根据式(6)，展开一阶状态方程:

其矩阵形式为:

其中，M为N×N－1的矩阵，N为采样点数。同样，根据式(7)，系统的观测方程为:

其矩阵形式为:

根据系统的状态方程(8)和观测方程(10)，通过最小二乘法，使得非一致相ε(n)以及测量误差和噪声ξ(n)的平方和最小［8］，求得未知量 ak(n)。

2 基于能量重心校正法的瞬时频率估计

2.1 瞬时频率原理

瞬时频率(IF)是非平稳信号分析中的一个重要概念，Ville于1948年提出的瞬时频率定义式是目前在学术界最为常用且得到了普遍认可［6，9］，即若实信号为s(t)=a(t)cosφ(t)，则 IF 定义为:

另一种等价形式为时频分布(TFD)的一阶矩，即:

2.2 能量重心离散频谱校正法

能量重心校正法是由三点卷积幅值校正法发展起来的。该方法可以同时对离散频谱的频率、幅值和相位进行校正［10］。能量重心法是通过利用各种窗函数离散频谱的能量重心无穷逼近坐标原点的特点进行离散频谱校正的。这种方法在分析平稳或准平稳信号时，具有较高的分析精度。

设Gi为加Hanning窗的谐波信号的功率谱第i条谱线值，Gk为主瓣内谱线最大值，k为幅值最大点对应的谱线号，采样频率为fs，采样点数为N，f0为主瓣中心，得到能量校正法校正频率的通用公式:

根据帕赛瓦定理，设能量恢复系数为Kt，则校正后的幅值为:

离散频率校正方法中的能量法对所有对称窗函数都适合，校正精度与窗函数有关，加Hanning窗时具有较高的校正精度。负频率成分和间隔较近的多频率成分产生的干涉现象所带来的误差对精度的影响小。不考虑信号中噪声的影响，是一种精度较高的近似校正方法。

2.3 基于能量重心法的瞬时频率估计(IFE)的原理

在STFT时频面上进行峰值搜索估计瞬时频率是目前获取旋转机械升降速阶段振动信号的瞬时频率的主要方法［6］。但是这种方法的估计精度依赖于频率分辨率，邻近频率成分的干扰也会对估计精度产生影响，存在时频集聚性不是很理想的问题。另外，在求取多分量信号的IF时，一般通过时频滤波将其它分量遮掩，提取某一分量的IF，依次递推，增加了算法的计算量。

离散频谱校正方法可以改善由于时域截断产生的能量泄漏的问题，提高分析精度。而旋转机械升降速过程的振动信号属于非平稳信号，不能直接使用离散频谱校正方法。因此，根据STFT算法原理，把短时间间隔的信号看成是准平稳信号，就能够利用离散频谱校正方法来校正信号的频谱。这样可以在频域中直接通过振动信号的频谱进行瞬时频率估计，并通过离散频谱校正方法提高分析精度，进而实现阶比跟踪。

比较式(13)和式(14)发现:式(14)是式(13)的近似离散表达形式，也就是说通过能量重心法校正获得的频率近似等于瞬时频率。基于能量重心法IFE的基本思路是根据假定的瞬时频率初始值进行整周期采样，同时做加窗的DFT;然后利用能量重心法对该段信号做离散频谱校正，得到新的瞬时频率;依次迭代重复，最后得到准确的瞬时频率值。该方法的误差主要来自信号的截断误差，所以采用整周期采样来减小误差，同时对信号加窗，以便提高精度和减小各相邻谱线之间的干扰。

下面给出基于能量重心法瞬时频率估计的计算步骤。

(1)对振动信号进行分段，数据长度为L，每M个点为一段，重叠长度为 l，共分为 m段，m∈［0，Num－1］，Num=int［L/(M －1)］=int(L/dN)，式中 int表示对计算结果取整;

(2)旋转机械升降速过程中总有一稳定转速(相对恒定)阶段，比如，最高转速为9000 r/min的降速阶段，对应的一阶最大转频为150 Hz。所以选择最高速阶段作为初始分析段，根据最大转频和该段的频谱成分设定初始转频f0，局部极大值搜索范围q和整周期采样系数K;

(3)根据当前m，计算当前分段数据的起始点位置:p=m×dN。在采样信号中截取M点数据{s(i)}，i∈［p，M －1+p］};

(4)根据f0和系数R计算出整周期采样点数N0，以分析段中点为中心重新选取N0点作为新的分析段，并对其做加窗的DFT，求得搜索范围q内最大频率谱线位置以及其左右邻近的n条谱线的功率值。其中，n的取值与加窗的类型有关，这里加Hanning窗，n=2;

(5)利用能量重心法对搜索到的峰值频率谱线进行频谱校正，得到校正频率f，令f为新的转频，重新计算出整周期采样所需的采样点数N;

(6)如果N=N0或者达到最大迭代次数j时，f为分析段中点所对应的频率估计值，否则，f0=f，转步骤(4);

(7)令m=m+1，返回(3)，直到m=Num－1。

最终根据迭代求得的参考轴转频f(m)，插值拟合后得到瞬时频率拟合值。

3 基于能量重心法瞬时频率估计的VKF-OT

从前面的讨论可以看出，通过能量重心法得到的参考轴转频近似等于瞬时频率，由此可以获得参考轴转速，进而实现阶比跟踪。本文方法的实现过程如图1所示。

获取旋转机械升降速阶段的振动信号，对该振动信号作基于能量重心法的瞬时频率估计，获得参考轴瞬时频率估计值;插值拟合后得到连续的参考轴瞬时频率，并对各个振动采样点进行阶比相位估计;然后，设置加权因子r，选取阶次p等参数，进行Vold-Kalman自适应滤波;最后得到所提取的第p阶信号分量。

图1 基于能量重心法的瞬时频率VKF-OT流程图Fig.1 Flow diagram of VKF-OT based on the improved IFE with energy centrobaric method

4 仿真与应用实例分析

4.1 仿真分析

为了验证算法，设计一个包含两个独立参考轴转速信号的仿真信号来模拟旋转机械升速阶段的振动信号。采样点数N=10 k，采样频率fs=1 kHz。仿真信号的时域波形如图2所示。

式中，η(n)为高斯白噪声［η(n)～N(0，1)］。此仿真信号由两组信号分量组成，第Ⅰ组包含第1阶、第3阶分量，基准频率从15 Hz到150 Hz线性变化，·Δt表示第n个采样时刻第Ⅰ转轴转过的相位角。第1、3阶比分量的幅值分别呈0到10和5到15线性变化。第Ⅱ组为添加了高斯白噪声的正弦信号，旋转频率f2(n)固定为 200 Hz，幅值固定为 10，·Δt表示第n个采样时刻第Ⅱ转轴转过的相位角。

采用基于能量重心法的瞬时频率估计方法对仿真的振动信号进行瞬时频率估计。分析中，加Hanning窗，n=2，相邻两段时移点数为500，结果如图3所示。

图2 仿真信号时域波形图Fig.2 Time waveform of synthetic signal

图3 基于能量重心法的瞬时频率估计Fig.3 IFE based on energy centrobaric method

图3中圆点表示采用基于能量重心法的瞬时频率估计得到的估计值，实线表示设计的参考轴Ⅰ的瞬时频率曲线，虚线为对估计值采用插值拟合后得到的拟合曲线。从图中可以看出，采用基于能量重心法得到的转频理论上近似于瞬时频率，与仿真信号的理想值非常吻合。

分别对基于STFT峰值搜索法和本文方法获得的瞬时频率估计值求相对于理想瞬时频率值的百分比误差。

式中，f(i)为瞬时频率的理想值，～f(i)为f(i)的拟合值，根据式(17)求得，基于STFT峰值搜索得到的瞬时频率估计值相对于理想瞬时频率值的百分比误差ζ1=0.37%，基于能量重心法的瞬时频率估计得到的拟合值相对于理想瞬时频率值的百分比误差ζ2=0.03%，对比结果表明，由本文提出的基于能量重心法IFE误差小，估计的瞬时频率逼近真实的瞬时频率值。取移动步长为1，分别用这两种方法计算得到瞬时频率估计值，取某一段的数据如图4所示。图中结果再次说明，本文方法连续性好，分析精度受频率分辨率的影响小。

图4 能量重心法IFE和STFT瞬时频率估计对比Fig.4 Comparison of IFE by energy centrobaric method and STFT

然后，采用基于能量重心法IFE获得的瞬时频率拟合值作为参考轴对应的瞬时频率，并求得相应的瞬时角位移信号，对振动信号进行VKF-OT分析。本方法中使用一阶自适应Vold-Kalman滤波器分离各个阶比分量。在图5中，图(a)、图(b)分别为第Ⅰ参考轴的第1、3阶比分量提取结果，图(c)为第Ⅱ参考轴的200 Hz常频分量提取结果。图6为两个轴各自阶比分量的幅值，可以看出由VKF-OT分析得到的幅值与设计值相符合。

为了能够清楚地看到提取结果与理想信号之间的差异，取第Ⅰ参考轴的第1阶分量的前1000个点与设计信号进行对比。如图7所示，实线表示用本文方法所提取的第1阶比分量，虚线表示理想信号，从图中可以看到，除了起始点处有一些偏差外，提取的阶比幅值都和设计值吻合。

图5 基于能量重心法瞬时频率估计的VKF-OT各阶比分量提取结果Fig.5 Extracted orders components using VKF-OT based on the improved IFE with energy centrobaric method

图6 VKF-OT分析得到的Ⅰ、Ⅱ轴各阶比分量幅值Fig.6 Amplitudes of orders about axes－ ⅠandⅡtracked by using VKF-OT

图7 第I参考轴的第1阶比分量与理想信号的前1000点对比结果Fig.7 Comparison between order－1 of axis－I and designed signal for n=1000

4.2 应用实例分析

下面对某一偏心直流电机在升速阶段引起的悬臂梁结构的振动信号进行测试。试验装置如图8所示，直流电机安装在简支梁振动台上，振动信号用YD－37型加速计拾取。采样频率为10 kHz，采样长度为75 k点。首先，通过基于能量重心离散频谱校正法的瞬时频率估计法提取振动信号的瞬时频率估计值，并通过三次样条拟合获得瞬时频率曲线。

图8 电机升速振动试验装置Fig.8 Test setup for vibration measurement during motor's run-up

结果如图9所示，圆点表示瞬时频率的估计值，实线表示三次样条插值拟合后的瞬时频率曲线。然后，利用参考轴瞬时转速与其瞬时频率的对应关系，求得参考轴的转速信号。最后，对原始振动信号进行VKFOT分析，提取第4阶比分量，结果如图10所示。分析结果表明，通过本文方法在提取振动信号瞬时频率的同时，能够在时域中直接提取所感兴趣的某阶比分量。

5 结论

图9 能量重心法瞬时频率估计的结果Fig.9 IFE based on energy centrobaric method

图10 偏心电机升速阶段振动信号第4阶比分量提取Fig.10 Extracted order-4 waveform of motor's vibration signal

本文分析了瞬时频率估计理论和离散频谱校正方法中的能量重心法，根据两者之间的理论联系，通过用能量重心法对信号的频谱做校正，得到了信号的瞬时频率，再采用Vold-Kalman阶比跟踪提取感兴趣的某阶分量，从而实现无转速计的旋转机械Vold-Kalman阶比跟踪。因为该方法根据获得的瞬时频率对原分段信号进行了整周期采样处理，使得频率分辨率能够自适应的变化，同时能量重心校正法的引入可以对误差进行了修正，提高了分析精度并减小了各相邻谱线之间的干扰。仿真分析和应用实例分析表明该方法的有效性，在精确提取参考轴转速信号的同时，提取了阶比分量，是对现有技术的一个有力补充。

［1］Vold H，Herlufson H，Mains M.Multi axle order tracking with the Vold-Kalman tracking filter［J］.Sound and Vibration Magazine，1997，13(5):30 －34.

［2］Pan M Ch，Lin Y F.Further exploration of Vold-Kalman filtering order tracking with shaft-speed information-I:Theoretical Part，numerical implementation and parameter investigations［J］.Mechanical Systems and Signal Processing，2006，20:1134－1154.

［3］Pan M Ch，Wu C X.Adaptive Vold-Kalman filtering order tracking［J］.Mechanical Systems and Signal Processing，2007，21:2957 －2969.

［4］孔庆鹏，宋开臣，陈 鹰.最小二乘自适应滤波旋转机械阶比跟踪研究［J］.浙江大学学报(工学版)，2003，40(9):1648－1651.

［5］孔庆鹏，刘敬彪，章雪挺，等.多轴机械自适应滤波交叉阶比跟踪［J］.农业机械学报，2009，40(1):213－216.

［6］郭 瑜，秦树人，汤宝平，等.基于瞬时频率估计的旋转机械阶比跟踪［J］.机械工程学报，2003，39(3):32－35.

［7］赵晓平，张令弥，郭勤涛.基于瞬时频率估计的自适应Vold-Kalman阶比跟踪研究［J］.振动与冲击，2008，27(12):112－116.

［8］Feldauer Ch，Holdrich R.Realisation of a Vold-Kalman Tracking Filter-A Least Square Problem［C］.Proceedings of the COST G-6 Conference on Digital Audio Effects(DAFX－410800)，Verona Italy，2000，December，7－9.

［9］Boashash B.Interpreting and estimating the instantaneous frequency of a signal—Part I:Fundamentals［J］.Proceedings of IEEE，1992，80(4):520－538.

［10］丁 康，江利旗.离散频谱的能量重心校正法［J］.振动工程学报，2001，14(3):354－358.