Search

In seismic data analysis, acoustic impedance may be used for predicting the porosity of a sedimentary layer with known lithology, where the influence on impedance from pore fluid is known. This method typically assumes that the stiffness of the rock frame is solely a function of porosity and that there is no variation in stiffness between rocks with the same composition and porosity. However, contact cementation can cause an increase in rock frame stiffness with no significant reduction in porosity. This means that rocks with similar porosities and varying degrees of cementation would have significantly different elastic moduli; thus, attempting to derive porosity from seismic impedance alone would incur errors. The primary goal of this study is to assess how significant this error is for porosity prediction. The secondary goal is to calculate the dry rock frame stiffness across the field and evaluate to what extent the time of hydrocarbon placement affects the lithification of the rock. To fulfill these objectives, we studied petrophysical logs and core data from the Dan field and derived seismic impedance, porosity and stiffness expressed in terms of Biot’s coefficient for 14 wells across the field. Trends defined between acoustic impedance, porosity, and Biot’s coefficient in the Maastrichtian unit of the Dan field showed that, for a given porosity and water saturation in the oil zone, the chalk located at the crest of the reservoir is softer than the chalk located at the flanks of the field.We also found that chalk is stiffer in the down-faulted north-western part of the field for a given porosity and water saturation in the oil zone. In the water zone, no difference was found. We speculate that contact cementation downflank in the northwest block caused an increase in the chalk stiffness without reducing the porosity, while earlier hydrocarbon filling of the south-eastern block and the crest of the structure delayed contact cementation between the chalk particles. This difference in contact cementation between rocks of the same porosity causes porosity prediction from acoustic impedance to vary by up to six porosity units.

Full-waveform inversion (FWI) is today a standard process for the inverse problem of seismic imaging. PDE-constrained optimization is used to determine unknown parameters in a wave equation that represent geophysical properties. The objective function measures the misfit between the observed data and the calculated synthetic data, and it has traditionally been the least-squares norm. In a sequence of papers, we introduced the Wasserstein metric from optimal transport as an alternative misfit function for mitigating the so-called cycle skipping, which is the trapping of the optimization process in local minima. In this paper, we first give a sharper theorem regarding the convexity of the Wasserstein metric as the objective function. We then focus on two new issues. One is the necessary normalization of turning seismic signals into probability measures such that the theory of optimal transport applies. The other, which is beyond cycle skipping, is the inversion for parameters below reflecting interfaces. For the first, we propose a class of normalizations and prove several favorable properties for this class. For the latter, we demonstrate that FWI using optimal transport can recover geophysical properties from domains where no seismic waves travel through. We finally illustrate these properties by the realistic application of imaging salt inclusions, which has been a significant challenge in exploration geophysics.

A novel approach of unconventional reservoir stress evaluation is proposed to enhance the accuracy of fracture development prediction. Differential horizontal stress ratio can reflect the stress characteristic of fractured reservoir. To calculate the differential horizontal stress ratio more intuitionistic and simpler, we re‐derive its formula with Poisson's ratio and fracture density. Meanwhile, a new azimuthal PP‐wave amplitude versus offset equation based on Poisson's ratio and fracture density was derived for inverting above elastic parameters directly. Therefore, two steps are needed to realize stress evaluation. First, amplitude‐versus‐azimuthal angle inversion in the Bayesian framework in constraint of prior information such as well log data or rock physics information is executed for elastic and fracture parameters using pre‐stack angle gathers of different azimuth angles. The second procedure is to estimate differential horizontal stress ratio with Poisson's ratio and fracture density. Finally, a real data set is studied to test the new approach and the result demonstrated that the estimated differential horizontal stress ratio can reflect the stress property and it agrees with geological law and the new drilled well interpretation. Therefore, we can conclude that the combination of newly derived differential horizontal stress ratio and azimuthal PP‐wave amplitude versus offset equation in this study provides an available method for estimating the differential horizontal stress ratio of unconventional reservoirs, and the new approach can offer reliable geophysical information for stress evaluation.

Seismic inversion is a central step within the seismic reservoir characterization workflow. In seismic inversion, seismic amplitudes are inverted to predict the spatial distribution of the subsurface elastic properties. In subsequent steps of the geomodelling workflow, the inverted models are then converted into rock properties and facies. The methodology to perform the seismic inversion does have an impact on the predictions about the spatial distribution of the rock properties of interest. Each method has advantages and limitations depending on data quality, availability and the objective of the study. This work compares the application of deterministic and stochastic elastic inversion in a challenging real dataset. The deterministic approach is based on a sparse spike approximation, while the stochastic one is a global geostatistical seismic amplitude‐versus‐angle inversion. Both methods are applied to a real dataset acquired over a producing field located offshore Greece. Both the complex geological setting, where one target is close to seismic resolution, and data quality affect the performance of the seismic inversion methods. The results of both inversion methods are compared in terms of facies probability using Bayesian classification over the inverted elastic models and by comparing predictions to direct measurements at a blind well location. This application example shows the ability of deterministic inversion to retrieve an inverted solution with broader geological predictions which matched gross features in the wells. On the other hand, geostatistcal inversion was able to predict thin continous sand bodies, which correspond to the finer details at and below the seismic resolution of the observed data.

Ill‐posedness is one of the most common and intractable issues that arise when solving geophysical inverse problems. Ill‐posedness could be induced by various factors such as noise, band‐limited intrinsic property of seismic data and inappropriate forward operators. Regularization has been proven to be an effective method widely accepted for mitigating the adverse effects of ill‐posedness. Aiming to improve the stability and fidelity of the pre‐stack seismic inversion process, we implement the inversion in a Bayesian framework, with a logarithmic absolute criterion taken as a likelihood function, and an l1−2‐norm metric as a priori constraint. Here, we exploit the linear approximation as the forward operator, and optimize the regularized misfit function by the alternating direction method of multipliers. Applications of the method to synthetic and real data sets yielded improved inversion results in terms of accuracy and resolution, and demonstrated the robustness of the method to noise.

