Rolling element bearing diagnosis based on probability box theory

Journal Pre-proof

Rolling Element Bearing Diagnosis Based on Probability Box Theory Hong Tang , Yi Du , Hong-Liang Dai PII: DOI: Reference:

S0307-904X(19)30671-7 https://doi.org/10.1016/j.apm.2019.10.068 APM 13136

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

13 January 2019 21 October 2019 30 October 2019

Please cite this article as: Hong Tang , Yi Du , Hong-Liang Dai , Rolling Element Bearing Diagnosis Based on Probability Box Theory, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.10.068

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights

(a) Decoupling of the bearing composite failure mechanism has been

avoided.

(b) The amount of statistical information rich in the raw signal has been

utilized.

(c) The sampling uncertainty caused by the average multi-segment signal and

the fluctuation of the signal with time have been solved.

(d) Subjective cognitive uncertainty and objective experimental uncertainty

of experimental data have been considered.

1

Title: Rolling Element Bearing Diagnosis Based on Probability Box Theory Authors: Hong Tang1, Yi Du*2 and Hong-Liang Dai1

Author affiliations: 1. State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha, 410082, China 2. City College, Kunming University of Science and Technology，Kunming，650051, China

Corresponding authors: Yi Du Tel.: +86 731 88664002 Fax: +86 731 88711911 Email: [email protected]

2

Abstract: Feature extraction leads to the loss of statistical information of raw data and ignores the sampling uncertainty and the fluctuations in the signal over time in mechanical fault diagnosis. In this paper, novel modeling methods for mechanical signals based on probability box theory were proposed to solve the above problem. First, the type of random distribution of the bearing signals were analyzed. Then, a Dempster-Shafer structure was obtained to establish a probability box model. To address the identification difficulty of the type of random distribution for the bearing signals, a second probability box model was established based on a vector consisting of features from the bearing signals. If the data are not found to follow a random distribution, a third modeling method based on the definition of probability boxes was proposed. The effectiveness and applicability of the three proposed models were compared with experimental data from rolling element bearings. The combination of probability box theory and mechanical fault diagnosis theory can open up a new research direction for mechanical fault diagnosis. Keywords: Compound fault detection; Feature extraction; Probability box theory; Rolling element bearing.

3

1 Introduction Bearings are an integral part of rotating machinery, and bearing fault diagnosis has been studied for many years [1]. The complexity of bearing faults is extremely nonlinear. The key to diagnosing bearing operating status is to extract overall feature information from faults [2]. There exist various classical methods for the signal processing of raw bearing signals. By using simulated bearing fault signals, Ho and Randall [3] presented an efficient application of self-adaptive noise cancellation in conjunction with envelope analysis to remove discrete frequency masking signals. For diagnostic reasons, the effects of different preprocessing techniques, such as high-pass or bandpass filtration, envelope detection (demodulation) and wavelet transform of the vibration signals, prior to feature extraction were investigated by Lou and Loparo [4]. The relationship between vibration and current frequencies caused by incipient bearing failures was described by Samanta and Al-Balushi [5] using artificial neural network-based fault diagnostics. Schoen et al. [6] took the initial step of investigating the efficacy of current monitoring for bearing fault detection. Fast Fourier transforms with envelope detection were shown to be efficient in detecting some types of bearing faults by Tse et al. [7]. In addition, some new nonlinear dynamics-based theories and methods have emerged for performing feature extraction of bearings [8-9]. An example of surface defects, preloads, and radial clearance was given by Cao and Xiao [10] using comprehensive dynamic model studies. Harsha [11] presented an analytical model to investigate the structural vibrations of a high-speed rotor supported by rolling element bearings. The nonlinear dynamic behavior of rolling element bearing systems, including surface defects, was analyzed by Rafsanjani [12]. Recent 4

theoretical developments have revealed that combinations of algorithms for feature extraction can make the information extracted with those algorithms complementary. Lei et al. [13] presented a novel method for fault diagnosis based on empirical mode decomposition that consisted of an improved distance evaluation technique and a combination of multiple adaptive neuro-fuzzy inference systems. Optimal resonance demodulation using a combination of a fast kurtogram for initial estimates and a genetic algorithm for final optimization was given by Zhang [14] using a model and algorithm to design the parameters. To develop a health assessment of the wheelset bearings of high-speed trains, a modeling method for the safe region concept was presented by Liu et al. [15]. The combinations of algorithms are diverse and can compensate for the defects of a single algorithm [16-17]. However, the methods in the above references are not without their own problems in the vast majority of cases, which can be categorized conclusively as follows: (1) Feature extraction from raw data leads to the loss of statistical information of the raw data. (2) Multisegment averaging results ignore the sampling uncertainty inherent in multisegment signals and ignore the fluctuations of the signal over time. (3) Subjective cognitive uncertainties and objective experimental uncertainties occur in experimental data, which limits fault analysis by traditional probability models. Probability box (p-box) theory is an effective solution for the above problems. The p-box is based on the upper and lower boundaries of probability and can unify the various deterministic and uncertain knowledge of real numbers, intervals, random distributions and Dempster-Shafer structures into the p-box model [18-21]. The p-box model incorporates stochastic algorithms of Bayesian inference and evidence theory and the artificial intelligence algorithms of neural networks, expert systems and 5

rough set theory [22-23]. From the theoretical basis and development of the p-box, it can be viewed as a branch of information fusion technology [24]. In 2014, Du et al. [25-27] presented the technical reports of implementing p-box in fault diagnosis. They reported that the p-box was effective and practical for use in fault diagnosis. Later, the same author and Tang et al. [28,29] investigated the use of p-box for bearing fault diagnosis, but a complete and independent modeling method was lacking in these technical reports. Three p-box modeling methods for compound faults of rolling element bearings are proposed in this paper. Subjective cognitive uncertainties and objective experimental uncertainties are taken into account by the p-box models using the envelope of the data, which avoids the drawback of feature extraction discarding the rich probability statistics in raw data and allows the uncertainty of the system to be fully reflected. The sampling uncertainty and the fluctuations in the signal over time can be calculated by converting the raw data into p-boxes and then using the p-boxes as a space-time registration framework. In this way, the p-boxes can obtain the range and random distribution of uncertain data, and the statistical information of a specified amount of data can be further applied. First, to analyze the type of random distribution of the raw data, a distribution-type p-box modeling method (DTPMM) is proposed. Then, dimensionless values of the raw data are extracted; to analyze the type of random distribution of the dimensionless values, a second, dimensionless p-box modeling method (DPMM) is proposed. Finally, a raw-data p-box modeling method (RDPMM) is proposed if verification of a random distribution of the data is not possible. The DPMM and RDPMM improve the applicability of the modeling methods and can be generalized to the modeling of other faults modeling. Analyzing the information of the p-boxes through the aggregated uncertainty measurement 6

method in different ways, the features of the different p-boxes are obtained. The correct recognition rate is used as the target, and then the different p-box modeling methods and traditional methods are compared to verify the validity of the modeling method.

2 P-box modeling of fault signals 2.1 Basic notions Let  F , F  denote the set of all nondecreasing functions F from the reals into [0,1] such that F  x   F  x   F  x  [30-31]. When the functions F

and F

circumscribe an imprecisely known probability distribution, we call  F , F  , specified by the pair of functions, the “probability box” or “p-box” for that distribution. This means that if  F , F  is a p-box for a random variable X whose distribution F is unknown except that it is within the p-box, then F  x  is a lower bound on F  x  , which is the (imprecisely known) probability that the random variable X is smaller than x [32]. Likewise, F  x  is an upper bound on the same probability. From a lower probability measure P for a random variable X , one can compute the upper and lower bounds on the distribution functions using

FX  x   1  P( X  x) ,

(1)

FX  x   P( X  x) .

(2)

Based on the comparison of the probabilistic characteristics of bearing signals, three kinds of p-box modeling methods are proposed by the author, and their relationships are shown in Fig. 1.

7

Bearing signal

Does it satisfy a random distribution?

Yes

DTPMM method

No

RDPMM method

No

Do the dimensionless values satisfy a random distribution?

Yes DPMM method

Fig. 1. Relationships among modeling methods.

The DTPMM is selected when the collected bearing signals satisfy a random distribution. Otherwise, the DTPMM is inapplicable, and then the DPMM is selected. However, if the features from the bearing signals do not satisfy a random distribution either, the DPMM is inapplicable, and then the RDPMM is selected. 2.2 The DTPMM In most cases, bearing signals collected from a sensor follow a random distribution because of measurement errors of the sensor, the different positions of the measurement and the variability of working conditions. A vibration acceleration signal is considered for the DTPMM. The signal is considered as the random variable X , which satisfies commonly used random distributions such as the Weibull

distribution, the lognormal distribution, the normal distribution, the Type I extreme value (Gumbel) distribution, the Type II extreme value (Frechet) distribution and the uniform exponential distribution. The random distribution parameter (1 ,..., n ) is an n-dimensional vector. Accordingly, a finite Dempster-Shafer structure on the real line  can be identified with its basic probability assignment [33], which corresponds to the mapping 8

m : 2  [0,1]

(3)

where m()  0 , m(ai )  pi for focal elements ai  ， i  1,2,..., n , and m( D)  0 whenever D  ai ， i  1,2,..., n , pi  0 and

p

i

1.

To illustrate the DTPMM modeling method, consider a normal distribution for random variable X . For simplicity, we also assume that the focal elements are closed intervals rather than more complicated sets. Hence, the  1 can be defined in left right an interval mean  [ ,  ] , and the  2 can be defined in an interval deviation

 [ lower , upper ] , where  left and  right refer to the left bound and right bound on the respective probabilities of the p-box, as do  lower and  upper . We define a Dempster-Shafer structure as a collection of pairs consisting of an interval of

([ left ,  right ],[ lower , upper ]) and a mass

m . The lower part of the left side of the

left p-box is determined by the CDF of the normal distribution with mean  and lower standard deviation  , which can be expressed as

  

left 1



 

 ,  1lower  , m1 ,  2left ,  2lower  , m2 ,...,  nleft2 ,  nlower 2  , mn 2

 .

The upper part of the left side is determined by the distribution with mean  left and upper standard deviation  by the following expression:

    left upper     left upper  left upper    n 1 ,  n 1  , mn 1  ,    n  2 ,  n  2  , mn  2  ,...,  n ,  n  , mn 2  2   2 2  2     2





  .  

Similarly, the lower part of the right side is determined by the distribution with mean

 right and standard deviation  lower , which can be expressed as

  

right 1

 



lower  ,  1lower  , m1 ,  2right ,  2lower  , m2 ,...,  nright 2 ,  n 2  , mn 2

9

 .

The upper part of the right side is determined by the distribution with mean  right and upper standard deviation  , namely,

    right upper     right upper  right upper    n 1 ,  n 1  , mn 1  ,    n  2 ,  n  2  , mn  2  ,...,  n ,  n  , mn 2  2   2 2  2    2 





 .  

The plausibility function Pl : 2  [0,1] corresponding to a Dempster-Shafer structure with a basic probability assignment is the sum of all masses associated with sets that overlap with or merely touch the set b   . Thus, FX  x  =Pl(b) 



m( a ) 

a a  b 



m(ai ) .

(4)

i ai b 

The belief function Bel : 2  [0,1] is the sum of all masses associated with sets that are subsets of b   so that

FX  x  =Bel(b)   m(a )  a a b

 m( a ) . i

(5)

i ai  b

Clearly, Bel(b)  Pl(b) . In fact, a Dempster-Shafer structure could also be identified with either of these functions. Given only one of m , Pl or Bel , one can compute the other two. Implementation of such a Dempster-Shafer structure on a computer would thus require storage for 3n floating-point numbers, one for each pi and two for each corresponding interval. However, since the system response must be treated as a black box, each parameter must be discretized in any case. Therefore, we propose to discretize each parameter in such a way that its upper and lower probabilities, Pli and Beli , are the plausibility and belief functions, respectively, of a random set (i , mi ) . Suppose N observations are made of ([ left ,  right ],[ lower , upper ]) , each of which results in an imprecise (nonspecific) measurement given by a set of values A . Two discretization 10

methods are proposed: the averaging discretization method (adm) and the outer discretization method (odm) [29], as shown in Fig. 2.

(b) (a) Fig. 2. Discretization of the p-box. (a) Average discretization method (adm), (b) outer discretization method (odm).

Consider Fig. 2(a) and (b). In the first step of both methods, the [0,1] ordinate intervals of Pli and Beli are both discretized into n subintervals of length

m j  0 ( j  1,..., n) . By definition, let m0  0 , and let Pli1 and Beli1 indicate the inverse functions of Pli and Beli , respectively. According to the adm, the jth focal element of (i , mi ) is the interval j 1  1  j 1 mj  m  1  Ai , j   Pli   ms   , Beli   ms  j   ，j  1,..., n 2  2   s 0  s 0  and its basic probability assignment is

mi ( Ai , j )  m j .

(6)

(7)

According to the odm, the jth focal element of (i , mi ) is the interval   j 1   j 1   Ai , j   Pli1   ms  , Beli1   ms   ，j  1,..., n  s 0   s 0   

(8)

Pli1 (0)  lim{ ] Pli1 ( x) ， Beli1 (1)  lim[ } Beli1 ( x)

