A computer simulation-based design optimization method for multi-level grouting in overburden.

Abstract

This invention relates to the field of computer simulation technology, and discloses a computer simulation-based method for optimizing multi-level grouting design in overburden. The method includes: constructing a three-dimensional geomechanical model of a goaf based on geological exploration data; calculating the development characteristics of each layer in the overburden strata of the goaf, and constructing the grouting-bearing layer based on these characteristics; analyzing fracture development characteristics and layer stress distribution; analyzing grouting simulation parameters for the grouting-bearing layer based on these characteristics; generating a multi-level layered grouting simulation scheme for the grouting-bearing layer; generating a digital trajectory of the grouting-bearing layer based on the multi-level layered grouting simulation scheme, and generating grouting control simulation commands according to the digital trajectory; performing multi-level grouting simulation operations on the grouting-bearing layer according to the grouting control simulation commands, and outputting the multi-level grouting simulation results. This invention can improve the accuracy of grouting levels in multi-level overburden design optimization.

Description

Technical Field

[0001] This invention relates to the field of computer simulation technology, and in particular to a computer simulation-based method for optimizing the design of multi-level grouting in overburden. Background Technology

[0002] Underground coal mining creates goaf areas, damaging the overlying strata and eventually affecting the surface, causing surface movement and deformation. This results in prolonged surface subsidence and makes the damaged land difficult to reuse. Through comprehensive analysis of the overlying strata damage in old goaf areas, and by reconstructing key bearing strata to support the overlying strata, the surface is no longer affected by residual surface deformation, thus freeing up land resources above the old goaf areas.

[0003] Currently, there are two main strategies for controlling surface subsidence. One approach focuses on mining technology, such as leaving coal pillars, controlling mining height, and backfilling to control the height of overburden damage and thus control surface subsidence. The other approach involves grouting to fill the overburden separation area, reducing the space for overburden migration and controlling surface subsidence. However, both methods have limitations. The first method wastes some coal resources or increases underground mining operations, resulting in limited economic benefits. The second method requires coordination with underground mining; the overburden movement cycle in goaf areas is long, and if the grouting layer has poor bearing capacity, it cannot sustain long-term surface subsidence control, thus its effectiveness in controlling residual surface deformation is limited. Both methods are based on controlling surface subsidence in the early stages of mining or initial mining areas, and their effectiveness in controlling surface subsidence in existing goaf areas is limited. Therefore, it is necessary to conduct multi-level grouting on key bearing layers for reconstruction based on the characteristics of the overburden strata in goaf areas to control surface subsidence. Summary of the Invention

[0004] This invention provides a computer simulation-based method for optimizing the design of multi-level grouting in overburden, the main purpose of which is to solve the problem of low grouting accuracy during the optimization of multi-level grouting design in overburden.

[0005] To achieve the above objectives, this invention provides a computer simulation-based method for optimizing multi-level grouting design in overburden, comprising:

[0006] Collect geological exploration data of the goaf area, digitize the geological exploration data, and construct a three-dimensional geomechanical model of the goaf area;

[0007] Finite element numerical simulation was performed on the three-dimensional geomechanical model to calculate the development characteristics of each layer in the overlying strata of the goaf, and the grouting bearing layer of the goaf was constructed based on the development characteristics.

[0008] Analyze the fracture development characteristics and stress distribution in the grouting bearing layer;

[0009] Based on the fracture development characteristics and the stress distribution of the strata, the grouting simulation parameters of the grouting bearing layer are analyzed, and the grouting simulation parameters are iteratively optimized based on the rheological characteristics of the grouting material.

[0010] A multi-level, layered grouting simulation scheme for the grouting bearing layer is dynamically generated using optimized grouting simulation parameters.

[0011] The digital trajectory of the grouting borehole in the grouting bearing layer is generated according to the multi-level layered grouting simulation scheme, and grouting control simulation instructions are generated according to the digital trajectory.

[0012] Perform multi-level grouting simulation operations on the grouting bearing layer according to the grouting control simulation command, and output multi-level grouting simulation results.

[0013] Optionally, the step of digitally processing the geological exploration data to construct a three-dimensional geomechanical model of the goaf includes:

[0014] Perform borehole lithology identification on the geological exploration data to generate lithological distribution data of the overlying strata in the goaf.

[0015] The fracture zone boundary was calibrated in the overlying strata of the goaf to obtain the coordinate range of the subsidence zone.

[0016] Spatial interpolation calculations are performed on the overburden strata based on the coordinate range of the subsidence zone to generate a density cloud map of the overburden strata structure surface;

[0017] Based on the density cloud map and the lithological distribution data, a block discrete network model with fault constraints is constructed on the goaf.

[0018] The block discrete network model is used as a three-dimensional geomechanical model of the goaf.

[0019] Optionally, the step of performing finite element numerical simulation on the three-dimensional geomechanical model to calculate the development characteristics of each layer in the overlying strata of the goaf includes:

[0020] Based on the geological exploration data, the rock mass mechanics parameters corresponding to the three-dimensional geomechanics model are loaded, and an elastoplastic constitutive model of the overlying strata of the goaf is constructed according to the rock mass mechanics parameters.

[0021] The three-dimensional geomechanical model is divided into grid cells based on the overlying elastoplastic constitutive model to obtain a grid cell group;

[0022] Multi-condition numerical simulations were performed on the grid cell group to obtain the development characteristic parameters of each layer in the overburden strata, and the development characteristics of each layer in the overburden strata were determined based on the development characteristic parameters.

[0023] Optionally, the step of constructing the grouting bearing layer of the goaf based on the development characteristics includes:

[0024] The developmental features are mapped onto the three-dimensional geomechanical model to output the stratigraphic feature matrix;

[0025] Based on preset layer conditions, the layers in the layer feature matrix are filtered to obtain the candidate area of ​​the bearing layer;

[0026] The boundary range of the grouting bearing layer is determined according to the preset algorithm for the height of the water-conducting fracture zone;

[0027] The grouting bearing layer position is determined in the candidate area of ​​the bearing layer according to the boundary range.

[0028] Optionally, the analysis of the fracture development characteristics and stress distribution in the grouting bearing layer includes:

[0029] Based on the stratum coordinates of the grouting bearing layer, the fracture density parameter, fracture orientation parameter, and stress concentration factor parameter are extracted from the stratum feature matrix.

[0030] A directed graph model of fracture connectivity is generated based on the fracture density parameter and the fracture orientation parameter.

[0031] The fracture development characteristics of the grouting bearing layer are determined based on the path trends in the directed graph model of fracture connectivity.

[0032] Based on the stress concentration factor parameters, a layer stress gradient cloud map is generated, and the layer stress distribution of the grouting bearing layer is determined according to the cloud map trend of the layer stress gradient cloud map.

[0033] Optionally, the step of analyzing the grouting simulation parameters of the grouting bearing layer based on the fracture development characteristics and the layer stress distribution includes:

[0034] The number of horizontal layers in the grouting bearing layer is determined based on the thickness of the grouting bearing layer and the grouting diffusion radius.

[0035] The grouting location for each layer in the horizontal stratification group is determined based on the fracture development characteristics and the layer stress distribution.

[0036] The fracture space volume of each layer in the number of horizontal stratification groups is calculated based on the fracture connectivity in the fracture development characteristics.

[0037] The grouting volume corresponding to the grouting location is calculated by the volume of the fracture space and the preset rock mass fragmentation coefficient, and the grouting pressure at the grouting location is determined according to the stress distribution of the strata.

[0038] The grouting volume, the grouting pressure, and the grouting location are determined as grouting simulation parameters.

[0039] Optionally, the iterative optimization of the grouting simulation parameters based on the rheological characteristics of the grouting material includes:

[0040] The dynamic viscosity function is extracted from the rheological characteristics based on the fracture development characteristics of the grouting bearing layer.

[0041] Based on the dynamic viscosity function and the layer stress distribution, a slurry diffusion control equation is generated.

[0042] Perform a finite difference solution operation on the slurry diffusion control equation and output the grouting coverage radius;

[0043] Calculate the pressure compensation gradient based on the deviation between the grouting coverage radius and the preset target coverage area;

