Fractal convergence properties of geophysical inversion

Abstract

Geophysical data, such as measurements of the Earth's gravitational or magnetic field, are routinely collected and studied to get information on the structure of the subsurface geology. The standard method of data analysis involves the least-squares inversion of the data with respect to various user-chosen model parameters (such as the geometry or density of the lithology). One serious problem with most inversion schemes is that they are liable to converge to local minima – that is they reach a set of model parameters whose geophysical response is a better fit to the observed data than the starting model, but they do not reach the set of parameters that would provide the best fit possible. The set of initial model parameters that converge to a particular minima of the misfit surface is studied here for some magnetic models, and is found to be a fractal when there are two minima available for the model to converge to, and at least two model parameters are inverted. The fractal dimension of the set is shown to be inversely proportional to the damping of the inversion process.

Introduction

Geophysical data, such as measurements of the strength of the magnetic or gravity field of the Earth, is collected routinely worldwide and is used to give information on the subsurface structure of the Earth. Common uses of geophysics include exploration for oil, gold, diamonds, groundwater, and other commodities. The Earth's magnetic field can be approximated by that of a dipole inclined at an angle of 11° from the axis of rotation and located at the centre of the Earth [1, p. 168]. Rocks within the crust that contain magnetic minerals possess induced magnetisation because of the presence of the geomagnetic field. Measurements of the strength of the magnetic field (generally made from aircraft or on the ground) therefore give information about the geology. Fig. 1 shows an example aeromagnetic dataset from South Africa.

Quantitative interpretation of data such as that shown in Fig. 1 is frequently performed by extracting profiles over features of interest and working with them. Geological structures possessing anomalous magnetisation are then approximated by polygonal shapes (see Fig. 2). This is known as 2D modelling, and is performed in preference to the 3D modelling necessary to model map or image data because the computational requirements are much less and because commercial 3D modelling software tends to be very expensive.

The geophysicist performing the interpretation will estimate some parameters of an initial geological model, such as its geometry, and compare the geophysical response of the model with the observed data (this is known as forward modelling). A process known as inversion is then generally used to modify the model such that its geophysical response becomes closer to the observed data. When the two match, the model is taken to be a possible representation of the subsurface geology. Due to ambiguities inherent in the geophysical techniques used, there are an infinite range of models that can be made to fit the data and it is normal for several different geophysical techniques to be applied in the same area, together with borehole and surface geological data, to constrain the possible solutions.

The inversion of most geophysical datasets is a difficult problem because of the non-linear response of the model to the parameters being inverted. This results in the inversion process being iterative, with the fit of the model response to the field data (hopefully) improving at each iteration. The change in the model parameters at each iteration can be obtained from the technique of Marquadt–Levenburg [2], [3], sometimes referred to as ridge regression [4];

$$\text{d}\text{P=(A}{}^{\mathrm{<mtext>T</mtext>}}\text{A+kI)}{}^{\mathrm{-1}}\text{A}{}^{\mathrm{<mtext>T</mtext>}}\text{e}$$

where dP is the m point parameter change matrix, A is termed the data kernel [5, p. 9], k is a constant and I is the identity matrix. e is the n point misfit matrix, i.e. e_i=O_i−C_i, where O_i is the ith observed data and C_i is the calculated response of the model at the same position. The factor k acts to damp the inversion, preventing the model parameters from assuming arbitrarily large values when A^TA is singular. It works by placing an upper limit on the parameter change vector. However, if the model response is a non-linear function of the parameter's values and hence the inversion process is iterative, then the greater the value of k, the slower the inversion will converge to the position of minimum misfit.