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

Abstract. The vibration signals of rotating machinery in fault conditions are non-stationary and nonlinear. For the non-stationary and nonlinear characteristics of fault vibration signals, a novel multifractal manifold (MFM) method based on detrended fluctuation analysis (DFA) is proposed. The proposed method consists of three steps. Firstly, calculate the multifractal fluctuation functions of signal series with an appropriate polynomial order, according to multifractal DFA method. Secondly, construct multifractal feature vector for each signal sample to reveal the nonlinear characteristics in different scales. Finally, implement manifold learning to reduce the dimension of multifractal feature vectors. The obtained low-dimensional MFM features can reveal the differences of signal samples from different fault patterns effectively, which are benefit for automatic pattern recognition and multiple fault diagnosis. The recognition performance of the proposed MFM method is verified by fault experiments of gearbox and rolling element bearing, which demonstrates the superiority of MFM method in rotating machinery fault diagnosis compared to other DFA-based methods.


Introduction
Rotating machinery is extensively applied in industrial production and usually runs in harsh environments.Therefore, in order to ensure the security and stability of industrial production, the corresponding fault diagnosis for rotating machinery are very essential.In fault conditions, the vibration signals of rotating machinery are non-stationary with nonlinear characteristics.In recent years, time-frequency analysis methods are applied to non-stationary signal processing such as wavelet transform [1,2], Wigner-Ville distribution [3,4], spectral kurtosis [5,6], empirical mode decomposition [7,8] and EMD derived methods [9].However, all these methods are based on signal decomposition or filtering and a comprehensive nonlinear characteristic of signal might be dispersed in several components, which make feature extraction difficult.
Detrended fluctuation analysis (DFA) [10,11] is proposed to reveal the long-range correlations of non-stationary and nonlinear time series and applied in rotating machinery fault diagnosis.Moura [12,13] applies fluctuation function vectors to represent signal features and reduces the vector dimension with linear technique PCA for gearbox fault diagnosis.However, the fluctuation function vectors are incomplete features, and PCA method neglects the nonlinear characteristics of fluctuation function vectors.Lin [14] applies the crossover characteristics of logarithmic curves based on DFA to represent the characteristics of signals for rotating machinery fault diagnosis.While the location and quantity of crossover points of different fault signals are different and they must be determined in advance, which is difficult to implement for automatic multiple fault pattern recognition.Moreover, the features segmented by crossover points still neglect the nonlinear characteristics of logarithmic curves in the vicinity of crossover points.
Further studies indicate that DFA method is only sensitive to the transient impulses of rotating machinery vibration signals.In order to comprehensively reveal the characteristics of signal, multifractal detrended fluctuation analysis (MFDFA) [15] is employed.MFDFA can represent the fluctuation status of time series in different scales.In order to overcome the disadvantages of nonlinear characteristic lack, the local consecutive slopes of MFDFA logarithmic curves are determined as multifractal features.For high-dimensional and nonlinear multifractal features, a global nonlinear dimensionality reduction method manifold learning [16][17][18] is applied to reduce the dimension of multifractal features and recognize signal patterns.
In this study, a multifractal manifold (MFM) method based on DFA is proposed.The performance of the proposed MFM method is verified by experimental data from fault rotating machinery components and compared with related DFA-based methods [12][13][14].The effects of different polynomial orders on diagnostic results are also discussed which are ignored in literatures [12][13][14].

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 , = 1, 2,…, , determine a cumulative series: (2) Divide the cumulative series ( ) into nonoverlapping segments with the same length , where ≡ int( / ).In many cases, the series length is not a multiple of the segment length , 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 2 segments are obtained.
(3) Implement the least-square algorithm to calculate the local trend of each segment: where is the order of polynomial, are polynomial coefficients and is the time of data point in the -th segment.Then the variance can be expressed as: when = 1, … , , and: when = + 1, … ,2 .( , ) is also known as the local root mean square (RMS) error of the th segment.
(4) For all segments, the -th order fluctuation function can be defined as: When = 2, the multifractal DFA degenerates to the standard DFA.Repeat step ( 2)-( 4) with different segment length and fluctuation function order , the trend of each segment in different scales can be eliminated and the local fractal characteristics can be revealed.
(5) The logarithmic relationships between ( ) and are analyzed with different order to determine the scaling behavior of the fluctuation functions.If the series is long-range power-law correlated [19], ( ) increases along with following a power-law: When = 2, ℎ( 2) is known as the Hurst exponent [15].Accordingly, ℎ( ) is the generalized Hurst exponent.When = 0, ( ) is expressed by a logarithmic averaging procedure: With different order , the fluctuation functions of signals can be represented more comprehensively and the differences of different signals can be highlighted in different scales.

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.As shown in Fig. 1, the -order local root mean square (qRMS, i.e., [ ( , )] ⁄ ) of fault bearing signal and normal bearing signal are calculated.With positive , [ ( , )] ⁄ are sensitive to the transient impulses in the signal (black dashed line).In contrast, with negative , [ ( , )] ⁄ are sensitive to local fluctuations with small magnitudes between adjacent impulses (purple dashed line).With = 2 and -2, [ ( , )] ⁄ are more sensitive to impulses and small local fluctuations than those with = 1 and -1.When = 0, [ ( , )] ⁄ is neutral to both of impulses and small local fluctuations.
The mono-fractal DFA method [12][13][14]  In multifractal spectrum [15], the generalized Hurst exponent ℎ( ) are corresponding to the slopes of linear fitting line of ( ) -logarithmic curves with different order .Since the non-stationary characteristics of rotating machinery vibration signals, the logarithmic relationships between ( ) and are nonlinear and cannot be represented by the slopes of overall linear fitting line.Thereby, Moura [12,13] selected ( ) with different 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 ( = 2) and the features lack some nonlinear characteristics in varying degrees.

