A Triangular Prism Spatial Interpolation Method for Mapping Geological Property Fields

Abstract

The spatial interpolation of property fields in 3D, such as the temperature, salinity, and organic content of ocean water, is an active area of research in the applied geosciences. Conventional interpolation methods have not adequately addressed anisotropy in these data. Thus, in our research we considered two interpolation methods based on a triangular prism volume element, as a triangular prism structure best represents directivity, to express the anisotropy inherent in geological property fields. A linear triangular prism interpolation is proposed for layered stratum that achieves a complete $C_0$ continuity based on the volume coordinates of the triangular prism. A triangular prism quadric interpolation (a unit function of a triangular prism spline with 15 nodes) is designed for a smooth transition between adjacent triangular prisms with approximately $C_1$ continuity, expressing the continuity of the entire model. We designed a specific model which accounts for the different spatial correlations in three dimensions. We evaluated the accuracy of our proposed linear and triangular prism quadric interpolation methods with traditional inverse distance weighting (IDW) and kriging interpolation approaches in comparative experiments. The results show that, in 3D geological modeling, the linear and quadric triangular prism interpolations more accurately represent the changes in the property values of the layered strata than the IDW and kriging interpolation methods. Furthermore, the triangular prism quadric interpolation algorithm with 15 nodes outperforms the other methods. This study of triangular prism interpolation algorithms has implications for the expression of data fields with 3D properties. Moreover, our novel approach will contribute to spatial attribute prediction and representation and is applicable to all 3D geographic information; for example, in studies of atmospheric circulation, ocean circulation, water temperature, salinity, and three-dimensional pollutant diffusion.

1. Introduction

Three-dimensional (3D) geological modeling (3DGM) is an area of interest and technological development, involving geology, geo-statistics, and computer science. Three-dimensional geometric shapes and the internal property distributions of geological bodies can be built with 3DGM techniques using the available geological modeling data and appropriate modeling methods. A geographic information system (GIS) platform is not only required to express the boundaries of geological bodies, but also the internal heterogeneous properties, such as the rock property fields (mechanical strength, porosity, etc.), the distribution of gas in coal beds and surrounding rocks, crustal stresses, temperature fields, grade changes in ore, and so on. A 3D property field is an abstract mathematical definition of the spatial distribution of geological properties. However, the existing GIS platforms can only describe the vector structures [1,2,3] of a geological model, and cannot achieve a heterogeneous expression of a 3D property field. Furthermore, due to the relatively high cost of underground sampling, sparsely distributed sampling points are often used to infer the distribution of the properties of a study area. As a result, 3D spatial interpolation in most geological models is undertaken using sparsely distributed sampling points.

Stratification is one manifestation of heterogeneity and is a critical element in geology. Materials with a layered structure always show large differences in property change rates between the horizontal and vertical directions. The change rate in the horizontal direction is relatively uniform and slow, while it is rapid in the vertical direction [4]. The mathematical expression of these heterogeneous properties benefits from the accurate calculation and analysis of geological property fields, such as organic matter measurements from bore samples. Improved methods for visualization and 3D modeling support advances in many practical applications, such as mining and oil and gas exploration.

Much research has been done on the interpolation of the property fields of geological bodies; Li has proposed a “volume function model” [5,6]. The volume function model complements traditional vector boundary and volume element models, splitting a geological body into adjacent and non-intersecting elements (such as tetrahedrons, triangular prisms or hexahedrons). A volume function can be built to mathematically express continuous change, the distribution and the inter-relationship of heterogeneous property fields. This interpolation method is faster and has a greater accuracy than IDW and kriging interpolation, as the heterogeneity of a geological body can be expressed by the volume function of the non-overlapping volume elements. Three-dimensional spatial interpolation using a fitted volume function ensures $C_0$, $C_1$ or $C_2$ continuity in the same layer rather than between different strata [5]. Interpolated property values using volume function satisfy positional continuity, first derivative continuity, or second derivative continuity, and are similar to curve and surface interpolation $C_0$, $C_1$ and $C_2$ continuity. Previous research has been done on the volume function and related fields. Zhou and Chen [7,8] used a tetrahedral region high-order volume function to fit values of a property to an entire stratum. Alfeld [9] proposed first-order smoothing ($C_1$) based on the Clough–Tocher tetrahedral interpolation method; Barnhill and Little [10] proposed a tetrahedral interpolation formula based on the Barnhill, Birkhof, and Gordon (BBG) transfinite interpolation method; Soh [11] successfully applied quadrilateral area coordinates to the thin plate element; Li [12] studied the interpolation discretization of two types of hexahedral volume coordinates based on mechanical 3D finite element calculation theory; Wang [13] used a piecewise spline polynomial that satisfied the continual condition in the finite element method. Bhattacharjee [14] proposed semantic kriging to blend the semantics of spatial features with an ordinary kriging method for prediction of an attribute. Experimentation was carried out with land surface temperature data of four major metropolitan cities in India. Anderson demonstrated that kriging produces a better estimation of a continuous surface of air temperature than IDW and spline [15]. A radial basis function in a neural network improved spatial interpolation performance more than traditional methods for soil properties in Kunming [16]. Karydas evaluated prediction maps created with interpolation for five common topsoil properties; namely, organic matter, total CaCO₃, and electric conductivity [17]. However, all these methods have a common drawback in that the directivity of their heterogeneous attributes cannot be effectively expressed.

