Sandstone grotto fracture water seepage risk assessment method based on numerical simulation

Abstract

The invention discloses a sandstone grotto fracture water seepage risk assessment method based on numerical simulation, and the method comprises the steps: S1, obtaining a three-dimensional fracture network and fracture data based on a fracture network model; s2, performing simulation calculation on the seepage conditions in the plurality of fractures to obtain a flow vector simulation value and a permeability vector simulation value of an outlet of each fracture; s3, establishing a fracture flow prediction model; s4, calculating the average contribution degree of each fracture permeability vector to the fracture outlet flow vector prediction value by adopting an LRP method; s5, determining a main seepage trajectory, and performing risk assessment by taking the fracture where the main seepage trajectory is located and / or the fracture with the maximum average contribution degree as a target fracture; according to the method, simulation operation is carried out on the basis of numerical simulation in combination with the neural network, then the fractures which have great influence on the grotto are selected from the fracture network to be analyzed through the LRP method, risk assessment is carried out, and the degree and influence of the fracture water seepage phenomenon in the grotto rock mass can be accurately assessed.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of grotto cultural relic protection, in particular to a sandstone grotto fissure water seepage risk assessment method based on numerical simulation. BACKGROUND

[0002] Under the long-term natural weathering and hydrogeological action, the fissure water seepage disease of sandstone grotto cultural relics seriously affects its stability and preservation condition. Fissure water seepage leads to problems such as rock softening, salt precipitation and weathering peeling, which seriously threatens the safety of grotto, and it is urgently needed to study, assess and prevent in advance.

[0003] The existing research and assessment methods mainly rely on traditional hydrogeological investigation and qualitative analysis, lack scientific and quantitative analysis means, and it is difficult to accurately identify the source and seepage trajectory of fissure water seepage disease, which leads to inaccurate analysis of the seepage mechanism of fissure water seepage and risk assessment. Therefore, a scientific, systematic and quantitative evaluation method is needed to clarify the seepage mechanism of fissure water seepage of sandstone grotto cultural relics, so as to improve the pertinence and effectiveness of grotto water seepage disease treatment. SUMMARY

[0004] In order to solve the above problems, the present application provides a sandstone grotto fissure water seepage disease assessment method based on numerical simulation, which can accurately identify the source and seepage trajectory of fissure water seepage disease, and provide technical support for grotto protection and water seepage disease treatment.

[0005] A sandstone grotto fissure water seepage risk assessment method based on numerical simulation, comprising the following steps:

[0006] S1, simulating a plurality of fissures in the sandstone grotto based on a fissure network model to obtain a three-dimensional fissure network and fissure data;

[0007] S2, according to the three-dimensional fissure network and fissure data, sequentially using constrained Delaunay triangulation method and finite volume method to simulate and calculate the seepage condition of the plurality of fissures, to obtain each fissure outlet flow vector simulation value and permeability vector simulation value;

[0008] S3, based on the fissure outlet flow vector simulation value and the permeability vector simulation value, a fissure flow prediction model is established; the fissure flow prediction model is used to represent the mapping relationship between the fissure permeability vector and the fissure outlet flow vector, the input of the fissure flow prediction model is the fissure permeability vector (the fissure permeability vector can be the true value measured, or the permeability vector simulation value), and the output is the fissure outlet flow vector prediction value;

[0009] S4, calculate the average contribution of each fracture permeability vector to the predicted value of fracture outlet flow vector by using the LRP method; then perform particle transport simulation operation on the water particles in the multiple fractures by using the dfnTrans method in the dfnWorks software to obtain the particle number in the fractures;

[0010] S5, determine the main seepage trajectory and the fracture where the main seepage trajectory is located according to the particle number and the average contribution, and take the fracture where the main seepage trajectory is located and / or the fracture with the largest average contribution of fracture permeability vector to the predicted value of fracture outlet flow vector as the target fracture; if the value of the seepage parameter in the target fracture exceeds the standard threshold value, it is considered that the selected range of sandstone grotto has high risk and preventive measures need to be taken, otherwise, it is considered that the sandstone grotto has low risk and preventive measures do not need to be taken.

[0011] Description: Through the above method, seepage data in fractures that are considered undetectable can be obtained by fracture seepage simulation in numerical simulation, providing more data for subsequent simulation and calculation. At the same time, the flow law of water in fractures can be simulated to evaluate the hydraulic characteristics of fractures. Particle transport simulation can visually display the motion trajectory and distribution of water particles in fractures, revealing the trajectory and diffusion law of water flow, which helps to understand the complexity of fracture network. Through the establishment of a neural network model and the use of the LRP method, key seepage trajectories in the fracture network can be identified and quantified based on numerical simulation, so that the flow mechanism of water in fractures can be more accurately understood, and potential risk areas can be identified. The above method not only helps to improve the understanding of fracture seepage behavior, but also provides a scientific basis for seepage control and risk management.

[0012] Further, the method for constructing the fracture network model in S1 comprises:

[0013] S1-1, randomly generate a plurality of simulated fracture data by using the Monte Carlo method;

[0014] S1-2, use the actual fracture data obtained by exploration and the simulated fracture data to perform geometric modeling by using the dfnWorks software;

[0015] S1-3, determine whether the fracture density in the geometric modeling meets the P 32 density requirement; if not, perform test correction on the simulated fracture data by using the FRAM method; then perform the steps of S1-1 and S1-2 in a loop until the fracture density in the geometric modeling meets the P 32 density requirement, and stop to obtain a three-dimensional fracture network; if yes, no test correction is needed for the simulated fracture data, and a three-dimensional fracture network is obtained.

[0016] Illustration: The model established by the above method can truly reflect the distribution and characteristics of the fissures in the sandstone grotto, improve the accuracy of the simulation, and help reveal the seepage law in the fissures and the flow mechanism of water in the fissures.

[0017] Further, the fissure data are actual fissure data and / or simulated fissure data, and the fissure data include fissure geometric shape data, spatial position data of the fissure, fissure permeability data, fluid property data, boundary condition and initial condition data. Specifically, the fissure geometric shape data include the length, width and depth of the fissure; the spatial position of the fissure (coordinates of the fissure in three-dimensional space, and spatial attitude information such as the strike, dip and dip angle of the fissure); the fissure permeability data (permeability of the fissure, and variation range of the permeability); the fluid property data (physical properties such as density and viscosity of the fluid, and chemical properties of the fluid); and the boundary condition and initial condition data (boundary shape of the fissure network, and fluid pressure or flow on the boundary).

[0018] Further, the step S3 of establishing the fissure flow prediction model comprises the following steps.

[0019] S3-1, obtaining a plurality of sets of fissure outlet flow vectors and permeability vectors as a data set;

[0020] S3-2, dividing the data set into a training set, a validation set and a test set, training the fissure flow prediction model constructed by using the training set; in the training process, the Adam algorithm is selected for optimization, and the mean square error is used as the loss function; and the neural network adopts a multi-branch feedforward neural network.

[0021] S3-3, after the training of the fissure flow prediction model is completed, the validation set is used for verification, the parameters of the fissure flow prediction model are adjusted according to the result of the validation set, the test set is used for testing the fissure flow prediction model and compared with the true value, so as to evaluate the model performance of the fissure flow prediction model.

