# A Heavy-Haul Railway Corrugation Diagnosis Method Based on WPD-ASTFT and SVM

Authors: Binghuan Xiao; Jinzhao Liu; Ziyuan Zhang
Journal: Shock and Vibration
 (2022)

---

## Abstract

Rail corrugation in heavy-haul railway increases the contact forces between the wheel and the rail and deteriorates the rail condition. Severe corrugation affects railway operational safety. Fast diagnosis techniques allow technical personnel to perform timely maintenance and repair, preventing the quick deterioration of rail corrugation. This paper presents a heavy-haul railway corrugation diagnosis method incorporating the time-frequency analysis with machine learning methods. First, the signal is decomposed into several subsignals by wavelet packet decomposition (WPD). The paper proposes an adaptive short-time Fourier transform (ASTFT) and performs the ASTFT on the subsignals to obtain the optimal resolution time-frequency distribution and compute the corresponding entropy. The dimensionality reduction based on mean entropy is then performed for the high-dimensional data. The training and testing samples are classified using Support Vector Machine (SVM). The adaptive short-time Fourier transform (ASTFT) is incorporated with the Renyi entropy and the particle swarm optimization algorithm, which achieves a better aggregation of the time-frequency distribution and reduces the computation cost. Finally, to assist the repair work and estimate the severity of the corrugation section, the corrugation index is proposed. The corrugation indices for the determined corrugation sections are calculated to measure the severity of the corrugation. Experimental studies performed on the axle-box vertical acceleration data collected from the heavy-haul comprehensive inspection train show that the method presented by this paper achieves higher accuracy when compared with conventional feature classification methods for time-frequency analysis. The accuracy of corrugation recognition for the presented method is 93%.


---

## Body

## 1. Introduction

Rail corrugation is an irregular wave-type wear of rail surface. It is a common type of rail wear in the heavy-haul railway. Rail corrugation increases the contact forces between wheel and rail. Since the heavy-haul railway has a higher axle load, abnormal wheel-rail force can damage the rail and vehicle components and cause contact fatigue. Severe corrugation can affect operational safety. Currently, there are no techniques to eliminate rail corrugation. Compared with techniques such as rail lubrication for mitigating friction coefficient and rail vibration absorber, rail grinding is considered as an effective maintenance technique to inhibit quick deterioration of corrugation [1]. Using the rail grinder and manually identifying the corrugation is the most direct maintenance method. However, due to the limitations of maintenance windows and the differences in the operating personnel’s proficiency leading to different assessments, manual identification results in low working efficiency. Identifying and assessing rail corrugation both effectively and in a timely manner becomes a key research topic for researchers and scholars.In 1999, Niu adopted the unsprung-mass acceleration signals of the Polish State Railway network to analyse the power spectrum for different corrugation areas [2]. In 2009, Hory et al. proposed a rail corrugation diagnosis method based on time-frequency analysis [3]. They used the ARCAP method to estimate the frequency of the corrugation and the amplitude of depth of the corrugation. In 2001, C. Mandriota et al. proposed a method to detect and identify corrugation based on the surface texture analysis [4]. In 2015, Li et al. proposed an automatic detection algorithm for corrugation based on the axle-box acceleration measurement for identifying the signature tunes of the wheel-track system [5]. In 2017, Kang et al. proposed a corrugation detection method based on laser imaging techniques [6]. They used multiple sensors in parallel at high sampling frequency to capture the rail profile. The region of interest in the captured image data is sent to the host computer to extract the rail profile data for detecting the rail corrugation. In 2018, Dong et al. proposed a method for rail corrugation detection based on wavelet packet energy entropy [7]. The method performed the four-layer wavelet packet decomposition of the axle-box vibration signal and computed the wavelet pocket energy entropy and wavelet energy of each node. Then, the rail corrugation condition is determined by comparative analysis. Li et al. used 3D structured light and wavelet analysis to detect the rail corrugation [8]. They used 3D structured light scanner scanning a rail segment and analysed the flatness of the rail surface; then, they analysed the corrugation using the time-frequency analysis. Zhu et al. proposed a detection method based on Hilbert-Huang transformation [9]. They established a vehicle-rail vertically coupled system and performed numerical simulation. They then performed the empirical mode decomposition and the Hilbert spectrum time-frequency analysis on the dynamic response data. In 2020, Zhang et al. proposed a corrugation detection method using the parameter optimization variational mode decomposition (VMD) method and combined it with the smooth pseudo-Wigner-Ville distribution (SPWVD) [10]. The decomposed signal is analysed using the SPWVD signal time-frequency analysis method to determine the wavelength and the location of the corrugation.The corrugation detection methods based on the image detection are susceptible to external environment factors, thus tending to cause missed detection. Also, using the time-frequency analysis to extract features from the data is feasible; however, due to the issue of the distribution resolution, this method often results in low detection accuracy. Therefore, we proposed a corrugation diagnosis method for heavy-haul railway based on the adaptive time-frequency analysis combined with machine learning methods. The proposed method uses the vertical axle-box acceleration for the corrugation diagnosis. The axle-box acceleration can reflect the surface condition of the rail because it is similar to the wheelset acceleration under the hypothesis of neglecting the stiffness of roller bearings. The method first uses wavelet packet decomposition (WPD) to decompose the acceleration signal into subsignals. Then, the method performs the adaptive short-time Fourier transform (ASTFT) on the subsignals to acquire a high-resolution time-frequency distribution and compute the entropy of each time-frequency distribution as the index for rail corrugation identification. Finally, after the dimensionality reduction to the feature data, the data is classified using a Support Vector Machine (SVM). The rail corrugation diagnosis method based on WPD-ASTFT and SVM improves the resolution of the time-frequency distribution via adaptive time-frequency analysis. Therefore, the method can achieve higher diagnostic accuracy.


