Research on Automatic Construction Method of Three-Dimensional Complex Fault Model

: Three-dimensional complex fault modeling is an important research topic in three-dimensional geological structure modeling. The automatic construction of complex fault models has research signiﬁcance and application value for basic geological theories, as well as engineering ﬁelds such as geological engineering, resource exploration, and digital mines. Complex fault structures, especially complex fault networks with multilevel branches, still require a large amount of manual participation in the characterization of fault transfer relationships. This paper proposes an automatic construction method for a three-dimensional complex fault model, including the generation and optimization of fault surfaces, automatic determination of the contact relationship between fault surfaces, and recording of the model. This method realizes the automatic construction of a three-dimensional complex fault model, reduces the manual interaction in model construction, improves the automation of fault model construction, and saves manual modeling time.


Introduction
Fault modeling is the basis of three-dimensional (3D) geological modeling and has important value in basic geological research and various geological science and technology applications. In areas with developed faults, there are often not only a large number of faults, but also complex faults with special structures, such as X-type faults, multilevel y-type faults, etc. It is very difficult to model faults in such areas. Currently, there are two main ways to construct fault models: the explicit way and the implicit way. Explicit modeling is constructed by geologists according to their own knowledge and experience in order to adjust the model through manual interactive operations, and the result is subject to the cognition of a modeler [1][2][3]. Implicit modeling adopts the mathematical solution method of scalar fields, but it requires the setting of boundary conditions and additional artificial data points to make the scheme feasible [2,4,5]. Godefroy et al. [6] and Grose et al. [7] introduced kinematics into implicit modeling and realized the parametric modeling of folds. Current implicit methods can only represent a limited number of unconformities owing to computational and storage limitations [3]. Explicit modeling can be flexible, but it is time-consuming due to numerous interactive operations. To free geologists from the complex interactive modeling work and promote the application of 3D geological modeling in related fields, the automatic construction of 3D complex fault models has become a major issue.
The construction of the fault model mainly includes four steps: data preparation, fault surface generation, fault-fault contact relationship processing, and model recording. In addition to the geological survey data, the data required for modeling often uses seismic data, logging data, and remote sensing images. These data are geologically interpreted to obtain geometric data such as points and lines that meet the requirements as input data for modeling [8,9]. Jessell et al. [10] realized the automatic deconstruction of geological maps and the automatic definition of fault contact relations; Wu, X. et al. [11,12] and Qi et al. [13] realized automatic fault interpretation of seismic data. As for the 3D fault surface, automatic generation and visualization technology has been relatively mature. The current mainstream method uses a triangular mesh model to finely describe the detailed characteristics of the fault surface [14,15]. The pillar grid and stair-step grid can also describe fault surfaces, but they will simplify the surface characteristics and reduce the model accuracy to improve the efficiency of numerical simulation [16][17][18][19][20]. Exploring and realizing the automatic processing of the fault-fault contact relationship are key issues for the automatic construction of 3D complex fault models. The processing of the faultfault contact relationship includes the determination of the intersection station, the major or minor fault, the termination of fault surfaces, and the truncation of the minor fault. Although the judgment of intersection has been automated, the determination of the major or minor fault, the processing of fault surface termination, and the truncation of minor faults still require manual interaction [21][22][23]. The determined result of the fault-fault contact relationship is mainly recorded by the key pillar, matrix, and binary tree methods. The pillar-based fault model uses a special pillar that is manually added to mark the intersection of fault surfaces, but it cannot handle multilevel branch faults. The matrix method can completely record the contact relationship of each fault, but it is not easy to edit intuitively with high information redundancy and high dependence on manual participation [24]. The fusion fault block method extends fault surfaces to divide the geological body into fault blocks and organizes the fault-fault contact relationships in the form of a binary tree, which can construct multilevel y-shaped faults [25,26]. However, the truncation operation still requires manual participation, and the accuracy of the model needs to be improved [27].
This study provides a new solution for the automatic processing of the fault-fault contact relationship of the complex fault model. Multisource data (seismic data, logging data, remote sensing images, etc.) were interpreted as labeled fault point data, and the knowledge of geologists was converted into point data that can be easily recognized by the computer through labels. With the labeled datasets, the modeler can automatically construct the same complex fault truncation model. Therefore, the shortcomings of explicit modeling, which are not easy to reproduce, can be fixed. The specific modeling process will be described in Section 2. The methods of model construction, minor fault truncation, and model recording are described in Section 3. After constructing a fine triangulation fault surface model with the fault point data, the intersecting relationship of the fault surfaces is automatically judged and recorded by the graph model of graph theory. Graph theory has some applications in the field of geology (such as fracture modeling [28][29][30], two-dimensional (2D) geological map recognition [10], the organization of geological relations [31,32], etc.). Compared with the fault model, the fracture model has a small scale and a simple surface structure. The application of the graph model in these two problems is also different. In the fracture network, the edges of the graph model represent the fracture surfaces, and the vertices represent the intersection points of the fracture surfaces. For large-scale objects such as faults, graph theory has more applications. For instance, in the recognition of 2D geological maps, the vertices and edges correspond to geometric points and edges [10,30]. In uncertainty modeling, the vertices represent the fault interpretation points, and the edges represent that the two fault points belong to the same fault [33]. In this study, each vertex represents a fault surface, and the edges represent the intersection relationship of the fault surfaces. The major and minor type of the intersecting faults are automatically updated according to the attributes of graph theory. Further, the minor faults are automatically truncated according to the truncation standard. The final result of the fault model is recorded in extensive markup language (XML), which not only makes it easy to analyze the structure of the fault-fault contact relationship, but also makes it easy to modify the type of major or minor fault. In short, the method in this study reduces the manual interaction in model construction, improves the automation of fault model construction, and saves manual modeling time.