[0022] Illustration: The above method can establish a neural network model, select the mean square error as the loss function and the Adam algorithm as the optimization method, effectively optimize the model parameters, and improve the prediction accuracy of the model. The above method sets the number of input layer nodes equal to the total number of fissures, and determines the number of branches according to the number of fissures intersecting with the rock wall, which can effectively capture the complexity of the fissure network. Through the combination of three hidden layers and ReLU activation function, the nonlinear expression ability of the model is ensured, and the calculation efficiency is maintained. The batch size of 4 realizes good parallel computing and memory management balance in the training process.

[0023] Further, the method for sequentially adopting the constrained Delaunay triangulation method and the finite volume method to simulate and calculate the seepage in the plurality of fissures to obtain the outlet flow vector and the permeability vector of each fissure comprises the following steps: first, the grid in the three-dimensional fissure network is divided into a triangular grid by using the constrained Delaunay triangulation method; then, the pressure value of the grid node is set as a variable in the solving; and then, the Voronoi polygon control body is established based on the triangular grid, and the finite volume method is used to solve the Voronoi polygon control body to obtain the outlet flow vector simulation value and the permeability vector simulation value of each fissure.

[0024] Description: The above method can complete the simulation and calculation of the seepage in the three-dimensional fissure network.

[0025] Further, in S4, the method for calculating the LRP method comprises the following steps:

[0026] S4-1, inputting a fissure permeability vector sample into the fissure flow prediction model to obtain a fissure outlet flow vector prediction value; wherein the fissure permeability vector sample comprises the permeability component of each fissure in the fissure network;

[0027] S4-2, for one fissure permeability vector sample, according to the alpha-beta rule, the fissure outlet flow vector prediction value is reversely distributed to the permeability component of each fissure in the fissure permeability vector of the fissure flow prediction model input layer layer by layer, to obtain a contribution degree vector c i of each component in the permeability vector to the fissure outlet flow vector prediction value.

[0028] S4-3, selecting n fissure permeability vector samples, repeating S3-1 and S3-2 n times to obtain n groups of contribution degree vectors; for each fissure i, the average contribution degree in the n fissure permeability vector samples is calculated by formula (1):

[0029]

[0030] In the formula, is the average contribution degree of fissure i, n is the number of fissure permeability vector samples, c i is the contribution degree vector of fissure i.

[0031] Description: The above method can quantitatively evaluate the relative importance of each fissure in the seepage process by inputting the fissure permeability vector sample multiple times and calculating the average contribution degree, which provides a scientific basis for identifying key fissures and effectively managing seepage risks.

[0032] Further, in S4, in the calculation of the alpha-beta rule, alpha is taken as 1 and beta is taken as 0, and the calculation formula is as follows:

[0033]

[0034] In the formula, is the crack outlet flow vector prediction value transmitted by the jth neuron of the l+1th layer to the ith neuron of the lth layer; x i is the output of the ith neuron of the lth layer; w ij is the output of the ith neuron of the lth layer; x i is the weight of the jth neuron of the l+1th layer; w is the contribution degree vector of each component of the permeability vector of the jth neuron of the l+1th layer to the crack outlet flow vector prediction value; u i ∈U (l) is the contribution degree vector of each component of the permeability vector of the jth neuron of the l+1th layer to the crack outlet flow vector prediction value; u i is the neuron belonging to the lth layer.

[0035] Further, the step S5 of determining the main seepage trajectory and the crack where the main seepage trajectory is located based on the particle number and the average contribution degree comprises:

[0036] S5-1, constructing a seepage topology graph based on the crack network model and the simulation calculation result of the seepage in the plurality of cracks in S2;

[0037] S5-2, determining the edge weight of the seepage topology graph according to the particle number and the average contribution degree;

[0038] S5-3, selecting the seepage trajectory with the minimum edge weight in the seepage topology graph, and taking the seepage trajectory as the main seepage trajectory.

[0039] Explanation: The above method can effectively identify the main seepage trajectory by constructing the seepage topology graph and determining the edge weight, and the identified main seepage trajectory helps to screen the target crack and take targeted seepage control measures on the target crack, thereby reducing the seepage risk.

[0040] Further, the edge weight of the seepage topology graph in S5-2 is calculated by the following formula (4):

[0041]

[0042] In the formula, r u is the average contribution degree of the crack corresponding to one of the position nodes u in the seepage topology graph; r v is the average contribution degree of the crack corresponding to one of the position nodes v in the seepage topology graph; N max is the maximum value of the particle number passing through the directed edge in the seepage topology graph; N min is the minimum value of the particle number passing through the directed edge in the seepage topology graph; N uv is the particle number passing through from one of the position nodes u to another position node v in the seepage topology graph; wuv The weight of the edge.

[0043] Furthermore, the seepage parameters mentioned in S5 include real-time seepage flow rate, seepage pressure, and permeability variation coefficient. The standard threshold for the real-time seepage flow rate is 0.8 times the design drainage capacity, the standard threshold for the seepage pressure is 0.7 times the compressive strength of the surrounding sandstone cave, and the standard threshold for the permeability variation coefficient is 0.5.

[0044] Note: The selection of the above seepage parameter values ​​and the setting of standard thresholds can facilitate the assessment and study of the risk of seepage damage in sandstone cave fissures.

[0045] The beneficial effects of this invention are:

[0046] This invention can accurately assess the degree and impact of fissure seepage in grotto rock masses. Based on numerical simulation combined with neural networks, it establishes the relationship between key seepage parameters in the fissure network. Then, using the LRP method, it selects fissures with a significant impact on the grotto from the fissure network for research and risk assessment. This allows for the quantitative expression of seepage conditions and risk levels within the grotto through specific numerical values. Through these research and assessment methods, researchers and grotto conservation workers can better understand the hydrogeological conditions of the grotto rock mass and the mechanism of fissure water transport, more accurately predict potential seepage problems, and thus formulate effective strategies for grotto protection and seepage repair. Attached Figure Description

[0047] Figure 1 This is an illustrative diagram of a multi-branch feedforward neural network according to an embodiment of the present invention. The diagram shows one main branch and two branches.

[0048] Figure 2 This is a schematic diagram of the analysis method flow according to an embodiment of the present invention;

[0049] Figure 3 This is an illustration of the flow velocity along the fracture intersection line in an embodiment of the present invention;

[0050] Figure 4 This is a geographical distribution map of Cave 168 of a certain mountain stone carving in an embodiment of the present invention;

[0051] Figure 5 This is a schematic diagram of the known cracks in Cave 168 of a certain mountain stone carving in an embodiment of the present invention;

[0052] Figure 6 This is the fissure network of Cave 168 of a certain mountain stone carving in this embodiment of the invention;

[0053] Figure 7 This is the distribution of the seepage field in Cave 168 of a certain mountain stone carving in this embodiment of the invention;

[0054] Figure 8 is a schematic diagram of a fissure water seepage source in an embodiment of the present application;

[0055] Figure 9 is a fissure seepage contribution degree of each fissure in Cave 168 of a certain stone carving in an embodiment of the present application;

[0056] Figure 10 is a main seepage track of a fissure network in Cave 168 of a certain stone carving in an embodiment of the present application. DETAILED DESCRIPTION

[0057] In order to further illustrate the manner of carrying out the present application and the effects achieved, the technical solutions of the present application will be described below in conjunction with experiments.

