Fault diagnosis of mechanical drives under non-stationary conditions based on manifold learning of kernel mapping

2.1. The definition of manifold learning

Manifold learning is to explore the intrinsic geometry information of dataset, i.e. to discover the inherent low-dimensional manifold embedded in the high-dimensional observation space, and the formal description is given as follows:

In a high-dimensional observation space, let $\mathbf{Y}$ be a $d$-dimensional domain contained in the Euclidean space and let $f:\mathbf{Y}\to R^d$be a smooth embedding. The objective of manifold learning is to recover $f$ and $\mathbf{Y}$ based on a given set of observed dataset $x_i$. The observed data arise as follows. Hidden data are generated randomly in $\mathbf{Y}$, and then mapped by $f$ to become the observed data, so $x_i=f\left(y_i\right)$.

In many cases it seems to be a reasonable assumption that the data lie on or close to a manifold, and the samples in different states also have different embedded manifolds. The embedded manifold is considered as a sub-manifold, which has a manifold structure itself. Let $\mathbf{N}$ be an $n$-dimensional manifold, and let $k$ be an integer such that $0\le k\le n$. A $k$-dimensional embedded sub-manifold of $\mathbf{N}$ is a subset $\mathbf{M}\subset \mathbf{N}$. Then, for different embedded manifolds, if there are no crosses each other, thus they can be unfolded with mapping transformation.

2.2. The difference degree of embedded manifolds

It was clear that if the manifold coordinate corresponding to the different fault classes has a large difference, then the fault pattern can be classified easily. So the evaluation of difference degree between embedded manifolds is proposed with Euclidean distance. Let original data $\mathbf{X}$, and manifold coordinate $\mathbf{Y}$, then the difference degree $J_D$ is defined as follows:

where ${\mathbf{H}}_{ij}=\left\{\begin{array}{cc}0& {\mathbf{Y}}_i,{\mathbf{Y}}_j\in p\\ 1& \text{others}\end{array}\right.$, and $\mathbf{H}$ represent the class lable, which is 0 for same elements and 1 for different elements. In order to improve the separability of embedded manifolds, it is necessary to increase the number of nonzero in function $J_D$.

2.3. The manifold learning for varying operation condition

At the viewpoint of geometry, the vibration data in the same state has the same geometric property in space distribution or topological structure. Its mapping points in the low-dimensional space can be distributed in an embedded manifold or in its neighborhoods. So the samples in different faults have different embedded manifolds. However, under the varying loads and speeds, the vibration patterns become ambiguous. Then for the dimension reduction, the difference degree of embedded mainfolds corresponding to vibration patterns maybe decrease, or overlap. The direct result is to cause erroneous judgement for new data. Therefore, in order to overcome the effect of varying operating conditions, the learning strategy of embedded manifold that can improve class separability become an important problem.

3.1. Background of kernel mapping

For machine learning, the kernel trick is a way of mapping observations from a general set into an inner product space, without ever having to compute the mapping explicitly. The trick to avoid the explicit mapping is to use learning algorithms that only require dot products between vectors in product space, and choose the mapping so that these high-dimensional dot products can be computed within the original space, by means of a kernel function. With respect to keeping the nonlinear topology structure, the feature space is projected to the Hilbert space by kernel mapping, and the class separability is also improved. Therefore, the kernel mapping strategy is transplanted to a manifold learning for the purpose of improving class separability.

3.2. Influence of the difference degree in kernel space

With the kernel mapping, the difference degrees of the embedded manifolds projected in Hilbert space also change. Let $\mathbf{X}$ denotes the analyzed data matrix, and the embedded manifold data matrix is $\mathbf{Y}$. The inner product matrix is defined by ${\mathbf{K}}_{ij}={\mathbf{Y}}_i\cdot {\mathbf{Y}}_j$, so the difference degree matrix is defined as follows:

where ${\mathbf{K}}_{ij}$ denotes a Gram matrix, and then the inner products between the vectors in feature space are alternatived with calculation of kernel function. The input and the output of manifold learning is related by an unknown matrix $\mathbf{V}$, and the relationship is defined by:

So the Eqs. (2) can be rewritten as:

where ${\mathbf{H}}_{ij}=\left\{\begin{array}{cc}0& \phi \left({\mathbf{X}}_i\right),\phi \left({\mathbf{X}}_i\right)\in p\\ 1& \text{others}\end{array}\right.$. Although the mapping matrix $\mathbf{V}$ is unexplicit, the class similarity in mapping $\phi \text{(}\mathbf{X}\text{)}$ can be improved. Hence the difference degree $J_D$ of the embedded manifold obtained by kernel mapping is larger than the traditional manifold learning.

3.3. Kernel mapping manifold learning for fault diagnosis

Considering the advantage of kernel mapping, the kernel mapping manifold learning (KMML) algorithm is also proposed to reveal the intrinsic structure which has better class separability for varying operation conditions. Therefore, the strategy of fault diagnosis based on KMML can be described as follows:

1) For the faults of mechanical drives, except for time features, the frequency spectrum of the signal may change, which signifies that new frequency components may appear. Thus, diagnosis feature set $\mathbf{X}=\text{{}{\mathbf{x}}_1,...,{\mathbf{x}}_N\text{}}$ is calculated in time domain or frequency domain.

2) The mapping kernel function $\phi$ is used to project the diagnostic feature space to Hilbert space $\mathbf{X}\subseteq R^D\to \mathbf{H}\subseteq R^N$.

