Method and device for constructing three-dimensional random geological model based on machine learning

Abstract

The present invention provides a method and device for constructing a three-dimensional random geological model based on machine learning. The construction method includes: constructing a geological model region and discretizing it into n voxels; constructing an initial geological model based on the label values ​​of known voxels determined by drilling data and the initial label values ​​of unknown voxels determined by a cubic learning domain; constructing a neighborhood system for each voxel and setting spatial constraints in 13 directions; obtaining a prediction result based on the label value of each voxel in the initial geological model according to a label probability prediction formula and an iterative update formula for the grain size coefficient; obtaining the edge probability of each unknown voxel and the predicted label value of each unknown voxel based on the prediction result; and updating the initial geological model based on the predicted label value of each unknown voxel to construct a three-dimensional random geological model. The three-dimensional random geological model constructed using the construction method of the present invention has accurate simulation results and can simultaneously characterize the uncertainty of the strata.

Description

Technical Field

[0001] The present invention belongs to the field of geotechnical engineering technology, and specifically relates to a method and device for constructing a three-dimensional random geological model based on machine learning. Background Art

[0002] A prerequisite for geotechnical reliability analysis is to fully quantify the sources of uncertainty in site conditions. Obtaining an accurate interpretation of the subsurface with quantified uncertainty is a critical and essential component of geotechnical system planning and design. However, soil, as a naturally occurring material, is inherently heterogeneous and stochastic, making accurate interpretation of site conditions particularly complex. Limited by the availability of field data, inferring the subsurface structure of a specific project site is already challenging, let alone quantifying its associated uncertainties.

[0003] Due to budget constraints and tight project timelines, typical geotechnical engineering projects typically only have a limited number of boreholes available. Consequently, subsurface structural analysis often relies on sparsely distributed borehole sampling data. Stratigraphic information at other locations often requires inference or simulation based on historical drilling data or relevant field observations. Due to the complex and limited understanding of geological formation mechanisms, coupled with a limited number of boreholes and insufficient test data, inferring subsurface structure faces significant uncertainty. This uncertainty not only impacts the accuracy of geological modeling but also poses potential risks to subsequent engineering design.

[0004] In recent years, machine learning-based models, such as Markov random fields (MRFs), coupled Markov chains (CMCs), and iterative convolutional XGBoost (IC-XGBoost), have been widely used in geotechnical engineering to more efficiently assess subsurface uncertainty. These models offer innovative technical solutions for addressing subsurface uncertainty. However, because these models primarily record and present geological data in two dimensions, a significant amount of spatial information is often overlooked during the modeling process. Summary of the Invention

[0005] The purpose of the present invention is to provide an efficient three-dimensional random geological model construction method based on Markov random field theory and Bayesian machine learning. The construction method uses 13 granularity coefficients to describe spatial correlation. The initial label values ​​of unknown voxels in the initial geological model that do not correspond to the drill hole position are determined according to the values ​​of known voxels in a cubic learning domain established with the unknown voxel as the center. Subsequently, the method is further refined through Bayesian machine learning, with fast iterative convergence. The constructed three-dimensional random geological model has accurate simulation results and can characterize the uncertainty of the stratum.

[0006] In a first aspect, the present invention provides a method for constructing a three-dimensional random geological model based on machine learning, the construction method comprising the following steps: step S10, constructing a three-dimensional geological model region based on the geological body of the modeling area, and discretizing the geological model region into n voxels, and numbering each voxel, wherein the modeling area includes multiple boreholes opened in the geological body; step S20, mapping the drilling position and drilling information of each borehole in the modeling area to the n voxels, completing label assignment for multiple known voxels corresponding to the drilling position, and determining the initial label value of the unknown voxel that does not correspond to the drilling position based on the label value of the known voxel in the cubic learning domain established with the unknown voxel as the center, and constructing an initial geological model based on the label value of the known voxel and the initial label value of the unknown voxel; step S30, constructing a neighborhood system for each unknown voxel based on the Markov random field, and setting spatial constraints in 13 directions for each unknown voxel, and at the same time using granularity coefficients β1~β 13 Describe the spatial constraint strength in each different direction, wherein the neighborhood system of each unknown voxel consists of other voxels that are adjacent to the unknown voxel and share nodes; step S40, based on the initial label value of each unknown voxel, the label value of each known voxel and the neighborhood system of each unknown voxel, obtain a prediction result according to the label probability prediction formula and the update iterative formula of the granularity coefficient, the prediction result includes M sampled label values ​​of each unknown voxel, M is a positive integer; step S50, based on the M sampled label values ​​of each unknown voxel, obtain the edge probability of each unknown voxel, and based on the edge probability of each unknown voxel, obtain the predicted label value of each unknown voxel; step S60, update the initial geological model based on the predicted label value of each unknown voxel to construct a three-dimensional random geological model.

[0007] In a specific embodiment, the construction method also includes step S70, obtaining the information entropy of each unknown voxel based on the calculation formula of the edge probability and information entropy of each unknown voxel, and obtaining the quantification result of the stratigraphic uncertainty corresponding to the three-dimensional random geological model based on the information entropy of each unknown voxel.

[0008] In a specific embodiment, in step S20, based on the label value of each known voxel, the initial label value of each unknown voxel is determined in sequence according to a preset order, wherein the step of determining the initial label value of any unknown voxel in the initial geological model includes: step (1), taking the unknown voxel i as the center voxel, constructing a cubic learning domain according to a preset size, and obtaining coordinate information of all voxels in the cubic learning domain, wherein the voxels in the cubic learning domain are divided into first voxels with known label values ​​and second voxels with unknown label values, wherein the coordinate information is represented by a coordinate index; step (2), obtaining the label value of each first voxel, based on each second voxel, The coordinate information and label value of a voxel are used to obtain a first covariance matrix between multiple first voxels with the same label value and a second covariance matrix between multiple first voxels with different label values; step (3), based on the first covariance matrix and the second covariance matrix, obtain the local measurement of the central voxel; step (4), based on the coordinate information and label value of each first voxel, the coordinate information of the central voxel and the local measurement of the central voxel, obtain the DANN distance between each first voxel and the central voxel; step (5), based on the DANN distance between each first voxel and the central voxel and the local probability calculation formula, obtain the initial label value of the central voxel.

[0009] In a specific embodiment, the step (5) includes: calculating the DANN distance harmonic mean for evaluating the possibility of selecting each label based on the DANN distance between the first voxel and the central voxel with the same label value; obtaining the label value probability distribution result of the central voxel based on the DANN distance harmonic mean for evaluating the possibility of selecting each label and the local probability calculation formula; and obtaining the initial label value of the central voxel by random sampling based on the label value probability distribution result of the central voxel.

[0010] In a specific embodiment, the formula for calculating the harmonic mean of the DANN distance used to evaluate the likelihood of selecting the label value m is:

[0011] ;

[0012] Where, is the harmonic mean distance from unknown voxel i to the first voxel g with the k nearest label value m. The DANN distance used to calculate the harmonic mean distance is based on Select, The matrix representing the DANN distance between each first voxel and the central voxel, x g is the label value of the first voxel g;

[0013] The local probability calculation formula is:

[0014] ;

[0015] Where, represents the local probability of selecting label m for unknown voxel i, It means that all possible labels in the first voxel are looped over L, where L is the set of label values.

[0016] In a specific embodiment, the local metric of the central voxel The calculation formula is:

[0017] ;

[0018] Where W is the first covariance matrix, B is the second covariance matrix, Is a tuning parameter, I is a multidimensional identity matrix, and its dimensions are consistent with W and B;

[0019] A matrix containing the DANN distances between each first voxel and the central voxel The calculation formula is:

[0020] ;

[0021] in, is the coordinate information of the first voxel g, is the coordinate information of the unknown voxel i, c, r, d represent the index in the X, Y, and Z directions respectively, represents the local metric of the central voxel and T represents the matrix transpose.

