A modeling method and prediction method for a reservoir rock fracture network expansion prediction model based on neural network

Abstract

The present invention provides a modeling method for a reservoir rock fracture network expansion prediction model based on a neural network, which specifically includes the following steps: step 1, obtaining original geological parameters, fracturing design parameters and natural fracture distribution parameters of each type of oil reservoir; generating multiple original geological parameter random values according to the original geological parameters, and using fracture combination parameters to generate multiple fracture combination parameter random values, and combining the original geological parameter random values and the fracture combination parameter random values to form a fracture network expansion simulation data set; step 2, obtaining fractal parameters of a fractal structure, combining the fracture network expansion simulation data set and the fractal parameters, constructing a feature data set based on the fractal structure, and dividing the feature data set into a training set, a validation set and a test set; step 3, establishing a prediction model based on a neural network; step 4, inputting the training set into the prediction model for training, and after the training is completed, inputting the validation set into the prediction model for verification, and inputting a test machine into the prediction model for testing.

Description

Technical Field

[0001] The present invention relates to the field of shale reservoir fracture prediction, and in particular to a modeling method and a prediction method of a reservoir rock fracture network expansion prediction model based on a neural network. Background Art

[0002] Hydraulic fracturing is a key technology for shale reservoir reconstruction. Fracture monitoring results show that shale reservoirs are highly developed with sedimentary bedding and natural fractures. Volumetric fracturing in shale reservoirs creates a complex fracture network that provides effective pathways for oil and gas flow. Rapidly and accurately predicting the expansion morphology of the complex fracture network formed by shale fracturing is crucial for the efficient development of shale oil and gas reservoirs.

[0003] Current research on fracture network expansion in shale reservoirs faces the following challenges: 1) Most existing studies are based on finite element analysis (FEM), which relies on accurate physical models and the solution of partial differential equations. This approach requires meshing and iterative solutions, resulting in a computationally complex and time-consuming solution process. Furthermore, FEM requires adjustments to the model and boundary conditions for different shale samples, requiring re-modeling and re-solving, resulting in relatively poor transferability. 2) Predicting the expansion of shale reservoir fracture networks involves the morphology of the fracture network. These rock fracture networks exhibit self-similar and self-affine fractal properties. Fractal structures embody the spatial distribution and multi-scale characteristics of fractures, and influence their connectivity and permeability, significantly impacting the paths and patterns of fracture expansion. However, traditional prediction methods often ignore the fractal structure when characterizing fracture network morphology, resulting in inaccurate predictions and impacting the assessment of fracture expansion behavior, permeability, and hydraulic fracturing effectiveness. Summary of the Invention

[0004] The main purpose of the present invention is to provide a modeling method and a prediction method of a reservoir rock fracture network expansion prediction model based on a neural network, so as to solve the above technical problems.

[0005] To achieve the above objectives, the present invention adopts a technical solution: a modeling method for a reservoir rock fracture network expansion prediction model based on a neural network, which specifically includes the following steps:

[0006] Step 1: Obtain original geological parameters and fracture combination parameters of each reservoir rock; generate multiple random values of original geological parameters based on the original geological parameters, and generate multiple random values of fracture combination parameters based on the original geological parameters and fracture combination parameters, and combine the random values of original geological parameters and fracture combination parameters to form a fracture network expansion simulation data set;

[0007] Step 2: Obtain fractal parameters of the fractal structure, combine the fracture network extension simulation data set and the fractal parameters, construct a feature data set based on the fractal structure, and divide the feature data set based on the fractal structure into a training set, a validation set, and a test set;

[0008] Step 3: Establish a prediction model based on neural network;

[0009] Step 4: Input the training set into the prediction model for training. After the training is completed, input the validation set into the prediction model for verification and input the test set into the prediction model for testing.

[0010] Preferably, generating a plurality of random values of fracture combination parameters according to the original geological parameters and the fracture combination parameters comprises the following steps:

[0011] Step 11: Setting initial parameters and boundary conditions: The initial parameters include geological parameters, shale Young's modulus and Poisson's ratio, initial fracture length and initial fracture inclination; the geological parameters include porosity and permeability; the boundary conditions include stress field and displacement conditions, and the initial conditions and initial parameters are input into the displacement discontinuity method model;

[0012] Step 12: Simulate the propagation of the crack by solving the displacement on the crack surface using the displacement discontinuity method: Calculate the stress and displacement on the crack surface and update the crack morphology at each time step. Through iterative calculation, obtain the crack morphology at different moments in the predetermined time series. Combine the obtained multiple groups of crack morphologies to obtain random values of the crack combination parameters. The crack morphology includes crack length and crack inclination.

