Structural plane three-dimensional modeling method based on trend control and occurrence constraint

Abstract

The invention belongs to the technical field of three-dimensional geologic modeling, and particularly relates to a three-dimensional modeling method for a structural plane based on strike control and occurrence constraint, which comprises the following steps of: extracting strata, fault strike lines and occurrence data based on a plane geologic map of a modeling target area, and extracting coordinates and azimuth angles of points on the strike lines through a curvature-constrained equal-interval sampling method; generating layer control points in batches by changing the sampling depth in combination with occurrence data; expanding a layer control point set by utilizing a stratum or fault layering point disclosed by a drill hole or a geological section, and further creating an irregular triangulation network of a structural plane; and extracting a structural surface intersection line by adopting a collision detection algorithm, reconstructing a stratum contact relation and a fault network space topological relation, finally obtaining a layer data file of the structural model, and completing data storage and three-dimensional visualization of the model. According to the method, the model of the structural plane such as the stratum and the fault can be quickly established, the spatial topological relation is reconstructed, and data file generation and three-dimensional visualization of the structural model are realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of three-dimensional geological modeling, and particularly relates to a three-dimensional modeling method for structural surfaces based on strike control and attitude constraints. Background Technique

[0002] The geological guarantee system is an information-based support technology for promoting the safe production of coal mines and is also an important part of the construction of intelligent mines. Three-dimensional geological modeling is a key technology in the construction of the mine geological guarantee system and is also the carrier and data foundation for the construction of transparent geology, digital twin, and intelligent management and control platforms under the intelligent mine strategy. Three-dimensional geological modeling uses computer technology to express the spatial form and internal attribute field of geological bodies in three-dimensional space. It can not only display the spatial structure of geological bodies but also use tools such as spatial information analysis technology and geostatistics to achieve the integration of multi-source heterogeneous data, process simulation, and analysis and prediction. It has been widely used in scenarios such as numerical simulation and disaster prevention.

[0003] Structural modeling mainly focuses on the creation of structural surfaces such as fault planes, stratigraphic planes, and unconformity surfaces. The modeling process takes the creation of three-dimensional surfaces as the core. Using known data to construct a structural model that conforms to actual geological features such as attitude descriptions, contact relationships, and fault networks is the most important link in three-dimensional geological modeling. Generally, it can be achieved through explicit modeling, implicit modeling, or a combination of both. Borehole data is the measured result of the stratigraphic structure and is one of the most reliable data sources for three-dimensional geological modeling. Establishing a layer triangulation network using borehole data is a common method for three-dimensional geological modeling. Boreholes with relatively uniform distribution can bring better structural surface generation effects and are therefore widely used in various three-dimensional geological modeling problems. Three-dimensional seismic can accurately reflect the spatial distribution of underground deep structural surfaces, and its interpretation results are also important data sources for three-dimensional geological modeling. However, due to its high upfront investment cost, it is mostly applied to models with complex structures, large burial depths, and high precision requirements. In addition, cross-sections can intuitively reflect the stratigraphic structure and its topological relationships and are also often used as important data sources for structural modeling. Using a certain number of cross-sections can also obtain a good three-dimensional geological model. Obviously, in structural modeling, borehole data, three-dimensional seismic interpretation results, geological cross-sections, etc. all play very important roles. The comprehensive use of one or more types of data with sufficient quantity and quality will determine the accuracy and reliability of structural modeling. Existing modeling methods and software are designed and developed based on this and can meet most structural modeling scenarios.

[0004] However, in actual work, one or more of the above-mentioned basic data may be lacking, or the distribution characteristics of the data are not ideal, making it difficult to meet the requirements of structural modeling, resulting in the difficulty in carrying out existing structural modeling techniques and even being unable to accurately describe the variation characteristics of the geological body space. Specifically, first, sparse borehole data often cannot effectively constrain the interpolation algorithm during implicit modeling, resulting in the easy drift and distortion of the structural surface and increased uncertainty. Second, the lack of sufficient profile data makes it difficult to generate the correct spatial topological relationship for the stratigraphic contact relationship and fault network. At the same time, if there is no corresponding 3D seismic interpretation data in the modeling target area, conventional modeling methods and modeling software cannot play a role due to the lack of necessary basic data. Moreover, when modeling special geological structures, the lack of data will further exacerbate the uncertainty of modeling, and existing modeling methods cannot even meet the requirements. For example, steeply inclined strata widely exist in the two wings of folds with strong tectonic compression. Among them, for example, steeply inclined coal seams are the main coal seams to be mined. However, due to the too large dip angle of the coal seams, it not only hinders the efficient mining of coal mines, but also because the exposure of boreholes to steeply inclined strata is extremely limited, conventional structural modeling techniques also face new challenges. In this case, generally, the original vertical boreholes can be converted into horizontal virtual boreholes by means of virtual boreholes, and structural modeling can be carried out by means of software deception. However, this method has high requirements for data preprocessing. When the geological structure is complex or the number of boreholes is large, a large amount of data collation is carried out manually, with a high error rate and low efficiency. And when the number of boreholes is limited, it is difficult to play an effective role, with strong limitations.

[0005] The planar geological map is a low-cost data source that is easy to obtain and directly reflects the regional geological characteristics. It integrates expert knowledge and experience, can not only provide a regional macro geological background for geological modeling, but also make up for the lack of fine modeling data such as boreholes, and is also widely used in 3D geological modeling. In the existing technology, a 3D modeling method for complex faults based on the planar geological map has been studied. In the existing technology, a research idea of comprehensively using low-cost data such as profile maps and geological maps for 3D geological modeling has also been proposed. Gao Shujuan, Zhou Liangchen and others have studied a 3D geological modeling method based on the planar geological map and framed by cross-section maps, and pointed out that using the planar geological map for 3D geological modeling is an effective solution when other geological data are lacking.

