# Fault Diagnosis of Oil Pumping Machine Retarder Based on Sound Texture-Vibration Entropy Characteristics and Gray Wolf Optimization-Support Vector Machine

Authors: Shutao Zhao; Ke Chang; Erxu Wang; Bo Li; Kedeng Wang; Qingquan Wu
Journal: Shock and Vibration
 (2020)

---

## Abstract

In order to diagnose the retarder faults of oil pumping machine accurately in complex environments and improve the generalization of the algorithm, a GWO-SVM fault diagnosis algorithm based on the combination of sound texture and vibration entropy characteristics was proposed. Firstly, the acquired sound signal was purified by band-pass filter, then generalized S-transform was developed to extract the box dimension, directivity, and contrast ratio, which reflect the characteristics of time-frequency spectrum, to construct three-dimensional texture eigenvectors. Secondly, theK parameter of variational mode decomposition (VMD) was reasonably selected by the energy method, and then the vibration signal was decomposed to get modal components, and the permutation entropy was obtained from modal components. Finally, joint eigenvectors were constructed and fed into SVM for learning. The gray wolf optimization (GWO) algorithm was used to optimize the parameters of the SVM model based on mixed kernel function, which reduces the impact of sensor frequency response, environmental noise, and load fluctuation disturbance on the accuracy of retarder fault diagnosis. The results showed that the GWO-SVM fault diagnosis method, which is based on the combination of sound texture and vibration entropy characteristics, makes full use of the complementary advantages of signal frequency band. And the overall diagnostic accuracy for the experimental samples reaches 100%, which has good generalization ability.


---

## Body

## 1. Introduction

Oil pumping machine relies on the up and down movement of the horsehead to complete the lifting of the crude oil from the wellbore. The retarder connecting the crank train is the key component of the power transmission. Due to friction and impact, oil pumping machine is prone to oil leakage, gear damage, bearing damage, belt wear, and other faults [1, 2]. Vibration signal characteristics can be used to diagnose mechanical failure of oil pumping machine retarder [3, 4]. But limited to the charge accumulating effect and coupling mode of the piezoelectric accelerometer, it is likely to cause detection failure from the continuous impact caused by the pitting and broken cog of the retarder. The reason is that the vibration sensor only detects a narrow frequency response range and is insensitive to high-frequency signal changes. The sound signal accompanying the operation of the oil pumping machine is homologous to the vibration signal and can be obtained by the noncontact electret film capacitive sensor, which can effectively compensate for the detection failure phenomenon caused by the vibration sensor band limitation.For nonstationary signals such as sound and vibration, there are mainly analysis methods such as dynamic time warping (DTW), wavelet transform (WT), empirical mode decomposition (EMD), and local mean decomposition (LMD) [5–8]. The DTW planning optimal path is prone to metamorphosis distortion, WT has energy leakage, and the two cannot adaptively decompose the signal. The EMD adaptive decomposition process is prone to over enveloping, end effect, and modal aliasing. LMD optimizes underfrequency, over envelopes, and other issues of EMD. However, LMD and EMD are essentially recursive decomposition methods, which cannot eliminate the endpoint effect, modal aliasing, and other issues fundamentally. Variational mode decomposition (VMD) is a new nonrecursive mode decomposition method, which avoids the modal aliasing problem caused by envelope error [9, 10]. The generalized S-transform is developed by short-time Fourier transform (STFT) and Wigner–Ville distribution (WVD), which makes up for the defects of STFT single resolution and WT phase and without cross-term interference problem. It introduces a frequency-dependent adjustment factor into the window function to adjust the time-frequency resolution [11].For the sound-vibration signal combination method, Zhao et al. [12] process sound and vibration signals via improving the ensemble empirical mode decomposition (EEMD) to calculate the two-dimensional spectral entropy. The diagnostic effect is improved. The literature [13] utilizes wavelet packet and feature entropy theory to extract the features of the collected sound and vibration signals, which improves the diagnostic accuracy. Zhang et al. [14] construct three-dimensional (3D) graph of sound and vibration and extract envelope of 3D graph and hierarchical eigenvectors based on shape, which provides a new idea for circuit breaker diagnosis. In [15], LMD is utilized to decompose sound and vibration signals and select appropriate PF component to obtain feature entropy as the eigenvectors. The diagnostic accuracy is improved. Although studies [12–15] have achieved some good results, the difference in sound and vibration signals is not considered. The feature extraction methods of the two signals are identical, which leads to misdiagnosis and poor generalization.It is a beneficial attempt to distinguish the defect of oil pumping machine retarder by the features combination of sound and vibration signals. The generalized S-transform is very sensitive to high-frequency impact signals and has high time-frequency resolution [16], which can accurately reflect the high-frequency impact characteristics of sound signals and rich time-frequency information. The load fluctuation of oil pumping machine’s up and down stroke causes the spectrum of the vibration signal of the retarder to be extremely complicated. The signal components include the combination of gear shaft rotation, gear meshing frequency, and so on [17]. Performing VMD method can effectively improve misdiagnosis and missed diagnosis caused by band aliasing. In this paper, the combined features of sound and vibration signals are extracted and combined with the GWO-SVM model for fault identification. The diagnosis process is shown in Figure 1.Figure 1
Diagnostic flowchart.AC144 piezoelectric acceleration sensor (0.6–10000 Hz) and NVL-AF-audio embedded waterproof (explosion-proof) high-fidelity pickup (20–20000 Hz) were acquired to collect CYJ10-3-48HB oil pumping machine retarder (CJH1100 × 73) multiple sets of sound and vibration sample signals under retarder oil leakage, gear pitting peeling, belt damage, and normal state to extract sound texture-vibration entropy characteristics, respectively.


## 2. Sound Signal Feature Extraction

