基于AR(p)型高阶模糊时间序列的磨削颤振预测方法

许小芬，陈换过，陈 培，许昆明

(浙江理工大学 浙江省机电产品可靠性技术研究重点实验室，浙江 杭州 310018)

磨削颤振是刀具和工件之间强烈作用而激发的一种振动，一般为再生型颤振.磨削颤振会降低零件加工表面的质量，加速刀具磨损，产生过大噪声，降低加工效率，加速机床损坏，严重时会导致加工的无法持续.对颤振进行预测，了解其变化规律，有利于及时采取控制颤振的相应措施，消除颤振对机械加工的不良影响.因此，针对磨削颤振的危害，寻求一种能及时有效预测磨削颤振的方法尤为重要.国内外许多学者曾对加工过程中颤振预测问题做过相关研究.Bediaga I等采用支持向量机的方法对高速铣削颤振进行了预测[1]；Gonzalez O等采用小波变换的方法对外圆磨削过程的颤振进行了预测[2]；Yao Z H等利用小波包分解方法构造颤振特征量，并用支持向量机对切削加工中颤振进行了预测[3]；Tansel I N等基于S变换获取振动信号特征值，采用模糊逻辑方法对车削加工中颤振进行了预测[4]；董新峰等提出一种基于经验模态分解复杂度和鉴别信息，对平面磨削颤振进行预测的方法[5].由于磨削颤振发生之际会产生幅值和频率均随时间变化的非平稳信号，且容易混入其他的干扰信号，因此上述大部分预测方法在具体应用时都会受到一定程度的限制.

基于模糊时间序列的预测方法在许多领域得到了应用.Chou H L等利用模糊时间序列预测模型对某一旅游景点人数进行了预测[6]；余文利等提出一种基于模糊C均值算法的模糊时间序列预测模型，对杭州市的气温进行了预测分析[7]；陈刚等针对模糊时间序列模型存在的论域划分问题，提出具有可调参数的模糊时间序列论域的非等分划分方法，来改进模糊聚类算法，对Alabama大学注册人数进行了预测[8]；蔺玉佩等利用模糊时间序列模型对股票指数走势进行了预测[9].但这些方法大都基于一阶模糊时间序列进行预测，预测结果不够理想.在模糊时间序列的实际建模过程中，只考虑相邻两个样本之间的联系，但忽略了不相邻样本之间可能存在的联系，使得预测结果不一定理想.本文采用AR(p)型高阶模糊时间序列的预测方法，考虑前期多个时刻的数据情况，提升样本数据之间的关联性，使得模型具有更强的预测性能.将AR(p)型高阶模糊时间序列的预测方法用于磨削颤振的预测，无需考虑磨削颤振信号的非平稳和非线性特点.首先，对磨削颤振信号进行小波降噪处理，提取降噪后信号的实时方差并将其作为反映颤振状态的特征量；然后引用本文方法对实时方差所构成的时间序列进行训练，选定AR(p)模型的阶数；最后，采取选定的模型进行拟合和预测.

1 AR(p)型高阶模糊时间序列

基于AR(自回归)变量的模糊时间序列的一般定义如下：

定义1：对任意一个固定的Y(t)(t=0,1,2,…)，设Y(t)⊂R，即为实数域的子集，在Y(t)上定义了一组模糊集fj(t)(j=1,2,…)，且F(t)={f1(t),f2(t),f3(t),…}，则称F(t)为定义在Y(t)上的模糊时间序列.

定义2：设R(t,t-1)为定义在F(t-1)到F(t)的模糊关系.对于t的任何值，如果R(t,t-1)不随t的变化而变化，则R(t,t-1)=R(t-1,t-2).此时F(t)叫作时不变模糊时间序列，否则叫作时变模糊时间序列.

定义3：若模糊时间序列F(t)只受F(t-1)的影响，则F(t-1)与F(t)之间的的模糊关系可表示为：

F(t-1)→F(t)

(1)

即被定义为一阶模糊AR(1)时间序列预测模型.在该模型中，时间序列仅由一个滞后时间序列决定.

定义4：若F(t)同时由多个滞后模糊时间序列F(t-1),F(t-2),…F(t-p)所确定，则F(t)与F(t-1),F(t-2),…F(t-p)之间的模糊关系可以表示为：

F(t-1),F(t-2),…F(t-p)→F(t)

(2)

式(2)即为第p阶模糊AR(p)时间序列预测模型.根据(2)式可以得出，在AR(p)模型中，时间序列受1阶、2阶、…p阶滞后时间序列的影响.

