# Intelligent Mechanical Fault Diagnosis Based on Multiwavelet Adaptive Threshold Denoising and MPSO

Authors: Hao Sun; Ke Li; Huaqing Wang; Peng Chen; Yi Cao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The condition diagnosis of rotating machinery depends largely on the feature analysis of vibration signals measured for the condition diagnosis. However, the signals measured from rotating machinery usually are nonstationary and nonlinear and contain noise. The useful fault features are hidden in the heavy background noise. In this paper, a novel fault diagnosis method for rotating machinery based on multiwavelet adaptive threshold denoising and mutation particle swarm optimization (MPSO) is proposed. Geronimo, Hardin, and Massopust (GHM) multiwavelet is employed for extracting weak fault features under background noise, and the method of adaptively selecting appropriate threshold for multiwavelet with energy ratio of multiwavelet coefficient is presented. The six nondimensional symptom parameters (SPs) in the frequency domain are defined to reflect the features of the vibration signals measured in each state. Detection index (DI) using statistical theory has been also defined to evaluate the sensitiveness of SP for condition diagnosis. MPSO algorithm with adaptive inertia weight adjustment and particle mutation is proposed for condition identification. MPSO algorithm effectively solves local optimum and premature convergence problems of conventional particle swarm optimization (PSO) algorithm. It can provide a more accurate estimate on fault diagnosis. Practical examples of fault diagnosis for rolling element bearings are given to verify the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Rolling element bearings are an important part of and widely used in rotating machinery. In practical application, bearing failures may cause the breakdown of equipment, and further, serious consequences may arise due to the failure. Thus, fault diagnosis and condition discrimination of bearings have an important significance for safe operation, guaranteeing production efficiency and reducing maintenance cost. Many reliability survey papers deal with failure statistics of rotating machinery subassemblies, focusing mainly on roller bearing because of their widespread use in industry [1–4]. Occurrence rate of bearing faults is very high in rotating machines, and other faults arising in rotation machines are often associated with bearing faults. In many instances, the accuracy of the instruments and devices used to monitor and control the rotation machines is highly dependent on the dynamic performance of bearings. Although fault diagnosis of rolling bearings is often artificially carried out using time or frequency analysis of vibration signals, there is a need for a reliable, fast automated diagnosis method.Vibration diagnosis is commonly used to detect the faults and identify the states in rotating machine. The condition diagnosis of rotating machinery depends largely on the feature analysis of vibration signals measured for the condition diagnosis because the signals carry dynamic information about the machine state [5–7]. However, feature extraction for fault diagnosis is difficult, because if the vibration signals are measured at an early stage of the machine failure, or at a location away from the fault part, the vibration signals contain strong noise. Stronger noise than the actual failure signal may lead to misrecognition of useful information for the condition diagnosis. Thus, it is important that the feature of the signal can be sensitively extracted at the state change of a machine.Wavelet transform (WT) is well known for its ability to focus on localized structures in time-frequency domain which has been widely used for fault diagnosis of rolling element bearings [8–10]. It has the local characteristic of time domain as well as frequency domain and its time-frequency window is changeable. In the processing of nonstationary signals it presents better performance than the traditional Fourier analysis. However, the measured signals often contain strong noise and the fault features are hidden in the background noise, and it is not the best way for WT to match the different fault features with a single wavelet and scaling functions, which will reduce the fault diagnosis accuracy. Multiwavelet transform is the new development of WT. It is constructed from translations and dilations of scaling and wavelet vector functions and has the predominant properties such as orthogonality, symmetry, compact support, and higher order vanishing moments. Multiwavelet transform decomposes the signal into subsignals of different frequency bands based on vector basis functions, via inner product principle. Because of the multiple scaling and wavelet basis functions, multiwavelet transform has predominant advantages in feature extraction of signals. Recently multiwavelet transform has been applied in fault diagnosis of rotating machinery as a powerful tool. In [11], multiwavelet system was introduced to diagnose gear faults. In [12, 13], multiwavelet lifting scheme was improved for compound faults separation and extraction. In [14], the undecimated multiwavelet was proposed for fault diagnosis of planetary gearboxes.PSO algorithm is a population based stochastic optimization technique developed by Kennedy and Eberhart in 1995 and inspired by social behavior of bird flocking or fish schooling [15]. In PSO algorithm, particles cooperate in finding good solutions for difficult discrete optimization problems. PSO algorithm has been applied to a variety of different problems, such as function optimization [16], scheduling [17], traveling salesman problem [18], neural network training [19, 20], and clustering task [21–23] which is the topic of interest in this paper. In recent years, PSO algorithm has been successfully applied in mechanical fault diagnosis; the domestic research on PSO fault diagnosis issues also has many articles reporting [24–27]. In [24], Bocaniala and Sa da Costa compared the time spent by PSO algorithm and genetic algorithm, testifying PSO algorithm with prominent superiority through fault diagnosis benchmark problem. In [25], Pan et al. used PSO algorithm to extract fault characteristics of rotation machinery. In [26, 27], PSO algorithm was used to diagnose gearbox fault. In this study, a clustering model is constructed by using an improved PSO called MPSO algorithm. It is used to classify the SPs calculated from the signals in each machine state for condition diagnosis, as well as obtaining their optimal clustering centers. According to these optimal clustering centers’ information, the conditions of the machine can be accurately identified.In order to extract the fault features of signals more effectively and identify mechanical condition more accurately, this paper proposes a novel fault diagnosis method for rotation machinery based on multiwavelet adaptive threshold denoising and MPSO algorithm. GHM multiwavelet is employed for extracting weak fault features under heavy background noise, and the method of adaptively selecting appropriate threshold values for multiwavelet with energy ratio of multiwavelet coefficient is presented. The six nondimensional SPs in the frequency domain are defined to reflect the features of the vibration signals measured in each state. DI using statistical theory has been also defined to evaluate the sensitiveness of SP for condition diagnosis. MPSO algorithm with adaptive inertia weight adjustment and particle mutation is proposed for condition identification. MPSO algorithm effectively solves local optimum and premature convergence problems of conventional particle swarm optimization (PSO) algorithm. It can provide a more accurate estimate on fault diagnosis. Practical examples of fault diagnosis for rolling element bearings are given to verify the effectiveness of the proposed method.


## 2. Feature Extraction by Multiwavelet Adaptive Threshold Denoising

### 2.1. Multiwavelet Theory