Analysis of amplitude variation with offset is an essential step for reservoir characterization. For an accurate reservoir characterization, the amplitude obtained with an isotropic assumption of the reservoir must be corrected for the anisotropic effects. The objective is seismic anisotropic amplitude correction in an effective medium, and, to this end, values and signs of anisotropic parameter differences (Δδ and Δe) across the reflection interfaces are needed. These parameters can be identified by seismic and well log data. A new technique for anisotropic amplitude correction was developed to modify amplitude changes in seismic data in transversely isotropic media with a vertical axis of symmetry. The results show that characteristics of pre‐stack seismic data, that is, amplitude variation with offset gradient, can be potentially related to the sign of anisotropic parameter differences (Δδ and Δe) between two layers of the reflection boundary. The proposed methodology is designed to attain a proper fit between modelled and observed amplitude variation with offset responses, after anisotropic correction, for all possible lithofacies at the reservoir boundary. We first estimate anisotropic parameters, that is, δ and e, away from the wells through Backus averaging of elastic properties resulted from the first pass of isotropic pre‐stack seismic inversion, on input data with no amplitude correction. Next, we estimate the anisotropic parameter differences at reflection interfaces (values and signs of Δδ and Δe). We then generate seismic angle gather data after anisotropic amplitude correction using Ruger's equation for the P‐P reflection coefficient. The second pass of isotropic pre‐stack seismic inversion is then performed on the amplitude‐corrected data, and elastic properties are estimated. Final outcome demonstrates how introduced methodology helps to reduce the uncertainty of elastic property prediction. Pre‐stack seismic inversion on amplitude‐corrected seismic data results in more accurate elastic property prediction than what can be obtained from non‐corrected data. Moreover, a new anisotropy attribute (ν) is presented for improvement of lithology identification.

Rock typing and flow unit detection are more challenging in clastic reservoirs with a uniform pore system. An integrated workflow based on well logs, inverted seismic data and rock physics models is proposed and developed to address such challenges. The proposed workflow supplies a plausible reservoir model for further investigation and adds extra information. Then, this workflow has been implemented in order to define different rock types and flow units in an oilfield in the Persian Gulf, where some of these difficulties have been observed. Here, rock physics models have the leading role in our proposed workflow by providing a diagnostic framework in which we successfully differentiate three rock types with variant characteristics on the given wells. Furthermore, permeability and porosity are calculated using the available rock physics models to define several flow units. Then, we extend our investigation to the entire reservoir by means of simultaneous inversion and rock physics models. The outcomes of the study suggest that in sediments with homogeneous pore size distribution, other reservoir properties such as shale content and cementation (which have distinct effects on the elastic domain) can be used to identify rock types and flow units. These reservoir properties have more physical insights for modelling purposes and can be distinguished on seismic cube using proper rock physics models. The results illustrate that the studied reservoir mainly consists of rock type B, which is unconsolidated sands and has the characteristics of a reservoir for subsequent fluid flow unit analysis. In this regard, rock type B has been divided into six fluid units in which the first detected flow unit is considered as the cleanest unit and has the highest reservoir process speed about 4800 to 5000 mD. Here, reservoir quality decreases from flow unit 1 to flow unit 6.

We compare two Monte Carlo inversions that aim to solve some of the main problems of dispersion curve inversion: deriving reliable uncertainty appraisals, determining the optimal model parameterization and avoiding entrapment in local minima of the misfit function. The first method is a transdimensional Markov chain Monte Carlo that considers as unknowns the number of model parameters, that is the locations of layer boundaries together with the Vs and the Vp/Vs ratio of each layer. A reversible‐jump Markov chain Monte Carlo algorithm is used to sample the variable‐dimension model space, while the adoption of a parallel tempering strategy and of a delayed rejection updating scheme improves the efficiency of the probabilistic sampling. The second approach is a Hamiltonian Monte Carlo inversion that considers the Vs, the Vp/Vs ratio and the thickness of each layer as unknowns, whereas the best model parameterization (number of layer) is determined by applying standard statistical tools to the outcomes of different inversions running with different model dimensionalities. This work has a mainly didactic perspective and, for this reason, we focus on synthetic examples in which only the fundamental mode is inverted. We perform what we call semi‐analytical and seismic inversion tests on 1D subsurface models. In the first case, the dispersion curves are directly computed from the considered model making use of the Haskell–Thomson method, while in the second case they are extracted from synthetic shot gathers. To validate the inversion outcomes, we analyse the estimated posterior models and we also perform a sensitivity analysis in which we compute the model resolution matrices, posterior covariance matrices and correlation matrices from the ensembles of sampled models. Our tests demonstrate that major benefit of the transdimensional inversion is its capability of providing a parsimonious solution that automatically adjusts the model dimensionality. The downside of this approach is that many models must be sampled to guarantee accurate posterior uncertainty. Differently, less sampled models are required by the Hamiltonian Monte Carlo algorithm, but its limits are the computational effort related to the Jacobian computation, and the multiple inversion runs needed to determine the optimal model parameterization.

We propose an optimized method to compute travel times for seismic inversion problems. It is a hybrid method combining several approaches to deal with travel time computation accuracy in unstructured meshes based on tetrahedral elementary cells. As in the linear travel time interpolation method, the proposed approach computes travel times using seismic ray paths. The method operates in two sequential steps: At a first stage, travel times are computed for all nodes of the mesh using a modified version of the shortest path method. The difference with the standard version is that additional secondary nodes (called tertiary nodes) are added temporarily around seismic sources in order to improve accuracy with a reasonable increase in computational cost. During the second step, the steepest travel time gradient method is used to trace back ray paths for each source–receiver pair. Travel times at each receiver are then recomputed using slowness values at the intersection points between the ray path and the traversed cells. A number of numerical tests with an array of different velocity models, mesh resolutions and mesh topologies have been carried out. These tests showed that an average relative error in the order of 0.1% can be achieved at a computational cost that is suitable for travel time inversion.

Tight oil siltstones are rocks with complex structure at pore scale and are characterized by low porosity and low permeability at macroscale. The production of tight oil siltstone reservoirs can be increased by hydraulic fracturing. For optimal fracking results, it is desirable to map the ability to fracture based on seismic data prior to fracturing. Brittleness is currently thought to be a key parameter for evaluating the ability to fracture. To link seismic information to the brittleness distribution, a rock physics model is required. Currently, there exists no commonly accepted rock physics model for tight oil siltstones. Based on the observed correlation between porosity and mineral composition and known microstructure of tight oil siltstone in Daqing oilfield of Songliao basin, we develop a rock physics model by combining the Voigt–Reuss–Hill average, self‐consistent approximation and differential effective medium theory. This rock physics model allows us to explore the dependence of the brittleness on porosity, mineral composition, microcrack volume fraction and microcrack aspect ratio. The results show that, as quartz content increases and feldspar content decreases, Young's modulus tends to increase and Poisson ratio decreases. This is taken as a signature of higher brittleness. Using well log data and seismic inversion results, we demonstrate the versatility of the rock physics template for brittleness prediction.

In this paper, we develop a general approach to integrating petrophysical models in three‐dimensional seismic full‐waveform inversion based on the Gramian constraints. In the framework of this approach, we present an example of the frequency‐domain P‐wave velocity inversion guided by an electrical conductivity model. In order to introduce a coupling between the two models, we minimize the corresponding Gramian functional, which is included in the Tikhonov parametric functional. We demonstrate that in the case of a single‐physics inversion guided by a model of different physical type, the general expressions of the Gramian functional and its gradients become simple and easy to program. We also prove that the Gramian functional has a non‐negative quadratic form, so it can be easily incorporated in a standard gradient‐based minimization scheme. The developed new approach of seismic inversion guided by the known petrophysical model has been validated by three‐dimensional inversion of synthetic seismic data generated for a realistic three‐dimensional model of the subsurface.

