Application of tentative variational mode decomposition in fault feature detection of rolling element bearing

Accepted Manuscript Application of tentative variational mode decomposition in fault feature detection of rolling element bearing Tingkai Gong, Xiaohui Yuan, Yanbin Yuan, Xiaohui Lei, Xu Wang PII: DOI: Reference:

S0263-2241(18)31136-9 https://doi.org/10.1016/j.measurement.2018.11.083 MEASUR 6126

To appear in:

Measurement

Received Date: Revised Date: Accepted Date:

7 June 2018 20 November 2018 26 November 2018

Please cite this article as: T. Gong, X. Yuan, Y. Yuan, X. Lei, X. Wang, Application of tentative variational mode decomposition in fault feature detection of rolling element bearing, Measurement (2018), doi: https://doi.org/ 10.1016/j.measurement.2018.11.083

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

Application of tentative variational mode decomposition in fault feature detection of rolling element bearing Tingkai Gong1, 2, Xiaohui Yuan1, 3,*, Yanbin Yuan4, Xiaohui Lei5, *, Xu Wang5 1. School of Hydropower and Information Engineering, Huazhong University of Science and Technology, Wuhan 430074, China 2. School of Aircraft Engineering, Nanchang Hangkong University, Nanchang 330063, China 3. Hubei Provincial Key Laboratory for Operation and Control of Cascaded Hydropower Station, China Three Gorges University, Yichang 443002, China 4. School of Resource and Environmental Engineering, Wuhan University of Technology, Wuhan 430070, China 5. State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China

* Corresponding author: Xiaohui Yuan ([email protected]); Xiaohui Lei ([email protected]) Abstract: In contrast to empirical mode decomposition (EMD), variational mode decomposition (VMD) has better robustness to noise and sampling. In the VMD, the mode number is a key parameter that VMD can be applied successfully. In order to enhance the reliability of the adaptive selection to this parameter, tentative variational mode decomposition (TVMD) is proposed in this study. On the other hand, when a real signal sampled in strong noise background is analyzed by TVMD, the mode mixing problem frequently happens, dynamic time warping (DTW) is thus adopted to reconstruct original signals well. The effectiveness of the proposed approach is verified by both simulation analysis and the vibration signals of bearings with an outer race, an inner race and a rolling element faults, respectively. The experimental results indicate that the proposed method can extract the bearing fault features and has more reliable ability for selecting the mode number. Compared with EMD, ensemble empirical mode decomposition (EEMD) and complete EEMD (CEEMD), the proposed method has superior performance in fault feature detection. 1

Keywords: variational mode decomposition; dynamic time warping; rolling element bearing; fault feature detection 1 Introduction Rolling element bearings play an important role in the transmission systems of machine equipment. In real situation, they are usually operating under certain stress or load during service, thus it easily results in the presence of bearing faults, which significantly affects the healthy condition of the machine equipment [1, 2]. In order to enhance the reliability during the useful life of a physical asset, the vibration-based condition monitoring is introduced in the bearing fault diagnosis area. However, conventional vibration based on fast frequency transform (FFT) shows its inability to process non-linear and non-stationary vibration signals produced by bearing defects. Therefore, it is necessary to exploit effective signal processing methods in engineering and theory applications [3]. Differing from common signal processing techniques, where the signal analysis is frequently conducted by convolution between original signals and a predefined basis function, Huang et.al [4] developed a signal decomposing scheme termed empirical mode decomposition (EMD), which can divide a signal into the sum of limited intrinsic mode functions (IMFs) through its dynamic features in time scale. Furthermore, these IMFs often show amplitude-modulation or frequency-modulation behavior, and can represent the intrinsic nature of the system. Consequently, EMD has strong adaptive decomposition capacity to time-varying signals [5-8]. Yet, the EMD highly relies on local extrema detection, thus resulting in the high-level sensitivity to noise and sampling frequency [9-11]. As a consequence, ensemble empirical mode decomposition [12] is exploited as an extension of EMD. It contains sifting an ensemble of additive white noise, and uses the mean as the eventual consequence. Because the finite-amplitude noise added allows the ensemble to extract all possible solutions in the sifting process, different sub-signals in raw signatures can collate into the proper IMFs and the barriers of EMD are handled effectively. Many studies [13-16] employ the new approach to fault diagnosis of 2

rotating machinery. However, the additive white noise produces a number of demerits, for example, the smear of the residue noise to signal reconstruction, the low-efficiency computation, the inconvenience in selecting the amplitude of the added noise and the number of the ensemble trials. A few attempts and improvements are proposed to deal with these problems [17-19]. Recently, Dragomiretskiy and Zosso [20] presented variational mode decomposition (VMD), which is entirely non-recursive decomposition method, and can synchronously estimate these decomposed modes with sparseness property in frequency-domain. The sparseness means that the obtained modes are very compact around their individual central frequencies, but the frequency components at the both sides of these central frequencies are sparse. It is a good behavior to divide a practical complex signal because in most cases the extracted feature frequencies are far away from the central frequencies. According to this advantage, many studies use VMD to detect mechanical faults [21- 24]. However, the main limitation of VMD is that the mode number κ requires a proper definition because it directly affects the performance of VMD. Generally, there are three ways to determine this parameter in real applications. The first one uses empirical knowledge [21-23, 25-26]. Clearly, it lacks the adaptiveness in complicated situations. The second way adopts optimization algorithms, such as particle swarm and fish swarm optimizations [24, 27]. Although they can obtain an appropriate one, it is inefficiency because it needs implementing considerable trials. The last way is adaptive [28-30], namely one or more indicators in time- or frequency-domain are used to evaluate the change of the quantity of the feature information in the decomposed modes, when VMD is performed by using the different two mode numbers respectively. If the change is positive, this parameter will be updated and the same process is conducted again. Otherwise, the update to κ is terminated immediately and an optimal κ is correspondingly achieved. In fact, it is based on an assumption, that is, the change of the feature information (i.e. the value of the calculated indicator) is strictly monotonic following the update of κ. However, due to the signal randomness, the monotonicity is not always satisfied in the 3

