Non-linear inversion using general measures of data misfit and model structure

We investigate the use of general, non- l ₂ measures of data misfit and model structure in the solution of the non-linear inverse problem. Of particular interest are robust measures of data misfit, and measures of model structure which enable piecewise-constant models to be constructed. General measures can be incorporated into traditional linearized, iterative solutions to the non-linear problem through the use of an iteratively reweighted least-squares (IRLS) algorithm. We show how such an algorithm can be used to solve the linear inverse problem when general measures of misfit and structure are considered. The magnetic stripe example of Parker (1994 ) is used as an illustration. This example also emphasizes the benefits of using a robust measure of misfit when outliers are present in the data. We then show how the IRLS algorithm can be used within a linearized, iterative solution to the non-linear problem. The relevant procedure contains two iterative loops which can be combined in a number of ways. We present two possibilities. The first involves a line search to determine the most appropriate value of the trade-off parameter and the complete solution, via the IRLS algorithm, of the linearized inverse problem for each value of the trade-off parameter. In the second approach, a schedule of prescribed values for the trade-off parameter is used and the iterations required by the IRLS algorithm are combined with those for the linearized, iterative inversion procedure. These two variations are then applied to the 1-D inversion of both synthetic and field time-domain electromagnetic data.     

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     

Time-domain electromagnetic migration in the solution of inverse problems

Time-domain electromagnetic (TDEM) migration is based on downward extrapolation of the observed field in reverse time. In fact, the migrated EM field is the solution of the boundary-value problem for the adjoint Maxwell's equations. The important question is how this imaging technique can be related to the solution of the geoelectrical inverse problem. In this paper we introduce a new formulation of the inverse problem, based on the minimization of the residual-field energy flow through the surface or profile of observations. We demonstrate that TDEM migration can be interpreted as the first step in the solution of this specially formulated TDEM inverse problem. However, in many practical situations this first step produces a very efficient approximation to the geoelectrical model, which makes electromagnetic migration so attractive for practical applications. We demonstrate the effectiveness of this approach in inverting synthetic and practical TDEM data.     

Reweighting strategies in seismic deconvolution

Reweighting strategies have been widely used to diminish the influence of outliers in inverse problems. In a similar fashion, they can be used to design the regularization term that must be incorporated to solve an inverse problem successfully. Zero-order quadratic regularization, or damped least squares (pre-whitening) is a common procedure used to regularize the deconvolution problem. This procedure entails the definition of a constant damping term which is used to control the roughness of the deconvolved trace. In this paper I examine two different regularization criteria that lead to an algorithm where the damping term is adapted to successfully retrieve a broad-band reflectivity.
Synthetic and field data examples are used to illustrate the ability of the algorithm to deconvolve seismic traces.     

Bayesian inversion with Markov chains—I. The magnetotelluricone-dimensional case

We use Monte Carlo Markov chains to solve the Bayesian MT inverse problem in layered situations. The domain under study is divided into homogeneous layers, and the model parameters are the conductivity of each layer. We use an a priori distribution of the parameters which favours smooth models. For each layer, the a priori and a posteriori distributions are digitized over a limited set of conductivity values.
  The Markov chain relies on updating the model parameters during successive scanning of the domain under study. For each step of the scanning, the conductivity is updated in one layer given the actual value of the conductivity in the other layers. Thus we designed an ergodic Markov chain, the invariant distribution of which is the a posteriori distribution of the parameters, provided the forward problem is completely solved at each step.
  We have estimated the a posteriori marginal probability distributions from the simulated successive values of the Markov chain. In addition, we give examples of complex magnetotelluric impedance inversion in tabular situations, for both synthetic models and field situations, and discuss the influence of the smoothing parameter.     

Generalized Born inversion of seismic reflection data

Summary . Born inverse methods give accurate and stable results when the source wavelet is impulsive. However, in many practical applications (reflection seismology) an impulsive source cannot be realized and the inversion needs to be generalized to include an arbitrary source function. In this paper, we present a Born solution to the seismic inverse problem which can accommodate an arbitrary source function and give accurate and stable results. It is shown that the form of the generalized inversion algorithm reduces to a Wiener shaping ***filter, which is solved efficiently using a Levinson recursion algorithm. Numerical examples of synthetic and real field data illustrate the validity of our method.     

Electromagnetic 2.5-dimensional forward modelling with a boundary integral formulation

An effective and accurate technique for the numerical solution of 2-D electromagnetic scattering problems with 3-D sources is presented. This solution introduces a set of the usual boundary integral equations and uses a scalar Green's function. In this scalar version, the unknowns of the problem are the boundary values of the longitudinal fields and their normal derivatives in the Fourier domain. A generalization of the usual boundary integral formulation enables us to handle a large class of models composed of piecewise homogeneous domains, including contiguous domains, multiply-connected domains and unbounded domains. This formulation involves the solution of a system of linear equations, and results in a significant saving in computation time in comparison with other rigorous methods.
  The requirements for the numerical implementation of this solution are described in detail. Numerical tests were carried out using the important example of electromagnetic tomography. The specific symmetry properties of the response function in this case are illustrated. Numerical accuracy is verified over a large frequency range, up to 1  MHz.     

