Method and device for constructing three-dimensional digital geological model

Abstract

The invention provides a method and a device for constructing a three-dimensional digital geological model. The construction method comprises the following steps: constructing a geologic model area range and discretizing the geologic model area range into n unit lattices; mapping the drilling data to n unit lattices to construct an initial geological model; constructing a neighborhood system of each unit lattice, and setting spatial constraints in seven directions; based on the label value of each unit lattice in the initial geologic model, obtaining a prediction result comprising M sampling label values of each unknown unit lattice according to a label probability prediction formula and an updating iteration formula of a granularity coefficient; based on the prediction result, obtaining the marginal probability of each unknown unit lattice and the prediction label value of each unknown unit lattice; and updating the initial geologic model based on the predicted label value of each unknown unit lattice, and constructing a three-dimensional digital geologic model. The three-dimensional digital geologic model constructed by adopting the construction method is accurate in simulation result and can represent the uncertainty of the stratum at the same time.

Description

Technical Field

[0001] The present invention belongs to the field of geotechnical engineering technology, and in particular relates to a method and device for constructing a three-dimensional digital geological model. 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 digital geological model construction method based on Markov random field theory and Bayesian machine learning. The construction method uses seven particle size coefficients to describe spatial correlation. The mean and standard deviation of each particle size coefficient initially obtained by iterative simulation of a two-dimensional geological profile are used as the prior distribution. Subsequently, Bayesian machine learning is used to further refine the method. The iterative convergence is fast, and the constructed three-dimensional digital 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 digital geological model, the construction method comprising the following steps: step S10, constructing a three-dimensional geological model region based on a geological body in a modeling area, and discretizing the geological model region into n unit lattices, and numbering each unit lattice, 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 unit lattices, completing initial assignment for the n unit lattices, and constructing an initial geological model, wherein the label values ​​of multiple known unit lattices corresponding to the drilling position in the initial geological model are determined based on the drilling information, and the initial label values ​​of other unknown unit lattices in the initial geological body model that do not correspond to the drilling position are randomly generated based on the total number of label values ​​determined based on the drilling information; step S30, constructing a neighborhood system for each unit lattice based on a Markov random field, and setting spatial constraints in seven directions for each unit lattice, and describing the spatial constraint strength in each different direction with granularity coefficients β1 to β7, wherein the neighborhood system of each unit lattice is composed of The unit lattice is composed of other unit lattices that are adjacent to the unit lattice and share nodes; step S40, based on the initial label value of each unknown unit lattice, the label value of each known unit lattice and the neighborhood system of each unit lattice, according to the label probability prediction formula and the update iterative formula of the granularity coefficient, obtaining a prediction result, the prediction result including M sampled label values ​​of each unknown unit lattice, M is a positive integer, wherein: the step of obtaining the prediction result includes: first, based on a two-dimensional geological profile, obtaining the mean and standard deviation of each granularity coefficient through multiple simulation iterations, and then, based on the mean and standard deviation of each granularity coefficient and the initial geological model, obtaining the prediction result through multiple iterations, the two-dimensional geological profile including a calculated profile of the YZ plane and a calculated profile of the XZ plane; step S50, based on the M sampled label values ​​of each unknown unit lattice, obtaining the edge probability of each unknown unit lattice, and based on the edge probability of each unknown unit lattice, obtaining the predicted label value of each unknown unit lattice; step S60, updating the initial geological model based on the predicted label value of each unknown unit lattice to construct a three-dimensional digital geological model.

[0007] In a second aspect, the present invention provides a device for constructing a three-dimensional digital geological model, the device comprising: a first construction module 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 unit lattices, and numbering each unit lattice, wherein the modeling region includes multiple boreholes opened in the geological body; a second construction module for mapping the drilling position and drilling information of each borehole in the modeling region to the n unit lattices, completing initial assignment of the n unit lattices, and constructing an initial geological model, wherein the label values ​​of multiple known unit lattices corresponding to the drilling position in the initial geological model are determined based on the drilling information, and the initial label values ​​of other unknown unit lattices in the initial geological body model that do not correspond to the drilling position are randomly generated based on the label values ​​of the multiple known unit lattices; a third construction module for constructing a neighborhood system for each unit lattice based on a Markov random field, and setting spatial constraints in seven directions for each unit lattice, while describing the spatial constraint strength in each different direction with granularity coefficients β1 to β7, wherein the neighborhood system of each unit lattice is composed of the unit lattice and the unit lattice. The method comprises the following steps: a first acquisition module for acquiring a prediction result based on the initial label value of each unknown unit cell, the label value of each known unit cell, and the neighborhood system of each unit cell, 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 unit cell, where M is a positive integer, and the step of acquiring the prediction result includes: first, based on a two-dimensional geological profile, acquiring a mean and a standard deviation of each granularity coefficient through multiple simulation iterations, and then, based on the mean and standard deviation of each granularity coefficient and an initial geological model, acquiring the prediction result through multiple iterations, wherein the two-dimensional geological profile includes a calculated profile of a YZ plane and a calculated profile of an XZ plane; a second acquisition module for acquiring an edge probability of each unknown unit cell based on the M sampled label values ​​of each unknown unit cell, and acquiring a predicted label value of each unknown unit cell based on the edge probability of each unknown unit cell; and an update module for updating the initial geological model based on the predicted label value of each unknown unit cell to construct a three-dimensional digital geological model.

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

[0009] 1. Based on Markov random field theory and Bayesian machine learning, the present invention provides an efficient three-dimensional digital geological model construction method. This construction method uses seven particle size coefficients to describe spatial correlation. The mean and standard deviation of each particle size coefficient obtained by iterative simulation of a two-dimensional geological profile are initially used as the prior distribution. Subsequently, Bayesian machine learning is used to further refine the accuracy and iterative convergence. The constructed three-dimensional digital geological model has accurate simulation results and can characterize the uncertainty of the strata. Compared with the two-dimensional model, the three-dimensional digital model constructed by the present invention is more intuitive, accurate and comprehensive in expressing spatial geological bodies. At the same time, the geological information in the three-dimensional digital model is richer, more realistic and reliable, facilitating more in-depth geological analysis.

[0010] 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.

[0011] 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

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

[0013] Figure 2 A schematic flow chart of the steps for obtaining the mean and standard deviation of each particle size coefficient through multiple simulation iterations based on a two-dimensional geological profile provided by the present invention;

[0014] Figure 3 A schematic flow chart of the steps for obtaining prediction results through multiple iterations based on the mean and standard deviation of each particle size coefficient and the initial geological model provided by the present invention;

[0015] Figure 4 A schematic diagram of a unit cell i and its neighborhood system provided by the present invention;

[0016] Figure 5 A diagram showing the spatial constraints set in seven directions for a unit cell provided by the present invention;

[0017] Figure 6 The three-dimensional geological model and drilling schematic diagram constructed in Example 1 of the present invention, wherein: Figure 6 (a) is a three-dimensional geological model diagram generated by the Gibbs sampler simulation in Example 1; Figure 6 (b) in the figure shows the drilling diagram of randomly selecting 5% of the columns as known drilling holes;

[0018] Figure 7 for Figure 6 (a) A top view of the three-dimensional geological model diagram;

[0019] Figure 8 For Figure 7 Schematic diagram of the calculated section of the YZ plane corresponding to the YZ projection line;

[0020] Figure 9 For Figure 7 Schematic diagram of the calculated cross section of the XZ plane corresponding to the XZ projection line;

[0021] Figure 10 for Figure 8 The trajectories of the posterior β1 to β4 on the calculated section of the YZ plane are shown;

