Fault feature extraction method for rolling bearing based on MVMD and complex Fourier transform

2.1. Multivariate modulated oscillations

To express conveniently a multi-oscillation signal, the amplitude modulation and frequency modulation signal is defined in the form of vector ${\left\{u_i\right(t\left)\right\}}_{\mathrm{i}=1}^C$ in the real number field:

where $a_i\left(t\right)$ and ${\varnothing }_i\left(t\right)$ are the amplitude function and the phase function, respectively, and $C$ is the number of channels. Hilbert operator is applied to Eq. (1), and the corresponding analytic vector of $u\left(t\right)$ is obtained. The analytic vector $u_{+}\left(t\right)$ is as follows [30]:

The amplitude function $a_i\left(t\right)$ and the phase function ${\varnothing }_i\left(t\right)$ can be obtained by the following modulo and phase operations:

The original real-value signal vector $u\left(t\right)$ can be obtained from the real part of $u_{+}\left(t\right):$

For vibration signals of rolling bearing defects collected by 2D or 3D sensors, so one or more components with the same frequency may be present. Therefore, a simple multivariate analytic signal model for $u_{+}\left(t\right)$ that assumes a single common component among all data channels is proposed:

2.2. Multivariate VMD

MVMD is an extension of the VMD method. For more details regarding the VMD, see [11]. As a generalization of the VMD algorithm in multidimensional space, the main goal of MVMD is to extract the predefined $K$ multivariate modulation oscillation signals $u_k\left(t\right)$ containing C channel data from the input data $u\left(t\right)$. MVMD is used to decompose a multivariate signal $u\left(t\right)$ into the sum of multiple modulated oscillation signals $u_k\left(t\right)$:

where $u_k\left(t\right)=[u_{k1}\left(t\right),u_{k2}\left(t\right),\dots ,u_{kC}\left(t\right)]$.$\mathrm{}K$ is the number of multivariable modulated oscillation signals.

The purpose of MVMD is to extract a set of multivariable modulation oscillations ${\left\{u_k\left(t\right)\right\}}_{k=1}^K$ from the input data. Additionally, the sum of bandwidths of the extracted sets is the smallest, and the original signals can be reconstructed accurately. To this end, Eq. (2) is used to represent an analytical vector, and the bandwidth of $u_k\left(t\right)$ can be estimated by the L2 norm of the gradient function of $u_{+}^k\left(t\right)$. Therefore, the multivariable function that needs to be optimized in MVMD is expressed by the following formula [30]:

The symbol ${\partial }_t$ represents the partial derivative operation on time.

Another characteristic of Eq. (8) is that a single frequency component ${\omega }_k$ is used for harmonic mixing of the entire vector $u_{+}^k\left(t\right)$. Therefore, the design provides the possibility to find the multivariate oscillation with a single common frequency component ${\omega }_k$ in all channels in $u_k\left(t\right)$, and this causes a multivariate modulation oscillation model defined in Eq. (6). The bandwidth of the modulation multivariable oscillation is estimated by shifting the single-side spectrum of all channels of $u_{+}^k\left(t\right)$ to ${\omega }_k$ and taking the Frobenius norm of the matrix. It is convenient to express $f$ as follows:

where $u_{+}^{k,c}\left(t\right)$represents an analytical modulation signal corresponding to channel number $C$ and mode number $k$. In contrast to the vector signal in Eq. (8), $u_{+}^{k,c}\left(t\right)$ is a single component complex value signal. The constrained optimization problem of MVMD is as follows:

where $\sum _k{u}_{k,c}\left(t\right)=x_c\left(t\right),c=\mathrm{1,2},\cdots ,C$, ${\{u}_{k,c}\left(t\right)\}$ and $\left\{{\omega }_k\right\}$represent the set of all $k$ modes and their central frequencies, respectively.

The corresponding augmented Lagrangian function becomes:

where $\alpha$ is the penalty factor.

In the MVMD, the alternating direction multiplier (ADMM) is used to solve the above unconstrained optimization problem. The ADMM method transforms the complex optimization problem in Eq. (11) into several simpler sub-optimization problems through Eq. (12-14):

