Attenuation characteristics analysis of seismic energy and its application to risk assessment in underground coal mines

Abstract Seismic energy analysis is the basis of coal burst risk assessment. The prerequisite for utilizing seismic energy to quantify coal burst hazards is the accurate understanding of seismic energy attenuation. Thereby the characteristics of seismic energy attenuation were first systematically investigated in this paper. Based on the seismic data and 24 coal burst records during the period of sixteen months, a risk assessment method including three indexes (Static Intensity Index (SII), Dynamic Intensity Index (DII), and Risk Assessment Index (RAI)) was derived from seismic energy attenuation. In consideration of the static and dynamic response, the superposition effects of seismic energy were proposed to improve the performance of risk assessment. Results show that residual seismic energy (RSE) is relatively correlated with seismic energy, source radius, and energy attenuation coefficient; the higher seismic energy, source radius, and lower energy attenuation coefficient would result in higher coal burst risk. The case study demonstrated that SII could effectively predict the damaged area, but its prediction efficiency of high-magnitude events (HMEs) is relatively low; DII could predict damaged areas and HMEs, but may lead to risk overestimation; and RAI could efficiently predict the damaged area and HMEs. The sensitivity of the monitoring system and seismicity activeness could be the two main restrictions of the proposed method, leading to low efficiency or even false risk assessment. This paper sheds light on seismic energy utilization to risk prediction in underground coal mines. Highlights The attenuation response characteristics of seismic energy to different factors including seismic energy, source radius, and attenuation coefficient were investigated. A new assessment method based on seismic energy attenuation calculation was proposed to assess coal burst risks. Three seismic-based indexes considering static and dynamic response and their superposition effects were proposed to assess coal burst risk. Sensitivity of monitoring system and seismicity activeness could be the two main restrictions leading to low accuracy or even false in risk assessing, which implies that the proposed method could also have false judgements in the place where rare seismicity is monitored.


