Focal Mechanism and Source Parameters Analysis of Mining-Induced Earthquakes Based on Relative Moment Tensor Inversion

Mining-induced earthquakes (MIEs) in underground coal mines have been a common phenomenon that easily triggers rock bursts, but the mechanism is not understood clearly. This research investigates the laws of focal mechanism and source parameters based on focal mechanism and source parameters analysis of MIEs in three frequent rock burst areas. The relative moment tensor inversion (MTI) method was introduced, and the way to construct the inversion matrix was modified. The minimum ray and source number conditions were calculated, and an optimized identification criterion for source rupture type was proposed. Results show that the geological structure, stress environment, and source horizon influence the focal mechanism. The tensile type sources can distribute in the roof and coal seam, while the shear types are primarily located in the coal seam. In the typical fold structure area, the difference in source rupture strength and stress adjustment between tensile and shear types is negligible, while the disturbance scale of tensile types is distinct. The shear types have higher apparent volume and seismic moment in the deep buried fault area but lower source energy. The apparent stress of the tensile types is higher than that of the shear types, representing that the stress concentration still exists in the roof after the MIEs, but the stress near the faults could be effectively released. In the high-stress roadway pillar area, the primary fracture of the coal pillar easily produces a continuous shear rupture along the dominant stress direction under the extrusion of the roof and floor. The source parameters (except apparent stress) of shear types are higher than tensile types and have higher dynamic risk. The results contribute to expanding the understanding of rock burst mechanisms and guide MIEs’ prevention.


Introduction
Mining-induced earthquakes (MIEs) describe the physical-mechanical response of coal and rock mass to the normal consequence of operation in underground coal mines [1,2]. As the mining depth increases, geological engineering conditions have become more complicated [3]. Subsequently, the term rock burst triggered by MIEs has been revealed as a common dynamic disaster in underground mining, which usually leads to fatalities, equipment trouble, and roadway damage with the release of lots of elastic energy [4][5][6][7][8]. The mechanism of rock burst is still one of the least understood and challenging tasks of

Inversion Method of the Relative Moment Tensor for the Focal Mechanism of Mining-Induced Earthquakes
Based on the arrangement of seismic monitoring systems in coal mines, the focal mechanism and source parameters of MIEs can be calculated to evaluate the dynamic hazard [5]. Along with the elaboration of the source theory, people have realized that the source mechanical state at the rupture time is more complex than simply measuring it by a single force or a pair of staggered forces [42]. Therefore, the MTI is introduced to explain the rupture behavior of sources. For example, the response characteristics and rupture types of key strata during periodic tests in a Chinese coal mine were identified via the MTI [35].