where $\tau$ is the time step.

Eq. (12) can be equivalent to the following formula:

One of these sub-problems is depicted in the updated diagram of the original VMD. In the MVMD, the relationship is updated using the following pattern [12] in the Fourier domain:

For the renewal of the central frequency of Eq. (13), considering that the last two terms of the Lagrangian function Eq. (11) do not depend on ${\omega }_k$, the related problems are simplified as follows:

By using Plancherel’s theorem of the inner product of the function in the time domain and frequency domain, the optimization process can be carried out more easily in the frequency domain. In the Fourier domain, the equivalence problem of Eq. (17) becomes:

By setting the first derivative of the quadratic function to zero and then performing some simple algebraic operations, the above sum of the quadratic function is minimized, and the following relations are obtained:

It can be seen from Eq. (19) that when updating ${\omega }_k$ of each mode during MVMD, the power spectrum contribution of all channels is considered.

2.3. Complex Fourier transform

Fourier transform is performed on the basis of a complex sequence signal $z\left(i\right)=x\left(i\right)+\mathrm{j}y\left(i\right)$ ($i$ is an integer and $i\in \left[1,N\right]$, $N$ is the sequence length,$\mathrm{}{\mathrm{j}}^2=-1$) to obtain $Z\left(I\right)$:

where FFT stands for fast Fourier transform.

Using the Fourier transform in the finite-length sequence of real and even functions on the $N$/2 symmetry and the Fourier transform of the pure imaginary odd-function finite-length sequence on the $N$/2 anti-symmetric and the linear properties of the Fourier transform, it is possible to get:

where $X\left(I\right)$ and $Y\left(I\right)$ are the Fourier transforms of the cosine signals $x\left(t\right)$ and $y\left(t\right)$ respectively, $I=$1, 2,…, $N$/2.

Because the amplitude of the Fourier transform of the complex sequence $z\left(i\right)$ is related to $\left|z\left(i\right)\right|$, if it is represented by a trigonometric function, it can be found that $\left|z\left(i\right)\right|$ is related to the initial phases of $x$ and $y$.

In order to study the relationship between CFT and FVS when the initial phase changes, let:

where $f=$ 50 Hz, $\phi =k\pi /8$, ($k=$ –8, –7,…, 0,…, 7, 8).

The sampling frequency $f_s$ is equal to 800 Hz. The characteristic frequency amplitude obtained by analyzing $z\left(t\right)$ with the phase difference of its imaginary and real parts based on the CFT and FVS methods, respectively, is shown in Fig. 1. In Fig. 1, when $\phi \in (-\pi ,0)$, the amplitude of the characteristic frequency based on CFT is larger. When $\phi \in (0,\pi )$, the amplitude of the conjugate characteristic frequency ($f_s-f$) based on CFT is larger. Therefore, the characteristic frequency of the complex signal can be selected according to the phase difference to characterize the fault.

Fig. 1Characteristic frequency amplitude obtained by analyzing zt with phase difference of its imaginary and real parts based on CFT and FVS methods

2.4. Proposed feature extraction method for determining rolling bearing defects

1) First, a 3D accelerometer mounted at rolling bearing housing is used to collect vibration signals from two channels. The collected horizontal signal $x\left(t\right)$ is used as the signal $u_1\left(t\right)$ of the first channel, and the collected vertical signal $y\left(t\right)$ is used as the signal $u_2\left(t\right)$ of the second channel. Then, a binary modulated vibration signal $u\left(t\right)$ is formed, where $u\left(t\right)=\left[\genfrac{}{}{0pt}{}{u_1\left(t\right)}{u_2\left(t\right)}\right]$.

2) MVMD is applied to decompose $u\left(t\right)$ into the sum of binary oscillation signals ${\{u}_k\left(t\right)\}$, $k=1,\mathrm{}2,\cdots K$. where $u_k\left(t\right)=[u_{k1}\left(t\right),u_{k2}\left(t\right)]$.