Inversion of seismic data using information from horizontal wells is often hampered by cumulative well-location errors. These errors can have a significant influence on the final subsurface model derived from the data. To achieve a proper data integration and arrive at correct uncertainty estimates, we formulate the problem in a fully probabilistic framework and present a numerical approach for improving subsurface imaging using uncertain well-log data and their uncertain locations as well as uncertain seismic data. The result is improved model error quantification in the seismic inversion process.

Seismic inverse modeling is a common method in reservoir prediction and it plays a vital role in the exploration and development of oil and gas. Conventional seismic inversion method is difficult to combine with complicated and abstract knowledge on geological mode and its uncertainty is difficult to be assessed. The paper proposes an inversion modeling method based on GAN consistent with geology, well logs, seismic data. GAN is a the most promising generation model algorithm that extracts spatial structure and abstract features of training images. The trained GAN can reproduce the models with specific mode. In our test, 1000 models were generated in 1 second. Based on the trained GAN after assessment, the optimal result of models can be calculated through Bayesian inversion frame. Results show that inversion models conform to observation data and have a low uncertainty under the premise of fast generation. This seismic inverse modeling method increases the efficiency and quality of inversion iteration. It is worthy of studying and applying in fusion of seismic data and geological knowledge., Comment: 22 pages,13 figures

In geophysical inverse problems, the posterior model can be analytically assessed only in case of linear forward operators, Gaussian, Gaussian mixture, or generalized Gaussian prior models, continuous model properties, and Gaussian‐distributed noise contaminating the observed data. For this reason, one of the major challenges of seismic inversion is to derive reliable uncertainty appraisals in cases of complex prior models, non‐linear forward operators and mixed discrete‐continuous model parameters. We present two amplitude versus angle inversion strategies for the joint estimation of elastic properties and litho‐fluid facies from pre‐stack seismic data in case of non‐parametric mixture prior distributions and non‐linear forward modellings. The first strategy is a two‐dimensional target‐oriented inversion that inverts the amplitude versus angle responses of the target reflections by adopting the single‐interface full Zoeppritz equations. The second is an interval‐oriented approach that inverts the pre‐stack seismic responses along a given time interval using a one‐dimensional convolutional forward modelling still based on the Zoeppritz equations. In both approaches, the model vector includes the facies sequence and the elastic properties of P‐wave velocity, S‐wave velocity and density. The distribution of the elastic properties at each common‐mid‐point location (for the target‐oriented approach) or at each time‐sample position (for the time‐interval approach) is assumed to be multimodal with as many modes as the number of litho‐fluid facies considered. In this context, an analytical expression of the posterior model is no more available. For this reason, we adopt a Markov chain Monte Carlo algorithm to numerically evaluate the posterior uncertainties. With the aim of speeding up the convergence of the probabilistic sampling, we adopt a specific recipe that includes multiple chains, a parallel tempering strategy, a delayed rejection updating scheme and hybridizes the standard Metropolis–Hasting algorithm with the more advanced differential evolution Markov chain method. For the lack of available field seismic data, we validate the two implemented algorithms by inverting synthetic seismic data derived on the basis of realistic subsurface models and actual well log data. The two approaches are also benchmarked against two analytical inversion approaches that assume Gaussian‐mixture‐distributed elastic parameters. The final predictions and the convergence analysis of the two implemented methods proved that our approaches retrieve reliable estimations and accurate uncertainties quantifications with a reasonable computational effort.

Reflection traveltime tomography has been used to describe subsurface velocity structures which, in practice, can be used as a background or initial model for pre‐stack depth migration or full waveform inversion. Conventional reflection traveltime tomography is performed by solving an optimization problem based on a ray‐tracing method. As a result, reflection traveltime tomography requires heavy computational efforts to carry out ray tracing and solve a large matrix equation. In addition, like most data‐domain tomography methods, reflection traveltime tomography depends on initial guesses and suffers from non‐uniqueness and uncertainty of solutions. In this research, we propose a deterministic ray‐based reflection traveltime tomography method by applying seismic interferometry. This method does not suffer from the non‐uniqueness problem and does not require a priori information on subsurface media. By adding a virtual layer (whose properties are known) on top of the real surface and applying convolution‐type interferometry, we approximately determine the stationary points (i.e., incident raypaths in the virtual layer). Then, we generate reflection points for a range of assumed velocities and estimate the velocity by considering the number of reflection points and the traveltime difference between the observed and calculated data. The reflection surface can then be recovered by using the estimated velocity. Once the first target layer is resolved, we can recover the whole media by recursively applying the same method to the lower layers. Numerical examples using surface seismic profile data for homogeneous‐layer (with a low‐velocity layer) and inhomogeneous‐layer models and real field data experiments on the Congo data set demonstrate that our method can successfully recover the velocities and depths of subsurface media without initial guesses. However, our method has some limitations for multi‐layer models because the method does not have sufficient reflection points for the deeper layers.

Transverse isotropy with a vertical axis of symmetry is a common form of anisotropy in sedimentary basins, and it has a significant influence on the seismic amplitude variation with offset. Although exact solutions and approximations of the PP‐wave reflection coefficient for the transversely isotropic media with vertical axis of symmetry have been explicitly studied, it is difficult to apply these equations to amplitude inversion, because more than three parameters need to be estimated, and such an inverse problem is highly ill‐posed. In this paper, we propose a seismic amplitude inversion method for the transversely isotropic media with a vertical axis of symmetry based on a modified approximation of the reflection coefficient. This new approximation consists of only three model parameters: attribute A, the impedance (vertical phase velocity multiplied by bulk density); attribute B, shear modulus proportional to an anellipticity parameter (Thomsen's parameter ε−δ); and attribute C, the approximate horizontal P‐wave phase velocity, which can be well estimated by using a Bayesian‐framework‐based inversion method. Using numerical tests we show that the derived approximation has similar accuracy to the existing linear approximation and much higher accuracy than isotropic approximations, especially at large angles of incidence and for strong anisotropy. The new inversion method is validated by using both synthetic data and field seismic data. We show that the inverted attributes are robust for shale‐gas reservoir characterization: the shale formation can be discriminated from surrounding formations by using the crossplot of the attributes A and C, and then the gas‐bearing shale can be identified through the combination of the attributes A and B. We then propose a rock‐physics‐based method and a stepwise‐inversion‐based method to estimate the P‐wave anisotropy parameter (Thomsen's parameter ε). The latter is more suitable when subsurface media are strongly heterogeneous. The stepwise inversion produces a stable and accurate Thomsen's parameter ε, which is proved by using both synthetic and field data.

