Modeling the Relationship between Vibration Features and Condition Parameters Using Relevance Vector Machines for Health Monitoring of Rolling Element Bearings under Varying Operation Conditions

Rotational speed and load usually change when rotating machinery works. Both this kind of changing operational conditions and machine fault could make the mechanical vibration characteristics change. Therefore, effective health monitoring method for rotating machinery must be able to adjust during the change of operational conditions. This paper presents an adaptive threshold model for the health monitoring of bearings under changing operational conditions. Relevance vector machines (RVMs) are used for regression of the relationships between the adaptive parameters of the threshold model and the statistical characteristics of vibration features. The adaptive threshold model is constructed based on these relationships. The health status of bearings can be indicated via detecting whether vibration features exceed the adaptive threshold. This method is validated on bearings running at changing speeds. The monitoring results show that this method is effective as long as the rotational speed is higher than a relative small value.


Introduction
Rolling element bearing is a fundamental component of rotating machinery.Its fault is one of the most important causes for the failure of rotating machinery.Health monitoring of bearings aims at finding faults at their early stage which allows maintenance actions to be taken in time to prevent machinery from the failures.Conventional approaches are based on the assumption that the effect of bearing faults can be considered as load fluctuations in the time associated with the impact of rolling elements when passing through the defective zone [1].These load fluctuations give periodic excitations of structural resonances as long as the rotational speed is constant.These periodic responses are named as the components of fault feature frequencies [1,2].Traditional frequency based methods have been developed to extract the components of fault feature frequencies.
Rotating machine usually operates under varying load or speed conditions.For example, the output power and speed of an engine have to change with a large range to complete the slide, climbing, and landing process of an aircraft.In the cases where a machine operates under varying speeds or loads, its dynamics and vibrations exhibit nonstationary characteristics.Both the amplitudes and the frequencies of the vibration change with time [1,3].Thus traditional frequency based methods are not applicable to such cases for health monitoring.
Great research efforts have been concentrated on nonstationary signal processing methods for health monitoring and fault diagnosis of rotating machinery.These methods include short time Fourier transform (STFT) [4], Wigner-Ville distribution (WVD), wavelet analysis [5][6][7], and empirical mode decomposition (EMD) [8,9].The STFT is to move a short time window along the vibration record and obtain the Fourier spectrum as a function of time shift.However, the frequency resolution is the reciprocal of the effective time window length; thus STFT could not give higher resolution both in time domain and frequency domain.The WVD obtains better local characteristics than STFT.But WVD is a nonlinear transform method and suffers from cross terms between the actual components.The wavelet analysis is to decompose the signal in terms of a family of "wavelets." The wavelet analysis possesses the multiscale character and gives a better time localization at high frequencies, and for that reason it has been widely applied to nonstationary signals analysis.But unfortunately, the wavelet analysis suffers from energy leakage in the vicinity of frequencies of interest due to the limited length of wavelet base function.The EMD decomposes the signal into a finite number of intrinsic mode functions (IMFs), which reflect the intrinsic of the signal.Therefore, the EMD method has been widely applied to nonstationary signal processing.However, the EMD also has problems, such as the end effect, the disturbances due to pseudo-information at low frequencies, and the IMF criterion.The order tracking method uses multiples of running speed, instead of absolute frequencies, as the frequency base.These multiples of running speed are named orders.Compared with other nonstationary signal processing methods, the order tracking method has advantages for analysis of speed-related vibration components, especially when rotational speed changes over a greater range [3,10].However, angle domain signals are used for order tracking.And to obtain the angle domain signals, either additional hardware is needed to sample the vibration at constant increments of shaft angle, or additional software is needed to resample the vibration from time domain signals.
These signal processing methods provide ways to analyze vibration signal under varying operation conditions.Features that are extracted by these methods are assumed to vary with health status.But these nonstationary signal processing methods can be unavailable in some cases.For example, conditions could change sharply in short time, and the amplitudes of features extracted could be changed with both machine fault and the varying operation conditions.Thus how to recognize whether the changes of features are caused by faults or the varying operation conditions is difficult.
This paper aims at solving the health monitoring problem of bearings under varying operation conditions with pattern recognition method.We consider root mean square (RMS)  R of the overall vibration, which can be computed easily, as a feature to indicate fault severity of bearings.RMS is a measurement of the effective energy content in a vibration signal.As vibration level generally increases when faults occur, RMS is well used for diagnostics and prognostics of rotating machines [11][12][13].The proposed method is based on the relation functions between vibration features and condition parameters.The relationship functions are obtained with relevance vector machines (RVMs) and are used as adaptive parameters of threshold model.Since the adaptive parameters are continuous functions of condition parameters, the threshold model obtained is available at any operation conditions.
Relevance vector machines (RVMs) and support vector machines (SVMs) have been widely used in condition monitoring and fault diagnosis/prognosis [14][15][16].But in previous works, RVMs/SVMs are mainly used for classification or regression.When RVMs/SVMs are used for classification, the inputs of RVMs/SVMs are vibration features and the outputs are pattern indexes.When RVMs/SVMs are used for regression, the inputs are vibration features and outputs are assessment of fault degree.The novelty of this paper is that RVMs are used to model the relationship between condition parameters and the statistics of vibration feature.The input of RVMs is condition parameter and the outputs are statistics of vibration feature.These statistics are used as adaptive parameters of fault monitoring model.And these adaptive parameters make the monitoring model adaptive and available under varying operation conditions.

