Multifractal manifold for rotating machinery fault diagnosis based on detrended fluctuation analysis

2.1. MFDFA method

In order to overcome the disadvantages of DFA and reveal the characteristics of signal comprehensively, multifractal theory is introduced to DFA, which can represent the local fractal feature of signal in different scales. The MFDFA method can be expressed as the following 5 steps [15]:

(1) For a time series $x_j$, $j=\mathrm{}$1, 2,…, $N$, determine a cumulative series:

(2) Divide the cumulative series $Y\left(i\right)$ into $N_s$ nonoverlapping segments with the same length$s$, where $N_s\equiv \mathrm{i}\mathrm{n}\mathrm{t}(N/s)$. In many cases, the series length $N$ is not a multiple of the segment length$s$, a small part of series might remain at the end. Then, the segmentation procedure is implemented again from the opposite end of the series to make full use of the remaining series. Accordingly, all of $2N_s$ segments are obtained.

(3) Implement the least-square algorithm to calculate the local trend $y_v$ of each segment:

where $m$ is the order of polynomial, $a_k$ are polynomial coefficients and $t$ is the time of data point in the $v$-th segment. Then the variance can be expressed as:

when $v=1,\dots ,N_s$, and:

when $v=N_s+1,\dots ,2N_s$. $F(v,s)$ is also known as the local root mean square (RMS) error of the $v$th segment.

(4) For all segments, the $q$-th order fluctuation function can be defined as:

When $q=\text{2}$, the multifractal DFA degenerates to the standard DFA. Repeat step (2)-(4) with different segment length $s$ and fluctuation function order $q$, the trend of each segment in different scales can be eliminated and the local fractal characteristics can be revealed.

(5) The logarithmic relationships between $F_q\left(s\right)$ and $s$ are analyzed with different order $q$ to determine the scaling behavior of the fluctuation functions. If the series $x_j$ is long-range power-law correlated [19], $F_q\left(s\right)$ increases along with $s$ following a power-law:

When $q=\text{2}$, $h\left(2\right)$ is known as the Hurst exponent [15]. Accordingly, $h\left(q\right)$ is the generalized Hurst exponent. When $q=0$, $F_q\left(s\right)$ is expressed by a logarithmic averaging procedure:

With different order $q$, the fluctuation functions of signals can be represented more comprehensively and the differences of different signals can be highlighted in different scales.

2.2. Fluctuation function order and polynomial order

In normal circumstances, mono-fractal series is very rare while multifractal series are widespread, and the vibration signals of rotating machinery in fault condition show significant multifractal characteristics.

Fig. 1Signals and corresponding q-order RMS

As shown in Fig. 1, the $q$-order local root mean square (qRMS, i.e., $\left[F^2\right(v,s){]}^{q/2}$) of fault bearing signal and normal bearing signal are calculated. With positive $q$, $\left[F^2\right(v,s){]}^{q/2}$ are sensitive to the transient impulses in the signal (black dashed line). In contrast, with negative $q$, $\left[F^2\right(v,s){]}^{q/2}$ are sensitive to local fluctuations with small magnitudes between adjacent impulses (purple dashed line). With $q=\text{2}$ and –2, $\left[F^2\right(v,s){]}^{q/2}$ are more sensitive to impulses and small local fluctuations than those with $q=\mathrm{}$1 and –1. When $q=0$, $\left[F^2\right(v,s){]}^{q/2}$ is neutral to both of impulses and small local fluctuations.

The mono-fractal DFA method [12-14] ($q=\text{2}$) can only detect the transient impulses in a specific scale, which is insufficient to reveal the characteristics of the signals. In contrast, MFDFA can embody fractal characteristics of impulses and small local fluctuations in different scales, which can represent the operating status of rotating machinery more comprehensively.