To handle the directivity of heterogeneous property changes, an appropriate body element needs to be selected as the base component unit for a geological body. A tetrahedron is suitable for complex strata modeling and can represent the heterogeneous property changes of a local area. However, it has several drawbacks for geological modeling. A geological model contains a large amount of data, and using a tetrahedron is time-consuming. Furthermore, a tetrahedron element cannot express the directivity of a geological body [5]. A hexahedron is also not appropriate, as its regular formation is not suitable for complex geological strata [12]. In this study, we selected the triangular prism as the body element for the modeling. A triangular prism can be decomposed into the basic tetrahedron elements and can fit the features of borehole data and the trends found in geological strata. It can also effectively describe regional and directional changes in the properties of a geological body.

This paper explores different triangular prism interpolation algorithms. Unlike conventional geostatistical methods, these methods are based on volume element interpolation. A new triangular prism linear interpolation method is proposed to express the properties of a layered stratum that achieves a complete $C_0$ continuity in a single triangular prism. A triangular prism quadric interpolation algorithm is designed for a smooth transition between adjacent triangular prisms, expressing the continuity of the whole model, reaching approximate $C_1$ continuity. We deduce the triangular prism linear interpolation formula based on the existing methods for tetrahedron and hexahedron volume coordinates [5,12]. The triangular prism quadric interpolation is designed by modifying the existing B-net method for tetrahedron domains [7,18,19]. In our research, the triangular prism element is used to express the 3D geometric shapes of geological bodies and structural mechanics, building on these past studies [20,21]. The related tetrahedron and hexahedron interpolation methods have been extensively researched [7,8,12]. Triangular prism interpolation, however, is a new solution for interpolating and predicting geological property fields based on sparse sample data.

The remainder of this paper is organized as follows. Section 2 summarizes several related geostatistical interpolation methods. Section 3 introduces the experimental data and algorithms. Section 4 describes the geostatistical model for 3D interpolation. Section 5 presents the results. In Section 6, the results are discussed, and we end with some concluding remarks.

2. Related Work

Interpolation issues are a focus of much work in the geostatistical field. Royse [22] argues that three-dimension modeling is developing in two directions: digital three-dimension interpolation and cognitive interpretation. Three-dimensional data interpolation is a heterogeneous expression of 3D properties. Scholars from different fields have done similar research on the 3D modeling of heterogeneous data fields; the main methods are the distance power inverse ratio, kriging interpolation, finite element and spline interpolation.

(1) Distance Power Inverse Ratio Interpolation

The distance power inverse ratio method was proposed by meteorologists and geologists. It is simple, and has a wide applicability. The closer the space points are, the stronger the similarity, and the farther the distance, the weaker the similarity. This indicates that the unknown point is similar to the values of the nearest $n$ known points. The property value $F(x,y,z)$ of the interpolation point is defined as the weighted average value of the known point property $f_k$, which treats the distance as inverse between the interpolation point and the known points, as shown in Equation (1).

$$F(x,y,z)={\displaystyle \sum _{k=1}^n{f}_k{W}_k(x,y,z),}{W}_k(x,y,z)=\frac{\frac{1}{d_k{(x,y,z)}^u}}{{\displaystyle \sum _{k=1}^n\frac{1}{d_k{(x,y,z)}^u}}}$$

(1)

where $W_k(x,y,z)$ is the weight of known points; $d_k(x,y,z)=\sqrt{{(x-x_k)}^2+{(y-y_k)}^2+{(z-z_k)}^2}$ is the Euclidean distance from $(x,y,z)$ to $(x_k,y_k,z_k)$; $u$ is the power exponent. The accuracy is acceptable when the sample points are adequate. In terms of its application, engineering practice shows that this interpolation is an effective method for simple-structure stratiform deposits; however, this method has some problems stemming from “outliers”, “decluttering” and “anisotropy” [23].

(2) Kriging Interpolation

Kriging is also called the spatial autocovariance best interpolation method, developed by the French geographical mathematician Georges Matheron and South African mining engineer D.G. Krige, and is a mature, unbiased, optimal interpolation method widely used in geo-statistics for the calculation of ore reserves, hydrology, meteorology, topography, and soil mapping. This geostatistical method provides an optimization strategy for spatial interpolation; in the interpolation process, the variables are determined dynamically according to set optimization criteria. Matheron and Krige focused on the determination of the weight coefficients, so that the interpolation function is in the best state: the best linear unbiased estimate for the variable values at the given point. Kriging is based on the known points to estimate the unknown points with a minimum prediction error. In addition, the corresponding minimum estimation variance can be obtained. The Kriging method however, is complex and requires a great deal of computation time. The kriging interpolation effect depends on the selection of variation functions in the theoretical model. However, a variogram appropriate for a geological model is not known in advance, thus the human factor in the kriging method is strong, as a user determines if the variogram function is reasonable or not. Some scholars [24,25,26] have extended the kriging method, but the kriging method is still a low-pass filter and cannot reconstruct local, high-frequency, or weak signals reflected by an original signal.

