基于瞬时频率估计与Vold-Kalman滤波的铣削颤振识别

汪晓姗， 彭志科， 陈是扦

(上海交通大学 机械系统与振动国家重点实验室，上海 200240)

颤振是金属加工中出现的自激振动，严重制约了加工精度和加工效率。随着机床加工柔性化的日益发展，要求加工能对不同的工件在不同的条件下进行，因此不能从根本上杜绝颤振。近年来，于英华等[1-2]很多国内外学者开始着力于研究颤振的在线监测和预报控制。

铣削颤振的形成是一个孕育的过程，在颤振早期，监测系统抢先把铣削过程即将发生颤振的预兆先通报给控制系统，是实现颤振预报的关键。各种各样的传感信号可以被用来采集颤振信号，比如动态铣削力信号[3-4]，加速度信号[5]，力矩信号[6]等。Kuljanic等[7]提出了利用多传感器来监测识别颤振，这种方法提高了颤振监测的精确性和鲁棒性；Tansel等[8]通过铣削位移信号的谐波分量来监测颤振；Tarng等[9]在发生颤振时切削力信号在的频谱上的窄带特性来预测颤振。

工程中，直接根据传感信号难以识别颤振，因此需要从采集来的信号中获得对故障敏感的特征量。时域中主要采用的是概率统计分析，如互相关系数[10]、标准差[11]等；频域中，迟玉伦等[12]基于通过识别铣削力的谱峰频率来监测颤振。然而，这些传统的方法只适合分析平稳信号，而且无法准确预报颤振的发生。对于成分复杂的颤振信号，需要将其分解成单分量再进行分析。Cao等[13]应用EMD分解铣削力信号，提取颤振的特征，但是无法克服EMD分解的模态混叠。该方法易产生虚假的信号分量，干扰故障特征的诊断；Zhang等[14]提出用VMD分解铣削力信号，但是VMD本质上通过自适应的频域滤波器组来提取信号，其分解铣削力信号得到的子分量没有物理意义。Vold-Kalman滤波方法[15]基于数据方程误差和最小化结构，有效地进行多分量信号分解。不仅可以避免EMD分解带来的模式混淆，而且可以有效的从复杂的颤振信号中提取单分量。Vold-Kalman滤波提取单信号首先需要根据时频脊线估计单分量信号的瞬时频率。经典的瞬时频率估计方法中，过零法和相位差分法较为简单，但是误差大。HHT继承了EMD的缺点。在参数化时频分析的基础上，Yang等[16]提出了基于优化频谱集中性指标，对信号瞬时频率参数估计的方法，该方法具有抗噪性强、精确度高等特点。

颤振产生时，信号的频率成分和能量分布发生改变。稳定铣削时，颤振的能量分布在噪声、主轴的转速频率及其谐波处。发生颤振时，振动的能量逐渐转移到颤振频率处。因此本文提出基于频域能量分布改变的颤振识别方法。针对复杂的铣削颤振信号，考虑Vold-Kalman滤波方法在单分量提取和优化频谱集中性指标在瞬时频率计算方面的独特优势，本文提出了基于两者的多分量信号分解方法。采用分解后的子信号在频域上的能量分布作为颤振识别的特征。多次铣削实验验证了本文所提出的方法的有效可行。

1 多分量信号分解

1.1 基于优化频谱集中性指标估计信号瞬时频率

该过程可以描述为，选取一个目标分量，采用合适的变换核参数构造对应的旋转算子，将原信号旋转变换到它的一个旋转域。为了在上述旋转域评估分量能量的集中性，提出频谱集中性指标。旋转域中能量越集中，集中性指标越大，可认为变换核参数与分量模型的相位参数越相近。

单分量解析信号可用多项式相位信号形式表示[17]

(1)

(2)

sd(t;C)=s(t)Φ(t;C)

(3)

sd(t;C)=a(t)exp(j(2πc0t+φ0))

(4)

此时能量集中于c0处， 为衡量频谱集中性，定义频谱集中性指标

(5)

式中：E(·)表示期望算子； F (·)表示傅里叶变换，以SCI最大为参数估计的指标。

(6)

本文，估计参数采用粒子群优化算法，通过找到解调信号的最大谱峰值对应的频率来估计初始频率参数

(7)

通过式(6)、式(7)估计瞬时频率的所有参数(c0,…,ck)。

1.2 Vold-Kalman滤波