(9)

where

x 0

x1

11

and its basic probability assignment satisfies Eq. (7). The adm produces a random set whose plausibility ‘averages’ are Pli and whose belief ‘averages’ are Beli . Therefore, the adm is very effective for calculating the expectation of a parameter or the expectation of its image through Eq. (6). On the other hand, the odm ensures that F  Pli and that F  Beli . The odm may severely overestimate Pli and underestimate Beli . 2.3 The DPMM As mentioned previously, bearing signals seen as a random variable X must satisfy a random distribution; otherwise, the DTPMM is inapplicable. Although applying the extracted features of bearing signals will lead to the loss of rich statistical information, the features can more intuitively reflect the nature of bearing signals, which consists of the compensation and support of the raw signal. In order to solve the above problems, the DPMM is proposed. Its specific algorithm is described by the following procedure. Bearing signals Fs  [1 , 2 ,

, mn ] can be transformed

into the matrix  11 12   21  22  Fs    ij   i1    m1  m 2

1n   2 n 

      mn 

(10)

where Fs denotes the bearing signals, ij (i  1,

, m),( j  1,

, n) denotes a value

of vibration acceleration, and the subscripts m and n denote sampling times and sampling frequency, respectively. The skewness features or the kurtosis features from each row vector of Fs are 12

then calculated. The skewness feature can be expressed as follows: 3 1 n  ij   i     n j 1 Si = , i  1, 1 n 3  ij n j 1

where Si (i  1,

,m

(11)

, m) is the skewness feature of the ith row vector of Fs , and  i is

the mean of the ith row vector. The equation for the kurtosis feature is 4 1 n       ij i  n  Ki = j 1 n , i  1, 1 4 ij n j 1

where Ki (i  1,

,m

(12)

, m) is the kurtosis feature of the ith row vector of Fs . A vector

comprising skewness features is obtained from all bearing signals as follows:

S vector =[S1 , S2 , ..., Sm ]T

(13)

The K vector can be expressed in terms of the kurtosis features as follows:

K vector =[K1 , K2 , ..., Km ]T

(14)

The feature vector with respect to the skewness features or kurtosis features can be seen as a random variable, which may satisfy a random distribution. Thus, the Dempster-Shafer structures for the DPMM are obtained and can be discretized to form the p-box. The DPMM is based on the DTPMM, where the random distributions related to the DPMM are limited to the extracted features. 2.4 The RDPMM Deterministic random distribution are essential for using the DTPMM and DPMM because they allows us to estimate (1 ,..., n ) and then discretize Dempster-Shafer structures based on the estimated (1 ,..., n ) to form the p-box.

13

However, in the majority of cases, bearing signals or their features do not satisfy a random distribution. To address this problem, the RDPMM based on the definition of the p-box is proposed when random distributions cannot be satisfied or verified. In the first step of the RDPMM, two vectors are made from the maxima and minima from each column vector of Fs , which can be expressed as column 1 i Fsmin  [ min ,...,  min ,

n i i 1 , ,  min ]T ,  min   min

(15)

column 1 i n i i 1 . Fsmax  [ max ,...,  max ,...,  max ]T ,  max   max

(16)

column column where Fsmin and Fsmax denote the vectors made by the minima and the

i i maxima, respectively, and  min (i  1,..., n) and  max (i  1,..., n) are the minima and

maxima from the ith column vector of Fs , respectively. A Dempster-Shafer structure is obtained for the RDPMM with the following expression:

([

1 min

1 ,  max ], m1 ),

i i n n , ([ min ,  max ], mi ),..., ([ min ,  max ], mn )

where mi (i  1,..., n) satisfies

mi ( Ai , j )  m j , j  1,..., n

(17)

Subintervals of equal length 1 n are used because the probability is 1 n for each value of the vibration acceleration, where n is the sampling frequency. Thus, Ai , j can be expressed through the adm or odm as follows:

 1 1  1 1   Ai , j  Pli 1   j    , Beli 1   j     ，j  1,..., n , 2  2   n n  or

14

(18)

 1  1  Ai , j   Pli 1   j  1  , Bel i 1   j    ，j  1,..., n . n  n  

(19)

Then, Eq. (17) can be simplified as follows: mi ( Ai , j ) 

1 . n

(20)

A discrete Dempster-Shafer structure is produced, and discrete samples on the upper and lower bounds are measured by the adm or odm, respectively, as follows:

1  ([ x1 ,  x1 ], ), n  where

1 1  , ([ x ,  xi ], ),..., ([ x ,  xn ], )  i n n n 

 x   x , and  x   x i

i

i

i 1

whenever  x   x . Thus, i i 1

approximated by discrete sampling on upper bound, and F

F

can be

can be made

approximated by discrete sampling on lower bound, which will therefore be an approximation to the p-box.

3 Feature vector of the p-box Analyzing the information of the p-box through the aggregated uncertainty measurement method in different ways, a feature of the p-box can obtained, which is used for identification. The feature of the p-box can be a single scalar value or interval. The aggregated uncertainty measurement methods for the p-box can be expressed as follows: Method one: Obtain the aggregated width to use basic probability assignment for all focal element intervals of weight, which can be expressed as n

e1   (mi  xi  xi )

(22)

i 1

where e1 is the first feature of the p-box, and xi and xi are the bound values of the focal element interval. The value of the aggregated width, which represents the 15

uncertainty area of the p-box, reflects the degree of distinguishing in a pattern recognition system, with typical p-box features representing the different classes to be distinguished, and a smaller value means a smaller uncertainty area in the p-box [34]. In cases where this is true, a smaller value indicates that the p-boxes are easier to distinguish in a pattern recognition system; otherwise it is difficult to distinguish p-boxes. Method two: Obtain the logarithm of the aggregated width to use basic probability assignment for all focal element intervals of weight, which can be expressed as n

e2   (mi  log 2 xi  xi )

(23)

i 1

where e2 is the second feature of the p-box. Method three: Obtain the logarithm of the aggregated width with 1 added to the argument of the logarithm to use basic probability assignment for all focal element intervals of weight, which can be expressed as n

e3   (mi  log 2 (1  xi  xi ))

(24)

i 1

where e3 is the third feature of the p-box. Method four: Obtain the lower and upper bounds of the p-box to use basic probability assignment for the aggregated interval bound values of weight, which can be expressed as n  n  e4    mi xi ,  mi xi  i 1  i 1 

(25)

where e4 is the fourth feature of the p-box, which is an interval value. Method five: Obtain aggregated bound values from the lower and upper bounds of the p-box under the conditional value of the aggregated assignment  , where  16

satisfies n 1

m i 1

i

n

    mi

(26a)

