Study of Rock Burst Risk Evolution in Front of Deep Longwall Panel Based on Passive Seismic Velocity Tomography

Key Laboratory of Deep Coal Resource Mining, Ministry of Education, China University of Mining and Technology, Xuzhou 221116, China School of Mines, China University of Mining and Technology, Xuzhou 221116, China Laboratory of Mine Earthquake Monitoring and Prevention, China University of Mining and Technology, Xuzhou 221116, China School of Energy and Safety, Anhui University of Science and Technology, Huainan 232001, China Shaanxi Zhengtong Coal Industry Co., Ltd., Xianyang 713600, China


Introduction
Due to long-term exploitation, coal resources at shallow depths are gradually exhausted, and mining depth of coal is continually increasing. Mining depth in overseas including Poland, Germany, Britain, Japan, and France had already exceeded 1000 m as early as in the 1980s, and now, it has reached 1500 m [1,2]. Currently, 47 coal mines in China have exceeded 1000 m [3]. With mining depth increasing, the stress state of surrounding rock is continuously deteriorating and rock burst occurs more frequently [4,5].
Various rock burst mechanisms have been put forward by different methods [6][7][8][9]. A universally acknowledged viewpoint is that rock burst is caused by the superposition of dynamic and static stress, as shown in Figure 1. The "contribu-tion rate" of static load and dynamic load for rock burst risk varied with mining depth. Due to high static stress in deep mining, a slight dynamic load increment caused by mininginduced tremor can make superimposed stress exceed the critical value and even lead to a rock burst. Prediction methods for stress field and rock burst risk in deep buried uncertain environments should be further investigated in the future [10].
Seismic velocity tomography, a new geophysical exploration method inferring the wave propagation velocity through structures, has been a novel measurement method for stress redistribution in underground coal mining. Superior to previous methods such as drilling bits [11], pressure sensor [12,13], electromagnetic emission [14], and acoustic emission [15], seismic velocity tomography can provide comprehensive and continuous stress redistribution by imaging the seismic wave velocity in coal and rock body. Regarding seismic sources, this method can be classified into two types, that is, active velocity tomography and passive velocity tomography [16,17]. Active velocity tomography has been conducted in geological structure and stress field detection and achieved good results [18][19][20][21][22]. However, due to extra labor and economic cost, the application of active velocity tomography is limited to some extent. Alternatively, passive seismic velocity tomography can rapidly and continuously present the stress redistribution during coal mining by using mining-induced seismicity as the sources seismic wave [23][24][25][26][27]. Passive seismic velocity tomography in existing literatures was mainly used to represent stress redistribution. The relationship between wave velocity change and rock burst risk was rarely inferred.
In this study, a quantitative assessment criterion for rock burst risk was established based on the relation between stress and wave velocity [28][29][30][31]. Rock burst risk and range of a deep longwall panel with folds and adjoining goaf were determined continuously and visually based on passive seismic velocity tomography. Moreover, the influence of pressure-relief measures on rock burst risk was analyzed based on tomography results.

Theory of Passive Seismic Velocity Tomography for Rock Burst Risk
During the transmission of seismic wave, P-wave is first monitored as it travels faster [23]. P-wave velocity variation with stress under different loading schemes has been conducted in laboratory [28], and it indicates the positive correlation between stress and P-wave velocity. P-wave velocity can reflect stress state and burst risk of rock and coal mass. A 3D mesh network of mining area is necessary for passive seismic velocity tomography, and mining area is divided into voxels in x, y, and z directions [25,31]. During the propagation, seismic rays will pass through voxels along the ray path from seismic source to sensors.
Suppose the ray path of the ith seismic wave is L i and the travel time is T i , then, the travel time of ith seismic wave from the source to the sensor is the integral of the slowness S (or the inverse of velocity), which can be expressed by Eq.
where VðX, Y, ZÞ is the velocity (m/s), L i is the ray path of the ith seismic wave (m), T i is the travel time (s), SðX, Y, ZÞ is the slowness (s/m), d ij is the distance of the ith ray in the jth voxel, N is the total number of rays, and M is the number of voxels. Generally, seismic event location and ray path are calculated using an initial velocity model [31]. Due to the unknown velocity, distance, and time in an individual voxel, thus, matrices expressed by Eq. (4) can be built with the voxel slowness, distance, and time. Then, the velocity can be determined by the following matrix [33].
where T is the travel time per ray matrix (1 × N), D is the distance per ray per voxel matrix (N × M), and P is the slowness per grid cell matrix (1 × M). Eventually, the key problem of seismic velocity tomography is to solve the slowness vector S. The most effective way to solve this problem is through an iterative process, and the simultaneous iterative reconstructive technique (SIRT) [16,29,34] is the famous and effective one, which is adopted in the paper.
Velocity of each voxel can be determined by seismic velocity tomography introduced above. The relationship between seismic velocity anomaly and stress coefficient was given in [30,31]. In deep coal mining, burst risk is primarily dominated by static stress. Accordingly, an assessment  2 Geofluids criterion for burst risk in Table 1 is built. Velocity anomaly A n is determined by Where V p is P-wave velocity in certain voxel, and V a p is the average velocity of the model. It should be noted that the zones with positive anomaly and negative anomaly are overstressed and pressure-relieved, respectively [31]. In this paper, the mining area was classified as None, Weak, Middle, and Strong burst risk zones according to rock burst risk assessment criterion in Table 1.