[0044] Based on the pressure compensation gradient iterative update of the grouting pressure parameters, when the coverage radius deviation of the preset number of iterations is less than or equal to the preset deviation threshold, the optimized grouting simulation parameters are output.

[0045] Optionally, the step of dynamically generating a multi-level layered grouting simulation scheme for the grouting bearing layer using optimized grouting simulation parameters includes:

[0046] The grouting bearing layer is divided into annular regions with increasing radial distances, and the grouting execution order is determined based on the annular regions.

[0047] The grouting sequence is determined based on the layer order of the grouting bearing layer.

[0048] A grouting sequence is generated based on the grouting distance execution order and the grouting up-down execution order;

[0049] A crack sealing simulation operation is performed on the grouting area of ​​the lower layer of the grouting bearing layer to construct a virtual boundary of the grout barrier layer;

[0050] Based on the virtual boundary of the grout barrier layer, a grout barrier layer simulation unit is constructed. When the strength of the grout barrier layer simulation unit meets the preset strength condition, a multi-level layered grouting simulation scheme is generated according to the grouting sequence and the optimized grouting simulation parameters.

[0051] Optionally, generating the digital trajectory of the grouting borehole in the grouting bearing layer according to the multi-level layered grouting simulation scheme includes:

[0052] The preliminary borehole axis of the grouting bearing layer is generated based on the grouting location in the multi-level layered grouting simulation scheme.

[0053] The drillability of the rock strata in the preliminary borehole axis is analyzed to obtain the target risk marker segment;

[0054] Based on the target risk marker segment and the preset gyro inclination data, the trajectory error compensation operation is performed on the preliminary borehole axis to obtain the borehole control point set;

[0055] The digital trajectory of the grouting bearing layer is generated based on the set of borehole control points.

[0056] Optionally, the step of performing multi-level grouting simulation operations on the grouting bearing layer according to the grouting control simulation command and outputting multi-level grouting simulation results includes:

[0057] The grouting control simulation command is used to monitor the orifice pressure and cumulative grouting volume of the grouting bearing layer in real time.

[0058] If the instantaneous rise rate of the orifice pressure exceeds a preset rise threshold and the cumulative grouting volume reaches a preset first grouting threshold, the stop grouting simulation command in the grouting control simulation command is triggered.

[0059] If the orifice pressure shows a decreasing trend within a preset time range or the cumulative grouting volume exceeds a preset second grouting threshold, the emergency positioning simulation command in the grouting control simulation command is triggered.

[0060] According to the emergency positioning simulation command, the weak surface of the overburden of the grouting bearing layer is scanned to obtain the grout leakage positioning point, and the sealing material injection simulation command is executed according to the grout leakage positioning point;

[0061] After the simulation command for injecting sealing material is executed, monitor whether the orifice pressure is stable. When the orifice pressure is stable, trigger the grouting control simulation command to resume grouting.

[0062] Grouting visualization diagrams are generated based on the stop grouting simulation command, the emergency positioning simulation command, and the resume grouting simulation command, respectively. The results of multi-level grouting simulation are determined through the grouting visualization diagrams.

[0063] This invention provides a comprehensive quantitative analysis and dynamic optimization of the entire process, from geological data acquisition to grouting simulation. It enables precise and controllable construction of multi-level grouting to reconstruct key bearing layers in goaf areas. The resulting new bearing layer effectively supports the overlying strata, significantly controls residual surface deformation, and liberates land resources above old goaf areas. Simultaneously, this process addresses technical problems in traditional methods, such as deviations caused by experience-based grouting layer selection, insufficient grouting or grout leakage due to unreasonable grouting simulation parameter design, and poor control of surface subsidence in old goaf areas. Through 3D modeling, finite element simulation, and rheological property optimization, it ensures the scientific nature of the grouting scheme and the stability of the bearing layer, improving the efficiency and reliability of mining-induced damage control. Therefore, the computer simulation-based multi-level grouting design optimization method proposed in this invention can solve the problem of low grouting accuracy during multi-level grouting design optimization. Attached Figure Description

[0064] Figure 1 This is a flowchart illustrating a computer simulation-based multi-level grouting design optimization method for overburden provided in an embodiment of the present invention.

[0065] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0066] It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0067] This application provides a computer-simulated multi-level grouting design optimization method for overburden. The execution entity of this computer-simulated multi-level grouting design optimization method includes, but is not limited to, at least one of the following electronic devices that can be configured to execute the method provided in this application: a server, a terminal, etc. In other words, the computer-simulated multi-level grouting design optimization method for overburden can be executed by software or hardware installed on a terminal device or a server device, and the software can be a blockchain platform. The server includes, but is not limited to, a single server, a server cluster, a cloud server, or a cloud server cluster. The server can be an independent server or a cloud server that provides basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communication, middleware services, domain name services, security services, content delivery networks (CDNs), and big data and artificial intelligence platforms.

[0068] Reference Figure 1The diagram shown is a flowchart illustrating a computer-simulated multi-level grouting design optimization method for overburden provided in an embodiment of the present invention. In this embodiment, the computer-simulated multi-level grouting design optimization method for overburden includes:

[0069] S1. Collect geological exploration data of the goaf area, digitize the geological exploration data, and construct a three-dimensional geomechanical model of the goaf area.

[0070] In this embodiment of the invention, the geological exploration data refers to a collection of geological information such as the distribution of rock strata, lithological characteristics, and fracture development status of the goaf area obtained through field measurement methods such as drilling and geophysical exploration.

[0071] In detail, geological exploration data of the goaf can be collected from a pre-stored storage area using computer statements with data scraping capabilities (such as Java statements, Python statements, etc.), where the storage area includes, but is not limited to, databases and blockchains.

[0072] Furthermore, in order to provide an accurate geological basis for subsequent overburden analysis and grouting scheme design, it is necessary to quantify and visualize the overburden characteristics of the goaf, so as to ensure that subsequent steps such as finite element numerical simulation and grouting layer determination have reliable data support.

[0073] In this embodiment of the invention, the three-dimensional geomechanical model refers to a three-dimensional digital model constructed based on geological exploration data of the goaf (such as core characteristics obtained by drilling, rock layer distribution information obtained by geophysical exploration, etc.). For example, in the scenario of an old goaf in a coal mine, the model can clearly present the alternating distribution of sandstone and mudstone, as well as the spatial location of fracture zones and bending subsidence zones, providing a visualized three-dimensional carrier for subsequent analysis of the stress state of the overburden.

[0074] In this embodiment of the invention, the step of digitally processing the geological exploration data to construct a three-dimensional geomechanical model of the goaf includes:

[0075] Perform borehole lithology identification on the geological exploration data to generate lithological distribution data of the overlying strata in the goaf.

[0076] The fracture zone boundary was calibrated in the overlying strata of the goaf to obtain the coordinate range of the subsidence zone.

[0077] Spatial interpolation calculations are performed on the overburden strata based on the coordinate range of the subsidence zone to generate a density cloud map of the overburden strata structure surface;

[0078] Based on the density cloud map and the lithological distribution data, a block discrete network model with fault constraints is constructed on the goaf.

[0079] The block discrete network model is used as a three-dimensional geomechanical model of the goaf.

[0080] In detail, laboratory analysis is conducted on borehole core samples from geological exploration data. X-ray diffraction is used to determine the mineral composition of the rock strata, and sonic logging data is combined to assess lithological hardness, generating lithological distribution data for the overlying strata of the goaf. The lithological identification results of each borehole are correlated with the borehole coordinates, and planar interpolation methods (such as inverse distance weighting) are used to draw planar and cross-sectional maps of the lithological distribution of the overlying strata in the goaf, clarifying the spatial distribution range of different lithologies (such as sandstone, mudstone, and limestone). For example, lithological identification data from 20 boreholes in a coal mine goaf shows that sandstone dominates within a depth of 100-150m, accounting for 70% of the stratum area, while mudstone dominates within a depth of 150-200m, accounting for 65%.