i 1

The aggregated bound value from the lower and upper bounds of the p-box can be expressed as n  n  e5     xi ,   xi  i 1  i 1 

(26b)

where e5 is the fifth feature of the p-box, which is an interval value. Method six: Obtain the contradiction interval statistics of the p-box; the aggregated uncertainty measurement results can be expressed as n

e6   mi  log 2 (c2'  c1' )

(27)

i 1

' ' where e6 is the sixth feature of the p-box, and c1 and c2 are the mean of the

probability statistics from the lower and upper bounds of the p-box, respectively. Thus, for a vector of the p-box, Eqs. (22) to (27) provide a way to calculate the features of the p-box. Accordingly, the feature vector of the p-box can be expressed in terms of the features as follows: e=  e1 , e2 , e3 , e4 , e5 , e6 

T

(28)

4 Experiment and discussion 4.1 Collecting the bearing signals In the present approach, the influence of various loads and speeds on the p-box leads to a shift in the lower and upper bounds of the p-box to negative or positive values on the abscissa, but the aggregated width of the p-box is unaffected because the 17

signal patterns remain unchanged. Thus, the diagnostic results are unaffected as long as the aggregated width remains unchanged. For simplicity, we can ignore the influence of different loads and speeds on the p-box to study the classification performance of the different p-box modeling methods. The vibration acceleration signals of rolling element bearings are used as the research object for these experiments. A kinematic diagram of the rolling element bearing fault simulation test bench is shown Fig. 3.

Fig. 3. A kinematic diagram of the bench.

The experimental setup, shown in Fig. 4, includes a motor, coupling, sensor, bearing seat, PXI-4472 data acquisition card and laptop. The signal acquisition system is an NI PXI-1042Q high performance acoustic vibration testing system. The sensor is a PCB M603C01 ICP accelerometer. The acceleration transducer is installed on the bearing seat by using a magnetic base to collect bearing signals. Bearing signals were collected through two channels of a dynamic signal test and analysis system.

18

Fig. 4. Experimental setup.

The experimental object is a 30305 SKF tapered roller bearing. The basic geometric parameters and fault characteristic frequencies are listed in Table 1. Table 1 Geometric parameters and fault characteristic frequencies of the 30305 SKF tapered roller bearing Outer Inner Pitch Roller Contact Number of diameter, diameter, diameter, diameter, angle, º balls mm mm mm mm 30305 62 25 44.6 28 13 9.06 Inner race frequency, Hz Outer race frequency, Hz Rotational frequency, Hz Rotation rate, r/min Bearing model

100.65

70.03

13.33

800

Photographs of bearing faults are shown in Fig. 5. The grooves, which are 0.5 mm deep and 0.5 mm wide, were processed in the inner race, the outer race and the rolling elements to simulate bearing failure by using wire cut electrical discharge machining technology [35,36].

(a) (b) (c) Fig. 5. Photographs of bearing faults. (a) Inner race fault, (b) Outer race fault, (c) Rolling element fault.

The main methodology adopted for the experiment can be categorized into five parts as follows: 19

(1) Connecting the experimental hardware as seen in Fig. 3, including assembly of the experimental rolling element bearing and the ICP accelerometer as well as connection of the dynamic signal test and analysis system. (2) Calibrating the ICP accelerometer with a YE5501 sensitivity calibrator (i.e., a handheld sensor). (3) Designing a collection and storage system for the vibration acceleration signal with LabVIEW software. (4) Starting up the experimental device and adjusting the motor speed to 800 r/min with a frequency converter until speed stabilization. (5) Setting the sampling frequency to 10240 Hz [36] according to Balderston [37], who noted that at high frequencies in the tens of kHz, measurable acceleration levels correspond to extremely small displacements, which could be accommodated in the clearance space between the surface asperities of a bearing race in its housing, even after fitting, and thus natural frequencies would not be greatly modified by the mounting. The sample time is then set to 60 seconds. The dynamic signal test and analysis system is subsequently started for the collection and storage of vibration acceleration signals. 4.2 The time domain waveforms and power spectra of the experimental data The time domain waveforms and power spectra of the partial bearing signals under eight conditions are shown in Fig. 6.

20

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

21

Fig. 6. The time domain waveforms and the power spectra of the experimental bearing acceleration signals. (a) Healthy bearing signal, (b) Inner race fault signal, (c) Outer race fault signal, (d) Rolling element fault signal, (e) Inner and outer race fault signals, (f) Inner race and rolling element fault signals, (g) Outer race and rolling element fault signals, (h) Inner race, outer race and rolling element fault signals.

Fig. 6(a) shows the time domain waveform and power spectrum of the healthy bearing. Fig. 6(b), (c) and (d) correspond to single fault signals of the rolling element bearing, and Fig. 6(e), (f), (g) and (h) correspond to compound fault signals. The basic time domain waveforms and power spectra are not generally informative (Fig. 6), partly because they are in the presence of strong masking signals from other machine components but also because the angle of the load from the radial plane varies with the position of each rolling element in the bearing as the ratio of the local radial to the axial load changes. Thus, each rolling element has a different effective rolling diameter and is trying to roll at a different speed, but the cage limits the deviation of the rolling elements from their mean position, thus causing some random slip. The resulting change in bearing frequencies is typically on the order of 1-2%, both as a deviation from the calculated value and as a random variation around the mean frequency. This random slip, while small, does result in a fundamental change in the character of the signal, which is why the effects of some preprocessing techniques, such as high-pass or bandpass filtration, envelope detection (demodulation) and wavelet transform of the vibration signals, prior to feature extraction are generally studied to remove the uncertainty (unknown transmission path) of bearing signals in traditional rolling element bearing diagnostics [38]. However, the uncertainty of bearing signals is the main characterization of bearing signals because it always exists and changes from the generation of healthy signals to the formation of fault signals. Collecting the uncertainty of bearing signals provides a new method for rolling 22

element bearing diagnostics. As expected, the uncertainty of the bearing signals changed over time and can be enveloped by a p-box, which avoids the decoupling process for compound fault signals [26]. 4.3 Analysis of experimental data The bearing signals were analyzed using the MATLAB statistics toolbox. The collected signals from the healthy bearing were verified to have a normal

Standard line of normal distribution Data of healthy bearing

0.999 0.997 0.99 0.98 0.95 0.90

Probability

Probability

distribution, as shown in Fig. 7(a).

0.75 0.50 0.25 0.10 0.05 0.02 0.01 0.003 0.001

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

0.999 0.997 0.99 0.98 0.95 0.90 0.75 0.50 0.25 0.10 0.05 0.02 0.01 0.003 0.001

-30

Data

Standard line of normal distribution Data of inner ring fault

-20

-10

0

10

