Residual useful life prediction of large-size low-speed slewing bearings – a data driven method

2.1. Ensemble empirical mode decomposition

As a self-adaptive data-driven method, Empirical Mode Decomposition (EMD) is very suitable for unstable signal processing. However, the non-linear components in non-stationary signals may result in different frequency scales projected on one IMF, and this phenomenon is also called mode mixing. To solve this problem, Wu and Huang [14] proposed the Ensemble Empirical Mode Decomposition method in 2009, which put white noise in EMD procedures to restrain mode mixing effect. For signal $a\left(t\right)$, the procedures of EEMD are as follows:

(I) Put a randomly generated white noise $s\left(t\right)$ with magnitude $z$ to raw signal $a\left(t\right)$ to generate a new data set $a_{new}\left(t\right)$:

(II) Use EMD to decompose $a_{new}\left(t\right)$ into $H$ IMFs, the $h$-th ($h\in \left[1,H\right]$) IMF is represented as $f_h\left(t\right)$.

(III) Repeat steps (I) and (II) $I$ times with different white noise series to obtain an ensemble of IMFs, the $h$th ($h\in \left[1,H\right]$) IMF after $i$th ($i\le I$) EMD is $f_h^i\left(t\right)$.

(IV) Obtain the $h$th IMF of EEMD by calculating the ensemble means of corresponding IMFs:

It should be noted that, the magnitude of white noise $z$ and the repeated times $I$ need to be determined during EEMD process. Chen et al. [15] pointed out that under a same noise level, $z$ is not a critical factor to calculation results, and the increase of $I$ can only slow down the calculation. Therefore, in this paper, $z$ is 0.2 and $I$ is 100.

2.2. Principal component analysis

PCA de-correlates variables through extraction of linear relationships between variables, and uses several unrelated components to represent all variables. For matrix $\mathbf{X}\left(m\times n\right)$, where m represents the number of samples and n stands for the number of variables, PCA decomposes it into principal component space and residual subspace:

where ${\mathbf{t}}_j$ ($m\times$1) is the $j$th principal component (PC), ${\mathbf{p}}_j$ ($n\times$1) stands for the $j$th load vector, and $\mathbf{E}$ represents the residual error. $\mathbf{X}$ can also be presented as:

where $\mathbf{T}\left(m\times r\right)$ is the principal component matrix, and $\mathbf{P}\left(n\times r\right)$ is the load matrix.

Therefore, the projection on residual space is:

To implement condition monitoring, most PCA based schemes use two typical statistical indices: Hotelling’s $T^2$ and squared prediction error (SPE), also known as $Q$ statistic. According to literature [16], SPE can better demonstrate the differences between new datasets and the original dataset. Meanwhile, it would be seen as abnormality if the overall trend of the vector $\mathbf{S}\mathbf{P}\mathbf{E}$ is greater than its threshold $Q_{\alpha }$, and they can be calculated by Eqs. (6), (7):

where $c_{\alpha }$ is the standard normal deviate corresponding to the upper ($1-\alpha$) percentile, ${\theta }_d={\sum }_{o=r+1}^n{\lambda }_o^d$ for $d=$1, 2, 3, $h_0=1-2{\theta }_1{\theta }_3/3{\theta }_2^2$, and ${\lambda }_o$ is the $o$th eigenvalue of the variance–covariance matrix of input data $\mathbf{X}$.

2.3. Least square – support vector machine

Proposed by Vapnik [17], SVM was developed to solve pattern recognition problems based on the structural risk minimization theory. With SVM, different data sets are firstly separated by hyperplanes as much as possible, and then an optimal hyperplane will be found to maximize the distances between all the hyperplanes. The normal vectors of the optimal hyperplane are called support vectors. In order to improve machine learning speed, LS-SVM [18] was then proposed to transfer a quadratic programming problem into a one-dimensional linear problem by replacing the insensitive loss function used in SVM with a squared error loss function.

For an $N$ sample training set {${\mathbf{x}}_k$, $y_k$}, $k=$1, 2,…, $N$, where ${\mathbf{x}}_k$ (${\mathbf{x}}_k\in {\mathbf{R}}^n$) is the $k$th input vector while $y_k$ ($y_k\in \mathbf{R}$) represents the output. The linear regression function of LS-SVM is:

where $\phi \left({\mathbf{x}}_k\right)$ stands for the mapping function in kernel space, $\mathbf{\omega }$ represents the normal vector of the hyperplane and $b$ is a constant. This linear problem can also be described as:

where $e$ is the training error, $\gamma$ is defined as a regularization parameter to adjust the penalty for training errors. By introducing Lagrange multipliers ${\mathbf{\alpha }}_k$ and optimal constraints, the LS-SVM prediction model can be described as:

where $K\left(\mathbf{x},{\mathbf{x}}_k\right)$ is the kernel function, which maps the sample space into a high-dimensional space. It should be noted that commonly used kernel functions are linear kernel function, polynomial kernel function and radial basic function (RBF).

4.1. Experimental setup and results.

To evaluate the performance of the proposed method, a home-made test platform is developed to implement an experiment. As shown in Fig. 3, G1, G2 and G3 are hydraulic rams to apply axial force, radial force and turnover movement to a tested slewing bearing, and G4 is the hydraulic motor to drive the slewing bearing. More information of the test platform can be found in Table 1.

Fig. 3Home-made test platform

As shown in Fig. 4 and Fig. 5, the type of the main tested slewing bearing is QNA-730-22, of which the inner ring is the rotatable ring and the outer ring is the fixed ring. The properties of the slewing bearing are shown in Table 2.

Fig. 4Improved test platform for validation of the proposed method

Given that the tested slewing bearing is inner geared, some improvements have been made to the test platform. With the frocks shown in Fig. 4, the inner ring (with gears) of the main tested slewing bearing is connected with the outer ring of the accompanied slewing bearing, and the outer (fixed) ring of the main tested slewing bearing is fixed with the upper flange. After which, as shown in Fig. 5, extreme loads are applied to the main tested slewing bearing and the experiment starts at a constant speed of 4 r/min.

Fig. 5Applied loads and the distribution of the four accelerometers

To implement condition monitoring, as shown in Fig. 5, four 495 mv/g Kistler-8395A010 accelerometers are mounted along the radial direction of the slewing bearing for vibration measurement, and the vibration signals are collected by NI PXIe-4492 at a sampling rate of 2 kHz. Besides, to assist monitoring the slewing bearing condition, the grease temperature, driven torque and rotational speed are also measured by NI PXI-6238 at a sampling rate of 10 Hz. What’s more, the main data acquisition (DAQ) system, NI PXIe-1062Q collects all the data from both PXIe-4492 and PXI-6238 through a LabVIEW interface. The schematic of the DAQ system is shown in Fig. 6.

Fig. 6DAQ system schematic diagram

The test stops when the slewing bearing gets stuck and comes to a critical failure after 11 days, and the final amount of the turning laps of the slewing bearing is 6.33×10⁴. Fig. 7(a) and Fig. 7(b) show the trends of the grease temperature and the driven torque respectively. It should be specially explained that, to study the slewing bearing condition when its fatigue life reaches the Chinese standard [19] of 30000 r, we suspended the test and maintained the tested slewing bearing on the 7th day, and that’s the reason that both the grease temperature and the driven torque decrease after the 7th day. Besides, Fig. 7 also illustrates that both the grease temperature and the driven torque increase dramatically after the 9th day until the slewing bearing comes to a critical failure. In addition, Fig. 8(a) shows the slip zone occurred in the outer ring raceway of the tested slewing bearing on the 7th day, while Fig. 8(b) and Fig. 8(c) show the damaged outer ring and the cracked balls respectively when the slewing bearing gets stuck.

Fig. 7Trends of both the grease temperature and the driven torque during whole test: a) trend of the grease temperature and b) trend of the driven torque

a)

b)

4.2. RUL prediction and discussions

The lowest fault frequency of the tested slewing bearing is 0.34 Hz according to the commonly used bearing fault frequency calculation equations [20], and 0.34 Hz will have one full cycle period in at least 2.9 seconds [8]. Besides, during the whole test, vibration signals acquired at late night may reduce the environmental interferences as much as possible. Therefore, as shown in Fig. 9, we choose 3 seconds late-night acquired vibration data (6144 samples) of each day for signal processing and RUL prediction.

Fig. 8Damaged slewing bearing parts

a) Slip zone

b) Damaged outer ring

c) Cracked balls

Fig. 9Measured life-cycle vibration signals

It can be observed from Fig. 9 that the details of all four signals are submerged under white noise and the energy of the periodical impulse is very low. For simplicity, only the de-noising procedures of ${\mathbf{a}}_1$ will be presented below and the de-noising effectiveness will be compared with another EEMD-PCA based de-noising method [6]. Firstly, signal ${\mathbf{a}}_1$ is divided into 11 segments according to the experimental days, and then all the signal segments are decomposed by EEMD respectively. The IMFs of the 5th signal segment are shown in Fig. 10 as an example.

Afterwards, the first IMF of both the 5th signal segment and the 1st signal segment are divided into matrix ${\mathbf{C}}_{11}$ and ${\mathbf{C}}_{15}$ respectively, and a PCA model is built on ${\mathbf{C}}_{11}$. Then, by projecting ${\mathbf{C}}_{15}$ onto the PCA model, $ER{R}_{15}$ can be obtained. Accordingly, the ERRs of all 12 IMFs of the 5th signal segment are calculated.

By repeating the above steps, the ERRs of all 13 IMFs of the remaining signal segments can also be calculated. Fig. 11 visualizes the waterfall profile of the ERR amplitude vs. experimental days (signal segments) of all 13 IMFs, in which the red dot indicates that the curve is greater than zero for the first time. As mentioned above, the amplitude of the ERR represents the degree of abnormality in the corresponding IMF, while the abnormality of the IMF is the result of some fault occurred in the tested slewing bearing. Therefore, it can be observed from Fig. 11 that some incipient abnormalities occur in the 8th-12th IMFs from the 2nd day to the 4th day, which indicates the run-in period of the slewing bearing. Two days later, abnormalities firstly occur in the 5th-7th IMFs and keep being more severe, and the ERR amplitudes of the 8th-12th IMFs keep increasing at the same time. This phenomenon stands for the performance degradation process of the slewing bearing, and different types of slewing bearing faults are projected onto different IMFs. After the 9th day, the ERR amplitudes of the 6th-9th IMFs increase rapidly, which demonstrates the critical failure of the slewing bearing. It also should be noted that, in accordance with our previous statement, the ERR amplitudes of the first four IMFs keep low and change little throughout the experiment because of white noise, which will result in neglect of these IMFs in the selection process.

Fig. 10IMF 1-13 of the 5th signal segment of a1

a) IMF 1-6 of the 5th signal segment

b) IMF 7-13 of the 5th signal segment

Fig. 11Trend of the ERR of each IMF in 11 days

Fig. 12Profile of the ERR amplitude vs. all 13 IMFs of each day

Nevertheless, Fig. 11 only exhibits the trends of the ERRs of each IMF throughout the experiment. It is still not clear enough for us to determine which IMFs to choose to reconstruct signal ${\mathbf{a}}_1$. Hence, take another look at Fig. 11, Fig. 12 demonstrates the waterfall profile of the ERR amplitude vs. all 13 IMFs of each day. It should be noted that the ERRs of each day are mapped to [–1, 1] respectively to clearly demonstrate which IMFs to be selected. For $\beta =$0.9, the selected IMFs are listed in Table 3, and the selected times of each IMF are shown in Fig. 13.

Fig. 13Selected times of each IMF

Therefore, the 5th, 6th, 8th, 11th and 12th IMFs are selected to reconstruct the life-cycle signal ${\mathbf{a}}_1$ to achieve de-noising, and all the four groups of de-noised vibration signals are shown in Fig. 14.

Fig. 14De-noised life-cycle vibration signals by the proposed method

As shown in Fig. 14, the life-cycle signals are divided into three stages according to the trends of the vibration signals throughout the experiment. The amplitudes of the signals are low and stable at stage I, and then increase gradually at the middle stage. Before the end of the experiment, the amplitudes of signals rise dramatically and some periodical impulses can be found at stage III.

Meanwhile, for comparison, another EEMD-PCA based de-noising method [6] (we define it as direct EEMD-PCA) is utilized to de-noise the four groups of vibration signals. The de-noised signals and the contribution of each IMF (principal component) of ${\mathbf{a}}_1$ are shown in Fig. 15(a) and Fig. 15(b) respectively. It can be seen from Fig. 15(b) that the contribution of the first five IMFs is very significant, and this will result in improper selection of the first few IMFs which contain large amount of white noise. Consequently, this method de-noises the raw signals to a certain extent, but the de-noised signals in Fig. 15(a) are still submerged by noise. Through the comparison, we can find out that, by selecting the most representative IMFs from life-cycle signals, the proposed de-noising method has a better de-noising effect on such low signal-to-noise vibration signals than direct EEMD-PCA based methods.

Fig. 15Direct EEMD-PCA de-noising

a) De-noised signals by direct EEMD-PCA

b) Contribution of each IMF (PC) of ${\mathbf{a}}_1$

