# Bispectrum Texture Feature Manifold for Feature Extraction in Rolling Bear Fault Diagnosis

Authors: Fei Wang; Liqing Fang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Effectively classify the fault types and the degradation degree of a rolling bearing is an important basis for accurate malfunction detection. A novel feature extract method - bispectrum image texture features manifold (BTM) of the rolling bearing vibration signal is proposed in this paper. The BTM method is realized by three main steps: bispectrum image analysis, texture feature construction and manifold feature dimensionality reduction. In this method, bispectrum analysis is employed to convert the mass vibration signals into bispectrum contour map, the typical texture features were extracted from the contour map by gray level co-occurrence matrix (GLCM), then the manifold dimensionality reduction method liner local tangent space alignment (LLTSA) is used to remove redundant information and reduce the dimension from the extracted texture features and obtain more meaningful low-dimensional information. Furthermore, the low-dimensional texture features were identified by support vector machine (SVM) which was optimized by genetic optimization algorithm (GA). The validity of BTM is confirmed by rolling bear experiments, the result show that the proposed feature extraction method can accurately distinguish different fault types and have a good performance to classify the degradation degree of inner race fault, outer race fault and rolling ball fault.


---

## Body

## 1. Introduction

The analysis of vibration signals is the most commonly used method for condition monitoring and fault diagnosis of rotating machinery, vibration signal contains a large number of state changes information in the process of mechanical operation, extracting fault features from the state information is the key technology of mechanical condition monitoring and diagnosis. Currently, mechanical vibration signal processing techniques are mainly based on the time-domain method, the frequency domain method and the time-frequency domain method. When the failure occurred, the vibration signal generated by rotating machinery contain complex components, including both non-gaussian distribution components (gear meshing frequency, rolling bearing characteristic frequency, or harmonic components) and gaussian distribution components (strong background noise). Therefore, these vibration signals presents nonlinear, non-stationary and non-gaussian, the stationary signal processing method such as FFT cannot meet the processing requirements, classical power spectrum analysis [1] and other time-frequency analysis methods (Short Time Fourier [2], Wigner-Ville distribution [3], Continuous Wavelet Transform) often lose amplitude of the signal phase information and can not handle non-minimum phase systems or non-Gaussian distributed signals [4]. In order to solve these problems, a good signal processing method is needed.High-order statistics (HOS), pointed out by Nikias, Mendel and Swami [5–7], have established a status as a sophisticated mathematical and signal processing tool for nonlinear system analysis. Higher order spectra are useful signal processing tools that have shown significant benefits over traditional power spectrum analysis, because HOS have unique properties of nonlinear system identification, Gaussian noise elimination and phase information retention [1, 8, 9]. The method of spectrum analysis is a nonlinear signal processing method based on higher-order statistics, the bispectrum is the third order spectrum and it results in a frequency-frequency-amplitude relationship which shows coupling effects between signals at different frequencies [10]. The bispectrum is easier to calculate and contains all the properties of HOS, at the same time it can completely suppress the gaussian noise theoretically and describe the nonlinear phase coupling of vibration signal which is closely related to the fault, because of these advantages it is widely used in non-stationary and non-Gaussian vibration signals processing. Reference [11] presents a novel gear wear monitoring method through a modulation signal bispectrum based motor current signal analysis, this method uses frequency changes to monitor fault conditions. Reference [12] proposed a method based on horizontal slice of cyclic bispectrum to analysis the vibration signals of rolling element bearings. Reference [13] used the technique of spectral phase analysis to explain the feature spectrum extracted from the gearbox fault data. In general, these papers only utilize the partial properties of bispectrum for feature extraction, for instance, frequency changes, horizontal slice or spectral phase. Preliminary research showed bispectrum contour map contains a lot of fault information, and the difference between the image of the bispectrum has a certain relationship with the fault status, more comprehensive fault information can be obtained by extracting the feature of spectrum image, now few papers have done research on the bispectrum image feature extraction of machinery.Texture is an inherent property of entities or scenes, which has the characteristics of brightness, color, shape, scale, etc. Texture analysis aiming to interpret and understand real-world visual patterns is an active and challenging research field [14]. Now, image texture-based feature extraction has been widely used in medical image recognition, facial recognition, finger printing and other fields [15–19], but few researches on mechanical fault recognition. Texture features are remarkable property in bispectrums, texture features of the bispectrum are very different under different fault conditions. Based on this, different fault states can be identified by extracting the quantized texture features. Gray-scale concurrence matrix (GLCM), also called spatial grayscale matrix, is a texture feature description algorithm proposed by Haralick in 1973 [20]. GLCM is a statistical theory-based algorithm that has a simple principle and a full description of image texture features. Using GLCM to extract the texture features in the bispectrum, the feature extraction of one-dimensional vibration signal is transformed into two-dimensional image texture feature extraction,with the increase of dimensionality, more fault information is obtained.As the GLCM has many characteristic parameters, the texture feature data set extracted from bispectrum are non-linear and may contain some redundant information, which is a challenge for classification and identification of fault features. Therefore, the feature dataset needs to be simplified. Manifold learning [21–24] is a new and effective method of non-linear dimensionality reduction, it has attracted more research attention in machinery, the purpose of manifold learning methods is to project the original high-dimensional data onto a lower dimension feature space with the local neighborhood structure preserved. The representative manifold learning methods contain linear embedding (LLE) [25], locality preserving projection (LPP) [26]. In this paper, a novel manifold learning method called liner local tangent space alignment (LLTSA) [27] is used for analyzing bispectrum texture feature set. Compared with PCA, LLE and LPP, LLTSA has a superior clustering performance and nonlinear complex information processing capabilities [27–29].In summary, combining bispectrum analysis with GLCM texture features can obtain better feature extraction effects, and the LLTSA method can further reduce the dimensions of feature vector, based on the above advantages, we propose a new fault feature extraction method based on bispectrum, GLCM and LLTSA. Firstly, the bispectrum method is used to obtain the contour map which contain signature information of vibration signal. Secondly, the texture feature from bispectrum contour map is extracted by GLCM, the texture features are marked as fault features. Thirdly, the LLTSA is applied to reduce the dimension for the fault features. Finally, SVM based on GA is established to achieve the fault diagnosis for machinery. The examples with new fault feature extraction method showed this method has a good performance in rolling bear fault diagnosis.The remainder of the paper is organized as following. In Section2, the theory and methods of bispectrum are reviewed briefly and extracting GLCM texture feature from bispectrum are described in detail. The proposed method is then verified by applications to bearing defect identification in Section 3, respectively. And finally the conclusion is given in Section 4.


