Coal mine three-dimensional geological modeling method based on machine learning and directional drilling multi-source data

Abstract

This invention discloses a three-dimensional geological modeling method for coal mines based on machine learning and multi-source data from directional boreholes, belonging to the field of three-dimensional geological modeling technology. The method includes: constructing a rock mass stress-damage constitutive relation library containing multiple physical state kernels, each physical state kernel having a mathematical expression describing the relationship between acoustic emission or microseismic response, resistivity change, and stress loading; acquiring microseismic events, identifying them as events reaching the damage kernels, analyzing their spatiotemporal distribution, using the constitutive relation library to reverse-drive the dynamic stress field, and determining the dynamic boundary of the stress shadow area; based on the dynamic boundary of the stress shadow area, virtually deploying computational nodes within it, and inferring virtual microseismic activity and virtual resistivity values ​​according to the constitutive relation library and the stress levels at corresponding points; this invention can achieve dynamic iterative updates of the three-dimensional geological model throughout the entire mining process, accurately reconstructing the entire process of rock mass damage evolution caused by mining-induced stress redistribution.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional geological modeling technology, and more specifically, to a method for three-dimensional geological modeling of coal mines based on machine learning and multi-source data from directional boreholes. Background Technology

[0002] Existing three-dimensional geological models for coal mines are mostly constructed based on static drilling and logging data before mining. They can only reflect the initial geological state before mining. However, as the mining face continues to advance, the stress field of the surrounding rock undergoes dynamic redistribution, triggering a series of mining-induced geological effects such as the evolution of the roof fracture field, the formation of exfoliation water, and the expansion of the floor failure depth. Static geological models cannot characterize the nonlinear dynamic changes under the influence of mining, resulting in a serious disconnect between the model and the actual geological conditions. This makes it impossible to provide effective support for disaster prevention and control during the mining process, and has become a problem that restricts the safe and efficient advancement of intelligent mining.

[0003] Current machine learning-based geological dynamic prediction schemes mostly use historical static data for model training. Their modeling logic is result-oriented and does not consider that the geological evolution caused by mining is a physical process in which stress, energy, and fluid continuously couple in porous media. They lack a mechanism for capturing and representing the driving forces of the process and cannot transform real-time multi-source monitoring data streams such as microseismic monitoring and resistivity monitoring into the dynamic evolution rules inherent in the model. This makes it difficult to achieve accurate time-series prediction of key disaster-causing parameters such as the height of water-conducting fracture zones under the superposition of mining stress. Especially in engineering scenarios where mining faces are continuously advancing and the roof collapses periodically, existing models are easily interfered with and misled by massive amounts of high-noise microseismic event data from mined-out areas behind them, which have no direct mapping relationship with the geological conditions ahead, in the process of pursuing the accuracy of predicting the height of water-conducting fracture zones in unmined areas ahead. This leads to problems such as decreased model generalization ability and prediction results deviating from reality, resulting in the technical dilemma of using chaotic historical data to predict future geological evolution. Summary of the Invention

[0004] To address the problems existing in the prior art, the purpose of this invention is to provide a three-dimensional geological modeling method for coal mines based on machine learning and multi-source data from directional boreholes. This method can achieve dynamic iterative updates of the three-dimensional geological model throughout the entire mining process, accurately reconstruct the entire process of rock mass damage evolution caused by the redistribution of mining stress, and provide technical support for intelligent mining by offering a geological model that is matched in real time with the actual site conditions.

[0005] To solve the above problems, the present invention adopts the following technical solution: A method for 3D geological modeling of coal mines based on machine learning and multi-source data from directional boreholes. This method includes: A rock mass stress-damage constitutive relation library containing multiple physical state kernels is constructed. Each physical state kernel has a mathematical expression describing the relationship between acoustic emission or microseismic response, resistivity change and stress loading. Microseismic events are acquired, identified as events reaching the damage core, their spatiotemporal distribution is analyzed, and the dynamic boundary of the stress shadow zone is determined by using the constitutive relation library to reverse the dynamic stress field. Based on the dynamic boundary of the stress shadow zone, virtual calculation nodes are deployed inside it, and virtual microseismic activity and virtual resistivity values ​​are deduced according to the constitutive relation library and the stress level of the corresponding points. By comparing the virtual resistivity value with the measured resistivity value, the difference is defined as the influence factor of water, and a corrected comprehensive geological state field is output through physical compensation. In the comprehensive geological state field, based on whether the rock mass has reached the damage core threshold, the cumulative amount of water influence factors, and the distance from the dynamic boundary of the stress shadow zone, the evolution potential of water-conducting fractures is calculated and a three-dimensional dynamic cloud map is formed. The evolution potential cloud map of water-conducting fractures is used to adjust the cutting path and advance speed of mining equipment.

[0006] Furthermore, a rock mass stress-damage constitutive relation library containing multiple physical state cores is constructed, including: Variable path loading and acoustic-electric synchronous response tests were performed on the rock core. The waveforms of acoustic emission events and complex impedance spectra under different loading paths were recorded. The inherent response characteristics determined by lithology were separated from the coupled signals by the waveform spectrum feature decoupling algorithm. Based on the inherent response characteristics, the stress values, damage degree and resistivity change values ​​measured in the rock core at different loading stages are projected into the potential energy field to construct a potential function space to describe the physical state and evolution path of the rock mass.

[0007] Furthermore, microseismic events are acquired, identified as events reaching the damage core, and their spatiotemporal distribution is analyzed, including: Waveform features are extracted from the original microseismic waveforms. The extracted waveform features are compared with the theoretical waveform features of the corresponding lithology in the potential function space during the damage core transition. Events that match the comparison are identified as pure damage core microseismic events. Based on pure damaged nuclear microseismic events, the spatiotemporal sequence and energy triggering relationship between events are analyzed, and a directed energy transfer network is constructed, with the edge of the network serving as the dynamic boundary of the stress shadow zone.

[0008] Furthermore, by utilizing the constitutive relation library to reverse the driving stress field, the dynamic boundary of the stress shadow region is determined, including: Based on the energy release nodes in the constructed stress release and transmission network, the failure conditions that satisfy the Mohr-Coulomb criterion at the nodes are taken as the first type of mechanical constraints, and the boundaries of the goaf and the exposed roadway are taken as the second type of mechanical constraints, thus constructing a set of mechanical boundary constraints. Based on the set of mechanical boundary constraints, a rock mass stress-damage constitutive relation library is constructed. With the stress release and transmission network as the skeleton, the stress values ​​of the nodes are recursively calculated along the network, and then a dynamic stress field is constructed by interpolation. Based on the dynamic stress field and rock mass stress-damage constitutive relation library, the ratio of the current stress level to the critical stress level of the damage core in the potential function space is calculated point by point as the stress over / under ratio. The outer boundary of the continuous region where the stress over / under ratio is less than 1 and greater than the preset threshold is tracked to obtain the dynamic boundary of the stress shadow area.

[0009] Furthermore, based on the dynamic boundary of the stress shadow zone, virtual computing nodes are deployed within it, including: Based on the dynamic boundary and dynamic stress field of the stress shadow area, the stress gradient tensor and stress excess / excess ratio of each point inside the stress shadow area are calculated. Regions where the stress gradient amplitude exceeds the threshold are identified and gradient streamline intersection points are extracted as first-class nodes. Points with stress excess / excess ratio greater than the preset first threshold and less than 1 are selected as second-class nodes. The two types of nodes are merged and assigned lithological identifiers and spatial coordinates to form a virtual node layout set. Based on the virtual node deployment set, the corresponding potential function space is retrieved from the rock mass stress-damage constitutive relation library according to the lithological identifier of each node. The initial virtual microseismic activity and virtual resistivity value are derived by taking the real-time stress value at each node as input. Then, spatial consistency correction is performed based on the energy transfer direction and correlation strength between nodes revealed by the stress release and transfer network, and the virtual microseismic activity and virtual resistivity cooperative field is output.

[0010] Furthermore, based on the constitutive relation library and the corresponding point stress levels, virtual microseismic activity and virtual resistivity values ​​are derived, including: For each virtual node, the dynamic stress field time history of the node since entering the dynamic boundary of the stress shadow zone is traced back, the stress loading path is extracted from it, the transition trajectory along the path is simulated in the potential function space, and the endpoint of the trajectory is taken as the real initial damage state of the virtual node. Based on the actual initial damage state and the real-time stress value of the node at the current moment, the preliminary virtual microseismic activity and virtual resistivity value are derived in the potential function space. Then, the actual energy release value of the adjacent nodes in the stress release and transmission network is introduced as a disturbance input for dynamic correction, and the corrected virtual microseismic activity and virtual resistivity value are output. The virtual cumulative released energy characterized by the modified virtual microseismic activity is compared with the measured cumulative released energy of the pure damage core microseismic event at the corresponding spatial location. The energy release deviation is calculated. When the deviation exceeds a preset threshold, the parameters of the potential function space are adjusted using the deviation.

[0011] Furthermore, by comparing the virtual resistivity values ​​with the measured resistivity values, the difference is defined as the influence factor of water, including: The calculated energy release deviation is introduced as a correction coefficient to correct the virtual resistivity value in the virtual microseismic activity and virtual resistivity synergistic field, thus obtaining the reference virtual resistivity field. The reference virtual resistivity field is compared with the measured resistivity field at multiple frequency points. The dispersion effect characteristics and induced polarization characteristics are extracted from the comparison results. The water influence factor vector is defined based on the dispersion effect characteristics and induced polarization characteristics.

[0012] Furthermore, the comprehensive geological state field, corrected through physical compensation, includes: Based on the defined water influence factor vector, the constructed dynamic stress field, and the constructed potential function space, the mechanical parameters of the corresponding lithology in the potential function space are dynamically weakened and corrected to generate a water-bearing damage correction field. Based on the water-bearing damage correction field, combined with the constructed dynamic stress field, the output virtual microseismic activity and virtual resistivity synergistic field, a dynamic comprehensive geological state field is constructed through state vector synthesis rules.

[0013] Furthermore, the evolution potential of the water-conducting fractures is calculated and a three-dimensional dynamic cloud map is generated, including: In the comprehensive geological state field, the ratio of the current stress level at each point to the critical stress level of the damage core is calculated point by point as the proximity of the damage core threshold, the cumulative amount of water influence factor, and the distance to the dynamic boundary of the stress shadow area. The initiation potential of water-conducting fractures is obtained through nonlinear coupling. Based on the initiation potential of water-conducting fractures, and constrained by the vector distribution of rock mass strength and water influence factors in the constructed comprehensive geological state field, the minimum energy dissipation path starting from each initiation potential point is searched, and the path superposition probability is used as the evolution potential of water-conducting fractures to form a three-dimensional dynamic cloud map.

