Numerical simulation method for estimating equivalent seepage of three-dimensional fracture network based on topological property

Abstract

The invention discloses a numerical simulation method for estimating equivalent seepage of a three-dimensional fracture network based on topological properties, which comprises the following steps: firstly, constructing a small three-dimensional fracture network sample set on the basis of related data of a fracture network in a research area, then mapping the sample fracture network into an equivalent two-dimensional undirected graph by using a graph theory, solving equivalent permeability, and calculating the equivalent seepage of the fracture network. Forming a training sample set; then, a maximum information index method is adopted to mine fracture network equivalent permeability main control factors, and the control relation between fracture network parameters and an underground flow field is determined; and finally, constructing a substitution model between the main control factors and the equivalent permeability based on a BPNN machine learning model, and realizing conversion from the target discrete fracture network to an equivalent medium. According to the method, a simulation scheme which reflects topological connectivity and is high in prediction precision and low in calculation cost is provided for numerical simulation of the fractured stratum underground flow field under the field scale.

Description

Technical Field

[0001] The invention relates to a numerical simulation method for a fracture network flow field in an underground stratum, and in particular to a numerical simulation method for estimating equivalent seepage of a three-dimensional fracture network based on topological properties. Background Art

[0002] Fracture networks are widely present in natural rock matrices, which can provide advantageous paths and water barriers for fluid flow and affect the distribution of underground flow fields. Therefore, for well-consolidated bedrock, obtaining the equivalent permeability of the fracture network in the underground strata relative to the lower pore permeability is the key to seepage simulation. However, most popular fracture equivalent permeability models rely on observed geometric features and ignore the topological structure of the fracture network. In fact, the fracture distribution has obvious anisotropy, and the topological structure has a strong control effect on seepage. The discrete fracture model can accurately describe the connectivity in the fracture network and the flow of fluid in the fractures, but when the fracture network is complex, the computational cost of the model is high. The equivalent permeable medium model is more efficient, but it does not consider the topological connectivity, thereby greatly reducing the complexity of the model, and the calculation results are highly biased. Summary of the invention

[0003] Purpose of the invention: The purpose of the present invention is to provide a numerical simulation method for estimating the equivalent seepage of a three-dimensional fracture network based on topological properties with high prediction accuracy and low computational cost for fractured strata at the site scale.

[0004] Technical solution: The numerical simulation method for estimating equivalent seepage of a three-dimensional fracture network based on topological properties of the present invention comprises:

[0005] (1) Collect basic data on fracture networks and basic hydrogeological data, and construct a target discrete three-dimensional fracture network based on the basic data on fracture networks;

[0006] (2) Construct a conceptual model of the underground flow field of the fracture network based on the target discrete three-dimensional fracture network and basic hydrogeological data;

[0007] (3) According to the geological conditions, hydrogeological conditions and the scale of the target area, the coverage of the fracture network parameters is determined, and then a number of small three-dimensional fracture networks are established by randomly generating fracture occurrences using the Monte Carlo method to form the initial training sample set;

[0008] (4) Use graph theory to map the sample fracture network into an equivalent two-dimensional undirected graph and solve the equivalent permeability of the sample fracture network;

[0009] (5) The maximum information index method is used to mine the main control parameters of equivalent permeability from several fracture network parameters, and a sample set with the main control parameters as input and equivalent permeability as output is constructed, and the input and output variables are normalized;

[0010] (6) Based on the BPNN machine learning model, a prediction model for the equivalent permeability of a three-dimensional fracture network is constructed. The target discrete fracture network is gridded, and the three-dimensional matrix composed of the main control parameters is used as input. The three-dimensional matrix composed of the equivalent permeability, i.e., the equivalent medium, is used as output. The hyperparameters of the model are adjusted to achieve high-precision input-output simulation. Finally, a verified numerical model for simulating the underground flow field of the fracture network is obtained.

[0011] Furthermore, in step (1), the basic data of the fracture network include fracture occurrence distribution characteristics, fracture diameter distribution characteristics, and fracture density distribution characteristics.

[0012] Furthermore, the basic hydrogeological data include structural geological maps and hydrogeological maps.

[0013] Furthermore, in step (3), the occurrence parameters of the small three-dimensional fracture network are as follows:

[0014] The crack diameter follows a power law distribution and is expressed as:

[0015] f(l)=βl -α

[0016] Where, f(l) is the probability density function of the crack diameter distribution; β is the minimum value l of the crack diameter. min and the maximum value l max The normalized index obtained; α is the power law index;

