Abstract

The oscillatory behavior of the center of mass (CoM) and the ground reaction forces (GRFs) of walking people can be successfully explained by a 2D spring-loaded inverted pendulum (SLIP) model. However, the application of the 2D model is just restricted to a two-dimensional plane as the model fails to take the GRFs in the lateral direction into consideration. In this article, we simulated the gait cycle with a nonlinear dynamic model—a three-dimensional bipedal walking model—that compensated for defects in the 2D model. An experiment was conducted to compare the simulation results with the experimental data, which revealed that the experimental data of the ground reaction forces were in good agreement with the results of numerical simulation. A correlation analysis was also conducted between several initial dynamic parameters of the model. Through an examination of the impact of 3D dynamic parameters on the peaks of GRFs in three directions, we found that the 3D parameters had a major effect on the lateral GRFs. These findings demonstrate that the characteristics of human walking can be interpreted from a simple spring-damper system.

1. Introduction

Many mathematical models for dynamic systems in the field of engineering science are formulated as nonlinear differential equations. Many researchers have applied quantity theory or other conservative theories to build differential equations to simulate and analyze nonlinear systems [13]. Human walking models derived mathematically are widely used in many engineering fields, such as gait analysis [411], human walking assistance [12, 13], motion stability control of robots [1418], and response prediction of civil structures induced by pedestrians [1921]. A variety of kinetic models, both passive and active ones, are developed from nonlinear differential equations to study the kinetic phenomenon in human walking. One of the typical passive models that deserves particular attention is the bipedal spring-loaded inverted pendulum (SLIP) model, in which the legs are equivalent to springs, and the GRF is generated through the alternation of the front and rear legs [9, 2224]. The control strategy to guarantee the stability of walking robots based on the SLIP model has been widely studied and successfully applied [2531]. Later, a damper was added by Kim et al. to the compliant leg in the classic walking model for a better match with the experimental data, and the result showed that the stiffness and damping ratio exhibited a linear relationship with walking speed [11]. Although the classic bipedal spring-mass model gives a good prediction of the trajectory of the centroid and the complex features of the ground reaction force (GRF), the problem lies in that the previous model overestimates the GRF and underestimates the contact time [32]. Moreover, the SILP model was confirmed to demonstrate the regularity of human walking dynamics concerning age [33]. Lim and Park [34] extended the classic SLIP model to a new one which considered the asymmetry of the ankle joints and proved that the dynamics of multijoint lower limbs during the walking process could be represented by mechanical systems. Li et al. [35] proposed an actuated dissipative spring-mass model that mimicked the actual pattern of human walking. Hao et al. [36] applied a hip torque actuator to the model proposed by Alexander [7], the main advantage of which resides in its reduced energy consumption and a further improvement of the stability. Kuo [37] extended the classic planar motions to a three-dimensional motion which made it possible for the model to tilt side to side, and several strategies were implemented to strengthen the stability of the motion. Roos and Dingwell [38] introduced a dynamic walking model that incorporated a lateral step controller to maintain the lateral stability. They also found that there was a nonlinear relationship between neuromuscular noise and the probability of the model falling. With the undamped inverted pendulum modeled by the Duffing equation as an example, Selvam et al. [39] obtained the existence, Hyers–Ulam stability, and Hyers–Ulam Mittag–Leffler stability of solutions for a discrete fractional Duffing equation with forcing term, and stability conditions for the pendulum were obtained.

A major limitation of current studies about the 2D walking model is that all models are restricted to a two-dimensional plane and they fail to include the 3D dynamic parameters [26, 40, 41]. From another perspective, the 2D model is defective as to its inability to calculate the lateral GRF in terms that is the cardinal feature of the oscillation behavior of human walking.

dTo compensate for the 2D model’s inability to incorporate lateral GRFs, we present a 3D walking model that embodies GRFs in three directions. The contributions of this paper are as follows: (1) we extended the classic two-dimensional model proposed by Kim and Park [11] to a three-dimensional model that enables a more comprehensive consideration of the GRF of human walking in three directions. (2) We proposed an approach to develop dynamic equations for a 3D walking model through the Lagrangian method. The simulation results showed that our proposed model agreed well with the experimental results. (3) We conducted a correlation analysis between the parameters of the 3D model and confirmed a correlation between some parameters. (4) An investigation was undertaken to study the influence of kinetic parameters on the peak value of GRF.

