CN108416841B - Automatic reconstruction method of multi-point statistical three-dimensional geological model based on local search strategy - Google Patents
Automatic reconstruction method of multi-point statistical three-dimensional geological model based on local search strategy Download PDFInfo
- Publication number
- CN108416841B CN108416841B CN201810060685.7A CN201810060685A CN108416841B CN 108416841 B CN108416841 B CN 108416841B CN 201810060685 A CN201810060685 A CN 201810060685A CN 108416841 B CN108416841 B CN 108416841B
- Authority
- CN
- China
- Prior art keywords
- sub
- probability
- grid
- simulation
- current
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- Computer Graphics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Remote Sensing (AREA)
- Geophysics And Detection Of Objects (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a multipoint statistic three-dimensional geological model automatic reconstruction method based on a local search strategy, which directly utilizes a two-dimensional horizontal and vertical geological profile to reconstruct a three-dimensional geological model, divides a three-dimensional modeling space into a plurality of sub-regions according to the spatial cross relationship of a known profile, and each sub-region is surrounded by N (N is more than or equal to 1 and less than or equal to 6) sub-sections; for unknown attribute positions distributed in each sub-region, acquiring N surrounding sub-sections according to the positions of the unknown attribute positions, customizing X, Y, Z data events in three directions according to set neighborhood parameters and known data points, and scanning the sub-sections in the corresponding directions by using the corresponding data events to acquire N conditional probability density functions; and then obtaining the joint probability density distribution by using a reasonable probability fusion strategy, and finally obtaining the attribute value of the position. The invention can be popularized and applied in various three-dimensional geological information systems, geographic information systems and geological modeling and simulation systems.
Description
Technical Field
The invention relates to the field of three-dimensional geological modeling and geological statistical simulation, in particular to a multipoint statistical three-dimensional geological model automatic reconstruction method based on a local search strategy.
Background
The three-dimensional geological modeling and visualization technology is one of core technologies in the field of geospatial information science, and is also a key and basic technology for the construction of three-dimensional geographic and geological information systems. The three-dimensional visual representation of the geological structure can provide more real and visual description of the underground geological phenomena and structures, and plays an important role in quantitative evaluation of underground water, petroleum and solid mineral resources. The spatial heterogeneity and connectivity of reservoirs play an important control role in the migration of subsurface fluids, which is crucial to quantifying and predicting the formation and distribution of subsurface resources. However, for a practical three-dimensional application, these properties are difficult to describe and model because the known data we can obtain is very sparse. Traditional two-point geostatistical methods have been used to reconstruct three-dimensional models of subsurface sedimentary structures and phase distributions. However, two-point geostatistical methods have difficulty capturing higher order statistics, and thus the anisotropic features and connectivity patterns are not well reproduced. Object-based methods are also inadequate for three-dimensional reconstruction of heterogeneous lithofacies due to the inability to parameterize complex geological structures and phenomena. The multi-Point Statistics (MPS) method is an important branch of the field of automatic construction and simulation of three-dimensional geological models, the development of the basic theory and the proposal of various algorithms provide technical support for three-dimensional reconstruction of anisotropic underground complex geological structures, and the MPS method is well applied to a plurality of geological fields such as reservoir simulation, hydrogeological modeling, porous medium reconstruction and the like.
The multipoint statistical method aims at describing spatially heterogeneous geometric features by extracting spatial patterns directly from training images rather than by using a variogram. The training image is a conceptual model abstracted from actual phenomena and plays a crucial role in stochastic simulation based on multipoint statistics. Theory of multi-point geostatistics since its introduction, various stochastic simulation and modeling methods based on this theory have been developed. In these methods, a training image is scanned using a fixed search template and the multipoint statistical patterns are stored in a tree or linked list data structure. Then, a conditional probability density function (cpdf) is calculated according to the current data event, and a simulation result of the current position is sequentially obtained. However, one notable contradiction is in the size of the data template and the reproducibility of the large-scale structures. The larger the data template, the larger the size of the data event, and therefore the fewer the number of repetitions of capturing the data event on the training image, which will make it difficult to infer the corresponding conditional probability density function. Conversely, when the data template is too small, the large scale structure of the training image will not be reproduced. In addition, searching for templates that are too large increases computational cost and memory consumption. The multiple-grid concept alleviates the above limitations, but the above problem is not completely eliminated due to the data template being fixed. As an approximate expression, the pattern distance may be used to represent the degree of match between the pattern in the training image and the data event. For the multipoint statistical method based on the probability density function, the use of the pattern distance can increase the matching success probability and reduce the number of pattern storage to some extent, but some complex detailed features can be filtered out due to the approximate expression.
A suitable, reliable training image is the basis of a stochastic simulation method based on multipoint statistics. For practical three-dimensional applications, the acquisition of three-dimensional training images is difficult, and the reliability thereof is difficult to evaluate. To overcome this difficulty, some methods attempt to reconstruct three-dimensional models of subsurface complex geological structures directly using low-dimensional data (e.g., boreholes, sections, outcrops, remote sensing images, etc.) rather than using complete three-dimensional training images. Such methods include two main types. One widely used method is Direct Sampling (DS) which directly takes an incomplete known data set as reference, obtains a multipoint statistical mode from the reference, and sequentially completes the reconstruction of the whole three-dimensional modeling space. The main advantage of this approach is that spatial heterogeneity in various directions is taken into account in the pattern matching between the data events and the known data. However, a fine, large-scale three-dimensional model will contain tens of millions of simulation meshes. To find a matching multi-point pattern from sparse known data, each grid would require millions of attempts to obtain the final result. Thus, for fine three-dimensional applications, this approach to search for matches on global data still suffers from computational efficiency. Another approach uses two-dimensional geological profiles as training images to fill the entire three-dimensional modeling space by performing a series of two-dimensional simulations. Although this method is highly efficient, spatial continuity of the simulation results is difficult to ensure.
From another perspective, it is a common workflow of geological workers to obtain two-dimensional geological data (e.g., geological plans, profiles, etc.) by means of sparse, scattered one-dimensional borehole, etc. sampling data, appropriate simulation methods, and expert knowledge and experience. It would then be logical and feasible to reconstruct a three-dimensional geological structure model using reliable simulation methods based on low-dimensional data. However, in practical modeling practice, due to geological phenomena and spatial heterogeneity of the process, the direct use of real geological profiles, remote sensing images and the like as training images for three-dimensional reconstruction of complex geological structures is still limited by non-stationarity characteristics of the real geological profiles, the remote sensing images and the like.
Probabilistic fusion is the basis for multivariate statistical analysis. A geologic phenomenon is often affected by factors from different aspects or angles, which requires a final joint probability distribution using a probabilistic fusion approach. The probability fusion method can be classified into an addition-based method and a multiplication-based method according to its mathematical properties. The Linear Pooling formula is a widely used method based on probability addition, and is mainly characterized by flexibility and simplicity. Multiplication methods include the Bordley formula and the Log-Linear Pooling formula. It has been shown that for the case of only two attributes, the two equations are equivalent, but the latter has better performance for non-binary events. The linear probability fusion method is used to combine several independent probability distributions to obtain a larger and more stable combined probability distribution, just like the logical operation "OR". While multiplication-based methods are more like the logical operation "AND", such methods are typically used to fuse probability distributions with significant correlation.
In summary, due to the limitation of the lack of a complete and reliable three-dimensional training image, the traditional random simulation method based on multipoint statistics is difficult to be well applied to the simulation and reconstruction of the actual three-dimensional geological structure. Therefore, reconstructing a three-dimensional model of the subsurface geological structure using existing low-dimensional data is a viable and amenable path for geological data processing. The existing multipoint statistical reconstruction method for performing global search matching on all known data is limited by the computational efficiency. The method of gradually filling the three-dimensional space domain by simulating the two-dimensional slice by using the multipoint statistic technology is difficult to ensure the spatial continuity of the reconstruction results in different directions of the whole space. In view of the above problems and challenges of the analytical multi-point statistical method in the three-dimensional reconstruction of complex geological structures, it is necessary to provide an efficient automatic reconstruction method of a multi-point statistical three-dimensional geological structure model based on low-dimensional data.
Disclosure of Invention
In view of this, the embodiment of the present invention provides an efficient method for automatically reconstructing a multi-point statistical three-dimensional geological model based on a local search strategy, which can characterize spatial heterogeneity, and specifically includes the following steps:
(1) loading data and distributing the known data into a simulation grid;
(2) according to the known position of each two-dimensional horizontal and vertical geological section in the simulation grid, the index of each two-dimensional horizontal and vertical geological section in the direction of X, Y, Z is obtained, and the prior probability P of all the sub-regions divided by the horizontal and vertical sections is calculated0；
(3) Determining multiple grids G, and processing each multiple grid in sequence until all multiple grids are processed; acquiring a current grid g, if the grid is effective, turning to the step (4), otherwise, turning to the step (13);
(4) defining a random simulation path for the current grid g so as to cover all unknown nodes on the current grid g;
(5) obtaining the current simulation position x from the random simulation path in sequence, if the current node exists, turning to the step (6), otherwise, turning to the step (12);
(6) determining X, Y, Z search neighborhoods N corresponding to three directions according to the position of the current node in the simulation gridX、NY、NZAnd corresponding data event devX、devY、devZ；
(7) Obtaining the indices { x ] of all N sub-sections around the sub-region where the current simulation position x is located0,x1}、{y0,y1}、{z0,z1N is more than or equal to 1 and less than or equal to 6;
(8) taking the sub-sections obtained in the step (7) as training images, and sequentially using data events devX、devY、devZScanning the sub-sections in the corresponding directions to obtain corresponding N conditional probability density functions { cpdf1,...,cpdfN}；
(9) From all the sub-sections, N conditional probability density functions { cpdf } corresponding to the current simulation position x are obtained1,...,cpdfNH, then the prior probability P obtained from step (2)0The prior probability p of the sub-region corresponding to the current simulation position x is obtained0；
(10) Fusing the obtained N conditional probability density functions { cpdf according to a probability fusion strategy1,...,cpdfNAnd a priori probability p0Thereby obtaining a joint probability distribution;
(11) randomly extracting a corresponding attribute from the final joint probability distribution, assigning the attribute to the current simulation position x, and turning to the step (5) to continue the attribute simulation of other unknown positions;
(12) all unknown nodes on the current grid g are simulated, and the step (3) is switched to simulate the next repeated grid;
(13) all the unknown positions on each multiple grid G are simulated, the simulation result is stored, and the reconstruction is finished;
the attribute values of all nodes on the whole simulation grid are completely obtained, so that a three-dimensional geological structure model which fills the whole modeling space is obtained.
Further, in the step (8), for the scanning process of any one of the sub-sections, the following steps (8-1) to (8-8) are performed:
(8-1) acquiring all nodes of the currently scanned sub-section according to the index of the sub-section in the whole simulation grid;
(8-2) setting a random scanning path and initializing a counter sum of the number of matched patterns to be 0;
(8-3) sequentially obtaining corresponding positions path in the sub-section from the random scanning pathiiIf ii is<If the path () × f is true, turning to the step (8-4), otherwise, turning to the step (9); where f is the proportion of a given scanned sub-section, i.e. allowing a portion of the sub-section to be scanned rather than the entirety;
(8-4) according to pathiiAcquiring a corresponding neighborhood N 'and a data event dev' at the position in the training image;
(8-5) calculating a pattern distance d { dev, dev '} between the data event dev from the simulation grid and the data event dev' from the training image using the following formula:
where k represents the number of nodes contained in the data event, Z (x)j) And Z (y)j) Respectively representing the attribute values of the central nodes corresponding to the data events in the simulation grid and the training image;
(8-6) if d { dev, dev' } t is true, turning to the step (8-7), otherwise, turning to the step (8-8); wherein t is a threshold value of the pattern distance corresponding to the current grid g, and if the current pattern distance is smaller than the threshold value, the current pattern distance is received;
(8-7) updating the conditional probability density function cpdf corresponding to the current sub-section according to the attribute value corresponding to the data event center node in the training image, and enabling sum + +;
(8-8) if sum > maximumN is true, turning to the step (9), otherwise, turning to the step (8-3) to continue scanning the current sub-section; where maximumN represents the maximum number of multi-point patterns to be acquired from the current training image.
Further, it is characterized byIn the step (8), when scanning the sub-section, the search neighborhood radius R of the multiple grids and the threshold t of the mode distance are allowed to be linearly decreased with the increase of the grid level; the neighborhood radius R and the threshold t on the first rescreen are given by the user, i.e. R ═ R0，t＝t0(ii) a Fixing the grid to be R1 and t 0 on the last grid; for any intermediate mesh, its neighborhood radius R and threshold t can be expressed as:
where m represents the total number of multiple grids and τ (1 ≦ τ ≦ m) represents the index of the current grid.
Further, in the step (8-8), a recommended value range of maximumN is [50,100], where maximumN represents the maximum number of multi-point patterns obtained from the current sub-section, and is also properly adjusted according to the scale of the model mesh and the complexity of the geological structure; therefore, the stability of the conditional probability density function cpdf is controlled by controlling the number of successfully matched modes, and unnecessary calculation time is avoided.
Further, in the step (10), the probability fusion strategy comprises the following two steps:
(10-1) fusing conditional probability density functions from two sub-sections parallel to each other in the same direction using Linear Pooling formula based on probability addition:
wherein A represents a probabilistic fused attribute class, PG(A) Representing the probability, P, of the property after fusioni(A) Representing probability of property, w, from any sub-sectioniIs the corresponding weight; at most two sections are arranged in the parallel direction of one subregion, and probability fusion is needed only when two sections exist; when m is 2, weight w1、w2The value of (d) is calculated by the following formula:
wherein d is1、d2Respectively representing the distances from the current simulation position to the two sub-sections;
(10-2) fusing probability distribution in three orthogonal directions processed in the step (10-1) and prior probability by using a Log-Linear Pooling formula based on probability multiplication:
wherein A represents the attribute category of probability fusion, n is the number of probability distribution participating in probability fusion, n is more than or equal to 1 and less than or equal to 3, namely the probability distribution in X, Y, Z three orthogonal directions obtained after the processing of step (10-1) at most, PG(A) Representing the probability, P, of the property after fusioni(A) Representing the probability of the attribute, P, of any sub-direction after the processing in step (10-1)0(A) Is the prior probability, w, of the propertyiWeights for the probabilities in each direction; weights w comprising a priori probabilities0Within, the sum must be 1, i.e.
Through the adjustment of each weight, the influence degree of the probability distribution in different directions on the final combined probability distribution is different, and the larger the weight in a certain direction is, the closer the final result is to the probability distribution condition; weight of prior probability
Namely when
When w0Is positive, the prior probability also participates in the probability fusion, and the weight w0The larger the effect, the correspondingly larger; when in use
When w0When the prior probability is 0, the regulation function fails; when in use
When w0The prior probability is a negative number, and the prior probability is a negative regulation effect, namely the final probability distribution deviates from the prior probability distribution; in the weight selection process, it is necessary to ensure that the sum of all weights is 1.
Compared with the prior art, the invention has the following beneficial effects:
(1) the invention uses low-dimensional data (such as two-dimensional section, one-dimensional drilling and other sampling data) as the basis, and uses the transverse and longitudinal sections to replace complete three-dimensional training images, thereby realizing the automatic reconstruction of a complex three-dimensional geological structure model, providing a new thought for the reconstruction from the low-dimensional geological data to the three-dimensional geological model, and simultaneously providing a new technology and a new method for the three-dimensional geological simulation and modeling based on multi-point statistics.
(2) According to the local search strategy provided by the invention, the acquisition of the multipoint statistical information is fixed to the local sub-section closest to the current simulation position, so that a multipoint mode too far away is prevented from being captured, and the reliability of the simulation result is ensured; and the search strategy of the local sub-section reduces the search times of the training image, thereby improving the simulation efficiency.
(3) The multiple grids and the parameter adjustment strategy of the search neighborhood adopted in the invention allow the thresholds of the radius and the mode distance of the search neighborhood to decrease with the increase of the grid level; on the initial multiple grids, the neighborhood radius is large, and the threshold value of the mode distance is also large, so that the acquisition of a large-scale structural mode is ensured; as the neighborhood radius and threshold decrease, the pattern matching will tend to be a perfect match, which will ensure the reconstruction of fine structural features; the strategy alleviates the problems of fixed search template and inapplicable size in the traditional random simulation method based on multipoint statistics to a certain extent.
(4) In the probability fusion strategy and the weight selection method, the final joint probability distribution comprises probability distributions of sub-sections from all directions; the probability fusion of the parallel sub-sections based on addition expands the capacity of probability distribution, so that the combined probability after fusion is more stable; probability fusion based on multiplication in the orthogonal direction ensures the final depiction of the multipoint modes in different directions, so that the simulation result shows good spatial continuity in different directions.
(5) The limitation on the maximum matching mode number in the invention not only ensures the acquisition of a stable conditional probability density function, but also avoids invalid multipoint statistical search and mode matching, thereby improving the efficiency of three-dimensional model reconstruction to a certain extent.
(6) The method can be popularized and applied to various three-dimensional geological information systems, geographic information systems, geological modeling and simulation systems and other software.
Drawings
FIG. 1 is a general flow of the method for automatically reconstructing a multi-point statistical three-dimensional geological model based on a local search strategy.
FIG. 2 is a diagram illustrating a local search strategy in an embodiment of the invention.
FIG. 3 is a diagram illustrating the search neighborhood, data events, and pattern matching process on adjacent subsections in an embodiment of the invention.
FIG. 4 is a diagram illustrating the relationship between multiple grids and search neighborhoods in an embodiment of the present invention.
FIG. 5 is a three-dimensional reference model, a staggered profile, and two reconstruction results of two different methods employed in embodiments of the invention: DS process and process of the invention.
FIG. 6 is a comparison of statistical properties of modeled results in an embodiment of the present invention, where each method includes 20 reconstructed results: FIG. 6(a) is a spatial connectivity curve of the reconstructed results of the two methods in X, Y, Z three directions; FIG. 6(b) is a plot of the overall variation function of the reconstruction results for both methods; fig. 6(c) is a comparison of the computational efficiency of the two methods as the grid scale increases.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, embodiments of the present invention will be further described with reference to the accompanying drawings.
The embodiment of the invention provides a multipoint statistic three-dimensional geological model automatic reconstruction method based on a local search strategy, which directly utilizes a two-dimensional horizontal and vertical geological profile to reconstruct a three-dimensional geological structure model and provides a local search strategy based on multipoint statistics; according to the strategy, a three-dimensional modeling space is divided into a plurality of sub-regions according to the spatial staggered relation of known profiles, and each sub-region is surrounded by N (N is more than or equal to 1 and less than or equal to 6) sub-sections; for the unknown attribute positions distributed in each sub-region, the multipoint statistical information is directly obtained from N sub-sections surrounding the sub-region: firstly, acquiring N surrounding sub-sections according to the positions of the sub-sections, customizing X, Y, Z data events in three directions according to set neighborhood parameters and known data points, and scanning the sub-sections in the corresponding directions by using the corresponding data events to obtain N conditional probability density functions; and then obtaining the joint probability density distribution by using a reasonable probability fusion strategy, and finally obtaining the attribute value of the position.
Referring to fig. 1, the method includes the following steps:
(1) loading data, and distributing known data points (discrete horizontal and vertical geological profiles, remote sensing images or drilling sampling data) into a simulation grid;
(2) according to the known positions of the two-dimensional horizontal and vertical geological profiles in the simulation grid, the indexes of the profiles in the direction of X, Y, Z in the simulation grid are obtained, and the prior probability P of all sub-regions divided by the horizontal and vertical profiles is calculated0；
Prior probability P0Is the probability distribution of the attributes of the known data points, since the local search strategy in the present invention only seeks to search several sub-sections around the current position to be simulated, the prior probability of each sub-region will also be obtained from these sub-sections around the sub-region.
(3) Determining multiple grids G, and processing each multiple grid in sequence until all multiple grids are processed; acquiring a current grid g, if the grid is effective, turning to the step (4), otherwise, turning to the step (13);
(4) defining a random simulation path for the current grid g so as to cover all unknown nodes on the current grid g;
(5) obtaining the current simulation position x from the random simulation path in sequence, if the current node exists, turning to the step (6), otherwise, turning to the step (12);
(6) determining X, Y, Z search neighborhoods N corresponding to three directions according to the position of the current node in the simulation gridX、NY、NZAnd corresponding data event devX、devY、devZ；
(7) Obtaining the index { x ] of all N (1 ≦ N ≦ 6) subsections surrounding the subregion where the current simulation position x is located0,x1}、{y0,y1}、{z0,z1}；
As shown in fig. 2, each sub-area is generally surrounded by 6 sub-sections, but sometimes the sub-sections are less than 6; while also allowing no profile in one or both directions, but at least in one direction;
(8) using the sub-sections as training images, and using data events dev in sequenceX、devY、devZScanning the sub-sections in the corresponding directions to obtain corresponding N conditional probability density functions { cpdf1,...,cpdfN}；
FIG. 3 is a schematic diagram of the neighborhood search and data events on a sub-section, and for the central node, the search window may be spread over the entire sub-section, while other nodes in the data event allow it to be grabbed from the neighboring sub-sections, which will ensure the spatial continuity of the simulation results around the section segmentation; fig. 4 is a diagram illustrating the relationship between a 3-fold grid and a search neighborhood.
For the scanning process of any one of the sub-sections, the following steps (8-1) to (8-8) are performed:
(8-1) acquiring all nodes of the currently scanned sub-section according to the index of the sub-section in the simulation grid;
(8-2) setting a random scanning path and initializing a counter sum of the number of matched patterns to be 0;
(8-3) sequentially obtaining corresponding positions path in the sub-section from the scanning pathiiIf ii is<If the path () × f is true, turning to the step (8-4), otherwise, turning to the step (9); where f is the proportion of a given scanned cross-section, i.e. allowing a portion but not all of the sub-sections to be scanned;
(8-4) according to pathiiAcquiring a corresponding neighborhood N 'and a data event dev' at the position in the training image;
(8-5) calculating a pattern distance d { dev, dev '} between the data event dev from the simulation grid and the data event dev' from the training image using the following formula:
where k represents the number of nodes contained in the data event, Z (x)j) And Z (y)j) Respectively representing the attribute values of the central nodes corresponding to the data events in the simulation grid and the training image;
(8-6) if d { dev, dev' } t is satisfied (t is a threshold value of the pattern distance corresponding to the current grid g), turning to the step (8-7), otherwise, turning to the step (8-8);
the search neighborhood radius R of the multiple grids and the threshold t of the mode distance are allowed to decrease linearly as the grid level increases; the neighborhood radius R and the threshold t on the first rescreen are given by the user, i.e. R ═ R0，t＝t0(ii) a In the example of implementation given in the present invention (shown in FIG. 5), an initial neighborhood radius R is given0Initial threshold t 300The final heavy grid is fixed to be R1 and t 0; in the middle of any mesh, the neighborhood radius R and the threshold t can be expressed as:
where m represents the total number of multiple grids and τ (1 ≦ τ ≦ m) represents the index of the current grid.
(8-7) updating the conditional probability density function cpdf corresponding to the current sub-section according to the attribute value corresponding to the data event center node in the training image, and enabling sum + +;
(8-8) if sum > maximumN is true, turning to the step (9), otherwise, turning to the step (8-3) to continue scanning the current sub-section; where maximumN denotes the maximum number of multi-point patterns obtained from the current sub-section, which is set to 80 in fig. 5; maximumN controls the stability of the conditional probability density function cpdf by controlling the number of successfully matched patterns; if maximumN is too small, the cpdf capacity is too small, and each probability distribution is unstable, so that the uncertainty of a simulation result is increased; with the increase of maximummN, the cpdf tends to be stable, so that the excessively large maximummN brings meaningless search and matching, and the calculation time is increased; one recommended value range is [50,100], and proper adjustment is made according to the scale of the model grid and the complexity of the geological structure;
(9) from P0The prior probability p of the sub-region corresponding to the current simulation position x is obtained0；
(10) Fusing the obtained N conditional probability density functions { cpdf according to a probability fusion strategy1,...,cpdfNAnd a priori probability p0Thereby obtaining a joint probability distribution;
the probability fusion strategy comprises the following two steps:
(10-1) fusing conditional probability density functions from two sub-sections parallel to each other in the same direction using Linear Pooling formula based on probability addition:
wherein A represents a probabilistic fused attribute class, PG(A) Representing the probability, P, of the property after fusioni(A) Representing probability of property, w, from any sub-sectioniIs the corresponding weight; in the invention, at most two sections are arranged in the parallel direction of one subregion, and probability fusion is needed only when two sections existCombining; when m is 2, weight w1、w2The value of (d) is calculated by the following formula:
wherein d is1、d2Respectively representing the distances from the current simulation position to the two sub-sections;
(10-2) fusing probability distribution and prior probability in three orthogonal directions after the first step of processing by using a Log-Linear Pooling formula based on probability multiplication:
where A denotes the attribute class of probability fusion, n is the number of probability distributions participating in probability fusion, PG(A) Representing the probability, P, of the property after fusioni(A) Representing the probability of the property, P, of any sub-direction after the first step of processing0(A) Is the prior probability, w, of the propertyiWeights for the probabilities in each direction; w is aiThe selection of the sub-section is determined by the complexity (importance) of the corresponding sub-section, if an upward profile structure is more complex and can provide more abundant multi-point patterns, a relatively larger corresponding weight is set, but it is necessary to ensure that
In the example of FIG. 5, let w be the same since the pattern complexity is similar in the three directions0＝w1＝w2＝w3＝0.25，w0The weights are corresponding to the prior probabilities.
(11) Randomly extracting a corresponding attribute from the final joint probability distribution, and assigning the attribute to the current simulation position x; continuing attribute simulation of other unknown positions by turning to the step (5);
(12) all unknown nodes on the current grid g are simulated, and the step (3) is switched to simulate the next repeated grid;
(13) and (4) completing the simulation of all the unknown positions on each multiple grid, storing the simulation result and finishing the reconstruction.
At this point, the attribute values of all nodes on the whole simulation grid are completely obtained, so that the three-dimensional geological structure model which fills the whole modeling space is obtained.
Fig. 5 shows an experimental example of the invention implemented to illustrate the superiority of the method, as a three-dimensional reference model, the interleaved profiles used (X, Y, Z two profiles in each direction), and two corresponding reconstructions of two different methods (DS method and method of the invention).
In order to highlight the comparative analysis of the statistical characteristics of the modeled results, 20 reconstructions were performed per method in this experimental case, and the statistical characteristics of the 20 output results are shown in fig. 6. The spatial connectivity is an important standard for measuring the reconstruction quality of the three-dimensional geological structure model, and fig. 6(a) is a spatial connectivity curve of the reconstruction results of the two methods in X, Y, Z three directions, so that it can be seen that the connectivity of the output model of the method in each direction is closer to that of a given three-dimensional reference model. Another measure of the anisotropy of the geological structure is the variability of the property, and fig. 6(b) is an overall variation function curve of the reconstruction results of the two methods, and the results also show that the method has stronger advantages in the reconstruction capability of the heterogeneity of the geological structure. Fig. 6(c) is a comparison curve of the computational efficiencies of the two methods with the increase of the grid scale, and it can be seen that the acceleration ratio of the computational efficiency of the method of the present invention is gradually increased with the increase of the number of simulation grids, so that the method has the capability of implementing the automatic reconstruction of a large-scale fine geological structure model.
The invention takes low-dimensional data (such as sampling data of a two-dimensional section, a one-dimensional drilling hole and the like) as the basis, and replaces a complete three-dimensional training image with a transverse section and a longitudinal section, thereby realizing the automatic reconstruction of a complex three-dimensional geological structure model, providing a new thought for the reconstruction from the low-dimensional geological data to the three-dimensional geological model, and simultaneously providing a new technology and a new method for the three-dimensional geological simulation and modeling based on multi-point statistics; according to the local search strategy provided by the invention, the acquisition of the multipoint statistical information is fixed to the local sub-section closest to the current simulation position, so that a multipoint mode too far away is prevented from being captured, and the reliability of the simulation result is ensured; the search strategy of the local sub-section reduces the search times of the training image, thereby improving the simulation efficiency; the multiple grids and the parameter adjustment strategy of the search neighborhood adopted in the invention allow the thresholds of the radius and the mode distance of the search neighborhood to decrease with the increase of the grid level; on the initial multiple grids, the neighborhood radius is large, and the threshold value of the mode distance is also large, so that the acquisition of a large-scale structural mode is ensured; as the neighborhood radius and threshold decrease, the pattern matching will tend to be a perfect match, which will ensure the reconstruction of fine structural features; the strategy reduces the problems that a search template is fixed and the size is not too large in the traditional random simulation method based on multipoint statistics to a certain extent; in the probability fusion strategy and the weight selection method, the final joint probability distribution comprises probability distributions of sub-sections from all directions; the probability fusion of the parallel sub-sections based on addition expands the capacity of probability distribution, so that the combined probability after fusion is more stable; probability fusion based on multiplication in the orthogonal direction ensures the final depiction of the multipoint modes in different directions, so that the simulation result shows good spatial continuity in different directions; the limitation on the number of the maximum matching modes ensures the acquisition of a stable conditional probability density function, avoids invalid multipoint statistical search and mode matching, and improves the efficiency of three-dimensional model reconstruction to a certain extent; the method can be popularized and applied to various three-dimensional geological information systems, geographic information systems, geological modeling and simulation systems and other software.
The features of the embodiments and embodiments described herein above may be combined with each other without conflict.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents, improvements and the like that fall within the spirit and principle of the present invention are intended to be included therein.
Claims (5)
1. The method for automatically reconstructing the multi-point statistical three-dimensional geological model based on the local search strategy is characterized by comprising the following steps of:
(1) loading data and distributing the known data into a simulation grid;
(2) according to the known position of each two-dimensional horizontal and vertical geological section in the simulation grid, the index of each two-dimensional horizontal and vertical geological section in the direction of X, Y, Z is obtained, and the prior probability P of all the sub-regions divided by the horizontal and vertical sections is calculated0；
(3) Determining multiple grids G, and processing each multiple grid in sequence until all multiple grids are processed; acquiring a current grid g, if the grid is effective, turning to the step (4), otherwise, turning to the step (13);
(4) defining a random simulation path for the current grid g so as to cover all unknown nodes on the current grid g;
(5) obtaining the current simulation position x from the random simulation path in sequence, if the current node exists, turning to the step (6), otherwise, turning to the step (12);
(6) determining X, Y, Z search neighborhoods N corresponding to three directions according to the position of the current node in the simulation gridX、NY、NZAnd corresponding data event devX、devY、devZ；
(7) Obtaining the indices { x ] of all N sub-sections around the sub-region where the current simulation position x is located0,x1}、{y0,y1}、{z0,z1N is more than or equal to 1 and less than or equal to 6;
(8) taking the sub-sections obtained in the step (7) as training images, and sequentially using data events devX、devY、devZScanning the sub-sections in the corresponding directions to obtain corresponding N conditional probability density functions { cpdf1,...,cpdfN}；
(9) From all the sub-sections, N conditional probability density functions { cpdf } corresponding to the current simulation position x are obtained1,...,cpdfNAnd then the prior probability obtained from step (2)P0The prior probability p of the sub-region corresponding to the current simulation position x is obtained0；
(10) Fusing the obtained N conditional probability density functions { cpdf according to a probability fusion strategy1,...,cpdfNAnd a priori probability p0Thereby obtaining a joint probability distribution;
(11) randomly extracting a corresponding attribute from the final joint probability distribution, assigning the attribute to the current simulation position x, and turning to the step (5) to continue the attribute simulation of other unknown positions;
(12) all unknown nodes on the current grid g are simulated, and the step (3) is switched to simulate the next repeated grid;
(13) all the unknown positions on each multiple grid G are simulated, the simulation result is stored, and the reconstruction is finished;
the attribute values of all nodes on the whole simulation grid are completely obtained, so that a three-dimensional geological structure model which fills the whole modeling space is obtained.
2. The local search strategy-based multipoint statistic three-dimensional geological model automatic reconstruction method according to claim 1, characterized in that in the step (8), for the scanning process of any one of the sub-sections, the following steps (8-1) to (8-8) are carried out:
(8-1) acquiring all nodes of the currently scanned sub-section according to the index of the sub-section in the whole simulation grid;
(8-2) setting a random scanning path and initializing a counter sum of the number of matched patterns to be 0;
(8-3) sequentially obtaining corresponding positions path in the sub-section from the random scanning pathiiIf ii is<If the path () × f is true, turning to the step (8-4), otherwise, turning to the step (9); where f is the proportion of a given scanned sub-section, i.e. allowing a portion of the sub-section to be scanned rather than the entirety;
(8-4) according to pathiiAcquiring a corresponding neighborhood N 'and a data event dev' at the position in the training image;
(8-5) calculating a pattern distance d { dev, dev '} between the data event dev from the simulation grid and the data event dev' from the training image using the following formula:
where k represents the number of nodes contained in the data event, Z (x)j) And Z (y)j) Respectively representing the attribute values of the central nodes corresponding to the data events in the simulation grid and the training image;
(8-6) if d { dev, dev' } t is true, turning to the step (8-7), otherwise, turning to the step (8-8); wherein t is a threshold value of the pattern distance corresponding to the current grid g, and if the current pattern distance is smaller than the threshold value, the current pattern distance is received;
(8-7) updating the conditional probability density function cpdf corresponding to the current sub-section according to the attribute value corresponding to the data event center node in the training image, and enabling sum + +;
(8-8) if sum > maximumN is true, turning to the step (9), otherwise, turning to the step (8-3) to continue scanning the current sub-section; where maximumN represents the maximum number of multi-point patterns to be acquired from the current training image.
3. The method for automatically reconstructing the multi-point statistical three-dimensional geological model based on the local search strategy as claimed in claim 2, wherein in the step (8), when scanning the sub-section, the search neighborhood radius R of the multiple grids and the threshold value t of the mode distance are allowed to be linearly decreased with the increase of the grid level; the neighborhood radius R and the threshold t on the first rescreen are given by the user, i.e. R ═ R0，t＝t0(ii) a Fixing the grid to be R1 and t 0 on the last grid; for any intermediate mesh, its neighborhood radius R and threshold t are expressed as:
where m represents the total number of multiple grids and τ (1 ≦ τ ≦ m) represents the index of the current grid.
4. The method for automatically reconstructing the multi-point statistical three-dimensional geological model based on the local search strategy as claimed in claim 2, wherein in the step (8-8), a value range of maximumN is [50,100], wherein maximumN represents the maximum number of multi-point patterns obtained from the current sub-section, and is adjusted moderately according to the scale of the model grid and the complexity of the geological structure; therefore, the stability of the conditional probability density function cpdf is controlled by controlling the number of successfully matched modes, and unnecessary calculation time is avoided.
5. The method for automatically reconstructing the multi-point statistical three-dimensional geological model based on the local search strategy as claimed in claim 2, wherein in the step (10), the probability fusion strategy comprises the following two steps:
(10-1) fusing conditional probability density functions from two sub-sections parallel to each other in the same direction using Linear Pooling formula based on probability addition:
wherein A represents a probabilistic fused attribute class, PG(A) Representing the probability, P, of the property after fusioni(A) Representing probability of property, w, from any sub-sectioniIs the corresponding weight; at most two sections are arranged in the parallel direction of one subregion, and probability fusion is needed only when two sections exist; when m is 2, weight w1The value of w is calculated by the following formula:
wherein d is1、d2Respectively representing the distances from the current simulation position to the two sub-sections;
(10-2) fusing probability distribution in three orthogonal directions processed in the step (10-1) and prior probability by using a Log-Linear Pooling formula based on probability multiplication:
wherein A represents the attribute category of probability fusion, n is the number of probability distribution participating in probability fusion, n is more than or equal to 1 and less than or equal to 3, namely the probability distribution in X, Y, Z three orthogonal directions obtained after the processing of step (10-1) at most, PG(A) Representing the probability, P, of the property after fusioni(A) Representing the probability of the attribute, P, of any sub-direction after the processing in step (10-1)0(A) Is the prior probability, w, of the propertyiWeights for the probabilities in each direction; weights w comprising a priori probabilities0Within, the sum must be 1, i.e.
Through the adjustment of each weight, the influence degree of the probability distribution in different directions on the final combined probability distribution is different, and the larger the weight in a certain direction is, the closer the final result is to the probability distribution condition; weight of prior probability
Namely when
When w0Is positive, the prior probability also participates in the probability fusion, and the weight w0The larger the effect, the correspondingly larger; when in use
When w0When the prior probability is 0, the regulation function fails; when in use
When the temperature of the water is higher than the set temperature,w0the prior probability is a negative number, and the prior probability is a negative regulation effect, namely the final probability distribution deviates from the prior probability distribution; in the weight selection process, it is necessary to ensure that the sum of all weights is 1.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810060685.7A CN108416841B (en) | 2018-01-22 | 2018-01-22 | Automatic reconstruction method of multi-point statistical three-dimensional geological model based on local search strategy |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810060685.7A CN108416841B (en) | 2018-01-22 | 2018-01-22 | Automatic reconstruction method of multi-point statistical three-dimensional geological model based on local search strategy |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108416841A CN108416841A (en) | 2018-08-17 |
| CN108416841B true CN108416841B (en) | 2021-05-28 |
Family
ID=63126030
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810060685.7A Active CN108416841B (en) | 2018-01-22 | 2018-01-22 | Automatic reconstruction method of multi-point statistical three-dimensional geological model based on local search strategy |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108416841B (en) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109671114A (en) * | 2018-12-10 | 2019-04-23 | 扬州大学 | River work measurement method based on three-dimensional reconstruction |
| CN110363299B (en) * | 2019-07-16 | 2023-04-07 | 成都理工大学 | Spatial case reasoning method for outcrop rock stratum layering |
| CN111415407B (en) * | 2020-03-27 | 2023-04-07 | 西北民族大学 | Method for improving performance of three-dimensional reconstruction image by adopting multi-template system |
| CN111709169B (en) * | 2020-05-29 | 2021-08-24 | 中国地质大学（武汉） | Multi-point geostatistical random simulation method based on conditional conduction probability |
| CN111739165B (en) * | 2020-06-15 | 2024-05-03 | 鞍钢集团矿业有限公司 | Method for locally updating three-dimensional grade model of ore body |
| CN111985800B (en) * | 2020-08-11 | 2022-11-11 | 武汉深能环保新沟垃圾发电有限公司 | Intelligent control method for garbage pool of garbage incineration power plant |
| CN113704371B (en) * | 2021-07-16 | 2024-06-28 | 重庆工商大学 | Method for adaptively detecting and dividing subareas in geographic information network |
| CN113537329B (en) * | 2021-07-30 | 2022-05-31 | 山西大学 | Method for rapidly estimating probability distribution of various ground objects position by position |
| CN114078183B (en) * | 2021-11-01 | 2023-06-20 | 清华大学 | Reconstruction method, device, equipment and medium of three-dimensional structure of porous medium |
| CN116958470B (en) * | 2023-07-25 | 2024-05-07 | 中山大学 | Geological modeling method and device integrating Markov chain and multipoint statistics |
| CN116912534B (en) * | 2023-09-14 | 2023-12-22 | 中国地质大学（武汉） | Method for identifying spatial structure of hydrothermal deposit ore-forming system by self-adaptive search matching |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101706956A (en) * | 2009-11-06 | 2010-05-12 | 上海第二工业大学 | Method for reconstructing image statistical information by multi-point geostatistics method |
| CN101706966A (en) * | 2009-11-06 | 2010-05-12 | 上海第二工业大学 | Method for three-dimensional reconstruction of porous medium on basis of two-dimensional images and multi-point statistical method |
| US8989447B2 (en) * | 2012-08-13 | 2015-03-24 | Texas Instruments Incorporated | Dynamic focus for computational imaging |
| CN107358654A (en) * | 2017-06-19 | 2017-11-17 | 中国地质大学（武汉） | Section reconstruct three-dimensional surface modeling method and system based on polygon morphing technology |
-
2018
- 2018-01-22 CN CN201810060685.7A patent/CN108416841B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101706956A (en) * | 2009-11-06 | 2010-05-12 | 上海第二工业大学 | Method for reconstructing image statistical information by multi-point geostatistics method |
| CN101706966A (en) * | 2009-11-06 | 2010-05-12 | 上海第二工业大学 | Method for three-dimensional reconstruction of porous medium on basis of two-dimensional images and multi-point statistical method |
| US8989447B2 (en) * | 2012-08-13 | 2015-03-24 | Texas Instruments Incorporated | Dynamic focus for computational imaging |
| CN107358654A (en) * | 2017-06-19 | 2017-11-17 | 中国地质大学（武汉） | Section reconstruct three-dimensional surface modeling method and system based on polygon morphing technology |
Non-Patent Citations (2)
| Title |
|---|
| The three-dimensional entity creation algorithm for surface reconstruction based on two-dimensional contour lines;Ning Chao 等;《 International Conference on Electrical Information and Mechatronics (ICEIM) 》;20111225;全文 * |
| 基于多点地质统计学岩相随机模拟研究;于明乐;《中国优秀硕士学位论文全文数据库基础科学辑（月刊 ）》;20170315;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108416841A (en) | 2018-08-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108416841B (en) | Automatic reconstruction method of multi-point statistical three-dimensional geological model based on local search strategy | |
| Shi et al. | Development of subsurface geological cross-section from limited site-specific boreholes and prior geological knowledge using iterative convolution XGBoost | |
| Tahmasebi et al. | Multiple-point geostatistical modeling based on the cross-correlation functions | |
| Arpat et al. | Conditional simulation with patterns | |
| Mariethoz et al. | The direct sampling method to perform multiple‐point geostatistical simulations | |
| Strebelle et al. | Reservoir modeling using multiple-point statistics | |
| WO2021147529A1 (en) | Multipoint geostatistical pre-stack inversion method based on updated theory of permanence of probability ratio | |
| US8775142B2 (en) | Stochastic downscaling algorithm and applications to geological model downscaling | |
| Suzuki et al. | History matching with an uncertain geological scenario | |
| US9140821B2 (en) | Ordered multipoint geostatistics simulation using non-symmetric search mask | |
| WO2017007924A1 (en) | Improved geobody continuity in geological models based on multiple point statistics | |
| WO2006023602A2 (en) | Method for making a reservoir facies model | |
| US20120059634A1 (en) | Method for optimizing the positioning of wells in an oil reservoir | |
| Chen et al. | Conditional multiple-point geostatistical simulation for unevenly distributed sample data | |
| US20230088307A1 (en) | Hierarchical Building and Conditioning of Geological Models with Machine Learning Parameterized Templates and Methods for Using the Same | |
| Gardet et al. | Pattern-based conditional simulation with a raster path: a few techniques to make it more efficient | |
| Deutsch et al. | A multidimensional scaling approach to enforce reproduction of transition probabilities in truncated plurigaussian simulation | |
| Tahmasebi | Structural adjustment for accurate conditioning in large-scale subsurface systems | |
| Liu et al. | A feature-enhanced MPS approach to reconstruct 3D deposit models using 2D geological cross sections: A case study in the Luodang Cu deposit, Southwestern China | |
| US9329289B2 (en) | Method of constructing a grid representative of a property distribution by conditional multipoint statistical simulation | |
| Koneshloo et al. | A workflow for static reservoir modeling guided by seismic data in a fluvial system | |
| Le Ravalec | Optimizing well placement with quality maps derived from multi-fidelity meta-models | |
| Alakeely et al. | Simulating oil and water production in reservoirs with generative deep learning | |
| Huang et al. | Theoretical generalization of Markov chain random field from potential function perspective | |
| WO2021160980A1 (en) | An improved micro-ct image acquisition method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |