Health condition monitoring of bearings based on multifractal spectrum feature with modified empirical mode decomposition-multifractal detrended fluctuation analysis

Abstract

Multifractal detrended fluctuation analysis (MFDFA) is proved to be a powerful tool for fault diagnosis of rotating machinery due to its ability to reveal multifractal structures hidden in nonstationary and nonlinear vibration signals. To overcome the discontinuity of the fitting scale-dependent trend and the poor adaptability of this algorithm, Empirical Mode Decomposition-Multifractal Detrended Fluctuation Analysis (EMD-MFDFA) is introduced. However, EMD-MFDFA runs into difficulties in reverse segmentation and the selection of the expected Intrinsic Mode Functions (IMFs). Aiming at solving these deficiencies, a Modified EMD-MFDFA (MEMD-MFDFA) approach with IMF selection strategy and Step-Moving Window (SMW) segmentation method is proposed in this paper. In MEMD-MFDFA, a metric for distinguishing deterministic and random components is established to select expected IMF components by scaling exponent. Meanwhile, SMW segmentation method is exploited to reduce the estimated errors caused by reverse segmentation. The robustness of the proposed method is investigated through comparing MEMD-MFDFA, MFDFA, and EMD-MFDFA by multifractality of simulated signals with different Signal-to-Noise Ratio (SNR). Furthermore, the proposed approach is applied to three bearing run-to-failure datasets containing three types of faults, and the results show that the multifeatures of the multifractal spectrum obtained by MEMD-MFDFA have the ability to simultaneously identify early fault and assess performance degradation of bearings.

Keywords

Empirical mode decomposition-multifractal detrended fluctuation analysis scaling exponent health condition monitoring step-moving window multifractal spectrum feature deterministic/random separation

Introduction

As one of the decisive parts responsible for the transmission of rotating machinery, the rolling bearing is more prone to failure under extreme working conditions or long operation time. Its abnormal operation or even failure is often the reason for the unplanned shutdown of rotating machinery.^1,2 Consequently, the development of a bearing health monitoring system is not only conducive to the early detection of the abnormal state of bearing in the operation process, but also provides valuable advice for predictive maintenance by assessing the severity of defect evolution.^3,4 With the purpose of monitoring health condition of bearings, it is necessary to extract effective information related to faults and construct degradation assessment indicators being capable of tracking dynamical behavior of the rolling bearing without prior knowledge of the kinematics.⁵

In practice, the establishment of health indicators mainly focuses on the extraction of statistical parameters in time domain and characteristics in frequency domain based on vibration signals.⁶ Sparsity measures including negative entropy, smoothness index, pq-mean, the ratio of Lp to Lq norm, kurtosis, and Gini index were reformulated with the ratio of different quasi-arithmetic means to monitor the health condition of bearings by Hou et al.⁷ Soualhi et al.⁸ extracted defective frequencies of different fault types from Hilbert marginal spectrum as health indicators to track the fault evolution of rolling bearings. Meng et al.⁹ constructed a health indicator on the basis of multifeatures of envelope spectrum to reflect the severity of bearings in real-time condition monitoring.

Since the fault information that is not obvious in time domain are able to be reflected in the frequency domain analysis, and vice versa, some scholars concentrated on the indicator construction based on features in time-frequency domain. The time-frequency analysis method is appropriate for revealing the time-varying features of frequency components from nonlinear and nonstationary vibration signals, which presents the partial characteristics in both time and frequency domains simultaneously.¹⁰ Although the high time-frequency resolution can be obtained using Wigner-Ville Distribution (WVD),¹¹ Short-Time Fourier Transform (STFT),¹² and Wavelet Transform (WT),¹³ the sufficient prior knowledge of the analyzed signal is required in advance. Empirical Mode Decomposition (EMD),¹⁴ as a parameter-free time-frequency method, can adaptively decompose a multi-component signal into multiple Intrinsic Mode Functions (IMFs) with fixed instantaneous frequencies and a residue, which indicates that EMD is a valid tool to deal with nonlinear and nonstationary signals.

As another feature extraction method for revealing the hidden fractal structures of vibration signals, Detrended Fluctuation Analysis (DFA) is also widely used in the processing of nonlinear signals.¹⁵ Because the inherent fractal structures of the bearing vibration signals under different health conditions are also different to some extent, it makes possible to reflect the performance degradation condition of bearings by measuring the fractal structures. However, since nonstationary and nonlinear vibration signals usually exhibit multifractal characteristics, only a mono-fractal characteristic obtained from DFA method would fail to adequately expose the nature of their dynamical mechanism.¹⁶ Although the crossover property of the scaling-law curve in DFA method can be used to demonstrate the multi-scale behaviors of complex signals and detect rotary machinery faults,^16,17 the DFA method only can reveal the total fractal characteristics of them, and lacks the depiction of their local fractal characteristics.¹⁸

In order to measure the multifractal characteristics of nonstationary time series from different local scales, the multifractal DFA (MFDFA) was proposed by Kantelhardt et al.¹⁹ Recently, MFDFA has gradually developed into a popular multifractal analysis method that can effectively uncover the multifractality of nonlinear vibration signals due to its simple execution. In view of the conspicuous difference in multifractality characteristics of vibration signals under different bearing health conditions due to different dynamic mechanisms, the extraction method based on multifractal features is applied to fault diagnosis and detection in practical application.^20,21 Furthermore, some studies focused on the application of adaptive signal decomposition methods to MFDFA, such as Local Characteristic-scale Decomposition (LCD),²² Intrinsic Time-scale Decomposition (ITD),²³ EMD,¹⁸ and Variational Mode Decomposition (VMD),²⁴ which are employed to decompose a signal into a series of sub-signals with different local scales, and improve the accuracy of fault diagnosis by extracting the inherent fractal features of this sub-signals. Although the effective behavior of multifractal characteristics in recognition of rotating machinery including diverse fault types and severity under different running conditions was verified by experimental results in the above researches, they cannot still fully explain the effectiveness of multifractal features for monitoring the health condition of the complex system in real time. Accordingly, Hu et al.²⁵ calculated multifractal characteristic parameters from vibration signals to identify the initial fault stage of the test bearing and further track its performance degradation trend during subsequent operation. Considering that vibration signals of running bearings are vulnerable to background noise in the incipient fault stage, Chen et al.²⁶ used the improved VMD algorithm to denoise the raw signal and carried out multifractal analysis on the preprocessed signal to obtain the improved multifractal features, and the results show that the extracted features are more sensitive to the incipient fault of rolling bearings. Even so, there are still some problems to be solved for the application of MFDFA in condition monitoring of bearings.^27–30 Firstly, the appropriate polynomial order needs to be chosen for the local scale-dependent trend. Secondly, scales are set empirically. Thirdly, new pseudo-fluctuation errors are produced on account of the discontinuity at the corresponding polynomial fitting trend of two adjacent segmentations. Fourthly, the order of the original time series is disrupted by reverse segmentation, which will lead to estimation errors of multifractal features.

In order to solve these problems of MFDFA, Ihlen²⁷ proposed Empirical Mode Decomposition-Multifractal Detrend Fluctuation Analysis (EMD-MFDFA), which exploited EMD to obtain smooth fitting scale-dependent trends defined by a group of IMFs with mono-component and specific time scale from a nonstationary and nonlinear signal for MFDFA. Unfortunately, the last problem of MFDFA mentioned above has not been solved. In addition, there may be pseudo-components in the IMF components acquired from EMD,³¹ which will have a negative impact on the derived multi-scale trends. Hence, Du et al.²⁸ developed a criterion to pick out the true components for fault diagnosis based on the correlation analysis between the profile of the original signal and IMF components derived by EMD in frequency domain. The results can be seen from the study of Du et al.²⁸ that the extracted multifractal features are suitable for identifying the abnormal state of the running bearing, but they cannot be used to track the fault evolution trend in the performance degradation stage. For a health indicator utilized to monitor bearing health condition, it needs to identify the early fault and simultaneously assess performance degradation of bearings. The poor performance of multifractal features based on this feature extraction method in bearing performance condition monitoring can be attributed to two reasons. On the one hand, the criterion used to distinguish the true components from IMFs is not established based on the impulse response characteristics of the fault bearing, and it is difficult to accurately distinguish the useful IMF components by the correlation criterion when the original signal is seriously contaminated by noise. On the other hand, the estimation errors of multifractal features caused by the reverse segmentation of MFDFA are overlooked.