This paper is organized as follows. Section 2 gives the kinetic formulation of the 3D walking model and provides a detailed description of the numerical simulation and the experiment protocols. In Section 3, the simulation, experimental results, and numerical analysis are displayed. In Section 4, we discuss the findings and the limitations of our study. Finally, we draw the conclusions in Section 5.

2. Materials and Methods

2.1. Bipedal 3D Inverted Pendulum Model

To describe the ground reaction forces (GRFs) during the human walking, we proposed a 3D inverted pendulum dynamic model based on the previous 2D model proposed by Kim and Park [11]. In the 3D model, the GRFs exerted by the foot against the ground are obtained from the swing of the center of mass (CoM) in a gait cycle. As shown in Figure 1, we introduce the angles and between limbs and the plumb line passing through the center of mass (CoM) to assess the vertical centroid oscillation in a gait cycle. The model consists of two legs with the same stiffness and damping , and the lengths of the leading and trailing leg are represented by and , respectively. and denote the angles between the projection of the two springs on the ground and the forward direction, respectively, to evaluate the lateral motion of the model during the walking process. For simplicity, we refer to the two angles and as the 3D angles and and as the 2D angles. It is worth noting that the direction of forward movement is self-defined and is fixed over a gait cycle which determines the initial angles and relevant parameters of the model. To construct the 3D model and demonstrate its lateral excursion, the semicircle curve of the classic 2D model is replaced by the hemispheric surface, with GRFs in three coordinate directions taken into account. The hemispheric surface with radius R = 0.3 m is used to simulate the foot (Figure 1). To ensure the transmission of force and moments from the CoM to the foot, the lower end of the spring and the spherical center of the hemispherical surface are rigidly tied. The upper end of the spring is hinged to the CoM, allowing the spring to rotate around the vertical axis freely in the three-dimensional space, with the vertical axis passing through the center of mass.

According to the Lagrange formulation (Appendix A) and space velocity relation (Appendix B), the motion equation of the 3D model for the single-support phase can be obtained as follows:

And the motion equation of the 3D model for the double-support phase can be obtained as follows:where , and represent the model mass, the radius of the foot hemisphere, 2D spring angle, spring length, the 3D angles, the 2D angular velocity, the change rate of spring length, the 3D angular velocity, slack length of spring, spring stiffness, and damping constant, respectively. The subscripts define the trailing and leading leg (Figure 1). It is worth noting that the three mutually independent variables , and can determine the double-support phase of the model due to the existence of geometric and the dynamic constraint.

2.2. Simulation Method

As shown in Figure 1, the positive direction of the x-axis is set as the initial direction of walking. With a similar pattern to that of the 2D model in [11, 22], the gait cycle in the 3D model is initiated from a single-support phase in which the 6 initial state parameters need to be determined. The parameters are , and the details of the parameters have been described before. During the simulation for the integration of the equations, the 3D model is initially positioned in a state when the point mass is at its apex during a gait cycle. The forward limb is positioned in front of the trailing limb at a predetermined angle . Before the contact of the leading leg with the ground, the motion equations (equation (1)) of the single-support phase are used to determine the dynamical behaviors of the model. It is worth noting that the angle of the 3D model is not initially zero, but a small value approaches zero, which is different from that in the 2D model to suit the need for convergence of the differential equation. In the numerical simulation, it is found that a stable result could be obtained when is in the range of 0.02 rad to 0.04 rad. Meanwhile, the initial angle and the corresponding angular velocity are selected through the Genetic Algorithm (GA) which will be described below.

When the model is transformed from the single-support phase to the double-support phase, the following geometric constraint needs to be satisfied:

Especially when the leading leg is just in contact with the ground, the following condition needs to be satisfied:where represents the touch-down angle. The calculation of the dynamic states of the model in the double-support phase requires the closed kinematic constraints (equation (B.4) in Appendix B) to determine the parameters related to the limb spring velocity , 2D angular velocity , and 3D angular velocity . If the trailing leg is lifted, and the point mass is at its apex, the trailing leg will be reset in front of the body with the predetermined impact angle. The model is instantaneously rotated along the plumb line passing through the body center, and the 3D angle and angular velocity after rotation are identical in value to those before the rotation but in the opposite direction. The process corresponds to the phase of “instant switch of legs” in Figure 2.