Introduction
Coal bursts have been recognized as some of the most hazardous issues that threaten coal mines and often cause fatalities and enormous property loss due to the release of strain energy (Cook 1965;Ortlepp and Stacey 1994;Dou et al. 2012;Cao et al. 2015;He et al. 2017). Recently, with the increase of mining depth and increasingly complex geological conditions, multiple casualties caused by coal bursts were recorded in many countries, including China, Poland, Australia, South Africa, USA, Russia, and Canada (Ortlepp 2005;Todd and Newsome 2014;Zhang et al. 2017;Keneti and Sainsbury 2018;Li et al. 2018;Chen et al. 2021). Therefore, more investigation should be conducted to obtain a comprehensive understanding of the inducing mechanism, quantitative risk assessment, and prevention of coal bursts.
Over past years, several hypotheses were put forward for the mechanism of coal burst, including strength (Romashov and Tsygankov 1993), rigidity (Chen et al. 1997;Ma et al. 2018), burst liability (Kidybi nski 1981;Cai et al. 2016), stability (Cook 1965;Linkov 1996) and energy theory (Hudson et al. 1972;Xu et al. 2017). Among these hypotheses, the inducing mechanism of coal bursts was attributed to static conditions. However, this hypothesis cannot explain the occurrence of serious coal bursts under the condition of shallow burial and low in-situ stress (e.g. a coal burst at the Kuangou coal mine, China, with the mining depth of 317 m (Lan et al. 2012)) on one hand, and fails to take into account dynamic loads, i.e. loading rate effects and mining seismicity on the other. To solve these problems, the theory of coupled static and dynamic stress was proposed by Dou et al. (2014) to interpret the coal burst mechanism. As is proposed, coal bursts would be initiated provided the total stress, that is, the superposition of the static and dynamic stress exceeds the critical stress. Moreover, Dou et al. (2014) first pointed out the mining-induced abutment stress should be regarded as static stress, and mining seismicity could cause dynamic stress. Consequently, assessment and prevention of coal bursts should take static stress, dynamic stress, and their superposition effects into consideration.
Several methods for dynamic disaster monitoring and prediction have been developed, such as drilling bits (Braeuner 1994;Gu et al. 2012), acoustic emission (Qi and Liu 1994;Dou and He 2001), borehole detection (Qu et al. 2011), borehole pressure monitoring (Zhang et al. 2014), electromagnetic emission Wang et al. 2011) and hydraulic support pressure , etc. These methods mentioned above meet the needs of mining operations to some extent, whereas they are not suitable for large-scale monitoring due to the shortcomings of localized monitoring range, interference with production. In practice, seismic monitoring has been proven to be the most convenient and effective tool to quantify seismicity and monitor coal burst risks for the following merits:continuous remote transmission, wide monitoring range, strong anti-interference ability, and no interference with production.
Based on seismic monitoring, considerable typical dynamic disaster prediction methods have been expended via different seismic processing techniques, including source parameters (Glazer 2018), source mechanism (Ortlepp 1997;Gibowicz and Lasocki 2001;Gibowicz 2009), velocity tomography , seismic clustering (Mason 1981;Si et al. 2015;Cao et al. 2016) and seismic energy and frequency (Urbancic et al. 1992;Cao et al. 2018;Si et al. 2020), etc. Meanwhile, several seismic indices have been proposed: b value, moment tensor, apparent stress/volume, fractal dimension, seismic energy, and frequency, etc. It is observed that seismic energy is one of the most fundamental and widely used indicators. Wang et al. (2016) investigated the evolution characteristics of seismic energy prior to the coal burst by means of energy density index. He et al. (2017) constructed the dynamic stress disturbance coefficient based on seismic energy to evaluate dynamic disturbance. Nevertheless, attenuation, an essential feature of seismic energy, was not considered in their research. Cai et al. (2018) developed a new coal burst forecasting method based on seismic strain energy. Recently, Cai et al. (2020) tried to estimate dynamic stress through seismic energy characteristics. The fitted attenuation coefficient was adopted to measure the attenuation of seismic energy, whereas the source size, an indispensable parameter in seismicity analysis, was not taken into account. Wang et al. (2020) developed the dynamic Load Index to assess the risk of dynamic loads caused by seismic events. The attenuation characteristics of seismic energy (including source size, fracture mode, and quality factor of rock mass) were investigated via theoretical analysis. However, in the process of assessing hazards, too many parameters in this theoretical model need to be field measured. Due to the complexity and variety of the mining environment, it is difficult to obtain the accurate measurement of these parameters. Therefore, it is of great necessity to develop a simple and reasonable method to describe attenuation characteristics of seismic energy and evaluate coal burst hazards. This paper proposes a new seismic-based method to assess coal burst hazards in underground coal mines. To refine the method, seismic data and damage records from a burst-prone longwall panel in a Chinese coal mine for sixteen months were used. The seismic energy attenuation characteristics were systematically investigated. Based on seismic energy attenuation calculation, three seismic-based indexes which considered the static response, dynamic response, and their superposition effects respectively were proposed to assess coal burst risk. The new method was validated through high-magnitude events and damaged records of a coal burst case. Combined with the monitoring record of coal burst in the coal mine, the proposed method and indexes were then reanalyzed.

Geology and mining conditions
Huating Coal Mine, a typical burst-prone coal mine, is located in the southeast of Gansu Province, China. The mining depth of No.5 coal seam in Huating Coal Mine is 400 $ 800 m, with the dip angle being between 1 $ 15 . Since the average thickness of coal seam is nearly 37 m, the coal seam is divided into three layers with an inclined slicing method to excavate from the top layer to bottom. Each layer is nearly 13 m thick and the longwall top coal caving (LTCC) method is adopted. During longwall panel retreat, the first mining height is 5 m, while the remaining 8 m coal seam depends on roof pressure to fall (seen in Figure 1(b)).
The case study site, longwall (LW) 250106-1, is the upper stratified panel of No.5 coal seam with a length of 2000 m and a width of 200 m. As depicted in Figure 1(a), the east of LW250106-1 is the goaf of LW250104-1, and the west LW250107-1. There are 6 m rib block coal pillars between longwall panels to isolate the goaf. The stratigraphic column of LW250106-1 is shown in Figure 1(b).

Seismic monitoring system
The seismic monitoring system in coal mines can help calculate the location, energy, and focal mechanism of seismic events, invert coal and rock stress state, and evaluate the coal burst risk. (Cao et al. 2015). Since 2008, The Huating Coal Mine was equipped with a seismic monitoring system 'SOS' developed by the Central Institute of Mining in Poland to continuously monitor seismicity related to mining activities. As shown in Figure 2, aiming to monitor the seismicity of LW250106-1, two types of geophones, namely entry geophones and remote geophones, were installed underground. The entry geophones 1#, 7#, 8#, 9# and 16# were installed in the longwall entries and moved regularly with the longwall retreat, and the other remote geophones were placed in the roadway at different heights 2 $ 5 km from the floor to cover the monitoring area. Owing to the short distance between the longwall seismicity active zones and entry geophones, these geophones can capture relatively clear waveforms, which means they can usually be used to locate the seismic events. However, the occurrence of high-magnitude events (HMEs) would lead to some geophones exceeding the limit , especially for entry geophones. In this case, it can result in an underestimate of seismic energy. But luckily, the waveforms of HMEs can usually be captured by remote geophones. Owing to the long-distance propagation of waves after attenuation, the remote geophones will not overrun. Thereby the remote geophones can be auxiliary to calculate more accurate source energy, particularly for HMEs.

