Digital-twin-based dynamic evaluation method for geological environment of mine goaf

Abstract

The present application relates to the technical field of mine safety and geological environment evaluation, and particularly relates to a kind of dynamic evaluation method of mine goaf geological environment based on digital twinning.The method comprises the following steps: constructing the digital twin model of mine goaf and designing data structure and virtual-real mapping rules, collecting multi-source heterogeneous data and performing spatio-temporal data indexing processing to form spatio-temporal correlation data stream; extracting microseismic travel time data from the spatio-temporal correlation data stream; performing parallel ray tracing and fast inversion on the microseismic travel time data to obtain the initial value of wave velocity distribution; using the initial value of wave velocity distribution to perform wave velocity field adaptive grid refinement on the microseismic travel time data to obtain a multi-resolution wave velocity field; analyzing the wave velocity variation rate of the multi-resolution wave velocity field to obtain the wave velocity variation degree tensor.The present application turns the abstract energy transfer process into a quantifiable data field for analysis, achieving accurate dynamic evaluation of the geological environment of the goaf.

Description

Technical Field

[0001] This invention relates to the field of mine safety and geological environment assessment technology, and in particular to a dynamic assessment method for the geological environment of mine goaf based on digital twins. Background Technology

[0002] Existing geological environment assessment technologies for mining goaf areas face challenges in integrating multi-source heterogeneous data. Data from different monitoring systems exhibits inconsistent formats and spatiotemporal benchmarks, resulting in severe "data silos." The lack of a unified spatiotemporal correlation mechanism leads to low data utilization efficiency. Traditional data processing methods struggle to achieve collaborative analysis of cross-source data and fail to uncover deep correlations between data. Traditional assessment methods also struggle to transition from static to dynamic perspectives. They rely on discrete time-point measurements, resulting in limited sampling frequency and sparse observation point distribution. Assessment models are based on static equilibrium, failing to characterize the dynamic evolution of the geological environment. Assessment results often lag behind actual geological changes, hindering early warning decision-making. Furthermore, stress redistribution and energy transfer processes are difficult to quantify. Traditional methods lack quantitative descriptions of stress and energy field evolution; energy transfer paths and patterns within rock masses are difficult to identify; and the spatiotemporal relationship between energy accumulation and release is ambiguous, making it impossible to accurately predict the development trend of potential instability areas.

[0003] In summary, existing technologies suffer from problems such as insufficient data fusion capabilities, statically outdated evaluation methods, and difficulty in quantifying energy transfer processes, which urgently need to be addressed. Summary of the Invention

[0004] Therefore, it is necessary to provide a dynamic evaluation method for the geological environment of mining goaf areas based on digital twins to solve at least one of the above-mentioned technical problems.

[0005] To achieve the above objectives, a dynamic evaluation method for the geological environment of mine goaf based on digital twins includes the following steps:

[0006] Step S1: Construct a digital twin model of the mining goaf area and design data structure and virtual-real mapping rules; collect multi-source heterogeneous data and perform spatiotemporal data indexing to form a spatiotemporal correlated data stream;

[0007] Step S2: Extract microseismic travel time data from the spatiotemporal correlated data stream; perform parallel ray tracing and fast inversion on the microseismic travel time data to obtain initial values ​​of wave velocity distribution; use the initial values ​​of wave velocity distribution to perform adaptive mesh refinement of the wave velocity field on the microseismic travel time data to obtain a multi-resolution wave velocity field; perform wave velocity change rate analysis on the multi-resolution wave velocity field to obtain the wave velocity variability tensor.

[0008] Step S3: Based on the wave velocity variability tensor, identify the high-speed and low-speed regions. Combine the microseismic energy release data from the spatiotemporal correlation data stream to calculate the regional energy accumulation to release ratio and obtain the energy balance index field. By tracking the spatiotemporal pattern of energy transfer from the high-speed region to the low-speed region through the energy balance index field, identify regions with high energy accumulation and low release, as well as stress shadow regions, and obtain the energy transfer flux field.

[0009] Step S4: Identify key structural planes and rock layer interfaces based on wave velocity variability tensor; analyze and evaluate the geological structural evolution characteristics of structural planes and rock layer interfaces to obtain a geological structural stability evaluation matrix;

[0010] Step S5: Integrate the wave velocity variability tensor, energy transfer flux field, and geological structure stability evaluation matrix to conduct a comprehensive evaluation of the geological environment and obtain the dynamic level spectrum of the geological environment.

[0011] This invention overcomes the "data silo" problem of multi-source heterogeneous data by constructing accurate 3D entity models and designing twin mapping rules that include a data cascading refresh mechanism. In particular, by constructing a composite spatiotemporal index structure and forming a spatiotemporally correlated data stream, it not only achieves real-time, high-fidelity mapping from the physical world to the digital world, but also ensures efficient retrieval and integration of all relevant data within the spatiotemporal neighborhood when critical events occur, providing a high-quality, strongly correlated data foundation for subsequent dynamic analysis.

[0012] By designing parallel ray tracing and fast inversion algorithms, the reconstruction of wave velocity fields, which traditionally takes several hours, is shortened to minutes, solving the fundamental problems of low computational efficiency and delayed evaluation results in traditional methods. Furthermore, by introducing an octree-based adaptive mesh refinement technique, computational resources are intelligently concentrated on key regions with drastic wave velocity changes without sacrificing overall computational efficiency, significantly improving spatial resolution and inversion accuracy. The resulting wave velocity variability tensor transforms static wave velocity parameters into a dynamic four-dimensional index capable of finely characterizing the stress state and damage evolution of the rock mass.

[0013] This research solves the problem of traditional methods struggling to quantify and assess stress redistribution and energy transfer processes. The energy balance index field can accurately identify high-risk areas such as "stress shadows" where energy is highly accumulated but insufficiently released. The energy transfer flux field visualizes and vectorizes the abstract energy transfer process, clearly revealing the path, direction, and intensity of energy transfer from stress concentration areas to potential instability areas. This provides a new physical dimension and direct quantitative evidence for understanding and predicting the formation mechanisms of dynamic disasters.

[0014] The evaluation process elevates from macroscopic displacement monitoring to a refined analysis of deformation patterns and integrity evolution of key geological structural surfaces. By introducing three-dimensional digital image correlation methods and strain tensor analysis, specific deformation modes such as tension, compression, and shear can be identified. In particular, by constructing a dynamic integrity index field that integrates mechanical deformation and energy convergence effects, a deeper and more accurate assessment of rock mass damage is achieved. The final geological structure stability evaluation matrix provides a comprehensive and quantitative diagnosis of the stability of each key component, making safety management more targeted.