The motor transmits power to the retarder through the belt. The retarder reduces the high-speed rotation of the motor through three-axis two-stage deceleration to the low-speed rotation of the crankshaft. The whole operation process is relatively complex, and it runs in the open air. As a result, the sound signals collected include the belt friction noise of oil pumping machine, the noise of the motor, wind noise, thunder noise, and human voice. Through spectrum analysis, it is not hard to find that the above noise is mostly low-frequency interference, and the frequency is concentrated between 0 and 10 kHz. So the frequency band noise of sound signals below 10 kHz and above 20 kHz is filtered by finite impulse response (FIR) band-pass filter. The sound signals in normal operation state between 10 and 20 kHz, before and after denoising, are shown in Figure2.Figure 2
Comparison of results before and after denoising.


### 2.1. Sound Signal Generalized S-Transform

The windowed Fourier transform for signalxt is performed:(1)Sf=∫−∞+∞xtwte−j2πftdt.The Gauss window function is scaledσ and translated τ, that is,(2)wτ−t=1σ2πe−τ−t2/2σ2.The time-frequency spectrum ofxt is obtained by introducing formula (2) into formula (1):(3)Sσ,τ=∫−∞+∞xt1σ2πe−τ−t2/2σ2e−j2πftdt.Setσf=1/f, get the S-transform as follows:(4)Sf,τ=∫−∞+∞f2πe−f2τ−t2/2xte−j2πftdt=∫−∞+∞ωf,τ−txte−j2πftdt,ωf,τ−t=f2πe−f2τ−t2/2,where f is the frequency, t is the time, τ is the positional parameter controlling the Gaussian window on the time axis t, and ωf,τ−t is the Gaussian window function, and the height and width of which vary with the frequency f. At the same frequency, the Gaussian window function of different signals is fixed. In order to improve the energy concentration of time-frequency, parameters αandβ are introduced. The basic structure is as follows:(5)Sf,τ=∫−∞+∞fαβ2πe−f2ατ−t2/2β2xte−j2πftdt=∫−∞+∞ωf,τ−txte−j2πftdt,where αandβ are adjustment factors and are generally positive. When β increases or α decreases, the Gaussian window width is stretched in the time domain and compressed in the frequency domain, and vice versa. The values of α and β are selected as 1 and 2, respectively. The time-frequency diagram of the generalized S-transform of the vibration signal in the four states is shown in Figure 3.Figure 3
Time-frequency diagram of generalized S-transform in four states: (a) normal state, (b) belt damage, (c) retarder oil leakage, and (d) gear pitting peeling.
(a)
(b)
(c)
(d)The time-frequency texture images, after generalized S-transform, embody different texture features under different faults. The box dimension, directivity, and contrast ratio which reflect the time-frequency texture features are extracted as eigenvectors to classify faults.


### 2.2. Box Dimension

The box dimension is particularly sensitive to the texture roughness of the image spatial distribution and can quantitatively characterize the distribution law of the time-frequency diagram. The calculation process is as follows: subblocks of sizeM × M are divided into L × L (1 < L ≤ M/2, L is an integer), make r = L/M, then each grid is a column L × L × L1 box, L1 = L × G/M, and G is gray level. Suppose the minimum and maximum values of the image grayscale in the (i, j) grid fall in the k-th and l-th boxes, respectively, then the number of boxes required for the (i, j) grid image is given by(6)nri,j=k−l+1.The number of boxes required to cover the entire image isNr:(7)Nr=∑nri,j.The box dimension can be calculated by the following formula:(8)D=limr⟶0logNrlog1/r.The logNr and log 1/r slopes are fitted by least-squares method, and the absolute value is the box dimension.


### 2.3. Directivity

Directivity describes the global characteristics of texture images and characterizes the trend of basic texture units and their arrangement tendency in all directions. Calculate the gradient vectorΔV of each pixel, whose modulus and direction are defined as follows:(9)ΔV=ΔP+ΔQ2,θ=tan−1ΔQΔP+π2,where ΔV is the gradient vector modulus; ΔP and ΔQ are the changes of ΔV in the horizontal and vertical directions, respectively; and θ is the angle of the gradient vector.Make a histogram and calculate the directivity by calculating the peak value of the gradient statistical histogram. The formula is as follows:(10)Fdir=1−r⋅np×∑p=1np∑φ∈wpϕ−ϕp2×HDϕ,where wp is the amplitude of the p-th peak, np is the number of histogram peaks, the subscript D of HD is the D-th peak of the histogram, and r is the normalization factor.


### 2.4. Contrast Ratio

The contrast ratio reflects the difference of the image in gray level. The larger the difference in gray value, the stronger the contrast ratio. When there is a significant peak in the gray value of 0 or 255, the deviation degree is measured by the kurtosisk4, and the calculation is as follows:(11)Fc=σhk4,where k4=μ4/σh4, μ4 is the 4th moment of the gray mean of the whole image and σh is the standard deviation of the image.Calculate the box dimension, directivity, and contrast ratio of the generalized S-transformed time-frequency diagram of the sound signal in the four states (belt breakage, retarder oil leakage, gear pitting peeling, and normal state), and record them asD, Fdir, and Fc, respectively. The sound texture features are shown in Table 1.Table 1
Sound texture features.
Operating state
Group
Sound texture
D
F
dir
Fc
Belt damage
1
2.5674
4.5897
1.2483
2
2.5661
4.5812
1.2429
Retarder oil leakage
1
3.8370
8.3845
2.2430
2
3.8367
8.3861
2.2423
Gear pitting peeling
1
5.4279
6.3344
3.3394
2
5.4286
6.3367
3.3323
Normal state
1
4.6391
5.4335
4.3428
2
4.6387
5.4345
4.3425The texture features of the abovementioned sound signal have antinoise and rotation invariance and can describe the local pattern and arrangement rule of the image. Longitudinal and horizontal comparison of the characteristics of the sound texture can identify the operating state of the retarder.


## 2.1. Sound Signal Generalized S-Transform

The windowed Fourier transform for signalxt is performed:(1)Sf=∫−∞+∞xtwte−j2πftdt.The Gauss window function is scaledσ and translated τ, that is,(2)wτ−t=1σ2πe−τ−t2/2σ2.The time-frequency spectrum ofxt is obtained by introducing formula (2) into formula (1):(3)Sσ,τ=∫−∞+∞xt1σ2πe−τ−t2/2σ2e−j2πftdt.Setσf=1/f, get the S-transform as follows:(4)Sf,τ=∫−∞+∞f2πe−f2τ−t2/2xte−j2πftdt=∫−∞+∞ωf,τ−txte−j2πftdt,ωf,τ−t=f2πe−f2τ−t2/2,where f is the frequency, t is the time, τ is the positional parameter controlling the Gaussian window on the time axis t, and ωf,τ−t is the Gaussian window function, and the height and width of which vary with the frequency f. At the same frequency, the Gaussian window function of different signals is fixed. In order to improve the energy concentration of time-frequency, parameters αandβ are introduced. The basic structure is as follows:(5)Sf,τ=∫−∞+∞fαβ2πe−f2ατ−t2/2β2xte−j2πftdt=∫−∞+∞ωf,τ−txte−j2πftdt,where αandβ are adjustment factors and are generally positive. When β increases or α decreases, the Gaussian window width is stretched in the time domain and compressed in the frequency domain, and vice versa. The values of α and β are selected as 1 and 2, respectively. The time-frequency diagram of the generalized S-transform of the vibration signal in the four states is shown in Figure 3.Figure 3
Time-frequency diagram of generalized S-transform in four states: (a) normal state, (b) belt damage, (c) retarder oil leakage, and (d) gear pitting peeling.
(a)
(b)
(c)
(d)The time-frequency texture images, after generalized S-transform, embody different texture features under different faults. The box dimension, directivity, and contrast ratio which reflect the time-frequency texture features are extracted as eigenvectors to classify faults.


## 2.2. Box Dimension

The box dimension is particularly sensitive to the texture roughness of the image spatial distribution and can quantitatively characterize the distribution law of the time-frequency diagram. The calculation process is as follows: subblocks of sizeM × M are divided into L × L (1 < L ≤ M/2, L is an integer), make r = L/M, then each grid is a column L × L × L1 box, L1 = L × G/M, and G is gray level. Suppose the minimum and maximum values of the image grayscale in the (i, j) grid fall in the k-th and l-th boxes, respectively, then the number of boxes required for the (i, j) grid image is given by(6)nri,j=k−l+1.The number of boxes required to cover the entire image isNr:(7)Nr=∑nri,j.The box dimension can be calculated by the following formula:(8)D=limr⟶0logNrlog1/r.The logNr and log 1/r slopes are fitted by least-squares method, and the absolute value is the box dimension.


## 2.3. Directivity

Directivity describes the global characteristics of texture images and characterizes the trend of basic texture units and their arrangement tendency in all directions. Calculate the gradient vectorΔV of each pixel, whose modulus and direction are defined as follows:(9)ΔV=ΔP+ΔQ2,θ=tan−1ΔQΔP+π2,where ΔV is the gradient vector modulus; ΔP and ΔQ are the changes of ΔV in the horizontal and vertical directions, respectively; and θ is the angle of the gradient vector.Make a histogram and calculate the directivity by calculating the peak value of the gradient statistical histogram. The formula is as follows:(10)Fdir=1−r⋅np×∑p=1np∑φ∈wpϕ−ϕp2×HDϕ,where wp is the amplitude of the p-th peak, np is the number of histogram peaks, the subscript D of HD is the D-th peak of the histogram, and r is the normalization factor.


## 2.4. Contrast Ratio

The contrast ratio reflects the difference of the image in gray level. The larger the difference in gray value, the stronger the contrast ratio. When there is a significant peak in the gray value of 0 or 255, the deviation degree is measured by the kurtosisk4, and the calculation is as follows:(11)Fc=σhk4,where k4=μ4/σh4, μ4 is the 4th moment of the gray mean of the whole image and σh is the standard deviation of the image.Calculate the box dimension, directivity, and contrast ratio of the generalized S-transformed time-frequency diagram of the sound signal in the four states (belt breakage, retarder oil leakage, gear pitting peeling, and normal state), and record them asD, Fdir, and Fc, respectively. The sound texture features are shown in Table 1.Table 1
Sound texture features.
Operating state
Group
Sound texture
D
F
dir
Fc
Belt damage
1
2.5674
4.5897
1.2483
2
2.5661
4.5812
1.2429
Retarder oil leakage
1
3.8370
8.3845
2.2430
2
3.8367
8.3861
2.2423
Gear pitting peeling
1
5.4279
6.3344
3.3394
2
5.4286
6.3367
3.3323
Normal state
1
4.6391
5.4335
4.3428
2
4.6387
5.4345
4.3425The texture features of the abovementioned sound signal have antinoise and rotation invariance and can describe the local pattern and arrangement rule of the image. Longitudinal and horizontal comparison of the characteristics of the sound texture can identify the operating state of the retarder.


## 3. Vibration Signal Feature Extraction

The energy spectrum of the vibration signal is concentrated within 10 kHz, and the vibration signals of the normal state and belt damage state of oil pumping machine are collected. The time domain waveform comparison is shown in Figure4.Figure 4
Comparison of normal and abnormal waveforms in time domain.


### 3.1. VMD Decomposition

#### 3.1.1. VMD Decomposition Principle

The VMD decomposition is mainly divided into two parts: the establishment and solution of the variational constraint problem. The following problems are solved for the vibration signal with the data lengthN in the oil production process:(12)minuk,ωk∑k∂tδt+jπt∗μkte−jωkt2s.t.∑kμk=f,where μp≔μ1,…,μp is the decomposed p modes and ωp≔ω1,…,ωp is the center frequency of the p modes.In order to solve the optimal solution, the quadratic penalty factorα and the Lagrangian operator λt are introduced to change the constrained variational problem into a nonbinding variational problem. The extended Lagrangian expression is as follows:(13)Lμk,ωk,λ=α∑k∂tδt+jπt∗μkte−jωkt22+ft−∑kμkt22+λt,ft−∑kμkt.The saddle point of equation (13) is solved by the ADMM method, so that μkn+1, ωkn+1, and λn+1 are continuously updated, and the modal component μk and center frequency ωk are solved as follows:(14)μ^kn+1ω=f^ω−∑i≠kμ^knω+λ^ω/21+2αω−ωk2,(15)ωkn+1=∫0∞ωμ^knωdω∫0∞μ^knωdω.The VMD steps are as follows:(1)
Initializeμ^k1, ωk1, λ^1, and n, let its initial value be 0, and set the decomposition modal number K to 2 (predecomposition optimization).(2)
Updateμk and wk according to equations (14) and (15), respectively.(3)
Updateλ^:(16)λ^nω+τf^ω−∑μ^kn+1⟶λ^n+1ω.(4)
If the following formula is satisfied, the iteration is stopped and the result is output; otherwise, it returns to step 2:(17)∑μ^kn+1−μ^kn22μ^kn22<ε.


#### 3.1.2. VMD ParameterK Optimization

In order to prevent the VMD from being decomposed, theK parameter is selected according to the energy conservation theory before and after decomposition. For the original retarder vibration signal sequence xi, the energy calculation formula is as follows:(18)E=∑i=1nx2in,where E represents the signal energy value and n is the sampling point. In order to characterize the energy difference before and after the VMD decomposition, the energy difference parameter ψ is defined and calculated as follows:(19)ψ=∑x=1KEx−EE×100%,where Ex corresponds to the energy of the x component, K is the number of components, and E is the original signal energy. The total energy of K components after VMD decomposition is equal to the energy of the original signal (i.e., the ideal value of ψ is 0), after many experiments and calculations, the trend of K is shown in Figure 5.Figure 5
Change trend ofK value.It can be seen from Figure5 that when K is greater than 6, the energy difference parameter λ is increased, and it can be judged that overdecomposition occurs. At this time, the K value at the turning point is the optimal decomposition mode number of the VMD. The time-frequency diagram obtained by decomposing the vibration signal in the belt damage state is shown in Figures 6 and 7.Figure 6
Time domain diagram of vibration signal VMD decomposition.Figure 7
Frequency domain diagram of vibration signal VMD decomposition.It can be seen from the frequency domain spectrum of Figure7 that the VMD decomposition of the vibration signal effectively improves the modal aliasing phenomenon and provides a strong support for accurate fault diagnosis.


### 3.2. Calculated Permutation Entropy

Permutation entropy can detect the sudden change of the signal, has strong antinoise ability and high time resolution, and is highly targeted to the nonstationary chaotic vibration signal under complex environment of oil production field.For the signal sequenceXi,i=1,2,…,n, phase space reconstruction is performed to obtain the following matrix:(20)x1x1+τ⋯x1+m−1τx2x2+τ⋯x2+m−1τxjxj+τ⋯xj+m−1τ⋯⋯xkxk+τ⋯xk+m−1τ,where j=1,2,…,k, k=n−m−1τ, and m and τ are the embedded dimension and delay time. The elements of each component (total of k) are arranged in ascending order according to the numerical value, and different symbol sequences are obtained according to the reconstructed component index. The m-dimensional phase space mapping symbol sequence has m! species, and the probability of occurrence of s different symbol sequences is P1,P2,…,Ps, respectively, calculated as follows:(21)Pem=−∑j=1sPjlnPj.NormalizePem:(22)0≤Pe=Pemlnm!≤1.After many experiments and comparative analysis, the embedding dimensionm is selected as 8 and the delay time τ is selected as 5, and the entropy of the six modes of the vibration signal is obtained as shown in Table 2.Table 2
Permutation entropy of modal components of vibration signals.
Operating state
Group
Vibration entropy
V
1
V
2
V
3
V
4
V
5
V
6
Belt damage
1
0.6139
0.6753
0.6214
0.6457
0.6390
0.6523
2
0.6138
0.6751
0.6211
0.6457
0.6391
0.6522
Retarder oil leakage
1
0.3269
0.3658
0.3569
0.3546
0.3426
0.3261
2
0.3265
0.3659
0.3565
0.3546
0.3425
0.3260
Gear pitting peeling
1
0.5689
0.5469
0.5218
0.5634
0.5849
0.5459
2
0.5687
0.5468
0.5217
0.5634
0.5847
0.5458
Normal state
1
0.1232
0.1526
0.1726
0.1521
0.1358
0.1189
2
0.1231
0.1523
0.1725
0.1523
0.1356
0.1189The entropy value is related to the degree of the sample sequence rule. From Table2, it can be seen that there are obvious differences in the permutation entropy of each state, the distinguishability is good, and the state information of the retarder can be effectively characterized.


## 3.1. VMD Decomposition

### 3.1.1. VMD Decomposition Principle

The VMD decomposition is mainly divided into two parts: the establishment and solution of the variational constraint problem. The following problems are solved for the vibration signal with the data lengthN in the oil production process:(12)minuk,ωk∑k∂tδt+jπt∗μkte−jωkt2s.t.∑kμk=f,where μp≔μ1,…,μp is the decomposed p modes and ωp≔ω1,…,ωp is the center frequency of the p modes.In order to solve the optimal solution, the quadratic penalty factorα and the Lagrangian operator λt are introduced to change the constrained variational problem into a nonbinding variational problem. The extended Lagrangian expression is as follows:(13)Lμk,ωk,λ=α∑k∂tδt+jπt∗μkte−jωkt22+ft−∑kμkt22+λt,ft−∑kμkt.The saddle point of equation (13) is solved by the ADMM method, so that μkn+1, ωkn+1, and λn+1 are continuously updated, and the modal component μk and center frequency ωk are solved as follows:(14)μ^kn+1ω=f^ω−∑i≠kμ^knω+λ^ω/21+2αω−ωk2,(15)ωkn+1=∫0∞ωμ^knωdω∫0∞μ^knωdω.The VMD steps are as follows:(1)
Initializeμ^k1, ωk1, λ^1, and n, let its initial value be 0, and set the decomposition modal number K to 2 (predecomposition optimization).(2)
Updateμk and wk according to equations (14) and (15), respectively.(3)
Updateλ^:(16)λ^nω+τf^ω−∑μ^kn+1⟶λ^n+1ω.(4)
If the following formula is satisfied, the iteration is stopped and the result is output; otherwise, it returns to step 2:(17)∑μ^kn+1−μ^kn22μ^kn22<ε.


