Wayside Bearing Fault Diagnosis Based on a Data-Driven Doppler Effect Eliminator and Transient Model Analysis

Abstract

A fault diagnosis strategy based on the wayside acoustic monitoring technique is investigated for locomotive bearing fault diagnosis. Inspired by the transient modeling analysis method based on correlation filtering analysis, a so-called Parametric-Mother-Doppler-Wavelet (PMDW) is constructed with six parameters, including a center characteristic frequency and five kinematic model parameters. A Doppler effect eliminator containing a PMDW generator, a correlation filtering analysis module, and a signal resampler is invented to eliminate the Doppler effect embedded in the acoustic signal of the recorded bearing. Through the Doppler effect eliminator, the five kinematic model parameters can be identified based on the signal itself. Then, the signal resampler is applied to eliminate the Doppler effect using the identified parameters. With the ability to detect early bearing faults, the transient model analysis method is employed to detect localized bearing faults after the embedded Doppler effect is eliminated. The effectiveness of the proposed fault diagnosis strategy is verified via simulation studies and applications to diagnose locomotive roller bearing defects.

1. Introduction

Bearing defects are the dominant type of fault in railway vehicles, which leads to serious accidents and significant costs for the rail transport industry [1]. Approximately 50 bearing-related derailments occur in the United States each year [2]. Thus, accurately and automatically detecting and diagnosing the existence and severity of these faults in the bearings are significant [3]. Wayside acoustic defective bearing detection techniques are based on the assumption that diagnostically relevant information is stored in the acoustic signal generated by the bearings of a passing vehicle. With the signal processing techniques, fault characteristic information can be extracted from the acoustic signals. A successful example is the wayside acoustic defective bearing detector (ADBD) system developed in the 1980s [4]. Compared with other systems, the ADBD system entails lower costs and can detect bearing defects earlier in the failure process or before overheating occurs, thereby allowing scheduled bearing maintenance [5].

However, the effectiveness of the wayside acoustic monitoring-based technique is decreased when vehicles pass by at high speeds. One of the problems caused by the high relative movement is the Doppler effect, as it can lead to obvious frequency shifts, frequency band expansion, and amplitude modulation for the recorded acoustic signal, which reduces the diagnostic performance [6].

Aside from the wayside acoustic monitoring system, the Doppler effect also widely exists in the signal from a moving acoustic source. For example, in the areas of underwater acoustic communication and acoustical holography for moving vehicles, the Doppler effect contained in the recorded acoustic signal is also a barrier that could significantly diminish the effectiveness of signal processing. Stojanovic et al. [7] and Johnson et al. [8] proposed a Doppler compensation method jointly based on phase synchronization and channel equalization. Yang and Wang [9] established the time-space relation among the measurement field, radiating field, and acoustic holography field; they also proposed a method based on the nonlinear mapping function between the sound source and the measured signal, in which the Doppler effect was eliminated. However, Doppler reduction methods for the wayside acoustic defective bearing detection method are rarely reported in current papers.

Dybała [10] proposed a disturbance-oriented dynamic signal resampling method based on the Hilbert Transform to eliminate the Doppler effect for wayside monitoring systems. A time-domain interpolation resampling (TIR) method is proposed by Liu [11] to remove the Doppler effect embedded in the acoustic signal. Shen [12] constructed a Doppler transient model that combines the Doppler transient model and parameter identification based on the Laplace wavelet and a spectrum correlation assessment to detect the locomotive bearing fault. However, in the above two methods, the source is assumed to move at a uniform velocity, and the sound speed is set at a constant value. In practice, the moving speed of the railway vehicle is unstable, and parameters such as the temperature and atmospheric pressure of the medium in which the wave propagates would definitely affect the value of the sound speed. These parameters could also have a significant effect on the results of the methods. What's more, all of the geometric parameters of the kinematic model are assumed to be known in advance in the proposed method.

In this paper, a fault diagnosis strategy is invented for locomotive bearing fault diagnosis based on the wayside acoustical monitoring technique. Through the proposed strategy, the Doppler effect embedded in the recorded bearing acoustic signal can be eliminated by a Doppler effect eliminator. In the Doppler effect eliminator, a so-called Parametric-Mother-Doppler-Wave (PMDW) is constructed based on the kinematic model parameters. The kinematic model parameters, including the moving speed and sound speed of the railway vehicle, as well as geometric parameters of the model are then identified via correlation filtering analysis. All of the parameters can be identified based on the signal itself. A time domain signal resampler is invented and employed to eliminate the Doppler effect using the identified kinematic model parameters. After the embedded Doppler effect is eliminated, the transient model parameters of the Doppler-free signal are identified to detect the localized bearing faults.

The rest of this paper is organized as follows: Section 2 introduces the proposed locomotive bearing fault diagnosis strategy. The construction of the PMDW, kinematic model parameter identification based on correlation filtering analysis, Doppler effect elimination based on the resampling method, and fault feature extraction based on transient model analysis are all introduced in this section. A simulation case study is provided in Section 3. An experimental verification test using defective locomotive roller bearings with outer race defect and inner race defect is discussed in Section 4. Finally, Section 5 presents the concluding remarks.

2. Proposed Wayside Bearing Fault Diagnosis Strategy Based on a Data-Driven Doppler Effect Eliminator and Transient Model Analysis

In the wayside acoustic bearing monitoring system, microphones are fixed by the wayside to record the acoustic signals emitted by the bearings of a passing vehicle. The basic kinematic model (Figure 1) involving a single moving acoustic source and a microphone is considered in this study. The source is moving along a straight line, and the microphone is placed perpendicular to the trail of the moving source. When the vehicle passes by, the microphone receives the spherically attenuated signals with the Doppler effect from the bearing source.