Seismic inversion is an important tool that transfers interface information of seismic data to formation information, which renders the seismic data easily understood by geologists or petroleum engineers. In this study, a novel multi‐trace basis‐pursuit inversion method based on the Bayesian theory is proposed to enhance the vertical resolution and overcome the lateral instability of inversion results between different traces occasionally seen in the traditional trace‐by‐trace basis‐pursuit inversion method. The Markov process is initially introduced to describe the relationship between adjacent seismic traces and their correlation, which we then close couple in the equation of our new inversion method. A recursive function is further derived to simplify the inversion process by considering the particularity of the coefficient matrix in the multi‐trace inversion equation. A series of numerical‐analysis and field data examples demonstrates that both the traditional and the new methods for P‐wave impedance inversion are helpful in enhancing the resolution of thin beds that are usually difficult to discern from original seismic profiles, thus highlighting the importance of acoustic‐impedance inversion for thin bed interpretation. Furthermore, in addition to yielding thin bed inversion results with enhanced lateral continuity and high vertical resolution, our proposed method is robust to noise and cannot be easily contaminated by it, which we verify using both synthetic and field data.

We use Legendre polynomials to reparameterize geophysical inversions solved through a particle swarm optimization. The subsurface model is expanded into series of Legendre polynomials that are used as basis functions. In this framework, the unknown parameters become the series of expansion coefficients associated with each polynomial. The aim of this peculiar parameterization is threefold: efficiently decreasing the number of unknowns, inherently imposing a 1D spatial correlation to the recovered subsurface model and searching for maximally decoupled parameters. The proposed approach is applied to two highly non-linear geophysical optimization problems: seismic-petrophysical inversion and 1D elastic full-waveform inversion. In this work, with the aim to maintain the discussion at a simple level, we limit the attention to synthetic seismic data. This strategy allows us to draw general conclusions about the suitability of this peculiar parameterization for solving geophysical problems. The results demonstrate that the proposed approach ensures fast convergence rates together with accurate and stable final model predictions. In particular, the proposed parameterization reveals to be effective in reducing the ill conditioning of the optimization problem and in circumventing the so-called curse-of-dimensionality issue. We also demonstrate that the implemented algorithm greatly outperforms the outcomes of the more standard approach to global inversion in which each subsurface parameter is considered as an independent unknown.

During the time taken for seismic data to be acquired, reservoir pressure may fluctuate as a consequence of field production and operational procedures and fluid fronts may move significantly. These variations prevent accurate quantitative measurement of the reservoir change using 4D seismic data. Modelling studies on the Norne field simulation model using acquisition data from ocean-bottom seismometer and towed streamer systems indicate that the pre-stack intra-survey reservoir fluctuations are important and cannot be neglected. Similarly, the time-lapse seismic image in the post-stack domain does not represent a difference between two states of the reservoir at a unique base and monitor time, but is a mixed version of reality that depends on the sequence and timing of seismic shooting. The outcome is a lack of accuracy in the measurement of reservoir changes using the resulting processed and stacked 4D seismic data. Even for perfect spatial repeatability between surveys, a spatially variant noise floor is still anticipated to remain. For our particular North Sea acquisition data, we find that towed streamer data are more affected than the ocean-bottom seismometer data. We think that this may be typical for towed streamers due to their restricted aperture compared to ocean-bottom seismometer acquisitions, even for a favourable time sequence of shooting and spatial repeatability. Importantly, the pressure signals on the near and far offset stacks commonly used in quantitative 4D seismic inversion are found to be inconsistent due to the acquisition timestamp. Saturation changes at the boundaries of fluid fronts appear to show a similar inconsistency across sub-stacks. We recommend that 4D data are shot in a consistent manner to optimize aerial time coverage, and that additionally, the timestamp of the acquisition should be used to optimize pre-stack quantitative reservoir analysis.

If there are some erratic data (e.g. outliers), which may arise from measurement error, or other reasons, in seismic data, the seismic deconvolution and inversion need to be implemented in a way that minimizes their effects. However, the deconvolution and inversion methods based on L2-norm misfit function are highly sensitive to these erratic seismic observations. As an alternative, L1-norm misfit functions are more robust and erratic-resistant. In order to find the solution of the inverse problem constrained by an L1-norm misfit function, an iteratively re-weighted least squares algorithm is used frequently. However, it is relatively time consuming. In this paper, we propose a new method based on the sparse signal representation theory. The overcomplete dictionary used for the sparse representation of seismic data with erratic data is composed of two bases: a wavelet basis used for representing the seismic data to implement deconvolution and a Dirac basis used for representing the erratic data. In addition, at the stage of seismic inversion after deconvolution, total variation and a priori model are used as the regularization constraint terms to estimate inversion results with a blocky and laterally continuous structure. The new method is successfully tested on the noisy synthetic seismic data with erratic data. Finally, the proposed method is performed on a real seismic data section, and the inversion results are reasonable, i.e. consistent with the geologic structure of the original seismic data. Compared to the conventional sparse deconvolution and inversion method, the proposed method not only eliminates the effect of outliers, but also has highly improved computational efficiency.

Seismic inversion has drawn the attention of researchers due to its capability of building an accurate earth model. Such a model will need to be discretised finely, and the dimensions of the inversion problem will be very high. In this paper, we propose an efficient differential evolution algorithm and apply it to high-dimensional seismic inversion. Our method takes into account the differences among individuals, which are disregarded in conventional differential evolution methods, resulting to a better balance between exploration and exploitation. We divide the entire population into three subpopulations and propose a novel mutation strategy with two phases. Furthermore, we optimise the crossover operator by applying the components having the best objective function values into the crossover operator. We embed this strategy into a cooperative coevolutionary differential evolution and propose a new differential evolution algorithm referred to as a differential evolution with subpopulations. Then, we apply our scheme to both synthetic and field data; the results of highdimensional seismic inversion have shown that the proposed differential evolution with subpopulations achieves faster convergence and a higher-quality solution for seismic inversion.

Seismic inversion plays an important role in reservoir modelling and characterization due to its potential for assessing the spatial distribution of the subsurface petro-elastic properties. Seismic AVA inversion methodologies allow to retrieve P-wave and S-wave velocities and density individually allowing a better characterization of existing litho-fluid facies. We present an iterative geostatistical seismic AVA inversion algorithm that inverts pre-stack seismic data, sorted by angle gather, directly for: density; P-wave; and S-wave velocity models. The proposed iterative geostatistical inverse procedure is based on the use of stochastic sequential simulation and co-simulation algorithms as the perturbation technique of the model parameter space; and the use of a genetic algorithm as a global optimizer to make the simulated elastic models converge from iteration to iteration. All the elastic models simulated during the iterative procedure honour the marginal prior distributions of P-wave velocity, S-wave velocity and density estimated from the available well-log data, and the corresponding joint distributions between density versus P-wave velocity and P-wave versus S-wave velocity. We successfully tested and implemented the proposed inversion procedure on a pre-stack synthetic dataset, built from a real reservoir, and on a real pre-stack seismic dataset acquired over a deep-water gas reservoir. In both cases the results show a good convergence between real and synthetic seismic and reliable high resolution elastic subsurface Earth models. This article is protected by copyright. All rights reserved