[0081] Specifically, the fracture zone refers to the area of ​​overlying strata in a goaf where numerous fractures have been created due to mining activities, but no collapse has occurred. The subsidence zone refers to the area of ​​overlying strata above the fracture zone where bending and subsidence have occurred. The boundary between the two needs to be determined through a combination of field measurements and theoretical calculations. The fracture zone boundary determination process includes: First, using geophysical methods (such as seismic wave reflection and acoustic wave transmission) to obtain the physical properties (such as wave velocity and resistivity) of the overlying strata. Because the wave velocity decreases in fracture-developed areas (e.g., the wave velocity in intact strata is 4000 m / s, while in fracture zones it drops to 2500 m / s), the wave velocity abrupt change interface is identified based on this, and the upper boundary of the fracture zone is preliminarily determined. Second, combining the "three zones"... The theoretical calculation formula for fracture zone height (such as the empirical formula for calculating fracture zone height based on mining height and lithology) is used to correct the geophysical exploration results. For example, if the mining height of a goaf is 5m, the fracture zone height is calculated to be 30m. Combined with the depth of the geophysical wave velocity change interface, the upper boundary of the fracture zone is finally determined. The coordinate range of the subsidence zone refers to the positional boundary of the subsidence zone in three-dimensional space. The upper and lower boundary coordinates of the subsidence zone are determined by measuring the surface subsidence points with a total station or GPS to invert the subsidence range of the rock strata.

[0082] Furthermore, the overburden strata data within the subsidence zone coordinate range are processed using the Kriging interpolation algorithm to generate a density contour map of the overburden strata structural surfaces. This contour map uses color gradients to represent the density of structural surfaces, with red areas representing densely packed areas. Based on the distribution of structural surfaces and lithological distribution data in the density contour map, finite element software is used to divide the overburden strata into block elements with fault constraints. The block discrete network model divides the overburden strata into several independent blocks, and the contact relationships between the blocks reflect the three-dimensional model of the rock strata structure. The construction process is as follows: First, the spatial distribution of structural planes is determined based on the density cloud map (e.g., the red high-density area corresponds to dense structural planes), and the structural planes are used as the dividing boundaries of the blocks; second, combined with lithological distribution data, the blocks within the same lithological region are made to have the same mechanical parameters (e.g., the elastic modulus of sandstone blocks is 25 GPa, and that of mudstone blocks is 8 GPa); finally, fault constraint conditions are introduced, that is, the blocks on both sides of the fault do not penetrate each other, and the relative displacement between the blocks must conform to the mechanical characteristics of the fault (e.g., the shear displacement coefficient along the fault strike is 0.2). Because the block discrete network model integrates lithological distribution (the basis of mechanical parameters), structural plane density (rock mass integrity characteristics), and fault constraints (boundary conditions), it fully covers the core information required for geomechanical analysis and can be directly used as a three-dimensional geomechanical model.

[0083] Furthermore, through multi-source data integration and numerical modeling, the geological characteristics of the goaf area were quantified and visualized, providing geometric boundaries and mechanical parameter carriers for finite element numerical simulation.

[0084] S2. Perform finite element numerical simulation on the three-dimensional geomechanical model to calculate the development characteristics of each layer of the overlying strata in the goaf, and construct the grouting bearing layer of the goaf based on the development characteristics.

[0085] In this embodiment of the invention, finite element numerical simulation refers to a numerical calculation method that discretizes a three-dimensional geomechanical model into a finite number of elements and simulates the stress state of the overburden strata by solving mechanical equilibrium equations; development characteristics refer to a set of parameters reflecting the stability of the overburden strata, such as fracture density, deformation, and stress concentration.

[0086] In this embodiment of the invention, the step of performing finite element numerical simulation on the three-dimensional geomechanical model to calculate the development characteristics of each layer in the overlying strata of the goaf includes:

[0087] Based on the geological exploration data, the rock mass mechanics parameters corresponding to the three-dimensional geomechanics model are loaded, and an elastoplastic constitutive model of the overlying strata of the goaf is constructed according to the rock mass mechanics parameters.

[0088] The three-dimensional geomechanical model is divided into grid cells based on the overlying elastoplastic constitutive model to obtain a grid cell group;

[0089] Multi-condition numerical simulations were performed on the grid cell group to obtain the development characteristic parameters of each layer in the overburden strata, and the development characteristics of each layer in the overburden strata were determined based on the development characteristic parameters.

[0090] In detail, rock mass mechanical parameters such as elastic modulus (e.g., 25 GPa for sandstone, 8 GPa for mudstone), Poisson's ratio (0.2-0.3), and internal friction angle (30°-40°) are extracted from geological exploration data and substituted into the Drucker-Prager elastoplastic constitutive equation to construct an overburden elastoplastic constitutive model. This model can realistically reflect the elastic deformation and plastic yielding characteristics of different lithological strata under stress. Based on the geometric characteristics of the strata reflected by the overburden elastoplastic constitutive model (e.g., the total overburden thickness is 100 m, divided into fracture zones, bending subsidence zones, etc., with each layer being 20-30 m thick), the basic scale of grid division is determined: fine grids are used in densely fractured areas (e.g., the red high-density areas identified by density cloud maps), with unit sizes set to 1 m × 1 m × 1 m, to accurately capture stress changes near fractures; in areas with homogeneous lithology, fine grids are used in areas with dense fracture development (e.g., the red high-density areas identified by density cloud maps), with unit sizes set to 1 m × 1 m × 1 m, to accurately capture stress changes near fractures; in areas with uniform .... For structurally intact areas (such as the blue low-density areas), a coarse mesh is used with a unit size of 5m×5m×5m to reduce computational load. Alternatively, based on the geometric characteristics of the overlying rock elastoplastic constitutive model, tetrahedral elements are used to divide the three-dimensional geomechanical model. The mesh size is set to 0.5m×0.5m in areas with dense fractures and 2m×2m in areas with intact rock strata, forming a mesh unit group containing 500,000 elements. The mesh unit group is a collection of several continuous and non-overlapping small units that discretize the three-dimensional geomechanical model. Each unit contains corresponding lithological information and mechanical parameters.

[0091] Specifically, multi-condition numerical simulation refers to setting different mining conditions (such as mining depth, goaf range, mining intensity, etc.) and calculating the mechanical response of the grid unit group using finite element software to comprehensively reflect the variation law of overburden strata under different scenarios. That is, multiple sets of simulation conditions are set: combined with the actual coal mining, three sets of conditions are set: mining depth of 300m and goaf range of 100m×50m; mining depth of 400m and goaf range of 150m×80m; mining depth of 500m and goaf range of 200m×100m. Each set of conditions simulates the overburden change process 100 days after mining. Then, the grid unit group is solved using finite element software: based on the overburden elastoplastic constitutive model, self-weight stress (calculated according to rock mass density of 2500kg / m³) and goaf boundary conditions (stress released to 0 within the goaf range) are applied to solve the stress, strain and displacement values ​​of each unit and output the development characteristic parameters of each unit. For example, at a mining depth of 300m, a certain unit is calculated to have a fracture density of 1.2 fractures / m², a settlement of 0.5m, a stress concentration factor of 1.3, and a plastic strain of 0.002. The average value of all unit parameters in the same stratum (such as a fracture zone) is then taken to obtain the development characteristic parameters of that stratum (such as an average fracture density of 1.5 fractures / m² and an average settlement of 0.8m in the fracture zone). The development characteristic parameters are quantitative indicators characterizing the state of the overburden strata, including fracture density, settlement, stress concentration factor, and plastic strain. The development characteristics are then determined based on these parameters: if the average fracture density of a certain stratum is >1 fracture / m² and the average plastic strain is >0.001, its development characteristic is judged as "the fractures are well developed and there is some plastic deformation"; if the average stress concentration factor is >1.5, it is judged as "significant stress concentration and poor stability".

[0092] Furthermore, by quantitatively analyzing the characteristics of overburden strata under different mining conditions, the limitations of relying on experience to judge the development status of strata were overcome, providing data support for the subsequent screening of grouting bearing strata. The development characteristic parameters are the core basis for screening candidate bearing strata areas, which is in line with the needs of coal mine overburden "three-zone" analysis.