## 2. Theoretical Background

### 2.1. Wavelet Packet Decomposition

Wavelet packet decomposition (WPD) resolves the problems in the wavelet decomposition that the resolution cannot be fine for both high and low-frequency bands [11]. The wavelet packet decomposition can perform more precise decomposition on the original signal for the entire frequency band; that is, it performs continuous decomposition for both high-frequency band and low-frequency band. Thereby, the wavelet packet decomposition can achieve a better frequency resolution for the original signal ft.Assume that the conjugate filterhn satisfies(1)∑hn−2khn−2l=δkl,∑hn=2,k,l∈ℝ.Letgk=−1kh1−k, and the coefficient equations of the wavelet packet decomposition are(2)W2pt=2∑khkWp2t−k,W2p+1t=2∑kgkWp2t−k.The reconstructed coefficient for the wavelet packet decomposition is(3)Wpt=2∑kht−2kW2pk+2∑kgt−2kW2p+1k.The signalft after p layer wavelet packet decomposition acquires N=2p subsignals and satisfies fkt,0≤k≤2p+1−1. The original signal after p layers wavelet packet decomposition can be expressed as(4)ft=∑k=2p−13·2p−1fkt.Through the wavelet packet decomposition, the vertical axle-box acceleration signal can be decomposed into different frequency bands, which can help extract the features of corrugation.


### 2.2. Support Vector Machine

Support Vector Machine classification is illustrated in Figure1. The two classes are separated by a red line in the figure, which is defined as the maximum margin separating the hyperplane. The equation for the hyperplane is expressed as(5)ωTx+b=0,where ω is the normal vector and x is the sample point. The distance from the sample points to the hyperplane is maximized. The distance from a sample point in the space x=x1,x2,…,xn to the hyperplane is(6)ωTx+bω,where(7)ω=ω12+⋯+ωn2.Figure 1
An illustration of the SVM for two classes.For classification of two-class data,(8)ωTx+bω≥d,y=1,ωTx+bω≤−d,y=−1.Rearranging in a compact form,(9)yωTx+b≥1.Therefore,(10)d=yωTx+bω.The Support Vector Machine satisfiesyωTx+b=1. To maximize the distance d, ω is minimized. Therefore, the optimization problem is formed as(11)minfω=12ω2,s.t.giω=1−yiωTxi+b≤0,i=1,2…,n.After solving the optimization problem, the classification result for each sample can be found by substituting the sample vector into the decision function. The decision function is(12)fx=signωTx+b,where the sign⋅ is the signum function. The SVM is used to classify the samples into two classes, that is, the corrugation sections and the normal sections.


## 2.1. Wavelet Packet Decomposition

Wavelet packet decomposition (WPD) resolves the problems in the wavelet decomposition that the resolution cannot be fine for both high and low-frequency bands [11]. The wavelet packet decomposition can perform more precise decomposition on the original signal for the entire frequency band; that is, it performs continuous decomposition for both high-frequency band and low-frequency band. Thereby, the wavelet packet decomposition can achieve a better frequency resolution for the original signal ft.Assume that the conjugate filterhn satisfies(1)∑hn−2khn−2l=δkl,∑hn=2,k,l∈ℝ.Letgk=−1kh1−k, and the coefficient equations of the wavelet packet decomposition are(2)W2pt=2∑khkWp2t−k,W2p+1t=2∑kgkWp2t−k.The reconstructed coefficient for the wavelet packet decomposition is(3)Wpt=2∑kht−2kW2pk+2∑kgt−2kW2p+1k.The signalft after p layer wavelet packet decomposition acquires N=2p subsignals and satisfies fkt,0≤k≤2p+1−1. The original signal after p layers wavelet packet decomposition can be expressed as(4)ft=∑k=2p−13·2p−1fkt.Through the wavelet packet decomposition, the vertical axle-box acceleration signal can be decomposed into different frequency bands, which can help extract the features of corrugation.


## 2.2. Support Vector Machine

Support Vector Machine classification is illustrated in Figure1. The two classes are separated by a red line in the figure, which is defined as the maximum margin separating the hyperplane. The equation for the hyperplane is expressed as(5)ωTx+b=0,where ω is the normal vector and x is the sample point. The distance from the sample points to the hyperplane is maximized. The distance from a sample point in the space x=x1,x2,…,xn to the hyperplane is(6)ωTx+bω,where(7)ω=ω12+⋯+ωn2.Figure 1
An illustration of the SVM for two classes.For classification of two-class data,(8)ωTx+bω≥d,y=1,ωTx+bω≤−d,y=−1.Rearranging in a compact form,(9)yωTx+b≥1.Therefore,(10)d=yωTx+bω.The Support Vector Machine satisfiesyωTx+b=1. To maximize the distance d, ω is minimized. Therefore, the optimization problem is formed as(11)minfω=12ω2,s.t.giω=1−yiωTxi+b≤0,i=1,2…,n.After solving the optimization problem, the classification result for each sample can be found by substituting the sample vector into the decision function. The decision function is(12)fx=signωTx+b,where the sign⋅ is the signum function. The SVM is used to classify the samples into two classes, that is, the corrugation sections and the normal sections.


