基于VMD改进算法的铣削颤振仿真识别

Simulation Identification of Milling Chatter Based on Improved VMD Algorithm

摘要: 颤振是机床铣削加工中的关键问题，严重影响加工效率与工件质量，因此早期颤振检测具有重要意义。本文提出一种基于优化变分模态分解与功率谱熵结合的铣削颤振识别方法。针对变分模态分解参数选择难题，采用麻雀搜索算法与最小包络熵相结合的自适应优化策略，实现参数高效寻优；基于分解信号能量比筛选本征模态分量，重构信号以去除噪声干扰；引入功率谱熵作为颤振识别指标，提取仿真信号特征并实现颤振状态的精准识别。实验结果表明，该方法能够有效检测颤振状态。

Abstract: Chatter is a critical issue in machine tool milling, seriously affecting machining efficiency and workpiece quality, so early chatter detection is of great significance. This paper proposes a milling chatter identification method based on optimized variational mode decomposition (VMD) and power spectrum entropy (PSE). To address the challenge of VMD parameters selection, an adaptive optimization strategy combining sparrow search algorithm (SSA) and minimum envelope entropy (EE) is adopted to achieve efficient parameters optimization. Based on the energy ratio of the decomposed signals, intrinsic mode functions (IMFs) are selected to reconstruct the signals and eliminate noise interference. PSE is introduced as a chatter identification index to extract characteristics of simulated signal and achieve accurate identification of chatter states. Experimental results demonstrate that the proposed method can effectively detect chatter.

1. 引言

颤振是机床加工过程中刀具和工件间产生的一种自激振动，是由于加工系统的刚度不足，以及铣削参数的不合理所导致的[1]。颤振会严重影响加工零件的质量和机床刀具的寿命。选择保守的铣削参数进行加工，虽然会减少颤振，但同时加工效率也会降低。因此，及时检测颤振是控制颤振，提高加工质量和效率的关键，在铣削加工中具有十分重要的意义。

颤振检测的关键是快速准确地识别颤振特征[2]。由于早期的颤振信号能量较小，且传感器采集的信号中含有较多噪声，无法直观有效地判断颤振是否发生，需要对原始信号分析处理才能获取颤振信息[3]。因此，信号处理是颤振检测的重要环节，已逐渐成为当前的研究热点。

目前，许多学者在此方面取得重要成果，同时也有很多不足。JING等[4]提出一种基于快速傅里叶变换(Fast Fourier Transform, FFT)和频谱向量的颤振识别方法，但颤振信号具有非线性和不稳定性，FFT难以准确检测颤振的发生。GAO等[5]提出一种基于连续小波变换的薄壁工件铣削颤振检测方法，然而小波变换(Wavelet Transform, WT)的高频系数分辨率较低，难以识别高频颤振。FU等[6]提出一种基于经验模态分解(Empirical Mode Decomposition, EMD)和希尔伯特–黄变换时频谱的颤振识别方法，但由于EMD模态函数频率分布不均匀，且存在模态混叠等问题，在实际应用中受到很大限制。LIU等[7]提出一种基于变分模态分解(Variational Mode Decomposition, VMD)和能量熵的颤振特征提取方法。虽然VMD具有诸多优点，但其分解层数k和惩罚因子α对颤振信号分解的效果有很大影响，需合理选择参数，而麻雀搜索算法(Sparrow Search Algorithm, SSA)能够准确寻得最优解。

因此，针对早期颤振状态难以检测，高频噪声干扰和VMD参数选择等问题，本文提出一种基于SSA-EE-VMD结合功率谱熵(Power Spectrum Entropy, PSE)的颤振识别方法，以实现准确提取颤振特征。

2. 算法理论

本文所提方法首先通过斜面铣削加工实验采集振动加速度信号，并进行FFT分析信号中颤振频率的能量分布特性。然后采用SSA结合包络熵(Envelope Entropy, EE)对VMD的分解参数进行寻优。再利用最佳参数进行VMD分解，并选取能量占比高的本征模态分量(Intrinsic Mode Functions, IMFs)进行信号重构以消除噪声干扰。最后提取信号中的PSE作为颤振识别指标。

2.1. 变分模态分解原理

VMD作为新兴的自适应信号处理方法，能够把信号分解成含有稀疏特性的分量，从而快速高效地提取非平稳信号中的特征。VMD比EMD具有更为严格的理论基础，其本质是构造和求解变分问题[8]。

VMD的具体分解过程如下：

步骤1：构造原始信号$X\left(t\right)$ 的变分问题，如公式(1)所示，其同时满足公式(2)的条件。