## 2. Technical Background

Bispectrum, GLCM, Texture analysis and LLTSA have been widely applied in many fields and are the basis of the BTM method introduced. The theoretical bases of these methods are briefly introduced in Sections2.1–2.3, and the process of the proposed method are presented in Section 2.4.


### 2.1. Bispectrum

From a theoretical point of view, bispectrum is one kind of high-order statistics (HOS), first the higher order statistics are introduced as follows:Setx(n) as discrete zero mean stationary random process, its kth cumulant is defined as(1)ckxτ1,…,τk-1=cum⁡xn,xn+τ1,…,xn+τk-1Ifckxτ1,…,τk-1 is absolutely summable, or ∑τ1=-∞∞⋯∑τk-1=-∞∞ckxτ1,…,τk-1<∞,Thekth order cumulant spectrum is defined as the discrete Fourier transform (DFT) of kth the order cumulant:(2)Skxω1,…,ωk-1=∑τ1=-∞∞⋯∑τk-1=-∞∞ckxτ1,…,τk-1exp⁡-j∑i=1k-1ωiτi<∞The higher order spectrum is the higher order cumulant spectrum described above.Whenk=3, the formula (2) is defined as the bispectrum, it can be considered as a two-dimensional Fourier transform of the third order statistics of the signal x(n):(3)Bω1,ω2=∑τ1=-∞∞∑τ2=-∞∞c3xτ1,τ2e-jτ1ω1+τ2ω2The bispectru has many advantages including: (a) it can minimize the effect of Gaussian (colored or white) noise; (b) it can extract the process information in term of deviations from Gaussianity; and (c) it is able to characterize non-linear properties of process in terms of the frequency variation and time variation [10, 30, 31].


### 2.2. Gray-Scale Concurrence Matrix Texture Feature Extraction

#### 2.2.1. Gray-Scale Concurrence Matrix (GLCM)

GLCM is a quantitative description method of the image, and it is a widely used texture statistical analysis and texture measurement techniques in the field of image processing. GLCM is essentially a matrix function of pixel distance and angle, by calculating the correlation between a certain distance in the image and the grayscale of two points in a certain direction, the comprehensive information of the image in the direction, the interval, the change range and the change speed of the image is reflected. The specific calculation process of GLCM is as follows:Starting from pixel position(x,y) with image gray i, statistics the symbiotic probability function P(i,j,d,θ) between (x,y) and (x+Dx,y+Dy), among them, the distance between (x,y) and (x+Dx,y+Dy) is d and the grayscale of (x+Dx,y+Dy) is j,(4)Pi,j,d,θ=x,y,x+Dx,y+Dy∣fx,y=i,fx+Dx,y+Dy=jIn the formula,i,j=0,1,2,…,N-1, N=256, Dx,Dy are position offset, d is the generate step value and θ is generate directions of GLCM, usually, θ = 0°, 45°, 90°, 135°,Before extracting the feature statistics of GLCM, we usually normalize the matrix according to the following formula:(5)gi,j=Pi,j∑i=0N-1∑j=0N-1Pi,jSummarize the above process, the gray level co-occurrence matrix is the joint probability distribution of two gray pixels with the distance D in the image.


#### 2.2.2. Texture Feature Parameters of GLCM

Texture analysis is a technology to obtain the texture feature parameters through image processing, it is the basic for analyzing local patterns and permutations of an image. In this paper, texture features are reflected by GLCM’s partial eigenvalues, this method essentially represents the texture features with conditional probability and is the representation of the grayscale correlation of adjacent pixels in the image.GLCM has 14 characteristic statistics and 9 relevant textural parameters will be used during the diagnosis:(1) Angular Second Moment. ASM is the uniformity of a image’s gray distribution. When the image gray distribution is more uniform, the ASM is larger; on the contrary, the ASM is smaller.(6)asm=∑i=0N-1∑j=0N-1g2i,j(2) Contrast. Contrast is a measure of the clarity of the texture. The deeper the texture in the image, the greater the contrast in the image.(7)con=∑i=0N-1∑j=0N-1i-j2gi,j(3) Correlation. Correlation is a measure of the degree of similarity of GLCM elements in a row or column direction and is the linear relationship of the grayscale of the image.(8)cor=∑i=0N-1∑j=0N-1i-uj-ugi,jσ2(4) Inverse Difference Moment (IDM). IDM measures the extent of local changes in the image texture and characterizes the degree of regularity of the texture. The more regular texture, the larger the inverse torque, and vice versa.(9)idm=∑i=0N-1∑j=0N-111+i-j2gi,j(5) Entropy. He represents the amount of information in the image, which is a measure of the randomness of the image content and can characterize the complexity or non-uniformity of the texture. When the image has no texture, entropy is 0, and when texture is full, entropy is the largest.(10)entro=-∑i=0N-1∑j=0N-1gi,jlog⁡gi,j(6) Max Probability. Max Probability reflects the maximum frequency of gray-level pairs in GLCM(11)mp=maxi,j⁡gi,j(7) Mean Value. The mean value characterizes the average gray level of the image(12)mean=∑i=0N-1∑j=0N-1igi,j(8) Variance. Variance characterizes the change of image gray level(13)var=∑i=0N-1∑j=0N-1i-u2gi,j(9) Variance of the Residuals. The more obvious the gray level contrast of the image pixels, the greater the variance of the gray level differences of the adjacent pixels of the image, the worse the variance value will increase.(14)vor=∑i=0N-1i-∑k=1N-1kgx+ykgx-yi


