一种基于变分模态分解和能量熵的铣削颤振识别方法

摘要

本发明公开了一种基于变分模态分解和能量熵的铣削颤振识别方法，属于机床加工颤振识别技术领域。其包括以下步骤：S1、建立VMD数学模型；S2、建立能量熵的数学模型；S3、进行分别代表稳定切削、微弱颤振及严重颤振三种切削状态的三组铣槽加工实验，通过测力仪来获取三组铣削力信号；S4、对三组铣削力信号进行FFT分析，证明三组铣削力信号分别代表加工处于稳定切削、微弱颤振及严重颤振状态；S5、通过基于峭度值的VMD参数自动选取方法确定VMD分解的最佳模态个数K和惩罚因子α；S6、求多个IMF的瞬时频率，确定铣削颤振特征频带；S7、采用锤击实验获取刀具模态；S8、基于能量熵对每个IMF提取颤振特征向量。本发明提高了VMD分解效果，实现了自动识别颤振。

说明书

技术领域

本发明属于机床加工颤振识别技术领域，涉及一种基于变分模态分解和能量熵的铣削颤振识别方法。

背景技术

为了提高材料切除率和减小切削力，高速铣削加工被广泛地应用在航空航天工业中。机床颤振源于切屑形成中的自激机理，刀具-工件系统的某个模态最初被切削力激励，在接近但不等于加工系统主结构频率处产生颤振。颤振会降低表面质量、生产效率，也会造成刀具的磨损。对于加工工艺系统的颤振问题，很多学者已经提出颤振稳定性预测、识别、抑制等方法，但对于机床操作者来说，颤振预测的方法很难得到有效实施，因此，为了实现零件的高效、高精加工，研究颤振识别技术变得重要。

智能加工和智能机床是智能制造的重要内容，研究铣削过程中的在线颤振识别技术也对实现智能制造具有重要意义。颤振特征提取是在线颤振识别的关键技术，而实现颤振的特征提取前提是准确确定颤振频带。早期颤振识别允许操作者干预加工过程和避免颤振对工艺系统的破坏。近年来，针对颤振识别的研究已经广泛的展开。

变分模态分解(VMD)是2014年由Dragomiretskiy等提出一种新的非递归式信号处理方法。VMD一经提出，便在不同的领域得到了应用。有文献指出，VMD方法更能准确提取特征。Wang利用VMD、经验小波变换(EWT)、集合平均模态分解(EEMD)、经验模态分解(EMD)等四种方法对频谱信号进行提取，结果表明VMD在信号分析的优势。VMD方法在特征提取中得到了广泛的应用，但是很少文献将VMD方法用于铣削过程中颤振特征的提取上。尽管VMD的理论基础扎实，分解效果较好，但VMD对信号进行分解时，模态个数K和惩罚因子的选择通常依据经验进行选择，VMD分解的效果不佳。

铣削过程中，伴随着能量的转移。稳定切削阶段，能量被消耗在主频及其谐频上，当切削状况发生恶化时，随着颤振频率幅值的增大，能量被消耗在颤振频带上。熵指的是体系的混乱的程度，能量熵是熵在能量域的延伸。当铣削过程发生颤振时，能量熵会增大。所以，基于能量熵随频率的变化来表明切削阶段是否发生颤振是可行和有效的方法。但目前通常选取整个信号的能量熵作为颤振特征，没有真正意义上提取到颤振频带和颤振特征向量。

综上所述，亟需一种自动选取VMD参数，基于能量熵来提取颤振特征向量的铣削颤振识别方法。

发明内容

(一)要解决的技术问题

为了解决现有技术的上述问题，本发明提供一种自动选取VMD参数，基于能量熵来提取颤振特征向量的铣削颤振识别方法。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

一种基于变分模态分解和能量熵的铣削颤振识别方法，包括以下步骤：

S1、建立VMD数学模型；

S2、建立能量熵的数学模型；

S3、进行分别代表稳定切削、微弱颤振及严重颤振三种切削状态的三组铣槽加工实验，通过测力仪来获取三组铣削力信号；

S4、对三组铣削力信号进行FFT分析，证明三组铣削力信号分别代表加工处于稳定切削、微弱颤振及严重颤振状态；