[0014] Furthermore, the evolution potential cloud map of water-conducting fractures is used to adjust the cutting path and advance speed of the mining equipment, including: Based on the three-dimensional dynamic cloud map of the evolution potential of the formed water-conducting fracture, three-dimensional gradient calculation is performed, gradient streamlines are tracked to identify the risk approach direction and the safety avoidance direction, and a risk tendency vector field covering the area in front of the mining is constructed. The cutting path is corrected based on the safe avoidance direction in the risk tendency vector field. At the same time, the allowable disturbance time window to reach the boundary of the high-risk zone is calculated based on the three-dimensional dynamic cloud map of the evolution potential of the formed water-conducting fracture, and the propulsion speed is adjusted accordingly.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: (1) This scheme adopts the technical means of constructing a rock mass stress-damage constitutive relation library and potential function space containing multi-physical state cores, and combining real-time multi-source monitoring data to invert the dynamic stress field of mining. This overcomes the technical problem that existing three-dimensional geological models of coal mines are mostly statically constructed and cannot characterize the nonlinear dynamic geological evolution of rock mass under the influence of mining, resulting in serious disconnect between the model and actual geological conditions and failing to provide effective support for disaster prevention and control. It realizes the dynamic iterative update of the three-dimensional geological model throughout the mining process, accurately restores the entire process of rock mass damage evolution caused by the redistribution of mining stress, and provides the technical effect of providing geological model support that matches the actual site in real time for intelligent mining.

[0016] (2) This scheme adopts the technical means of identifying pure damage core microseismic events to construct a directed energy transfer network, delineating the dynamic boundary of the stress shadow area and setting up virtual computing nodes inside it. This overcomes the technical problems of existing machine learning-based geological prediction schemes being easily interfered with by the massive high noise of monitoring data in the goaf area and having no direct mapping relationship with the geological conditions of the unmined area ahead, resulting in a decrease in the model's generalization ability and deviation of the prediction results from reality. It effectively shields invalid interference data and accurately limits the calculation and prediction range to the critical evolution area of ​​the rock mass affected by mining, greatly improving the model's anti-interference ability and the accuracy of geological evolution prediction in the unmined area.

[0017] (3) This scheme adopts the technical means of constructing a risk trend vector field based on the three-dimensional dynamic cloud map of the evolution potential of water-conducting fractures, correcting the cutting path along the safe avoidance direction, and controlling the advance speed in combination with the allowable disturbance time window. It overcomes the technical problem that the existing coal mine geological modeling and intelligent mining control links are disconnected and the geological evolution prediction results cannot be converted into executable real-time mining control commands. It realizes the closed-loop linkage between mining-induced water inrush risk prediction and early warning and intelligent mining control. While effectively avoiding the risk of water inrush through water-conducting fractures, it also takes into account the working face recovery efficiency, and provides full-process technical support for intelligent, safe and efficient mining of coal mines. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 This is a flowchart of the coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes, as described in this invention. Detailed Implementation

[0020] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0021] Please see Figure 1 A method for three-dimensional geological modeling of coal mines based on machine learning and multi-source data from directional boreholes. This method includes the following steps: Step 1: Construct a rock mass stress-damage constitutive relation library containing multiple physical state kernels. Each physical state kernel has a mathematical expression describing the relationship between acoustic emission or microseismic response, resistivity change, and stress loading. The specific operations are as follows: During the stress loading process under mining disturbance, rock masses simultaneously experience internal damage accumulation, acoustic emission or microseismic event responses, and resistivity characteristic changes. These three phenomena are intrinsically linked to lithology. Each physical state core corresponds to the physical state of a rock mass of a specific lithology at a particular stress loading stage, including independent physical state cores corresponding to the elastic deformation stage, plastic yielding stage, microcrack initiation stage, fracture propagation stage, and macroscopic failure stage. Each physical state core is accompanied by a corresponding mathematical expression. This mathematical expression uses the stress loading on the rock mass as the input independent variable and the acoustic emission or microseismic response characteristic parameters and resistivity change characteristic parameters generated by the rock mass under the corresponding stress state as the output dependent variables. The construction logic of the expression... The system follows the fundamental physical laws of stress-strain-damage in rock mechanics, and matches the measured mapping relationship between the acoustic and electrical responses of corresponding rock masses in laboratory tests and stress loading. During the construction of the constitutive relation library, for various rock masses commonly found in underground coal mines, including sandstone, mudstone, limestone, and coal, multiple sets of physical state kernels and their corresponding mathematical expressions were constructed and fitted and calibrated. All physical state kernels and their corresponding mathematical expressions were stored according to lithology classification, forming a callable and matchable rock mass stress-damage constitutive relation library. Subsequently, based on the lithological information of the rock mass in the target area, the physical state kernels and mathematical expressions of the corresponding lithology can be directly retrieved to complete the quantitative calculation and inversion of the acoustic and electrical response characteristics under the corresponding stress level.

[0022] Step 1 also includes the following steps: Step 11: Perform variable path loading and acoustic-electric synchronous response tests on the rock core, record the acoustic emission event waveforms and complex impedance spectra under different loading paths, and separate the inherent response characteristics determined by the lithology from the coupled signals using a waveform spectrum feature decoupling algorithm. The specific operation is as follows: During the implementation process, undisturbed rock cores are first taken from the corresponding strata of the target coal mine's mining area. These cores are then prepared into standard-sized specimens according to rock mechanics testing standards. Multiple sets of parallel specimens are prepared for the same lithology to eliminate testing errors caused by specimen dispersion. Variable-path loading, conforming to the stress evolution characteristics of underground coal mining, is applied to the prepared core specimens. The loading paths cover various loading modes that match the dynamic changes in surrounding rock stress during underground mining, such as uniaxial stepwise loading, cyclic loading and unloading, and triaxial graded loading. Acoustic emission signal acquisition and complex impedance testing are carried out simultaneously throughout the loading process. Acoustic emission signal acquisition is achieved using acoustic emission sensors placed on the specimen surface, continuously recording the complete waveform data of acoustic emission events generated throughout the loading process. Complex impedance testing is performed using electrodes placed at both ends of the specimen, employing a multi-frequency sweep test mode to continuously record complex impedance spectrum data at different test frequencies corresponding to different loading stages. All acquired data is synchronized with the real-time stress values ​​and loading duration timestamps during the loading process to ensure spatiotemporal consistency of the data.

[0023] The acquired acoustic emission waveform and complex impedance spectrum signals simultaneously contain inherent response characteristics determined by the mineral composition, pore structure, and mechanical properties of the lithology itself, as well as non-inherent interference characteristics caused by vibration of the loading equipment, electromagnetic interference of the test circuit, and friction effects at the sample end. A waveform spectrum feature decoupling algorithm is needed to separate these two types of features. During the implementation of the decoupling algorithm, the acquired acoustic emission waveform and complex impedance spectrum data are first preprocessed. The preprocessing operations cover filtering and denoising, trend term elimination, and baseline correction to obtain standardized time-domain and frequency-domain data. Then, feature decoupling is performed on the preprocessed data. The extracted features cover the amplitude, rise time, duration, energy, dominant frequency, and centroid frequency of the acoustic emission waveform, as well as the real part, imaginary part, phase, amplitude spectrum, and phase spectrum of the complex impedance spectrum. Then, based on the common patterns of feature data measured by multiple parallel samples of the same lithology under the same loading conditions, and the differences in feature data measured by different lithology samples under the same loading conditions, the features are decoupled by independent component analysis. The feature components that only vary with lithology and do not vary with loading path and interference conditions are separated from the coupled feature set. These feature components are the inherent response features of the corresponding lithology and rock mass.

[0024] Step 12: Based on the inherent response characteristics, the stress values, damage degree, and resistivity changes measured at different loading stages of the rock core are projected into the potential energy field to construct a potential function space to describe the physical state and evolution path of the rock mass. The specific operations are as follows: For each set of core samples that have completed testing, the cumulative acoustic emission energy calibration method is used to quantify the degree of rock mass damage at the corresponding loading stage based on the synchronously acquired acoustic emission event data. During the calculation, the total cumulative acoustic emission energy when the sample reaches macroscopic failure is used as the benchmark value. The ratio of the cumulative acoustic emission energy corresponding to a certain loading stage to the benchmark value is used as the degree of rock mass damage at that loading stage. At the same time, for the same loading stage, the resistivity change value relative to the initial state is calculated synchronously. The resistivity change value is characterized by the relative change of the resistivity value corresponding to the characteristic frequency in the complex impedance spectrum measured at that loading stage to the resistivity value corresponding to the same characteristic frequency in the initial unloaded state of the sample. Based on the inherent response characteristics of the corresponding lithology obtained by separation, a basic coordinate system of multidimensional potential energy field is constructed. The dimensions of the coordinate system cover the stress dimension, damage dimension, resistivity dimension, and characteristic dimension corresponding to the inherent response characteristics. Each dimension is set with a corresponding quantitative scale, and the scale range covers the full range of change of the corresponding parameters in the rock mass from the intact state to the macroscopic failure state.

[0025] The stress values ​​measured at different loading stages of the rock core, the quantified damage degree, the calculated resistivity change, and the inherent response characteristic components extracted at the corresponding loading stages are used as a set of multidimensional state parameters. These parameters are projected into the aforementioned multidimensional potential energy field coordinate system to obtain the state points of the rock mass at that loading stage. All state points corresponding to the loading stages form a continuous state evolution trajectory in the potential energy field, which corresponds to the physical state evolution process of the rock mass under the corresponding loading path. Based on the state evolution trajectories obtained from multiple sets of different loading paths for the same lithology, the potential function of the rock mass in the potential energy field is fitted. The physical meaning of the potential function is to characterize the potential energy of the rock mass at a certain state point in the potential energy field. The state evolution always follows the principle of potential energy minimization, that is, the state change of the rock mass during the stress loading process always proceeds along the direction of the fastest decrease in potential energy in the potential energy field. The input parameters of the potential function are the multidimensional state parameters of the rock mass, and the output parameters are the potential energy values ​​of the corresponding state points. Finally, based on the fitted potential function, the constructed multidimensional potential energy field coordinate system, and the state points and evolution trajectories of the corresponding lithology and rock mass at different stress loading stages, a complete potential function space is formed. This potential function space can completely describe all reachable physical states of the corresponding lithology and rock mass, as well as the evolution path of the physical state of the rock mass under different stress loading paths. At the same time, based on the current state parameters of the rock mass, its state change trend in subsequent stress loading processes can be deduced.

[0026] In some embodiments of the present invention, step 2 is further included: acquiring microseismic events, identifying them as events that reach the damage core, analyzing their spatiotemporal distribution, and using the constitutive relation library to reverse the dynamic stress field to determine the dynamic boundary of the stress shadow region. The specific operations are as follows: During the advancement of the mining face, the stress redistribution within the surrounding rock can trigger rock mass damage and fracturing, generating microseismic events. The spatiotemporal distribution and energy characteristics of these microseismic events directly correspond to the stress state and damage evolution process within the surrounding rock. First, effective events directly related to irreversible rock mass damage transitions are screened from the real-time acquired raw microseismic data, while noise interference events unrelated to rock mass damage are eliminated. Then, a network structure representing stress transmission paths is constructed based on the spatiotemporal correlation patterns of the effective events. A mechanical constraint system is built by combining basic rock mechanics principles with the actual engineering boundaries underground. Based on this, a three-dimensional dynamic stress field updated synchronously with the mining progress is obtained through inversion. Finally, the dynamic boundary of the stress shadow zone is delineated based on the stress field distribution and the critical threshold of rock mass damage. This limits the scope of subsequent calculations and analyses to the area where mining stress has already had an impact but the rock mass has not yet undergone macroscopic damage, avoiding interference from invalid data from mined-out areas in the subsequent prediction process and achieving a dynamic and accurate definition of the mining impact range.

[0027] Step 2 also includes the following steps: Step 21: Extract waveform features from the original microseismic waveform. Compare the extracted waveform features with the theoretical waveform features of the corresponding lithology in the potential function space during damage core transition. Events that match the comparison are identified as pure damage core microseismic events. The specific operation is as follows: Based on the raw microseismic waveform data collected in real time by microseismic monitoring stations pre-deployed underground in the coal mine, the deployment of microseismic monitoring stations covers the roadways and chambers around the mining face, and can continuously collect microseismic waveform data generated in the surrounding rock throughout the mining process. All collected waveform data have accurate timestamps and station reception information required for positioning calculation. In the implementation process, the raw waveform data corresponding to a single microseismic event is first preprocessed. The preprocessing operations include filtering and denoising, baseline correction and direct wave picking to obtain standardized and effective waveform data. Then, feature extraction is performed on the standardized waveform data. The extracted waveform features include waveform amplitude, rise time, duration, event energy, dominant frequency, centroid frequency and spectral distribution characteristics. The extracted features are consistent with the acoustic emission waveform feature types collected in the indoor test in step 11 to ensure feature dimension matching in the subsequent comparison process.

[0028] Before feature comparison, the three-dimensional spatial coordinates of the microseismic event are calculated using a microseismic location algorithm. The lithological information of the corresponding underground formation is matched with these coordinates. The theoretical waveform features of the corresponding lithology during damage core transition are retrieved from the potential function space constructed in step 12. The theoretical waveform features corresponding to damage core transition are the standard waveform features corresponding to the inherent response features of the corresponding lithology rock mass when irreversible damage abruptly occurs, which were separated in step 11. The features extracted from the measured microseismic waveforms are matched with the retrieved theoretical waveform features. During the matching calculation, each feature item is standardized and normalized to eliminate matching errors caused by differences in the dimensions of different features. When the calculated similarity value is within the pre-calibrated qualified threshold range, the two sets of features are considered to be consistent, and the corresponding microseismic event is identified as a pure damage core microseismic event. For microseismic events with similarity values ​​exceeding the qualified threshold range, they are identified as interference events unrelated to rock mass damage core transition and are removed, thus completing the screening of valid events.

[0029] Step 22: Based on pure damage core microseismic events, analyze the spatiotemporal sequence and energy triggering relationship between events, construct a directed energy transfer network, and use the edge of this network as the dynamic boundary of the stress shadow zone. The specific operation is as follows: Based on all the pure damage core microseismic events identified in step 21, the correlation analysis between events and the construction of the energy transfer network are completed. First, all pure damage core microseismic events are sequentially ordered according to their timestamps to clarify the temporal relationship of all events. Then, the energy triggering relationship between events is analyzed based on the ordering results. The determination of the energy triggering relationship is based on the stress wave transmission law of rock mass. For two microseismic events with a temporal relationship, the straight-line distance between the spatial coordinates of the two events is calculated first. Combined with the stress wave propagation velocity of the corresponding rock mass, the theoretical time required for the stress wave to travel from the location of the first event to the location of the second event is calculated. When the difference between the timestamp of the second event and the timestamp of the first event is within a reasonable time window corresponding to the theoretical transmission time, the energy triggering relationship is determined. When the energy release of two events conforms to the energy attenuation law of stress transmission in the rock mass, an energy triggering relationship is identified between the two events, with the earlier event being the triggering event and the later event being the triggered event. Based on the energy triggering relationships between all events, a directed energy transfer network is constructed. The nodes in the network correspond to the three-dimensional spatial location of each pure damage core microseismic event, and the directed edges in the network correspond to the connection between two nodes with an energy triggering relationship. The direction of the directed edges points from the node corresponding to the triggering event to the node corresponding to the triggered event, and the weight of the directed edges corresponds to the energy transfer intensity between the two events. After the network is constructed, the outermost nodes of the entire directed energy transfer network are extracted, and the closed contour line formed by connecting the outermost nodes in sequence is used as the initial dynamic boundary of the stress shadow zone.

[0030] Step 23: Based on the energy release nodes in the constructed stress release and transmission network, the failure conditions satisfying the Mohr-Coulomb criterion at the nodes are taken as the first type of mechanical constraints, and the boundaries of the goaf and the exposed roadways are taken as the second type of mechanical constraints. A set of mechanical boundary constraints is constructed, and the specific operations are as follows: In step 22, the energy release nodes in the directed energy transfer network correspond to the pure damage-dependent probability events identified in step 21. The rock mass at each node has undergone irreversible shear failure, and its stress state must satisfy the Mohr-Coulomb criterion for judging rock shear failure in rock mechanics, serving as the first type of mechanical constraint. The derivation logic of the Mohr-Coulomb criterion is based on the physical nature of rock mass shear failure. The condition for rock mass to undergo shear failure is that the shear stress on the shear surface reaches the shear strength of the rock mass itself. The shear strength of the rock mass consists of two parts: cohesion and frictional strength corresponding to the internal friction angle. This criterion quantifies the critical stress condition for rock mass shear failure through the relationship between the maximum and minimum principal stresses, and is a widely used failure judgment criterion in the field of rock mechanics. To obtain the critical maximum principal stress σ1 corresponding to the rock mass satisfying the Mohr-Coulomb shear failure condition, based on the constructed first type of mechanical constraint, the minimum principal stress σ3 at rock mass failure, the internal friction angle φ of the corresponding lithology and rock mass, and the cohesion c of the corresponding lithology and rock mass are required as input parameters. The solution is completed through step-by-step decomposition calculation, as follows: First, calculate the first intermediate result: based on the input internal friction angle φ, calculate the sum of 45° and φ / 2, then calculate the tangent value corresponding to the sum of the angles, square the tangent value, and multiply it by the input minimum principal stress σ3 to obtain the first calculation result; Then, the second intermediate result is calculated: based on the input internal friction angle φ, the sum of 45° and φ / 2 is calculated, the tangent value corresponding to the sum of the angles is obtained, and the tangent value is multiplied by twice the rock mass cohesion c to obtain the second calculation result; Finally, the target critical maximum principal stress σ1 is obtained by summing the first and second calculation results obtained above. The critical maximum principal stress σ1 required for the rock mass to undergo shear failure is finally obtained.

[0031] This formula, as the basic expression of the first type of mechanical constraint, defines the mechanical relationship that the stress state of each energy release node must satisfy. The second type of mechanical constraint is based on the boundary of the goaf formed underground and the boundary of the exposed roadway. The boundary of the goaf and the exposed roadway belongs to the free surface of the rock mass. The normal stress of the free surface is zero, and the tangential stress also meets the mechanical boundary conditions of the free surface. This serves as the fixed boundary constraint for stress field inversion. The first type of mechanical constraint and the second type of mechanical constraint are integrated, and corresponding mechanical boundary conditions are assigned to each constraint node and constraint boundary to form a complete set of mechanical boundary constraints.

[0032] Step 24: Based on the set of mechanical boundary constraints and the constructed rock mass stress-damage constitutive relation library, using the stress release and transfer network as the framework, the stress values ​​of the nodes are recursively calculated along the network, and then a dynamic stress field is constructed through interpolation. The specific operations are as follows: Based on the mechanical boundary constraint set constructed in step 23 and the rock mass stress-damage constitutive relation library constructed in step 1, known constraint nodes close to the working face and goaf boundary are first extracted from the mechanical boundary constraint set. The stress state of these nodes satisfies the Mohr-Coulomb failure condition corresponding to the first type of mechanical constraint and can be used as the starting node for stress recursion calculation. Then, along the directed edge of the stress release and transmission network, the stress value of adjacent nodes is calculated one by one according to the direction of energy transmission. During the recursion calculation, for each node to be calculated, the lithological information of the corresponding strata at the node location is matched first. The physical state kernel and the corresponding mathematical expression of the corresponding lithology are retrieved from the rock mass stress-damage constitutive relation library. Combined with the energy value released by the pure damage kernel microseismic event at the node, the stress level corresponding to the rock mass damage failure at the node is inferred by the mathematical expression of the physical state kernel. At the same time, the limiting conditions of the mechanical boundary constraint set are strictly followed during the calculation process to ensure that the calculation results of each node satisfy the corresponding mechanical constraints.

[0033] After calculating the stress values ​​of all network nodes, the stress values ​​of all nodes are used as known control points in three-dimensional space. Spatial interpolation is used to interpolate the stress values ​​of all spatial points in the entire target calculation area. During the interpolation process, the differences in mechanical parameters of different lithological strata are considered. The interpolation results of different strata are corrected by region. Finally, a three-dimensional dynamic stress field covering the entire mining influence area is obtained. This dynamic stress field can be updated in real time as the mining face advances and new pure damage core microseismic events are identified, accurately characterizing the dynamic redistribution process of the surrounding rock stress field under the influence of mining.

[0034] Step 25: Based on the dynamic stress field and the rock mass stress-damage constitutive relation library, calculate the ratio of the current stress level to the critical stress level of the damage core in the potential function space point by point. Use this ratio as the stress excess / deficit ratio. Track the outer boundary of continuous regions where the stress excess / deficit ratio is less than 1 and greater than a preset threshold to obtain the dynamic boundary of the stress shadow zone. The specific operation is as follows: First, for each spatial calculation point within the coverage area of ​​the three-dimensional dynamic stress field, the lithological information of the corresponding location is matched. The critical stress level of the damage core of the corresponding lithology is retrieved from the potential function space constructed in step 12. The critical stress level of the damage core is the minimum principal stress threshold when the rock mass of the corresponding lithology undergoes irreversible damage core transition. This threshold is calibrated through the indoor test in step 11. For each spatial calculation point, the ratio of the current stress level of the point to the critical stress level of the corresponding lithology damage core is calculated. This ratio is defined as the stress excess / excess ratio. The value of the stress excess / excess ratio can directly characterize the relative relationship between the stress level of the rock mass at that point and the critical damage threshold. When the stress excess / excess ratio is greater than or equal to 1, the current stress level of the point has reached or exceeded the critical stress of the damage core, and the rock mass has undergone irreversible damage. When the stress excess / excess ratio is less than 1, the current stress level of the point has not reached the critical stress of the damage core, and the rock mass has not yet undergone irreversible damage. Irreversible damage mutations occur; based on the geological conditions, lithological characteristics, and mining disturbance patterns of the target mine, a lower limit threshold for the stress excess / deficit ratio is pre-calibrated. Then, all spatial points within the three-dimensional calculation area are traversed and calculated to select all spatial points with a stress excess / deficit ratio less than 1 and greater than the preset lower limit threshold. The continuous region formed by these points is extracted, and the outermost closed boundary of this continuous region is traced. This outer boundary is the final dynamic boundary of the stress shadow area. The continuous region corresponding to the stress shadow area is the region where mining stress has had a significant impact, the rock mass stress level is close to the critical damage threshold, and damage evolution and fracture propagation are highly likely to occur during subsequent mining. By delineating this boundary, the subsequent virtual node layout and state inference work can be completely limited to this region, effectively avoiding the interference of invalid data from damaged goaf areas and stable areas unaffected by mining on the subsequent calculation process.

[0035] In some embodiments of the present invention, step 3 is further included: based on the dynamic boundary of the stress shadow zone, virtual computing nodes are deployed within it, and virtual microseismic activity and virtual resistivity values ​​are deduced according to the constitutive relation library and the corresponding point stress levels. The specific operation is as follows: During the continuous advancement of the mining face, the rock mass within the stress shadow zone is in a critical state of damage evolution. Its subsequent state changes directly determine the initiation and propagation of water-conducting fractures. However, there is insufficient measured monitoring data in this area to directly characterize the internal state of the rock mass. Therefore, it is necessary to achieve a refined quantitative characterization of the physical state of the rock mass in this area through the deployment of virtual nodes and theoretical deduction. The implementation of this step follows a complete logic from node deployment to state deduction and then to closed-loop correction. First, based on the distribution characteristics of the stress field, virtual computing nodes are deployed at key locations within the stress shadow zone to ensure that computing resources are concentrated in highly sensitive areas of damage evolution. Then, combined with the rock mass stress-damage constitutive relation library and stress loading history, the virtual microseismic activity and virtual resistivity value corresponding to each virtual node are deduced. At the same time, the deduction results are corrected based on the microseismic data measured downhole, and the relevant parameters of the potential function space are optimized simultaneously to ensure that the deduction results are consistent with the mechanical response characteristics of the actual rock mass downhole. Finally, a coordinated distribution field that can completely characterize the acoustic and electrical response characteristics of the rock mass within the stress shadow zone is formed.

[0036] Step 3 also includes the following steps: Step 31: Based on the dynamic boundary and dynamic stress field of the stress shadow area, calculate the stress gradient tensor and stress excess / deficit ratio at each point within the stress shadow area. Identify regions where the stress gradient amplitude exceeds a threshold and extract the intersection points of gradient streamlines as the first type of nodes. Select points with a stress excess / deficit ratio greater than a preset first threshold and less than 1 as the second type of nodes. Merge the two types of nodes and assign them lithological identifiers and spatial coordinates to form a virtual node layout set. The specific operations are as follows: First, all spatial calculation points within the stress shadow area are traversed. Based on the spatial distribution data of the three-dimensional dynamic stress field, the stress gradient tensor and stress excess / excess ratio of each spatial point are calculated. The stress gradient tensor is obtained by taking the partial derivatives of the three principal stress components in the three-dimensional dynamic stress field along the three-dimensional spatial coordinate axes. The magnitude of the tensor characterizes the degree of spatial change in the stress state at that point. The stress excess / excess ratio is calculated using the same method as in step 25, which is the ratio of the current stress level at that point to the critical stress level of the corresponding lithological damage core. Based on the lithological characteristics of the target mine and historical mining monitoring data, the threshold for judging the stress gradient amplitude is pre-calibrated. All spatial points within the stress shadow area are traversed to identify areas where the stress gradient amplitude exceeds the preset threshold. Within these areas, stress gradient streamlines are drawn along the direction of the fastest stress growth. The intersection points of the stress gradient streamlines are extracted as the first type of nodes. These nodes correspond to key locations of stress concentration and are key points for subsequent rock mass damage evolution.

[0037] A first threshold for the stress over / under ratio is preset. The value of this threshold is between the lower limit threshold of the stress over / under ratio in step 25 and the value 1. All spatial points within the stress shadow area are traversed, and spatial points with a stress over / under ratio greater than the preset first threshold and less than 1 are selected as the second type of nodes. The rock mass stress level corresponding to this type of node is already very close to the critical stress of the damage core, and irreversible damage transition is about to occur. It needs to be included in the key calculation scope. The selected first type of nodes and second type of nodes are merged, and nodes with duplicate spatial coordinates are removed. Each merged node is matched with a corresponding three-dimensional spatial coordinate. At the same time, the lithological information of the downhole strata corresponding to the spatial coordinate position is matched and assigned a corresponding lithological identifier, and finally a complete set of virtual nodes is formed.

[0038] Step 32: Based on the virtual node deployment set, retrieve the corresponding potential function space from the rock mass stress-damage constitutive relation library according to the lithological identifier of each node. Use the real-time stress value at each node as input to deduce the initial virtual microseismic activity and virtual resistivity value. Then, perform spatial consistency correction based on the energy transfer direction and correlation strength between nodes revealed by the stress release and transfer network, and output the virtual microseismic activity and virtual resistivity cooperative field. The specific operation is as follows: First, traverse each virtual node in the virtual node deployment set. Based on the lithology identifier attached to the node, retrieve the potential function space of the corresponding lithology and the mathematical expression of the corresponding physical state kernel from the rock mass stress-damage constitutive relation library constructed in step 1. Using the real-time stress value at the current moment of the node's spatial location obtained from step 2 as the input independent variable, substitute it into the retrieved physical state kernel mathematical expression. Combining the mapping relationship between the stress state in the potential function space and the acoustic and electrical response characteristics of the rock mass, deduce the initial virtual microseismic activity and the initial virtual resistivity value corresponding to the node. The virtual microseismic activity is used to quantify the probability of a microseismic event occurring at the node under the current stress state, the event frequency per unit time, and the expected energy released. The virtual resistivity value is used to characterize the theoretical resistivity parameters of the rock mass corresponding to the node under the current stress and damage state.

[0039] After completing the initial value derivation for all nodes, based on the energy transfer direction and correlation strength between nodes revealed by the stress release and transfer network constructed in step 22, spatial consistency correction is performed on the initial derivation results. During the correction process, the adjacent nodes of each virtual node in the stress release and transfer network, as well as the correlation strength weights corresponding to energy transfer between nodes, are first determined. Combining the initial derivation results and correlation strength weights of adjacent nodes, the initial virtual microseismic activity and initial virtual resistivity values ​​of that node are smoothly corrected to ensure that the parameter changes of adjacent nodes conform to the mechanical distribution law of the continuous rock mass, eliminating the spatial parameter abrupt change problem caused by isolated calculations. After completing the correction of all nodes, using the corrected parameters of all virtual nodes as spatial control points, a three-dimensional spatial interpolation method is used for global interpolation calculation, ultimately obtaining the virtual microseismic activity and virtual resistivity synergistic field covering the entire stress shadow area.

[0040] Step 33: For each virtual node, trace back the dynamic stress field time history of the node since entering the dynamic boundary of the stress shadow zone, extract the stress loading path, simulate the transition trajectory along the path in the potential function space, and take the endpoint of the trajectory as the true initial damage state of the virtual node. The specific operation is as follows: As the mining face advances, the dynamic boundary of the stress shadow zone continuously expands. Each spatial point within the stress shadow zone gradually moves from a stable area unaffected by mining activity into the stress shadow zone as the face progresses. Its stress state exhibits continuous changes over time. The damage accumulation process of the rock mass is directly related to the entire stress loading path; the current stress value alone cannot accurately characterize its true initial damage state. During implementation, for each virtual node in the centralized virtual node deployment, the dynamic stress field data corresponding to all time steps from the moment the node entered the dynamic boundary of the stress shadow zone to the current calculation time is traced back based on the node's three-dimensional spatial coordinates. The three-dimensional principal stress value corresponding to each time step of the node is extracted, forming a complete value for that node. The stress loading path is fully reproduced, representing the complete stress change process experienced by the rock mass after entering the mining influence zone. The extracted stress loading path is substituted into the potential function space of the corresponding lithology of the node, and the transition trajectory of the rock mass state along the loading path is simulated in the potential function space. The simulation process follows the principle of potential energy minimization in rock mechanics, that is, the state transition of the rock mass in each stress change process proceeds along the direction of the fastest decrease in potential energy in the potential function space. The simulation process starts from the complete initial state of the rock mass corresponding to the moment when the node enters the stress shadow zone, and completes the simulation calculation of the state transition step by step along the stress loading path. Finally, the endpoint of the transition trajectory is obtained. The state parameters in the potential function space corresponding to the endpoint are the real initial damage state of the virtual node at the current calculation moment.

[0041] Step 34: Based on the actual initial damage state and the real-time stress value of the node at the current moment, the preliminary virtual microseismic activity and virtual resistivity value are derived in the potential function space. Then, the actual energy release value of the adjacent nodes in the stress release and transmission network is introduced as a disturbance input for dynamic correction, and the corrected virtual microseismic activity and virtual resistivity value are output. The specific operation is as follows: Using the true initial damage state of each virtual node obtained in step 33 as the starting point, and combining the real-time stress value at the current spatial location of the node, the state change process of the rock mass at the node from the true initial damage state under the current real-time stress is deduced in the potential function space of the corresponding lithology, thus obtaining the final damage evolution result of the node under the current stress state. Based on this damage evolution result, and combining the mathematical expression of the corresponding physical state kernel in the rock mass stress-damage constitutive relation library constructed in step 1, the preliminary virtual microseismic activity and preliminary virtual resistivity value corresponding to the node are calculated. After completing the preliminary value calculation, the true energy release value of adjacent nodes in the stress release and transmission network is introduced as a disturbance input to dynamically correct the preliminary calculation results. The values ​​are derived from the measured energy data of the pure damage core microseismic events identified in step 21. Stress disturbances and damage evolution at adjacent locations in the rock mass have a direct mutual influence. The actual energy release that has occurred at adjacent nodes will change the stress environment around the virtual node, thereby affecting its damage evolution process and acoustic-electric response characteristics. In the correction process, the disturbance weight corresponding to the actual energy release value of each adjacent node is first determined according to the energy transfer direction and correlation strength between the virtual node and adjacent nodes in the stress release and transmission network. Then, the weighted disturbance value is substituted into the preliminary calculation results to dynamically correct the preliminary virtual microseismic activity and preliminary virtual resistivity values, eliminate the calculation deviation caused by damage disturbances that have occurred in adjacent areas, and finally output the corrected virtual microseismic activity and virtual resistivity values.

[0042] Step 35: Compare the virtual cumulative released energy characterized by the corrected virtual microseismic activity with the measured cumulative released energy of the corresponding spatially pure damage core microseismic event, calculate the energy release deviation, and when the deviation exceeds a preset threshold, adjust the parameters of the potential function space using this deviation. The specific operation is as follows: First, for each virtual node in the virtual node deployment set, based on the corrected virtual microseismic activity output in step 34, calculate the virtual cumulative released energy at the spatial location corresponding to that node. The calculation process is the sum of the expected released energy represented by the virtual microseismic activity at all time steps from the moment the node enters the stress influence zone to the current calculation moment. Simultaneously, acquire the measured energy values ​​of all pure damage core microseismic events identified in step 21 at the three-dimensional spatial location of the node, and sum all the measured energy values ​​to obtain the measured cumulative released energy at that spatial location. Energy; Based on the virtual cumulative released energy and the measured cumulative released energy, the energy release deviation is calculated. The energy release deviation is used to quantify the relative deviation between the theoretical deduction results and the downhole measured results. The energy release deviation D needs to be calculated. This calculation process uses two physical quantities determined in the aforementioned steps as input parameters. One is the virtual cumulative released energy Ev represented by the corrected virtual microseismic activity at the corresponding spatial location. The other is the measured cumulative released energy Em obtained from the statistics of pure damage core microseismic events at the corresponding spatial location. The specific step-by-step derivation and calculation process is as follows: The first step is to calculate the absolute deviation of energy release. This involves calculating the difference between the input virtual cumulative energy release Ev and the measured cumulative energy release Em, and then taking the absolute value of the difference to obtain the absolute deviation of energy release. This step eliminates the interference of the positive or negative direction of the deviation on the evaluation results, retaining only the deviation between the theoretical deduction and the field measurement results. This ensures that the final calculated deviation can uniformly and objectively represent the degree of deviation of the deduction results. The second step is to determine the baseline value for deviation assessment. The measured cumulative released energy Em is used as the baseline value for this deviation assessment. This baseline value corresponds to the actual energy release result of irreversible damage to the rock mass. Using this as a baseline can eliminate the influence of the difference in energy release magnitude of the rock mass at different spatial locations on the deviation assessment result, and ensure that the deviation of different calculation points has uniform horizontal comparability.

[0043] The third step is to summarize and calculate the energy release deviation D. The absolute deviation of energy release calculated in the first step is used as the numerator, and the measured cumulative energy release benchmark value determined in the second step is used as the denominator. A division operation is then performed to finally obtain the energy release deviation D for the corresponding spatial location.

[0044] The derivation of this calculation is based on the fundamental principles of error analysis. Using the measured cumulative energy release as a benchmark, the relative deviation is obtained by calculating the ratio of the absolute deviation between the virtual and measured values ​​to the benchmark value. This relative deviation eliminates the influence of differences in energy release magnitude at different spatial locations, ensuring the objectivity and universality of the deviation assessment. Based on the accuracy of the monitoring data and the lithological characteristics of the target mine, a threshold for the energy release deviation is pre-set. When the calculated energy release deviation exceeds this preset threshold, it indicates a deviation between the parameters within the current parameter range and the mechanical response characteristics under actual underground conditions. The parameters within the parameter range need to be adjusted using this deviation. The adjustment process aims to minimize the energy release deviation, iteratively optimizing the correlation coefficients of the parameter range until the deviation between the derived virtual cumulative energy release and the measured cumulative energy release within the adjusted parameter range falls within the preset threshold range. This process achieves dynamic self-correction of the constitutive model, ensuring the matching degree between subsequent derivation results and the actual underground geological conditions.

[0045] In some embodiments of the present invention, step 4 is further included: comparing the virtual resistivity value with the measured resistivity value, defining the difference as the water influence factor, and outputting a corrected comprehensive geological state field through physical compensation. The specific operation is as follows: The resistivity variation of rock mass is jointly controlled by two independent factors: one is the change in pore structure caused by stress loading and internal damage evolution, and the other is the state of groundwater in the pores and fractures of the rock mass. The virtual resistivity value derived in step 3 is generated based on the stress-damage constitutive relation library of rock mass and only reflects the effect of stress and damage on rock mass resistivity, without including the additional influence of groundwater. The measured resistivity field collected in the well is the result of the combined effect of stress, damage and water. The implementation logic of this step is to compare the two types of resistivity fields, isolate the resistivity difference caused only by water, quantify this difference as the influence factor of water, and then perform physical compensation correction on the rock mass mechanical parameters based on the weakening effect of groundwater on the mechanical properties of the rock mass. Finally, the multi-dimensional physical field data are integrated to construct a dynamic comprehensive geological state field. This state field can completely characterize the stress state, damage degree and water-bearing characteristics of the rock mass under the influence of mining.

[0046] Step 4 also includes the following steps: Step 41: Introduce the calculated energy release deviation as a correction coefficient to correct the virtual resistivity value in the virtual microseismic activity and virtual resistivity combined field, thus obtaining the reference virtual resistivity field. The specific operation is as follows: The virtual resistivity value in the virtual microseismic activity and virtual resistivity co-field generated in step 3 is based on the extrapolation results of the potential function space. The energy release deviation calculated in step 35 can quantify the relative deviation between the potential function space extrapolation results and the actual damage state of the downhole rock mass. This deviation will be directly transmitted to the calculation process of the virtual resistivity value. Therefore, the virtual resistivity value needs to be corrected first to eliminate the systematic error brought about by the model extrapolation and obtain a benchmark virtual resistivity field that only reflects the stress and damage effects. During implementation, for each spatial calculation point in the stress shadow zone, the energy release deviation corresponding to that point is matched first. The energy release deviation is used as a correction coefficient to correct the virtual resistivity value of the corresponding point in the virtual resistivity co-field point by point. The correction amplitude is positively correlated with the value of the energy release deviation. The larger the value of the energy release deviation, the larger the correction amplitude of the corresponding virtual resistivity value. The correction process follows a linear mapping law between the degree of rock mass damage and the change in resistivity. The higher the degree of rock mass damage and the higher the degree of internal fracture development, the greater the corresponding resistivity change. The two show a stable linear correlation. The energy release deviation directly quantifies the relative deviation of the damage degree estimation. Therefore, the relative deviation of resistivity estimation is linearly correlated with the energy release deviation. Based on this linear relationship, a correction formula is constructed. In order to eliminate the virtual resistivity systematic error caused by the deviation of rock mass damage estimation, a benchmark virtual resistivity value ρb that only reflects the effect of stress and damage on rock mass resistivity is obtained. Three physical quantities determined in the aforementioned steps are required as input parameters: the virtual resistivity value ρv to be corrected at the corresponding spatial location, the linear proportionality coefficient k of the same lithology calibrated by indoor core tests, and the energy release deviation D calculated in step 35 at the corresponding spatial location. The specific step-by-step derivation and calculation process is as follows: The first step is to calculate the core term of the deviation correction. The input linear scaling factor k is multiplied by the energy release deviation D at the corresponding spatial location to obtain the core term of the deviation correction. This term is used to quantify the degree of damage deviation caused by the model extrapolation and the corresponding resistivity relative correction magnitude, which provides the basis for the subsequent calculation of the comprehensive correction coefficient. The second step is to calculate the comprehensive correction coefficient. Using constant 1 as the base value, subtract the deviation correction core term obtained in the first step to obtain the comprehensive correction coefficient. This coefficient is a dimensionless multiplier value, which fully represents the overall correction multiplier required for the virtual resistivity value at the corresponding spatial location, and can offset the systematic error brought about by the model derivation. The third step is to summarize and calculate the baseline virtual resistivity value ρb, and then use the input virtual resistivity value ρv to be corrected as the basis for calculation. Multiply it with the comprehensive correction coefficient calculated in the second step to finally obtain the baseline virtual resistivity value ρb after eliminating system errors.

[0047] The derivation of this formula is based on the linear mapping relationship between rock mass damage and resistivity, as well as the basic law of error propagation. The deviation in the derivation of the degree of rock mass damage will be linearly propagated to the derivation result of resistivity. By measuring the deviation through energy release and constructing a linear correction term, the systematic error brought about by the model derivation can be eliminated, ensuring that the corrected resistivity value only reflects the effect of stress and damage. After the virtual resistivity values ​​of all spatial points are corrected, a reference virtual resistivity field covering the entire stress shadow area is formed.

[0048] Step 42: Compare the reference virtual resistivity field with the measured resistivity field at multiple frequency points. Extract dispersion effect features and induced polarization features from the comparison results. Define the water influence factor vector based on the dispersion effect features and induced polarization features. The specific operations are as follows: The underground measured resistivity field was acquired through a mine multi-frequency resistivity monitoring system. Sensors of the system were deployed in directional boreholes and roadways around the mining face. Data acquisition was completed using a multi-frequency sweep test mode, obtaining three-dimensional spatial resistivity distribution data at different test frequencies, forming a multi-frequency measured resistivity field. Specifically, first, the three-dimensional spatial registration of the reference virtual resistivity field and the measured resistivity field was completed, ensuring complete alignment of the spatial coordinate systems of the two fields and a one-to-one correspondence of the three-dimensional coordinates of each spatial calculation point, eliminating comparison errors caused by spatial position deviations. After spatial registration, for each spatial calculation point, the reference virtual resistivity value and the measured resistivity value were compared at all test frequency points to obtain the resistivity difference sequence corresponding to each frequency point. The resistivity difference obtained from the multi-frequency comparison was then analyzed. In the resistivity difference sequence, dispersion effect characteristics and induced polarization characteristics are extracted. The dispersion effect characteristics are the amplitude and variation law of resistivity value with test frequency, and the induced polarization characteristics are the parameters related to the attenuation of polarization potential generated by the rock mass under the action of alternating electric field, covering charging rate, time constant and frequency correlation coefficient. The resistivity change caused by stress and damage alone will not produce obvious dispersion effect and induced polarization effect. The changes of the two types of characteristics are completely determined by the occurrence state of groundwater in the rock mass. Based on the extracted dispersion effect characteristics and induced polarization characteristics, a water influence factor vector is defined. The multiple dimensions of the water influence factor vector correspond to the rock mass water content, the degree of fissure water filling and the groundwater mineralization, respectively. The magnitude of the vector characterizes the comprehensive influence of water on the rock mass resistivity, and the direction of the vector characterizes the dominant type of water influence.

[0049] Step 43: Based on the defined water influence factor vector, the constructed dynamic stress field, and the constructed potential function space, the mechanical parameters of the corresponding lithology in the potential function space are dynamically weakened and corrected to generate a water-bearing damage correction field. The specific operations are as follows: The presence of groundwater softens rock masses, reducing their cohesion, internal friction angle, elastic modulus, and Poisson's ratio, thus accelerating the damage evolution process. The intensity of this softening is directly related to the magnitude of the water influence factor vector; a larger magnitude indicates a greater weakening of the rock mass's mechanical parameters. During implementation, for each spatial calculation point within the stress shadow zone, the water influence factor vector, the current dynamic stress field's three-dimensional principal stress value, and the potential function space of the corresponding lithology are first matched. The initial mechanical parameters of the corresponding lithology and rock mass are then retrieved from the potential function space. The initial mechanical parameters are those obtained from the indoor core test in step 11 under dry rock conditions. Based on the modulus of the water influence factor vector and combined with the indoor water softening test results of the corresponding lithology and rock mass, the weakening coefficient corresponding to each mechanical parameter is determined. The value of the weakening coefficient is negatively correlated with the modulus of the water influence factor vector; the larger the modulus of the water influence factor vector, the smaller the value of the weakening coefficient, and the greater the correction range of the corresponding mechanical parameter. Using the determined weakening coefficient, the initial mechanical parameters of the corresponding lithology are dynamically weakened and corrected point by point to obtain the corrected mechanical parameters of each spatial point under water-bearing conditions. Based on the corrected mechanical parameters of all spatial points within the stress shadow area and combined with the spatial distribution characteristics of the dynamic stress field, a water-bearing damage correction field covering the entire stress shadow area is generated. This correction field fully characterizes the weakening influence of groundwater occurrence on the mechanical properties and damage evolution process of the rock mass.