[0022] Figure 11 for Figure 9 The trajectories of the posterior β2 and posterior β5 to β7 on the calculated section of the XZ plane are shown;

[0023] Figure 12 The three-dimensional digital 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 digital geological model are shown, wherein: Figure 12 (a) is the three-dimensional digital geological model map constructed in Example 1, Figure 12 (b) is the quantification result diagram of stratigraphic uncertainty corresponding to the three-dimensional digital geological model;

[0024] Figure 13 This is a schematic diagram of 25 drill holes provided in Example 2, wherein: Figure 13 (a) is the plan view of the positions of 25 drill holes. Figure 13 (b) is the geological information map of 25 drill holes;

[0025] Figure 14 The three-dimensional digital geological model constructed in Example 2 and the stratum uncertainty quantification result diagram corresponding to the three-dimensional digital geological model, wherein: Figure 14 (a) is the three-dimensional digital geological model map constructed in Example 2. Figure 14 (b) is the quantification result diagram of stratigraphic uncertainty corresponding to the three-dimensional digital geological model;

[0026] Figure 15 A module diagram of a device for constructing a three-dimensional digital geological model provided in another embodiment of the present invention. DETAILED DESCRIPTION

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

[0028] 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 digital geological model is provided. Compared with traditional two-dimensional models, the proposed three-dimensional digital model is more intuitive, accurate, and comprehensive in its representation of spatial geological volumes. Furthermore, the geological information contained in the three-dimensional model is richer, more authentic, and more reliable, facilitating more in-depth geological analysis.

[0029] See also Figures 1 to 3 , the construction method provided by the present invention comprises the following steps:

[0030] 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 unit cells, and number each unit cell. The modeling region includes a plurality of boreholes opened in the geological body.

[0031] In the present invention, each unit cell is in the shape of a cube.

[0032] 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.

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

[0034] In the present invention, each unit cell is numbered to construct a set S, wherein the set S = {i|i = 1, 2, ..., n}, where each i represents the index number of a specific unit cell, and the set includes the indices of all unit cells in the study area.

[0035] Step S20: Mapping the drilling position and drilling information of each borehole in the modeling area to the n unit lattices, completing initial value assignment for the n unit lattices, and constructing an initial geological model, wherein the label values ​​of multiple known unit lattices corresponding to the drilling positions in the initial geological model are determined based on the drilling information, and the initial label values ​​of other unknown unit lattices in the initial geological body model that do not correspond to the drilling positions are randomly generated based on the total number of label values ​​determined based on the drilling information.

[0036] In the present invention, a known unit lattice can be understood as a unit lattice with a known rock and soil type. The rock and soil type of the known unit lattice is determined by drilling information and is known and no longer needs to be predicted. The unit lattice that does not correspond to the drilling position is an unknown unit lattice, and predicting its rock and soil type is the purpose of the present invention.

[0037] 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 unit cells with the same label value are regarded as the same type of soil.

[0038] 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 unit lattice is one of them.

[0039] In the present invention, the label value of the rock and soil type is represented by a set L, L = {1, 2, 3, ..., m}, L represents the set of all possible soil type label values, in which there are m types of soil types.

[0040] For ease of understanding, let us take an example and assume that the geotechnical types include seven types: 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.

[0041] Construct a label field x corresponding to n unit cells in the order of numbering from small to large, wherein the label field x = (x1, x2, x3, ..., x n ), each x in the label field i All values ​​are taken from the set L.

[0042] It should be noted that determining the rock and soil type / label value of a known unit cell corresponding to a drilling position through drilling information is an existing technology, and the present invention does not improve on this.

[0043] In the present invention, the initial label value of each unknown unit cell is randomly generated based on the total number of label values. For ease of understanding, let's assume that the total number of label values ​​is 7, i.e., the label values ​​are 1, 2, 3, 4, 5, 6, and 7. The initial label value of each unknown unit cell is any one of 1 to 7, and each label value has equal probability of being generated.

[0044] Step S30: construct a neighborhood system for each unit cell based on the Markov random field, set spatial constraints in seven directions for each unit cell, and describe the strength of the spatial constraints in each different direction with granularity coefficients β1 to β7, wherein the neighborhood system of each unit cell is composed of other unit cells that are adjacent to the unit cell and share nodes with the unit cell.

[0045] In the present invention, the number of other unit cells adjacent to the unit cell and sharing a node is 26.

[0046] In the present invention, the neighborhood system is defined as in It represents the set of all unit cells adjacent to unit cell i and sharing nodes, thus forming a local neighborhood system centered on unit cell i; S is a set of n cell numbers.

[0047] In the present invention, the first of the seven directions is parallel to the X-axis, and its particle size coefficient is represented by β1; the second direction is parallel to the Z-axis, and its particle size coefficient is represented by β2; the third direction is located on the same straight line as the minor diagonal of a plane parallel to both the X-axis and the Z-axis, and its particle size coefficient is represented by β3; the fourth direction is located on the same straight line as the main diagonal of a plane parallel to both the X-axis and the Z-axis, and its particle size coefficient is represented by β4; the fifth direction is parallel to the Y-axis, and its particle size coefficient is represented by β5; the sixth direction is located on the same straight line as the minor diagonal of a plane parallel to both the Y-axis and the Z-axis, and its particle size coefficient is represented by β6; the seventh direction is located on the same straight line as the main diagonal of a plane parallel to both the Y-axis and the Z-axis, and its particle size coefficient is represented by β7.

[0048] It should be noted that, in the present invention, the main diagonal line of the rectangle is defined as running from the upper left to the lower right with a negative slope, and the secondary diagonal line is defined as running from the upper right to the lower left with a positive slope.

[0049] Please refer to Figure 4 and Figure 5 , Figure 4 Schematic diagram of the unit cell i and its neighborhood system provided by the present invention, Figure 5 The spatial constraints in seven directions are shown.

[0050] Step S40: Based on the initial label value of each unknown unit cell, the label value of each known unit cell, and the neighborhood system of each unit cell, 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 unit cell, where M is a positive integer.

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

[0052] 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, as shown in formulas (1) to (3), specifically as follows:

[0053]

[0054] in, is the local conditional probability distribution of lattice i, represents the local energy, x i represents a specific label at lattice i, x i ' represents any possible label at lattice i, looping over L, where L is the set of label values.

[0055] In the present invention, local energy is calculated using formula (2) and formula (3):

[0056]

[0057] represents the local energy, V i,j (x i ,x j ) is the potential function, represents the neighborhood system of unit cell i, x i represents a specific label at lattice i, x j represents a specific label of other lattices j in the neighborhood system of lattice i, β d ∈β={β1,β2,β3,β4,β5,β6,β7}.

[0058] It can be understood that when the iterative operation is performed below, in the first iteration, a specific label of unit cell i and unit cell j is the initial label (the lattice belongs to an unknown unit cell), and in the second iteration and subsequent iterations, a specific label of unit cell i and unit cell j is the sampled label value obtained in the previous iteration, and in each iteration, the label value of the known unit cell is the label value determined based on the drilling information.

[0059] In the present invention, the update iterative formula of the granularity coefficient includes a likelihood function calculation formula and an objective function calculation formula, as shown in formula (4) and formula (5), where:

[0060] The likelihood function calculation formula is:

[0061]

[0062] The objective function calculation formula is:

[0063] log(target)=log(Prior(β))+log(L(x BH |β,x unknown )) (5).

[0064] In formula (4) and formula (5), β represents the particle size coefficient, Prior(β) represents the prior distribution of β, and x unknown Indicates the label value of the unknown unit cell, x BH Represents the label value of a known unit cell, represents the neighborhood system of unit cell i.