[0022] In a specific embodiment, the label probability prediction formula includes a local conditional probability distribution calculation formula, a local energy calculation formula, and a potential function calculation formula, wherein:

[0023] The local conditional probability distribution calculation formula is:

[0024] ;

[0025] The local energy calculation formula is:

[0026] ;

[0027] The potential function calculation formula is:

[0028] ;

[0029] In the label probability prediction formula, is the local conditional probability distribution at the unknown voxel i; q is and neighborhood system energy The weight coefficient of It is a potential function defined solely on the unknown voxel i, which represents the preference for selecting different labels on the unknown voxel i; represents local energy; x i represents a specific label at the unknown voxel i, Represents any possible label at the unknown voxel i circulating on L, where L is the set of label values; is the potential function, represents the neighborhood system of unknown voxel i; β d ∈β={β1,β2,β3,β4,β5,β6,β7,β8,β9,β 10 , β 11 , β 12 , β 13}, x j Represents a specific label of other voxels j in the neighborhood system of unknown voxel i;

[0030] The update iterative formula of the granularity coefficient includes a likelihood function calculation formula and an objective function calculation formula, wherein:

[0031] The likelihood function calculation formula is:

[0032] ;

[0033] The objective function calculation formula is:

[0034] ;

[0035] In the update iterative formula of the granularity coefficient, β represents the granularity coefficient, represents the prior distribution of β, x unknown Represents the sampled label value of the unknown voxel, x BH represents the label value of a known voxel, represents the neighborhood system of unknown voxel i, It represents the probability of unknown voxel i taking each label value based on the voxel label value and β value in the neighborhood system of unknown voxel i.

[0036] In a specific embodiment, the step S40 includes: step S401, based on the preset β1~β 13 Mean and standard deviation, establish the relationship between β1~β 13 Corresponding prior multivariate Gaussian distribution; Step S402, based on the label value of the known voxel, the sampled label value of the unknown voxel obtained by the last iteration update and the prior multivariate Gaussian distribution, according to the update iteration formula of the granularity coefficient, obtain β1~β 13 The sampling value of the current iteration is β1~β 13 The corresponding objective function value is greater than β1~β obtained in the last iteration13 Corresponding objective function value; Step S403, according to the β1~β 13 The sampling value and the label probability prediction formula are used to calculate the local conditional probability distribution of each unknown voxel obtained in the current iteration, and according to the local conditional probability distribution of each unknown voxel, the sampling label value of each unknown voxel obtained in the current iteration is obtained; step S404, repeating steps S402 and S403 multiple times until β1~β 13 Converges and in β1~β 13 After convergence, continue to repeat step S402 and step S403 for a preset number of M times, and set β1~β 13 After convergence, the M sampled label values ​​of each unknown voxel obtained by M iterations are used as the prediction results.

[0037] In a specific embodiment, the step S402 includes: step (a), using the Metropolis-Hastings method to sample from the prior multivariate Gaussian distribution to obtain β1~β 13 The current sampling value of step (b), based on the label value of the known voxel, the sampling label value of the unknown voxel obtained by the last iteration, the prior multivariate Gaussian distribution, and β1~β 13 The current sampling value of β1~β is obtained according to the likelihood function calculation formula and the objective function calculation formula. 13 The current objective function value corresponding to the current sampling value of the current objective function; Step (c), the current objective function value and the β1~β obtained in the last iteration 13 The corresponding objective function values ​​are compared. When the current objective function value is greater than β1~β obtained in the previous iteration, 13 When the corresponding objective function value is 13 The current sampling value is used as the β1~β obtained in the current iteration 13 The sampling value of the current objective function is less than or equal to the β1~β obtained in the previous iteration. 13 When the corresponding objective function value is reached, return to step (a) and resample.

[0038] In a second aspect, the present invention provides a device for constructing a three-dimensional random geological model based on machine learning, the device comprising: a first construction module for constructing a three-dimensional geological model region based on a geological body in a modeling area, and discretizing the geological model region into n voxels, and numbering each voxel, wherein the modeling area includes a plurality of boreholes opened in the geological body; a second construction module for mapping the drilling position and drilling information of each borehole in the modeling area to the n voxels, completing label assignment for a plurality of known voxels corresponding to the drilling position, and determining the initial label value of the unknown voxel that does not correspond to the drilling position based on the label value of the known voxel in a cubic learning domain established with the unknown voxel as the center, and constructing an initial geological model based on the label value of the known voxel and the initial label value of the unknown voxel; a third construction module for constructing a neighborhood system for each unknown voxel based on a Markov random field, and setting spatial constraints in 13 directions for each unknown voxel, and at the same time using granularity coefficients β1~β 13 Describe the spatial constraint strength in each different direction, wherein the neighborhood system of each unknown voxel is composed of other voxels adjacent to the unknown voxel and sharing nodes; a first acquisition module is used to obtain a prediction result based on the initial label value of each unknown voxel, the label value of each known voxel and the neighborhood system of each unknown voxel according to the label probability prediction formula and the update iterative formula of the granularity coefficient, wherein the prediction result includes M sampled label values ​​of each unknown voxel, where M is a positive integer; a second acquisition module is used to obtain the edge probability of each unknown voxel based on the M sampled label values ​​of each unknown voxel, and obtain the predicted label value of each unknown voxel based on the edge probability of each unknown voxel; an update module is used to update the initial geological model based on the predicted label value of each unknown voxel to construct a three-dimensional random geological model.

[0039] The beneficial effects of the present invention include at least:

[0040] 1. The present invention is based on Markov random field theory and Bayesian machine learning, and provides a method for constructing a three-dimensional random geological model based on machine learning. First, an initial geological model is constructed, wherein the label value of the voxels in the initial geological model is determined based on the drilling position and drilling information, and the initial label value of the unknown voxels in the initial geological model is determined according to the label value of the known voxels in the cubic learning domain established with the unknown voxel as the center; then a neighborhood system of each voxel is constructed, and 13-directional spatial constraints are set for each voxel, including all neighbors of the unknown voxel; based on the initial label value of each unknown voxel, the label value of each known voxel and the neighborhood system of each unknown voxel, a prediction result is obtained according to the label probability prediction formula and the update iterative formula of the granularity coefficient, and the prediction result includes M sampling label values ​​are set for each unknown voxel, where M is a positive integer; finally, the label value of the unknown voxel is determined according to the prediction result, and the initial geological model is updated based on the label value of the unknown voxel to obtain a three-dimensional random geological model; compared with the existing technology, on the one hand, spatial constraints in 13 directions are set, which include all neighbors of the unknown voxel, making the simulation results more accurate; on the other hand, the initial label value of the unknown voxel is determined based on the label value of the known voxel in the cubic learning domain, so that it is further accurate and the iterative convergence is fast in the subsequent Bayesian machine learning; at the same time, compared with the two-dimensional model, the three-dimensional random geological model constructed by the present invention is more intuitive, accurate and comprehensive in the expression of spatial geological bodies, and the geological information in the three-dimensional random model is richer, more real and reliable, which facilitates more in-depth geological analysis.

[0041] 2. The construction method provided by the present invention has a wide range of potential application value, especially in the field of geotechnical engineering, which is specifically reflected in the following aspects: (1) Improving the accuracy of stratum prediction: By combining Bayesian learning with the Markov random field model, the accuracy and reliability of stratum prediction can be significantly improved; (2) Supporting complex geotechnical engineering projects: This method is of great significance for the risk assessment and management of complex geotechnical engineering projects such as foundation treatment, slope reinforcement, and tunnel construction, and helps to reduce uncertainty and risk during the engineering design and implementation stages; (3) Enhancing risk control capabilities: By quantitatively analyzing stratum uncertainty, this method provides a more accurate tool for risk control of engineering projects, which can improve construction safety and sustainability.

[0042] In addition to the above-described objects, features and advantages, the present invention has other objects, features and advantages. The present invention will be further described in detail below with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 A schematic flow chart of the steps of a method for constructing a three-dimensional random geological model provided by one embodiment of the present invention;

