Local coordinate weight reconstruction for rolling bearing fault diagnosis

2.1. Local coordinate computation

Let sample matrix $\mathbf{X}=\left({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\right)\in {\mathbf{R}}^{D\times N}$, $D$ is the number of sample dimensions and $N$ is the number of samples. Sample ${\mathbf{x}}_i$ and its $k$ nearest samples (including ${\mathbf{x}}_i$) constitute local neighborhood matrix ${\mathbf{X}}_i=({\mathbf{x}}_{i1},{\mathbf{x}}_{i2},\dots ,{\mathbf{x}}_{ik})\in {\mathbf{R}}^{D\times k}$, each neighborhood exists a local tangent space ${\mathbf{Q}}_i\in {\mathbf{R}}^{D\times d}$$(d<D)$ consisting of standard orthogonal basis vectors, and ${\mathbf{\Theta }}_i=\left({\mathbf{\theta }}_{i1},{\mathbf{\theta }}_{i2},\dots ,{\mathbf{\theta }}_{ik}\right)\in {\mathbf{R}}^{d\times k}$ is the projection of ${\mathbf{X}}_i$ onto ${\mathbf{Q}}_i$. For the main geometric structure information within neighborhood, minimize the sum of square of distance between ${\mathbf{X}}_i$ and ${\mathbf{\Theta }}_i$:

where $\overline{{\mathbf{X}}_i}=1/k{\mathbf{X}}_i\mathbf{e}{\mathbf{e}}^T$ is a mean matrix of ${\mathbf{X}}_i$, $\mathbf{e}$ is a unit column matrix of length $k$, and ${\mathbf{I}}_d$ is a $d$ by $d$ unit matrix.

Apply singular value decomposition to ${\mathbf{X}}_i-\overline{{\mathbf{X}}_i}$:

${\mathbf{Q}}_i{\mathbf{\Theta }}_i$, ${\mathbf{Q}}_i$, ${\mathbf{\Theta }}_i$ are computed as below:

where $\mathbf{U}\in {\mathbf{R}}^{D\times D}$ is a left singular vector set, $\mathbf{V}\in {\mathbf{R}}^{k\times k}$ is a right singular vector set, $\sum \in {\mathbf{R}}^{D\times k}$ is a singular value diagonal matrix, ${\sum }_d\in {\mathbf{R}}^{d\times d}$ is a diagonal matrix with $d$ maximum singular values in descending order, ${\mathbf{U}}_d\in {\mathbf{R}}^{D\times d}$ is the corresponding left singular vector set and ${\mathbf{V}}_d\in {\mathbf{R}}^{k\times d}$ is right singular vector set.

Thus, ${\mathbf{\Theta }}_i$ holds the most important geometric structure information in ${\mathbf{X}}_i$.

2.2. Global coordinate construction

Suppose ${\mathbf{T}}_i=\left({\tau }_{i1},{\tau }_{i2},\dots ,{\tau }_{ik}\right)\in {\mathbf{R}}^{d\times k}$ is the local low-dimensional coordinate of ${\mathbf{X}}_i$. Build the affine transformation between ${\mathbf{T}}_i$ and ${\mathbf{\Theta }}_i$:

The local reconstruction error ${\mathbf{E}}_i$ is written as:

where ${\mathbf{L}}_i\in {\mathbf{R}}^{d\times d}$ is a local affine matrix, $\overline{{\mathbf{T}}_i}=1/k{\mathbf{T}}_i\mathbf{e}{\mathbf{e}}^T$ is a mean matrix of ${\mathbf{T}}_i$, ${\mathbf{E}}_i\in {\mathbf{R}}^{d\times k}$ is a local reconstruction error matrix.

However, owing to the different contribution of samples to their manifold structure, a novel LCWR based on weight coefficient is proposed to reduce permutation error and reconstruct local coordinates more accurately. According to LCWR, the closer a sample is to its manifold, the larger its weight coefficient is. Likewise, the farther a sample is from its manifold, the smaller its weight coefficient is. So an exponential function of geodesic distance between a sample and the center point of its neighborhood, reflecting the proximity to its manifold structure, is defined as weight coefficient, namely:

where ${\omega }_{ij}$ denotes the weight coefficients of the $j$th nearest neighbor $x_{ij}$ in ${\mathbf{X}}_i$, $G_{ij}$ and ${\sigma }_i$ denote the geodesic distance from $x_{ij}$ to the center of ${\mathbf{X}}_i$ and the mean square error of ${\mathbf{X}}_i$ respectively, and $\beta$ is adjustment parameter.

Then ${\mathbf{E}}_i$ is rewritten as:

where ${\mathbf{W}}_i=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\sqrt{{\omega }_{i1}},\sqrt{{\omega }_{i2}},\cdots ,\sqrt{{\omega }_{ik}})\in {\mathbf{R}}^{k\times k}$is a weight coefficient matrix, ${\mathbf{H}}_i=(\mathbf{I}-1/k\mathbf{e}{\mathbf{e}}^T){\mathbf{W}}_i\in {\mathbf{R}}^{k\times k}$.

