Rolling bearing compound fault diagnosis based on spatiotemporal intrinsic mode decomposition

2.1. Fast independent component analysis

Independent component analysis is a method to solve the problem of blind source separation. It separates signals based on signal independence and can obtain the separation matrix when the source signal is unknown. It assumes that the source signals are statistically independent, the source signals are non-Gaussian distribution, and the mixed matrix $B$ is orthogonal. $X=BS$ [20], where $X$ represents the mixed signal matrix, $B$ represents the mixed matrix, and $S$ represents the signal matrix after unmixing. If $S$ follows a Gaussian distribution, then a unique $S$ cannot be reduced. If $S$ is not Gaussian, then $S=WX$ can be calculated by calculating $W$, where $W$ represents the separation matrix. The central limit theorem states the conditions under which a large number of random variables approximate normal distribution [21, 22]. The central limit theorem is the main theoretical basis of the objective function [23]. By extension, it can be known that, compared with the source signal, the mixed signal received by the sensor is more gauss, and conversely, the non-gauss is weak. Therefore, the separation of mixed signals can be completed by performing non-Gaussian maximization operation on mixed signals. The independent component analysis algorithm [24] computes the $B$ and $S$ matrices using kurtosis or negative entropy by maximizing the non-Gaussian property of the source signal. Fast independent component analysis finds the direction of maximizing non-Gaussian measurement through projection pursuit [25].

The new signal is obtained by preprocessing the observed signal. $w$ is the row vector in the unmixing matrix $W$. The FastICA algorithm maximizes the non-Gaussian property of a vector $w^{T}z$ by a fixed point iterative algorithm. The most common evaluation function for non-Gaussianness is the approximation of negative entropy. The objective function of the FastICA algorithm is:

where $G$ is any non-quadratic function and $v$ is a Gaussian random variable whose mean value and unit error function is zero. $z$ means $z=B^{T}w$. Get the best value by maximizing in $J\left(w\right)$. When ${‖w‖}^2=1$, the negative entropy $J\left(w\right)$ is the highest, that $E\{G\left(w^{T}z\right)$ is the maximum. Through the Lagrange multiplier algorithm, the objective function of the fixed point algorithm is [26]:

where $\beta$ is a constant, so the optimization problem is equivalent to:

where the function $\mathrm{g}$ is the derivative of $G$. The solution to the equation 3 is similar to the Newtonian iteration, that is, leaving the left side of the equation unchanged:

The gradient of $F\left(w\right)$ is:

where $T$ for transpose. Since $S$ are independent sources, it can be assumed that the variance of $s_i$ is consistent. So the covariance of $S$ is equal to $I$ [27]:

The final iteration form of the fast ICA algorithm can be obtained by approximating the equation:

The weight vector $w$ is normalized after each iteration:

2.2. Nonlinear matching pursuit

Nonlinear matching pursuit (NMP) can be used to separate intrinsic mode function (IMF) from mixed signals instead of empirical mode decomposition. Assuming that the signal $x\left(t\right)$ can be decomposed into a finite number of IMFs, the following problems can be solved:

where $s_j$ is assumed to contain the following interwave frequency modulated signals. $\theta \left(t\right)$ can be known to be monotone by the invertible function theorem, and IMF $s\left(\theta \right)$ is expressed by the function of $\theta$:

If the $a\left(\theta \right)$ and ${\theta }^{\text{'}}\left(\theta \right)={\left.\frac{d\theta }{dt}\right|}_{\theta }$ are smoother than $\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)$, then $s\left(\theta \right)$ is defined as an interwave frequency modulated signal. That is to say, the $a\left(\theta \right)$ and $\theta \text{'}\left(\theta \right)$ are in the following set:

Among them $L_{\theta }=\left[\frac{\theta \left(t_1\right)-\theta \left(t_0\right)}{2\pi }\right]$, $\lambda =1/2$.

A signal that physically contains an interwave frequency modulation component can roughly correspond to the solution of the second-order differential equation:

where the $b\left(t\right)$ and $c\left(t\right)$ are the sufficiently smooth functions.

To sum up, the IMF dictionary with an interwave frequency modulation component is:

when the sum of $a_j{cos}_j$ is equal to the original signal $x$, the quantity $M$ is minimum. For noisy signals the equality is replaced by inequality${‖x-\sum _{j=1}^M{a}_j{cos}_j‖}_2\le \delta$.

To solve the problem of NMP method minimization, an alternate scheme is proposed, that is, fixed $\theta$ minimization $a$ and then fixed $a$ updating $\theta$. In the minimization calculation, projecting $a\left(\theta \right)$ onto the surface $V\left(\theta \right)$ is equivalent to applying a low-pass filter to the $\theta$ coordinates. The guess value of the initial phase function is generally $2\pi ft$.

2.3. Fast spectral kurtosis

Kurtosis is a dimensionless parameter index in the time domain which is very sensitive to the instantaneous characteristics of the signal. Kurtosis coefficient refers to the impact pulse generated at the defect of the working face when the working surface of rotating parts such as bearings is faulty. It is often used to detect the strength of the impact component in the vibration signal of rolling bearings [25]:

Spectral kurtosis (SK) is used to characterize the kurtosis value of each spectral line of the signal, so as to find the non-stationary components of the signal and their positions in the frequency domain [20]:

where $E\left(x\right)$ is the mathematical expectation; $H(t,f)$ is the complex envelope of the original signal at frequency $f$. In order to obtain the transient components in the signal, the kurtosis value of each frequency band should be calculated to find the frequency band with the maximum kurtosis. In order to reduce the calculation time of spectral kurtosis significantly, fast spectral kurtosis calculation is carried out by using the 1/3 – binary tree filter banks to realize the fast calculation of spectral kurtosis of each sub-band, and the kurtosis value is represented according to the depth of the color graph [28].

2.4. Spatiotemporal intrinsic mode decomposition methods

The spatiotemporal intrinsic mode decomposition method establishes an over-complete dictionary base through empirical mode decomposition as the constraint condition of nonlinear matching pursuit method, which can effectively separate the impact component from the vibration signal of hybrid bearing. The STIMD method integrates the idea of nonlinear matching pursuit, so that it can detect the frequency characteristic component of mixed signal near the initial phase function, so it can still effectively separate various mixed components in the case of high noise. The STIMD method combines the advantages of fast independent component analysis and nonlinear matching pursuit. It can separate the multi-dimensional signal blind source to obtain the component signal generated by each vibration source and analyze the component signal to get the fault characteristics of the rolling bearing vibration signal. By nonlinear modelling of complex multidimensional mixed signals and constructing the eigenmode function dictionary conforming to the characteristics of bearing vibration signals, the spatiotemporal intrinsic mode decomposition method can effectively extract the impact signal components of bearing vibration from the mixed signals. By analyzing the envelope spectrum of the component, the characteristic frequency of the fault can be seen and the fault type can be determined.

Fig. 1Flow chart of STIMD

It is particularly important to select the initial phase function of the initial parameter, which can directly affect the accuracy of the solution of the following minimization problem. The selection of the initial phase function can be based on the fault frequency band with high kurtosis found by the fast spectral kurtosis method, and then the center frequency of the fault band can be determined. A suitable initial phase function is selected for nonlinear matching pursuit of the newly formed signal to obtain the source signal component. In the field of actual mechanical fault diagnosis, the initial phase function can be expressed as $2\pi ft$, $f$ is the center frequency determined by the fast spectral kurtosis method. The specific flow chart is shown in Fig. 1.

4.1. Test equipment

In order to further verify the practicability of the proposed method, the vibration signals collected in the actual experiment is analyzed. Relevant experimental research is conducted in Guilin University of Electronic Science and Technology with the help of corresponding author Yanxue Wang in 2021.The MFS-MG mechanical fault comprehensive simulation test bench is used to conduct the test. The test system consists of a driving motor, a speed monitor, an acceleration sensor, a hand regulator, a flexible coupling, a load gearbox and an experimental bearing. The test bench is shown in Fig. 9. ER-12K type rolling bearings were selected for the test bearing, whose dimensions were shown in Table 1. Bearings with compound faults were installed near the driving motor and multi-channel vibration information was collected by piezoelectric acceleration sensors. The signal acquisition position was shown in Fig. 10. The faulty bearing is shown in Fig. 11.

Fig. 8Relationship between SNR and kurtosis

a) Simulation signal 1 and component 1 of the two methods

b) Simulation signal 2 and component 2 of the two methods

