A method for three-dimensional modeling of structural planes based on strike control and attitude constraint

Abstract

The present application belongs to the technical field of three-dimensional geological modeling, and particularly relates to a three-dimensional modeling method of structural planes based on strike control and occurrence constraint. The method extracts strata, fault strike lines and occurrence data based on a plane geological map of a modeling target area, extracts coordinates and azimuth angles of points on the strike lines through an equal-interval sampling method with curvature constraint, generates layer control points in batches through transformation of sampling depth in combination with the occurrence data, expands the layer control point set by using strata or fault stratification points exposed by drilling or geological profiles, and further creates an irregular triangular net of structural planes. The method extracts intersection lines of the structural planes by using a collision detection algorithm, reconstructs strata contact relations and spatial topological relations of fault networks, finally obtains layer data files of the structural model, and completes data storage and three-dimensional visualization of the model. The present application can quickly establish a structural plane model of strata, faults and the like and reconstruct spatial topological relations, and realize generation of data files of the structural model and three-dimensional visualization.

Description

Technical Field

[0001] This invention belongs to the field of three-dimensional geological modeling technology, specifically relating to a three-dimensional modeling method for structural surfaces based on strike control and attitude constraints. Background Technology

[0002] Geological support systems are information technology pillars for promoting safe coal mine production and an important component of smart mine construction. Three-dimensional geological modeling is a key technology in the construction of mine geological support systems and serves as the carrier and data foundation for building transparent geology, digital twins, and intelligent management platforms under the intelligent mine strategy. 3D geological modeling uses computer technology to express the spatial morphology and internal attribute fields of geological bodies in three-dimensional space. It can not only display the spatial structure of geological bodies but also utilize spatial information analysis techniques and geostatistics tools to achieve the integration of diverse heterogeneous data, process simulation, and analysis and prediction. It has been widely applied in numerical simulation, disaster prevention, and other scenarios.

[0003] Structural modeling primarily focuses on creating structural surfaces such as fault planes, bedding planes, and unconformities. The modeling process centers on the creation of 3D surfaces. Constructing a structural model that conforms to actual geological features such as occurrence descriptions, contact relationships, and fault networks using known data is the most crucial step in 3D geological modeling. This can generally be achieved through explicit modeling, implicit modeling, or a combination of both. Borehole data, representing measured stratigraphic structures, is one of the most reliable data sources for 3D geological modeling. Establishing bedding plane triangulations using borehole data is a common method, and relatively uniformly distributed boreholes can produce better structural surface generation results, thus it is widely used in various 3D geological modeling problems. 3D seismic data, due to its accurate representation of the spatial distribution of deep subsurface structural surfaces, is also an important data source for 3D geological modeling. However, due to its high initial investment cost, it is mostly used in models with complex structures, large burial depths, and high accuracy requirements. Furthermore, cross-sectional views can visually reflect stratigraphic structures and their topological relationships and are also frequently used as important data sources for structural modeling. Using a sufficient number of cross-sectional views can also yield a good 3D geological model. Clearly, borehole data, 3D seismic interpretation results, and geological profiles all play a very important role in structural modeling. 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 principle, and can meet the needs of most structural modeling scenarios.

[0004] However, in real-world work, situations arise where one or more of the aforementioned basic data are lacking, or the data distribution characteristics are unsatisfactory, making it difficult to meet the needs of structural modeling. This hinders the implementation of existing structural modeling techniques and makes it impossible to accurately describe the spatial variation characteristics of geological bodies. Specifically, firstly, sparse borehole data often fails to effectively constrain interpolation algorithms during implicit modeling, leading to easy drift and distortion of structural surfaces, increasing uncertainty. Secondly, the lack of sufficient profile data makes it difficult to generate correct spatial topological relationships for stratigraphic contact relationships and fault networks. Furthermore, if the target area for modeling lacks corresponding 3D seismic interpretation data, conventional modeling methods and software become ineffective due to the lack of necessary basic data. Moreover, when modeling specific geological structures, the lack of data further exacerbates the uncertainty of modeling, and existing modeling methods may even be insufficient. For example, steeply dipping strata are widely present on the flanks of folds subjected to intense tectonic compression. For instance, steeply dipping coal seams may be the main mining seams, but their large dip angles hinder efficient coal mining, and the extremely limited exposure of steeply dipping strata by boreholes presents new challenges for conventional structural modeling techniques. In such cases, virtual drilling can be used to convert the original vertical boreholes into horizontal virtual boreholes, using software deception for structural modeling. However, this method requires high data preprocessing. When the geological structure is complex or the number of boreholes is large, manual data processing is necessary, which results in a high error rate and low efficiency. When the number of boreholes is limited, it is difficult to play an effective role, and it has strong limitations.

[0005] Planar geological maps are readily available, low-cost data sources that directly reflect regional geological characteristics. They incorporate expert knowledge and experience, providing both a macro-geological background for geological modeling and compensating for the lack of detailed modeling data such as borehole data. They are also widely used in 3D geological modeling. Existing technologies have researched 3D modeling methods for complex faults based on planar geological maps. Furthermore, existing technologies have proposed research approaches that comprehensively utilize low-cost data such as cross-sectional maps and geological maps for 3D geological modeling. Gao Shijuan, Zhou Liangchen, and others have studied a 3D geological modeling method based on planar geological maps and using map-cut cross-sections as a framework, pointing out that using planar geological maps for 3D geological modeling is an effective solution when other geological data is scarce.

