A novel method for self-adaptive feature extraction using scaling crossover characteristics of signals and combining with LS-SVM for multi-fault diagnosis of gearbox

2.1. Describe detrended fluctuation analysis (DFA)

Since the seminal study on DFA was proposed by Peng et al. [29], this technique has been widely used in time series analysis. The DFA aims to improve the evaluation of correlations in a time series by eliminating trends in the data. The method is described as follows.

1) Consider a signal $x\left(i\right)$, where $i=\text{1,}\text{…}\text{,}N$, and $N$ is the length of the signal. An integrated time series $\tilde{x}\left(i\right)$ of $x\left(i\right)$ is obtained as:

where$⟨x⟩$ is the mean of $x\left(i\right)$ given as:

2) Split series $\tilde{x}\left(i\right)$ into $N_s=\left[N/s\right]$ non-overlapping windows with$s$ data points. $\left[\cdot \right]$ denotes the round numbers. If $N$ is not divisible by $s$, there are some remaining $r=N-N_{s}s$ values at the end of the series $\tilde{x}\left(i\right)$. To solve this problem, the other $N_s$ segments are acquired but starting from $\tilde{x}\left(r+1\right)$. In this way, $\text{2}{N}_s$ windows of $s$ data points are obtained.

3) Let $v$ be the index of the $\text{2}{N}_s$ window. For each one of these windows it takes the polynomial $y_v\left(i\right)$ of degree $m$ to fit the data in the window where $i$ is the data index. Then the $\text{2}{N}_s$ local variance is obtained as Eq. (3). If $N$ is divided by$s$, it simply repeats the $N_s$ values of $F^2\left(v,s\right)$:

where $v=\text{1,}\text{…}\text{,}\text{2}{N}_s$, must be the same for every step of this technique, determining the detrend polynomial order of analysis.

The averaged root mean square value, i.e. fluctuation function $F\left(s\right)$ for the entire time series is given as:

4) Repeated all the above steps 1-3 for several values of $s$ in order to construct the relation between fluctuation function $F\left(s\right)$ and time scale $s$. The relation is presented in the form as:

By a double logarithmic operation, the parameter $\mathrm{\alpha }$ which is called the Hurst exponent or scale exponent can be evaluated. Simultaneously the relation between the fluctuation function and the time scale can be illustrated graphically by the scaling-law curve in a log-log plot. The value of $\alpha =\text{0.5}$ indicates that there are no correlations or only short-term correlations [17]. If $\alpha >\text{0.5}$, the data are long-term correlated. With the increasing of $\mathrm{\alpha }$ the persistent long-range correlation of data is enhanced. The case of $\alpha <\text{0.5}$ corresponds to long-term anti-correlations, meaning that large fluctuation is more possible to be followed by small values.

In fact, the complex fractal time series generally do not exhibit mono-scaling behavior characterized by a single scaling exponent. Multi-scaling behaviors are actually very common in natural phenomena. An example for vibration signal of multi-fault gearbox and its scaling-law curve analyzed by DFA is shown in Fig. 1. There are two obvious crossover points which are shown as the square points described in Fig. 1(b). As a result, the two crossover points can be employed to divide the whole scaling-law curve into three different scale regions, in each of which a single scale exponent can be estimated. However, it hardly determines the scale exponent without a theoretical model for multi-scaling behaviors. If the determinations of these regions are subjective, the results derived from different analysis of the same series can not be consistent. Therefore, an independent criterion to estimate the optimal scale region of scaling-law curve is developed in the Subsection 2.2.

Fig. 1Example for the vibration data of muti-fault gearbox with multi-scaling behaviors

a) Time domain signal

b) Scaling-law curve obtained by DFA

2.2. The method of uncovering optimal scale intervals

The optimal scale regions of time series with crossover characteristics can be objectively and efficiently estimated by using the method of Quasi-Monte Carlo. The Quasi-Monte Carlo method is described as follows.

1) Suppose that $x\left(i\right)$ is a time series which has $M$ given scales $s\in \left\{s_1,s_2,\dots ,s_M\right\}$ where the fluctuation function $F\left(s\right)\in \left\{F\left(s_k\right)\right\}$, $\left(\text{1}\le k\le M\right)$ is calculated.