$\underset{u_k,{\omega }_k}{\min }=\left\{{{\displaystyle \sum _k\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)u_k\left(t\right)\right]{\text{e}}^{-j{\omega }_{k}t}\Vert }}_2\right\}$ (1)

$s.t.{\displaystyle \sum }_k{u}_k=X\left(t\right)$ (2)

式中，$k$ 为分解的模态数目，$u_k\left(t\right)$ 为分解的第$k$ 个模态分量，${\omega }_k$ 为估计的中心频率，${\partial }_t$ 为原始信号$X\left(t\right)$ 关于$t$ 的偏导数，$\delta \left(t\right)$ 为狄拉克函数。

步骤2：用拉格朗日乘子$\gamma \left(t\right)$ 和二次惩罚因子$\alpha$ 构造公式的增广拉格朗日函数，如公式(3)所示。

$\begin{array}{l}L\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\gamma \right)=\alpha {{\displaystyle \sum _k\Vert \partial t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}{u}_k\left(t\right)\right)\right]{\text{e}}^{-j{\omega }_{k}t}\Vert }}_2^2+{\Vert X\left(t\right)-{\displaystyle \sum _k{u}_k\left(t\right)}\Vert }_2^2\\ +\left(\gamma \left(t\right),X\left(t\right)-{\displaystyle \sum _k{u}_k\left(t\right)}\right)\end{array}$ (3)

步骤3：用交替方向乘子法对公式(3)迭代求解。

步骤4：重复步骤3过程，直至结果满足公式(4)，此时$u_k\left(t\right)$ 即为所求模态分量。

${\displaystyle \sum _z\frac{{\Vert u_k^{z+1}-u_k^z\Vert }_2^2}{{\Vert u_k^z\Vert }_2^2}}<\epsilon$ (4)

式中，$z$ 为步骤3的迭代次数，$\epsilon$ 为收敛误差。

2.2. 麻雀搜索算法寻优

2.2.1. 麻雀搜索算法原理

SSA算法是一种模拟麻雀觅食行为的启发式优化算法，具有全局搜索能力强、适应性广、收敛速度快和寻优精度高等特点[9]。

SSA将麻雀划分成发现者、加入者与警戒者，位置更新分别如公式(5) (6) (7)所示。

发现者位置：

