Geophysical Inverse Problems Research Articles

Recursive ant colony optimization: a new technique for the estimation of function parameters from geophysical field data

We introduce Recursive Ant Colony Optimization (RACO), a new technique for the estimation of function parameters from field data obtained from various geophysical surveys. RACO is a modified form of the ant colony method (ACO) and can be used to determine the best solution for various geophysical problems. As RACO is an extension of ACO, it simulates the social behaviour of ants, optimizing their path from the nest to the food source. RACO applies ACO recursively, introducing an additional term ‘depth’ that decides the extent of the recursion. The results of each depth contribute towards the construction of models for the following depth and the range of values for each parameter is condensed around the actual solution. The algorithm is tested on a two parameter mathematical function and to further test its efficiency and stability on real world problems, we applied it on a few geophysical inversion problems. We show its application on synthetic data sets for the self‐potential anomaly (noise free and 10% noise) caused in the case of an inclined sheet like a body buried inside the earth, synthetic 1D vertical electrical sounding (VES) geo‐electric data sets for three‐layer and five‐layer earth models, a real 1D VES data set obtained from the Tangasol region, India and P‐wave logs for two sites located in the Krishna‐Godavari (KG) Basin, India. The results generated by RACO inversions are found to be in good agreement with previous studies conducted using these data sets.

Bayesian inversion of marine CSEM data with a trans-dimensional self parametrizing algorithm

The posterior distribution of earth models that fit observed geophysical data convey information on the uncertainty with which they are resolved. From another perspective, the non-uniqueness inherent in most geophysical inverse problems of interest can be quantified by examining the posterior model distribution converged upon by a Bayesian inversion. In this work we apply a reversible jump Markov chain Monte Carlo method to sample the posterior model distribution for the anisotropic 1-D seafloor conductivity constrained by marine controlled source electromagnetic data. Unlike conventional gradient based inversion approaches, our algorithm does not require any subjective choice of regularization parameter, and it is self parametrizing and trans-dimensional in that the number of interfaces with a resistivity contrast at depth is variable, as are their positions. A synthetic example demonstrates how the algorithm can be used to appraise the resolution capabilities of various electromagnetic field components for mapping a thin resistive reservoir buried beneath anisotropic conductive sediments. A second example applies the method to survey data collected over the Pluto gas field on the Northwest Australian shelf. A benefit of our Bayesian approach is that subsets of the posterior model probabilities can be selected to test various hypotheses about the model structure, without requiring further inversions. As examples, the subset of model probabilities can be viewed for models only containing a certain number of layers, or for models where resistive layers are present between a certain interval as suggested by other geological constraints such as seismic stratigraphy or nearby well logs.

Efficient methods for large‐scale linear inversion using a geostatistical approach

In geophysical inverse problems, such as estimating the unknown parameter field from noisy observations of dependent quantities, e.g., hydraulic conductivity from head observations, stochastic Bayesian and geostatistical approaches are frequently used. To obtain best estimates and conditional realizations it is required to perform several matrix‐matrix computations involving the covariance matrix of the discretized field of the parameters. In realistic three‐dimensional fields that are finely discretized, these operations as performed in conventional algorithms become extremely expensive and even prohibitive in terms of memory and computational requirements. Using Hierarchical Matrices, we show how to reduce the complexity of forming approximate matrix‐vector products involving the Covariance matrices in log linear complexity for an arbitrary distribution of points and a wide variety of generalized covariance functions. The resulting system of equations is solved iteratively using a matrix‐free Krylov subspace approach. Furthermore, we show how to generate unconditional realizations using an approximation to the square root of the covariance matrix using Chebyshev matrix polynomials and use the above to generate conditional realizations. We demonstrate the efficiency of our method on a few standard test problems, such as interpolation from noisy observations and contaminant source identification.

The levenberg-marquardt method for approximation of solutions of irregular operator equations