## 3. WPD-ASTFT and SVM Diagnosis Method

We proposed the WPD-ASTFT and SVM diagnosis method to determine the rail corrugation. The diagnosis flow is presented in Figure2. For the short-time Fourier transform, only a fixed window length can be selected for the signal with multifrequency components at a one-time interval, which causes the insufficiency of the time-frequency distribution resolution of the signal. We first used WPD to decompose the components with large frequency differences into different subsignals so that for different subsignals we can choose different window lengths for STFT at the same time interval. We then performed the ASTFT for each subsignal. The optimal window length is obtained by calculating the window length that achieves the best aggregation of the time-frequency distribution. We used Renyi entropy to measure the aggregation of time-frequency distribution; the smaller the entropy, the better aggregation of the time-frequency distribution. The Renyi entropy is also used as the fitness function in the particle swarm optimization (PSO) algorithm. The PSO algorithm is computationally faster than other optimization algorithms, such as the genetic algorithm, when solving the optimization problem, so the PSO algorithm is used when solving the optimal window length. Each subsignal is transformed through ASTFT to obtain the time-frequency distribution corresponding to its lowest Renyi entropy, and this Renyi entropy value is used as the classification feature. The mean feature dimension reduction method is proposed to reduce the high-dimensional Renyi entropy value obtained by each segment of the signal after the WPD-ASTFT. Finally, the inspection data are classified into normal and corrugation classes using the SVM classifier. For the corrugation sections, we further calculated the corrugation index to determine the degree of corrugation.Figure 2
Flowchart of the corrugation detection method.


### 3.1. Adaptive Short-Time Fourier Transform

The short-time Fourier transform (STFT) is a conventional time-frequency analysis method, and its window length has a significant impact on the resolution of the time-frequency distribution. To avoid manual adjustment of the window length and ensure that the signal has a better resolution time-frequency distribution after STFT, we proposed an adaptive short-time Fourier transform (ASTFT) combined with the Renyi entropy and particle swarm optimization. We transformed each subsignal using ASTFT and calculated the Renyi entropy of the time-frequency distribution under the corresponding optimal window length.For a series of signals containingn information, y1,y2,y3⋯yn, the probability for each information is py1,py2,py3⋯pyn. They form a system S:(13)S=YP=y1y2⋯ynpy1py2⋯pyn.The Renyi entropy for this system is(14)HY=11−qln∑i=1npyiq,where q>0 and q≠1. Usually, q=2.In this paper, the Renyi entropy is adopted as the criterion to evaluate the concentration of the time-frequency representation. To get the best time-frequency resolution after STFT, the optimal window length for STFT is needed for the subsignal. Therefore, the Renyi entropy of the time-frequency distributionSTFTt,f obtained by STFT with different window lengths is calculated using the following equation:(15)Eς=11−qln∫0T∫0∞STFTb,ξ2qdξdb∫0T∫0∞STFTb,ξ2dξdbq,where b and ξ are the dummy variables representing the time and frequency, respectively, in the STFT.Particle swarm optimization (PSO) is a computational optimization method and is a population iteration process in which particles search for the optimal particles in the search space. Compared with the genetic algorithm, PSO is simple to implement and easy to adjust parameters and has fast convergence. In D-dimensional space, the position for a particlei is xi=xi1,xi2,…,xiD, substituting the position xi into the fitness function fx for solving the fitness value. The velocity for a particle i is vi=vi1,vi2,…,viD. Thus, the velocity update equation for the particle i in the d th dimension is(16)vidk=wvidk−1+c1r1pbestid−xidk−1+c2r2gbstd−xidk−1.The position update equation for the particlei in the d th dimension is(17)xidk=xidk−1+vidk−1,where vidk is the velocity vector component in the d th dimension for particle i in the k th iteration, xidk is the position vector component in the d th dimension for particle i in the k th iteration, c1 and c2 are the acceleration coefficients for cognitive and social components, respectively, to adjust the learning step, r1 and r2 are random numbers between 0 and 1, and w is the nonnegative inertia weight to adjust the search space. The particle i obtains the optimal position pbesti=pi1,pi2,…,piD, and the population obtains the optimal position gbest=g1,g2,…,gD.The PSO algorithm flow is shown in Figure3.Figure 3
The particle swarm optimization for finding the best window length for STFT.The different window lengths for the subsignals are set as the particles for optimization, and the Renyi entropy for the time-frequency distribution after the STFT is set as the fitness function. The Renyi entropy is computed at the minimal fitness value.The signal after wavelet packet decomposition obtainsN subsignals. The frequency for each subsignal has certain differences. The optimal window length for the STFT is found using the PSO. For the signal xt and Gaussian window g, the STFT is(18)Gjt,f=∫xjξg∗ξ−te−i2πftdξ,where j=1,2,…,N, gt is the optimal window length, and f is the frequency.After the adaptive short-time Fourier transform, the time-frequency distribution for the minimal Renyi entropy is obtained for theN subsignals under the corresponding optimal window length. Also, the Renyi entropy value in the transform can be used as the resolution criterion and the feature of the data in the classification.


### 3.2. Feature Extraction and Dimensionality Reduction

We propose a dimensionality reduction method to extract the features of the N-column entropy values and to reduce the dimension of the entropy data. Dimensionality reduction of data is based on the mean difference of columnwise entropy values between the two classes of data. The data columns with the entropy values that differ greatly between the two classes are selected as the data features. The computation procedure is as follows:(1)
Computing the mean valueComputing the columnwise means of the entropy values for the acceleration data of the normal rail segments,(19)mi=∑j=1JNjiJ,wherei=1,2,…,N is the index of columns, J is the number of normal segments, and Nj×i is the entropy matrix.Computing the columnwise means of the entropy values for the acceleration of the corrugation rail segments,(20)mi′=∑k=1KCkiK,wherei=1,2,…,N is the index of columns, K is the number of corrugation segments, and Ck×i is the entropy matrix.(2)
Computing the entropy value differencesiSubtracting by elements of the columnwise means between the normal and corrugation segment and taking the absolute difference,(21)si=mi−mi′.(3)
Computing the average values¯ of si as the characteristic entropy,(22)s¯=∑i=1NsiN.(4)
Finding the location wheresi>s¯,i=1,2,…,N and assuming the indices as x,y,…,z,1≤x,y,z≤N, all the data columns with indices x,y,…,z are kept.


### 3.3. WPD-ASTFT and SVM Diagnosis Method

The steps for the heavy-haul railway corrugation diagnosis method based on WPD-ASTFT and SVM are described as follows:(1)
Dividing the axle-box vertical acceleration signal of the heavy-haul railway into units by 50 meters(2)
Performing the wavelet packet decomposition (WPD) for each unit signal to get several subsignals(3)
Performing the adaptive short-time Fourier transform (ASTFT) for every subsignal in each unit and computing the Renyi entropy corresponding to each time-frequency distribution(4)
Performing the dimensionality reduction for the entropy values(5)
Using the data to train the SVM classifier and making the predictionsThe schematic for the diagnosis procedure is shown in Figure4.Figure 4
The schematic for the WPD-ASTFT diagnosis method.


### 3.4. Corrugation Index

After the corrugation sections are determined using the SVM classifier, the severity of the corrugation sections is further determined using the corrugation index. Based on the simulation, the amplitude of the acceleration vibration is related to both the corrugation depth and the wavelength. When the wavelength is fixed, the amplitude of the vibration increases with the increased depth. With the increase of the corrugation wavelength, the vibration amplitude reaches a maximum value and then decreases with the increase of the wavelength. Therefore, in a real situation, it is difficult to determine the corrugation depth only from the acceleration amplitude. Since the change of speed in the heavy-duty railway is small and the wavelength of corrugation is concentrated in a small range of wavelength from 200 mm to 300 mm, the severity of the corrugation could be determined from the aspect of the vibration energy.For the determined corrugation sections, the moving RMS of the axle-box vertical acceleration is calculated. Then, the corrugation index is calculated (i.e., the mean over time of the RMS for each section) to measure the severity of corrugation.The moving RMS is calculated using(23)Rt=∑t=rr+Kf2tK,where K is the forward window length.The corrugation index is calculated using(24)C=∑t=0TRtT,where T is the duration of the corrugation section.


## 3.1. Adaptive Short-Time Fourier Transform

The short-time Fourier transform (STFT) is a conventional time-frequency analysis method, and its window length has a significant impact on the resolution of the time-frequency distribution. To avoid manual adjustment of the window length and ensure that the signal has a better resolution time-frequency distribution after STFT, we proposed an adaptive short-time Fourier transform (ASTFT) combined with the Renyi entropy and particle swarm optimization. We transformed each subsignal using ASTFT and calculated the Renyi entropy of the time-frequency distribution under the corresponding optimal window length.For a series of signals containingn information, y1,y2,y3⋯yn, the probability for each information is py1,py2,py3⋯pyn. They form a system S:(13)S=YP=y1y2⋯ynpy1py2⋯pyn.The Renyi entropy for this system is(14)HY=11−qln∑i=1npyiq,where q>0 and q≠1. Usually, q=2.In this paper, the Renyi entropy is adopted as the criterion to evaluate the concentration of the time-frequency representation. To get the best time-frequency resolution after STFT, the optimal window length for STFT is needed for the subsignal. Therefore, the Renyi entropy of the time-frequency distributionSTFTt,f obtained by STFT with different window lengths is calculated using the following equation:(15)Eς=11−qln∫0T∫0∞STFTb,ξ2qdξdb∫0T∫0∞STFTb,ξ2dξdbq,where b and ξ are the dummy variables representing the time and frequency, respectively, in the STFT.Particle swarm optimization (PSO) is a computational optimization method and is a population iteration process in which particles search for the optimal particles in the search space. Compared with the genetic algorithm, PSO is simple to implement and easy to adjust parameters and has fast convergence. In D-dimensional space, the position for a particlei is xi=xi1,xi2,…,xiD, substituting the position xi into the fitness function fx for solving the fitness value. The velocity for a particle i is vi=vi1,vi2,…,viD. Thus, the velocity update equation for the particle i in the d th dimension is(16)vidk=wvidk−1+c1r1pbestid−xidk−1+c2r2gbstd−xidk−1.The position update equation for the particlei in the d th dimension is(17)xidk=xidk−1+vidk−1,where vidk is the velocity vector component in the d th dimension for particle i in the k th iteration, xidk is the position vector component in the d th dimension for particle i in the k th iteration, c1 and c2 are the acceleration coefficients for cognitive and social components, respectively, to adjust the learning step, r1 and r2 are random numbers between 0 and 1, and w is the nonnegative inertia weight to adjust the search space. The particle i obtains the optimal position pbesti=pi1,pi2,…,piD, and the population obtains the optimal position gbest=g1,g2,…,gD.The PSO algorithm flow is shown in Figure3.Figure 3
The particle swarm optimization for finding the best window length for STFT.The different window lengths for the subsignals are set as the particles for optimization, and the Renyi entropy for the time-frequency distribution after the STFT is set as the fitness function. The Renyi entropy is computed at the minimal fitness value.The signal after wavelet packet decomposition obtainsN subsignals. The frequency for each subsignal has certain differences. The optimal window length for the STFT is found using the PSO. For the signal xt and Gaussian window g, the STFT is(18)Gjt,f=∫xjξg∗ξ−te−i2πftdξ,where j=1,2,…,N, gt is the optimal window length, and f is the frequency.After the adaptive short-time Fourier transform, the time-frequency distribution for the minimal Renyi entropy is obtained for theN subsignals under the corresponding optimal window length. Also, the Renyi entropy value in the transform can be used as the resolution criterion and the feature of the data in the classification.


## 3.2. Feature Extraction and Dimensionality Reduction

We propose a dimensionality reduction method to extract the features of the N-column entropy values and to reduce the dimension of the entropy data. Dimensionality reduction of data is based on the mean difference of columnwise entropy values between the two classes of data. The data columns with the entropy values that differ greatly between the two classes are selected as the data features. The computation procedure is as follows:(1)
Computing the mean valueComputing the columnwise means of the entropy values for the acceleration data of the normal rail segments,(19)mi=∑j=1JNjiJ,wherei=1,2,…,N is the index of columns, J is the number of normal segments, and Nj×i is the entropy matrix.Computing the columnwise means of the entropy values for the acceleration of the corrugation rail segments,(20)mi′=∑k=1KCkiK,wherei=1,2,…,N is the index of columns, K is the number of corrugation segments, and Ck×i is the entropy matrix.(2)
Computing the entropy value differencesiSubtracting by elements of the columnwise means between the normal and corrugation segment and taking the absolute difference,(21)si=mi−mi′.(3)
Computing the average values¯ of si as the characteristic entropy,(22)s¯=∑i=1NsiN.(4)
Finding the location wheresi>s¯,i=1,2,…,N and assuming the indices as x,y,…,z,1≤x,y,z≤N, all the data columns with indices x,y,…,z are kept.


## 3.3. WPD-ASTFT and SVM Diagnosis Method

The steps for the heavy-haul railway corrugation diagnosis method based on WPD-ASTFT and SVM are described as follows:(1)
Dividing the axle-box vertical acceleration signal of the heavy-haul railway into units by 50 meters(2)
Performing the wavelet packet decomposition (WPD) for each unit signal to get several subsignals(3)
Performing the adaptive short-time Fourier transform (ASTFT) for every subsignal in each unit and computing the Renyi entropy corresponding to each time-frequency distribution(4)
Performing the dimensionality reduction for the entropy values(5)
Using the data to train the SVM classifier and making the predictionsThe schematic for the diagnosis procedure is shown in Figure4.Figure 4
The schematic for the WPD-ASTFT diagnosis method.


## 3.4. Corrugation Index

After the corrugation sections are determined using the SVM classifier, the severity of the corrugation sections is further determined using the corrugation index. Based on the simulation, the amplitude of the acceleration vibration is related to both the corrugation depth and the wavelength. When the wavelength is fixed, the amplitude of the vibration increases with the increased depth. With the increase of the corrugation wavelength, the vibration amplitude reaches a maximum value and then decreases with the increase of the wavelength. Therefore, in a real situation, it is difficult to determine the corrugation depth only from the acceleration amplitude. Since the change of speed in the heavy-duty railway is small and the wavelength of corrugation is concentrated in a small range of wavelength from 200 mm to 300 mm, the severity of the corrugation could be determined from the aspect of the vibration energy.For the determined corrugation sections, the moving RMS of the axle-box vertical acceleration is calculated. Then, the corrugation index is calculated (i.e., the mean over time of the RMS for each section) to measure the severity of corrugation.The moving RMS is calculated using(23)Rt=∑t=rr+Kf2tK,where K is the forward window length.The corrugation index is calculated using(24)C=∑t=0TRtT,where T is the duration of the corrugation section.


## 4. Results and Experimental Verification

The experimental signal is the axle-box vertical acceleration signal collected by the Chinese heavy-haul railway comprehensive inspection train. The heavy-haul railway inspection train is shown in Figure5(a), and the axle-box acceleration sensor is shown in Figure 5(b). The sensor is a biaxial accelerometer and can measure the acceleration in the vertical direction and the axial direction. The sampling rate for the acceleration signal is 2000 Hz. The axle-box acceleration signal is divided into 50-meter units. One unit axle-box acceleration signal is decomposed into several subsignals with different frequency bands after the wavelet packet decomposition, and the original signal can be reconstructed using the subsignals xt=x1t+⋯+xNt, where N is the number of subsignals. In this paper, 8 subsignals are obtained using the three-layer wavelet packet decomposition, and the 6th, 7th, and 8th subsignals are superimposed as the final 6th subsignal. This processing can reduce the computation time and fully decompose the axle-box acceleration signal.Figure 5
Photo of the inspection train and the sensor location. (a) The heavy-haul railway comprehensive inspection train. (b) The installed location of the axle-box acceleration sensor.
(a)(b)


### 4.1. Adaptive Time-Frequency Feature Analysis and Feature Extraction for the Acceleration Signal