[0017] The crack diameter is expressed as:

[0018]

[0019] Where R is a random number;

[0020] The occurrence of fractures follows the Fisher distribution:

[0021]

[0022] Among them, θ is the inclination; θ0 is the average inclination; φ is the dip; φ0 is the average dip; κ is the discrete degree parameter; sinh is the hyperbolic sine function;

[0023] The crack opening λ is between 0.1 mm and 1 mm, which conforms to the uniform random distribution; it is assumed that the cracks in each crack network have the same opening.

[0024] Furthermore, in step (4), the sample fracture network is mapped into an equivalent two-dimensional undirected graph using graph theory, including:

[0025] Considering the fractures as nodes V and the intersections between fractures as edges E, the three-dimensional fracture network is converted into a two-dimensional undirected graph G(V, E) as follows:

[0026] The three-dimensional fracture network F consists of N fractures f, which can be expressed as:

[0027] F={f i}, i = 1, ..., N

[0028] Define a bijective mapping Φ:

[0029] Φ:f i →v i

[0030] If the two cracks f i With f j intersect, Then there is an edge in E connecting the corresponding vertices:

[0031]

[0032] Among them, (v i , v j )∈E represents vertex v i and v j The edge between ij , edge ij The weight w ij Set to vertex v i and v j The length of the intersection between:

[0033] w ij =Lengthf i ∩f j

[0034] Similarly, considering the flow direction, the inlet plane x0 is regarded as the source node s, and the outlet plane x L Considered as the target node t, when the crack intersects the entry and exit boundaries:

[0035]

[0036]

[0037] Among them, e si is the edge between the source node and the non-source node, e it is the edge between the non-source sink node and the target node; thus, the three-dimensional fracture network is converted into a two-dimensional undirected graph;

[0038] When the three-dimensional fracture network is sparse, the two-dimensional undirected graph will be divided into several parts. Only the fracture subnetwork corresponding to the connected subgraph containing the source node and the target node can effectively connect the inlet plane and the outlet plane to induce percolation. Therefore, the connected subgraph without the source node and the target node can be ignored. That is, the two-dimensional undirected graph G can be pruned to the union of all connected subgraphs containing the source node and the target node except the edge nodes, that is, the subset G' of the two-dimensional undirected graph G. G' has a unique corresponding subnetwork F' of the three-dimensional fracture network F.

[0039] Furthermore, the equivalent permeability of the sample fracture network is solved, including: solving the equivalent permeability of the fracture network flow direction on the equivalent subgraph G' of the subnetwork F', and the solution is based on the following assumptions:

[0040] i. The rock matrix is ​​incompressible and impermeable;

[0041] ii. The fluid is incompressible;

[0042] iii. The two planes of the crack are parallel and smooth;

[0043] iv. Fluid flow follows Darcy’s law and exhibits laminar behavior.

[0044] Furthermore, assuming that there are M nodes in the two-dimensional undirected graph G, the nodes represent the intersections between cracks, and node i is adjacent to node j. According to the cubic law, the total flow Q between node i and node j is j It is expressed as:

[0045]

[0046] Where k is the equivalent permeability; g is the gravitational acceleration; μ is the dynamic viscosity coefficient of water flow; I is the hydraulic gradient; λ ij is the opening between crack i and crack j; b is the crack width; H i is the upstream head, H j is the downstream head, L ij is the Euclidean distance between node i and node j;

[0047] According to the law of conservation of mass, the flow at node i is expressed as:

[0048]

[0049] The regional hydraulic head distribution of the sample fracture network is expressed as:

[0050]

[0051] Where ΔH ij is the upstream and downstream head difference; w ijis the equivalent hydraulic conductivity between nodes i and j, and is also the weight of the edge connecting nodes i and j. When there is no connection between nodes i and j, w ij The value is 0;

[0052] From this, we can obtain the adjacency matrix A of the two-dimensional undirected graph G and the vertex degree matrix D, which are defined as

[0053] The equivalent permeability k of the entire sample fracture network area is expressed as:

[0054]

[0055] Among them, L zone is the domain length of the sample fracture network, Q si is the inflow, H s is the inlet head, H t It is the water head.

[0056] Furthermore, in step (5), the fracture network parameters include topological and geometric parameters, and the topological parameters include:

[0057] degree centrality, which represents the number of rifts that intersect a rift;

[0058] Eigenvector centrality is an iterative index of degree centrality, indicating the importance of the fracture in the entire fracture network;

[0059] clustering coefficient, representing the diversity of fracture network pathways;

