Extraction method of composite fault features of gear transmission system based on demodulation of NMD and Teager energy operators

Abstract

In order to effectively identify and extract the composite fault characteristics of the gear transmission system, a composite fault diagnosis method combining nonlinear mode decomposition (NMD) and Teager energy operator demodulation is proposed. Because the envelope demodulation of Hilbert transform has the disadvantages of large amount of calculation and end effect, it uses Teager energy operator to solve the problem of large amount of calculation, and NMD solves the problem that the fault signal features is not easy to extract under the mode aliasing. First, the NMD method is used to decompose the fault simulation signal, and the nonlinear modal component with practical physical significance is obtained. Secondly, the Teager energy operator demodulation is carried out for the nonlinear modal components, and the demodulation results are analyzed to verify the feasibility of the method. Then, the method is applied to the composite fault diagnosis of gear pitting wear in gear transmission system, and the characteristic frequency obtained from the test data is compared with the calculated characteristic frequency. The comparative analysis shows that the method can separate the fault characteristic frequencies of large and small gears. The comparative analysis with EMD and EEMD methods in simulation signal analysis and experimental research shows that this method is superior.

Keywords

Nonlinear mode decomposition composite fault of gear transmission system feature extraction Teager energy operator demodulation

Introduction

Gearbox is an integral part of most mechanical products. In order to prevent accidental failure and possible life loss, the fault diagnosis of gearbox has received great attention.^1,2 The vibration signal of gearbox is usually complex, non-stationary and nonlinear.³ When a composite fault occurs, multiple faults will interfere with each other and be coupled with each other.⁴ In recent years, time-frequency analysis methods have been widely used in composite fault diagnosis. Wang et al.⁵ proposed a composite fault diagnosis method for Gearbox Based on improved integrated local mean decomposition. Through the analysis of simulation and experimental signals, it can be seen that the method effectively overcomes the mode mixing phenomenon and reduces the reconstruction error. At the same time, the method avoids the appearance of pseudo components and reduces the calculation amount. It is verified that the improved integrated local mean decomposition method is an effective method to extract composite fault features. Isham et al.⁶ proposed the VMDEA and ELM combined fan gearbox fault diagnosis method. This study provides a good solution for determining optimized VMD parameters to decompose vibration signals, and also provides a more effective fault diagnosis method for wind turbine gearboxes. Li et al.⁷ proposed a gearbox composite fault diagnosis method based on SWD-AVDIF. By analyzing the simulation signals and experimental signals of the gearbox composite fault, the results showed that the method could effectively separate the composite fault signals and strengthen the fault characteristics. Yang et al.⁸ put forward a gear fault diagnosis method combining EEMD and Choi Williams. Through the analysis of simulation and experimental signals, the adaptability and effectiveness of this method in gearbox fault diagnosis were verified. The traditional time-domain or frequency-domain analysis method is only suitable for linear stationary signals. The vibration signals of composite faults are mostly non-linear and non-stationary signals. Using only time-domain or frequency-domain analysis methods often fails to effectively extract fault features.⁹ Therefore, a time-frequency analysis method should be used in feature extraction of composite faults. Windowed Fourier transform (WFT) and wavelet transform (WT) are more commonly used time-frequency analysis methods in signal processing. Both are suitable for processing nonlinear, non-stationary signals.^10,11 The difference between them is that they analyze the form of the components that appear in the signal. WFT has linear frequency resolution, that is, it is distinguished according to the frequency difference of each component. WT has logarithmic frequency resolution and is distinguished according to the component frequency ratio. According to Heisenberg’s uncertainty principle, for any kind of time-frequency transformation, the frequency resolution and time resolution cannot reach the optimum at the same time. The time resolution of the WT will increase as the frequency increases. The time resolution of the WFT is the same at all frequencies. So we merge the two, Using WT to get high temporal resolution, Use WFT for high frequency resolution. Huang first proposed a new time-frequency analysis method-empirical mode decomposition (EMD). It can decompose complex signals into finite intrinsic mode functions (IMF), which are suitable for nonlinear and non-stationary signals.¹² Since this method came out, it is widely used in fault diagnosis of mechanical equipment.^13,14 Wu and Huang¹⁵ also proposed an EEMD method, adding white noise to the signal and then performing EMD decomposition. However, the signals decomposed by EMD and EEMD are prone to produce pseudo components. The non-sinusoidal complex components generated by the same excitation source will be decomposed into multiple simple components, which will increase the number of components. In addition, they are sensitive to noise and their robustness is not ideal. In view of their disadvantages, Iatsenko et al.¹⁶ proposed a new time-frequency analysis method—nonlinear mode decomposition (NMD), which can decompose the given signal into oscillation with practical physical significance and decompose it into nonlinear mode component (NM). In the process of decomposition, the advantages of WT and WFT are integrated, the accuracy of extraction and reconstruction of main signals is improved, and the noise robustness is very strong. At present, the envelope demodulation method used in many domestic and foreign literatures of fault diagnosis is Hilbert transform envelope demodulation.^17,18 But there are also endpoint effects in Hilbert transform, which will affect the whole signal analysis results.¹⁹ Hilbert has a lot of calculations, the calculation of the Teager energy operator is relatively small. The calculation amount of Teager energy operator demodulation is mainly 4N times real number multiplication, equivalent to N times complex multiplication; The calculation amount of Teager energy operator demodulation is about 1/logN of Hilbert demodulation method.²⁰ When N = 1024, it is 1/10. Teager energy operator can enhance the transient characteristics of the signal and is suitable for detecting the impact component in the signal. This method has high time resolution, and it has good adaptive ability to the instantaneous change of signal. It has low computational complexity and high efficiency.²¹

In this paper, a composite fault feature extraction method combining nonlinear mode decomposition (NMD) and Teager energy operator is proposed. The composite fault signal is decomposed by NMD, and then the obtained nonlinear modal component (NM) is demodulated by Teager energy operator, and the fault type is determined by comparing the obtained fault characteristic frequency with the calculated fault characteristic frequency.

The second part of this paper introduces the basic theory of the proposed method. The third part analyzes the simulation signal. The fourth part compares and analyzes the effectiveness of the three methods for composite fault feature extraction. The fifth part expounds the conclusion of this paper.

The proposed method

Basic theory of NMD

Under a certain excitation, nonlinear system will produce response, which can be regarded as NM component. It is assumed that the signal consists of a series of NM components and noise.

s (t) = \sum_{i} c_{i} (t) + η (t)

(1)

Among them, $η (t)$ is noise, and $c_{i} (t)$ is NM component. $c_{i} (t)$ can be expressed as a series of FM and AM harmonics:

c (t) = A (t) v [φ (t)] = \sum_{h} A_{h} \cos φ_{h}

(2)

r_{h} = A_{h} \cos φ_{h}

(3)

Where $v [φ (t)]$ period is 2 π. $r_{h} (t)$ is the h-th harmonic of $c_{i} (t)$ . $r_{1} (t)$ is the main harmonic. The rest are subharmonics (h = 2, 3, 4 …).

The purpose of NMD decomposition is to decompose the NM signal. The characteristic parameters of each harmonic $r_{h} (t)$ are obtained. Amplitude is $A_{h} (t)$ , phase is $φ_{h} (t)$ , frequency is $v_{h} (t) = φ_{h}^{'} (t)$ .

To construct WT and WFT:

(1) Wavelet transform the original signal. Given the signal s(t), WT can be expressed as:

\begin{matrix} w_{s} (ω, t) \equiv \int_{- \infty}^{\infty} s^{+} (μ) ψ^{*} (\frac{ω (μ - t)}{ω_{ψ}}) \frac{ω d_{μ}}{ω_{ψ}} \\ = \frac{1}{2 π} \int_{0}^{\infty} e^{i ε t} \hat{s} (ε) {\hat{ψ}}^{*} (\frac{ω_{ψ} ε}{ω}) d ε \end{matrix}

(4)

In the formula, $ω$ is the frequency, $μ$ and $ε$ are the other forms of time t in the time and frequency domains respectively; $ω_{ψ}$ is the peak frequency of wavelet; $ψ_{t}$ is a normal wavelet basis function; $\hat{ψ}$ and $ψ^{*}$ are the Fourier transform and complex conjugate function of $ψ (t)$ respectively.

(2) Extraction time-frequency ridge $ω_{p}$

A curve consisting of a series of ridge points is called a ridge curve. Represents the peak frequency of each moment. In the paper, the ridge extraction method proposed by Iatsenko et al.²² is adopted. The ridge $r_{1} (t)$ of the main harmonic $ω_{p} (t)$ is a sequence of peaks that maximizes the function path. For wavelet transform (WT), the path function can be expressed as:

S [ω_{p} (t)] \equiv \sum_{n = 1}^{N} [In | G_{s} (ω_{p} (t_{n}), t_{n}) | - \frac{1}{2} {(\frac{Δ ω_{p} (t_{n}) - < Δ In ω_{p} >}{std | Δ In ω_{p} |})}^{2} - \frac{1}{2} {(\frac{ω_{p} (t_{n}) - < In ω_{p} >}{std [In ω_{p}]})}^{2}]

(5)

Where, $G_{s} (ω, t)$ is the formula represented by equation (4). $ω_{p} (t_{n})$ is the expression of ridge $ω_{p}$ when $t = t_{n}$ , $t_{n}$ is the n-th peak; $< f (t) >$ and $std [f (t)]$ represent the average and standard deviation of $f (t)$ respectively.

(3) Main harmonic reconstruction

After the ridgeline curve is extracted, the main harmonic $r_{1} (t)$ needs to be reconstructed. There are two methods to reconstruct the main harmonic: the direct method and the ridge method. Because the ridge method has good robustness. Therefore, we first use ridge method to reconstruct the main harmonic. For wavelet transform, the ridge method can be expressed as:

\begin{matrix} v^{r} (t) = ω_{p} (t) \exp (δ In v_{d} (t)) \\ A^{r} (t) \exp (i φ (t)) = \frac{2 ω_{s} (ω_{p} (t), t)}{{\hat{ψ}}^{*} (\frac{ω_{ψ} v (t)}{ω_{p} (t)})} \end{matrix}}

(6)

(4) WT or WFT selection criteria

The paper combines the advantages of WT and WFT. It enables NMD to adaptively select the most suitable type of time-frequency transformation. The empirical formula can be expressed as:

R = (1 + V [\partial_{t} v (t), v (t)])^{- 1} + (1 + V [\partial_{t} A (t), v (t)])^{- 1}

(7)

When R < 1, WFT is used, and WT is used otherwise.

(5) Determination of reconstruction method

Since there are two reconstruction methods, in order to make the reconstructed harmonics more accurate and reasonable. According to different time-frequency analysis methods and two reconstruction methods, the final reconstruction method is determined.

(1) When WFT is the optimal time-frequency transform, the given signal G(s) is transformed by WFT:

\begin{matrix} G_{s} (ω, t) \equiv \int_{- \infty}^{+ \infty} s^{+} (μ) g (μ - t) \exp (- i ω (μ - t)) d_{μ} \\ = \frac{1}{2 π} \int_{0}^{\infty} \exp (i ε t) \hat{s} (ε) \hat{g} (ω - ε) d_{ε} \end{matrix}

(8)

Where, $g (t)$ is the Gaussian window function. The time-frequency ridge can be extracted according to equation (5). For WFT, according to the difference between the two, it is only necessary to change the logarithmic frequency scale $In ω_{p}$ in equation (5) to $ω_{p}$ . Then the parameters of $A^{(d, r)} (t)$ , $φ^{(d, r)} (t)$ , $v^{(d, r)} (t)$ , respectively, using the direct method and ridge line method refactoring (d represent direct method, r represent direct method):

\begin{matrix} v^{r} (t) = ω_{p} (t) + δ v_{d} (t) \\ A^{r} (t) e^{i φ (t)} = \frac{2 G_{s} (ω_{p} (t), t)}{\hat{g} [ω_{p} (t) - v (t)]} \end{matrix}}

(9)

\begin{matrix} A^{d} (t) e^{i φ (t)} = C_{g}^{- 1} \int_{ω_{- (t)}}^{ω_{+ (t)}} G_{s} (ω, t) d_{ω} \\ C_{g} \equiv \frac{1}{2} \int_{- \infty}^{\infty} \hat{g} (ε) d ε = π g (0) \\ v^{d} (t) = Re [\frac{\int_{ω_{- (t)}}^{ω_{+ (t)}} ω G_{s} (ω, t) d_{ω}}{\int_{ω_{- (t)}}^{ω_{+ (t)}} G_{s} (ω, t) d_{ω}} - \bar{ω_{g}}] \\ \bar{ω_{g}} = \frac{C_{g}^{- 1}}{2} \int_{- \infty}^{+ \infty} ε \hat{g} (ε) d ε \end{matrix}}

(10)

Where, δ is the discrete correction coefficient.

First, calculate the weighted average ${\tilde{A}}^{(d, r)}$ , ${\tilde{φ}}^{(d, r)}$ , ${\tilde{v}}^{(d, r)}$ . The reconstruction method is determined according to the inconsistency of the reconstruction method. For harmonic parameters, the inconsistency can be expressed as:

\begin{matrix} ε_{A}^{(d, r)} \equiv k_{A}^{(d, r)} \sqrt{< {[{\hat{A}}^{(d, r)} (t) - A^{(d, r)} (t)]}^{2} >} \\ ε_{φ}^{(d, r)} \equiv k_{φ}^{(d, r)} \sqrt{1 - {| < \exp (i ({\tilde{φ}}^{(d, r)} (t) - {\tilde{φ}}^{(d, r)} (t))) > |}^{2}} \\ ε_{v}^{(d, r)} \equiv k_{v}^{(d, r)} \sqrt{< {({\tilde{v}}^{(d, r)} (t) - v^{(d, r)} (t))}^{2}} > \end{matrix}}

(11)

Where $k_{A, φ, v}^{(d, r)}$ is the coefficient of adjustment performance, $k_{A, φ, v}^{(d)}$ = ${3.4.2}$ , $k_{A, φ, v}^{(r)}$ =1.

(2) According to formula (11), if $ε_{A}^{d} < ε_{A}^{r}$ , the direct method is used, otherwise the ridge method is used. According to the above, WT is set as the best, the direct method is used to reconstruct the harmonic parameters.

\begin{matrix} A^{d} (t) e^{i φ (t)} = C_{ψ}^{- 1} \int_{ω_{-} (t)}^{ω_{+} (t)} ω_{s} (ω, t) \frac{d_{ω}}{ω} \\ C_{ψ} \equiv \frac{1}{2} \int_{0}^{\infty} {\hat{ψ}}^{*} (ε) \frac{d_{ε}}{ε} \\ v^{d} (t) = R_{e} \frac{D_{ψ}^{- 1} \int_{ω_{-} (t)}^{ω_{+} (t)} ω W_{s} (ω, t) \frac{d_{ω}}{ω}}{C_{ψ}^{- 1} \int_{ω_{-} (t)}^{ω_{+} (t)} W_{s} (ω, t) \frac{d_{ω}}{ω}} \\ D_{ψ} \equiv \frac{ω_{ψ}}{2} \int_{0}^{\infty} \frac{1}{ε} {\hat{ψ}}^{*} (ε) \frac{d_{ε}}{ε} \end{matrix}}

(12)

(6) Extraction of sub-harmonics

After obtaining the main harmonic $r_{1} (t)$ , its instantaneous frequency is $v^{(1)} (t)$ . For the sub-harmonic $r_{h} (t)$ , its ridge line $ω_{p}^{(h)}$ is the closest to the peak sequence of $h v^{1} (t)$ in the direction of amplitude increase. After finding the sub-harmonic ridge line $ω_{p}^{(h)}$ , the direct method and the ridge method are used to reconstruct the sub-harmonics $r_{h} (t)$ . Then, the optimal sub-harmonic parameter is determined by equation (11).

(7) Identification of true harmonics

The amplitude, phase, and frequency of true harmonics should be related to the main harmonic. After extracting the primary harmonic $r_{1} (t)$ and the secondary harmonic $r_{h} (t)$ , we need to determine if they are true harmonics. By using the method of time shift instead of data, the independent zero hypothesis test of the main harmonic $r_{1} (t)$ and the sub-harmonic $r_{h} (t)$ is carried out, so as to identify the authenticity of the harmonic.²³ The alternative parameter of the main harmonic $r_{1} (t)$ can be expressed as:

\begin{matrix} A_{d}^{(1)} (τ) = A^{(1)} (τ - \frac{Δ T_{d}}{2}) \\ φ_{d}^{(1)} (τ) = φ^{(1)} (τ - \frac{Δ T_{d}}{2}) \\ v_{d}^{(1)} (τ) = v^{(1)} (τ - \frac{Δ T_{d}}{2}) \\ τ = {t_{i = 1 + \frac{M}{2}, \dots, N - M / 2}} \\ Δ T_{d = 1, 2, 3 \dots N_{d}} = M (1 - \frac{2 d}{N_{d}}) / (2 f_{s}) \end{matrix}}

(13)

Where, $N$ is the sampling length, $f_{s}$ is the sampling frequency, $d$ is the number of time shift points of the substitute data and $M$ is the maximum time shift of the substitute data.

Alternative parameter process of subharmonic $r_{h} (t)$ :

(1) Move its signal forward in the time domain by $T_{d} / 2$ ; $G (ω, τ + Δ T_{d} / 2)$ or $W (ω, τ + Δ T_{d} / 2)$ .

(2) Use step (6) to select the optimal substitute harmonic parameters $A_{d}^{(h)}$ , $φ_{d}^{(h)}$ , $v_{d}^{(h)}$ .

(3) Introduce the statistical feature q_(A,v)^((h)) to measure the degree of consistency between the primary and secondary harmonics. Its dependence can be expressed as:

\begin{matrix} q_{A}^{(h)} \equiv \exp {- \frac{\sqrt{< {[A^{(h)} (t) < A^{(1)} (t) > - A^{(1)} (t) < A^{(h)} (t) >]}^{2} >}}{< A^{(1)} (t) A^{(h)} (t) >}} \\ q_{φ}^{(h)} \equiv a | < \exp {i [φ^{(h)} (t) - h φ^{(1)} (t)]} > | \\ q_{v}^{(h)} \equiv \exp {- \frac{\sqrt{< {[v^{(h)} (t) - h v^{(1)} (t)]}^{2} >}}{< v^{(h)} (t) >}} \end{matrix}}

(14)

Consistency coefficient can be expressed as:

ρ^{(h)} (ω_{A,} ω_{φ,} ω_{v,}) = (q_{A}^{(h)})^{ω_{A}} (q_{φ}^{(h)})^{ω_{φ}} (q_{v}^{(h)})^{ω_{v}}

(15)

The parameter $ω_{A, φ, v}^{(h)}$ is the weight of each consistency coefficient $q_{A, φ, v}^{(h)}$ . $ρ_{0}^{(h)}$ is the consistency coefficient when $Δ T_{0} = 0$ .

(4) In order to reduce the error caused by harmonic interference, the previous true harmonic must be subtracted and then harmonic identification (it is also applicable to the main harmonic. When identifying the sub-harmonic $r_{h} (t)$ , the main harmonic $r_{1} (t)$ should be subtracted).

(8) Noise detection.

The extracted true harmonics are added to obtain the NM component, which is then subtracted from the original signal. Until there is only noise in the remaining signal. The surrogate test method with Fourier transform is used to detect whether the residual signal contains noise. The statistical characteristic quantity D of the alternative data is expressed as:

\begin{matrix} D (α_{A}, α_{v}) \equiv α_{A} Q [\hat{A} (ε)] + α_{v} Q [\hat{v} (ε)] \\ Q [f (x)] \equiv - \int \frac{{| f (x) |}^{2}}{\int^{{| f (x) |}^{2}} d_{x}} In \frac{{| f (x) |}^{2}}{\int^{{| f (x) |}^{2}} d_{x}} d_{x} \end{matrix}}

(16)

Where, $Q [\hat{A} (ε)]$ and $Q [\hat{v} (ε)]$ are the discrete spectral entropy of amplitude and frequency respectively. The inspection steps are as follows:

(1) Calculate the statistical characteristic quantity $D_{0} (α_{A}, α_{v})$ .

(2) Create N_s Fourier transform alternatives, and calculate the corresponding $D_{s = 1, 2, \dots, N_{s}} (α_{A}, α_{v})$ .

(3) If the confidence level of $D_{s} > D_{0}$ is greater than or equal to the preset confidence level, the noise test zero hypothesis is rejected. D(1,1), D(0,1), and D (1,0) are used as discriminant statistics to test three different ones $(α_{A}, α_{v})$ . If at least one of them negates the zero hypothesis, it can be considered that the signal is not noise, and the test should be continued. The NMD flowchart is as follows (Figure 1):

Figure 1.

NMD flow chart.

Basic theory of Teager energy operator

The Teager energy operator is a non-linear operator, and the energy operator method is developed on the signal form of equation (17):

x (t) = a (t) \cos [φ (t)]

(17)

An energy operator is defined as follows:

Ψ (x) = {\overset{\cdot}{x}}^{2} - x \overset{\cdot\cdot}{x}

(18)

$\overset{\cdot}{x}$ and $\overset{\cdot\cdot}{x}$ are the first reciprocal and second derivative of x(t), respectively. There are:

\begin{matrix} Ψ [x (t)] = [a (t) \cos [φ (t)]^{'}]^{2} + a^{2} (t) φ^{″} (t) * \frac{\sin [2 φ (t)]}{2} \\ + \cos^{2} [φ (t)] Ψ [a (t)] \end{matrix}

(19)

In general, the variation of the modulated signal is much slower than that of the carrier signal, so the amplitude-frequency value can be treated approximately as a constant. Thus there are:

Ψ [a (t)] \approx 0, φ^{″} (t) \approx 0

(20)

φ [x (t)] \approx {[a (t) φ^{'} (t)]}^{2} = a^{2} (t) ω^{2} (t)

(21)

In the same way:

Ψ [x (t)] \approx a^{2} (t) ω^{4} (t)

(22)

The instantaneous amplitude and phase of the modulated signal are:

| a (t) | \approx Ψ [x (t)] / \sqrt{Ψ [x^{'} (t)]}

(23)

ω (t) \approx \sqrt{Ψ [x^{'} (t)] / Ψ [x (t)]}

(24)

Therefore, this method can be used to demodulate and analyze the signal.

In this paper, NMD method and Teager energy operator demodulation method are combined to perform NMD decomposition for composite fault signals, and then Teager energy operator demodulation for decomposed signals. At the same time, it is compared with EMD and EEMD. The specific fault diagnosis process is shown in Figure 2.

Figure 2.

Flow chart of composite fault diagnosis.

Simulation analysis

This paper compares NMD with EMD and EEMD to verify the advantages of NMD. In the paper, the sum of two Am and FM signals and one random signal is selected for simulation analysis. The simulated signals and components are as follows:

\begin{matrix} x (t) = x_{1} (t) + x_{2} (t) + 0.1 randn \\ x_{1} (t) = (1.5 + 0.5 \sin 26 π t) * \cos 600 π t \\ x_{2} (t) = (1 + \sin 20 π t) * \cos 200 π t \end{matrix}}

(25)

Among them, $randn$ is random noise of N [0,1], the sampling time is 1 s, and the sampling frequency is 1024 Hz. The simulation signal and each component is shown in the Figure 3:

Figure 3.

Simulated signal and two component time-domain waveforms.

Next, the NMD, EMD, and EEMD decompositions are performed on the simulation signals. Their decomposition results are shown in Figures 4 to 6 respectively.

Figure 4.

NMD simulation signal decomposition results.

Figure 5.

EMD simulation signal decomposition results.

Figure 6.

EEMD simulation signal decomposition results.

It can be seen from the decomposition results that NMD only decomposes the simulated signal into two NM components. There are no redundant other components. The decomposed NM1 and NM2 correspond to the simulation signals X1 and X2 respectively, and are basically consistent. However, after EMD decomposition, we can actually get more components, six of which are IMF1–IMF6, among which IMF4–IMF6 is basically meaningless. After EEMD decomposition, more components can be obtained, the first six are IMF1–IMF6, among which IMF4–IMF6 is basically a pseudo component, which has no practical significance. Therefore, from the perspective of the decomposition effect of simulation signals, the decomposition effect of NMD is indeed better than that of EMD and EEMD.

Next, the simulated components obtained from NMD, EMD and EEMD decomposition are demodulated by Teager energy operator respectively, obtain the NM simulation component demodulation spectrum (as shown in Figure 7), the IMF simulation component demodulation spectrum after EMD decomposition (as shown in Figure 8) and the IMF simulation component demodulation spectrum after EEMD decomposition (as shown in Figure 9). From the demodulation results of the NM simulation component and the demodulation results of the IMF simulation component decomposed by EMD and EEMD, it can be seen that the NM simulation component can be extracted to 1 and 2 times the characteristic frequency of the simulated signal fault, respectively. The IMF1 simulation component of EMD decomposition can be extracted to 1 frequency and 2 frequency, but the IMF2 simulation component of EMD decomposition can only be extracted to 1 frequency. The IMF1 simulation component of EEMD decomposition can also extract 1 and 2 times of frequency, but the IMF2 simulation component of EEMD decomposition can only extract 1 times of frequency. Although the demodulation results of EMD and EEMD simulation signal components can also get fault characteristic frequency, the number of peak points is less, and the phenomenon of mode aliasing occurs in the demodulation spectrum of EMD and EEMD simulation components. Therefore, the demodulation results show that the decomposition effect of NMD is indeed better than EMD and EEMD, and it can effectively suppress the phenomenon of modal aliasing, the decomposition and demodulation results of the simulated signal illustrate the advantages of the NMD method. Therefore, this paper uses the decomposition method of NMD.

Figure 7.

Demodulation spectrum of NM simulation component: (a) demodulation spectrum of NM1 simulation component and (b) demodulation spectrum of NM2 simulation component.

Figure 8.

Demodulation spectrum of IMF simulation component decomposed by EMD: (a) demodulation spectrum of IMF1 simulation component decomposed by EMD and (b) demodulation spectrum of IMF2 simulation component decomposed by EMD.

Figure 9.

Demodulation spectrum of IMF simulation component decomposed by EEMD: (a) demodulation spectrum of IMF1 simulation component decomposed by EEMD and(b) demodulation spectrum of IMF1 simulation component decomposed by EEMD.

Experimental analysis of gearbox composite fault diagnosis

NMD, EMD, and EEMD methods are compared and applied to the complex fault diagnosis of gear-gear transmission system respectively, so as to verify the effectiveness of this method in practical application. The equipment model used in the experiment is HFXC-I. The experimental equipment is shown in Figure 10. By manual setting, it is set up as a composite failure of pitting of big gear and wear of small gear. The acceleration sensor is used in the experiment, and the experiment is carried out under the load of 0.2 A, the sampling time is 2 S, the sampling frequency is 5120 Hz, and the measured speed is n = 825 r/min. The number of teeth of the small gear and the large gear is $Z_{1} = 55$ and $Z_{2} = 75$ , The pinion shaft is the driving shaft, By calculation, the fault characteristic frequency of the large gear is $f_{0}$ =10.08 Hz, The fault characteristic frequency of the pinion is $f_{r}$ =13.75 Hz. The calculation formula of fault characteristic frequency is shown in equations (26) and (27).

f_{r} = \frac{n}{60}

(26)

f_{0} = \frac{f_{r} * Z_{1}}{Z_{2}}

(27)

Figure 10.

Experimental equipment diagram.

The original composite fault signals were decomposed by EMD, EEMD and NMD methods, and the first two signals obtained by the decomposition of EMD and EEMD were taken, namely, IMF1 and IMF2, The results are shown in Figures 11 and 12, NMD decomposition results in two NM components as shown in Figure 13. Then the signals IMF1 and IMF2 obtained from the decomposition of EMD and EEMD were demodulated by Teager energy operator respectively and the demodulation results of the EMD decomposed signal is shown in Figure 14, and the demodulation results of the EEMD decomposed signal is shown in Figure 15. Then the two NM components obtained by NMD decomposition are demodulated by Teager energy operator, and the demodulation results is shown in Figure 16.

Figure 11.

Time domain diagram of IMF components decomposed by EMD: (a) time domain diagram of IMF1 components decomposed by EMD and (b) time domain diagram of IMF2 components decomposed by EMD.

Figure 12.

Time domain diagram of IMF components decomposed by EEMD: (a) time domain diagram of IMF1 components decomposed by EEMD and (b) time domain diagram of IMF2 components decomposed by EEMD.

Figure 13.

NM component time domain waveform: (a) NM1 component time domain waveform and (b) NM2 component time domain waveform.

Figure 14.

Demodulation spectrum of IMF components decomposed by EMD: (a) demodulation spectrum of IMF1 component of EMD decomposition and (b) demodulation spectrum of IMF2 component of EMD decomposition.

Figure 15.

Demodulated spectrum of IMF component of EEMD decomposition: (a) demodulation spectrum of IMF1 component of EEMD decomposition and (b) demodulation spectrum of IMF1 component of EEMD decomposition.

Figure 16.

NM component demodulation spectrum: (a) NM1 component demodulation spectrum and (b) NM2 component demodulation diagram

As shown in Figure 14(a), the demodulation results of the IMF1 component after EMD decomposition only extracted the 1-fold, 2-fold, and 3-fold frequency of the characteristic frequency of the large gear. As shown in Figure 15(a), the demodulation spectrum of IMF1 component decomposed by EEMD extracts the 1-fold frequency, 2-fold frequency, and 3-fold frequency of large gear fault feature frequency, and the 2-fold frequency of small gear fault feature frequency. As shown in Figure 16(a), the NM1 demodulation spectrum only extracted the 1-fold, 2-fold, and 3-fold frequency of the characteristic frequency of the large gear. The demodulation results show that EMD and NMD can extract the 1-fold, 2-fold, and 3-fold frequency of the fault feature frequency of large gear. Although EEMD can also extract the 1-fold, 2-fold, and 3-fold frequency of the fault feature frequency of large gear, it also extracts the 2-fold frequency of the fault feature frequency of small gear, which shows that EEMD cannot extract the fault feature frequency of large gear alone. As shown in Figure 14(b), the demodulation results of imf2 component after EMD decomposition extract the 1-fold and 2-fold frequency of pinion fault feature frequency, and extract the 1-fold, 2-fold, and 3-fold frequency of large gear fault feature frequency. As shown in Figure 15(b), the demodulation spectrum of IMF2 component after EEMD decomposition only extracted the 1-fold, 2-fold, and 3-fold the fault characteristic frequency of pinion. As shown in Figure 16(b), the NM2 demodulation spectrum only extracted the 1-fold, 2-fold, and 3-fold frequency of the characteristic frequency of the pinion. The demodulation results show that the EEMD and NMD methods can extract the fault characteristic frequency of the pinion. The effect is better than the EMD method.

Therefore, through the comparative analysis of the demodulation results, we can know that the EMD method can extract the fault feature frequency of the large gear separately, but it can’t extract the fault feature frequency of the small gear separately. The EEMD method can extract the fault characteristic frequency of the small gear, but it cannot extract the fault characteristic frequency of the large gear. It shows that the EMD and EEMD methods cannot effectively extract the fault characteristic frequency of the large and small gears separately. NMD method can not only extract the fault feature frequency of small gear, but also the fault feature frequency of big gear. The superiority of NMD method is illustrated. The NMD method can separate the fault characteristic frequency of the gears separately. It shows that the NMD method is better than the EMD and EEMD methods. Therefore, NMD method is suitable for composite fault diagnosis of gear transmission system, and the effect is good.

Conclusion

The structure of gear transmission system is complex, and it is easy to be affected by noise and other external factors, so it is not easy to extract fault feature frequency and identify fault type. Therefore, this paper proposes a method of Composite Fault Feature Extraction Based on the combination of nonlinear mode decomposition (NMD) and Teager energy operator.

Through comparative analysis with EMD and EEMD methods, it can be seen that the composite fault feature extraction method combining the nonlinear mode decomposition (NMD) and Teager energy operator proposed in this paper effectively suppresses the modal aliasing phenomenon.

Through the analysis of simulation signals and experimental data, we can see that the method proposed in this paper can effectively extract the composite fault characteristics and accurately identify the fault type. The validity of the proposed method is verified.

The proposed method has strong noise robustness and can effectively eliminate modal aliasing. This paper can effectively identify and extract the composite fault features of gear pitting wear. However, the correctness and accuracy of composite fault diagnosis of other mechanical equipment need to be further explored.

Footnotes

Acknowledgements

The authors gratefully acknowledge the support of project funded by the Science and Technology Research Projects of Education Department of Liaoning Province of China (LG201921), Liaoning BaiQianWan Talents Program (2020921031), Natural Science Foundation of Liaoning Province of China (2020-MS-216), State Key Laboratory of Mechanical Transmissions (SKLMT-KFKT-201605), and China Postdoctoral Science Foundation (2017M610496).

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Jingyue Wang

References

Singh

Parey

. Gearbox fault diagnosis under fluctuating load conditions with independent angular re-sampling technique, continuous wavelet transform and multilayer perceptron neural network. IET Sci Meas Technol 2017; 11(2): 220–225.

Yang

Shen

Wang

. Multi-Fault diagnosis of gearbox based on improved multipoint optimal minimum entropy deconvolution. Entropy 2018; 20(8): 611.

Zhao

Zhu

, et al. Application of local mean decomposition in composite fault diagnosis of rotating machinery. Noise Vib Control 2012; 5: 144–149.

Zhang

. Energy operator demodulating of optimal resonance components for the composite faults diagnosis of gearboxes. Meas Sci Technol 2015; 26(11): 115003.

Wang

Cai

, et al. Application of an improved ensemble local mean decomposition method for gearbox composite fault diagnosis. Complexity 2019; 2019(7): 1–17.

Isham

Leong

Lim

, et al. Intelligent wind turbine gearbox diagnosis using VMDEA and ELM. Wind Energy 2019; 22(6): 813–833.

Cheng

Huang

, et al. A gearbox compound fault diagnosis method based on SWD-AVDIF. Noise Vib Control 2019; 39(01): 171–176.

Yang

Sheng

Jing

, et al. Gear fault diagnosis based on EEMD and Choi Williams distribution. Mech Trans 2019; 43(04): 120–124.

Yang

Zeng

Cheng

. Research on the scale decomposition method of local features and its component criteria. China Mech Eng 2013; 02: 195–201.

10.

Liang

. Research on windowed Fourier transform harmonic detection algorithm and its interpolation improvement. J Electr Power Sci Tech 2007; 03: 40–44.

11.

Wang

Zhang

. Fault diagnosis of rolling bearings based on wavelet packet transform. China Mech Eng 2012; 23(3): 295–298.

12.

Huang

Shen

Long

, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc A 1998; 454(1971): 903–995.

13.

Cheng

Tang

, et al. Local rub-impact fault dianosis of the rotor systems based on EMD. Mech Mach Theory 2009; 44: 784–791.

14.

Yang

Suh

. Interpretation of crack-induced rotor non-linear response using instantaneous frequency. Mech Syst Signal Process 2004; 18(3): 491–513.

15.

Huang

. Ensemble empirical mode decomposition: a noise-assisted data analysis method. Adv Adapt Data Anal 2009; 01(01): 1–41.

16.

Iatsenko

Mcclintock

Ste-Fanovska

. Nonlinear mode decomposition: a nosise-roubst adaptive decomposition method. Phys Rev E 2015; 9(3): 032916.

17.

Feng

Qin

. Fault diagnosis of planetary gearbox based on Hilbert vibration decomposition and high-order energy operator. Vib Shock 2016; 35(5): 47–54.

18.

Gengel

Pikovsky

. Phase demodulation with iterative Hilbert transform embeddings. Signal Process 2019; 165: 115–127.

19.

Shi

Ding

, et al. An improved EMD method and its application in fault diagnosis of rotating machinery. J Vib Shock 2013; 32(4): 185–190.

20.

Liu

Chen

. Energy operator demodulation method and its application in mechanical signal demodulation. J Mech Eng 1998; 34(5): 85–90.

21.

Wang

Feng

Hao

. Research on fault diagnosis of rolling bearings based on teager energy operator. Vib Shock 2012; 31(2): 1–5.

22.

Iatsenko

Mcclintock

Stefanovska

. Extraction of instantaneous frequencies from ridges in time–frequency representations of signals. Signal Process 2016; 125: 290–303.

23.

Zhang

Tang

Deng

. Multi-component signal instantaneous parameter estimation and application based on wavelet ridge. J Mech Eng 2014; 20(10): 1–8.