Quantitative Detection of Vertical Track Irregularities under Non-Stationary Conditions with Variable Vehicle Speed

Track irregularities directly affect the quality and safety of railway vehicle operations. Quantitative detection and real-time monitoring of track irregularities are of great importance. However, due to the frequent variable vehicle speed, vehicle operation is a typical non-stationary process. The traditional signal analysis methods are unsuitable for non-stationary processes, making the quantitative detection of the wavelength and amplitude of track irregularities difficult. To solve the above problems, this paper proposes a quantitative detection method of track irregularities under non-stationary conditions with variable vehicle speed by order tracking analysis for the first time. Firstly, a simplified wheel–rail dynamic model is established to derive the quantitative relationship between the axle-box vertical vibration and the track vertical irregularities. Secondly, the Simpson double integration method is proposed to calculate the axle-box vertical displacement based on the axle-box vertical acceleration, and the process error is optimized. Thirdly, based on the order tracking analysis theory, the angular domain resampling is performed on the axle-box vertical displacement time-domain signal in combination with the wheel rotation speed signals, and the quantitative detection of the track irregularities is achieved. Finally, the proposed method is validated based on simulation and field test analysis cases. We provide theoretical support and method reference for the quantitative detection method of track irregularities.


Introduction
Track irregularity is one of the main reasons for vibration in railway vehicle operation [1].When these vibrations are unacceptable, several problems will follow, such as riding comfort reduction, running stability decline, wheel-rail noise pollution, fatigue strength failure of components, etc. [2][3][4][5].Therefore, track irregularity should be controlled and managed as much as possible to create a better vehicle operating environment.
The current hot research topic of track irregularity can be summarized as the following aspects: 1  ⃝ Track infrastructure evolution-due to the coupling effect of temperature, humidity, load, and other factors, the railway infrastructure undergoes the evolution of physical and chemical properties.The track infrastructure presents overall or local cumulative deformation, which affects the track irregularity state [6][7][8]. 2 ⃝ Track-vehicle interaction-taking the track and vehicle as a whole system, the long-term interaction between the track and the vehicle will affect the track geometry.The deterioration of track irregularity will make the vehicle system produce a more significant response, accelerating the deterioration of track and vehicle components [9][10][11]. 3 ⃝ Track management and maintenance-however, before the causing mechanisms of the problem mentioned above are thoroughly clarified, and some basic solutions are proposed, scientifically managing all kinds of railway (high-speed, heavy-haul, urban subways, etc.) operating lines and Sensors 2024, 24, 3804 2 of 30 maintaining a safe, stable, and reliable railway vehicle service environment is a crucial problem [12,13].
The track management and maintenance prerequisite is to obtain actual track irregularity state data.The selection of measurement methods greatly affects the authenticity of the results. 1  ⃝ Static measurement-is directly measured by hand or using some small devices to obtain the track irregularity state data for a short distance.Since the force acting on the track by the measuring device is relatively small, the elastic deformation of the rail, bed, roadbed, etc., is negligible.The measured track irregularity data is geometric irregularity at this time, which is called static irregularity. 2 ⃝ Dynamic measurement-however, when the vehicles are running on the track, the enormous dynamic load can cause comprehensive elastic deformation of the rail, bed, and roadbed, called dynamic irregularity.Strictly speaking, compared with static irregularity, dynamic irregularity is more reflective of the actual vehicle working environment because it includes both the geometric irregularity before loading and track elastic deformation after loading [14][15][16].Obviously, dynamic irregularity can not be obtained by direct measurement.Currently, two main indirect methods for detecting dynamic vertical track irregularity are as follows: the chord measurement method and the inertial reference method. 1 ⃝ The chord measurement method-takes the fixed chord length as the measurement reference and equal interval sampling is used to obtain the chord measurement value of each sampling point.Its mathematical calculation model is as follows [17]: where k is the mileage position of the sample point, y(k) is the chord measurement value at the mileage position k, p is the half chord length, and h(k), h(k − p), h(k + p) are the vertical irregularity amplitudes of the track at the mileage positions k, k − p, k + p.
Performing the discrete Fourier transform to Equation (1), the transfer function H(ω) can be expressed as: where ω is the spatial angular frequency (ω = 2π/λ), and λ is the track irregularity wavelength).
As shown in Equation ( 2), H(ω) is the transfer function related to the spatial angular frequency, independent of the vehicle running speed and time.Therefore, there is no nonstationary problem.However, the value of the transfer function H(ω) is varied with ω in the range [0, 2] and cos(ωp) is in the range [−1, 1] [18].Especially when the value of the transfer function H(ω) is close to 0, it is easy to be disturbed by the noise signal, and the track irregularity state near its corresponding wavelength is difficult to detect, leading to a lack of stability in the measurement accuracy.In addition, the chord measurement method is only applied to a specific rail inspection vehicle.With the increasing growth of rail transportation tasks, the free time of rail inspection vehicles for line inspection is insufficient, and it is impossible to achieve real-time monitoring.Hence, the chord measurement method's shortcomings in the economy and efficiency are gradually exposed [19][20][21].

2
⃝ The axle-box on-board measurement method-due to the small elastic deformation caused by the interaction between the wheel and rail, the wheel and rail can be considered rigid bodies, and there is no suspension system between the wheelset and axle-box.All the track irregularity can be reflected in the axle-box vibration signal [22][23][24][25].Therefore, evaluating the track irregularity state by collecting axle-box vibration information online has become a hot research topic in recent years.Kawasaki et al. obtained the acceleration vibration features of an axle-box by establishing a rail-vehicle system simulation model and then judged the track irregularity state through the features identification [26]; M. Bocciolon et al. developed a real-time track condition monitoring technique by calculating the mean square value of axle-box vibration acceleration to analyze the track irregularity state [27]; Yang et al.
Sensors 2024, 24, 3804 3 of 30 applied wavelet transform theory to eliminate the axle-box vibration acceleration trend term, performing the double integration to the acceleration signal and evaluating the track irregularity state based on the obtained axle-box displacement signal [28]; Ning et al. combined the improved empirical mode decomposition (EMD) with Cohen's kernel, and applied it to the analysis of the axle-box non-stationary vibration signal to achieve online monitoring of the track irregularity fault location [29]; Zhang et al. designed a refined system transfer function identification test method by combining the iterative analysis of simulation and test to reproduce the track irregularity based on the axle-box vibration acceleration [30]; Sun et al. used the inertial reference method and digital filtering method to analyze the axle-box vibration acceleration of high-speed trains to achieve online monitoring of track irregularity [17].However, the above studies were limited in solving the fundamental problem of non-stationary vibration signals, and have not been able to achieve quantitative detection of the essential parameters of the track irregularity.
This paper proposes a complete set of quantitative detection methods for track irregularities under variable vehicle speed conditions.In Section 1, a brief overview of the current research about track irregularities is given.In Section 2, a simplified wheelrail dynamics model is established, and the quantitative relationship between axle-box vertical vibration frequency, amplitude, and track irregularity is demonstrated by theoretical derivation.In Section 3, the axle-box vertical vibration displacement calculation method based on Simpson time-domain double integration is proposed, and the improved Hodrick-Prescott decomposition method is used to optimize the process error.Section 4 proposes the quantitative detection method of track vertical irregularity under vehicle variable speed conditions based on the order tracking analysis theory.In Sections 5 and 6, the correctness of the above quantitative detection method is verified based on simulation and field test cases, respectively.The conclusion is given in Section 7. The overall technical route is shown in Figure 1.
ulation model and then judged the track irregularity state through the features identifica tion [26]; M. Bocciolon et al. developed a real-time track condition monitoring techniqu by calculating the mean square value of axle-box vibration acceleration to analyze th track irregularity state [27]; Yang et al. applied wavelet transform theory to eliminate th axle-box vibration acceleration trend term, performing the double integration to the acce eration signal and evaluating the track irregularity state based on the obtained axle-bo displacement signal [28]; Ning et al. combined the improved empirical mode decompos tion (EMD) with Cohen's kernel, and applied it to the analysis of the axle-box non-station ary vibration signal to achieve online monitoring of the track irregularity fault locatio [29]; Zhang et al. designed a refined system transfer function identification test metho by combining the iterative analysis of simulation and test to reproduce the track irregu larity based on the axle-box vibration acceleration [30]; Sun et al. used the inertial refe ence method and digital filtering method to analyze the axle-box vibration acceleration o high-speed trains to achieve online monitoring of track irregularity [17].However, th above studies were limited in solving the fundamental problem of non-stationary vibra tion signals, and have not been able to achieve quantitative detection of the essential pa rameters of the track irregularity.
This paper proposes a complete set of quantitative detection methods for track irreg ularities under variable vehicle speed conditions.In Section 1, a brief overview of the cu rent research about track irregularities is given.In Section 2, a simplified wheel-rail dy namics model is established, and the quantitative relationship between axle-box vertica vibration frequency, amplitude, and track irregularity is demonstrated by theoretical de ivation.In Section 3, the axle-box vertical vibration displacement calculation metho based on Simpson time-domain double integration is proposed, and the improved Ho drick-Prescott decomposition method is used to optimize the process error.Section 4 pro poses the quantitative detection method of track vertical irregularity under vehicle varia ble speed conditions based on the order tracking analysis theory.In Sections 5 and 6, th correctness of the above quantitative detection method is verified based on simulation an field test cases, respectively.The conclusion is given in Section 7. The overall technica route is shown in Figure 1.

The Quantitative Relationship between the Axle-Box Vertical Vibration and Vertical Track Irregularity
To derive the quantitative relationship between the axle-box's vertical vibration an vertical track irregularity, a simplified vertical wheel-rail dynamics model is established as shown in Figure 2. Due to the minimal elastic deformation caused by the interactio between the wheel and rail, the dynamic model can be considered a rigid body [31,32 Because there is no suspension system between the wheelset and axle-box, the track irreg ularity characteristics can be reflected in the vibration signal of the axle-box.