With the purpose of addressing these two problems and extracting more sufficient information related to multifractal characteristics hidden in vibration signals for incipient fault identification and performance degradation assessment, a Modified EMD-MFDFA (MEMD-MFDFA) method is proposed based on the IMF selection strategy and Step-Moving Window (SMW) segmentation method. To begin with, a metric for selecting the expected IMF components to reconstruct signal is established in this paper by applying the scaling exponent, which is an effective tool in DFA method to measure the determinism and randomness without relying on comparison and reference to the original noisy signal. Furthermore, the SMW algorithm is introduced to substitute for the original piecewise method to avoid unnecessary errors caused by the order disruption of vibration signals. Additionally, the feasibility of MEMD-MFDFA method is verified numerically and experimentally through comparing with other multifractal analysis methods.

The remainder of this paper is organized as follows. In section “Theoretical background,” the theoretical background with respect to DFA and EMD is described concisely. In section “Modified EMD-based MFDFA algorithm,” a metric for distinguishing deterministic and random IMF components is developed based on the scaling exponent. Meanwhile, MEMD-MFDFA on the basis of IMF selection strategy and SMW segmentation method is introduced. In section “Analysis of simulated multifractal signal” and section “Application of proposed method in run-to-failure lifetime testing,” the developed MEMD-MFDFA method is exploited to analyze the simulated multifractal signal and experimental bearing vibration signals in the run-to-failure experiment, respectively. Finally, some further discussions and conclusions are given in section “Discussion” and section “Conclusion.”

Theoretical background

Detrended fluctuation analysis method

As a property hidden in the time series, long-range correlation is used to reflect the dependence strength between the time series and its time-shifted version. Compared with the Hurst exponent,³² the scaling exponent based on DFA is more appropriate for detecting and correctly quantifying the long-range correlation of nonstationary time series.³³

Considering a given time-series signal x(t) $(t = 1,2, \dots, N)$ of length N, the calculation of the scaling exponent h based on DFA³⁴ is conducted by the following procedure.

Step 1: Obtain the cumulative form Y(t) of the time-series signal x(t) by means of equation (1)

Y (t) = \sum_{k = 1}^{t} [x (k) - \bar{x}]

(1)

where

\bar{x}

denotes the mean of x(t).

Step 2: Separate the cumulative time series Y(t) into $N_{s} = [N / s]$ sections of equivalent scale s. The polynomial function is applied to fit for each section, of which the corresponding local trend Y _s(t) can be determined.

Step 3: Define the fluctuation function for different scale s as equation (2)

F (s) = \sqrt{\frac{1}{N} \sum_{t = 1}^{N} {[Y (t) - Y_{s} (t)]}^{2}}

(2)

The relationship between $F (s)$ and s is defined to determine the power-law behavior of the time series, which is analyzed by calculating the slope h of the $F (s)$ versus s in log-log scale as equations (3) and (4)

F (s) \propto s^{h}

(3)

\ln (F (s)) = h \ln (s) + C

(4)

where the slope h is called the scaling exponent in DFA method.

As reported in the study of Dong et al.,³⁵ the correlation of the analytic time series x(t) is closely related to the relationship between $F (s)$ and s. When h locates the range from 0 to 0.5, the time series always exhibits short-range correlation. While $0.5 < h < 1$ , temporal correlation of the time series is persistent. In the case of $h > 1$ , the time series is characterized by strong non-stationarity. Therefore, the smoothness of the analytic time series can be quantified through the scaling exponent h. Specifically, the smaller the value of h is, the rougher the time series is, and vice versa.¹⁷

Empirical mode decomposition algorithm

EMD is a data-driven adaptive multi-component signal decomposition tool introduced by Huang et al.,¹⁴ which is capable of decomposing nonlinear and nonstationary signals into IMFs and a residue so that their sum is equal to the original signal. In EMD, IMFs need to satisfy two conditions:³⁶ First of all, the number of local extreme points and over-zero points of the whole data are required to be equal to or at most one difference from each other. The second is that the average of the two envelopes estimated by the local maxima and minima is supposed to be zero. After the sifting process,³¹ the all obtained IMFs $y_{i} (t)$ and the residue $r (t)$ can be summed up to reconstruct the raw signal $x (t)$ as equation (5)

x (t) = \sum_{i = 1}^{n} y_{i} (t) + r (t)

(5)

The Hilbert transform is performed on each IMF to obtain the instantaneous amplitude and instantaneous frequency. In the Hilbert transform theory,³⁷ a complex-valued analytic signal can be constructed from a real signal by an exponential term and its real envelope. Therefore, a complex-valued analytic signal $y_{A i} (t)$ of an IMF $y_{i} (t)$ is described by equation (6)

y_{A i} (t) = A_{i} (t) e^{j ϕ_{i} (t)}

(6)

where

A_{i} (t)

and

ϕ_{i} (t)

stands for the instantaneous amplitude and the phase of

y_{A i} (t)

given as equations (7) and (8)

A_{i} (t) = \sqrt{y_{i}^{2} (t) + {\tilde{y}}_{i}^{2} (t)}

(7)

ϕ_{i} (t) = arctg (\frac{{\tilde{y}}_{i} (t)}{y_{i} (t)})

(8)

where

{\tilde{y}}_{i} (t)

denotes the Hilbert transform of

y_{i} (t)

The instantaneous frequency corresponding to each IMF is defined by equation (9)

ω_{i} (t) = \frac{d ϕ_{i} (t)}{d t}

(9)

Modified empirical mode decomposition-multifractal detrended fluctuation analysis algorithm

Overlap segmentation based on step-moving window

In the studies of temperature series,³⁸ it is considered that points on a temperature series usually have a better correlation with points close to them. Likewise, the relationship would hold for points in the time series. Therefore, in time series data processing, the dependence of data at different time points is very important.³⁹

As depicted in Figure 1(a), $2 N_{s} = 2 [N / s]$ non-overlapping sections of equivalent length s denoted as $ε (v, s) (v = 1,2, \dots, 2 N_{s})$ are derived by segmenting the scale-dependence fluctuation $ε_{s} (t)$ with a length of N in opposite order, respectively, where notation $[\cdot]$ signifies the downward integral function. This segmentation method ignores the fact that the points on one local segment also have a close relationship with those on another local segment next to them. It is more problematic especially that points located on both sides of the local segment may not be related to the points in the same local segment owing to the long distance. It is natural to expect to change the segmentation by introducing SMW so that the points at these particular positions on one segment are in the middle of another segment. The detailed procedure with respect to SMW algorithm shown in Figure 1(b) is described as follows.

Figure 1.

Schematics of different segmentation methods. (a) Non-overlapping division in opposite order, (b) Overlap division with SMW.

Step 1: Take the time scale s as the width of the moving window.

Step 2: Set the step length of forward movement as $S L = 1$ for the moving window.

Step 3: The fixed-width moving window is applied on the scale-dependence fluctuation $ε_{s} (t)$ to acquire $N - s + 1$ overlapping segments with length s.

The correlation is enhanced by placing the adjacent points on the same segment of time series, which helps to obtain a reliable local fluctuation of time series and extract a more accurate trend. The segmenting method based on the SMW algorithm avoids the divisibility between the length N of time series and the time scale s, so it makes full use of the information provided by the time series. Additionally, there is not order disruption of the original time series caused by reverse segmentation.

Evaluation indicator for intrinsic mode function selection