### 3.1.2. VMD ParameterK Optimization

In order to prevent the VMD from being decomposed, theK parameter is selected according to the energy conservation theory before and after decomposition. For the original retarder vibration signal sequence xi, the energy calculation formula is as follows:(18)E=∑i=1nx2in,where E represents the signal energy value and n is the sampling point. In order to characterize the energy difference before and after the VMD decomposition, the energy difference parameter ψ is defined and calculated as follows:(19)ψ=∑x=1KEx−EE×100%,where Ex corresponds to the energy of the x component, K is the number of components, and E is the original signal energy. The total energy of K components after VMD decomposition is equal to the energy of the original signal (i.e., the ideal value of ψ is 0), after many experiments and calculations, the trend of K is shown in Figure 5.Figure 5
Change trend ofK value.It can be seen from Figure5 that when K is greater than 6, the energy difference parameter λ is increased, and it can be judged that overdecomposition occurs. At this time, the K value at the turning point is the optimal decomposition mode number of the VMD. The time-frequency diagram obtained by decomposing the vibration signal in the belt damage state is shown in Figures 6 and 7.Figure 6
Time domain diagram of vibration signal VMD decomposition.Figure 7
Frequency domain diagram of vibration signal VMD decomposition.It can be seen from the frequency domain spectrum of Figure7 that the VMD decomposition of the vibration signal effectively improves the modal aliasing phenomenon and provides a strong support for accurate fault diagnosis.


## 3.1.1. VMD Decomposition Principle

The VMD decomposition is mainly divided into two parts: the establishment and solution of the variational constraint problem. The following problems are solved for the vibration signal with the data lengthN in the oil production process:(12)minuk,ωk∑k∂tδt+jπt∗μkte−jωkt2s.t.∑kμk=f,where μp≔μ1,…,μp is the decomposed p modes and ωp≔ω1,…,ωp is the center frequency of the p modes.In order to solve the optimal solution, the quadratic penalty factorα and the Lagrangian operator λt are introduced to change the constrained variational problem into a nonbinding variational problem. The extended Lagrangian expression is as follows:(13)Lμk,ωk,λ=α∑k∂tδt+jπt∗μkte−jωkt22+ft−∑kμkt22+λt,ft−∑kμkt.The saddle point of equation (13) is solved by the ADMM method, so that μkn+1, ωkn+1, and λn+1 are continuously updated, and the modal component μk and center frequency ωk are solved as follows:(14)μ^kn+1ω=f^ω−∑i≠kμ^knω+λ^ω/21+2αω−ωk2,(15)ωkn+1=∫0∞ωμ^knωdω∫0∞μ^knωdω.The VMD steps are as follows:(1)
Initializeμ^k1, ωk1, λ^1, and n, let its initial value be 0, and set the decomposition modal number K to 2 (predecomposition optimization).(2)
Updateμk and wk according to equations (14) and (15), respectively.(3)
Updateλ^:(16)λ^nω+τf^ω−∑μ^kn+1⟶λ^n+1ω.(4)
If the following formula is satisfied, the iteration is stopped and the result is output; otherwise, it returns to step 2:(17)∑μ^kn+1−μ^kn22μ^kn22<ε.


## 3.1.2. VMD ParameterK Optimization

In order to prevent the VMD from being decomposed, theK parameter is selected according to the energy conservation theory before and after decomposition. For the original retarder vibration signal sequence xi, the energy calculation formula is as follows:(18)E=∑i=1nx2in,where E represents the signal energy value and n is the sampling point. In order to characterize the energy difference before and after the VMD decomposition, the energy difference parameter ψ is defined and calculated as follows:(19)ψ=∑x=1KEx−EE×100%,where Ex corresponds to the energy of the x component, K is the number of components, and E is the original signal energy. The total energy of K components after VMD decomposition is equal to the energy of the original signal (i.e., the ideal value of ψ is 0), after many experiments and calculations, the trend of K is shown in Figure 5.Figure 5
Change trend ofK value.It can be seen from Figure5 that when K is greater than 6, the energy difference parameter λ is increased, and it can be judged that overdecomposition occurs. At this time, the K value at the turning point is the optimal decomposition mode number of the VMD. The time-frequency diagram obtained by decomposing the vibration signal in the belt damage state is shown in Figures 6 and 7.Figure 6
Time domain diagram of vibration signal VMD decomposition.Figure 7
Frequency domain diagram of vibration signal VMD decomposition.It can be seen from the frequency domain spectrum of Figure7 that the VMD decomposition of the vibration signal effectively improves the modal aliasing phenomenon and provides a strong support for accurate fault diagnosis.


## 3.2. Calculated Permutation Entropy