[0065] The step of obtaining the prediction result includes: first, based on the two-dimensional geological profile, obtaining the mean and standard deviation of each particle size coefficient through multiple simulation iterations, and then based on the mean and standard deviation of each particle size coefficient and the initial geological model, obtaining the prediction result through multiple iterations, wherein the two-dimensional geological profile includes the calculated profile of the YZ plane and the calculated profile of the XZ plane.

[0066] In an optional embodiment, the step of obtaining the mean and standard deviation of each particle size coefficient through multiple simulation iterations based on the two-dimensional geological profile includes:

[0067] Step S41: Based on the drilling positions of multiple boreholes in the modeling area, select the estimated section of the YZ plane, and map the drilling information of multiple first boreholes on the estimated section of the YZ plane and multiple second boreholes adjacent to the estimated section of the YZ plane to the estimated section of the YZ plane to obtain the calculated section of the YZ plane.

[0068] In the present invention, a top view of the geological body in the modeling area is first obtained; then the distribution of the boreholes is checked from the top view; and then an estimated section of the YZ plane is selected based on the distribution of the boreholes.

[0069] In the present invention, the estimated section of the YZ plane is selected by those skilled in the art based on experience. In principle, the criterion is that the sum of the number of drill holes on the section and the number of adjacent drill holes is larger.

[0070] In the present invention, the method for determining the second borehole adjacent to the estimated profile of the YZ plane is: the distance between the second borehole and the estimated profile of the YZ plane is within a preset distance range, and it is theoretically believed that the label values ​​of each lattice within the preset distance range and located in the same geological layer are the same.

[0071] When constructing the model, those skilled in the art may determine the second borehole based on experience.

[0072] Step S42: Based on the calculated profile of the YZ plane, β1 to β4 are first converged through multiple iterations. After β1 to β4 converge, the first preset number of iterations K are continued to be performed to obtain K sample values ​​of β1 to β4 obtained by iterating K times after β1 to β4 converge.

[0073] It can be understood that the simulation iterative operation is performed on the calculation section of the YZ plane. Among the seven constraint directions set by the present invention, only the first direction, the second direction, the third direction and the fourth direction are on the calculation section of the YZ plane. Therefore, when the simulation iterative operation is performed based on the calculation section of the surface, the values ​​of β1 to β4 can be determined.

[0074] This step includes:

[0075] Step 4.2.1. Obtain the number of each first unit cell and the number of each second unit cell on the calculated section of the YZ plane, where the first unit cell is a unit cell of an unknown rock and soil type, and the second unit cell is a unit cell of a known rock and soil type.

[0076] Step 4.2.2: Based on the preset means and standard deviations of β1 to β4, establish a first priori multivariate Gaussian distribution corresponding to β1 to β4.

[0077] In the present invention, the preset means of β1 to β4 are all 1 (i.e., β1 is 1, β2 is 1, β3 is 1, and β4 is 1), and the standard deviations of β1 to β4 are all 100 (i.e., the standard deviation of β1 is 100, the standard deviation of β2 is 100, the standard deviation of β3 is 100, and the standard deviation of β4 is 100).

[0078] In the present invention, the standard deviation is set to be relatively large, which can make the range of the first priori multivariate Gaussian distribution relatively large, so that reasonable posterior values ​​can be sampled.

[0079] Step 4.2.3, based on the initial label value of each first unit lattice, the label value of each second unit lattice and the neighborhood system of each unit lattice, perform iterative operation according to the label probability prediction formula and the update iteration formula of the granularity coefficient, first through multiple iterations until β1~β4 converge, and after β1~β4 converge, continue to iterate the first preset number of times K times, and obtain K sampling values ​​of β1~β4 obtained by iterating K times after β1~β4 converge.

[0080] In an optional implementation, K is greater than or equal to 100.

[0081] In the present invention, the convergence of β1 to β4 can be understood as β1 to β4 reaching stability, that is, the absolute value of the difference between β1 to β4 obtained twice before and after is less than or equal to a preset tolerance threshold.

[0082] In the present invention, judging whether β1 to β4 converge is a prior art and will not be described in detail here.

[0083] For ease of understanding, let's take an example. Assume that in step (3), a total of 1000 iterations are performed. When the number of iterations is 200, β1 to β4 converge, which is equivalent to K = 800. That is, after the iteration is completed, 800 sample values ​​of β1 to β4 obtained from 201 iterations to 1000 iterations are obtained.

[0084] In the present invention, each complete iteration includes obtaining the sampled values ​​of β1 to β4 and the sampled label value of each first unit cell.

[0085] It should be noted that before β1~β4 converge, the sampling values ​​of β1~β4 are discarded. After β1~β4 converge, the sampling values ​​of β1~β4 are used to calculate the mean and variance of β1~β4, and the sampling label value of the first unit cell obtained in each iteration is used to update β1~β4 for the next use.

[0086] In an optional implementation, any iterative update of β1 to β4 includes:

[0087] Step (a): use the Metropolis-Hastings method to sample from the first prior multivariate Gaussian distribution to obtain the current sampling values ​​of β1 to β4.

[0088] Step (d) is to obtain the current first objective function value corresponding to the current sampling value of β1 to β4 based on the label values ​​of multiple second unit cells, the sampling label value of the first unit cell obtained from the last iterative update, the first prior multivariate Gaussian distribution and the current sampling values ​​of β1 to β4 according to the likelihood function calculation formula and the objective function calculation formula.

[0089] It should be noted that when performing the first iterative update of β1 to β4, the label value of the first unit lattice used is the initial label value of the first unit lattice. When performing the second iterative update of β1 to β4, the label value of the first unit lattice used is the sampled label value of the first unit lattice obtained in the first iterative update.

[0090] In the present invention, the likelihood function calculation formula is shown in the above formula (4), and the objective function calculation formula is shown in the above formula (5), which will not be repeated here.

[0091] Step (c): compare the current first objective function value with the first objective function values ​​corresponding to β1 to β4 obtained in the previous iteration. When the current first objective function value is greater than the first objective function values ​​corresponding to β1 to β4 obtained in the previous iteration, use the current sampling values ​​of β1 to β4 as the sampling values ​​of β1 to β4 obtained in the current iteration; when the current first objective function value is less than or equal to the first objective function values ​​corresponding to β1 to β4 obtained in the previous iteration, return to step (a) for resampling.

[0092] It should be noted that in the present invention, the first objective function value corresponding to the sampling values ​​of β1~β4 obtained in the current iteration needs to be greater than the first objective function value corresponding to the sampling values ​​of β1~β4 obtained in the previous iteration, that is, each sampling is performed based on the conditional posterior distribution.

[0093] Step (d) obtains the local conditional probability distribution of the first unit cell updated in the current iteration based on the sampled values ​​of β1 to β4 obtained in the current iteration and the label probability prediction formula, and obtains the sampled label value of each first unit cell obtained in the current iteration based on the local conditional probability distribution of each first unit cell. The label probability prediction formula is formula (1) to formula (3) above and will not be repeated here.

[0094] In the present invention, a sample label value is randomly acquired according to the local conditional probability distribution of each first unit cell.

[0095] For ease of understanding, the following example illustrates how to obtain the sampling label value for the first unit cell: assuming that the types of rock and soil include A, B, and C, the local conditional probability distribution for the first unit cell 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 first unit cell, 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.