Vold-Kalman滤波方法可以从复杂多分量信号中提取有效的频率成分。与传统的滤波方法相比，该方法避免了由时域至频域变换带来的相位偏差，可以将复杂多分量信号分解为单分量。Vold-Kalman滤波的中心频率能根据瞬时频率进行自适应调节，因此能够有效分离在时频域内邻近甚至交叉的信号分量[18]。

调制信号可以表示为

(8)

式中：Ak(t)为第k个分量的幅值包络；Θk(t)为载波信号；k为阶次。

(9)

振幅包络线是一个低频率调制的载波信号，可以表示为一个低阶多项式。对于离散信号

(10)

Ak(n-1)-2Ak(n)+Ak(n+1)=εk(n)

(11)

对于实际信号，Ak(0)=0。 当n=1时，Ak(0)-2Ak(1)+Ak(2)=εk(1)； 当n=N时， 其中是信号的长度。 式(11)的矩阵的形式

(12)

矩阵形式为

MA=ε

(13)

实际测试信号

(14)

为噪声或误差项，式(14)可以写成矩阵形式

X-BA=δ

(15)

优化频谱集中性指标， 估计信号的瞬时频率ωk(n)， 从而得到B， 根据(13)和(16)得到各信号分量的幅值包络矩阵A， 由矩阵C和A可以重构各信号分量

S=AC

(16)

信号的非线性调频分量分解的步骤如下：

(1) 优化频谱集中性指标，估计分量的瞬时频率；

(2) 根据估计到的瞬时频率应用Vold-Kalman滤波器滤出分量ui；

(3) 从信号中减去复原的信号分量作为初始信号，重复以上步骤，得到信号的各个分量。

2 特征提取

为了颤振预报及时并且精确，正确提取特征是关键。所选择的特征参数要充分反映铣削颤振过程的本质，同时还必须考虑信号采集和数据处理的简易可行，使监测系统可以在几秒时间内完成计算工作，向控制系统发出预报信号。

频带熵、信息熵、近似熵等的引用使信号复杂性表征成为可能，但是它们都存在缺点：在单一尺度上对时间序列进行分析、需要比较自身数据段且参数改变时结果一致性差。因此，本文采用多分量信号分解和熵结合的方法。在稳定铣削时，除了主轴转速频率及其谐波处，其他频率处具有能量小、分布相对平均、不确定的特点。当发生颤振时，除了主轴转速频率及其谐波处，在相应的频带内就会出现相应共振频率，能量便会集中在此频率带内，使能量分布的不确定性减少。因此可以通过计算分解得到的各子信号的能量分布以检测颤振。信号经过分解可以得到n个单分量u1,u2,…,un， 计算出分别对应的能量

(17)

假设残余分量可以忽略，由于本文的分解方法具有正交性，n个子信号的能量之应该等于原始信号的能量。由于各个子信号包含不同的频率成分和能量，从而形成在铣削力信号在频域的能量分布。由此引入能量熵铣削过程中的振动信号包括：由于铣刀旋转引起的周期信号部分、由于再生效应引起的颤振和由于系统噪声等引起的随机扰动部分。从实验信号中提取颤振分量对分析颤振有重要意义。本文提出的颤振监测方法如图1(b)所示。首先，动态铣削力信号经过Vold-Kalman滤波器去除由于主轴旋转产生的信号分量。然后优化频谱集中性指标对信号进行瞬时频率估计。Vold-Kalman滤波将复杂的多分量信号分解为单分量，从时域中提取对分析颤振有意义的信号分量。最后计算分量信号的能量熵作为识别颤振的特征。

(18)

3 实验设置

实验装置图1(a)所示，所有的铣削实验都是在数控铣削机床上完成的，动态铣削力信号测量系统由Kistler 9272测力计、Kistler 5070A电荷放大器、5697数据采集卡和PC组成，采样频率为20 000 Hz。工件材料为400 mm←100 mm←100 mm的铝合金， 刀具为四齿硬质合金球头铣刀， 半径4 mm， 悬伸长度为40 mm。 颤振的产生与铣削参数有直接的关系， 因此可以通过改变主轴转速n、进给速度f、轴向铣削深度ap获取感兴趣的实验数据。

本文进行了5次实验，实验条件如表1所示。 Test1、Test2、Test3、Test4构成第一组实验； Test5分别在稳定、过渡、颤振时采样，这3段采样信号构成第二组实验。

图1 实验装置及方法流程图Fig.1 The experimental platform and the flow chart