2) Define a number of $\delta$ data points. In practice, the parameter $\delta$ must be less than the length of minimum scale region. A detailed discussion with the selection of parameter $\delta$ is given at the end of this subsection.

3) Calculate the logarithmic series $L_s\in \left\{\mathrm{l}\mathrm{o}\mathrm{g}\left(s_k\right)\right\}$ and $L_F\in \left\{\mathrm{l}\mathrm{o}\mathrm{g}\left(F\left(s_k\right)\right)\right\}$, $\left(\text{1}\le k\le M\right)$.

4) Define a matrix $r_{M\times M}$ of default values equal to 0.

5) Compute all non-zero elements of $r_{M\times M}$ according as:

where $R^2\left(s_i,s_j\right)$ is the coefficient of determination with the linear fit of $L_F$ versus $L_s$ between the logarithmic time scale $\mathrm{l}\mathrm{o}\mathrm{g}\left(s_i\right)$ and $\mathrm{l}\mathrm{o}\mathrm{g}\left(s_j\right)$, $\text{1}\le i\le (M-\delta +\text{1})$ and $i+\delta -\text{1}\le j\le M$.

6) Sort all non-zero values $r_{ij}$ in decreasing order while keeping a record of their original subindices. The first element of this list will then provide the first optimal fitting scale interval. If there are repeated values in $r_{M\times M}$, sorting requires an extra criterion that the longest interval length is selected. The slope of linear fit is regard as scale exponent$\mathrm{\alpha }$ of the scale interval.

7) Eliminate searched interval in step 6. Then, repeat above steps 2 to 6. The next optimal region can be searched in remaining interval until all linear regions are uncovered.

This algorithm is called as Quasi-Monte Carlo since it actually performs all possible fits at least $\delta$ data points between scales $s_1$ and $s_M$. In order to more clearly clarify how to operate the method for multi-fault vibration signals, there are two projects used to illustrate the operational process subsequently.

A simple artificial series is performed by this method. The 126 equally spaced points $x_i$ are generated in interval of [1, 3.5]. Then, the Eq. (7) is calculated.

where ${\sigma }_{1i}$, ${\sigma }_{2i}$, ${\sigma }_{3i}$ are generated as a set of uniform random numbers of variances ${\sigma }_1$, ${\sigma }_2$, ${\sigma }_3$ used as the additive noise. Then the artificial fluctuation functions $F\left(s_i\right)$ and scales $s_i$ are generated with series $x_i$ and $y_i$ according to:

Consider two kinds of situations. One is that the relationship of variances for random numbers ${\sigma }_{1i}$, ${\sigma }_{2i}$, ${\sigma }_{3i}$, is ${\sigma }_1\le {\sigma }_2\le {\sigma }_3$ and the other one is ${\sigma }_2\le {\sigma }_1\le {\sigma }_3$. The $\delta$ data points are set as 18. The analysis results of artificial series by the method of uncovering optimal fitting scale intervals are shown in Fig. 2 and Fig. 3 respectively. It can be seen that that the fitted curves are perfectly coincident with original data depicted in Fig. 2(b) and Fig. 3(b). The slopes of two examples are (0.5464, 0.2021, 0.9377) and (0.2002, 0.5505, 0.9228) which approximate to the genuine value (0.55, 0.20, 0.90). Moreover, the order of the obtained linear region is in accordance with the sizes of ${\sigma }_1$, ${\sigma }_2$, ${\sigma }_3$. As seen in Fig. 2(c) and Fig. 3(c), the coefficients of determination approach to one. It demonstrates that the method of uncovering optimal scale intervals is reliable.

It should be remarked that the key of the method for uncovering optimal scale intervals is to set up the parameter $\delta$. As mentioned in paper [26], the computing time $\mathrm{t}$ for this implementation might be approximated by:

where $k$ is constant. $\#C$ is the number of active cores. This inverse proportionality is due to the fact that elements in the matrix $r_{M\times M}$ are independent of each other, therefore enabling the possibility to split the processing task.

Fig. 2Analysis results of series with noise

a) Series with noise ${\sigma }_1\le {\sigma }_2\le {\sigma }_3$

b) Fitting result by the method of uncovering optimal fitting scale intervals

c) Coefficient of determination of the linear fit