### 2.3. Liner Local Tangent Space Alignment (LLTSA)

Linear local tangent space alignment (LLTSA) is a novel linear dimensionality reduction algorithm. It uses the tangent space in the neighborhood of a data point to represent the local geometry, and then aligns those local tangent spaces in the low-dimensional space which is linearly mapped from the raw high-dimensional space. The method can be viewed as a linear approximation of the nonlinear local tangent space alignment algorithm.The procedure of LLTSA algorithm is summarized as follows:Given a data setX=[x1,…,xN]Step 1 (principal component analysis (PCA) projection).
The PCA method is used to map the data set to the subjective subspace.APCA represents the transfer matrix of PCA, at the same time, for conveniently presentation, X is applied to represent the PCA subspace data set.Step 2 (determine the neighbor region).
The K-nearest neighbor method (KNN) is used to search for the data points’ neighbor region, that means the distance matrix of all data points is constructed by using Euclidean distance, and then thek similar points xj of the data points xi are found by analyzing the distance matrix.Step 3 (calculate local information).
Matrix ofd feature vectors corresponding to thed largest eigenvalue of determined by XiHk(Xi=[xi1,…,xik]) is calculated. Among them, Hk=I-eeT/kStep 4 (constructing arrangement matrix).
Construct matrixB by local accumulation:(15)BIi,Ii←BIi,Ii+WiWiTi=1,…,N
InitializeB=0, Ii=ii,…,ik denote the index set of k nearest neighbors for xi, Wi=Hki-ViViT(i=1,…,N).Step 5 (calculate mapping).
Calculate the eigenvalues and eigenvectors of the generalized feature problem:(16)XHNBHNXTα=λXHNXTα
The corresponding characteristic vector solution with the eigenvalueλ1,λ2,…,λn(λ1<λ2⋯<λn) is α1,α2,…,αd, ALLTSA=(α1,α2,…,αd). The ultimate transformation matrix is as follows: A=APCA×ALLTSA, X→Y=ATXORGHN.


### 2.4. Technological Process of the Proposed Method

Motived by the advantages of bispectrum analysis, the GLCM texture feature and LLTSA for feature extraction is proposed, the flowchart of this method for bearing fault diagnosis is depicted in Figure1, the work steps are as follows:Figure 1
Flow chart of the proposed method.(1) Conduct the bispectral analysis of the original vibration signal and get the contour map which reflecting the spectrum characteristics of each fault status.(2) Convert the each map to GLCM and calculate 9 typical texture features in4 directions (0∘, 45∘, 90∘, 135∘).(3) Process the texture features by LLTSA and get the low-dimensional characteristics of the essential manifold. In this step, the amplitude of target dimensionn and neighborhood parameters ε should be determined by the specific data.(4) The texture features processed in Step(3) are entered into the SVM for training and testing.Kernel functions of SVM is radial basis function (RBF), penalty parameterγ and kernel function parameters σ is optimized by genetic optimization algorithm (GA).(5) Get fault classification accuracy under different fault conditions.


## 2.1. Bispectrum

From a theoretical point of view, bispectrum is one kind of high-order statistics (HOS), first the higher order statistics are introduced as follows:Setx(n) as discrete zero mean stationary random process, its kth cumulant is defined as(1)ckxτ1,…,τk-1=cum⁡xn,xn+τ1,…,xn+τk-1Ifckxτ1,…,τk-1 is absolutely summable, or ∑τ1=-∞∞⋯∑τk-1=-∞∞ckxτ1,…,τk-1<∞,Thekth order cumulant spectrum is defined as the discrete Fourier transform (DFT) of kth the order cumulant:(2)Skxω1,…,ωk-1=∑τ1=-∞∞⋯∑τk-1=-∞∞ckxτ1,…,τk-1exp⁡-j∑i=1k-1ωiτi<∞The higher order spectrum is the higher order cumulant spectrum described above.Whenk=3, the formula (2) is defined as the bispectrum, it can be considered as a two-dimensional Fourier transform of the third order statistics of the signal x(n):(3)Bω1,ω2=∑τ1=-∞∞∑τ2=-∞∞c3xτ1,τ2e-jτ1ω1+τ2ω2The bispectru has many advantages including: (a) it can minimize the effect of Gaussian (colored or white) noise; (b) it can extract the process information in term of deviations from Gaussianity; and (c) it is able to characterize non-linear properties of process in terms of the frequency variation and time variation [10, 30, 31].


## 2.2. Gray-Scale Concurrence Matrix Texture Feature Extraction

### 2.2.1. Gray-Scale Concurrence Matrix (GLCM)

GLCM is a quantitative description method of the image, and it is a widely used texture statistical analysis and texture measurement techniques in the field of image processing. GLCM is essentially a matrix function of pixel distance and angle, by calculating the correlation between a certain distance in the image and the grayscale of two points in a certain direction, the comprehensive information of the image in the direction, the interval, the change range and the change speed of the image is reflected. The specific calculation process of GLCM is as follows:Starting from pixel position(x,y) with image gray i, statistics the symbiotic probability function P(i,j,d,θ) between (x,y) and (x+Dx,y+Dy), among them, the distance between (x,y) and (x+Dx,y+Dy) is d and the grayscale of (x+Dx,y+Dy) is j,(4)Pi,j,d,θ=x,y,x+Dx,y+Dy∣fx,y=i,fx+Dx,y+Dy=jIn the formula,i,j=0,1,2,…,N-1, N=256, Dx,Dy are position offset, d is the generate step value and θ is generate directions of GLCM, usually, θ = 0°, 45°, 90°, 135°,Before extracting the feature statistics of GLCM, we usually normalize the matrix according to the following formula:(5)gi,j=Pi,j∑i=0N-1∑j=0N-1Pi,jSummarize the above process, the gray level co-occurrence matrix is the joint probability distribution of two gray pixels with the distance D in the image.


### 2.2.2. Texture Feature Parameters of GLCM