[0093] In this embodiment of the invention, the grouting bearing layer refers to the target layer that has the conditions for grouting and can support the overlying rock layer.

[0094] In this embodiment of the invention, the step of constructing the grouting bearing layer of the goaf based on the development characteristics includes:

[0095] The developmental features are mapped onto the three-dimensional geomechanical model to output the stratigraphic feature matrix;

[0096] Based on preset layer conditions, the layers in the layer feature matrix are filtered to obtain the candidate area of ​​the bearing layer;

[0097] The boundary range of the grouting bearing layer is determined according to the preset algorithm for the height of the water-conducting fracture zone;

[0098] The grouting bearing layer position is determined in the candidate area of ​​the bearing layer according to the boundary range.

[0099] In detail, each developmental characteristic parameter (such as a fracture density of 1.2 fractures / m² and a stress concentration factor of 1.3 in a certain stratum) is bound to the coordinates of the corresponding block unit in the three-dimensional geomechanical model. Data gridding is used to divide the continuous three-dimensional space into 10m×10m×5m grid units. The average value of the developmental characteristic parameter in each grid unit is calculated (such as the fracture densities of 1.1, 1.3 and 1.2 fractures / m² in 3 block units in a certain grid, with an average value of 1.2 fractures / m²). A stratum characteristic matrix is ​​constructed with grid unit coordinates as rows and developmental characteristic parameters as columns. The preset stratigraphic conditions are screening criteria set according to the functional requirements of reconstructing key bearing strata. These mainly include the degree of fracture development, bearing capacity potential, and stability requirements, such as "the reconstructed key bearing strata must have a bearing capacity of 30~45MPa" and "the target area can be selected from fracture zones and flexural subsidence zones." For example, preset stratigraphic conditions might be: fracture density of 0.8-2.0 fractures / m² (ensuring grout penetration and sufficient space for cementation), stress concentration factor ≤1.5 (avoiding strata easily damaged due to excessive stress), and lithology of hard sandstone or limestone (the target strata are mostly hard sandstone). The stratigraphic feature matrix is ​​then traversed, and each grid cell is checked to see if it meets all preset conditions. For example, a grid cell with a fracture density of 1.5 fractures / m², a stress concentration factor of 1.2, and sandstone lithology is considered to meet the conditions; if another cell has a fracture density of 0.6 fractures / m² (below the lower limit), it is excluded. All grid cells that meet the conditions are spatially aggregated to form a continuous region, which is the candidate area for bearing strata.

[0100] Specifically, the water-conducting fracture zone refers to the area of ​​rock strata in the overlying strata of the goaf that has fractures due to mining and may conduct water. Its height is the vertical distance from the top of the goaf to the upper limit of the fracture zone. The spatial range of the fracture zone is quantitatively determined based on the water fracture zone height prediction method. The boundary range of the grouting bearing layer refers to the upper and lower limit coordinates of the bearing layer in three-dimensional space, used to constrain the spatial position of the final layer. For example, in a coal mine with a mining height of 5m, the upper limit of the water-conducting fracture zone is 54.7m from the top of the goaf. Combining the goaf burial depth, the absolute elevation is determined as follows: if the elevation of the goaf top is -400m, then the upper limit elevation of the water-conducting fracture zone is -400 + 54.7 ≈ -345.3m. Based on the condition that the target layer can be selected 30m above the fracture zone, the lower boundary of the grouting bearing layer is determined to be -345.3m + 10m = -335.3m (avoiding the core area of ​​the fracture zone), and the upper boundary is -335.3m. +20m = -315.3m (ensuring sufficient layer thickness) forms the boundary range; compare the spatial coordinates of the boundary range with the candidate area of ​​the bearing layer, and extract the part of the candidate area located within this vertical range. If the extracted area is a continuous 2000㎡ (meets engineering requirements) and the internal development characteristic parameters all meet the preset conditions (such as average crack density of 1.3 cracks / m², stress concentration factor of 1.1), it is determined to be valid, and the area is identified as the grouting bearing layer, and its three-dimensional coordinate range is output.

[0101] Furthermore, it can accurately pinpoint the grouting target layer that meets the requirements, solving the problem of "insufficient bearing capacity" or "ineffective grouting" caused by selecting layers based solely on experience. By superimposing the candidate area and the boundary range, the layer can be accurately located, providing a clear spatial carrier for the subsequent design of grouting simulation parameters.

[0102] S3. Analyze the fracture development characteristics and stress distribution of the grouting bearing layer.

[0103] In this embodiment of the invention, the fracture development characteristics refer to the spatial characteristics of fracture distribution density, orientation, and connectivity in the grouting bearing layer; the layer stress distribution refers to the magnitude and trend of stress values ​​within the layer.

[0104] In this embodiment of the invention, the analysis of the fracture development characteristics and stress distribution in the grouting bearing layer includes:

[0105] Based on the stratum coordinates of the grouting bearing layer, the fracture density parameter, fracture orientation parameter, and stress concentration factor parameter are extracted from the stratum feature matrix.

[0106] A directed graph model of fracture connectivity is generated based on the fracture density parameter and the fracture orientation parameter.

[0107] The fracture development characteristics of the grouting bearing layer are determined based on the path trends in the directed graph model of fracture connectivity.

[0108] Based on the stress concentration factor parameters, a layer stress gradient cloud map is generated, and the layer stress distribution of the grouting bearing layer is determined according to the cloud map trend of the layer stress gradient cloud map.

[0109] In detail, the fracture density parameter refers to the number of fractures per unit area in the grouting bearing layer, reflecting the density of fractures; the fracture orientation parameter refers to the angle between the fracture extension direction and the due north direction, characterizing the spatial extension trend of the fractures; the stress concentration factor parameter refers to the ratio of the local stress value to the average stress value within the layer, used to measure the degree of stress concentration, that is, to define the three-dimensional coordinate range of the grouting bearing layer, and to select entries whose coordinates fall within this range from the layer feature matrix; from these entries, the corresponding fracture density parameter, fracture orientation parameter, and stress concentration factor parameter are extracted. For example, in a grouting layer of a coal mine, the parameters corresponding to coordinates (100, 50, -320) are: fracture density 1.5 fractures / m², orientation 35°, and stress concentration factor 1.2.

[0110] Specifically, the directed graph model of fracture connectivity is a graphical model used to describe the connection relationship between fractures. It uses the endpoints of fractures as nodes and the extension direction of fractures as directed edges (the arrow direction is consistent with the fracture direction). The node density is determined based on the extracted fracture density parameters (higher density means denser nodes). The direction of the directed edges is determined according to the fracture orientation parameters (e.g., for a fracture with a 35° orientation, the arrow of the edge points in the 35° direction). Finally, a directed graph reflecting the interconnection of fractures is constructed. The extension direction and connectivity of most paths (e.g., over 70%) in the directed graph are analyzed. If the path mainly extends along the 30° direction and the connection rate between nodes exceeds 60%, the fracture development characteristic of the grouting bearing layer is determined to be "mainly oriented at 30° with good connectivity"; if the path directions are disordered and the connection rate is less than 30%, it is determined to be "dispersed orientations with poor connectivity". A stratigraphic stress gradient cloud map is a graphic representation of stress changes within a stratigraphic layer using color gradients. Typically, blue represents low stress (stress concentration factor < 1.0), green represents moderate stress (1.0 ≤ stress concentration factor ≤ 1.3), and red represents high stress (stress concentration factor > 1.3). The generation and analysis steps are as follows: the extracted stress concentration factor parameters are mapped onto a two-dimensional plane using coordinates, and an interpolation algorithm is used to generate the cloud map. Based on the color distribution trend in the cloud map (e.g., red areas are concentrated in the middle of the stratigraphic layer, and blue areas are distributed at the edges), the stratigraphic stress distribution characteristics are determined (e.g., "stress concentration in the middle, and gentler stress at the edges"). The overburden stress distribution then affects the grouting effect.

[0111] Furthermore, based on the characteristics of fracture development and the stress distribution in the strata, the internal structure of the grouting strata can be clearly determined, providing a basis for the design of subsequent grouting simulation parameters.

