Big data-based coal mine rock pressure intelligent early warning analysis method and system

Abstract

The application discloses a coal mine rock pressure intelligent early warning analysis method and system based on big data, and particularly relates to the technical field of coal mine safety monitoring. The method collects radar interference phase and navigation satellite observation data in a mining area, removes errors through inter-satellite epoch double-difference solution and Bayesian interpolation, constructs a three-dimensional ground surface displacement rate matrix, synchronously analyzes underground microseismic and electromagnetic detection sequences, extracts a view volume cluster density matrix and a resistivity dynamic gradient matrix, fuses the two matrices to perform grid tomography clustering to intercept a fracture boundary, calibrates overburden three-zone interface coordinates and outputs a dynamic thickness of a curved subsidence zone, substitutes into an elastic limit load equation to calculate a dynamic stress threshold, extracts a deformation partial derivative and an energy increment to splice a hysteresis characteristic vector, solves a node stress accumulation index through a logarithmic logic mapping function, and selects an out-of-limit node to output a three-dimensional early warning coordinate instruction. The system is used to implement the above method.

Description

Technical Field

[0001] This invention relates to the field of coal mine safety monitoring technology, specifically to a method and system for intelligent early warning and analysis of coal mine rockburst based on big data. Background Technology

[0002] During deep coal mining, as the goaf expands, the original stress balance of the overlying strata is disrupted, easily leading to large-scale surface subsidence and rockbursts caused by the rapid release of energy from deep rock masses. Comprehensive analysis of multi-source heterogeneous monitoring data is crucial for understanding the geological structure evolution of mining areas and for early disaster assessment. This typically requires combining surface deformation characteristics with underground rock mass stress changes to predict potential dynamic disaster risks in the mining area.

[0003] However, in actual monitoring of mine stress disasters, different mining geological conditions and subsidence gradients can interfere with the phase calculation of surface interferometric radar. Large-gradient surface deformation can easily lead to severe decoherence in the interferogram, and atmospheric tropospheric delay errors caused by local micro-meteorological changes in the mining area can directly affect the accuracy of deformation monitoring. Existing early warning schemes usually process surface displacement data and underground stress monitoring sequences independently, and do not adequately consider the spatial correlation evolution mechanism between the lag of macroscopic surface deformation and the energy release of microscopic fractures in underground rock masses. This makes it difficult for the data fusion process to accurately map the dynamic topological evolution of overburden caving zones, fracture zones, and bending subsidence zones in goaf areas. Facing the multi-temporal and spatial dimensions of rock mass stress transmission and energy accumulation processes, it is difficult to achieve accurate critical threshold determination, resulting in limitations or even misjudgments in actual disaster early warning results. Therefore, it is still necessary to provide a big data-based intelligent early warning analysis method and system for coal mine rockburst to improve the accuracy and timeliness of dynamic disaster prediction in mining areas. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the embodiments of the present invention provide a method and system for intelligent early warning analysis of coal mine rockburst based on big data to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a big data-based intelligent early warning and analysis method for coal mine rockburst, comprising the following steps: S1, collecting radar interferometric phase map sequences and navigation satellite epoch observation data in the mining area; performing inter-satellite epoch double difference calculation on the navigation satellite epoch observation data to extract tropospheric water vapor delay residuals; mapping the tropospheric water vapor delay residuals to the radar interferometric phase map sequences in the mining area using an empirical Bayesian interpolation algorithm to perform error stripping; and outputting a three-dimensional continuous surface displacement rate matrix; S2, simultaneously acquiring underground microseismic waveform signals and transient electromagnetic detection sequences; analyzing the underground microseismic waveform signals to extract the microseismic apparent volume cluster density matrix; and performing time-series evolution gradient inversion on the transient electromagnetic detection sequences to output a dynamic resistivity gradient matrix; S3, fusing microseismic data... Spatial grid tomographic clustering is performed on the apparent volume cluster density matrix and resistivity dynamic gradient matrix. The attribute abrupt change surface of the resistivity dynamic gradient matrix is ​​extracted to cut the rupture boundary of the microseismic apparent volume cluster density matrix. The three-dimensional interface coordinates of the caving zone, fracture zone and bending subsidence zone are calibrated and the dynamic thickness scalar of the bending subsidence zone is output. S4. The dynamic thickness scalar of the bending subsidence zone is substituted into the elastic limit load equation to calculate the dynamic stress critical threshold. The spatiotemporal partial derivative of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates is extracted and concatenated with the energy increment of the microseismic apparent volume cluster density matrix to form a hysteresis feature vector. The hysteresis feature vector is input into the logarithmic logic mapping function to output the nodal stress accumulation index. Grid nodes with nodal stress accumulation index greater than the dynamic stress critical threshold are screened and three-dimensional early warning coordinate instructions are output.

[0006] In a preferred embodiment, the specific process of acquiring the phase map sequence of radar interferometry in the mining area and the epoch observation data of navigation satellites, and performing inter-ephemeral double-difference calculation on the navigation satellite epoch observation data to extract the tropospheric water vapor delay residual is as follows: Analyze the phase map sequence of radar interferometry in the mining area to extract single-view complex images and construct a spatial registration data stack; analyze the epoch observation data of navigation satellites to extract the dual-frequency carrier phase observation sequence and pseudorange observation parameters; fuse the dual-frequency carrier phase observation sequence and pseudorange observation parameters to construct the inter-station and inter-satellite double-difference observation equations; perform ionospheric dispersion term and satellite clock error term separation calculation on the inter-station and inter-satellite double-difference observation equations to extract the line-of-sight tropospheric delay component; separate the dry delay parameter in the line-of-sight tropospheric delay component to extract the wet delay parameter and output the tropospheric water vapor delay residual.

[0007] In a preferred embodiment, the specific process of mapping the tropospheric water vapor delay residual to the radar interferometric phase map sequence of the mining area using the empirical Bayesian interpolation algorithm to perform error stripping and output a three-dimensional continuous surface displacement rate matrix is ​​as follows: extracting the spatial distribution reference coordinates of the tropospheric water vapor delay residual and fitting a spatial semivariogram; substituting the spatial semivariogram into the empirical Bayesian interpolation algorithm to generate a continuous atmospheric delay phase surface that matches the spatial resolution of the radar interferometric phase map sequence of the mining area; introducing the continuous atmospheric delay phase surface into the spatial registration data stack constructed based on the radar interferometric phase map sequence of the mining area and performing interferometric phase compensation calculation to extract a pure deformation phase map sequence; performing minimum cost flow phase unwrapping and three-dimensional spatial reference projection transformation on the pure deformation phase map sequence to output a three-dimensional continuous surface displacement rate matrix.

[0008] In a preferred embodiment, the specific process of simultaneously acquiring downhole microseismic waveform signals and transient electromagnetic detection sequences, analyzing the downhole microseismic waveform signals to extract the microseismic apparent volume cluster density matrix, and performing time-series evolution gradient inversion on the transient electromagnetic detection sequences to output the resistivity dynamic gradient matrix is ​​as follows: The downhole microseismic waveform signals are subjected to phase initiation picking calculation to extract the source spatial coordinates and radiated energy release parameters; the source spatial coordinates and radiated energy release parameters are mapped to a preset three-dimensional discrete grid, and three-dimensional spatial grid integration is performed to output the microseismic apparent volume cluster density matrix; the transient electromagnetic detection sequences are subjected to transient attenuation curve fitting inversion to generate a three-dimensional apparent resistivity matrix sequence; the three-dimensional apparent resistivity matrix sequence is subjected to spatiotemporal partial derivative analysis along the time evolution axis to extract the local resistivity change rate parameter and assemble it into a resistivity dynamic gradient matrix.

