基于自適應動態無偏最小二乘支持向量機的刀具磨損預測建模

肖鵬飛 張超勇 羅 敏 林文文

1.華中科技大學機械科學與工程學院，武漢，4300742.寧波大學機械工程與力學學院，寧波，315211

0 引言

刀具在切削過程中會逐漸磨損退化直至失效，刀具在初期磨損、正常磨損到嚴重磨損直至失效過程中將不同程度地降低產品質量，甚至會損壞機床，影響加工系統運行，造成重大損失[1]。統計研究表明:采用刀具監測技術可以有效縮短停機時間，顯著提高生產率和機床利用率；實時監測刀具狀態不但延長了刀具自身的使用時間，而且還能有效防范刀具失效造成的工件報廢和機床故障，進而削減成本[2]。

根據測量手段的不同，刀具狀態監測方法主要分為直接法和間接法。直接法較常見的有接觸測量法和機器視覺圖像測量法[3],直接法因監測過程的高成本和不協調而在加工領域的應用受到限制。間接法通過利用刀具磨損間接相關的信號進行監測，可在不中斷加工過程的情況下實現在線監測，非常適合自動加工系統，受到越來越多的關注。切削加工進程中的許多特征都可用于判別刀具的狀態，如刀具和主軸振動所發出的聲音和振幅信號，切削力、切削溫度、主軸功率或扭矩的實時變化，以及聲發射信號等。在間接監測方法中，刀具磨損監測系統的組成主要包括研究對象(刀具)、加工條件、傳感器檢測、信號處理、特征信號提取與選擇以及磨損量預測等模塊和環節。

張臣等[4]研究表明，與前向神經網絡方法相比，采用最小二乘支持向量機(least square support vector machine，LSSVM)構建監測模型所得到的監測結果更加準確。VIJAYAKUMAR[5]在標準支持向量機結構風險的置信范圍中添加了b2/λ2(λ>0)項，提高學習速度的同時，對泛化能力的影響不大。由于訓練樣本數量有限，以及滑動時間窗長度和監測模型不能自適應調整和更新等因素，傳統基于機器學習的刀具磨損預測模型存在精度和效率較低等問題，本文在滿足預測精度要求的基礎上引入了滑動時間窗長度自適應確定算法，提高了數據集長度的自適應性，構建的動態無偏最小二乘支持向量機(dynamic non-bias least square support vector machine，DNLSSVM)在設計上消除偏置項，增加訓練樣本動態增減過程，有效利用核擴展矩陣對稱正定性質，簡化動態學習過程中Lagrange乘子的求解，從而提高了建模效率。

1 自適應動態無偏最小二乘支持向量機

1.1 建立初始預測模型

首先利用滑動時間窗內的初始訓練數據集建立初始預測模型[6]：設長度為l，則初始訓練數據集的形式為{(xi,yi)}(i=1,2,…,l)，其中輸入數據xi∈Rn，輸出數據yi∈R。在最小二乘支持向量機結構風險的置信范圍中添加b2/(2λ2)(λ>0)項，該模型的目標函數和約束條件為

(1)

i=1,2,…,l

令ω′=(ω;b/λ)，則式(1)變成

(2)

i=1,2,…,l

式中,φ(·)為非線性映射函數，它將數據集從輸入空間映射到特征空間；b為線性擬合模型的偏置項，參與決定超平面到原點之間的距離；ξ為擬合誤差變量；C為針對擬合誤差的懲罰因子，C越大，則目標函數對擬合誤差的容忍度越??；λ為為了消除偏置b而引入的參數，取較大值會更近似原模型，但值太大又會影響學習過程中的收斂速度和穩定性。

為快速求解上述問題，采用Lagrange乘子法，建立Lagrange函數并利用KKT條件把約束優化問題轉變成函數優化問題，解得

(3)

式中，αi為Lagrange乘子。

對于i=1,2,…,l，消去ω′、ξi，并將式(3)寫成矩陣形式，即

(K+λ2E+C-1I)α=Y

(4)

Kij=(φ(xi)·φ(xj))=k(xi,xj)K=[Kij]

Y=(y1,y2,…,yl)Tα=(α1,α2,…,αl)T

式中，E為l階全一方陣；I為l階單位方陣。

初始預測模型的形式為

(5)

由式(5)可以看出，通過引入參數λ消除了回歸函數中的偏置項。