[0096] In the present invention, the local conditional probability distribution result of each first unit cell is a result obtained by predicting the label values ​​of the unit cells adjacent to it, and the adjacent unit cells are determined based on the set constraint directions. In this embodiment, when performing simulation iterative operations based on the YZ profile of the two-dimensional geological model, the number of unit cells adjacent to the predicted unit cell is 8, and the 8 adjacent unit cells are distributed in the first direction, second direction, third direction and fourth direction set for the predicted unit cell.

[0097] Step S43: Calculate the mean and standard deviation of β1 to β4 based on the K sample values ​​of β1 to β4.

[0098] In the present invention, the mean and standard deviation of the particle size coefficients β1 to β4 refer to the mean and standard deviation of the particle size coefficient β1, the mean and standard deviation of the particle size coefficient β2, the mean and standard deviation of the particle size coefficient β3, and the mean and standard deviation of the particle size coefficient β4.

[0099] In the present invention, the mean value of the particle size coefficient is the arithmetic mean.

[0100] Step S44: Based on the drilling positions of multiple boreholes in the modeling area, select the estimated section of the XZ plane, and map the drilling information of multiple third boreholes on the estimated section of the XZ plane and multiple fourth boreholes adjacent to the estimated section of the XZ plane to the estimated section of the XZ plane to obtain the calculated section of the XZ plane.

[0101] In the present invention, a top view of the geological body in the modeling area is first obtained; then the distribution of the boreholes is checked from the top view; and then an estimated section of the XZ plane is selected based on the distribution of the boreholes.

[0102] In the present invention, the estimated section of the XZ plane is selected by those skilled in the art based on experience. In principle, the criterion is that the sum of the number of drill holes on the section and the number of adjacent drill holes is larger.

[0103] In the present invention, the method for determining the fourth borehole adjacent to the estimated profile of the XZ plane is: the distance between the fourth borehole and the estimated profile of the XZ plane is within a preset distance range, and it is theoretically believed that the label values ​​of each lattice within the preset distance range and located in the same geological layer are the same.

[0104] When constructing the model, those skilled in the art may determine the fourth borehole based on experience.

[0105] Step S45: Based on the calculated profile of the XZ plane and the mean and standard deviation of the particle size coefficient β2, β5~β7 are first converged through multiple iterations. After β5~β7 converge, continue to iterate for a second preset number of times L to obtain L sampling values ​​of β5~β7 obtained by iterating L times after β5~β7 converge.

[0106] It can be understood that in step S45, the simulation iterative operation is performed based on the calculated section of the XZ plane. Among the seven constraint directions set by the present invention, only the second direction, the fifth direction, the sixth direction and the seventh direction are on the calculated section of the XZ plane, and the mean and standard deviation of β2 corresponding to the second direction have been determined in step S42. Therefore, when the simulation iterative operation is performed based on the calculated section of the XZ plane, the values ​​of β5 to β7 are determined.

[0107] This step includes:

[0108] Step 4.5.1. Obtain the number of each third unit cell and the number of each fourth unit cell on the calculated section of the XZ plane, wherein the third unit cell is a unit cell of an unknown rock and soil type, and the fourth unit cell is a unit cell of a known rock and soil type.

[0109] Step 4.5.2: Based on the preset mean and standard deviation parameters of β5 to β7 and the mean and standard deviation of β2, a second priori multivariate Gaussian distribution corresponding to β2 and β5 to β7 is established.

[0110] In the present invention, the preset mean values ​​of β5 to β7 are all 1 (i.e., β5 is 1, β6 is 1, and β7 is 1), and the preset standard deviations of β5 to β7 are all 100 (i.e., the standard deviation of β5 is 100, the standard deviation of β6 is 100, and the standard deviation of β7 is 100).

[0111] In the present invention, the standard deviations of β5 to β7 are set relatively large, which can make the range of the first priori multivariate Gaussian distribution relatively large, so that reasonable posterior values ​​can be sampled.

[0112] The mean and standard deviation of β2 are determined in step S43 and have converged.

[0113] Step 4.5.3, based on the second prior multivariate Gaussian distribution, the initial label value of each third unit lattice, the label value of each fourth unit lattice and the neighborhood system of each unit lattice, perform iterative operation according to the label probability prediction formula and the updated iterative formula of the granularity coefficient, first through multiple iterations until β5~β7 converge, and after β5~β7 converge, continue to iterate the second preset number of times L times, and obtain K sampling values ​​of β5~β7 obtained by iterating L times after β5~β7 converge.

[0114] In the present invention, L is greater than or equal to 100.

[0115] In the present invention, the convergence of β5 to β7 can be understood as β5 to β7 reaching stability, that is, the absolute value of the difference between β5 to β7 obtained twice before and after is less than or equal to a preset tolerance threshold.

[0116] In the present invention, judging whether β5 to β7 are converged belongs to the prior art and will not be described in detail here.

[0117] For ease of understanding, let's take an example. Assume that in step 4.5.3, a total of 1000 iterations are performed. When the number of iterations is 100, β5~β7 converge, which is equivalent to L=900. That is, after the iteration is completed, 900 sample values ​​of β5~β7 obtained from 101 iterations to 1000 iterations are obtained.

[0118] In the present invention, each complete iteration includes obtaining the sampled values ​​of β5 to β7 and the sampled label values ​​of each third unit cell.

[0119] It should be noted that before β5~β7 converge, the sampling values ​​of β5~β7 are discarded. After β5~β7 converge, the sampling values ​​of β5~β7 are used to calculate the mean and variance of β5~β7, and the sampling label value of the third unit cell obtained in each iteration is used to update β5~β7 for the next use.

[0120] In an optional implementation, any iterative update of β5 to β7 includes:

[0121] Step (a): use the Metropolis-Hastings method to sample from the second prior multivariate Gaussian distribution to obtain the current sampling values ​​of β2 and β5~β7.

[0122] Step (b) is based on the label values ​​of multiple fourth unit cells, the sampling label values ​​of the third unit cell obtained from the last iterative update, the second prior multivariate distribution, and the current sampling values ​​of β2 and β5 to β7, and according to the likelihood function calculation formula and the objective function calculation formula, obtain the current second objective function value corresponding to the current sampling value of β5 to β7.

[0123] It should be noted that when performing the first iterative update of β5 to β7, the label value of the third unit lattice used is the initial label value of the third unit lattice; when performing the second iterative update of β5 to β7, the label value of the third unit lattice used is the sampled label value of the third unit lattice obtained in the first iterative update; when performing the third iterative update of β5 to β7, the label value of the third unit lattice used is the sampled label value of the third unit lattice obtained in the second iterative update... and so on.

[0124] In the present invention, the likelihood function calculation formula is shown in the above formula (4), and the objective function calculation formula is shown in the above formula (5), which will not be repeated here.

[0125] Step (c): compare the current second objective function value with the second objective function value corresponding to β5~β7 obtained in the previous iteration. When the current second objective function value is greater than the second objective function value corresponding to β5~β7 obtained in the previous iteration, use the current sampling value of β5~β7 as the sampling value of β5~β7 obtained in the current iteration; when the current second objective function value is less than or equal to the second objective function value corresponding to β5~β7 obtained in the previous iteration, return to step (a) for resampling.

[0126] Step (d): Based on the sampling values ​​of β2 and β5~β7 obtained in the current iteration and the label probability prediction formula, obtain the local conditional probability distribution of the third unit lattice obtained by the current iteration update, and based on the local conditional probability distribution of each third unit lattice, obtain the sampling label value of each third unit lattice obtained in the current iteration.

[0127] The label probability prediction formula is the above formula (1) to formula (3), which will not be repeated here.