[0009] In a preferred embodiment, the logical process of performing spatial grid tomography clustering by fusing the microseismic apparent volume cluster density matrix and the resistivity dynamic gradient matrix is ​​as follows: extract the density feature parameters of the microseismic apparent volume cluster density matrix and the gradient variation feature parameters of the resistivity dynamic gradient matrix; map the density feature parameters and gradient variation feature parameters to a preset unified three-dimensional spatial grid framework to perform grid node feature dimension expansion and splicing to output a multi-dimensional grid feature tensor; import the multi-dimensional grid feature tensor into the spatial density tomography engine to perform spatial tomography scanning to extract high-dimensional density core points; and divide the three-dimensional cluster set with similar physical response properties by connecting the neighboring grid nodes according to the spatial distribution topology of the high-dimensional density core points.

[0010] In a preferred embodiment, the specific process of extracting the attribute mutation surface of the resistivity dynamic gradient matrix, truncating the rupture boundary of the microseismic apparent volume cluster density matrix, calibrating the three-dimensional interface coordinates of the caving zone, fracture zone, and flexural subsidence zone, and outputting the dynamic thickness scalar of the flexural subsidence zone is as follows: Traversing the edge clusters belonging to resistivity response variation in the three-dimensional cluster set, extracting the longitudinal gradient maxima connected mesh to generate the attribute mutation surface; analyzing the clusters belonging to the core of the microseismic event cluster in the three-dimensional cluster set, extracting the high-frequency energy release dense boundary to generate the rupture boundary; extracting the rupture boundary and topological intersection node set by truncating the rupture boundary in the three-dimensional spatial projection domain through the attribute mutation surface; fitting the three-dimensional interface coordinates of the caving zone, fracture zone, and flexural subsidence zone in layers according to the elevation distribution characteristics of the topological intersection node set; extracting the longitudinal elevation difference between the top plate spatial nodes and the bottom plate spatial nodes of the flexural subsidence zone in the three-dimensional interface coordinates and outputting the dynamic thickness scalar of the flexural subsidence zone.

[0011] In a preferred embodiment, the specific process of substituting the dynamic thickness scalar of the bent subsidence zone into the elastic limit load equation to calculate the dynamic stress critical threshold, and extracting the spatiotemporal partial derivatives of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates and the energy increment of the microseismic apparent volume cluster density matrix to splice into a hysteresis feature vector is as follows: The dynamic thickness scalar of the bent subsidence zone and pre-set rock mass physical parameters are introduced into the analytical elastic limit load equation to perform mechanical limit state solution and output the dynamic stress critical threshold; the deformation time series of the corresponding spatial projection grid nodes are extracted from the three-dimensional interface coordinates in the three-dimensional surface continuous displacement rate matrix, and the spatiotemporal partial derivative parameters are extracted by spatiotemporal partial derivative analysis; the energy increment parameters are extracted from the microseismic apparent volume cluster density matrix within the synchronous time window; the time axis alignment transformation and dimension normalization processing are performed on the deformation spatiotemporal partial derivative parameters and the energy increment parameters, and then spliced ​​together to generate the hysteresis feature vector.

[0012] In a preferred embodiment, the specific process of inputting the hysteresis feature vector into the logarithmic logic mapping function to output the node stress accumulation index, and filtering out grid nodes with node stress accumulation indices greater than the dynamic stress critical threshold and outputting a three-dimensional early warning coordinate command is as follows: extract the nonlinear transfer term of the hysteresis feature vector input into the logarithmic logic mapping function and perform spatiotemporal displacement-stress coupling solution; traverse the three-dimensional spatial grid frame to map the spatiotemporal displacement-stress coupling solution results and output the node stress accumulation index; extract the node stress accumulation index and perform a full-space grid extreme value comparison and judgment with the dynamic stress critical threshold; filter out node stress accumulation indices higher than the dynamic stress critical threshold to lock the target over-limit grid nodes, extract the spatial position coordinate parameters of the target over-limit grid nodes, encapsulate and output the three-dimensional early warning coordinate command.

[0013] The big data-based intelligent early warning and analysis system for coal mine rockburst is used to execute the aforementioned big data-based intelligent early warning and analysis method for coal mine rockburst. It includes: a surface deformation reconstruction module, used to collect radar interferometry phase map sequences and navigation satellite epoch observation data from the mining area; performing inter-satellite epoch double-difference calculation on the navigation satellite epoch observation data to extract tropospheric water vapor delay residuals; mapping the tropospheric water vapor delay residuals to the radar interferometry phase map sequence from the mining area using an empirical Bayesian interpolation algorithm to perform error stripping; and outputting a three-dimensional continuous surface displacement rate matrix; an underground feature analysis module, used to simultaneously acquire underground microseismic waveform signals and transient electromagnetic detection sequences; analyzing the underground microseismic waveform signals to extract the microseismic apparent volume cluster density matrix; performing time-series evolution gradient inversion on the transient electromagnetic detection sequences to output a dynamic resistivity gradient matrix; and spatial clustering. The inversion module is used to fuse the microseismic apparent volume cluster density matrix and the resistivity dynamic gradient matrix to perform spatial grid tomographic clustering, extract the attribute abrupt change surface of the resistivity dynamic gradient matrix to extract the rupture boundary of the microseismic apparent volume cluster density matrix, calibrate the three-dimensional interface coordinates of the caving zone, fracture zone and bending subsidence zone and output the dynamic thickness scalar of the bending subsidence zone; the intelligent inference and early warning module is used to substitute the dynamic thickness scalar of the bending subsidence zone into the elastic limit load equation to calculate the dynamic stress critical threshold, extract the spatiotemporal partial derivative of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates and concatenate it with the energy increment of the microseismic apparent volume cluster density matrix to form a hysteresis feature vector, input the hysteresis feature vector into the logarithmic logic mapping function to output the nodal stress accumulation index, filter the grid nodes whose nodal stress accumulation index is greater than the dynamic stress critical threshold and output the three-dimensional early warning coordinate command.

[0014] The technical effects and advantages of this invention are as follows: (1) A method and system for intelligent early warning analysis of coal mine rockburst based on big data, in the stage of acquiring monitoring characteristics, performs atmospheric delay error stripping on navigation satellite and radar interferometric measurement data through inter-satellite epoch double difference calculation and empirical Bayesian interpolation algorithm, and outputs a high-precision three-dimensional continuous surface displacement rate matrix; simultaneously analyzes the underground microseismic waveform signal and transient electromagnetic detection sequence, extracts the microseismic apparent volume cluster density matrix respectively, and performs time-series evolution gradient inversion to generate a resistivity dynamic gradient matrix. As a result, it can improve the problem of interferometric phase being easily affected by local meteorological interference and large gradient deformation incoherence in the existing monitoring scheme, realize the synchronous and accurate reconstruction of the macroscopic high-resolution spatial deformation field of the surface and the underground microscopic rock mass fracture stress field, and improve the analysis quality of multi-source heterogeneous monitoring data.

[0015] (2) A method and system for intelligent early warning analysis of coal mine rockburst based on big data integrates the microseismic apparent volume cluster density matrix and resistivity dynamic gradient matrix to perform spatial grid tomographic clustering during the early warning simulation stage. It calibrates the three-dimensional interface coordinates of the caving zone, fracture zone, and bending subsidence zone and extracts the dynamic thickness scalar of the bending subsidence zone. It substitutes this into the elastic limit load equation to generate the dynamic stress critical threshold. Furthermore, it spatiotemporal partial derivatives of the three-dimensional surface continuous displacement rate matrix and energy increments are concatenated into a hysteresis feature vector, which is then input into a logarithmic logic mapping function for nodal stress accumulation index comparison. This changes the traditional mechanism of single-point monitoring of independent physical quantities and static fixed threshold early warning. It couples the macroscopic hysteresis deformation of the surface with the stress transmission of underground rock strata in the full spatial grid dimension, improving the accuracy of spatial coordinate positioning and graded early warning of rockburst disasters.

[0016] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0017] Figure 1 This is a flowchart of the intelligent early warning and analysis method for coal mine rockburst based on big data, as described in this invention. Figure 2 This is a schematic diagram of the cross-section of the three-zone interface of the overlying rock and the microseismic fracture boundary in an embodiment of the present invention; Figure 3 This is a schematic diagram comparing the nodal stress accumulation index with the dynamic stress critical threshold warning in an embodiment of the present invention; Figure 4 This is a flowchart of the intelligent early warning and analysis system for coal mine rockburst based on big data, according to the present invention. Detailed Implementation

[0018] This application's embodiments address the problems of low correlation between surface and underground data and susceptibility to environmental noise interference leading to inaccurate disaster prediction results in existing mine pressure monitoring processes by providing a big data-based intelligent early warning and analysis method and system for coal mine rockburst.

[0019] The overall approach of the scheme in this application embodiment is as follows: Error stripping is performed on the surface radar interferometry and navigation satellite observation sequences of the mining area to reconstruct the continuous surface deformation field; density and gradient evolution matrices are extracted by simultaneously analyzing the downhole microseismic and transient electromagnetic data; grid clustering is performed on the downhole data to calibrate the three-zone boundaries of the overburden and solve the dynamic critical threshold; then, the surface deformation partial derivatives and downhole energy increment are extracted and spliced ​​into a feature vector; the stress accumulation state of the grid nodes is calculated through a logical mapping function; and nodes exceeding the limit are screened and early warning coordinate commands are output.

[0020] Example 1; please refer to Figure 1This invention provides a technical solution: a big data-based intelligent early warning and analysis method for coal mine rockburst, comprising the following steps: S1, collecting radar interferometric phase map sequences and navigation satellite epoch observation data in the mining area; performing inter-satellite epoch double difference calculation on the navigation satellite epoch observation data to extract tropospheric water vapor delay residuals; mapping the tropospheric water vapor delay residuals to the radar interferometric phase map sequences in the mining area using an empirical Bayesian interpolation algorithm to perform error stripping; and outputting a three-dimensional continuous surface displacement rate matrix; S2, simultaneously acquiring underground microseismic waveform signals and transient electromagnetic detection sequences; analyzing the underground microseismic waveform signals to extract the microseismic apparent volume cluster density matrix; and performing time-series evolution gradient inversion on the transient electromagnetic detection sequences to output a dynamic resistivity gradient matrix; S3, fusing the microseismic apparent volume... Spatial grid tomographic clustering is performed using the cluster density matrix and resistivity dynamic gradient matrix. The attribute abrupt change surface of the resistivity dynamic gradient matrix is ​​extracted to truncate the rupture boundary of the microseismic apparent volume cluster density matrix. The three-dimensional interface coordinates of the caving zone, fracture zone, and bending subsidence zone are calibrated, and the dynamic thickness scalar of the bending subsidence zone is output. S4. The dynamic thickness scalar of the bending subsidence zone is substituted into the elastic limit load equation to calculate the dynamic stress critical threshold. The spatiotemporal partial derivative of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates is extracted and concatenated with the energy increment of the microseismic apparent volume cluster density matrix to form a hysteresis feature vector. The hysteresis feature vector is input into the logarithmic logic mapping function to output the nodal stress accumulation index. Grid nodes with nodal stress accumulation index greater than the dynamic stress critical threshold are screened and three-dimensional early warning coordinate instructions are output.

[0021] In this implementation plan, step S1 involves preprocessing and spatially fusing multi-source deformation monitoring data of the mining area surface. Specifically, the system simultaneously receives interferometric phase map sequences from radar satellites and epoch files output from ground navigation satellite receiving stations. During this processing stage, inter-satellite epoch double-difference calculation is a numerical processing logic that eliminates common errors such as satellite clock bias and receiver clock bias. The system compares the phase differences between different satellites at different observation time points to separate the signal delay caused by uneven water vapor distribution in the atmosphere. Subsequently, the system calls the empirical Bayesian interpolation algorithm to transform this discrete station-level water vapor delay residual into continuous spatial surface data. This algorithm combines the spatial autocorrelation properties of the data itself with the geographical topological distribution characteristics of known measuring points for unbiased estimation interpolation. The system substitutes the generated continuous error surface into the initial radar interferometric phase map sequence to perform subtraction compensation operations, eliminating atmospheric phase disturbances caused by local micro-meteorological conditions, and calculating the three-dimensional continuous surface displacement rate matrix covering the target mining area.

[0022] Step S2 focuses on the extraction and reconstruction of characteristic parameters of the physical field monitoring of deep rock masses in the mine. The system extracts dynamic feedback parameters during the stress-induced fracturing process of the rock mass through a sensor network deployed in the roadway. The microseismic apparent volume cluster density matrix is ​​a spatial grid data structure constructed by the system based on the three-dimensional spatial coordinates and energy release magnitude of microseismic events. It is used to characterize the density and frequency of microfractures in deep rocks within a specific area. The transient electromagnetic detection sequence records the electromagnetic response evolution data of rock resistivity over mining time. The system performs temporal evolution gradient inversion on the transient electromagnetic detection sequence, that is, extracts the apparent resistivity difference data of nodes at the same spatial location along the time axis, and filters out the dynamic resistivity gradient matrix that reflects changes in water content or conductivity changes during fracture propagation within the rock strata, thus establishing a multi-dimensional feature set for the underlying geological microphysical field.

[0023] Step S3 relies on the downhole physical field feature matrix obtained in the previous step to perform parameter stitching and mechanical boundary delineation in three-dimensional space. The system maps the microseismic cluster density features and resistivity dynamic gradient features under the same time window to a unified three-dimensional grid coordinate system and performs spatial grid tomographic clustering. This is a discrete grid partitioning mechanism based on unsupervised learning. The system divides the space into three-dimensional clusters with different stress attributes according to the vector similarity of two sets of physical parameters on the grid nodes. The system locks the attribute mutation surface with a step change in value in the resistivity dynamic gradient data and uses the spatial projection range of the surface to intercept the spatial boundary of dense microseismic fracture energy release. Through the above topological intersection processing, the system delineates the geometric boundaries of the rock mass collapse area, the delamination fracture development area, and the overall bending area of ​​the overlying strata in the grid model, calibrates the three-dimensional interface coordinate parameters of the caving zone, fracture zone, and bending subsidence zone, and calculates the elevation difference between the top plate grid and the bottom plate grid of the uppermost bending subsidence zone to generate the dynamic thickness scalar of the bending subsidence zone.

[0024] Step S4 integrates macroscopic surface deformation data with microscopic rock layer structure parameters from the well to perform a joint assessment of stress risk levels. The system inputs the dynamic thickness scalar of the bending subsidence zone calculated in step S3 into the elastic limit load equation, and calculates the ultimate stress state parameters that the current overburden spatial structure can withstand based on cantilever beam mechanics logic, using these as adjustable dynamic stress critical thresholds. The system extracts the spatiotemporal partial derivative parameters of the corresponding interface coordinates in the three-dimensional continuous surface displacement rate matrix to obtain the acceleration measure of surface subsidence and the spatial deformation gradient, and combines this with the current energy increment value in the microseismic apparent volume cluster density matrix to assemble a hysteresis feature vector. The hysteresis feature vector reflects the time delay difference and spatial displacement difference in the process of deep fracture stress propagating to the surface. The system substitutes the vector into the logarithmic logic mapping function to perform nonlinear integral derivation, calculates the nodal stress accumulation index representing the current stress concentration degree of a specific grid rock block, traverses and compares the index with the dynamic stress critical threshold in the three-dimensional grid space, filters out abnormal grid nodes that exceed the threshold limit and extracts their coordinate parameters, and encapsulates and outputs a three-dimensional early warning coordinate command.

[0025] Specifically, the process of collecting phase map sequences from radar interferometry in the mining area and epoch observation data from navigation satellites, and extracting tropospheric water vapor delay residuals by performing inter-satellite double-difference calculations on the navigation satellite epoch observation data is as follows: The phase map sequences from radar interferometry in the mining area are analyzed to extract single-view complex images and construct a spatial registration data stack; the epoch observation data from navigation satellites are analyzed to extract dual-frequency carrier phase observation sequences and pseudorange observation parameters; the dual-frequency carrier phase observation sequences and pseudorange observation parameters are fused to construct inter-station and inter-satellite double-difference observation equations; the ionospheric dispersion term and satellite clock error term are separated from the inter-station and inter-satellite double-difference observation equations to extract the line-of-sight tropospheric delay components; the dry delay parameters in the line-of-sight tropospheric delay components are separated to extract the wet delay parameters and output the tropospheric water vapor delay residuals.

[0026] In this implementation plan, the primary and secondary images of the radar interferometric phase map sequence in the mining area are first spatially registered and resampled to generate a single-view complex image containing complex amplitude and phase information. A unified time reference is then used to construct a spatially registered data stack. Simultaneously, the system reads epoch observation data from navigation satellites and extracts the frequency bands... and The system obtains the dual-frequency carrier phase observation sequence and the corresponding pseudorange observation parameters. Based on this, the system selects the ground reference station u and the mining area mobile monitoring station v, as well as the navigation reference satellite s and the target satellite k, to construct the inter-station and inter-satellite double-difference observation equation, the mathematical expression of which is: ;in, This represents the double-difference carrier phase observation value; This represents the double-difference geometric spatial distance between the station and the satellite; This represents the double-difference line-of-sight tropospheric delay component of the target extraction; Indicates the delayed component of the double-difference ionosphere; The wavelength parameter representing the carrier wave transmitted by the navigation satellite; This represents the integer ambiguity parameter of the double difference; This represents the thermal noise error term within the observation system. In the solution process, since the differential operator has already physically eliminated satellite clock errors and receiver clock errors, the system further employs a dual-frequency, ionosphere-free linear combination to eliminate dispersion effects. Fixed fuzziness Then the line-of-sight tropospheric delay component can be separated. For the dry and wet media contained in this delay component, the system calls the static atmospheric pressure model to calculate the zenith dry delay reference, and introduces an elevation angle mapping function to separate the wet delay parameter. The formula is as follows: ;in, This represents the final output tropospheric water vapor delay residual; This represents the zenith dry delay parameter determined based on meteorological experience data; Indicates the corresponding satellite elevation angle Dry atmospheric geometric projection mapping function; This represents the corresponding geometric projection mapping function for the moist atmosphere.

[0027] Specifically, the process of mapping the tropospheric water vapor delay residual to the radar interferometric phase map sequence of the mining area using the empirical Bayesian interpolation algorithm to perform error stripping and output a three-dimensional continuous surface displacement rate matrix is ​​as follows: Extracting the spatial distribution reference coordinates of the tropospheric water vapor delay residual and fitting a spatial semivariogram; substituting the spatial semivariogram into the empirical Bayesian interpolation algorithm to generate a continuous atmospheric delay phase surface that matches the spatial resolution of the radar interferometric phase map sequence of the mining area; introducing the continuous atmospheric delay phase surface into the spatial registration data stack constructed based on the radar interferometric phase map sequence of the mining area and performing interferometric phase compensation calculation to extract a pure deformation phase map sequence; performing minimum cost flow phase unwrapping and three-dimensional spatial reference projection transformation on the pure deformation phase map sequence to output a three-dimensional continuous surface displacement rate matrix.

[0028] In this implementation scheme, the system extracts the tropospheric water vapor delay residuals output by each of the above-mentioned monitoring stations and their corresponding three-dimensional spatial reference coordinates, divides the distance step size within the projected coordinate system of the mining area, and fits the spatial semivariogram function: ;in, This indicates that the spatial distance lag is The estimated value of the semivariance variation over time; Indicates spatial spacing as The total number of valid station data point pairs; Represents the discrete spatial coordinate vector of the g-th navigation satellite station; Indicates coordinate position The tropospheric water vapor delay residual value at the location. After obtaining the variation pattern, it is substituted into the empirical Bayesian interpolation algorithm engine for any pixel coordinate to be estimated in the spatial registration data stack. Generate continuous atmospheric delay phase values ​​matching spatial resolution. Where Z represents the total number of known stations participating in the interpolation calculation within the defined search neighborhood; This represents the coordinates of the z-th known station; These represent the adaptive spatial interpolation weight coefficients. To ensure the minimum interpolation variance, the adaptive spatial interpolation weight coefficients... The method for determining this is to solve the unbiased linear covariance equation system: Under strict adherence Under spatial constraints, combined with Lagrange multipliers Calculate the spatial weight allocation for each station. In the equation Characteristic coordinates and The semivariance between the values. The interpolated continuous atmospheric delay phase surface. Pixel-by-pixel mapping is performed in the spatial registration data stack to perform interferometric phase compensation calculation: ;in, Indicates radar azimuth direction With distance Pure deformation phase at two-dimensional matrix coordinates; This represents the original radar interferometric entanglement phase; This represents the emission wavelength of the radar sensor. Finally, the system constructs a flow network consisting of pixel nodes and connected edges, and calls the minimum cost flow algorithm to... The phase residual span in the network is subjected to global optimal unentanglement to eliminate... The whole-cycle fuzzy constraint is used, and the external digital elevation model of the mining area is used to perform a three-dimensional spatial projection transformation from radar slant range to geographic vertical direction, and the final three-dimensional continuous displacement rate matrix of the ground surface is output.

[0029] Specifically, the process of simultaneously acquiring downhole microseismic waveform signals and transient electromagnetic detection sequences, analyzing the downhole microseismic waveform signals to extract the microseismic apparent volume cluster density matrix, and performing time-series evolution gradient inversion on the transient electromagnetic detection sequences to output the resistivity dynamic gradient matrix is ​​as follows: The downhole microseismic waveform signals are processed by phase initiation picking to extract the source spatial coordinates and radiated energy release parameters; the source spatial coordinates and radiated energy release parameters are mapped to a preset three-dimensional discrete grid, and three-dimensional spatial grid integration is performed to output the microseismic apparent volume cluster density matrix; the transient electromagnetic detection sequences are processed by transient attenuation curve fitting inversion to generate a three-dimensional apparent resistivity matrix sequence; and the three-dimensional apparent resistivity matrix sequence is processed along the time evolution axis by spatiotemporal partial derivative analysis to extract the local resistivity change rate parameter and assemble it into a resistivity dynamic gradient matrix.

[0030] In this implementation scheme, microseismic waveform signals and transient electromagnetic secondary field induced voltage signals returned by monitoring nodes deployed downhole are acquired. For the microseismic waveform signals, the system employs a long-short time window energy ratio algorithm to identify the initial motion to time sequences of P-waves and S-waves. It then combines the travel-time residual equations from multiple observation nodes to perform nonlinear location calculations, extracting the source spatial coordinates of independent microseismic events. And the radiation energy release parameters based on source spectrum analysis ;in, These represent the three-dimensional spatial components of a single microseismic event in an independent local coordinate system downhole. Subsequently, the system constructs a pre-defined three-dimensional discrete mesh framework with a constant step size, targeting any discrete node within the mesh. Set the spatial search radius Perform three-dimensional spatial mesh integration calculations to output the microseismic apparent volume cluster density parameter. ;in, Represents grid nodes The cluster density value at the location; W represents the search radius of the falling space. The total number of microseismic events within the area; This represents the parameter indicating the radiative energy release of the w-th microseismic event; This represents the static shear modulus of the underlying rock mass medium in the well. Indicates the rock mass fracture strain parameter; This represents the Euclidean spatial distance between the earthquake source and the corresponding grid node. The energy space attenuation function coefficients are used to iterate through all nodes and output the microseismic apparent volume cluster density matrix. Simultaneously, the system analyzes the transient electromagnetic detection sequence, performs transient attenuation curve fitting on the secondary field induced voltage curves intercepted at different delay time windows, and calls a transient field inversion algorithm based on damped least squares to map the voltage attenuation rate to apparent resistivity values ​​varying with formation depth, generating a three-dimensional apparent resistivity matrix sequence arranged according to the detection period. ;in, This represents the nth transient electromagnetic detection time epoch. The system applies the central difference operator to the three-dimensional apparent resistivity matrix sequence along the time evolution axis to perform spatiotemporal partial derivative analysis, calculating the local resistivity change rate parameter. ;in, The transient gradient value representing the location of a grid node; This parameter represents the fixed time interval between two adjacent electromagnetic detection cycles. The system arranges the local resistivity change rate parameters of all grid nodes according to the spatial topological dimension and assembles them into a dynamic resistivity gradient matrix.