[0058] The present application provides a sandstone grotto fissure water seepage risk assessment method based on numerical simulation, which can accurately identify the source and main seepage track of fissure water seepage disease. The present application trains a multi-branch feedforward neural network by using a large amount of numerical simulation data to learn the corresponding relationship between the fissure permeability vector and the outlet flow vector, and evaluates the seepage contribution degree of each fissure through the LRP algorithm. The edge weight of the seepage topology graph is calculated in combination with the particle number between fissures, and the main seepage track is identified through the shortest path search method. Through the identification of the main seepage track and the above research, the seepage mechanism and risk degree in the fissure can be clearly reflected. The embodiments of the present application will provide certain theoretical support for the accurate treatment of Dazu Stone Carving fissure water seepage disease, and will also provide a beneficial reference for the research of water seepage mechanism of similar grottoes in other regions.

[0059] In the prior art, the sandstone grotto fissure water seepage data is usually obtained through investigation and detection methods. However, the existing exploration technology cannot directly observe the unexposed fissures in the mountain, and needs to generate a random unproven fissure network based on statistical parameters for numerical simulation to obtain a large amount of water seepage data for subsequent research. Therefore, the present application has the problem of low accuracy. Therefore, in the embodiments of the present application, investigation and exploration and numerical simulation are used to obtain the outlet flow data and permeability data of the fissures. See the following embodiments for details:

[0060] Embodiment 1: A sandstone grotto fissure water seepage risk assessment method based on numerical simulation, comprising the following steps:

[0061] First, the outlet flow data and permeability data of the fissures, as well as other hydrogeological data, are obtained through field investigation to facilitate researchers to conduct research;

[0062] 1) Hydrogeological and fissure investigation;

[0063] Geophysical exploration and drilling techniques were used to investigate and statistically analyze lithological distribution, stratigraphic structure, special geological structures, aquifer distribution and types, and fracture characteristics. Geometric elements of fractures (attitude, trace length, and aperture) in rock outcrops were measured. Based on drilling data, core composition, sedimentary sequence, attitude thickness, and variation sequence were analyzed to determine lithological characteristics. Permeability coefficients of different fractures were calculated through pumping and pressure tests. Environmental isotope tracing and hydrochemical methods were applied to track the flow of water in fractured rock masses, analyzing the sources and pathways of water flow in the grottoes. The fracture levels, connectivity, and direct responses to hydrological processes in seepage channels were investigated, along with the hydrological significance of fractures at different levels.

[0064] 2) Investigation of the characteristics of seepage diseases;

[0065] A survey was conducted on historical and existing seepage hazards in the study area. These hazards were classified based on frequency, location, and current condition, and the relationships between different types were analyzed. Continuous monitoring was performed at selected seepage points within the study area, with data recorded at least three times daily at each point. Monitoring equipment was installed at each point, and the flow rate was used to characterize the amount of seepage. Rainfall was also monitored using rain gauges to analyze the response mechanism between seepage volume and rainfall. Chemical parameters such as pH, conductivity, and ion concentration were measured in the seepage water. XRD diffraction analysis was used to investigate the composition of crystalline salts at the locations of seepage hazards, elucidating the chemical characteristics of seepage hazards in the study area.

[0066] 3) Information obtained from the survey; such as Figure 4 As shown, a certain mountain carving is located near the summit of Beishan Mountain in a certain area. The landform type is a flat, eroded, denuded low mountain, with an overall narrow valley and deep hill shape, and an average elevation of 500m. The carving is located on a near-vertical rock wall in the middle and lower part of the north slope of the NEE-trending hill. It is carved on a steep sandstone cliff that extends north and south in a crescent shape, with a height of 3-7m and a top elevation of 506-519m. The valleys in the carving area are not well developed. There are no surface water bodies in the study area. The groundwater consists of pore water from the Quaternary residual slope deposits and fissure water from the bedrock. The pore water is stored in the silty clay layer (mixed with sand and gravel) above the carving wall, while the fissure water is found in the sandstone fissure network. The groundwater in the carving area is entirely dependent on precipitation. Precipitation infiltrates along the soil layer, replenishing the underlying bedrock fissure network. However, due to the obstruction of the mudstone impermeable base beneath the sandstone body, the fissure water drains out on the surface of the carving wall, resulting in fissure seepage damage. Among them, Cave 168 has the most severe and typical water seepage, with dimensions of 3.1m wide and 3.3m high, and is the simulated area in the embodiment of this invention.

[0067] 4) Data: Fracture statistical parameters are based on 590 fracture data obtained from hydrogeological investigation in Beishan area, and the measurement indexes include dip, dip angle, trace length, opening degree, etc. The permeability coefficient of fractured rock mass is calculated through borehole water injection experiment, the experiment covers the thickness of fractured sandstone mass, and water seepage mainly occurs in fractures, so the obtained permeability coefficient can effectively represent the permeability of fractured rock mass.

[0068] In combination with the above, the steps of the embodiments of the present application are as follows:

[0069] Before establishing the fracture network model, a hydrogeological conceptual model is first established for preliminary analysis;

[0070] ①Establish a hydrogeological conceptual model:

[0071] According to the seepage supply range of the stone carving and the regional hydrogeological conditions, the model range is determined, and the hydraulic characteristics are generalized to determine whether the underground water flow is stable flow, Darcy flow and saturated water. In combination with the hydrogeological conditions of the 168 cave area and the field water pressure test results, the fracture roughness is evaluated and the permeability coefficient is determined. The boundary conditions of the generalized conceptual model are determined, including the boundary type (constant water head, constant flow, mixed boundary) and the water head and flow characteristics. The aquifer is generalized to determine the position, top and bottom plate elevation, hydraulic parameters and control fracture distribution of each aquifer, and the hydraulic connection between the aquifers is analyzed. The source and sink items are generalized, the distribution characteristics of the supply (precipitation, lateral runoff of groundwater) and discharge (seepage point, spring outcrop point) are investigated, and the supply source and composition are analyzed by using water chemistry and isotope.

[0072] As shown in Figure 5 , the permeability of the sandstone matrix can be ignored, and the seepage medium is only the fractures of the bedrock. The model range is the fracture network within the 6x20x10m sandstone body of the 168 cave. The upper boundary is generalized as a constant water head boundary, the stone carving wall is a discharge boundary, and the other boundaries are generalized as water-resistant boundaries.

[0073] S1, simulate a plurality of fractures in the sandstone cave based on the fracture network model to obtain a three-dimensional fracture network and fracture data;

[0074] The construction method of the fracture network model comprises:

[0075] S1-1, a plurality of simulated fracture data are randomly generated by using the Monte Carlo method;

[0076] S1-2, the actual fracture data obtained by exploration and the simulated fracture data are used to perform geometric modeling by using the dfnWorks software;

[0077] S1-3, it is judged whether the fracture density in the geometric modeling satisfies P 32If not, the density requirement is met, and the simulated fracture data is verified and corrected by the FRAM method; then the steps of S1-1 and S1-2 are executed in a loop until the fracture density in the geometric modeling meets the P 32 If yes, the three-dimensional fracture network is obtained without verifying and correcting the simulated fracture data. Figure 6

[0078] Specifically, ② a fracture network model is established: existing exploration means cannot observe all the unexposed fractures in the mountain, and the model generates the fracture network by using the deterministic fracture (actual fracture) and the Monte Carlo method simulation (simulated fracture). For the exposed fractures that have been explored, the measured tendency, dip angle, trace length, and opening are generated; for a large number of unmeasured fractures inside the rock mass, the Monte Carlo method is used to generate the fractures by randomly sampling the parameters such as occurrence, radius, and opening.