Permutation entropy can detect the sudden change of the signal, has strong antinoise ability and high time resolution, and is highly targeted to the nonstationary chaotic vibration signal under complex environment of oil production field.For the signal sequenceXi,i=1,2,…,n, phase space reconstruction is performed to obtain the following matrix:(20)x1x1+τ⋯x1+m−1τx2x2+τ⋯x2+m−1τxjxj+τ⋯xj+m−1τ⋯⋯xkxk+τ⋯xk+m−1τ,where j=1,2,…,k, k=n−m−1τ, and m and τ are the embedded dimension and delay time. The elements of each component (total of k) are arranged in ascending order according to the numerical value, and different symbol sequences are obtained according to the reconstructed component index. The m-dimensional phase space mapping symbol sequence has m! species, and the probability of occurrence of s different symbol sequences is P1,P2,…,Ps, respectively, calculated as follows:(21)Pem=−∑j=1sPjlnPj.NormalizePem:(22)0≤Pe=Pemlnm!≤1.After many experiments and comparative analysis, the embedding dimensionm is selected as 8 and the delay time τ is selected as 5, and the entropy of the six modes of the vibration signal is obtained as shown in Table 2.Table 2
Permutation entropy of modal components of vibration signals.
Operating state
Group
Vibration entropy
V
1
V
2
V
3
V
4
V
5
V
6
Belt damage
1
0.6139
0.6753
0.6214
0.6457
0.6390
0.6523
2
0.6138
0.6751
0.6211
0.6457
0.6391
0.6522
Retarder oil leakage
1
0.3269
0.3658
0.3569
0.3546
0.3426
0.3261
2
0.3265
0.3659
0.3565
0.3546
0.3425
0.3260
Gear pitting peeling
1
0.5689
0.5469
0.5218
0.5634
0.5849
0.5459
2
0.5687
0.5468
0.5217
0.5634
0.5847
0.5458
Normal state
1
0.1232
0.1526
0.1726
0.1521
0.1358
0.1189
2
0.1231
0.1523
0.1725
0.1523
0.1356
0.1189The entropy value is related to the degree of the sample sequence rule. From Table2, it can be seen that there are obvious differences in the permutation entropy of each state, the distinguishability is good, and the state information of the retarder can be effectively characterized.


## 4. Experiment and Result Analysis

### 4.1. GWO-SVM Diagnostic Model

SVM is especially suitable for small sample fault diagnosis, but the selection of penalty factorC and kernel function directly affect the classification performance of SVM. Considering that the sample in this paper is a multisensor data feature set, the sensor frequency response, environmental noise, and load output fluctuation interference have a great influence on the diagnostic accuracy, and GA, PSO, and other algorithms are easy to fall into the local optimal solution, so a mixed kernel function is constructed and the gray wolf optimization (GWO) algorithm is introduced to find the best penalty factor and mixing coefficient.It is necessary to combine the different kernel functions to obtain a kernel function with strong promotion ability and learning ability and good extrapolation ability. The Taylor-Kernel with moderate decreasing (T-KMOD) function satisfies the zero-point near-descent criterion, with good flexibility, fast convergence, and good locality. Polynomial functions have high classification accuracy and strong generalization ability. Therefore, the two are combined in the following form:(23)SM=λSpoly+1−λST−KMOF,where ST−KMOD=T−KMODx,x′=L∑i=1nγ/x−x′2+σ2k∗i, Spoly=x⋅x′+Cq, λ∈0,1 C≥0; x⋅x′∈Rn, L > 0, L is used to control the value of the kernel function at 0. σandγ are used to control the width and convergence rate of the kernel function, respectively, and both i and n are positive integers.GWO has a fast convergence speed and a simple structure, and it is easier to achieve optimal classification. The mathematical model is as follows:(24)B=C⋅Xpt−Xt,Xt+1=Xpt−A⋅B,where t represents the current number of iterations; Xpt is the prey position vector, and A and C are coefficient vectors. A and C are as follows:(25)A=2a⋅r1−a,C=2r2,where a is the convergence factor, satisfying a∈0,2, and r1 and r2 are random vectors in [0, 1].60 sets of samples were collected for each state, 40 sets were used for training, and 20 sets were used for testing. Initialize the population, calculate the fitness value of each wolf, select the first three optimal fitness values, determine the gray wolf rank, and update the gray wolf position, head wolf, coefficient vector, and other parameters until the SVM parameters are optimal. The algorithm flowchart is shown in Figure8:Figure 8
Gray wolf optimization algorithm flowchart.


### 4.2. Experimental Result

In the experiment, the MVP-6000 acquisition card of ADLINK Company is used, whose sampling rate was set to 40 kHz. The sound sensor is placed about 50 cm away from the sound source, and the vibration sensor is adsorbed on the surface of the vibration body of the retarder. The acquisition card is equipped with IEPE constant current source, ±10 V voltage range, 24 bit resolution, and 110 dB dynamic range.The sound texture feature extracted by the generalized S-transform and the vibration entropy feature calculated by the VMD decomposition are combined to construct a joint eigenvector matrix and sent to the SVM for training. GWO is used to optimize the penalty factorC and the mixing coefficient λ in order to improve the SVM classification performance. The number of iterations is set to 100, and the optimal parameter values are 2.5693 and 0.17, respectively. The convergence of the algorithm is shown in Figure 9:Figure 9
Convergence of GWO algorithm.The sound-vibration joint eigenvectors are constructed with a combination of Tables1 and 2. Set the normal state sample characteristic value label to 1 (1–20 groups), the belt damage is 2 (21–40 groups), the retarder oil leakage is 3 (41–60 groups), and the gear pitting peeling is 4 (61–80 groups). Some test data are shown in Table 3.Table 3
Partial test samples and diagnostic results.
Joint eigenvector
Diagnosis label
Diagnosis result
D
F
dir
Fc
V1
V2
V3
V4
V5
V6
4.6393
5.4331
4.3423
0.1233
0.1525
0.1728
0.1523
0.1360
0.1190
1
Normal state
4.6385
5.4348
4.3422
0.1229
0.1523
0.1723
0.1523
0.1358
0.1193
1
2.5673
4.5897
1.2485
0.6137
0.6752
0.6213
0.6459
0.6391
0.6527
2
Belt damage
2.5665
4.5813
1.2430
0.6139
0.6753
0.6209
0.6458
0.6392
0.6523
2
3.8371
8.3847
2.2431
0.3271
0.3661
0.3570
0.3541
0.3428
0.3265
3
Retarder oil leakage
3.8365
8.3859
2.2425
0.3267
0.3660
0.3566
0.3542
0.3427
0.3261
3
5.4276
6.3346
3.3395
0.5692
0.5468
0.5219
0.5633
0.5851
0.5463
4
Gear pitting peeling
5.4282
6.3361
3.3325
0.5685
0.5471
0.5219
0.5639
0.5849
0.5465
4The diagnosis results are as follows, in which Figures10 and 11 show the fault identification results before and after the model optimization.Figure 10
Fault identification results of the original SVM model.Figure 11
Fault identification results of the GWO-SVM model.According to the results of Figure10, two groups of normal state samples were misjudged as belt damage state, two groups of retarder oil leakage state samples were misjudged as belt damage state, and one group of gear pitting peeling state samples was misjudged as retarder oil leakage state. The identification accuracy is 93.75%. The results of the GWO-SVM model are shown in Figure 11. All 80 groups of test samples are classified correctly, and the identification accuracy is 100%. Although the input characteristics are the same, the diagnostic accuracy differs significantly. Compared with SVM, the accuracy of GWO-SVM is improved by 6.25%. The reason is that GWO-SVM diagnostic model can find the appropriate parameters of SVM classifier through GWO, which makes full use of the classification advantages of SVM in constructing optimal hyperplane, and thus, the accuracy of fault diagnosis is improved. At the same time, the important influence of the accurate selection of penalty factor C and kernel function on the accuracy of state identification is verified. Compared with the SVM model, GWO-SVM can more accurately reflect the type of retarder defects.


