Dynamic Evaluation of the Degradation Process of Vibration Performance for Machine Tool Spindle Bearings

Real-time condition monitoring and fault diagnosis of spindle bearings are critical to the normal operation of the matching machine tool. In this work, considering the interference of random factors, the uncertainty of the vibration performance maintaining reliability (VPMR) is introduced for machine tool spindle bearings (MTSB). The maximum entropy method and Poisson counting principle are combined to solve the variation probability, so as to accurately characterize the degradation process of the optimal vibration performance state (OVPS) for MTSB. The dynamic mean uncertainty calculated using the least-squares method by polynomial fitting, fused into the grey bootstrap maximum entropy method, is utilized to evaluate the random fluctuation state of OVPS. Then, the VPMR is calculated, which is used to dynamically evaluate the failure degree of accuracy for MTSB. The results show that the maximum relative errors between the estimated true value and the actual value of the VPMR are 6.55% and 9.91%, and appropriate remedial measures should be taken before 6773 min and 5134 min for the MTSB in Case 1 and Case 2, respectively, so as to avoid serious safety accidents that are caused by the failure of OVPS.


Introduction
Precision, as one of the most important performance indicators of spindle bearings, has a significant impact on the normal operation of machine tools. Super-precision spindle bearings (SPSB) refer to spindle bearings with low vibration, a wide speed range, a high rotational accuracy, low heating, high system rigidity, and low noise. SPSB works at the OVPS, which is the basis for a machine tool to achieve optimal performance. The failure of MTSB are not only due to fatigue. A variety of performance failures, such as stuck, sintering, plastic deformation, crack, or fracture may occur before fatigue failure [1][2][3]. The probability distribution of these failure modes is unknown, and the characteristic data of the spindle bearings are scarce. In particular, the non-linear dynamic contact and collision between the parts inside the spindle bearings, the non-linear damage, viscoustemperature and viscous-pressure characteristics of lubricating medium, and the accuracy loss, all present uncertain and non-linear characteristics [4][5][6][7]. These uncertainties bring new challenges to analyze the degradation trend of vibration performance. Therefore, it is of great significance to evaluate the degradation process of vibration performance dynamically before the MTSB fail.
When analyzing the degradation process of bearings, the mathematical statistical methods generally assume that the amount of data are limited, and then extract characteristic parameters of bearings [8][9][10]. Based on the vibration signal collected, Ye et al. [11] used the maximum entropy method to calculate the PDF of bearings, which was regarded as degradation characteristics. Then, the Bayesian method, the bootstrap method and the used the maximum entropy method to calculate the PDF of bearings, which was regarded as degradation characteristics. Then, the Bayesian method, the bootstrap method and the maximum entropy method were fused to evaluate the reliability of bearings in service. Its performance degradation and reliability evaluation process are shown in Figure 1. However, due to the differences of experimental data samples caused by random factors and interference factors in the experimental process, or the incompleteness of information of data samples, it is difficult to represent the entire distribution range of bearing life accurately. For two sets of bearings with the same batch and model, the prediction results may also be different, that is, the generalization ability and stability of the model cannot be guaranteed. In order to reduce the impact of random and interfering factors on the analysis results, many experts and scholars at home and abroad preprocess the experimental data. Huang et al. [12] fused empirical mode decomposition method and Hilbert transform method to analyze the time-frequency spectrum of signals, making up for the shortcomings of traditional statistical theory in analyzing nonlinear and non-stationary data. In order to avoid multiple spline curve fitting in the process of empirical mode decomposition, Feldman [13] proposed a time-varying vibration decomposition method based on Hilbert transform. In order to not affect the accuracy of condition monitoring and fault diagnosis, Guo et al. [14] proposed a signal compression method based on integrated empirical mode decomposition, which adaptively decomposes the vibration signal into signal components with different frequency bands. These preprocessing methods may reduce the influence of random and interference factors, while also filtering out some useful information, resulting in "distortion" of the analysis results, which in turn leads to the inability to evaluate the degradation process of bearings effectively. However, due to the differences of experimental data samples caused by random factors and interference factors in the experimental process, or the incompleteness of information of data samples, it is difficult to represent the entire distribution range of bearing life accurately. For two sets of bearings with the same batch and model, the prediction results may also be different, that is, the generalization ability and stability of the model cannot be guaranteed. In order to reduce the impact of random and interfering factors on the analysis results, many experts and scholars at home and abroad preprocess the experimental data. Huang et al. [12] fused empirical mode decomposition method and Hilbert transform method to analyze the time-frequency spectrum of signals, making up for the shortcomings of traditional statistical theory in analyzing nonlinear and non-stationary data. In order to avoid multiple spline curve fitting in the process of empirical mode decomposition, Feldman [13] proposed a time-varying vibration decomposition method based on Hilbert transform. In order to not affect the accuracy of condition monitoring and fault diagnosis, Guo et al. [14] proposed a signal compression method based on integrated empirical mode decomposition, which adaptively decomposes the vibration signal into signal components with different frequency bands. These preprocessing methods may reduce the influence of random and interference factors, while also filtering out some useful information, resulting in "distortion" of the analysis results, which in turn leads to the inability to evaluate the degradation process of bearings effectively.
Traditional machine learning and deep learning, as methods widely used in data preprocessing, have their own advantages and disadvantages. In the field of traditional machine learning, SVM has a low error rate of generalization, a low computational overhead, and it is easy to interpret its results [15,16]. However, SVM is too sensitive to parameters of Sensors 2023, 23, 5325 3 of 27 kernel functions. The idea of KNN is simple, with no assumptions about data input, but its disadvantage is that the computational complexity is high [17]. The logistic regression method has low computational cost, and it is easy to understand and implement [18]. However, it is prone to underfitting and its classification accuracy is not high. The decision tree is insensitive to missing intermediate values, and can handle irrelevant feature data [19]. It is better than KNN in understanding the inherent meaning of the data. Its disadvantage is that it is prone to overfitting and the construction process is time-consuming.
Deep learning and traditional machine learning are similar in data preprocessing. The core difference lies in the feature extraction process, where deep learning does not require manual extraction and the extraction process is performed by a machine [20][21][22][23]. Although deep learning can learn the features of patterns automatically and achieve high recognition accuracy, the prerequisite is that amounts of data are provided. For a limited amount of data, deep learning algorithms cannot estimate the laws of data without bias. The weight parameters of CNN are less than those of DNN connected fully, which makes the training speed of CNN model faster and is not prone to overfitting [24]. Meanwhile, CNN requires less data to train. The disadvantage is that it requires large sample size and parameter adjustment, and its physical meaning is also unclear. RNN cannot solve the problem of long-term dependence, while LSTM implements temporal memory function and prevents gradient disappearance. At the same time, LSTM can better handle the tasks of time series than CNN. In addition, it also solves the long-term dependency problem [25][26][27]. However, the model structure of LSTM is complex relatively, and its training process is more timeconsuming than CNN. The characteristics of RNNs determine that they cannot parallelize data well. It is also difficult for LSTM to handle longer data sequences.
Vibration data collected through the experiment vary with time, which can be regarded as a non-linear time series. From the beginning of service to the failure of accuracy, the vibration performance varies continuously, forming multiple time series. The intrinsic series refers to the time series with OPS. According to the time series of MTSB, the degradation process of accuracy is defined as a Poisson process, and the degradation probability is taken as a parameter of the counting process [28]. Influenced by many factors, the degradation process of bearings has the characteristics of non-linear dynamics [29,30]. For each time series of accuracy, the degradation probability relative to the intrinsic series also has the characteristics of non-linearity and diversity. This causes a dynamic change in information. In particular, the degradation probability functions of accuracy vary with time and environmental factors. In view of this, a polynomial is used to fit the parameters, and the data samples of the degradation probability of OVPS are obtained for each time series. Since it is difficult to obtain enough original information of degradation probability in a short time, the grey bootstrap method is used to generate a large set of sample data [31].
Information fusion technology and the real-time update method can improve the prediction accuracy, but most of the existing fusion methods focus on the fusion of homogeneous information collected by the sensors with same type. They rarely consider the fusion of heterogeneous information, especially the fusion and prediction of event single value data and waveform data monitored, which needs further in-depth study [32,33]. In addition, the current real-time update methods basically take Bayesian method as the theoretical framework. They can only apply the observed data to update the priori probability distribution, but cannot utilize other available information such as moments of parameters or functions of moments.
On the basis of this, a dynamic evaluation model is established to study the performance degradation process for MTSB. This model considers the interference of random factors during the service process of MTSB. First, the vibration signal are processed into segments to obtain the time series with OVPS based on rolling average method. The maximum entropy method is used to calculate the PDF of intrinsic sequence. Based on the Poisson counting principle, the variation probabilities of OVPS of each time series are calculated. Considering the interference of random factors, the small data sample of variation probability is obtained for each time series by polynomial fitting method. The segments to obtain the time series with OVPS based on rolling average method. Th imum entropy method is used to calculate the PDF of intrinsic sequence. Based Poisson counting principle, the variation probabilities of OVPS of each time series culated. Considering the interference of random factors, the small data sample of va probability is obtained for each time series by polynomial fitting method. The gre strap method (namely, GBM (1,1)) is applied to generate enough data samples of va probability, so as to obtain the estimated true value, estimated interval and dynam certainty of variation probability. The PMR and PMRR of OVPS of MTSB are calc according to Poisson process theory. Experimental verification is utilized to verify t proposed model can be used effectively to estimate the failure degree of OVPS to m and analyze the performance degradation process of MTSB.
The flow diagram of proposed method is shown in Figure 2.

