Failure Threshold Determination of Rolling Element Bearings Using Vibration Fluctuation and Failure Modes

: One of the challenges in predicting the remaining useful life (RUL) of rolling element bearings (REBs) is determining a proper failure threshold (FT). In the literature, the FT is usually assumed to be a constant value of an extracted feature from the vibration signals. In this study, a degradation indicator was extracted to describe damage to REBs by applying principal component analysis (PCA) to their run-to-failure data. The relationship between this degradation indicator and the vibration peak was represented through a joint probability distribution using statistical copula models. The FT was proposed as a probability distribution based on the ﬂuctuation increase in the vibration trend. A set of run-to-failure tests was conducted. Applying the proposed method to this dataset led to various FTs for the different failure modes that occurred. It is shown that, for inner race degradation, a higher FT can be assumed than for rolling element degradation. This could help extend the lives of REBs regarding the degrading elements. A dataset for an industrial machine was also analyzed and it is shown that the proposed model estimated a reasonable and proper FT in an actual case study.


Introduction
Rolling element bearings (REBs) are one of the most common mechanical components of rotary machines. The failure of these components is the main reason for 45% to 55% of machines breakdown [1]. This is why researchers have paid more attention to predicting the remaining useful life (RUL) of REBs over the past two decades. The RUL of REBs is defined as the remaining time until a vibration feature reaches the predetermined failure threshold (FT) [2], which is often assumed to be a constant value in the literature. Therefore, the accuracy of the RUL prediction depends on selecting a proper FT. Furthermore, the determination of FTs for various vibration features is a challenging issue. In general, to achieve protective or predictive goals in preventing the breakdown of rotary machines, the FT is determined experimentally, based on the manufacturers recommendation or using the recommendations of standard codes (e.g., ISO 10816-3 [3]). The recommendations in ISO 10816-3 and other similar codes correspond to different fault types in machines in the frequency range of 10 to 1000 Hz. They are thus not specifically for the failure of REBs. Therefore, employing a suitable and proper FT for REBs leads to the more reliable and economical operation of rotary machines.
Reviewing the literature, three different approaches can be considered for defining the FT: (1) as a constant value of a vibration feature, (2) as a probability distribution of a vibration feature or (3) as a specific level of physical damage. In most research studies on the RUL prediction of REBs, the FT is defined as a constant value of a vibration feature (e.g., the root mean square (RMS) or the time signal entropy) [4][5][6][7][8]. This approach is the most convenient method for this issue. For instance, in the PRONOSTIA dataset [9], which is a set of run-to-failure tests of REBs conducted by the FEMTO-ST institute in 2011, the constant value of the vibration peak (20 g) was employed as the stop criterion of the tests. After the publication of this dataset at the PHM 2012 conference, many researchers used it to validate their studies [10][11][12]. Most of them have used the test stop criterion as the FT in their RUL prediction problems.
In the second approach, the FT is defined as a probability distribution of a vibration feature. This approach has been widely used in reliability literature. In this regard, Wang and Coit [13] studied the effect of different FT distributions on equipment reliability. Nystad et al. [14] used the gamma distribution of the degradation indicator for the FT, the parameters of which were determined by experts, to estimate the probability of chock valves' RUL accurately. In another study, Yu and Fuh [15] studied the failure of an element in the crack propagation problem. They proposed a probability distribution of the crack size to define a failure threshold in the fatigue test results.
Using a specific level of physical damage as a FT is less common than the previous approaches. Another dataset that has been employed in the RUL prediction literature is that of the REB run-to-failure tests conducted by the Center for Intelligent Maintenance Systems (IMS) at the University of Cincinnati [16]. In this dataset, the test stopping criterion was defined as a specific contaminant level in the REB lubrication feedback system. Reaching this level led to the activation of an electrical switch that stopped the test running. The dataset included three run-to-failure experiments on four similar parallel REBs mounted on a single shaft. These experiments stopped with the same physical damage criterion, which was explained earlier. However, investigating the results of the vibration records revealed a significant difference between the levels of various vibration features at the end of these three tests. This means that the stopping criterion in this dataset employed led to extremely different results from the first approach, which was explained earlier. Some researchers [17,18] have used the stopping criteria of these tests as the FT in their studies on the IMS dataset. However, others [19,20] have employed a constant value of the vibration features as the FT in IMS tests. Therefore, the end of life (EoL) may differ from the end of the test in these studies.
Reviewing the literature, it can be seen that the details of defining a proper FT were less extensively considered in the RUL prediction of REBs. As the second class of approaches mentioned in the previous paragraphs deal with the uncertainty of the FT, a proposal of a probability distribution for the FT is accordingly desired. In this paper, a procedure to define the FT for REBs is proposed. The extracted feature from the principal component analysis (PCA) is employed to describe the defect growth in the elements of REBs. The copula method is used to relate the fluctuation of the vibration features to the defect growth and selection of a proper FT. The procedure proposed in this paper is applied to experimental run-to-failure data, including accelerated life tests in the lab and the actual field data of an industrial machine. For the accelerated life test dataset, three different scenarios for types of failure that could occur in the tests are planned and studied. The results show that different FTs can be used for the various probable failure modes in REBs. Furthermore, applying the proposed method to the actual field data shows that it provides a proper FT suggestion in the industrial case.