Methodology
2.1.Adaptive Threshold Model.We select Gaussian threshold model (GTM) as the health monitoring model.The GTM is derived from Chebyshev's inequality [17].For a feature  that has Gaussian distribution, if its mean and standard deviation are  and , respectively, the following Chebyshev's inequality holds: ( The inequality means that the probability of  ∈ ( − ,  + ) is bigger than 1 −  −2 . is called threshold tolerance.The above inequality is called three-sigma criterion when  = 3.This model is also effective sometimes for features that are not Gaussian as long as the threshold tolerance is well optimized.According to Chebyshev's inequality, for a given threshold tolerance , the threshold interval of vibration RMS  R can be constructed as  = [ R −  R ,  R +  R ].As we are concerned about the upper bond only, a decision function can be obtained as ( R ) =  R +  R −  R , which allows the monitoring to be implemented explicitly.( R ) is considered as the health index that indicates the health status of bearings; in particular, ( R ) ≥ 0 means that the vibration level  R is normal and the bearing is healthy, whereas ( R ) < 0 means that  R is abnormal and the bearing is faulty.
Here only rotational speed is addressed to simplify the problem.As the varying rotational speeds can also result in the change of vibration energy level, the adaptive GTM should be the function of both RMS  R and the rotational speed .Here we define the detecting function of the adaptive GTM as where  R () and  R () are, respectively, mean and standard deviation of the RMS at rotational speed  and  R () +  R () is the threshold of vibration level at rotational speed .
In order to apply the GTM to data at any rotational speeds, (2) should be a continuous function of , which can be ensured as long as  R () and  R () are continuous functions.Two RVMs are used to fit  R () and  R ().The first RVM, which is trained on samples {(),  R () |  = 1, 2, . . ., }, gives the rotational-speed-varying parameters  R ().() is the rotational speed of the th sample. R () is the mean of vibration RMS at rotational speed (). is the number of training samples.The second RVM, which is trained on samples {(),  R ()}, gives the rotational-speedvarying parameters  R (). R () is the standard deviation of vibration RMS at rotational speed ().These training samples are extracted from tests of normal bearings.

Estimation of Rotational Speed and Vibration Level.
The raw rotational speed signals   collected with the optical encoder are series of pulses.One pulse is generated for each revolution.We search the raw data for points that satisfy   () ≤ 0 and   ( + 1) ≥ 0. Then the zero-crossing time   can be estimated with linear interpolation by where () and ( + 1) are the time instants for   () and   ( + 1), respectively.Then the speed   at time   can be estimated by   = 60/(  −  −1 ), in which   is the th zero-crossing time and  −1 is the ( − 1)th zero-crossing time.
The vibration RMS  R  at time   is estimated with the vibration points that are sampled in the time interval ( −1 ,   ].

Statistical Characteristics Estimation of RMS.
We suppose {  } and { R  } to be the rotational speed series and the RMS series that are obtained from the th test.It should be noted that the series {  } of different tests are different.In order to compute the statistical parameters of RMS at different rotational speeds, we need to resample { R  } to make them share a common rotational speed series {() |  = 1, 2, . . ., }.The series () has equal speed increments.
The corresponding RMS series { R  ()} of the th test can be yielded with linear interpolation.We search the series {  } for the points that satisfy (  − ())( (+1) − ()) ≤ 0. The corresponding RMS can be estimated by The RMS series { R  ()} yielded shares a common rotational speed coordinate.
We carry out  times of tests and yield  series { R  () |  = 1, 2, . . ., }.The mean series and the standard deviation series of RMS can be yielded by (5) . ., },  detectors will be generated and these detectors are only applicable to  R whose corresponding rotational speed is ().Such detectors are not practical.In order to solve this problem, we introduce two RVMs to generate the detector of (2).The first RVM is trained on the set {(),  R ()} and gives continuous function  R ().The second one is trained on the set {(),  R ()} and gives continuous function  R ().Putting  R () and  R () into (2) yields the adaptive GTM, which is applicable at any rotational speed.A brief introduction about RVM regression will be given in Section 2.4.