[0013] Step 13: Combine the generated original geological parameter random values and the fracture combination parameter random values.

[0014] Preferably, before executing step 2, the fracture network extension simulation dataset is processed as follows:

[0015] Step 14: Map the combination obtained in step 13 from high dimension to low dimensional feature space through principal component analysis;

[0016] Step 15: construct an advanced feature dataset based on the data obtained in step 14. The advanced feature dataset is the fracture network extension simulation dataset in step 2.

[0017] Preferably, step 12 specifically includes the following steps:

[0018] Step 121: Based on the initial parameters, boundary conditions, and basic principles of elasticity, the stress boundary integral equation is used to solve the displacement discontinuity on the crack surface:

[0019]

[0020] in, is the initial ground stress field, σ ij (x) is the total stress field after considering the effect of cracks, Γ is the crack surface, ∫ Γ G ijkl (x, y) is the Green's function, i, j, k, and l represent the components of the tensor in the Green's function, i and j represent the displacement or stress component direction of the observation point, that is, the shear stress acting on the i-th surface along the j-th direction, k and l represent the direction of the force or stress applied by the source point, x represents the corresponding position of the calculated stress point, y represents the point on the integration path, Δu is the displacement on the crack surface, and the process of solving Δu is as follows:

[0021] The crack surface is discretized into N sub-elements, and the displacement on the sub-element is Δu q , through the above discretization, for each crack element p, the stress boundary integral equation can be transformed into a linear algebraic equation:

[0022]

[0023] Among them, p and q are both discrete crack units, p represents the observation point or the specified discrete position for calculating stress; q represents the subunit index on the crack surface, which is used to traverse all subunits. is the total stress field at discrete point p, Δu q is the displacement on the qth subunit, It represents the Green's function at the discrete point p, describing the influence of unit q on the stress at p. By summing all N crack subunits q, it represents the cumulative stress contribution of each subunit at p.

[0024] When solving formula (2) in the iterative process, firstly calculate Δu q Make a guess as the initial value and calculate the And the corresponding energy norm E(u), the energy norm can be defined as:

[0025]

[0026] in, is the strain tensor, which is expressed by the displacement Δu q Substitute into the geometric equation to calculate:

[0027]

[0028] : represents the double inner product operation between tensors, Ω is the volume domain of shale, ρ is the density of shale volume force, which affects the overall force balance and energy distribution;

[0029] By adjusting Δu qThe value of is to minimize the energy norm E(u), so as to gradually make Δu q and Converge to the true value;

[0030] Step 122: According to the obtained Δu q and Update the shape of the crack. The specific steps are as follows:

[0031] Step 1221, calculate the stress intensity factor, which is divided into the cracking intensity factor K I and slip intensity factor K II , respectively expressed as:

[0032]

[0033]

[0034] Among them, σ and τ are normal stress and shear stress respectively, and normal stress and shear stress are The weight, is the crack length, Y, Y' is a geometric factor related to the crack shape and loading conditions;

[0035] Step 1222: Calculate the capacity release rate G, which is calculated using the following formula:

[0036]

[0037] Where E' is the equivalent elastic modulus. For plane stress, E'=E, and for plane strain, ν is Poisson's ratio;

[0038] Step 1223: According to the crack extension criterion, whether the crack extends depends on whether G reaches or exceeds the fracture toughness G of the material. c , if G ≥ G c , the crack will expand;

[0039] Step 1224: When the crack is determined to have expanded, the crack morphology is updated; otherwise, this step is not performed. Updating the crack morphology includes the following steps:

[0040] 1) Add a certain amount of stress to get the new stress field σ i ' j (x), and let G = Gc;

[0041] 2) The new stress field σ i ' j (x) and the new G are substituted into formulas (5), (6), and (7) to obtain the new crack length a new and the crack inclination angle θ,

[0042]

[0043] 3) According to the obtained a new and θ to update the fracture morphology, wherein the fracture morphology includes fracture length and fracture inclination, and the fracture morphology obtained each time is saved;

[0044] Step 123, return to step 121, calculate the next sub-unit until all sub-units are calculated, and save the multiple groups of fracture morphologies obtained by simulation as data files, which are random values of fracture combination parameters.

[0045] Preferably, step 14 specifically includes the following steps:

[0046] Step 141: Standardize the data sample so that its mean is 0 and its variance is 1;

[0047] Step 142: Calculate the covariance matrix C:

[0048]

[0049] Find the eigenvectors and eigenvalues of the covariance matrix C, that is, solve the equation Cv = λv, where λ is the eigenvalue vector, v is the corresponding eigenvector, n represents the number of data samples, x i is the feature vector of the i-th data sample, is the average vector of all data samples in the data sample set;

