A New Method of Two-stage Planetary Gearbox Fault Detection Based on Multi-Sensor Information Fusion

Abstract

Due to their high transmission ratio, high load carrying capacity and small size, planetary gears are widely used in the transmission systems of wind turbines. The planetary gearbox is the core of the transmission system of a wind turbine, but because of its special structure and complex internal and external excitation, the vibration signal spectrum shows strong nonlinearity, asymmetry and time variation, which brings great trouble to planetary gear fault diagnosis. The traditional time-frequency analysis technology is insufficient in the condition monitoring and fault diagnosis of wind turbines. For this reason, we propose a new method of planetary gearbox fault diagnosis based on Compressive sensing, Two-dimensional variational mode decomposition (2D-VMD) and full-vector spectrum technology. Firstly, the nonlinear reconstruction and noise reduction of the signal is carried out by using compressed sensing, and then the signal with multiple degrees of freedom is adaptively decomposed into multiple sets of characteristic scale components by using 2D-VMD. Then, Rényi entropy is used as the optimization index of 2D-VMD analysis performance to extract the effective target intrinsic mode function (IMF) component, reconstruct the dynamics signal in the planetary gearbox, and improve the signal-to-noise ratio. Then, using the full-vector spectrum technique, the homologous information collected by numerous sensors is data layer fused in the spatial domain and the time domain to increase the comprehensiveness and certainty of the fault information. Finally, the Teager–Kaiser energy operator is used to demodulate the potential low-frequency dynamics frequency characteristics from the high-frequency domain and detect the fault characteristic frequency. Furthermore, the correctness and validity of the method are verified by the fault test signal of the planetary gearbox.

1. Introduction

In recent years, environmental pollution has aroused the attention of various countries, and clean energy such as wind energy and solar energy has become the most important renewable energy in the world. With the expansion of the installation and operation scale of global wind turbines, the maintenance of wind turbines has become one of the most important research topics in the field of wind power. The wind turbine unit is in a complex environment with strong random wind load for a long time in operation, and the internal and external excitation is complex, which easily causes the failure of key parts of the transmission system, and then leads to the shutdown of the whole unit. Therefore, research on the structural health monitoring (SHM) of key components of wind turbines plays an important role in ensuring the safe and stable operation of the unit [1]. Structural health monitoring (SHM) can effectively avoid early failures of the wind turbine drive system and reduce maintenance and repair costs [2,3].

The planetary gearbox is the core of the wind turbine drive system. It has the advantages of compact structure, high transmission ratio, and small volume. It is mainly composed of planetary gears, planet carriers, sun gears, ring gears, boxes, and bearings. The planetary gearbox has several groups of gears engaged at the same time during transmission. Compared with the fixed shaft gear transmission, the planetary gearbox has special structure and dynamic characteristics, which also leads to its complex dynamics response, thus increasing the spectral complexity of the vibration signal [4,5,6]. The frequency spectrum of vibration signals of a planetary gearbox has typical asymmetry and nonlinearity due to the rotation of planetary gears [7,8,9]. The literature [10,11] summarizes the research methods of wind turbine fault detection proposed in recent years. The research shows that the traditional time-frequency analysis technology has some shortcomings in the condition monitoring and fault detection of wind turbines.

In view of these shortcomings, in recent years, researchers in the field of rotating machinery fault detection have proposed many new methods, such as: Feng et al. [4] proposed a joint amplitude frequency demodulation analysis method based on ensemble empirical mode decomposition (EEMD) and an energy separation algorithm. Zhang et al. [12] combined a blind deconvolution algorithm and a vibration model to detect seed faults in a helicopter planetary gearbox. Lei et al. [13,14] proposed four kinds of planetary gearbox fault detection parameters: root mean square of the filter signal, normalization of the positive amplitude of the difference spectrum between an unknown signal and a healthy signal, carrier series based on difference spectrum, and cumulative amplitude based on energy ratio. The key task of planetary gearbox fault detection is to determine the gear characteristic frequency [15]. Feng et al. [16] used the adaptive optimal kernel method to extract the frequency components and time-varying characteristics of the non-stationary signals of the planetary gearbox. Cheng et al. [17] proposed a planetary gear fault detection method based on ensemble empirical mode decomposition entropy feature fusion. Li et al. [18] proposed a fault detection method for the planetary gearbox based on improved multi-scale symbolic dynamic entropy and minimum-redundancy maximum-relevancy.

