1 Introduction

Wheel out-of-roundness, especially wheel polygonal wear, is a very common type of wheel defect for a railway vehicle especially for a locomotive [1,2,3,4]. Wheel polygonization is a phenomenon that the wheel profile periodically deviates from an ideal round shape (i.e., periodic radial deviation) along circumference [5]. The defect will seriously deteriorate wheel-rail contact behavior, increase wheel-rail vertical force and even cause strong impact load to locomotive components and infrastructure [6, 7]. As a result, wheel polygonization accelerates failure of the locomotive and track components and thus seriously threatens the locomotive running safety and reliability [6]. In particular, with increasing in the running speed and axle load of modern rail transit equipment, the harms of the wheel polygonization behave more obvious. Therefore, researches of the wheel polygonal wear have always been hot topics in railway engineering.

In the past decades, most of studies have been carried out on formation mechanism of the polygonal wear. For instance, Jin et al. [8] investigated the root causes of polygonal wear of metro train wheels and found that the formation of the ninth-order polygonal wear is related to the first bending resonance of the wheelset. Tao et al. [9] found the wheelset bending vibration excited by rail corrugation with a specific wavelength can cause higher-order (i.e., the 12th–14th orders) wheel polygonal wear for metro trains. There are also studies on the mechanism of wheel polygonal wear of electric locomotives [10] and high-speed trains [11], and most of the researches show that the polygonal wear is closely related to wheelset vibrations (bending or self-excited vibration). On the other hand, many studies focus on the influence of polygonal wear on vehicle and infrastructure components. For example, Wang et al. [12, 13] investigated the influence of the polygonal wear on dynamic performance of gearbox housing and axle-box bearing of high-speed trains. The results show that wheel polygonization can increase vibrations and forces of gearbox and bearing, and thus accelerates degeneration of the components. Chen et al. [14] studied the effect of the polygonal wear on subgrade performance by constructing a vehicle–track–subgrade coupled dynamics model, and showed that vibration and stress of the subgrade are mainly affected by lower-order polygonal wear. It is worth noting that some researchers have investigated vibration characteristics induced by polygonal wheels. For instance, Yang et al. [15] compared amplitudes, root-mean-square (RMS) and peak-to-valley values of axle box accelerations of heavy-haul locomotives with and without wheel polygonal wear, and the results indicate that the polygonization significantly increases the locomotive vibration.

Due to the severe harms of polygonal wear, condition monitoring and fault detection of wheel polygonization have aroused considerable attention. Song et al. [16] employed piezoelectric strain sensors installed on the rail to measure wheel-rail forces caused by a passing vehicle, and then identified polygonal wheels according to magnitudes of the impact forces. Similarly, Wei et al. [17] defined a wheel condition index based on track strain data measured by fiber Bragg grating sensors, and developed a real-time wheel condition monitoring system based on the proposed index. Given that the changes of the track strain caused by wheel minor defects may be unobvious and thus can be easily affected by wheel load effect, Liu et al. [18] employed a Bayesian blind source separation method to extract defect-sensitive components from raw strain data. To improve wheel fault detection performance, Ni et al. [19] extracted cumulative distribution function values and frequency points of strain data as features for training, and then applied a sparse Bayesian learning method to achieve wheel defect detection. Note that the above methods all rely on track-side sensor data collected in a certain railway section, and thus cannot monitor wheel conditions in the whole vehicle running lines. Of particular interest to us are fault detection methods using on-board sensor data such as vehicle vibration [20, 21]. For instance, Song et al. [22] identified wheel polygonization by detecting a passing frequency of a polygonal wheel from an axle-box acceleration, and they also investigated acceleration characteristics induced by wheel polygonal wear with different harmonic orders and amplitudes. This method assumes that the vehicle is running under a constant speed condition. Sun et al. [23] applied an angle domain synchronous averaging technique to enhance fault-related components of axle-box acceleration signals, and then constructed a health indicator to reflect roughness levels of wheel polygonal wear. However, this method requires a speed coder to measure the vehicle speed, which inevitably increases hardware costs. In addition, it cannot directly identify amplitudes of the polygonal wear.

A railway vehicle often suffers from non-stationary running conditions such as traction and braking, resulting in time-varying passing frequencies of a polygonal wheel. Traditional frequency spectrum-based fault detection methods cannot be applied for the time-varying fault feature extraction. Time–frequency (TF) analysis is a powerful tool to address this issue. Various TF transforms have been developed to generate TF representations (TFRs) of non-stationary signals such as short-time Fourier transform (STFT), continuous wavelet transform, TF reassignment and synchrosqueezing transform [24,25,26,27], some of which have already been widely applied in machine fault diagnosis [28,29,30]. However, most of the methods assume that a signal is quasi-stationary and thus cannot accurately analyze strongly modulated non-stationary signals. It is worth noting that adaptive signal decomposition methods like empirical mode decomposition, empirical wavelet transform, variational mode decomposition, can be combined with Hilbert transform to generate improved TFRs [31,32,33]. Nevertheless, these methods are mainly designed for narrow-band signals and thus may suffer from over-decomposition issue for wide-band non-stationary signals. To deal with this issue, Chen et al. developed an adaptive chirp mode decomposition (ACMD) framework and proposed various methods based on the unified framework [34,35,36,37,38]. The ACMD is able to accurately estimate instantaneous frequencies (IFs) and instantaneous amplitudes (IAs) of non-stationary signals and thus generates high-quality TFRs of signals. However, the output results of an ACMD are sensitive to inputted initial IFs and its applicability in wheel polygonization detection has not been investigated.