[0006] In view of this, the present invention provides a three-dimensional modeling method for structural planes based on strike control and attitude constraint. In actual work, a solution for structural modeling is proposed through analysis and summary, which is based on low-cost plane geological maps or structural outline maps to obtain basic data. By using the layer control point generation algorithm based on strike control and attitude constraint and the spatial topological relationship reconstruction algorithm, computer programs are used to batch generate structural plane control point data, quickly establish structural plane models such as strata and faults, and reconstruct the spatial topological relationship, so as to realize the generation of data files of the structural model and three-dimensional visualization. The feasibility of the method is verified through application examples, providing a solution for three-dimensional geological modeling in sparse data areas. Summary of the Invention

[0007] The object of the present invention is to provide a three-dimensional modeling method for structural planes based on strike control and attitude constraint. By using the layer control point generation algorithm based on strike control and attitude constraint and the spatial topological relationship reconstruction algorithm, computer programs are used to batch generate structural plane control point data, quickly establish structural plane models such as strata and faults, and reconstruct the spatial topological relationship, so as to realize the generation of data files of the structural model and three-dimensional visualization. The feasibility of the method is verified through application examples, providing a solution for three-dimensional geological modeling in sparse data areas.

[0008] The technical solutions adopted by the present invention are specifically as follows:

[0009] A three-dimensional modeling method for structural planes based on strike control and attitude constraint, comprising the following steps:

[0010] Step 1: Extract the strike lines and attitude data of strata and faults from the plane geological map of the modeling target area. Through the equidistant sampling method with curvature constraint, extract the coordinates and azimuth angles of points on the strike line; combined with the attitude data, batch generate layer control points by changing the sampling depth; then use the stratification points of strata or faults exposed by limited boreholes or geological profiles to expand the layer control point set, and then create an irregular triangular network of the structural plane;

[0011] Step 101: Use GIS software to extract the strike line of the layer from the plane geological map of the study area. The layer is the stratigraphic interface or fault plane. Perform equidistant sampling on the strike line and extract the coordinate set {P i} = {X i , Y i , Z0, θ i}, where i = 0, 1,..., n; Z0 represents the initial elevation value or depth value of the sampling of the layer object, which can be obtained by extracting the elevation of the strike line and the DEM. The corresponding initial depth can be set for modeling in the depth domain. θ i represents the azimuth angle of the strike line at this point, and organize to obtain the data table of the strike line;

[0012] Meanwhile, curvature constraints are used to control the morphology of sampling points and thin out the data. The principle is as follows: The B-spline curve equation passing through the strike sampling points can be quickly obtained through data fitting, and let it be y = f(x); taking one end of the strike line as the starting point, every three adjacent points are used to calculate the curvature, and traverse to the other end point of the curve; P i (x i ,y i ) The curvature at the point can be calculated by formula (1);

[0013]

[0014] Thus, the curvature set {κ1, κ2, …, κ n-2} of the strike line is obtained;

[0015] Step 102: Taking the strike line as the top boundary line of the layer, a right triangle is constructed in three-dimensional space for the strike line sampling points according to the layer dip and dip angle. Let the layer dip angle be α, and the projection point of the sampling point P i at the depth h is denoted as O. According to the trigonometric function relationship of plane projection, the coordinates (X′ i , Y′ i ) of P′ i at the depth h are related to the layer dip and strike azimuth angle; when the layer dip is southward, the coordinates (X′ i , Y′ i ) of P′ i can be calculated by formulas (2) and (3) respectively. When the formation dip is northward, the coordinates (X′ i , Y′ i ) of P′ i are calculated by formulas (4) and (5) respectively;

[0016]

[0017] Based on this, the control point set {P′ i} of the formation surface at the depth h is obtained = {X′ i , Y′ i , Z h}, i = 0, 1, …, n, Z h represents the elevation corresponding to the depth h at P′ i . If depth-domain modeling is adopted, it can be directly converted to the depth Z h = h; by controlling the strike line sampling interval and setting different sampling depths, the layer control point set {P i , P′ i , P″ i , …} can be obtained; based on this, with the constraints of the layer occurrence and sampling depth, a large number of layer control points can be generated quickly in batches; at the same time, the sampling depth The positive and negative value changes can flexibly adjust the sampling range of the layer control points;

[0018] Step 103: Insertion and correction of borehole stratification points; Extract the layer stratification points from the borehole and insert them into the control point set {P i , P′ i , P″ i , …} created in Step 102, and encrypt and correct the layer control point set.

[0019] Step 104: Store the layer control point data, where the Type field records the layer type represented by the control point set {P i , P′ i , P″ i , …}′; Horizon represents the ground layer, Fault represents the fault, the Name field records the layer name, and the PointID field records the index number of the control point; within the same layer, PointID is a unique value, and X, Y, and Z respectively record the coordinate values of the control point.

[0020] Step 105: Creation of layer triangulation network; After obtaining the layer control point set, directly perform Delaunay triangulation to generate the TIN of the layer; Save the vertex coordinates and triangle connection order in the TIN in text format.

[0021] Step 2: Use the collision detection algorithm to extract the intersection lines of the structural surfaces, reconstruct the stratigraphic contact relationship and the spatial topological relationship of the fault network, and finally obtain the layer data file of the structural model and complete the data storage and three-dimensional visualization of the model.

[0022] After establishing the layer TIN according to Step 1 in the said Step 2, construct the stratigraphic contact relationship and the fault network topological relationship. The method principle is as follows:

[0023] Step 201: Use the collision detection algorithm to screen the intersecting triangles. Use the method of axis-aligned bounding box collision detection to extract the intersection lines of the TINs of two intersecting layers; Enclose the triangle completely with a cube with edges parallel to the coordinate axes, and use a recursive method to form a multi-way tree structure to implement the collision detection algorithm;

[0024] Step 202: Extraction of intersection surface lines: Since the layers are all constructed using triangulation networks, after screening by the axis-aligned bounding box collision detection method, the extraction of the intersection lines is transformed into a triangle intersection operation; When the bounding boxes of two triangles overlap, it is considered that the triangles intersect, and the sides of one triangle are respectively judged with the other triangle to obtain the intersection points; Use the Möller–Trumbore algorithm to calculate the intersection lines. The principle is as follows:

[0025] Suppose the plane is represented by the normal vector N and a point p' on the plane of the triangle. Then, for any point p on the plane, the following formula (6) is satisfied:

[0026] p: (p - p')·N = 0 (6)

[0027] A ray is represented by its starting point O and direction vector Let t represent the ray length ratio. Then, the ray r satisfies the following formula (7):

[0028]

[0029] Let the three vertices of the triangle be P0, P1, and P2. According to the rules of barycentric coordinate interpolation, the intersection of the ray and any point on the plane can be expressed by the following formula (8):

[0030]

[0031] In the formula, b1, b2, and 1 - b1 - b2 represent the barycentric coordinates of any point on the plane of the triangle; P0, P1, and P2 represent the coordinates of the three vertices of the triangle. There are only three unknowns, t, b1, and b2, in formula (8). When t, b1, b2, and 1 - b1 - b2 are all greater than or equal to 0, the ray intersects with the triangle.

[0032] Accordingly, traverse the vertices of the two triangles respectively. A total of six calculation operations are required for each pair of triangles. When t ∈ [0, 1], it indicates that there is an intersection between the side of the current triangle and the other triangle. Record the obtained intersection points and their adjacent intersection points. If there is no adjacent point, it is the end point of the intersection line. Form an intersection point linked list and trace all the intersection points in sequence to obtain the intersection lines of the two levels. All the intersection points are marked with a unique PointID. While saving the point coordinates, set the Neighbors field to save the adjacent points before and after the current point. The end point has only one adjacent point. According to the tracing order of the intersection points, the first recorded is the starting point, and the last recorded is the end point of the intersection line.

[0033] Step 203: Spatial topology reconstruction and layer clipping; Based on the TIN data created above, including the layer control point coordinates and triangle vertex indices, construct and clip the spatial topology relationships such as the formation contact relationships and fault networks, which are transformed into the problem of dividing the layer control point set by the intersection lines and reconstructing the triangulation network after division.

[0034] First, set the master-slave relationship between the two surfaces according to the stratigraphic column or fault cutting relationship, and then generate the bounding box using the intersection point set. At this time, for the slave surface; when reconstructing the stratigraphic contact relationship, if the strata are in an erosional relationship, the master surface is the unconformity surface and the slave surface is the older strata; when constructing the fault network, it is determined by the sequence or importance of fault development; when the fault cuts the strata, the control point set with the fault plane as the master and the strata as the slave will be divided into three parts: inside the bounding box and on both sides of the bounding box; subsequently, according to the coordinate range, the slave surface control points outside the bounding box are divided into two groups, A and B, and any point is selected from them; assume a point P is selected from group A, and a point pair is formed with the slave surface control point Q inside the intersection line bounding box. Again, use the Möller–Trumbore algorithm to determine whether the line segment PQ intersects the master surface. If the intersection point exists, Q should be classified into group B, otherwise it belongs to group A; after the classification is completed, according to the erosional relationship, retain the data points on one side of the slave surface points to complete the surface clipping;

[0035] Using the intersection line detection and surface clipping algorithm, pairwise judgment is performed by traversing the TIN file. While extracting the intersection line, the topological reconstruction is completed.

[0036] Step 204: TIN repair; Add the intersection line points to the clipped slave surface point set, and use the Delaunay triangulation algorithm to locally repair the TIN and update the surface data file.

[0037] The technical effects achieved by the present invention are:

[0038] The present invention utilizes strike control and attitude constraints, and proposes an algorithm for batch generation of structural surface control points, and on this basis, forms an algorithm flow for three-dimensional surface creation. This method can make full use of the limited geological data in the modeling target area, quickly generate surface control point data, realize structural modeling, and solve the pain point of being unable to model due to insufficient data.

[0039] Based on mature computer graphics basic models such as collision detection, the present invention realizes the algorithm for reconstructing the spatial topological relationship of structural surfaces, and forms the entire process from data generation, storage to three-dimensional visualization of the model. The method flow is simple, easy to program and implement, and considers the compatibility and expansion of data formats, which is convenient for data sharing.

[0040] The application of the present invention example shows that relying on the attitude information and structural features extracted from the plane geological map, the method of the present invention effectively solves the problems such as the special structure of steeply dipping strata in the target area, the lack of basic data such as boreholes, and the inability to realize structural modeling by conventional methods, verifying the feasibility and reliability of the method. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 is the flow chart of the method of the present invention;

[0042] Figure 2 It is a schematic diagram of curvature calculation of the present invention;

[0043] Figure 3 It is a schematic diagram of calculating layer sampling points under the occurrence constraint of the present invention;

[0044] Figure 4 It is a schematic diagram of the principle of calculating sampling points of the present invention;

[0045] Figure 5 It is the encryption of borehole layer data and the correction of layer control points of the present invention;

[0046] Figure 6 It is a schematic diagram of topological clipping of the present invention;

[0047] Figure 7 It is a simplified structure diagram of Well Wei-2 Minefield and the acquisition results of main fault plane data in the first example of the present invention;

[0048] Figure 8 It is a 3D structure model of the first mining area of Well Wei-2 Minefield in the first example of the present invention;

[0049] Figure 9 It is the occurrence of bedrock strata in the modeling target area in the second example of the present invention;

[0050] Figure 10 It is a geological section (partial) of the modeling target area in the second example of the present invention;

[0051] Figure 11 It is the modeling process and results in the steeply inclined strata area in the second example of the present invention. (a) is the extraction result of sampling points of the strike line in the modeling target area; (b) is the sampling result of the occurrence-constrained layer control points proposed in this paper; (c) is the steeply inclined strata and fault structural planes rendered based on the layer control points; (d) is the extraction of the intersection line between two layers; (e) and (f) are the tracking results of the layer intersection line and the structural model after topological reconstruction and layer clipping respectively. Specific Embodiments

[0052] In order to make the purpose and advantages of the present invention clearer, the present invention will be specifically described below in conjunction with embodiments. It should be understood that the following text is only used to describe one or several specific implementation manners of the present invention, and does not strictly limit the scope of protection of the specific claims of the present invention.

[0053] Embodiment 1:

[0054] A 3D modeling method for structural planes based on strike control and occurrence constraint, comprising the following steps:

[0055] Step 1: Extract formation, fault strike lines and attitude data from the planar geological map of the modeling target area. By using the equally spaced sampling method with curvature constraint, extract the coordinates and azimuth angles of the points on the strike lines. Combine the attitude data and batch generate layer control points by varying the sampling depth. Then, use the formation or fault bedding points exposed by limited boreholes or geological profiles to expand the layer control point set, and further create an irregular triangular network of the structural surface.

[0056] Step 2: Based on the ideas of computational geometry and computer graphics, adopt a collision detection algorithm to extract the intersection lines of the structural surfaces, and reconstruct the spatial topological relationships such as formation contact relationships and fault networks. Finally, obtain the layer data file of the structural model, and complete the data storage and three-dimensional visualization of the model.

[0057] In the present invention, when constructing a three-dimensional surface in a computer, a triangle is the basic unit. By connecting 3 adjacent points to form a triangle and then connecting the triangles to form an irregular surface, the creation of a Triangulated Irregular Network (TIN) is an essential step in the modeling of the structural surface. Geological modeling based on TIN can adjust the size and quantity of triangles according to the complexity of the modeling target. While retaining the topological relationship, it is easy to handle complex geological structures. Based on the core idea of constructing the TIN of the structural surface, this invention uses the planar geological map as the basic data, proposes an automatic batch generation algorithm for the control points of the formation surface and fault surface aiming at the problem of insufficient basic modeling data, and designs a reconstruction algorithm for the topological relationship of the structural surface on this basis, forming a method flow for generating and modeling the structural surface data in the case of lack or unsatisfactory data of the modeling basic data, and considering the generality of the data format and the storage and interaction problems. The method flow of this invention is as Figure 1 shown. The core work is as follows: First, extract formation, fault strike lines and attitude data from the planar geological map of the modeling target area. By using the equally spaced sampling method with curvature constraint, extract the coordinates and azimuth angles of the points on the strike lines. Combine the attitude data and batch generate layer control points by varying the sampling depth. Then, use the formation or fault bedding points exposed by limited boreholes or geological profiles to expand the layer control point set, and further create an irregular triangular network of the structural surface. Then, based on the ideas of computational geometry and computer graphics, adopt a collision detection algorithm to extract the intersection lines of the structural surfaces, and reconstruct the spatial topological relationships such as formation contact relationships and fault networks. Finally, obtain the layer data file of the structural model, and complete the data storage and three-dimensional visualization of the model.

[0058] Preferably, in view of the fact that existing geological modeling software is difficult to accurately and effectively construct a three-dimensional model in the case of lack of necessary modeling data such as boreholes, this paper proposes a three-dimensional modeling algorithm for the structural surface based on strike control and attitude constraint. The specific steps of step 1 are as follows:

[0059] Step 101: Extract the strike lines of the horizons from the planar geological map of the study area with the aid of GIS software. The horizons are stratigraphic interfaces or fault planes. Sample the strike lines at equal intervals and extract the coordinate set {P i} = {X i , Y i , Z0, θ i}, where i = 0, 1, …, n; Z0 represents the initial elevation value or depth value of the sampling of the horizon object, which can be obtained by extracting the elevation from the strike line and the DEM. The corresponding initial depth can be set in the depth domain modeling. θ i represents the azimuth angle of the strike line at this point. Organize to obtain the data table of the strike line;

[0060] Meanwhile, to avoid sharp changes in the strike during sampling or computational redundancy caused by overly dense data points, curvature constraints are used to control the shape of the sampling points and thin the data. The principle is as follows: The B-spline curve equation passing through the strike sampling points can be quickly obtained by data fitting. Let it be y = f(x); Take one end of the strike line as the starting point, and use every three adjacent points to calculate the curvature, traversing to the other end point of the curve; For Figure 2 the P i (x i , y i ) point shown, the curvature can be calculated by formula (1);

[0061]

[0062] Thus, the curvature set {κ1, κ2, …, κ n-2} of the strike line is obtained; Set a reasonable threshold for r in formula (1), then the change trend of the strike line can be controlled by the curvature, achieving the avoidance of sharp changes in the formation strike. At the same time, the data points can be thinned in the area where the strike changes little to reduce the computational amount in the subsequent steps;

[0063] Step 102: Use the strike line as the top boundary line of the horizon, and construct right-angled triangles in three-dimensional space for the strike line sampling points according to the dip and dip angle of the horizon; As Figure 3 shown; Let the dip angle of the horizon be α, and the projection point of the sampling point P i at depth h is denoted as O. According to the trigonometric function relationship of plane projection, the coordinates (X′ i , Y′ i ) of P′ i at depth h are related to the dip and strike azimuth angle of the horizon; When the dip of the horizon is southward, as Figure 4 -a and Figure 4 -b shown, the coordinates (X′ i , Y′ i , Y′ i) can be calculated respectively using formulas (2) and (3). When the formation dip is northward, as shown in Figure 4 -c and Figure 4 -d, the coordinates (X′ i , Y′ i ) of P′ i are calculated respectively using formulas (4) and (5);

[0064]

[0065] Based on this, the control point set of the formation surface at depth h is obtained i = 0, 1, …, n, Z h represents the elevation corresponding to the depth h at P′ i . If depth-domain modeling is adopted, it can be directly converted to depth Z h = h; By controlling the sampling interval of the strike line and setting different sampling depths, the control point set of the formation surface {P i , P′ i , P″ i , …} can be obtained; Based on this, with the constraints of the formation attitude and sampling depth, a large number of control points of the formation surface can be generated quickly in batches; At the same time, the positive and negative value changes of the sampling depth can flexibly adjust the sampling range of the control points of the formation surface;

[0066] Step 103: Insertion and correction of borehole stratification points; The stratification points of the formation made based on the borehole columnar diagram are control point data with relatively higher accuracy. Even with a small number of borehole stratification data, they can still correct the formation surface; Therefore, extract the formation stratification points from the boreholes and insert them into the control point set {P i , P′ i , P″ i , …} created in step 102 to encrypt and correct the control point set of the formation surface.

[0067] Step 104: Store the control point data of the formation surface, where the Type field records the formation surface type represented by the control point set {P i , P′ i , P″ i , …}′; Horizon represents the formation surface, Fault represents the fault, the Name field records the formation surface name, and the PointID field records the index number of the control point; In the same formation surface, PointID is a unique value, and X, Y, and Z record the coordinate values of the control point respectively; In practical applications, each formation surface can be saved separately in a file; The data file storage format shown in Table 1 can be directly used in general 3D geological modeling software; In addition, with the supplement of borehole stratification data or seismic interpretation data, organizing the data according to the format of Table 1 can quickly realize the expansion and update of the control points of the formation surface.

[0068] Table 1 Data Structure of Layer Control Point Set

[0069]

[0070] Step 105: Creation of layer triangulation network; After obtaining the layer control point set, Delaunay triangulation can be directly performed to generate the TIN of the layer; The vertex coordinates and triangle connection order in the TIN are saved in text format, and the file format is shown in Table 2; According to the layer point data saved in Table 2, it can be directly rendered into a three-dimensional layer, or imported as point cloud data into other three-dimensional geological modeling software to achieve data interaction.

[0071] Table 2 Saving Format of Layer TIN Data

[0072]

[0073] Preferably, in structural modeling, the determination of the intersection line of layers is crucial. Taking the intersection line of strata and faults as an example, the general method is to encrypt the strata data, and obtain the intersection line according to the influence range and throw of the fault, which is called the overall interpolation method of the fault. It requires manual intervention to specify the influence range of the layer, and the algorithm is complex and has poor adaptability. After establishing the layer TIN in step 2 according to step 1, the contact relationship of the strata and the topological relationship of the fault network are constructed, and the method principle is as follows:

[0074] Step 201: Use the collision detection algorithm to screen intersecting triangles. Using the method of axis-aligned bounding box (AABB) collision detection, the intersection line of the TINs of two intersecting layers is extracted; The axis-aligned bounding box is based on the ideas of computational geometry and computer graphics. Taking the triangle as the smallest unit, a cube with edges parallel to the coordinate axes is constructed to completely enclose the triangle, and a multi-tree structure is formed by recursion to implement the collision detection algorithm. It only needs to judge whether the projections of the two bounding boxes on the coordinate axes coincide to judge whether the triangles intersect; This method has a simple structure, less storage space occupancy, and low algorithm complexity, and can quickly screen the intersecting triangles.

[0075] Step 202: Extraction of the intersection surface line: Since the layers are all constructed by triangulation networks, after being screened by the axis-aligned bounding box collision detection method, the extraction of the intersection line is transformed into a triangle intersection operation; When the bounding boxes of two triangles overlap, it is regarded as the triangles intersecting, and the sides of one triangle are respectively judged with the other triangle to obtain the intersection points; The Möller–Trumbore algorithm is used to calculate the intersection line, and the principle is as follows:

[0076] Let the plane be represented by the normal vector N and a point p' on the plane of the triangle, then any point p on the plane satisfies formula (6):

[0077] p: (p - p')·N = 0(6)

[0078] Using the starting point O and the direction vector to represent the ray, and t represents the ray length ratio, then the ray r satisfies formula (7):

[0079]

[0080] Let the three vertices of the triangle be P0, P1, P2. According to the rules of barycentric coordinate interpolation, the intersection of the ray and any point on the plane can be expressed as formula (8):

[0081]

[0082] In the formula, b1, b2, 1 - b1 - b2 represent the barycentric coordinates of any point on the triangle plane; P0, P1, P2 represent the coordinates of the three vertices of the triangle; there are only three unknowns t, b1, b2 in formula (8). When t, b1, b2, 1 - b1 - b2 are all greater than or equal to 0, the ray intersects with the triangle;

[0083] Accordingly, traverse the vertices of the two triangles respectively. A total of six calculation operations are required for each pair of triangles. When t ∈ [0, 1], it means that there is an intersection between the side of the current triangle and the other triangle; record the obtained intersection points and their adjacent intersection points; if there is no adjacent point, it is the end point of the intersection line; form an intersection point linked list, and trace all the intersection points in sequence to obtain the intersection lines of the two layers; the storage format of the intersection point data is shown in Table 3. All intersection points are marked by a unique PointID. While saving the point coordinates, set the Neighbors field to save the front and back adjacent points of the current point. The end point has only one adjacent point; according to the tracing order of the intersection points, the first recorded is the starting point, and the last recorded is the end point of the intersection line.

[0084] Table 3 Storage format of intersection line points

[0085]

[0086] Step 203: Spatial topology reconstruction and layer clipping; based on the TIN data created above; including the construction and layer clipping of spatial topology relationships such as layer control point coordinates, triangle vertex indices, stratigraphic contact relationships, and fault networks, which are converted into the problem of dividing the layer control point set by the intersection line and reconstructing the triangulation network after division;

[0087] First, the master-slave relationship of the two levels is set according to the stratigraphic column or fault cutting relationship, and then the bounding box is generated using the intersection point set. At this time, the slave surface; when reconstructing the stratigraphic contact relationship, if the strata are in an erosion relationship, the master surface is the unconformity surface and the slave surface is the older strata; when constructing the fault network, it is determined by the order or priority of the fault development; when the fault cuts the strata, the fault surface is the master and the stratigraphic surface is the slave) The control point set will be divided into three parts: the inside of the bounding box and the sides of the bounding box; then, according to the coordinate range, the slave surface control points outside the bounding box are divided into two groups AB, and any point is selected from them; assuming that a point P is selected from group A, such as Figure 6 As shown; the control point Q of the slave surface inside the intersection line bounding box forms a point pair, and the Mohler-Tremblay algorithm is used again to determine whether the line PQ has an intersection with the main surface. If the intersection exists, it should be classified into group B with point Q, otherwise it should be classified into group A; after the classification is completed, according to the erosion relationship, the data points on one side of the slave surface point are retained to complete the surface clipping;

[0088] By using the intersection detection and layer clipping algorithm, the TIN file is traversed to make pairwise judgments, and the topology reconstruction is completed while the intersection extraction is completed.

[0089] Step 204: TIN repair; adding the intersection points to the clipped secondary surface point set, using the Delaunay triangulation algorithm to perform local repair on the TIN, and updating the layer data file.

[0090] Example 1:

[0091] The structural morphology of the Weizhou mining area in the Ningdong coalfield is mainly controlled by the Weizhou syncline, and the secondary structure is mainly faults. The Weier mining area is located in the central and southern part of the Weizhou mining area, that is, the southern part of the eastern wing of the Weizhou syncline. It is a monocline structure with a stratum dip of 5 to 35 degrees to the west. The stratum gradually slows down from east to west. There are two groups of large-scale and large-scale oblique faults in the mining area, the NW (mainly compressional reverse faults) and the NNE (mainly tensional normal faults). There are more than 40 proven faults in the mining area, including 11 large faults with long extension lengths. Table 4 lists the occurrence of these faults.

[0092] Table 4 Fault occurrence in Weier Mine Field

[0093]

[0094]

[0095] Figure 7 (a) shows the main structure and coal seam boundary of Weier Mine Field. Figure 7 (b) Using the proposed construction surface control point batch generation algorithm, based on Figure 7 The fault strike line of (a) and the fault strike data of Table 4 generate the main fault plane control points of the well field. Figure 7(c) Fault structure planes generated for the control points at the usage level. It is worth mentioning that during the batch generation of fault plane control points, considering the characteristics of the minefield elevation and calculating the layer points sampling based on the strike line, the sampling depth h has corresponding changes. Therefore, in Figure 7 (a) and Figure 7 (b), the strike line is shown as the initial sampling position near the middle of the layer.

[0096] In the Wei-2 minefield, the coal seams are relatively developed. The coal-bearing strata are the Lower Permian Shanxi Formation (P1s) and the Carboniferous - Permian Taiyuan Formation (C2-P1t). Among them, the main exploitable coal seams in the Taiyuan Formation are mainly concentrated in the second and third sections of the Taiyuan Formation, with a total of 6 main exploitable seams (seams 12, 14, 15, 16, 17, 20), and 3 exploitable coal seams in the Shanxi Formation (seams 2, 3, 4). Using the method proposed in this paper and combining the coalfield geological borehole data, the structural plane data of the exploitable coal seams are generated by the same principle.

[0097] To serve the selection of well locations for coalbed methane wells in the coordinated development of coal and coalbed methane in the minefield, a 3D structural modeling is carried out for the early mining area of the Wei-2 minefield. After obtaining the above data, a 3D geological model of the early mining area of the Wei-2 minefield is further constructed to display the spatial distribution characteristics of the geological structure in the mining area. Figure 8 (a) shows the boundary of the early mining area of the minefield and the distribution of boreholes, Figure 8 (b) shows the structural model of the early mining area, where the spatial position relationship between the Shanxi Formation and Taiyuan Formation coal seams and the main faults in the minefield is clearly and intuitively displayed.

[0098] Example 2:

[0099] The structure of another modeling target area in the actual work of the present invention is relatively special. The bedrock strata in the area are strongly compressed by the structure, showing steeply inclined monoclinic strata, which are unconformably in contact with the overlying loose strata. The plane geological sketch of the target area is as Figure 9 shown. The bedrock strata and faults are all in the northeast strike. The cross-section ([[]]END]] Figure 10 ) shows that the dip angle of the bedrock strata > 80 degrees, showing a nearly vertical state.

