Virtual sensing for gearbox condition monitoring based on extreme learning machine

2.1. Kernel principal component analysis

Principal component analysis (PCA) allows linear dimensionality reduction. However, if the data has more complicated structures which cannot be well represented in a linear subspace, traditional PCA will not be very helpful. Fortunately, kernel principal component analysis (KPCA) allows us to generalize traditional PCA to nonlinear dimensionality reduction, which is a nonlinear version of PCA and has been widely used for feature selection and fusion applications. The key idea of KPCA is to define a nonlinear transformation $\varphi \left(\right)$ which transforms the sample data into a high-dimensional data space. Then each data point $X_i$ is projected to a point $\varphi \left(X_i\right)$. Next, we can perform traditional PCA in the new feature space [15]. It transforms a set of observations of possible correlated variables into a set of uncorrelated variables called principal components. The first principal component has the largest variance, and each succeeding principal component has comparative lower variance orthogonal to the preceding principal components. The first several principal components can represent the original data with minimal mean squared approximation error, and thus KPCA can be used in dimensionality reduction.

Mathematically, given a set of input vectors ${\left(X_i\left(1\right),X_i\left(2\right),\dots ,X_i\left(m\right)\right)}^T$, $i=$1, 2,…, $p$, the sample data $X_i$ is mapped into $\varphi \left(X_i\right)$ via the nonlinear kernel function $\varphi \left(\right)$, i.e. With the assumption of centered data $\frac{1}{p}\sum _{i=1}^p\varphi \left(X_i\right)=0$, the principal components are obtained by solving eigenvalue problem in KPCA:

where $C$ is the sample covariance matrix of $\varphi \left(X_i\right)$, ${\lambda }_i$ is one of the eigenvalues of covariance matrix $C$, and $u_i$ is the corresponding eigenvector. The covariance matrix is constructed as:

Define a Gram matrix $K$ with its elements as:

where $x_i$ and $x_j$ are the sample vectors. Assuming $k\left(\right)$ is a symmetric kernel function, the dot production in Eq. (3) can be replaced by a kernel function $k\left(\right)$ based on the Mercer’s theorem. Since the data points need to be centered in the feature space, the centered kernel matrix $\tilde{K}$ is defined as [16]:

where $A_p$ is a $p\times p$ matrix with ${\left(A_p\right)}_{ij}=1/p$. The eigenvalue Eq. (1) can be rewritten as:

Then the $n$th kernel principal component is readily obtained by projecting the observations in the direction of the $n$th eigenvector [15]:

Since the number of eigenvectors is the same as the sample size in KPCA, it can deal with nonlinear problems which cannot be solved by PCA. By calculating the accumulated contribution rate (e.g., 95 %), the number of the most significant principal components can be selected for dimensionality reduction:

2.2. Support vector regression model

Support vector regression (SVR) is the term used when Support vector machines (SVMs) are used to solve nonlinear regression estimation problems [17], which is based on statistics learning theory and has already been widely applied in most fields and made great results. Comparing with other data mining techniques such as artificial neural networks (ANN), it reveals good generalization capability and needs less training samples [18]. SVR transforms the original feature space into a higher dimensional space to determine an optimal hyperplane by maximizing the separation distances among the classes. Given an input training data set $z\in \chi$, the transformed higher dimensional feature space can be obtained as:

where $\varphi$ is the transformation function. A hyperplane $f\left(z’\right)=$ 0 can be formulated as [18]:

where $\mathbf{\tau }$ is a $n$-dimensional vector and $b$ is a scalar. The vector $\mathbf{\tau }$ and scalar $b$ are used to define the position of the separating hyperplane. The hyperplane is built to maximize the distance $D$ among the closest classes through the following optimization:

where $y_i$ is the class labeler. For example, it is labeled as {–1, 1} for two classes. Taking into account the noise with slack variables ${\xi }_i$ and error penalty $C$, Eq. (10) can be rewritten as [18]:

The hyperplane can be determined as the following sign function ($sgn\left(t\right)=$ 1 for $t\ge$ 0, and $sgn\left(t\right)=$ –1 for $t<$ 0). The linear decision function is given by:

where $\alpha$ is the Lagrange multiplier. The hyperplane function can be determined by kernel function $K\left(z_i,z_j\right)={\varphi }^T\left(z_i\right)\varphi \left(z_j\right)$ by computing the inner products without specifying the explicit form of the transformation function. Different kernels can be formulated such as linear, polynomial, Gaussian RBF, and Sigmoid kernel functions. Accordingly, the associated decision function for regression analysis is expressed as [19]:

2.3. Extreme learning machine model

ELM [20–22] was originally proposed for the single hidden-layer feedforward neural networks (SLFNs) and was then extended to the generalized SLFNs where the hidden layer need not be neuron alike. ELMs have both universal approximation and classification capabilities; and they can also build a direct link between multiple theories (specifically, ridge regression, optimization, neural network generalization performance, linear system stability, and matrix theory). Consequently, ELMs, which can be biologically inspired, offer significant advantages such as fast learning speed, ease of implementation, and minimal human intervention. They thus have strong potential as a viable alternative technique for large-scale computing and machine learning.

Compared with other machine learning methods including back-propagation (BP) and SVMs, the ELM methods have the following advantages. Firstly, since the support vectors obtained by SVM are much larger than the required hidden neurons in ELM, the testing time spent SVMs for the same testing data set is much longer than the ELM. In addition, ELM only needs to set the number of hidden neurons, and it does not need to adjust the input weights and the bias of hidden neurons during the implementation of the algorithm, which leads to producing a unique optimal solution. However, it should be noted that in order to obtain as good performance as possible for SVM, long time effort has been made to find the appropriate parameters for SVM. That means, after trained and deployed the ELM may react to new observations much faster than SVMs in such real application.

The ELM for SLFNs shows that hidden nodes can be randomly generated. The input data is mapped to $L$-dimensional ELM random feature space, and the network output is:

where $\mathbf{\beta }={\left[{\beta }_1,{\beta }_2,\dots ,{\beta }_L\right]}^T$ is the vector of the output weights between the hidden layer of $L$ nodes and the output node and $\mathbf{h}\left(x\right)=\left[h_1\left(x\right),h_2\left(x\right),\dots ,h_L\left(x\right)\right]$ is the output (row) vector of the hidden layer with respect to the input $x$, and $h_i\left(x\right)$ is the output of the $i$th hidden node. $\mathbf{h}\left(x\right)$ actually, maps the data from the $d$-dimensional input space to the $L$-dimensional hidden-layer feature space (ELM feature space) $H$, and thus, $\mathbf{h}\left(x\right)$ is indeed a feature mapping.

Given $N$ training samples $\left\{\right(x_i,t_i){\}}_{i=1}^N$, the ELM can resolve the following learning problem:

where $\mathbf{T}={\left[t_1,t_2,\dots ,t_N\right]}^T$ are target labels, and $\mathbf{H}={\left[{\mathbf{h}}^T\left(x_1\right),{\mathbf{h}}^T\left(x_2\right),\dots ,{\mathbf{h}}^T\left(x_N\right)\right]}^T$.

Different from traditional learning algorithms, ELM tends to reach not only the smallest training error but also the smallest norm of output weights. According to Bartlett’s theory [23], for feedforward neural networks reaching smaller training error, the smaller the norms of weights are, the better generalization performance the networks tend to have. ELM is to minimize the training error as well as the norm of the output weights [20, 21].

Minimize:

The minimal norm least square method instead of the standard optimization method was used for calculating the output weights $\mathbf{\beta }$ in the original implementation of ELM [20, 21]:

where $H^{†}$ is the Moore-Penrose generalized inverse of matrix $\mathbf{H}$ [24, 25]. Different methods can be used to calculate the Moore–Penrose generalized inverse of a matrix: orthogonal projection method, orthogonalization method, iterative method, and singular value decomposition (SVD) [25].

The theoretical background of KPCA and machine learning methods discussed here, including SVR and ELM, forms the basis of virtual gearbox condition sensing model formulated in the next section.

3.1. Data acquisition and feature extraction

