Rolling Bearing Reliability Assessment via Kernel Principal Component Analysis and Weibull Proportional Hazard Model

Reliability assessment is a critical consideration in equipment engineering project. Successful reliability assessment, which is dependent on selecting features that accurately reflect performance degradation as the inputs of the assessment model, allows for the proactive maintenance of equipment. In this paper, a novel method based on kernel principal component analysis (KPCA) and Weibull proportional hazards model (WPHM) is proposed to assess the reliability of rolling bearings. A high relative feature set is constructed by selecting the effective features through extracting the time domain, frequency domain, and time-frequency domain features over the bearing’s life cycle data. The kernel principal components which can accurately reflect the performance degradation process are obtained by KPCA and then input as the covariates of WPHM to assess the reliability. An example was conducted to validate the proposed method. The differences in manufacturing, installation, and working conditions of the same type of bearings during reliability assessment are reduced after extracting relative features, which enhances the practicability and stability of the proposed method.


Introduction
The rolling bearing is one of the most important components of rotating machinery [1], and its running state directly influences the health of the entire system of equipment [2].It is also very easily damageable, as it not only supports load but also permits relative motion [3,4].Effective rolling bearing maintenance strategies can not only reduce the amount of downtime and cost of maintenance, but also ensure the normal operation of the equipment [5,6].Accurate reliability assessment is essential for making predictive maintenance decisions based on the real-time status of equipment.
Reliability assessment comes with two key challenges: the construction of an appropriate reliability assessment model and the selection of features which can accurately reflect performance degradation process.Reliability assessment based on real-time equipment conditions has become a popular research topic in recent years [7], and advancements in information technology and artificial intelligence have brought about a number of valuable contributions to the literature.For example, the proportional hazard model (PHM), first introduced by the British statistician Cox [8] in 1972, is a powerful statistical analysis methodology.It is an important statistical regression model based on lifetime data and condition monitoring data and has been successfully used for reliability assessment in accelerated life testing.Ding et al. [9] used the Weibull proportional hazard model (WPHM) to assess the reliability of a rolling bearing in real time.Liao et al. [10] used logistic regression model and PHM to assess the reliability of an individual unit.Zhang et al. [11] used a mixed WPHM to predict the failure of a mechanical system with multiple failure modes.
The WPHM is a well-established mathematical model.However, when it is applied in real equipment life prediction, it is problematic as far as covariates selection, setting reliability threshold, trend prediction, and other issues.In terms of covariates selection, most previous studies concern direct time domain statistical analysis where one or more time 2 Shock and Vibration domain features are selected to build a reliability assessment model.However, one single feature or features based on one single domain cannot accurately reflect the performance degradation process and thus impact the overall accuracy of reliability assessment.Although time domain, frequency domain, and time-frequency domain features can comprehensively reflect the performance degradation process of bearings over their entire service lifetime, excessive parameters lead to data redundancy.Further, selecting more covariates of the WPHM makes the parameter estimation process more challenging.The vibration signals of faulty machinery are generally nonstationary and nonlinear under complicated operating conditions [12,13].Therefore, it is crucial to select the features through nonlinearly reducing dimensionality and removing redundant features.
Kernel principal component analysis (KPCA), first proposed by Schölkopf et al. [14,15], is generalized principal component analysis (PCA) that is applied to nonlinear cases by nonlinearly mapping input samples to a higher dimensional feature space before performing PCA per usual [16].It has been successfully utilized in process monitoring and fault diagnosis applications.Lee et al. [17] developed a new nonlinear process monitoring technique based on KPCA.Jiang et al. [18] proposed a fault diagnosis approach based on KPCA and multiclass classifiers of a support vector machine.He et al. [19] used the low-dimensional principal component representations from the statistical features of measured signals to characterize and monitor gearbox conditions.Su et al. [20] used a Euclidian distance discriminating approach to distinguish bearing fault data by adopting the first seven principal components as inputs.
In order to overcome the weakness for the selection of WPHM covariates, this paper proposes a novel method for assessing the reliability of rolling bearings based on KPCA and WPHM.The novelty of this research is in improving the covariates selection method of WPHM, which has considerable value in practical application.A high relative feature set is constructed by selecting the effective features through extracting the time domain, frequency domain, and timefrequency domain features over the bearing's life cycle data.Then the first three kernel principal components (KPCs), which can accurately reflect the performance degradation process through KPCA, are selected as WPHM covariates to assess the reliability.The feasibility and effectiveness of the method were validated using bearing's life cycle data, and it can provide important basis for equipment proactive maintenance.The differences in manufacturing, installation, and working conditions of the same type of bearings during reliability assessment are reduced after extracting relative features, which enhances the practicability and stability of the proposed method compared to traditional assessment techniques.It enriches the theory of covariates selection and is more emphasis on application innovation.
The remainder of this paper is organized as follows.Section 2 presents the fundamental theories of KPCA and WPHM.Section 3 presents the proposed method for reliability assessment in detail.Section 4 discusses the features extraction method used to reflect the bearing performance degradation process through KPCA.The case study we conducted to validate the proposed method reported in Section 5, and conclusions are given in Section 6.