[0100] Due to the steeply inclined characteristics of the strata, the exposure of the bedrock in the borehole data is extremely limited. Generally, there are very few multiple strata stratification points in the same borehole. When using conventional modeling software for modeling, it is difficult to form effective constraints. In implicit modeling, there are strata distortions and holes, and the unconformable contact between the bedrock strata and the overlying sedimentary strata cannot generate the correct spatial topological relationship. In addition, since there is only one available geological cross-section in the modeling target area, explicit modeling is also difficult to play a role.

[0101] Using the algorithm for generating control points of the structural plane proposed by the present invention, based on Figure 9For the trend line of the modeling target area, generate equally spaced sampling points, such as Figure 11 (a). Using step 2) of Section 2.2, obtain the layer control points at different depths, such as Figure 11 (b). After generating the layer data file, import it into the modeling software for rendering. The effect after local magnification of the layer is as shown in Figure 11 (c). Overlay the bottom surface of the loose layer in the modeling target area, and perform TIN layer intersection line tracking and unconformity surface topological relationship reconstruction. The effects are respectively as shown in Figure 11 (d) and Figure 11 (e). After completing the layer clipping, the final structural surface model of the modeling target area is as shown in Figure 11 (f).

[0102] In summary, in view of the lack of basic modeling data such as boreholes, cross-sections, and seismic interpretations in 3D geological modeling work, and the difficulty of forming effective modeling constraints by conventional methods, this paper proposes a convenient structural modeling method using low-cost and easily accessible data such as plane geological maps and tectonic outline maps.