Coal bursts
Twenty-four coal bursts recorded during LW250106-1 retreat between 27 July 2016 and 10 November 2017 were used for this analysis. The twenty-four coal bursts with different seismic energy are described in Figure 3, and their key features are listed in Table 1. The seismic energy-triggered coal bursts varied from 1.89E þ 04 J to 1.35E þ 06 J, and most of them were located on the side of goaf adjacent to LW250401-1. The damaged length of coal burst extended from 5 m to 120 m, which was characterized by instantaneous floor heave, rib convergence, and roof fall. The shortest distance from the damaged area to coal burst events (NDDACBs) varied from 13 m to 920 m. Figure 4 shows the percentage distribution of coal bursts of different seismic energy levels and NDDACBs. As can be seen, the seismic energy of coal bursts between 10 4 J $10 5 J, 10 5 J $10 6 J, and 10 6 J $10 7 J accounted for 37.5%, 54.17%, and 8.33%, respectively, which means that coal bursts could be probably triggered when the seismic energy is larger than 10 4 J. As is shown in Figure 4  concluded that the closer HMEs are to the entry, the greater risk there is to trigger coal bursts. Figure 5 shows the distribution of damaged areas and the shortest distance from the damaged area to panel position. There were nearly three-quarters (18/24) of damaged areas located in tail entry, which means that the goaf side entry of the panel carries a higher risk. Furthermore, almost four-fifths (20/25) of damaged areas were located within in the 300 m range ahead of the panel. It signifies that the 300 m range ahead of the panel is high burst areas due to high abutment stress and strong disturbance seismic activities.  3. The statistics-based analysis of seismic energy attenuation characteristics

Description of seismic energy attenuation
A seismic event can be defined as a volume of a sudden inelastic strain of rock that radiates perceptible seismic waves (Alexey 2013). When a seismic event occurs, a stress drop and energy release happen on the ruptured plane. Particularly in underground mining, the fracture size of HMEs may be up to hundreds of meters (Cao 2009). So it is essential to consider the influence of fracture size on seismic events. The seismic wave propagation is characterized by sustaining inherence attenuation, geometric diffusion attenuation, and scattering effects. . In consideration of the effects of seismic propagation, the seismic waves propagating seismic energy E ij from a mining-induced seismic source epicenter j to point i could be expressed as follows in exponential attenuation form (Cao 2009):  where E j denotes the seismic energy of event j, R ij refers to the distance between event j and point i; and g is the attenuation coefficient used to describe the attenuation effects of three modes as mentioned above. Moreover, a seismic source is assumed to be characterized by globosity with a circular rupture surface. For further simplification, it is assumed that the seismic energy remains constant within fracture radius (source radius) r 0 and decays from the sphere surface. Therefore, the released strain energy from seismic sources is derived as follows: In this study, the radius of apparent volume, named apparent radius, is used to estimate the source radius on the ground that apparent volume is a scalar estimation of the volume of co-seismic inelastic deformation that radiates seismic waves, depending on the seismic moment and seismic energy which are independent of the source model (Mendecki 1996). The apparent volume is calculated as follows: where M 0 denotes seismic moment, G is shear modulus; and E s refers to seismic energy. The apparent radius is calculated as:

Attenuation coefficient distribution of different seismic energy
As is assumed, the seismic energy E can be completely converted into kinetic energy when a particle with mass m at some distance from the seismic source is accelerated to velocity v. It can be expressed as follows (Roberts and Brummer 1988;Daehnke and Roberts 2001): In addition, due to the attenuation of the seismic waves propagating process, the vibration velocity of the particle will attenuate correspondingly. However, owing to the complex underground environment, it is difficult to distinguish the attenuation effects of three modes as referred above . Hence this attenuation is also usually expressed in exponential form: (Kijko 1994): Combining Eqs.
(2), (5), and (6), the following can be inferred: In underground mining, the vibration velocity peak of a particle is also called peak particle velocity (PPV) , which can be measured by capturing the amplitude of a series of seismic wave signals by a geophone. Figure 6 depicts an actual capture process of PPV for a certain seismic event.
Furthermore, according to the relationship between the distance and the corresponding PPV, the PPV amplitude attenuation coefficient can be computed by fitting the data ) (seen in Figure 7). As shown in Figure 7, the PPV amplitude attenuation coefficient b of this seismic event is fitted as 0.0004713, thereby the energy attenuation coefficient g is twice as b, reaching up to 0.0009426. In this way, a total of 11578 seismic events during LW250106-1 retreat were introduced to fit b and further calculate g. Figure 8 donates the PPV attenuation coefficient statistic of different seismic energy levels, which vary from 2.7857 Â 10 À4 to 6.9839 Â 10 À4 . There is a negative correlation between b and seismic energy levels. It suggests that the higher the seismic energy is, the less attenuation there would be during seismic wave propagating. In other words, the HMEs will cause more coal burst risk not only because of the high seismic energy but less attenuation as well. In this sense, it should be noted that attenuation of seismic energy is a critical issue when assessing the coal burst risk or investigating the mechanism of seismic disaster. Table 2 shows the calculated g based on the PPV attenuation coefficient of different seismic energy levels, which will be used to further calculate the coal risk assessment indexes in the later section.

Attenuation response characteristics of seismic energy to different factors
Because seismic energy is the most fundamental parameter of seismic monitoring, it is widely used to evaluate the dynamic hazard. For example, in Poland, scholars put forward the criterion for quantitative evaluation of coal burst risk through seismic energy , which included: the maximum energy recorded by the seismic system in 24 hours (E max ) and the total seismic energy recorded by the seismic system every 5 m forward of the longwall face. According to these two indexes, the coal burst hazard level can be divided into none, weak, medium, and strong risk, respectively. Particularly the principles indicate that the monitoring area is a low hazard when E max > 10 4 J. However, it may not be feasible to measure the disturbance of seismic events accurately. As seen in Eq. (7), the PPV and seismic energy will be attenuating, which means that the residual seismic energy (RSE) which represents the energy of the actual action on the coal and rock mass after attenuation could be remarkably different. Therefore, to assess the burst risk caused by seismic events, it is necessary to investigate the characteristics of RSE.   Seismic energy level (J) Energy attenuation coefficient g 10 2 $10 3 J 1.3968 Â 10 -3 10 3 $10 4 J 1.2046 Â 10 -3 10 4 $10 5 J 9.1698 Â 10 -4 10 5 $10 6 J 7.2674 Â 10 -4 10 6 $10 7 J 5.5714 Â 10 -4 In this section, the attenuation response characteristics of seismic energy to different factors, including seismic energy, source radius, and energy attenuation coefficient, were investigated via the control variable method, which is helpful to understand the burst risk caused by RSE. Figure 9(a) shows the RSE distribution of different seismic energy level when g ¼ 0.0001, r 0 ¼ 10 m. (In the context aiming to investigate the effect of a single factor on RSE, the values of other factors are relatively conservative, e.g. g equal to 0.0001 and r 0 equal to 10 m rather than strictly referring to section 3.2.) The RSE with different seismic energy has experienced a similar attenuation process versus distance between the point and seismic spherical surface, but the magnitude of RSE is distinct-different. For example, if the seismic energy equals 1E7 J, the RSE will decrease to 1E4 J when propagating nearly 300 m. However, if the seismic energy equals 1E6 J, the seismic waves will only propagate about 100 m when the RSE decreases to 1E4 J. Hypothetically, the RSE ¼ 1E4 J could be regarded as the minimum risk threshold according to the Polish assessment principles as mentioned early, and the predicted area with risk is also totally distinct. Therefore, it suggests that RSE, instead of seismic energy, is a superior indictor for evaluating the coal burst risk or researching the burst mechanism. Figure 9(b) shows the RSE distribution of different source radius when g ¼ 0.0001, and seismic energy ¼ 10 6 J. The RSE with different source radius also displayed similar attenuation laws. Meanwhile, there is a positive correlation between RSE and source radius. Owing to the complexity of the mechanism of seismic sources in underground mining, the rupture types, and mechanisms of the sources are obviously distinct with similar seismic energy (Wojtecki et al. 2016). Therefore, it can be inferred that the larger the source radius is, the higher coal burst risk there is for the source with the same seismic energy. Figure 9(c) shows the RSE distribution of different energy attenuation coefficients when r 0 ¼ 10 m, and seismic energy ¼ 10 6 J. The RSE magnitude of different energy attenuation coefficients varies slightly for a certain distance between the point and seismic spherical surface. And the change of energy attenuation coefficient exerts a weak impact on the RSE. However, in fact, the change of source energy and source radius has a great influence on RSE, and in this case, different energy attenuation coefficients also exerts a great influence on energy. Moreover, as shown in Table 2, the reason why the energy attenuation coefficients of different seismic energy levels have inconspicuous differences may be that the seismic waves propagating in underground mines are faced with a similar geological environment, and accordingly they will have similar frequency and quality factors .
According to the analysis of attenuation response characteristics of RSE to different seismic energy, source radius, and energy attenuation coefficient, the following results can be obtained: The RSE increases with the increase of seismic energy and source radius, which means higher seismic energy and source radius will lead to higher coal burst risk. The RSE, rather than seismic energy, is a more precise indicator for evaluating the coal burst risk. Figure 9. The RSE distribution characteristics of different energy attenuation coefficients, source radiuses, or energy attenuation coefficients experimented by control variable method. (a) The RSE distribution of different seismic energy levels, when g ¼ 0.0001, r 0 ¼ 10m, (b) The RSE distribution of different source radiuses, when g ¼ 0.0001, and seismic energy ¼ 10 6 J, and(c) The RSE distribution of different energy attenuation coefficients, when r 0 ¼ 10m, and seismic energy ¼ 10 6 J.
Aiming to calculate the RSE accurately, a reasonable source radius and energy attenuation coefficient must be determined. In this paper, the apparent radius and statistical energy attenuation coefficient of different source energy levels are adopted, as mentioned in section 3.1.