[0015] By establishing a scientific and systematic comprehensive evaluation process, an objective, dynamic, and forward-looking assessment of the geological environment of goaf areas was achieved. A weight optimization method combining the analytic hierarchy process (AHP) and entropy weighting avoided subjective arbitrariness in the evaluation. A Markov state transition matrix was innovatively introduced to determine the overall evolutionary stage of the system, and dynamic correction coefficients were designed accordingly. This ensures that the final evaluation results not only reflect the current state but also include predictions of future development trends. The generated dynamic geological environment hierarchy provides an intuitive, multi-dimensional, and predictive scientific basis for mine safety management.

[0016] Therefore, this method achieves unified management and correlation analysis of multi-source data by constructing a digital twin model and a spatiotemporal index structure; it applies parallel computing and adaptive grid technology to improve the reconstruction speed of the wave velocity field to a near real-time level, enabling dynamic monitoring of the geological environment; and it innovatively introduces the concepts of wave velocity variability tensor and energy transfer flux field to concretize the abstract energy transfer process into a quantifiable and analyzable data field, thereby achieving accurate dynamic evaluation of the geological environment of the goaf and providing a scientific basis for mine safety management. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating the steps of a dynamic evaluation method for the geological environment of a mine goaf based on digital twins. Detailed Implementation

[0018] The objectives, features, and advantages of this invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings.

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

[0020] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0021] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0022] To achieve the above objectives, please refer to Figure 1 This invention provides a method for dynamic evaluation of the geological environment of mine goaf based on digital twins, comprising the following steps:

[0023] Step S1: Construct a digital twin model of the mining goaf area and design data structure and virtual-real mapping rules; collect multi-source heterogeneous data and perform spatiotemporal data indexing to form a spatiotemporal correlated data stream;

[0024] In this embodiment of the invention, a three-dimensional NURBS solid model of the goaf area is constructed based on AutoCAD DXF format drawings and geological data, and sensor coordinates are accurately labeled in the model. Next, an XML-based twin mapping rule set is designed, defining the pairing relationship between physical equipment and the digital model, the data update cycle, and a data cascading refresh mechanism triggered by high-energy events. Then, multi-source monitoring data is collected via TCP / IP protocol, uniformly converted to JSON standard format using a data adapter, and filtered and denoised. A composite spatiotemporal index structure is constructed using GeoHash and B+ trees. Finally, a data processing pipeline updates the model attributes with real-time data according to the mapping rules, and uses a spatiotemporal correlation algorithm to query and bundle multi-source data within a specific spatiotemporal neighborhood, forming a spatiotemporal correlated data stream.

[0025] Step S2: Extract microseismic travel time data from the spatiotemporal correlated data stream; perform parallel ray tracing and fast inversion on the microseismic travel time data to obtain initial values ​​of wave velocity distribution; use the initial values ​​of wave velocity distribution to perform adaptive mesh refinement of the wave velocity field on the microseismic travel time data to obtain a multi-resolution wave velocity field; perform wave velocity change rate analysis on the multi-resolution wave velocity field to obtain the wave velocity variability tensor.

[0026] In this embodiment of the invention, microseismic travel time data is extracted from a spatiotemporally correlated data stream. A region decomposition strategy is employed to grid the computational domain, and a distributed conjugate gradient iterative solver based on MPI is used to perform parallel ray tracing and fast inversion to obtain initial values ​​for wave velocity distribution. Next, key regions are identified based on wave velocity gradients and ray densities, and the grids of these regions are recursively refined using an octree structure. Local wave velocity re-inversion is performed on the refined multi-level grids, and cubic spline interpolation is used to smooth the boundaries of grids at different resolutions, resulting in a multi-resolution wave velocity field. Finally, a time series sequence of the wave velocity field is constructed using a sliding time window, the wave velocity change rate between adjacent windows is calculated, and the absolute wave velocity value, change rate, and spatial gradient are combined into a four-dimensional wave velocity variability tensor.

[0027] Step S3: Based on the wave velocity variability tensor, identify the high-speed and low-speed regions. Combine the microseismic energy release data from the spatiotemporal correlation data stream to calculate the regional energy accumulation to release ratio and obtain the energy balance index field. By tracking the spatiotemporal pattern of energy transfer from the high-speed region to the low-speed region through the energy balance index field, identify regions with high energy accumulation and low release, as well as stress shadow regions, and obtain the energy transfer flux field.

[0028] In this embodiment of the invention, based on ±15% of the background wave velocity in the goaf as a threshold, high-velocity and low-velocity regions are identified in the wave velocity variability tensor, forming a wave velocity anomaly region map. Then, based on rock physics relationships, the elastic strain energy accumulation intensity of the high-velocity region is estimated, and its ratio to the microseismic energy release density is calculated. After logarithmic transformation and normalization, an energy balance index field is obtained. By analyzing the gradient of this index field and microseismic activity, stress shadow zones with high energy accumulation and low activity are identified. Utilizing... The pathfinding algorithm constructs a network of energy transfer channels connecting high-speed and low-speed regions. Finally, the negative gradient of the energy balance exponential field is projected onto these channels to construct an energy transfer flux field that quantifies the direction and intensity of energy flow.

[0029] Step S4: Identify key structural planes and rock layer interfaces based on wave velocity variability tensor; analyze and evaluate the geological structural evolution characteristics of structural planes and rock layer interfaces to obtain a geological structural stability evaluation matrix;

[0030] In this embodiment of the invention, the spatial gradient of the wave velocity variability tensor is calculated using the three-dimensional Sobel operator. Key structural surfaces are identified by combining this with areas of dense microseismic activity, and then segmented using the watershed algorithm. Next, the displacement of the structural surfaces within a continuous time window is tracked using the three-dimensional digital image correlation (3D-DIC) method, forming a structural displacement field sequence. Then, the Green-Lagrange strain tensor is obtained by calculating the spatial derivative of the displacement field, and the tensile and compressive deformation modes of the structural surfaces are analyzed accordingly. Normalized principal strains and energy are weighted and combined with flux to construct a dynamic integrity index field. Finally, an instability risk assessment is performed on the comprehensive integrity index, displacement rate, and acceleration of key load-bearing components, and all quantitative evaluation results are summarized into a geological structure stability evaluation matrix.

[0031] Step S5: Integrate the wave velocity variability tensor, energy transfer flux field, and geological structure stability evaluation matrix to conduct a comprehensive evaluation of the geological environment and obtain the dynamic level spectrum of the geological environment;

[0032] In this embodiment of the invention, data from wave velocity variability tensors, energy transfer flux fields, and geological structure stability evaluation matrices are integrated into a unified three-dimensional grid using spatial interpolation and Z-score normalization to form comprehensive evaluation data. Then, a multi-dimensional evaluation index system encompassing structural stability, energy balance, and stress rationality is constructed, and the weights of each index are determined using a combination of the Analytic Hierarchy Process (AHP) and entropy weighting. The geological environment state of the goaf is partitioned and graded using a weighted scoring method and K-means clustering algorithm. The historical evolution of the state partitions is analyzed by constructing a Markov state transition matrix to determine the overall evolutionary stage of the system. Finally, a development trend correction coefficient is set according to the evolutionary stage to dynamically adjust the static comprehensive score, generating a dynamic geological environment level spectrum.

[0033] Preferably, step S1 includes the following steps:

[0034] Step S11: Obtain and establish a basic structural model of the goaf based on the geological exploration data and geometric structure of the goaf, wherein the key monitoring points of the sensors and the data acquisition area are marked in the basic structural model of the goaf.

[0035] Step S12: Design the twin mapping rule set for the basic structural model of the goaf;

[0036] Step S13: Collect and standardize multi-source data of the goaf area according to the twin mapping rule set to obtain a standardized monitoring dataset;

[0037] Step S14: Construct a spatiotemporal index structure for a standardized monitoring dataset and a basic structural model of the goaf area;

[0038] Step S15: Map the standardized monitoring dataset to the basic structure model of the goaf area according to the twin mapping rule set, and then use the spatiotemporal index structure to identify the spatiotemporal correlation between different data sources to form a spatiotemporal correlated data stream.

[0039] In this embodiment of the invention, AutoCAD DXF format design drawings and borehole columnar geological exploration data, including goaf roadways, working faces, faults, and rock strata interfaces, are acquired. The data is then converted into a three-dimensional solid model of the goaf using a 3D modeling program. The geological structures are represented using non-uniform rational B-spline (NURBS) surfaces, and the roadway structures are represented using scanned volumes. Simultaneously, the precise three-dimensional coordinates (x, y, z) of the microseismic sensors, stress gauges, and displacement gauges deployed on-site are loaded into the model, and each monitoring point is assigned a unique numerical identifier, forming a basic structural model of the goaf that includes the geometric structure and the layout of the monitoring points.

[0040] A twin mapping rule set was designed. This rule set is stored in XML file format and contains a series of mapping rules. Each rule defines a pairing relationship between a physical monitoring device and its corresponding numerical identifier in the goaf foundation structure model. The rules explicitly specify the data update cycle; for example, the stress gauge data update frequency is set to 1Hz, while microseismic data is updated event-driven. The rule set also includes a data cascading mapping mechanism. For example, when the energy detected by a microseismic sensor exceeds a preset threshold E_0, this mechanism is automatically triggered, forcibly refreshing the data of all stress gauges and displacement gauges within a 50-meter radius of that sensor, ensuring data synchronization in relevant areas during high-energy events.

[0041] Data acquisition and standardization are performed based on a twin mapping rule set. The system connects to the mine microseismic monitoring network and fiber optic stress monitoring system via TCP / IP protocol. A data acquisition adapter is developed that listens to the data port in real time, uniformly parsing and converting received raw binary data streams of different formats into JSON format data containing timestamps (accurate to milliseconds), unique device identifiers, 3D coordinates, monitoring values, and units. For the acquired raw data, a medium-range filter is applied to remove noise points, and the 3-sigma criterion is used to identify and remove outliers exceeding the normal fluctuation range, forming a standardized monitoring dataset.

[0042] A spatiotemporal index structure is constructed. For the spatial dimension, the GeoHash encoding algorithm is used to convert the three-dimensional coordinates of monitoring points and model grid cells into one-dimensional string codes, and a hash index is built based on this. For the temporal dimension, a B+ tree index is constructed using the timestamps of standardized monitoring data as keys. The two indexes are combined to form a composite spatiotemporal index structure, where the leaf nodes of the hash index point to the root node of a B+ tree. This B+ tree organizes the data from all historical moments at that spatial location, thus enabling efficient composite queries of data within a specified time range and spatial region.

[0043] Perform real-time data mapping and correlation fusion. Initiate a data processing pipeline that continuously acquires data from a standardized monitoring dataset. Based on the twin mapping rule set, update the monitoring value of each data record to the corresponding numerical identifier attribute in the goaf foundation structure model. Simultaneously, initiate a spatiotemporal correlation algorithm. For example, when a microseismic event data is received, the algorithm uses the event's occurrence time T and source location P as the center, and utilizes a spatiotemporal index structure to query all stress gauge and displacement gauge data within a time window [T-5s, T+5s] and a spherical radius of 50 meters. Package the queried multi-source data with the microseismic event data into a single data object, forming an element of the spatiotemporally correlated data stream, achieving dynamic correlation of data from different sources under event-driven conditions.

[0044] Preferably, step S2, which involves parallel ray tracing and fast inversion of the microseismic travel time data, includes:

[0045] Computational domain gridding and initial model construction were performed on the microseismic travel time data to obtain the initial wave velocity model for each region;

[0046] Parallel ray path tracing and sensitivity matrix construction are performed on the initial wave velocity model of the partition to obtain the local sensitivity matrix;

[0047] By combining microseismic travel time data, zonal initial wave velocity models, and local sensitivity matrices, a global sparse linear equation system is constructed.

[0048] The wave speed correction vector is obtained by solving the global sparse linear equation system using distributed conjugate gradient iteration.

[0049] The initial value of the wave velocity distribution is obtained by updating and judging the wave velocity based on the wave velocity correction vector.

[0050] In this embodiment of the invention, the microseismic travel time data is used for computational domain gridding and initial model construction. The goaf and its affected area are divided into a uniform Cartesian grid of 10m × 10m × 10m. Based on historical rock sample test data, a uniform initial P-wave velocity value, such as 4500 m / s, is assigned to all grid cells. Subsequently, the entire computational domain grid is divided into N subdomains along the X, Y, and Z directions using a domain decomposition method, where N is the number of processor cores in the parallel computing cluster. Each subdomain and its internal initial velocity value constitute a partitioned initial velocity model.

[0051] Parallel ray path tracing and sensitivity matrix construction are performed on the partitioned initial wave velocity model. On each processor core, the shortest path ray tracing algorithm is executed independently on the partitioned initial wave velocity model it is responsible for. This algorithm calculates the ray path from the microseismic source to the receiving sensor. For rays crossing subdomain boundaries, the position and direction information of the ray's exit point are exchanged between adjacent processor cores via a message passing interface (MPI). Simultaneously, the propagation distance of each ray within each grid cell is calculated; this distance is the element A_ij of the sensitivity matrix A, representing the sensitivity of the i-th ray to the slowness change of the j-th grid cell. Finally, each core generates a local sensitivity matrix containing only information from the grid cells within its subdomain.

[0052] By combining microseismic travel time data, zonal initial wave velocity models, and local sensitivity matrices, a global sparse linear equation system is constructed. This equation system takes the form Am = d, where m is the grid slowness (the reciprocal of wave velocity) correction vector to be determined, and d is the travel time residual vector, whose elements are the differences between observed travel times and theoretical travel times calculated based on the current wave velocity model. The global sensitivity matrix A is a logical combination of all local sensitivity matrices; since each ray passes through only a very small number of grids, this matrix is ​​highly sparse.