adaptive definition process. It can be proved in the two studies [31, 32], where the used indicators such as kurtosis and signal-to-noise ratio show the non-monotonicity in optimizing their individual parameters. On the other hand, in order to avoid the over-decomposition of VMD, a threshold value defined experimentally is frequently introduced to compare with the difference between the indicators obtained by implementing VMD with two different mode numbers. If the comparison is positive, the update of the mode number is conducted again. Otherwise, it will be finished. Obviously, the experimental definition way indicates that it is hard to choose the threshold in practical applications. In this case, if using the one-time comparison is to determine whether the update is carried out, it is inadequate still because of the signal randomness. Clearly, these problems will result in the low reliability in defining the mode number adaptively. In order to improve this situation, tentative variational mode decomposition is proposed. Additionally, dynamic time warping is applied to handle the mixing problem between the obtained modes. As a result, the proposed approach can isolate a multi-component signal in strong noise background more effectively. The structure of this paper is organized as follows. The principle of VMD is briefly described in Section 2. Section 3 gives the proposed method, followed by a simulation analysis to evaluate its effectiveness. In Section 4, the experimental results are demonstrated by using the proposed method to extract fault features from vibration signals of the bearings with an inner race, an outer race and a rolling element faults, respectively. Eventually, conclusions are summarized in Section 5. 2 Principle of variational mode decomposition The objective of VMD is to isolate a complex signal into a series of sub-signals (i.e. mode) with specific sparseness property. The sparseness embodies on the frequency bandwidths of the modes, namely the frequency spectrum of each mode is compact around a center frequency. The problem is handled by solving the following constrained variational problem.

4

  j   min   t   (t )    mk (t )  e jwk t mk , wk  t     k K

subject to

m

k

k 1

   2 2

(1)

 f

where κ is the number of the modes; mκ and μκ strand for the κth mode and its center frequency, respectively. The constraint can be solved by using a quadratic penalty and Lagrangian multipliers. Correspondingly, the augmented Lagrangian can be written as: 2

 j   L(mk , wk ,  )     t   (t )    mk (t )  e jwk t t  k  

2

 f (t )   mk (t )   (t ), f (t )   mk (t ) k

2

2

(2)

k

where α accounts for the balance parameter of the data-fidelity. The solution of Eq. (1) is transformed into searching the saddle point of Eq. (2) by applying alternate direction method of multipliers (ADMM). In the process, mκ and ωκ are continually optimized until optimal values are obtained. The steps of VMD are expressed as follows: Step 1: Initialize m1 , 1 ,  1 , and iterations N starting from 0; Step 2: Update mκ and ωκ as Eqs. (3) and (4) respectively; mkn +1 ( w) 

f ( w)   i  k mi ( w)  1  2 ( w  wk ) 2



n +1 k

w

 ( w) 2

 w m (w) dw   m (w) dw 2

k

0



0

(3)

2

(4)

k

Step 3: Update β as Eq. (5), then n = n+1;

 n1   n   ( f  i mi )

(5)

Step 4: Justify convergence as Ineq. (6) for a given ε. If n ≤ N and the convergence is not satisfied, repeat Step 2. Otherwise, mκ and ωκ are achieved.

m

n 1 k

k

 mkn

2 2



(6)

3 Details of the proposed method 3.1 Tentative variational mode decomposition When an unsuitable κ is used, it leads to the over- or under- decomposition of 5

VMD. The under-decomposition means that the useful signal components in the obtained modes will be enhanced following the update of κ. Conversely, if they will be decreased, it is over-decomposition. In fact, such a case is also considered as a type of the over-decomposition. Namely, when VMD is conducted by using an updated κ, the interested signal components are strengthened, but the performance is not so apparent with inefficient computation, thus the situation is inadvisable in the actual applications. In order to determine the mode number adaptively, a common process is expressed as blow. Firstly, the mode number κ is initialized, and then VMD is implemented. Subsequently, the used indicators of all modes are calculated, and the optimal one Kn is achieved. The same steps are performed when κ = κ +1, and the new optimal Km is generated. Next, the difference λ between Kn and Km is calculated as: λ = Km - Kn

(7)

If λ < 0, it is over-decomposition and VMD is terminated. On the contrary, VMD is continuously implemented by using a new κ until the calculated λ is less than 0. As mentioned previously, if λ is a faint increase, it is also over-decomposition. In order to avoid this case, Ineq. (8) is considered as: λ > ζ

(8)

where the positive parameter ζ stands for the increase intensity of λ, that is, if λ is greater than ζ, it shows that the interested signal components are enhanced significantly , and VMD should be carried out again by a new κ. Otherwise, it will be ended. The above process basically describes the adaptive definition of the mode number. In general, either Ineq.(7) or Ineq. (8) is used to complete the adaptive process, but they all result in the weak reliability because of the signal randomness. Tentative VMD is thus considered. In the novel method, the update of the mode number will be determined by the tentative operations, i.e. multiple comparisons as Ineq. (8). the details of TVMD are given as: (1) If λ < 0, the update to κ is terminated and an optimal one is achieved. 6

