Automatic Decision Method of Optimum Symptom Parameters and Frequency Bands for Intelligent Machinery Diagnosis: Application to Condition Diagnosis of Centrifugal Pump System

Abstract

A combination method of statistical filtering (SF) and the ant colony optimization (ACO) is proposed for automatic decision of optimum symptom parameters and frequency bands for machinery diagnosis. The noise of vibration signals is canceled by using SF. Similarity factor I_pq is defined to evaluate the filtering performance; the significance level α is optimized by genetic algorithms (GA). The optimum symptom parameters in different frequency bands, by which the states of rotating machinery can be sensitively distinguished, are automatically and sequentially selected by ACO based on the Mahalanobis distance between different machine states. Finally, the Mahalanobis distance is used to identify failure types based on the sequential diagnostic method. The new method proposed in this paper has been used to diagnose a centrifugal pump system for faults which often occur in the pump, such as impeller unbalance, shaft misalignment, and cavitation. The verification results of the condition diagnosis for a centrifugal pump show that the new method has good performance.

1. Introduction

When the vibration signal is used for condition diagnosis of rotating machinery and measured at an early stage of machine failure or at a distant location from the faulty part of the machine, the optimum decision of symptom parameters and distinction of the fault types are difficult, because the failure signal is strongly contaminated by noise. It is important to cancel the noise from the vibration signal as cleanly as possible in order to increase the sensitivity of the condition diagnosis [1–3]. For noise cancelling, many methods have been proposed, for example, band pass filter [4], adaptive filter [5], Wiener filter [6], Kalman filter [7], and so forth. However in the field of machinery diagnosis, these methods cannot always be applied to failure signal extraction, due to the following adverse conditions. Firstly, in the case of the band pass filter, the wide band noise cannot be cancelled. Secondly, when applying the adaptive filter for noise cancelling, the reference noise must be simultaneously measured with the signal. Simultaneous measurement of the reference noise is not easily realized in most cases of fault diagnosis. Finally, the noise cannot be effectively removed by the Wiener filter and Kalman filter if noise and signal do not follow the normal distribution. In order to overcome the problems above, we propose a new extraction method of the fault signal from the vibration signal measured in the abnormal state of a machine by using genetic algorithms (GA) and statistical information.

Rotating machines are used in a variety of applications. To increase the availability and safety of them, continuous and effective condition monitoring is required. Up to now, a lot of intelligent diagnosis methods were presented for this purpose. Particularly artificial neural networks, fuzzy logic, and failure tree were employed for condition monitoring and diagnosis of rotating machine [8, 9]. When an intelligent diagnosing method is used for diagnosis of plant machinery, symptom parameters are required to express the information indicated by the signal measured for diagnosing machinery faults [10]. However at present, there are few accepted methods for deciding the optimum combination of symptom parameters, by which many machine faults can be sensitively and precisely detected. In many cases, this has been done merely by trial and error.

As one of the intelligent algorithms, ant colony optimization (ACO) is tending to develop towards pattern recognition and parameter optimization domain and gaining popularity. ACO is used for optimally clustering $N$ objects into $K$ clusters based on Euclidean distance and compared with other stochastic algorithms in [11]. Li et al. and Song et al. used ACO as clustering system to identify the states of roller bearing [12, 13]. Taher Niknam combined fuzzy adaptive particle swarm optimization (PSO), ACO, and k-means algorithms to find better cluster partition. ACO is used to make sure that the global best position is unique for every particle [14]. Zhang et al. used ACO to optimize the parameter of SVM applied to machinery fault diagnosis [15]. And a hybrid method based on ACO and ANNs to address feature selection in the inductive learning of the neural network pattern classifier used for medical diagnosis is proposed by Sivagaminathan and Ramakrishnan. ANNs are used as the objective function of ACO, but, as a disadvantage, this maybe leads to a slow running speed of the program [16]. In my past study, the input symptom parameters of the clustering system are chosen by distinction index [17], while the ACO clustering system is used to identify the types of the signal [13]. However, research applying ACO to the selection of optimum symptom parameters based on Mahalanobis distance is still undiscovered. In this paper, we propose a new method by which the optimum combination of symptom parameters can be automatically and sequentially decided by using the ACO and Mahalanobis distance.

To sum up, in this paper the combination method of statistical filtering with GA and the ant colony optimization with Mahalanobis distance is proposed for deciding optimum symptom parameters and frequency bands automatically, and the Mahalanobis distance is used as the objective function for selecting the symptom parameter. On the other hand, in order to adapt to different failure types, the signals in different frequency bands are extracted by band pass filter before using the statistical filter (SF). The method proposed in this paper was used to diagnose a centrifugal pump system for distinguishing faults which often occur in the pump, such as impeller unbalance, shaft misalignment, and cavitation, as verifying example.

The procedure of the condition diagnosis using the proposed method is follows:

(1)

measuring vibration signals in each state;

(2)