[0060] Betweenness centrality, which represents the number of cracks in the shortest path of the network, plays an important role in controlling percolation;

[0061] The maximum flow is the feasible flow with the maximum flow in a network with capacity, and is an important topological feature of percolation;

[0062] Anisotropy index f A , the expression is:

[0063]

[0064] Where N is the number of cracks, λ is the crack opening, i is the unit vector of the flow direction, and n is the crack normal vector;

[0065] Geometric parameters include the number of cracks, crack opening, and crack area per unit volume P. 32 and unit body size.

[0066] Furthermore, in step (5), the maximum information index method is used to mine the main control parameters, including:

[0067] Assume that the data set has Γ data points, draw the scatter plot of the two attributes in two-dimensional space, and divide it into m×n grids along the horizontal and vertical axes, m×n<B, B=Γ 0.6 ; The grid with coordinates (x0, y0) contains γ data points, and the frequency P(x, y) = γ / Γ; Based on this, the maximum grid resolution of the normalized mutual information is calculated as the MIC metric:

[0068]

[0069] Among them, X and Y are two variables whose correlation needs to be calculated;

[0070] The value range of MIC is |0, 1|. When the two variables are independent of each other, MIC approaches 0. On the contrary, when the two variables are strongly correlated, MIC approaches 1.

[0071] Furthermore, in step (6), the BPNN machine learning model has a three-layer structure, including an input layer, a hidden layer, and an output layer, and the Tan-sigmoid and nonlinear LM algorithms are set as transfer functions and training functions, respectively.

[0072] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: the present invention determines the main control parameters of the equivalent permeability of the fracture network through the maximum information index method, and constructs an estimation model of the equivalent permeability of the three-dimensional fracture network based on the BPNN machine learning model, which comprehensively improves the high computational cost of the discrete fracture model and the shortcomings of the equivalent permeability medium model that does not consider the topological connectivity. The fracture network can be efficiently and reliably converted into the corresponding equivalent medium, and provides a simulation scheme with high prediction accuracy and low computational cost that reflects the topological connectivity for the numerical simulation of the underground flow field of the fracture strata at the site scale. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 It is a flowchart of a numerical simulation method for estimating equivalent seepage of a three-dimensional fracture network based on topological properties provided by an embodiment of the present invention;

[0074] Figure 2 It is a schematic diagram of the structure of the BPNN machine learning model in an embodiment of the present invention;

[0075] Figure 3 is a schematic diagram of the correlation analysis results between the fracture network parameters and the equivalent permeability in an embodiment of the present invention;

[0076] Figure 4 (a) is a target discrete crack network diagram in an embodiment of the present invention, Figure 4 (b) is an estimated equivalent permeable medium diagram in an embodiment of the present invention;

[0077] Figure 5(a) is a top view of the target discrete fracture network in the vertical direction ranging from 0 to 15 meters in an embodiment of the present invention, Figure 5 (b) is a top view of the estimated equivalent permeable medium in the vertical direction ranging from 0 to 15 meters in an embodiment of the present invention;

[0078] Figure 6 (a) is a top view of the target discrete fracture network underground flow field in the vertical direction ranging from 0 to 15 meters in an embodiment of the present invention, Figure 6 (b) is a top view of the estimated equivalent permeable medium underground flow field in the vertical direction ranging from 0 to 15 meters in an embodiment of the present invention;

[0079] Figure 7 It is a top view of a target discrete fracture network, an estimated equivalent permeable medium, and an estimated equivalent permeable medium underground flow field streamline in a vertical direction ranging from 0 to 15 meters in an embodiment of the present invention. DETAILED DESCRIPTION

[0080] The present invention will be further described below in conjunction with the accompanying drawings.

[0081] like Figure 1 As shown, an embodiment of the present invention provides a numerical simulation method for estimating equivalent seepage of a three-dimensional fracture network based on topological properties, comprising the following steps:

[0082] (1) Collect basic data on fracture networks and basic data on hydrogeology, and construct a target discrete three-dimensional fracture network based on the basic data on fracture networks. The basic data on fracture networks include the distribution characteristics of fracture occurrence, fracture diameter, and fracture density, and the basic data on hydrogeology include structural geological maps and hydrogeological maps.

[0083] (2) A conceptual model of the underground flow field of the fracture network is constructed based on the target discrete three-dimensional fracture network and basic hydrogeological data.

[0084] (3) According to the geological conditions, hydrogeological conditions and the scale of the target area, the coverage of the fracture network parameters is determined, and then a number of small three-dimensional fracture networks are established by randomly generating fracture occurrences using the Monte Carlo method to form the initial training sample set;

