**1. Introduction**

With the vigorous development of mineral pipeline transportation technology, highpressure diaphragm pump operation, maintenance, and fault monitoring have become a concern. The high-pressure diaphragm pump provides the power required in the process of mineral transportation [1], and most of its failures are caused by one of the core components of the pump check valve. The check valve is a directional element used to control the feeding and discharging, so that the conveying medium flows in one direction and cannot flow back. Generally, the average price of each diaphragm pump exceeds CNY20 million, and the daily delivery of pulp exceeds 30,000 tons. If the stroke coefficient of the high-pressure diaphragm pump is 50 R/min, the feed and discharge check valves in a normal operating day need to act in a reciprocating manner 72,000 times. Therefore, the check valve is the most frequent and fault-prone component in a high-pressure diaphragm pump. Local faults and unplanned shutdowns in its operation state easily cause equipment pressure and flow fluctuations, resulting in equipment vibration and damage to the mineral transmission pipeline. To sum up, monitoring the operation state of the check valve is very important for the safe, stable, and efficient operation of mineral pipeline transportation. By detecting the vibration signal of the check valve and extracting the characteristics that can characterize the operation state for diagnosis and identification, the engineers and technicians can learn the operation state of the equipment in time and significantly improve the ability of operation and maintenance management. However, in the actual detection process, due to the complex field background noise, there is often strong impulse interference and random noise, making it difficult to distinguish the useful signal and noise interference, affecting the effective extraction of fault features.

Traditional signal processing methods are mostly based on Fourier transform, such as short-time Fourier transform, Wigner–Ville, power spectrum analysis, etc. [2–4]. Although

**Citation:** Yang, J.; Zhou, C. A Fault Feature Extraction Method Based on LMD and Wavelet Packet Denoising. *Coatings* **2022**, *12*, 156. https:// doi.org/ 10.3390/coatings12020156

Academic Editors: Ke Feng, Jinde Zheng and Qing Ni

Received: 23 December 2021 Accepted: 24 January 2022 Published: 27 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the above methods have been popularized in many fields of noise reduction, there are still some defects: they cannot decompose the signal adaptively. The vibration signal of the check valve is a typical non-stationary and nonlinear signal. Before processing and analyzing it, it is necessary to obtain the local properties of the signal in the time domain. Through the ongoing efforts of scholars in related fields, some breakthroughs have been made. The empirical mode decomposition (EMD) proposed by Huang et al. [5] can adaptively decompose the signal into a series of intrinsic mode functions (IMF). However, this method has problems such as the endpoint effect, mode aliasing, etc. In recent years, the local mean decomposition (LMD) [6] proposed by Smith is a new time-frequency analysis method, which can adaptively decompose the signal to be processed into several product functions (PF) with high to low frequency and has been gradually applied. Compared with traditional methods, LMD has the advantages of better adaptability and time-frequency aggregation, solves the problems of over the envelope and under the envelope of EMD, has fewer iterations, and effectively improves the accuracy of signal decomposition. Dong Linlu et al. [7] combined the respective advantages of LMD and the singular-value decomposition method to denoise the noisy micro-vibration signal. The results show that the proposed method can better remove the high-frequency noise of the micro-vibration signal, which lays a foundation for further analysis of vibration signals. To solve the influence of complex noise and measuring point position in the monitoring process, Wang Haijun et al. [8] proposed a processing method combining a data fusion algorithm and LMD. The results show that the combined method is better than the single digital filtering method, and the vibration monitoring of a hydropower plant is realized better. Wang Zhijian et al. [9] introduced the mask signal method to process the component signals screened after LMD decomposition, which effectively suppressed the noise and weakened the mode aliasing phenomenon in the process of signal decomposition. Through the signal's LMD decomposition, the complete time-frequency distribution information can be obtained, and it is convenient to extract the signal characteristics. At the same time, the latest applications of the LMD method in other fields have been reported. Gupta P et al. [10] proposed the method of combining LMD and an artificial neural network to complete the tool chatter feature extraction in the turning process. From the experimental analysis, it can be seen that the obtained safe cutting zone is important and can limit chatter. Liao L et al. [11] completed the extraction of low-frequency noise under the background of aerodynamic noise by using the improved LMD method. Through the wind turbine experiment, the results show that the signal noise extraction effect is good, and the calculation efficiency of the algorithm is good. To solve the problematic identification of welding quality, Huang y et al. [12] proposed a classification algorithm combining LMD and the deep belief network (DBN). Experimental results show that this method has a higher recognition rate than traditional classification methods. Lee CY et al. [13] introduced LMD, wavelet packet decomposition (WPD), and other ways to solve the problems existing in the rotor fault diagnosis of an induction motor. The results show that the model can reduce the dimension of the original data and has high robustness. However, due to noise interference, the extracted PF component will be distorted to a certain extent, which affects the interpretation of the useful components of the signal. Using wavelet packet to process the fault signal can effectively overcome the problem of low-frequency resolution in the high-frequency band of the traditional wavelet transform [14] and accurately describe the signal information, which is conducive to the effective elimination of noise information.

However, due to the bad working environment and changeable working conditions of the high-pressure diaphragm pump, some noise will still remain in the signal component decomposed by the time-frequency analysis method. If fault feature extraction and analysis are carried out directly, it will have a certain impact on fault diagnosis results. Therefore, it is also necessary to reduce the noise of the signal. Wavelet packet transform is based on wavelet transform. Compared with wavelet transform, wavelet packet transform has higher-frequency resolution in high-frequency and low-frequency bands and has good adaptability. In recent years, wavelet transform has been widely used in the noise reduction

of non-stationary signals. Kuanfang he et al. [15] proposed a wavelet packet denoising method for welding acoustic emission signals. The experimental results show that the proposed method can effectively process the acoustic emission signals of welding cracks. Wang X et al. [16] proposed a method based on optimized variational mode decomposition (VMD) and wavelet packet threshold denoising and applied it to remove strong white noise signals. The results show that this method retains the practical components of the signal well. Sun W et al. [17] used a wavelet packet to remove the noise in the collected bearing signal and then combined it with the LMD method to extract the fault feature. The results show that this method can effectively extract fault features.

Based on the above analysis, a fault feature extraction method based on local mean decomposition (LMD) and wavelet packet denoising is proposed in this paper. Firstly, the vibration signal of the check valve was decomposed by LMD, and the K-L divergence of each AM-FM component was obtained. The practical component signal was selected for reconstruction through the comparative analysis with the set threshold. Then, wavelet packet transform was used to denoise the reconstructed signal to reduce the noise interference. Finally, Hilbert envelope spectrum analysis asper formed on the denoised signal to extract the fault signal features. The effectiveness of the proposed method was verified by analyzing the actual fault data of the check valve.

The rest of the paper is organized as follows: In Section 2, we introduce the basic principle of the LMD algorithm and K-L divergence and describe the process of PF component selection. In Sections 3 and 4, we introduce the basic principles of Hilbert envelope demodulation and wavelet packet denoising. In Section 5, we describe the modeling process of check valve fault feature extraction based on LMD and wavelet packet denoising. Section 6 is experimental analysis. The effectiveness of the method proposed in this paper was verified by the collected actual fault data of the check valve. Section 7 concludes the paper.

#### **2. LMD Algorithm Based on K-L Divergence**

## *2.1. LMD Algorithm*

LMD is a time-frequency analysis method proposed by Smith. Its principle is to decompose the signal into product function components of different frequencies and a margin, and each component is obtained by multiplying a pure FM signal and envelope signal. If the original signal is *x*(*t*), the decomposition step can be described in the following form.

(1) Calculate all local extreme points of the original signal *x*(*t*), and then calculate the average value *mi* of adjacent local extreme points *ni* and *ni+*1:

$$m\_i = \frac{n\_i + n\_{i+1}}{2} \tag{1}$$

The moving average method is used to smooth the line composed of *mi*, and the local mean function *m*11(*t*) can be obtained.

(2) The envelope estimation value *ai* is calculated according to the adjacent extreme points *ni* and *ni+*1:

$$n\_i = \frac{|n\_i - n\_{i-1}|}{2} \tag{2}$$

Similarly, the envelope estimation function *a*11(*t*) is obtained by the moving average method.

(3) The local mean function *m*11(*t*) is separated from the original signal *x*(*t*).

$$h\_{11}(t) = \mathbf{x}(t) - m\_{11}(t) \tag{3}$$

(4) Divide *h*11(*t*) by *a*11(*t*) to obtain FM signal *s*11(*t*)

$$s\_{11}(t) = h\_{11}(t) / a\_{11}(t) \tag{4}$$

When *a*12(*t*) = 1, *s*11(*t*) is a standard pure FM signal, and when *a*12(*t*) is not equal to 1, *s*11(*t*) is taken as the original data, and the above processes repeated. Until *<sup>a</sup>*1(*n*+<sup>1</sup>)(*t*) = 1 is met, the pure FM signal *<sup>s</sup>*1*n*(*t*) is calculated, as shown in the following formula:

$$\begin{cases} \begin{array}{l} h\_{11}(t) = \mathbf{x}(t) - m\_{11}(t) \\ h\_{12}(t) = \mathbf{x}(t) - m\_{12}(t) \\ \vdots \\ h\_{1n}(t) = \mathbf{x}(t) - m\_{1n}(t) \end{array} \\\\ \begin{cases} s\_{11}(t) = h\_{11}(t) / a\_{11}(t) \\ s\_{12}(t) = h\_{12}(t) / a\_{12}(t) \\ \vdots \\ s\_{1n}(t) = h\_{1n}(t) / a\_{1n}(t) \end{cases} \end{cases} \tag{6}$$

The conditions for iteration termination are as follows:

$$\lim\_{n \to \infty} a\_{1n}(t) = 1 \tag{7}$$

(5) Calculate the envelope signal:

$$a\_1(t) = a\_{11}(t)a\_{12}(t) \cdots a\_{1n}(t) = \prod\_{q=1}^{n} a\_{1q}(t) \tag{8}$$

(6) The product function can be obtained by multiplying the pure FM signal *<sup>s</sup>*1*n*(*t*) and the envelope signal *<sup>a</sup>*1(*t*)

$$PF\_1(t) = a\_1(t)s\_{1n}(t)\tag{9}$$

(7) Separate *PF*1(*t*) from *x*(*t*) and find the signal *<sup>u</sup>*1(*t*). Then, repeat the above steps *k* times. When *uk*(*t*) becomes a monotone function, the operation is terminated.

$$\begin{cases} u\_1(t) = x(t) - PF\_1(t) \\\\ u\_2(t) = u\_1(t) - PF\_2(t) \\ \vdots \\ u\_k(t) = u\_{k-1}(t) - PF\_k(t) \end{cases} \tag{10}$$

The following formula can be obtained by combining all PF components and recombining.

$$\mathbf{x}(t) = \sum\_{p=1}^{k} PF\_p(t) + u\_k(t) \tag{11}$$

From the above description, it can be seen that the LMD is an adaptive decomposition method based on the local extreme points of the signal itself. The decomposition process is a multi-cycle iterative process.

#### *2.2. PF Component Selection Based on K-L Divergence*

Several PF components obtained by LMD have a different correlation with the original signal. To effectively screen component signals, appropriate screening methods must be adopted. K-L divergence is called directional divergence, describing the difference between the two probability distributions [18]. Let the two probability distributions be *p*1(*x*) and *p*2(*x*), then the K-L distance is defined as follows:

$$\delta(p\_1, p\_2) = \int p\_1(\mathbf{x}) \log \frac{p\_1(\mathbf{x})}{p\_2(\mathbf{x})} d\mathbf{x} \tag{12}$$

The K-L divergence values of *p*1(*x*) and *p*2(*x*) are:

$$D(p\_1, p\_2) = \delta(p\_1, p\_2) + \delta(p\_2, p\_1) \tag{13}$$

The vibration signals collected in the equipment characterize the vibration of the equipment. Therefore, assuming that the probabilities between the two signals are *p*1(*x*) and *p*2(*x*), the closer the vibration between them, the closer *p*1(*x*) and *p*2(*x*) will be, and the smaller the K-L divergence value will be.

Since several PF components are obtained through LMD, the solution of K-L divergence for the original signal *X* = {*<sup>x</sup>*1, *x*2 ... *xn*} and component *Y* = {*y*1, *y*2 ... *yn*} decomposed by the LMD algorithm is as follows:

(1) Solve the probability distribution of the above two signals and let the function *p*1(*x*) be the kernel density estimation of signal *X*:

$$p\_1(\mathbf{x}) = \frac{1}{nh} k \sum\_{i=1}^{n} \left[ \frac{\mathbf{x}\_i - \mathbf{x}}{h} \right], \mathbf{x} \in \mathcal{R} \tag{14}$$

where *k*(·) is the kernel function, which usually uses the Gaussian kernel function, which is *k*(*u*) = √12*π e*<sup>−</sup>*<sup>u</sup>*2/2, and *h* is a given positive number. Similarly, *p*2(*x*) can be obtained in the above way.

(2) The K-L distances *X* = {*<sup>x</sup>*1, *x*2 ... *xn*} and *Y* = {*y*1, *y*2 ... *yn*} are obtained by substituting the probability distributions of signal *<sup>δ</sup>*(*p*1, *p*2) and signal *<sup>δ</sup>*(*p*2, *p*1) into Formula (12).

(3) By substituting the calculated *<sup>δ</sup>*(*p*1, *p*2) and *<sup>δ</sup>*(*p*2, *p*1) into Formula (13), the corresponding K-L divergence can be obtained.

(4) Normalize all calculated K-L divergence values.

(5) Filter according to the preset threshold. If the component is less than the threshold, it will be regarded as the component containing an obvious fault signal. According to the data obtained in this study, the set threshold is 0.03.

#### **3. Hilbert Envelope Demodulation**

The Hilbert envelope demodulation method mainly converts the actual signal into an analytical signal by Hilbert transform and then takes the modulus to obtain its envelope. Assuming that the signal *x*(*t*) has amplitude modulation envelope *A*(*t*) and phase modulation function *ϕ*(*t*), its expression is as follows:

$$x(t) = A(t)\cos(2\pi ft + q(t))\tag{15}$$

The Hilbert transform of signal *x*(*t*) can be an approximately 90◦ phase shift of *x*(*t*):

$$
\hat{x}(t) = A(t)\sin(2\pi ft + \varphi(t))\tag{16}
$$

Construct the analytical signal *Z*(*t*) so that:

$$Z(t) = \mathbf{x}(t) + j\mathbf{\hat{x}} = A(t)e^{j\varphi(t)}\tag{17}$$

By calculating the modulus of the above formula, the following formula can be obtained:

$$|Z(t)| = A(t) \tag{18}$$

In the above formula, |*Z*(*t*)| is the envelope of the signal.

#### **4. Principle of Wavelet Packet Denoising**

As one of the more detailed signal analysis methods, the advantage of wavelet packet transform is that it can adaptively select the corresponding frequency band and divide the frequency band at multiple levels [19,20]. For a one-dimensional signal containing noise, it can be described as the following formula:

$$S(i) = f(i) + \mathfrak{a}\beta(i) \text{ ( $i = 0, 1, 2, \dots, n - 1$ )}\tag{19}$$

where *f*(*i*) is the actual signal, *S*(*i*) is the signal containing noise, and *β*(*i*) is noise. The process of wavelet packet denoising is as follows:

(1) Determine the corresponding wavelet base and decomposition level, and then start wavelet packet decomposition;

(2) According to the established entropy standard, the optimal tree is obtained, and the optimal wavelet basis is determined;

(3) The appropriate threshold is selected, and the high-frequency coefficients with different decomposition scales are quantified at the same time;

(4) The signal reconstruction operation is carried out according to the decomposition coefficients and quantization coefficients of the N-th layer wavelet packet.

The decomposition structure of wavelet packet is shown in Figure 1 [21].

**Figure 1.** Decomposition structure of wavelet packet.

#### **5. Fault Feature Extraction Model of Check Valve Based on Local Mean Decomposition Wavelet Packet Denoising**

Based on the above theoretical analysis, a check valve fault feature extraction method based on local mean decomposition wavelet packet denoising is proposed in this paper. The research outline of the proposed method is shown in Figure 2, and the method implementation process is shown in Figure 3. The specific steps are as follows:

**Figure 2.** Research outline of the proposed method.

**Figure 3.** Fault feature extraction process of check valve based on LMD and wavelet packet denoising.

(1) The experimental failure data of wear breakdown of check valve are selected, and then the above samples are decomposed into several PF components by LMD decomposition.

(2) According to the calculation method of K-L divergence, the K-L divergence between the original signal and each PF component obtained by decomposition is calculated. Then, all K-L divergence values are normalized.

(3) According to the set threshold, the calculated K-L divergence value is compared with it, and then any PF component less than the threshold is filtered.

(4) The filtered PF components are reconstructed, and then wavelet packet denoising is carried out.

(5) Hilbert envelope demodulation is performed on the denoised signal.