### 4.3. Diagnostic Effect Comparison and Verification

#### 4.3.1. Sound-Vibration Joint Characteristic Method Verification

Contrast the diagnosis effect of sound signal, vibration signal, sound-vibration combined, as shown in Figure12. For the experimental sample data, the diagnostic accuracy of sound characteristics and vibration characteristics were 91% and 94%, respectively, and the sound-vibration combined diagnosis accuracy rate reached 100%. Therefore, based on the combined characteristics of sound and vibration, the state information of the retarder can be fully reflected, the extracted eigenvectors are complementary, and the diagnostic effect is improved.Figure 12
Comparison of three eigenvectors.


#### 4.3.2. Verification of Generalization Performance

Because the source and structure of the data are different in the actual oilfield operation, it is necessary to classify the fault data of the same type and different characteristics. In the generalization experiment, the sampling rate was changed from 40 kHz to 30 kHz, and the PCB357B21 type vibration sensor and the WM-025N type pickup were replaced, and the sensor placement position was changed. The diagnosis result is shown in Figure13.Figure 13
Accuracy comparison in four states.It can be seen from Figure13 that the overall diagnostic accuracy of the optimized SVM model still reaches 97.8% in the case of changes in acquisition parameter settings, sensor types, and positions, which is much higher than the unoptimized model, indicating that the optimized model has stronger adaptability to fresh samples and better generalization ability.


## 4.1. GWO-SVM Diagnostic Model

SVM is especially suitable for small sample fault diagnosis, but the selection of penalty factorC and kernel function directly affect the classification performance of SVM. Considering that the sample in this paper is a multisensor data feature set, the sensor frequency response, environmental noise, and load output fluctuation interference have a great influence on the diagnostic accuracy, and GA, PSO, and other algorithms are easy to fall into the local optimal solution, so a mixed kernel function is constructed and the gray wolf optimization (GWO) algorithm is introduced to find the best penalty factor and mixing coefficient.It is necessary to combine the different kernel functions to obtain a kernel function with strong promotion ability and learning ability and good extrapolation ability. The Taylor-Kernel with moderate decreasing (T-KMOD) function satisfies the zero-point near-descent criterion, with good flexibility, fast convergence, and good locality. Polynomial functions have high classification accuracy and strong generalization ability. Therefore, the two are combined in the following form:(23)SM=λSpoly+1−λST−KMOF,where ST−KMOD=T−KMODx,x′=L∑i=1nγ/x−x′2+σ2k∗i, Spoly=x⋅x′+Cq, λ∈0,1 C≥0; x⋅x′∈Rn, L > 0, L is used to control the value of the kernel function at 0. σandγ are used to control the width and convergence rate of the kernel function, respectively, and both i and n are positive integers.GWO has a fast convergence speed and a simple structure, and it is easier to achieve optimal classification. The mathematical model is as follows:(24)B=C⋅Xpt−Xt,Xt+1=Xpt−A⋅B,where t represents the current number of iterations; Xpt is the prey position vector, and A and C are coefficient vectors. A and C are as follows:(25)A=2a⋅r1−a,C=2r2,where a is the convergence factor, satisfying a∈0,2, and r1 and r2 are random vectors in [0, 1].60 sets of samples were collected for each state, 40 sets were used for training, and 20 sets were used for testing. Initialize the population, calculate the fitness value of each wolf, select the first three optimal fitness values, determine the gray wolf rank, and update the gray wolf position, head wolf, coefficient vector, and other parameters until the SVM parameters are optimal. The algorithm flowchart is shown in Figure8:Figure 8
Gray wolf optimization algorithm flowchart.


## 4.2. Experimental Result

