基于贝叶斯理论的刀具剩余寿命预测

Prediction of Remaining Tool Life Based on Bayesian Theory

摘要: 刀具健康状态的监测以及刀具剩余寿命的预测，对于企业与厂家的降本增效来说异常重要。针对单一通道预测所带来的数据预测不精确，包含的退化信息过少，可靠性低等问题，提出一种基于隐马尔可夫模型及多通道融合的刀具剩余寿命预测方法。首先通过计算多通道信号所提取特征的时间序列与对应时间向量的斯皮尔曼等级相关系数对特征时序做单调性排序，其次筛选出单调性优异的特征代入隐马尔可夫模型进行训练，并构建得到健康因子作为观测数据，最后借助贝叶斯理论和马尔科夫链蒙特卡洛采样估计退化模型，实现实时在线不断更新的退化模型模拟真实的退化进程。本文用PHM2010公开数据挑战赛中三槽球头硬质合金铣刀切削不锈钢过程的磨损全寿命数据验证了方法的有效性。

Abstract: The monitoring of tool health status and the prediction of remaining tool life are exceptionally im-portant for cost reduction and efficiency of companies and manufacturers. To address the problems of inaccurate data prediction, low reliability due to the inclusion of too little degradation infor-mation caused by single channel prediction, a tool remaining life prediction method based on Hid-den Markov Model and multi-channel fusion is proposed. Firstly, the time series of the extracted features of the multi-channel signals are ranked monotonically by calculating the Spearman rank correlation coefficients of the corresponding time vectors, secondly, the features with excellent monotonicity are selected and substituted into the hidden Markov model for training, and the health factors are constructed as the observed data, and finally the degradation model is estimated with the help of Bayesian theory and Markov Chain Monte Carlo sampling to achieve real-time online The continuously updated degradation model simulates the real degradation process. This paper validates the validity of the method using wear full life data from the PHM2010 Open Data Challenge for a three-flute ball-tipped carbide milling cutter cutting stainless steel.

1. 引言

刀具作为直接与工件接触的部分，其性能的好坏直接影响工件的优劣。而在机械加工的过程中不可避免地存在刀具磨损的现象。严重的刀具磨损对加工效率和加工质量均会产生较大影响，同时也关系到加工成本，在精密超精密切削过程中，刀具的微量磨损会引起零件的精度误差 [1] [2] 。要想获得一把刀具的磨损量，从而达到刀具剩余寿命预测的目的 [3] ，一般都是对刀具进行大量实验，进而得以分析该刀具的剩余使用寿命。因此，本文对刀具所进行的实时剩余寿命预测对保障生产安全，企业大幅降本增效有着重大的意义。

随着机械加工向集成化，机械化，复杂化方向的发展，传统的刀具状态监测技术无法满足生产加工中复杂多变的需求。近年来，随着智能传感技术的大力发展，人工智能的广泛运用，研究人员借助诸如机器学习，人工运维等方法对工件的剩余寿命进行了高效的预测。目前，常用的方法诸如：支持向量机(support vector machine, SVM)并结合多种优化算法对模型或者流程改进的方法。Celestino [4] 等人使用自动回归综合移动平均模型(ARIMA)来预先估计预测变量的值，并将其作为支持向量回归模型(SVM)的输入，新方法远远超过了VARMA模型所提供的预测能力。梅飞等人 [5] 基于融合粒子群的模糊聚类与支持向量机结合，实现设备在故障样本数据量较小情况下建模，解决断路器故障分类诊断问题。除了支持向量机的方法，Yang等人 [6] 结合长短时记忆(LSTM)法与递归神经网络实现基于大数据的旋转机械设备故障智能分类。Li和Liu [7] 提出了一个改进的HMM，构建了一个危险模型来描述时间变化和条件适应的状态转换概率，并使用正向算法在线预测RUL。Akhilesh Kumar [8] 等人建立了基于隐藏马尔科夫模型和多项式回归的刀具自主诊断和预测模型。但是，上述成果大多基于单通道信号进行研究，对于本文的刀具剩余寿命预测而言，采用单一信号进行预测易受外界环境干扰，同时基于单一传感器的剩余寿命预测方法存在数据利用率低，可靠性较低的问题。因此，多通道的信息融合和巧妙应用就显得尤为重要。对于多通道问题，李乃鹏等 [9] 提出一种基于多传感器数据融合模型的剩余寿命预测方法，该算法首先根据传感器在RUL预测中的单独性能进行优先排序，然后根据其综合性能选择一个最佳的传感器组。