[0103] Using strike control and attitude constraints, this invention proposes an algorithm for batch generation of structural surface control points, and based on this, forms an algorithmic process for 3D layer creation. This method can make full use of the limited geological data in the modeling target area, quickly generate layer control point data, realize structural modeling, and solve the pain point of being unable to model due to insufficient data.

[0104] Based on mature computer graphics basic models such as collision detection, this invention realizes the algorithm for reconstructing the spatial topological relationship of structural surfaces, and forms the entire process from data generation, storage to 3D visualization of the model. The method process is simple, easy to program and implement, and considers the compatibility and expansion of data formats, facilitating data sharing.

[0105] The application examples of this invention show that relying on the attitude information and structural features extracted from the plane geological map, using the method of this paper effectively solves problems such as the special structure of steeply inclined strata in the target area, the lack of basic data such as boreholes, and the inability to realize structural modeling by conventional methods, verifying the feasibility and reliability of the method.

[0106] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention. The structures, devices, and operation methods not specifically described and explained in the present invention are implemented according to the conventional means in the art without special description and limitation.

Claims

1. A three-dimensional modeling method for structural surfaces based on strike control and occurrence constraints, characterized in that: The following steps are involved: Step 1: Extract the stratigraphic and fault strike lines and occurrence data based on the planar geological map of the target area, and extract the coordinates and azimuth of the points on the strike line through the curvature-constrained equidistant sampling method; combine the occurrence data and generate layer control points in batches by changing the sampling depth; then use the stratigraphic or fault stratification points revealed by the borehole or geological profile to expand and constrain the layer control points, establish the structural surface control point set, and then create the irregular triangulated network of the structural surface; Step 2: Use collision detection algorithm to extract the intersection lines of structural surfaces, reconstruct the contact relationship of strata and the spatial topological relationship of fault networks, and finally obtain the layer data file of the structural model, and complete the data storage and three-dimensional visualization of the model.

2. A three-dimensional modeling method for structural surface based on strike control and occurrence constraint according to claim 1, characterized in that: The step 1 specifically comprises the following steps: Step 101: Extract the strike line of the layer from the planar geological map of the study area with the help of GIS software. The layer is the stratigraphic interface or fault plane. Sample the strike line at equal intervals and extract the coordinate set of the sampling points {P i }={X i ,Y i ,Z0,θ i }, where i = 0, 1, ..., n; Z0 represents the initial elevation value or depth value of the layer object sampling, which can be extracted through the trend line and DEM, and the corresponding initial depth can be set in the depth domain modeling, θ i Indicates the azimuth of the trend line at the point, and the data table of the trend line is obtained by sorting; At the same time, curvature constraints are used to control the morphology of sampling points and perform data thinning. The principle is as follows: the equation of the B-spline curve passing through the trend sampling points can be quickly obtained through data fitting, and it is set as y = f(x); taking one end of the trend line as the starting point, every three adjacent points are used to calculate the curvature and traverse to the other end of the curve; P i (x i ,y i ) point, the curvature can be calculated by formula (1); Thus, we can obtain the curvature set of the trend line {κ1,κ2,…,κ n-2 }; Step 102: Take the strike line as the top boundary line of the layer, and construct a right triangle in three-dimensional space for the strike line sampling points according to the layer inclination and dip angle. Let the layer dip angle be α, and the sampling point P i The projection point at depth h is denoted as O. According to the trigonometric function relationship of plane projection, P′ at depth h i The coordinates (X′ i ,Y′ i ) is related to the dip and strike azimuth of the layer; when the dip of the layer is southward, P′ i The coordinates (X′ i ,Y′ i ) can be calculated using formulas (2) and (3) respectively. When the stratum dips to the north, P′ i The coordinates (X′ i ,Y′ i ) are calculated using formulas (4) and (5) respectively; Based on this, we can get the control point set {P′ i }={X′ i ,Y′ i ,Z h },i=0,1,…,n,Z h Represents P′ i The elevation corresponding to the depth h at the position can be directly converted to the depth Z if depth domain modeling is used. h = h; By controlling the sampling interval of the strike line and setting different sampling depths, the layer control point set {P i ,P′ i ,P″ i ,…}; Based on this, the layer control points can be generated quickly in batches with the help of the constraints of layer occurrence and sampling depth; at the same time, the positive and negative value changes of the sampling depth h can flexibly adjust the sampling range of the layer control points; Step 103: Insert and correct the drilling layer points; insert the layer points extracted from the drilling into the control point set {P i ,P′ i ,P″ i ,…}, encrypt and correct the layer control point set.

