Fitting of the ground level is the most important step in multi-source geologic modeling. Whether the nodes are added by Delaunay subdivision method or by grid hierarchy, their elevation values need to be further solved. Therefore, interpolation methods must be used to accurately and smoothly process the initial ground level formed by limited data point information. Interpolation is the process of calculating unknown values based on known data, which results in the formation of a continuously distributed data field (Spragueetal., 2005). Numerous interpolation methods exist at present, but the interpolation methods suitable for the expression of 3D data field are still limited (Hu et al., 2007), the commonly used ones are Inverse-Distance Weighting (IDW) and Ordinary Kriging, among which the distance-weighted inverse-distance method belongs to a kind of deterministic difference method, Lu et al. ( 2008) extended the method, according to the number of samples and distribution density characteristics, so that it can determine the value of its parameters according to the characteristics of the samples; while the Kriging method belongs to a kind of deterministic difference method, Jessell (2001) carried out an in-depth study of it, based on which a kind of interpolation method of the potential-field (potential-field) is proposed, which is is able to handle discontinuous data fields in the presence of faults. Here, only these two methods are discussed.

5.3.3.1 Distance-weighted inverse method

The distance-weighted inverse method is one of the most commonly used methods for interpolating geologic data. It was first proposed by meteorologists and geologists and later improved by D. Shepard, so the method is called Shepard's method.

The basic idea of the inverse distance weighting method is to define the interpolation function f(P) as a weighted average of the function fk of each data point, which is considered to be the closest to the point to be interpolated with a number of known sampling points to the point to be interpolated with the largest contribution to the value of the point to be interpolated, and its contribution is inversely proportional to a certain power of the distance.

The basic principle of the inverse distance weighting method can be expressed by the following equation:

Research on three-dimensional geologic modeling and application of digital subspace and engineering

In which: f(P) is the estimated value of the point to be interpolated, P; fi is the sample value of the ith (i=1,..., n) known sampling point, Pi; di is the generalized distance between the ith sample point, Pi, and the point to be interpolated. interpolated point P; v is the power of the distance, which significantly affects the result of interpolation, and its selection criterion is the minimum mean absolute error. The results of related studies show that the higher the power, the more smoothing effect the interpolation result has (Lu et al., 2008). At a power index of 2, not only can the interpolation result be more satisfactory, but also has the advantage of being easy to calculate. In practical applications, the inverse distance squared method is generally used to find the value to be estimated.

The traditional inverse distance weighting method is very simple and natural, but when it is applied to the actual interpolation of geological features, it has the following obvious defects:

(1) When the number of data points is very large, the computation of f(P) will become very huge, and the huge amount of computation may even lead to the method becoming unachievable.

(2) The method only considers the distance from Pi to P, but not its direction. In fact, it is not sufficient to consider only the magnitude of the distance; there are known discrete points that are closer to the point to be interpolated, but its influence on the point to be interpolated may be shielded by other points.

(3) In the neighborhood of the known sampling point Pi, the error calculated due to di ≈ 0 will become very sensitive, especially when the two term forms are dominant and opposite in sign.

Therefore, the traditional inverse distance weighting method must be improved in practical application. Considering the spatial correlation of geologic feature data, two constraints of geologic body structure and influence distance can be attached to the inverse distance weighting method in order to improve the rationality and accuracy of spatial geometric or attribute data interpolation. The specific improvements are as follows:

(1) Geological body geometry constraints. In the process of spatial feature interpolation, only samples within the same stratigraphic lithology geology are selected when choosing the original sample data affecting the value to be estimated; samples that are not within the same stratigraphic lithology geology are not used even if they are very close. In other words, the spatial unit is divided according to the geologic structure, and the sample points are selected only within the same stratigraphic lithologic geology.

(2) Selection conditions for neighboring sample points. When choosing the neighboring points to be interpolated, three principles can be considered: one is the distance principle, i.e., according to the characteristics of the geological attribute data, a distance r is given, and the sample points outside the distance have no influence on the estimation of the points to be interpolated; the second is the number of points principle, i.e., given a data m, the m nearest sample points to be interpolated are used for the estimation; and the third is the use of the Voronoi diagram to find out the neighboring points of the points to be interpolated.

(3) Establish a data point index table to improve the search efficiency of the sample points around the to-be-interpolated point, thus greatly reducing the amount of computation of the large data volume team.

5.3.3.2 Ordinary Kriging Interpolation Method

Let the study area be A, and the regionalization variable (i.e., the physical attribute variable to be studied) be {Z(x)∈A}, with x denoting the spatial location. the attribute value of Z(x) at the sample point xi (i=1, 2, ..., n) (or the primary realization of the regionalization variable ) is Z(xi)(i=1, 2, ..., n), then according to the principle of ordinary Kriging interpolation, the estimated value of the attribute value Z(x0) at the unsampled point x0 is the weighted sum of the attribute values of the n known sampling points, i.e.

Research on 3D Geological Modeling and Application of Digital Subsurface and Engineering

Wherein, λi(i=1, 2 , ..., n) is the weighting coefficient to be sought.

It is assumed that the regionalized variable Z(x) satisfies the second-order smoothness assumption throughout the study area:

(1) The mathematical expectation of Z(x) exists and is equal to a constant: E[Z(x)] = m (constant).

(2) The covariance Cov(xi, xj) of Z(x) exists and is related only to the relative position between two points. Or the eigenhypothesis is satisfied:

(3) E[Z(xi)-Z(xj)]=0.

(4) The incremental covariance exists and is smooth: Var[Z(xi)-Z(xj)]=E[Z(xi)-Z(xj)]2.

Based on the requirement of unbiasedness: E[Z*(x0)]=E[Z(x0)].

The derivation leads to .

Minimize the estimation variance under the unbiased condition, i.e.

min{ Var[Z*(x0)-2μ( (λi-1))}, where μ is the Lagrange multiplier.

A system of equations for solving the weight coefficients λii(=1, 2, ..., n) can be obtained:

Research on 3D Geological Modeling and Application of Digital Subsurface Space and Engineering

After solving the plurality of the weight coefficients λii(=1, 2, ..., n), then the value of the attribute at the sampling point x0, Z*( x0).

The covariance Cov(xi, xj) in the above system of equations for solving the weight coefficients λii(=1, 2, ..., n), if expressed by the variational function γ(xi, xj), is of the form:

Research on three-dimensional geologic modeling and application of digital subsurface space and engineering

The definition of the variational function is as follows. Digital Subsurface Space and Engineering 3D Geological Modeling and Applied Research

The variance obtained by Kringing interpolation is:

Digital Subsurface Space and Engineering 3D Geological Modeling and Applied Research

or

Digital Subsurface Space and Engineering 3D Geological Modeling and Applied Research

A. Kringing interpolation in the Variational function

Variational function is the basis of Kringing interpolation method interpolation. In interpolation it is necessary to first determine the variance function of the regionalized variable under study. Assuming that the area under study is A, and there is a regionalized variable Z(x) in area A, and its one time sampling at position xi(i=1, 2,..., N) is Z(xi)(i=1, 2,..., N), the variational function of Z(x) is defined as:

Digital subsurface space and engineering three-dimensional geological modeling and Applied Research

The spatial variability of a spatial variable is the property of how this variable varies with location in space. The variability function reflects the spatial variability from different perspectives through its own structure and its parameters, and the process of determining the variability function is a process of analyzing the structure of spatial variability.

Let h be a vector with mode r=|h| and direction a. If there exist Nh pairs of observation data points separated by the vector h, the experimental variability function γ*(h) corresponding to the vector h in the direction of a can be expressed in the following form:

Research on three-dimensional geologic modeling and application of digital subsurface space and engineering

In which, z(xi+h) and z( xi) are the observations located on points xi+h and xi(i=1, 2, ..., Nh), respectively.

B. Theoretical model of the variational function

After obtaining the values of the experimental variational function, it is necessary to select the theoretical model of the variational function, and then fit the parameters of the selected theoretical model of the variational function, which is called "structural analysis".

The parameters of the theoretical model of the variational function generally include: range (generally denoted by a), abutment (sill, generally denoted by C(0)), arch height (generally denoted by C), and nugget constant (generally denoted by C0), as shown in Figure 5.9.

Figure 5.9 Theoretical variational function

Variational range a represents the demarcation line of the transition from the state of spatial correlation (|h| < a) to the state of non-existence of correlation (|h|>a); the discontinuity of the variational function at the origin is known as the Nugget effect, and the corresponding constant, C0= (h), is known as Nugget's constant; the abutment, C(0), has a covariance function: C(h) = C(0) - γ (h) the a priori variance of the second-order smooth regionalized variable Z(x): Var{Z(x)}=C(0)=γ(h); the arch height C is the difference between the abutment C(0) and the nugget of gold constant C0 in the variational function: C=C(0)-C0.

C. Classification of theoretical models of variational functions

The theoretical models of variational functions are generally divided into two major categories: those with abutment values and those without abutment values. Variational function theory models with abutment values include spherical models, exponential models, Gaussian models, etc. (Figure 5.10). The most commonly used is the spherical model. The formula of the spherical model is:

Research on three-dimensional geological modeling and application of digital underground space and engineering

Figure 5.10 Spherical model, exponential model, Gaussian model in the variability function

D. Reflection of the variability function on the structure of spatial variability

The variability function, as a kind of statistical tool for quantitative description of spatial variability, reflects spatial variability from different perspectives through its own structure and its parameters, it reflects the structure of spatial variability from different angles. Using the variability function, the continuity of spatial variables, correlation, influence range of variables, scale effect, discontinuity at the origin, anisotropy and other elements can be described.

E. Variational function theoretical model parameter fitting

Variational function theoretical model parameter fitting is the use of the original sampling point data or experimental variational function values of the selected theoretical model parameters are estimated in a specific way. The fitting method is generally hand-fitting.

Manual fitting is based on the experimental variance function values, on the one hand, through the observation of the experimental variance function graph; on the other hand, the necessary analysis of the regionalized variables under study, the use of visual observation to determine the variance function model parameters, and repeated cross-validation of the parameters, and ultimately determine the model parameters. The general process of its fitting is as follows:

(1) Firstly, the necessary analyses of the regionalized variables under study are carried out in terms of structure, background, etc., and the theoretical model of the variance function is determined by combining the experience of experts.

(2) The experimental variational function scatterplot is utilized to determine the block gold constants, the abutment values, the variational ranges, the anisotropy angles, and the anisotropy ratios in the variational function parameters.

(3) Cross-validation.