Abstract
In this paper, a fault diagnosis method based on symmetric polar coordinate image and Fuzzy C-Means clustering algorithm is proposed to solve the problem that the fault signal of axial piston pump is not intuitive under the time-domain waveform diagram. In this paper, the sampled vibration signals of axial piston pump were denoised firstly by the combination of ensemble empirical mode decomposition and Pearson correlation coefficient. Secondly, the data, after noise reduction, was converted into images, called snowflake images, according to symmetric polar coordinate method. Different fault types of axial piston pump can be identified by observing the snowflake images. After that, in order to evaluate the research results objectively, the obtained images were converted into Gray-Level Cooccurrence Matrixes. Their multiple eigenvalues were extracted, and the eigenvectors consisting of multiple eigenvalues were classified by Fuzzy C-Means clustering algorithm. Finally, according to the accuracy of classification results, the feasibility of applying the symmetric polar coordinate method to axial piston pump fault diagnosis has been validated.
1. Introduction
With their outstanding advantages, such as light weight, great power-mass ratio, flexible control, and fast response speed, hydraulic systems have received extremely high attention and extensive applications in industry, agriculture, and national defence [1]. The hydraulic pump, as the power component of the hydraulic system, converts the mechanical energy, provided by the prime mover, into the pressure energy of the working medium. A working mechanism of hydraulic system can be driven only by pressure energy of the working medium. Thus, hydraulic pump is also called the heart of hydraulic system [2, 3].
Once the hydraulic pump fails, it will affect the entire system. Sometimes, the faults even lead to terrible safety accidents. In addition, many developed countries around the world have put forward the concept of intelligent production; and the fault diagnosis technology of equipment is one of the important contents [4, 5]. Therefore, research on fault diagnosis technology for hydraulic pumps is particularly important for equipment safety and intellectualization [6, 7]. In this paper, a new algorithm based on symmetric polar coordinate method and Fuzzy C-Means (FCM) clustering was proposed for fault diagnosis of axial piston pump. The method can project the time-domain vibration signals into the polar coordinate through the symmetrical polar coordinate. Then the snowflake images were generated according to the mirror symmetry plane rotation angle φ and angle magnification factor . After that, the snowflake images were transformed into Gray-Level Cooccurrence Matrix (GLCM), which is easy to calculate by computer. The eigenvalues of the matrix were extracted. Finally, the FCM algorithm was used to cluster the eigenvalues to achieve the purpose of fault diagnosis. In addition, the vibration signals of the axial piston pump always have noise interference. This paper adopted a method, combination of ensemble empirical mode decomposition (EEMD) and Pearson correlation coefficient, to denoise the signal.
The difference between the method described in this paper and the traditional methods is that the symmetric polar coordinate is employed. This method has the advantages of small computation and symmetrical distribution. Hence, compared with the traditional methods that generate the time-domain waveform and frequency domain waveform, the image generated by the symmetric polar coordinate can reflect the tiny difference of the signals more clearly. Because of the advantages of the symmetric polar coordinate, the diagnosis rate based on it is higher.
Four types of vibration signals of axial piston pumps have been sampled, such as swash plate wear, loose slipper, sliding shoe wear, and normal operation, through experiments. The diagnosis methods introduced in this paper are used to analyse these types of signals. The operation finally gets the FCM clustering result. The analysis results show that the method proposed in this paper has a high accuracy rate for the fault diagnosis of the axial piston pump.
2. Related Theories and Methods
2.1. Noise Reduction Algorithm Combined EEMD with Pearson Correlation Coefficient
Axial plunger pump is a typical rotating machine [8, 9]. When a certain part of it is worn or cracked, it will generate some abnormal vibration signals through the periodic rotation of the pump [10, 11]. However, the installation of vibration sensors cannot alter the space structure of the pump. The vibration sensors hence can only be installed on the shell of pump. In the process of vibration signals transmitting to the sensors, noise will inevitably be mixed into the signals and reduces the signal-noise ratio [12, 13]. Therefore, the sampled original signals are nonstationary and nonlinear signals [14, 15].
Because of the characteristics of the vibration signals, being nonstationary and nonlinear, this paper adopted the method of EEMD combined with Pearson correlation coefficient to denoise the original vibration signals.
2.1.1. Empirical Mode Decomposition
In 1998, American scientist Norden E. Huang proposed a new method, Hilbert-Huang Transform (HHT), for processing nonstationary signals. Empirical mode decomposition (EMD) was also introduced firstly as an important part of this method. The EMD algorithm flow chart is shown in Figure 1.
According to the flow chart shown in Figure 1, the concrete steps of the EMD algorithm are as follows: Step 1: Parameters are initialized, and all local extremes are computed from the signal x(t). Step 2: The upper envelope E1(t) and lower envelope E2(t) of signal x(t) are constructed by cubic splines; then the mean of the envelopes mi(t) is calculated. Step 3: The mean mi(t) is subtracted from the single x(t); then the ith component hi(t) is gained. Step 4: If the component hi(t) is in accordance with the conditions of Intrinsic Mode Function (IMF), hi(t) will be taken as the ith IMF component ci(t), the difference between x(t) and hi(t) is denoted as the residual r(t) and i is added with 1. Otherwise, x(t) will be set as hi(t), and steps 1–3 should be repeated. Step 5: If the residual r(t) is monotonous, the decomposition will be stopped. Otherwise, x(t) will be set as r(t), and steps 1–4 will be repeated.
Finally, the original signal x(t) can be expressed aswhere i = 1,2,3, …, N; N is the number of IMF components obtained by decomposition [16].
Compared with previous signal processing methods, EMD can decompose signals without setting any basis functions. It has advantages of intuitiveness, directness, and posterior and self-adaptation. Based on these advantages, EMD can be used theoretically to decompose various signals, including nonstationary and nonlinear signals. Therefore, EMD was applied to various engineering fields as soon as it was proposed.
However, EMD has end effects and modal aliasing. These problems affect the quality and performance of decomposition. Because of these defects, Wu and Huang proposed EEMD in 2009 [17].
2.1.2. Ensemble Empirical Mode Decomposition
The EEMD is based on EMD, and it can effectively avoid modal aliasing. Its principle of decomposition is the following: an original signal is added with multiple groups of Gaussian white noise with zero mean. Then the EMD algorithm is executed on the processed signal, and the signal will be automatically decomposed into different frequency bands. Because the Gaussian white noise average value is zero, the white noise can be eliminated from the signal through averaging operation and restore to the original signal [18]. The EEMD algorithm flow chart is shown in Figure 2.
According to the flow chart shown in Figure 2, the concrete steps of the EEMD algorithm are as follows: Step 1: The EMD execution times M and the amplitude coefficient of Gaussian white noise a are initialized, respectively, and i is set as 1. Step 2: Gaussian white noise ni(t), with a zero-mean value and a constant standard deviation, is added to the original signal x(t) many times to obtain a new signal xi(t): where ni(t) is the ith time Gaussian white noise sequence added. Step 3: EMD is performed on xi(t) and several IMF components cij(t) and a residual ri(t) are obtained, where cij(t) is the jth IMF obtained by EMD after adding the ith Gaussian white noise to the signal x(t). ri(t) is the residual after the ith EMD. Step 4: If i < M, i will be added with 1 and steps 2 and 3 will be repeated until i = M. Step 5: The average of all IMF components cj(t) and residual r(t) are obtained after M times EMD:where i = 1,2,3, …, M, j = 1,2,3, …, N, and N is the number of IMF components. cj(t) is the ith IMF component of EEMD. r(t) is the residual of EEMD.
The original signal can be reconstructed with multiple cj(t) and a r(t) [19]:
2.1.3. Pearson Correlation Coefficient
The correlation coefficient was first proposed by the statistician Carl Pearson in the 19th century. He established the maximum likelihood method, based on the correlation and regression statistical concepts previously proposed by Galton, Weldon, and others. He used the correlation coefficient r to represent the correlation degree of bivariate normal distribution. It should be noted that the Pearson correlation coefficient is only one type of correlation coefficients. The correlation coefficients described below refer to the Pearson correlation coefficient. The formula is as follows [20, 21]:where Cov(X, Y) is the covariance of X and Y; Var|X| and Var|Y| are the variances of X and Y, respectively. The value range of the correlation coefficient r is from −1 to 1. If r > 0, X and Y will be positively correlated. If r < 0, X and Y will be negatively correlated. If r = 1, X and Y will be identical. If r = 0, X and Y will have zero correlation. If r = −1, X will be equal to minus Y.
2.1.4. Simulation Signal-Noise Reduction
In order to verify the effectiveness of the noise reduction algorithm combined EEMD with Pearson correlation coefficient, a set of simulated signals are constructed and processed with this algorithm. Their sampling frequency is 10 kHz, sampling time is 2 s, and sampling number is 20000. The mathematical expression of the original signal x(t) is as follows:where x1(t) is a Sine signal and x2(t) is an amplitude-modulated signal.
White Gaussian noise of 8 dB was added to the original signal x(t). The original signal x(t), the noise signal, and the synthetic signal are shown in Figure 3.
The synthetic signal decomposed by EEMD is shown in Figure 4(a), and each component is observed in the frequency domain, as shown in Figure 4(b).
(a)
(b)
The Pearson correlation coefficients between each component and the original signal are calculated, and the results are shown in Table 1.
The signal is reconstructed through the two components with the greatest correlation coefficients, as shown in Figure 5.
In Figure 5, the reconstructed signal is basically consistent with the original simulation signal. After calculation, the correlation coefficient between the reduced noise signal and the original signal is 99.55%, which proves that this method can effectively reduce the noise in the noised signal.
2.2. Symmetrical Polar Coordinate Algorithm
The symmetrical polar coordinate algorithm transforms the sampled time-domain signal into polar coordinate and expresses it in the form of image. This image is called a snowflake image. Because of its symmetry, snowflake images can well show the differences between each other. In addition, they are more intuitive than time-domain waveform graphs to show the difference between different fault types.
The basic principle of the symmetrical polar coordinate algorithm is as follows:
The amplitude of signal at time i is x(i) and at time i + l it is x(i + l). It can be converted to a point in polar coordinate by the following formulas:where r(i) is the polar coordinate radius; xmax is the maximum amplitude in the signal; xmin is the minimum amplitude in the signal; α(i) is the rotation angle in the counter-clockwise direction from the mirror symmetry plane; φ is the rotation angle of the mirror symmetry plane; k is the angle magnification factor; β(i) represents the rotation angle in the clockwise direction from the mirror symmetry plane [22]. The physical quantities are represented in polar coordinate as shown in Figure 6.
The size of φ is inversely proportional to the number of mirror planes. If φ is too large, the number of petals in the snowflake image will be too small and the information contained by the graphics will be less. But φ cannot be too small. If it is too small, the number of petals will be too many, even overlapping. It will lead to the graphics being too chaotic to find the characteristics [23]. Usually φ is set as 60°, and the resulting mirror plane angles are 0°, 60°, 120°, 180°, 240°, and 300°. These six mirror planes are evenly distributed in polar coordinate to form a six-petal snowflake image. The value of l is proportional to the width of petals of the snowflake images; generally, 3∼10 is better. The value of k represents the maximum angle that half of the petals can cover. The selection of the value of k will directly affect the degree of overlap between the petals. Generally, 20°∼60° is better [24].
2.3. Gray-Level Cooccurrence Matrix and Its Eigenvalues
The statistical method of Gray-Level Cooccurrence Matrix (GLCM) was proposed by Haralick et al. in the early 1970s. It is a universal image analysis method for images that have texture information in the spatial distribution relationship between pixels. The GLCM is generated as follows: take one pixel (x, y) and another pixel (x + dx, y + dy) in the gray image. dx is the deviation between two pixels in x direction. dy is similar to dx. The gray values of these two pixels are and , respectively. The relationship of pixels is shown in Figure 7.
The calculation formula of the probability is as follows:where i = 0, 1, 2, …, N − 1, j = 0, 1, 2, …, N − 1, and N is the gray level of the image; x, y are the horizontal and vertical values of coordinate of a certain pixel in the image, x = 0, 1, 2, …, Lx − 1, y = 0, 1, 2, …, Ly − 1, Lx, Ly are the number of rows and columns of image pixels, respectively; δ is the distance between two pixels; and θ is the angle between the connecting line of two pixels and the horizontal line [25].
There are N2 combinations of and . Arranging the probability of each combination into a square matrix is the GLCM. The structure of the matrix is shown in Figure 8.
Because the deviation values dx and dy could take different values, the GLCM can be obtained under different position relationships. Generally, the generation direction θ of the GLCM takes four directions (0°, 45°, 90°, and 135°), as shown in Figure 9. Different generated directions reflect the texture features of different directions of the image, and the GLCM obtained from different parameters is also different [26].
Since the GLCM cannot directly reflect the texture of the image, some statistics based on the matrix are usually used as classification features. R Haralick et al. proposed a total of 14 eigenvalues of GLCM. In this paper, four commonly and effectively used statistical features are employed: Angular Second Moment (ASM), Contrast (Con), Correlation (Cor), and Homogeneity (Hom).
2.3.1. Angular Second Moment
The ASM of the GLCM is also called energy. This feature value is the sum of the squares of each matrix element. It reflects the uniformity of the texture distribution of an image and the thickness of the texture. The formula is
2.3.2. Contrast
Con reflects the sharpness of the image and the gray-level difference of the texture. The greater the gray-level difference, the greater edge the contrast value. The formula is
2.3.3. Correlation
Cor reflects the degree of similarity of grayscale images in the row and column directions. The higher the degree of similarity, the greater the autocorrelation. The formulas are
2.3.4. Homogeneity
Hom reflects the degree of local changes in the image texture. When the local area differs greatly, the value is large and vice versa. The formula is
2.4. Fuzzy C-Means Clustering Algorithm
In most clustering problems, the samples in the data set cannot be clearly classified into a certain category. It is even wrong to specify a sample to a specific category. Lotfi Zadeh, an American automatic control expert, is the “Father of Fuzzy Sets.” He proposed fuzzy set theory and fuzzy logic to solve these problems in the clustering process in the middle of the 20th century. These theories were quickly applied in many fields and the emergence of Fuzzy C-means (FCM) clustering algorithm also credits to this. The FCM algorithm does not give a certain limit to the category, also called cluster, like the K-means clustering algorithm. But there is a weight for each sample and category, which shows the membership of the sample to the category. The sum of the memberships of all samples to all categories is 1. Compared with the weights given by statistical methods, this method can better avoid the difficulty of selecting statistical models and give a natural and nonprobabilistic classification result [27, 28].
The core of FCM is minimizing the objective function Jm, the sum of squares of errors. The formula is as follow:where N is the total number of samples; C is the number of categories; m is the weighted index number; uij is the membership of sample xi to category j; and cj represents the center of category j.
The flow of the FCM is continuously iterating the membership degree uij and the category center cj to make the objective function reach the best. The calculation formulas for these two values are as follows [29, 30]:
3. Experimental System and Fault Data Sampling
In order to obtain the data required for this paper, our research team has built a hydraulic pump failure simulation test bench. The hydraulic schematic diagram is shown in Figure 7.
In Figure 10, the tested plunger pump 10 was connected to the motor. There was also a vane pump 3 to provide sufficient hydraulic oil to the tested plunger pump. A direct-acting relief valve 23 was used to ensure the stable pressure of the oil source. The pilot-operated proportional relief valve 21 and pilot-operated relief valve 22 were switched by the two-position three-way electromagnetic directional valve 19. It could establish the different working pressures of the tested plunger pump. There were three vibration sensors 11, which were, respectively, fixed on the plunger pump housing along the radial horizontal direction x, radial vertical direction y, and axial direction z. Most importantly, this test device could sample the front’s and the rear’s pressure information of the tested plunger pump and the motor speed information at the same time. The basic parameters of some main components are shown in Table 2. The picture of the test bench is shown in Figure 11.
The software of data acquisition used was NI’s LabVIEW, and the acquisition card was a USB6221 data acquisition card. The acquisition system can guarantee a 20 kHz acquisition rate per channel. The front panel of the acquisition program is shown in Figure 12.
Before the experiment, the corresponding faulty parts had been prepared. They included the grinding swash plate and the sliding shoe and pulling the plunger and the sliding shoe artificially. In the experiment, the normal equipment operated under a pressure of 10 MPa, the computer sampled vibration signals in the x, y, and z directions, and the sampling frequency was 20 kHz. After that, the previously prepared faulty parts were replaced one by one into the normal plunger pump. Then, the vibration signals under the faults were sampled under the same conditions.
The three-direction vibration signals of the axial piston pump sampled in the normal state are shown in Figure 13.
The amplitude of the vibration signals in the z-axis direction is significantly higher than that in the other two directions as shown in Figure 13. It means that the z-axis direction is the most sensitive to the vibration of the piston pump. For this reason, this article used the z-axis direction vibration signals to diagnose the faults of the axial piston pump.
4. Data Processing and Diagnostic Analysis
4.1. Noise Reduction
After calculation, each circle corresponds to 800 data points, so the data length is set as 1000 to include a complete turn. In this paper, MATLAB was used to perform all calculations. Before the noise reducing Fast Fourier Transform (FFT) algorithm was used on the various vibration signals sampled from the axial piston pump and frequency spectrum diagrams were generated. These diagrams reflected the distribution of each signal in the frequency domain, as shown in Figure 14.
(a)
(b)
(c)
(d)
In this paper, the swash plate wear vibration signals were taken as an example to explain the process of noise reduction by the method combined EEMD with the Pearson correlation coefficient. The noise reduction processes for other conditions are similar.
First, the sampled vibration signals were decomposed by the EEMD algorithm and obtained each IMF component signal, as shown in Figure 15.
(a)
(b)
Second, the FFT was applied to each component to observe the distribution of each component in the frequency domain, as shown in Figure 16.
(a)
(b)
According to the observation of Figure 16, the first and second components distributed widely in the frequency domain. They were considered to be noise components.
Third, the Pearson correlation coefficients between the original signal and the components are calculated, except the first and second IMF components. The results are shown in Table 3.
Finally, according to the correlation coefficients between each IMF component and the original signal, the five components with the largest correlation coefficients were selected to reconstruct the signal and reduce noise. The frequency spectrum diagram of swash plate wear signal after noise reduction is shown in Figure 17.
Compared with the original vibration signal, the noise in the high frequency part of the signal has been better eliminated after noise reduction, as shown in Figure 17. It could be considered that the noise reduction effect was obvious.
Considering the change of signals themselves in the process of signal sampled, the distribution of each group of signals on the IMF components might not be completely consistent. For this reason, certain fixed IMF components were not suitable for all signals’ reconstruction. Therefore, in the noise reduction process of each signal, the Pearson correlation coefficients between each IMF component and the original signal need to be calculated; and the five components, except the first and second, with the largest correlation coefficients were selected to reconstruct the signals. This could reduce the loss of useful information of each signal.
4.2. The Snowflake Images Generated
80 sets of vibration signals were obtained after noise reduction, and each state has 20 sets equally. These signals were substituted into the polar coordinate algorithm and 80 snowflake images were obtained. In this paper, φ was set as 60°, and the resulting mirror plane angles were 0°, 60°, 120°, 180°, 240°, and 300°; l was 4; and k was 30°. Figure 18 shows the snowflake images corresponding to various conditions.
(a)
(b)
(c)
(d)
It is intuitive to show the difference between the snowflake images of normal state and various fault states, as shown in Figure 18. Among them, (1) the snowflake image’s petals of the normal state were in the shape of a long water drop. They were thick and not curved and were evenly distributed on the entire circumference. (2) The image’s petals of the swash plate wear like short and thick water drop. Each two petals between the two mirror planes were closer, but on both sides of a single mirror plane they were relatively distant. The centroid of petals was farther from the center of the circle. (3) The image’s petals of the loose slipper were curved water drop. Each two petals on both sides of a single mirror plane were close first and then separated as the distance increased in the radial direction. A single petal was slender near the center, thicker away from the center, and with more divergent points at the end. (4) The image’s petals of the sliding shoe wear were crescent-shaped. The petals were arc-shaped near the mirror plane and flat away from the mirror plane. A single petal was thin at both ends and thick in the middle, with fewer divergent points.
4.3. Gray-Level Cooccurrence Matrix Generated and Feature Extraction
In normal circumstances, the gray level of a gray image is generally 256 levels, from 0 to 255. However, in the calculation process, 256 levels will produce a tremendous amount of computation. For example, the snowflakes in 4.2 have 469 × 469 pixels. If a computer uses 256 gray levels and operates this picture, it is going to calculate 1.4 × 1010 times. On the premise of ensuring the image texture as much as possible, the number of operations can be reduced by reducing the gray level. Usually, the gray level is compressed to 16 or 8 to reduce the size of the GLCM. In this paper, the gray scale is set to 16 levels.
Using the GLCM algorithm, the snowflake images generated in Section 4.2 were converted into matrixes in the four directions of 0°, 45°, 90°, and 135°. There were 80 matrixes generated in each state, and a total of 320 matrixes were gotten. Then, the respective eigenvalues like ASM, Con, Cor, and Hom were calculated for each matrix. Finally, the four eigenvalues of each state were averaged as the eigenvalue benchmark of the GLCM of this state. The average results are shown in Table 4.
4.4. Fuzzy C-Means Clustering Results
In this paper, 80 samples were used to train algorithm. These 80 samples were 20 normal samples, 20 swash plate wear samples, 20 loose slipper samples, and 20 sliding shoe wear samples. Similarly, there were 80 samples as a test set. They were in identical condition but were different in data. The numbers of samples and eigenvalues contained in the training set and test set are shown in Table 5.
These eigenvalues of test samples were put into the FCM algorithm for calculation, and the classification results are shown in Figure 19. When drawing the classification results, the dimensionality was reduced to realize its drawing in the three-dimensional space [31].
After being clustered by the FCM algorithm, the four states of the axial piston pump, the normal state and the three types of failure states, were clearly separated, which is shown in Figure 19. The accuracy rate of the classification results is shown in Table 6.
In addition, this paper has used EMD to decompose the same data. The three IMF components with the greatest correlation with the original signal were selected. The energy feature values of these components were extracted. Then the FCM clustering algorithm was used to classify these feature values. The clustering result is shown in Figure 20.
Compared with the classification results shown in Figure 19, the classification of the four types of states in Figure 20 is not clear. The classification accuracy rate is also lower, as shown in Table 7. Thus, the effectiveness and superiority of the method proposed in this paper are proved.
5. Conclusions
In this paper, a fault diagnosis algorithm for axial piston pump was proposed which was primarily based on symmetrical polar coordinate image and FCM. First, the noise reduction algorithm, combined EEMD with Pearson correlation coefficients, was employed to preprocess the original signals. Second, the symmetrical polar coordinate algorithm converted the processed samples into snowflake images. Third, the snowflake images were transformed to GLCM. Fourth, the eigenvalues corresponding to the sample could be gotten by the eigenvalue of GLCM algorithm, and the samples of eigenvalues were gained. Finally, the FCM algorithm performed clustering on samples according to the clustering center and completed the fault diagnosis. Through the above content, the conclusions can be drawn:(1)Compared with time-domain waveform diagram, the state’s snowflake images, drawn by the symmetrical polar coordinate algorithm, could reflect the difference of the fault and normal types of the axial piston pump more intuitively. Because people’s eyes are more sensitive to symmetrical patterns, the type of fault can be identified directly by observing the snowflake images.(2)Due to the images having a high degree of differentiation for each state, these GLCM eigenvalues samples, derived from the image, represent the characteristics of the failure state more effectively, and the result of FCM clustering was exact. After statistics, the comprehensive accuracy rate of the fault diagnosis algorithm proposed in this paper was 98.75%.(3)Compared with the EMD method of extracting energy eigenvalues, the accuracy of the method proposed in this paper was significantly higher, and the accuracy rate of the EMD method was only 92.5%. The superiority of the method proposed in this paper has been proved.
The method proposed in this paper is feasible and is more effective than other departed methods in fault diagnosis. This is a research topic worthy of further study. Moreover, this fault diagnosis algorithm can be extended to other rotating machinery. We will continue to study this area in the future.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grants nos. 51875498 and 51475405) and Key Project of Natural Science Foundation of Hebei Province, China (Grants nos. E2018203339 and F2020203058).