[0085] The occurrence parameters of the small 3D fracture network are as follows:

[0086] The crack diameter follows a power law distribution and is expressed as:

[0087] f(l)=βl -α (1)

[0088] Where f(l) is the probability density function of the crack diameter distribution; β is the minimum value l of the crack diameter. min and the maximum value l maxThe normalized exponent obtained; α is the power law exponent, set to 3.0;

[0089] The crack diameter is expressed as:

[0090]

[0091] Where R is a random number;

[0092] The occurrence of fractures follows the Fisher distribution:

[0093]

[0094] Among them, θ is the inclination; θ0 is the average inclination; φ is the dip; φ0 is the average dip; κ is the discrete degree parameter, which takes the value as 0; sinh is the hyperbolic sine function;

[0095] The fracture opening λ is between 0.1 mm and 1 mm, which is consistent with the uniform random distribution. The fracture opening differs from the fracture size by more than two orders of magnitude, and the effect of the size change of the fracture opening on the seepage is so small that it can be ignored. Therefore, the present invention does not consider the change of the fracture aperture, and assumes that the fractures in each fracture network have the same opening.

[0096] Based on the above parameter ranges and value selection methods, the Monte Carlo method is used to generate the initial training sample set.

[0097] (4) Use graph theory to map the sample fracture network into an equivalent two-dimensional undirected graph and solve the equivalent permeability of the sample fracture network;

[0098] (4-1) Use graph theory to map the sample fracture network into an equivalent two-dimensional undirected graph

[0099] Considering the fractures as nodes V and the intersections between fractures as edges E, the three-dimensional fracture network is converted into a two-dimensional undirected graph G(V, E) as follows:

[0100] The three-dimensional fracture network F consists of N fractures f, which can be expressed as:

[0101] F={f i}, i = 1, ..., N

[0102] Define a bijective mapping Φ:

[0103] Φ:f i →v i

[0104] If the two cracks f i With f j intersect, Then there is an edge in E connecting the corresponding vertices:

[0105]

[0106] Among them, (v i , v j )∈E represents vertex v i and v j The edge between ij , edge ij The weight w ij Set to vertex v i and v j The length of the intersection between:

[0107] w ij =Lengthf i ∩f j

[0108] Similarly, considering the flow direction, the inlet plane x0 is regarded as the source node s, and the outlet plane x L is regarded as the target node t. When the crack intersects the entry and exit boundaries,

[0109]

[0110]

[0111] Among them, e si is the edge between the source node and the non-source node, e it is the edge between the non-source node and the target node. So far, the three-dimensional fracture network is converted into a two-dimensional undirected graph.

[0112] When the three-dimensional fracture network is sparse, the two-dimensional undirected graph will be divided into several parts. Only the fracture subnetwork corresponding to the connected subgraph containing the source node and the target node can effectively connect the inlet plane and the outlet plane to induce percolation. Therefore, the connected subgraph without the source node and the target node can be ignored. That is, the two-dimensional undirected graph G can be pruned to the union of all connected subgraphs containing the source node and the target node except the edge nodes, that is, the subset G' of the two-dimensional undirected graph G. G' has a unique corresponding subnetwork F' of the three-dimensional fracture network F.

[0113] (4-2) Solving the equivalent permeability of the sample fracture network

[0114] The equivalent permeability of the fracture network flow direction is solved on the equivalent subgraph G' of the subnetwork F', and the solution is based on the following assumptions:

[0115] i. The rock matrix is ​​incompressible and impermeable;

[0116] ii. The fluid is incompressible;

[0117] iii. The two planes of the crack are parallel and smooth;

[0118] iv. Fluid flow follows Darcy’s law and exhibits laminar behavior.

[0119] Assume that there are M nodes in the two-dimensional undirected graph G, and the nodes represent the intersections between cracks. Node i is adjacent to node j. According to the cubic law, the total flow Q between node i and node j is j It is expressed as:

[0120]

[0121] Where k is the equivalent permeability; g is the gravitational acceleration; μ is the dynamic viscosity coefficient of water flow; I is the hydraulic gradient; λ ij is the opening between crack i and crack j; b is the crack width; H i is the upstream head, H j is the downstream head, L ij is the Euclidean distance between node i and node j;

[0122] According to the law of conservation of mass, the flow at node i is expressed as:

[0123]

[0124] The regional hydraulic head distribution of the sample fracture network is expressed as:

[0125]

