Seismic source analysis inversion system based on mine earthquake monitoring

Abstract

The invention relates to the technical field of seismic source analysis and inversion, and discloses a seismic source analysis and inversion system based on mine earthquake monitoring, which is characterized in that energy parameters are generated through waveform feature extraction, mine earthquake event identification and accurate positioning are realized by using a space-time aggregation rule, and seismic source mechanism parameters are inverted in combination with a seismic wave initial motion parameter and amplitude relationship. On the basis, a regional stress field initial inversion model is constructed, a stress gradient is calculated in combination with a working face track and stratum information, a principal stress rotation composite index is determined through interlayer mechanical difference and a spatial change rate, and self-adaptive deformable partitions are generated. And finally, secondary stress field inversion is carried out in the subareas, the corrected maximum principal stress direction is output, high-precision dynamic analysis of a mine earthquake mechanism and stress distribution is realized, and the mine disaster prediction and early warning capability is remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of source analysis and inversion technology, and more specifically, to a source analysis and inversion system based on mine seismic monitoring. Background Technology

[0002] In deep mining, the roof typically consists of multiple layers of rock, characterized by alternating hard and soft layers, uneven thickness, and significant differences in mechanical properties. As the working face advances, the original tectonic stress and mining-induced stress superimpose, causing the stress state of the roof to exhibit complex changes both spatially and temporally. Especially near the goaf and within the advanced influence zone, the principal stress direction often becomes unstable, exhibiting a clear dynamic rotation. This phenomenon is even more pronounced under composite roof conditions, with the mechanical differences between different layers further exacerbating the uncertainty of the principal stress direction.

[0003] Existing focal mechanism analysis and stress field inversion methods generally rely on fixed partitions, assuming that the principal stress directions are relatively consistent within a certain spatial range. However, when the principal stress direction of the roof changes rapidly over a short distance or within a short time, traditional methods will classify multiple events with significantly different orientations into the same partition, leading to obvious conflicts between mechanism solutions. Since the inversion process depends on the statistical results of these partitions, the final stress direction often deviates systematically from the actual situation, thus affecting the accurate assessment of mine dynamic disaster risks.

[0004] The fundamental reasons for the aforementioned deviations are threefold: First, mining disturbances cause the principal stress direction to rotate rapidly along the working face advance direction; second, differences in the mechanical properties between layers of the composite roof cause inconsistent changes in the principal stress direction at different layers, creating a gradient effect in the vertical direction; and third, existing methods fail to adjust for dynamic rotation, still using fixed spatial divisions. Consequently, the source mechanism zoning and stress inversion results inevitably become distorted under complex stress environments, making it difficult to meet the engineering requirements for dynamic monitoring of the mine stress field. Summary of the Invention

[0005] This invention provides a source analysis and inversion system based on mine seismic monitoring, which solves the technical problems mentioned in the background art.

[0006] This invention provides a source analysis and inversion system based on mine seismic monitoring, comprising: The waveform feature extraction module acquires continuous waveforms from mine seismic monitoring and generates continuous parameters characterizing waveform energy changes. The event identification and localization module generates mine seismic events based on the spatiotemporal aggregation patterns of continuous parameters, and determines the spatial location and occurrence time of the mine seismic events; The source parameter inversion module uses the initial motion parameters and amplitude relationships of seismic waves to inversely derive the source mechanism parameters of mining earthquake events. The initial stress field inversion module, based on the source mechanism parameters and the spatiotemporal distribution of mining earthquake events, inverts the direction of the maximum principal stress in the regional stress field; The stress gradient calculation module determines the first spatial direction by combining the mining trajectory of the working face and the second spatial direction by combining the stratigraphic layering information, and calculates the spatial variation rate of the maximum principal stress direction in the first and second spatial directions. The coupling index calculation module obtains the interlayer mechanical difference index based on the mechanical parameters and thickness information of each layer of the composite top plate, and determines the principal stress rotation composite index by combining the two spatial change rates and the interlayer mechanical difference index. The deformable partitioning module generates adaptive deformable partitions along the direction of maximum principal stress based on the principal stress rotation composite index, and allocates and weights mine seismic events according to their spatial position relationship with the partition centerline. The stress field fine inversion module performs secondary stress field inversion based on weighted seismic events within the deformable partition, and outputs the corrected direction of the maximum principal stress.

[0007] Furthermore, continuous waveforms from mine seismic monitoring are acquired, and continuous parameters characterizing waveform energy changes are generated, including: The three-component raw waveform sequences collected from each station are subjected to mean removal and bandpass filtering to obtain standardized waveform sequences; The short-time average energy continuity, long-time average energy continuity, and energy ratio continuity are calculated based on standardized waveform sequences, as follows: The short-time average energy continuity is a time series obtained by averaging the absolute values ​​of a standardized waveform sequence over a first preset time window. The long-term average energy continuity is a time series obtained by averaging the absolute values ​​of a standardized waveform sequence over a second preset time window. The length of the second preset time window is greater than the length of the first preset time window; The energy ratio continuous quantity is a time series obtained by comparing the short-time average energy continuous quantity with the long-time average energy continuous quantity at the same time point; The short-time average energy continuous quantity, long-time average energy continuous quantity, and energy ratio continuous quantity are aggregated according to the station number to form a continuous parameter set.

[0008] Furthermore, based on the spatiotemporal aggregation law of continuous parameters, mine seismic events are generated, and the spatial location and occurrence time of these events are determined, including: The spatial coordinates of each station are obtained to form a station spatial location set, and the parameters containing the seismic wave propagation velocity of each rock layer are obtained to form a propagation velocity parameter set. Based on the station spatial location set and the propagation velocity parameter set, the straight-line distance from the candidate source point to each station is calculated, and then the straight-line distance is divided by the seismic wave propagation velocity of the corresponding rock layer to obtain the theoretical propagation time. The energy ratio continuous parameters of each station are time-aligned according to the theoretical propagation time and summed to obtain the spatiotemporal energy aggregation function. Local maxima are searched within the three-dimensional neighborhood of the spatiotemporal energy aggregation function, with each local maximum corresponding to a seismic event. The time coordinates of the local maximums are used as the initial values ​​for the occurrence time of the seismic event, and the spatial coordinates are used as the initial values ​​for the spatial location of the seismic event. Based on the initial values ​​of occurrence time and spatial location, the time corresponding to the maximum value in the energy ratio continuity of each station is extracted as the arrival time observation within a fixed window centered on the theoretical propagation time; Construct the arrival residual for each station; where the arrival residual is the sum of the arrival observation time and the theoretical propagation time; by minimizing the sum of squared arrival residuals for all stations, the occurrence time and spatial location of the seismic event are obtained. Select any two earthquake events from the set of earthquake events whose time interval is less than a first set value or whose spatial distance is less than a second set value to form an earthquake event pair; calculate the difference between the arrival observation difference and the theoretical propagation time difference for each station for the earthquake event pair, as the double-difference residual; by minimizing the sum of squares of the double-difference residuals of all stations, solve for the spatial position correction of each earthquake event; add the corresponding spatial position correction to the spatial position to obtain the updated spatial position of the earthquake event.

[0009] Furthermore, by utilizing the initial motion parameters and amplitude relationships of seismic waves, the focal mechanism parameters of mining-related seismic events can be derived inversely, including: Based on the occurrence time and updated spatial location of the seismic event, the theoretical propagation time of the seismic waves to each station is calculated. In the standardized waveform sequence, a fixed length P-wave window and a S-wave window are cut out with the theoretical propagation time as the center. The positive and negative polarities of the first wave of the P-wave are extracted from the P-wave window as the initial motion polarity of the P-wave, and the peak amplitudes are extracted from the P-wave window and the S-wave window as the P-wave amplitude and the S-wave amplitude, respectively. The propagation distance is calculated based on the updated spatial location of the seismic event and the spatial location of the station. The medium attenuation coefficient is calculated based on the quality factor, seismic wave frequency and propagation velocity in the propagation velocity parameter set. The medium attenuation coefficient is inversely proportional to the quality factor, directly proportional to the frequency and inversely proportional to the propagation velocity. Geometric diffusion correction is performed on the P-wave and S-wave amplitudes using the propagation distance, with the correction method being the P-wave amplitude or S-wave amplitude divided by the propagation distance. Attenuation correction is then performed on the corrected P-wave or S-wave amplitudes using the medium attenuation coefficient and the propagation distance, with the correction method being the P-wave amplitude or S-wave amplitude multiplied by the exponential attenuation term of the product of the medium attenuation coefficient and the propagation distance. The amplitude ratio observation is obtained by comparing the S-wave amplitude after the second correction with the P-wave amplitude. An objective function for source mechanism inversion is constructed, which includes: a P-wave first motion polarity consistency term and an amplitude ratio residual term. The P-wave first motion polarity consistency term is the sum of the number of inconsistencies between the observed P-wave first motion polarity at each station and the P-wave first motion polarity predicted by the candidate source mechanism; the amplitude ratio residual term is the sum of the squares of the logarithmic differences between the observed amplitude ratio at each station and the amplitude ratio predicted by the candidate source mechanism. The objective function for source mechanism inversion is minimized using a grid search algorithm. The parameters of the grid search are: strike angle from 0 to 360 degrees, dip angle from 0 to 90 degrees, and slip angle from -90 degrees to +90 degrees, with a preset step size. The strike angle, dip angle, and slip angle of the seismic event are obtained by solving the algorithm. The minimum value of the objective function is used as the sum of squared residuals as a quality index of the source mechanism parameters. The strike angle, dip angle, slip angle, and sum of squared residuals are aggregated to form a set of source mechanism parameters.

[0010] Furthermore, based on the focal mechanism parameters and the spatiotemporal distribution of the mine-induced seismic event, the directions of the maximum principal stresses in the regional stress field are obtained through inversion, including: Multiple inversion windows are defined; the temporal and spatial ranges of the inversion windows are set based on the density of seismic events, with the temporal range not exceeding ten times the sampling period of the stations and the spatial range not less than the average spacing of the station distribution; events that occur within the temporal range and whose updated spatial location is within the aforementioned spatial range are selected from the seismic event set to form a subset of inversion events; For each seismic event in the inversion event subset, the fault plane normal vector and the slip direction vector are derived through three-dimensional geometric relationships based on the corresponding strike angle, dip angle, and slip angle. The fault plane normal vector is perpendicular to the fault plane, and the slip direction vector points along the fault plane towards the slip direction. Based on the sum of squared residuals of the seismic events, weighting coefficients are generated through an exponentially monotonically decreasing function. The exponentially monotonically decreasing function is the weighting coefficient being the power of the product of the negative sum of squared residuals of the natural exponential function and the preset scale parameter. An objective function is constructed using the stress tensor as the unknown. The objective function is the sum of the angle mismatch terms for each seismic event, weighted by a weighting coefficient. The angle mismatch term is the square of the angle between the sliding direction vector and the shear stress vector of the event. The shear stress vector is the remaining part after subtracting the normal component of the normal vector of the fault plane from the stress tensor. The objective function is minimized using the least squares optimization algorithm, and the stress tensor is obtained by solving it. The stress tensor is decomposed into three eigenvalues ​​and three corresponding eigenvectors. The eigenvector corresponding to the largest eigenvalue is selected as the direction of the maximum principal stress in the inversion window. The direction of the maximum principal stress is bound to the midpoint of the time range and the geometric center of the spatial range of the inversion window. The maximum principal stress direction field is formed by traversing all inversion windows.

[0011] Furthermore, by combining the mining trajectory of the working face to determine the first spatial direction and combining the stratigraphic stratification information to determine the second spatial direction, the spatial variation rate of the direction of maximum principal stress in the first and second spatial directions is calculated, including: For the geometric center of the spatial range of the inversion window, the nearest target trajectory point is searched in the mining trajectory set of the working face by calculating the Euclidean distance; the tangential vector is calculated by the coordinate difference between the target trajectory point and the previous and next adjacent trajectory points, and the first spatial direction is obtained by dividing the tangential vector by the modulus. The first spatial direction is the tangential direction of the working face advance. For the geometric center of the spatial range of the inversion window, the corresponding local bedding plane is extracted from the stratigraphic information set. The local bedding plane is described by the coordinates of multiple discrete points. The least squares fitting method is used to fit the discrete point coordinates to obtain the plane equation of the local bedding plane. The normal vector is determined based on the coefficients of the plane equation. The normal vector is divided by the modulus to obtain the second spatial direction, which is perpendicular to the local bedding plane. The angle between the projection of the direction of the maximum principal stress onto the horizontal plane and the north direction is taken as the direction angle, which is defined by a 180-degree period. Angle dewinding is applied to the direction angle sequence to eliminate jumps when crossing zero or 180 degrees. The step size of the inversion window in the first spatial direction and the step size in the second spatial direction are read. The step size is equal to the Euclidean distance between the geometric centers of the spatial extents of adjacent inversion windows in the corresponding directions. The rate of change of the orientation angle sequence is calculated using the central difference method. The rate of change of the inversion window is equal to the difference in orientation angle between the forward adjacent inversion window and the reverse adjacent inversion window divided by twice the step size. For the first and last inversion windows in all inversion windows, forward difference or backward difference is used for calculation. The rate of change of the first inversion window is the difference in orientation angle between the second inversion window and the first inversion window divided by the step size, and the rate of change of the last inversion window is the difference in orientation angle between the last inversion window and the second to last inversion window divided by the step size. The first spatial rate of change sequence and the second spatial rate of change sequence are formed.

[0012] Furthermore, based on the mechanical parameters and thickness information of each layer of the composite roof, an interlayer mechanical difference index is obtained. Combining the two spatial variation rates and the interlayer mechanical difference index, the principal stress rotation composite index is determined, including: Within the vertical range corresponding to the geometric center of the spatial range of each inversion window, the thickness and equivalent Young's modulus of each rock layer within the vertical range are read; the proportion of the thickness of each rock layer to the total thickness of all rock layers is calculated to obtain the proportional weight; the weighted average of the equivalent Young's modulus is calculated based on the proportional weight; where the weighted average is equal to the sum of the products of the equivalent Young's modulus of each rock layer and the corresponding proportional weight; the weighted standard deviation of the equivalent Young's modulus is calculated; where the weighted standard deviation is equal to the square root of the sum of the squares of the differences between the equivalent Young's modulus and the weighted average and the corresponding proportional weight; the ratio of the weighted standard deviation to the weighted average is used as the interlayer mechanical difference index, which is used to characterize the degree of dispersion of the interlayer mechanical properties of the composite roof; For each inversion window, the first spatial rate of change is multiplied by the first spatial step size to obtain the standard first spatial change; the second spatial rate of change is multiplied by the second spatial step size to obtain the standard second spatial change. The sum of squares of the standard first spatial variation and the standard second spatial variation is calculated, and the square root of the result is taken to obtain the comprehensive rotational intensity. The comprehensive rotational intensity is multiplied by the interlayer mechanical difference index to obtain the principal stress rotational composite index. The principal stress rotational composite index sequence is formed by traversing all inversion windows.

[0013] Furthermore, based on the principal stress rotation composite index, adaptive deformable partitions are generated along the direction of maximum principal stress, and seismic events are allocated and weighted according to their spatial relationship with the partition centerline, including: For each inversion window, read the direction of the maximum principal stress corresponding to the inversion window, determine the first spatial step size and the second spatial step size of the inversion window, and take the larger value as the reference length. Using the geometric center of the spatial range of the inversion window as the midpoint, a center line segment with a length equal to the reference length is generated along the direction of the maximum principal stress. The coordinates of the two endpoints of the center line segment are obtained by moving the midpoint in the direction of the maximum principal stress and the opposite direction of the maximum principal stress by half the reference length, respectively. Multiply the reference length by (1 + principal stress rotational composite index) to obtain the partition width; The deformable partition is defined as: a three-dimensional tubular neighborhood consisting of all spatial points whose shortest radial distance from the center line segment does not exceed the width of the partition; where the shortest radial distance refers to the minimum value of the straight-line distances from the spatial points to any point on the center line segment; Events occurring within the inversion window time range are selected from the set of seismic events to form a subset of events within the window. For each seismic event, the shortest radial distance from the updated spatial location of the seismic event to the center line segment is calculated. Using the partition width as the scale parameter, an exponential decay function is used to calculate the weighting coefficients. The weighting coefficient is equal to the square of the ratio of the negative shortest radial distance of the natural exponential function to the partition width. If a seismic event falls into multiple deformable partitions, it is normalized and allocated according to the weight of the seismic event in each deformable partition based on the weighting coefficient of the corresponding deformable partition. That is, the weight of each deformable partition is equal to the corresponding weighting coefficient divided by the sum of the weighting coefficients of all deformable partitions it falls into.

[0014] Furthermore, within the deformable partition, a secondary stress field inversion is performed based on the weighted seismic events, outputting the corrected directions of the maximum principal stresses, including: For each deformable partition, read all seismic events within the deformable partition to form an event list; For each mine tremor event in the event list, extract the weighting coefficients based on the sum of squared residuals; And the weighting coefficients obtained based on the shortest straight-line distance; multiply the weighting coefficients and the weighting coefficients to obtain the comprehensive weighting coefficients; bind the comprehensive weighting coefficients to the fault plane normal vector and the slip direction vector of the seismic event to form a set of regional events; A partition objective function is constructed with the stress tensor as the unknown. The partition objective function is the sum of the angle mismatch terms of all mining seismic events in the deformable partition after being weighted by a comprehensive weighting coefficient. The partition objective function is minimized using a quasi-Newton optimization algorithm to obtain the stress tensor of the deformable partition. The convergence criterion is set as the relative difference between the objective function values ​​of two adjacent iterations being less than a preset minimum value. When the convergence criterion is met, the iteration stops and the current stress tensor is output. The current stress tensor is decomposed into three eigenvalues ​​and three corresponding eigenvectors. The eigenvector corresponding to the largest eigenvalue is selected as the corrected maximum principal stress direction of the deformable partition. The maximum principal stress directions of all deformable partitions are collected to form the corrected maximum principal stress direction field.

[0015] The beneficial effects of this invention include: by constructing a complete source analysis and inversion system, it organically combines mine seismic waveform feature extraction, event identification and location, source parameter inversion, initial stress field inversion, stress gradient calculation, principal stress rotation composite index evaluation, and deformable zoning with secondary inversion, achieving a full-process analysis of mine seismic events from their occurrence mechanism to stress field distribution. This invention not only significantly improves the accuracy of mine seismic event identification and location but also dynamically reflects the changing patterns of principal stress directions under complex geological conditions, effectively overcoming the shortcomings of traditional methods in terms of insufficient accuracy and poor adaptability. It provides a scientific and reliable basis for mine disaster early warning, rock strata stability analysis, and safe production decision-making. Attached Figure Description

[0016] Figure 1 This is a block diagram of the source analysis and inversion system based on mine seismic monitoring of the present invention. Detailed Implementation

[0017] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0018] like Figure 1 As shown, the source analysis and inversion system based on mine seismic monitoring includes: The waveform feature extraction module acquires continuous waveforms from mine seismic monitoring and generates continuous parameters characterizing waveform energy changes. The event identification and localization module generates mine seismic events based on the spatiotemporal aggregation patterns of continuous parameters, and determines the spatial location and occurrence time of the mine seismic events; The source parameter inversion module uses the initial motion parameters and amplitude relationships of seismic waves to inversely derive the source mechanism parameters of mining earthquake events. The initial stress field inversion module, based on the source mechanism parameters and the spatiotemporal distribution of mining earthquake events, inverts the direction of the maximum principal stress in the regional stress field; The stress gradient calculation module determines the first spatial direction by combining the mining trajectory of the working face and the second spatial direction by combining the stratigraphic layering information, and calculates the spatial variation rate of the maximum principal stress direction in the first and second spatial directions. The coupling index calculation module obtains the interlayer mechanical difference index based on the mechanical parameters and thickness information of each layer of the composite top plate, and determines the principal stress rotation composite index by combining the two spatial change rates and the interlayer mechanical difference index. The deformable partitioning module generates adaptive deformable partitions along the direction of maximum principal stress based on the principal stress rotation composite index, and allocates and weights mine seismic events according to their spatial position relationship with the partition centerline. The stress field fine inversion module performs secondary stress field inversion based on weighted seismic events within the deformable partition, and outputs the corrected direction of the maximum principal stress.

