Geophysical inversion with a neighbourhood algorithm—I. Searching a parameter space

This paper presents a new derivative-free search method for finding models of acceptable data fit in a multidimensional parameter space. It falls into the same class of method as simulated annealing and genetic algorithms, which are commonly used for global optimization problems. The objective here is to find an ensemble of models that preferentially sample the good data-fitting regions of parameter space, rather than seeking a single optimal model. (A related paper deals with the quantitative appraisal of the ensemble.)
  The new search algorithm makes use of the geometrical constructs known as Voronoi cells to derive the search in parameter space. These are nearest neighbour regions defined under a suitable distance norm. The algorithm is conceptually simple, requires just two 'tuning parameters', and makes use of only the rank of a data fit criterion rather than the numerical value. In this way all difficulties associated with the scaling of a data misfit function are avoided, and any combination of data fit criteria can be used. It is also shown how Voronoi cells can be used to enhance any existing direct search algorithm, by intermittently replacing the forward modelling calculations with nearest neighbour calculations.
  The new direct search algorithm is illustrated with an application to a synthetic problem involving the inversion of receiver functions for crustal seismic structure. This is known to be a non-linear problem, where linearized inversion techniques suffer from a strong dependence on the starting solution. It is shown that the new algorithm produces a sophisticated type of 'self-adaptive' search behaviour, which to our knowledge has not been demonstrated in any previous technique of this kind.     

Three-dimensional massively parallel electromagnetic inversion—I. Theory

An iterative solution to the non-linear 3-D electromagnetic inverse problem is obtained by successive linearized model updates using the method of conjugate gradients. Full wave equation modelling for controlled sources is employed to compute model sensitivities and predicted data in the frequency domain with an efficient 3-D finite-difference algorithm. Necessity dictates that the inverse be underdetermined, since realistic reconstructions require the solution for tens of thousands of parameters. In addition, large-scale 3-D forward modelling is required and this can easily involve the solution of over several million electric field unknowns per solve. A massively parallel computing platform has therefore been utilized to obtain reasonable execution times, and results are given for the 1840-node Intel Paragon. The solution is demonstrated with a synthetic example with added Gaussian noise, where the data were produced from an integral equation forward-modelling code, and is different from the finite difference code embedded in the inversion algorithm     

Stable Iterative Methods for the Inversion of Geophysical Data

Summary. Interpretation of earth electrical measurements can often be assisted by inversion, which is a non-linear model-fitting problem in these cases. Iterative methods are normally used, and the solution is defined by' best fit'in the sense of generalized least-squares.
The inverse problems we describe are ill-posed. That is, small changes in the data can lead to large changes in both the solution and in the iterative process that finds the solution. Through an analysis of the problem, based on local linearization, we define a class of methods that stabilize the iteration, and provide a robust solution. These methods are seen as generalizations of the well-known Singular Value Truncation and Marquardt Methods of iterative inversion.
Here, and in a companion paper, we give examples illustrating the successful application of the method to ill-posed problems relating to the resistivity of the Earth.     

Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data

A tomographic inversion technique that inverts traveltimes to obtain a model of the subsurface in terms of velocities and interfaces is presented. It uses a combination of refraction, wide-angle reflection and normal-incidence data, it simultaneously inverts for velocities and interface depths, and it is able to quantify the errors and trade-offs in the final model. The technique uses an iterative linearized approach to the non-linear traveltime inversion problem. The subsurface is represented as a set of layers separated by interfaces, across which the velocity may be discontinuous. Within each layer the velocity varies in two dimensions and has a continuous first derivative. Rays are traced in this medium using a technique based on ray perturbation theory, and two-point ray tracing is avoided by interpolating the traveltimes to the receivers from a roughly equidistant fan of rays. The calculated traveltimes are inverted by simultaneously minimizing the misfit between the data and calculated traveltimes, and the roughness of the model. This 'smoothing regularization' stabilizes the solution of the inverse problem. In practice, the first iterations are performed with a high level of smoothing. As the inversion proceeds, the level of smoothing is gradually reduced until the traveltime residual is at the estimated level of noise in the data. At this point, a minimum-feature solution is obtained, which should contain only those features discernible over the noise.
The technique is tested on a synthetic data set, demonstrating its accuracy and stability and also illustrating the desirability of including a large number of different ray types in an inversion.     

Some aspects of non-linearity in inversion

Summary. A direct attempt is made to examine the influence of non-linearity, in the dependence of data on parameters, on model resolution. An extension of the Marquardt approach to the non-linear case leads to bounded domains in discrete parameter space for a specified level of data misfit. The simplest class of non-linear influences on these domains may be simulated by modified linear terms via perturbation theory. Asymmetric non-linear effects are considered using an ensemble of attempts at the solution of the inverse problem.     

Approximate depth averages of electrical conductivity from surface magnetotelluric data

This paper presents a simple non-linear method of magnetotelluric inversion that accounts for the computation of depth averages of the electrical conductivity profile of the Earth. The method is not exact but it still preserves the non-linear character of the magnetotelluric inverse problem. The basic formula for the averages is derived from the well-known conductance equation, but instead of following the tradition of solving directly for conductivity, a solution is sought in terras of spatial averages of the conductivity distribution. Formulas for the variance and the resolution are then readily derived. In terms of Backus-Gilbert theory for linear appraisal, it is possible to inspect the classical trade-off curves between variance and resolution, but instead of resorting to linearized iterative methods the curves can be computed analytically. The stability of the averages naturally depends on their variance but this can be controlled at will. In general, the better the resolution the worse the variance. For the case of optimal resolution and worst variance, the formula for the averages reduces to the well-known Niblett-Bostick transformation. This explains why the transformation is unstable for noisy data. In this respect, the computation of averages leads naturally to a stable version of the Niblett-Bostick transformation. The performance of the method is illustrated with numerical experiments and applications to field data. These validate the formula as an approximate but useful tool for making inferences about the deep conductivity profile of the Earth, using no information or assumption other than the surface geophysical measurements.     

Non-linear altimetric geoid inversion for lithospheric elastic thickness and crustal density

The inversion of high-resolution geoid anomaly maps derived from satellite altimetry should allow one to retrieve the lithospheric elastic thickness, T _e , and crustal density, _c . Indeed, the bending of a lithospheric plate under the load of a seamount depends on both parameters, and the associated geoid anomaly is correspondingly dependent on the two parameters. The difference between the observed and modelled geoid signatures is estimated by a cost function, J , of the two variables, T _e and _c . We show that this cost function forms a valley structure along which many local minima appear, the global minimum of J corresponding to the true values of the lithospheric parameters. Classical gradient methods fail to find this global minimum because they converge to the first local minimum of J encountered, so that the final parameter estimate strongly depends on the starting pair of values ( T _e ,   _c ). We here implement a non-linear optimization algorithm to recover these two parameters from altimetry data. We demonstrate from the inversion of synthetic data that this approach ensures robust estimates of T _e and _c by activating two search phases alternately: a gradient phase to find a local minimum of J , and a tunnelling phase through high values of the cost function. The accuracy of the solution can be improved by a search in an iteratively restricted parameter subspace. Applying our non-linear inversion to the Great Meteor Seamount geoid data, we further show that the inverse problem is intrinsically ill-posed. As a consequence, minute geoid (or gravity) data errors can induce large changes in any recovery of lithospheric elastic thickness and crustal density.     

Groups, algebras, and the non-linearity of geophysical inverse problems

Mathematical methods from the theory of continuous groups are used to determine whether a non-linear inverse problem, in the form of a functional, can be transformed into a linear inverse problem. If such transformations exist they can be constructed from the solutions of a linear system of differential equations. An illustration of the methodology is given by the linearization of the functional relating basement topography to observed surface gravity. The linearized inversion of gravity data for basement topography is applied to observations from Yucca Mountain, Nevada. A 2.0 km step in the basement to the west of Yucca Mountain, corresponding to the Bare Mountain fault, matches the Bouguer gravity anomaly. The resolution and uncertainty associated with the estimates of basement topography indicate that the structure directly beneath the gravity line is well constrained.     

Error propagation in non-linear delay-time tomography

Delay-time tomography can be either linearized or non-linear. In the case of linearized tomography, an error due to the linearization is introduced. If the tomography is performed in a non-linear fashion, the theory used is more accurate from the physical point of view, but if the data have a statistical error, a noise bias in the model is introduced due to the non-linear propagation of errors. We investigate the error propagation of a weakly non-linear delay-time tomography example using second-order perturbation theory. This enables us to compare the linearization error with the noise bias. We show explicitly that the question of whether a non-linear inversion methods leads to a better estimation of the model parameters than a linearized method is dependent on the signal-to-noisc ratio. We also show that, in cases of poor data quality, a linearized inversion method leads to a better estimation of the model parameters than a non-linear method.     

Inversion of time domain three-dimensional electromagnetic data

We present a general formulation for inverting time domain electromagnetic data to recover a 3-D distribution of electrical conductivity. The forward problem is solved using finite volume methods in the spatial domain and an implicit method (Backward Euler) in the time domain. A modified Gauss–Newton strategy is employed to solve the inverse problem. The modifications include the use of a quasi-Newton method to generate a pre-conditioner for the perturbed system, and implementing an iterative Tikhonov approach in the solution to the inverse problem. In addition, we show how the size of the inverse problem can be reduced through a corrective source procedure. The same procedure can correct for discretization errors that inevidably arise. We also show how the inverse problem can be efficiently carried out even when the decay time for the conductor is significantly larger than the repetition time of the transmitter wave form. This requires a second processor to carry an additional forward modelling. Our inversion algorithm is general and is applicable for any electromagnetic field  ( E , H , d B / dt )  measured in the air, on the ground, or in boreholes, and from an arbitrary grounded or ungrounded source. Three synthetic examples illustrate the basic functionality of the algorithm, and a result from a field example shows applicability in a larger-scale field example.     

A general inverse method for modelling extensional sedimentary basins

A two-dimensional inverse model for extracting the spatial and temporal variation of strain rate from extensional sedimentary basins is presented and applied. This model is a generalization of a one-dimensional algorithm which minimizes the misfit between predicted and observed patterns of basin subsidence. Our calculations include the effects of two-dimensional conduction and advection of heat as well as flexural rigidity. More importantly, we make no prior assumptions about the duration, number or intensity of rifting periods. Instead, the distribution of strain rate is permitted to vary smoothly through space and time until the subsidence misfit has been minimized. We have applied this inversion algorithm to extensional sedimentary basins in a variety of geological settings. Basin stratigraphy can be accurately fitted and the resultant spatiotemporal distributions of strain rate are corroborated by independent information about the number and duration of rifting episodes. Perhaps surprisingly, the smallest misfits are achieved with flexural rigidities close to zero. Spatiotemporal strain rate distributions will help to constrain the dynamical evolution of thinning continental lithosphere. The strain rate pattern governs the heat-flow history and so two-dimensional inversion can be used to construct accurate maturation models. Finally, our inversion algorithm is a stepping stone towards a generalized three-dimensional implementation.     

Smooth inversion for ground surface temperature histories: estimating the optimum regularization parameter by generalized cross-validation

The inversion of recent borehole temperatures has proved to be a successful tool to determine ancient ground surface temperature histories. To take into account heterogeneity of thermal properties and their non-linear dependence on temperature itself, a versatile 1-D inversion technique based on a finite-difference approach has been developed. Regularization of the generally ill-posed problem is obtained by an appropriate version of Tikhonov regularization of variable order. In this approach, a regularization parameter has to be determined, representing a trade-off between data fit and model smoothness. We propose to select this parameter by generalized cross-validation. The resulting technique is employed in case studies from the Kola ultradeep drilling site, and another borehole from northeastern Poland. Comparing the results from both sites corroborates the hypothesis that subglacial ground surface temperatures as met in Kola often are much higher than the ones in areas exposed to atmospheric conditions (Poland).     

Three-dimensional magnetotelluric inversion using non-linear conjugate gradients

We have formulated a 3-D inverse solution for the magnetotelluric (MT) problem using the non-linear conjugate gradient method. Finite difference methods are used to compute predicted data efficiently and objective functional gradients. Only six forward modelling applications per frequency are typically required to produce the model update at each iteration. This efficiency is achieved by incorporating a simple line search procedure that calls for a sufficient reduction in the objective functional, instead of an exact determination of its minimum along a given descent direction. Additional efficiencies in the scheme are sought by incorporating preconditioning to accelerate solution convergence. Even with these efficiencies, the solution's realism and complexity are still limited by the speed and memory of serial processors. To overcome this barrier, the scheme has been implemented on a parallel computing platform where tens to thousands of processors operate on the problem simultaneously. The inversion scheme is tested by inverting data produced with a forward modelling code algorithmically different from that employed in the inversion algorithm. This check provides independent verification of the scheme since the two forward modelling algorithms are prone to different types of numerical error.     

New advances in three-dimensional controlled-source electromagnetic inversion

New techniques for improving both the computational and imaging performance of the three-dimensional (3-D) electromagnetic inverse problem are presented. A non-linear conjugate gradient algorithm is the framework of the inversion scheme. Full wave equation modelling for controlled sources is utilized for data simulation along with an efficient gradient computation approach for the model update. Improving the modelling efficiency of the 3-D finite difference (FD) method involves the separation of the potentially large modelling mesh, defining the set of model parameters, from the computational FD meshes used for field simulation. Grid spacings and thus overall grid sizes can be reduced and optimized according to source frequencies and source–receiver offsets of a given input data set. Further computational efficiency is obtained by combining different levels of parallelization. While the parallel scheme allows for an arbitrarily large number of parallel tasks, the relative amount of message passing is kept constant. Image enhancement is achieved by model parameter transformation functions, which enforce bounded conductivity parameters and thus prevent parameter overshoots. Further, a remedy for treating distorted data within the inversion process is presented. Data distortions simulated here include positioning errors and a highly conductive overburden, hiding the desired target signal. The methods are demonstrated using both synthetic and field data.     

Two efficient algorithms for iterative linearized inversion of seismic waveform data

We present two equivalent algorithms for iterative linearized waveform inversion for 3-D Earth structure with respect to an arbitrary 3-D starting model; one is a matrix formulation, and the second is a wavefield formulation. Both algorithms require the computation of accurate synthetic seismograms, but neither requires that any particular method be used to compute the synthetics. The matrix formulation is equivalent to our previously published algorithm (Hara, Tsuboi & Geller 1991), but requires less than 10 per cent of the CPU time of the previous algorithm. The wavefield algorithm is equivalent to that of Tarantola (1986) and Mora (1987), but appears to be substantially more efficient.     

A method of obtaining a velocity-depth envelope from wide-angle seismic data

Summary. Due to the non-uniqueness of traveltime inversion of seismic data, it is more appropriate to determine a velocity-depth ( v-z ) envelope, rather than just a v-z function. Several methods of obtaining a v-z envelope by extremal inversion have been proposed, all of which invert the data primarily from either x-p , or T-p , or both domains. These extremal inversion methods may be divided into two groups: linear extremal and non-linear extremal. There is some debate whether the linearized perturbation techniques should be applied to the inherently non-linear problem of traveltime inversion. We have obtained a v-z envelope by extremal inversion in T-p with the constraint that the inversion paths also satisfy x-p observations. Thus we use data jointly in r-p and x-p , and yet avoid the linearity assumptions.
This joint, non-linear extremal inversion method has been applied to obtain a v-z envelope down to a depth of about 30 km in the Baltimore Canyon trough using x-t data from an Expanding Spread Profile acquired during the LASE project. We have found that the area enclosed by the v-z envelope is reduced by about 15 per cent using x-p control on the T-p inversion paths, compared to the inversion without x-p control.     

Inversion for interface structure using teleseismic traveltime residuals

An inversion method is presented for the reconstruction of interface geometry between two or more crustal layers from teleseismic traveltime residuals. The method is applied to 2-D models consisting of continuous interfaces separating constant-velocity layers. The forward problem of determining ray paths and traveltimes between incident wave fronts below the structure and receivers located on the Earth's surface is solved by an efficient and robust shooting method. A conjugate gradient method is employed to solve the inverse problem of minimizing a least-squares type objective function based on the difference between observed and calculated traveltimes. Teleseismic data do not accurately constrain average vertical structure, so a priori information in the form of layer velocities and average layer thicknesses is required. Synthetic tests show that the method can be used to reconstruct interface geometry accurately, even in the presence of data noise. Tests also show that, if layer velocities and initial interface positions are poorly chosen, lateral structure is still recoverable. The inversion method was applied to previously published teleseismic data recorded by an in-line array of portable seismographs that traversed the northern margin of the Musgrave Block, central Australia. The solution based on interface parametrization is consistent with models given by other studies that used the same data but different methods, most notably the standard tomographic approach that inverts for velocity rather than interface structure.     

Resolution analysis for discrete systems

summary . Treatments of geophysical inverse problems have tended to polarize into approaches intended to generate models either described by piecewise continuous functions or with some prior discretization. The two approaches are here developed in parallel, and the ideas of a trade-off between the anticipated error and the attainable level of detail in the model estimate are extended to the discrete case, either with even or uneven discretization.
An alternative approach to specifying the potential resolution of a model is to establish upper and lower bounds on parameter values. Linear programming methods are extended to determine bounds which allow for subjective limits on parameter values.
For a non-linear system the possible resolution may be investigated by estimation procedures based on the full set of successful solutions obtained by Monte-Carlo inversion.     

Optimized and robust experimental design: a non-linear application to EM sounding

Pragmatic experimental design requires objective consideration of several classes of information including the survey goals, the range of expected Earth responses, acquisition costs, instrumental capabilities, experimental conditions and logistics. In this study we consider the ramifications of maximizing model parameter resolution through non-linear experimental design. Global optimization theory is employed to examine and rank different EM sounding survey designs in terms of model resolution as defined by linearized inverse theory. By studying both theoretically optimal and heuristic experimental survey configurations for various quantities of data, it is shown that design optimization is critical for minimizing model variance estimates, and is particularly important when the inverse problem becomes nearly underdetermined. We introduce the concept of robustness so that survey designs are relatively immune to the presence of potential bias errors in important data. Bias may arise during practical measurement, or from designing a survey using an appropriate model.