During operation process of mechanical equipment, online measurements can acquire the in-process parameters such as accelerometer and tachometer signals reflecting gearbox conditions. Due to low signal to noise ratio (SNR), it is usually difficult to model the relationship between raw measurement and gearbox condition. To tackle this problem, effective feature extraction techniques are usually performed to reduce data dimension without losing gearbox defect signatures. In this study, 21 features indictors or condition indictors (CIs) from time domain, frequency domain, and time-frequency domain are investigated. Time domain methods involve statistical features such as root mean square (RMS), Kurtosis (KT), crest factor (CF), peak-to-peak value (P2P). RMS is a measure for the magnitude of a varying quantity. It is also related to the energy of a signal. Kurtosis indicates the spikiness of a signal. Features from the frequency domain provide another perspective of gearbox mechanical components’ conditions, and reveal information that are not included in the statistical domain. In frequency domain, spectral skewness and spectral kurtosis are extracted, where $S\left(f_i\right)$ is the power spectrum density obtained using the Welch method. In time-frequency domain, wavelet transform can be used for signal denoising and feature extraction. The wavelet coefficient with higher energy is selected which is related to the bearing defect characteristic frequencies. Thus, the energy of the selected wavelet coefficient is also extracted as a feature.

3.2. Feature selection and fusion

There are an overwhelming number of features extracted from the raw measurements. In general, the extracted features can be viewed as a high-dimensional multivariate matrix composed of several feature vectors. It is not feasible to input the above matrix to a virtual sensing model without dimension reduction because of the curse of dimensionality and the high correlation between vectors. For improved computational efficiency in virtual sensing model, a proper feature selection and fusion strategy is needed to lower the dimension of a feature space.

By implementing feature selection and fusion algorithms, the complexity of modeling process could be reduced and new feature vectors are reconstructed. Different representative dimensional reduction techniques are investigated for feature selection and fusion, including KPCA, PCA and LPP algorithms. Generally, it is difficult to determine which feature is more sensitive to gearbox conditions. The goal of feature selection and fusion is to preserve as much of the relevant information as possible by removing redundant or irrelevant information in acquired sensory signals. The top ranked features of these three schemes (e.g. KPCA, PCA, LPP, etc.) are then selected and fused into the computational model to infer the actual quantities of gearbox mechanical components’ conditions. Afterwards, their performance is evaluated by comparing the predicted gearbox conditions with actual offline measurement in terms of model accuracy.

3.3. ELM based virtual sensing model

Given the complex relationship between fused features and actual quantities of gearbox mechanical components’ conditions, it is difficult to describe it in an explicit analytic form. By exploiting the underlying structure of data measurements, the paper utilizes a novel machine learning tool-ELM model, which has significant advantages such as fast learning speed, ease of implementation, and minimal human intervention compared with ANN and SVM in real applications. Thus, in this paper, the ELM model is used to investigate the dependency between fused features and gearbox conditions.

As previously mentioned, the parameters of ELM model mainly involve a tunable parameter of hidden neurons compared with SVR algorithm. In general, the proper hidden neurons number is obtained through a lot of artificial experiments or cross validation (CV) method. In this paper, the CV method is adopted to find the optimal hidden neurons number. Additionally, in order to validate the effectiveness of ELM method, SVR based virtual sensing scheme is investigated and compared in terms of model accuracy. And the selection of parameters and kernel functions in SVR model is determined using grid search algorithm whose evaluation of model performance follows leave-one-out cross-validation method.

4.1. Experimental setup and data collection

In this chapter, data from a spiral bevel gear case study conducted on the spiral bevel gear test facility are used to validate the presented method. And the experimental analyses of different virtual sensing schemes on three data sets are carried out in MATLAB 7.11.0 environment running in Intel(R) Xeon(R) E5-2650 v2, 2.60-GHZ CPU with 64-GB RAM.

Vibration data from experiments performed on the spiral bevel gear test facility was reprocessed for this analysis. A more detailed description of the test rig and test procedure is given in [2]. The rig is used to quantify the performance of gear materials, gear tooth design and lubrication additives on the fatigue strength of gears. During this testing, vibration condition indictors (CIs) and oil debris monitoring were used to detect pitting damage on spiral bevel gears.

The tests consisted of running the gears under load through a “back to back” closed loop torque regenerative system. Accelerometers were installed on the right and left side of the gearbox per Fig. 3. Vibration data was collected once per minute using a sampling rate of 100 kHz for 2 seconds duration. Shaft speed was measured by an optical sensor once per each gear shaft revolution, generating time synchronous averages (TSA) on the gear shaft (36 teeth). The pinion, on which the damage occurred, has 12 teeth. The tests were performed for a specific number of hours or until surface fatigue occurs. In this paper, data collected from three gears were used for performance comparison of virtual sensing schemes based on different feature selection and fusion methods. At test completion, destructive pitting could be observed on the teeth of the pinions (see Fig. 3).