In the experiment, the MVP-6000 acquisition card of ADLINK Company is used, whose sampling rate was set to 40 kHz. The sound sensor is placed about 50 cm away from the sound source, and the vibration sensor is adsorbed on the surface of the vibration body of the retarder. The acquisition card is equipped with IEPE constant current source, ±10 V voltage range, 24 bit resolution, and 110 dB dynamic range.The sound texture feature extracted by the generalized S-transform and the vibration entropy feature calculated by the VMD decomposition are combined to construct a joint eigenvector matrix and sent to the SVM for training. GWO is used to optimize the penalty factorC and the mixing coefficient λ in order to improve the SVM classification performance. The number of iterations is set to 100, and the optimal parameter values are 2.5693 and 0.17, respectively. The convergence of the algorithm is shown in Figure 9:Figure 9
Convergence of GWO algorithm.The sound-vibration joint eigenvectors are constructed with a combination of Tables1 and 2. Set the normal state sample characteristic value label to 1 (1–20 groups), the belt damage is 2 (21–40 groups), the retarder oil leakage is 3 (41–60 groups), and the gear pitting peeling is 4 (61–80 groups). Some test data are shown in Table 3.Table 3
Partial test samples and diagnostic results.
Joint eigenvector
Diagnosis label
Diagnosis result
D
F
dir
Fc
V1
V2
V3
V4
V5
V6
4.6393
5.4331
4.3423
0.1233
0.1525
0.1728
0.1523
0.1360
0.1190
1
Normal state
4.6385
5.4348
4.3422
0.1229
0.1523
0.1723
0.1523
0.1358
0.1193
1
2.5673
4.5897
1.2485
0.6137
0.6752
0.6213
0.6459
0.6391
0.6527
2
Belt damage
2.5665
4.5813
1.2430
0.6139
0.6753
0.6209
0.6458
0.6392
0.6523
2
3.8371
8.3847
2.2431
0.3271
0.3661
0.3570
0.3541
0.3428
0.3265
3
Retarder oil leakage
3.8365
8.3859
2.2425
0.3267
0.3660
0.3566
0.3542
0.3427
0.3261
3
5.4276
6.3346
3.3395
0.5692
0.5468
0.5219
0.5633
0.5851
0.5463
4
Gear pitting peeling
5.4282
6.3361
3.3325
0.5685
0.5471
0.5219
0.5639
0.5849
0.5465
4The diagnosis results are as follows, in which Figures10 and 11 show the fault identification results before and after the model optimization.Figure 10
Fault identification results of the original SVM model.Figure 11
Fault identification results of the GWO-SVM model.According to the results of Figure10, two groups of normal state samples were misjudged as belt damage state, two groups of retarder oil leakage state samples were misjudged as belt damage state, and one group of gear pitting peeling state samples was misjudged as retarder oil leakage state. The identification accuracy is 93.75%. The results of the GWO-SVM model are shown in Figure 11. All 80 groups of test samples are classified correctly, and the identification accuracy is 100%. Although the input characteristics are the same, the diagnostic accuracy differs significantly. Compared with SVM, the accuracy of GWO-SVM is improved by 6.25%. The reason is that GWO-SVM diagnostic model can find the appropriate parameters of SVM classifier through GWO, which makes full use of the classification advantages of SVM in constructing optimal hyperplane, and thus, the accuracy of fault diagnosis is improved. At the same time, the important influence of the accurate selection of penalty factor C and kernel function on the accuracy of state identification is verified. Compared with the SVM model, GWO-SVM can more accurately reflect the type of retarder defects.


## 4.3. Diagnostic Effect Comparison and Verification

### 4.3.1. Sound-Vibration Joint Characteristic Method Verification

Contrast the diagnosis effect of sound signal, vibration signal, sound-vibration combined, as shown in Figure12. For the experimental sample data, the diagnostic accuracy of sound characteristics and vibration characteristics were 91% and 94%, respectively, and the sound-vibration combined diagnosis accuracy rate reached 100%. Therefore, based on the combined characteristics of sound and vibration, the state information of the retarder can be fully reflected, the extracted eigenvectors are complementary, and the diagnostic effect is improved.Figure 12
Comparison of three eigenvectors.


### 4.3.2. Verification of Generalization Performance

Because the source and structure of the data are different in the actual oilfield operation, it is necessary to classify the fault data of the same type and different characteristics. In the generalization experiment, the sampling rate was changed from 40 kHz to 30 kHz, and the PCB357B21 type vibration sensor and the WM-025N type pickup were replaced, and the sensor placement position was changed. The diagnosis result is shown in Figure13.Figure 13
Accuracy comparison in four states.It can be seen from Figure13 that the overall diagnostic accuracy of the optimized SVM model still reaches 97.8% in the case of changes in acquisition parameter settings, sensor types, and positions, which is much higher than the unoptimized model, indicating that the optimized model has stronger adaptability to fresh samples and better generalization ability.


## 4.3.1. Sound-Vibration Joint Characteristic Method Verification

Contrast the diagnosis effect of sound signal, vibration signal, sound-vibration combined, as shown in Figure12. For the experimental sample data, the diagnostic accuracy of sound characteristics and vibration characteristics were 91% and 94%, respectively, and the sound-vibration combined diagnosis accuracy rate reached 100%. Therefore, based on the combined characteristics of sound and vibration, the state information of the retarder can be fully reflected, the extracted eigenvectors are complementary, and the diagnostic effect is improved.Figure 12
Comparison of three eigenvectors.


## 4.3.2. Verification of Generalization Performance

Because the source and structure of the data are different in the actual oilfield operation, it is necessary to classify the fault data of the same type and different characteristics. In the generalization experiment, the sampling rate was changed from 40 kHz to 30 kHz, and the PCB357B21 type vibration sensor and the WM-025N type pickup were replaced, and the sensor placement position was changed. The diagnosis result is shown in Figure13.Figure 13
Accuracy comparison in four states.It can be seen from Figure13 that the overall diagnostic accuracy of the optimized SVM model still reaches 97.8% in the case of changes in acquisition parameter settings, sensor types, and positions, which is much higher than the unoptimized model, indicating that the optimized model has stronger adaptability to fresh samples and better generalization ability.


## 5. Conclusion

Defect identification of related retarder of beam oil pumping machine has always been a technical problem in the state monitoring of distributed oil production wells. The combination of sound texture and vibration entropy characteristics and the GWO-SVM classification algorithm, proposed in this paper, can effectively and accurately diagnose the field faults under the complementary frequency band. The main contributions and novels of the proposed method are summarized as follows:(1)
A fault diagnosis method based on the complementarity combination of sound-vibration signals is proposed for retarder equipment of oil pumping machine, which improves the accuracy of fault identification on the basis of nonmissing detection of retarder defects.(2)
For sound signals, the box dimension, directivity, and contrast ratio of time-frequency diagram are calculated after generalized S-transform to construct the sound texture features. For vibration signals, parameterK of the VMD method is selected by the energy method, and the permutation entropy of modal components is obtained to construct the vibration entropy characteristics. The combination of the two effectively characterizes comprehensive information on various types of fault samples.(3)
The diagnosis model of GWO-SVM proposed in this paper, which optimizes the penalty factor and mixing coefficient, greatly improves the accuracy and generalization of SVM identification. It has broad application value in fault diagnosis of oil well machinery and equipment.



---

*Source: 2709384-2020-05-11.xml*