Automatic Construction Process of Complex Fault Model
The automatic construction of a 3D complex fault model includes three steps: the construction of the 3D fault surface model, the calculation of the fault contact relationship, and the description of the complex fault model (Figure 1). The contact relationship of faults is automatically calculated, the truncated fault model is automatically generated, and the results of the fault model construction can be reused after recording.

Data Preparation
Multisource data can be interpreted into fault point data, fault sticks, and fault polygons. If there are enough fault point data, fault sticks and fault polygons can be transported to a pointset. Therefore, the fault data in this paper were interpreted as labeled fault point data. Fault points belonging to the same fault have the same label. To avoid the misunderstanding of truncation, fault point data on two sides of one fault will be marked with different labels.

3D Fault Surface Model Construction
The triangulation method of the two-dimensional plane is very mature [12], so the 3D fault point data are first projected to the two-dimensional (2D) plane (a least squares fitting plane illustrated by Figure 2) one by one and triangulated to generate a 2D triangulation plane. Then, according to the mapping rules, the two-dimensional triangulation model is mapped onto a 3D triangulated surface mesh, and mesh optimization is performed to finely describe the fault surface shape.

Determination of Fault-Fault Contact Relationship
The contact relationship of fault surfaces includes disjoint ( Figure 3a) and intersecting relationships, in which an intersecting relationship contains intersection stopping ( Figure 3b) and intersection cutting ( Figure 3c). To determine the fault-fault contact relationships, it is necessary to determine whether the faults intersect first. There is a threshold about fault-fault distance between two disjoint faults. If the distance is within the threshold, one fault's surface will be extended to make an intersection as shown in Figure 3b [34], and the contact relationship is deemed to be intersection stopping. According to recorded intersection status, the intersection line of the two intersecting faults is calculated. The location of the intersection line on the fault surface will decide whether to extend the intersection line to truncate one of the fault surfaces (details in Section 3.2).