dividing the vibration signals into three bands from low frequency to high frequency, as 0 Hz~1000 Hz, 1000 Hz~5000 Hz, and 5000 Hz~25000 Hz, respectively;

(3)

automatically extracting the feature spectra in each frequency band by using spectrum statistical filter and genetic algorithms (GA);

(4)

calculating the value of symptom parameters in time domain and frequency domain;

(5)

automatically and sequentially deciding the optimum symptom parameters and frequency bands for distinguishing each state by using ant colony optimization (ACO) and Mahalanobis distance;

(6)

verifying the effectiveness of the proposed method with sequential diagnosis method and new vibration signals measured in each state.

In Section 2, the principle of the proposed method is described. Section 3 describes the experimental setup and procedures. After that, the specific experimental process and results are presented in Section 4. Finally, Section 5 states conclusions.

2. Principle of the Proposed Method

2.1. The Flowchart of the Proposed Method

The flowchart of the method proposed in this study is shown as Figure 1. This study used band pass filter, SF, GA, ACO, and the Mahalanobis distance for fault diagnosis through statistical feature extracted from vibration signals of normal and abnormal states. The vibration signals are acquired by accelerometers set on the examined part of the rotating machine.

Figure 1

Flowchart of the study.

Fault characteristics of structural failures, such as unbalance, looseness, and shaft misalignment, mostly focus on low and middle frequency areas. The characteristic frequency area of cavitation is more widely distributed than that of other faults. Resulting signals of failure brought about by self-excited and impact vibration can be found in high frequency area, such as bearing flaw and shock and rubbing. For this reason, this paper introduced a feature extraction method called multiband analysis.

Feature extraction was firstly employed by band pass filter. The signals induced by four different states were, respectively, divided into three frequency ranges from low to high, as 0 Hz~1000 Hz, 1000 Hz~5000 Hz, and 5000 Hz~25000 Hz.

Then the optimum symptom parameters in different frequency bands, by which the states of rotating machinery can be sensitively distinguished, are automatically and subsequently selected by ACO based on the Mahalanobis distance between different states.

Finally, the proposed method was used to diagnose the condition of a centrifugal pump for verifying its effectiveness.

2.2. Principle of Spectrum Statistics Filter (SF) Method

SF is based on hypothesis testing in the frequency domain to eliminate the identical component between the reference signal and the primary signal. The basic principle is originated from significant difference testing based on statistics [18, 19]. A theoretical flow chart of the statistical filter method using GA is shown in Figure 2.

Figure 2

The theoretical flow chart of SF.

The primary signal $s(t)$ measured for diagnosis is contaminated by noise. The reference signal $n(t)$ is measured in normal state beforehand. Firstly, the spectra $n(f)$ and $s(f)$ are obtained by fast Fourier transformation. Then, the average ( ${\mu }_1(f)$ , ${\mu }_2(f)$ ) and variance ( ${\sigma }_1^2(f)$ , ${\sigma }_2^2(f)$ ) are going through statistical analysis. The two null hypotheses are set as follows:

$$\begin{array}{c}\text{Hypothesis}\hspace{0.17em}1:\hspace{0.17em}{\mu }_1\left(f\right)={\mu }_2\left(f\right),\\ \text{Hypothesis}\hspace{0.17em}2:{\sigma }_1^2\left(f\right)={\sigma }_2^2\left(f\right).\end{array}$$

(1)

Hypothesis testing based on statistics states with the level of significance $\alpha$ the following: if the assumption is logically impossible, the original hypothesis is denied. Therefore, the alternative hypotheses are ${\mu }_1(f)\ne {\mu }_2(f)$ and ${\sigma }_1^2(f)\ne {\sigma }_2^2\left(f\right)$ . If both null hypotheses would prove to be received, this means that at the frequency $f$ the spectrum component of the primary signal $s(t)$ would be similar to that of the reference signal $n(t)$ . The component at the frequency $f$ does not contain fault information and will be cancelled. Then the spectrum $s(f)$ of the primary signal $s(t)$ can be divided into estimated fault signal $s^{*}(f)$ and estimated noise signal $n^{*}(f)$ . Similarity factor $I_{pq}$ is defined to evaluate the similar degree between the estimated noise signal $n^{*}(t)$ and the reference signal $n(t)$ . That is to say, if the estimated noise signal $n^{*}(f)$ is optimally extracted, the estimated fault signal $s^{*}(f)$ will be also obtained. The optimum extraction of the estimated fault signal $s^{*}(f)$ is decided by the optimum level of significance $\alpha$ . The optimum $\alpha$ will be automatically decided by GA. The estimated noise signal $n^{*}(t)$ in the time domain is obtained by inverse fast Fourier transformation of the estimated noise signal $n^{*}(f)$ .