3) With the k-nearest neighbor algorithm (KNN), determine $k$ nearest neighbors $\phi \left({\mathbf{x}}_i^j\right)$ of $\phi \left({\mathbf{x}}_i\right)$, $j=$ 1, 2, …, $k$, and set:

4) Let the minimal weight:

5) Let $\mathbf{Q}{\left(\phi \right)}_i$ is the $d$ left singular vectors of matrix $\phi {\left(\mathbf{x}\right)}_i-\phi \left({\stackrel{-}{\mathbf{x}}}_i\right){\mathbf{e}}^{\text{T}}$ corresponding to its $d$ largest singular values, and the local coordinates $\mathbf{\Theta }$ is given by $\mathbf{\Theta }{\left(\phi \right)}_i$ defined as:

6) Calculate the alignment matrix $\mathbf{B}\left(\phi \right)=\mathbf{S}\mathbf{W}{\mathbf{W}}^{\text{T}}{\mathbf{S}}^{\text{T}}$, where $\mathbf{S}=[{\mathbf{S}}_1,{\mathbf{S}}_2,\dots {\mathbf{S}}_N]$ and ${\mathbf{S}}_i$ is the selection matrix that $\text{[}\phi \text{(}{\mathbf{x}}_1\text{)},...,\phi \text{(}{\mathbf{x}}_N\text{)}]{\mathbf{S}}_i=\text{[}\phi \text{(}{\mathbf{x}}^1\text{)},...,\phi \text{(}{\mathbf{x}}^k\text{)]}$, and $\mathbf{W}=\text{diag(}{\mathbf{W}}_{\text{1}},\dots ,{\mathbf{W}}_N\text{)}$ with:

7) The low-dimensional embedding $\mathbf{T}$ is given by the $d$ eigenvectors of the alignment matrix $\mathbf{B}$, corresponding to the 2nd to $d$+1st smallest eigenvalues of $\mathbf{B}$.

Finally, for the new testing samples ${\mathbf{X}}_{new}=\text{{}{\mathbf{x}}_{N+1},{\mathbf{x}}_{N+2},\dots ,{\mathbf{x}}_{N+M}\text{}}$, the new embedding can be also obtained as follows:

(1) Calculate the mapping $\phi \left({\mathbf{x}}_{\mathbf{N}+\mathbf{i}}\right)$ and its neighborhood matrix.

(2) Reconstruct weight ${\mathbf{W}}_{ij}$ and search the $k$ embedding ${\mathbf{t}}_j$ in matrix $\mathbf{T}$ according to the neighborhood of $\phi \text{(}{\mathbf{x}}_{N+i}\text{)}$.