[0050] Step 44: Based on the water-bearing damage correction field, and combining the constructed dynamic stress field, the output virtual microseismic activity, and the virtual resistivity synergistic field, a dynamic comprehensive geological state field is constructed through state vector synthesis rules. The specific operations are as follows: The dynamic integrated geological state field is a complete quantitative representation of the multi-physical field state of the rock mass under the influence of mining. It needs to integrate physical information from multiple dimensions such as stress state, damage evolution, and water occurrence. First, for each spatial calculation point in the stress shadow area, the corrected mechanical parameters of the water-bearing damage correction field, the three-dimensional principal stress value of the dynamic stress field, and the microseismic activity parameters and resistivity parameters in the virtual microseismic activity and virtual resistivity synergistic field are extracted. Each parameter is treated as an independent state component to construct the multi-dimensional state vector corresponding to that point. Based on the pre-set state vector synthesis rules, the multi-dimensional state vector of each spatial point is synthesized. The synthesis rules follow the basic physical laws of rock mechanics and hydrogeology to ensure that the synthesized state parameters can fully reflect the comprehensive geological state of the rock mass at that point, covering stress level, damage degree, water-bearing state, and mechanical properties. After the state vector synthesis of all spatial points in the stress shadow area is completed, a dynamic integrated geological state field covering the entire stress shadow area is formed. This state field can be dynamically updated in real time as the mining face advances and the on-site monitoring data is updated, accurately representing the dynamic evolution process of the rock mass geological state under the influence of mining.

[0051] In some embodiments of the present invention, step 5 is further included: in the comprehensive geological state field, based on whether the rock mass has reached the damage core threshold, the cumulative amount of water influence factors, and the distance to the dynamic boundary of the stress shadow zone, the evolution potential of water-conducting fractures is calculated and a three-dimensional dynamic cloud map is formed. The specific operation is as follows: Following the dynamic integrated geological state field constructed in step 4, the quantitative calculation of the evolution potential of water-conducting fractures and the construction of a three-dimensional dynamic cloud map are completed. The initiation and expansion of water-conducting fractures are the direct causes of water inrush disasters during underground mining in coal mines. Their evolution process is jointly controlled by the degree of rock mass damage, the state of groundwater occurrence, and the evolution process of mining stress. The complete logic from quantifying the initiation risk to deducing the expansion trend is as follows: First, based on the multi-dimensional parameters in the integrated geological state field, the possibility of water-conducting fracture initiation in the rock mass is quantified point by point. Then, based on the basic principles of rock mass fracture mechanics, the expansion path and distribution probability after fracture initiation are deduced. Finally, the spatial distribution characteristics of water-conducting fracture evolution are transformed into a three-dimensional dynamic cloud map, which intuitively presents the evolution risk distribution of water-conducting fractures in the mining area.

[0052] Step 5 also includes the following steps: Step 51: In the comprehensive geological state field, calculate the ratio of the current stress level to the critical stress level of the damage core at each point as the proximity of the damage core threshold, the cumulative amount of water influence factor, and the distance to the dynamic boundary of the stress shadow zone. Obtain the initiation potential of water-conducting fractures through nonlinear coupling. The specific operation is as follows: First, each spatial calculation point within the stress shadow area covered by the comprehensive geological state field is traversed, and the values ​​of three key parameters are calculated point by point. The first parameter is the damage core threshold proximity, which is calculated as the ratio of the current stress level at that point to the critical stress level of the corresponding lithology damage core. The current stress level is taken from the maximum principal stress value of the three-dimensional principal stress corresponding to that point in the comprehensive geological state field. The critical stress level of the damage core is taken from the critical stress threshold for irreversible damage transition of the corresponding lithology rock mass in the potential function space constructed in step 1. The closer the damage core threshold proximity value is to 1, the closer the stress level of the rock mass at that point is to the critical damage state, and the higher the probability of crack initiation. The second parameter is the cumulative amount of water influence factor, which is calculated as the cumulative amount of water influence factor at that point. The first parameter is the cumulative value of the water influence factor vector magnitude corresponding to all time steps from the moment of entering the dynamic boundary of the stress shadow zone to the current calculation time. The water influence factor vector is taken from the calculation result of step 42. The larger the cumulative value of the water influence factor, the longer and deeper the rock mass at this point is softened by groundwater, the greater the weakening of the rock mass's mechanical properties, and the easier it is for cracks to initiate. The third parameter is the distance from the dynamic boundary of the stress shadow zone. It is calculated as the shortest straight-line distance from the three-dimensional spatial coordinates of this point to the dynamic boundary of the stress shadow zone delineated in step 2. The smaller the value of this distance, the closer the point is to the front edge area of ​​the mining stress influence, the greater the space for stress growth as the mining face advances, and the higher the potential risk of crack initiation.

[0053] After point-by-point calculation of the three key parameters, a nonlinear coupling function is used to perform coupled calculations of the three parameters, obtaining the initiation potential of water-conducting fractures at each spatial point. The construction of the nonlinear coupling function follows the basic laws of rock mass damage evolution and water-rock interaction. There is a mutually reinforcing positive feedback effect among the three parameters. The increase in the degree of rock mass damage will intensify the infiltration and softening effect of groundwater. The softening effect of groundwater will further reduce the critical stress threshold of rock mass damage, accelerating the stress-driven fracture initiation process. To quantify the probability of initial initiation of water-conducting fractures in the rock mass at the corresponding spatial location within the stress shadow zone, a source for subsequent fracture propagation path search is constructed. The point criterion requires calculating the initiation potential Pi of water-conducting fractures. This calculation process is based on the parameters determined in the aforementioned steps, including dimensionless weighting coefficients α, β, and γ calibrated from indoor core tests and field measurements in the mine; the damage core threshold proximity Cd, the cumulative water influence factor Cw, and the normalized distance Cl from the dynamic boundary of the stress shadow zone calculated point by point according to the corresponding spatial location; and the nonlinear coupling power coefficient n calibrated from the corresponding lithological damage evolution test. Simultaneously, a fixed mathematical constant, the natural constant e, is used. The calculation process strictly follows the mathematical operation priority and is executed step by step in conjunction with the physical laws of rock mass damage and water-rock coupling. The specific process is as follows: The first step is to calculate the weighted contribution value of each disaster-causing driving factor by performing three independent multiplication operations to complete the weight matching of each influencing factor: First, multiply the stress-damage driving term weight coefficient α with the damage kernel threshold proximity Cd at the corresponding spatial location to obtain the weighted contribution value of stress approaching the critical damage state to fracture initiation; Second, multiply the water-rock softening driving term weight coefficient β with the cumulative influence factor Cw of water at the corresponding spatial location to obtain the weighted contribution value of groundwater softening to fracture initiation; Third, multiply the mining front driving term weight coefficient γ with the distance normalized value Cl at the corresponding spatial location to obtain the weighted contribution value of mining stress growth potential to fracture initiation. This step achieves standardized weight matching of the three types of disaster-causing factors and eliminates the coupling bias caused by differences in the dimensions and magnitudes of different parameters. The second step is to calculate the comprehensive index of linear coupling of multiple disaster-causing factors. The three sets of weighted contribution values ​​obtained in the first step are summed to obtain the comprehensive index after linear superposition of the three types of disaster-causing factors. This step integrates the comprehensive influence of three core factors, namely stress driving, water-rock interaction, and mining front, on fracture initiation, and provides a basis for subsequent nonlinear coupling calculations. The third step is to calculate the nonlinear coupling amplification term of the multi-factor positive feedback effect. The linear coupling comprehensive index obtained in the second step is subjected to power operation according to the calibrated nonlinear coupling power coefficient n to obtain the nonlinear coupling amplification term. This step is used to characterize the positive feedback coupling effect between the three types of disaster-causing factors. That is, the increase in the degree of rock mass damage will exacerbate the groundwater infiltration and softening, and the groundwater softening will further reduce the critical stress of rock mass damage and accelerate the initiation of stress-driven cracks, which is in complete agreement with the nonlinear physical law of rock mass damage evolution. The fourth step is to determine the power term of the exponential function. The negative value of the nonlinear coupling amplification term obtained in the third step is taken to obtain the power term of the exponential function with the natural constant e as the base. This step provides basic parameters that conform to the exponential distribution law for subsequent probabilistic mapping. The fifth step is to calculate the baseline index value of the crack initiation risk probability. Using the natural constant e as the base and the power term obtained in the fourth step as the exponent, the exponentiation operation is performed to obtain the baseline index value of the risk probability. This operation is based on the exponential probability distribution model, which maps the comprehensive index after nonlinear coupling to the numerical range of 0 to 1, ensuring that the calculation results have a unified probabilistic quantitative scale. The sixth step is to summarize and calculate the final water-conducting fracture initiation potential Pi. Using a constant of 1 as the baseline value, the risk probability baseline index value obtained in the fifth step is subtracted to obtain the water-conducting fracture initiation potential Pi at the corresponding spatial location. The final calculation result is in the range of 0 to 1. The larger the value, the higher the probability of initial initiation of water-conducting fractures in the rock mass at that location, providing a clear source criterion for subsequent fracture propagation path search.

[0054] The derivation of this formula is based on a risk probability model of rock mass fracture initiation. The risk of fracture initiation increases non-linearly with the increase of three key parameters. An exponential probability distribution function is used to map the coupled comprehensive index to the range of 0 to 1, ensuring that the numerical value of the initiation potential has a unified quantitative scale. The introduction of the power term is used to characterize the positive feedback coupling effect between the three parameters, which conforms to the non-linear evolution characteristics of rock mass damage and water-rock interaction. Among them, α corresponds to the contribution of mining-induced stress to the rock mass when it approaches the damage criticality. The harder the rock, the stronger the brittleness, and the higher the geostress level, the higher the value of α. β The values ​​of β correspond to the contribution of groundwater to the weakening of rock mass mechanical properties and water pressure-induced fracturing. The stronger the water-bearing capacity of the aquifer and the more significant the softening and disintegration characteristics of the rock mass when exposed to water, the higher the value of β. The values ​​of γ correspond to the contribution of mining advancement to the stress growth potential. The faster the working face advances, the greater the mining height, and the more severe the mining pressure manifestation, the higher the value of γ. In hard rock and highly water-bearing working faces, the values ​​of the three coefficients are α=0.45, β=0.45, and γ=0.10, respectively. In soft rock and shallow buried high-intensity mining working faces, the values ​​of the three coefficients are α=0.25, β=0.45, and γ=0.30, respectively.