An ill-posed problem is considered in the form of a nonlinear operator equation with a discontinuous inverse operator. It is known that in investigating a high convergence of the methods of the type of Levenberg-Marquardt (LM) method, one is forced to impose very severe constraints on the problem operator. In the suggested article the LM method convergence is set up not for the initial problem, but for the Tikhonov-regularized equation. This makes it possible to construct a stable Fejer algorithm for approximation of the solution of the initial irregular problem at the conventional, comparatively nonburdensome conditions on the operator. The developed method is tested on the solution of an inverse problem of geophysics.

Computational recipes for electromagnetic inverse problems

SUMMARY The Jacobian of the non-linear mapping from model parameters to observations is a key component in all gradient-based inversion methods, including variants on Gauss–Newton and non-linear conjugate gradients. Here, we develop a general mathematical framework for Jacobian computations arising in electromagnetic (EM) geophysical inverse problems. Our analysis, which is based on the discrete formulation of the forward problem, divides computations into components (data functionals, forward and adjoint solvers, model parameter mappings), and clarifies dependencies among these elements within realistic numerical inversion codes. To be concrete, we focus much of the specific discussion on 2-D and 3-D magnetotelluric (MT) inverse problems, but our analysis is applicable to a wide range of active and passive source EM methods. The general theory developed here provides the basis for development of a modular system of computer codes for inversion of EM geophysical data, which we summarize at the end of the paper.

A Study of the Three‐Term Non‐Gaussian Pre‐Stack Inversion Method

AbstractIn view of the wide use of the “Gaussian distribution hypothesis of noise” in geophysical inversion problems, we study the non‐Gaussian distribution characteristics of the noise in pre‐stack seismic data. Upon this analysis, the concept of non‐Gaussian inversion of pre‐stack three‐term inversion is put forward, in which the mixed‐norm is proposed for suppressing both the Gaussian and non‐Gaussian noises. In view of the non‐derivable feature of the objective function, the Powell algorithm is used to solve the objective function. Model test and real seismic data inversion results show the correctness and reliability of the algorithm.

Crosshole traveltime tomography using particle swarm optimization: A near-surface field example

Particle swarm optimization (PSO) is a relatively new global optimization approach inspired by the social behavior of bird flocking and fish schooling. Although this approach has proven to provide excellent convergence rates in different optimization problems, it has seldom been applied to inverse geophysical problems. Until today, published geophysical applications mainly focus on finding an optimum solution for simple, 1D inverse problems. We have applied PSO-based optimization strategies to reconstruct 2D P-wave velocity fields from crosshole traveltime data sets. Our inversion strategy also includes generating and analyzing a representative ensemble of acceptable models, which allows us to appraise uncertainty and nonuniqueness issues. The potential of our strategy was tested on field data collected at a well-constrained test site in Horstwalde, Germany. At this field site, the shallow subsurface mainly consists of sand- and gravel-dominated glaciofluvial sediments, which, as known from several boreholes and other geophysical experiments, exhibit some well-defined layering at the scale of our crosshole seismic data. Thus, we have implemented a flexible, layer-based model parameterization, which, compared with standard cell-based parameterizations, allows for significantly reducing the number of unknown model parameters and for efficiently implementing a priori model constraints. Comparing the 2D velocity fields resulting from our PSO strategy to independent borehole and direct-push data illustrated the benefits of choosing an efficient global optimization approach. These include a straightforward and understandable appraisal of nonuniqueness issues as well as the possibility of an improved and also more objective interpretation.

S- and R-approximations of the Earth's surface topography

The linear integral representation method has been developed from the Backus and Gilbert [G. Backus and F. Gilbert, Numerical application of formalism for geophysical inverse problems, Geophys. J. R. Astron. Soc. 13 (1967), pp. 247–276; G. Backus and F. Gilbert, The resolving power of gross Earth data, Geophys. J. R. Astron. Soc. 16 (1968), pp. 169–205] method for solving problems, in which the number of data is sufficiently large; at the same time the spatial distribution of some physical properties is represented by a large number of the parameters, or even is an unknown function of spatial coordinates. The linear integral representation method enables us to obtain the solution using discrete and approximately given data. The two versions of the linear integral representation method are applied to the determination of a function describing the Earth's surface topography. Results of numerical experiments and examples of a real topography are presented. The methods proposed are applicable to the solution of various geophysical, geomorphological and geodetic problems.

