The core idea of the regular voxel splitting model is based on the regular voxels. In the modeling process, with the filling of the voxel in the geological body, the positional relationship between the voxel and the stratum interface is judged in real time, and the voxel is split differently. That is, when the voxel is cut by the stratum interface, it will be divided into several irregular entities; and if the voxel lies completely within the interface of the two layers, its regular shape will be retained. The experimental results show that this method can compensate for the deficiency of the regular voxel in geological modeling, so that it can express the geological surface smoothly and describe the boundary position accurately.

The modeling method using the regular voxel split model is shown in figure 2:

*A. Voxel classification*

As shown in figure 3, when the voxel is filled from the top down inside the bounding box, the voxel can be divided into the following three categories according to the positional relationship between the voxel and the formation interface.

1) Complete removal (type I): Because the bounding box is larger than the actual boundary, the type I voxel generated at the beginning of filling or during the filling to the bottom does not belong to any geological body.

2) Need to split (type Ⅱ): This kind of voxel is passed through the stratigraphic interface; therefore, it needs to be split, and the entities formed after splitting belong to different geological bodies.

3) Complete retention (type Ⅲ): The voxels exist completely between two stratigraphic interfaces.

When using the regular voxel splitting model for modeling, accurately classifying all voxels and assigning geological attributes to each voxel is the key to ensuring the reliability of the model. To classify voxels, this paper uses regular grid control to interpolate at each grid point of the stratigraphic interface, and the voxels and stratigraphic points will have a positional relationship, as shown in figure 4. When the voxel is traversed by the stratum, there will be vector points on the corresponding edge. We only need to judge whether there is a stratum point on the edge, and the three types of voxels can be judged. In addition, the basic data of modeling all have geographic coordinates, and all the models are constructed are in the same coordinate system. Using the coordinates as a reference can easily judge the spatial positional relationship of voxels.

As shown in figure 4, the four edges of the voxel are respectively denoted as *e*1-*e*4. First, the two endpoint coordinates *P*1u (*x*, *y*, *z*u) and *P*1d (*x*, *y*, *z*d) of edge *e*1 are calculated from the coordinates of the center point of the voxel. Then, the vector points of the stratum interface at the corresponding lattice points are retrieved by abscissa *x* and ordinate *y*, which are marked as *P*1~*P*n (n=6 in this article), and the coordinates are recorded as (*x*, *y*, *z*u). Afterwards, by comparing the sizes of *z*u, *z*d and *z*n, the positional relationship between edge *e*1 and the formation interface point is judged, and the positions of the other three edges can be judged using the same method. If there is a vector point on the edge of the voxel, the voxel needs to be split. For example, there are vector points on edges *e*1 and *e*4 of voxel ① in the figure. If there are no vector points on the four edges of the voxel, the voxel may be removed or retained. Voxel ② in the figure should be retained because it is between two strata.

*B. Voxel split realization*

Even if they are all split voxels, due to the uncertainty of the stratum interface, the stratum will have different directions and angles when passing through the voxels. In this paper, according to the positional relationship between the upper and lower adjacent voxels and the vector points of the stratigraphic interface, the split voxels are divided into five types: "4 type", "3-1 type", "1-3 type", "2-2 adjacent type" and "2-2 crossing type". Figure 5 takes the "2-2 adjacent type" as an example to illustrate the voxel splitting process.

The flow of the entire splitting process is shown in figure 6. In the end, the voxel is split into two irregular entities. Entity ① and entity ② belong to the same voxel, but they are separated by a certain stratum; therefore, that they belong to different geological bodies. In addition, *F*1 and *F*2 shown in the figure fall on *e*3 and *e*4, respectively (reflected in figure 4). However, because of the uncertainty, it is possible to fall on other edges in other voxels, so locating the region (index) of the edge of *F*1 and using it as the benchmark to solve the model can ensure the versatility of the algorithm.

In the figure, *V*1~*V*8 represent the vertexes of the upper regular voxel, which are obtained by calculating the coordinates of the center of the voxel. *F*1~*F*4 are the vector points belonging to the same stratum interface, and the coordinates have been calculated during grid interpolation. Points *O*, *P* and *Q* are the unknown points that need to be solved, which are the points of the intersection between *F*1*F*3, *F*2*F*3, *F*1*F*4 and the lower surface of the voxel. In this paper, these points are solved by the vector method in the implementation of the program.

*Data structure design of the voxel splitting model*

To make the data organization, management and sharing of the split model more convenient, the model was divided into four basic elements in this paper: points, triangles, regular voxels and irregular voxels. Then, the data structure was designed to realize data storage and analysis.

(1) Point

In the process of constructing the geological model, two types of point data are produced. One type of point data is the vector points of the stratigraphic interface obtained by the interpolation of boreholes and sections. As shown in figure 7 (a), when designing the data structure, it is necessary to record the sequence number of the stratigraphic interface by using the layerId field in addition to recording the coordinate information (*x*, *y*, *z*). The second type of point data is the vertices of the irregular voxels generated during the split. As shown in figure 7 (b), in order to facilitate model reconstruction and visualization, the orderNum field records the number of each vertex in the voxel and binds the voxel through the vexelId field.

(2) Triangle

To achieve different visualization requirements, the organization of data should meet the requirements of the geological entity model and geological surface model at the same time. As shown in figure 8 (a), when the stratum surface passes through two upper and lower adjacent voxels, the four vector points of the same stratum controlled by the regular grid will form two triangular faces on the cutting plane of the voxel. Each triangular plane is composed of three formation vector points called LayerPoints, and all triangular planes of the same formation can be combined into a formation surface model. Figure 8 (b) is a triangular data structure. In addition to recording three vertices, the layerId field is used to identify the interface of the strata to which it belongs.

(3) Regular voxel

The regular voxel is the most common data model in the modeling process. Its characteristics are that the length and width of the horizontal direction are controlled by the grid, and the height is strictly unified. When visualizing the model, the position of each vertex and surface can be calculated from the central coordinates. As shown in figure 9, in order to save storage space, only the central coordinates (*x*, *y*, *z*) of the regular voxel are recorded when storing the regular voxel, and the geological body to which the voxel belongs is recorded using the geoBody field.

(4) irregular voxel

An irregular voxel is the most complex data model, and the entity shapes formed by the 5 splitting situations are also different. Since the vertex position and each surface of the entity cannot be solved directly, it is necessary to record each vertex and surface at the same time when designing the data structure of irregular voxels. Previously, when designing the data structure of points, the vertices of irregular voxels were stored in VexelVector, and the identification ID of the voxel was bound; therefore, when reading an irregular voxel, all the vertices of the voxel can be retrieved through their IDs. Figure 10 shows the data structure of irregular voxels. Json is a general data exchange format that can be used to store every surface of irregular voxels, and each vertex on the surface is recorded counterclockwise.

*D. Model realization*

The model implementation process is the parsing process of the data structure. How to calculate the shapes, positions and attributes of voxels through the reading of the database and the analysis of the data structure is the core task of the model. As shown in figure 11, in the experiment of this article, voxels are divided into two types: regular voxels (RegularVexel) and irregular voxels (SplitVexel). It is easy to analyze a regular voxel, and the positions of the 8 vertices and 6 surfaces can be calculated directly from its central coordinates. When analyzing irregular voxels, we should first retrieve the coordinates and numbers of each vertex according to the identification ID and obtain the vertices contained in each surface by parsing the faceJson field. Then, it is necessary to complete the rendering of each surface of the irregular voxel according to the order and coordinates of the vertices. Finally, the texture (color) is set for each voxel through the geological attribute information in the geoBody field.