Spatial interpolation

Global interpolation method

1) Boundary interpolation method

The boundary interpolation method assumes that any significant changes occur at the boundary, and the changes within the boundary are uniform and homogeneous, ie, the same in all directions. This conceptual model is often used for soil and landscape mapping, and can be used to define other soil and landscape features by defining “homogeneous” soil units and landscape maps.

The simplest statistical model for boundary interpolation is the standard variance analysis (ANOVAR) model.

In the formula, Z is in x.₀ The attribute value of the location, Mu is the total average value, alpha_k it is the difference between the average value of class k and the average value of class mu, and the error of class E is the average error between classes (noise).

The model assumes that the attribute values of each class k are normally distributed; the mean of each class k (μ+α_k) is estimated from an independent sample set and is assumed to be spatially independent; the average error ε between classes assumes that all classes are identical.

The index to evaluate the classification effect is that the smaller the ratio is, the better the classification effect is, for the variance between classes and for the total variance. F test can be used to test the significance of classification results.

In essence, the theoretical assumption of boundary interpolation method is:

• Attribute value Z varies randomly in “patches” or landscape units, but not regularly.

• All “patches” of the same category have the same class variance (noise).

• All attribute values are normal distribution.

• All spatial variations occur at the boundary and are sudden rather than gradual.

When using boundary interpolation, careful consideration should be given to whether the data source conforms to these theoretical assumptions.

2) Trend surface analysis

The continuous change of some geographical attributes in space can be described by a smooth mathematical plane. The idea is to fit a smooth mathematical plane equation with the data of known sampling points, and then calculate the data of points without measurement value according to the equation, according to the relationship between the attribute data of sampling points and geographical coordinates, the method of smoothing mathematical plane equation by multiple regression analysis is called trend surface analysis. Its theoretical assumption is geographic coordinates.(x, y) It is an independent variable, the attribute value Z is also an independent variable and is normal distribution, and the regression error is also an independent variable independent of location.

Polynomial regression analysis is the simplest method to describe the characteristics of long-distance gradual change. The basic idea of polynomial regression is to use polynomials to represent lines and surfaces, and to fit data points according to the principle of least square method. The choice of line or surface polynomials depends on whether the data is one-dimensional or two-dimensional.

Use a simple example to illustrate that the eigenvalue z in a geographic or environmental survey is along a section at x₁, x₂…x_nSampling, if the z-value increases linearly with increasing x, the long-term variation of the eigenvalue can be calculated using one of the following regression equations:

Where b₀, b₁is the regression coefficient, and ε is the normal distribution residual (noise) independent of x.

However, in many cases, the quadratic polynomial is needed to fit the change in a more complex way rather than in a linear function.

For the two-dimensional case, the surface polynomials obtained by the multivariate regression analysis of XY coordinates are as follows:

The first three forms are:

Plane

Oblique plane

Two degree surface

Where P is the number of times of the trend surface equation.

It is the minimum number of terms of trend surface polynomials under normal conditions. Zero degree polynomial is a plane with one item number; first degree polynomial is an oblique plane with three items number; quadric surface has six items number and cubic trend surface has 10 items number.

The calculation factor b:sub:`i`is a standard multivariate regression problem. The advantages of trend surface analysis are very easy to understand, at least in terms of calculations. In most cases, low-order polynomials can be used for fitting, but it is more difficult to assign complex polynomial to explicit physical meaning.

The trend surface is a smoothing function, it is difficult to pass the original data point. Unless the data point is small and the trend surface number is high, the surface just passes through the original data point, so the trend surface analysis is an approximate interpolation method. In fact, the most effective application of trend surface is to reveal the largest deviation from the general trend in the region, so the main purpose of trend surface analysis is to use trend surface analysis to remove some macroscopic features from the data before using a local interpolation method, instead of directly using it for spatial interpolation.

As with multiple regression analysis, F distribution can be used to test the fitting degree of trend surface, and its test statistics are as follows:

Where U is the sum of the squares of the regression, Q is the sum of the squares of the residuals (the sum of the remaining squares), p is the number of polynomial terms (excluding the constant term b₀), and n is the number of data points used. When F>F_a, the trend surface is significant, otherwise it is not significant.

3) Transform function interpolation

The overall spatial interpolation based on empirical equations of one or more spatial parameters is also a frequently used spatial interpolation method, which is called a transformation function, the following is a description of a research example.