Mathematical Models
Based on the maximum entropy method and the Poisson counting principle, the variation probability of the OVPS is calculated for MTSB. The small data sample of variation probability is obtained by polynomial fitting for each time series. The estimated true value curve and the upper-and lower-bound curves of variation probability are obtained by fusing the grey bootstrap method into the maximum entropy method. The VPMR is calculated based on the Poisson process, which is used to dynamically evaluate the failure degree of OVPS for MTSB.

Calculating Variation Probability of OVPS
During the service period of the spindle bearings, r time sequences are obtained by periodic sampling vibration acceleration signals of spindle bearings. The time series collected after initially wearing is considered to be the intrinsic sequence, which is marked as the first time series and expressed by the vector X 1 as where x(k) is kth performance data in the intrinsic sequence; k is the order number of performance data in intrinsic sequence, k = 1, 2, 3, . . . , N; N is the total number of performance data in the intrinsic sequence.
With the change of time variable, vibration acceleration data are collected continuously, and the nth time series vector X n is obtained.
where x n (k) is the kth performance data of the nth time series; n is the order number of time series, n = 1, 2, 3, . . . , r.
The maximum entropy method can make the optimal estimation for the unknown probability distribution with minimal subjective bias. Lagrange multipliers are introduced in the process of solving probability distributions, so the problem of solving probability distribution is transformed into the problem of solving Lagrange multipliers. For the convenience of description, the continuous variable x is used to express the discrete variable x(k) in the intrinsic sequence.
According to the maximum entropy method, the probability density function should meet the condition that the value of entropy function H(x) is the maximum.
where f (x) is the probability density function of continuous variable x; lnf (x) is the logarithm of the probability density function f (x); S is the feasible domain of the performance random variable x, S = [S 1 , S 2 ]; S 1 is the lower-bound value of the feasible domain, and S 2 is the upper-bound value of the feasible domain. Then, the Lagrange multiplier method is used to solve this problem by adjusting the probability density function f (x) to maximize the value of entropy function H(x) [11]. Assume H is the Lagrange function.
where i is the order number of origin moment, i = 1, 2, . . . , j, usually j = 5; m i is the ith order origin moment, m 0 = 1; x i is the coefficient of the function f (x); c i is the (i + 1)th Lagrange multiplier and c 0 is the first Lagrange multiplier, i = 1, 2, . . . , j. The probability density function f (x) of data samples can be expressed as With The other j Lagrange multipliers should satisfy that To ensure the convergence of solution, the original data interval is mapped to interval [−e, e] by the substitution of variable. The sample data are divided into ξ groups in the incremental order, and the histogram can be drawn. At the same time, the values z µ and frequency Γ µ can be calculated, and u = 2, 3, . . . , ξ + 1. Then, the histogram is extended into (ξ + 2) group and let Γ 1 = Γ ξ +2 . Let where ψ is the variable to be transformed, ψ ∈ [−e, e]; a and b are the mapping parameters; e has a value of 2.71828.
The mapping parameters a and b can be calculated based on dx = dψ/a.
Therefore, the probability density function f (x) can be transformed into Set a significant level α, α ∈ [0, 1], so the confidence level P is given by The maximum entropy estimated interval is [X L , X U ] for a given confidence level P, and the lower-bound value X L should meet that X L = X α/2 . And The upper-bound value should meet that X U = X (1−α)/2 . And According to Equation (16), the maximum entropy estimated interval [X L1 , X U1 ] can be calculated for the intrinsic time series, where X L1 is the lower-bound value and X U1 is the upper-bound value of the maximum entropy estimated interval of the intrinsic time series.
Based on the Poisson counting principle, record the number N n that performance data of the nth time series is outside the estimated interval [X L1 , X U1 ]. Then, the variation frequency λ n can be given as Equation (17) for the nth time series. With where N n1 is the number showing that performance data are less than X L1 for the nth time series; N n2 is the number showing that performance data are more than X U1 for the nth time series; n is the sequence number of time series, n = 1, 2, 3, . . . , r.