The examples in Figure6 show the acceleration signal of axle-box in normal heavy-haul railway section and in corrugation section. Figure 6(a) shows the acceleration signal in the normal rail section for the mileage from K401 + 050 to K401 + 100. Figure 6(b) shows the acceleration signal in the corrugation rail section for the mileage from K387 + 150 to K387 + 200. After three-layer wavelet packet decomposition and superposition, the original signal for each section is decomposed into 6 subsignals. The subsignals for the two sections are presented in Figure 7.Figure 6
The plots of the axle-box acceleration signal for (a) the normal rail sections and (b) the corrugation rail section.
(a)(b)Figure 7
The plots of the decomposed subsignals of the acceleration signal for (a) the normal rail sections and (b) the corrugation rail section.
(a)(b)Adaptive short-time Fourier transform is applied to the subsignals of the normal rail section and the corrugation section. The results for one normal section and one corrugation section are presented as an example in Figure8 for illustration. Comparing the time-frequency distribution of the 1st subsignal between the corrugation section and normal section, it can be found that the aggregation of the frequency for the time-frequency distribution of the axle-box acceleration is significantly higher in the corrugation section than in the normal section. The Renyi entropy values corresponding to the time-frequency distribution for the example are calculated according to equation (15), and the values are presented in Table 1.Figure 8
The plots of the time-frequency distribution for the six subsignals for (a) the normal rail sections and (b) the corrugation rail section.
(a)(b)Table 1
The Renyi entropy values for a normal section and a corrugation section.
Sample sectionsSubsignal 1Subsignal 2Subsignal 3Subsignal 4Subsignal 5Subsignal 6Normal section13.4317.1819.3118.6417.7319.77Corrugation section15.7416.8518.5818.3317.7419.25The Renyi entropy corresponding to the 6 subsignals is obtained by performing the WPD-ASTFT to the acceleration signals for all the sections. In Table2, column 1 to column 75 show the Renyi entropy for the acceleration signal in the normal section, and column 76 to column 100 show the Renyi entropy for the corrugation section. The mean of the Renyi entropy value and the characteristic entropy are calculated according to equations (19) to (22). The mean of the entropy values for the normal section and corrugation section are tabulated in Table 3, and the characteristic entropy is calculated as s¯=0.47.Table 2
The Renyi entropy values for the experimental data.
SamplesSubsignal 1Subsignal 2Subsignal 3Subsignal 4Subsignal 5Subsignal 6Sample 113.7916.3418.5518.0317.7519.10Sample 213.7016.1819.2018.4817.4119.40Sample 313.8916.7016.6115.9717.3218.24……Sample 7413.8915.7415.4816.9017.8518.70Sample 7513.6618.0418.9019.3718.1619.56Sample 7615.5716.5018.5218.0617.7619.49Sample 7715.7416.8518.5818.3317.7419.25Sample 7815.0516.2418.3217.8916.7419.06……Sample 9915.3916.4718.3317.8617.7019.00Sample 10015.9117.2118.4718.1317.4219.14Table 3
The columnwise means of the Renyi entropy and the difference of the means for the normal section and corrugation section.
Subsignal 1Subsignal 2Subsignal 3Subsignal 4Subsignal 5Subsignal 6Mean for normal sectionmi15.4816.7118.3118.1017.3719.02Mean for corrugation sectionmi′13.6116.7018.8218.3517.4719.11Difference of the means1.870.010.510.250.100.08The means of the entropy values for each subsignal of the normal section and corrugation section are plotted in Figure9(a). The differences in the mean values between the normal section and corrugation section are plotted in Figure 9(b). We used a threshold on the difference of the mean entropy to select the features of the entropy. Using the average of the differences 0.47 as the threshold, the entropies of the 1st subsignal and the 3rd subsignal are kept as the first feature entropy and the second feature entropy, respectively.Figure 9
(a) The plot of the columnwise mean values of the entropy for each subsignal. (b) The difference of the means between the normal section and the corrugation section.
(a)(b)


### 4.2. Classification Results

The first feature entropy and the second feature entropy obtained after the dimensionality reduction are used to train the SVM classifier. The sample 1 to sample 75 in Table2 represent the normal section and are labelled as “1,” whereas sample 76 to sample 100 in Table 2 represent the corrugation section and are labelled as “−1.” The training samples are numbered from 1 to 100 and are used to train the SVM classifier. The two classes are separated by the hyperplane as shown in Figure 10. The testing results of the classification of the normal section and corrugation section using the SVM classifier are plotted in Figure 11. The actual labels and prediction labels of the test data are shown in Table 4. For the 15 sections, 14 sections are predicted accurately, with a test accuracy as high as 93.33%.Figure 10
The training results for the Support Vector Machine classier.Figure 11
The prediction results using the trained Support Vector Machine classifier.Table 4
The corrugation index for the detected corrugation section in the testing samples.
Section12345Corrugation index1.751.751.41.591.43


### 4.3. On-Site Review

The rail section K199 + 350 to K199 + 400 is diagnosed as having corrugation using the WPD-ASTFT and SVM-based corrugation diagnostic method. This conclusion is consistent with the actual situation. Also, the time-frequency distribution with the minimum entropy of the subsignal is concentrated at 61 Hz, with a train speed of 65 kilometres per hour, which corresponds to a corrugation wavelength of 296 millimetres. The rail surface condition during the on-site review in this section is shown in Figure12. The surface roughness of the rail measured by the digital roughness gauge is plotted as Figure 13(a). Using the measured data to compute the power spectrum shown in Figure 13(b), the calculated corrugation wavelength is 297 millimetres.Figure 12
The surface condition of the rail for the mileage section K199 + 357.Figure 13
(a) The surface roughness of the rail surface. (b) The power spectrum.
(a)(b)


### 4.4. Corrugation Index

After the corrugation sections are filtered out by the SVM classifier, the corrugation index for each detected section is calculated to determine the severity of the corrugation problem. The corrugation indices for the testing sections are tabulated in Table4. The maximum value of the corrugation index is 1.75 corresponding to Section 1 and Section 2. The actual corrugation situation for Section 1 and Section 2 is shown in Figures 14(a) and 14(b), respectively. It can be seen from the pictures that the corrugation problem is serious for those two sections, and the depth for the two sections all has exceeded 1 mm.Figure 14
The actual corrugation situation for the corrugation sections in the testing samples. The pictures (a)–(e) show the corrugation situation for the testing Section1 to Section 5.
(a)(b)(c)(d)(e)


### 4.5. Method Comparison

We compared our method with another method using a CEEMD-STFT time-frequency analysis [12] to compute the classification features and repeated the analysis using the same data. The accuracy for the CEEMD-STFT method is 86.67%, whereas the accuracy for the proposed WPD-ASTFT method is 93.33%. The results using the CEEMD-STFT method are shown in Figure 15, and the prediction results’ comparison is tabulated in Table 5. This illustrates the advantage of the adaptive time-frequency analysis algorithms when computing the classification features.Figure 15
The prediction result of the SVM classifier using CEEMD-STFT method.Table 5
The prediction results for the testing samples using the conventional method and the presented method.
Testing samples123456789101112131415Actual label1111111111−1−1−1−1−1Predicted label using CEEMD-STFT method11111111−11−1−1−11−1Predicted label using WPD-ASTFT and SVM11111111−11−1−1−1−1−1When using the WPD-ASTFT method to calculate the entropy of the time-frequency distribution, the time-frequency distribution is obtained by an adaptive process under the corresponding optimal window length instead of selecting a fixed window length. Also, the optimal window length corresponding to different frequency components is different, which makes the time-frequency distribution of each signal more aggregated than the result of the CEEMD-STFT methods. Therefore, when calculating the entropy value, the advantage of the WPD-ASTFT method can improve the differentiation in the feature extraction for the normal rail section and the corrugation rail section.


## 4.1. Adaptive Time-Frequency Feature Analysis and Feature Extraction for the Acceleration Signal

The examples in Figure6 show the acceleration signal of axle-box in normal heavy-haul railway section and in corrugation section. Figure 6(a) shows the acceleration signal in the normal rail section for the mileage from K401 + 050 to K401 + 100. Figure 6(b) shows the acceleration signal in the corrugation rail section for the mileage from K387 + 150 to K387 + 200. After three-layer wavelet packet decomposition and superposition, the original signal for each section is decomposed into 6 subsignals. The subsignals for the two sections are presented in Figure 7.Figure 6
The plots of the axle-box acceleration signal for (a) the normal rail sections and (b) the corrugation rail section.
(a)(b)Figure 7
The plots of the decomposed subsignals of the acceleration signal for (a) the normal rail sections and (b) the corrugation rail section.
(a)(b)Adaptive short-time Fourier transform is applied to the subsignals of the normal rail section and the corrugation section. The results for one normal section and one corrugation section are presented as an example in Figure8 for illustration. Comparing the time-frequency distribution of the 1st subsignal between the corrugation section and normal section, it can be found that the aggregation of the frequency for the time-frequency distribution of the axle-box acceleration is significantly higher in the corrugation section than in the normal section. The Renyi entropy values corresponding to the time-frequency distribution for the example are calculated according to equation (15), and the values are presented in Table 1.Figure 8
The plots of the time-frequency distribution for the six subsignals for (a) the normal rail sections and (b) the corrugation rail section.
(a)(b)Table 1
The Renyi entropy values for a normal section and a corrugation section.
Sample sectionsSubsignal 1Subsignal 2Subsignal 3Subsignal 4Subsignal 5Subsignal 6Normal section13.4317.1819.3118.6417.7319.77Corrugation section15.7416.8518.5818.3317.7419.25The Renyi entropy corresponding to the 6 subsignals is obtained by performing the WPD-ASTFT to the acceleration signals for all the sections. In Table2, column 1 to column 75 show the Renyi entropy for the acceleration signal in the normal section, and column 76 to column 100 show the Renyi entropy for the corrugation section. The mean of the Renyi entropy value and the characteristic entropy are calculated according to equations (19) to (22). The mean of the entropy values for the normal section and corrugation section are tabulated in Table 3, and the characteristic entropy is calculated as s¯=0.47.Table 2
The Renyi entropy values for the experimental data.
SamplesSubsignal 1Subsignal 2Subsignal 3Subsignal 4Subsignal 5Subsignal 6Sample 113.7916.3418.5518.0317.7519.10Sample 213.7016.1819.2018.4817.4119.40Sample 313.8916.7016.6115.9717.3218.24……Sample 7413.8915.7415.4816.9017.8518.70Sample 7513.6618.0418.9019.3718.1619.56Sample 7615.5716.5018.5218.0617.7619.49Sample 7715.7416.8518.5818.3317.7419.25Sample 7815.0516.2418.3217.8916.7419.06……Sample 9915.3916.4718.3317.8617.7019.00Sample 10015.9117.2118.4718.1317.4219.14Table 3
The columnwise means of the Renyi entropy and the difference of the means for the normal section and corrugation section.
Subsignal 1Subsignal 2Subsignal 3Subsignal 4Subsignal 5Subsignal 6Mean for normal sectionmi15.4816.7118.3118.1017.3719.02Mean for corrugation sectionmi′13.6116.7018.8218.3517.4719.11Difference of the means1.870.010.510.250.100.08The means of the entropy values for each subsignal of the normal section and corrugation section are plotted in Figure9(a). The differences in the mean values between the normal section and corrugation section are plotted in Figure 9(b). We used a threshold on the difference of the mean entropy to select the features of the entropy. Using the average of the differences 0.47 as the threshold, the entropies of the 1st subsignal and the 3rd subsignal are kept as the first feature entropy and the second feature entropy, respectively.Figure 9
(a) The plot of the columnwise mean values of the entropy for each subsignal. (b) The difference of the means between the normal section and the corrugation section.
(a)(b)