Proposing a New Failure Threshold
In this section, different features extracted from the vibration signals are studied. Then, by using PCA, the first principal component is extracted as a degradation feature. Copula models are used to investigate the statistical relationship between the first principal component and the other popular features. The new proposed criterion for the FT is introduced based on this relationship. This proposed criterion will be applied to the experimental test results in the next sections.

Feature Extraction
In this paper, fifteen popular features (from the time and frequency domains) [21,22] are employed to study the degradation of REBs. These features and their formulas are introduced in Table 1. Each of these features can describe degradation in different ways and they may reveal various points of information about the hidden, growing defects in a degrading REB. Therefore, all of them are considered for analysis to extract a proper degradation indicator in the next subsection. Table 1. Definition of signal features in time and frequency domains.

Feature Definition Feature Definition
Peak amplitude RMS(x(n)) 6 Peak-to-peak max(x(n)) − min(x(n)) Margin factor Root mean square frequency

PCA Algorithm
The features introduced in the previous subsection were employed as the inputs of the PCA algorithm. This algorithm uses the matrix characteristics to create a new feature that describes the maximum data variance. This feature is known as the first principal component (PC 1 ). The PC 1 can be used as a degradation indicator in prognostic problems. Some researchers have assumed a linear relationship between PC 1 and the size of the growing defect [23].
The PCA algorithm is an orthogonal linear transform used to reduce the dimensions of data. This algorithm creates new features based on their importance (eigenvalues) so that the primary features have more data variance and better illustrate the nature of data. Assume X n × L = [x 1 x 2 x 3... x L ] represents the input data, where L is the number of observations in an n-dimensional space of which n is the number of rows (features). The PCA algorithm converts X data intoX m × L = [x 1x2x3...xL ] via Equation (1): where m is smaller than n and W T m × n is the transformation matrix. The algorithm then uses eigenvalues and eigenvectors, as in Equation (2), to calculate matrix W, where λ is the eigenvalue, V i is the eigenvector and C x is the covariance matrix of the X matrix that is calculated through Equation (3). As a result, the covariance matrix of the new components (new features) is diagonal and these components lack any correlation. Finally, the PCA algorithm outputs the eigenvalues in descending order (Equation (4)). The first component of the new features (PC 1 ) that corresponds to the largest eigenvalue explains the most data variance and is as assumed to be the best indicator of data. In the rest of this paper, this component is called the degradation indicator.

Copula Statistical Models
In the next step, Pearson's correlation coefficient (PCC) was applied to analyze the correlation between PC 1 and the time domain features. The feature with the highest correlation was selected for the next step of the analysis. Then, copula models were employed to indicate the statistical relationship between PC 1 and the selected feature. Hence, using the copula model to create a joint probability distribution, it is possible to obtain a conditional probability distribution for the selected feature and PC 1 . In other words, the probability distribution of the time feature can be determined for different defect sizes.
The copula models are mathematical functions used to create a joint probability distribution of two random variables, such as X and Y, in the form H (X, Y) = P (X ≤ x, Y ≤ y). Sklar's theory [24] claims that if H is a cumulative joint distribution function and F and G are univariate marginal distributions, there is a C copula model where H (X, Y) = C (F(X), G(Y)). This copula is unique if F and G are continuous. Similarly, if H is joint distribution, u = F(X) and v = G(Y), the copula model is obtained by equation C (u, v) = H (F −1 (u), G −1 (v)). One of the most common copula families in the literature is the Gaussian copula, which is defined as follows: where Φ Σ is the cumulative joint normal distribution function with a mean of zero, Σ is the covariance matrix and φ −1 is the inverse indicator of the standard Gaussian distribution. Moreover, the conditional probability distribution is determined using the copula model as follows: In this study, eight parametric models of copula families, with details listed in Table 2, were employed [25]. Moreover, the model parameters were updated using the Bayesian method. The root mean square error (RMSE) and Nash-Sutcliffe efficiency (NSE) criteria were also used to assess the goodness of model fit.
In an ideal model, we would have RMSE = 0 and NSE = 1. In these equations, y i and ∼ y i represent the model values and the observation values, respectively. Table 2. Mathematical description of copulas [25].

Mathematical Description Parameter Range
Gaussian

