Bearing Fault Diagnosis Based on Optimal Time-Frequency Representation Method⁎

5th IFAC International Conference on 5th IFAC International Conference on Intelligent Control andConference Automationon Sciences 5th IFAC IFAC International International 5th Intelligent Control andConference Automationon Sciences Available online at www.sciencedirect.com Belfast, Northern Ireland, August 21-23, 2019 5th IFAC International Conference on Intelligent Control and Automation Sciences Intelligent ControlIreland, and Automation Sciences Belfast, Northern August 21-23, 2019 Intelligent Control and Automation Sciences Belfast, Northern Ireland, August 21-23, 2019 Belfast, Northern Ireland, August 21-23, 2019 Belfast, Northern Ireland, August 21-23, 2019

ScienceDirect

IFAC PapersOnLine 52-11 (2019) 194–199

Bearing Bearing Fault Fault Diagnosis Diagnosis Based Based on on Optimal Optimal ⋆⋆ Bearing Fault Diagnosis Based on Optimal Time-Frequency Representation Bearing Fault Diagnosis Based onMethod Optimal Time-Frequency Representation Method Time-Frequency Representation Method ⋆⋆ Time-Frequency Representation Method Israel Ruiz Quinde, Jorge Chuya Sumba, ∗∗

{ { { { {

Israel Ruiz Quinde, Jorge Chuya Sumba, ∗∗ ∗ Israel Quinde, Jorge Chuya Sumba, Luis Ochoa, Antonio Vallejo Guevara, ∗ Israel Ruiz Ruiz Quinde, Jorge Jr. Chuya Sumba, Luis Escajeda Escajeda Ochoa, Antonio Jr. Vallejo Guevara, ∗ ∗ ∗ Israel Ruiz Quinde, Jorge Chuya Sumba, Luis Escajeda Ochoa, Antonio Jr. Vallejo Guevara, Ruben Morales-Menendez ∗ Luis Escajeda Ochoa, Antonio Jr. Vallejo Guevara, ∗∗ Ruben Morales-Menendez ∗ Luis Escajeda Ochoa, Antonio Jr. Vallejo Guevara, ∗ Ruben Morales-Menendez Ruben Morales-Menendez ∗ ∗ Ruben Morales-Menendez o gico de Monterrey, Monterrey NL 64,489 ∗ Tecnol´ ogico de Monterrey, Monterrey NL 64,489 M´ M´eexico xico ∗ Tecnol´ ∗ Tecnol´ o gico de Monterrey, Monterrey NL 64,489 M´eexico xicormm A00824089, A00824119, A00823938 } @itesm.mx, { avallejo, Tecnol´ o gico de Monterrey, Monterrey NL 64,489 M´ A00824089, A00824119, A00823938 } @itesm.mx, { avallejo, ∗ Tecnol´ oA00824119, gico de Monterrey, Monterrey NL 64,489 M´exicormm A00824089, A00823938 } @itesm.mx, { avallejo, @tec.mx A00824089, A00824119, A00823938 } @itesm.mx, { avallejo, rmm @tec.mx} @itesm.mx, { avallejo, rmm A00824089, A00824119, A00823938 rmm @tec.mx @tec.mx @tec.mx

}} }} }

