Axial Piston Pump Fault Diagnosis Method Based on Symmetrical Polar Coordinate Image and Fuzzy C-Means Clustering Algorithm

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 ﬁrstly by the combination of ensemble empirical mode decomposition and Pearson correlation coeﬃcient. Secondly, the data, after noise reduction, was converted into images, called snowﬂake images, according to symmetric polar coordinate method. Diﬀerent fault types of axial piston pump can be identiﬁed by observing the snowﬂake 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 classiﬁed by Fuzzy C-Means clustering algorithm. Finally, according to the accuracy of classiﬁcation results, the feasibility of applying the symmetric polar coordinate method to axial piston pump fault diagnosis has been validated.


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]. e 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. us, 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]. erefore, 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. e method can project the time-domain vibration signals into the polar coordinate through the symmetrical polar coordinate. en the snowflake images were generated according to the mirror symmetry plane rotation angle φ and angle magnification factor k. After that, the snowflake images were transformed into Gray-Level Cooccurrence Matrix (GLCM), which is easy to calculate by computer.
e 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.
is paper adopted a method, combination of ensemble empirical mode decomposition (EEMD) and Pearson correlation coefficient, to denoise the signal. e difference between the method described in this paper and the traditional methods is that the symmetric polar coordinate is employed. is 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. e diagnosis methods introduced in this paper are used to analyse these types of signals. e operation finally gets the FCM clustering result. e analysis results show that the method proposed in this paper has a high accuracy rate for the fault diagnosis of the axial piston pump.

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. e 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 signalnoise ratio [12,13]. erefore, 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.

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. e 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: e upper envelope E 1 (t) and lower envelope E 2 (t) of signal x(t) are constructed by cubic splines; then the mean of the envelopes m i (t) is calculated.
Step 3: e mean m i (t) is subtracted from the single Step 4: If the component h i (t) is in accordance with the conditions of Intrinsic Mode Function (IMF), h i (t) will be taken as the ith IMF component c i (t), the difference between x(t) and h i (t) is denoted as the residual r(t) and i is added with 1. Otherwise, x(t) will be set as h i (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 as where 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, Input signal x (t), i = 1 Fit the upper envelope E 1 (t) and lower envelope E 2 (t) and find the mean m i (t) The ith component h i (t) is the difference between x (t) and m i (t) and c i (t) is the residual amount, denoted as r (t) x (t) = r (t)  erefore, EMD was applied to various engineering fields as soon as it was proposed.
However, EMD has end effects and modal aliasing. ese problems affect the quality and performance of decomposition. Because of these defects, Wu and Huang proposed EEMD in 2009 [17].

Ensemble Empirical Mode Decomposition.
e 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. en 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]. e 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: e 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 n i (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 x i (t): where n i (t) is the ith time Gaussian white noise sequence added.
Step 3: EMD is performed on x i (t) and several IMF components c ij (t) and a residual r i (t) are obtained, where c ij (t) is the jth IMF obtained by EMD after adding the ith Gaussian white noise to the signal x(t). r i (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: e average of all IMF components c j (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. c j (t) is the ith IMF component of EEMD. r(t) is the residual of EEMD. e original signal can be reconstructed with multiple c j (t) and a r(t) [19]:

Pearson Correlation Coefficient.
e 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. e correlation coefficients described below refer to the Pearson correlation coefficient. e formula is as follows [20,21]: Shock and Vibration where Cov(X, Y) is the covariance of X and Y; Var|X| and Var|Y| are the variances of X and Y, respectively. e 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.

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. eir sampling frequency is 10 kHz, sampling time is 2 s, and sampling number is 20000. e mathematical expression of the original signal x(t) is as follows: where x 1 (t) is a Sine signal and x 2 (t) is an amplitudemodulated signal. White Gaussian noise of 8 dB was added to the original signal x(t). e original signal x(t), the noise signal, and the synthetic signal are shown in Figure 3. e 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 e Pearson correlation coefficients between each component and the original signal are calculated, and the results are shown in Table 1.
e 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.

Symmetrical Polar Coordinate Algorithm.
e symmetrical polar coordinate algorithm transforms the sampled time-domain signal into polar coordinate and expresses it in the form of image. is 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. e basic principle of the symmetrical polar coordinate algorithm is as follows: e 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 P(r(i), α(i), β(i)) in polar coordinate by the following formulas: where r(i) is the polar coordinate radius; x max is the maximum amplitude in the signal; x min is the minimum amplitude in the signal; α(i) is the rotation angle in the counterclockwise 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]. e physical quantities are represented in polar coordinate as shown in Figure 6. e 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°. ese six mirror planes are evenly distributed in polar coordinate to form a six-petal snowflake image. e value of l is proportional to the width of petals of the snowflake images; generally, 3∼10 is better. e value of k represents the maximum angle that half of the petals can cover. e selection of the value of k will directly affect the degree of overlap between the petals. Generally, 20°∼60°is better [24].

Gray-Level Cooccurrence Matrix and Its Eigenvalues.
e 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. e GLCM is generated as follows: take one pixel (x, y) and another pixel (x + d x , y + d y ) in the gray image. d x is the deviation between two pixels in x direction. d y is similar to d x . e gray values of these two pixels are g i and g j , respectively. e relationship of pixels is shown in Figure 7. e calculation formula of the probability P(g i , g j , δ, θ) is as follows:    Shock and Vibration 5 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, . . ., L x − 1 , y � 0, 1, 2, . . ., L y − 1 , L x , L y 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].
ere are N 2 combinations of g i and g j . Arranging the probability of each combination into a square matrix is the GLCM. e structure of the matrix is shown in Figure 8.
Because the deviation values d x and d y 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).