S5、通过基于峭度值的VMD参数自动选取方法确定VMD分解的最佳模态个数K和惩罚因子α；

S6、求解多个IMF的瞬时频率并确定铣削颤振特征频带；

S7、采用锤击实验来获取刀具模态；

S8、采用基于能量熵的方法对每个IMF提取颤振特征向量。

根据本发明，所述步骤S1包括以下步骤：

S11、定义IMF为一个调幅调频信号；

S12、构造变分问题；

S13、得到每个IMF的解析信号，并得到其单边频谱；

S14、将每个IMF的频谱调制到相应的基频带；

S15、通过平方范数估计出各IMF的带宽；

S16、求取约束变分模型的最优解；

S17、采用ADMM解决步骤S12的变分问题；

S18、解得中心频率的更新方法。

根据本发明，所述步骤S2包括以下步骤：

S21、求得每个IMF的能量；

S22、基于香农熵定义VMD的IMF的能量熵。

根据本发明，所述步骤S5包括以下步骤：

S51、引入峭度；

S52、将峭度引入铣削颤振；

S53、选定惩罚因子α及模态个数K的范围及步长；

S54、在每种K和α的组合下，利用VMD对铣削力信号进行多次分解，将每次VMD获得的n个IMF相加，得到该次VMD后的重组信号；

S55、计算多个重组信号的峭度值，并保存最大的峭度值和记录此时的惩罚因子α和模态个数K。

根据本发明，所述步骤S8具体为：

S81、结合FFT分析和瞬时频率，对各个IMF进行判断，判断各个IMF是否为颤振频带，若是，选最大能量熵的IMF作为颤振特征向量，否则此IMF不能作为颤振特征向量。

根据本发明，所述步骤S54中的重组信号个数为380。

根据本发明，所述步骤S55中的最大的峭度值为4744，取最大的峭度值时的模态个数K为13，惩罚因子α为4200。

(三)有益效果

本发明可以基于峭度值自动选取VMD参数，解决了现有采用VMD分解时获得的模态个数K和惩罚因子α需凭经验设定，导致VMD分解效果不佳的问题，大大提高了VMD分解的效果；在最佳模态个数K和惩罚因子α的参数组合下对信号进行VMD分解得到各本征模态函数(IMF)，当铣削过程中发生颤振时，能量向颤振频带转移，为准确确定颤振频带，采用基于能量熵对每个IMF提取颤振特征向量的方法，真正意义上提取到了颤振特征向量，不仅理论基础扎实，而且可以实现自动识别颤振，准确有效地提取颤振特征，为铣削颤振特征提取提供了新的手段，对实现智能加工有重要意义。

附图说明

图1为本发明具体实施方式中基于VMD和能量熵的铣削颤振识别方法的流程图；

图2为三组时域和频域切削力分析图，其中，(a)第1组切削参数下；(b)第2组切削参数下；(c)第3组切削参数下；

图3为不同K和α组合下的峭度三维图；

图4为第三组信号VMD后的各IMF的瞬时频率；

图5为三组信号在K＝13、α＝13号在时的各IMF的能量及能量熵值。

具体实施方式

为了更好的解释本发明，以便于理解，下面通过具体实施方式，对本发明作详细描述。

具体实施方式

一种基于VMD和能量熵的铣削颤振识别方法，如图1所示，包括以下步骤：

S1：建立VMD数学模型。

本实施方式中，VMD基于变分问题来求解本征模态函数(IMF)。IMF是具有一定宽度的频带。变分问题是求泛函的极值问题。为解决变分问题，利用乘法算子交替方向法(ADMM)不断更新各IMF及其中心频率，限制条件是其带宽之和为最小。然后将解得的IMF解调到相应的基频带，最终提取各个IMF及相应的中心频率。VMD具有两个重要的分解参数：模态个数K和惩罚因子α。VMD算法的目标是通过构造并求解约束变分问题，将原始信号分解为K个IMF分量。

S11、定义IMF为一个调幅调频信号。

本实施方式中，定义IMF为一个调幅调频信号，其表达式为：

其中，A_k(t)为u_k(t)的瞬时幅值，w_k(t)为r_k(t)的瞬时频率：

其中，A_k(t)及w_k(t)相对于相位来说是缓变的，即在[t-δ，t+δ]的间隔范围内u_k(t)可以看作是一个幅值为A_k(t)、频率为w_k(t)的谐波信号。

