CN111968231A

CN111968231A - Three-dimensional stratum modeling method based on geological map section

Info

Publication number: CN111968231A
Application number: CN202010819474.4A
Authority: CN
Inventors: 简兴祥; 王丽坤; 钟红梅; 李宏涛; 罗海金; 袁宏; 马婵华; 姚华龙; 叶娇珑
Original assignee: Sichuan Institute Of Nuclear Industry Geology Southwest Institute Of Nuclear Industry Geology; Sichuan Jinhe Mining Co ltd; Chengdu Univeristy of Technology
Current assignee: Sichuan Institute Of Nuclear Industry Geology Southwest Institute Of Nuclear Industry Geology; Sichuan Jinhe Mining Co ltd; Chengdu Univeristy of Technology
Priority date: 2020-08-14
Filing date: 2020-08-14
Publication date: 2020-11-20
Anticipated expiration: 2040-08-14
Also published as: CN111968231B

Abstract

The invention discloses a three-dimensional stratum modeling method based on a geological map section, which comprises the following steps: s1, editing the section of the graph to complete the setting of basic parameters required by modeling; s2, carrying out discrete sampling on the stratum boundary and the inner and outer regions to obtain the discrete distribution of stratum attributes in space; s3, carrying out interpolation calculation on the discrete points of the stratum three-dimensional space to generate a three-dimensional space stratum attribute data field; s4, extracting an isosurface from the three-dimensional space stratum attribute data field to obtain a three-dimensional contour surface of the stratum; s5, repeating the steps S2-S4 to obtain the three-dimensional contour surfaces of all the stratums; and S6, spatially splicing the three-dimensional contour surfaces of all the stratums to complete three-dimensional stratum modeling. The method avoids the problems of complex contour matching, bifurcation, corresponding triangular net connection and the like in a two-dimensional contour reconstruction method, and can obtain a three-dimensional vector geological model with higher precision; the method is suitable for stratum modeling under any complex condition, and has high modeling efficiency and wide application range.

Description

Three-dimensional stratum modeling method based on geological map section

Technical Field

The invention relates to the fields of geology, mineral products and engineering, in particular to a three-dimensional stratum modeling method based on a geological map cutting section.

Background

Three-dimensional geological modeling is always a technical hotspot and difficult problem for research in engineering, geology and mineral industry. The geological structure shape and the rock physical property parameters under the limited geological data recovery site play a very important role in engineering geological analysis, oil gas and mineral exploration and development. Due to the diversity and the multi-solution of geological data and the limitation of data acquisition means, the three-dimensional geological modeling has high difficulty and complex realization process. The three-dimensional geological modeling by using a geological map cut section is a common modeling mode.

The figure cutting section modeling is to draw a series of parallel (or approximate parallel) cutting lines on a geological map, divide a three-dimensional geological target into a series of two-dimensional geological sections by adopting a geological analysis method, and then construct a three-dimensional target geological body model by adopting a three-dimensional reconstruction method. The technology of modeling by using a cutting profile mainly has two modes, one mode is a three-dimensional reconstruction technology based on a two-dimensional contour; one is a three-dimensional mesh modeling technique that rasterizes a two-dimensional profile.

The basic idea of the three-dimensional surface reconstruction technology based on the two-dimensional contour is to intercept a series of two-dimensional parallel sections from a three-dimensional object, then extract the boundary contour line of each section aiming at a specific target, and perform matching and corresponding triangular connection on adjacent contour lines in a three-dimensional space, thereby completing the modeling of the surface of the three-dimensional object. The gridding modeling based on the two-dimensional section is to grid the two-dimensional section according to the attribute classification of the geologic body and then map the two-dimensional section into a three-dimensional space to generate a three-dimensional gridding geologic model. Compared with a gridded three-dimensional geological model, the three-dimensional surface reconstruction method based on the two-dimensional contour is a three-dimensional vector data model, has high model precision and good visualization effect, and is the most common three-dimensional modeling technology.

In the research of reconstructing a three-dimensional surface based on parallel contour lines, only the problem of contour line reconstruction on two adjacent parallel planes is generally considered. Since the reconstruction problem of multi-layer contours can be simplified to the superposition of reconstructed surfaces of a plurality of adjacent parallel contours, the connection of contours of adjacent sections is the key to complete modeling. The connection of contour lines mainly comprises three problems of contour correspondence, contour splicing and bifurcation, wherein the problems of contour correspondence and bifurcation are difficult problems. The existing algorithms have many limitations on adjacent contour lines, such as the requirement of small space between adjacent tangent planes, high coverage between corresponding contour lines, similar shape, etc., while the algorithms related to the branching problem focus on solving the problem of correct connection between a contour line and a plurality of contour lines.