[0031] Specifically, the logical process of performing spatial grid tomography clustering by integrating the microseismic apparent volume cluster density matrix and the resistivity dynamic gradient matrix is ​​as follows: extract the density feature parameters of the microseismic apparent volume cluster density matrix and the gradient variation feature parameters of the resistivity dynamic gradient matrix; map the density feature parameters and gradient variation feature parameters to a preset unified three-dimensional spatial grid framework, perform grid node feature dimension expansion and splicing to output a multi-dimensional grid feature tensor; import the multi-dimensional grid feature tensor into the spatial density tomography engine to perform spatial tomography scanning to extract high-dimensional density core points; and divide the three-dimensional clusters with similar physical response properties by connecting the neighboring grid nodes according to the spatial distribution topology of the high-dimensional density core points.

[0032] In this implementation scheme, the apparent volumetric density matrix and resistivity dynamic gradient matrix of microseismic activity are read, and each spatial grid element within the matrix is ​​traversed to extract density characteristic parameters representing the degree of fracture density within the rock mass. and gradient variation characteristic parameters characterizing the development and connectivity of rock strata fractures. Considering the spatial resolution differences in acquiring the two heterogeneous physical field data, the system sets up a pre-defined unified three-dimensional spatial mesh framework with a unified coordinate origin and subdivision size. Trilinear interpolation is then used to resample the two feature parameters into this framework. For any mesh node within this framework... The system performs grid node feature dimension expansion and splicing on the one-dimensional density feature parameter and one-dimensional gradient mutation feature parameter mapped to the node, constructing a multi-dimensional grid feature tensor containing multiple physical properties. ;in, This indicates that the sequence number of the three-dimensional spatial axes within the preset unified three-dimensional spatial mesh framework is... The node position; This represents the two-dimensional feature column vector generated by concatenation; and These represent the dimensional normalized weight coefficients for the density and gradient dimensions, respectively. These coefficients are determined by using the sum of the maximum absolute value of the corresponding feature within the entire grid and a preset minimum zero-prevention constant as the divisor, performing a reciprocal operation to eliminate differences in the physical dimensions of multiple sources. Subsequently, the system imports the multidimensional grid feature tensors of all nodes into the spatial density tomography engine. A tomographic clustering core algorithm based on multidimensional feature spatial distance metrics is then used to perform a spatial tomographic scan, calculating the radius of each grid node within its feature space neighborhood. The frequency of neighboring nodes within the range. When a specific grid node... The frequency of neighboring nodes is greater than or equal to the judgment threshold. At that time, it is extracted as a high-dimensional density core point, and the judgment threshold is... The method for determining the dimension is to calculate the product coefficient of the feature tensor of the grid, that is, to take the feature dimension value multiplied by an integer multiple of two. For all extracted high-dimensional density core points, the system examines their relative Euclidean distance in three-dimensional geometric space. According to the spatial distribution topology, for core points whose physical distance is less than the geometric neighborhood radius parameter, graph structure edge connectivity is performed, and cluster expansion is carried out along the connected edges by absorbing neighboring grid nodes with achievable density. Finally, the global grid is divided into multiple independent three-dimensional cluster sets. ;in, This represents the total set of clusters output. denoted as b, representing the b-th independent three-dimensional cluster with similar physical response properties; B represents the total number of discrete clusters automatically divided by spatial tomography clustering.

[0033] Specifically, the process of extracting the attribute mutation surface of the resistivity dynamic gradient matrix, truncating the rupture boundary of the microseismic apparent volume cluster density matrix, calibrating the three-dimensional interface coordinates of the caving zone, fracture zone, and flexural subsidence zone, and outputting the dynamic thickness scalar of the flexural subsidence zone is as follows: Traverse the edge clusters belonging to resistivity response variation in the three-dimensional cluster set, extract the longitudinal gradient maxima connected mesh to generate the attribute mutation surface; analyze the clusters belonging to the core of the microseismic event cluster in the three-dimensional cluster set, extract the high-frequency energy release dense boundary to generate the rupture boundary; extract the rupture boundary in the three-dimensional spatial projection domain through the attribute mutation surface to extract the set of topologically intersecting nodes; fit the three-dimensional interface coordinates of the caving zone, fracture zone, and flexural subsidence zone in layers according to the elevation distribution characteristics of the topologically intersecting node set; extract the longitudinal elevation difference between the top plate spatial nodes and the bottom plate spatial nodes of the flexural subsidence zone in the three-dimensional interface coordinates and output the dynamic thickness scalar of the flexural subsidence zone.

[0034] In this implementation scheme, the edge clusters marked as resistivity response variations in the aforementioned generated three-dimensional cluster set are traversed. For each vertical grid cylinder space within this cluster, discrete grid nodes whose vertical resistivity dynamic gradient values ​​reach local maxima are extracted. The extracted nodes are then topologically connected in the horizontal spatial domain to generate a property mutation surface characterizing the drastic evolution of water content or fractures. Simultaneously, the system analyzes the clusters identified as the core of microseismic event clusters in the three-dimensional cluster set, extracts boundary grid nodes with high-frequency energy release densities exceeding the statistically preset order of magnitude, and generates fracture boundaries reflecting the concentrated fracture profiles of the rock mass. The system utilizes attribute mutation surfaces. Within the three-dimensional spatial projection domain, the fracture boundary Perform Boolean intersection operations to extract the set of topologically intersecting nodes that simultaneously satisfy the dual physical characteristics of resistivity abrupt change and microseismic high-density rupture. Based on the set of intersecting nodes in the topology. Based on the elevation distribution characteristics of each node, the system uses a kernel density estimation algorithm to generate elevation probability density curves. The three minimum probability stationary points of these curves in the vertical distribution are extracted as the stratified elevation critical planes. From bottom to top, the upper limit elevation surface of the lowest caving zone is then fitted sequentially. Upper limit elevation surface of the intermediate fracture zone and the upper limit elevation surface of the topmost curved and sunken zone. The system completes the calibration of the three-dimensional interface coordinates. It also extracts the elevation surface corresponding to the top spatial nodes of each independent grid column within the coverage area of ​​the curved subsidence zone. Elevation surface corresponding to the spatial nodes of the base plate Perform longitudinal elevation difference calculation node by node, and output dynamic thickness scalar of the bending subsidence zone. Combining Figure 2 As shown, the execution logic of spatial tomographic clustering and interface labeling in this embodiment can be intuitively verified. Figure 2 The horizontal axis represents the spatial horizontal grid sequence, and the vertical axis represents the vertical elevation nodes; the circular scattered points in the figure map the discrete microseismic event cluster nodes in the microseismic apparent volume cluster density matrix. By extracting the attribute abrupt change surface to truncate the microseismic rupture boundary, the system fits three interfaces in the vertical space: the bottom dashed line represents the upper limit of the caving zone. The dashed line in the middle represents the upper limit of the fracture zone. The topmost dotted line represents the upper limit of the curved and sunken zone. The microseismic cluster nodes are mainly concentrated between the upper limit of the caving zone and the upper limit of the fracture zone, confirming the technical characteristic of dense rock fractures in this area; while the longitudinal intercept between the upper limit of the flexural subsidence zone and the upper limit of the fracture zone directly represents the dynamic thickness scalar of the flexural subsidence zone output by the system. The physical length provides a reliable geometric parameter basis for subsequent calculation of the dynamic stress threshold.

[0035] Specifically, the process of substituting the dynamic thickness scalar of the bent subsidence zone into the elastic limit load equation to calculate the dynamic stress critical threshold, and extracting the spatiotemporal partial derivatives of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates and the energy increment of the microseismic apparent volume cluster density matrix to form a hysteresis feature vector is as follows: The dynamic thickness scalar of the bent subsidence zone and pre-set rock mass physical parameters are introduced into the analytical elastic limit load equation to perform mechanical limit state solution and output the dynamic stress critical threshold; the deformation time series of the corresponding spatial projection grid nodes are extracted from the three-dimensional interface coordinates in the three-dimensional surface continuous displacement rate matrix, and the spatiotemporal partial derivative parameters are extracted by spatiotemporal partial derivative analysis; the energy increment parameters are extracted from the microseismic apparent volume cluster density matrix within the synchronous time window; and the time axis alignment transformation and dimension normalization processing are performed on the deformation spatiotemporal partial derivative parameters and the energy increment parameters to splice and assemble them to generate the hysteresis feature vector.

[0036] In this implementation scheme, the elastic limit load equation, which incorporates the force balance mechanism of the cantilever beam, is analyzed, and the dynamic thickness scalar of the bending subsidence zone is obtained from the solution. Substituting the pre-set rock mass physical parameters into the equation, the mechanical limit state solution is executed, and the dynamic stress critical threshold is output. ;in, This parameter represents the ultimate tensile strength of the overburden rock beneath the mine foundation. This represents the exposed span parameter of the goaf roof, calculated in real time as mining progresses. Based on the calibrated three-dimensional interface coordinates, the system extracts the deformation time series vector of the corresponding spatial projection grid nodes from the three-dimensional continuous surface displacement rate matrix. The deformation-space-time partial derivative parameters are extracted analytically by performing second-order space-time partial derivative analysis on it. ;in, This indicates the observation time point in the current warning cycle; This represents the Laplace differential operator that measures the spatial deformation gradient. Simultaneously, the system analyzes the microseismic apparent volume cluster density matrix within the current time window, extracts the microseismic energy release difference between two adjacent observation epochs, and generates the energy increment parameter. The system employs a zero-mean standard deviation normalization method for the deformation spatiotemporal partial derivative parameters. and energy increment parameter The process involves performing a time series alignment transformation along the time evolution axis, and then concatenating and assembling the data to generate a hysteresis feature vector. ;in, and These represent the partial derivatives and energy parameters after removing the dimensionless effect, respectively; and The weight of the vector dimension feature concatenation is determined by statistically analyzing the variance contribution rate of two types of physical parameters in the historical measured mine pressure manifestation samples, and using the ratio of the variance contribution rate as the corresponding weight coefficient distribution value.

[0037] Specifically, the process of inputting the hysteresis feature vector into the logarithmic logic mapping function to output the node stress accumulation index, and then filtering out grid nodes with node stress accumulation indices greater than the dynamic stress critical threshold and outputting 3D early warning coordinate commands is as follows: Extract the nonlinear transfer term of the hysteresis feature vector input into the logarithmic logic mapping function and perform spatiotemporal displacement-stress coupling solution; traverse the 3D spatial grid frame to map the spatiotemporal displacement-stress coupling solution results and output the node stress accumulation index; extract the node stress accumulation index and perform a full-space grid extreme value comparison and judgment; filter out node stress accumulation indices higher than the dynamic stress critical threshold to lock the target over-limit grid nodes, extract the spatial position coordinate parameters of the target over-limit grid nodes, encapsulate them, and output the 3D early warning coordinate commands.

[0038] In this implementation scheme, hysteresis feature vectors are extracted. The input is then fed into the nonlinear transit term of the logarithmic logic mapping function to perform a spatiotemporal displacement-stress coupling solution. The solution logic conforms to the formula: ;in, This represents the nodal stress accumulation index in the solution output; The Euclidean L2 norm of the hysteresis eigenvector is calculated. This represents the mapping scale amplification factor set based on the extreme values ​​of past disaster stresses in the mining area; This represents the shape adjustment parameter controlling the slope of the nonlinear stress response in rock mechanics. The system traverses the entire three-dimensional spatial mesh framework, mapping the spatiotemporal displacement-stress coupling solution results to the corresponding meshes in the three-dimensional spatial coordinate system, completing the reconstruction output of the global node stress accumulation index. For any node within the spatial mesh framework, the system extracts its nodal stress accumulation index. The dynamic stress critical threshold obtained from the previous deduction Perform a full-space mesh extremum comparison and determination. The system filters out those that satisfy the condition that the exponential amplitude is greater than the dynamic stress critical threshold. The discrete grid nodes are identified and identified as target grid nodes at risk of rockburst initiation. For each identified target grid node, the system extracts its three-dimensional coordinate parameters of absolute spatial location in the coordinate system. This coordinate parameter sequence, along with the accumulation index exceeding the limit, is then encapsulated into structured data. A three-dimensional early warning coordinate command is output to the coal mine safety monitoring terminal, driving on-site anti-rockburst drilling and directional pressure relief intervention operations. Combined with... Figure 3 As shown, the dynamic early warning mechanism of multi-criticality determination in this system is further clarified. Figure 3The horizontal axis represents the observation time sequence, and the vertical axis represents the calculated characteristic values. The solid line in the figure, showing a smooth downward trend, represents the dynamic stress critical threshold. Its downward trend reflects the progress made by the working face (exposed span) The increasing load (accelerating load) reflects the mechanical and physical law that the bearing capacity of the overburden gradually decreases according to the elastic limit load equation. The sharply fluctuating broken line in the figure represents the nodal stress accumulation index output by the logarithmic logic mapping function. During the full-space grid extreme value comparison and determination process, when the broken line value representing the actual accumulated stress penetrates upward through the smooth threshold line representing the compressive strength limit (as shown by the pin annotation with the limit-crossing mark in the figure), it indicates that the rock mass stress of the corresponding grid has exceeded the elastic limit stage. The system will accurately lock these limit-crossing nodes and trigger a three-dimensional early warning coordinate command. This dynamic comparison mechanism effectively avoids the problem of easy underreporting by traditional static thresholds and achieves intelligent early warning that matches the actual spatiotemporal evolution process of the mine.

[0039] Example 2; please refer to Figure 4 A big data-based intelligent early warning and analysis system for coal mine rockburst, used to execute the big data-based intelligent early warning and analysis method for coal mine rockburst described in the embodiments, includes: a surface deformation reconstruction module, used to collect radar interferometric phase map sequences and navigation satellite epoch observation data in the mining area, perform inter-satellite epoch double difference calculation on the navigation satellite epoch observation data to extract tropospheric water vapor delay residuals, map the tropospheric water vapor delay residuals to the radar interferometric phase map sequences in the mining area using an empirical Bayesian interpolation algorithm to perform error stripping, and output a three-dimensional continuous surface displacement rate matrix; an underground feature analysis module, used to simultaneously acquire underground microseismic waveform signals and transient electromagnetic detection sequences, analyze the underground microseismic waveform signals to extract the microseismic apparent volume cluster density matrix, and perform time-series evolution gradient inversion on the transient electromagnetic detection sequences to output a dynamic resistivity gradient matrix; and a spatial... The clustering inversion module is used to perform spatial grid tomographic clustering by fusing the microseismic apparent volume cluster density matrix and the resistivity dynamic gradient matrix. It extracts the attribute mutation surface of the resistivity dynamic gradient matrix to cut the rupture boundary of the microseismic apparent volume cluster density matrix, calibrates the three-dimensional interface coordinates of the caving zone, fracture zone, and bending subsidence zone, and outputs the dynamic thickness scalar of the bending subsidence zone. The intelligent inference and early warning module is used to substitute the dynamic thickness scalar of the bending subsidence zone into the elastic limit load equation to calculate the dynamic stress critical threshold. It extracts the spatiotemporal partial derivative of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates and concatenates it with the energy increment of the microseismic apparent volume cluster density matrix to form a hysteresis feature vector. The hysteresis feature vector is input into the logarithmic logic mapping function to output the node stress accumulation index. The grid nodes with the node stress accumulation index greater than the dynamic stress critical threshold are screened and three-dimensional early warning coordinate instructions are output.