Static and dynamic response and its induced coal burst process in underground mining
In underground mining, the extraction of solid coal leads to stress redistribution around the longwall panels, as described in Figure 10. The areas with abnormally high vertical stress are called abutment zones, and the abnormal vertical stress is known as abutment stress. Along the direction of the panel advance, there will be a mostly quasi-static loading process, and the coal mass close to the excavation boundary is fractured, and the vertical stress corresponds to the post-peak state post-peak residual stress (DE). It can be inferred that vertical stress can be divided into three stages from deeper solid coal to excavation boundary, namely, elastic stage (AB), the pre-peak plastic deformation stage (BD), and the post-peak strain-softening stage (DE). In this process, the primary fissure in coal mass will also propagate, converge, and connect to form macro-fractures or even slip-zones. Not only can it lead to the cumulative damage of coal, but also generate lots of seismic waves. As referred to in the literature by Cai et al. (2020), the relationship between damage level and distributional vertical stress is shown in Figure 10. Moreover, the continued mining process may induce roof fracturing or fault activation and associated seismic events, as depicted in Figure 10. Especially for vast layers of sandstone fracturing, it may generate seismic events up to 10 6 J . It could generate additional dynamic loads following the law of exponential attenuation ahead of the longwall panel. In this process, the coal burst is more likely to occur when the total stress is superimposed by static stress (i.e. abutment stress) and dynamic stress exceeds the critical stress for coal burst initiation (Dou et al. 2014;Li et al. 2015). The conceptual model can be established as follows (He et al. 2017): where r s denotes the static stress of coal and rock mass; r d is the dynamic stress cause by roof fracturing, fault slipping etc.; and r c refers to the critical stress condition of coal burst initiation. As can be known from the above statement, coal burst risk assessment or prediction needs to be carried out from three aspects: static stress, dynamic stress, and their superposition effect. However, as stress is a vector, it cannot simply be superimposed numerically. In addition, the dynamic stress field induced by seismic waves is extraordinarily complex. Hence there are few effective methods to carry out the inversion of the dynamic stress field. Fortunately, some strategies have been proposed to solve this problem. For example, Wang et al. (2020) proposed that passive velocity tomography (PVT) could be used to evaluate risks induced by static stress concentrations, and Dynamic Load Index (DLI) could express dynamic load disturbance intensity caused by seismic events quantitatively. Furthermore, by combining PVT and DLI, a dynamic and static load evaluation index could be constructed, which effectively assesses coal burst risks with the superposition effect of static and dynamic stress taken into consideration. Therefore, combined with static and dynamic load, an index can be feasibly constructed to indicate the degree of burst risks comprehensively.