[0050] Step 143: Select the largest m eigenvalues in λ: λ1≥λ2≥…≥λ m The corresponding eigenvectors constitute the transformation matrix W=[v1,v2,…v m ],

[0051] Let Z = XW,

[0052] X is the data matrix of the data sample, each row of X is a data sample, each column is a feature, Z is the principal component matrix, each row of Z is a sample, and each column is a principal component.

[0053] Preferably, step 15 adopts one of the following two methods:

[0054] 1) Standardization: Scale the features to the same scale, where z is one of the columns of the principal component matrix Z, z' is the standardized result, μ is the mean of the column, and σ is the standard deviation of the column.

[0055]

[0056] 2) Normalization: Scale the features to the interval [0,1], where z is one column of the principal component matrix Z and z' is the normalized result.

[0057]

[0058] After this step, the feature dataset is obtained:

[0059]

[0060] in It is a set of feature data. is the number of eigenvectors contained in the set, F n is the nth eigenvector Its dimensional components Represents the high-level features after principal component decomposition and standardization / normalization, corresponding to the i-th standardized / normalized principal component z' n ,i ,i∈(1,m).

[0061] Preferably, in step 2, a feature data set based on a fractal structure is constructed and the fractal dimension is calculated by a box counting method, and the steps are as follows:

[0062] Step 21. Cover the crack structure with grids of different sizes and count the number of boxes N(∈) covering the crack for each grid size, where ∈ is the side length of the grid. Plot the double logarithm logN(ε) and log(1 / ε). The slope of the graph is the fractal dimension D, which can be calculated using the formula:

[0063]

[0064] Step 21: Calculate the fractal length distribution. First, collect the crack length L data and plot the corresponding probability distribution graph and log P(L) versus log L. P(L) is the ratio of the number of cracks of a certain length to the total number of cracks. The slope of the length distribution graph and the width distribution graph is calculated to be the length distribution index α.

[0065] Step 22: adding the obtained fractal dimension D and the length distribution index α as additional features to the fracture network extension simulation data set obtained in step 1 to obtain a feature data set based on the fractal structure;

[0066] Step 23: Obtain the crack network extension morphology distribution data at multiple prediction moments as corresponding prediction labels, add the prediction labels to the fractal structure-based feature data set obtained in step 22 to obtain a final feature data set, save the final feature data set in a data frame format (DataFrame), and export it as a CSV file;

[0067] Step 24: The final feature dataset obtained in step 23 is randomly divided into a training set, a validation set, and a test set in a ratio of 6:2:2. The training set is used for model training, the validation set is used for verification during model training, and the test set is used to finally evaluate the generalization ability of the model.

[0068] Preferably, the step three specifically includes the following steps:

[0069] Step 31: Set the first layer as a convolutional layer, which is used to capture the relationship between geological parameters or the pattern and distribution characteristics of fractures in space;

[0070] Step 32: The second layer is provided with a flattening layer, which is used to flatten the multi-dimensional tensor data output by the convolutional layer into a one-dimensional vector;

[0071] Step 33: The third layer is provided with a multilayer perceptron, which is used to extract global features from the one-dimensional vector output by the flattening layer and capture the relationship between geological parameters, fracture combination parameters and analytical structural data. Two fully connected layers are provided, and the number of neurons in each layer can be adjusted according to data complexity and model requirements. The multilayer perceptron composed of two fully connected layers is used to extract its global features and capture the relationship between geological parameter features, fracture combination parameters and analytical structural data.

[0072] Step 34. Select the activation function ReLU.

[0073] Compared with the prior art, the present invention has the following beneficial effects:

[0074] The prediction model established by the present invention has high prediction accuracy, strong generalization ability, and is applicable to different types of oil reservoirs; the prediction model only takes a few seconds to predict the expansion of the fracture network, overcoming the defect that the traditional fracture network expansion simulation method takes several hours or even longer to simulate the fracture network expansion process, greatly reducing the time cost and improving the efficiency of predicting the expansion of the shale reservoir fracturing fracture network. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 is a flow chart of the modeling method according to the present invention. DETAILED DESCRIPTION

[0076] The following description is intended to disclose the present invention so that those skilled in the art can implement the present invention. The preferred embodiments described below are merely examples, and those skilled in the art may conceive of other obvious variations.

[0077] like Figure 1 As shown, a modeling method for a shale reservoir fracture network expansion prediction model based on a neural network specifically includes the following steps:

[0078] Step 1: Obtain the Required Data: Analyze and organize field data to obtain the original geological parameters and fracture combination parameters for each reservoir type. Based on the original geological parameters, geostatistical methods are used to generate multiple random values for the original geological parameters. Simulate the fracture network expansion using a displacement discontinuity model to generate multiple random values for the fracture combination parameters. The random values for the original geological parameters and the fracture combination parameters are combined to form a fracture network expansion simulation dataset. The original geological parameters include porosity distribution, permeability distribution, and original geostress field, while the fracture combination parameters include fracture length and inclination.

[0079] Specifically, geostatistical methods such as Kriging or Monte Carlo simulation can effectively simulate the spatial distribution and variation of geological parameters. In the present invention, Monte Carlo simulation is used to generate random values of the original geological parameters. The original geological parameters are assumed to conform to a normal distribution, with a set mean and standard deviation. For example, the mean of porosity is 0.2 and the standard deviation is 0.05; the mean of permeability is 100 and the standard deviation is 20. Monte Carlo simulation can generate random values for these geological parameters.

[0080] The initial parameters and boundary conditions of the displacement discontinuity method model are initialized based on the original geological parameters and fracture combination parameters. The fracture combination parameters are output in each simulation. The simulation process is repeated to generate multiple fracture combination parameter samples. The fracture combination parameters reflect the geometric characteristics and spatial distribution of fractures and are important inputs for simulating the propagation morphology of fracture networks.

[0081] Specifically, the displacement discontinuity method is used to simulate the propagation of the fracture network to generate random values of multiple fracture combination parameters, including the following steps:

[0082] Step 11: Setting initial parameters and boundary conditions: The initial parameters include geological parameters, shale Young's modulus and Poisson's ratio, initial fracture length and initial fracture inclination; geological parameters include porosity, permeability, etc.; the boundary conditions include stress field and displacement conditions. The initial parameters and boundary conditions are input into the displacement discontinuity method model to initialize the fracture distribution in the simulation area;

[0083] Step 12: Use the displacement discontinuity method to solve the discontinuous displacement on the crack surface to simulate the expansion of the crack: calculate the stress and displacement on the crack surface and update the crack morphology at each time step. The crack morphology includes the crack length and crack direction. Through iterative calculation, the distribution results of the crack network expansion morphology at different times in the predetermined time series can be obtained. The specific steps include the following:

[0084] Step 121: Calculate the discontinuous displacement and total stress field on the crack surface based on the initial parameters and boundary conditions. According to the basic principles of elasticity, the stress boundary integral equation can be used to solve the displacement discontinuity on the crack surface:

[0085]

[0086] in, is the initial ground stress field, σ ij (x) is the total stress field after considering the effect of cracks, Γ is the crack surface, ∫ Γ G ijkl (x, y) is the Green's function, which can represent the displacement of a point in the elastic body in response to the unit force of another point. i, j, k, and l represent the components of the tensor in the Green's function. i and j represent the displacement or stress component direction of the observation point, that is, the shear stress acting on the i-th surface along the j-th direction. k and l represent the direction of the force or stress applied by the source point. x represents the corresponding position of the calculated stress point, and y represents the point on the integral path. The Green's function describes the response of the observation point in the i and j directions when the force is applied in the direction of the source point k and l. Δu is the displacement on the crack surface. Since σ ij (x) is also an unknown quantity that needs to be solved, so when solving the actual situation of discontinuous displacement Δu k (y) is solved indirectly through discretization and iterative methods, as follows:

[0087] The crack surface is discretized into N sub-units. Assume that the discontinuous displacement on each sub-unit is a constant or a linear function, and the displacement on the sub-unit is Δu q Through the above discretization, for each crack element p, the stress boundary integral equation can be transformed into a linear algebraic equation:

[0088]

[0089] Among them, p and q are both discrete sub-units, p represents the specific discrete position of the observation point or calculated stress; q represents the sub-unit index on the crack surface, which is used to traverse all sub-units. is the total stress field at subunit p, Δu q is the displacement on the qth subunit, represents the Green's function at the discrete point p, which describes the influence of subunit q on the stress at p. By summing all N crack subunits q, it represents the cumulative stress contribution of each subunit at p;

[0090] When solving formula (2) in the iterative process, firstly calculate Δu q Make an initial guess and calculate the value based on the initial guess And the corresponding energy norm E(u), the energy norm can be defined as:

[0091]

[0092] in, is the strain tensor, which is expressed by the displacement Δu q Substitute into the geometric equation to calculate:

[0093]

[0094] : represents the double inner product operation between tensors, Ω is the volume domain of shale, ρ is the density of shale body force (here only gravity is included), which affects the overall force balance and energy distribution, and x and y are the x-direction and y-direction in the coordinate system, respectively.

[0095] According to the principle of least action, when there is a deviation between the estimated value of the displacement field and the actual value, the energy state of the system will deviate from the minimum value. Therefore, the energy norm is optimized as the objective function. By minimizing the energy norm, the displacement field can be adjusted to make it closer to the actual physical state. The optimization method using gradient descent is used to adjust Δu q To minimize the energy norm E(u), so as to gradually make Δu q and Converges to the true value.

[0096] Step 122: According to the obtained Δu q and Update the crack morphology. Based on fracture mechanics theory, this process involves calculating the stress intensity factors (SIFs) at the crack tip and using the energy release rate (ERR) or crack growth criterion to predict the crack behavior. The specific steps are as follows:

[0097] Step 1221: Calculate the stress intensity factor. The stress intensity factor is an important parameter that describes the tendency of crack propagation. The stress intensity factors of mode I (cracking) and mode II (slip) are expressed as:

[0098]

[0099]

[0100] Among them, σ and τ are normal stress and shear stress respectively, and normal stress and shear stress are The normal stress will cause the crack to crack, and the shear stress will cause the crack to slip. The crack length of the entire crack surface at the crack unit p, Y, Y' is a geometric factor related to the crack shape and loading conditions. The geometric factor depends on the normal stress and shear stress on the crack surface, the geometric characteristics of the loaded object, and the position of the load. When specifically calculated, it is generally calculated through finite element simulation;

[0101] Step 1222: Calculate the energy release rate G, which is the energy released per unit crack area growth. For plane stress and plane strain conditions, G can be calculated using the following formula:

[0102]

[0103] Where E' is the equivalent elastic modulus. For plane stress problems, E' = E, and for plane strain problems, we have

[0104] ν is Poisson's ratio;

[0105] Step 1223: According to the crack extension criterion, whether the crack extends depends on whether G reaches or exceeds the fracture toughness G of the material. c , if G ≥ G c , the crack will expand; the direction of crack expansion is usually determined by the Maximum Energy Release Rate Criterion.

[0106] Step 1224: When the crack extension (G≥G c ), update the crack morphology, otherwise this step is not performed; updating the crack morphology includes the following steps:

[0107] 1) Add a certain amount of stress to get the new stress field σ i ' j (x), and let G = Gc;

[0108] 2) The new stress field σ i ' j (x) and the new G are substituted into formulas (5), (6), and (7) to obtain the new crack length a new and the crack inclination angle θ,

[0109]

[0110] 3) According to the obtained a new The fracture morphology is updated by using θ and θ, wherein the fracture morphology includes the fracture length and the fracture inclination, and the fracture morphology obtained each time is saved.

[0111] Step 123: Advance the time step and return to step 121 to perform calculations on the next subunit until the end time of the predetermined time series is reached. The multiple groups of simulated fracture morphologies are saved as data files, i.e., random values of fracture combination parameters. The maximum time series corresponding to each fracture surface is N. In practice, multiple fracture surfaces can be simulated.

[0112] Step 13: Combine the generated original geological parameter random values and the fracture combination parameter random values to form a complete fracture network expansion simulation data sample set. Each data sample in the data sample set includes a set of original geological parameters (such as porosity and permeability) and corresponding fracture combination parameters (fracture length and fracture inclination). Through this combination, a multidimensional feature vector is created, which covers both geological characteristics and fracture characteristics.

[0113] Step 14: Use principal component analysis (PCA) to map high-dimensional data samples to a low-dimensional feature space, retaining the main information and reducing feature redundancy:

[0114] Step 141: Standardize the data sample so that its mean is 0 and its variance is 1;

[0115] Step 142: Calculate the covariance matrix C:

[0116]

[0117] Find the eigenvectors and eigenvalues of the covariance matrix C, that is, solve the equation Cv = λv, where λ is the eigenvalue vector, v is the corresponding eigenvector, n represents the number of data samples, x i is the feature vector of the i-th data sample, is the average vector of all data samples in the data sample set;

[0118] Step 143: Select the largest m eigenvalues in λ: λ1≥λ2≥…≥λ m The corresponding eigenvectors constitute the transformation matrix W=[v1,v2,…v m ],

[0119] Let Z = XW,

[0120] X is the data matrix of the data sample, each row of X is a data sample, each column is a feature, Z is the principal component matrix, each row of Z is a sample, and each column is a principal component.

[0121] Step 15: Use statistics to construct more effective high-level features and scale multiple samples of the same feature, that is, each column of the principal component matrix Z. The following two methods can be used:

[0122] 1) Standardization: Scale the features to the same scale, where z is one of the columns of the principal component matrix Z, z' is the standardized result, μ is the mean of the column, and σ is the standard deviation of the column.

[0123]