Multiparameter inversion for pre-stack seismic data plays a significant role in quantitative estimation of subsurface petrophysical properties. However, it remains a complicated problem due to the non-unique results and unstable nature of the processing; the pre-stack seismic inversion problem is ill-posed and band-limited. Combining the full Zoeppritz equation and additional assumptions with edge-preserving regularisation can help to alleviate these problems. To achieve this, we developed an inversion method by constructing a new objective function that includes edge-preserving regularisation and soft constraints based on anisotropic Markov random fields and is intended especially for layered formations. We applied a fast simulated annealing algorithm to solve the nonlinear optimisation problem. The method directly obtains reflectivity RPP values using the full Zoeppritz equation instead of its approximations and effectively controls the stability of the multiparameter inversion by assuming a sectionally constant S- and P-wave velocity ratio and using the generalised Gardner equation. We substituted the inverted parameters, i.e., the P-wave velocity, the fitting deviation of S-wave velocity, and the density were inverted instead of the P-wave velocity, the S-wave velocity, and the density, and the generalised Gardner equation was applied as a constraint. Test results on two-dimensional synthetic data indicated that our substitution obtained improved results for multiparameter inversion. The inverted results could be improved by utilising high-order anisotropic Markov random field neighbourhoods at early stages and low-order anisotropicMarkov random field neighbourhoods in the later stages. Moreover, for layered formations, using a large horizontal weighting coefficient can preserve the lateral continuity of layers, and using a small vertical weighting coefficient allows for large longitudinal gradients of the interlayers. The inverted results of the field data revealed more detailed information about the layers and matched the logging curves at the wells acceptably over most parts of the curves.

Full waveform inversion for reection events is limited by its linearized update re-quirements given by a process equivalent to migration. Unless the background velocity model is reasonably accurate, the resulting gradient can have an inaccurate update direction leading the inversion to converge what we refer to as local minima of the objective function. In our approach, we consider mild lateral variation in the model, and thus, use a gradient given by the oriented time-domain imaging method. Specifically, we apply the oriented time-domain imaging on the data residual to obtain the geometrical features of the velocity perturbation. After updating the model in the time domain, we convert the perturbation from the time domain to depth using the average velocity. Considering density is constant, we can expand the conventional 1D impedance inversion method to 2D or 3D velocity inversion within the process of full waveform inversion. This method is not only capable of inverting for velocity, but it is also capable of retrieving anisotropic parameters relying on linearized representations of the reection response. To eliminate the cross-talk artifacts between different parameters, we utilize what we consider being an optimal parametrization for this step. To do so, we extend the prestack time-domain migration image in incident angle dimension to incorporate angular dependence needed by the multiparameter inversion. For simple models, this approach provides an efficient and stable way to do full waveform inversion or modified seismic inversion and makes the anisotropic inversion more practicable. The proposed method still needs kinematically accurate initial models since it only recovers the high-wavenumber part as conventional full waveform inversion method does. Results on synthetic data of isotropic and anisotropic cases illustrate the benefits and limitations of this method. This article is protected by copyright. All rights reserved

Seismic inversions produce seismic models, which are 3-dimensional (3D) images of wave velocity of the entire planet retrieved by fitting seismic measurements made on records of past earthquakes or other seismic events. Computing power of the TeraFlop era, along with the dataflow from new, very dense, seismic arrays, has led to a new generation of 3D seismic Earth models with an unprecedented level of resolution. Here we compare two recent models of western United States from the Dynamic North America (DNA) seismic imaging effort. The two models only differ in the wave propagation that was used for their inversion: one is based on ray theory (RT), and the other on finite frequency (FF). We evaluate the two models using an independent numerical method and statistical tests. We show that they differ in how they produce seismic signals from a subset of earthquakes that were used in the original inversion and were recorded on the US array. This is especially true for measurements done in the Yellowstone area which has a large negative seismic anomaly. This result is of importance for seismologists who have been debating on the practical benefit of using FF in ill-posed Earth inversions. Model evaluation, such as the one reported here, represents an opportunity for collaboration between geophysical and statistical communities. More opportunities should arise with the upcoming Exascale era, which will provide enough computational power to explore together several sources of errors in models with thousands of parameters, opening the way of uncertainty quantification of seismic models.

Rock fractures are of great practical importance to petroleum reservoir engineering because they provide pathways for fluid flow, especially in reservoirs with low matrix permeability, where they constitute the primary flow conduits. Understanding the spatial distribution of natural fracture networks is thus key to optimizing production. The impact of fracture systems on fluid flow patterns can be predicted using discrete fracture network models, which allow not only the 6 independent components of the second-rank permeability tensor to be estimated, but also the 21 independent components of the fully anisotropic fourth-rank elastic stiffness tensor, from which the elastic and seismic properties of the fractured rock medium can be predicted. As they are stochastically generated, discrete fracture network realizations are inherently non-unique. It is thus important to constrain their construction, so as to reduce their range of variability and hence the uncertainty of fractured rock properties derived from them. This paper presents the underlying theory and implementation of a method for constructing a geologically realistic discrete fracture network, constrained by seismic amplitude variation with offset and azimuth data. Several different formulations are described, depending on the type of seismic data and prior geologic information available, and the relative strengths and weaknesses of each approach are compared. Potential applications of the method are numerous, including the prediction of fluid flow, elastic and seismic properties of fractured reservoirs, model-based inversion of seismic amplitude variation with offset and azimuth data, and the optimal placement and orientation of infill wells to maximize production. This article is protected by copyright. All rights reserved

Geostatistical seismic inversion methods are routinely used in reservoir characterization studies because of their potential to infer the spatial distribution of the petro-elastic properties of interest (e.g., density, elastic and acoustic impedances) along with the associated spatial uncertainty. Within the geostatistical seismic inversion framework, the retrieved inverse elastic models are conditioned by a global probability distribution function and global spatial continuity model as estimated from the available well-log data for the entire inversion grid. However, the spatial distribution of the real subsurface elastic properties is complex, heterogeneous and in many cases non-stationary since they directly depend on the subsurface geology, i.e., the spatial distribution of the facies of interest. In these complex geological settings, application of a single distribution function and spatial continuity model is not enough to properly model the natural variability of the elastic properties of interest. In this study, we propose a 3D geostatistical inversion technique that is able to incorporate the reservoir's heterogeneities. This method uses a traditional geostatistical seismic inversion conditioned by local multi-distribution functions and spatial continuity models under non-stationary conditions. The procedure of the proposed methodology is based on a zonation criterion along the vertical direction of the reservoir grid. Each zone can be defined by conventional seismic interpretation, with the identification of the main seismic units and significant variations of seismic amplitudes. The proposed method was applied to a highly non-stationary synthetic seismic dataset with different levels of noise. The results of this work clearly show the advantages of the proposed method against conventional geostatistical seismic inversion procedures. It is important to highlight the impact of this technique in terms of a higher convergence between real and inverted reflection seismic data and the more realistic approximation towards the real subsurface geology comparing with traditional techniques. This article is protected by copyright. All rights reserved

