## 1 Introduction  
Underground rock masses are inevitably in a certain geological and tectonic environment, and are subject to initial stresses such as gravitational stress, tectonic stress, temperature stress, etc. The initial stresses have a strong effect on the mechanical behavior of the underground rock masses and the stability of the underground engineering [1, 2, 3]. Meanwhile, the drill-and-blast method is the most widely used technique for tunnel excavation and underground mining. In this process, the initial stressed underground rock masses are under strong dynamic disturbances, and the underground structure can be damaged by the stress waves generated during the excavation process. Hence, it is of great practical significance to study the propagation of stress waves in underground rock masses under initial stress for the optimal design of underground rock mass blasting excavation parameters and the dynamic stability analysis of underground engineering. Rock masses contain various types of discontinuous interfaces such as joints and fractures and so on, which have a noticeable effect on the mechanical response of the rock mass [4, 5]. Discontinuous interfaces in natural rock mass are usually distributed in groups, such as a stratified rock mass, and it is particularly important to study the propagation of blasting stress waves in the layered rock mass and to monitor the vibrations generated during the blasting of the stratified rock masses [6]. Intensive studies have been conducted to investigate the propagation of stress waves across jointed rock masses via various theoretical and experimental methods. In terms of theoretical research, the displacement discontinuity model (DDM) proposed by [PERSON] [7] has been widely applied to study stress wave propagation through a jointed rock mass [8; 9]. The DDM was also combined with other analysis methods, e.g., the method of characteristics (MC) [10], the scattering matrix method [11], and the time-domain recursive method (TDRM) [12] to study the stress waves passing through linear and nonlinear joints [13; 14], one single joint and a set of parallel joints [15; 16], and even intersecting rock joints [17]. By experimental means to date, the split Hopkinson pressure bar (SHPB) apparatus has been mainly used to study stress wave propagation across rock masses [18; 19]. In contrast, the numerical simulation method is an economical and feasible alternative to survey the stress wave propagation across a jointed rock mass. Based on the discrete element method (DEM) proposed by [PERSON] [20], the universal distinct element code (UDEC) has been widely used to calculate the propagation problems of stress waves in a jointed rock mass [21; 22; 23]. Furthermore, other numerical methods and software have been adopted to solve the problems involving the stress wave propagation in a rock mass, e.g., the particle manifold method (PMM) [24; 25], the numerical manifold method (NMM) [26; 27], the particle flow code (PFC) [28; 29], and the three-dimensional element code (3 DEC) [30]. However, the above theoretical, experimental, and numerical methods have mainly focused on the effect of the parameters of the joints, e.g., joint stiffness, joint spacing, joint number, and the parameters of the stress wave, e.g., waveform, amplitude, frequency, and the incident angle of stress wave on the stress wave propagation pattern, and have proposed that the attenuation of the stress wave only occurs at the joints, while it has been assumed that the intact rock is elastic. Few works have been conducted that investigate the effect of initial stresses on the stress wave propagation in the jointed rock mass, and studies considering the initial damage of intact rocks in the rock mass under initial stresses are much rarer. On the other hand, besides discontinuous interfaces, intact rocks are the other part of the rock mass.