Proposing a New Criterion for the FT
REBs are usually the weak and failure-prone elements of rotary machines. Various mechanisms may lead to the failure of these elements [26]. If an REB is designed accurately, adequately lubricated and installed, continuously aligned during operation and is kept away from moisture, corrosive materials and excessive loading, the only mechanism that can lead to the failure of the REB is rolling contact fatigue (RCF) [27]. The failure caused by RCF can include soft failure and hard failure stages. In the soft failure stage, degradation occurs gradually over a long period of time and the REB does not entirely lose its efficiency. However, in the hard failure stage, the REB undergoes sudden impacts that can impose catastrophic damage on the equipment and the machine completely loses its function [4]. When the hard failure stage begins, fast growth is observed in the vibration level trend ( Figure 1). Guidelines highly recommended terminating the operation of a machine when its vibration level passes a specific limit. Manufacturers and some standard codes [3] suggest values for this limit based on the type and properties of the machines. This limit is employed to define a high-risk zone for the equipment. Moreover, the transformation of degradation from soft failure to hard failure corresponds to the changing vibrations from the stable stage to the unstable stage. The beginning the unstable stage can be interpreted as entering the high-risk zone. Therefore, if one can predict this transition point (from soft failure to hard failure), it can be used as the FT and define the end of life for the machine. Furthermore, instead of considering a static FT, predicting a probability distribution of the FT can enhance the reliability estimation. Guidelines highly recommended terminating the operation of a machine when its vibration level passes a specific limit. Manufacturers and some standard codes [3] suggest values for this limit based on the type and properties of the machines. This limit is employed to define a high-risk zone for the equipment. Moreover, the transformation of degradation from soft failure to hard failure corresponds to the changing vibrations from the stable stage to the unstable stage. The beginning the unstable stage can be interpreted as entering the high-risk zone. Therefore, if one can predict this transition point (from soft failure to hard failure), it can be used as the FT and define the end of life for the machine. Furthermore, instead of considering a static FT, predicting a probability distribution of the FT can enhance the reliability estimation.
There is no accurate and unanimous definition for the physical point of failure for a machine. Therefore, the recommendations of experts in the relevant field can be helpful in considering a limit for the high-risk region and decreasing failure possibility. In this paper, it is proposed that the high fluctuation occurrence in the vibration trend can be used as the FT in order to prevent machine operation in the high-risk zone. High fluctuation is assumed here to be an amount of 2 dB or ±25% variation of the vibration feature.
As mentioned in Section 2.2, PC 1 , which was proportional to the defect growth, was used as the degradation indicator. To provide a new suggestion for the FT, the relationship between vibration features and PC 1 must be studied. With this aim, a conditional probability distribution of the vibration feature given variable PC 1 was obtained through a proper copula model.
The conditional probability distribution function, in which 95% of the data are scattered in the range of ±25% around the mode of distribution, was used as the indicator of the initiation of high fluctuation in the vibration trend in this study. In other words, this distribution, H (F | PC 1 = pc 1 ), should satisfy the equation below: In Equation (9), F is the vibration feature most highly correlated with PC 1 , which was described in the previous subsections. The parameter f m is the mode of H (F | PC 1 = pc 1 ). Since H was assumed to correspond to high fluctuation initiation, it represents the FT probability distribution and f m is suggested as an FT.
The procedure proposed in this paper to estimate the FT of REBs is summarized in the flowchart depicted in Figure 2. First, different vibration features in the time and frequency domains were extracted from the run-to-failure history. These features were the inputs to the PCA algorithm that produced a feature as a degradation indicator (PC 1 ). Afterward, the vibration feature most highly correlated with PC 1 was selected. By using the copula models, a suitable joint probability distribution was fitted for the selected feature and PC 1 . Finally, the FT probability distribution was obtained based on the proposed criterion. There is no accurate and unanimous definition for the physical point of failure for a machine. Therefore, the recommendations of experts in the relevant field can be helpful in considering a limit for the high-risk region and decreasing failure possibility. In this paper, it is proposed that the high fluctuation occurrence in the vibration trend can be used as the FT in order to prevent machine operation in the high-risk zone. High fluctuation is assumed here to be an amount of 2 dB or ±25% variation of the vibration feature.
As mentioned in Section 2.2, PC1, which was proportional to the defect growth, was used as the degradation indicator. To provide a new suggestion for the FT, the relationship between vibration features and PC1 must be studied. With this aim, a conditional probability distribution of the vibration feature given variable PC1 was obtained through a proper copula model.
The conditional probability distribution function, in which 95% of the data are scattered in the range of ±25% around the mode of distribution, was used as the indicator of the initiation of high fluctuation in the vibration trend in this study. In other words, this distribution, H (F | PC1 = pc 1 ), should satisfy the equation below: In Equation (9), F is the vibration feature most highly correlated with PC1, which was described in the previous subsections. The parameter fm is the mode of H (F | PC1 = pc 1 ). Since H was assumed to correspond to high fluctuation initiation, it represents the FT probability distribution and fm is suggested as an FT.
The procedure proposed in this paper to estimate the FT of REBs is summarized in the flowchart depicted in Figure 2. First, different vibration features in the time and frequency domains were extracted from the run-to-failure history. These features were the inputs to the PCA algorithm that produced a feature as a degradation indicator (PC1). Afterward, the vibration feature most highly correlated with PC1 was selected. By using the copula models, a suitable joint probability distribution was fitted for the selected feature and PC1. Finally, the FT probability distribution was obtained based on the proposed criterion.