20

Data

(a) (b) Fig. 7. Results of the verification of normal distributions of the experimental data. (a) Healthy bearing data, (b) Inner race fault data.

The data from the healthy bearing basically coincide with the standard line of normal

distribution

(red

dotted

line),

as

shown

in

Fig.

7(a).

The

Kolmogorov-Smirnov test of normality was used to accurately illustrate the normal distribution of the data from the healthy bearing; the Kolmogorov-Smirnov test finds the degree of plausibility of obtaining the points of a sample distribution from a reference normal distribution [39]. The result show that the agreement between the data from the healthy bearing and the standard normal distribution is 91.23%, which means the healthy bearing data satisfy a normal distribution. For the 23

sake of comparison, the experimental results for the data from the inner race fault is shown in Fig. 7(b), which does not satisfy normal distribution since the data do not coincide with the standard line of normal distribution. Accordingly, as indicated in Section 2.2, the p-box of the healthy bearing can be obtained by the DTPMM because the data from the healthy bearing satisfy a normal distribution (Fig. 7(a)), as shown in Fig. 8. 1.0

N(μright, σupper)

Cumulative probability

N(μleft, σupper) 0.8

Pl-healthy Bel-healthy

0.6

0.4

0.2

N(μright, σlower)

N(μleft, σlower)

0.0 -8

-6

-4

-2

0

2

4

6

8

Interval

Fig. 8. P-box of the healthy bearing signals using the DTPMM.

The uncertainty of the data is included in the p-box discretized by the adm (Fig. 8), where the number of focal elements is 1000. The p-box of the healthy bearing differs from the traditional Monte Carlo method in forming a single CDF (Fig. 8), which includes F and F of the p-box calculated with Eqs. (4) and (5), where F and F are a stepwise approximation to Pli and Beli , respectively. The upper part of the left side of the p-box is determined by the CDF of the normal distribution with mean left and standard deviation  upper . The lower part of the left side is determined by the distribution with mean left and standard deviation  lower . The right side of the p-box similarly involves the right value for the mean but both values for the standard deviation. The remaining bearing signals (i.e., the 24

bearing signals of the remaining seven conditions) were assessed to determine whether they satisfied random distributions, was assessed, but the experimental results showed that they did not. Thus, the DTPMM is inapplicable in these cases. 4.4 P-boxes modeling of experimental data As mentioned above, the bearing data did not satisfy a random distribution in the majority of cases except for the data from the healthy bearing. According to Sections 2.3 and 2.4, on the one hand, the fault features could be analyzed to determine if they satisfy a random distribution, and the type of random distribution could be obtained to establish a p-box. On the other hand, the raw data could be used to establish a p-box, which avoids needing to verify satisfaction of a random distribution. We first consider the DPMM modeling method. In the first step of the method, the kurtosis features are obtained from the data from the inner race fault calculated with Eq. (12), and its normal distribution can be assessed. The kurtosis data basically coincide with the standard line of normal distribution (red dotted line), as shown in Fig. 9. The agreement between the kurtosis data and the standard normal distribution is only 89.87% with the Kolmogorov-Smirnov test because of the sudden decrease in the amount of data caused by feature extraction; these results can be improved by increasing the amount of the data from the inner race fault. The skewness features can be obtained with Eq. (11). The normal distribution of the skewness data is similar to that of the kurtosis data.

25

0.99 0.98

Standard line of normal distribution

0.95

Data of kurtosis

Probability

0.90 0.75

0.50

0.25 0.10 0.05 0.02 0.01 -0.4

-0.2

0

0.2

0.4

Data

Fig. 9. Normal distribution verification for the kurtosis data.

Similar to the fault features calculated from the data from the inner race fault and the verification of a normal distribution, the kurtosis features of the raw data from the remaining bearing conditions can also be obtained (calculated with Eq. (12)); then, the random distribution of these features can be verified to establish a p-box based on the DPMM, where the result for the random distribution of the kurtosis features satisfies a normal distribution. Accordingly, the p-boxes of the different conditions discretized by the adm were established based on the kurtosis features as shown in Fig. 10, where H stands for healthy bearing, IR stands for inner race fault, OR stands for outer race fault, RE stands for rolling element fault, IOR stands for inner and outer race faults, IRRE stands for inner race and rolling element faults, ORRE stands for outer race and rolling element faults, and IORRE stands for inner race, outer race and rolling element faults.

26

(a) P-box of H, (b) P-box of IR, (c) P-box of OR, (d) P-box of RE, (e) P-box of IOR, (f) P-box of IRRE, (g) P-box of ORRE, (h) P-box of IORRE.

Pl Pl Bel

Bel

Pl

Bel

Be l

Pl

(h)

(g)

Pl

1.0

0.8

0.6

Pl Pl

(e)

Bel

Bel

P-box es

Cumulative probability

(f)

Bel

(d)

Pl

(c)

Bel

0.4

(b)

0.2

(a)

0.0 0.5

1.0

1.5

2.0

2.5

Kurtosis Fig. 10. Kurtosis p-box modeling results.

‘Pl’ and ‘Bel’ refer to the plausibility upper and belief lower bounds of the p-boxes, respectively (Fig. 10), calculated with Eqs. (4) and (5). The number of focal elements is 1000. Fig. 10(a), (b), (c) and (d) depict p-boxes of the single fault signals (i.e., the p-boxes of the healthy bearing, inner race fault, outer race fault and rolling element fault), where the aggregated width for the p-box of the healthy bearing can be seen to be smaller than that of the other p-boxes, which means that the p-box of the healthy bearing can be easily distinguished in a pattern recognition system, because the healthy bearing signals envelop a minimal amount of uncertainty. Fig. 10(e), (f), (g) and (h) show the p-boxes of compound fault signals, where the aggregated width for the p-box of the inner and outer race faults is similar to the aggregated width for the p-box of the inner race and rolling element faults (see Figs. 10(e) and (f)). It is worth emphasizing that the aggregated width for the p-box of the inner race and rolling element faults (Fig. 10(f)) is smaller than the p-box of the rolling element fault (Fig. 10(d)) since there are signals between the inner fault data and rolling element fault data that cancel each other out in the p-box of the inner 27

and rolling element faults (i.e., the p-box of the compound fault signal). P-boxes for the different conditions discretized by the adm were also established based on the skewness features (calculated with Eq. (11)), as shown in Fig. 11, where the abscissa indicates the skewness value and the ordinate indicates the cumulative

Cumulative probability

probability density values. The number of focal elements is 1000. 0.99 0.66 0.33 0.00 0.99 0.66 0.33 0.00 0.99 0.66 0.33 0.00 0.99 0.66 0.33 0.00 0.99 0.66 0.33 0.00 0.99 0.66 0.33 0.00 0.99 0.66 0.33 0.00 0.99 0.66 0.33 0.00

