Abstract
A realtime tool health assessment has a significant implication on reliable machining operations. This paper proposes a health assessment method for tools in milling machine using the MahalanobisTaguchi system (MTS) based on wavelet packet transformation and autoregression. In this method, the nonlinear and nonstationary vibration signal from the milling process is first decomposed using wavelet packet transforms. Then, an autoregression (AR) model is constructed for each coefficient of the reconstructed signal, and the parameters as well as variance of the remnants of each AR model are employed to form the initial feature matrix. Singular values of this feature matrix are obtained through singular value decomposition, at which point MTS is employed. In this study, MTS provides (1) a computational scheme based on the Mahalanobis distance for obtaining the health index of a tool, and (2) Taguchi methods to extract the key features and reduce the redundant ones. Finally, the performance and effectiveness of the proposed method are validated by vibration signals acquired from the milling machining process.
1. Introduction
In a modern machining system, tool health assessment systems are necessary to ensure high quality production and prevent downtime of mechanical tools as a result of catastrophic failures [1]. Tool health status directly affects the surface finish and dimensional accuracy of a product. Thus, tool health assessment is essential in machining operations [2].
One challenge in effective tool health assessment is selecting a proper feature space that can reflect comprehensive performance [3]. Nowadays, the tool vibration signals usually display strong nonlinear, nonGaussian and nonstationary features because of the influence of the nonlinear factors including clearance, friction, stiffness and so on, so it is difficult to extract a proper feature space which can accurately recognize the working condition of the tools.
The Autoregressive (AR) model is a practical method with simple calculations. It can extract the fault feature of the vibration signals effectively without a requirement of constructing a mathematical model and studying the fault mechanism of the system [4, 5]. Baillie and Mathew used the coefficients of AR model as the input features of an artificial neural network to realize correct classification of the bearing faults [6]. Ocak and Loparo also used the coefficients of AR model as features and adopted hidden Markov model (HMM) to classify different kinds of bearing abnormalities [7]. Weixiang Zhao and Cristina E. Davis demonstrated the distinct advantages of an AR model based feature extraction method for the classification of the chromatography data with time shifts [8]. These applications prove that the autoregressive model is effective to characterize the features of the vibration signal, and it is demonstrated that the coefficients of AR model are very sensitive to condition changes. However, as a model to characterize the stationary time series, estimation of the AR parameters may not be available for nonlinear and nonstationary signals [9]. Aim at this problem, in this paper, the feature vectors of the tool vibration signal are extracted with the combination of wavelet packet transform (WPT) and AR model. Compared with common timefrequency analysis methods, such as WignerVille distributions (WVD), the short time Fourier transformation (STFT), the wavelet transform (WT), etc., the WPT has particular advantages for characterizing signals with a complete multilevel decomposition in both time and frequency domains. Meanwhile, the WPT has particular advantages in extracting features that combine nonlinear and nonstationary characteristics, which is crucial for nonlinear and nonstationary vibration signal [4]. Thus, the WPT and AR models are integrated to ensure better efficacy in the fault feature extraction. Since the singular values are the intrinsic features of a feature vector matrix and have good stability [10], the SVD technique is further developed in this study to decompose the feature vector matrices and obtain singular values.
After extracting feature vectors, a large amount of data can be collected for processes analysis, but another challenge of realizing effective tool health assessment is how to deal effectively with the vast amount of data and build an intelligent assessment model.
The MahalanobisTaguchi system (MTS), introduced by G. Taguchi, is a prediction and diagnosis method which is used to analyze data patterns in multivariate systems [11]. It is a pattern information technology, which has been used in different applications to make quantitative decisions by constructing a multivariate measurement scale. In the MTS approach, Mahalanobis distance (MD) is used to measure the degree of abnormality of patterns, whereas the Taguchi methods which utilize orthogonal arrays (OA), and signaltonoise ratio (S/N), are used to evaluate the accuracy of predictions, and to optimize the system [12].The benefits of MTS as a pattern recognition tool can be summarized as follows [13, 14]: (1) It is a robust methodology that is insensitive to variations in multidimensional systems; (2) It can handle many different types of data sets and effectively consolidates the data into a useful metric; (3) Implementation of MTS requires limited knowledge of statistics; (4) It relies typically on simple arithmetic, contextual knowledge, and intuition; (5) It employs the Taguchi methods to determine the key features and reduce the rest; (6) Its success has been demonstrated in various practical applications. In this study, MTS is employed to provide a computational scheme based on the MD, and the key features are extracted using Taguchi methods.
2. Methodology
2.1. Feature extraction based on WPTAR
Wavelet packet transform is a powerful tool that has been applied extensively for feature extraction. The wavelet packet method is an expansion of classical wavelet transform that presents more possibilities for signal processing. WPT splits the signal utilizing both lowfrequency components and the highfrequency components, which can divide the full frequency band of signal into multilevel, so that each band’s signal contain more elaborate information about the original signal. This flexibility of a rich collection of abundant information with arbitrary timefrequency resolution allows extraction of features that combine nonstationary and stationary characteristics. For n levels decomposition, the WPT produces 2n different sets of coefficients. These coefficients are feature representations of the original signal in different wavelet packet bases [15, 16].
In this paper, a 3level WPT is utilized for feature extraction, and 8 different sets of coefficients are calculated. For each set of coefficient, an AR model is constructed and the parameters of the AR models are extracted to form the feature vector.
An AR model is simply a linear regression of the current value of the series against one or more prior values of the series. The AR model of an order $p$ is defined as:
where, the $n$th value $x\left(n\right)$ can be predicted by its previous $p$ values:$\mathrm{}{x}_{n1}\text{,}\text{}{x}_{n2}\text{,}\text{}\text{...,}\text{}{x}_{np}$, $a\left(k\right)$ are the AR coefficients and $e\left(n\right)$ is a Gaussian whitenoise series with zero mean and variance ${\sigma}^{2}$. The goal of an AR model is to estimate the AR coefficients that can fit the original data as much as possible through an optimization process [8]. The advantage of the AR model is that unlike a moving average model (and the moving average part of ARMA), its parameters can be determined by solving a linear set of equations. Generally, an AR model requires far fewer coefficients than the corresponding moving average model, and hence, it is more efficient in feature extraction [17].
2.2. MahalanobisTaguchi system
2.2.1. Mahalanobis distance (MD)
MD, introduced by P. C. Mahalanobis in 1936 [18], is a multivariate generalized measure used to determine the distance of a data point to the mean of a group. The advantage of the MD is that it is sensitive to the intervariable changes in the reference data. Therefore, it has traditionally been used to classify observations into different groups and diagnoses. It is calculated as follows [14, 18]:
1) Calculate the mean for each characteristic in the normal operation data set as:
2) Then, calculate the standard deviation ${s}_{i}$ for each characteristic ($i=$1, 2, 3,…):
3) Next, normalize each characteristic, form the normalized data matrix (${Z}_{ij}$), and take its transpose ${Z}_{ij}^{T}$:
4) Form the correlation matrix $C$ for normalized data. Calculate the matrix $C$ elements ${c}_{ij}$ as follows:
5) Finally, calculate MD as:
where ${C}^{1}$ is the inverse of the correlation matrix $C\text{,}$ which contains correlation coefficients between the variables, ${Z}_{ij}$ is the column vector of standardized variables, $M{D}_{j}$ is the Mahalanobis distance for the $j$th observation, and $k$ is the number of characteristics.
2.2.2. Taguchi methods
Taguchi methods are employed to extract key features. An orthogonal array (OA) is a table that lists a combination of characteristics, thereby enabling the effects of the presence or absence of a characteristic to be tested. One example of OA (${L}_{9}\left({2}^{8}\right)$) is as shown in Table 1.
Table 1Taguchi orthogonal array
Level  Factor  
$E$  1  2  3  4  5  6  7  8 
1  1  1  1  1  1  1  1  1 
2  1  1  1  1  1  2  2  2 
3  1  1  2  2  2  1  1  1 
4  1  2  1  2  2  1  2  2 
5  1  2  2  1  2  2  1  2 
6  1  2  2  2  1  2  2  1 
7  2  1  2  2  1  1  2  2 
8  2  1  2  1  2  2  2  1 
In MTS, characteristics in the OA have two levels. Level 1 represents the presence of a characteristic, whereas Level 2 represents the absence of a characteristic. In abnormal cases, MD values are calculated using the combination of characteristics dictated by the OA, and the largerthebetter signaltonoise ratio is calculated as follows [12]:
where ${\eta}_{q}$ is the signaltonoise ratio for the $q$th row of the OA, and $n$ is the number of observations under consideration. Afterwards, an average signaltonoise ratio at Level 1 (${t}_{1}$) and Level 2 (${t}_{2}$) of each characteristic is obtained. If ${t}_{1}{t}_{2}>\text{0}$, then this characteristic is useful; otherwise, it is useless for diagnosis.
2.3. CV calculation
This method adopts confidence value ($CV$), which was proposed by Qiu [19], as an index for evaluating the degree of health condition. A normalized function based on sigmoid function is proposed as follows to convert the $MD$ obtained by MTS to $CV$:
In Eq. (8), ${c}_{0}$ is the scale parameter determined by the $MD$ and the $CV$ value under normal state. The CVs under normal conditions, defined as $C{V}_{pre}$, can be denoted as:
where $mean\left(M{D}_{normal}\right)$ is the average of MDs under normal conditions. By converting the aforementioned equation, we can denote ${c}_{0}$ as follows:
2.4. Health assessment method based on WPTAR and MTS
As mentioned, the proposed health assessment scheme can be summarized as shown in Fig. 1.
Fig. 1Scheme of the proposed health assessment method
The procedure can be divided into two main parts:
1) WPTAR based feature extraction
The tool vibration signal can be decomposed to ${2}^{n}$ subbands with different frequencies by $n$level WPT. Coefficients of the reconstructed signal can be defined as ${c}_{1}\left(t\right)\text{,}\text{}{c}_{2}\left(t\right)\text{,...,}\text{}{c}_{j}\left(t\right)$ ($j={2}^{n}$), each of which implies different characteristic scale information. Therefore, the characteristics of the tool vibration signal can be obtained by extracting the characteristics of ${c}_{1}\left(t\right)\text{,}\text{}{c}_{2}\left(t\right)\text{,...,}\text{}{c}_{j}\left(t\right)$. The AR model for each coefficient ${c}_{i}\left(t\right)$ is constructed as:
where ${\phi}_{ik}$ ($k=$ 1, 2,…, $p$) and $p$, respectively, are the $k$th model parameters and model order of the AR model of the ith coefficient; ${e}_{i}\left(t\right)$ is the remnant of the model, a white noise sequence with zero mean value and variance ${\sigma}_{i}^{2}$ [10]. As the parameters ${\phi}_{ik}$ ($k=$1, 2,…, $p$) can reflect the inherent characteristics of a tool vibrating system and the variance of the remnant ${\sigma}_{i}^{2}$ is tightly related with the output characteristics of the system, ${\phi}_{ik}$ ($k=$1, 2,…, $p$) and ${\sigma}_{i}^{2}$ can be chosen as feature vectors ${A}_{i}=[{\phi}_{i1},{\phi}_{i2},\dots ,{\phi}_{ip},{\sigma}_{i}^{2}]$ to identify the tool condition in the milling machine [5]. Thereafter, the singular value decomposition (SVD) method is used to decompose the feature vector matrices and obtain singular values [20].
2) Health assessment
After the feature extraction, the MDs between the testing data and the benchmark space composed of the normal training data are determined as follows:
where ${R}^{1}$ is the inverse of the correlation matrix $R$, which contains correlation coefficients between the variables, ${x}_{j}$ is the column vector of standardized variables, $M{D}_{j}$ is the $MD$ for the $j$th observation, and $k$ is the number of characteristics.
Then, the Taguchi method is employed to optimize the $n$dimensional feature vector. The orthogonal table is established, and the largerthebetter signaltonoise ratio is calculated. By comparing the average signaltonoise ratio between Level 1 (${t}_{1}$) and Level 2 (${t}_{2}$) of each characteristic, the principal characteristics are extracted and the redundant features are reduced. With the identified key features, the MDs are recalculated for health assessment and finally converted into CVs, thereby achieving the health condition assessment.
3. Case study
The data in this study are sampled by a vibration sensor that measured runs on a milling machine. Fig. 2(a) shows the mounting position of the vibration sensor. Tool wear occurs in different forms. In our experiments, we measured the flank wear VB as a generally accepted parameter for evaluating tool wear, as shown in Fig. 2(b). Two experiments with cut depth of 1.5 mm feed rate of 0.25 mm/rev, and cut depth of 0.75 mm feed rate of 0.25 mm/rev were chosen as example 1 and 2. Cast iron is used as the material. Here, we mainly focus on example 1.
Fig. 2a) Mounting position of the vibration sensor, b) Measure of flank wear VB
a)
b)
3.1. Feature extraction
In example 1, a total of 40 samples have been acquired under normal conditions. The first 30 normal samples are used to construct the normal space and the others are used for testing. Among the samples, 10 have a VB of 0.18 mm, and these samples are also acquired for the testing process. Fig. 3 illustrates a sample of WPT results of the normal vibration signal.
Fig. 3A sample of WPT results of the normal tool signal
Table 2Extracted singular values for model training
Number  State  1  2  3  4  5  6  7  8 
1  Normal condition  12.6596  11.6476  2.0088  1.4113  0.8029  0.2152  0.1703  0.1328 
2  12.6885  11.6316  1.9789  1.3542  0.7955  0.2255  0.1741  0.1356  
3  12.6753  11.6344  1.9908  1.3871  0.8048  0.2168  0.1707  0.1353  
4  12.6696  11.6170  1.9811  1.3794  0.8038  0.2169  0.1782  0.1333  
5  12.6602  11.6247  1.9926  1.3935  0.8071  0.2120  0.1749  0.1238  
6  12.7062  11.6544  1.9851  1.3911  0.7982  0.2188  0.1701  0.1367  
7  12.6904  11.6415  1.9842  1.3958  0.8019  0.2133  0.1681  0.1425  
8  12.7053  11.6178  1.9692  1.3776  0.8057  0.2307  0.1792  0.1370 
The AR model is constructed for each coefficient of the reconstructed signal ${c}_{1}\left(t\right)\text{,}\text{}{c}_{2}\left(t\right)\text{,...}\text{}\text{,}\text{}{c}_{8}\left(t\right)$. In this case, the order, m, is determined as 7, so the AR parameters ${\phi}_{ik}$ ($k=$1, 2,… , 7) and ${\sigma}_{i}^{2}$ are chosen to form the initial feature matrix. Then, SVD is employed to obtain an 8dimensional feature vector of singular values. Some of the extracted singular values are listed in Tables 2 and 3.
Table 3Extracted characteristic vectors for testing
Number  State  1  2  3  4  5  6  7  8 
1  Normal condition  12.7162  11.6418  1.9826  1.3946  0.8055  0.2281  0.1677  0.1303 
2  12.7171  11.6276  1.9775  1.3630  0.8071  0.2273  0.1749  0.1199  
3  12.7272  11.6938  1.9952  1.3934  0.8066  0.2262  0.1592  0.1271  
4  12.7631  11.6370  1.9657  1.3776  0.8092  0.2355  0.1642  0.1317  
1  Worn condition  12.6851  11.7671  1.8333  1.3133  0.7581  0.4003  0.1844  0.1161 
2  12.6662  11.7570  1.8504  1.3368  0.7604  0.4170  0.1797  0.1168  
3  12.6434  11.7149  1.8185  1.3063  0.7584  0.3785  0.1872  0.1193  
4  12.6591  11.7237  1.8676  1.3518  0.7586  0.3786  0.1816  0.1210 
3.2. Health assessment based on MTS
3.2.1. MD calculation
After feature extraction, the $MD$ values of the testing samples under normal and worn conditions are calculated. Table 4 presents the experimental results. These results prove that the MD method is highly effective for fault diagnosis.
Table 4Mean, minimum, and maximum of MDs
Normal condition  Worn condition  
Mean  2.6454  337.1591 
Minmax  1.01535.3637  252.4657408.8422 
3.2.2. Optimization based on the Taguchi method
The extracted feature vector is constructed by 8dimensional singular values. Therefore, the ${L}_{9}\left({2}^{8}\right)$ OA is employed to identify the useless characteristics. We used 30 samples of normal conditions ($VB=$0 mm) and 10 samples of worn conditions ($VB=$0.18 mm) in the Taguchi experiment. The results are reported in Table 5.
Table 5Results of Taguchi method
1  2  3  4  5  6  7  8  SNR  
1  1  1  1  1  1  1  1  1  21.24 
2  1  1  1  1  1  2  2  2  20.84 
3  1  1  2  2  2  1  1  1  15.75 
4  1  2  1  2  2  1  2  2  21.86 
5  1  2  2  1  2  2  1  2  15.94 
6  1  2  2  2  1  2  2  1  17.92 
7  2  1  2  2  1  1  2  2  19.23 
8  2  1  2  1  2  2  2  1  13.14 
9  2  1  1  2  2  2  1  2  13.98 
${t}_{1}$  18.93  17.36  19.48  17.79  19.81  19.52  16.73  17.01  
${t}_{2}$  15.45  18.57  16.40  17.75  16.13  16.36  18.60  18.37  
${t}_{1}{t}_{2}$  3.48  –1.21  3.08  0.05  3.68  3.15  –1.87  –1.35  
Useful  Y  N  Y  Y  Y  Y  N  N 
According to the rules of the Taguchi method, the 2nd, 7th, and 8th singular values should be excluded as redundant factors. With the identified key features, the $MD$ values for all of the testing datasets were recalculated. Table 6 presents a comparison between the initial and optimized MDs for each dataset. The result confirms that the Taguchi methods are highly effective in characteristic optimization.
Table 6Comparison between initial and optimized MDs
Normal condition  
Initial  Optimized  
Mean  2.6454  2.0733 
Min–max  1.01535.3637  0.82843.5922 
Worn condition  
Initial  Optimized  
Mean  337.1591  130.8750 
Minmax  252.4657408.8422  112.5295155.6077 
3.2.3. Health assessment
A total of 80 samples from example 1 and 280 samples from example 2 are chosen to verify the effectiveness of the proposed method. In example 1, every 10 samples are in the same $VB$. The $MD$ values and CVs of each group are calculated, as shown in Fig. 4.
Blue curve in Fig. 5 indicates the health degradation of the tool. The first 20 samples with CVs above 0.9 can be considered as “normal stage”, which shows that the tool is in good health condition. Thereafter, the MDs increase while the CVs decrease, and the tool is considered to be at “danger stage”. Once the $CV$ drops below the threshold of 0.7, the tool is considered to be at “failure stage”.
Fig. 4MDs and CVs of the tools in example 1
a)
b)
In example 2, every 40 samples are in the same $VB$. The $MD$ values and CVs of each group are shown in Fig. 5. The blue curve shows the same trend as in example 1. The first 40 samples are in the “normal stage”, with CVs above 0.9. As the tool health degrading, the CVs fall at the same time, which can be considered as “danger stage”. The last 80 samples are in the “failure stage”, with CVs bellow 0.7.
Thus, the health condition of the tool can be coincidentally represented by the $CV$. The experiment results demonstrate that the proposed health assessment method based on WPTAR and MTS performs efficiently.
Fig. 5MDs and CVs of the tools in example 2
a)
b)
4. Conclusions
This paper presents a health assessment method for tools in milling machine using MTS based on WPTAR. Targeting nonlinear and nonstationary signals from the milling process, the WPTAR method is employed in feature extraction and forms the initial feature matrix. Singular values of this feature matrix are obtained via SVD. Then, MTS is employed for health assessment, which has been verified by vibration signals obtained from the milling machining process. Generally speaking, the proposed approach has two benefits. First, it extracts the suitable features that evidently reflect the health state of tools because the combination of WPT and AR prevents drawbacks in processing nonlinear and nonstationary signals. Second, it is a reliable realtime multivariate analytical method that offers a systematic way to determine the principal characteristics through MTS.
References