After de-noising, two time-domain features: root mean square (RMS) values and impulse factors of the four groups of vibration signals are calculated and shown in Fig. 16(a) and Fig. 16(b) respectively. It can be seen from Fig. 16(a) that the overall trend of the RMS for each vibration signal is in accordance with a general performance degradation trend of equipment, while the details differences are caused by some random factors such as geometric errors, processes and loads. Fig. 16(b) illustrates the potential faults development in different positions of the slewing bearing. On one hand, since the accelerometer ${\mathbf{a}}_1$ and ${\mathbf{a}}_3$ are mounted near the load points, their peak values of the impulse factors are quite high, which indicates severe damages of the slewing bearing around sensor ${\mathbf{a}}_1$ and ${\mathbf{a}}_3$. On the other hand, the peak values of the impact factor of ${\mathbf{a}}_1$ appear more and earlier than ${\mathbf{a}}_3$, which demonstrates that the slewing bearing raceway around ${\mathbf{a}}_1$ damages earlier and more severe because of the soft area. It is worthy to be found that all the phenomena observed from Fig. 16 are consistent with the conclusions drawn in our earlier work [21]. In a word, compared with one dimensional signal, multi-dimensional vibration signals can explain the performance degradation trend of the slewing bearing more thoroughly and precisely.

Fig. 16Two time-domain features

a) RMSs of the four de-noised signals

b) Impulse factors of the four de-noised signals

Nevertheless, the number of time-domain, frequency-domain and time-frequency-domain features of the four groups of vibration signals can reach up to 120 [10], which brings great challenges to RUL prediction modeling. To avoid dimension disasters, C-SPE of the four groups of signals and its RMS are calculated and shown in Fig. 17.

Fig. 17C-SPE and its RMS

a) C-SPE of the four de-noised signals

b) RMS of the C-SPE

It can be seen from Fig. 17 that the trend of the RMS of C-SPE is exactly in accordance with the performance degradation process of the slewing bearing observed above. The first three days reflect the run-in period, while the normal service period is demonstrated from the 3rd day to the 7th day. After maintenance, the condition of the slewing bearing improves a little on the 8th day, and soon after that, the condition becomes worse dramatically since the 9th day until a critical failure occurs. Therefore, C-SPE can be an effective feature to reduce the dimensions of the vibration signals and meanwhile accurately assess the performance degradation process as much as possible.

Given that the C-SPE in frequency domain has no practical significance, only ten dimensionless time-domain features of the C-SPE are calculated and shown in Fig. 18(a), which will be used as the inputs in RUL prediction modeling. At the same time, the same ten features of the four groups of de-noised signals are also calculated for comparison. For simplicity, only the features of ${\mathbf{a}}_1$ are shown in Fig. 18(b).

Fig. 18Ten time-domain features of both C-SPE and a1

a) Features of C-SPE

b) Features of ${\mathbf{a}}_1$

Afterwards, to establish an RUL prediction model (C-SPE based model) with LS-SVM, ten time-domain features of the C-SPE and the residual useful life (expressed in revolutions) are introduced as the inputs and the outputs respectively. It should be noted that the kernel function used in the LS-SVM is the RBF, and the optimal regularization parameter $\gamma$ and the kernel bandwidth ${\sigma }^2$ in the RBF are determined by cross validation. In addition, another RUL prediction model (signal based model) is built based on the 40 time-domain features of the four groups of signals. Both of the established models are shown in Fig. 19.

Fig. 19LS-SVM based RUL prediction models

a) C-SPE based RUL prediction model: $\gamma =$ 103.22, ${\sigma }^2=$ 0.71, Time = 10.215 s, Accuracy = 99.40 %

b) C-SPE based RUL prediction model: $\gamma =$ 75.14, ${\sigma }^2=$ 0.93, Time = 68.466 s, Accuracy = 94.30 %

It can be observed from Fig. 19 that the C-SPE based RUL prediction model has a shorter calculation time and a higher accuracy than the RUL prediction model built based on the features of signals. In fact, the C-SPE itself is an indication of the differences of signals in different time periods, therefore it can be seen from Fig. 17 and Fig. 18(a) that the C-SPE and most of its time-domain features keep increasing during the entire life cycle. On the contrary, as shown in Fig. 18(b), there are some fluctuations in many of the time-domain features of signals. Besides, the number of the input vectors of the C-SPE based model is only 10 while the number of the input vectors of the signal based model is up to 40. As a result, the signal based model is less accurate and more time consuming. It should be noted that the signal based model may achieve a higher performance if the input vectors are carefully selected from the 40 features. However, the selection procedure itself requires human involvement and thus is still inefficient.