[0079] Specific data includes: fracture geometric data (including the length, width, and depth of the fracture); spatial position of the fracture (coordinates of the fracture in the three-dimensional space, and spatial attitude information such as the strike, tendency, and dip angle, and opening of the fracture); permeability data of the fracture (permeability of the fracture, and variation range of the permeability); fluid property data (physical properties such as density and viscosity of the fluid, and chemical properties of the fluid); boundary condition and initial condition data (boundary shape of the fracture network, and fluid pressure or flow on the boundary).

[0080] Based on the fracture investigation statistical data, the Monte Carlo method is used to randomly generate the unexposed fractures in the mountain, and corresponding parameters are assigned by grouping. The fracture is initially assumed to be an ellipse, the center of which is located at the origin, and the major and minor axes are distributed along the X and Y axes. The X-axis direction axis length 2a is sampled according to the logarithmic normal distribution, and the Y-axis direction axis length 2b is determined according to the ratio set by the user. In the embodiment of the present application, the axis length ratio is set to 1, that is, the Baecher disc model.

[0081] Subsequently, the ellipse circumference is calculated according to the axis length, and the circumference is equally divided according to the vertex number specified by the user, and the vertex number is set to 6 in the embodiment of the present application; the ellipse circumference coordinates can be described by a parametric equation group (see formula (1-2)), and t can be regarded as the radian turned by a certain point on the circumference from the positive direction of the X axis as the starting point. Set s as the arc length of the ellipse circumference, according to formula (1-3), the coordinates of each vertex can be calculated according to formula (1-2) by integrating from the (a, 0) point counterclockwise with the equally divided arc length as the step. At this point, an elliptical fracture has been approximated to a polygon.

[0082]

[0083] In the formula, a is the half-axis length in the X-axis direction, b is the half-axis length in the Y-axis direction, s is the arc length of the ellipse circumference, and t can be regarded as the starting point with the positive direction of the X axis.​

[0084] The fracture unit normal vector is initially along the Z-axis. A new normal vector is generated based on Fisher distribution sampling. The old normal vector is rotated to the new normal vector using the cross product direction of the old and new normal vectors as the rotation axis. At the same time, the coordinates of the fracture vertices are transformed accordingly so that the fracture orientation conforms to the random sampling results.

[0085] When the center point of a crack is located outside the model area, it may partially extend into the model. Ignoring this situation will result in a lower crack density near the boundary. Therefore, the sampling range of the crack center point coordinates should extend beyond the model area. In this embodiment of the invention, the sampling range of the center point coordinates is set to a cuboid extending 2m outward from the model, and x, y, and z coordinates are randomly generated in a uniform distribution. Then, the crack vertices are translated and transformed to form random cracks generated by the Monte Carlo method.

[0086] To facilitate mesh generation, dfnWorks uses the FRAM method to filter cracks and sets a minimum geometric element length h (0.05m in this embodiment). First, the crack portion extending beyond the model boundary is trimmed to ensure it remains a polygon. Then, starting from a vertex, all vertices are traversed clockwise / counterclockwise. If the distance between adjacent vertices is less than h, the vertex is deleted, and adjacent vertices are connected to form a new edge. This process continues until the traversal ends or the number of vertices is less than 3; in the latter case, the crack is rejected.

[0087] The selected fractures are further subjected to intersection tests. If the length of the intersection line, the distance from the endpoint to the fracture vertex, the spacing between the intersection lines is less than h or the included angle is less than 45°, the fracture is rejected; otherwise, it is accepted and included in the fracture network.

[0088] The above crack generation process will be executed cyclically, with dfnWorks using P 32 Density is the criterion for determining the termination of random generation of a certain group of fractures. 32 Density refers to the ratio of the sum of the areas of the two walls of a fracture to the volume of the simulated region; it is an indicator of fracture density. When P of all groups... 32 When all density values ​​are met, the crack generation cycle terminates.

[0089] The results of the investigation and exploration are as follows Figure 5 As shown, a total of 8 fissures have been identified in Cave 168. These fissures were generated based on measured parameters. Among them, 5 layer fissures cover the model area in the horizontal direction, and the parameters of the remaining 3 fissures are shown in Table 1.

[0090] Table 1. Parameters of three identified steeply dipping fractures in Cave 1168

[0091]

[0092] Unidentified fissures within the grotto were randomly generated using the Monte Carlo method of this invention. Fissure orientations were sampled according to a Fisher distribution, diameters and apertures according to a log-normal distribution, and center point coordinates according to a uniform distribution. The corresponding parameters are shown in Table 2.

[0093] Table 2168 Fracturation Network Parameters

[0094]

[0095] like Figure 6 As shown, P32 represents the density at the termination of each group of random cracks, calculated using a simplified statistical window method, with the P32 of deterministic cracks subtracted to avoid redundant calculations. Since the exposed cracks on the wall surface are already determined, random cracks are set to not intersect with the stone carving wall surface. Ultimately, only the crack clusters connecting the upper boundary and the stone carving wall surface are retained, totaling 129 cracks.

[0096] S2. Based on the three-dimensional fracture network and fracture data, the constrained Delaunay triangulation method and the finite volume method are used in sequence to simulate the seepage in the multiple fractures, and the simulated values ​​of the flow rate vector and permeability vector at the outlet of each fracture are obtained.

[0097] Specifically, the open-source software LaGriT was used, and constrained Delaunay triangulation was applied to divide the 3D fracture network into triangular meshes. The pressure values ​​of the mesh nodes after triangulation were used as the variables to be solved. A Thiessen polygon control volume was established based on the triangular mesh, and the open-source software PFLOTRAN was used to numerically solve the problem using the finite volume method. The governing equations used in PFLOTRAN are shown in formula (1-4), which are essentially Richard equations:

[0098]

[0099] In the formula, η is porosity; s is saturation; η is the molar density of water, in kmol / m³. 3 ; = represents flow velocity, in m / s; Q w Source and sink terms, unit kmol / (m 3 ·s); k is the saturated permeability, in meters. 2 ;k r ρ is a dimensionless coefficient that is a function of relative permeability and saturation; μ is the dynamic viscosity coefficient, in Pa·s; P is the pressure, in Pa; ρ is the density of water, in kg / m³. 3 As mentioned in the conceptual model, the fracture is generalized as a smooth parallel plate. Let b be the fracture aperture, and the permeability k (i.e., the simulated value of the permeability vector) be calculated according to the cubic law, as shown in formula (1-5); the permeability of all grids within the same fracture takes the same value:

[0100]

[0101] Since the numerical solution does not directly obtain the flow velocity distribution, only the normal flow of each edge of the Thiessen polygon control body is obtained, in order to describe the flow field, the flow velocity of each node of the triangular mesh is reconstructed, and four flow velocity vectors are generated on the same node on the intersection line, each of which represents the flow velocity on one side of the intersection line. Flow velocity distribution reconstruction: PFLOTRAN does not directly obtain the flow velocity distribution in the numerical solution, only the normal flow of each edge of the Thiessen polygon control body is obtained, that is, the flow from one node to another adjacent node. In order to describe the flow field, the flow velocity of each node of the triangular mesh must be reconstructed; in the flow velocity reconstruction stage, it is assumed that the flow velocity in each control body is constant, as shown in formula (1-6);