Kaya B., Oysu C., Ertunc H. M. 2011 Forcetorque based online tool wear estimation system for CNC milling of Inconel 718 using neural networks. Advances in Engineering Software, Vol. 42, 2011, p. 7684.

Koshy P., Dumitrescu P., Ziada Y. Novel methods for rapid assessment of tool performance in milling. International Journal of Machine Tools & Manufacture, Vol. 44, 2004, p. 15991605.

Pan Y., Chen J., Li X. Bearing performance degradation assessment based on lifting wavelet packet decomposition and fuzzy cmeans. Mechanical Systems and Signal Processing, Vol. 24, 2010, p. 559566.

Li Z. X., Yan X. P., Yuan C. Q., et al. Virtual prototype and experimental research on gear multifault diagnosis using waveletautoregressive model and principal component analysis method. Mechanical Systems and Signal Processing, Vol. 25, 2011, p. 25892607.

Cheng J. S., Yu D. J., Yang Y. A fault diagnosis approach for roller bearings based on EMD method and AR model. Mechanical Systems and Signal Processing, Vol. 20, 2006, p. 350362.

Baillie D. C., Mathew J. A comparison of autoregressive modeling techniques for fault diagnosis of rolling element bearings. Mechanical Systems and Signal Processing, Vol. 10, Issue 1, 1996, p. 117.