状态n/(r·min-1)f/(mm·min-1)ap/mm第一组Test1稳定4 2005806.6Test2稳定3 6002502.6Test3颤振4 0005806.6Test4颤振4 0006406.7第二组Test5过渡4 0005806.3

4 结果分析

对于第一组实验，分别对Test1、Test2、Test3、Test4进行多分量信号分解，计算子信号的能量在频域上的分布。对比分析在不同实验条件下，稳定和颤振的子信号的能量熵。

图2(a)为Test5采集到的进给方向动态铣削力的时域波形，Test5经历了从稳定铣削过渡到颤振的过程，实验前0.6 s铣削稳定，0.6 s～0.8 s信号振幅急剧增大， 0.8 s开始颤振。 分别自0.2 s、0.8 s、2 s起， 往后截取5 000个采样点长的信号记作Sample1(稳定)、Sample2(过渡)、Sample3(颤振)，作为第二组实验。对比Sample1、Sample2、Sample3经过信号分解得到的子信号的能量分布，分析在实验条件不变的情况下，随着铣削从稳定过渡到颤振，子信号在频域的能量熵的变化。

图2 Test5及各采样点时域波形Fig.2 Time domain waveform of test5 and sampling points

4.1 数据分析

第一组原始实验数据的时域波形和频谱如图 3所示，Test1(图3(a))、Test2(图3(b))铣削稳定，频谱由主轴转速频率及其谐波、噪声组成；改变铣削条件，动态铣削力振幅急剧增加，发生颤振(图3(c))，频谱上的出现能量远远高于其他频率的颤振峰值；继续增加铣削深度，颤振的更加剧烈，如图3(d)。

图3 Test1、Test2、Test3、Test4的波形和FFTFig.3 Vibration signal in time domain andfrequency domain of Tests

第二组原始实验数据的频谱如图4(a)、图4(b)、图4(c)所示，f表示主轴转速频率，fc表示颤振频率。铣削稳定时(图4(a))，频谱由主轴转速及其谐波成分和分布不确定的噪声组成；过渡状态(图4(b))频谱除了主轴转速谐波成分和噪声，还出现了颤振分量和它的各阶谐波成份，此时工件表面还没有明显的颤振；随着颤振程度增加(图4(c))，颤振频率及其谐波对应的峰值也急剧增长，信号的能量除了集中在主轴转速频率处，都汇集在fc及其谐波上，且颤振峰值还再继续增加。

在分解信号、分析子信号能量频域分布之前，为了避免主轴转速及其谐波分量对颤振识别产生的干扰，先估计这些分量的瞬时频率，用Vold-Kalman滤波提取此分量信号并从原始信号中剔除这一部分。对预处理过的铣削信号进行FFT，稳定铣削时频谱(图4(d))只剩下在频率域上分布均匀且峰值较低的随机振动；过渡到颤振时(图4(e)、图4(f))， 频谱在fc及2倍的fc处出现能量很高的峰值，除此之外还有能量比较低的噪声。

通常动态铣削力信号具有非线性、非平稳等特征，直接进行傅里叶变换无法揭示频率分量随时间的变化规律；短时傅里叶变换、小波变换等能够建立信号在时间和频率上的分布，然而它们是先验性的自适应能力差。并且它们只能用来离线判断铣削状态，无法运用在实际中。预报铣削颤振要准确、及时，一旦出现颤振发生的预兆，计算机可以识别并且迅速的启动控制系统，调整转速、铣削用量等，将颤振控制在孕育阶段。故障振动信号的能量分布会随着颤振发生相应的改变，故采用分解后子信号在频域分布的熵值的范围可用来判断是否颤振。

图4 Sample1、Sample2、Sample3滤波前后的FFTFig.4 FFT of Sample1、Sample2、Sample3 beforeand after filtering

首先，Vold-Kalman滤波器波除去信号的转速频率及其谐波成分。优化频谱集中性指标来估计发生颤振时预处理过的信号的瞬时频率参数。用估计到的参数将非线性调频信号解调为平稳信号。最后Vold-Kalman滤波器提取解调单分量信号。本文对实验信号精细分解，按照分解后单分量子信号能量从高到低，得到这些分量集中度较好的频谱表示。由于主要的故障信息集中在前几个子信号中，剩余信号能量远小于原始信号总能量，可以忽略，因此本文选用了前8个子信号，8个子信号的能量之和应该恒等于原始信号的总能量。由于篇幅限制，本文仅展示预处理后的第二组信号的Sample1、Sample3的分解结果。Sample1的经过分解后的各子信号按能量从高到低如图5(a)所示，均为随机分布的噪声，能量分布相对平均和不确定；Sample3经过分解后的各子信号能量从高到低如图5(b)所示，其中u1，u2能量远高于其他各子信号能量，信号频谱的幅值范围均在0～4。u1的峰值正好在铣削系统的固有频率附近，证明本实验出现的是典型的再生颤振。u2为颤振频率的2倍频。剩余各阶子信号是被齿频调制的颤振分量或噪声，其时域波形的振幅小，能量明显低于u1和u2。

