3D Modeling Method for Dome Structure Using Digital Geological Map and DEM

Geological maps have wide coverage with low acquisition difficulty. When other geological survey data are scarce, they are a valuable source of geological structure information for geological modeling. However, for structures with large deformation, geological map information has difficulty meeting the requirement of its 3D geological modeling. Therefore, this paper takes the dome structure as an example to explore a 3D modeling method based on geological maps, DEM and related geological knowledge. The method includes: (1) adaptively calculating the attitude of points on the stratigraphic boundaries; (2) inferring and generating the bottom boundary of the model from the attitude data of the boundary points; (3) generating the model interface constrained by Bézier curves based on the bottom boundary; (4) generating the top and bottom surfaces of the stratum; and (5) stitching each surface of the geological body to generate the final dome model. Case studies of the dome in Wulongshan in China and the Richat structure in Mauritania show that this method can build a solid model of the dome based only on geological maps and DEM data, whose morphological features are basically consistent with those embodied in the section view or the model generated by traditional methods.


Introduction
Geological phenomena and structures exist in three dimensions (3D). The solution of deep earth problems requires 3D visualization and geological analysis as the basis, and the premise of all this is 3D geological modeling [1,2]. Three dimensional (3D) geological modeling not only has very important application value in many fields such as urban planning, engineering construction, oil and gas storage, digital mines, etc., but also has a certain significance in the research of geological phenomenon interpretation, geological disaster prediction, geological environment assessment, and the quantitative simulation of geological effects [3,4]. The 3D geological modeling methods can be basically divided according to the data sources of modeling and the application of geological knowledge, as shown in Table 1. Mallet, 2002;Wu, et al., 2005;Hao, et al., 2021 [6,12,13] application mode of geological knowledge Explicit application In the modeling process, the knowledge of geology is explicitly used to construct the geological interface The geological data are limited, the geological structure type is clear and the modeling area is relatively simple  [19][20][21][22] According to the source of modeling data, the modeling method for geological bodies mainly includes borehole-based modeling, cross-section based modeling, geological-mapbased modeling and multi-source data modeling [2]. (1) Borehole-based modeling is a method of directly fitting and generating stratum solid models from borehole data [5]. This method is relatively mature and has a high degree of automation. However, due to the limitation of the number of boreholes, this method is mainly suitable for loose layer and key exploration mining area modeling; (2) Cross-section-based modeling is a method of generating a certain number of 2-D geological sections based on geological exploration data, and then generating a 3D geological model based on the constraints of the geological section [7][8][9]. However, the generation of the geological section requires a lot of survey work and manual editing by experts, which means large investment, so it is also suitable for loose layer and key exploration mining area modeling; (3) The modeling based on the geological map uses the information on the plane geological map drawn by the direct geological survey data to generate a 3D geological model. Because the plane geological map combines the results of geological field survey work and the knowledge of geological experts, it reveals the information of geological structures in the region. In addition, the stratiform structure of formations is constrained by surface boundaries and attitude. On the premise of a lack of other geological data, it is an effective solution for constructing 3D regional geological model by using geological maps [10]; (4) The modeling method with multi-source data integration is a method of fusing geological data from multiple sources such as boreholes and geological sections for 3D modeling [6,12,13]. This method has advantages in the construction of complex geological bodies, and has higher accuracy.
However, it is difficult to obtain different types of data with the same volume at the same time [23]. In addition, different methods of acquiring modeling data and data standards will inevitably lead to certain information conflicts, which will affect the accuracy of the model to a certain extent [24,25], so it is also suitable for loose layer and key exploration mining area modeling.
The expression and application of geological knowledge is a powerful supplement to geological modeling. It can improve the accuracy and automation of modeling by making reasonably constraints on the geological interface shape according to the characteristics of the geological body [26,27]. By increasing the application of geoscience knowledge in the process of geological modeling, on the one hand, when the data are so sparse that it is difficult to build a reasonable model, specific geological knowledge can reduce the dependence on modeling data to a certain extent; on the other hand, knowledge of geological experts can be integrated in the modeling process to improve the construction accuracy of the model [16,28,29]. In addition, different types of geological bodies have different genesis and tectonic forms, whose corresponding geoscience knowledges and modeling methods are also different. Therefore, it is an important basis for knowledgedriven modeling to effectively construct a geological knowledge base of a specific geological body type for knowledge application.
The application of knowledge in the modeling process can be divided into explicit application and implicit application. Explicit application uses an explicit definition of each object in the model and directly obtains the coordinates of key nodes, line segments and patches to form the surface of the strata [30]. More important, this method is based on geological semantic used to define the rules of geology, and these rules are used to describe different geological events and the correlation between these events, and then this correlation is used to restrict the construction of geological model [29]. In order words, this method explicitly considers the specific geological dimension of the model. For example, Perrin [16] proposed a knowledge-driven approach for SEM that shares throughout the workflow with the geological interpretation related to this model. Hou et al. [17] used the geological knowledge contained in the geological map to analyze the constraint relationship between strata and faults, and proposed a method to construct a 3D model of complex faults based on the planar geological map. In addition, Fernández et al. [14,15] and Vidal-Royo et al. [18] carried out research on obtaining dip domain based on geological mapping and outcrop interpretations, which simplifies geometries to volumes in which bedding attitude is constant, and proved the effectiveness and high accuracy of this method in experiments on Ainsa basin and Pico del Águila anticline, respectively. The method can well meet the resolution required of the reconstruction on the irregular geometry of the syncline or anticline. Implicit application uses an implicit definition of geological interfaces, which are defined as the iso-surfaces of one or several scalar fields in 3D space [30][31][32][33][34]. For example, De Kemp and Sprague [19] fits the stratigraphic boundary based on the Bézier curve, and generates structural zones based on the attitude data to upgrade the planar geological map. Laurent et al. [22] integrated the foliation data on fold structure in the numerical framework, which effectively improved the accuracy of fold structure modeling. When facing a modeling area with limited geological survey data for which the geological structure type is clear and relatively single, the explicit application of knowledge modeling method is more feasible; when facing a modeling area with relatively sufficient geological survey data which has relatively complex geological structure types, the implicit application of the knowledge modeling method is more feasible. However, when facing geological structures with complicated morphological ones, there are certain limitations in using explicit or implicit modeling based on geological maps alone [30]. Explicit modeling is prone to generating a considerable number of topological errors, and the degree of automation is difficult to improve, requiring a lot of manual interaction or repair measures. Further, the mathematical curve or surface model constructed by implicit modeling is too complicated to fit the actual shape of the geological interface.
The dome structure originates from the rock deformation caused by the movement inside the earth [35]. The structure is a kind of anticline with very significant non-cylindricity, whose rock layer is inclined outward and closed around and its geometric shape resembles an inverted grain bowl [36]. When eroded, the older strata in the core of the dome structure are exposed, which is roug y in the form of several strata including others on the geological map, that is, the overlying strata surround the underlying strata [20]. At the same time, the dome structure is a classic oil and gas mass closure, which finds it easy to accumulate oil and gas and mineral resources [36]. Therefore, the 3D modeling of dome structures has certain geological research significance and application value. Its geological interface is curved, and more survey data are needed to control the lateral shape. However, as a kind of bedrock structure, the dome also lacks sufficient survey data, so it is necessary to apply geoscience knowledge to reduce the data dependence of its 3D modeling. As the main geological survey results, the geological map has a wide coverage, a large amount of information, and contains considerable professional knowledge of geology, which is very suitable as a knowledge base for modeling applications. Therefore, for the regional geological structure such as the dome structure where the geological survey data are insufficient, the knowledge-driven modeling method based on the geological structure knowledge on the plane geological map can increase the feasibility of constructing a reasonable model [37]. Moreover, considering that the geological type of the modeling target in this paper is clear, and specific geoscience knowledge can be used to constrain the model, this paper mainly adopts explicit means in the application of geoscience knowledge. At the same time, in order to control the shape of the model more easily, knowledge is also implicitly used in some parts, such as generating side lines constrained by Bézier curves.
Further, for the dome structure strata with curved stratum boundaries and variable attitude, this paper mainly needs to use the knowledge of geology to solve the following three problems: (1) the attitude on the maps are scarce, which find it difficult to truly reflect the actual attitude of the stratum boundary; (2) separately deducing the stratigraphic interface from top edge points is likely to generate topology problems, such as self-intersection or burr of the stratigraphic interface; (3) simply deducing the stratigraphic side surface from the attitude of the stratigraphic boundary cannot reflect the true shape of the force and distortion in different parts of the stratum side.
The purpose of this paper is-taking the dome structure as an example based on the plane geological map and DEM-to study a method of using geoscience knowledge to make up for the lack of data and reasonably construct the 3D geological entity model in areas with sparse geological survey data. This paper will solve the problems above through an improved calculation method of attitude, the generation and optimization method of bottom boundary, and the refined modeling method on stratum side surface constrained by the Bézier curve. The organization of this paper is as follows: Section 2 introduces the methodology, Section 3 presents the experimental results, Section 4 presents the discussion, and Section 5 presents the conclusions and future work.