[0112] S4. Analyze the grouting simulation parameters of the grouting bearing layer based on the fracture development characteristics and the layer stress distribution, and iteratively optimize the grouting simulation parameters based on the rheological characteristics of the grouting material.

[0113] In this embodiment of the invention, the grouting simulation parameters refer to the key indicators that need to be controlled during the grouting process, including grouting volume, grouting pressure, grouting location, etc.

[0114] In this embodiment of the invention, the step of analyzing the grouting simulation parameters of the grouting bearing layer based on the fracture development characteristics and the layer stress distribution includes:

[0115] The number of horizontal layers in the grouting bearing layer is determined based on the thickness of the grouting bearing layer and the grouting diffusion radius.

[0116] The grouting location for each layer in the horizontal stratification group is determined based on the fracture development characteristics and the layer stress distribution.

[0117] The fracture space volume of each layer in the number of horizontal stratification groups is calculated based on the fracture connectivity in the fracture development characteristics.

[0118] The grouting volume corresponding to the grouting location is calculated by the volume of the fracture space and the preset rock mass fragmentation coefficient, and the grouting pressure at the grouting location is determined according to the stress distribution of the strata.

[0119] The grouting volume, the grouting pressure, and the grouting location are determined as grouting simulation parameters.

[0120] In detail, the number of horizontal stratification groups refers to the number of horizontal grouting layers that divide the grouting bearing layer along the vertical direction. This is used to achieve multi-level segmented grouting and ensure that the grout uniformly covers the target layer. The grouting diffusion radius refers to the maximum horizontal distance that the grout can naturally diffuse in the rock stratum fissures. It is determined by the characteristics of the grouting material (such as viscosity) and the permeability of the rock stratum. The thickness of the grouting bearing layer is extracted from the geological data, and the grouting diffusion radius is determined through laboratory grouting simulation tests. The number of horizontal stratification groups can be determined according to the formula: number of horizontal stratification groups = layer thickness ÷ (2 × grouting diffusion radius). Grouting location refers to the specific spatial coordinates of grouting holes set in each horizontal layer. It needs to be determined in combination with the characteristics of fracture development (such as fracture density and connectivity) and the stress distribution of the stratum (such as stress concentration) to improve grout filling efficiency and avoid the risk of grout leakage in areas with excessive stress. That is, within the horizontal layer, areas with a fracture density ≥1.2 fractures / m² and good connectivity are preferred, as grout can easily diffuse in such areas and form effective cementation. In combination with the stress distribution of the stratum, high stress areas with stress concentration factor >1.5 are avoided (to avoid rock fracture due to pressure superposition during grouting), and medium-low stress areas with stress concentration factor ≤1.3 are selected. Grouting holes are arranged in a 5m × 5m grid within the candidate area to ensure that the diffusion range of each grouting hole can cover the surrounding densely fractured area. Fracture connectivity refers to the degree to which fractures are interconnected within a horizontal stratum; fracture space volume refers to the total volume of grout that can be contained in all connected fractures within that stratum, and is a fundamental parameter for calculating grouting volume. In the statistical model, the proportion of connected fractures to the total fractures is used. The geometric parameters of the fractures within that stratum are extracted through a three-dimensional geomechanical model. The volume of a single fracture is calculated according to the formula "volume of a single fracture = length × average width × average height" and then summed. Therefore, fracture space volume = total fracture volume × fracture connectivity.

[0121] Specifically, the preset rock mass fragmentation coefficient refers to the ratio of the volume of the rock mass after fragmentation to its original volume, used to compensate for the volume change caused by rock mass compression; the grouting volume refers to the volume of grout required to fill the fracture space, i.e., grouting volume = fracture space volume × preset rock mass fragmentation coefficient. The grouting pressure refers to the pressure applied during grouting, which needs to be matched with the stress of the stratum to avoid grout leakage or insufficient filling. If the stress distribution of the stratum is a low-stress zone (stress concentration factor ≤ 1.0), a lower pressure (e.g., 0.5-1 MPa) is used to avoid excessive grout diffusion leading to grout leakage (the grouting pressure of the following stratum is lower, used to seal fractures); if it is a medium-stress zone (1.0 < stress concentration factor ≤ 1.3), a conventional pressure (e.g., 1.5-2 MPa) is used to ensure sufficient grout filling; if it is a high-stress zone (stress concentration factor > 1.3), a higher pressure (e.g., 2-3 MPa) is used to overcome the rock mass stress to ensure grout penetration (meeting the requirement of "pressurized grouting, different pressures for different strata"). Then, the grouting location coordinates (e.g., X=100, Y=50, Z=-320), corresponding grouting volume (108m³), and grouting pressure (1.8MPa) of each horizontal layer are bound together to form a structured parameter table to determine the grouting simulation parameters.

[0122] Furthermore, by determining precise grouting simulation parameters through quantitative calculations, the problems of "insufficient grouting" or "excessive waste" caused by reliance on experience in grouting are solved, meeting the requirements of precise pressurized grouting reinforcement. The grouting simulation parameters provide a data foundation for multi-level grouting sequence and process design, forming a technical chain of feature analysis, parameter calculation, and scheme implementation.

[0123] In this embodiment of the invention, the rheological characteristics of the grouting material refer to its flow and deformation characteristics under external force, specifically the variation law of viscosity with shear rate, time or temperature. By dynamically optimizing the grouting pressure, the problem of insufficient or excessive grouting caused by fixed parameters is solved, making the grout coverage more accurate.

[0124] In this embodiment of the invention, the iterative optimization of the grouting simulation parameters based on the rheological characteristics of the grouting material includes:

[0125] The dynamic viscosity function is extracted from the rheological characteristics based on the fracture development characteristics of the grouting bearing layer.

[0126] Based on the dynamic viscosity function and the layer stress distribution, a slurry diffusion control equation is generated.

[0127] Perform a finite difference solution operation on the slurry diffusion control equation and output the grouting coverage radius;

[0128] Calculate the pressure compensation gradient based on the deviation between the grouting coverage radius and the preset target coverage area;

[0129] Based on the pressure compensation gradient iterative update of the grouting pressure parameters, when the coverage radius deviation of the preset number of iterations is less than or equal to the preset deviation threshold, the optimized grouting simulation parameters are output.

[0130] In detail, the dynamic viscosity function is a mathematical expression describing the change in viscosity of grouting materials with shear rate. Rheological tests were conducted on the grouting materials (according to a preset ratio: 30% fly ash, 40% 1-10mm gangue, 20% cement, and 3% accelerator) using a rotational viscometer to determine different shear rates (e.g., 10s). -1 50s -1 100s -1 The viscosity value at a given point is used to obtain the rheological curve. If the fracture density at the grouting layer is 1.5 fractures / m² and the connectivity is 0.7 (moderate connectivity), then the interval that meets the moderate shear thinning characteristics is selected from the rheological curve. The data of this interval is fitted using the least squares method to obtain the dynamic viscosity function. ,in denoted as shear rate.

[0131] Specifically, the grout diffusion control equation is a mathematical model describing the relationship between the grout diffusion distance in the fracture and the grouting pressure, material viscosity, and layer stress. Therefore, the grout diffusion control equation is: ,in This represents the slurry diffusion distance. For grouting pressure, Let be the grout diffusion velocity, and 0.02 be an empirical coefficient. The higher the pressure, the lower the stress, and the lower the viscosity, the longer the diffusion distance. Therefore, the continuous governing equations are discretized into algebraic equations on grid nodes. A numerical method for obtaining the grout diffusion distance at each node is used through iterative calculation. The grouting coverage radius refers to the maximum horizontal distance the grout can diffuse under a certain grouting pressure, and is a key indicator for evaluating the grouting effect. The governing equations are then discretized into a difference scheme. ,in For the first Step diffusion distance, For the first The diffusion rate is calculated step by step, and the diffusion distance of each grid node is obtained by iterative calculation for 500 steps (simulating 500s diffusion). The maximum value is taken as the grouting coverage radius.