Most of the above methods are for the case of measuring a single-stage planetary transmission system with a single sensor. However, planetary gear transmission systems generally consist of multi-stage planetary transmissions. At the same time, multiple sets of planetary gears and sun gears along with planetary gears and ring gears are meshed simultaneously. The vibrations caused by these meshing pairs are mutually coupled, superimposed, and neutralized, which greatly affects the decoupling of fault characteristics, increases the complexity of the frequency spectrum, and further increases the difficulty of fault detection [19]. Therefore, it is necessary to study the method specially used to detect faults in a planetary gearbox with multi-stage planetary rows. According to the characteristics of the vibration signal of the planetary gear transmission system, using a multi-sensor data fusion method is a feasible method. Multi-sensor fusion has always been one of the key areas of data analysis [20], and it is an effective method to solve the fault detection problem of mechanical systems [21]. Multi-sensor data fusion enhances the effectiveness and success rate of fault monitoring by providing complementary information [22]. The information collected by multiple sensors simultaneously can overcome the problems of uncertainty and low accuracy in fault recognition due to the single-sided information of a single sensor, and also overcome the complexity of environmental noise and external complicated excitation, thereby improving the robustness and reliability of fault detection [23]. Multi-sensor information fusion can be divided into three grades [24,25], in which data level fusion directly and synchronously processes multiple sensors’ data from the tested object. This level of data fusion contains the most direct and reliable information of the measured object and can provide the most accurate troubleshooting results [26]. However, in the research on fault monitoring and detection of rotating machinery based on multiple sensors, data-level fusion methods mostly rely on artificial intelligence methods such as neural networks and support vector machines [27,28,29], and there are few research works on multi-sensor information fusion in the field of data layers [30,31].

The fault detection method based on multi-sensor data fusion needs to solve three problems in practical application [32]; these are: (1) how to remove noise effectively, (2) how to improve the fault feature, that is, the resolution of the signal, and (3) how to fuse the data of multiple sensors scientifically and effectively on the data layer. Next, this manuscript elaborates on the above three issues.

Empirical Mode Decomposition (EMD) and Local mean decomposition (LMD) are classical signal detecting methods proposed in recent years. However, research shows that the EMD and LMD methods also have the following shortcomings in signal detection and feature extraction: (1) The robustness of anti-noise is weak [33]; (2) the inefficient ability of decoupling frequency components [34]; (3) multivariate data cannot be processed at the same time; and (4) the mode mixing problem [35]. These shortcomings greatly affect the accuracy of signal processing [36].

In 2014, Dragomiretskiy and Zosso [37] proposed a new multi-resolution modal decomposition method called variational mode decomposition (VMD). Since it was published, it has been widely considered by scholars in various fields, and the method has been applied to medical treatment [38], image processing [39], financial analysis [40], electrocardiogram (ECG) [41] and other fields. The results show that VMD has an accurate mathematical model and theoretical basis in theory, and is superior to EMD and LMD [42] in signal separation and noise robustness. In the aspect of rotating machinery fault detection, researchers have carried out in-depth research. Zhang et al. [43] used VMD to detect defect signals in different positions of rolling bearings in multistage centrifugal pumps. Li et al. [44] proposed a new rolling bearing fault detection method based on variational mode decomposition and improved the Kernel Extremum Learning Machine. Li et al. [45] proposed an independent VMD method based on correlation analysis to adaptively extract weak fault features and compound fault features of wheelset bearings. An and Tang [46] proposed a bearing fault detection method based on energy distribution variational mode decomposition. However, the above research is still based on single-sensor information collection. Dragomiretskiy extended VMD to two-dimensional areas and proposed two-dimensional variational mode decomposition [47] in 2014, which makes simultaneous analysis and noise reduction of multivariate data possible.

The application of multiple sensors inevitably brings another problem, that is, the synchronous processing of massive monitoring data. The on-line condition monitoring system of the wind turbine functions by using multiple sensors installed in different key positions of the unit to collect real-time vibration acceleration, velocity, stress and strain and other parameters synchronously. Long-cycle real-time acquisition will generate a huge amount of data [48]. The huge amounts of data used for data transmission, storage, and real-time processing pose significant challenges [49]. At the same time, big data brings problems such as noise accumulation [50]. Therefore, research on multi-sensor data noise reduction methods has attracted researchers in various fields. Compressed sensing (CS) is a new signal acquisition technology, also known as compressed sampling [51]. As a new sampling theory, it develops the sparse characteristics of the signal, obtains discrete samples of the signal by random sampling, and then reconstructs the signal perfectly by a nonlinear reconstruction algorithm [52]. The recovered signal is sparse and robust at the same time [53]. Compressed sensing technology is widely used in signal denoising and feature extraction in various professional fields. Du et al. [54] applied compressed sensing theory to wireless sensor networks, and proposed a method of data sampling compression and reconstruction based on distributed compressed sensing theory. Guo [55] applied compressible sensing theory to the data recovery problem of practical sensor networks. Chen et al. [56] proposed a new sparse extraction of impulses by an adaptive dictionary scheme to extract impulse components.

Based on the above analysis, in order to solve the problem of planetary gearbox fault detection with multi-stage planetary rows, this paper proposes a new method of planetary gear fault detection with multi-sensor compressed sensing, 2D-VMD and full-vector spectrum theory. The paper is arranged as follows: in the preface, we review the recently proposed planetary gear fault detection method, and explain the necessity of studying the planetary gearbox fault detection method that is specially used to detect faults in a planetary gearbox with multi-stage planetary rows. In Section 2, we introduce compressed sensing, 2D-VMD, Rényi entropy and full-vector spectrum theory, and study the sparsity and noise robustness of compressed sensing in detail. A new method of planetary gear fault detection based on multi-sensor compressed sensing, 2D-VMD and full-vector spectrum is proposed. In Section 3, the synchronous acquisition of multivariate planetary gear tooth failure test data is carried out by using the planetary gear fault detection test-bed, and the spectrum characteristics of planetary gear tooth failure are studied in detail by using the proposed method. The fault features are extracted accurately, and the correctness and effectiveness of the method are proved. The conclusions are presented in Section 4.