本文针对定义4的预测模型提出了具体算法.

步骤1:将论域划分为若干个模糊区间.

U为时间序列x(t)的论域，将论域划分为b个子区间，即U={u1,u2,…ub}.

论域U的数据范围由时间序列x(t)所决定.根据合适的区间长度将U划分为子区间ui.其中Xmax和Xmin分别为给定时间序列x(t)的最大值和最小值，给定两个任意值D1和D2，则论域U可表示为如下闭区间：

U=[Xmin-D1,Xmax+D2]

(3)

步骤2：定义模糊集，计算隶属度值.

模糊集是根据U和ui的时间序列来定义的，可表示为：

Ai=fAi(u1)/u1+fAi(u2)/u2+…+fAi(uq)/uq,fori=1,2,…,q

(4)

(5)

步骤3：将时间序列模糊化.

子区间ui是根据具体观测值确定的.最高隶属度对应的子区间确定了模糊集Ai.该观测值的模糊值为Ai.

步骤4：确定模糊逻辑关系，并建立AR(p)模型的模糊逻辑组关系表.

这一步骤中的模糊逻辑关系由式(2)定义的AR(p)模型确定.如：对于AR(2)模型中模糊逻辑关系，有：A2,A3→A3，A2,A3→A3,A2,A3→A5,则其模糊逻辑组关系就可以表示为A2,A3→A3,A3,A5.

步骤5：对时间序列进行模糊预测.

若滞后变量为：F(t-p)=Ai,F(t-(p-1))=Ar,…,F(t-1)=As且F(t)=Aj，则3种可能的情形如下：

情形1：若在模糊逻辑组关系表中存在Ai,Ar,…,As→Aj,Aj，则此时模糊预测值为Aj；

情形2：若在模糊组逻辑关系表中存在Ai,Ar,…,As→Aj,…Aj,Ak,…,Ak,A1,…,A1，且这里有a个Aj，b个Ak，c个A1，则此时的模糊预测值为Aj,…,Aj,Ak,…,Ak,A1,…,A1；

情形3：若Ai,Ar,…,As在模糊组逻辑关系表中并不存在对应的模糊规则，则此时的模糊预测值为Ai,Ar,…As.

步骤6：去模糊化预测.

将集中化方法用于去模糊化.当满足步骤5中情形1且模糊预测值为Aj时，去模糊化预测值为具有最高隶属度的某个模糊值Aj所对应的模糊区间uj的中间值.对于情形2和情形3，可通过计算加权平均数的方法来计算去模糊化预测值，即：

(6)

此时，具有最高隶属度的模糊值Aj,…,Aj,Ak,…,Ak,A1,…,A1所对应的模糊区间uj,uk,…,u1的中间值分别为mj,mk,…,m1.

2 磨削颤振预测的流程及实验方法

基于AR(p)型模糊时间序列预测模型的颤振预测流程如图1所示.

图1 颤振预测的流程

本文采用杭机股份有限公司KD4020X16型数控动梁龙门导轨磨床进行实验.该磨床的主轴转速为0～5 000 r/min，砂轮转速为1 400 r/min，进给速度为25 m/min，采用切入式磨削方式.磨削颤振实验平台主要由计算机、IEPE(Integral Electronic Piezoelectric)加速度传感器、CompactDAQ-9178动态信号采集系统等部分组成.

将8个IEPE加速度传感器分别设置在磨头、横梁和立柱上，具体测点位置分布在磨头的X、Y、Z向和电机立柱的X、Y、Z向以及横梁的X、Y向.

在实验过程中，通过动态信号采集系统采集8路振动加速度信号，并使用系统内的信号处理软件对信号进行预处理和降噪[10].对比振动比较强烈的测点信号后，选取从磨头X向测得的磨削颤振信号进行实验.采样频率为3 200 Hz.

3 磨削颤振特征量的提取

取一个包含完整颤振发展过程的振动信号(所取信号历时40 s，其中：0～8.5 s为平稳磨削阶段；8.5～16 s为颤振过渡阶段；16～40 s为磨削颤振阶段)，计算振动信号的实时方差.磨削颤振信号的时域如图2所示，磨削颤振信号的实时方差如图3所示.

图2 磨削颤振信号的时域

图3 磨削颤振信号的实时方差

由图3可看出：在平稳磨削阶段，磨削颤振信号的实时方差较小，且在接近于0的较小值附近变化；在颤振过渡阶段，实时方差急剧增大；到了磨削颤振阶段，实时方差在较大值附近平稳波动.显然，实时方差能够较好地反映磨削状态变化，故本文以它作为磨削颤振的特征量.

4 AR(p)型高阶模糊时间序列预测模型的建立

以一段完整的从磨削平稳状态到磨削颤振发生后的时间序列为元素，求取这段时间序列对应信号的实时方差，构成实时方差时间序列{X(t)}，共有128 600个数据点作为训练集，建立AR(p)型高阶模糊时间序列模型.分别取颤振过渡阶段磨削颤振信号实时方差时间序列的128个和690个数据点，作为两组测试集(分别为测试集1和测试集2)，为预测结果的验证做准备.

(1) 在应用程序中，用训练集进行训练，确定模型的参数，得到一组拟合值.首先，对时间序列{X(t)}进行论域划分，且每次训练均尝试不同长度的区间，从100个增加到2 000个，按每次增加50个模糊区间的方式进行训练；然后，采用三角形隶属度函数计算隶属度，确定每个原始数据的模糊集，从而，确定相邻数据之间的模糊逻辑关系(在这一算法步骤中，每次训练都需要尝试不同阶数的模型，本文采用p=2～5阶模型)；最后，根据算法的模糊规则，确定原始数据的拟合值.根据这些不同的长度区间和不同的模型阶数，按照本文提出的算法，重复步骤1到步骤6的计算过程，由此得到多组拟合值，并可计算每组拟合值的均方根误差(RMSE).在拟合值组别中，具有最小RMSE值的模型被确定为AR(p)型模糊时间序列的最佳模型.

(2) 在应用程序中，根据本文提出算法的步骤2到步骤4确定模糊逻辑关系.通过查找由训练过程得到的模糊逻辑关系表，对测试集进行去模糊化预测.计算测试集中平均绝对误差(MAE)、RMSE、平均百分比误差(MAPE)和决定系数(R2).通过这4个指标即可确定所提出方法的预测性能.

从p取值不同和论域划分方式不同的多种情形，对模型的参数进行反复调试，当p=2且选取1 000个模糊区间时，模型预测误差降到了最小，从而选取p=2且模糊区间个数为1 000的情形，相应地，AR(2)模糊时间序列被确定为预测模型.

5 预测结果及其分析

采用选取的模型，预测测试集1和测试集2中数据，并计算各误差指标，以检验所建立预测模型是否满足要求.其中，测试集1中真实值(又称原始值)取从颤振过渡阶段将要结束到颤振阶段约2.7 s内的原始方差；测试集2中原始值取颤振过渡阶段约2.2 s内的原始方差.所得测试集的预测值与原始值对比如图4所示.经计算，可得表1所示各误差指标结果.

图4 实时方差的预测值与真实值(即原始值)对比

测试集中数据/个MAEMAPERMSER21280.0 0180.0 1900.0 0240.9 9986900.0 0190.1 1530.0 0640.9 996

由图4(a)可看出，测试集实时方差的原始值和预测值总体变化趋势相同，可见预测效果良好.图4(a)和(b)曲线均表明，尽管测试集不同,但实时方差的预测值与原始值的总体偏差仍较小.

对于表1中4个误差指标来说，MAE、MAPE和RMSE的值越小，表明误差越小，模型效果越好；而R2的值越接近1，表明模型对数据拟合越好；越接近0，表明模型对数据拟合越差(一般而言，若R2>0.4，则认为模型拟合效果好).

在颤振过渡阶段，根据本文提出的预测模型，在当前时刻前期数据的基础上可预测当前时刻的实时方差.作为特征量的实时方差能很好地反映信号的变化情况.若实时方差在某时间段逐渐增大到一个相对平稳的较大值，则可以判定颤振发生，颤振发展过程已经处于颤振阶段.

6 结 论

(1) 根据磨床磨削过程中信号的特点，提取实时方差作为磨削颤振的特征量.在磨削过程中，实时方差随着磨削状态的变化而变化.在平稳磨削状态下，实时方差在一个接近于0的较小值附近波动；在磨削颤振过渡阶段，实时方差逐渐增大，呈现上升趋势；在磨削颤振阶段，实时方差在一个较大值附近平稳波动.在磨损加工中，可根据实时方差进行颤振状态的识别.

(2) 本文提出了一种基于AR(p)型高阶模糊时间序列颤振预测模型.通过该模型，以完整的颤振发展过程的信号特征量时间序列作为训练集，可得到测试集的预测值，在保证良好预测性能的基础上实现了较长时间段的预测.