The central idea of the regular voxel splitting model is based on 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 divided differently. That is, when the voxel is cut by the stratum interface, it will be split into several irregular entities; and if the voxel lies completely inside 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 regular voxel in geological modeling, so that it can express the geological surface smoothly and accurately describe the position of the boundary.
The modeling method using the regular voxel splitting model is shown in Fig. 2.
Classification of voxels
As shown in Figure 3, when the voxel is filled from top to bottom inside the bounding box, the voxel can be divided into the following three categories based on the positional relationship between the voxel and the interface of training.

(1)
Complete removal (type I): Because the bounding box is larger than the actual boundary, the type I voxel generated at the start of filling or when filling down does not belong to any geologic body.

(2)
Need to split (type II): This type of voxel has passed through the stratigraphic interface; therefore, it must be divided and the entities formed after the division belong to different geological bodies.

(3)
Complete retention (Type III): Voxels exist completely between two stratigraphic interfaces.
When using the regular voxel splitting model for modeling, accurately classifying all voxels and assigning geologic attributes to each voxel is key to ensuring model reliability. To classify the voxels, this article uses a 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 Fig. 4. When the voxel is crossed 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 all three types of voxels can be judged. In addition, the basic modeling data all have geographic coordinates, and all models built are in the same coordinate system. Using coordinates as a reference can easily judge the spatial position relationship of voxels.
As shown in Figure 4, the four edges of the voxel are respectively denoted by e_{1}–e_{4}. First, the two coordinates of the endpoints P_{1u} (X, there, z_{you}) and P_{1d} (X, there, z_{D}) edge e_{1} are calculated from the coordinates of the center point of the voxel. Then vector points from stratum interface to corresponding lattice points are retrieved by abscissa X and orderly therewhich are marked as P_{1} ~P_{not}(n=6 in this article), and the coordinates are recorded as ( X, there, z_{you}). Then, comparing the sizes of z_{you} , z_{D}and z_{not} the positional relationship between the edge e_{1} and the forming 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 must be split. For example, there are vector points on the edges e_{1} and e_{4} of voxels ① in the figure. If there are no vector points on all four edges of the voxel, the voxel can be deleted or kept. Voxel ② in the figure must be kept because it is between two strata.
Realization of the split voxel
Even though 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 article, 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 “2–2 adjacent type” as an example to illustrate the process of dividing voxels.
The flow of the whole dividing process is shown in Fig. 6. At 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 Fig. 4). However, due to uncertainty, it is possible to stumble upon other edges in other voxels, thus locating the (hint) region of the edge of F_{1} and using it as a reference to solve the model can guarantee the versatility of the algorithm.
In the figure, V_{1} ~ V_{8}represent the vertices 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 were calculated during grid interpolation. Points O, P and Q are the unknown points to solve, which are the intersection points betweenF_{1}F_{3}, F_{2}F_{3}, F_{1}F_{4}and the bottom surface of the voxel. In this article, these points are resolved by the vector method in the program implementation.
Design of the voxel splitting model data structure
To make it easier to organize, manage, and share split model data, the model has been divided into four basic elements in this article: points, triangles, regular voxels, and irregular voxels. Next, the data structure was designed to realize data storage and analysis.
Indicate
In the process of building the geological model, two types of point data are produced. One type of point data are vector points of the stratigraphic interface obtained by the interpolation of boreholes and sections. As shown in Figure 7a, when designing the data structure, it is necessary to record the sequence number of the stratigraphic interface using the layerId field in addition to recording the coordinate information (X, there, z). The second type of point data are the vertices of the irregular voxels generated during the split. As shown in Figure 7b, in order to facilitate the reconstruction and visualization of the model, the orderNum field records the number of each vertex in the voxel and links the voxel via the vexelId field.
Triangle
To achieve different visualization requirements, the data organization must meet the requirements of the geological feature model and the geological surface model at the same time. As shown in Figure 8a, when the stratum surface crosses two adjacent upper and lower voxels, the four vector points of the same stratum controlled by the regular grid will form two triangular faces on the voxel cutting plane. Each triangular plane is made up of three training vector points called LayerPoints, and all triangular planes of the same formation can be combined into a training surface model. Figure 8b is a triangular data structure. In addition to registering three vertices, the layerId field is used to identify the interface of the strata to which it belongs.
Regular voxel
The regular voxel is the most common data pattern 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 viewing the model, the position of each vertex and surface can be calculated from the center coordinates. As shown in Fig. 9, in order to save storage space, only the center coordinates (X,there,z) of the regular voxel are saved when storing the regular voxel, and the geologic body to which the voxel belongs is saved using the geoBody field.
irregular voxel
An irregular voxel is the most complex data pattern, and the feature shapes formed by the 5 splitting situations are also different. Since the position of the vertex and each surface of the feature cannot be solved directly, it is necessary to register each vertex and each surface at the same time when designing the irregular voxel data structure. Previously, when designing the point data structure, the vertices of irregular voxels were stored in VexelVector and the voxel identification ID was bound. therefore, when reading an irregular voxel, all the vertices of the voxel can be retrieved thanks to their identifiers. Figure 10 shows the data structure of irregular voxels. Json is a general data interchange format that can be used to store every irregular voxel surface, and every vertex of the surface is recorded counterclockwise.
Realization of the model
The model implementation process is the data structure analysis process. How to calculate the shapes, positions and attributes of voxels through reading the database and analyzing the data structure is the main 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 parsing irregular voxels, we first need to retrieve the coordinates and numbers of each vertex based on 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 the coordinates of the vertices. Finally, the texture (color) is defined for each voxel via the geologic attribute information in the geoBody field.