[0132] Furthermore, the preset target coverage area refers to the minimum coverage radius set to ensure complete filling of the crack. The difference between the actual coverage radius and the target radius is calculated as the deviation. The pressure compensation gradient is determined by the increase in grouting pressure required to compensate for the unit deviation. For example, if the pressure-coverage sensitivity of the grouting material is "for every 0.1 MPa increase in pressure, the coverage radius increases by 0.2 m", the deviation is 0.5 m, so the pressure compensation gradient is 0.5 MPa / m, and the total compensation pressure is 0.25 MPa. Iterative update refers to the cyclical process of adjusting the grouting pressure, recalculating the coverage radius, and verifying the deviation multiple times. The preset number of iterations (e.g., 3 times) and deviation threshold (e.g., 0.1 m) are to balance calculation accuracy and efficiency and ensure parameter optimization convergence.

[0133] Furthermore, by dynamically optimizing the grouting pressure, the problem of insufficient or excessive grouting caused by fixed parameters is solved, making the grout coverage more precise.

[0134] S5. Dynamically generate a multi-level layered grouting simulation scheme for the grouting bearing layer using the optimized grouting simulation parameters.

[0135] In this embodiment of the invention, the multi-level layered grouting simulation scheme is a construction guidance document that includes the grouting sequence and parameters for each step (location, pressure, quantity, and material ratio).

[0136] In this embodiment of the invention, the step of dynamically generating a multi-level layered grouting simulation scheme for the grouting bearing layer using optimized grouting simulation parameters includes:

[0137] The grouting bearing layer is divided into annular regions with increasing radial distances, and the grouting execution order is determined based on the annular regions.

[0138] The grouting sequence is determined based on the layer order of the grouting bearing layer.

[0139] A grouting sequence is generated based on the grouting distance execution order and the grouting up-down execution order;

[0140] A crack sealing simulation operation is performed on the grouting area of ​​the lower layer of the grouting bearing layer to construct a virtual boundary of the grout barrier layer;

[0141] Based on the virtual boundary of the grout barrier layer, a grout barrier layer simulation unit is constructed. When the strength of the grout barrier layer simulation unit meets the preset strength condition, a multi-level layered grouting simulation scheme is generated according to the grouting sequence and the optimized grouting simulation parameters.

[0142] In detail, the annular region is a concentric ring-shaped area centered on the grouting main well, divided by increasing radial distance. This is used to achieve uniform diffusion and coverage of the grout. The radial distance refers to the horizontal straight-line distance between the boundary of each annular region and the main well. Based on the location of the main well, n annular regions with increasing radial distances are divided according to the grout diffusion radius and the layer area. The division is based on the principle that "the width of each ring should not exceed twice the diffusion radius" (ensuring no overlap between adjacent ring grouting areas). The grouting execution order is clearly defined as "Ring 1 → Ring 2 → Ring 3," thus filling the distant fractures first to avoid excessive grout accumulation towards the near end and reduce the risk of grout leakage. The layer sequence refers to the vertical arrangement of the grouting bearing layers according to the vertical depth (lower layer, middle layer, upper layer), determined by the Z-coordinate in the three-dimensional geomechanical model (the greater the depth, the lower the layer). The grouting sequence refers to the order in which grouting is performed from bottom to top according to the grouting layers. The grouting sequence is determined as "lower layer → middle layer → upper layer," based on the principle that grouting the lower layer can first seal the cracks below to form a grout-isolating layer, preventing grout from leaking into the old goaf collapse zone below. The grouting sequence is a specific construction step sequence table that combines the near-far grouting sequence with the upper-lower grouting sequence. It is used to guide the logical sequence of multi-level layered grouting, such as near-far: Ring 1 → Ring 2 → Ring 3; upper-lower: Lower layer → Middle layer → Upper layer, combining these to generate the grouting sequence.

[0143] Specifically, the crack sealing simulation operation refers to predicting the effect of grout solidification and crack sealing after grouting in the lower layer using numerical simulation software. This involves inputting parameters for the lower layer grouting area: grouting pressure 0.5 MPa, grouting volume 80 m³, grout setting time 6 hours, and crack density 1.5 cracks / m². The simulation software then calculates the grout diffusion path and the strength distribution after solidification: the grout fills the first to third ring cracks in the lower layer under a pressure of 0.5 MPa. After 6 hours, the stone mass strength reached 3MPa (meeting the strength requirements of the grout barrier simulation unit). The boundary of the closed area formed was Z=-330--320m, X=70-130m, Y=30-90m. This is the virtual boundary of the grout barrier simulation unit. The lower layer grouting was carried out according to the virtual boundary of the grout barrier simulation unit: grouting was carried out in the 1st-3rd rings of the lower layer according to the sequence steps 1-3, using quick-setting material (20% cement content, 3% quick-setting agent). The grout diffusion was monitored in real time (deviation from the simulated path ≤0.5m). After 6 hours, samples were taken through boreholes for testing. The strength of the grout mass was measured, and the result was 3.2 MPa (≥3 MPa, meeting the preset conditions). The formation of the grout barrier simulation unit was confirmed, and a multi-level layered grouting simulation scheme was generated: with the grouting sequence as the framework, the parameters of each step were bound: lower layer: pressure 0.5 MPa, material quick-setting type (3% quick-setting agent), volume 80-100 m³; middle layer: pressure 2 MPa, material conventional type (1% quick-setting agent), volume 95-110 m³; upper layer: pressure 2.5 MPa, material high-strength type (30% cement admixture), volume 90-105 m³.

[0144] For example, the multi-level grouting method involves segmented and layered grouting, with the main grouting well as a reference, grouting proceeds sequentially from far to near and from bottom to top. The initial grouting material in the lower layer needs to quickly seal fractures and form a grout-isolating layer simulation unit, depending on the intended use. The grouting pressure is low, and the grouting material requires a fast setting rate. After the lower layer forms a stable grout-isolating layer simulation unit through grouting from far to near, the grouting of the next layer begins, proceeding from far to near. The grouting pressure is increased to the conventional grouting pressure, and the grouting material ratio is adjusted. After the grouting of this layer is completed, the grouting of the next layer begins, proceeding from far to near, until all layers are grouted.

[0145] Furthermore, the grouting sequence is based on layered orderly grouting and grout barrier simulation unit protection, annular region division and layer sequence. The simulation and construction of grout barrier simulation unit ensures that the sequence can be executed safely. The final generated scheme is a comprehensive application of all the aforementioned optimized parameters.

[0146] S6. Generate the digital trajectory of the grouting borehole in the grouting bearing layer according to the multi-level layered grouting simulation scheme, and generate grouting control simulation instructions according to the digital trajectory.

[0147] In this embodiment of the invention, the digital trajectory refers to a smooth curve fitted based on a set of control points, which is used to guide directional drilling construction and ensure that the grouting pipe accurately reaches each grouting position.

[0148] In this embodiment of the invention, generating the digital trajectory of the grouting borehole in the grouting bearing layer according to the multi-level layered grouting simulation scheme includes:

[0149] The preliminary borehole axis of the grouting bearing layer is generated based on the grouting location in the multi-level layered grouting simulation scheme.

[0150] The drillability of the rock strata in the preliminary borehole axis is analyzed to obtain the target risk marker segment;

[0151] Based on the target risk marker segment and the preset gyro inclination data, the trajectory error compensation operation is performed on the preliminary borehole axis to obtain the borehole control point set;

[0152] The digital trajectory of the grouting bearing layer is generated based on the set of borehole control points.

[0153] In detail, the grouting location refers to the three-dimensional coordinates of each grouting point specified in the multi-level layered grouting simulation scheme, including the coordinates of the annular area in the horizontal direction and the layer coordinates in the vertical direction; the preliminary borehole axis refers to the straight path from the ground drilling station to the grouting location, which is used to initially plan the direction of the borehole, that is, to determine the coordinates of the ground drilling station as the starting point of the borehole, and for each grouting location, to calculate the path from the starting point to the ending point using the spatial straight line equation to obtain the angular parameters of the preliminary borehole axis (such as a depression angle of 75° and an azimuth angle of 30°).