The Quantitative Relationship between the Axle-Box Vertical Vibration and Vertical Track Irregularity
To derive the quantitative relationship between the axle-box's vertical vibration and vertical track irregularity, a simplified vertical wheel-rail dynamics model is established, as shown in Figure 2. Due to the minimal elastic deformation caused by the interaction between the wheel and rail, the dynamic model can be considered a rigid body [31,32].Because there is no suspension system between the wheelset and axle-box, the track irregularity characteristics can be reflected in the vibration signal of the axle-box.In Figure 2, v is the vehicle running speed, m1 is the unsprun the steel rail, k1 is the elastic contact stiffness between the wheel stiffness of the rail, c2 is the vertical damping of the rail, z1 is the the wheelset/axle-box, z2 is the vertical displacement of the steel displacement excitation.The directions to the right and down ar The motion differential equation of the above dynamic mod When the initial displacement and velocity are both 0, the La tion (3) can be written as: The irregularity excitation can be seen as a linear superposit waves with different wavelengths, amplitudes, and phases.The ir of a single harmonic, as shown in Figure 3, can be expressed as a where a is the wave depth of the irregularity and ω is the time an In Figure 2, v is the vehicle running speed, m 1 is the unsprung mass, m 2 is the mass of the steel rail, k 1 is the elastic contact stiffness between the wheel and rail, k 2 is the vertical stiffness of the rail, c 2 is the vertical damping of the rail, z 1 is the vertical displacement of the wheelset/axle-box, z 2 is the vertical displacement of the steel rail, and u is the vertical displacement excitation.The directions to the right and down are positive here.
The motion differential equation of the above dynamic model is expressed as: When the initial displacement and velocity are both 0, the Laplace transform of Equation (3) can be written as: The irregularity excitation can be seen as a linear superposition of multiple harmonic waves with different wavelengths, amplitudes, and phases.The irregularity excitation u(t) of a single harmonic, as shown in Figure 3, can be expressed as a cosine function: where a is the wave depth of the irregularity and ω is the time angular frequency.In Figure 2, v is the vehicle running speed, m1 is the unsprung the steel rail, k1 is the elastic contact stiffness between the wheel a stiffness of the rail, c2 is the vertical damping of the rail, z1 is the the wheelset/axle-box, z2 is the vertical displacement of the steel r displacement excitation.The directions to the right and down are The motion differential equation of the above dynamic mode When the initial displacement and velocity are both 0, the Lap tion (3) can be written as: The irregularity excitation can be seen as a linear superpositi waves with different wavelengths, amplitudes, and phases.The ir of a single harmonic, as shown in Figure 3, can be expressed as a where a is the wave depth of the irregularity and ω is the time an  By the Laplace transformation, Equation ( 5) can be rewritten as: Sensors 2024, 24, 3804 5 of 30 Then, to find the transfer function of the responses z 1 and z 2 corresponding to the irregularity excitation u, assume the following relationship: where z 1 (s) and z 2 (s) are the Laplace transforms corresponding to z 1 (t) and z 2 (t), and H z1 (s) and H z2 (s) are the transfer functions corresponding to the irregularity excitation u(s) in response to z 1 (s) and z 2 (s).By combining Equations ( 4) and ( 7), the transfer function H z1 (s) can be solved: Comparing the oscillation characteristic terms in Equations ( 12) and ( 13), the amplitude of the vertical acceleration of the axle-box is ω 2 times the vertical displacement with a phase difference π.According to the vertical displacement z 1 (t) of the axle-box and the frequency response function H z1 (jω), the vertical track irregularity excitation can be calculated, as shown in Equations ( 7) and (8).With the vehicle running speed v known, the vertical vibration frequency of the axle-box can be used to estimate the wavelength of the track irregularity.The vertical vibration amplitude D of the axle-box acceleration or displacement can be used to estimate the wave depth of the track irregularity.

Time-Domain Track Vertical Irregularity Displacement Excitation
Equation ( 12) can be written as: ..
where D is the amplitude of axle-box vertical vibration acceleration, which is a constant factor related to ω 1 , ω 2 , c 2 , k 1 , k 2 , and ω. ψ is the phase and Q 0 is the constant term.Then, the vertical vibration velocity and displacement of the axle-box can be obtained by the first and second times integrals of Equation ( 14), expressed as: .
Comparing the oscillation characteristic term in Equations ( 15) and ( 16), the amplitude of the vertical acceleration of the axle-box is ω 2 times the vertical displacement with a phase difference π.This characteristic is entirely consistent with the conclusion in Section 2, which indicates that it is theoretically possible to use the double integration method to calculate the vertical oscillation displacement of the axle-box.Then, combining the frequency response function H z1 (jω), the track irregularity excitation can be calculated.

Improved Simpson Time-Domain Integration Method
There are two methods for calculating displacement by using integrated acceleration: the frequency-domain method and the time-domain method.However, the railway vehicle speed in operation frequently changes.The vibration signal is the typical non-stationary signal.The traditional discrete Fourier transform will lead to frequency ambiguity in the frequency domain, leading to serious estimation errors.Therefore, the frequency domain integration method is unsuitable for analyzing and processing this situation.However, the time-domain integration can avoid frequency distortion, frequency ambiguity, and other problems caused by discrete Fourier transform [33,34].The Simpson integral formula can be expressed as: where a and b represent the lower and upper limits of the integration interval, respectively.Substituting the discrete acceleration time series a(n), (n = 1, 2, . .., N 0 ) into Equation ( 18) can obtain the velocity time series v(n), (n = 1, 2, 3, . .., N 0 − 1).
where ∆t is the time interval.
Furthermore, the integral remainder error of Simpson formula can be expressed as: Equation (20) shows that the residual error of the Simpson formula is related to (b − a).The absolute value of the error decreases as the integration range [a, b] between adjacent points decreases.Therefore, we propose to add a nonlinear interpolation to the initial signal before integration to reduce the time interval of the signal and, thus, improve the integration accuracy.This idea is validated in Section 3.3.

Trend Term Decomposition and Elimination
In addition to the integration errors, the errors caused by the trend term should be given more attention.Since the initial value of the acceleration signal is affected by uncertain factors during the acquisition process, the signal is often accompanied by phenomena such as zero drift and noise interference.This error information accumulates and amplifies during integration, affecting the calculation results.Removing these trend terms is a critical issue in time-domain integration operations.
The main trend term elimination methods are polynomial fitting, CEEMDAN empirical mode decomposition, and Hodrick-Prescott trend decomposition.The research on polynomial fitting and CEEMDAN empirical mode decomposition is quite mature, and their theories and methods can be seen in references [35][36][37], which do not need to be elaborated further.This paper will focus on the theory and application of the Hodrick-Prescott trend decomposition method and compare it with the above two methods to verify its correctness.
Assuming there is a trend term g(n) and a residual term c(n) in the signal a(n): where a(n) represents the initial signal, g(n) represents the trend term, and c(n) represents the residual term, which contains the fluctuation and noise terms.Then, establish a loss function, which is expressed as: where λ is the smoothing factor.As shown in Equation ( 22), the first term on the right is essentially the principle of least squares, characterizing the closeness between the initial signal a(n) and the trend term g(n).The second term on the right can be viewed as a second-order difference equation for g(n), i.e., ∆g(n) − ∆g(n − 1) = [g(n + 1) − g(n)] − [g(n) − g(n − 1)], characterizing the smoothness of the desired trend term g(n).When the loss function M takes the minimum value, it means that the resulting trend term g(n) is both closest to the trend of the initial signal a(n) and sufficiently smooth.
Calculate the first order derivatives of Equation (22) to g(1); g(2); . .., respectively, and set them equal to 0: It can be expressed as matrix form: where A is the initial signal vector, G is the trend term vector, C is the residual term vector, and F is the coefficient matrix, expressed as Equation (25).
Therefore, the formula for the trend term vector G can be expressed as: where I is the unit vector.
As shown in Equation ( 26), the trend term g(n) can be separated from the initial signal a(n) as long as the smoothing factor λ takes an appropriate value.When the loss function M is a number very close to 0, the optimal smoothing factor λ can be expressed as: where σ c 2 is the variance of the residual term c(n) and σ 2 ∆ 2 g is the variance of the secondorder difference of the trend term g(n).
Since the trend term g(n) and the residual term c(n) are unknown before the calculation, the smoothing factor λ is also unknown.Here, we propose an optimal smoothing factor λ search method, as the following technical route: (a) consider the smoothing factor as a coefficient to be determined and set the initial value as λ j (j = 0); (b) substitute λ j into Equation (26) to separate the trend term g j (n) and the residual term c j (n); (c) calculate the variance σ 2 ∆ 2 g of the second-order difference of the trend term g j (n) and the variance σ cj 2 of the residual term cj(n), respectively; (d) calculate the equivalent error ε j , (e) reset λ j = λ j + j∆λ (j = 1, 2, . ..) and repeat the steps (a) to (d); (f) find λ j corresponding to the smallest equivalent error ε j , which is the optimal smoothing factor λ opt .
In addition, it is worth noting that the optimal smoothing factor λ opt determined by the above technical route must be based on certain assumptions, namely that the trend term g(n) and the residual term c(n) are subordinate to independent normal distributions.Meanwhile, the second-order differential equation for g(n) in Equation (22) shows that the determination of the smoothing factor is closely related to the sampling time interval.Therefore, it has been proposed that the smoothing factor should be corrected according to the sampling frequency (or sampling interval).
However, the correction method is controversial, and the formula has no exact derivation and proof [38][39][40].In this paper, a series of simulated signals were cyclically tested.The fitted trend curves were analyzed for different values of the smoothing factor.It was proposed that the fitted trend curve for a vibration signal with a sampling frequency of 5120 Hz and λ = 1.47 × 10 12 has the smallest error with the simulated trend curve, as shown in Figure 4.

Validation and Analysis
To verify the correctness of the above time-domain integra decomposition method, assume a simulation vibration acceleratio above methods to calculate the vibration displacement.The sim perimposed by five signals.a1(t), a2(t), and a3(t) represent the vibra ated by short, medium, and long wave periodic track irregularitie ( ) ( ) ( ) ( ) ( ) ( ) Finally, add a Gaussian white noise signal a5(t) into the signa its variance is 100.Perform the analogue-to-digital conversion to obtain the discrete signal a(n).The sampling time is 2 s and the sa Hz.
The three trend term decomposition methods are used for th Figure 5.All three methods as a whole can fit the trend term wel mial fit decomposition method has a poor fit superiority at the bre tation function (e.g., 0 s, 1.0 s).The CEEMDAN empirical modal has serious trend errors at the end, caused by the inherent defect decomposition class of methods (i.e., the endpoint error effect).Prescott decomposition method has the most stable overall perfo

