The Method and Practice of Constructing 3D Geological Model from Coalfield Exploration 2D Maps ()
1. Introduction
The application of 3D geological modeling in mineral resource exploration and development is increasing, and its role is becoming increasingly significant. With the proposal and practice of the concept of intelligent mining and transparent mining, ore mining is transitioning from traditional mining mode to green, fine, digital, and transparent direction [1] [2] [3] . 3D geological modeling is an inevitable choice for coal exploration. The result of coalfield geological exploration is a geological report. Previous coalfield geological exploration work has accumulated a large number of coalfield geological reports. By the data of China National Geological Data Center, China has accumulated 20,789 coalfield exploration reports [4] . Among these, more than 4051 reports were submitted after 1990 [4] , the maps are stored in MapGIS/ AutoCAD format. According to the current coal geological exploration specifications in China [5] , 3D geological modeling is not mandatory requirement. So, in these coal exploration reports, the spatial expression form for the coal seams and their surrounding rocks are 2D maps. These electronic 2D maps are excellent data sources for constructing 3D geological models of coal field exploration areas. In the traditional geological reports, the spatial morphology and mining technical conditions of ore seams are based on 1) topographic and geological maps, 2) seam floor contour and reserve calculation maps, 3) seam thickness contour maps, 4) exploration line profiles, and 5) borehole columns maps. These maps all are two dimensional. The mine design, construction, and production departments need to invest a lot of manpower and resources for transporting the 2D maps to the 3D model, and then it is possble to add subsequent accessories for the intelligent mine. If the geologic prospecting party, as a service provider for ore mining, can provide both traditional geological reports and the digital 3D geological model, it will be welcomed by the mine owners. This will increase the added value of geological exploration results, broaden the service width and depth of the exploration team. How to build 3D model from these 2D maps has been studying for a long time, many methods have been proposed [6] - [14] . But, the 3D modeling for complex and relatively complex coal fields has not yet achieved a fundamental breakthrough. People need to accumulate little by little progress in order to achieve a leap from quantitative to qualitative change.
The author’s team has developed a software of “Coal Field Exploration 3D Modeling Platform” using VC++2010. Although there are still many functions that need to be supplemented and improved, it is basically available. This paper introduces the principles, experience, and insights used by the authors in the development of this software, with a focus on the principles and algorithms for constructing a 3D geological model using contour lines of seam floor. Semi-automation construction patch boundary ring and reverse fault duplicate area dealing method are originally created by the authors. The authors hope these works can promote the application of 3D geological modeling in coal exploration.
2. Geometry Expression of 3D Geological Surface
The geological surfaces are geometric abstraction of surfaces with geological significance, such as ground surfaces, stratum unconformity surfaces, rock/stratum interfaces, ore bed interfaces, fault surfaces, collapse column boundary interfaces, etc. The key to three-dimensional geological modeling in geological exploration is construction of geological surfaces. In metal ore and oil/gas exploration, the expression of non-uniform attribute fields within ore bodies has highly valued. In the coalfield exploration, people pay more attention to the spatial morphology and interrelationships of coal/rock layers, faults, collapse columns and so on. Coal is a typical layered mineral deposit structure, which is relatively stable horizontally compared with other mineral deposits. For 3D modeling of coal field, the bifurcation, merger, pinch out and tectonics expression of coal seams are the most important [15] . Therefore, the construction of geological surfaces plays a crucial role in the three-dimensional geological modeling of coal fields. In computer 3D graphics, geological surfaces belong to irregular free surface and are difficult to accurately express using an analytical expression. In the field of 3D geographic information systems and 3D geological modeling, there are three main modes to express geoscience surfaces [16] .
2.1. Contour Mode
Contour refers to the curve formed by adjacent points of equal height on a surface. Contour mode is a traditional form of geological surface expression. Before the emergence of computers, the geoscience community had always expressed geoscience surfaces using contours. The drawing method of expressing mountains with closed Contour line appeared on the Fangmatan wooden board map unearthed in Tianshui, Gansu Province at the end of the Warring States Period (239 BC) and on the silk Topographic map unearthed in Mawangdui, Changsha City at the Western Han Dynasty (186 BC) [17] [18] . However, this painting method has not been passed down and promoted. In 1728, Krucky of the Netherlands first used the isobath method to represent the depth and bed condition of rivers, and later applied it to represent the depth of the sea; In 1791, DuPontriel in France drew the first contour map; Subsequently, contour lines gradually began to be used in surveying topographic maps. The contour mapping method is a revolution in geological mapping. From then on, people can express the three-dimensional morphology of natural surfaces such as terrain, stratigraphic interfaces, and fault planes in a precise and measurable manner on a 2D map, as shown in Figure 1. Geological engineers can modify the shape of surfaces by modifying the shape of contour lines. So, it has unique advantages in coal exploration mapping. The contour drawing method also has some drawbacks, that is, when the surface is upright, the contour lines overlap with each other, and when the surface is inverted from positive to negative, the contour lines intersect with each other chaos. In a computer 3D drawing environment, although people can “see” the 3D shape of a surface by drawing 3D contour curves, it does not have the feeling of “surface continuity” and cannot achieve “color gradients” effect, as shown in Figure 2. Therefore, in the context of computer 3D geological modeling, the contour representation method ultimately gives way to the TIN method.
2.2. TIN (Triangulated Irregular Network) Mode
After the 1960s, computers began to be used for graphic display and rendering. People had found that the most effective way to display a 3D surface with a computer is to discretize the surface into many small triangle flat patches for display. The surface is represented by a collection of many small triangular flat patches embedded in an irregular triangular network. When displaying surfaces on a computer, setting different colors for the vertices of each small triangular
![]()
Figure 1. Contour line of coal seam 3 floor in a coal exploration area (2D view).
![]()
Figure 2. In a 3D viewport, contour lines can reflect the undulation of the coal seam floor surface (3D view).
patch based on their elevation can achieve the color gradien effect, making people’s “feeling” of the 3D shape of the surface more realistic. TIN has the advantages that neither contour nor regular rectangular grid has. Therefore, TIN is the standard expression of free surface in computer graphics, 3D geographic information system and 3D geological modeling.
2.3. DEM (Digital Elevation Matrix) Mode
In the field of surveying and geographic information, it is generally believed that DEM is the abbreviation of Digital Elevation Model. The digital elevation model is a digital simulation of geospatial surfaces such as the ground, by limited terrain elevation data (i.e. the digital expression of terrain surface morphology) [19] . According to the definition, digital elevation models should not only include regular square grid network, but also include contour and TIN. The essence of the grid network method is to represent the elevation of a surface in the form of an ordered numerical array. The DEM data provided by the surveying and mapping industry is a digital elevation matrix data of a certain resolution. So, DEM, as an abbreviation for Digital Elevation Matrix, is more appropriate. DEM data is a type of grid data in which data points are equidistant in the abscissa and ordinate directions. It is not necessary to record the horizontal coordinates of each grid point, but only one point horizontal coordinate of lower left corner (or upper left corner) and the grid spacing in abscissa and ordinate direction. It is easy to calculate the abscissa and ordinate coordinates of each grid point. So, for every grid point only needs to record its elevation. For data points without elevation values, use −9999 instead. Dividing each rectangular grid of DEM data into two isosceles triangles along the diagonal, and then it can be converted into regular TIN data. The regular TIN is a special case of irregular TIN. DEM is suitable for expressing geological surfaces with huge amounts of data such as terrain, but it is not suitable for expressing geological surfaces that may appear nearly upright or inverted.
3. Method of Fault Lines Obtaining Elevation
In 3D geological modeling by 2D map, it is required to assign elevation to all linear features such as contour lines of coal seam floor, fault lines, boundary lines of collapse columns, boundary lines of coal free areas, and boundary line of the work area, as shown in Figure 3. This means that these 2D curves need to be upgraded to 3D curves in order to participate in the surface construction of coal seam floor. The contour lines of the coal seam floor can be accurately assigned elevation value based on the digital marks next to the contour line. But, fault line, boundary line of collapse column, boundary line of coal free area, boundary line of work area, have no elevation in the map, and their elevations are different in different parts of the same line. How to correctly assign elevation values to them? Different 3D modeling software has different strategy. The strategy adopted by this software will be explained. Below, the fault line is taken as example to explain other 2D curve features than contour how to get elevation after contour lines having been assigned elevation.
3.1. Fitting a 3D Curve by Intersection Points of the Fault Line with Contour Lines
The fitting steps are: 1) obtaining the elevation from the endpoint of the contour line started or ended at the fault line, as shown in Figure 4; 2) Fit a 3D spline curve using these intersection points which obtained elevation from contour; 3) Assign elevations to other knots of the fault line based on the fitted 3D spline
![]()
Figure 3. Fault lines and workarea boundary lines need to been assigned elevation (2D view).
![]()
Figure 4. Fault line get elevation from contour line endpoint (2D view).
curve. The more endpoints the contour line intersects with the fault line, the more accurate the fitted fault line elevation is.
3.2. Fitting a Refer-Surface with the Contour Lines
In some cases, a fault line may not have been encountered any contour line, so another method should be used, which is to first triangulate based on the contour lines of the coal seam floor as the limiting line, and fit a rough reference surface—“refer-surface”. Since the fault line, the boundary line of the collapse column, and the boundary line of the coal free area are not used to build “refer-surface”, the fitted “refer-surface” is very rough. But, with this “refer-surface” as a reference, the elevation of each knot the fault line can be calculated by horizontal abscissa and ordinate coordinates of the knots.
3.3. Modifying the Elevation of the Fault Line on a Vertical 3D Profile
The elevation accuracy of fault line which obtained by the first two methods may not be satisfactory and require editing, as shown in Figure 5. It can be modified in 3D environment, by generating a vertical profile that curves along the fault lines, then move the knots in the vertical profile, where knots only can move in vertical direction, can not move in horizontal direction, as shown in Figure 6.
3.4. Modifying the Elevation of the Fault Line on a 2D Section
Although it is possible to modify the elevation of a fault line in the 3D environment, its operation is not convenient enough, and people are still more accustomed to editing 3D curves in a 2D environment. So, another method is to project the 3D curved vertical profile onto a 2D vertical section. It is a reverse operation to translate 2D profile to 3D profile [11] . Modifying the elevation of knots of the fault line on the 2D section and then mapping the modifying into the 3D profile.
![]()
Figure 5. Fault lines and workarea boundary lines have been gotten elevation from contour line endpoint (3D view).
![]()
Figure 6. Along fault line generate a vertical profile (3D view).
By comprehensively applying the above four methods, the fault line can be obtained basically satisfactory elevation accuracy. The accuracy of 3D modeling requires continuous iteration to approach accuracy. After the following semi-automatic partition triangulation, the step 2) “refer-surface” should be replaced by the partition triangulation surface. Then recalculating the elevation of each node on the intersection line of the sections will improve the modeling accuracy has been improved.
4. Semi-Automatic Partitioning Triangulation
A basic principle for expressing surfaces through irregular triangulation is that spatially separated surface patches must be triangulated independently. The surface patches that are divided by boundary elements such as the fault line, the boundary line of the collapsed column, and the boundary line of the coal free area, must be triangulated independently. If the coal seam floor surface is completely cut off by a fault line, it must be divided into independent patches. If the fault line does not completely cut off the surface, it should be treated as a limited ring within the patch. Fully automatic triangulation is a goal pursued by people, but its implementation is quite difficult. This software take a method by user click to select line features to construct a surface patch boundary ring, and then automatically processes various line features that fall into the partition boundary ring, thereby achieving semi-automatic partition triangulation [20] . The specific steps are as follows.
4.1. Semi-Automatic Construction of Partition Boundary Ring
Users need to determine in advance how many partitions the entire area needs to be divided into based on the cutting situation of faults, collapsed column boundary lines and coal free boundary lines. Each partition is a surface patch, which is boundared by a ring consisted of fault line, collaped column boundary line and coal free zone boundary line. Then, for each partition, the user clicks to select the line features (including fault line, collapsed column boundary line, coal free boundary line, and the boundary line of the work area) consisted of the partition boundary ring. When clicking in anticlockwise order, the software will automatically cut and connect the line feature, and finally form the boundary ring of the surface patch. When clicked, select one layer at a time as an editable layer (note: not to edit or modify this layer here, but only to set the layer in the editing state as the selected layer), forming a closed boundary ring, as shown in Figure 7.
After the partition boundary ring is closed, the line features that fall into the partition boundary ring should be treated as limiting lines (as shown in Figure 8) or limiting rings in triagulation procedure.
4.2. The Construction of a Fault Ring
The software will scans and identifies whether the fault line ring formed by the fault lines that fallen into the partition boundary ring is a normal fault ring or a reverse fault ring, and then treated separately.
A fault ring is formed by closing the fault lines on the upper and lower walls of a same fault. The algorithm for generating a fault ring from the fault line is as follows: 1) Scan all fault lines, record all the fault lines that fall into the boundary ring of the patch, and generate a set of fault lines to be matched; 2) Check the set of the fault lines one by one, match each fault line with all other fault line in the set. If one endpoint meets, check the other endpoint; If the other end also coincides, the intersection line connecting these two sections in a clockwise direction forms a fault ring; If the two endpoints on the other end do not coincide, check whether the other end of the intersection line of the two sections intersects with the patch boundary ring. If so, connecting the intersection line of these two sections clockwise with the intersecting boundary segment to form a fault ring. 3) For each fault ring formed, the two fault lines are removed from the set of intersection lines.
![]()
Figure 7. Click sequence of new patch boundary ring (2D view).
![]()
Figure 8. New surface patch 4 triangulation (2D view).
4.3. Determination a Boundary Ring Is Normal Fault Ring or Reverse Fault Ring
Based on whether exist contour lines of the coal seam floor between the fault rings, it can be determinated that the fault ring is a normal fault ring or a reverse fault ring. If exist one or more contour line in the fault ring, it indicates that the fault ring area is a duplicate area of coal seams, and the fault is a reverse fault (Figure 9, left); If not, it is a normal fault (Figure 9, right).
4.4. Algorithm for Processing Normal Fault Rings
The interior of the normal fault ring is a blank area with missing coal seams, so triangulation is relatively simple. The processing of normal fault ring is to organize the sampling point of the normal fault ring in a clockwise direction as a inner boundary ring in triangulation. As shown in Figure 10.
The treatment method for boundary rings of collapsed columns and coal free areas is the same as that for normal fault rings.
The processing algorithm of reverse fault rings is very complex, it will be discussed detailly in next Chapter 5.
![]()
Figure 9. Reverse fault (left) and normal fault (right) (2D view).
![]()
Figure 10. Normal fault ring in the partition is used as triangulation limit ring clockwise (2D view).
5. Algorithm for Processing Reverse Fault Rings
Triangulation of coal seam floor surface cut by reverse faults has been bothering software developers of 3D geological modeling for a long time. The difficulty is that exist a duplicate area of coal seams. The author proposes a method for process the reverse fault ring duplicate area [21] . It can be roughly summarized as next four steps.
5.1. Generate an Outer Expanding Ring of Reverse Fault Ring
Firstly, an outer expanding ring is generated by off reverse fault ring a distance, such as 10 meters. The procedure is shown as Figure 11.
5.2. Dividing “Duplicate Area” from the Patch
The outer expanding ring separates the reverse fault area from other part of the patch which is normal area without the duplicate area. For the normal area, it can be triangulation normally, as shown in Figure 12.
5.3. Duplicate Area Being Divided Tow Sub-Patches and Triangulation Separately
The duplicate area will be divided to two side, left side and right side. The two sides left side and right side all are contain diclicasd area. Left side and right side can be triangulated severally. As shown in Figure 13, Figure 14.
![]()
Figure 11. Reverse fault ring generate outer encloure ring (2D view). (a) fault line 96,97 (b) fault lines form a reverse fault ring (c) expanding from reverse ring to form outer enclosure ring.
![]()
Figure 12. Normal area is triangulation normally (2D view).
![]()
Figure 13. (a) the left side of the outer ring and right side of the reverse fault ring combine form a left ring, (b) the left side of the reverse fault ring and the right side of the outer ring combine to form a right ring, (c) triangulation left ring, (d) triangulation right ring, (e) put left ring TIN and right ring TIN together (2D view).
5.4. Removing Expanding Ring and Suturing the Duplicate Area with Mather Sub-Patch
Remove the outer ring of the duplicate area and non original data with topological adjacency, locally reconstruct the triangular network, and generate the final set of research area triangular networks. As shown in Figure 15
![]()
Figure 14. Expanding ring connecting normal area TIN and overlap area TIN (2D view).
![]()
Figure 15. After remove expanding ring (note: among fault ring, the left TIN and right TIN are overlapped each other) (2D view).
6. Extraction Method for Ridge and Valley Lines
Ridge lines and valley lines are invisible line elements that can be extracted from contour lines, similar to the lines of ridges and valleys in topographic maps on a surface. In the geographic information community, ridge lines and valley lines are collectively referred to as geographic lines. When constructing triangulation on a surface, there is a principle that the edges of the triangle should not cross ridges and valleys. If the ridge and valley lines can be automatically extracted after triangulation, and the ridge and valley lines are also used as the limiting lines for triangulation, the quality of the triangulation network will be improved. Our team has studied a center of gravity method based on constrained Delaunay triangulation to extract mountaintop or valley lines. The steps are as follows.
6.1. Preliminary Delaunay Triangulation
Using the contour line’s own range as the outer boundary, Delaunay triangulation is carried out internally. Observing the results of the triangulation, the triangulated triangles can be divided into three types: Type I, where two edges in the triangle are constraint edges, that is, two edges are on the contour line; Type II, in a triangle, one edge is a constraint edge, that is, only one edge is on the contour line; Class III, unconstrained edges in a triangle, meaning that all three edges of the triangle are not on the contour line, as shown in Figure 16.
6.2. The Center of Gravity Method Is Used to Extract Geographical Lines
In the constructed Delaunay triangulation network, for type I triangles, the midpoint of the unconstrained edge in the triangle is connected to the triangle vertex relative to the unconstrained edge; For type II triangles, the midpoint of the two unconstrained edges of the triangle is connected; For type III triangles, the midpoint of the three sides of the triangle and the center of gravity of the triangle are connected. Finally, connect these connected line segments according to a certain topological relationship to form a whole, which is the desired central axis, as shown in Figure 17.
![]()
Figure 16. Three types of triangles, the solid line is a segment of Contour line (2D view). (a) Type I triangle; (b) Type II triangle; (c) Type III triangle.
![]()
Figure 17. The center axis acquired by gravity method (2D view).
6.3. Add Geolines to Reconstruct 3D Surfaces
Add the extracted geolines in Figure 17 as constraints and reconstruct the triangulation network, as shown in Figure 18. Extracting geolines and using them as constraint lines will make the surface triangulation network more optimized.
7. Construction Fault Surface from Fault Lines
Fault surface reconstruction is an important part of three-dimensional geological modeling. When there are multiple coal seam floor contour maps in the exploration area from top to bottom, the intersection lines of the same fault plane and different coal seam floor contour lines can be used as the skeleton line. Triangulation can be carried out between the skeleton lines to construct the fault plane and extend it downward and upward to a certain height.
Figure 19 and Figure 20 are examples of constructing fault planes by intersecting the upper and lower coal seams of a certain coal mine in Shanxi.
![]()
Figure 18. The triangles after reconstructing with feature points from terrainlie in it (2D view).
![]()
Figure 19. Preparing construct fault surface by cross section lines of different upper and lower coal seams (3D view).
8. Construction Stratum and Coal Seam Body
8.1. Construction Stratum Body
The construction method of stratigraphic body is to use the boundary rings of the top interface and bottom interfaces of the stratigraphic body for constrained triangulation, forming a triangular network of the lateral torus of the stratigraphic body, as shown in Figure 21. After coloring the lateral torus, it is shown in Figure 22. After closed top, bottom, and side torus surfaces, it is shown in Figure 23.
8.2. Construction Coal Seam Body
Generally, there is no contour map of coal seam roof in the geological report of coal field exploration, but a contour map of coal seam thickness is provided, and there is basically no fault information on the contour map of coal seam thickness. Therefore, for the construction of coal seam body, the zoning boundary of the
![]()
Figure 20. Fault surface constructed by fault curve of upper and lower coal seams (3D view).
![]()
Figure 21. Triangulation of boundary rings on the top and bottom surface of the stratum to construct the lateral torus of the stratum (3D view).
![]()
Figure 22. Triangle coloring of lateral torus (3D view).
![]()
Figure 23. Top, bottom and side Annulus Sealing Formation Body (3D view).
coal seam floor surface can be used to overlay the coal seam thickness as the coal seam roof surface, and then the construction method is the same as that of the stratum body. The bifurcation, merger and pinch out expression of coal seams are the most important in the construction of coal seams 3D modelling [14] .
9. Conclusions
This paper research methods of construction 3D geological model from 2D coal maps. The paper gives new ideals in the next problems:
1) There are three modes to express a geological surface: contour, TIN and DEM. Contour is the most suit for geologist modify surface shape; but TIN is the most suit for geological surface 3D modeling; DEM as Digital Elevation Matrix, is not suit for expressing geological surface.
2) Fault line get height by meet point with contour and refer surface and can be edited in 3D Profile and 2D section.
3) Semi-automatic partition method by clicking select boundary line features in counter-clockwise direction sequence way is proposed.
4) Automatic method to deal with fault ring fallen partitioning boundary ring is proposed. It can automatic recognize normal fault ring and reverse fault ring.
5) A method of automatically processing reverse fault ring is proposed. By adding an outer expanding ring, duplicate area will be separated with normal area of the patch, then dividing the duplicate area as two sub-patches which are single layer.
6) Algorithm for extracting ridge and valley line is proposed.
7) Construction fault surface from different coal layer fault lines is proposed.
8) Construction method for stratum body and coal body is given. The coal top surface is made by coal bottum surface + coal thickness.
NOTES
*First author.
#Corresponding author.