S12、构造变分问题。

本实施方式中，VMD算法中将变分问题描述为：在约束条件为各IMF之和等于输入信号f下，寻求k个IMF使得每个IMF的估计带宽之和最小，将每个IMF的频谱调制到相应的基频带。

S13、得到每个IMF的解析信号，并得到其单边频谱。

本实施方式中，通过Hilbert变换，得到每个IMF的解析信号，并得到其单边频谱：

S14、将每个IMF的频谱调制到相应的基频带。

本实施方式中，通过加入指数项调整各IMF估计的中心频率，将每个IMF的频谱调制到相应的基频带：

S15、通过平方范数估计出各IMF的带宽。

本实施方式中，通过平方范数来估计出各IMF的带宽，如果将原始信号y(t)分解为K个IMF分量，则对应的约束变分模型表达式如下：

其中，{u_k}＝{u₁,....,u_k}代表VMD分解得到的K个IMF分量；{w_k}＝{w₁,....,w_k}表示各IMF分量的频率中心。

S16、求取约束变分模型的最优解。

本实施方式中，为求取约束变分模型的最优解，VMD通过引入二次惩罚因子α和Lagrange算子λ(t)，将待求解的约束性变分问题转变为非约束性变分问题，其中二次惩罚因子α可保证信号的重构精度。

S17、采用乘法算子交替方向法(ADMM)解决步骤S12的变分问题。

本实施方式中，采用ADMM解决步骤S12的变分问题，通过交替更新λⁿ⁺¹寻求扩展Lagrange表达式的解。得到待求解的二次优化问题解为：

S18、解得中心频率的更新方法。

本实施方式中，解得中心频率的更新方法：

其中，相当于当前剩余量的维纳滤波，为当前IMF功率谱的中心。

S2、建立能量熵的数学模型。

本实施方式中，熵(entropy)指的是体系的混乱的程度，能量熵是熵在能量域的延伸。当铣削过程发生颤振时，能量熵会增大。所以，基于能量熵和频率的变化来表明切削阶段是否发生颤振是可行和有效的方法。

S21、求得每个IMF的能量。

本实施方式中，设VMD分解后的各IMF分别为u₁(t)、u₂(t)、…u_N(t)，代表从低到高的频带，则每个IMF的能量如下：

其中，R_i分别代表VMD的IMF的能量，由于各IMF的正交性，故原信号的能量可以用各IMF的能量之和表示。

S22、基于香农熵定义VMD的IMF的能量熵。

本实施方式中，基于香农熵定义VMD的IMF的能量熵为：

T_i＝R_i/R表示每个VMD分解后的IMF的能量占整个信号能量的百分比。

其中，R＝{R₁、R₂、…R_N}构成了VMD分解部分的能量在频域上的分布。

S3、进行分别代表稳定切削、微弱颤振及严重颤振三种切削状态的三组铣槽加工实验，通过测力仪来获取三组铣削力信号。

本实施方式中，进行三组铣槽加工实验，通过测力仪来获取铣削过程中铣削力信号。三组铣槽加工实验分别代表三种不同切削状态：稳定切削、轻微颤振、严重颤振。

本实施方式中，采用沈阳中捷机床厂的四轴数控机床TH5650，铣刀直径D＝10mm，齿数N＝2，铝合金方块尺寸为100mm×50mm×50mm。采用切削力信号来识别颤振，测力仪Kistler 9257B采集x、y、z三向力信号，整个切削过程是在干铣削条件下进行。测力仪用固定板固定在工作台上，工件通过虎钳夹持。设置测力仪采样频率fs＝2000Hz。其加工参数如表1所示：

表1 切削参数

S4、为说明三组力信号的频率组成，对三组铣削力信号进行快速傅里叶(FFT)分析，证明三组铣削力信号分别代表加工处于稳定切削、微弱颤振及严重颤振状态。

本实施方式中，三组力信号及FFT分析如图2所示：第1组切削参数下，力的幅值在区间[-150，150]上，第2组切削参数下，力的幅值在区间[-250，250]上，第3组切削参数下，力的幅值在区间[-250，400]上，可知随轴向切深的增加，时域切削力的幅值增大。

