基于HMM和優化的PF的數控轉臺精度衰退模型

王 剛， 陳 捷,2， 洪榮晶,2， 王 華,2

(1.南京工業大學 機械與動力工程學院，南京 211800；2. 江蘇省工業裝備數字制造及控制技術重點實驗室，南京 211800)

數控轉臺是高檔數控機床的重要功能部件，數控轉臺的分度精度和重復定位精度對數控機床的加工精度具有重要影響[1]。在數控機床中轉臺處于加工區域，工作環境惡劣，容易導致精度的過快衰減。已有的測量精度方法過于繁瑣，且在很多情況下限于環境因素難以測量。因此，找到一種以描述數控轉臺精度衰減為核心的建模方法不失為一種新的嘗試。目前，建模思路主要分為兩種[2]，一種是根據數控轉臺的物理結構建立數學模型；Yu等[3]分析了一種大型數控轉臺的電液耦合的動態特性，并且認為數控轉臺蝸輪傳動系統的綜合反向間隙是影響分度精度和動態特性的關鍵因素之一；然而，建立一個能準確且全面地反映數控轉臺的精度變化的數學模型十分困難；另一種是通過智能算法在振動信號等數據驅動的基礎上建立模型[4]，與數學模型相比智能算法不需要分解數控轉臺的結構，而是通過振動信號等數據作為觀察值，對蘊含在數據中的反應數控轉臺狀態變化的信息進行發掘，從而轉化為對數控轉臺狀態的評價。Benkedjouh等[5]使用小波變換進行降噪和提取特征向量，運用支持向量機(Support Vector Machine， SVM)建立剩余壽命預測模型。Lu等[6]在大型回轉支承的研究中使用小波能量熵(Wavelet Energy Entropy， WEE)和檢驗模態分解法(Empirical Mode Decomposition， EMD)進行降噪和特征向量提取，使用主成分分析法(Principal Component Analysis， PCA)對多維特征向量進行降維，最后運用最小二乘支持向量機建立了衰退預測模型。然而，數控轉臺精度變化的趨勢信息在轉臺的振動信號中往往非常微弱，很容易被周圍環境產生的白噪聲掩蓋。Matej等[7]在研究回轉支承時運用EEMD-MSPCA(Ensemble Empirical Mode Decomposition Multi-scale Principal Component Analysis)進行降噪并提取特征向量，為解決信號微弱問題提供了思路。在此基礎上Feng等[8]進行了改進，建立了EEMD-PCA算法，克服了其在降噪時經驗參數難以獲得的缺點。除此之外，數控轉臺的驅動形式、摩擦特性、系統剛度、阻尼因素、消隙預載等眾多因素都會影響到其最終的動靜態分度精度以及其精度保持性[9]，其中某些因素難以直接觀察，因此需要一些更具有針對性的算法。Wang等[10]在研究中將隱馬爾科夫(Hidden Markov Model, HMM)方法引入到齒輪箱的壽命預測中，比較有效地識別出某些難以觀察到的狀態。張西寧等[11]將主分量分析法(PCA)與隱馬爾科夫(HMM)結合，用到對軸承的診斷檢測中。

針對以上的問題，結合數控轉臺串聯系統的結構特點，利用隱馬爾科夫(HMM)算法既可以討論可觀察狀態又可以討論隱藏狀態的特點，提出用振動信號訓練過的隱馬爾科夫(HMM)模型進行精度衰退診斷并建立健康狀態指標，以粒子群算法(Particle Swarm Optimization, PSO)優化過的粒子濾波(Particle Filtering, PF)算法作為預測模型計算數控轉臺精度衰退的剩余壽命，并將優化前后的結果以及其他優化方法的結果進行對比分析。同時，在數控轉臺精度衰退加速壽命試驗過程中使用激光干涉儀定期測量轉臺的重復定位精度，以驗證該模型的有效性。

1 構建數控轉臺精度衰退預測模型

1.1 隱馬爾科夫(HMM)算法

隱馬爾科夫算法是建立在馬爾科夫假設、不動性假設和輸出獨立性假設的基礎上。其核心思想是通過可觀察到的狀態來推斷出隱藏狀態。對HMM有下面五個參數來表述[12]：馬爾科夫鏈的狀態數N；每種狀態對應的觀察值的數目M；狀態轉移概率矩A={αij}；在狀態j時的觀察概率分布B={βj(k)}；初始狀態概率分布v={vi}。HMM的表達式可以簡寫成λ=(A,B,v)。