Validation and Analysis
To verify the correctness of the above time-domain integration method and trend decomposition method, assume a simulation vibration acceleration signal a(t) and use the above methods to calculate the vibration displacement.The simulation signal a(t) is superimposed by five signals.a 1 (t), a 2 (t), and a 3 (t) represent the vibration acceleration generated by short, medium, and long wave periodic track irregularities features, respectively: where a 1 (t) = 30•cos(2π•50t + 60) and a 2 (t) = 10•cos(2π•6t + 120); a 3 (t) = 2•cos(2π•2t + 201).Then, to express the non-stationary characteristics, assume the signal a(t) contains an initial trend term a 4 (t): Finally, add a Gaussian white noise signal a 5 (t) into the signal a(t); its mean is 0, and its variance is 100.Perform the analogue-to-digital conversion to the analogue signal and obtain the discrete signal a(n).The sampling time is 2 s and the sampling frequency is 500 Hz.
The three trend term decomposition methods are used for the analysis, as shown in Figure 5.All three methods as a whole can fit the trend term well.However, the polynomial fit decomposition method has a poor fit superiority at the breakpoints of the segmentation function (e.g., 0 s, 1.0 s).The CEEMDAN empirical modal decomposition method has serious trend errors at the end, caused by the inherent defect of the empirical modal decomposition class of methods (i.e., the endpoint error effect).However, the Hodrick-Prescott decomposition method has the most stable overall performance.This is because the Hodrick-Prescott decomposition method contains both the least squares fitting idea, which makes the fitted trend term closest to the initial signal, and the fluctuation term constraint, which makes the fitted trend term smooth enough.
has serious trend errors at the end, caused by the inherent defect of the empirical modal decomposition class of methods (i.e., the endpoint error effect).However, the Hodrick-Prescott decomposition method has the most stable overall performance.This is because the Hodrick-Prescott decomposition method contains both the least squares fitting idea, which makes the fitted trend term closest to the initial signal, and the fluctuation term constraint, which makes the fitted trend term smooth enough.As shown in Figure 6, both the single-segment and segmented trend curves show good adaptability.Therefore, the Hodrick-Prescott decomposition method is also used in the subsequent integration calculation for trend elimination of the curves before and after integration.Then, as shown in Equations ( 17)-( 19), the Simpson double integration method can be used to calculate the axle-box vertical displacement based on the acceleration signal.It is necessary to perform three trend term eliminations on the signal before the first integration, before the second integration, and after the second integration.Taking the simulation acceleration signal a(t) as an example, the results obtained by the first and second direct Simpson integrations are shown in Figure 7 for the velocity signal and the displacement signal, respectively.The time-domain integration process produces the trend terms, and the velocity and displacement signals after removing the trend term have a reduced error with the theoretical signal.As shown in Figure 6, both the single-segment and segmented trend curves show good adaptability.Therefore, the Hodrick-Prescott decomposition method is also used in the subsequent integration calculation for trend elimination of the curves before and after integration.As shown in Figure 6, both the single-segment and segmented trend curve good adaptability.Therefore, the Hodrick-Prescott decomposition method is also the subsequent integration calculation for trend elimination of the curves before an integration.Then, as shown in Equations ( 17)-( 19), the Simpson double integration meth be used to calculate the axle-box vertical displacement based on the acceleration si is necessary to perform three trend term eliminations on the signal before the first i tion, before the second integration, and after the second integration.Taking the sim acceleration signal a(t) as an example, the results obtained by the first and second Simpson integrations are shown in Figure 7 for the velocity signal and the displa signal, respectively.The time-domain integration process produces the trend term the velocity and displacement signals after removing the trend term have a reduce with the theoretical signal.Then, as shown in Equations ( 17)-( 19), the Simpson double integration method can be used to calculate the axle-box vertical displacement based on the acceleration signal.It is necessary to perform three trend term eliminations on the signal before the first integration, before the second integration, and after the second integration.Taking the simulation acceleration signal a(t) as an example, the results obtained by the first and second direct Simpson integrations are shown in Figure 7 for the velocity signal and the displacement signal, respectively.The time-domain integration process produces the trend terms, and the velocity and displacement signals after removing the trend term have a reduced error with the theoretical signal.
tion, before the second integration, and after the second integration.Taking the simulation acceleration signal a(t) as an example, the results obtained by the first and second direct Simpson integrations are shown in Figure 7 for the velocity signal and the displacement signal, respectively.The time-domain integration process produces the trend terms, and the velocity and displacement signals after removing the trend term have a reduced error with the theoretical signal.Meanwhile, since the error is related to the size of the integration region [a, b] between adjacent points, theoretically, the error is smaller when [a, b] is smaller.We propose Meanwhile, since the error is related to the size of the integration region [a, b] between adjacent points, theoretically, the error is smaller when [a, b] is smaller.We propose to reduce the sampling interval of the signal by adding a spline interpolation step before the integration and then improve the accuracy of the time-domain integration.As shown in Figure 8, the Simpson integration error can indeed be effectively reduced by reducing the sampling time interval of the signal, improving the accuracy of the time-domain integration.
Sensors 2024, 24, 3804 11 of 30 to reduce the sampling interval of the signal by adding a spline interpolation step before the integration and then improve the accuracy of the time-domain integration.As shown in Figure 8, the Simpson integration error can indeed be effectively reduced by reducing the sampling time interval of the signal, improving the accuracy of the time-domain integration.To quantitatively evaluate the accuracy of the interpolated and uninterpolated methods before time-domain integration, the degree of agreement between the theoretical displacement signal and the displacement signal obtained by integration is reflected by calculating the correlation index R.The closer the value of R is to 1, the higher the accuracy of the integration calculation [41].The correlation index R is 0.7343 for the displacements calculated directly without interpolation before integration, while the correlation index R is 0.9210 for the displacements calculated after interpolation before integration, and the correlation is significantly improved.
In summary, it is theoretically possible to calculate the axle-box vertical displacement signal by double integration of the axle-box vibration acceleration.Considering that the measured signal is non-stationary, the time-domain integration method should be used, and the effects of integration error and trend term error should be focused on.Finally, the time-domain track vertical irregularity displacement excitation can be calculated according to the axle-box vertical displacement and the frequency response function H(ω).

Quantitative Detection of Vertical Track Irregularities Based on the Order Tracking Analysis Theory
Railway vehicles always operate at variable speeds, causing the axle-box vertical vibration acceleration to become the typical non-stationary signal.Hence, the time-domain track vertical irregularity displacement excitation calculated according to the method in Section 3 is the non-stationary signal.Traditional signal processing methods cannot detect the wavelength and amplitude of track irregularities from non-stationary signals.But, if these time-domain signals can be converted into the spatial domain signal, the problem of non-stationary signals will be completely solved.Considering the wheel as the rotating mechanical component, we propose the order tracking analysis theory for the non-stationary signal processing for the first time.This section will focus on the detailed derivation and verification of applying this theory in the quantitative detection of vertical track ir- To quantitatively evaluate the accuracy of the interpolated and uninterpolated methods before time-domain integration, the degree of agreement between the theoretical displacement signal and the displacement signal obtained by integration is reflected by calculating the correlation index R.The closer the value of R is to 1, the higher the accuracy of the integration calculation [41].The correlation index R is 0.7343 for the displacements calculated directly without interpolation before integration, while the correlation index R is 0.9210 for the displacements calculated after interpolation before integration, and the correlation is significantly improved.
In summary, it is theoretically possible to calculate the axle-box vertical displacement signal by double integration of the axle-box vibration acceleration.Considering that the measured signal is non-stationary, the time-domain integration method should be used, and the effects of integration error and trend term error should be focused on.Finally, the time-domain track vertical irregularity displacement excitation can be calculated according to the axle-box vertical displacement and the frequency response function H(ω).

Quantitative Detection of Vertical Track Irregularities Based on the Order Tracking Analysis Theory
Railway vehicles always operate at variable speeds, causing the axle-box vertical vibration acceleration to become the typical non-stationary signal.Hence, the time-domain track vertical irregularity displacement excitation calculated according to the method in Section 3 is the non-stationary signal.Traditional signal processing methods cannot detect the wavelength and amplitude of track irregularities from non-stationary signals.
But, if these time-domain signals can be converted into the spatial domain signal, the problem of non-stationary signals will be completely solved.Considering the wheel as the rotating mechanical component, we propose the order tracking analysis theory for the non-stationary signal processing for the first time.This section will focus on the detailed derivation and verification of applying this theory in the quantitative detection of vertical track irregularities.

Introduction and Derivation of Order Tracking Analysis Theory
Order tracking analysis is an effective method for extracting fault features under variable speed conditions.The key is to measure and analyze the wheel rotational speed or angular velocity, then resample the time-domain signal non-uniformly by using an equal-angle sample method.Therefore, the non-stationary signal in the time-domain is transformed into a stationary signal in the angular domain (spatial domain).Finally, the resampled angular-domain signal is analyzed by order spectrum to detect the fault features of the components quantitatively.
The order is a signal characterization parameter related to wheel rotational speed, expressed as: where o is the order (dimensionless), f is the vibration frequency, R is the rotational speed of the wheel (revolutions per minute, rpm), and R = 60•ω/2π; ω is the angular velocity of the wheel.
As shown in Equation ( 30), the order tracking analysis is essentially a kind of frequency analysis.However, the frequency information is not fixed; it changes with the rotational speed change.The larger the order o, the higher the corresponding vibration frequency f, which reflects the vibration characteristics of short wave track irregularities.When the order o is smaller or even the fractional (a number less than 1), the corresponding vibration frequency f is lower, which reflects the vibration characteristics of long-wave vertical track irregularities.
A simple case is used to illustrate how to convert the time-domain track vertical irregularity displacement excitation into a spatial signal in the angular domain, and finally, the order spectrum analysis is performed.Let a track vertical irregularity displacement excitation be z(t), as shown in Equation ( 31) and Figure 9.The sampling frequency is 1024 Hz and the total sampling length is 10 s.
Sensors 2024, 24, 3804 12 of 30 resampled angular-domain signal is analyzed by order spectrum to detect the fault features of the components quantitatively.The order is a signal characterization parameter related to wheel rotational speed, expressed as: where o is the order (dimensionless), f is the vibration frequency, R is the rotational speed of the wheel (revolutions per minute, rpm), and R = 60•ω/2π; ω is the angular velocity of the wheel.
As shown in Equation ( 30), the order tracking analysis is essentially a kind of frequency analysis.However, the frequency information is not fixed; it changes with the rotational speed change.The larger the order o, the higher the corresponding vibration frequency f, which reflects the vibration characteristics of short wave track irregularities.When the order o is smaller or even the fractional (a number less than 1), the corresponding vibration frequency f is lower, which reflects the vibration characteristics of long-wave vertical track irregularities.
A simple case is used to illustrate how to convert the time-domain track vertical irregularity displacement excitation into a spatial signal in the angular domain, and finally, the order spectrum analysis is performed.Let a track vertical irregularity displacement excitation be z(t), as shown in Equation (31) and Figure 9.The sampling frequency is 1024 Hz and the total sampling length is 10 s. ) The signal z(t) contains two vibration features z1(t) and z2(t), having the characteristic of increasing vibration frequency with time, which is typical of a non-stationary signal.These vibration features in Figure 9 are similar to the axle-box vertical irregularity displacement excitation generated when rail vehicles accelerate on tracks.z1(t) and z2(t) correspond to the wavelength of l1 and l2 of the two kinds of periodic vertical track irregularity.If the discrete Fourier transform is directly applied to estimate the irregularity characteristics (wavelength, amplitude) from the signal, z(t) will cause serious errors, even incorrect results.As shown in Figure 10, it exhibits severe frequency ambiguity with characteristic frequencies ranging from 0 to 500 Hz.The wavelength and amplitude of irregularities are completely unrecognizable.As shown in Figure 11, When the time-frequency analysis method is used, the change in signal frequency with time is analyzed by a shorttime Fourier transform.It exhibits the trend of frequency variation in the time-frequency plot.However, the magnitude of the track irregularity is reflected by the depth of the color and does not quantitatively express the degree of track irregularity.In addition, the frequency ridges are also very blurred (thick), which will cause significant errors in estimating the frequency (wavelength).Therefore, traditional frequency domain analysis meth- The signal z(t) contains two vibration features z 1 (t) and z 2 (t), having the characteristic of increasing vibration frequency with time, which is typical of a non-stationary signal.These vibration features in Figure 9 are similar to the axle-box vertical irregularity displacement excitation generated when rail vehicles accelerate on tracks.z 1 (t) and z 2 (t) correspond to the wavelength of l 1 and l 2 of the two kinds of periodic vertical track irregularity.
If the discrete Fourier transform is directly applied to estimate the irregularity characteristics (wavelength, amplitude) from the signal, z(t) will cause serious errors, even incorrect results.As shown in Figure 10, it exhibits severe frequency ambiguity with characteristic frequencies ranging from 0 to 500 Hz.The wavelength and amplitude of irregularities are completely unrecognizable.As shown in Figure 11, When the time-frequency analysis method is used, the change in signal frequency with time is analyzed by a shorttime Fourier transform.It exhibits the trend of frequency variation in the time-frequency plot.However, the magnitude of the track irregularity is reflected by the depth of the color and does not quantitatively express the degree of track irregularity.In addition, the frequency ridges are also very blurred (thick), which will cause significant errors in estimating the frequency (wavelength).Therefore, traditional frequency domain analysis methods are unsuitable for the above non-stationary signals.However, the above problems will be very clear when the order tracking ana used.Basically, the vertical displacement excitation frequency is related to the speed and the wavelength, i.e., f = v/l.Since the vertical track irregularity l is fix unchanged, the ultimate influence on the vibration frequency f is the operating speed v. Due to the mutual conversion relationship between the vehicle running s and the wheel rotation speed R, the complete derivation process is as follows: where f is the vertical displacement excitation frequency (axle-box vertical vibrat quency), v is the vehicle running speed (wheel longitudinal speed), ω is the wheel a velocity, r is the wheel radius, l is the vertical track irregularity wavelength, C is th circumference, o is the order (the ratio of the wheel circumference C to the track irregularity wavelength l, C/l), and R is the wheel rotation speed.
Therefore, according to the simulated signal in Equation ( 31), its wheel rot  However, the above problems will be very clear when the order tracking an used.Basically, the vertical displacement excitation frequency is related to the speed and the wavelength, i.e., f = v/l.Since the vertical track irregularity l is fi unchanged, the ultimate influence on the vibration frequency f is the operating speed v. Due to the mutual conversion relationship between the vehicle running and the wheel rotation speed R, the complete derivation process is as follows: where f is the vertical displacement excitation frequency (axle-box vertical vibra quency), v is the vehicle running speed (wheel longitudinal speed), ω is the wheel velocity, r is the wheel radius, l is the vertical track irregularity wavelength, C is th circumference, o is the order (the ratio of the wheel circumference C to the track irregularity wavelength l, C/l), and R is the wheel rotation speed.
Therefore, according to the simulated signal in Equation ( 31), its wheel ro However, the above problems will be very clear when the order tracking analysis is used.Basically, the vertical displacement excitation frequency is related to the vehicle speed and the wavelength, i.e., f = v/l.Since the vertical track irregularity l is fixed and unchanged, the ultimate influence on the vibration frequency f is the operating vehicle speed v. Due to the mutual conversion relationship between the vehicle running speed v and the wheel rotation speed R, the complete derivation process is as follows: where f is the vertical displacement excitation frequency (axle-box vertical vibration frequency), v is the vehicle running speed (wheel longitudinal speed), ω is the wheel angular velocity, r is the wheel radius, l is the vertical track irregularity wavelength, C is the wheel circumference, o is the order (the ratio of the wheel circumference C to the track vertical irregularity wavelength l, C/l), and R is the wheel rotation speed.
Therefore, according to the simulated signal in Equation (31), its wheel rotational speed can be calculated: Since the vibration response characteristics of z 1 (t) and z 2 (t) are generated by the same wheel in operation, the wheel rotational speed must be the same, i.e., R(t) = R 1 (t) = R 2 (t).Different vibration characteristics actually correspond to different track irregularity characteristics.As shown in Equations ( 33) and (34), o 1 = o 2 /4, which is expressed as the wavelength l 1 is four times l 2 , i.e., l 1 = 4l 2 .In this case, assuming the wheel circumference C = 2l 2 , there are orders o 1 = 1/2, o 2 = 2. Therefore, the signals z 1 (t) and z 2 (t) essentially represent the physical meaning of the vibration response generated by a wheel traveling on a track with two periods of irregular characteristics at a speed of R(t).The wavelength of these two periods of irregularity is one-half of the wheel circumference, i.e., l 2 = C/2; one is two times the wheel circumference, i.e., l 1 = 2C.The wheel rotational speed is shown in Figure 12.
Sensors 2024, 24, 3804 14 of 30 two periods of irregularity is one-half of the wheel circumference, i.e., l2 = C/2; one is two times the wheel circumference, i.e., l1 = 2C.The wheel rotational speed is shown in Figure 12.

Resample the Time-Domain Signal Non-Uniformly in the Equal-Angle Sample Method
Once the wheel rotation speed time-domain signal is obtained, the original time signal can be resampled non-uniformly in the equal-angle sample method.Before that, we summarize some parameter relationships for order tracking analysis.When the length of the sampled data is an integer multiple P of the arc length of the wheel rotation, and the number of sample points for each revolution of the wheel is M, the maximum analysis order Omax in the order domain is M/2, and the order analysis resolution ΔO is 1/P.In other words, M is the sampling rate and M × P is the sampling length.The larger the value of M, the wider the order domain of the analysis.The larger the value of P, the higher the accuracy of the order domain of the analysis.
Since the order o1 = 1/2, o2 = 2 in this case, according to Nyquist's sampling theorem, the sampling points number for each wheel revolution should be M = 2•max(o1, o2) = 4.Meanwhile, like the traditional frequency domain analysis, the order domain analysis also exists the phenomenon of signal aliasing, caused by the sampling interval Δθ being too large or the sampling rate M being too small, which will cause serious estimation errors.There are two ways to solve the signal aliasing problem.① An anti-aliasing filter-using the low-pass filter to clear the frequency components larger than half of the sampling frequency can avoid signal aliasing problems.However the actual filter hard has the ideal filter characteristics, and it is easy to introduce artificial errors.② Higher sampling rateincreasing the sampling rate or decreasing the sampling interval can solve the signal aliasing problems.In this case, since M = (2.56~4)•max(o1,o2), let M = 8 and Δθ = 2π/8 = π/4.Finally, considering the order analysis resolution ΔO = 1/2 can guarantee that the order spectrum does not suffer from spectrum value leakage.Therefore, P = 2 n .In this case, let P = 2, which means the signal must be collected for two revolutions of the wheel.
After resampling the z(t) signal, the angular domain signal is shown in Figure 13.The horizontal coordinate is the wheel rotation phase angle, Δθ = π/4.The vertical coordinate is the vertical irregularity displacement excitation whose amplitude remains unchanged

Resample the Time-Domain Signal Non-Uniformly in the Equal-Angle Sample Method
Once the wheel rotation speed time-domain signal is obtained, the original time signal can be resampled non-uniformly in the equal-angle sample method.Before that, we summarize some parameter relationships for order tracking analysis.When the length of the sampled data is an integer multiple P of the arc length of the wheel rotation, and the number of sample points for each revolution of the wheel is M, the maximum analysis order O max in the order domain is M/2, and the order analysis resolution ∆O is 1/P.In other words, M is the sampling rate and M × P is the sampling length.The larger the value of M, the wider the order domain of the analysis.The larger the value of P, the higher the accuracy of the order domain of the analysis.
Since the order o 1 = 1/2, o 2 = 2 in this case, according to Nyquist's sampling theorem, the sampling points number for each wheel revolution should be M = 2•max(o 1 , o 2 ) = 4.Meanwhile, like the traditional frequency domain analysis, the order domain analysis also exists the phenomenon of signal aliasing, caused by the sampling interval ∆θ being too large or the sampling rate M being too small, which will cause serious estimation errors.There are two ways to solve the signal aliasing problem. 1  ⃝ An anti-aliasing filter-using the low-pass filter to clear the frequency components larger than half of the sampling frequency can avoid signal aliasing problems.However the actual filter hard has the ideal filter characteristics, and it is easy to introduce artificial errors. 2  ⃝ Higher sampling rate-increasing the sampling rate or decreasing the sampling interval can solve the signal aliasing problems.In this case, since M = (2.56~4)•max(o 1 , o 2 ), let M = 8 and ∆θ = 2π/8 = π/4.Finally, considering the order analysis resolution ∆O = 1/2 can guarantee that the order spectrum does not suffer from spectrum value leakage.Therefore, P = 2 n .In this case, let P = 2, which means the signal must be collected for two revolutions of the wheel.
After resampling the z(t) signal, the angular domain signal is shown in Figure 13.The horizontal coordinate is the wheel rotation phase angle, ∆θ = π/4.The vertical coordinate is the vertical irregularity displacement excitation whose amplitude remains unchanged at 10 µm and 5 µm.Compared with the z(t) signal in Figure 9, it shows good stationarity after converting the time-domain signal to the angular domain signal.

Order Spectrum Analysis and Quantitative Detection
Then, the discrete Fourier transform is performed to obtain the order spectrum, as shown in Figure 14.The horizontal coordinate is the order that reflects the ratio of the wheel circumference C to the track vertical irregularity wavelength, as shown in Equation (32).The vertical coordinate is the amplitude.Compared with Figures 10 and 11, the amplitude of the corresponding vibration can be clearly identified with 0 error, proving the advantage of order tracking analysis to the non-stationary signal whose vibration frequency is related to the wheel rotational speed, as shown in Figures 13 and 14.In practical engineering, since the exact wavelength of the track irregularity is not known at the beginning, it is difficult to ensure that the collected signal corresponds to an integer multiple of the period irregularity, causing the signals to be non-periodic.If the discrete Fourier transform is directly performed on these non-periodic signals, it will cause energy leakage in the order spectrum.For this case, like the traditional frequency domain analysis, the window function (Hanning) can be used to the resampled angular domain signal, and then correcting the spectrum value of the discrete Fourier transformed spectrum by multiplying 2.0, the desired result can be obtained.Based on the order information, it can be converted into a spatial domain signal of track vertical irregularity: where u is the vertical track irregularity, x is the longitudinal distance of the track, l1 and l2 are the track irregularity wavelength, C is the wheel circumference, and o1 and o2 are the orders.
The curves in Figure 13 could be smoother, mainly because the M and P determined

Order Spectrum Analysis and Quantitative Detection
Then, the discrete Fourier transform is performed to obtain the order spectrum, as shown in Figure 14.The horizontal coordinate is the order that reflects the ratio of the wheel circumference C to the track vertical irregularity wavelength, as shown in Equation (32).The vertical coordinate is the amplitude.Compared with Figures 10 and 11, the amplitude of the corresponding vibration can be clearly identified with 0 error, proving the advantage of order tracking analysis to the non-stationary signal whose vibration frequency is related to the wheel rotational speed, as shown in Figures 13 and 14.In practical engineering, since the exact wavelength of the track irregularity is not known at the beginning, it is difficult to ensure that the collected signal corresponds to an integer multiple of the period irregularity, causing the signals to be non-periodic.If the discrete Fourier transform is directly performed on these non-periodic signals, it will cause energy leakage in the order spectrum.For this case, like the traditional frequency domain analysis, the window function (Hanning) can be used to the resampled angular domain signal, and then correcting the spectrum value of the discrete Fourier transformed spectrum by multiplying 2.0, the desired result can be obtained.

Order Spectrum Analysis and Quantitative Detection
Then, the discrete Fourier transform is performed to obtain the order spectrum, as shown in Figure 14.The horizontal coordinate is the order that reflects the ratio of the wheel circumference C to the track vertical irregularity wavelength, as shown in Equation (32).The vertical coordinate is the amplitude.Compared with Figures 10 and 11, the amplitude of the corresponding vibration can be clearly identified with 0 error, proving the advantage of order tracking analysis to the non-stationary signal whose vibration frequency is related to the wheel rotational speed, as shown in Figures 13 and 14.In practical engineering, since the exact wavelength of the track irregularity is not known at the beginning, it is difficult to ensure that the collected signal corresponds to an integer multiple of the period irregularity, causing the signals to be non-periodic.If the discrete Fourier transform is directly performed on these non-periodic signals, it will cause energy leakage in the order spectrum.For this case, like the traditional frequency domain analysis, the window function (Hanning) can be used to the resampled angular domain signal, and then correcting the spectrum value of the discrete Fourier transformed spectrum by multiplying 2.0, the desired result can be obtained.Based on the order information, it can be converted into a spatial domain signal of track vertical irregularity: where u is the vertical track irregularity, x is the longitudinal distance of the track, l1 and l2 are the track irregularity wavelength, C is the wheel circumference, and o1 and o2 are the orders.
The curves in Figure 13 could be smoother, mainly because the M and P determined Based on the order information, it can be converted into a spatial domain signal of track vertical irregularity: where u is the vertical track irregularity, x is the longitudinal distance of the track, l 1 and l 2 are the track irregularity wavelength, C is the wheel circumference, and o 1 and o 2 are the orders.The curves in Figure 13 could be smoother, mainly because the M and P determined in the above cases are smaller enough, which satisfies the analysis accuracy requirement.Compared with Figures 13 and 14, as shown in Figures 15-18, with increasing M, the sampling point number of each wheel revolution will increase, and the sampling interval ∆θ will decrease, resulting in a smoother curve.With the increase of P, the number of sampled rotations will increase, resulting in a longer signal.Therefore, to detect short-wave track irregularity, we should try to make M larger.Increasing the P value will achieve a better order analysis resolution, which will help to distinguish the two types of track irregularity with similar wavelengths in the order spectrum.
Sensors 2024, 24, 3804 16 of 3 rotations will increase, resulting in a longer signal.Therefore, to detect short-wave trac irregularity, we should try to make M larger.Increasing the P value will achieve a bette order analysis resolution, which will help to distinguish the two types of track irregularit with similar wavelengths in the order spectrum.Sensors 2024, 24, 3804 16 of 3 rotations will increase, resulting in a longer signal.Therefore, to detect short-wave track irregularity, we should try to make M larger.Increasing the P value will achieve a bette order analysis resolution, which will help to distinguish the two types of track irregularity with similar wavelengths in the order spectrum.Sensors 2024, 24, 3804 16 of rotations will increase, resulting in a longer signal.Therefore, to detect short-wave trac irregularity, we should try to make M larger.Increasing the P value will achieve a bette order analysis resolution, which will help to distinguish the two types of track irregularit with similar wavelengths in the order spectrum.

Simulation Case Analysis and Verification
To verify the accuracy and reliability of the above methods, this section conducts simulation analysis and verification by establishing a vehicle system dynamics model.

Track-Vehicle Dynamic Model
Due to the weak coupling between the vertical and lateral dynamics of rail vehicles, the modeling can be carried out independently to reduce the complexity of the problem.The vehicle vertical dynamics model is shown in Figure 19.

Simulation Case Analysis and Verification
To verify the accuracy and reliability of the above methods, this section conducts simulation analysis and verification by establishing a vehicle system dynamics model.

Track-Vehicle Dynamic Model
Due to the weak coupling between the vertical and lateral dynamics of rail vehicles, the modeling can be carried out independently to reduce the complexity of the problem.The vehicle vertical dynamics model is shown in Figure 19.where Mb is the body mass, Mt is the bogie frame mass, Jb is the body nodding inertia, Jt is the bogie frame nodding inertia, ks is the secondary spring vertical rigidity, kp is the primary spring vertical rigidity, cs is the secondary spring vertical damping, and cp is the primary spring vertical damping.
As shown in Table 1, a series of track vertical irregularities are designed and are input as the excitation in the above model.The wheel circumference is 2639 mm.The vehicle runs on a 1000 m straight track with an initial speed of 10.8 km/h.The axle-box vertical vibration acceleration and the vehicle running speed are output to estimate the track irregularities.The sampling frequency is 1000 Hz, and the sampling time is 190 s, as shown in Figures 20 and 21.where M b is the body mass, M t is the bogie frame mass, J b is the body nodding inertia, J t is the bogie frame nodding inertia, k s is the secondary spring vertical rigidity, k p is the primary spring vertical rigidity, c s is the secondary spring vertical damping, and c p is the primary spring vertical damping.
As shown in Table 1, a series of track vertical irregularities are designed and are input as the excitation in the above model.The wheel circumference is 2639 mm.The vehicle runs on a 1000 m straight track with an initial speed of 10.8 km/h.The axle-box vertical vibration acceleration and the vehicle running speed are output to estimate the track irregularities.The sampling frequency is 1000 Hz, and the sampling time is 190 s, as shown in Figures 20 and 21.

Signal Analysis and Pre-Processing
The axle-box vertical vibration displacement time-domain signal can be obtained by time-domain double integration and trend term elimination, according to the method in Section 3, as shown in Figure 22.Since the low vehicle speed in this validation scheme, the transfer relationship between the axle-box vertical displacement and the track vertical irregularity displacement excitation is close to 1. Therefore, the axle-box vertical displacement basically follows the track vertical irregularity at this time.However, as shown in Figure 21, the vehicle speed frequently changes during operation.The axle-box vertical displacement signal in Figure 22 is non-stationary.As shown in Figures 10 and 11, traditional signal processing methods cannot quantitatively detect the wavelength and amplitude of track irregularities from non-stationary signals.However, the time signal in Figure 22 can be resampled non-uniformly in the equal-angle sample method and converted into the angular (spatial) domain signal by using the proposed method in Section 4. The quantitative detection of track irregularities can be achieved.According to the analysis method in Section 4, it is known that the estimation of the wheel rotational speed R is the key to converting the time-domain signal into the angular domain signal.The relationship between the wheel rotational speed R and the vehicle running speed v is shown in Equation ( 32), and the wheel rotational speed is calculated, as shown in Figure 23.

Signal Analysis and Pre-Processing
The axle-box vertical vibration displacement time-domain signal can be obtained by time-domain double integration and trend term elimination, according to the method in Section 3, as shown in Figure 22.Since the low vehicle speed in this validation scheme, the transfer relationship between the axle-box vertical displacement and the track vertical irregularity displacement excitation is close to 1. Therefore, the axle-box vertical displacement basically follows the track vertical irregularity at this time.

Signal Analysis and Pre-Processing
The axle-box vertical vibration displacement time-domain signal can be obtained by time-domain double integration and trend term elimination, according to the method in Section 3, as shown in Figure 22.Since the low vehicle speed in this validation scheme, the transfer relationship between the axle-box vertical displacement and the track vertical irregularity displacement excitation is close to 1. Therefore, the axle-box vertical displacement basically follows the track vertical irregularity at this time.However, as shown in Figure 21, the vehicle speed frequently changes during operation.The axle-box vertical displacement signal in Figure 22 is non-stationary.As shown in Figures 10 and 11, traditional signal processing methods cannot quantitatively detect the wavelength and amplitude of track irregularities from non-stationary signals.However, the time signal in Figure 22 can be resampled non-uniformly in the equal-angle sample method and converted into the angular (spatial) domain signal by using the proposed method in Section 4. The quantitative detection of track irregularities can be achieved.According to the analysis method in Section 4, it is known that the estimation of the wheel rotational speed R is the key to converting the time-domain signal into the angular domain signal.The relationship between the wheel rotational speed R and the vehicle running speed v is shown in Equation ( 32), and the wheel rotational speed is calculated, as shown in Figure 23.However, as shown in Figure 21, the vehicle speed frequently changes during operation.The axle-box vertical displacement signal in Figure 22 is non-stationary.As shown in Figures 10 and 11, traditional signal processing methods cannot quantitatively detect the wavelength and amplitude of track irregularities from non-stationary signals.However, the time signal in Figure 22 can be resampled non-uniformly in the equal-angle sample method and converted into the angular (spatial) domain signal by using the proposed method in Section 4. The quantitative detection of track irregularities can be achieved.According to the analysis method in Section 4, it is known that the estimation of the wheel rotational speed R is the key to converting the time-domain signal into the angular domain signal.
The relationship between the wheel rotational speed R and the vehicle running speed v is shown in Equation (32), and the wheel rotational speed is calculated, as shown in Figure 23.
the wavelength and amplitude of track irregularities from non-stationary signals.However, the time signal in Figure 22 can be resampled non-uniformly in the equal-angle sample method and converted into the angular (spatial) domain signal by using the proposed method in Section 4. The quantitative detection of track irregularities can be achieved.According to the analysis method in 4, it is known that the estimation of the wheel rotational speed R is the key to converting the time-domain signal into the angular domain signal.The relationship between the wheel rotational speed R and the vehicle running speed v is shown in Equation (32), and the wheel rotational speed is calculated, as shown in Figure 23.Based on the improved Simpson integration method and trend item decomposition method in Section 3, the wheel rotation phase angle time-domain signal can be calculated as shown in Equation ( 36).Then, determine the time t n (n = 0, 1, 2, . ..) according to the wheel rotation phase angle (θ n = n • ∆θ, n = 0, 1, 2, . ..).Finally, we sequentially extract the corresponding values from the axle-box vertical displacement time-domain signal (Figure 22) based on the time t n (n = 0, 1, 2, . ..) to obtain the resampled signal.
where θ n is the wheel rotation phase angle time series, and R(n) is the wheel rotational speed time series.
In the above case study, since the wavelength characteristics of track irregularity, as shown in Table 1, combined with the wheel circumference and vehicle running speed, it is necessary for the order spectrum analysis to accurately express the characteristics of track irregularity when ∆O ≤ 0.01 and M ≥ 4 × 52.779 (i.e., P ≥ 50, M ≥ 211.11).Therefore, let P = 100 and M = 240 here to ensure sufficient accuracy.That is, collect the vibration signal corresponding to 100 revolutions of the wheel with resampled at equal angular intervals for 240 points per revolution (i.e., ∆θ = 2π/M = π/120).However, all the above signals are discrete (non-continuous) and cannot be expressed by specific functions.When resampling by equal angles, the time t n (n = 0, 1, 2, . ..) may not exist or between two discrete time points, causing the sampling failure.Therefore, this paper proposes a nonlinear interpolation method to obtain an accurate value between two discrete points as follows.
Assuming a cubic piecewise polynomial function S(t) that passes through t 1 , t 2 , t 3 , and t 4 , to calculate the S(t n ) (t 2 < t n < t 3 ), the analytical expression of the function S(t) needs to be solved.Therefore, 12 conditional formulas must be needed to determine the 12 polynomial coefficients in total: 1  ⃝ the S(t 1 ) = B 1 , S(t 2 ) = B 2 , S(t 3 ) = B 3 , and S(t 4 ) = B 4 are already known by simulated output signal; 2  ⃝ based on the fact that S(t) has a continuous second derivative and satisfies continuity conditions at t 2 and t 3 , i.e., S(t j−0 ) = S(t j+0 ), .S(t j−0 ) = .S(t j+0 ), .. S(t j−0 ) = .. S(t j+0 ), j = 2, 3; 3  ⃝ there are two endpoint boundary conditions are already known by simulated output signal, i.e., .
We solve the cubic piecewise polynomial S(t) based on the above 12 conditional formulas, and substitute t = t n into it to obtain S(t n ).The resampled angular domain signal of the axle-box vertical vibration displacement is shown in Figure 24, which is the spatial domain signal and does not have non-stationary characteristics.
; ③ there are two endpoint boundary conditions are already known by simulated output signal, i.e., . We solve the cubic piecewise polynomial S(t) based on the above 12 conditional formulas, and substitute t = tn into it to obtain S(tn).The resampled angular domain signal of the axle-box vertical vibration displacement is shown in Figure 24, which is the spatial domain signal and does not have non-stationary characteristics.

Order Spectrum Analysis
We perform the discrete Fourier transform to the angle-domain signal in Figure 24 to obtain the order spectrum.As shown in Figure 25, the horizontal axle represents the order, which is expressed as the ratio of the wheel circumference to the track irregularity wavelength.The vertical axle represents the amplitude of the track irregularity.Compared with the real amplitude in Table 1, the wavelength estimation error of track irregularity is less than 2.31%, and the amplitude estimation error of the track irregularity is less than 20%, as shown in Table 2.

Order Spectrum Analysis
We perform the discrete Fourier transform to the angle-dom obtain the order spectrum.As shown in Figure 25, the horizontal a which is expressed as the ratio of the wheel circumference to the length.The vertical axle represents the amplitude of the track irre the real amplitude in Table 1, the wavelength estimation error of than 2.31%, and the amplitude estimation error of the track irreg as shown in Table 2.As shown in Table 2, the accuracy of the order analysis in estimating the short and medium wave irregularity is significantly higher than the long wave irregularity.The main reasons are as follows: 1  ⃝ since the long wave irregularity generally corresponds to low frequency signals, which are easily susceptible to interference from trend terms.Because the trend term can not be completely eliminated, even a small residual error will affect the estimation accuracy; 2  ⃝ due to the wide variety of track irregularity wavelengths, it is impossible to ensure the intercepted signal is an integral multiple of the period of corresponding wavelength irregularity.There must be obvious truncation errors when analyzing the track irregularity.Overall, compared to the traditional frequency domain analysis method, the order analysis method for identifying track irregularities under nonstationary conditions has significant advantages.Section 6 will further analyze and discuss its application in practical engineering.

Field Test Case Analysis and Verification
Based on the analysis and derivation in Sections 3 and 4, using the axle-box vertical vibration acceleration and vehicle running speed to estimate the track vertical irregularity state (wavelength, amplitude) is feasible.This section will conduct analysis and research to apply the above series of methods in practical engineering.Since the measured signals often have severe test errors and outliers, it is to perform the required pre-processing to the test signals, such as zero drift el singular value elimination, electronic noise elimination, etc.In this paper, ba signal processing toolbox that comes with MATLAB R2023b, the Savitzky G smoothing function is called to pre-process the initial signal to eliminate the zero trend terms due to real-time testing.As shown in Figure 27, the axle-box vibra eration signal collected from the down-process shows severe zero-frequency st rors in the local time domain, and the difference before and after signal proces nificant, showing the importance of signal pre-processing.Since the measured signals often have severe test errors and outliers, it is necessary to perform the required pre-processing to the test signals, such as zero drift elimination, singular value elimination, electronic noise elimination, etc.In this paper, based on the signal processing toolbox that comes with MATLAB R2023b, the Savitzky Golay filter smoothing function is called to pre-process the initial signal to eliminate the zero-drift and trend terms due to real-time testing.As shown in Figure 27, the axle-box vibration acceleration signal collected from the down-process shows severe zero-frequency stepwise errors in the local time domain, and the difference before and after signal processing is significant, showing the importance of signal pre-processing.

Calculation of Time-Domain Track Irregularity
smoothing function is called to pre-process the initial signal to eliminate the zerotrend terms due to real-time testing.As shown in Figure 27, the axle-box vibratio eration signal collected from the down-process shows severe zero-frequency step rors in the local time domain, and the difference before and after signal processin nificant, showing the importance of signal pre-processing.According to the proposed method in Section 3, we performed the time-dom ble integration and trend term elimination to calculate the axle-box vibration d ment.Then, according to the transfer relationship H(jω) between the axle-box ver placement and the track irregularity excitation in Equation ( 8), the track irregula placement excitation u(t) can be obtained, as shown in Figure 28.According to the proposed method in Section 3, we performed the time-domain double integration and trend term elimination to calculate the axle-box vibration displacement.Then, according to the transfer relationship H(jω) between the axle-box vertical displacement and the track irregularity excitation in Equation ( 8), the track irregularity displacement excitation u(t) can be obtained, as shown in Figure 28.As shown in Figure 28, there is no apparent short-wave irregularity, which mainly shows height variation throughout the trajectory.This variation can be seen as a low-frequency trend term.The trend and fluctuation terms can be obtained separately by applying the Hodrick-Prescott decomposition method in Section 3. Since the trend term is very close to Figure 28, it is omitted here.The fluctuation term is shown in Figure 29.Compared the axle-box acceleration signal in Figure 27 with the track irregularity displacement excitation in Figure 29, it is found that these two are similar in amplitude and frequency characteristics, with a phase difference of π, proving the correctness of the above calculation method.

Resampling of Time-Domain Displacement Excitation Signal
Due to the vehicle running at variable speed, the time-domain signal in Figure 29 is non-stationary.However, the traditional signal processing method cannot quantitatively detect the wavelength and amplitude of track irregularities from it.According to the method proposed in Section 4, it is necessary to further resample the time-domain dis- As shown in Figure 28, there is no apparent short-wave irregularity, which mainly shows height variation throughout the trajectory.This variation can be seen as a lowtrend term.The trend and fluctuation terms can be obtained separately by applying the Hodrick-Prescott decomposition method in Section 3. Since the trend term is very close to Figure 28, it is omitted here.The fluctuation term is shown in Figure 29.Compared the axle-box acceleration signal in Figure 27 with the track irregularity displacement excitation in Figure 29, it is found that these two are similar in amplitude and frequency characteristics, with a phase difference of π, proving the correctness of the above calculation method.As shown in Figure 28, there is no apparent short-wave irregularity, which mainly shows height variation throughout the trajectory.This variation can be seen as a low-frequency trend term.The trend and fluctuation terms can be obtained separately by applying the Hodrick-Prescott decomposition method in Section 3. Since the trend term is very close to Figure 28, it is omitted here.The fluctuation term is shown in Figure 29.Compared the axle-box acceleration signal in Figure 27 with the track irregularity displacement excitation in Figure 29, it is found that these two are similar in amplitude and frequency characteristics, with a phase difference of π, proving the correctness of the above calculation method.

Resampling of Time-Domain Displacement Excitation Signal
Due to the vehicle running at variable speed, the time-domain signal in Figure 29 is non-stationary.However, the traditional signal processing method cannot quantitatively detect the wavelength and amplitude of track irregularities from it.According to the method proposed in Section 4, it is necessary to further resample the time-domain dis-