Static intensity index based on damage mechanics
As described in Figure 10, the static stress redistribution and dynamic disturbance result in the damage of coal mass, and the damage level gradually increases to a critical damage point D F as the static stress evolves from pre-peak to post-peak state. Based on continuum damage mechanics (Lemaitre 1985), the static stress state transition can be described as follows: where r and e represent static stress and strain respectively; E denotes Elastic Modulus; and D refers to damage variable. In field implementation, the strain of coal mass can hardly ever be measured immediately. However, Benioff (1951) pointed out that the elastic strain rebound increment was proportional to the square root of the energy released after the earthquakes, which was later named as Benioff strain. Therefore, the Benioff strain can be used to estimate the accumulation strain of coal mass under loading. Furthermore, the research region could be discretized into grids and damage variable could be constructed by Benioff strain as well, which can be expressed as follows: where N denotes the number of seismic events; e ei represents the cumulative Benioff strain of the ith grid caused by all the seismic events; k is the ratio coefficient. Then e f can be calculated as follows (Cai et al. 2018): Where e ef denotes the maximum value of e e ; and the critical damage D F ¼0.95 corresponds to the ultimate damage state of coal and rock mass .
Furthermore, the seismic-based static stress can be expressed by substituting Eqs. (10) and (11) into Eq. (9): Where r si is the static stress of the ith grid. The stress levels at different positions in the research area may vary greatly, leading to the stress values of different nodes differing greatly and reducing visualization when cloud images are drawn. To overcome this obstacle, a new index, named 'Static Intensity Index' (SII), is introduced for better static stress visualization based on the normalization method. SII is defined as: Where SII i is the SII value at point i, r smax and r smin are the maximum and minimum value of r si :

Dynamic Intensity index based on energy concentration
The earthquake energy is initially stored in the medium as elastic strain within a certain volume. The energy is converted into kinetic energy of the rock mass around the fault after sliding along the fault plane. (Storcheus 2011). As there is no essential difference between seismic events and earthquakes, the seismic energy after attenuation (i.e. RSE) can also be stored as elastic strain in the coal and rock mass before damage, resulting in dynamic stress concentration. It will induce a coal burst when the relation in equation (8) is satisfied somewhere. Consequently, it is indispensable to evaluate the dynamic stress concentration caused by seismic events.
The conventional analytical method illustrates dynamic stress concentration in different energy levels with circles, which are convenient for intuitive and qualitative analysis, but not for quantitative judgment . Moreover, this method does not take the attenuation effect of energy into account. Therefore, the energy accumulation is proposed as follows: Where r di denotes the accumulated energy of the ith grid; and E ij is the seismic waves propagating seismic energy from a seismic source epicenter j to point i (i.e. RSE at the ith grid).
To improve visualization, the energy accumulation index named 'Dynamic Intensity Index' (DII) is normalized as: Where DII i is the DII value at point i, r dmax and r dmin are the maximum and minimum value of r di :