Calculating Estimated Truth Value and Estimated Interval of Variation Probability
The least-squares method is used to fit the variation probability with different order polynomials. The fitting effect depends on the correlation coefficient R 2 . The closer the correlation coefficient R 2 is to 1, the better the polynomial fitting effect. If the correlation coefficient R 2 is less than 0.8, the fitting effect is worse, so the polynomial will not be used in the later analysis process.
The polynomial function G q (λ) is given by G q (λ) = p q0 λ 0 + p q1 λ 1 + p q2 λ 2 . . . + p qγ λ γ + . . . + p qq λ q ; γ = 0, 1, 2, . . . , q; q = 1, 2, . . . , 6; where G q (λ) is the qth order polynomial; q is the order number of polynomial function; p qγ is the coefficient of the power function λ γ . According to the above polynomial function models, small data samples of variation probability of OVPS are obtained for each time series. Y(n) = (y n (1), y n (2), . . . , y n (6)) = (y n (u)); u = 1, 2, . . . , 6;n = 1, 2, . . . , r; (20) where Y(n) is the data sample of variation probability of OVPS for the n time series; y n (u) is the uth data in the variation probability data sample for the n time series. Using the bootstrap method, B bootstrap re-sampling samples of size z, namely the bootstrap re-sampling samples V Bootstrap , can be obtained by an equiprobable sampling as where V β is the βth bootstrap re-sampling sample, β = 1, 2, . . . , B; B is the times of the bootstrap re-sampling, and also the number of bootstrap samples, with where v β (Θ) is the Θth data in the βth bootstrap re-sampling sample.
According to the grey model GM (1,1) [34], suppose the first-order accumulated generating operator (1-AGO) of V β is given by The grey generated model can be described as the differential equation as follows: where u is the time variable, and c 1 and c 2 are the coefficients to be estimated. Use the increment to replace the differential, viz., where ∆u is equal to the unit interval, 1. Furthermore, assume the generated vector of the mean series as follows The least-squares solution with the initial condition y ξβ (1) = v ξβ (1) is given bŷ where the coefficients c 1 and c 2 are as follows According to the inverse AGO [34], the βth generated data are expressed aŝ Therefore, B generated data for the vibration performance can be obtained as Y B = w 1 , w 2 , . . . , w β , . . . , w B = ν 1 (z + 1),ν 2 (z + 1), . . . ,ν β (z + 1), . . . ,ν B (z + 1) (32) where w β is the βth generated data. The maximum entropy method is used to calculate the probability density function of the generated sample Y B . According to the probability density function, the true value and upper-and lower-bound values are estimated for the data sample of variation probability of each time series.
The probability density function of variation probability data sample can be calculated as The estimated true value of the variability probability is obtained as Set a significant level and let α ∈ (0, 1). The maximum entropy estimated interval is given as where λ nL is the lower-bound value and λ nU is the upper-bound value of the variation probability data sample for the nth time series.

Evaluation of the Uncertainty of Variation Probability
The fluctuation range of variation probability can be expressed as where U λn is the estimated uncertainty of variation probability, namely, the instantaneous uncertainty at the confidence level P. The number η is calculated for variation probability sample data more than the upperbound value λ nU based on the Poisson counting process, and the reliability of evaluation result is defined as where P R is the reliability of the polynomials fitting effect using the least-squares method. Define where U mean is the dynamic average uncertainty; | PR=100% stands for that the calculation process is under the condition of P R = 100% [9].

Evaluation of PMR and PMRR
Any counting process can be described by the Poisson process [28] as where ξ stands for the time variable with ξ = 1, 2, 3, . . . , and ξ ≥ 1; e is the number of occurring failure events with e = 0, 1, 2, 3, . . . , namely, the serious variation in working condition that may cause the OVPS failure; Q is the probability of failure events occurring e times. Thus, the reliability R for events occurring can be obtained by the Poisson process. When solving for the PMR of bearings working at the OVPS, let e = 0, viz., which indicates that R is the probability before the OVPS fails. Let ξ = 1, which indicates that R is the PMR of time series in current time. According to Equation (40), PMR can be expressed as where R(λ n ) is a function of the variation probability λ n . The upper-and lower-bound values of variation probability are applied into the reliability Equation (41), so the estimated true values R 0 and upper-and lower-bound intervals [R L , R U ] are gained for performance maintaining reliability during the time intervals corresponding to the performance time series. The range of variation probability λ n is [0, 1]. 0 represents the OVPS of bearings without any variation, which is an ideal state and the reliability is 100% for bearing working at the OVPS. 1 represents that the OVPS of rolling bearings fail completely and the performance state is very unreliable. Therefore, if the value of λ nL is less than 0, let λ nL = 0 artificially in the process of solving the performance maintaining reliability.
According to the concept of relative error in measurement theory, performance maintaining relative reliability (PMRR) d(λ n ) of MTSB is obtained to characterize the failure degree of OVPS [11].
where R(λ 1 ) is the PMR for the intrinsic series of MTSB, R(λ n ) is the PMR for the nth time series of MTSB; d(λ n ) is the PMRR for the nth time series of MTSB. The basic classification principle of degree failure of OVPS for MTSB is as follows: (1) If d(λ n ) is not less than 0%, which shows that the PMR during this period is not less than PMR of OVPS, and it cannot deny that the performance has reached its optimal state; otherwise, it can deny that the performance has achieved its optimal state. (2) When d(λ n ) is less than 0%, if the absolute value of d(λ n ) is in (0%, 15%], this indicates that the error between the evaluation value and the optimum value is very small. If the absolute value of d(λ n ) is in (15%, 30%], this indicates that the error between the evaluation value and the optimum value is gradually increasing. If the absolute value of d(λ n ) is greater than 30%, this indicates that the error between the evaluation value and the optimum value is very large.
Based on that, the degree failure of OVPS for MTSB is divided into S1, S2, S3, S4 for a total of four levels: S1: If d(λ n ) ≥ 0%, it indicates the performance states of MTSB reaches the optimum and has almost no failure possibility.
S2: If d(λ n ) ∈ [−15%, 0%), it indicates the performance states of MTSB is normal, and the degree failure of OVPS is very small. S3: If d(λ n ) ∈ [−30%, −15%), it indicates the performance states of MTSB is gradually becoming worse, and the degree failure of OVPS is gradually increasing.
S4: If d(λ n ) < −30%, it indicates the performance states of MTSB is worse, and the degree failure of OVPS is very large.
The negative value of d(λ n ) indicates the performance states has attenuation, namely, PMR currently is less than PMR during the optimum period, and the positive value indicates no attenuation. The smaller d(λ n ) is, the worse the performance states of MTSB is, and the larger the degree failure of OVPS is. Therefore, the time corresponding to d(λ n ) = −30%, is the critical time where OVPS becomes poor. Taking corresponding measures before the critical time can avoid serious safety accidents that are caused by the failure of OVPS.