Resampling of Time-Domain Displacement Excitation Signal
Due to the vehicle running at variable speed, the time-domain signal in Figure 29 is non-stationary.However, the traditional signal processing method cannot quantitatively detect the wavelength and amplitude of track irregularities from it.According to the method proposed in Section 4, it is necessary to further resample the time-domain displacement excitation non-uniformly in the equal-angle sample method based on the wheel angular velocity signal and finally achieve quantitative detection of track irregularity through order spectrum analysis.Hence, estimating the wheel angular velocity (rotational speed) signal is critical for equal angle resampling and subsequent order analysis.There are two methods proposed for obtaining wheel angular velocity signals: Based on vehicle running speed signal, rail vehicles usually have their own speed sensors, which can directly derive the real-time vehicle running speed signal from the control system.Then, the wheel angular velocity signals can be converted according to Equation (32), as shown in Figure 30.However, the vehicle running speed derived from the control system is essentially the overall operating speed of the vehicle.Due to the non-rigid connection between vehicles, the overall speed of the vehicle may not represent the real-time angular velocity of a certain wheel, and minor differences may also lead to significant estimation errors.Based on key-phase signal, to actually obtain a more realistic real-time wheel rotation angular velocity, the measurement of key-phase signals can be carried out.The keyphase signal is a continuous pulse signal that occurs quantitatively with each wheel revolution and is measured by using the direct beam photoelectric sensor, which is attached to the bottom of the axle-box and aligns the direction of the light source with the wheel as shown in Figure 31.We then attach a reflective plate to the outside of the wheel and align it with the light source emitted by the photoelectric sensor.Finally, we connect the data line of the photoelectric sensor to the signal acquisition box in Figure 26b to collect the continuous photoelectric pulse signal generated during wheel operation.The sampling frequency is essential for the collection of key phase signals.In this case study, the typical operating speed of the subway vehicle is 65 km/h, and the wheel diameter is 0.74 m.The time for the wheels to rotate once is 0.15 s.To ensure the photoelectric sensor can clearly detect the key-phase signals, a higher sampling frequency must be satisfied, as shown in Equation (37).The sampling frequency should be at least 3420 Hz.In this case study, the sampling frequency of the key-phase signal is set to 5120 Hz, which is consistent with the vibration acceleration sampling frequency.It is convenient for later signal processing.The collected key-phase signal is shown in Figure 32.Based on key-phase signal, to actually obtain a more realistic real-time wheel rotation angular velocity, the measurement of key-phase signals can be carried out.The key-phase signal is a continuous pulse signal that occurs quantitatively with each wheel revolution and is measured by using the direct beam photoelectric sensor, which is attached to the bottom of the axle-box and aligns the direction of the light source with the wheel as shown in Figure 31.We then attach a reflective plate to the outside of the wheel and align it with the light source emitted by the photoelectric sensor.Finally, we connect the data line of the photoelectric sensor to the signal acquisition box in Figure 26b to collect the continuous photoelectric pulse signal generated during wheel operation.Based on key-phase signal, to actually obtain a more realistic real-time wheel rotation angular velocity, the measurement of key-phase signals can be carried out.The keyphase signal is a continuous pulse signal that occurs quantitatively with each wheel revolution and is measured by using the direct beam photoelectric sensor, which is attached to the bottom of the axle-box and aligns the direction of the light source with the wheel as shown in Figure 31.We then attach a reflective plate to the outside of the wheel and align it with the light source emitted by the photoelectric sensor.Finally, we connect the data line of the photoelectric sensor to the signal acquisition box in Figure 26b to collect the continuous photoelectric pulse signal generated during wheel operation.The sampling frequency is essential for the collection of key phase signals.In this case study, the typical operating speed of the subway vehicle is 65 km/h, and the wheel diameter is 0.74 m.The time for the wheels to rotate once is 0.15 s.To ensure the photoelectric sensor can clearly detect the key-phase signals, a higher sampling frequency must be satisfied, as shown in Equation (37).The sampling frequency should be at least 3420 Hz.In this case study, the sampling frequency of the key-phase signal is set to 5120 Hz, which is consistent with the vibration acceleration sampling frequency.It is convenient for later signal processing.The collected key-phase signal is shown in Figure 32.The sampling frequency is essential for the collection of key phase signals.In this case study, the typical operating speed of the subway vehicle is 65 km/h, and the wheel diameter is 0.74 m.The time for the wheels to rotate once is 0.15 s.To ensure the photoelectric sensor can clearly detect the key-phase signals, a higher sampling frequency must be satisfied, as shown in Equation (37).The sampling frequency should be at least 3420 Hz.In this case study, the sampling frequency of the key-phase signal is set to 5120 Hz, which is consistent with the vibration acceleration sampling frequency.It is convenient for later signal processing.The collected key-phase signal is shown in Figure 32.
where f s is the sampling frequency, T min is the minimum rotation period of the wheel, and N is the minimum number of sampling points for one cycle of wheel rotation.When N ≥ 500, it can ensure that the signal is almost undistorted.The key-phase signal is characterized as a rectangular waveform.The interval between two adjacent rising edges (or falling edges) represents one wheel rotation period, as T1 and T2, as shown in Figure 32.The pulse widths δ1 and δ2 represent the movement time length of the light point emitted by the photoelectric sensor on the reflective white plate.It can be found that the above two adjacent rotation periods T and pulse width δ are not equal, indicating that the rotation speed of this wheel is non-uniform.Therefore, the axle-box acceleration collected during this period is typically non-stationary.
Then, using the time point of the rising edge, the wheel rotation speed can be calculated as follows: where t1_r and t2_r represents the time point of two adjacent rising edges, T1 = t2_r-t1_r, k is the number of pulses per wheel revolution, k = 1; and r1_r is the wheel rotation speed at time t1_r.
As shown in Equation (38), this calculation method takes the average speed within the time interval [t1_r, t2_r] as the instantaneous velocity at time t1_r.Similarly, the time point of the falling edge also can be used as follows: ( ) Finally, we obtain the instantaneous wheel rotation speed at the time point (t1_r + t1_f)/2 by calculating the average value of r1_r and r1_f.According to the conversion relationship between wheel rotation speed and wheel angular velocity in Equation (32), the real-time wheel rotation angular velocity signal calculated, as shown in Figure 33.The key-phase signal is characterized as a rectangular waveform.The interval between two adjacent rising edges (or falling edges) represents one wheel rotation period, as T 1 and T 2 , as shown in Figure 32.The pulse widths δ 1 and δ 2 represent the movement time length of the light point emitted by the photoelectric sensor on the reflective white plate.It can be found that the above two adjacent rotation periods T and pulse width δ are not equal, indicating that the rotation speed of this wheel is non-uniform.Therefore, the axle-box acceleration collected during this period is typically non-stationary.
Then, using the time point of the rising edge, the wheel rotation speed can be calculated as follows: where t 1_r and t 2_r represents the time point of two adjacent rising edges, T 1 = t 2_r − t 1_r , k is the number of pulses per wheel revolution, k = 1; and r 1_r is the wheel rotation speed at time t 1_r .
As shown in Equation (38), this calculation method takes the average speed within the time interval [t 1_r , t 2_r ] as the instantaneous velocity at time t 1_r .Similarly, the time point of the falling edge also can be used as follows: Finally, we obtain the instantaneous wheel rotation speed at the time point (t 1_r + t 1_f )/2 by calculating the average value of r 1_r and r 1_f .According to the conversion relationship between wheel rotation speed and wheel angular velocity in Equation (32), the real-time wheel rotation angular velocity signal calculated, as shown in Figure 33.
( ) Finally, we obtain the instantaneous wheel rotation speed at the time point (t1_r + t1_f)/2 by calculating the average value of r1_r and r1_f.According to the conversion relationship between wheel rotation speed and wheel angular velocity in Equation (32), the real-time wheel rotation angular velocity signal calculated, as shown in Figure 33.As shown in Figure 33, the calculated real-time wheel angular velocity based on the photoelectric key-phase signal is almost identical to that based on vehicle running speed with only some minor fluctuations, indicating the calculation method is theoretically wholly correct.
Then, we perform the time-domain integration to the wheel rotation angular velocity signal (Figure 33) to calculate the phase angle time-domain signal, as shown in Figure 34.The accumulated phase angle of the wheel rotation gradually increases with time until the vehicle finally stops.Meanwhile, the phase angle curve does not increase uniformly, which is related to the frequent changes in the vehicle running speed.As shown in Figure 33, the calculated real-time wheel angular velocity based on the photoelectric key-phase signal is almost identical to that based on vehicle running speed with only some minor fluctuations, indicating the calculation method is theoretically wholly correct.
Then, we perform the time-domain integration to the wheel rotation angular velocity signal (Figure 33) to calculate the phase angle time-domain signal, as shown in Figure 34.The accumulated phase angle of the wheel rotation gradually increases with time until the vehicle finally stops.Meanwhile, the phase angle curve does not increase uniformly, which is related to the frequent changes in the vehicle running speed.According to the relationship between the order o, wheel circumference C, and track irregularity wavelength l (i.e., o = C/l), when the minimum wavelength of concern was 0.05 m, the highest order should be omax = C/lmin = 46.5.As shown in Section 4.2, the order spectrum analysis must at least satisfy ΔO ≤ 0.01, M ≥ 2 × 46.5, i.e., P ≥ 50, M ≥ 93.Meanwhile, a larger M must be considered to avoid the signal aliasing phenomenon and let M = 240 ≥ 5 × 46.5 here.Then, P is expressed as the number of wheel revolutions.Since the resampling process is necessarily performed on the entire sampled signal, the value of P essentially depends on the length of the original sampled signal.The longer the sampling time, the greater the number of wheel revolutions, the greater the value of P, which can be expressed as: where   x represents the largest integer less than x, and θ is the phase angle, as the maxi- mum value in Figure 34.
In this case study, the phase angle of the wheel rotation until the vehicle finally stops is 57,931.05rad.The P≈9220 indicates that the total wheel rotation during the time course in Figure 33 is approximately 9220 revolutions.Hence, the resolution of the order spectrum is ΔO = 1/P = 0.00010846 ≤ 0.01, which can satisfy the above estimation accuracy requirements.According to the resampling setting conditions above, the time-domain track irregularity displacement excitation in Figure 28 is resampled in the equal-angle sample method (Δθ = 2π/240).As shown in Figure 35, the equal-angle sample result is expressed as non-uniform sampling in the time domain.When the vehicle speed is low, the sampling points are sparse; otherwise, the sampling points are more intensive and can be expressed in the form of an angular domain signal, as shown in Figure 36.The horizontal coordinate is the phase angle of the wheel, and the vertical coordinate is the vertical According to the relationship between the order o, wheel circumference C, and track irregularity wavelength l (i.e., o = C/l), when the minimum wavelength of concern was 0.05 m, the highest order should be o max = C/l min = 46.5.As shown in Section 4.2, the order spectrum analysis must at least satisfy ∆O ≤ 0.01, M ≥ 2 × 46.5, i.e., P ≥ 50, M ≥ 93.Meanwhile, a larger M must be considered to avoid the signal aliasing phenomenon and let M = 240 ≥ 5 × 46.5 here.Then, P is expressed as the number of wheel revolutions.Since the resampling process is necessarily performed on the entire sampled signal, the value of P essentially depends on the length of the original sampled signal.The longer the sampling time, the greater the number of wheel revolutions, the greater the value of P, which can be expressed as: where ⌊x⌋ represents the largest integer less than x, and θ is the phase angle, as the maximum value in Figure 34.
In this case study, the phase angle of the wheel rotation until the vehicle finally stops is 57,931.05rad.The P ≈ 9220 indicates that the total wheel rotation during the time course in Figure 33 is approximately 9220 revolutions.Hence, the resolution of the order spectrum is ∆O = 1/P = 0.00010846 ≤ 0.01, which can satisfy the above estimation accuracy requirements.According to the resampling setting conditions above, the time-domain track irregularity displacement excitation in Figure 28 is resampled in the equal-angle sample method (∆θ = 2π/240).As shown in Figure 35, the equal-angle sample result is expressed as non-uniform sampling in the time domain.When the vehicle speed is low, the sampling points are sparse; otherwise, the sampling points are more intensive and can be expressed in the form of an angular domain signal, as shown in Figure 36.The horizontal coordinate is the phase angle of the wheel, and the vertical coordinate is the vertical irregularity amplitude of the track.The fundamental difference between the two types of expression forms is that the track irregularity in the time domain is a non-stationary signal, while in the spatial domain, it is a stationary signal.

Order Spectrum Analysis and Quantitative Detection
The track vertical irregularity in Figure 36 is a stationary signal in the spatial domain, whose wavelength and amplitude can be detected by discrete Fourier transform.But before that, to avoid the spectrum leakage problem caused by truncation error, refer to the traditional spectrum analysis method, add the Hanning window to it, and then perform spectrum value compensation correction (multiplied by 2.0) for the spectrum result after the discrete Fourier transform.The calculated order analysis spectrum is shown in Figure 37.
Based on the relationship between order, track vertical irregularity wavelength, and wheel circumference, these irregularities with a wavelength less than 116.25 mm correspond to the order over 20th in Figure 37 and have minimal amplitude, indicating the track has a good condition in short-wavelength irregularity wear.The lower order region (less than 1) represents medium and long wave irregularities, whose wavelengths are greater than the wheel's circumference, with relatively stable overall and no obvious periodic irregularity feature.Hence, this subway track line in service for about 2 years has good quality and is mainly characterized by random irregularity.

Order Spectrum Analysis and Quantitative Detection
The track vertical irregularity in Figure 36 is a stationary signal in the spatial domain, whose wavelength and amplitude can be detected by discrete Fourier transform.But before that, to avoid the spectrum leakage problem caused by truncation error, refer to the traditional spectrum analysis method, add the Hanning window to it, and then perform spectrum value compensation correction (multiplied by 2.0) for the spectrum result after the discrete Fourier transform.The calculated order analysis spectrum is shown in Figure 37.
Based on the relationship between order, track vertical irregularity wavelength, and wheel circumference, these irregularities with a wavelength less than 116.25 mm correspond to the order over 20th in Figure 37 and have minimal amplitude, indicating the track has a good condition in short-wavelength irregularity wear.The lower order region (less than 1) represents medium and long wave irregularities, whose wavelengths are greater than the wheel's circumference, with relatively stable overall and no obvious periodic irregularity feature.Hence, this subway track line in service for about 2 years has good quality and is mainly characterized by random irregularity.

Order Spectrum Analysis and Quantitative Detection
The track vertical irregularity in Figure 36 is a stationary signal in the spatial domain, whose wavelength and amplitude can be detected by discrete Fourier transform.But before that, to avoid the spectrum leakage problem caused by truncation error, refer to the traditional spectrum analysis method, add the Hanning window to it, and then perform spectrum value compensation correction (multiplied by 2.0) for the spectrum result after the discrete Fourier transform.The calculated order analysis spectrum is shown in Figure 37.

Order Spectrum Analysis and Quantitative Detection
The track vertical irregularity in Figure 36 is a stationary signal in the spatial domain, whose wavelength and amplitude can be detected by discrete Fourier transform.But before that, to avoid the spectrum leakage problem caused by truncation error, refer to the traditional spectrum analysis method, add the Hanning window to it, and then perform spectrum value compensation correction (multiplied by 2.0) for the spectrum result after the discrete Fourier transform.The calculated order analysis spectrum is shown in Figure 37.
Based on the relationship between order, track vertical irregularity wavelength, and wheel circumference, these irregularities with a wavelength less than 116.25 mm correspond to the order over 20th in Figure 37 and have minimal amplitude, indicating the track has a good condition in short-wavelength irregularity wear.The lower order region (less than 1) represents medium and long wave irregularities, whose wavelengths are greater than the wheel's circumference, with relatively stable overall and no obvious periodic irregularity feature.Hence, this subway track line in service for about 2 years has good quality and is mainly characterized by random irregularity.Based on the relationship between order, track vertical irregularity wavelength, and wheel circumference, these irregularities with a wavelength less than 116.25 mm correspond to the order over 20th in Figure 37 and have minimal amplitude, indicating the track has a good condition in short-wavelength irregularity wear.The lower order region (less than 1) represents medium and long wave irregularities, whose wavelengths are greater than the wheel's circumference, with relatively stable overall and no obvious periodic irregularity feature.Hence, this subway track line in service for about 2 years has good quality and is mainly characterized by random irregularity.
It is worth noting that, as shown in Figure 37, the amplitude of irregularities over the 20th order is minimal and has almost no effect on the overall evaluation of track irregularities.It is unnecessary to take M = 240 when resampling the time-domain displacement excitation.In the following analysis cases, compared with when M = 100, the number of sample points for one wheel revolution is reduced from 240 to 100, significantly improving the resampling efficiency and reducing the computation time from about 32 h to about 16 h.
To verify the accuracy of the quantitative detection method, another set of test data (upprocess) was used to calculate the track vertical irregularity of the above subway line again.Due to the difference in operating speed, the axle-box vibration acceleration and key-phase signal generated by the same vehicle running on the same track will also be completely different from the above signals (Figures 27 and 32).If the estimated track irregularity results based on two completely different sets of test data do not differ significantly, it can prove the accuracy and reliability of the proposed method.As shown in Figure 38, the similarity between the estimation results of track irregularity based on the two test results is extremely high, especially the peaks at orders 6, 7, 8, 9, 10, 11, 12, 13, and 14 are almost identical, which fully proves the accuracy and reliability of the proposed method for identifying vertical periodic track irregularity in this paper.It is worth noting that, as shown in Figure 37, the amplitude of irregularities over the 20th order is minimal and has almost no effect on the overall evaluation of track irregularities.It is unnecessary to take M = 240 when resampling the time-domain displacement excitation.In the following analysis cases, compared with when M = 100, the number of sample points for one wheel revolution is reduced from 240 to 100, significantly improving the resampling efficiency and reducing the computation time from about 32 h to about 16 h.
To verify the accuracy of the quantitative detection method, another set of test data (up-process) was used to calculate the track vertical irregularity of the above subway line again.Due to the difference in operating speed, the axle-box vibration acceleration and key-phase signal generated by the same vehicle running on the same track will also be completely different from the above signals (Figures 27 and 32).If the estimated track irregularity results based on two completely different sets of test data do not differ significantly, it can prove the accuracy and reliability of the proposed method.As shown in Figure 38, the similarity between the estimation results of track irregularity based on the two test results is extremely high, especially the peaks at orders 6, 7, 8, 9, 10, 11, 12, 13, and 14 are almost identical, which fully proves the accuracy and reliability of the proposed method for identifying vertical periodic track irregularity in this paper.The reasons for the detection error in Figure 38 are as follows.① During the testing process in this case study, the axle-box vibration signal in Figure 27 and the photoelectric key-phase signal in Figure 32 are collected separately, not absolutely synchronously.These signals are mainly processed manually to align the time starting points of the acceleration and key phase signals, whose errors present in these processes will affect the accuracy of the subsequent resampling.② The quality of the vibration data collected in this field test is not well.As shown in Figure 27, there was severe step-type zero drift.The introduction of human interference through the moving average processing of the signal will lead to signal distortion.③ There are some differences in the testing interval.The upward process includes a transition section for exiting the station, which is not included in the downward process.Due to the collection without stopping throughout the process, it is difficult to remove this section accurately.However, overall, the detection results based on the two sets of test data have minimal differences, proving the proposed method's reliability and accuracy.

Conclusions
Quantitative detection of track vertical irregularity directly affects the quality and safety of vehicle operation.In this paper, starting from the establishment of a simplified wheel-rail dynamics model, the quantitative relationship between track vertical irregularity excitation and axle-box vertical vibration is theoretically derived in the whole pro- The reasons for the detection error in Figure 38 are as follows. 1  ⃝ During the testing process in this case study, the axle-box vibration signal in Figure 27 and the photoelectric key-phase signal in Figure 32 are collected separately, not absolutely synchronously.These signals are mainly processed manually to align the time starting points of the acceleration and key phase signals, whose errors present in these processes will affect the accuracy of the subsequent resampling. 2  ⃝ The quality of the vibration data collected in this field test is not well.As shown in Figure 27, there was severe step-type zero drift.The introduction of human interference through the moving average processing of the signal will lead to signal distortion. 3 ⃝ There are some differences in the testing interval.The upward process includes a transition section for exiting the station, which is not included in the downward process.Due to the collection without stopping throughout the process, it is difficult to remove this section accurately.However, overall, the detection results based on the two sets of test data have minimal differences, proving the proposed method's reliability and accuracy.

Conclusions
Quantitative detection of track vertical irregularity directly affects the quality and safety of vehicle operation.In this paper, starting from the establishment of a simplified wheel-rail dynamics model, the quantitative relationship between track vertical irregularity excitation and axle-box vertical vibration is theoretically derived in the whole process, and the systematic estimation method and process of track vertical irregularity based on axlebox vertical vibration acceleration is proposed, which mainly involves the following aspects.⃝ The double time-domain integration of axle-box vertical vibration acceleration by using Simpson's time-domain integration method is proposed to obtain the axle-box vertical vibration displacement signal, and the calculation error of the method is optimized and improved. 2 ⃝ A method based on the Hodrick-Prescott decomposition is proposed to eliminate the trend term generated in the time-domain integration process.Meanwhile, the optimal value of the smoothing factor λ is derived theoretically, and the determination process and the suggested value are proposed; 3  ⃝ The method and principle of order ratio analysis are derived, analyzed and verified, and finally, the method of estimating the track vertical irregularity based on the axle-box non-stationary vibration signal under the vehicle variable speed operation condition is proposed, which solves the problems of signal non-stationarity, time-domain and space domain conversion, equal angle resampling, and speed estimation based on the photoelectric key signal.
But it is worth pointing out that: 1 ⃝ the effectiveness of the method proposed in this paper relies on the accurate acquisition of the frequency response function H(ω) in practice.For the different track types and vehicle system parameters, the frequency response function curves are quite different. 2 ⃝ There is a P2 resonance effect between the vehicle and the track system, indicating it is better to avoid the first wheel-rail resonance peak with low passing frequency to estimate the track vertical irregularity, especially when |H(ω)| ≈ 1.

3
⃝ The lower the operating speed of the vehicle, the wider the wavelength range of track irregularity detection.However, a speed that is too low reduces the detection efficiency.This is a problem that needs to be considered as a compromise. 4 ⃝ The effect of vehicle lateral motion on the detection results of vertical track irregularities is simplified without considering the weak coupling effect of vehicle lateral motion and vertical motion.It has to be further optimized for practical applications.

Figure 1 .
Figure 1.The overall technical route.

Figure 1 .
Figure 1.The overall technical route.

Figure 4 .
Figure 4. Searching for the optimal smoothing factor.

Figure 4 .
Figure 4. Searching for the optimal smoothing factor.

Figure 8 .
Figure 8. Compare the calculation results of interpolation and non-interpolation before integration.(a) Speed signal; (b) Displacement signal.

Figure 8 .
Figure 8. Compare the calculation results of interpolation and non-interpolation before integration.(a) Speed signal; (b) Displacement signal.

Figure 19 .
Figure 19.Track-vehicle vertical dynamic model.whereMb is the body mass, Mt is the bogie frame mass, Jb is the body nodding inertia, Jt is the bogie frame nodding inertia, ks is the secondary spring vertical rigidity, kp is the primary spring vertical rigidity, cs is the secondary spring vertical damping, and cp is the primary spring vertical damping.

Figure 19 .
Figure 19.Track-vehicle vertical dynamic model.whereM b is the body mass, M t is the bogie frame mass, J b is the body nodding inertia, J t is the bogie frame nodding inertia, k s is the secondary spring vertical rigidity, k p is the primary spring vertical rigidity, c s is the secondary spring vertical damping, and c p is the primary spring vertical damping.

Figure 25 .
Figure 25.Order analysis spectrum of track irregularity.

Figure 25 .
Figure 25.Order analysis spectrum of track irregularity.
Taking a Chinese subway vehicle as an example, we conducted the axle-box vibration acceleration tests to detect the track irregularity.As shown in Figure 26, glue the acceleration sensor (type: J13510, Shanghai Beizhi Electronic Technology Co., Ltd.(Shanghai, China)) to the axle-box.When the vehicle operates under AW0 (empty load) condition and ATO (automatic train operation) mode, the axle-box vibration acceleration signal is transmitted to the signal acquisition cabinet (type: INV3020D, China Orient Institute of Noise & Vibration, Beijing, China) via the data line.The sampling frequency is 5120 Hz, and the sampling duration is for full process collection.Sensors 2024, 24, 3804 6.1.Calculation of Time-Domain Track Irregularity Taking a Chinese subway vehicle as an example, we conducted the axle-box acceleration tests to detect the track irregularity.As shown in Figure 26, glue t ation sensor (type: J13510, Shanghai Beizhi Electronic Technology Co., Ltd.China)) to the axle-box.When the vehicle operates under AW0 (empty load) con ATO (automatic train operation) mode, the axle-box vibration acceleration sign mitted to the signal acquisition cabinet (type: INV3020D, China Orient Institut & Vibration, Beijing, China) via the data line.The sampling frequency is 5120 H sampling duration is for full process collection.

Figure 29 .
Figure 29.The fluctuation term of track irregularity displacement excitation (down-process).

Figure 29 .
Figure 29.The fluctuation term of track irregularity displacement excitation (down-process).

Figure 29 .
Figure 29.The fluctuation term of track irregularity displacement excitation (down-process).

Figure 31 .
Figure 31.Schematic diagram of the photoelectric sensor installation.

Figure 31 .
Figure 31.Schematic diagram of the photoelectric sensor installation.

Figure 31 .
Figure 31.Schematic diagram of the photoelectric sensor installation.

Figure 35 . 30 Figure 35 .
Figure 35.Equal angle resampling points (resampling points are shown at intervals) in the time domain.

Figure 38 .
Figure 38.Comparative analysis of detection results for track irregularities.

Figure 38 .
Figure 38.Comparative analysis of detection results for track irregularities.

Table 2 .
Detection error analysis of track irregularities.

Table 2 .
Detection error analysis of track irregularities.