## 2 Semi-Analytical Solution for the Coupled Track-Ground-Tree Grid System
### Fundamental Solutions to the Dynamic Response of a Stratified Half-Space  
A stratified ground system composed of \(Q+1\) layers, where the lowermost layer is idealized as a homogeneous elastic half-space, is assumed to be isotropic and linearly elastic, as shown in Figure 1. The dynamic response of this layered medium under harmonic point loading is evaluated using the thin-layer method (TLM) [52]. The wave propagation behavior in the homogeneous elastic domain is governed by the following frequency-domain equations [53]:  
\[(\lambda+2\mu)\
abla\
abla\cdot\mathbf{u}-\mu\
abla\times(\
abla\times \mathbf{u})+\rho\mathbf{b}=\rho\mathbf{\ddot{u}} \tag{1}\]  
where \(\mathbf{u}=[u_{x},u_{y},u_{z}]^{\mathrm{T}}\) denotes the responses in three directions, \(\mathbf{b}=[b_{x},b_{y},b_{z}]^{\mathrm{T}}\) denotes the force vector applied, \(\lambda\) and \(\mu\) denote the Lame constants, and \(\rho\) is the density of the elastic mass. The derivation of the solution in the frequency-wavenumber domain necessitates the use of two types of Fourier transforms [54]:  
\[ \begin{array}{l}\hat{f}(\omega)=\int\limits_{-\infty}^{\infty}f(t)e^{-i \omega t}\mathrm{d}t,f(t)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\hat{f}( \omega)e^{i\omega t}\mathrm{d}\omega\\ \overline{\hat{f}}(k_{x},k_{y})=\int\limits_{-\infty}^{\infty}\int\limits_{- \infty}^{\infty}f(x,y)e^{i(k_{x}x+k_{y}y)}\mathrm{d}x\mathrm{d}y,f(x,y)=\frac {1}{4\pi^{2}}\int\limits_{-\infty}^{\infty}\overline{\hat{f}}(k_{x},k_{y})e^{- i(k_{x}x+k_{y}y)}\mathrm{d}k_{x}\mathrm{d}k_{y}\end{array}  \tag{2}\]  
where the superscript symbols \"\(\cdots\)\", \"\(\cdots\), and \"\(\cdots\)\" represent the Fourier transformations performed over time \(t\) and spatial coordinates \(x\) and \(y\), respectively. Following the procedure outlined by [PERSON] [55] and based on the integration of Equations (1) and (2) together with the discretization of the ground system into multiple thin layers, as shown in Figure 2, the governing equation in the frequency-double wavenumber domain is derived as follows:  
\[\mathbf{\widehat{\mathbf{P}}}=\left[k_{x}^{2}\mathbf{A}_{xx}+k_{x}k_{y} \mathbf{A}_{xy}+k_{y}^{2}\mathbf{A}_{yy}+i(k_{x}\mathbf{B}_{x}+k_{y}\mathbf{B }_{y})+(\mathbf{G}-\omega^{2}\mathbf{M})\right]\mathbf{\widehat{\mathbf{U}}} \tag{3}\]  
A notable strength of the TLM is its ability to analytically perform at least a single inverse Fourier transform during the derivation process of the fundamental solution in the spatial domain. Such a process is carried out using the modal superposition method [53], in which the response is expressed as a combination of basic functions. Accordingly, the coefficient matrices in Equation (3) are initially assembled following the conventional finite element (FE) procedure to establish the system of equations. Figure 2: Schematic diagram of discretization of the half-space into thin layers. (DOFs) are then reorganized by grouping the horizontal (\(x\)), horizontal (\(y\)), and vertical (\(z\)) components together. Finally, the equation for eigenmodes is expressed as follows:  
\[[k^{2}\mathbf{A}+k\mathbf{B}+(\mathbf{G}-\omega^{2}\mathbf{M})]\widehat{ \overline{\Psi}}=0 \tag{4}\]  
where \(k\) represents the characteristic value and \(\widehat{\overline{\Psi}}\) represents an undetermined eigenmode. The specific structures of the matrices in Equation (4) are detailed as follows:  
\[\mathbf{A} =\left[ \begin{array}{ccc}\mathbf{A_{1}}&\mathbf{O}&\mathbf{O}\\ \mathbf{O}&\mathbf{A_{2}}&\mathbf{O}\\ \mathbf{O}&\mathbf{O}&\mathbf{A_{3}}\end{array} \right]^{\mathbf{T}},\mathbf{B } =\left[ \begin{array}{ccc}\mathbf{B_{1}}&\mathbf{O}&\mathbf{O}\\ \mathbf{O}&\mathbf{B_{2}}&\mathbf{O}\\ \mathbf{O}&\mathbf{O}&\mathbf{B_{3}}\end{array} \right]^{\mathbf{T}}, \tag{5}\] \[\mathbf{G} =\left[ \begin{array}{ccc}\mathbf{G_{1}}&\mathbf{O}&\mathbf{O} \\ \mathbf{O}&\mathbf{G_{2}}&\mathbf{O}\\ \mathbf{O}&\mathbf{O}&\mathbf{G_{3}}\end{array} \right]^{\mathbf{T}},\mathbf{M } =\left[