Fundamental Theory
2.1.Kernel Principal Component Analysis.KPCA essentially works by nonlinearly mapping input samples to a highly dimensional feature space F and applying a linear PCA to the transformed signals.KPCA performs nonlinear data processing more effectively than PCA.
In KPCA, a set of multidimensional signals x  ,  = 1, 2, . . ., , is mapped to Φ(x  ),  = 1, 2, . . ., , by nonlinear mapping Φ : R → F. Assume Φ(x  ),  = 1, 2, . . ., , has been mean-centered.PCA is performed by finding the eigenvalues  > 0 and eigenvectors a ∈ F satisfying V = S  V, where the sample covariance matrix of Φ(x  ) is Substituting (1) into the eigenvector equation yields The eigenvectors can be expanded as follows: where   is correlation coefficient.Substituting (3) into (2) yields and K is symmetric matrix, where Equation ( 4) can then be written as a = Ka. 1 ≥  2 ≥ ⋅ ⋅ ⋅ ≥   ≥ ⋅ ⋅ ⋅ ≥   are the corresponding eigenvalues of K, and a 1 , a 2 , . . ., a  , . . ., a  are the eigenvectors of K.If   is the minimum eigenvalue (a nonzero number), normalized eigenvectors can be obtained successfully: Finally, the principal components for testing examples {y 1 , y 2 , . . ., y  } can be calculated as follows: Shock and Vibration 3 where The above algorithm is based on the assumption that Φ(x  ) is mean-centered, but the assumption is suitable in general.The mapping data to be mean-centered is expressed as follows: where l  is the  ×  unit matrix of the coefficient which is 1/.The cumulative contribution rate (CCR) is utilized to determine number  of principal components.
The cumulative contribution rate (CCR) threshold is referenced by [21,22].This threshold could be set at 85%, 90%, or 95%.In general, once the CCR exceeds 85%, the first principal components contain most information of the original feature set, so this paper is set at 85%.
There are three common types of kernel functions: polynomial kernel function, radial basis Gaussian (RBG) kernel function, and neural network kernel function.The transformation matrices of the RBG kernel function have positive definiteness and a wide convergence field.It only contains one parameter, and the calculation process is relatively simple [23].The RBG kernel function (x, y) = exp{−‖x − y‖ 2 /2 2 } is utilized here, where  2 is a parameter related to the kernel width. 2 can be obtained by optimizing the parameter of the kernel function [24].