在使用HMM建立觀察序列的模型時通常使用三種算法來解決三個基本的問題[13]：前向-后向算法解決評估問題，即已知觀察序列和概率模型，計算觀察序列的概率值P(O|λ)。Viterbi算法解決解碼問題，即已知觀察序列和概率模型，找到一個最可能的隱藏狀態序列S。Baum-Welch算法解決學習問題，即通過已知的觀察序列建立最有可能的HMM模型。

本文使用對精度變化敏感的若干統計特征量組成觀察矩陣OB作為觀察值序列，人為設置了四個隱藏狀態，訓練一個HMM模型。再使用訓練過的模型進行精度衰退的早期診斷。與以往使用HMM算法進行狀態識別不同的是，本文在進行早期診斷的過程中，使用HMM算法將數控轉臺的正常精度狀態與后面變化的精度狀態作比較，得出的是各個精度狀態下的精度與正常狀態下的精度的匹配程度，由此建立數控轉臺精度衰退的健康狀態指標H。

1.2 粒子濾波算法(PF)

(1)

在建立粒子濾波算法預測數控轉臺精度變化的過程中，首先要建立其狀態空間模型方程。在已有的健康指標H的基礎上，使用高斯函數對H進行非線性回歸擬合，得到擬合的狀態空間模型[15]

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

1.3 粒子群(PSO)算法優化粒子濾波算法

粒子濾波算法雖然在理論上可以表示任意形式的概率分布，但是依然存在一系列的問題[16]，上文中的三個狀態空間的參數ak,bk,ck在算法預測前就需要確定其初始值，現有的文獻中確定這個初始值時雖然可以利用粒子濾波的狀態跟蹤能力通過觀察數據跟蹤獲得，但獲得的最優參數只是一個范圍，具體值的確定依然需要依靠經驗。本文提出使用粒子群(PSO)算法找出這三個參數的最優值，達到優化粒子濾波算法的預測能力的目的。

PSO可以描述為：對d維空間中隨機初始化的一個粒子群的速度和位置進行迭代，每次迭代之后，計算其自身搜索到的歷史最優解pbest和全部粒子的歷史最優解gbest，并由適應度函數對每個粒子進行優化評價。迭代過程中每個粒子的速度和位置公式為

v=wi×vp+c1×r1×(pbest-x)+c2×r2×(gbest-x)

(10)

x=x+vp

(11)

式中：wi為慣性權重，本文取為0.7[17]；c1為粒子自身路徑中的歷史最值權重，一般取2；c2為群體最優值權重，一般取2；r1,r2為[0,1]內隨機數。

本文中的適應度函數為

(12)

式中：y為觀測值；yp為預測值；在經過300次迭代之后，粒子群算法停止優化。具體的優化過程如表1所示。

表1 PSO優化粒子濾波算法

2 轉臺精度衰退預測實現過程

轉臺精度衰退預測的具體實現過程如下：

(1) 信號降噪與重構，使用EEMD-PCA方法對原始數據降噪并進行信號重構[18]，通過使用聚合經驗模態法(EEMD)將原始振動信號分解為若干組IMF(Intrinsic Mode Function)分量，再對這些組的IMF分量進行PCA，從而篩選出其中包含低頻有效信息的IMF，再對篩選過的IMF分量重構成矩陣Ej，從而得到反映進度衰退的特征信號，同時起到降噪的目的。

(2) 選擇統計特征量，組成特征觀察值矩陣。通過Baum-Welch算法在已知的觀察序列的基礎上建立最有可能的HMM模型。本文通過觀察序列OB10建立一個HMM模型λ10，代入測試數據計算得到精度衰退的對數似然概率曲線并診斷是否發生衰退，在此基礎上建立健康狀態指標H。本文選擇了峭度(KU)、偏度(SK)、裕度指標(CL)、峰值指標(CF)、脈沖指標(IM)和波形指標(SI)等六個統計特征量。

由這些統計特征量組成的觀察值矩陣為

(13)

式中：OBj=[OBj1,OBj2, …,OBjK]的j對應下文的第j組數據。

定義描述轉臺精度衰退的健康狀態指標為

H=P(OBj|λ10)

(14)

式中：λ10為用健康狀態下的觀察值矩陣OB10訓練的HMM模型。