The scaling exponent h of multifractal simulation signal under different SNR that arises from the p-model algorithm in the study of Meneveau et al.⁴⁰ is calculated as shown in Figure 2. It can be concluded from the constant scaling exponent that the influence of the small amount of noise on the long-range correlation of the original signal is very small when the SNR is large. However, with the continuous reduction of the SNR, the scaling exponent h decreases continuously due to the destruction of the long-range correlation of the signal caused by the gradually strong noise.

Figure 2.

Scaling exponent h of multifractal simulation signal with different SNR.

With the assistance of EMD, signals with SNR of 60–80 dB and 10–30 dB are decomposed into multiple IMF components with specific time scales, of which scaling exponent h are computed as depicted in Figure 3. From Figure 3, it can be observed that the trends of the scaling exponent h at both ends are steady, which implies that the IMF components in these two regions have a good correlation, respectively. Meanwhile, each of the IMF components located in the ascending tendency of the scaling exponent h shows varying degrees of correlation with other components. By comparing Figures 3(a) and (b), it is obvious that the similar IMF components with smaller scaling exponent h are more sensitive to background noise and exhibit much difference in long-range correlation with the variance of noise intensity. On the contrary, in the case of different SNR, the similar IMF components with larger scaling exponent h still maintain basically unchanged long-range correlation, which means that the information on the deterministic original signal is retained in the IMF components with larger scaling exponent. Furthermore, it can be reasonably inferred that the random noise is more concentrated in the IMF components with smaller scaling exponent h. As a result, the deterministic and random components could be effectively distinguished by means of the scaling exponent h.

Figure 3.

Scaling exponent h of IMF components from multifractal simulation signal with different SNR in the area of: (a) 60–80 dB, (b) 10–30 dB.

In the signal processing for recovering the diagnostic information, the components with different properties in the vibration signal of the bearing are supposed to be separated in time. A more universal splitting at present is between deterministic and random components.⁴¹ The repeated shocks caused by a local defect in rolling bearings do not exactly possess periodicity. The reason is that the rolling elements are not quite uniformly placed in the clearance of the cage and the rotating speed of the cage may be unstable due to slipping in the process of operation. Therefore, it is of great significance for fault diagnosis to separate the bearing signals with random impact characteristics from the deterministic signals produced by other components.⁴²

For the randomly oscillating IMFs in the vibration signals, they will be regarded as useful components in the multifractal analysis of bearing vibration signals on the basis of MEMD-MFDFA. At the same time, the deterministic IMF components with weak oscillation are discarded as noise to enhance fault signals in terms of multifractal characteristics.

In the denoising method, based on the adaptive decomposition of signals, the most critical step is to provide an appropriate threshold for the selected indicator to recover the contaminated signals. According to the above analysis, the maximum scaling exponent corresponding to the obtained IMF components is almost not affected by noise, and can be used as a reference value to judge the effectiveness of the IMF component. Thus, the threshold $θ$ on the strength of maximum scaling exponent is defined as equation (10) in this paper

θ = η \max (h_{i}) (i = 1,2, \dots, n)

(10)

where

h_{i}

represents the scaling exponent of the i-th IMF and

η

denotes the ratio factor in the range of 0–1.

The larger the ratio factor $η$ , the greater the proportion of the deterministic components in the raw signal.

Modified empirical mode decomposition-multifractal detrended fluctuation analysis algorithm

The proposed MEMD-MFDFA consists of five stages as follows.

Stage 1 Adaptive decomposition based on EMD

Suppose that x(t) $(t = 1,2, \dots, N)$ is a time-series signal of length N. For the sake of eliminating the constant offset over different time scales s that occurs in the time series, a cumulative series Y(t) is constructed as equation (1).⁴³

EMD is employed to decompose the cumulative series Y(t) into several IMFs $y_{i} (t) (t = 1,2, \dots, n)$ with specific physical significance and a residual component r(t) that reveals the average trend hidden in the time series according to the sifting process.¹⁴ Theoretically speaking, all obtained IMFs and the residue r(t) can be summed up to obtain the reconstructed cumulative series $Y_{r} (t)$ as equation (5).

Then, the scaling exponent h corresponding to each IMF component is compared with the threshold $θ$ . In accordance with the metric for IMF selection, $Y_{r} (t)$ can also be expressed by equation (11)

\begin{array}{l} Y_{r} (t) = \sum_{j \in A} y_{j} (t) + \sum_{k \in B} y_{k} (t) + r (t) = \sum_{j \in A} y_{j} (t) + \tilde{r} (t) \\ A = {A | h_{A} \leq θ}, B = 1 - A \end{array}

(11)

where

h_{A}

denotes the scaling exponent of noise-free IMF components.

\tilde{r} (t)

is considered as the new residue.

Stage 2 Extraction of scale-dependent fluctuation

Since the time-series signal is decomposed in terms of the local characteristic time scale through EMD, a group of IMFs with mono-component and specific time scale are acquired. It should be noted that the time scale s of each IMF component can reflect its own oscillation pattern, which has a reciprocal relationship with the mean of instantaneous frequency $ω (t)$ of IMFs as defined by equation (12)²⁷

s = \frac{1}{ω}

(12)

where

\bar{ω}

represents the mean of

ω (t)

In order to discriminate $y_{j} (t)$ mentioned above, $d_{s} (t)$ is used to represent the oscillation component with time scale s. Then the time scale sequence $s = [s_{1}, s_{2}, \dots, s_{j}]$ of noise free IMFs is formed. The combination of IMFs with a larger time scale than $s_{j}$ and the residue is defined to the scale-dependence trend ${\hat{c}}_{s} (t)$ for each scale s as equation (13)

{\hat{c}}_{s} (t) = \sum_{s > s_{j}} d_{s} (t) + \tilde{r} (t)

(13)

Subsequently, the scale-dependence fluctuation $ε_{s} (t)$ at corresponding time scale s could be defined as equation (14)

ε_{s} (t) = Y_{r} (t) - {\hat{c}}_{s} (t)

(14)

Stage 3 Construction of detrended fluctuation function

The scale-dependence fluctuation $ε_{s} (t)$ is divided into $N - s + 1$ overlapping sections of equivalent length s. Consequently, a total of $N - s + 1$ sections with length s denoted as $ε (v + i - 1)$ are derived. The square of the detrended fluctuation function $F (v, s)$ is calculated for each of $ε (v + i - 1)$ as equation (15)

F^{2} (v, s) = \frac{1}{s} \sum_{i = 1}^{s} {[ε (v + i - 1)]}^{2} (1 \leq v \leq N - s + 1)

(15)

The fluctuation function $F_{q} (s)$ of the q-th order can be calculated by averaging over all the segments as described in equation (16)

F_{q} (s) = {\begin{cases} {(\frac{1}{N - s + 1} \sum_{v = 1}^{N - s + 1} {[F^{2} (v, s)]}^{\frac{q}{2}})}^{\frac{1}{q}} (q \neq 0) \\ \exp [\frac{1}{2 [N - s + 1]} \sum_{v = 1}^{N - s + 1} \ln F^{2} (v, s)] (q = 0) \end{cases}

(16)

Stage 4 Calculation of multifractal spectrum

The logarithmic relationships between $F_{q} (s)$ and s can be analyzed to reveal the scaling behavior of $F_{q} (s)$ for different values of q given as equation (17)

F_{q} (s) \sim s^{h (q)}

(17)

where

h (q)

denotes the generalized Hurst exponent.

Here the least square method is applied on the log-log graphs $F_{q} (s)$ versus s to obtain its slope regarded as the $h (q)$ in the multifractal analysis. Afterward, the corresponding multifractal scaling exponent $τ (q)$ can be calculated based on equation (18) for each q

τ (q) = q h (q) - 1

(18)

With the assistance of Legendre transform, two parameters of singularity exponent

α

and multifractal spectrum

f (α)

reflecting the multifractal characteristics of time series can be expressed by

τ (q)

as equation (19) and (20)

α (q) = τ^{'} (q) = h (q) + q h^{'} (q)

(19)

f (α) = q α - τ (q) = q [α - H (q)] + 1

(20)

In the fault diagnosis of rolling bearings based on multifractal characteristics, the complex fractals of the bearing vibration signal can be divided into small regions with different degrees of singularity by multifractal method, and the fault states of the bearing can be found by analyzing their long-range correlation. Singularity exponent $α$ reflects the non-uniformity strength of the signal in the local probability measure distribution. Multifractal spectrum $f (α)$ denotes the fractal dimensions of segments with the same local roughness strength, which describes the distribution of some fractal exponents hidden in the vibration signal.

Analysis of simulated multifractal signal

Aiming at demonstrating the effectiveness of the proposed method for multifractal characteristic analysis, the simulated noisy multifractal signal generated through the p-model in the study of Meneveau et al.⁴⁰ is analyzed by the proposed method in this section. The multifractal signal can be constructed by means of the equal division of unit interval based on this multifractal model. The obtained two sections with the same length will be given different weights, that is, p ₁ and p ₂. Then the same step above is repeated for each of the sections. Finally, when the number k of repetition is equal to m, the generated multifractal time series $P (k)$ would be composed of $N = 2^{m}$ pieces of equal sections, which is defined as equation (21)

P (k) = p_{1}^{(m - k)} p_{2}^{k} (k = 0,1, \dots, m)

(21)

where m denotes the maximum number of repetition.

p_{1} = 1 - p_{1}

and

0.5 < p_{2} < 1

Herein, we use $p_{1} = 0.3$ , $p_{1} = 0.7$ , and $m = 16$ in the p-model. In theory, the generalized Hurst exponent $h (q)$ and the multifractal scaling exponent $τ (q)$ corresponding to this p-model could be expressed as equations (22) and (23)

h (q) = \frac{1}{q} - \frac{\ln [p_{1}^{q} + p_{2}^{q}]}{q \ln 2}

(22)

τ (q) = - \frac{\ln [p_{1}^{q} + p_{2}^{q}]}{\ln 2}

(23)

Similarly, the singularity exponent $α (q)$ and multifractal spectrum $f (α)$ of this p-model could be calculated according to equations (24) and (25)

α (q) = - \frac{p_{1}^{q} \ln p_{1} + p_{2}^{q} \ln p_{2}}{\ln 2 (p_{1}^{q} + p_{2}^{q})}

(24)

f (α) = - \frac{q p_{1}^{q} \ln p_{1} + q p_{2}^{q} \ln p_{2}}{(p_{1}^{q} + p_{2}^{q}) \ln 2} + \frac{\ln (p_{1}^{q} + p_{2}^{q})}{\ln 2}

(25)

The theoretical values of four multifractal characteristic parameters could be utilized to compare the results of several MFDFA-based employed in multifractal simulation signals. Meanwhile, two evaluation indicators are employed to appraise the performance of these methods, which is beneficial to avoid the subjective judgment in the comparison of analysis results. The average absolute error $δ$ between the estimated value $E (i)$ and the theoretical value $T (i)$ can be used as an evaluation indicator to represent the accuracy of the multifractal analysis method as equation (26)

δ = \frac{1}{n} \sum_{i = 1}^{n} | E (i) - T (i) |

(26)

In addition, the variance $σ$ corresponding to the errors of all estimated points is calculated to quantify the degree of deviation from the average absolute error, which can be used as another indicator to evaluate the stability of the analysis method as equation (27)

σ = \frac{1}{n} \sum_{i = 1}^{n} {[E (i) - T (i) - δ]}^{2}

(27)

Considering that the IMFs in the ascending stage of the scaling exponent h may contain a large number of deterministic and/or random components, it is necessary to obtain more sufficient and effective multifractal characteristics through qualitative evaluation and separation of them. As these IMFs with uncertain components are continuously mixed with the random components in the increasing order of scaling exponent, the average absolute errors between the estimated and theoretical $α (q)$ are calculated for the multifractal signals corresponding to SNR of 30–40 dB, which is shown in Figure 4. As depicted in Figure 4, the average absolute error $δ$ for each of the diverse noise intensity distinctly reveals a crossover point at the same modal number. Consequently, these IMF components could be roughly split into two categories using the crossover point. In addition, it is found that the scaling exponent $h_{cross}$ corresponding to the crossover points of these results are very close in value, so the ratio factor $η$ related to the threshold setting can be determined by their ratios to the maximum scaling exponent as described by equation (28)

η > \frac{h_{cross}}{\max (h_{i})} (i = 1,2, \dots, n)

(28)

Figure 4.

Relationship of average absolute error $δ$ and number of uncertain modes under different SNR.

The ratio between $h_{cross}$ and the maximum scaling exponent for multifractal simulated signals with different SNR is shown in Figure 5. It can be observed in Figure 5 that the golden ratio 0.618 is greater than the ratio between $h_{cross}$ and the maximum scaling exponent, so it is used as the ratio factor $η$ in this paper.

Figure 5.

Ratio between $h_{cross}$ and the maximum scaling exponent for multifractal simulated signals with different SNR.

MFDFA, EMD-MFDFA, and the proposed method are also applied to estimate $h (q)$ , $τ (q)$ , , $f (α)$ and $α (q)$ of the noisy multifractal signal with SNR of 20 dB for comparing their performance. Different order $γ$ in the range of [1,3] of the polynomial function in MFDFA is used in this paper, which are denoted as MFDFA1, MFDFA2, and MFDFA3. Additionally, the exponent order q of these five MFDFA-based methods is set as −10 to 10 in this paper. The estimated values of the four multifractal characteristic parameters obtained by five MFDFA-based methods are compared with their corresponding theoretical values as shown in Figure 6. It can be seen in Figures 6(a) and (b) that the results of these five approaches are very close to the theoretical calculation for the positive q. However, when q is equal to a negative value, except for the MEMD-MFDFA method, the estimated generalized Hurst exponent $h (q)$ and the multifractal scaling exponent $τ (q)$ of the other four approaches are seriously inconsistent with the theoretical value. As shown in Figures 6(a) and (c), the relationship between $τ (q)$ and q is nonlinear and the unimodal convex curve derived through theoretical computation appears in the multifractal spectrum, which confirms the multifractal characteristic of the analyzed signal. Different from the MFDFA1, MFDFA2, and MFDFA3 easily affected by noise, the proposed method is more capable of suppressing the interference of white Gaussian noise to reveal the multifractal characteristics emerged in the contaminated signal. Although EMD-MFDFA can do them as shown in Figures 6(a) and (c), its accuracy in estimating the multifractal parameters seems to be poor. It can be observed from Table 1 that the average absolute error $δ$ of the estimated $h (q)$ , $τ (q)$ , $f (α)$ , and $α (q)$ is equal to 1.24, 0.23, 0.22, and 0.27, which is obviously inferior to 0.13, 0.05, 0.05 and 0.04 of the proposed method. A similar conclusion can still be acquired in the comparison of the results with respect to another evaluation indicator $σ$ .

Figure 6.

Relationships of (a) $τ (q) \sim q$ , (b) $h (q) \sim q$ , and (c) of noisy p-model time series with $S N R = 20 dB$ .

Table 1.

Evaluation indicators of multifractal parameters for noisy p-model time series with $S N R = 20 dB$ .

	$τ (q)$		$h (q)$		$α (q)$		$f (α)$
	$δ$	$σ$	$δ$	$σ$	$δ$	$σ$	$δ$	$σ$
MFDFA1	0.44	0.14	2.51	8.73	0.55	0.21	0.48	0.16
MFDFA2	0.44	0.14	2.42	8.33	0.52	0.20	0.37	0.12
MFDFA3	0.42	0.13	2.34	7.97	0.50	0.19	0.31	0.11
EMD-MFDFA	0.23	0.05	1.24	2.61	0.27	0.07	0.22	0.05
Proposed method	0.05	1.34e−4	0.13	0.014	0.04	6.23e−4	0.05	2.80e−3
Proposed method without SMW	0.06	1.78e−4	0.15	0.018	0.09	1.09e−3	0.08	4.74e−3

Additionally, the proposed method without SMW is also employed to estimate $h (q)$ , $τ (q)$ , $f (α)$ , and $α (q)$ of this noisy multifractal signal and the results are presented in Figure 6. Compared with the estimated values of these multifractal parameters obtained by other existing approaches, the computation results of the proposed method without SMW are closer to the theoretical values, which suggests that the proposed IMF selection strategy based on the scaling exponent has the ability to identify and remove the pseudo-components from the noisy signal, so that it is more effective than other existing methods in measuring the multifractal characteristics of analyzed signal under strong noise. Meanwhile, making comparison between the results of the proposed method with or without SMW shown in Figure 6, the multifractal parameters calculated by the proposed method with SMW are in better accordance with the theoretical values, especially in the calculation of multifractal spectrum of Figure 6(c). As reported in Table 1, the average absolute error $δ$ and the variance $σ$ between the estimated $α (q)$ according to the proposed method with SMW and its theoretical value are 0.04 and 6.23e-4, respectively, while those of the proposed method without SMW are about 0.09 and 1.09e-3, respectively. The average absolute error $δ$ and the variance $σ$ for $α (q)$ are reduced by 55.56% and 42.84% in terms of the proposed method with SMW compared to that without SMW, and there are also some improvements in $h (q)$ , $τ (q)$ , and $f (α)$ . Consequently, it implies that SMW segmentation method also has enormous influence on the performance of the proposed method.

To further verify the anti-noise performance of the proposed MEMD-MFDFA method, the multifractal scaling exponent $τ (q)$ is taken as an example. The capability of the proposed method in extracting multifractal features under different SNR is compared with other four methods by calculating their corresponding average absolute error $δ$ and the variance $σ$ , respectively. The changing trends representing the average absolute error $δ$ and the variance $σ$ with different SNR are displayed in Figure 7. As can be seen from Figure 7, both $δ$ and $σ$ of MFDFA2 and MFDFA3 increase significantly when SNR is greater than 70 dB, which indicates that the analysis results of these two methods are more susceptible to the influence of noise. Although the values of $δ$ and $σ$ acquired by MFDFA1 with the SNR between 55 and 80 dB are smaller than those of MFDFA2 and MFDFA3, its average absolute error and variance are higher than 0.22 and 0.0321, respectively, when the SNR is less than 50 dB. Thus, choosing the appropriate order of polynomials may obtain more reliable multifractal features for the polynomial-based multifractal analysis method.

Figure 7.

Evaluation indicators of analysis results: (a) Average absolute error $δ$ (b) Variance $σ$ .

Furthermore, it can be seen from Figure 7 that with the decrease of SNR, both average absolute error and variance of EMD-MFDFA are increasing. In comparison with EMD-MFDFA, both the variation trends of average absolute error and variance of the proposed method are steady. In addition, the resultant $δ$ of EMD-MFDFA are larger than 0.1 when $S N R < 45 dB$ , while the average absolute error $δ$ of the proposed method is just over 0.1 in the case of the lowest SNR. Similarly, the advantages of the proposed method can also be reflected in the variance of another indicator and the same conclusion can be drawn. It is obvious that the IMF selection strategy based on the scaling exponent is conductive to the accurate extraction of multifractal features.

Application of proposed method in run-to-failure lifetime testing

Experimental set-up

Three run-to-failure datasets for bearings containing outer ring fault, inner ring fault and ball fault, respectively, were provided by the Center for Intelligent Maintenance Systems (IMS) of Rexnord Corp.,⁴⁴ which are used to verify the effectiveness of multifractal characteristics extracted through the proposed method in reflecting the natural degradation of bearings under different health conditions. On the accelerated fatigue test bench shown in Figure 8, four Rexnord ZA-2115 double-row bearings, which were uniformly mounted on the same shaft with the constant rotation speed of 2000 rpm driven by an AC motor, were subjected to a 6000 lbs radial force generated by a spring mechanism. The vibration signals with data length of 20,480 points were collected every 10 min in terms of PCB 353B33 High Sensitivity Quartz ICP accelerometers with sampling frequency at 20 kHz. The detailed description of the dataset corresponding to the two cases studied in the presented work is shown in Table 2.

Figure 8.

Test bench of IMS. (a) Experimental system, (b) System diagram.

Table 2.

Description of dataset from different faulty bearing.

Case	Dataset	Fault type	Number of data groups
Case 1	Bearing 3-Ch 5	Inner race fault	2156
Case 1	Bearing 4-Ch 8	Ball fault	2156
Case 2	Bearing 1-Ch 1	Outer race fault	984

Experiment Results

The accelerated fatigue test of the bearing continued for nearly a week until an obvious failure occurred at the outer ring of Bearing 1-Ch 1, of which the vibration signal is displayed in Figure 9. Based on the previous work,⁹ Growth Rate of Real-time Mahalanobis Distance with Cumulative Sum (GRRMD-CUMSUM) was proved to be able to effectively distinguish the different health states during the entire operation by calculating the rate of performance degradation of the bearing based on the features of envelope spectrum. Consequently, this indicator is used to reflect the degradation process throughout the life cycle of Bearing 1-Ch 1.

Figure 9.

Vibration signal of Bearing 1-Ch 1 in time domain.

As depicted in Figure 10, four phases including normal phase, incipient fault phase, severe fault phase, and failure phase are highlighted for Bearing 1-Ch 1 according to the change trend of GRRMD-CUMSUM. Twenty data groups are then randomly selected from each of identified degradation stages to carry out multifractal analysis. Subsequently, these bearing vibration signals are analyzed in the light of the proposed method. As a result, the calculated multifractal spectrums are displayed in Figure 11(a). It can be seen from Figure 11(a) that for negative q, the multifractal spectrum corresponding to the failure stage is distinguished from the rest but that corresponding to the other three stages cannot well distinguish different health conditions during the performance degradation process of the bearing, which indicates that it is not significant to extract the features from the multifractal spectrum at negative q for the condition monitoring. The multifractal spectrums in the range of $q < 0$ at the failure stage have a better distinction from those at other stages, which may be helpful for the detection of bearing failure.

Figure 10.

GRRMD-CUMSUM of test data.

Figure 11.

(a) Complete multifractal spectrums, (b) Partial multifractal spectrums at positive q of different performance degradation stages based on proposed method.

Partial multifractal spectrums at positive q of different performance degradation stages based on the proposed method are shown in Figure 11(b). Normal, incipient fault and severe fault conditions can be effectively identified in accordance with the location differences of the partial multifractal spectrums as depicted in Figure 11(b). Different from the multifractal features extracted from bearing vibration signals in the first three stages, it can also be observed in Figure 11(b) that the multifractal spectrums corresponding to the failure stage have serious dispersion, and even overlap significantly with those of severe fault stage. The irregular contact between the fragmented outer ring raceway surface and other components in the failure stage results in the instability of multifractal features. Therefore, the horizontal and vertical coordinates of each point on the partial multifractal spectrum at positive q, namely, the singularity exponent $α (q)$ and multifractal spectrum $f (α)$ , can be, respectively, used as indicators to characterize the performance degradation trend of bearings.

$f (α)$ and $α (q)$ with different $q \in [0, 10]$ are, respectively, derived for each data group of Bearing 1-Ch 1 on the basis of their multifractal spectrums. In addition, to compare with the results in the case of positive q, $q = - 1$ , as one of the cases of negative q, is also used to study the change trend of its corresponding singularity exponent $α (q)$ and multifractal spectrum $f (α)$ in the process of bearing fault evolution, which are presented in Figure 12.

Figure 12.

Extraction results of multifractal features: (a) Change trend of singularity exponent $α (q)$ under different exponent order q, (b) Zoom-in view of figure (a), (c) Change trend of multifractal spectrum $f (α)$ under different exponent order q in the range of −1–4, (d) Change trend of multifractal spectrum $f (α)$ under different exponent order q in the range of 5–10, and (e) Zoom-in view of figure (d).

It can be seen from the change trend of $α (q)$ in Figure 12(a) that the performance degradation evolution of the bearing can be segmented into four different stages, which is similar to that as depicted in Figure 10.

For stage I, the changes of $α (q)$ are very slight, which are consistent with the normal condition of the healthy bearing. Thus, the bearing can be considered to be in healthy condition based on the variance trend of $α (q)$ . The period of the normal condition has lasted for 88.67 h, which exposes that for the whole life cycle of a normal bearing, the long running period of the bearing is in this stage.

For stage II, when a small defect is generated, the trend of $α (q)$ slightly decreases at the beginning of 88.83 h. As the defect size increases, the trend of $α (q)$ exhibits decreasing with the running time. Because the fault features are often masked by the strong background noise in the initial fault stage, it is difficult to confirm the start time of this stage.

For stage III, the trailing edge of a small defect would become smoother under the mutual effects of abrasive wear, spalling, and over-rolling behavior,⁴⁵ which may be cause the healing phenomenon in the performance degradation of bearings. The amplitude of $α (q)$ rises and falls irregularly until the running bearing enters into the stage IV.

For stage IV, the bearing continues to operate at a rapidly increasing rate of performance degradation in the failure stage. Although the trend of $α (q)$ presents slight oscillation during this period, its overall trend is still decreasing, which is able to reflect the severity of performance degradation of the bearing at this stage.

Additionally, some results are provided as follows on the basis of Figures 12(a) and (b).

1) When q starts to increase from 1, the trend of $α (q)$ can manifest a stronger response when the bearing performance changes, that is, the global change tendency is more conspicuous, which helps to classify the performance degradation stage.

2) When q increases to 4, the larger q is not much assistance to improve the classification results of the bearing performance degradation stage.

3) There are several local fluctuations highlighted by red circles in the same position of the trends of $α (q)$ , which are depicted in Figure 12(b). The larger q is, the larger the peak at these local fluctuations is.

4) Relative to the singularity exponents $α (q)$ at positive q, $α (q)$ in the case of $q = - 1$ does not have the ability to identify the first three performance degradation stages of the bearing.

According to equation (20), when q is equal to 0, the multifractal spectrum $f (α)$ is 1, thus the multifractal spectrum of q being close to 0 is taken as the value of $f (α)$ in the case of $q = 0$ in this paper. In addition, due to the obvious difference in the amplitudes of $f (α)$ under different q, “Min-Max scaling” as a normalization technique is used to map $f (α)$ to a specific interval [0,1].⁴⁶ As a result, the trend of $f (α)$ with respect to accumulated running time is derived to monitor the performance degradation conditions of Bearing 1-Ch 1 as depicted in Figures 12(c) and (d).

As shown in Figure 12(c), when q is in [−1, 2], the trends of multifractal spectrum $f (α)$ calculated at each performance degradation stage are obvious, which is conducive to easily identifying different health conditions of the test bearing. However, the variance trend of $f (α)$ in the case of larger q becomes inconspicuous attributing to its violent fluctuation, which especially leads to the inability to accurately estimate the entry time of the incipient fault stage. Thus, taking $q = 5$ as an example, the local fluctuation of $f (α)$ can be eliminated by calculating the mean value of every 30 neighborhoods, so as to obtain the expected trend, which is shown in Figure 12(d). As can be observed from Figure 12(d), when the trend transits from stage I to stage II, the value of $f (α)$ increases from 0.31 to 0.35, which indicates that the initial fault of the running bearing occurs. The obvious fluctuation of the performance degradation trend in the stage III may be caused by the fact that the defect edge surface is alternately smooth or rough under frequent contact with the rolling element. Until the bearing enters the failure stage, its performance degradation trend increases rapidly. From the zoom-in view of the $f (α)$ trends at q in [5, 10] as depicted in Figure 12(e), with the increase of q, the trend of $f (α)$ fluctuates more violently in the local area, which means that it is more difficult to reflect the health condition of running bearings according to their trends.

In addition, the other two run-to-failure datasets containing different bearing fault types are also utilized to verify the performance of multifractal features in condition monitoring of the running bearing. The singularity exponent $α (q)$ at $q = 1$ or $q = 2$ are taken as examples to track the evolution trends of inner race defect in Bearing 3-Ch 5 and ball defect in Bearing 4-Ch 8. Their results are displayed in Figure 13 and Figure 14. It can be observed in Figure 13(a) and Figure 14(a) that the $α (q)$ present obvious downward trends when the running bearings enter into the initial fault stage. The downward trends of $α (q)$ do not continue in the whole time. With the change of race surface, the trends of $α (q)$ rise or fall irregularly until they drop again rapidly and continuously, which indicates the complete failure of the bearing. Hence, the characteristic parameters on the multifractal spectrum obtained by the proposed method can be used to simultaneously identify early fault and assess the bearing performance degradation.

Figure 13.

(a) Trend analysis, (b) Condition monitoring (after 200 h of operation) of singularity exponent $α (q)$ at $q = 0$ for Bearing 3-Ch 5.

Figure 14.

(a) Trend analysis, (b) Condition monitoring of singularity exponent $α (q)$ at $q = 1$ for Bearing 4-Ch 8.

Discussion

Both MFDFA1 and EMD-MFDFA are also employed to extract multifractal features by the same dataset of Bearing 1-Ch 1. The multifractal spectrum obtained under different bearing health conditions can be utilized to evaluate the ability of the multifractal features extracted by these analysis methods in distinguishing the different performance degradation conditions of bearings.

Similar to the above experiment, the 20 data groups randomly selected from each of identified degradation stages are analyzed by MFDFA1 and the results are shown in Figure 15. As shown in Figure 15(a), the multifractal spectrums at negative q of the bearing under four conditions overlap each other and are difficult to be distinguished, which indicates that multifractal characteristic parameters at negative q cannot reveal the differences in the dynamic behaviors of the bearing under different conditions. From the partial multifractal spectrum at positive q shown in Figure 15(b), it can be seen that the normal condition can be well distinguished from other abnormal conditions, but there is a slight overlap between incipient and severe fault stages.

Figure 15.

(a) Complete multifractal spectrums, (b) Partial multifractal spectrums at positive q of different performance degradation stages based on MFDFA1.

Also, these selected vibration signals are analyzed by EMD-MFDFA and the results are displayed in Figure 16. As shown in Figure 16(a), for negative q, it is hard to distinguish the four stages of the running bearing by their multifractal spectrums, which is similar to the results obtained from MFDFA1. Compared with the proposed method and MFDFA1, Figure 16(b) shows that the partial multifractal spectrum at positive q calculated by EMD-MFDFA under different bearing health conditions have more serious overlap with each other. Nevertheless, the results of EMD-MFDFA also manifest that the partial multifractal spectrum at positive q is more suitable than that at negative q to reflect the fault evolution of bearings in the whole life cycle, which is consistent with the results of the proposed method and MFDFA1.

Figure 16.

(a) Complete multifractal spectrums, (b) Partial multifractal spectrums at positive q of different performance degradation stages based on EMD-MFDFA.