经过预处理的动态铣削力信号经过分解后得到多个分量，且各个分量的瞬时频率在频谱上不相交，子信号的能量随分解阶数的增加而降低。在信噪比低的情况下，基于频谱集中性指标估计分量的瞬时频率较精确，抗噪性强。即使对于能量较小的单分量的瞬时频率参数的估计也比较准确。

图5 本文方法分解结果Fig.5 Decomposition results of this method

4.2 颤振识别

随着颤振的发生，振动的能量发生转移，所以振动的概率分布也发生了变化，这些都可以用熵来进行判断。由于各阶子信号含不同的频率成分且具有不同的能量，分别求它们的能量E={E1,E2,…,En}。首先根据式(18)将各子信号的能量归一化。分析第一、二组实验经过分解后的能量分布，分别如图6(a)，图6(b)所示。对于第一组实验，将分解得到的子信号归一化，u1～u8的能量从高到低，分布相对均匀(图6(a))，熵值较大；而Test3能量集中在u1和u2，所以熵值相对较小；Test4颤振更加剧烈，随着颤振的程度增加，u3～u8所占能量几乎为零，所以计算得到的熵值最小。计算Test1、Test2、Test3、Test4的能量熵结果如图7(a)所示。对于第二组实验，短暂的稳定铣削过程中，子信号能量分布较为均匀。随着颤振的出现，能量急剧转移到u1和u2，颤振剧烈时，u1和u2汇集了几乎全部的能量。根据式(19)，计算Sample1、Sample2、Sample3的能量熵结果如图7(b)所示，可以看出铣削从稳定过渡到颤振，熵的值急剧减小。第二组的分析更加证明了本文所提出的颤振识别方法的可靠性。

图6 第一组、第二组归一化后各子信号能量分布Fig.6 Normalized energy ration of the sub-signals

稳定铣削时，预处理后的铣削力信号经分解，得到的子信号能量分布较平均和不确定，能量熵较高。当出现颤振时，在相应的频带内就会出现相应的颤振频率。此时，能量集中在颤振频率上，使能量分布的不确定性减小，从而熵减小。通过第一组、第二组实验的分析，可见本文提出的方法可以有效识别颤振。然而由于实验条件的限制，本文没有给出发生颤振时能量熵的阈值。

图8所示为Test5从稳定铣削到颤振的过程，可见在0.6 s颤振发生时，能量熵明显变小，证实了能量熵可以作为监测颤振的指标。本文的不足之处在于，由于实验前后时间跨度大，无法提供工件表面加工质量前后对比照片。

图7 第一组、第二组能量熵Fig.7 Energy entropy of the tests

图8 Test5能量熵随铣削力的变化Fig.8 The vibration signal and the relatedenergy entropy in Test 5

5 结 论

本文提出了一种简单有效的铣削颤振识别方法。动态铣削力信号具有不平稳、频率时变、噪声大等特点，分解信号前先经过Vold-Kalman滤波器去除信号中的主轴转速及其谐波分量，再通过优化频谱集中性指标来估计信号瞬时频率参数，用估计到的参数将非线性调频信号解调为平稳信号，最后用Vold-Kalman滤波器提取单分量信号，并得到能量从高到低的各阶子信号。该信号分解方法抗噪性好，对实验信号分解精细，并得到这些分量集中度较好的频谱表示。随着铣削条件变化，经过分解得到的各子信号，其能量分布也呈现变化。因此，本文提出表征能量分布特征的能量熵作为颤振识别的特征。在稳定铣削时，实验信号进行预处理后，其分解得到的各阶子信号主要为噪声，因此能量熵小。随着颤振越来越剧烈，振动的能量集中到含有颤振信息的子信号上，因此能量分布从分散变得集中起来，能量熵变大。经过实验验证，该方法简单有效，尤其在颤振发生的早期识别精度高，是实现颤振在线实时控制的关键。