Case 1
This is a strength lifetime test on the machine tool spindle bearings, which is conducted at a radial load of 4.58 kN, an axial load of 4.17 kN, and a motor speed 6000 r/min. The vibration acceleration sensor has a measuring range of ±2000 g and a resolution of 0.0001 g. The data-acquisition rate is 20 KHz. The test machine, the bearing used and the vibration data in this case are exactly the same as those of Case 1 in Reference [11]. The vibration signals are automatically collected by the computer control system, as shown in Figure 3.
As shown in Figure 3, based on the rolling average method, the performance degradation stages are divided as shown in Table 1.
From the OVPS of the bearing, there are 7321 vibration signals. Among them, the 473rd to 7472nd data points are divided into 10 segments at 700 data intervals. The 7473rd to 7793rd data are separated into one segment, that is, the 11th segment.
It is worth noting that if the vibration data are divided into fewer segments, the specific variation process of characteristic parameter cannot be obtained. If the vibration data are divided into a large number of segments, it will cause large computational workload and complex calculation process. At the same time, each segment of data contain less information, which makes the calculation error of the probability density function larger. Therefore, the vibration data are divided into ten segments usually in the analysis process. Sensors 2023, 23, x FOR PEER REVIEW 11 of 32

Data Point Degradation Stage 1st to 472nd
Initial wear stage 473rd to 2574th Optimal performance state 2575th to 6659th Normal wear stage 6660th to 7446th Degeneration stage 7447th to 7793rd Deterioration stage From the OVPS of the bearing, there are 7321 vibration signals. Among them, the 473rd to 7472nd data points are divided into 10 segments at 700 data intervals. The 7473rd to 7793rd data are separated into one segment, that is, the 11th segment.
It is worth noting that if the vibration data are divided into fewer segments, the specific variation process of characteristic parameter cannot be obtained. If the vibration data are divided into a large number of segments, it will cause large computational workload and complex calculation process. At the same time, each segment of data contain less information, which makes the calculation error of the probability density function larger. Therefore, the vibration data are divided into ten segments usually in the analysis process.     Assume that the significance level α is 0.01 and the confidence level is P = 99%, based on the Equations (13)- (16), the maximum entropy estimated interval of the intrinsic series is [2.4867, 4.5127] m·s −2 .
Based on the Equations (17) and (18), the numbers are calculated, respectively, that the 11 data samples fall outside the maximum entropy estimated interval [X L1 , X U1 ] of intrinsic time series. According to the Poisson counting principle, the variation frequencies [λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 , λ 7 , λ 8 , λ 9 , λ 10 , λ 11 ] are as shown in Figure 5. As shown in Figure 5, the performance variation probability is non-linear and uncertain-relative to the intrinsic time series. Before the corresponding time interval of the third time series, the variation probability of the MTSB vibration performance is almost zero. Between the corresponding time interval of the third time series and the fifth time series, the variation probability has an increased trend. Between the corresponding time interval of the fifth time series and the eighth time series, the variation probability has a decreasing trend. Between the corresponding time interval of the eighth time series and the eleventh time series, the variation probability has a rapidly increasing trend.

Estimated Truth Value and Estimated Interval of Variation Probability of MTSB
The least-squares method is used to fit the variation probability with second-, third-, fourth-, fifth-and sixth-order polynomials. Based on the Equation (19), the fitting results are shown in Tables 2 and 3, and Figure 6.   As shown in Figure 5, the performance variation probability is non-linear and uncertain-relative to the intrinsic time series. Before the corresponding time interval of the third time series, the variation probability of the MTSB vibration performance is almost zero. Between the corresponding time interval of the third time series and the fifth time series, the variation probability has an increased trend. Between the corresponding time interval of the fifth time series and the eighth time series, the variation probability has a decreasing trend. Between the corresponding time interval of the eighth time series and the eleventh time series, the variation probability has a rapidly increasing trend.

Estimated Truth Value and Estimated Interval of Variation Probability of MTSB
The least-squares method is used to fit the variation probability with second-, third-, fourth-, fifth-and sixth-order polynomials. Based on the Equation (19), the fitting results are shown in Tables 2 and 3, and Figure 6. Second order G 2 (λ n ) = 0.0155λ n 2 − 0.1162λ n + 0.2065 0.7863 Third order G 3 (λ n ) = 0.0048λ n 3 − 0.0712λ n 2 + 0.3181λ n − 0.3193 0.9390 Fourth order G 4 (λ n ) = 0.0009λ n 4 − 0.0162λ n 3 + 0.0958λ n 2 − 0.1751λ n + 0.0900 0.9725 Fifth order G 5 (λ n ) = −0.0002λ n 5 + 0.0076λ n 4 − 0.0900λ n 3 + 0.4533λ n 2 − 0.9059λ n + 0.5576 0.9860 Sixth order G 6 (λ n ) = −0.00008λ n 6 + 0.0027λ n 5 − 0.0335λ n 4 + 0.1897λ n 3 − 0.4996λ n 2 + 0.5821λ n − 0.2323 0.9955  The closer the correlation coefficient R 2 is to 1, the better the polynomial fitting. As shown in Table 4, the correlation coefficient R 2 for the second-order polynomial is less than 0.8. This shows that the fitting is worse than other higher-order polynomials. Therefore, the second-order polynomial is not used in further analysis. The closer the correlation coefficient R 2 is to 1, the better the polynomial fitting. As shown in Table 4, the correlation coefficient R 2 for the second-order polynomial is less than 0.8. This shows that the fitting is worse than other higher-order polynomials. Therefore, the second-order polynomial is not used in further analysis.  Figure 6 shows that the fitted variation probability is almost identical with the actual variation probability for the fourth, tenth and eleventh time series. The fitting was performed using the least-square method by the third-order polynomial. With fittings using the least-square method by the fourth-order polynomial, the fitted variation probability values are basically identical with the actual values for the second, fourth, sixth, seventh and ninth time series. The fitted variation probability values are almost identical with the actual variation probability values for the third, sixth, ninth and eleventh time series when fitting was performed using the least-square method by the fifth-order polynomial. When fitting was performed using the least-square method by the sixth-order polynomial, the fitted variation probability values are basically identical with the actual values for the first, third, fourth, fifth, seventh, eighth and ninth time series. In short, the above four polynomials have their own advantages in fitting. Therefore, relevant information should be fully fused to effectively monitor the evolution process of rolling bearing vibration performance.