As shown in Figure 17(a), four different health conditions of bearings can be discriminated from each other in the trend of singularity exponent $α (q)$ obtained from MFDFA1. However, according to the trend of $α (q)$ for EMD-MFDFA in Figure 17(b), the initial fault stage is indistinguishable from the rest and the other three stages can only be roughly distinguishable from each other, which indicates that EMD-MFDFA cannot reveal the multifractality contained in the bearing vibration signals due to the interference of heavy background noise in the early fault stage. The performance degradation trends of running bearings described by multifractal spectrum $f (α)$ with different exponent order q are also shown in Figure 18 and Figure 19. It can be seen in Figure 18 that according to the performance degradation trends characterized by $f (α)$ of MFDFA1, the severe fault stage and the failure stage of the bearing can be easily distinguished from the other two stages, but the division between the early fault stage and the normal stage is hard to determine in Figure 18. The $f (α)$ of EMD-MFDFA corresponding to q in [-1, 4] cannot be utilized to determine the specific time owing to severe fluctuation when the running bearing enters the incipient fault stage as depicted in Figure 19(a). It can also be observed from Figure 19(b) that the ambiguous trend of $f (α)$ corresponding to q = 5 is not conducive to dividing the health conditions of the bearing.

Figure 17.

Trend of singularity exponent $α (q)$ under different exponent order q by (a) MFDFA1, (b) EMD-MFDFA.

Figure 18.

Trend of multifractal spectrum $f (α)$ obtained by MFDFA1 under different exponent order q in the range of: (a) −1–4, (b) 5–10.

Figure 19.

Trend of multifractal spectrum $f (α)$ obtained by EMD-MFDFA under different exponent order q in the range of: (a) −1–4, (b) 5–10.

The normal stage is well known as the comparatively steady state of the whole life cycle of the running bearing, which is significantly important to recognize the anomalies due to the occurrence of small defect or crack area. To further prove the superiority of the proposed method in monitoring the health condition of bearing in stage I, MFDFA1 and EMD-MFDFA are applied to Bearing 3-Ch 5 and Bearing 4-Ch 8 for comparison and the results are shown in Figures 20 and 21. In accordance with the analysis of the health conditions of Bearing 3-Ch 5 and Bearing 4-Ch 8 by MFDFA1 and EMD- MFDFA, it can be observed in the stage I that the singularity exponent $α (q)$ shows a downward/upward trend at the beginning of bearing operation as depicted in Figures 20 and 21, which is in keeping with the results of the proposed method shown in Figures 13(a) and 14(a). In fact, from the perspective of manufacturing process differences and rolling contact fatigue, stage I corresponding to the normal stage of bearings can also be divided into running-in stage and steady-state stage.⁴⁵ In the running-in stage, the increase of material strength and microyield stress caused by work hardening will produce residual stress, which will reduce the mechanical properties of the metal surface of the component.⁴⁷ The surface hardness continues to increase until the work hardening reaches saturation, which represents that the bearings enter into steady-state stage. In the steady-state stage, the bearings are expected to run in a stable manner under the normal operation condition. Through the comparison between Figure 20 and Figure 13(a), the singularity exponent $α (q)$ obtained by EMD-MFDFA and the proposed method exhibit smaller variation in the steady-state stage of Bearing 3-Ch 5. Nevertheless, for Bearing 4-Ch 8, condition assessment curves, respectively, indicated by $α (q)$ of MFDFA1, EMD-MFDFA and the proposed method show varying degrees of fluctuation in steady-state stage as shown in Figure 21 and Figure 14(a). From the fault location of Bearing 4-Ch 8, it can be inferred that the rolling body of this bearing is more sensitive to contact stress than other parts in loading zone, so its surface is more prone to local microstructure deformation. In fact, there is a difference in surface machining accuracy among the rolling bodies. As these rolling bodies interact with raceway and cage, the contact stresses between the interaction surfaces would be diverse from each other, which leads to different degrees of local microstructure deformation on their surfaces in the frequent collision with raceway and cage. At the same time, due to the autorotation of the rolling body, microscopic shape on its surface in contact with raceway or cage will affect the amplitude of the impact impulse. The simultaneous action of these complex factors will make the running health condition of Bearing 4-Ch 8 fluctuate slightly in the steady-state stage. Even so, compared with the other two methods, $α (q)$ of the proposed method has smaller variation and is more suitable for bearing condition monitoring.

Figure 20.

Trend analysis for stage I of singularity exponent $α (q)$ at $q = 0$ for Bearing 3-Ch 5 with (a) MFDFA1, (b) EMD-MFDFA.

Figure 21.

Trend analysis for stage I of singularity exponent $α (q)$ at $q = 1$ for Bearing 4-Ch 8 with (a) MFDFA1, (b) EMD-MFDFA.

It can be found that both multifractal spectrums in Figure 11 and Figure 15 present anomalies at both ends. It indicates that both the proposed method and MFDFA1 seem to destroy the multifractal structure of the original signal. However, according to the previous analysis results, it can be known that for the partial multifractal spectrums at positive q, the corresponding multifractal spectrum characteristic parameters $α (q)$ and $f (α)$ exhibit increasingly violent oscillations when q gradually deviates from 0. Consequently, the variance $σ_{1}$ of the multifractal characteristic parameters obtained by the proposed method, MFDFA1 and EMD-MFDFA are used to quantify their respective volatility, which are shown in Figure 22. It can be seen from Figure 22 that the q corresponding to the least volatile $α (q)$ and $f (α)$ extracted by these three methods lies in [−2, 2]. On the contrary, the trend of $σ_{1}$ indicates that the multifractal characteristic parameters at both ends of the multifractal spectrum are more volatile. The reason leading to much volatility is that they may be more susceptible to the effects of background noise and slight variation of operating conditions. Therefore, it is acceptable to sacrifice the structure in this part of the multifractal spectrum to enhance the performance of the characteristic parameters at other locations for distinguishing the different health conditions of the bearing.

Figure 22.

Fluctuation evaluation of trend curves of various multifractal spectral characteristic parameters calculated by proposed method, MFDFA1 and EMD-MFDFA: (a) $α (q)$ , (b) $f (α)$ .

Conclusion

In this paper, an approach of MEMD-MFDFA based on the IMF selection strategy and SMW segmentation method is proposed to extract multifractal features for the early fault identification and performance degradation assessment of rolling bearings. Several simulated and experimental cases are given to verify the correctness of the proposed method. The following conclusions could be drawn.

(1) The developed MEMD-MFDFA method is superior to MFDFA and EMD-MFDFA in anti-noise by investigating the simulated multifractal signal produced from the p-model algorithm.

(2) Multifractal characteristic parameters on the multifractal spectrum at positive q of vibration signals can be employed to reflect the performance degradation trend of running bearings.

(3) With the increase of q, the fluctuation of the corresponding multifractal characteristic parameters becomes stronger. Moreover, the appropriate characteristic parameters of the multifractal spectrum in the middle region can be determined according to their fluctuation degree, which provides a new idea for the establishment of an effective performance degradation indicator based on the multifractal spectrum.

(4) Compared with MFDFA and EMD-MFDFA, the multifractal characteristic parameters extracted by MEMD-MFDFA are proved to be more effective in early fault identification and performance degradation assessment according to the experimental results.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (Grant no. 51765034) and Excellent Graduate Student “Innovation Star“ project of Education Department of Gansu Province (Project no. 2021CXZX-446). This work was also supported by Science and Technology Projects of Gansu Province (Project no. 21JR7RA305) and Key Laboratory of Cloud Computing of Gansu Province (Gansu Computing Center).

ORCID iD

Changfeng Yan

References

Lee

Zhao

, et al. Prognostics and health management design for rotary machinery systems-reviews, methodology and applications. Mech Syst Sig Process 2014; 42(1-2): 314–334.

Wang

Qiao

, et al. Fractional envelope analysis for rolling element bearing weak fault feature extraction. IEEE/CAA J Autom Sin 2017; 4(2): 353–360.

Cerrada

Sanchez

R-V

, et al. A review on data-driven fault severity assessment in rolling bearings. Mech Syst Sig Process 2018; 99: 169–196.

Miao

Azarian

, et al. Health monitoring of cooling fan bearings based on wavelet filter. Mech Syst Sig Process 2015; 64-65: 149–161.

Peeters

Antoni

Helsen

. Blind filters based on envelope spectrum sparsity indicators for bearing and gear vibration-based condition monitoring. Mech Syst Sig Process 2020; 138: 106556.