[0044] Figure 2 for Figure 1 The sub-step flow chart of step S40 is shown;

[0045] Figure 3 A schematic diagram of an unknown voxel i and its neighborhood system provided by the present invention;

[0046] Figure 4 Schematic diagram of an unknown voxel i and its learning domain provided by the present invention;

[0047] Figure 5 This is a diagram showing the spatial constraints set in 13 directions for an unknown voxel i provided by the present invention, wherein: Figure 5 (a) Schematic diagram of the five intersecting planes corresponding to the spatial constraints; Figure 5 (b) A diagram showing the spatial constraints of the unknown voxel i in 13 directions based on five intersecting planes and coordinate axes.

[0048] Figure 6 A module diagram of a device for constructing a three-dimensional random geological model provided by another embodiment of the present invention;

[0049] Figure 7 The three-dimensional geological model and drilling schematic diagram constructed in Example 1 of the present invention, wherein: Figure 7 (a) is a three-dimensional geological body model diagram generated in Example 1; Figure 7 (b) in the figure shows the drilling diagram where 5% of the columns are randomly selected as known drilling holes;

[0050] Figure 8 The initial geological model map constructed in Example 1 of the present invention;

[0051] Figure 9 The three-dimensional random geological model constructed by the construction method of the present invention in Example 1 and the stratum uncertainty quantification result diagram corresponding to the three-dimensional random geological model are shown, wherein: Figure 9 (a) is the three-dimensional random geological model diagram constructed in Example 1. Figure 9 (b) is the quantification result diagram of stratigraphic uncertainty corresponding to the three-dimensional random geological model;

[0052] Figure 10 A three-dimensional schematic diagram of 33 drill holes provided in Example 2;

[0053] Figure 11 for Figure 10 A top view of the 33 drill holes shown;

[0054] Figure 12 This is a bar chart summarizing the accuracy of the cross-validation of 33 drill holes in Example 2. DETAILED DESCRIPTION

[0055] The embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0056] According to a first aspect of the present invention, based on Markov random field (MRF) theory and Bayesian machine learning, a method for constructing a three-dimensional random geological model using machine learning is provided. Compared with traditional two-dimensional models, the proposed three-dimensional random model is more intuitive, accurate, and comprehensive in representing spatial geological volumes. Furthermore, the geological information contained in the three-dimensional model is richer, more accurate, and more reliable, facilitating more in-depth geological analysis.

[0057] See also Figures 1 to 5 The method for constructing a three-dimensional random geological model based on machine learning provided by the present invention comprises the following steps:

[0058] Step S10: construct a three-dimensional geological model region based on the geological body of the modeling area, discretize the geological model region into n voxels, and number each voxel. The modeling region includes a plurality of boreholes opened in the geological body.

[0059] In the present invention, each voxel is in the shape of a cube.

[0060] In the present invention, the geological model region may be in the shape of a cuboid, or in the shape of an undulating hillside, etc., which is determined according to the actual shape of the geological body and is not limited by the present invention.

[0061] In the present invention, a three-dimensional geological model area is constructed based on a geological body of a modeling area. A plurality of sparse boreholes are opened on the geological body. The rock and soil type of the voxel corresponding to the borehole position can be determined through the borehole position and the borehole information.

[0062] In the present invention, each voxel is numbered to construct a set S, where S={s|s=1, 2, ...s, n}, wherein each element in each set S represents the index number of a specific voxel, and a total of n voxels are included. The set contains the indexes of all voxels in the study area.

[0063] Step S20: Map the drilling position and drilling information of each borehole in the modeling area to the n voxels, complete label assignment for multiple known voxels corresponding to the drilling position, and determine the initial label value of the unknown voxel that does not correspond to the drilling position based on the label value of the known voxel in the cubic learning domain established with the unknown voxel as the center; and construct an initial geological model based on the label value of the known voxel and the initial label value of the unknown voxel.

[0064] In the present invention, known voxels can be understood as voxels with known geotechnical types, and the geotechnical types of known voxels are determined by drilling information. The label values ​​of known voxels are determined and no longer need to be predicted; while voxels that do not correspond to drilling positions are unknown voxels, and predicting their geotechnical types is the purpose of the present invention.

[0065] In the present invention, the total number of label values ​​(ie, the number of types of label values) is the same as the number of rock and soil types. Different rock and soil types have different label values, and voxels with the same label value are considered to be the same soil type.

[0066] It can be understood that the number of geotechnical types is determined based on the drilling information of all boreholes. For the sake of ease of understanding, for example, assuming that the number of boreholes on the geological body is 10, and based on the number of 10 boreholes, it is known that there are 7 geotechnical types, then the geological body in the modeling area is considered to be composed of these 7 geotechnical types, and the geotechnical type of the unknown voxel is one of them.

[0067] In the present invention, the label value of the geotechnical type is represented by a set L, L = {1, 2, 3, ..., l}, where L represents the set of all possible soil type label values, and in this set, there are l types of soil types.

[0068] For ease of understanding, let's take an example and assume that the geotechnical types include rockfill, muck, medium sand, silty clay, fully weathered mixed granite (FWMG), strongly weathered mixed granite (SWMG) and moderately weathered mixed granite (MWMG). Then m=7, and the geotechnical type labels consist of 1, 2, 3, 4, 5, 6 and 7, with each label value representing a geotechnical type.

[0069] Construct a label field x corresponding to n voxels in ascending order of number, where x = (x1, x2, x3, ..., x n ), each x in the label field i All values ​​are taken from the set L.

[0070] It should be noted that determining the geotechnical type / label value of a known voxel corresponding to a borehole location through borehole information is an existing technology, and the present invention does not improve on this.

[0071] In an optional embodiment, in step S20, based on the label value of each known voxel, the initial label value of each unknown voxel is determined in sequence according to a preset order, wherein the step of determining the initial label value of any unknown voxel in the initial geological model includes:

[0072] Step (1): With the unknown voxel i as the center voxel, a cubic learning domain is constructed according to a preset size, and the coordinate information of all voxels in the cubic learning domain is obtained. The voxels in the cubic learning domain are divided into a first voxel with a known label value and a second voxel with an unknown label value, wherein the coordinate information is represented by a coordinate index.

[0073] In the present invention, the length of the cubic learning domain in the X direction is the same as the length of the three-dimensional geological model area in the X direction, the length of the cubic learning domain in the Y direction is the same as the length of the three-dimensional geological model area in the Y direction, and the ratio of the length of the cubic learning domain in the Z direction to the side length of the voxel is 6.

[0074] In the present invention, the coordinate index includes indexes in the X direction, the Y direction, and the Z direction.

[0075] Please refer to Figure 4 , Figure 4 The vertex 11 in the upper left corner is the origin of the physical space coordinates. The first row is row 0, and the indexes increase in the descending rows. The first column is column 0, and the indexes increase in the right columns. The depth of the first row is 0, and the indexes increase inward. The index of each voxel in the X direction, Y direction, and Z direction is 0 or other positive integers.

[0076] In the present invention, the index in the X direction corresponds to the column in the X direction, the index in the Y direction corresponds to the column in the Y direction, and the index in the Z direction corresponds to the row.

[0077] See also Figure 4 , Figure 4 The smallest cube 12 in represents the unknown voxel i / center voxel, Figure 4 The cube 13 shown represents the learning domain of the unknown voxel i / center voxel, Figure 4 The number 14 indicates drilling.

[0078] It can be understood that, in the present invention, the central voxel can be understood as any unknown voxel whose label value is to be determined.

[0079] Step (2): obtaining a label value of each first voxel, and based on the coordinate information and label value of each first voxel, obtaining a first covariance matrix between multiple first voxels with the same label value and a second covariance matrix between multiple first voxels with different label values.

[0080] In this embodiment, the calculation formula of the first covariance matrix is:

[0081] .