对刀具而言常用监测信号有切削力、振动、声发射、电流等。从原始监测信号中提取的特征有些可能与刀具磨损退化的相关性较小，另一方面，受外部环境、工况波动等随机因素的影响，信号与刀具退化的相关性会出现波动。而对多通道信号提取的特征进行融合能够减小对某单一信号的依赖，即使某监测信号出现异常也能通过其他信号构建稳定的退化指示量，避免了依赖单一信号而造成的预测失稳，大大提高预测模型的鲁棒性。综上，本文提出一种基于隐马尔可夫模型及多通道融合的刀具剩余寿命预测的方法，能够规避不利因素，实时迭代退化模型逼近刀具实际磨损退化趋势，同时对每个时刻的剩余寿命进行实时预测。

本文首先通过计算多通道信号所提取特征的时间序列与对应时间向量的斯皮尔曼等级相关系数对特征时序做单调性排序，其次筛选出单调性优异的特征代入隐马尔可夫模型进行训练，并构建得到健康因子作为观测数据，最后借助贝叶斯理论和马尔科夫链蒙特卡洛采样估计退化模型，实现实时在线不断更新的退化模型模拟真实的退化进程。

2. 信息融合及健康指标构建

2.1. 多通道特征提取和筛选

刀具剩余寿命整体预测框架如图1所示。首先对每个通道提取包括std，peak2peak，RMS，PV，AP，RP，在内的六种时域特征。根据Spearman等级相关系数对所提特征进行排序，过滤掉一些单调性分数高于预定阈值的特征。

Figure 1. Framework diagram of the proposed forecasting methodology

图1. 预测方法框架图

${\sigma }_s=\sqrt{\frac{1}{n-1}}{\displaystyle {\sum }_{i=1}^n{\left(s_i-{\mu }_s\right)}^2}$ (1)

$F_{pp}=\max s-\min \left(\text{s}\right)$ (2)

$F_{rms}=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{s}_i^2}}$ (3)

$\text{PV}=\max {\displaystyle {\sum }_{i=1}^n\left|s_i\right|}$ (4)

$\text{AP}=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n\left|s_i\right|}$ (5)

$\text{RP}={\left(\frac{1}{n}{\displaystyle {\sum }_{i=1}^n\sqrt{\left|s\right|}}\right)}^2$ (6)

2.2. 基于多链马尔科夫过程的多通道融合健康指标构建

在剩余寿命预测问题中，健康因子是用来描述系统退化趋势的定量指标，可分为物理量 [10] (直接监测数据)和虚拟量(间接监测数据)，具体取决于监测数据是直接还是间接 [11] 。例如，在电池健康状态评估和剩余寿命预测中 [12] ，锂离子电池的剩余使用寿命(RUL)是由放电表面温度预测的，便于实时测量和在线估算。在机械设备的剩余寿命预测中，由于空间可及性和退化机制的隐蔽性，很难直接观察到机械设备运行状态的变化，通常只能通过监测振动、力和声发射等信号间接反映设备的退化情况。

由于受到外部环境中随机因素的影响，单通道监测信号会出现异常值波动，这样构筑的健康因子会导致预测结果严重不稳定。在本文中，我们提出了一种特征选择和多通道融合的组合，通过将选定的特征代入隐马尔可夫模型，并将其作为贝叶斯-MCMC退化模型的参数估计的输入观测值来构建健康因子。

2.3. 信号特征的提取