Angular Second Moment.
e ASM of the GLCM is also called energy.
is 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. e formula is     Shock and Vibration

Correlation.
Cor reflects the degree of similarity of grayscale images in the row and column directions. e higher the degree of similarity, the greater the autocorrelation. e formulas are

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. e formula is

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. ese theories were quickly applied in many fields and the emergence of Fuzzy C-means (FCM) clustering algorithm also credits to this. e 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. e 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]. e core of FCM is minimizing the objective function J m , the sum of squares of errors. e formula is as follow: where N is the total number of samples; C is the number of categories; m is the weighted index number; u ij is the membership of sample x i to category j; and c j represents the center of category j. e flow of the FCM is continuously iterating the membership degree u ij and the category center c j to make the objective function reach the best. e calculation formulas for these two values are as follows [29,30]:

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. e hydraulic schematic diagram is shown in Figure 7.
In Figure 10, the tested plunger pump 10 was connected to the motor.
ere was also a vane pump 3 to provide sufficient hydraulic oil to the tested plunger pump. A directacting relief valve 23 was used to ensure the stable pressure of the oil source. e pilot-operated proportional relief valve 21 and pilot-operated relief valve 22 were switched by the twoposition three-way electromagnetic directional valve 19. It could establish the different working pressures of the tested plunger pump. ere 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. e basic parameters of some main components are shown in Table 2. e picture of the test bench is shown in Figure 11. Shock and Vibration e software of data acquisition used was NI's Lab-VIEW, and the acquisition card was a USB6221 data acquisition card.
e acquisition system can guarantee a 20 kHz acquisition rate per channel. e front panel of the acquisition program is shown in Figure 12.
Before the experiment, the corresponding faulty parts had been prepared. ey 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. en, the vibration signals under the faults were sampled under the same conditions.    8 Shock and Vibration e three-direction vibration signals of the axial piston pump sampled in the normal state are shown in Figure 13. e 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.

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. ese diagrams reflected the distribution of each signal in the frequency domain, as shown in Figure 14.
In this paper, the swash plate wear vibration signals were taken as an example to explain the process of noise reduction  First, the sampled vibration signals were decomposed by the EEMD algorithm and obtained each IMF component signal, as shown in Figure 15.  Second, the FFT was applied to each component to observe the distribution of each component in the frequency domain, as shown in Figure 16.
According to the observation of Figure 16, the first and second components distributed widely in the frequency domain. ey were considered to be noise components.
ird, the Pearson correlation coefficients between the original signal and the components are calculated, except the first and second IMF components. e 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. e 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. erefore, 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. is could reduce the loss of useful information of each signal.

e Snowflake Images
Generated. 80 sets of vibration signals were obtained after noise reduction, and each state has 20 sets equally. ese 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.   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. ey were thick and not curved and were evenly distributed on the entire circumference. (2) e 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. e centroid of petals was farther from the center of the circle. (3) e 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) e image's petals of the sliding shoe wear were crescent-shaped. e 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.

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 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°.
ere were 80 matrixes generated in each state, and a total of 320 matrixes were gotten. en, 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. e average results are shown in Table 4.

Fuzzy C-Means Clustering
Results. In this paper, 80 samples were used to train algorithm. ese 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. ey were in identical condition but were different in data. e numbers of samples and eigenvalues contained in the training set and test set are shown in Table 5.
ese 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. e accuracy rate of the classification results is shown in Table 6.   Figure 19: Clustering results of symmetrical polar coordinate combined with GLCM.  In addition, this paper has used EMD to decompose the same data. e three IMF components with the greatest correlation with the original signal were selected. e energy feature values of these components were extracted. en the FCM clustering algorithm was used to classify these feature values. e clustering result is shown in Figure 20.  Table 7.
us, the effectiveness and superiority of the method proposed in this paper are proved.

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.
ird, 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. rough 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%. e superiority of the method proposed in this paper has been proved.
e method proposed in this paper is feasible and is more effective than other departed methods in fault diagnosis.
is 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
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.