本实施方式中，实验中采用n＝1000r/min，故主频f＝n/60＝16.67Hz，而刀齿通过频率f_tooth＝2×f＝33.34Hz，频域图中用等间隔的浅细实线来表示主频及其倍频。

本实施方式中，频谱图中可以看出，当在第1组切削参数下即轴向切深a_p＝0.5mm时，系统较大幅值的频率主要集中在刀齿通过频率及其倍频上，表明此切削过程处于稳定切削；当在第2组切削参数下即轴向切深a_p＝0.7mm时，系统较大幅值的频率还是出现在主频的倍频处，但也出现了其他较大幅值的频率，图中黑色圆圈放大后显示，较大幅值的频率点没有和浅细实线重合，说明系统的能量向这些频率点转移，系统出现微弱颤振；当在第3组切削参数下即轴向切深a_p＝1mm时，系统较大幅值的频率不再出现在主频的倍频处，图中黑色圆圈放大后显示，较大幅值的频率点没有和浅细实线重合，说明系统的能量主要集中在这些频率点，系统出现严重颤振。

S5、通过基于峭度值的VMD参数自动选取方法确定VMD分解的最佳模态个数K和惩罚因子α。

本实施方式中，当铣削过程中发生颤振时，信号中大幅值的概率密度增加，信号幅值的分布偏离正态分布，峭度值增大，故把峭度值作为VMD分解参数选取的指标。

S51、引入峭度。

本实施方式中，峭度(Kurtosis)是反映信号分布特性的数值统计量，是归一化的4阶中心矩，对于离散数据，峭度定义为：

其中，x_i为信号值，为信号均值，N为采样长度，σ_i为标准差。

S52、将峭度引入铣削颤振。

本实施方式中，铣削过程中，颤振的发生伴随着信号幅值的增大。峭度值越大，说明信号幅值增长较大，表明在铣削过程中发生了颤振。

S53、选定惩罚因子α及模态个数K的范围及步长。

本实施方式中，选定惩罚因子α及模态个数K的范围及步长，设[S_α，L_α]和T_α分别表示α的范围和步长，[S_K，L_K]和T_K分别表示K的范围和步长，则影响参数α和K共有种组合形式。

本实施方式中，对第3组力信号采用VMD分解，选定惩罚因子α的[S_α，L_α]和T_α分别为[200，2000]和100，模态个数K的[S_K，L_K]和T_K分别为[1，20]和1。

S54、在每种K和α的组合下，利用VMD对铣削力信号进行多次分解，将每次VMD获得的n个IMF相加，得到该次VMD后的重组信号。

本实施方式中，共有19x20＝380种组合，对信号进行380次VMD分解。VMD分解时，如果模态个数K＝n(其中n＝1，....20)，则可以获得n个IMF。如果把n个IMF相加，得到的和就是VMD分解后的重组信号。

S55、计算380个重组信号的峭度值，并保存最大的峭度值和记录此时的惩罚因子α和模态个数K。

本实施方式中，对380个重组信号计算峭度值，通过MATLAB画出380个K和α组合下的峭度值，如图3所示，从图3中可以看出，最大峭度为4744，此时模态个数K＝13，惩罚因子α＝4200。峭度值最大的重组信号说明信号幅值的分布严重偏离正态分布，最能体现颤振的影响，故此组合(K，α)＝(13，4200)为最佳分解组合。

S6、在最佳模态个数K＝13和惩罚因子α＝4200组合下对力信号进行VMD分解，求解13个IMF的瞬时频率并确定铣削颤振特征频带。

本实施方式中，因为模态个数K＝13，所以共获得13个IMF。

本实施方式中，瞬时频率不仅能定位各频率发生的时间，也能表明各IMF的频率组成。为表明各IMF的频率组成，对每个IMF求瞬时频率。

本实施方式中，瞬时频率的求解过程如下：

设解析信号为：

则瞬时幅度为：

A(t)＝±|x(t)| (12)

则瞬时相位为：

其中，对瞬时相位求一阶导数为瞬时频率：