2. Theory of Method

The application of multiple sensors brings the problem of massive data, using compressed sensing technology to collect information at a sampling frequency far lower than Nyquist sampling frequency. Compared with the data before processing, the nonlinear reconstruction algorithm perfectly reconstructs the signal. The processed signal has both sparsity and robustness, avoiding the redundancy of the collected information. The algorithm is as follows:

Suppose a one-dimensional discrete time signal $X$ [50]:

$$X=(x_1,x_2\cdots x_N{)}^T$$

(1)

The sparse basis matrix $\psi$ is used to calculate the sparse time signal $X$:

$$\psi ={[{\phi }_1,{\phi }_2\cdots {\phi }_N]}^T$$

(2)

$$X={\displaystyle \sum {\phi }_i}{s}_i$$

(3)

where ${\phi }_i$ is the i-th of N orthogonal basis vectors in the space, and the corresponding $s_i$ is the coefficient of the corresponding orthogonal basis vector.

Select the measurement matrix $\Phi$ and calculate the observation value $y$:

$$\Phi =[{\varphi }_1,{\varphi }_2\dots {\varphi }_N]$$

(4)

$$y=\Phi X$$

(5)

where ${\varphi }_1,{\varphi }_2\dots {\varphi }_N$ is the column vector of the measurement matrix $\Phi$.

The sensing matrix $\Theta$ is as follows:

$$\Theta =\Phi \psi$$

(6)

If there is a constant $\delta$, so that for any K sparse vector $x({\Vert x\Vert }_0\le K)$, it satisfies the restricted equidistant characteristic, and then the sparse signal can be effectively and completely recovered [57]:

$$(1-\delta ){\Vert x\Vert }_2^2\le {\Vert \Phi x\Vert }_2^2\le (1+\delta ){\Vert x\Vert }_2^2$$

(7)

In each iteration, the column vector having the largest projection coefficient in the sensing matrix is selected by the residual to form the recovery matrix $\stackrel{\sim }{\Phi }$ [58], and the column vector having the largest projection coefficient in the sensing matrix must be zeroed each time it is updated:

$${\varpi }_j=\left|\langle {\gamma }_{i-1},{\theta }_i\rangle \right|$$

(8)

$${\lambda }_i=\arg \max {\varpi }_j$$

(9)

where ${\gamma }_{i-1}$ is the residual of the i-1th iteration, ${\theta }_i$ is the column vector of the sensing matrix $\Theta$, ${\lambda }_i$ is the index of the maximum projection coefficient b, ${\lambda }_i$ is the indexes of the maximum value of the projection coefficient ${\varpi }_j$, and $\stackrel{\sim }{\Phi }$ is the set of column vectors having the largest projection coefficient in each iterative.

Solve the Moore–Penrose pseudo-inverse of the recovery matrix, and then based on Modified Cholesky decomposition, obtain the inverse $Q^{-1}$ of the symmetric matrix $Q$:

$$Q^{-1}={(L^{-1})}^T{D}^{-1}{L}^{-1}$$

(10)

The reconstructed signal ${\stackrel{\sim }{x}}_i$ obtained after each iteration is:

$${\stackrel{\sim }{x}}_i=Q^{-1}{\varpi }_j$$

(11)

The residuals that need to be iterated next time are updated according to the recovery matrix $\stackrel{\sim }{\Phi }$ and the reconstruction vector ${\stackrel{\sim }{x}}_i$ in each iteration until the final requirement is reached, and the final reconstructed signal $\stackrel{\sim }{x}$ is obtained:

$${\gamma }_i=y-\stackrel{\sim }{\Phi }{\stackrel{\sim }{x}}_i$$

(12)

We hope to find a method that can perform a frequency division operation on signals with multiple degrees of freedom at the same time, in order to achieve the filtering effect, and can maintain the coupling relationship between the original signals to the greatest extent. Therefore, we use the recently proposed 2D-VMD technique to adaptively decompose signals with three degrees of freedom into multiple sets of feature scale components. Then, Rényi entropy is used as the optimization index of 2D-VMD analysis performance to measure the correlation and coupling between the characteristic scale components of various degrees of freedom in different states. The signal is reconstructed using the extracted effective target IMF component to improve the signal resolution and signal-to-noise ratio.

The essence of 2D-VMD is a non-recursive and completely adaptive decomposition of the original function to obtain a specified number of modal components. Let the expression of the synthesized two-degree-of-freedom directional vibration signal sequence be f(x), which is equivalent to solving the following optimization problem [47]:

$$\underset{m_K:{ℝ}^2\to ℝ,{\omega }_K\in {ℝ}^2}{\min }\{{{\displaystyle \sum _K{\beta }_K\Vert \nabla \left[m_{AS,K}\left(x\right)e^{-j\langle {\omega }_K,x\rangle }\right]\Vert }}_2^2\}s.t.\forall x\in {ℝ}^2:{\displaystyle \sum _K{m}_K(x)=f(x)}$$

(13)