P-box of H

(a)

Pl

P-box of IR

Bel

(b)

Pl

Bel

P-box of OR

(c)

Pl

Bel

P-box of RE

(d) Pl

Bel

P-box of IOR

(e)

Pl Bel

P-box of IRRE

(f)

Pl Bel

P-box of ORRE

(g)

Pl Bel

P-box of IORRE

(h)

Pl Bel

-2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

2.5

Skewness

Fig. 11. Skewness p-box modeling results. (a) P-box of H, (b) P-box of IR, (c) P-box of OR, (d) P-box of RE, (e) P-box of IOR, (f) P-box of IRRE, (g) P-box of ORRE, (h) P-box of IORRE.

The aggregated width of the p-box of the healthy bearing can be seen to be smaller than that of the other p-boxes (see Fig. 11), but it is larger than the aggregated width for the kurtosis p-box of the healthy bearing (Fig. 10(a)). Similarly, the aggregated widths of the skewness p-boxes for the other models are larger than the aggregated widths of the respective kurtosis p-boxes. For the sake of comparison, the aggregated widths of the kurtosis p-boxes (see Fig. 10) and the skewness p-boxes (see Fig. 11) were calculated with Eq. (22) and are shown in Fig. 12.

28

0.05

0.3581

0.3187

0.2368 0.1611

0.1707

0.1131

0.1742

0.15 0.10

0.2101

0.20

0.2312 0.2731

0.25

0.251 0.2921

0.30

0.0413

Aggregated width

0.35

0.3796

0.4007

Kurtosis p-boxes Skewness p-boxes

0.40

0.4012

0.45

