The stiffness matrix given in [29] can be adopted for the dynamic stiffness matrix of the 3D layered foundation [ M P , S V , S H ] . If the uniform distributed vertical load moving along the y-axis direction has a distribution length of 2Δ along the x-axis direction, the density is q ₀, which can be expressed as q ( x , y , t ) = q 0 δ ( y − c t ) ϕ ( x ) , (1) where δ is the Dirac function; ϕ ( x ) is the distribution function, when | x | ≤ Δ , ϕ ( x ) = 1 , otherwise ϕ ( x ) = 0 ; c is the speed of the moving load. The Fourier transform of Eq. 1 can be obtained as q ≅ ( k x , k y , ω ) = 1 8 π 3 ∫ − ∞ ∞ ∫ − ∞ ∞ ∫ - Δ Δ q ⋅ e exp ( i ω t − i k x x − i k y y ) d x d y d t = q 0 ⁡ sin ( k x Δ ) 4 π 2 k x c δ ( k y − ω c ) , (2)

Here, k x and k y are the wavenumbers in the x and y directions, and ω is the circular frequency. In the layered elastic foundation, the external force and displacement on the interface of each soil layer satisfy the equilibrium equation [30] as { w ≅ x 0 , w ≅ y 0 , i w ≅ z 0 , ⋯ , w ≅ x N , w ≅ y N , i w ≅ z N } = [ M P , S V , S H ] − 1 { 0,0 , i q ≅ , ⋯ , 0,0,0 } . (3)

Among them, { w ≅ x 0 , w ≅ y 0 , i w ≅ z 0 , ⋯ , w ≅ x N , w ≅ y N , i w ≅ z N } is the displacement amplitude vector at the interface of the soil layer, and { 0,0 , i q ≅ , ⋯ , 0,0,0 } is the external force vector at the interface of the soil layer. Considering that there is only the vertical load on the surface, only the third element has a value, and the other elements are zero. [ M P , S V , S H ] is the 3D stiffness matrix of the layered half-space. Substituting Eq. 2 into Eq. 3, the displacement on the interface of each soil layer can be obtained. Subsequently, the amplitude of the upward and downward waves in each soil layer can be obtained from the relationship between the amplitude of the upward and downward waves in each soil layer and the displacement at the interface of each soil layer. Finally, the displacement { w ≅ x f , w ≅ y f , w ≅ z f } T and stress { T ≅ x f , T ≅ y f , T ≅ z f } T of any point in the frequency-wavenumber domain can be obtained from the relationship between the displacement and stress and the amplitude of the up and down waves.

For the scattered wavefield caused by the in-filled trench, let [ G u D ( x , y , z ) ] , [ G t D ( x , y , z ) ] be the displacement and stress Green’s function matrixes in the outer half-space of the in-filled trench and [ G u S ( x , y , z ) ] , [ G t S ( x , y , z ) ] be the displacement and stress Green’s function matrixes in the inner half-space of the in-filled trench, respectively [31]; then, the displacement and stress of any point outside and inside the in-filled trench can be expressed as { w ˜ x g D , w ˜ y g D , w ˜ z g D } T = [ G u D ( x , y , z ) ] { p ˜ 1 x , p ˜ 1 y , p ˜ 1 z } T , (4a) { T ˜ x g D , T ˜ y g D , T ˜ z g D } T = [ G t D ( x , y , z ) ] { p ˜ 1 x , p ˜ 1 y , p ˜ 1 z } T , (4b) { w ˜ x g S , w ˜ y g S , w ˜ z g S } T = [ G u S ( x , y , z ) ] { p ˜ 2 x , p ˜ 2 y , p ˜ 2 z } T , (4c) { T ˜ x g S , T ˜ y g S , T ˜ z g S } T = [ G t S ( x , y , z ) ] { p ˜ 2 x , p ˜ 2 y , p ˜ 2 z } T . (4d)