Introducing Two Bearing Datasets
In this section, two experimental run-to-failure datasets are introduced to study the effectiveness of the method proposed in this paper. Dataset 1 corresponds to seven accelerated life tests of an REB in the laboratory. Dataset 2 corresponds to three run-to-failure tests of an REB in an industrial machine. The details of these two datasets will be introduced in the following subsections.

Dataset 1: Accelerated Life Test Dataset
A set of run-to-failure experiments were designed and conducted for the purpose of this research. The platform of these experiments is shown in Figure 3. Using this platform, an accelerometer was installed vertically on the housing of the test REB. The test REB was

Introducing Two Bearing Datasets
In this section, two experimental run-to-failure datasets are introduced to study the effectiveness of the method proposed in this paper. Dataset 1 corresponds to seven accelerated life tests of an REB in the laboratory. Dataset 2 corresponds to three run-tofailure tests of an REB in an industrial machine. The details of these two datasets will be introduced in the following subsections.

Dataset 1: Accelerated Life Test Dataset
A set of run-to-failure experiments were designed and conducted for the purpose of this research. The platform of these experiments is shown in Figure 3. Using this platform, an accelerometer was installed vertically on the housing of the test REB. The test REB was a 6907 deep groove single-row bearing (55 mm OD, 35 mm ID). The accelerated life experiments were conducted under constant operating conditions of 2000 rpm speed and 9000 N radial load. The sampling frequency for the accelerometer was 25.6 kHz. In these experiments, reaching an amplitude of 30 g in the acceleration signals was the test stop criterion.  In total, seven run-to-failure experiments were conducted. Table 3 shows the final useful life results and the observed failure types in these experiments. At the end of the experiments, the observed failure in four tests related to the rolling elements and for the other three tests related to the inner race ( Figure 4). Rolling element 6498 7 Rolling element 5546  In total, seven run-to-failure experiments were conducted. Table 3 shows the final useful life results and the observed failure types in these experiments. At the end of the experiments, the observed failure in four tests related to the rolling elements and for the other three tests related to the inner race ( Figure 4).  In total, seven run-to-failure experiments were conducted. Table 3 shows the final useful life results and the observed failure types in these experiments. At the end of the experiments, the observed failure in four tests related to the rolling elements and for the other three tests related to the inner race ( Figure 4). Rolling element 6498 7 Rolling element 5546

Dataset 2: Industrial REB Dataset
As well as the accelerated life test dataset, the industrial data for a cylindrical roller bearing NU 2309 ECP (45 mm ID, 100 mm OD) of a sanding machine in a wood-based products manufacturing company was studied. These data were provided by Behravesh Vibration Engineering Company. The offline recorded vibration of this REB from March 2011 to June 2016 was provided for this research. In total, thirty-eight acceleration measurements in the frequency range of 10 to 8000 Hz were collected for this period. The trends for the RMS and the peak of the recorded signals in the mentioned period of time are depicted in Figure 5. During this period, three REB failures were detected and the machine was repaired. These three points are illustrated with a star (*) in the graph. These measurements were considered as three run-to-failure tests and their use in the study is described in the next section.

Dataset 2: Industrial REB Dataset
As well as the accelerated life test dataset, the industrial data for a cylindrical roller bearing NU 2309 ECP (45 mm ID, 100 mm OD) of a sanding machine in a wood-based products manufacturing company was studied. These data were provided by Behravesh Vibration Engineering Company. The offline recorded vibration of this REB from March 2011 to June 2016 was provided for this research. In total, thirty-eight acceleration measurements in the frequency range of 10 to 8000 Hz were collected for this period. The trends for the RMS and the peak of the recorded signals in the mentioned period of time are depicted in Figure 5. During this period, three REB failures were detected and the machine was repaired. These three points are illustrated with a star (*) in the graph. These measurements were considered as three run-to-failure tests and their use in the study is described in the next section.

Results
This section describes the application of the proposed method for FT determination to the two datasets introduced previously.