[0019] In one embodiment of the present invention, acquiring continuous waveforms from mine seismic monitoring and generating continuous parameters characterizing waveform energy changes includes: The three-component raw waveform sequences collected from each station are subjected to mean removal and bandpass filtering to obtain standardized waveform sequences; In detail, the original waveforms undergo standardization processing. The three-component original waveform sequences collected by each monitoring station must first undergo mean removal and bandpass filtering. Mean removal eliminates the constant DC component in the waveform to avoid baseline interference in energy calculations; bandpass filtering preserves the effective frequency range of the seismic signal, ultimately resulting in a stable standardized waveform sequence that ensures the comparability of waveforms from different stations and at different times.

[0020] The short-time average energy continuity, long-time average energy continuity, and energy ratio continuity are calculated based on standardized waveform sequences, as follows: The short-time average energy continuity is a time series obtained by averaging the absolute values ​​of a standardized waveform sequence over a first preset time window. In detail, a first preset time window is selected (used to capture instantaneous energy changes), and the absolute values ​​of the standardized waveform sequence are averaged over time within this time window to obtain a short-time average energy sequence that changes continuously with time. Its function is to quickly respond to sudden increases in waveform energy (such as energy surges during mine tremors) and sensitively capture short-term fluctuations in the signal.

[0021] The long-term average energy continuity is a time series obtained by averaging the absolute values ​​of a standardized waveform sequence over a second preset time window. In detail, a second preset time window is selected, and the absolute values ​​of the standardized waveform sequence are averaged over time within this time window to obtain a long-term average energy sequence. Its function is to characterize the long-term background energy state of the signal.

[0022] The length of the second preset time window is greater than the length of the first preset time window; The energy ratio continuous quantity is a time series obtained by comparing the short-time average energy continuous quantity with the long-time average energy continuous quantity at the same time point; In detail, the energy ratio sequence is obtained by comparing the short-term average continuous energy with the long-term average continuous energy at the same time point. Since short-term energy is sensitive to instantaneous changes and long-term energy reflects the background, this ratio can amplify the relative energy changes during a seismic event (e.g., when short-term energy increases sharply while long-term energy is stable, the ratio increases significantly), effectively highlighting the energy characteristics of potential seismic events.

[0023] The short-time average energy continuous quantity, long-time average energy continuous quantity, and energy ratio continuous quantity are aggregated according to the station number to form a continuous parameter set.

[0024] In one embodiment of the present invention, generating mine seismic events based on the spatiotemporal aggregation law of continuous parameters and determining the spatial location and occurrence time of the mine seismic events includes: The spatial coordinates of each station are obtained to form a station spatial location set, and the parameters containing the seismic wave propagation velocity of each rock layer are obtained to form a propagation velocity parameter set; In detail, the spatial coordinates of all monitoring stations are collected to form a set of station spatial locations, which is used to calculate the spatial distance between the seismic source and the station; at the same time, the seismic wave propagation velocity of each rock layer is collected to form a set of propagation velocity parameters.

[0025] Based on the set of station spatial locations and the set of propagation velocity parameters, the straight-line distance from the candidate source point to each station is calculated, and then the straight-line distance is divided by the seismic wave propagation velocity of the corresponding rock layer to obtain the theoretical propagation time. The spatiotemporal energy aggregation function is obtained by summing the continuous parameters of energy ratio of each station after time alignment according to the theoretical propagation time. In detail, for each candidate source point, the straight-line distance from it to each station is calculated based on the station spatial location set. This distance is then divided by the seismic wave propagation velocity of the corresponding rock layer to obtain the theoretical propagation time of the seismic wave from the candidate point to each station. Since the distances between different stations and the source vary, the wave arrival times differ. Therefore, the energy ratio continuity parameter of each station needs to be time-aligned according to the theoretical propagation time. This involves adjusting the time axis of the energy sequence of each station so that the energy signals theoretically generated by the same source coincide in time. The aligned energy ratio parameters are then summed to obtain the spatiotemporal energy aggregation function. The local maxima of the spatiotemporal energy aggregation function correspond to the regions of highest energy concentration, indicating the spatiotemporal location of the potential source.

[0026] Based on the initial values ​​of occurrence time and spatial location, the time corresponding to the maximum value in the energy ratio continuity of each station is extracted as the arrival time observation within a fixed window centered on the theoretical propagation time; In detail, the time and spatial coordinates of the local maximum point in the spatiotemporal energy aggregation function are used as the initial values ​​of the occurrence time and spatial location of the seismic event. Based on this, within a fixed window centered on the theoretical propagation time (the window size is set according to the time error range of wave propagation) in the continuous energy ratio of each station, the time corresponding to the energy maximum is extracted as the actual arrival time of the seismic wave recorded by that station (i.e., arrival time observation).

[0027] Construct the arrival residual for each station; where the arrival residual is the sum of the arrival observation time and the theoretical propagation time; by minimizing the sum of squared arrival residuals for all stations, the occurrence time and spatial location of the seismic event are obtained. In detail, the arrival residual for each station is constructed: this is the arrival observation minus the sum of the initial occurrence time and the theoretical propagation time. This residual reflects the deviation between the theoretical prediction and the actual observation. By minimizing the sum of squared arrival residuals for all stations, the occurrence time and spatial location that minimizes the deviation are determined. This process is based on the least squares principle, iteratively adjusting the time and location parameters to achieve the highest overall agreement between the theoretical propagation time and the arrival observation, thus obtaining a more accurate event occurrence time and spatial location.

[0028] Select any two earthquake events from the set of earthquake events whose time interval is less than a first set value or whose spatial distance is less than a second set value to form an earthquake event pair; calculate the difference between the arrival observation difference and the theoretical propagation time difference for each station for the earthquake event pair, as the double-difference residual; by minimizing the sum of squares of the double-difference residuals of all stations, solve for the spatial position correction of each earthquake event; add the corresponding spatial position correction to the spatial position to obtain the updated spatial position of the earthquake event.

[0029] In detail, due to the complex geological conditions, the location of a single event may be affected by velocity model errors, thus requiring further optimization. From the set of seismic events, events with time intervals less than a first preset value (e.g., events occurring consecutively within a short period) or spatial distances less than a second preset value (e.g., spatially adjacent events) are selected to form seismic event pairs. The difference between the observed arrival time difference (the difference in arrival times of the two events recorded by the two stations) and the theoretical propagation time difference (the difference in theoretical propagation times from the two events to the station) for each station is calculated as the double-difference residual. The double-difference residual can effectively eliminate the combined influence of velocity model errors (similar to the propagation path errors of the same station for the two events, which can be canceled out in the difference). By minimizing the sum of squared double-difference residuals for all stations, the spatial location correction for each event is obtained. This correction is then added to the spatial location obtained in step four to finally obtain the updated spatial location.

[0030] In one embodiment of the present invention, the focal mechanism parameters of a mining seismic event are derived by inversely using the initial motion parameters and amplitude relationship of seismic waves, including: Based on the occurrence time and updated spatial location of the mine earthquake event, the theoretical propagation time of the seismic waves to each station is calculated; In detail, based on the occurrence time and updated spatial location of the seismic event, and combined with the propagation velocity parameter set, the theoretical propagation time of seismic waves to each station is recalculated. This ensures that the subsequent wave window selection can accurately cover the effective signal segments of both P-waves and S-waves.

[0031] In the standardized waveform sequence, a fixed-length longitudinal wave window and a transverse wave window are cut out with the theoretical propagation time as the center; In detail, within the standardized waveform sequence, fixed-length longitudinal wave windows and transverse wave windows are extracted, centered on the theoretical propagation time. The longitudinal wave, with its faster propagation speed, arrives first, while the transverse wave, with its slower speed, arrives later. Separate window settings avoid signal interference between the two waves, ensuring accurate extraction of their respective characteristic parameters.

[0032] The positive and negative polarities of the first wave of the P-wave are extracted from the P-wave window as the initial motion polarity of the P-wave, and the peak amplitudes are extracted from the P-wave window and the S-wave window as the P-wave amplitude and the S-wave amplitude, respectively. In detail, the positive and negative polarities of the initial P-wave are extracted from the P-wave window as the initial polarity of the P-wave. The initial vibration direction of the particles at the source is directly related to the stress state of the source. The peak amplitudes are extracted from the P-wave and S-wave windows respectively as the P-wave amplitude and S-wave amplitude. The magnitude of the amplitude reflects the energy intensity of the wave when it propagates to the station, but it needs further correction due to the influence of the propagation path.

[0033] The propagation distance is calculated based on the updated spatial location of the seismic event and the spatial location of the station. The medium attenuation coefficient is calculated based on the quality factor, seismic wave frequency and propagation velocity in the propagation velocity parameter set. Among them, the medium attenuation coefficient is inversely proportional to the quality factor, directly proportional to the frequency, and inversely proportional to the propagation speed; Geometric diffusion correction is performed on the P-wave and S-wave amplitudes using the propagation distance. The correction method is to divide the P-wave amplitude or S-wave amplitude by the propagation distance. The amplitude of the longitudinal wave or the transverse wave is attenuated by using the medium attenuation coefficient and the propagation distance. The correction method is to multiply the amplitude of the longitudinal wave or the transverse wave by the exponential attenuation term of the product of the medium attenuation coefficient and the propagation distance. In detail, the propagation distance is calculated based on the updated spatial location of the seismic event and the spatial location of the seismic station. The medium attenuation coefficient is calculated based on the quality factor of the seismic wave frequency and propagation velocity from the propagation velocity parameter set. The medium attenuation coefficient is inversely proportional to the quality factor (the lower the quality factor, the stronger the medium absorption), directly proportional to the frequency (high-frequency waves attenuate faster), and inversely proportional to the propagation velocity (low-speed waves interact more frequently in the medium). Geometric diffusion correction is achieved by dividing the P-wave amplitude or S-wave amplitude by the propagation distance, eliminating the natural energy attenuation caused by spherical diffusion during wave propagation. Attenuation correction is achieved by multiplying the P-wave amplitude or S-wave amplitude by the exponential attenuation term of the product of the medium attenuation coefficient and the propagation distance, eliminating energy loss caused by medium absorption. After the secondary correction, the amplitude is closer to the original radiation characteristics at the source.

[0034] The amplitude ratio observation is obtained by comparing the second-corrected shear wave amplitude with the corrected longitudinal wave amplitude. In detail, the amplitude ratio observation is obtained by comparing the second-corrected shear wave amplitude with the corrected p-wave amplitude. This ratio is not affected by the absolute value of the total source energy, but is only related to the source radiation mode, which is determined by the source mechanism parameters.

[0035] An objective function for source mechanism inversion is constructed, which includes: a P-wave first motion polarity consistency term and an amplitude ratio residual term. The P-wave first motion polarity consistency term is the sum of the number of inconsistencies between the observed P-wave first motion polarity at each station and the P-wave first motion polarity predicted by the candidate source mechanism; the amplitude ratio residual term is the sum of the squares of the logarithmic differences between the observed amplitude ratio at each station and the amplitude ratio predicted by the candidate source mechanism. In detail, the objective function comprises two terms: a P-wave first motion polarity consistency term and an amplitude ratio residual term. The P-wave first motion polarity consistency term is the sum of the number of discrepancies between the observed P-wave first motion polarity at each station and the P-wave first motion polarity predicted by the candidate source mechanism, used to measure the degree of polarity matching. The amplitude ratio residual term is the sum of the squared logarithmic differences between the observed amplitude ratio at each station and the amplitude ratio predicted by the candidate source mechanism, used to measure the degree of agreement in amplitude ratio. The combination of these two terms provides a comprehensive assessment of the overall matching degree between the candidate mechanism and the observed data.

[0036] The objective function for source mechanism inversion is minimized using a grid search algorithm. The parameters of the grid search are: strike angle from 0 to 360 degrees, dip angle from 0 to 90 degrees, and slip angle from -90 degrees to +90 degrees, with a preset step size. The strike angle, dip angle, and slip angle of the seismic event are obtained by solving the algorithm. The minimum value of the objective function is used as the sum of squared residuals as a quality index of the source mechanism parameters. The strike angle, dip angle, slip angle, and sum of squared residuals are aggregated to form a set of source mechanism parameters.