Here, the subscripts “u” and “t”, respectively, represent the displacement and stress Green’s function; the subscript “G” represents the dynamic response caused by the moving uniform distributed line load; the superscripts “D” and “S”, respectively, represent the layered elastic half-space outside the in-filled trench and the elastic soil layer inside the in-filled trench. { p ˜ 1 x , p ˜ 1 y , p ˜ 1 z } T and { p ˜ 2 x , p ˜ 2 y , p ˜ 2 z } T are virtual distributed line load.

The continuous conditions of stress and displacement on the interface s of the in-filled trench can be expressed as ∫ S [ W ( s ) ] T ( [ T ˜ x g D ( s ) T ˜ y g D ( s ) T ˜ z g D ( s ) ] + [ T ˜ x f ( s ) T ˜ y f ( s ) T ˜ z f ( s ) ] ) d s = ∫ s [ W ( s ) ] T [ T ˜ x g S ( s ) T ˜ y g S ( s ) T ˜ z g S ( s ) ] d s , (5a) ∫ S [ W ( s ) ] T ( [ w ˜ x g D ( s ) w ˜ y g D ( s ) w ˜ z g D ( s ) ] + [ w ˜ x f ( s ) w ˜ y f ( s ) w ˜ z f ( s ) ] ) d s = ∫ s [ W ( s ) ] T [ w ˜ x g S ( s ) w ˜ y g S ( s ) w ˜ z g S ( s ) ] d s . (5b)

Here, T ˜ x g D ( s ) , T ˜ y g D ( s ) , T ˜ z g D ( s ) , w ˜ x g D ( s ) , w ˜ y g D ( s ) , and w ˜ z g D ( s ) are the stress and displacement along the x-, y-, and z-coordinate directions generated by the layered elastic half-space out of the in-filled trench acting on each element of the in-filled trench boundary. T ˜ x g S ( s ) , T ˜ y g S ( s ) , T ˜ z g S ( s ) , w ˜ x g S ( s ) , w ˜ y g S ( s ) , and w ˜ z g S ( s ) are the stress and displacement along the x-, y-, and z-coordinate directions generated by the layered elastic half-space inner of the in-filled trench acting on each element of the in-filled trench boundary. T ˜ x f ( s ) , T ˜ y f ( s ) , T ˜ z f ( s ) , w ˜ x f ( s ) , w ˜ y f ( s ) , and w ˜ z f ( s ) are the stress and displacement along the x-, y-, and z-coordinate directions generated by the free field on each element at the boundary of the in-filled trench. [ W ( s ) ] is the weight function matrix, which can be taken as the unit matrix, so that the integral can be performed independently on each element. Substituting Eq. 4 into Eq. 5 can be obtained [ T p L ] { p ˜ 1 x , p ˜ 1 y , p ˜ 1 z } T + [ T f ] = [ T p T ] { p ˜ 2 x , p ˜ 2 y , p ˜ 2 z } T , (6a) [ V p L ] { p ˜ 1 x , p ˜ 1 y , p ˜ 1 z } T + [ V f ] = [ V p T ] { p ˜ 2 x , p ˜ 2 y , p ˜ 2 z } T , (6b) where [ T p L ] = ∫ S [ W ( s ) ] T [ G t D ( s ) ] d s , (7) [ T p L ] = ∫ S [ W ( s ) ] T [ G t D ( s ) ] d s , (8) [ T f ] = ∫ S [ W ( s ) ] T [ T ˜ x f ( s ) , T ˜ y f ( s ) , T ˜ z f ( s ) ] T d s , (9) [ V p L ] = ∫ S [ W ( s ) ] T [ G u D ( s ) ] d s , (10) [ V p T ] = ∫ S [ W ( s ) ] T [ G u S ( s ) ] d s , (11) [ V f ] = ∫ S [ W ( s ) ] T [ w ˜ x f ( s ) , w ˜ y f ( s ) , w ˜ z f ( s ) ] T d s . (12)

