Collection of articles inspired by the Mathematics of Planet Earth 2013 programme

In order to support the Mathematics of Planet Earth 2013 programme mpe2013.org, Inverse Problems is making all its related content published since 2010 available free until the end of 2013. These articles report significant research developments related to the themes of the programme and will be a useful resource to researchers in the mathematical, physical and biological sciences.

The authors of some of the articles have provided us with short 'Insight' summaries of the advances described in their articles, explaining how their research can help in understanding a wide range of phenomena that take place on our planet and in dealing with the ever-evolving challenges that life on the planet faces. We hope that you enjoy reading these Insights and that you find this collection of articles useful and stimulating.

Papers with Insights

An inverse radiative transfer model of the vegetation canopy based on automatic differentiation M Vossbeck, M Clerici, T Kaminski, T Lavergne, B Pinty and R Giering

Insight

All papers Show article list


Blind deconvolution of seismograms regularized via minimum support

A A Royer et al 2012 Inverse Problems 28 125010

The separation of earthquake source signature and propagation effects (the Earth's 'Green's function') that encode a seismogram is a challenging problem in seismology. The task of separating these two effects is called blind deconvolution. By considering seismograms of multiple earthquakes from similar locations recorded at a given station and that therefore share the same Green's function, we may write a linear relation in the time domain ui(t)*sj(t) − uj(t)*si(t) = 0, where ui(t) is the seismogram for the ith source and sj(t) is the jth unknown source. The symbol * represents the convolution operator. From two or more seismograms, we obtain a homogeneous linear system where the unknowns are the sources. This system is subject to a scaling constraint to deliver a non-trivial solution. Since source durations are not known a priori and must be determined, we augment our system by introducing the source durations as unknowns and we solve the combined system (sources and source durations) using separation of variables. Our solution is derived using direct linear inversion to recover the sources and Newton's method to recover source durations. This method is tested using two sets of synthetic seismograms created by convolution of (i) random Gaussian source-time functions and (ii) band-limited sources with a simplified Green's function and signal to noise levels up to 10% with encouraging results.

A wave-equation-based Kirchhoff operator

Fons ten Kroode 2012 Inverse Problems 28 115013

In this paper, I will study a Kirchhoff-type integral, which can be seen as a linear operator mapping angle–azimuth-dependent reflection coefficients along a reflector into reflection data for the acoustic wave equation. I will show that a minor adaptation of a construction of angle–azimuth-dependent images as proposed by Sava and Fomel leads to a left inverse of this operator, which maps primary reflection data to angle–azimuth-dependent reflection coefficients. The new construction naturally leads to a reformulation of the Kirchhoff operator, acting on space-shift-extended images, which can be implemented completely in terms of the fundamental solutions of the wave equation. I will study the composition of this new wave-equation-based Kirchhoff operator with an operator forming space-shift-extended images from data. I will show that these operators are partial inverses of each other, with their compositions being pseudo-differential operators that reconstruct suitably microlocalized versions of primary reflection data and extended images focused at space-shift zero.

Data inversion in coupled subsurface flow and geomechanics models

Marco A Iglesias and Dennis McLaughlin 2012 Inverse Problems 28 115009

We present an inverse modeling approach to estimate petrophysical and elastic properties of the subsurface. The aim is to use the fully coupled geomechanics-flow model of Girault et al (2011 Math. Models Methods Appl. Sci. 21 169–213) to jointly invert surface deformation and pressure data from wells. We use a functional-analytic framework to construct a forward operator (parameter-to-output map) that arises from the geomechanics-flow model of Girault et al. Then, we follow a deterministic approach to pose the inverse problem of finding parameter estimates from measurements of the output of the forward operator. We prove that this inverse problem is ill-posed in the sense of stability. The inverse problem is then regularized with the implementation of the Newton-conjugate gradient (CG) algorithm of Hanke (1997 Numer. Funct. Anal. Optim. 18 18–971). For a consistent application of the Newton-CG scheme, we establish the differentiability of the forward map and characterize the adjoint of its linearization. We provide assumptions under which the theory of Hanke ensures convergence and regularizing properties of the Newton-CG scheme. These properties are verified in our numerical experiments. In addition, our synthetic experiments display the capabilities of the proposed inverse approach to estimate parameters of the subsurface by means of data inversion. In particular, the added value of measurements of surface deformation in the estimation of absolute permeability is quantified with respect to the standard history matching approach of inverting production data with flow models. The proposed methodology can be potentially used to invert satellite geodetic data (e.g. InSAR and GPS) in combination with production data for optimal monitoring and characterization of the subsurface.