[0040] In this implementation scheme, the surface deformation reconstruction module serves as the hardware abstraction layer or logical operation node for the system's surface monitoring data processing, directly interacting with external receiving antennas and radar satellite databases. Upon receiving the mine area radar interferometric phase map sequence and navigation satellite epoch observation data, this module invokes its internal double-difference solution engine to perform inter-satellite double-difference solution on the navigation satellite epoch observation data, extracting the tropospheric water vapor delay residual. Subsequently, the module initiates its embedded empirical Bayesian interpolation algorithm program to expand the extracted tropospheric water vapor delay residual in the spatial domain and directly map it to the mine area radar interferometric phase map sequence, completing the error stripping operation for the interferometric phase. After error filtering, the module outputs a three-dimensional continuous surface displacement rate matrix to the system's internal bus, storing it in the underlying shared memory for downstream modules to access.

[0041] The downhole feature analysis module is responsible for feature transformation and data dimensionality reduction of the underlying multiphysics field, and establishes direct communication connections with various underlying sensors deployed in the working face and roadways. Within a preset operating cycle, this module synchronously acquires downhole microseismic waveform signals and transient electromagnetic detection sequences, and initiates parallel analysis threads. In the microseismic signal processing channel, the module analyzes the three-dimensional coordinates and energy parameters of the incoming waveform signals to generate a microseismic apparent volume cluster density matrix. In the electromagnetic signal processing channel, the module processes the voltage decay data of the transient electromagnetic detection sequence epoch-by-epoch, performs time-series evolution gradient inversion operations along the time axis, and outputs a resistivity dynamic gradient matrix. After these two matrices are generated, the module uniformly applies a timestamp and aligns their spatial coordinates before pushing them into the system's data transmission link.

[0042] The spatial clustering inversion module is the core computational unit for the system to perform spatial topology reconstruction of multi-source heterogeneous data. This module reads the microseismic apparent volume cluster density matrix and resistivity dynamic gradient matrix via a link, loads them into a preset 3D spatial grid framework, and performs spatial grid tomographic clustering. After clustering, the module traverses the bottom-level grid data, extracts the attribute abrupt change surface that reaches the extreme value in the resistivity dynamic gradient matrix, and controls this surface within the system projection space to intercept the fracture boundary of the microseismic apparent volume cluster density matrix. Based on the geometric distribution of the intersecting regions, the module calibrates the 3D interface coordinates of the overlying caving zone, fracture zone, and flexural subsidence zone in the grid model. Then, the module extracts the dynamic thickness scalar of the flexural subsidence zone based on the vertical elevation difference of the interface, and passes it forward as a key geometric parameter.

[0043] The intelligent simulation and early warning module undertakes the final mechanical criticality determination and terminal command distribution tasks of the system. This module receives the dynamic thickness scalar of the bending subsidence zone, directly substitutes it into the stored elastic limit load equation to perform calculations, and determines the dynamic stress critical threshold under the current overburden structure state. Simultaneously, based on the three-dimensional interface coordinates, the module extracts the deformation spatiotemporal partial derivatives of the corresponding nodes in the three-dimensional continuous surface displacement rate matrix, and extracts the energy increment from the microseismic apparent volume cluster density matrix. These two types of parameters are then concatenated and assembled into a hysteresis feature vector in the underlying vector space. Subsequently, the module inputs the hysteresis feature vector into a logarithmic logic mapping function to perform nonlinear solutions, outputting the nodal stress accumulation index of the global grid. The module uses an internal comparator to traverse and filter grid nodes with stress accumulation indices greater than the dynamic stress critical threshold, extracts their spatial parameters, and encapsulates and outputs three-dimensional early warning coordinate commands.

[0044] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A method for intelligent early warning and analysis of coal mine rockburst based on big data, characterized in that, Includes the following steps: S1. Collect the phase map sequence of radar interferometry in the mining area and the epoch observation data of navigation satellites. Perform inter-satellite double difference calculation on the epoch observation data of navigation satellites to extract the tropospheric water vapor delay residual. Map the tropospheric water vapor delay residual to the phase map sequence of radar interferometry in the mining area using the empirical Bayesian interpolation algorithm to perform error stripping and output the three-dimensional continuous displacement rate matrix of the Earth's surface. S2. Simultaneously acquire downhole microseismic waveform signals and transient electromagnetic detection sequences, analyze downhole microseismic waveform signals to extract the microseismic apparent volume cluster density matrix, and perform time-series evolution gradient inversion on transient electromagnetic detection sequences to output the resistivity dynamic gradient matrix. S3. Perform spatial grid tomographic clustering by fusing the microseismic apparent volume cluster density matrix and the resistivity dynamic gradient matrix, extract the attribute mutation surface of the resistivity dynamic gradient matrix to cut the rupture boundary of the microseismic apparent volume cluster density matrix, calibrate the three-dimensional interface coordinates of the caving zone, fracture zone and bending subsidence zone and output the dynamic thickness scalar of the bending subsidence zone. S4. Substitute the dynamic thickness scalar of the bending subsidence zone into the elastic limit load equation to calculate the dynamic stress critical threshold. Extract the spatiotemporal partial derivative of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates and the energy increment of the microseismic apparent volume cluster density matrix, and concatenate them into a hysteresis feature vector. Input the hysteresis feature vector into the logarithmic logic mapping function to output the node stress accumulation index. Filter the grid nodes whose node stress accumulation index is greater than the dynamic stress critical threshold and output the three-dimensional early warning coordinate command.

2. The intelligent early warning and analysis method for coal mine rockburst based on big data as described in claim 1, characterized in that: The specific process of collecting phase map sequences from radar interferometry in the mining area and epoch observation data from navigation satellites, and extracting tropospheric water vapor delay residuals by performing inter-satellite double-difference calculations on the navigation satellite epoch observation data is as follows: The single-view complex image was extracted from the phase map sequence of radar interferometry in the mining area and a spatial registration data stack was constructed. The dual-frequency carrier phase observation sequence and pseudorange observation parameters were extracted from the epoch observation data of navigation satellites. Inter-station and inter-satellite double-difference observation equations are constructed by integrating dual-frequency carrier phase observation sequences and pseudorange observation parameters; The ionospheric dispersion term and satellite clock error term were separated and solved for the inter-station and inter-satellite double-difference observation equations to extract the line-of-sight tropospheric delay component. Extract the wet delay parameter from the dry delay parameter in the tropospheric delay component of the separated line of sight, and output the tropospheric water vapor delay residual.

3. The intelligent early warning and analysis method for coal mine rockburst based on big data as described in claim 2, characterized in that: The specific process of mapping the tropospheric water vapor delay residual to the phase map sequence of radar interferometry in the mining area using the empirical Bayesian interpolation algorithm, performing error stripping, and outputting a three-dimensional continuous surface displacement rate matrix is ​​as follows: Extract the spatial distribution baseline coordinates of the tropospheric water vapor delay residual and fit the spatial semivariogram function; Substituting the spatial semivariogram into the empirical Bayesian interpolation algorithm, a continuous atmospheric delay phase surface is generated that matches the spatial resolution of the phase map sequence of radar interferometry in the mining area. In the spatial registration data stack constructed based on the phase map sequence of radar interferometry in the mining area, a continuous atmospheric delay phase surface is introduced to perform interferometric phase compensation calculation to extract a pure deformation phase map sequence; Perform minimum cost flow phase unwrapping and 3D spatial reference projection transformation on the pure deformation phase map sequence to output a 3D continuous displacement rate matrix of the Earth's surface.