Heavy metal pollution in alluvial plains is related to several important factors, of which the distance from the source (river) and the elevation are the most important. Generally, the coarse sediment carrying heavy metals is deposited on the river beach, and the fine sediment carrying heavy metals is deposited in low-lying areas which are easy to be flooded during the flood season. Where the frequency of floods is low, the sediment particles carrying heavy metals are relatively small and are lightly polluted. Since the distance and elevation from the river are relatively easy to obtain spatial variables, the spatial interpolation between the contents of various heavy metals and their empirical equations can be used to improve the prediction of heavy metal pollution. The form of the regression equation in this case is as follows:

The formula is Z (x), a certain heavy metal content (ppm), b.₀ … B_n it’s a regression coefficient, p₁ … P_n is an independent spatial variable, in this case p₁ it is the distance factor from the river, p₂ it is an elevation factor.

This regression model is often called a conversion function and can be calculated by most GIS software. The transfer function can be applied to other independent variables such as temperature, elevation, rainfall, and distance from the sea and vegetation can be combined into a function of excess water content, the geographic location and its attributes can be combined into as many regressions as possible into the required regression model and then spatially interpolated. However, it should be noted that the physical meaning of the regression model must be clear, it should also be noted that all regression transfer functions belong to approximate spatial interpolation.

Global interpolation methods usually use standard statistical methods such as difference analysis and regression equation, and the calculation is relatively simple. Many other methods can also be used for global spatial interpolation, such as Fourier series and wavelet transform, especially for remote sensing image analysis, but they require a large amount of data.

Local interpolation method

Local interpolation only uses adjacent data points to estimate the value of unknown points, including several steps:

1. Define a neighborhood or search range;

2. Search for data points that fall within this neighborhood;

3. Choosing mathematical functions to express the spatial variations of these finite points;

4. assign values to data points that fall on regular grid cells. Repeat this step until all points on the grid are assigned.

The following aspects should be paid attention to when using local interpolation method: the interpolation function used; the size, shape and direction of neighborhood; the number of data points; whether the distribution of data points is regular or irregular.

1) Nearest Neighbor Method: Tyson Polygon Method

Thiessen polygon (also known as Dirichlet or Voronoi polygon) uses an extreme boundary interpolation method, using only the nearest single point for region interpolation. Tyson polygon divides the region into sub-regions according to the location of data points. Each sub-region contains a data point. The distance between each sub-region and its data points is less than any distance to other data points, and the data points in each sub-region are assigned. Connections connecting all data points form Delaunay triangle, which has the same topological structure as irregular triangle network TIN.

In GIS and geographic analysis, Tyson polygon is often used for fast assignment. In fact, an implicit assumption of Tyson polygon is that the nearest meteorological station is used for meteorological data at any location. In fact, unless there are enough meteorological stations, this assumption is inappropriate, because precipitation, air pressure, temperature and other phenomena are continuously changing. The results obtained by Tyson polygon interpolation method only occur on the boundary, and are homogeneous and unchanged within the boundary.

2) Moving average interpolation method: distance reciprocal interpolation

The distance reciprocal interpolation method combines the advantages of the adjacent point method of Tyson polygon and the gradient method of trend surface analysis, which assumes that the attribute value of the unknown point xsub:0is the distance weighted average of all data points in the local neighborhood. The distance reciprocal interpolation method is one of the weighted moving average methods, the weighted moving average method is calculated as follows:

In the formula, the weight coefficients are calculated by functions, requiring the reciprocal or negative exponential form at that time. The most common form is the reciprocal distance weighting function, which is as follows:

In the formula, x_j for unknown points, x_i for known data points.

The simplest form of the weighted moving average formula is called linear interpolation, fhe formula is as follows:

The distance reciprocal interpolation method is the most common method for GIS software to generate raster layers based on point data. The calculated value of the distance reciprocal method is susceptible to the cluster of data points. The calculation result often shows a “duck egg” distribution pattern with isolated point data significantly higher than the surrounding data points, which can be improved by dynamically modifying the search criteria during the interpolation process.

3) Interpolation method of spline function

Before the computer was used to fit curves and data points, the plotter used a flexible curve gauge to fit smooth curves piecewise. The piecewise curves drawn by this flexible curve plan are called splines. The data points matched with splines are called pile points, and the position of the control curve is controlled by the pile points when drawing the curve. The curve drawn by the curve gauge is mathematically described by a piecewise cubic polynomial function with continuous first-order and second-order continuous derivatives at its junction.

Spline function is a mathematical equation equivalent to flexible curve gauge in mathematics with a piecewise function. Only a few points can be fitted once, and at the same time, the continuity of the connection of curve segments can be ensured. This means that the spline function can modify the registration of a few data points without recalculating the whole curve. Trend surface analysis can not do this (Fig. 16).

Figure 8-16: Local Characteristics of Spline Functions

（a：当二次样条曲线的一个点位置变化时，只需要重新计算四段曲线；

B: When one point of a spline curve changes, only two segments of the curve need to be recalculated.

The general piecewise polynomial p(x)is defined as:

P (x) = psili (x)

I < x < x_i+1 (iTunes 1, 2, 3… , kmur1)

( j=0, 1, 2, …, r-1; i=1, 2, … , k-1 )

x:sub:`1`, …, x:sub:`k-1`The interval x:sub:`0`, x:sub:`k`Divided into ksubintervals, these split points are called “breakpoints”, and the points on the curve with these xvalues are called “sections”. The function p:sub:`i` (x)is a polynomial less than or equal to m times. The rentry is used to indicate the constraints of the spline function:

r=0, no constraint;

r=1, the function is continuous and has no constraints on its derivative;

r=m-1, interval [x:sub:`0`x:sub:`k`]can be represented by a polynomial;

R = m, with the most constraints.

The splines of m=1, 2, and 3 are the first, second, and third spline functions respectively, and the derivatives are 0th order, 1st order, and 2nd order derivatives. There must be one order at each node of the quadratic spline function. For continuous derivatives, each node of the cubic spline function must have a second-order continuous derivative at the beginning. * The simple spline of R = m* has only K + m* degrees of freedom, r = M = 3* has special significance because it is a cubic polynomial, which is called a spline function for the first time. The term “cubic spline” is used in a three-dimensional case where surface interpolation is performed instead of curve interpolation.

Since the range of discrete subintervals is wide, it may be a digitized curve, alculating simple splines in this range will cause certain mathematical problems, therefore, in practical applications, *B*splines are used—a special kind Spline function.

BSplines are often used for digitizing lines that are smoothed before display, such as soils, various boundaries on geological maps, and traditional cartography always wants to draw smoother curves. However, there are some problems with polygon boundary smoothing with Bsplines, especially the calculation of polygon area and perimeter, the result will be different from that before smoothing.

To sum up, spline function is piecewise function, only a few data points are used at a time, so the interpolation speed is fast. Compared with trend surface analysis and moving average method, spline function retains local variation characteristics. Both linear and surface spline functions are visually satisfactory. Some shortcomings of spline function are that the error of spline interpolation can not be estimated directly, at the same time, the problem to be solved in practice is the definition of spline blocks and how to assemble these “blocks” into complex surfaces in three-dimensional space without introducing any abnormal phenomena in the original surface.

4) Space autocovariance best interpolation method: Kriging interpolation

Several interpolation methods introduced in the previous section have not been well solved, such as the control parameters of trend surface analysis and the weights of reciprocal distance interpolation methods, which have a great impact on the results. These problems include:

• The number of average data points needs to be calculated.

• How to determine the neighborhood size, direction and shape of search data points?

• Is there a better method for estimating weight coefficients than calculating simple distance functions?

• Error problems related to interpolation.

To solve these problems, Georges Matheron, a French geomathematician, and D. G. Krige, a South African mining engineer, have studied an optimal interpolation method for mine exploration. This method is widely used in groundwater simulation, soil mapping and other fields, and has become an important part of GIS software geostatistical interpolation. This method fully absorbs the idea of geostatistics and considers that any attributes that change continuously in space are very irregular. It can not be simulated by simple smooth mathematical functions, but can be described appropriately by random surface. The spatial attributes of this continuity change are called “regional variables”, which can describe the index variables of image pressure, elevation and other continuity changes. This method of spatial interpolation using geostatistics is called Kriging interpolation. Geostatistical method provides an optimization strategy for spatial interpolation, i.e. the dynamic determination of variables according to some optimization criterion function in the interpolation process. The interpolation methods studied by Matheron, Krige and others focus on the determination of weight coefficients, so that the interpolation function is in the best state, that is, to provide the best linear unbiased estimation of the variable value at a given point.

The regional variable theory of Kriging interpolation assumes that the spatial variations of any variable can be expressed as the sum of the following three main components (Fig. 8-17):

1. Structural components related to constant mean or trend;

2. Random variables related to spatial change, i.e. regional variables;

3. Spatially independent random noise or residual error terms.

Figure 8-17: The theory of regional variables divides complex spatial changes into three parts

1. the average characteristics of terrain;

2. irregular spatial variations;

3. random and local variations”

   In addition, xis a position in 1D, 2D or 3D space. The value of the variable zat xcan be calculated as follows:

In the formula, m (x) is a deterministic function describing the structural components of Z (x); a random variable related to spatial variation, i.e. a regional variable; and a residual error term with zero mean value and spatial independent Gauss noise term.

The first step in the Kriging method is to determine the appropriate m(x)function. In the simplest case, m(x)is equal to the average of the sampling area, and the distance vector his separated. The mathematical expectation between two points x, x+his equal to zero:

Where z(x), z(x+h) is the value of the random variable zat x, x+h, and also assumes that the variance between the two points is only Related to distance h, then:

In the formula, it is a kind of function, called semi-variance function.

The two inherent assumptions of regional variable theory are the stability and variability of differences. Once the structural components are determined, the remaining differences belong to homogeneous changes, and the differences between different locations are only functions of distance. In this way, the formula for calculating regional variables can be written in the following form:

The estimating formulas of semi-variance are as follows:

In the formula, nis the number of sample point pairs with distance h(nto point), and the sampling interval his also called delay. The map corresponding to his called the “half variance graph”. Figures 8-18 show a typical half-variance plot.

Semi-variance is the first step in quantitatively describing regional changes, which provides useful information for spatial interpolation and optimization of sampling schemes. In order to obtain the semi-variance graph, we must first obtain the theoretical model of fitting the semi-variance, in the semi-variance model, we can get the following results:

1. The part of the curve with the larger value of delay his horizontal, the horizontal part of the curve becomes “Sill”. Note that there is no spatial correlation of data points in this range of delays because all variances do not change with distance.

2. The delay range of curve from low value to beam is called Range. The range is the most important part of the semivariogram because it describes how space-related differences vary with distance. The closer the distance is within the range of the range, the more similar features are. The range provides a way to determine the size of the window for the moving weighted averaging method. Obviously, the distance between the data point and the unknown point is greater than the range of the range, indicating that the data point is too far away from the unknown point and has no effect on the interpolation.

3. The fitting model in the graph does not pass through the origin, but is intersected with the coordinate axis in the positive direction.

Figure 8-18: Semi-variance Diagram

According to the semi-variance formula, when h = 0, it must be zero. The positive value appearing in the model is the estimated value of residual error, which is space-independent noise. Called “Nugget” variance, it is a combination of observation errors and spatial variations with small distance intervals.

When there are obvious ranges and beams, and the kernel variance is also important, but the numerical value is not too large, the spherical model can be used to fit the semi-variance (Figure 8-19-1). The formula is:

0 < h < a

h >= a

h=0

In the formula, A is a variable, h is a delay, c₀ +c₁ is the beam. Generally, the fitting effect of spherical model is ideal.

The exponential model (Figure 8-19-2) can be used for fitting if there are obvious kernel variances and beams without gradual variation.

If the kernel variance is small relative to the random variation related to the spatial variation, it is better to use curved curves, such as the Gauss curve: (Figure 8-19-4)

If the spatial variation varies with the range, but there are no beams, the linear model (Figure 8-19-3) can be used to fit:

In formula B is the slope of the line. Linear models are also used when the range of variation is much larger than the desired interpolation range.

Previous discussions assume that the surface features are identical in all directions, but in many cases, spatial variations have obvious directivity, which requires the use of models with different parameters to fit semi-variance maps.

Figure 8-19: Variable variance diagrams

The important use of the fitted semi-variance graph is to determine the weight factor needed for local interpolation. The process of determination is similar to the weighted moving interpolation method, but it is calculated not by a fixed function, but by the statistical analysis principle of the semi-variance graph of the sampled data. Namely:

The selection of weights should be unbiased and the variance of the estimation is less than that of other linear combinations of observations.

The minimum variance is:

The minimum variance can be obtained only if the following formula is true:

式中，是z在采样点 x_i ,

In the formula, it is the semivariance of Z between sampling points x sub:i, x sub:j and sampling points sub: x sub: J *and unknown points sub:xand sub: 0 *which can be obtained from the semi-variance diagram of the fitted model. The amount of Lagrangian operator needed to calculate the minimum variance.

This method is called conventional Kriging interpolation. It is an exact interpolation model, and the interpolation or the best local mean is consistent with the data points. In cartography, regular grids with thinner sampling intervals are often used for interpolation, which can be converted into contour maps by the methods mentioned above. Similarly, the estimated error, also known as Kriging variance, can be used to map the reliability of interpolation results throughout the study area.

An example of Kriging interpolation: (Figure 8-20)

Figure 8-20: An example of Kriging interpolation

Calculate the weight of the kriging method of the unknown point I0 in the graph (x=5, y = 5), and use the spherical model to fit the semivariance, parameter c₀ = 2.5 , c:sub :1 = 7.5; range a = 10.0. The data of these 5 sampling points are:

The solution process is represented by a matrix as follows:

Where A is the semivariance matrix between the data points, b is the semivariogram vector between each data point and the unknown point, λis the weight coefficient vector to be calculated, Φis the solution equation Lagrange operator.

Firstly, the distance matrix between the five data points is calculated.

*I 1 2 3 4 5*
*I1 0.0 5.099 9.899 5.000 3.162*
*I2 5.099 0.0 6.325 3.606 4.472*
*I3 9.899 6.325 0.0 5.0 7.211*
*I4 5.0 3.606 5.0 0.0 2.236*
*I5 3.162 4.472 7.211 2.236 0.0 *

And the distance vectors between them and unknown points:

*I I0*
*I1 4.243*
*I2 2.808 *
*I3 5.657 *
*I4 1.0*
*I5 2.0*

By introducing the above values into the spherical model, the corresponding semi-variance (A, b):

*A= I 1 2 3 4 5 6*
*1 2.500 7.739 9.999 7.656 5.939 1.000*
*2 7.739 2.500 8.677 6.381 7.196 1.000*
*3 9.999 8.677 2.500 7.656 9.206 1.000*
*4 7.656 6.381 7.656 2.500 4.936 1.000*
*5 5.939 7.196 9.206 4.936 2.500 1.000*
*6 1.000 1.000 1.000 1.000 1.000 0.000*
*b = I 0 *
*1 7.151*
*2 5.597*
*3 8.815*
*4 3.621*
*5 4.720*
*6 1.000*

Note that the extra 6th and 6th columns are to ensure that the sum of the weights is 1.

By calculating the inverse matrix of A, the following results are obtained:

A:sup:`-1`= I 1 2 3 4 5 6

1 -0.172 0.050 0.022 -0.026 0.126 0.273

2 0.050 -0.167 0.032 0.077 0.007 0.207

3 0.022 0.032 -0.111 0.066 -0.010 0.357

4 -0.026 0.077 0.066 -0.307 0.190 0.030

5 0.126 0.007 -0.010 0.190 -0.313 0.134

6 0.273 0.207 0.357 0.003 0.134 -6.873

Thus, the weight λis:

I weight coefficient distance

1 0.0175 4.423

2 0.2281 2.828

3 -0.0891 ∑ = 1 5.657

4 0.6437 1.000

5 0.1998 2.000

6 0.1182 φ

After interpolating the unknown point, the numerical value is as follows:

z(I0) =0.0175*3 + 0.2281*4 - 0.0891*2 + 0.6437*4 + 0.1998*6

= 4.566

The estimated variance is:

σ_e_2 = [0.0175* 7.151+ 0.2281* 5.597 - 0.0891* 8.815

+ 0.6437 *3.621 + .1998 *4.720 ] +φ

= 3.890 + 0.1182

= 4.008

The purpose of the Kriging interpolation method is to provide a method for determining the optimal weight coefficient and to describe the error information. Since the interpolation value of the Kriging point model (conventional kriging model) is related to the capacity of the original sample, when the sample is small, the interpolation result graph with simple point conventional Kriging interpolation will have obvious unevenness. This shortcoming can be overcome by modifying the Kriging equation to estimate the average within the sub-block B,this method is called block Kriging interpolation, which is suitable for estimating the average of a given area test cell or interpolating a regular grid of a given grid size.

The mean value of Z in subblock B is:

In the formula,S_Bis the area of sub-block B. z(x) still uses the same estimation formula:

The minimum variance is still:

However, the conditions for its establishment are as follows:

The variance of block Kriging interpolation is usually smaller than that of point Kriging interpolation, so the smooth interpolation surface generated will not produce concave and convex phenomena of point model.