Geophysical parametrization and interpolation of irregular data using natural neighbours

Summary An approach is presented for interpolating a property of the Earth (for example temperature or seismic velocity) specified at a series of ‘reference’ points with arbitrary distribution in two or three dimensions. The method makes use of some powerful algorithms from the field of computational geometry to efficiently partition the medium into ‘Delaunay’ triangles (in 2-D) or tetrahedra (in 3-D) constructed around the irregularly spaced reference points. The field can then be smoothly interpolated anywhere in the medium using a method known as natural-neighbour interpolation. This method has the following useful properties: (1) the original function values are recovered exactly at the reference points; (2) the interpolation is entirely local (every point is only influenced by its natural-neighbour nodes); and (3) the derivatives of the interpolated function are continuous everywhere except at the reference points. In addition, the ability to handle highly irregular distributions of nodes means that large variations in the scale-lengths of the interpolated function can be represented easily. These properties make the procedure ideally suited for ‘gridding’ of irregularly spaced geophysical data, or as the basis of parametrization in inverse problems such as seismic tomography. We have extended the theory to produce expressions for the derivatives of the interpolated function. These may be calculated efficiently by modifying an existing algorithm which calculates the interpolated function using only local information. Full details of the theory and numerical algorithms are given. The new theory for function and derivative interpolation has applications to a range of geophysical interpolation and parametrization problems. In addition, it shows much promise when used as the basis of a finite-element procedure for numerical solution of partial differential equations.