4. The intelligent early warning and analysis method for coal mine rockburst based on big data as described in claim 1, characterized in that: The specific process of simultaneously acquiring downhole microseismic waveform signals and transient electromagnetic detection sequences, analyzing the downhole microseismic waveform signals to extract the microseismic apparent volume cluster density matrix, and performing time-series evolution gradient inversion on the transient electromagnetic detection sequences to output the dynamic resistivity gradient matrix is ​​as follows: Perform phase initiation picking and calculation on downhole microseismic waveform signals to extract source spatial coordinates and radiated energy release parameters; The spatial coordinates of the seismic source and the parameters of radiated energy release are mapped to a preset three-dimensional discrete grid. Three-dimensional spatial grid integration is performed to output the apparent volume cluster density matrix of the microseismic event. A three-dimensional apparent resistivity matrix sequence is generated by fitting and inverting the transient attenuation curve of the transient electromagnetic detection sequence. The local resistivity rate parameter is extracted by performing spatiotemporal partial derivative analysis on the three-dimensional apparent resistivity matrix sequence along the time evolution axis and assembled into a dynamic resistivity gradient matrix.

5. The intelligent early warning and analysis method for coal mine rockburst based on big data as described in claim 4, characterized in that: The logical process of performing spatial grid tomographic clustering by fusing the microseismic apparent volume cluster density matrix and the resistivity dynamic gradient matrix is ​​as follows: Extract the density characteristic parameters of the apparent volume cluster density matrix of microseismic events and the gradient variation characteristic parameters of the resistivity dynamic gradient matrix; The density feature parameter and gradient variation feature parameter are mapped to a preset unified three-dimensional spatial mesh framework, and the feature dimension of the mesh nodes is expanded and spliced ​​to output a multi-dimensional mesh feature tensor. Import the multidimensional mesh feature tensor into the spatial density tomography engine to perform spatial tomography scans and extract high-dimensional density core points; Based on the spatial distribution topology of high-dimensional density core points, neighboring grid nodes are connected to divide a set of three-dimensional clusters with similar physical response properties.

6. The intelligent early warning and analysis method for coal mine rockburst based on big data as described in claim 5, characterized in that: The specific process of extracting the attribute abrupt change surface of the resistivity dynamic gradient matrix, truncating the microseismic apparent volume cluster density matrix rupture boundary, calibrating the three-dimensional interface coordinates of the caving zone, fracture zone, and flexural subsidence zone, and outputting the dynamic thickness scalar of the flexural subsidence zone is as follows: Traverse the edge clusters belonging to resistivity response variation in the three-dimensional cluster set, extract the vertical gradient maxima connected mesh to generate attribute mutation surfaces; Analyze the clusters belonging to the core of the microseismic event cluster in the three-dimensional cluster set, and extract the high-frequency energy release dense boundary to generate the rupture boundary; The set of topologically intersecting nodes is extracted by intercepting the broken boundary in the three-dimensional spatial projection domain using a surface with abrupt property changes. Based on the elevation distribution characteristics of the topologically intersecting node set, the three-dimensional interface coordinates of the caving zone, fissure zone, and bending subsidence zone are fitted in layers. Extract the spatial nodes of the top and bottom plates of the curved subsidence zone from the three-dimensional interface coordinates, perform longitudinal elevation difference calculation, and output the dynamic thickness scalar of the curved subsidence zone.

7. The intelligent early warning and analysis method for coal mine rockburst based on big data as described in claim 6, characterized in that: The specific process of substituting the dynamic thickness scalar of the bent subsidence zone into the elastic limit load equation to calculate the dynamic stress critical threshold, and extracting the spatiotemporal partial derivatives of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates and concatenating the energy increments of the microseismic apparent volume cluster density matrix into a hysteresis eigenvector is as follows: The analytical elastic limit load equation introduces a dynamic thickness scalar of the bending subsidence zone and pre-set rock mass physical parameters to perform mechanical limit state solution and output the dynamic stress critical threshold. Based on the three-dimensional interface coordinates, the deformation time series of the corresponding spatial projection grid nodes is extracted from the three-dimensional continuous displacement rate matrix on the Earth's surface, and the spatiotemporal partial derivative parameters of deformation are extracted by spatiotemporal partial derivative analysis. Energy increment parameters are extracted by analyzing the microseismic apparent volume cluster density matrix within the synchronous time window. Perform time axis alignment transformation and dimension normalization on the deformation spatiotemporal partial derivative parameter and the energy increment parameter, and then assemble them to generate hysteresis feature vectors.

8. The intelligent early warning and analysis method for coal mine rockburst based on big data according to claim 7, characterized in that: The specific process of inputting the hysteresis feature vector into the log-logarithmic mapping function to output the node stress accumulation index, filtering out mesh nodes whose node stress accumulation index is greater than the dynamic stress critical threshold, and outputting the three-dimensional early warning coordinate command is as follows: Extract the hysteresis feature vector and input it into the nonlinear transit term of the log-logarithmic mapping function to perform spatiotemporal displacement-stress coupling solution; The solution of spatiotemporal displacement-stress coupling is mapped by traversing the three-dimensional spatial mesh frame, and the output is the nodal stress accumulation index. Extract the nodal stress accumulation index and perform a full-space mesh extreme value comparison and determination based on the dynamic stress critical threshold; Filter nodes with stress accumulation indices exceeding the dynamic stress critical threshold to lock the target out-of-limit grid nodes, extract the spatial location coordinate parameters of the target out-of-limit grid nodes, encapsulate and output the three-dimensional early warning coordinate command.

9. A big data-based intelligent early warning and analysis system for coal mine rockburst, used to execute the big data-based intelligent early warning and analysis method for coal mine rockburst as described in any one of claims 1-8, characterized in that, include: The surface deformation reconstruction module is used to collect phase map sequences of radar interferometry in the mining area and epoch observation data of navigation satellites. It performs inter-epoch double difference calculation on the navigation satellite epoch observation data to extract the tropospheric water vapor delay residual. It then uses the empirical Bayesian interpolation algorithm to map the tropospheric water vapor delay residual to the phase map sequence of radar interferometry in the mining area to perform error stripping and outputs a three-dimensional continuous surface displacement rate matrix. The downhole feature analysis module is used to simultaneously acquire downhole microseismic waveform signals and transient electromagnetic detection sequences, analyze downhole microseismic waveform signals to extract the microseismic apparent volume cluster density matrix, and perform time-series evolution gradient inversion on transient electromagnetic detection sequences to output the resistivity dynamic gradient matrix. The spatial clustering inversion module is used to fuse the microseismic apparent volume cluster density matrix and the resistivity dynamic gradient matrix to perform spatial grid tomographic clustering, extract the attribute change surface of the resistivity dynamic gradient matrix to cut the rupture boundary of the microseismic apparent volume cluster density matrix, calibrate the three-dimensional interface coordinates of the caving zone, fracture zone and bending subsidence zone and output the dynamic thickness scalar of the bending subsidence zone. The intelligent simulation and early warning module is used to substitute the dynamic thickness scalar of the bending subsidence zone into the elastic limit load equation to calculate the dynamic stress critical threshold. It extracts the spatiotemporal partial derivative of the three-dimensional surface continuous displacement rate matrix in the three-dimensional interface coordinates and the energy increment of the microseismic apparent volume cluster density matrix, and concatenates them into a hysteresis feature vector. The hysteresis feature vector is input into the logarithmic logic mapping function to output the node stress accumulation index. The module then filters out grid nodes whose node stress accumulation index is greater than the dynamic stress critical threshold and outputs three-dimensional early warning coordinate instructions.

Images