Risk assessment index linked with static and dynamic loading effects
Based on the superposition principle of dynamic and static stress of coal burst (Eq. (1)), the 'Risk Assessment Index' (RAI) integrating SII and DII is defined as follows: Where RAI i denotes the RAI value at point i, w 1 and w 2 are weights of SII i and DII i , respectively, and w 1 þ w 2 ¼ 1: In this section, the bias difference method is used to determine the weights: Where the prefix 'std' represents the standard deviation of the corresponding index. On this basis, if the standard deviation of the index is higher, it means that the index can reflect the situation more realistically, and thus it will occupy a higher weight. According to Eqs. (13), (14), and (16), the 'Static Intensity Index' (SII), 'Dynamic Intensity Index' (DII), or Risk Assessment Index (RAI) are distributed within the range from 0 to 1. The none, weak, moderate and strong risk of coal bursts can be classified quantitatively by the three indexes based on the index values ranging from 0 to 0.25, 0.25 to 0.50, 0.50 to 0.75, and 0.75 to 1. According to the various risk levels of coal bursts, corresponding control measures could be taken as summarized in Table 3. Figure 11 shows an example for calculating cumulative Benioff strain or RSE at any point on the grids. First, to plot the contour maps of RAI, the research area is discretized into grids and N seismic events are assumed to be detected in a given period. Second, the cumulative Benioff strain or RSE from all the events to point i can be derived by using the attenuation law. Third, every point on the grids can be implemented by this calculation process. Figure 12 shows the process for calculation and cloud map drawing of Static Intensity Index (SII), Dynamic Intensity Index (DII) and Risk Assessment Index (RAI). The seismic parameters (x, y, z, E j , r 0 , and g) are used in this process. First, the cumulative Benioff strain or RSE is calculated for each grid as described in Figure 11. Subsequently, each grid's SII, DII, and RAI index values can be obtained by Eqs. (10)-(18), and pixel interpolation is used to plot the contour map accordingly. Table 3. Relationship between the aforesaid indexes and coal burst risk levels and prevention measures. (after Dou and He (2007) Figure 13 shows the evolution of seismic event frequency, total daily energy, and maximum daily energy prior to the coal burst on 2017/11/01 ('11.01' coal burst). The figure shows that seismic event frequency and total daily energy remain a gradual upward trend during this period, and maximum daily energy exceeds 1E4 J in most cases. It indicates that the longwall has undergone noticeable changes in seismic energy release, thus enhancing coal burst risk. Figure 14 shows the distribution of seismic events in LW250106-1 in a 10-day time window before the coal burst which occurred on 2017/11/01. At this stage, the stress level of the longwall increases sharply as the longwall gradually advances toward the synclinal axis (seen Figure 1(a) and Figure 10(a)-(c)). Simultaneously, the seismic events gradually gather near the Face position, and the proportion of HMEs increases apparently. It also indicates a relatively high coal burst risk. However, due to the superposition and tanglesome distribution of seismic events, it can only be judged that the risk area is qualitatively based on abnormal concentration of seismic events; quantitative and accurate identification of risk areas cannot be obtained.

Risk-based assessment
To investigate the prediction efficiency of SII, DII and RAI for the HMEs and coal burst damage, the coal burst that occurred on 2017/11/01 ('11.01' coal burst) was reanalyzed based on the seismic data of a month. In this process, the SII and DII were first calculated based on seismic data in the 10-day time window and the 3-day time window, respectively. Second, combining SII and DII, the RAI was obtained by Eqs. (16)-(18). Third, the HMEs during the next ten days and the final coal burst (i.e. '11.01' coal burst) were used to evaluate the prediction efficiency of the three indexes.      (or damaged area) within the next 10 days. The figures show that the predicted risk areas based on SII were always located on the side adjacent to the goaf ahead of the panel. It means that the lateral overhanging roof of LW250401-1 goaf and the abnormal activities of overlying rock in goaf increase stress level and decrease the entry stability in the goaf side. Furthermore, the prediction accuracy for HMEs (Ratio of the number of HMEs in the above risk areas to the number not in the areas) was used to measure prediction efficiency. As shown in Figure 15(a)-(b), the prediction accuracy of HMEs was 43% and 69%, respectively. And Figure 15(c) shows that the damaged area falls into the strong risk area and the coal burst falls into the moderate risk area. It can be concluded that the SII can effectively predict the damaged area, but the prediction efficiency of HMEs is relatively low. Figure 16(a)-(c) depict the predicted risk areas based on DII and event distribution (or damaged area) within the next 10 days. As shown in Figure 16(a)-(b), the prediction accuracy of HMEs was 95% and 88%, respectively. In addition, Figure  16(c) shows that the damaged area fell into the strong risk area and the coal burst fell into the moderate risk area. It can be seen that the DII represents a pretty high prediction accuracy for the HMEs and damaged area. The predicted moderate risk area covered most areas within the 300 m range ahead of the longwall face. These areas were strongly disturbed by dynamic load caused by roof movement, crack propagation, etc. However, the risk area predicted by DII is so wide as to overestimate the coal burst risk, thus reducing the prediction efficiency to some extent, which is not conducive to the implementation of corresponding prevention and control measures as shown in Table 3. Figure 17(a)-(c) show the predicted risk areas based on RAI and event distribution (or damaged area) within the next 10 days. As shown in Figure 17(a)-(b), the prediction accuracy was 67% and 88%, respectively. In addition, Figure 17(c) indicates that the location of coal burst and damaged area was predicted effectively. The result suggests that the RAI not only improves the prediction efficiency of HMEs and risk area but also overcomes the shortcoming of overestimating risk areas, which also verifies Figure 18. Prediction accuracy of HME and damaged area with magnitude of completeness variation of studied coal burst cases. the necessity of coal burst assessment from both static and dynamic aspects. Furthermore, a series of controls such as coal mass pressure relief blasting, floor distressing borehole, and roof pre-blasting should be implemented on the predicted areas to eliminate or reduce the coal burst risk.

Discussion
As the RAI is fully dependent on seismicity, the seismic data pattern and quality would control the outcome. For daily assessing results, its quality should be quantified and the influencing factors should be figured out. The prediction accuracy of damaged area (P damage ) and HMEs (P HME ) is proposed here to evaluate result quality, which is expressed as follows: Where S RAI>0:5 is the moderate and above-predicted risk area based on RAI, S recorded denotes recorded damaged area; N RAI>0:5 represents the number of HMEs falling into moderate or above risk zone; and N recorded refers to the total number of HMEs recorded in the statistical time window. In addition, the magnitude of completeness m c , the minimum magnitude that seismic event can be reliably recorded, can be used to evaluate the sensitivity of seismic monitoring system. The lower m c indicates that more complete seismic events could be captured, and the calculation method of m c could be referred to the research of Wang et al. (2021). The m c was calculated to analyze the prediction results auxiliary in the context.
The P damage and P HME of 24 coal burst cases were reanalyzed in the same manner as in section 5 as shown in Figure 18. In the figure, the blank column means no HME occurred during the statistical period, and the following analysis results were obtained: The m c fluctuated greatly in different statistical periods. The P HME was nearly up to 60% when the mc was at a low level. Except two cases occurring in 2017/03/15 and 2017/09/03, there are extreme variations of the P damage , either 0% or 100%.
Seismic signals are recorded in spatial-temporal heterogeneous underground surroundings due to the interference of mining activities under complex geological conditions. Though geophones are regularly installed along with the excavation, the seismic networks are still limited by the layout of coal mines. As shown in Figure 18, the m c fluctuated greatly in different statistical periods, which means that great changes in seismic data integrity took place in different periods of excavation. It is noted that P HME is higher when m c is lower, which indicates good integrity of seismic data. Meanwhile, there are extreme variations of the P damage , either 0% or 100%. In most cases, the damaged area can be predicted by RAI when the m c is low. In the false assessments, either the m c is very high, or the damaged area is far from the face position. In both cases, there is no enough seismic data to analyze, which in turn leads to false assessments due to the restricted sensitivity of the monitoring system or low seismicity activeness. This demonstrates that the risk assessing method is applicable in a zone with high seismicity activeness. In those low seismicity activeness areas, some additional monitoring methods, such as drilling bits, electromagnetic radiation, acoustic emission, borehole detection, borehole pressure monitoring, and active velocity tomography, could be adopted to detect and monitor potential risks.

Conclusions
Seismic events occur with excavation in underground coal mines. It is of great significance to understand the attenuation characteristics of seismic energy for the mechanism of coal burst and its prediction. This paper systematically investigates seismic energy attenuation characteristics based on reliable statistics. Combined with the static and dynamic effects of underground mining, a new risk assessment method based on seismic energy attenuation calculation was proposed to assess coal burst risk. The effectiveness of this method was verified by back analysis based on the data from LW250106-1. The following conclusions can be drawn: 1. The seismic energy attenuation coefficient ranges from 2.7857 Â 10 À4 to 6.9839 Â 10 À4 and is negatively related to seismic energy level. This indicates that HMEs will cause more coal burst risk not only because of the high seismic energy but less attenuation as well. 2. The RSE is relatively correlated with seismic energy, source radius and energy attenuation coefficient. The higher seismic energy, source radius, and lower energy attenuation coefficient will result in higher coal burst risk. 3. The proposed SII based on damage mechanics could effectively predict the damaged area, but the prediction efficiency of HMEs is relatively low. The proposed DII based on energy concentration could availably predict damaged areas and HMEs, but it may lead to overestimation of coal burst risk. Combined with the advantages of SII and DII, RAI can efficiently predict the damaged area and HMEs. The case study demonstrates that the average prediction accuracy of HMEs based on RAI was up to 67% and 88%, respectively. 4. Sensitivity of monitoring system and seismicity activeness could be the two main restrictions leading to low accuracy or even false risk assessment, which implies that the proposed method could also have false judgments in the place where rare seismicity is monitored. The diversified auxiliary monitoring methods should be implemented in low seismicity activeness areas.

Availability of data and material
Not applicable.

Code availability
Not applicable.

Disclosure statement
No potential conflict of interest was reported by the authors.