3) Hilbert operations are respectively performed on $u_{k1}$ and $u_{k2}$, and the envelope function $a_{k1}$ of $u_{k1}$ and the envelope function $a_{k2}$ of $u_{k2}$ are obtained respectively.

4) Let $a_k=a_{k1}+\mathrm{j}{a}_{k2}$, ${\mathrm{j}}^2=-1$. Then, Eq. (21) is applied to $a_k$ to obtain a complex Fourier spectrum, and to select the fault feature frequency according to the phase difference between the imaginary part and the real part of $a_k$.

The key point of the whole procedure is $K$ selection. Generally, the vibration signals caused by rolling bearing defects have four components: fundamental frequencies, defect frequencies, natural vibration frequencies and ultrasonic frequency. Therefore, when the MVMD method is used to analyze vibration signals caused by rolling bearing defects, $K$ is set to 4.

3.1. Experimental Setup

In this paper, the test bench shown in [31] is used. As described in [31], the experiment rig consists of a rotor-bearings system driven by a 1 HP triple-phase asynchronous motor, as shown in Fig. 2(a). Multiple rotating speeds, which fluctuation range is almost 60 rpm, can be monitored by a motor speed controller. The vibration data were obtained from 3D-accelerometer mounted on a rotor-bearing system. Different kinds of bearing faults are simulated by replacing bearing 1 with a rolling element defective bearing or compound defective bearing having roller defect and outer race defect. The defect on the outer race and ball of bearings is shown in Fig. 2(c) and (d). All channels of vibration data from the 3D-accelerometer were acquired simultaneously using a LMS SCADAS Mobile data acquisition system in Fig. 2(b) at a sampling rate of 25.6 KHz. The bearing model is MB ER-16K, which has 9 rollers, 38.5 mm pitch diameter, 7.9 mm rolling body diameter, and 0 deg. contact angle. Ball pass frequency outer (BPFO) is 3.572. Ball Spin Frequency (BSF) is 2.322. The motor speed is 2700 rpm, and the rotation frequency $f_n=$45 Hz. So, the ball fault frequency $f_b$ is 104 Hz, and the outer race fault frequency $f_o$ is 161 Hz. Ball defect of 1.2 mm wide and 0.5 mm deep was made by wire cutting method. Outer race defect of 1.2 mm wide and 1.75 mm deep was made by wire cutting method.

Fig. 2View of experimental setup of rotor-bearing system: 1 – variable speed controller; 2 – motor; 3 – coupling; 4 – 3D-accelerometer; 5 – bearing 1; 6 – rotor disk; 7 – shaft; 8 – healthy bearing; 9 – ball defect; 10 – outer race defect; 11 – LMS SCADAS mobile

a) Machine fault simulator (MFS)

b) Data acquisition system

c) Oter race defect

d) Ball defect

3.2. Rolling element defect

In Fig. 2, rolling bearings are mainly under the action of the radial force excited by the movement of rotor disk and the motor, so this paper only analyzes the radial vibration signals. If the collected vibration signal in the horizontal direction is x and the vibration signal in the vertical direction is $y$, then a binary modulation signal $u$ is equal to $u=\left[\genfrac{}{}{0pt}{}{x}{y}\right]$. Let $K=$4, $\alpha =$2000 and $\tau =$ 0, the MVMD method is used to decompose $u$ into the sum of a series of IMFs $u_k$ ($k=$ 1, 2, 3, 4), where $u_k=[u_{k1},u_{k2}]$. $u_{k1}$ represents the horizontal component of the vibration signal $x$, and $u_{k2}$ represents the vertical component of the vibration signal $y$. The time-domain waveforms of $x$, $y$ and $u_k$ are shown in Fig. 3. The evolution of IMFs center-frequencies is shown in Fig. 4.

Fig. 3Time-domain waveforms of rolling element defect signals x, y and uk obtained by MVMD