If the number of contours on the same layer exceeds one, the correspondence of the contours (which contours on different layers belong to the same surface) must be determined. If the distance between layers is extremely small, the corresponding problem can be solved by the coincidence degree of the contour lines^[1](ii) a However, when the distance is smallWhen the coincidence degree is large and the correct result cannot be obtained, the corresponding problem becomes a problem which is difficult to solve. Even if the correspondence problem is solved, the branching problem needs to be solved when splitting (whether one contour line on one layer and a plurality of contour lines on the other layer belong to the same surface is determined by checking and analyzing the contour line correspondence problem, and one contour line and the plurality of contour lines are connected into an object surface through a triangular patch by a branching structure). Many documents discuss the branch problem and the corresponding problem, but no satisfactory solution has been found yet^[2][3][4]. Barequet proposes an optimized triangulation strategy based on a dynamic programming algorithm, namely a BPLI (Barequtps PIecewise Linear interpolation) method^[5][6]. For two-layer contours of the input, the method does not directly address the branching problem and the corresponding problem, but tries to join them into a three-dimensional surface that does not appear to self-intersect. The method can obtain the subdivision result with correct topological structure for the complex surface, thereby being capable of ingeniously bypassing the corresponding problem and the branch problem.

How to determine a desired combination from among a plurality of acceptable surface combinations is currently performed by: method for maximizing reconstructed three-dimensional body volume^[7]Surface area minimization method^[8]Side length shortening method^[9]Shortest diagonal method^[10]Synchronous advancing method of adjacent contour lines^[11]Construction algorithm based on domain knowledge^[12]Cutting and sewing method^[13]And the like. For solving the problem of reconstruction of adjacent contour lines with large difference of shapes and vertex numbers, the methods have different limitations in the aspects of adaptability, connection effect and efficiency.

He jin nation^[14](2003) The BPLI is improved, tried and researched, and a contour line segment matching algorithm and a new algorithm for solving the space polygon triangulation are provided, so that a degradation area is eliminated, and the matching effect is improved; qiweili wine^[15](2007) A method of inserting a new contour line in the intermediate layer is proposed, which can solve the branch-to-branch problem under simple conditions; what is afraid^[16](2010) On the basis of minimum surface areaOn the basis of the maximum volume method, an Optimal Path (OP) algorithm is provided, and optimization attempt is performed on contour connection; chengxiao tea^[17](2011) And the method for evolving the polygon curve solves the problems of vertex matching and connection of adjacent contour lines by simplifying the polygon and extracting the characteristic points of the polygon. However, this method is also directed to the connection problem corresponding to a single contour, and a long and narrow triangle may appear for a more complex contour; graceful thread^[18](2015) The method introduces that the reconstruction of the three-dimensional shape of the geologic body with a single contour is realized by a Discrete Smooth Interpolation (DSI) algorithm by means of GOCAD software; yangyang^[19](2015) And the like, a mode of inserting transition contour lines between contour lines is provided so as to improve the modeling quality. Jingyongbing^[20](2016) And the contour line veneering algorithm is provided, wherein a wavelet filter bank is utilized to simplify contour lines into low-resolution contour lines similar to the current situation of original contour lines, then the low-resolution contour lines are optimally veneered, and then connecting triangles between the contour lines are constructed. The algorithm reduces the complexity of the directed loop graph algorithm to a certain extent, but for adjacent contour lines with large differences, similar veneering is difficult to obtain through filtering.

Due to the fact that the modeling topological structure of the BPLI method is complex, for complex contour structures, particularly for complex contour data with large slice distances, intersected three-dimensional surface models are likely to occur, and correct three-dimensional model results cannot be obtained. These problems can be avoided by using a rasterized modeling approach.

The rasterization modeling implementation technology is simple, only the two-dimensional section needs to be rasterized, then a three-dimensional space grid model is generated through a three-dimensional interpolation algorithm, and finally the required stratum grid is extracted from the grid model. The three-dimensional geologic body in the grid model is expressed by adopting a block model, and the three-dimensional geologic body is dispersed into a regular geometric body in space (such as a hexahedral model, a tetrahedral model, a prismatic model and the like). The modeling method does not need to consider complex topological structure and triangular connection in space, is simple to realize, and is commonly used for mineral modeling, reserve estimation and the like. However, the method needs a large number of two-dimensional profile sections, the workload of the graph cut sections is large, and the visualization effect of the three-dimensional model generated by the method is poor.