[0082] In the above formula: W represents the first covariance matrix, is the coordinate of the first voxel with label value m in The covariance matrix in , yes A subset of , containing only known voxels with the same label m; express exist The proportion of represents the cubic learning domain of unknown voxel i, yes contains only a subset of the first voxels.

[0083] In the present invention, The calculation formula is:

[0084] .

[0085] In the present invention, the calculation formula of the second covariance matrix is:

[0086] .

[0087] In the above formula, B represents the second covariance matrix, and They are and The center of mass, express exist The ratio in , T represents the matrix transpose.

[0088] Step (3): obtaining a local metric of the central voxel based on the first covariance matrix and the second covariance matrix.

[0089] In the present invention, the local metric of the central voxel The calculation formula is:

[0090] .

[0091] Where W is the first covariance matrix, B is the second covariance matrix, Is a tuning parameter, I is a multidimensional identity matrix, and its dimensions are consistent with W and B. You can measure Adjusts from an infinite strip to an ellipsoid, avoiding using known voxels that are far from the center voxel.

[0092] It can be understood that, in the present invention, the calculated local metric of the central voxel is a matrix.

[0093] Step (4): Based on the coordinate information and label value of each first voxel, the coordinate information of the central voxel, and the local metric of the central voxel, the DANN distance between each first voxel and the central voxel is obtained.

[0094] In the present invention, the matrix including the DANN distance between each first voxel and the central voxel is The calculation formula is:

[0095] .

[0096] in, is the coordinate information of the first voxel g, is the coordinate information of the unknown voxel i, c, r, d represent the index in the X, Y, and Z directions respectively, represents the local metric of the central voxel and T represents the matrix transpose.

[0097] In the present invention, V g is a matrix. If the number of voxels is known to be 30, then V g It is a 3-row by 30-column matrix, where each column represents the coordinates of a known voxel.

[0098] Accordingly, is a matrix of multiple rows and one column, where the elements of each row represent the DANN distance between a known voxel and the central voxel. That is, if the number of known voxels is 30, then A matrix of 30 rows and 1 column.

[0099] Step (5): Based on the DANN distance between each first voxel and the central voxel and the local probability calculation formula, the initial label value of the central voxel is obtained.

[0100] It can be understood that the initial label value of the central voxel obtained in step (5) is the initial label value of the unknown voxel i. If the number of unknown voxels is 2000, it is necessary to repeat steps (1) to (5) a total of 2000 times to obtain the initial label value of each unknown voxel. In the present invention, the unknown voxel i specifically refers to the unknown voxel for which the initial label value is to be obtained, which is any one of all the unknown voxels.

[0101] In the present invention, the initial label value of each unknown voxel is determined in sequence according to the number of the unknown voxel, that is, steps (1) to (5) are repeated multiple times, and the initial label value of each unknown voxel is determined in sequence according to the number order from small to large.

[0102] It should be noted that when the unknown voxel is located at the edge of the geological model area, the constructed cubic learning domain does not contain any voxels outside the geological model area, and only the known voxel labels within the learning domain within the geological model area are used when calculating the initial label value.

[0103] Furthermore, step (5) includes:

[0104] Step ①: Based on the DANN distance between the first voxel with the same label value and the central voxel, calculate the DANN distance harmonic mean for evaluating the possibility of selecting each label.

[0105] In the present invention, the formula for calculating the harmonic mean of the DANN distance used to evaluate the likelihood of selecting the label value m is:

[0106] .

[0107] is the harmonic mean distance from unknown voxel i to the first voxel g with the k nearest label value m. The DANN distance used to calculate the harmonic mean distance is based on Select, is the matrix representing the DANN distance between each first voxel and the central voxel, x g is the label value of the first voxel g.

[0108] In the present invention, the k value is determined based on experience, preferably, k=5.

[0109] For ease of understanding, assume that the matrix is ​​30 rows × 1 column, where the DANN distances between the 8 voxels with label value 1 and the unknown voxel i correspond to the matrices From the 1st to the 8th row of , when k=5, the 5 smallest DANN distance values ​​are selected from the data in the 1st to the 8th row to calculate the DANN distance harmonic mean to evaluate the possibility of selecting the label value 1.

[0110] Step ②: Based on the DANN distance harmonic mean used to evaluate the possibility of selecting each label and the local probability calculation formula, the label value probability distribution result of the central voxel is obtained.

[0111] In the present invention, the local probability of selecting label m for unknown voxel i is It can be calculated by the following formula:

[0112] .

[0113] Where, It means that all possible labels in the first voxel are looped over L, where L is the set of label values.

[0114] In the present invention, the local probability calculation formula can be used to obtain the local probability of the central voxel at any label value. After summarizing the local probabilities of the central voxel at each label value, the label value probability distribution result of the central voxel can be obtained.

[0115] Step ③: Based on the probability distribution result of the label value of the central voxel, the initial label value of the central voxel is obtained by random sampling.

[0116] It can be understood that when random sampling is performed based on the probability distribution of label values, the label value with the highest probability has the highest probability of being drawn.

[0117] For ease of understanding, the following example illustrates how to obtain the initial label value for an unknown voxel: assuming the known label values ​​include 1, 2, and 3, the resulting probability distribution of the label value for the unknown voxel is: label value 1, corresponding probability 0.5; label value 2, corresponding probability 0.3; label value 3, corresponding probability 0.2. Based on the probability distribution of the label value for the unknown voxel, a label value is randomly extracted to obtain the initial label value. This label value can be 1, 2, or 3. It can be understood that the label value 1 has the highest probability and is relatively more likely to be randomly extracted.

[0118] Step S30: construct a neighborhood system for each unknown voxel based on the Markov random field, and set spatial constraints in 13 directions for each unknown voxel, and use granularity coefficients β1~β 13 Describes the strength of spatial constraints in each different direction, where the neighborhood system of each unknown voxel consists of other voxels that are adjacent to the unknown voxel and share nodes with it.

[0119] See also Figure 3 , Figure 3 The schematic diagram of the unknown voxel i and its neighborhood system provided by the present invention is as follows: Figure 3 It can be seen that the number of other voxels adjacent to the unknown voxel i and sharing nodes is 26.

[0120] In the present invention, the neighborhood system is defined as ,in represents the set of all voxels that are adjacent to the unknown voxel i and share nodes, thereby forming a local neighborhood system centered on the unknown voxel i; S is a set of n voxel numbers.

[0121] In the present invention, the spatial constraints correspond to 5 intersecting planes in 13 independent directions. These intersecting planes contain all neighbors of the unknown voxel i, that is, the other 26 voxels that are adjacent to the unknown voxel i and share nodes, such as Figure 5 (a). In each 2D voxel grid, there are four spatial constraint parameters that affect four independent directions: 0, π / 2, π / 4, and 3π / 4. Therefore, in a model containing 5 intersecting planes, after excluding overlapping parameters, there are a total of 13 parameters, as shown in Figure 5 (b), that is, β = [β1, β2, β3, β4, β5, β6, β7, β8, β9, β 10 , β 11 , β 12 , β 13], these parameters are called particle size coefficients.

[0122] Specifically: the particle size coefficients corresponding to the four independent directions of the first plane parallel to the X-axis and the Z-axis are β1, β2, β3, and β4, wherein β1 is the 0 degree direction, β2 is the π / 2 direction, β3 is the π / 4 direction, and β3 is the 3π / 4 direction; the particle size coefficients corresponding to the four independent directions of the second plane parallel to the Y-axis and the Z-axis are β2, β5, β6, and β7, wherein β2 is the 0 degree direction, β5 is the π / 2 direction, β6 is the π / 4 direction, and β7 is the 3π / 4 direction; the particle size coefficients corresponding to the four independent directions of the third plane parallel to the X-axis and the Y-axis are β1, β5, β8, and β9, wherein β1 is the 0 degree direction, β5 is the π / 2 direction, β8 is the π / 4 direction, and β9 is the 3π / 4 direction; the particle size coefficients corresponding to the four independent directions of the fourth plane at an angle of 45 degrees to the first plane in the clockwise direction are β2, β9, β1. 10 , β 11 , where β2 is the 0 degree direction, β9 is the π / 2 direction, and β 10 is the π / 4 direction, β 11 The particle size coefficients corresponding to the four independent directions of the fifth plane with an angle of 45 degrees to the first plane in the counterclockwise direction are β2, β8, β 12 , β 13 , where β2 is the 0 degree direction, β8 is the π / 2 direction, and β 12 is the π / 4 direction, β 13 It is in the 3π / 4 direction.