## 4.2. Classification Results

The first feature entropy and the second feature entropy obtained after the dimensionality reduction are used to train the SVM classifier. The sample 1 to sample 75 in Table2 represent the normal section and are labelled as “1,” whereas sample 76 to sample 100 in Table 2 represent the corrugation section and are labelled as “−1.” The training samples are numbered from 1 to 100 and are used to train the SVM classifier. The two classes are separated by the hyperplane as shown in Figure 10. The testing results of the classification of the normal section and corrugation section using the SVM classifier are plotted in Figure 11. The actual labels and prediction labels of the test data are shown in Table 4. For the 15 sections, 14 sections are predicted accurately, with a test accuracy as high as 93.33%.Figure 10
The training results for the Support Vector Machine classier.Figure 11
The prediction results using the trained Support Vector Machine classifier.Table 4
The corrugation index for the detected corrugation section in the testing samples.
Section12345Corrugation index1.751.751.41.591.43


## 4.3. On-Site Review

The rail section K199 + 350 to K199 + 400 is diagnosed as having corrugation using the WPD-ASTFT and SVM-based corrugation diagnostic method. This conclusion is consistent with the actual situation. Also, the time-frequency distribution with the minimum entropy of the subsignal is concentrated at 61 Hz, with a train speed of 65 kilometres per hour, which corresponds to a corrugation wavelength of 296 millimetres. The rail surface condition during the on-site review in this section is shown in Figure12. The surface roughness of the rail measured by the digital roughness gauge is plotted as Figure 13(a). Using the measured data to compute the power spectrum shown in Figure 13(b), the calculated corrugation wavelength is 297 millimetres.Figure 12
The surface condition of the rail for the mileage section K199 + 357.Figure 13
(a) The surface roughness of the rail surface. (b) The power spectrum.
(a)(b)


## 4.4. Corrugation Index

After the corrugation sections are filtered out by the SVM classifier, the corrugation index for each detected section is calculated to determine the severity of the corrugation problem. The corrugation indices for the testing sections are tabulated in Table4. The maximum value of the corrugation index is 1.75 corresponding to Section 1 and Section 2. The actual corrugation situation for Section 1 and Section 2 is shown in Figures 14(a) and 14(b), respectively. It can be seen from the pictures that the corrugation problem is serious for those two sections, and the depth for the two sections all has exceeded 1 mm.Figure 14
The actual corrugation situation for the corrugation sections in the testing samples. The pictures (a)–(e) show the corrugation situation for the testing Section1 to Section 5.
(a)(b)(c)(d)(e)


## 4.5. Method Comparison

We compared our method with another method using a CEEMD-STFT time-frequency analysis [12] to compute the classification features and repeated the analysis using the same data. The accuracy for the CEEMD-STFT method is 86.67%, whereas the accuracy for the proposed WPD-ASTFT method is 93.33%. The results using the CEEMD-STFT method are shown in Figure 15, and the prediction results’ comparison is tabulated in Table 5. This illustrates the advantage of the adaptive time-frequency analysis algorithms when computing the classification features.Figure 15
The prediction result of the SVM classifier using CEEMD-STFT method.Table 5
The prediction results for the testing samples using the conventional method and the presented method.
Testing samples123456789101112131415Actual label1111111111−1−1−1−1−1Predicted label using CEEMD-STFT method11111111−11−1−1−11−1Predicted label using WPD-ASTFT and SVM11111111−11−1−1−1−1−1When using the WPD-ASTFT method to calculate the entropy of the time-frequency distribution, the time-frequency distribution is obtained by an adaptive process under the corresponding optimal window length instead of selecting a fixed window length. Also, the optimal window length corresponding to different frequency components is different, which makes the time-frequency distribution of each signal more aggregated than the result of the CEEMD-STFT methods. Therefore, when calculating the entropy value, the advantage of the WPD-ASTFT method can improve the differentiation in the feature extraction for the normal rail section and the corrugation rail section.


## 5. Conclusions

Combined with the time-frequency analysis method and machine learning techniques, the paper proposes a diagnosis method of rail corrugation based on WPD-ASTFT and SVM. The method is applied to diagnose the rail corrugation in the heavy-haul railway, and the reliability of the method is verified by comparing the predicted results with the actual situation. By comparing the results from the methods of extracting classification features by the CEEMD-STFT time-frequency analysis, the accuracy of the prediction results presented by this paper is up to 93.33%, which is better than the CEEMD-STFT methods. The severity level of the detected corrugation sections can be further measured by the corrugation indices. The results can effectively guide the maintenance and repair work to inhibit the deterioration of rail corrugation [13–18].



---

*Source: 8370796-2022-01-29.xml*