Fix ${\mathbf{T}}_i$ and minimize ${\mathbf{E}}_i$ to preserve as much local information as possible, namely:

where ${\mathbf{L}}_i\in {\mathbf{R}}^{d\times d}$, ${\mathbf{\Phi }}_i={\mathbf{\Theta }}_i{\mathbf{W}}_i\in {\mathbf{R}}^{d\times k}$.

Substitute Eq. (9) into Eq. (7) and obtain ${\mathbf{E}}_i$:

where ${\mathbf{B}}_i={\mathbf{H}}_i-{\mathbf{H}}_i{\mathbf{\Phi }}_i^T\times inv\left({\mathbf{\Phi }}_i{\mathbf{\Phi }}_i^T\right){\mathbf{\Phi }}_i\in {\mathbf{R}}^{k\times k}$.

Minimize the sum of all neighborhoods reconstruction errors to obtain the global low-dimensional coordinate $\mathbf{T}$:

where $\mathbf{T}=\left({\tau }_{1,}{\tau }_{2,}\dots ,{\tau }_N\right)\in R^{d\times N}$ is the global low-dimensional coordinate of $\mathbf{X}$, ${\mathbf{T}}_i=\mathbf{T}{\mathbf{R}}_i$, ${\mathbf{R}}_i=\mathbf{R}\left(:,{\mathbf{I}}_i\right)\in {\mathbf{R}}^{N\times k}$is a selection matrix, $\mathbf{R}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathrm{o}\mathrm{n}\mathrm{e}\left(1,N\right)\right)$, ${\mathbf{I}}_i=\left(i_1,i_1,\dots ,i_k\right)$ is the subscript of $k$ nearest neighbors of $x_i$.

It is equal to solve differential equation:

where $\mathbf{M}={\mathbf{M}}^T={\sum }_{i=1}^N{\mathbf{R}}_i{\mathbf{B}}_i{\mathbf{B}}_i^T{\mathbf{R}}_i^T\in {\mathbf{R}}^{N\times N}$.

Therefore, the optimal solution of ${\mathbf{T}}^T$ is composed of eigenvectors corresponding to the 2nd to the ($d+1$)th smallest eigenvalues of $\mathbf{M}$. ${\mathbf{T}}^T=\left(u_2,u_3,\dots ,u_{d+1}\right)$ or $\mathbf{T}={\left(u_2,u_3,\dots ,u_{d+1}\right)}^T\in {\mathbf{R}}^{d\times N}$ is an orthogonal global low-dimensional coordinate mapping matrix of nonlinear manifold in $\mathbf{X}$, where $u_j\in {\mathbf{R}}^{1\times N}$ is an eigenvector corresponding to the $j$th eigenvalue of $\mathbf{M}$.

LCWR is summarized as follows:

(1) Look for neighborhood. Local neighborhood ${\mathbf{X}}_i$ of each sample ${\mathbf{x}}_i$ ($i=\mathrm{1,2},\dots ,N$) is found by k nearest neighbors.

(2) Extract local coordinates. Projection matrix ${\mathbf{Q}}_i$ and local low-dimensional coordinate matrix ${\mathbf{\Theta }}_i$ of each neighborhood ${\mathbf{X}}_i$ are obtained according to Eq. (2) and Eq. (3), respectively.

(3) Construct global coordinates. Global low-dimensional coordinate matrix $\mathbf{T}$ is conducted from ${\mathbf{\Theta }}_i$ reconstruction by weight coefficient matrix ${\mathbf{W}}_i$ and expressed as Eq. (14).

The flow chart of LCWR is shown in Fig. 1.

3.1. High-dimensional feature construction

As shown in Table 1, 12 time-domain statistical indicators and 8 frequency-domain statistical indicators are selected to constitute a 20-dimensional sample to characterize the bearing state. An original sample signal is decomposed into eight sub-band signals by three-layer db8 wavelet packet decomposition, and the ratio of the energy of each sub-band to the total energy of all sub-bands is taken as the frequency domain indicator. That is, $e{}_j{}^{\mathrm{}}\mathrm{}=E_j/E$, $E={\sum }_{j=1}^8{E}_j$, $E_j$ is the energy of sub-band signals. Thus, a high-dimensional feature matrix $\mathbf{X}\in {\mathbf{R}}^{192\times 20}$ is created.

Fig. 1Procedure of LCWR

Fig. 2LCWR based fault diagnosis

3.2. Low-dimensional feature extraction

According to Fig. 2, some of the low-dimensional feature samples extracted by LCWR are used as training samples to train support vector machine (SVM) while the others as test samples to be recognized by trained SVM. When using LCWR to extract bearing state characteristics, three parameters such as neighbor number $k$, dimension $d$ and adjustment parameter $\beta$ need to be optimized. Because the recognition rate can be regarded as a function of three parameters $k$, $d$ and $\beta$, these parameters interact to determine the recognition rate. By changing these parameters in a certain range and the corresponding recognition rate obtained, it is proved that there exist optimization values of parameters with the peak recognition rate. The trend of recognition rate with respect to a single parameter variable while the other two parameters fixed is shown in Fig. 3-5, respectively. Because of different parameters, the trend of recognition rate is also different from each other.