During a gait cycle, it is assumed that the change in the CoM and gait is symmetrical, which is consistent with the symmetric features of human walking. Based on the symmetric features and the consideration of computational efficiency, half of the gait cycle of the 3D model is calculated which is expressed as the single-double-single process in Figure 2. By iterating the initial conditions to be consistent with the end conditions, a stable GRFs curve can be finally obtained. The vertical GRF deserves more attention than the transverse and longitudinal components of the GRF for the vertical GRF is assumed to better reflect the oscillation characteristics of human walking. To optimize the initial parameter values, the Genetic Algorithm (GA) Solver of MATLAB® is utilized to identify the optimal solution for a specified nonlinear system with the fitness function set as the root mean square error between the data of the simulation of vertical GRF and the experimental ones. To prevent blind selection of parameters from causing nonconvergence of the solution, the parameter range is obtained from trial calculations by reference to the results of [22], and the final parameter range is as shown in Table 1. The parameters of GA are given as follows: the population size is set to 100, the mutation rate to 0.01, and the maximum number of generations for each run to 50.

2.3. Experiment

To validate the model proposed for 3D walking, the experimental data were collected to compare with the simulation results. Thirteen young and healthy subjects (body mass of 67.4 ± 12.2 kg and height of 1.67 ± 0.12 m) with no gait disorders participated in the experiment. The subjects walked on a walkway equipped with three force plates with piezoelectric 3-component force sensors. The participants were instructed to keep pace with a metronome placed nearby, which provided a fixed frequency each time. During the experiment, each participant was guided to walk at the frequency of the metronome, three times at each frequency. To get the ground reaction force (GRF) of each participant, three six-dimensional force plates (JP6060, CN) were used to measure the real-time GRF at a sample rate of 100 Hz. To ensure that the participants would cross the force plate at a comfortable and stable pace, the pre- and post-transition parts of the walkway made of wood were placed before and after the force plate (Figure 3). An AHRS system (AH100B, CN) was used to measure the acceleration and direction angle of CoM during the walking process at a sample rate of 100 Hz, as shown in Figure 4. The AHRS system could accurately and stably perform real-time 3D dynamic measurements with a Kalman filter. The velocity of the CoM could be calculated by integrating the acceleration data collected by the AHRS system.

3. Results

3.1. The Experimental Results Fit Well with the Numerical Simulation of the 3D Walking Model

To minimize the differentiation of the GRF due to the variation of participants in mass, the GRF is standardized by dividing the ground reaction force by the bodyweight as an evaluation criterion (Figure 5). With the initial parameters optimized by the genetic algorithm (GA) method, the simulation results are consistent with the experimental results well, with the root mean square error of for vertical forces. As can be seen from Figures 5(a) and 5(c), for GRF in the vertical and longitudinal direction, the simulation fits the data well in peak and valley values. There is a minor difference between the observed and simulated data in the valley values of the lateral GRF, which can be identified from Figure 5(b).

3.2. Effect of Model Dynamics Parameters on GRF Shape

The effects of spring stiffness and gait frequency on the GRF simulation results are shown in Figure 6(a). As demonstrated in Figure 6(a), the greater stiffness and gait frequency bring higher peaks, shorter gait cycles, and a more distinct gap between the peak and valley values but have little impact on the lateral and longitudinal forces. However, there is a significant effect of the 3D dynamics parameters on lateral forces, which is indicated in Figure 6(b). The results in Figure 6(b) reveal that the lateral GRF exhibits an increase in peaks with the increase of 3D angles and a decrease of 3D angular velocity. The effect of other parameters on lateral forces will be discussed later in Figure 7.

3.3. Correlations between the Initial Dynamic Parameters of the 3D Model

To explore the potential relationship between various parameters of the model and features of human gait, the kinetic parameter is defined as , and is obtained by experiment (Figure 8). A significantly positive correlation (r = 0.793, , r = correlation coefficient) between walking velocity and gait frequency is observed (Figure 9(a)).