(3) Finite Element Interpolation

Geostatistical methods interpolate subsurface property fields by modeling voxels [27], expressing a continuous geological field using discrete adjacent volume elements such as tetrahedrons, triangular prisms, or hexahedrons. A volume element is only a property value, and is very small; thus, if expressed reasonably and effectively, it can represent the heterogeneous properties in the real body. However, with the gradual division of the volume element, the storage space becomes huge, and the amount of data will expand three-fold, while the computation speed is greatly reduced. However, the approximation function of a volume element is expressed by the node of each element, and its interpolation function in an unknown field effectively and efficiently expresses internal features of the data. Li, Cen [12,28] and others studied the interpolation methods for two kinds of hexahedral volume coordinates from the point of view of mechanical three-dimensional finite element calculation. The tetrahedron unstructured model is suitable for complex structural models and finite element simulations based on irregular discrete points [29]. In geology, a triangular prism can be decomposed into the basic tetrahedron elements to fit the features of borehole data and the trends of the strata. It can also effectively describe the regional and directional changes of the properties of a geological body.

(4) Spline Interpolation

In the numerical analysis, Spline interpolation is a form of interpolation where the interpolant is a special type of a piecewise polynomial called a spline. Spline interpolation is often preferred over polynomial interpolation because the interpolation error can be reduced even when using low degree polynomials for the spline [13]. Spline interpolation avoids the problem of Runge’s phenomenon, as oscillation can occur between points when interpolating using high degree polynomials. The main advantages of spline interpolation are its stability and calculation simplicity. Sets of linear equations used to solve to construct splines are very well-conditioned; therefore, the polynomial coefficients are calculated precisely.

In the existing body of work, however, little attention has been paid to methods that describe the regional and directional changes of the properties of a geological body. To fill this gap, we propose two interpolation methods based on a triangular prism volume element, as a triangular prism structure best represents directivity, to express geological property fields.

3. Materials and Methods

3.1. Materials

There are many types of original data that can be used for 3DGM, including borehole data, contour line data, fault line data, and profile data. This paper employs borehole data from a reliable source. The original geology three-dimensional assistant engineer (G3DA) text file contains the original document for modeling.

Firstly, the data including the 3D coordinates and the data of four stratum layers were obtained from 365 borehole points from a mining area in Qinghai Province, China. After establishing the connected relationships of the borehole points, a triangular surface model of the point data was constructed by triangulation, confirming the digital elevation map (DEM) of the four stratum layers. We plotted and drew the triangular prisms using the geometric relationships of points, lines and faces in DEM. The triangular prism model is shown in Figure 1.

To validate the following interpolation algorithms, we simulated the property fields according to the change characteristics; propagating the property value, however, along the vertical direction is faster than in the horizontal direction. Hence, different damping coefficients between the horizontal and vertical directions were used in order to propagate the property outwards. The property type is organic matter content. Assuming that this property propagates along the horizontal and then to the vertical, the propagation progress can be expressed as:

$$v_1=P_a/(1.0+d_1*r_1)$$

(2)

where $P_a$ denotes the maximum property value, $d_1$ denotes the largest distance between the sample point of the triangular prism and the maximum property point along the horizontal direction, and the corresponding minimum property value is $v_1$. $r_1$ denotes the horizontal change rate, given by:

$$r_1=(P_a-P_i)/(P_i*d_1)$$

(3)

where $P_i$ equates to $v_1$.

$v_2$ is the property value of the corresponding $v_1$ with distance $d_2$ in the vertical direction. Its vertical change rate is $r_2$, and expressed as follows:

$$v_2=v_1/(1.0+d_2*r_2)$$

(4)

where $r_2$ denotes the vertical change rate, given by:

$$\begin{array}{l}{r}_{21}=(P_a/P_{ai}-1.0)/d_{21}\\ r_{22}=(P_i/P_{ii}-1.0)/d_{22}\\ r_2=(r_{21}+r_{22})/2\end{array}$$

(5)

where $P_a$, $P_i$ are the maximum and minimum property values in the horizontal direction, respectively. $P_{ai}$ is given the property value of the corresponding $P_a$ with distance $d_{21}$ in the vertical direction, and $P_{ii}$ is the same as $P_{ai}$. Figure 1a is the original triangular prism model. The simulated triangular prism property field is shown in Figure 1b.

3.2. Triangular Prism Linear Interpolation

Tetrahedron volume coordinates have been studied in the fields of mathematics and mechanical engineering for years. Li [5,6] undertook in-depth research on this topic from a geographic information perspective. Much of the existing research had addressed hexahedron volume coordinates [12]. In this paper, we deduce the triangular prism volume coordinate formulas by studying the ideas and methods of the related work discussed in the introduction. In this section, we describe the linear interpolation algorithm based on triangular prism volume coordinates.

3.2.1. Triangular Prism Volume Coordinates

(1) Generalized Triangular Prism (GTP) Shapes