[0128] For ease of understanding, the following example illustrates how to obtain the sampling label value for the third unit cell: Assuming that the types of rock and soil include A, B, and C, the local conditional probability distribution for the third unit cell is: Soil type A, corresponding probability 0.6; Soil type B, corresponding probability 0.3; Soil type C, corresponding probability 0.1. Based on the local component probability distribution of the first unit cell, 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.

[0129] In the present invention, the local conditional probability distribution result of each third unit cell is a result obtained by predicting the label values ​​of the unit cells adjacent to it, and the adjacent unit cells are determined based on the set constraint directions. In this embodiment, when performing simulation iterative operations based on the calculation profile of the XZ plane, the number of unit cells adjacent to the predicted unit cell is 8, and the 8 adjacent unit cells are distributed in the second direction, fifth direction, sixth direction and seventh direction set for the predicted unit cell.

[0130] It should be noted that in step 4.5.3, when performing iterative calculations based on the calculated profile of the XZ plane, β2 is sampled based on the calculation needs of the iterative process. Its sampling value does not need to be used to calculate the mean and standard value of β2. This round of iterative process is only for obtaining the mean and standard deviation of β5 to β7.

[0131] Step S46: Calculate the mean and standard deviation of β5 to β7 based on the L sample values ​​of β5 to β7.

[0132] In the present invention, the mean and standard deviation of the particle size coefficients β5 to β7 refer to the mean and standard deviation of the particle size coefficient β5, the mean and standard deviation of the particle size coefficient β6, and the mean and standard deviation of the particle size coefficient β7.

[0133] In the present invention, the mean value of the particle size coefficient is the arithmetic mean.

[0134] In an optional embodiment, the step of obtaining the prediction result through multiple iterations based on the mean and standard deviation of each particle size coefficient and the initial geological model includes:

[0135] Step S401: Use the mean and standard deviation of each particle size coefficient as parameters of a priori multivariate Gaussian distribution to establish a third priori multivariate Gaussian distribution corresponding to β1 to β7.

[0136] In the present invention, the mean and standard deviation of each particle size coefficient obtained by iterative simulation based on the two-dimensional geological profile are used as parameters of the prior multivariate Gaussian distribution to establish a third prior multivariate Gaussian distribution, which can accelerate the convergence speed of the three-dimensional iterative operation.

[0137] Step S402: Based on the label value of the known lattice unit, the sampling label value of the unknown unit lattice obtained in the last iteration, and the third prior multivariate Gaussian distribution, according to the updated iteration formula of the granularity coefficient, the sampling values ​​of β1 to β7 obtained in the current iteration are obtained, and the third objective function value corresponding to β1 to β7 obtained in the current iteration is greater than the third objective function value corresponding to β1 to β7 obtained in the last iteration.

[0138] The step S402 includes:

[0139] Step (1): sampling from the third prior distribution using the Metropolis-Hastings method to obtain current sampling values ​​of β1 to β7.

[0140] Step (2): Based on the label value of the known lattice unit, the sampling label value of the unknown lattice unit obtained by the last iterative update, the third prior multivariate Gaussian distribution, and the current sampling values ​​of β1 to β7, the current third objective function value corresponding to the current sampling value of β1 to β7 is obtained according to the likelihood function calculation formula and the objective function calculation formula.

[0141] It should be noted that when performing the first iterative update of β1 to β7, the label value of the unknown unit lattice used is the initial label value of the unknown unit lattice; when performing the second iterative update of β1 to β7, the label value of the unknown unit lattice used is the sampled label value of the unknown unit lattice obtained in the first iterative update; when performing the third iterative update of β1 to β7, the label value of the unknown unit lattice used is the sampled label value of the unknown unit lattice obtained in the second iterative update... and so on.

[0142] In the present invention, the likelihood function calculation formula is shown in the above formula (4), and the objective function calculation formula is shown in the above formula (5), which will not be repeated here.

[0143] Step (3): compare the current third objective function value with the third objective function values ​​corresponding to β1 to β7 obtained in the previous iteration. When the current third objective function value is greater than the third objective function values ​​corresponding to β1 to β7 obtained in the previous iteration, use the current sampling values ​​of β1 to β7 as the sampling values ​​of β1 to β7 obtained in the current iteration; when the current third objective function value is less than or equal to the third objective function values ​​corresponding to β1 to β7 obtained in the previous iteration, return to step (1) and resample.

[0144] Step S403: Calculate the local conditional probability distribution of each unknown lattice unit obtained in the current iteration based on the sampling values ​​of β1 to β7 obtained in the current iteration and the label probability prediction formula, and obtain the sampling label value of each unknown unit lattice obtained in the current iteration based on the local conditional probability distribution of each unknown unit lattice.

[0145] The label probability prediction formula is the above formula (1) to formula (3), which will not be repeated here.

[0146] The method for obtaining the sampling label value of each unknown unit cell is the same as the method for obtaining the sampling label value of the first unit cell and the sampling label value of the third unit cell above. Both are randomly obtained and related to the conditional probability distribution results, which will not be repeated here.

[0147] Step S404, repeat step S402 and step S403 multiple times until β1~β7 converge, and after β1~β7 converge, continue to repeat step S402 and step S403 a third preset number M times, and use the M sampling label values ​​of each unknown unit cell obtained by iterating M times after β1~β7 converge as the prediction result.

[0148] In the present invention, the convergence of β1 to β7 can be understood as β1 to β7 reaching stability, that is, the absolute value of the difference between β1 to β7 obtained twice before and after is less than or equal to a preset tolerance threshold.

[0149] In the present invention, judging whether β1 to β7 are converged belongs to the prior art and will not be described in detail here.

[0150] For ease of understanding, let's take an example and assume that in step S404, a total of 200 iterations are performed. When the number of iterations is 40, β1 to β7 converge, which is equivalent to M = 160. That is, after the iteration is completed, the sampling label values ​​of 160 unknown unit cells obtained from 41 iterations to 200 iterations are obtained.

[0151] In the present invention, each complete iteration includes obtaining the sampled values ​​of β1 to β7 and the sampled label value of each unknown unit cell.

[0152] It should be noted that before β1~β7 converge, the sampling label values ​​of the unknown unit lattice are discarded. After β1~β7 converge, the sampling label values ​​of the unknown unit lattice are used to calculate the edge probability, and the sampling values ​​of β1~β7 obtained in each iteration are used to update the sampling label values ​​of the unknown unit lattice in the next iteration and to determine whether it converges.

[0153] Step S50 : Based on the M sampled label values ​​of each unknown unit cell, obtain the edge probability of each unknown unit cell, and based on the edge probability of each unknown unit cell, obtain the predicted label value of each unknown unit cell.

[0154] Step (1): based on the M sampled label values ​​of each unknown unit cell, count the number of each label value in the M sampled label values ​​of each unknown unit cell.

[0155] For ease of understanding, let's take an example and assume that the number of sampling label values ​​of the unknown unit cell is 50 (i.e., M=50), and the label values ​​corresponding to the soil types include three types: 1, 2, and 3. Then calculate the number of sampling label values ​​with label values ​​of 1, 2, and 3 respectively.

[0156] Step (2) calculates the ratio of the number of each label value of each unknown unit cell to the total number M of sampled label values ​​to obtain the probability of each label value.

[0157] For ease of understanding, let's take an example and assume that the number of sampled label values ​​with a label value of 1 is 25, the number of sampled label values ​​with a label value of 2 is 15, and the number of sampled label values ​​with a label value of 3 is 10. Then, 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, is calculated respectively 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%.

[0158] Step (3) summarizes the probability of each unknown unit cell at each label value to obtain the edge probability of each unknown unit cell.

