Inversion Of Seismic Reflection Data Research Articles

3D Seismic to absolute acoustic impedance inversion of the Lower Devonian Tadrart Sandstone in Ghadames Basin, NW Libya

The aim of this study is to map 3D acoustic impedance from 3D post-stacked seismic data, and then predict 3D porosity values from the acoustic impedance results for the Lower Devonian Tadrart reservoir sandstone using seismic inversion. Theoretically, seismic inversion is the process of transforming seismic reflection data into qualitative rock properties such as acoustic impedance and porosity, which is also very important for reservoir evaluation. Seismic inversion can be performed on pre or post stacked seismic data, it can be applied when the conventional seismic interpretation may become misleading under certain conditions. Most of oil and gas companies use the seismic inversion to improve the seismic interpretation by removing the side lobes and tuning effects from seismic data, to improve the estimation of rock properties, and to increase resolution and reliability. Inversion of seismic reflection data for various lithological and petro physical attributes is broadly used for reservoir characterization and hydrocarbons detection. Rock property related attributes are easier to interpret than the seismic reflectivity, which is related to boundaries between zones of contrasting acoustic impedance properties. Broadly, the commonalities between all impedance type properties are in their relations to the values measured from the seismic traces. The fundamental problem is the lack of low-frequency information in the seismic data leading to many uncertainties in the solution. The absolute acoustic impedance inversion is the method that has been applied in this study, and it is applied in Hamra field, northwest Ghadames basin, the approach of this work by extrapolating the well logs information into the seismic properties, this is, in turn, better estimations of reservoir properties such as porosity and an additional benefit that the interpretation efficiency is greatly improved.

Integration of Borehole Data in Geophysical Inversion Using Fuzzy Clustering

Borehole data is critical information to constrain geophysical inversions. Borehole information is usually used as a petrophysical constraint of the inversion scheme. Therefore, borehole information is only valuable if the borehole features are the same as the physical parameters of geophysical models. In fact, we may have many drilled holes, but the physical parameter of the geophysical model is only available in a few holes. Thus, the question is how we can exploit other borehole features to assist the inversion. We present an application of fuzzy clustering to incorporate multiple borehole features such as lithological, assay and wireline logs in the geophysical inversion. The integration of this extra information assist inversion process to build a model that fit surface geophysical data and simultaneously honour the prior information of borehole data. We apply this approach to the case study over the Kevitsa deposit within the Kevitsa ultramafic intrusion in northern Finland. The inversion of seismic reflection data with assistance from borehole information produces a more geologically interpretable image than seismic reflection data. The integration of both petrophysical and spatial attribute of the borehole enable the inversion to build a better model than only petrophysical constraint.

Optimised spectral merge of the background model in seismic inversion

The inversion of seismic reflection data to absolute impedance generates low-frequency deviations around the true impedance if the frequency content of the background impedance model does not merge seamlessly into the spectrum of the inverted seismic data. We present a systematic method of selecting a background model that minimises the mismatch between the background model and the relative impedance obtained by inverting the seismic data at wells. At each well a set of well-log relative impedances is formed by passing the impedance log through a set of zero-phase high-pass filters. The corresponding background models are constructed by passing the impedance log through the complementary zero-phase low-pass filters and a set of seismic relative impedances is computed by inverting the seismic data using these background models. If the inverted seismic data is to merge perfectly with the background model, it should correspond at the well to the well-log relative impedance. This correspondence is the basis of a procedure for finding the optimum combination of background model and inverted seismic data.It is difficult to predict the low-frequency content of inverted seismic data. These low frequencies are affected by the uncertainties in (1) measuring the low-frequency response of the seismic wavelet and (2) knowing how inversion protects the signal-to-noise ratio at low frequencies. Uncertainty (1) becomes acute for broadband seismic data; the low-frequency phase is especially difficult to estimate. Moreover we show that a mismatch of low-frequency phase is a serious source of inversion artefacts. We also show that relative impedance can estimate the low-frequency phase where a well tie cannot. Consequently we include a low-frequency phase shift, applied to the seismic relative impedances, in the search for the best spectral merge. The background models are specified by a low-cut corner frequency and the phase shifts by a phase intercept at zero frequency. A scan of the misfit between the well-log and seismic relative impedances as a function of corner frequency and phase intercept then reveals the most consistent combination of the background and relative impedances. The method is illustrated by application to a synthetic example and a broad-band data set.