[0055] Step 52: Based on the initiation potential of water-conducting fractures, and constrained by the vector distribution of rock mass strength and water influence factors in the constructed comprehensive geological state field, search for the minimum energy dissipation path starting from each initiation potential point. The path superposition probability is used as the evolution potential of water-conducting fractures to form a three-dimensional dynamic cloud map. The specific operations are as follows: Based on the water-conducting fracture initiation potential calculated in step 51, the water-conducting fracture evolution potential is calculated and a three-dimensional dynamic cloud map is constructed. The water-conducting fracture evolution potential is used to characterize the probability of directional propagation of fractures that have already initiated in the rock mass. Its calculation process follows the principle of minimum energy dissipation in rock fracture mechanics, that is, the propagation of fractures in the rock mass always proceeds along the path that requires the least energy to propagate per unit length. First, the spatial point where the water-conducting fracture initiation potential calculated in step 51 exceeds a preset threshold is taken as the starting point of fracture initiation, that is, the source point of fracture propagation. The preset threshold is determined according to the target... The geological conditions and water hazard prevention requirements of the mine have been calibrated. The rock mass strength and water influence factor vector distribution in the comprehensive geological state field constructed in step 4 are used as constraints for fracture propagation path search. The rock mass strength is the shear and tensile strength of the rock mass after water-weakening correction, taken from the water-damaged correction field generated in step 43. The water influence factor vector distribution is taken from the calculation results in step 42. The lower the rock mass strength, the larger the magnitude of the water influence factor vector, the less energy is required for fracture propagation along that location, and the easier it is to become a fracture propagation path. For each fracture initiation source point, all possible fracture propagation paths originating from that source point are searched in three-dimensional space. The total energy dissipation value corresponding to each path is calculated, and the path with the smallest total energy dissipation value is selected as the minimum energy dissipation path corresponding to that source point. The energy dissipation value per unit length at each point on the path is jointly determined by the rock mass strength, the magnitude of the water influence factor vector, and the current stress level at that point. The energy dissipation value per unit length is positively correlated with the rock mass strength and negatively correlated with the magnitude of the water influence factor vector.

[0056] After completing the search for the minimum energy dissipation path corresponding to all fracture initiation points, all paths are superimposed in three-dimensional space. The number of times each spatial calculation point is traversed by the minimum energy dissipation path is calculated point by point. The ratio of this number to the total number of paths is taken as the path superposition probability of that point. This path superposition probability is the water-conducting fracture evolution potential of the corresponding spatial point. The larger the value of the water-conducting fracture evolution potential, the higher the probability that the point will become a water-conducting fracture propagation path, and the greater the risk of water-conducting channel connection. After completing the calculation of the water-conducting fracture evolution potential of all spatial points in the stress shadow zone, based on the evolution potential values ​​of all points in three-dimensional space, a continuous color-gradient mapping method is used to convert the numerical value of the evolution potential into the corresponding color feature, constructing a three-dimensional dynamic cloud map covering the entire impact area in front of the mining. This cloud map can be updated synchronously with the advancement of the mining face and the real-time update of the comprehensive geological state field, intuitively presenting the spatial distribution characteristics and evolution trend of water-conducting fractures.

[0057] In some embodiments of the present invention, step 6 is further included, which involves using the water-conducting fracture evolution potential cloud map to adjust the cutting path and advance speed of the mining equipment. The specific operation is as follows: The system receives the 3D dynamic cloud map of the evolution potential of water-conducting fractures generated in step 5, and performs dynamic optimization and adjustment of the operating parameters of the mining equipment to achieve closed-loop linkage between geological hazard prediction and early warning and intelligent mining control. The preceding steps completed the dynamic inversion of the geological state of the rock mass under the influence of mining and the quantitative prediction of the evolution risk of water-conducting fractures. The generated 3D dynamic cloud map fully characterizes the spatial distribution and risk level of water-conducting fracture development in the mining front area. This quantitative risk result is transformed into directly executable control commands for the mining equipment. Through active correction of the cutting path and dynamic control of the advance speed, the system avoids water inrush disasters caused by the penetration of water-conducting fractures, while ensuring the continuity and efficiency of the working face mining operation. The complete logic from risk vector characterization to control parameter optimization first transforms the scalar evolution potential cloud map into a vector field that can characterize the risk direction, clarifying the risk transmission direction and safe avoidance path. Then, based on the vector field information, the system completes the correction of the cutting trajectory and the matching control of the advance speed. All adjustment actions are synchronously iterated with the real-time update of the water-conducting fracture evolution potential cloud map, ensuring that the mining operation is always within a safe and controllable range.

[0058] Step 6 also includes the following steps: Step 61: Based on the three-dimensional dynamic cloud map of the evolution potential of the formed water-conducting fracture, perform three-dimensional gradient calculation, trace the gradient streamlines to identify the risk approach direction and the safety avoidance direction, and construct a risk tendency vector field covering the area in front of the mining. The specific operations are as follows: First, a global three-dimensional gradient calculation is performed on the three-dimensional scalar field. The gradient calculation is obtained by taking the partial derivatives of the water-conducting fracture evolution potential value at each spatial point along the three orthogonal coordinate axes of the three-dimensional spatial coordinate system. The calculation result for each spatial point is a three-dimensional gradient vector. The magnitude of the gradient vector characterizes the spatial rate of change of the water-conducting fracture evolution potential at that point. The positive direction of the gradient vector points towards the direction of the fastest increase in the evolution potential value, which is defined as the risk approach direction. The negative direction of the gradient vector points towards the direction of the fastest decrease in the evolution potential value, which is defined as the safety avoidance direction. After completing the global gradient vector calculation, starting from all calculation points on the current coal face of the mining face, along each starting point... The gradient vector corresponding to a point is used for streamline tracing in the positive direction. During the tracing process, the gradient vector extends point by point along the direction of the gradient vector in three-dimensional space to form continuous gradient streamlines. Each gradient streamline completely represents the risk transmission path and evolution trend starting from the coal face. Based on the gradient vectors of all spatial points in the entire domain and the gradient streamlines obtained by tracing, a risk trend vector field covering the entire mining area in front of the mining is constructed. The vector direction of each spatial point in the vector field corresponds to the risk approach direction at that location, and the vector magnitude corresponds to the rate of risk growth. This vector field is synchronously iterated with the update of the three-dimensional dynamic cloud map of the water-conducting fracture evolution potential, and always remains consistent with the geological risk evolution state under the current mining process.

[0059] The three-dimensional gradient is calculated as follows: the gradient vector of the hydroconductive fracture evolution potential at a certain spatial point is equal to the partial derivative of the evolution potential value along the x-axis multiplied by the x-axis unit vector, plus the partial derivative of the evolution potential value along the y-axis multiplied by the y-axis unit vector, plus the partial derivative of the evolution potential value along the z-axis multiplied by the z-axis unit vector.

[0060] The calculation logic is abstracted from the basic mathematical definition of gradient in vector analysis. The gradient of a scalar field is a vector whose direction is the direction of the maximum rate of change of the scalar field at that point, and its magnitude is the value of that maximum rate of change. The three-dimensional dynamic cloud map of the evolution potential of water-conducting fractures is a typical three-dimensional spatial scalar field. Through gradient calculation, the scalar form of risk value can be transformed into the vector form of risk direction information, clarifying the specific direction of risk transmission and safety avoidance. All parameters involved in the calculation process have clear physical meanings. Among them, the gradient vector represents the direction and rate of change of the evolution potential of water-conducting fractures at a certain point in three-dimensional space, the evolution potential value represents the risk level of water-conducting fracture development at the corresponding spatial point, x, y, and z represent the three orthogonal coordinate axes of the three-dimensional spatial coordinate system, the unit vectors correspond to the standard directions of the three coordinate axes, and the partial derivatives represent the rate of change of the evolution potential value along the corresponding coordinate axis.

[0061] Step 62: Based on the safe avoidance direction in the risk tendency vector field, the cutting path is corrected. At the same time, the allowable disturbance time window to reach the boundary of the high-risk zone is calculated according to the three-dimensional dynamic cloud map of the evolution potential of the formed water-conducting fracture, and the propulsion speed is adjusted accordingly. The specific operation is as follows: During the correction of the cutting path, the safe avoidance direction corresponding to each point in the area to be cut in the working face is first extracted from the risk tendency vector field. Based on the initial cutting path preset by the mining equipment, the cutting trajectory is corrected segment by segment along the safe avoidance direction. The initial cutting path is preset based on the working face design mining parameters. During the correction process, the trajectory direction is adjusted in two dimensions: cutting height and cutting depth, so that the corrected cutting path always avoids the risk area with a high value of water-conducting fracture evolution potential, while ensuring the continuity and smoothness of the cutting trajectory and meeting the operating conditions requirements of the mining equipment. The corrected cutting path can retain a safe isolation coal-rock pillar that meets the requirements of mine water hazard prevention and control specifications between the boundary of the goaf formed by the cutting operation and the high-risk area. During the control of the advance speed, based on the three-dimensional dynamic cloud map of the evolution potential of water-conducting fractures, the evolution potential threshold corresponding to the high-risk area is pre-calibrated according to the geological conditions of the target mine and the requirements for water hazard prevention and control. Continuous spatial areas with evolution potential values ​​exceeding the threshold are selected as high-risk areas. The boundary of the high-risk area is determined on the side closest to the current mining face. The shortest straight-line distance from the coal wall of the current working face to the boundary of the high-risk area is calculated. Combining the time effect of rock mass damage evolution and the cumulative effect of mining disturbance, the allowable disturbance time window to reach the boundary of the high-risk area is calculated. The allowable disturbance time window is the longest working time that the working face can advance to near the boundary of the high-risk area without triggering the penetration of water-conducting fractures in the high-risk area.

[0062] To quantify the maximum safe operating time of the mining face without triggering the penetration of water-conducting fractures in high-risk areas and ensuring the safety of water hazard prevention and control, and to provide a precise quantitative basis for the dynamic control of the face advance speed, it is necessary to calculate the allowable disturbance time window T. This calculation process is based on four physical quantities determined or pre-calibrated in the aforementioned steps as input parameters: the rock mass damage disturbance coefficient η, calibrated by indoor core mining damage tests and on-site mine pressure measurement data; the shortest spatial straight-line distance D from the current coal face to the boundary of the high-risk area, calculated based on the three-dimensional dynamic cloud map of the water-conducting fracture evolution potential; the safety isolation distance S, pre-set according to the geological conditions of the target mine and the water hazard prevention and control specifications; and the reference advance speed vr, pre-determined in the working face mining design. The specific step-by-step derivation and calculation process strictly follows the mathematical operation priority and is executed in combination with the physical laws of mining-induced rock mass damage evolution, as follows: The first step is to calculate the effective safe distance for the working face to advance. This involves subtracting the preset safety isolation distance S from the shortest straight-line distance D from the current working face coal wall to the boundary of the high-risk zone. The physical meaning of this step is that the safety isolation distance S is the minimum reserved rock pillar distance to ensure that mining disturbances will not be transmitted to the high-risk zone and to avoid triggering the accelerated evolution of water-conducting fractures. By eliminating the safety reserved sections that must be retained through difference calculation, the maximum effective spatial distance that the working face can advance without exceeding the safety limit is obtained, providing a spatial basis for subsequent time calculations. The second step is to calculate the total allowable advance distance after the rock mass damage characteristics are corrected. The effective safe distance calculated in the first step is multiplied by the input rock mass damage disturbance coefficient η to obtain the total allowable advance distance after the rock mass damage characteristics are corrected. The physical meaning of this step is that the rock mass damage disturbance coefficient η is a correction coefficient that matches the damage evolution rate and disturbance effect diffusion characteristics of the corresponding rock mass after mining disturbance. It can eliminate the deviation between the theoretical spatial distance and the actual mining disturbance response of the rock mass on site, ensure that the calculation results are consistent with the real damage evolution law of the underground rock mass, and avoid the theoretical calculation being too conservative or exceeding the safety threshold. The third step involves summarizing and calculating the final allowable disturbance time window T. The corrected allowable total advance distance calculated in the second step is used as the numerator, and the input working face design reference advance speed vr is used as the denominator. A division operation is then performed to obtain the allowable disturbance time window T under the corresponding mining conditions. The physical meaning of this calculation result is the maximum time the working face can operate continuously and safely under the current geological risk state and mining design parameters, without triggering the risk of water inrush through water-conducting fractures in high-risk areas. This provides a direct and implementable quantitative criterion for the dynamic control of the subsequent working face advance speed.