[0124] 2) Normalization: Scale the features to the interval [0,1], where z is one column of the principal component matrix Z and z' is the normalized result.

[0125]

[0126] After this step, the feature dataset is obtained:

[0127]

[0128] in It is a set of feature data. is the number of eigenvectors contained in the set, F n is the nth eigenvector Its dimensional components Represents the high-level features after principal component decomposition and standardization / normalization, corresponding to the i-th standardized / normalized principal component z' n,i ,i∈(1,m).

[0129] Step 2: Construct a feature dataset based on fractal structures. Fractal structures are a set of complex geometric forms inherent to fractures, which help describe and simulate the formation and propagation trends of faults and natural fractures. Because rock fracture networks are self-similar and self-affine, fractal structure parameters more comprehensively describe the complexity and multi-scale characteristics of fracture networks. This allows for a more accurate description of complex three-dimensional fracture network morphology, thereby improving the ability to predict the mechanical behavior and permeability characteristics of fracture networks.

[0130] Considering the scale invariance of fractal structures, i.e., the statistical characteristics of parameters remain unchanged at different scales, it is possible to calculate them using statistical methods: the fractal dimension is calculated using the box counting method, and the fractal length distribution index is extracted by statistically analyzing the scale invariance of the crack length distribution. The box counting method is a commonly used method for calculating fractal dimension, particularly suitable for crack structures in two-dimensional and three-dimensional space. The steps are as follows:

[0131] Here are the steps:

[0132] Step 21. Cover the crack structure with grids of different sizes and count the number of boxes N(∈) covering the crack for each grid size, where ∈ is the side length of the grid. Plot the double logarithm logN(ε) and log(1 / ε). The slope of the graph is the fractal dimension D, which can be calculated using the formula:

[0133]

[0134] Step 21: Calculate the fractal length distribution. First, count the lengths L of the cracks and plot the corresponding probability distribution graph and the graph of log P(L) versus log Ld. P(L) is the ratio of the number of cracks within a specified length to the total number of cracks. The slope of the calculated length distribution graph is the length distribution index α.

[0135] Step 22: adding the obtained fractal dimension D and the length distribution index α as additional features to the feature data set obtained in step 1 to obtain a new feature data set;

[0136] Step 23: Obtain the crack network extension morphology distribution data at multiple prediction moments as corresponding prediction labels, add the prediction labels to the feature data set obtained in step 22 to obtain a new feature data set, and save the newly obtained feature data set in a data frame format (DataFrame) and export it as a CSV file;

[0137] Step 24: The feature dataset obtained in step 23 is randomly divided into three non-overlapping datasets: a training set, a validation set, and a test set in a ratio of 6:2:2. See Table 1. The training set is used for model training to help the model learn data features; the validation set is used for verification during model training, to adjust model parameters and evaluate model performance; and the test set is used to ultimately evaluate the generalization ability of the model to ensure that the model performs well on unseen data.

[0138] Table 1 Division of feature datasets

[0139]

[0140] Step 3: Establish a prediction model based on neural network:

[0141] Step 31: Set the first layer as a convolutional layer, which is used to capture the relationship between geological parameters or the pattern and distribution characteristics of fractures in space;

[0142] Step 32: The second layer is provided with a flattening layer, which is used to flatten the multi-dimensional tensor data output by the convolutional layer into a one-dimensional vector;

[0143] Step 33: The third layer is provided with a multi-layer perceptron (MLP): used to extract global features from the one-dimensional vector output by the flattening layer and capture the relationship between geological parameters, fracture combination parameters and analytical structural data. Two fully connected layers are provided, and the number of neurons in each layer can be adjusted according to data complexity and model requirements. The multi-layer perceptron composed of two fully connected layers is used to extract its global features and capture the relationship between geological parameter features, fracture combination parameters and analytical structural data.

[0144] Step 34. Select the activation function ReLU.

[0145] Step 4: Input the training set obtained in step 2 into the prediction model, and use the validation set for verification and the test machine to test the prediction results.

[0146] The above shows and describes the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions merely illustrate the principles of the present invention. Various changes and modifications may be made to the present invention without departing from the spirit and scope of the present invention. Such changes and modifications are intended to fall within the scope of the present invention. The scope of protection claimed by the present invention is defined by the appended claims and their equivalents.

Claims

