Neural network-based reservoir rock fracture network extension prediction model modeling method and prediction method

Abstract

The invention provides a modeling method of a reservoir rock fracture network extension prediction model based on a neural network. The modeling method specifically comprises the following steps that 1, original geological parameters, fracturing design parameters and natural fracture distribution parameters of each type of oil reservoir are obtained; generating a plurality of original geological parameter random values according to the original geological parameters, generating a plurality of crack combination parameter random values by using the crack combination parameters, and combining the original geological parameter random values and the crack combination parameter random values to form a crack network expansion simulation data set; step 2, acquiring fractal parameters of the fractal structure, constructing a feature data set based on the fractal structure in combination with the fracture network expansion simulation data set and the fractal parameters, and dividing the feature data set into a training set, a verification set and a test set; step 3, establishing a prediction model based on a neural network; and step 4, inputting the training set into the prediction model for training, inputting the verification set into the prediction model for verification after training is completed, and inputting the 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 sedimentary bedding and natural fractures are highly developed in shale reservoirs. The complex fracture network formed by volume fracturing in shale reservoirs provides an effective channel for oil and gas flow. Rapidly and accurately predicting the expansion morphology of the complex fracture network formed by shale reservoir fracturing is of great significance to the efficient development of shale oil and gas reservoirs.

[0003] In the current research on the expansion of shale reservoir fracture networks, there are problems such as the following: 1) Most existing studies are based on finite element analysis. This method relies on accurate physical models and solving partial differential equations, and requires meshing and iterative solutions. Therefore, the solution process is computationally complex and time-consuming. In addition, FEM needs to adjust the model and boundary conditions for different shale samples to re-model and solve, and the transferability is relatively poor; 2) The prediction of shale reservoir fracture network expansion involves the morphology of the fracture network. These rock fracture networks have self-similar and self-affine fractal characteristics. The fractal structure contains the spatial distribution and multi-scale characteristics of the fracture, and affects its connectivity and permeability, which has an important impact on the fracture expansion path and pattern. However, traditional prediction methods often ignore the fractal structure when characterizing the fracture network morphology, resulting in inaccurate prediction results, affecting the evaluation of fracture expansion behavior, permeability and hydraulic fracturing effect. 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-mentioned technical problems.

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

[0006] Step 1, obtaining original geological parameters and fracture combination parameters of each reservoir rock; generating multiple random values ​​of original geological parameters according to the original geological parameters, and generating multiple random values ​​of fracture combination parameters according to the original geological parameters and the fracture combination parameters, and combining the random values ​​of the original geological parameters and the random values ​​of the fracture combination parameters to form a fracture network expansion simulation data set;

[0007] Step 2: obtaining fractal parameters of the fractal structure, combining the fracture network extension simulation data set and the fractal parameters, constructing a feature data set based on the fractal structure, and dividing 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, and 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 model of the displacement discontinuity method;

[0012] Step 12, solving the displacement on the crack surface by the displacement discontinuity method to simulate the expansion of the crack: in each time step, the stress and displacement on the crack surface are calculated and the crack morphology is updated. Through iterative calculation, the crack morphology at different times in the predetermined time series is obtained, and the obtained multiple groups of crack morphologies are combined to obtain the random value of the crack combination parameter. The crack morphology includes the crack length and the 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 data set is processed as follows:

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

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

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

[0018] Step 121, based on the initial parameters and boundary conditions and the basic principles of elasticity, the displacement discontinuity on the crack surface is solved using the stress boundary integral equation:

[0019]

[0020] in, is the initial geostress 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, 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 integral 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-units, 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:

[0022]

[0023] Among them, p and q are both discrete crack units, p represents the specified discrete position of the observation point or calculated 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, represents the Green's function at 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, Δ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, expressed by the displacement Δu q Substitute the geometric equation into the calculation:

[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 the slip intensity factor K II , respectively expressed as:

[0032]

[0033]

[0034] Among them, σ and τ are normal stress and shear stress, respectively. 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, where G 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 it is determined that the crack has 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 a new stress field σ i ' j (x), and let G = Gc;

[0041] 2) The new stress field σ i ' j Substitute (x) and the new G 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 eigenvector and eigenvalue 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 ith 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 form the transformation matrix W = [v1, v2, …v m ],