Wangdan pill^[21](2012) The method comprises the steps of discussing the utilization of a two-dimensional section, rasterizing, mapping into a three-dimensional space and generating a three-dimensional grid model by using an octree; zhou Liang Chen^[22](2013) The method combines cubic spline fitting and manual editing to process complex structures such as folds and the like by adopting encrypted occurrence data of a research area so as to improve the quality of a grid model; zhang Bao Yi^[23](2017) The process of constructing a three-dimensional mesh model by gridding the map-cut sections using a large number of parallel map-cut sections; plum fruit juice^[24](2108) The method comprises the steps of generating a distance field of a two-dimensional profile by using a distance function, generating a three-dimensional attribute field by using a B-spline interpolation method, and generating the contour surface of the geologic body by using a method of extracting a physical equivalent surface by using a Marching cube algorithm on the basis. The method avoids the problem of corresponding connection of complex contour lines, and is suitable for geological modeling under various complex conditions. However, the method is difficult to find an accurate boundary value through a three-dimensional physical equivalent surface generated by the distance field, so that the boundary error of the constructed three-dimensional model is large, and the result is inaccurate. Second, because the distance field of a profile is a uniformly subdivided mesh, the amount of mesh data for the profile can be large. If the three-dimensional attribute field is calculated by adopting global interpolation, a large amount of computer memory is consumed. The distance matrix is calculated by dividing the grid according to 200 multiplied by 200, 20 sections and performing inverse distance weighted interpolation calculation, and the memory consumption of the computer can reach 2384G

Although this situation can be avoided by using the local interpolation algorithm, the local interpolation method (or the linear interpolation method) has obvious limitations. The local interpolation method does not take into account the spatial overall relevance of the data distribution. Under the condition that two-dimensional contour lines are sparse, the modeling result is fractured between the contour lines, and the spatial distribution is discontinuous.

In summary, in the existing method for performing three-dimensional modeling based on a two-dimensional profile, a three-dimensional surface reconstruction method based on a two-dimensional profile can obtain a high-quality three-dimensional geological model. But the realization technology is complex, and the requirement on two-dimensional contour lines is high. The corresponding contour lines require simple convex polygons, and non-simple polygons need polygon simplification processing, so that the modeling workload is increased. However, for the connection problems of complex multi-contour correspondence, large difference, triangle intersection, long and narrow triangle and the like of the contour lines far away, the prior art still has no good solution. The cross-section grid modeling method based on rasterization can avoid the problems and is simple in implementation technology. However, because the rasterized mesh model is influenced by subdivision, the model precision is low, the visualization effect is poor, and the application has certain limitations. And the method for modeling by extracting the isosurface from the three-dimensional attribute field has less research. The existing method has the problems of low profile sampling precision, large sampling data volume, difficult interpolation of a three-dimensional attribute field, inaccurate isosurface extraction and the like.

Reference documents:

[1]Meyers D,Skinner S.Surfaces from Contours[J].ACM Transactions on Graphics.1992,11(3):228～258.

[2]Klein R,Schilling A,StraBer W.Reconstruction and Simplification of Surfaces from Contours[J].Graphical Models.2000,62:429～443.

[3]Barequet G,Shapiro D,Tal A.Multilevel Sensitive Reconstruction of Polyhedral Surfaces from Parallel Slices[J].The Visual Computer.2000,16:116～133.

[4]Levin D.Multidimensional Reconstruction by Set 2valued Approximation[J].IMA J Numer Anal.1986,6:173～184.

[5]Barequet G,SharirM.Piecewise-linear Interpolation between Polygonal Slices[C].Proc10th Ann ACM Symp onComputational Geometry(SoCG).1994.

[6]Barequet G,SharirM.Piecewise-linear Interpolation between Polygonal Slices[C].ComputerVision ImageUnderstanding.1996.63:251～272.

[7]Keppel E.Approximating complex surfaces by triangulationof contour lines[J].IBM Journal of Research and Development.1975,19(1):2-11.

[8]Fuchs H,Kedem Z M,Uselton S P.Optimal surfaces reconstruction from planar contours[J].Communication of the A CM.1977,20(10):693-702.

[9]Christiansen H N,Sederberg T W.Conversion of complex contour line definitions into polygonal element mosaics[J].Computer Graphics.1978,12(3):187-192.

[10] liu just, Ware, Deng Lin. Algorithm discussion of layered three-dimensional body surface reconstruction [ J ]. University of Oncology, Nature Science Edition,2003,30(5):537- & 540 (Liu Gang, Hu Yuan2lai, Deng Lin. discussion on the understanding of the treatment 32D surface [ J ]. Journal of Chengdu University of Technology: Science & Technology Edition,2003,30(5):537 540.).

[11] Three-dimensional data field visualization of Tangjiang [ M ]. Beijing, Qinghua university Press.2000.

[12] Wanengde, Li Yan, Bayubin, the three-dimensional visual modeling technology of ore body [ J ]. university of northeast, proceedings of Nature sciences, 2005,26(9): 890-.

[13] A new multi-contour line reconstruction three-dimensional body algorithm of Marhon, Guomeing, a cutting-stitching method [ J ]. university of northeast proceedings.2007, 1.Vol 28). No.1.

[14] Helian, Charbins, a novel algorithm for reconstructing a three-dimensional surface from two-dimensional parallel contour lines based on BPLI [ J ]. the university of Beijing (Nature science edition): 2003(03)

[15] Qiweili, qinxingqiang, wangxi, songliping, an algorithmic study [ C ] to reconstruct a three-dimensional surface based on two-dimensional parallel contour lines, the 18 th academic conference on computer technology and applications (CACIS) nationwide, 2007.8, p891-895.