The Uncertainty of Variation Probability of MTSB
The variation probability of the first time series is taken as an example here. According to the grey bootstrap method, the number of bootstrap re-sampling is taken to be B = 20,000, and confidence level is taken to be P = 95%. Based on the Equations (20)- (32), the grey bootstrap sample Y Bootstrap is obtained by sampling the sample data of fitting the variation probability values using the above four polynomials, as shown in Figure 7.

The Uncertainty of Variation Probability of MTSB
The variation probability of the first time series is taken as an example here. According to the grey bootstrap method, the number of bootstrap re-sampling is taken to be B = 20,000, and confidence level is taken to be P = 95%. Based on the Equations (20)-(32), the grey bootstrap sample YBootstrap is obtained by sampling the sample data of fitting the variation probability values using the above four polynomials, as shown in Figure 7. The probability density function of the grey bootstrap sample Y Bootstrap is shown in Figure 8.  The grey bootstrap maximum entropy method is used to obtain the weighted average values, upper-and lower-bound values for the 11 time series under the condition that the confidence level is 0.05, as shown in Figure 9. The grey bootstrap maximum entropy method is used to obtain the weighted average values, upper-and lower-bound values for the 11 time series under the condition that the confidence level is 0.05, as shown in Figure 9.  Figure 9 shows that the values of variation probability calculated using the method in the Reference [34] (The second method) are larger than those calculated using this method. The variation probability has reached 1 for the 10th time series using the second method, which is less in line with the actual degradation process of OVPS of MTSB. Because the interference of random factors during the service process of MTSB was not taken into account.
Based on the Poisson counting process, the total number is calculated as η = 0 for variation probability sample data falling outside the upper-bound value λnU. That is, the reliability of the evaluation results is PR = 100%, which satisfies PR > P.
The estimated uncertainties are calculated for the variation probability of rolling bearing vibration performance, as shown in Table 4.
Based on the Equations (36)-(38), the dynamic average Umean is calculated for the variation probability of time series as Umean = 1.7515/11 = 0.1775.

PMR and PMRR of MTSB
According to the calculated results of variation probability, the estimated true values R0 and upper-and lower-bound intervals [RL, RU] are found as R0 = exp(−λn0), RL = exp(−λnU), RU = exp(−λnL). The range of variation probability λ is [0, 1]. Here, 0 represents the OVPS of bearings without any variation, which is an ideal state with the reliability of bearing performance at 100%. On the other hand, 1 represents that the OVPS of spindle bearings fail completely and the performance is very unreliable. Therefore, if the value of λL is less than 0, let λL = 0 (artificially) in the process of solving the performance maintaining reliability. Based on the Equations (39)-(41), the dynamic evaluation results of the VPMR are shown in Figure 10.  Figure 9 shows that the values of variation probability calculated using the method in the Reference [34] (The second method) are larger than those calculated using this method. The variation probability has reached 1 for the 10th time series using the second method, which is less in line with the actual degradation process of OVPS of MTSB. Because the interference of random factors during the service process of MTSB was not taken into account.
Based on the Poisson counting process, the total number is calculated as η = 0 for variation probability sample data falling outside the upper-bound value λ nU . That is, the reliability of the evaluation results is P R = 100%, which satisfies P R > P.
The estimated uncertainties are calculated for the variation probability of rolling bearing vibration performance, as shown in Table 4.
Based on the Equations (36)-(38), the dynamic average U mean is calculated for the variation probability of time series as U mean = 1.7515/11 = 0.1775.

PMR and PMRR of MTSB
According to the calculated results of variation probability, the estimated true values R 0 and upper-and lower-bound intervals [R L , R U ] are found as R 0 = exp(−λ n0 ), R L = exp(−λ nU ), R U = exp(−λ nL ). The range of variation probability λ is [0, 1]. Here, 0 represents the OVPS of bearings without any variation, which is an ideal state with the reliability of bearing performance at 100%. On the other hand, 1 represents that the OVPS of spindle bearings fail completely and the performance is very unreliable. Therefore, if the value of λ L is less than 0, let λ L = 0 (artificially) in the process of solving the performance maintaining reliability. Based on the Equations (39)-(41), the dynamic evaluation results of the VPMR are shown in Figure 10.
As shown in Figure 10, the values of PMR remain unchanged before the time point corresponding to the third time series. During the period corresponding to the third time series to the fifth time series, the values of PMR have a decreasing trend. The value of PMR reaches the minimum, that is, 81.17442% for the fifth time series. The values of PMR have an increasing trend during the period corresponding to the fifth time series to the eighth time series. The values of PMR have a rapidly decreasing trend during the period corresponding to the eighth time series to the eleventh time series. Figure 10 shows that the values of PMR calculated using the second method are smaller than those calculated using this method. As shown in Figure 10, the values of PMR remain unchanged before the time point corresponding to the third time series. During the period corresponding to the third time series to the fifth time series, the values of PMR have a decreasing trend. The value of PMR reaches the minimum, that is, 81.17442% for the fifth time series. The values of PMR have an increasing trend during the period corresponding to the fifth time series to the eighth time series. The values of PMR have a rapidly decreasing trend during the period corresponding to the eighth time series to the eleventh time series. Figure 10 shows that the values of PMR calculated using the second method are smaller than those calculated using this method.
The estimated true value curve coincides with the actual value curve for the PMR of spindle bearings, but the estimated true value curve is smoother. Moreover, the upperand lower-bound curves envelope the actual value curves fully for the PMR, which verifies the rationality of the evaluation model again.
In order to visualize the difference between the estimated true value and the actual value of the PMR, the error bars of PMR are calculated, as shown in Figure 11.  The estimated true value curve coincides with the actual value curve for the PMR of spindle bearings, but the estimated true value curve is smoother. Moreover, the upper-and lower-bound curves envelope the actual value curves fully for the PMR, which verifies the rationality of the evaluation model again.
In order to visualize the difference between the estimated true value and the actual value of the PMR, the error bars of PMR are calculated, as shown in Figure 11.  Figure 10 shows that the values of PMR calculated using the second method are smaller than those calculated using this method.
The estimated true value curve coincides with the actual value curve for the PMR of spindle bearings, but the estimated true value curve is smoother. Moreover, the upperand lower-bound curves envelope the actual value curves fully for the PMR, which verifies the rationality of the evaluation model again.
In order to visualize the difference between the estimated true value and the actual value of the PMR, the error bars of PMR are calculated, as shown in Figure 11.  As can be seen from Figure 11, the maximum error appears in the eighth time series, but is only 6.55%, and the errors were lower than 1.00% for the first, fourth, sixth, ninth, and tenth time series, which shows that the analysis results of the proposed method have good consistency. By analyzing the vibration data in Figure 3, it can be seen that the vibration performance has a sudden change or a large fluctuation during the period corresponding to these time series. Therefore, it is difficult to accurately estimate the values of PMR during the corresponding period.
According to Equation (42), the dynamic evaluation results of the PMRR are shown in Figure 12. and tenth time series, which shows that the analysis results of the proposed method have good consistency. By analyzing the vibration data in Figure 3, it can be seen that the vibration performance has a sudden change or a large fluctuation during the period corresponding to these time series. Therefore, it is difficult to accurately estimate the values of PMR during the corresponding period.
According to Equation (42), the dynamic evaluation results of the PMRR are shown in Figure 12.
, which shows that the performance states of MTSB is gradually becoming worse, and the degree failure of OVPS is gradually increasing; d(λ10), d(λ11) < −30%, which shows that the performance states of MTSB is worse, and the degree failure of OVPS is very large. Thus, taking appropriate remedial measures are necessary steps before 6773 min, which can avoid serious safety accidents that are caused by the failure of OVPS. Figure 12 also shows that the values of PMRR calculated using the second method are smaller than those calculated using this method. The values of PMRR, d(λ7), d(λ8), d(λ9), d(λ10), d(λ11) < −30%, which shows that the degree failure of OVPS is very large after the corresponding time period of the 7th time series. However, this is less in line with the actual degradation process of OVPS of MTSB. Because the interference of random factors during the service process of MTSB was not taken into account.

