Abstract: The numerical simulation is critical for the seismic analysis of underground structures. However, its accuracy and efficiency are affected by the factors such as input of ground motion, boundary effect, model scale, which leads to the incompatibility among computational scale, time and accuracy. How to efficiently and accurately simulate the seismic response of the underground structure-ground interaction system is still an open question. A novel hybrid boundary element-finite element method in the framework of the domain reduction method (DRM) is proposed to efficiently simulate the seismic response characteristics of the subsurface structure-strata system. First, the overall subsurface structure-stratum model is divided into the inner domain sub-model of the near-field stratum structure and the outer domain sub-model of the far-field stratum, in which the displacement continuity at the inner-outer domain coupling boundary is ensured by overlapping nodes. Second, the non-free field vibration characteristics under the influences of free field and topography in the outer domain are solved by the boundary element method, and the dynamic response in the outer domain is converted into the equivalent seismic load by the DRM. The method greatly reduces the size of the outer domain, ensures the reasonable input of ground vibration, and realizes the rapid parametric analysis of the seismic response of subsurface structures in the inner domain. Further, two typical cases are designed to test the reliability and efficiency of the method. In the case of a two-line tunnel without the influences of topography, the accuracy of the method is validated by comparison with the reference solution. In another case of a two-line tunnel under the influences of terrain conditions, the numerical results show that compared with that of the remote boundary method and the traditional viscoelastic method, the computational cost is reduced by about 72% and 58%, respectively, and the computational scale is reduced by about 97% and 83%, respectively. In addition, the proposed method can be extended to the dynamic response analysis of subsurface structures under the action of oblique incident ground shaking.