Multiwavelet consists of wavelet function vectorΨ and a function vector Φ is called multiscaling function. They are denoted as follows [28]:
(1)ϕ=[ϕ1,ϕ2,…,ϕr]T,ψ=[ψ1,ψ2,…,ψr]T.
For a multiresolution of multiplicity r>1.Similar to scalar wavelet,Ψ and Φ satisfy the two-scale matrix refinement equations:
(2)ϕ(t)=∑k=0NHKϕ(2t-k),ψ(t)=∑k=0NGKϕ(2t-k),
where k∈Z and Z is the set of integers. HK and GK are low-pass and high-pass matrix filter banks, respectively.Letcj-1=[c1,j-1,…,cr,j-1]T be the vector low frequency coefficients and let dj-1=[d1,j-1,…,dr,j-1]T be the vector high frequency coefficient; multiwavelet decomposition and composition are denoted as follows:
(3)cj-1,n=∑kHk-2ncj,k,dj-1,n=∑kGk-2ncj,k,cj,k=∑nHk-2n*cj,k+Gk-2n*dj-1,n,
where * means the complex conjugate transpose.Figure1 shows decomposition and reconstruction of multiwavelet.Decomposition and reconstruction of multiwavelet.
(a)
Decomposition
(b)
ReconstructionGHM multiwavelet constructed by Geronimo, Hardin, and Massopust is one of the most important multiwavelet systems with two pairs of scaling and wavelet functions and has the superior properties of short support, symmetry, orthogonality, and second approximation order [29]. Because of the excellent properties, GHM multiwavelets are adopted in this study. The multiscaling functions and multiwavelet functions of GHM multiwavelets are presented in Figure 2. The dilation and wavelet equations for GHM multiwavelet have four coefficients as follows:
(4)[ϕ1ϕ2]=∑kHk|ϕ1(2t-k)ϕ2(2t-k)|,H0=|310225-240-320|,H1=|3100924012|,H2=|009240-320|,H3=|00-2400|,[ψ1ψ2]=∑kGk|ψ1(2t-k)ψ2(2t-k)|,G0=|-240-320-120-3220|,G1=|9240-129200|,G2=|9240-320-920-3220|,G3=|-2400-1200|.Multiscaling functions and multiwavelet functions of GHM. (a) Scaling functionΦ1, (b) scaling function Φ2, (c) multiwavelet function Ψ1, and (d) multiwavelet function Ψ2.
(a)
(b)
(c)
(d)In view of the matrix filter banks, preprocessing is necessary to translate one stream input signal into two streams. Some preprocessing of preprocessing for multiwavelets has been proposed, such as repeated-row preprocessing, balanced multiwavelet, and prefilter methods [30]. Different preprocessing methods will produce different effect on performances of multiwavelets. It is a fundamental problem for each multiwavelet function to choose an appropriate preprocessing method for specific applications. In this study, the preprocessing method of repeated-row preferable for GHM multiwavelet is adopted and given as follows [31]:
(5)sn=|xncxn|,
where xn is original signal, sn is the signal after preprocessing, and c=2.


### 2.2. Multiwavelet Adaptive Threshold Denoising

Similar to single wavelet, multiwavelet denoising depends largely on the threshold denoising. The effect of threshold denoising depends on the selection of thresholds. A variety of threshold choosing methods can be mainly divided into two categories: global thresholding and level-dependent thresholding. The former chooses a single value ofλ to be applied globally to all empirical wavelet coefficients, while the latter chooses different threshold value λι for each wavelet level. However, it is difficult to choose appropriate threshold values for different wavelet coefficient. A large threshold value cuts too many coefficients, resulting in the loss of useful information. Conversely, a too small threshold value will leave much noise. In this study, the method of adaptive selecting appropriate threshold values for multiwavelet denoising based on comparison of noise energy in different levels is proposed. According to the noise levels of wavelet coefficients, the adaptive threshold value is determined by energy ratio.The energy spectrum of multiwavelet coefficient is denoted as follows:(6)Ej=∑i=1rEij=∑i=1r∑n=1m|di(n)|2,
where Ej represents multiwavelet coefficient total energy in the jth layer; Eij represents the multiwavelet coefficient energy of the i-dimensional in the jth layer; dj(n) is the multiwavelet coefficient in the jth layer after r-dimensional multiwavelet decomposition; and r is the number of dimensions of multiwavelet coefficient.Energy ratio of multiwavelet coefficient can be obtained as follows:(7)p=EijEj.According to noise level of multiwavelet coefficients, the threshold values of each multiwavelet coefficient can be adaptively obtained as follows:(8)μ=p×λij,λij=M×2ln⁡⁡(n)0.6745,
where M is the median absolute value of multiwavelet coefficient; n is signal length.In conclusion, the processing steps of multiwavelet adaptive threshold denoising are summarized as follows.(1)
Preprocess the original signal to transform it into two streams by the method of repeated-row preferable.(2)
Decompose two stream signals using multiwavelets.(3)
Threshold values are adaptively determined by energy ratio of the wavelet coefficients.(4)
Threshold the wavelet coefficients.(5)
Reconstruct the threshold wavelet coefficients and the scale coefficients.(6)
Postprocess the two stream results to get the denoising signal.In order to test effectiveness of multiwavelet adaptive threshold denoising proposed in this paper, a simulation experiment is designed as follows.The simulation signal is composed of a periodic impulse component and white Gaussian noise to simulate a bearing fault. The periodic impulse signal with the period of 0.01 s is expressed as(9)x(t)=x0e-ξωntsinωn1-ξ2t,
where ξ is damp coefficient; ωn denotes natural frequency; x0 indicates displacement constant. The shock impulse signal is displayed in Figure 3(a). In this case, ξ=0.1, ωn=3 kHz, x0=5, and sampling frequency and sampling points are 20 kHz and 4096, respectively. The simulation signal is shown in Figure 3(b), the signal has a low signal-noise ratio (SNR), and no useful features can be seen in the dynamic signal in the time domain.The simulation signal in the time domain: (a) the shock impulse signal; (b) the noisy signal.
(a)
(b)The noisy signal is processed using GHM multiwavelet adaptive threshold denoising, GHM multiwavelet neighboring coefficient denoising, and Daubechies 2 (db2) wavelet threshold denoising, respectively. Type of thresholding used is soft thresholding, and decomposition level is four. Denoised signal’s performance is evaluated based on mean square error (MSE) and SNR. Figure4 shows the denoising results using different wavelet denoising techniques. The SNR and MSE of different wavelet denoising techniques are calculated, as shown in Table 1. The results indicate that the method of GHM multiwavelet adaptive threshold denoising has the maximum SNR and the minimum MSE, which means the method proposed in this study can effectively extract the defect-induced shock impulses and eliminate much noise from the simulation signal.Table 1
SNR and MSE of different wavelet denoising techniques.
Multiwavelet adaptive threshold denoising
Multiwavelet neighboring coefficient denoising
db2 wavelet threshold denoising
SNR
14.516
11.098
9.667
MSE
0.212
0.235
0.306The denoising results using different wavelet denoising techniques. (a) Multiwavelet adaptive threshold denoising, (b) multiwavelet neighboring coefficient denoising, and (c) db2 wavelet threshold denoising.
(a)
(b)
(c)


## 2.1. Multiwavelet Theory