We extend the frequency- and angle-dependent poroelastic reflectivity to systematically analyse the characteristic of seismic waveforms for highly attenuating reservoir rocks. It is found that the mesoscopic fluid pressure diffusion can significantly affect the root-mean-square amplitude, frequency content, and phase signatures of seismic waveforms. We loosely group the seismic amplitude-versus-angle and -frequency characteristics into three classes under different geological circumstances: (i) for Class-I amplitude-versus-angle and -frequency, which corresponds to well-compacted reservoirs having Class-I amplitude-versus-offset characteristic, the root-mean-square amplitude at near offset is boosted at high frequency, whereas seismic energy at far offset is concentrated at low frequency; (ii) for Class-II amplitude-versus-angle and -frequency, which corresponds to moderately compacted reservoirs having Class-II amplitude-versus-offset characteristic, the weak seismic amplitude might exhibit a phase-reversal trend, hence distorting both the seismic waveform and energy distribution; (iii) for Class-III amplitude-versus-angle and -frequency, which corresponds to unconsolidated reservoir having Class-III amplitude-versus-offset characteristic, the mesoscopic fluid flow does not exercise an appreciable effect on the seismic waveforms, but there exists a non-negligible amplitude decay compared with the elastic seismic responses based on the Zoeppritz equation.

Seismic wavelet estimation is an important part of seismic data processing and interpretation, whose reliability is directly related to the results of deconvolution and inversion. The methods for seismic wavelet estimation can be classified into two basic types: deterministic and statistical. By combining the deterministic spectral coherence method and the statistical skewness attribute method, the amplitude and phase of the time-varying wavelet are estimated separately. There is no assumption on the wavelet's time-independent nature or the phase characteristic. The advantage of this method is the ability to estimate time-varying phase. Phase-only corrections can then be applied by means of a time-varying phase rotation. Alternatively, amplitude and phase deconvolution can be achieved to enhance the resolution and support the fine reservoir prediction and description. We illustrate the method with both synthetic and real data examples. Synthetic examples certify its feasibility while real data example demonstrates the ability to estimate the time-varying characteristic of wavelets.

We present a new inversion method to estimate, from prestack seismic data, blocky P- and S-wave velocity and density images and the associated sparse reflectivity levels. The method uses the three-term Aki and Richards approximation to linearise the seismic inversion problem. To this end, we adopt a weighted mixed l2, 1-norm that promotes structured forms of sparsity, thus leading to blocky solutions in time. In addition, our algorithm incorporates a covariance or scale matrix to simultaneously constrain P- and S-wave velocities and density. This a priori information is obtained by nearby well-log data. We also include a term containing a low-frequency background model. The l2, 1 mixed norm leads to a convex objective function that can be minimised using proximal algorithms. In particular, we use the fast iterative shrinkage-thresholding algorithm. A key advantage of this algorithm is that it only requires matrix–vector multiplications and no direct matrix inversion. The latter makes our algorithm numerically stable, easy to apply, and economical in terms of computational cost. Tests on synthetic and field data show that the proposed method, contrarily to conventional l2- or l1-norm regularised solutions, is able to provide consistent blocky and/or sparse estimators of P- and S-wave velocities and density from a noisy and limited number of observations.

This paper presents a comparison between subsurface impedance models derived from different deterministic and geostatistical seismic inversion methodologies applied to a challenging synthetic dataset. Geostatistical seismic inversion methodologies nowadays are common place in both industry and academia, contrasting with traditional deterministic seismic inversion methodologies that are becoming less used as part of the geo-modelling workflow. While the first set of techniques allows the simultaneous inference of the best-fit inverse model along with the spatial uncertainty of the subsurface elastic property of interest, the second family of inverse methodology has proven results in correctly predicting the subsurface elastic properties of interest with comparatively less computational cost. We present herein the results of a benchmark study performed over a realistic three-dimensional non-stationary synthetic dataset in order to assess the performance and convergence of different deterministic and geostatistical seismic inverse methodologies. We also compare and discuss the impact of the inversion parameterisation over the exploration of the model parameter space. The results show that the chosen seismic inversion methodology should always be dependent on the type and quantity of the available data, both seismic and well-log, and the complexity of the geological environment versus the assumptions behind each inversion technique. The assessment of the model parameter space shows that the initial guess of traditional deterministic seismic inversion methodologies is of high importance since it will determine the location of the best-fit inverse solution.

Seismic inversion with prestack seismic data such as amplitude variation with offsets (AVO) inversion is an important tool in the estimation of elastic parameters for predicting lithology and discriminating fluid in conventional or unconventional hydrocarbon reservoirs. The product of Young's modulus and density (Young's impedance, YI) and the product of Poisson ratio and density (Poisson ratio impedance, PI) show great potential in lithology prediction and fluid discrimination of unconventional resources such as shale gas or oil. The high quality requirements for prestack data in density inversion render the estimation of YI and PI arduous and inaccurate with a conventional prestack inversion approach. In this study, a direct AVO inversion approach is proposed to estimate YI, PI, and density directly from P-wave seismic data. The linearized P-wave reflectivity approximate equation in terms of YI, PI, and density is initially derived. Five models, including four typical AVO classes, are utilized to verify the accuracy of the derived linearized P-wave reflectivity equation in comparison with the exact P-wave reflectivity equation and the frequently used linearized reflectivity approximate equation involving P- and S-wave velocities and density. Parameter sensitivity analysis illustrates that YI and PI can reasonably be estimated from P-wave reflectivity if a decorrelation scheme is utilized in the inversion algorithm. In addition, a pragmatic AVO inversion using a Bayesian scheme is suggested for the direct inversion of YI and PI from prestack seismic data. Synthetic and field data examples demonstrate the feasibility of the proposed inversion approach in the estimation of YI and PI and show the potential of this approach in fluid discrimination.