Texture analysis is a technology to obtain the texture feature parameters through image processing, it is the basic for analyzing local patterns and permutations of an image. In this paper, texture features are reflected by GLCM’s partial eigenvalues, this method essentially represents the texture features with conditional probability and is the representation of the grayscale correlation of adjacent pixels in the image.GLCM has 14 characteristic statistics and 9 relevant textural parameters will be used during the diagnosis:(1) Angular Second Moment. ASM is the uniformity of a image’s gray distribution. When the image gray distribution is more uniform, the ASM is larger; on the contrary, the ASM is smaller.(6)asm=∑i=0N-1∑j=0N-1g2i,j(2) Contrast. Contrast is a measure of the clarity of the texture. The deeper the texture in the image, the greater the contrast in the image.(7)con=∑i=0N-1∑j=0N-1i-j2gi,j(3) Correlation. Correlation is a measure of the degree of similarity of GLCM elements in a row or column direction and is the linear relationship of the grayscale of the image.(8)cor=∑i=0N-1∑j=0N-1i-uj-ugi,jσ2(4) Inverse Difference Moment (IDM). IDM measures the extent of local changes in the image texture and characterizes the degree of regularity of the texture. The more regular texture, the larger the inverse torque, and vice versa.(9)idm=∑i=0N-1∑j=0N-111+i-j2gi,j(5) Entropy. He represents the amount of information in the image, which is a measure of the randomness of the image content and can characterize the complexity or non-uniformity of the texture. When the image has no texture, entropy is 0, and when texture is full, entropy is the largest.(10)entro=-∑i=0N-1∑j=0N-1gi,jlog⁡gi,j(6) Max Probability. Max Probability reflects the maximum frequency of gray-level pairs in GLCM(11)mp=maxi,j⁡gi,j(7) Mean Value. The mean value characterizes the average gray level of the image(12)mean=∑i=0N-1∑j=0N-1igi,j(8) Variance. Variance characterizes the change of image gray level(13)var=∑i=0N-1∑j=0N-1i-u2gi,j(9) Variance of the Residuals. The more obvious the gray level contrast of the image pixels, the greater the variance of the gray level differences of the adjacent pixels of the image, the worse the variance value will increase.(14)vor=∑i=0N-1i-∑k=1N-1kgx+ykgx-yi


## 2.2.1. Gray-Scale Concurrence Matrix (GLCM)

GLCM is a quantitative description method of the image, and it is a widely used texture statistical analysis and texture measurement techniques in the field of image processing. GLCM is essentially a matrix function of pixel distance and angle, by calculating the correlation between a certain distance in the image and the grayscale of two points in a certain direction, the comprehensive information of the image in the direction, the interval, the change range and the change speed of the image is reflected. The specific calculation process of GLCM is as follows:Starting from pixel position(x,y) with image gray i, statistics the symbiotic probability function P(i,j,d,θ) between (x,y) and (x+Dx,y+Dy), among them, the distance between (x,y) and (x+Dx,y+Dy) is d and the grayscale of (x+Dx,y+Dy) is j,(4)Pi,j,d,θ=x,y,x+Dx,y+Dy∣fx,y=i,fx+Dx,y+Dy=jIn the formula,i,j=0,1,2,…,N-1, N=256, Dx,Dy are position offset, d is the generate step value and θ is generate directions of GLCM, usually, θ = 0°, 45°, 90°, 135°,Before extracting the feature statistics of GLCM, we usually normalize the matrix according to the following formula:(5)gi,j=Pi,j∑i=0N-1∑j=0N-1Pi,jSummarize the above process, the gray level co-occurrence matrix is the joint probability distribution of two gray pixels with the distance D in the image.


## 2.2.2. Texture Feature Parameters of GLCM

Texture analysis is a technology to obtain the texture feature parameters through image processing, it is the basic for analyzing local patterns and permutations of an image. In this paper, texture features are reflected by GLCM’s partial eigenvalues, this method essentially represents the texture features with conditional probability and is the representation of the grayscale correlation of adjacent pixels in the image.GLCM has 14 characteristic statistics and 9 relevant textural parameters will be used during the diagnosis:(1) Angular Second Moment. ASM is the uniformity of a image’s gray distribution. When the image gray distribution is more uniform, the ASM is larger; on the contrary, the ASM is smaller.(6)asm=∑i=0N-1∑j=0N-1g2i,j(2) Contrast. Contrast is a measure of the clarity of the texture. The deeper the texture in the image, the greater the contrast in the image.(7)con=∑i=0N-1∑j=0N-1i-j2gi,j(3) Correlation. Correlation is a measure of the degree of similarity of GLCM elements in a row or column direction and is the linear relationship of the grayscale of the image.(8)cor=∑i=0N-1∑j=0N-1i-uj-ugi,jσ2(4) Inverse Difference Moment (IDM). IDM measures the extent of local changes in the image texture and characterizes the degree of regularity of the texture. The more regular texture, the larger the inverse torque, and vice versa.(9)idm=∑i=0N-1∑j=0N-111+i-j2gi,j(5) Entropy. He represents the amount of information in the image, which is a measure of the randomness of the image content and can characterize the complexity or non-uniformity of the texture. When the image has no texture, entropy is 0, and when texture is full, entropy is the largest.(10)entro=-∑i=0N-1∑j=0N-1gi,jlog⁡gi,j(6) Max Probability. Max Probability reflects the maximum frequency of gray-level pairs in GLCM(11)mp=maxi,j⁡gi,j(7) Mean Value. The mean value characterizes the average gray level of the image(12)mean=∑i=0N-1∑j=0N-1igi,j(8) Variance. Variance characterizes the change of image gray level(13)var=∑i=0N-1∑j=0N-1i-u2gi,j(9) Variance of the Residuals. The more obvious the gray level contrast of the image pixels, the greater the variance of the gray level differences of the adjacent pixels of the image, the worse the variance value will increase.(14)vor=∑i=0N-1i-∑k=1N-1kgx+ykgx-yi


## 2.3. Liner Local Tangent Space Alignment (LLTSA)