Relevance Vector Machines for Regression.
RVMs models were originally proposed by Tipping [18,19].RVMs share the functional form with the popular support vector machines (SVMs).RVMs exploit a probabilistic Bayesian learning framework and utilize dramatically fewer kernel functions than comparable SVMs.RVMs have a number of additional advantages.These advantages include the benefits of probabilistic predictions, automatic estimation of parameters, and the facility to utilize arbitrary kernel functions.In this section, a brief outline of the RVMs theory is present for applying it to the bearing datasets.The details of RVM can be found in [18,19].Given a dataset of training pairs {x  ,   }  =1 , where x  is the input feature vector,   is the target output of x  , and  is the number of samples.RVMs make predictions based on the model of ( 6) that shares a functional form with SVMs.Consider where w = [ 1 , . . .,   ] T is the weight vector, K = [(x  , x 1 ), . . ., (x  , x  )] T is the kernel function vector, and  0 is the offset.Gaussian kernel function is generally selected.Consider where  is the basis width of the kernel.For a test sample x  whose output target   is unknown, the trained RVMs give ỹ as the predicted value of   .  is class label in cases where RVMs are used for classification, and   is real value in cases where RVMs are used for regression.
In order to overcome overfitting, RVMs adopt a Bayesian perspective, and zero-mean Gaussian prior distribution is used to introduce an individual hyperparameter   on each weight   .The prior distribution of weights then can be written as where  = [ 0 , . . .,   ].The posterior distribution over the weights is then obtained from Bayes' rule: with Σ = (Φ T BΦ + A) −1 ,  = ΣΦ T By, where A = diag( 0 ,  1 , . . .,   ) and B =  −2 I  .By integrating out the weights, the marginal likelihood can be obtained for the hyperparameters: Tipping released a Matlab software package named Sparse-Bayes Version 2.0.The software computes the posterior distribution over the hyperparameters  and returns their most-probable values by maximizing the marginal likelihood function equation (11).As many of the weights are mathematically set to zero, a lower number of kernel functions are used to construct the predictive model.This results in a sparse model, which can generalize well and is fast to compute.The pairs of data (x  ,   ) with nonzero weights are called relevance vectors.2.5.Discussions.The proposed health monitoring method concerns the selection of two parameters, in particular the threshold tolerance  in (2) and the basis width  in Gaussian kernel equation (7).The threshold tolerance  is determinative to balance false negative and false positive.According to (2), the larger the  is, the wider the threshold interval is and hence the smaller the probability that the normal samples exceed the threshold intervals, which means that the system has fewer false negatives.Conversely, the smaller the  is, the narrower the threshold interval is and the smaller the probability that the novel samples are located in the threshold intervals, which means that the system is more sensitive.So  can be set with cross-validation as follows: (1) initialize  to be 3 or any value; (2) if monitoring results of historical data show unacceptable false negative alarms, increase ; (3) on the contrary, if monitoring results of historical data show unacceptable false positive rate, reduce ; (4) false negative rate and false positive rate cannot be reduced simultaneously, and a balance between them should be reached [17].
The basis width  is a determinative parameter for RVMs.The larger the basis width is, the sparser the RVMs are and fewer points are relevance vectors.On the contrary, the smaller the basis width is, the more the relevance vectors will be selected, and the RVMs have worse generalization.Choosing of basis width is experiential.Cross-validation results can be referred to for choosing of basis width.