Coal elasticity is one of the important mechanical characteristics reflecting the material composition and structure of coal. Under the constraint of drilling and well logging, applying seismic exploration to predict elastic characteristics of the coal seams would be of important engineering significance for coal mining and CBM reservoir evaluation and development. In this regard, ultrasonic experiment of coal samples is the basis to realize the seismic inversion of physical properties of coal seam. In this paper, based on 30 pieces of six kinds of metamorphic raw coals collected from six mining areas of Yima, Fukang, Huainan, Pingdingshan, Hebi and Jiaozuo in China, laboratory ultrasonic measurements of these coal samples were respectively conducted in three directions of strike, dip and perpendicular to the beddings of coal seam under room temperature and pressure. The experimental results show:the samples' velocities of both pressure wave (P-wave) and shear wave (S-wave) decrease in turn in the three directions; there are obvious anisotropies of velocities, and the average anisotropy of P-wave velocity is stronger than S-wave's. Moreover, larger differences exist among both quality factors and modules of elasticity in the three directions, and the S-wave quality factors are greater than P-wave's; with exception of Poisson's ratio, modules of elasticity of coal are less than those common sedimentary rocks. Through the experiment and analysis, it can be demonstrated that Gardener and Castagna formulas are not suitable to represent relations of coal elastic parameters in China's coal fields, and two empirical formulas with higher precision are statistically given.

In this paper we propose Russell fluid factor direct estimation method based on stochastic seismic inversion. It is a Monte Carlo based strategy for non-linear inversion, which can effectively integrate the highfrequency information of well-logging data and have a higher resolution. And the method is formulated in a Bayesian framework. Firstly, we can calculate the Russell fluid factor using well-logging data and get the a priori information of fluid factor through sequential Gaussian simulation (SGS). Then we construct the likelihood function. Finally, we apply Metropolis algorithm in order to obtain an exhaustive description of the posteriori probability density. In this paper, we use the sequential Gaussian simulation (SGS) in a new implementation way, which can improve the computation speed. According to the numerical calculations, we can conclude that the final results match the model and real well-logging data well and have a higher resolution. In addition, Russell fluid factor we inverted is a sensitive indicator for reservoir fluid identification.

Seismic inversion is generally implemented with certain optimization algorithm. However, the inverse operator estimation algorithm proposed in this study is to perform the inversion of data matrix directly under the hypothesis that the inverse mapping exists in the empirically constrained subspaces. The key point of the proposed approach is to search those subspaces instead of searching for the solution indirectly as optimization algorithms do and it's more efficient. AVO/AVA (amplitude variation with offset or angle) inversion is widely utilized in exploration geophysics, and the inversion process is restricted by the quality of seismic data. L1 norm is applied in the construction of the kernel function of inversion by combining the constraint from initial models, which is helpful in enhancing the efficiency and stability of the inversion. Model and field data examples indicate that the proposed AVO inversion algorithm based on inverse operator estimation is more accurate and reliable.

We explore the link between basin modelling and seismic inversion by applying different rock physics models. This study uses the E-Dragon II data in the Gulf of Mexico. To investigate the impact of different rock physics models on the link between basin modelling and seismic inversion, we first model relationships between seismic velocities and both (1) porosity and (2) effective stress for well-log data using published rock physics models. Then, we build 1D basin models to predict seismic velocities derived from basin modelling with different rock physics models, in a comparison with average sonic velocities measured in the wells. Finally, we examine how basin modelling outputs can be used to aid seismic inversion by providing constraints for the background low-frequency model. For this, we run different scenarios of inverting near angle partial stack seismic data into elastic impedances to test the impact of the background model on the quality of the inversion results. The results of the study suggest that the link between basin modelling and seismic technology is a two-way interaction in terms of potential applications, and the key to refine it is establishing a rock physics models that properly describes changes in seismic signatures reflecting changes in rock properties.

The main difficulty in seismic full waveform inversion (FWI) is the strong nonlinearity, which is caused by the complexity of seismic wave propagation. Different components of elastic parameters result in different characteristics in seismic data. Meanwhile, different inversion quality is required in different stages of exploration or exploitation. So, there is no need to pursue matching all the seismic information during the inversion. Some problems can be solved by matching part of seismic information and the strong nonlinearity can also be avoided by this way. According to this consideration, a generalized FWI strategy and method based on seismic data subsets is presented. For different seismic data subsets, the gradients of the misfits have the same form and can be calculated by two times of seismic wave forward modeling just as traditional FWI. The only difference in calculating the gradients is the adjoint source. During the waveform inversion based on seismic data subsets, the responses of seismic data to different scales of perturbation on different media parameters should be analyzed intensively. Different seismic data subsets should be used in different stages of full waveform inversion. And then the residual of this seismic data subset is back-projected along the reasonable sub-kernels to decide where and which component of the model parameter need to be updated. As examples, envelope and reflection data subsets are used in FWI with synthetic and practical data to prove the validity and effectiveness of our presented FWI strategy and method. Especially, in the absence of low frequency contents, reasonable misfits can be used to recover the background compressional and shear velocity models using these FWI methods.

We obtain a three-dimensional (3-D) density structure in the crust and the upper-mantle beneath Northeast China through a joint inversion of gravity and P-wave travel time data. Compared to inversions of gravity or seismic data alone, the joint inversion can not only increase the vertical resolution of the gravity inversion, but also improve the reliability of the seismic inversion. Our inversion shows two results. First, the residual density anomalies in the crust and the upper-mantle correlate with the tectonics quite well: the orogen corresponds to low density whereas the basin corresponds to high density. Second, the areas beneath volcanoes have low density. One explanation is that the partial dehydrated Pacific slab was stagnant in the process of subducting into the mantle transition zone, and the hot material in the upper-mantle generated the upwelling magma and caused the volcano eruption.

We present a structural smoothing regularization scheme in the context of inversion of marine controlled-source electromagnetic data. The regularizing hypothesis is that the electrical parameters have a structure similar to that of the elastic parameters observed from seismic data. The regularization is split into three steps. First, we ensure that our inversion grid conforms with the geometry derived from seismic. Second, we use a seismic stratigraphic attribute to define a spatially varying regularization strength. Third, we use an indexing strategy on the inversion grid to define smoothing along the seismic geometry. Enforcing such regularization in the inversion will encourage an inversion result that is more intuitive for the interpreter to deal with. However, the interpreter should also be aware of the bias introduced by using seismic data for regularization. We illustrate the method using one synthetic example and one field data example. The results show how the regularization works and that it clearly enforces the structure derived from seismic data. From the field data example we find that the inversion result improves when the structural smoothing regularization is employed. Including the broadside data improves the inversion results even more, due to a better balancing between the sensitivities for the horizontal and vertical resistivities.