[0102]

[0103] In the formula, G is an nxd matrix, n is the number of control body edges, and each row is the normal area vector of one edge of the control body; is a d-dimensional flow velocity vector, and is an unknown quantity; is an n-dimensional vector, and each component is the flow of one edge of the control body (i.e. the outlet flow vector simulation value of the fracture).

[0104] The flow velocity on the fracture plane is a two-dimensional vector, d=2. The control body is a planar polygon, and the minimum number of edges is 3. n>d, the equation is overdetermined, and the least square method is used to solve the optimal solution See formula (1-7), and the solution is shown in formula (1-8).

[0105]

[0106] The above method for solving flow velocity is applicable to the case where the control body inside the fracture or part of the control body edge is located on the constant head boundary. When part of the control body edge is on the constant flow boundary, the additional constraint of the boundary condition should be considered, and the linear constraint least square method is used to solve the optimal solution See formula (1-10). Normally, nb

[0107]

[0108] In the formula, B is an nbxd matrix, each row is the normal area vector of the control body edge coinciding with the constant flow boundary; is an nb-dimensional vector, and each component is the flow of the control body edge coinciding with the constant flow boundary. The flow velocity information of the node on the intersection line of the fracture is crucial to describe the flow splitting of groundwater at the intersection of the fracture.

[0109] To describe the direction of groundwater movement on both sides of the intersection line, the central node control volume on the intersection line is divided into two planar polygon control volumes, and if two fractures intersect, the intersection line node corresponds to four control volumes. Each control volume reconstructs the node flow rate according to the method described above, so the same node on the intersection line will generate four flow rate vectors, as shown in FIG. 6, each of which represents the flow rate on one side of the fracture intersection line. Figure 3

[0110] The following describes the direction of groundwater movement on both sides of the intersection line:

[0111] Control volume division: the central node control volume on the intersection line is divided into two planar polygon control volumes. Therefore, on the intersection line of two fractures, each node corresponds to four control volumes.

[0112] Flow rate reconstruction: each control volume reconstructs the flow rate at the node according to the method described above.

[0113] Flow rate vector generation: the same node on the intersection line will generate four flow rate vectors, each corresponding to the flow rate on one side of the fracture intersection line, indicating the direction of groundwater movement (see FIG. 6); Figure 3

[0114] Opening setting: considering the roughness of the fracture and the weakening of the permeability of the filling, the fracture opening is set in the way of reducing the mechanical opening to the hydraulic equivalent opening. The mechanical opening of the proven fracture is used, and the mechanical opening of the other fractures is randomly generated according to the logarithmic normal distribution grouping, and the mean and standard deviation of the random fracture opening meet the statistical characteristics of the millimeter-level fracture.

[0115] Average permeability coefficient: according to the average permeability coefficient obtained by the water injection test, the water injection test mainly reflects the lateral permeability of the fractured rock mass. When testing the reduction coefficient, based on the generated fracture network, the back boundary (parallel to the stone wall surface) is set to 1.8×10 5 Pa (corresponding to 3m water head on the bedrock surface), the stone wall surface is 1×10 5 Pa, and the remaining boundaries are water-resistant boundaries. Select multiple reduction coefficients for trial calculation, uniformly reduce all fracture mechanical openings (mechanical opening × reduction coefficient), and calculate the equivalent permeability coefficient in the horizontal direction through numerical simulation, and select the reduction coefficient closest to the measured value. Since the fracture opening is randomly generated, multiple numerical simulations are required to take the average value as the final equivalent permeability coefficient. Tests show that when the number of simulations reaches 20 times, the average equivalent permeability coefficient tends to be stable, so the average value of 20 simulations is taken as the final calculation result.

[0116] Boundary conditions: define the boundary conditions in the form of specified absolute pressure. The boundary of the fracture that does not appear to seep water on the stone wall surface is set as a water-resistant boundary. ​​

[0117] Numerical simulation strategy: Since the mechanical aperture of part of the fractures in the fracture network is randomly generated, the seepage field has randomness under the same boundary conditions. To weaken this influence and explore the general law of the seepage field of the fracture network, the average value of multiple numerical simulations is used as the final result in the embodiment of the application. After each numerical solution and flow velocity reconstruction, the pressure value and flow velocity vector of each node are obtained. The grid is not changed, the node position is fixed, only the pressure head is taken as the arithmetic average, and the final head is calculated according to formula (1-11). The flow velocity only concerns the size, and the flow velocity size of the intersection node of the fractures in a single simulation is calculated according to the arithmetic average of the flow velocity vector modules in each direction, and the final result is the average value of multiple simulation of the flow velocity size of each node.

[0118]

[0119] In the formula, h is the total head, the unit is m; z is the coordinate of node z, the unit is m; n is the number of repetitions; Pi is the pressure value of node i simulation, the unit is Pa. In order to determine the appropriate number of repetitions, the number of repetitions is set to 100, 150 respectively and the average is taken for multiple times.

[0120] The particle transport simulation includes the following (1)-(5);

[0121] (1) Flow velocity interpolation: In order to calculate the motion trajectory of the particle, the flow velocity of any position on the fracture should be clear. The flow velocity of each node on the reconstructed triangular mesh is reconstructed, and the standard barycentric interpolation method is used to interpolate the flow velocity at any position in the triangular mesh to obtain the interpolated flow velocity.

[0122] (2) Time step calculation: The time step needs to be known to update the particle position, and the barycentric interpolation method is also used to calculate the time step. The control body area of each node of the unstructured mesh is different, and this method can obtain a non-uniform time step suitable for the area of each control body and meet the stability criterion.

[0123] (3) New position calculation: After the flow velocity and time step of the present position are calculated by interpolation, the first-order prediction correction method is used to update the new position of the particle. For each time step, the new position is obtained by prediction and correction.

[0124] (4) Selection of downstream fracture at intersection: When the tracer particle moves to the triangular mesh beside the fracture intersection, if the distance between the particle and the intersection is less than the moving distance of the last time step, whether to cross the intersection, if it is detected to pass, the particle stops on the intersection. The fracture intersection is the common edge of the triangular meshes of multiple intersecting fractures, and when the particle moves to the fracture intersection, there will be multiple possible downstream triangular meshes. The probability of the particle entering the downstream mesh is proportional to the flow into the mesh. Then, the position is updated along the fracture where the selected triangular mesh is located.

[0125] (5) Particle number, initial position and outlet: the total trace length of the intersection of each fracture on the upper boundary of the different section fracture network was calculated, and particles were injected at equal intervals on the intersection of each fracture on the upper boundary at intervals of 0.2 m. From the initial position on the upper boundary, the above steps of updating the position of the tracer particle were repeated until the particle flowed out of the stone wall surface outlet.

[0126] Simulation results: based on the fracture network, saturated steady flow simulation and particle transport simulation were performed. When the fourth system overburden was filled with more water, the maximum thickness of the aquifer was 3 m, the upper boundary pressure was set to 1.3 x 10 5 Pa, and the stone wall surface pressure was 1.0 x 10 5 Pa.