When the kinetic parameters are identified by the GA to achieve a better fit between the simulations and the experimental results, a significant positive linear correlation can be observed between gait frequency and spring stiffness (r = 0.597, ) (Figure 9(b)) and initial 2D angular velocity (r = 0.338, ) (Figure 9(c)). Meanwhile, there is a significant negative correlation between the impact angle and gait frequency (r = −0.469, ), as illustrated in Figure 9(d).

Figure 10(a) shows the impact angle versus the initial 2D angle, with a correlation coefficient of 0.394 (). In contrast, a negative correlation is observed between impact angle and initial 3D angle (r = −0.306, ) (Figure 10(b)) and impact angle and damping ratio (r = −0.372, ) (Figure 10(c)).

Figures 9 and 10 show the correlation coefficient of the linear fit results for some of the dynamic parameters. As shown in Figure 11, a scatter plot matrix is designed to detect possible correlations among the initial parameters of the 3D model optimized by the genetic algorithm.

A scatter plot matrix is an effective multivariate visualization tool that plots the correlation between all variables. The nondiagonal part of the ith row and jth column represents the scatter plot of the ith variable against the jth variable, which is intended to describe the correlation between these two variables. If the two variables are positively correlated, the dependent variable y increases as the independent variable x becomes larger. If the two variables are negatively correlated, the dependent variable y becomes smaller as the independent variable x becomes larger. If the two variables are uncorrelated, the dependent variable y does not change as the independent variable x changes. In Figure 11, the solid line represents a least square fit of the data to a straight line, and the positive or negative sign of its slope determines whether the two variables are positively or negatively correlated. The thick dashed line indicates a smooth nonparametric curve under a smoother loess (the loess span 0.75) and the fine dashed line indicates the boundaries of the 95% confidence interval of the smoothed curve. The box plots of the dynamic parameters of the 3D model are distributed diagonally. From Figure 11, it can be drawn that all parameters are correlated to some extent.

To illustrate the level of correlation between the optimized initial variables, the linear fit results of the 10 parameters are used to generate a heat map (Figure 12). Correlation analysis in the heat map (Figure 12) shows that walking velocity is positively correlated with gait frequency (r = 0.79, ), step length (r = 0.63, ), and leg stiffness (r = 0.33, ), but has a negative correlation with impact angle (r = −0.40, ). Among the initial parameters, only the impact angle and initial 3D angular velocity are correlated with the damping ratio (r = −0.37 and r = 0.39, ). The impact angle is correlated with initial variables, except 2D and 3D angular velocity and step length. It is worth noting that in addition to its positive correlation with velocity, gait frequency is also positively correlated with leg stiffness (r = 0.6, ) and initial 2D angular velocity (r = 0.34, ), and there is a relatively significant negative correlation with impact angle (r = −0.47, ). There is a weak positive correlation between initial 2D angle and initial 3D angle (r = 0.32, ) and a weak negative correlation between the initial 2D angular velocity and initial 3D angular velocity (r = −0.38, ).

3.4. Effect of Initial Dynamic Parameters on GRFs
3.4.1. Effect of Initial Dynamic Parameters on GRF in the Longitudinal Direction

Kinetics studies are conducted to investigate the effect of various initial conditions on GRF in the longitudinal, vertical, and lateral direction by setting the initial walking velocity as and initial 2D angular velocity as for estimation.

The impacts that initial 3D angular velocity has on GRF in the longitudinal direction are grouped here into two phases (Figure 13). Phase I corresponds to the case where initial 3D angular velocity and peak forward forces show a negative correlation. In Phase II, with the increasing 3D angular velocity, the GRF values first increase till reaching a peak value and then decrease, which is exhibited in Figure 13(a) (), Figure 13(b) (), and Figure 13(c) (). It is manifested that the peak values tend to increase with increasing initial 2D angle . The effects of different initial factors on GRF peaks in the longitudinal direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.

3.4.2. Effect of Initial Dynamic Parameters on GRF in the Vertical Direction