Figure 3 shows the preparatory process of the statistics filter using GA. The statistics filter has been explained in theory, and it will be described in more detail during signal analysis. An example would be the case of an unbalanced signal. Measured normal signal and unbalanced signal can be divided into $M$ parts. Spectrum analysis is performed for each part, which generates a set of spectra. The probability distribution at each frequency of normal signal and unbalanced signal is obtained from the set, as shown in Figure 2, where $F_{\text{nor}}(f_1)$ and $F_{\text{nor}}(f_2)$ are the probability distribution of the normal signal at the frequencies $f_1$ and $f_2$ . $F_{\text{unb}}(f_1)$ and $F_{\text{unb}}(f_2)$ are the probability distribution of the unbalanced signal at the frequencies $f_1$ and $f_2$ .

Figure 3

The preparatory process of statistics filter.

In order to determine whether there is significant difference between the spectrums of normal and diagnosis state at the frequency $f_i$ , hypothesis testing shown in formulae (1) is applied. The specific process of hypothesis testing is shown as Figure 4. If according to the result of hypothesis testing there is significant difference at the frequency $f_1$ between normal and diagnosis states, the spectrum at this frequency will remain. On the contrary, the spectrum at this frequency will be removed. The estimated fault spectrum $s^{*}(f)$ after statistics filter is shown in Figure 4. Moreover, the filtered spectrums $s^{*}(f)$ will be reconstructed into time domain by using inverse fast Fourier transformation (IFFT) to obtain the estimated fault signal $s^{*}(t)$ .

Figure 4

The specific process of hypothesis testing.

In the study, the coefficient optimization model is based on the basic theory of GA. A high value of the objective function is an indicator of the system's success. Then the optimal value of significance level α is automatically found.

Population size of GA is $m$ , and gene length is $n$ . The genes’ list is randomly obtained at first as shown in Table 1. The genes are encoded as formula (2). Here, $\text{gene}(k,i)$ is encoded in $i$ th codon of $k$ th gene. Then the steps, namely, mutation and crossover, are carried out. In addition, gene iteration times are 3, mutation rate is 0.5, and crossover rate is 0.3. Consider

$$\begin{array}{cc}{\alpha }_k=0.01\cdot \left({\displaystyle \sum _{i=i}^n\text{gene}\left(K,i\right)\cdot 2^{n-i}}\right),& k=1\hspace{0.17em}to\hspace{0.17em}m\end{array}.$$

(2)

Table 1

Genes’ list and the encoding process.

$\text{gene}(k,i)$	$i=\mathrm{1}$	$i=\mathrm{2}$	$i=\mathrm{3}$	$i=\mathrm{4}$	$i=\mathrm{5}$	$i=\mathrm{6}$	${\alpha }_{\mathrm{k}}$
$k=1$	1	0	1	1	1	1	${\mathrm{\alpha }}_{\mathrm{1}}=\mathrm{0.47}$
$k=2$	1	1	0	1	1	0	${\mathrm{\alpha }}_{\mathrm{2}}=\mathrm{0.54}$
$k=3$	0	1	1	1	1	0	${\mathrm{\alpha }}_{\mathrm{3}}=\mathrm{0.3}$
$k=4$	1	1	0	0	0	0	${\mathrm{\alpha }}_{\mathrm{4}}=\mathrm{0.48}$
$k=5$	0	0	0	1	1	0	${\mathrm{\alpha }}_{\mathrm{5}}=\mathrm{0.06}$
$k=6$	0	1	1	1	0	1	${\mathrm{\alpha }}_{\mathrm{6}}=\mathrm{0.29}$

Similarity factor $I_{pq}$ is defined as the evaluation factor to represent the similar degree of extracted noise $n^{*}(t)$ and reference signal $n(t)$ . The greater the similarity existing between them, the higher the accuracy of the statistical test. This paper takes the similarity factor $I_{pq}$ as the objective function to evaluate the similar degree between the estimated noise signal $n^{*}(t)$ and the reference signal $n(t)$ and sets up the GA model to optimize the significance level $\alpha$ . The similarity factor $I_{pq}$ is defined as follows:

$$I_{pq}={\displaystyle \sum _{i=1}^K\left|\log \left(\frac{q_i}{q_i^{\text{'}}}\right)\right|}\cdot K^{-1}+{\displaystyle \sum _{i=1}^K\left|\log \left(\frac{p_i}{p_i^{\text{'}}}\right)\right|\cdot K^{-1}},$$

(3) where $p_i=(\sqrt{\mathrm{2}\pi }·q_i)/{\sigma }_v$ , $q_i=({\sigma }_v/\mathrm{2}\pi {\sigma }_x)\exp (-n_i^2/\mathrm{2}{\sigma }_x^2)$ , ${\sigma }_x^2={\int }_{\mathrm{0}}^{\infty }{F}_n^2(f_k)d{f}_k$ , and ${\sigma }_v^2=(\mathrm{2}\pi f_k{)}^{\mathrm{3}}/\mathrm{3}+{\sigma }_x^2$ .