[0053] A distributed conjugate gradient iterative solution is performed on the global sparse linear equation system. During the iteration process, the core computation of the product of the matrix and vector (A×p) is performed in a distributed manner. Each processor core only needs to store its corresponding local luminosity matrix and a partial component of the vector p. After completing the multiplication locally, all local results are aggregated through a global reduction operation (MPI_Allreduce) to obtain the global product result. This process avoids storing the entire global matrix on a single node.

[0054] Wave velocity updates and convergence checks are performed based on the wave velocity correction vector. After obtaining the slowness correction vector m in each iteration, the slowness s_old of each grid cell is updated: s_new = s_old + m, and converted into a new wave velocity v_new = 1 / s_new. The root mean square (RMS) value of the travel time residuals of all rays is calculated. When the rate of change of the RMS value between two consecutive iterations is less than 1%, the calculation is considered converged, and the iteration terminates. The wave velocity field obtained at this time is the initial value of the wave velocity distribution.

[0055] Preferably, step S2, which uses the initial value of the wave velocity distribution to perform adaptive mesh refinement of the microseismic travel time data, includes:

[0056] Based on the initial value of wave velocity distribution, key areas of microseismic travel time data are comprehensively identified and marked to obtain a key area marking map.

[0057] Based on the key region marking map, octree-driven mesh recursive processing is performed to obtain an octree multi-level mesh structure.

[0058] By using an octree multi-level grid structure to refine the local wave velocity re-inversion of microseismic travel time data, a non-smoothed multi-resolution wave velocity field is obtained.

[0059] By traversing the octree multi-level mesh structure, a multi-level mesh boundary smoothing process is performed on the non-smooth multi-resolution wave velocity field to obtain the multi-resolution wave velocity field.

[0060] In this embodiment of the invention, key regions are comprehensively identified and marked based on the initial wave velocity distribution value of the microseismic travel time data. The spatial gradient of the initial wave velocity distribution value at each grid cell is calculated, and cells with gradient values ​​greater than a preset threshold G_0 (e.g., 50 (m / s) / m) are marked as "high gradient regions". Simultaneously, the ray density of microseismic events within each grid cell is statistically analyzed, and cells with ray densities greater than a preset threshold D_0 (e.g., 10 rays / cell) are marked as "high density regions". Grid cells that simultaneously meet both the "high gradient region" and "high density region" conditions are identified as key regions requiring further refinement, and a key region marking map corresponding to the wave velocity field spatial grid is generated, where key region cells are marked as 1, and non-key region cells are marked as 0.

[0061] An octree-driven recursive meshing process is performed based on the critical region marker map. Starting with a root cube cell encompassing the entire computational domain, it is checked whether this cell intersects with any critical region cell in the critical region marker map. If they intersect, the cube cell is recursively divided into eight equal-sized sub-cube cells. This checking process is repeated for each newly generated sub-cell until the sub-cell size reaches a preset minimum resolution (e.g., 2.5m) or no longer intersects with any critical region cell. This process ultimately generates an octree multilevel mesh structure with a higher mesh density in critical regions and maintains the initial coarse mesh density in non-critical regions.

[0062] Local wave velocity re-inversion of microseismic travel time data in refined regions was performed using an octree multilevel grid structure. The inversion problem was reconstructed only for the refined regions within the octree multilevel grid structure. All small grid cells within the refined regions were treated as new unknowns, and a new, smaller set of linear equations was established using microseismic ray travel time data passing through these regions. The conjugate gradient method was used to solve these equations, obtaining a higher-resolution wave velocity distribution within the refined regions. This result was then combined with the wave velocity values ​​from the unrefined regions to form a non-smoothed multi-resolution wave velocity field with different resolutions in different regions.

[0063] A multi-level grid boundary smoothing process is performed on the non-smooth multi-resolution wave velocity field by traversing an octree multi-level grid structure. All interfaces between coarse and fine grids in the octree structure are identified. At these interfaces, a cubic spline interpolation function is applied to the wave velocity values ​​on the fine grid side, forcing their boundary values ​​and their first derivatives to be continuous with the wave velocity values ​​on the coarse grid side. This boundary smoothing process eliminates pseudo-discontinuities in wave velocity caused by abrupt grid changes between different resolution regions, thus obtaining a final, spatially continuous multi-resolution wave velocity field.

[0064] Preferably, step S2, which involves analyzing the rate of change of wave velocity in the multi-resolution wave velocity field, includes:

[0065] A time window sequence of wave velocity fields is constructed from the multi-resolution wave velocity fields to obtain the wave velocity field time series sequence;

[0066] Wave velocity rate of change is calculated and anomalies are identified from the time series of wave velocity fields to obtain the wave velocity rate of change field.

[0067] The wave velocity variability tensor is constructed and mapped from the wave velocity rate of change field and the multi-resolution wave velocity field to obtain the wave velocity variability tensor.

[0068] In this embodiment of the invention, a time-window wave velocity field sequence is constructed for the multi-resolution wave velocity field. A time window of 24 hours and a sliding strategy with a step size of 1 hour are set. The multi-resolution wave velocity field calculated within each time window is taken as a data slice, and the time at the midpoint of the window is used as its timestamp. All these wave velocity field data slices with timestamps are stored in a time-indexed data structure in chronological order to form a continuous wave velocity field time series sequence, which records the complete process of the wave velocity field evolution over time.

[0069] Wave velocity change rate calculation and anomaly identification are performed on the time series of wave velocity fields. For each grid cell in space, the wave velocity values ​​V(t_1) and V(t_2) at two adjacent time windows (e.g., t_1 and t_2) are extracted. The wave velocity change rate of this cell is calculated using the formula (V(t_2)-V(t_1)) / (t_2-t_1). After calculation for all cells in the entire space, an instantaneous wave velocity change rate field is formed. Subsequently, the mean μ and standard deviation σ of this change rate field in the local spatial neighborhood are calculated, and regions with absolute change rate values ​​exceeding μ+3σ are identified as wave velocity change anomaly areas.

[0070] A wave velocity variability tensor was constructed and mapped for the wave velocity rate of change field and the multi-resolution wave velocity field. A four-dimensional (x, y, z, t) tensor data structure was constructed, where the tensor element at any spatiotemporal point (x, y, z, t) is a composite data volume. This data volume contains three core components: 1) the absolute wave velocity value V at that point; 2) the scalar wave velocity rate of change dV / dt at that point; and 3) the spatial gradient vector ▽(dV / dt) of the wave velocity rate of change field at that point. The direction of this gradient vector indicates the spatial direction of the most drastic wave velocity change. The resulting wave velocity variability tensor, with the magnitude and directionality of its components, is used to directly reflect the concentration of stress in the rock mass and the dynamic characteristics of damage evolution.

[0071] Preferably, step S3 includes the following steps:

[0072] Step S31: Based on the wave velocity variability tensor and the preset wave velocity anomaly threshold, identify the high-speed region above the threshold and the low-speed region below the threshold to form a wave velocity anomaly region map.

[0073] Step S32: Extract the energy information and spatial location of microseismic events from the spatiotemporal correlation data stream, calculate the cumulative energy released by microseismic events per unit volume, and obtain the energy release density field;

[0074] Step S33: Calculate the energy accumulation and release ratio based on the wave velocity anomaly region map and energy release density field to obtain the energy balance index field;

[0075] Step S34: Analyze the gradient distribution characteristics in the energy balance index field, determine the potential stress shadow area, calculate the geometric characteristic parameters and energy accumulation rate of the stress shadow area, and form a stress shadow area characteristic map;

[0076] Step S35: Based on the spatial distribution of the wave velocity anomaly region map, determine the energy transfer direction and construct an energy transfer channel network;

[0077] Step S36: Construct a flux vector field based on the energy transfer channel network, establish a developmental correlation analysis between the flux vector field and the stress shadow area characteristic spectrum, and obtain the energy transfer flux field.

[0078] In this embodiment of the invention, the background reference value of the wave velocity in the goaf is first calculated. This value is the arithmetic mean of the wave velocities across the entire evaluation area. The wave velocity anomaly threshold is set as follows: ±15%. Traverse the wave velocity components in the wave velocity variability tensor, and assign wave velocities higher than... The grid cells are marked as high-speed regions, and the wave velocity is lower than that of the high-speed regions. The grid cells are marked as low-velocity regions. A 3×3×3 morphological opening kernel is applied to the marked 3D grid data to remove isolated anomalous cells. Finally, a 3D wave velocity anomaly map is generated, classifying and marking high-velocity and low-velocity regions.

[0079] Filter all microseismic event records with timestamps from the spatiotemporal correlated data stream within the past 24 hours. Extract the source energy of each event. and three-dimensional coordinates The entire evaluation space is divided into a three-dimensional grid consistent with the wave velocity field. For each grid cell... The energy of all microseismic events falling within this unit is accumulated, i.e. Divide the total energy by the volume of the grid cells. The energy release density of the unit was obtained. The energy release density values ​​of all grid cells together constitute the energy release density field.

[0080] Based on rock physics relationships, using formulas Estimate the accumulated elastic strain energy intensity within each grid cell in the high-speed region. ,in This represents the current wave velocity of the unit. This is a constant determined by the rock density and elastic modulus. For each grid cell, its energy accumulation intensity is calculated. With energy release density ratio ,in It is a very small positive number (e.g.) To avoid the denominator being zero, the ratio of all grid cells is calculated. By performing logarithmic transformation and min-max normalization, it is mapped to the interval [0,1] to form a dimensionless energy balance exponential field.

[0081] The spatial gradient of the energy balance index field is calculated. Connected regions with an index value greater than 0.9 and local microseismic frequencies below a preset background value of 10% are identified; these regions are designated as potential stress shadow areas. A three-dimensional connected component labeling algorithm is used to segment these regions. For each segmented stress shadow area, its volume, geometric center coordinates, and the temporal rate of change of the average energy balance index over the past three time windows are calculated; this rate of change represents the energy accumulation rate. These parameters collectively constitute a stress shadow area feature map.

[0082] The geometric center of each high-speed and low-speed region in the wave velocity anomaly region map is considered a node in the graph. The energy transfer direction is defined as from the high-speed region to the low-speed region. For each pair of high-speed-low-speed nodes, an operation is performed on the 3D mesh. The pathfinding algorithm finds the optimal path connecting two points. The path cost function of the algorithm is set to the reciprocal of the wave velocity value of the grid cell, that is, it prioritizes passing through regions with lower wave velocities. The set of all found optimal paths is constructed into a directed acyclic graph, i.e., an energy transfer channel network.

[0083] The negative gradient of the energy balance index field is calculated to obtain the energy potential gradient field, which points in the direction of the fastest decrease in the energy index. Then, within each grid cell covered by the energy transfer channel network, the energy potential gradient vector is projected onto the local tangent direction of the channel containing that cell to obtain the direction and magnitude of the flux vector. This vector field is the instantaneous energy flux vector field. By analyzing the changes in the flux field within a continuous time window, especially the spatiotemporal correspondence between the positional evolution of the source and sink points in the flux field and the expansion or contraction of the region in the stress shadow region feature map, the energy transfer flux field is finally established and output.

[0084] Preferably, step S33 includes:

[0085] Wave velocity data from the wave velocity anomaly region map is extracted, and the elastic strain energy accumulation intensity is estimated to obtain the elastic strain energy accumulation field.

[0086] The energy components of the elastic strain energy accumulation field and the energy release density field are spatially registered to obtain the registered energy component data pairs.

[0087] The point-by-point energy accumulation and release ratio of the registered energy component data pairs is calculated to obtain the original energy ratio field.

[0088] The original energy ratio field is constructed by energy balance index normalization to obtain the energy balance index field.

[0089] In this embodiment of the invention, wave velocity data is extracted from the wave velocity anomaly region map to estimate the accumulated elastic strain energy. For each grid cell marked as a high-speed region in the wave velocity anomaly region map, its... Wave velocity value Based on rock physics theory, through formulas Estimate the accumulated elastic strain energy of this element. ,in The density of the rock mass is taken as 2700 kg / m³. The background wave velocity is a pre-calculated reference value. For mesh elements in the non-high-speed region, the accumulated elastic strain energy intensity... Set to 0. This sets all grid cells to 0. The values ​​are combined to form a three-dimensional elastic strain energy accumulation field corresponding to the original mesh.

[0090] Spatial registration of the energy components of the elastic strain energy accumulation field and the energy release density field is performed. It is confirmed that both the elastic strain energy accumulation field and the energy release density field use the same Cartesian coordinate system and mesh partitioning structure. All mesh elements are traversed, and for each spatial coordinate... The mesh elements are extracted, and their accumulated elastic strain energy is also extracted. and energy release density Treat these two values ​​as a data pair. and the coordinates of the cell The data are correlated to form a set of registered energy component data pairs that contain energy component information of all grid cells.

[0091] The point-by-point energy accumulation-release ratio is calculated for the registered energy component data pairs. For each grid cell's registered energy component data pair... Calculate the ratio of energy accumulation to energy release. The calculation formula is: ,in It is a value A constant is used to prevent the denominator from being zero. This calculation is performed on all mesh cells to generate a new three-dimensional data field, namely the original energy ratio field, where the value of each cell is the calculated ratio R.

[0092] The original energy ratio field is constructed using energy balance index normalization. First, the natural logarithm of all values ​​in the original energy ratio field R is taken, yielding... Adding 1 ensures that the parameter of the logarithmic function is positive. Then, find... Maximum value in the field and minimum value Apply the minimum-maximum normalization formula. , will each unit The values ​​are linearly mapped to the interval [0,1] to obtain the final energy balance index. All grid cells These values ​​together constitute the energy balance index field.