Multiwavelet consists of wavelet function vectorΨ and a function vector Φ is called multiscaling function. They are denoted as follows [28]:
(1)ϕ=[ϕ1,ϕ2,…,ϕr]T,ψ=[ψ1,ψ2,…,ψr]T.
For a multiresolution of multiplicity r>1.Similar to scalar wavelet,Ψ and Φ satisfy the two-scale matrix refinement equations:
(2)ϕ(t)=∑k=0NHKϕ(2t-k),ψ(t)=∑k=0NGKϕ(2t-k),
where k∈Z and Z is the set of integers. HK and GK are low-pass and high-pass matrix filter banks, respectively.Letcj-1=[c1,j-1,…,cr,j-1]T be the vector low frequency coefficients and let dj-1=[d1,j-1,…,dr,j-1]T be the vector high frequency coefficient; multiwavelet decomposition and composition are denoted as follows:
(3)cj-1,n=∑kHk-2ncj,k,dj-1,n=∑kGk-2ncj,k,cj,k=∑nHk-2n*cj,k+Gk-2n*dj-1,n,
where * means the complex conjugate transpose.Figure1 shows decomposition and reconstruction of multiwavelet.Decomposition and reconstruction of multiwavelet.
(a)
Decomposition
(b)
ReconstructionGHM multiwavelet constructed by Geronimo, Hardin, and Massopust is one of the most important multiwavelet systems with two pairs of scaling and wavelet functions and has the superior properties of short support, symmetry, orthogonality, and second approximation order [29]. Because of the excellent properties, GHM multiwavelets are adopted in this study. The multiscaling functions and multiwavelet functions of GHM multiwavelets are presented in Figure 2. The dilation and wavelet equations for GHM multiwavelet have four coefficients as follows:
(4)[ϕ1ϕ2]=∑kHk|ϕ1(2t-k)ϕ2(2t-k)|,H0=|310225-240-320|,H1=|3100924012|,H2=|009240-320|,H3=|00-2400|,[ψ1ψ2]=∑kGk|ψ1(2t-k)ψ2(2t-k)|,G0=|-240-320-120-3220|,G1=|9240-129200|,G2=|9240-320-920-3220|,G3=|-2400-1200|.Multiscaling functions and multiwavelet functions of GHM. (a) Scaling functionΦ1, (b) scaling function Φ2, (c) multiwavelet function Ψ1, and (d) multiwavelet function Ψ2.
(a)
(b)
(c)
(d)In view of the matrix filter banks, preprocessing is necessary to translate one stream input signal into two streams. Some preprocessing of preprocessing for multiwavelets has been proposed, such as repeated-row preprocessing, balanced multiwavelet, and prefilter methods [30]. Different preprocessing methods will produce different effect on performances of multiwavelets. It is a fundamental problem for each multiwavelet function to choose an appropriate preprocessing method for specific applications. In this study, the preprocessing method of repeated-row preferable for GHM multiwavelet is adopted and given as follows [31]:
(5)sn=|xncxn|,
where xn is original signal, sn is the signal after preprocessing, and c=2.


## 2.2. Multiwavelet Adaptive Threshold Denoising

Similar to single wavelet, multiwavelet denoising depends largely on the threshold denoising. The effect of threshold denoising depends on the selection of thresholds. A variety of threshold choosing methods can be mainly divided into two categories: global thresholding and level-dependent thresholding. The former chooses a single value ofλ to be applied globally to all empirical wavelet coefficients, while the latter chooses different threshold value λι for each wavelet level. However, it is difficult to choose appropriate threshold values for different wavelet coefficient. A large threshold value cuts too many coefficients, resulting in the loss of useful information. Conversely, a too small threshold value will leave much noise. In this study, the method of adaptive selecting appropriate threshold values for multiwavelet denoising based on comparison of noise energy in different levels is proposed. According to the noise levels of wavelet coefficients, the adaptive threshold value is determined by energy ratio.The energy spectrum of multiwavelet coefficient is denoted as follows:(6)Ej=∑i=1rEij=∑i=1r∑n=1m|di(n)|2,
where Ej represents multiwavelet coefficient total energy in the jth layer; Eij represents the multiwavelet coefficient energy of the i-dimensional in the jth layer; dj(n) is the multiwavelet coefficient in the jth layer after r-dimensional multiwavelet decomposition; and r is the number of dimensions of multiwavelet coefficient.Energy ratio of multiwavelet coefficient can be obtained as follows:(7)p=EijEj.According to noise level of multiwavelet coefficients, the threshold values of each multiwavelet coefficient can be adaptively obtained as follows:(8)μ=p×λij,λij=M×2ln⁡⁡(n)0.6745,
where M is the median absolute value of multiwavelet coefficient; n is signal length.In conclusion, the processing steps of multiwavelet adaptive threshold denoising are summarized as follows.(1)
Preprocess the original signal to transform it into two streams by the method of repeated-row preferable.(2)
Decompose two stream signals using multiwavelets.(3)
Threshold values are adaptively determined by energy ratio of the wavelet coefficients.(4)
Threshold the wavelet coefficients.(5)
Reconstruct the threshold wavelet coefficients and the scale coefficients.(6)
Postprocess the two stream results to get the denoising signal.In order to test effectiveness of multiwavelet adaptive threshold denoising proposed in this paper, a simulation experiment is designed as follows.The simulation signal is composed of a periodic impulse component and white Gaussian noise to simulate a bearing fault. The periodic impulse signal with the period of 0.01 s is expressed as(9)x(t)=x0e-ξωntsinωn1-ξ2t,
where ξ is damp coefficient; ωn denotes natural frequency; x0 indicates displacement constant. The shock impulse signal is displayed in Figure 3(a). In this case, ξ=0.1, ωn=3 kHz, x0=5, and sampling frequency and sampling points are 20 kHz and 4096, respectively. The simulation signal is shown in Figure 3(b), the signal has a low signal-noise ratio (SNR), and no useful features can be seen in the dynamic signal in the time domain.The simulation signal in the time domain: (a) the shock impulse signal; (b) the noisy signal.
(a)
(b)The noisy signal is processed using GHM multiwavelet adaptive threshold denoising, GHM multiwavelet neighboring coefficient denoising, and Daubechies 2 (db2) wavelet threshold denoising, respectively. Type of thresholding used is soft thresholding, and decomposition level is four. Denoised signal’s performance is evaluated based on mean square error (MSE) and SNR. Figure4 shows the denoising results using different wavelet denoising techniques. The SNR and MSE of different wavelet denoising techniques are calculated, as shown in Table 1. The results indicate that the method of GHM multiwavelet adaptive threshold denoising has the maximum SNR and the minimum MSE, which means the method proposed in this study can effectively extract the defect-induced shock impulses and eliminate much noise from the simulation signal.Table 1
SNR and MSE of different wavelet denoising techniques.
Multiwavelet adaptive threshold denoising
Multiwavelet neighboring coefficient denoising
db2 wavelet threshold denoising
SNR
14.516
11.098
9.667
MSE
0.212
0.235
0.306The denoising results using different wavelet denoising techniques. (a) Multiwavelet adaptive threshold denoising, (b) multiwavelet neighboring coefficient denoising, and (c) db2 wavelet threshold denoising.
(a)
(b)
(c)


## 3. Symptom Parameters for Fault Diagnosis and Sensitivity Evaluation

### 3.1. Symptom Parameters for Fault Diagnosis

When developing intelligent condition diagnosis system by computer, symptom parameters (SPs) are required to express the information indicated by a signal measured for diagnosing machinery faults. A good symptom parameter can correctly reflect states and the condition trend of plant machinery [32–34]. Many symptom parameters have been defined in the pattern recognition field. Here, six SPs in the frequency domain, commonly used for the fault diagnosis of plant machinery, are considered.Frequency-domain skewness:(10)P1=∑i=1I(fi-f-)3·F(fi)σ3I.Frequency-domain kurtosis:(11)P2=∑i=1I(fi-f-)4·F(fi)σ4·I.Mean frequency that wave shape cross the mean of time-domain signal:(12)P3=∑i=1Ifi4·F(fi)∑i=1Ifi2·F(fi).Stabilization factor of wave shape:(13)P4=∑i=1Ifi2·F(fi)∑i=1IF(fi)∑i=1Ifi4·F(fi).Sum of the squares of the power spectrum:(14)P5=∑i=1IF(fi).Square root of the sum of the squares of the power spectrum:(15)P6=∑i=1IF2(fi),
where I is the number of spectrum lines, fi is frequency, and from 0 Hz to the maximum analysis frequency, F(fi) is the power spectrum value at frequency fi, and i=1~I. f- is mean value of the analysis frequency, and f-=(∑i=1Ifi·F(fi))/∑i=1IF(fi); σ is standard deviation, and σ=(∑i=1I(fi-f-)2·F(fi))/I.


### 3.2. Detection Index

For automatic diagnosis, SPs are needed that can sensitively distinguish the fault types. In order to evaluate the sensitivity of a SP for distinguishing two states, such as a normal or an abnormal state, DI is defined as follows.Supposing thatx1 and x2 are the SP values calculated from the signals measured in state 1 and state 2, respectively, their average value and standard deviation are μ and σ. The DI is calculated by
(16)DI=|μ1-μ2|σ1+σ2.The distinction rate (DR) is defined as(17)DR=1-12π∫-∞-DIexp⁡(-μ22)dμ.It is obvious that the larger the value of the DI, the larger the value of the DR will be and, therefore, the better the SP will be. Thus, the DIcan be used as the index of the quality to evaluate the distinguishing sensitivity of the SP.The number of symptom parameters used for diagnosis and fault types areM and N, respectively; the synthetic detection index (SDI) is defined as follows:
(18)SDI=∑i=1N-1∑j=i+1N∑k=1M|μik-μjk|σik2+σjk2.


## 3.1. Symptom Parameters for Fault Diagnosis

When developing intelligent condition diagnosis system by computer, symptom parameters (SPs) are required to express the information indicated by a signal measured for diagnosing machinery faults. A good symptom parameter can correctly reflect states and the condition trend of plant machinery [32–34]. Many symptom parameters have been defined in the pattern recognition field. Here, six SPs in the frequency domain, commonly used for the fault diagnosis of plant machinery, are considered.Frequency-domain skewness:(10)P1=∑i=1I(fi-f-)3·F(fi)σ3I.Frequency-domain kurtosis:(11)P2=∑i=1I(fi-f-)4·F(fi)σ4·I.Mean frequency that wave shape cross the mean of time-domain signal:(12)P3=∑i=1Ifi4·F(fi)∑i=1Ifi2·F(fi).Stabilization factor of wave shape:(13)P4=∑i=1Ifi2·F(fi)∑i=1IF(fi)∑i=1Ifi4·F(fi).Sum of the squares of the power spectrum:(14)P5=∑i=1IF(fi).Square root of the sum of the squares of the power spectrum:(15)P6=∑i=1IF2(fi),
where I is the number of spectrum lines, fi is frequency, and from 0 Hz to the maximum analysis frequency, F(fi) is the power spectrum value at frequency fi, and i=1~I. f- is mean value of the analysis frequency, and f-=(∑i=1Ifi·F(fi))/∑i=1IF(fi); σ is standard deviation, and σ=(∑i=1I(fi-f-)2·F(fi))/I.


## 3.2. Detection Index

For automatic diagnosis, SPs are needed that can sensitively distinguish the fault types. In order to evaluate the sensitivity of a SP for distinguishing two states, such as a normal or an abnormal state, DI is defined as follows.Supposing thatx1 and x2 are the SP values calculated from the signals measured in state 1 and state 2, respectively, their average value and standard deviation are μ and σ. The DI is calculated by
(16)DI=|μ1-μ2|σ1+σ2.The distinction rate (DR) is defined as(17)DR=1-12π∫-∞-DIexp⁡(-μ22)dμ.It is obvious that the larger the value of the DI, the larger the value of the DR will be and, therefore, the better the SP will be. Thus, the DIcan be used as the index of the quality to evaluate the distinguishing sensitivity of the SP.The number of symptom parameters used for diagnosis and fault types areM and N, respectively; the synthetic detection index (SDI) is defined as follows:
(18)SDI=∑i=1N-1∑j=i+1N∑k=1M|μik-μjk|σik2+σjk2.


## 4. MPSO for Condition Diagnosis

### 4.1. Brief of PSO

PSO algorithm is based on the groups, and according to the environmental fitness, individual in groups will be moved to the good region. The algorithm evaluates the optimal result by using evolutionary fitness function of group, and each particle in the algorithm has a fitness value determined by the fitness function; two properties of position and speed that are used to show the position and moving speed of the current articles in the solving space, by the fitness function value corresponding to particle position coordinate, determine the performance of particles. In PSO algorithm, each particle adjusts its position according to its own experience and according to the experience of a neighboring particle, making use of the best position encountered by itself and its neighbor.In theR-dimensional search space, the i particle’s space position is defined as follows:
(19)P(i)·location=[Xi1,Xi2,…,XiR].The velocity of particlei is defined as follows:
(20)P(i)·velocity=[Vi1,Vi2,…,ViR].The best previous position of particlei is defined as follows:
(21)P(i)·best=[Pi1,Pi2,…,PiR].The best position among all particles experienced is defined as follows:(22)g(i)·best=[gi1,gi2,…,giR].The particle updates the position and velocity according to the following equations:(23)P(i)·velocity(t+1)=ωP(i)·velocity(t)+η1r1[P(i)·best(t)-P(i)·location(t)]+η2r2[g(i)best(t)-P(i)·location(t)],P(i)·location(t+1)=P(i)·location(t)+P(i)·velocity(t+1),
where r1 and r2 are the random numbers within (0,1) and η1 and η2 are the acceleration which constants the control of how far a particle moves in a single generation. The inertia weight ω controls the previous velocity of particle, and it is defined as follows:
(24)ω=0.5+rand2,
where rand is random generated number between 0 and 1.Although PSO algorithm is easy to realize, the method is easy to trap into local optimum. Shi and Eberhart proposed a linearly decreasing weight particle swarm optimization (WPSO) of which a linearly decreasing inertia factor was introduced into the velocity of the updated equation from the original PSO [35, 36]. The performance of WPSO is significantly improved over the original PSO because WPSO balances out the global and local search abilities of the swarm effectively. The equation for the linearly decreased weight is defined as follows:
(25)ωl=ωmax⁡-iteration×ωmax⁡-ωmin⁡iterationmax⁡,
where ωmax⁡ is 1, ωmin⁡ is 0.1, and iterationmax⁡ is the maximum number of the allowed iterations.The velocity of the updated equation for WPSO is defined as follows:(26)P(i)·velocity(t+1)=ωlP(i)·velocity(t)+η1r1[P(i)·best(t)-P(i)·location(t)]+η2r2[g(i)best(t)-P(i)·location(t)].


