Time-varying fault feature extraction of rolling bearing via time–frequency sparsity

As a traditional method of processing signals, spectrum analysis can diagnose and locate machine faults with the help of the characteristic frequency [1, 2]. It has been widely used in fault identification of rolling bearings [3, 4]. However, traditional spectrum analysis only applies to the stable working condition with fixed speed. The rotation frequency and fault characteristic frequency (FCF) of each component will often change with rotation speed in engineering practice. Therefore, traditional spectrum analysis fails to solve this problem [5]. For accurate fault identification and early warning, it is necessary to eliminate the influence of speed and load change. Time–frequency analysis (TFA) has close ties with the tendency for the frequency characteristics of a signal to change with time. Thus, TFA can be used as a powerful tool to analyze time-varying fault characteristics under non-stationary working conditions [6, 7].

Presently, there are advantages and disadvantages to the commonly used TFA methods. Traditional TFA technology has various defects, such as the Heisenberg uncertainty principle, cross terms and mode aliasing [8–10]. These defects seriously interfere with representations of the signal features, leading to difficulty in determining the time-varying rule. Thetwo typical methods are short-time Fourier transform (STFT) and empirical mode decomposition (EMD) [11, 12]. In theory, STFT is based on the inner product and on the idea of fixing sliding windows. Therefore, the choice of window function and window width has an influence on the correctness of the analysis results. EMD belongs to the algorithm of signal adaptive decomposition. Nevertheless, EMD and its upgraded versions such as ensemble EMD still have some problems, such as a lack of theoretical basis and mode aliasing [13]. Subsequently, variational mode decomposition and empirical wavelet transform have been successively proposed to achieve multi-component signal decomposition [14, 15]. However, parameter selection is still inevitable. Daubechie et al proposed the synchrosqueezed wavelet transform (SWT) [16], which plays an important role in improving the energy concentration of the time–frequency (TF) plane both theoretically and practically. In this method, the complex spectrum of wavelet transform is squeezed and rearranged along the frequency axis to obtain a more concentrated TF representation. SWT has expanded into the field of multivariable signal processing [17]. However, this method is poor in noise immunity and requires an accurate estimation of the instantaneous frequency (IF). Recently, based on the theoretical framework of the synchrosqueezing transform (SST), several TFA methods, such as STFT-based SST (FSST), second-order STFT-based SST (FSST2), four-order STFT-based SST (FSST4) [18], synchroextracting transform (SET) [19], synchroextracting chirplet transform [20] and high-order SST [21], have been put forward to address these problems. These synchrosqueezing-based methods are all post-processing processors and do not change the qualitative behavior of the TF representation. Objectively, the actual analysis result of this kind of method still depends on standard STFT or wavelet transform (WT). Based on the above descriptions, optimal STFT plays a significance role in time-varying feature identification from the TF plane, especially under the circumstances of outliers and strong interference. To address this common problem, this paper introduces high-order STFT which is based on the high-order derivative of amplitude and phase of the signal [22]. In particular, considering the calculation cost and the feasibility of practical application, we derive a second-order short-time Fourier transform (SOSTFT) to obtain a sparse TF coefficient, which is based on second-order local estimation of the phase function.

In contrast with the commonly used post-processing method, this paper focuses on obtaining the optimal STFT coefficients. The commonly improved STFT methods include adaptive STFT [11, 23] and multi-resolution STFT [24]. Usually, these methods work well under certain assumptions. Nevertheless, significant reconstruction errors will appear if these methods are used to deal with signals with strong amplitude modulation and frequency modulation (AM–FM). In particular, the signals collected from industrial equipment often contain noise and are in non-stationary working conditions, which often leads to unsatisfactory results. Recently, Zhu has proposed two sound methods to reconstruct multi-component signals [25]. The methods are aimed at eliminating the influence of strong noise and outliers with a better loss function. However, it cannot precisely reconstruct strong modulation signals under the condition of loud noise.