本实施方式中，当IMF的瞬时频率中如果出现了FFT(图2)中的颤振频率点，则此IMF是颤振频带。从图4中可以看出，u1和u2是包含2f和4f的频带，u3-u7属于主频的谐振频带，原信号频谱中颤振频率主要集中在[600-800Hz]，而u8-u12的频带范围在[600-800Hz]，可知u8-u12是颤振特征频带，包含颤振频率。此方法虽能确定颤振特征频带，但不利于自动提取颤振特征向量。

S7、为进一步确定颤振的发生，依据颤振理论，采用锤击实验来获取刀具模态。

本实施方式中，通过FFT分析来表明铣削过程中是否发生颤振，缺乏理论依据。颤振理论指出：机床颤振源于切屑形成中的自激机理，刀具-工件系统的某个模态最初被切削力激励，在接近但不等于加工系统主结构频率处产生颤振。因此采用锤击实验为获得加工系统主结构频率。如果力频谱中存在接近但不等于加工系统主结构频率的频率，说明铣削过程发生了颤振。

本实施方式中，刀具总长75mm，刀具悬长45mm，冲击力锤型号为B&K8206-002，灵敏度为0.23mV/N，最大冲击力22N，谐振频率722kHz，加速度传感器型号为877a500，灵敏度为10.00mV/g，用Cutpro动态信号采集系统来处理信号，USB信号采集器型号为T/O-9233，设定采集频率范围是50-5000Hz。

本实施方式中，将加速度传感器先后固定在刀具和工件上进行锤击实验来获取刀具和工件的传递函数。当固定在刀具上时，传感器的位置要尽量靠近刀尖部分。

本实施方式中，通过力锤模态实验结果中如表2所示：

表2 刀具模态参数

刀具参数固有频率阻尼比刚度一阶模态695.50.05467.9451×10⁶二阶模态9260.02484.6542×10⁶

本实施方式中，刀具的第一模态为695.5Hz，而第二组和第三组力信号的频谱图(FFT)中，存在着较大幅值的695.5Hz频率。依据颤振理论，说明第二组和第三组力信号中包含颤振信息。

S8、为实现自动提取颤振特征向量，采用基于能量熵的方法对每个IMF提取颤振特征向量。

本实施方式中，针对结合瞬时频率和FFT方法不利于自动提取颤振特征向量，提出基于能量熵的颤振特征向量提取方法。铣削过程处于稳定切削时，能量被消耗在机床主结构的主频及其谐频处；当有颤振发生时，能量会向颤振频带转移。

本实施方式中，从图5中可以看出，能量及能量熵的变化趋势是一致的。当铣削过程处于稳定切削时，能量及能量熵的峰值出现在u2上，而u2对应着刀齿通过频率所在频带；当轻微颤振时，能量及能量熵的峰值出现在u10上，但u1、u2也具有较大的能量及能量熵，说明能量已从主频及其倍频处转移到颤振频带；当严重颤振时，能量及能量熵的前5个峰值出现在u8、u9、u10、u11、u12上，而u1、u2有较小的能量及能量熵，说明能量值主要集中在颤振频带。

本实施方式中，u8、u9、u10、u11、u12都可以作为颤振特征向量。为减少颤振特征向量的数量和实现自动提取颤振特征向量，选取所有能量熵中具有最高峰值的u10作为颤振特征向量，其在三组切削实验下u10的能量熵值如表3所示：

表3 三组切削实验下u10的能量熵

切削状态稳定切削(实验1)轻微颤振(实验2)严重颤振(实验3)能量熵0.05490.3550.359

但需要注意的是：在颤振产生初期，能量主要还是集中在主频及其倍频处，所以最大能量熵的IMF不一定是颤振频带，需要排除IMF不属于主频带及其谐频带。

本实施方式中，当铣削实验处于稳定切削时，u10的能量熵较小，仅为0.0549；而当铣削中发生轻微颤振时，u10的能量熵增加到0.355；当铣削中发生严重颤振时，u10的能量熵略有增加。结果表明，提出的基于VMD和能量熵结合的颤振特征提取方法能准确有效提取颤振特征，为铣削颤振特征提取提供了新的手段。

以上结合具体实施方式描述了本发明的技术原理。这些描述只是为了解释本发明的原理，而不能以任何方式解释为对本发明保护范围的限制。基于此处的解释，本领域技术人员不需要付出创造性劳动即可联想到本发明的其它具体实施方式，这些方式都将落入本发明保护范围之内。