Ocak H., Loparo K. A. HMMbased fault detection and diagnosis scheme for rolling element bearings. Journal of Vibration and Acoustics, Vol. 127, 2005, p. 299366.

Zhao W. X., Davis C. E. Autoregressive model based feature extraction method for time shifted chromatography data. Chemometrics and Intelligent Laboratory Systems, Vol. 96, 2009, p. 252257.

Wang G. F., Liu C., Cui Y. H. Clustering diagnosis of rolling element bearing fault based on integrated Autoregressive/Autoregressive Conditional Heteroscedasticity model. Journal of Sound and Vibration, Vol. 331, 2012, p. 43794387.

Wang Y. J., Kang S. Q., Jiang Y. C., et al. Classification of fault location and the degree of performance degradation of a rolling bearing based on an improved hyperspherestructured multiclass support vector machine. Mechanical Systems and Signal Processing, Vol. 29, 2012, p. 404414.

Reséndiz E., RullFlores C. A. MahalanobisTaguchi system applied to variable selection in automotive pedals components using Gompertz binary particle swarm optimization. Expert Systems with Applications, Vol. 40, 2013, p. 23612365.

Cudney E. A., Hong J., Jugulum R., et al. An evaluation of MahalanobisTaguchi System and neural network for multivariate pattern recognition. Journal of Industrial and Systems Engineering, Vol. 1, Issue 2, 2007, p. 139150.