A fast automatic multichannel blind seismic inversion for high-resolution impedance recovery

The inversion of seismic reflection data for acoustic impedance (AI) is a common and accepted method for the interpretation of poststack seismic data. The original mathematical problem is nonlinear due to the nonlinearity of relation between the AI and the earth reflectivity series and also the insufficiency of information about the source wavelet. Furthermore, the problem is ill posed due to the lack of low and high frequencies in the data. We have developed a multichannel blind inversion and solved it for obtaining the AI model and the wavelet directly from seismic reflection data. We found a solution to the overall problem by alternating between two subproblems, corresponding to the AI and wavelet recovery. Having an estimation of the wavelet, the algorithm directly inverts multichannel data for a high-resolution AI model, having blocky structures in the sense of total variation (TV) regularization, while satisfying a priori low-frequency information. The advantages of the split Bregman technique and the discrete cosine transform are used to build a fast and efficient algorithm for solving the nonlinear impedance inversion with TV regularization. Having an estimate of the AI model, the wavelet is updated by restricting it to have a sparse representation in a wavelet transformation domain while predicting the observed data. Numerical tests using simulated 2D data obtained from the benchmark Marmousi model and also 2D field data confirm that the proposed algorithm stably generates accurate estimates of complex AI models and complicated wavelets, simultaneously.

Full waveform inversion of seismic reflection data from the Forsmark planned repository for spent nuclear fuel, eastern central Sweden

The Swedish Nuclear Fuel and Waste Management Company (SKB) has been carrying out extensive studies at the planned repository for spent nuclear fuel at the Forsmark site in the eastern part of central Sweden since 2002. Identification of subhorizontal to gently dipping seismic reflections is especially important since these may represent transport routes for radionuclides. Studies have shown that such reflections can be generated by water filled fracture zones that have a lower velocity than the surrounding bedrock. Lithological changes, that is, mafic sills, may also be responsible for reflections in some cases. At the Forsmark site, it is difficult to distinguish fracture zones from mafic sills in the standard reflection seismic processed sections. However, since mafic sills usually have a positive velocity contrast with the background velocity field compared to fractures zones that have a negative one, the two possibilities could be differentiated if we could reconstruct the underground velocity field. Seismic full waveform inversion has the potential to perform this reconstruction, allowing us to discriminate between fractures zones and mafic sills. In this study, we apply a 2-D waveform inversion code on crooked line data sets acquired at the Forsmark site. This implies we are dealing with a 3-D geometry. We handle this problem by applying 3-D to 2-D coordinate projections. First, we perform a synthetic benchmark test with a similar geometry to that of the projected real data. We test both amplitude and phase inversion and phase only inversion on the synthetic data. The results show that the phase only inversion has fewer artefacts and is more stable. After successful application on the synthetic data, we apply the phase only waveform inversion on the real data. The resulting velocity fields show more details compared with the starting model based on first arrival traveltime tomography. Time domain synthetic data sets generated from the final velocity fields show a good match between the early arriving phase of the real data and the synthetic data sets. Checkerboard tests are also applied on the final models for resolution checks and indicate good resolution down to 200 m in the central parts of the profiles. Based on a comparison between the Kirchhoff pre-stack depth migrated sections and the final waveform inversion results we conclude that the reflections investigated most likely are generated at low-velocity fracture zones.

Seismic sparse-layer reflectivity inversion using basis pursuit decomposition