Linear local tangent space alignment (LLTSA) is a novel linear dimensionality reduction algorithm. It uses the tangent space in the neighborhood of a data point to represent the local geometry, and then aligns those local tangent spaces in the low-dimensional space which is linearly mapped from the raw high-dimensional space. The method can be viewed as a linear approximation of the nonlinear local tangent space alignment algorithm.The procedure of LLTSA algorithm is summarized as follows:Given a data setX=[x1,…,xN]Step 1 (principal component analysis (PCA) projection).
The PCA method is used to map the data set to the subjective subspace.APCA represents the transfer matrix of PCA, at the same time, for conveniently presentation, X is applied to represent the PCA subspace data set.Step 2 (determine the neighbor region).
The K-nearest neighbor method (KNN) is used to search for the data points’ neighbor region, that means the distance matrix of all data points is constructed by using Euclidean distance, and then thek similar points xj of the data points xi are found by analyzing the distance matrix.Step 3 (calculate local information).
Matrix ofd feature vectors corresponding to thed largest eigenvalue of determined by XiHk(Xi=[xi1,…,xik]) is calculated. Among them, Hk=I-eeT/kStep 4 (constructing arrangement matrix).
Construct matrixB by local accumulation:(15)BIi,Ii←BIi,Ii+WiWiTi=1,…,N
InitializeB=0, Ii=ii,…,ik denote the index set of k nearest neighbors for xi, Wi=Hki-ViViT(i=1,…,N).Step 5 (calculate mapping).
Calculate the eigenvalues and eigenvectors of the generalized feature problem:(16)XHNBHNXTα=λXHNXTα
The corresponding characteristic vector solution with the eigenvalueλ1,λ2,…,λn(λ1<λ2⋯<λn) is α1,α2,…,αd, ALLTSA=(α1,α2,…,αd). The ultimate transformation matrix is as follows: A=APCA×ALLTSA, X→Y=ATXORGHN.


## 2.4. Technological Process of the Proposed Method

Motived by the advantages of bispectrum analysis, the GLCM texture feature and LLTSA for feature extraction is proposed, the flowchart of this method for bearing fault diagnosis is depicted in Figure1, the work steps are as follows:Figure 1
Flow chart of the proposed method.(1) Conduct the bispectral analysis of the original vibration signal and get the contour map which reflecting the spectrum characteristics of each fault status.(2) Convert the each map to GLCM and calculate 9 typical texture features in4 directions (0∘, 45∘, 90∘, 135∘).(3) Process the texture features by LLTSA and get the low-dimensional characteristics of the essential manifold. In this step, the amplitude of target dimensionn and neighborhood parameters ε should be determined by the specific data.(4) The texture features processed in Step(3) are entered into the SVM for training and testing.Kernel functions of SVM is radial basis function (RBF), penalty parameterγ and kernel function parameters σ is optimized by genetic optimization algorithm (GA).(5) Get fault classification accuracy under different fault conditions.


## 3. Experimental Data Analysis

In this section, rolling bearing fault data is processed according to the proposed diagnosis flow chart procedure to verify the validity of the proposed method. The test range includes the identification of different fault types and degradation degrees of rolling bearing include outer race fault, inner race fault with rolling elements fault.


### 3.1. Acquire Fault Data of Rolling Bearing

Rolling bearing experiment data are offered by Rolling Bearing Data Center of Case Western Reserve University. The rolling bearing test-bed is showed in Figure2. The test bearings are 6205-2RSJEM SKF deep groove ball bearings. Single point fault was arranged to the test bearings using electro-discharge machining with fault diameters of 7 mils, 14mils, 21 mils, 28 mils (1mils=0.001inch= 0.00254cm). The data in the paper is the drive end accelerometer data. The vibration acceleration signals from four kinds rolling bearings including healthy bearing (H) and bearings with outer race fault (OF), inner race fault (IF) and rolling elements fault (BF) were collected under the speed 1797 r/min of the rotor, with the sampling frequency 12 kHz and each group of data contains 2000 points. The rolling bearing fault data were showed in Table 1.Table 1
The rolling bearing fault data used in this paper.
Rolling bearing type Rotor speed Sampling frequency H OF (inch) IF (inch) BF (inch) 6205-2RSJEM SKF 1797(r/min) 12000Hz Normaldata 0.007 0.007 0.007 0.014 0.014 0.014 0.021 0.021 0.021 0.028 0.028Figure 2
The rolling bearing test-bed.


### 3.2. Bispectrum Analysis of Rolling Bearing

Figure3 is the rolling bearing vibration signals of variety fault conditions, due to the presence of a large number of noise and other complex components, the morphological of each signal is diversity, it can not directly distinguish the various forms of signal fault by analyzing the vibration signal. The vibration signals in all fault states were processed by bispectrum analysis, bispectral contour maps from the analysis result were shown in Figure 4, the maps of each fault state are quite different in texture characteristics and each image’s feature is unique. All kinds of fault conditions would be identified while accurately extract the texture features of each image.Figure 3
The rolling bearing vibration signal of healthy, outer race fault (OF), inner race fault (IF) and ball fault (BF).Figure 4
The rolling bearing bispectral contour map of healthy, outer race fault (OF), inner race fault (IF) and ball fault (BF).


### 3.3. Image Texture Feature Extraction

Nine kinds of texture feature parameters of 4 directions (0°, 45°, 90°, 135°) under each fault condition are calculated by GLCM and a 36-dimensional texture feature matrix is established for each fault conductions. To illustrate the distribution of texture features with different fault status, Figure5 shows part of the texture feature parameters’ distribution under different fault conditions, it can be known from the figure that these texture parameters have better classification performance, but in some states, the values of the parameters are relatively stable, but in some fault conductions there exists serious distribution alias. Using part of the texture feature parameters can distinguish the fault state to a limited extent, so it needs to calculate the texture feature parameters of multiple directions and get a better distribution effect.Figure 5
Part texture feature parameters distribution in 0° direction with different fault conductions. (a) Angular Second Moment (b) Variance (c) Correlation (d) Contrast.
(a) (b) (c) (d)