Three-dimensional inversion without blocks

Summary. We propose a method for solving non-linear inverse problems in the case where the unknown is a function of the spatial coordinates and the data set is discrete and finite. The method is based on a generalized leastsquares criterion, it is defined directly for non-linear problems (without previous linearization of the forward problem), and in the particular linear case it gives the same results (although slightly more general) than the Backus & Gilbert approach. As an example, we apply the method to the three-dimensional seismic velocity inverse problem, using as data the arrival times of seismic waves. The following paper (Nercessian et al .) shows some esults obtained using the present method.     

An algebraic formulation of geophysical inverse problems

Many geophysical inverse problems derive from governing partial differential equations with unknown coefficients. Alternatively, inverse problems often arise from integral equations associated with a Green's function solution to a governing differential equation. In their discrete form such equations reduce to systems of polynomial equations, known as algebraic equations. Using techniques from computational algebra one can address questions of the existence of solutions to such equations as well as the uniqueness of the solutions. The techniques are enumerative and exhaustive, requiring a finite number of computer operations. For example, calculating a bound to the total number of solutions reduces to computing the dimension of a linear vector space. The solution set itself may be constructed through the solution of an eigenvalue problem. The techniques are applied to a set of synthetic magnetotelluric values generated by conductivity variations within a layer. We find that the estimation of the conductivity and the electric field in the subsurface, based upon single-frequency magnetotelluric field values, is equivalent to a linear inverse problem. The techniques are also illustrated by an application to a magnetotelluric data set gathered at Battle Mountain, Nevada. Surface observations of the electric ( E _y ) and magnetic ( H _x ) fields are used to construct a model of subsurface electrical structure. Using techniques for algebraic equations it is shown that solutions exist, and that the set of solutions is finite. The total number of solutions is bounded above at 134 217 728. A numerical solution of the algebraic equations generates a conductivity structure in accordance with the current geological model for the area.     

Multichannel appraisal deconvolution

Summary The single channel scalar deconvolution method presented by Oldenburg has been extended to include N channels of data and vector models of the form     ( t ) = ( m ₁( t ), m ₂( t ), …, m _α( t ))^T. The solution has its foundation in the linear inverse theory of Backus & Gilbert and is effected by computing a set of N filters, which, when convolved with the data, yield unique averages of one of the scalar functions of the model. Those averages are the summation of the scalar model convolved with a primary averaging function plus contamination from secondary averaging functions convolved with other model components. It is shown how a set of suitably selected weights can annihilate these secondary averaging functions and thereby greatly simplify the interpretation. The computations are efficiently carried out in the frequency domain and require the inversion of an N × N Hermitian matrix at each frequency. As a type example, we have shown how the time varying elements of a seismic moment tensor might be computed from a set of seismograms.     

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.     

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.     

A general correspondence principle for time-domain electromagnetic wave and diffusion fields

A general correspondence principle is presented that relates any time-domain electromagnetic diffusion field to an electromagnetic wavefield in a 'corresponding' configuration. The principle applies to arbitrarily inhomogeneous and anisotropic media and arbitrary transmitters and receivers. For the correspondence between the two types of electromagnetic fields to hold, the electric conductivity in the diffusive case and the permittivity in the wavefield case should have the same spatial variation, while the permeability distributions in space in the two cases are to be identical. Essential steps in the derivation of the correspondence principle are the use of the time Laplace transformation of causal signals, taken at real, positive values of the transform parameter, the Schouten-Van der Pol theorem in the theory of the Laplace transformation, and the reliance upon Lerch's theorem of the uniqueness of the interrelation between causal field quantities and their time-Laplace-transform representations at real, positive values of the transform parameter. Correspondence is then established between the tensorial Green's functions in the two cases, where the Green's functions are the point-receiver responses (either electric or magnetic field) to point-transmitter excitations (either electric- or magnetic-current source).
Through the correspondence principle, all transient electromagnetic wavefields (where losses are neglected) have as a counterpart a transient diffusive electromagnetic field (where the electric displacement current is neglected). The interrelation yields the tool to compare quantitatively the potentialities of the two types of fields in transient electromagnetic geophysical prospecting.
Finally, a general medium-parameter scaling law for time-domain electromagnetic wavefields is presented.     

Modelling of electromagnetic fields in thin heterogeneous layers with application to field generation by volcanoes—theory and example

When interpreting electromagnetic fields observed at the Earth's surface in a realistic geophysical environment it is often necessary to pay special attention to the effects caused by inhomogeneities of the subsurface sedimentary and/or water layer and by inhomogeneities of the Earth's crust. The inhomogeneities of the Earth's crust are expected to be especially important when the electromagnetic field is generated by a source located in a magma chamber of a volcano. The simulation of such effects can be carried out using generalized thin-sheet models, which were independently introduced by Dmitriev (1969 ) and Ranganayaki & Madden (1980 ). In the first part of the paper, a system of integral equations is derived for the horizontal current that flows in the subsurface inhomogeneous conductive layer and for the vertical current crossing the inhomogeneous resistive layer representing the Earth's mantle. The terms relating to the finite thickness of the laterally inhomogeneous part of the model are retained in the equations. This only marginally complicates the equations, whilst allowing for a significant expansion of the approximation limits.
  The system of integral equations is solved using the iterative dissipative method developed by the authors in the period from 1978 to 1988. The method can be applied to the simulation of the electromagnetic field in an arbitrary inhomogeneous medium that dissipates the electromagnetic energy. When considered on a finite numerical grid, the integral equations are reduced to a system of linear equations that possess the same contraction properties as the original equations. As a result, the rate at which the iterative-perturbation sequence converges to the solution remains independent of the numerical grid used for the calculations. In contrast to previous publications on the method, aspects of the algorithm implementation that guarantee its effectiveness and robustness are discussed here.     

Numerical modelling and inversion of travel times of seismic body waves in inhomogeneous anisotropic media

Summary. Two approaches to travel-time computations in laterally inhomogeneous anisotropic media are presented. The first method is based on ray tracing in an anisotropic inhomogeneous medium, the second on the linearization procedure. The linearization procedure, which can be applied to inhomogeneous, slightly anisotropic media, does not require ray tracing in an anisotropic medium. Applications of linearized equations to the solutions of direct and inverse kinematic problems are discussed. A program package to perform the linearized computations for rather general 2-D laterally inhomogeneous layered structures is described and a numerical example is presented.     

Calculating the two-dimensional magnetotelluric Jacobian in finite elements using reciprocity

To speed up the calculation of the field Jacobian for 2-D magnetoteliuric inversion using finite elements, the principle of electromagnetic reciprocity is applied. The governing relationship for the Jacobian of the field along strike is obtained by differentiating the Helmholtz equation with respect to the resistivity of each region in the finite-element mesh. The result is a similar Helmholtz equation for the Jacobian, with new sources distributed over all nodes within the parameter medium. However, according to the principle of electromagnetic reciprocity, the roles of sources and receivers are interchangeable. Utilizing reciprocity, the field values obtained from the original forward problem and for new unit sources imposed at the receivers are then utilized in the calculation of the Jacobian by simple multiplication and summation with finite-element terms at each rectangle in the mesh. For the auxiliary (across-strike) fields, the Jacobian terms are obtained by solving source vectors loaded with parabola coefficients used in the approximation to Maxwell's equations. Jacobian terms for the apparent resistivity ( p _a), the impedance phase (φ) and the vertical magnetic field ( K _zy) are then calculated utilizing the parallel- and auxiliary-field Jacobians. Comparison of Jacobian values obtained from reciprocity calculations and by differencing two forward solutions show that the reciprocity method is accurate and can be used to decrease the number of calculations required to obtain sensitivities by one to two orders of magnitude.     

Interpretation of local two-dimensional electromagnetic anomalies by formalized trial procedure

Summary. An iterative algorithm is presented to be used in the search for the shape of a 2-D local deep geoelectric inhomogeneity lying in a layered medium; an anomaly having been identified in the usual way by observing an alternating time-harmonic electromagnetic field along the surface of the Earth. The normal section parameters (conductivities and thicknesses) and the excess electrical conductivity (inside inhomogeneity) are assumed to be known. The shape of the inhomogeneity is determined by means of a misfit functional minimization technique. A gradient minimization algorithm is constructed and Tikhonov's regularization scheme is applied to achieve stability of the solution. The effectiveness of such an approach is demonstrated by model calculations and by the interpretation of the Carpathian geomagnetic anomaly. Finally, a brief discussion of the problems of the practical application of this formalized trial procedure is presented. Because of the lack of reliable estimates of the excess conductivity, it is proposed to consider a family of models selected for the set of probable values of model parameters. This family can be treated as a generalized solution of the interpretation problem.     

Numerical inversion of deformation caused by pressure sources: application to Mount Etna (Italy)

The interpretation of geodetic data in volcanic areas is usually based on analytical deformation models. Although numerical finite element (FE) modelling allows realistic features such as topography and crustal heterogeneities to be included, the technique is not computationally convenient for solving inverse problems using classical methods. In this paper, we develop a general tool to perform inversions of geodetic data by means of 3-D FE models. The forward model is a library of numerical displacement solutions, where each entry of the library is the surface displacement due to a single stress component applied to an element of the grid. The final solution is a weighted combination of the six stress components applied to a single element-source. The pre-computed forward models are implemented in a global search algorithm, followed by an appraisal of the sampled solutions. After providing extended testing, we apply the method to model the 1993–1997 inflation phase at Mt Etna, documented by GPS and EDM measurements. We consider four different forward libraries, computed in models characterized by homogeneous/heterogeneous medium and flat/topographic free surface. Our results suggest that the elastic heterogeneities of the medium can significantly alter the position of the inferred source, while the topography has minor effect.     

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.