[0123] In order to simplify the model, the model parameters are divided into three groups, where β a = [β1, β5, β8, β9] represents the spatial constraints between unknown voxel i and its horizontally adjacent voxels; β b =[β2] represents the spatial constraint between unknown voxel i and its vertically adjacent voxels; β c =[β3,β4,β6,β7,β 10 , β 11 , β 12 , β 13 ] represents the spatial constraint between the unknown voxel i and its diagonally adjacent voxels (in the π / 4 direction). Spatial constraints can be equal in the same direction, so this simplification is feasible. Accordingly, the granularity coefficient is expressed as: T ={β a , β b , β c}.

[0124] Step S40: Based on the initial label value of each unknown voxel, the label value of each known voxel, and the neighborhood system of each unknown voxel, a prediction result is obtained according to the label probability prediction formula and the update iterative formula of the granularity coefficient. The prediction result includes M sampled label values ​​of each unknown voxel, where M is a positive integer.

[0125] In the present invention, M is greater than or equal to 100.

[0126] In the present invention, the label probability prediction formula includes a local conditional probability distribution calculation formula, a local energy calculation formula, and a potential function calculation formula, wherein:

[0127] The local conditional probability distribution calculation formula is:

[0128] .

[0129] The local energy calculation formula is:

[0130] .

[0131] The potential function calculation formula is:

[0132] .

[0133] In the label probability prediction formula, is the local conditional probability distribution at the unknown voxel i; q is and neighborhood system energy The weight coefficient of It is a potential function defined solely on the unknown voxel i, which represents the preference for selecting different labels on the unknown voxel i; represents local energy; x i represents a specific label at the unknown voxel i, Represents any possible label at the unknown voxel i circulating on L, where L is the set of label values; is the potential function, represents the neighborhood system of unknown voxel i; β d ∈β={β1,β2,β3,β4,β5,β6,β7,β8,β9,β 10 , β 11 , β 12 , β 13}, x j Represents a specific label of other voxels j in the neighborhood system of unknown voxel i.

[0134] In the present invention, HMD is integrated into local energy to maintain The local label preference is characterized and the neighborhood system is further updated through random simulation.

[0135] Preferably, q=4 in the local conditional probability distribution calculation formula.

[0136] The present invention also includes the steps of setting multiple q values ​​and obtaining the optimal q value among the multiple q values. When performing the q value test, multiple original three-dimensional geological bodies with complete geological information are used, some drill holes are extracted from them, and simulation is performed. The geological information finally obtained is compared with the original geological information, and the q value with the highest accuracy is taken as the optimal value.

[0137] In the present invention, the update iterative formula of the granularity coefficient includes a likelihood function calculation formula and an objective function calculation formula, wherein:

[0138] The likelihood function calculation formula is:

[0139] .

[0140] The objective function calculation formula is:

[0141] .

[0142] In the update iterative formula of the granularity coefficient, β represents the granularity coefficient, It represents the prior distribution of β, x unknown Represents the label value of the unknown voxel, x BH represents the label value of a known voxel, represents the neighborhood system of unknown voxel i, It represents the probability of unknown voxel i taking each label value based on the voxel label value and β value in the neighborhood system of unknown voxel i.

[0143] In an optional implementation, step S40 includes:

[0144] Step S401: Based on the preset β1~β 13 Mean and standard deviation, establish the relationship between β1~β 13 The corresponding prior multivariate Gaussian distribution.