The hydrodynamic characterization of the epikarst, the shallow part of the unsaturated zone in karstic systems, has always been challenging for geophysical methods. This work investigates the feasibility of coupling time-lapse refraction seismic data with petrophysical and hydrologic models for the quantitative determination of water storage and residence time at shallow depth in carbonate rocks. The Biot–Gassmann fluid substitution model describing the seismic velocity variations with water saturation at low frequencies needs to be modified for this lithology. I propose to include a saturation-dependent rock-frame weakening to take into account water–rock interactions. A Bayesian inversion workflow is presented to estimate the water content from seismic velocities measured at variable saturations. The procedure is tested first with already published laboratory measurements on core samples, and the results show that it is possible to estimate the water content and its uncertainty. The validated procedure is then applied to a time-lapse seismic study to locate and quantify seasonal water storage at shallow depth along a seismic profile. The residence time of the water in the shallow layers is estimated by coupling the time-lapse seismic measurements with rainfall chronicles, simple flow equations, and the petrophysical model. The daily water input computed from the chronicles is used to constraint the inversion of seismic velocities for the daily saturation state and the hydrodynamic parameters of the flow model. The workflow is applied to a real monitoring case, and the results show that the average residence time of the water in the epikarst is generally around three months, but it is only 18 days near an infiltration pathway. During the winter season, the residence times are three times shorter in response to the increase in the effective rainfall.

Fluid depletion within a compacting reservoir can lead to significant stress and strain changes and potentially severe geomechanical issues, both inside and outside the reservoir. We extend previous research of time-lapse seismic interpretation by incorporating synthetic near-offset and full-offset common-midpoint reflection data using anisotropic ray tracing to investigate uncertainties in time-lapse seismic observations. The time-lapse seismic simulations use dynamic elasticity models built from hydro-geomechanical simulation output and a stress-dependent rock physics model. The reservoir model is a conceptual two-fault graben reservoir, where we allow the fault fluid-flow transmissibility to vary from high to low to simulate non-compartmentalized and compartmentalized reservoirs, respectively. The results indicate time-lapse seismic amplitude changes and travel-time shifts can be used to qualitatively identify reservoir compartmentalization. Due to the high repeatability and good quality of the time-lapse synthetic dataset, the estimated travel-time shifts and amplitude changes for near-offset data match the true model subsurface changes with minimal errors. A 1D velocity–strain relation was used to estimate the vertical velocity change for the reservoir bottom interface by applying zero-offset time shifts from both the near-offset and full-offset measurements. For near-offset data, the estimated P-wave velocity changes were within 10% of the true value. However, for full-offset data, time-lapse attributes are quantitatively reliable using standard time-lapse seismic methods when an updated velocity model is used rather than the baseline model.

In contrast to the use of marine seismic reflection techniques for reservoir-scale applications, where seismic inversion for quantitative sediment analysis is common, shallow-water, very-high-resolution seismic reflection data are seldom used for more than qualitative reflection interpretation. Here, a quantitative analysis of very-high-resolution marine seismic reflection profiles from a shallowwater (

In certain seismic data processing and interpretation tasks such as spiking deconvolution, tuning analysis, impedance inversion, and spectral decomposition, it is commonly assumed that the vertical direction is normal to reflectors. This assumption is false in the case of dipping layers and may therefore lead to inaccurate results. To overcome this limitation, we propose a coordinate system in which geometry follows the shape of each reflector and the vertical direction corresponds to normal reflectivity. We call this coordinate system stratigraphic coordinates. We develop a constructive algorithm that transfers seismic images into the stratigraphic coordinate system. The algorithm consists of two steps. First, local slopes of seismic events are estimated by plane‐wave destruction; then structural information is spread along the estimated local slopes, and horizons are picked everywhere in the seismic volume by the predictive‐painting algorithm. These picked horizons represent level sets of the first axis of the stratigraphic coordinate system. Next, an upwind finite‐difference scheme is used to find the two other axes, which are perpendicular to the first axis, by solving the appropriate gradient equations. After seismic data are transformed into stratigraphic coordinates, seismic horizons should appear flat, and seismic traces should represent the direction normal to the reflectors. Immediate applications of the stratigraphic coordinate system are in seismic image flattening and spectral decomposition. Synthetic and real data examples demonstrate the effectiveness of stratigraphic coordinates.

Seismic conditioning of static reservoir model properties such as porosity and lithology has traditionally been faced as a solution of an inverse problem. Dynamic reservoir model properties have been constrained by time-lapse seismic data. Here, we propose a methodology to jointly estimate rock properties (such as porosity) and dynamic property changes (such as pressure and saturation changes) from time-lapse seismic data. The methodology is based on a full Bayesian approach to seismic inversion and can be divided into two steps. First we estimate the conditional probability of elastic properties and their relative changes; then we estimate the posterior probability of rock properties and dynamic property changes. We apply the proposed methodology to a synthetic reservoir study where we have created a synthetic seismic survey for a real dynamic reservoir model including pre-production and production scenarios. The final result is a set of point-wise probability distributions that allow us to predict the most probable reservoir models at each time step and to evaluate the associated uncertainty. Finally we also show an application to real field data from the Norwegian Sea, where we estimate changes in gas saturation and pressure from time-lapse seismic amplitude differences. The inverted results show the hydrocarbon displacement at the times of two repeated seismic surveys.

Full-waveform inversion is re-emerging as a powerful data-fitting procedure for quantitative seismic imaging of the subsurface from wide-azimuth seismic data. This method is suitable to build high-resolution velocity models provided that the targeted area is sampled by both diving waves and reflected waves. However, the conventional formulation of full-waveform inversion prevents the reconstruction of the small wavenumber components of the velocity model when the subsurface is sampled by reflected waves only. This typically occurs as the depth becomes significant with respect to the length of the receiver array. This study first aims to highlight the limits of the conventional form of full-waveform inversion when applied to seismic reflection data, through a simple canonical example of seismic imaging and to propose a new inversion workflow that overcomes these limitations. The governing idea is to decompose the subsurface model as a background part, which we seek to update and a singular part that corresponds to some prior knowledge of the reflectivity. Forcing this scale uncoupling in the full-waveform inversion formalism brings out the transmitted wavepaths that connect the sources and receivers to the reflectors in the sensitivity kernel of the full-waveform inversion, which is otherwise dominated by the migration impulse responses formed by the correlation of the downgoing direct wavefields coming from the shot and receiver positions. This transmission regime makes full-waveform inversion amenable to the update of the long-to-intermediate wavelengths of the background model from the wide scattering-angle information. However, we show that this prior knowledge of the reflectivity does not prevent the use of a suitable misfit measurement based on cross-correlation, to avoid cycleskipping issues as well as a suitable inversion domain as the pseudo-depth domain that allows us to preserve the invariant property of the zero-offset time. This latter feature is useful to avoid updating the reflectivity information at each non-linear iteration of the full-waveform inversion, hence considerably reducing the computational cost of the entire workflow. Prior information of the reflectivity in the full-waveform inversion formalism, a robust misfit function that prevents cycle-skipping issues and a suitable inversion domain that preserves the seismic invariant are the three key ingredients that should ensure well-posedness and computational efficiency of fullwaveform inversion algorithms for seismic reflection data.