[16] Wu bird based on cross-section contour line human-computer interaction three-dimensional geologic body modeling [ J ], geophysical prospecting technique, 2010,32(4): P433-436.

[17] The method is characterized by comprising the following steps of ageing green, enabling jasper, Zhangguojian and Dingyu, a matching algorithm [ J ] of adjacent sections of complex ore bodies, college newspaper of northeast China (Nature science edition), 2011,32(4): 579-.

[18] The geologic body three-dimensional modeling based on a planar geological map [ J ] geological prospecting discourse, 2015,30(4): 594-.

[19] Yang, Pan 25035, Wu Zhen Yu, etc. a new contour line three-dimensional geological surface reconstruction method [ J ] in the scientific bulletin of terrestrial information, 2015, 17(3):253-259.

[20] A method for reconstructing the two-dimensional contour line of the three-dimensional model of complex ore body in Jingyongbing Wang Gong Zhong Sun light [ J ] metal mine 2016,285(11):124-127.

[21] Wandan three-dimensional geological modeling method based on planar geological map studies [ D ].2012, master thesis of university in Nanjing.

[22] Zhouqing, Linbingxian, Wangdan, etc. three-dimensional geologic body modeling method research of planar geological map [ J ]. scientific academic on Earth information, 2013,15(1):46-54.

[23] Zhang Bao Yi, Yangli, Chen laugh, etc. three-dimensional modeling and resource evaluation of regional mineralizing geological body based on map-cutting geological section, taking manganese ore in southwest area of Guixi as an example [ J ] the university of Jilin, university of Earth's republic of China (geoscience edition), 2017,47(3): 933-.

[24] Li Chen, three-dimensional complex geological modeling technology based on profile study [ D ]. Master thesis of university of Chengdu Physician, 2018.

[25] Xiayanhua, Baishiwei, level set-based complex three-dimensional stratigraphic model modeling and application research in underground engineering [ J ]. geotechnics, 2012,33(5): 1445-.

[26] Lie-range geological profiling based on MT algorithm [ J ] digital communications 2014,41(4):35-41.

[27] Handsome, a kind of improved MC three-dimensional reconstruction algorithm [ J ]. Chinese digital medicine, 2016,11(3):83-86.

[28] Wang, Li Ruimin, moving cube algorithm face ambiguity problem research [ J ] software engineering, 2017, 20(9):6-9.

[29] Wang. three-dimensional surface modeling in a section data field [ J ] software engineering, 2018,21(3):10-15.

Disclosure of Invention

Aiming at the problems of the existing modeling method, the invention provides a three-dimensional stratum modeling method based on a geological map section, which combines the rasterization modeling with the surface reconstruction of a geological model, thereby not only avoiding the problems of complex contour correspondence, triangulation network connection and the like, but also obtaining a three-dimensional vector geological model with higher precision; and the problems of sampling precision and data quantity in the modeling process, calculation of a three-dimensional data field by a global interpolation method, the precision of isosurface extraction and the like are solved.

In order to achieve the purpose, the technical scheme adopted by the invention is as follows:

a three-dimensional stratum modeling method based on geological map cutting profiles comprises the following steps:

s1, editing the section of the graph to complete the setting of basic parameters required by modeling; the process of editing the graph cut section comprises the following steps: (1) converting the plane coordinates of the section into three-dimensional coordinates of a geological space by calibrating a space positioning point; (2) obtaining formation parameter information through formation attribute filling; (3) obtaining the contour boundary of the stratum through stratum boundary tracking;

s2, carrying out discrete sampling on the stratum boundary and the inner and outer regions, respectively setting attribute values of respective sampling points inside the stratum, on the stratum boundary and outside the stratum, converting all the sampling points of the graph cut section into space three-dimensional coordinate points, and combining the space three-dimensional coordinate points into a space discrete point set;

s3, carrying out interpolation calculation on the three-dimensional space discrete points to generate a three-dimensional space stratum attribute data field;

s4, extracting an isosurface from the three-dimensional stratum attribute data field to obtain a three-dimensional contour surface of the stratum;

s5, repeating the steps S2-S4 to obtain the three-dimensional contour surfaces of all the stratums, and then executing the step S6;

and S6, spatially splicing the three-dimensional contour surfaces of all the stratums to complete three-dimensional stratum modeling.

Specifically, in step S1, the formation is filled by using a fast flooding filling algorithm based on a linear queue, and boundary points of the filling algorithm are recorded, so as to obtain a formation polygon and obtain formation parameter information.

Specifically, in step S1, for the formation boundary points obtained by tracking, projection processing is adopted to eliminate the boundary point tracking error due to the influence of the raster image display resolution, and ensure the consistency between the formation contour boundary obtained by tracking and the original formation contour boundary.

Further, the projection processing comprises the following steps:

(1) simplifying the stratum contour line, deleting redundant points on the contour straight line, and then forming a contour line segment { L ] by every two points_kN, where N is the number of contour line segments;