For non-stationary rotating machinery vibration signals, there are significant nonlinear characteristics, especially in the range of transient impulse waveform. The local fluctuations depend on the fitting result of local trend. The fitting effects of fault bearing signal with different polynomial order $m$ are shown in Fig. 2. The cumulative series are divided into 3 regions according to waveform. The mean of local RMS in each region are shown in Table 1.

Fig. 2Local RMS with different polynomial order

When $m=$ 1, the local trend are linear and the means of RMS in transient impulse regions 1 and 3 are much larger than those with higher polynomial orders. With the polynomial order increasing, the means of RMS are reduced significantly. When $m$ increases to 5, the RMS reduction becomes small.

Transient impulses are widespread in rotating machinery fault signals, similar but different fault patterns might contain similar transient impulses. Therefore, when linear local trend ($m=\text{1}$) is powerless to highlight the subtle differences among similar fault patterns, a higher polynomial order is necessary.

3.1. Multifractal feature based on DFA

In multifractal spectrum [15], the generalized Hurst exponent $h\left(q\right)$ are corresponding to the slopes of linear fitting line of $F_q\left(s\right)$-$s$ logarithmic curves with different order $q$. Since the non-stationary characteristics of rotating machinery vibration signals, the logarithmic relationships between $F_q\left(s\right)$ and $s$ are nonlinear and cannot be represented by the slopes of overall linear fitting line. Thereby, Moura [12, 13] selected $F\left(s\right)$ with different $s$ as feature vector for each vibration signal and reduce the dimension of feature vectors with PCA to recognize signal patterns; Lin [14] used segmented slopes of logarithmic curves according to crossover characteristics to recognize signal patterns. The 2 methods are both based on mono-fractal DFA ($q=\text{2}$) and the features lack some nonlinear characteristics in varying degrees.

In order to represent the nonlinear characteristics of $F_q\left(s\right)$-$s$ logarithmic curves, the multifractal feature based on DFA can be defined as:

$h(q,i)$ are the slopes of adjacent $F_q\left(s\right)$ and $s$ in logarithmic curves. Where the segment length $s=[s_1,s_2,\dots ,s_K]$, $q=\mathrm{}$[–2, –1, 0, 1, 2] and $i=\mathrm{}$1, 2,…, $K-1$. Thereby, a $h(q,i)$ matrix is obtained. The $h(q,i)$ matrix can be converted to a multifractal feature vector $mf$ ($\in R^D$) by connecting each row:

3.2. Manifold learning

The multifractal feature vectors $mf$are high-dimensional features with nonlinear characteristics. In order to recognize the patterns of high-dimensional nonlinear features, manifold learning [16, 17] is applied. Manifold learning is a global nonlinear dimensionality reduction method, which is mainly used in pattern and image recognition.

LTSA [18] is an effective manifold learning technique for nonlinear dimensionality reduction. For $L$ samples from different signal patterns, a multifractal feature matrix $MF$ can be expressed as:

The feature dimension $D$ of multifractal feature matrix can be reduced to $d$ ($d<D$) by LTSA as following steps [18]:

(1) Extracting local information.

Determine the $k$ nearest neighborhoods for each $mf$ vector through the Euclidean distance:

Centralize the $m{f}_i^k$ matrix by removing the mean of corresponding $k$nearest neighborhoods and then calculate the $d$ largest right singular vectors of the centralized $m{f}_i^k$ matrix to obtain the orthogonal basis $V_i$:

(2) Constructing alignment matrix.

Determine the 0-1 selection matrix $S_i$:

Calculate the correlation matrix $W_i$:

where $e_k$ is a column vector of $k$ ones.

Construct the alignment matrix $B$:

(3) Aligning global coordinates.

Calculate the $d+1$ smallest eigenvectors of matrix $B$, and then the 2nd to the $d+1$th smallest eigenvalues are selected to constitute the global coordinates$MFM$ ($\in R^{d\times L}$). Thence, the $d$-dimension multifractal manifold features $MFM$ can be expressed as:

where $u_{i+1}\in R^L$ ($i=\mathrm{}$1, 2,…, $d$) is corresponding to the $i$th dimensionality reduction results.

The low-dimensional multifractal manifold features can reveal the differences of different signal patterns which are beneficial for pattern recognition and fault diagnosis.

Combining the advantages of MFDFA in time series nonlinear characteristics extraction and manifold learning in nonlinear dimensionality reduction, a novel multifractal manifold method is proposed. The proposed MFM method contains 3 steps:

(1) Calculate the fluctuation functions of series according to MFDFA with an appropriate polynomial order.

(2) Construct nonlinear multifractal feature vector $mf$ for each signal series according to Eqs. (9-10).

(3) Reduce the dimension of multifractal feature vectors through manifold learning for pattern recognition.

4.1. Case I: Gearbox fault

The signal samples are acquired from a single cylindrical gear reducer. The signal patterns are normal, driven gear pitting, driving gear wear and mixed fault (driving gear wear and driven gear pitting). The sampling frequency is 5120 Hz. Each sample contains 5120 sampling points. The parameters of gears are shown in Table 2.

The signal of each pattern is shown in Fig. 3(a). The multifractal $F\left(s\right)$-$s$ logarithmic curves and multifractal features $h(q,i)$ are shown in Fig. 3(b) ($m=\text{3}$).

Fig. 3a) Gearbox signals, b) multifractal logarithmic curves

a)

b)

The clustering results of MFM and MF-PC with different polynomial order $m$ are shown in Fig. 4. The first 3 low-dimensional features of MFM ($m=\text{3}$) reveal the differences of signal samples, the samples from each pattern have the same features and all samples are distributed at 4 points, as shown in Fig. 4(a), (b); when $m=\text{1}$, the samples from wear and mixed fault patterns are not completely distinguished, as shown in Fig. 4(c), (d). As shown in Fig. 4(e), (f), the signal samples are distinguished from each other by the first 3 low-dimensional features of MF-PC. However, the distribution variance of simples with $m=\text{3}$ is smaller than that with $m=\text{1}$.

The $F\left(s\right)$-$s$ logarithmic curves based on mono-fractal DFA are calculated. For each signal pattern, the curves of 5 samples are plotted in Fig. 5. The curves are divided into 3 regions according to the crossover points and the slopes of fitting lines in each region are H1, H2 and H3.

The clustering results of crossover characteristics and DFA-PC are shown in Fig. 6. Most signal samples are distinguished from each other; only several samples are difficult to identify.

The classification results of SVM are shown in Table 3. All samples from different patterns are classified by 4 methods with high accuracies, while only MF-PC method and the MFM method with $m=$ 3 classify samples completely correct.

Fig. 4Process results

a) MFM features with $m=$ 3

b) MFM with $m=$ 3

c) MFM features with $m=$ 1

d) MFM with $m=$ 1

e) MF-PC with $m=\text{3}$

f) MF-PC with $m=\text{1}$

Fig. 5Mono-fractal logarithmic curves of gearbox signals

a)$m=\text{3}$

b)$m=\text{1}$

Fig. 6Process results

a) Crossover characteristics with $m=$ 3

b) Crossover characteristics with $m=$ 1

c) DFA-PC with $m=$ 3

d) DFA-PC with $m=$ 1

4.2. Case II: bearing rolling element with multi defect severities

Bearing data with different defect severities is obtained from the CWRU bearing data center [22]. All signal samples are acquired from normal bearing (type SKF-6205) and bearings with defects on the surface of rolling element and the defect sizes are 0.007, 0.014, 0.021 and 0.028 inch. The rotating speed of drive shaft is 1730 revolutions per minute (RPM) and the sampling frequency is 12000 Hz. Each sample contains 4800 sampling points.

The signal of each pattern is shown in Fig. 7(a). The multifractal $F\left(s\right)-s$ logarithmic curves and multifractal features $h(q,i)$ are shown in Fig. 7(b) ($m=$ 4).

Fig. 7Signals of bearings

a) Signals

b) Multifractal logarithmic curves

The clustering results of MFM and MF-PC are shown in Fig. 8. 4 MFM features ($m=$ 4) reveal the differences of signal samples; in order to represent these features in 3-dimension, 4 MFM features are reduced to 3 by PCA and the corresponding features are MFM-PCs, the samples from 5 patterns are distributed at 5 points, as shown in Fig. 8(a), (b); when $m=\text{1}$, the samples from fault patterns of defect sizes 0.007 and 0.021 inch are recognized as a same pattern incorrectly, as shown in Fig. 8(c), (d). Simultaneously, the samples from fault patterns of defect sizes 0.007 and 0.021 inch patterns are still difficult to be distinguished by MF-PC method, as shown in Fig. 8(e), (f).

The $F\left(s\right)$-$s$ logarithmic curves based on mono-fractal DFA are calculated. For each signal pattern, the curves of 5 samples are plotted in Fig. 9. The curves are divided into 3 regions according to the crossover points and the slopes of fitting lines in each region are H1, H2 and H3.

The clustering results of crossover characteristics and DFA-PC are shown in Fig. 10. 5 signal patterns are distinguished from each other with a higher polynomial order $m=\text{4}$, as shown in Fig. 10(a), (c); while the samples from fault patterns of defect sizes 0.007 and 0.021 inch are unrecognizable when $m=\text{1}$, as shown in Fig. 10(b), (d).

The classification results of SVM are shown in Table 4. When $m=\text{1}$, there are different quantities of misclassified samples in the result of each method and the error samples are mainly from fault patterns of defect sizes 0.007 and 0.021 inch. With a higher polynomial order $m=\text{4}$, the classification accuracies of all 4 methods are significantly improved.

Fig. 8Process results

a) MFM features with $m=$ 4

b) MFM with $m=$ 4

c) MFM features with $m=$ 1

d) MFM with $m=$ 1

e) MF-PC with $m=\text{4}$

f) MF-PC with $m=\text{1}$

Fig. 9Mono-fractal logarithmic curves of bearing signals

a)$m=$ 4

b)$m=\text{1}$

Fig. 10Process results

a) Crossover characteristics with $m=$ 4

b) Crossover characteristics with $m=$ 1

c) DFA-PC with $m=$ 4

d) DFA-PC with $m=$ 1

4.3. Case III: bearing consecutive extended defect

Self-aligning ball bearings (type HRB-1209) with different ranges of pitting defects and normal bearing are installed to the test end and tested individually, as shown in Fig. 11. The depths of artificially machined pitting defects are 0-0.5 mm. The rotating frequency of motor is 15 Hz and the signal sampling frequency is 25600 Hz. Each sample contains 12800 sampling points. The size ranges of bearing defects are shown in Table 5.

Fig. 11Experimental apparatus

a) Experimental platform

b) Bearings with defects

The signal of each pattern is shown in Fig. 12(a). The multifractal $F\left(s\right)$-$s$ logarithmic curves and multifractal features $h(q,i)$ are shown in Fig. 12(b) ($m=$ 3).

Fig. 12Signals of bearings

a) Signals

b) Multifractal logarithmic curves

The clustering results of MFM and MF-PC are shown in Fig. 13. The first 2 low-dimensional features of MFM ($m=$ 3) reveal the differences of signal samples, the samples from each pattern have the same features and all samples are distributed at 3 points, as shown in Fig. 13(a), (b); when $m=$ 1, the samples from local defect and extended defect fault patterns are not completely distinguished, as shown in Fig. 13(c), (d). The first 2 low-dimensional features of MF-PC are applied to identify the signal samples; the mixing phenomenon of samples is obvious, as shown in Fig. 13(e), (f).