[0154] Specifically, rock strata drillability refers to the ease with which a rock stratum can be penetrated by a borehole, usually expressed as a drillability grade (1-10, with higher grades indicating greater difficulty in drilling). From the lithological distribution data of the three-dimensional geomechanical model, the rock strata types and corresponding drillability grades traversed by the initial borehole axis are extracted. For example, the axis in the Z=-100-150m section is limestone (drillability grade 7), and the Z=-200-250m section is densely fractured sandstone (fracture density 2.5 fractures / m²). These sections are marked and defined as target risk marker sections. For these target risk marker sections, a preset compensation coefficient is used: when the drillability grade is 7, the deviation angle is adjusted by 0.5° every 10m to avoid excessive deviation. Combining gyro tilt data (e.g., if the actual deviation angle is 3°, exceeding the allowable error ±2°), the compensation amount is calculated: a 0.3° adjustment in the opposite direction is required to keep the deviation within 2°. Every 50m along the compensated axis... Set a control point, record its three-dimensional coordinates to form a set of borehole control points, ensure that the lines connecting the points are smooth, and then fit the control points in the set of borehole control points to obtain a digital trajectory.

[0155] Furthermore, the grouting control simulation commands are a set of specific instructions guiding drilling and grouting operations, including drilling parameters (angle, speed) and grouting simulation parameters (pressure, volume, material ratio). For example, the drilling stage commands specify the drilling angle, drilling speed, and drill bit type for each control point; the grouting stage commands specify the following according to the grouting sequence: the lower layer, first ring, grouting pressure is 0.5 MPa, grouting volume is 80 m³, and the material is quick-setting (3% quick-setting agent); the middle layer, second ring, pressure is 2 MPa, volume is 95 m³, and the material is conventional (1% quick-setting agent); the monitoring commands monitor the borehole pressure (deviation ≤ 0.1 MPa) and grouting volume (cumulative error ≤ 5%) in real time. If the pressure drops sharply (e.g., < 0.3 MPa), the "stop grouting + seal weak surfaces" command is triggered.

[0156] S7. Perform multi-level grouting simulation operation on the grouting bearing layer according to the grouting control simulation command, and output the multi-level grouting simulation results.

[0157] In this embodiment of the invention, the multi-level grouting simulation result refers to a comprehensive evaluation of the grouting effect, such as coverage area, filling density, and whether the strength requirements of the bearing layer are met.

[0158] In this embodiment of the invention, the step of performing multi-level grouting simulation operation on the grouting bearing layer according to the grouting control simulation command and outputting multi-level grouting simulation results includes:

[0159] The grouting control simulation command is used to monitor the orifice pressure and cumulative grouting volume of the grouting bearing layer in real time.

[0160] If the instantaneous rise rate of the orifice pressure exceeds a preset rise threshold and the cumulative grouting volume reaches a preset first grouting threshold, the stop grouting simulation command in the grouting control simulation command is triggered.

[0161] If the orifice pressure shows a decreasing trend within a preset time range or the cumulative grouting volume exceeds a preset second grouting threshold, the emergency positioning simulation command in the grouting control simulation command is triggered.

[0162] According to the emergency positioning simulation command, the weak surface of the overburden of the grouting bearing layer is scanned to obtain the grout leakage positioning point, and the sealing material injection simulation command is executed according to the grout leakage positioning point;

[0163] After the simulation command for injecting sealing material is executed, monitor whether the orifice pressure is stable. When the orifice pressure is stable, trigger the grouting control simulation command to resume grouting.

[0164] Grouting visualization diagrams are generated based on the stop grouting simulation command, the emergency positioning simulation command, and the resume grouting simulation command, respectively. The results of multi-level grouting simulation are determined through the grouting visualization diagrams.

[0165] In detail, the orifice pressure refers to the grout pressure at the outlet of the grouting pipe, reflecting the diffusion resistance of the grout in the fracture; the cumulative grouting volume refers to the total volume of grout injected from the start of grouting to the current moment, used to determine the degree of fracture filling, i.e., pressure sensors and electromagnetic flowmeters are installed at the orifice of the grouting pipe to collect pressure and flow data in real time. The instantaneous rise rate of the orifice pressure refers to the change in pressure per unit time, reflecting the saturation degree of fracture filling by grout. The preset rise threshold is a critical value set according to the solidification characteristics of the grout. The first grouting threshold refers to the theoretically calculated amount of grout required to fill the fracture. When the instantaneous rise rate of the orifice pressure exceeds the preset rise threshold and the cumulative grouting volume reaches the preset first grouting threshold, it is determined that the fracture has been filled, the control system automatically sends a stop grouting simulation command, shuts down the grouting pump, and records the current status (e.g., "10:30 Stop grouting, pressure 2.0MPa, volume 108m³").

[0166] Specifically, the downward trend of orifice pressure refers to a continuous decrease in pressure within a preset time range, which may be due to grout leakage caused by the penetration of fractures. The preset second grouting threshold is the maximum allowable grouting volume. Exceeding this value indicates that the grout has not been effectively filled (possibly due to grout leakage). The control system then triggers an emergency positioning simulation command, automatically activating the ground-penetrating radar (detection accuracy 0.5m) to scan the grouting bearing layer, focusing on identifying areas with dense fractures. The ground-penetrating radar scan shows abnormal reflection signals in area A (dense fractures). Combined with borehole television observation, the grout leakage location point is determined. The weak surface of the overburden refers to the area with dense fractures and low strength in the overburden layer, which is the main channel for grout leakage. The grout leakage location point refers to the specific three-dimensional coordinates of grout leakage in the weak surface. The control system then sends a command to inject sealing material into the location point through the branch grouting pipe and monitors whether the orifice pressure is stable after the sealing material injection simulation command is executed. That is, if the orifice pressure fluctuation within a preset time is less than the preset fluctuation threshold (e.g., stable at 1.7-1.8MPa), it indicates that the sealing is effective and the fracture leakage channel has been closed. The command to resume grouting simulation refers to the command to restart the grouting operation. The control system sends a resumption command, the grouting pump restarts, and grouting continues according to the original parameters. The current cumulative volume is recorded (85+5=90m³ has been injected, and 18m³ remains).

[0167] For example, the ground directional drilling grouting method is a pressurized grouting method, which requires real-time monitoring of grouting pressure and grouting volume. If the borehole pressure rises rapidly and the grouting volume reaches the preset value, grouting is stopped; if the borehole pressure remains stable, grouting continues; if the borehole pressure drops rapidly or the grouting volume exceeds the preset value, grout leakage may occur, and grouting must be stopped immediately. In case of grout leakage, the leakage location must be located, the weak surface of the overburden rock must be sealed, and normal grouting can then proceed; in case of grout leakage, borehole observation must be conducted to locate the leakage location, and grouting can continue after the leakage location is sealed.

[0168] Furthermore, grouting visualization refers to grout diffusion cloud maps, pressure distribution contour maps, and command execution time axes generated by computer software. Visualization maps are generated according to different commands, such as the stop command map: showing the pressure surge point (10:30), the cumulative volume of 108m³, and the grout diffusion cloud map covering the target area (red indicates the filling area); the emergency command map: marking the grout runoff location points (97m, 57m, -325m), and showing the pressure change curves before and after the sealing; the recovery command map: showing the pressure stable section (10:45-10:50) and the remaining grout diffusion path.

[0169] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the present invention.

[0170] Therefore, the embodiments should be regarded as exemplary and non-limiting in all respects. The scope of the invention is not limited to the foregoing description, and all variations within the meaning and scope of equivalents falling within the protection scope are intended to be included in the invention.

[0171] The embodiments of this application can acquire and process relevant data based on artificial intelligence technology. Artificial intelligence (AI) refers to the theories, methods, technologies, and application systems that use digital computers or machines controlled by digital computers to simulate, extend, and expand human intelligence, perceive the environment, acquire knowledge, and use that knowledge to obtain optimal results.

[0172] Furthermore, it is clear that the word "comprising" does not exclude other units or steps, and the singular does not exclude the plural. Multiple units or systems stated in a system claim may also be implemented by a single unit or system through software or hardware. The terms "first," "second," etc., are used to indicate names and do not indicate any specific order.