Weibull Proportional Hazard Model.
The PHM builds a mathematical relationship between the feature parameters of the equipment running status and the reliability.According to the feature parameters of the real-time operation, PHM can get the device hazard rate in its current state, to assess the current reliability of the equipment.The hazard rate at time  is expressed as follows: where ℎ 0 () is the baseline hazard rate dependent on the service time, z  is a row vector composed of monitoring values at time  that is time-dependent, and  is a column vector composed of the regression parameters corresponding to the monitoring variables.In the PHM, z  is regarded as a vector of covariates that increases or decreases the system hazard rate proportionally; its coefficient vector  defines the influence of the monitoring variables on the failure process.The Weibull distribution is frequently used to model the failure time of mechanical systems.The hazard rate function of the Weibull distribution is commonly selected as the baseline hazard rate of the PHM.The hazard rate for the twoparameter Weibull distribution is written as follows: where  > 0 and  > 0 are the shape and scale parameter of the Weibull distribution, respectively.The PHM with the Weibull baseline function is called the WPHM, the hazard function of which is defined as follows: According to the principle of reliability analysis [25], reliability and failure probability density can be estimated as follows: The key of using WPHM to assess the operating status of equipment is to estimate unknown parameters according to the feature data and time data of the real-time status.The maximum likelihood method is commonly applied to estimate unknown WPHM parameters.In practice, a mechanical system may be run until it fails but may be repaired prior to failure.The lifetime data usually contains failure times and suspension times to reflect this.To properly account for both types of data, the likelihood function where the covariates are time-dependent is defined as follows: where  indexes the failure times,  indexes the suspension times,  is the number of failure samples, and  is the number of suspension samples.By substituting ( 14) into ( 15), the likelihood function can be rewritten as follows: where  indexes both the failure times and the suspension times.The log-likelihood function is In the above equations, the covariates of WPHM are timedependent.When the covariates only relate to the current time (i.e., they are non-time-dependent), the reliability and the failure probability density can be, respectively, rewritten as follows: Therefore, by substituting ( 18) into ( 17), the log-likelihood function can be rewritten as follows: By setting the partial derivatives of ( 17) or ( 19) with respect to the parameters , , and  equal to zero, β, η, and γ can be obtained via Newton iterative method.Unless the initial value is suitable, the numerical solution is very difficult to obtain.
With increase in the number of covariates, the complexity of the maximum likelihood estimation increases substantially.Therefore, the Nelder-Mead iterative algorithm [26,27] is applied to estimate these mixed parameters.

Proposed Method
A flowchart of the proposed method is shown in Figure 1.
The reliability assessment process takes place in a stepwise manner: (1) Select effective features that comprehensively reflect the performance degradation process from the time domain, frequency domain, and time-frequency domain features of training bearings data to compose the feature vector.
(2) Build a high relative feature set by extracting samples from lifetime data of training bearings.
(3) Obtain KPCs and KPCs mapping from KPCA for the high training lifetime relative feature set, and select the first KPCs with CCR exceeding 85%.
(4) Build the high relative feature set from lifetime data of the test bearing; obtain KPCs of the test bearing through KPC mapping of the training bearings, and then verify whether the KPCs of the test bearing reflect the performance degradation process.(5) Take the KPCs of the training bearings as the WPHM covariates to estimate the WPHM parameters.(6) Take the KPCs of the test bearing as the WPHM covariates to assess the test bearing reliability.

Performance Degradation Process of Rolling Bearing
4.1.Experimental Setup.The rolling bearing life cycle test data used in this paper was provided by the Center for Intelligent Maintenance Systems (IMS), University of Cincinnati [28].The experimental data sets were generated from bearing run-to-failure tests under constant load conditions on the specially designed test rig as shown in Figure 2.
There are four test double row bearings (Rexnord ZA-2115) on one shaft of the bearing test rig.The shaft is driven by an AC motor and coupled with rub belts.A radial load of 6,000 lbs was added to the shaft and bearings by a spring mechanism; the rotation speed was kept constant at 2,000 rpm during the experiment.A magnetic plug is installed in the oil feedback pipe to collect debris from the oil as evidence of bearing degradation.The test was stopped when the accumulated debris adhered to the magnetic plug exceeded a certain level, causing an electrical switch to close.Vibration data were collected every 20 minutes with a National Instruments DAQCard-6062E data acquisition card (data sampling rate 20 kHz and data length 20,480 points).Data collection was conducted in the National Instruments LabVIEW program.Table 1 summarizes the test results.Figure 3 shows the components of the failure bearing.
Bearing 3 (Test 1) is used as test bearing and the other bearings (Test 1 and Test 2) are used as training bearings.