[0159] When the probability of label value 1 of the unknown unit cell is 50%, the probability of label value 2 is 30%, and the probability of label value 3 is 20%, the edge probabilities of the unknown unit cell are: label value 1, 50%; label value 2, 30%, and label value 3, 20%.

[0160] Step (4): Based on the edge probability of each unknown unit cell, a predicted label value of each unknown unit cell is obtained.

[0161] In the present invention, the predicted label value of the unknown unit cell has the maximum label probability.

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

[0163] Step S60: updating the initial geological model based on the predicted label value of each unknown unit cell to construct a three-dimensional digital geological model.

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

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

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

[0167]

[0168] Among them, E(x i ) is information entropy, p (m) (x i ) represents the marginal probability of soil type m in unit cell 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 unit cell i, that is, the higher the uncertainty of the unit cell. Therefore, information entropy provides an effective measurement tool for quantifying uncertainty without assuming any specific functional form or statistical model for probability distribution.

[0169] By calculating the entropy value for each unit cell, 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 unit cell uncertainty but also lays the foundation for further geological analysis.

[0170] Example 1

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

[0172] First, a Gibbs sampler is used to simulate and generate a three-dimensional geological model, such as Figure 6 As shown in (a), the modeling size is 30×60×100 (number of grids), and the parameters used in the simulation are set to β = {4.0, 0.1, 0.1, 0.1, 2.0, 0.1, 0.1}. 5% of the columns are randomly selected from the 3D geological case as known drill holes (i.e., 150 drill holes), as shown in Figure 6 (b) is used to estimate the stratigraphic information and stratigraphic uncertainty of the entire area.

[0173] It should be noted that the 30 in the modeling size can be understood as 30 grids in the X direction, 60 grids in the Z direction, and 100 grids in the Y direction. Assuming that each grid is a cubic lattice, the side length of each grid can be 100m, 200m, etc., depending on the actual geological body settings.

[0174] Other schematic diagrams involved in the process of constructing a 3D digital geological model are detailed in Figures 7 to 12 , Figure 7 for Figure 6 (a) A top view of the 3D geological model; Figure 8 For Figure 7 Schematic diagram of the YZ section corresponding to the YZ projection line; Figure 9 For Figure 7 Schematic diagram of the XZ section corresponding to the XZ projection line; Figure 10 for Figure 8 The trajectories of the posterior β1 to β4 on the YZ profile are shown; Figure 11for Figure 9 The trajectories of posterior β2 and posterior β5 to β7 on the XZ section are shown; Figure 12 (a) is the three-dimensional digital geological model map constructed in Example 1, Figure 12 (b) in the figure is the quantification result diagram of stratigraphic uncertainty corresponding to the three-dimensional digital geological model.

[0175] from Figure 10 It can be seen that this round of iterative simulation was performed for 1000 iterations, and the chain of β1 to β4 gradually converged after about the 200th iteration. The posterior estimation results from the 200th iteration to the end of the chain are summarized in Table 1. Figure 11 It can be seen that this round of iterative simulation was performed for a total of 1000 iterations, and the chain of β5 to β7 gradually converged after about the 100th iteration. The posterior estimation results from the 100th iteration to the end of the chain are summarized in Table 2.

[0176] Table 1. Posterior β of YZ plane

[0177]

[0178] Table 2 Posterior β of XZ plane

[0179]

[0180] Example 2

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

[0182] To further validate the performance of the geological model constructed in this paper in actual engineering, we studied deep foundation pit drilling data in Shenzhen. This deep foundation pit drilling data was obtained from a CMC model-based study by Jiang et al. (Jiang QH, Zhang JZ, Zhang DM, et al. Simulation of geological uncertainty based on an improved three-dimensional coupled Markov chain model [J]. Engineering Geology, 2024, 340: 107647).

[0183] A total of 25 borehole data were collected in this project, and seven types of geotechnical soils were found in the boreholes: rockfill, silt, medium sand, silty clay, fully weathered mixed granite (FWMG), strongly weathered mixed granite (SWMG), and moderately weathered mixed granite (MWMG). The plane positions of the 25 boreholes are as follows: Figure 13 As shown in (a), the geological information in the borehole is as follows Figure 13 (b) shown.

[0184] The present invention adopts a cross-validation method to perform geological simulation using different borehole combinations, that is, one of the boreholes is extracted in turn as a validation borehole, and the other 24 boreholes are used as test boreholes to simulate the stratum. The simulation results of the cross-validation method are shown in Table 3. As can be seen from Table 3, the overall performance of the geological model provided by the present invention in this engineering example is very excellent. The accuracy of most boreholes is above 0.90, among which BH7, BH14, BH15, BH19, BH20, BH21, and BH25 have an accuracy of up to 0.96, indicating that the model has high reliability in most cases. At the same time, the accuracy of BH4 (0.75), BH17 (0.72), and BH23 (0.71) is significantly lower than that of other boreholes, which may be due to the complex geological conditions in these areas or the uneven distribution of training data, resulting in a decrease in the generalization ability of the model. Except for individual boreholes, the model performed stably in 24 cross-validations, with an average accuracy of approximately 0.91, verifying the robustness of the model.

[0185] Using BH24 as a verification borehole, the three-dimensional digital geological model is constructed as shown in the figure below. Figure 14 As shown in (a), the quantitative results of the formation uncertainty corresponding to the three-dimensional digital geological model are shown in Figure 14 (b) is shown. Figure 14 It can be seen that the geology of this area is basically distributed in a horizontal state, and the uncertainty of the simulation results is mainly distributed at the junction of strata.

[0186] Table 3 Cross-validation simulation results

[0187]

[0188] In order to compare with the study of Jiang et al., the same borehole combination (i.e., BH14, BH15, BH16, BH17, BH18, and BH19 as verification boreholes and other boreholes as test boreholes) was used for simulation, and the simulation results were compared with those of Jiang et al., as shown in Table 4.

[0189] Table 4 Comparison of research results

[0190]

[0191] As can be seen from Table 4, the accuracy of BH14 to BH19 of the present invention is generally higher than the results of the prior art Jiang et al., which are 0.96, 0.96, 0.95, 0.71, 0.90, and 0.96, respectively. Among them, the accuracy of BH17 is significantly lower, which is consistent with the cross-validation results and may be related to the geological anomaly in the area where the borehole is located.

[0192] See also Figure 15 According to a second aspect of the present invention, the present invention further provides a three-dimensional digital geological model construction device 100, the construction device 100 comprising:

[0193] The first construction module 101 is used to construct a three-dimensional geological model region based on the geological body of the modeling area, and discretize the geological model region into n unit cells and number each unit cell. The modeling region includes a plurality of boreholes opened in the geological body.

[0194] The second construction module 102 is used to map the drilling position and drilling information of each borehole in the modeling area to the n unit cells, complete initial value assignment for the n unit cells, and construct an initial geological model, wherein the label values ​​of multiple known unit cells corresponding to the drilling positions in the initial geological model are determined based on the drilling information, and the initial label values ​​of other unknown unit cells in the initial geological body model that do not correspond to the drilling positions are randomly generated based on the label values ​​of the multiple known unit cells.

[0195] The third construction module 103 is used to construct a neighborhood system for each unit cell based on the Markov random field, and set spatial constraints in seven directions for each unit cell, while describing the strength of the spatial constraints in each different direction with granularity coefficients β1 to β7, wherein the neighborhood system of each unit cell is composed of other unit cells that are adjacent to the unit cell and share nodes with the unit cell.

[0196] The first acquisition module 104 is used to obtain a prediction result based on the initial label value of each unknown unit cell, the label value of each known unit cell, and the neighborhood system of each unit cell, 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 unit cell, where M is a positive integer, wherein: the step of obtaining the prediction result includes: first, based on a two-dimensional geological profile, obtaining the mean and standard deviation of each granularity coefficient through multiple simulation iterations, and then, based on the mean and standard deviation of each granularity coefficient and the initial geological model, obtaining the prediction result through multiple iterations, wherein the two-dimensional geological profile includes a calculated profile of the YZ plane and a calculated profile of the XZ plane.

[0197] The second acquisition module 105 is configured to acquire an edge probability of each unknown unit cell based on the M sampled label values ​​of each unknown unit cell, and acquire a predicted label value of each unknown unit cell based on the edge probability of each unknown unit cell.

[0198] The updating module 106 is configured to update the initial geological model based on the predicted label value of each unknown unit cell to construct a three-dimensional digital geological model.

[0199] 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 unit lattice based on the calculation formula of the edge probability and information entropy of each unknown unit lattice, and obtain the formation uncertainty quantification result corresponding to the three-dimensional digital geological model based on the information entropy of each unknown unit lattice.

[0200] 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.

[0201] 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 digital geological model, 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 unit cells, and numbering each unit cell, 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 unit cells, completing initial value assignment for the n unit cells, and constructing an initial geological model, wherein label values ​​of a plurality of known unit cells corresponding to the drilling positions in the initial geological model are determined based on the drilling information, and initial label values ​​of unknown unit cells in the initial geological body model that do not correspond to the drilling positions are randomly generated based on the total number of label values ​​determined based on all the drilling information; Step S30: constructing a neighborhood system for each unit cell based on the Markov random field, and setting spatial constraints in seven directions for each unit cell. At the same time, the strength of the spatial constraints in each different direction is described by granularity coefficients β1 to β7. The neighborhood system of each unit cell consists of other unit cells that are adjacent to the unit cell and share nodes with it. Step S40: Based on the initial label value of each unknown unit cell, the label value of each known unit cell, and the neighborhood system of each unit cell, 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 unit cell, where M is a positive integer, wherein: The step of obtaining the prediction result includes: first, based on a two-dimensional geological profile, obtaining a mean value and a standard deviation of each particle size coefficient through multiple simulation iterations, and then, based on the mean value and the standard deviation of each particle size coefficient and an initial geological model, obtaining the prediction result through multiple iterations, wherein the two-dimensional geological profile includes a calculated profile of a YZ plane and a calculated profile of an XZ plane; Step S50: obtaining an edge probability of each unknown unit cell based on the M sampled label values ​​of each unknown unit cell, and obtaining a predicted label value of each unknown unit cell based on the edge probability of each unknown unit cell; Step S60: updating the initial geological model based on the predicted label value of each unknown unit cell to construct a three-dimensional digital geological model.

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

3. The method for constructing a three-dimensional digital geological model according to claim 1, wherein: The label probability prediction formula includes a local conditional probability distribution calculation formula, a local energy calculation formula, and a potential function calculation formula, as shown in formulas (1) to (3), where: In formula (1) to formula (3), is the local conditional probability distribution at lattice i; represents local energy; x i represents a specific label at lattice i, x i ' represents any possible label at lattice i, looping on L, where L is the set of label values; V i,j (x i ,x j ) is the potential function, represents the neighborhood system of unit cell i; β d ∈β={β1,β2,β3,β4,β5,β6,β7},x j represents a specific label of other lattices j in the neighborhood system of lattice 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: log(target)=log(Prior(β))+log(L(x BH |b,x unknown )) (5); In formula (4) and formula (5), β represents the particle size coefficient, Prior(β) represents the prior distribution of β, and x unknown Represents the sampled label value of the unknown unit cell, x BH Represents the label value of a known unit cell, represents the neighborhood system of unit cell i.

4. The method for constructing a three-dimensional digital geological model according to claim 3, characterized in that: The step of obtaining the prediction result through multiple iterations based on the mean and standard deviation of each particle size coefficient and the initial geological model includes: Step S401: Using the mean and standard deviation of each particle size coefficient as parameters of a priori multivariate Gaussian distribution, a third priori multivariate Gaussian distribution corresponding to β1 to β7 is established; Step S402: Based on the label value of the known lattice unit, the sampled label value of the unknown unit lattice obtained in the previous iteration, and the third priori multivariate Gaussian distribution, according to the update iteration formula of the granularity coefficient, obtain the sampled values ​​of β1 to β7 obtained in the current iteration, and the third objective function value corresponding to β1 to β7 obtained in the current iteration is greater than the third objective function value corresponding to β1 to β7 obtained in the previous iteration; Step S403: Calculate the local conditional probability distribution of each unknown lattice unit obtained in the current iteration based on the sampled values ​​of β1 to β7 obtained in the current iteration and the label probability prediction formula, and obtain the sampled label value of each unknown unit cell obtained in the current iteration based on the local conditional probability distribution of each unknown unit cell; Step S404, repeat step S402 and step S403 multiple times until β1~β7 converge, and after β1~β7 converge, continue to repeat step S402 and step S403 a third preset number M times, and use the M sampling label values ​​of each unknown unit cell obtained by iterating M times after β1~β7 converge as the prediction result.

5. The method for constructing a three-dimensional digital geological model according to claim 4, characterized in that: The step S402 includes: Step (1), sampling from the third prior distribution using the Metropolis-Hastings method to obtain current sampled values ​​of β1 to β7; Step (2): based on the label value of the known lattice unit, the sampled label value of the unknown lattice unit obtained by the last iterative update, the third prior multivariate Gaussian distribution, and the current sampled values ​​of β1 to β7, according to the likelihood function calculation formula and the objective function calculation formula, obtain the current third objective function value corresponding to the current sampled values ​​of β1 to β7; Step (3): compare the current third objective function value with the third objective function values ​​corresponding to β1 to β7 obtained in the previous iteration. When the current third objective function value is greater than the third objective function values ​​corresponding to β1 to β7 obtained in the previous iteration, use the current sampling values ​​of β1 to β7 as the sampling values ​​of β1 to β7 obtained in the current iteration; when the current third objective function value is less than or equal to the third objective function values ​​corresponding to β1 to β7 obtained in the previous iteration, return to step (1) and resample.

6. The method for constructing a three-dimensional digital geological model according to any one of claims 3 to 5, characterized in that: The first of the seven directions is parallel to the X-axis, and its particle size coefficient is represented by β1; the second direction is parallel to the Z-axis, and its particle size coefficient is represented by β2; the third direction is on the same straight line as the main diagonal of a plane parallel to both the X-axis and the Z-axis, and its particle size coefficient is represented by β3; the fourth direction is on the same straight line as the secondary diagonal of a plane parallel to both the X-axis and the Z-axis, and its particle size coefficient is represented by β4; the fifth direction is parallel to the Y-axis, and its particle size coefficient is represented by β5; the sixth direction is on the same straight line as the secondary diagonal of a plane parallel to both the Y-axis and the Z-axis, and its particle size coefficient is represented by β6; the seventh direction is on the same straight line as the main diagonal of a plane parallel to both the Y-axis and the Z-axis, and its particle size coefficient is represented by β7.