[0127] 500 tracer particles were injected at equal intervals on the intersection of each fracture on the upper boundary, and the position of the particle was updated using the first-order prediction correction method. In the simulation, the fracture opening was reduced to the hydraulic equivalent opening according to the average permeability coefficient (8.5 x 10 -5 cm / s) obtained from the borehole water injection experiment. Considering the randomness of the opening, the average value of the simulation water head and flow rate was used as the final result. Tests showed that when the simulation times reached 150 times, the average value tended to be stable.

[0128] S3, based on the fracture outlet flow vector simulation value and the permeability vector simulation value, a fracture flow prediction model is established; the fracture flow prediction model is used to represent the mapping relationship between the fracture permeability vector and the fracture outlet flow vector, and the input of the fracture flow prediction model is the fracture permeability vector (the fracture permeability vector can be a real value measured or a permeability vector simulation value), and the output is the fracture outlet flow vector prediction value;

[0129] S3 includes:

[0130] S3-1, a plurality of sets of fracture outlet flow vector simulation values and permeability vector simulation values are obtained as a data set;

[0131] S3-2, the data set is divided into a training set, a validation set and a test set, and the training set is used to train the constructed fracture flow prediction model; in the training process, the Adam algorithm is selected for optimization, and the mean square error is used as the loss function; the neural network uses a multi-branch feedforward neural network; wherein the number of input layer nodes of the multi-branch feedforward neural network is equal to the total number of fractures in the fracture network, and the number of branches is the number of fractures intersecting the rock wall; the hidden layer of the main network and each branch network in the multi-branch feedforward neural network adopts a three-layer structure, the activation function selects a rectified linear unit, and the batch size is 4;

[0132] S3-3, after the fracture flow prediction model is trained, the validation set is used for verification, the parameters of the fracture flow prediction model are adjusted according to the result of the validation set, the test set is used for testing the fracture flow prediction model and compared with the true value, so as to evaluate the model performance of the fracture flow prediction model.

[0133] Specifically, the training method and result are as follows (1)-(4):

[0134] (1) Data set: the fracture opening is randomly set and not reduced, and other parameters and boundary conditions remain unchanged. 5000-7000 numerical simulations are performed, each simulation generates a sample, including a fracture permeability vector (feature) and a fracture intersection flow vector (label). The data set is divided into a training set of about 3000, a validation set and a test set of 1000 each.

[0135] (2) Multi-branch feedforward neural network structure: a multi-branch feedforward neural network is used to map the fracture permeability to the outlet flow. The network contains one main branch and nine branches, each with two fully connected hidden layers (width same as input layer), input layer width is the number of fractures (129), branch output layer is linear layer, width 1, no activation function.

[0136] (3) Training method: small batch gradient descent method (batch size 4) is used for sample normalization. The loss function is mean square error, the optimizer is ADAM, the initial learning rate is 0.00016, and the learning rate is halved when the validation set loss does not decrease for 50 consecutive times. To prevent overfitting, an early stopping strategy is used, and the training is terminated when the validation set loss does not decrease for 150 consecutive times, and the optimal model is saved.

[0137] (4) Training result: the training is terminated at the 223th cycle, and the minimum validation set loss is 0.00015. The test set evaluation result shows that the average relative error of flow regression is 10%, the regression coefficient of each fracture outlet flow is shown in Table 3, and the average coefficient of determination is 0.86. The neural network regression effect of each fracture outlet flow is good.

[0138] Table 3-168 each fracture outlet flow regression coefficient of determination

[0139]

[0140] S4, the LRP method is used to calculate the average contribution of each fracture permeability vector to the fracture outlet flow vector prediction value; the dfnTrans method in the dfnWorks software is used for particle transport simulation operation of water particles in multiple fractures, to obtain the particle number in the fracture;

[0141] In S4, the method of using LRP each fracture permeability vector to calculate the average contribution of the fracture outlet flow vector includes:

[0142] S4-1, input a fracture permeability vector sample into the fracture flow prediction model to obtain a fracture outlet flow vector prediction value; wherein the fracture permeability vector sample contains the permeability component of each fracture in the fracture network;

[0143] S4-2, for a fracture permeability vector sample, according to the α-β rule, the fracture outlet flow vector prediction value is inversely distributed to the permeability component of each fracture of the fracture permeability vector of the input layer of the fracture flow prediction model layer by layer, to obtain the contribution degree vector c i of each component in the permeability vector to the fracture outlet flow vector prediction value;

[0144] S4-3, select n fracture permeability vector samples, repeat S3-1 and S3-2 n times to obtain n groups of contribution degree vectors; for each fracture i, the average contribution degree in the n fracture permeability vector samples is calculated by formula (1):

[0145]

[0146] In the formula, is the average contribution degree of fracture i, n is the number of fracture permeability vector samples, c i is the contribution degree vector of fracture i.

[0147] In the calculation of the α-β rule, α=1 and β=0 are taken, and the calculation formula is as follows formula (3):

[0148]

[0149] In the formula, is the fracture outlet flow vector prediction value transmitted by the jth neuron of the l+1 layer to the ith neuron of the l layer; x i is the output of the ith neuron of the l layer; w ij is the weight of x i to the jth neuron of the l+1 layer; is the contribution degree vector of each component in the permeability vector of the jth neuron of the l+1 layer to the fracture outlet flow vector prediction value; u i ∈U (l) is u i belongs to the neuron of the l layer.

[0150] It should be understood that the LRP (Layer-wise Relevance Propagation) algorithm is a method of eXplainable AI (explainable artificial intelligence) for explaining the neural network calculation process. The algorithm can calculate the contribution of the input layer elements to the output layer. In the embodiments of the present application, the alpha-beta rule (alpha=1, beta=0) is used for information back propagation.

[0151] In this calculation, the LRP method is used to back propagate the fracture outlet flow information, quantitatively calculate the contribution of the fracture permeability to the outlet flow, and evaluate the importance of individual fractures in the fracture network. The method includes the following steps:

[0152] Accumulate data sets: use the fracture seepage network model, keep the boundary conditions unchanged. Randomly sample the fracture opening according to the lognormal distribution, perform seepage numerical simulation, and calculate the fracture flow of the outlet boundary. The fracture permeability vector (input) and the outlet flow vector (output) form a data pair. Repeat the simulation and accumulate enough data.

[0153] Train the neural network: divide the data set into training set, validation set and test set. Train the neural network to learn the mapping relationship between the fracture permeability vector and the outlet flow vector, which is a generalization model of complex numerical simulation.

[0154] LRP back propagation calculation: for each sample in the data set, input the fracture permeability vector for forward calculation to get the outlet flow prediction value. The LRP algorithm is used for back propagation to calculate the contribution of each fracture. The contribution of all samples is arithmetically averaged to obtain the final contribution of each fracture to the outlet flow.

[0155] S5, according to the particle number and the average contribution, determine the main seepage trajectory and the fracture where the main seepage trajectory is located, and take the fracture where the main seepage trajectory is located and / or the fracture with the largest average contribution of the fracture permeability vector to the fracture outlet flow vector prediction value as the target fracture; when the value of the seepage parameter in the target fracture exceeds the standard threshold value, it is considered that the selected range of sandstone grotto has high risk and preventive measures need to be taken, otherwise, it is considered that the sandstone grotto has low risk and preventive measures do not need to be taken.

[0156] The determination of the main seepage trajectory and the fracture where the main seepage trajectory is located based on the particle number and the average contribution includes:

[0157] S5-1, based on the fracture network model and the seepage simulation results of the plurality of fractures in S2, a seepage topology graph is constructed;