Calculating FT for Dataset 1 (Accelerated Life Test Dataset)
Three different scenarios were designed and considered in order to analyze the experimental data for this section. The results of the scenarios are presented as three models. In the first model, the features of all the data, without considering the failure modes, were given to the PCA algorithm to obtain the degradation indicator PC1. Then, a joint probability distribution was fitted between PC1 and the best correlated feature. With this distribution in hand, the approach proposed in Section 2.4 was applied to determine the FT. In the next two models, the experimental data were separated in order to be analyzed with regard to their occurred failure mode. In the second model, PC1 was extracted similarly to the previous model. However, data classified according to the failure mode was used to extract the joint probability distribution separately. Thus, considering separated joint distributions led to specific FTs for each failure mode. In the third model, a specific PC1 was extracted for each failure mode. Then, the FT was estimated similarly to the second model.
The reason for introducing the second and third models was to include consideration of the failure mode in the FT estimation. Model 1 did not include consideration of the failure mode at all. Model 2 included consideration of the failure mode only to generate the copula models. However, Model 3 considered the failure mode in the PCA algorithm as well as in the generation of the copula models.

Results
This section describes the application of the proposed method for FT determination to the two datasets introduced previously.

Calculating FT for Dataset 1 (Accelerated Life Test Dataset)
Three different scenarios were designed and considered in order to analyze the experimental data for this section. The results of the scenarios are presented as three models. In the first model, the features of all the data, without considering the failure modes, were given to the PCA algorithm to obtain the degradation indicator PC 1 . Then, a joint probability distribution was fitted between PC 1 and the best correlated feature. With this distribution in hand, the approach proposed in Section 2.4 was applied to determine the FT. In the next two models, the experimental data were separated in order to be analyzed with regard to their occurred failure mode. In the second model, PC 1 was extracted similarly to the previous model. However, data classified according to the failure mode was used to extract the joint probability distribution separately. Thus, considering separated joint distributions led to specific FTs for each failure mode. In the third model, a specific PC 1 was extracted for each failure mode. Then, the FT was estimated similarly to the second model.
The reason for introducing the second and third models was to include consideration of the failure mode in the FT estimation. Model 1 did not include consideration of the failure mode at all. Model 2 included consideration of the failure mode only to generate the copula models. However, Model 3 considered the failure mode in the PCA algorithm as well as in the generation of the copula models.

Model 1: Proposing the FT without Considering the Failure Modes in PCA and Copulas
The time and frequency domain features (Table 1) were extracted from the run-tofailure test data. These features were given to the PCA algorithm as inputs and PC 1 was calculated as the output of this algorithm. The five features most highly correlated with PC 1 are specified in Table 4. It can be seen that the peak of the time signal had the highest correlation with PC 1 . Next, using the statistical copula models (described in Table 2), the Gumbel copula was selected as the best to create the joint probability distribution function. The criteria of the model fitness evaluation (RMSE and NSE) and corresponding parameter values are listed in Table 5. Thus, the Gumbel copula was employed to provide the joint probability distribution of the peak and PC 1 (Figure 6). The time and frequency domain features (Table 1) were extracted from the run-tofailure test data. These features were given to the PCA algorithm as inputs and PC1 was calculated as the output of this algorithm. The five features most highly correlated with PC1 are specified in Table 4. It can be seen that the peak of the time signal had the highest correlation with PC1. Next, using the statistical copula models (described in Table 2), the Gumbel copula was selected as the best to create the joint probability distribution function. The criteria of the model fitness evaluation (RMSE and NSE) and corresponding parameter values are listed in Table 5. Thus, the Gumbel copula was employed to provide the joint probability distribution of the peak and PC1 ( Figure 6). Equation (6) can be used to obtain the conditional probability distribution function of the peak given PC1. According to the definition proposed for the FT in this paper, high vibration fluctuation occurred when PC1 was 13.5. The estimated distribution of the FT using the proposed method is represented in Figure 7. Moreover, the mode of distribution was 3.5 g. This value can be used as the constant value for the FT in predictions for the RUL problem. Equation (6) can be used to obtain the conditional probability distribution function of the peak given PC 1 . According to the definition proposed for the FT in this paper, high vibration fluctuation occurred when PC 1 was 13.5. The estimated distribution of the FT using the proposed method is represented in Figure 7. Moreover, the mode of distribution was 3.5 g. This value can be used as the constant value for the FT in predictions for the RUL problem.