Monitoring Model Construction.
The proposed method was applied to a test rig of machinery fault simulation, as shown in Figure 1.In order to compute the mean series { R ()} and the standard deviation series { R ()} as training samples, 20 tests of healthy bearings were carried out on the test rig.The rotational speed interval of the motor is set to be [200∼1600] rpm.The vibration signals and rotational speed signals were collected synchronously.And vibration RMS and rotational speeds were computed for each revolution.Figure 2 shows the rotational speeds and the vibration RMS values of a healthy bearing obtained by a test.It can be seen from Figure 2 that although the bearing is healthy, vibration energy level changes greatly at different rotational speeds.At the same time, the fault status may only lead to a little change in system vibration response.Thus it is difficult to identify whether the change of vibration energy is caused by faults or varying operation conditions.As mentioned in Section 2.2.3, different tests generate different rotational speed series.In order to compute the statistical mean and standard deviation of RMS at different rotational speeds, a new rotational speed series with equal speed increments is created first, and then the RMS series corresponding to the new speed series is estimated for each test.As the rotational speeds range from 200 rpm to 1600 rpm, the new speed series will be {(),  = 1, 2, . . ., } = {203, 209, 215, . . ., 1397} when we set  = 200.The RMS series of each test corresponding to the rotational speed series can be estimated with the linear interpolation method introduced in Section 2.2.3.
Selection of the threshold tolerance  using cross-validation requires fault samples for computation of false positive rates.In order to generate fault samples, tests of healthy bearings and bearings with seeded inner race defects and outer race defects were carried out, as shown in Figure 1.The bearings with seeded defects were installed in the position of bearing I.The vibration RMS values and rotational speeds of testing samples were also computed with each circle the machine rotates.
20 RMS series { R  (),  = 1, 2, . . ., 20,  = 1, 2, . . ., 200} were generated from 20 tests of healthy bearings.Then the mean series { R ()} and standard deviation series { R ()} were estimated with (5).{(),  R ()} and {(),  R ()} are training samples of RVMs.The parameter  R () of the adaptive GTM was yielded by training a RVM with {(),  R ()}, and  R () was yielded by training another RVM with {(),  R ()}.In order to validate the regression results and to select the basis width  with cross-validation, we need to generate testing samples of RVMs.We generated a new rotational speed series {()} = {206, 212, 218, . . ., 1394}, and the corresponding mean series { R ()} and standard deviation series { R ()} were yielded in the same way as above.{(),  R ()} and {(),  R ()} are testing samples of RVMs.Selection of the basis width  is vital for RVM regression.For cross-validation of , we define training error  1 and testing error  2 of RVM, respectively, as where    is the real target of a training sample, ỹ  is the estimated output given by RVM,    is the real target of a testing sample, and ỹ  is the estimated output given by RVM.The cross-validation results of the first RVM, which is used for regression of  R (), are shown in Figure 3. Figure 3(a) illustrates the training error and testing error with varying basis width in the interval [2,70].It can be seen that both the training error and the testing error decrease dramatically with the increase of  at the beginning.If  becomes bigger than 6, the training error and the testing error turn out to increase slowly with the increase of . Figure 3(b) illustrates the number of relevance vectors.It can be seen that the number of relevance vectors decreases with the increase of .We can infer from Figure 3 that when  takes values in the interval [30, 70], both the testing error and the training error are smaller than 0.09, and the model is relatively sparse as the number of relevance vectors is smaller than 38.We set the basis width to be  = 48; in such case, both the training error and the testing error are 0.06, and the number of relevance vectors is 24.The cross-validation results of the second RVM, which is used for regression of  R (), are shown in Figure 4. Similarly, we can infer from Figure 4 that  can take values in the interval [20,70].Once  = 28 is chosen, both the training error and the testing error are 0.06, and the number of relevance vectors is 34.
The regression and test results of the parameters  R () and  R () are, respectively, shown in Figures 5 and 6.It can be seen from Figure 5 that vibration RMS values have a general uptrend when rotational speeds increase.Vibration level reaches the local peak when the rotational speeds are around 4 particular values, 257, 527, 1037, and 1349 rpm.The regression precision of  R () is high as few testing samples or training samples are far away from the regression curve.It can be seen from Figure 6 that the standard deviation of vibration RMS is also higher when the rotational speeds are around the 4 particular values.The adaptive parameter  R () is also well fitted.There are only several testing samples and training samples that are far away from the regression curve, and these samples are around the 4 particular rotational speeds.
We let the threshold tolerance  change in the interval [0.5, 10] and apply the corresponding detectors to testing samples.The false negative rates and the false positive rates yielded by detectors with different threshold tolerance are shown in Figure 7.It can be seen that false negative rate decreases with the increase of .The false negative rate is close to 0.04 when  increases to 3.5.The false positive rate increases with the increase of .The false positive rate is 0.21 when  = 3.5.As high false negative rate is generally more intolerable than high false positive rate for field engineers, a relative bigger  = 3.5 is an optimized selection.Figure 9 shows the monitoring results of a healthy bearing during Test 19.Test 19 lasted 59.1 seconds, during which the rotational speeds increased from 200 rpm to 1400 rpm.The number of the extracted samples is 669.There is only a sample whose RMS value exceeds its corresponding thresholds as shown in Figure 9(a).The sample is at around 282 rpm.The health index of the sample is negative as shown in Figure 9(b).This sample is false negative.The false negative rate is 1/669 ≈ 0.002.