Given the high relative speed between the railway vehicle and the microphone, the recorded signal is distorted by the Doppler effect, which causes the signal frequency to shift and the frequency band to expand. This condition is a barrier to further analysis, especially for the methods based on frequency domain analysis.

Two steps are implemented for the proposed strategy to detect the localized faults of the bearing of the fast moving railway vehicle. In the first step, the embedded Doppler effect is eliminated. The frequential structure disturbance is eliminated, and the amplitude is demodulated. In the next step, the transient model parameters of the Doppler-free signal are identified to detect the localized bearing faults. A flowchart of the proposed strategy that includes two signal processing modules is presented in Figure 2. In the following subsections, the two signal processing modules are discussed in detail.

2.1. Doppler Effect Eliminator

As shown in Figure 1, we assume that during the source point movement from A to B, a length of the acoustic Doppler-shifted signal, X_dop(t), is recorded by the microphone. The kinematic model (Figure 1) can be defined by the following kinematic model parameter set:

$${\gamma }_k=\left[c,r,X_0,V_0,\left[a\right]\right]$$

(1)

where c denotes the sound speed in the air, r is the distance from the microphone to the trail of the moving source point, X₀ stands for the distance between starting point A and point O′, and point O′ stands for the point where the source point is closest to the microphone during the pass-by movement. V₀ stands for the initial velocity. [a] is a vector containing m elements that stand for the 1 to m-order derivative of the source moving speed V, which is equal to the following:

$$\left[a\right]=\left[\frac{dV}{dt},\frac{d^{2}V}{d{t}^2},\mathrm{\dots },\frac{d^{m}V}{d{t}^m}\right]$$

(2)

When the different structures of the passing vehicles are considered, the parameter r is non-constant. The speed of sound propagation c is unstable because of the various atmospheric environments. The parameters X₀, V₀ and [a] are also difficult to measure accurately.

A data-driven Doppler effect eliminator is invented in this paper. Through this eliminator, all of the kinematic model parameters can be identified based on the recorded signal itself. A flowchart of the eliminator is shown in Figure 3. First, the kinematic model parameter set, γ_k, is identified by correlation filtering analysis between the input Doppler-shifted signal and the investigated PMDW generated by the PMDW generator. After the kinematic model parameters are identified, the Doppler-shifted signal is resampled by a signal resampler.

In the following subsection, the correlation filtering analysis method and the two signal processing modules (i.e., PMDW generator and signal resampler) are discussed in detail.

2.1.1. Parametric Mother-Doppler-Wave (PMDW) Generator

The construction of the PMDW is introduced in this section. As illustrated in Figure 1, we assume that the source point keeps emitting a pure harmonic sound wave during the movement from A to B. Let the harmonic sound wave be written as S_e = sin(2πf_ct_e), 0 < t_e < T. This wave can be sampled as a discrete amplitude sequence as Equation (3) at a sampling frequency of f_s. Thus:

$$S_e\left(n\right)=[\begin{array}{cccc}\text{sin}\left(2\pi f_c{t}_e\left(1\right)\right),& \text{sin}\left(2\pi f_c{t}_e\left(2\right)\right)& \mathrm{\dots }& \text{sin}\left(2\pi f_c{t}_e\left(N\right)\right)\end{array}]$$

(3)

that can be called “emit-amplitude-vector” and with the time sequence:

$$t_e\left(n\right)=[\begin{array}{cccc}0,& 1/\text{fs}& \mathrm{\dots }& \left(\text{N}-1\right)/\text{fs}\end{array}]$$

(4)

that can be called “emit-time-vector”. N equals the length of the signal. An amplitude weight of S_e(i) is emitted from the source point at a certain emit time t_e(i) and this amplitude weight arrives at the microphone at time t_r(i) where the t_r can be called “receive-time-vector”. Thus, we can determine that t_r is equal to the following:

$$t_r\left(n\right)=t_e\left(n\right)+R\left(n\right)/c$$

(5)

where c denotes sound speed, and R(n) stands for the distance between the source point to the microphone during the pass-by movement, which can be calculated by the following equation:

$$R\left(n\right)=\sqrt{{\left[X_0-X\left(n\right)\right]}^2+r^2}$$

(6)

where the X(n) denotes the displacement of the source point to the initial point A and it can be calculated as shown in Appendix A.1.

According to the Morse acoustic theory [13], the received amplitude weight can be described by the following equation:

$$S_r=\frac{r}{R{(1-Mcos\theta )}^2}\times S_e$$

(7)

where θ stands for the deflection angle of the source moving between the direction of source velocity and the microphone direction, and M = V/c stands for the mach number of the source point velocity. Thus, the “receive-amplitude-vector” can be calculated by the following equation:

$$S_r\left(n\right)=\frac{r}{R\left(n\right)\times {\left[1-\text{V}\left(n\right)/c\times cos\theta (n)\right]}^2}\times S_e\left(n\right)$$

(8)

where V(n) denotes the velocity of the source point equal to Appendix A.2. The cosθ(n) can be calculated by the following:

$$\mathit{cos}\theta (n)=X(n)/R(n)$$

(9)

The detailed procedure of the PMDW construction is presented in Figure 4. A corresponding example is shown in Figure 5a–d to illustrate the procedure more clearly. The additional parameter, f_c, stands for the center frequency, and the PMDW model parameter set is as follows:

$${\gamma }_{\Psi }=\left[f_c,c,r,X_0,V_0,\left[a\right]\right]$$

(10)

These parameters belong to the subsets S_fc, S_c, S_r, S_X0, S_V0, and S_a, which are shown as follows:

$$\begin{array}{l}{S}_{fc},S_c,S_r,S_{X0},S_{V0}\subset R^{+},\hfill \\ S_a\subset R\hfill \end{array}$$

(11)