In Fig. 3(a) and (b), the original vibration signals $x$ and $y$ have no obvious amplitude modulation characteristics, and the fault feature information is submerged by noise. In Fig. 3(g) and Fig. 3(h), the $u_{31}$ and $u_{32}$ contain very obvious modulation characteristics, but the amplitude of $u_{32}$ is larger than that of $u_{31}$. In Fig. 4, the central frequency of IMFs fluctuated greatly at the early stage of evolution, but it entered a stable state after 24 evolutions.

Fig. 4Evolution of IMFk center-frequencies ωk

To illustrate the frequency band structure of a vibration signal excited due to the bearing defect, the reason for taking $K=$4, and the advantage of using the MVMD to analyze multi-channel signals, the Fourier spectrums of the original signals and the Fourier spectrums of their decomposition results based on MVMD are shown in Fig. 5. As generally accepted, the Fourier spectrums of the rolling bearing fault vibration signals contain four subbands. Therefore, it is appropriate to set $K$ to 4 when using MVMD to decompose the rolling bearing fault vibration signals. In addition, $u_{k1}$ and $u_{k2}$ have similar spectral structures and are arranged in order of frequency from low to high, and there are no modal aliasing and modal splitting between the subbands of $u_{k1}$ and the subbands of $u_{k2}$. The fault features obtained by using the MVMD method are decomposed into the same layer, which is beneficial to the subsequent fault feature fusion. Hilbert operations are respectively performed on $u_{k1}$ and $u_{k2}$, and the envelope functions $a_{k1}\mathrm{}$and $a_{k2}\mathrm{}$are obtained. The Fourier spectrums of $a_{k1}\mathrm{}$and $a_{k2}$ are shown in Fig. 6.

Fig. 5Fourier spectrums of rolling element defect signals x, y and uk obtained by MVMD

In Fig. 6, the rolling element fault feature frequency $f_b=$101 Hz is in the mid-frequency subband 3. In addition, the amplitude of $a_{31}$ is much larger than that of $a_{32}$.

Let $a_3=a_{31}+\mathrm{j}{a}_{32}$, then the phase difference spectrum will be between $a_{32}$ and $a_{31}$. So it is possible to depict the complex Fourier spectrum of $a_3$ and full vector spectrum of $a_3$ as shown in Fig. 7. Because the phase difference between $a_{32}$ and $a_{31}$ is between –$\pi$ and 0, the amplitude of the feature frequency $f_b$, which is 3.563, in the complex Fourier spectrum of $a_3$ is larger than that of the conjugate frequency ($f_s$-$f_b$), which is 3.456. According to Fig. 6 and Fig. 7, it can be seen that the amplitude of $f_b$ is the largest in the complex Fourier spectrums, i.e. 0.3563, followed by 0.3506 in the full vector spectrum, and it is the smallest in the Fourier spectrum of $a_{32}$, i.e. 0.3223.

Fig. 6Fourier spectrums of ak1 and ak2

Fig. 7Phase difference spectrum, complex Fourier spectrum and full vector spectrum of a3

The Fourier spectrums of the first three-order IMFs $c_{k1}$ and $c_{k2}$ obtained by the BEMD method are shown in Fig. 8. The Fourier spectrums of first three-order product functions (PFs) $p{f}_{k1}$ and $p{f}_{k2}$ obtained by the CLMD method are shown in Fig. 9. The envelope spectrums of the IMFs and PFs are shown in Fig. 10 and Fig. 11, respectively.

Obviously, the IMFs obtained by the BEMD method or the PFs obtained by the CLMD method induce modal aliasing and mode splitting, which will cause the extracted fault features to be indistinct. For example, Fig. 8(a) shows that the mid-frequency subband is almost submerged by the dominant high-frequency subband, and mode splitting occurs in the mid-frequency subband. As a result, the amplitude of the rolling element fault feature frequency $2f_b$ in Fig. 10(a) is very small and almost submerged. Fig. 8(b) demonstrates mode aliasing in $c_{12}$, and the mid-frequency subbands containing fault features are decomposed into $c_{12}$ and $c_{22}$, obviously the mode splitting occurs. As a result, the amplitude of the fault feature frequency $f_b$ in Fig. 10(b) and Fig. 10(d) is reduced. when applying the CLMD method to analyze rolling element defects, the similar situation as described above occurs, as shown in Fig. 11.

Fig. 8Fourier spectrum of first three-order IMFs ck1 and ck2 obtained by BEMD method

Fig. 9Fourier spectrum of first three-order PFs pfk1 and pfk2 obtained by CLMD method

In Fig. 10 and Fig. 11, the fault features on the same layer are inconsistent. This means the information fusion becomes increasingly challenging. Comparing the envelope spectrums of IMFs and PFs, the amplitude of the fault characteristic frequency of the vibration signal in the vertical direction is obviously larger than that in the horizontal direction. It is very likely to make a wrong diagnosis based on only the signal feature in the horizontal direction.

Fig. 10Envelope spectrums of ck1 and ck2 based on BEMD

The full vector envelope spectrums of $c_1$ and $c_2$ are obtained by applying the FVS to fuse the envelope signals of $c_{k1}$ and $c_{k2}$ ($k=$1, 2) as shown in Fig. 12. The complex envelope spectrums of $p{f}_1$ are obtained by applying the CFT to fuse the envelope signals of $p{f}_{11}$ and $p{f}_{12}$ as shown in Fig. 13. The maximum amplitudes of the rolling element fault feature frequency $f_b$ obtained by the above methods are shown in Table 1. The amplitude of $f_b$ obtained by the MVMD method is much larger than that obtained by the BEMD or CLMD methods. The amplitude obtained by the proposed MVMD+CFT method is larger than that obtained by the MVMD+FVS method.

Fig. 11Envelope spectrums of pfk1 and pfk2 based on CLMD

Fig. 12Full vector envelope spectrums of c1 and c2 based on BEMD

Fig. 13Complex Fourier spectrum of envelope signal of pf1 based on CLMD

3.3. Composite defects

The time-domain waveform and Fourier spectrum of the composite fault vibration signals with rolling element defects and outer race defects are shown in Fig. 14. Similar to rolling element defects, the spectrums of composite fault vibration signals mainly contain 4 subbands.

Let $K=$4, $\alpha =$2000 and $\tau =$0. Then the MVMD is applied to decompose the composite fault signals into a series of $u_{k1}$ and $u_{k2}$ ($k=$1, 2, 3, 4). Their Fourier spectrums are shown in Fig. 15.

It is clearly seen from Fig. 14 and Fig. 15 that $u_{11}$-$u_{41}$ and $u_{12}$-$u_{42}$ come from subband 1, subband 2, subband 3 and subband 4 in the spectrum of signals $x$ and $y$, respectively, bearing no modal aliasing or modal splitting in the spectrums of $u_{k1}$ and $u_{k2}$. This property of MVMD provides the possibility to determine where the fault signature is located. In Fig. 15, the frequency centers of the subbands in the spectrums of $u_{k1}$ and $u_{k2}$ are approximately the same that is very beneficial to the subsequent fault feature information fusion.

Fig. 14Time-domain waveform and Fourier spectrum of composite fault vibration signal

Fig. 15Fourier spectrums of uk1 and uk2 obtained by MVMD

The HT is applied to $u_{k1}$ and $u_{k2}$ respectively, and the envelope spectrums of $u_{k1}$ and $u_{k2}$ are shown in Fig. 16. The rolling element defect feature frequency $f_b=$101 Hz is still in the mid-frequency subband 3. The outer race defect feature frequency $f_o=$160 Hz is in low-frequency subband 1, subband 2 and in high-frequency subband 4. Compared with Fig. 6, it can be seen that when only the rolling element defect occurs, the fault feature amplitude in the vertical direction is larger than that in the horizontal direction. But when a composite defect occurs, the situation changes cardinally to the reverse.

Let $a_k=a_{k1}+j{a}_{k2}$, where $a_{k1}$ and $a_{k2}$ are the envelope signals of $u_{k1}$ and $u_{k2}$ respectively. Then the phase difference spectrums of $a_{k1}$ and $a_{k2}$ shown in Fig. 17, and the complex Fourier spectrums and the full vector spectrums of $a_k$ shown in Fig. 18 and in Fig. 19 are actual.

In Fig. 18, the larger amplitude is selected from the characteristic frequency and conjugate characteristic frequency to reveal the characteristics of composite defects according to the phase difference between $a_{k1}$ and $a_{k2}$. Taking $a_3$ and $a_4$ as two examples, In Fig. 17(c) and (d), the phase difference ${\phi }_3$ between $a_{32}$ and $a_{31}$ at the ball defect characteristic frequency $f_b$ is 6.481° and the phase difference ${\phi }_4$ between $a_{42}$ and $a_{41}$ at the outer race defect characteristic frequency $f_o$ is 82.31°, which belongs to (0, $\pi$). According to 2.3, when the phase difference is between (0, $\pi$), the magnitude of the conjugate characteristic frequency is larger than that of the characteristic frequency. In Fig. 18, the magnitude of the ball defect characteristic frequency $f_b$ is 0.2053, while the magnitude of the conjugate characteristic frequency $f_b$ is 0.2228. The magnitude of the outer race fault characteristic frequency $f_o$ is 0.238, while the magnitude of the conjugate characteristic frequency ($f_s$-$f_o$) is 0.5825. In Fig. 19, the amplitude of $f_o$ in the full vector spectrum of $a_4$ is only 0.4102. In Fig. 16, the amplitude of $f_o$ in the envelope spectrum of $u_{42}$ is larger than that of $u_{41}$, i.e. 0.4094, but much smaller than that in the full vector spectrums or in the complex Fourier spectrums.

Fig. 16Envelope spectrums of uk1 and uk2

Fig. 17Phase difference spectrums of ak1 and ak2

The BEMD is used to decompose the composite defect signals into a series of IMFs $c_{k1}$ and $c_{k2}$, where $k=$1, 2, 3, 4. The Fourier specturms of the first four $c_{k1}$ and $c_{k2}$ are shown in Fig. 20. Their envelope specturms are shown in Fig. 21. Obviously, there are modal aliasing and modal splitting in the spectrum of IMFs shown in Fig. 20. The envelope spectrums of $c_{k1}$ and $c_{k2}$ have significant differences as shown in Fig. 21. The outer race fault feature frequency $f_o$ is in the envelope spectrums of $c_{11}$ and $c_{41}$, and the rolling element fault feature frequency $f_b$ is in $c_{21}$. In the envelope spectrum of $c_{k2}$, each of them has an outer race fault feature frequency $f_o$, but without rolling element fault feature frequency $f_b$. Obviously, the diagnosis based only on the envelope spectrum of $c_{k2}$ is easy to induce errors and inaccuracies.

Fig. 18Complex Fourier spectrums of ak based on MVMD

Fig. 19Full vector spectrums of ak based on MVMD

Fig. 20Fourier spectrums of first four ck1 and ck2 obtained by BEMD

Fig. 21Envelope spectrums of first four ck1 and ck2 obtained by BEMD

The full vector envelope spectrums of the IMFs $c_k$ obtained by BEMD are shown in Fig. 22. The complex envelope spectrums of $p{f}_k$ obtained by CLMD are shown in Fig. 23. In Fig. 22 and Fig. 23, the rolling element fault features are small and difficult to be identified.

The maximum amplitudes of the compound defect feature frequencies $f_b$ and $f_o$ obtained by the above methods are shown in Table 2. It can be seen from Table 2 that even if compound defects are analyzed, the amplitude of defect features extracted by the proposed method is clearer.

Fig. 22Full vector envelope spectrums of IMFs ck obtained by BEMD

Using fast spectral kurtosis based on 4th order statistics, it is possible demodulate the rolling element defect feature frequency $f_b$ and its frequency multiplication from the signal $x$, and the outer race defect feature frequency $f_o$ and its frequency multiplication from the signal $y$. The demodulation result is shown in Fig. 24. Obviously, only the rolling element defect is diagnosed according to the envelope spectrum of $x$, and the only outer race defect is diagnosed according to the envelope spectrum of $y$.

Fig. 23Complex envelope spectrums of first two pf obtained by CLMD