[0006] In view of this, this invention provides a 3D modeling method for structural surfaces based on strike control and attitude constraints. Through analysis and summarization in practical work, it proposes a solution for structural modeling based on low-cost, fundamental data obtained from planar geological maps or structural outline maps. It employs a layer control point generation algorithm based on strike control and attitude constraints, as well as a spatial topology reconstruction algorithm. Using a computer program, it generates structural surface control point data in batches, rapidly establishing stratigraphic, fault, and other structural surface models and reconstructing spatial topological relationships, thus achieving data file generation and 3D visualization of the structural model. Application examples verify the feasibility of the method, providing a solution for 3D geological modeling in sparse data areas. Summary of the Invention

[0007] The purpose of this invention is to provide a 3D modeling method for structural surfaces based on strike control and attitude constraints. This method employs a strike control point generation algorithm and a spatial topology reconstruction algorithm based on strike control and attitude constraints. It utilizes a computer program to batch generate structural surface control point data, rapidly establishing structural surface models such as strata and faults, and reconstructing spatial topology relationships. This enables the generation of data files for the structural models and 3D visualization. Application examples verify the feasibility of the method, providing a solution for 3D geological modeling in sparse data areas.

[0008] The specific technical solution adopted by this invention is as follows:

[0009] A 3D modeling method for structural surfaces based on orientation control and attitude constraints includes the following steps:

[0010] Step 1: Extract stratigraphic, fault strike lines, and attitude data from the planar geological map of the target area. Extract the coordinates and azimuths of points on the strike lines using a curvature-constrained equal-interval sampling method. Combine the attitude data with batch generation of bedding control points by changing the sampling depth. Then, expand the bedding control point set using stratigraphic or fault layering points revealed by limited boreholes or geological profiles, thereby creating an irregular triangular network of structural surfaces.

[0011] Step 101: Using GIS software, extract the strike lines of the strata from the planar geological map of the study area. The strata are stratigraphic interfaces or fault planes. Sampling is performed on the strike lines at equal intervals, and the coordinate set {P} of the sampling points is extracted. i}={X i ,Y i ,Z0,θ i}, where i = 0, 1, ..., n; Z0 represents the initial elevation or depth value of the sampled layer object, which can be extracted by the orientation line and the DEM. In depth domain modeling, the corresponding initial depth can be set, θ i This represents the azimuth angle of the directional line at that point, and the data table of the directional line is obtained by organizing the data.

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

[0013]

[0014] This yields the curvature set of the directional line {κ1,κ2,…,κ}. n-2};

[0015] Step 102: Using the strike line as the top boundary line of the layer, construct right triangles in three-dimensional space for the sampling points along the strike line based on the layer dip 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 relationships of plane projection, P′ at depth h is... i coordinates (X′) i ,Y′ i The angle of P′ is related to the dip and strike azimuth of the stratum; when the stratum dips southward, P′ i coordinates (X′) i ,Y′ i P′ can be calculated using formulas (2) and (3) respectively, when the strata dip northward. i coordinates (X′) i ,Y′ i The results were obtained using formulas (4) and (5) respectively.

[0016]

[0017] Based on this, the set of control points {P′} at depth h is obtained. i}={X′ i ,Y′ i Z h}, i = 0, 1, ..., n, Z h P′ represents i The elevation corresponding to depth h at a given location can be directly converted to depth Z if depth domain modeling is used. h =h; By controlling the sampling interval of the traverse lines and setting different sampling depths, the set of layer control points {P} can be obtained. i ,P′ i ,P″ i Based on this, by leveraging the constraints of bedding attitude and sampling depth, bedding control points can be generated rapidly in batches; simultaneously, the sampling depth... The positive and negative values ​​can be varied to flexibly adjust the range of layer control point sampling;

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

[0019] Step 104: Store the layer control point data, where the Type field records the control point set {P} i ,P′ i ,P″ i The type of the layer is represented by ,…}′; Horizon represents the ground plane, 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 record the coordinate values ​​of the control point respectively.

[0020] Step 105: Create the layer triangulation 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.

[0021] Step 2: Use collision detection algorithm to extract structural surface intersections, reconstruct stratigraphic contact relationships and fault network spatial topology, and finally obtain the layer data file of the structural model, and complete the data storage and 3D visualization of the model.

[0022] Step 2, based on the establishment of the layer TIN in Step 1, constructs the contact relationships and fault network topology of the strata. The principle of the method is as follows:

[0023] Step 201: Use a collision detection algorithm to filter intersecting triangles. The intersection line of the TINs on two intersecting layers is extracted using an axially aligned bounding box collision detection method. The collision detection algorithm is implemented by constructing a cube with edges parallel to the coordinate axes to completely enclose the triangles and using a recursive multi-branch tree structure.

[0024] Step 202: Intersection Line Extraction: Since all layers are constructed using triangular meshes, after filtering using the axially aligned bounding box collision detection method, the extracted intersection lines are transformed into triangle intersection operations. When the bounding boxes of two triangles overlap, they are considered to intersect. The edges of one triangle are compared with the other triangle to obtain the intersection point. The Moeller-Trumble algorithm is used to calculate the intersection line, and the principle is as follows:

[0025] Let the plane be represented by the normal vector N and a point p′ on the plane containing the triangle. Then any point p on the plane satisfies formula (6):

[0026] p:(pp′)·N=0(6)

[0027] Using the starting point O and the direction vector Let denot be a ray, and t be the proportion of the ray length. Then, ray r satisfies 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 a ray and any point on the plane can be expressed as formula (8):

[0030]

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

[0032] Accordingly, the vertices of the two triangles are traversed separately, requiring six calculation operations for each pair of triangles. When t∈[0,1], it indicates that the edge of the current triangle intersects with the other triangle. The obtained intersection points and their adjacent intersection points are recorded. If there are no adjacent points, they are the endpoints of the intersection line. An intersection point linked list is formed, and the intersection line of the two levels can be obtained by tracing all intersection points in sequence. All intersection points are marked with a unique PointID. While saving the point coordinates, the Neighbors field is set to save the previous and next adjacent points of the current point. The endpoint has only one adjacent point. According to the tracing order of the intersection points, the first one recorded is the starting point, and the last one recorded is the endpoint of the intersection line.

[0033] Step 203: Spatial topology reconstruction and layer clipping; Based on the TIN data created above; Constructing and clipping spatial topological relationships such as layer control point coordinates and triangle vertex indices, stratigraphic contact relationships, fault networks, etc., and transforming it into the problem of dividing the layer control point set by intersection lines and reconstructing the triangular network after division.

[0034] First, the master-slave relationship between two surfaces is established based on the stratigraphic column or fault cutting relationship. Then, a bounding box is generated using the intersection point set, which is the slave surface. When reconstructing the stratigraphic contact relationship, if the stratigraphy is in an erosion relationship, the master surface is the unconformity, and the slave surface is the older stratigraphy. When constructing the fault network, the order or priority of fault development is used to determine the slave surface. When a fault cuts the stratigraphy, the control point set with the fault surface as the master and the stratigraphic surface as the slave will be divided into three parts: inside the bounding box and on both sides of the bounding box. Subsequently, the slave surface control points outside the bounding box are divided into two groups, A and B, based on the coordinate range, and one point is randomly selected from them. Assuming a point P is selected from group A, it forms a point pair with the slave surface control point Q inside the intersection bounding box. The Moeller-Trumble algorithm is used again to determine whether the line connecting PQ intersects with the master surface. If an 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, the data points on one side of the slave surface point are retained according to the erosion relationship to complete the surface trimming.

[0035] By using intersection detection and layer clipping algorithms, the TIN file is traversed to make pairwise judgments, thereby completing the extraction of intersection lines and the reconstruction of the topology.

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

[0037] The technical effects achieved by this invention are as follows:

[0038] This invention utilizes strike control and attitude constraints to propose a batch generation algorithm for structural surface control points, and based on this, forms an algorithmic flow for creating 3D layers. This method can fully utilize the limited geological data of the modeling target area to quickly generate layer control point data, achieving structural modeling and solving the pain point of being unable to model due to insufficient data.

[0039] This invention, based on mature computer graphics fundamentals such as collision detection, implements an algorithm for reconstructing the spatial topological relationships of constructed surfaces, forming a complete process from data generation and storage to 3D model visualization. The method is simple, easy to program, and considers data format compatibility and extensibility, facilitating data sharing.

[0040] The application examples of this invention demonstrate that, relying on the occurrence information and structural features extracted from the planar geological map, the method presented in this paper effectively solves the problems of special steeply dipping strata structures in the target area, lack of basic data such as borehole data, and the inability of conventional methods to achieve structural modeling, thus verifying the feasibility and reliability of the method. Attached Figure Description

[0041] Figure 1 This is a flowchart of the method of the present invention;

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

[0043] Figure 3 This is a schematic diagram of the calculation of the layer sampling points under the attitude constraint of the present invention;

[0044] Figure 4 This is a schematic diagram illustrating the principle of sampling point calculation in this invention;

[0045] Figure 5 These are the drilling layer data encryption and correction layer control points of this invention;

[0046] Figure 6 This is a schematic diagram of the topology trimming of the present invention;

[0047] Figure 7 This is the structural diagram of the Weier coalfield and the data acquisition results of the main fault planes in Example 1 of this invention;

[0048] Figure 8 This is the three-dimensional structural model of the Weier mining area in Example 1 of this invention;

[0049] Figure 9 This refers to the bedrock strata attitude of the target area modeled in Example 2 of this invention;

[0050] Figure 10 This is a partial geological profile of the target area modeled in Example 2 of this invention;

[0051] Figure 11 This is the modeling process and results of a steeply dipping stratum in Example 2 of this invention. (a) is the result of sampling points extracted from the strike line of the target area; (b) is the result of sampling points for the attitude-constrained strata control points proposed in this paper; (c) is the steeply dipping strata and fault structure surface rendered based on the strata control points; (d) is the extraction of the intersection line of the two strata; (e) and (f) are the results of strata intersection line tracing and the structural model after topological reconstruction and strata trimming, respectively. Detailed Implementation

[0052] To make the objectives and advantages of this invention clearer, the invention will be specifically described below with reference to embodiments. It should be understood that the following text is merely used to describe one or more specific embodiments of the invention and does not strictly limit the scope of protection specifically claimed by the invention.

[0053] Example 1:

[0054] A 3D modeling method for structural surfaces based on orientation control and attitude constraints includes the following steps:

[0055] Step 1: Extract stratigraphic, fault strike lines, and attitude data from the planar geological map of the target area. Extract the coordinates and azimuths of points on the strike lines using a curvature-constrained equal-interval sampling method. Combine the attitude data with batch generation of bedding control points by changing the sampling depth. Then, expand the bedding control point set using stratigraphic or fault layering points revealed by limited boreholes or geological profiles, thereby creating an irregular triangular network of structural surfaces.

[0056] Step 2: Based on computational geometry and computer graphics, a collision detection algorithm is used to extract the intersection lines of structural surfaces, and the stratigraphic contact relationships and fault network spatial topology are reconstructed. Finally, the layer data file of the structural model is obtained, and the data storage and 3D visualization of the model are completed.

[0057] In this invention, triangles are the basic unit when constructing 3D surfaces in a computer. Connecting three adjacent points forms triangles, and these triangles, in turn, form irregular surfaces. Therefore, creating a Triangulated Irregular Network (TIN) is essential for structural surface modeling. TIN-based geological modeling can adjust the size and number of triangles according to the complexity of the modeling target, preserving topological relationships while easily handling complex geological structures. This invention is based on this core idea of ​​structural surface TIN construction. Using planar geological maps as basic data, it proposes an automatic batch generation algorithm for control points of strata and fault planes to address the problem of insufficient basic modeling data. Based on this, a structural surface topology reconstruction algorithm is designed, forming a method for generating and modeling structural surface data in situations where basic modeling data is lacking or unsatisfactory. The method also considers the universality of data formats and storage interaction issues. The method flow of this invention is as follows: Figure 1 As shown, the core work is as follows: First, based on the planar geological map of the target area, stratigraphic, fault strike lines, and attitude data are extracted. Using a curvature-constrained, equally spaced sampling method, the coordinates and azimuths of points on the strike lines are extracted. Combined with the attitude data, layer control points are generated in batches by changing the sampling depth. Then, using the stratigraphic or fault layering points revealed by limited boreholes or geological profiles, the set of layer control points is expanded, thereby creating an irregular triangular network of structural surfaces. Next, based on computational geometry and computer graphics concepts, a collision detection algorithm is used to extract the intersection lines of structural surfaces and reconstruct spatial topological relationships such as stratigraphic contact relationships and fault networks. Finally, the layer data file of the structural model is obtained, and the model's data storage and 3D visualization are completed.

[0058] Preferably, given that existing geological modeling software struggles to accurately and effectively construct 3D models when necessary modeling data such as borehole data is lacking, this paper proposes a 3D modeling algorithm for structural surfaces based on strike control and attitude constraints. Step 1 specifically includes the following steps:

[0059] Step 101: Using GIS software, extract the strike lines of the strata from the planar geological map of the study area. The strata are stratigraphic interfaces or fault planes. Sampling is performed on the strike lines at equal intervals, and the coordinate set {P} of the sampling points is extracted. i}={X i ,Y i ,Z0,θ i}, where i = 0, 1, ..., n; Z0 represents the initial elevation or depth value of the sampled layer object, which can be extracted by the orientation line and the DEM. In depth domain modeling, the corresponding initial depth can be set, θ i This represents the azimuth angle of the directional line at that point, and the data table of the directional line is obtained by organizing the data.

[0060] Meanwhile, to avoid abrupt changes in direction during sampling, or computational redundancy due to excessively dense data points, curvature constraints are used to control the shape of the sampling points and thin out the data. The principle is as follows: the equation of the B-spline curve passing through the sampling points can be quickly obtained through data fitting, let it be y = f(x); taking one end of the direction line as the starting point, every three adjacent points are used to calculate the curvature, traversing to the other end of the curve; for Figure 2 P shown i (x i ,y i The curvature at point () can be calculated using formula (1);

[0061]

[0062] This yields the curvature set of the directional line {κ1,κ2,…,κ}. n-2}; By setting a reasonable threshold for r in equation (1), the trend of change of the strike line can be controlled by curvature, so as to avoid the drastic change of the strike of the strata. At the same time, in the area where the strike change is small, the data points can be thinned to reduce the amount of calculation in the subsequent steps.

[0063] Step 102: Using the strike line as the top boundary line of the layer, construct right-angled triangles in three-dimensional space for the sampling points along the strike line based on the layer's dip and dip angle; for example... Figure 3 As shown; let the dip angle of the plane be α, and the sampling point P i The projection point at depth h is denoted as O. According to the trigonometric relationships of plane projection, P′ at depth h is... i coordinates (X′) i ,Y′ i It is related to the dip and azimuth of the stratum; when the stratum dips southward, such as Figure 4 -a and Figure 4 As shown in -b, P′ i coordinates (X′) i ,Y′ iThe results can be calculated using formulas (2) and (3) respectively. When the strata dip northward, such as Figure 4 -c and Figure 4 As shown in -d, P′ i coordinates (X′) i ,Y′ i The results were obtained using formulas (4) and (5) respectively.

[0064]

[0065] Based on this, the set of control points at depth h of the ground stratum is obtained. i = 0, 1, ..., n, Z h P′ represents i The elevation corresponding to depth h at a given location can be directly converted to depth Z if depth domain modeling is used. h =h; By controlling the sampling interval of the traverse lines and setting different sampling depths, the set of layer control points {P} can be obtained. i ,P′ i ,P″ i Based on this, by leveraging the constraints of bedding attitude and sampling depth, bedding control points can be generated rapidly in batches; simultaneously, the sampling depth... The positive and negative values ​​can be varied to flexibly adjust the range of layer control point sampling;

[0066] Step 103: Insertion and Correction of Borehole Stratigraphic Points; Stratigraphic stratification points created based on borehole columnar sections are relatively high-precision control point data. Even a small number of borehole stratification data can still play a role in correcting the stratigraphic layers. Therefore, the stratigraphic layers extracted from the boreholes are inserted into the control point set {P} created in Step 102. i ,P′ i ,P″ i ...}, to encrypt and correct the set of layer control points.

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

[0068] Table 1. Data Structure of Level Control Point Set

[0069]

[0070] Step 105: Create the layer triangulation 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, as shown in Table 2; the layer point data saved in Table 2 can be directly rendered as a 3D layer, or it can be imported as point cloud data into other 3D geological modeling software to achieve data interaction.

[0071] Table 2. Level TIN Data Storage Format

[0072]

[0073] Preferably, in the construction modeling, the determination of the intersection line of the strata is crucial. Taking the intersection line of strata and faults as an example, the general method is to obtain the intersection line by densifying the strata data and based on the fault's influence range and elevation difference, which is the overall fault interpolation method. This method requires manual intervention to specify the influence range of the strata, and the algorithm is complex and has poor adaptability. Step 2, based on the establishment of the strata TIN in step 1, constructs the contact relationship of the strata and the fault network topology. The principle of the method is as follows:

[0074] Step 201: Use a collision detection algorithm to filter intersecting triangles. The Axis-aligned bounding box (AABB) collision detection method is used to extract the intersection line of the TINs of two intersecting layers. Based on computational geometry and computer graphics, the AABB uses triangles as the smallest unit. It completely encloses the triangles by constructing cubes with edges parallel to the coordinate axes. A recursive multi-branch tree structure is used to implement the collision detection algorithm. It only needs to determine whether the projections of the two bounding boxes on the coordinate axes coincide to determine if the triangles intersect. This method is simple to construct, requires little storage space, and has low algorithm complexity, enabling rapid filtering of intersecting triangles.

[0075] Step 202: Intersection Line Extraction: Since all layers are constructed using triangular meshes, after filtering using the axially aligned bounding box collision detection method, the extracted intersection lines are transformed into triangle intersection operations. When the bounding boxes of two triangles overlap, they are considered to intersect. The edges of one triangle are compared with the other triangle to obtain the intersection point. The Moeller-Trumble 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 containing the triangle. Then any point p on the plane satisfies formula (6):

[0077] p:(pp′)·N=0(6)

[0078] Using the starting point O and the direction vector Let denot be a ray, and t be the proportion of the ray length. Then, ray r satisfies formula (7):

[0079]

[0080] Let the three vertices of the triangle be P0, P1, and P2. According to the rules of barycentric coordinate interpolation, the intersection of a 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 centroid coordinates of any point on the plane of the triangle; P0, P1, P2 represent the coordinates of the three vertices of the triangle; in formula (8), there are only three unknowns: t, b1, b2. When t, b1, b2, 1-b1-b2 are all greater than or equal to 0, the ray intersects the triangle.

[0083] Accordingly, the vertices of the two triangles are traversed separately, requiring six calculation operations for each pair of triangles. When t∈[0,1], it indicates that the edge of the current triangle intersects with the other triangle. The obtained intersection points and their adjacent intersection points are recorded. If there are no adjacent points, they are the endpoints of the intersection line. An intersection point linked list is formed, and the intersection lines of the two levels can be obtained by tracing all intersection points in sequence. The intersection point data storage format is shown in Table 3. All intersection points are marked with a unique PointID. While saving the point coordinates, the Neighbors field is set to save the previous and next adjacent points of the current point. The endpoint has only one adjacent point. According to the tracking order of the intersection points, the first one recorded is the starting point, and the last one recorded is the endpoint of the intersection line.

[0084] Table 3. Storage format of intersection points

[0085]

[0086] Step 203: Spatial topology reconstruction and layer clipping; Based on the TIN data created above; Constructing and clipping spatial topological relationships such as layer control point coordinates and triangle vertex indices, stratigraphic contact relationships, fault networks, etc., and transforming it into the problem of dividing the layer control point set by intersection lines and reconstructing the triangular network after division.

[0087] First, based on the stratigraphic column or fault cutting relationship, the master-slave relationship between the two surfaces is established. Then, a bounding box is generated using the intersection point set, where the surface is the slave surface. When reconstructing the stratigraphic contact relationship, if the stratigraphy is in an erosion relationship, the master surface is the unconformity, and the slave surface is the older strata. When constructing the fault network, the order or priority of fault development is used to determine the control point set (fault surface is master, and stratigraphic surface is slave) will be divided into three parts: the inside of the bounding box and the two sides of the bounding box. Subsequently, based on the coordinate range, the slave surface control points outside the bounding box are divided into two groups, A and B, and one point is randomly selected from them. Assuming a point P is selected from group A, such as... Figure 6 As shown; the control point Q inside the bounding box of the intersection line forms a point pair, and the Moeller-Trenbre algorithm is used again to determine whether the connecting line PQ intersects 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 surface point are retained to complete the surface trimming.

[0088] By using intersection detection and layer clipping algorithms, the TIN file is traversed to make pairwise judgments, thereby completing the extraction of intersection lines and the reconstruction of the topology.

[0089] Step 204: TIN Repair; Add the intersection points to the clipped set of surface points, use the Delaunay triangulation algorithm to locally repair the TIN, and update 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, with secondary structures dominated by faults. The Weier mining area is located in the south-central part of the Weizhou mining area, specifically in the southern part of the eastern wing of the Weizhou Syncline. It is a monocline structure with a westward dip of 5–35°, and the strata gradually flatten from east to west. The mining area contains two sets of relatively large and numerous oblique faults: a north-northwest fault (mainly compressive reverse faults) and a north-northeast fault (mainly extensional normal faults). More than 40 faults have been identified within the mining area, including 11 large faults with considerable extension lengths. Table 4 lists the attitudes of these faults.

[0092] Table 4 Fault occurrence in the Weierjing field

[0093]

[0094]

[0095] Figure 7 (a) shows the main geological structure and coal seam boundaries of the Weier coalfield. Figure 7 (b) To utilize the proposed algorithm for batch generation of control points for structural surfaces, based on... Figure 7 (a) shows the fault strike line and the fault attitude data in Table 4, which are used to generate the main fault plane control points for the well field. Figure 7(c) Fault structural surfaces generated using bedding plane control points. It is worth noting that during the batch generation of fault bedding plane control points, the sampling depth h was adjusted accordingly when calculating bedding plane points based on the strike line, taking into account the characteristics of the mine field elevation. Therefore, in Figure 7 (a) and Figure 7 In (b), the orientation line is shown as the initial sampling position near the middle of the plane.

[0096] The Weier coalfield has relatively well-developed coal seams, with coal-bearing strata consisting of the Lower Permian Shanxi Formation (P1s) and the Carboniferous-Permian Taiyuan Formation (C2-P1t). The mineable coal seams in the Taiyuan Formation are mainly concentrated in the second and third sections, with a total of six mineable seams (coal seams 12, 14, 15, 16, 17, and 20). The Shanxi Formation has three mineable coal seams (coal seams 2, 3, and 4). Using the method proposed in this paper, combined with coalfield geological borehole data, structural surface data of the mineable coal seams were generated based on the same principle.

[0097] To support the selection of coalbed methane well locations in the research on the coordinated development of coal and coalbed methane within the mining area, a three-dimensional structural model was constructed for the initial mining area of ​​the Wei'er mining area. After obtaining the aforementioned data, a three-dimensional geological model of the initial mining area of ​​the Wei'er mining area was further constructed, showcasing the spatial distribution characteristics of the geological structure of the mining area. Figure 8 (a) shows the boundaries of the early mining area and the distribution of boreholes in the well field. Figure 8 (b) shows the structural model of the early mining area, in which the spatial relationship between the Shanxi Formation and Taiyuan Formation coal seams and the main faults in the mining area is clearly and intuitively displayed.

[0098] Example 2:

[0099] Another target area in the actual modeling work of this invention has a rather unique geological structure. The bedrock strata in this area are subjected to strong tectonic compression, exhibiting steeply dipping monocline strata in unconformable contact with the overlying loose layers. A simplified geological plan of the target area is shown below. Figure 9 As shown, the bedrock strata and faults all trend northeast, as shown in the cross-section diagram. Figure 10 The data shows that the bedrock strata dip at an angle greater than 80 degrees, and are in a nearly vertical state.

[0100] Due to the steep dip of the strata, the exposure of bedrock in borehole data is extremely limited. Multiple stratigraphic stratification points are rarely found within the same borehole. When using conventional modeling software, it is difficult to create effective constraints. Implicit modeling results in stratigraphic distortion and voids, and the unconformity between the bedrock and overlying sedimentary strata fails to generate correct spatial topological relationships. Furthermore, since only one geological profile is available for the target area, explicit modeling is also ineffective.

[0101] Using the construction surface control point generation algorithm proposed in this invention, based on Figure 9The modeling target area's orientation line generates equally spaced sampling points, such as... Figure 11 (a) Using step 2) in section 2.2, layer control points at different depths are obtained, such as... Figure 11 (b) After generating the layer data file, import it into the modeling software for rendering. The effect after zooming in on a specific layer is as follows. Figure 11 (c) The bottom surface of the loose layer in the target area of ​​the superimposed modeling region is used to perform TIN plane intersection line tracing and unconformity topology reconstruction, with the results shown below. Figure 11 (d) and Figure 11 As shown in (e), the final structural surface model of the target area after layer trimming is as follows: Figure 11 (f).

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

[0103] This invention utilizes strike control and attitude constraints to propose a batch generation algorithm for structural surface control points, and based on this, forms an algorithmic flow for creating 3D layers. This method can fully utilize the limited geological data of the modeling target area to quickly generate layer control point data, achieving structural modeling and solving the pain point of being unable to model due to insufficient data.

[0104] This invention, based on mature computer graphics fundamentals such as collision detection, implements an algorithm for reconstructing the spatial topological relationships of constructed surfaces, forming a complete process from data generation and storage to 3D model visualization. The method is simple, easy to program, and considers data format compatibility and extensibility, facilitating data sharing.

[0105] The application examples of this invention demonstrate that, relying on the occurrence information and structural features extracted from the planar geological map, the method presented in this paper effectively solves the problems of special steeply dipping strata structures in the target area, lack of basic data such as borehole data, and the inability of conventional methods to achieve structural modeling, thus verifying the feasibility and reliability of the method.

[0106] The above description is merely a preferred embodiment of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention. Structures, devices, and operating methods not specifically described or explained in this invention are implemented according to conventional methods in the art unless otherwise specified or limited.

Claims

1. A three-dimensional modeling method for structural surfaces based on orientation control and attitude constraints, characterized in that: Includes the following steps: Step 1: Extract stratigraphic strike lines, fault strike lines, and attitude data from the planar geological map of the target area. Extract the coordinates and azimuths of points on the strike lines using a curvature-constrained equal-interval sampling method. Combine the attitude data with batch generation of bedding control points by changing the sampling depth. Then, expand and constrain the bedding control points using stratigraphic or fault layering points revealed by boreholes or geological profiles to establish a set of structural surface control points, thereby creating an irregular triangular network of structural surfaces. Step 2: Use collision detection algorithm to extract structural surface intersections, reconstruct stratigraphic contact relationships and fault network spatial topology, and finally obtain the layer data file of the structural model, and complete the data storage and 3D visualization of the model; Step 2 specifically involves rapidly filtering intersecting triangles using axially aligned bounding boxes and multi-way tree collision detection. The Moeller-Trenbre algorithm is then used to find the intersections of edges and faces of overlapping triangles, extracting intersection lines and organizing the intersection point list using adjacency relationships. Furthermore, based on the master-slave relationship of stratigraphic contact or fault cutting, the control points are divided into three parts using bounding boxes of intersection points. Point affixing is determined by ray intersection, and one side of the points is retained to complete the layer trimming, simultaneously achieving spatial topology reconstruction of the stratigraphic contact and fault networks. Finally, the intersection points are added to the trimmed point set, and Delaunay triangulation is used for local TIN repair and to update the layer data file, thereby obtaining the layer data required for model construction.

2. The three-dimensional modeling method for structural surfaces based on orientation control and attitude constraints according to claim 1, characterized in that: Step 1 specifically includes the following steps: Step 101: Using GIS software, extract the strike lines of the strata from the planar geological map of the study area. The strata are stratigraphic interfaces or fault planes. Sample the strike lines at equal intervals and extract the coordinate set of the sampling points. ,in ; The initial elevation or depth value of the sampled layer object is represented by the elevation extracted from the DEM using the orientation line, and the corresponding initial depth is set in the depth domain modeling. This represents the azimuth angle of the directional line at that point, and the data table of the directional line is obtained by organizing the data. Simultaneously, curvature constraints are used to control the shape of sampling points and thin out the data. The principle is as follows: the equation of the B-spline curve passing through the sampling points is quickly obtained through data fitting, let it be... Starting from one end of the traverse line, use every three adjacent points to calculate the curvature, and iterate to the other end of the curve. The curvature at the point is calculated using formula (1); （1） This yields the curvature set of the directional lines. ; Step 102: Using the strike line as the top boundary line of the layer, construct right triangles in three-dimensional space for the sampling points along the strike line based on the layer dip and dip angle. Let the layer dip angle be... Sampling points In depth The projection point at that location is denoted as Based on the trigonometric relationships of plane projection, depth Place coordinates It is related to the dip and azimuth of the stratum; when the stratum dips southward, coordinates The results were obtained using formulas (2) and (3) respectively, when the strata dip northward. coordinates The results were obtained using formulas (4) and (5) respectively. （2） （3） （4） （5） Based on this, depth is obtained Set of control points at the ground level , express depth at The corresponding elevation, if modeled using the depth domain, is directly converted into depth. By controlling the sampling interval of the traverse lines and setting different sampling depths, the set of layer control points can be obtained. Therefore, by leveraging the constraints of bedding attitude and sampling depth, bedding control points can be generated rapidly in batches; simultaneously, the sampling depth... The positive and negative values ​​can be varied to flexibly adjust the sampling range of the layer control points; Step 103: Insertion and correction of borehole layer points; Extract the layer points from the borehole and insert them into the control point set created in step 102. The control point set at the level is encrypted and corrected.

3. The method for three-dimensional modeling of structural surfaces based on orientation control and attitude constraints according to claim 2, characterized in that: Step 1 further includes: Step 104: Store the layer control point data, where the Type field records the control point set. The types of strata represented are: Horizon for strata, Fault for faults, Name for strata name, and PointID for control point index number. Within the same stratum, PointID is a unique value. Record the coordinates of the control points respectively.

4. The three-dimensional modeling method for structural surfaces based on orientation control and attitude constraints according to claim 3, characterized in that: Step 1 also includes the following steps: Step 105: Create the layer triangulation 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 three-dimensional modeling method for structural surfaces based on orientation control and attitude constraints according to claim 4, characterized in that: Step 2, based on the establishment of the layer TIN in Step 1, constructs the contact relationships and fault network topology of the strata. The principle of the method is as follows: Step 201: Use a collision detection algorithm to filter intersecting triangles; use the axially aligned bounding box collision detection method to extract the intersection line of the TINs of two intersecting layers; construct a cube with edges parallel to the coordinate axes to completely enclose the triangle, and use a recursive method to construct a multi-branch tree structure to implement the collision detection algorithm. Step 202: Intersection Line Extraction: Since all layers are constructed using triangular meshes, after filtering using the axially aligned bounding box collision detection method, the extracted intersection lines are transformed into triangle intersection operations. When the bounding boxes of two triangles overlap, they are considered to intersect. The edges of one triangle are compared with the other triangle to obtain the intersection point. The Moeller-Trumble algorithm is used to calculate the intersection line, and the principle is as follows: Let's use the normal vector. and a point on the plane of the triangle Let the plane be defined, then any point on the plane is defined. Satisfying formula (6): （6） Use the starting point and direction vector Represents rays, Indicating the proportion of ray length, then the ray... Satisfying formula (7): （7） Let the three vertices of the triangle be... , , According to the rules of barycentric coordinate interpolation, the intersection of a ray and any point on the plane can be expressed as formula (8): （8） In the formula, , , Represents the centroid coordinates of any point on the plane of the triangle; , , Represents the coordinates of the three vertices of the triangle; only in equation (8) are the coordinates of the three vertices of the triangle. , , Three unknowns, when , , , When all values ​​are greater than or equal to 0, the ray intersects the triangle; Therefore, the vertices of the two triangles are traversed separately, requiring a total of six calculation operations for each pair of triangles. When the current triangle has an intersection point with another triangle, the intersection point and its adjacent intersection points are recorded; if there are no adjacent intersection points, the intersection point is the endpoint of the line of intersection. A linked list of intersection points is formed, and the intersection lines of the two levels can be obtained by tracing all intersection points in sequence. All intersection points are marked with a unique PointID. While saving the point coordinates, the Neighbors field is set to save the previous and next adjacent points of the current point. The endpoint has only one adjacent point. According to the tracing order of the intersection points, the first one recorded is the starting point, and the last one recorded is the endpoint of the intersection line.

6. The three-dimensional modeling method for structural surfaces based on orientation control and attitude constraints according to claim 5, characterized in that: Step 2 also includes the following steps: Step 203: Spatial topology reconstruction and layer clipping; Based on the TIN data created above; Constructing and clipping spatial topological relationships such as layer control point coordinates and triangle vertex indices, stratigraphic contact relationships, fault networks, etc., and transforming it into the problem of dividing the layer control point set by intersection lines and reconstructing the triangular network after division. First, establish the master-slave relationship between two surfaces based on the stratigraphic column or fault cutting relationship. Then, generate a bounding box using the set of intersection points. At this point, the secondary surface is the primary surface. When reconstructing stratigraphic contact relationships, if the strata are in an erosion relationship, the master surface is the unconformity, and the secondary surface is the older strata. When constructing the fault network, the order or priority of fault development is used to determine the primary and secondary strata. When a fault cuts through strata, the control point set with the fault surface as the master and the stratigraphic surface as the secondary will be divided into three parts: the inside of the bounding box and the two sides of the bounding box. Subsequently, based on the coordinate range, the control points of the secondary surface outside the bounding box are divided into two groups, A and B, and one point is randomly selected from them. Assuming a point is selected from group A... Control points of the surface inside the bounding box of the intersection line. Forming point pairs, the Moeller-Trumpé algorithm is used again to determine whether the connecting line PQ intersects with the principal surface; if an intersection exists, then... Points should be classified into group B, otherwise they should be classified into group A. After classifying the points, retain the data points from one side of the surface according to the erosion relationship to complete the surface trimming. By using intersection detection and layer clipping algorithms, the TIN file is traversed to make pairwise judgments, thereby completing the extraction of intersection lines and the reconstruction of the topology.

7. The three-dimensional modeling method for structural surfaces based on orientation control and attitude constraints according to claim 6, characterized in that: Step 2 further includes the following steps: Step 204: TIN repair; add the intersection points to the clipped set of surface points, use the Delaunay triangulation algorithm to locally repair the TIN, and update the layer data file.

Images