7. The method for constructing a three-dimensional digital geological model according to claim 6, characterized in that: The step of obtaining the mean and standard deviation of each particle size coefficient through multiple simulation iterations based on the two-dimensional geological profile includes: Step S41: Based on the drilling positions of multiple boreholes in the modeling area, an estimated section of the YZ plane is selected, and drilling information of multiple first boreholes on the estimated section of the YZ plane and multiple second boreholes adjacent to the estimated section of the YZ plane are mapped to the estimated section of the YZ plane to obtain a calculated section of the YZ plane; Step S42: Based on the calculated profile of the YZ plane, first converge β1 to β4 through multiple iterations, and after β1 to β4 converge, continue to iterate for a first preset number of times K times to obtain K sample values ​​of β1 to β4 obtained by iterating K times after β1 to β4 converge; Step S43: Calculate the mean and standard deviation of β1 to β4 based on the K sample values ​​of β1 to β4; Step S44: Based on the drilling positions of the multiple boreholes in the modeling area, an estimated cross-section of the XZ plane is selected, and drilling information of the multiple third boreholes on the estimated cross-section of the XZ plane and the multiple fourth boreholes adjacent to the estimated cross-section of the XZ plane are mapped to the estimated cross-section of the XZ plane to obtain a calculated cross-section of the XZ plane; Step S45: Based on the calculated profile of the XZ plane and the mean and standard deviation of β2, first converge β5-β7 through multiple iterations, and after β5-β7 converge, continue to iterate a second preset number L times to obtain L sample values ​​of β5-β7 obtained by iterating L times after β5-β7 converge; Step S46: Calculate the mean and standard deviation of β5 to β7 based on the L sample values ​​of β5 to β7.

8. The method for constructing a three-dimensional digital geological model according to claim 7, characterized in that: Step S42 includes: Obtaining a number of each first unit cell and a number of each second unit cell on the calculated section of the YZ plane, wherein the first unit cell is a unit cell of an unknown rock and soil type, and the second unit cell is a unit cell of a known rock and soil type; Based on the preset means and standard deviations of β1 to β4, a first priori multivariate Gaussian distribution corresponding to β1 to β4 is established; Based on the initial label value of each first unit cell, the label value of each second unit cell, and the neighborhood system of each unit cell, an iterative operation is performed according to the label probability prediction formula and the update iteration formula of the granularity coefficient, first through multiple iterations until β1 to β4 converge, and after β1 to β4 converge, continue to iterate a first preset number of times K times, and obtain K sample values ​​of β1 to β4 obtained by iterating K times after β1 to β4 converge, wherein any iterative update of β1 to β4 includes: Step (a), using the Metropolis-Hastings method to sample from the first prior multivariate Gaussian distribution to obtain current sampled values ​​of β1 to β4; Step (d), based on the label values ​​of the plurality of second unit cells, the sampled label values ​​of the first unit cell obtained by the last iterative update, the first prior multivariate Gaussian distribution, and the current sampled values ​​of β1 to β4, obtaining a current first objective function value corresponding to the current sampled values ​​of β1 to β4 according to a likelihood function calculation formula and an objective function calculation formula; Step (c), comparing the current first objective function value with the first objective function values ​​corresponding to β1 to β4 obtained in the previous iteration; when the current first objective function value is greater than the first objective function values ​​corresponding to β1 to β4 obtained in the previous iteration, using the current sampled values ​​of β1 to β4 as the sampled values ​​of β1 to β4 obtained in the current iteration; when the current first objective function value is less than or equal to the first objective function values ​​corresponding to β1 to β4 obtained in the previous iteration, returning to step (a) for resampling; Step (d): Based on the sampling values ​​of β1 to β4 obtained in the current iteration and the label probability prediction formula, obtain the local conditional probability distribution of the first unit lattice obtained by the current iteration update, and based on the local conditional probability distribution of each first unit lattice, obtain the sampling label value of each first unit lattice obtained in the current iteration.

9. The method for constructing a three-dimensional digital geological model according to claim 8, characterized in that: Step S45 includes: Obtaining the number of each third unit cell and the number of each fourth unit cell on the calculated section of the XZ plane, wherein the third unit cell is a unit cell of an unknown rock and soil type, and the fourth unit cell is a unit cell of a known rock and soil type; Based on the preset mean and standard deviation parameters of β5-β7 and the mean and standard deviation of β2, a second prior multivariate Gaussian distribution corresponding to β2 and β5-β7 is established; Based on the second priori multivariate Gaussian distribution, the initial label value of each third unit cell, the label value of each fourth unit cell, and the neighborhood system of each unit cell, an iterative operation is performed according to the label probability prediction formula and the update iteration formula of the granularity coefficient, first through multiple iterations until β5 to β7 converge, and after β5 to β7 converge, continue to iterate a second preset number of times L times, and obtain K sample values ​​of β5 to β7 obtained by iterating L times after β5 to β7 converge, wherein any iterative update of β5 to β7 includes: Step (a), sampling from the second prior multivariate Gaussian distribution using the Metropolis-Hastings method to obtain current sampled values ​​of β2 and β5-β7; Step (d), based on the label values ​​of the plurality of fourth unit cells, the sampled label values ​​of the third unit cell obtained by the last iterative update, the second prior multivariate distribution, and the current sampled values ​​of β2 and β5-β7, obtaining a current second objective function value corresponding to the current sampled values ​​of β5-β7 according to the likelihood function calculation formula and the objective function calculation formula; Step (c), comparing the current second objective function value with the second objective function values ​​corresponding to β5 to β7 obtained in the previous iteration; when the current second objective function value is greater than the second objective function values ​​corresponding to β5 to β7 obtained in the previous iteration, using the current sampled values ​​of β5 to β7 as the sampled values ​​of β5 to β7 obtained in the current iteration; when the current second objective function value is less than or equal to the second objective function values ​​corresponding to β5 to β7 obtained in the previous iteration, returning to step (a) for resampling; Step (d): Based on the sampling values ​​of β2 and β5~β7 obtained in the current iteration and the label probability prediction formula, obtain the local conditional probability distribution of the third unit lattice obtained by the current iteration update, and based on the local conditional probability distribution of each third unit lattice, obtain the sampling label value of each third unit lattice obtained in the current iteration.

10. A device for constructing a three-dimensional digital geological model, 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 unit cells, and number each unit cell, 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 unit cells, perform initial value assignment on the n unit cells, and construct an initial geological model, wherein label values ​​of a plurality of known unit cells corresponding to the drilling position in the initial geological model are determined based on the drilling information, and initial label values ​​of other unknown unit cells in the initial geological body model that do not correspond to the drilling position are randomly generated based on the label values ​​of the plurality of known unit cells; The third building block is used to construct a neighborhood system for each unit cell based on a Markov random field, and set spatial constraints in seven directions for each unit cell. The granularity coefficients β1 to β7 are used to describe the strength of the spatial constraints in each different direction. The neighborhood system of each unit cell consists of other unit cells that are adjacent to the unit cell 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 unit cell, the label value of each known unit cell, and the neighborhood system of each unit cell, 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 unit cell, where M is a positive integer, wherein: The step of obtaining the prediction result includes: first, based on a two-dimensional geological profile, obtaining a mean value and a standard deviation of each particle size coefficient through multiple simulation iterations, and then, based on the mean value and the standard deviation of each particle size coefficient and an initial geological model, obtaining the prediction result through multiple iterations, wherein the two-dimensional geological profile includes a calculated profile of a YZ plane and a calculated profile of an XZ plane; A second acquisition module is used to acquire an edge probability of each unknown unit cell based on the M sampled label values ​​of each unknown unit cell, and acquire a predicted label value of each unknown unit cell based on the edge probability of each unknown unit cell; The updating module is used to update the initial geological model based on the predicted label value of each unknown unit cell to construct a three-dimensional digital geological model.

Images