Machado de Azevedo

Araujo

Bouchonneau

. A review of wind turbine bearing condition monitoring: State of the art and challenges. Renew Sust Energ Rev 2016; 56: 368–379.

Hou

Wang

Xia

, et al. Investigations on quasi-arithmetic means for machine condition monitoring. Mech Syst Sig Process 2021; 151: 107451.

Soualhi

Medjaher

Zerhouni

. Bearing health monitoring based on hilbert-huang transform, support vector machine, and regression. IEEE Trans Instrum Meas 2015; 64(1): 52–62.

Meng

Yan

Chen

, et al. Health indicator of bearing constructed by RMS-CUMSUM and GRRMD-CUMSUM with multi-features of envelope spectrum. IEEE Trans Instrum Meas 2021; 70: 1–16.

10.

Feng

Liang

Chu

. Recent advances in time-frequency analysis methods for machinery fault diagnosis: A review with application examples. Mech Syst Sig Process 2013; 38(1): 165–205.

11.

Nassser

Mohammadi

. Condition monitoring using wavelet transform and fuzzy logic by vibration signals. Life Sci J 2012; 9(4): 5680–5685.

12.

Miao

. Deep learning based approach for bearing fault diagnosis. IEEE Trans Ind Appl 2017; 53(3): 3057–3065.

13.

Tong

Cao

Han

, et al. A fault diagnosis approach for rolling element bearings based on RSGWPT-LCD bilayer screening and extreme learning machine. IEEE Access . 2017; 5: 5515–5530.

14.

Huang

Shen

Long

, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc Roy Soc Lond A Math Phys Eng Sci 1998; 454: 903–995.

15.

de Moura

Souto

Silva

, et al. Evaluation of principal component analysis and neural network performance for bearing fault diagnosis from vibration signal processed by RS and DF analyses. Mech Syst Sig Process 2011; 25(5): 1765–1772.

16.

Jiang

Wang

. Study on nature of crossover phenomena with application to gearbox fault diagnosis. Mech Syst Sig Process 2017; 83: 272–295.

17.

Lin

Chen

. A novel method for feature extraction using crossover characteristics of nonlinear data and its application to fault diagnosis of rotary machinery. Mech Syst Sig Process 2014; 48(1–2): 174–187.

18.

Liu

Zhang

Cheng

, et al. Fault diagnosis of gearbox using empirical mode decomposition and multi-fractal detrended cross-correlation analysis. J Sound Vib 2016; 385: 350–371.

19.

Kantelhardt

Zschiegner

Koscielny-Bunde

, et al. Multifractal detrended fluctuation analysis of nonstationary time series. Phys A Stat Mech Appl 2002; 316(1–4): 87–114.

20.

Lin

Chen

. Fault diagnosis of rolling bearings based on multifractal detrended fluctuation analysis and Mahalanobis distance criterion. Mech Syst Sig Process 2013; 38(2): 515–533.

21.

Xiong

Zhang

, et al. A fault diagnosis method for rolling bearings based on feature fusion of multifractal detrended fluctuation analysis and alpha stable distribution. Shock Vib . 2016; 2016: 1–12.

22.

Liu

Wang

. Rolling bearing fault diagnosis based on LCD-TEO and multifractal detrended fluctuation analysis. Mech Syst Sig Process 2015; 60-61: 273–288.

23.

Yuan

Peng

, et al. Rolling bearing fault diagnosis based on adaptive smooth ITD and MF-DFA method. J Low Freq Noise V A 2020; 39(4): 968–986.

24.

Liu

Wang

, et al. Feature extraction method based on VMD and MFDFA for fault diagnosis of reciprocating compressor valve. J Vibroeng 2017; 19(8): 6007–6020.

25.

Zhang

Liang

. Dynamic degradation observer for bearing fault by MTS-SOM system. Mech Syst Sig Process 2013; 36(2): 385–400.

26.

Chen

Yan

Meng

, et al. A novel two-stage feature extraction based on multi-fractal spectrum with CI-VMD for early fault diagnosis of bearing. In: Proceedings of IEEE International Conference Prognostics and System Health Management, Jinan, China, 23–25 October, 2020, 2020, pp. 261–269.

27.

Ihlen

EAF

. Multifractal analyses of response time series: a comparative study. Beh Res Meth 2013; 45(4): 928–945.

28.

Kang

Pecht

. Fault diagnosis using adaptive multifractal detrended fluctuation analysis. IEEE Trans Ind Electron 2020; 67(3): 2272–2282.

29.

G-F

Zhou

W-X

. Detrending moving average algorithm for multifractals. Phys Rev E 2010; 82(1): 011136.

30.

Schumann

Kantelhardt

. Multifractal moving average analysis and test of multifractal model with tuned correlations. Phys A Stat Mech Appl 2011; 390(14): 2637–2654.

31.

Guo

Gong

, et al. Multifractal characterization of mechanical vibration signals through improved empirical mode decomposition-based detrended fluctuation analysis. Proc Inst Mech Eng C: J Mech Eng Sci 2017; 231(22): 4139–4149.

32.

Mohanty

Gupta

Raju

. Hurst based vibro-acoustic feature extraction of bearing using EMD and VMD. Measurement . 2018; 117: 200–220.

33.

Mert

Akan

. Detrended fluctuation thresholding for empirical mode decomposition based denoising. Digit Signal Process 2014; 32: 48–56.

34.

Peng

Buldyrev

Havlin

, et al. Mosaic organization of DNA nucleotides. Phys Rev E 1994; 49(2): 1685–1689.

35.

Dong

Jia

, et al. Non-iterative denoising algorithm for mechanical vibration signal using spectral graph wavelet transform and detrended fluctuation analysis. Mech Syst Sig Process 2021; 149.

36.

Rilling

Flandrin

. One or two frequencies: the empirical mode decomposition answers. IEEE Trans Signal Process 2007; 56(1): 85–95.

37.

Boashash

. Estimating and interpreting the instantaneous frequency of a signal. II. Algorithms and applications. IEEE Proc . 1992; 80(4): 540–568.

38.

Zhou

Leung

. Multifractal temporally weighted detrended fluctuation analysis and its application in the analysis of scaling behavior in temperature series. J Stat Mech-Theory E 2010; 2010(06): P06021.

39.

Que

Jin

. Remaining useful life prediction for bearings based on a gated recurrent unit. IEEE Trans Instrum Meas 2021; 70: 3511411.

40.

Meneveau

Sreenivasan

. Simple multifractal cascade model for fully developed turbulence. Phys Rev Lett 1987; 59(13): 1424.

41.

Randall

Sawalhi

Coats

. A comparison of methods for separation of deterministic and random signals. Int J Cond Monit 2011; 1(1): 11–19.

42.

Peeters

Guillaume

Helsen

. A comparison of cepstral editing methods as signal pre-processing techniques for vibration-based bearing fault detection. Mech Syst Sig Process 2017; 91: 354–381.

43.

Kantelhardt

Koscielny-Bunde

Rego

, et al. Detecting long-range correlations with detrended fluctuation analysis. Phys A Stat Mech Appl 2001; 295(3): 441–454.

44.

Lee

Qiu

, et al. Rexnord Technical Services: Bearing Data Set [Online]. 2007. Available at: http://ti.arc.nasa.gov/project/prognostic-data-repository (Accessed 13 November 2020).

45.

El-Thalji

Jantunen

. A descriptive model of wear evolution in rolling bearings. Eng Fail Anal 2014; 45: 204–224.

46.

Ren

, Remaining useful life estimation of rolling bearings based on sparse representation. In: Proceeding of 7th International Conference Mechanical Aerospace Engineering, London, UK, 18–20 July, 2016, 2016, pp. 209–213.

47.

Arakere

Subhash

. Work hardening response of M50-NiL case hardened bearing steel during shakedown in rolling contact fatigue. Mater Sci Technol 2012; 28(1): 34–38.