经过前期调研，对于长度为n的信号 $s=\left[s_1,s_2,\cdots ,s_n\right]$ ，本文在提取特征时选用常见的六个时域特征：标准差，峰峰值，有效值，峰值，平均幅值，均方根。对N个通道在T个时间单位中的时序数据，分别提取表中所列特征的时间序列，形成特征矩阵X^T*M，其中M = 6*N。用Fi{1:T}表示X^T*M的第i列特征时序。

2.4. 特征筛选

因为刀具退化过程随时间是严格单调的，因此需要筛选出单调性好的特征以表现出退化趋势。斯皮尔曼等级相关系数用于评估两个变量之间的关系程度，因此需要计算每个特征时序与对应的时间矢量的斯皮尔曼等级系数作为该种特征的单调性评分。对于 $F_{\left\{1:T\right\}}^i$ ，其计算公式为

$S\left(F_{\left\{1:T\right\}}^i\right)=1-\frac{6\times {\displaystyle {\sum }_{k=1}^T{r}_k^2}}{T\left(T^2-1\right)}$ (7)

式中r_k为矢量R中元素，

$R=rank\left(F_{\left\{1:T\right\}}^i\right)-rank\left(t_{\left\{1:T\right\}}\right)$ (8)

式中， $rank\left(F_{\left\{1:T\right\}}^i\right)$ 为 $F_{\left\{1:T\right\}}^i$ 的排序。

对特征时序按上述方法进行筛选。由于各通道信号所提取的特征数据的量级有很大差别，因此需要进行标准化处理，将不同的数据转到相同的量级上。

2.5. 构造隐马尔可夫模型

经过对多通道所提取的特征信号的筛查分析可知，被测刀具在前10次切削中完全处于正常磨损状态，因此选取前10次的刀具切削数据构造隐马尔可夫模型，利用隐马尔可夫模型对极少的刀具正常磨损数据进行分析和建模，并将生成的结果构造健康因子，便于代入后续的寿命预测与分析。

通过隐马尔可夫模型构建的健康因子包含了过多的杂波，为了增加预测的稳定性，对输出的健康因子进行三次一阶指数平滑处理，如式所示，其中，a为平滑系数， $y_t^{\left(1\right)}{y}_t^2{y}_t^{\left(3\right)}$ 分别为第一二三次的平滑值。经过三次平滑后的值所构成的时间序列，即 $Y_{\{1}{}_{:T\}}=\left[y_1,\cdots ,y_t,\cdots ,y_T\right]$ 为最终构建的健康因子。