[0145] For normal sedimentary stratigraphic structures (i.e., flat ground layers), their distribution usually exhibits horizontal stratification. If expressed using the spatial constraint or granularity coefficient of this model, this means that the coefficient β related to the horizontal direction a will be significantly larger than the coefficients in the other two directions (i.e., β b , β c ). In addition, the model parameters were extensively tested, not only for manually delineated strata, but also for different β T By extracting a large number of virtual known boreholes and T Prior value (i.e. βT The results show that no matter β T The generated strata always obtain satisfactory results regardless of the MCMC tracking performance. It is particularly noteworthy that when the borehole data are sparse, as long as the prior information indicates that the horizontal β a Significantly larger than the vertical β b and β c Therefore, for sites with a limited number of known boreholes, if geological knowledge confirms that the site is a stationary field, β can be effectively used. a , β b and β c The relationship between setting β T The prior value of β a Set to much higher than β b and β c , in order to conduct a more informed random simulation. Accordingly, in this stationary field case, the mean of the prior parameter is set to β T = {5.0, 0.1, 0.1}, indicating that the stratigraphic distribution tends to be horizontally stratified, while the standard deviation σ = 0.1 means that β can fluctuate within a small range during the sampling process.

[0146] In the present invention, β a The mean of β1, β5, β8, and β9 is 5; the mean of β2 is 0.1; c The mean of β3, β4, β6, β7, β 10 , β 11 , β 12 , β 13 The mean value of β1~β 13 The standard deviation of β1 is 0.1, the standard deviation of β2 is 0.1, the standard deviation of β3 is 0.1, the standard deviation of β4 is 0.1... 13 The standard deviation is 0.1.

[0147] Step S402: Based on the label value of the known voxel, the sampled label value of the unknown voxel obtained in the previous iteration and the prior multivariate Gaussian distribution, according to the update iteration formula of the granularity coefficient, obtain β1~β 13 The sampling value of the current iteration is β1~β 13 The corresponding objective function value is greater than β1~β obtained in the last iteration 13 The corresponding objective function value.

[0148] The step S402 includes:

[0149] Step (a) uses the Metropolis-Hastings method to sample from the prior multivariate Gaussian distribution to obtain β1~β 13 The current sample value of .

[0150] Step (b), based on the label value of the known voxel, the sampled label value of the unknown voxel obtained from the last iteration, the prior multivariate Gaussian distribution, and β1~β 13 The current sampling value of β1~β is obtained according to the likelihood function calculation formula and the objective function calculation formula. 13 The current objective function value corresponding to the current sampling value of .

[0151] It should be noted that, when conducting β1~β 13 In the first iteration update, the label value of the unknown voxel used is the initial label value of the unknown voxel; 13 In the second iteration update, the label value of the unknown voxel used is the sampled label value of the unknown voxel obtained in the first iteration update; 13 During the third iterative update, the label value of the unknown voxel used is the sampled label value of the unknown voxel obtained by the second iterative update... and so on.

[0152] Step (c) Compare the current objective function value and the β1~β obtained in the last iteration 13 The corresponding objective function values ​​are compared. When the current objective function value is greater than β1~β obtained in the previous iteration, 13 When the corresponding objective function value is 13 The current sampling value is used as the β1~β obtained in the current iteration 13 The sampling value of the current objective function is less than or equal to the β1~β obtained in the previous iteration. 13 When the corresponding objective function value is reached, return to step (a) and resample.

[0153] Step S403: β1~β2 obtained in the current iteration 13 The sampling value and the label probability prediction formula are used to calculate the local conditional probability distribution of each unknown voxel obtained in the current iteration, and according to the local conditional probability distribution of each unknown voxel, the sampling label value of each unknown voxel obtained in the current iteration is obtained.

[0154] In the present invention, a sampled label value is randomly obtained according to the local conditional probability distribution of each unknown voxel.

[0155] For ease of understanding, the following example illustrates how to obtain sampling label values ​​for unknown voxels: assuming that the types of rock and soil include A, B, and C, the local conditional probability distribution of the unknown voxels is: soil type A, corresponding probability 0.4; soil type B, corresponding probability 0.5; soil type C, corresponding probability 0.1. Based on the local component probability distribution of the unknown voxels, a label value is randomly extracted to obtain a sampling label value. This label value can be the label value corresponding to soil type A, soil type B, or soil type C. It can be understood that soil type B has the highest probability and is relatively more likely to be randomly extracted.

[0156] Step S404: Repeat steps S402 and S403 multiple times until β1~β 13 Converges and in β1~β 13 After convergence, continue to repeat step S402 and step S403 for a preset number of M times, and set β1~β 13 After convergence, the M sampled label values ​​of each unknown voxel obtained by M iterations are used as the prediction results.

[0157] In the present invention, β1~β 13 Convergence can be understood as β1~β 13 Reaching stability, that is, the β1~β obtained twice before and after 13 The absolute value of the difference is less than or equal to the preset tolerance threshold.

[0158] In the present invention, it is determined that β1~β 13 Whether convergence occurs belongs to the existing technology and will not be elaborated here.

[0159] For ease of understanding, let's take an example. Assume that in step S404, a total of 200 iterations are performed. When the number of iterations is 40, β1~β 13 Convergence is equivalent to M=160, that is, after the iteration is completed, the sampling label values ​​of 160 unknown voxels obtained from 41 iterations to 200 iterations are obtained.

[0160] In the present invention, each complete iteration includes obtaining β1~β 13 The sampled values ​​of and the sampled label values ​​of each unknown voxel.

[0161] It should be noted that in β1~β 13 Before convergence, the sampled label values ​​of the unknown voxels are discarded. 13 After convergence, the sampled label values ​​of the unknown voxels are used to calculate the edge probability, and the β1~β 13 The sampling value is used to update the sampling label value of the unknown voxel in the next iteration and to judge whether it has converged.

[0162] It can be understood that when performing iterative operations, in the first iteration, a specific label of the unknown voxel i and voxel j is the initial label, and in the second iteration and subsequent iterations, a specific label of the unknown voxel i and voxel j is the sampled label value obtained in the previous iteration, and in each iteration, the label value of the known voxel is the label value determined based on the drilling information.

[0163] Step S50 : Based on the M sampled label values ​​of each unknown voxel, obtain the marginal probability of each unknown voxel, and based on the marginal probability of each unknown voxel, obtain the predicted label value of each unknown voxel.

[0164] Step 5.1: Based on the M sampled label values ​​of each unknown voxel, count the number of each label value in the M sampled label values ​​of each unknown voxel.

[0165] For ease of understanding, let's take an example. Suppose the number of sampled label values ​​of the unknown voxel is 50 (i.e., M=50), and the label values ​​corresponding to the soil types include 1, 2, and 3, a total of three types. Then calculate the number of sampled label values ​​with label values ​​of 1, 2, and 3 respectively.

[0166] Step 5.2: Calculate the ratio of the number of each label value of each unknown voxel to the total number of sampled label values ​​M to obtain the probability of each label value.

[0167] For ease of understanding, let's take an example. Assume that there are 25 sampled label values ​​with a label value of 1, 15 sampled label values ​​with a label value of 2, and 10 sampled label values ​​with a label value of 3. Then, calculate the ratio of the number of sampled label values ​​with label values ​​of 1, 2, and 3 to the total number of sampled label values, 50, to obtain the probability of each label value. The probability of label value 1 is 50%, the probability of label value 2 is 30%, and the probability of label value 3 is 20%.

[0168] Step 5.3: Summarize the probability of each unknown unit cell at each label value to obtain the edge probability of each unknown voxel.

[0169] When the probability of the label value 1 of the unknown voxel is 50%, the probability of the label value 2 is 30%, and the probability of the label value 3 is 20%, then the marginal probabilities of the unknown voxel are: label value 1, 50%; label value 2, 30%, label value 3, 20%.

[0170] Step 5.4: Based on the marginal probability of each unknown voxel, obtain the predicted label value of each unknown voxel.

[0171] In the present invention, the predicted label value of the unknown voxel has the maximum label probability.

[0172] For ease of understanding, let's take an example. Assuming there are three types of geotechnical soil, there are 3 corresponding label values. The marginal probability may be: the probability of geotechnical soil type A (label value 1) is 0.5, the probability of geotechnical soil type B (label value 2) is 0.3, and the probability of geotechnical soil type C (label value 3) is 0.2. Then the predicted label value of the voxel is 1, and its geotechnical soil type is A. That is, the label value with the maximum probability is used as the predicted label value of the unknown voxel.

[0173] Step S60: updating the initial geological model based on the predicted label value of each unknown voxel to construct a three-dimensional random geological model.

[0174] It can be understood that the predicted label value of each unknown voxel can determine the geotechnical type of each unknown voxel. In this way, after the geotechnical type of each voxel is determined, a three-dimensional random geological model can be constructed.

[0175] Step S70: Based on the calculation formula of the edge probability and information entropy of each unknown voxel, obtain the information entropy of each unknown voxel, and obtain the stratum uncertainty quantification result corresponding to the three-dimensional random geological model according to the information entropy of each unknown voxel.

[0176] In the present invention, the calculation formula of information entropy is:

[0177] .

[0178] Among them, E(x i ) is information entropy, p (m) (x i ) represents the marginal probability of soil type m in unknown voxel i, and L is the set of label values ​​of soil type. i ) value, the more difficult it is to determine the soil type of the unknown voxel i, that is, the higher the uncertainty of the voxel. Therefore, information entropy provides an effective measurement tool for quantifying uncertainty without assuming any specific functional form or statistical model for probability distribution.

[0179] By calculating the entropy value for each voxel, an entropy distribution map for the entire study area can be generated. This entropy distribution map can reveal areas of high uncertainty in the stratigraphic estimate, helping to better understand the model's performance. In areas of high uncertainty, the calculated information entropy provides an intuitive tool for assessing and visualizing stratigraphic estimate quality. For example, if drillhole data is available, this data can be compared with the model estimate to assess model performance and identify areas with insufficient information or complex geological conditions. Therefore, information entropy not only provides a concise and effective means of quantifying voxel uncertainty but also lays the foundation for further geological analysis.

[0180] It should be noted that, in the present invention, i represents an unknown voxel, j represents other voxels in the neighborhood system of the unknown voxel i, and g represents a known voxel in the cubic learning domain of the unknown voxel i.

[0181] See also Figure 6 According to a second aspect of the present invention, a three-dimensional random geological model construction device 100 is provided. The construction device 100 includes: a first construction module 101 for constructing a three-dimensional geological model region based on a geological body in a modeling region, and discretizing the geological model region into n voxels, and numbering each voxel, wherein the modeling region includes a plurality of boreholes opened in the geological body; a second construction module 102 for mapping the drilling position and drilling information of each borehole in the modeling region to the n voxels, completing label assignment for a plurality of known voxels corresponding to the drilling position, and determining the initial label value of the unknown voxel that does not correspond to the drilling position based on the label value of the known voxel in the cubic learning domain established with the unknown voxel as the center, and constructing an initial geological model based on the label value of the known voxel and the initial label value of the unknown voxel; a third construction module 103 for constructing a neighborhood system for each unknown voxel based on a Markov random field, and setting spatial constraints in 13 directions for each unknown voxel, and at the same time using granularity coefficients β1~β 13 Describe the spatial constraint strength in each different direction, wherein the neighborhood system of each unknown voxel consists of other voxels that are adjacent to the unknown voxel and share nodes; a first acquisition module 104 is used to obtain a prediction result based on the initial label value of each unknown voxel, the label value of each known voxel, and the neighborhood system of each unknown voxel according to a label probability prediction formula and an update iterative formula of the granularity coefficient, wherein the prediction result includes M sampled label values ​​of each unknown voxel, where M is a positive integer; a second acquisition module 105 is used to obtain the edge probability of each unknown voxel based on the M sampled label values ​​of each unknown voxel, and obtain the predicted label value of each unknown voxel based on the edge probability of each unknown voxel; an updating module 106 is used to update the initial geological model based on the predicted label value of each unknown voxel to construct a three-dimensional random geological model.