Case 2
The test machine and the bearing used in this case are exactly the same as those of Case 1. The vibration data are shown in Figure 13 by changing the test conditions of the motor to a speed of 4000 r/min, an axial load of 4.17 kN, and a radial load of 4.58 kN.  Figure 12 shows that the actual value of PMRR, d(λ 1 ), d(λ 3 ) ≥ 0%, which shows that the performance states of MTSB reaches the optimum and has almost no failure possibility at the corresponding stage; d(λ 2 ), d(λ 7 ), d(λ 8 ) ∈ [−15%, 0%), which shows that the performance states of MTSB is normal, and the degree failure of OVPS is very small; d(λ 4 ), d(λ 5 ), d(λ 6 ), d(λ 9 ) ∈ [−30%, −15%), which shows that the performance states of MTSB is gradually becoming worse, and the degree failure of OVPS is gradually increasing; d(λ 10 ), d(λ 11 ) < −30%, which shows that the performance states of MTSB is worse, and the degree failure of OVPS is very large. Thus, taking appropriate remedial measures are necessary steps before 6773 min, which can avoid serious safety accidents that are caused by the failure of OVPS. Figure 12 also shows that the values of PMRR calculated using the second method are smaller than those calculated using this method. The values of PMRR, d(λ 7 ), d(λ 8 ), d(λ 9 ), d(λ 10 ), d(λ 11 ) < −30%, which shows that the degree failure of OVPS is very large after the corresponding time period of the 7th time series. However, this is less in line with the actual degradation process of OVPS of MTSB. Because the interference of random factors during the service process of MTSB was not taken into account.

Case 2
The test machine and the bearing used in this case are exactly the same as those of Case 1. The vibration data are shown in Figure 13 by changing the test conditions of the motor to a speed of 4000 r/min, an axial load of 4.17 kN, and a radial load of 4.58 kN. From the OVPS of the bearing, there are 5883 vibration signals. Among them, the 334th to 5933rd data points are divided into seven segments at 800 data intervals. The 5934th to 6216th data are separated into one segment, that is, the 8th segment.  From the OVPS of the bearing, there are 5883 vibration signals. Among them, the 334th to 5933rd data points are divided into seven segments at 800 data intervals. The 5934th to 6216th data are separated into one segment, that is, the 8th segment. The probability density estimated truth function f 1 (x) is calculated as shown in Figure 14. Based on the Equations (17) and (18), the numbers are calculated, respectively, that the 8 data samples fall outside the maximum entropy estimated interval of intrinsic time series. According to the Poisson counting principle, the variation frequencies [λ1, λ2, λ3, λ4, λ5, λ6, λ7, λ8, λ9, λ10, λ11] = [0, 0.0136, 0.0738, 0.2450, 0.4563, 0.4370, 0.4213, 0.9963, 1], as shown in Figure 15. Based on the Equations (17) and (18), the numbers are calculated, respectively, that the 8 data samples fall outside the maximum entropy estimated interval of intrinsic time series. According to the Poisson counting principle, the variation frequencies [λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 , λ 7 , λ 8 , λ 9 , λ 10 , λ 11 ] = [0, 0.0136, 0.0738, 0.2450, 0.4563, 0.4370, 0.4213, 0.9963, 1], as shown in Figure 15. As shown in Figure 15, the performance variation probability is non-linear and uncertain-relative to the intrinsic time series. Before the corresponding time interval of the second time series, the variation probability of the bearings vibration performance is almost zero. Between the corresponding time interval of the second time series and the fifth As shown in Figure 15, the performance variation probability is non-linear and uncertain-relative to the intrinsic time series. Before the corresponding time interval of the second time series, the variation probability of the bearings' vibration performance is almost zero. Between the corresponding time interval of the second time series and the fifth time series, the variation probability has an increased trend. Between the corresponding time interval of the fifth time series and the sixth time series, the variation probability has a decreasing trend. Between the corresponding time interval of the sixth time series and the seventh time series, the variation probability has a rapidly increasing trend.