[0158] S5-2, according to the particle number and the average contribution, the edge weight of the seepage topology graph is determined;

[0159] S5-3, selecting a seepage trace with the minimum edge weight in the seepage topology graph as the main seepage track.

[0160] The edge weight of the seepage topology graph is calculated by formula (4) as follows:

[0161]

[0162] In the formula, r u is the average contribution degree of a fissure corresponding to one of the position nodes u in the seepage topology graph; r v is the average contribution degree of a fissure corresponding to one of the position nodes v in the seepage topology graph; N max is the maximum value of the number of particles passing through the directed edge in the seepage topology graph; N min is the minimum value of the number of particles passing through the directed edge in the seepage topology graph; N uv is the number of particles passing through from one of the position nodes u to another position node v in the seepage topology graph; w uv is the weight of the edge.

[0163] The seepage parameters include real-time seepage flow, seepage pressure and seepage rate variation coefficient, the standard threshold value of the real-time seepage flow is 0.8 times the design drainage capacity, the standard threshold value of the seepage pressure is 0.7 times the compressive strength of the sandstone cave surrounding rock, and the standard threshold value of the seepage rate variation coefficient is 0.5.

[0164] Specifically, the main seepage track search result: through 150 times of particle transport simulation, a large number of traces are accumulated to establish a seepage topology graph of the overall fissure network. The nodes in the figure represent fissures, the directed edges represent particle motion directions, and the edge weights are determined according to the number of particles and the seepage contribution of the fissures. The NetworkX is used to identify the main seepage track by searching the shortest simple path.

[0165] Results and analysis: the average water head distribution of 150 times of seepage simulation of the fissure network of 168 cave is shown in Figure 7 (a), the overall water head presents a gradually decreasing trend from top to bottom and from back to front, which is consistent with the direction of water flow. The local seepage layer presents a relatively large flow rate, which indicates that it plays a leading role in seepage, which is consistent with the serious water damage of the cave. The average flow rate distribution is shown in Figure 7 (b), the groundwater flow rate generally presents a trend of becoming larger as it is closer to the rock wall surface, which is consistent with the change of the equipotential line density. The overall seepage is slow, and the groundwater flow rate in most fissures is lower than 0.0097 m / s, and the minimum flow rate is 7.4×10 -11m / s. Given the numerical simulation boundary conditions are conducive to seepage, the flow rate in the actual fracture network should be lower, and groundwater in some fractures may be stagnant. Accordingly, it can be reasonably speculated that when the bedrock top surface is recharged, groundwater near the wall side will quickly drain; groundwater in the deeper fractures inside the mountain will be stored in the mountain, thus forming slow and continuous recharge to the fractures on the wall side.

[0166] Fracture outlet seepage source;

[0167] In Cave 168, No. 5 layer fracture and J96, J97 fractures represent three typical fracture outlets, respectively. Particle trajectory distribution reflects groundwater migration trajectory. Based on particle transport simulation, only the particle trajectories flowing out of the specified outlet are counted, the number of particles passing through the fracture is recorded, and the arithmetic mean of 150 simulations is taken as the final result. Delete the average number of particles less than 1, and the remaining fracture network constitutes the catchment network corresponding to the fracture outlet. Particle motion trajectory presents obvious characteristics. In J96, J97 and J98 steep fractures, particle trajectories in deep rock mass are mainly vertical downward from the top, while fractures near the wall are inclined downward to the outlet. The particles in the layer fracture move from inside to outside, and the trajectory is generally perpendicular to the stone wall. The particle trajectory is dense at the intersection of large steep fractures. The distribution of particle trajectories is consistent with the distribution of water head and flow rate, indicating that its density can intuitively reflect the seepage trajectory and flow of groundwater. The seepage channels of different fracture outlets are different, and the seepage sources are different: J96 fracture is mainly supplied by the top layer fracture; J97 fracture is derived from its upper fracture and No. 4 layer fracture channel; No. 5 layer fracture is relatively scattered, and is greatly affected by large steep fractures that cut through multiple layers of fractures.

[0168] Fracture seepage contribution evaluation;

[0169] The seepage contribution of each fracture calculated by LRP algorithm is linearly mapped to [0, 10], [10, 1], and the results are shown in Figure 9 J97 fracture has the largest contribution, followed by No. 5 layer fracture. The seepage contribution of steep fractures is basically proportional to their development scale, and the contribution of small random fractures is generally low. The seepage contribution of layer fractures changes regularly with height, and the contribution of No. 2, 4 and 5 layer fractures directly discharged on the wall is significantly higher than that of other layers, and increases with the decrease of height, while the contribution of No. 1 and 3 layer fractures shows the opposite trend.

[0170] Main seepage trajectory evaluation;

[0171] The information of the shortest five main seepage trajectories of the overall fracture network in Cave 168 is shown in Table 4, and the specific location is shown in Figure 10The fastest seepage track is supplied by the bedrock top surface, and directly discharged to the wall surface through J97 and J96 fissures, in which J97 is more advantageous, and the length of the two tracks is much smaller than other channels. The second is the seepage track formed by the combination of J97 and No. 4 and No. 5 layer fissures. Although the discharge of J97 and J96 fissures alone is obviously advantageous, the seepage capacity between them is not as good as that of the combination of layer fissures. Overall, the main seepage track of the 168 cave section is composed of fissures directly connecting the bedrock top surface and the stone carving wall surface, and steep fissures vertically penetrating multiple layer fissures.

[0172] Table 4 - Main seepage track of 168 cave fissure network

[0173]

[0174]

[0175] In the embodiment of the present application, the distribution of the seepage field and the distribution of the seepage track of the 168 cave fissure network are obtained by numerical simulation.

[0176] The results show that the groundwater in the 168 cave mainly migrates obliquely downward along the bedrock top surface and is discharged in the grotto area. The groundwater near the wall surface is quickly discharged, while the groundwater in the deep fissures of the mountain body is stored and slowly supplies the wall surface fissures. The seepage channels of different fissure outlets are different, and the sources of seepage water are different: the J96 fissure is mainly supplied by the top layer fissure; the J97 fissure is derived from the upper fissure and the No. 4 layer fissure channel; the No. 5 layer fissure is relatively dispersed and is greatly affected by the large steep fissure cutting through multiple layers. The contribution of different fissures to the seepage network is obviously different, among which the J97 fissure contributes the most, followed by the No. 5 layer fissure. The seepage contribution of steep fissures is roughly proportional to their development scale, while small-scale random fissures contribute less. In the fissure network of the cave section, J96 and J97 fissures are the fastest seepage track, and the advantage is far superior to other tracks. The risk detection results of J96 and J97 fissures show that their seepage indexes do not exceed the standard, indicating that the overall risk of the grotto is low, and no protective measures need to be taken.

Claims