[0182] In an optional embodiment, the construction device also includes: an uncertainty quantification module 107, which is used to obtain the information entropy of each unknown voxel based on the calculation formula of the edge probability and information entropy of each unknown voxel, and obtain the stratigraphic uncertainty quantification result corresponding to the three-dimensional random geological model based on the information entropy of each unknown voxel.

[0183] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the working process of the construction device described above can refer to the contents of the aforementioned method embodiment and will not be repeated here.

[0184] Example 1

[0185] Construct a three-dimensional random geological model according to the construction method described above

[0186] 1.1 Constructing the original stratum

[0187] A Figure 7 The synthetic formation shown in (a) is used as the original formation: 1) Four planes are first generated using a Gaussian random field to divide the modeling domain into several regions; 2) These regions are then filled with three different labels. The modeling domain dimensions are 100 × 50 × 100 (number of voxels). 0.5% of the columns in the original formation domain are randomly selected as known boreholes (i.e., 50 boreholes), as shown in the following example: Figure 7 (b) is used to estimate the stratigraphic uncertainty for the entire region.

[0188] It should be noted that the 100 in the modeling size can be understood as 100 voxels in the X direction, 50 voxels in the Z direction, and 100 voxels in the Y direction. Assuming that each voxel is a cubic voxel, the side length of each voxel can be 100m, 200m, etc., depending on the actual geological body setting.

[0189] 1.2 Initial formation configuration for sampling