(2) If 0 < λ < ζ, updating κ will not be ended immediately. Now, the tentative behavior is introduced. Namely, the mode number is still updated, and VMD using the new value is implemented. In the tentative process, if λ < 0, updating the mode number is terminated and a corresponding optimal κ is obtained. If the number of the tentative operation is greater than the predefined one Ω and λ is still less than ζ, the update is finished. If λ > ζ, Kn is replaced by the new optimal one. Besides, the tentative number is set to be 0, and the tentative process is conducted again. (3) If λ > ζ, the mode number is updated again until the aforementioned two cases happen. In fact, TVMD aims primarily to the situation, i.e. 0 < λ < ζ. In such a situation, the decomposition of VMD is maybe excessive-decomposition. Therefore, it will be proved by the tentative process. If λ > ζ in the tentative process, it demonstrates that that conclusion is pseudo, and the update of the mode number will be conducted until the corresponding stop conditions are satisfied. Accordingly, TVMD can select the mode number adaptively and has high reliability in processing a complex real signal. 3.2 Dynamic time warping Under ideal conditions, VMD can separate the signal components in a signal into different modes. However, it is hard to complete the objective in strong noise background. In order to reconstruct original signals more effectively, dynamic time warping is used in this study. DTW [33] is originally developed for spoke word recognition in different speaking speed. Let X be a sequence with the length m, X={x1, x2, …, xm}, and Y be a sequence with the length n, Y={y1, y2 ,…, yn}. A m-by-n path matrix is created, where the (ith, jth) element is the distance d(xi, yj) the two points xi and yj (Typically the Euclidean distance is considered, so d(xi, yj) = (xi - yj)2). Each element in the path matrix (i, j) stands for a alignment between xi and yj. Therefore, the warping path Z is a continuous set of these matrix elements defining a map between the two 7

time sequences X and Y. In the warping path Z, the q-th element is given as zq= (i, j)q. As a result, it has: Z  z1 , z2 , z3 ,

zq

zQ max(m, n)  Q  m  n 1

(9)

At the same time, the warping path has three constraints: (1) Boundary constraints: z1=(1, 1) and zQ=(m, n). It indicates that the warping path starts and terminates at diagonally opposite corner elements of the path matrix. (2) Continuity: if zq= (c, d), zq-1= (c’, d’), where c-c’≤1 and d-d’≤1. This means that the selectable steps in the warping path are limited at adjacent cells. (3) Monotonicity: if zq= (c, d), zq-1= (c’,d’), where c-c’≥0 and d-d’≥0. It allows the cells in Z to be monotonic spacing in time. In fact, many warping paths can meet the above constraints, but only one that can minimize the warping cost is given as: 

DTW ( X , Y )  min 



Q

z q 1

q

  

(10)

It can be determined by using dynamic programming for estimating the following recurrence, which gives the accumulative distance v(i, j) as the d(xi, yj) and the minimum of the accumulative distances of neighbor cells is:

 (i, j ) d ( xi , y j )  min{ (i 1, j 1), (i, j 1), v(i 1, j )}

(11)

The details on DTW, interested readers can refer in [33-35]. DTW has two advantages. First, the two time sequences can be the same or different lengths. Second, it has good robustness even if the different experimental parameters are adopted separately in the two time series. Due to the application of DTW, the similarity between the modes obtained by TVMD is evaluated effectively. However, the optimal κ often changes for different signals. Therefore, the number of principal modes Γshould be considered. Based on an overall consideration, it is given as:

     8

(12)

where



stands for the operator of round toward minus infinity.

The procedure of selecting the principal modes is given as follows. First, the optimal κ is gained by TVMD, the indicator of each mode is calculated, and the best mode Λ that has an optimal indicator is found. Based on Eq. (12), the main modes should include (Γ-1) ones except for the mode Λ. Then the similarity between Λ and the others is measured by DTW, and the (Γ-1) modes that are very similar to Λ (i.e. the accumulative distances between them and Λ are smaller) are obtained. Eventually, the reconstructed signals are achieved through averaging all principal modes. 3.3 Procedure of the proposed method Except that the mode number is determined by TVMD, several related parameters need to be specified in advance. The first one is the balancing index of the data-fidelity constraint α. Wang et.al [36] found that if α takes a larger value, it results in the loss of useful information. Besides, when VMD is employed to pick the transient signatures produced by mechanical faults, this parameter should be given a smaller one. Accordingly, α is set to be 2500 in this study. Center frequencies of all modes can be initialized by two ways, uniform distribution and zero. If the uniformly spaced distribution is utilized, VMD works essentially as an expansion of short time Fourier transform and can extract oscillatory components. Instead, when zero is applied, VMD shows a nature of wavelet-packet-like decomposition and has good behavior in impulsive feature extraction. Consequently, the latter one is adopted. Next, time-step of the dual ascent is defined as zero, which means that the Lagrangian multiplier is closed effectively. Additionally, the parameter ζ in Ineq.(8) needs to be selected. When the impulsiveness in the raw vibration signals is not evident, a smaller value should be proposed. If the transient features are weak and ε chooses a large value, it generates over-decomposition of TVMD easily. Therefore, ζ is set to 1 in this study. At the same time, the tentative number Ω is defined as 2. Finally, an indicator that is sensitive to impulses produced by

bearing

faults

should

be

defined. 9

In

the

fault

diagnosis

field,

Lempel–ZivComplexity [37], Approximate Entropy [38], Permutation Entropy [39] and kurtosis are frequently used. However, it is low efficiency in computation when the first three indicators are used to process the longer samples, thus kurtosis is adopted. According to the initializations mentioned above, the proposed method can be performed as the following Steps, and its flowchart is showed in Fig. 1. Start Initialize κ , ζ , Ω, α, center frequencies of all modes, time-step of the dual ascent Conduct VMD using κ and κ +1, calculate the kurtosis of all modes at each operation, then find the two max ones: Kn and Km Calculate λ= Km – Kn κ =κ+1

λ > 0

No

Yes Ω =0 K n= K m

Yes

λ > ζ No

Obtain the optimal κ κ = κ - 1

Ω = Ω+1

No

Ω> 2 Yes Obtain the optimal κ

Calculate the distances by using DTW Select principal modes

Obtain the rebuilt signals by averaging the selected modes End

Fig. 1. Flowchart of the proposed method.

Step1: Initialize these parameters: κ = 3, ζ= 1, Ω = 0, α = 2500, Center frequencies of all modes and time-step of the dual ascent are initialized as zero; Step2: Conduct TVMD using κ and κ+1 respectively, calculate the kurtosis values of all modes at each operation, find the two maximal ones: Kn and Km, and 10

calculate λ = Km - Kn. Step3: If λ < 0, the mode number = κ -1, TVMD is finished, and go to Step5. If 0< λ <ζ, κ = κ +1, and Ω = Ω +1, go to Step 2. If λ > ε, Ω = 0, and κ = κ +1, go to Step2. Step4: If Ω > 2; the current κ is the optimal one, and TVMD ends; Step5: Determinate the best mode Λ with max kurtosis, calculate the distance between it and the other modes, and select principal modes; Step6: Obtain the reconstructed signals through averaging the principal modes. 3.4 Simulation analysis In order to evaluate the proposed approach, a simulated signal is written as: M

x(t )  1.2cos(60 t )  1.5cos(130 t )   3e10t sin(10 t )  x3

(13)

i 1

where M stands for the repeated frequency of impulse and is set to 32. Because VMD is robust to noise, the standard deviation of Gaussians white noise x3 takes the bigger value 3. The sampling frequency is 1024 Hz, and the data length is 1024 points. Fig. 2(a) plots the time waveform of the simulation signals. (a)

Amplitude

20 10 0 -10 -20

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Time (s)

(b) 0.4

Amplitude

0.3

65 Hz 30 Hz

0.2 0.1 0

0

64

128

192

256

320

384

448

512

Frequency (Hz)

Fig. 2. The simulated signals: (a) original waveform, (b) frequency spectrum.

It can be noticed in Fig. 2(b) that the impulsive features are completely inhibited by the harmonic components and the strong noisy signatures, thus the proposed method is employed. After the improved VMD is performed, the mode number is 11

defined as 5. Figs. 3(a) and 3(b) present the time waveforms and frequency spectra of the five modes, respectively. In order to emphasize the impulsive features in these modes, then DWT is applied. Fig. 4(a) presents the selected modes, i.e. mode1 and mode3. After they are reconstructed, it is observed that the two impulsive feature frequencies 32Hz and 64Hz are evident in the envelope spectrum in Fig. 4(b). The results indicate that the proposed method can extract the impulsive features in the strong noise background. (a)

(b) A

1

32Hz

0.5 0

0

64

128

192

256

320

384

448

512

192

256

320

384

448

512

192

256

320

384

448

512

256

320

384

448

512

192 256 320 Frequency (Hz)

384

448

512

A

0.4

56Hz

0.2 0

0

64

128

A

1

64Hz

0.5 0

0

64

128

A

0.4

158Hz

0.2 0

0

64

128

192

A

0.4

220Hz

0.2 0

0

64

128

Fig. 3. Results of TVMD: (a) five modes, (b) corresponding frequency spectra. (b)

Amplitude

2 1 0 -1 -2

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Time (s) 0.15

64 Hz

32 Hz Amplitude

(a)

0.1

0.05

0

0

32

64

96

128

160

192

224

Frequency (Hz)

Fig.4. (a) principal modes, (b) the reconstructed signals and its envelope spectrum.

12

256

4 Applications In this section, practical vibration signals collected from faulty bearings are analyzed for further validating the efficiency of the proposed method. They are provided by the bearing data center of Case Western Reserve University (CWRU) [40]. As illustrated in Fig. 5, the test rig is composed mainly of a 2 HP motor, a torque transducer, a dynamometer, and various control units.

Fig. 5. The bearing tested rig.

In the experiment，the types of the used bearings include SKF 6203 and 6205, and their individual structural parameters are listed in Table 1. Single-point faults were introduced into experiment using electro-discharge machining. Raw vibration signals were sampled with the sampling frequency of 12 kHz, and the sampling time is 0.5s. Table 1. Structural parameters of the tested bearings.

Type

Pitch diameter (mm)

Outerside diameter

Inside diameter

Contact

(mm)

(mm)

angle

Number of rolling elements

SKF 6205

39

52

25

0

9

SKF 6203

28.5

40

17

0

8

4.1 Fault feature detection for an inner race fault A fault with 0.18 mm diameter was introduced on the inner race of bearing SKF 6203 installed at fan end, and the shaft rotating speed is 1730 rpm. By mathematical calculations, the ball passing frequency inner race (BPFI) is about 142.6 Hz. The vibration signals are achieved through an accelerometer attached to the motor 13

supporting base plate. Because the sensor is far from the fault and the signal attenuation is very fast in transmission path, the impulsive features are almost absent in Fig. 6(a). (a)

(b)

Fig. 6. Original signals for an inner race fault: (a) time waveform, (b) frequency spectrum. (b) 1

0

0.5

-0.05

A

0.05

0 0

0.1

0.2

0.3

0.4

0.5 0.4

0

0.2

-0.05

A

0.05

0 0

0.1

0.2

0.3

0.4

0.5 0.4

0

0.2

-0.05

A

0.05

0 0

0.1

0.2

0.3

0.4

0.5 0.4

0

0.2

-0.1

A

0.1

0 0

0.1

0.2

0.3

0.4

0.5

0.5

0.4

0

0.2

-0.5

A

Mode5

Mode4

Mode3

Mode2

Mode1

(a)

0 0

0.1

0.2

0.3

0.4

0.5

Time (s)

120Hz 0

1000

2000

3000

4000

5000

6000

2000

3000

4000

5000

6000

3000

4000

5000

6000

3000

4000

5000

6000

5000

6000

412Hz 0

1000

1003Hz 0

1000

2000

1389Hz 0

1000

2000

3639Hz 0

1000

2000

3000 4000 Frequency (Hz)

Fig.7. Results of TVMD for an inner race fault: (a) five modes, (b) corresponding frequency spectra.

Subsequently, TVMD is conducted, and the corresponding results are showed in Fig.7. The tentative process is terminated at the second time, which means that the calculated kurtosis is not obviously enhanced by TVMD, thus the mode number is 5. Then DTW is applied to choosing optimal modes of reconstructed signals. Fig. 8(a) presents the main modes mode3 and mode2 that are averaged as the 14

reconstructed signatures. Fig. 8(b) illustrates its envelope spectrum, showing that the two feature frequencies BPFI and 2BPFI are detectable. The results demonstrate that the improved method can extract the impulsive features in the heavy noise background. EMD is also used to process the inner race fault for comparison. In general, first six modes comprise more useful information. Fig. 9(a) shows the FFT spectra of the first six modes obtained by EMD. It can be noticed that the spectra of the latter modes present evident overlap, indicating that the fault features are decomposed into different modes. Consequently, it can be noticed in Fig. 9(b) that EMD fails to extract the fault feature frequencies. (a)

(b) Amplitude (m s -2 )

0.04 0.02 0 -0.02 -0.04

0

0.1

0.2

0.3

0.4

0.5

Time (s)

Amplitude

0.1

BPFI 0.05

0

2BPFI

0

100

200

300

400

500

600

Frequency (Hz)

Fig. 8. (a) principal modes, (b) the reconstructed signals and its envelope spectrum.

15

(b) A

3000

4000

5000

6000 A

2000

0

1000

2000

3000

4000

5000

6000 A

0.4 0.2 0

1000

0

0

1000

1000

2000

2000

3000

3000

4000

4000

5000

5000

6000

0.2 0.1 0 0.2 0.1 0

6000

A A

0.2 0.1 0

0

0

1000

1000

2000

2000

3000

4000

3000 4000 Frequency (Hz)

5000

5000

0.1 0.05 0 0.2 0.1 0

0.5 0

0.1 0.05 0

A

A

0.2 0.1 0

0

A

A

1 0.5 0

0.04 0.02 0

A

A

0.4 0.2 0

A

(a)

6000

6000

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300 400 Frequency (Hz)

500

600

Fig. 9. Results of EMD: (a) Frequency spectra of first six modes, (b) corresponding envelope spectra. (a)

(b)

Fig. 10. Original vibration signals for an outer race fault: (a) time waveform, (b) frequency spectrum.

4.2 Fault feature detection for an outer race fault The type of the bearings studied in the next two cases is SKF 6205. First, an outer race fault at drive end is analyzed. The fault point is 0.53 mm in diameter, and the shaft rotating speed is 1750 rpm. According to theoretical computations, the ball passing frequency outer race (BPFO) approaches 104.6 Hz. The original signals are sampled by the accelerometer installed on fan end and displayed in Fig. 10(a), where the modulation phenomena in amplitude is clear and the impulsive signals are 16

heavily smeared by noise. After TVMD is implemented, the samples are isolated into six modes illustrated in Fig. 11(a). In the tentative process, the number mode is defined as 3, 4, 5, 6 and 7, respectively. Correspondingly, the calculated maximal kurtosis value at each implementation of TVMD is 11.19, 11.96, 12.31, 17.38 and 17.16, respectively. It is observed that the obtained kurtosis has relatively good improvement when κ = 5. Thus, the tentative number is set to be 0 again, and then the new tentative process is performed again. Eventually, the optimal kurtosis 17.38 is obtained when κ = 6. Clearly, if the tentative operations are not conducted, the mode number will be defined by 4. In this case, it will result in the loss of the fault features. Therefore, it demonstrates that TVMD does have high reliability in this application. (a)

Mode3

0.1 0 -0.1 0.2 0 -0.2

A

Mode1 Mode2

0.5 0 -0.5

Mode4

(b) 0.05

0

0

0.1

0.1

0.1

0.2

0.3

0.2

0.3

0.2

0.3

0.2

0.3

0.2

0.3

0.4

0.4

0.4

0.4

0.4

0.4

0.5 A

0

0.1

0.3

0.4 0.2 0

A

0

0.1

0.2

0.4 0.2 0

A

0

0.1

0.4 0.2 0

0.5

0.5

0.5 A

0.5 0 -0.5

0

0.1 0.05 0

1 0.5 0

A

Mode5

0.1 0 -0.1

Mode6

0

0.2 0.1 0

0.5

0.5

Time (s)

156Hz 0

1000

2000

3000

4000

5000

6000

2000

3000

4000

5000

6000

2000

3000

4000

5000

6000

3000

4000

5000

6000

3000

4000

5000

6000

5000

6000

568Hz 0

1000

685Hz 0

1000

0

1000

1080Hz 2000

1400Hz 0

1000

2000

2730Hz 0

1000

2000

3000 4000 Frequency (Hz)

Fig. 11. Results of TVMD: (a) six modes, (b) corresponding frequency spectra.

17

(b) Amplitude (m s -2 )

(a)

0.2 0.1 0 -0.1 -0.2

0

0.1

0.2

0.3

0.4

0.5

Time (s) 0.2

Amplitude

0.15

BPFO

0.1

2BPFO 0.05 0

0

100

200

300

400

500

600

Frequency (Hz)

Fig. 12. (a) principal modes, (b) the reconstructed signals and its envelope spectrum.

After the similarity between the six modes is estimated by dynamic time warping, the principal modes comprise mode5 and mode6 that are illustrated in Fig. 12(a), and the reconstructed signal and its envelope spectrum are correspondingly showed in Fig. 12(b). In this envelope spectrum, the fault features BPFO and 2BPFO are obvious, and sidebands spaced at the shaft rotating frequency are detected. This means that the proposed approach detect the bearing damage successfully. On the other hand, complete ensemble empirical mode decomposition [41] has better robustness to noise than EMD and is utilized for comparison. Figs. 13(a) and 13(b) separately display the frequency and envelope spectra of the first six modes produced by CEEMD. It can be found that the fault characteristics are absent in Fig. 13(b). Therefore, it indicates that the CEEMD method can not diagnose the bearing fault.

18

(a)

0

1000

1000

2000

2000

2000

3000

3000

3000

4000

4000

4000

4000

3000 4000 Frequency (Hz)

5000

5000

5000

5000

5000

5000

A

0

1000

2000

3000

4000

A

0.2 0.1 0

0

1000

2000

3000

0.2 0.1 0

A

A

0.2 0.1 0

0

1000

2000

0.2 0.1 0

A

A

0.2 0.1 0

0

1000

0.2 0.1 0

A

A

0.2 0.1 0

0

0.2 0.1 0

0.2 0.1 0

A

A A

0.2 0.1 0

A

(b) 0.2 0.1 0

0.4 0.2 0

6000

6000

6000

6000

6000

6000

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300 Frequenc (Hz)

400

500

600

Fig. 13. Results of CEEMD: (a) frequency spectra of first six modes, (b) corresponding envelope spectrum.

4.3 Fault feature detection for an rolling element fault Compared with the inner and outer race faults, the rolling element one is diagnosed more difficultly. A main reason is that the fault point is always in motion status following the rolling elements rotating, thus the impulsive signals are subject to the interference from noise. Generally speaking, there are two shocks in a basic period, and the even harmonics of the feature frequency, i.e. ball-spin frequency (BSF), are quite often dominated. The related experimental parameters are: the diameter of the rolling element fault is 0.36 mm and the shaft speed is 1750 rpm. Correspondingly, the characteristic frequency BSF is about 68.7 Hz. Figs. 14(a) and 14(b) separately present the time waveform and frequency spectrum of the vibration signals.

19

(a)

(b)

Fig. 14. Original vibration signals for a rolling element fault: (a) time waveform, (b) frequency spectrum. (b) 0.1

1

0

0.5

0.1

0.2

0.3

0.4

0

0.5

0.1

0.4

0

0.2

0

1000

2000

3000

4000

5000

6000

2000

3000

4000

5000

6000

3000 4000 Frequency (Hz)

5000

6000

628Hz

-0.1

Mode3

0

A

Mode2

-0.1

A

360Hz

0

0.1

0.2

0.3

0.4

0

0.5

0.1

1

0

0.5

-0.1

A

Mode1

(a)

0

0.1

0.2

0.3

0.4

0.5

Time (s)

0

0

1000

1402Hz

0

1000

2000

Fig. 15. Results of TVMD: (a) three modes, (b) corresponding frequency spectra.

By the analysis of TVMD for the raw vibrations, it is found that the optimal mode number is 3. Figs. 15(a) and 15(b) show the time waveforms and frequency spectra of the three modes, respectively. According to Eq. (12), only mode2 illustrated in Fig. 16(a) acts as the reconstructed signals. In Fig. 16(b), the fault feature BSF and its second harmonic can be extracted. The results indicate that the proposed method is still efficient to diagnose the rolling element defect.

20

(a) Amplitude (m s -2 )

0.1 0.05 0 -0.05 -0.1

0

0.1

0.2

0.3

0.4

0.5

Time (s)

(b)

Amplitude

0.1 BSF 2BSF

0.05

0

0

100

200

300

400

500

600

Frequency (Hz)

Fig. 16. The reconstructed signals: (a) time waveform, (b) envelope spectrum.

In order to prove the superiority of the proposed approach, ensemble empirical mode decomposition (EEMD) is also applied. As an improved one, EEMD inherits the dyadic filter bank property of EMD. However, it does not well show this merit when decomposing the vibration signals, for instance, the frequency spectra of the latter four modes in Fig. 17(a) exhibit the obvious overlap behavior. As a result, the fault features are not exhibited in the corresponding envelope spectra. This means that the proposed method has stronger signal decomposition ability than EEMD.

(a)

(b)

Fig.17. Results of EEMD: (a) frequency spectra of first six modes, (b) corresponding envelope spectra. 21

4.4 Discussion In this study, EMD, EEMD and CEEMD are considered to process the vibration signals of the three faulty bearings, even if the latter two methods can improve the limitations of EMD effectively, they all do not obtain the satisfied diagnosis performance. The main reason is the feature of the decomposition procedure of the EMD-based methods, i.e. sifting local extrema. In strong noise background, the sifted extrema sometimes deviate the nature of original signals easily. Besides, backward error can not be corrected in the subsequent sifting procedures, thus generating the sensitivity of the three methods to noise and sampling. Thereby, they can not extract these bearing fault features. Differently, VMD shows a particular behavior, i.e. the frequency-domain sparseness under constraint in time-domain. In other words, it is a posterior signal decomposition method. Because the sampling process for a signal can not avoid the interference from background noise and sampling frequency to the signal-to-noise ration of the signal, VMD decompose a given signal by such consideration based on its whole frequency spectrum, namely allowing the sum of the basebands of all modes to be minimized by optimizing an objective function. The minimization makes the structure of the spectra of all modes simpler, that is, the main components in these spectra concentrate at their individual center frequencies. It can be known that if the frequency components of a signal can show the compacting feature in frequency domain, the signal itself (or the obtained modes) should be succinct in time domain because the frequency components of noise generally has wide range and seriously affects such compact performance in frequency domain. In other words, the decomposition of VMD shows the Winner filter denoising performance [20]. Besides, the feature frequencies are usually smaller than the center frequencies of the estimated modes, and the compact performance of VMD in frequency domain enables the frequency components around the feature frequencies to be sparser. The factors mentioned above works together, thus VMD can more effectively extract the fault features in strong noise background. TVMD is an improvement based on VMD, 22

enhancing the reliability in this decomposition process. At the same time, applying DWT can deal with the mode mixing problem effectively. As a result, the proposed method obtains a superior diagnosis to the three bearing faults when compared with the EMD-based methods. 5 Conclusions In order to detect fault features of faulty bearings, the combination of TVMD and DTW is presented in this study. The proposed method is evaluated by a simulated signal and the vibration signals sampled from the bearings separately with an outer, an inner and a rolling element faults. The results show that the proposed method can diagnose the three bearing faults. Conclusions are given as follows: (a) TVMD can determine the mode number adaptively and has high reliability in heavy noise background. (b) DTW is an effective way to solve the mode mixing problem of TVMD (c) In comparisons with EMD, EEMD and CEEMD, the proposed method has superior performance in extracting the bearing fault features. Acknowledgments This work was supported by National Natural Science Foundation of China (No. U1765201, No. 41571514) and Fundamental Research Funds for the Central Universities (No. 2017KFYXJJ204). References [1] Z. Chen, X. Yuan, B. Ji, P. Wang, H. Tian, Design of a fractional order PID

controller for hydraulic turbine regulating system using chaotic non-dominated sorting genetic algorithm II, Energy Conversion and Management 84 (2014) 390-404. [2] J. Liang, X. Yuan, Y. Yuan, Z. Chen, Y. Li, Nonlinear dynamic analysis and

robust controller design for Francis hydraulic turbine regulating system with a straight-tube surge tank, Mech. Syst. Signal Process. 85 (2017) 927-946. [3] X. Yuan, B. Ji, Y. Yuan, R. Ikram, X. Zhang, Y. Huang, An efficient chaos

embedded hybrid approach for hydro-thermal unit commitment problem, Energy Conversion and Management 91 (2015) 225-237. 23

[4] N.E Huang, Z. Shen, S.R. Long, M. Wu, H. Shih, Q. Zheng, N. Yen, C. Tung, H.

Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 454 (1971) (1998) 903-995. [5] Y. Lv, R. Yuan, G. Song, Multivariate empirical mode decomposition and its

application to fault diagnosis of rolling bearing, Mech. Syst. Signal Process. 81 (2016) 219-234. [6] J. Dybala, R. Zimroz, Rolling bearing diagnosing method based on empirical

mode decomposition of machine vibration signal Appl. Acoust. 77 (2014) 195-203. [7] M. Amarnath, I.P. Krishna, Empirical mode decomposition of acoustic signals

for diagnosis of faults in gears and rolling element bearings, IET. Sci. Mea. Technol. 6 (2012) 279-287. [8] M. Li, F. Li, B. Jing, H. Bai, H. Li, G. Meng, Multi-fault diagnosis of rotor

system based on differential-based empirical mode decomposition, J. Vib. Control. 21(9) (2015) 1821-1837. [9] G. Rilling, P. Flandrin, Sampling effects on the empirical mode decomposition,

Adv. Adapt. Data Anal. 1(2009) 43-59. [10] X. Yuan, Q. Tan, X. Lei, Y. Yuan, X. Wu, Wind power prediction using hybrid

autoregressive fractionally integrated moving average and least square support vector machine, Energy 129 (2017) 122-137. [11] G. Rilling, P. Flandrin, One or two frequencies? The empirical mode

decomposition answers, IEEE. T. Signal. Process. 56(2008) 85-95. [12] Z. Wu, N.E. Huang, Ensemble empirical mode decomposition: a noise-assisted

data analysis method, Adv. Adapt. Data Anal. 1(2009) 1-41. [13] S. Wang, N. Zhang, L. Wu , Y. Wang, Wind speed forecasting based on the

hybrid ensemble empirical mode decomposition and GA-BP neural network method, Renew. Energy 94 (2016) 629-636 [14] X. Zhang, Y. Liang, J. Zhou, Y. Zang, A novel bearing fault diagnosis model

integrated permutation entropy, ensemble empirical mode decomposition and optimized, Measurement 69 (2015) 164-179. [15] L. Meng, J. Xiang, Y. Wang, Y. Jiang, H. Gao, A hybrid fault diagnosis method

using morphological filter–translation invariant wavelet and improved ensemble 24

empirical mode decomposition ,Mech. Syst. Signal Process.50 (2015)101-115. [16] Y. Lei, N. Li, J. Lin, S. Wang, Fault diagnosis of rotating machinery based on