As Figure 14 shows, the initial 3D angular velocity has a similar impact on vertical forces. An increase in the 3D angular velocity results in the decrease of vertical forces when the initial angle is small in Figure 14(a) (), Figure 14(b) (), and Figure 14(c) (). In contrast, the initial 3D angular velocity has little effect on the peak values when the initial angle is large, as shown in Figure 14(b) () and Figure 14(c) (). However, the trend is different when , and the vertical values first increase and then decrease as grows when and remain basically unchanged when (Figure 14(a)). As evident from Figure 14, with a fixed initial 3D angle , the larger the , the less influence the initial 3D angular velocity will have on the vertical peak values. The effects of different initial factors on GRF peaks in the vertical direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.

3.4.3. Effect of Initial Dynamic Parameters on GRF in the Lateral Direction

Figure 7(a) reveals that an increase in the initial 3D angular velocity results in a decrease of peak values of GRF in the lateral direction while the decreasing initial shows an opposite effect. The initial 3D angle is expected to have a negligible impact on the GRF values, which is illustrated in Figure 7(b). The effects of different initial factors on GRF peaks in lateral direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.

4. Discussion

The application of the conventional SLIP model is restricted to a two-dimensional plane, and the model can only predict the trajectory of CoM and the corresponding values of GRF in the vertical and longitudinal directions [11, 22] in the human walking process. In the conventional model, a curved foot is employed to mimic the pressure center excursions [8, 11, 22, 42]. In this article, we introduce the spatial angle and the corresponding angular velocity to construct a 3D model, with the foot treated as a hemispheric surface (Figure 1). Similar to the traditional SLIP model, the proposed model simulates bipedal walking through a simple mass, spring, and damper system, with legs equivalent to springs and dampers. The difference between the original and the proposed model lies in the fact that the proposed one can simulate the gait process and take the 3D dynamic parameters into account.

With the initial parameters optimized by the genetic algorithm (GA) method, the simulation results are consistent with the experimental results of the GRF well, which validates the 3D model. Due to the introduction of 3D kinematic parameters, the model proposed in this paper can generate GRF forces in three directions, rather than a calculation of the forces just in two directions conducted in previous models [810, 22, 24, 32, 4345]. The initial dynamics parameters of the 3D model have a significant effect on the shape and peak of the GRF, which is in line with those results presented in previous studies [22, 34, 46].

To thoroughly explore the relationship between several dynamic parameters, the correlation analysis is carried out with a linear correlation (Pearson’s r) test. The results show that there is a correlation between some of the initial parameters (Figure 12). A positive significant correlation is found between walking velocity and gait frequency and walking velocity and step length, a pattern similar to previous studies [33, 46, 47]. The results reveal a positive association between the frequency and stiffness and initial 2D angular velocity and a negative association with the impact angle. There also exists a negative and significant correlation between impact angle and frequency, velocity, stiffness, and damping ratio, in line with [22]. The initial 2D angle is significantly positively correlated with the impact angle and initial 3D angle; 2D angular velocity is negatively correlated with 3D angular velocity. A moderate positive correlation is found between damping and 3D angular velocity; a moderate negative correlation is found between impact angle and initial 3D angle. Regarding leg stiffness , our research suggests that there is a positive correlation between walking speed and stiffness, which is consistent with the results of [33, 48]. However, the finding that the damping ratio is not significantly associated with walking speed does not support the previous research [33]. The possible reason is that the spring constant taken by [33] is in the range of 15–30 kN/m, but we take it in the range of 17.5–22.5 kN/m. Contrary to the previous study [11], our study shows that there is a correlation between the initial angle, initial angular velocity, and damping with walking speed. The existence of 3D angles and 3D angular velocities may explain such an inconsistency.

A detailed parametric study is conducted to analyze the impact of initial dynamic parameter values on the GRF peaks. The results show that both 2D and 3D parameters have a nonnegligible effect on the peak of the GRF, but relative to the GRF in the vertical and longitudinal direction, the 3D parameters have the most significant effect on the GRF in the lateral direction. The possible explanation for the sensitivity of the lateral GRF to the 3D angular velocity may be as follows. The increase in absolute values of 3D angular velocity means that under the same impact angle and time for leg switching, the 3D angle will increase when the foot against the ground. Since the 3D angle is small and the lateral force is the product of the sine of the 3D angle and the spring force in the horizontal projection, the lateral force shows a positive correlation with the 3D angle and a negative correlation with 3D angular velocity (note: the 3D angular velocity is a negative value).