From the recognition rate about the nearest neighbor number $k$ in Fig. 3, there is the maximum rate 96.43 % at $k=$ 8. The role of $k$ on recognition rate is carried out by influencing the intrinsic geometry structure of high-dimensional samples, the close relationship between similar samples and the nonlinear dimensionality reduction ability of LCWR. If $k$ is too small, LCWR can’t maintain the intrinsic geometry of high-dimensional samples and close association between similar samples. If $k$ is too large, it weakens the nonlinear dimensionality reduction capability of LCWR. Hence, the low-dimensional manifold structure hiding in the high-dimensional samples can be found to the greatest extent and achieve the maximum rate at the optimum value $k=$ 8. However, due to the comprehensive effect of different factors affected by $k$, the recognition rate fluctuates and there are multiple turnover points rather than a monotonous trend.

Dimension $d$ affects the recognition rate by mining the sensitive features of high-dimensional samples in the neighborhood and eliminating redundant and interference components. From Fig. 4, it can be seen that the maximum recognition rate is 96.43 % at $d=$ 3, because an appropriate $d$ makes similar samples have approximate low-dimensional features, leading to better clustering effect and improvement in recognition rate. Otherwise, LCWR can’t fully mine the sensitive features from the high-dimensional samples in neighborhood if $d$ is too small or the low-dimensional features contain redundancy and interference if $d$ is too large. Likewise, due to different factors, the recognition rate has multiple turnover points.

Adjustment parameter $\beta$ affects the recognition rate by changing the degree of clustering and global geometry retention. In relationship between $\beta$ and recognition rate in Fig. 5, if $\beta$ is too small, the proximity of samples is low and the clustering is obvious, but the degree of retention of the global geometric structure is poor. If $\beta$ is too large, the global geometric structure can be improved but the clustering reduced. These factors cause reduction in recognition rate. As a result, there is an optimal $\beta =$ 0.1 where the recognition rate reaches the maximum 96.43 %.

Fig. 3Recognition rate with neighbor k

Fig. 4Recognition rate with dimension d

Fig. 5Recognition rate with parameter β

The three-dimensional sensitive features of bearing high-dimensional characteristic samples reduced by LCWR according to the optimized values of $k=$ 8, $d=$ 3 and $\beta =$ 0.1 are shown in Fig. 6. After dimension reduction, four kinds of bearing samples have no intersection and overlap. Each of them has its own clustering center. The clustering effect and global geometric structure are objective, displaying a distinct manifold structure.

Fig. 6Embedded 3 dimensions by LCWR (k= 8, β= 0.1)

3.3. Dimensionality reduction effect analysis

LCWR is compared with LTSA, locally linear embedding (LLE) and principal component analysis (PCA) to verify its dimensionality reduction effect. The dimensionality reduction results of LTSA, LLE and PCA are shown in Fig. 7, Fig. 8, and Fig. 9, respectively. Generally speaking, the reduced dimensionality samples via these methods have different degrees of intersection and overlap, poor clustering within class, lack of clustering centers. It is difficult to mine the essential characteristics of the bearing state and the differences between classes. Although LTSA and LLE find the manifold structure of high-dimensional samples, they are unable to expand the gaps between dissimilar samples in neighborhood. PCA belongs to one of linear statistical distributions without considering the local structure of the samples. It makes the intraclass aggregation poor and the differences between classes unclear, which fails to reveal the non-linear manifold structure of the bearing state as shown in Fig. 9.

Fig. 7Embedded 3 dimensions by LTSA (k= 9)

Fig. 8Embedded 3 dimensions by LLE (k= 5)

Combining its weight coefficient with local coordinate’s permutation, LCWR enhances the intraclass aggregation and the differences between classes, overcomes the shortcomings that LTSA and LLE can’t enlarge the gaps between dissimilar samples in neighborhood, simplifies the dimension while retaining the low-dimensional principal characteristics of high-dimensional samples, accurately reflects the relationship between signal characteristics and the bearing state, and effectively distinguishes four kinds of bearing states. As shown in Table 2, the features extracted by various methods are sent to SVM and the recognition rate of LCWR reaches the highest 96.43 % despite a little more time to run LCWR than LTSA, LLE and PCA as shown in Table 3. Therefore, in contrast to other dimensionality reduction methods, LCWR can achieve higher accuracy and prove its effectiveness.

Fig. 9Embedded 3 dimensions by PCA

Besides, it can be founded that recognition rate of manifold dimensionality reduction using LCWR, LTSA, LLE and PCA (all greater than 90 %) is higher than that of the non-dimensionality reduction method (only 84.72 %). It is further proved that these manifold learning methods can filter redundancy and interference of the high-dimensional features and extract the intrinsic low-dimensional manifold characteristics, which can significantly improve the recognition rate of the bearing state as shown in Table 2. Meanwhile, these manifold learning methods except LCWR consume less time and get better recognition efficiency.