[0051] Let Z = XW,

[0052] X is the data matrix of data samples, 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 of the columns of the principal component matrix Z, and z' is the normalized result.

[0057]

[0058] After this step, the feature data set 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, step 2, constructing a feature data set based on a fractal structure, and calculating the fractal dimension by a box counting method, 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 cracks under each grid size, where ∈ is the side length of the grid; draw a double logarithm logN(ε) and log(1 / ε) graph, and the slope in the graph is the fractal dimension D, that is, the calculation formula of the fractal dimension is:

[0063]

[0064] Step 21, calculating the fractal length distribution, firstly, counting the length L data of the cracks, and drawing the corresponding probability distribution diagram and log P(L) versus log L, where P(L) is the ratio of the number of cracks of a certain length to the total number of cracks, and calculating the slope of the length distribution diagram and the width distribution diagram is 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, obtaining the crack network extension morphology distribution data at multiple prediction moments as corresponding prediction labels, adding the prediction labels to the feature data set based on the fractal structure obtained in step 22, obtaining the final feature data set, and saving the final feature data set in a data frame format (DataFrame), and exporting it as a CSV file;

[0067] Step 24: The final feature data set 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, setting the first layer as a convolutional layer, the convolutional layer is used to capture the relationship between geological parameters or the pattern and distribution characteristics of cracks in space;

[0070] Step 32: The second layer is provided with a flattening layer, where the flattening layer 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, and two fully connected layers are set. 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 prediction efficiency of the shale reservoir fracturing fracture network expansion. 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 used to disclose the present invention so that those skilled in the art can implement the present invention. The preferred embodiments described below are only examples, and those skilled in the art may think of other obvious variations.

[0077] like Figure 1 As shown, a modeling method of 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: sort out and analyze the data on site to obtain the original geological parameters and fracture combination parameters of each type of oil reservoir; generate multiple random values ​​of original geological parameters based on the original geological parameters and use geostatistical methods, and use the displacement discontinuity method model to simulate the fracture network expansion and generate multiple random values ​​of fracture combination parameters. Combine the original geological parameter random values ​​and the fracture combination parameter random values ​​to form a fracture network expansion simulation data set. The original geological parameters include porosity distribution, permeability distribution, and original geostress field, and the fracture combination parameters include fracture length and fracture 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 original geological parameters, assuming that the original geological parameters conform to normal distribution, and setting the 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. Random values ​​of these geological parameters can be generated by Monte Carlo simulation.

[0080] According to the original geological parameters and fracture combination parameters, the initial parameters and boundary conditions of the displacement discontinuity method model are initialized. In each simulation, the fracture combination parameters are used as output. 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 extension morphology of fracture networks.

[0081] Specifically, the displacement discontinuity method is used to simulate the expansion 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, and the geological parameters include porosity, permeability, etc.; the boundary conditions include stress field and displacement conditions, and 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 in 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. Specifically, the following steps are included:

[0084] Step 121, according to the initial parameters and boundary conditions, the discontinuous displacement and total stress field on the crack surface are calculated. According to the basic principles of elastic mechanics, the displacement discontinuity on the crack surface can be solved by using the stress boundary integral equation:

[0085]

[0086] in, is the initial geostress 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 a unit force at 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 a 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 fracture 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 discrete point p, describing 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, Δu q Make an initial guess and calculate the And the corresponding energy norm E(u), the energy norm can be defined as:

[0091]

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

[0093]

[0094] : represents the double inner product operation between tensors, Ω is the volume domain of shale, ρ is the density of shale volume force (here only gravity is included), which affects the overall force balance and energy distribution, and x and y are the x and y directions 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 real 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. According to 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 extension criteria to predict the behavior of the crack. The specific steps are as follows:

[0097] Step 1221, calculate the stress intensity factor. The stress intensity factor is an important parameter to describe the crack extension tendency. The stress intensity factors of mode I (cracking) and mode II (slipping) are respectively expressed as:

[0098]

[0099]

[0100] Among them, σ and τ are normal stress and shear stress, respectively. 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, the position of the load, etc. When calculating specifically, it is generally calculated by finite element simulation;

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

[0102]

[0103] Where E' is the equivalent elastic modulus. For plane stress problems, E' = E. 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 executed; updating the crack morphology includes the following steps:

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

[0108] 2) The new stress field σ i ' j Substitute (x) and the new G 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 fracture length and fracture inclination, and the fracture morphology obtained each time is saved.

[0111] Step 123, advance the time step, return to step 121, calculate the next subunit until the end time of the predetermined time series is reached, save the multiple groups of fracture morphologies obtained by simulation as data files, that is, the random values ​​of the fracture combination parameters, the maximum time series corresponding to each fracture surface is N, and multiple fracture surfaces can be simulated in practice;

[0112] Step 13, combining 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, permeability) and corresponding fracture combination parameters (fracture length, fracture inclination), and through this combination, a multidimensional feature vector is created, which covers both geological characteristics and fracture characteristics;

[0113] Step 14: Map the high-dimensional data samples to the low-dimensional feature space through principal component analysis (PCA), retain the main information and reduce 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 eigenvector and eigenvalue 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 ith 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 form the transformation matrix W = [v1, v2, …v m ],

[0119] Let Z = XW,

[0120] X is the data matrix of data samples, 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: With the help of statistics, more effective high-level features are constructed. The scale of multiple samples of the same feature, that is, each column of the principal component matrix Z, can be scaled by the following two methods:

[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 of the columns of the principal component matrix Z, and z' is the normalized result.

[0125]

[0126] After this step, the feature data set 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 data set based on fractal structure. Fractal structure is a set of complex geometric forms inherent in fractures, which helps to describe and simulate the formation and expansion trends of faults and natural fractures. Based on the self-similar and self-affine properties of rock fracture networks, fractal structure parameters more comprehensively describe the complexity and multi-scale characteristics of fracture networks, and can more accurately describe the morphology of complex three-dimensional fracture networks, thereby improving the prediction ability of fracture network mechanical behavior and permeability characteristics.

[0130] Considering the scale invariance of fractal structures, that is, the statistical characteristics of parameters remain unchanged at different scales, it can be calculated by statistical methods: the fractal dimension is calculated by 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 dimensions, especially 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 cracks under each grid size, where ∈ is the side length of the grid; draw a double logarithm logN(ε) and log(1 / ε) graph, and the slope in the graph is the fractal dimension D, that is, the calculation formula of the fractal dimension is:

[0133]

[0134] Step 21, calculate the fractal length distribution. First, count the length L of the cracks, and draw the corresponding probability distribution diagram and the image of log P(L) versus log Ld. P(L) is the ratio of the number of cracks within the specified length to the total number of cracks. The slope of the calculated length distribution diagram 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, obtaining the crack network extension morphology distribution data at multiple prediction moments as corresponding prediction labels, adding the prediction labels to the feature data set obtained in step 22 to obtain a new feature data set, and saving the newly obtained feature data set in a data frame format (DataFrame), and exporting it as a CSV file;

[0137] Step 24: The feature data set obtained in step 23 is randomly divided into three non-overlapping data sets: training set, validation set, and test set according to 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; the test set is used to finally evaluate the generalization ability of the model to ensure that the model performs well on unseen data;

[0138] Table 1 Division of feature data sets

[0139]

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

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

[0142] Step 32: The second layer is provided with a flattening layer, where the flattening layer 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 flattened layer and capture the relationship between geological parameters, fracture combination parameters and analytical structural data, set two fully connected layers, the number of neurons in each layer can be adjusted according to data complexity and model requirements, and use a multi-layer perceptron composed of two fully connected layers 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 check 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 only describe the principles of the present invention. The present invention may be subject to various changes and improvements without departing from the spirit and scope of the present invention. These changes and improvements fall within the scope of the present invention. The scope of protection claimed by the present invention is defined by the attached 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, obtaining original geological parameters and fracture combination parameters of each reservoir rock; generating multiple random values ​​of original geological parameters according to the original geological parameters, and generating multiple random values ​​of fracture combination parameters according to the original geological parameters and the fracture combination parameters, and combining the random values ​​of the original geological parameters and the random values ​​of the fracture combination parameters to form a fracture network expansion simulation data set; Step 2: obtaining fractal parameters of the fractal structure, combining the fracture network extension simulation data set and the fractal parameters, constructing a feature data set based on the fractal structure, and dividing 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.

2. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 1 is characterized in that: 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, and 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 model of the displacement discontinuity method; Step 12, solving the displacement on the crack surface by the displacement discontinuity method to simulate the expansion of the crack: in each time step, the stress and displacement on the crack surface are calculated and the crack morphology is updated. Through iterative calculation, the crack morphology at different times in the predetermined time series is obtained, and the obtained multiple groups of crack morphologies are combined to obtain the random value of the crack combination parameter. The crack morphology includes the crack length and the crack inclination. Step 13: Combine the generated original geological parameter random values ​​and the fracture combination parameter random values.

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

4. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 2 is characterized in that: Step 12 specifically includes the following steps: Step 121, based on the initial parameters and boundary conditions and the basic principles of elasticity, the displacement discontinuity on the crack surface is solved using the stress boundary integral equation: in, is the initial geostress 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, 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 integral 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-units, 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: 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, 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; When solving formula (2) in the iterative process, firstly, Δ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: Where ò(u) is the strain tensor, which is expressed by the displacement Δu q Substitute the geometric equation into the calculation: 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, and x and y are the x and y directions 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 the slip intensity factor K II , respectively expressed as: Among them, σ and τ are normal stress and shear stress, respectively. The weight, is the crack length, Y, Y' is a geometric factor related to the crack shape and loading conditions; Step 1222: Calculate the capacity release rate G, where 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 it is determined that the crack has 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 Substitute (x) and the new G 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.

5. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 3 is 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 eigenvector and eigenvalue 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 ith 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 form the transformation matrix W = [v1, v2, …v m ], Let Z = XW, X is the data matrix of data samples, 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.

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: 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 of the columns of the principal component matrix Z, and z' is the normalized result. After this step, the feature data set 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).

7. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 5, characterized in that: Step 2: Construct a feature data set based on fractal structure and calculate the fractal dimension by 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 cracks under each grid size, where ∈ is the side length of the grid; draw a double logarithm logN(ε) and log(1 / ε) graph, and the slope in the graph is the fractal dimension D, that is, the calculation formula of the fractal dimension is: Step 22, calculating the fractal length distribution, firstly, counting the length L data of the cracks, and drawing the corresponding probability distribution graph and the image of log P(L) versus log L, where P(L) is the ratio of the number of cracks of a specified length to the total number of cracks, and the slope of the calculated 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, obtaining the crack network extension morphology distribution data at multiple prediction moments as corresponding prediction labels, adding the prediction labels to the feature data set based on the fractal structure obtained in step 22, obtaining the final feature data set, and saving the final feature data set in a data frame format (DataFrame), and exporting it as a CSV file; Step 25: The final feature data set 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.

8. The modeling method of a reservoir rock fracture network expansion prediction model based on a neural network according to claim 7, characterized in that: The step three specifically includes the following steps: Step 31, setting the first layer as a convolutional layer, the convolutional layer is used to capture the relationship between geological parameters or the pattern and distribution characteristics of cracks in space; Step 32: The second layer is provided with a flattening layer, where the flattening layer 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 analysis structure data. Two fully connected layers are set, 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. Capture the relationship between geological parameter characteristics, fracture combination parameters and analytical structural data; Step 34. Select the activation function ReLU.

9. A method for predicting reservoir rock fracture network expansion based on neural network, characterized in that: The specific steps include: Step 1: Obtaining 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 by the modeling method described in any one of claims 1 to 8 to obtain the fracture morphology.

Images