Model 2: Proposing the FT by Considering the Failure Modes in Copulas
In this model, PC1 was obtained using a similar approach to that described for the previous model. However, in the current model, the data classified by the failure modes were used to extract the joint probability distribution separately. Thus, having separated joint distributions resulted in a specific FT for each failure mode.
The experimental data used in this study were classified into two groups according to the observed failure mode (inner race failure or rolling element failure) at the end of each experiment (Table 3). Following the analysis of the correlation between vibration features with PC1, the peak of the time signal had the highest correlation for both failure modes (Table 6). Therefore, using the statistical copula models (Table 2), Gumbel and Gaussian copulas were chosen as the best to create the joint distribution functions for inner race and rolling element failure modes, respectively. The criteria of the model fitness evaluation and corresponding parameter values are listed inTable 7.

Model 2: Proposing the FT by Considering the Failure Modes in Copulas
In this model, PC 1 was obtained using a similar approach to that described for the previous model. However, in the current model, the data classified by the failure modes were used to extract the joint probability distribution separately. Thus, having separated joint distributions resulted in a specific FT for each failure mode.
The experimental data used in this study were classified into two groups according to the observed failure mode (inner race failure or rolling element failure) at the end of each experiment (Table 3). Following the analysis of the correlation between vibration features with PC 1 , the peak of the time signal had the highest correlation for both failure modes (Table 6). Therefore, using the statistical copula models (Table 2), Gumbel and Gaussian copulas were chosen as the best to create the joint distribution functions for inner race and rolling element failure modes, respectively. The criteria of the model fitness evaluation and corresponding parameter values are listed inTable 7.
Finally, the conditional probability distribution function of the peak given PC 1 was obtained using Equation (6). According to the definition of the FT proposed in this paper, the high fluctuation of vibration occurred for the rolling element degradation and inner race degradation when PC 1 was equal to 15.1 and 35, respectively. Furthermore, the conditional distributions of the peak for the mentioned values of PC 1 were extracted from the corresponding joint probability distributions. These distributions, which were used as the FT distributions, are depicted in Figure 8. Additionally, the mode of each distribution can be used as the FT value for the degradation of each component (3.1 g for the rolling element and 8.7 g for the inner race).  Finally, the conditional probability distribution function of the peak given PC1 was obtained using Equation (6). According to the definition of the FT proposed in this paper, the high fluctuation of vibration occurred for the rolling element degradation and inner race degradation when PC1 was equal to 15.1 and 35, respectively. Furthermore, the conditional distributions of the peak for the mentioned values of PC1 were extracted from the corresponding joint probability distributions. These distributions, which were used as the FT distributions, are depicted inFigure 8. Additionally, the mode of each distribution can be used as the FT value for the degradation of each component (3.1 g for the rolling element and 8.7 g for the inner race). In this model, the classification of accelerated life tests was also considered in the PCA step. The data based on the rolling element or inner race failure modes were categorized and given to the PCA algorithm and PC1 was obtained for each failure mode separately. Then, considering the results of the correlation analysis displayed inTable 8, the peak feature for both failure modes was specified as the highest correlated feature. Comparisons of the goodness-of-fit parameters for different copulas are shown inTable 9. In light of the results shown in this table, the Gumbel copula was chosen as the best fitted copula model for both failure modes.  In this model, the classification of accelerated life tests was also considered in the PCA step. The data based on the rolling element or inner race failure modes were categorized and given to the PCA algorithm and PC 1 was obtained for each failure mode separately. Then, considering the results of the correlation analysis displayed inTable 8, the peak feature for both failure modes was specified as the highest correlated feature. Comparisons of the goodness-of-fit parameters for different copulas are shown inTable 9. In light of the results shown in this table, the Gumbel copula was chosen as the best fitted copula model for both failure modes.  Finally, by calculating the conditional probability distribution function of the peak given PC 1 , the FTs were specified for rolling element degradation and inner race degradation when PC 1 was equal to 11.2 and 40, respectively. Similar to the results reported for model 2, the FT distributions for the degradation of each failure mode are depicted in Figure 9. The mode values of the FT distributions were 5 g and 14.5 g for the rolling element degradation and the inner race degradation, respectively. These values can be considered as constant FTs.
Appl. Sci. 2021, 11 Finally, by calculating the conditional probability distribution function of the peak given PC1, the FTs were specified for rolling element degradation and inner race degradation when PC1 was equal to 11.2 and 40, respectively. Similar to the results reported for model 2, the FT distributions for the degradation of each failure mode are depicted in Figure 9. The mode values of the FT distributions were 5 g and 14.5 g for the rolling element degradation and the inner race degradation, respectively. These values can be considered as constant FTs.

Calculating FT for Dataset 2 (Industrial REB Dataset)
In order to explore the significance of FTs in industrial applications, the approach presented in this study was applied to dataset 2. Since this dataset included limited sets of run-to-failure data (three sets in total), only model 1 was investigated., Furthermore, models 2 and 3 could not be employed since the failure modes were not clear in this dataset. As can be seen in Table 10, the peak of the time signal had the highest correlation with the PC1 resulting from the PCA algorithm. By comparing various copula families, the Gumbel copula was obtained as the most accurate one to model the joint probability distribution of data (Table 11).  Table 11. Comparing the copula models' goodness-of-fit parameters.