A basis pursuit inversion of seismic reflection data for reflection coefficients is introduced as an alternative method of incorporating a priori information in the seismic inversion process. The inversion is accomplished by building a dictionary of functions representing reflectivity patterns and constituting the seismic trace as a superposition of these patterns. Basis pursuit decomposition finds a sparse number of reflection responses that sum to form the seismic trace. When the dictionary of functions is chosen to be a wedge-model of reflection coefficient pairs convolved with the seismic wavelet, the resulting reflectivity inversion is a sparse-layer inversion, rather than a sparse-spike inversion. Synthetic tests suggest that a sparse-layer inversion using basis pursuit can better resolve thin beds than a comparable sparse-spike inversion. Application to field data indicates that sparse-layer inversion results in the potentially improved detectability and resolution of some thin layers and reveals apparent stratigraphic features that are not readily seen on conventional seismic sections.

Elastic full waveform inversion of multi‐component OBC seismic data

ABSTRACTElastic full waveform inversion of seismic reflection data represents a data‐driven form of analysis leading to quantification of sub‐surface parameters in depth. In previous studies attention has been given to P‐wave data recorded in the marine environment, using either acoustic or elastic inversion schemes. In this paper we exploit both P‐waves and mode‐converted S‐waves in the marine environment in the inversion for both P‐ and S‐wave velocities by using wide‐angle, multi‐component, ocean‐bottom cable seismic data. An elastic waveform inversion scheme operating in the time domain was used, allowing accurate modelling of the full wavefield, including the elastic amplitude variation with offset response of reflected arrivals and mode‐converted events. A series of one‐ and two‐dimensional synthetic examples are presented, demonstrating the ability to invert for and thereby to quantify both P‐ and S‐wave velocities for different velocity models. In particular, for more realistic low velocity models, including a typically soft seabed, an effective strategy for inversion is proposed to exploit both P‐ and mode‐converted PS‐waves. Whilst P‐wave events are exploited for inversion for P‐wave velocity, examples show the contribution of both P‐ and PS‐waves to the successful recovery of S‐wave velocity.

The adjoint method in seismology: I. Theory

The adjoint method is a powerful tool in geophysics that permits the computation of the exact first derivative of a physical observable or an associated objective function with respect to its parameters. Typical observables are displacement fields and flow patterns at the surface. Possible geophysical parameters are density, viscosity or elastic coefficients. When the observable can be modelled by solving a differential equation, the computation of the derivative only consists in solving the forward problem and its adjoint problem. Therefore, the adjoint method is far more efficient than any finite difference approximation. Here we present a mathematical formalism that generalises the derivation of the adjoint problem. In order to connect to work by Tarantola [Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics 49 (8), 1259–1266] we first give a derivation of the adjoint equations for the scalar wave equation in two dimensions. As objective function we choose the L 2 distance between the modelled wave field and real data. In this case the adjoint problem coincides with the original forward problem, the only difference being that the adjoint field satisfies terminal rather than initial conditions. A numerical example in two dimensions demonstrates that the adjoint field focusses near the location of a parameter perturbation at the same time when the original wavefront reaches that location. Based on this simple example, we introduce a generalised formalism for the adjoint method. It is independent of the existence of Green’s functions and their spatio-temporal reciprocity relations. Moreover, the formalism applies to non-linear equations such as the Navier–Stokes equations. This may become important in mantle flow reconstructions. The source term of the adjoint equations depends only on the specific objective function. Choosing the objective function to coincide with the observable itself allows us the computation of Jacobians, i.e., the derivative of the observable with respect to the model parameters. To demonstrate the consistency of our formalism with earlier analyses, we consider the anisotropic elastic wave equation with attenuation, which is of major interest in seismology.

Evidence for fault weakness and fluid flow within an active low-angle normal fault.

Determining the composition and physical properties of shallow-dipping, active normal faults (dips < 35 degrees with respect to the horizontal) is important for understanding how such faults slip under low resolved shear stress and accommodate significant extension of the crust and lithosphere. Seismic reflection images and earthquake source parameters show that a magnitude 6.2 earthquake occurred at about 5 km depth on or close to a normal fault with a dip of 25-30 degrees located ahead of a propagating spreading centre in the Woodlark basin. Here we present results from a genetic algorithm inversion of seismic reflection data, which shows that the fault at 4-5 km depth contains a 33-m-thick layer with seismic velocities of about 4.3 km s(-1), which we interpret to be composed of serpentinite fault gouge. Isolated zones exhibit velocities as low as approximately 1.7 km s(-1) with high porosities, which we suggest are maintained by high fluid pressures. We propose that hydrothermal fluid flow, possibly driven by a deep magmatic heat source, and high extensional stresses ahead of the ridge tip have created conditions for fault weakness and strain localization on the low-angle normal fault.

Detailed structure of the top of the melt body beneath the East Pacific Rise at 9°40′N from waveform inversion of seismic reflection data

We have applied waveform inversion to multichannel seismic reflection data collected at the East Pacific Rise at 9°40′N in order to determine the precise velocity structure of the magma body causing the axial magma chamber reflection. Our analysis supports the idea of a molten sill as previously suggested from forward modeling of seismic data from this location. Our inverted solution has a 30‐m‐thick sill with a P wave seismic velocity of 2.6 km s−1. Although not well constrained by the data we believe that the S wave velocity in the sill is not significantly different from 0.0 km s−1. The low P‐ and S wave velocities in the sill imply that it contains less than 30% crystals. The molten sill is underlain by a velocity gradient in which the P wave velocity increases from 2.6 to 3.5 km s−1 over a vertical distance of 50‐m. The shape of our velocity‐depth profile implies that accretion of material to the roof of the sill is minor compared to accretion to the floor. The underlying velocity gradient zone may represent crystal settling under gravity. We suggest that only material from the 30‐m‐thick layer can erupt.

Velocity structure of a gas hydrate reflector at Ocean Drilling Program site 889 from a global seismic waveform inversion

Strong subhorizontal reflectors in the upper few hundred meters of oceanic sediments have been imaged by seismic reflection profiling adjacent to many continental margins. These reflectors are known as “bottom simulating reflectors” (BSRs) as they run roughly parallel to the seabed. BSRs generally have been interpreted as representing the base of a methane gas hydrate stability field. A series of boreholes was drilled in October–November 1992 during leg 146 of the Ocean Drilling Program to examine fluid discharge and its relationship to the BSR, in the Cascadia accretionary wedge. Independently, we performed a global waveform inversion of seismic reflection data collected near site 889. We use a multistep strategy for global waveform inversion of seismic reflection data. The strategy is based on the Bayesian inference theory, in which the a posteriori probability distribution is determined by the a priori probability distribution and the likelihood function relating the data with the model. Three different likelihood functions have been defined relating different parts of the data with different wavelengths of model parameter variation. First, global nonlinear search algorithms are used to estimate the large‐scale features of one‐dimensional velocity variations; then a nonlinear iterative search is used to retrieve the velocity structure in fine detail. At each stage, bounds on the model parameters are estimated. The strategy has been implemented in the slowness and intercept‐time domain. We find that the BSR corresponds to a 15–20 m zone where the velocity drops from about 1.8 km/s to a minimum of 1.4 km/s, suggesting the presence of free methane gas in the low‐velocity zone, and a 15–20 m transition zone where the velocity gradually increases to a velocity of the surrounding sediments. Thus the total thickness of the low‐velocity zone is about 30–40 m. A small or negligible velocity increase just above the BSR, suggests that the hydrate fills less than 10% of the pore space in the sediments. Our result was confirmed by the drilling. The velocity derived using a vertical seismic profile and sonic log from the drill hole agree closely with our result. Waveform inversion can recover detailed and well‐constrained velocity variations from a high‐amplitude target such as a BSR; combined with the detailed interpretation of these velocity variations which the results from site 889 will provide, the inversion approach is a powerful tool for analysis of BSRs elsewhere.

Seismic mapping of the marmousi data set with the common reflecting element method (CRE method)