Site Characteristics of Selected
Longwall Panel There are two inferior key strata and one main key stratum above coal seam, as shown in Figure 2 No.4 coal seam is of strong burst propensity, and the rock strata in the floor and roof are both of weak burst propensity. The maximum principal stress reaches 38.2-44.8 MPa, and the average uniaxial compressive strength of coal is 19.3 MPa; hence, the stress concentration coefficient is 1.98-2.32, which indicates its high burst risk [35].

Microseismic Data Acquisition and Processing.
A microseismic monitoring system called "SOS," manufactured by Central Mining Institute of Poland, was installed in the coal mine, and the maximum locating errors are 20 m in horizontal direction and 30 m in vertical direction, respectively. Seis-mic monitoring network in July 2019 is shown in Figure 3. The system has been optimized three times during the panel 204 retreating to ensure the monitoring accuracy. During April 2019 to March 2020, over 3800 mining-induced tremors in No.2 mining district were recorded by SOS. The first arrival time of P-wave in each sensor was calibrated manually till the error between the calculated and theoretical values is less than 20 ms. Some typical tremors in panel 204 are illustrated in Figure 4(b). Seismic wave attenuates by power function as Eq. (6). Seismic waves with higher energy attenuate slower and travel farther, which is consistent with previous studies [36,37]. Consequently, the quantity and accuracy of microseismic events can satisfy the tomography.
Where L is the distance between tremor and sensor, V 0,max and V 0 ðLÞ are the particle velocity at tremor and sensor, respectively, and λ is the attenuation coefficient.

Inversion Parameters.
Tomography area is 2200 m long, 2560 m wide, and 360 m high, which is divided into voxels by 20 m in x and y directions and 30 m in z direction, respectively, as shown in Figure 3. More ray paths enable a higher accuracy in the computation [38]. Thus, only tremors recorded by more than six sensors are considered as seismic sources, and the voxels with more than 10 rays are considered reliable. To improve the efficiency and accuracy of inversion and source locating, SIRT was adopted to recalculate the seismic location, and the slowness in each voxel along the seismic rays was modified by iteration till the threshold value was reached. To start the first iteration and reduce indeterminacy, P-wave initial velocity its range was assumed 4.48 km/s and 3.5 km/s-6.5 km/s, respectively, which was obtained from the P-velocity data in Figure 5.