Fig. 3a) The bevel gear test rig and bevel gears in [2], b) damaged spiral bevel gear in experiment 1, c) damaged spiral bevel gear in experiment 3

a)

b)

c)

4.2. Data processing

TSA data was processed with gear CI algorithms presented in [26] and [27] to compute the following CIs:

(1) TSA: RMS, Kurtosis (KT), Peak-to-Peak (P2P), Crest Factor (CF); (2) Residual RMS, KT, P2P, CF; (3) Energy Operator RMS, KT; (4) Energy Ratio; (5) FM0; (6) Sideband Level factor; (7) Narrowband (NB) RMS, KT, CF; (8) Amplitude Modulation (AM) RMS, KT; (9) Derivative AM KT; (10) Frequency Modulation (FM) RMS, KT. However, not all the CIs generated from TSA data were good candidates for virtual sensing model. For the purpose of prognostics, one is interested in selecting the CIs that have shown a good trending correlation. In order to select the best CIs, correlation coefficients of the CIs with the time index were computed. The following 6 CIs with correlation coefficients over 0.5 were selected for virtual sensing model: (1) residual RMS, (2) energy operator RMS, (3) FM0, (4) narrowband kurtosis, (5) amplitude modulation kurtosis, and (6) frequency modulation RMS. Next, different dimensional reduction techniques including KPCA, PCA and LPP are performed to exploit the 6 features by fusing the extracted 21 features. In all feature selection and fusion algorithms, the reduced dimensions are set: $k=$4 by calculating the accumulated contribution rate (e.g. 95 %), and kernel functions are optimized using grid search algorithm, respectively.

4.3. Performance evaluation

The selected and fused features obtained by KPCA, PCA and LPP are fed into the ELM model to infer the gearbox mechanical components’ conditions from the in-process vibration signals. And its parameter, hidden neurons number, is optimized by a 10-fold cross validation process. Afterwards, the performance of the proposed virtual sensing scheme is compared with SVR based virtual sensing scheme, in which two hyperparameters, the cost parameter $C$ and the Gaussian kernel parameter $\gamma$, are selected using grid search method in the cross validation process to prevent overfitting. A total of three sets of gearbox life test data (e.g. Y1, Y2, and Y3, etc.) are available. Take dataset Y1 as example. The predicted spiral bevel gear ODM using different virtual sensing schemes including KPCA, PCA and LLP algorithms with SVR and ELM models are illustrated in Figs. 4-6, respectively. To compare their performance, the actual spiral bevel gear ODM measured offline is also included. It is found that the predicted spiral bevel gear ODM by the virtual sensing model follows the trend of the truth data well.

Fig. 4Performance comparison of different virtual sensing schemes based on PCA method

Fig. 5Performance comparison of different virtual sensing schemes based on LPP method

Fig. 6Performance comparison of different virtual sensing schemes based on KPCA method

To quantitatively evaluate the performance of the proposed virtual sensing model, different criteria are investigated including Pearson Correlation coefficient (PCC), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE). PCC is a statistical measure of independence of two or more random variables which is defined as:

where $x$ is the actual spiral bevel gear ODM, and $\widehat{x}$ is the predicted spiral bevel gear ODM using the virtual sensing model. The model with the highest correlation coefficient is considered as the one of the best. RMSE is defined as the square root of the average of the square of all difference between estimated spiral bevel gear ODM $\widehat{x}$ and actual spiral bevel gear ODM $x$:

MAE is defined as the mean of all absolute difference between estimated spiral bevel gear ODM $\widehat{x}$ and actual spiral bevel gear ODM $x$:

MAPE is defined as the mean of all absolute percentage differences between estimated spiral bevel gear ODM $\widehat{x}$ and actual spiral bevel gear ODM $x$:

Fig. 7Performance comparison of different virtual sensing schemes using different criteria

Next, three different virtual sensing schemes with SVR and ELM models are quantitatively evaluated according to different criteria including PCC, RMSE, MAE, and MAPE. The performance of these three virtual sensing schemes is compared as shown in Fig. 7. In general, the larger the PCC value, the better the model performance, while the less the RMSE/MAE/MAPE value, the better the model performance. From the perspective of estimation accuracy, it is found from Table 1 that the predicted gearbox conditions using KPCA model follows the trend of the truth data best compared with other feature selection and fusion methods. In addition, ELM based virtual sensing scheme can obtain better prediction performance in comparison with SVR based virtual sensing scheme.