# Research on the Simulation of Wheelset Response Characteristic Identification of Railway Fastener Loosening

Authors: Wenbai Zhang; Lele Peng; Shubin Zheng; Xun Guo; Yuling Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Rail fastener is a crucial component equipment to ensure the safe operation of the train, and it is very paramount to detect the loose state of the fastener. In this paper, the vertical vibration acceleration signal of wheelset is taken as the research object, and the loose state of fastener is identified by separating and calculating the key IMF energy entropy. Firstly, based on the finite element theory and the principle of multibody dynamics, the rigid-flexible coupling simulation model of vehicle track is established. Then, the vertical vibration acceleration signals of the wheelset under the speed of 200 km/h are obtained by setting the different loosening degrees of the fastener. Finally, we use optimized HHT to process signals, and the orthogonal empirical mode decomposition method (OEMD) is proposed to optimize the orthogonality of the intrinsic mode function, to eliminate the IMF component having poor correlation with the original signal; the Hilbert time spectrum and information entropy theory are combined to calculate the energy entropy of the key IMF, and the HHT energy entropy evaluation algorithm of the vertical acceleration response signal of the train wheelset is proposed. The simulation results show that the HHT energy entropy of 100% fastener looseness is less than 25%, 50%, and 75%, decreasing trend. The algorithm can recognize the looseness of track fastener through the experiment under different working conditions.


---

## Body

## 1. Introduction

At present, the methods of detecting fasteners for railway lines in China mainly include manual inspection, track inspection car, and computer-aided visual inspection. When the vehicle is running on the rail, the interaction between wheel and rail will happen, and when the fastener is damaged or missing, the dynamic parameters of the structure will change to some extent [1]. These changes will be responded to by the vibration signals of some vehicle components in some forms [2, 3].Liu et al. [4] used ANSYS/ LS-DYNA simulation software to establish the vertical coupling vibration model of vehicle-ballastless track-subgrade system, which verified that the sudden change of fastener stiffness had tiny effect on the vertical acceleration of bogie and had obvious effect on the wheelset vibration acceleration; Zhao and Tan [5] established a simulation model of metro vehicle-track flexible body and compared and analyzed the maximum dynamic responses of vehicle-track system under different conditions of fastener failure and different speeds. The research shows that fastener failure has a certain impact on the vibration response of track; Huang [6] analyzed the orthogonality of each intrinsic mode function by numerical simulation and applied the improved HHT method to the damage identification of large structure system; Zhang [7] combined HHT with various theories and applied it to identify rail vibration signals with different fastener looseness under moving load impact. The results show that the method can reflect the change of rail fastener looseness to a certain extent.To sum up, there is no way to judge the loose state of fastener by separating the vibration acceleration of wheelset [8, 9]. This paper presents a method to identify whether the fastener is loose from the response signal of the wheelset before and after the change of the fastener state. It is of great practical significance to identify the fastener.


## 2. Vehicle-Track Rigid-Flexible Coupling Modeling and Dynamic Simulation

### 2.1. Flexible Track Modeling

In order to obtain the wheelset vibration signal, an accurate vehicle-track rigid-flexible coupling model needs to be established [10]. At first, the track system is considered as a flexible body on the basis of the existing rigid vehicle-track model [11, 12]. The slab ballastless track system is shown in Figure 1. The finite element model of rail and rail slab is established in ANSYS. In the multibody dynamics software SIMPACK, force element is used to simulate CA mortar and fastener and assembled into the whole ballastless track system.Figure 1
Slab ballastless track system, which consists of rail (black), fasteners (red), rail slab (green), elastic layer (purple), and concrete subgrade (blue).


#### 2.1.1. Track System Dynamics Model Based on ANSYS

When the rail is modeled, the rail is regarded as a continuous elastic beam. China’s standard 60 kg/section rail parameters are used for modeling and analysis. Because the number of marker points in SIMPACK is limited and if there are too many master DOF (degrees of freedom), it will lead to the incomputable situation in the SIMPACK; therefore, in order to ensure the calculation efficiency, and combined with the purpose of this paper, a 120-meter rail model in ANSYS is established. The distance between two fasteners is 0.6 m, and this length is divided into 20 small units. A total of 201 master DOF nodes are selected. The discrete model of the rail is shown in Figure2(a). The solid model of the rail slab is shown in Figure 2(b).Figure 2
Finite element model of rail and rail slab: (a) model of rail, the distance between two fasteners is 0.6 m, and this length is divided into 20 small units; the total of main DOF nodes is 201. (b) Model of rail slab, continuous elastomer, having an elastic modulus of 3.6e10 Pa, a Poisson’s ratio of 0.1, and a density of 2500 kg/m3, is the solid-45 eight-node-space entity unit. (a) Rail. (b) Rail slab.
(a)(b)


#### 2.1.2. Dynamics Analysis on Rail Subsystem

In this paper, the substructure modal analysis of rail and rail slab needs to reduce the main degree of freedom of the model, and it is reduced by Guyan reduction method in ANSYS. The substructure model of rail and rail slab in this paper is shown in Figure3.Figure 3
Substructure model of rail and rail slab: (a) rail substructure model. (b) Rail board substructure model; the purple triangle mark in the figure is the location of the selected master DOF node; the total number of master DOF should be greater than the order of the subsequent modal analysis, or equal to twice. The predicted deformation direction of rail and rail slab is selected as the master DOF. The location of the load and constraint shall be selected as the main degree of freedom, master DOF is evenly distributed as far as possible, and the nodes at both ends of rail and rail slab and the nodes at fastener positions are selected as the main nodes of degrees of freedom. (a) Rail. (b) Rail slab.
(a)(b)The first 100 modes of the rail and the top 20 modes of the rail slab are calculated, respectively. When selecting the model of rail and rail slab, it is indispensable to accurately reflect the vibration characteristics of rail and rail slab and consider the weight of each model [13]. When the train passes, the vertical vibration is the main vibration of rail, and Figure 4 shows the main mode shapes of the rail and rail slab.Figure 4
Main vibration model of rail and rail slab; vibration model of rail is shown in the 1st, 2nd, and 55th order; vibration model of rail slab is shown in the 1st, 2nd, and 16th order; the order is determined by the length of the model.


### 2.2. Vehicle-Track Rigid-Flexible Coupling Model

#### 2.2.1. FEMBS (Finite Element Multibody System) Interface Program

The basic process flow of flexible data transmission between ANSYS and SIMPACK is shown in Figure5. After modal analysis of rail and rail slab, ANSYS generates geometric model file (.CDB), mass stiffness matrix file (.Sub), and modal mode file (.RST). The interface of FEMBS (finite element multibody system) converts the flexible body data to the standard code (SID) in the format of ASCII readable by SIMPACK, that is, inputting the characteristics of flexible body of rail and rail slab into motion equation to generate SID file. The input information for FEMBS includes DOF and coordinates of nodes, mass attribute, translational and rotational vibration, modal mass matrix, stiffness matrix, damping matrix, and geometric stiffening matrix for initial loads.Figure 5
ANSYS and SIMPACK data transmission flow chart; the left is the flow of ANSYS preprocessing and the right is the flow of SIMPACK/ FEMBS preprocessing.


#### 2.2.2. Realization of Rigid-Flexible Coupling

When assembling the track system in SIMPACK, the fastener and the CA mortar are modeled by linear force element. The flexible rails and the rigid wheels are connected with Hertz spring; however, the contact conditions of wheel and rail in SIMPACK must be rigid wheelset and rigid track and longitudinal relative displacement cannot occur between wheel and rail, and the track model is flexible; the rail does not move with the wheel, so to define a virtual rail between each wheel and flexible rail and define a moving marker point allow the wheel rail force to be transmitted downward to the flexible rail system.The virtual rail body is an object with zero mass and zero-moment of inertia, which only acts as a connection in the dynamic system, and does not affect the other connection structures, but the mass and the moment of inertia of the object in SIMPACK cannot be zero, so the relevant parameters of the virtual rail body are fetched as small as possible, which is assumed to be 1.0e − 6 in this paper, whose impact on wheel rail force transmission is negligible.In order to ensure that each virtual rail body moves longitudinally along the rail with its corresponding wheelset, the moving marker points on the virtual rail body and the flexible rail are required to define the restraint to hold them, and the wheelset is assembled to the virtual rail body through the hinges, so that the virtual rail body is equivalent to a part of the rail and the vehicle can run along the rail. The wheel-rail rigid-flexible coupling model is shown in Figure6(a). In simulation, the wheel rail force is calculated through Hertz contact between the virtual rail body and rigid wheel. Data exchange between the virtual rail body and the flexible rail is done through the deformation coordination condition and the force balance condition [14, 15]. According to the above method, the vehicle and the flexible track are assembled together in the SIMPACK. The established vehicle-track rigid-flexible coupling model is shown in Figure 6(b).Figure 6
Coupling mode: (a) wheel (blue), hertz contract (green), virtual rail (red), constraint and force element (purple), move marker (yellow). (b) There is a motor vehicle with two bogies in the coupling model. (a) Wheel-rail coupling model. (b) Vehicle-track rigid-flexible coupling model.
(a)(b)


#### 2.2.3. Model Verification When Fasteners Are Not Loose

Before simulation, firstly we should verify whether the initial state of vehicle-track rigid-flexible coupling model is balanced, which is similar to the verification of rigid body vehicle model. In SIMPACK, we verify the model by calculating the nominal force maximum acceleration value of initial state. The maximum acceleration is calculated to be6.811∗10−6 m/s2, which is less than the evaluation order of magnitude 10−4 m/s2, so it can be determined that the built model is balanced in the initial state and the coupling model is correct.In [16], the actual measured values of the vehicle vibration acceleration when the high-speed vehicle runs through under different conditions of track irregularity are given, and the simulation calculation is carried out with the VICT simulation software under the same conditions. In this paper, the same test conditions are set in the model to carry out the simulation calculation of vehicle dynamics, and the simulation results are compared with the above results. See Table 1 for comparison. We can see that the simulation results of this paper are not much different from the test results and the VICT simulation results in [17]. The reason for the slight difference is because the vehicle simulation parameters of this paper are different from the other article, and the vertical correctness of the simulation model is verified, which lays a foundation for the subsequent simulation.Table 1
The comparison of max value from experiment [16], VICT simulation [17], and current coupled model.
Test conditionsLongitudinal irregularityWavelengthλ (m)1012122424Wave deptha (mm)10991620Test speedv (km/h)160135150160160Max measured value in the reference (g)0.120.060.080.120.13Max value of the VICT simulation (g)0.1040.0780.0850.0960.120Value of simulation in the current model (g)0.1080.0670.0790.1100.111


### 2.3. Vehicle Dynamics Simulation

#### 2.3.1. SIMPACK Simulation Excitation Model

SIMPACK can generate excitation by means of inputting the coefficient of the track spectral density formula and the method of converting the measured data. The software can be set up for two types of excitations: one is track-related; the other is rail-related. Due to the fact that there is no unified track spectrum standard in our country, the low-interference track spectrum of German high-speed railway is adopted in this paper. The vertical track irregularities generated in SIMPACK are shown in Figure7.Figure 7
Track vertical irregularity. In (a) and (b), the first figure is the vertical irregularity imposed on the track, the second figure is the irregularity under the speed excitation, and the third figure is the irregularity under the acceleration excitation. (a) Left track irregularity. (b) Right track irregularity.
(a)(b)


#### 2.3.2. Simulation Conditions

When two fasteners are loosened on the same section of the same track, and two fasteners are loosened continuously, the impact on the vehicle and track system is relatively large [18, 19]. Therefore, three typical simulation conditions are set up for the looseness of the fastener, as shown in the diagram. The cross line is the fastener loosening position in the picture, and there is one loose fastener on the first condition, there are two loose fasteners on the same section of the track in the second condition, and there are two loose fasteners on the same rail in the third condition.Selecting the fastener classic stiffness of 50 kN/mm, the vehicle passes through the fastener in different degrees of looseness at 200 km/h speed (not loose, 25% loose, 50% loose, 75% loose, and completely loose), the looseness of the fastener is simulated by changing the stiffness value of the fastener force element in the SIMPACK, the corresponding stiffness value of the fastener loosening degree, and the equivalent stiffness value of the mortar unit length as shown in Table2, and the looseness degree of fastener is the same in the three working conditions, shown in Figure 8. The vertical acceleration of wheelset and vehicle body is calculated by simulation.Table 2
The corresponding stiffness value of the loosening degree of the fastener and the equivalent stiffness value of the length of the mortar unit.
Degrees of fastener stiffness looseness (%)Vertical, horizontal, and longitudinal stiffness value of the fastener (kN/mm)Vertical, horizontal, and longitudinal damping value of the fastener (kN·s/mm)Equivalent stiffness of CA mortar unit length (N/m3)Equivalent damping of CA mortar unit length (N·s/m)05075/60/601.25∗1093.46∗1042537.550257512.51000Figure 8
Loose fastener condition. The red square and cross is the position of loose fastener: one loose fastener in condition 1, two opposite in condition 2, and two consecutive in condition 3. The looseness of fastener could be set separately. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)Figure 9
Wheelset vertical acceleration signal with different degrees of fastener looseness in condition 1, looseness 0% (black), 25% (brown), 50% (green), 75% (purple), and 100% (blue).Figure 10
Wheelset vertical acceleration with different degrees of fastener looseness in condition 2.Figure 11
Wheelset vertical acceleration with different degrees of fastener looseness in condition 3.


#### 2.3.3. Simulation Results of Wheelset Vertical Acceleration

(1)
Wheelset vertical acceleration in condition 1The vehicle response is the same before the vehicle arrives at the loosened part of the fastener [20, 21], so only the vertical acceleration response of the 0.6 s∼1.7 s is given. The vertical acceleration comparison diagram of the wheelset in different degrees of fastener loosening is given, as shown in Figure 9; the acceleration at about 0.72 s∼0.76 s starts changing; with the increase of the degree of fastener looseness, the vibration signal also changes obviously. The vertical acceleration of the wheelset with loosened fasteners has a variation of 1 m/s2 to 6 m/s2, when the fastener is completely loose, the change reaches the maximum. Wheelset vertical acceleration comparison of five looseness degrees in condition 2 and condition 3 is shown in Figures 10 and 11. The looseness of the fastener will increase the vibration and displacement of the rail, which is equivalent to the influence of the irregularity on the train in the running line.(2)
Wheelset vertical acceleration in condition 2(3)
Wheelset vertical acceleration in condition 3


## 2.1. Flexible Track Modeling

In order to obtain the wheelset vibration signal, an accurate vehicle-track rigid-flexible coupling model needs to be established [10]. At first, the track system is considered as a flexible body on the basis of the existing rigid vehicle-track model [11, 12]. The slab ballastless track system is shown in Figure 1. The finite element model of rail and rail slab is established in ANSYS. In the multibody dynamics software SIMPACK, force element is used to simulate CA mortar and fastener and assembled into the whole ballastless track system.Figure 1
Slab ballastless track system, which consists of rail (black), fasteners (red), rail slab (green), elastic layer (purple), and concrete subgrade (blue).


### 2.1.1. Track System Dynamics Model Based on ANSYS

When the rail is modeled, the rail is regarded as a continuous elastic beam. China’s standard 60 kg/section rail parameters are used for modeling and analysis. Because the number of marker points in SIMPACK is limited and if there are too many master DOF (degrees of freedom), it will lead to the incomputable situation in the SIMPACK; therefore, in order to ensure the calculation efficiency, and combined with the purpose of this paper, a 120-meter rail model in ANSYS is established. The distance between two fasteners is 0.6 m, and this length is divided into 20 small units. A total of 201 master DOF nodes are selected. The discrete model of the rail is shown in Figure2(a). The solid model of the rail slab is shown in Figure 2(b).Figure 2
Finite element model of rail and rail slab: (a) model of rail, the distance between two fasteners is 0.6 m, and this length is divided into 20 small units; the total of main DOF nodes is 201. (b) Model of rail slab, continuous elastomer, having an elastic modulus of 3.6e10 Pa, a Poisson’s ratio of 0.1, and a density of 2500 kg/m3, is the solid-45 eight-node-space entity unit. (a) Rail. (b) Rail slab.
(a)(b)


### 2.1.2. Dynamics Analysis on Rail Subsystem

In this paper, the substructure modal analysis of rail and rail slab needs to reduce the main degree of freedom of the model, and it is reduced by Guyan reduction method in ANSYS. The substructure model of rail and rail slab in this paper is shown in Figure3.Figure 3
Substructure model of rail and rail slab: (a) rail substructure model. (b) Rail board substructure model; the purple triangle mark in the figure is the location of the selected master DOF node; the total number of master DOF should be greater than the order of the subsequent modal analysis, or equal to twice. The predicted deformation direction of rail and rail slab is selected as the master DOF. The location of the load and constraint shall be selected as the main degree of freedom, master DOF is evenly distributed as far as possible, and the nodes at both ends of rail and rail slab and the nodes at fastener positions are selected as the main nodes of degrees of freedom. (a) Rail. (b) Rail slab.
(a)(b)The first 100 modes of the rail and the top 20 modes of the rail slab are calculated, respectively. When selecting the model of rail and rail slab, it is indispensable to accurately reflect the vibration characteristics of rail and rail slab and consider the weight of each model [13]. When the train passes, the vertical vibration is the main vibration of rail, and Figure 4 shows the main mode shapes of the rail and rail slab.Figure 4
Main vibration model of rail and rail slab; vibration model of rail is shown in the 1st, 2nd, and 55th order; vibration model of rail slab is shown in the 1st, 2nd, and 16th order; the order is determined by the length of the model.


## 2.1.1. Track System Dynamics Model Based on ANSYS

When the rail is modeled, the rail is regarded as a continuous elastic beam. China’s standard 60 kg/section rail parameters are used for modeling and analysis. Because the number of marker points in SIMPACK is limited and if there are too many master DOF (degrees of freedom), it will lead to the incomputable situation in the SIMPACK; therefore, in order to ensure the calculation efficiency, and combined with the purpose of this paper, a 120-meter rail model in ANSYS is established. The distance between two fasteners is 0.6 m, and this length is divided into 20 small units. A total of 201 master DOF nodes are selected. The discrete model of the rail is shown in Figure2(a). The solid model of the rail slab is shown in Figure 2(b).Figure 2
Finite element model of rail and rail slab: (a) model of rail, the distance between two fasteners is 0.6 m, and this length is divided into 20 small units; the total of main DOF nodes is 201. (b) Model of rail slab, continuous elastomer, having an elastic modulus of 3.6e10 Pa, a Poisson’s ratio of 0.1, and a density of 2500 kg/m3, is the solid-45 eight-node-space entity unit. (a) Rail. (b) Rail slab.
(a)(b)


## 2.1.2. Dynamics Analysis on Rail Subsystem

In this paper, the substructure modal analysis of rail and rail slab needs to reduce the main degree of freedom of the model, and it is reduced by Guyan reduction method in ANSYS. The substructure model of rail and rail slab in this paper is shown in Figure3.Figure 3
Substructure model of rail and rail slab: (a) rail substructure model. (b) Rail board substructure model; the purple triangle mark in the figure is the location of the selected master DOF node; the total number of master DOF should be greater than the order of the subsequent modal analysis, or equal to twice. The predicted deformation direction of rail and rail slab is selected as the master DOF. The location of the load and constraint shall be selected as the main degree of freedom, master DOF is evenly distributed as far as possible, and the nodes at both ends of rail and rail slab and the nodes at fastener positions are selected as the main nodes of degrees of freedom. (a) Rail. (b) Rail slab.
(a)(b)The first 100 modes of the rail and the top 20 modes of the rail slab are calculated, respectively. When selecting the model of rail and rail slab, it is indispensable to accurately reflect the vibration characteristics of rail and rail slab and consider the weight of each model [13]. When the train passes, the vertical vibration is the main vibration of rail, and Figure 4 shows the main mode shapes of the rail and rail slab.Figure 4
Main vibration model of rail and rail slab; vibration model of rail is shown in the 1st, 2nd, and 55th order; vibration model of rail slab is shown in the 1st, 2nd, and 16th order; the order is determined by the length of the model.


## 2.2. Vehicle-Track Rigid-Flexible Coupling Model

### 2.2.1. FEMBS (Finite Element Multibody System) Interface Program

The basic process flow of flexible data transmission between ANSYS and SIMPACK is shown in Figure5. After modal analysis of rail and rail slab, ANSYS generates geometric model file (.CDB), mass stiffness matrix file (.Sub), and modal mode file (.RST). The interface of FEMBS (finite element multibody system) converts the flexible body data to the standard code (SID) in the format of ASCII readable by SIMPACK, that is, inputting the characteristics of flexible body of rail and rail slab into motion equation to generate SID file. The input information for FEMBS includes DOF and coordinates of nodes, mass attribute, translational and rotational vibration, modal mass matrix, stiffness matrix, damping matrix, and geometric stiffening matrix for initial loads.Figure 5
ANSYS and SIMPACK data transmission flow chart; the left is the flow of ANSYS preprocessing and the right is the flow of SIMPACK/ FEMBS preprocessing.


### 2.2.2. Realization of Rigid-Flexible Coupling

When assembling the track system in SIMPACK, the fastener and the CA mortar are modeled by linear force element. The flexible rails and the rigid wheels are connected with Hertz spring; however, the contact conditions of wheel and rail in SIMPACK must be rigid wheelset and rigid track and longitudinal relative displacement cannot occur between wheel and rail, and the track model is flexible; the rail does not move with the wheel, so to define a virtual rail between each wheel and flexible rail and define a moving marker point allow the wheel rail force to be transmitted downward to the flexible rail system.The virtual rail body is an object with zero mass and zero-moment of inertia, which only acts as a connection in the dynamic system, and does not affect the other connection structures, but the mass and the moment of inertia of the object in SIMPACK cannot be zero, so the relevant parameters of the virtual rail body are fetched as small as possible, which is assumed to be 1.0e − 6 in this paper, whose impact on wheel rail force transmission is negligible.In order to ensure that each virtual rail body moves longitudinally along the rail with its corresponding wheelset, the moving marker points on the virtual rail body and the flexible rail are required to define the restraint to hold them, and the wheelset is assembled to the virtual rail body through the hinges, so that the virtual rail body is equivalent to a part of the rail and the vehicle can run along the rail. The wheel-rail rigid-flexible coupling model is shown in Figure6(a). In simulation, the wheel rail force is calculated through Hertz contact between the virtual rail body and rigid wheel. Data exchange between the virtual rail body and the flexible rail is done through the deformation coordination condition and the force balance condition [14, 15]. According to the above method, the vehicle and the flexible track are assembled together in the SIMPACK. The established vehicle-track rigid-flexible coupling model is shown in Figure 6(b).Figure 6
Coupling mode: (a) wheel (blue), hertz contract (green), virtual rail (red), constraint and force element (purple), move marker (yellow). (b) There is a motor vehicle with two bogies in the coupling model. (a) Wheel-rail coupling model. (b) Vehicle-track rigid-flexible coupling model.
(a)(b)


### 2.2.3. Model Verification When Fasteners Are Not Loose

Before simulation, firstly we should verify whether the initial state of vehicle-track rigid-flexible coupling model is balanced, which is similar to the verification of rigid body vehicle model. In SIMPACK, we verify the model by calculating the nominal force maximum acceleration value of initial state. The maximum acceleration is calculated to be6.811∗10−6 m/s2, which is less than the evaluation order of magnitude 10−4 m/s2, so it can be determined that the built model is balanced in the initial state and the coupling model is correct.In [16], the actual measured values of the vehicle vibration acceleration when the high-speed vehicle runs through under different conditions of track irregularity are given, and the simulation calculation is carried out with the VICT simulation software under the same conditions. In this paper, the same test conditions are set in the model to carry out the simulation calculation of vehicle dynamics, and the simulation results are compared with the above results. See Table 1 for comparison. We can see that the simulation results of this paper are not much different from the test results and the VICT simulation results in [17]. The reason for the slight difference is because the vehicle simulation parameters of this paper are different from the other article, and the vertical correctness of the simulation model is verified, which lays a foundation for the subsequent simulation.Table 1
The comparison of max value from experiment [16], VICT simulation [17], and current coupled model.
Test conditionsLongitudinal irregularityWavelengthλ (m)1012122424Wave deptha (mm)10991620Test speedv (km/h)160135150160160Max measured value in the reference (g)0.120.060.080.120.13Max value of the VICT simulation (g)0.1040.0780.0850.0960.120Value of simulation in the current model (g)0.1080.0670.0790.1100.111


## 2.2.1. FEMBS (Finite Element Multibody System) Interface Program

The basic process flow of flexible data transmission between ANSYS and SIMPACK is shown in Figure5. After modal analysis of rail and rail slab, ANSYS generates geometric model file (.CDB), mass stiffness matrix file (.Sub), and modal mode file (.RST). The interface of FEMBS (finite element multibody system) converts the flexible body data to the standard code (SID) in the format of ASCII readable by SIMPACK, that is, inputting the characteristics of flexible body of rail and rail slab into motion equation to generate SID file. The input information for FEMBS includes DOF and coordinates of nodes, mass attribute, translational and rotational vibration, modal mass matrix, stiffness matrix, damping matrix, and geometric stiffening matrix for initial loads.Figure 5
ANSYS and SIMPACK data transmission flow chart; the left is the flow of ANSYS preprocessing and the right is the flow of SIMPACK/ FEMBS preprocessing.


## 2.2.2. Realization of Rigid-Flexible Coupling

When assembling the track system in SIMPACK, the fastener and the CA mortar are modeled by linear force element. The flexible rails and the rigid wheels are connected with Hertz spring; however, the contact conditions of wheel and rail in SIMPACK must be rigid wheelset and rigid track and longitudinal relative displacement cannot occur between wheel and rail, and the track model is flexible; the rail does not move with the wheel, so to define a virtual rail between each wheel and flexible rail and define a moving marker point allow the wheel rail force to be transmitted downward to the flexible rail system.The virtual rail body is an object with zero mass and zero-moment of inertia, which only acts as a connection in the dynamic system, and does not affect the other connection structures, but the mass and the moment of inertia of the object in SIMPACK cannot be zero, so the relevant parameters of the virtual rail body are fetched as small as possible, which is assumed to be 1.0e − 6 in this paper, whose impact on wheel rail force transmission is negligible.In order to ensure that each virtual rail body moves longitudinally along the rail with its corresponding wheelset, the moving marker points on the virtual rail body and the flexible rail are required to define the restraint to hold them, and the wheelset is assembled to the virtual rail body through the hinges, so that the virtual rail body is equivalent to a part of the rail and the vehicle can run along the rail. The wheel-rail rigid-flexible coupling model is shown in Figure6(a). In simulation, the wheel rail force is calculated through Hertz contact between the virtual rail body and rigid wheel. Data exchange between the virtual rail body and the flexible rail is done through the deformation coordination condition and the force balance condition [14, 15]. According to the above method, the vehicle and the flexible track are assembled together in the SIMPACK. The established vehicle-track rigid-flexible coupling model is shown in Figure 6(b).Figure 6
Coupling mode: (a) wheel (blue), hertz contract (green), virtual rail (red), constraint and force element (purple), move marker (yellow). (b) There is a motor vehicle with two bogies in the coupling model. (a) Wheel-rail coupling model. (b) Vehicle-track rigid-flexible coupling model.
(a)(b)


## 2.2.3. Model Verification When Fasteners Are Not Loose

Before simulation, firstly we should verify whether the initial state of vehicle-track rigid-flexible coupling model is balanced, which is similar to the verification of rigid body vehicle model. In SIMPACK, we verify the model by calculating the nominal force maximum acceleration value of initial state. The maximum acceleration is calculated to be6.811∗10−6 m/s2, which is less than the evaluation order of magnitude 10−4 m/s2, so it can be determined that the built model is balanced in the initial state and the coupling model is correct.In [16], the actual measured values of the vehicle vibration acceleration when the high-speed vehicle runs through under different conditions of track irregularity are given, and the simulation calculation is carried out with the VICT simulation software under the same conditions. In this paper, the same test conditions are set in the model to carry out the simulation calculation of vehicle dynamics, and the simulation results are compared with the above results. See Table 1 for comparison. We can see that the simulation results of this paper are not much different from the test results and the VICT simulation results in [17]. The reason for the slight difference is because the vehicle simulation parameters of this paper are different from the other article, and the vertical correctness of the simulation model is verified, which lays a foundation for the subsequent simulation.Table 1
The comparison of max value from experiment [16], VICT simulation [17], and current coupled model.
Test conditionsLongitudinal irregularityWavelengthλ (m)1012122424Wave deptha (mm)10991620Test speedv (km/h)160135150160160Max measured value in the reference (g)0.120.060.080.120.13Max value of the VICT simulation (g)0.1040.0780.0850.0960.120Value of simulation in the current model (g)0.1080.0670.0790.1100.111


## 2.3. Vehicle Dynamics Simulation

### 2.3.1. SIMPACK Simulation Excitation Model

SIMPACK can generate excitation by means of inputting the coefficient of the track spectral density formula and the method of converting the measured data. The software can be set up for two types of excitations: one is track-related; the other is rail-related. Due to the fact that there is no unified track spectrum standard in our country, the low-interference track spectrum of German high-speed railway is adopted in this paper. The vertical track irregularities generated in SIMPACK are shown in Figure7.Figure 7
Track vertical irregularity. In (a) and (b), the first figure is the vertical irregularity imposed on the track, the second figure is the irregularity under the speed excitation, and the third figure is the irregularity under the acceleration excitation. (a) Left track irregularity. (b) Right track irregularity.
(a)(b)


### 2.3.2. Simulation Conditions

When two fasteners are loosened on the same section of the same track, and two fasteners are loosened continuously, the impact on the vehicle and track system is relatively large [18, 19]. Therefore, three typical simulation conditions are set up for the looseness of the fastener, as shown in the diagram. The cross line is the fastener loosening position in the picture, and there is one loose fastener on the first condition, there are two loose fasteners on the same section of the track in the second condition, and there are two loose fasteners on the same rail in the third condition.Selecting the fastener classic stiffness of 50 kN/mm, the vehicle passes through the fastener in different degrees of looseness at 200 km/h speed (not loose, 25% loose, 50% loose, 75% loose, and completely loose), the looseness of the fastener is simulated by changing the stiffness value of the fastener force element in the SIMPACK, the corresponding stiffness value of the fastener loosening degree, and the equivalent stiffness value of the mortar unit length as shown in Table2, and the looseness degree of fastener is the same in the three working conditions, shown in Figure 8. The vertical acceleration of wheelset and vehicle body is calculated by simulation.Table 2
The corresponding stiffness value of the loosening degree of the fastener and the equivalent stiffness value of the length of the mortar unit.
Degrees of fastener stiffness looseness (%)Vertical, horizontal, and longitudinal stiffness value of the fastener (kN/mm)Vertical, horizontal, and longitudinal damping value of the fastener (kN·s/mm)Equivalent stiffness of CA mortar unit length (N/m3)Equivalent damping of CA mortar unit length (N·s/m)05075/60/601.25∗1093.46∗1042537.550257512.51000Figure 8
Loose fastener condition. The red square and cross is the position of loose fastener: one loose fastener in condition 1, two opposite in condition 2, and two consecutive in condition 3. The looseness of fastener could be set separately. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)Figure 9
Wheelset vertical acceleration signal with different degrees of fastener looseness in condition 1, looseness 0% (black), 25% (brown), 50% (green), 75% (purple), and 100% (blue).Figure 10
Wheelset vertical acceleration with different degrees of fastener looseness in condition 2.Figure 11
Wheelset vertical acceleration with different degrees of fastener looseness in condition 3.


### 2.3.3. Simulation Results of Wheelset Vertical Acceleration

(1)
Wheelset vertical acceleration in condition 1The vehicle response is the same before the vehicle arrives at the loosened part of the fastener [20, 21], so only the vertical acceleration response of the 0.6 s∼1.7 s is given. The vertical acceleration comparison diagram of the wheelset in different degrees of fastener loosening is given, as shown in Figure 9; the acceleration at about 0.72 s∼0.76 s starts changing; with the increase of the degree of fastener looseness, the vibration signal also changes obviously. The vertical acceleration of the wheelset with loosened fasteners has a variation of 1 m/s2 to 6 m/s2, when the fastener is completely loose, the change reaches the maximum. Wheelset vertical acceleration comparison of five looseness degrees in condition 2 and condition 3 is shown in Figures 10 and 11. The looseness of the fastener will increase the vibration and displacement of the rail, which is equivalent to the influence of the irregularity on the train in the running line.(2)
Wheelset vertical acceleration in condition 2(3)
Wheelset vertical acceleration in condition 3


## 2.3.1. SIMPACK Simulation Excitation Model

SIMPACK can generate excitation by means of inputting the coefficient of the track spectral density formula and the method of converting the measured data. The software can be set up for two types of excitations: one is track-related; the other is rail-related. Due to the fact that there is no unified track spectrum standard in our country, the low-interference track spectrum of German high-speed railway is adopted in this paper. The vertical track irregularities generated in SIMPACK are shown in Figure7.Figure 7
Track vertical irregularity. In (a) and (b), the first figure is the vertical irregularity imposed on the track, the second figure is the irregularity under the speed excitation, and the third figure is the irregularity under the acceleration excitation. (a) Left track irregularity. (b) Right track irregularity.
(a)(b)


## 2.3.2. Simulation Conditions

When two fasteners are loosened on the same section of the same track, and two fasteners are loosened continuously, the impact on the vehicle and track system is relatively large [18, 19]. Therefore, three typical simulation conditions are set up for the looseness of the fastener, as shown in the diagram. The cross line is the fastener loosening position in the picture, and there is one loose fastener on the first condition, there are two loose fasteners on the same section of the track in the second condition, and there are two loose fasteners on the same rail in the third condition.Selecting the fastener classic stiffness of 50 kN/mm, the vehicle passes through the fastener in different degrees of looseness at 200 km/h speed (not loose, 25% loose, 50% loose, 75% loose, and completely loose), the looseness of the fastener is simulated by changing the stiffness value of the fastener force element in the SIMPACK, the corresponding stiffness value of the fastener loosening degree, and the equivalent stiffness value of the mortar unit length as shown in Table2, and the looseness degree of fastener is the same in the three working conditions, shown in Figure 8. The vertical acceleration of wheelset and vehicle body is calculated by simulation.Table 2
The corresponding stiffness value of the loosening degree of the fastener and the equivalent stiffness value of the length of the mortar unit.
Degrees of fastener stiffness looseness (%)Vertical, horizontal, and longitudinal stiffness value of the fastener (kN/mm)Vertical, horizontal, and longitudinal damping value of the fastener (kN·s/mm)Equivalent stiffness of CA mortar unit length (N/m3)Equivalent damping of CA mortar unit length (N·s/m)05075/60/601.25∗1093.46∗1042537.550257512.51000Figure 8
Loose fastener condition. The red square and cross is the position of loose fastener: one loose fastener in condition 1, two opposite in condition 2, and two consecutive in condition 3. The looseness of fastener could be set separately. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)Figure 9
Wheelset vertical acceleration signal with different degrees of fastener looseness in condition 1, looseness 0% (black), 25% (brown), 50% (green), 75% (purple), and 100% (blue).Figure 10
Wheelset vertical acceleration with different degrees of fastener looseness in condition 2.Figure 11
Wheelset vertical acceleration with different degrees of fastener looseness in condition 3.


## 2.3.3. Simulation Results of Wheelset Vertical Acceleration

(1)
Wheelset vertical acceleration in condition 1The vehicle response is the same before the vehicle arrives at the loosened part of the fastener [20, 21], so only the vertical acceleration response of the 0.6 s∼1.7 s is given. The vertical acceleration comparison diagram of the wheelset in different degrees of fastener loosening is given, as shown in Figure 9; the acceleration at about 0.72 s∼0.76 s starts changing; with the increase of the degree of fastener looseness, the vibration signal also changes obviously. The vertical acceleration of the wheelset with loosened fasteners has a variation of 1 m/s2 to 6 m/s2, when the fastener is completely loose, the change reaches the maximum. Wheelset vertical acceleration comparison of five looseness degrees in condition 2 and condition 3 is shown in Figures 10 and 11. The looseness of the fastener will increase the vibration and displacement of the rail, which is equivalent to the influence of the irregularity on the train in the running line.(2)
Wheelset vertical acceleration in condition 2(3)
Wheelset vertical acceleration in condition 3


## 3. Identification of Fastener Loosening Feature Based on Response Signal

### 3.1. Identification of Fasteners State in Varying Degrees

HHT has good reliability in dealing with nonstationary and nonlinear signals, detecting structural faults. However, it is not based on a complete theory and has energy leakage problems. In order to solve this problem, this paper optimizes the algorithm according to the theory of orthogonality and correlation and combines the improved Hilbert Huang algorithm with the theory of information entropy [22, 23], which is applied to the identification of rail fastener looseness.


#### 3.1.1. Hilbert–Huang Transform

The crux of HHT algorithm is EMD decomposition. The result of EMD decomposition directly affects the accuracy of subsequent signal processing, so it is necessary to ensure that the IMF component of the EMD decomposition should possess completeness and orthogonality, so that no energy leakage in the decomposition can be guaranteed [24, 25]. From the practical point of view, Huang et al. consider that there exists orthogonality between all IMFs decomposed by EMD, but there is no rigorous theoretical derivation proving that they are rigorously orthogonal in the overall situation[29, 30].Huang proposed two orthogonality indicators to measure the orthogonality between IMF components, namely, the overall orthogonality index (IOT) and the orthogonality index between two components (IOjk) [31]. When all the IMF components are rigorously orthogonal to each other, the IOT and IOjk should be zero.(1)IOT=∑j=1n+1∑k=1k≠jn+1∑i=1Nimfjiimfki∑i=1Nxi2,IOjk=∑i=1Nimfjiimfki∑i=1Nimfji2+imfki2.In the formula,xi is the ith original signal, imfji is the jth IMF component of xi, imfk is the kth IMF component of xi, and j is not equal to k. In addition, the energy can be used to measure the degree of orthogonality between IMFs. The energy of the original signal X(t) = [x1, x2, …, xn] is(2)Ex=∫0TX2tdt=∑i=1Nxi2.The energy of each IMF is(3)Ej=∫0Timfj2tdt=∑i=1Nimfji2t,j=1,…,n+1.In the formula, imfj (t) is the jth IMF component; if all the IMFs are strictly orthogonal to each other, Ex = Ej and the leaking energy between the IMFs is zero, that is,(4)EIOT=∑j=1n+1Ejt=Ex,where EIOT is the energy of IOT, total energy after decomposition which is constant.(5)Ejk=∫0Timfjtimfktdt=∑i=1Nimfjiimfki=0,j,k=1,…,n+1;j≠k.Ejk is the leaking energy between imfj and imfk. From the processing steps of EMD decomposition, we can see that, in the actual decomposition process, the mean value curve is obtained by fitting the approximate value, and the decomposed IMF is a part of the original signal, so the different IMF components are approximately orthogonal [29, 30]. Huang’s theory proves that the EMD orthogonality is in the order of 10−3 to 10−2, and the IMFs are not strictly orthogonal, resulting in the phenomenon of energy leakage and modal aliasing during signal analysis, which will bring errors in the later analysis, making the result inaccurate. In order to eliminate energy leakage, it is necessary to ensure that the orthonormal mode functions are strictly orthogonal to each other [31]. Therefore, the traditional EMD algorithm is improved.


#### 3.1.2. OEMD Algorithm

In order to ensure the strict orthogonality between the IMF components after EMD decomposition [32, 33], this paper orthogonalizes the IMF components; the basic process is shown as follows:(1)
The original signalX(t) is decomposed into the form of the sum of multiple IMF components ci(t) (i = 1, 2, …, n) and a residual component rn (t) through EMD:(6)Xt=∑j=1ncjt+rnt.(2)
Assume imf1 (t) = c1 (t) and it is the first orthogonalized IMF component of the original signal, which is the highest frequency component of the IMF.(3)
In the process of decomposition to obtain the second IMF componentc2 (t), we can see that it is impossible to ensure the orthogonality between c2 (t) and imf1 (t). Therefore, we need to remove imf1 (t) from c2 (t); that is,(7)imf2t=c2t−β21imf1t,imf2 (t) is the second orthogonalized IMF component of the original signal X (t), and β21 is the orthogonalization coefficient between c2 (t) and imf1 (t). Each side of equations (2)–(7) is multiplied by imf1 (t) so as to integrate t. Since imf1 (t) and imf2 (t) are orthogonal, it can be concluded that(8)∫0Timf1timf2tdt=∫0Tc2timf1tdt−β21∫0Timf12tdt=0,(9)β21=c2⟶Timf1⟶imf1⟶Timf1⟶.The numerator and denominator in equations (2)–(9) are inner product operations of two vectors.(4)
By analogy, thej + 1th orthogonal component, imfj+1 (t), of the original signal can be obtained by eliminating the imfi (t) (i = 1, 2, …, j) components in the j + 1th IMF component cj+1 (t). That is,(10)imfj+1t=cj+1t−∑i=1jβj+1,iimfit.Each side of the above equation is multiplied by imfi (t) (i ≤ j), so as to integrate t. Owing to the fact that imfi (t) and imfj+1 (t) are orthogonal, it can be concluded that(11)βj+1,i=cj+1⟶Timfi⟶imfi⟶Timfi⟶.After the above calculation, the original signalX (t) is decomposed into the form of the sum of the orthogonalized IMF components and the residual components, that is,(12)Xt=c1t+c2t+⋯+cnt+rnt=imf1t+imf2t+β21imf1t+imf3t+β31imf1t+β32c2t+⋯+imfnt+βn1imf1t+βn2imf2t+⋯+βn,n−1imfn−1t+rnt=1+β21+β31+β41+⋯+βn1imf1t+1+β32+β42+⋯+βn2imf2t+⋯+1+βn,n−1imfn−1t+imfnt+rnt=imf1∗t+imf2∗t+⋯+imfn∗t+rnt=∑j=1nimfj∗t+rnt=∑j=1najimfjt+rnt.In equations (2)–(12), aj=∑i=jnβi,j,j=1,2,…,n,βi,j=1i=j. From the above calculation, we can see that there is a strict orthogonality in imfj (t), so the linear transformation on imfj (t) will not change its orthogonality, and imfj∗t is also strictly orthogonal.The integration of the intrinsic mode components after the decomposition of the above EMD is called the orthogonal empirical mode decomposition (OEMD). According to the different orthogonal sequence of IMF components, the OEMD can be divided into three kinds. The way to orthogonalize the high frequency IMF component before orthogonalizing the low frequency IMF component is called OEMD1, the way to deal with IMF in reverse order is called OEMD2, and the way to perform orthogonalization from any IMF is called OEMD3. The essence of OEMD is still EMD; only the IMF components decomposed by EMD are orthogonalized and then reorganized to achieve strict orthogonality of each IMF component.


### 3.2. OEMD Decomposition of the Vertical Acceleration Signals of Wheelset

The EMD is used to decompose the vertical acceleration signal (Figure9) of the wheelset without loosening of the fastener under the working condition 1 and obtain 7 IMF components and 1 residual component, as shown in Figure 12. The same signals are processed by OEMD1 and OEMD2, respectively; the results of 7 IMF components and 1 residual component can be obtained.Figure 12
EMD decomposition results of simulation signal (no loose fasteners).The overall orthogonal index obtained by the decomposition of EMD is 0.7546, which can be calculated by formula (1). The overall orthogonal indices obtained by OEMD1 and OEMD2 are 0.006463 and 0.007347, respectively. It can be seen that the overall orthogonal index obtained by the OEMD is 2 orders of magnitude higher than that of the overall EMD. There is little difference in accuracy of the orthogonality of the two orthogonal EMD algorithms. Since the 8th IMF is a residual component with a very small value, it has no orthogonality with the first 7 components.The orthogonality between the 7 IMF components after the decomposition of EMD, OEMD1, and OEMD2 is calculated, respectively, by formula (1); the results of the calculation are shown in Tables 3 and 4. The upper triangular data of Table 3 are the orthogonal indexes between IMF components obtained by EMD decomposition, and the lower triangular data are the orthogonal indexes between IMF components obtained by OEMD1 decomposition. The upper triangular data of Table 4 are the orthogonal indexes obtained by EMD decomposition, and the lower triangular data are the orthogonal indexes obtained by OEMD2 decomposition.Table 3
Orthogonal index of each component in EMD and OEMD2.
IMF1234567810.55.60e − 29.82e − 31.17e − 26.07e − 32.16e − 31.80e − 23.63e − 225.68e − 190.56.29e − 26.65e − 31.13e − 21.83e − 27.31e − 21.59e − 232.45e − 178.16e − 170.53.80e − 21.92e − 21.02e − 21.11e − 24.79e − 342.16e − 179.92e − 181.48e − 160.55.84e − 22.36e − 36.89e − 33.48e − 353.31e − 182.11e − 171.91e − 171.22e − 160.56.19e − 42.80e − 34.35e − 365.56e − 182.65e − 181.75e − 174.80e − 181.18e − 170.58.86e − 24.76e − 374.90e − 187.17e − 181.78e − 179.65e − 183.10e − 171.93e − 160.51.53e − 283.63e − 21.94e − 27.18e − 33.62e − 32.73e − 44.38e − 21.08e − 20.5Table 4
Orthogonal index of each component in EMD and OEMD1.
IMF1234567810.58.86e − 22.80e − 36.89e − 31.11e − 27.31e − 21.80e − 21.70e − 224.76e − 170.56.19e − 42.36e − 31.02e − 21.83e − 22.16e − 31.07e − 233.15e − 182.76e − 170.55.84e − 21.92e − 21.13e − 26.07e − 31.03e − 343.10e − 171.07e − 174.52e − 170.53.80e − 26.65e − 31.17e − 23.37e − 452.27e − 171.17e − 173.00e − 178.13e − 180.56.29e − 29.82e − 31.89e − 361.64e − 181.24e − 173.46e − 171.37e − 172.66e − 160.55.60e − 21.81e − 273.87e − 201.97e − 182.60e − 184.42e − 175.67e − 171.23e − 160.59.93e − 281.53e − 22.80e − 33.90e − 33.95e − 33.71e − 41.65e − 24.25e − 20.5It can be seen from Tables3 and 4 that the orthogonality index between IMFs by the traditional EMD decomposition is up to 10−4, and the IMF orthogonal index by orthogonalization can reach 10−19, and the accuracy is improved by 15 orders of magnitude, which is less than 10−16, the effective magnitude that the computer calculates in orthogonal index. The result shows that the IMF obtained by using the orthogonalized algorithm has strict orthogonality.Table 5
Energy index of each component in OEMD1 and EMD.
MethodEXE1E2E3E4E5E6E7E8EIOTErrorEMD2.63e54.6e31.5e41.1e51.5e54.6e41.3e44.3e31.2e33.44e530.8%OEMD12.63e57.9e39.4e49.5e45.3e44.1e34.0e31.6e31.2e32.61e50.76%In addition, the energy index can also be used to analyze the orthogonality of IMF. Based on the total energy of EMD, OEMD1, and OEMD2 which can be drawn from formulas (2) and (3), and IMF-decomposed energy shown in Tables 5 and 6, it can be seen that the energy sum of IMF components obtained by EMD has a larger error than that of the original signal, which is 30.8%, and the energy sums of IMF components obtained by the decomposed algorithm of OEMD1 and OEMD2 have the errors of 0.76% and 1.14% against that of the original signal, so the application of OEMD1 and OEMD2 can effectively reduce the energy error, and the effect is obvious.Table 6
Energy index of each component in OEMD2 and EMD.
MethodEXE1E2E3E4E5E6E7E8EIOTErrorEMD2.63e54.3e31.3e44.6e41.5e51.1e51.5e44.6e31.2e33.44e530.8%OEMD22.63e52.6e34.2e32.6e41.1e51.0e51.2e44.3e31.2e32.60e51.14%It is known from Table5 that the energy of this signal is mainly concentrated in the first 5 orders, accounting for 96.05% of the total energy. Therefore, OEMD3 is used to process the signal in the orthogonal order of c5-c4-c3-c2-c1-c6-c7-c8. The decomposed IMFs are shown in Figure 13. The overall orthogonal index is 0.001121 drawn from formula (1)–IOT. The orthogonal index of each component is calculated according to formula (1)–IOTjk, as shown in Table 7; the upper triangular data of this table are the orthogonal indexes between each IMF component obtained by EMD decomposition, and the lower triangular data are the orthogonal indexes by OEMD3 decomposition. From the table, it can be seen that the application of OEMD3 algorithm improves the IMF orthogonal order of magnitude as well as the first two algorithms. This algorithm ensures the strict orthogonality of each IMF component that is decomposed.Figure 13
OEMD3 decomposes IMF components (no loose fasteners).Table 7
Orthogonality index of each component in EMD and OEMD3.
IMF1234567810.53.80e − 21.97e − 21.02e − 21.11e − 26.29e − 29.82e − 34.79e − 321.05e − 170.55.84e − 22.36e − 36.89e − 36.65e − 31.17e − 23.48e − 334.42e − 171.93e − 160.56.19e − 42.80e − 31.13e − 26.07e − 34.35e − 341.23e − 188.65e − 188.68e − 180.58.86e − 21.83e − 32.15e − 34.76e − 354.38e − 177.20e − 185.32e − 171.39e − 160.57.31e − 31.80e − 41.53e − 361.62e − 161.38e − 171.10e − 173.53e − 191.87e − 180.55.6e − 21.59e − 274.78e − 188.10e − 175.82e − 186.86e − 183.79e − 181.20e − 160.53.63e − 284.79e − 33.26e34.75e34.68e31.15e21.65e24.25e20.5According to equations (2) and (3), the energy values of IMFs obtained by the decomposition of OEMD3 algorithm are shown in Table 8. The table shows that the energy error between IMFs decomposed by this algorithm is 0.38%, which is smaller than that of the first two algorithms. The energy errors between the IMF components reduce the energy leakage effectively.Table 8
The energy values of IMFs decomposed by OEMD3 algorithm.
MethodEXE1E2E3E4E5E6E7E8EIOTError (%)EMD2.63e54.6e41.5e51.1e51.5e44.6e31.3e44.3e31.2e33.39e529.0OEMD32.63e52.3e41.1e59.5e41.1e44.5e31.5e44.2e31.2e32.64e50.38To sum up, from the calculation of the wheelset vertical acceleration signals without fastener loosening, it can be seen that the orthogonal index of IMFs decomposed by OEMD1, OEMD2, and OEMD3 has a great improvement, and the strict orthogonality between the components is guaranteed compared with the orthogonal index of the IMFs decomposed by EMD. According to the results of energy index and orthogonality index calculated by the three algorithms, the total energy index error of each IMF component decomposed by OEMD3 and the error of each component energy index are less than the other two algorithms. Therefore, this algorithm is used to decompose the vertical acceleration of the wheelset with different degrees of fasteners loosening in condition 1, condition 2, and condition 3. The wheelset vertical accelerations decomposed by EMD and OEMD3 with fasteners loosening 50% and 100% in condition 1 are shown in Figure14. It can be seen that, in condition 1, there are some certain differences in the IMF obtained from OEMD3 with the vertical acceleration signals that have different degrees of loose fasteners. In particular, there is one more orthogonal IMF component when the fastener is completely loosened than in the other loosening degrees. This shows that the different degrees of fastener loosening have different effects on the frequency components of the vertical acceleration signal; the difference is most obvious when the fastener is completely loosened, so the time-frequency analysis is considered in the follow-up analysis to identify the loosening features of the fastener.Figure 14
Decomposition results of EMD and OEMD3: (a, b) fastener is loosened by 50% in condition 1. (c, d) Fastener is loosened by 50% in condition 1. (a) EMD decomposition (50%). (b) OEMD3 decomposition (50%). (c) EMD decomposition (100%). (d) OEMD3 decomposition (100%).
(a)(b)(c)(d)


### 3.3. Orthogonal IMF Selection Based on Correlation

The above analysis shows that the OEMD can well guarantee the strict orthogonality between the IMF components and ensure that the energy does not leak, but it can be seen from Figure13 that some interfering signals are introduced by IMFs. These interfering signals have an impact on the accuracy of the subsequent analysis and should be eliminated. Orthogonalized IMF should have a good correlation with the original signal, so this paper uses the correlation coefficient to select a large correlation IMF. Setting the threshold to 0.5, when the correlation coefficient of IMF and the original signal is greater than 0.5, it is believed that the reliability of the signal is relatively higher. The correlation coefficients between the IMF components and the original signals decomposed by the OEMD3 in the first condition when the fastener is not loosed are shown in Table 9. It can be seen from the table that only IMF2 and IMF3 satisfy the conditions, and the rest of IMF components and residual components are discarded without any computing.Table 9
Correlation coefficient between IMF and original signal.
IMFCorrelation coefficient10.184220.566130.623540.388150.026460.100370.058080.083990.0123


### 3.4. Fastener Looseness Signal Analysis Based on Hilbert–Huang Transform

Hilbert transform is applied to the IMF signal processed by OEMD3 to obtain the corresponding Hilbert time spectrum, and the loose feature is extracted by quantitative analysis. Time-frequency spectrum can reflect the distribution of energy in time and frequency and can reflect the change of signal amplitude with time and frequency. However, due to the complexity of wheelset vertical acceleration signal components and Hilbert spectrum, it requires a lot of experience to directly identify the feature information from Hilbert spectrum. In the time and frequency plane, the difference of different degrees of loosening signals is reflected in the different energy distribution in the corresponding area, that is, the uniformity of energy distribution in the same area. Therefore, the information entropy theory is applied to the quantitative recognition of Hilbert time spectrum, which simplifies the complexity of the algorithm in the calculation speed and arduousness.


### 3.5. Fastener Looseness Feature Recognition Based on Energy Entropy

In the information theory, supposing the sample space of discrete random variableX is S = {x1, x1, …, xn} and the probability of random variable X = xj is Pj, then the information quantity is(13)Ixj=IX=xj=log1pj=−logpj.The amount of information describes the relationship between the probability of occurrence and information content. The mean value of information quantityI (xj) in S is the information entropy of X, written as(14)HX=−∑j=1Npjlogpj.In this paper, the information entropy and the energy distribution of Hilbert time-frequency spectrum are combined to calculate the entropy of fastener looseness in each condition. This method is called HHT energy entropy. The specific method is as follows.The energy of the entire Hilbert spectrum plane is denoted asA, which is equally divided into N equal areas, the energy in each small area is denoted as Wi(i = 1, 2, …, N), and the energy of each small area is normalized; thus, qi = Wi/A. According to the formula of information entropy, the formula of HHT energy entropy of wheel acceleration signal is(15)Eq=−∑i=1Nqilnqi.According to the basic properties of energy entropy of HHT, the value ofE (q) can reflect the uniformity of the energy in the time-frequency plane. The more uniform the energy distribution in the plane, the larger the value of E (q); otherwise, E (q) becomes smaller.According to the principle of correlation, the IMF components which are closely related to the original signal are selected for the Hilbert frequency spectrum analysis. Finally, the time-frequency spectrum is divided into 170 equal portions, and the entropy of HHT energy in time-frequency spectrum is calculated. The specific process is shown in Figure15.Figure 15
Calculation of the Hilbert energy entropy process. Firstly, varying vertical acceleration signals are gotten and then decomposed by OEMD3. Secondly, IMFS is selected based on correlation principle, which is calculated by Hilbert transform. Lastly, energy entropy theory is used to compute the energy in different looseness conditions.The values of HHT energy entropy and their trends under various conditions are shown in Table10 and Figure 16. It can be seen that the HHT energy entropy of the first condition (one loose fastener) is higher than that of condition 2 and condition 3 (two loose fasteners); in case of 100% looseness, the energy entropy of condition 1 is 4.2672, and those of condition 2 and condition 3 are 4.2043 and 4.2339, respectively. For the same working condition, the HHT energy entropy is less than that without looseness; for example, in working condition 1, the energy entropy decreases from 4.5678 to 4.2672 as looseness increases. Therefore, the HTT energy entropy of wheelset vibration signals decreases with the increase of fastener loosening degree. When the fastener is completely loose, the HHT energy entropy is the smallest. This is because the stiffness of the fastener affects the vibration of various components of the vehicle and rail system. As the rigidity decreases, the rail vibration and the rail displacement will change, just like higher irregularity’s impact on a vehicle in the running line, which will affect the vibration of the vehicle system, make the time-frequency distribution of the vibration signal inhomogeneous, and cause the value of energy entropy to drop. Therefore, calculating HHT energy entropy of wheelset vertical acceleration can reflect the looseness of fastener to a certain extent, and it can provide reference for the identification of fastener looseness.Table 10
Energy entropy of fastener loosening in various working conditions.
Level of fastener looseness0%25%50%75%100%Energy entropy (condition 1)4.56784.50324.46724.39864.2672Energy entropy (condition 2)4.56784.47914.36274.32764.2043Energy entropy (condition 3)4.56784.43564.39144.30224.2339Figure 16
Energy entropy trend of vertical vibration of wheelset. Condition 1 (black square), condition 2 (red dot), and condition 3 (blue triangle); the looseness is from 0% to 100%, increasing gradually at an interval of 25%.


## 3.1. Identification of Fasteners State in Varying Degrees

HHT has good reliability in dealing with nonstationary and nonlinear signals, detecting structural faults. However, it is not based on a complete theory and has energy leakage problems. In order to solve this problem, this paper optimizes the algorithm according to the theory of orthogonality and correlation and combines the improved Hilbert Huang algorithm with the theory of information entropy [22, 23], which is applied to the identification of rail fastener looseness.


### 3.1.1. Hilbert–Huang Transform

The crux of HHT algorithm is EMD decomposition. The result of EMD decomposition directly affects the accuracy of subsequent signal processing, so it is necessary to ensure that the IMF component of the EMD decomposition should possess completeness and orthogonality, so that no energy leakage in the decomposition can be guaranteed [24, 25]. From the practical point of view, Huang et al. consider that there exists orthogonality between all IMFs decomposed by EMD, but there is no rigorous theoretical derivation proving that they are rigorously orthogonal in the overall situation[29, 30].Huang proposed two orthogonality indicators to measure the orthogonality between IMF components, namely, the overall orthogonality index (IOT) and the orthogonality index between two components (IOjk) [31]. When all the IMF components are rigorously orthogonal to each other, the IOT and IOjk should be zero.(1)IOT=∑j=1n+1∑k=1k≠jn+1∑i=1Nimfjiimfki∑i=1Nxi2,IOjk=∑i=1Nimfjiimfki∑i=1Nimfji2+imfki2.In the formula,xi is the ith original signal, imfji is the jth IMF component of xi, imfk is the kth IMF component of xi, and j is not equal to k. In addition, the energy can be used to measure the degree of orthogonality between IMFs. The energy of the original signal X(t) = [x1, x2, …, xn] is(2)Ex=∫0TX2tdt=∑i=1Nxi2.The energy of each IMF is(3)Ej=∫0Timfj2tdt=∑i=1Nimfji2t,j=1,…,n+1.In the formula, imfj (t) is the jth IMF component; if all the IMFs are strictly orthogonal to each other, Ex = Ej and the leaking energy between the IMFs is zero, that is,(4)EIOT=∑j=1n+1Ejt=Ex,where EIOT is the energy of IOT, total energy after decomposition which is constant.(5)Ejk=∫0Timfjtimfktdt=∑i=1Nimfjiimfki=0,j,k=1,…,n+1;j≠k.Ejk is the leaking energy between imfj and imfk. From the processing steps of EMD decomposition, we can see that, in the actual decomposition process, the mean value curve is obtained by fitting the approximate value, and the decomposed IMF is a part of the original signal, so the different IMF components are approximately orthogonal [29, 30]. Huang’s theory proves that the EMD orthogonality is in the order of 10−3 to 10−2, and the IMFs are not strictly orthogonal, resulting in the phenomenon of energy leakage and modal aliasing during signal analysis, which will bring errors in the later analysis, making the result inaccurate. In order to eliminate energy leakage, it is necessary to ensure that the orthonormal mode functions are strictly orthogonal to each other [31]. Therefore, the traditional EMD algorithm is improved.


### 3.1.2. OEMD Algorithm

In order to ensure the strict orthogonality between the IMF components after EMD decomposition [32, 33], this paper orthogonalizes the IMF components; the basic process is shown as follows:(1)
The original signalX(t) is decomposed into the form of the sum of multiple IMF components ci(t) (i = 1, 2, …, n) and a residual component rn (t) through EMD:(6)Xt=∑j=1ncjt+rnt.(2)
Assume imf1 (t) = c1 (t) and it is the first orthogonalized IMF component of the original signal, which is the highest frequency component of the IMF.(3)
In the process of decomposition to obtain the second IMF componentc2 (t), we can see that it is impossible to ensure the orthogonality between c2 (t) and imf1 (t). Therefore, we need to remove imf1 (t) from c2 (t); that is,(7)imf2t=c2t−β21imf1t,imf2 (t) is the second orthogonalized IMF component of the original signal X (t), and β21 is the orthogonalization coefficient between c2 (t) and imf1 (t). Each side of equations (2)–(7) is multiplied by imf1 (t) so as to integrate t. Since imf1 (t) and imf2 (t) are orthogonal, it can be concluded that(8)∫0Timf1timf2tdt=∫0Tc2timf1tdt−β21∫0Timf12tdt=0,(9)β21=c2⟶Timf1⟶imf1⟶Timf1⟶.The numerator and denominator in equations (2)–(9) are inner product operations of two vectors.(4)
By analogy, thej + 1th orthogonal component, imfj+1 (t), of the original signal can be obtained by eliminating the imfi (t) (i = 1, 2, …, j) components in the j + 1th IMF component cj+1 (t). That is,(10)imfj+1t=cj+1t−∑i=1jβj+1,iimfit.Each side of the above equation is multiplied by imfi (t) (i ≤ j), so as to integrate t. Owing to the fact that imfi (t) and imfj+1 (t) are orthogonal, it can be concluded that(11)βj+1,i=cj+1⟶Timfi⟶imfi⟶Timfi⟶.After the above calculation, the original signalX (t) is decomposed into the form of the sum of the orthogonalized IMF components and the residual components, that is,(12)Xt=c1t+c2t+⋯+cnt+rnt=imf1t+imf2t+β21imf1t+imf3t+β31imf1t+β32c2t+⋯+imfnt+βn1imf1t+βn2imf2t+⋯+βn,n−1imfn−1t+rnt=1+β21+β31+β41+⋯+βn1imf1t+1+β32+β42+⋯+βn2imf2t+⋯+1+βn,n−1imfn−1t+imfnt+rnt=imf1∗t+imf2∗t+⋯+imfn∗t+rnt=∑j=1nimfj∗t+rnt=∑j=1najimfjt+rnt.In equations (2)–(12), aj=∑i=jnβi,j,j=1,2,…,n,βi,j=1i=j. From the above calculation, we can see that there is a strict orthogonality in imfj (t), so the linear transformation on imfj (t) will not change its orthogonality, and imfj∗t is also strictly orthogonal.The integration of the intrinsic mode components after the decomposition of the above EMD is called the orthogonal empirical mode decomposition (OEMD). According to the different orthogonal sequence of IMF components, the OEMD can be divided into three kinds. The way to orthogonalize the high frequency IMF component before orthogonalizing the low frequency IMF component is called OEMD1, the way to deal with IMF in reverse order is called OEMD2, and the way to perform orthogonalization from any IMF is called OEMD3. The essence of OEMD is still EMD; only the IMF components decomposed by EMD are orthogonalized and then reorganized to achieve strict orthogonality of each IMF component.


## 3.1.1. Hilbert–Huang Transform

The crux of HHT algorithm is EMD decomposition. The result of EMD decomposition directly affects the accuracy of subsequent signal processing, so it is necessary to ensure that the IMF component of the EMD decomposition should possess completeness and orthogonality, so that no energy leakage in the decomposition can be guaranteed [24, 25]. From the practical point of view, Huang et al. consider that there exists orthogonality between all IMFs decomposed by EMD, but there is no rigorous theoretical derivation proving that they are rigorously orthogonal in the overall situation[29, 30].Huang proposed two orthogonality indicators to measure the orthogonality between IMF components, namely, the overall orthogonality index (IOT) and the orthogonality index between two components (IOjk) [31]. When all the IMF components are rigorously orthogonal to each other, the IOT and IOjk should be zero.(1)IOT=∑j=1n+1∑k=1k≠jn+1∑i=1Nimfjiimfki∑i=1Nxi2,IOjk=∑i=1Nimfjiimfki∑i=1Nimfji2+imfki2.In the formula,xi is the ith original signal, imfji is the jth IMF component of xi, imfk is the kth IMF component of xi, and j is not equal to k. In addition, the energy can be used to measure the degree of orthogonality between IMFs. The energy of the original signal X(t) = [x1, x2, …, xn] is(2)Ex=∫0TX2tdt=∑i=1Nxi2.The energy of each IMF is(3)Ej=∫0Timfj2tdt=∑i=1Nimfji2t,j=1,…,n+1.In the formula, imfj (t) is the jth IMF component; if all the IMFs are strictly orthogonal to each other, Ex = Ej and the leaking energy between the IMFs is zero, that is,(4)EIOT=∑j=1n+1Ejt=Ex,where EIOT is the energy of IOT, total energy after decomposition which is constant.(5)Ejk=∫0Timfjtimfktdt=∑i=1Nimfjiimfki=0,j,k=1,…,n+1;j≠k.Ejk is the leaking energy between imfj and imfk. From the processing steps of EMD decomposition, we can see that, in the actual decomposition process, the mean value curve is obtained by fitting the approximate value, and the decomposed IMF is a part of the original signal, so the different IMF components are approximately orthogonal [29, 30]. Huang’s theory proves that the EMD orthogonality is in the order of 10−3 to 10−2, and the IMFs are not strictly orthogonal, resulting in the phenomenon of energy leakage and modal aliasing during signal analysis, which will bring errors in the later analysis, making the result inaccurate. In order to eliminate energy leakage, it is necessary to ensure that the orthonormal mode functions are strictly orthogonal to each other [31]. Therefore, the traditional EMD algorithm is improved.


## 3.1.2. OEMD Algorithm

In order to ensure the strict orthogonality between the IMF components after EMD decomposition [32, 33], this paper orthogonalizes the IMF components; the basic process is shown as follows:(1)
The original signalX(t) is decomposed into the form of the sum of multiple IMF components ci(t) (i = 1, 2, …, n) and a residual component rn (t) through EMD:(6)Xt=∑j=1ncjt+rnt.(2)
Assume imf1 (t) = c1 (t) and it is the first orthogonalized IMF component of the original signal, which is the highest frequency component of the IMF.(3)
In the process of decomposition to obtain the second IMF componentc2 (t), we can see that it is impossible to ensure the orthogonality between c2 (t) and imf1 (t). Therefore, we need to remove imf1 (t) from c2 (t); that is,(7)imf2t=c2t−β21imf1t,imf2 (t) is the second orthogonalized IMF component of the original signal X (t), and β21 is the orthogonalization coefficient between c2 (t) and imf1 (t). Each side of equations (2)–(7) is multiplied by imf1 (t) so as to integrate t. Since imf1 (t) and imf2 (t) are orthogonal, it can be concluded that(8)∫0Timf1timf2tdt=∫0Tc2timf1tdt−β21∫0Timf12tdt=0,(9)β21=c2⟶Timf1⟶imf1⟶Timf1⟶.The numerator and denominator in equations (2)–(9) are inner product operations of two vectors.(4)
By analogy, thej + 1th orthogonal component, imfj+1 (t), of the original signal can be obtained by eliminating the imfi (t) (i = 1, 2, …, j) components in the j + 1th IMF component cj+1 (t). That is,(10)imfj+1t=cj+1t−∑i=1jβj+1,iimfit.Each side of the above equation is multiplied by imfi (t) (i ≤ j), so as to integrate t. Owing to the fact that imfi (t) and imfj+1 (t) are orthogonal, it can be concluded that(11)βj+1,i=cj+1⟶Timfi⟶imfi⟶Timfi⟶.After the above calculation, the original signalX (t) is decomposed into the form of the sum of the orthogonalized IMF components and the residual components, that is,(12)Xt=c1t+c2t+⋯+cnt+rnt=imf1t+imf2t+β21imf1t+imf3t+β31imf1t+β32c2t+⋯+imfnt+βn1imf1t+βn2imf2t+⋯+βn,n−1imfn−1t+rnt=1+β21+β31+β41+⋯+βn1imf1t+1+β32+β42+⋯+βn2imf2t+⋯+1+βn,n−1imfn−1t+imfnt+rnt=imf1∗t+imf2∗t+⋯+imfn∗t+rnt=∑j=1nimfj∗t+rnt=∑j=1najimfjt+rnt.In equations (2)–(12), aj=∑i=jnβi,j,j=1,2,…,n,βi,j=1i=j. From the above calculation, we can see that there is a strict orthogonality in imfj (t), so the linear transformation on imfj (t) will not change its orthogonality, and imfj∗t is also strictly orthogonal.The integration of the intrinsic mode components after the decomposition of the above EMD is called the orthogonal empirical mode decomposition (OEMD). According to the different orthogonal sequence of IMF components, the OEMD can be divided into three kinds. The way to orthogonalize the high frequency IMF component before orthogonalizing the low frequency IMF component is called OEMD1, the way to deal with IMF in reverse order is called OEMD2, and the way to perform orthogonalization from any IMF is called OEMD3. The essence of OEMD is still EMD; only the IMF components decomposed by EMD are orthogonalized and then reorganized to achieve strict orthogonality of each IMF component.


## 3.2. OEMD Decomposition of the Vertical Acceleration Signals of Wheelset

The EMD is used to decompose the vertical acceleration signal (Figure9) of the wheelset without loosening of the fastener under the working condition 1 and obtain 7 IMF components and 1 residual component, as shown in Figure 12. The same signals are processed by OEMD1 and OEMD2, respectively; the results of 7 IMF components and 1 residual component can be obtained.Figure 12
EMD decomposition results of simulation signal (no loose fasteners).The overall orthogonal index obtained by the decomposition of EMD is 0.7546, which can be calculated by formula (1). The overall orthogonal indices obtained by OEMD1 and OEMD2 are 0.006463 and 0.007347, respectively. It can be seen that the overall orthogonal index obtained by the OEMD is 2 orders of magnitude higher than that of the overall EMD. There is little difference in accuracy of the orthogonality of the two orthogonal EMD algorithms. Since the 8th IMF is a residual component with a very small value, it has no orthogonality with the first 7 components.The orthogonality between the 7 IMF components after the decomposition of EMD, OEMD1, and OEMD2 is calculated, respectively, by formula (1); the results of the calculation are shown in Tables 3 and 4. The upper triangular data of Table 3 are the orthogonal indexes between IMF components obtained by EMD decomposition, and the lower triangular data are the orthogonal indexes between IMF components obtained by OEMD1 decomposition. The upper triangular data of Table 4 are the orthogonal indexes obtained by EMD decomposition, and the lower triangular data are the orthogonal indexes obtained by OEMD2 decomposition.Table 3
Orthogonal index of each component in EMD and OEMD2.
IMF1234567810.55.60e − 29.82e − 31.17e − 26.07e − 32.16e − 31.80e − 23.63e − 225.68e − 190.56.29e − 26.65e − 31.13e − 21.83e − 27.31e − 21.59e − 232.45e − 178.16e − 170.53.80e − 21.92e − 21.02e − 21.11e − 24.79e − 342.16e − 179.92e − 181.48e − 160.55.84e − 22.36e − 36.89e − 33.48e − 353.31e − 182.11e − 171.91e − 171.22e − 160.56.19e − 42.80e − 34.35e − 365.56e − 182.65e − 181.75e − 174.80e − 181.18e − 170.58.86e − 24.76e − 374.90e − 187.17e − 181.78e − 179.65e − 183.10e − 171.93e − 160.51.53e − 283.63e − 21.94e − 27.18e − 33.62e − 32.73e − 44.38e − 21.08e − 20.5Table 4
Orthogonal index of each component in EMD and OEMD1.
IMF1234567810.58.86e − 22.80e − 36.89e − 31.11e − 27.31e − 21.80e − 21.70e − 224.76e − 170.56.19e − 42.36e − 31.02e − 21.83e − 22.16e − 31.07e − 233.15e − 182.76e − 170.55.84e − 21.92e − 21.13e − 26.07e − 31.03e − 343.10e − 171.07e − 174.52e − 170.53.80e − 26.65e − 31.17e − 23.37e − 452.27e − 171.17e − 173.00e − 178.13e − 180.56.29e − 29.82e − 31.89e − 361.64e − 181.24e − 173.46e − 171.37e − 172.66e − 160.55.60e − 21.81e − 273.87e − 201.97e − 182.60e − 184.42e − 175.67e − 171.23e − 160.59.93e − 281.53e − 22.80e − 33.90e − 33.95e − 33.71e − 41.65e − 24.25e − 20.5It can be seen from Tables3 and 4 that the orthogonality index between IMFs by the traditional EMD decomposition is up to 10−4, and the IMF orthogonal index by orthogonalization can reach 10−19, and the accuracy is improved by 15 orders of magnitude, which is less than 10−16, the effective magnitude that the computer calculates in orthogonal index. The result shows that the IMF obtained by using the orthogonalized algorithm has strict orthogonality.Table 5
Energy index of each component in OEMD1 and EMD.
MethodEXE1E2E3E4E5E6E7E8EIOTErrorEMD2.63e54.6e31.5e41.1e51.5e54.6e41.3e44.3e31.2e33.44e530.8%OEMD12.63e57.9e39.4e49.5e45.3e44.1e34.0e31.6e31.2e32.61e50.76%In addition, the energy index can also be used to analyze the orthogonality of IMF. Based on the total energy of EMD, OEMD1, and OEMD2 which can be drawn from formulas (2) and (3), and IMF-decomposed energy shown in Tables 5 and 6, it can be seen that the energy sum of IMF components obtained by EMD has a larger error than that of the original signal, which is 30.8%, and the energy sums of IMF components obtained by the decomposed algorithm of OEMD1 and OEMD2 have the errors of 0.76% and 1.14% against that of the original signal, so the application of OEMD1 and OEMD2 can effectively reduce the energy error, and the effect is obvious.Table 6
Energy index of each component in OEMD2 and EMD.
MethodEXE1E2E3E4E5E6E7E8EIOTErrorEMD2.63e54.3e31.3e44.6e41.5e51.1e51.5e44.6e31.2e33.44e530.8%OEMD22.63e52.6e34.2e32.6e41.1e51.0e51.2e44.3e31.2e32.60e51.14%It is known from Table5 that the energy of this signal is mainly concentrated in the first 5 orders, accounting for 96.05% of the total energy. Therefore, OEMD3 is used to process the signal in the orthogonal order of c5-c4-c3-c2-c1-c6-c7-c8. The decomposed IMFs are shown in Figure 13. The overall orthogonal index is 0.001121 drawn from formula (1)–IOT. The orthogonal index of each component is calculated according to formula (1)–IOTjk, as shown in Table 7; the upper triangular data of this table are the orthogonal indexes between each IMF component obtained by EMD decomposition, and the lower triangular data are the orthogonal indexes by OEMD3 decomposition. From the table, it can be seen that the application of OEMD3 algorithm improves the IMF orthogonal order of magnitude as well as the first two algorithms. This algorithm ensures the strict orthogonality of each IMF component that is decomposed.Figure 13
OEMD3 decomposes IMF components (no loose fasteners).Table 7
Orthogonality index of each component in EMD and OEMD3.
IMF1234567810.53.80e − 21.97e − 21.02e − 21.11e − 26.29e − 29.82e − 34.79e − 321.05e − 170.55.84e − 22.36e − 36.89e − 36.65e − 31.17e − 23.48e − 334.42e − 171.93e − 160.56.19e − 42.80e − 31.13e − 26.07e − 34.35e − 341.23e − 188.65e − 188.68e − 180.58.86e − 21.83e − 32.15e − 34.76e − 354.38e − 177.20e − 185.32e − 171.39e − 160.57.31e − 31.80e − 41.53e − 361.62e − 161.38e − 171.10e − 173.53e − 191.87e − 180.55.6e − 21.59e − 274.78e − 188.10e − 175.82e − 186.86e − 183.79e − 181.20e − 160.53.63e − 284.79e − 33.26e34.75e34.68e31.15e21.65e24.25e20.5According to equations (2) and (3), the energy values of IMFs obtained by the decomposition of OEMD3 algorithm are shown in Table 8. The table shows that the energy error between IMFs decomposed by this algorithm is 0.38%, which is smaller than that of the first two algorithms. The energy errors between the IMF components reduce the energy leakage effectively.Table 8
The energy values of IMFs decomposed by OEMD3 algorithm.
MethodEXE1E2E3E4E5E6E7E8EIOTError (%)EMD2.63e54.6e41.5e51.1e51.5e44.6e31.3e44.3e31.2e33.39e529.0OEMD32.63e52.3e41.1e59.5e41.1e44.5e31.5e44.2e31.2e32.64e50.38To sum up, from the calculation of the wheelset vertical acceleration signals without fastener loosening, it can be seen that the orthogonal index of IMFs decomposed by OEMD1, OEMD2, and OEMD3 has a great improvement, and the strict orthogonality between the components is guaranteed compared with the orthogonal index of the IMFs decomposed by EMD. According to the results of energy index and orthogonality index calculated by the three algorithms, the total energy index error of each IMF component decomposed by OEMD3 and the error of each component energy index are less than the other two algorithms. Therefore, this algorithm is used to decompose the vertical acceleration of the wheelset with different degrees of fasteners loosening in condition 1, condition 2, and condition 3. The wheelset vertical accelerations decomposed by EMD and OEMD3 with fasteners loosening 50% and 100% in condition 1 are shown in Figure14. It can be seen that, in condition 1, there are some certain differences in the IMF obtained from OEMD3 with the vertical acceleration signals that have different degrees of loose fasteners. In particular, there is one more orthogonal IMF component when the fastener is completely loosened than in the other loosening degrees. This shows that the different degrees of fastener loosening have different effects on the frequency components of the vertical acceleration signal; the difference is most obvious when the fastener is completely loosened, so the time-frequency analysis is considered in the follow-up analysis to identify the loosening features of the fastener.Figure 14
Decomposition results of EMD and OEMD3: (a, b) fastener is loosened by 50% in condition 1. (c, d) Fastener is loosened by 50% in condition 1. (a) EMD decomposition (50%). (b) OEMD3 decomposition (50%). (c) EMD decomposition (100%). (d) OEMD3 decomposition (100%).
(a)(b)(c)(d)


## 3.3. Orthogonal IMF Selection Based on Correlation

The above analysis shows that the OEMD can well guarantee the strict orthogonality between the IMF components and ensure that the energy does not leak, but it can be seen from Figure13 that some interfering signals are introduced by IMFs. These interfering signals have an impact on the accuracy of the subsequent analysis and should be eliminated. Orthogonalized IMF should have a good correlation with the original signal, so this paper uses the correlation coefficient to select a large correlation IMF. Setting the threshold to 0.5, when the correlation coefficient of IMF and the original signal is greater than 0.5, it is believed that the reliability of the signal is relatively higher. The correlation coefficients between the IMF components and the original signals decomposed by the OEMD3 in the first condition when the fastener is not loosed are shown in Table 9. It can be seen from the table that only IMF2 and IMF3 satisfy the conditions, and the rest of IMF components and residual components are discarded without any computing.Table 9
Correlation coefficient between IMF and original signal.
IMFCorrelation coefficient10.184220.566130.623540.388150.026460.100370.058080.083990.0123


## 3.4. Fastener Looseness Signal Analysis Based on Hilbert–Huang Transform

Hilbert transform is applied to the IMF signal processed by OEMD3 to obtain the corresponding Hilbert time spectrum, and the loose feature is extracted by quantitative analysis. Time-frequency spectrum can reflect the distribution of energy in time and frequency and can reflect the change of signal amplitude with time and frequency. However, due to the complexity of wheelset vertical acceleration signal components and Hilbert spectrum, it requires a lot of experience to directly identify the feature information from Hilbert spectrum. In the time and frequency plane, the difference of different degrees of loosening signals is reflected in the different energy distribution in the corresponding area, that is, the uniformity of energy distribution in the same area. Therefore, the information entropy theory is applied to the quantitative recognition of Hilbert time spectrum, which simplifies the complexity of the algorithm in the calculation speed and arduousness.


## 3.5. Fastener Looseness Feature Recognition Based on Energy Entropy

In the information theory, supposing the sample space of discrete random variableX is S = {x1, x1, …, xn} and the probability of random variable X = xj is Pj, then the information quantity is(13)Ixj=IX=xj=log1pj=−logpj.The amount of information describes the relationship between the probability of occurrence and information content. The mean value of information quantityI (xj) in S is the information entropy of X, written as(14)HX=−∑j=1Npjlogpj.In this paper, the information entropy and the energy distribution of Hilbert time-frequency spectrum are combined to calculate the entropy of fastener looseness in each condition. This method is called HHT energy entropy. The specific method is as follows.The energy of the entire Hilbert spectrum plane is denoted asA, which is equally divided into N equal areas, the energy in each small area is denoted as Wi(i = 1, 2, …, N), and the energy of each small area is normalized; thus, qi = Wi/A. According to the formula of information entropy, the formula of HHT energy entropy of wheel acceleration signal is(15)Eq=−∑i=1Nqilnqi.According to the basic properties of energy entropy of HHT, the value ofE (q) can reflect the uniformity of the energy in the time-frequency plane. The more uniform the energy distribution in the plane, the larger the value of E (q); otherwise, E (q) becomes smaller.According to the principle of correlation, the IMF components which are closely related to the original signal are selected for the Hilbert frequency spectrum analysis. Finally, the time-frequency spectrum is divided into 170 equal portions, and the entropy of HHT energy in time-frequency spectrum is calculated. The specific process is shown in Figure15.Figure 15
Calculation of the Hilbert energy entropy process. Firstly, varying vertical acceleration signals are gotten and then decomposed by OEMD3. Secondly, IMFS is selected based on correlation principle, which is calculated by Hilbert transform. Lastly, energy entropy theory is used to compute the energy in different looseness conditions.The values of HHT energy entropy and their trends under various conditions are shown in Table10 and Figure 16. It can be seen that the HHT energy entropy of the first condition (one loose fastener) is higher than that of condition 2 and condition 3 (two loose fasteners); in case of 100% looseness, the energy entropy of condition 1 is 4.2672, and those of condition 2 and condition 3 are 4.2043 and 4.2339, respectively. For the same working condition, the HHT energy entropy is less than that without looseness; for example, in working condition 1, the energy entropy decreases from 4.5678 to 4.2672 as looseness increases. Therefore, the HTT energy entropy of wheelset vibration signals decreases with the increase of fastener loosening degree. When the fastener is completely loose, the HHT energy entropy is the smallest. This is because the stiffness of the fastener affects the vibration of various components of the vehicle and rail system. As the rigidity decreases, the rail vibration and the rail displacement will change, just like higher irregularity’s impact on a vehicle in the running line, which will affect the vibration of the vehicle system, make the time-frequency distribution of the vibration signal inhomogeneous, and cause the value of energy entropy to drop. Therefore, calculating HHT energy entropy of wheelset vertical acceleration can reflect the looseness of fastener to a certain extent, and it can provide reference for the identification of fastener looseness.Table 10
Energy entropy of fastener loosening in various working conditions.
Level of fastener looseness0%25%50%75%100%Energy entropy (condition 1)4.56784.50324.46724.39864.2672Energy entropy (condition 2)4.56784.47914.36274.32764.2043Energy entropy (condition 3)4.56784.43564.39144.30224.2339Figure 16
Energy entropy trend of vertical vibration of wheelset. Condition 1 (black square), condition 2 (red dot), and condition 3 (blue triangle); the looseness is from 0% to 100%, increasing gradually at an interval of 25%.


## 4. Conclusions

This paper starts with the idea of identifying the loosening of the fastener from the vibration response signal of the vehicle, and the traditional vehicle-track rigid-body model is improved. The vehicle-track rigid-flexible coupling model is established by using multibody dynamics software SIMPACK and finite element analysis software ANSYS. The rigid-flexible coupling model simulates the vibration response signals of the wheelsets and the vehicle body in three operating conditions. Combined with improved HHT algorithm and the information entropy theory, the energy entropy is applied to the recognition of the vertical acceleration signal of wheelset with fastener loosening; a certain degree of recognition effects had been achieved. The main conclusions of this paper are listed as follows:(1)
Firstly, the finite element model of rail and rail slab is established in ANSYS, and the model is discretized. The substructure modal analysis of rail and rail slab is carried out, and the modeling parameters are obtained. In order to solve the shortcoming that SIMPACK wheel rail module can only realize the physical contact between rigid body and rigid body, a virtual rail without mass and moment of inertia is added between the wheel and the rail, and the transfer of force is realized based on the balance of force and deformation compatibility condition, the fastener and CA mortar are simulated by force element, and the flexible track and rigid vehicle are assembled in the SIMPACK, and a vehicle-track rigid-flexible coupling model is established. Finally, the vertical accuracy of the model is verified by comparing the measured acceleration of quasi-high speed vehicles with different track irregularities.(2)
Three typical lines with loosened fasteners are set up and the five degrees of fastener loosening are set up. Based on the model built in this paper, the vehicle runs on the rails that the fasteners are loosened in varying degrees and under the effect of German low-interference spectrum at the speed of 200 km/h, and the corresponding vertical acceleration of wheelset and the vibration response of vehicle vertical acceleration are obtained. Finally, the vertical acceleration of the vehicle wheel with obvious vibration response characteristics is selected for the subsequent analysis.(3)
The Hilbert–Huang Transform is applied to the analysis of wheelset vibration signals. In order to solve the problems of low orthogonality and energy leakage in the first step EMD of HHT, the EMD is replaced by the orthogonal empirical mode decomposition in processing the vertical acceleration signal of wheelset. Based on this, the correlation principle is applied to eliminate the components whose correlation with the original signal in the decomposition results is weak. Owing to the fact that Hilbert time-frequency spectrum of vertical acceleration signals under different working conditions is complicated, the algorithm of Hilbert energy-spectrum entropy is proposed by combining the theory of HHT time-frequency spectrum with the theory of information entropy and is applied to the analysis and identification of wheelset vertical acceleration signals in different conditions. The results show that the method of energy entropy can be used to analyze the vertical vibration signals of the wheelset under different degrees of fastener looseness, and it can realize the identification of the looseness of rail fasteners. The simulation results show that the HHT energy entropy of the wheelset vertical vibration signals decreases with the increase of the degree of fastener loosening; the algorithm can recognize the loose state of rail fastener.



---

*Source: 4518624-2020-09-19.xml*