令H=K+λ2E+C-1I(λ>0,C>0)，式(4)簡化為Hα=Y，稱H為核擴展矩陣，可以證明其為對稱正定。H為對稱正定矩陣，能夠進行Cholesky分解，可以被唯一地分解為H=UTU，其中U為上三角陣。由UTUα=Y，設P=Uα，則UTP=Y。初始預測模型的Lagrange乘子向量α的求取公式為

(6)

(7)

式中，Pi、αi分別為P和α的第i個分量;uik為矩陣U的第i行、第k列元素。

如果通過求H的逆矩陣進而求解Lagrange乘子，則H維度很大時計算復雜，而利用上述公式可以簡化計算過程。

1.2 動態更新模型

設t時刻滑動時間窗中有l個樣本[7]，則該時刻訓練數據集的形式為{(xi,yi)}(i=t+1,t+2,…,t+l)。t+1時刻進來一個新樣本(xt+l+1,yt+l+1)，同時淘汰時間窗中該時刻最早的舊樣本(xt+1,yt+1)。

在線訓練算法中，核函數矩陣K、訓練樣本輸出Y和待求的Lagrange 乘子α都是t的函數，表示如下：

Y(t)=(yt+1,yt+2,…,yt+l)T

(8)

α(t)=(αt+1,αt+2,…,αt+l)T

(9)

Kij(t)=k(xi,xj)

(10)

(11)

H(t)為動態正定矩陣，能夠進行Cholesky分解。由

K(t)=

(12)

相應地

H(t)=

(13)

由t時刻的學習結果，令H(t)=U(t)TU(t)。t+1時刻增加新樣本(xt+l+1,yt+l+1)后，對應H(t)有

(14)

V(t+1)=[k(xt+l+1,xt+1)+λ2

k(xt+l+1,xt+2)+λ2…k(xt+l+1,xt+l)+λ2]T

v(t+1)=k(xt+l+1,xt+l+1)+λ2+C-1

(15)

1.3 自適應計算訓練數據集長度算法

一般情況下，用于儲存訓練數據的滑動時間窗過短可能會導致信息不全和模型精度不高，過長又會導致過擬合和建模速度慢，因此本文采用自適應的時間窗長度選取算法，綜合構建自適應動態無偏最小二乘支持向量機(adaptive dynamic non-bias least square support vector machine，ADNLSSVM)。

損失函數Qt-1的計算公式如下:

(16)

令Qt-2=Qt-1/l，則損失函數相對下降量

(17)

滑動時間窗的調整需要預先設定預測誤差閾值θ和目標函數相對下降量閾值ε，在初始時間窗中，不斷添加新的樣本建立動態模型，并預測下一樣本輸出值，在模型預測誤差小于閾值θ的情況下，如果連續迭代m次目標函數相對下降量Δt-1均不超過閾值ε，則滿足迭代終止條件[6]，如圖1所示。

圖1 滑動時間窗自適應選取算法流程圖Fig.1 Algorithm of adaptively selecting the lengthof sliding time window

2 數據分析與處理

2.1 數據分析

本文采用美國紐約預測與健康管理學會(PHM)2010年高速數控機床刀具健康預測競賽開放數據中的銑削實驗數據[8]進行試驗，試驗的系統結構和主要設備如圖2所示。銑削試驗用到的主要設備和試驗條件見表1。

圖2 試驗系統結構和設備Fig.2 Structure and devices of the experimental system

記每一次走刀過程為108 mm長度的端面銑，由于端面銑加工材料為正方形，刀具每一次走刀時間間隔相等，因此可以用走刀數來代替切削時間，每次走刀后采用LEICA MZ12顯微鏡測量刀具實際的后刀面磨損量作為刀具磨損結果，一次走刀為一個監測過程。

表1 銑削試驗主要設備和試驗條件Tab.1 Main equipments and machining conditions ofmilling experiment

試驗釆用6把獨立的銑刀(C1號銑刀～C6號銑刀)進行刀具全壽命試驗。下面以C4組銑削試驗作為研究對象，該組數據記錄了315次走刀試驗的振動信號和對應三列磨損數據值，分別表示球頭銑刀三個切削刃的磨損量，本文采用三個切削刃磨損量的均值表示銑刀磨損量大小。圖3繪出了刀具三個切削刃磨損均值的變化曲線。

圖3 刀具磨損曲線Fig.3 The curve of tool wear amount

由圖3可以發現，1～19次走刀曲線較陡，磨損量上升較快后趨于平緩，該階段內銑刀處于初期磨損階段；20～176次走刀銑刀磨損量上升過程相對較慢，處于正常磨損階段；177次走刀以后刀具磨損速度急劇增加，進入急劇磨損階段，磨損量超過閾值便磨鈍失效了[9]。