Calculating FT for Dataset 2 (Industrial REB Dataset)
In order to explore the significance of FTs in industrial applications, the approach presented in this study was applied to dataset 2. Since this dataset included limited sets of run-to-failure data (three sets in total), only model 1 was investigated. Furthermore, models 2 and 3 could not be employed since the failure modes were not clear in this dataset. As can be seen in Table 10, the peak of the time signal had the highest correlation with the PC 1 resulting from the PCA algorithm. By comparing various copula families, the Gumbel copula was obtained as the most accurate one to model the joint probability distribution of data (Table 11).
The probability distribution of the FT was estimated through the aforementioned criterion ( Figure 10). This distribution showed that the mode value of 10.5 g could be used as a proper FT value in this industrial REB. The probability distribution of the FT was estimated through the aforementioned criterion ( Figure 10). This distribution showed that the mode value of 10.5 g could be used as a proper FT value in this industrial REB.

Discussion on the Results of Dataset 1
As previously mentioned, the experimental accelerated life test data included seven run-to-failure tests, classified into two groups according to the type of final observed failure (Table 3). This subsection describes the application of the results of the FT estimation by the three models in Section 4.1 to these run-to-failure tests in order to compare their outputs.
The graphs in Figure 11 illustrate the trends for the vibration peak feature for three tests corresponding to inner race failure. In these graphs, three proposed values for the FT using the three models described in Section 4.1 are shown. The EoL estimation of each model is depicted by a star (*) on the trend. As can be seen, the estimated EoLs of the first and second models were close together and represent more conservative failure points than the third model. However, all of the models prevented unstable vibration. Also, the third model provided a more extended lifetime.
Similarly, the trends for the vibration peak in the four tests with the rolling element failure are shown in Figure 12. In this figure, it can be seen that the first and second models estimated the same EoL for all four tests. The third model also estimated a similar EoL in three tests out of four. However, in one test (test 7), the EoL estimated by model 3 was somewhat later than the two other models.
The results of the three FT estimation models with regard to useful life for the seven Figure 10. Estimated FT distribution for dataset 2.

Discussion on the Results of Dataset 1
As previously mentioned, the experimental accelerated life test data included seven run-to-failure tests, classified into two groups according to the type of final observed failure (Table 3). This subsection describes the application of the results of the FT estimation by the three models in Section 4.1 to these run-to-failure tests in order to compare their outputs.
The graphs in Figure 11 illustrate the trends for the vibration peak feature for three tests corresponding to inner race failure. In these graphs, three proposed values for the FT using the three models described in Section 4.1 are shown. The EoL estimation of each model is depicted by a star (*) on the trend. As can be seen, the estimated EoLs of the first and second models were close together and represent more conservative failure points than the third model. However, all of the models prevented unstable vibration. Also, the third model provided a more extended lifetime.  Figure 11. Comparison of the FT models in the tests for inner race faults. Figure 11. Comparison of the FT models in the tests for inner race faults.

Test 3 Test 5
Similarly, the trends for the vibration peak in the four tests with the rolling element failure are shown in Figure 12. In this figure, it can be seen that the first and second models estimated the same EoL for all four tests. The third model also estimated a similar EoL in three tests out of four. However, in one test (test 7), the EoL estimated by model 3 was somewhat later than the two other models.
The results of the three FT estimation models with regard to useful life for the seven tests are given in Table 12. The parameter ∆ is used to show the difference in useful life between the proposed models and the predetermined test stop criterion. The reader should note that this difference does not represent the error of the models presented in this paper since the predetermined criterion is just an assumption for stopping the tests and not the ideal. Figure 11. Comparison of the FT models in the tests for inner race faults.

Test 3 Test 5
Appl. Sci. 2021, 11, x FOR PEER REVIEW 15 of 18 Test 6 Test 7 Figure 12. Comparison of the FT models in the tests for rolling element faults. The highest value of ∆ was obtained in test 2. Comparing the vibration trends from all seven tests reveals the different behavior of the REB degradation in test 2. The vibration trend started to grow gradually from the beginning of its life. Furthermore, this REB had a relatively short life. The lack of a healthy stage in this test (unlike in the others) led to the high value of ∆.
Comparing the degradation trends for two groups of data shows that the degradation data included a soft failure stage when the growing defect was on the inner race. However, when the defect was on the rolling element, the soft failure stage did not exist. So, sudden growth and hard failure occurred after a healthy stage. The physical reason behind this behavior difference can be explained by the stress states in the rolling element and inner race. The smaller curvature radius in the rolling element (compared to the inner race) led to more severe stress and resulted in quicker degradation. Another explanation  The highest value of ∆ was obtained in test 2. Comparing the vibration trends from all seven tests reveals the different behavior of the REB degradation in test 2. The vibration trend started to grow gradually from the beginning of its life. Furthermore, this REB had a relatively short life. The lack of a healthy stage in this test (unlike in the others) led to the high value of ∆.
Comparing the degradation trends for two groups of data shows that the degradation data included a soft failure stage when the growing defect was on the inner race. However, when the defect was on the rolling element, the soft failure stage did not exist. So, sudden growth and hard failure occurred after a healthy stage. The physical reason behind this behavior difference can be explained by the stress states in the rolling element and inner race. The smaller curvature radius in the rolling element (compared to the inner race) led to more severe stress and resulted in quicker degradation. Another explanation for this behavior could be the sharp edge flattening of spall regions that occurs in the inner race after removing particles from the surface. However, this flattening does not occur in the spall region of rolling elements [28].
In the third model proposed in this paper, the type of growing defect was considered in detail. Its output was therefore more accurate in estimating a proper EoL in the tests. However, model 1 proposed an EoL regardless of the type of failure and its output was more conservative than the others.

Discussion on the Results of Dataset 2
The peak trend for the industrial data on the three REB lifespans is depicted in a single plot in Figure 13. The total lifespans and the fluctuations of the vibrations in the three REBs can be compared in this figure. The FT extracted by analyzing these data with the approach proposed in this paper is depicted by a red dashed line in the figure.  Keeping in mind that there is no suggested FT for the acceleration of REBs in the guidelines, comparing the output of the proposed method with the actual decision of the condition monitoring (CM) group, which is based on their own experiences with the machine, shows that the proposed FT seems reasonable.
It should be noted that, since these data were gathered from an actual industrial field, the failure point in this dataset does not correspond to catastrophic failure and repairing activities were conducted to prevent such failure, based on CM recommendations. However, in the accelerated life tests which were previously discussed, the final failure referred to more severe and controlled damage under laboratory conditions.
In industrial applications, many of the machines are under the offline CM program. The available record during the degradation period is thus very limited. Since the goodness of fit in copula model estimation directly depends on the amount of available data, this limited data could affect the performance of the proposed algorithm.

Conclusions
In this paper, a method was suggested to determine the FT of REBs using vibration data. In this method, PCA was used to produce PC1 as a degradation growth indicator from vibration features. Moreover, copula models were applied to reveal the relationship between the vibration features and PC1. With the best-fitted copula model in hand, the FT Keeping in mind that there is no suggested FT for the acceleration of REBs in the guidelines, comparing the output of the proposed method with the actual decision of the condition monitoring (CM) group, which is based on their own experiences with the machine, shows that the proposed FT seems reasonable.
It should be noted that, since these data were gathered from an actual industrial field, the failure point in this dataset does not correspond to catastrophic failure and repairing activities were conducted to prevent such failure, based on CM recommendations. However, in the accelerated life tests which were previously discussed, the final failure referred to more severe and controlled damage under laboratory conditions.
In industrial applications, many of the machines are under the offline CM program. The available record during the degradation period is thus very limited. Since the goodness of fit in copula model estimation directly depends on the amount of available data, this limited data could affect the performance of the proposed algorithm.

Conclusions
In this paper, a method was suggested to determine the FT of REBs using vibration data. In this method, PCA was used to produce PC 1 as a degradation growth indicator from vibration features. Moreover, copula models were applied to reveal the relationship between the vibration features and PC 1 . With the best-fitted copula model in hand, the FT distribution was defined as a probability distribution based on a 25% vibration fluctuation criterion and the mode of this distribution was used as a constant value for the FT. The proposed method was applied to one accelerated life test dataset as well as an industrial dataset. Applying the method to the accelerated life tests led to the development of three different models for the failure mode. The third model, which used the separated data of each failure mode for PCA and copula analysis, provided more accurate FTs that took the type of failure into account. It was shown that the estimated FT for the rolling element was less than that for the inner race. this was because of prompt growth of degradation after the initiation of the defect in the rolling element. Additionally, using model 1 to estimate the FT regardless of the failure mode led to a more conservative FT. Applying the method to the industrial data showed that the proposed method provided reasonable FTs for actual field data too. Due to the limitations of industrial data, only model 1 was used to estimate the FTs for this dataset.
FT selection plays a significant role in the RUL prediction of REBs. Since the method proposed in this paper estimates FTs using the high fluctuations in the degradation trend, it can provide a more accurate RUL prediction. Also, the provided FT distribution can be used in the reliability analysis of rotary machines. Furthermore, the method proposed in this paper can be used to estimate FTs under different machine operating conditions, which is what we plan to investigate in our future work.