3. The method for three-dimensional modeling of a structural surface based on strike control and occurrence constraint according to claim 2, characterized in that: The step 1 also includes: Step 104: Store the layer control point data, where the Type field records the control point set {P i ,P′ i ,P″ i ,…}′ represents the layer type; Horizon represents the stratigraphic layer, Fault represents the fault, the Name field records the layer name, and the PointID field records the index number of the control point; in the same layer, PointID is a unique value, and X, Y, and Z respectively record the coordinate values ​​of the control point.

4. The method for three-dimensional modeling of a structural surface based on strike control and occurrence constraint according to claim 3, characterized in that: The step 1 also includes the following steps: Step 105: Create a layer triangulated network; after obtaining the layer control point set, directly perform Delaunay triangulation to generate the layer TIN; save the vertex coordinates and triangle connection order in the TIN in text format.

5. The method for three-dimensional modeling of a structural surface based on strike control and occurrence constraint according to claim 4, characterized in that: After the layer TIN is established according to step 1, step 2 constructs the contact relationship of the strata and the topological relationship of the fault network. The principle of the method is as follows: Step 201: Use the collision detection algorithm to screen the intersecting triangles. Use the axially aligned bounding box collision detection method to extract the intersection lines of the TINs of the two intersecting planes; construct a cube with edges parallel to the coordinate axis to completely enclose the triangle, and use a recursive method to form a multi-branch tree structure to implement the collision detection algorithm; Step 202: Extraction of intersection lines: Since all layers are constructed using triangulated meshes, after being screened by the axially aligned bounding box collision detection method, the extracted intersection lines are converted into triangle intersection operations; when the bounding boxes of two triangles overlap, the triangles are considered to intersect, and the edges of one triangle are judged with the other triangle to obtain the intersection points; the intersection lines are calculated using the Mohler-Tremblay algorithm, and the principle is as follows: Assume that the plane is represented by the normal vector N and a point p′ on the plane where the triangle is located. Then any point p on the plane satisfies formula (6): p:(pp′)·N=0(6) With starting point O and direction vector represents the ray, t represents the ray length ratio, then the ray r satisfies formula (7): Assume that the three vertices of the triangle are P0, P1, and P2. According to the rule of barycentric coordinate interpolation, the intersection of the ray and any point on the plane can be expressed as formula (8): In the formula, b1, b2, 1-b1-b2 represent the centroid coordinates of any point on the triangle plane; P0, P1, P2 represent the coordinates of the three vertices of the triangle; in formula (8), there are only three unknown quantities: t, b1, and b2. When t, b1, b2, and 1-b1-b2 are all greater than or equal to 0, the ray intersects the triangle; Based on this, the vertices of the two triangles are traversed separately, and each pair of triangles requires six calculation operations. When t∈[0,1], it means that the edge of the current triangle has an intersection with another triangle; the obtained intersection point and its adjacent intersection points are recorded; if there is no adjacent point, it is the endpoint of the intersection line; A linked list of intersections is formed, and the intersection lines of the two levels can be obtained by tracing all the intersections in sequence; all intersections are marked with a unique PointID. While saving the point coordinates, the Neighbors field is set to save the front and back neighboring points of the current point. The endpoint has only one neighboring point; according to the order of tracing the intersections, the first one recorded is the starting point, and the last one recorded is the end point of the intersection.