Fig. 9Bearing fault experiment platform

In the test process, the motor frequency $f_r=$29.87 Hz, sampling frequency $f_s=$25.6 kHz, analysis sampling number 25600, that is, analysis sampling time 1s. According to the basic dimension parameters of bearings, the frequency information of motor rotation and Eq. (16) and (17), the inner ring fault frequency is $f_i=$ 147.81 Hz, and the outer ring fault frequency is $f_o=$91.15 Hz:

where, $Z$ is the number of rolling bodies, and $d$ is the diameter of rolling bodies; $D$ is the pitch circle diameter and $\theta$ is the contact Angle.

Fig. 10Signal acquisition position of the experiment platform

Fig. 11Compound fault bearing

4.2. Bearing compound fault diagnosis

The bearing compound fault collected by the test equipment is a signal containing an inner ring fault and outer ring fault. And the time-domain waveform is shown in Fig. 12.

Fig. 12Time domain diagram of the signal of bearing composite fault

a) Time domain diagram of channel 1

b) Time domain diagram of channel 2

There is no obvious impact characteristic in the time domain vibration signal, and the fault type cannot be judged. Rapid spectral kurtosis analysis is performed on the bearing composite vibration signals collected in the experiment, as shown in Fig. 13.

Fig. 13Fast spectral kurtosis of experimental signal

Two fault bands with high energy are obtained by fast spectral kurtosis analysis and then filtered and analyzed by the STIMD method. The result is shown in Fig. 14.

Fig. 14Decomposition result waveforms of STIMD

a) Waveform diagram of outer ring fault component

b) Waveform of inner ring fault component

It is not difficult to find periodic impact signals in the time domain waveform in Fig. 14. The envelope spectrum analysis of the fault component obtained by the STIMD method is carried out, as shown in Fig. 15. $f_i$ is the inner ring fault frequency, and $f_o$ is the outer ring fault frequency.

Fig. 15Decomposing component envelope diagram of STIMD

a) Envelope diagram of outer ring fault component

b) Inner ring fault component envelope diagram

The existence of inner and outer ring faults can be determined from the envelope spectrum of components, and the effective separation and extraction of fault types can be realized. From Fig. 15(a), the outer ring failure frequency with a large value can be analyzed; from Fig. 15(b), the inner ring failure frequency and its double frequency can be analyzed. It is shown that the inner and outer ring faults can be decomposed and separated from the multi-channel bearing vibration signals by the STIMD method which selects the initial phase function by fast spectral kurtosis optimization, and the bearing compound fault diagnosis can be realized.

When the same test data is used to analyze the bearing vibration compound fault signals through the fast ICA method, the result is shown in Fig. 16.

The time domain signal component does not have an obvious impact component. The envelope analysis is carried out on the time domain component decomposed by the fast ICA method, and the envelope spectrum as shown in Fig. 17 is obtained. $f_i$ is the inner ring fault frequency.

Fig. 16Time domain diagram of fast ICA

a) Waveform diagram of component 1

b) Waveform diagram of component 2

Fig. 17Envelope spectrum of results of fast ICA

a) Envelope spectrum of component 1

b) Envelope spectrum of component 2

In the envelope spectrum of fast independent component analysis, only the inner circle fault can be detected, but the outer circle fault cannot be distinguished. Therefore, the composite bearing fault diagnosis method applied in the test cannot realize the separation and judgment of the fault frequency of the inner ring and the outer ring.

In conclusion, the spatiotemporal intrinsic mode decomposition method selects the key parameter of the initial phase function through fast spectral kurtosis optimization, which can effectively solve the problem that the fast independent component analysis method cannot effectively judge and separate the inner and outer ring faults in bearing compound fault diagnosis.