The inversion of seismic reflection data collected across a laterally heterogeneous Earth is an important objective in reflection seismology. The purpose of this paper is to illustrate the potential of the common reflecting element (CRE) imaging method, a new topological technique that involves the stacking of multi-offset reflection seismic data. The CRE method has four advantages: 1. (1) it is based on ray theory and, thus, is an elementary geometrical method, 2. (2) it only requires knowledge of seismic velocities at a “reference level” (near-surface) as a-priori information, 3. (3) it delivers wavefront parameters (radius of curvature and angle of incidence) that can be used to formulate new inversion methods, and 4. (4) normal-moveout stretching is avoided. The local dip of common subsurface reflecting elements (CRE's) is used as a criterion for choosing traces to be stacked in order to approximate “true” zero-offset sections. Since the method includes a prestack coherency analysis it may be regarded as a data-adaptive stacking technique. In the present study, the CRE method is applied to a complex synthetic data set (Marmousi model) to investigate its applicability under “ideal field conditions”. Our results show that the CRE algorithm, which works automatically, is stable and yields reasonable images. Information on the kinematic properties of the stacked reflected wavefronts is also produced.

Velocity Structure of a Gas Hydrate Reflector

Seismic reflection profiles across many continental margins have imaged bottom-simulating reflectors (BSRs) parallel to the seabed; these are often interpreted as the base of a zone in which methane hydrate "ice" is stable. Waveform inversion of seismic reflection data can be used to estimate from seismic data worldwide the velocity structure of a BSR and its thickness. A test of this method at a drill site of the Ocean Drilling Program predicts that sediment pores beneath the BSR contain free methane for approximately 30 meters. The hydrate and underlying gas represent a large global reservoir of methane, which may have economic importance and may influence global climate.

INVERSION OF POST‐STACK SEISMIC DATA USING SIMULATED ANNEALING1

AbstractModel‐based inversion of seismic reflection data is a global optimization problem when prior information is sparse. We investigate the use of an efficient, global, stochastic optimization method, that of simulated annealing, for determining the two‐way traveltimes and the reflection coefficients.We exploit the advantage of an ensemble approach to the inversion of full‐scale target zones on 2D seismic sections.In our ensemble approach, several copies of the model‐algorithm system are run in parallel. In this way, estimation of true ensemble statistics for the process is made possible, and improved annealing schedules can be produced.It is shown that the method can produce reliable results efficiently in the 2D case, even when prior information is sparse.

A simultaneous inversion for background velocity and impedance maps

Full‐waveform inversion of seismic reflection data is highly nonlinear because of the irregular form of the function measuring the misfit between the observed and the synthetic data. Since the nonlinearity results mainly from the parameters describing seismic velocities, an alternative to the full nonlinear inversion is to have an inversion method which remains nonlinear with respect to velocities but linear with respect to impedance contrasts. The traditional approach is to decouple the nonlinear and linear parts by first estimating the background velocity from traveltimes, using either traveltime inversion or velocity analysis, and then estimating impedance contrasts from waveforms, using either waveform inversion or conventional migration. A more sophisticated strategy is to obtain both the subsurface background velocities and impedance contrasts simultaneously by using a single least‐squares norm waveform‐fit criterion. The background velocity that adequately represents the gross features of the medium is parameterized using a sparse grid, whereas the impedance contrasts use a dense grid. For each updated velocity model, the impedance contrasts are computed using a linearized inversion algorithm. For a 1-D velocity background, it is very efficient to perform inversion in the f-k domain by using the WKBJ and Born approximations. The method performs well both with synthetic and field data.

A strategy for nonlinear elastic inversion of seismic reflection data

The problem of interpretation of seismic reflection data can be posed with sufficient generality using the concepts of inverse theory. In its roughest formulation, the inverse problem consists of obtaining the Earth model for which the predicted data best fit the observed data. If an adequate forward model is used, this best model will give the best images of the Earth’s interior. Three parameters are needed for describing a perfectly elastic, isotropic, Earth: the density ρ(x) and the Lamé parameters λ(x) and μ(x), or the density ρ(x) and the P-wave and S-wave velocities α(x) and β(x). The choice of parameters is not neutral, in the sense that although theoretically equivalent, if they are not adequately chosen the numerical algorithms in the inversion can be inefficient. In the long (spatial) wavelengths of the model, adequate parameters are the P-wave and S-wave velocities, while in the short (spatial) wavelengths, P-wave impedance, S-wave impedance, and density are adequate. The problem of inversion of waveforms is highly nonlinear for the long wavelengths of the velocities, while it is reasonably linear for the short wavelengths of the impedances and density. Furthermore, this parameterization defines a highly hierarchical problem: the long wavelengths of the P-wave velocity and short wavelengths of the P-wave impedance are much more important parameters than their counterparts for S-waves (in terms of interpreting observed amplitudes), and the latter are much more important than the density. This suggests solving the general inverse problem (which must involve all the parameters) by first optimizing for the P-wave velocity and impedance, then optimizing for the S-wave velocity and impedance, and finally optimizing for density. The first part of the problem of obtaining the long wavelengths of the P-wave velocity and the short wavelengths of the P-wave impedance is similar to the problem solved by present industrial practice (for accurate data interpretation through velocity analysis and “prestack migration”). In fact, the method proposed here produces (as a byproduct) a generalization to the elastic case of the equations of “prestack acoustic migration.” Once an adequate model of the long wavelengths of the P-wave velocity and of the short wavelengths of the P-wave impedance has been obtained, the data residuals should essentially contain information on S-waves (essentially P-S and S-P converted waves). Once the corresponding model of S-wave velocity (long wavelengths) and S-wave impedance (short wavelengths) has been obtained, and if the remaining residuals still contain information, an optimization for density should be performed (the short wavelengths of impedances do not give independent information on density and velocity independently). Because the problem is nonlinear, the whole process should be iterated to convergence; however, the information from each parameter should be independent enough for an interesting first solution.

MIGRATION STRATEGY*

ABSTRACTThe geophysical data processor today has on offer a great variety of tools for the inversion of seismic reflection data to estimate geological structure. The major subset of these comprises migration procedures, which span a wide range of sophistication and cost in terms both of computation time and manual effort on the part of interpreters and processing staff. The choice of an over‐powerful process can be very wasteful, but on the other hand too naive a migration procedure can lead to wrong interpretations which are much more costly still.Complete inversion procedures which aim to delineate all changes in rock densities and elastic properties in the subsurface are still in the imaginative stages of research. Not even the most sophisticated migration procedure in current use with real data, however, provides a complete inversion, but all depend in some measure on prior knowledge of the velocity structure of the section of the earth traversed by the seismic energy. Such knowledge may be very approximate at first, but each inversion should, through the skill of the interpreter, allow him to revise his velocity model and, up to some limit imposed by the quality and ambiguity of the original data, to improve the next inversion. Paradoxically, he can often be helped by using forward modeling procedures to check the implications of his ideas in the data domain, both in deciding how to update the velocity model and in selecting the most appropriate migration process to use next.We review here the currently available toolkit of migration and modeling processes and make suggestions as to how each process can fit into a learning strategy which can improve the interpretation as economically as possible and in as many iterative steps as the complexity of the earth's velocity structure makes necessary. An example is shown of the strategy being used in a complex overthrust region.The authors wish to thank the Chairman and Board of Directors of BP Exploration Co. for permission to publish this paper, and also make acknowledgment to our colleagues whose labours in research and development have made available to our use many of the essential tools required to implement the strategies we describe.

Inversion of seismic reflection data in the acoustic approximation

The nonlinear inverse problem for seismic reflection data is solved in the acoustic approximation. The method is based on the generalized least‐squares criterion, and it can handle errors in the data set and a priori information on the model. Multiply reflected energy is naturally taken into account, as well as refracted energy or surface waves. The inverse problem can be solved using an iterative algorithm which gives, at each iteration, updated values of bulk modulus, density, and time source function. Each step of the iterative algorithm essentially consists of a forward propagation of the actual sources in the current model and a forward propagation (backward in time) of the data residuals. The correlation at each point of the space of the two fields thus obtained yields the corrections of the bulk modulus and density models. This shows, in particular, that the general solution of the inverse problem can be attained by methods strongly related to the methods of migration of unstacked data, and commercially competitive with them.