[0093] Preferably, step S36 includes:

[0094] Calculate the energy potential gradient field of the energy balance exponential field;

[0095] Construct a constrained energy flow field based on the energy potential gradient field and the energy transfer channel network;

[0096] The instantaneous energy flux vector field is obtained by synthesizing the constrained energy flow direction field and the energy potential gradient field.

[0097] A flux dynamic evolution characteristic map is constructed based on the instantaneous energy flux vector field;

[0098] Based on the flux dynamic evolution characteristic map, energy transfer correlation analysis is performed on the instantaneous energy flux vector field and stress shadow region characteristic map to obtain the energy transfer flux field.

[0099] In this embodiment of the invention, the energy potential gradient field of the energy balance exponential field is calculated. For a three-dimensional energy balance exponential field... The spatial gradient at the center of each grid cell is calculated using the central difference method. Specifically, in The gradient component in the direction is ,in This is the grid spacing. Calculate similarly. and Components of direction and These three components are combined into a gradient vector. The negative gradient vector Defined as an energy potential gradient field, the direction of each vector in this field indicates the direction of the fastest exponential decrease in energy, i.e. the potential driving direction of energy transfer.

[0100] A constrained energy flow field is constructed based on the energy potential gradient field and the energy transfer channel network. Each grid cell in the energy transfer channel network is traversed. For each cell, the local tangent vector of its corresponding channel is obtained. Then, the energy potential gradient vector of this unit is... Projected onto the tangent vector Above, a projection vector consistent with the channel direction is obtained. The projection vector The magnitude of represents the driving force component of energy transfer along the channel, and its direction is constrained along the channel path. The projection vector of all channel elements. The set constitutes a constrained energy flow field.

[0101] The constrained energy flow field and the energy potential gradient field are synthesized to obtain the instantaneous energy flux vector field. For each vector in the constrained energy flow field... Its size is multiplied by the energy balance index value of that grid cell. The flux intensity at that point is obtained, i.e. .vector This is the instantaneous energy flux vector at that point, with its direction indicating the direction of energy flow at that point and its magnitude indicating the intensity of the energy flow. The F vectors of all grid cells together constitute the instantaneous energy flux vector field.

[0102] Construct a dynamic evolution characteristic map of flux based on the instantaneous energy flux vector field. Calculate the divergence of the instantaneous energy flux vector field. Regions with positive divergence are defined as energy source regions, indicating that energy flows out of these regions; regions with negative divergence are defined as energy sink regions, indicating that energy converges into these regions. By analyzing the flux field divergence maps over three consecutive time windows, the spatial location, shape, and intensity changes of the source and sink regions are tracked, forming a flux dynamic evolution characteristic map describing the evolution of the flux field structure.

[0103] Energy transfer correlation analysis was performed on the instantaneous energy flux vector field and stress shadow region feature map based on the flux dynamic evolution characteristic map. The energy sink region in the flux dynamic evolution characteristic map was spatiotemporally superimposed and compared with the stress shadow region in the stress shadow region feature map. If the spatial location of a stress shadow region highly overlaps with a continuously increasing energy sink region, then the energy transfer process is confirmed as the main cause of energy accumulation in the stress shadow region. This instantaneous energy flux vector field, confirmed by the correlation analysis and possessing clear physical meaning, is defined as the final energy transfer flux field.

[0104] Preferably, step S4 includes the following steps:

[0105] Step S41: Calculate the gradient distribution of the wave velocity variability tensor in space, and combine it with the spatiotemporal correlation data stream to identify and segment key structural surfaces, thereby obtaining a structural surface feature library;

[0106] Step S42: Establish time-series tracking identifiers based on the structural surface feature library, and calculate the structural surface displacement by combining the time series of wave velocity variability tensor to obtain the structural displacement field sequence;

[0107] Step S43: Calculate the strain tensor distribution of the structural surface based on the structural displacement field sequence, perform rock mass deformation mode analysis, and obtain the deformation mode spectrum.

[0108] Step S44: Identify high strain energy regions based on deformation mode spectrum and energy transfer flux field, perform dynamic assessment of rock mass integrity, and obtain the integrity dynamic index field;

[0109] Step S45: Conduct an instability risk assessment of key parts of the goaf based on the integrity dynamic index field to obtain an instability risk assessment map;

[0110] Step S46: Conduct a comprehensive evaluation of the stability of the geological structure based on the instability risk assessment map to obtain the geological structure stability evaluation matrix.

[0111] In this embodiment of the invention, a three-dimensional Sobel operator is applied to the wave velocity component in the wave velocity variability tensor to calculate its spatial gradient magnitude in each grid cell. Regions with gradient magnitudes greater than a preset threshold are identified as potential structural surfaces. Simultaneously, the spatial distribution of microseismic events is extracted from the spatiotemporal correlated data stream. Regions with densely distributed microseismic events are spatially superimposed with high-gradient regions, and regions with high overlap are identified as key structural surfaces. A three-dimensional watershed segmentation algorithm is used to segment the identified regions, extracting the independent geometric morphology of each structural surface. Each structural surface is assigned a unique identifier, and its geometry is stored in a triangular mesh model format, forming a structural surface feature library.

[0112] Establish a time-series tracking identifier for each structural surface in the structural surface feature library. For Time and The wave velocity variability tensor at time 1 is used to repeat the identification process in step S41 to obtain the structural surfaces at two time points. Using the three-dimensional digital image correlation (3D-DIC) method, the geometric feature points on the structural surfaces at the two time points are matched to calculate the position of each feature point. The three-dimensional displacement vector over time. Interpolation is performed on the displacement vectors of all feature points on the structural surface to construct a continuous displacement field covering the entire structural surface. The displacement fields calculated at each time step are organized in chronological order to form a structural displacement field sequence.

[0113] The strain tensor distribution of the structural surface is calculated based on the structural displacement field sequence. For the displacement field at each time step, a Green-Lagrange strain tensor is constructed by calculating the first-order partial derivative of the displacement vector with respect to spatial coordinates. Eigenvalue decomposition is performed on this tensor to obtain the three principal strains and their directions. Based on the sign and magnitude of the principal strains, the structural surface is divided into a tensile deformation region (positive principal strain), a compressive deformation region (negative principal strain), and a shear deformation region (opposite signs of principal strains). This partitioning information is then superimposed on the structural surface model using color coding to generate a deformation mode spectrum.

[0114] Spatially register the deformation mode spectrum with the energy transfer flux field. Identify the regions with the top 10% absolute principal strain values ​​in the deformation mode spectrum (high strain region) and the regions with the top 10% energy convergence intensity (negative divergence) in the energy transfer flux field (high energy sink region). Design the rock mass integrity index. ,in For the normalized maximum principal strain, For normalized energy inflow flux, and Weighting coefficients (e.g.) , ). Calculate the integrity index of each structural surface element to form a dynamic integrity index field.