It must be acknowledged that our study has some limitations. First, the 3D model results in computational complexity for its adding of 3D variables entails the restriction of variables range. Second, the 3D model cannot perform a complete gait cycle continuously and spontaneously if compared to the model proposed by Lim and Park [48]. Third, the objective function of the genetic algorithm is set to be the root mean square error between the simulation results and the experimental data of the vertical GRF. Despite the advantages of simplified calculations, it must be conceded that the optimization approach may omit some correlation features between the parameters.

5. Conclusions

In this study, we proposed a three-dimensional bipedal walking model that could generate the repeated gait cycle with the simulated results of GRF in three directions fitting with the experimental data well. The correlation analysis between the initial dynamic parameters of the model was conducted. We examined the impact of 3D dynamics parameters on the peaks of GRF in three directions and found that the 3D parameters had a major effect on the lateral GRF. The results of our study implied that the dynamic characteristics presented by the bipedal walking pedestrian could be simulated by a spring-damping system.

Appendix

A. Lagrangian Formulation of Equations of Motion

The equations of motion of the bipedal 3D inverted pendulum model are derived from the Lagrangian formulation. Lagrange’s equation of motion is shown in the following:where represents a degree of the system and indicates the corresponding nonconservative force. The total potential energy of the system consists of the elastic potential energy of the spring and the gravitational potential energy of the center of mass:where denotes the spring stiffness of the leg, represents the slack length of leg, represents the length of the trailing leg, represents the length of the leading leg, indicates the angle between the limb and the plumb line passing through the center of mass (CoM), and indicates the corresponding velocity. represents the radius of the foot, represents the mass of the body, and represents the acceleration of gravity. For the convenience of equation derivation, the velocity of the center of mass is decomposed into three mutually perpendicular velocity components as , , and . The direction vector describes the direction of leg extension, describes the direction which is perpendicular to the plane defined by the leg and vertical plumb line, and is perpendicular to the direction of the leg and in the plane defined by the leg and vertical plumb line (Figure 15).

The total kinetic energy of the system can be given by

Lagrangian, , of the system can be obtained from equations (A.2) and (A.4) as follows:

The system has three mutually independent degrees of freedom, , , and , respectively, and their relationship can be obtained from the spatial velocity relationship, as shown in (B.4) in Appendix B. The incremental changes in as a function of , , and can be described from the Jacobian terms and can be given by

The Lagrangian equations of motion are given by

For a small unit of time, the work done by the damping force can be written as

According to the equations (B.5)–(B.7) in Appendix B, the damping forces , , and are given by

Substituting equations (A.6)–(A.8) and equations (A.13)–(A.15) into equations (A.9)–(A.12) results in the equations of motion of the system. For the convenience of calculation, this dynamic equation can be written in the matrix form as follows:

The bipedal phase of the gait can be expressed by equations (A.16)–(A.19). If not considering , , and in equation (A.16), the kinetic equations for the single-step phase can be obtained as

It is worth noting that the inertia matrix in equation (A.16) is not a symmetric positive definite matrix; this is because the walking model has kinematic constraints in the bipedal phase; the state of motion of the system requires only three variables as to be expressed. And the other three variables can be expressed by through equations (B.1)–(B.13).

B. Calculation of the Kinematic Constraints of the Motion

During the bipedal phase of a gait cycle, the velocity of the center of mass can be calculated in two ways. The velocity of the mass can not only be described by the variables , , and , but also by the variables , , and . Decomposing the velocity of the center of mass along the three directions X, Y, and Z, the equations can be obtained as follows:

The equation for velocity along the X-direction:

The equation for velocity along the Y-direction:

The equation for velocity along the Z-direction:

From equations (B.1)–(B.3), the relationship between , , and and , , and can be derived as follows:

The elements of the matrix are expressed separately for the readability of the formulas as follows:

Data Availability

The authors confirm that the data supporting the findings of this study are available within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors thank HY Liang (Soochow University) for revising the manuscript and polishing the language. This study was supported by the National Natural Science Foundation of China (Grant No. 51278106), Project of Science and Technology Research and Development Program of China Railway Corporation (Grant No. P2018G049), and Key R & D Project of Jiangsu Province (Grant No. BE2019107).