### 3.4. Texture Feature Matrix Dimension Reduction Using LLTSA

LLTSA is used to reduce the dimensionality of the extracted texture features, to investigate the dimensionality reduction effectiveness of this method, the 36-dimensional feature reduces to 3 dimensions by LLTSA, the neighborhood parametersε is set to 12. The distribution of 3 dimensions features (FE1, FE2, FE3) under different fault status are plotted in Figure 6. From the result, we can find LLTSA has good feature classification performance on the whole, Figures 6(a) and 6(c) have a better performance in fault sample clustering and can clearly distinguish different fault status, due to the influence of aliasing, the classification effect of fault status in Figures 6(b) and 6(d) are not very well, this shows that the effect of using unified parameter settings for fault classification is not very well. To achieve a good dimension reduction effect, it needs to find the best parameters reduction dimension n and neighborhood parameter ε.Figure 6
Dimension reduction results of LLTSA in different fault conductions.
(a) (b) (c) (d)In order to study the relationship between classification accuracy and parameter selection, LIBSVM is used for classification and identification, and genetic algorithm(GA) is used to optimize SVM parameters, the population size is set to 20, the maximum number of iterations is 200, and the search range of the SVM penalty parameterC and the kernel function parameter σ are both 0 to 500; 100 data sets are used for each fault condition, among them,20 data sets are used to train the GA-SVM, 80 data sets are used for testing. Figure 7 shows the relationship between the classification accuracy and the two parameter settings in distinguish health, inner race fault, outer race fault and rolling element fault, the highest recognition rate corresponding to the best parameters. through this way we can get the optimal parameter settings in fault recognition for other fault recognition situations. Table 2 shows the rolling bearing fault diagnosis classification results under different fault conditions by GA-SVM. The total fault recognition rate is more than 95%, in some cases the recognition rate is close to 100%. On the one hand, this method can accurately identify different fault types, on the other hand, it has excellent diagnostic capability for the degree of fault. The experiment results show the method proposed has a good performance in fault diagnosis of rolling bearingsTable 2
The classification accuracy using the proposed method.
Fault type Accuracy (%) Mean (%) optimal parameter Fault type Accuracy (%) Mean (%) optimal parameter Healthy 100 99.69 n = 19  ε=22 Healthy 100 99.75 n = 13  ε=18 IF(0.007) 100 IF(0.007) 98.75 OF(0.007) 100 IF(0.014) 100 BF(0.007) 98.75 IF(0.021) 100 IF(0.028) 100  Fault type accuracy (%) Mean (%) optimal parameter Fault type accuracy (%) Mean (%) optimal parameter  Healthy 100 95.6 n = 7  ε=26 Healthy 100 90.5 n = 9  ε=23 OF(0.007) 93.75 BF(0.007) 100 OF(0.014) 91.25 BF(0.014) 87.5 OF(0.021) 97.5 BF(0.021) 82.5 BF(0.028) 82.5Figure 7
The relationship of the fault classification accuracy vs. parametern and ε.In order to show the advantages of LLTSA in data fusion and dimensionality reduction, The principal component analysis (PCA), Locally Linear Embedding (LLE), Locality Preserving Projections (LPP) and Local Tangent Space Alignment (LTSA) methods are used for data dimension reduction and the parameter settings are the same as the LLTSA method. Case 1 is the classification accuracy for health, inner race fault, outer race fault and rolling element fault, case 2 is inner race fault severity identification accuracy, case 3 is outer race fault severity identification accuracy, case 4 is rolling element fault severity identification accuracy, the result of fault identification accuracy under different cases are shown in Table4.Table3 shows that the LLTSA method achieves the highest fault recognition accuracy in four cases, the manifold learning algorithm including LEE, LPP, LTSA and LLTSA have obvious advantages compared with the linear dimensionality reduction method PCA as the extracted high-dimensional texture feature data present locality and nonlinear. Based on the LTSA algorithm and the linear block method, the LLTSA is proposed to solve the problem which LTSA is not suitable for dealing with high curvature distribution and sparse heterogeneous distribution data and it has better dimensionality reduction performance than other manifold learning methods. The final recognition result also shows that the feature extraction method proposed in this paper can accurately identify the fault location and the fault degeneration degree of the rolling bearing.Table 3
Identification accuracy of four cases under different dimension reduction algorithms.
identification accuracy (%) PCA LLE LPP LTSA LLTSA Case 1 79.68 98.44 95.63 92.19 99.69  Case 2 74.25 98.25 99.75 97.75 99.75  Case 3 76.25 81.88 88.44 91.25 95.6  Case 4 23.2 79.2 80 80.8 90.5Table 4
EEMD-energy feature result of different fault types.
Fault V 1 V 2 V 3 V 4 V 5 Healthy 0.7107 0.0349 0.1639 0.0655 0.0250  IF(0.007) 0.5167 0.3980 0.0657 0.0159 0.0038  IF(0.014) 0.5199 0.4487 0.0153 0.0087 0.0074  IF(0.021) 0.5060 0.4702 0.0079 0.0120 0.0038  IF(0.028) 0.5314 0.3278 0.0585 0.0803 0.0020  BF(0.007) 0.5089 0.4653 0.0109 0.0090 0.0059  BF(0.014) 0.5111 0.4572 0.0071 0.0062 0.0184  BF(0.021) 0.5005 0.4773 0.0100 0.0082 0.0041  BF(0.028) 0.5005 0.4927 0.0051 0.0012 0.0001  OF(0.007) 0.4980 0.4902 0.0074 0.003 0.002  OF(0.014) 0.5086 0.4357 0.0171 0.0244 0.0141


### 3.5. Comparison with EEMD Feature Extraction Method