Constructing
High Relative Feature Set.More than 50 features of time domain, frequency domain, and time-frequency domain [29,30] were extracted from the life cycle data of seven training bearings.To ensure that the first three KPCs contain as much useful information as possible, it is necessary to minimize the dimension and ensure the validity of the information of each dimension prior to using KPCA.First, each feature was drawn to value-time figures by lifetime data; next, the features that did not reflect the degradation process (e.g., mean, skewness) were eliminated; the similar features with poor performance degradation process were then removed by comparison, such as empirical mode decomposition (EMD) normalized energy spectrum; finally, 11 features that can comprehensively reflect different degradation stage were selected through the comparison.They are as follows: (1) Time domain: root mean square (RMS), kurtosis, peak-peak value, and peak factor.(2) Frequency domain: spectrum mean, spectrum variance, and spectrum RMS.
Even the same type of bearings differ due to differences in manufacturing, installation, and working conditions; thus there are differences among them even in the same work period.Take time domain features as an example; for bearings 1-8, the time domain features for stable trend of normal work period were selected and averaged as shown in Figure 4.
Figure 4 shows where there was a sizable difference in time domain features from the stable trend of all eight bearings.As shown in Figure 4(a), the average RMS of bearing 1 in the normal working period was 0.154 while that of bearing 5 was 0.077.The features should be standardized to reduce their influence on the assessment results.First, the features Shock and Vibration   for stable trend of normal work period are averaged as the standard values; then the ratio of original features to the standard values is calculated to obtain relative features (standardized features).Relative features can be defined as follows: where () is the relative feature, () is the original feature, and  is the mean value of normal work period.The advantages of relative features are discussed further in Section 4.3.
For the seven training bearings, each had 100 samples that can reflect the process of the lifetime (at a total of 700 samples).A high training lifetime relative feature set of 700 × 11 was then composed accordingly.For the test bearing, a total of 2,152 whole life cycle samples were used to obtain the 2152 × 11 high test lifetime relative feature set shown in Figure 5, where In order to effectively verify the following analysis, according to the features of bearing 3, the degradation process was divided into five stages as shown in Table 2. To show the   features of each degradation stage, Figure 5 was redrawn as Figure 6. Figure 5 shows the overall trend and details of the bearing degradation process, while Figure 6 shows the differences in the features of the bearing degradation process at different stages.Many previous researchers have used RMS and kurtosis as PHM covariates.Figure 6(a) shows that although RMS can distinguish the normal working stage, early failure stage, medium wear state, and severe wear stage, it is merely a time domain feature and contains less information than principal components.Figure 6(c) shows that kurtosis is only sensitive to early failure.It fluctuates greatly with time and varies in a large range in the early failure stage.Moreover, kurtosis in medium wear stage and severe wear stage is slightly larger than the normal working stage.Therefore, kurtosis, which is unable to reflect the performance degradation process, is suitable as an important indicator of early failure warning rather than as a covariate.Most importantly, each feature is inevitably subject to large or small fluctuations in the whole life cycle.In this paper, the KPCA method can effectively solve the above problems, and the details are shown in the following section.3.
The CCR of the first three principal components of high training original feature set was lower than high training relative feature set.In other words, the dimension reduction effect of high training original feature set was lower than high training relative feature set.This fully verifies the advantages of relative features.Hence high training relative feature set was selected for subsequent analysis.
The first three principal components of high test relative feature set can be obtained by KPCA.To verify the result of KPCA, the first three KPCs were projected onto threedimensional coordinate system as shown in Figure 7(a), and the first two KPCs were projected onto two-dimensional coordinate system as shown in Figure 7(b).
As shown in Table 3, the contribution rate of the first kernel principal component (KPC1) reached 63.24%, suggesting that it contains most of the high test relative feature set information.Figure 7 shows that the first principal component well reflects the degradation trend of test bearing.The second principal component (KPC2) and third principal component (KPC3) contribution rates were 17.47% and 8.49%, respectively, indicating that they contained little information.As time progressed, the data points formed an obvious trend direction with a fairly smooth change process in general.The   major degradation trend is KPC1 positive direction, KPC2 negative direction, and KPC3 negative direction.The normal period was concentrated within a small area.The early fault period deviated from the main direction and developed into KPC1 positive direction, KPC2 positive direction, and KPC3 positive direction.The early fault period was easy to distinguish from other periods.The recovery period reverted to the main direction.The medium wear stage and severe wear stage continued to rapidly develop along the main direction.Because we choose the kurtosis value which is sensitive in the early failure stage, the kernel principal component degradation trend direction was offset in early failure stage.
By comparison, the first three principal components were projected as shown in Figure 8(a), and the first two KPCs are projected onto two-dimensional coordinate system as shown in Figure 8(b).The contribution rate of the first principal component (PC1) was 63.85%, while the early fault data point distribution range was large enough to overlap with the medium wear stage and severe wear stage.The stages were difficult to distinguish from each other in Figure 8, and trend is rougher and with larger offset compared to Figure 7.The bearing performance degradation process overall trend is more obvious in Figure 7  , so the offset is large (Figure 8).KPCA, which exploits nonlinear analysis, results in a well-offset early fault stage in Figure 7 while the change trends are smoother and offset is smaller than Figure 8; in the medium and severe wear stages (33.49-34.17/d), Figure 7 data points show a clearer change trend overall than in Figure 8.
In conclusion, the first three KPCs contain most time domain, frequency domain, and time-frequency domain information as well as nonlinear components.Between the sample points, the offset is relatively small and the trend is obvious.These first three KPCs, which fully represent the bearing performance degradation process, can be utilized as covariates to establish a WPHM model that is highly stable and reliable.