(3) Combine ${\mathbf{W}}_{ij}$ with embedding ${\mathbf{t}}_j$ to calculate the new embedding which is defined as ${\mathbf{t}}_{new}^i=\sum _j{\mathbf{W}}_{ij}{\mathbf{t}}_j$, and update the optimal embedding $\mathbf{T}$ with $\mathbf{T}=\text{(}\mathbf{T}\text{,}{\mathbf{t}}_{new}^i\text{)}$.

(4) The KNN algorithm is adopted to classify in embedding $\mathbf{T}$.

5.1. Description of the simulator

The experiments were conducted on Spectra Quest’s Gearbox Dynamics Simulator, as shown in Fig. 5(a). The objective is to relate the development of the feature component changes in varying operation conditions. The gearbox is driven by a 3 HP motor with a speed range of 0-4000 rpm. The speed inverter can be programmed, and the load provided by the brake was controlled by a current source. The modular design of the gearbox allows one to simulate various faults, either individually or jointly in a totally controlled environment. Fig. 5(b) illustrates a picture taken for the inside of the gearbox and a schematic of the two-stage parallel gearbox layout.

Two gears were tested in this study: (1) a health gear, (2) a large-chipped gear. The defected gears shown in Fig. 6 are the input gear (the 32-tooth one). For gear configuration, two operation conditions (varying load, varying load and speed simultaneously) were tested. For each test, accelerometer was mounted on the bearing seat. The vibration data was collected by a Noise & Vibration Analyzer – OR35, where the sampling frequency is 10 kHz.

Fig. 5Gearbox dynamics simulator

(a) Experiment table structure

(b) Two-stage parallel gear transmission

Fig. 6Faulted gears for the chipped tooth

5.2. The experiment of varying rotation speed

For the rotation speed of simulator, the cosine function is adopted as the model of varying speed which changes from 10 Hz to 20 Hz in 65 seconds. Therefore, the speed value per second is calculated according to the speed-time curve shown in Fig. 7, for the time resolution is 1 second and speed resolution is 0.1 Hz in the software of motor PC-control program. Finally, the rotation speed is programmed by the 65 discrete values of speed. Fig. 8 shows the waveform of the signal obtained from the gearbox in two states with varying rotation speed. It can obviously be seen that the amplitude of the waveform increases with the increase of rotation speed. However, the vibration level in Fig. 8(b) is the same as that in Fig. 8(a) with the effect of the varying speed.

Fig. 7Relationship of the speed and the time

Fig. 8Waveform obtained from two operation states in varying speed

(a) Normal state

(b) Chipped tooth state

For such non-stationary signal, the time-frequency representation is shown in Fig. 9. However, it is difficult to visually identify the fault state from normal operating data. Thus in order to distinguish gearbox operation state, an appropriate dimension reduction tool is needed to extract the low-dimensional embedding.

Fig. 9Time-frequency representation of the signal in two operation states

(a) Normal state

(b) Chipped tooth state

The vibration data is divided by slide Hanning window whose window width is 10000, and the degree of overlapping is 50 % (to ensure the continuous samples). Thus the set of training data consisting of 119 samples is constructed for each state (normal and chipped). The frequency spectra are divided into 20 equal bands, and then the 20 spectral energy features are taken as the statistical features to characterize the running state pattern.

By the influence of the rotation speed, the Pseudo degree of freedom also increases. Therefore, in order to visually identify the fault, the intrinsic dimension is set to 2. The low-dimensional embedding results obtained by LTSA and KMML are compared, where the kernel function is gauss function with variance 0.01. Fig. 10(a) shows the two-dimensional representation of normal and chipped states. It embraces an overlap around the higher speeds. This may results in false diagnosis in the visual coordinate space. The main reason is that, in the higher speeds, the features of frequency band energy of chipped tooth are similar to that in the normal state, thus due to the influence of space reduction, the feature set of the two operation states can emerge overlap partly.

Fig. 10Two-dimensional representation for different learning algorithm

(a) LTSA

(b) KMML

Fig. 10(b) shows the two-dimensional representation of the feature set obtained by KMML method. It can obviously be seen that the embedded manifold coordinate has a significant difference between normal and chipped-tooth states. Due to kernel mapping, the class separability can be enhanced. Even though the embedding value increase with the rotation speed increases in embedding element ${\mathbf{T}}_1$coordinate direction, the embedded manifold of the two running conditions also has no overlaps. It is very convenient to identify the fault pattern in the varying speed condition. Finally, the 33 test samples are classified correctly.

