TSCK guided parameter convex optimization tunable Q -factor wavelet transform and its application in wheelset bearing fault diagnosis

Abstract

Wheelset bearing is a typical vulnerable structural component in high-speed trains and heavy haul vehicles. In addition to the typical nonlinear and nonstationary characteristics, the vibration signal of wheelset bearing also contains track subgrade vibration and transmission path coupling interference components. To solve this problem, this paper proposes a new feature extraction method for wheelset bearing faults. This method constructs the Teager energy spectrum correlation kurtosis, which is purposely sensitive to periodic fault impulse components, as the objective function. The Q-factor and redundancy of tunable Q-factor wavelet transform are selected by using the parameter convex optimization method, which makes the signal decomposition have better sparsity, so as to extract fault information accurately. Simulated analysis, experimental signal analysis of QPZZ-II test-bed, and experimental signal analysis of wheelset bearing test-bed show that the proposed method can suppress the influence of nonperiodic transient impulse components, harmonic components, and noise components in the signal and accurately extract the periodic impact characteristics of bearings.

Introduction

Wheelset bearings are widely used in railway trains, subways, and heavy loaded vehicles with large equipment, and their operation status is related to the safety of vehicle operation. Under heavy load operation, especially under uneven heavy load conditions, it is very easy to cause wheelset bearing damage. When the vehicle runs under wheelset bearing damage, it is easy to introduce subgrade vibration, which will further damage the mechanical structure of the vehicle, and in serious cases, it will lead to seizure or even overturning. Therefore, it is of great significance to adopt necessary analysis means to enable the wheelset bearing to realize monitoring and diagnosis during vehicle operation.^1–5

Vibration signal is an effective medium to analyze the running state of rotating machinery equipment, especially bearings. Many scholars have carried out a lot of research work on the nonstationary and nonlinear characteristics of bearing vibration signals and fault feature extraction methods.^6–10 In order to effectively suppress noise and signal modulation and accurately identify wheelset bearing fault impact, Li et al. ¹¹ proposed a method using morphological signal and image-processing technology. This method mainly includes two aspects: a new double cross-correlation algorithm for noise reduction and an improved image-processing algorithm for highlighting fault features. Traditional boundary determination methods cannot evaluate the capacity of fault information contained in the decomposition frequency band. To solve this problem, Zhang et al. ¹² proposed an adaptive boundary determination method based on time delayed kurtosis (TDK), which uses the boundary determined by TDK to construct Meyer wavelet to extract fault-related resonance bands, and the proposed method was applied to wheelset bearing fault diagnosis. The parameter dictionary obtained by the classical correlation-filtering algorithm is difficult to match the potential fault influence characteristics of the wheelset bearing vibration signal. To solve this problem, Deng et al. ¹³ proposed a fault impact-matching algorithm (FIMA) combined with orthogonal match pursuit. The core of FIMA mainly includes two stages: partial segmentation and global analysis. The correlation function and kurtosis are used to comprehensively evaluate the local and global matching degree of Laplace wavelet and signal. In view of the problem that traditional sparse representation is difficult to ideally extract fault features when the signal contains strong interference components, Xing et al. ¹⁴ proposed a new method of fault feature extraction for wheel bearing based on wavelet sparse representation and adaptive local iterative filtering. The adaptive local iterative filtering of this method effectively reduces the influence of interference components and is helpful to extract sparse pulses. Wavelet sparse representation uses the global optimal solution of L1 regular optimization in sparse coding to extract the inherent features of faults in the wavelet domain. Based on the developed autocorrelation kurtosis and the proposed spectral segmentation method, Zhang et al. ¹⁵ proposed a new wheel bearing fault detection method called adaptive autocorrelation kurtosis. The autocorrelation kurtosis of this method aims to reduce the interference of irrelevant components and improve the signal-to-noise ratio. In view of the poor adaptability of time-frequency analysis method under multiple working conditions, Wei et al. ¹⁶ proposed a bearing fault detection method combining empirical mode decomposition and adaptive time-varying parameter short-time Fourier synchronous compression transform, which solved the problem of signal adaptation under multiple working conditions. Aiming at the influence of complex transmission path and strong background noise on the vibration signal of early fault of rolling bearing, Xiong et al. ¹⁷ proposed a diagnosis method of parameter adaptive multipoint optimal minimum entropy deconvolution adjustment combined with dynamic mode decomposition. The analysis of analog signals, public database signals, and fan bearing signals shows that the proposed method can successfully extract the early fault features of rolling bearings. To solve the problem of regularization parameter selection of sparse denoising algorithm, Sun et al. ¹⁸ proposed an adaptive regularization parameter selection strategy, which directly measures the importance of fault features in sparse denoising signals. The minimum angle regression algorithm with high computational efficiency is used to obtain the entire solution path, and the feature probability index is designed to measure the importance of fault features in the reconstructed signal envelope demodulation spectrum. This method is not limited by the noise probability distribution. In order to solve the problem that the system parameters need to be matched in real time when stochastic resonance (SR) processes nonstationary signals, Wang et al. ¹⁹ studied the adaptive frequency shift SR of tristable nonlinear systems. The variation of potential well depth, bifurcation phenomenon, and steady state switching process caused by system parameters are studied in order to control system parameters and determine steady state conditions. Moreover, SR is applied to fault diagnosis of rolling bearings under variable speed conditions and strong noise background, and weak fault features are effectively extracted. In order to effectively analyze the fault signals of nonlinear and nonstationary rolling bearings, Li et al. ²⁰ studied derivative and enhanced discrete analysis wavelet algorithms for rolling bearing fault detection. The morphological spectrum of fault signal is obtained through multi-scale morphological open operation. Yazdanianasr et al. ²¹ constructed an improved envelope spectrum via feature optimization-gram (IESFOgram) method, which can accurately diagnose the wear bearing running at low speed. Amin et al. ²² present a machine-learning framework for early damage detection in gearboxes based on the cyclostationary and Kurtogram analysis of sensor data. Li et al. ²³ proposed a fault diagnosis method based on weighted variational mode decomposition (VMD) and cyclic spectrum slice energy (CSSE). The signal is decomposed into intrinsic mode functions (IMFs) through VMD, and the amount of information contained in each IMF is measured using sparsity. Then, the reconstructed signal is analyzed using the advantage that CSSE can accurately mediate fault information. Feng et al. ²⁴ provides a vibration-based wear mechanism identification procedure, and then the wear evolution is tracked using an indicator of vibration cyclostationarity. Ni et al. ²⁵ proposed a fault information guided VMD (FIVMD) method to extract the weak bearing repetitive transient for the problem of parameter selection of VMD. Two nested statistical models based on fault cyclic information are used to approximate the number of modes in combination with the statistical threshold of a specific significance level. Then ratio of fault characteristic amplitude is defined and used to identify the optimal bandwidth control parameters, so as to extract the maximum fault information.

In addition to nonlinear and nonstationary characteristics, the vibration signals of wheelset bearings are coupled with subgrade vibration impulses, transmission path coupling interference, and other factors. To solve these problems, this paper proposes a fault feature extraction method based on Teager energy spectrum correlation kurtosis (TSCK) guided parameter convex optimization tunable Q-factor wavelet transform (TQWT). The proposed method constructs the TSCK, which is particularly sensitive to the periodic impulse components of the wheelset bearing fault, as the objective function. Combining with the parameter convex optimization method, the key parameter combination of TQWT can be found, thus realizing the optimal sparse decomposition and feature extraction of the signal. The simulated and experimental analyses show that the method can extract the periodic impulse components of bearing fault in the presence of nonperiodic transient impulses, harmonic, and noise interference. The outline of this paper is as follows. Section “Theoretical framework of TSCK guided parameter convex optimization TQWT” introduces the theory and process of the proposed algorithm. Section “Simulated analysis” verifies the diagnosis effect of the proposed algorithm on bearings under different working conditions.

Theoretical framework of TSCK guided parameter convex optimization TQWT

Tunable Q-factor wavelet transform

TQWT is a wavelet transform method designed from the perspective of frequency domain filtering. ²⁶ Its implementation process depends on a decomposition filter bank and a reconstruction filter bank. Compared with the constant Q-factor of traditional wavelet transform, TQWT selects the quality factor according to the oscillation characteristics of the signal to be analyzed. The purpose of signal decomposition and reconstruction can be achieved by iteratively applying the two dual channel filter banks. The three-stage decomposition and reconstruction is taken as an example, and the corresponding dual channel filter banks are shown in Figure 1.

Figure 1.

The corresponding dual channel filter banks of (a) decomposition process and (b) reconstruction process.

${\mathrm{H}}_1\left(\omega \right)$ and ${\mathrm{H}}_0\left(\omega \right)$ are the frequency responses of low-pass and high-pass filters, respectively, and ${\mathrm{H}}_1^{*}\left(\omega \right)$ and ${\mathrm{H}}_0^{*}\left(\omega \right)$ are the complex conjugate of ${\mathrm{H}}_1\left(\omega \right)$ and ${\mathrm{H}}_0\left(\omega \right)$ , respectively. ${\mathrm{c}}^j$ and ${\mathrm{w}}^j$ are the scale coefficients and wavelet coefficients obtained from the jth layer wavelet transform, respectively.

Low-pass filter response and high-pass filter response in TQWT are defined by Equations (1) and (2), respectively.

$${\mathrm{H}}_0\left(\omega \right)=\{\begin{array}{cc}1,\hfill & \left|\omega \right|\le (1-\beta )\hfill \\ \multicolumn{1}{c}{\theta \left[\frac{\omega +(\beta -1)\pi }{\alpha +\beta -1}\right],}& (1-\beta )\pi \le \omega <\alpha \pi \hfill \\ \multicolumn{1}{c}{0,}& \alpha \pi \le \left|\omega \right|<\pi \hfill \end{array}$$

(1)

$${\mathrm{H}}_1\left(\omega \right)=\{\begin{array}{cc}0,\hfill & \left|\omega \right|\le (1-\beta )\hfill \\ \multicolumn{1}{c}{\theta \left[\frac{\alpha \pi -\omega }{\alpha +\beta -1}\right],}& (1-\beta )\pi \le \omega <\alpha \pi \hfill \\ \multicolumn{1}{c}{1,}& \alpha \pi \le \left|\omega \right|<\pi \hfill \end{array}$$

(2)

where $\alpha$ is the low-pass scale factor, $\beta$ is the high-pass scale factor, $\theta$ is the transition interval design function, ${\mathrm{H}}_1\left(\omega \right)$ and ${\mathrm{H}}_0\left(\omega \right)$ are the frequency responses of low-pass and high-pass filters, respectively. The maximum number of decomposition layers of TQWT can be obtained by Equation (3).

$$J_{max}=\lfloor \frac{\log (\beta \mathrm{Ns}/8)}{\log (1/\alpha )}\rfloor$$

(3)

where N is the length of the signal to be analyzed. The scale factors $\alpha$ and $\beta$ can be substituted into Equation (3) to obtain Equation (4).

$$J_{max}=\lfloor \frac{\log (Ns/4(Q+1))}{\log [(Q+1)/(Q+1-2/r)]}\rfloor$$

(4)

Teager energy spectrum correlation kurtosis

Teager energy operator (TEO) is a nonlinear operator, which amplifies the instantaneous energy component of the signal through nonlinear combination of the instantaneous value and its differential of the vibration signal, thus enhancing the impact characteristics of the signal. ²⁷ It is suitable for detecting the impact component in the signal, and the envelope demodulation effect of TEO has higher time resolution than Hilbert transform.

For a continuous-time signal $x\left(t\right)$ , the calculation process of the TEO can be expressed by Equation (5).

$$\psi \left[x\left(t\right)\right]={\left[\stackrel{·}{x}\left(t\right)\right]}^2-x\left(t\right)\stackrel{··}{x}\left(t\right)$$

(5)

where $\stackrel{·}{x}\left(t\right)$ and $\stackrel{··}{\mathrm{x}}\left(t\right)$ are the first and the second derivatives of signal x relative to time t. The amplitude modulation and frequency modulation process can be represented by Equation (6).

$$s\left(t\right)=a\left(t\right)\cos \left[\phi \left(t\right)\right]$$

(6)

where $a\left(t\right)$ and $\phi \left(t\right)$ are the instantaneous amplitude and instantaneous phase of $s\left(t\right)$ , respectively. The instantaneous frequency of the signal $s\left(t\right)$ is $\omega =\stackrel{·}{\phi }\left(t\right)$ . The energy operator of $s\left(t\right)$ can be expressed as Equation (7).

$$\begin{array}{c}\multicolumn{1}{c}{\psi \left[s\left(t\right)\right]={\left[a\left(t\right)\stackrel{·}{\phi }\left(t\right)\right]}^2+a^2\left(t\right)\stackrel{··}{\phi }\left(t\right)\times \sin \left[2\phi \left(t\right)\right]/2+{\cos }^2\left[\phi \left(t\right)\right]}\\ \multicolumn{1}{c}{[{\stackrel{·}{a}}^2\left(t\right)-a\left(t\right)\stackrel{··}{a}\left(t\right)]}\end{array}$$

(7)

The change of the carrier signal is faster than that of the modulated signal, so the instantaneous amplitude and frequency of the modulated signal are relatively slow compared with the high-frequency carrier, and can be approximately equal to a constant. Therefore, Equation (8) can be obtained.

$$\{\begin{array}{c}\stackrel{·}{a}\left(t\right)\approx 0\hfill \\ \multicolumn{1}{c}{\stackrel{··}{a}\left(t\right)\approx 0}\\ \multicolumn{1}{c}{\stackrel{··}{\phi }\left(t\right)\approx 0}\end{array}$$

(8)

Equation (7) can be simplified to Equation (9).

$$\psi \left[s\left(t\right)\right]\approx {\left[a\left(t\right)\stackrel{·}{\phi }\left(t\right)\right]}^2=a^2\left(t\right){\omega }^2\left(t\right)$$

(9)

The energy operator of the first-order differential of $\stackrel{·}{s}\left(t\right)$ can be expressed as Equation (10).

$$\psi \left[\stackrel{·}{s}\left(t\right)\right]\approx a^2\left(t\right){\omega }^4\left(t\right)$$

(10)

Equations (11) and (12) can be obtained from Equation (10).

$${\sin }^2\left(\omega \right)=\frac{1-\cos \left(2\omega \right)}{2}=\frac{\psi \left[\stackrel{·}{s}\left(n\right)\right]}{\psi \left[s\left(n\right)\right]}$$

(11)

$$a{\left(t\right)}^2=\frac{{\psi }^2\left[\stackrel{·}{s}\left(n\right)\right]}{\psi \left[s\left(n\right)\right]}$$

(12)

Note that if $x\left(n\right)$ has been obtained from sampling a continuous-time signal, then the actual frequency $f\left(n\right)$ (in Hz) can be found from $\omega$ as $f=\omega /2\pi$ . The instantaneous frequency and amplitude of signal $s\left(t\right)$ can be obtained, respectively, by Equations (13) and (14).

$$f\left(n\right)=\frac{1}{4\pi }\arccos \{1-\frac{\psi [s(n+1)-s(n-1)]}{2\psi \left[s\left(n\right)\right]}\}$$

(13)

$$a\left(n\right)=\frac{2\psi \left[s\left(n\right)\right]}{\sqrt{\psi [s(n+1)-s(n-1)]}}$$

(14)

Furthermore, Fourier transform is used to process the signal processed by TEO, as shown in Equation (15).

$$F\left(\psi \right)=\underset{0}{\overset{T}{\int }}\psi \left(s\left(t\right)\right)e^{-\mathrm{jwt}}\mathrm{d}t$$

(15)

The autocorrelation coefficient describes the relationship of the same signal at different times, which can be expressed by Equation (16).

$$R_x\left(\tau \right)=\frac{1}{T}\underset{T\to \infty }{lim}\underset{0}{\overset{T}{\int }}x\left(t\right)x(t+\tau )\mathrm{d}t$$

(16)

where $\tau$ reflects the time difference between $x\left(t\right)$ and $x(t+\tau )$ .

The independent variable frequency of Teager energy spectrum is regarded as a generalized time variable. Combine Equation (15) with Equation (16) to obtain the autocorrelation function of Teager energy spectrum, as shown in Equation (17).

$$R_F\left(\tau \right)=\frac{1}{T}\underset{T\to \infty }{lim}\underset{0}{\overset{T}{\int }}F\left(f\right)F(f+\tau )\mathrm{d}t$$

(17)

The process of Equation (17) can enhance the frequency components with coupling relationship in Teager energy spectrum and suppress the frequency components without coupling relationship.

The autocorrelation function of Teager energy spectrum often contains a large trend term, which can be removed by polynomial fitting. The autocorrelation function R in Equation (17) can be expressed by higher-order polynomial as shown in Equation (18).

$${\stackrel{\wedge }{R}}_k=a_0+a_{1}k+a_2{k}^2+···+a_m{k}^m(k=1,2,3,...,n)$$

(18)

where R is the data to be fitted, a is the undetermined coefficient, m is the order, and n is the data length

All sampled data can be calculated as shown in matrix relationship Equation (19).

$$\mathrm{\Delta }=\mathrm{KA}$$

(19)

where $\mathrm{\Delta }={[{\stackrel{\wedge }{R}}_1{\stackrel{\wedge }{R}}_2···{\stackrel{\wedge }{R}}_n]}^T$ , $K={\left[\begin{array}{cccc}1\hfill & 1^2\hfill & ···\hfill & 1^m\hfill \\ \multicolumn{1}{c}{1}& 2^2\hfill & ···\hfill & 2^m\hfill \\ \multicolumn{1}{c}{···}& ···\hfill & ···\hfill & ···\hfill \\ \multicolumn{1}{c}{1}& n^2\hfill & ···\hfill & n^m\hfill \end{array}\right]}^T$ , $A={[a_0{a}_1···a_m]}^T$

The undetermined coefficient $a_i(i=0,1,...,m)$ of function ${\stackrel{\wedge }{R}}_k$ is obtained to minimize the sum of squares of errors between function ${\stackrel{\wedge }{R}}_k$ and sampling data $R_k$ , which can be expressed as Equations (20) and (21).

$$min\left(E\right)=min\left[{(\stackrel{\wedge }{R}-R)}^T(\stackrel{\wedge }{R}-R)\right]$$

(20)

$$\begin{array}{c}\multicolumn{1}{c}{E={(\stackrel{\wedge }{R}-R)}^T(\stackrel{\wedge }{R}-R)={(\mathrm{KA}-R)}^T(\mathrm{KA}-R)}\\ \multicolumn{1}{c}{=A^T{K}^T\mathrm{KA}-A^T{K}^{T}R-\mathrm{AK}{R}^T+R^{T}R}\end{array}$$

(21)

Calculate the partial derivative of E to $a_i(i=0,1,...,m)$ , as shown in Equation (22).

$$\frac{\partial E}{\partial A}=2(K^T\mathrm{KA}-K^{T}R)$$

(22)

When the partial derivative is 0, the sum of squares of the errors is taken as the minimum, and the polynomial coefficients can be expressed by Equation (23).

$$A={\left(K^{T}K\right)}^{-1}{K}^{T}R$$

(23)

The result after removing the trend term can be expressed as shown in Equation (24).

$$Y=R_k-{\stackrel{\wedge }{R}}_k=R_k-(a_0+a_{1}k+a_2{k}^2+···+a_m{k}^m)$$

(24)

In order to verify the robustness of the proposed TSCK to noise compared with Hilbert envelope spectrum correlation kurtosis (HSCK), Hilbert envelope spectrum dispersion entropy (HSDE), Dispersion entropy (DE), and Shannon entropy (SE),^28–31 different levels of noise are added to the periodic pulse analog signal to simulate different signal-to-noise ratio (SNR) states. The mathematical model of the analog signal can be expressed by Equation (25). ³²

$$S_n\left(t\right)=\sum _{i=1}^{n}A{e}^{-\xi {[t-q_i(t)/f_{\mathrm{outer}}]}^2}·\sin (2\pi f_{\mathrm{o}}t)+n(t)$$

(25)

where A is the amplitude of periodic signal, $\xi$ is the damping of the periodic signal, $f_{\mathrm{o}}$ is the natural frequency of the system, $f_{\mathrm{outer}}$ is the characteristic frequency, and $n(t)$ is noise of different degrees (The SNR range is (−18 to −4 dB)).

The time-domain waveform, envelope spectrum, and Teager energy spectrum of the simulated signals are shown in Figure 2(a)–(c). The variation trend of five different indicators under different SNR is shown in Figure 2(d). It can be seen from the figure that the slope of DE and HSDE curves is low, indicating low sensitivity and resolution to different SNR. SE has good sensitivity and resolution at high SNR stage (greater than −8 dB), but has poor stability at low SNR stage. TSCK and HSCK have good sensitivity and stability, and TSCK is better than HSCK.

Figure 2.

Simulated signal under different SNR: (a) time-domain waveform, (b) envelope spectrum, (c) Teager energy spectrum, and (d) variation trend of five different indicators under different SNR.

The greater the Q, the stronger the oscillation. The selection of the Q-factor determines the selection of the final filtering band, and the appropriate band evaluation index is the key to ensure the final filtering effect. The parameter r is the redundancy, which determines the coincidence degree between wavelet frequency responses. The greater the r, the higher the coincidence degree of frequency responses. Therefore, reasonable selection of parameter combination Q and r of TQWT is of great significance for effective separation of the periodic impulse components of wheelset bearing faults. Considering that TSCK has the special sensitivity to periodic impulse components and the ability to suppress nonperiodic transient impulse components and harmonic interference, TSCK is constructed as an objective function, and the parameter optimization of TQWT is realized by the parameter convex optimization method. The process is that the signal to be processed is decomposed by cyclic TQWT with different parameter combinations, the TSCKs of the subband after each decomposition is calculated, the maximum TSCK value is selected to represent the decomposition result, the characteristic matrix is constructed with the maximum TSCK values of each result of cyclic TQWT decomposition. Then, parameter convex optimization is used to obtain discrete points in the connected space, and 3D interpolation smooth fitting is performed to obtain the distribution of the maximum TSCK value under different parameter combinations. Finally, the distribution of the maximum TSCK values is projected to the Q and r parameter planes and meshed to obtain the optimal TQWT parameter combination at the maximum point.

The proposed diagnosis scheme

The flowchart of the proposed method is shown in Figure 3.

1. Parametric TQWT ([Q r]) is used to decompose the original signal, in which parameters Q and r are selected within a specific interval such as [1 10].

2. Each subband after TQWT decomposition is demodulated by the TEO, and the spectrum of each demodulated signal is analyzed to obtain the Teager energy spectrum.

3. The spectral independent variable of Teager energy spectrum is regarded as a generalized time quantity for autocorrelation analysis, and the autocorrelation function is obtained.

4. The trend term of autocorrelation coefficient of Teager energy spectrum is removed by polynomial fitting method to obtain a smooth autocorrelation coefficient sequence, and the kurtosis of the autocorrelation coefficient is calculated.

5. Combining steps (2), (3), and (4), the TSCK of each subband signal under parametric combination ([Q r]) can be obtained, and the optimal TQWT parameter combination ([Q^*r^*]) can be obtained by using the parameter convex optimization visualization method, so as to obtain the optimal decomposition mode.

6. The optimal subband interval is obtained under the optimal TQWT decomposition mode, and the fault feature frequency is extracted to diagnose the fault type.

Figure 3.

Flowchart of the proposed algorithm.

Simulated analysis

When monitoring the running state of the wheelset bearing, the sensors are generally arranged on the wheelset bearing seat or box. The signal components received by the sensors may be complex. In addition to the periodic impulse components of the bearing fault, they may also contain a variety of interference components, such as the nonperiodic transient impulse components introduced by the subgrade vibration, the noise components and harmonic components introduced by the transmission path. The schematic diagram of the generation, transmission, and superposition of each signal component is shown in Figure 4.

Figure 4.

Schematic diagram of generation, transmission, and superposition of each signal components.

The mathematical model shown in Equation (26) is used to construct the simulated signal, ³³ which mainly includes three parts, namely, periodic fault impulse components $\left(y_1\right)$ , nonperiodic transient impulse components $\left(y_2\right)$ and harmonic components $\left(y_3\right)$ and noise components $\left(n\left(t\right)\right)$ .

$$\{\begin{array}{c}{y}_1=\sum _{i=1}^{n}A{e}^{-\xi {[t-q_i(t)/f_{\mathrm{outer}}]}^2}·\sin (2\pi f_{o}t)\hfill \\ \multicolumn{1}{c}{y_2=A_1{e}^{[-b_1\times (t-t_1)]}·\sin [2\pi f_1·(t-t_1)]}\\ \multicolumn{1}{c}{+A_2{e}^{[-b_2\times (t-t_2)]}·\sin [2\pi f_2·(t-t_2)]}\\ \multicolumn{1}{c}{y_3=A_3·\sin (2\pi f_{\mathrm{mesh}}t+{\beta }_1)(1+A_4\sin (2\pi f_{h}t+{\beta }_2))}\\ \multicolumn{1}{c}{y_s=y_1+y_2+y_3+n(t)}\end{array}$$

(26)

where A, $\xi$ , $f_{\mathrm{outer}}$ , and $f_{\mathrm{o}}$ are the vibration amplitude, damping ratio, fault characteristic frequency, and natural frequency of periodic fault impulse components, respectively. $A_1/A_2$ , $f_1/f_2$ , and $b_1/b_2$ are the vibration amplitude, natural frequency, and damping ratio of nonperiodic transient impulse components, respectively. $A_3$ , $f_{\mathrm{mesh}}$ , and $f_h$ are the amplitude, carrier frequency, and modulation frequency of harmonic signal, respectively. The parameters of the simulated signal are listed in Table 1. The time-domain waveform of the simulated signals is shown in Figure 5.

Table 1.

Parameters of the simulated signal.

Parameter type	Parameter values (Hz)
Sampling frequency	8192
Natural frequency fo	2100
Natural frequency f1 = f2	3000
Carrier frequency fmesh	360
Modulation frequency fh	15
Fault frequency fouter	50

Figure 5.

Time-domain waveform of the simulated signals.

The calculation process of Kurtogram defined in the time domain is to decompose the original signal into a series of subband signals through tree filtering, and then calculate the kurtosis value of the time-domain sequence of each subband, and normalize it to obtain the kurtosis distribution of the subbands. The optimal subband determined by this method is located in the twelfth subband of the fourth layer (concentrated near 3000 Hz) (Figure 6(a)). The time-domain waveform and envelope spectrum of this subband signal are shown in Figure 6(b) and (c), respectively. The signal components of this subband are mainly nonperiodic transient impulse components, rather than the periodic impulse components to be extracted. The calculation process of Kurtogram defined in the frequency domain is to analyze the envelope spectrum of each subband after the tree decomposition, and then calculate the kurtosis value of the envelope spectrum sequence separately, and normalize it to obtain the kurtosis distribution of the subbands. The optimal subband located by this method is located in the second subband of the fourth layer (concentrated near 360 Hz) (Figure 7(a)). The time-domain waveform and envelope spectrum of this subband are shown in Figure 7(b) and (c), respectively. The signal components of this subband are mainly harmonic signal components. It can be seen that the Kurtogram algorithm defined by the frequency domain index will be affected by the harmonic components (which can be regarded as nonperiodic transient impulse components in the frequency domain), so it is unable to effectively identify the fault components.

Figure 6.

Analysis results of the simulated signals using the Kurtogram algorithm defined in time domain: (a) the optimal band to be located, (b) waveform of the selected subband, and (c) envelope spectrum of (b).

Figure 7.

Analysis results of the simulated signals using the Kurtogram algorithm defined in frequency domain: (a) the optimal band to be located, (b) waveform of the selected subband, and (c) envelope spectrum of (b).

The method proposed in this paper is used for analysis. The process is that the simulated signals is decomposed by cyclic TQWT with different parameter combinations (the double parameter optimization interval of Q and r is 1 to 10), the TSCKs of the subband after each decomposition is calculated, the maximum TSCK value is selected to represent the decomposition result, and the characteristic matrix is constructed with the maximum TSCK values of each result of cyclic TQWT decomposition. Then, parameter convex optimization is used to obtain discrete points in the connected space, and 3D interpolation smooth fitting is performed to obtain the distribution of the maximum TSCK value under different parameter combinations. Finally, the distribution of the maximum TSCK values is projected to the Q and r parameter planes and meshed to obtain the optimal TQWT parameter combination at the maximum point. The results are shown in Figure 8. Figure 8(a) shows the maximum TSCK discrete data distribution of each subband after TQWT decomposition under different [Q r] parameter combinations. Figure 8(b) shows the result of 3D interpolation smooth fitting. Figure 8(c) shows the projection of Figure 8(b) on the 2D dimensional plane, and Figure 8(d) shows the grid division results of the parameter plane for the convenience of parameter selection. From Figure 8(d), it can be clearly seen that the optimal combination of TQWT decomposition parameters [Q r] that make the maximum value of TSCK index region is [8 2].

Figure 8.

Analysis results of the simulated signals by the proposed method: (a) the maximum TSCK distribution of each subband after TQWT decomposition under different [Q r], (b) the result of 3D interpolation smooth fitting, (c) the projection on the 2D dimension plane of (b), and (d) grid generation result of parameter plane.

When the TQWT parameter combination is [8 2], the TSCK value of the sixth subband signal is the largest, and its frequency response results are shown in Figure 9(a). The frequency corresponding to the selected subband is concentrated around 2100 Hz, that is, the subband range of the periodic fault impulse components. The time-domain waveform and envelope spectrum of this subband are shown in Figure 9(b) and (c). There are obvious periodic impulse components in the time-domain waveform. The fault characteristic frequency and its multiple frequency can be seen in the envelope spectrum. Therefore, the proposed method can accurately extract the periodic impulse components.

Figure 9.

(a) Frequency responses with parameter combination [8 2], (b) time-domain waveform of the selected subband, and (c) envelope spectrum of (b).

Experimental analysis

Case 1: Experimental signals with nonperiodic transient impulse components

In order to further verify the effectiveness of the proposed method, experimental fault signals are used for analysis. The experimental equipment is the QPZZ-II simulated experimental platform. The structure of the test-bed is shown in Figure 10(a). The model of the bearing was 6205E, and electrical discharge machining is used to process faults on the inner, as shown in Figure 10(b). The entire test-bed is driven by a drive motor $(N=1470\mathrm{rpm})$ to drive the coupling, shaft, and bearing to rotate. The left bearing seat is installed with normal bearings, the right bearing seat is installed with faulty bearings, the long shaft is installed with eddy current displacement sensors, and the right bearing seat is installed with vibration acceleration sensors in the horizontal and vertical directions.

Figure 10.

(a) Bearing fault experimental equipment of QPZZ-II, and (b) 6205E bearing with inner ring fault.

The theoretical fault characteristic frequency (as shown in Table 2) of the bearing can be calculated by using Equation (27) based on the bearing structural dimension parameters (as shown in Table 3) and rotating speed (1470 rpm).

$$\{\begin{array}{c}{f}_i=\frac{Z}{2}(1+\frac{d}{D}\cos \alpha )\frac{N}{60}=132.6\mathrm{Hz}\hfill \\ \multicolumn{1}{c}{f_{\mathrm{o}}=\frac{Z}{2}(1-\frac{d}{D}\cos \alpha )\frac{N}{60}=87.9\mathrm{Hz}}\\ \multicolumn{1}{c}{f_{\mathrm{ball}}=\frac{D}{d}(1-{\left(\frac{d}{D}\right)}^2{\cos }^2\alpha )\frac{N}{60}=116\mathrm{Hz}}\end{array}$$

(27)

Table 2.

Theoretical fault characteristic frequencies of 6205E bearing.

Parameter type	Parameter values
Inner ring fault (fi)	132.6
Outer ring fault (fo)	87.9
Ball fault (fball)	116

Table 3.

Structural dimension parameters of 6205E bearing.

Parameter type (6205E)	Parameter values
Number of balls (Z)	9
Ball diameter (d) (mm)	7.9
Pitch diameter (D) (mm)	39
Contact angle (α) (°)	0

In the first group of experiments, the motor speed is 1470 rpm, and the sampling frequency is 12,800 Hz. The fault type is inner race fault, and the corresponding bearing structure parameters are shown in Table 3. The time-domain waveform of the experimental signals is shown in Figure 11. It can be seen that there is a large amplitude of nonperiodic transient impulse components between 0.3 and 0.4 s of the time-domain waveform.

Figure 11.

Time-domain waveform of the first group of experimental signals.

Kurtogram algorithm defined by kurtosis index in time domain is used to analyze the first group of experimental signals, and the results are shown in Figure 12. The optimal location area is located in the second subband of the third layer. The time-domain waveform and envelope spectrum of the subband signals are shown in Figure 12(b) and (c), respectively. The signal component of this subband is mainly the nonperiodic transient impulse components, rather than periodic impulse components caused by bearing failure in the experimental signals.

Figure 12.

Analysis results of the first group experimental signals using the Kurtogram algorithm defined in time domain: (a) the optimal band to be located, (b) waveform of the selected subband, and (c) envelope spectrum of (b).

The method proposed in this paper is used for analysis. TQWT is used to decompose the first group of experimental signals, and the TSCK of the decomposed subband signal is used to convex optimize the parameters Q and r of TQWT to find the optimal parameter combination. The double parameter optimization interval of Q and r is 1 to 10. The results are shown in Figure 13. Figure 13(a) shows the maximum TSCK discrete data distribution of each subband after TQWT decomposition under different [Q r] parameter combinations. Figure 13(b) shows the result of 3D interpolation smooth fitting. Figure 13(c) shows the projection of Figure 13(b) on the 2D dimensional plane, and Figure 13(d) shows the grid division results of the parameter plane for the convenience of parameter selection. From Figure 13(d), it can be clearly seen that the optimal combination of TQWT decomposition parameters [Q r] that make the maximum value of TSCK index region is [10 6].

Figure 13.

Analysis results of the first group of experimental signals by the proposed method: (a) the maximum TSCK distribution of each subband after TQWT decomposition under different [Q r], (b) the result of 3D interpolation smooth fitting, (c) the projection on the 2D dimension plane of (b), and (d) grid generation result of parameter plane.

When the TQWT parameter combination is [10 6], the TSCK value of the fifth subband signal is the largest, and its frequency band response results are shown in Figure 14(a). The time-domain waveform and envelope spectrum of this subband are shown in Figure 14(b) and (c). There are obvious periodic impulse components in the time-domain waveform. The fault characteristic frequency and its multiple frequency can be seen in the envelope spectrum. Therefore, the proposed method can accurately extract the periodic impulse components.

Figure 14.

(a) Frequency responses with parameter combination [10 6], (b) time-domain waveform of the selected subband, and (c) envelope spectrum of (b).

Case 2: Experimental signals with harmonic components

The second group of experiments is also carried out on the test bench in Figure 10(a). The experimental bearing is replaced with N205. The theoretical fault characteristic frequency (as shown in Table 4) of the bearing can be calculated based on the bearing structural dimension parameters (as shown in Table 5) and rotating speed. The fault type is outer ring fault (Figure 15). The motor speed was 1440 rpm, and the sampling frequency is 12,800 Hz. The time-domain waveform of the experimental signals is shown in Figure 16, from which a strong harmonic modulation component can be seen.

Table 4.

Theoretical fault characteristic frequencies of N205 bearing.

Parameter type	Parameter values
Inner ring fault (fi)	172.1
Outer ring fault (fo)	115.9
Ball fault (fball)	118.5

Table 5.

Structural dimension parameters of N205 bearing.

Parameter type (N205)	Parameter values
Number of balls	12
Ball diameter (mm)	7.5
Pitch diameter (mm)	38.5
Contact angle (°)	0

Figure 15.

N205 bearing with outer ring fault.

Figure 16.

Time-domain waveform of the second group of experimental signals.

Kurtogram algorithm defined by kurtosis index of frequency domain envelope spectrum is used to analyze the second group of experimental signals, and the results are shown in Figure 17. The optimal subband located by this method is the first subband of the fourth layer. The time-domain waveform and envelope spectrum of this subband are shown in Figure 17(b) and (c), respectively. The signal components of this subband are mainly harmonic signal components, which cannot effectively extract the fault characteristic frequency.

Figure 17.

Analysis results of the second group experimental signals using the Kurtogram algorithm defined in frequency domain: (a) the optimal band to be located, (b) waveform of the selected subband, and (c) envelope spectrum of (b).

The second group of experimental signals is analyzed using the method proposed in this paper. TSCK is used to guide the decomposition parameters Q and r of TQWT to obtain the optimal combination of signal decomposition parameters. The parameter optimization range of Q and r is 1 to 10, and the optimization results are shown in Figure 18. Figure 18(a) shows the maximum TSCK discrete data distribution of each subband after TQWT decomposition under different [Q r] parameter combinations. Figure 18(b) shows the result of 3D interpolation smooth fitting. Figure 18(c) shows the projection of Figure 18(b) on the 2D dimensional plane, and Figure 18(d) shows the grid division results of the parameter plane for the convenience of parameter selection. From Figure 18(d), it can be clearly seen that the optimal combination of TQWT decomposition parameters [Q r] that make the maximum value of TSCK index region is [4 7].

Figure 18.

Analysis results of the second group of experimental signals by the proposed method: (a) the maximum TSCK distribution of each subband after TQWT decomposition under different [Q r], (b) the result of 3D interpolation smooth fitting, (c) the projection on the 2D dimension plane of (b), and (d) grid generation result of parameter plane.

The TQWT with the combination of Q and r values as [4 7] decomposes the second group of experimental signals. Its frequency response results are shown in Figure 19(a). The TSCK value of the fifth subband is the largest. The time-domain waveform and envelope spectrum of the subband are shown in Figure 19(b) and (c). The time-domain signals contain periodic impulse components, and the fault characteristic frequency and its double frequency components can be extracted from the envelope spectrum.

Figure 19.

(a) Frequency responses with parameter combination [4 7], (b) time-domain waveform of the selected subband, and (c) envelope spectrum of (b).

Signal analysis of wheelset bearing experiment

In order to verify the effectiveness of the proposed method on the wheelset bearing experiment signals, the experiment is carried out on the experimental platform as shown in Figure 20(a). The bearing with inner ring failure is shown in Figure 20(b). Vibration acceleration sensor is installed on the wheelset bearing seat. The sampling frequency is 25.6 kHz. The shaft rotational speed is set to 465 rpm ( $f_r=7.8\mathrm{Hz}$ ). The structural dimension parameters of the bearing are shown in Table 6. The theoretical fault characteristic frequencies of the bearing (Table 7) can be calculated from the above bearing parameters and working conditions. The fault characteristic frequencies of bearing inner ring is $f_i=88.2\mathrm{Hz}$ . The time-domain waveform of wheelset bearing inner race fault is shown in Figure 21.

Figure 20.

(a) Wheelset experiment platform, and (b) bearing with inner ring fault.

Table 6.

Structural dimension parameters of the wheelset bearing.

Parameter type	Parameter values
Number of balls	20
Ball diameter (mm)	24.74
Pitch diameter (mm)	176.29
Contact angle (deg)	8.833

Table 7.

Theoretical fault characteristic frequencies of the wheelset bearing.

Parameter type	Parameter values
Inner ring fault	88.2
Outer ring fault	66.8
Ball fault	27.1

Figure 21.

Time-domain waveform of the wheelset bearing signals.

The method proposed in this paper is used to analyze the experimental signal of wheelset bearing. TSCK is used to guide the decomposition parameters Q and r of TQWT to obtain the optimal combination of signal decomposition parameters. The parameter optimization range of Q and r is 1 to 10, and the optimization results are shown in Figure 22. Figure 22(a) shows the maximum TSCK discrete data distribution of each subband after TQWT decomposition under different [Qr] parameter combinations. Figure 22(b) shows the result of 3D interpolation smooth fitting. Figure 22(c) shows the projection of Figure 22(b) on the 2D dimensional plane, and Figure 22(d) shows the grid division results of the parameter plane for the convenience of parameter selection. From Figure 22(d), it can be clearly seen that the optimal combination of TQWT decomposition parameters [Qr] that make the maximum value of TSCK index region is [10 2].

Figure 22.

Analysis results of the wheelset bearing signals by the proposed method: (a) the maximum TSCK distribution of each subband after TQWT decomposition under different [Q r], (b) the result of 3D interpolation smooth fitting, (c) the projection on the 2D dimension plane of (b), and (d) grid generation result of parameter plane.

The combination of Q and r parameters of TQWT is set as [10 2] to analyze the experimental signal of wheelset bearing. Its frequency response results are shown in Figure 23(a). The TSCK value of the sixth subband is the largest. The time-domain waveform and envelope spectrum of the subband are shown in Figure 23(b) and (c). The time-domain signals contain periodic impulse components, and the fault characteristic frequency and its double frequency components can be extracted from the envelope spectrum.

Figure 23.

(a) Frequency responses with parameter combination [10 2], (b) time-domain waveform of the selected subband, and (c) envelope spectrum of (b).

Conclusions

The running state of the wheelset bearing is related to the running safety of the vehicle, and the vibration signal is a typical medium in the wheelset bearing monitoring. To solve the problem that it is difficult to extract the fault features of wheel bearing vibration signals due to the coupling interference of nonlinear, nonstationary, subgrade vibration and transmission path, this paper proposes a fault feature extraction method based on TSCK guided parameter convex optimization TQWT. This method constructs the TSCK which is particularly sensitive to the periodic impulse components of wheelset bearing fault as the objective function and combines with the parameter convex optimization method to find the key parameter combination (Q^*r^*) of TQWT, so as to achieve the goal of optimal sparse decomposition and feature extraction of the signals. Through the analysis of simulated signals and experimental signals, the advantages of the proposed method in dealing with the interference of nonperiodic transient impulse components, harmonic components, and noise components are verified.