### 4.2. MPSO

Although WPSO algorithm improved conventional PSO to a certain extent, it cannot adapt to all of complex practical problems. The main reasons can be explained as follows. (1) The inertia weight of conventional WPSO algorithm is monotone decreasing, and adjustment ability of WPSO algorithm is limited. If particles cannot find optimal point in the initial stage of the algorithm, WPSO algorithm is easy to trap into local optimum with the decrease of the inertia weight. (2) With increasing iterations, particle diversity of WPSO algorithm decreases; it causes deterioration of global search ability; WPSO algorithm is also easy to trap into local optimum and premature convergence.To improve global search ability and adjustment ability of conventional PSO algorithm and prevent local optimum and premature convergence problems, MPSO algorithm with adaptive inertia weight adjustment and particle mutation is proposed in this paper.


#### 4.2.1. Adaptive Inertia Weight

Define change rate of fitness value:(27)R=|f(t+5)-f(t)||f(t)|,
where f(t) is optimum fitness value of the tth iteration; f(t+5) is optimum fitness value of the (t+5)th iteration; R indicates change rate of fitness value in five iterations.According to the variation ofR, the inertia weight ω adaptively adjusts as follows:
(28)ω={k1+0.5q,R>0.05,k2+0.5q,R≤0.05,
where q is a random number with a uniform probability within 0~1; k1 and k2 are parameters; k1>k2; the choice of k1 and k2 is determined experimentally; here k1=0.5 and k2=0.2. When R>0.05, the algorithm is in the exploration stage, and a large ω is beneficial to the algorithm’s convergence. When R≤0.05, the algorithm is in the development stage, and a small ω is beneficial to searching optimum point.


#### 4.2.2. Particle Mutation

To increase particle diversity of PSO algorithm, the method of particle mutation is proposed. In the operation process of PSO algorithm, if the best position among all particlesgbest does not change in a long time, some particles are mutated according to a certain probability. The execution process of the mutation for PSO is as follows.(1)
All particles are arranged in ascending order according to the values of the fitness function.(2)
Them (m>1) particles with smaller fitness functions are selected.(3)
Random datari  {i=1,2…,m} for selected particles are produced automatically.(4)
A weightPm is set, and 0.1<P<0.5.(5)
Pm is compared with ri, if Pm>ri, and then the particle’s space position is updated by using (29).(6)
Steps (3)–(5) are looped until the space position ofm particles are updated:
(29)xijt+1=xijt(1+0.5η),where η is random data that obeys Gaussian(0,1) distribution.


### 4.3. Fitness Function of MPSO for Condition Diagnosis

Assume thatN is the sample set of vibration signals measured in m different states; the length of N is n, N={s1,s2,…,sn}. Every sample signal has t identified symptoms (in this paper, the symptoms are P1–P6). Then, the clustering analysis is to divide n sample data into m states, such that the fitness function F shown in (30) is minimized:
(30)min⁡F=∑j=1m∑i=1n∑k=1taij∥Sik-Xjk∥2,(31)Xjk=∑i=1naijSik∑i=1naij(j=1,2,…,m;k=1,2,…,t),(32)aij={1,ifSi∈statej0,ifSi∉statej(i=1,2,…,n;j=1,2,…,m).


## 4.1. Brief of PSO

PSO algorithm is based on the groups, and according to the environmental fitness, individual in groups will be moved to the good region. The algorithm evaluates the optimal result by using evolutionary fitness function of group, and each particle in the algorithm has a fitness value determined by the fitness function; two properties of position and speed that are used to show the position and moving speed of the current articles in the solving space, by the fitness function value corresponding to particle position coordinate, determine the performance of particles. In PSO algorithm, each particle adjusts its position according to its own experience and according to the experience of a neighboring particle, making use of the best position encountered by itself and its neighbor.In theR-dimensional search space, the i particle’s space position is defined as follows:
(19)P(i)·location=[Xi1,Xi2,…,XiR].The velocity of particlei is defined as follows:
(20)P(i)·velocity=[Vi1,Vi2,…,ViR].The best previous position of particlei is defined as follows:
(21)P(i)·best=[Pi1,Pi2,…,PiR].The best position among all particles experienced is defined as follows:(22)g(i)·best=[gi1,gi2,…,giR].The particle updates the position and velocity according to the following equations:(23)P(i)·velocity(t+1)=ωP(i)·velocity(t)+η1r1[P(i)·best(t)-P(i)·location(t)]+η2r2[g(i)best(t)-P(i)·location(t)],P(i)·location(t+1)=P(i)·location(t)+P(i)·velocity(t+1),
where r1 and r2 are the random numbers within (0,1) and η1 and η2 are the acceleration which constants the control of how far a particle moves in a single generation. The inertia weight ω controls the previous velocity of particle, and it is defined as follows:
(24)ω=0.5+rand2,
where rand is random generated number between 0 and 1.Although PSO algorithm is easy to realize, the method is easy to trap into local optimum. Shi and Eberhart proposed a linearly decreasing weight particle swarm optimization (WPSO) of which a linearly decreasing inertia factor was introduced into the velocity of the updated equation from the original PSO [35, 36]. The performance of WPSO is significantly improved over the original PSO because WPSO balances out the global and local search abilities of the swarm effectively. The equation for the linearly decreased weight is defined as follows:
(25)ωl=ωmax⁡-iteration×ωmax⁡-ωmin⁡iterationmax⁡,
where ωmax⁡ is 1, ωmin⁡ is 0.1, and iterationmax⁡ is the maximum number of the allowed iterations.The velocity of the updated equation for WPSO is defined as follows:(26)P(i)·velocity(t+1)=ωlP(i)·velocity(t)+η1r1[P(i)·best(t)-P(i)·location(t)]+η2r2[g(i)best(t)-P(i)·location(t)].


## 4.2. MPSO

Although WPSO algorithm improved conventional PSO to a certain extent, it cannot adapt to all of complex practical problems. The main reasons can be explained as follows. (1) The inertia weight of conventional WPSO algorithm is monotone decreasing, and adjustment ability of WPSO algorithm is limited. If particles cannot find optimal point in the initial stage of the algorithm, WPSO algorithm is easy to trap into local optimum with the decrease of the inertia weight. (2) With increasing iterations, particle diversity of WPSO algorithm decreases; it causes deterioration of global search ability; WPSO algorithm is also easy to trap into local optimum and premature convergence.To improve global search ability and adjustment ability of conventional PSO algorithm and prevent local optimum and premature convergence problems, MPSO algorithm with adaptive inertia weight adjustment and particle mutation is proposed in this paper.


### 4.2.1. Adaptive Inertia Weight

Define change rate of fitness value:(27)R=|f(t+5)-f(t)||f(t)|,
where f(t) is optimum fitness value of the tth iteration; f(t+5) is optimum fitness value of the (t+5)th iteration; R indicates change rate of fitness value in five iterations.According to the variation ofR, the inertia weight ω adaptively adjusts as follows:
(28)ω={k1+0.5q,R>0.05,k2+0.5q,R≤0.05,
where q is a random number with a uniform probability within 0~1; k1 and k2 are parameters; k1>k2; the choice of k1 and k2 is determined experimentally; here k1=0.5 and k2=0.2. When R>0.05, the algorithm is in the exploration stage, and a large ω is beneficial to the algorithm’s convergence. When R≤0.05, the algorithm is in the development stage, and a small ω is beneficial to searching optimum point.


### 4.2.2. Particle Mutation

To increase particle diversity of PSO algorithm, the method of particle mutation is proposed. In the operation process of PSO algorithm, if the best position among all particlesgbest does not change in a long time, some particles are mutated according to a certain probability. The execution process of the mutation for PSO is as follows.(1)
All particles are arranged in ascending order according to the values of the fitness function.(2)
Them (m>1) particles with smaller fitness functions are selected.(3)
Random datari  {i=1,2…,m} for selected particles are produced automatically.(4)
A weightPm is set, and 0.1<P<0.5.(5)
Pm is compared with ri, if Pm>ri, and then the particle’s space position is updated by using (29).(6)
Steps (3)–(5) are looped until the space position ofm particles are updated:
(29)xijt+1=xijt(1+0.5η),where η is random data that obeys Gaussian(0,1) distribution.


## 4.2.1. Adaptive Inertia Weight

Define change rate of fitness value:(27)R=|f(t+5)-f(t)||f(t)|,
where f(t) is optimum fitness value of the tth iteration; f(t+5) is optimum fitness value of the (t+5)th iteration; R indicates change rate of fitness value in five iterations.According to the variation ofR, the inertia weight ω adaptively adjusts as follows:
(28)ω={k1+0.5q,R>0.05,k2+0.5q,R≤0.05,
where q is a random number with a uniform probability within 0~1; k1 and k2 are parameters; k1>k2; the choice of k1 and k2 is determined experimentally; here k1=0.5 and k2=0.2. When R>0.05, the algorithm is in the exploration stage, and a large ω is beneficial to the algorithm’s convergence. When R≤0.05, the algorithm is in the development stage, and a small ω is beneficial to searching optimum point.


## 4.2.2. Particle Mutation

To increase particle diversity of PSO algorithm, the method of particle mutation is proposed. In the operation process of PSO algorithm, if the best position among all particlesgbest does not change in a long time, some particles are mutated according to a certain probability. The execution process of the mutation for PSO is as follows.(1)
All particles are arranged in ascending order according to the values of the fitness function.(2)
Them (m>1) particles with smaller fitness functions are selected.(3)
Random datari  {i=1,2…,m} for selected particles are produced automatically.(4)
A weightPm is set, and 0.1<P<0.5.(5)
Pm is compared with ri, if Pm>ri, and then the particle’s space position is updated by using (29).(6)
Steps (3)–(5) are looped until the space position ofm particles are updated:
(29)xijt+1=xijt(1+0.5η),where η is random data that obeys Gaussian(0,1) distribution.


## 4.3. Fitness Function of MPSO for Condition Diagnosis

Assume thatN is the sample set of vibration signals measured in m different states; the length of N is n, N={s1,s2,…,sn}. Every sample signal has t identified symptoms (in this paper, the symptoms are P1–P6). Then, the clustering analysis is to divide n sample data into m states, such that the fitness function F shown in (30) is minimized:
(30)min⁡F=∑j=1m∑i=1n∑k=1taij∥Sik-Xjk∥2,(31)Xjk=∑i=1naijSik∑i=1naij(j=1,2,…,m;k=1,2,…,t),(32)aij={1,ifSi∈statej0,ifSi∉statej(i=1,2,…,n;j=1,2,…,m).


## 5. Diagnosis and Application

In this section, the application of condition diagnosis for a rolling bearing is shown to verify that the method proposed in this paper is effective.


### 5.1. Experimental System

Figure5 shows the experimental system for a roller bearing fault diagnosis test. The most commonly occurring faults in a roller element bearing are the outer-race defect, the inner-race defect, and the roller element defect. These fault bearings are shown in Figure 6 and were created artificially using a wire-cutting machine. In this work an accelerometer (PCB MA352A60) was used to measure the vibration signals of the vertical direction in the normal, the outer-race defect, the inner-race defect, and the roller element defect states, respectively. The original vibration signals in each state are measured at a constant speed (800 rpm), and a 150 kg load was also transported on the rotating shaft by the loading equipment (RCS2-RA13R) while the vibration signals were being measured. The sampling frequency of the signal measurement was 50 kHz, and the sampling time was 20 s. All of the data was divided to 100 parts; 40 parts were used to train diagnosis system; other parts were used for condition identification test. Spectrum values at frequency Figures 7(a), 8(a), 9(a), and 10(a) show the original vibration signal in each state, and Figures 7(b), 8(b), 9(b), and 10(b) show the multiwavelet adaptive threshold denoising results of the vibration signal in each state.Figure 5
Experimental system for bearing fault diagnosis.Bearing defects. (a) Outer-race defect. (b) Inner-race defect. (c) Roller defect.
(a)
(b)
(c)The vibration signal in normal state: (a) the original vibration signal, (b) after multiwavelet adaptive threshold denoising, and (c) Fourier spectrum of the denoised signal.
(a)
(b)
(c)The vibration signal in outer-race defect state: (a) original vibration signal, (b) after multiwavelet adaptive threshold denoising, and (c) Fourier spectrum of the denoised signal.
(a)
(b)
(c)The vibration signal in inner-race defect state: (a) original vibration signal, (b) after multiwavelet adaptive threshold denoising, and (c) Fourier spectrum of the denoised signal.
(a)
(b)
(c)The vibration signal in roller defect state: (a) original vibration signal, (b) after multiwavelet adaptive threshold denoising, and (c) Fourier spectrum of the denoised signal.
(a)
(b)
(c)


### 5.2. Diagnosis by the Proposed Method

The main procedure for fault diagnosis using GHM multiwavelet adaptive threshold denoising and MPSO algorithm is shown in Figure11 and explained as follows.(1)
Vibration signals are measured in each known state.(2)
Weak fault feature is extracted by using GHM multiwavelet adaptive threshold denoising.(3)
SPs are calculated using (10)–(15).(4)
The highly sensitive SPs are selected for condition diagnosis by DI.(5)
MPSO algorithm is trained with SPs selected by DI, and the optimal clustering centers are obtained.(6)
Condition of the bearing can be diagnosed by the trained MPSO algorithm and SPs.Figure 11
Flowchart for the condition diagnosis by the method proposed in this study.In this study, the good SPs which have high sensitivity for distinguishing each fault state of the bearing are selected by the method of DI. As an example, Table2 lists parts of DIs of SPs. The maximum value (50.49) of SDI is obtained in the case of the combination of P1, P5, and P6, and when P1, P5, and P6, are used for distinguishing each state separately, the DIs are larger than 1.41, and all of the DRs are larger than 92.1%. Therefore, P1, P5, and P6 have high sensitivity for distinguishing each fault state of the bearing.Table 2
DIs of each SP.
P
1
P
2
P
3
P
4
P
5
P
6
DIN:O
5.733
2.607
0.953
13.973
6.920
3.287
DIN:I
1.947
1.467
0.740
1.593
3.387
1.840
DIN:R
4.540
3.580
0.513
0.707
2.287
2.820
DIO:I
3.793
0.467
0.587
1.367
2.413
1.680
DIO:R
2.007
0.693
1.040
1.533
1.567
1.606
DII:R
1.560
0.467
0.813
1.087
1.687
1.413In this study, MPSO automatically obtains the optimal clustering centers according to the classification of the sample data information. The purpose of training MPSO is the acquisition of optimum clustering centers. The SPs selected by DI were input into MPSO. MPSO converged to the optimum clustering centers. In the training process of MPSO, at first, the sample data are classified into the normal, the outer-race defect, the inner-race defect, and the roller element defect randomly. The fitness values and the clustering centers are calculated by (30) and (31). With increasing iterations, the speed and position (classification of the sample data) of the particle are updated incessantly, and according to the classification of the sample data information, the clustering centers are also updated. Finally, the optimal clustering centers with a minimum fitness value are calculated.To explain the effectiveness of MPSO algorithm, a comparison is made among MPSO, WPSO, and PSO algorithms. The optimal clustering centers of each state are obtained by MPSO, WPSO, and PSO algorithms, respectively. Particle number and iteration number of the three algorithms are 50 and 1000, respectively. The clustering centers obtained by each method are shown in Tables3, 4, and 5. Figure 12 shows the fitness curve of PSO, WPSO, and MPSO algorithms. It is obvious that MPSO algorithm has the minimum fitness value; namely, the optimal clustering centers obtained by MPSO algorithm are the most accurate. MPSO algorithm has stronger capacity of searching optimal solution.Table 3
Clustering centers obtained by MSPO.
Machinery condition
Clustering centers
P
1
P
5
P
6
Normal
0.378
2.882
0.076
Outer-race defect
0.433
6.094
0.571
Inner-race defect
0.118
4.051
0.187
Roller element defect
0.729
5.246
0.022Table 4
Clustering centers obtained by SPO.
Machinery condition
Clustering centers
P
1
P
5
P
6
Normal
0.371
2.922
0.069
Outer-race defect
0.583
6.081
0.558
Inner-race defect
0.082
3.825
0.175
Roller element defect
0.579
5.325
0.020Table 5
Clustering centers obtained by WSPO.
Machinery Condition
Clustering Centers
P
1
P
5
P
6
Normal
0.395
2.831
0.073
Outer-race defect
0.558
6.105
0.583
Inner-race defect
0.155
4.228
0.223
Roller element defect
0.796
5.045
0.021Figure 12
The fitness curve of PSO, WPSO, and MPSO algorithms.After training MPSO algorithm, to verify the diagnostic capability of the proposed method in this paper, the test data measured in each known state that had not been used to train MPSO algorithm were used. When inputting the test data into the trained MPSO algorithm, MPSO classified the test data according to the information of the optimum clustering centers shown in Table3 and correctly and quickly output identification results. As an example, some diagnosis results are listed in Table 6. We also identified condition of the rolling bearing using PSO and WPSO algorithms, respectively, and some identification results are shown in Tables 7 and 8. The comparison of diagnostic capability of each method is shown in Figure 13. Viewing the overall diagnostic results, diagnostic accuracy of each state using MPSO algorithm is 100%, 88.3%, 86.7%, and 81.7%, respectively; they are the largest in three methods. The method proposed in this study provides a more accurate estimate in the case of the rolling bearing faults diagnosis. These results verified the efficiency of the intelligent diagnosis method using multiwavelet adaptive threshold denoising and MPSO proposed in this paper.Table 6
Diagnosis result using proposed method.
Machinery condition
Number of training data
Number of test data
Number of  correct results
Diagnostic accuracy (%)
Normal
40
60
60
100%
Outer-race defect
40
60
53
88.3%
Inner-race defect
40
60
52
86.7%
Roller element defect
40
60
49
81.7%Table 7
Diagnostic result using PSO.
Machinery condition
Number of training data
Number of test data
Number of  correct results
Diagnostic accuracy (%)
Normal
40
60
52
86.7%
Outer-race defect
40
60
43
71.6%
Inner-race defect
40
60
41
68.3%
Roller element defect
40
60
43
71.6%Table 8
Diagnosis result using WPSO.
Machinery condition
Number of training data
Number of test data
Number of  correct results
Diagnostic accuracy (%)
Normal
40
60
60
100%
Outer-race defect
40
60
49
81.7%
Inner-race defect
40
60
45
75%
Roller element defect
40
60
46
76.7%Figure 13
Comparison of diagnostic capability.


## 5.1. Experimental System

Figure5 shows the experimental system for a roller bearing fault diagnosis test. The most commonly occurring faults in a roller element bearing are the outer-race defect, the inner-race defect, and the roller element defect. These fault bearings are shown in Figure 6 and were created artificially using a wire-cutting machine. In this work an accelerometer (PCB MA352A60) was used to measure the vibration signals of the vertical direction in the normal, the outer-race defect, the inner-race defect, and the roller element defect states, respectively. The original vibration signals in each state are measured at a constant speed (800 rpm), and a 150 kg load was also transported on the rotating shaft by the loading equipment (RCS2-RA13R) while the vibration signals were being measured. The sampling frequency of the signal measurement was 50 kHz, and the sampling time was 20 s. All of the data was divided to 100 parts; 40 parts were used to train diagnosis system; other parts were used for condition identification test. Spectrum values at frequency Figures 7(a), 8(a), 9(a), and 10(a) show the original vibration signal in each state, and Figures 7(b), 8(b), 9(b), and 10(b) show the multiwavelet adaptive threshold denoising results of the vibration signal in each state.Figure 5
Experimental system for bearing fault diagnosis.Bearing defects. (a) Outer-race defect. (b) Inner-race defect. (c) Roller defect.
(a)
(b)
(c)The vibration signal in normal state: (a) the original vibration signal, (b) after multiwavelet adaptive threshold denoising, and (c) Fourier spectrum of the denoised signal.
(a)
(b)
(c)The vibration signal in outer-race defect state: (a) original vibration signal, (b) after multiwavelet adaptive threshold denoising, and (c) Fourier spectrum of the denoised signal.
(a)
(b)
(c)The vibration signal in inner-race defect state: (a) original vibration signal, (b) after multiwavelet adaptive threshold denoising, and (c) Fourier spectrum of the denoised signal.
(a)
(b)
(c)The vibration signal in roller defect state: (a) original vibration signal, (b) after multiwavelet adaptive threshold denoising, and (c) Fourier spectrum of the denoised signal.
(a)
(b)
(c)


## 5.2. Diagnosis by the Proposed Method

The main procedure for fault diagnosis using GHM multiwavelet adaptive threshold denoising and MPSO algorithm is shown in Figure11 and explained as follows.(1)
Vibration signals are measured in each known state.(2)
Weak fault feature is extracted by using GHM multiwavelet adaptive threshold denoising.(3)
SPs are calculated using (10)–(15).(4)
The highly sensitive SPs are selected for condition diagnosis by DI.(5)
MPSO algorithm is trained with SPs selected by DI, and the optimal clustering centers are obtained.(6)
Condition of the bearing can be diagnosed by the trained MPSO algorithm and SPs.Figure 11
Flowchart for the condition diagnosis by the method proposed in this study.In this study, the good SPs which have high sensitivity for distinguishing each fault state of the bearing are selected by the method of DI. As an example, Table2 lists parts of DIs of SPs. The maximum value (50.49) of SDI is obtained in the case of the combination of P1, P5, and P6, and when P1, P5, and P6, are used for distinguishing each state separately, the DIs are larger than 1.41, and all of the DRs are larger than 92.1%. Therefore, P1, P5, and P6 have high sensitivity for distinguishing each fault state of the bearing.Table 2
DIs of each SP.
P
1
P
2
P
3
P
4
P
5
P
6
DIN:O
5.733
2.607
0.953
13.973
6.920
3.287
DIN:I
1.947
1.467
0.740
1.593
3.387
1.840
DIN:R
4.540
3.580
0.513
0.707
2.287
2.820
DIO:I
3.793
0.467
0.587
1.367
2.413
1.680
DIO:R
2.007
0.693
1.040
1.533
1.567
1.606
DII:R
1.560
0.467
0.813
1.087
1.687
1.413In this study, MPSO automatically obtains the optimal clustering centers according to the classification of the sample data information. The purpose of training MPSO is the acquisition of optimum clustering centers. The SPs selected by DI were input into MPSO. MPSO converged to the optimum clustering centers. In the training process of MPSO, at first, the sample data are classified into the normal, the outer-race defect, the inner-race defect, and the roller element defect randomly. The fitness values and the clustering centers are calculated by (30) and (31). With increasing iterations, the speed and position (classification of the sample data) of the particle are updated incessantly, and according to the classification of the sample data information, the clustering centers are also updated. Finally, the optimal clustering centers with a minimum fitness value are calculated.To explain the effectiveness of MPSO algorithm, a comparison is made among MPSO, WPSO, and PSO algorithms. The optimal clustering centers of each state are obtained by MPSO, WPSO, and PSO algorithms, respectively. Particle number and iteration number of the three algorithms are 50 and 1000, respectively. The clustering centers obtained by each method are shown in Tables3, 4, and 5. Figure 12 shows the fitness curve of PSO, WPSO, and MPSO algorithms. It is obvious that MPSO algorithm has the minimum fitness value; namely, the optimal clustering centers obtained by MPSO algorithm are the most accurate. MPSO algorithm has stronger capacity of searching optimal solution.Table 3
Clustering centers obtained by MSPO.
Machinery condition
Clustering centers
P
1
P
5
P
6
Normal
0.378
2.882
0.076
Outer-race defect
0.433
6.094
0.571
Inner-race defect
0.118
4.051
0.187
Roller element defect
0.729
5.246
0.022Table 4
Clustering centers obtained by SPO.
Machinery condition
Clustering centers
P
1
P
5
P
6
Normal
0.371
2.922
0.069
Outer-race defect
0.583
6.081
0.558
Inner-race defect
0.082
3.825
0.175
Roller element defect
0.579
5.325
0.020Table 5
Clustering centers obtained by WSPO.
Machinery Condition
Clustering Centers
P
1
P
5
P
6
Normal
0.395
2.831
0.073
Outer-race defect
0.558
6.105
0.583
Inner-race defect
0.155
4.228
0.223
Roller element defect
0.796
5.045
0.021Figure 12
The fitness curve of PSO, WPSO, and MPSO algorithms.After training MPSO algorithm, to verify the diagnostic capability of the proposed method in this paper, the test data measured in each known state that had not been used to train MPSO algorithm were used. When inputting the test data into the trained MPSO algorithm, MPSO classified the test data according to the information of the optimum clustering centers shown in Table3 and correctly and quickly output identification results. As an example, some diagnosis results are listed in Table 6. We also identified condition of the rolling bearing using PSO and WPSO algorithms, respectively, and some identification results are shown in Tables 7 and 8. The comparison of diagnostic capability of each method is shown in Figure 13. Viewing the overall diagnostic results, diagnostic accuracy of each state using MPSO algorithm is 100%, 88.3%, 86.7%, and 81.7%, respectively; they are the largest in three methods. The method proposed in this study provides a more accurate estimate in the case of the rolling bearing faults diagnosis. These results verified the efficiency of the intelligent diagnosis method using multiwavelet adaptive threshold denoising and MPSO proposed in this paper.Table 6
Diagnosis result using proposed method.
Machinery condition
Number of training data
Number of test data
Number of  correct results
Diagnostic accuracy (%)
Normal
40
60
60
100%
Outer-race defect
40
60
53
88.3%
Inner-race defect
40
60
52
86.7%
Roller element defect
40
60
49
81.7%Table 7
Diagnostic result using PSO.
Machinery condition
Number of training data
Number of test data
Number of  correct results
Diagnostic accuracy (%)
Normal
40
60
52
86.7%
Outer-race defect
40
60
43
71.6%
Inner-race defect
40
60
41
68.3%
Roller element defect
40
60
43
71.6%Table 8
Diagnosis result using WPSO.
Machinery condition
Number of training data
Number of test data
Number of  correct results
Diagnostic accuracy (%)
Normal
40
60
60
100%
Outer-race defect
40
60
49
81.7%
Inner-race defect
40
60
45
75%
Roller element defect
40
60
46
76.7%Figure 13
Comparison of diagnostic capability.


## 6. Conclusions

In order to diagnose faults of rotation machinery at an early stage, this paper proposed a novel intelligent condition diagnosis method using multiwavelet adaptive threshold denoising and MPSO to detect faults and distinguish fault types at an early stage. The main conclusions are summarized as follows.(1)
The method of multiwavelet adaptive threshold denoising was presented for extracting weak fault features under background noise. It could adaptively select appropriate threshold for multiwavelet with energy ratio of multiwavelet coefficient. The simulation experiment verified that the method of multiwavelet adaptive threshold denoising can effectively extract fault features and eliminate much noise from the noisy signal.(2)
The six SPs in the frequency domain were defined for reflecting the features of vibration signals measured in each state. DI using statistical theory had been also defined to evaluate the applicability of the SPs for the condition diagnosis measured in each state. DI could be used to indicate the fitness of a SP for condition identification.(3)
MPSO algorithm with adaptive inertia weight adjustment and particle mutation was proposed for condition identification. MPSO algorithm was used to classify the SPs calculated from the signals in each machine state for condition diagnosis, as well as obtaining their optimal clustering centers. According to these optimal clustering centers’ information, the conditions of rotation machinery could be accurately identified. MPSO algorithm effectively solved local optimum and premature convergence problems of conventional PSO algorithm and raised diagnostic accuracy.(4)
Practical example of condition diagnosis for a rolling bearing verified that the method proposed in this paper was effective. Moreover, a comparison was also made among MPSO, WPSO, and PSO algorithms. The diagnostic results show that MPSO algorithm could provide a more accurate estimate in the case of the rolling bearing faults diagnosis.



---

*Source: 142795-2014-07-22.xml*