A Dual‐Parameter Regularization Method for Electrical Conductivity Imaging

AbstractThis paper presents a dual‐parameter regularization method with a term that has a second order regularization operator. The optimal value of the regularization parameter is determined by applying the L‐curve criterion, the discrepancy principle and generalized cross‐validation. The validity and superiority of the proposed method is verified by numerical simulations of the theoretical model with truncated singular value decomposition, conjugate gradient, and standard‐form Tikhonov regularizations. The results are of very high precision and stable against random data noise. Finally, we apply the proposed method to geophysical inversion problems in electrical conductivity imaging. The imaging results suggest that the proposed method improves the imaging quality and fidelity. This further validates that, when it is applied to geophysical inversion problems, the proposed method improves the stability of inversion, reduces the possible multiplicity of solutions, and enhances computation speed.

Efficient stochastic estimation of the model resolution matrix diagonal and generalized cross–validation for large geophysical inverse problems

[1] In recent years, larger geophysical data sets and novel model parameterizations have dramatically increased both the data and model space dimensions of many inverse problems. Because of their relatively low computational expense, trade–off curve corner estimation for choosing regularized models and “checkerboard” tests for evaluating model resolution are commonly applied, despite their limitations. We present and demonstrate a low–cost method for accurately estimating the diagonal elements of the model resolution matrix and for implementing generalized cross–validation (GCV) for optimal regularization parameter selection. The ability to estimate the diagonal of the resolution matrix and GCV function thus facilitates the introduction of additional tools for diagonal resolution analysis and regularization evaluation, even for very large inverse problems, with storage and computational costs comparable to those required for obtaining model solutions. We demonstrate the method using a Tikhonov regularized teleseismic body wave velocity inversion example with approximately 260,000 model parameters, where we validate randomly selected Rm diagonal elements against explicitly calculated values and compare GCV-estimated regularized model results to those obtained through traditional methods.

Frequency‐domain acoustic‐elastic coupled waveform inversion using the Gauss‐Newton conjugate gradient method

ABSTRACTWe developed a frequency‐domain acoustic‐elastic coupled waveform inversion based on the Gauss‐Newton conjugate gradient method. Despite the use of a high‐performance computer system and a state‐of‐the‐art parallel computation algorithm, it remained computationally prohibitive to calculate the approximate Hessian explicitly for a large‐scale inverse problem. Therefore, we adopted the conjugate gradient least‐squares algorithm, which is frequently used for geophysical inverse problems, to implement the Gauss‐Newton method so that the approximate Hessian is calculated implicitly. Thus, there was no need to store the Hessian matrix. By simultaneously back‐propagating multi‐components consisting of the pressure and displacements, we could efficiently extract information on the subsurface structures. To verify our algorithm, we applied it to synthetic data sets generated from the Marmousi‐2 model and the modified SEG/EAGE salt model. We also extended our algorithm to the ocean‐bottom cable environment and verified it using ocean‐bottom cable data generated from the Marmousi‐2 model. With the assumption of a hard seafloor, we recovered both the P‐wave velocity of complicated subsurface structures as well as the S‐wave velocity. Although the inversion of the S‐wave velocity is not feasible for the high Poisson's ratios used to simulate a soft seafloor, several strategies exist to treat this problem. Our example using multi‐component data showed some promise in mitigating the soft seafloor effect. However, this issue still remains open.

Fast solution of geophysical inversion using adaptive mesh, space-filling curves and wavelet compression

Summary Many geophysical inverse problems involve large and dense coefficient matrices that often exceed the limitations of physical memory in commonly available computers. The repeated multiplications of such matrices to vectors during processing or inversion require an immense amount of computing power. These two factors pose a significant challenge to solving large-scale inverse problems in practice and can render many realistic problems intractable. To overcome these limitations, we develop a new computational approach for this class of problems by combining an adaptive quadtree or octree model discretization and wavelet transforms on reordered parameter sets. The adaptive mesh discretizes the model region according to the required resolutions based on localized anomalies. Hilbert space-filling curves and similar ordering of the reduced parameter set then enable a higher compression of the coefficient matrix by forming its sparse representation in the 1-D wavelet domain. This combination can reduce the storage requirement by 100 to 1000 times and, therefore, also speeds up the computation during the processing stage by the same factor. As a result, problems can now be solved that were computationally prohibitive. We present the algorithm and illustrate its effectiveness with an example from equivalent source construction in potential-field processing.