0.00 (a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure legends

Fig. 12. Kurtosis p-boxes vs. skewness p-boxes in terms of aggregated width.

The aggregated widths of the kurtosis p-boxes are smaller than those of the skewness p-boxes (Fig. 12), which means that the kurtosis p-boxes are more conducive to a pattern recognition system. This is because the kurtosis is a greater-order moment of a probability function than the skewness, so the kurtosis p-box has more certainty, and therefore its aggregated width is smaller. As mentioned previously, raw data are used to establish p-boxes based on the RDPMM. Using Eqs. (18) and (20) to discretize the Dempster-Shafer structure for the RDPMM, eight p-boxes were generated, as shown in Fig. 13 (calculated with Eqs. (1) and (2)), where the abscissa indicates the interval value of a random distribution and the ordinate indicates the cumulative probability density value. The number of focal elements is 10240. To better present each p-box, we cropped each p-box data to the interval values between -2 and 2.

29

1.0

1.0

(b)

0.8

0.6

Pl-H Bel-H Pl-IR Bel-IR Pl-OR Bel-OR Pl-RE Bel-RE

0.4

0.2

0.0 -2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

Cumulative probability

Cumulative probability

(a) 0.8

0.6 Pl-IOR Bel-IOR

0.4

Pl-IRRE Bel-IRRE Pl-ORRE Bel-ORRE

0.2

Pl-IORRE Bel-IORRE 2.0

0.0 -2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

Interval

Interval

Fig. 13. RDPMM p-box modeling results. (a) P-boxes of single faults and (b) P-boxes of compound faults.

The aggregated widths of the p-boxes based on the RDPMM (Fig. 13) are smaller than obtained from the DPMM because the feature extraction of raw data was not preformed prior to performing the RDPMM, so the integrity of raw data has been maintained, increasing the certainty of the p-boxes. However, the higher integrity of the raw data leads to an increase in the amount of processed data, resulting in greater computational costs and time. 4.5 Feature extraction from different p-boxes It is worth emphasizing that the kurtosis p-boxes are more conducive to a pattern recognition system than the skewness p-boxes based on the analysis from the previous sections. For this calculation, a conditional value for the aggregated assignment  of 0.95 was used in Eq. (26) [34]. Only the features from the kurtosis p-boxes were extracted for the DPMM and were calculated with Eqs. (22) to (27); these features are listed in Table 2.

30

Table 2 Features of the kurtosis p-boxes (fs=10240 Hz) Bearing status

Method 1

Method 2

Method 3

Method 4

Method 5

Method 6

H

0.0413

-4.6258

0.0584

[0.0085, 0.0498]

[0.0038, 0.0562]

0.8950

IR

0.1131

-2.2565

0.2795

[0.0831, 0.2982]

[0.5124, 0.7920]

1.4719

OR

0.2510

-2.1018

0.3184

[0.0129, 0.2381]

[0.3331, 0.7186]

1.1543

RE

0.2312

-2.1494

0.2985

[0.1443, 0.3755]

[0.5513, 0.8506]

1.3022

IOR

0.1707

-2.6162

0.2257

[0.0056, 0.1763]

[0.2522, 0.4934]

1.1208

IRRE

0.1611

-0.8248

0.6587

[0.2384, 0.3496]

[0.2062, 1.0291]

0.3429

ORRE

0.3187

-1.7621

0.3921

[0.0058, 0.3129]

[0.3972, 0.8905]

1.0611

IORRE

0.2368

-2.2278

0.3008

[0.0908, 0.1460]

[0.1166, 0.5028]

0.9832

Vectors of the different kurtosis p-box features were obtained (Table 2). Methods 1, 2, 3 and 6 show the single scalar features of different p-boxes, where method 1 results in the aggregated widths of the different p-boxes. Methods 4 and 5 show the interval features of the different p-boxes. Similarly, the features of the p-boxes obtained from the RDPMM were extracted as listed in Table 3. Table 3 Features of p-boxes obtained from the RDPMM (fs=10240 Hz) Bearing status

Method 1

Method 2

Method 3

Method 4

Method 5

Method 6

H

0.0060

-7.4355

0.0086

[0.2116, 0.2176]

[0.1999, 0.2078]

1.2308

IR

0.0430

-4.5644

0.0607

[0.3342, 0.3772]

[0.2783, 0.3320]

0.6282

OR

0.1435

-2.8466

0.1926

[1.5208, 1.6643]

[1.3060, 1.4989]

0.7964

RE

0.0824

-3.7121

0.1135

[1.0571, 1.1395]

[0.8305, 0.9569]

1.7471

IOR

0.0974

-3.3884

0.1337

[1.3548, 1.4522]

[1.0771, 1.1971]

1.8112

IRRE

0.0550

-4.2693

0.0769

[0.9064, 0.9614]

[0.7235, 0.8034]

2.0429

ORRE

0.1647

-2.6460

0.2190

[1.0759, 1.2407]

[0.8247, 1.0450]

0.8285

IORRE

0.1103

-3.2212

0.1505

[1.7569, 1.8673]

[1.4807, 1.6248]

1.6037

Contrasting the values of method 1 obtained from the kurtosis p-boxes and the RDPMM (Tables 2 and 3, respectively), the aggregated widths of the p-boxes obtained from the RDPMM are smaller than those of the kurtosis p-boxes, which means that the p-boxes obtained from the RDPMM are more conducive to a pattern recognition system than the kurtosis p-boxes. The effectiveness of the aggregated width is demonstrated as follows. 31

4.6 Pattern recognition with different methods The correct recognition rate is taken as the target. The kurtosis p-boxes and the p-boxes obtained by the RDPMM were quantitatively compared. A total of 200 p-boxes were obtained from the bearing signals from each condition. Thus, 1600 p-boxes were obtained for all eight conditions. The feature vector sets from the 1600 p-boxes

were

calculated

with

Eq. (28).

Before

training

the

Python

environment-based support vector machine (SVM) classifier, the data were divided into training, validation and testing sets. The number of feature vectors for each set can be represented as {960, 320, 320}. A sigmoid function was selected as the kernel function, and the parameters C and gamma were set to [1e3, 5e3, 1e4, 5e4, 1e5] and [1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 1e-1], respectively. The parameters C and gamma were automatically optimized. Using 10-fold cross-validation, the correct recognition rates of the kurtosis p-boxes and the p-boxes obtained by the RDPMM are listed in Table 5, where the number of the samples, the correct cases of classification, and the total correct classification rates are represented by N, CC and TCR, respectively. Table 5 Classification performance of the kurtosis and RDPMM p-boxes for various sets State

N

H IR OR RE IOR IRRE ORRE

200 200 200 200 200 200 200

IORRE

200

Method

RDPMM

CC 200 200 185 194 175 198 177

TCR

95.2%

194

Method

Kurtosis p-boxes

CC 187 176 173 159 155 166 161

TCR

83.4%

157

The classification performances of the kurtosis and RDPMM p-boxes for various sets are given in Table 5. As expected, the experimental results show that the p-boxes obtained by the RDPMM result in the greater score, because the RDPMM maintains the integrity of the raw data, increasing the certainty of the p-boxes. 32

Time-domain statistical features have been widely used in traditional feature extraction methods. Samanta and Al-Balushi [40] studied an ANN-based fault diagnosis of rolling element bearings using five time-domain statistical features—the root mean square, variance, skewness, kurtosis and normalized sixth central moment—to classify the status of the machine in terms healthy or faulty bearings. They reported that the features extracted directly from the measured vibration data were suitable for machine diagnostics. Kankar et al. [41] considered five cases (i.e., H, IR, OR, RE and IORRE) for the bearing condition, and 14 features from the bearing signals were used for the study. Then, they observed that the classification accuracy of the SVM was better than that of the ANN. To verify the performances of the p-boxes, the results obtained from the SVM for the H, IR, OR, RE and IORRE p-boxes were chosen and compared with Ref. [41]. The results are displayed in Table 6. Table 6 Classification performance of the current methods and the method from Ref. [41] Methods

RDPMM p-boxes

Kurtosis p-boxes

Ref. [41]

Correctly classified rate

97.3%

85.2%

73.9726%

The experimental results show that the p-boxes obtained by the RDPMM and the kurtosis p-boxes obtain better classification accuracy rates than Ref. [41] in Table 6. This may be because feature extraction leads to a loss of the statistical information from raw data and ignores the sampling uncertainty and the fluctuations in the signal over time in mechanical fault diagnosis.

5 Conclusions and future work Novel modeling methods for the fault signals of a rolling element bearing have been proposed based on p-box theory in this paper. It is not necessary to decouple the individual signals from compound fault signals. An SVM was used as the

33

classification algorithm to evaluate the classification performance of the three methods for eight fault data sets. In the bearing fault simulation experiments, the random distributions of the eight different bearing signals were verified, and the experimental results showed that the healthy bearing signals and the features of all the bearing signals satisfy a normal distribution, which means that the features were extracted from the bearing signals to model the p-boxes in the vast majority of cases. Consequently, the DPMM was more efficient than the DTPMM. For features that satisfied a random distribution, the kurtosis p-boxes and skewness p-boxes were obtained based on the DPMM, and their aggregated width were calculated. The results showed that the aggregated widths for the kurtosis p-boxes were smaller than those of the skewness p-boxes. Therefore, the kurtosis p-box is more conducive to a pattern recognition system. To quantitatively compare the classification performance of kurtosis p-boxes, RDPMM p-boxes and traditional vectors on various sets, an SVM was used as the classification algorithm. The experimental results showed that the total correct classification rates for the RDPMM and kurtosis p-boxes and traditional vectors were 95.2%, 83.4% and 73.9726%, respectively. Thus, the RDPMM yields the highest classification accuracy rate. The good classification accuracy rate of the p-box modeling method suggests that the investigation is pointing in the right direction. However, the present modeling methods were only used to classify different fault types of rolling element bearings. An approach based on the p-box could be used to investigate the benefits of a nonuniform discretization step aimed at increasing the accuracy of the (lower) tails. However, approaches based on the 34

p-box have inherent shortcomings, such that the aggregated with of the p-box may not be optimal for each p-box model. This may necessitate further p-box modeling methods and may be included in the machine test and condition monitoring program since a small aggregated width of the p-box is required for the pattern recognition system. This issue will be the subject of further study.

Acknowledgements The research is supported by the Creative Research Groups Foundation of China (51621004), National Natural Science Foundation of China (51365020, 51467007), the Special Fund Project of the Hunan Provincial Civil-Military Integration Industry Development ([2018]23), Changsha Bureau of Science and Technology (KQ1701030) and Hunan Provincial Natural Science Foundation of China (2017JJ2044).

References [1] R. B. Randall, J. Antoni, Rolling element bearing diagnostics-A tutorial, Mech Syst Signal Pr, 25 (2011) 485-520. [2] H. Cao, L. Niu, S. Xi, X. Chen, Mechanical model development of rolling bearing-rotor systems: A review, Mech Syst Signal Pr, 102 (2018) 37-58. [3] D. Ho, R.B. Randall, Optimisation of bearing diagnostic techniques using simulated and actual bearing fault signals, Mech Syst Signal Pr, 14 (2000) 763-788. [4] X. S. Lou, K. A. Loparo, Bearing fault diagnosis based on wavelet transform and fuzzy inference, Mech Syst Signal Pr, 18 (2004) 1077-1095. [5] B. Samanta, K. R. Al-Balushi, Artificial neural network based fault diagnostics of rolling element bearings using time-domain features, Mech Syst Signal Pr, 17 (2003) 317-328. [6] R. R. Schoen, T. G. Habetler, F. Kamran, R. G. Bartheld, Motor bearing damage detection using stator current monitoring, Ieee T Ind Appl, 31 (1995) 1274-1279. 35

[7] P. W. Tse, Y. H. Peng, R. Yam, Wavelet analysis and envelope detection for rolling element bearing fault diagnosis-Their effectiveness and flexibilities, J Vib Acoust, 123 (2001) 303-310. [8] P. Gao, L. Hou, R. Yang et al., Local defect modelling and nonlinear dynamic analysis for the inter-shaft bearing in a dual-rotor system, Appl Math Model, 68 (2019) 29-47. [9] K. Worden, R. J. Barthorpe, E. J. Cross et al., On evolutionary system identification with applications to nonlinear benchmarks, Mech Syst Signal Pr, 112 (2018) 194-232. [10] M. Cao, J. Xiao, A comprehensive dynamic model of double-row spherical roller bearing-model development and case studies on surface defects, preloads, and radial clearance, Mech Syst Signal Pr, 22 (2008) 467-489. [11] S. P. Harsha, Nonlinear dynamic analysis of a high-speed rotor supported by rolling element bearings, J Sound Vib, 290 (2006) 65-100. [12] A. Rafsanjani, S. Abbasion, A. Farshidianfar et al., Nonlinear dynamic modeling of surface defects in rolling element bearing systems, J Sound Vib, 319 (2009) 1150-1174. [13] Y. Lei, Z. He, Y. Zi et al., Fault diagnosis of rotating machinery based on multiple ANFIS combination with GAS, Mech Syst Signal Pr, 21 (2007) 2280-2294. [14] Y. X. Zhang, R. B. Randall, Rolling element bearing fault diagnosis based on the combination of genetic algorithms and fast kurtogram, Mech Syst Signal Pr, 23 (2009) 1509-1517. [15] Z. Liu, J. Kang, X. Zhao et al., Modeling of the safe region based on support vector data description for health assessment of wheelset bearings, Appl Math Model. 73 (2019) 19-39. [16] C. Sobie, C. Freitas, M. Nicolai, Simulation-driven machine learning: Bearing fault classification, Mech Syst Signal Pr, 99 (2018) 403-419. [17] J. Zheng, H. Pan, S. Yang et al., Generalized composite multiscale permutation entropy and Laplacian score based rolling bearing fault diagnosis, Mech Syst Signal Pr, 99 (2018) 229-243. [18] M. Beer, S. Ferson, Fuzzy probability inengineering analyses, in: B. Ayyub(Ed.), Proceedings of the First International Conference on Vulnerability and Risk Analysis and Management (ICVRAM 2011) and the Fifth International Symposium on Uncertainty Modeling and Analysis (ISUMA2011), 11-13 April 2011, University of Maryland, ASCE, Reston, VA, USA, 2011, pp. 53-61. [19] M. Beer, S. Ferson, V. Kreinovich, Imprecise probabilities in engineering analyses. Mech 36

Syst Signal Pr, 37 (2013) 4-29. [20] D. Dubois, H. Prade, Possibility theory, probability theory and multiple-valued logics: a clarification, Ann Math Artif Intell, 32 (2001) 35-66. [21] S. Ferson, J. G. Hajagos, Arithmetic with uncertain numbers: rigorous and (often) best possible answers, Reliab Eng Syst Saf, 85 (1-3) (2004) 135-152. [22] Savic. Neural generation of uncertainty reliability functions bounded by belief and plausibility frontiers. ECSR, (2005) 1757-1762. [23] N. Walker. A first course in fuzzy logic, Boca Raton, Florida: CRC Press, (2006). [24] D. Lee, N. H. Kim, H. S. Kim, Validation and updating in a large automotive vibro-acoustic model using a P-box in the frequency domain, Struct Multidiscip O, 54 (2016) 1485-1508. [25] J. Ding, Y. Du, Q. Wang et al., Fault Diagnosis of Power Transformer based on Probability-box Theory, MCE, (2014) 419-422. [26] A. Liu, D. Yi, W. He et al., Gear fault diagnosis correlation analysis based on probability box theory, MCE, (2014) 423-426. [27] Du Y , Ding J , Liu A . A probability box modeling method of dimensionless mechanical fault features in time domain, MCE, (2014) 427-430. [28] T. Hong, D. Yi, D. Jiaman et al., Compound fault diagnosis based on probability box theory, ICMIC, (2017) 336-341. [29] H, Tang, Y, Du, J, Ding et al., An improved probability box discretization method for mechanical fault diagnosis, CCC, (2017) 11011-11016. [31] R. R. Yager, J. Kacprzyk, M. Fedrizzi, Advances in the Dempster-Shafer Theory of Evidence, John Wiley & Sons, New York, NY, USA, 1994. [30] K. Sentz, S. Ferson, “Combination of evidence in Dempster-Shafer theory,” SAND 2002-0835, SandiaNational Laboratories, Albuquerque, NM, USA, 2002. [32] W. L. Oberkampf, C. J. Roy, Verification and validation in scientific computing. Cambridge: Cambridge university press, 2010. [33] D. Dubois, H. Prade, Random sets and fuzzy interval analysis, Fuzzy Sets Syst 42 (1991) 87-101. [34] G. J. Klir, Uncertainty and information: Foundations of generalized information theory. 37

Wiley-IEEE Press (2005). [35] X. A. Yan, M. P. Jia, Application of CSA-VMD and optimal scale morphological slice bispectrum in enhancing outer race fault detection of rolling element bearings, Mech Syst Signal Pr, 122 (2019) 56-86. [36] T. Lin, G. Chen, W. Ouyang et al., Hyper-spherical distance discrimination: A novel data description method for aero-engine rolling bearing fault detection, Mech Syst Signal Pr, 109 (2018) 330-351. [37] H. L. Balderston, The detection of incipient failure in bearings, Material Evaluation 27 (June) (1969) 121–128. [38] K. Worden, W. J. Staszewski, J. J. Hensman, Natural computing for mechanical systems research: A tutorial overview. Mech Syst Signal Pr, 25 (1) (2011) 4-111. [39] H. W. Lilliefors, On the Kolmogorov-Smirnov test for normality with mean and variance unknown, J Am Stat Assoc, 62 (318) (1967) 399-402. [40] B. Samanta, K. R. Al-Balushi, Artificial neural network based fault diagnostics of rolling element bearings using time-domain features, Mech Syst Signal Pr, 17 (2003) 317-328. [41] P. K. Kankar, S. C. Sharma, S. P. Harsha, Fault diagnosis of ball bearings using machine learning methods, Expert Syst Appl, 38 (3) (2011) 1876-1886.

38