[0115] Based on the goaf design drawings, key bearing coal pillars and roof structures are identified. The integrity dynamic index field, displacement rate, and displacement acceleration of these structures are extracted. An instability risk index is then established. ,in For displacement rate, For displacement acceleration, This is an empirical coefficient. Based on the risk index. The values ​​are set to three levels of risk thresholds, classifying the risk status of key components into three levels: "stable," "critical," and "high-risk." The classification results are visualized to form an instability risk assessment map.

[0116] The evaluation results of the instability risk assessment map are integrated into a geological structure stability evaluation matrix. Each row of this matrix represents a key structural component (e.g., "Coal Pillar No. 1"), and each column represents an evaluation index. Column items include: structural component number, current integrity index, average displacement rate, instability risk level, and an evolution trend prediction ("stable," "deteriorating," or "improving") calculated based on the rate of change of data over the past 12 hours. This matrix provides a comprehensive evaluation conclusion on the overall geological structural stability of the goaf using a combination of quantitative and qualitative methods.

[0117] Preferably, step S5 includes the following steps:

[0118] Step S51: Integrate and standardize multi-source evaluation data for the wave velocity variability tensor, energy transfer flux field, and geological structure stability evaluation matrix to obtain comprehensive evaluation data;

[0119] Step S52: Construct multi-dimensional evaluation indicators and optimize their weights for the comprehensive evaluation data to obtain the evaluation indicator weight matrix;

[0120] Step S53: Based on the evaluation index weight matrix, the comprehensive evaluation data is divided into geological environment state zones and levels to obtain an environmental state zoning map;

[0121] Step S54: Determine the environmental evolution stage and extract features from the comprehensive evaluation data according to the environmental state partition map, to obtain an environmental evolution feature spectrum;

[0122] Step S55: Identify key influencing factors and analyze their contribution degrees according to the environmental evolution feature spectrum, to obtain an influencing factor contribution spectrum;

[0123] Step S56: Conduct a comprehensive assessment of the dynamic grade of the geological environment according to the influencing factor contribution spectrum, to obtain a dynamic grade spectrum of the geological environment.

[0124] In the embodiment of the present invention, multi-source evaluation data integration and standardization are carried out. The data of wave velocity variability tensor and energy transfer flux field are mapped onto a unified three-dimensional grid with a resolution of 5m×5m×5m through spatial interpolation method. For the data in the geological structure stability evaluation matrix, the three-dimensional grid cells occupied by the associated key structural parts are assigned the corresponding index values in the matrix. The Z-score standardization method is applied to the multi-source data (such as wave velocity change rate, energy flux magnitude, integrity index, etc.) in all grid cells, that is, each index is subtracted by its global mean and then divided by its global standard deviation to eliminate the influence of dimension, and comprehensive evaluation data is formed.

[0125] Multi-dimensional evaluation index construction and weight optimization are carried out. Three first-level indexes are constructed from the comprehensive evaluation data: structural stability, energy balance, and stress rationality. Each first-level index has second-level indexes under it. For example, the structural stability has the normalized integrity index and displacement rate under it. The analytic hierarchy process (AHP) is used to construct a judgment matrix, and 5 domain experts score the relative importance of the indexes pairwise, and the initial subjective weight is calculated. Then the entropy weight method is used to calculate the objective weight according to the data dispersion degree of each index in the entire evaluation area. The final weight is calculated by the formula to obtain an evaluation index weight matrix.

[0126] The geological environment state is partitioned and graded according to the evaluation index weight matrix for the comprehensive evaluation data. For each grid cell, the standardized index values of each item are multiplied by the corresponding weights in the evaluation index weight matrix and summed to obtain a comprehensive evaluation score S. Five levels of score thresholds are set: S>0.8 is safe and stable, 0.6<S≤0.8 is basically stable, 0.4<S≤0.6 is critically stable, 0.2<S≤0.4 is unstable, and S≤0.2 is severely unstable. The K-means clustering algorithm is used, with the index vectors of each grid cell as features, to cluster the units with similar geological environment states into different state regions, and the corresponding grades are assigned according to the average score of the units in the region, and an environmental state partition map is generated.

[0127] Environmental evolution stages were determined and features extracted. An environmental state partitioning map was analyzed across 24 consecutive time points, and a Markov state transition matrix was constructed. This matrix records the probability of a grid cell transitioning from one state level to another. Evolutionary stages were defined based on the transition probability characteristics: if the probability of transitioning to an unstable state is significantly higher than the probability of transitioning to a stable state, it is considered a period of dynamic instability; if the probabilities of bidirectional transitions are roughly equal, it is considered an adjustment transition period; if the vast majority of cells remain at a stable or basically stable state, it is considered a period of stable adaptation. The dominant transition paths and average residence time for each stage were extracted to form an environmental evolution feature spectrum.

[0128] Key influencing factors were identified and their contributions analyzed. For each region classified as "unstable" or "severely unstable" in the environmental status zoning map, the weighted score (indicator value × weight) of each secondary evaluation indicator within it was calculated. The indicator with the highest weighted score was identified as the dominant influencing factor for the current state of that region. The contribution of this indicator was defined as the percentage of its weighted score to the sum of the weighted scores of all indicators in that region. The dominant factors and their contributions for all unstable regions were summarized to form an influencing factor contribution spectrum.

[0129] A comprehensive assessment of the dynamic level of the geological environment is conducted. Based on the static comprehensive evaluation score S calculated in step S53, a development trend correction coefficient determined by the environmental evolution characteristic spectrum is introduced. If a region is in a period of "dynamic instability," then If it is in the "adjustment and transition period", then If it is in a "stable adaptation period", then Calculate the final dynamic rating score. This dynamic rating score integrates the current state and future evolution trend, and is combined with the ratings of three primary indicators—structural stability, energy balance, and stress rationality—to form the final dynamic geological environment rating spectrum.

