CN112016048B - Geological attribute prediction method based on multi-point geological statistics - Google Patents
Geological attribute prediction method based on multi-point geological statistics Download PDFInfo
- Publication number
- CN112016048B CN112016048B CN202010740254.2A CN202010740254A CN112016048B CN 112016048 B CN112016048 B CN 112016048B CN 202010740254 A CN202010740254 A CN 202010740254A CN 112016048 B CN112016048 B CN 112016048B
- Authority
- CN
- China
- Prior art keywords
- geological
- point
- points
- query
- adjacent
- 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
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/22—Matching criteria, e.g. proximity measures
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Evolutionary Biology (AREA)
- Mathematical Physics (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Operations Research (AREA)
- Probability & Statistics with Applications (AREA)
- Artificial Intelligence (AREA)
- Algebra (AREA)
- Evolutionary Computation (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention relates to a geological attribute prediction method based on multipoint geological statistics, which comprises the following steps: importing observation data of each actual geological exploration point, wherein the observation data at least comprises three-dimensional coordinates and corresponding geological attributes of each actual geological exploration point; inputting the three-dimensional coordinates of the query points, and judging the spatial adjacency relation graph of the query points according to a Voronoi algorithm or an envelope circle algorithm based on Euclidean distance; confirming adjacent points of the query points, wherein the adjacent points are actual geological exploration points falling into the spatial adjacent relation graph; calculating the probability distribution of the three-dimensional geological attribute at each adjacent point; weighting and summing the probability distribution of all the adjacent points by adopting a weighting algorithm and a weight function to obtain a three-dimensional geological attribute probability distribution function at the query point; and determining the maximum likelihood estimation of the three-dimensional geological attribute at the query point according to the probability distribution function, and obtaining the geological attribute of the input point. The method has simple principle and convenient modification and maintenance, and is convenient for front-line employees to use.
Description
Technical Field
The invention relates to the technical field of three-dimensional geological modeling, in particular to a geological attribute prediction method based on multi-point geological statistics.
Background
The three-dimensional geological information is a kind of space correlation variable, and the classic space correlation variable estimation method mainly comprises a variation function method, a Krigin estimation method and variants of the two methods. The variance function is a main tool for quantitatively describing the spatial correlation of variables, and is generally represented by a curve, and the main parameters describing the curve include a variable range, a lump value, a base value, and the like. The variogram can only describe statistical properties between two points, and it is therefore difficult to characterize a geologic volume with a complex spatial structure and geometry. Common theoretical variogram models include a spherical model, an exponential model, a gaussian model, and the like. Kriging estimation is a method of performing local estimation. It provides the best estimator of the mean of the regionalized variable over a local region, i.e. the optimal (least estimate variance), unbiased (mathematical expectation of estimation error is zero) estimate.
There are some professional geological analysis business programs that develop relevant functional modules and secondary development custom interfaces for this purpose, but these analytical modeling works are highly professional and heavy, and are a serious challenge for most front-line practitioners.
Disclosure of Invention
The invention aims to provide a geological attribute prediction method based on multi-point geological statistics, which has the advantages of simple principle, convenience in modification and maintenance and convenience for front-line employees.
In order to realize the purpose of the invention, the technical scheme is as follows:
a geological attribute prediction method based on multi-point geological statistics comprises the following steps:
s1: importing observation data of each actual geological exploration point, wherein the observation data at least comprises a three-dimensional coordinate of each actual geological exploration point and a geological attribute corresponding to the three-dimensional coordinate;
s2: inputting three-dimensional coordinates of a query point, and judging a spatial adjacency relation boundary of the query point according to a Voronoi algorithm or an envelope circle algorithm based on Euclidean distance;
s3, confirming adjacent points of the query points, wherein the adjacent points are actual geological exploration points falling into the space adjacent relation boundary;
s4, calculating the probability distribution of the three-dimensional geological attribute at each adjacent point;
s5, carrying out weighted summation on the probability distribution of all adjacent points by adopting a weighting algorithm and a weight function to obtain a three-dimensional geological attribute probability distribution function at the query point;
and S6, determining the maximum likelihood estimation of the three-dimensional geological attribute at the query point according to the probability distribution function in the step S5, and obtaining the geological attribute of the input point.
Further, in step S1, observation data of each actual geological exploration point is obtained by means of drilling geophysical prospecting.
Further, in step S2, the three-dimensional coordinates of the actual geological survey point are projected onto the XY horizontal plane to form projected points, and then the adjacency relationship between the projected points is compared in the XY horizontal plane.
Further, in step S4, a search interval is set around the adjacent point to be calculated, and the probability distribution of the geological property of the adjacent point in the search interval is calculated.
Further, in step S4, a search radius σ is set, the elevation at the query point is z, and the search interval is [ z- σ, z + σ ]]Calculating the probability distribution P of the geological attribute in the adjacent point search interval ij In which P is ij ＝∑h j /2σ，∑h j Is adjacent point elevation interval [ z-sigma, z + sigma ]]The sum of the thicknesses of soil layers corresponding to similar geologies in the interior, i is any natural number from 1 to n, n is the total number of adjacent points, and j is the search interval [ z-sigma, z + sigma ] at the adjacent points]Different geological category codes.
Further, in step S4, when the search interval exceeds the borehole revealing interval, the geological class code of the excess portion is identified as "null".
Further, in step S5, an inverse distance weighting algorithm is adopted, and the weight function is w i ＝D/d i Wherein d is i The projection distance from the ith adjacent point to the query point on the XY horizontal plane,
probability distribution P of different geological classes at a query point j ＝∑w i P ij 。
Further, in step S6, the largest is selectedProbability distribution P j And outputting the corresponding geological category code to obtain the geological attribute of the query point.
Compared with the prior art, the invention has the following beneficial effects:
the invention not only comprehensively considers the spatial distribution relation of geological exploration data to ensure that the prediction result is reliable and credible, but also adopts a modular design, and the implementation method of each step can be independently called and replaced by functions, is convenient for modification and maintenance, can meet the service requirements of general engineering technicians, and can be used for testing and developing front-edge research methods by deep-funding researchers.
Drawings
FIG. 1 is a schematic flow chart of a geological property prediction method based on multi-point geostatistics according to an embodiment of the invention;
FIG. 2 is a schematic diagram of a Voronoi algorithm applied to determine the adjacency relationship in a data point plane according to an embodiment of the present invention;
FIG. 3 is a schematic diagram of a probability distribution calculation of geological borehole soil layering according to an embodiment of the present invention.
Description of the reference numerals:
10. actual geological survey point, 20 query point, 30 abutment point, 40 box, 50 circle, 60 top, 70 ground, 80 aperture.
Detailed Description
The embodiments of the invention are described in detail below with reference to the drawings:
as shown in fig. 1 to 3, a geological property prediction method based on multipoint geostatistics includes the following steps:
s1: importing observation data of each actual geological exploration point, wherein the observation data at least comprises a three-dimensional coordinate of each actual geological exploration point and a geological attribute corresponding to the three-dimensional coordinate;
s2: inputting the three-dimensional coordinates of the query points, and judging the spatial adjacency relation boundary of the query points according to a Voronoi algorithm or an envelope circle algorithm based on Euclidean distance;
s3, confirming an adjacent point 30 of the query point 20, wherein the adjacent point 30 is an actual geological exploration point 10 falling into the space adjacent relation boundary;
s4, calculating the probability distribution of the three-dimensional geological attribute at each adjacent point 30;
s5, carrying out weighted summation on probability distribution of all the adjacent points 30 by adopting a weighting algorithm and a weight function to obtain a three-dimensional geological attribute probability distribution function at the query point 20;
and S6, determining the maximum likelihood estimation of the three-dimensional geological attribute at the query point 20 according to the probability distribution function in the step S5, and obtaining the geological attribute of the input point.
In the embodiment, the soil texture prediction is based on geological drilling, the input data are three-dimensional coordinates of a drilling hole opening 80, observed soil stratification data and three-dimensional coordinates of the query point 20, and the output data are maximum likelihood estimation of soil texture categories at the query point 20. The method can also be based on the soil property prediction of the deep well according to the actual requirement.
In step S1, the data in the drilling exploration point floor plan map is recorded in the format of table 1, the drilling observation information in the drilling histogram is recorded in the format of table 2, and a code is added to the geological formation according to the figure number and the formation name (when the combination of the figure number and the formation name is completely consistent, the geological formation is regarded as the same geology), and the geological formation is imported into an analysis program, and the analysis program can be written in different programming languages, such as Python, MATLAB and the like, according to actual needs. The top 60, bottom 70, and orifice 80 positions in table 2 are shown in fig. 3.
TABLE 1 borehole survey Point numbering and Orifice coordinates
TABLE 2 geological drilling observations
In step S2, different thiessen polygon maps are formed according to different query points, as shown in fig. 2. Due to the equipartition characteristic of the Thiessen polygon on the space subdivision, the adjacency relation of the drilling exploration points in the space can be analyzed by adopting a Voronoi function, so that the adjacency drilling point distribution of any drilling point is obtained, as shown in FIG. 2, the spatial adjacency relation boundary judged by Voronoi is a block 40, and the points in the block 40 are the Voronoi adjacency points 30 of the query point 20, wherein the total number is 7; the space adjacency relation boundary judged by the Euclidean distance enveloping circle algorithm is a circle 50, and the points in the circle 50 are Euclidean adjacency points 30 of the query point 20, and the number of the points is 10.
In step S3, a certain adjacent drilling point of the query point 20 is selected, as shown in fig. 3, the search radius σ is set, and the elevation at the query point 20 is z, so that only the search interval [ z- σ, z + σ ] of the adjacent drilling point is examined]The distribution situation of geological categories in the hole can determine the soil property probability distribution P at the drilling point ij . When the search interval exceeds the drill hole revealing interval, the soil property type of the exceeding portion is specified as "empty" and is marked as "null". For the vertical drilling hole shown in the figure 3, the soil layer information is inquired from top to bottom, and the sum sigma h of the soil layer thicknesses corresponding to the same geology is used j Divided by the total length of the search interval, 2 σ, then P ij ＝∑h j I is any natural number from 1 to n, n is the total number of the adjacent points 30, and j is the search range [ z- σ, z + σ ] at the adjacent points 30]Different geological category codes are included. As shown in fig. 2, if a Voronoi algorithm is adopted, i is any one of natural numbers 1 to 7; if an envelope circle algorithm based on Euclidean distance is adopted, i is any natural number from 1 to 10. Table 3 shows the geological distribution probability in the search interval at a certain adjacency point 30, where the geological class of the adjacency point 30 in the corresponding search interval has 3 types, and the corresponding geological class codes are a, B, and C, respectively.
TABLE 3 probability of geological distribution within search interval at certain adjacency 30
In step S4, an inverse distance weighting algorithm is selected, and the weight function is w i ＝D/d i Wherein d is i The projection distance from the ith adjacent point 30 to the query point 20,
the probability distribution P at query point 20 may be denoted P i ＝∑w i P ij . As can be seen from table 4, the geological categories of all the adjacent points 30 in the corresponding search interval have 5 categories in total, and the corresponding geological category codes are a, B, C, D, null, respectively.
TABLE 4 probability of geological distribution at query point 20
In step S5, a geological category with the highest probability is selected and output, and the geology is the maximum likelihood estimation of the soil texture state at the query point 20. For the data in Table 4, the output geological code is "B", the geological name is [ Lengend _3, name _2], and the corresponding probability is 0.36.
The three-dimensional geological attribute rapid prediction method based on multi-point geological statistics obtains the maximum likelihood estimation of the three-dimensional geological attribute at the designated point on the region to be modeled by analyzing the adjacent observation data of the region to be modeled and combining the spatial distribution characteristics of the data points. The method not only comprehensively considers the spatial distribution relation of geological exploration data to ensure that a prediction result is reliable and credible, but also adopts a modular design, and the implementation method of each step can be independently called and replaced by a function, so that the modification and maintenance are convenient, the service requirements of general engineering technicians can be met, and the advanced research method can be tested and developed by deep researchers.
The technical features of the embodiments described above may be arbitrarily combined, and for the sake of brevity, all possible combinations of the technical features in the embodiments described above are not described, but should be considered as being within the scope of the present specification as long as there is no contradiction between the combinations of the technical features.
The above-mentioned embodiments only express several embodiments of the present invention, and the description thereof is more specific and detailed, but not construed as limiting the scope of the invention. It should be noted that, for a person skilled in the art, several variations and modifications can be made without departing from the inventive concept, which falls within the scope of the present invention. Therefore, the protection scope of the present patent should be subject to the appended claims.
Claims (8)
1. A geological attribute prediction method based on multipoint geological statistics is characterized by comprising the following steps:
s1: importing observation data of each actual geological exploration point, wherein the observation data at least comprises a three-dimensional coordinate of each actual geological exploration point and a geological attribute corresponding to the three-dimensional coordinate;
s2: inputting the three-dimensional coordinates of the query points, and judging the spatial adjacency relation boundary of the query points according to a Voronoi algorithm or an envelope circle algorithm based on Euclidean distance;
s3, confirming adjacent points of the query points, wherein the adjacent points are actual geological exploration points falling into the space adjacent relation boundary;
s4, calculating the probability distribution of the three-dimensional geological attribute at each adjacent point;
s5, carrying out weighted summation on the probability distribution of all adjacent points by adopting a weighting algorithm and a weight function to obtain a three-dimensional geological attribute probability distribution function at the query point;
and S6, determining the maximum likelihood estimation of the three-dimensional geological attribute at the query point according to the probability distribution function in the step S5, and obtaining the geological attribute of the input point.
2. The multipoint-geostatistical-based geologic property prediction method of claim 1, wherein in step S1, the observation data of each actual geologic prospecting point is obtained by means of drilling geophysical prospecting.
3. The multipoint-geostatistical-based geologic property prediction method of claim 2, wherein in step S2, the three-dimensional coordinates of the actual geologic survey point are projected onto an XY horizontal plane to form projected points, and then the adjacency relations between the projected points are compared within the XY horizontal plane.
4. The multipoint-geostatistical-based geological property prediction method according to claim 2, wherein in step S4, a search interval is set around the adjacent point to be calculated, and the probability distribution of the geological property of the adjacent point in the search interval is calculated.
5. The multipoint-geostatistical-based geological property prediction method according to claim 4, wherein a retrieval radius σ is set, the elevation at the query point is z, and the retrieval interval is [ z- σ, z + σ ]]Calculating probability distribution P of geological attribute in adjacent point search interval ij In which P is ij ＝∑h j /2σ，∑h j Search interval [ z-sigma, z + sigma ] for neighbor points]The sum of the thicknesses of soil layers corresponding to similar geology in the interior, i is any natural number from 1 to n, n is the total number of adjacent points, and j is the search interval [ z-sigma, z + sigma ] at the adjacent points]Different geological category codes.
6. The method of predicting geologic properties based on multi-point geostatistics of claim 5, wherein when the search interval exceeds the borehole reveal interval, the exceeding part of the geologic type code is identified as "null".
7. The method for predicting the geological properties based on the multi-point geostatistics of claim 5, characterized in that in the step S5, an inverse distance weighting algorithm is adopted, and the weight function is w i ＝D/d i Wherein d is i The projection distance from the ith adjacent point to the query point on the XY horizontal plane,
different geological categories at query pointsProbability distribution P of j ＝∑w i P ij 。
8. The multipoint-geostatistical-based geologic property prediction method of claim 7, wherein in step S6, the largest probability distribution P is chosen j And outputting the corresponding geological category code to obtain the geological attribute of the query point.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010740254.2A CN112016048B (en) | 2020-07-28 | 2020-07-28 | Geological attribute prediction method based on multi-point geological statistics |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010740254.2A CN112016048B (en) | 2020-07-28 | 2020-07-28 | Geological attribute prediction method based on multi-point geological statistics |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112016048A CN112016048A (en) | 2020-12-01 |
| CN112016048B true CN112016048B (en) | 2022-10-18 |
Family
ID=73499980
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010740254.2A Active CN112016048B (en) | 2020-07-28 | 2020-07-28 | Geological attribute prediction method based on multi-point geological statistics |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112016048B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113537329B (en) * | 2021-07-30 | 2022-05-31 | 山西大学 | Method for rapidly estimating probability distribution of various ground objects position by position |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3121622A1 (en) * | 2015-07-24 | 2017-01-25 | Bergen Teknologioverforing AS | Method of predicting parameters of a geological formation |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20200033495A1 (en) * | 2018-07-30 | 2020-01-30 | Cgg Services Sas | Method and apparatus for obtaining statistical measures of geological properties values related to lateral wells using seismic-derived maps |
| CN109633781B (en) * | 2018-08-16 | 2021-03-02 | 清能艾科（深圳）能源技术有限公司 | Geological property acquisition method and device, electronic equipment and storage medium |
| CN111311746B (en) * | 2019-12-09 | 2023-09-08 | 中交广州航道局有限公司 | Intelligent three-dimensional geological modeling method based on drilling data |
-
2020
- 2020-07-28 CN CN202010740254.2A patent/CN112016048B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3121622A1 (en) * | 2015-07-24 | 2017-01-25 | Bergen Teknologioverforing AS | Method of predicting parameters of a geological formation |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112016048A (en) | 2020-12-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Mariethoz et al. | Reconstruction of incomplete data sets or images using direct sampling | |
| Alnahwi et al. | Mineralogical composition and total organic carbon quantification using x-ray fluorescence data from the Upper Cretaceous Eagle Ford Group in southern Texas | |
| AU2018211194A1 (en) | Context based geo-seismic object identification | |
| US8234073B2 (en) | System and method for estimating geological architecture of a geologic volume | |
| AU2021208450B2 (en) | Correlation of multiple wells using subsurface representation | |
| US9146903B2 (en) | Method of using spatially independent subsets of data to calculate vertical trend curve uncertainty of spatially correlated reservoir data | |
| Eskandari et al. | Reservoir Modelling of Complex Geological Systems--A Multiple-Point Perspective | |
| WO2014099201A1 (en) | Geophysical modeling of subsurface volumes based on horizon extraction | |
| KR20210150917A (en) | Method and apparatus for estimating lithofacies by learning well log | |
| CN112016048B (en) | Geological attribute prediction method based on multi-point geological statistics | |
| CN108595749A (en) | A kind of resource reserve appraisal procedure using variation function single direction structural analysis | |
| CN114036829A (en) | Geological profile generation method, system, equipment and storage medium | |
| Prades | Geostatistics and clustering for geochemical data analysis | |
| CN112083144B (en) | Fault on-off prediction method and device, computer equipment and storage medium | |
| CN112800518B (en) | Stratum surface model correction method utilizing adjacent stratum cross-correlation constraint | |
| CN114707597A (en) | River facies tight sandstone reservoir complex lithofacies intelligent identification method and system | |
| US20220035070A1 (en) | Determination of representative wells to characterize subsurface regions | |
| CN114581556A (en) | Digital map filling method in regional geological survey | |
| CN112987091A (en) | Reservoir detection method and device, electronic equipment and storage medium | |
| Caers | Geostatistics: from pattern recognition to pattern reproduction | |
| CN118194162B (en) | Method, system, electronic equipment and storage medium for locating mining target area based on multivariate data | |
| Rawlinson et al. | Gaussian process modeling of well logs | |
| La Pointe et al. | Geological discrete fracture network model for the Laxemar site. Site Descriptive Modelling. SDM-Site Laxemar | |
| Silva Maureira | Enhanced geologic modeling with data-driven training images for improved resources and recoverable reserves | |
| CN117407841B (en) | Shale layer seam prediction method based on optimization integration algorithm |
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 |