an adaptive ensemble empirical mode decomposition, Sensors 13(12) (2013) 16950-16964. [17] X. Yuan, Z. Chen, Y. Yuan, Y. Huang, Design of fuzzy sliding mode controller

for hydraulic turbine regulating system via input state feedback linearization method, Energy 93 (2015) 173-187. [18] J. Zhang, R. Yan, R. Gao, Z. Feng, Performance enhancement of ensemble

empirical mode decomposition, Mech. Syst. Signal Process 24 (2010)2104-2123 [19] J. Zheng, J. Cheng, Y. Yang, Partly ensemble empirical mode decomposition:

an improved noise-assisted method for eliminating mode mixing, Signal Process. 96 (2014) 362-374. [20] K. Dragomiretskiy, D. Zosso, Variational mode decomposition, IEEE. T. Signal

Process. 62 (2014) 531-544. [21] K.K. Gupta, K.S. Raju, Bearing fault analysis using variational mode

decomposition, 9th International Conference on Industrial and Information Systems (ICIIS), Gwalior, India, Dec. 15-17, 2014, pp. 1-6. [22] S. Zhang, Y. Wang, S. He, Z. Jiang, Bearing fault diagnosis based on variational

mode decomposition and total variation denoising, Meas. Sci. Technol. 27(7) (2016) DOI: 10.1088/0957-0233/27/7/075101 [23] Y. Wang, R.Markert, J. Xiang, W. Zheng, Research on variational mode

decomposition and its application in detecting rub-impact fault of the rotor system, Mech. Syst. Signal. Process. 60 (2015) 243-251. [24] C. Yi, Y. Lv, Z. Dang, A fault diagnosis scheme for rolling bearing based on

particle swarm optimization in variational mode decomposition. Shock Vib. 2016 (2016). [25] Z. Lv, B. Tang, Y. Zhou, C. Zhou, A novel method for mechanical fault

diagnosis based on variational mode decomposition and multikernel support vector machine, Shock Vib. (2016) DOI: 10.1155/2016/3196465. [26] Z. Li, Y. Jiang, Q. Guo, C. Hu, Z. Peng, Multi-dimensional variational mode

decomposition for bearing-crack detection in wind turbines with large driving-speed variations, Renew. Energ. 116 (2018) 55-73. [27] J. Zhu, C. Wang, Z. Hu, F. Kong, X. Liu, Adaptive variational mode

decomposition based on artificial fish swarm algorithm for fault diagnosis of 25

rolling bearings, Proceedings of the Institution of Mechanical Engineers, Part C: P. I. Mech. Eng. C-J. Mec. 231(4) (2017) 635-654. [28] Y. Liu, G. Yang, M. Li, H. Yin, Variational mode decomposition denoising

combined the detrended fluctuation analysis, Signal Process. 125 (2016) 349-364. [29] Z. Li, J. Chen, Y. Zi, J. Pan, Independence-oriented VMD to identify fault

feature for wheel set bearing fault diagnosis of high speed locomotive, Mech. Syst. Signal. Process. 85 (2017) 512-529. [30] J. Liang, Z. Liu, H. Wang, X. Dong, Adaptive variational mode decomposition

method for signal processing based on mode characteristic, Mech. Syst. Signal. Process. 107 (2018) 53-77. [31] Y. Dong, M. Liao, X. Zhang, F. Wang, Faults diagnosis of rolling element

bearing based on modified morphological method, Mech. Syst. Signal Process. 25(4) (2011) 1276-1286. [32] A.S. Raj, N. Murali, Early classification of bearing faults using morphological

operators and fuzzy inference, IEEE T. Ind. Electron. 60(2) (2013) 567-574. [33] H. Sakoe, S. Chiba, Dynamic programming algorithm optimization for spoken

word recognition, IEEE .T. Acoust. Speech. Signal Process. 26 (1978) 43-49. [34] J. Tang, H. Cheng, Y. Zhao, H. Guo, Structured dynamic time warping for

continuous hand trajectory gesture recognition, Pattern Recogn. 80 (2018) 21-31. [35] D.J. Berndt, J. Clifford, Using dynamic time warping to find patterns in time

series, KDD. workshop. 10 (1994) 359-370. [36] Y.X. Wang, R. Markert, Filter bank property of variational mode decomposition

and its applications, Signal Process. 120 (2016) 509–521. [37] A. Lempel, J. Ziv, On the complexity of finite sequences, IEEE Trans. Inf.

Theory. 22 (1976) 75-81. [38] F. Kaspar, H.G. Schuster, Easily calculable measure for the complexity of

spatiotemporal patterns, Phys. Rev. 36 (1987) 842–848. [39] C. Bandt, B. Pompe, Permutation entropy: a natural complexity measure for

time series, Phys. Rev. Lett. 88 (2002) [40] Case Western Reserve University Bearings Vibration Dataset available:

http://csegroups.case.edu/bearingdatacenter/home (Accessed: October 2013). [41] M. Torres, M. Colominas, G. Schlotthauer, P. Flandrin, A complete ensemble 26

empirical mode decomposition with adaptive noise, 2011 IEEE International Conference on Acoustics, Speech and Signal Processing, Prague Congress Ctr, Prague, Czech Republic, May 22-27, 2011, pp. 4144-4147.

● TVMD is proposed to detect fault feature of rolling element bearing. ● TVMD has good adaptivity and reliability in heavy noise signal. ● Selection of principal mode based DTW is sensitive to bearing fault feature. ● The proposed TVMD has superior behavior in fault feature detection.

27