6. The method for three-dimensional modeling of a structural surface based on strike control and occurrence constraint according to claim 5, characterized in that: The step 2 further comprises the following steps: Step 203: spatial topology reconstruction and layer clipping; based on the TIN data created above; including layer control point coordinates and triangle vertex indexes, construction and layer clipping of spatial topological relationships such as stratum contact relationships and fault networks, conversion into the division of layer control point sets by intersection lines, and reconstruction of the triangulated network after division; First, the master-slave relationship of the two planes is set according to the stratigraphic column or fault cutting relationship, and then the bounding box is generated by the intersection point set, which is the slave surface at this time; when the stratigraphic contact relationship is reconstructed, if the strata are in an erosion relationship, the master surface is the unconformity surface and the slave surface is the older stratum; when the fault network is constructed, it is determined by the order or priority of the fault development; when the fault cuts the stratum, the control point set with the fault plane as the master stratigraphic plane and the slave plane as the slave will be divided into three parts: the inside of the bounding box and the two sides of the bounding box; then, according to the coordinate range, the slave surface control points outside the bounding box are divided into two groups AB, and any point is selected from them; suppose a point P is taken from group A, and it forms a point pair with the slave surface control point Q inside the intersection bounding box, and the Moller-Tremblay algorithm is used again to determine whether the connecting line PQ has an intersection with the master surface; if the intersection exists, it should be classified into group B with point Q, otherwise it should be classified into group A; after the classification is completed, according to the erosion relationship, the data points on one side of the slave surface point are retained to complete the surface clipping; By using the intersection detection and layer clipping algorithm, the TIN file is traversed to make pairwise judgments, and the topology reconstruction is completed while the intersection extraction is completed.

7. The method for three-dimensional modeling of a structural surface based on strike control and occurrence constraint according to claim 6, characterized in that: The step 2 also includes the following steps: Step 204: TIN repair; adding the intersection points to the clipped secondary surface point set, using the Delaunay triangulation algorithm to locally repair the TIN, and updating the layer data file.

Images