Health assessment of rotary machinery based on integrated feature selection and Gaussian mixed model

2.1. Empirical mode decomposition

The method of EMD is developed from the assumption that any signal consists of different simple intrinsic modes of oscillation. Each signal can be decomposed into a number of IMFs, each of which must satisfy the following definition [17]: (1) In the whole data set, the number of extreme crossings should either be equal or differ by at most one to the number of zero crossings; (2) At any point, zero is the mean value of the envelope defined by the local maxima and the envelope defined by the local minima.

An IMF represents a simple oscillatory mode compared with the simple harmonic function. With the definition, any signal $x\left(t\right)$ can be decomposed as follows [18]:

Step 1: Identify all the local extrema, and then connect all the local maxima by a cubic spline line as the upper envelope.

Step 2: Repeat the procedure for the local minima to produce the lower envelope. The upper and lower envelopes should cover all the data between them.

Step 3: The mean of upper and low envelope value is designated as $m_1\left(t\right)$; and the difference between the signal $x\left(t\right)$ and $m_1\left(t\right)$ is the first component, $h_1\left(t\right)$, i.e.:

Ideally, if $h_1\left(t\right)$ is an IMF, then $h_1\left(t\right)$ is the first component of $x\left(t\right)$.

Step 4: If $h_1\left(t\right)$ is not an IMF, $h_1\left(t\right)$ is treated as the original signal and repeat step 1, 2, 3; then:

After repeated sifting, i.e. up to $k$ times, $h_{1k}\left(t\right)$ becomes an IMF, that is:

Then, it is designated as:

The first IMF component from the original data $c_1\left(t\right)$ should contain the finest scale or the shortest period component of the signal.

Step 5: Separate $c_1\left(t\right)$ from $x\left(t\right)$, we could get:

$r_1\left(t\right)$ is treated as the original data and repeat the above processes; therefore, the second IMF component $c_2\left(t\right)$ of $x\left(t\right)$ could be got. Let’s repeat the process as described above for $n$ times, then $n$-IMFs of signal $x\left(t\right)$ could be got. Then:

The decomposition process can be stopped when $r_n\left(t\right)$ becomes a monotonic function from which no more IMFs can be extracted.

With the definition, any signal can be decomposed [19]. We finally obtain:

Thus, one can achieve a decomposition of the signal into $n$-empirical modes and a residue $r_n\left(t\right)$, which is the mean trend of $x\left(t\right)$.

EMD decomposes signal into several IMFs one by one iteratively. The new IMF $c_n\left(t\right)$ is decomposed from the previous residue $r_{n-1}\left(t\right)$. Each IMF $c_1\left(t\right)$, $c_2\left(t\right)$,…, $c_n\left(t\right)$ contains lower-frequency oscillations than the prior-extracted one. With the higher frequency parts are separated from signal $x\left(t\right)$, each residue $r_1\left(t\right)$, $r_2\left(t\right)$,…,$\mathrm{}{r}_n\left(t\right)$ is more smooth than the prior one. Until $r_n\left(t\right)$ becomes a monotonic function from which no more IMFs can be extracted, the process of EMD can be stopped. Therefore, the final residue $r_n\left(t\right)$ represents the central tendency of signal $x\left(t\right)$, and is usually a monotonic trend rather than white noise.

EMD has been proved as a fast, effective, self-adaptive method for nonlinear and non-stationary time series analysis, but there is still a drawback, which is worth noting: the end effect, appearing distortion at the end of signal in the decomposition process. In this paper，the mirror periodic extending method is employed to solve this problem.

2.2. Singular value decomposition

SVD is a powerful and effective feature extraction tool in linear algebra, and has been used for fault diagnosis in rotating machinery [21]. The SVD is defined as follows:

$X$ denotes an $m\text{×}n$ matrix with $m\ge n$. According to the SVD theorem, the matrix $X$ can be decomposed in the form given by [22]:

where $U\left(m\text{×}m\right)$ and $V\left(n\text{×}n\right)$ are orthogonal matrices, and $\sigma$ is an $m\text{×}n$ diagonal matrix of singular values (${\sigma }_{ij}\ne 0$, if $i=j$ and ${\sigma }_{11}\ge {\sigma }_{22}\ge \cdots \ge 0$). The columns of the orthogonal matrix $U$ and $V$ are called the left and right singular vectors, respectively. An important property of $U$ and $V$ is that they are mutually orthogonal. The singular values represent the importance of individual singular vector in the matrix composition. In other words, a singular vector that corresponds to larger singular values contains more information about the structure of the pattern embedded in the matrix than the other singular vector.