[0173] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A computer simulation-based method for optimizing multi-level grouting design in overburden rock, characterized in that, The method includes: Geological exploration data of the goaf area is collected, and the data is digitized to construct a three-dimensional geomechanical model of the goaf area. This includes: performing borehole lithology identification on the geological exploration data to generate lithological distribution data of the overlying strata in the goaf area; calibrating the fracture zone boundaries of the overlying strata in the goaf area to obtain the coordinate range of the subsidence zone; performing spatial interpolation calculations on the overlying strata based on the coordinate range of the subsidence zone to generate a density cloud map of the overlying strata structure; constructing a block discrete network model of the goaf area with fault constraints based on the density cloud map and the lithological distribution data; and using the block discrete network model as the three-dimensional geomechanical model of the goaf area. Finite element numerical simulation was performed on the three-dimensional geomechanical model to calculate the development characteristics of each layer in the overlying strata of the goaf, and the grouting bearing layer of the goaf was constructed based on the development characteristics. Analyze the fracture development characteristics and stress distribution in the grouting bearing layer; Based on the fracture development characteristics and the stress distribution of the strata, the grouting simulation parameters of the grouting bearing layer are analyzed, and the grouting simulation parameters are iteratively optimized based on the rheological characteristics of the grouting material. A multi-level layered grouting simulation scheme is dynamically generated for the grouting bearing layer using optimized grouting simulation parameters. This includes: dividing the grouting bearing layer into annular regions with increasing radial distances, and determining the grouting execution order based on the annular regions; determining the grouting execution order based on the layer order of the grouting bearing layer; generating a grouting sequence based on the grouting execution order and the grouting execution order; performing a crack sealing simulation operation on the lower grouting area of ​​the grouting bearing layer to construct a virtual boundary of the grout barrier layer; constructing a grout barrier layer simulation unit based on the virtual boundary of the grout barrier layer; and generating a multi-level layered grouting simulation scheme according to the grouting sequence and optimized grouting simulation parameters when the strength of the grout barrier layer simulation unit meets a preset strength condition. The digital trajectory of the grouting borehole in the grouting bearing layer is generated according to the multi-level layered grouting simulation scheme, and grouting control simulation instructions are generated according to the digital trajectory. Perform multi-level grouting simulation operations on the grouting bearing layer according to the grouting control simulation command, and output multi-level grouting simulation results.

2. The method for designing and optimizing multi-level grouting in overburden based on computer simulation as described in claim 1, characterized in that, The finite element numerical simulation of the three-dimensional geomechanical model, calculating the development characteristics of each layer in the overlying strata of the goaf, includes: Based on the geological exploration data, the rock mass mechanics parameters corresponding to the three-dimensional geomechanics model are loaded, and an elastoplastic constitutive model of the overlying strata of the goaf is constructed according to the rock mass mechanics parameters. The three-dimensional geomechanical model is divided into grid cells based on the overlying elastoplastic constitutive model to obtain a grid cell group; Multi-condition numerical simulations were performed on the grid cell group to obtain the development characteristic parameters of each layer in the overburden strata, and the development characteristics of each layer in the overburden strata were determined based on the development characteristic parameters.

3. The method for designing and optimizing multi-level grouting in overburden based on computer simulation as described in claim 1, characterized in that, The construction of the grouting bearing layer in the goaf based on the development characteristics includes: The developmental features are mapped onto the three-dimensional geomechanical model to output the stratigraphic feature matrix; Based on preset layer conditions, the layers in the layer feature matrix are filtered to obtain the candidate area of ​​the bearing layer; The boundary range of the grouting bearing layer is determined according to the preset algorithm for the height of the water-conducting fracture zone; The grouting bearing layer position is determined in the candidate area of ​​the bearing layer according to the boundary range.

4. The method for designing and optimizing multi-level grouting in overburden based on computer simulation as described in claim 3, characterized in that, The analysis of the fracture development characteristics and stress distribution in the grouting bearing layer includes: Based on the stratum coordinates of the grouting bearing layer, the fracture density parameter, fracture orientation parameter, and stress concentration factor parameter are extracted from the stratum feature matrix. A directed graph model of fracture connectivity is generated based on the fracture density parameter and the fracture orientation parameter. The fracture development characteristics of the grouting bearing layer are determined based on the path trends in the directed graph model of fracture connectivity. Based on the stress concentration factor parameters, a layer stress gradient cloud map is generated, and the layer stress distribution of the grouting bearing layer is determined according to the cloud map trend of the layer stress gradient cloud map.

5. The method for designing and optimizing multi-level grouting in overburden based on computer simulation as described in claim 1, characterized in that, The analysis of grouting simulation parameters for the grouting bearing layer based on the fracture development characteristics and the layer stress distribution includes: The number of horizontal layers in the grouting bearing layer is determined based on the thickness of the grouting bearing layer and the grouting diffusion radius. The grouting location for each layer in the horizontal stratification group is determined based on the fracture development characteristics and the layer stress distribution. The fracture space volume of each layer in the number of horizontal stratification groups is calculated based on the fracture connectivity in the fracture development characteristics. The grouting volume corresponding to the grouting location is calculated by the volume of the fracture space and the preset rock mass fragmentation coefficient, and the grouting pressure at the grouting location is determined according to the stress distribution of the strata. The grouting volume, the grouting pressure, and the grouting location are determined as grouting simulation parameters.

6. The method for designing and optimizing multi-level grouting in overburden based on computer simulation as described in claim 1, characterized in that, The iterative optimization of the grouting simulation parameters based on the rheological characteristics of the grouting material includes: The dynamic viscosity function is extracted from the rheological characteristics based on the fracture development characteristics of the grouting bearing layer. Based on the dynamic viscosity function and the layer stress distribution, a slurry diffusion control equation is generated. Perform a finite difference solution operation on the slurry diffusion control equation and output the grouting coverage radius; Calculate the pressure compensation gradient based on the deviation between the grouting coverage radius and the preset target coverage area; Based on the pressure compensation gradient iterative update of the grouting pressure parameters, when the coverage radius deviation of the preset number of iterations is less than or equal to the preset deviation threshold, the optimized grouting simulation parameters are output.

7. The method for designing and optimizing multi-level grouting in overburden based on computer simulation as described in claim 1, characterized in that, The step of generating the digital trajectory of the grouting borehole in the grouting bearing layer according to the multi-level layered grouting simulation scheme includes: The preliminary borehole axis of the grouting bearing layer is generated based on the grouting location in the multi-level layered grouting simulation scheme. The drillability of the rock strata in the preliminary borehole axis is analyzed to obtain the target risk marker segment; Based on the target risk marker segment and the preset gyro inclination data, the trajectory error compensation operation is performed on the preliminary borehole axis to obtain the borehole control point set; The digital trajectory of the grouting bearing layer is generated based on the set of borehole control points.

8. The method for designing and optimizing multi-level grouting in overburden based on computer simulation as described in claim 1, characterized in that, The step of performing multi-level grouting simulation operations on the grouting bearing layer according to the grouting control simulation command, and outputting multi-level grouting simulation results, includes: The grouting control simulation command is used to monitor the orifice pressure and cumulative grouting volume of the grouting bearing layer in real time. If the instantaneous rise rate of the orifice pressure exceeds a preset rise threshold and the cumulative grouting volume reaches a preset first grouting threshold, the stop grouting simulation command in the grouting control simulation command is triggered. If the orifice pressure shows a decreasing trend within a preset time range or the cumulative grouting volume exceeds a preset second grouting threshold, the emergency positioning simulation command in the grouting control simulation command is triggered. According to the emergency positioning simulation command, the weak surface of the overburden of the grouting bearing layer is scanned to obtain the grout leakage positioning point, and the sealing material injection simulation command is executed according to the grout leakage positioning point; After the simulation command for injecting sealing material is executed, monitor whether the orifice pressure is stable. When the orifice pressure is stable, trigger the grouting control simulation command to resume grouting. Grouting visualization diagrams are generated based on the stop grouting simulation command, the emergency positioning simulation command, and the resume grouting simulation command, respectively. The results of multi-level grouting simulation are determined through the grouting visualization diagrams.

Images