$\{\begin{array}{l}{y}_t^{\left(1\right)}=\alpha x_t+\left(1-\alpha \right)y_{t-1}^1,y_0^{\left(1\right)}=x^0\hfill \\ y_t^{\left(2\right)}=\alpha y_t^{\left(1\right)}+\left(1-\alpha \right)y_{t-1}^2,y_0^{\left(2\right)}=y_0^{\left(1\right)}\hfill \\ y_t^{\left(1\right)}=\alpha y_t^{\left(2\right)}+\left(1-\alpha \right)y_{t-1},y_0=y_0^{\left(2\right)}\hfill \end{array}$ (9)

3. 基于贝叶斯-MCMC的参数后验联合概率密度计算及采样方法的基本原理

通过前期对刀具特性的考察与研究，本文采用指数退化模型z = −1 + a * exp (b * t)描述刀具退化过程。指数退化模型能良好刻画系统退化趋势，被包括轴承、刀具等在内的机械设备RUL预测广泛采用 [13] [14] 。指数模型中由两个未知参数θ = [a, b]控制退化轨迹的趋势，准确估计a和b的值即为退化模型精确与否的关键。而贝叶斯推断则能够不断融合获取的时序观测数据即健康因子，实时更新θ的联合后验概率分布并从中采集ns个样本，将采集的样本代入退化模型中沿时间外推直至超过预先设定的阈值，得到剩余寿命的概率分布。而贝叶斯后验概率密度通常没有封闭解析解，因此结合MCMC采样算法，基于贝叶斯后验概率密度数值计算公式从后验概率密度中获得待估参数的样本。

在初始时刻待估参数θ = [a, b]中a和b的先验概率分布分别设定为服从参数为μ_a,σ_a的对数正态和参数为μ_b,σ_b的正态分布。初始时刻假设a,b相互独立，则其联合概率密度的先验分布f_pr如式(10)所示。参数μ_aσ_bμ_aσ_b值的指定遵循以下规律：如果对指定μ_a,μ_b的值有较高置信度，则σ_a,σ_b应取较小的值，反之，σ_a,σ_b应指定较大的值以保证对a,b的估计主要取决于观测数据而非a,b的先验分布。

$f_{pr}\left(\theta |{\mu }_a,{\sigma }_a,{\mu }_b,{\sigma }_b\right)=\frac{1}{a{\sigma }_a\sqrt{2\pi }}\exp \left(\frac{-\left(\ln a-{\mu }_a\right)}{2{\sigma }_a}\right)\times \frac{1}{{\sigma }_b\sqrt{2\pi }}\exp \left(\frac{-{\left(b-{\mu }_b\right)}^2}{2{\sigma }_b}\right)$ (10)

在基于贝叶斯理论的预测理论背景中，似然函数的意义为：在t时刻，当待估参数θ取其可能取值 ${\theta }^i=\left[a^i,b^i\right]$ 时，理论退化值 $z_t^i\left(a_t^i,b_t^i\right)$ 与当前观测数据y_t的相似程度，记为 $f_{li}\left(y_t|z_t^i\left({\theta }^i\right)\right)$

$f_{li}\left(y_t|z_t^i\left({\theta }^i\right)\right)=\frac{1}{{\sigma }_s\sqrt{2\pi }}\exp \left(-\frac{{\left(y_t-z_t^i\left(a_t^i,b_t^i\right)\right)}^2}{2{\sigma }_s{}^2}\right)$ (11)

可以看出，似然函数同时是时间和θ的函数：对于固定时刻t它是θ的函数，而当θ取其可能取值 ${\theta }^i=\left[a^i,b^i\right]$ 时则是时间的函数。在获取到T个时序观测数据 $Y_{\{1:T\}}$ 后，θ的后验联合概率分布为给定 $Y_{\{1:T\}}$ 时θ的条件概率，记为 $f_{po}\left(\theta |Y_{\{1:T\}}\right)$ ，根据贝叶斯理论

$f_{po}\left(\theta |Y_{\left\{1:T\right\}}\right)=\frac{f_{pr}\left(\theta \right)\times f_{li}\left(Y_{\left\{1:T\right\}}|z_T\left(\theta \right)\right)}{{{\displaystyle \int }}^{\text{}}{f}_{pr}\left(\theta \right)\times f_{li}\left(Y_{\left\{1:T\right\}}|z_T\left(\theta \right)\right)\text{d}\theta }$ (12)

$F_{li}\left(Y_{\{1:T\}}|z_T\left(\theta \right)\right)$ 为在获取T个观测数据时的总体似然函数，为每一时刻似然函数的连乘。基于贝叶斯-MCMC的参数后验联合概率密度计算方法总结为初始化参数 ${\theta }^0=\left[a^0,b^0\right]$ ，令 ${\theta }^{old}={\theta }^0$ ，代入计算 $f_{po}\left({\theta }^{old}|Y_{\{1:T\}}\right)$ ，for $i=1:n_s$ ，

$f_{po}\left({\theta }^i|Y_{\left\{1:T\right\}}\right)\propto \frac{1}{a^i{\sigma }_a\sqrt{2\pi }}\exp \left(\frac{-\left(\ln a^i-{\mu }_a\right)}{2{\sigma }_a}\right)\times \frac{1}{{\sigma }_b\sqrt{2\pi }}\exp \left(\frac{-{\left(b^i-{\mu }_b\right)}^2}{2{\sigma }_b}\right)\times {\displaystyle {\prod }_{j=1}^T\exp \left(-\frac{{\left(y_j-z_j^i\left({\theta }^i\right)\right)}^2}{2{\sigma }^2}\right)}$ (13)

从建议的对称分布中采样得到参数 ${\theta }^{new}$ 。将 ${\theta }^{new}$ 代入式(13)中，计算 $f_{po}\left({\theta }^{new}|Y_{\{1:T\}}\right)$ 。 if $u<\min \left\{1,f_{po}\left({\theta }^{new}|Y_{\{1:T\}}\right)/f_{po}\left({\theta }^{old}|Y_{\{1:T\}}\right)\right\}$ ，则接受 ${\theta }^{new}$ ，并将 ${\theta }^{new}$ 与 $f_{po}\left({\theta }^{new}|Y_{\{1:T\}}\right)$ 的值赋给 ${\theta }^{old}$ 与 $f_{po}\left({\theta }^{old}|Y_{\{1:T\}}\right)$ ，进入下一循环，Else，拒绝 ${\theta }^{new}$ ，接受 ${\theta }^{old}$ ，最终获得T时刻n_s个参数 $\theta$ 的估计值 $\hat{\theta }$ 。在T时刻刀具RUL的分布由时间外推形成的退化轨迹刻画。

4. 实验验证

4.1. 实验数据介绍

本文所提方法采用PHM2010数据挑战赛的公开数据进行验证 [15] [16] [17] [18] ，该组数据可以清晰展现刀具退化趋势，便于研究人员进行数据分析，被多国学者引用 [16] [17] [18] [19] 。试验环境如图2所示，在高速数控机床Roders Tech RFM760上采用三槽球头硬质合金铣刀进行铣削试验，被切削工件材料

Figure 2. Experimental environment and raw data presentation

图2. 实验环境和原始数据展示

为硬度HRC52不锈钢，主轴转速10,400 r/min，进给率1555 mm/min，径向切深0.2 mm，轴向切深0.125 mm。在切削过程中，以50 kHz采样速率采集切削力信号(三个方向)、加速度信号(三个方向)及声发射信号共7个通道的信号，每完成一次工件表面铣削，用莱卡MZ12型号显微镜离线测量刀具磨损量，直至现场专业人员认为磨损量超出正常范围时停止试验，共完成315次切削。共使用6把刀具进行了6组切削试验。本案例以一次切削过程所采集的信号为一个时间单位进行单一通道特征提取，共315个时间单位。以4号(#4)刀为例，图2展示了切削力和振动信号的全寿命数据，可以看出，随切削次数增加导致的磨损加剧，切削力和加速度信号随时间有明显变化。

4.2. 特征提取及优选

由于在前期仅获取少量数据的情况下估计退化模型的参数十分困难，所估计的参数不确定性大，对剩余寿命预测意义甚微，另一方面，在给定切削参数及切削材料的情况下，刀具寿命有大致范围，且刀具前期磨损退化较稳定，在早期对刀具开展剩余寿命预测工程意义不大。因此，本案例从第150次切削开始预测刀具剩余寿命。前150次切削所采集的原始数据作为训练数据为特征筛选及构建健康因子提供基础和依据。仍以第四把刀为例，针对每一次切削，每个通道构建六个特征，共构建42个特征的时间序列，形成特征矩阵X^{150 × 42}。

用上文所述的斯皮尔曼等级相关系数计算每个特征时序的单调性得分。以#4刀为例，结果如图3所示，共选出x方向切削力信号的标准差，峰峰值，有效值，峰值，平均幅值，均方根；y方向切削力信号的标准差，峰峰值，有效值，峰值，平均幅值，均方根；z方向切削力信号的标准差，峰峰值，有效值，峰值，平均幅值，均方根；x方向加速度信号的峰峰值和峰值；z方向加速度信号的标准差，峰峰值，有效值，峰值，平均幅值，均方根在内的26种特征。如图4所示，每个通道都提取包括std，peak2peak，RMS，PV，AP，RP，在内的六种时域特征，因此对于不同刀具，依据上述方法筛选出的特征种类和数量存在差异，且能够自适应地筛选出单调性好的特征。本例中，#3，#5，#6刀具分别提取出36，31，36种特征。

Figure 3. Monotonicity feature for all channels of #4 tool extraction feature timing

图3. #4刀具在所有通道下提取的特征图

Figure 4. Feature maps extracted under all channels

图4. 在所有通道下提取的特征图

Figure 5. Health indicator of #4 tool

图5. #4号刀具的健康因子

4.3. 健康指标构建

对筛选出的特征代入隐马尔可夫模型进行训练，并得到所构建的健康因子，如图5所示。同时，健康因子将作为贝叶斯-MCMC计算框架的观测值对指数退化模型参数实时更新，并在每时刻对刀具剩余寿命迭代预测。贝叶斯-MCMC计算框架参数设定如下：μ_a = 1, σ_a = 100, μ_b = 1, σ_b =100, a₀ = 4, b₀ = 0.1。MCMC采样中的建议分布为均匀分布，采样宽度范围为[−0.5, 0.5]。待估参数a,b初始时刻的先验概率密度函数中，由于对均值μ_a = 1，μ_b = 1置信度较低，因此指定了较大的方差σ_a = 100,σ_b = 100，在这种情况下a, b的估计值就主要取决于观测数据而非人为指定的先验知识，在对a, b没有很好先验知识的情况下降低人为设定的先验分布对参数估计的影响。

4.4. 刀具剩余寿命预测

以四号刀为例，其全寿命数据的验证结果如图6所示，图6中给出了RUL真实值、预测值及95%置信区间曲线。

Figure 6. #4 Tool 150 cuts prediction chart

图6. #4号刀具150次切削预测图

4.5. 结果分析与总结

从图6可以发现，在切削后期RUL预测值有些偏离RUL真实值，但一直处于“欠预测”状态，即保守估计了刀具可用剩余切削次数，该情况只会导致提前换刀，不会对产品加工精度造成影响，因此，相比“过预测”(可用剩余切削次数的预测值大于实际值)，“欠预测”是真实场景中较为理想的情况。

引入两个常用预测算法评价指标：均方根误差(Root mean square error, RMSE)、平均绝对误差(Mean absolute error, MAE)来评价预测结果。

RMSE与MAE计算如式所示，其中er_t为t时刻剩余寿命预测值与真实值之差。从第200次切削开始计算RMSE及MAE，结果如图7所示。

$RMSE=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{t=1}^{n}e{r}_t^2}}$ (14)