The $F\left(s\right)-s$logarithmic curves based on mono-fractal DFA are calculated. For each signal pattern, the curves of 5 samples are plotted in Fig. 14. The curves are divided into 2 regions according to the crossover point and the slopes of fitting lines in each region are H1 and H2.

The curves from 3 patterns are very similar in both of $F\left(s\right)$ and segmented slopes. Therefore, the clustering results of the 2 methods based on DFA are unsatisfactory, as shown in Fig. 15.

The classification results of SVM are shown in Table 6. With a higher polynomial order, the classification accuracies are improved. The signal samples are correctly classified by the proposed MFM method; however, the signal samples are still hardly identified by the 2 DFA-based methods and MF-PC method.

Fig. 13Process results

a) MFM features with $m=$ 3

b) MFM with $m=$ 3

c) MFM features with $m=$ 1

d) MFM with $m=$ 1

e) MF-PC with $m=$ 3

f) MF-PC with $m=\text{1}$

Fig. 14Mono-fractal logarithmic curves of bearing signals

a)$m=\text{3}$

b)$m=\text{1}$

Fig. 15Process results

a) Crossover characteristics with $m=$ 3

b) Crossover characteristics with $m=$ 1

c) DFA-PC with $m=$ 3

d) DFA-PC with $m=$ 1

4.4. Case IV: bearing discrete extended defect

Cylindrical roller bearings (type TMB-N209M) with different number of defects and normal bearing are also tested on the experiment platform in Fig. 11(a). The sizes of artificially machined defects are 0.5 mm in depth and 0.6 mm in width. The rotating frequency of motor is 15 Hz and the signal sampling frequency is 25600 Hz. Each sample contains 12800 sampling points. The amount and feature of bearing defects are shown in Table 7 and Fig. 16.

Fig. 16Bearings with discrete defects

The signal of each pattern is shown in Fig. 17(a). The multifractal $F\left(s\right)$-$s$ logarithmic curves and multifractal features $h(q,i)$ are shown in Fig. 17(b) ($m=$ 4).

Fig. 17Signals of bearings

a) Signals

b) Multifractal logarithmic curves

The clustering results of MFM and MF-PC are shown in Fig. 18. The first 3 low-dimensional features of MFM ($m=$ 4) reveal the differences of signal samples, the samples from each pattern have the same features and all samples are distributed at 4 points, as shown in Fig. 18(a), (b); when $m=$ 1, the samples from 3 fault patterns are not distinguished from each other accurately, as shown in Fig. 18(c), (d). The clustering results of MF-PC with different polynomial orders are similar and the samples from normal and single defect patterns are mixed, as shown in Fig. 18(e), (f).

The $F\left(s\right)$-$s$ logarithmic curves based on mono-fractal DFA are calculated. For each signal pattern, the curves of 5 samples are plotted in Fig. 19. The curves are divided into 3 regions according to the crossover points and the slopes of fitting lines in each region are H1, H2 and H3.

The segmented slopes of curves from 3 fault patterns are similar, and $F\left(s\right)$ of samples from normal and single defect patterns are similar. Therefore, some samples from different patterns are mixed in the clustering results of the 2 DFA-based methods, as shown in Fig. 20.

The classification results of SVM are shown in Table 8. With a higher polynomial order, the classification accuracies are improved. The signal samples are correctly classified by the proposed MFM method; however, the misclassified samples still exist in the results of the 2 DFA-based methods and MF-PC method.

Fig. 18Process results

a) MFM features with $m=$ 4

b) MFM with $m=$ 4

c) MFM features with $m=$ 1

d) MFM with $m=$ 1

e) MF-PC with $m=$ 4

f) MF-PC with $m=\text{1}$

Fig. 19Mono-fractal logarithmic curves of bearing signals

a)$m=\text{4}$

b)$m=\text{1}$

Fig. 20Process results

a) Crossover characteristics with $m=$ 3

b) Crossover characteristics with $m=$ 1

c) DFA-PC with $m=$ 3

d) DFA-PC with $m=$ 1