[0130] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0131] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A digital twin-based dynamic evaluation method for a geological environment of a mine goaf, characterized in that, The method comprises the following steps: Step S1: Constructing a digital twin model of a mine goaf and designing data structure and virtual-real mapping rules, collecting multi-source heterogeneous data and performing spatio-temporal data indexing processing to form a spatio-temporal correlation data stream; Step S2: Performing calculation domain gridding and initial model construction on the microseismic travel time data to obtain a partitioned initial wave velocity model; Performing parallel ray path tracking and sensitivity matrix construction on the partitioned initial wave velocity model to obtain a local sensitivity matrix; combining the microseismic travel time data, the partitioned initial wave velocity model and the local sensitivity matrix to construct a global sparse linear equation set; Performing distributed conjugate gradient iteration solving on the global sparse linear equation set to obtain a wave velocity correction vector; performing wave velocity updating and convergence judgment according to the wave velocity correction vector to obtain a wave velocity distribution initial value; Performing key area comprehensive identification and marking on the microseismic travel time data according to the wave velocity distribution initial value to obtain a key area marking map; performing octree-driven grid recursive processing according to the key area marking map to obtain an octree multi-level grid structure; performing refined area local wave velocity re-inversion on the microseismic travel time data by using the octree multi-level grid structure to obtain a non-smooth multi-resolution wave velocity field; performing multi-level grid boundary smooth transition processing on the non-smooth multi-resolution wave velocity field by traversing the octree multi-level grid structure to obtain a multi-resolution wave velocity field; Performing time window wave velocity field sequence construction on the multi-resolution wave velocity field to obtain a wave velocity field time sequence; performing wave velocity change rate calculation and anomaly identification on the wave velocity field time sequence to obtain a wave velocity change rate field; Performing wave velocity variation degree tensor construction and mapping on the wave velocity change rate field and the multi-resolution wave velocity field to obtain a wave velocity variation degree tensor; Step S3 comprises the following steps: Step S31: Identifying high-speed areas higher than a threshold value and low-speed areas lower than the threshold value based on the wave velocity variation degree tensor and a preset wave velocity anomaly threshold value to form a wave velocity anomaly area map; Step S32: Extracting energy information and spatial positions of microseismic events from the spatio-temporal correlation data stream, calculating cumulative energy release of microseismic events per unit volume to obtain an energy release density field; Step S33: Performing energy accumulation and release ratio calculation according to the wave velocity anomaly area map and the energy release density field to obtain an energy balance index field; Step S34: Analyzing gradient distribution characteristics in the energy balance index field to determine a potential stress shadow area, calculating geometric characteristic parameters and energy accumulation rate of the stress shadow area to form a stress shadow area characteristic map; Step S35: Determining an energy transfer direction based on the spatial distribution of the wave velocity anomaly area map, constructing an energy transfer channel network; Step S36: Constructing a flux vector field according to the energy transfer channel network, establishing development correlation analysis between the flux vector field and the stress shadow area characteristic map to obtain an energy transfer flux field; Step S4: Identifying key structural planes and rock stratum interfaces based on the wave velocity variation degree tensor; analyzing and evaluating geological structure evolution characteristics of the structural planes and the rock stratum interfaces to obtain a geological structure stability evaluation matrix; Step S5: Integrating the wave velocity variation degree tensor, the energy transfer flux field and the geological structure stability evaluation matrix to perform comprehensive geological environment evaluation to obtain a geological environment dynamic grade spectrum.

2. The digital-twin-based dynamic evaluation method of the geological environment of a mine goaf according to claim 1, characterized in that, Step S1 comprises the following steps: Step S11: Obtain and establish a goaf basic structure model according to goaf geological exploration data and geometric structure, wherein the goaf basic structure model is marked with key monitoring points of sensors and data collection areas; Step S12: Design a twin mapping rule set of the goaf basic structure model; Step S13: Collect and standardize multi-source data of the goaf according to the twin mapping rule set, to obtain a standardized monitoring data set; Step S14: Construct a space-time index structure of the standardized monitoring data set and the goaf basic structure model; Step S15: Map the standardized monitoring data set to the goaf basic structure model according to the twin mapping rule set, and identify the space-time correlation between different data sources by using the space-time index structure, to form a space-time correlation data stream.

3. The digital-twin-based dynamic evaluation method of the geological environment of a mine goaf according to claim 1, characterized in that, Step S33 includes: Extract the wave velocity data of the wave velocity anomaly area atlas, perform elastic strain energy accumulation strength estimation, and obtain an elastic strain energy accumulation field; Perform energy component spatial registration on the elastic strain energy accumulation field and the energy release density field, to obtain a registered energy component data pair; Perform point-by-point energy accumulation and release ratio calculation on the registered energy component data pair, to obtain an original energy ratio field; Perform energy balance index normalization construction on the original energy ratio field, to obtain an energy balance index field.

4. The digital-twin-based dynamic evaluation method of the geological environment of a mine goaf according to claim 1, characterized in that, Step S36 includes: Calculate an energy potential gradient field of the energy balance index field; Construct a constraint energy flow direction field according to the energy potential gradient field and the energy transfer channel network; Synthesize the constraint energy flow direction field and the energy potential gradient field, to obtain an instantaneous energy flux vector field; Construct a flux dynamic evolution feature map according to the instantaneous energy flux vector field; Perform energy transfer correlation analysis on the instantaneous energy flux vector field and the stress shadow feature map according to the flux dynamic evolution feature map, to obtain an energy transfer flux field.

5. The digital-twin-based dynamic evaluation method of the geological environment of a mine goaf according to claim 1, characterized in that, Step S4 includes the following steps: Step S41: Calculate the gradient distribution of the wave velocity variation degree tensor in space, identify and segment key structural surfaces in combination with the space-time correlation data stream, to obtain a structural surface feature library; Step S42: Establish a time sequence tracking identifier based on the structural surface feature library, and calculate the displacement of the structural surface in combination with the time sequence of the wave velocity variation degree tensor, to obtain a structural displacement field sequence; Step S43: Calculate the strain tensor distribution of the structural surface based on the structural displacement field sequence, perform rock mass deformation mode analysis, and obtain a deformation mode spectrum; Step S44: Identify a high strain energy area based on the deformation mode spectrum and the energy transfer flux field, perform dynamic evaluation of rock mass integrity, and obtain an integrity dynamic index field; Step S45: Perform key part instability risk evaluation of the goaf according to the integrity dynamic index field, to obtain an instability risk evaluation map; Step S46: Perform comprehensive evaluation of the stability of the geological structure according to the instability risk evaluation map, to obtain a geological structure stability evaluation matrix.

6. The digital-twin-based dynamic evaluation method of the geological environment of a mine goaf according to claim 1, characterized in that, Step S5 includes the following steps: Step S51: Integrate and standardize multi-source evaluation data of the wave velocity variation degree tensor, the energy transfer flux field and the geological structure stability evaluation matrix, to obtain comprehensive evaluation data; Step S52: Construct multi-dimensional evaluation indexes and optimize weights, to obtain an evaluation index weight matrix; Step S53: According to the evaluation index weight matrix, the geological environment state partition and classification are performed on the comprehensive evaluation data to obtain an environment state partition map; Step S54: According to the environment state partition map, the environment evolution stage determination and feature extraction are performed on the comprehensive evaluation data to obtain an environment evolution feature spectrum; Step S55: According to the environment evolution feature spectrum, the key influence factor identification and contribution analysis are performed to obtain an influence factor contribution spectrum; Step S56: According to the influence factor contribution spectrum, the geological environment dynamic grade comprehensive evaluation is performed to obtain a geological environment dynamic grade spectrum.

Images