$\text{MAE}=\frac{1}{n}{\displaystyle {\sum }_{t=1}^n\left|e{r}_t\right|}$ (15)

预测值与实际值偏差较大的问题可通过延迟“首次预测时刻”ts以获取更多观测数据得到缓解，因为获取的观测数据越多，对退化模型参数的估计会越准确，当把ts推迟至200时，RUL预测曲线得到明显改善，如图7所示。RMSE值从15.8080降为14.0258，MSE值从12.6565降为11.9783。即使预测从第200次切削开始，也可以提前准确估计出全寿命的三分之一所能获得的剩余刀具寿命，从而为及时更换刀具提供理论依据，在保证工件加工精度和实现精确的刀具管理的同时，降低了成本，提高了效率。

Figure 7. #4 Tool 200 cuts prediction chart

图7. #4号刀具200次切削预测图

5. 结论

本文提出一种基于隐马尔可夫模型及多通道信号融合的刀具剩余寿命预测方法。能够从多通道信号提取的特征中自适应筛选出单调性高的特征，通过隐马尔可夫模型构建表征磨损退化趋势的健康因子。结合贝叶斯-MCMC退化模型参数后验概率分布计算及采样框架，对每个时刻的剩余寿命进行迭代估计。本文对多通道信息的融合利用可以精准把握刀具退化趋势并且基于贝叶斯理论对刀具退化过程进行实时监测，精准获取剩余寿命。整个过程避免了基于深度学习方法需要依赖大量全寿命数据来训练预测模型，且模型对新预测任务适应性有限的问题，简单易行，鲁棒性强，具有较好的工程实用性。