EEMD method was first introduced by Huang et al [32], it is a self -adaptive decomposition method and have a good performance for feature extraction in the machinery vibration signals. In order to illustrate the advantages of the feature extraction methods mentioned in this paper, we compared it with the EEMD-energy feature method, this method extracting energy values of different frequency bands as feature values for fault diagnosis. Firstly, the bearing vibration signal is decomposed to several IMFs by EEMD, then we selected the first five IMF components based on the principle of correlation coefficient and calculated the energy characteristics. The method for calculating the energy of each IMF is as follows:(17)Ei=∑i=1kfit2E=∑i=1nEiVi=EiEwhere Ei is the energy ofith IMF, E is the sum of energies of n IMFs and Vi represents the percentage of energy ofith IMF in the whole signal energyE. Figure 8 is the EEMD result for inner fault bearing vibration signals, Table 4 is the five IMF energy feature value of different faults.Figure 8
Results of EEMD for inner fault (7mils) bearing vibration signals.The energy values of the five extracted IMFs are used as feature vectors for pattern recognition by GA-SVM, 20 sets of data are used as training samples, and 80 sets of data are used as test samples. The recognition result are showed in Table5. Table 5 showed the EEMD-energy method have a good performance in distinguish fault types compared with distinguish the fault degree, the proposed method is superior to the EEMD-energy method in distinguish the fault types and fault degrees.Table 5
Fault identification comparison result.
identification accuracy EEMD- energy method Proposed method Case 1 99.37% 99.69%  Case 2 81.87% 99.75%  Case 3 92.00% 95.60%  Case 4 81.87% 90.50%


## 3.1. Acquire Fault Data of Rolling Bearing

Rolling bearing experiment data are offered by Rolling Bearing Data Center of Case Western Reserve University. The rolling bearing test-bed is showed in Figure2. The test bearings are 6205-2RSJEM SKF deep groove ball bearings. Single point fault was arranged to the test bearings using electro-discharge machining with fault diameters of 7 mils, 14mils, 21 mils, 28 mils (1mils=0.001inch= 0.00254cm). The data in the paper is the drive end accelerometer data. The vibration acceleration signals from four kinds rolling bearings including healthy bearing (H) and bearings with outer race fault (OF), inner race fault (IF) and rolling elements fault (BF) were collected under the speed 1797 r/min of the rotor, with the sampling frequency 12 kHz and each group of data contains 2000 points. The rolling bearing fault data were showed in Table 1.Table 1
The rolling bearing fault data used in this paper.
Rolling bearing type Rotor speed Sampling frequency H OF (inch) IF (inch) BF (inch) 6205-2RSJEM SKF 1797(r/min) 12000Hz Normaldata 0.007 0.007 0.007 0.014 0.014 0.014 0.021 0.021 0.021 0.028 0.028Figure 2
The rolling bearing test-bed.


## 3.2. Bispectrum Analysis of Rolling Bearing

Figure3 is the rolling bearing vibration signals of variety fault conditions, due to the presence of a large number of noise and other complex components, the morphological of each signal is diversity, it can not directly distinguish the various forms of signal fault by analyzing the vibration signal. The vibration signals in all fault states were processed by bispectrum analysis, bispectral contour maps from the analysis result were shown in Figure 4, the maps of each fault state are quite different in texture characteristics and each image’s feature is unique. All kinds of fault conditions would be identified while accurately extract the texture features of each image.Figure 3
The rolling bearing vibration signal of healthy, outer race fault (OF), inner race fault (IF) and ball fault (BF).Figure 4
The rolling bearing bispectral contour map of healthy, outer race fault (OF), inner race fault (IF) and ball fault (BF).


## 3.3. Image Texture Feature Extraction

Nine kinds of texture feature parameters of 4 directions (0°, 45°, 90°, 135°) under each fault condition are calculated by GLCM and a 36-dimensional texture feature matrix is established for each fault conductions. To illustrate the distribution of texture features with different fault status, Figure5 shows part of the texture feature parameters’ distribution under different fault conditions, it can be known from the figure that these texture parameters have better classification performance, but in some states, the values of the parameters are relatively stable, but in some fault conductions there exists serious distribution alias. Using part of the texture feature parameters can distinguish the fault state to a limited extent, so it needs to calculate the texture feature parameters of multiple directions and get a better distribution effect.Figure 5
Part texture feature parameters distribution in 0° direction with different fault conductions. (a) Angular Second Moment (b) Variance (c) Correlation (d) Contrast.
(a) (b) (c) (d)


## 3.4. Texture Feature Matrix Dimension Reduction Using LLTSA