Sequence Number n of Time Series
Variation Probability The closer the correlation coefficient R 2 is to 1, the better the polynomial fitting. As shown in Table 5, the correlation coefficients R 2 for all polynomials are more than 0.8. Therefore, the first-order, second-order, third-order, fourth-order, fifth-order, and sixthorder polynomials will be used in further analysis. Figure 16 shows that the fitted variation probability is almost identical with the actual variation probability for the second and fifth time series. The fitting was performed using the least-square method by the first-order polynomial. With fittings using the least-square method by the second-order polynomial, the fitted variation probability value is basically identical with the actual value for the fourth time series. The fitted variation probability values are almost identical with the actual variation probability values for the first, second and third time series when fitting was performed using the least-square method by the third-order polynomial. When fitting was performed using the least-square method by the fourth-order polynomial, the fitted variation probability value is basically identical with the actual value for the third time series. With fittings using the least-square method by the fifth-order polynomial, the fitted variation probability values are basically identical with the actual values for the fourth, seventh and eighth time series. The fitted variation probability values are almost identical with the actual variation probability values for the first, fifth, sixth, seventh and eighth time series when fitting was performed using the least-square method by the sixth-order polynomial. In short, the above six polynomials have their own advantages in fitting. Therefore, relevant information should be fully fused to effectively monitor the evolution process of rolling bearing vibration performance.

Uncertainty of Variation Probability of MTSB (Case 2)
The variation probability of the first time series is taken as an example here. According to the grey bootstrap method, the number of bootstrap re-sampling is taken to be B = 20,000, and confidence level is taken to be P = 95%. Based on the Equations (20)- (32), the grey bootstrap sample Y Bootstrap is obtained by sampling the sample data of fitting the variation probability values using the above six polynomials, as shown in Figure 17. and third time series when fitting was performed using the least-square method by the third-order polynomial. When fitting was performed using the least-square method by the fourth-order polynomial, the fitted variation probability value is basically identical with the actual value for the third time series. With fittings using the least-square method by the fifth-order polynomial, the fitted variation probability values are basically identical with the actual values for the fourth, seventh and eighth time series. The fitted variation probability values are almost identical with the actual variation probability values for the first, fifth, sixth, seventh and eighth time series when fitting was performed using the least-square method by the sixth-order polynomial. In short, the above six polynomials have their own advantages in fitting. Therefore, relevant information should be fully fused to effectively monitor the evolution process of rolling bearing vibration performance.

Uncertainty of Variation Probability of MTSB (Case 2)
The variation probability of the first time series is taken as an example here. According to the grey bootstrap method, the number of bootstrap re-sampling is taken to be B = 20,000, and confidence level is taken to be P = 95%. Based on the Equations (20)- (32), the grey bootstrap sample YBootstrap is obtained by sampling the sample data of fitting the variation probability values using the above six polynomials, as shown in Figure 17. The probability density function of the grey bootstrap sample YBootstrap is shown in Figure 18. The probability density function of the grey bootstrap sample Y Bootstrap is shown in Figure 18. The grey bootstrap maximum entropy method is used to obtain the weighted average values, upper-and lower-bound values for the 8 time series under the condition that the confidence level is 0.05, as shown in Figure 19. The grey bootstrap maximum entropy method is used to obtain the weighted average values, upper-and lower-bound values for the 8 time series under the condition that the confidence level is 0.05, as shown in Figure 19.  Figure 19 shows that the values of variation probability calculated using the method in the Reference [34] (The second method) are larger than those calculated using this method. The variation probability has reached 1 for the 7th time series using the second method, which is less in line with the actual degradation process of OVPS of MTSB. Because the interference of random factors during the service process of MTSB was not taken into account.
Based on the Poisson counting process, the total number is calculated as η = 0 for variation probability sample data falling outside the estimated upper-bound curve λnU. That is, the reliability of the evaluation results is PR = 100%, which satisfies PR > P.
The estimated uncertainty is calculated for the variation probability of rolling bearing vibration performance, as shown in Table 7.   Figure 19 shows that the values of variation probability calculated using the method in the Reference [34] (The second method) are larger than those calculated using this method. The variation probability has reached 1 for the 7th time series using the second method, which is less in line with the actual degradation process of OVPS of MTSB. Because the interference of random factors during the service process of MTSB was not taken into account.
Based on the Poisson counting process, the total number is calculated as η = 0 for variation probability sample data falling outside the estimated upper-bound curve λ nU . That is, the reliability of the evaluation results is P R = 100%, which satisfies P R > P.
The estimated uncertainty is calculated for the variation probability of rolling bearing vibration performance, as shown in Table 7. According to the calculated results of variation probability, the estimated true values R 0 and upper-and lower-bound intervals [R L , R U ] are found as R 0 = exp(−λ n0 ), R L = exp(−λ nU ), R U = exp(−λ nL ). The range of variation probability λ is [0, 1]. Here, 0 represents the OVPS of bearings without any variation, which is an ideal state with the reliability of bearing performance at 100%. On the other hand, 1 represents that the OVPS of spindle bearings fail completely and the performance is very unreliable. Therefore, if the value of λ L is less than 0, let λ L = 0 (artificially) in the process of solving the performance maintaining reliability. Based on the Equations (39)-(41), the dynamic evaluation results of the VPMR are shown in Figure 20. As shown in Figure 20, the values of VPMR remain unchanged before the time point corresponding to the second time series. During the period corresponding to the second time series to the fifth time series, the values of VPMR have a decreasing trend. The value of VPMR reaches the minimum, that is, 63.36544% for the fifth time series. The values of VPMR have an increasing trend during the period corresponding to the fifth time series to the sixth time series. The values of VPMR have a rapidly decreasing trend during the period corresponding to the sixth time series to the seventh time series. Figure 20 shows that the values of PMR calculated using the second method are smaller than those calculated using this method.
The estimated true value curve coincides with the actual value curve for the VPMR of spindle bearings, but the estimated true value curve is smoother. Moreover, the upperand lower-bound curves envelope the actual value curves fully for the VPMR, which verifies the rationality of the evaluation model again.
In order to visualize the difference between the estimated true value and the actual value of the VPMR, the error bars of PMR are calculated, as shown in Figure 21.  As shown in Figure 20, the values of VPMR remain unchanged before the time point corresponding to the second time series. During the period corresponding to the second time series to the fifth time series, the values of VPMR have a decreasing trend. The value of VPMR reaches the minimum, that is, 63.36544% for the fifth time series. The values of VPMR have an increasing trend during the period corresponding to the fifth time series to the sixth time series. The values of VPMR have a rapidly decreasing trend during the period corresponding to the sixth time series to the seventh time series. Figure 20 shows that the values of PMR calculated using the second method are smaller than those calculated using this method.
The estimated true value curve coincides with the actual value curve for the VPMR of spindle bearings, but the estimated true value curve is smoother. Moreover, the upper-and lower-bound curves envelope the actual value curves fully for the VPMR, which verifies the rationality of the evaluation model again.
In order to visualize the difference between the estimated true value and the actual value of the VPMR, the error bars of PMR are calculated, as shown in Figure 21. As can be seen from Figure 21, the maximum error appears in the sixth time series, but is only 9.91%, and the errors were lower than 0.05% for the first, second, fourth and eighth time series, which shows that the analysis results of the proposed method have good consistency. By analyzing the vibration data in Figure 13, it can be seen that the vibration performance has a sudden change or a large fluctuation during the period corresponding to these time series. Therefore, it is difficult to accurately estimate the values of VPMR during the corresponding period.
According to Equation (42), the dynamic evaluation results of the PMRR are shown in Figure 22. As can be seen from Figure 21, the maximum error appears in the sixth time series, but is only 9.91%, and the errors were lower than 0.05% for the first, second, fourth and eighth time series, which shows that the analysis results of the proposed method have good consistency. By analyzing the vibration data in Figure 13, it can be seen that the vibration performance has a sudden change or a large fluctuation during the period corresponding to these time series. Therefore, it is difficult to accurately estimate the values of VPMR during the corresponding period.
According to Equation (42), the dynamic evaluation results of the PMRR are shown in Figure 22.  Figure 22 shows that the actual value of PMRR, d(λ1) ≥ 0%, which shows that the performance states of MTSB reaches the optimum and has almost no failure possibility at the corresponding stage; d(λ2), d(λ3) ∈ [−15%, 0%), which shows that the performance states of MTSB is normal, and the degree failure of OVPS is very small; d(λ4), d(λ5), d(λ6) ∈ [−30%, −15%), which shows that the performance states of MTSB is gradually becoming worse, and the degree failure of OVPS is gradually increasing; d(λ7), d(λ8) < −30%, which shows that the performance states of MTSB is worse, and the degree failure of OVPS is very large. Thus, taking appropriate remedial measures are necessary steps before 5134 min, which can avoid serious safety accidents that are caused by the failure of OVPS. Figure 22 also shows that the values of PMRR calculated using the second method are smaller than those calculated using this method. The values of PMRR, d(λ4), d(λ5), d(λ6), d(λ7), d(λ8) < −30%, which shows that the degree failure of OVPS is very large after the corresponding time period of the 4th time series. However, this is less in line with the actual degradation process of OVPS of MTSB. Because the interference of random factors  Figure 22 shows that the actual value of PMRR, d(λ 1 ) ≥ 0%, which shows that the performance states of MTSB reaches the optimum and has almost no failure possibility at the corresponding stage; d(λ 2 ), d(λ 3 ) ∈ [−15%, 0%), which shows that the performance states of MTSB is normal, and the degree failure of OVPS is very small; d(λ 4 ), d(λ 5 ), d(λ 6 ) ∈ [−30%, −15%), which shows that the performance states of MTSB is gradually becoming worse, and the degree failure of OVPS is gradually increasing; d(λ 7 ), d(λ 8 ) < −30%, which shows that the performance states of MTSB is worse, and the degree failure of OVPS is very large. Thus, taking appropriate remedial measures are necessary steps before 5134 min, which can avoid serious safety accidents that are caused by the failure of OVPS. Figure 22 also shows that the values of PMRR calculated using the second method are smaller than those calculated using this method. The values of PMRR, d(λ 4 ), d(λ 5 ), d(λ 6 ), d(λ 7 ), d(λ 8 ) < −30%, which shows that the degree failure of OVPS is very large after the corresponding time period of the 4th time series. However, this is less in line with the actual degradation process of OVPS of MTSB. Because the interference of random factors during the service process of machine tool spindle bearings (MTSB) was not taken into account.
In summary, the model realizes on-line monitoring of the degradation process of OVPS, which can give timely feedback so as to take preventive and remedial measures before the failure of OVPS for MTSB. Compared to the method in the Reference [8], the vibration threshold does not need to be set manually. Compared to other AI prediction methods, it does not require training on vibration signals, and the adjustment process of parameter is also omitted. In addition, there is no requirement for the length of time series. It is worth noting that, if there is an amount of missing data in the process of signal acquisition, the evaluation results of this approach may be inaccurate, and it cannot diagnose which component of bearings the fault occurred on.