$x_{i,j}\left(t+1\right)=\{\begin{array}{l}{x}_{i,j}\left(t\right)\cdot \exp \left(\frac{-i}{Q\cdot T}\right),R_2<ST\hfill \\ x_{i,j}\left(t\right)+M\cdot W,R_2\ge ST\hfill \end{array}$ (5)

式中，$t$ 为实时迭代次数，$T$ 为最大迭代次数，$x_{i,j}\left(t\right)$ 为第$j$ 维的第$i$ 只麻雀的位置，$Q$ 为$\left[0,1\right]$ 的随机值，$M$ 为正态分布的随机值，$W$ 为元素均为1的矩阵，$R_2$ 为预警值，$ST$ 为安全值。

加入者位置：

$x_{i,j}\left(t+1\right)=\{\begin{array}{l}M\cdot \exp \left(\frac{x_w\left(t\right)-x_{i,j}\left(t\right)}{i^2}\right),i>0.5g\hfill \\ x_{bj}\left(t+1\right)+\left|x_{i,j}\left(t\right)-x_{bj}\left(t+1\right)\right|\cdot A^{+}\cdot W,i\le 0.5g\hfill \end{array}$ (6)

式中，$x_w\left(t\right)$ 为第$t$ 次迭代的最差位置，$x_{bj}\left(t+1\right)$ 为第$t+1$ 次迭代的最佳位置，$A^{+}$ 为元素1或−1的矩阵，$g$ 为麻雀数量。

警戒者位置：

$x_{i,j}\left(t+1\right)=\{\begin{array}{l}{x}_{bj}\left(t\right)+\beta \left|x_{i,j}\left(t\right)-x_{bj}\left(t\right)\right|,f_i>f_{\text{g}}\hfill \\ x_{i,j}\left(t\right)+D\left(\frac{\left|x_{i,j}\left(t\right)-x_w\left(t\right)\right|}{\left(f_i-f_w\right)+e}\right),f_i=f_{\text{g}}\hfill \end{array}$ (7)

式中，$\beta$ 为步长限定参数，$D$ 为麻雀搜索方向的随机值，$e$ 为防止分母为0的常数，$f_i$ 为麻雀实时适应度值，$f_{\text{g}}$ 和$f_w$ 分别为实时最佳和最差适应度值。

2.2.2. 适应度函数及优化步骤

SSA算法需根据适应度函数值的大小来评判最优参数。而EE可以有效反应信号的稀疏特性，熵值越小则稀疏性越强，也表示信号中的颤振信息越充足[10]。因此，以EE作为SSA的适应度函数，其计算过程如公式(8)。

$\{\begin{array}{l}{E}_p=-{\displaystyle \sum }_{j=1}^J{p}_j\log \left(p_j\right)\hfill \\ p_j=\frac{a\left(j\right)}{{{\displaystyle \sum }}_{j=1}^{J}a\left(j\right)}\hfill \\ a\left(j\right)=\sqrt{x^2\left(j\right)+{\widehat{x}}^2\left(j\right)}\hfill \end{array}$ (8)

式中，$E_p$ 为$p_j$ 的包络熵，$J$ 为信号个数，$p_j$ 为$a\left(j\right)$ 的归一化表示形式，$a\left(j\right)$ 为零均值信号，$x\left(j\right)$ 为包络信号，$\widehat{x}\left(j\right)$ 由$x\left(j\right)$ 希尔伯特变换得到的信号。

为了能够搜索到VMD参数的最优解，以IMFs的平均包络熵(Average Envelope Entropy, AEE)最小为全局优化目标，AEE值的计算过程如公式(9)所示。

$\text{minAEE}=\text{min}\frac{1}{A}{\displaystyle \sum }_{i=1}^A{E}_{pi}$ (9)

式中，$A$ 为本征模态分量个数，$E_{pi}$ 为各本征模态分量的包络熵值。

SSA优化VMD关键参数的具体步骤为：

步骤1：初始化各参数。如参数k和α的取值范围、发现者和警戒者的比例、示警阈值、麻雀数量、迭代次数等；

步骤2：计算寻优过程的适应度函数即平均包络熵的值；

步骤3：计算预警值并更新发现者、加入者、警戒者位置；

步骤4：根据更新得到的适应度函数值来确定麻雀最优位置；

步骤5：判断是否满足停止条件，满足后输出VMD的最优参数组合。

2.3. 功率谱熵

PSE是香农信息熵在频域方面的扩展，是衡量信号频谱分布的指标[11]。其基于信号在各频率上的功率分布来计算熵值，具体步骤如下：

步骤1：计算信号的功率谱。对于长度为$N$ 的信号段$s=\left[s\left(1\right),s\left(2\right),\cdots ,s\left(N\right)\right]$ ，其功率谱可表示为公式(10)。

$P_s\left(w_i\right)=\frac{1}{2\pi N}{\left|s\left(w_i\right)\right|}^2$ (10)

式中，$s\left(w_i\right)$ 为信号$s$ 的第$i$ 个FFT，$w_i$ 为第$i$ 个对应频率。

步骤2：标准化信号功率谱。将信号功率谱归一化为概率分布，其概率密度函数如公式(11)所示。

$P_i=\frac{P_s\left(w_i\right)}{{{\displaystyle \sum }}_{i=1}^N{P}_s\left(w_i\right)},i=1,2,\cdots ,N$ (11)

式中，$P_i$ 是$P_s\left(w_i\right)$ 的概率密度。

步骤3：计算功率谱熵的值。利用熵的定义公式(12)计算PSE的值$H$ 。

$H=-{\displaystyle \sum }_{i=1}^N{P}_i\ln P_i$ (12)

将结果用因子$lnN$ 进行归一化处理，如公式(13)所示。

$E=\frac{H}{\ln N}=\frac{-{{\displaystyle \sum }}_{i=1}^N{P}_i\ln P_i}{\ln N}$ (13)

归一化后的PSE值$E$ 是在$\left[0,1\right]$ 内的一个没有量纲的指标，其值接近于1表示信号频谱分布均匀且随机，接近于0表示频谱分布集中且规律。

3. 实验研究

3.1. 铣削加工实验

3.1.1. 铣削实验设置

铣削加工实验系统如图1所示，主要由数控机床、加工工件、切削刀具、振动加速度传感器、数据采集仪和电脑组成。铣削实验选在一台FANUC ROBODRILL数控铣床上进行；实验工件采用长宽为140 × 130 mm、高度为62~70 mm线性递增的铝合金7085斜面工件；实验刀具选用四槽整体硬质合金立米尔斯铣刀，详细参数如表1所示。铣削实验类型为轴向变切深加工，采用干式逆铣。

Table 1. Experimental tool parameters

表1. 实验刀具参数

刀具齿数	刀具直径(mm)	刀具长度(mm)	螺旋角(˚)	刀齿长度(mm)
4	8	85	45	25

Figure 1. Milling experiment system

图1. 铣削实验系统

传感器和工件的布局如图1(a)所示，加工件被夹具平整的夹紧在工作台上，铣刀被机床主轴固定，悬伸长度为35 mm，传感器按图示方向安装在主轴底部用于采集铣削过程中的主轴振动信号，X方向设为铣刀的进给方向，传感器型号为Kistler-5429291型三向加速度传感器，采样频率设定为5000 Hz。图1(b)是动态信号采集分析系统，由DHDAS-5922D型数据采集仪和电脑组成，用于存储传感器采集的实验数据。铣削实验的采集现场如图1(c)所示，数控机床进行高速铣削加工时，传感器采集主轴加速度信号，由数据采集仪传输到笔记本电脑端，再通过DHDAS动态信号采集分析软件进行存储。

3.1.2. 铣削加工参数

在颤振稳定性理论中，进给速度只与静态切削厚度有关，对切屑载荷再生机制没有影响，颤振主要受到主轴转速、轴向切削深度、加工系统的动态特性等因素的影响[12]。因此，为得到稳定铣削、轻微颤振、严重颤振3种不同加工状态的加速度时域信号，保持进给速度和径向切削深度恒定，主轴转速设为6000 rpm、8000 rpm、10000 rpm三种，轴向切削深度由0 mm逐渐递增到8 mm，进行多组铣削加工实验。从中选出三组明显发生颤振现象的加速度信号，切削参数如表2所示。

Table 2. Experimental machining parameter

表2. 实验加工参数

序号	主轴转速(rpm)	轴向切深(mm)	径向切深(mm)	进给速度(mm/min)
1	6000	0~8	2	600
2	8000	0~8	2	600
3	10,000	0~8	2	600

3.2. 铣削状态研究

3.2.1. 加工状态划分

选取第2组的加工信号为例进行分析。根据其原始加速度信号的时域波形变化粗略划分铣削状态，如图2所示。

Figure 2. Original acceleration signal

图2. 原始加速度信号

图2中，a部分处于稳定铣削状态，b部分处于轻微颤振状态，c部分处于严重颤振状态。从这三部分各随机选取一段0.2 s的子信号进行对比分析，结果如图3所示。

Figure 3. Time domain diagram of stable milling, slight flutter and severe flutter

图3. 稳定铣削、轻微颤振、严重颤振状态的时域图

图3(a)~(c)分别是稳定铣削、轻微颤振、严重颤振状态的子信号时域图，其中稳定状态下的加速度信号幅值较小且平稳；轻微颤振状态下的信号幅值增大但较为平缓；严重颤振状态下的信号幅值波动剧烈，且呈现非线性、非平稳性特征。

3.2.2. 频率特性分析

为验证时域特征的准确性，将3种状态的子信号分别进行FFT分析，结果如图4所示。

Figure 4. Frequency domain diagram of stable milling, slight flutter and severe flutter

图4. 稳定铣削、轻微颤振、严重颤振状态的频域图

图4(a)~(c)分别是稳定铣削、轻微颤振、严重颤振状态的子信号频域图，其中稳定状态的信号频谱整体幅值较小，频率分布简单，主要集中在主轴旋转频率及其倍频上；轻微颤振状态的信号频谱幅值增大，频率分布逐渐复杂，除主轴倍频之外还出现幅值较小的颤振频率1575 Hz；严重颤振状态的信号频谱均匀，频率分布更加复杂，且出现多个幅值较大的颤振频率。

4. 结果分析

4.1. 基于SSA-EE-VMD的信号分解

通过时频谱分析虽然可以检测到颤振，但是无法及时识别早期颤振，在工业铣削加工中并不实用。为准确提取颤振特征，选取第2组加速度信号进行SSA-EE-VMD算法的分解处理。首先采用SSA对VMD的分解层数k和惩罚因子α进行寻优，设置麻雀种群个数为100，迭代次数50，参数的寻优过程如图5所示。

Figure 5. SSA optimization process for VMD parameters

图5. SSA优化VMD参数过程

根据图5中的AEE值随迭代次数变化曲线可知，在第28次迭代后寻得最小AEE值7.3809，对应的最优参数：分解层数k为6，惩罚因子α为1315，此后迭代曲线逐渐收敛。采用优化后的VMD分解信号，所得的时域和频域图如图6所示。

由图6可知，严重铣削子信号被分解为6个IMF，颤振频率1575 Hz和2084 Hz分别被分离到IMF4和IMF6频带中，且分解后的各个频带之间界限清晰，未出现模态混叠现象。因此，所提SSA-EE-VMD算法分解效果较好，成功将颤振频率单独分离出来，避免了模态混叠。

4.2. 基于PSE的仿真信号颤振特征提取

颤振的特征可以由频率和能量分布变化所体现。稳定铣削时，能量主要集中在主轴旋转频率及谐波上，随着颤振的发生，能量逐渐向颤振频率上转移。因此，将第2组原始加速度信号进行SSA-EE-VMD分解后，计算各IMF的能量占比，如图7所示，并选取前4个能量比较高的IMFs进行信号的重构以消除噪声干扰，重构后的仿真信号如图8所示。再提取仿真信号的PSE作为加工特征指标进行颤振状态识别，并对其进行归一化处理，结果如图9所示。

Figure 6. Time domain and frequency domain diagrams of SSA-EE-VMD

图6. SSA-EE-VMD分解的时域、频域图

Figure 7. Normalized energy proportion of IMF components

图7. IMF分量的归一化能量占比

Figure 8. Reconstruct simulation signal and original signal

图8. 重构仿真信号与原始信号

Figure 9. Normalized power spectral entropy of simulation signal

图9. 重构仿真信号的归一化功率谱熵

PSE本质上是衡量信号中频率成分分布特性的指标，其值越高表示频率分量越分散，反之频率分量越集中。由于不同加工状态下的信号频率分布情况不同，PSE能够从频域角度反应它们的信号特征。在稳定铣削状态下，频谱分布较为简单，频率分量比较集中，此时功率谱熵较低。当颤振开始出现，频谱分布逐渐复杂，频率分量变得分散，导致功率谱熵增高。随着颤振严重程度的增加，频谱分布开始混乱，频率分量更加分散，功率谱熵进一步增大。由图8可知，在1.5秒前后，归一化的PSE值小幅上升，到第4 s左右，PSE值又逐渐下降，再从第5秒开始，PSE值又急剧增大，最终趋于稳定。因此可知，第2组铣削实验的加工过程是在1.5 s左右时从稳定状态变为轻微颤振状态，从第4 s开始回归稳定铣削，然后再从第5 s开始逐渐过渡到严重颤振状态。此外根据铣削状态变化对应的时间与PSE值可知，归一化PSE值低于0.3可以认为加工过程处于稳定铣削状态，当PSE值处于0.3~0.6之间为轻微颤振，高于0.6时为严重颤振。

4.3. 不同参数下的算法验证

为验证所提算法的通用性，选取剩余两组铣削加工信号进行算法验证。通过SSA分别寻得第1组与第3组信号的最佳分解参数(k, α)为(5, 1670)和(3, 1092)，再进行分解重构，并提取归一化PSE值，结果如图10所示。

由图10(a)可知，通过第1组铣削信号的时域波形可以判断，该组实验经历了稳定加工–轻微颤振–严重颤振–稳定加工–轻微颤振–严重颤振的加工变化过程，同时PSE值在1.5 s时小幅上升进入0.3~0.6区间，在第6s时急剧增大超过0.6，随后逐渐下降到0.3以下，再在9 s时回归0.3~0.6区间，最终在12.5 s后又超过0.6，PSE对应的加工状态符合判断。同理，第3组铣削实验经历了稳定加工–轻微颤振–稳定加工–轻微颤振–严重颤振的加工过程，PSE值也符合加工状态。因此，两组实验的铣削状态变化不仅与加工信号的时域波形变化对应，同时也与不同状态的PSE值区间相对应，由此验证了所提方法的通用性。

(a) 第1组实验结果

(b) 第3组实验结果

Figure 10. The original signal, reconstructed signal, and normalized power spectral entropy

图10. 原始信号、重构信号和归一化功率谱熵

5. 结论

由于颤振现象对机床铣削系统的影响巨大，而传统的时频分析难以实时地识别加工状态。因此，提出一种基于SSA-EE-VMD与PSE结合的颤振识别方法。本文首先进行轴向变切深铣削实验以采集加工振动信号，通过对加速度信号的时频分析得出颤振频率的能量分布特性。然后采用SSA基于EE进行寻优，解决VMD的分解参数选择难题。再利用能量比筛选IMFs重构信号以消除噪声干扰。最后，提取PSE作为颤振特征指标识别加工状态，并进行不同参数下的算法验证。

研究结果表明，通过时频分析得出，频率分布随加工状态的改变，由主轴旋转频率向颤振频率转移。采用SSA-EE找到的最佳分解参数，能够使VMD清晰地分离颤振频率，避免模态混叠。此外PSE可以从频率分量角度作为颤振特征指标识别不同加工状态。由此验证了本文所提算法的有效性，并为铣削加工系统的颤振监测技术提供新的途径。