[0126] Where ΔH ij is the upstream and downstream head difference; w ij is the equivalent hydraulic conductivity between nodes i and j, and is also the weight of the edge connecting nodes i and j. When there is no connection between nodes i and j, w ij The value is 0;

[0127] From this, we can obtain the adjacency matrix A of the two-dimensional undirected graph G and the vertex degree matrix D, which are defined as

[0128] The equivalent permeability k of the entire sample fracture network area is expressed as:

[0129]

[0130] Among them, L zone is the domain length of the sample fracture network, Q si is the inflow, H s is the inlet head, H t It is the water head.

[0131] (5) The maximum information index method is used to mine the main control parameters (or main control parameters) of equivalent permeability from several fracture network parameters, and a sample set with the main control parameters as input and equivalent permeability as output is constructed, and the input and output variables are normalized;

[0132] (5-1) Fracture network parameters include topological and geometric parameters

[0133] Several important topological parameters are listed below:

[0134] degree centrality, which represents the number of rifts that intersect a rift;

[0135] Eigenvector centrality is an iterative index of degree centrality, indicating the importance of the fracture in the entire fracture network;

[0136] clustering coefficient, representing the diversity of fracture network pathways;

[0137] Betweenness centrality, which represents the number of cracks in the shortest path of the network, plays an important role in controlling percolation;

[0138] The maximum flow is the feasible flow with the maximum flow in a network with capacity, and is an important topological feature of percolation;

[0139] The above topological parameters only represent the properties of a single fracture, so the present invention takes the global average value and variance for them to reflect the seepage capacity and heterogeneity of the fracture network.

[0140] The anisotropy index fA is expressed as:

[0141]

[0142] Where i is the unit vector of the flow direction and n is the crack normal vector.

[0143] Geometric parameters include the number of cracks, crack opening, and crack area per unit volume P. 32 and unit body size.

[0144] (5-2) Use the maximum information index method to mine the main control parameters of the above topological parameters and a small number of geometric parameters

[0145] The principle of the maximum information index method is as follows:

[0146] Assume that the data set has Γ data points, draw the scatter plot of the two attributes in two-dimensional space, and divide it into m×n grids along the horizontal and vertical axes, m×n<B, B=Γ 0.6 ; The grid with coordinates (x0, y0) contains γ data points, and the frequency P(x, y) = γ / Γ; Based on this, the normalized maximum grid resolution of mutual information is calculated as the MIC metric, that is:

[0147]

[0148] Among them, X and Y are two variables whose correlation needs to be calculated.

[0149] The value range of MIC is |0, 1|. When the two variables are independent of each other, MIC approaches 0. On the contrary, when the two variables are strongly correlated, MIC approaches 1.

[0150] (6) A prediction model for the equivalent permeability of a three-dimensional fracture network is constructed based on the BPNN machine learning model. The target discrete fracture network in step (1) is gridded, and the three-dimensional matrix composed of the main control parameters is used as input, and the three-dimensional matrix composed of the equivalent permeability, i.e., the equivalent medium, is used as output. The hyperparameters of the model are adjusted to achieve high-precision input-output simulation, and finally a verified numerical model for simulating the underground flow field of the fracture network is obtained.

[0151] like Figure 2 As shown, the BPNN machine learning model includes an input layer, a hidden layer and an output layer, each layer contains multiple neurons, and the neurons of adjacent layers are interconnected through weighted connections, w and V are connection coefficients. The input layer serves as the data inlet, accepts input x, and obtains data h after passing it to the hidden layer. The hidden layer calculates the data h and generates output y. After minimizing the loss through the training function, the connection coefficient in the network is adjusted by back propagation until the loss reaches the expected level or the training reaches the expected number of times. The BPNN machine learning model established by the present invention has three layers, and the Tan-sigmoid and nonlinear LM (Levenberg-Marquardt) algorithms are set as the transfer function and the training function respectively.

[0152] The BPNN machine learning model training uses MATLAB software, and all sample analysis, spatial discretization and numerical simulation in the present invention are performed based on MATLAB software.

[0153] The numerical simulation method for estimating equivalent seepage of a three-dimensional fracture network based on topological properties provided in an embodiment of the present invention adopts a BPNN machine learning model to estimate the equivalent permeable medium containing topological properties in the stratum, which greatly reduces the computational load of solving the underground flow field model by a discrete fracture network model, and provides an efficient and reliable simulation solution for fracture medium seepage simulation under actual site conditions, and has broad application prospects.

[0154] A specific example is given below.