Manifold learning
The multifractal feature vectors 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 samples from different signal patterns, a multifractal feature matrix can be expressed as: The feature dimension of multifractal feature matrix can be reduced to ( < ) by LTSA as following steps [18]: (1) Extracting local information.Determine the nearest neighborhoods for each vector through the Euclidean distance: Centralize the matrix by removing the mean of corresponding nearest neighborhoods and then calculate the largest right singular vectors of the centralized matrix to obtain the orthogonal basis : (2) Constructing alignment matrix.
Determine the 0-1 selection matrix : Calculate the correlation matrix : where is a column vector of ones.Construct the alignment matrix : (3) Aligning global coordinates.
Calculate the + 1 smallest eigenvectors of matrix , and then the 2nd to the + 1th smallest eigenvalues are selected to constitute the global coordinates (∈ × ).Thence, the -dimension multifractal manifold features can be expressed as: where ) is corresponding to the 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.
(3) Reduce the dimension of multifractal feature vectors through manifold learning for pattern recognition.

Experimental verification
The proposed MFM method is tested by 4 cases of rotating machinery fault.For comparison, several related method are also tested.The corresponding methods are shown as follows: (1) MFM: the proposed method, the low-dimensional features are MFMs.
(2) Crossover characteristics [14]: the segmented slopes of ( )logarithmic curve based on mono-fractal DFA are applied to represent the characteristics of signal series, which are H1, H2 and H3 (polynomial order = 1).
(3) DFA-PC [12,13]: the mono-fractal fluctuation functions ( ) with different are selected as feature vectors to represent the characteristics of signal series, the dimension of vectors are reduced by PCA.The low-dimensional features are DFA-PCs (polynomial order = 1).
(4) MF-PC: based on DFA-PC method, the multifractal fluctuation functions are selected as feature vectors and the dimension of vectors are also reduced by PCA.The low-dimensional features are MF-PCs.
To test the performance of each method, clustering results of signal samples are obtained.50 samples from each pattern are divided into two parts randomly: 20 train samples and 30 test samples, and then test samples are classified by support vector machine (SVM) [20,21].

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 The clustering results of MFM and MF-PC with different polynomial order are shown in Fig. 4. The first 3 low-dimensional features of MFM ( = 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 = 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 = 3 is smaller than that with = 1.
The ( )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 = 3 classify samples completely correct.

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 ( ) − logarithmic curves and multifractal features ℎ( , ) are shown in Fig. 7(b) ( = 4).The clustering results of MFM and MF-PC are shown in Fig. 8. 4 MFM features ( = 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 = 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 ( )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 = 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 = 1, as shown in Fig. 10(b), (d).
The classification results of SVM are shown in Table 4.When = 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 = 4, the classification accuracies of all 4 methods are significantly improved.

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.The clustering results of MFM and MF-PC are shown in Fig. 13.The first 2 low-dimensional features of MFM ( = 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 = 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 ( ) − 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 ( ) 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.

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.The signal of each pattern is shown in Fig. 17 The ( )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 ( ) 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.

Discussions
For case III, the clustering results of MFM with = 3 and = 1 are shown in Fig. 13(b), (d  = 1 have achieved high classification accuracies, therefore, the improvement of method based on multifractal DFA is not noticeable. In case II, the classification results of 4 methods with high polynomial order are more accurate than those with polynomial order = 1.While in case III and IV, the DFA-based methods [12][13][14] cannot identify the samples accurately since the ( ) and segmented slopes based on mono-fractal logarithmic curves are similar, as shown in Fig. 14 and Fig. 19; due to the nonlinear characteristics of high-dimensional features, even the multifractal theory is introduced into DFA-PC method the recognition performance of MF-PC method is not improved obviously compared to DFA-PC method [12,13].
In 4 cases, the proposed MFM method with higher polynomial orders ( ≥ 2) has a superior performance in signal pattern recognition, which is benefit for rotating machinery fault diagnosis.

Conclusions
In this study, a multifractal manifold method based on detrended fluctuation analysis is proposed to reveal the nonlinear characteristics of vibration signals.The fault diagnosis performance of the MFM method is compared with other DFA-based methods and verified by fault experiments of gearbox and rolling element bearing.The experiment and verification results indicate that a higher polynomial order can improve the recognition results of signal patterns, which is ignored in DFA-based methods [12][13][14].The low-dimensional MFM features can accurately represent the patterns of different signal samples, which are benefit for automatic pattern recognition and multiple fault diagnosis.The comparison results indicate that the MFM method can overcome the disadvantages of DFA-based methods and has a superior performance in signal pattern recognition.YI FENG, BAOCHUN LU, DENGFENG ZHANG

17 .
(a).The multifractal ( )logarithmic curves and multifractal features ℎ( , ) are shown in Fig.17(b) ( = 4).a) Signals b) Multifractal logarithmic curves Fig. Signals of bearings The clustering results of MFM and MF-PC are shown in Fig. 18.The first 3 low-dimensional features of MFM ( = 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 = 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).

1 Fig. 18 . 1 Fig. 19 . 1 Fig. 20 .
features with = 4 b) MFM with = 4 c) MFM features with = 1 d) MFM with = 1 e) MF-PC with = 4 f) MF-PC with = Process results a) = 4 b) = Mono-fractal logarithmic curves of bearing signals characteristics with = 3 b) Crossover characteristics with = 1 c) DFA-PC with = 3 d) DFA-PC with = Process results ), and the clustering result with = 2 is shown in Fig. 21.For case IV, the clustering results of MFM with = 4 and = 1 are shown in Fig. 18(b), (d), and the clustering results with = 3 and = 2 are shown in Fig.22(a), (b).The comparisons show that the mixed samples are gradually separated from each other with the increase of polynomial order.In 4 cases, the recognition performances of 4 methods are improved in varying degrees with higher polynomial orders.

Fig. 21 .
Fig. 21.MFM of case III with = 2As the mono-fractal logarithmic curves of case I shown in Fig.5, the logarithmic curves of different patterns are easily identified and the DFA-based methods with linear polynomial order

2 Fig. 22 .
a) MFM with = 3 b) MFM with = Process results of case IV MANIFOLD FOR ROTATING MACHINERY FAULT DIAGNOSIS BASED ON DETRENDED FLUCTUATION ANALYSIS.

Table 3 .
Classification results

Table 4 .
Classification results

Table 6 .
Classification results

Table 7 .
Parameters of bearings

Table 8 .
Classification results