Linear Inverse Theory Research Articles

Bayesian linearized inversion for petrophysical and pore-connectivity parameters with seismic elastic data of carbonate reservoirs

Abstract Carbonate reservoirs are important targets for promoting the oil and gas reserve exploration and production in China. However, such reservoirs usually contain developed complex pore structures, which heavily affect the precision in seismic prediction of petrophysical parameters. As one of the most important parameters to characterize reservoir rock, pore-related parameters can not only describe the pore structure, but also be used to evaluate the oil/gas-bearing capabilities of potential reservoirs. The conventional rock-physics models (e.g. Gassmann's model) are formulated assuming fully connected pores, which is unable to accurately capture the geometrical complexity in real rocks. To characterize the influences of multiple pores on the elastic properties, this work presents a rock-physics modeling method for carbonates, wherein the percentage composition of connected pores is equivalently quantified as the pore-connectivity factor. The method treats the pore-connectivity factor as an objective variable to characterize the spatial variations of pore structure. Specifically, the method combines the differential equivalent medium theory and Gassmann's model, and derives a linearized forward operator to quantitatively link porosity, water saturation, and pore-connectivity factor to seismic elastic parameters. According to the Bayesian linear inverse theory, the simultaneous estimation of petrophysical and pore-connectivity parameters is achieved. To characterize the statistical variations with multiple lithofacies, the Gaussian mixture model is employed to quantify the prior distribution of the objective variables. The posterior distribution of the objective variables is analytically expressed with the linearized forward operator. Numerical experiments show that the accuracy of the proposed method in predicting elastic parameters is improved. Compared with the conventional Xu–White model and the varying pore aspect-ratio method, the accuracy of predicted P-wave velocity increases by 10.29% and 1.33%, respectively, and the predicted S-wave velocity increases by 6.44% and 0.03%, in terms of correlation coefficient. The application to the field data validates the effectiveness of the method, wherein the porosity and water saturation results help indicating the spatial distribution of potential reservoirs.

Waveform Imaging Based on Linear Forward Representations for Scalar Wave Seismic Data

The current reverse-time migration, which is based on wave equations for imaging wavefields, employs an imaging formula derived from Claerbout’s imaging principle. This imaging formula is only valid for plane waves with small incident angles on the perfectly flat reflecting surface. However, the complexity of seismic wave propagation may lead to situations that do not meet this requirement. Therefore, this paper divides the subsurface into local scattering and reflecting bodies. It proposes linear forward representations for scattering and reflection data based on perturbations in the physical parameters and wave impedance, respectively. To further describe the effect on the reflecting body boundary, the local reflection coefficient is defined and the linear forward representation for the reflection data based on it is obtained. After that, the proposed linear forward representations are used as the forward equations for the linear inverse of the seismic data, and the seismic data waveform imaging method is developed based on linear inversion theory. At the same time, the specific waveform imaging calculation formulas for scalar wave scattering data and scalar wave reflection data are provided and validated via numerical experiments. Compared with the current reverse-time migration, waveform migration not only has the correct phase and higher resolution in theory but also does not increase the computational complexity. To some extent, it improves the deficiencies of the current structural imaging and provides a basis for subsurface lithological imaging.

Stochastic inversion of geophysical data by a conditional variational autoencoder

Recovering geologically realistic physical property models by geophysical inversion is a long-standing challenge. Generative neural networks offer a promising path to meet this challenge because they can produce spatially complex models that exhibit the characteristics of a set of training models, even when those characteristics are difficult to quantify. In the context of geophysical inversion, these characteristics may include faults, layers, and sharp contacts between rock units. Here, we develop a framework for incorporating prior geologic knowledge into geophysical inversions using conditional variational autoencoders (CVAEs). We train a CVAE to reconstruct training density models while honoring relative gravity data. Once trained, the decoder network of the CVAE inverts gravity data. The inputs to the decoder are observed gravity data and a set of latent variables that are sampled from a standard normal distribution. The decoder maps from the observed data and latent variables to density models such that the resulting models are consistent with the training models and the input data. Consequently, the inversion fits the observed data and incorporates the information embedded in the training models. The decoder can produce many inverted models instantaneously, sampling an approximation to the posterior model distribution efficiently. We find that the latent variables correspond to independent interpretable ways in which the model can vary while still honoring the observed data. We draw a connection to linear inverse theory, positing that the latent variables are analogous to the principal components of the local posterior model covariance.

High-resolution topography from planetary images and laser altimetry

Mapping landforms on the Moon is of great interest and importance for future human settlements and resources exploration. One of the first steps is to map the topography in great detail and resolution. However, data from the Lunar Orbiter Laser Altimeter (LOLA) provide low-resolution elevation maps in comparison to the size of detailed geological features. To improve resolution, we developed a new method to upscale topographic maps to a higher resolution using images from the Lunar Reconnaissance Orbiter Camera (LROC). Our method exploits the relation between topographic gradients and degrees of shading of incoming sunlight. In contrast to earlier published methods, our approach is based on probabilistic, linear inverse theory, and its computational efficiency is very high due to its formulation through the Sylvester Equation. The method operates on multiple images and incorporates albedo variations. A further advantage of the method is that we avoid/reduce the use of arbitrary tuning parameters through a probabilistic formulation where all weighting of data and model parameters is based on prior information about data uncertainties and reasonable bounds on the model. Our results increase the resolution of the topography from ∼60 ​m per pixel to 0.9 ​m per pixel, bringing it to the same pixel resolution as the optical images from LROC, allowing in some cases detection of craters as small as ∼3 ​m of diameter. We estimate uncertainties of the topographic model due to noise in the images, and in the low-resolution (LOLA) model. • Mapping planetary landforms in detailed resolution is hampered by the poor resolution of Laser Altimeter instruments. • Method computes high-resolution topography from the relation between topographic slopes and degrees of shading of sunlight. • The method uses probabilistic inverse theory. It is computationally efficient and avoids arbitrary tuning parameters. • We obtain DEMs of 0.9 m/pixel resolution (same as used optical images), showing surface features as small as ~3m diameter.

A mixed, unified forward/inverse framework for earthquake problems: fault implementation and coseismic slip estimate

SUMMARY We introduce a new finite-element (FE) based computational framework to solve forward and inverse elastic deformation problems for earthquake faulting via the adjoint method. Based on two advanced computational libraries, FEniCS and hIPPYlib for the forward and inverse problems, respectively, this framework is flexible, transparent and easily extensible. We represent a fault discontinuity through a mixed FE elasticity formulation, which approximates the stress with higher order accuracy and exposes the prescribed slip explicitly in the variational form without using conventional split node and decomposition discrete approaches. This also allows the first order optimality condition, that is the vanishing of the gradient, to be expressed in continuous form, which leads to consistent discretizations of all field variables, including the slip. We show comparisons with the standard, pure displacement formulation and a model containing an in-plane mode II crack, whose slip is prescribed via the split node technique. We demonstrate the potential of this new computational framework by performing a linear coseismic slip inversion through adjoint-based optimization methods, without requiring computation of elastic Green’s functions. Specifically, we consider a penalized least squares formulation, which in a Bayesian setting—under the assumption of Gaussian noise and prior—reflects the negative log of the posterior distribution. The comparison of the inversion results with a standard, linear inverse theory approach based on Okada’s solutions shows analogous results. Preliminary uncertainties are estimated via eigenvalue analysis of the Hessian of the penalized least squares objective function. Our implementation is fully open-source and Jupyter notebooks to reproduce our results are provided. The extension to a fully Bayesian framework for detailed uncertainty quantification and non-linear inversions, including for heterogeneous media earthquake problems, will be analysed in a forthcoming paper.

Inversion-based imaging: from LSRTM to FWI imaging

Summary Least Squares RTM (LSRTM) is a powerful inversion-based imaging algorithm which minimizes the data misfit between observed seismic recordings and forward-modelled synthetic data. The algorithm, which can be implemented in either data or image domains, carries a fundamental limitation because it is based on a linear inversion theory which cannot accommodate velocity refinement as part of its model update process. Successful application of LSRTM therefore requires highly accurate velocity information, and if the velocity model is in significant error, modelled events will not be aligned kinematically with the observed data, and the algorithm will tend to produce unsatisfactory results. FWI is another inversion-based algorithm that enjoys widespread industry use. Unlike LSRTM, FWI poses its inverse problem within a non-linear framework whereby it updates the velocity model and associated wave paths throughout its iterative process, gradually aligning modelled events with observed events. With the recent convergence of FWI and LSRTM methodologies, FWI is not only being used as a velocity update tool, but also as a direct imaging tool, thereby achieving two key imaging goals, namely refining the velocity model and deriving a better-quality seismic image. The latter process, which is known as ‘FWI imaging’, has recently been gaining a lot of industry attention as it offers the possibility of high-quality imaging along with workflow simplification. In this article we will compare and contrast LSRTM and FWI. We conclude that the process of generating the FWI-imaging essentially amounts to nonlinear, data-domain inversion. This recognition facilitates a ready comparison against the data-domain form of LSRTM, the latter being a linear, data-domain inversion.

Angle‐domain imaging of relative impedance perturbation

ABSTRACTImpedance is a physical parameter that plays an important role in seismic data processing and interpretation. A relative impedance perturbation (the ratio of the impedance perturbation and the impedance for the background models) imaging method in depth domain based on the reflection wave equation is proposed. Under the small perturbation assumption, primary wave and high‐frequency approximation condition, a linear propagation equation of the primary reflection waves based on the relative impedance perturbation was first derived. On this basis, we further derived the imaging formula of the relative impedance perturbation using a linear inversion theory. Then, the source–receiver bidirectional illumination compensation was used to improve the image quality of the subsurface structures. The image result obtained by this method can be used to estimate the relative impedance perturbation. In the angle domain, the extracted near‐angle‐domain image gather with amplitude compensation can estimate the relative impedance perturbation, and the far‐angle image gather provides the estimation of the relative velocity perturbation (the ratio of the velocity perturbation and the background velocity). Finally, several numerical tests demonstrate the effectiveness of the method.

Prediction of Favorable Carbonate Reservoirs under Extremely Thick Salts via Poststack Facies-Controlled and Prestack Zoeppritz Equation Inversions in the Santos Basin of Brazil

Subsalt carbonate rocks in Brazil’s deepwater region have broad prospects for oil and gas exploration and development. Due to the low-frequency bandwidth of the seismic data and the poor signal quality for this kind of reservoir target, there is a demand for accurate seismic prediction methods. We employ the facies-controlled inversion using a low-pass filter matrix to ensure the accuracy of the low frequency and to improve the robustness of the inversion results. We integrated the concept of adaptive regularization constraint of the Zoeppritz equation into the generalized linear inversion theory framework, which overcomes the shortcomings of the approximate equation. Making full use of the large angle prestack seismic information, Zoeppritz equation inversion improves the accuracy of the inversion results. The application of this method in carbonate reservoirs under extremely thick salts in the Santos Basin of Brazil indicates the feasibility and practicality of the proposed integrated prediction method.

Gaussian process models—I. A framework for probabilistic continuous inverse theory

SUMMARYWe develop a theoretical framework for framing and solving probabilistic linear(ized) inverse problems in function spaces. This is built on the statistical theory of Gaussian Processes, and allows results to be obtained independent of any basis, avoiding any difficulties associated with the fidelity of representation that can be achieved. We show that the results of Backus–Gilbert theory can be fully understood within our framework, although there is not an exact equivalence due to fundamental differences of philosophy between the two approaches. Nevertheless, our work can be seen to unify several strands of linear inverse theory, and connects it to a large body of work in machine learning. We illustrate the application of our theory using a simple example, involving determination of Earth’s radial density structure.

Reflected wave least squares reverse time migration with angle illumination compensation

High-quality seismic data imaging plays an important role in the lithological interpretation of subsurface structures. However, high-quality imaging remains a challenging task. Based on the linear inversion theory of reflected wave equations, this paper proposes reflected wave least squares reverse time migration with angle illumination compensation to better balance the amplitude of seismic imaging. We use the reflected wave migration equation to unify forward and backward propagation, which helps to obtain an image with correct phase and symmetric waveform. Under the assumption that the spectrum of seismic wavefield remains unchanged, the Poynting vector method is used to efficiently calculate the propagation direction of seismic waveform and seismic illumination in the angle domain. During iteration, angle-domain illumination is used as a preconditioner to compensate for the amplitude of the iterated gradient terms based on the angle value. In this manner, we can enhance the imaging energy of steeply inclined structures. To improve the stability of linear inversion, the spatial derivative of the image is used as a regularized constraint term. Numerical tests show that the proposed method can suppress imaging noise as well as improve resolution and amplitude fidelity of the images. Furthermore, the inversed result can be used to estimate underground reflectivity, which is important for the further development of seismic inversion technology.

The diminishing residual matrices method for surface-consistent statics — A review

The diminishing residual matrices (DRM) method can be used to surface-consistently decompose individual trace statics into source and receiver components. The statics to be decomposed may either be first-arrival times after the application of linear moveout associated with a consistent refractor as used in refraction statics or residual statics obtained by crosscorrelating individual traces with corresponding model traces (known as pilot traces) at the same common-midpoint (CMP) location. The DRM method is an iterative process like the well-known Gauss-Seidel (GS) method, but it uses only source and receiver terms. The DRM method differs from the GS method in that half of the average common shot and receiver terms are subtracted simultaneously from the observations at each iteration. DRM makes the under-constrained statics problem a constrained one by implicitly adding a new constraint, the equality of the contribution of shots and receivers to the solution. The average of the shot statics and the average of the receiver statics are equal in the DRM solution. The solution has the smallest difference between shot and receiver statics profiles when the number of shots and the number of receivers in the data are equal. In this case, it is also the smallest norm solution. The DRM method can be derived from the well-known simultaneous iterative reconstruction technique. Simple numerical tests as well as results obtained with a synthetic data set containing only the field statics verify that the DRM solution is the same as the linear inverse theory solution. Both algorithms can solve for the long-wavelength component of the statics if the individual picks contain them. Yet DRM method is much faster. Application of the method to the normal moveout-corrected CMP gathers on a 3D land survey for residual statics calculation found that pick-decompose-apply-stack stages of the DRM method need to be iterated. These iterations are needed because of time and waveform distortions of the pilot traces due to the individual trace statics. The distortions lessen at every external DRM iteration.

Three-dimensional inversion technique in ocean and seabed acoustics using the parabolic equation method

A three-dimensional parabolic equation (PE) and perturbation approach is used to invert for the range-dependent geoacoustic characteristics of the seabed. The model assumes that the sound speed profile is the superposition of a known range-independent profile and an unknown depth and range-dependent perturbation. Using a Green’s function approach, the total measured pressure field in the water column is decomposed into a background field, which is due to the range-independent profile, and a scattered field, which is due to the range-dependent perturbation. When the Born approximation is applied to the resulting integral equation, it can be solved for the range-dependent profile using linear inverse theory. For simplicity, the sound speed profile in the water column is assumed to be known, and the range-dependent perturbation is added to the index of refraction n(x,y,z), rather than the sound speed profile c(x,y,z). The method is implemented in both Cartesian (x,y,z) and cylindrical (r,θ,z) coordinates with the forward propagation of the field in x and r, respectively. Synthetic and experimental data are used to demonstrate the validity of the method. Keywords: Three-dimensional parabolic equation method, geoacoustic inversion, range-dependent sound speed profile, linear inversion, Born approximation, and Green’s functions.

Zoeppritz equation-based prestack inversion and its application in fluid identification

Prestack seismic inversion methods adopt approximations of the Zoeppritz equations to describe the relation between reflection coefficients and P-wave velocity, S-wave velocity, and density. However, the error in these approximations increases with increasing angle of incidence and variation of the elastic parameters, which increases the number of inversion solutions and minimizes the inversion accuracy. In this study, we explore a method for solving the reflection coefficients by using the Zoeppritz equations. To increase the accuracy of prestack inversion, the simultaneous inversion of P-wave velocity, S-wave velocity, and density by using prestack large-angle seismic data is proposed based on generalized linear inversion theory. Moreover, we reduce the ill posedness and increase the convergence of prestack inversion by using the regularization constraint damping factor and the conjugate gradient algorithm. The proposed prestack inversion method uses prestack large-angle seismic data to obtain accurate seismic elastic parameters that conform to prestack seismic data and are consistent with logging data from wells.

Evaluating the utility of gravity gradient tensor components

Gravity gradiometry allows individual components and combinations of components to be used in interpretation. Knowledge of the information content of different components and their combinations is therefore crucial to their effectiveness, and a quantitative rating of information level is needed to guide the choice. To this end, I use linear inverse theory to examine the relationship between the different tensor components and combinations thereof and the model parameters to be determined. The model used is a rectangular prism, characterized by seven parameters: the prism location [Formula: see text], [Formula: see text]; its width [Formula: see text] and breadth [Formula: see text]; the density [Formula: see text]; the depth to top [Formula: see text]; and thickness [Formula: see text]. Varying these values allows a variety of body shapes, e.g., blocks, plates, dykes, and rods, to be considered. The Jacobian matrix, which relates parameters and their associated gravity response, clarifies the importance and stability of model parameters in the presence of data errors. In general, for single tensor components and combinations, the progression from well to poorly determined parameters follows the trend of [Formula: see text], [Formula: see text], [Formula: see text], [Formula: see text], [Formula: see text], [Formula: see text], to [Formula: see text]. Ranking the estimated model errors from a range of models showed that data sets consisting of concatenated components produced the smallest parameter errors. For data sets comprising combined tensor components, the invariants of the tensor produced the smallest parameter errors. Of the single tensor components, [Formula: see text] gave the best performance overall, but those single components with strong directional sensitivity can produce some individual parameters with smaller estimated errors (e.g., [Formula: see text] and [Formula: see text] estimated from [Formula: see text]).

Evaluating the utility of gravity gradient tensor components

Gravity gradiometry offers multiple single components and possible combinations of components to be used in interpretation. Knowledge of the information content of components and their combinations is therefore crucial to their effectiveness and so a quantitative rating of information level is needed to guide the choice. To this end we use linear inverse theory to examine the relationship between the different tensor components and combinations thereof and the model parameters to be determined. The model used is a simple prism, characterized by seven parameters: the prism location, xc, yc, its width w and breadth b, the density I?, the depth to top z, and thickness t. Varying these values allows a wide variety of body shapes, e.g., blocks, plates, dykes, rods, to be considered. The Jacobian matrix, which relates parameters and their associated gravity response, clarifies the importance and stability of model parameters in the presence of data errors. In general, for single tensor components and combinations, the progression from well- to poorly-determined parameters follows the trend of I?, xc, yc, w, b, z to t. Ranking the estimated model errors from a range of models shows that data sets consisting of concatenated components produce the smallest parameter errors. For data sets comprising combined tensor components, the invariants I1 and I2 produce the smallest model errors. Of the single tensor components, Tzz gives the best performance overall, but those single components with strong directional sensitivity can produce some individual parameters with smaller estimated errors (e.g., w and xc estimated from Txz).

Estimating the electron energy distribution during ionospheric modification from spectrographic airglow measurements

The electron energy distribution during an Fregion ionospheric modification experiment at the HAARP facility near Gakona, Alaska, is inferred from spectrographic airglow emission data. Emission lines at 630.0, 557.7, and 844.6 nm are considered along with the absence of detectable emissions at 427.8 nm. Estimating the electron energy distribution function from the airglow data is a problem in classical linear inverse theory. We describe an augmented version of the method of Backus and Gilbert which we use to invert the data. The method optimizes the model resolution, the precision of the mapping between the actual electron energy distribution and its estimate. Here, the method has also been augmented so as to limit the model prediction error. Model estimates of the suprathermal electron energy distribution versus energy and altitude are incorporated in the inverse problem formulation as representer functions. Our methodology indicates a heater‐induced electron energy distribution with a broad peak near 5 eV that decreases approximately exponentially by 30 dB between 5–50 eV.

Evaluation of Joint Inversion Algorithm with Logarithmic Treatment for Resistivity Log

The joint inversion of dual laterolog and dual-induction log was evaluated from the linear inversion theory. The generalized inversion method was used for the inverse calculation. A couple of parameters in the inversion algorithm were treated using logarithmic method. The abnormal states resulting from the logarithmic treatment on equations were analyzed. The effects of both the data-resolution matrix and model-resolution matrix on convergence of inversion process were discussed. It indicated that the logarithm treatment can not only weaken the abnormal state obviously, but also can reduce the spread of the data-resolution matrix. The improved inversion algorithm is efficient to enhance the stability and to improve the behaviors of convergence during the inversion process.

A Relaxed Matrix Inversion Method for Retrieving Water Constituent Concentrations in Case II Waters: The Case of Lake Kasumigaura, Japan

The matrix inversion method (MIM) is an effective algorithm for estimating water constituent concentrations in case II waters. To apply this method, appropriate and accurate specific inherent optical properties (SIOPs) for each constituent in water are essential. However, many routine observations of lake water quality do not in fact provide SIOPs, thus limiting the application of the MIM. In this paper, an alternative MIM method based on linear matrix inversion theory was proposed to relax the requirement of SIOPs measurement. For this, so-called ESIOPs (Estimated SIOPs) were first derived by an unusual application of MIM based on adequate calibration samples; then the water constituent concentrations for the whole study area were retrieved by the standard application of MIM based on the derived ESIOPs. For each calibration sample, measurement of the reflectance spectrum and corresponding water constituent concentrations, which can be obtained from periodical satellite data and routine field surveys, is required. The performance of the proposed method was evaluated using the simulation data from Hydrolight and three MEdium Resolution Imaging Spectrometer Instrument (MERIS) images. The results showed that this method yielded satisfactory estimations of the water constituent concentrations for the noise-contaminated simulation data sets. For MERIS data in our study area (Lake Kasumigaura, Japan), the average bias (mean normalized bias or MNB) and relative random uncertainty (normalized root mean square error, or NRMS) were in the range of -11.2% to 3.4% and 4.8% to 29.7% for each water constituent concentration. These findings imply that the algorithm proposed in this study is theoretically reasonable and practically applicable.

Seabed shear modulus profile inversion using surface gravity (water) wave-induced bottom motion

Summary. The motion of a layered seabed induced by propagating gravity water waves is modelled by the coupling matrix for a massless incompressible elastic bed according to Yamamoto. An amplitude inversion scheme is developed to extract the bottom shear modulus profile from the motion of the seabed at a point on the bed surface using the linear inverse theory combined with the Yamamoto theory of wave-seabed interaction. Numerical tests using synthetic data without noise confirmed that the inversion is unique and consistent. The inversion technique is applied to a set of field data from the Mississippi River Delta. A good agreement is obtained between the inverted profile and the direct measurements.

LINEAR INVERSE ANALYSIS OF STRESS AND SHAPE CONTROL OF STRUCTURES WITH MEMBERS OF VARIABLE LENGTH

There are many structures that are controlled their shapes and/or stress distribution by being actively varied the length of their members. Some of the structures with vibration control system, so-called variable geometry trusses and many of the tension structures are examples of such structures. One of the difficulties in the design and realization of such structural systems is how to grasp the relationship between the member length change, and their shapes and stress distribution. In this paper a simple and effective linear inverse theory for the control of such structures is described. Then its validity is proved by showing some numerical examples and the comparison with the results of real scale tests of a tensegrity frame.