2.2 銑削振動信號的小波降噪

首先對信號進行零均值化和去除趨勢項的預處理。

依據分類情況，為了增強信號的對比度，在三類磨損狀態中分別選取第1、第158和第315次的三向振動信號進行分析，其對應的磨損量分別為24.2 μm、95.3 μm、203.1 μm。對這三次走刀過程中各方向(X、Y、Z)銑削振動信號進行頻譜分析，按照小波包分解理論[10]，濾除高頻噪聲，其中,X方向降噪前后頻域圖見圖4。

(a)降噪前

(b)降噪后圖4 X方向降噪前后分類頻域圖Fig.4 The frequency-domain magnitude before andafter denoising in X direction

其他方向銑削信號的降噪處理采用相同的方法完成。

2.3 特征提取

運用時域分析、頻域分析及小波時頻域分析方法對切削過程的切削振動信號進行分析處理，把振動信號中能夠直接有效反映刀具態變化的監測特征抽取并選擇出來。

2.3.1時域特征提取

(1)時域有量綱特征。有量綱特征參數以數值形式直觀反映振動信號的時域統計特性，計算簡單快捷，便于在線實時處理，對刀具磨損和破損等嚴重切削故障較為敏感，是對振動信號最簡單的時域描述，如均值、標準差、均方根值等統計量。

(2)時域量綱一特征。時域有量綱特征是統計量，需要較多的歷史數據，而且受加工材料、切削參數等工況條件的影響，沒有統一的衡量標準；量綱一特征能在一定程度上克服有量綱特征的這些缺點；兩個具有相同量綱的量的比值組成量綱一特征參數，這些參數對信號的幅值和頻率變化不敏感，對故障缺陷敏感[8]，如峰值因子、峭度指標、偏度指標等。

圖5繪出了銑削振動信號(X方向)的時域有量綱和量綱一參數與刀具磨損量間的變化關系。由圖5可知：振動信號的均值在零附近波動，逐漸接近于零，直觀上與磨損量無相關性(相關系數為0.129 0)；銑削振動信號的標準差與均方根特征相關性好(相關系數分別為0.922 4、0.921 9)，隨著磨損量的增加而遞增；量綱一參數中，偏度指標與磨損量的相關性好(相關系數為0.716 2)，其他兩個指標與刀具磨損量基本不相關(相關系數分別為0.344 5、0.441 1)。

(a)有量綱

(b)量綱一圖5 時域有量綱和量綱一參數變化圖Fig.5 The variation curves of features of dimension and non-dimension in time-domain

2.3.2頻域特征提取

圖6 頻域特征參數變化圖Fig.6 The variation curves of features in frequency-domain

信號的功率譜能很好地揭示刀具運行過程中切削振動頻率成分變化的動態信息，被廣泛應用于刀具磨損狀態監測。常用的譜特征獲取方法有相關圖法功率譜估計、周期圖法功率譜估計。頻域特征參數有頻段能量、重心頻率、頻率方差和均方頻率等。本文釆用周期圖功率譜估計法對銑削監測信號進行功率譜估計，然后計算其頻域特征參數。圖6所示為X向切削振動信號的頻域特征參數隨刀具磨損量變化情況。由圖6可知：銑削振動頻率特征參數隨刀具磨損量的增大總體呈上升趨勢，但具有波動性。

2.3.3時頻域提取

本節分析和提取小波包頻帶能量與小波熵特征。先對釆樣信號進行頻譜分析以確定小波包分解的層次。分析得知：不同磨損狀態下信號頻譜規律性較強，頻譜圖中各峰值頻率呈倍頻關系；進行小波包分解時，頻帶分辨率應不大于518 Hz。由采樣頻率50 kHz，可以確定進行6層小波包分解。小波包分解層次可以依據以下公式求得[10]：

fmin=Fs/2n+1

(18)

式中，Fs為采樣頻率;n為層數;fmin為最小頻段。

由于銑削振動頻率主要集中在1～10 kHz范圍內，故應取前32頻帶(1～12.51kHz)能量作為監測特征。

圖7列舉了銑削X方向振動信號小波包分解的1～4頻帶能量。由圖7可知：銑削振動信號的部分頻帶能量特征相關性好，小波包頻帶能量熵能反映刀具磨損狀態的變化。銑削振動小波能譜熵與磨損量的相關性比較好。

2.4 特征選擇