Given an initialized model parameter set of γ_Ψ = [f_c, c, r, X₀, V₀, [a]], a PMDW can be constructed through the following steps:

(1)

Discrete sine wave generation. The center frequency of the sine wave is f_c. The discrete wave is an amplitude vector, S_e, (Equation (3)) with a time vector, t_e, (Equation (4)). An example is illustrated in Figure 5a.

(2)

Amplitude modulation. Through the “Doppler-amplitude-modulator” implemented by Equation (8), a receive-amplitude-vector, S_r, is obtained. The procedure is illustrated in Figure 5b.

(3)

Amplitude vector rearrangement and curve fitting. S_r obtained in Step 2 is rearranged by matching with the non-linear receive-time-vector, t_r, which is calculated through the receive-time-calculator implemented by Equation (5). A new fitting curve, χ^p, can then be obtained with t_r and S_r. The procedure is illustrated in Figure 5c.

(4)

Signal resampling. The amplitude vector of a PMDW is then obtained by resampling the curve, χ^p, obtained in Step 3 with the delayed-time-vector, t_d, which is equal to t_e + R₀/c. R₀ represents the distance between the starting point A and the microphone.

Two PMDWs are generated through the preceding steps. The PMDW shown in Figure 6a with its STFT spectrum in Figure 6b is generated with the following model parameter set: [1500, 340, 1, 4, 30, [0]]. The PMDW shown in Figure 6c with its STFT spectrum in Figure 6d is generated with the following parameter set: [1500, 340, 1, 4, 30, [100]].

2.1.2. Improved Correlation Filtering Analysis

Correlation filtering analysis aims to measure the strength and direction of the linear relationship between two signals by calculating the Pearson's correlation coefficient between them [14]. S.B. Wang et al. [15] and D. Wang, et al. [16] employed the correlation filtering analysis method to identify the transient model parameters that were then used as rotating machine fault detection features.

In the current paper, the kinematic model parameter set described by Equation (1) is identified by an improved correlation filtering analysis between the constructed PMDWs and the recorded Doppler-shifted signal. The PMDW is first transformed into a discrete analytic wavelet through the method investigated by Marple [17]. The Pearson's correlation coefficient between the recorded Doppler-shifted signal and the real part, as well as the imaginary part of the analytic wavelet, are then calculated. Finally, the square root of the two obtained coefficients is employed as a criterion of the inherent linear relationship between the two signals. Through this improved correlation filtering analysis method, the consideration of the initial phase of the two signals could be ignored, as the imaginary part of the analytic wavelet is a 90-degree phase shift signal of the PMDW.

Given an input Doppler-shifted signal, X_dop(n), and a parametric wavelet, Ψ(n), we can implement the improved correlation filtering analysis through the following procedure:

(1)

The discrete-time analytic signal of the parametric wavelet is obtained:

$$A*\left[{\Psi }_{\gamma }\left(n\right)\right]={\Psi }_{\gamma }\left(n\right)+H\left[{\Psi }_{\gamma }\left(n\right)\right]$$

(12)

(2)

The Pearson's correlation coefficients are calculated as follows:

$${\rho }_{X_{\mathit{\text{dop}}}\left(n\right),\text{R}\left(n\right)}=\frac{\sum _{i=1}^N\left(X_{\mathit{\text{dop}}}\left(i\right)-{\overline{X}}_{\mathit{\text{dop}}}\right)\left(\text{R}\left(i\right)-\overline{\text{R}}\right)}{\sqrt{\sum _{i=1}^N{\left(X_{\mathit{\text{dop}}}\left(i\right)-{\overline{X}}_{\mathit{\text{dop}}}\right)}^2}\sqrt{\sum _{i=1}^N{\left(\text{R}\left(i\right)-\overline{\text{R}}\right)}^2}}$$

(13)

$${\rho }_{X_{\mathit{\text{dop}}}\left(n\right),I\left(n\right)}=\frac{\sum _{i=1}^N\left(X_{\mathit{\text{dop}}}\left(i\right)-{\overline{X}}_{\mathit{\text{dop}}}\right)\left(I\left(i\right)-\overline{\text{I}}\right)}{\sqrt{\sum _{i=1}^N{\left(X_{\mathit{\text{dop}}}\left(i\right)-{\overline{X}}_{\mathit{\text{dop}}}\right)}^2}\sqrt{\sum _{i=1}^N{\left(I\left(i\right)-\overline{\text{I}}\right)}^2}}$$

(14)

where R(n) is the real part of A* [Ψ_γ(n)], and I(n) is the imaginary part of A* [Ψ_γ(n)], X̄_dop stands for the mean of X_dop, R̄ and Ī stand for the mean of R and I respectively.

(3)

The criterion of the inherent linear relationship between x(n) and Ψ(n) is as follows:

$${\rho }_{X_{\mathit{\text{dop}}}\left(n\right),{\Psi }_{\gamma }\left(\text{n}\right)}=\sqrt{{{\rho }^2}_{X_{\mathit{\text{dop}}}\left(n\right),\text{R}\left(\text{n}\right)}+{{\rho }^2}_{X_{\mathit{\text{dop}}}\left(n\right),\text{I}\left(\text{n}\right)}}$$

(15)

2.1.3. Signal Re-Sampler

After the kinematic model parameters are identified, the Doppler-shifted signal and the identified kinematic model parameters are inputted into the signal resampler, through which the embedded Doppler effect can be clearly eliminated.

If the identified kinematics model parameters are as follows:

$${{\gamma }^{\mathit{\text{opt}}}}_k=\left[c^{\mathit{\text{opt}}},r^{\mathit{\text{opt}}},{X_0}^{\mathit{\text{opt}}},{V_0}^{\mathit{\text{opt}}},\left[a^{\mathit{\text{opt}}}\right]\right]$$