In summary, quantitative detection of harmonic orders and amplitudes of railway wheel polygonal wear using vehicle vibration accelerations is a challenging task which has not been properly solved by existing methods. In particular, non-stationary running conditions of a vehicle make the issue much more difficult. Utilizing the good advantages of the ACMD for non-stationary signal processing, a novel method for quantitative fault detection of locomotive wheel polygonization under non-stationary conditions is proposed in this paper. Firstly, the ACMD is combined with a TF ridge detection method to accurately estimate a gear meshing frequency from motor vibration acceleration signal and therefore a wheel rotating frequency can be obtained according to gear transmission relation. At this stage, the IF initialization issue of ACMD is addressed by TF ridge detection. Then, a signal resampling technique is applied to vibration acceleration of axle box to remove the non-stationary effect of the signal and therefore harmonic orders of polygonal wear can be accurately identified by order analysis. Finally, some polygonal-wear-related harmonics of the axle-box acceleration signal are extracted by ACMD and then an inertial algorithm is applied to the extracted harmonics to estimate amplitudes of the polygonal wear. The proposed method achieves the wheel rotating frequency estimation and fault detection only using locomotive vibration signals (i.e., free of speed coder).

The structure of the paper is organized as follows. In Sect. 2, the signal model is briefly described at first and then the principle of ACMD is introduced. Section 3 details the proposed wheel polygonization detection method. Effectiveness of the proposed method is demonstrated by dynamics simulation in Sect. 4 and field test in Sect. 5. Conclusions are drawn in Sect. 6.

2 Theoretical background

2.1 Signal model

When a locomotive undergoes non-stationary running conditions, its vibration responses will exhibit time-varying amplitude modulated (AM) and frequency modulated (FM) features. In addition, practical vibration responses contain multiple components induced by different excitation sources. Therefore, a multi-component AM-FM model is often employed to characterize the locomotive vibration signal as

$$g\left( t \right) = \sum\limits_{p = 1}^{P} {g_{p} \left( t \right)} = \sum\limits_{p = 1}^{P} {A_{p} \left( t \right)\cos \left[ {2{\uppi }\int_{0}^{t} {f_{p} \left( \upsilon \right){\text{d}}\upsilon } + \varphi_{p} } \right]},$$
(1)

where P denotes the number of signal components, and \(A_{p} \left( t \right)\), \(f_{p} \left( t \right)\) and \(\varphi_{p}\) stand for IA, IF and initial phase of the p-th signal component, respectively. It is worth noting that the vibration response induced by wheel polygonal wear often shows a harmonic structure where the IFs of the signal components are integer multiples of wheel rotating frequency. Namely, Eq. (1) can be rewritten as

$$g\left( t \right) = \sum\limits_{p = 1}^{P} {A_{p} \left( t \right)\cos \left[ {2{\uppi }o_{{p}} \int_{0}^{t} {f_{{\text{w}}} \left( \upsilon \right){\text{d}}\upsilon } + \varphi_{p} } \right]},$$
(2)

where \(f_{p} \left( t \right) = o_{p} f_{{\text{w}}} \left( t \right)\), \(f_{{\text{w}}} \left( t \right)\) is wheel rotating frequency, and \(o_{p}\) is an integer called harmonic order of the polygonal wear. It can be seen that, to effectively detect the wheel polygonization, one should accurately estimate \(f_{{\text{w}}} \left( t \right)\) in advance.

2.2 Adaptive chirp mode decomposition

The ACMD is designed to effectively estimate the IFs and IAs of a multi-component AM-FM signal (also called chirp signal) and thus achieves reconstruction and separation of the signal components. The ACMD can accurately analyze strongly time-varying modulated signals. The basic idea of the ACMD is to use a demodulation technique to reduce modulation degree of an AM-FM signal. Specifically, Eq. (1) can be rewritten into a frequency demodulated form as

$$g\left( t \right) = \sum\limits_{p = 1}^{P} {\left\{ {x_{p} \left( t \right)\cos \left[ {2{\uppi }\int_{0}^{t} {\tilde{f}_{p}^{{\text{d}}} \left( \upsilon \right)} {\text{d}}\upsilon } \right] + y_{p} \left( t \right)\sin \left[ {2{\uppi }\int_{0}^{t} {\tilde{f}_{p}^{{\text{d}}} \left( \upsilon \right)} {\text{d}}\upsilon } \right]} \right\}}$$
(3)

with