1. A method for assessing the risk of water seepage through fissures in sandstone grottoes based on numerical simulation, characterized in that, Includes the following steps: S1. Simulate multiple fissures in the sandstone cave based on the fissure network model to obtain a three-dimensional fissure network and fissure data; S2. Based on the three-dimensional fracture network and fracture data, the constrained Delaunay triangulation method and the finite volume method are used in sequence to simulate the seepage in the multiple fractures, and the simulated values ​​of the flow rate vector and permeability vector at the outlet of each fracture are obtained. S3. Based on the simulated values ​​of the fracture outlet flow vector and the simulated values ​​of the permeability vector, a fracture flow prediction model is established. The fracture flow prediction model is used to characterize the mapping relationship between the fracture permeability vector and the fracture outlet flow vector. The input of the fracture flow prediction model is the fracture permeability vector, and the output is the predicted value of the fracture outlet flow vector. S4. The LRP method is used to calculate the average contribution of each fracture permeability vector to the predicted value of the fracture outlet flow vector; then the dfnTrans method in dfnWorks software is used to perform particle transport simulation calculations on the water particles seeping in multiple fractures to obtain the number of particles in the fractures. S5. Based on the number of particles and the average contribution, determine the main seepage trajectory and the fracture where the main seepage trajectory is located. The fracture where the main seepage trajectory is located and / or the fracture with the largest average contribution of the fracture permeability vector to the predicted value of the fracture outlet flow vector are taken as the target fracture. When the value of the seepage parameter in the target fracture exceeds the standard threshold, the sandstone cave within the selected range is considered to have a high risk and preventive measures need to be taken. Otherwise, the sandstone cave is considered to have a low risk and no preventive measures are required.

2. The method for assessing the risk of seepage through fissures in sandstone grottoes based on numerical simulation as described in claim 1, characterized in that, The method for constructing the fracture network model described in S1 includes: S1-1. Multiple simulated fracture data are randomly generated using the Monte Carlo method; S1-2. Using the actual fracture data obtained from the exploration and the simulated fracture data, perform geometric modeling using dfnWorks software; S1-3. Determine whether the fracture density in the geometric modeling satisfies P. 32 Density requirement; if not, use the FRAM method to verify the simulated fracture data; then repeat steps S1-1 and S1-2 until the fracture density in the geometric modeling meets P. 32 Stop after reaching the required density, and you will obtain the three-dimensional fracture network and the verified simulated fracture data; if so, there is no need to verify or correct the simulated fracture data, and you will obtain the three-dimensional fracture network.

3. The method for assessing the risk of water seepage in sandstone grottoes based on numerical simulation as described in claim 2, characterized in that, The fracture data refers to actual fracture data and / or simulated fracture data / verified simulated fracture data. The fracture data includes fracture geometry data, fracture spatial location data, fracture permeability data, fluid property data, boundary conditions, and initial conditions data.

4. The method for assessing the risk of water seepage in sandstone grottoes based on numerical simulation as described in claim 1, characterized in that, The establishment of the fracture flow prediction model described in S3 includes: S3-1. Obtain multiple sets of simulated values ​​of fracture outlet flow rate and permeability vector as a dataset; S3-2. The dataset is divided into a training set, a validation set, and a test set. The constructed fracture flow prediction model is trained using the training set. During the training process, the Adam algorithm is selected for optimization, and the mean squared error is used as the loss function. The neural network is a multi-branch feedforward neural network. S3-3. After the fracture flow prediction model is trained, it is validated using a validation set. The parameters of the fracture flow prediction model are adjusted based on the results of the validation set. Then, the fracture flow prediction model is tested using a test set and compared with the true values ​​to evaluate the model performance.

5. The method for assessing the risk of seepage through fissures in sandstone grottoes based on numerical simulation as described in claim 1, characterized in that, The method of sequentially using constrained Delaunay triangulation and finite volume method to simulate the seepage in the multiple fractures and obtain the flow rate vector and permeability vector at each fracture outlet includes: firstly, using constrained Delaunay triangulation to divide the mesh in the three-dimensional fracture network into triangular meshes; then, setting the pressure value of the mesh node as a variable in the solution; then, establishing a Thiessen polygon control volume based on the triangular mesh; and finally, using the finite volume method to solve the Thiessen polygon control volume to obtain the simulated values ​​of the flow rate vector and permeability vector at each fracture outlet.

6. The method for assessing the risk of seepage through fissures in sandstone grottoes based on numerical simulation as described in claim 1, characterized in that, In S4, the method of averaging the contribution of each LRP fracture permeability vector to the predicted fracture outlet flow vector includes: S4-1. Input a fracture permeability vector sample into the fracture flow prediction model to obtain the predicted value of the fracture outlet flow vector; wherein, a fracture permeability vector sample contains the permeability component of each fracture in the fracture network; S4-2. For a single fracture permeability vector sample, following the α-β rule, the predicted fracture outlet flow rate vector is inversely distributed layer by layer to the permeability component of each fracture in the fracture permeability vector of the input layer of the fracture flow rate prediction model, thus obtaining the contribution vector c of each component in the permeability vector to the predicted fracture outlet flow rate vector. i ; S4-3. Select n fracture permeability vector samples, repeat S3-1 and S3-2 n times to obtain n sets of contribution vectors; for each fracture i, the average contribution in the n fracture permeability vector samples is calculated by formula (1): In the formula, Let n be the average contribution of fracture i, n be the number of fracture permeability vector samples, and c be the average contribution of fracture i. i This is the contribution vector of crack i.

7. The method for assessing the risk of seepage through fissures in sandstone grottoes based on numerical simulation as described in claim 6, characterized in that, In S4, when calculating the α-β rule, taking α = 1 and β = 0, the calculation formula is as follows (3): In the formula, x is the predicted value of the slit exit flow vector transmitted from the j-th neuron in layer l+1 to the i-th neuron in layer l; i The output of the i-th neuron in layer l; w ij For linear computation of fully connected layers, x i The weights of the j-th neuron in layer l+1; u is the contribution vector of each component in the permeability vector of the (l+1)th neuron j to the predicted value of the fracture outlet flow vector; i ∈U (l) For u i Neurons belonging to layer l.

8. The method for assessing the risk of water seepage in sandstone grottoes based on numerical simulation as described in claim 2, characterized in that, The determination of the main seepage trajectory and the fracture where the main seepage trajectory is located based on the number of particles and the average contribution rate, as described in S5, includes: S5-1. Based on the fracture network model and the seepage simulation results in the multiple fractures in S2, construct a seepage topology map; S5-2. Determine the edge weights of the seepage topology graph based on the number of particles and the average contribution. S5-3. Select the seepage trajectory line that minimizes the edge weight in the seepage topology graph, and use the seepage trajectory line as the main seepage trajectory.

9. The method for assessing the risk of water seepage in sandstone grottoes based on numerical simulation as described in claim 8, characterized in that, The edge weights of the seepage topology graph in S5-2 are calculated using the following formula (4): In the formula, r u The average contribution of a fracture to a node u in the seepage topology diagram. r v N represents the average contribution of a fracture to a node v in the seepage topology diagram. max N represents the maximum number of particles passing through a directed edge in the seepage topology graph. min N represents the minimum number of particles passing through a directed edge in the seepage topology graph. uv w represents the number of particles that pass through from one location node u to another location node v in the seepage topology graph. uv The weight of the edge.

10. The method for assessing the risk of water seepage in sandstone grottoes based on numerical simulation as described in claim 1, characterized in that, The seepage parameters mentioned in S5 include real-time seepage flow rate, seepage pressure, and permeability variation coefficient. The standard threshold for the real-time seepage flow rate is 0.8 times the design drainage capacity, the standard threshold for the seepage pressure is 0.7 times the compressive strength of the surrounding sandstone cave, and the standard threshold for the permeability variation coefficient is 0.5.

Images