An inverse diffusivity problem for the helium production–diffusion equation

Gang Bao and Xiang Xu 2012 Inverse Problems 28 085002

Thermochronology is a technique for the extraction of information about the thermal history of rocks. Such information is crucial for determining the depth below the surface at which rocks were located at a given time (Bao G et al 2011 Commun. Comput. Phys. 9 129). Mathematically, extracting the time–temperature path can be formulated as an inverse diffusivity problem for the helium production–diffusion equation which is the underlying process of thermochronology. In this paper, to reconstruct the diffusivity which depends on space only and accounts for the structural information of rocks, a local Hölder conditional stability is obtained by a Carleman estimate. A uniqueness result is also proven for extracting the thermal history, i.e. identifying the time-dependant part of the diffusion coefficient, provided that it is analytical with respect to time. Numerical examples are presented to illustrate the validity and effectiveness of the proposed regularization scheme.

Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations

M Andrle and A El Badia 2012 Inverse Problems 28 075009

We consider the inverse problem of identifying multiple moving pollution sources in a linear advection–dispersion–reaction equation. Although we consider the specific application of pollution source identification in surface waters or atmospheric media, this problem has many other important applications in ecological and diffusive systems. We establish an identifiability result using observations on a non-empty subset of the domain boundary and develop an identification method by reformulating the inverse problem into a minimization problem. Finally, we provide numerical results to support the theoretical results.

Sparse regularization of inverse gravimetry—case study: spatial and temporal mass variations in South America

D Fischer and V Michel 2012 Inverse Problems 28 065012

Sparse regularization has recently experienced high popularity in the inverse problems community. In this paper, we show that a sparse regularization technique can also be developed for linear geophysical tomography problems. For this purpose, we adapt a known matching pursuit algorithm. The main theoretical features (existence, stability, and convergence) of the new method are given. We also show further properties of some trial functions which we use. Moreover, the algorithm is applied to a static and a monthly varying gravitational field of South America which yields spatial and temporal variations in the mass distribution. The new approach represents essential progress in comparison to a corresponding wavelet method, which is not flexible enough for the use of heterogeneous data, and a respective spline method, where the resolution cannot exceed approximately 104 basis functions due to experienced numerical problems with the ill-conditioned and dense matrix. The novel sparse regularization technique does not require homogeneous data and is not limited in the number of basis functions due to its iterative algorithm.

On the inverse problems for the coupled continuum pipe flow model for flows in karst aquifers

Shuai Lu et al 2012 Inverse Problems 28 065003

We investigate two inverse problems for the coupled continuum pipe flow (CCPF) model which describes fluid flows in karst aquifers. After generalizing the well posedness of the forward problem to the anisotropic exchange rate case which is a space-dependent variable, we present the uniqueness of this parameter by measuring the Cauchy data. The uniqueness of the geometry of the conduit by the Cauchy data is verified as well. These results enhance the practicality of the CCPF model.

Joint inversion approaches for geophysical electromagnetic and elastic full-waveform data

A Abubakar et al 2012 Inverse Problems 28 055016

We present joint inversion approaches for integrating controlled source electromagnetic data and seismic full-waveform data for geophysical applications. The first approach is the joint petrophysical inversion carried out by reconstructing petrophysical parameters such as porosity and saturations instead of the usual geophysical parameters such as resistivity, seismic velocities and mass density. This approach utilizes the strong correlation between the electromagnetic and seismic parameters through the petrophysical relationships. Another approach that does not require the a priori petrophysical correlation is the joint structural inversion method. In this approach, the inversion is carried out by employing a regularization function for enforcing the structural similarity between the resistivity and the seismic velocities and the mass density. In this method, we employ the cross-gradient function, which has been shown on many occasions to be quite effective. By using a time-lapse reservoir monitoring example, we show that both joint inversion approaches produce results that are superior to those obtained by disjointed inversions. Hence, these joint inversions have great potential to be the next-generation tools for reservoir characterization and monitoring.

Estimation of position and intensity of a pollutant source in channel flow using transmittance functions

T Maalej et al 2012 Inverse Problems 28 055010

This paper considers the one-dimension linear advection diffusion equation with constant coefficients for a transient pollutant point source. Theoretical transmittances deduced from the analytical solution of the transport equation (state-space), that is, transfer functions between concentrations, have been calculated. The convolution product between transmittances obtained by Laplace inversion and a reference concentration for channel flow application allows us to estimate the parameters of the model (velocity and diffusion coefficient). A nonlinear fixed point method based on least-squares estimation based on noisy concentration responses and using truncated singular value decomposition regularization is implemented to determinate gaseous contaminant source locations and strengths.

A boundary perturbation method for recovering interface shapes in layered media

Alison Malcolm and David P Nicholls 2011 Inverse Problems 27 095009

The scattering of linear acoustic radiation by a periodic layered structure is a fundamental model in the geosciences as it closely approximates the propagation of pressure waves in the earth's crust. In this contribution, the authors describe new algorithms for (1) the forward problem of prescribing incident radiation and, given the known structure, determining the scattered field, and (2) the inverse problem of approximating the form of the structure given prescribed incident radiation and measured scattered data. Each of these algorithms is based upon a novel statement of the problem in terms of boundary integral operators (Dirichlet–Neumann operators), and a boundary perturbation algorithm (the method of operator expansions) for their evaluation. Detailed formulas and numerical simulations are presented to demonstrate the utility of these new approaches.

Bayesian matched-field geoacoustic inversion

Stan E Dosso and Jan Dettmer 2011 Inverse Problems 27 055009

This paper describes a Bayesian approach to matched-field inversion (MFI) of ocean acoustic data for seabed geoacoustic properties. In a Bayesian formulation, the unknown environmental and experimental parameters are considered random variables constrained by noisy data and prior information, and the goal is to interpret the multi-dimensional posterior probability density (PPD). The PPD is typically characterized in terms of point estimates, marginal distributions, and posterior correlations (or joint statistics). Computing these requires numerical optimization and integration of the PPD, which are carried out efficiently here using adaptive hybrid optimization and Metropolis–Hastings sampling in principal-component space, respectively. Likelihood and misfit functions for multi-frequency MFI with incomplete source spectral information are derived based on the assumption of complex Gaussian-distributed data errors with covariance matrices estimated from residual analysis; posterior statistical tests are applied to validate these estimates and assumptions. Model selection is carried out by applying the Bayesian information criterion to determine the simplest seabed parameterization consistent with the resolving power of the data. Bayesian MFI is illustrated for shallow-water acoustic data measured in the Mediterranean Sea.

Inverse bifurcation problems for nonlinear Sturm–Liouville problems

Tetsutaro Shibata 2011 Inverse Problems 27 055003

We consider the inverse nonlinear eigenvalue problems for the equation which is motivated biologically by the problem of population dynamics. It is assumed that f(u) is an unknown nonlinear term. Under the standard growth conditions on f, for any given α > 0, there exists a unique solution of the equation with ||uα||q = α, where || ⋅ ||q denotes usual Lq-norm. This curve λ(q, α) is called the Lq-bifurcation curve. We show that (i) the unknown nonlinear term is determined as f(u) = up (p > 1) nearly exponentially for u ≫ 1 under the reasonable condition on λ(q, α) for α ≫ 1, and (ii) by variational method, the unknown nonlinear term is determined uniquely for u ⩾ 0 under the additional conditions on f when q = 2.

A 2D nonlinear inversion of well-seismic data

Ludovic Métivier et al 2011 Inverse Problems 27 055005

Well-seismic data such as vertical seismic profiles are supposed to provide detailed information about the elastic properties of the subsurface at the vicinity of the well. Heterogeneity of sedimentary terrains can lead to far from negligible multiple scattering, one of the manifestations of the nonlinearity involved in the mapping between elastic parameters and seismic data. We present a 2D extension of an existing 1D nonlinear inversion technique in the context of acoustic wave propagation. In the case of a subsurface with gentle lateral variations, we propose a regularization technique which aims at ensuring the stability of the inversion in a context where the recorded seismic waves provide a very poor illumination of the subsurface. We deal with a huge size inverse problem. Special care has been taken for its numerical solution, regarding both the choice of the algorithms and the implementation on a cluster-based supercomputer. Our tests on synthetic data show the effectiveness of our regularization. They also show that our efforts in accounting for the nonlinearities are rewarded by an exceptional seismic resolution at distances of about 100 m from the well. They also show that the result is not very sensitive to errors in the estimation of the velocity distribution, as far as these errors remain realistic in the context of a medium with gentle lateral variations.

Level-set techniques for facies identification in reservoir modeling

Marco A Iglesias and Dennis McLaughlin 2011 Inverse Problems 27 035008

In this paper we investigate the application of level-set techniques for facies identification in reservoir models. The identification of facies is a geometrical inverse ill-posed problem that we formulate in terms of shape optimization. The goal is to find a region (a geologic facies) that minimizes the misfit between predicted and measured data from an oil–water reservoir. In order to address the shape optimization problem, we present a novel application of the level-set iterative framework developed by Burger in (2002 Interfaces Free Bound. 5 301–29; 2004 Inverse Problems 20 259–82) for inverse obstacle problems. The optimization is constrained by (the reservoir model) a nonlinear large-scale system of PDEs that describes the reservoir dynamics. We reformulate this reservoir model in a weak (integral) form whose shape derivative can be formally computed from standard results of shape calculus. At each iteration of the scheme, the current estimate of the shape derivative is utilized to define a velocity in the level-set equation. The proper selection of this velocity ensures that the new shape decreases the cost functional. We present results of facies identification where the velocity is computed with the gradient-based (GB) approach of Burger (2002) and the Levenberg–Marquardt (LM) technique of Burger (2004). While an adjoint formulation allows the straightforward application of the GB approach, the LM technique requires the computation of the large-scale Karush–Kuhn–Tucker system that arises at each iteration of the scheme. We efficiently solve this system by means of the representer method. We present some synthetic experiments to show and compare the capabilities and limitations of the proposed implementations of level-set techniques for the identification of geologic facies.

Identification of moving pointwise sources in an advection–dispersion–reaction equation

M Andrle et al 2011 Inverse Problems 27 025007

We consider the inverse problem of identifying a moving source in a linear advection–dispersion–reaction equation. The main application, but not the only one, is the identification of an environmental pollution source in a river. An identifiability result is established and an identification method proposed using measurement records at two locations, one upstream and the other downstream from the source. Finally, numerical simulations are performed to assess the identification process.

Reconstruction of river bed topography from free surface data using a direct numerical approach in one-dimensional shallow water flow

A F Gessese et al 2011 Inverse Problems 27 025001