2.3. Principal component analysis

Principal Component Analysis (PCA) has been called one of the most valuable results from applied linear algebra [5]. PCA is used abundantly in several forms of analysis, because it is a simple, non-parametric method of extracting relevant information from confusing data sets. With minimal additional effort PCA provides a roadmap for how to reduce a complex data set to a lower dimension to reveal the sometimes hidden, simplified structure that often underlie it [8]. It is a way of identifying patterns in data, and expressing the data in such a way as to highlight their similarities and differences. Since patterns in data can be hard to find in data of high dimension, where the luxury of graphical representation is not available, PCA is a powerful tool for analyzing data.

The purpose of PCA is to identify the structure behind a multivariable stochastic observation in order to obtain a compact description of it. Clusters are more likely to be distinguished by projection in a high-variance direction than in a low variance one. Hence, the dimensionality reduction in the PCA could retain most of the intrinsic information in the data [9]. In PCA, the $k$th principal component direction is along an eigen vector direction belonging to the $k$th largest eigen value of the covariance matrix of the input data. An advantage in using PCA to preprocessing data is a reduction in data dimension that eliminates redundant information and allows simplifying the classification process [10].

2.4. Health assessment based on the Gaussian mixed model

This paper regards health assessment as evaluating the overlap between the current observed features and those observed during the normal condition. This overlap adopts Confidence Values (CVs) as an index for evaluating the severity of a health condition. This index was introduced by Jay Lee [9, 23] and has been accepted by some researchers [24-26]. The Confidence Value lies between 0 and 1, with 1 being the normal state. In this paper, the GMM is employed to assess the health state by using only normal condition datasets.

The GMM, as shown in Fig. 1, is a statistical model that uses a weighted sum of probability density functions of multiple Gaussian distributions to depict the distribution of a vector in the probability space [24]. This method can theoretically be applied to approximate an arbitrary feature distribution within an arbitrary accuracy [23].

Fig. 1Structure of a GMM with M mixture components

The mathematic model of GMM is parameterized by mixture weights $\alpha$, mean vectors $\overrightarrow{\mu }$, and the covariance $\sum$ from all $M$ mixture components as follows:

GMM can be given as:

where $p\left(x\right|{\theta }_i)$ is a single Gaussian function, $x$ is the $d$­dimensional feature vector, $M$ is the number of mixtures, and ${\theta }_i$ represents the parameters of the $i$-th Gaussian function, including mean vector ${\overrightarrow{\mu }}_i$ and covariance matrix ${\sigma }_i$.

For health assessment, given the training feature vectors from normal machine conditions, the goal of model training is to estimate the parameters of the GMM built under normal machine conditions [24]. Several methods can be used to estimate these parameters. Maximum likelihood (ML) estimation, which is the most popular and well-established method, is used in the current work. The ML parameters can be estimated iteratively by using the expectation maximization algorithm. The main steps as follows [25-28]:

1) (E step) Using Bayes’ theorem, calculate the posterior probability of model $i$, given a data vector $x_n$ and current model parameter $\theta$:

2) (M step) Maximum likelihood (ML) re­estimation of model parameters is given by:

Repeating the E and M steps, and the calculation converges to a stable solution that represents the maximum likelihood solution of the problem, thus getting converged mean, covariance matrix and prior vectors.

When using GMM to approximate the feature distribution, how to properly choose the number of the mixture is critical, or the problems of underfitting oroverfitting will arise. Bayesian Information Criterion (BIC) is introduced to address this issue [26].

BIC defines a log­likelihood function and penalty term, and its mathematical description given GMM is depicted as follows:

where $H_j$ is the $j$th candidate model, $D$ is the training feature, $In\left(P\right(D|{\widehat{\theta }}_j,H_j))$ is the maximized log­likelihood of the $j$th candidate model, $k$ is the number of the parameters to be estimated and $n$ is the size of the feature. The candidate GMM with the lowest BIC score can be deemed as the best model.

After establishing the GMMs for normal condition and monitoring condition separately, the Confidence Value ($CV$) can be calculated by the overlap of the GMMs as follows:

where $g_1\left(x\right)$ and $g_2\left(x\right)$ are GMMs that represent the normal and current monitoring conditions, respectively. This paper uses only normal condition datasets for overlap calculation, where each overlap value between the clustered normal dataset and online measurement data can be considered Confidence Value ($CV$), as well as for tracking health degradation status.