LLTSA is used to reduce the dimensionality of the extracted texture features, to investigate the dimensionality reduction effectiveness of this method, the 36-dimensional feature reduces to 3 dimensions by LLTSA, the neighborhood parametersε is set to 12. The distribution of 3 dimensions features (FE1, FE2, FE3) under different fault status are plotted in Figure 6. From the result, we can find LLTSA has good feature classification performance on the whole, Figures 6(a) and 6(c) have a better performance in fault sample clustering and can clearly distinguish different fault status, due to the influence of aliasing, the classification effect of fault status in Figures 6(b) and 6(d) are not very well, this shows that the effect of using unified parameter settings for fault classification is not very well. To achieve a good dimension reduction effect, it needs to find the best parameters reduction dimension n and neighborhood parameter ε.Figure 6
Dimension reduction results of LLTSA in different fault conductions.
(a) (b) (c) (d)In order to study the relationship between classification accuracy and parameter selection, LIBSVM is used for classification and identification, and genetic algorithm(GA) is used to optimize SVM parameters, the population size is set to 20, the maximum number of iterations is 200, and the search range of the SVM penalty parameterC and the kernel function parameter σ are both 0 to 500; 100 data sets are used for each fault condition, among them,20 data sets are used to train the GA-SVM, 80 data sets are used for testing. Figure 7 shows the relationship between the classification accuracy and the two parameter settings in distinguish health, inner race fault, outer race fault and rolling element fault, the highest recognition rate corresponding to the best parameters. through this way we can get the optimal parameter settings in fault recognition for other fault recognition situations. Table 2 shows the rolling bearing fault diagnosis classification results under different fault conditions by GA-SVM. The total fault recognition rate is more than 95%, in some cases the recognition rate is close to 100%. On the one hand, this method can accurately identify different fault types, on the other hand, it has excellent diagnostic capability for the degree of fault. The experiment results show the method proposed has a good performance in fault diagnosis of rolling bearingsTable 2
The classification accuracy using the proposed method.
Fault type Accuracy (%) Mean (%) optimal parameter Fault type Accuracy (%) Mean (%) optimal parameter Healthy 100 99.69 n = 19  ε=22 Healthy 100 99.75 n = 13  ε=18 IF(0.007) 100 IF(0.007) 98.75 OF(0.007) 100 IF(0.014) 100 BF(0.007) 98.75 IF(0.021) 100 IF(0.028) 100  Fault type accuracy (%) Mean (%) optimal parameter Fault type accuracy (%) Mean (%) optimal parameter  Healthy 100 95.6 n = 7  ε=26 Healthy 100 90.5 n = 9  ε=23 OF(0.007) 93.75 BF(0.007) 100 OF(0.014) 91.25 BF(0.014) 87.5 OF(0.021) 97.5 BF(0.021) 82.5 BF(0.028) 82.5Figure 7
The relationship of the fault classification accuracy vs. parametern and ε.In order to show the advantages of LLTSA in data fusion and dimensionality reduction, The principal component analysis (PCA), Locally Linear Embedding (LLE), Locality Preserving Projections (LPP) and Local Tangent Space Alignment (LTSA) methods are used for data dimension reduction and the parameter settings are the same as the LLTSA method. Case 1 is the classification accuracy for health, inner race fault, outer race fault and rolling element fault, case 2 is inner race fault severity identification accuracy, case 3 is outer race fault severity identification accuracy, case 4 is rolling element fault severity identification accuracy, the result of fault identification accuracy under different cases are shown in Table4.Table3 shows that the LLTSA method achieves the highest fault recognition accuracy in four cases, the manifold learning algorithm including LEE, LPP, LTSA and LLTSA have obvious advantages compared with the linear dimensionality reduction method PCA as the extracted high-dimensional texture feature data present locality and nonlinear. Based on the LTSA algorithm and the linear block method, the LLTSA is proposed to solve the problem which LTSA is not suitable for dealing with high curvature distribution and sparse heterogeneous distribution data and it has better dimensionality reduction performance than other manifold learning methods. The final recognition result also shows that the feature extraction method proposed in this paper can accurately identify the fault location and the fault degeneration degree of the rolling bearing.Table 3
Identification accuracy of four cases under different dimension reduction algorithms.
identification accuracy (%) PCA LLE LPP LTSA LLTSA Case 1 79.68 98.44 95.63 92.19 99.69  Case 2 74.25 98.25 99.75 97.75 99.75  Case 3 76.25 81.88 88.44 91.25 95.6  Case 4 23.2 79.2 80 80.8 90.5Table 4
EEMD-energy feature result of different fault types.
Fault V 1 V 2 V 3 V 4 V 5 Healthy 0.7107 0.0349 0.1639 0.0655 0.0250  IF(0.007) 0.5167 0.3980 0.0657 0.0159 0.0038  IF(0.014) 0.5199 0.4487 0.0153 0.0087 0.0074  IF(0.021) 0.5060 0.4702 0.0079 0.0120 0.0038  IF(0.028) 0.5314 0.3278 0.0585 0.0803 0.0020  BF(0.007) 0.5089 0.4653 0.0109 0.0090 0.0059  BF(0.014) 0.5111 0.4572 0.0071 0.0062 0.0184  BF(0.021) 0.5005 0.4773 0.0100 0.0082 0.0041  BF(0.028) 0.5005 0.4927 0.0051 0.0012 0.0001  OF(0.007) 0.4980 0.4902 0.0074 0.003 0.002  OF(0.014) 0.5086 0.4357 0.0171 0.0244 0.0141


## 3.5. Comparison with EEMD Feature Extraction Method

EEMD method was first introduced by Huang et al [32], it is a self -adaptive decomposition method and have a good performance for feature extraction in the machinery vibration signals. In order to illustrate the advantages of the feature extraction methods mentioned in this paper, we compared it with the EEMD-energy feature method, this method extracting energy values of different frequency bands as feature values for fault diagnosis. Firstly, the bearing vibration signal is decomposed to several IMFs by EEMD, then we selected the first five IMF components based on the principle of correlation coefficient and calculated the energy characteristics. The method for calculating the energy of each IMF is as follows:(17)Ei=∑i=1kfit2E=∑i=1nEiVi=EiEwhere Ei is the energy ofith IMF, E is the sum of energies of n IMFs and Vi represents the percentage of energy ofith IMF in the whole signal energyE. Figure 8 is the EEMD result for inner fault bearing vibration signals, Table 4 is the five IMF energy feature value of different faults.Figure 8
Results of EEMD for inner fault (7mils) bearing vibration signals.The energy values of the five extracted IMFs are used as feature vectors for pattern recognition by GA-SVM, 20 sets of data are used as training samples, and 80 sets of data are used as test samples. The recognition result are showed in Table5. Table 5 showed the EEMD-energy method have a good performance in distinguish fault types compared with distinguish the fault degree, the proposed method is superior to the EEMD-energy method in distinguish the fault types and fault degrees.Table 5
Fault identification comparison result.
identification accuracy EEMD- energy method Proposed method Case 1 99.37% 99.69%  Case 2 81.87% 99.75%  Case 3 92.00% 95.60%  Case 4 81.87% 90.50%


## 4. Conclusions

In this paper, vibration image recognition technology is applied to fault diagnosis, it combines bispectrum analysis, GLCM texture features and LLTSA data dimension reduction method. This method is an exploratory study of fault feature extraction method. Comparison to other manifold learning dimensionality reduction method and EEMD-energy feature extraction method showed the advantages of this method: this method has a better performance in distinguish faults locations and the fault degeneration degree of rolling bearing. But the method still have problems that need to be solved: LLTSA’s parameters determination takes a lot of calculation time, which is not conducive to condition monitoring and fault diagnosis. Using adaptive methods to select suitable parametern and ε is the focus of the next step in the research.



---

*Source: 3805729-2019-02-26.xml*