Bearing Reliability Assessment
The first three KPCs were selected as the covariates of WPHM to assess the reliability.By substituting the failure and suspension data of KPCs from the high training relative feature set into ( 17) and ( 19), β, η, and γ can be obtained as shown in Table 4.The KPCs of the high test relative feature set are plugged into the WPHM with non-time-dependent covariates to calculate the reliability, as shown in Figure 9.
The reliability of the normal working stage remained stable between 0.99 and 0.90; the reliability of the early failure stage fluctuated between 0.90 and 0.75; the reliability of the healing stage fluctuated between 0.87 and 0.84; the reliability of the medium wear state gradually fell to 0.50 from 0.75; and the reliability of the severe wear stage fell rapidly from 0.50 to 0.45.The reliability variable accurately reflects the state of the test bearing.When reliability drops to 0.90, the bearing requires attention.When reliability drops to 0.75, it urgently requires attention, and a maintenance plan should be enacted immediately.When the reliability falls to 0.50, the equipment must be stopped to avoid an accident.
By contrast, The KPCs of the high test relative feature set were plugged into the WPHM with time-dependent covariates to calculate the reliability as shown in Figure 10.
The reliability decreased at a generally steady rate in the normal working stage and then began to decrease at a quicker pace in the early failure stage.The reliability decline rate increased sharply in the medium wear stage and severe wear stage.Therefore, the reliability assessment of WPHM with time-dependent covariates can reflect not only the degradation process of the bearing, but also the stage in its life cycle.
Considering the influence of the historical data, the reliability assessment in WPHM with time-dependent covariates is highly stable and credible.Therefore, it is better suited to maintaining important equipment.However, reliability as reflected in WPHM with non-time-dependent covariates is only related to the present time without considering the historical data; the complexity of calculation decreases, so it is more suitable for normal equipment.More importantly, it is more suitable for the reliability assessment of the bearings which lack historical data or have been repaired.Time-dependent covariates or non-time-dependent covariates can be selected according to actual conditions.Whether the covariates are time-dependent can be determined by the actual needs of the WPHM.In conclusion, the results indicate that this method can accurately assess reliability and timely provide effective maintenance decisions.

Conclusion
In this study, KPCs based on KPCA were successfully used as WPHM covariates to assess the reliability of rolling bearings.Based on the relative multiple features, KPCs can sufficiently describe the bearing performance degradation process.KPCs as WPHM covariates provide accurate reliability to support timely maintenance decisions.The relative features also enhance the practicability and stability of the proposed method compared to traditional assessment techniques.

Figure 1 :
Figure 1: The flowchart of the proposed method.
Time domain-peak factor

Figure 4 :
Figure 4: Mean time domain features in normal work period.
(a)-(d) are the time domain features; (e)-(g) are the frequency domain features; (h)-(i) are the normalized energy spectrum of wavelet packet 3,7 bands; (j)-(k) are the sample entropy of wavelet packet 3,7 bands; and (l) denotes whether the vibration data was collected in the given test period.The discontinuous points indicate a lack of data collection.
The time of data collection

Figure 5 :
Figure 5: High test lifetime relative feature data (overall degradation process).
The time of data collection

Figure 6 :
Figure 6: High test lifetime relative feature data (individual degradation stages).

Figure 7 :Figure 8 :
Figure 7: Kernel principal component projection: (a) the first three Kernel principal components; (b) 1st kernel principal component and 2nd Kernel principal component.
as well.For example, in the early fault stage (31.05-33.34/d),changes in peak-to-peak value, kurtosis value, and peak factor change range are large and fluctuate with time (Figures 5(b)-5(d))