where $m_{AS,K}$ is the analytical signal obtained from the modal component $m_K$ according to ${\omega }_K$; ${\beta }_K$ is the weight coefficient of the mode; ${\Vert \cdot \Vert }_2^2={{\displaystyle {\int }_{{ℝ}^2}\left|\cdot \left(X\right)\right|}}^{2}dx$ can be understood as a norm in a two-dimensional case.

Lagrange multiplier and quadratic penalty function are introduced to further constrain the above problems:

$$\begin{array}{c}L(\{m_K\},\{{\omega }_K\},\lambda )={{\displaystyle \sum _K{\beta }_K\Vert \nabla [m_{AS,K}(x)e^{-j\langle {\omega }_K,x\rangle }]\Vert }}_2^2\\ +{\Vert f(x)-{\displaystyle \sum _K{m}_K(x)}\Vert }_2^2+\langle \lambda (x),f(x)-{\displaystyle \sum _K{m}_K(x)}\rangle \end{array}$$

(14)

Among them, $\lambda :{ℝ}^2\to ℝ$ is a Lagrangian multiplier.

The problem is further transformed into an unconstrained saddle point problem by the augmented Lagrangian equation defined above.

$$m_K:{ℝ}^n\to \stackrel{\min }{ℝ},{\omega }_K\in {ℝ}^n\lambda :\stackrel{\max }{{ℝ}^n}\to ℝ\begin{array}{c}L(\{m_K\},\{{\omega }_k\},\lambda )\\ \end{array}$$

(15)

The Lagrange multipliers are combined with the quadratic penalty terms:

$$\begin{array}{c}L(\{m_K\},\{{\omega }_K\},\lambda )={\displaystyle \sum _K{\beta }_K}{\Vert \nabla [m_{AS,K}(x)e^{-j\langle {\omega }_K,x\rangle }]\Vert }_2^2\\ +{\Vert f(x)-{\displaystyle \sum _K{m}_K(x)+\frac{\lambda (x)}{2}}\Vert }_2^2-{\Vert \frac{\lambda (x)}{2}\Vert }_2^2\end{array}$$

(16)

In each iteration, the modal components, center frequency and Lagrangian multipliers in the updated state are taken as the saddle points of the augmented Lagrangian equation to solve the optimal data [59,60,61]:

$$m_K^{p+1}\leftarrow \begin{array}{c}\arg \min \\ m_K:{ℝ}^n\to ℝ\end{array}L(\{m_{i<K}^{p+1}\},m_K,\{m_{i>K}^p\},\{{\omega }_i^p\},{\lambda }^p)$$

(17)

$${\omega }_K^{p+1}\leftarrow \begin{array}{c}\arg \min \\ {\omega }_K\in {ℝ}^n\end{array}L(\{m_i^{p+1}\},\{m_{i<K}^{p+1}\},{\omega }_K,\{{\omega }_{i>K}^p\},{\lambda }^p)$$

(18)

$${\lambda }^{p+1}\leftarrow {\lambda }^p+\tau (f-{\displaystyle \sum _K{m}_K^{p+1}})$$

(19)

where $\tau \ge 0$.

After eliminating the first variation of the modal component and reaching the optimal condition, solve the modal component $\stackrel{\wedge }{m_K}$ to obtain the Wiener filter update:

$${\stackrel{\wedge }{m}}_K^{p+1}(\omega )=(\stackrel{\wedge }{f(\omega )}-{\displaystyle \sum _{i\ne K}\stackrel{\wedge }{m_i}}(\omega )+\frac{\stackrel{\wedge }{\lambda }(\omega )}{2})\frac{1}{1+2{\beta }_K{\left|\omega -{\omega }_K\right|}^2},\forall \omega \in {\Omega }_K$$

(20)

where ${\Omega }_K$ is the frequency domain half space.

The center frequency ${\omega }_K$ in the update is:

$${\omega }_K^{p+1}=\underset{{\omega }_K\in {ℝ}^2}{\arg \min }\{{\beta }_K{\Vert \nabla [m_{AS,K}(x)e^{-j\langle {\omega }_K,x\rangle }]\Vert }_2^2$$

(21)

The maximum Lagrange multiplier $\lambda$ is:

$${\lambda }^{p+1}(x)={\lambda }^p(x)+\tau (f(x)-{\displaystyle \sum _K{m}_K^{p+1}(x))}$$

(22)

The above processing reduces the influence of some external noise on the fault feature information to some extent, but the energy of the fault feature information will inevitably leak due to the filtering processing. However, we can use the IMF component with the smallest energy leakage to complete the fault diagnosis of many typical nonlinear dynamic systems. Therefore, screening the IMF component with the smallest amount of energy leakage is the key to optimizing 2D-VMD technology. As an effective tool to quantify the correlation and coupling between systems, Rényi entropy is well applied in fault diagnosis and other fields. The more characteristic fault information contained in the IMF component in the fault state, the smaller the value of Rényi entropy, and the smaller the correlation and coupling degree between the fault and the healthy system. Therefore, the extraction of the effective target IMF component is completed by calculating the Rényi entropy between each IMF component in the fault state and the healthy planetary gearbox signal, thereby achieving the purpose of improving the signal-to-noise ratio of the signal and optimizing the targets of 2D-VMD. The Rényi entropy algorithm is as follows:

Uncorrelated vectors are orthogonalized using Schmidt orthogonalization.

$$|{\psi }_1\rangle ={\{\mathsf{\Upsilon }\}}_1$$

(23)

$$|{\psi }_k\rangle ={\{\mathsf{\Upsilon }\}}_k-\frac{\langle |{\psi }_k\rangle ,{\{\mathsf{\Upsilon }\}}_{k-1}\rangle }{\langle {\{\mathsf{\Upsilon }\}}_{k-1},{\{\mathsf{\Upsilon }\}}_{k-1}\rangle }{\{\mathsf{\Upsilon }\}}_{k-1}$$

(24)

In the formula, ${\{\mathsf{\Upsilon }\}}_k$ is a dynamics signal in the form of a column vector.

$$\rho ={\displaystyle \sum _k{p}_k|{\psi }_k\rangle }\langle {\psi }_k|$$

(25)

$${\displaystyle \sum _k{p}_k=1,0<p_k<1}$$

(26)

where $\rho$ is the density matrix; $p_k$ is the probability value of each power system (here,$p_1=p_2=1/2$).

Rényi Entropy is:

$$S_{\alpha }(\rho )=\frac{1}{1-\alpha }\log tr({\rho }^{\alpha }),\alpha \in {\mathbb{Z}}_{+} \text{and} \alpha \ne 1$$

(27)

where $\alpha$ is the scale parameter of the Rényi entropy.

The following is a full-vector synthesis of the vibration acceleration signals of the multi-stage planetary gearbox in both the x and y directions. The algorithm is as follows [62]:

Let the vibration signal of the rotor in two mutually perpendicular directions be {x_i}, {y_i}, i = 1, 2 …, N/2-1, complex sequence {z_i} = {x_i} + j{y_i}. The Fourier transform is:

$$\{Z_i\}=\{X_i\}+j\{Y_i\}$$

(28)

Available:

$$\{\begin{array}{c}{M}_{ci}=\frac{1}{2}[X_i+X_{N-i}]\\ M_{ai}=\frac{1}{2}[Y_i-Y_{N-i}]\end{array},\{\begin{array}{c}{N}_{ci}=\frac{1}{2}[Y_i+Y_{N-i}]\\ N_{ai}=\frac{1}{2}[X_i-X_{N-i}]\end{array}$$

(29)

In the above formula, i = 1, 2,…, N/2 − 1; $M_{ci}=M_i\cos {\varphi }_{xi}$; ${\varphi }_{xi}$ and $M_i$ are the phase angle and amplitude of x-direction signals, respectively. $M_{ai}=M_i\sin {\varphi }_{xi}$, $N_{xi}=N_i\cos {\varphi }_{yi}$, $N_{ai}=M_i\sin {\varphi }_{yi}$; ${\varphi }_{yi}$ and $N_{xi}$ are the phase angle and amplitude of y-direction signals, respectively. Ellipse parameters:

$$\{\begin{array}{c}{R}_{ci}=\frac{1}{2N}[\left|Z_i\right|+\left|Z_{N-i}\right|]\\ R_{si}=\frac{1}{2N}[\left|Z_i\right|-\left|Z_{N-i}\right|]\\ tg2{\beta }_i=\frac{Y_i{X}_{N-i}-X_i{Y}_{N-i}}{X_i{X}_{N-i}+Y_i{Y}_{N-i}}\\ tg{\varphi }_i=\frac{Y_i+Y_{N-i}}{X_i+X_{N-i}}\end{array},i=1,2\cdots ,N/2-1$$

(30)

In the formula, ${\beta }_i$ is the angle between the main vibration vector and the x-axis; $R_{ci}$ is the main vibration vector (i.e., the elliptical long axis), $R_{si}$ is the secondary vibration vector (i.e., the elliptical short axis); ${\varphi }_i$ is the initial phase angle of the ellipse.

We hope to demodulate the characteristic frequency of a low-frequency fault with higher accuracy by a better demodulation method, so as to effectively avoid the windowing effect brought by the traditional method. The energy operator has a good time resolution for the instantaneous change of any discrete-time signal, and can effectively detect the transient components in the signal. The detection signal has been decomposed into a certain number of single-component signals by 2D-VMD in the above, which lays the foundation for the use of energy operators.

The Teager–Kaiser energy operator can measure the instantaneous energy change of non-linear and unstable signals. For any signal h (t), the Teager–Kaiser energy operator Ψ is defined as [63,64,65]:

$$\psi \left[h\left(t\right)\right]={\left[\dot{h}\left(t\right)\right]}^2-h\left(t\right)\ddot{h}\left(t\right)$$

(31)

where $\dot{h}\left(t\right)$ and $\ddot{h}\left(t\right)$ are the first and second derivatives of the signal h(t), respectively. For discrete time signal h(t), the Teager–Kaiser energy operator is:

$$\psi \left[h(n)\right]={\left[h(n)\right]}^2-h(n-1)h(n+1)$$

(32)

The Teager–Kaiser energy operator can separate the instantaneous amplitude and instantaneous frequency of a single-component AM-FM signal, that is [66,67]:

$$\left|a(t)\right|=\frac{{\psi }_c\left[h(t)\right]}{\sqrt{{\psi }_c\left[\dot{h}(t)\right]}}$$

(33)

$$w_i(t)=\sqrt{\frac{{\psi }_c\left[\dot{h}(t)\right]}{{\psi }_c\left[h(t)\right]}}$$

(34)

The planetary gearbox has a special structure and complex internal and external excitation, which makes the vibration signal spectrum show strong nonlinear, asymmetric and time-varying characteristics. In this paper, a new fault detection method for the planetary gearbox is proposed by combining compressed sensing, 2D-VMD, full-vector spectrum technology and energy operator demodulation. The specific process is as follows:

1. Signal extraction: firstly, several vibration acceleration sensors are arranged in the key position of the planetary gearbox, and the state signals of the planetary gearbox are collected synchronously.

2. Signal preprocessing: The original signals collected by multiple sensors are preprocessed, and the signals that can be sparsely represented are subsampled with compressed sensing technology. Then, the atoms in the sensing matrix that have the greatest correlation with the residual are selected through greedy iteration to achieve signal vector approximation and complete signal deposition.

3. State information enhancement: the noise reduction signal is decomposed adaptively into multiple sets of feature scale components by the 2D-VMD technique, and Rényi entropy is used to measure the linear correlation of information as an optimization index of 2D-VMD performance, improve the signal-to-noise ratio, enhance the state information, and extract the fault features contained in the effective target IMF component.

4. Vibration signal fusion: using the full-vector spectrum technique, the homologous information collected by multiple sensors is data layer fused in the spatial domain and the time domain to increase the comprehensiveness and certainty of the fault information.

5. Fault detection: the envelope spectrum of amplitude modulation and frequency modulation signal is demodulated by the energy operator, and then the frequency characteristic of low-frequency signal is restored by Fourier transform.

The flow chart of the proposed method is shown in Figure 1.

3. Experimental Study

3.1. Introduction to the Experiment

Compared with fixed shaft gear drive and single-stage planetary gear drive systems, the dynamics response of the multi-stage planetary gear drive system is characterized by high complexity, high nonlinearity, complex internal and external excitation, and strong unpredictability of the operation environment. In order to verify the effectiveness of the proposed method in the multi-stage planetary gearbox fault detection, the multi-stage planetary gearbox fault diagnosis experimental platform provided by spectraquest company is used for experimental verification. The experimental platform mainly includes: a two-stage planetary gear transmission system, editable servo motor, magnetic powder brake, data acquisition system, etc. [17,18], as shown in Figure 2. The test bench can simulate typical faults such as broken teeth, tooth surface wear, pitting and missing teeth of key parts such as planetary gears, sun gears, ring gears, etc., and adopt a three-way acceleration sensor to collect x, y, and z in the state of broken gears. The acceleration signal is in three directions, the sensor arrangement is shown in Figure 3, the sampling frequency is 12,800 Hz, and the sampling time is 10 s [67,68,69]. The broken tooth fault is set at the second stage solar gear of the planetary gearbox, and the input frequency of the servo motor is 40 Hz. The gearbox parameters are shown in

Table 1. Planetary gearbox parameters.

	Gear	Sun Wheel	Planetary Gears (Number)	Gear Ring
Tooth Number	First stage	28	36(4)	100
	Second stage	20	40(3)	100

.

The local characteristic frequency algorithm of the planetary gearbox is as follows [70]:

$$\{\begin{array}{c}{f}_m=f_c{Z}_r=(f_s^{(r)}-f_c)Z_s\\ f_s=\frac{f_m}{Z_s}N\\ f_p=2\frac{f_m}{Z_P}\\ f_r=\frac{f_m}{Z_r}N\end{array}$$

(35)

where $Z_r$ and $Z_s$ are the number of teeth of the gear ring and sun wheel; $f_m$ is the meshing frequency; $f_s^{(r)}$ is the rotation frequency of the solar wheel; $f_c$ is the rotation frequency of the planetary frame; $f_s$ is the characteristic frequency of local faults of the solar wheel; $f_p$ is the local fault characteristic frequency of the planetary gear; $f_r$ is the characteristic frequency of local faults of the gear ring; $Z_s$ is the number of teeth of the solar gear; $N$ is the number of planetary wheels; $Z_P$ is the number of planetary gear teeth; $Z_r$ is the number of tooth rings. See the

Table 2. Characteristic frequency of the planetary gearbox.

Meshing Frequency (Hz)	Absolute Rotation Frequency (Hz)		Local Fault Characteristic Frequency (Hz)
	$\mathbf{Sun}\mathbf{Wheel}{\mathit{f}}_{\mathit{a}}$	$\mathbf{Planet}\mathbf{Carrier}{\mathit{f}}_{\mathit{c}}$	$\mathbf{Sun}\mathbf{Wheel}{\mathit{f}}_{\mathit{s}}$	$\mathbf{Planet}\mathbf{Gear}{\mathit{f}}_{\mathit{p}}$	$\mathbf{Gear}\mathbf{Ring}{\mathit{f}}_{\mathit{r}}$
145	6.67	1.45	20.84	8.11	5.83

below for specific parameters.

3.2. Experimental Analysis

The time-domain diagram and spectrogram in the x-direction of the planetary gearbox in a healthy state are shown in Figure 4. The time-domain diagram and spectrum in the x-direction of the planetary gearbox in the broken state are shown in Figure 5.

After the vibration signal of the planetary gearbox in the state of broken teeth is processed by the compressed sensing technology, the time domain diagram and frequency spectrum diagram of the x-direction signals are as shown in Figure 6. Compared with the original frequency domain signal, the peak frequency and sideband of noise reduction signal are more concentrated in 1200–2000 Hz. It can be seen that the compressed sensing technology can improve the accuracy of the feature information range, strengthen the broken tooth state signal of the planetary gearbox in the high-frequency domain, effectively reduce the noise interference, and improve the signal-to-noise ratio of the experimental signal.

2D-VMD analysis was carried out on the x and y vibration signals processed by compressed sensing, and several IMF components were obtained adaptively. The spectrum of each IMF component is shown in Figure 7. In the full frequency domain, the Lagrangian multiplier and the quadratic penalty function limit the center frequency position and the width of the surrounding sidebands. The non-stationary features are also clearly picked up while effectively decomposing the signal. The center frequencies of the six IMF components are approximately: 120, 1950, 1318, 2600, 3713, and 6060 Hz. In the actual low signal-to-noise ratio condition, there is a phenomenon of superposition and coupling between vibration signals of different amplitudes, and rich harmonic components are embedded in the frequency spectrum, resulting in the presentation of fault spectrum characteristics in IMF being affected [35]. Therefore, it is necessary to further deal with the IMF components obtained by decomposition, that is, to use Rényi entropy to solve the above problems.

Rényi entropy reflects the nonlinear degree and coupling degree between different systems. Therefore, a method of extracting effective target IMF components based on Rényi entropy is proposed. Rényi entropy is used to be the optimization index of 2D-VMD performance by measuring the coupling degree of information, so as to extract the fault features contained in the effective target IMF component, solve the problem of the interference of mixed mode to the fault feature frequency, and further improve the signal-to-noise ratio. Rényi entropy can quantify the coupling degree between the fault system and the healthy system. The smaller the entropy is, the less the coupling degree between the two systems is, and the richer the fault information contained in the fault system is. The Rényi entropy with the parameter value $\alpha =4$ between the IMF component of each broken tooth state of the planetary gearbox and the healthy signal is calculated, respectively, and the calculated results are shown in Figure 8. In Figure 8a, IMF₁ > IMF₄ > IMF₂ > IMF₆ > IMF₃ > IMF₅; In Figure 8b, IMF₁ > IMF₄ > IMF₂ > IMF₆ > IMF₃ > IMF₅. As shown in the figure, the Rényi entropy between the IMF₁ component and healthy signal is the largest under the condition of planetary gear teeth broken, indicating that the coupling relationship between the IMF₁ information feature and normal signal is the largest, and the fault feature information is relatively the least. The Rényi entropy corresponding to IMF₃, IMF₅ and IMF₆ is small, indicating that there is a strong nonlinear relationship between the broken tooth IMF component and the healthy signal, and the fault information is more abundant. The components of IMF₃, IMF₅ and IMF₆ are in the high-frequency domain. After data fusion, the influence of noise can be reduced and the workload of data calculation can also be reduced. So, IMF₃, IMF₅ and IMF₆ are selected to extract and process the follow-up fault information, which provides the basis for the follow-up fault diagnosis.

Multi-channel homologous information gathering realizes the superposition and complementarity of the information of each degree of freedom, and fully reflects the integrity of the original information. Figure 9 shows the spectrum diagram of the broken tooth state of the planetary gearbox and the local magnification of part of the meshing frequency obtained by multi-channel acquisition and full-vector spectrum processing. As shown in Figure 9a, in the spectrum diagram, the sun gear in the broken tooth state will have a peak value at its engagement frequency and its frequency doubling, and will generate periodic impact signal (f_s = 20.84 Hz). Due to the influence of manufacturing error, there is a certain difference in the impact characteristics between the solar wheel and the three planetary gears. The local fault characteristic frequency of the solar wheel will become f_s/3 = 6.947 Hz and there will be a peak value at its frequency doubling. At the same time, based on the modulation effect of the transmission effect of the system on the vibration signal, the peak value will appear at the absolute rotation frequency f_c and its frequency doubling of the planetary frame. Based on the above factors, the peak value will appear at the frequency of $k{f}_m\pm m{f}_c\pm (m\pm n)f_s/3$ Hz in the frequency spectrum of the multi-stage planetary gearbox when the gear is broken. The side band is formed around the meshing frequency and its doubling frequency, and its amplitude and width are larger than the healthy signal [70]. As shown in Figure 9b, peak values of 1.422 × 10⁻³, 2.396 × 10⁻³, 8.519 × 10⁻⁴, 1.950 × 10⁻³ and 3.453 × 10⁻³ appeared at positions 9f_m + 7f_c – 5f_s = 1211 Hz, 9f_m – 3f_c – 3f_s = 1237 Hz, 9f_m + 2f_c + f_s/3 = 1314 Hz, 9f_m + f_c + 5f_s/3 = 1340 Hz, 9f_m + 3f_c + 4 = 1357 Hz and so on, respectively. In Figure 9c, peak values of 8.331 × 10⁻⁴, 9.986 × 10⁻⁴, 1.856 × 10⁻⁴, 4.919 × 10⁻⁴, 1.280 × 10⁻³ appeared at positions 11f_m + 3f_c – 2f_s = 1557 Hz, 11f_m – 2f_c – 2f_s/3 = 1577 Hz, 11f_m + f_c = 1597 Hz, 11f_m + f_c + 5f_s3 = 1630 Hz, 11f_m + f_c + 2f_s = 1637 Hz and so on, respectively. There are peaks at 9f_m + f_s/3 = 1311 Hz, 9f_m – 2f_c = 1301 Hz of Figure 9b and 11f_m – 2f_c = 1591 Hz, 11f_m + 4f_s/3 = 1622 Hz and so on of Figure 9c. However, most of the peaks are located at mf_c + (m+n)f_s/3. The above experimental results verify the feasibility of extracting effective target IMF components based on Rényi entropy, further optimize the decomposition performance of 2D-VMD, and prove the accuracy of the proposed method as well as the ability of fault detection in complex internal excitation environments.

In Figure 9, the partial fault frequency information of the broken tooth sun gear is analyzed. It can be found that the feature information is all distributed in a relatively high-frequency domain range. Because the characteristic frequency (f_s/3 = 6.947 Hz) of the local fault solar wheel, the absolute rotation frequency (f_c = 1.45 Hz) of the planetary rack, and the absolute rotation frequency (f_a = 6.67 Hz) of the solar wheel are too small, they will be affected by external complex excitation, internal excitation, noise and other signals, and the modulation will disappear in the low-frequency range, which has a great impact on fault diagnosis. Therefore, we use the Teager–Kaiser energy operator method to demodulate the potential low-frequency broken-tooth fault frequency modulated into the high-frequency region. The result is shown in Figure 10. It can be seen from Figure 10 that there are peaks at the local fault characteristic frequency of the sun wheel (nf_s/3), the absolute rotation frequency of the sun wheel (mf_a), and the absolute rotation frequency (kf_c) of the planet carrier (such as 5f_c = 7.5Hz, 8f_c = 12 Hz, 2f_a = 13.5 Hz, 5f_a = 33.5 Hz, 2f_s = 42 Hz, and 3f_s = 62 Hz), and the amplitude is higher than the healthy signal. There are a large number of side bands around the above frequency positions, and most of the peaks are located at the combination of planetary frame rotating frequency and local fault characteristic frequency of the solar wheel [70] (such as f_s − 3f_c = 16.5 Hz, f_s + f_c = 22.5 Hz and f_s + 3f_c = 23 Hz near 12f_c). The presence of 2f_s − 2f_c = 39 Hz, 4f_s/3 + 14f_c = 48 Hz and 2f_s + 5f_c = 49 Hz near 2f_s further illustrates the existence of the local fault of the solar wheel, and they also prove the correctness and practicability of the proposed method.

4. Conclusions

We propose a new method for planetary gearbox fault detection based on compressed sensing, 2D-VMD and full-vector spectrum theory. The conclusions are as follows:

1. In view of the complexity and low signal-to-noise ratio of the vibration signal extracted from the multi-stage planetary gearbox, compressed sensing is used to denoise the signal, which proves that compressed sensing has strong robustness and noise reduction ability for noise signals.

2. In order to realize the synchronous noise reduction of multi-degree-of-freedom signals, we propose a new idea of multi-sensor data fusion, which extends one-dimensional variational mode decomposition to multi-dimensional. The 2D-VMD method with non-recursive mode and strong mathematical theory is used to process the planetary gearbox fault signal, and Rényi entropy is used to be the optimization index of 2D-VMD performance by measuring the coupling degree of information, so as to extract the fault features contained in the effective target IMF component, and the feasibility of the method is proved by analyzing and showing the fault characteristic frequency of the solar wheel in the planetary gearbox. The proposed method realizes the synchronous processing of multivariate vibration signals and can efficiently and comprehensively extract the effective spectrum characteristics of planetary gearbox fault state signals.

3. Using the full-vector spectrum technique, the homologous information collected by multiple sensors is data layer fused in the spatial domain and the time domain to increase the comprehensiveness and certainty of the fault information.

4. When the multi-stage planetary gearbox is in the state of broken teeth, there are peak values at the meshing frequency, at its doubling frequency, and at $k{f}_m\pm m{f}_c\pm (m\pm n)f_s/3$, and the side frequency band is formed at the same time, while the width and amplitude are greater than the healthy signal. Potential low-frequency features need to be demodulated from the high-frequency domain using energy operators, and the energy operator has higher resolution than envelope demodulation. When the multi-stage planetary gearbox is in local fault, the peak values will appear at the frequency multiplication (nf_s/3) of the sun wheel’s local fault characteristic frequency, at the frequency multiplication (mf_c) of the planet carrier’s absolute rotation frequency, and at the frequency multiplication (kf_a) of the sun wheel’s absolute rotation frequency; the peak values appear at the combined position of the three and form the side band at the same time.