(2) selecting a tracked stratum boundary point P_i(i ═ 1,2, …, M), the set of distances { D) from that point to all contour segments is calculated_ik}(k＝1，2，..N)；

(3) From the set of distances { D_ikFind the line segment L of the minimum distance in_jDividing the boundary point P_iProjected to line segment L_jTo obtain P_ijA 1 is to P_iReplacement by P_ijAs new stratigraphic boundary points;

(4) repeating the steps (2) and (3), and projecting all the stratum boundary points onto the stratum contour lines;

(5) and simplifying the contour line of the projected stratum boundary, deleting redundant points on the straight line, and finishing the tracing of the stratum boundary.

Specifically, the step S2 includes the steps of:

s2-1, setting a sampling interval step, carrying out boundary point sampling along the boundary of the stratum contour line to obtain a stratum boundary point set B1, and setting attribute values of all points in the set to be 1;

s2-2, setting a subdivision grid M multiplied by N with a specified size, performing two-dimensional grid subdivision on a graph cutting section, and recording the coordinates of the center point of each grid to obtain a subdivision point set B2;

s2-3, carrying out stratum inner and outer judgment on the points in the set B2, if the points are in a stratum polygon, setting the attribute value of the point to be 1, and otherwise, setting the attribute value to be 0;

s2-4, performing space coordinate conversion on all the points in the sets B1 and B2, and converting the plane coordinate points into geological space coordinate points;

s2-6, putting the converted spatial coordinate points and the attribute values thereof into a scattered point data set {

p

_i1,2, …, N, after obtaining formation sampleAnd (4) finishing sampling according to the discrete space point data.

Preferably, in step S3, a three-dimensional discrete point spatial interpolation is performed by using an inverse multivariate quadratic radial basis function as an interpolation function and using a radial basis function interpolation algorithm to generate a three-dimensional spatial formation property data field; the radial basis function form is as follows:

wherein Z₀Is the interpolated point value to be estimated, w_iIs the weight of the corresponding point, | | p_i-p₀Is p_iAnd p₀The distance between them;

the interpolation calculation process is as follows:

(1) calculate distance matrix { D_ij}(i＝1,2,…,N,j＝1,2,…,N)，

(2) Calculating an inverse distance matrix D by adopting a Gaussian elimination method and combining an OpenCL multithreading acceleration technology^-1；

(3) A matrix W of weight coefficients is calculated,

Z_jis the jth discrete point P_jA formation property value of;

(4) calculating field values of properties, field values at arbitrary positions in three-dimensional space

Preferably, in step S4, the Marching Cubes algorithm is used to extract the iso-surface from the three-dimensional property data field of the formation as the three-dimensional surface contour of the formation, and the process is as follows:

(1) selecting 8 adjacent points from the three-dimensional grid data to form a hexahedron;

(2) giving an isosurface value V0, judging the state value Vi of each node in the hexahedron, and if Vi > V0, setting the state of the node as 1; otherwise, setting the value to 0; 1,2, …, 8;

(3) according to the node state sequence, binary bytes are combined in sequence to obtain the index value T of the hexahedron unit configuration table_i(0≤T_i≤255)；

(4) Index value T according to configuration_iAnd retrieving a corresponding edge node state combination from the configuration table, and calculating the intersection point coordinate P (x, y, z) of the isosurface with the edge node states of 0 and 1 on the edge, wherein the calculation formula is as follows:

p(x，y，z)＝p_1(x，y，z)+(v-v1)/(v2-v1)*(p_2(x，_y，z)-p_1(x，y，z))

in the formula, p (x, y, z) is the coordinate of the intersection point to be calculated, and v is the isosurface value to be calculated; p1, p2 are two-node coordinates of edges, v1 and v2 are two-node attribute values;

(5) indexing the value T according to the configuration table_iRetrieving corresponding triangle vertexes and connection sequences from the configuration table, and combining the triangle vertexes and the connection sequences into a closed triangle surface;

(6) and (5) traversing all the grid units, and repeating the steps (1) to (5) to obtain all the triangular surfaces of the surface outline of the target stratum.

Specifically, the step S6 includes the steps of:

(1) combining all vertexes, and performing repeated point elimination to obtain a new vertex set;

(2) and modifying the triangular patch indexes of all models according to the new vertex set to complete the three-dimensional stratum modeling.

Compared with the prior art, the invention has the following technical effects:

1. the invention can be suitable for stratum or geologic body modeling under any complex conditions, and has wider application range compared with the existing modeling method based on two-dimensional contour lines, because:

(1) the existing three-dimensional modeling method based on two-dimensional contour lines mostly adopts a modeling method for connecting the surfaces of the adjacent contour lines in space, and the requirements on the parallel or approximately parallel of the section planes of the graph cut are higher, and the requirements on the manufacture of the section planes of the graph cut are higher. The invention has no special requirement on the manufacture of the graph cut section, can adopt the graph cut section with any direction and any angle for modeling, and can even adopt the graph cut section with non-straight line (broken line).

(2) Because the global interpolation algorithm is adopted to carry out interpolation calculation on the discrete sampling points of the stratum, the invention has no requirement on the space (or density) of the two-dimensional map cutting section. The section can be a sparse contour line as long as the section can control the space form of the three-dimensional geologic body of the model to be built. In addition, the inverse multivariate quadratic radial basis function is used as an interpolation function, and a radial basis function interpolation algorithm is adopted to carry out three-dimensional discrete point spatial interpolation, so that better model smoothness and spatial continuity can be obtained.

2. The modeling efficiency of the invention is higher, which is mainly expressed in that:

(1) the method can be used for modeling by adopting the graph-cut section in any direction, thereby greatly reducing the workload of the graph-cut section and controlling the precision of the model more easily.

(2) The existing contour line of the graph cut section is hardly required to be subjected to additional work of splitting, simplifying, editing, modifying and the like, and only stratum filling and attribute assignment are required.

(3) A fast flooding algorithm based on linear queues is adopted, a kernel function based on threads is created, and stratum boundaries are filled in a multithreading mode, so that the filling efficiency is improved.

3. The modeling precision of the invention is higher, which is mainly expressed in that:

(1) since the filling process is affected by the display resolution of the contour lines, the higher the graphic display resolution, the closer the tracked boundary points are to the true bed boundaries. In the process of filling and stratum boundary detection, the invention performs projection processing on the boundary points and the stratum contour lines, automatically converts the tracked stratum boundary points into contour line boundary points, so as to realize the consistency of the stratum tracking boundary points and the original stratum contour lines, reduce the calculation errors of the stratum boundary points, improve the accuracy of model boundary detection and solve the problem of the dependence of the edge detection on the resolution of display pixels.

(2) The method has the advantages that two sampling modes of boundary sampling and grid sampling are combined, so that the precision of the stratum original model is guaranteed, overlarge sampling data volume is avoided, a three-dimensional global interpolation algorithm can be used for generating a three-dimensional attribute data field, an isosurface is extracted from the stratum three-dimensional attribute data field by combining a Marchang cube algorithm to serve as a three-dimensional surface profile of the stratum, the precision of surface splicing of the follow-up stratum three-dimensional profile is fully guaranteed, and a modeling result is more reasonable and accurate.

In conclusion, the invention improves the stratum boundary tracking algorithm and improves the stratum boundary tracking precision based on the thought of the three-dimensional stratum attribute field; by adopting a special section sampling mode, the accuracy of stratum boundary sampling is guaranteed, overlarge sampling data volume is avoided, and a guarantee is provided for the application of a three-dimensional global interpolation algorithm; by adopting the radial basis function and the GPU-based multithreading calculation technology, the spatial continuity and the smoothness of the stratum model are ensured, and the calculation efficiency is high.

Drawings

FIG. 1 is a schematic flow chart of the present invention.

FIG. 2 is a schematic cross-sectional view of a figure cut in an embodiment of the present invention.

FIG. 3 is a schematic diagram of a cross-sectional view taken along a straight line and a cross-sectional view taken along a broken line in accordance with an embodiment of the present invention.

FIG. 4 is a schematic illustration of the formation fill in an embodiment of the present invention.

FIG. 5 is a diagram illustrating the result of tracking the formation boundary points using a flood fill algorithm according to an embodiment of the present invention.

FIG. 6 is a diagram illustrating the results of projection of formation boundary points and elimination of redundant points according to an embodiment of the present invention.

FIG. 7 is a diagram illustrating the results of stratigraphic boundary tracking, boundary point projection, and redundant point elimination in accordance with an embodiment of the present invention.

FIG. 8 is a schematic diagram illustrating a relationship between plane coordinates and space coordinates of a graph cut section in an embodiment of the present invention.

FIG. 9 is a diagram illustrating the discrete sampling results of the stratigraphic boundary and the grid in accordance with an embodiment of the present invention.

FIG. 10 is a schematic diagram of a discrete sampling result of a three-dimensional space stratum of a graph cut section in an embodiment of the invention.

FIG. 11 is a schematic diagram illustrating a result of a formation property data field generated by three-dimensional interpolation of radial basis functions according to an embodiment of the present invention.

FIG. 12 is a three-dimensional grid diagram of a stratigraphic boundary with an attribute value of 1 according to an embodiment of the present invention.

Fig. 13 is a schematic diagram illustrating a result of calculating a three-dimensional contour surface of a formation by using a Marching Cubes algorithm according to an embodiment of the present invention.

FIG. 14 is a diagram illustrating the results of the three-dimensional geological model after stitching according to an embodiment of the present invention.

Detailed Description

The invention provides a three-dimensional stratum modeling method based on a geological map section, which mainly utilizes the geological map section or other two-dimensional contour lines to carry out three-dimensional geological model modeling, the modeling process is shown as figure 1 and mainly comprises the following steps:

s1, editing the section of the graph to complete the setting of basic parameters required by modeling; the process of editing the graph cut section comprises the following steps: (1) positioning a section space; (2) filling stratum properties; (3) tracing the stratum boundary; wherein: converting the plane coordinates of the section into three-dimensional coordinates of a geological space by calibrating a space positioning point; obtaining formation parameter information through formation attribute filling; obtaining the contour boundary of the stratum through stratum boundary tracking;

The present invention will be further described with reference to the following description and examples, which include but are not limited to the following examples.

Examples

Firstly, editing the section of the graph

FIG. 2 shows a schematic cut-away view supporting either vector graphics format (e.g., DXF) or raster graphics format (JPG, GIF, BMP, etc.). The picture cutting section editing comprises the following steps:

(1) loading a DXF graph of a cutting section, selecting four corner points (at least 2 positioning points are needed) on the graph of the cutting section, setting three-dimensional space coordinates of 4 corner points, and mapping the graph of the cutting section to a space coordinate system. The present invention supports a cut profile of a straight or a folded profile. The straight section means a straight line for forming a drawing section, and the broken section means a broken line for forming a drawing section. The space positioning point corresponding to the folding section is marked on the turning point of the folding line, as shown in fig. 3;

(2) selecting a stratum, filling the stratum by adopting a flooding filling algorithm, recording stratum boundary points of the stratum to obtain a stratum polygon, and setting the name and attribute value (stratum unique identification number) of the stratum, as shown in figure 4;

(3) because the stratum boundary point tracked by the flooding filling algorithm and the actual stratum contour boundary have an error (as shown in fig. 5), the stratum boundary point needs to be subjected to projection processing to obtain an accurate stratum boundary, and the projection processing steps are as follows:

a. and simplifying the stratum contour line, and deleting redundant points on the contour straight line. Then, each two points are formed into a contour line segment { L_kN, where k is 1,2,. N), and N is the number of contour segments;

b. selecting a tracked stratum boundary point P_i(i ═ 1,2, …, M), the set of distances { D) from that point to all contour segments is calculated_ik}；

c. From the set of distances { D_ikFind the line segment L of the minimum distance in_kDividing the boundary point P_iProjected to line segment L_kTo obtain P_ikA 1 is to P_iReplacement by P_ikAs a new groundLayer boundary points;

d. repeating the steps b and c, and projecting all the stratum boundary points onto the stratum contour lines;

e. simplifying the contour line of the boundary of the projected stratum, deleting redundant points on the line, and completing the tracing of the boundary of the stratum (as shown in figures 6 and 7).

Two, discrete sampling and combination into a space discrete point set

(1) Setting an interval step, sampling boundary points along the boundary of the stratum contour line to obtain a stratum boundary point set B1, and setting attribute values of all points in the set to be 1;

(2) setting a subdivision grid M multiplied by N with a specified size (for example, M is 100, N is 80), performing two-dimensional grid subdivision on a map-cutting section, and recording coordinates of a central point of each grid to obtain a subdivision point set B2;

(3) carrying out stratum inside and outside judgment on the points in the set B2, wherein if the points are in a stratum polygon, the attribute value of the points is 1, otherwise, the attribute value is 0;

(4) performing space coordinate conversion on all the points in the sets B1 and B2, and converting the plane coordinate points into geological space coordinate points; the conversion process is as follows (as shown in fig. 8):

wherein, P (x, y) is the coordinate of the plane point to be converted, and P' (x, y, z) is the space coordinate after conversion; p1, P2 is the plane coordinate of the positioning point interval where the point P is located, V1(x, y, z), V2(x, y, z) is the space coordinate corresponding to P1 and P2; y1 is the minimum value of the positioning point in the Y direction of the plane, Y2 is the maximum value of the positioning point in the Y direction of the plane; zmin and Zmax are corresponding minimum elevation and maximum elevation values in space;

the calculation formula for converting the plane point P (x, y) coordinates to the space coordinates P' (x, y, z) is as follows:

p′_x＝V1_x+(p_x-p1_x)*(V2_x-V1_x)/(p_2x-p1_x)

p′_y＝V1_y+(p_x-p1_x)*(V2_y-V1_y)/(p_2x-p1_x)

p′_z＝Z_min+(Z_max-Z_min)*(p_y-Y1)/(Y2-Y1)

the distribution of the sampling points after the sampling of the section of the graph is shown in FIG. 9;

(6) the converted spatial coordinates and their attribute values (0 or 1) are put into the scattered point data set S to obtain the scattered spatial point data after the formation sampling (as shown in fig. 10).

Thirdly, generating a three-dimensional space stratum attribute data field

And performing three-dimensional discrete point spatial interpolation by adopting a radial basis function interpolation algorithm to generate a three-dimensional spatial stratum attribute data field. The radial basis function form is as follows:

wherein Z₀Is the interpolated point value to be estimated, w_iIs the weight of the corresponding point, | | p_i-p₀Is p_iAnd p₀The distance between them. The basis function is an inverse multivariate quadratic function, and the function form is as follows:

the interpolation calculation process is as follows:

(1) calculate distance matrix { D_ij}(i＝1,2,…,N,j＝1,2,…,N)，

(3) A matrix W of weight coefficients is calculated,

Z_jis the jth discrete point P_jA formation property value of;

The three-dimensional attribute data field after the three-dimensional interpolation is completed is shown in fig. 11, in which the grid distribution with attribute value of 1 is shown in fig. 12.

Fourthly, obtaining the three-dimensional contour surface of the stratum

And extracting an isosurface with the attribute value of 1 from the stratum three-dimensional attribute data field by adopting a Marching Cubes algorithm to be used as a three-dimensional surface profile of the stratum, and completing the three-dimensional modeling of the stratum, as shown in FIG. 13.

The process of extracting the isosurface by the Marching Cubes algorithm is as follows:

p(x，y，z)＝p_1(x，y，z)+(v-v1)/(v2-v1)*(p_2(x，y，z)-p_1(x，y，z))

in the formula, p (x, y, z) is the coordinate of the intersection point to be calculated, and v is the isosurface value to be calculated; p1, p2 are two-node coordinates of an edge, v1 and v2 are two-node attribute values (the coordinates meet a vector algorithm); in the embodiment, the HashMap linear table is adopted to store the edge intersections, so that the intersections of the isosurface and each edge are not repeatedly calculated and stored;

And fifthly, repeating the process to create three-dimensional surface models of all the stratums.

Sixthly, completing three-dimensional stratum modeling

Splicing all stratum models spatially, wherein the splicing process is as follows:

(2) the triangular patch indices of all models are modified according to the new set of vertices.

Thus, a three-dimensional model is generated and output in a three-dimensional graphic format (such as obj or ply), thereby completing the three-dimensional formation modeling, as shown in fig. 14.

The method has reasonable and rigorous design, avoids the problems of complex contour correspondence, triangulation network connection and the like, and can obtain a high-precision three-dimensional vector geological model. The method well breaks through the limitation of the prior art, conforms to the trend of scientific and technological development, and well matches the requirements of the current stage engineering on three-dimensional geological modeling. Therefore, compared with the prior art, the invention has obvious technical progress, and has outstanding substantive characteristics and remarkable progress.

The above-mentioned embodiment is only one of the preferred embodiments of the present invention, and should not be used to limit the scope of the present invention, and all the technical problems solved by the present invention should be consistent with the present invention, if they are not substantially modified or retouched in the spirit and concept of the present invention.

Claims

1. A three-dimensional stratum modeling method based on geological map cutting profiles is characterized by comprising the following steps:

2. The method of claim 1, wherein in step S1, the stratum is filled by a fast flooding filling algorithm based on linear queues, and boundary points of the stratum are recorded to obtain a stratum polygon, and stratum parameter information is obtained.

3. The method as claimed in claim 1, wherein in step S1, projection processing is used to obtain the boundary points of the formation by tracking, so as to eliminate the tracking error of the boundary points due to the influence of the display resolution of the raster image and ensure the consistency between the boundary of the contour of the tracked formation and the contour of the original formation.

4. The method of claim 3, wherein the projection process comprises the steps of:

(1) simplifying the stratum contour line, deleting redundant points on the contour straight line, and then forming a contour by every two pointsLine segment { L_kN, where k is 1,2,. N), and N is the number of contour segments;

5. The geological map cut-section-based three-dimensional stratum modeling method according to any of claims 1-4, wherein the step S2 comprises the following steps:

s2-6, putting the converted spatial coordinate points and the attribute values thereof into a scattered point data set { p_iAnd (i is 1,2, …, N), obtaining discrete spatial point data after stratum sampling, and completing sampling.

6. The method according to claim 1, wherein in step S3, an inverse multivariate quadratic radial basis function is used as an interpolation function, and a radial basis function interpolation algorithm is used to perform three-dimensional discrete point spatial interpolation to generate a three-dimensional spatial formation property data field; the radial basis function form is as follows:

the interpolation calculation process is as follows:

(1) calculate distance matrix { D_ij}(i＝1,2,…,N,j＝1,2,…,N)，

(3) A matrix W of weight coefficients is calculated,

Z_jis the jth discrete point P_jA formation property value of;

7. The method for modeling three-dimensional stratum based on geological map cutting section as claimed in claim 1, wherein in step S4, Marching Cubes algorithm is used to extract isosurface from stratum three-dimensional attribute data field as three-dimensional surface contour of stratum, the process is as follows:

p(x，y，z)＝p_1(x，y，z)+(v-v1)/(v2-v1)*(p_2(x，y，z)-p_1(x，y，z))

8. The method for modeling three-dimensional stratum based on geological map cutting section according to claim 1, characterized in that said step S6 includes the following steps: