Fault feature extraction for rolling element bearings based on multi-scale morphological filter and frequency-weighted energy operator

2.1. Basic operator of morphology filter

Mathematical morphology was first proposed by Matheron and Serra as an image processing tool [22], which can simplify the amount of data in the images and removes irrelevant information while retaining important features. Then, in 1987, morphology filter was used as a nonlinear data processing method for 1-D time series [23].

There are four basic operators when using morphology filter, including, erosion, dilation, opening and closing, meanwhile the opening and closing operators are composed of erosion and dilation operators. Let $f\left(n\right)$ be the 1-D discrete signal that is the function over a domain $F=\left\{\mathrm{0,1},2,\dots ,N-1\right\}$. Let $g\left(n\right)$ be the SE and the discrete function over a domain $G=\left\{\mathrm{0,1},2,\dots ,M-1\right\}$. Both $N$ and $M$ are integers, $N>M$. Thus, four basic operators can be defined as:

Dilation:

Erosion:

Opening:

Closing:

where in Eq. (1)-(4), $\mathrm{\Theta }$, $\oplus$,$\circ$ and $•$ represent erosion, dilation, opening and closing operator respectively while max and min represent maximum and minimum values. In theory, the erosion operator can smooth negative impulses and reduce positive ones while the dilation operator has an opposite effect. The closing operator can be used to preserve positive impulses and eliminate negative ones. In contrast, the opening operator is able to extract negative impulses.

2.2. The definition of ACDIF

ACDIF [16] is constructed by combining four basic operators. At first, a combination of dilation operator and closing operator is defined called Dilation-Closing operator and Closing-Dilation operator:

Then, a combination of erosion operator and opening operator is defined called Erosion-Opening operator and Opening-Erosion operator:

Thus, the ACDIF can be defined as:

where the $F_{CD-EO}$ and $F_{DC-EO}$can be defined as:

By using ACDIF, it can not only suppress noise, but also improve the impact characteristics of signals

2.3. Multi-scale morphology filter

The main idea of multi-scale morphology filter is to analyze the signal at different scales by using different SEs. The denoising ability of large-scale SE is relatively strong, and it can reflect the contour features of signals, but some detail information may be neglected. A small-scale SE can preserve the detail features of the signal, but the ability of denoising is relatively weak. Thus, multi-scale morphology filter seems to be more proper than single-scale morphology filter when applied to complicated bearing fault signals.

Let $\lambda$ be the scale, the SE at scale $\lambda$ is defined as:

Then the multi-scale ACDIF at scale $\lambda$ can be expressed as:

The multi-scale filtered results are usually obtained by averaging the filtered results at all scales or the sum of weighted filtered results at all scales. In this paper, the weighted results are calculated which can be expressed as:

where ${\lambda }_1$,${\lambda }_{max}$ denotes the minimum and maximum SE scale, $w_i$ represents weight coefficient of the filtered result at scale $i$ and $N_{\lambda }$ denotes the number of all selected SE scales.

2.4. The process of adaptive selection of weight coefficients

SE is also an important operator for morphology filter, which has a great influence on the filtered results. SE is mainly distinguished in shape, amplitude and length, commonly used shapes of SEs include cosine, semicircle, triangle, flat and their combination. Due to the focus of many studies are the optimal selection of SE scale, to improve the computational efficiency, flat SE with a magnitude of zero has been widely used. Recently, a new SE called double-dot is proposed by Chen [24], which is proved to be better than flat SE in fault feature extraction. Therefore, in this study, the double-dot SE is used instead of the commonly used flat SE. For double-dot SE, when the scale $\lambda =$ 1, the SE is {1,0,1}; when $\lambda =$ 2, the SE is {1,0,0,1}; when $\lambda =$ 3, the SE is {1,0,0,0,1}, and so on.

To adaptively select the weight coefficients of SE scale, a newly proposed H-PSO-SCAC method is used. The H-PSO-SCAC is proposed to solve the problem existing in the traditional particle swarm optimization (PSO) method, such as premature convergence and easily trapped in the local optimum solution, which has been proved to achieve good results [25].

Assuming that the size of the population is $N$ and the dimension of the search space is $D$, the position of the $i$th particle can be represented by a vector as $X_i^d=$ ($x_{i1}$, $x_{i2}$, …, $x_{iD}$). The velocity of a particle $i$ is denoted by vector $V_i^d=$ ($v_{i1}$, $v_{i2}$, …, $v_{iD}$). The optimal position vector of $i$th particle is represented by $Pbest=$ ($Pbes{t}_i^1$, $Pbes{t}_i^2$, …, $Pbes{t}_i^D$), called personal best position and the best position of the population is denoted as $gbest=$ ($gbes{t}^1$, $gbes{t}^2$, …, $gbes{t}_i^D$). The position $X_i$ and velocity $V_i$ are initialized randomly and updates to the D-dimension of the $i$ particle. The process of the standard PSO is expressed as:

where $c_1$ and $c_2$ are the cognitive component and the social component respectively. $\omega$ denotes the inertia weight and $r_1$, $r_2$ are random number between 0 to 1.

In H-PSO-SCAC, the modification is as follows:

First, $c_1$ and $c_2$ are defined as:

where $M_j$ represents the $j$th iteration and $M_{max}$ denotes the maximum number of iteration.

Second, $\omega$ is defined as:

where $x_k\in \left(\mathrm{0,1}\right)$, $c$ denotes a random number between 0 to 4, and $k$ is the current iteration number.

Third, the population initialization method is changed which can increase the opportunities of reaching the global optimal solution by fifty percent [25]:

The initial population is randomly chosen as $p$($M=$ 0)$=\left\{x_{ij}\right\}$, $i=$1, 2, …, $N$ and $j=$1, 2, …, $D$, then, the reverse population $p\text{'}\left(M=0\right)=\left\{x_{ij}^{\text{'}}\right\}$ is calculated as:

where $x_{max,j}$ and $x_{min,j}$ are population position $x_i$ at $j$ dimension’s max value and min value respectively. Finally, the smaller half of $x_{ij}$ and $x_{ij}^{\text{'}}$ are selected and combined as the initial position.

Forth, the search form described in Eq. (16) is modified as:

where $w_{ij}$,$w_{ij}^{\text{'}}$ and $\psi$ are defined as:

where $u$ is the mean fitness value in the first iteration, $iter$ is the current iteration and $f\left(j\right)$ is the fitness of $j$th particle.

As the correlation kurtosis (CK) can evaluate the impact characteristics of a signal, which means the larger CK value corresponds to the more obvious impact components in the signal. Hence, the maximum CK is chosen as the fitness function in the optimization process which can be expressed as:

where $x_n$ represents the discrete signal, $N$ denotes the length of the signal, $T$ denotes the period of interest, $M$ represents the order of shift. CK overcomes the disadvantage that kurtosis cannot reflect specific signal characteristics. When CK is applied to fault diagnosis of rolling element bearings, the $T$ is set to be failure cycle. As a result, when a localized defect occurs, the impulse signal corresponds to the bearing fault cycle has the maximum CK, while the CK of other signals is relatively small.

5.1. Simulation signal

To verify the effectiveness of the proposed fault feature extraction method in this paper, a simulation signal is established and the model of simulation signal is:

In Eq. (31), the value of the rotating frequency $f_r$ is 42 Hz, $A_i$ is the amplitude modulation which cycle is 1/$f_r$, and $f_r=$ 42 Hz. $f_n=$ 3200 is the natural frequency of the system and $r=$ 0.05 is the damping coefficient, $t_i=$ 0.01 $T$ represents the delay caused by rolling element slippage in $i$th cycle. The repetition period of the impulses is 1/185 s which means the fault characteristic frequency is 185 Hz. The harmonic $B\left(t\right)$ is applied to simulate the interference components existing in the signal while the $f_1$ and $f_2$ denotes 50 Hz and 100 Hz, respectively. The signal $n\left(t\right)$ is used to simulate the random noise which can be applied by MATLAB using $randn$(1, $n$), while $n$ equal to the length of the signal. The sampling frequency in simulation is 32768Hz while the length of the signal is 1 second.

Fig. 2 describes the simulation signal $x\left(t\right)$ and its FFT spectrum. It can be seen from Fig. 2(a) that periodic impulses are submerged by strong background noise which can hardly be observed. Due to the influence of noise and interference components 50 Hz and 100 Hz, it is difficult to find the fault characteristic frequency and its harmonics in Fig. 2(b). This means some other method is needed for the fault feature extraction of simulation signal.

Fig. 2Simulation signal: a) time domain waveform, b) FFT spectrum of a)

a)

b)

5.2. Simulation analysis

Then the proposed method is applied to the simulation signal. As is shown in Ref. [29] that the optimal SE scale is closely related to both the sampling frequency of the signal and the points in a fault cycle. To get a better preference of multi-scale morphological filter, the maximum SE scale should between 0.2 to 0.3 times of the points in a fault cycle and 0.2 times is chosen in this study. In the simulation, the sampling frequency is 32768 while the fault characteristic frequency is 185 Hz, so the points in a fault cycle is about 177 and the maximum SE scale is 35. What’s more, as the adjacent SE scale has the similar filtering result [16], the SE scales are set to:

Hence, 12 SEs are used in the multi-scale ACDIF. Fig. 3(a) shows the convergence curve, the optimal fitness value is 3.03×10^-14 after 12th iteration while the optimal weight coefficients are:

After the acquiring the optimal weight coefficients, the filtered result of multi-scale ACDIF is obtained which is shown in Fig. 3(b). Fig. 3(c) shows the demodulation result of the filtered signal by FWEO where the rotating frequency 42 Hz and its harmonic 84 Hz, the fault characteristic frequency 185 with two of its harmonics 370 Hz and 555 Hz are prominent. Moreover, the sidebands modulated by rotating frequency are able to be observed. As a result, the extracted fault features indicate the existence of a defect on the inner race.

Fig. 3Results of the proposed method in simulation: a) the convergence curve, b) the final filtered result, c) FWEO spectrum of b)

a)

b)

c)

Fig. 4Results of optimal single-scale ACDIF in simulation: a) relationship between SE scales and TEK, b) FWEO spectrum of optimal filtered result

a)

b)

By comparison, optimal SE scale is selected using the method proposed in Ref. [16] by teager energy kurtosis (TKE), then the optimal single-scale ACDIF will still be demodulated by FWEO and the results are shown in Fig. 4. Fig. 4(a) shows the relationship between the SE scales 1-35 and the TEK. The SE scale of 35 corresponds to the maximum TEK will then be used in optimal scale ACDIF. Fig. 4(b) displays the FWEO spectrum of the optimal filtered result where the fault characteristic frequency is prominent but its harmonics can hardly observed which is no better than the result in Fig. 3(c). Hence, the proposed method may be more suitable for the fault feature extraction in this simulation.

6.1. Experimental data

In this study, a bearing test rig is used to further evaluate the effectiveness of the proposed method for fault feature extraction of rolling element bearing. The test rig showed in Fig. 5(a) consists several main components which are bearing support structure, main shaft, experimental bearing, lubricating oil system, servo-driven motor, radial loading device, axial loading device and control system.

Fig. 5Experimental setup: a) the test rig, b) bearing defect on the inner race

a)

b)

The experimental bearing type is 6010, the outer and inner diameter of the bearing are 80 mm and 50 mm respectively, thus the pitch diameter is $D=$ 65 mm. The diameter and number of rolling elements are $d=$ 9 mm and $Z=$ 13. The contact angle is $\alpha =$ 0°.

The experimental signals are collected by the B&K vibration test system where the models of the acceleration sensor, the signal conditioning and acquisition module are B&K4354B-004 and B&K3053-B-120, respectively. Two conditions are simulated in the experiment when collecting the vibration signals: (1) normal bearing condition (signal 1) (2) bearing with a defeat on the inner race generated by laser wire-electrode cutting (shown in Fig. 5(b), signal 2).

To make the fault feature extraction more difficult, the measuring point is chosen on the outside of bearing support structure (12 o’clock, shown in Fig. 5(a)), the driving motor is set to rotate at a speed of $f_r=$ 3000 r/min. The sampling frequency of the data acquisition device is set to 65536 Hz. The fault characteristic frequencies for different localized defects can be calculated by the following equations:

where $f_i$, $f_o$ and $f_b$ represent the ball pass frequency of the inner race (BPFI), the ball pass frequency of the outer race (BPFO) and the ball spin frequency (BSF), respectively. Therefore, the theoretical BPFI, BPFO and BSF are 370 Hz, 280 Hz and 177.09 Hz. The time domain waveform of the collected vibration signals (all the experimental signals are collected in the radial direction) and their FFT spectrum are shown in Fig. 6.

Fig. 6(a) and (b) shows the time domain waveform of signal 1 and its FFT spectrum, respectively. Fig. 6(c) shows the time domain waveform of signal 2 with the length of 32768 points, where periodic impulses are submerged by strong background noise and can hardly recognized. The FFT spectrum of signal 2 is shown in Fig. 6(d) which displays a complicated frequency map containing many frequency components. The BPFI cannot be detected easily, hence, judging whether the bearing is damaged or not by the FFT spectrum of the original signal is difficult.

Fig. 6Experimental bearing vibration signal: a) time domain waveform of signal 1, b) FFT spectrum of a), c) time domain waveform of signal 2, d) FFT spectrum of c)

a)

b)

c)

d)

6.2. Fault feature extraction based on proposed method

Then the proposed method is applied to the experimental signals. For bearings with inner race failure, as the sampling frequency and BPFI are 65536 Hz and 370 Hz respectively in the experiment, so the points in a fault cycle is about 177 and the maximum SE scale is 35 (0.2 times of the points in a fault cycle). To avoid selecting the adjacent SE scales, the SE scales are set to:

which is similar to the simulation. So, 12 SEs are used in the multi-scale ACDIF and the H-PSO-SCAC is used to select the optimal weight coefficients. Similarly, for bearing with outer race failure, fifteen SE scales are chosen which are:

for bearing with rolling element failure, totally 25 SE scales are selected which are:

For signal 1, as the selection SE scales are closely related to the fault characteristic frequency, for a bearing with no defects, the proper SE scales may be hard to choose. So, all three range of SE scales mentioned above are used, results are shown in Fig. 7, where the filtered results are shown in (a), (c) and (e), corresponding FWEO spectrum are shown in (b), (d) and (f). As can be seen from Fig. 7(b), (d) and (f), only the spectral lines of 50 Hz and 100 Hz are prominent which represent the rotating frequency and its harmonic (mainly caused by the misalignment existing in the test rig). The BPFI, BPFO and BSF are unable to be observed which is consistent with the actual working condition of the rolling element bearing.

Fig. 7Processing results of signal1 using the proposed method: a) the filtered result using e¨1, b) FWEO spectrum of a), c) the filtered result using e¨2, d) FWEO spectrum of c), e) the filtered result using e¨3, f) FWEO spectrum of e)

a)

b)

c)

d)

e)

f)

Subsequently, signal 2 is analysed. As depicts in Fig. 8(a), based on the convergence curve, the optimal fitness value is achieved after 8th iteration which value is 3.03×10^-14. Thus, the corresponding best position of the population at 8th iteration is considered as the optimal weight coefficients which is: [0.815, 0.802, 0.720, 0.185, 0.042, 0.282, 0.052, 0.013, 0.125, 0.022, 0.376, 0.370]$.$

Then the final filtered result of multi-scale ACDIF is obtained which is shown in Fig. 8(b). The FWEO spectrum of the final filtered signal is shown in Fig. 8(c) where the rotating frequency 50 Hz and its harmonic 100 Hz, the frequency component 368 Hz (close to the theoretical value of BPFI, the difference may cause by size error of bearing or the slippage of rolling elements) with two of its harmonics 736 Hz and 1104 Hz are prominent. Moreover, the sidebands modulated by rotating frequency around the BPFI and its harmonics are able to be observed. All the features above clearly indicates the existence of a defect on the inner race.

Fig. 8Processing results of signal 2 using the proposed method: a) the convergence curve, b) the final filtered result, c) FWEO spectrum of b)

a)

b)

c)

Fig. 9Processing results of signal 2 using optimal signal-scale ACDIF: a) relationship between chosen SE scales and TEK, b) FWEO spectrum using the SE scale 35, c) FWEO spectrum using the SE scale 6

a)

b)

c)

Two methods are used to make comparisons when analysing signal 2, one is similar to the method used in simulation, where the TEK of filtered results at scales 1-35 are calculated and the results are shown in Fig. 9(a). Fig. 9(a) depicts when using the SE scale of 35, the maximum TEK can be achieved. So, the optimal-scale ACDIF is employed by using the SE length of 35, then the FWEO is applied to the filtered signal to extract the fault features of bearing. Fig. 9(b) shows prominent rotating frequency 50 Hz and its harmonic 100 Hz, the BPFI 368 Hz along with one of its harmonic 736 Hz, but the frequency line corresponds to 1104 Hz are unable to observed. Moreover, as the TEK at scale 6 is prominent compared with the adjacent scales, the FWEO spectrum using the SE scale 6 is also calculated and the result is shown in Fig. 9(c) where the BPFI and its harmonics can be clearly observed, which surely indicates the existence of a fault on the inner race. Hence, when the method in Ref. [16] is applied to the experimental signal in this study, using the SE scale corresponds to the maximum TEK may not achieve the optimal result.

Another method used for comparison is the Fast Kurtogram proposed in [29] and the results are shown in Fig. 10. In Fig. 10(a), the central frequency and bandwidth of the optimal analysis frequency band can be determined as 30208 Hz and 1024 Hz whose square envelope spectrum is shown in Fig. 10(b). In Fig. 10(b), the BPFI 368 Hz and its harmonic 736 Hz are able to be observed but not prominent while the interference frequency components are obvious. Thus, the Fast Kurtogram cannot achieve the ideal results when analysing the signal 2, which further verifies that the proposed method is effective in impulse components extracting from the vibration signal with heavy background noise.

Fig. 10Analysing results of signal 2 using the fast Kurtogram: a) figure of fast Kurtogram, b) square envelope spectrum of signal within optimal analysis frequency band

a)

b)