Recording of Complex Fault Model
Graph theory was utilized in this paper to record complex fault models. The vertices represent the fault surface models, and the edges represent the contact relationship of two fault models. If there is an edge between two vertices, it means that the fault surfaces represented by these two vertices intersect. The XML can flexibly record the graph structure [35], and it is convenient to give abundant information to the elements in the graph. The XML can also describe the complex fault model and record its information, which includes details regarding the fault surface and the fault contact relationship. The output of the model result will be a triangular mesh and graph network about fault-fault contact relationships.

Automatic Construction Method of Complex Fault Model
At present, the construction of complex fault models, especially when calculating the complex contact relationships of faults, requires significant manual participation. The key point in realizing the automatic construction of complex fault models is to allow the computer to automatically calculate the fault contact relationship and determine the truncation form. This paper proposes a set of complex fault model construction methods based on triangulation network fault surfaces, including (1) the automatic generation method of fault surfaces, (2) the automatic calculation method of fault contact relationships, and (3) complex fault model recording methods.

Automatic Generation Method of Fault Surface
The automatic construction process of the 3D fault surface model is as follows ( Figure 4):

1.
Calculate the centroid coordinates and 2D fitting plane normal vector (w 1 , w 2 , w 3 ) of the 3D fault point dataset. The plane equation (S : w 1 x + w 2 y + w 3 z + w 4 = 0) can be calculated according to the centroid coordinates and normal vector. Project the 3D points' set to the 2D surface one by one and perform triangulation to generate a 2D triangular mesh surface.

2.
Based on the one-to-one correspondence between the points in the 3D space and the 2D space, the topological relationship of the two-dimensional triangulation model is mapped to the 3D space to form a 3D fault surface.

3.
Optimize the fault surface model, reduce the deformed triangular grid unit, and ensure the quality of the grid.