Methodology
This paper intends to study and realize an effective 3D modeling method of a dome structure based on DEM and digital geological maps. The modeling process mainly involves the following steps (   In reality, the dome structure may be affected by various geological processes and may have examples of various forms. Generally, affected by relatively single uplift-erosion event, the strata of the dome structure are in the form of several rings containing each other, which is the main form of this method. If the dome structure presents other forms due to later geological events, it needs to be restored to this form in a certain way. The situations that need to be dealt with are listed as follows: (1) if the fault causes serious displacement of the original geological structure, it is necessary to judge the feasibility of restoring the structure to the previous damage according to expert experience and restore it, and then use the method in this paper for modeling; (2) for the situation that magmatic rock or loose layer affects the complete exposure of bedrock, it is also necessary to be restored according to expert experience before modeling. Moreover, in order to evaluate the accuracy of the modeling cases in this paper, this paper selects the research area where the survey data can support the modeling of traditional modeling methods. Then, this paper constructs the model only based on sparse input data and compares it with the model constructed by traditional methods to analyze the accuracy.
The output of the solid model constructed in this paper is a surface model, that is, a closed space surface formed by TIN. A closed space surface represents a continuous stratum body. In order to construct this surface model, this paper firstly generates the exposed surface of the dome stratum, the upper and lower interface (that is, the inner and outer sides) and the artificially constrained bottom surface, and then stitches the surfaces to obtain the solid model of the dome structure.

Input Data
Strata: vector face data of the strata, in ESRI shapefile format, including the information of structural type and coordinates of all discrete points on the stratum boundary; DEM: digital elevation model, a kind of image in tif format, from which it can obtain contour lines, vector linear data of elevation contour, containing elevation information, in ESRI shapefile format; Attitude on map/measured attitude: the point data obtained after vectorization of the attitude elements marked on the plane geological map, in ESRI shapefile format, reflecting the attitude data measured at the marked position.

Method Parameters
Bottom elevation: the bottom elevation of the model set by the user is generally set according to the characteristics of the geological body. This value should not be set too small; otherwise, it will lead to low accuracy of the model. Furthermore, the maximum value should not be larger than the lowest part of the stratum boundary; Burr angle threshold: an angle value set by the user according to the characteristics of the modeling object. When identifying the burr of the bottom edge, if the included angle of the bottom edge is smaller than this value, it is necessary to further judge whether the included angle is a burr. Specifically, setting this value requires experts to comprehensively consider the structural deformation, erosion degree and geological survey data of the modeling object.

Attitude Calculation
The bottom boundary, which constrains the side and bottom of the stratum, is the key factor that determines the overall modeling effect. Additionally, the basis and prerequisite for the inference of the stratigraphic bottom boundary is sufficient information on the attitude of stratigraphic boundary. However, by only relying on the small amount of measured attitude information on the geological map, it is difficult to meet the high requirements for the quantity and accuracy of stratigraphic attitude information.
For this reason, based on the existing methods, this paper has developed an adaptive method for calculating the attitude of each point on the stratigraphic boundary. The main steps include ( Figure 2): (1) based on the five geological attitude rules under different conditions proposed by Xu [11], assign the measured attitude to corresponding stratigraphic boundary points; (2) extract the intersection points of the stratigraphic boundaries and the contour lines, and calculate the stratigraphic attitude according to the three-point or four-point method, and then insert intersection points on the stratigraphic boundaries; (3) calculate the attitudes of other boundary points based on linear interpolation, and correct abnormal attitudes.

Measured Attitude Conduction
The measured attitude points are distributed within the exposed range of each stratum in the plane geological map, indicating the attitude information of the stratum at the location. Ideally, each area with discontinuous attitude should have its attitude information. However, due to factors such as weathering of rock in the field, difficulty in distinguishing between cleavages and lineal planes, unexposed rock formations, constraints of geological survey fund and time, etc., the measurement of attitude is restricted.
According to the principle of ground stacking, if there is no major tectonic deformation in the later stage, excluding a certain degree of inclination, the interfaces of the sedimentary rock are basically parallel, and the measured attitudes of the outcropping parts are basically the same as those of the interfaces. Generally, both of the difference in the dip direction and angle between the several measured attitudes do not exceed 5° [38]. Further, the mean value of the measured attitudes can be conducted to to the whole boundary of the stratum. Furthermore, if several measured attitude data points of a series of mutually conformable contact strata are basically consistent, the mean value of the measured attitudes can be conducted to all the boundaries of these strata. If a certain degree of tectonic deformation has occurred, the measured attitude values of these strata may be quite different. However, due to the conformable contact relationship between the strata, there is still a certain correlation between the measured attitudes and the strata on the direction of the dip line. Thus, the measured attitudes can be conducted along the dip line to the intersections of the boundaries and the dip lines. Similarly, the measured attitudes on several strata that are parallel unconformable can also be conducted to the intersections of the boundaries and the dip lines. (Xu Feng, 2014). The specific method is as follows: firstly, between the strata with conformable contact relationship or parallel unconformable relationship, pass the point with measured attitude and draw a straight line along the dip direction; secondly, assign the measured attitude to the intersection of the straight line and each stratigraphic boundary; then, the attitude of point P1 can be assigned to the intersection points P2, P3, and P4 ( Figure 3). Besides stratigraphic attitude points, geological maps also contain attitude information of fault planes. A fault cuts a stratum and causes the stratum to break, forming a new stratum boundary, where the layer morphology is controlled by the fault and not transmitted to its neighbors, so the boundary here can directly give the fault attitude. As shown in Figure 4, there are faults passing through the formation, that is, two red lines in the figure. The attitude of the stratum boundary is the same as that of the fault, and the attitude of the fault interface can be directly assigned to the corresponding boundary points.

Adaptively Indirect Calculation of Attitude
Traditional indirect calculation methods of attitude commonly use the four-or threepoint method. Among them, the four-point method, namely the adjacent contour method, calculates the attitude of the strata from the four intersections of two contour lines with unequal elevations and the selected stratigraphic boundary [39,40]. The principle of the three-point method is to determine the shape of the plane by finding three points that are coplanar and non-collinear [41]. The four-point method uses more than three points, which can reduce the attitude calculation error caused by point error to a certain extent through adjustment, so that the accuracy of the four-point method is higher than that of the three-point method. However, the conditions of the four-point method are harsher than that of the three-point method. Therefore, in this paper, when several calculation points meet the calculation conditions of the four-point method and the three-point method, the four-point method is preferred to calculate the attitude, otherwise the three-point method should be considered. (Table 2). Both methods are suitable for calculating the attitude of stable rock stratum. Even for strata with severe tectonic deformation, the calculated attitude value can also reflect the overall attitude trend to some extent. Each small part of the stratum can be regarded as the plane with the same attitude. Therefore, based on several adjacent data points, three and four-point method can be used to carry out the adaptive calculation of the attitude of each part of the stratigraphic boundary in this paper. Therefore, this paper adopts the three-point method and the four-point method to carry out the adaptive calculation of the attitude of each part of the stratigraphic boundary. The specific steps are as follows: (1) extract the intersections of the stratigraphic boundaries and the contour lines as the calculation points; (2) take four consecutive calculation points on the same geological boundary, and judge whether the four-point method can be used to calculate attitude of the four points, and if so, calculate the attitude based on the four-point method, and go to step (4), otherwise go to step (3); (3) judge whether the three-point method can be used on the first three points of the four points, and if so, calculate; (4) judge whether there are four consecutive calculation points that are still not used for attitude calculation, if so, go to step (2). Table 2. Calculation conditions of three-point and four-points methods.

Three-Points Method
Four-Points Method

Recognition and Correction of Abnormal Attitude
As a kind of anticline structure, the dome is inclined and closed [36]. In other words, its dip direction around the boundary generally point outward relative to the stratigraphic boundary. However, the calculated attitude based on the measured attitude conduction, indirect method calculation, attitude interpolation or other related steps may have an inward dip direction, and this is an abnormal attitude for the dome structure. How to effectively correct these abnormal attitudes is the key issue to be solved in this section. This paper proposes a method for correcting the attitude of data points by boundary curvature. This method, which satisfies the law of general boundary attitude on dome structure, can make the dip direction toward the outside of the stratum according to the local morphology of the boundary. Therefore, for abnormal situations that face inward, the inclination can be corrected based on the curvature method ( Figure 5). The dip direction of the boundary point P j points to inside and is judged to be abnormal. At this time, based on the curvature method, the dip direction of P j can be changed to point from the center OP j of the circumcircle to P j .

Bottom Boundary Generation and Optimization
Near the stratum surfaces, the shapes of stratigraphic interfaces are constrained by the boundaries and their attitudes [38]. In this paper, each corresponding bottom point is calculated from each stratum boundary point, and then the bottom boundary is obtained by connecting each bottom boundary point in turn. Therefore, using the position and attitude of the point on the stratum boundary, the bottom boundary point can be reasonably generated.
However, due to the approximate elliptical shape and different attitudes of the dome structures, the generated initial bottom boundary often has problems of self-intersections and burrs. Thus, optimizing the initial bottom boundary, in other words, eliminating the self-intersection and burrs on the initial bottom boundary, is the key to reasonably generating the bottom boundary and the most important process that affects the 3D modeling accuracy of the dome structure. These two aspects are discussed in detail in the two subsections, respectively.

Treatment of Bottom Side Self-Intersecting
According to the characteristic that monotonic chains never intersect themselves, Yang [42] proposed a fast algorithm for judging the self-intersection of polylines based on the monotonic chain of computational geometry and an improved parallel line scanning algorithm. According to this algorithm, based on the monotonicity of the abscissa or ordinate of the adjacent boundary points ( Figure 6): firstly, divide the bottom edge into several monotonic chains; then, by pairwise judging whether the monotonic chains intersect, identify potential self-intersecting points; and finally, based on any self-intersecting point, divide the bottom edge into two parts, and then identify and eliminate the bottom edge "self-intersecting ring" according to the area of the enclosed polygon. Compared with the conventional method that finds the intersection of the line segments on bottom in pairs, this method has low complexity and high execution efficiency.

Treatment of Bottom Burrs
The bottom edge burr treated in this paper refers to the burr on the bottom edge calculated by a group of continuous data points on the top boundary without burr. It may cause the interface to form a long and narrow bending structure with sharp inward depression or outward protrusion in the longitudinal direction. The bottom burr is shown in Figure 7a. Generally, the unexposed conformable contact stratum layer will not form burr without experiencing severe structural deformation. However, when the difference of the dip direction between two adjacent boundary points on top is large, burrs may appear on the bottom edge. In addition, after the process of eliminating self-intersection on the bottom edge, burrs may also appear. The existence of burrs due to the bottom boundary generation will cause the generated stratum interface not be smooth enough, which will have a greater impact on the quality of the constructed 3D dome model. It should be pointed out that those burrs (thin valleys, canyons, etc.) produced by erosion generally appear on the exposed stratum erosion surface, on which this method will not be used.
For this reason, a curve point extraction method is proposed in this paper and is based on the included angle of adjacent segments to straighten the burr line segment and generate a smooth bottom boundary. The specific idea of this processing method is as follows (Figure 7): firstly, search for the boundary point with an angle less than a certain threshold on the bottom boundary, which can be judged as a "burr"; secondly, connect two adjacent boundary points of the current boundary point to fill the included the angle, and judge whether the burr is completely eliminated according to the included angle threshold. If not eliminated, repeat this operation until all burrs are eliminated.

Interface Generation under Bézier Constraint
Using a straight line connecting the corresponding points of the top and bottom boundaries, the formation side with straight side lines can be quickly generated. However, due to the different forces in the formation process, the side line of the dome structure is not a simple straight line determined by the boundary point and its attitude, but a smooth curve with a certain radian. Therefore, based on the corresponding boundary points of the top and bottom, a Bézier curve should be generated and discretized into several continuous line segments to obtain a relatively smooth side curve [43][44][45].
The steps based on Bézier constraints are mainly as follows: (1) generating control points of a Bézier curve, based on the points on the top boundary and the parameters of the Bézier curve set by the user with reference to the actual shape of the dome; (2) generating the Bézier based on the control points curve with parameterized expression, which will be discretized into t continuous line segments immediately by the user-specified number of pieces t; (3) connecting the i-th point on each side line in turn (i = 1, 2, . . . , t − 1) to obtain t − 1 closed curves; (4) Using the method in Section 2.2 to eliminate the self-intersection and burrs of the closed curve, and store the points on the corresponding curve in the side point set; (5) constructing a Delaunay triangulation network as the side surface of stratum based on the points of top and bottom boundaries and the side points.
The parameters can be set reasonably according to the borehole or section data of the experimental area in order to make the model of the dome structure more in line with the actual shape. See Section 4.1 for the specific method. The parameters that affect the generation of the Bézier curve are shown in Table 3. The bending coefficient of a control point is the ratio of the dip angle of the next control point to its dip angle, and the number of bending coefficients of the B é zier curve is n-1 t

Number of Pieces Discretize the Bézier side line into t segments
In the case of n = 3, correspondingly, there are four control points and two bending coefficients, and the generation process of the side curve is as follows: (1) Obtain an upper boundary point P 1 (x 1 , y 1 , z 1 ) of the stratum, that is, the first control point, and its attitude is α 1 (ρ, θ 1 ), where ρ 1 is the dip direction and θ 1 is the dip angle; (2) Calculate the difference ∆H between x 1 and bottom elevation H; (3) Calculate the position of the second control point P 2 (x 2 , y 2 , z 2 ) on the side line according to Equation (1), while the attitude of P 2 is α 2 (ρ, θ 2 ) = (ρ, m 1 × θ 1 ); Among them, ∆H is the difference between z 1 and the bottom elevation set by the user, and i is control point number, i ∈ [2,n+1]; (4) In the same way as in step (3), obtain the third control point P 3 (x 3 , y 3 , z 3 ) and the fourth control point P 4 (x 4 , y 4 , z 4 );

Bottom Generation
The modeling process of the stratum bottom surface can be boiled down to establishing a boundary constraint triangulation network based on inner and outer boundaries of the bottom [47]. Generally, in areas with small fluctuations, switching to a two-dimensional network construction method has less impact on the horizontal structure of the triangulation network, and because the dimension of the network construction process is lower, the computational complexity is much lower than that of the 3D method. Based on the above facts, this paper firstly extracts the inner and outer boundary points from a stratigraphic element and projects them on the horizontal plane. Secondly, it connects the boundary points in turn to obtain the inner and outer rings of the bottom polygon, and generates a two-dimensional bottom triangle based on the Delaunay rule. Finally, assign the elevation value of the bottom specified by the user to all the points of the triangle network to obtain the bottom model.

Top Generation
The top modeling process is similar to the bottom. Firstly, establish an initial triangulation based on the inner and outer boundaries of the strata; secondly, according to the sampling interval, extract sampling points from the contour lines within the exposure range of the corresponding strata, and inserts them into the initial triangulation under Delaunay's criteria; finally, give the elevation value of the top surface to all the points of the triangle network, and then obtains the top surface model.

Model Stitching
After generating the interface, top, and bottom surfaces of each layer in the dome structure, these surfaces need to be further stitched into solid models of each stratum. The essence of model stitching is to merge points with the same name on different faces of the same entity.
In the model stitching process (Table 4), if there is an underlying stratum s (i+1) in the stratum s i , the outer interface of s (i+1) needs to be regarded as the inner interface of s i . Then, based on the inner and outer interfaces, stitch top and bottom surfaces. Otherwise, only the outer interface and the top and bottom surfaces need to be stitched together. In addition, the attribute information of each stratum model is obtained from the corresponding stratum elements of the plane geological map, and stored in the corresponding JSON format model attribute file together with the model ID, which can effectively support the associated query between the stratum 3D model and attribute information.

Data and Experimental Platform
In this paper, the DEM and vector stratum and attitude point data of digital geological map are used for dome modeling. All algorithms were implemented using GDAL (Geospatial Data Abstraction Library, www.gdal.org, accessed on 28 March 2018) and compiled using the Microsoft Visual C# 2012 compiler, accessed on 26 November 2012.

Study Area and Dataset
This experimental area is the Wulongshan dome developed in the strata of Paleozoic in Zigui County, Hubei Province, China. The Wulongshan dome is located in the core area of the Yangtze Craton in South China, between the Huangling dome and the main body of the Xianglongshan anticline. It is a typical dome structure, which has continuous sedimentary strata and an abundant fossil record. This area has always been a hotspot for studying the frontier issues of geoscience, such as the growth and evolution of the Early Precambrian continental crust of the South China Yangtze Craton, and intracontinental compression, uplift and extensional tectonics of South China in the Mesozoic-Cenozoic era [48]. The structure has six strata from inside to outside, namely, Baota Formation, Nanjinguan Formation (including Honghuayuan Formation, Dawan Formation, Guniutan Formation, etc.), Longmaxi Formation of Ordovician, Shamao Formation, Luoraping Formation of Silurian, and Paomagang Formation of Cretaceous. In addition, the main structure is basically not covered by volcanic rocks and loose layers, which is conducive to the implementation of the modeling method in this article.
The source of the modeling data is the 1:200,000 scale geological map and the DEM data with 30-meter resolution (Figures 9 and 10). Both of them use the coordinate system of WGS84 and the projection of Web Mercator. There are two attitude points measured in the modeling area on the geological map and there are two faults, one of which corresponds to a part of the stratigraphic boundary. Therefore, it can be used for assigning the attitude of the corresponding boundary.

Modeling Performance
The internal stratigraphic contact relationship is parallel unconformity [49], so the attitude information of the geological map can be assigned to the corresponding strata boundary positions using the method proposed in Section 2.1.1 After calculating the attitudes of boundary points, the attitude information of each layer in the model is obtained and shown in Figure 11. In this case, on each side line, the number of pieces is set to 2 and the bending coefficient is 1.5. Without expressing the faults of this area, the generated model is shown in Figures 12 and 13.   Deng [49] drew a geological section of the Wulongshan dome based on seismic data (Figure 14a) after studying the structure of the Huangling dome and neighboring areas. The Xianglongshan anticline is mainly controlled by the nearly North-South compressive tectonic stress; that is, the stratum on the South side is formed by wedging the stratum on the North side of the anticline along the Silurian shale. The dip of the middle part and core part of the structure is small, while the dip of the two wings is large. The North wing of the anticline as a whole has a smaller uplift over the South wing, which suffers more denudation, especially where the F-F' section line passes through. The vertical section model (Figure 14c) is obtained by cutting the Wulongshan dome structure model generated in this case. By comparing it with the former section (Figure 14b) obtained in the same location, it turns out that their shapes are basically consistent. It indicates that the 3D model generated in this case can reflect the actual form of the dome structure to a certain extent. In addition, considering that the deeper the depth is, the lower the accuracy of the stratum level deduction will be, this experiment set the bottom elevation as 300 m.

Study Area and Dataset
The Richat structure is located in the Western part of the Sahara Desert in the Islamic Republic of Mauritania ( Figure 15). The structure is composed of several groups of concentric rings with a huge width. The entire area has a diameter of 40 km and an altitude of about 400 m. The current academic community believes that the multiple groups of prominent ring-shaped single-sided mountains in the Richard structure originate from the alkaline magmatic rock intruding into the rock plate during the Cretaceous period. Influenced by the different proportions of quartz in the rock formations, the alternating hard and soft rock layers at the top gradually rise, and then the dome structure is formed under the action of differential erosion [50]. From the core to the edge, the diagenetic age of the strata ranges from Late Proterozoic to Ordovician. The stratum dip is about 10-20°. Its central area is composed of siliceous breccia, covering an area of at least 30 km in diameter [51,52]. The geological map used in this experiment, based on the 1:1,000,000 geological map of Mauritania in the PRISM II project [53], is obtained by restoring its bedrock boundary after eliminating lava effluents and loose layers in the modeling area (Figure 15a). Further, the DEM data source is a digital product from ALOS PALSAR (search.asf.alaska.edu, accessed on 7 May 2022) with a resolution of 30 m (Figure 15b). Both of them use the coordinate system of WGS84 and the projection of Web Mercator.

Modeling Performance
In this case, the six strata of the Proterozoic era exposed in the core of the dome are selected for modeling, whose codes are At2, Ahb, Aht, Ahz, Ah, Te from the inside to the outside. Except for At2 formed in the Mesoproterozoic era, the rest of the strata all belong to Neoproterozoic, with continuous geological ages. However, since the original geological map did not mark the actual attitude information, the attitudes used in the modeling process were mainly calculated by the four-point and three-point methods. On each side line, the number of pieces is set to 5, and the bending coefficients are 1.2 and 0.8, respectively. The final model effect is shown in Figure 16. Matton and Jébrak [54], in the study of structural evolution and its genetic mechanism of Richat structure, based on high-resolution seismic data, exploration well location data and geological survey data, mapped the possible perspective interfaces of Richat structure, which showed the East-West internal structure (Figure 17a). The vertical section model of Richat structure obtained in this case (Figure 17b) is basically consistent with the model of Figure 17a, which proves that the model generated in this case can reflect the actual dome structure to a certain extent. In this paper, the adaptive three and four-points method is used to calculate the attitude of rock strata. In this method, the interface of fold is regarded as a curved surface composed of multiple micro planes, on which the calculated attitudes have a higher resolution. However, the attitude calculation method based on a few data points is more sensitive to local irregularity. Specifically, due to structural deformation or point location measurement error, there may be a large difference between the calculated attitude of adjacent points, which may lead to the need to deal with more spatial topology problems later. In the face of this local irregularity, Fernández [55] proposed two methods, planar region and moment of inertia analysis, to obtain average surface attitudes from points belonging to the surface, which have higher robustness in the face of local irregularity.

Effects of Different Side Line Bending Coefficients
In reality, the side line of the dome structure is not a straight line. According to the actual development, the shape of the side line may be convex or concave, or convex and concave at the same side. Furthermore, it is difficult to confirm the shape of these sidelines based only on the surface attitude data. Although the borehole and geological section data of dome structures are often missing or relatively sparse, and are not enough to be directly used as the basic data for modeling, these are the most objective data that reflect the regional geological structure and stratum distribution. Therefore, they are important reference data for reasonable determination of modeling parameters, and for quality verification and optimization of the model. Among them, the borehole data includes the location information of the borehole, the stacking sequence and thickness information of the drilled stratum, while the geological section data include the position information of the section and the stratum development form at the section line. When generating the side line of the Bézier curve, according to the borehole or section data, by setting the control point parameters reasonably and controlling the spatial position of several parts of the stratum interface, it can effectively correct and optimize the thickness and spatial distribution of the stratum, and improve the modeling accuracy of the side. Specifically, it is to determine the spatial position of several parts of the stratum interface according to the borehole and section data of the modeling area, and set appropriate Bézier curve parameters to control the generated stratum layer model close to the actual shape.
The comparative experimental results with different numbers of control points on each sideline and bending coefficients are presented Figure 18 below. By comparing the deviation between the generated and the measured side line, the fitting degree of the side line can be judged. In general, the more control points, the better they fit.

Treatment of Faults on the Dome Structure
In order to express the faults developed in the dome structure area, it is necessary to perform fault modeling based on the original dome structure model through the cutting and suture processing. The specific modeling method is as follows: firstly, based on the line data and attitude of fault, generate a slice-shaped 3D model; secondly, obtain the dome structure with its fault by the subtraction of 3D Boolean operations between the fault plane model and the dome model ( Figure 19). Since the stratigraphic boundary used in the dome modeling has already reflected the stratigraphic movement caused by fault, the 3D model established can also reflect it. In addition, when there are multiple faults in the modeling area, fault modeling should be carried out one by one according to the sequence of fault development, from old to new. For several intersecting faults, it should be judged whether they belong to the relationship of cutting through or cutting off, and then appropriately clip the fault models based on Boolean operation.

Interface Penetration Treatment between Strata
According to the spatial topology rule, in the absence of later structural transformation, on a series of conformable strata the inner side of a dome should not penetrate its outer side and the thickness of the strata is basically the same. Thus, fora series of conformable strata if the dip angle of the overlying stratum is greater than that of the underlying one, there may be a problem that the overlying stratum is penetrated by the underlying stratum.
This problem can be solved according to the principle of ground stacking. The specific method is to optimize the boundary line of the overlying stratum at the penetration part by extrapolating the boundary line of a certain thickness based on that of the underlying stratum, according to the principle of consistency of stratum thickness ( Figure 20). Because the outer stratum of the geological entity structure is often covered by plenty of loose layers, the boundary of the exposed stratum is not the real stratigraphic boundary. The boundary of loose layers may be misjudged as an abnormal attitude, which affects the calculation of the attitude there, thus reducing the quality of modeling. Therefore, in the actual modeling process, two treatment methods are needed. One is to avoid the use of loose layer boundaries as the input data of the three-point or four-point methods; the other is to use bedrock geological maps as the input data for modeling as much as possible. Compared with the regional geological map containing loose layers, the bedrock geological map is closer to the real shape of the bedrock structure. For example, in the Lingyanshan dome structure modeling in Nanjing (Figure 21), a regional geological map containing loose layers is used for modeling, which causes some attitudes of the loose layer boundary to be abnormal. When the bedrock geological map is used for modeling, there is no such misjudgment.

Model Accuracy Validation
The accuracy of the model is an important index to evaluate the quality of the model and the availability of modeling methods. The accuracy evaluation is generally carried out by comparing with various geological survey data in the modeling area, or models constructed by other modeling methods based on these data. However, this method is mainly used to improve the feasibility of modeling of area with sparse geological survey data. For this kind of modeling area, the geological survey data are difficult to meet the requirements of large-scale model accuracy evaluation, and the model constructed by datadriven modeling method also has certain limitations. Therefore, this paper puts forward two possible ideas: (1) calculate the deviation of the corresponding part of the model based on a small amount of geological survey data such as borehole and section. The smaller the deviation, the higher the accuracy of the model; (2) the rationality of the model can be evaluated by using geoscience knowledge, such as judging whether the side of the dome stratum extends outward and whether the thickness of the conformable stratum is basically the same. It should be pointed out that in the same modeling area, even if a large amount of geological knowledge is used to build a 3D model based on sparse data, its accuracy is generally difficult to exceed the model built with the support of rich data. The two validation ideas proposed in this paper can evaluate the accuracy of the constructed model to a certain extent, which can basically meet the needs of sparse data-based modeling.

Applicability of the Method
The dome structure modeling method based on a geological map in this paper is a modeling method that obtains and makes full use of the stratigraphic boundary and attitude, to derive the shape of stratum interface on the basic of considering the geometric characteristics of the dome structure. The data that this method rely on are fewer and more easily available. For common geological structures developed in bedrock, such as fault structures, fold structures, monoclinic structures, horizontal structures, etc., the stratum layer morphology has a certain correlation with the attitude on the surface, while the attitudes measured on these structures are few. For these geological structures, the geological map and the measured attitude combined with the calculation method of stratigraphic attitude proposed in this paper can generally meet the needs of layer morphology deduction. Therefore, the attitude calculation method and the 3D modeling method described in this paper have a certain applicability to this kind of modeling of geological structures.
The shape of the tectonic basin on the geological map is similar to that of the dome, which is roughly in the form of several interlocking strata, but its internal shape characteristics are opposite those of the dome structure. As a syncline structure with strong non-cylindricity, it resembles a placed grain bowl. The strata incline inward and the stratigraphic boundaries are closed. The geological interface within the geological boundary is a continuous curved surface [36]. In most syncline structures, the geological interface within the geological boundary is continuous curved surface. Therefore, after the downward deduction based on the attitude of various places on the surface, there is a problem of stitching the stratigraphic interfaces deduced by the strata boundaries. It is difficult to directly determine the specific shape of the inner interface. Especially at the turning point, it is difficult to directly determine whether the curvature is large or small and the specific form it presents based on the plane geological map. This requires a comprehensive analysis in combination with regional boreholes or cross-sections and the development characteristics of local folds. Therefore, the method in this paper is difficult to directly apply to the 3D modeling of basin structure.

Conclusions
This paper studies the method of constructing a 3D solid model of the dome structure using relevant geological knowledge based on DEM and geological maps, which mainly includes three sub methods. Among them, the attitude calculation method based on stratum boundary points can more accurately obtain the attitudes of the dome stratum boundary and can effectively identify and correct abnormal attitudes, which can solve the problem of sparse attitude and precisely control the 3D morphology of the stratum; the self-intersection and burr processing of the bottom boundary can effectively eliminate the self-intersecting surface and the narrow gap on the side surface of the stratum; and the modeling method of the side surface constrained by the Bézier curve can effectively control the morphology of the stratum side surface. This paper selects the Wulongshan dome and Richat structure as the experimental area. The experimental results show that the 3D model of the dome structure constructed in this paper using geological maps is basically in line with the real development form of the experimental area. Therefore, in the absence of borehole data and section data, the method of this paper can effectively model the dome structure. However, the method in this paper is susceptible to loose layer coverage and volcanic intrusion in the study area. Moreover, since the accuracy of the deduction of stratigraphic inference decreases as the depth increases, this method is mainly suitable for modeling shallow parts of geological bodies.
The method of this paper is helpful for research on the morphology and geological characteristics of dome structures and promotes the exploration and utilization of natural resources in the region where a dome structure is developed. In addition, the basic modeling idea of the 3D modeling method of a dome structure based on a geological map and DEM in this paper also has a certain reference value for the 3D modeling of other geological structures, such as fault structure, horizontal structure, fold structure and so on, based on sparse geological survey data.