From Eq. 6, the virtual moving uniform distributed line loads { p ˜ 1 x , p ˜ 1 y , p ˜ 1 z } T and { p ˜ 2 x , p ˜ 2 y , p ˜ 2 z } T applied on the boundary of the filling trench can be obtained, and substituting them into Eq. 4 can obtain the scattered wave field. Finally, by combining the free wave field and the scattered wavefield, the total displacement of any point in the layered elastic half-space can be obtained { u ˜ x , u ˜ y , u ˜ z } T = { G u x ( x , y , z , ω ) , G u y ( x , y , z , ω ) , G u z ( x , y , z , ω ) } T q ˜ ( y , ω ) . (13)

Among them, G u x ( x , y , z , ω ) , G u y ( x , y , z , ω ) , and G u z ( x , y , z , ω ) , respectively, represent the dynamic flexibility coefficient of the layered elastic foundation with in-filled trench under the unit moving uniform distributed line load with a width in the x, y, and z directions; q ˜ ( y , ω ) is the vertical uniformly distributed line load; and { w ˜ x , w ˜ y , w ˜ z } T is the displacement vector generated by the elastic foundation under moving load. The displacement, velocity, and acceleration of any point in the time domain can be obtained by the inverse Fourier transform as shown in the following: { w x ( x , y , z , t ) w y ( x , y , z , t ) w z ( x , y , z , t ) } = ∫ − ∞ ∞ { w ˜ x ( x , y , z , ω ) w ˜ y ( x , y , z , ω ) w ˜ z ( x , y , z , ω ) } e i ω t d ω , (14a) { w ˙ x ( x , y , z , t ) w ˙ y ( x , y , z , t ) w ˙ z ( x , y , z , t ) } = ∫ − ∞ ∞ i ω { w ˜ x ( x , y , z , ω ) w ˜ y ( x , y , z , ω ) w ˜ z ( x , y , z , ω ) } e i ω t d ω , (14b) { w ¨ x ( x , y , z , t ) w ¨ y ( x , y , z , t ) w ¨ z ( x , y , z , t ) } = ∫ − ∞ ∞ ( − ω 2 ) { w ˜ x ( x , y , z , ω ) w ˜ y ( x , y , z , ω ) w ˜ z ( x , y , z , ω ) } e i ω t d ω . (14c)

Assuming that the interaction pressure density between the track and the foundation is q ( y , t ) , the vibration equation of the track’s vertical displacement Ω can be described as [32]: EI ∂ 4 Ω ∂ y 4 + M ∂ 2 Ω ∂ t 2 = 2 Δ q ( y , t ) + ∑ n = 1 M p n ( y − c t + L n ) . (15)

Here, ∑ n = 1 M p n (y-ct + L _n ) is the wheel load of the train, and its specific parameters can be found in [32]. Equation 15 is converted into the frequency-wavenumber domain to obtain ( EI k y 4 − M ω 2 ) Ω ¯ ˜ ( k y , ω ) = 2 Δ q ¯ ˜ ( k y , ω ) + ∑ n = 1 M p ¯ ˜ n ( k y , ω ) . (16)

Considering that the displacement of each point of the track in the transverse direction is Ω , the continuous condition of the track center displacement and the corresponding foundation point is q ¯ ˜ ( k y , ω ) G w ( 0,0,0 , ω ) = Ω ¯ ˜ . (17)

Here, G w ( 0,0,0 , ω ) is the dynamic flexibility of the layered foundation with an in-filled trench under the action of the moving load in the z-direction. It expresses the vertical displacement of the midpoint of the element due to the uniform distributed vertical load acting on the element, which can be obtained by Eq. 13. Combining Eqs 16 and 17 can obtain the interaction force between the track and the layered foundation and the track displacement in the frequency-wavenumber domain. Then, the inverse Fourier transform is used to transform the result into the frequency-space domain, and then the displacement response of any position of the elastic foundation in the frequency-space domain can be obtained by Eq. 13. Finally, the displacement, velocity, and acceleration responses of any position of the elastic foundation in the time-space domain are obtained from Eq. 14.