Abstract
To overcome the shortcomings that the early fault characteristics of rolling bearing are not easy to be extracted and the identification accuracy is not high enough, a novel collaborative diagnosis method is presented combined with VMD and LSSVM for incipient faults of rolling bearing. First, the basic concept of VMD was introduced in detail, and then, the adaptive selection principle of parameter K in VMD was constructed by instantaneous frequency mean. Furthermore, we used Lagrangian polynomial and Euclidean norm to verify the value of K accurately. Secondly, we proposed a classification algorithm based on PSO-optimized LSSVM. Meanwhile, the flowchart of the classification algorithm of fault modes may be also designed. Third, the experiment shows that the presented algorithm in this paper is effective by using the existing failure data provided by the laboratory of Guangdong Petrochemical Research Institute. Finally, some conclusions and application prospects were discussed.
1. Introduction
In recent years, the machinery has become more high-speed, intelligentized, and complicated with the development of the modern industrialization. As we all know, the rotating machinery is the cornerstone of transportation, power electronics, and manufacturing. So, how to guarantee the security of whole rotating machinery systems is very important in the industrial field. In the actual industrial scenario, the engineers and researchers have noticed that the safety of the bearings is often one of the critical joints, which ensures the global safety of whole rotating machinery [1]. Therefore, it is essential to detect and assess the performance of the running state of the bearings. The traditional fault diagnosis methods that judge and evaluate the running state of the bearing are operated or implemented by observing the frequency of the vibration signal. The skeleton of these methods consists of just three steps: signal processing, feature extraction, and fault pattern recognition. In most realistic scenarios, signal processing is often used as the preparing work for the feature extraction. Of course, the feature extraction is also used as the prepared work for the fault pattern recognition because the classification accuracy of the fault modes is the final objective in fault diagnosis of the bearing, so the signal processing and feature extraction are often integrated to analyze the vibration signals of the bearing. And, how to implement them becomes critical.
For the question raised above, the scholars have presented and constructed some models such as Empirical Mode Decomposition (EMD), Wavelet Transform (WT), Local Mean Value Decomposition (LMVD), and Variational Mode Decomposition (VMD) in references [2–4]. The experiment results show that these methods may acquire the most of the valuable information in the specified scenarios. Unfortunately, almost all these methods have some shortcomings. For example, the EMD and LMD have the phenomenon such as modal aliasing and endpoint effect. The WT needs to select the wavelet base and decomposition scale because the finite length may cause inaccurate decomposition of complete components. In the VMD, if the parameter K is wrongly selected, the phenomenon such as overdecomposition or underdecomposition will appear. To overcome these shortcomings, some improved algorithms have been presented such as Simplistic Geometric Mode Decomposition (SGMD), Adaptive Chirped Mode Decomposition (ACMD), and New Spectral Analysis Methods (NSAM) in [5–7]. Especially, to address the shortcoming of the EMD, the study in [8] has constructed the EMD envelope correction method using B-spline interpolation and base spline. These methods may alleviate the modal aliasing problem of high-frequency signals. Meanwhile, to optimally select the parameter K of VMD, the genetic variation sample group, kurtosis criterion variational mode decomposition, and self-organizing mapping (SOM) neural network have been adopted to adaptively determine the optimal value of the parameter K in [9–12]. To verify the effectiveness of these new methods, some experiment examples have been used to simulate in [13–17]. The simulated results showed that these improved models may solve the shortcomings to a certain extent. For practical application, the constant improvement of the existing methods is the goal of the engineers and scholars. Thus, we will treat the problem in this paper.
On the contrary, in the view of fault diagnosis, to get the accurate classification of fault modes is the other main objective of the bearing fault diagnosis. In fact, an excellent pattern recognition method of the fault modes has an important influence for the final diagnosis accuracy. Based on this objective, support vector machine (SVM), least-squares’ support vector machine (LSSVM), BP neural network (BPNN), fuzzy logic (FM), and other methods have been successfully applied in [18–23]. And, then, these fault pattern recognition methods have been widely used in different industrial environments. Further, some improved methods of the fault pattern recognition were studied in [24, 25]. For example, the double support-vector machine and smooth iterative online-support tensor algorithm are proposed to improve the performance of the traditional support vector machine in [26, 27]. The least-squares’ ground projection method of the double support-vector machine reduces the diagnostic error in [28]. Meanwhile, to optimize the penalty factor C and kernel parameter of LSSVM, some new algorithms such as the Moth-flame Optimization (MFO), the von Neumann Topology Whale Optimization Algorithm (VNWOA), Quantum Particle Swarm (QPS), and Chaotic Antlion Algorithm (CAA) were introduced to implement the optimization operation for enhancing the precision of fault diagnosis in [29–34]. The experiments have verified the performance of these presented algorithms. However, the global searching ability of these algorithms is weak in the real industrial environment. So, searching the improved pattern recognition method to enhance the global searching ability and improve the classification accuracy of fault modes is another concern in our paper.
Based on the two points mentioned above, an improved fault diagnosis of the bearing will be presented combined with the VMD algorithm based on instantaneous frequency optimization and particle swarm optimization least-squares’ support vector machine in our paper. The rest of this paper is arranged as follows. In Section 2, the adaptive selection principle of K value in the VMD algorithm is given in detail. In Section 3, the least-squares’ support vector machine classification model for particle swarm optimization (PSO) is established, and the concrete flowchart of the fault diagnosis process is designed and analyzed. In Section 4, some simulated examples were used to verify the effectiveness of our algorithm through the existing failure data provided by the laboratory of Guangdong Petrochemical Research Institute. Finally, some conclusions are summarized in Section 5.
2. Adaptive Selection Principle of Parameter K in the VMD Algorithm
2.1. The Basic Concept of the VMD Decomposition Principle
The intrinsic mode function (IMF) is defined as an FM and AM signal by VMD decomposition and is expressed as follows:where expresses the instantaneous amplitude, is the instantaneous frequency, and K represents the number of signal components after decomposition.
Suppose the original signal is a multicomponent signal, which is composed of the K IMF component with limited bandwidth, and the central frequency of each IMF is . To determine the bandwidth of each mode, the following steps are used to obtain it:(1)Analytic signals of modal functions are obtained, and Hilbert transformation is performed for each modal function :(2)Mix the estimated center frequency of each modal analytic signal. The spectrum of each modal is modulated to the corresponding baseband as follows:(3)Calculate the square norm of the gradient of the above demodulation signal, and estimate the bandwidth of each modal component. The constraint variational model is established as follows:where represents the K IMF components obtained by decomposition and represents the center frequency of each component.
In order to solve the above constraint variational model, the quadratic penalty factor α and Lagrangian multiplication operator are introduced, where the quadratic penalty factor can guarantee the reconstruction accuracy of the signal in the presence of Gaussian noise and keeps the constraint conditions strict. The expanded Lagrangian expression is as follows:
The multipliers’ alternating direction algorithm is used to update the IMF and its center frequency, and the saddle point of formula (4) is the optimal solution of the original problem. All IMF in the frequency domain can be obtained by the following formula:where is the current residual quantity and is the result of Wiener filtering. The new IMF power-spectrum centers in the algorithm are as follows:where is the power spectrum center.
The above process is the adaptive decomposition process of VMD. From the decomposition principle, it can be known that VMD can well avoid the endpoint effect and modal confusion. But, from the perspective of the actual decomposition process, the VMD algorithm loses the ability to decompose signals independently, which needs to preset the value of K. And, the reasonableness of the K value determines the signal decomposition accuracy of VMD. If the K value is estimated according to the existing observation method, that is, observing the center frequency differentiation of the signal component, the better the center frequency differentiation is, the better the selection of the K value is, and there is no overdecomposition and underdecomposition. However, there is a large error in this method, which makes it difficult to guarantee the decomposition accuracy of the signal and also affects the final classification accuracy. Therefore, this paper proposes a method to optimize the K value of VMD by using instantaneous frequency, which can make use of the difference of instantaneous frequency between signal components to measure the advantage of the K value.
2.2. K Value Estimation of the VMD Component Based on Instantaneous Frequency
If the K value is set too high, the decomposition number will be too large, and then, the component will be fragmenting, especially at high frequency, and the average instantaneous frequency will decrease. If the K value is set too low, the signal will not be completely decomposed, and the superiority of the signal component cannot be reflected. Therefore, the original signal may be decomposed by VMD once the K value was traversed from 2 to 10. And then, the mean values of instantaneous frequencies may be calculated under different K values, and the line graph may be also drawn. Lagrange polynomials were used to fit the discrete points, and the polynomial coefficients under different K values were extracted to construct the coefficient vector, and then, the Euclidean norm of the vector was calculated. The smaller the norm was, the smoother the fitting instantaneous frequency curve was and the better the value was.
The definition of instantaneous frequency is as follows:where is a single-valued function of time t, that is, a single-component signal on frequency. The analytic signal of instantaneous frequency is defined as follows:where is the Hilbert transform of x(t), z(t) is the analytic signal of , is the module of the signal , and is the phase of the signal, which is expressed as . The instantaneous frequency multiplying the integral of the density function over the entire time axis is the average frequency of the signal. Through the Fourier transform of the analytic signal in formula (9), we can get the following formula:
According to the principle of the stationary phase, the integral of equation (10) has a maximum value at the frequency , which needs to meet the condition , namely, . This conclusion indicates that the energy of nonstationary signals is mainly concentrated at the instantaneous frequency. This conclusion indicates that the instantaneous frequency plays a very important role in the recognition, detection, estimation, and modeling of signals, and it can also be used as the evaluation index of VMD decomposition signals. Therefore, on the basis of the original VMD, the original signal is decomposed into different signal components, and the decomposed number is the K value. Then, the average instantaneous frequency under different mode numbers of K from 2 to 10 is calculated to judge the trend of the corresponding line graph. The flatter the trend is, the better the corresponding K value is, so as to realize the optimization of the parameters of the VMD algorithm. However, this method may be misjudged to some extent. In order to measure the instantaneous frequency change more accurately and select the optimal K value, corresponding methods should be adopted to obtain numerical results.
2.3. The Superiority Distinction of K Values Based on Lagrange Polynomials
After calculating the average instantaneous frequency of each component, we need to adopt an index to measure the variation trend of the instantaneous average frequency, which can avoid the error caused by subjective judgment. By fitting the mean instantaneous frequency, the Lagrangian polynomials can be calculated, the vector norm of their coefficients can be compared, and the merits and disadvantages of the K value can be evaluated.
For any point in the interpolation node , make a polynomial of degree n, which satisfies the following formula:
The basic function of Lagrange interpolation is as , and the node is presented as . So, is a polynomial with n null points. Therefore,where is the n-order basic interpolation polynomial or n-order Lagrangian interpolation basis function on n + 1 interpolation nodes. Using the n-order basic interpolation polynomial, the n-order Lagrange polynomial satisfying the interpolation condition can be written as follows:
The average instantaneous frequency of different components is taken as the discrete point of calculating Lagrangian polynomials. After obtaining the simplest form of the Lagrangian polynomial by calculation, the coefficients of the polynomial are extracted and constructed into a vector, and the Euclidean distance of the vector with different K values is calculated. For the coefficient vector , the Euclidean distance of the vector is .
3. Classification Algorithm-Based Least-Squares’ Support Vector Machine with Particle Swarm Optimization
3.1. Basic Concept of LSSVM
The LSSVM is an improved algorithm of the support vector machine; however, as a binary classifier, its core idea remains unchanged, that is, to find a hyperplane that optimizes classification and maximizes the gap between classifications, so as to improve the credibility of classification. The difference between LSSVM and SVM is that the construction of the objective function of LSSVM is through the binomials for the error factor, and at the same time, constraints are equally constraint, and in terms of solving the optimization problem, because the LSSVM is the constraint equation form, the solution is the system of linear equations; to a certain extent, it reduced the difficulty of the algorithm and raised the solving speed, and these advantages make it different from other improvement on the SVM algorithm. The basic principle of this method is described as follows.
The sample of training data can be expressed as , is the input vector of the ith sample, is the target value of the ith sample, and is the number of training samples. In special space, the LSSVM model can be expressed aswhere is the mapping function of nonlinear transformation, which maps the input sample data to the high-dimensional feature space. W is the weight vector, and B is the offset. The objective function of least-squares’ support vector machines is described as
Type is the error variable, and ϒ > 0 is the penalty coefficient. For the simplicity of analyzing, the Lagrangian function is designed as follows:where is the Lagrange multiplier. In the real operation, the KKT optimal condition is used to calculate , , , and . So, the following system of linear equations should be obtained:
In equation (17), I is the identity matrix. According to the Mercer condition, the kernel function can be written as
After and b can be obtained from equations (18) and (19), the nonlinear function of LSSVM can be obtained as follows:
3.2. PSO Parameter Optimization for LSSVM
LSSVM requires two parameters to be tuned: gam and sig2, where gam is the regularization parameter, which determines the minimization and smoothness of the adaptation error, and sig2 is the parameter of the RBF function. PSO optimizes two parameters of LSSVM, gam and sig2, to find the optimal combination of parameters, so as to improve the classification accuracy. The general optimization steps are as follows:(1)Initializing the various parameters of the PSO algorithm, such as population size, learning factor, the maximum number of iterations, initial position, and the velocity of particles.(2)Respectively, in the LSSVM predictive learning sample of each particle vector, get the prediction error of the current position value of the particle, which is used as the fitness value for each particle. Then, the current fitness value of each particle is compared with the best fitness value of the particle itself. If there are many, the current position of the particle is taken as the best position of the particle.(3)The adaptive value of the optimal position of each particle was compared with the adaptive value of the optimal position of the population. If it is better, the optimal position of the particle is regarded as the optimal position of the population.(4)Use formulas (21) and (22) to update the particle velocity and position: where V is the particle speed, X is the position of the current particle, Rand () is a random number between (0, 1), and C1 and C2 are the learning factors, usually C1=C2 = 1.5.(5)Check the result of optimization (maximum number of iterations or expected accuracy) is met or not. If so, the optimization is completed and the optimal solution is found. Otherwise, go to Step (2) and continue the search.
3.3. Rolling Bearing Fault Diagnosis Steps
The acceleration sensor is used to collect four state signals of the rolling bearing, which are normal, bearing external crack, bearing internal crack, and bearing wear. 10 groups of data of each state signal are collected. Take the normal state as an example, they were normal 1, normal 2, …, normal 10. Set the signal period to 1024, which means that, in a file such as “normal 1” with 10240 pieces of data, it is divided into 1024∗10 groups. Based on the above analysis, in order to ensure the realization of fault diagnosis and classification, the classification algorithm flow chart can be designed as follows:(1)Traverse the K value of VMD (K value is the number of original signals decomposed into different components), input the first group of data of each state of the original signal collected into the VMD algorithm, and get K components under different K values (K = 2, 3, ..., 10).(2)Calculate the average instantaneous frequency corresponding to different components of K, draw a line chart, and estimate the K value through the trend of the line chart.(3)In order to further verify the pros and cons of K, the different components of the average instantaneous frequency may be used as computing Lagrange polynomial of discrete points, and then the most simplified forms of Lagrange polynomial can be computed. So, the extraction of polynomial coefficients may be constructed as a vector. In fact, the vector may be demonstrated and calculated under different K values, coefficient of Euclidean distance, and judge norm size to determine the optimal values of K.(4)Set the optimal K value as the mode number that VMD needs to decompose. Decompose the 10 groups of data of each state and extract the time-domain features to form the feature set.(5)The parameters of gamand sig2 of the LSSVM algorithm were optimized by the PSO algorithm.(6)The obtained data set is input into the LSSVM classification algorithm, which is divided into training data and test data. The parameters of the model are updated with the training data, and the test data is input into the trained model to obtain the diagnosis results of fault pattern recognition;
The corresponding flowchart is shown in Figure 1.
4. Classification Experiment
To test the validity and rationality of the algorithm, some test data of the bearing provided by Guangdong Key Laboratory of Petrochemical Equipment Fault Diagnosis was used to simulate and experiment. This data set included the acceleration changes with four different states: normal, bearing internal crack, bearing external crack, and bearing wear. The bearing damage and data acquisition platform are shown in Figure 2.
(a)
(b)
(c)
Figures 2(a)–2(c), respectively, represent bearing internal crack, bearing external crack, and bearing fault data acquisition platform. The acceleration sensor was used for data collection, with the collection period T = 1024. The collected fault data was divided into 10 groups according to the period, and the four different bearing states were divided into 40 groups.
Since these data are the most original vibration signal data, it is difficult to extract subsequent features without processing, so VMD is used to preprocess vibration signals. In order to select the optimal decomposed mode number K, first select 1 group from the 10 groups of data of each bearing state to input VMD, traverse K values from 2 to 10, calculate the instantaneous frequency mean, and get the corresponding broken line chart, as shown in Figure 3.
(a)
(b)
(c)
(d)
From Figure 3, a rough estimate of K may be obtained. Noticing that four kinds of condition is the most gentle the most ideal when K = 2, there is no high frequency under the intermittent and suddenly curved because the original signal only is decomposed into two components. Because the result do not conform to the actual, K = 2 is not as the objects of choice. Thus, it can be estimated that the optimal value K in the normal state is 3, the optimal value K in the bearing wear state is 6, the optimal value K in the bearing internal crack state is 4, and the optimal value K in the bearing external crack state is 5. However, such estimation may lead to wrong choices when the difference between the broken lines is not large. Therefore, it is necessary to choose an index to accurately judge the advantages and disadvantages of the K value. In this paper, Lagrange polynomials are proposed to be established, and instantaneous frequencies under each K value are used as discrete points to calculate Lagrange polynomials. After obtaining the simplest polynomials, coefficients are extracted and corresponding coefficient vectors are calculated. The smaller the Euclidean norm is, the better the K value is. The normal state data are selected here for experimental verification.
Table 1 shows the average instantaneous frequency corresponding to different K values in the normal state of the bearing, and Table 2 shows the corresponding Euclidian norm. From Table 2, it can be seen that the norm is the smallest when K = 2, but because it is not consistent with the actual situation, the value of K is excluded as 2. Therefore, it can be known that when the K value is 3, the corresponding norm is the smallest, and the optimal modal component number in the normal state of the bearing is 3, which is also consistent with the estimated result of the line graph of the K value above.
Setting the optimal K = 3 as the number of modes that VMD needed to decompose, 10 groups of data in the normal state were decomposed in a cycle to obtain the spectrum diagram and time-domain feature set of the modal components after VMD decomposition.
Figure 4 represents the spectrum diagram of VMD decomposition when K = 3 under the normal state. From the spectrum diagram obtained, VMD avoids the defects of modal aliasing and endpoint effect of decomposition methods such as EMD. Figure 5 shows the characteristic signals extracted from the signal components under the normal state. In this paper, 16 time-domain indexes are used to reflect the features. The three states of wear, internal crack, and external crack are also obtained through these steps, and then, these feature data are put into an Excel sheet to form a feature set. The feature set is divided into training data and test data and input into the LSSVM toolbox for fault pattern recognition.
This paper made three contrast figures of fault diagnosis precision. They are, respectively, as follows: VMD was optimized and unoptimized, LSSVM was optimized and unoptimized, as well as the condition of the VMD and LSSVM was optimized and unoptimized. Figure 6(a) shows that the fault diagnosis accuracy of optimized VMD is 91.5%. Figure 6(b) shows that the diagnostic accuracy of unoptimized VMD is 88.3333%, and the contrast figure from this group that can validate the proposed VMD optimization method is effective; Figure 7(a) shows that the fault diagnosis accuracy of optimized LSSVM is 91.8333%. Compared with the result of 88.3333% without optimization, the optimization of LSSVM also has the effect of improving the accuracy. When VMD and LSSVM were optimized, it further improved the accuracy of fault diagnosis, as shown in Figure 8(a), as the accuracy was 92%. Let the optimized VMD be abbreviated as P-VMD, and the optimized LSSVM is abbreviated as P-LSSVM. Table 3 lists and illustrates the fault diagnosis accuracy of different algorithms and our algorithm.
(a)
(b)
(a)
(b)
(a)
(b)
From the comparison of three sets of results, we can clearly see that the proposed method in this paper based on instantaneous frequency optimization of the VMD fault diagnosis method is effective.
5. Conclusions
In this paper, the K value optimization problem of the variational modal decomposition algorithm is studied. Considering that the mode number K of VMD needs to be selected according to prior knowledge, improper selection will lead to overdecomposition or underdecomposition so that useful characteristic data cannot be extracted, ultimately leading to the problem of low accuracy of fault diagnosis. In this paper, the instantaneous frequency is used to find the optimal K value of VMD decomposition. Finally, the LSSVM model optimized by particle swarm optimization is combined to carry out fault pattern recognition. The results show the following:(1)The advantage of measuring the value of K by the change of the instantaneous frequency of the signal component after VMD decomposition is more accurate and simple than the previous observation method to judge the value of K, which can avoid overdecomposition and underdecomposition.(2)The optimized VMD decomposition algorithm can better reflect the characteristic parameters of vibration signals, which make subsequent feature extraction easier and helps to improve the diagnostic accuracy. As shown in the final experimental results, the accuracy of the optimized VMD is nearly 4% higher than that of the unoptimized results, indicating the effectiveness of this method.(3)The use of the PSO-LSSVM classification model for fault diagnosis can further improve the accuracy of the final diagnosis. This conclusion can be verified by Figures 7 and 8 in the final experimental results.
It can be seen that a joint fault diagnosis method based on optimized VMD and LSSVM proposed in this paper improves the accuracy of fault diagnosis.
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
The authors would like to thank the Guangdong Institute of Petrochemical Technology for providing the data set of the rolling bearing in this paper. This work was supported in part by the National Natural Science Foundation of P.R. China under Grant 61663008 and Chongqing Technology Innovation and Application Special Key Project under Grant cstc2019jscx-mbdxX0015.