Abstract: Abstract: Abstract: Wigner-Ville is probably probably the the most most used used non-linear non-linear time-frequency time-frequency distridistriAbstract: Distribution Wigner-Ville Distribution (WVD) (WVD) is Abstract: Wigner-Ville Distribution (WVD) is used non-linear time-frequency distribution for signal signal processing in fault fault diagnosis,the duemost to the the advantages of excellent resolution resolution Wigner-Ville Distribution (WVD) is probably probably the most usedadvantages non-linear of time-frequency distribution for processing in diagnosis, due to excellent Wigner-Ville Distribution (WVD) is probably the most used non-linear time-frequency distribution for signal processing in fault diagnosis, due to the advantages of excellent resolution and localization in time-frequency domain. However, the presence of cross terms when they are bution for signal processing in fault diagnosis, due to the advantages of excellent resolution and localization in time-frequency domain. However, the presence of cross terms when they are bution for signalin processing in fault diagnosis, due the to the advantages of excellent resolution and localization in time-frequency domain. However, the presence of cross terms when they are applied to multicomponent signals can give misleading interpretations. A methodology based and localization time-frequency domain. However, presence of cross terms when they are applied to multicomponent signalsdomain. can giveHowever, misleading methodology based and localization in time-frequency the interpretations. presence of more crossA terms when they are applied to multicomponent signals can give misleading interpretations. A methodology based on Local Mean Decomposition (LMD) and WVD is proposed proposed to get get reliable bearing fault applied to multicomponent signals canand give misleading interpretations. Areliable methodology based on Local Mean Decomposition (LMD) WVD is to more bearing fault applied to multicomponent signals can give misleading interpretations. A methodology based on Local Mean Decomposition (LMD) and WVD is proposed to get more reliable bearing fault diagnosis based on vibration signals. Kullback-Leibler Divergence (KLD) guides the selection of on Local Mean (LMD)Kullback-Leibler and WVD is proposed to get moreguides reliable fault diagnosis basedDecomposition on vibration signals. Divergence (KLD) thebearing selection of on Local Mean Decomposition (LMD) and WVD is proposed to get more reliable bearing fault diagnosis based on vibration signals. Kullback-Leibler Divergence (KLD) guides the selection of the optimal frequency band with the most relevant information about the fault. Early results diagnosis based on vibration signals. Kullback-Leibler Divergence (KLD) guides the selection of the optimal frequency band with theKullback-Leibler most relevant information fault. results diagnosis based on vibration signals. Divergence about (KLD)the guides theEarly selection of the optimal frequency band with the most relevant information about the fault. Early results based on experimental data show successful diagnosis. the optimal frequency band with the most relevant information about the fault. Early results based on experimental data show successful diagnosis. the optimal frequency band with the most relevant information about the fault. Early results based on experimental data show successful diagnosis. based on experimental data show successful diagnosis. © 2019,on IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. based experimental show successful diagnosis. Keywords: Wigner-Villedata Distribution, Fault Diagnosis, Diagnosis, Bearing Spindles. Spindles. Keywords: Wigner-Ville Distribution, Fault Bearing Keywords: Keywords: Wigner-Ville Wigner-Ville Distribution, Distribution, Fault Fault Diagnosis, Diagnosis, Bearing Bearing Spindles. Spindles. Keywords: Wigner-Ville Distribution, Fault Diagnosis, Bearing Spindles. 1. INTRODUCTION Table 1. INTRODUCTION Table 1. 1. Definition Definition of of acronyms acronyms 1. Table 1. 1. Definition Definition of of acronyms acronyms 1. INTRODUCTION INTRODUCTION Table 1. INTRODUCTION 1. Definition of acronyms Acronym Table Definition The conditions conditions under under which which rotating rotating machines machines work, work, play play Acronym Definition The Acronym Acronym Definition Definition The conditions under which machines work, AM Amplitude Modulation an important role in machinery machinery diagnosis. There is an anplay inTheimportant conditionsrole under which rotating rotating machines work, play AM Amplitude Modulation an in diagnosis. There is inAcronym Definition The conditions under which rotating machines work, play ANN Artificial Neural Network AM Amplitude Modulation an important role in machinery diagnosis. There is an increased requirement for effective/efficient techniques that AM Amplitude Modulation an important role in for machinery diagnosis. techniques There is anthat inANN Artificial Neural Network creased requirement effective/efficient AOK Adaptive Optimal Kernel AM Amplitude Modulation an important role in in for machinery diagnosis. Therethe is an inANN Artificial Neural Network creased requirement effective/efficient techniques that monitor machines real time. This could allow operANN Artificial Neural Network AOK Adaptive Optimal Kernel creased requirement for effective/efficient techniques that monitor machines in real time. This could allow the operCWRU Case Western Reserve University ANN Artificial Neural Network AOK Adaptive Optimal Kernel creased requirement for effective/efficient techniques that AOK Adaptive Optimal Kernel monitor machines in real time. This could allow the operator to schedule maintenance before the defect ends in a CWRU Case Western Reserve University monitor machines in real time. This could allow the operator to schedule maintenance before the defect ends in a EK Excess of Optimal Kurtosis AOK Adaptive Kernel CWRU Case Western Reserve University monitor machines in real time. This could allow the operCWRU Case Western Reserve University ator to schedule maintenance before the defect ends in a EK Excess of Kurtosis full damage of the machine and therefore an unexpected ator damage to schedule maintenance before the defect ends in a full of the machine and therefore an unexpected EMD Empirical Mode Decomposition CWRU Case Western Reserve University EK Excess of Kurtosis EK Excess of Kurtosis ator to schedule maintenance before the defect ends in a EMD Empirical Mode Decomposition full damage of the machine and therefore an unexpected production downtime with high maintenance costs. Spinfull damagedowntime of the machine and maintenance therefore an costs. unexpected FM Frequency Modulation EK Excess of Kurtosis production with high SpinEMD Empirical Mode Decomposition EMD Empirical Mode Decomposition FM Frequency Modulation full damage of the machine and therefore an unexpected production maintenance costs. dle bearingsdowntime are the the with most high critical and vulnerable vulnerable comproduction downtime with high maintenance costs. SpinSpinHSM High SpeedMode Machining EMD Empirical Decomposition FM Frequency Modulation dle bearings are most critical and comFM Frequency HSM High SpeedModulation Machining production downtime with high maintenance costs. Spindle bearings are the most critical and vulnerable components in High Speed Machining (HSM). Friction, load dle bearings are the most critical and vulnerable comIR Inner Race Modulation FM Frequency HSM High Speed Machining ponents in High Speed Machining (HSM). Friction, load HSM High Speed Machining IR Inner Race dle bearings are Speed the actuating most critical and vulnerable components in High Machining Friction, load forces and vibrations over(HSM). bearings can produce produce KLD Kullback - Leibler Divergence HSM High Speed Machining IR Inner Race ponents in vibrations High Speed Machining (HSM). Friction, load IR Inner Race forces and actuating over bearings can KLD Kullback - Leibler Divergence ponents in vibrations High Speed Machining (HSM). Friction, load forces and actuating over bearings can produce LMD Local Mean Decomposition IR Inner Race wear, fatigue and impending cracks on these. A survey for KLD Kullback Leibler Divergence forces and vibrations actuating over bearings can produce KLD Kullback - Leibler Divergence LMD Local Mean Decomposition wear, fatigue and impending cracks these. A survey for forces andmotors vibrations actuating overon bearings can produce OI Orthogonality Index KLD Kullback - Leibler Divergence LMD Local Mean Mean Decomposition wear, fatigue and impending cracks on these. A survey for induction reports that about 40 % of their failures LMD Local Decomposition wear, fatigue and impending cracks on these. A survey for OI Orthogonality Index induction motors reports that about 40 % of their failures OR Outer Race LMD Local Mean Decomposition wear, fatigue and impending cracks on these. A survey for OI Orthogonality Index induction motors reports that about 40 % of their failures are relatedmotors to bearings, bearings, Hyers et al. al. (2006). (2006). OI Orthogonality Index OR Outer Race induction reportsHyers that about 40 % of their failures are related to et PF Profuct Function OI Orthogonality Index OR Outer Race Race induction motors reportsHyers that about 40 % of their failures OR Outer are related to bearings, et al. (2006). PF Profuct Function are related to bearings, Hyers et al. (2006). Condition based monitoring of rotating machines can be RE Rolling Elements OR Outer Race PF Profuct Function are related to bearings, Hyers et al. (2006). PF Profuct Elements Function Condition based monitoring of rotating machines can be RE Rolling RMS Root Mean Square PF Profuct Function Condition monitoring of machines be RE Rolling Elements achieved bybased the effective effective analysis of vibrations, vibrations, wearcan debris Conditionby based monitoring of rotating rotating machines can be RE Rolling Elements RMS Root Mean Square achieved the analysis of wear debris Condition based monitoring of rotating machines can be SET Synchro Extracting RE Rolling Elements RMS Root Mean Square Transform achieved by the effective analysis of vibrations, wear debris in oil, acoustic emission, etc, which are recorded from the RMS Root Mean Square achieved by the emission, effective analysis of vibrations, wear debris SET Synchro - Extracting Transform in oil, acoustic etc, which are recorded from the SPWVD Smoothed Pseudo - Ville Distribution RMS Root Mean SquareWigner SET Synchro -- Extracting Transform achieved by the emission, effective analysis of vibrations, wear debris in oil, acoustic etc, which are recorded from the SET Synchro Extracting Transform machine. SPWVD Smoothed Pseudo Wigner - Ville Distribution in oil, acoustic emission, etc, which are recorded from the machine. TFR Time Frequency Representation SET Synchro Extracting Transform SPWVD Smoothed Pseudo Wigner Ville Distribution Distribution in oil, acoustic emission, etc, which are recorded from the SPWVD Smoothed Pseudo Wigner -- Ville machine. TFR Time - Frequency Representation machine. uLBP Uniform Local Binary Patterns SPWVD Smoothed Pseudo Wigner - Ville Distribution The limitation which which posses posses linear linear Time-Frequency Time-Frequency RepRepTFR Time -- Frequency Representation TFR Time Frequency Representation machine. uLBP Uniform Local Binary Patterns The limitation VMD Variational Mode Decomposition TFR Time - Frequency Representation uLBP Uniform Local Binary Patterns The which posses linear Time-Frequency Representation (TFR) such as Short Short Time Fourier Transform Transform uLBP Uniform Local Binary Patterns The limitation limitation which posses linear Time-Frequency RepVMD Variational Mode Decomposition resentation (TFR) such as Time Fourier WT Wavelet Transform uLBP Uniform Local Binary Patterns The limitation which posses linear Time-Frequency RepVMD Variational Mode Decomposition resentation (TFR) such as Short Time Fourier Transform (STFT) and Wavelet Transform (WT), regarding to the VMD Variational Mode Decomposition resentation (TFR) such as Short Time Fourier Transform WT Wavelet Transform (STFT) and(TFR) Wavelet Transform (WT), regarding to the WVD Wigner -Transform Ville Distribution VMD Variational Mode Decomposition WT Wavelet Transform resentation such as Short Time Fourier Transform WT Wavelet (STFT) and Wavelet Transform (WT), regarding to the selection of a suitable window size for effective signal analWVD Wigner - Ville Distribution (STFT) and Wavelet window Transform regarding the selection of a suitable size (WT), for effective signalto analWT Wavelet WVD Wigner Ville (STFT) Wavelet Transform (WT), regarding the WVD Wigner --Transform Ville Distribution Distribution selection of aa suitable size for analysis, has and aroused the window interest of researchers researchers insignal the to appliselection of suitable window size for effective effective signal analysis, has aroused the interest of in the appliWVD Wigner - Ville Distribution selection of a suitable window size for effective signal analysis, has aroused the interest of researchers in the application of non-linear distributions. Quadratic distributions, ysis, has aroused the interest of researchers in the application of non-linear distributions. Quadratic distributions, ysis, aroused the interest researchers in the application non-linear distributions. Quadratic distributions, such asof WVD offers very highofresolution resolution regarding good Du et al. (2016) proposed a method based on the smoothed cationhas ofWVD non-linear distributions. Quadratic distributions, such as offers very high regarding good Du et al. (2016) proposed a method based on the smoothed cation of non-linear distributions. Quadratic distributions, such as WVD offers very high resolution regarding good location of the signal in the time-frequency plane. HowHowsuch as WVD offers very high resolution regarding good version Du et aa method based on WVD to the of compolocation of the signal in the time-frequency plane. Du et al. al.of (2016) proposed method based on the the smoothed smoothed version of(2016) WVDproposed to obtain obtain the TFR TFR of resultant resultant composuch as every WVD offers high resolution good nents location of signal in the time-frequency plane. ever, as non-linear distribution, their regarding representations Du et al. (2016) proposed aDecomposition method based on the smoothed location of the the signalvery indistribution, the time-frequency plane. HowHowversion of WVD to obtain the TFR of resultant compoin Empirical Mode (EMD) method, ever, as every non-linear their representations version of WVD to obtain the TFR of resultant components in Empirical Mode Decomposition (EMD) method, location of the signal in the time-frequency plane. However, non-linear distribution, their representations present interference terms, Cohen (1989). (1989). To address this this and version ofEmpirical WVD to obtain the TFR ofthe resultant compoever, as as every every non-linear distribution, their To representations nents in Mode Decomposition (EMD) method, by intuition, define the zone with greatest lumipresent interference terms, Cohen address nents in Empirical Mode Decomposition (EMD) method, and byinintuition, define the zone with the greatest lumiever, as every non-linear distribution, their To representations present interference terms, Cohen (1989). address this issue, there are some techniques that combined with WVD nents Empirical Mode Decomposition (EMD) method, present interference terms, Cohen (1989). To address this and by intuition, define the zone with the greatest luminosity in the energy distribution as the resonance band. issue, there are some techniques that combined with WVD and by intuition, define the zone with the greatest luminosity in the energy distribution as the resonance band. present interference terms, Cohen (1989). To address this issue, there are some techniques that combined with WVD increases its effectiveness. byin intuition, define the zone with the greatest lumiissue, there are some techniques that combined with WVD and nosity in the energy distribution as the resonance band. Other signal decomposition method commonly used in increases its effectiveness. nosity the energy distribution as the resonance band. Other signal decomposition method commonly used in issue, there some techniques that combined with WVD nosity increases its effectiveness. in the energy distribution as the resonance band. increases itsare effectiveness. Other signal decomposition method commonly used in fault diagnosis is Local Mean Decomposition (LMD), Xie Other signal decomposition method commonly used in fault diagnosis is Local Mean Decomposition (LMD), Xie increases its effectiveness. Other signal In decomposition method commonly usedXie in fault diagnosis is Local Mean Decomposition (LMD), et al. (2012). Liu et al. (2017), important contributions fault diagnosis is Local Mean Decomposition (LMD), Xie et al. diagnosis (2012). InisLiu et al. (2017), important contributions ⋆ Authors thank Tecnol´ fault Local Mean Decomposition (LMD), Xie et al. (2012). In Liu et al. (2017), important contributions were exposed regarding setting parameters that must be o gico de Monterrey for its support. ⋆ Authors thank Tecnol´ et al. exposed (2012). In Liu et al.setting (2017),parameters important that contributions were regarding must be ogico de Monterrey for its support. ⋆ et al. exposed (2012). In Liu et al.setting (2017),parameters important that contributions ⋆ Authors were exposed regarding setting parameters that must be be Authors thank thank Tecnol´ Tecnol´ ogico gico de de Monterrey Monterrey for for its its support. support. were regarding must o ⋆ Authors thank Tecnol´ were exposed regarding setting parameters that must be gico de Monterrey for its 2405-8963 © 2019, IFAC o(International Federation of support. Automatic Control) Hosting by Elsevier Ltd. All rights reserved.