1. A modeling method for a reservoir rock fracture network expansion prediction model based on a neural network, characterized in that: The specific steps include: Step 1: Obtain original geological parameters and fracture combination parameters of each reservoir rock; generate multiple random values of original geological parameters based on the original geological parameters, and generate multiple random values of fracture combination parameters based on the original geological parameters and fracture combination parameters, and combine the random values of original geological parameters and fracture combination parameters to form a fracture network expansion simulation data set; Step 2: Obtain fractal parameters of the fractal structure, combine the fracture network extension simulation data set and the fractal parameters, construct a feature data set based on the fractal structure, and divide the feature data set based on the fractal structure into a training set, a validation set, and a test set; Step 3: Establish a prediction model based on neural network; Step 4: Input the training set into the prediction model for training. After the training is completed, input the validation set into the prediction model for verification and input the test set into the prediction model for testing. Generating multiple random values of fracture combination parameters according to original geological parameters and fracture combination parameters includes the following steps: Step 11: Setting initial parameters and boundary conditions: The initial parameters include geological parameters, shale Young's modulus and Poisson's ratio, initial fracture length and initial fracture inclination; the geological parameters include porosity and permeability; the boundary conditions include stress field and displacement conditions, and the boundary conditions and initial parameters are input into the displacement discontinuity method model; Step 12: Simulate the propagation of the crack by solving the displacement on the crack surface using the displacement discontinuity method: Calculate the stress and displacement on the crack surface and update the crack morphology at each time step. Through iterative calculation, obtain the crack morphology at different moments in the predetermined time series. Combine the obtained multiple groups of crack morphologies to obtain random values of the crack combination parameters. The crack morphology includes crack length and crack inclination. Step 13: combining the generated original geological parameter random values and the fracture combination parameter random values; Step 12 specifically includes the following steps: Step 121: Based on the initial parameters, boundary conditions, and basic principles of elasticity, the stress boundary integral equation is used to solve the displacement discontinuity on the crack surface: in, is the initial ground stress field, σ ij (x) is the total stress field after considering the effect of cracks, Γ is the crack surface, ∫ Γ G ijkl (x, y) is the Green's function, i, j, k, and l represent the components of the tensor in the Green's function, i and j represent the displacement or stress component direction of the observation point, that is, the shear stress acting on the i-th surface along the j-th direction, k and l represent the direction of the force or stress applied by the source point, x represents the corresponding position of the calculated stress point, y represents the point on the integration path, Δu is the displacement on the crack surface, and the process of solving Δu is as follows: The crack surface is discretized into N sub-elements, and the displacement on the sub-element is Δu q , through the above discretization, for each crack element p, the stress boundary integral equation can be transformed into a linear algebraic equation: Among them, p and q are both discrete crack units, p represents the observation point or the specified discrete position for calculating stress; q represents the subunit index on the crack surface, which is used to traverse all subunits. is the total stress field at subunit p, Δu q is the displacement on the qth subunit, represents the Green's function at subunit p, which describes the influence of unit q on the stress at p. By summing all N crack subunits q, it represents the cumulative stress contribution of each subunit at p; When solving formula (2) in the iterative process, firstly calculate Δu q Make a guess as the initial value and calculate the And the corresponding energy norm E(u), the energy norm is defined as: Where ò(u) is the strain tensor, which is expressed by the displacement Δu q Substitute into the geometric equation to calculate: represents the double inner product operation between tensors, Ω is the volume domain of shale, ρ is the density of shale volume force, which affects the overall force balance and energy distribution, x and y are the x-direction and y-direction in the coordinate system respectively; By adjusting Δu q The value of is to minimize the energy norm E(u), so as to gradually make Δu q and Converge to the true value; Step 122: According to the obtained Δu q and Update the shape of the crack. The specific steps are as follows: Step 1221, calculate the stress intensity factor, which is divided into the cracking intensity factor K I and slip intensity factor K II , respectively expressed as: Among them, σ and τ are normal stress and shear stress respectively, and normal stress and shear stress are The weight, is the crack length, Y and Y are geometric factors related to the crack shape and loading conditions; Step 1222: Calculate the energy release rate G. G is calculated using the following formula: Where E′ is the equivalent elastic modulus. For plane stress, E′=E, and for plane strain, ν is Poisson's ratio; Step 1223: According to the crack extension criterion, whether the crack extends depends on whether G reaches or exceeds the fracture toughness G of the material. c , if G ≥ G c , the crack will expand; Step 1224: When the crack is determined to have expanded, the crack morphology is updated; otherwise, this step is not performed. Updating the crack morphology includes the following steps: 1) Add a certain amount of stress to get the new stress field σ′ ij (x), and let G = Gc; 2) The new stress field σ′ ij (x) and the new G are substituted into formulas (5), (6), and (7) to obtain the new crack length a new and the crack inclination angle θ, 3) According to the obtained a new and θ to update the fracture morphology, wherein the fracture morphology includes fracture length and fracture inclination, and the fracture morphology obtained each time is saved; Step 123, return to step 121, calculate the next sub-unit until all sub-units are calculated, and save the multiple groups of fracture morphologies obtained by simulation as data files, which are random values of fracture combination parameters.

2. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 1, characterized in that: Before executing step 2, the fracture network extension simulation dataset is processed as follows: Step 14: Map the combination obtained in step 13 from high dimension to low dimensional feature space through principal component analysis; Step 15: construct an advanced feature dataset based on the data obtained in step 14. The advanced feature dataset is the fracture network extension simulation dataset in step 2.

3. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 2, characterized in that: Step 14 specifically includes the following steps: Step 141: Standardize the data sample so that its mean is 0 and its variance is 1; Step 142: Calculate the covariance matrix C: Find the eigenvectors and eigenvalues of the covariance matrix C, that is, solve the equation Cv = λv, where λ is the eigenvalue vector, v is the corresponding eigenvector, n represents the number of data samples, x i is the feature vector of the i-th data sample, is the average vector of all data samples in the data sample set; Step 143: Select the largest m eigenvalues in λ: λ1≥λ2≥…≥λ m The corresponding eigenvectors constitute the transformation matrix W=[v1,v2,…v m ], Let Z = XW, X is the data matrix of the data sample, each row of X is a data sample, each column is a feature, Z is the principal component matrix, each row of Z is a sample, and each column is a principal component.

4. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 3, characterized in that: Step 15 uses one of the following two methods: 1) Standardization: Scale the features to the same scale, where z is one of the columns of the principal component matrix Z, z' is the standardized result, μ is the mean of the column, and σ is the standard deviation of the column. 2) Normalization: Scale the features to the interval [0,1], where z is one column of the principal component matrix Z and z' is the normalized result. After this step, the feature dataset is obtained: in It is a set of feature data. is the number of eigenvectors contained in the set, F n is the nth eigenvector z' n , its dimensional components Represents the high-level features after principal component decomposition and standardization / normalization, corresponding to the i-th standardized / normalized principal component z′ n,i ,i∈(1,m).

5. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 4, characterized in that: Step 2: Construct a feature dataset based on fractal structure and calculate the fractal dimension using the box counting method. The steps are as follows: Step 21. Cover the crack structure with grids of different sizes and count the number of boxes N(∈) covering the crack for each grid size, where ∈ is the side length of the grid. Plot the double logarithm logN(ε) and log(1 / ε). The slope of the graph is the fractal dimension D, which can be calculated using the formula: Step 22: Calculate the fractal length distribution. First, collect the crack length L data and draw the corresponding probability distribution graph and the graph of log P(L) versus log L. P(L) is the ratio of the number of cracks of a specified length to the total number of cracks. The slope of the length distribution graph is the length distribution index α. Step 23: adding the obtained fractal dimension D and the length distribution index α as additional features to the fracture network extension simulation data set obtained in step 1 to obtain a feature data set based on the fractal structure; Step 24: Obtain the crack network extension morphology distribution data at multiple prediction moments as corresponding prediction labels, add the prediction labels to the feature data set based on the fractal structure obtained in step 22 to obtain a final feature data set, save the final feature data set in a data frame format, and export it as a CSV file; Step 25: The final feature dataset obtained in step 23 is randomly divided into a training set, a validation set, and a test set in a ratio of 6:2:

2. The training set is used for model training, the validation set is used for verification during model training, and the test set is used to finally evaluate the generalization ability of the model.

6. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 5, characterized in that: The step three specifically includes the following steps: Step 31: Set the first layer as a convolutional layer, which is used to capture the relationship between geological parameters or the pattern and distribution characteristics of fractures in space; Step 32: The second layer is provided with a flattening layer, which is used to flatten the multi-dimensional tensor data output by the convolutional layer into a one-dimensional vector; Step 33: The third layer is provided with a multilayer perceptron, which is used to extract global features from the one-dimensional vector output by the flattening layer and capture the relationship between geological parameters, fracture combination parameters and analytical structural data. Two fully connected layers are provided, and the number of neurons in each layer is adjusted according to the data complexity and model requirements. The multilayer perceptron composed of two fully connected layers is used to extract its global features. Capture the relationship between geological parameter characteristics, fracture combination parameters and analytical structural data; Step 34. Select the activation function ReLU.

7. A method for predicting reservoir rock fracture network expansion based on neural network, characterized in that: The specific steps include: Step 1: Obtain the original geological parameters and natural fracture distribution parameters of the reservoir rock to be predicted; Step 2: Input the original geological parameters and natural fracture distribution parameters into the model established and trained using the modeling method described in any one of claims 1 to 6 to obtain the fracture morphology.

Images