[0190] Based on the formation information determined by the known borehole, the formation information of the unknown voxels in the original formation is determined according to the method for obtaining the initial label value of the unknown voxel provided by the present invention. The initial formation configuration is as follows: Figure 8 As shown, this configuration is consistent with the original profile ( Figure 7 (a) has some similarity, but there is more noise in the initial formation configuration.

[0191] 1.3 Simulation Results

[0192] In this embodiment, β a The mean of β1, β5, β8, and β9 is 5; the mean of β2 is 0.1; c The mean of β3, β4, β6, β7, β 10 , β 11 , β 12 , β 13 The mean value of β1~β 13 The standard deviation of β1 is 0.1, the standard deviation of β2 is 0.1, the standard deviation of β3 is 0.1, the standard deviation of β4 is 0.1... 13 The standard deviation is 0.1.

[0193] The MAP of the formation is simulated using the label values ​​determined by known boreholes and the prior distribution established by the above mean and standard deviation, such as Figure 9 (a) The accuracy of MAP estimation reaches 0.93, indicating that the proposed 3D MRF model can effectively infer the stratigraphic formation and quantify its uncertainty using sparse borehole data with only 0.5% known voxels. Figure 9 (b) shows the corresponding information entropy (IE), which reflects the uncertainty level of each voxel.

[0194] Example 2

[0195] The simulation effect of the geological model constructed by the construction method of the present invention is verified based on experimental engineering

[0196] At a certain engineering site, a total of 33 boreholes (BH1~BH33) were collected, and the drilling diagram is shown as follows: Figure 10 As shown, the drilling plan is as follows Figure 11 shown.

[0197] The present invention adopts the cross-validation method to use different borehole combinations for geological simulation, that is, one of the boreholes is taken as the verification borehole and the other 32 boreholes are taken as the test boreholes to simulate the stratum. The simulation results of the cross-validation method are shown in Table 1, and the accuracy distribution results are shown in Table 1. Figure 12 As shown in Table 1, the geological model provided by the present invention performs very well in this engineering example. The accuracy of most boreholes is above 0.80, and the accuracy of BH07, BH20, BH21, BH23, and BH29 is as high as 0.9, indicating that the model has high reliability in most cases.

[0198] Table 1 Cross-validation simulation results

[0199]

[0200] The above content is a further detailed description of the present invention in conjunction with specific preferred embodiments, and the specific implementation of the present invention should not be considered to be limited to these descriptions. For those skilled in the art of the present invention, several simple deductions and substitutions can be made without departing from the concept of the present invention, and all of these should be considered to fall within the scope of protection of the present invention.

Claims

1. A method for constructing a three-dimensional random geological model based on machine learning, characterized in that: The construction method comprises the following steps: Step S10: constructing a three-dimensional geological model region based on the geological body of the modeling region, discretizing the geological model region into n voxels, and numbering each voxel, wherein the modeling region includes a plurality of boreholes opened in the geological body; Step S20: mapping the drilling position and drilling information of each borehole in the modeling area to the n voxels, assigning labels to a plurality of known voxels corresponding to the drilling positions, and determining initial label values ​​of unknown voxels that do not correspond to the drilling positions based on the label values ​​of known voxels in a cubic learning domain established with the unknown voxels as the center; and constructing an initial geological model based on the label values ​​of the known voxels and the initial label values ​​of the unknown voxels. In step S20, based on the label value of each known voxel, the initial label value of each unknown voxel is determined in sequence according to a preset order, wherein the step of determining the initial label value of any unknown voxel in the initial geological model includes: Step (1): With the unknown voxel i as the center voxel, a cubic learning domain is constructed according to a preset size, and coordinate information of all voxels in the cubic learning domain is obtained. The voxels in the cubic learning domain are divided into a first voxel with a known label value and a second voxel with an unknown label value, wherein the coordinate information is represented by a coordinate index; Step (2), obtaining a label value of each first voxel, and based on the coordinate information and label value of each first voxel, obtaining a first covariance matrix between a plurality of first voxels with the same label value and a second covariance matrix between a plurality of first voxels with different label values; Step (3), obtaining a local metric of the central voxel based on the first covariance matrix and the second covariance matrix; Step (4), based on the coordinate information and label value of each first voxel, the coordinate information of the central voxel, and the local metric of the central voxel, obtaining the DANN distance between each first voxel and the central voxel; Step (5), based on the DANN distance between each first voxel and the central voxel and the local probability calculation formula, obtain the initial label value of the central voxel; Step S30: construct a neighborhood system for each unknown voxel based on the Markov random field, and set spatial constraints in 13 directions for each unknown voxel, and use granularity coefficients β1~β 13 Describe the spatial constraint strength in each different direction, where the neighborhood system of each unknown voxel consists of other voxels that are adjacent to the unknown voxel and share nodes with it; Step S40: Based on the initial label value of each unknown voxel, the label value of each known voxel, and the neighborhood system of each unknown voxel, a prediction result is obtained according to a label probability prediction formula and an update iterative formula of a granularity coefficient, wherein the prediction result includes M sampled label values ​​of each unknown voxel, where M is a positive integer. Step S50: obtaining the edge probability of each unknown voxel based on the M sampled label values ​​of each unknown voxel, and obtaining the predicted label value of each unknown voxel based on the edge probability of each unknown voxel; Step S60: updating the initial geological model based on the predicted label value of each unknown voxel to construct a three-dimensional random geological model.

2. The method for constructing a three-dimensional random geological model based on machine learning according to claim 1, characterized in that: The construction method further includes step S70, obtaining the information entropy of each unknown voxel based on a calculation formula of the edge probability and information entropy of each unknown voxel, and obtaining a quantification result of the formation uncertainty corresponding to the three-dimensional random geological model based on the information entropy of each unknown voxel.

3. The method for constructing a three-dimensional random geological model based on machine learning according to claim 1, characterized in that: The step (5) comprises: Based on the DANN distance between the first voxel with the same label value and the center voxel, the harmonic mean of the DANN distance is calculated to evaluate the likelihood of selecting each label; Obtaining a probability distribution result of the label value of the central voxel based on the DANN distance harmonic mean and the local probability calculation formula for evaluating the possibility of selecting each label; Based on the probability distribution result of the label value of the central voxel, the initial label value of the central voxel is obtained by random sampling.

4. The method for constructing a three-dimensional random geological model based on machine learning according to claim 3, characterized in that: The formula for calculating the harmonic mean of the DANN distances used to assess the likelihood of selecting a label value m is: ； Where, is the harmonic mean distance from unknown voxel i to the first voxel g with the k nearest label value m. The DANN distance used to calculate the harmonic mean distance is based on Select, The matrix representing the DANN distance between each first voxel and the central voxel, x g is the label value of the first voxel g; The local probability calculation formula is: ； Where, represents the local probability of selecting label m for unknown voxel i, It means that all possible labels in the first voxel are looped over L, where L is the set of label values.

5. The method for constructing a three-dimensional random geological model based on machine learning according to claim 4, characterized in that: Local metric of the central voxel The calculation formula is: ； Where W is the first covariance matrix, B is the second covariance matrix, Is a tuning parameter, I is a multidimensional identity matrix, and its dimensions are consistent with W and B; A matrix containing the DANN distances between each first voxel and the central voxel The calculation formula is: ； in, is the coordinate information of the first voxel g, is the coordinate information of the unknown voxel i, c, r, d represent the index in the X, Y, and Z directions respectively, represents the local metric of the central voxel and T represents the matrix transpose.

6. The method for constructing a three-dimensional random geological model based on machine learning according to any one of claims 1 to 5, characterized in that: The label probability prediction formula includes the local conditional probability distribution calculation formula, the local energy calculation formula and the potential function calculation formula, where: The local conditional probability distribution calculation formula is: ； The local energy calculation formula is: ； The potential function calculation formula is: ； In the label probability prediction formula, is the local conditional probability distribution at the unknown voxel i; q is and neighborhood system energy The weight coefficient of It is a potential function defined solely on the unknown voxel i, which represents the preference for selecting different labels on the unknown voxel i; represents local energy; x i represents a specific label at the unknown voxel i, Represents any possible label at the unknown voxel i circulating on L, where L is the set of label values; is the potential function, represents the neighborhood system of unknown voxel i; β d ∈β={β1,β2,β3,β4,β5,β6,β7,β8,β9,β 10 , β 11 , β 12 , β 13 }, x j Represents a specific label of other voxels j in the neighborhood system of unknown voxel i; The update iterative formula of the granularity coefficient includes a likelihood function calculation formula and an objective function calculation formula, wherein: The likelihood function calculation formula is: ； The objective function calculation formula is: ； In the update iterative formula of the granularity coefficient, β represents the granularity coefficient, represents the prior distribution of β, x unknown Represents the sampled label value of the unknown voxel, x BH represents the label value of a known voxel, represents the neighborhood system of unknown voxel i, It represents the probability of unknown voxel i taking each label value based on the voxel label value and β value in the neighborhood system of unknown voxel i.

7. The method for constructing a three-dimensional random geological model based on machine learning according to claim 6, characterized in that: The step S40 includes: Step S401: Based on the preset β1~β 13 Mean and standard deviation, establish the relationship between β1~β 13 The corresponding prior multivariate Gaussian distribution; Step S402: Based on the label value of the known voxel, the sampled label value of the unknown voxel obtained in the previous iteration and the prior multivariate Gaussian distribution, according to the update iteration formula of the granularity coefficient, obtain β1~β 13 The sampling value of the current iteration is β1~β 13 The corresponding objective function value is greater than β1~β obtained in the previous iteration 13 The corresponding objective function value; Step S403: β1~β2 obtained in the current iteration 13 The sampling value and the label probability prediction formula are used to calculate the local conditional probability distribution of each unknown voxel obtained in the current iteration, and according to the local conditional probability distribution of each unknown voxel, obtain the sampling label value of each unknown voxel obtained in the current iteration; Step S404: Repeat steps S402 and S403 multiple times until β1~β 13 Converges and in β1~β 13 After convergence, continue to repeat step S402 and step S403 for a preset number of M times, and set β1~β 13 After convergence, the M sampled label values ​​of each unknown voxel obtained by M iterations are used as the prediction results.

8. The method for constructing a three-dimensional random geological model based on machine learning according to claim 7, characterized in that: The step S402 includes: Step (a) uses the Metropolis-Hastings method to sample from the prior multivariate Gaussian distribution to obtain β1~β 13 The current sample value of Step (b), based on the label value of the known voxel, the sampled label value of the unknown voxel obtained from the last iteration, the prior multivariate Gaussian distribution, and β1~β 13 The current sampling value of β1~β is obtained according to the likelihood function calculation formula and the objective function calculation formula. 13 The current objective function value corresponding to the current sampling value of ; Step (c) Compare the current objective function value and the β1~β obtained in the last iteration 13 The corresponding objective function values ​​are compared. When the current objective function value is greater than β1~β obtained in the previous iteration, 13 When the corresponding objective function value is 13 The current sampling value is used as the β1~β obtained in the current iteration 13 The sampling value of the current objective function is less than or equal to the β1~β obtained in the previous iteration. 13 When the corresponding objective function value is reached, return to step (a) and resample.

9. A device for constructing a three-dimensional random geological model based on machine learning, characterized in that: The construction device comprises: A first construction module is configured to construct a three-dimensional geological model region based on a geological body in a modeling region, discretize the geological model region into n voxels, and number each voxel, wherein the modeling region includes a plurality of boreholes opened in the geological body; a second construction module, configured to map the drilling position and drilling information of each borehole in the modeling area to the n voxels, assign labels to a plurality of known voxels corresponding to the drilling positions, determine initial label values ​​of unknown voxels that do not correspond to the drilling positions based on the label values ​​of known voxels in a cubic learning domain established with the unknown voxels as the center, and construct an initial geological model based on the label values ​​of the known voxels and the initial label values ​​of the unknown voxels; The second construction module is further configured to determine, based on the label value of each known voxel, an initial label value of each unknown voxel in sequence according to a preset order, wherein the step of determining the initial label value of any unknown voxel in the initial geological model includes: Step (1): With the unknown voxel i as the center voxel, a cubic learning domain is constructed according to a preset size, and coordinate information of all voxels in the cubic learning domain is obtained. The voxels in the cubic learning domain are divided into a first voxel with a known label value and a second voxel with an unknown label value, wherein the coordinate information is represented by a coordinate index; Step (2), obtaining a label value of each first voxel, and based on the coordinate information and label value of each first voxel, obtaining a first covariance matrix between a plurality of first voxels with the same label value and a second covariance matrix between a plurality of first voxels with different label values; Step (3), obtaining a local metric of the central voxel based on the first covariance matrix and the second covariance matrix; Step (4), based on the coordinate information and label value of each first voxel, the coordinate information of the central voxel, and the local metric of the central voxel, obtaining the DANN distance between each first voxel and the central voxel; Step (5), based on the DANN distance between each first voxel and the central voxel and the local probability calculation formula, obtain the initial label value of the central voxel; The third building block is used to construct a neighborhood system for each unknown voxel based on the Markov random field, and set spatial constraints in 13 directions for each unknown voxel, while using the granularity coefficients β1~β 13 Describe the spatial constraint strength in each different direction, where the neighborhood system of each unknown voxel consists of other voxels that are adjacent to the unknown voxel and share nodes with it; A first acquisition module is configured to acquire a prediction result based on the initial label value of each unknown voxel, the label value of each known voxel, and the neighborhood system of each unknown voxel according to a label probability prediction formula and an update iterative formula of a granularity coefficient, wherein the prediction result includes M sampled label values ​​of each unknown voxel, where M is a positive integer; a second acquisition module, configured to acquire an edge probability of each unknown voxel based on the M sampled label values ​​of each unknown voxel, and acquire a predicted label value of each unknown voxel based on the edge probability of each unknown voxel; The updating module is used to update the initial geological model based on the predicted label value of each unknown voxel to construct a three-dimensional random geological model.

Images