特征與磨損量之間的相關程度不盡相同并且特征之間也可能有某種關系，這會造成不必要的計算負擔和冗余。本文通過相關性分析即分析特征和磨損量間的關系緊密程度，從而識別并挑選出與磨損量相關程度更高的特征作為模型輸入，提高建模效果。相關性研究方法包括灰色關聯度分析法和相關系數法。

2.4.1相關系數法

用相關系數ρxy衡量兩個數據序列X和Y之間的線性相關性，計算公式如下：

(a)部分頻帶能量

(b)小波能譜熵圖7 部分頻帶能量和小波能譜熵變化圖Fig.7 The variation curves of 4-band energy and the wavelet energy entropy

(19)

經過計算，得到時域、頻域上各個特征量與磨損量間的相關系數，如表2所示。

根據表2，以|ρxy|>0.9為標準，可以篩選出一部分特征量，在X、Y、Z方向上滿足條件的分別有4、2、2個特征量，故總共有8個特征量。

經過計算，X、Y、Z各個方向振動信號小波包分解后時頻域上前32個頻帶能量和能量熵與磨損量間的部分相關系數,如表3所示。

根據表3，以|ρxy|>0.9為標準，可以篩選出一部分特征量，在X、Y、Z方向上時頻域滿足條件的分別有15、13、12個特征量，故總共40個特征量。

匯總時域、頻域和時頻域上的特征量，總共48個特征量，對這些特征量進行特征編號，如表4所示。

2.4.2灰色關聯度分析

灰色信息是指外延可以確定，但是內涵不確定的、隱含而未公開的信息。兩個系統中的子系統或因素間存在隨時間或不同對象而變化的關聯性，它們之間的相似或相異性即灰色關聯度，被用作衡量它們之間的關聯程度。

分析之前，首先進行數據標準化處理，即量綱一化。設參考數列X*為磨損量序列，其他數列Xi(i=1,2,…)為特征量序列，計算Xi與X*的關聯系數ζi(t)。令αi(t)=|Xi(t)-X*(t)|，t表示序列中元素順序，t∈{1,2,…,N},則有

式中，ρ∈(0,1),一般取值0.5。

表2 時域和頻域特征量相關系數表Tab.2 The correlation coefficients of features in time and frequency domain

表3 時頻域32個頻帶能量和能量熵的部分相關系數Tab.3 Parts of correlation coefficients of 32-band energy and the wavelet energy entropy in time domain

表4 特征匯總編號表Tab.4 The sequencing number of features

由關聯系數得到Xi與X*的關聯度

(20)

將48個編號后的特征量進行灰色關聯度分析。圖8為特征量關聯度降序排列圖。由圖8可知：以關聯度0.65為分界，可以淘汰掉10個特征(編號10～13，16，29，30，32，41，43)，則總共還剩余38個特征量。

圖8 灰色關聯度分析排序圖Fig.8 Grey relational grade values and their sorting sequence

3 建模與仿真

在64位計算機MATLAB 2015b平臺上編輯程序，選用高斯徑向基核函數，用粒子群優化算法得到組合參數懲罰因子C=810.4、核函數參數σ2=1.2和引入參數λ=26.8。以第四組數據中的前60個數據作為訓練數據，后225個數據作為測試數據，分別采用傳統在線LSSVM和DNLSSVM進行225次在線預測，比較在不同滑動時間窗長度條件下，兩種算法的運行時間。等間隔選取10個不同時間窗長度，對應運行時間如圖9所示。

圖9 算法的測試運行時間Fig.9 The test running time of DNLSSVM algorithom

由圖9分析得出：隨著滑動時間窗長度的增大，兩種算法的測試運行時間都增長，DNLSSVM算法運行時間增長相對平穩和緩慢；一般情況下滑動時間窗長度越長，DNLSSVM的效率優勢越明顯。

為了平衡訓練數據量和建模時間，利用滑動時間窗長度自適應調整算法，可以求取滑動時間窗長度，目標函數相對下降量隨滑動時間窗長度變化如圖10所示。由圖10可知，在預測精度滿足要求θ≤1.5 μm的條件下，選取閾值ε=0.03，可以確定滑動時間窗長度為23。

圖10 目標函數相對下降量隨滑動時間窗長度變化圖Fig.10 The variation curve of the relative decrementof object function

輸出預測仿真擬合效果,如圖11所示。預測過程性能用平均相對誤差r和平均絕對誤差a衡量，5次試驗后兩種算法的預測誤差見表5，計算公式如下：

由圖9、圖11以及表5可知，ADNLSSVM算法建模具有較好的建模效率和擬合精度。

4 結論與展望

(1)ADNLSSVM算法通過在設計上消除偏置項，增加訓練樣本動態增減過程，能夠有效利用核擴展矩陣對稱正定性質以及歷史數據，提高模型的適應性，簡化動態學習過程中Lagrange乘子的求解，因此該方法和傳統在線LSSVM建模方法相比具有更高的建模效率，并且在滑動時間窗長度越長時優勢表現越明顯。

(2)自適應窗口長度確定算法兼顧了時間窗長度和擬合精度要求。與在線LSSVM算法相比，在不影響模型泛化能力的情況下，具有滿足要求的擬合精度。用平均相對誤差和平均絕對誤差作為衡量模型精度的評價指標，試驗結果表明建立的模型具有相對更好的擬合精度。

(a)DNLSSVM算法

(b)LSSVM算法圖11 模型擬合效果圖Fig.11 Fitting performance of the prediction model

實驗次數ADNLSSVM算法傳統在線LSSVM算法a(μm)r(%)a(μm)r(%)11.16630.931.24220.9821.11140.891.11190.8931.02120.881.12700.9041.10160.871.11630.8951.11190.891.13690.91

(3)模型動態更新過程中需要利用上一階段的實際磨損量，可以結合一些較為成熟的無接觸直接測量刀具磨損量的技術比如機器視覺檢測等技術進行狀態在線跟蹤和預測。

參考文獻：

[1] 莊子杰.基于聲發射和振動法的刀具磨損狀態檢測研究[D]. 上海：上海交通大學，2009.

ZHUANG Zijie. Tool Wear Condition Monition Based on AE and Vibration[D]. Shanghai：Shanghai Jiao Tong University,2009.

[2] TONSBOFF H K. Developments and Trends in Monitoring and Control of Machining Processes[J].Annual of the CIRP,1988,37(2):611-622.

[3] 楊吟飛,李亮,何寧.一種新的刀具磨損面圖像邊界提取方法[J].南京航空航天大學學報,2007,39(6):786-789.

YANG Yinfei, LI Liang, HE Ning. Extraction Method for Image Boundary Based on Tool Wear[J]. Journal of Nanjing University of Aeronautics and Astronautics,2007,39(6):786-789.

[4] 張臣,周來水,安魯陵,等.球頭銑刀刀具磨損建模與誤差補償[J]. 機械工程學報,2008,44(2):207-212.

ZHANG Chen, ZHOU Laishui, AN Luling, et al.Modeling and Wear-induced Error Compensation of Ball-end Milling Cutter Wear[J]. Chinese Journal of Mechanical Engineering,2008,44(2):207-212.

[5] VIJAYAKUMAR S, WU S. Sequential Support Vector Classifiers and Regression[C]// International Conference on Soft Computing. Genoa,1999:610-619.

[6] 蔡艷寧, 汪洪橋, 葉雪梅.復雜系統支持向量機建模與故障預報[M].北京：國防工業出版社,2015:4.CAI Yanning, WANG Hongqiao, YE Xuemei. Support Vecto Machine Modeling and Fault Prediction for Complex System[M]. Beijing: National Defence Industry Press,2015:4.

[7] 張浩然,汪曉東.回歸最小二乘支持向量機的增量和在線式學習算法[J].計算機學報,2006,29(3):394-406.

ZHANG Haoran, WANG Xiaodong. Incremental and Online Learning Algorithm for Regression Least Squares Support Vector Machine[J]. Chinese Journal of Computers,2006,29(3):394-406.

[8] 李威霖，傅攀，張爾卿.基于粒子群優化LS-SVM車刀磨損量識別技術研究[J].計算機應用研究,2014,31(4):1094-1097.

LI Weilin, FU Pan, ZHANG Erqing. Application of Particle Swarm Optimization Least Square Support Machine in Tool Wear Monitoring[J]. Application Research of Computers,2014,31(4):1094-1097.

[9] 謝厚正,黃民.基于振動測試的數控機床刀具磨損監測方法[J].儀表技術與傳感器,2013(2):73-76.

XIE Houzheng, HUANG Min. Research of Numerical Control Machine Tools Wear Monitoring Method Based on Vibration Test[J]. Instrument Technique and Sensor, 2013(2):73-76.

[10] 叢龍慧,韓玉杰.小波分析在刀具磨損狀態檢測中的應用[J].林業機械與土木工程,2009,37(4):51-52.

CONG Longhui, HAN Yulin. Application of Wavelet Analysis in Tool Wear Status Detection. Forestry Machinery & Woodworking Equipment, 2009, 37(4):51-52.