Magnetic inverse modelling of a dike using the artificial neural network approach

ABSTRACTArtificial neural systems have been used in a variety of problems in the fields of science and engineering. Here we describe a study of the application of neural networks in solving some geophysical inverse problems. In particular, we try to estimate the depth of dikes using magnetic data and a three‐layer feed forward neural network. The network is trained by synthetic data as input and output. For forward neural network training we use the back‐propagation algorithm. Results indicate that forward neural networks, if adequately trained, can predict a reasonably accurate depth for dikes. The proposed method was applied to magnetic data over the Darmian Iron field in Iran. Results were compared to real values from well data and proved the good performance of the trained neural network in predicting the dike’s depth.

A MULTIPLICATIVE REGULARIZED GAUSS-NEWTON ALGORITHM AND ITS APPLICATION TO THE JOINT INVERSION OF INDUCTION LOGGING AND NEAR-BOREHOLE PRESSURE MEASUREMENTS

Due to the ill-posed nature of nonlinear inverse problems of borehole geophysics, a parameterization approach is necessary when the available measurement data are limited and measurements are only carried out from sparse transmitter-receiver positions (limited data diversity). A potential remedy is the joint inversion of multi- physics measurements. A parametric inversion approach has desirable attributes for multi-physics measurements with difierent resolutions. It provides a ∞exible framework to put the sensitivities of multi- physics multi-resolution measurements on equal footing. In addition, the number of unknown model parameters to be inverted is rendered tractable with parameterization. Consequently, a Gauss-Newton based inversion algorithm taking advantage of the Hessian information can be advantageously employed over inversion approaches that rely only on gradient information. We describe a new dual-physics parametric joint-inversion algorithm to estimate near-borehole ∞uid permeability and porosity distributions of rock formations from ∞uid-∞ow and electromagnetic measurements. In order to accommodate the cases in which the measurements are redundant or lack sensitivity with respect

The 2-D magnetotelluric inverse problem solved with optimization

SUMMARY The practical 2-D magnetotelluric inverse problem seeks to determine the shallow-Earth conductivity structure using finite and uncertain data collected on the ground surface. We present an approach based on using PLTMG (Piecewise Linear Triangular MultiGrid), a special-purpose code for optimization with second-order partial differential equation (PDE) constraints. At each frequency, the electromagnetic field and conductivity are treated as unknowns in an optimization problem in which the data misfit is minimized subject to constraints that include Maxwell's equations and the boundary conditions. Within this framework it is straightforward to accommodate upper and lower bounds or other conditions on the conductivity. In addition, as the underlying inverse problem is ill-posed, constraints may be used to apply various kinds of regularization. We discuss some of the advantages and difficulties associated with using PDE-constrained optimization as the basis for solving large-scale nonlinear geophysical inverse problems. Combined transverse electric and transverse magnetic complex admittances from the COPROD2 data are inverted. First, we invert penalizing size and roughness giving solutions that are similar to those found previously. In a second example, conventional regularization is replaced by a technique that imposes upper and lower bounds on the model. In both examples the data misfit is better than that obtained previously, without any increase in model complexity.

Method of some inverse geophysical problem solution by using weakly coupled connected distributed computing systems

In this work a possible model for organizing a grid-based application that performs the solution of several inverse geophysical problems is described. As an example, we consider the problem of determining the parameters of seismic anisotropy in the Earth's mantle by the inversion of seismic waveforms. It is shown that this class of problems is reduced to the tabulation of a complicated multidimensional function. In this approach, the calculation at each point in a definition interval is calculated independently, so this is ideally appropriate for calculations that use a loosely connected distributed computing infrastructure.

A unified transformation function for lower and upper bounding constraints on model parameters in electrical and electromagnetic inversion