The generalized triangular prism (GTP) is a 3D element, as illustrated in Figure 2a. The GTP element is degenerated when there is a fault in the strata or a narrow distance between two-layer boundaries [20], as shown in Figure 2b,c.

(2) The Node and Surface Serial Numbers of the Triangular Prism Element

The triangular prism element node and surface serial numbers when the serial numbers of the points at the surface follow the “right-hand screw rule” [12] are shown in Figure 3.

(3) The Basic Definition of Triangular Prism Volume Coordinates

The convex pentahedron consists of quadrilateral lateral surfaces and two triangular bottom surfaces, as shown in Figure 3. Point P is given randomly and has volume coordinates that are similar to the tetrahedral and hexahedron volume coordinates [5,12]:

$$L_i=\frac{V_i}{V}\ge 0(i=1,2,\dots ,5)$$

(6)

where $V_i$ is the volume of the rectangular pyramid or tetrahedron comprised by point P and the $i$-th surface, and $V$ is the volume of the triangular prism. Therefore, these volume coordinates satisfy $L_1+L_2+L_3+L_4+L_5=1$, expressed as $(u,v,w,r,s)$. The position of any point on the triangular prism can be determined by its volume coordinates. The volume coordinates are independent of the shape of the triangular prism and interpolated point position in a global coordinate system.

A random triangular prism can be split into three adjacent and disjoint tetrahedrons without degeneration; thus, the volume of the triangular prism can be expressed as the summation of the three tetrahedron volumes:

$$V_{Tri-prism}=V_{\mathsf{\Delta }4651}+V_{\mathsf{\Delta }1652}+V_{\mathsf{\Delta }1623}$$

(7)

where $V_{\mathsf{\Delta }4651}=\frac{1}{6}(\overrightarrow{45}*\overrightarrow{46})\cdot \overrightarrow{41}$.

3.2.2. Interpolation Principle

In the following discussion, the property value of arbitrary interpolation point P is deduced by the coordinates and property values of six vertex points of the triangular prism. We call the method “triangular prism linear interpolation”. The basic idea of this algorithm is illustrated in Figure 4. The line perpendicular to the vertical lines of edges 1–4, 2–5, and 3–6 passes through point P. Points 7–9 are perpendicular feet. Points 10 and 11 are formed by the vertical edge line that is perpendicular to the horizontal plane through point P intersecting surface 4 and 5. From Figure 4, it can be observed that the distance from the interpolation point to the surface affects the triangular prism volume coordinates. The greater the distance from interpolation point P to surface 1, the greater the proportion the rectangular pyramid occupies in the total volume of the triangular prism. This also means that the larger the value of $u$, the closer the distance to point 9 at edges 3–6, which means that the property value of point P will more likely be affected by point 9. Similar procedures can be used to explain that the values of points 8, 7, 11, and 10 will have more impact on point P.

Because the values of points 7–11 are not given, we need to use interpolation techniques to estimate them. The basic idea of this algorithm is as follows. The property values of points 7–9 can be obtained according to linear interpolation, as shown in Equation (8).