(3) 粒子濾波算法預測精度趨勢。在經過粒子群算法優化得到最優的起始點參數ak,bk,ck的值后，使用粒子濾波算法對精度衰退趨勢和剩余壽命進行預測，如表2所示。

表2 標準粒子濾波算法

其中，

(15)

訓練的粒子數定為2 000，σ1，σ2，σ3，σ4分別為8，0.007, 0.005，9。粒子濾波算法根據狀態空間模型對精度狀態進行迭代更新，輸出每一步的狀態估計值，在達到精度的閾值之后，根據精度的PDF分布以及預測失效點的位置計算數控轉臺的精度剩余壽命(Remaining Useful Life, RUL)。

轉臺精度衰退評估和剩余壽命計算的模型方案如圖1所示。

圖1 數控轉臺精度衰退模型方案Fig.1 Model of precision decline of NC rotary table

3 數控轉臺精度衰退加速壽命試驗

3.1 試驗設備簡介

本次的試驗對象為國內某機床附件公司生產的TK13250E型數控立臥轉臺。工作臺最高轉速為11.1 r/min，定位精度為30"，重復定位精度為8"(調整前)，經間隙補償后重復定位精度可達6"。如圖2(a)所示，該數控轉臺由西門子s7-200 224XP型PLC控制的步進電機驅動，該PLC(Programmable Logic Controller)同時控制一個加載液壓缸，振動加速度信號由PXI數據采集箱采集。加速度傳感器1安裝在蝸桿的軸線延長線位置上，加速度傳感器2安裝在電機傳動齒輪的鏈接法蘭上。由于蝸輪副的狀態變化對轉臺的定位精度具有顯著影響，所以本文中利用加速度傳感器1采集的精度衰退全壽命信號進行分析。本次試驗為轉臺的變載荷加速壽命試驗，在試驗的前200 h以100%的載荷運行轉臺，后520 h以200%載荷運行轉臺，采樣頻率為2 048 Hz。將樣本分成了87組，其中以第10組樣本為訓練樣本，其余為測試樣本，每組122 886個樣本。

為了驗證該模型的有效性，在實驗中使用激光干涉儀定期測量了轉臺的定位精度。實驗設備為雷尼紹XL-80激光干涉儀系統，該設備回轉測量時分度精度為±5 μm/m，重復精度為1 μm/m。如圖2(b)所示，使用激光干涉儀每60 h測一次共測得12組數據。每次測量時需要將轉臺平放在操作臺上，安裝測量儀器。試驗時，數控轉臺以10°為步長，正反交替轉兩圈獲得一次測量值，每次試驗測三次并取平均值。整個試驗系統的原理圖如圖2(c)所示。

圖2 數控轉臺精度衰退加速壽命試驗系統Fig.2 Accelerated life test system ofaccuracy decline of NC rotary table

3.2 試驗結果分析

通過對比分析試驗采集的87組原始數據，在第10組～第30組數據的振動幅值十分平穩，從第39組數據開始幅值出現較大的增加，并且出現規律性的變化，圖3(a)中對第10組、第39組和第50組數據進行了對比，從中可以看出，振動信號的幅值在全壽命周期中躍升了三個數量級。

圖3 第10組、第39組和第50組原始數據以及EEMD-PCA方法重構前后的第50組數據Fig.3 Raw data at file number 10、39 and 50 and Fiftieth sets of data before and after reconstruction by EEMD-PCA

獲取了數據后使用EEMD-PCA方法對原始數據降噪并進行信號重構，圖3(b)中以第50組數據為例，給出了重構前后的數據。從中可以看出，信號降噪和重構之后出現了比較明顯的能反映精度變化的特征。同時，圖4給出了表2中的統計特征量在代入重構信號時的計算結果，這些特征量對轉臺精度衰退的變化比較敏感，驗證了表2的合理性。

將測試樣本輸入到經過第10組數據訓練的HMM模型中，計算出一組對數似然概率值，該值反映了測試樣本與訓練樣本的匹配度，如圖5所示。選用第10組數據作為訓練樣本主要是考慮到轉臺在運轉初期有一個磨合的過程，其初期的精度不是轉臺運轉的最佳狀態，經過對比前后的數據，認為第10組數據可以代表正常精度狀態下的振動信號，并且由于本次試驗的采樣頻率較高，所以一組數據的信息量足夠大。

圖4 統計特征量計算結果Fig.4 Calculation results of statistical characteristic quantities

從圖中可以看出，在第39組數據時，匹配度有一個比較大的下降，可以判斷轉臺精度從此時開始逐漸下降，本文將從第39組開始計算出的概率值作為健康狀態指標H。同時對比圖5和圖6，可以看出在用激光干涉儀測量結果中第8次測量與第9次測量之間，轉臺的精度超過了設計允許的誤差精度值8"，時間是在第480 h～第540 h，對應HMM計算的似然函數是在第67組，因此可以認為函數值在-100為精度失效的閾值。

這里需要說明的是，圖6是利用激光干涉儀在數控轉臺精度衰退全壽命試驗中測量得到的，它反映了整個精度衰退的過程的全貌，而在轉臺實際運轉工作過程中，利用激光干涉儀獲得的僅僅是某一個或幾個點的精度值，它難以獲得對整個運轉周期完整的描述，只能判斷是否精度失效，而無法反映出轉臺當前的精度狀態。本文提出的算法模型通過對轉臺實際運轉過程中監測獲得的歷史數據進行分析，可以得到轉臺精度衰退的趨勢，描繪出轉臺所處的精度狀態，并且可以以某一個精度未失效的點為起始點，來預測轉臺精度的變化趨勢，判斷出精度的剩余壽命，由此判斷是否要對轉臺進行必要的調整。

圖5 數控轉臺精度衰退健康指標曲線Fig.5 Health indicators curve of accuracy decline of NC rotary table

圖6 激光干涉儀測量重復定位精度值Fig.6 Measurement accuracy of repeated positioning by laser interferometer

圖7 三種算法預測結果圖Fig.7 Prediction results of three algorithms

方法平均絕對誤差均方誤差計算時間/sPSO優化初始參數0．5750．88943．67未優化的粒子濾波1．4753．53110．33PSO優化粒子集1．3173．28339．89

表4給出以不同信號分組作為預測起始點的情況下三種算法的預測結果。其中，從第50組數據開始預測時，結合圖7可以看出，三種方法的不確定度為50%時的預測結果分別是71、75和60，則可得到其剩余壽命(RUL)分別為21、25和10。通過激光干涉儀測量精度變化獲得的實際的精度剩余壽命為17，差值分別為4、8和7。綜合對比不同預測起始點的預測結果后認為本文提出的方法更加接近實際的結果。圖8給出了本文提出的優化粒子濾波初始參數的方法在不同預測起始點的情況下，預測結果的PDF的分布曲線，可以看出起始點越大其分布越集中，說明其不確定度在變小，與實際情況是一致的。

表4 不同起始點剩余壽命預測結果

圖8 PSO優化粒子濾波算法預測結果的PDF圖Fig.8 PDF results of PF initial parameters optimizedby PSO

4 結 論

為了解決數控轉臺的精度衰退缺乏有效的描述方法且現有測量儀器只是測量工具而無法實現狀態評估等問題。本文引入數據驅動的方法，提出數控轉臺重復定位精度衰退趨勢預測模型，并得出如下的結論：

(1)本文區別于以往的狀態識別，將HMM算法以新的角度運用于數控轉臺的精度衰退的研究中，結合優化后的粒子濾波算法，建立了一種預測轉臺精度衰退趨勢的模型。在以第50組為預測起始點時，預測剩余壽命為21，實際剩余壽命為17，相差4，說明預測結果準確度較高。

(2)在使用粒子濾波算法的過程中，對粒子濾波算法的預測點初始參數使用PSO算法進行了優化，找出了最優的結果，并與未優化的粒子濾波算法以及用PSO優化粒子濾波粒子集方法的結果進行了對比分析，結果表明，雖然PSO優化初始參數的程序計算過程需要更長的時間，但是以第50組為預測起始點時，其預測結果與實際結果最接近，相差4，其他兩種方法分別相差8和7。其他預測起始點時也是相同的情況。

[ 1 ] YU C J, HUANG X D, FANG C G, et al. Research on damping and vibration characteristic of the large and precision NC rotary table[J]. Advanced Materials Research, 2010, 97/98/99/100/101: 1216-1222.

[ 2 ] ?VOKELJ M, ZUPAN S, PREBIL I. Non-linear multivariate and multiscale monitoring and signal denoising strategy using kernel principal component analysis combined with ensemble empirical mode decomposition method [J]. Mechanical Systems & Signal Processing, 2011, 25(7): 2631-2653.

[ 3 ] YU C J, HUANG X D, FANG C G. Research on dynamic characteristics of NC rotary table considering leakage factors of its hydrostatic guideway[J]. Journal of Mechanical Engineering Science, 2012, 226(11): 2674-2685.

[ 4 ] HENG A, ZHANG S, TAN A C C, et al. Rotating machinery prognostics: state of the art, challenges and opportunities [J]. Mechanical Systems & Signal Processing, 2009, 23(3): 724-739.

[ 5 ] BENKEDJOUH T, MEDJAHER K, ZERHOUNI N. Remaining useful life estimation based on nonlinear feature reduction and support vector regression [J]. Engineering Applications of Artificial Intelligence, 2013, 26(7): 1751-1760.

[ 6 ] LU C, CHEN J, HONG R, et al. Degradation trend estimation of slewing bearing based on LSSVM model [J]. Mechanical Systems & Signal Processing, 2016, 76/77:

353-366.

[ 7 ] MATEJ Z, SAMO Z, IVAN P. Non-linear multivariate and multiscale monitoring and signal denoising strategy using kernel principal component analysis combined with ensemble empirical mode decomposition method[J]. Mechanical Systems and Signal Processing, 2011, 25(7): 2631-2653.

[ 8 ] FENG Y, HUANG X, HONG R, et al. Residual useful life prediction of large-size low-speed slewing bearings-a data driven method[J]. Journal of Vibroengineering, 2015, 17(8): 4164-4179.

[ 9 ] 于春健.多物理場耦合的大重型數控回轉工作臺動靜態精度研究[D]. 南京：南京工業大學，2012.

[10] WANG D, MIAO Q, ZHOU Q, et al. An intelligent prognostic system for gear performance degradation assessment and remaining useful life estimation [J]. Journal of Vibration & Acoustics, 2015, 137(2): 12.

[11] 張西寧, 雷威, 李兵. 主分量分析和隱馬爾科夫模型結合的軸承監測診斷方法[J]. 西安交通大學學報,2017, 51(6): 1-7.

ZHANG Xining, LEI Wei, LI Bing. Bearing fault detection and diagnosis method based on principal component analysis and hidden markov model [J].Journal of Xi’an Jiaotong University, 2017, 51(6)： 1-7.

[12] 余龍華,王宏,鐘洪聲. 基于隱馬爾科夫模型的人臉識別[J]. 計算機技術與發展,2012, 22(2): 25-28.

YU Longhua，WANG Hong，ZHONG Hongsheng. Face recognition based on hidden markov model [J].Computer Technology and Development, 2012, 22(2): 25-28.

[13] RABINER L R. A tutorial on hidden Markov models and selected applications in speech recognition [J]. Readings in Speech Recognition, 1990, 77(2): 267-296.

[14] ARULAMPALAM M S, MASKELL S, GORDON N. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[J]. IEEE Transactions on Signal Processing, 2002, 50(2): 174-188.

[15] 王法勝, 魯明羽, 趙清杰, 等. 粒子濾波算法[J]. 計算機學報,2014,37(8): 1679-1694.

WANG Fasheng, LU Mingyu, ZHAO Qingjie, et al. Particle filtering algorithm [J]. Chinese Journal of Computer, 2014, 37(8): 1679-1694.

[16] 方正, 佟國峰, 徐心和. 粒子群優化粒子濾波方法[J]. 控制與決策,2007, 22(3): 273-277.

FANG Zheng, TONG Guofeng, XU Xinhe. Particle swarm optimized particle filter [J].Control and Decision, 2007, 22(3): 273-277.

[17] 王維博, 林川, 鄭永康. 粒子群算法中參數的實驗與分析[J]. 西華大學學報(自然科學版),2008, 27(1): 76-80.

WANG Weibo, LIN Chuan, ZHENG Yongkang. Experiment and analysis of parameters in particle swarm optimization [J]. Journal of Xi Hua University (Natural Science), 2008, 27(1): 76-80.

[18] 封楊, 黃筱調, 陳捷, 等. 大型回轉支承非平穩振動信號的EEMD-PCA降噪方法[J]. 南京工業大學學報(自然科學版), 2015, 37(3): 61-66.

FENG Yang, HUANG Xiaodiao, CHEN Jie, et al. EEMD-PCA based denosing method for non-stationary vibration signals of large-size slewing bearings [J].Journal of Nanjing Tech University(Natural Science) ,2015, 37(3):61-66.