Hu J. Q., Zhang L. B., Liang W. Dynamic degradation observer for bearing fault by MTSSOM system. Mechanical Systems and Signal Processing, Vol. 36, 2013, p. 385400.

Soylemezoglu A., Jagannathan S., Saygin C. MahalanobisTaguchi System as a multisensor based decision making prognostics tool for centrifugal pump failures. IEEE Transactions on Reliability, Vol. 60, Issue 4, 2011, p. 864878.

Yusuff A. A., Fei C., Jimoha A. A., Munda J. L. Fault location in a series compensated transmission line based on wavelet packet decomposition and support vector regression. Electric Power Systems Research, Vol. 81, 2011, p. 12581265.

Ekici S., Yildirim S., Poyraz M. Energy and entropybased feature extraction for locating fault on transmission lines by using neural network and wavelet packet decomposition. Expert Systems with Applications, Vol. 34, 2008, p. 29372944.

Wang X. Y., Makis V. Autoregressive modelbased gear shaft fault diagnosis using the KolmogorovSmirnov test. Journal of Sound and Vibration, Vol. 327, 2009, p. 413423.

Mahalanobis P. C. On the generalized distance in statistics. Proceedings, National Institute of Science of India, Vol. 21, 1936, p. 4955.

Qiu H., Lee J., Lin J., Yu G. Robust performance degradation assessment methods for enhanced rolling element bearing prognostics. Advanced Engineering Informatics, Vol. 17, 2003, p. 127140.

Feiyun C., Jin C., Guangming D., Fagang Z. Shorttime matrix series based singular value decomposition for rolling bearing fault diagnosis. Mechanical Systems and Signal Processing, Vol. 34, 2013, p. 218230.
About this article
This study is supported by the National Natural Science Foundation of China (Grant Nos. 61074083, 51105019), as well as the Technology Foundation Program of National Defense (Grant No. Z132013B002).