Purpose of Focal Mechanism Inversion via MTI
The purpose of focal mechanism inversion is that people hope that the stress state at the time of source rupture can be grasped to analyze the mechanism of source rupture better. When the earthquake wavelength is much smaller than the distance between the source and the station in space and the period of earthquake wave is much longer than the rise time of the earthquake wave in temporal, the earthquake sources can be regarded as point sources [41]. For a point source, in the Cartesian coordinate system, the moment tensor of sources M can be expressed as: where M ij (i, j = 1, 2, 3) represents the value of the moment tensor and the first and second subscript represents the axial direction of the force couples and the direction in which it points. The nine different force couples of the moment tensor are depicted in Figure 1. Due to the moment tensor meeting the condition that angular momentum is conserved (e.g., M ij = M ji ), the M has only six independent elements [43]. The nine different force couples that make up the moment tensor components. (After Shearer [43]).
The relationship between the monitored displacement u n ij,k (x, t), Green's function G n ij,k (x, x 0 , t), and the moment tensor M jk can be expressed as Equation (2), which will be used to invert the moment tensor of the point source. u n ij,k (x, t) = G n ij,k (x, x 0 , t)M jk * S(t) i, j, k, n = 1, 2, 3 where the u n ij,k (x, t) means the displacement of the i-th source caused by the k-th seismic phase in the n-th direction of the j-th station, here k = 1, 2, or 3 indicates the waves excited from the P, SH, and SV seismic phase, respectively, and n = 1, 2, or 3 represents three orthogonal directions; G n ij,k (x, x 0 , t) represents the relationship between the unit force at time t and point x 0 and the displacement produced at point x; the symbol " * " is the convolution operation, and the S(t) represents the time function of source. It can be noticed that S(t) is replaced as an impulse function δ(t) due to the synchrosource hypothesis (i.e., the period of earthquake wave is much longer than the rise time of the earthquake wave in temporal.). Especially for MIEs (i.e., the medium and small earthquakes), the earthquake wave propagation time is approximately zero, so the impulse function δ(t) is nearly equal to 1 [44]. Therefore, Equation (2) for MIEs can be simplified as: u n ij,k (x, t) = G n ij,k (x, x 0 , t)M jk i, j, k, n = 1, 2, 3.
Furthermore, Equation (3) can be expressed in matrix form as: where u is the far-field displacement sequence matrix obtained from different stations, which means that only the displacement sequence obtained from stations far enough from the source Equation (4) can be constructed, and the minimum displacement condition is proposed as more than 500 m in underground mining [45]; G is the theoretical Green's function related to the velocity model of formation structure. It is widely accepted that the low-frequency spectral level Ω 0 can be regarded as the displacement u i , which can be calculated as [46]: where S D2 and S V2 are the squared spectral displacement and velocity integrals for both Pand S-waves. In addition, the direction of Ω 0 is judged according to the relative location between MIEs and stations, which can be expressed as [47]: Positive,The sensor located above the source, and the waveform initially moved up Negative,The sensor located above the source, and the waveform initially moved down Positive,The sensor located below the source, and the waveform initially moved down Negative,The sensor located below the source, and the waveform initially moved up (6) Therefore, the goal of MTI is acquiring the M according to Equation (4) and decomposing it to further explain the source rupture mechanism. The MTI can be divided into absolute, relative, and mixed inversion methods according to the treatment of Green's function [48]. Figure 2 describes the schematic diagram of the absolute MTI method. In this method, Green's function for each event is calculated based on the relationship between source location, station layout, and attenuation model of medium [41]. Hence, the inversion accuracy of the absolute MTI method mainly depends on Green's function. Although the hybrid moment tensor method reduces the influence of noise through repeated iteration, the inversion theory is the same as that of the absolute method. However, due to the complex geological environments and possible lateral inhomogeneities, it is always difficult to calculate Green's function accurately [48]. The relative MTI method was adopted to solve this problem since Green's function for each event was not calculated in the process.  Figure 3 describes the schematic diagram of the relative MTI method. In this method, they assume that the sources of the cluster experience a similar ray-path or propagation process, and the Green's function of the sources from the same source region is common. Generally, a reference event whose radiation pattern is known is used to estimate Green's function [49][50][51]. However, it is challenging to find a perfectly proper reference. Luckily, the relative MTI method proposed by Dahm [52] can overcome this shortcoming, i.e., the reference source is not needed anymore. Green's functions are simplified in his method as linear wave propagation terms and radiation patterns based on the ray-theoretical limit. The linear wave propagation terms are eliminated when radiation patterns are known-thereby avoiding the explicit calculation of Green's function. It is emphasized that the relative MTI method has two significant advantages over the absolute one: (1) it avoids solving complex Green's function and greatly improves the inversion accuracy; (2) the inversion efficiency is greatly improved by using source cluster compared with single-source inversion in absolute MTI. However, when introducing the relative MTI in the source mechanism of MIEs, several difficulties must be solved: (1) how to select the appropriate source cluster is necessary since MIEs are randomly distributed in temporal and space, i.e., the minimum ray and source number conditions should be calculated; (2) the construction of the inversion matrix needs to be optimized under the background of the stations in coal mine move along with the excavation in contrast to Dahm's; (3) the constructed of the inversion matrix, which may be a large sparse matrix, needs an effective solution method. The detailed description of the relative method inversion matrix construction and solution method can be seen in Appendices A and B, respectively.

Focal Mechanism Solution of MIEs
The focal mechanism solution is a prerequisite for analyzing the seismic source rupture mechanism by moment tensor theory. In detail, the rupture type and the occurrence of the rupture surface (i.e., azimuth, dip and slip angle) can be determined by decomposing the moment tensor. In addition, many applications of moment tensor theory are derived from these parameters. For instance, Šílený and Milev [53] distinguished roof caving, roof rupture, coal pillar instability, and low subduction reverse fault type fracture modes by inverting the focal mechanism of five mining-induced earthquakes in the Driefontein gold mine in South Africa. Moreover, several typical applications of MTI of mining-induced earthquakes have been developed, including the inversion of the rupture propagation process in hydraulic fracturing [53,54], the guidance of hydraulic fracturing parameter design [55], mining-induced stress field inversion, [56,57] etc. In this process, the effectiveness of inversion analysis depends on the suitability of the focal mechanism solution method.

Solution of Source Rupture Occurrence
According to the source rupture theory, a general dislocation model (seen in Figure 4) can describe the source rupture occurrence. The relationship between moment tensor component and rupture surface occurrence is as follows: M 11 = −κ cos α(sin 2φ sin θ cos γ + sin 2 φ sin 2θ sin γ) + sin α(1 + 2κ sin 2 φ sin 2 θ) M 12 = +κ cos α(cos 2φ sin θ cos γ + 0.5 sin 2φ sin 2θ sin γ) −κ sin α sin 2φ sin 2 θ M 13 = −κ cos α(cos φ cos θ cos γ + sin φ cos 2θ sin γ) +κ sin α sin φ sin 2θ M 22 = +κ cos α(sin 2φ sin θ cos γ − cos 2 φ sin 2θ sin γ) + sin α(1 + 2κ cos 2 φ sin 2 θ) M 23 = −κ cos α(sin φ cos θ cos γ − cos φ cos 2θ sin γ) −κ sin α cos φ sin 2θ M 33 = +κ cos α sin 2θ sin γ + sin α(1 + 2κ cos 2 θ) where M ij = M ij /M 0 (i, j = 1 to 3); M ij represents the moment tensor without the isotropic components; M 0 represents the scalar seismic moment (formula for the calculation is shown in Appendix A); κ represents the ratio of Lame constant at the rupture surface (i.e., κ = λ/µ); λ and µ represent the first and second Lame constant, respectively; φ, θ, and γ represent the strike, dip, and slip angle, respectively; α represents the dislocation angle.  [59]. N represents the north direction; D t and D s represent the shear and vertical displacement vector, respectively; φ, θ, and γ represent the strike, dip, and slip angle, respectively; α represents the dislocation angle. In this fault model, the initial fault gap is not taken into account. Equation (7) contains six nonlinear equations with five unknowns. Therefore, the problem of solving these equations can be transformed into searching for the best position in the space to make it closest to the target in the five-dimensional space. Meanwhile, the five unknowns have a clear range, which makes the range of five-dimensional space limited, and reduces the difficulty of solving the problem and the possibility of multiple solutions to a certain extent. Therefore, Particle swarm optimization (PSO) [58] was adopted here to solve this question.

Identification of Source Rupture Type
For identification of the rupture type, there are two widely accepted approaches: proportion of moment tensor components-and dislocation angle of rupture surface (seen in Figure 4)-based. For the former, the moment tensor component could be divided into a double-couple (M DC ), isotropic component (M ISO ), and a compensated linear vector dipole component (M CLVD ), whose proportion adds up to 100% (i.e., M DC % + M ISO % + M CLVD % = 100%) [43]. Consequently, [60] proposed to quantify the rupture type by the proportion of shear components (M DC %), in which M DC % ≥ 60% is defined as shear type rupture, M DC % ≤ 40% is defined as tensile type rupture, and 40% < M DC %<60% is a mixed-type rupture. For the latter, the dislocation angle of the rupture surface was proposed by Vavrycuk [61] as the basis for judging the rupture type. However, Ohtsu and Vavrycuk proposed the corresponding discriminant criteria for acoustic emission events and natural earthquakes, respectively. Due to the complexity and diversity of MIEs, if these criteria are applied to the classification of rupture type of MIEs, there may be a bias in the accurate identification of rupture type. Therefore, which parameters and ranges apply to the identification of MIEs' rupture types should be further investigated.
In this section, to obtain the appropriate criterion of source rupture type for MIEs, a parametric discussion about the proportion of shear ratio in moment tensor components (M DC %) and dislocation angle of rupture surface was conducted.
In source rupture theory, both the moment tensor components and occurrence of the rupture surface are related to the Lame constant of the medium, where the Lame constant is the material parameter characterizing the deformation of the medium in the strain-stress relationship (the detailed information is shown in Equation (7)). Therefore, it is indispensable to discuss the relationship between the M DC % and dislocation angle (α) under the different proportions of the Lame constant (i.e., different µ/λ). Figure 5 depicts the relationship between M DC % and cos α for different µ/λ. The following results can be obtained:

•
There is a positive correlation between the M DC % and cos α, which means that the two indexes are consistent in identifying the source rupture type, but both need a reasonable range of identification.

•
There is a nonlinear correlation between them, indicating that the shear displacement is still dominant when the M DC % is relatively low (M DC % ∈ [20,40]). It reminds us that there may be a bias on those sources characterizing low M DC % but high shear displacement.

•
Due to the nonlinear relationship between M DC % and cos α the specific interval of both M DC % and cos α should change with µ/λ.
It follows from the above statement that when the research object is a natural tectonic earthquake or rock specimen with pre-existed cracks, it is relatively efficient to identify the rupture type by M DC %. However, MIEs may have more complex rupture mechanisms, including shear, tension, compression, implosion, and other fractures [62]. It can be noticed that the Lame constants must be taken into account when identifying source rupture types based on not only M DC % but dislocation angle. The Lame constants are too abstract to apply to the rock mechanism. Therefore, the uniaxial compressive, tensile, and shear strength widely used in rock mechanics are used to replace Lame constants [45]. For extended-source events (dislocation angle larger than zero), the dislocation forms are shear dislocation and tensile dislocation, and the corresponding detection criteria for dislocation angle are as follows: For the shrinking-source event (the dislocation angle is less than zero), the dislocation form is shear dislocation and compression dislocation, and the corresponding detection standard of the dislocation angle is as follows: where F n , F s , and F t are the sources' uniaxial compressive, tensile, and shear strength. Suppose the mechanical properties of sources are not measured. In that case, the hypothetical relationship among them (i.e., the compressive strength of coal and rock mass is three times the shear strength and ten times the tensile strength) is proposed [45]. Consequently, the source rupture criterion combined with dislocation angle and Lame constants can be expressed as: It is emphasized that if the mechanical properties of rock near the source have been tested, the criteria for identification of rupture types for MIEs should be constructed according to Equations (8) and (9). Figure 6 shows the flowchart for implementing the method to apply. First, we should confirm the sources to be retrieved and the corresponding station information, including the three-dimensional coordinates of the sources, the station coordinates, and the wave velocity near the stations. Second, we should judge whether the sources satisfy the far-field condition. If meeting the far-field condition, the corresponding ray between source and station is valid; if not, those rays should be eliminated. Third, we should judge whether the screened conditions of source and ray numbers (seen in Equation (27)) are satisfied.

Guidelines on the Implementation
Provided the two conditions are met, the inversion matrix can be constructed; if not, we should increase the number of sources involved in the inversion return to the first step. After constructing the inversion matrix, we can solve it by Equations (28)-(31) and further decompose the moment tensor. Furthermore, the rupture occurrence (using Equation (7)), source rupture types (using Equation (10)), and beach balls can be obtained.

Response Laws of Focal Mechanism and Source Parameters (RLFMSPs) in Frequent Rock Burst Area
The focal mechanism can reveal the type and source of rupture occurrence. However, in the process of monitoring hazards induced by MIEs, the intensity information of source rupture, the scale of disturbance, and the level of stress change should be further understood to conduct targeted prevention. In addition, to further reveal the mechanism of MIEs and improve the ability to identify catastrophic events, the relationship between focal mechanism and source parameters should be analyzed in detail. Therefore, in this section, we compare and analyze the response laws of focal mechanism and source parameters (abbreviated as 'RLFMSPs' in this paper) in three different frequent rock burst mining areas in China. Therein, the Brune model [63] was used to calculate the source parameters, including scalar seismic moment, source radius, apparent volume, stress drop, and apparent stress. The physical meanings and formula for the calculation are listed in Appendix C.

A Case of RLFMSPs in Fold Structure Area
Huating Coal Mine, a typical burst-prone coal mine, is located in the southeast of Gansu Province, China. The first case study site, longwall (LW) 250105−1, is the upper stratified panel of No.5 coal seam with a length of 2000 m and a width of 200 m. There are 6 m rib block coal pillars between longwall panels to isolate the goaf. The mining depth of No.5 coal seam in LW250105−1 is 500~700 m, and the dip angle is between 1~5 • . Figure 7a shows that the heterogeneous fold structure leads to a complex stress environment and higher local horizontal stress. The west of LW250105−1 is the goaf of LW250103−1 (seen in Figure 7c), and the east is solid coal mass. There are 6 m rib block coal pillars between longwall panels to isolate the goaf, and the stratigraphic column of LW250105−1 is shown in Figure 7b. longwall panels to isolate the goaf, and the stratigraphic column of LW250105-1 is shown in Figure 7b.
Since 2008, The Huating Coal Mine has been equipped with a seismic monitoring system "SOS" developed by the Central Institute of Mining in Poland to monitor seismicity related to mining activities continuously. As shown in Figure 7c, aiming to monitor the seismicity of LW250105-1, two types of stations are installed underground, i.e., entry stations and remote stations. The entry stations 1#, 2#, 7#, 13#, and 16# were installed in the longwall entries and moved regularly with the longwall retreat, and the other remote stations were placed in the roadway at different heights 2~5 km from the floor to cover the monitoring area. It should be emphasized that as a similar seismic system is used to monitor MIEs in the mining area in the later section, more details about the system and layouts are not mentioned.   Since 2008, The Huating Coal Mine has been equipped with a seismic monitoring system "SOS" developed by the Central Institute of Mining in Poland to monitor seismicity related to mining activities continuously. As shown in Figure 7c, aiming to monitor the seismicity of LW250105−1, two types of stations are installed underground, i.e., entry stations and remote stations. The entry stations 1#, 2#, 7#, 13#, and 16# were installed in the longwall entries and moved regularly with the longwall retreat, and the other remote stations were placed in the roadway at different heights 2~5 km from the floor to cover the monitoring area. It should be emphasized that as a similar seismic system is used to monitor MIEs in the mining area in the later section, more details about the system and layouts are not mentioned.
There were 34 strong MIEs that occurred in LW250105−1 during the retreat from March 2014 to May 2015. The typical seismic waves of strong tensile type MIEs in LW250105−1 are depicted in Figure 8. We can see that the strong MIE caused several channels (such as channels 1#, 2#, 13#, and 16#) to exceed the limit, and the waveform detected by other stations far from the source lasted longer and attenuated slowly, which indicates that the strong MIE is characterized by strong high disturbance. In addition, the signal has poor P, SH, and SV amplitudes in some panels, e.g., panels 1#, 13#, and 16# in Figure 7. However, in the moment tensor inversion, these panels are not used, and only panels 2#, 3#, 4#, 5#, 6#, and 8# were used. Figure 9a shows the source location and focal mechanism solution of 34 strong MIEs during the LW250105−1 retreat. There were 20 tensile type ruptures and 14 shear-type ruptures, indicating that the tensile-type ruptures were dominant among the 34 MIEs investigated. Moreover, the tensile rupture is mainly located in the roof and coal seam, while the shear rupture is primarily located near the coal seam. It can be seen that the separation and tensile rupture of the high and low roofs are the main reasons for the strong tensile type MIEs, and the buckling and shear failure of the floor under high-level tectonic stress is the main factor of the strong shear-type MIEs triggered by complex fold structure and mining disturbance.  Figure 9b depicts the source rupture occurrence of 34 strong MIEs. The dip angle of most sources with tensile type rupture is less than 30 • , while the dip angle of those with shear type rupture is more than 30 • , which indicates that strong MIEs are closely related to the rupture type. In addition, the dip angle of the rupture plane located on the floor is more than 60 degrees, and the tensile type sources on the roof have more rupture dip angles than those located on the floor, which indicates that the occurrence of the rupture plane of the strong MIEs is closely related to the source horizon. Figure 9c displays the focal mechanism and response characteristics of the focal parameters of 34 strong MIEs. There are a few differences in source energy between the tensile-type and shear-type ruptures, which indicates almost no difference between the twotype ruptures in source intensity. Although both the tensile-type and shear-type sources have relatively close source radiuses, the apparent volume of the tensile types is 1.28 times that of the shear types. It showed us that the tensile types in the fold structure area have a larger disturbance scale. The apparent stress difference between tensile-type and shear-type ruptures is negligible, while the stress drop of shear types is about 10% higher than that of tensile types. It can be concluded that the rupture type of strong MIEs has little influence on the stress adjustment at the source.

A Case of RLFMSPs in Deep Buried Fault Area
The second case study site, LW3302, is located in the middle of Shandong Province, China. The LW3302 is fairly deep at about 1200 m underground with high in situ stress. In addition, as shown in Figure 10, a fault with a drop ranging from 1.4 m to 3.0 m and a dip of 60 degrees extended across the middle of the LW3302. During the LW3302 retreat between June and July 2015, 35 strong MIEs were detected near the LW3302, and seven rock bursts were induced in this process. The rock burst distribution is shown in Figure 10. It can be inferred that the existence of the fault not only led to the frequent MIEs and induced higher dynamic disaster but only high in situ stress. Figure 11a shows the source location and focal mechanism solution of 35 strong MIEs during the LW3302 retreat between June and July 2014. The sources are mainly distributed in coal seams or low-level roofs on elevation, and parts are located on the high-level roof. Among them, the sources located in the roof are judged as tensile-type ruptures, while the sources that occurred in coal are mainly shear-type ruptures. It shows that when mining near the deep buried fault area, the strong mining disturbance of the longwall is easy to induce fault activation, and the coal and rock mass near the fault is easier to produce slip rupture, and the shear component is dominant in the moment tensor. Meanwhile, fault activation can easily lead to roof breaking and instability characterized by tensile type rupture.   It can be seen that the rupture azimuth of the sources is about 150 • or 300 • and approximately parallel to the mining direction of the LW3302. In addition, the dip angle of most sources is less than 40 • , and only a few source dip angles are larger than 60 • , which indicates that the occurrence of MIEs is closely related to the mining direction of longwall face and fault orientation during mining in deep fault areas. Figure 11c shows the focal mechanism and response characteristics of focal parameters of 35 strong MIEs. The source energy of tensile-type sources is significantly higher than that of shear-type sources, while the seismic moment of shear-type sources is higher than that of tensile-type sources. Compared with the shear-type sources, the source energy radiating from sources in the roof will experience less attenuation and lead to higher dynamic risk. Furthermore, as the coal and rock mass in the fault area is relatively fractured, the fault will produce larger shear dislocation, which results in the seismic moment of shear types being higher than that of tensile types.

A Case of RLFMSPs in High-Stress Roadway Pillar Area
The third case study site, the roadway pillar area, is located in the south-middle of Shanxi Province, China. The operation target is the No.4 coal seam with an average thickness of 9.43 m and a buried depth of 800~1000 m. Due to the large buried depth, fold structure, and large goaf area around the roadways, special static stress is applied on the roadway pillars and leads to obvious stress concentration. Fifty strong MIEs were monitored from August to October 2018, and there were five rock bursts being induced (seen in Figure 12). As seen in Figures 12 and 13a, the coal pillars of the roadway bear a high abutment pressure due to the goaf area around the roadway. In addition, the repeated disturbance of the longwall mining leads to the development of primary fractures in the coal pillars, leading to a decrease in the stability of the coal pillars. Figure 13a also shows the source location and focal mechanism solution of 50 strong MIEs in the roadway area of the Gaojiapu Coal Mine from August to October 2018. In the 50 strong MIEs, shear-type sources accounted for 46 times, and tensile-type sources only accounted for 4 times. In addition, the horizontal projection of the five rock bursts is located in the coal pillar area of the roadways, vertically located in the coal seam or adjacent to the coal seam. It is suggested that the coal pillar is prone to shear instability along the pre-existing primary fracture when subjected to high-stress loading. As seen in Figure 13b, the azimuth of 50 strong MIEs was mostly concentrated between 120~150 • and 310~330 • , and dip angles were mostly less than 60 • . Combined with the geological condition of the roadways, it was seen that coal pillars were easily subjected to continuous shear failure along the dominant direction of stress under the extrusion of roof and floor under the high static load stress (large buried depth, isolated coal pillar, etc.) in coal pillar area and dense layout of roadways. Moreover, as shown in Figure 13c, the average value of source parameters except for the apparent stress of shear-type sources are higher than that of tensile types. It can be concluded that the rock burst risk of shear-type MIEs in the high-stress roadway is significantly higher than that of tensile types.
It follows from the above statement that there is a significant difference in the focal mechanism and source parameters, where the geological structure, stress environment, and focal layer have a prominent influence. These response characteristics provide a basis for the accurate identification of strong MIEs and dynamic risk prediction in underground coal mines.

Conclusions
Mining-induced earthquakes are one of the most common dynamic phenomena in underground coal mines. It is of great significance to understand the response characteristics of focal mechanisms and source parameters of the MIEs to predict and control dynamic disasters. An inversion method of relative moment tensor was proposed to obtain the focal mechanism solution. An inversion matrix was constructed based on the source group, and the identification method of source rupture type was optimized. The response laws of focal mechanism and source parameters of mining-induced earthquakes in three frequent rock burst areas were analyzed. There was a significant difference between MIEs in different geological and engineering backgrounds. Three main conclusions were drawn as follows:

•
In the typical fold structure area, the tensile-type sources are located in the roof and coal seam, while the shear types are primarily located near the coal seam. In addition, there is a close relationship between the occurrence of the rupture surface and the source elevation. The difference in source rupture strength and stress adjustment between tensile types and shear types is negligible, while the disturbance scale of tensile types is distinct.

•
In the deep buried fault area, the rupture of sources in the roof and coal are dominant tensile and shear types, respectively. The tensile types have higher source energy but lower seismic moment than shear types. The apparent volume of shear types is about three times for tensile types. The apparent stress of the tensile types is higher than that of the shear types, representing that the stress concentration still exists in the roof after the MIEs, but the stress near the faults could be effectively released.

•
In the high-stress roadway pillar area, the primary fracture of coal pillars easily produces continuous shear ruptures along the dominant stress direction under the extrusion of roofs and floors, which leads to a dominated shear-type rupture. The source parameters (except apparent stress) of shear types are higher than that of tensile types and show a higher dynamic risk. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.
The basic principle of relative MTI is that the linear wave propagation terms can be eliminated. For example, provided that two sources in the same region are simultaneously monitored when each source is observed by m stations, and they have the same "I-term", thereby "I-term" can be eliminated [52]: Equation (A7) can be further converted into an equivalent form: . . . . . .
To avoid a homogeneous system, we adopted an additional linear condition from Dahm [52], namely, When the whole inversion process is not limited to using two sources in the source cluster but selecting n sources observed by m stations, the inversion process can be further evolved into the following form: where: u n p1 a 11,q u n p1 a 12,q · · · u n p1 a 16,q u n p2 a 21,q u n p2 a 22,q · · · u n p2 a 26,q . . . . . . . . . . . . u n pm a m1,q u n pm a m2,q · · · u n pm a m6,q       p, q = 1, 2, 3 . . . n Equation (A10) can be converted to the matrix form as: R = GS . Where R represents a composite column matrix consisting of row zeros and a non-zero number; G represents the matrix describing the ray propagation effect between sources and stations; S represents the column matrix composed of the moment tensor components in the source cluster.
It follows from Equation (A1) to Equation (A11) that all sources in the source cluster must be monitored by several common stations during the construction of the matrix G. However, since the relatively small scale of MIEs and the complex mining surroundings, even for a fixed source cluster, it is difficult to ensure that the fixed numbered stations can record all sources. To solve this problem, we improve the way to construct matrix G to make it easier to construct under the background of uncertain fixed stations. Therefore, if the way to construct matrix G of Equation (A11) is transformed from sources to stations, Equation (A12) can be obtained: where G m represents the moment tensor inversion component matrix of all sources received by the m-th station with n m (n m −1) 2 lines and 6ncolumns, here n m represents the total number of sources received by the m-th station. In detail, the k-th line is shown as Equation (A13): αqk a 11,q u n αqk a 12,q · · · u n αqk a 16,q (A13) where q indicates the station number; α and β represent the source numbers received by station q, which is the common number within the whole source cluster; 0 qα , 0 qαβ and 0 qβ are single-row zero matrices, whose occupied lines are 6(α − 1), 6(α + β − 1), and 6(ω − β), respectively, herein ω indicates the total number of sources; similar to Equation (A13), B αβ represents the relative moment tensor relationship between two sources.

Inversion Matrix Solution of the Relative Method
In the relative MTI process, the selection of source number and ray number is the prerequisite for constructing an inversion matrix aiming to make sure the inversion matrix is solvable. Therefore, in this section, the minimum station and ray number conditions are derived according to the matrix solution requirements and the station layout background of coal mines. In addition, since the constructed matrix is likely to be a large sparse matrix, a suitable solution method is introduced here.
(1) Calculation of Ray and Source Number Conditions Provided that q stations and θ sources are involved in the inversion process; the number of sources detected by each station is x i (i = 1, 2, 3, · · · q;x i ≥ 2); the total number of rays received by q stations is x i > 4θ, according to the common P-wave first arrival localization method). Since θ sources have 6θ moment tensor components (i.e., 6θ unknowns), the line number of the constructed matrix according to Equation (A12) must be greater than 6θ, which can be expressed as: In principle, when the source and ray number meet the conditions described in Equation (A14), the minimum amount of data required by the relative MTI method can be confirmed. However, the ray distribution may be heterogeneous due to the complex underground surroundings, and the reliability of the monitoring system is relatively low for the stations with less ray propagation. Therefore, the accuracy of matrix construction should be further improved. Aiming to reduce the deviation of inversion caused by errors of individual stations, the ray number of a single station should be satisfied [45]: At the same time, the ray number of a single station must not be greater than the source number: Combining Equations (A14)-(A16), the minimum ray and source number conditions can be obtained: (2) Inversion matrix solution method Since the inversion matrix R = GS (G with R rows and C lines) is a large sparse matrix that has more rows than columns, it is extremely difficult to solve it by inverting it, and a flexible transformation is necessary. By multiplying the left and right sides of G transpose, the R= GS can be transformed as: Therefore, the solving inverse matrix of G −1 can be transformed into solving inverse matrix of the square matrix G T G. Furthermore, the singular value decomposition method can be used to solve the inverse problem. According to Lanczos's decomposition theory [64], matrix G can be decomposed into the following form: where V is a C × C dimensional matrix composed of orthogonal eigenvectors of a square matrix GG T ; E is a square matrix with a dimension of C × C, which consists of the eigenvalues of G; V is a C × C matrix consisting of orthogonal eigenvectors of a square matrix GG T . Since V is an orthogonal matrix (i.e., V T = V −1 ), thereby Equation (A19) can be further simplified as: Therefore, the inverse matrix problem is transformed into the eigenvector and eigenvalue problem of the square matrix G T G, which is expressed as: Appendix C Table A1. Common source parameters and physical meanings.

Source Parameters Parameter Meaning Formula for the Calculation Parameter Information
Corner frequency Preliminarily judge the size of the mining-induced earthquake scale; the larger the scale, the richer the low-frequency components of the spectrum.
The frequency corresponding to the intersection of the progressive high frequency trend and low frequency level of the source amplitude spectrum in logarithmic coordinates. ρ = medium density in source area; c = wave velocity in source area; R = the distance between source and station; Ω 0 = level of spectrum at low frequencies; F = radiation pattern (P-waves = 0.52, S-waves = 0.63); k = constant for the Brune model