3.1. Experimental setup

In this paper, two datasets from different bearing test rigs were used to verify the proposed method. The first dataset was generated from a test rig in the Case Western Reserve University Bearing Data Center (http://csegroups.case.edu/bearingdatacenter/pages/download-data-file). In this experiment, as Table 1 shows, 6205-2RS JEM SKF deep-groove ball bearings were tested. Three common types of bearing faults were tested: inner race fault, outer race fault, and rolling element fault. Under four states, including the normal condition, the vibration signal was acquired with a stabilized motor speed of 1797 r/min and a sampling rate of 12 kHz.

The second dataset was generated from bearing run-to-failure tests performed under a constant load condition by the National Science Foundation Industry and University Cooperative Research Program (NSF I/UCR） Center on Intelligent Maintenance Systems (IMS) [29]. As Table 2 shows, Rexnord ZA-2115 double row bearing was selected in the run-to-failure test. The rotation speed was kept constant at 2000 rpm, and 6000 lbs (2730 kg) of radial load was placed onto the shaft and bearing via a spring mechanism [29]. Vibration data were collected every 20 min at a sampling rate of 20 kHz. The test was continued until an inner race fault in the bearing was found. Both datasets were analyzed to verify the proposed method.

3.2. Feature extraction

Firstly, the vibration signal from two datasets was decomposed into $n$-empirical modes that contain lower-frequency oscillations than the previously extracted ones and a residue $r_n\left(t\right)$, as shown in Fig. 2. Information on bearing health condition is considered in the high-frequency bands. Therefore, in this paper, $n$ was designated as 7. The first seven IMFs and a residue $r_7\left(t\right)$ were extracted to construct the feature matrix for the following SVD process. Then, the constructed feature matrix was decomposed to obtain an 8-dimensional feature vector of singular values. After that, PCA was employed to reduce a complex data set to a lower dimension to reveal the hidden, simplified information. In this paper, 8 dimensions were reduced to 2 dimensions.

Fig. 2A sample of EMD decomposition results of normal condition datasets

3.3. Health assessment by using only normal condition data

To verify this proposed method, both datasets from different test rigs were used. From the first dataset, a sample of normal condition data were employed to construct $g_1\left(x\right)$ in Eq. (17) and train the GMM model. Then, another 4 samples of normal condition data and 2 samples of data acquired under each fault mode were used to test the trained GMM model.

As shown in Table 3, Confidence Values (CVs) of normal condition from the first dataset were near 1, while Confidence Values (CVs) of all fault modes were close to 0. This result verified that the GMM is very notable for health assessment.

From the second dataset, with the first 400 min collected datasets considered normal condition data, the Confidence Values (CVs) from the beginning to the end of the life testing of the bearing were calculated by using Eq. (17), as shown in Fig. 3.

Fig. 3Confidence Values obtained from the second dataset

To explain the advantages of the proposed method, we compared it with the result of the reference [29], which employed the Self Organizing Map (SOM) to process the same dataset for performance degradation assessment. The result is shown in Fig. 4. It is found that reference [29] uses the value of minimum quantization error (MQE) to reflect the health state. The larger the value of MQE, the worse the health state is. If there is no fault data, it is unavaiable to know how large the MQE is means the failure appearance. In this situation, the method of reference [29] is useless for condition-based predictive maintenance. Meanwhile, the proposed method in this paper can be used without fault datasets in real applications, since the Confidence Values (CVs) gradually tend to 0. Besides, by comparison, it is also implied that the proposed method shows more sensitive to small initial changes in health conditions, and can detect defects at the incipient stage.

Fig. 4. Health assessment results from reference [29]

The blue curve in Fig. 3 represents the detected health degradation. In the first 900 min, the bearing was in good health and considered to be in the “normal stage”. From approximately 900 min to 5400 min, initial defects appeared and propagated, which decreased the Confidence Values (CVs); this is considered the “defect stage”. From approximately 5400 min to 6400 min, the defects worsened considerably, and the Confidence Values (CVs) decreased rapidly; this condition is considered the “danger stage”. A few minutes after 6400 min, fault occurred, and this condition is considered the “failure stage”. In Fig. 3, several thresholds of 0.8, 0.4, and 0.1 were selected. These thresholds were used to determine the health scale of the bearing, which is significant for condition-based predictive maintenance.