Since many geophysical inverse problems are ill-posed, implementation of model parameter constraints is important in reducing solution ambiguity. This paper presents a unified transformation function for enforcing lower and upper bounds on model parameters to restrict model updates during the inversion process such that unrealistic results are suppressed. The transformation is realized using a logarithmic or inverse hyperbolic tangent function on model parameters, and contains the conventional transformations as special cases. It is straightforward to recast the cost functional gradient in terms of the transformed variable. The interval between lower and upper bounds reflects the reliability of a priori information, which may be obtained from well logging and surface geological surveys. Tests on a synthetic crosswell electromagnetic example reveal that the use of bound constraints is effective in deriving valid physical results from inversion.

Magnetotellurics and radio-wave interference sounding

The plane harmonic electromagnetic fields are considered in the theory of magnetotelluric methods in the range of frequencies from 0.0001 Hz to 20 kHz. These fields are natural by their origin and contain information on the depths from tens of meters up to 100 km and more. The magnetotelluric soundings, which use the fields of radio stations, expand the frequency band almost up to 1 MHz and make it possible to study the depths from the first few meters. The method of radio-wave interference sounding supplements geoelectric prospecting on plane waves into the range of even higher frequencies (up to 100 MHz). In this case, the conduction and displacement currents become comparable, which makes it possible to distinguish objects both by their electrical conductivity and by their dielectric permittivity. For the two-layered model of a medium, there exist simple kinematic methods of data interpretation of a radio-interferometry sounding. Within multilayer, and especially horizontally heterogeneous, media, methods for solving equations of electrodynamics and inverse problems of geophysics are required. In the present paper, the foundations of the theory of radio-interferometry sounding, the methodology, its role in geoelectric prospecting, and the opportunities for the solution of geological problems are discussed.

Comparative appraisal of linear inverse models constructed via distinctive parameterizations (comparing distinctly inverted models)

Most geophysical inverse problems deal with models of the continuous Earth structure. The classical Backus‐Gilbert's theory demonstrates that the resolvable model variation is the truth viewed through the resolution kernel that is woven by the data kernels. The actual numerical resolution amounts to the inversion of the Gram matrix formed by the inner products among data kernels. However, due to the usually sizable amount of data constraints and/or imperfection of the forward theory, the practical implementations are usually tackled through certain a priori finite parameterizations based on rather arbitrary choices of bases such as spatial voxels, splines, spherical harmonics, or spherical wavelets. The cross assessment on the consistency among inverse models parameterized or regularized differently has long been downplayed. It is shown in this study that straightforward conversions among different model representations also enable the direct conversions of the Gram matrices. This leads to significant flexibility in formulating the forward data rule in one representation and then carrying out the actual inversion in an alternate domain. Furthermore, it is also fairly easy to convert both the model covariance and resolution matrices across different representations. These conversions thus enable direct assessments across inverse models obtained via different parameterizations and different regularization schemes. An example utilizing preliminary results of an experiment of ambient noise tomography of a plate boundary region of complex tectonics for the northeast coast and offshore area of Taiwan is shown to demonstrate such comparisons.

Study of multiscale global optimization based on parameter space partition

Inverse problems in geophysics are usually described as data misfit minimization problems, which are difficult to solve because of various mathematical features, such as multi-parameters, nonlinearity and ill-posedness. Local optimization based on function gradient can not guarantee to find out globally optimal solutions, unless a starting point is sufficiently close to the solution. Some global optimization methods based on stochastic searching mechanisms converge in the limit to a globally optimal solution with probability 1. However, finding the global optimum of a complex function is still a great challenge and practically impossible for some problems so far. This work develops a multiscale deterministic global optimization method which divides definition space into sub-domains. Each of these sub-domains contains the same local optimal solution. Local optimization methods and attraction field searching algorithms are combined to determine the attraction basin near the local solution at different function smoothness scales. With Multiscale Parameter Space Partition method, all attraction fields are to be determined after finite steps of parameter space partition, which can prevent redundant searching near the known local solutions. Numerical examples demonstrate the efficiency, global searching ability and stability of this method.