Copyright © 2019 IFAC 194 Copyright 2019 responsibility IFAC 194Control. Peer review©under of International Federation of Automatic Copyright © 2019 IFAC 194 Copyright © 2019 IFAC 194 10.1016/j.ifacol.2019.09.140 Copyright © 2019 IFAC 194

2019 IFAC ICONS Israel Ruiz Quinde et al. / IFAC PapersOnLine 52-11 (2019) 194–199 Belfast, Northern Ireland, August 21-23, 2019

considered in LMD to obtain reliable results in fault diagnosis. Ma et al. (2018) introduced a method based on Variational Mode Decomposition (VMD) and the Synchro-Extracting Transform (SET) for improving time-frequency analysis and estimating accurately the instantaneous frequency of signals of rolling bearing; experimental results showed an error of 0.98 %. Chen et al. (2017) proposed a method for fault diagnosis using texture features extracted from TFR of signals, by using Adaptive Optimal Kernel (AOK ) function and uniform Local Binary Patterns (uLBP).

4. The signal s11 (n) may not be a pure FM after the first decomposition, thus repeat procedures from step 1 to 3 for r times until the Smith criterion, eqn (8), is fulfilled and a purely FM signal s1n (t) is obtained. This is called the sifting process: h11 (n) = x(n) − m11 (n), h12 (n) = s11 (t) − m12 (n), .. . h1r (n) = s1(r−1) (n) − m1r (n).