Fig. 3Analysis result of series with noise

a) Series with noise ${\sigma }_2\le {\sigma }_1\le {\sigma }_3$

b) Fitting result by the method of uncovering optimal fitting scale intervals

c) Coefficient of determination of the linear fit

According to Eq. (9), the large parameter $\delta$ is conducive to reducing computing time. Moreover, the reasonable large parameter $\delta$ makes more obvious statistical sense for the linear fit. However, the parameter $\delta$ should be less than the data points of the shortest scale regions. Because the parameter $\delta$ is too large it could cause a poor statistical sense. As delineated in Eq. (7), the data points of the scale region are 42. Therefore, the $\delta$ data points are set as 18 which is a tradeoff between computing time and statistical sense.

2.3. The method of LS-SVM

As an improvement of Vapnik’s standard SVM [30], LS-SVM classifier [31] leads to solve a linear system instead of a quadratic optimization method. Its objective function is formulated as:

where $z_i\in \{-\text{1,}\text{1}\}$ stands for the $\mathrm{i}$th desired output. $\phi \left(y_i\right)$ denotes a nonlinear mapping of the $i$th input sample of the training data set $\{y_i,z_i{\}}_{i=1}^n$ from the primal space to the feature space. $C$ is viewed as a form of a “regularization” parameter which controls the tradeoff between the complexity of the machine and the number of non-separable points. $\mathbf{\omega }$ is a vector which determines the orientation of a separating hyperplane. ${\xi }_i$ is the slack variables.

In general, $\mathbf{\omega }$ may potentially become infinite dimensional so that a separating hyper-plane does not exist. Hence the constrain is relaxed by introducing slack variables ${\xi }_i\ge 0\text{.}$ The corresponding Lagrangian equation is built as:

where ${\beta }_i$ ($i=\text{1,}\text{…}\text{,}n$) is the Lagrange multiplier.

Then support vectors are got by Karush-Kuhn-Tucker (KKT) condition. The classifier in the dual space takes the form as Eq. (12). The detailed procedure can consult in paper [32]:

where $K\left(y,y_i\right)$ is called the kernel function that must satisfy Mercer condition.

2.4. Frame of the novel method for feature extraction using crossover characteristics

As a result, the method of uncovering optimal fitting scale intervals can be employed to divide the whole scaling-law curve into several different scaling regions, in each of which a single scale exponent can be estimated. Taking advantage of the natural crossover characteristics of the scaling-law curve, a novel method for self-adaptive feature extraction using crossover characteristics is proposed in this study. The whole procedure of the proposed novel method can be decomposed as the following four steps:

1) Analyze the vibration signal of multi-fault gearbox by DFA and obtain the scaling-law curve described in Subsection 2.1.

2) Segment the entire scaling-law curve into several different scaling regions by the method of uncovering optimal scale intervals demonstrated in Subsection 2.2.

3) Extract the scale exponents in each scale region as the feature parameters and utilize them to characterize the fault types.

4) Use the method of LS-SVM to classify the feature parameters obtained in the third step.

4.1. Extract the feature parameters from optimal fitting scale intervals

The vital step of the proposed method is the extraction of feature parameters. In this subsection, the method of uncovering optimal scale intervals in Subsection 2.2 is used to analyze the vibration data of multi-fault gearbox. As suggested in paper [16], the degree$\mathrm{m}$ of DFA introduced in Subsection 2.1 which can be obtained by attempting different orders is finally set as one. Additionally, the range for data points$\mathrm{s}$ of non-overlapping windows is from 10 to 3000, i.e. the interval of logarithmic $s$ ($\mathrm{l}\mathrm{o}\mathrm{g}\left(s\right)$) is set to range from 1 to 3, described in paper [17]. In our study, under the same assumptions, the degree$\mathrm{m}$ of DFA is 1 and the range for the data points of non-overlapping windows is from 10 to 3162, i.e. the interval of $\mathrm{l}\mathrm{o}\mathrm{g}\left(s\right)$ is uniformly set to range from 1 to 3.5 with the spacing of 0.02 to observe the analysis results in the larger windows. As depicted in Fig. 1(b) and Fig. 6(a)-(f), the scaling-law curves obtained by DFA for all fault types of gearbox clearly exhibit two crossover points. The scaling-law curves can be roughly divided into three scaling regions by using these crossover points where the $x$-coordinates of two crossover points are close to $\mathrm{l}\mathrm{o}\mathrm{g}\left(s_{25}\right)$ and $\mathrm{l}\mathrm{o}\mathrm{g}\left(s_{100}\right)=\text{3}$ respectively. Therefore, the ranges of three scale regions approximate to $\left[s_1,s_{25}\right]$, $\left[s_{25},s_{100}\right]$, and $\left[s_{101},s_{126}\right]$ respectively. As displayed in Fig. 6(a)-(f), the first and third linear regions have a shorter length than the second region. The shortest length of these linear intervals approximates to 25. To uncovering the optimal scale intervals, the $\delta$ data points are set as 18 which is a tradeoff between computing time and statistical sense as explained in Subsection 2.2.

Fig. 7Coefficients of determination (R2) of the linear fitting under load-1 with rotational speed 880 rpm

a) Type 1

b) Type 2

c) Type 3

d) Type 4

e) Type 5

f) Type 6

As shown in Fig. 6, the black heavy lines represent the optimal scale intervals fitted by the criterion based on a Quasi-Monte Carlo algorithm. Fig. 6(a)-(f) manifest that the optimal fitted lines can completely depict the outline of scaling-law curves. It can be seen from Fig. 7(a)-(f) that the coefficients of determination before the non-overlapping window of $\mathrm{l}\mathrm{o}\mathrm{g}\left(s_{100}\right)$ are more apparent than the remaining parts. That is, the linearity of first and second linear regions $\left[s_1,s_{25}\right]$, $\left[s_{25},s_{100}\right]$ has a better performance than the one of the third region $\left[s_{101},s_{126}\right]$. As a consequence, the order of the obtained linear regions may be $\left[s_1,s_{25}\right]$, $\left[s_{25},s_{100}\right]$ and $\left[s_{101},s_{126}\right]$ or $\left[s_{25},s_{100}\right]$, $\left[s_1,s_{25}\right]$ and $\left[s_{101},s_{126}\right]$ which are consistent with analysis in Subsection 2.2. Lastly, three scale exponents, in turn, are estimated from these linear regions by least square fitting. Afterwards three scale exponents are determined, the sequence of scale exponents should be rearranged to conform the original order $\left[s_1,s_{25}\right]$, $\left[s_{25},s_{100}\right]$ and $\left[s_{101},s_{126}\right]$ for uniformly forming a standard vector employed in Subsection 4.2 as the input samples. Table 2 reports the three scale exponents for six fault types of gearbox under load-1 with rotational speed of 880 rpm. As reported in Table 2, the scale exponents in the first two scaling regions are less than 0.5 and hence the corresponding data exhibit extremely strong anti-persistent long-range correlations while the opposite situation is shown in third linear region. In next subsectioon, we will present the discriminability of the features vector $V_1$ which is expressed as:

where $SE1$, $SE2$ and $SE3$ mean the three scale exponents in sequence.

4.2. Compare the classified peformance of proposed method with statistical parameters

In general, it is difficult to identify fault types by directly observing features. In this section, LS-SVM are used to realize the recognition for multi-fault gearbox.

As described in Subsection 2.3, the primary principle of LS-SVM which is an intelligent learning algorithm and derived from the work of Vapnik et al. [30, 35] is to map the input from the primal space to the feature space by kernel function shown in Eq. (12). The selection of kernel function $K\left(y,y_i\right)$ will affect the learning ability and generalization of LS-SVM. Since the Gaussian Radial Basis Function (RBF) kernel always have a superior performance than the other kinds of kernel functions in many practical applications [36], the RBF kernel which is shown in Eq. (13) is adopted in our research. As well, the RBF kernel is a better choice than other kernels like polynomial kernel because it has lesser hyper-parameters and so the problem becomes less computationally intensive. In addition, the LS-SVM is designed intrinsically for binary classification which is not suitable for multi-class of multi-fault diagnosis. However, in our paper, a six-class classification problem should be implemented by extending the binary classification. Many approaches have been considered for multi-class classification issue, such as “one-against-one”, “one-against-all”, and directed acyclic graph SVM [37]. “one-against-one” is to break down the multi-class problem into a number of smaller binary problems. “one-against-all” considers all classes at once and solves the multi-class problem in one step. In general, “one-against-all” is computationally more expensive to solve a multi-class problem than “one-against-one”. Hsu and Lin [38] discussed the characteristics of these methods and pointed out that the “one-against-one” method is more suitable for practical application than other methods. In this study, we select the “one-against-one” to identify the different fault types:

where $\sigma$ is the bandwidth of the Gaussian RBF kernel function which controls the degree of non-linearity of the model.

In order to make fair and objective evaluation for each type of multi-fault gearbox, two key parameters in LS-SVM, $C$ and $\sigma$, which play an important role in the classification performance for LS-SVM are optimized by grid searching [39]. The search scopes of parameters $C$ and $\sigma$ are between 0 and 2 with spacing of 0.1. As mentioned in Section 3, 44 data samples for each operating condition of one fault type are collected from the surface of gearbox. 44 samples are divided into training samples and testing samples, respectively. 15 samples are selected as training samples and the remaining 29 samples are regarded as testing samples. Therefore, a total number of training samples are 90 and testing samples are 174 for six fault types of gearbox under one working condition. The objective function of grid searching method is the average percentage of correct classification, i.e. the ratio of the total number of training samples correctly classified to the total number of training samples.

As displayed in Fig. 8, five of the training samples of $V_1$ are misclassified, yielding the maximum average percentage of classification accuracy (CA) which is calculated as Eq. (15) is 94.4 % under load-1 with rotational speed of 880 rpm when $C$ and $\sigma$ are equal to 1.1 and 1.6 respectively:

Fig. 8Optimal training results are obtained by grid searching in LS-SVM

Fig. 9Classification results in testing samples by LS-SVM

To illustrate the process concretely, the classification results of $V1$ in testing phases are drawn in Fig. 9. It can be seen from Fig. 9 that 16 of the testing samples are misclassified, yielding the success rates 90.8 %. Among them, 2 samples of type-1, 3 samples of type-3, 3 samples of type-5 and 8 samples of type-6 are misclassified respectively. For brevity, the CA of testing samples obtained by the proposed method for all types of multi-fault data are listed in Tables 3-4.

Additionally, compared with the proposed method, the vector $V_2$ of statistical parameters expressed as Eq. (16) is employed to examine the same multi-fault data which have been investigated by proposed method:

where $P$, $K$ and $F$ denote pulse factor, kurtosis and form factor, respectively.

The main parameters in LS-SVM are set in the same way as ones in our proposed method. The analysis results of testing samples about vector $V_2$ are demonstrated in Table 5 and Table 6.

As illustrated in Table 3 and Table 5, the classified performances of our proposed method are obviously over the performances of statistical parameters in all fault types of gearbox except the slightly less CA for the load-1 of type-5 and load-2 of type-6. Similarly, the comparisons of Table 4 and Table 6 clearly indicate that the performances of proposed method are also over the performances of statistical parameters at the absence of external load. Especially, zero correct classification crops up at type-4 and type-6 of 880 rpm shown in Table 6. As observed in Table 3, the classified performance of heavy load (load-1) is better than light load (load-2). Table 4 presents that the working condition with highly rotational speed has a slightly helpful for improving CA. The contrasts of Table 3 and Table 4 imply that the performances of the absence of external load are slightly poor. The classified results of statistical parameters are unstable and lowly accurate depicted in Table 6.

Successively, to further improve the classified performance of proposed method, one idea that vector $V_3$ which is form as Eq. (17) (the combination of three scale exponents $V_1$ and three statistical parameters $V_2$) is attempted to used to characterize the types of multi-fault gearbox:

where $V_3$ represent the combination of vector $V_1$ and vector $V_2$.

The main parameters in LS-SVM are set in the same way as ones in the proposed method. The classified performances of vector $V_3$ are presented in Table 7 and Table 8. From the comparison of Table 3 and Table 7, the combination of these parameters has a better performance than the single vector $V_1$. Identically, Table 4 and Table 8 state that the performances of novel idea are superior to the single vector $V_1$ under working condition of the absence external load. As reported in Table 7, it must be point out that the CA is 100 % for all fault types under load-1 with rotational speed of 880 rpm.