[0063] All parameters involved in the calculation process have clear physical meanings. Among them, the permissible disturbance time window represents the maximum safe operating time without triggering the penetration of water-conducting fractures in the high-risk area; the shortest spatial straight-line distance represents the spatial interval between the current working face and the boundary of the high-risk area; the safe isolation distance represents the pre-set minimum interval distance that meets the requirements for water hazard prevention and control; the reference advance speed represents the standard advance rate designed for the working face; and the rock mass damage disturbance coefficient is determined by the damage evolution characteristics of the corresponding lithology and rock mass and the intensity of mining disturbance, and is calibrated through indoor tests and field measurement data. The working face advance speed is dynamically controlled based on the length of the permissible disturbance time window. When the permissible disturbance time window is short, the advance speed of the mining equipment is reduced to decrease the mining disturbance intensity corresponding to a single cutting cycle, thus delaying the evolution process of rock mass damage and water-conducting fractures inside the surrounding rock. When the permissible disturbance time window is long, the advance speed can be maintained or increased within the range that meets safety requirements to ensure the working face recovery efficiency. The correction of the cutting path and the regulation of the advance speed are executed synchronously. After each complete cutting cycle, the control parameters are iteratively optimized based on the updated three-dimensional dynamic cloud map of the evolution potential of water-conducting fractures and the risk trend vector field, so as to realize the closed-loop dynamic control of mining operations and geological risk prevention and control.

[0064] The above description is merely a preferred embodiment of the present invention; however, the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and its improved concepts, should be covered within the scope of protection of the present invention.

Claims

1. A method for three-dimensional geological modeling of coal mines based on machine learning and multi-source data from directional boreholes, characterized in that: include: A rock mass stress-damage constitutive relation library containing multiple physical state kernels is constructed. Each physical state kernel has a mathematical expression describing the relationship between acoustic emission or microseismic response, resistivity change and stress loading. Microseismic events are acquired, identified as events reaching the damage core, their spatiotemporal distribution is analyzed, and the dynamic boundary of the stress shadow zone is determined by using the constitutive relation library to reverse the dynamic stress field. Based on the dynamic boundary of the stress shadow zone, virtual calculation nodes are deployed inside it, and virtual microseismic activity and virtual resistivity values ​​are deduced according to the constitutive relation library and the stress level of the corresponding points. By comparing the virtual resistivity value with the measured resistivity value, the difference is defined as the influence factor of water, and a corrected comprehensive geological state field is output through physical compensation. In the comprehensive geological state field, based on whether the rock mass has reached the damage core threshold, the cumulative amount of water influence factors, and the distance from the dynamic boundary of the stress shadow zone, the evolution potential of water-conducting fractures is calculated and a three-dimensional dynamic cloud map is formed. The evolution potential cloud map of water-conducting fractures is used to adjust the cutting path and advance speed of mining equipment.

2. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 1, characterized in that, Construct a rock mass stress-damage constitutive relation library containing multiple physical state cores, including: Variable path loading and acoustic-electric synchronous response tests were performed on the rock core. The waveforms of acoustic emission events and complex impedance spectra under different loading paths were recorded. The inherent response characteristics determined by lithology were separated from the coupled signals by the waveform spectrum feature decoupling algorithm. Based on the inherent response characteristics, the stress values, damage degree and resistivity change values ​​measured in the rock core at different loading stages are projected into the potential energy field to construct a potential function space to describe the physical state and evolution path of the rock mass.

3. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 2, characterized in that, Acquire microseismic events, identify those reaching the damage core, and analyze their spatiotemporal distribution, including: Waveform features are extracted from the original microseismic waveforms. The extracted waveform features are compared with the theoretical waveform features of the corresponding lithology in the potential function space during the damage core transition. Events that match the comparison are identified as pure damage core microseismic events. Based on pure damaged nuclear microseismic events, the spatiotemporal sequence and energy triggering relationship between events are analyzed, and a directed energy transfer network is constructed, with the edge of the network serving as the dynamic boundary of the stress shadow zone.

4. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 3, characterized in that, Using the constitutive relation library to reverse the driving stress field, the dynamic boundary of the stress shadow region is determined, including: Based on the energy release nodes in the constructed stress release and transmission network, the failure conditions that satisfy the Mohr-Coulomb criterion at the nodes are taken as the first type of mechanical constraints, and the boundaries of the goaf and the exposed roadway are taken as the second type of mechanical constraints, thus constructing a set of mechanical boundary constraints. Based on the set of mechanical boundary constraints, a rock mass stress-damage constitutive relation library is constructed. With the stress release and transmission network as the skeleton, the stress values ​​of the nodes are recursively calculated along the network, and then a dynamic stress field is constructed by interpolation. Based on the dynamic stress field and rock mass stress-damage constitutive relation library, the ratio of the current stress level to the critical stress level of the damage core in the potential function space is calculated point by point as the stress over / under ratio. The outer boundary of the continuous region where the stress over / under ratio is less than 1 and greater than the preset threshold is tracked to obtain the dynamic boundary of the stress shadow area.

5. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 4, characterized in that, Based on the dynamic boundary of the stress shadow zone, virtual computing nodes are deployed within it, including: Based on the dynamic boundary and dynamic stress field of the stress shadow area, the stress gradient tensor and stress excess / excess ratio of each point inside the stress shadow area are calculated. Regions where the stress gradient amplitude exceeds the threshold are identified and gradient streamline intersection points are extracted as first-class nodes. Points with stress excess / excess ratio greater than the preset first threshold and less than 1 are selected as second-class nodes. The two types of nodes are merged and assigned lithological identifiers and spatial coordinates to form a virtual node layout set. Based on the virtual node deployment set, the corresponding potential function space is retrieved from the rock mass stress-damage constitutive relation library according to the lithological identifier of each node. The initial virtual microseismic activity and virtual resistivity value are derived by taking the real-time stress value at each node as input. Then, spatial consistency correction is performed based on the energy transfer direction and correlation strength between nodes revealed by the stress release and transfer network, and the virtual microseismic activity and virtual resistivity cooperative field is output.

6. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 5, characterized in that, Based on the constitutive relation library and the corresponding point stress levels, virtual microseismic activity and virtual resistivity values ​​are derived, including: For each virtual node, the dynamic stress field time history of the node since entering the dynamic boundary of the stress shadow zone is traced back, the stress loading path is extracted from it, the transition trajectory along the path is simulated in the potential function space, and the endpoint of the trajectory is taken as the real initial damage state of the virtual node. Based on the actual initial damage state and the real-time stress value of the node at the current moment, the preliminary virtual microseismic activity and virtual resistivity value are derived in the potential function space. Then, the actual energy release value of the adjacent nodes in the stress release and transmission network is introduced as a disturbance input for dynamic correction, and the corrected virtual microseismic activity and virtual resistivity value are output. The virtual cumulative released energy characterized by the modified virtual microseismic activity is compared with the measured cumulative released energy of the pure damage core microseismic event at the corresponding spatial location. The energy release deviation is calculated. When the deviation exceeds a preset threshold, the parameters of the potential function space are adjusted using the deviation.

7. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 6, characterized in that, By comparing virtual resistivity values ​​with measured resistivity values, the difference is defined as the influence factor of water, including: The calculated energy release deviation is introduced as a correction coefficient to correct the virtual resistivity value in the virtual microseismic activity and virtual resistivity synergistic field, thus obtaining the reference virtual resistivity field. The reference virtual resistivity field is compared with the measured resistivity field at multiple frequency points. The dispersion effect characteristics and induced polarization characteristics are extracted from the comparison results. The water influence factor vector is defined based on the dispersion effect characteristics and induced polarization characteristics.

8. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 7, characterized in that, The comprehensive geological state field, corrected through physical compensation, includes: Based on the defined water influence factor vector, the constructed dynamic stress field, and the constructed potential function space, the mechanical parameters of the corresponding lithology in the potential function space are dynamically weakened and corrected to generate a water-bearing damage correction field. Based on the water-bearing damage correction field, combined with the constructed dynamic stress field, the output virtual microseismic activity and virtual resistivity synergistic field, a dynamic comprehensive geological state field is constructed through state vector synthesis rules.

9. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 8, characterized in that, Calculate the evolution potential of water-conducting fractures and generate a three-dimensional dynamic cloud map, including: In the comprehensive geological state field, the ratio of the current stress level at each point to the critical stress level of the damage core is calculated point by point as the proximity of the damage core threshold, the cumulative amount of water influence factor, and the distance to the dynamic boundary of the stress shadow area. The initiation potential of water-conducting fractures is obtained through nonlinear coupling. Based on the initiation potential of water-conducting fractures, and constrained by the vector distribution of rock mass strength and water influence factors in the constructed comprehensive geological state field, the minimum energy dissipation path starting from each initiation potential point is searched, and the path superposition probability is used as the evolution potential of water-conducting fractures to form a three-dimensional dynamic cloud map.

10. The coal mine three-dimensional geological modeling method based on machine learning and multi-source data from directional boreholes according to claim 9, characterized in that, The evolution potential cloud map of water-conducting fractures is used to adjust the cutting path and advance speed of mining equipment, including: Based on the three-dimensional dynamic cloud map of the evolution potential of the formed water-conducting fracture, three-dimensional gradient calculation is performed, gradient streamlines are tracked to identify the risk approach direction and the safety avoidance direction, and a risk tendency vector field covering the area in front of the mining is constructed. The cutting path is corrected based on the safe avoidance direction in the risk tendency vector field. At the same time, the allowable disturbance time window to reach the boundary of the high-risk zone is calculated based on the three-dimensional dynamic cloud map of the evolution potential of the formed water-conducting fracture, and the propulsion speed is adjusted accordingly.

Images