[0155] This example conducts numerical simulation to estimate equivalent seepage of three-dimensional fracture networks under actual site conditions, taking into account the topological properties of the fracture network. It mainly includes two parts: 1) mining the main controlling factors of equivalent permeability of fracture networks based on the maximum information index method; 2) an equivalent permeability estimation model based on the BPNN machine learning model.

[0156] 1) Mining the main controlling factors of equivalent permeability of fracture networks based on the maximum information index method

[0157] In this example, four groups of typical occurrences were selected to prepare the sample set, and the specific discrete fracture network parameters are shown in Table 1. In scenario 1, the fracture diameter distribution is relatively discrete; in scenario 2, the fracture diameter distribution is relatively uniform; in scenario 3, there is a group of conjugate fractures; in scenario 4, the fracture orientation conforms to a uniform random distribution, to emphasize the universality of the present invention.

[0158] The unit cell side length of the simulation area is the maximum crack diameter l max As the size of the unit body changes, the fracture diameter changes accordingly. Different simulation scenarios are set accordingly to solve the equivalent permeability of the fracture network in the main flow direction and establish a database. In this example, a large amount of data is obtained by controlling the number of fractures between 5 and 4000 and the unit body size between 0.55 and 100 meters. After removing outliers using the k-means method, a database with a capacity of 5044 is obtained.

[0159] Table 1 Discrete fracture network parameters of sample set

[0160]

[0161] Among them, Δφ0 is the difference range between the dip angle and the average dip angle, and Δθ0 is the difference range between the inclination and the average inclination.

[0162] The maximum information index method was used to conduct a correlation analysis on the fracture network parameters and equivalent permeability in the sample set. The specific analysis results are as follows: Figure 3 As shown, Figure 3 Where lnk represents the logarithm of equivalent permeability (the order of magnitude of equivalent permeability is usually small, at least 1E-20, so its logarithm is used in analysis and prediction). Fracture networks of different occurrences have the same main control parameters, and the change of unit size has no obvious effect on the correlation of parameters. In different scenarios, there is a good correlation between geometric parameters, and the correlation of topological parameters is quite different. Fracture aperture, P 32 The variance of the clustering coefficient is consistent with the characteristics of strong correlation with the equivalent permeability and weak correlation with other parameters, and is determined as the main control parameter of the estimation model to achieve streamlined model operation and accurate results.

[0163] 2) Equivalent permeability estimation model based on BPNN machine learning model

[0164] This example takes the synthetic fracture network as an example, and the fracture network parameters are shown in Table 2. The BPNN machine learning model used in this example has a maximum training number of 1000 times, a learning rate of 0.01, and a training target minimum error of 0.00001. The number of input layer nodes is 3, and the number of hidden layer nodes is 5. The sample set area is a cube with a side length of 15 meters. The number of discrete fractures ranges from 16 to 336. According to the fracture survey results, 1050 small three-dimensional fracture networks are generated. The data set is divided into a training set and a test set in a ratio of 8:2. The training set is used to train the BPNN machine learning model, and the test set is used to verify the performance of the BPNN machine learning model. The target area is a cube with a side length of 75 meters, containing 1000 discrete fractures. Based on the maximum diameter of the fracture, the target area is discretized into 125 sub-areas of 5×5×5 along the x, y, and z axes.

[0165] Table 2 Example model fracture network parameters

[0166]

[0167] By adjusting the number of training set samples, R 2 , RMSE and MAPE to obtain the optimal training set capacity of the estimation model, and the results are shown in Table 3. In all cases, the estimation effect of the model is very good, and the prediction results have a good fit with the calculated values, which shows that the method of the present invention has high accuracy. When the data set capacity is 1050, the estimation model performs best, and the fitting effect is slightly improved compared with the estimation results in other scenarios.

[0168] Table 3 Reliability index of the estimated model

[0169]

