Radial basis and spline functions

<< Click to Display Table of Contents >>

Navigation:  Surface and Field Analysis > Deterministic Interpolation Methods >

Radial basis and spline functions

Radial basis interpolation is the name given to a large family of exact interpolators. In many ways the methods applied are similar to those used in geostatistical interpolation, but without the benefit of prior analysis of variograms. On the other hand they do not make any assumptions regarding the input data points (other than they are not co-linear) and provide excellent interpolators for a wide range of data. See also Section 8.3.2, Radial basis function networks, for an alternative (neural network) approach to modeling with radial basis functions.

A great deal of research has been conducted into the quality of these interpolators, across many disciplines. For terrain modeling and Earth Sciences generally the so-called multiquadric function has been found particularly effective, as have thin plate splines. The simplest variant of this method, without smoothing (i.e. as an exact interpolator) can be viewed as a weighted linear function of distance (or inverse distance) from grid point to data point, plus a “bias” factor, m. The model is of the form:

or the equivalent model, using the untransformed data values and data weights λi:

In these expressions zp is the estimated value for the surface at grid point p; φ(ri) is the radial basis function selected, with ri being the radial distance from point p to the ith data point; the weights wi and λi and the bias value m (or Lagrangian multiplier) are estimated from the data points. This requires solving a system of n linear equations. Using the second of the two models above the procedure is then essentially the same as for Ordinary Kriging (see further, Section 6.7.2, Kriging interpolation). For clarity we outline this latter procedure below as a series of steps using matrix notation:

Compute the nxn matrix, D, of inter-point distances between all (x,y) pairs in the source dataset (or a selected subset of these)

Apply the chosen radial basis function, φ(), to each distance in D, to produce a new array Φ

Augment Φ with a unit column vector and a unit row vector, plus a single entry 0 in position (n+1),(n+1). Call this augmented matrix A — see below for an illustration

Compute the column vector r of distances from the grid point, p, to each of the source data points used to create D

Apply the chosen radial basis function to each distance in r, to produce a column vector φ and then create the (n+1) column vector c as φ plus a single 1 as the last entry

Compute the matrix product b=A-1c. This provides the set of n weights to be used in the calculation of the estimated value at p, plus the Lagrangian value m using the linear equation cited earlier in this Section. The value m is not used in calculating the estimate at p in this formulation

In matrix form the system of linear equations being solved is of the form:

i.e. Ab=c hence b=A-1c, assuming that A is non-singular.

A variety of different radial basis functions may be used. A selection of those commonly available follows. Where more than one functional form is cited by separate GIS software suppliers the name of the package has been included. The parameter c in these expressions (not to be confused with the vector c above), which may be 0, determines the amount of smoothing. For an example implementation using a number of alternative Radial Basis functions on track (transect-like) data see Carlson and Foley (1991, 1992). The GRASS function that performs bi-linear or bicubic spline interpolation (incorporated into QGIS, SEXTANTE plugin) incorporates a modification of standard spline interpolation by the addition of a regularization component. This addresses problems of instability with some problems, i.e. where the matrix A above is non-singular but ill-conditioned (such problems can result in very variable surface interpolations even when input data values only vary by a small amount). The procedure, known as Tikhonov regularization, results in smoothing of the interpolation, depending on the value of the smoothing parameter selected.


If the application of this function is restricted to some range, a, and we write h=r/a, then

is approximately linear over the range [c,1] and approximately c over the range [0,c].

Inverse multiquadric:

Thin plate spline:


Natural cubic spline:

Spline with tension:

where I0() is the modified Bessel function and γ=0.5771… is Euler’s constant. As noted earlier in this Guide (Table 1‑3) the modified Bessel function of order 0 is given by:

Completely regularized spline function:

where E1() is the exponential integral function and γ is Euler’s constant, as before. The exponential integral function is given by:

In each of the above models the smoothing parameter c remains to be determined. ArcGIS Geostatistical Analyst seeks to optimize c by computing the RMSE of the prediction versus the data point values using simple cross-validation. Surfer selects a default value based on the diagonal length of the MBR of the data points and the number of points this MBR contains. This may then be modified to meet the user’s requirements.