River bed topography is of paramount importance for the study of fluvial hydraulics, flood prediction and river flow monitoring. It is therefore important to develop fast, easy-to-implement and cost-effective methods to determine underwater river topography. This paper presents a new one-step approach for reconstructing the river bed topography from known free surface data. This problem corresponds to the inverse of the classical hydrodynamic problem where the shallow water equations provide the free surface profile for a given river bed. We show in this work that instead of treating this inverse problem in the traditional partial differential equation (PDE)-constrained optimization framework, we can conveniently rearrange the governing equations for the direct problem to obtain an explicit PDE for the inverse problem. This leads to a direct solution of the inverse problem. An interesting consequence of the analysis is that the equations governing the forward and inverse problems have a very similar form, and the same discretization technique, based on an upwind conservative numerical scheme, can be used. The proposed methodology is successfully tested on a range of benchmark problems for noisy and noiseless free surface data. It was found that this solution approach creates very little amplification of noise.

Illumination analysis of wave-equation imaging with curvelets

Shen Wang et al 2010 Inverse Problems 26 115013

We present a comprehensive framework for wave-equation illumination analysis and introduce a target-oriented illumination correction that simultaneously accounts for a limited acquisition aperture and, locally, compensates for the so-called normal operator in inverse scattering to yield a partial reconstruction of reflectivity or reflection coefficient, while minimizing (orientation-dependent) phase distortions and artifacts. The corrections are nested while the most significant correction is typically the one associated with the limited acquisition aperture. To carry out the analysis, we make use of higher dimensional curvelets, which provide the means of extracting directional information, and introduce associated matrix representations for the component operators, including the tapers associated with the acquisition aperture, that make up wave-equation imaging; we essentially exploit the properties of these matrices. Curvelets can be viewed as 'fat', localized plane waves and hence form a natural candidate to generalize geophysical diffraction tomography, which is at the basis of our approach. Our approach admits the formation of caustics and, hence, is valid in complex wave-speed models, though these need to be known for the illumination compensation to be effective.

Effective solution of nonlinear subsurface flow inverse problems in sparse bases

Lianlin Li and Behnam Jafarpour 2010 Inverse Problems 26 105016

Identification of spatially variable hydraulic rock properties such as permeability and porosity is necessary for accurate prediction of fluid flow displacement in subsurface environments. The estimation of these properties from dynamic flow data usually involves solving a highly underdetermined nonlinear inverse problem in which a limited set of measurements is combined with prior knowledge to estimate the unknown parameters. The overwhelming number of unknowns, relative to available data, leads to many parameter combinations that explain the data equally well, but fail to predict the flow behavior in the reservoir. To improve solution non-uniqueness and numerical stability, additional information is typically incorporated into the solution procedure. One way to regularize underdetermined inverse problems is to demand certain structural properties from the solution that are usually derived from the physics of the problem. In this paper, we exploit the compact representation of spatially correlated geologic formations in sparse bases to facilitate the reconstruction of rock hydraulic property distributions from flow measurements that are nonlinearly related to unknown parameters. By formulating the solution in a compressive basis such as the wavelet or Fourier, we show that minimizing a data misfit cost function augmented with an additive or multiplicative regularization term that promotes sparse solutions, the reconstruction results can be improved. Convergence to a relevant sparse solution is adaptively carried out through an iteratively reweighting algorithm in the transform domain. We evaluate the performance of our inversion algorithm using a set of two-phase waterflooding experiments in an oil reservoir where nonlinear dynamic flow data are integrated to infer the spatial distribution of rock permeability.

An inverse radiative transfer model of the vegetation canopy based on automatic differentiation

M Voßbeck et al 2010 Inverse Problems 26 095003

This paper presents an inverse model of radiation transfer processes occurring in the solar domain in vegetation plant canopies. It uses a gradient method to minimize the misfit between model simulation and observed radiant fluxes plus the deviation from prior information on the unknown model parameters. The second derivative of the misfit approximates uncertainty ranges for the estimated model parameters. In a second step, uncertainties are propagated from parameters to simulated radiant fluxes via the model's first derivative. All derivative information is provided by a highly efficient code generated via automatic differentiation of the radiative transfer code. The paper further derives and evaluates an approach for avoiding secondary minima of the misfit. The approach exploits the smooth dependence of the solution on the observations, and relies on a database of solutions for a discretized version of the observation space.

