基于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)型高階模糊時間序列顫振預測模型.通過該模型，以完整的顫振發展過程的信號特征量時間序列作為訓練集，可得到測試集的預測值，在保證良好預測性能的基礎上實現了較長時間段的預測.