[0037] In detail, the parameter ranges for the grid search are set as follows: strike angle 0 to 360 degrees (the direction of the fault's projection onto the horizontal plane), dip angle 0 to 90 degrees (the angle between the fault and the horizontal plane), and slip angle -90 degrees to +90 degrees (the angle between the relative slip directions of the two fault blocks and the strike). The step size is a preset angle (balancing computational efficiency and accuracy). By traversing all candidate parameter combinations, the corresponding objective function value is calculated, and the parameter combination that minimizes the objective function is selected as the optimal solution, yielding the strike, dip, and slip angles of the seismic event. The minimum value of the objective function is used as the residual sum of squares, reflecting the reliability of the inversion results (the smaller the value, the higher the matching degree). Finally, the strike, dip, slip angles, and residual sums of squares are aggregated to form the focal mechanism parameter set.

[0038] It should be noted that the quality factor represents the energy attenuation characteristics of seismic waves propagating in a medium, reflecting the strength of the medium's absorption of seismic wave energy. Obtaining the quality factor involves: conducting acoustic or seismic wave propagation experiments on rock samples of typical lithologies in the mining area through laboratory measurements, monitoring the energy attenuation pattern of the waves propagating in the samples, and then calculating the corresponding lithology's quality factor.

[0039] In one embodiment of the present invention, based on the source mechanism parameters and the spatiotemporal distribution of mining seismic events, the direction of the maximum principal stress in the regional stress field is obtained by inversion, including: Multiple inversion windows are defined; the temporal and spatial ranges of the inversion windows are set based on the density of seismic events, with the temporal range not exceeding ten times the sampling period of the stations and the spatial range not less than the average spacing of the station distribution. In detail, the inversion window is an analytical unit designed to aggregate seismic events within a specific spatiotemporal range. Its temporal and spatial ranges are determined based on the density of seismic events: the temporal range does not exceed ten times the station sampling period to ensure the temporal correlation of events within the window (avoiding excessive stress field changes due to excessively long time spans); the spatial range is not less than the average spacing of the stations to ensure sufficient station coverage within the window and guarantee spatial representativeness. A reasonable window size is crucial for balancing inversion stability (data volume) and resolution (spatiotemporal detail).

[0040] From the set of mining tremor events, select events whose occurrence time is within the time range and whose updated spatial location is within the aforementioned spatial range to form an inversion event subset; In detail, events that occurred within the inversion window time range and whose updated spatial location was within the window space range were selected from the set of mining seismic events to form a subset of inversion events. These events collectively reflect the stress state within the spatiotemporal region corresponding to the window.

[0041] For each seismic event in the inverted event subset, the fault plane normal vector and the slip direction vector are derived through three-dimensional geometric relationships based on the corresponding strike angle, dip angle and slip angle; wherein, the fault plane normal vector is perpendicular to the fault plane, and the slip direction vector points along the fault plane in the slip direction. In detail, for each seismic event in the inverted event subset, based on its focal mechanism parameters (strike angle, dip angle, slip angle), two key vectors are derived through three-dimensional geometric relationships: the fault plane normal vector (perpendicular to the fault plane, pointing towards the hanging wall or footwall) and the slip direction vector (along the tangent direction of the fault plane, pointing towards the direction of relative slip between the two sides of the fault). These two vectors quantitatively describe the spatial attitude and motion of the fault, serving as a bridge connecting the focal mechanism and the stress field. Theoretically, the fault slip direction should be consistent with the direction of shear stress on the fault plane.

[0042] In detail, the strike angle is the clockwise angle between the intersection of the fault plane and the horizontal plane (the strike line) and its deviation from true north, ranging from 0 degrees to 360 degrees. For example, a strike angle of 90 degrees indicates that the strike line is in the direction of true east.

[0043] In detail, the dip angle is the angle between the fault plane and the horizontal plane, measured downwards from the strike line, ranging from 0 degrees to 90 degrees. A dip angle of 90 degrees indicates that the fault plane is vertical, and 0 degrees indicates that the fault plane is horizontal.

[0044] In detail, the slip angle is the angle between the relative slip directions of the two sides of a fault and the strike line. It is measured along the fault plane: the slip direction is positive (right-handed) if it is clockwise and negative (left-handed) if it is counterclockwise, with a range of -90 degrees to 90 degrees.

[0045] The normal vector must be perpendicular to the fault plane, and its direction is determined by both the strike angle and dip angle, as derived below: Determine the dip direction of the fault plane: The dip direction is the direction of rotation of the strike line 90 degrees clockwise (i.e., the direction in which the fault plane is tilted). If the strike angle is φ, then the horizontal angle of the dip direction is φ + 90 degrees (subtract 360 degrees if it exceeds 360 degrees).

[0046] Calculate the horizontal component of the normal vector: The projection of the normal vector onto the horizontal plane is perpendicular to the dip direction (since the normal is perpendicular to the fault plane), so the angle of the horizontal component is φ + 180 degrees (opposite to the dip). The magnitude of the horizontal component is determined by the dip angle. If the normal vector is a unit vector, then the magnitude of the horizontal component is sin(dip angle) (the larger the dip angle, the steeper the fault plane, and the higher the proportion of the horizontal component).

[0047] Calculate the perpendicular component of the normal vector: the magnitude of the perpendicular component (in the z-axis direction) is cos(tilt angle), and its direction is upward (following the right-hand rule: four fingers along the direction, palm facing the inclination, and thumb in the normal direction).

[0048] Composite normal vector components: x component (due north): sin(tilt angle) × cos(φ + 180 degrees) = -sin(tilt angle) × cosφ; y-component (due east direction): sin(angle of inclination) × sin(φ + 180 degrees) = -sin(angle of inclination) × sinφ; z-component (vertical direction): cos(tilt angle); That is, the normal vector is [-sinδ・cosφ,-sinδ・sinφ,cosδ]; where δ is the dip angle and φ is the strike angle; The slip direction vector must be along the tangent direction of the fault plane and is determined by the strike angle, dip angle, and slip angle. The derivation process is as follows: Determine the two reference directions within the fault plane: Direction of travel: along the horizontal direction of the travel line, with vector [cosφ,sinφ,0].

[0049] Dip direction: the direction of inclination along the fault plane (perpendicular to the strike direction), and the vector is [-sinφ・cosδ,cosφ・cosδ,-sinδ] (composed of the horizontal dip direction and the vertical downward component).

[0050] The sliding angle is converted into components of two reference directions: the sliding angle ψ represents the angle between the sliding direction and the heading direction, therefore the sliding direction vector can be synthesized by decomposing the heading direction vector and the dipping direction vector by angle. The component along the direction of travel: cosψ (when the sliding angle is 0, the sliding direction is consistent with the direction of travel).

[0051] The component along the direction of inclination: sinψ (when the sliding angle is 90°, the sliding direction is downward along the direction of inclination).

[0052] Synthesize the sliding direction vector: x component: cosφ・cosψ+(-sinφ・cosδ)・sinψ; y component: sinφ・cosψ+(cosφ・cosδ)・sinψ; z-component: 0・cosψ+(-sinδ)・sinψ; That is, the sliding direction vector d = [cosφ・cosψ-sinφ・cosδ・sinψ,sinφ・cosψ+cosφ・cosδ・sinψ,-sinδ・sinψ]; where ψ is the sliding angle.

[0053] Based on the residual sum of squares of a mining earthquake event, weighting coefficients are generated through an exponentially monotonically decreasing function. The exponentially monotonically decreasing function is the weighting coefficients being the product of the negative residual sum of squares of the natural exponential function and a preset scale parameter raised to the power of the product. In detail, based on the residual sum of squares of seismic events (the minimum value of the objective function in focal mechanism inversion, reflecting the reliability of focal parameters), weighting coefficients are generated through an exponentially monotonically decreasing function. Specifically, the weighting coefficients are equal to the power of the product of the negative residual sum of squares of the natural exponential function and the preset scale parameter. The smaller the residual sum of squares, the more reliable the focal mechanism parameters, and the larger the corresponding weighting coefficients, allowing reliable events to play a greater role in the inversion and reducing the interference of unreliable data.

[0054] An objective function is constructed using the stress tensor as the unknown. The objective function is the sum of the angle mismatch terms for each seismic event, weighted by a weighting coefficient. The angle mismatch term is the square of the angle between the sliding direction vector and the shear stress vector of the event. The shear stress vector is the remaining part after subtracting the normal component of the normal vector of the fault plane from the stress tensor. The objective function is minimized using the least squares optimization algorithm, and the stress tensor is obtained by solving it. In detail, the stress tensor is a second-order tensor describing the stress state at a point in space, containing six independent components. The objective function is defined as the weighted sum of the angle mismatch terms for each seismic event: the angle mismatch term is the square of the angle between the slip direction vector and the shear stress vector of the event (the smaller the angle, the more consistent the theoretical shear stress direction is with the actual slip direction); the shear stress vector is obtained by subtracting the normal component of the stress tensor from the normal vector of the fault plane (i.e., the shear component of the stress on the fault plane, which is the direct force driving fault slip). The objective function constructed with the stress tensor as the unknown is to minimize the overall deviation between the theoretical shear stress direction and the actual slip direction for all seismic events.

[0055] In detail, the components of the stress tensor are iteratively adjusted using a least-squares optimization algorithm to minimize the objective function constructed with the stress tensor as the unknown, thus obtaining the stress tensor that best matches the subset of inverted events. Based on the mechanical principle that the fault slip direction is consistent with the shear stress direction, the stress state is inferred through data fitting.

[0056] The stress tensor is decomposed into three eigenvalues ​​and three corresponding eigenvectors. The eigenvector corresponding to the largest eigenvalue is selected as the direction of the maximum principal stress in the inversion window. The direction of the maximum principal stress is bound to the midpoint of the time range and the geometric center of the spatial range of the inversion window. The maximum principal stress direction field is formed by traversing all inversion windows.

[0057] In detail, the stress tensor is decomposed into three eigenvalues ​​(reflecting the magnitude of the principal stresses) and three corresponding eigenvectors (reflecting the directions of the principal stresses). The eigenvectors are mutually perpendicular and correspond to the three principal stress directions, respectively. The eigenvector corresponding to the largest eigenvalue is the direction of the maximum principal stress in the inversion window (the direction of the strongest stress in the region).

[0058] In one embodiment of the present invention, a first spatial direction is determined by combining the mining trajectory of the working face, a second spatial direction is determined by combining the stratigraphic stratification information, and the spatial variation rate of the direction of maximum principal stress in the first and second spatial directions is calculated, including: For the geometric center of the spatial range of the inversion window, the nearest target trajectory point is searched in the mining trajectory set of the working face by calculating the Euclidean distance; The tangential vector is calculated by the coordinate difference between the target trajectory point and the preceding and following trajectory points adjacent to the target trajectory point. The first spatial direction is obtained by dividing the tangential vector by the magnitude. The first spatial direction is the tangential direction of the working face advancement. In detail, the first spatial direction corresponds to the tangential direction of the working face advancement and needs to be calculated in conjunction with the working face mining trajectory: For the geometric center of the spatial range of each inversion window (representing the spatial position of the inversion window), the nearest target trajectory point (i.e., the point on the mining trajectory closest to the window center) is searched in the working face mining trajectory set by calculating the Euclidean distance (the straight-line distance between two points); using the coordinate difference between this target trajectory point and its adjacent previous and next trajectory points, the tangential vector (describing the extension direction of the trajectory at that point) is calculated; the tangential vector is normalized by dividing it by its magnitude (the length of the vector) to obtain the first spatial direction in unit vector form. The first spatial direction reflects the instantaneous direction of the working face advancement and is the benchmark for analyzing the changes in mining stress along the advancement direction.

[0059] For the geometric center of the spatial range of the inversion window, the corresponding local bedding plane is extracted from the stratigraphic layer information set. The local bedding plane is described by the coordinates of multiple discrete points. The least squares fitting method is used to fit the coordinates of discrete points to the plane, and the plane equation of the local bedding plane is obtained. The normal vector is determined based on the coefficients of the plane equation. The second spatial direction is obtained by dividing the normal vector by the magnitude. The second spatial direction is perpendicular to the local stratification plane. In detail, the second spatial direction, perpendicular to the local bedding plane, needs to be calculated in conjunction with stratigraphic information: For the geometric center of the inversion window's spatial range, the corresponding local bedding plane (a rock stratum described by multiple discrete point coordinates) is extracted from the stratigraphic information set; the least squares fitting method is used to fit the coordinates of these discrete points to a plane. By minimizing the sum of the squared distances from each point to the fitted plane, the plane equation of the local bedding plane is obtained (in the form ax + by + cz + d = 0); based on the coefficients (a, b, c) of the plane equation, the normal vector (a vector perpendicular to the plane) is determined, and the normal vector is normalized by dividing it by its magnitude to obtain the second spatial direction in unit vector form. The second spatial direction, perpendicular to the rock stratum, is used to analyze the variation of stress direction along the direction perpendicular to the bedding, adapting to the layered structural characteristics of the composite roof.

[0060] The angle between the projection of the direction of the maximum principal stress onto the horizontal plane and the north direction is taken as the direction angle, which is defined by a 180-degree period. Angle dewinding is applied to the direction angle sequence to eliminate jumps when crossing zero or 180 degrees. In detail, the angle between the projection of the maximum principal stress direction onto the horizontal plane and the north direction is defined as the orientation angle. The orientation angle is defined in a 180-degree cycle (0 to 180 degrees) (because the stress direction has bidirectional equivalence, such as 0 degrees and 180 degrees representing the same direction). Since the orientation angle will jump when crossing 0 degrees or 180 degrees (e.g., from 170 degrees to 10 degrees, it should actually be a continuous change of -160 degrees), the orientation angle sequence needs to be de-wrapped. By adding or subtracting integer multiples of 180 degrees, the jump is corrected into a continuous angle change, ensuring that the orientation angle sequence is numerically continuous.

[0061] Read the step size of the inversion window in the first spatial direction and the step size in the second spatial direction. The step size is equal to the Euclidean distance between the geometric centers of the spatial extents of adjacent inversion windows in the corresponding directions. In detail, the step size of the inversion window is read in the first and second spatial directions: the step size is equal to the Euclidean distance between the geometric centers of the spatial extents of adjacent inversion windows in the corresponding directions, that is, the straight-line distance between the centers of two adjacent windows along the first (or second) spatial direction. The step size reflects the spatial sampling interval and is the basic scale parameter for calculating the rate of change.

[0062] The rate of change of the orientation angle sequence is calculated using the central difference method; where the rate of change is the rate of change of the inversion window, which is equal to the difference in orientation angle between the forward adjacent inversion window and the reverse adjacent inversion window divided by twice the step size. For the first and last inversion windows in all inversion windows, forward difference or backward difference is used for calculation; that is, the rate of change of the first inversion window is the difference in the direction angle between the second inversion window and the first inversion window divided by the step size, and the rate of change of the last inversion window is the difference in the direction angle between the last inversion window and the second to last inversion window divided by the step size; thus forming the first spatial rate of change sequence and the second spatial rate of change sequence.

[0063] In detail, the rate of change of the orientation angle sequence in two spatial directions is calculated using the finite difference method: For the window in the middle of all inversion windows, the central difference method is used, that is, the rate of change of the window is equal to the difference in orientation angle between its forward adjacent window and its reverse adjacent window divided by twice the step size (using the average of the data on both sides to improve accuracy); for the first and last inversion windows, since there is a lack of adjacent data on one side, forward difference or backward difference is used respectively. The rate of change of the first window is the difference in orientation angle between the second window and the first window divided by the step size, and the rate of change of the last window is the difference in orientation angle between the last window and the second-to-last window divided by the step size. Through the above calculations, the first spatial rate of change sequence (the rate of change of stress direction along the working face advancement direction) and the second spatial rate of change sequence (the rate of change of stress direction along the direction perpendicular to the bedding) are finally formed.

[0064] It should be noted that the stratigraphic layer information set is a dataset used to describe the layered structure of underground rock strata in a mine. It contains the spatial distribution characteristics of each rock stratum, and may specifically include information such as the vertical range of different rock strata, the coordinates of discrete points of the bedding planes, the rock stratum type and physical properties.

[0065] It should be noted that local bedding planes refer to the bedding planes within a small-scale region corresponding to the geometric center of the inversion window's spatial range. They are stratified structural interfaces formed by differences in material composition, grain size, etc., during the deposition or formation of rock strata. Their spatial morphology is characterized by multiple discrete point coordinates (such as the three-dimensional coordinates of sampling points on the bedding plane). These discrete points reflect the undulation or tilt characteristics of the bedding planes within this local region.

[0066] It should be noted that the mathematical equation of a plane in space can be expressed as ax + by + cz + d = 0, where a, b, c, and d are constants, and a, b, and c are not all zero simultaneously. The geometric meaning of this equation is that any point (x, y, z) on the plane satisfies this equation. For a plane, the normal vector is the vector perpendicular to the plane, and its direction can be directly determined by the coefficients of the plane's equation. The coefficients a, b, and c correspond to the components of the normal vector along the x, y, and z coordinate axes, respectively; that is, the normal vector is (a, b, c).

[0067] In one embodiment of the present invention, an interlayer mechanical difference index is obtained based on the mechanical parameters and thickness information of each layer of the composite top plate, and the principal stress rotation composite index is determined by combining two spatial change rates and the interlayer mechanical difference index, including: Within the vertical range corresponding to the geometric center of the spatial range of each inversion window, read the thickness and equivalent Young's modulus of each rock layer within the vertical range; In detail, for the geometric center of each inversion window spatial range, determine its corresponding vertical range (i.e., the rock strata within a certain depth above and below the center), and read the thickness (the size of the rock strata in the vertical direction) and equivalent Young's modulus (a parameter characterizing the rock strata's ability to resist elastic deformation; the larger the value, the higher the stiffness) of each rock strata within this vertical range.

[0068] The proportional weight is obtained by calculating the proportion of the thickness of each rock layer to the total thickness of all rock layers. In detail, for all rock strata within the vertical range, the proportion of each rock stratum's thickness to the total thickness (the sum of the thicknesses of all rock strata) is calculated to obtain the proportional weight of that rock stratum. The proportional weight reflects the relative proportion of each rock stratum within the vertical range, with thicker strata having a greater weight to ensure that their mechanical properties occupy a more important position in the overall analysis.

[0069] The weighted average of the equivalent Young's modulus is calculated based on the proportional weights; where the weighted average is equal to the sum of the products of the equivalent Young's modulus of each rock layer and the corresponding proportional weight. In detail, the equivalent Young's modulus of each rock layer is multiplied by its corresponding proportional weight, and then all products are summed to obtain the weighted average of the equivalent Young's modulus. This value comprehensively reflects the overall average stiffness level of the composite roof and serves as a benchmark for measuring interlayer differences.

[0070] Calculate the weighted standard deviation of the equivalent Young's modulus; where the weighted standard deviation is equal to the square root of the sum of the squares of the differences between the equivalent Young's modulus and the weighted average of each rock layer and the corresponding proportional weights. In detail, the difference between the equivalent Young's modulus of each rock layer and the weighted average is calculated first. The square of the difference is then multiplied by the proportional weight of that rock layer. All the results are then summed, and finally the square root is taken to obtain the weighted standard deviation. The weighted standard deviation quantifies the degree of dispersion of the equivalent Young's modulus of each rock layer relative to the average value. The larger the value, the more significant the difference in stiffness between layers.

[0071] The ratio of the weighted standard deviation to the weighted mean is used as the interlayer mechanical difference index, which is used to characterize the dispersion of the interlayer mechanical properties of the composite roof. In detail, the ratio of the weighted standard deviation to the weighted mean is defined as the interlayer mechanical difference index. The interlayer mechanical difference index eliminates the influence of absolute values ​​and more intuitively characterizes the degree of inhomogeneity of the interlayer mechanical properties of the composite roof. A larger interlayer mechanical difference index indicates a more significant difference in stiffness between different rock strata and a more dispersed interlayer mechanical property.

[0072] For each inversion window, the first spatial rate of change is multiplied by the first spatial step size to obtain the standard first spatial change; the second spatial rate of change is multiplied by the second spatial step size to obtain the standard second spatial change. In detail, for each inversion window, the first spatial rate of change (the rate of change of the principal stress direction along the working face advancement direction) is multiplied by the first spatial step size (the distance between the centers of adjacent windows along this direction) to obtain the standard first spatial rate of change, which reflects the actual change in the principal stress direction within a unit step size along the working face advancement direction; similarly, the second spatial rate of change (the rate of change of the principal stress direction along the direction perpendicular to the bedding) is multiplied by the second spatial step size to obtain the standard second spatial rate of change, which reflects the change in the stress direction along the direction perpendicular to the bedding.

[0073] The sum of squares of the standard first space variation and the standard second space variation is calculated, and then the square root of the result is taken to obtain the comprehensive rotational intensity. In detail, the sum of squares of the standard first spatial variation and the standard second spatial variation is calculated, and then the square root of the result is taken to obtain the comprehensive rotational intensity. This value is based on the principle of vector composition, combining the stress rotation amplitudes in two key directions to quantify the overall rotational intensity of the principal stress directions in space.

[0074] Multiplying the comprehensive rotational intensity with the interlayer mechanical difference index yields the principal stress rotational composite index; traversing all inversion windows forms a sequence of principal stress rotational composite indices.

[0075] In detail, the principal stress rotation composite index is obtained by multiplying the comprehensive rotation intensity with the interlayer mechanical difference index. This index couples the rotation intensity along the principal stress direction with the interlayer mechanical difference of the composite roof, effectively characterizing their synergistic effect. A larger principal stress rotation composite index indicates more intense rotation along the principal stress direction in rock strata with significant differences in mechanical properties, leading to higher potential stress concentration and seismic risks. After traversing all inversion windows, a principal stress rotation composite index sequence is formed, fully presenting the spatial distribution characteristics of this coupling effect.

[0076] In one embodiment of the present invention, based on the principal stress rotational composite index, an adaptive deformable partition is generated along the direction of the maximum principal stress, and mine seismic events are allocated and weighted according to their spatial positional relationship with the partition centerline, including: For each inversion window, read the direction of the maximum principal stress corresponding to the inversion window, determine the first spatial step size and the second spatial step size of the inversion window, and take the larger value as the reference length. In detail, for each inversion window, the direction of the maximum principal stress (as the reference direction for partition extension) is read, and the first spatial step size (sampling interval along the working face advancement direction) and the second spatial step size (sampling interval along the direction perpendicular to the bedding plane) of the inversion window are extracted. The larger of the two values ​​is taken as the reference length. The selection of the reference length ensures that the partition has sufficient spatial coverage and balances the differences in sampling density in different directions.

[0077] Using the geometric center of the spatial range of the inversion window as the midpoint, a center line segment with a length equal to the reference length is generated along the direction of the maximum principal stress. The coordinates of the two endpoints of the center line segment are obtained by moving the midpoint in the direction of the maximum principal stress and the opposite direction of the maximum principal stress by half the reference length, respectively. In detail, a center line segment with a length equal to the reference length is generated along the direction of the maximum principal stress, with the geometric center of the inversion window spatial range as the midpoint. The coordinates of the two endpoints of the line segment are obtained by moving the midpoint along the direction of the maximum principal stress and the opposite direction by half the reference length, respectively. The center line segment serves as the central axis of the partition, and its direction is consistent with the direction of the maximum principal stress, ensuring that the partition extends along the dominant stress direction and adapts to the spatial distribution characteristics of the stress field.

[0078] Multiply the reference length by (1 + principal stress rotational composite index) to obtain the partition width; In detail, the width of the partition is determined by the reference length and the principal stress rotation composite index, specifically by multiplying the reference length by (1 plus the principal stress rotation composite index). The larger the principal stress rotation composite index (indicating intense stress rotation and significant interlayer mechanical differences), the wider the partition. This is because the stress field changes in such areas are complex, requiring a larger spatial range to gather relevant seismic events and avoid missing key events due to an overly narrow partition.

[0079] The deformable partition is defined as: a three-dimensional tubular neighborhood consisting of all spatial points whose shortest radial distance from the center line segment does not exceed the width of the partition; where the shortest radial distance refers to the minimum value of the straight-line distances from the spatial points to any point on the center line segment; In detail, the deformable partition is defined as a three-dimensional tubular neighborhood, containing all spatial points whose shortest radial distance from the centerline segment does not exceed the width of the partition. The shortest radial distance refers to the minimum straight-line distance from a spatial point to any point on the centerline segment, i.e., the perpendicular distance from the point to the segment (if the perpendicular falls on the segment) or the distance from the point to the endpoint of the segment (if the perpendicular extends beyond the segment's boundaries). The tubular neighborhood design allows the partition to extend along the centerline segment, and its width is dynamically adjustable, flexibly adapting to spatial variations in the stress field.

[0080] Events whose occurrence time falls within the time range of the inversion window are selected from the set of mining tremor events to form a subset of events within the window; In detail, events that occurred within the inversion window time range are selected from the set of seismic events, forming a subset of events within the window. These events are matched with the time range of the current inversion window, ensuring that the events used for zonal analysis are correlated with the corresponding spatiotemporal stress field characteristics.

[0081] For each seismic event, calculate the shortest radial distance from the updated spatial location of the seismic event to the centerline segment; In detail, for each seismic event in the event subset within the window, the shortest radial distance from its updated spatial location to the centerline segment is calculated. This distance quantifies the spatial proximity of the event to the partition's central axis; the smaller the distance, the closer the event is to the stress field characteristics represented by that partition.

[0082] Using the partition width as the scale parameter, the weighting coefficients are calculated using an exponential decay function; where the weighting coefficients are equal to the square of the ratio of the negative shortest radial distance of the natural exponential function to the partition width. In detail, using the partition width as the scale parameter, an exponential decay function is employed to calculate the weighting coefficient for each event. Specifically, the weighting coefficient is equal to the square of the negative of the natural exponential function (the ratio of the shortest radial distance to the partition width). This function's characteristic ensures that events closer to the center line segment have a larger weighting coefficient (approaching 1), while the coefficient decays exponentially (approaching 0) as the distance increases. This achieves spatial attenuation of the event's contribution to the partition, highlighting the dominant role of events near the center.

[0083] If a seismic event falls into multiple deformable zones simultaneously, the weight of the seismic event in each deformable zone is normalized and allocated according to the weighting coefficient of the corresponding deformable zone. That is, the weight of each deformable zone is equal to the corresponding weighting coefficient divided by the sum of the weighting coefficients of all the deformable zones it falls into.

[0084] In detail, if the spatial location of a seismic event falls into multiple deformable partitions simultaneously (since partitions may overlap), its weights need to be normalized: the weight assigned to each partition is equal to the weighting coefficient of the event in that partition divided by the sum of the weighting coefficients of all partitions it falls into. This allocation method ensures that the total weight of a single event is 1, avoids double counting, and ensures that the event's contribution in different partitions matches its spatial correlation.

[0085] In one embodiment of the present invention, a secondary stress field inversion is performed within a deformable partition based on weighted seismic events, and the corrected direction of the maximum principal stress is output, including: For each deformable partition, read all seismic events within the deformable partition to form an event list; In detail, for each deformable zone, all seismic events contained within that zone are read and compiled into an event list. These events were selected through a prior allocation and weighting process and are closely related to the stress field characteristics represented by the deformable zone.

[0086] For each mine tremor event in the event list, extract the weighting coefficients based on the sum of squared residuals; And the weighting coefficients obtained based on the shortest straight-line distance; multiply the weighting coefficients and the weighting coefficients to obtain the comprehensive weighting coefficients; bind the comprehensive weighting coefficients to the fault plane normal vector and the slip direction vector of the seismic event to form a set of regional events; In detail, for each seismic event in the event list, two types of weighting coefficients are extracted: one is a weighting coefficient based on the sum of squared residuals of the source mechanism parameters (reflecting the reliability of the source mechanism parameters; the smaller the residual, the larger the coefficient); the other is a weighted coefficient based on the shortest radial distance from the event to the center line of the partition (reflecting the spatial correlation between the event and the partition; the closer the distance, the larger the coefficient). These two types of coefficients are multiplied to obtain a comprehensive weighting coefficient. This comprehensive weighting coefficient integrates the event's reliability and spatial correlation, allowing more reliable events that are spatially closer to the partition core to play a greater role in the inversion, improving the stability and relevance of the inversion results. The comprehensive weighting coefficient is then bound to the fault plane normal vector and slip direction vector of the event to form a set of partitioned events containing mechanical characteristics and weighting information.

[0087] A partition objective function is constructed with the stress tensor as the unknown. The partition objective function is the sum of the angle mismatch terms of all mining seismic events in the deformable partition after being weighted by a comprehensive weighting coefficient. Using the stress tensor as the unknown, a zoning objective function is constructed: the zoning objective function is the sum of the angle mismatch terms of all seismic events within the deformable zone, weighted by a comprehensive weighting coefficient. The angle mismatch term is the square of the angle between the event's slip direction vector and the shear stress vector (the portion remaining after subtracting the normal component from the stress tensor's application to the fault plane's normal vector), used to measure the deviation between the theoretical shear stress direction and the actual fault slip direction. Through comprehensive weighting, events with high contribution (reliable and close to the zone's core) dominate the objective function.

[0088] The partition objective function is minimized using a quasi-Newton optimization algorithm to obtain the stress tensor of the deformable partition. The convergence criterion is set as the relative difference between the objective function values ​​of two adjacent iterations being less than a preset minimum value. When the convergence criterion is met, the iteration stops and the current stress tensor is output. In detail, a quasi-Newton optimization algorithm is used to iteratively adjust the components of the stress tensor to minimize the partitioned objective function. The quasi-Newton algorithm, by constructing a quadratic approximation model of the objective function, converges quickly without calculating the second derivative, making it suitable for handling such nonlinear optimization problems. The convergence criterion is set as follows: the relative difference between the objective function values ​​of two adjacent iterations is less than a preset minimum. When this criterion is met, it indicates that the objective function is close to its minimum, and the stress tensor tends to stabilize. At this point, iteration stops, and the current stress tensor is output, ensuring a balance between accuracy and efficiency in the inversion results.

[0089] The current stress tensor is decomposed into three eigenvalues ​​and three corresponding eigenvectors. The eigenvector corresponding to the largest eigenvalue is selected as the corrected maximum principal stress direction of the deformable partition. The maximum principal stress directions of all deformable partitions are collected to form the corrected maximum principal stress direction field.

[0090] In detail, the current stress tensor obtained from the solution is subjected to eigenvalue decomposition, yielding three eigenvalues ​​(reflecting the magnitude of the principal stresses) and three corresponding eigenvectors (reflecting the directions of the principal stresses). The eigenvectors are mutually perpendicular, each corresponding to one of the three principal stress directions. The eigenvector corresponding to the largest eigenvalue is the corrected maximum principal stress direction for this deformable partition. The corrected maximum principal stress direction is obtained through weighted event fine-grained inversion at the partition scale, resulting in higher spatial resolution and accuracy compared to previous regional inversion results.

[0091] It should be noted that the corrected direction of the maximum principal stress is a high-precision stress field characteristic quantity.

[0092] The corrected maximum principal stress direction overcomes the limitations of traditional fixed zoning methods, solving the problem of stress field inversion distortion caused by rapid rotation of principal stress under composite roof. Through deformable zoning and secondary inversion design, the interference of stress rotation and interlayer mechanical differences is eliminated, reflecting the true stress direction distribution of the mine more accurately than the initial inversion results.

[0093] In detail, the effect of the corrected direction of the maximum principal stress is as follows: It reflects the changes in the principal stress direction in different time and space regions (especially complex areas such as composite roof and mining-affected areas), and helps to understand the characteristics of stress accumulation, transfer and rotation under mining disturbance.

[0094] Abnormal changes in the principal stress direction (such as rapid rotation or large deflection) are often precursors to rock instability or mine tremors. The corrected direction field can more accurately capture these anomalies.

[0095] By clearly identifying the distribution of high stress directions, the working face advance direction, support strength, or mining pace can be adjusted accordingly to avoid strong disturbances in stress concentration areas (such as areas where the angle between the direction of maximum principal stress and the working face advance direction is too large), thereby reducing the risk of mine tremors.

[0096] By combining the corrected principal stress direction with the bedding structure of the composite roof, the deformation trend of different rock strata under stress (such as whether shear failure occurs along the principal stress direction) can be analyzed, providing a mechanical reference for roof management and tunnel maintenance.

[0097] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A source analysis and inversion system based on mine seismic monitoring, characterized in that, include: The waveform feature extraction module acquires continuous waveforms from mine seismic monitoring and generates continuous parameters characterizing waveform energy changes. The event identification and localization module generates mine seismic events based on the spatiotemporal aggregation patterns of continuous parameters, and determines the spatial location and occurrence time of the mine seismic events; The source parameter inversion module uses the initial motion parameters and amplitude relationships of seismic waves to inversely derive the source mechanism parameters of mining earthquake events. The initial stress field inversion module, based on the source mechanism parameters and the spatiotemporal distribution of mining earthquake events, inverts the direction of the maximum principal stress in the regional stress field; The stress gradient calculation module determines the first spatial direction by combining the mining trajectory of the working face and the second spatial direction by combining the stratigraphic layering information, and calculates the spatial variation rate of the maximum principal stress direction in the first and second spatial directions. The coupling index calculation module obtains the interlayer mechanical difference index based on the mechanical parameters and thickness information of each layer of the composite top plate, and determines the principal stress rotation composite index by combining the two spatial change rates and the interlayer mechanical difference index. The deformable partitioning module generates adaptive deformable partitions along the direction of maximum principal stress based on the principal stress rotation composite index, and allocates and weights mine seismic events according to their spatial position relationship with the partition centerline. The stress field fine inversion module performs secondary stress field inversion based on weighted seismic events within the deformable partition, and outputs the corrected direction of the maximum principal stress.

2. The source analysis and inversion system based on mine seismic monitoring according to claim 1, characterized in that, Acquire continuous waveforms from mine seismic monitoring and generate continuous parameters characterizing waveform energy changes, including: The three-component raw waveform sequences collected from each station are subjected to mean removal and bandpass filtering to obtain standardized waveform sequences; The short-time average energy continuity, long-time average energy continuity, and energy ratio continuity are calculated based on standardized waveform sequences, as follows: The short-time average energy continuity is a time series obtained by averaging the absolute values ​​of a standardized waveform sequence over a first preset time window. The long-term average energy continuity is a time series obtained by averaging the absolute values ​​of a standardized waveform sequence over a second preset time window. The length of the second preset time window is greater than the length of the first preset time window; The energy ratio continuous quantity is a time series obtained by comparing the short-time average energy continuous quantity with the long-time average energy continuous quantity at the same time point; The short-time average energy continuous quantity, long-time average energy continuous quantity, and energy ratio continuous quantity are aggregated according to the station number to form a continuous parameter set.

3. The source analysis and inversion system based on mine seismic monitoring according to claim 2, characterized in that, Mining seismic events are generated based on the spatiotemporal aggregation law of continuous parameters, and the spatial location and occurrence time of the seismic events are determined, including: The spatial coordinates of each station are obtained to form a station spatial location set, and the parameters containing the seismic wave propagation velocity of each rock layer are obtained to form a propagation velocity parameter set. Based on the station spatial location set and the propagation velocity parameter set, the straight-line distance from the candidate source point to each station is calculated, and then the straight-line distance is divided by the seismic wave propagation velocity of the corresponding rock layer to obtain the theoretical propagation time. The energy ratio continuous parameters of each station are time-aligned according to the theoretical propagation time and summed to obtain the spatiotemporal energy aggregation function. Local maxima are searched within the three-dimensional neighborhood of the spatiotemporal energy aggregation function, with each local maximum corresponding to a seismic event. The time coordinates of the local maximums are used as the initial values ​​for the occurrence time of the seismic event, and the spatial coordinates are used as the initial values ​​for the spatial location of the seismic event. Based on the initial values ​​of occurrence time and spatial location, the time corresponding to the maximum value in the energy ratio continuity of each station is extracted as the arrival time observation within a fixed window centered on the theoretical propagation time; Construct the arrival residual for each station; where the arrival residual is the sum of the arrival observation time and the theoretical propagation time; by minimizing the sum of squared arrival residuals for all stations, the occurrence time and spatial location of the seismic event are obtained. Select any two earthquake events from the set of earthquake events whose time interval is less than a first set value or whose spatial distance is less than a second set value to form an earthquake event pair; calculate the difference between the arrival observation difference and the theoretical propagation time difference for each station for the earthquake event pair, as the double-difference residual; by minimizing the sum of squares of the double-difference residuals of all stations, solve for the spatial position correction of each earthquake event; add the corresponding spatial position correction to the spatial position to obtain the updated spatial position of the earthquake event.

4. The source analysis and inversion system based on mine seismic monitoring according to claim 3, characterized in that, By utilizing the initial motion parameters and amplitude relationships of seismic waves, the focal mechanism parameters of mining-related seismic events can be derived inversely, including: Based on the occurrence time and updated spatial location of the seismic event, the theoretical propagation time of the seismic waves to each station is calculated. In the standardized waveform sequence, a fixed length P-wave window and a S-wave window are cut out with the theoretical propagation time as the center. The positive and negative polarities of the first wave of the P-wave are extracted from the P-wave window as the initial motion polarity of the P-wave, and the peak amplitudes are extracted from the P-wave window and the S-wave window as the P-wave amplitude and the S-wave amplitude, respectively. The propagation distance is calculated based on the updated spatial location of the seismic event and the spatial location of the station. The medium attenuation coefficient is calculated based on the quality factor, seismic wave frequency and propagation velocity in the propagation velocity parameter set. The medium attenuation coefficient is inversely proportional to the quality factor, directly proportional to the frequency and inversely proportional to the propagation velocity. Geometric diffusion correction is performed on the P-wave and S-wave amplitudes using the propagation distance, with the correction method being the P-wave amplitude or S-wave amplitude divided by the propagation distance. Attenuation correction is then performed on the corrected P-wave or S-wave amplitudes using the medium attenuation coefficient and the propagation distance, with the correction method being the P-wave amplitude or S-wave amplitude multiplied by the exponential attenuation term of the product of the medium attenuation coefficient and the propagation distance. The amplitude ratio observation is obtained by comparing the S-wave amplitude after the second correction with the P-wave amplitude. An objective function for source mechanism inversion is constructed, which includes: a P-wave first motion polarity consistency term and an amplitude ratio residual term. The P-wave first motion polarity consistency term is the sum of the number of inconsistencies between the observed P-wave first motion polarity at each station and the P-wave first motion polarity predicted by the candidate source mechanism; the amplitude ratio residual term is the sum of the squares of the logarithmic differences between the observed amplitude ratio at each station and the amplitude ratio predicted by the candidate source mechanism. The objective function for source mechanism inversion is minimized using a grid search algorithm. The parameters of the grid search are: strike angle from 0 to 360 degrees, dip angle from 0 to 90 degrees, and slip angle from -90 degrees to +90 degrees, with a preset step size. The strike angle, dip angle, and slip angle of the seismic event are obtained by solving the algorithm. The minimum value of the objective function is used as the sum of squared residuals as a quality index of the source mechanism parameters. The strike angle, dip angle, slip angle, and sum of squared residuals are aggregated to form a set of source mechanism parameters.

5. The source analysis and inversion system based on mine seismic monitoring according to claim 4, characterized in that, Based on the focal mechanism parameters and the spatiotemporal distribution of mining-induced seismic events, the directions of the maximum principal stresses in the regional stress field are obtained through inversion, including: Multiple inversion windows are defined; the temporal and spatial ranges of the inversion windows are set based on the density of seismic events, with the temporal range not exceeding ten times the sampling period of the stations and the spatial range not less than the average spacing of the station distribution; events that occur within the temporal range and whose updated spatial location is within the aforementioned spatial range are selected from the seismic event set to form a subset of inversion events; For each seismic event in the inversion event subset, the fault plane normal vector and the slip direction vector are derived through three-dimensional geometric relationships based on the corresponding strike angle, dip angle, and slip angle. The fault plane normal vector is perpendicular to the fault plane, and the slip direction vector points along the fault plane towards the slip direction. Based on the sum of squared residuals of the seismic events, weighting coefficients are generated through an exponentially monotonically decreasing function. The exponentially monotonically decreasing function is the weighting coefficient being the power of the product of the negative sum of squared residuals of the natural exponential function and the preset scale parameter. An objective function is constructed using the stress tensor as the unknown. The objective function is the sum of the angle mismatch terms for each seismic event, weighted by a weighting coefficient. The angle mismatch term is the square of the angle between the sliding direction vector and the shear stress vector of the event. The shear stress vector is the remaining part after subtracting the normal component of the normal vector of the fault plane from the stress tensor. The objective function is minimized using the least squares optimization algorithm, and the stress tensor is obtained by solving it. The stress tensor is decomposed into three eigenvalues ​​and three corresponding eigenvectors. The eigenvector corresponding to the largest eigenvalue is selected as the direction of the maximum principal stress in the inversion window. The direction of the maximum principal stress is bound to the midpoint of the time range and the geometric center of the spatial range of the inversion window. The maximum principal stress direction field is formed by traversing all inversion windows.

6. The source analysis and inversion system based on mine seismic monitoring according to claim 5, characterized in that, The first spatial direction is determined by combining the mining trajectory of the working face, and the second spatial direction is determined by combining the stratigraphic stratification information. The spatial variation rate of the direction of maximum principal stress in the first and second spatial directions is calculated, including: For the geometric center of the spatial range of the inversion window, the nearest target trajectory point is searched in the mining trajectory set of the working face by calculating the Euclidean distance; the tangential vector is calculated by the coordinate difference between the target trajectory point and the previous and next adjacent trajectory points, and the first spatial direction is obtained by dividing the tangential vector by the modulus. The first spatial direction is the tangential direction of the working face advance. For the geometric center of the spatial range of the inversion window, the corresponding local bedding plane is extracted from the stratigraphic information set. The local bedding plane is described by the coordinates of multiple discrete points. The least squares fitting method is used to fit the discrete point coordinates to obtain the plane equation of the local bedding plane. The normal vector is determined based on the coefficients of the plane equation. The normal vector is divided by the modulus to obtain the second spatial direction, which is perpendicular to the local bedding plane. The angle between the projection of the direction of the maximum principal stress onto the horizontal plane and the north direction is taken as the direction angle, which is defined by a 180-degree period. Angle dewinding is applied to the direction angle sequence to eliminate jumps when crossing zero or 180 degrees. The step size of the inversion window in the first spatial direction and the step size in the second spatial direction are read. The step size is equal to the Euclidean distance between the geometric centers of the spatial extents of adjacent inversion windows in the corresponding directions. The rate of change of the orientation angle sequence is calculated using the central difference method. The rate of change of the inversion window is equal to the difference in orientation angle between the forward adjacent inversion window and the reverse adjacent inversion window divided by twice the step size. For the first and last inversion windows in all inversion windows, forward difference or backward difference is used for calculation. The rate of change of the first inversion window is the difference in orientation angle between the second inversion window and the first inversion window divided by the step size, and the rate of change of the last inversion window is the difference in orientation angle between the last inversion window and the second to last inversion window divided by the step size. The first spatial rate of change sequence and the second spatial rate of change sequence are formed.

7. The source analysis and inversion system based on mine seismic monitoring according to claim 6, characterized in that, Based on the mechanical parameters and thickness information of each layer of the composite roof, the interlayer mechanical difference index is obtained. Combining the two spatial variation rates and the interlayer mechanical difference index, the principal stress rotation composite index is determined, including: Within the vertical range corresponding to the geometric center of the spatial range of each inversion window, the thickness and equivalent Young's modulus of each rock layer within the vertical range are read; the proportion of the thickness of each rock layer to the total thickness of all rock layers is calculated to obtain the proportional weight; the weighted average of the equivalent Young's modulus is calculated based on the proportional weight; where the weighted average is equal to the sum of the products of the equivalent Young's modulus of each rock layer and the corresponding proportional weight; the weighted standard deviation of the equivalent Young's modulus is calculated; where the weighted standard deviation is equal to the square root of the sum of the squares of the differences between the equivalent Young's modulus and the weighted average and the corresponding proportional weight; the ratio of the weighted standard deviation to the weighted average is used as the interlayer mechanical difference index, which is used to characterize the degree of dispersion of the interlayer mechanical properties of the composite roof; For each inversion window, the first spatial rate of change is multiplied by the first spatial step size to obtain the standard first spatial change; the second spatial rate of change is multiplied by the second spatial step size to obtain the standard second spatial change. The sum of squares of the standard first spatial variation and the standard second spatial variation is calculated, and the square root of the result is taken to obtain the comprehensive rotational intensity. The comprehensive rotational intensity is multiplied by the interlayer mechanical difference index to obtain the principal stress rotational composite index. The principal stress rotational composite index sequence is formed by traversing all inversion windows.

8. The source analysis and inversion system based on mine seismic monitoring according to claim 7, characterized in that, Based on the principal stress rotation composite index, adaptive deformable partitions are generated along the direction of maximum principal stress, and seismic events are allocated and weighted according to their spatial relationship with the partition centerline, including: For each inversion window, read the direction of the maximum principal stress corresponding to the inversion window, determine the first spatial step size and the second spatial step size of the inversion window, and take the larger value as the reference length. Using the geometric center of the spatial range of the inversion window as the midpoint, a center line segment with a length equal to the reference length is generated along the direction of the maximum principal stress. The coordinates of the two endpoints of the center line segment are obtained by moving the midpoint in the direction of the maximum principal stress and the opposite direction of the maximum principal stress by half the reference length, respectively. Multiply the reference length by (1 + principal stress rotational composite index) to obtain the partition width; The deformable partition is defined as: a three-dimensional tubular neighborhood consisting of all spatial points whose shortest radial distance from the center line segment does not exceed the width of the partition; where the shortest radial distance refers to the minimum value of the straight-line distances from the spatial points to any point on the center line segment; Events occurring within the inversion window time range are selected from the set of seismic events to form a subset of events within the window. For each seismic event, the shortest radial distance from the updated spatial location of the seismic event to the center line segment is calculated. Using the partition width as the scale parameter, an exponential decay function is used to calculate the weighting coefficients. The weighting coefficient is equal to the square of the ratio of the negative shortest radial distance of the natural exponential function to the partition width. If a seismic event falls into multiple deformable partitions, it is normalized and allocated according to the weight of the seismic event in each deformable partition based on the weighting coefficient of the corresponding deformable partition. That is, the weight of each deformable partition is equal to the corresponding weighting coefficient divided by the sum of the weighting coefficients of all deformable partitions it falls into.

9. The source analysis and inversion system based on mine seismic monitoring according to claim 8, characterized in that, Within the deformable partition, a secondary stress field inversion is performed based on weighted seismic events, outputting the corrected directions of the maximum principal stresses, including: For each deformable partition, read all seismic events within the deformable partition to form an event list; For each mine tremor event in the event list, extract the weighting coefficients based on the sum of squared residuals; And the weighting coefficients obtained based on the shortest straight-line distance; multiply the weighting coefficients and the weighting coefficients to obtain the comprehensive weighting coefficients; bind the comprehensive weighting coefficients to the fault plane normal vector and the slip direction vector of the seismic event to form a set of regional events; A partition objective function is constructed with the stress tensor as the unknown. The partition objective function is the sum of the angle mismatch terms of all mining seismic events in the deformable partition after being weighted by a comprehensive weighting coefficient. The partition objective function is minimized using a quasi-Newton optimization algorithm to obtain the stress tensor of the deformable partition. The convergence criterion is set as the relative difference between the objective function values ​​of two adjacent iterations being less than a preset minimum value. When the convergence criterion is met, the iteration stops and the current stress tensor is output. The current stress tensor is decomposed into three eigenvalues ​​and three corresponding eigenvectors. The eigenvector corresponding to the largest eigenvalue is selected as the corrected maximum principal stress direction of the deformable partition. The maximum principal stress directions of all deformable partitions are collected to form the corrected maximum principal stress direction field.

Images