(6)

s11 (n) = h11 (n)/a11 (n), s12 (n) = h12 (n)/a12 (n), .. . s1r (n) = h1r (n)/a1r (n).

(7)

lim a1r (n) = 1

(8)

This paper is organized as follows. Section 2 presents the theoretical background. Section 3 shows the methodology. Section 4 discusses experimental results. Finally, section 5 concludes the research work. Tables 1 and 2 summarize the definitions of used acronyms and variables. 2. THEORETICAL BACKGROUND

n→∞

Time-frequency analysis methods such as LMD and WVD emerged as a powerful tool to address these problems. 2.1 Rolling Bearings Some possible defects in bearings are: Inner/Outer Races (IR, OR), Rolling Elements (RE) and Cage defects. When RE strike with a local fault on races, a shock caused by impact is introduced exciting high-frequency resonances of the whole structure, Smith and Randall (2015). The characteristic frequencies are computed according to the dimensions of the bearings as: ORf = nα (1 − β) ) ( 2 REf = αD 1 − (β) /d

195

IRf = nα (1 + β)

(1)

Cagef = α (1 − β)

(2)

where the subscripts in hir (n), mir (n), sir (n), and air (n) refer to the ith PF and j th sifting process. 5. Compute the corresponding envelope signal a1 (n) of s1n (n): r ∏ a1 (n) = a11 (n)a12 (n)...a1r (n) = a1k (n) (9) k=1

6. Calculate the first PF:

(10) PF1 (n) = a1 (n)s1n (n) 7. The next PF are formed using a residual signal u(n), the whole process is repeated q-times until u(n) is a constant or contains no more oscillations. u1 (n) = x(n) − PF1 (n), u2 (n) = u1 (n) − PF2 (n), .. . uq (n) = uq−1 (n) − PFq (n).

where α = fr /2 and β = (d/D)cosϕ.

(11)

P Fq (n) and their envelopes, aq (n), can be used for performing the fault diagnosis in high/low frequencies.

2.2 LMD LMD uses an iterative approach to decompose multicomponent Amplitude-Frequency Modulation (AM-FM) signals into a series of Product Functions (PF), preserving their local features, Smith (2005). The LMD algorithm is: 1. For a discrete time signal x(n) of length N , find the local extremes, compute the raw version of local mean, m0 (n), and local magnitude, a0 (n), by the mean and the difference of two successive extrema nl and nl+1 in x(n) respectively: m0 (n) = [x(nl ) + x(nl+1 )]/2;

l ∈ 1, 2, ...N

(3)

a0 (n) = [x(nl ) − x(nl+1 )]/2; l ∈ 1, 2, ...N (4) Due to m0 (n) and a0 (n) signals are stepped waves, they are smoothed by using a moving averaging and new m(n) and a(n) are formed. 2. Compute a zero-mean signal: h11 (n) = x(n) − m11 (n) 3. Get the estimated FM signal s11 (t): s11 (n) = h11 (n)/a11 (n)

(5) 195

2.3 WVD For a signal x(t) the WVD is defined as: ∫ ∞ ( ) ( ) WVDx (t, f ) =

−∞

x t−

τ 2

x∗ t +

τ 2

e−j2πf τ dτ

(12)

WVD is the Fourier Transform (FT) of the local autocovariance or autocorrelation function of x(t). It is a quadratic energy distribution due to the product of the analyzed signal with its complex conjugate x∗ , eqn (12).

To suppress cross terms, Smoothed Pseudo WVD (SPWVD) has been proposed, Boashash (2016). This uses independent windows to smooth WVD in time and frequency. ∫ ∞ ( ) ( ) SPVWDx (t, f ) =

g(t)H(f )x t +

−∞

τ 2

x∗ t −

τ 2

e−j2πf τ dτ

(13)

Smoothed version of WVD can be applied to identify changes in energy signal and thus to detect the fault components with a good localization in frequency, Ibrahim and Albarbar (2011).

2019 IFAC ICONS 196 Israel Ruiz Quinde et al. / IFAC PapersOnLine 52-11 (2019) 194–199 Belfast, Northern Ireland, August 21-23, 2019

2.4 KLD Relative entropy or KLD measures the distinguish ability of two probability distributions: p and q, Rold´an (2014): D [p(x) ∥ q(y)] =

∫

∞

p(x)ln

−∞

p(x) dx q(y)

(14)

KLD is always a positive value. KLD can be used as an indicator of the correlation between each PF and the original signal. The smaller KLD value, a high correlation between the PF and original signal, Si et al. (2017). Table 2. Definition of variables Variable

Definition

ϕ θ a(t) a1 (t) a11 (t) Am (t) Ax (t) d D f fr g(t) h11 (t) H(f ) Lg Lh m(t) m11 (t) n N OI p(x), q(x) P Fn (t) Pi q, i s11 (t) u1 (t), u2 (t)...uq (t) x(t)

Contact angle measures from the center line of the RE and bearing axis Threshold for LM D (0.005) Local magnitude (after smoothing) First AM signal Initial envelope estimation Average amplitude of the fault component Amplitude of estimated cross terms Number of the RE Pitch circle diameter Objective function Shaft Speed of rotating machine Window function in time domain First zero-mean signal Window function in frequency domain Length of g(t) expressed in # of samples in time Length of H(f ) expressed in # of frequency points Local mean (after smoothing) First computation of local mean Discrete time Length of x(t) and P F expressed in # of samples Orthogonal Index Probability distributions of x and y nth P F Percent of energy of each TFR block Indexes of P F Initial F M signal estimation Residual signals in LM D Raw signal

OI = ∑Ns

1

n=1

x2 (n)

q ∑ i ∑ N ∑

i=1 j=1 n=1

P Fi (n) × P Fj (n)

(16)

To verify the effect of proposed objective function, eqn (15), this was tested with experimental data. The modified stopping criteria is adopted to decompose vibration signals and minimize mode mixing between resultant PF. 2.6 Smoothing the WVD The effect of cross terms in WVD can be reduced by applying the proper window function (size and type). The maximum shaft speed of the experiments provides the highest values of fault frequencies and therefore the minimum length of window function, that has to be used to avoid over-smoothing in the WVD. Sometimes energy of the fault components is dispersed in the sidebands which appear as a result of modulation effects that occurs when fault frequencies are modulated by bearing housings natural frequency, Sinha (2015) . In case of IR signals, sidebands are spaced at fr regarding IR component and its harmonics. In RE signals, sidebands are spaced at Cagef , eqn (2), regarding RE component and its harmonics. OR does not show sidebands, Smith and Randall (2015). To select the proper window length in frequency, IR and RE signals were analyzed to avoid suppressing these sidebands. The fault frequencies in number of samples were computed for each baring condition using using fault frequencies. A signal length of 2,048 points and a frequency bin of 2.93 Hz, Table 3. Table 3. Window functions in frequency. Shaft Speed in RPM

Shaft Speed (Hz)

1,797

29.35

OR 107

Fault Frequency (Hz) IR+fr RE +Cagef 196 156

Samples in frequency OR IR RE 37 69 55

2.5 LMD Optimal Parameters An important limitation that LMD technique can suffer is the mode mixing problem. The eqn (8) describes the ideal condition which must be satisfied to stop the envelope estimation in LMD; but, it can result a high computational cost. The objective function, Liu et al. (2017), to stop the sifting procedure, considers two parameters, the Root Mean Square (RMS ) and the Excess of Kurtosis (EK ) computed over a zero-baseline envelope signal, that is obtained by subtracting 1 from a1n (t). A small variation of the aforementioned objective function is proposed, eqn (15), with a threshold (θ = 0.005). f = RM S(x(n)) × EK(x(n)),

(15)

where x(n) is the analyzed signal. The Orthogonality Index (OI) verifies if obtained PF are mutually orthogonal, and it is defined as:

196

Fig. 1. Window functions in frequency with (Lh). Figure 1 exemplifies a smoothed version of WVD computed for an envelope of the PF1 extracted from a experimental signal, with a shaft speed of 1,730 RPM. Hanning windows in frequency of sizes (Lh) defined in Table 3, were

2019 IFAC ICONS Israel Ruiz Quinde et al. / IFAC PapersOnLine 52-11 (2019) 194–199 Belfast, Northern Ireland, August 21-23, 2019

testing. Sidebands in IR signal were completely suppressed when Lh of 37, Fig. 1 (top left), and partially suppressed when Lh of 55 is used, Fig. 1 (top right). However, they were preserved when Lh of 69, Fig. 1 (bottom left). The same behavior ocurred for RE signal, Fig. 1 (bottom right). Therefore, Lh of 69 was select as the window length in frequency for the remainder of analysis. The window type was selected considering the cross terms to fault components (mainlobes) ratio (XMR). Cross terms were estimated, Sicic and Boashash (2001), the main goal was to minimize the XMR criteria, eqn (17) to reduce the cross terms: XMR = Ax (t)/Am (t) (17) Different types of window functions (Rectangular, Hamming, Hanning, Blackman, Barlett and Papoluis) were tested using IR signal in Fig. 1 and the window type with the smallest XMR was selected, Table 4. Table 4. XMR for each window type Criteria XMR

3.1 Signal decomposition (Step: 1) Figure 3 exemplifies the LMD procedure with an OR signal. Vibration signal is decomposed in a series of PF, Fig. 3 (bottom left). Each PF is composed by an envelope signal and a purely FM signal computed using eqns (7) and (9), Fig. 3 (top). The effect of spectrum segmentation achieved at the end of the sifting procedure is shown in Fig. 3 (bottom). Mode mixing between PF is minimized using the proposed objective function f , eqn (15). 3.2 Selection of the optimal (Steps: 2, 3) The probability density distribution (pdf ) of original signal x(t) and each PF were compared using KLD, eqn (14). Table 5 shows t PF1 and PF2 have the smallest values of KLD; therefore, these two components are added and the result is chosen for the remainder of the analysis. Table 5. Selection of PF based on KLD. Bearing Condition

Window function types Rectangular

Hanning

Hamming

Blackman

Barlett

Papoulis

0.14192

0.03941

0.04921

0.03946

0.07417

0.06076

197

OR fault

KLD values

Selected PF

PF1

PF2

PF3

PF4

PF5

0.0075

0.2056

0.2898

0.3661

0.4958

PF1 + PF2

3.3 Time-frequency analysis (Step: 4)

3. METHODOLOGY Figure 2 shows the proposed methodology for bearing fault diagnosis.

Time-frequency analysis is carried out by SPWVD by establishing independent sliding window functions in time and frequency. Hanning windows with 117 samples (Lg) was applied in time considering the lowest value of shaft speed of the experiments (1,730 RPM), because it provides the higher values of time duration of fault components. Lh of 69 was set as the window function applied in frequency, Table 3. The whole time-frequency plane obtained by SPWVD is then divided into N time-frequency blocks of equal area with an energy magnitude of wi (i = 1, 2, 3, , N ). Each wi is normalized, eqn (18), Marie et al. (2018). Let m × m the size of each block, m is varied using values of the form 2n , with n = 1,2,...,7, i.e., TFR is divided into blocks of size: 4 × 4, 8 × 8, 16 × 16, 32 × 32, 64 × 64, 128 × 128 and 256 × 256. N ∑ pi = wi / wi (18) i=1

The fluctuation of energy in SPWVD can be analyzed based on the block size to obtain specific features from each signal, with a minimum interference cross terms. 3.4 Feature extraction method (Steps: 5, 6) TFR of signals of healthy and faulty bearings are different, representative features can be extracted from them. Shanon Entropy (SE) is included to obtain a measure of the change of energy distribution of a signal in timefrequency domain, eqn (19) (using the segmented TFR plane): N ∑ SE = − pi log2 pi (19) i=0

Fig. 2. Proposed methodology for bearing fault diagnosis. 197

Table 6 shows the feature extraction results for 7 signals which correspond to each bearing condition. Each row

2019 IFAC ICONS 198 Israel Ruiz Quinde et al. / IFAC PapersOnLine 52-11 (2019) 194–199 Belfast, Northern Ireland, August 21-23, 2019

represents SE values for each TFR segmentation. To increase the effectiveness of the method, 5 mean peaks was also extracted from each signal spectrum and finally feature vectors with SE and frequencies values are formed. SE measures the changes in energy signal while mean peaks of frequencies give information about the location of fault components. Table 6. SE feature vectors

17.68 17.96 17.46 17.52 17.27 17.93 17.01

16.02 16.31 15.82 15.85 15.62 16.27 15.33

16 × 16

TFR division 32 × 32 64 × 64

128 × 128

256 × 256

14.19 14.47 13.99 13.99 13.81 14.43 13.54

SE values 12.26 10.30 12.52 10.55 12.11 10.25 12.08 10.18 11.89 9.96 12.49 10.52 11.67 9.76

8.33 8.57 8.35 8.28 8.05 8.54 7.84

6.39 6.62 6.41 6.41 6.14 6.58 5.96

3.5 Fault Classifier (Step: 7)

Fig. 3. LMD results applied to OR signal, (q = 5).

An Artificial Neural Network (ANN) was used as fault classifier to test the proposed methodology. Drive End (DE) and Fan End (FE) acceleration signals from CWRU were classified in Table 7. They were grouped according to bearing condition and faults sizes in inches (0.007, 0.021). A set of 500 signals for each condition were generated using rectangular windows of 50% overlap, every signal has 2,048 points. Table 7. Classification of experimental data. Class 1 2 3 4

Bearing condition IR (0.007-Initial) IR (0.021-Several) OR (0.007-Initial) OR (0.021-Several)

Class 5 6 7

Bearing condition RE (0.007-Initial) RE (0.021-Several) Normal

Optimal PF

Magnitude of PF3(f)

Bearing condition 1 2 3 4 5 6 7

8×8

3000

Analyzed Signals

4×4

2500

2000

1500

1000

500

0 1

2

PF

3

0.018

Fault frequency = 107.3 Hz

0.016 0.014 0.012 0.01 0.008 0.006 0.004 0.002 0 0

200

400

600

800

1000

Frequency (Hz)

Fig. 4. Optimal PF method (left). Spectrum of PF3 (right).

4. RESULTS The proposed methodology was evaluated with experimental dataset. The maximum number of PF defined for LMD was q = 5, because it was enough to segment the whole signal spectrum. Automatic optimal PF selection was tested with signals in the training group. Figure 4(left) shows that components 1 and 2 are the most frequently selected, because their power spectra are located near to the bearing resonance band. A random signal of this group (OR, class of 4) was analyzed via FT of its envelope AM3, proving that the algorithm can still solve the fault frequency (107 Hz) in the PF3 bandwidth, Fig. 4(right). The method sometimes chooses PF3 for this reason. The envelope of the optimal PF, can be selected for further analysis by WVD to have a better illustration about fault characteristic frequencies on time-frequency domain. Figure 5 confirms the existence of peaks coincident with the OR, IR and RE fault frequencies, computed by eqns (1)-(2), when 3 signals of these bearing conditions are analyzed. The performance of the classifier can be visualized as a scattering plot, also in the confusion matrix, Fig. 6(top). 198

Fig. 5. Visual diagnosis based on proposed methodology for OR, IR and RE fault signals. The obtained performance is shown in Table 8; EMD combined with Hilbert-Huang Transform (HHT ) is included for comparison purposes. The procedure of selecting the optimal bandwidth based on KLD was applied to optimal Intrinsic Mode Functions (IMF ) generated by EMD and the feature extraction method was applied to their TFR obtained by HHT. Each result is the average of five tests. RE class achieved the lowest performance. A RE signal (shaft speed of 1,797 RPM) was analyzed. The marginal in frequency of the WVD computed for the envelope of selected PF is plotted in Fig 7. Dominant peaks coincident

2019 IFAC ICONS Israel Ruiz Quinde et al. / IFAC PapersOnLine 52-11 (2019) 194–199 Belfast, Northern Ireland, August 21-23, 2019

in their TFR. Taking advantage of WVD, feature vectors of short length were enough to represent the whole TFR. The method ends to be efficient because optimal bandwidth is analyzed instead of the whole signal spectrum. The method uses a simple feature extraction approach to reduce the computational cost using the TFR as an array of energy values instead of plotting the WVD for image processing. The average results for classification are higher than other TFR proposals such as EMD + HHT. According to the obtained results by the ANN, the proposed methodology effectively separates different bearing conditions with an average accuracy of 98.5 %.

Data classification

60

t-SNE dimension 2

40 20 0 -20 -40 -40

-20

0

20

199

REFERENCES

40

t-SNE dimension 1

Fig. 6. Features classification based on ANN Table 8. Classification results (rate of false negatives) Classification accuracy of diferent bearing condition Average Method IR(1) IR(2) OR(3) OR(4) RE(5) RE(6) Normal(7) accuracy Proposed 100% 97.90% 97.20% 98.70% 99.40% 96.00% 100% 98.50% EMD + HHT 87.70% 86.70% 87.20% 91.40% 77.30% 87.70% 97.80% 87.80%

with OR (107 Hz) and IR (162 Hz) frequencies, this could lead to misleading classification.

Fig. 7. A RE signal (other faults) 5. CONCLUSIONS A time-frequency methodology by combining LMD and WVD for bearing fault diagnosis is proposed. Fault frequency components of vibration signals were highlighted 199

Boashash, B. (2016). Advanced Time-Frequency Signal and System Analysis (Chapters 4 & 7). Academic Press, 2nd Ed, Oxford. Chen, H., Wang, J., Li, J., and Tang, B. (2017). A Texture-based Rolling Bearing Fault Diagnosis Scheme using Adaptive Optimal Kernel Time Frequency Representation and Uniform Local Binary Patterns. Measurement Science and Technology, 28. Cohen, L. (1989). Time-Frequency Distributions-A Review. Proc. of the IEEE, 77(7), 941–981. Du, W., Wang, Z., Gong, X., Wang, L., and Luo, G. (2016). Optimum IMFs Selection based Envelope Analysis of Bearing Fault Diagnosis in Plunger Pump. Shock and Vib., 2016, 1–8. Hyers, R., McGowana, J., K. Sullivan, J.M., and Syrett, B. (2006). Condition Monitoring and Prognosis of Utility Scale Wind Turbines. Energy Materials, 1(3), 187–203. Ibrahim, G.R. and Albarbar, A. (2011). Comparison between WignerVille Distribution- and Empirical Mode Decomposition Vibration-based Techniques for Helical Gearbox Monitoring. Proc. of the Institution of Mech. Eng., Part C: J. of Mech. Eng. Science, 225(8), 1833–46. Liu, Z., Zuo, M., Jin, Y., Pan, D., and Qin, Y. (2017). Improved Local Mean Decomposition for Modulation Information Mining and its Application to Machinery Fault Diagnosis. J. of Sound and Vib., 397, 266 – 281. Ma, Z., Ruan, W., Chen, M., and Li, X. (2018). An Improved Time-Frequency Analysis Method for Instantaneous Frequency Estimation of Rolling Bearing. Shock and Vib, (8710190). Marie, T., Han, D., An, B., and Li, J. (2018). A Research on FiberOptic Vibration Pattern Recognition based on Time-Frequency Characteristics. Advances in Mech. Eng., 10(12). ´ (2014). Dissipation and Kullback-Leibler Divergence, 37– Rold´ an, E. 59. Springer Int. Publishing. Si, L., Wang, Z., Tan, C., and Liu, X. (2017). Vibration-based Signal Analysis for Shearer Cutting Status Recognition based on Local Mean Decomposition and Fuzzy C-means Clustering. Applied Sciences, 7(2), 164. Sicic, V. and Boashash, B. (2001). Parameter Selection for Optimizing Time-Frequency Distributions and Measurements of TimeFrequency Characteristics of Non-stationary Signals. In IEEE Int. Conf. on Acoustics, Speech, and Signal Proc, 3557–3560. Sinha, J. (2015). Analysis, Instruments, and Signal Processing. CRC Press, 1st Ed. Smith, J. (2005). The Local Mean Decomposition and its Application to EEG Perception Data. J. of the Royal Society Interface, 2(5), 443–454. Smith, W. and Randall, R. (2015). Rolling Element Bearing Diagnostics using the Case Western Reserve University Data: A Benchmark Study. Mech. Syst. and Signal Proc., 64, 100–131. Xie, P., Yang, Y., Jiang, G., Du, Y., and Li, X. (2012). A New Fault Detection and Diagnosis Method based on Wigner-Ville Spectrum Entropy for the Rolling Bearing. Applied Mech. and Materials, 197, 346–350.