Rock Burst Risk Evolution with Longwall Panel
Retreating. Passive seismic velocity tomography in the mining area was carried out regularly with 2 times per month since May 2019. Thereinto, 8 tomography results of panel 204 are illustrated in Figure 6. Mining-induced tremors in the following one month or half one month are plotted simultaneously to verify the burst risk zones highlighted by tomography. In Figure 6, the velocity anomaly (A n ) in green, yellow, and red area is 0.05-0.15, 0.15-0.25, and >0.25, and the corresponding burst risk is weak, middle, and strong, respectively, based on the rock burst risk criterion given in Table 1.
Notably, panel 205 was in production since November 2019, and tremors were induced in the goaf of panel 204 and abutment area ahead of panel 205, which would bring about seismic rays and high burst risk zones across the goaf. The distance between panels 204 and 205 always exceeds 700 m. Hence, the extraction of panel 205 will not interfere with this study. To illustrate burst risk in front of panel 204 clearly, burst risk zones in the goaf were eliminated artificially.
As clearly shown in Figure 6, the general evolution of rock burst risk and range in front of the face line is closely 3 Geofluids related to the retreating of panel 204. High burst risk zone scale ahead of the face line in the strike, including middle and strong burst risk zones, is illustrated in Figure 7. Due to increasingly broadened goaf behind panel 204, goaf of panel 203, and the folds, high burst risk zone scale in front of face line first increases from 60 m, reaches the peak of 410 m, and then decreases to 80 m gradually. The peak is reached when panel located in 1# syncline shaft area in October 2019 (as shown in Figure 6(e)). Burst risk is the highest at that time as well. When the face line is close to the crossheading, high burst risk zones distribute along the crossheading and further intersect with the burst risk zones in 1# syncline area, which indicates the high stress level of the zones adjacent to the crossheading.
In terms of high burst risk along the inclination of panel 204, the range of burst risk zones on the ventilation roadway side is always larger than that on the haulage roadway side. However, the evolution of high burst risk shows distinct zoning features. Hence, the mining process of panel 204 can be divided into 3 phases along the strike, as shown in Figure 8.   Moreover, the distribution of mining-induced tremors over 1E4 J in the following mining period (one month or half one month after tomography) are well corresponding to high burst risk zones determined by the previous tomography above.      9 Geofluids microseism activity in Figure 11. In addition, from May to June 2019, high working resistance zones are mainly located on the haulage roadway side, and it is identical to the high burst risk zones in Figures 6(a) and 6(b). From September to November 2019, high working resistance zones are mainly distributed on the ventilation roadway side, and it is identical to the high burst risk zones in Figures 6(c)-6(f). From January to March 2020, high working resistance zones are mainly on the haulage roadway side, which is identical to the high burst risk zones in Figures 6(g) and 6(h).
In summary, the field drilling bit results, rock burst occurrence, microseism activity, and working resistance of hydraulic supports are well identical to the tomography results, which indicate the feasibility and accuracy of passive seismic velocity tomography for rock burst risk.   10 Geofluids in the intensified pressure-relief zones, and the microseismic activity decreased obviously (as shown in Figure 11), which indicates that pressure-relief measures have a significant influence on rock burst risk.
To investigate the influence of pressure-relief measures on rock burst risk visually and quantitatively, pressurerelief measures and zones in panel 204 from September 2018 to May 2019 were analyzed, and three passive tomographies were conducted in early May 2019.
Big-diameter borehole is the primary pressure prerelief method. Secondary big-diameter boreholes or coal blasting boreholes will be implemented again at interval of primary boreholes if the risk is not effectively eliminated. Bigdiameter borehole is 153 mm in diameter, 25 m in depth, and 1 m in interval. Coal blasting borehole is 42 mm in diameter, 10 m in depth, 5 m in interval, 1.2 m from the floor, and charged with 3 kg explosive, as shown in Figure 13.
As illustrated in Figure 14 respectively. Then, the corresponding high burst risk zones are dispelled obviously, as shown in Figure 14(b). From 7 May to 12 May, 47 and 30 boreholes were drilled at haulage and ventilation roadway, respectively. Similarly, the corresponding burst risk zones are eliminated as well, especially for the zones on the ventilation roadway side, as shown in Figure 14(c).
Generally, with production, stress field is adjusted; high stress and burst risk zones are shifted closer to the face line in the strike direction and to haulage or ventilation roadway in the inclination direction. Stress level and burst risk increase simultaneously. However, burst risk zones on the ventilation roadway side, contrarily, are being away from the face line. Burst risk eliminated zones match well with the pressure-relieved area, which indicates that pressurerelief measures and intensity in ventilation roadway can decrease burst risk effectively. However, high burst risk in zones 200 m in front of the face line on the haulage roadway side is not eliminated substantially, which is further proved by the future tremors in May 2019 as shown in Figure 6(a). Hence, more pressure-relief measures, such as coal blasting or deep presplitting blasting in roof, should be adopted till the risk is eliminated thoroughly. Therefore, pressure-relief measures have a distinct influence on stress field and burst risk of panel, which can account for the none or weak burst risk zones closely ahead of the face line in Figure 6. And moreover, passive seismic velocity tomography before and after the pressure-relief measures implemented can be a novel and effective method to assess pressure-relief effect.

Conclusion
Based on passive seismic velocity tomography results and rock burst risk assessment criterion, the burst risk and range evolution of a deep longwall panel were determined, and the influence of pressure-relief measures on burst risk was analyzed. The main conclusions are as follows: (1) Seismic wave propagation velocity in rock and coal mass is in positive correlation with the stress level. Seismic wave velocity distribution of rock and coal mass can be used to assess the rock burst risk in deep coal mining. And accordingly, a rock burst risk assessment criterion with velocity analogy was built (2) Passive tomography results indicate that due to the tectonic stress, abutment pressure, and mining layout, rock burst risk and range of the deep panel firstly increase, then decrease, and reach the peak at the 1# syncline shaft area during panel retreating. High burst risk zones in the inclination of panel show distinct zoning features. When panel approaching the crossheading, high burst risk zones distribute along the crossheading and further intersect with the burst risk zones in 1# syncline shaft area (3) High burst risk zones identified by passive seismic velocity tomography are well correlated with drilling bits results, rock burst records, microseism activity, and the working resistance of hydraulic supports, which indicates the practicability and accuracy of passive tomography for rock burst risk in deep coal mining (4) Pressure-relief measures and mining layout of panel have a distinct influence on rock burst risk, which can be assessed by passive seismic velocity tomography. Pressure-relief measures and intensity in different mining phases should be optimized timely based on the tomography results

Data Availability
The figures and tables used to support the findings of this study are included in the article.

Conflicts of Interest
The authors declare no conflict of interest.