[0170] like Figures 4 to 6 As shown in Figure 2, (a) is the discrete fracture model and (b) is the estimated equivalent medium model. Figure 4 It can be seen that the estimation model converts the three-dimensional discrete fracture model of the site into an equivalent medium model. Figure 5 The discrete fracture model in the range of 0 to 15 m on the Z axis is described, which contains a total of 187 discrete fractures. The equivalent permeability is consistent with the distribution trend of the discrete fractures. At position (1,2), the mesh has a higher fracture density and exhibits a higher equivalent permeability. In contrast, the meshes at positions (1,3), (1,4), (5,2), and (5,3) show lower fracture density and lower equivalent permeability. Figure 6The steady-state seepage fields of the discrete fracture network model and the estimated equivalent medium model are shown respectively. Although the spatial discretization of the equivalent medium masks the strong heterogeneity of the fractures to a certain extent, it effectively captures the flow resistance and diversion phenomena caused by the uncertain fracture distribution at y = 15m and y = 45m, and the pressure front presents a distinct concave and convex shape. Figure 7 The flow field in the equivalent medium is shown. By comparing the discrete fractures and the equivalent medium, it is shown that the streamline trend is consistent with the fracture distribution. The fracture distribution is highly uneven near the inlet and the streamlines are concentrated around the grid (1,2) position. In the middle of the seepage path, the fracture distribution is more uniform and the streamlines are more dispersed. The increase in the unevenness of the fracture distribution at the outlet causes the streamlines to move to the area with denser fractures, indicating the important role of the fracture distribution for the seepage path, and further proves the good ability of the simulation method of the present invention to restore the topological connectivity of the fracture network.

[0171] In general, the trained BPNN machine learning model can effectively predict the equivalent medium distribution of the three-dimensional fracture network, more realistically reflect the heterogeneity and topological connectivity of rock fractures, and achieve a high prediction accuracy for the spatial distribution of underground seepage. Therefore, this prediction model can be used for the three-dimensional visualization of the underground flow field of the three-dimensional fracture network under actual site conditions.

Claims

1. A numerical simulation method for estimating equivalent seepage in a three-dimensional fracture network based on topological properties, characterized in that: include: (1) Collect basic data on fracture networks and basic hydrogeological data, and construct a target discrete three-dimensional fracture network based on the basic data on fracture networks; (2) Construct a conceptual model of the underground flow field of the fracture network based on the target discrete three-dimensional fracture network and basic hydrogeological data; (3) According to the geological conditions, hydrogeological conditions and the scale of the target area, the coverage of the fracture network parameters is determined, and then a number of small three-dimensional fracture networks are established by randomly generating fracture occurrences using the Monte Carlo method to form the initial training sample set; (4) Use graph theory to map the sample fracture network into an equivalent two-dimensional undirected graph and solve the equivalent permeability of the sample fracture network; (5) The maximum information index method is used to mine the main control parameters of equivalent permeability from several fracture network parameters, and a sample set with the main control parameters as input and equivalent permeability as output is constructed, and the input and output variables are normalized; (6) Based on the BPNN machine learning model, a prediction model for the equivalent permeability of a three-dimensional fracture network is constructed. The target discrete fracture network is gridded, and the three-dimensional matrix composed of the main control parameters is used as input. The three-dimensional matrix composed of the equivalent permeability, i.e., the equivalent medium, is used as output. The hyperparameters of the model are adjusted to achieve high-precision input-output simulation. Finally, a verified numerical model for simulating the underground flow field of the fracture network is obtained.

2. The numerical simulation method according to claim 1, characterized in that: In step (1), the basic data of the fracture network include fracture occurrence distribution characteristics, fracture diameter distribution characteristics, and fracture density distribution characteristics.

3. The numerical simulation method according to claim 2, characterized in that: Basic hydrogeological data include structural geological maps and hydrogeological maps.

4. The numerical simulation method according to claim 1, characterized in that: In step (3), the occurrence parameters of the small three-dimensional fracture network are as follows: The crack diameter l follows a power law distribution, expressed as: f(l)＝βl -α Where, f(l) is the probability density function of the crack diameter distribution; β is the minimum value l of the crack diameter. min and the maximum value l max The normalized index obtained; α is the power law index; The crack diameter l is expressed as: Where R is a random number; The occurrence of fractures follows the Fisher distribution: Among them, θ is the inclination; θ0 is the average inclination; φ is the dip; φ0 is the average dip; κ is the discrete degree parameter; sinh is the hyperbolic sine function; The crack opening λ is between 0.1 mm and 1 mm, which conforms to the uniform random distribution; it is assumed that the cracks in each crack network have the same opening.