$$\{\begin{array}{c}{m}_7=m_4+\frac{(z_7-z_4)}{(z_1-z_4)}*(m_1-m_4)\\ \begin{array}{l}{m}_8=m_5+\frac{(z_8-z_5)}{(z_2-z_5)}*(m_2-m_5)\\ m_9=m_6+\frac{(z_9-z_6)}{(z_3-z_6)}*(m_3-m_6)\end{array}\end{array}$$

(8)

Meanwhile, the property values of points 10 and 11 can be determined by linear interpolation based on the triangle area coordinates [5]. $(u_1,v_1,w_1)$ are the area coordinates of point 10 at surface 4. The property value of point 10 can then be represented as follows:

$$m_{10}=m_1*u_1+m_2*v_1+m_3*w_1$$

(9)

A similar procedure can be used to estimate the property value of point 11.

Point P can be in different positions in the triangular prism. Therefore, we must consider different situations:

(1) The property value function of point P is as shown below when it is located in the inside of the triangular prism:

$$m_p=u\ast m_9+v*m_8+w*m_7+r*m_{11}+s*m_{10}$$

(10)

(2) The property value function of point P is as shown below when it overlaps point 1, where its coordinates are $(0,v,0,0,s)$:

$$m_p=m_1*(v+s)=m_1*1=m_1$$

(11)

(3) When point P is on edge 1–4, as shown in Figure 5a, its volume coordinates are $(0,0,w,r,s)$. Point P overlaps point 7. Therefore, the property value of point P is:

$$m_p=m_1*s+m_4*r$$

(12)

The result shows that a property value is only related to the corners of the edges when the interpolation point is located on the edge of the triangular prism.

(4) When point P is located on surface 4, as shown on Figure 5b, its volume coordinates are $(u,v,w,0,s)$. The property value of point P strongly depends on the property values of points 1–3, and its area coordinates are $(u_1,v_1,w_1)$ at surface 4. Therefore, its property value is expressed by the triangular area coordinates function, as follows:

$$m_p=m_1*u_1+m_2*v_1+m_3*w_1$$

(13)

(5) When point P is located on surface 1, its volume coordinates are $(0,v,w,r,s)$, as shown in Figure 5c. Similarly, the property value of point P is correlated to the property value of points 7, 8, 10, and 11. Therefore, the property value of interpolation point P is:

$$m_p=v*m_8+w*m_7+r*m_{11}+s*m_{10}$$

(14)

(6) When the triangular prism degenerates to a pyramid, as shown in Figure 5d, the pyramid consists of point P and surface 2, and the calculation method of the property value of point P is the same as (1).

(7) When the triangular prism degenerates to a tetrahedron, as shown in Figure 5e, surface 3 overlaps surface 4. It forms four tetrahedrons with surfaces 1, 2, 4, and 5, and its volume coordinates are $(u,v,0,r,s)$. Therefore, the calculation method for the property value of point P is the same as the linear interpolation of tetrahedron volume coordinates [5].

3.3. Triangular Prism Quadric Interpolation

The triangular prism linear interpolation can achieve $C_0$ continuity, but there are still some problems that need to be resolved. For example, the triangular prism linear interpolation algorithm cannot express the smooth transition of interior properties; therefore, it cannot easily convey the smooth transitions in a geological body. To overcome this weakness, we must find a way to represent the directivity and smooth transition of geological properties. To realize this goal, a triangular prism quadric interpolation approach with $C_1$ continuity is selected for further research.

Since a triangular prism with six points is not simplex, it is difficult to obtain the fitted smooth quadric interpolation by the volume coordinates of the triangular prism listed. Li and Cen [12] studied the volume coordinates of the hexahedron, but they did not upgrade the volume coordinate interpolation method to a higher order.

In mathematics, splines are widely used in finite element methods, as they satisfy certain continuity criteria of piecewise polynomials [18,21]. Li [21] proposed a quadrangle spline element based on triangle area coordinates and the B-net method. These 2D splines are highly accurate and insensitive to mesh distortion [30]. To upgrade a triangular prism element from 2D elements, internal smoothing can be achieved by constructing a triangular prism spline element by the B-net method [21].

3.3.1. Interpolation Principle

Bernstein–Bézier tetrahedron polynomials, the B-net method, and the splice method between tetrahedral domains are the theoretical bases of the quadric continuity for the internal triangular prism element [19,31,32,33]. The Bernstein–Bézier tetrahedron is an extension of the Bernstein–Bézier triangular surface in 3D space. It is popular in mathematics and engineering as it solves mechanical problems efficiently [34]. The B-net method is used to create a multivariate spline based on triangles or tetrahedrons. It originated from work on Bernstein polynomials for interpolating barycentric coordinates for triangular areas or the tetrahedral volumes [21].

This section focuses on how to apply these methods to the quadric continuity of the interior triangular prism element and build the 15 spline bases in a triangular prism quadric interpolation. Given that the property values of the six points are $g_1,g_2,g_3,g_4,g_5,g_6$ in the triangular prism shown in Figure 6b, then the triangular prism can be split into three sub-tetrahedrons, which are ${\mathsf{\Delta }}_1=p_6{p}_1{p}_4{p}_5$, ${\mathsf{\Delta }}_2=p_6{p}_1{p}_5{p}_2$, and ${\mathsf{\Delta }}_3=p_6{p}_1{p}_2{p}_3$.

$P_7,P_8,\cdots P_{18}$, $g_7,g_8,\cdots g_{18}$ are the coordinates and property values at the midpoints of the edges, and are calculated by averaging the coordinates and property endpoints for each edge. Each tetrahedron has 10 domain points according to the quadratic polynomial, as estimated using the tetrahedral volume coordinates and the B-net method, shown in Figure 6a. Hence, there are a total of 18 domain points placed in the prism, as shown in Figure 6b. The related B-net coefficients are $b_1,b_2,\dots ,b_{18}$.

In order to make spline functions of the $C_1$ continuity at points $p_{16},p_{17},p_{18}$, 15 interpolated spline bases can be obtained according to the smoothing cofactor method [18]. This means that the B-net coefficients of points $p_{16},p_{17},p_{18}$ are confirmed by these 15 points $p_1,p_2\dots ,p_{15}$. This splice method between the tetrahedral domains and the triangular prism interior meets the $C_1$ continuity.

The calculation processes of the B-net coefficients of the 15 interpolated spline bases in the triangular prism are listed below, and the constants are defined as:

$$\begin{array}{l}{\delta }_1=Det(4,5,2,6),{\delta }_2=-Det(1,5,2,6),{\delta }_3=Det(1,4,2,6)\\ {\delta }_4=-Det(1,4,5,6),{\delta }_5=Det(1,4,5,2)\\ {\eta }_1=Det(5,6,3,1),{\eta }_2=-Det(2,6,3,1),{\eta }_3=Det(2,5,3,1)\\ {\eta }_4=-Det(2,5,6,1),{\eta }_5=Det(2,5,6,3)\end{array}$$

(15)

The Bézier factors of $b_{16},b_{17},b_{18}$ are gained from $b_1,b_2,\dots ,b_{15}$:

$${(b_{16},b_{17},b_{18})}^T=-A_1^{-1}{A}_2{(b_1,b_2,\dots ,b_{15})}^T$$

(16)

$$A_1=\left(\begin{array}{ccc}{\delta }_3& 0& {\delta }_5\\ 0& {\eta }_3& 0\\ {\eta }_2& 0& {\eta }_3\end{array}\right),A_2=\left(\begin{array}{ccccccccccccccc}{\delta }_1& 0& 0& 0& 0& 0& {\delta }_2& 0& 0& {\delta }_4& 0& 0& 0& 0& 0\\ 0& {\eta }_1& 0& 0& 0& 0& 0& {\eta }_2& 0& {\eta }_5& {\eta }_4& 0& 0& 0& 0\\ {\eta }_5& 0& 0& 0& 0& 0& 0& 0& 0& {\eta }_1& 0& {\eta }_4& 0& 0& 0\end{array}\right)$$

(17)

By solving for unknown values of $b_1\dots b_{18}$, we can gain 15 linear independent vectors $b^{(i)}={\left\{b_j^{(i)}\right\}}_{j=1}^{18}(i=1,\dots 15)$; therefore, the B-net coefficients are gained, which form a 15 × 18 matrix:

$$\left(\begin{array}{c}{b}^{(1)}\\ b^{(2)}\\ ⋮\\ b^{(15)}\end{array}\right)=(I|A)$$

(18)

where $I$ is a 15 × 15 identity matrix, and $A=-{(A_1^{-1}{A}_2)}^T$ forms a 15 × 3 submatrix. Therefore, 15 spline bases $L_i(i=1,\dots ,15)$ restricted on the three sub-tetrahedrons can be obtained by the B-net method from the corresponding Bezier coefficients. The 15 spline functions $L_i(i=1,\dots ,15)$ satisfy ${\displaystyle \sum _{\mathrm{i}=1}^{15}{L}_i=1}$. By a simple linear transformation, we obtain another set of interpolated spline bases corresponding to the nodes $P_{\mathrm{i}}=(x_i,y_i,z_i)(i=1,\dots ,15)$ as follows:

$$\begin{array}{l}{N}_1=L_1-\frac{1}{2}(L_7+L_{10}+L_{12})N_2=L_2-\frac{1}{2}(L_8+L_{10}+L_{11})\\ N_3=L_3-\frac{1}{2}(L_9+L_{11}+L_{12})N_4=L_4-\frac{1}{2}(L_7+L_{13}+L_{15})\\ N_5=L_5-\frac{1}{2}(L_8+L_{13}+L_{14})N_6=L_6-\frac{1}{2}(L_9+L_{14}+L_{15})\\ N_7=2L_7,N_8=2L_8,\dots ,N_{15}=2L_{15}\end{array}$$

(19)

The 15 triangular prism spline bases are constructed by $N_1,N_2,\dots N_{15}$. The expression of the property value of the interpolation point is written as:

$$G(x,y,z)={\displaystyle \sum _{i=1}^{15}{g}_i{N}_i(x,y,z)}$$

(20)

3.3.2. Adjustment of Property Values

The triangular prism interpolation described above only possesses second-order completeness ($C_1$ continuity) in a single triangular prism element. However, $C_0$ continuity is obtained between neighboring triangular prisms. A geological model is structured by large numbers of triangular prism elements, and the smooth transition between adjacent triangular prism elements can express the continuity of the whole model. Therefore, we have to consider the transition between adjacent triangular prism elements in different directions. We adjust the property value of each quadric control point of the triangular prism in the horizontal and vertical and thus achieve approximate $C_1$ continuity.

(1) Horizontal Adjustments

To solve the smooth transition of the property values of adjacent triangular prisms in the horizontal direction, we select the $C_1$ continuity condition of the triangle patch based on the triangle area coordinates as the adjustment method, as the change rate of the property value is smaller along the horizontal layer direction. As shown in Figure 7a, point P of the triangular prism adjoins to triangular prisms 2, 4, 5, and 6. Therefore, the area coordinates of point P are $u_i,v_i,w_i(i=1,2\cdots ,n)$ at contiguous triangle surfaces. The property value $m_{pi}(i=1,2\dots n)$ can be obtained according to Equation (21):

$$m_{p1}=m_{20}*u_1+m_{21}*v_1+m_{22}*w_1$$

(21)

where $i=1$, $m_{20},m_{21},m_{22}$ are the property values of the corners of triangular prism 2.

It is difficult to achieve complete $C_1$ continuity in adjacent triangular prisms, as a simple triangular prism has many neighboring triangular prisms and requires many adjustments. The $C_1$ continuity condition only aims at one-time adjustment. Considering the geological application, the final property value is calculated by a weighted average of all the continuity adjustments, expressed as:

$$m_p=\frac{1}{n}{\displaystyle \sum _{i=1}^n{m}_{pi}}$$

(22)

where n is the number of adjustment times.

(2) Vertical Adjustments

The change rate of the property values along the vertical layer is quite large. Therefore, when calculating the property value of midpoint $Q$ at edge $AB$, it is necessary to consider the contribution of edge points $A,B,C,D$ to point $Q$, as shown in Figure 7b. In a geological model, triangular prisms are not always straight triangular prisms, but the coordinates $x,y$ at the nodes on the vertical layer are very close. Subsequently, the property value of the second control point in the vertical direction depends on the height changes. We can determine the cubic spline interpolation of the heights and property values, finally attaining the property value of the second control point.

The property value of an interpolation point can be gained by the adjusted property values and the corresponding basis functions of the 15 nodes.

4. Geostatistical Model of 3D Interpolation

Anisotropy is common phenomenon in geoscience applications [4]; therefore, it is necessary to design a specific model that accounts for spatial correlations in three dimensions. The variogram is a basic visualization tool to represent spatial correlations in geostatistics. In terms of a regionalized variable, a variogram reflects the spatial variation characteristics of properties, especially the structural characteristics. We can explore geological phenomena by fitting the variogram model. The variogram guarantees an effective spatial prediction in interpolation methods, most notably in kriging [35].

The variogram is estimated according to the coordinates and properties of a set of sampling points in geological space, as shown in Equation (23).

$$r\left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle \sum _{i=1}^{N\left(h\right)}{\left[m\left(x_i+h\right)-m\left(x_i\right)\right]}^2}$$

(23)

where $N\left(h\right)$ is the number of pairs of sample data from a distance of $\left|h\right|$. $m\left(x_i+h\right)$ and $m\left(x_i\right)$ are the property values located in $x_i+h$ and $x_i$$\left(i=1,2,\cdot \cdot \cdot ,N\left(h\right)\right)$. A variogram is a variation map which plots change in $r\left(h\right)$ together with the changes in $h$, as shown in Figure 8.

Figure 8. Variogram.

Figure 8. Variogram.

where $C+C_0$ is the still, $C_0$ is nugget, $C$ is the arch height, and $a$ is the range.

When interpolating 3D spatial properties, we use the eighth-sphere method to search all the near neighbors of the interpolation point with an appropriate search radius [35]. Given the heterogeneity and geometrical anisotropy of 3D geological property fields, the direction of the variable range will be an ellipsoid, as an ellipsoid generalizes the direction and size of the autocorrelation structure of the data represented by the variogram. In our experimental data, organic matter displays a uniform and slow change in the horizontal direction, while rapid change occurs in the vertical direction. We infer that the horizontal direction is the minimum correlation and the vertical direction is maximum correlation. Variograms of organic matter are shown in Figure 9. We find that the range of the vertical direction (the direction of the maximum) is shorter than the horizontal direction (the direction of the minimum). This indicates that organic matter changes faster with stronger anisotropy and heterogeneity in the vertical direction. In addition, the still is larger, indicating that the magnitude of the change in organic matter is greater and anisotropy stronger in the vertical direction than in the horizontal direction. Our experimental results are consistent with the statistical theory [36]. In kriging interpolation, we use the fitting variogram model to obtain weights for the sampling points and interpolation point properties.

5. Results and Discussion

The experimental design is similar to that used by Zhou and Chen [7,8]. The performance of the interpolation tested algorithms was evaluated by the profile results, average errors, and standard errors. We used the cross-validation method to assess the accuracy and stability of the model.

5.1. Profile Results

The profile is not only an effective method to express the changes of morphology and properties of the volume elements, but also an effective mean for visualization and inspection of 3D interpolation results. In order to test the performance of the interpolation algorithms, we simulated the property profile, as shown in Figure 10. Horizontal profiles and fence diagrams show the profiles of the non-uniform property changes. The profile results from the interpolation algorithms are shown in Figure 11 and Figure 12.

In these profiles, the quadric triangular prism interpolation results seen in Figure 11d and Figure 12d attest to the faster and more accurate performance of the proposed algorithm compared to the other three algorithms tested. These results are closest to the real values and show the most detail. The quadric triangular prism interpolation expresses, at a high level of detail, heterogeneity in the horizontal and vertical directions for the organic matter found in the samples. The second-best algorithm is the linear triangular prism interpolation algorithm, as shown in Figure 11c and Figure 12c, which reveals both overall and local property changes. Figure 11b and Figure 12b show that kriging interpolation expresses the overall situation. The IDW method is the slowest and the least accurate; the interpolation results are sparse with large differences between the interpolated results and measured changes in organic matter content, as shown in Figure 11a and Figure 12a. In contrast to the kriging and IDW interpolation results, the linear triangular prism and quadric triangular prism interpolation algorithms are less time-consuming, with more distinct local property values for organic matter.

5.2. Statistical Errors

In order to evaluate the four interpolation algorithms more precisely, we calculated the error, average error, and standard deviation error of the property values. The error distributions for the interpolation results (shown in the Figure 11 horizontal profiles) are displayed in Figure 13. There are 21,873 interpolation points, and the unit of organic matter error is (‱). Therefore, the unit of the vertical coordinates is (‱), and the unit of the horizontal coordinates is the number of points.

$p_i$ is the interpolated property value, $m_i$ is the simulated data, and n is the number of samples.

The error formula is:

$$\mathsf{\Delta }v=p_i-m_i(i=1,2,\dots ,n)$$

(24)

The average error formula is:

$$v={\displaystyle \sum _{i=1}^n(\left|p_i-m_i\right|)}/n$$

(25)

The standard deviation of the error formula is:

$$\sigma =\sqrt{{\displaystyle \sum _{i=1}^n{(\left|p_i-m_i\right|-v)}^2/n}}$$

(26)

In Figure 13, considering stability, the closer to the zero-horizontal line, the more stable the performance of the algorithm. We can see that the triangular prism quadric interpolation algorithm achieved the greatest stability, as shown in Figure 13d. Figure 13c show that the error range of the triangular prism linear interpolation is larger than the quadric interpolation; however, most of the errors are close to the zero-horizontal line. Figure 13a,b show that IDW and kriging perform are less accurate, as suggested by the error distribution.

For a more quantitative comparison of the four algorithms, the average errors and standard errors of all the algorithms are listed in

Table 1. Average error and standard error statistics.

Interpolation Algorithm	Sample Number	Average Error (‱)	Standard Error (‱)
Inverse distance weighting	21,873	12.66	8.91
Kriging	21,873	11.36	5.81
Linear triangular prism	21,873	11.68	5.29
Quadric triangular prism	21,873	8.99	3.32

.

5.3. Cross-Validation Method

Cross-validation can determine the accuracy and stability of a model. All known property points were divided into a training set and testing set with the cross-validation method for prediction, using the interpolation methods discussed in previous sections. The precision of the interpolation algorithms was compared using the correlation coefficient (R²) of the new interpolation values to the sample values. As shown in Figure 14a–d, for 37 sampling points, the quadric triangular prism interpolation (R² = 0.9703) is nearest the actual value and the points in a scatter plot are concentrated nearest the zero-error line, reflecting an unbiased valuation. The linear triangular prism interpolation (R² = 0.8443) is followed by kriging (R² = 0.8457), and the IDW method (R² = 0.8078) is last.

5.4. Results Analysis

Kriging is effective at whole geological body interpolation as it is a linear unbiased method that accounts for the range of influence of the interpolation point. It is a common method in geoscience, although the variation theory is complex. In our experiment, the kriging interpolation process requires more time, as it searches the nearest point from the all sample points. In contrast, quadric triangular prism interpolation is fast; it only needs to find the triangular prism for the interpolation points and the adjacent triangular prisms. This method focuses on building the spline bases of each triangular prism volume element and creating a smooth transition between adjacent triangular prisms. Therefore, quadric triangular prism interpolation has an acceptable level of precision and better represents the anisotropy of layered stratum than kriging and IDW. The IDW method is simple and has a wide applicability, as there is no need to adjust the method according to the characteristics of the data. The accuracy is acceptable when the sample points are adequate; however, the method does not consider irregular or heterogeneous fields of attributes. Linear triangular prism interpolation achieves complete $C_0$ continuity in a single triangular prism and offers computational simplicity. In addition, this interpolation method expresses local property changes and is less precise than kriging when considering the overall trend of change in a property field.

6. Conclusions

Most current 3DGM techniques focus on geometric and visualized models to explore geological data. A geometric model however, can only satisfy basic requirements. Obtaining an accurate and precise model of property fields in a geological body permits more advanced exploration for greater insights into a phenomenon of interest. Most of the interpolation methods generate large errors, as they apply general 3D statistical interpolation methods, such as kriging and IDW, and ignore the property characteristics in a geological model.

This paper compensates for the shortcomings of the traditional interpolation methods algorithmically, and proposes a function that expresses the heterogeneous distribution of properties in geological bodies: a spatial interpolation function based on the GTP. With this idea, we can build both a vector model and a geometric model. In addition, we can also express the internal change in the properties of a geological body. The experimental profile results, statistical errors, and cross-validation are used to analyze the accuracy of interpolation methods.

The experiments undertaken in this study show that the triangular prism quadric interpolation algorithm is more than comparable to kriging and IDW [26,37,38], and is actually superior. The spline elements of the 15 nodes on the triangular prism show a high accuracy and stability and are insensitive to mesh distortion. In terms of its application, it reflects both overall trends and local change in a property field, achieving smoothness when reaching an adjacent body. In addition, the triangular prism quadric interpolation algorithm also has directivity, expressing properties along both the horizontal and vertical layers. Some other interpolation algorithms, such as the 13-node pyramid element and the 21-node hexahedral element, have been successfully applied in elasticity mechanics [21]. The linear triangular prism interpolation algorithm does not perform as well as the interpolation algorithm tested; however, it is easy to use, has a high efficiency, and expresses local changes in properties. Kriging is a traditional low-pass filter interpolation algorithm [26], and expresses a property better than the linear triangular prism overall; however, this approach cannot reflect local information found in the original data. IDW interpolation has a “cluster” effect with no directivity; therefore, it gives the least satisfying results.

While we argue that the triangular prism interpolation algorithm is an effective and efficient method for expressing the heterogeneity of a property field, its accuracy and reliability have nevertheless not yet been fully established when considering restricted conditions such as folds, faults, and so on. Furthermore, in this paper, we achieve complete $C_0$ and $C_1$ continuity in a single triangular prism, but the splicing method for the adjacent triangular prisms yields only an approximate $C_1$ continuity. These issues will be addressed in our future research.