$$\left\{ \begin{gathered} x_{p} \left( t \right) = A_{p} \left( t \right)\cos \left\{ {2{\uppi }\int_{0}^{t} {\left[ {f_{p} \left( \upsilon \right) - \tilde{f}_{p}^{{\text{d}}} \left( \upsilon \right)} \right]} {\text{d}}\upsilon + \varphi_{p} } \right\} \hfill \\ y_{p} \left( t \right) = - A_{p} \left( t \right)\sin \left\{ {2{\uppi }\int_{0}^{t} {\left[ {f_{p} \left( \upsilon \right) - \tilde{f}_{p}^{{\text{d}}} \left( \upsilon \right)} \right]} {\text{d}}\upsilon + \varphi_{p} } \right\} \hfill \\ \end{gathered} \right. ,$$
(4)

where \(\tilde{f}_{p}^{{\text{d}}} \left( t \right)\) is a demodulation frequency, and \(x_{p} \left( t \right)\) and \(y_{p} \left( t \right)\) represent two demodulated signals. Obviously, if \(\tilde{f}_{p}^{{\text{d}}} \left( t \right) = f_{p} \left( t \right)\), the FM effect of the signals in Eq. (4) will be completely removed, resulting in purely AM signals with the narrowest bandwidth. Therefore, the ACMD estimates IFs and reconstructs signal components by minimizing bandwidths of the demodulated signals [34, 35] as

$$\begin{aligned} &\mathop {\min }\limits_{{\left\{ {x_{p} \left( t \right)} \right\},\left\{ {y_{p} \left( t \right)} \right\},\left\{ {\tilde{f}_{p}^{{\text{d}}} \left( t \right)} \right\}}} \left\{ {{\mathcal{L}}_{\tau } \left( {\left\{ {x_{p} \left( t \right)} \right\},\left\{ {y_{p} \left( t \right)} \right\},\left\{ {\tilde{f}_{p}^{{\text{d}}} \left( t \right)} \right\}} \right)} \right\} \hfill \\ &= \mathop {\min }\limits_{{\left\{ {x_{p} \left( t \right)} \right\},\left\{ {y_{p} \left( t \right)} \right\},\left\{ {\tilde{f}_{p}^{{\text{d}}} \left( t \right)} \right\}}} \left\{ {\sum\limits_{p = 1}^{P} {\left[ {\left\| {x_{p}^{\prime \prime } \left( t \right)} \right\|_{2}^{2} + \left\| {y_{p}^{\prime \prime } \left( t \right)} \right\|_{2}^{2} } \right]} + \tau \left\| {g\left( t \right) - \sum\limits_{p = 1}^{P} {g_{p} \left( t \right)} } \right\|_{2}^{2} } \right\}, \hfill \\ &{\text{ s}}{\text{.t}}{.}\,\,g_{p} \left( t \right) = x_{p} \left( t \right)\cos \left[ {2{\uppi }\int_{0}^{t} {\tilde{f}_{p}^{{\text{d}}} \left( \upsilon \right)} {\text{d}}\upsilon } \right] + y_{p} \left( t \right)\sin \left[ {2{\uppi }\int_{0}^{t} {\tilde{f}_{p}^{{\text{d}}} \left( \upsilon \right)} {\text{d}}\upsilon } \right], \hfill \\ \end{aligned}$$
(5)

where \(\left\| \cdot \right\|_{2}\) denotes \(l_{2}\) norm, \(\cdot^{\prime\prime}\) stands for the second derivative, \(\left\| {x_{p}^{\prime \prime } \left( t \right)} \right\|_{2}^{2}\) and \(\left\| {y_{p}^{\prime \prime } \left( t \right)} \right\|_{2}^{2}\) are used for assessing signal bandwidths, and \(\tau > 0\) is a weighting coefficient.

To efficiently solve the optimization problem in Eq. (5), an iterative algorithm which alternately updates the IFs and demodulated signals are available in [34, 35]. Main procedures of the algorithm include: (1) input initial IFs; (2) calculate demodulated signals with the inputted IFs; (3) update the IFs using phase information of the demodulated signals; (4) repeat steps (2) and (3) until the algorithm converges. Assuming that the finally estimated IFs and demodulated signals by ACMD are denoted by \(\tilde{f}_{p} \left( t \right)\) and \(\tilde{x}_{p} \left( t \right)\), \(\tilde{y}_{p} \left( t \right)\), for \(p = 1,2,\, \ldots ,\,P\), respectively, then each signal component can be reconstructed as

$$\tilde{g}_{p} \left( t \right) = \tilde{x}_{p} \left( t \right)\cos \left[ {2{\uppi }\int_{0}^{t} {\tilde{f}_{p} \left( \upsilon \right)} {\text{d}}\upsilon } \right] + \tilde{y}_{p} \left( t \right)\sin \left[ {2{\uppi }\int_{0}^{t} {\tilde{f}_{p} \left( \upsilon \right)} {\text{d}}\upsilon } \right].$$
(6)

Moreover, a high-resolution TFR of the multi-component signal can be constructed as

$${\text{ACMD}}\left( {t,f} \right) = \sum\limits_{p = 1}^{P} {\tilde{A}_{p} \left( t \right)\delta \left[ {f - \tilde{f}_{p} \left( t \right)} \right]},$$
(7)

where \(\tilde{A}_{p} \left( t \right) = \sqrt {\tilde{x}_{p}^{2} \left( t \right) + \tilde{y}_{p}^{2} \left( t \right)}\) is estimated IA, and \(\delta \left( \cdot \right)\) denotes Dirac delta function.

3 Wheel polygonization detection using ACMD

The proposed wheel polygonization detection method based on ACMD is presented in this section. The method mainly includes three steps: (1) wheel rotating frequency estimation; (2) signal resampling and order analysis; (3) polygonal wear amplitude estimation.

3.1 Wheel rotating frequency estimation

As shown in Eq. (2), accurately estimating the wheel rotating frequency \(f_{{\text{w}}} \left( t \right)\) is crucial for detection of harmonic orders of the polygonal wear. Most existing methods employ a speed coder to measure the vehicle speed for wheel rotating frequency estimation [23]. However, due to the limited space and unpredictable safety concerns, it may not be allowed to install a speed measuring device in a railway vehicle. To address this issue, a novel rotating frequency estimation method based on vibration signal analysis is proposed in this paper. Note that it is difficult to directly detect wheel rotating frequency from vibration signal of a vehicle with a healthy wheel shaft. Therefore, it is suggested to extract other dominant frequencies and then obtain the rotating frequency according to transmission relation. It is known that gear transmission systems are widely used in a locomotive. Due to the continuous meshing operation, there exists a distinct meshing frequency component in vibration response of the gear transmission system of the locomotive. Therefore, we propose to estimate the gear meshing frequency \(f_{{\text{m}}} \left( t \right)\) at first and then calculate the wheel rotating frequency as

$$f_{{\text{w}}} \left( t \right){ = }\frac{{f_{{\text{m}}} \left( t \right)}}{{Z_{{{\text{shaft}}}} }},$$
(8)

where \(Z_{{{\text{shaft}}}}\) denotes the number of the teeth of the gear fixedly connected with the wheel shaft.

To accurately estimate the gear meshing frequency, the ACMD is employed to analyze vibration acceleration signals of motor which is directly connected with gearbox and is less influenced by wheel-rail excitations [39, 40]. Since IF initialization plays an important role in analysis results of ACMD, it is necessary to obtain a good estimate of an initial IF. To this end, a TF ridge detection technique is utilized to deal with the IF initialization issue. Firstly, a TFR of the vibration signal is obtained by a proper TF transform. For simplicity, STFT of the signal \(g\left( t \right)\) is calculated as

$${\text{STFT}}\left( {t,f} \right) = \int_{ - \infty }^{ + \infty } {g\left( \upsilon \right)h^{ * } \left( {\upsilon - t} \right){\text{e}}^{{ - {\text{j}}2{\uppi }f\upsilon }} {\text{d}}\upsilon },$$
(9)

where \(\text{j}^{2} = - 1\), \(h\left( t \right)\) is a nonnegative, symmetric and real window, and * stands for complex conjugate. Ridge curve of the TFR is often deemed as a good estimate of the IF and can be detected by solving [41, 42]

$$\tilde{f}_{{\text{m}}}^{{\text{r}}} \left( t \right) = \mathop {\arg \max }\limits_{\varOmega \left( t \right)} \left\{ {\sum\limits_{i = 0}^{N - 1} {\left| {{\text{STFT}}\left( {t_{i} ,\varOmega \left( {t_{i} } \right)} \right)} \right|^{2} } - \lambda \sum\limits_{i = 1}^{N - 1} {\left| {\varOmega \left( {t_{i} } \right) - \varOmega \left( {t_{i - 1} } \right)} \right|^{2} } } \right\},$$
(10)

where \(t = t_{0} , t_{1},\, \ldots ,\,t_{N - 1}\) denote sampling time instants, N is the number of samples, \(\varOmega \left( t \right)\) stands for a set of all the TF paths from \(t = t_{0}\) to \(t = t_{N - 1}\), and \(\lambda\) is a penalty factor. Eq. (10) implies that a desired ridge curve should pass through TF points with large magnitudes and there is no significant jump phenomenon between adjacent points (i.e., the curve is smooth enough).

The accuracy of IF estimation by TF ridge detection is dominated by resolution of the TFR. Due to the limitation of traditional TF methods in analyzing strongly modulated signals, it is difficult to get a high-accuracy estimation of the IF by ridge detection. Therefore, the ACMD is further applied to refine IF estimation. Firstly, with the coarsely estimated gear meshing frequency by ridge detection, i.e., \(\tilde{f}_{{\text{m}}}^{{\text{r}}} \left( t \right)\), the demodulated signals can be obtained as (see Eq. (5); since ACMD is only employed to estimate the gear meshing frequency at this stage, the number of signal components is simply set to 1, i.e., \(P = 1\))

$$\left\{ {\tilde{x}_{{\text{m}}} \left( t \right),\tilde{y}_{{\text{m}}} \left( t \right)} \right\} = \mathop {\arg \min }\limits_{{x_{{\text{m}}} \left( t \right),y_{{\text{m}}} \left( t \right)}} \left\{ {{\mathcal{L}}_{\tau } \left( {x_{{\text{m}}} \left( t \right),y_{{\text{m}}} \left( t \right),\tilde{f}_{{\text{m}}}^{{\text{r}}} \left( t \right)} \right)} \right\},$$
(11)

where \(x_{{\text{m}}} \left( t \right)\) and \(y_{{\text{m}}} \left( t \right)\) denote the demodulated signals corresponding to gear meshing frequency. Then, according to Eq. (4), phase information of the demodulated signals can be used for calculating IF deviation as

$$\Delta f_{{\text{m}}} \left( t \right) = - \frac{1}{{2{\uppi }}}\frac{{\text{d}}}{{{\text{d}}t}}\left\{ {\arctan \left[ {\frac{{\tilde{y}_{{\text{m}}} \left( t \right)}}{{\tilde{x}_{{\text{m}}} \left( t \right)}}} \right]} \right\} = \frac{1}{{2{\uppi }}}\left\{ {\frac{{\tilde{y}_{{\text{m}}} \left( t \right) \cdot \tilde{x}_{{\text{m}}}^{\prime } \left( t \right) - \tilde{x}_{{\text{m}}} \left( t \right) \cdot \tilde{y}_{{\text{m}}}^{\prime } \left( t \right)}}{{\tilde{x}_{{\text{m}}}^{2} \left( t \right) + \tilde{y}_{{\text{m}}}^{2} \left( t \right)}}} \right\},$$
(12)

where \(\cdot^{\prime}\) denotes derivative. Finally, estimation of the gear meshing frequency can be improved by error compensation as

$$\tilde{f}_{{\text{m}}} \left( t \right){ = }\tilde{f}_{{\text{m}}}^{{\text{r}}} \left( t \right){ + }\ell \left\{ {\Delta f_{{\text{m}}} \left( t \right)} \right\},$$
(13)

where \(\ell \left\{ \cdot \right\}\) represents a low-pass filtering operator [35] for reducing noise interference.

3.2 Signal resampling and order analysis

With the method proposed in Sect. 3.1, the gear meshing frequency \(f_{{\text{m}}} \left( t \right)\) can be accurately estimated and a wheel rotating frequency \(f_{{\text{w}}} \left( t \right)\) can be calculated according to Eq. (8). The obtained \(f_{{\text{w}}} \left( t \right)\) is used as a reference frequency for signal resampling and order analysis. At this stage, a vibration acceleration signal of axle box is analyzed since the signal is sensitive to the wheel polygonization. Firstly, an angular variable is defined as \(s = \phi_{{\text{w}}} \left( t \right) = \int_{0}^{t} {f_{{\text{w}}} \left( \upsilon \right){\text{d}}\upsilon }\) and therefore Eq. (2) can be rewritten as

$$g\left( t \right) = \sum\limits_{p = 1}^{P} {A_{p} \left( t \right)\cos \left[ {2{\uppi }o_{p} \phi_{{\text{w}}} \left( t \right) + \varphi_{p} } \right]}.$$
(14)

Then, if the time variable is expressed with the angular variable as \(t = \phi_{{\text{w}}}^{ - 1} \left( s \right)\) where \(\phi_{{\text{w}}}^{ - 1}\) denotes inverse function, Eq. (14) can be mapped into angular domain as

$$g_{\phi } \left( s \right) = \sum\limits_{p = 1}^{P} {A_{p}^{\phi } \left( s \right)\cos \left[ {2{\uppi }o_{p} s + \varphi_{p} } \right]},$$
(15)

where \(g_{\phi } \left( s \right){ = }g\left( {\phi_{{\text{w}}}^{ - 1} \left( s \right)} \right)\) and \(A_{p}^{\phi } \left( s \right){ = }A_{p} \left( {\phi_{{\text{w}}}^{ - 1} \left( s \right)} \right)\). It can be seen that the angular-domain signal in Eq. (15) is free of FM effect and exhibits constant frequencies at \(o_{p}\) for \(p = 1, 2, \cdots ,P\). Therefore, it is easy to identify harmonic orders (i.e., \(o_{p}\)) by applying Fourier transform to the signal. In practice, the nonlinear mapping operation (i.e., Eq. (14) to Eq. (15)) is achieved by signal resampling or data interpolation [43, 44]. Namely, the time-domain signal \(g\left( t \right)\) is interpolated from discrete time points \(\left\{ {t_{n} } \right\}_{n = 0, 1, \cdots ,N - 1}\) to discrete angles \(\left\{ {s_{l} } \right\}_{l = 0, 1, \cdots ,L - 1}\) as

$$\begin{gathered} g_{\phi } \left( {s_{l} } \right) = {\text{Interpolate}}\left( {\phi_{{\text{w}}} \left( {t_{n} } \right),g\left( {t_{n} } \right),s_{l} } \right), \hfill \\ \, n = 0, 1, \cdots , \, N - 1, \, l = 0, 1, \cdots , \, L - 1, \hfill \\ \end{gathered}$$
(16)

where \(s_{l} = {{l\phi_{{\text{w}}} \left( {t_{N - 1} } \right)} \mathord{\left/ {\vphantom {{l\phi_{{\text{w}}} \left( {t_{N - 1} } \right)} {\left( {L - 1} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {L - 1} \right)}}\) for \(l = 0, 1,\, \ldots ,\,L - 1\) denote uniformly discretized samples of angular variable \(s\), and \(L \ge N\) represents the number of angular samples.

3.3 Polygonal wear amplitude estimation

Apart from harmonic order, harmonic amplitude of the polygonal wear is also a key factor which directly determines whether a wheel should be repaired. In reality, the polygonal wear amplitude is detected through static measurement of wheel roughness levels using specialized instruments. To the best of our knowledge, on-board polygonal wear amplitude detection methods based on vehicle vibration signals have hardly been reported yet. Note that the wheel polygonal wear along circumferential direction can be regarded as a type of irregularity whose dynamic effect on a vehicle is similar to that of rail corrugation. Therefore, we propose to estimate the polygonal wear amplitude by combining the ACMD with track irregularity detection methods. In this paper, a well-known irregularity measurement method called inertial algorithm is employed due to its easy implementation and wide applicability. As illustrated in Fig. 1, the basic principle of inertial algorithm is that the quadratic integral of acceleration gives displacement [45, 46]. Therefore, the track irregularity can be measured by processing the acceleration signal of axle box as

$$d\left( t \right) = \iint {a\left( t \right)}{\text{d}}t{\text{d}}t,$$
(17)

where \(a\left( t \right)\) denotes acceleration signal. Note that the integral operation in Eq. (17) may introduce an undesired signal trend which should be removed in a post-processing step.

Fig. 1
figure 1

Illustration of the principle of inertial algorithm

Full size image

It is worth noting that apart from harmonics induced by polygonal wear, the acceleration signal of axle box contains other interference components such as those caused by random track irregularities. To correctly estimate polygonal wear amplitudes, the ACMD is employed to extract the harmonics to be interested before the inertial algorithm is performed. Procedures of the proposed amplitude estimation method are listed as follows:

  1. (1)

    Construct initial IFs for the harmonics of the acceleration signal as \(\tilde{o}_{p} \tilde{f}_{{\text{w}}} \left( t \right)\), \(p = 1, 2, \cdots ,P\), according to an improved ACMD in [36], where \(\tilde{f}_{{\text{w}}} \left( t \right)\) and \(\tilde{o}_{p}\) denote the estimated rotating frequency and harmonic orders in Sects. 3.1 and 3.2, respectively.

  2. (2)

    Reconstruct harmonics by ACMD according to the obtained initial IFs and then summate these harmonics as

    $$\tilde{g}\left( t \right) = \sum\limits_{p = 1}^{P} {\tilde{g}_{p} \left( t \right)}.$$
    (18)
  3. (3)

    Calculate quadratic integral of the reconstructed signal in Eq. (18) to obtain a displacement signal \(\tilde{d}\left( t \right)\). Next, least squares spline fitting is applied to \(\tilde{d}\left( t \right)\) to obtain a signal trend \(\tilde{r}\left( t \right)\) which is then subtracted to get the irregularity caused by polygonal wear as

    $$\tilde{Z}\left( t \right) = \tilde{d}\left( t \right) - \tilde{r}\left( t \right).$$
    (19)
  4. (4)

    Perform the order analysis approach in Sect. 3.2 to the irregularity \(\tilde{Z}\left( t \right)\) to estimate amplitudes of the harmonics of the wheel polygonal wear.

Note that it is not possible to extract all the harmonics for polygonal wear identification. Therefore, only dominant harmonics (e.g., higher-order harmonics) in order spectrum (in step 1) are reconstructed since they cause stronger vibration to a locomotive. The whole flow chart of the proposed wheel polygonal wear detection method is presented in Fig. 2. It can be seen that, by synthetically analyzing vibration acceleration signals of motor and axle box, the proposed method can estimate the wheel rotating frequency and quantificationally detect the harmonic orders and their amplitudes of wheel polygonal wear without using speed coder.

Fig. 2
figure 2

Flow chart of the proposed wheel polygonization detection method

Full size image

4 Dynamics simulation

4.1 Dynamics model

In this section, the proposed polygonal wear detection method is validated by dynamics simulation. To this end, a locomotive-track coupled dynamics model considering gear transmissions [47,48,49] is introduced as shown in Fig. 3. The model fully takes into account rail-wheel interaction and thus can well simulate practical locomotive vibration responses.

Fig. 3
figure 3

Locomotive-track coupled dynamics model considering gear transmissions: a model and b gear transmission systems in one bogie

Full size image

The dynamics model contains three sub-models including locomotive, track and gear transmission sub-models. In the locomotive sub-model, a car body is linked with two bogie frames via secondary suspension (\(K_{{{\text{sz}}}}\)–\(C_{{{\text{sz}}}}\)) while each bogie frame is supported by two wheelsets through primary suspension (\(K_{{{\text{pz}}}}\)–\(C_{{{\text{pz}}}}\)). In each bogie, there are two gearboxes which are rigidly connected with two traction motors, respectively. In the dynamics model, most of the locomotive components are regarded as rigid bodies represented by mass (\(M_{q}\)) and rotational inertia (\(J_{q}\)), and each component has three degrees of freedom, i.e., vertical (\(Z_{q}\)), rotational (\(\beta_{q}\)) and longitudinal (\(X_{q}\)) motions, where \(q = {\text{c}},{\text{ m, w}}\) denote the car body, motor and wheelset, respectively. The track subsystem consists of rail, rail pads, sleepers, ballasts and subgrade. The stiffness and damping of the rail pads, ballasts and subgrade are represented as \(K_{{{\text{p}}i}}\)–\(C_{{{\text{p}}i}}\), \(K_{{{\text{b}}i}}\)–\(C_{{{\text{b}}i}}\), and \(K_{{{\text{f}}i}}\)–\(C_{{{\text{f}}i}}\), respectively, where i denotes sleeper number. As for the gear transmission subsystem, the motor with gearbox is suspended through the bearings (\(K_{{{\text{br}}}}\)–\(C_{{{\text{br}}}}\)) mounted on the wheel axle and the hanger rod (\(K_{{{\text{ms}}}}\)–\(C_{{{\text{ms}}}}\)) linked with bogie frame. The pinion is connected with rotor of the motor and the gear is mounted on the wheel axle. The power of the traction motor is transmitted from the pinion to gear via their teeth engagement described by the stiffness and damping elements (\(K_{{\text{m}}}\)–\(C_{{\text{m}}}\)) along the line of action (LOA). More details about the model and parameter setting can be found in [47,48,49].

4.2 Response analysis

The dynamics model in Fig. 3 is employed to calculate locomotive vibration responses caused by wheel polygonal wear. Firstly, the irregularity induced by the polygonal wear with P harmonics is used as an excitation of the model as

$$Z_{0} \left( t \right) = \sum\limits_{p = 1}^{P} {\varOmega_{p} \cos \left[ {\frac{{o_{p} u\left( t \right)}}{{R_{{\text{w}}} }} + \psi_{p} } \right]},$$
(20)

where \(\varOmega_{p}\), \(o_{p}\) and \(\psi_{p}\) stand for amplitude, harmonic order and initial phase of the polygonal wear, respectively, \(u\left( t \right)\) denotes running distance function of locomotive (i.e., the integral of locomotive speed), and \(R_{{\text{w}}}\) is wheel radius (the radius is set to 0.625 m in the simulation). In this section, the polygonal wear with three harmonics is considered and the parameters are provided in Table 1. Figure 4 illustrates the polygonal wear and its order analysis result. To simulate the random interference of track irregularity, the American fifth-grade track spectrum is employed to generate random excitation of the dynamics model. It is assumed that the locomotive is running under a traction condition and the speed is gradually increased to 80 km/h.

Table 1 Simulated wheel polygonal wear with three harmonics
Full size table
Fig. 4
figure 4

Wheel polygonal wear in the simulation case: a polygonal wear described in the polar coordinate and b order analysis of the polygonal wear

Full size image

The simulated vibration acceleration signals of axle box and motor based on the dynamics model are provided in Fig. 5. The gear meshing frequency and frequency contents induced by wheel polygonal wear can be observed in the TFRs by STFT. However, due to poor resolution of the TFR, details of the wheel polygonization are not available. Next, the proposed fault detection method is applied. Firstly, as shown in Fig. 6, the gear meshing frequency is estimated and then the wheel rotating frequency can be calculated according to Eq. (8), where the number of the gear teeth is set to \(Z_{{{\text{shaft}}}} { = }120\) for the locomotive. The mean-square errors of the estimated gear meshing frequencies by ridge detection and ACMD are 6.43 and − 5.73 dB, respectively. It can be seen that the IF estimation results by ridge detection exhibit staircase effects due to limited resolution of the TFR (see Fig. 6a, b). After frequency refinement by ACMD, the IF estimation accuracy can be significantly improved.

Fig. 5
figure 5

Simulated vibration signals of the locomotive with wheel polygonization: a and b show waveform and STFT of the acceleration signal of axle box, respectively; c and d illustrate the signal of motor

Full size image
Fig. 6
figure 6

IF estimation results from the simulated signal in Fig. 5(d): a and b show the estimated gear meshing frequency and wheel rotating frequency by ridge detection, respectively; c and d show the refined results by ACMD

Full size image

With the estimated wheel rotating frequency, the signal resampling technique is applied to acceleration signal of axle box and then the order spectrum is obtained and compared with original spectrum as shown in Fig. 7. Due to the time-varying frequency contents, Fourier spectrum cannot reveal effective fault information of the polygonal wear. On the contrary, since the non-stationary effect of the signal is fully removed by resampling technique, harmonic orders of the polygonal wear can be clearly identified in order spectrum (see Fig. 7b). In addition, the order spectrum based on acceleration signal shows much better resolution than that in Fig. 4b since the duration of the acceleration signal is long enough.

Fig. 7
figure 7

Comparison of spectrums of the simulated vibration signals of axle box: a original Fourier spectrum; b order spectrum after signal resampling with the estimated wheel rotating frequency

Full size image

Then, polygonal wear amplitude estimation is carried out. Firstly, based on the order analysis results, the harmonics induced by polygonal wear can be extracted by ACMD as shown in Fig. 8. An improved TFR of these harmonics by ACMD (see Eq. (7)) is also provided and compared with that by STFT as shown in Fig. 9. It can be observed that STFT cannot resolve these harmonics due to the poor resolution and interferences from track irregularity. The TFR obtained by ACMD can clearly uncover time-varying features of the harmonics. Next, inertial algorithm is applied to the extracted harmonics and the analysis results are shown in Fig. 10. It is indicated that the integral operation will lead to a complex signal trend with a large magnitude and therefore signal detrending is necessary for detection of the irregularity caused by polygonal wear. After detrending, radial deviations (or irregularities) caused by polygonal wear of different wheel rotation cycles are clearly recovered as illustrated in Fig. 10b. Finally, amplitudes of harmonics of the polygonal wear are estimated by performing order analysis to the recovered radial deviation of a cycle (suggested to use average data of different rotation cycles at lower-speed stage to reduce the influence of track irregularities) as shown in Fig. 10c. One can find that the estimated harmonic amplitudes are very close to the true ones showing the effectiveness of the proposed method for quantitative detection of locomotive wheel polygonal wear.

Fig. 8
figure 8

Harmonic extraction by ACMD for the simulated vibration signals of axle box: a extracted harmonics; b sum of the extracted harmonics

Full size image
Fig. 9
figure 9

Comparison of TFRs of the simulated vibration signals of axle box: a STFT and b ACMD

Full size image
Fig. 10
figure 10

Polygonal wear amplitude estimation by inertial algorithm in the simulation case: a obtained wheel radial deviation by quadratic integral of the reconstructed harmonics; b radial deviation after detrending; c order analysis of the reconstructed polygonal wear

Full size image

5 Practical application

In this section, the proposed fault detection method is applied to vibration signal analysis of an in-service locomotive to further demonstrate the effectiveness of the method.

5.1 Field test

The locomotive operated on the Chinese heavy-haul railway from Zhongwei to Yuci. Vibration alarm phenomenon occurred several times for on-board monitoring system of the locomotive. To investigate root cause of the alarm phenomenon, a research group from Train and Track Research Institute (in State Key Laboratory of Traction Power) carried out a thorough field test for the locomotive [15]. Firstly, wheel roughness levels around circumference were measured by Müller-BBM instrument to assess the wheel polygonal wear as shown in Fig. 11. The measuring accuracy of the instrument is 0.1 μm with a sampling precision of 1 mm. The roughness measurement was achieved by manually twirling wheels. In addition, vibration accelerations of key components such as axle box and motor were also measured during operation process of the locomotive as shown in Fig. 12. To determine whether the root cause of the vibration alarm is the wheel wear, two round tests before and after wheel re-profiling were conducted, respectively. The wheel re-profiling was achieved by wheel lathe to reduce radial deviation caused by polygonal wear as shown in Fig. 13.

Fig. 11
figure 11

Field measurement of wheel polygonal wear of the locomotive [15]: a wheel polygonal wear; b polygonal wear measurement by Müller-BBM instrument

Full size image
Fig. 12
figure 12

Field measurement of vibration of the locomotive: a the locomotive; b schematic diagram of measurement positions; c Measurement positions of axle box; d measurement positions of motor

Full size image
Fig. 13
figure 13

Wheel re-profiling with wheel lathe: a schematic diagram and b wheel lathe [15]

Full size image

5.2 Detection results

Firstly, the measured data of the locomotive before wheel re-profiling is analyzed. The measured wheel roughness data and its analysis result are provided in Fig. 14. It clearly shows that the wheel suffers from high-order polygonal wear (i.e., the 17th–24th orders). Then, the proposed method is applied to locomotive vibration data to test whether the polygonal wear information can be detected. In this test, the locomotive operates under a traction condition. Figure 15 illustrates the measured vibration signals and their TF analysis results. It can be seen that the dominant frequency contents of the vibration signal of axle box are induced by polygonal wear while those of the motor vibration signal are mainly related to gear meshing operation. Therefore, the proposed method is applied to the motor vibration signal to estimate the gear meshing frequency and wheel rotating frequency as shown in Fig. 16, where the number of gear teeth \(Z_{{{\text{shaft}}}}\) is 120. After the rotating frequency is obtained, signal resampling and order analysis are performed to the vibration signal of axle box as illustrated in Fig. 17. Compared with the raw Fourier spectrum (see Fig. 17a), the developed order spectrum (see Fig. 17b) correctly reveals all the harmonic orders of the polygonal wear and it even shows much better resolution than the analysis results based on the measured wheel roughness in Fig. 14b. TFRs of the polygonal-wear-related harmonics are shown in Fig. 18. It can be seen that ACMD successfully resolves and extracts all the harmonics and thus generates a high-resolution TFR (see Fig. 18b; due to the limited space, instead of waveforms, only TFRs are provided). These extracted harmonics can be further used for wear amplitude estimation.

Fig. 14
figure 14

Measured wheel polygonal wear of the locomotive: a polygonal wear described in the polar coordinate; b order analysis of the polygonal wear

Full size image
Fig. 15
figure 15

Measured vibration signals of the locomotive with wheel polygonization: a and b show waveform and STFT of the acceleration signal of axle box, respectively; c and d illustrate the signal of motor

Full size image
Fig. 16
figure 16

Estimated gear meshing frequency and wheel rotating frequency by ACMD from the signal in Fig. 15d: a gear meshing frequency and b wheel rotating frequency

Full size image
Fig. 17
figure 17

Comparison of spectrums of the measured vibration signals of axle box: a original Fourier spectrum; b order spectrum after signal resampling with the estimated wheel rotating frequency

Full size image
Fig. 18
figure 18

Comparison of TFRs of the measured vibration signals of axle box: a STFT and b ACMD

Full size image

For comparison, wheel roughness and locomotive vibration data after wheel re-profiling are analyzed. Figure 19 shows the measured roughness data. It can be seen that the wheel lathe does not fully eliminate wheel polygonization although the roughness magnitudes and the number of harmonics are significantly reduced (compared with that in Fig. 14). In this case, the dominant harmonic orders of the polygonal wear are 17 and 24. To show the effectiveness of the method in different working conditions, vibration signals of the locomotive during braking process are analyzed as illustrated in Fig. 20. It can be found that magnitudes of the vibration acceleration signals and energy of the polygonal-wear-induced components are decreased evidently. It implies that the vibration alarm phenomenon of the locomotive is mainly attributed to the wheel polygonal wear. Then, the proposed method is applied to demonstrate its usefulness in detecting such a slight fault. IF estimation and order analysis results based on the vibration signals are given in Figs. 21 and 22, respectively. It shows that the harmonic orders of the polygonal wear are clearly identified (see Fig. 22b). Comparison of the TFRs obtained by STFT and ACMD is illustrated in Fig. 23 indicating higher resolution of ACMD in separating closely spaced harmonics induced by wheel polygonal wear.

Fig. 19
figure 19

Measured wheel polygonal wear of the locomotive after wheel re-profiling: a polygonal wear described in the polar coordinate; b order analysis of the polygonal wear

Full size image
Fig. 20
figure 20

Measured vibration signals of the locomotive after wheel re-profiling: a and b show waveform and STFT of the acceleration signal of axle box, respectively; c and d illustrate the signal of motor

Full size image
Fig. 21
figure 21

Estimated gear meshing frequency and wheel rotating frequency by ACMD from the signal in Fig. 20d: a gear meshing frequency and b wheel rotating frequency

Full size image
Fig. 22
figure 22

Comparison of spectrums of the measured vibration signals of axle box after wheel re-profiling: a original Fourier spectrum; b order spectrum after signal resampling with the estimated wheel rotating frequency

Full size image
Fig. 23
figure 23

Comparison of TFRs of the measured vibration signals of axle box after wheel re-profiling: a STFT and b ACMD

Full size image

Finally, inertial algorithm is applied to the harmonics extracted by ACMD for polygonal wear amplitude estimation. The amplitude estimation results before and after wheel re-profiling are compared in Fig. 24. In both cases, the estimated polygonal wear amplitudes based on acceleration signals are close to those measured by instrument as shown in Fig. 24. The results indicate that the proposed method is effective for detection of wheel polygonal wear of different levels.

Fig. 24
figure 24

Estimated polygonal wear amplitudes by inertial algorithm for the experimental vibration signals: a before wheel re-profiling and b after wheel re-profiling

Full size image

6 Conclusions

In this paper, a novel quantitative detection method for wheel polygonal wear under non-stationary conditions using locomotive vibration signals has been developed. Firstly, TF ridge detection has been combined with ACMD to accurately estimate the time-varying wheel rotating frequency from motor vibration signal without using speed coder. Next, the order analysis technique based on signal resampling has been applied to vibration signal of axle box for harmonic order detection. Finally, polygonal wear amplitudes have been estimated by integrating the ACMD with inertial algorithm. Dynamics simulation based on a locomotive-track coupled model considering gear transmissions has been carried out showing that the proposed method can achieve high-accuracy detection of both harmonic orders and amplitudes of wheel polygonal wear. Field tests of a locomotive with wheel polygonization before and after wheel re-profiling have also been carried out indicating that the proposed method can effectively detect higher-order polygonal wear of different levels.

Our future work will focus on the integration of ACMD with advanced machine learning (ML) methods [50] for automatic detection of wheel polygonal wear of different kinds of railway vehicles. For example, in initialization stage, ML methods can be used to automatically estimate initial IFs from TF images while in decision stage, ML methods can be applied to the signal modes extracted by ACMD for automatic identification of harmonic orders and amplitudes of polygonal wear.