5. The numerical simulation method according to claim 1, characterized in that: In step (4), the sample fracture network is mapped into an equivalent two-dimensional undirected graph using graph theory, including: Considering the fractures as nodes V and the intersections between fractures as edges E, the three-dimensional fracture network is converted into a two-dimensional undirected graph G(V,E) as follows: The three-dimensional fracture network F consists of N fractures f, which can be expressed as: F={f i },i=1,……,N Define a bijective mapping Φ: F:f i →v i If the two cracks f i With f j intersect, Then there is an edge in E connecting the corresponding vertices: Among them, (v i ,v j )∈E represents vertex v i and v j The edge between ij , edge ij The weight w ij Set to vertex v i and v j The length of the intersection between: w ij =Lengthf i ∩f j Similarly, considering the flow direction, the inlet plane x0 is regarded as the source node s, and the outlet plane x L Considered as the target node t, when the crack intersects the entry and exit boundaries: Among them, e si is the edge between the source node and the non-source sink node, e it is the edge between the non-source sink node and the target node; thus, the three-dimensional fracture network is converted into a two-dimensional undirected graph; When the three-dimensional fracture network is sparse, the two-dimensional undirected graph will be divided into several parts. Only the fracture subnetwork corresponding to the connected subgraph containing the source node and the target node can effectively connect the inlet plane and the outlet plane to induce percolation. Therefore, the connected subgraph without the source node and the target node can be ignored. That is, the two-dimensional undirected graph G can be pruned to the union of all connected subgraphs containing the source node and the target node except the edge nodes, that is, the subset G' of the two-dimensional undirected graph G. G' has a unique corresponding subnetwork F' of the three-dimensional fracture network F.

6. The numerical simulation method according to claim 5, characterized in that: Solve the equivalent permeability of the sample fracture network, including: solve the equivalent permeability of the fracture network flow direction on the equivalent subgraph G' of the subnetwork F', and the solution is based on the following assumptions: i. The rock matrix is ​​incompressible and impermeable; ii. The fluid is incompressible; iii. The two planes of the crack are parallel and smooth; iv. Fluid flow follows Darcy’s law and exhibits laminar behavior.

7. The numerical simulation method according to claim 6, characterized in that: Assume that there are M nodes in the two-dimensional undirected graph G, and the nodes represent the intersections between cracks. Node i is adjacent to node j. According to the cubic law, the total flow Q between node i and node j is ij It is expressed as: Where k is the equivalent permeability; g is the gravitational acceleration; μ is the dynamic viscosity coefficient of water flow; I is the hydraulic gradient; λ ij is the opening between crack i and crack j; b is the crack width; H i is the upstream head, H j is the downstream head, L ij is the Euclidean distance between node i and node j; According to the law of conservation of mass, the flow at node i is expressed as: The regional hydraulic head distribution of the sample fracture network is expressed as: Where ΔH ij is the upstream and downstream head difference; w ij is the equivalent hydraulic conductivity between node i and node j, and is also the weight of the edge connecting node i and node j. When there is no connection between node i and node j, w ij The value is 0; From this, we can obtain the adjacency matrix A of the two-dimensional undirected graph G and the vertex degree matrix D, which are defined as The equivalent permeability k of the entire sample fracture network area is expressed as: Among them, L zone is the domain length of the sample fracture network, Q si is the inflow, H s is the inlet head, H t It is the water head.

8. The numerical simulation method according to claim 1, characterized in that: In step (5), the fracture network parameters include topological and geometric parameters, and the topological parameters include: degree centrality, which represents the number of rifts that intersect a rift; Eigenvector centrality is an iterative index of degree centrality, indicating the importance of the fracture in the entire fracture network; clustering coefficient, representing the diversity of fracture network pathways; Betweenness centrality, which represents the number of cracks in the shortest path of the network, plays an important role in controlling percolation; The maximum flow is the feasible flow with the maximum flow in a network with capacity, and is an important topological feature of percolation; Anisotropy index f A , the expression is: Where N is the number of cracks, λ is the crack opening, i is the unit vector of the flow direction, and n is the crack normal vector; Geometric parameters include the number of cracks, crack opening, and crack area per unit volume P. 32 and unit body size.

9. The numerical simulation method according to claim 8, characterized in that: In step (5), the maximum information index method is used to mine the main control parameters, including: Assume that the data set has Γ data points, draw the scatter plot of the two attributes in two-dimensional space, and divide it into m×n grids along the horizontal and vertical axes, m×n<B, B=Γ 0.6 ; The grid with coordinates (x0, y0) contains γ data points, and the frequency P(x, y) = γ / Γ; Based on this, the normalized maximum grid resolution of mutual information is calculated as the MIC metric: Among them, X and Y are two variables whose correlation needs to be calculated; The value range of MIC is [0,1]. When the two variables are independent of each other, MIC approaches 0. On the contrary, when the two variables are strongly correlated, MIC approaches 1.

10. The numerical simulation method according to claim 1, characterized in that: In step (6), the BPNN machine learning model has a three-layer structure, including an input layer, a hidden layer, and an output layer, and the Tan-sigmoid and nonlinear LM algorithms are set as transfer functions and training functions, respectively.

Images