Monitoring Results.
Figure 10 shows the monitoring results of a faulty bearing during Test 14.The rotational speeds decreased from 1271 rpm to 200 rpm in 55.7 seconds.It can be seen that  there are lots of samples that take RMS values smaller than their corresponding thresholds.The health indexes of these samples are positive as shown in Figure 10(b).These samples are false positives.It should be noted that most of these false positives are in case where the rotational speeds are lower than 570 rpm.But there are only 2 false positives when the rotational speeds are higher than 570 rpm.These 2 false positives are at around 990 rpm.The total number of the samples whose speeds are higher than 570 rpm is 521.Thus the false positive rate is 2/521 ≈ 0.004 when rotational speeds take values in the interval [570 1400] rpm.
Figure 11 shows the monitoring results of a faulty bearing during Test 16.The rotational speeds of the tested bearing increased from 200 rpm to 1400 rpm in 68.2 seconds.It can be seen that there are lots of false positives where the rotational speeds are lower than 558 rpm.Most samples exceed their thresholds when the rotational speeds are higher than 558 rpm except 2 samples at around 1009 rpm.The total number of the samples whose speeds are higher than 558 rpm is 502.Thus the false positive rate is 2/502 ≈ 0.004 when rotational speeds take values in the interval [558 1400] rpm.
The above monitoring results show that the proposed method has small false negatives and false positives are focused on the low rotational speed cases.Monitoring results of other tests show similar results.It is safe to reach the conclusion according to all the monitoring results that the proposed method is effective for the bearings of the test rig when the rotational speeds are beyond 600 rpm.

Discussions.
As introduced in Section 1, the effects of bearing faults can be mainly considered as additive load fluctuations that are caused by the impact of rolling elements

Figure 1 :
Figure 1: Test rig of machinery fault simulation and bearings with defects.

Figure 2 :
Figure 2: Varying of rotational speeds and vibration level of a healthy bearing: (a) rotational speeds, (b) vibration RMS.

Figure 3 :Figure 4 :
Figure 3: Regression characteristics of  R () with different basis width: (a) regression error; (b) number of relevance vectors.

Figure 8
shows the monitoring results of a healthy bearing during Test 18.The rotational speeds of the tested bearing decreased from 1400 rpm to 200 rpm in 20.0 seconds.The total number of the extracted samples is 215.

Figure 8 (
a) describes the vibration RMS values and their thresholds at different time and rotational speeds.

Figure 8 (
b) shows the health indexes at different rotational speeds.There are 8 samples whose RMS values exceed their thresholds.Two of these 8 samples are at around 471 rpm, another 5 are at the interval [941 993] rpm, and the other one is at round 1202 rpm.These 8 samples, whose health indexes are negative as shown in Figure 8(b), are false negatives.The false negative rate is 8/215 ≈ 0.037.

Figure 5 :
Figure 5: Regression and test results of parameter  R (): (a) regression results; (b) testing results.

Figure 6 :
Figure 6: Regression and test results of parameter  R (): (a) regression results; (b) testing results.

Figure 7 :Figure 8 :Figure 9 :
Figure 7: Change of test error at different threshold tolerance.

Figure 10 :Figure 11 :
Figure 10: The monitoring results of a faulty bearing during Test 14.
Data Preparation 2.2.1.Raw Data Collection.We adjust the rotational speed of healthy bearings in an interval [ min ,  max ], specifically, we first raise the rotational speed from  min to  max and then lower it from  max to  min .During the operation, vibration signals V  and rotational speed signals   are sampled synchronously, and the subscript  denotes that the signals are yielded during the th test,  = 1, 2, . . ., .Vibration signals are measured with accelerometers, and the rotational speed signals are measured with an optical encoder.V  and   are the raw data required in the proposed method.