Conclusions
By means of the Poisson process theory and taking variation probability as a time variable, the proposed model can effectively realize the dynamic evaluation for the degradation process of the OVPS for MTSB. This provides a technical and theoretical basis for on-line health detection and fault diagnosis of spindle bearings.

•
The variation probability, obtained using the maximum entropy method and the Poisson counting principle, can accurately describe the degradation information and evolution process of the OVPS of MTSB.

•
Considering the interference of random factors, the least-squares method by polynomial fitting, fused into the grey bootstrap maximum entropy method, can be used to calculate the dynamic mean uncertainty, so as to evaluate the random fluctuation state of OVPS.

•
The results show the maximum relative errors between the estimated true value and the actual value of the PMR are 6.55% and 9.91% for the MTSB in the two studied cases: Case 1 and Case 2, respectively. Appropriate remedial measures should be taken before 6773min and 5134 min for the MTSB in the two studied cases: Case 1 and Case 2, respectively, which can avoid serious safety accidents caused by the failure of OVPS. Institutional Review Board Statement: The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest:
The authors declare that there is no conflict of interest regarding the publication of this paper.

Nomenclature
x(k) kth performance data in the intrinsic sequence. k order number of performance data in intrinsic sequence. N total number of performance data in the intrinsic sequence.  α significant level. P confidence level.

N n1
Number showing that performance data are less than X L1 for the nth time series.

N n2
Number showing that performance data are more than X U1 for the nth time series. G q (λ) qth order polynomial. q order number of polynomial function. p qγ coefficient of the power function λ γ . Y(n) data sample of variation probability of OVPS for the n time series. y n (u) uth data in the variation probability data sample for the n time series. V β βth bootstrap re-sampling sample. B times of the bootstrap re-sampling. v β (Θ) Θth data in the βth bootstrap re-sampling sample. u time variable. c 1 ; c 2 coefficients to be estimated. λ nL ; λ nU lower-bound value and upper-bound value of the variation probability data sample sample for the nth time series. U λn estimated uncertainty of variation probability. P R reliability of the polynomials fitting effect using least-squares method. U mean dynamic average uncertainty. | PR = 100% calculation process is under the condition of P R = 100%. e number of occurring failure events. Q probability of failure events occurring e times. R(λ n ) function of the variation probability λ n .