Amplitude calculations for 3D Gaussian beam migration using complex-valued traveltimes

Norman Bleistein and Samuel H Gray 2010 Inverse Problems 26 085017

Gaussian beams are often used to represent Green's functions in three-dimensional Kirchhoff-type true-amplitude migrations because such migrations made using Gaussian beams yield superior images to similar migrations using classical ray-theoretic Green's functions. Typically, the integrand of a migration formula consists of two Green's functions—each describing propagation to the image point—one from the source and the other from the receiver position. The use of Gaussian beams to represent each of these Green's functions in 3D introduces two additional double integrals when compared to a Kirchhoff migration using ray-theoretic Green's functions, thereby adding a significant computational burden. Hill (2001 Geophysics 66 1240–50) proposed a method for reducing those four integrals to two, compromising slightly on the full potential quality of the Gaussian beam representations for the sake of more efficient computation. That approach requires a two-dimensional steepest descent analysis for the asymptotic evaluation of a double integral. The method requires evaluation of the complex traveltimes of the Gaussian beams as well as the amplitudes of the integrands at the determined saddle points. In addition, it is necessary to evaluate the determinant of a certain (Hessian) matrix of second derivatives. Hill (2001 Geophysics 66 1240–50) did not report on this last part; thus, his proposed migration formula is kinematically correct but lacks correct amplitude behavior. In this paper, we derive a formula for that Hessian matrix in terms of dynamic ray tracing quantities. We also show in a simple example how the integral that we analyze here arises in a true amplitude migration formula.

Application of a two-and-a-half dimensional model-based algorithm to crosswell electromagnetic data inversion

Maokun Li et al 2010 Inverse Problems 26 074013

In this paper, we apply a model-based inversion scheme for the interpretation of the crosswell electromagnetic data. In this approach, we use open and closed polygons to parameterize the unknown configuration. The parameters that define these polygons are then inverted for by minimizing the data misfit cost function. Compared with the pixel-based inversion approach, the model-based inversion uses only a few number of parameters; hence, it is more efficient. Furthermore, with sufficient sensitivity in the data, the model-based approach can provide quantitative estimates of the inverted parameters such as the conductivity. The model-based inversion also provides a convenient way to incorporate a priori information from other independent measurements such as seismic, gravity and well logs.

Guided Bayesian optimal experimental design

M R Khodja et al 2010 Inverse Problems 26 055008

A Bayesian methodology is described for designing experiments or surveys that will optimally complement all previously available information. This methodology uses strong prior information to linearize the problem, and to guide the design toward maximally reducing forecast uncertainties in the interpretation of the future experiment. The prior information could possibly be correlated among model parameters or the observation noise. With no prior information this approach reduces to the fast recursive implementation of the D-optimality criterion. Synthetic geophysical tomography examples are used to illustrate the benefits of this approach. The dependence of the design on the model representation is also discussed.

A comparison of seismic velocity inversion methods for layered acoustics

T van Leeuwen and W A Mulder 2010 Inverse Problems 26 015008

In seismic imaging, one tries to infer the medium properties of the subsurface from seismic reflection data. These data are the result of an active source experiment, where an explosive source and an array of receivers are placed at the surface. Because of the absence of low frequencies in the data, the corresponding inverse problem is strongly nonlinear in the slowly varying component of the velocity. The least-squares misfit functional typically exhibits local minima and has a small basin of attraction. The usual approach to fitting the data in a least-squares sense by employing a gradient-based optimization method will therefore most likely result in a wrong velocity model. In the geophysical community, this problem has long been recognized and alternative formulations of the inverse problem have been developed. We review several of these formulations and analyse the sensitivity to the error in the smooth velocity component. This analysis is carried out for laterally homogeneous velocities using an asymptotic solution of the wave equation. The analysis suggests that formulations which are geared towards fitting the phases of the data, rather than the amplitudes, have smooth corresponding misfit functionals with a large basin of attraction.