Recovery of the acoustic impedance from seismic data is one of the most important steps in the reflection seismology (Oldenburg et al., 1983). Seismic inversion is a technique for recovering the acoustic impedance of the subsurface from the data measured at the surface known as seismic data (Russell, 1988). There are different methods for seismic inversion which are classified into two main groups: deterministic methods and stochastic (geostatistical) methods. The deterministic methods of seismic inversion are based on the minimization of the difference between two terms. The first term is the synthetic seismogram computed from the convolution of the reflectivity series (obtained from an initial estimated acoustic impedance model) with an extracted wavelet, and the second one is the real seismic profile (Francis, 2005). Due to the band-limited nature of the seismic data, the results of deterministic methods are smooth maps of acoustic impedance and may be far from of the reservoir facts. Accordingly, deterministic inversion methods are not capable of modeling heterogeneities of the reservoir and are not therefore appropriate for reservoir characterization (for example, the volumetric calculations of oil in place and the simulation of fluid flow in the reservoir) (Francis, 2006). The stochastic methods of seismic inversion are based on the sampling of the probability density functions (PDF’s). In other words, PDF’s of acoustic impedance are computed in different locations constraining well data and the variograms of the wells. The sampling of the resulted PDF’s is performed using a Monte Carlo method based on a random number generator (Francis, 2006). Consequently, stochastic inversion methods, unlike deterministic ones generate multiple realizations of acoustic impedance with the same probability. Furthermore, stochastic inversion produces high-resolution maps of the acoustic impedance because the spatial continuity models (variograms) control the frequency content of stochastic inversion results (Azevedo and Soares, 2017).

The stochastic inversion of seismic data was first introduced by Hass and Dubrule (1994). They considered well data as primary variable and implemented a sequential Gaussian simulation (SGS) technique in a random path of the trace locations. Then, in each trace location, a local optimization was performed. In another way, in each trace location, the synthetic trace computed from the simulated acoustic impedance model was compared with the real seismic trace, and the model of acoustic impedance was accepted only if the correlation coefficient between the compared terms were sufficient. Debeye et al. (1996) modified this technique using simulated annealing algorithm. Buland and Omre (2003) developed a new Bayesian technique for amplitude-variation-with-offset (AVO) inversion. They analytically computed posterior distributions of P-wave velocity, S-wave velocity, and density. In addition, they expressed that the analytical form of the posterior distribution deduces a fast stochastic inversion from the computational point of view. Buland et al. (2003) developed stochastic AVO inversion in the Fourier domain. They demonstrated that applying the inversion in this domain is computationally much faster than doing it in the time domain, but it has two limitations. The first restriction is the regular sampling of data and model in the time and space dimension, and the second one is the non-stationarity of the covariance function. Escobar et al. (2006) presented a method to overcome these two limitations. They used AVO inversion in a Bayesian framework according to the work of Buland and Omre (2003) and computed the log-Gaussian joint posterior distribution of the acoustic impedance and shear impedance. Then, an SGS technique was used to decompose the global PDF to local PDF’s in each trace location, and multiple realizations were generated by the sampling of the local PDF’s. Ansari et al. (2014) developed an approach for geostatistical inversion based on spectral geostatistical simulation. They applied their method to an Iranian gas-field in the Persian Gulf basin and concluded that the geostatistical seismic inversion can recover the low-frequency components. Sabeti et al. (2017) investigated the influence of the non-stationary patterns on geostatistical seismic inversion and showed that, in the complex geological scenarios, geostatistical inversion for non-stationary patterns provides a better result compared to conventional methods.

One of the main advantages of stochastic inversion is its ability to recover the thin beds. Widess (1973) defined the thin bed as one the thickness of which is less than one-eighth of predominant wavelength calculated using the thin bed velocity. Sancevero et al. (2005) applied both stochastic and deterministic inversions to Campos basin in Brazil. They concluded that the deterministic method reproduces thick bodies, while stochastic inversion shows details and heterogeneities and provides the possibility to evaluate the uncertainties associated with the results. McCrank et al. (2009) used geostatistical inversion for the recovery of thinly bedded Ardley Coals and showed that the geostatistical method, as opposed to the deterministic one, is able to map heterogeneities and assess uncertainties. Ziyisyian and Vocaturo (2015) applied stochastic inversion to San Jorge basin in Argentina. They found out the capability of stochastic inversion to provide high detailed impedance simulation for modeling individual thin sands. Yakovleva et al. (2016) discussed a comparison between deterministic and stochastic seismic inversion methods. They applied an AVO geostatistical inversion to seismic data to find out the thin coal beds in a complex gas saturated thin beds. Finally, they concluded that reservoir characterization would be highly efficient when using the geostatistical inversion method. Chen et al. (2017) combined statistical rock physics and geostatistical seismic inversion to predict thin reservoir sandstone. Then, they used Bayes theory to get facies distribution within the reservoir. Their workflow showed a successful prediction of the thin reservoir which is consistent with well log data and geology. Ghosh et al. (2018) presented a workflow for reservoir characterization incorporating rock physics and stochastic inversion. They used rock physics model for shear velocity estimation and lithofacies differentiation; then, they applied stochastic inversion to seismic data for mapping lithofacies distribution. The result of the stochastic inversion, as opposed to deterministic results, delineated the lateral changes of lithofacies and provided a high resolution of lithofacies map.

Common methods of stochastic inversion which are based on geostatistical simulation techniques suffer from the limitation of slowness because the new grids will be conditioned to previous estimated locations. Bayesian principle provides the possibility to calculate the global PDF in an analytical form, which saves the time. Hence, in this research, a Bayesian linearized AVO inversion is applied to an oilfield in the Persian Gulf. The sampling is performed after the decomposition of global PDF to a number of local PDF’s by the SGS technique, and multiple realizations of elastic properties are generated; finally, the capability of this method in detecting thin beds is discussed.

2. Methodology

Aki and Richards (1980) presented an approximation to the Zoeppritz equations by assuming weak contrast in the elastic properties of two mediums on the both sides of the interface. They showed that PP reflection coefficient (P-wave reflection coefficient from a P-wave incident) changes with the incident angle as given below:

where,  is the incident angle (or reflection angle); , , and  are the differences of the P-wave velocities, S-wave velocities, and densities on the both sides of the reflector respectively; , , and  are the averages of the P-wave velocities, S-wave velocities, and densities over the reflecting interface respectively. Coefficients , , and  are defined as:

and

Acoustic impedance and shear impedance are obtained by the multiplication of density and corresponding velocity. Fatti et al. (1994) rewrote Equation 1 in terms of acoustic impedance (I) and shear impedances (J):

where,  and  are PP and PS reflection coefficients at the normal incident. The third coefficient in Equation 5 is written as follows:

Fatti et al. (1994) expressed that Equation 6 approaches zero at low and moderate angles. The normal incident coefficients in Equation 5 can be written logarithmically, so Equation 5 is changed to:

Acoustic impedance and shear impedance are assumed to be log-Gaussian. Also, differences in these parameters are Gaussian due to the linearity of the differentiation operator (Buland and Omre, 2003).  is defined as a discrete matrix comprising of the acoustic impedance and shear impedance at the column i of layers:

Reflection coefficients for a set of n incident angles are computed as reads (Escobar, 2006; Buland and Omre, 2003):

where,  is the difference of ;  (the matrix of reflectivity series) and  (the matrix of coefficients) are defined by:

and

The inversion problem is solved in a Bayesian framework. Bayesian principle is expressed as:

where,  is the posterior distribution of elastic properties conditioned by the observed seismic data;  is the PDF of seismic data given that the elastic parameters of m  is known as likelihood function;  is the prior distribution of elastic properties, and  is the unconditional PDF of seismic data obtained by (Buland and Omre, 2003):

It is assumed that the PDF of the prior model m is log-Gaussian and can be written as (Escobar et al., 2006):

where,  and  are the mean and covariance function of the prior model respectively; the superscript T denotes the transpose of the matrix.

Escobar et al. (2006) assumed the seismic noise is a random field in space and time dimensions. They calculated the likelihood function by:

where,  is a matrix containing both coefficients of  and components of the extracted wavelet;  gives the synthetic trace in the trace location i and the incident angle of ;  is the noise covariance matrix, which can be estimated from the well data and the seismic trace in the well locations (Buland and Omre, 2003).

The conditional Gaussian posterior distribution of elastic properties is obtained by combining Equations 11 to 14. It is difficult to sample the global posterior PDF directly due to its high multidimensionality. Therefore, an SGS technique is employed to decompose the global PDF to the local PDF’s in each trace location conditioned by the previously simulated traces as written below (Escobar et al., 2006):

where,  is the global PDF;  is the local PDF at trace 1 conditioned by seismic trace 1;  is the local PDF at trace 2 conditioned by the simulated trace 1 and seismic trace 2, and  is the local PDF at trace N conditioned by simulated traces 1 to N-1 and seismic trace N.

Figure 1 shows the workflow of the stochastic seismic inversion described above. The proposed method of stochastic inversion is based on the Bayesian principle to produce the probability density function of the elastic properties such as acoustic impedance and shear impedance. To reach this goal, a prior model of elastic properties, a likelihood function, and the unconditional PDF of seismic data are needed. The first term taking part in Bayesian interference is the prior model, which is constructed using well logs and some horizons. Geostatistical techniques like kriging methods are appropriate for estimating an initial model of elastic properties. The spatial continuity models (variograms) of elastic properties are required for geostatistical interpolation and are calculated from known values in well locations. A kriging method is applied between horizons using variograms, and elastic properties measured for well logs and the prior models of acoustic impedance and shear impedance are generated. The second term participating in Bayesian interference is likelihood function which is calculated according to Equation 14 and is dependent on seismic data, prior models, and a seismic wavelet. Extracting the wavelet from seismic is required to calculate the likelihood function. The third term, the unconditional PDF of seismic data, is obtained using Equation 12. Combining Equations 11 to 14 provides the joint posterior distribution of elastic properties. The sampling of the high dimensional PDF is difficult, so it is better to decompose the global PDF to a number of local PDF’s in each trace location. SGS technique is a geostatistical simulation which is performed to generate the local PDF’s from a global PDF; the next step is to sample the PDF’s. It is conducted in a random seed, and different models of elastic properties called realizations are generated. Multiple realizations of acoustic impedance and shear impedance can be provided by sampling the PDF’s. Finally, the assessment of uncertainties is possible by calculating the variances of the elastic properties in the different locations.

3. Field data example

The stochastic method of inversion in a Bayesian framework was applied to one of the Persian Gulf oilfields located in the west of the Persian Gulf, south of Iran. This field produces oil from an Oligo-Miocene sandstone formation. Uniform and continuous sedimentation with little tectonic activity occurred from Cambrian to early Cretaceous when the platform experienced compressive tectonic stresses from north, northeast, and east in the areas of the Taurus, Zagros, and Oman mountain ranges. These Alpine orogenic movements resulted in gentle north-south trending structures. Shallow marine deposition and subsidence continued in the Persian Gulf until the intense thrusting and folding in the Taurus-Zagros-Oman orogenic belt during late Tertiary. The studied field owes its present structural configuration to the superposition of the resulting northwest-southwest folding over the earlier north-south trend. Sedimentological study of the field describes the reservoir formation as about 250 feet of unconsolidated sands with some thin carbonate beds embedded within the reservoir. The sands are very finely grained with minor amounts of feldspar. The reservoir is not homogeneous, and the presence of thin carbonate beds within the reservoir in some areas of the field will certainly reduce the vertical permeability locally. Fluid flow simulation requires a model including details and heterogeneities of the reservoir. Because of these heterogeneities and requiring an appropriate model for fluid flow simulation, this field has been selected to be studied.

Figure 1

The workflow of stochastic seismic inversion in a Bayesian framework.

There are three wells in the studied area which contain petrophysical logs such as P-wave velocity, S-wave velocity, gamma ray, density etc. Check shot corrections were used to correct the P-wave velocity logs in all well locations. S-wave velocity has been measured in one of the wells and is estimated in the other two wells using Greenberg-Castagna’s equation, which has the following form:

where,  and  are S-wave velocity and P-wave velocity logs respectively; coefficients , , and  have various values for different minerals or structures; they are respecified in an iterative process to estimate the S-wave velocity log.

The south portion of the seismic survey was used in the inversion. Seismic data quality plays an important role in data analysis and directly affects the results of inversion process. The data used in this study are 40-fold prestack seismic data with the dominant frequency of 35 Hz. The quality of the seismic data of the reservoir is superior because the reservoir is not deep, and the amplitude will be relatively preserved as opposed to deep reservoirs. However, processing steps should be applied to seismic data to make them interpretable. For AVO inversion, amplitude recovery steps are very important and should be performed carefully to show the true changes of amplitude with offset. Processing steps, including statics corrections, noise attenuation, deconvolution, demultipling, amplitude recovery, NMO corrections, and velocity analysis were applied to the seismic data. In the amplitude recovery step, the effects of geometrical spreading and adsorption were compensated to recover the true amplitudes of seismic data as much as possible. The velocity model obtained in the velocity analysis step was used to transform seismic data from the offset-time domain to the angle-time domain. The range of the angles was between 0 and 35 degrees, so it was suitable to classify it into two parts: near-angle part and middle-angle part. Statistical wavelets were extracted from seismic traces for two ranges of angles: 0 to 15 degrees (near-angle) and 15 to 35 degrees (middle-angle). The results of wavelet extraction showed that frequency content of near angle wavelet is a little more than that of the middle angle one. After wavelet extraction, correlating synthetic seismic traces and real seismic traces were performed in the well locations by shifting the P-wave velocity logs in the time dimension.

The prior model of acoustic impedance and shear impedance was constructed using corrected P-wave velocity logs, S-wave velocity logs, and density logs by means of the kriging interpolation method. The constructed prior model and seismic data took part in a Bayesian framework and reproduced the joint posterior distribution of the elastic properties. Then, an SGS technique was performed to decompose the global PDF to a number of local PDF’s. The sampling of the local PDF’s was performed in a random seed, and two hundred realizations of the acoustic impedance and shear impedance were generated. Also, a deterministic method of seismic inversion was applied to the seismic data to compare its result with the stochastic inversion results. The mean and variance of the elastic properties were calculated to evaluate the uncertainties in the different locations.

4. Results and discussion

As mentioned above, two hundred realizations of the acoustic impedance and shear impedance were generated according to the workflow shown in Figure 1. The crossline of the studied field containing the well with original S-wave velocity log is chosen to present the results of the research.

Figures 2a-2d represent the results of stochastic inversion for acoustic impedance recovery. Figure 2a shows the mean of two hundred generated realizations of the acoustic impedance. Figures 2b and 2c illustrate three single realizations of the acoustic impedance specified by their numbers, namely realizations 35, 100, and 175, respectively. Since the averaging process acts as a low-pass filter, the mean of the realizations is a smoothed map which is similar to the result of deterministic inversion and does not meet the expectations of modeling the heterogeneities of the reservoirs. Unlike the smoothness of the mean of the realizations, every single realization illustrates a high-frequency map of acoustic impedance showing the details in the reservoir because variogram models control the frequency content of such results. The reservoir sands have been separated from the surrounding carbonate layers by low acoustic impedance. Three realizations in Figure 2 are relatively different in high-frequency contents, but they are the same in the background or low-frequency contents. Furthermore, they have the same probability of occurrence within the reservoir. Some thin beds are specified by high acoustic impedance in Figures 2b and 2c; their structures are different in each realization. Some of them are continuous in a realization, while they are discontinuous in another realization. All realizations have a feature in common; they have distinguished the thin beds whether continuous or discontinuous.

Figure 2

Results of stochastic inversion of seismic data for acoustic impedance recovery; a) mean of two hundred realizations of acoustic impedance; b) realization 35 of acoustic impedance; c) realization 100 of acoustic impedance; and d) realization 175 of acoustic impedance; well location is specified by a dark direct line, and the reservoir is between 715 and 790 ms in the well location. Two carbonate layers surrounded the sandstone reservoir are specified by high values of the acoustic impedance; the mean of the acoustic impedance realizations is a smoothed map, but single realizations show more details and are closer to reservoir facts.

Figures 3a-3d illustrate the results of shear impedance recovery from the seismic data. Figure 3a shows the mean of the shear impedance realizations, and Figures 3b-3d represent the numbers 35, 100, and 175 of the shear impedance realizations. As the mean of acoustic impedance, the shear impedance mean contains low-frequency components and is not able to extract the details and heterogeneities of the reservoir, but single realizations are closer to reservoir facts and include high-frequency components in the background of a smoothed map. Shear impedance realizations have the same probability to occur within the reservoir. Mean and single realizations show that shear impedance values are low in the reservoir area and increase at the top and bottom of the reservoir. The mean of all the realizations shows the changes of impedance on a large scale, but single realizations demonstrate both overall and local changes of shear impedance which are capable of recognizing thin beds.

Figure 3

Results of stochastic inversion of seismic data for shear impedance recovery; a) mean of two hundred realizations of shear impedance; b) realization 35 of shear impedance; c) realization 100 of shear impedance; and d) realization 175 of shear impedance; well location is specified by a dark direct line; according to (c) and (d), single realizations are able to model the heterogeneities of the reservoir, but the mean of realizations does not contain the details of the reservoirs.

All the single realizations have the same probability of occurrence in the reservoir. Because of the high resolution of the single realizations they can model the expected heterogeneities in the reservoir. Vertical and horizontal permeabilities are two important factors affecting fluid flow in the reservoir and are directly connected to the heterogeneities of the reservoir. Therefore, it is essential to model spars and thin beds embedded within the main reservoir layer. Due to modeling the heterogeneities of the reservoirs, stochastic seismic inversion is appropriate for calculating oil in place volume, reservoir connectivities, and fluid flow simulation in the reservoir.

All the realizations of the elastic properties reproduce the seismic section with high correlation coefficients to real data. Figure 4 shows the comparison of the stochastic inversion results and measured values in the well location. Four variables, including the acoustic impedance, shear impedance, the ratio of P-wave velocity to S-wave velocity, and the synthetic trace are compared between the well logs and the stochastic inversion results in the reservoir interval. The purple curves represent the results of stochastic inversion (realization 100 for both acoustic impedance and shear impedance), but the green curves show the measured properties by well logging. In the well location, the stochastic inversion has recovered elastic properties efficiently, and the synthetic trace from the well logging is highly correlated to synthetic trace from the stochastic inversion result.

Figure 4

Variance of the elastic properties obtained from two hundred realizations; a) variance of acoustic impedance, and b) variance of shear impedance; uncertainty increases in the locations far from the well and approaches its minimum in the well location.

One of the main advantages of stochastic inversion is its ability to evaluate uncertainty. Uncertainties associated with the results of stochastic inversion are obtained by calculating the variance of the elastic properties in each grid. All of two hundred realizations are used to evaluate the uncertainties within the reservoir. Figure 5 illustrates the uncertainties related to the elastic properties in the whole inverted area. Figure 5a is the variance of the acoustic impedance, and Figure 5b shows the shear impedance variance. By increasing distance from the well, the uncertainty increases in both acoustic impedance and shear impedance.

Figure 5

Analysis of the result of stochastic inversion in the well location and in the reservoir interval; the purple curves are the result of the stochastic inversion (from realization 100), while the green curves are measured well logs; compared curves are the acoustic impedance, shear impedance, the ratio of P-wave velocity to S-wave velocity, and the synthetic trace from the left to right; stochastic inversion has recovered the elastic properties efficiently in the well location.

Figure 6

Specification of the thin beds on the well log data; P-wave and S-wave velocities increase in carbonate thin beds and are highlighted by a green color.

Figure 6 shows some thin carbonate beds on the well log data specified by increasing P-wave and S-wave velocities. Four thin beds are highlighted by a green color, and their recovery is expected from the inversion process. A deterministic method of inversion was applied to seismic data, and its result is compared with the stochastic one for detecting the thin beds within the reservoir. Figure 7 illustrates the results of both stochastic and deterministic seismic inversion and shows their capability in manifesting the subsurface structures. Figures 7a and 7b show the results of stochastic inversion (realizations 35 and 100 of the acoustic impedance), and Figure 7d represents the result of the deterministic inversion. Stochastic inversion results show that there are four thin beds within the reservoir which are specified by high values of acoustic impedance. The recovered thin beds are at 722, 735, 767, and 773 ms and have the thicknesses of 5.7, 4.3, 3.5 and 4.5 m in the well location, respectively. The structure of one thin bed is not the same in both realizations, but the important thing is recovering this thin bed in both realizations. Comparison of these two methods shows the main differences between them; the results of the stochastic inversion demonstrate its capability in the recovery of thin beds, but the smoothness of deterministic inversion result has removed the details and just shows the changes of the acoustic impedance on a large scale. Thin beds are specified by blue ellipses in Figures 7b and 7c.

Figure 7

Comparison of stochastic and deterministic inversion results for detecting thin beds; a) realization 35 of acoustic impedance; b) realization 100 of acoustic impedance; and c) deterministic result of acoustic impedance; well location is specified by a dark direct line; stochastic inversion is able to recover the thin beds, but deterministic inversion shows a large-scale map of acoustic impedance; thin beds have been specified by blue ellipses in both (b) and (c).

5. Conclusions

Stochastic inversion was applied to one of the oilfields in the Persian Gulf, and two hundred realizations of acoustic impedance and shear impedance were generated using a Bayesian framework and an SGS technique. The results of stochastic inversion are high resolution maps of acoustic impedance and provide the possibility to evaluate the uncertainties within the reservoir. Uncertainty assessment shows that increasing the distance from the well causes an increase in the variance of the elastic properties. Moreover, stochastic inversion is capable of recovering the thin beds from 3.5 to 5.7 m, but deterministic inversion result is a low frequency map of elastic properties and does not contain the details of the reservoir. Thus, stochastic inversion is an appropriate method for modeling the heterogeneities of the reservoir and provides better results for reservoir characterization.

Acknowledgments

The authors are thankful to CGG Company for providing the academic license of Hampson-Russell software for Amirkabir University of Technology.

Nomenclatures