Inspired by the idea that the TF transformation coefficient is the local maximum of the signal energy in the TF plane, we can innovatively transform the problem of time-varying feature extraction into a typical convex optimization problem. In other words, our goal is to get a sparse approximate representation in the TF plane. Convex optimization, a subfield of mathematical optimization, studies the problem of minimizing convex functions defined in convex sets. Due to the sound properties of convex optimization, we can obtain the approximate solution for many non-convex optimization problems in our daily life [26, 27]. This paper proposes a robust approach called second-order time–frequency sparse representation (SOTFSR), which defines SOSTFT as the matrix operator and the optimal TF coefficients as an objective function of the convex optimization problem, in contrast to extant TFA methods like SST. Based on the Huber function and minimax concave (MC) penalty, this paper introduces a multivariate generalized minimax concave (GMC) penalty to guarantee the strict convexity of the regularized linear least-squares cost function [28–30]. Therefore, the obtainment of STFT coefficients of interest through our proposed SOTFSR method is critical for subsequent feature extraction and mode reconstruction.

The rest of this paper is arranged as follows: in section 2, we will give a theoretical introduction to SOSTFT and SOTFSR. In section 3, numerical analyses of simulation signals with multi-components will be presented to verify the feasibility of the proposed method in dealing with strong frequency modulation signals. In addition, the results provided by STFT, SST and SET will be used for comparison. In section 4, we will analyze gravitational wave signals and vibration signals of mechanical equipment with inner faults under variable speed, which will further demonstrate the feasibility of the proposed method. Our conclusions will be given in section 5.

Zoom In Zoom Out Reset image size

2.1. Second-order short-time Fourier transform

The equation of a multi-component AM–FM signal is defined as follows:

$$\begin{equation}f\,(t) = \sum\limits_{k = 1}^K\,\,\, {{f_k}(t)} = \sum\limits_{k = 1}^K {{a_k}(t)} {e^{\,\,j{\phi _k}(t)}}\end{equation} \tag{ 1 }$$

where $K$ is a positive integer representing the number of AM–FM components, ${a_k}(t)$ and ${\phi _k}(t)$ are the instantaneous amplitude and instantaneous phase function of the $k{\text{th}}$ component, respectively, and ${\phi _k}^{\prime} (t)$ is determined as the IF. If we set $g(t)$ as a real window function in the Schwartz class, the STFT of signal $f\,(t)$ is expressed as:

$$\begin{equation}{S_{\text{f}}}(t,w) = \int_{ - \infty }^{ + \infty } {f\,(\mu )g(\mu - t){e^{ - jw(\mu - t)}}{\text{d}}} \mu. \end{equation} \tag{ 2 }$$

The first-order approximation STFT of signal $f\,(t)$ can be written as follows:

$$\begin{equation}{S_{\text{f}}}(t,w) \approx \sum\limits_{k = 1}^K {{a_k}(t)} {e^{j{\phi _k}(t)}}\hat g(w - {\phi _k}^{\prime} (t))\end{equation} \tag{ 3 }$$

where $\hat g(w)$ is the Fourier transform of window function $g(t)$. From equation (3), the original signal $f\,(t)$ can be easily obtained as:

$$\begin{equation}{f_k}(t) \approx \frac{{{S_{\text{f}}}(t,{\phi _k}^{\prime} (t))}}{{\hat g(0)}}.\end{equation} \tag{ 4 }$$

The modified STFT of f(t) in equation (2) can be rewritten as:

$$\begin{equation}{S_{\text{f}}}(t,w) = \int_{ - \infty }^{ + \infty } {f\,(t + \mu )g(\mu ){e^{ - jw\mu }}{\text{d}}} \mu. \end{equation} \tag{ 5 }$$

Using the Taylor series approximation, the complex-valued function $\ f\,(t + \mu )$ can be estimated as:

$$\begin{equation}f\,(t + \mu ) = f\,(t){e^{j\phi ^{\prime}(t)\mu }} \times \left[ {1 + \sum\limits_{n = 1}^\infty {\frac{1}{{n!}}{\rho _n}\left( t \right){\mu ^n}} } \right]\end{equation} \tag{ 6 }$$

where ${\rho _n}\left( t \right) = {\left. {\frac{{{d^n}}}{{d{\mu ^n}}}\left[ {\frac{{f\left( {t + \mu } \right){e^{ - j\phi ^{\prime}(t)\mu }}}}{{f\left( t \right)}}} \right]} \right|_{\mu = 0}}$. Substituting equation (6) into (5), we can obtain the high-order form of STFT. Details are as follows:

$$\begin{equation}{S_{\text{f}}}(t,w) = f\,(t) \times \left[ {\hat g(w - \phi ^{\prime}(t)) + \sum\limits_{n = 1}^\infty {\frac{{{j^n}}}{{n!}}{\rho _n}\left( t \right){{\hat g}^{\left( n \right)}}(w - \phi ^{\prime}(t))} } \right]\end{equation} \tag{ 7 }$$

where ${\overset{\frown}{g}}^{(n)}(w)$ denotes the nth-order derivative of $\stackrel{\frown}{g}(w)$ and $n = 1,2, \cdots ,N$. Equation (7) can obtain the Nth-order approximation of STFT and be expressed as follows:

$$\begin{equation}{S_{\text{f}}}(t,w) \approx {G_N}\left( {t,w} \right)F_N^T(t)\end{equation} \tag{ 8 }$$

where ${G_N}\left( {t,w} \right)$ is a feature analysis matrix and $F_N^T(t)$ is an original signal matrix.

In further detail, ${G_N}\left( {t,w} \right)$ and $F_N^T(t)$ are respectively calculated as:

$$\begin{align} {G_N}\left( {t,w} \right)\ =\ &[\hat g(w - \phi ^{\prime}(t)),\ j\hat g^{\prime}(w - \phi ^{\prime}(t)),\cdots,\nonumber\\ &\frac{{{j^N}}}{{N!}}{\hat g^{\left( N \right)}}(w - \phi ^{\prime}(t))] \end{align} \tag{ 9 }$$

$$\begin{equation}F_N^T(t) = [\,f\,(t),f\,(t){\rho _1}\left( t \right),\cdots,\,f\,(t){\rho _N}\left( t \right)].\end{equation} \tag{ 10 }$$

Then, the nth-order derivative of ${S_{\text{f}}}(t,w)$ can be estimated as:

$$\begin{equation}\frac{{{\partial ^n}}}{{\partial {w^n}}}{S_{\text{f}}}(t,w) \approx [{\hat g^{(n)}}(w - \phi ^{\prime}(t)),\cdots,\frac{{{j^N}}}{{N!}}{\hat g^{\left( {n + N} \right)}}(w - \phi ^{\prime}(t))]F_N^T(t).\end{equation} \tag{ 11 }$$

Combining (8) and (11), and setting $w = \phi ^{\prime}(t)$, $F_N^T(t)$ can be approximately calculated as:

$$\begin{equation}F_N^T(t) \approx {\left[ \begin{array}{ccc} \hat g(0)&\cdots&\frac{{{j^N}}}{{N!}}{{\hat g}^{\left( N \right)}}(0)\\[3pt] \hat g^{\prime}(0)&\cdots&\frac{{{j^N}}}{{N!}}{{\hat g}^{\left( {N{\text{ + 1}}} \right)}}(0)\\[3pt] \vdots&\vdots&\vdots\\ {{\hat g}^{(N)}}(0)&\cdots&\frac{{{j^N}}}{{N!}}{{\hat g}^{\left( {2N} \right)}}(0) \end{array} \right]^{{\text{ - }}1}}\left[ \begin{array}{c} \left.{S_{\text{f}}}(t,w) \right|_{w = \phi ^{\prime}(t)}\\[3pt] \left.\frac{\partial }{{\partial w}}{S_{\text{f}}}(t,w) \right|_{w = \phi ^{\prime}(t)}\\[3pt] \vdots\\ \left.{S_{\text{f}}}(t,w) \right|_{w = \phi ^{\prime}(t)} \end{array} \right].\end{equation} \tag{ 12 }$$

For simplicity, equation (12) can be expressed as a matrix product:

$$\begin{equation}F_N^T(t) \approx \bar G_N^{ - 1}{S_N}(t).\end{equation} \tag{ 13 }$$

From equation (13), the original signal can be reconstructed by the equation $f\,(t) \approx \left[ {\bar G_N^{ - 1}(1,1),\ldots,\bar G_N^{ - 1}(1,N + 1)} \right]{S_N}(t)$.

Concretely speaking, when $N = 2$, the second-order approximation of STFT for $f\,(t)$, i.e. SOSTFT, can be obtained from equation (8). It is derived to strike a balance between reducing computing costs and maintaining good TF analysis performance.

2.2. Second-order time–frequency sparse representation

Many signal-processing techniques can be regarded as a problem of the sparse approximation solution in the TF domain. Usually, the convex optimization scheme can be used to seek a sparse approximate solution to a system of linear equations. In this section, the SOTFSR method proposed herein introduces a novel non-convex penalty for sparse-regularized linear least squares, in order to strictly keep the minimized convexity of the least-squares cost function.

There is no doubt that standard STFT is an important theoretical basis for TFA. The STFT values of the analyzed signal are given at the sample points regarded as the approximate points of the time-varying fault features. Thus, the proposed SOTFSR can be promoted in the TF sparsity and transformed to solve the following cost function based on convex optimization:

$$\begin{equation}{c^ * } = \,\mathrm{arg}\, \mathop {\,\mathrm{min}\, }\limits_c \left\{ {G(c) = \frac{1}{2}\left\| {\,\,f - Ac} \right\|_2^2 + \lambda {\varphi _B}\left( c \right)} \right\}\end{equation} \tag{ 14 }$$

where $f \in {R^N}$ is an actual measured vibration signal and $c \in {R^{M \times N}}$ is the STFT coefficients. In terms of the sparse transformation, the matrix $A \in {R^{M \times N}}$ represents the inverse SOSTFT operator defined in equation (8) and $A{A^T} = I$. The non-convex penalty ${\varphi _B}$ is a sparsity-promoting penalty function, which is parameterized by a matrix B. The convexity of $G(c)$ largely depends on B. It also should be pointed out that the choice of B will depend on A. l₁ norm is a common penalty to introduce sparse representations, but it may produce local suboptimal solutions. This paper introduces a multivariate GMC penalty, which is based on the generalized Huber function [28–30].

Firstly, the Huber function $s(c)$ is defined as:

$$\begin{equation}s(c) = \left\{ \begin{gathered} \frac{1}{2}{c^2}\;\;\;\;\;\;\;\;\;\left| c \right| \leq 1 \hfill \\ \left| c \right| - \frac{1}{2}{} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} \left| c \right| \geq 1. \hfill \\ \end{gathered} \right.\end{equation} \tag{ 15 }$$

The MC penalty function is defined as:

$$\begin{equation}\phi (c) = \left\{ \begin{gathered} \left| c \right| - \frac{1}{2}{c^2}\;\;\;\;\;\;\;\;\;\left| c \right| \leq 1 \hfill \\ \;\;\;\;\frac{1}{2}{} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} \;\;\;\;\;\;\;\;{} \left| c \right| \geq 1. \hfill \\ \end{gathered} \right.\end{equation} \tag{ 16 }$$

Accordingly, the MC penalty can be expressed as:

$$\begin{equation}{} \phi (c){} = {} {} \left| c \right| - s(c). \end{equation} \tag{ 17 }$$

The above equation for the MC penalty will be used to generalize the GMC penalty. Then, we will define scaled versions of the Huber function and MC penalty, respectively. The scaled Huber function is defined as:

$$\begin{equation}{s_b}(c) =\left\{\begin{array}{ll} \frac{1}{2}{b^2}{c^2}&\left| c \right|\leq1/{b^2} \\[4pt] \left| c \right| - \frac{1}{{2{b^2}}}&\left| c \right|\geq1/{b^2}. \end{array}\right.\end{equation} \tag{ 18 }$$

Similar to equation (17), the scaled MC penalty function is defined as:

$$\begin{equation}{} {\phi _b}(c){} = {} {} \left| c \right| - {s_b}(c). \end{equation} \tag{ 19 }$$

Letting $B \in {R^{M \times N}}$, we define the generalized Huber function as:

$$\begin{equation}{S_B}(c) = \,\mathrm{inf}\, \left\{ {{{\left\| v \right\|}_1} + \frac{1}{2}\left\| {B(c - v)} \right\|_2^2} \right\}.\end{equation} \tag{ 20 }$$

Subsequently, we introduce a multivariate generalization of the MC penalty (17), which is generated by using the l₁ norm and the generalized Huber function defined in equation (20). The GMC penalty function is defined as:

$$\begin{equation}{\varphi _B}\left( c \right) = {\left\| c \right\|_1} - {S_B}(c)\end{equation} \tag{ 21 }$$

where ${S_B}(c)$ corresponds to the generalized Huber function. Matrix B is used to accurately describe the convexity of G(c), and its selection is related to matrix A, namely ${B^T}B \le 1/\lambda \left( {{A^T}A} \right)$.

Given the inverse SOSTFT coefficient matrix A, we can simply set:

$B = \sqrt {\gamma /\lambda } A$ $0 \leq \gamma \leq 1$(22) where the parameter $\gamma$ controls the convexity of the penalty ${\varphi _B}\left( c \right)$. When $\gamma = 0$, the penalty ${\varphi _B}\left( x \right)$ becomes the typical l₁ norm. According to our experiment, better results can be obtained if the value of $\gamma$ ranges from 0.5 to 0.8.

The minimum value of equation (14) can be solved by the forward/backward algorithm, which can be rewritten as a saddle point problem as shown below:

$$\begin{equation}({c^{{\text{opt}}}},{v^{{\text{opt}}}}) = \,\mathrm{arg}\, \;\mathop {\,\mathrm{min}\, }\limits_{c \in {R^N}} \;\mathop {\,\mathrm{max}\, }\limits_{v \in {R^N}} \;G(c,v)\end{equation} \tag{ 23 }$$

where $G(c,v) = \frac{1}{2}\left\| {y - Ac} \right\|_2^2 + \lambda {\left\| c \right\|_1} - \lambda {\left\| v \right\|_1} - \frac{\gamma }{2}\left\| {A(c - v)} \right\|_2^2$ is the saddle function and the matrix ${c^{{\text{opt}}}} \in {R^{M \times N}}$ is the optimal STFT coefficients by the proposed SOTFSR method. In addition, $v \in \left\{ {0,x - 1,x + 1} \right\}$ represents the auxiliary variable. When the local optimum value of equation (23) is obtained, the one-dimensional signal component of interest x can be obtained by $x = A{c^{{\text{opt}}}}$.

Based on this, assuming that the number of modes K is known, we can achieve accurate estimation of the fault characteristic component in the method described in [31]. This technique heavily relies on the minimization of the following energy function:

$$\begin{equation}\begin{aligned}{E_{\text{f}}}\left( {{\varphi _1}, \cdots ,{\varphi _K}} \right) = &\sum\limits_{k = 1}^K { - \int_R {{{\left| {T{F_{\text{f}}}({\varphi _k}(t),t)} \right|}^2}} } {\text{d}}t\\ & + \int_R {{\lambda _{\text{1}}}{{\varphi ^{\prime}}_k}{{(t)}^2}{\text{d}}t + {\lambda _{\text{2}}}{{\varphi {^{\prime\prime}}}_k}{{(t)}^2}{\text{d}}t} \end{aligned}\end{equation} \tag{ 24 }$$

where $\left( {{\varphi _1}, \cdots ,{\varphi _K}} \right)$ is the estimation of different modes with time-varying fault characteristics. $T{F_{\text{f}}}$ is one of the TF representations given by ${c^{{\text{opt}}}}$, which is provided in equation (23). ${\lambda _{\text{1}}}$ and ${\lambda _{\text{2}}}$ are regularization parameters which can maximize the balance between smoothness of ${\varphi _K}$ and energy. In general, the values of ${\lambda _{\text{1}}}$ and ${\lambda _{\text{2}}}$ have little effect on the final analysis results.

To verify the validity of the proposed method, numerical simulation analysis is firstly applied to multi-component signals. The ultimate goal is to accurately extract time-varying features and to reconstruct signals. The proposed method based on the convex optimization is employed to deal with noise signals. The analog simulation signal is expressed in the following form:

$$\begin{equation}{s_1}(t) = \,\mathrm{exp}\, ( - 0.005{t^2})\,\mathrm{sin}\, (2\pi (75t + \,\mathrm{cos}\, ({\text{2}}\pi t)))\end{equation} \tag{ 25 }$$

$$\begin{equation}{s_2}(t) = (1 - 0.01{t^2})\,\mathrm{cos}\, (2\pi (20t + 3{t^2}))\end{equation} \tag{ 26 }$$

where $f\,(t) = {s_1}(t) + {s_2}(t)$ is a typical multi-component AM–FM signal. The IF corresponding to ${s_1}(t)$ and ${s_2}(t)$ can be expressed as $I{{F}_1} = 75 + 2\pi \,\mathrm{sin}\, (2\pi t)$ and $I{{F}_2} = 20 + 6t$. To test noise immunity, Gaussian white noise is added to the signal $f\,(t)$ and the signal-to-noise ratio (SNR) is calculated as 3 dB.

The time- and frequency-domain waveforms of multi-component signals are respectively plotted in figure 2. Judging from this, traditional signal-processing methods, like fast Fourier transform (FFT), cannot effectively reflect the rule that frequency changes with time. Meanwhile, noise will interfere with identification of multi-component signals. The ideal time-varying feature relevant to multi-component signals is shown in figure 3. The red and blue TF curves respectively correspond to IF₁ and IF₂. It clearly shows the law of a mode with sinusoidal phase and a linear chirp. Subsequently, figure 4 shows TF analysis results obtained by four different methods, including STFT, SST and SET. Due to the influence of noise, the results provided by STFT are ambiguous in the TF plane in figure 4(a). SST and SET are typical TF rearrangement algorithms based on IF estimation. The calculation results are shown in figures 4(b) and (c). The results show that SST has low energy aggregation and SET shows a dramatic improvement in performance. However, irrelevant signal components still exist in figure 4(c). Compared with SST and SET, the proposed SOTFSR has clear advantages in sharpening TF representations in figure 4(d).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Renyi entropy is a means of quantitative evaluation of information uncertainty. The more random the signal is, the bigger the uncertainty and the corresponding entropy are, and vice versa. In the field of TF analysis, more concentrated TF energy means that the uncertainty is less and the Renyi entropy value is smaller. To quantitatively evaluate the effectiveness of TFA, we choose Renyi entropy as an assessment criterion [19]. The results obtained in different TF methods are listed in table 1. The proposed method clearly has a smaller Renyi entropy, which indicates that this method is better in TF energy aggregation.

To demonstrate the robustness to noise of the proposed method, we provide the SNR obtained by the three methods in different SNR simulation signals. Figure 5 represents the SNR changing curves obtained in these three methods. The vertical coordinate denotes the SNR after signal reconstruction. A larger SNR value in the vertical coordinate of figure 5 means a better performance in feature extraction. From the comparison shown in figure 5, we can conclude that the proposed method is superior in signal reconstruction.

Zoom In Zoom Out Reset image size

Gravitational waves are one of the most important predictions of general relativity [32]. After nearly a century of efforts, GW150914 was the first gravitational wave signalsuccessfully detected on September 14th 2015. The received signal proves that the prediction of general relativity is right. There are two black hole systems in the universe which may be combined into a larger one. We will use the method proposed in this paper to analyze gravitational wave signals. It is aimed at recognizing the law that gravitational wave intensity changes with time and frequency, which is observed by a Laser Interferometer Gravitational-Wave Observatory detector.

The observed initial Hanford signal has a length of 4096 samples and the corresponding time-domain waveform is shown in figure 6(a). Then, the TFA by means of STFT is plotted in figure 6(b). It is observed that the TF blur seriously affects the resolution of feature recognition. Then, gravitational wave signals are analyzed by an advanced method FSST2, and by the proposed method. The main reason for choosing FSST2 as the comparison method is that its computational complexity is acceptable. The analysis results are shown in figure 7. Clearly, the result provided by the proposed SOTFSR in figure 7(b) has fewer cross-interference phenomena and clearer TF representations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Then, the time-varying feature extraction is applied to figure 7 and the corresponding result is shown in figure 8. Theoretically, FSST4 can obtain more accurate estimations of the instantaneous frequencies of the modes making up the signal by using higher-order approximations for both the amplitude and phase. However, FSST4 will make calculation more complex, influencing its practical application in engineering. The TF representations generated by FSST4 are shown in figure 9. In the field of TF analysis, smaller Renyi entropy indicates the concentrated TF representation is better [19]. Through calculation, the Renyi entropy corresponding to figures 8 and 9 is 13.2889, 11.1412 and 12.0375, which demonstrates that the proposed method has satisfactory performance. We believe that the proposed method can detect the changing feature more effectively. The extracted feature is highlighted in green and manganese purple in the TF plane. Based on the time-varying feature extraction, the mode reconstruction is achieved. For the purpose of evaluation, we consider the difference between the reconstructed signal and the numerical relativity waveform obtained through independent calculation based on [32]. On the basis of the Baumgarte–Shapiro–Shibata–Nakamura formulation of Einstein's equations, Campanelli presents a new algorithm for obtaining the numerical relativity waveform [33]. Thus, numerical relativity waveforms can be considered as theoretical true values. The reconstruction error of the two TF methods is expressed in figure 10. It is certain that the proposed SOTFSR in this paper has a smaller reconstruction error and better performance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To verify the feasibility of the approach, this study conducts vibration signal analysis under time-varying rotational speed [34, 35]. This is aimed at achieving high-resolution representation of fault characteristic components in the TF plane under variable conditions. The data are collected by an ICP accelerometer placed on the faulty bearing of a SpectraQuest machinery fault simulator (MFS-PK5M) and the specific experimental equipment is shown in figure 11. A motor is used to drive the shaft and the rotational speed is adjusted by an AC drive. Two ER16 K ball bearings are located at the end of the shaft. It should be pointed out that the left bearing is in good condition while the right one requires replacement, i.e. is faulty. Table 2 clearly shows the structure parameters of the rolling bearings under study. The fault characteristic coefficient (FCC) corresponding to the outer race fault (FCC_o) and the inner race fault (FCC_I) can be calculated as follows [36].

$$\begin{equation}{\textrm{} FC{C_O}} = \frac{{{n_b}}}{2}\left( {1 - \frac{d}{D}\,\mathrm{cos}\, \varphi } \right)\end{equation} \tag{ 27 }$$

$$\begin{equation}{\textrm{} FC{C_I}} = \frac{{{n_b}}}{2}\left( {1{\text{ + }}\frac{d}{D}\,\mathrm{cos}\, \varphi } \right)\end{equation} \tag{ 28 }$$

Zoom In Zoom Out Reset image size

where n_b is the number of rolling elements, d is the diameter of the rolling element, D is the pitch diameter of the bearing, and $\varphi$ is the angle of the load from the radial plane. It is clear that FCC (FCC_o or FCC_I) is the ratio of the FCF to the shaft rotational frequency (SRF). Under time-varying speed conditions, the FCC can be used to identify the failure mode because it is determined by the bearing structural parameters and also immune to velocity variation. By identifying the instantaneous fault characteristic frequency (IFCF) and the instantaneous shaft rotational frequency (ISRF) from the extracted TF curves, the average curve-to-curve ratio can be determined by theoretical calculation. By comparing the actual calculated curve-to-curve ratio with the FCC determined by the bearing structural parameters, the bearing faults can be diagnosed under the condition of time-varying speed without using a tachometer to measure the rotational speed.

Moreover, an encoder, which is regarded as a tool for comparative analysis, is employed to measure the shaft rotational speed. The sampling frequency of the analyzed signals is set as 200 kHz at 10 s intervals. Down-sampling must be performed because the high sampling rate results in high data volume. In order to verify that this method is effective in many areas, we collect the vibration signals from a faulty bearing with an inner race defect and decrease the rotational speed from 25.1 to 13.1 Hz.

The down-sampled bearing vibration signal is shown in figure 12. From the calculation results of FFT operation plotted in figure 12(b), we cannot inspect the FCF component. Envelope spectrum analysis is a common method of vibration signal analysis of mechanical equipment and the corresponding results are shown in figure 13. Unfortunately, we still cannot effectively identify the fault characteristics, either from the time-domain waveform or spectrum analysis results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results provided by the envelope spectrum of the signal after 1 k high-pass filter are plotted in figure 14. The FCF, its harmonics, and the SRF cannot be identified in the spectrum analysis results. Generally, bearing faults can be diagnosed through identifying the FCF in the frequency spectrum of the vibration signals. However, under the time-varying rotational speed, the fault characteristics of bearings are often non-stationary. Therefore, traditional spectral analysis methods based on Fourier transform will be ineffective.

Zoom In Zoom Out Reset image size

Parametric TFA is an important technical means of realizing fault diagnosis of mechanical equipment under variable working conditions. To date, the traditional TF method STFT and its advanced version, second-order STFT-based SST (FSST2), still plays an important role in TF representations of multi-component signals in noisy environments. At the same time, STFT and FSST2 are chosen as analysis tools of comparison, so as to reduce calculation. The results generated by STFT and FSST2 are shown in figure 15.

Zoom In Zoom Out Reset image size

Due to the defects of the theory, a blur of TF occurs in the method of STFT, with a high level of noise in the TF plane displayed in figure 15(a). Therefore, it is necessary to improve the TF resolution of STFT. According to figure 15(b), in contrast to STFT, FSST2, which is based on the SST framework, is better at dealing with non-stationary vibration signals. However, the time-varying feature corresponding to the FCC is indistinct and the low-frequency speed signal is not recognized. In other words, the performance of FSST2 still cannot meet the analysis requirements in sharpening TF representation. Next, we adopt a new method, SOTFSR, which uses the second-order STFT framework and convex optimization scheme, with the results plotted in figure 16. Fortunately, FCCs and rotational speed signals are successfully extracted by the proposed method. Compared with figures 15 and 16, we can conclude that the proposed approach is better at fault feature extraction of complex multicomponent vibration signals.

Zoom In Zoom Out Reset image size

The number of feature curves that require extracting is set as 4. Subsequently, figure 17 shows the extraction of time-varying feature curves by the proposed SOTFSR. Through the extracted time-varying features, it is estimated that the bottom TF curve and the second-lowest TF curve respectively stand for ISRF and IFCF. By theoretical calculation, the average curve-to-curve ratio is determined as 5.47, which matches the FCC_I. The two highest curves in figure 17 are related to the FCFs twice and three times (2FCC_I and 3FCC_I). For clarity, figure 18(a) demonstrates the estimated ISRF and the identified IFCF. Moreover, the shaft encoder signal can provide the measured ISRF. Figure 18(b) shows a comparison between the estimated ISRF and the measured ISRF. The average error, |estimated ISRF–measured ISRF|/measured ISRF, is close to 2.41%. Thus, based on the proposed method of SOTFSR, the conclusion can be drawn that the identified failure mode belongs to the inner ring fault, which is consistent with the actual situation. Additionally, the estimated ISRF is so accurate that we can analyze orders without a tachometer in future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Time-varying fault feature extraction is an important technique of TF analysis and signal reconstruction. However, its performance relies on signal modulation and noise interference. To resolve these problems, we use GMC penalty functions with the help of a convex optimization to promote the TF sparsity of the coefficients in the SOSTFT domain. The main contents of this paper are summarized as follows: (1) we introduce a SOSTFT which is based on second-order local estimation of the phase function. The aim of SOSTFT is to obtain a high-resolution TF representation and set a TF transform operator in the following convex optimization scheme; (2) subsequently, the work of time-varying feature extraction can be innovatively transformed to a typical convex optimization problem. GMC penalty is introduced to maintain the convexity of the objective function so as to obtain the STFT coefficient of interest in the TF domain. On this basis, accurate time-varying feature identification can be accomplished by the proposed method of SOTFSR; (3) to verify the feasibility of this approach, we carry out the numerical simulation of multi-component signals from a faulty bearing under variable speed. It is worth mentioning that the extracted feature corresponding to ISRF can be used for order-tracking analysis in the process of dealing with practical engineering signals. Consequently, the result proves that the proposed SOTFSR has clear superiority in dealing with vibration signals of mechanical equipment.

This work was supported by the National Natural Science Foundation of China under Grant No. 51805382.

The authors declare that there are no conflict of interests regarding the publication of this article.