Automatic Determination Method of Fault-Fault Contact Relationship
The two intersecting faults are called major and minor faults, respectively. The major fault will not be truncated, whereas if the minor fault needs to be truncated, further judgment is required [21]. The automatic determination of the fault-fault contact relationship first needs to determine whether the fault surfaces intersect. Two faults that have a distance within the threshold, as discussed in Section 2.2, will also be regarded as intersecting surfaces. After determination of the intersection, a complex fault model is generated according to the following process (   The first basis for determining whether the fault is a major or minor fault is the formation time of the fault, and the major fault has a later formation time than the minor fault does. Abbott et al. used geometric elements as the basis for judging the major and minor faults, such as the number of fault points, the degree of extension of the fault surface, and the number of fault points on both sides of the intersection line of the fault surface [21]. In the simple case where two faults intersect, geometric elements are feasible as a basis, but, for multilevel complex fault models, this method cannot be used to make effective judgments. In addition, the number of faults associated with a certain fault should also be considered to ensure that the faults in the model are as interconnected as possible.

Calculation of Intersection & Extension Line of Fault Surface
The intersection line of the intersecting fault surface is formed by the intersection point produced by the intersection of the triangular patches. The intersection line of the two intersecting fault surfaces was calculated according to the following steps: 1.
Record the intersecting triangle patches in each fault surface.

2.
Select one of the fault surfaces and calculate the intersection points on each triangle plane recorded in Step 1. One triangular patch may intersect with multiple triangular patches, resulting in multiple intersection points ( Figure 6).

3.
Based on the topological relationship between the triangular patches, all the intersection points on the fault surface are related to generate the intersection line.    1 and F 2 are two intersecting fault surfaces in space, and it is known that F 2 is the major fault and F 1 is the minor fault. S 1 is a reference surface, which is actually the xOy plane paralleled with the horizontal plane, and S 2 is the parallel plane that fits the 2D plane of F 2 . F 1 and F 2 intersect to produce intersection points to form a line segment L = {J 1 , J 2 , J 3 , J 4 , J 5 }, which cannot divide F 1 into two parts. Therefore, the intersection line needs to be extended to divide F 1 ; since the process of extension is the same on both sides of the line, the following example explains the process in detail for only one side.
Find the triangular patch t 1 that contains the point J 1 (x 1 , y 1 , z 1 ) in F 1 according to the geometric topological relationship and add the vertices into the point dataset of F 2 . Calculate the new fit plane F 1 2 :w 1 x + w 2 y + w 3 z − (w 1 x 1 + w 2 y 1 + w 3 z 1 ) = 0 based on the normal vector (w 1 , w 2 , w 3 ) of F 2 with the new point dataset and the point J 1 (x 1 , y 1 , z 1 ). The intersection point J 1 will be formed by the triangular patch t 1 and F 1 2 . J 1 is a point on the extending line of the intersection. Add another three vertices of the triangular patch t 2 into the point set of F 2 and fit a new plane F 2 2 . Calculate the normal vector (nw 1 , nw 2 , nw 3 ) of the new plane, F 2 2 , which can form a new plane with point J 1 . The triangular patch contains J 1 , and the new plane can create a new intersection point J 2 , . . . . . . , and so on. Finally, the extended line is traced and generated until the edge of F 1 . Figure 8 shows that F 1 is truncated into two parts by F 2 using the extension line segment set L = {J 1 , J 2 , J 3 } and intersection line segment set L = {J 1 , J 2 , J 3 , J 4 , J 5 }. According to the rules of truncation described in Section 2.2, the right part pointed by the arrow will be removed, and the left part will be kept.

Updating Geometrical Topological Relations with Ear Clipping Algorithm
The fault surface is represented by a triangulated irregular net (TIN). However, the calculated intersection points of the two intersecting fault surfaces destroyed the geometrical topological structure of the original triangulation fault surface. To ensure the consistency of the fault surface data, the ear-clipping algorithm was used to triangulate the triangular patches that have the intersection and extension lines of the fault surface [36,37]. The ear-clipping algorithm includes two steps: building a closed polygon with a triangle and segmenting the triangle. For triangles with intersections, a closed polygon is constructed with intersections and triangle vertices sorted based on the coordinates. The polygon is divided into multiple triangles, and the topology of the fault surface is updated.
Two fault surfaces intersect where a triangle t(A 1 B 1 C 1 ) on a certain fault surface and a triangle on another fault surface produce intersection points. The sorted intersection points constitute a set, J = {J 1 , J 2 , · · · , J n }, and a closed polygon Pc is constructed in three cases: 1.
Both  After the closed polygon Pc is constructed, traverse and delete the points on the closed polygon to generate new triangles until all new triangles are generated by segmentation.
Taking the triangle in Figure 10a as an example; there are four intersection points J 1 , J 2 , J 3 , J 4 in a triangle t(A 1 B 1 C 1 ), and the intersection points at both ends are within the triangle, forming a closed polygon P = {A 1 , J 1 , J 2 , J 3 , J 4 , J 4 , J 3 , J 2 , J 1 , A 1 , B 1 , C 1 , A 1 }. After deleting the points J 1 , J 2 , J 3 , J 4 , J 4 , J 3 , J 2 , J 1 one by one, divide the polygon P into new triangles Figure 10b. closed polygon to generate new triangles until all new triangles are generated by segmentation.

as shown in
Taking the triangle in Figure10a as an example; there are four intersection points   After the closed polygon Pc is constructed, traverse and delete the points on the closed polygon to generate new triangles until all new triangles are generated by segmentation.
Taking the triangle in Figure10a as an example; there are four intersection points

Judging the Truncated Part
When the relationships between the intersecting faults are defined, the minor fault is appropriately truncated. Taking the intersection line and its extension as the boundary, the minor fault can be divided into two parts: the upper and lower parts, and the truncated part is separated from the fault model. Usually, the smaller spread part is truncated. There are some labels: The upper part is truncated as 'up'; the lower part is truncated as 'down'; and not truncated as 'none.' Obviously, the lower part of F 2 in Figure 12a will be truncated and marked as 'down' in the relation XML. However, if the truncated part of a certain fault is known from prior geological knowledge, then the known part is truncated directly without judgment. For example, it is known that the upper part of F 1 in Figure 12a is truncated by F 3 , and the upper part after the truncation is separated, marked as 'up', with the black track on F 3 as truncation line. are some labels: The upper part is truncated as 'up'; the lower part is truncated as 'down'; and not truncated as 'none.' Obviously, the lower part of 2 F in Figure 12a will be truncated and marked as 'down' in the relation XML. However, if the truncated part of a certain fault is known from prior geological knowledge, then the known part is truncated directly without judgment. For example, it is known that the upper part of 1 F in Figure  12a is truncated by 3 F , and the upper part after the truncation is separated, marked as 'up', with the black track on 3 F as truncation line.

Complex fault Model Recording Methods
The 3D complex fault model needs to record the fault surface information and the fault-fault contact relationship. The fault surface is regarded as the entity model, and the fault surface contact relationship is the relationship model of the entities. The relationship model records the intersection information of a fault. XML is an effective tool for building entity and relationship models [35], which can accurately describe entity-related relationships. This paper proposes a method of recording complex fault models based on the XML, which can accurately describe the contact relationship between fault surfaces and effectively record complex fault information. The complex fault model recording method records information on fault surfaces, fault-fault contact relationships, and complex fault models in sequence. , and each entity represents a fault surface. Table 1 shows the record template of the fault surface entity, which contains three types of elements: the name of the fault represented by the fault ID, the geometric elements of the fault composed of vertices and triangles, and the time of fault formation.

Complex Fault Model Recording Methods
The 3D complex fault model needs to record the fault surface information and the faultfault contact relationship. The fault surface is regarded as the entity model, and the fault surface contact relationship is the relationship model of the entities. The relationship model records the intersection information of a fault. XML is an effective tool for building entity and relationship models [35], which can accurately describe entity-related relationships. This paper proposes a method of recording complex fault models based on the XML, which can accurately describe the contact relationship between fault surfaces and effectively record complex fault information. The complex fault model recording method records information on fault surfaces, fault-fault contact relationships, and complex fault models in sequence.

Record the Fault Surfaces
Construct a fault entity model to record the fault information. The set of fault surface entities E contains n entities, E = {e 1 , e 2 , . . . , e n }, and each entity represents a fault surface. Table 1 shows the record template of the fault surface entity, which contains three types of elements: the name of the fault represented by the fault ID, the geometric elements of the fault composed of vertices and triangles, and the time of fault formation. The fault entity model record template is represented by the XML structure diagram shown in Figure 13. The fault entity model is expressed as an example in xml as follows: <FaultsEntity> <FaultName>F1</FaultName> <FaultGeom> <Vertex>1721</Vertex> <Triangle>3384</Triangle> </FaultGeom> <FaultTime>Time</FaultTime> </FaultsEntity>

Record the Fault-Fault Contact Relationship
The fault-fault contact relationship was recorded using a relationship model. The set R of fault-fault contact relationships contains m relationships, R = {r 1 , r 2 , . . . , r m }, each of r m = e i , e j , which represents the intersection of two fault entities e i and e j . If there is no intersection relationship between the fault entities, no record is required. Table 2 is a record template of the fault relationship model, including four types of elements: fault relationship name, two intersecting fault IDs, major fault ID, and truncation type of minor faults. The fault relationship model record template is represented by an XML structure diagram, as shown in Figure 14.

Record the Complex Fault Model
The complex fault model is composed of a fault entities model and a contact relationship model between the fault entities. Graph theory is a mathematical model used to describe the relationship between things [38], in which undirected graphs are suitable for describing the structure of complex fault models. Using graphs to represent a complex fault model can ignore the geometric details of the fault model and only pay attention to the topological relationship between faults, making model adjustment more intuitive [33,39].
In the mathematical model, the graph is represented as: This is a two-tuple, in which V = {v 1 , v 2 , . . . , v m , v n } represents the point elements in the graph, and E = {v 1 v 2 , v 2 v 3 , . . . , v m v n } is a set of edges composed of points that have correlations in V [33].
The fault entities in the complex fault model are equivalent to the point elements in the graph, and the contact relationship between the two faults is equivalent to the edges in the graph (Figure 15). Some of the concepts in graph theory can also correspond to the problems studied in complex fault models. The degree of a point element D indicates the number of edges adjacent to the point, that is, the number of faults that intersect with the fault entity. The disconnected state in the graph theory indicates that there is no intersection relationship between the faults, and the circle indicates the self-truncating phenomenon in the fault, and the isolated point indicates that the fault has no intersection relationship in the model. In Figure 15, the degree of the point F 1 is 3, and the fault entity F 1 intersects with three faults, F 3 , F 4 , and F 5 , and the points F 1 , F 2 , and F 3 form a circle in the graph, indicating that, F 1 , F 2 , and F 3 are the self-truncating faults that intersect each other in pairs. The complex fault model associates the fault entities to the relationship model and is represented by the XML structure in Figure 16. When generating a complex fault model, the labels of the major and minor faults in the fault-fault contact relationship may be changed. The change rules are as follows: 1.
Remove isolated fault entity points and fault entity points with only two fault entities connected.

2.
Process the remaining connected fault entities sequentially. The degree of the fault entity F i is D i in the fault-fault contact relationship; r m = e i = F i , e j = F j if D i > D j , then F i is the major fault and F j is the minor fault; otherwise, it remains unchanged.

3.
If a circle is present, a self-truncating fault exists. The judgment of major and minor faults of self-truncating faults is too complicated, so they are not truncated, and its truncation label is marked as 'none.'

Examples of 3D Complex Fault Model Construction
The method proposed in this paper can solve the automatic construction of complex faults including X-type faults, multilevel y-type or λ-type faults, and self-interrupted faults. Figure 17 is a complex fault model composed of seven faults with multilevel and selftruncating characteristics. Each fault is regarded as a fault entity e i = F i , and the contact relationship between the fault entities is regarded as a relationship entity r m . The modeling process is as follows: 1.
Generate fault-fault contact relationship entities.

2.
Generate a graph model of a complex fault structure.
The point represents the fault entity, the edge represents the fault-fault contact relationship entity, and the graph model is constructed based on the fault-fault contact relationship entity shown in Figure 18.

3.
Update the label of the major fault and minor fault.
The major fault in the fault relationship entity r m was recorded as r m Ma. In the model in Figure 17, r 1 Ma = e 1 , r 2 Ma = e 1 , r 3 Ma = e 1 , r 4 Ma = e 1 , r 5 Ma = e 2 , r 6 Ma = e 4 , and r 7 Ma = e 6 . Update the major fault label according to Figure 18. Calculate the degree D of each node in Figure 18 and determine the existence of the circle C = {F 1 F 6 , F 6 F 7 , F 7 F 1 } in the graph shown in Figure 18.
After the update, the major fault in r 5 changed from e 2 to e 3 (Table 3), and the truncation label in r 3 , r 4 , r 7 changed to 'none.' Table 3. Results of updating the major fault labels.

Fault Entity
New Major Fault

Discussion and Conclusions
Based on the labeled fault point data, this study proposes and realizes the key process of automatic construction of complex fault network model, using 3D TIN fault surface model. This method has achieved good application effects in one area of China (Figure 21), replacing at least 1-2 days of manual interaction with a few minutes of computer automation. This method has the following advantages:

1.
It can automatically establish complex fault type such as x-shaped faults, multilevel y-shaped or λ-shaped faults, and self-truncating faults.

2.
It can automatically determine the fault-fault contact relationship of complex faults and independently truncate the minor faults, which greatly reduces the degree of manual participation in the model building process. 3.
XML was used to record complex fault models' information. It can be converted with a graph structure to facilitate the analysis of the contact relationship between faults in a visual environment. Prior geological knowledge can be expressed in XML and directly participate in model construction, simplifying modeling steps and improving modeling efficiency. Going forward, the standard for automatic judgment of fault type and truncation part for particularly complex fault models will still need to be further explored.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The XML record result examples of the complex fault model in Figure 15