(16)

and if the input Doppler-shifted signal is S_r^dop(n) with a length of N^dop points, then the variables similar to those mentioned in Section 2.1.1 can be calculated as follows:

Emit-time-vector:

$${t_e}^{\mathit{\text{dop}}}\left(n\right)=[\begin{array}{cccc}0,& 1/\text{fs}& \mathrm{\dots }& \left({\text{N}}^{\mathit{\text{dop}}}-1\right)/\text{fs}\end{array}]$$

(17)

where N^dop equals to the length of the input Doppler-shifted signal.

Receive-time-vector:

$${t_r}^{\mathit{\text{dop}}}\left(n\right)={t_e}^{\mathit{\text{dop}}}\left(n\right)+R^{\mathit{\text{dop}}}\left(n\right)/c$$

(18)

where R^dop (n) equals:

$$R^{\mathit{\text{dop}}}\left(n\right)=\sqrt{{\left[{X_0}^{\mathit{\text{opt}}}-X^{\mathit{\text{dop}}}\left(n\right)\right]}^2+r^{op{t}^2}}$$

(19)

where X^dop(n) equals to Appendix A.3.

According to Equation (7), the “emit-amplitude-vector” can be calculated by the following equation:

$${S_e}^{\mathit{\text{dop}}}\left(n\right)={S_r}^{\mathit{\text{dop}}}\left(n\right)\times {\left[\frac{r^{\mathit{\text{opt}}}}{R^{\mathit{\text{dop}}}\left(n\right)\times {\left[1-{\text{V}}^{\mathit{\text{dop}}}\left(n\right)/{\text{c}}^{\mathit{\text{opt}}}\times cos{\theta }^{\mathit{\text{dop}}}(\text{n})\right]}^2}\right]}^{-1}$$

(20)

where V^dop(n) equals to Appendix A.4.

and cosθ^dop (n) equals:

$$\mathit{cos}{\theta }^{\mathit{\text{dop}}}\left(n\right)=X^{\mathit{\text{dop}}}\left(n\right)/R^{\mathit{\text{dop}}}\left(n\right)$$

(21)

The detailed procedure of the signal resampler is shown in Figure 7. A corresponding example is shown in Figure 5e–h to illustrate the procedure more clearly.

Given a certain identified kinematic model parameter set of γ_k = [c^opt, r^opt, X₀^opt, V₀^opt, [a^opt]], the procedure of the resampler can be described as follows:

(1)

Curve fitting. The amplitude vector of the input Doppler-shifted signal is matched with the delayed-time-vector, t_d^dop, which equals [R₀^dop/c^opt, R₀^dop/c^opt + 1/f_s, …, R₀ ^dop/c^opt + (N^dop-1)/f_s]. A fitting curve, χ^d, can then be obtained by fitting the amplitude of the Doppler-shifted signal with t_d^dop. An example is illustrated in Figure 5e.

(2)

Signal resampling. The receive-amplitude-vector, S_r^dop(n), is then obtained by resampling the curve, χ^d, obtained in Step 1 with the receive-time-vector, t_r^dop, shown in Equation (18). This procedure is illustrated in Figure 5f.

(3)

Amplitude vector rearrangement. The receive-amplitude-vector, S_r(n), obtained in Step 2 is rearranged by matching with the linear emit-time-vector, t_e^dop, shown in Equation (17). This procedure is illustrated in Figure 5g.

(4)

Amplitude demodulation. The Doppler-free signal can be obtained by demodulating the rearranged amplitude vector in Step 2 through the “Doppler-amplitude-demodulator” implemented by Equation (20). This procedure is illustrated in Figure 5h.

2.2. Transient Model Analysis

The Doppler effect is eliminated through the aforementioned data-driven Doppler effect eliminator. Conventional fault feature extraction methods can be employed to analyze the Doppler-free signal, extract features, and then make a maintenance decision. In the past decades, numerous methods have been proposed to extract features for bearing signals, such as time-domain analysis [18], frequency-domain analysis [19], time-frequency-domain analysis [20–23], envelope spectrum [24,25], wavelet transform [26–29], empirical mode decomposition [30–32], and manifold learning [33].

A method of transient modeling by wavelet and parameter identification based on correlation filtering is first introduced and applied on bearing fault diagnosis by Wang et al. [15]. Whenever a bearing suffers a localized fault, the transients with a potential cyclic characteristic are generated by the rollers striking the localized fault. This phenomenon is an early bearing fault feature. The extraction of the transients is therefore beneficial to identify an early bearing fault [16]. The method is employed in this section to detect the localized defect of the locomotive bearing. A flowchart of this method is presented in Figure 8. It follows the steps of transient model construction, parameter identification through the correlation filtering analysis method, and bearing fault type identification through the recognized impact periods. Each step is discussed in detail in the following subsections.

2.2.1. Transient Model Construction Based on the Laplace Wavelet

The rolling element bearing typically consists of an inner race, an outer race, a number of rolling elements, and a cage. Once a localized fault is formed on the surface of the inner or outer race, a transient with an exponential decay is generated by the roller striking the localized fault.

The Laplace wavelet, a single-sided damped exponential function formulated as the impulse response of a single mode system, is highly similar to the waveform feature commonly encountered in bearing fault signal detection tasks.

The results reported in reference [15] show that the real part of the complex Laplace wavelet was the most sensitive wavelet to the transients generated by the localized bearing faults. The formula of the real part of the Laplace wavelet is as follows:

$${\psi }_{\text{Laplace}}(f,\zeta ,\tau ,t)=\{\begin{array}{l}{e}^{-\zeta /\sqrt{1-{\zeta }^2}2\pi f\times (t-\tau )}cos(2\pi f\times (t-\tau )),\tau \le t\le \tau +W\\ 0,\mathit{\text{else}}\end{array}$$

(22)

where W is the temporal range, f is the discrete frequency, ζ is the discrete damping coefficient, and τ is the delay time. These parameters belong to subsets F, Z, and T_d, which are shown as follows:

$$\begin{array}{c}F\subset R^{+}\\ Z\subset (R^{+}\cap [0,1))\\ T\subset R^{+}\end{array}$$

(23)

The speed variation has been removed during the procedure of Doppler effect elimination, however, the train of transients is not a strict periodic phenomenon when considering the “jitter” [34] during the operation. When considering this ‘jitter’, a different model should be used [35] and the traditional methods need to be modified [34]. The studies in this paper are based on the hypothesis that the train of transients is periodic after the speed variation has been removed. Then a periodic multi-transient model based on the Laplace wavelet is therefore constructed to simulate the waveform characteristics by introducing parameter T, as follows:

$$\chi \left(t\right)=\sum _m{\psi }_{\text{Laplace}}(f,\zeta ,\tau ,t-mT)$$

(24)

Figure 9 illustrates the single and periodic Laplace wavelet transient models, respectively.

2.2.2. Locomotive Bearing Fault Detection by Transient Model Parameters Identification

If the surface of the outer race of the bearing suffers a single defect based on the bearing geometries and rotation speed, f_r, the ball passing frequency over the outer race defect (BPFO) can be calculated as follows:

$$\mathit{\text{BPFO}}=\frac{1}{2}(1-\frac{d}{D_m}cos\alpha )f_{n}Z$$

(25)

where d and D_m represent the diameter of the rolling elements and the pitch diameter, respectively. α denotes the contact angle, f_n denotes the rotational frequency, and Z is the number of rolling elements.

Similarly, if a single defect occurs on the surface of the inner race of the bearing, the ball pass frequency over the inner race defect (BPFI) can be obtained by the following:

$$\mathit{\text{BPFI}}=\frac{1}{2}(1+\frac{d}{D_m}cos\alpha )f_{n}Z$$

(26)

Every time the rolling element passes through the defect, periodic impulses are created with time interval, Δt, as follows:

$$\mathrm{\Delta }t=\frac{1}{\mathit{\text{BPFO}}}$$

(27)

or:

$$\mathrm{\Delta }t=\frac{1}{\mathit{\text{BPFI}}}$$

(28)

After the Doppler effect embedded in the acoustic signal of the bearing is eliminated. The time interval can be identified through the improved correlation filtering analysis introduced in Section 2.1.2 between the periodic multi-transient model shown in Equation (24) and the Doppler-free signal of the bearing. The identified impact period in the transient model is the related bearing fault impact interval. The fault type can be determined by referring to the calculated theoretical fault-related impact intervals.

3. Simulation Case Study

In this section, a simulated Doppler-shifted bearing signal is analyzed to verify the effectiveness of the investigated diagnosis strategy. The source signal of the bearing without the Doppler effect can be described as follows:

$$X\left(t\right)=e^{-{\zeta }_0/\sqrt{1-{{\zeta }_0}^2}2\pi f_{0}mod(\text{t},1/40)}cos(2\pi f_{0}t)+n\left(t\right),t=0:1/50000:9999/50000$$

(29)

where the damping ratio, ζ₀, and frequency, f₀, are set at 0.05 and 1000 Hz, respectively. The impact interval embedded in the simulated signal is 0.025 s. The number of data points is 10,000, and the sampling frequency is 50 kHz.

A randomly distributed noise, n(t), is added to the simulated signal. The waveform of the simulated signal without noise is illustrated in Figure 10a. The polluted signal is shown in Figure 10b with its FFT spectrum in Figure 10c.

Through the PMDW generator introduced in Section 2.1.1, a sine wave can be embedded with the Doppler effect. If the sine wave is replaced with the aforementioned simulated signal, a Doppler-shifted bearing signal can be obtained with the same procedure and with the following kinematic parameters: c = 340 m/s, r = 1 m, X₀ = 4 m, V₀ = 30 m/s, [a] = [a₁] = 40 m/s². The wave form and FFT spectrum of the simulated Doppler-shifted bearing signal are illustrated in Figure 10d,e, respectively.

The proposed diagnosis strategy is then applied to the Doppler-shifted bearing signal. The Doppler-shifted bearing signal is first inputted to the Doppler effect eliminator introduced in Section 2.1 to eliminate the embedded Doppler effect. Through the eliminator, the kinematic parameters are first identified through the correlation filtering analysis introduced in Section 2.1.2 between the input Doppler-shifted signal and the PMDWs constructed with parameter subsets described by Equation (11).

The selection of the parameter subsets is crucial. On the one hand, the larger interval range and the smaller step of the parameter subset obtain a more accurate result. On the other hand, the larger interval range and the smaller step of the parameter subset cost excessive computation and decrease the efficiency of the method. When both accuracy and efficiency are considered, the subset, S_fc, is determined as [800:10:1200] by inspecting the FFT spectrum of the Doppler-shifted signal. The sound speed is set as S_c = [320:1:360]. The parameter sets S_r, S_X0, and S_V0 are set as [0.5:0.1:1.5], [2:0.1:6] and [20:0.1:40], respectively. Only the first-order derivative of the moving speed, a₁, is considered, and S_a is set as [35:0.5:45].

A searching grid of the model parameters is constructed based on the aforementioned six parameter subsets. Once a group of parameters is determined, the parameters are inputted to the PMDW generator introduced in Section 2.1.1 to generate a PMDW. The correlation filtering analysis introduced in Section 2.1.2 is then performed between the PMDW and the simulated Doppler-shifted signal. Figure 11 shows the maximal correlation coefficients for the different elements from a specified model parameter subset. When the parameter set of the PMDW is determined as

Table 1. Kinematic model parameters subset and the optimal kinematic model parameters of the simulated signal.

Items	fc [Hz]	c [m/s]	r [m]	x0 [m]	v0 [m/s]	a0 [m/s2]
Range	[800,1200]	[320,360]	[0.5,1.5]	[2,6]	[20,40]	[35,45]
Step	10	1	0.1	0.1	0.1	0.5
Optimal value	1000	340	1	4	30	40
True value	1000	340	1	4	30	40
Error	0%	0%	0%	0%	0%	0%

, the maximal correlation coefficient between the PMDW and the simulated Doppler-shifted signal can be obtained.

The identified kinematic model parameters and the Doppler-shifted bearing signal are then inputted to the signal resampler introduced in Section 2.1.3 to eliminate the embedded Doppler effect. First, the amplitude vector of the input Doppler-shifted signal is matched with the delayed-time-vector, t_d^dop = [R₀^dop/c^opt, R₀^dop/c^opt + 1/f_s,…, R₀^dop/c^opt + (N^dop−1)/f_s], where $R_0=\sqrt{r^{op{t}^2}+{X_0}^{op{t}^2}}=\sqrt{1^2+4^2}=4.1231$, and a fitting curve, χ^d, can be obtained by fitting the amplitude vector of the input Doppler-shifted signal with t_d^dop. Through the receive-time-calculator described by Equation (18), the receive-time-vector, t_r^dop, can be calculated with the emit-time-vector, t_e^dop = [0,1/f_s,…, (10000−1)/f_s]. The two time vectors are shown in Figure 12a. After resampling the fitting curve, χ^d, by t_r^dop, we obtain the receive-amplitude-vector, S_r^dop(n). Finally, the obtained S_r^dop(n) is rearranged by matching with t_e^dop and then demodulated by Equation (20). The curve of the demodulation weights is indicated in Figure 12b.

The wave form of the obtained Doppler-free signal is plotted as an overlay on the original bearing signal, shown in Figure 13a, with their FFT spectrums in Figure 13b. The correlation coefficient between the obtained Doppler-free signal and the original bearing signal is 0.9695, which indicates that the Doppler effect is effectively eliminated.

After the Doppler effect is eliminated, the transient model analysis method introduced in Section 2.2 is applied to the Doppler-free signal to detect the defect. The transient model is constructed according to Equation (24). Its parameters require optimization from the sets T, F, and Z. The range of T, F, and Z are set as [500/50,000:1/50,000:2000/50,000], [800:10:1200], and { [0.005:0.001:0.03] ∪ [0.04:0.01:0.1] ∪ [0.2:0.1:0.9]}, respectively. The subset of Z is non-uniform to provide higher resolution at lower damping ratio values so that the efficiency of the method can be retained. The grid of the model parameters is constructed according to F and Z for each element from set T. When a group of parameters is determined, the transient model is constructed according to the procedures introduced in Section 2.2.1. The correlation filtering analysis introduced in Section 2.1.2 is then performed between the transient model and the input Doppler-free bearing signal. Figure 14a shows the maximal correlation coefficients for the different elements from set T. The optimal parameters f_c = 1000 and ζ = 0.05 when T is set at 0.025 s.

As a comparison, the simulated signal before the Doppler effect elimination is also analyzed via the transient model analysis method with the same model parameter subsets. The maximal correlation coefficients for the different elements from set T are shown in Figure 14b. The maximal correlation coefficient is obtained when the impact period is 0.028 s, which is an incorrect value. Comparison of the transmit model analysis results between the signals before and after Doppler effect elimination is shown in

Table 2. Comparison of the transmit model analysis results between the signals before and after Doppler effect elimination.

Items	Coef	fres	ζ	T
Before Doppler effect elimination	0.195	1060	0.019	0.028
After Doppler effect elimination	0.780	1000	0.050	0.025

. Thus, performing the introduced Doppler effect elimination method before the conventional transient model analysis is necessary.

4. Application in Locomotive Bearings Fault Diagnosis

4.1. Experimental Setup and Data Acquisition

Two separate experiments were implemented to achieve the Doppler-shifted acoustic signals from defective locomotive roller bearings and verify the effectiveness of the proposed method. The obtained Doppler-shifted bearing acoustic signals are then analyzed by the proposed diagnosis strategy.

In the first experiment, the acoustic signal of a locomotive roller bearing with a single localized defect was acquired by the test bench in Figure 15a. The test bench consists of a drive motor, two supporting pillow blocks (mounted with normal bearings), a testing bearing loaded on the outer race, and a mechanical radial loader. A 4944-A type of microphone (B&K Company, Nærum, Denmark) was mounted adjacent to the outer race of the tested bearing for acoustic signal measurement. The advanced data acquisition system (DAS) provided by National Instruments (Austin, TX, US) was used to record the data. The specifications of the testing bearing are listed in

Table 3. Specifications of the testing bearing.

Type	NJ(P)3226XI
Diameter of the outer race	250 mm
Diameter of the inner race	130 mm
Pitch diameter (D)	190 mm
Diameter of the roller (d)	32 mm
Number of the roller (z)	14

. The bearings were tested under the radial load of 3 t, the rotation speed of the motor was set at 1430 r/min, and the sampling frequency was set at 50 kHz.

The second experiment can be referred to the model illustrated in Figure 1. As shown in Figure 15b, the experiment consists of a moving vehicle (mounted with the acoustic source), a same microphone, and a same DAS as those in the first experiment. When the car passes by the microphone, a Doppler-shifted signal of the testing bearing can be received and recorded. In this experiment, the car was passing by at an accelerated speed. The sampling frequency was also set at 50 kHz.

A single artificial crack with 0.18 mm width was set by the wire-electrode cutting machine on the outer and inner race, as shown in Figure 16. The Doppler-shifted bearing signals were then acquired by the aforementioned experimental steps. In the subsequent sections, the Doppler-shifted signals with the out-race defect (size: 0.18 mm) and inner-race defect (size: 0.18 mm) are analyzed.

4.2. Results

The wave form of the Doppler-shifted bearing signal with the out-race defect is shown in Figure 17a with a length of 12,000 points. As the original recorded signal is highly distorted by the low-frequency component and our study shows that the low-frequency component will reduce the effectiveness of the following correlation analysis. So a fourth-order Butterworth band-pass filter is employed to reduce this distortion. And our study also shows that if the band of the filter is too narrow, for example just around the structural resonance band, the effectiveness of the following correlation analysis will be reduced as the waveform of the signal will be damaged by the filter. So, the filter band is set as 50 Hz to 5 kHz. The wave form of the filtered signal is presented in Figure 17b with its FFT spectrum in Figure 17c.

The proposed diagnosis strategy is then applied to the filtered Doppler-shifted bearing signal. The signal is first inputted to the Doppler effect eliminator introduced in Section 2.1 to identify the kinematic model parameters and eliminate the embedded Doppler effect. During the kinematic parameter identification procedure, the subset S_fc is determined as [800:10:1600] by inspecting the FFT spectrum in Figure 17c. The sound speed is set as S_c = [320:1:360], determined by the value of 343.2 m/s that stands for the sound speed in standard air atmosphere.

The value of sound speed c is determined by the properties of the medium of air. In the Earth's atmosphere, the chief factor affecting the sound speed is the temperature. In practice, it can be estimated through the value calculated by the following practical formula:

$$c_{\mathit{\text{air}}}=331.3+0.606T$$

(30)

The parameter sets S_r and S_X0 are set as [1.5:0.1:2.5] and [2:0.1:6], respectively, based on a previous coarse manual measurement. As the momentum of the car is large and the time duration of the pass-by movement is extremely short, only the first-order derivative of the moving speed, a₁, is considered. S_V0 and S_a is set as [20:0.1:40] and [0:1:10], respectively, by considering the acceleration performance of the testing car. In parctice, the range of initial speed can also be estimated by the length of the recorded signal and the dimension of the vehicle.

PMDWs are then generated with different parameters from the previously determined parameter subsets by the PMDW generator introduced in Section 2.1.1 to generate a PMDW. The correlation filtering analysis introduced in Section 2.1.2 is then performed between the generated PMDW and the Doppler-shifted signal. Figure 18 shows the maximal correlation coefficients for the different elements from a specified parameter subset. When the parameter set of the PMDW is determined as those in

Table 4. Kinematic model parameter subsets and the optimal kinematic model parameters of the bearing signal with the out-race defect.

Items	fc [Hz]	c [m/s]	r [m]	x0 [m]	v0 [Hz]	a0 [Hz2]
Range	[800:1600]	[320:360]	[1.5:2.5]	[2:6]	[20:40]	[0:10]
Step	10	1	0.1	0.1	0.1	1
Optimal value	1250	339	2	3.9	30.5	3

, the maximal correlation coefficient between the PMDW and the Doppler-shifted signal can be obtained.

The signal resampler introduced in Section 2.1.3 is then employed to eliminate the embedded Doppler effect with the identified kinematic parameters in Table 4. The amplitude vector of the input Doppler-shifted signal is matched with the delayed-time-vector, t_d ^dop, which equals [R₀^dop/c^opt, R₀^dop/c^opt + 1/f_s,…, R₀^dop/c^opt + (N^dop−1)/f_s], where R₀^dop equals $\sqrt{r^{op{t}^2}+{X_0}^{op{t}^2}}=\sqrt{2^2+{3.9}^2}=4.3829$, and a fitting curve, χ^d, can be obtained by fitting the amplitude vector of the input Doppler-shifted signal with t_d ^dop. The receive-time-vector, t_r^dop, can be calculated from the emit-time-vector, t_e^dop = [0,1/f_s,…, (12000−1)/f_s], through the receive-time-calculator described by Equation (18). The two time vectors are shown in Figure 19a. The next step is resampling the fitting curve, χ^d, by t_r^dop to obtain the receive-amplitude-vector, S_r^dop(n). After the amplitude vector rearrangement by matching with t_e^dop and the amplitude demodulation procedure implemented by Equation (20), the curve of the demodulation weights can be obtained (Figure 19b).

The wave form of the obtained Doppler-free signal is shown in Figure 17d, with its FFT spectrum in Figure 17e. Compared with the shape of the FFT spectrum in Figure 17c, the problems of frequency shift and frequency band expansion are clearly solved. When Figure 17d is compared with Figure 17b, the amplitude is also clearly demodulated.

The transient model analysis method introduced in Section 2.2 is then applied to detect the characteristic interval of the Doppler-free fault signal. A periodical transient model with parameters adjustable using Equation (24) is constructed from the sets T = [50/50000:1/50000:600/50000], F = [800:5:2400], and Z = { [0.005:0.001:0.03] ∪ [0.04:0.01:0.1] ∪ [0.2:0.1:0.9]}. When a group of parameters is determined, the transient model is constructed according to the procedures introduced in Section 2.2.1. The correlation filtering analysis introduced in Section 2.1.2 is then performed between the transient model and the input Doppler-free bearing signal.

The outer race characteristic frequency is 138.74 Hz as calculated by Equation (25), and the periodical impact interval is 0.0072 s as calculated by Equation (27). Figure 20a reports the maximal correlation coefficients for the different elements from set T. The optimal impact period is 0.0072 s, which is equal to the real bearing fault-related impact interval. The optimal transient model is shown in Figure 20c.

The signal before the Doppler effect elimination is also analyzed via transient model analysis method with the same model parameter subsets. The maximal correlation coefficients for the different elements from set T are shown in Figure 20d. The maximal correlation coefficient is obtained when the impact period is 0.0083 s, which is incorrect.

The obtained Doppler-shifted bearing signal (Figure 21a) with an inner-race defect is also analyzed to further confirm the effectiveness of the proposed fault diagnosis strategy. Thus, the signal was also first pre-processed with a four-order Butterworth band-pass filter (band: 50 Hz to 5 kHz). The wave form of the filtered signal is shown in Figure 21b, with its FFT spectrum in Figure 21c.

The proposed Doppler effect eliminator is first applied to identify the kinematic model parameters and eliminate the embedded Doppler effect. The subset S_fc is determined as [1400:10:2000] by inspecting the FFT spectrum in Figure 21c. The parameter sets S_c, S_r, S_X0, S_V0, and S_a, are set as [320:1:360], [1.5:0.1:2.5], [2:0.1:6], [20:0.1:40], and [0:1:10], respectively. PMDWs are then generated with different parameters from the determined parameter subsets by the PMDW generator introduced in Section 2.1.1 to generate a PMDW. The correlation filtering analysis introduced in Section 2.1.2 is then performed between the generated PMDWs and the simulated Doppler-shifted signal. Figure 22 shows the maximal correlation coefficients for the different elements from a specified parameter subset. When the parameter set of the PMDW is determined as those in

Table 5. Kinematics model parameters subset and the optimal kinematics model parameters of the bearing signal with the inner-race defect.

Items	fc [Hz]	c [m/s]	r [m]	x0 [m]	v0 [m/s]	a0 [m/s2]
Range	[1400:2000]	[320:360]	[1.5:2.5]	[2:6]	[20,40]	[0,10]
Step	10	1	0.1	0.1	0.1	1
Optimal value	1750	342	1.9	4.1	32	5

, the maximal correlation coefficient between the PMDW and the simulated Doppler-shifted signal can be obtained.

The embedded Doppler effect is then eliminated through the signal resampler introduced in Section 2.1.3. First, by matching the amplitude vector of the input Doppler-shifted signal with the delayed-time-vector, t_d^dop, which equals [R₀^dop/c^opt, R₀^dop/c^opt + 1/f_s,…, R₀^dop/c^opt + (N^dop−1)/f_s], where R₀^dop equals $\sqrt{r^{op{t}^2}+{X_0}^{op{t}^2}}=\sqrt{{1.9}^2+{4.1}^2}=4.5188$, a fitting curve, χ^d, can be obtained. Second, the receive-time-vector, t_r^dop, is calculated from the emit-time-vector, t_e^dop = [0,1/f_s,…, (12000−1)/f_s], through the receive-time-calculator described by Equation (18). The two time vectors are shown in Figure 23a. After resampling the fitting curve, χ^d, by t_r^dop, the receive-amplitude-vector, S_r^dop(n), is obtained. Finally, the amplitude vector is rearranged by matching with t_e^dop and performing the amplitude demodulation procedure through Equation (20). The curve of the demodulation weights is shown in Figure 23b.

The wave form of the obtained Doppler-free signal is illustrated in Figure 21d with its FFT spectrum in Figure 21e. Compared with the shape of the FFT spectrum in Figure 21c and the wave form in Figure 21b, the problems of frequency shift, frequency band expansion, and amplitude modulation are clearly solved.

The transient model analysis method introduced in Section 2.2 is then applied to detect the characteristic interval embedded in the Doppler-free fault signal. A periodical transient model with parameters adjustable using Equation (24) is constructed from the sets T = [50/50000:1/50000:600/50000], F = [1200:5:2200], and Z = { [0.005:0.001:0.03] ∪ [0.04:0.01:0.1] ∪ [0.2:0.1:0.9]}. Figure 22a shows the maximal correlation coefficients for the different elements from set T.

The periodical impact interval is 0.0052 s, calculated by Equations (29) and (31). Figure 24a shows that the optimal impact period equals the real bearing fault-related impact interval. The optimal transient model is shown in Figure 24c.

The Doppler-shifted signal is then directly analyzed by the transient model analysis method with the same model parameter subsets. The maximal correlation coefficients for the different elements from set T are shown in Figure 24d. The maximal correlation coefficient is obtained when the impact period is 0.0068 s, which is an incorrect value.

5. Conclusions

In this paper, a fault diagnosis strategy based on the wayside acoustic monitoring technique is invented for locomotive bearing fault diagnosis. A parametric wavelet called PMDW is introduced and employed to identify the kinematic model parameters based on correlation analysis. A time domain signal resampler is introduced and employed to eliminate the embedded Doppler effect in the recorded bearing acoustic signal. The transient model analysis method is also employed to detect the localized bearing faults after the Doppler effect is eliminated. One of the best benefits of the proposed strategy is that all the kinematic model parameters, including the sound speed and the moving speed of the vehicle, as well as the geometric parameters of the model, can be identified based on the recorded signal itself. Thus, the proposed strategy overcomes the difficulties of kinematic model parameter measurement and is adjustable to different types of passing vehicles. Besides, the embedded Doppler effect can be eliminated through the proposed strategy, paving the way for the conventional invented signal processing methods and feature extraction methods. The performance of the proposed strategy has been evaluated by both simulated and practical Doppler-shifted bearing signals carrying fault information. Given the merits revealed in this study, the proposed fault diagnosis strategy can be widely used in wayside health monitoring systems, particularly in situations when vehicles pass by at high moving speeds and kinematic model parameters are difficult to estimate. The proposed data-driven Doppler effect eliminator is also hopeful to be used in other areas such as acoustic communication techniques and sound field holography for moving vehicles.