5.3. The experiment of variable rotating speed and load

As discussed above, the proposed method is an effective classification for varying speed. However, in practice, except for the varying speed, load of gearbox system also changes by the influence of many factors, such as non-uniformity in material which can result in the load fluctuation of spindle. Then the speed increase can result in the amplitudes increasing of features. So it is difficult to identify this content as loading-related or fault-related. Therefore further experiments with varying loads need to be done.

In the process of varying speeds and loads, loading force is controlled by a current source using the magnetic brake. Fig. 11 shows the varying curve of rotation speed and load. It can be seen that the current values of load force change from 0 A to 40 A in 130 seconds, where it is large in slow speeds and vice versa. Fig. 12 shows the waveform of the signal obtained from the two running states, where sampling frequency is 6400 Hz. Each plot depicts the changes of waveform in varying speed and load condition. It can obviously be seen that the amplitude of vibration waveform increases with the increase of rotation speed. The vibration level in Fig. 12(b) is also the same as that in Fig. 12(a).

Fig. 11Relationship between speed-time and current-time

(a) Speed curve

(b) Current curve

Fig. 12Waveform of signals obtained from two operation states in varying speeds and loads

(a) Normal state

(b) Chipped tooth state

It should be noted that classic diagnostic features usually depend on the value of the external load. Fig. 13 shows the diagnostic features [16] (standard deviation, kurtosis, root mean square, absolute mean, peak factor, waveform factor, and margin factor) of 200 samples of varying speeds and loads in two operating states. The vibration level of each variable in normal (in Fig. 13(a)) is compared with that of chipped fault (in Fig. 13(b)). It was clear that considering the vibration level of the variable itself, the fault pattern cannot be detected or identified.

Considering the effect of speed and load, the vibration data is divided by sliding window, where the window width is 6400, and degree of overlapping 50 % is used. The analysis sample number for each operation state is set to 240, and then each sample is transformed by the 'db3' wavelet which the scale $a=$ 1, 2, …, 20. Thus in the discrete 20 scales, the features of sample are calculated as the sum of absolute of the wavelet coefficients and mathematically is defined by $E_a=\sum _{i=1}^n\left|\mathbf{W}\right({\tau }_i,a\left)\right|(a=\mathrm{1,2},...,20)$ where $\mathbf{W}$ are wavelet coefficients. After that, a total of 480 samples including two states were considered for the 20 features, thus the high-dimensional feature space $\mathbf{X}\in R^{N\times m}$, where $N=$ 480 and $m=$ 20.

Fig. 13Feature value obtained from two operation states in varying speeds and loads

(a) Normal state

(b) Chipped tooth state

The two-dimensional representation of the feature set obtained by KMML is shown in Fig. 14(a). It can obviously be seen that, due to increasing the degree of freedom of embedding in the condition of varying speed and load, the embedded manifolds captured by KMML form two operation states have a larger overlap. This means that it is difficult to have a good recognition effect. Therefore, the three-dimensional embedded manifold is also extracted and shown in Fig. 14(b). It is very convenient to identify that the overlap of two embedding distribution is much lower. The 20 test samples are adopted to be classified and the classification error is only 5 %.

Obviously, in embedding space, the misjudgment of two states only occurs at overlapping region for the similar the embedding coordinates, such as the samples of normal and chipped in 11.9 Hz speed. The corresponding original 20-dimensional features are shown in Fig. 15. The number 1, 2, …, 20 of the horizontal axis in the graph respectively stands for different scales, and vertical axis is the feature value under scales. It is found that the feature values under normal and chipped states are almost the same representation. Therefore, in the condition of varying speed and load, except for the samples in overlapping region, other samples can be classified correctly for the larger different of embedding coordinates.

Fig. 14Representation for varying operation conditions with different embedding spaces

(a) Two-dimensional embedding

(b) Three-dimensional embedding

Fig. 15Frequency domain feature in the condition of speed is 11.9 Hz