$I_{pq}$ is often used as the factor which reflects the statistical similarity of the data set. The vibration signalis divided into $K$ sections along the longitudinal axis. $n_i$ is the boundary value of each section. $p_i$ is the probability density distribution of the vibration signal of the $i$ th section. $q_i$ is the frequency density distribution of the vibration signal of the $i$ th section. And for discrete signals, $I_{pq}$ is always calculated with the approximation equations shown in this paper by using ${\sigma }_x$ and ${\sigma }_v$ , where $f_k$ is half of the sampling frequency. The transformation principle from passage frequency to frequency spectrum is introduced in [20].

2.3. Principle of ACO

In the early 1990s, ACO was introduced by Dorigo and his colleagues as a novel nature-inspired metaheuristic solution for TSP problems. The TSP problem can be described as the problem of finding the minimal distance in a closed tour which visits each town once. In the following parts, we will simply explain about the basic mathematical theory on ACO in solving TSP [21].

The medium used to communicate between ants is pheromone. In the feeding process, the moving ants lay some pheromone on the ground to mark the path. When other ants encounter a previously passed trail, they will detect the pheromone and decide with high probability to follow it and thus add their own pheromone on the trail. The process is thus characterized by a positive feedback loop, where the probability with which an ant chooses a path increases with the number of ants that chose the same path in the preceding steps [22]. After arriving at the food source, ants return to the nest following the same route.

Pheromone on the previously passed route will evaporate along with time. After each round, the route may change depending on the evaporations of the pheromone and the pheromone applied by the new ants. Pheromone evaporation mechanism makes the ants ignore the poor path selected before. As a result, the algorithm converging to local optimization early can be avoided. $\rho$ is an evaporation coefficient. When the evaporation rate is set to 1, there is no pheromone evaporation, and it not easy to get convergence. But setting $\rho$ too low is prone to get a local best answer.

The intensity of pheromone on path- $ij$ at time $t+1$ is given by [21]

$${\tau }_{ij}\left(t+1\right)=\rho \cdot {\tau }_{ij}\left(t\right)+\Delta {\tau }_{ij}\left(t,t+1\right).$$

(4)

The transition probability for the $k$ th ant from town $i$ to town $j$ is defined as follows [21]:

$$p_{ji}^k={\displaystyle \sum _{l=1}^n\frac{{\left[{\tau }_{ij}\right]}^{\alpha }{\left[{\eta }_{ij}\right]}^{\beta }}{{\left[{\tau }_{il}\right]}^{\alpha }{\left[{\eta }_{il}\right]}^{\beta }}}.$$

(5)

The field of algorithms of the ant system can be divided into three categories determined by the way the trail is updated: ant-cycle, ant-quantity, and ant-density algorithms; the formulas are given by formulae (6)–(8) [21]. In ant-cycle model, ants lay their pheromone in the end of the round, but the other two models updated the pheromone after each step. Ant-cycle is widely used, and the other two models have been gradually eliminated.

$\Delta {\tau }_{ij}^k$ is the quantity per unit of length of pheromone laid on edge ( $i,j$ ) by the $k$ th ant between time $t$ and time $t+1$ ; $Q$ is a constant; $d_{ij}$ is the distance between $i$ and $j$ ; $L^k$ is the tour length of the $k$ th ant:

ANT-quantity:

$$\Delta {\tau }_{ij}^k\left(t,t+1\right)=\{\begin{array}{l}\frac{Q_1}{d_{ij}},\\ 0,\end{array}$$

(6)

ANT-density:

$$\Delta {\tau }_{ij}^k\left(t,t+1\right)=\{\begin{array}{l}{Q}_2,\\ 0,\end{array}$$

(7)

ANT-cycle:

$$\Delta {\tau }_{ij}^k\left(t,t+n\right)=\{\begin{array}{l}\frac{Q_3}{L^K},\\ 0.\end{array}$$

(8)

In this study, the evaporation coefficient $\rho$ is set as 0.15, information intensity factor $\alpha$ is set as 0.1, and inspired factor $\beta$ is set as 0.01. The ANT-cycle algorithms model is used to update the pheromone in each trail.

2.4. Mahalanobis Distance

Clustering analysis is often based on the measurement of the distance between objects. The Euclidean distance and the Mahalanobis distance are the most commonly used distances. The Euclidean distance metric assumes that each feature of data point is equally important. On the other hand, the Mahalanobis distance takes into account the correlation of the data [23]. It is calculated using the inverse of the variance-covariance matrix of the reference data. These will be illustrated with a simple example in two dimensions. $X_i=(x_{i1},x_{i2})$ and $Y_i=(y_{i1},y_{i2})$ are data sets in two dimensions. The data set $X$ is the reference data set that contains useful information. The Euclidian distance towards the center of the data set $X$ can be calculated for each object in the set $Y$ as

$$\begin{array}{cc}E{D}_i=\sqrt{{\left(y_{i1}-{\bar{x}}_1\right)}^2+{\left(y_{i2}-{\bar{x}}_2\right)}^2,}& \text{for}\end{array}\hspace{0.17em}i=1\hspace{0.17em}to\hspace{0.17em}n,$$

(9) where $\overline{x_{\mathrm{1}}}$ and $\overline{x_{\mathrm{2}}}$ are the mean of variables $x_{\mathrm{1}}$ and $x_{\mathrm{2}}$ . The Mahalanobis distance can be calculated as

$$\begin{array}{cc}M{D}_i=\sqrt{\left(Y_i-\bar{X}\right)C_x^{-1}{\left(Y_i-\bar{X}\right)}^{\text{'}},}& \text{for}\end{array}\hspace{0.17em}i=1\hspace{0.17em}to\hspace{0.17em}n,$$

(10) where the covariance matrix of the data set $X$ is equal to $C_x$ and $C_x^{-\mathrm{1}}$ is the inverse matrix of it. However, the computation of the matrix $C_x^{-\mathrm{1}}$ can cause problems. When the data set matrix is a singular matrix, the covariance matrix cannot be inverted. In this paper, when the matrix is singular, the Mahalanobis distance of this object is set as a small value. In the symptom parameter selection process, the Mahalanobis distance is used as the objective function. The symptom parameters combination with the maximum Mahalanobis distance will be selected.

3. Brief Review of the Experimental Method

3.1. Experiment Device

Rotating machinery can be found widespread in industrial applications, which operates by means of bearings and other rotating parts [24–26]. It is well known that faults of rotating parts would cause invalidations and breakdowns of a machine. Consequently, the detection of their defects at early stage turns to be increasingly important. Centrifugal pump, one kind of rotating machinery, is used to transport fluids by the conversion energy from kinetic energy of a motor to hydrodynamic energy of the fluid flow. It is not only used as water upgrade machine but also used in food industry, plants, agriculture, oil and gas industry, paper and pulp industry, and so on. The key elements in centrifugal pump are bearing, seal, and impeller. These are exactly rotating parts and directly affect the working performance of the centrifugal pump. On another hand, the faults that are presented in pump systems are usually cavitation, leakage, impeller unbalance, and shaft misalignment. These would cause a series of troubles such as abnormal noise, high vibration, deterioration of the hydraulic performance, and pitting or erosion [27–29].

In this work, a centrifugal pump system for condition diagnosis is shown in Figure 5(a). The flow rate of pump can also be adjusted through the valve control system. The series of the motor is SF-JRO. It is employed to drive the pump through a coupling, and the rotating speed can be varied through control panel. The experimental pump made by the HONDA Company is type HAS, having output of 3.7 kW and capacity of 7.5 m³/h [30]. Five accelerometers are used to measure vibration signals for fault detection as shown in Figure 5(b). Two sensors are placed at the pump inlet, two sensors are placed at the pump outlet and the other one is placed at the pump housing. The radial direction signals of inlet are used to diagnose fault in this paper.

Figure 5

The experiment device: (a) the centrifugal pump system and (b) the position of the sensors.

The signals are measured at the rotation speed of 3000 rpm (rotating frequency is 50 Hz). The sampling frequency is 50 KHz, and the sampling time is about 10 seconds. Defects have been artificially introduced. The types of measured signals are the normal state and three typical faults, namely, cavitation, impeller unbalance, and shaft misalignment.

At the beginning of the experiment, the valves in the suction and delivery lines were fully opened. The valve in suction line was set at position 425. The flow in delivery line was 19 m³/h. In this condition, pump was in the best operating situation, which can be regarded as normal state. Then, the opening valve in the suction line was made to keep turning down slowly as a few bubbles appeared. In this status, the condition of the pump turned into cavitation state. With the valve in suction line having been turned down, the degree of cavitation was gradually aggravated. When the valve in suction line was set at position 375, cavitation signal used in this paper was measured. In the case of impeller unbalance, the damage areas caused by the impeller are 100 mm². The shaft misalignment signal has the deviation of 1.0 mm.

3.2. Symptom Parameters Used in the Paper

This paper used 13 symptom parameters in time domain and frequency domain. When time-domain symptom parameters and frequency-domain symptom parameters are used together, they must be nondimension parameters.

Five of them are time-domain symptom parameters, shown as formulae (11), and the remaining eight are frequency-domain symptom parameters as formulae (12) and (13). $P_{f6}$ to $P_{f11}$ are defined in paper [10] and are used to diagnose bearing fault. They are sensitive to feature information in high frequency band. Therefore they are more suitable for spot fault of a bearing, such as inner-race flaw, outer-race flaw, and rolling element flaw. $P_{f12}$ and $P_{f13}$ are proposed in this paper for structure failure such as unbalance and misalignment, whose fault feature information is mostly focused on low frequency band:

$$\begin{array}{cc}{P}_{t1}& =K_v=\frac{(1/N){\displaystyle \sum }_{i=1}^N{x}_i^4}{x_{\text{rms}}^4}(\text{Kurtosis}\hspace{0.17em}\text{vlaue})\\ P_{t2}& =C=\frac{x_p}{{\text{x}}_{\text{rms}}}(\text{Crest}\hspace{0.17em}\text{factor})\\ P_{t3}& =\frac{x_p}{x^{\text{'}}}(\text{Impulse}\hspace{0.17em}\text{index})\\ P_{t4}& =L=\frac{x_p}{x_r}(\text{Clearance factor}).\\ P_{t5}& =k=\frac{x_{\text{rms}}}{x\text{'}}(\text{Shape}\hspace{0.17em}\text{Factor}).\end{array}$$

(11) $x_i$ is the value of the vibration signal of point $i$ . $x_p$ is peak value, which is commonly used feature time-domain symptom parameters. $x_{\text{rms}}$ is effective value. $x_r$ is root-mean-square value. $x^{\prime }$ is the average absolute value. Consider

$$P_{f6}={\left(\frac{{\displaystyle {\sum }_{i=1}^N{f}_i^{2}S\left(f_i\right)}}{{\displaystyle {\sum }_{i=1}^{N}S\left(f_i\right)}}\right)}^{1/2}$$

(Average characteristic frequency factor)

$$P_{f7}={\left(\frac{{\displaystyle {\sum }_{i=1}^N{f}_i^4\cdot S\left(f_i\right)}}{\left({\displaystyle {\sum }_{i=1}^N{f}_i^2\cdot S\left(f_i\right)}\right)}\right)}^{1/2}$$

(Average frequency in unit time)

$$P_{f8}=\frac{{\displaystyle {\sum }_{i=1}^N{f}_i^2\cdot S\left(f_i\right)}}{{\left({\displaystyle {\sum }_{i=1}^{N}S\left(f_i\right){\displaystyle {\sum }_{i=1}^N{f}_i^4\cdot S\left(f_i\right)}}\right)}^{1/2}}$$

(Waveform stability index)

$$P_{f9}=\frac{{\displaystyle {\sum }_{i=1}^N{\left(\left|f_i-\bar{f}\right|\right)}^{1/2}}\cdot S\left(f_i\right)}{\left(\sqrt{\sigma }\cdot N\right)}$$

(Crest factor of frequency spectrum)

$$P_{f10}=\frac{\sigma }{\bar{f}}$$

(Change rate of frequency spectrum)

$$\begin{array}{c}{P}_{f11}=\frac{{\displaystyle {\sum }_{i=1}^N{\left(f_i-\bar{f}\right)}^4}\cdot S\left(f_i\right)}{{\sigma }^4\cdot N}\\ \text{(Stridency of frequency spectrum) }.\end{array}$$

(12)

$f_i$ is the frequency value of point $i$ . $S(f_i)$ is the spectrum value in frequency $f_i$ , and $N$ is the total number of the points in the spectrum; $\bar{f}$ and $\sigma$ are defined as follows: $\bar{f}={\sum }_{i=\mathrm{1}}^N{f}_i·S\left(f_i\right)/{\sum }_{i=\mathrm{1}}^{N}S\left(f_i\right)$ , $\sigma =\sqrt{{\sum }_{i=\mathrm{1}}^N{\left(f_i-\bar{f}\right)}^{\mathrm{2}}·S\left(f_i\right)/(N-\mathrm{1})}$ ,

$$\begin{array}{l}{P}_{f12}=\sqrt[N]{{\displaystyle \prod _{i=1}^{N}S\left(f_i\right)}}\\ \text{(Geometric average of frequency spectrum)}\\ P_{f13}=\sqrt[M]{{\displaystyle \prod _{j=1}^{M}S\left(f_j\right)}\hspace{0.17em}}\hspace{0.17em}\text{(Geometric average of harmonics) }.\end{array}$$

(13)

$f_j$ is the frequency value around rotating frequency or its harmonic. The range of $M$ is the exact frequency value ±5 Hz. $S(f_i)$ is the spectrum value in frequency $f_j$ , and total number of harmonic is set as 5. These two symptom parameters are aiming for unbalance and misalignment fault diagnosis. Fault characters always show as having higher peak around rotating frequency and its harmonics in spectrum.

4. Experiments Process

In the experiments, vibration signals in four states are measured, which are normal (N), impellers unbalance (U), shaft misalignment (M), and cavitation (C). The spectrums of raw signals are shown as Figure 6. Feature components for diagnosing each state are not obvious enough in these spectra.

Figure 6

The spectrums of signals before processing. (a) Normal, (b) unbalance (c) misalignment, and (d) cavitation.

4.1. Band Pass Filter

It has shown that the feature frequencies of unbalance and misalignment are mostly focused on low and middle frequency bands. And the feature frequencies of cavitation may distribute widely in all frequency bands. Therefore, band pass filters of three frequency ranges were processed, respectively. The three frequency bands are 0 Hz~1000 Hz, 1000 Hz~5000 Hz, and 5000 Hz~25000 Hz. Spectra of the signals in different bands are shown as Figures 7, 8, and 9. After band pass filter was applied, the feature information obtained from the spectra was also insufficient and unclear.

Figure 7

The spectrum of signals in the low frequency band after band pass filter. (a) Normal, (b) unbalance, (c) misalignment, and (d) cavitation.

Figure 8

The spectrum of signals in the middle frequency band after band pass filter. (a) Normal, (b) unbalance (c), misalignment, and (d) cavitation.

Figure 9

The spectrum of signals in the high frequency band after band pass filter. (a) Normal, (b) unbalance, (c) misalignment, and (d) cavitation.

4.2. Statistic Filter (SF)

SF method requires lots of data. The larger the amount of data used, the better the performance of the filtering. But the more the data that has to be processed, the more the computing time that will be spent. This is one of the difficulties of the experiments. In this study, the computing time and filtering effect have been considered. Each part of the statistics filter data contains about 65536 points. And the total number of the parts is 6.

The coefficient optimization model is based on the basic theory of GA as shown in Section 2, with the population size 6, gene length 6, an iteration of 3, mutation rate 0.5, crossover rate 0.3, and random crossover point. High fitness is the indicator of the success of the system. The optimal value of significance level α is found as Table 2 in each state and frequency band.

Table 2

Significance level $\alpha$ .

	Unb	Mis	Cav
0 Hz–1000 Hz	0.21	0.04	0.09
1000 Hz–5000 Hz	0.06	0.02	0.17
5000 Hz–25000 Hz	0.04	0.17	0.04

13 nondimension symptom parameters are calculated after SF, as formulae (11) to (13). $P_{t1}$ and $P_{t5}$ are time-domain symptom parameters; $P_{f6}$ to $P_{f11}$ are frequency-domain symptom parameters. $P_{f12}$ and $P_{f13}$ are specially proposed for unbalance and misalignment, whose fault information is more focused on low frequency band and only applied on the first frequency band.

4.3. Decision of Optimum Symptom Parameters and Frequency Bands Using ACO

The failure types that have to be identified are multiple, and it is hard to find one group parameter that can identify all of those types at once. However, the symptom parameters for identification of two types are relatively easy to find. Therefore, the automatic and sequential selection method of optimum symptom parameters using the ACO and Mahalanobis distance for sequential diagnosis is proposed in this paper. In the sequential diagnosis, we should distinguish only two states in one step [8].

The selecting model of optimum symptom parameters was inspired by the ACO system that used to solve the TSP. The current model is suitable under the condition that the number of selected parameters is known. But it can be improved to a variable number condition in a follow-up study. The symptom parameters selection course is considered as travel between cities. The number of the directed line segments is considered as the number of cities. And each directed line segment has its pheromone matrix whose dimension is the number of all symptom parameters. The pheromone matrix is set as a tiny number during initialization. Firstly, the machine ants randomly select symptom parameters out of the list that included all the parameters in this frequency band. Then the pheromone matrix is updated according to the Mahalanobis distance. And the symptom parameters with high pheromone are selected with larger probability during the following process. Finally, the sensitive symptom parameters are automatically selected by ACO.

The results of selection in different frequency band are shown in Table 3. The proportional comparison figure of Mahalanobis distance between each frequency band is shown as Figure 10. The results are the same as expected beforehand. The defect information of unbalance and misalignment is mostly focused on low frequency, and cavitation leads to a relatively wide distribution in the frequency domain. $P_{f12}$ and $P_{f13}$ are showing good performance for distinguishing unbalance and misalignment.

Table 3

Selection results of different condition in different frequency band.

Distinguish states	Frequency band	Selected symptom parameters (number)	Mahalanobis distance
N and U	0 Hz–1000 Hz (optimum)	4, 5, 12, 13	$1.27E+03$
	1000 Hz–5000 Hz	4, 7, 8, 12	$2.80E+02$
	5000 Hz–25000 Hz	6, 7, 11, 12	$4.99E+02$
N and M	0 Hz–1000 Hz (optimum)	4, 5, 12, 13	$5.99E+02$
	1000 Hz–5000 Hz	2, 3, 4, 12	$5.26E+01$
	5000 Hz–25000 Hz	2, 3, 4, 9	$1.35E+01$
N and C	0 Hz–1000 Hz	4, 5, 12, 13	$1.69E+03$
	1000 Hz–5000 Hz	1, 7, 8, 10	$1.34E+02$
	5000 Hz–25000 Hz (optimum)	3, 6, 7, 8	$2.79E+03$
U and M	0 Hz–1000 Hz	2, 7, 8, 10	$6.86E+01$
	1000 Hz–5000 Hz (optimum)	1, 2, 6, 8	$1.46E+03$
	5000 Hz–25000 Hz	5, 6, 7, 8	$1.20E+02$
U and C	0 Hz–1000 Hz	2, 7, 8, 10	$6.00E+01$
	1000 Hz–5000 Hz (optimum)	8, 9, 10, 11	$1.93E+02$
	5000 Hz–25000 Hz	6, 7, 8, 9	$1.92E+02$
M and C	0 Hz–1000 Hz	7, 9, 10, 11	$1.40E+02$
	1000 Hz–5000 Hz	6, 9, 10, 12	$5.12E+02$
	5000 Hz–25000 Hz (optimum)	4, 7, 8, 11	$4.45E+03$

Figure 10

Mahalanobis distance (MD) proportional compared in different condition and different frequency band.

4.4. Results and Discussion of Sequential Diagnosis

As shown in Figure 11, the first step of the sequential diagnosis can be used to distinguish the normal state (N) from impellers unbalance (U), shaft misalignment (M), and cavitation (C) using the sensitive symptom parameters selected by ACO. Normal state and unbalance are identified by symptom parameters [4, 5, 12, 13] in the low frequency band. Normal state and misalignment are also coincidentally identified by symptom parameters [4, 5, 12, 13] in low frequency band. Normal operation and cavitation are identified by symptom parameters [3, 6–8] in high frequency band as shown in Table 2.

Figure 11

Flowchart of sequential diagnosis for pump in the paper.

According to the result of first step, the second step can be used to distinguish the impellers unbalance (U) from shaft misalignment (M) and cavitation (C). Then the last step is used to distinguish shaft misalignment (M) and cavitation (C). This is the most difficult condition; if the state is identified as [U, M], [U, C], or [M, C] in the first step, then the second step is taken to judge the exact one of the two states which confirmed in the first step. The flowchart of sequential diagnosis in the paper is shown in Figure 11.

The practicality of the proposed method is verified by applying it to testing signals. The tests have shown that the algorithm is capable of detecting four different states in the centrifugal pump. The testing signal of cavitation is taken for an example. Eight vibration signal samples are used for the diagnosis. The type of testing signal will be identified as the state that the majority of sample signals belong to. In the first step, the diagnosis result is N, M, and C, as shown in Table 4. According to Figure 11, in this case, the states are identified as M and C. And the second step distinguishes between M and C. The diagnosis result of the second step is shown in Table 5. Finally, the state is identified as C.

Table 4

First step of sequential diagnosis.

Number	States identified based on Mahalanobis distance
	Between N and U		Between N and M		Between N and C
1	96.27	109.70	1101.65	24.68	3180.53	5.18
2	108.66	90.59	3753.48	16.98	2616.61	2.16
3	90.73	115.90	1638.37	21.56	2863.59	1.49
4	102.46	103.74	4005.46	19.57	2558.91	2.66
5	86.33	117.30	1658.46	20.47	2911.03	2.66
6	114.92	92.39	3993.00	13.15	2867.03	4.15
7	99.01	103.89	1425.83	5.83	2790.35	0.73
8	87.20	117.82	1876.30	22.81	2918.85	3.63
Diagnose by majority diagnostic method	Identified as N 6 times		Identified as M 8 times		Identified as C 8 times
	N		M		C

Table 5

Second step of sequential diagnosis.

Number	States identified based on Mahalanobis distance
	Between M and C
1	507.84	1.55
2	537.01	1.78
3	508.65	1.15
4	512.47	1.88
5	496.12	1.73
6	527.30	1.96
7	526.35	0.91
8	503.17	3.29
Diagnose by majority diagnostic method	Identified as C 8 times
	C

5. Conclusions

This paper presents a method combining statistical filter (SF) with GA, the ant colony optimization (ACO), and Mahalanobis distance for automatic decision of optimum symptom parameters and frequency bands, which can be effectively used in machinery diagnosis. In order to prove the effectiveness of the methods, it is applied to diagnose the condition of a centrifugal pump system. The considered fault types are cavitation, impeller unbalance, and shaft misalignment, which often occur in pump system. Because the feature information of different fault types focuses on different frequency bands, the band pass filter is processed to extract fault information roughly. The three bands are 0–1000 Hz as low frequency, 1000–5000 Hz as middle frequency, and above 5000 Hz as high frequency. The noise of vibration signals is optimally canceled by using the statistic filter (SF). In this step the similarity factor $I_{pq}$ is defined to evaluate the effectiveness of the filter and taken as the object function of GA to optimize its significance level. After SF application, 13 nondimensional symptom parameters are used to reflect the characteristic of the signals, 5 symptom parameters in the time domain and 8 symptom parameters in the frequency domain. Then, the selecting model of symptom parameters using ACO, which evolved from the TSP, is applied to automatically decide the optimum symptom parameters and frequency bands. Mahalanobis distance between different states is used as the objective function during the procedure of symptom parameters selection. Finally, the sequential diagnostic method is employed to identify the type of the machine state. The verification results show that ACO had good performance in feature selection, and the proposed signal analysis method in different frequency band is proven effective for centrifugal pump fault diagnosis. The sequential diagnosis method also shows a good performance as expected beforehand.