Application of different classification methods for litho-fluid facies prediction: a case study from the offshore Nile Delta

Abstract

In this work we test four classification methods for litho-fluid facies identification in a clastic reservoir located in the offshore Nile Delta. The ultimate goal of this study is to find an optimal classification method for the area under examination. The geologic context of the investigated area allows us to consider three different facies in the classification: shales, brine sands and gas sands. The depth at which the reservoir zone is located (2300–2700 m) produces a significant overlap of the P- and S-wave impedances of brine sands and gas sands that makes discrimination between these two litho-fluid classes particularly problematic. The classification is performed on the feature space defined by the elastic properties that are derived from recorded reflection seismic data by means of amplitude versus angle Bayesian inversion. As classification methods we test both deterministic and probabilistic approaches: the quadratic discriminant analysis and the neural network methods belong to the first group, whereas the standard Bayesian approach and the Bayesian approach that includes a 1D Markov chain a priori model to constrain the vertical continuity of litho-fluid facies belong to the second group. The ability of each method to discriminate the different facies is evaluated both on synthetic seismic data (computed on the basis of available borehole information) and on field seismic data. The outcomes of each classification method are compared with the known facies profile derived from well log data and the goodness of the results is quantitatively evaluated using the so-called confusion matrix. The results show that all methods return vertical facies profiles in which the main reservoir zone is correctly identified. However, the consideration of as much prior information as possible in the classification process is the winning choice for deriving a reliable and physically plausible predicted facies profile.

litho-fluid facies identification, classification methods, AVA inversion

1 Introduction

Seismic reservoir characterization exploits recorded seismic reflection data to infer the petrophysical rock properties and the litho-fluid classes around reservoir zones. Facies classification is one of the key modeling components in reservoir characterization, because the distribution of rock properties in the static reservoir model depends on the outcomes of this classification procedure. The classification process is often based on seismic attributes (i.e. P-wave impedance, I_p; S-wave impedance, I_s; density, ρ) estimated by inverting pre-stack seismic data. One of the most used inversion approaches is the amplitude versus angle (AVA) method (Ostrander 1984, Rutherford and Williams 1989, Mazzotti 1990, Mazzotti and Zamboni 2003, Aleardi 2015a, Aleardi and Tognarelli 2016) that derives the elastic properties of the subsurface by exploiting the variations in the reflected amplitudes with increasing incidence angle. Other approaches that are being presently researched include full waveform inversion, which considers the full information content of the seismic data to derive high-resolution subsurface models (Bachrach et al2004, Sears et al2008, Vigh and Starr 2008, Sajeva et al2016, Aleardi and Mazzotti 2017).

For each investigated case, different litho-fluid classes are defined according to specific knowledge about the investigated area, which can be derived from a preliminary geologic interpretation integrated with core or well log data. In the optimal case, these data represent all the possible litho-fluid facies expected in the area. Otherwise, many methods based on rock physics modeling can be used to mimic all the admissible geologic scenarios in the area. For example, the Gassmann fluid substitution can be used to simulate the effect on the elastic properties of different saturation conditions and different fluid types (Avseth et al2005). These considerations make it clear that the available information used to define the characteristic of each facies plays a central role in facies classification because unreliable and erroneous information can lead to erroneous and biased predictions.

Many methods have been developed to classify the elastic attributes in different litho-fluid classes. In this work we test and compare four litho-fluid facies identification methods on a clastic reservoir located in the offshore Nile Delta. The targets of the investigation are gas-bearing sands at the depth range of 2300–2700 m that pertain to stacked turbidite channel systems of Plio-Pleistocene age. The reservoir covers an area of approximately 100 km² and mainly consists of clean-sand layers (with negligible cementation and low clay content) interbedded with laminated, non-permeable shales. Additional information about the geologic setting of the offshore Nile Delta area can be found in Aal et al2000, Samuel et al2003, Roberts et al2005, Vandré et al2007, Cross et al2009, and Sharaf et al2015.

More than 4000 well log samples from four wells drilled in the investigated area and through the reservoir zone were used to separate the different facies in terms of petrophysical properties (water saturation, shaliness and porosity) and elastic properties (seismic impedances and density). The classification was performed on the feature space defined by the P- and S-wave impedances (I_p and I_s, respectively) that were derived from the observed seismic data by means of a Bayesian linearized AVA inversion (Buland and Omre 2003). The analyzed case was particularly challenging due to the significant overlap between the elastic properties of brine and gas sands. This peculiarity can be explained taking into account the deep depth interval where the reservoir zone is located (2300–2700 m). In fact, it is well known (Avseth et al2003) that the elastic properties of sand and shale are more and more similar as the burial depth increases. The depth increase tends to hide the effect of different fluid saturations on the elastic properties and thus makes the discrimination between different saturation conditions more problematic.

The classification methods we consider are well known in seismic reservoir characterization and have been extensively applied worldwide for lithology and fluid predictions (some examples and additional references can be found in Avseth et al2005). These methods can be conveniently divided into two main categories: deterministic and probabilistic approaches. The quadratic discriminant analysis (DA) and the neural network (NN) approach belong to the first group, whereas the standard Bayesian (SB) approach and the Bayesian approach that includes a downward 1D Markov chain a priori model (MC) to constrain the vertical continuity of litho-fluid facies belong to the second group. The two deterministic methods return classifications in which each time sample is classified into one class or into another class without giving any indication about the probability that the sample effectively belongs to that class. Among the deterministic approaches, the quadratic DA (Avseth et al2005) considers that the different classes are divided by quadratic discriminant surfaces in the feature space, while the NN method is able to discriminate classes divided by highly non-linear discriminant surfaces. This difference allows us to check whether the assumption of quadratic discriminant surfaces yields a reliable classification in the investigated area or if more complicated non-linear surfaces are needed.

The two Bayesian classification methods exploit the seismic likelihood function and a set of a priori information (derived from well log data) to produce a posteriori probability that describes the probability that each given sample belongs to a particular litho-fluid class. The first Bayesian method considered is what we call the SB approach in which only the overall proportion of facies in the target interval is given as a priori information (Avseth et al2005). If we consider a 1D vertical profile, this approach classifies each given input sample independently from the adjacent classified samples. In the second Bayesian approach (MC) a 1D Markov chain a priori model (in the form of a transition probability matrix) is given as additional a priori information to constrain the vertical continuity of the litho-fluid facies (Larsen et al2006).

The paper starts with a brief reminder of the Bayesian linear AVA inversion and on the four classification methods we used. Then, we analyze in more detail the available well log data and the elastic properties of each facies. A synthetic inversion based on actual well log data is performed to check the reliability, the drawbacks and the applicability of each classification method in an optimal case with noise-free seismic data. Finally, we discuss the results of facies classification on field data.

2 The methods

In this section we first introduce the linear Bayesian AVA inversion that we apply to derive the P-wave and S-wave impedances, and the density from seismic data. Then, the four classification methods we tested are described. Additional references are given for the interested reader. In the following discussion m is the vector containing the natural logarithm of elastic properties (m = [ln(I_p), ln(I_s), ln(ρ)]^T), d represents the observed seismic data, and f is the vector expressing the three facies we consider (f = [shale, brine sand, gas sand]^T).

2.1 Bayesian linear AVA inversion

Buland and Omre (2003) proposed a linear AVA inversion method in which they assumed log-Gaussian distributed elastic properties (density, seismic velocities or seismic impedances) and Gaussian distributed noise in the seismic data with zero mean and a covariance matrix Σ_e. They also defined a Gaussian a priori model for the natural logarithm of the elastic properties (p(m)) with expectation μ_m and covariance Σ_m. These considerations result in the final a posteriori distribution p(m∣d) being Gaussian, with analytical expressions for the a posteriori expectation μ_m∣d and covariance Σ_m∣d:

$${\mathit{\mu }}_{\mathbf{m}\mid \mathbf{d}}={\mathit{\mu }}_{\mathbf{m}}+{\displaystyle {\sum }_{\mathbf{m}}{\mathbf{G}}^{\mathbf{T}}{\left(\mathbf{G}{\sum }_{\mathbf{m}}{\mathbf{G}}^{\mathbf{T}}+{\sum }_{\mathbf{e}}\right)}^{-1}(\mathbf{d}-\mathbf{G}{\mathit{\mu }}_{\mathbf{m}})}$$

1a

$${\sum }_{\mathbf{m}\mid \mathbf{d}}={\sum }_{\mathbf{m}}+{\displaystyle {\sum }_{\mathbf{m}}{\mathbf{G}}^{\mathbf{T}}{\left(\mathbf{G}{\sum }_{\mathbf{m}}{\mathbf{G}}^T+{\sum }_{\mathbf{e}}\right)}^{-1}\mathbf{G}{\sum }_{\mathbf{m}}}$$

1b

where G is the linear matrix operator that relates the observed seismic data to the natural logarithm of the elastic properties. The a priori assumptions together with equations (1a) and (1b), determine that the a posteriori probability distribution (PPD) can be defined as:

$$p(\mathbf{m}\mid \mathbf{d})=N\left({\mathit{\mu }}_{\mathbf{m}\mid \mathbf{d}},{\sum }_{\mathbf{m}\mid \mathbf{d}}\right)$$

2

where N indicates the normal Gaussian distribution.

2.2 Discriminant analysis (DA)

Linear or quadratic DA is a very common method used in pattern recognition or machine learning to separate a set of data in different classes or categories basing on a set of measured properties or attributes (Davis and Sampson 1986). This method assumes that the input properties of each class follow a Gaussian distribution and use a training dataset on which the statistical characteristics (in terms of mean vector and covariance matrix) of this distribution are defined. In the case of linear DA, the covariance matrices for all the classes are assumed to be identical, whereas quadratic DA considers a different covariance matrix for each class. If the covarance matrices are identical the discriminant surfaces in the feature space are linear, whereas they become quadratic when the covariance matrices are different.

In this work, we use the quadratic DA and we classify each sample on the basis of the Mahalanobis distance, which in our case can be computed as follows:

$${\mathbf{d}}_z^{\mathrm{f}}={({\mathbf{m}}_z-{\mathit{\mu }}_{\mathbf{m}}^{\mathrm{f}})}^T{\left({\sum }_{\mathbf{m}}^{\mathrm{f}}\right)}^{-1}({\mathbf{m}}_z-{\mathit{\mu }}_{\mathbf{m}}^{\mathrm{f}})$$

3

where m is the vector containing the input sample to be classified, that is the estimated elastic properties at a given vertical position z; μ_m and Σ_m are the mean vector and the covariance matrix of the elastic properties in each facies, while the superscript f indicates that these statistical properties are facies dependent. In our case, μ_m and Σ_m are computed from available well log data. After determining the ${\mathbf{d}}_z^{\mathrm{f}}$ values pertaining to each sample, these are classified into the litho-fluid class resulting in the minimum Mahalanobis distance. Application of the Mahalanobis distance to the litho-fluid facies classification can be found in Mukerji et al2001, Avseth and Mukerji 2002, Avseth et al2005 and Del Monte et al2011.

2.3 Neural networks (NN)

NN are a class of algorithms inspired by the central nervous system of animals that are used to estimate a function (a relation) that transforms an input ensemble of data (e.g. measured properties) in a given set of output. The output can be another set of properties or even a set of classes or categories in which the input data has been divided. Once this relation is defined on a training dataset, it can be used to classify or transform any other set of input data. The NN can be simply presented as a system of interconnected neurons exchanging messages between each other. The connections between the different neurons are defined by numeric weights that are iteratively tuned during the learning process of the network. This approach can be very useful in facies classification when the discriminant surfaces in the feature are highly non-linear (Avseth et al2005). Among the many implementations of NN, we use the popular multi-layer feed-forward NN, in which the weights for each neuron are adjusted using the back-propagation method (Van der Baan and Jutten 2000). The architecture of this network consists of one input layer, one output layer, and one or more hidden layers. The nodes (also called neurons) of each layer are connected by weights, and the nodes of the hidden layers can be non-linear or linear activation functions. We consider one hidden layer with ten neurons and a sigmoid transfer function (the function that connects the different neurons). Figure 1 gives a schematic representation of the network used in this work. Unfortunately, there is no unique and fast rule for deciding the number of hidden layers and the number of neurons in each layer, and in practice this choice is made by a trial and error procedure that tries to balance the computational effort, the convergence of the algorithm and the network performance. It is known that NN are prone to the overfitting problem (Hampson et al2001) and in this work the cross-validation method is used to attenuate this drawback (Geisser 1993). In particular, when training multi-layer networks, the general practice is to divide the data input to the network into three subsets: the training, the validation, and the test sets. The training set is used for computing the gradient and updating the network characteristics, while the error on the validation set is monitored during the training process. The validation and the training set errors normally decrease during the initial phase of training. However, when the network begins to overfit the data, the error on the validation set begins to rise. The test set error is not used during training, but it is used to compare different models. It is also useful to plot the test set error during the training process. If the error in the test set reaches a minimum at a significantly different iteration number than the validation set error, this might indicate a poor division of the dataset. In our work, 15% of the entire ensemble of available well log samples constitutes the validation set, 15% is used as the test set, while the remaining 70% is used as the training set. The final network will be the one determined at the minimum of the validation set error. Examples of NN applications to solve geophysical problems can be found in Avseth and Mukerji 2002, Aleardi 2015b, and Aleardi and Ciabarri 2017.

Figure 1.

Schematic representation of the network used in this work.

2.4 Standard Bayesian classification (SB)

The classification problem can be conveniently cast in a Bayesian framework (Tarantola 2005). This approach, different from the classical DA and NN approach, returns not only a deterministic solution (that is the facies predicted for each input sample), but also the probability that a given input sample belongs to a particular class. Assuming that the properties in each class are Gaussian distributed, the conditional probability of class f given a vector m is:

$$p(\mathbf{f}\mid \mathbf{m})={\displaystyle \frac{p(\mathbf{m}\mid \mathbf{f})p(\mathbf{f})}{p(\mathbf{m})}}$$

4

where $p(\mathbf{f}\mid \mathbf{m})$ is the likelihood function, p(f) is the a priori distribution of the overall proportions of the facies that can be provided by geological analysis of facies profiles derived from wells in the investigated zone, and p(m) is the unconditional probability of the properties values across all the facies. The p(m) distribution in practical applications can be simply viewed as a normalizing constant. The likelihood function can be computed as follows:

$$p(\mathbf{m}\mid \mathrm{f})\propto \exp \left(-{\displaystyle \frac{1}{2}}{(\mathbf{m}-{\mathit{\mu }}_{\mathbf{m}}^{\mathrm{f}})}^T{\left({\sum }_{\mathbf{m}}^{\mathrm{f}}\right)}^{-1}(\mathbf{m}-{\mathit{\mu }}_{\mathbf{m}}^{\mathrm{f}})\right).$$

5

In equation (5), the term in the brackets is similar to the Mahalanobis distance shown in equation (3). However, note that the Bayesian classification as given in equation (4) additionally considers the a priori probability of observing a particular class. Therefore, if a given input sample has the same likelihood of belonging to different classes, it is assigned to the most probable class, that is the class with the highest p(f). Application of this method to facies prediction can be found in Avseth and Mukerji (2002) and Avseth et al (2005).

2.5 Bayesian classification with 1D Markov chain a priori model (MC)

This classification method was proposed by Larsen et al (2006), to which we refer the reader for a detailed description. They developed a 1D facies classification method in a Bayesian framework in which a stationary downward MC is used to constrain the continuity of litho-fluid facies along a vertical profile. The a priori vertical continuity of facies is defined by the transition matrix that gives the probability of having a facies x at vertical position z after observing a facies y at position z-1. This matrix can be determined from the available well log data and from geological knowledge about the investigated area. Different from the previously described SB approach, in which each sample along the vertical profile is classified independently from the classified samples above or below it, this method classifies each sample taking into account the classifications of the adjacent samples. The introduction of a vertically coupled a priori Markov model, or in other words, the introduction of information about the vertical continuity of litho-fluid classes, constitutes additional a priori information that has the effect of stabilizing the classification process. This method can be implemented by means of a downward/upward recursive algorithm that was described in detail by Larsen et al (2006). The target a posteriori probability expressing the probability of a facies f given the estimated elastic properties m at a given vertical position z can be expressed as:

$$p({\mathbf{f}}_z\mid {\mathbf{m}}_z)={\displaystyle \prod _{z}p({\mathbf{f}}_z\mid {\mathbf{f}}_{z-1})}{\displaystyle \int p({\mathbf{m}}_z\mid {\mathbf{f}}_z)p({\mathbf{m}}_z\mid \mathbf{d})\mathrm{d}\mathbf{m}}$$

6

where the probability $p({\mathbf{f}}_z\mid {\mathbf{f}}_{z-1})$ can be obtained from the downward Markov chain transition matrix of litho-fluid classes estimated on actual well log classification.

3 Log data analysis and distribution of elastic properties in each facies

Given the geological setting of the investigated area, three litho-fluid facies were defined: shale, brine sand and gas sand. Figure 2 displays the logged elastic properties (figures 2(a)–(c)) and the interpreted petrophysical properties (figures 2(d)–(f)) of a single well that reached the reservoir interval. This figure shows that the reservoir consists of a sequence of three/four main sand layers interbedded with laminated non-permeable shales. Figure 2(g) shows the corresponding vertical facies profile derived on the basis of the interpreted petrophysical properties and of the chosen cut-off values:

• shale: shaliness >50%, porosity <10%;

• brine sand: shaliness <50%, porosity >10%;

• gas sand: shaliness <50%, porosity >10%, water saturation <50%.

Figure 2.

Example of well log data around the reservoir interval. (a)–(f) refer to V_p, V_s, density, water saturation, porosity and shaliness. In (g) the vertical facies profile derived from the petrophysical curves of porosity, shaliness and water saturation is shown. Black, yellow and red code the shale, brine sand and gas sand, respectively.

Tests carried out with small variations in the cut-off values gave negligible differences in the results.

Cross-plots showing the influence of the petrophysical properties of water saturation, porosity and shaliness on the natural logarithm of P- and S-wave impedances, are represented in figure 3, whereas figure 4 contains cross-plots that clearly display the facies dependency of the I_p and I_s values. The three facies we considered represent the major litho-fluid classes found in the investigated levels. Although defining more classes might be theoretically possible, the inability to recognize litho-fluid classes in finer detail (e.g. discriminating shaly sand or low porosity sand) from their elastic properties has discouraged a finer classification. In figure 4 note that I_p progressively decreases passing from shale, to brine sand and to gas sand. The lower I_p value for the gas sand with respect to the brine sand, is related to the decrease in bulk modulus that occurs when gas replaces brine in the pore space and to the lower density value of gas with respect to brine. Conversely, passing from shale to sand I_s increases due to the lower shear modulus that characterizes the shale (figure 4(a)). If we limit our attention to the sands, the S-wave impedance remains more or less constant. This fact is related to the null effect of the saturating fluid on the shear modulus and to the density decrease that occurs when brine is replaced by gas. These two contrasting processes (increasing V_s and density decrease) determine that the net I_s value remains constant when passing from total brine to gas saturation. Figure 4(b) shows the distribution of elastic properties for the three facies projected onto the ln(I_p)–ln(density) plane. It is clear that the density shows a considerable decrease moving from shale, to brine sand, and to gas sand. The different I_p, I_s and density values of the shale with respect to the sand can be obviously ascribed to the different mineralogical composition and geometric texture, whereas the significant overlap that characterizes the distribution of I_p and I_s for gas sand and brine sand can be ascribed to the depth range where the reservoir is located. In fact, it is well known (Avseth et al2003) that the increasing of burial depth tends to mask the fluid effect, and makes the discrimination between different saturation conditions problematic. For this reason, the case under examination can be considered a challenging test for any classification method that aims to identify the litho-fluid facies from the elastic property values.

Figure 3.

(a)–(c) Rock physics templates derived from actual well log data and showing the I_p and I_s changes with water saturation, porosity and shaliness, respectively.

Figure 4.

(a) and (b) Cross-plots showing the distribution of the natural logarithm of P-wave, S-wave impedances and density for each facies. The data were extracted from the available borehole data around the target interval and were classified basing on the petrophysical properties derived information evaluation analysis.

The DA and the two Bayesian classification methods we consider assume that the properties used to classify each input sample follow a Gaussian distribution in each class. To verify whether this basic assumption is met in the case under examination, we used the so-called Gaussian probability plot. This plot looks close to a straight line if the data are approximately normally distributed. Deviations from a straight line suggest departures from normality. In figure 5 we show the Gaussian probability plots derived on the basis of well log data extracted around the target interval. We observe that, despite some minor curvature particularly visible for the shale, the Gaussian assumption can be considered acceptable for each facies.

Figure 5.

Gaussian probability plot for the ln(I_p), ln(I_s) and ln(density) for each facies. (a)–(c) Shale, brine sand and gas sand, respectively. The red lines represent the expected alignment in the case of a theoretical Gaussian distribution, whereas the blue circles represent the distribution of the natural logarithm of logged properties. Note that the distribution of the actual data is very close to the Gaussian distribution. Minor deviations from the normal distribution characterize the elastic properties of the shale.

Based on the results shown in figures 4 and 5, we can approximate the overall distribution of the natural logarithm of the elastic properties in the target area to a Gaussian mixture distribution (Gallop 2006). This is a linear combination of different normal distributions that in our case can be written as:

$$\begin{array}{ccc}p(\ln (I_{\mathrm{p}},I_{\mathrm{s}},\rho ))\hfill & =\hfill & {\displaystyle \sum _{k=1}^n{\lambda }_k}{p}_k(\ln (I_{\mathrm{p}},I_{\mathrm{s}},\rho ))\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{k=1}^n{\lambda }_k}N\left({\mathit{\mu }}_{\mathbf{m}}^{\mathbf{k}},{\sum }_{\mathbf{m}}^{\mathbf{k}}\right)\hfill \end{array}$$

7

where n indicates the number of components of the mixture that in our application coincides with the number of litho-fluid facies we consider; N identifies the Gaussian distribution, λ_k represents the weights of the components (i.e., the proportion of each the facies), while μ_m and Σ_m are the mean vector and the covariance matrix. To derive the statistical properties of the Gaussian mixture model expressed in equation (7) we use the expectation maximization algorithm (Hastie et al2002). The resulting distribution is shown in figure 6, together with the 1D marginal distributions projected along the ln(I_p), ln(I_s) and ln(ρ) axes. Again, the significant overlap between the I_p and I_s values associated with the brine sand and gas sand is clearly visible.

Figure 6.

The Gaussian mixture distribution for the natural logarithm of the elastic properties around the target interval. (a) 2D joint distribution of ln(I_p) and ln(I_s) and the corresponding 1D marginal distributions projected along the ln(I_p) and ln(I_s) axes. (b) 2D joint distribution and the corresponding projections along the ln(I_p) and ln(ρ) axes. Note the significant overlap between the I_p and I_s values associated with the brine sand and gas sand.

4 Inversion and classification of synthetic seismic data

To check the practical applicability of the four analyzed methods for litho-fluid facies classification in the target area, we show the classification results obtained on the elastic properties estimated by inverting synthetic seismic data. The synthetic data were computed on depth models derived from the well log recordings already shown in figure 2, and by means of 1D convolutional forward modeling. A 50 Hz Ricker wavelet was considered as the source signature. Figure 7(a) shows the synthetic CMP gather in which the offset has been converted to incidence angles, whereas figures 7(b)–(d) illustrate the results of the Bayesian AVA inversion. In figure 7(a), the strong reflection at 2.46 s indicates the target gas sand. The decreases in I_p and density and the increases in I_s that characterize the gas sand interval generate the typical class III AVA anomaly (Castagna and Swan 1997) clearly visible in the synthetic seismogram. In figures 7(b)–(d) the blue curves depict the actual well log data resampled at the seismic sample interval, the red curves illustrate the maximum a posteriori (MAP) solution, and the gray curves are Monte Carlo realizations drawn from the a posteriori distribution p(m∣d) (see equation (2)). Each Monte Carlo realization represents a possible solution of the Bayesian AVA inversion. As expected, the uncertainties increase passing from I_p, to I_s and to density. However, note that, although with a different resolution, the predicted elastic properties closely match the true ones. Figure 7(e) shows the true vertical facies profile determined from the logged petrophysical property values.

Figure 7.

Results of the synthetic AVA inversion based on actual well log measurements. The synthetic CMP, P-impedance, S-impedance and density are represented in (a) to (d), respectively. For impedances and density, the x-axes are drawn with the same scale to better illustrate the increase in uncertainty that occurs passing from the I_p, to the I_s and to the density estimates. The blue curves represent the true elastic properties, the red the MAP solution, whereas the gray curves show Monte Carlo samples drawn from the PPD. (e) The true vertical facies profile derived from the actual well log measurements. The black, yellow and red colors identify shale, brine sand and gas sand, respectively. Note the clear class III AVA anomaly that corresponds to the gas sand interval.

As previously discussed, the Bayesian AVA inversion performed according to Buland and Omre (2003) returns the Gaussian probability a posteriori distribution of the natural logarithm of the elastic properties (I_p, I_s and density). However, we do not consider the density estimates in the classification process because the linear AVA inversion cannot retrieve reliable information on this parameter with realistic noise levels. Moreover, before estimating the elastic properties of each facies from well log data, the borehole data are filtered by applying the Backus averaging method (Backus 1962). This allows us to take into consideration in the classification procedure the different vertical resolutions of borehole measurements and of the elastic properties estimated from seismic data.

In the following tests, to qualitatively evaluate the ambiguity related to the different classification algorithms, we use as inputs for the classification methods different Monte Carlo realizations drawn from the PPD p(m∣d). In addition, for the two Bayesian classification approaches we also show the PPD representing the probability of occurrence of each facies at each given vertical position. The final results are represented in figure 8. Comparing the outcomes of the DA and the NN approaches we observe that they return very similar results, thus demonstrating that in this particular case the discriminant surfaces in the feature space (the 2D space containing the natural logarithm of I_p and I_s) can be conveniently approximated to quadratic surfaces. It can be noted that the gas sand layer and most of the brine sand layers are correctly classified, although some shale layers are erroneously classified as brine sand or gas sand. These misclassifications can be ascribed to the overlap between the elastic properties of each facies. Also, note that the vertical facies profiles predicted by the DA and NN methods show some unphysical characteristics, with a brine sand layer placed just above a gas sand interval. This is obviously an unphysical situation, with the gas less dense than brine. If we analyze the SB method, which considers the additional information on the a priori occurrence of each facies, we observe that the predicted vertical facies profile is slightly closer to the true one, but some unphysical vertical successions of facies are still present. These unphysical characteristics are not present in the Bayesian classification with a 1D Markov chain prior model in which we have imposed a null probability for the transition from a brine sand to a gas sand layer (table 1).

Figure 8.

Results of facies classification obtained with the results of the linear AVA inversion of synthetic seismic data. (a)–(d) The results for the DA, NN, SB and MC, respectively. In each part, the true vertical facies profile determined from well log data is shown on the left. The other profiles represent the results obtained by considering as input to the classification process five Monte Carlo realizations drawn from the a posteriori distribution p(m∣d), and the MAP solution of the AVA inversion. For the Bayesian classifications, the PPD pertaining to the MAP solution are also shown.

Table 1.

The transition probability matrix derived from well log data for the considered case. See the text for additional explanations.

	Shale	Brine sand	Gas sand
Shale	0.9587	0.0393	0.0020
Brine sand	0.0685	0.9315	0
Gas sand	0.0050	0.0500	0.9450

The stationary transition probability matrix derived for the analyzed case is shown in table 1. Rows correspond to shale, brine sand, and gas sand at generic vertical position i, whereas columns refer to shale, brine sand, and gas sand at vertical position i-1. Note that the downward transition from brine sand to gas sand is impossible, being associated with a null transition probability. This additional a priori information yields a predicted facies profile that is very close to the true one. In this case (figure 8(d)) the majority of sand and shale intervals are correctly predicted, thus demonstrating that when the elastic properties of each facies overlap each other, it is crucial to consider as much a priori information as possible to better constrain the final result.

The quantitative assessment of the classification results can be done by using contingency analysis tools. In particular, in figure 9 the so-called ‘confusion matrices’ are shown, which are computed on the facies profiles predicted when the MAP solution of the AVA inversion is the input of the classification process. In a confusion matrix, each column represents the instances in a predicted class, while each row represents the instances in an actual class. Along the main diagonal, we observe the percentage of samples belonging to a litho-fluid class that have been correctly classified in that class, whereas the off-diagonal terms indicate the percentage of samples that have been misclassified. The confusion matrices for the DA and the NN approaches are very similar, with about 90% of shale, 80% of brine sands and 65% of gas sands correctly predicted. As expected, the overlap between I_p and I_s values produces a significant number of actual gas sand samples that were erroneously classified as brine sand (35%, approximately). The confusion matrix resulting from the SB approach is similar to those derived from the NN and DA. The great advantage in this case is offered by the a posteriori probability, which indicates the probability of the occurrence of each facies at each time sample. Finally, the MC method produces a confusion matrix with 100% of gas sand correctly classified and negligible misclassification errors. Figure 9 further confirms that any a priori information and/or additional constraints, assume a crucial role in litho-fluid facies prediction in the investigated zone.

Figure 9.

(a)–(d) The confusion matrices pertaining to the DA, the NN, the SB, and the MC.

5 Inversion and classification of field seismic data

The acquired seismic data have been processed in an AVA friendly mode (Mazzotti and Ravagnan 1995), that is particular attention was paid to preserve the true amplitude of the reflections. However, the strong attenuation of high frequencies produced by several gas clouds occurring in the shallow layers, determines a very low dominant frequency (around 15 Hz) at the depth of the target level. This characteristic, in addition to the overlap between the elastic properties of the different facies, makes the classification procedure very challenging for any classification method. The discussion about the field data inversion starts by describing the results on a single CMP gather where well control is available to validate the results. This CMP is the nearest to the well that has already been considered in the previous synthetic inversion. Figure 10 illustrates the results of the Bayesian AVA inversion for the considered CMP gather (figure 10(a)). Despite the very low resolution, a clear class III AVA anomaly is visible at the target level (around 2.46 s). In figures 10(b)–(d) the true elastic properties resampled at the seismic sampling interval are shown by the blue curves, whereas the green curves show the true properties up-scaled to the seismic frequency band. The red curves represent the MAP solutions, while the gray lines are Monte Carlo realizations derived from the a posteriori distribution. As observed in the previous synthetic example, the uncertainties increase passing from the impedances to the density estimates. However, the true up-scaled elastic properties (green curves) show a fair match with the estimated ones (red curves) and, more importantly, they lie inside the range defined by the Monte Carlo realizations.

Figure 10.

Example of AVA inversion of field data performed on the CMP location nearest to the well. The synthetic CMP, P-impedance, S-impedance and density are represented from (a) to (d), respectively. For impedances and density, the x-axes are drawn with the same scale to better illustrate the increase of uncertainty that occurs passing from the I_p, to the I_s and to the density estimates. In blue are represented the true elastic properties resampled at the seismic sampling interval, in red the MAP solution, in green the true properties up-scaled to the seismic bandwidth, whereas the gray curves show Monte Carlo realizations drawn from the PPD.

In figure 11 we show the facies profiles predicted when the results of AVA inversion shown in figure 10 are classified. Similar to figure 8, the 1D true facies profile determined from well log data is compared to the 1D vertical facies profiles obtained by the four classification methods. We observe that all methods correctly identify the target gas sand around 2.46 s. In addition, all the methods identify both the sand interval just below the gas layer and the thick sand interval between 2.62 and 2.64 s but the NN method interprets this layer as gas saturated instead of brine saturated. Moreover, all the four methods tend to overestimate the thickness of the layers and are not able to correctly identify the finely layered shale–brine sand sequence. This low resolution affecting the classified facies profiles is obviously imputable to the very low resolution of the seismic data. The Mahalanobis distance strongly overestimates the occurrence of brine saturated layers, whereas the NN method is prone to overestimate the gas sand intervals. Once again, the a priori information inserted into the Bayesian classifications make the probabilistic approaches more reliable than the DA and NN methods. In particular, the a priori information on the occurrence of each facies around the target interval makes the proportions of facies predicted by the SB very close to the true one, although a brine sand interval is erroneously positioned at the top of the gas sand interval. The further inclusion of the stationary vertical transition probability matrix as a priori information, makes the classification given by the MC method very close to the true vertical facies profile with a physically plausible transition between the different facies.

Figure 11.

Results of facies classification on the CMP gather extracted from the 3D seismic cube nearest to the well previously considered in the synthetic test. (a)–(d) The results for the DA, NN, SB and MC, respectively. In each part, the true profile determined from the actual petrophysical properties is shown on the left. The other profiles represent the classification obtained on five Monte Carlo realizations drawn from the p(m∣d) distribution, and on the MAP solution of AVA inversion. For the Bayesian classifications, the PPD pertaining to the MAP solution is also shown.

These considerations are summarized in the information brought by the confusion matrices (figure 12). Figure 12(a) clearly shows that the DA overestimates the occurrence of brine sand layers and underestimates the gas-saturated intervals. In particular, this method tends to interpret 40% of actual shale intervals as brine sand layers, whereas 25% of gas intervals are misclassified as brine sands. The NN (figure 12(b)) tends to interpret 50% of brine saturated sands as gas saturated, whereas 20% of shales are interpreted as brine sands. The misclassifications errors are lower for the Bayesian methods and particularly for the MC approach, which returns a confusion matrix with very low off-diagonal terms values. In particular, this approach produces a confusion matrix with 100% of gas sand correctly classified and negligible misclassification errors that can be primarily attributed to the different resolution between well log information and the seismic data and to noise contaminating the seismic data.

Figure 12.

(a)–(d) The confusion matrices pertaining to the DA, NN, SB, and MC, respectively.

We now describe the elastic properties estimated by the Bayesian AVA inversion along an in-line section extracted from the 3D seismic volume. For simplicity, in the 3D inversion each CMP gather is inverted separately, thus overlooking the spatial correlation of elastic attributes. However, the vertical correlation between elastic properties is taken into account in the inversion by the a priori covariance matrix of elastic properties (Buland and Omre 2003). We point out that the lateral continuity of our results is mainly related to the lateral correlation of seismic data that is imposed by the migration operator. In fact, this operator can be considered a spatial filter with correlation length associated with the Fresnel zone. The lateral continuity of the results is clearly shown in figure 13 in which the white dashed lines identify the spatial location of the well already considered. Note that the lateral extensions of the decreases in I_p and density and increases in I_s that are intersected by the well at 2.46 s. This time interval is associated with the gas-saturated reservoir zone.

Figure 13.

MAP solution of AVA inversion represented along on an in-line section extracted from the 3D seismic cube. (a)–(c) I_p, I_s and density estimates, respectively. The white dashed lines identify the well position. Note the decrease in I_p and density and the increase in I_s occurring at 2.46 s that correspond to the target gas sand interval.

Figure 14 shows the facies predicted along the in-line section of figure 13. These sections show for each CMP position the 1D vertical facies profile estimated from the MAP solution of the AVA inversion. Once again, note the overestimation of brine sand intervals and gas sand intervals produced by the quadratic DA and the NN approach, respectively. The importance of including a priori information clearly stands out when comparing the facies sections of figures 14(a) and (b) with the results of the Bayesian classifications of figures 14(c) and (d). In detail, the poorer performance of the NN and SB methods in the field test with respect to the synthetic example can be mainly ascribed to the narrow bandwidth of the available seismic data that yields predicted elastic properties with low vertical resolution. This characteristic, together with the overlap between the elastic characteristics of the different facies, makes even more crucial the inclusion of a priori information, or constraints in the classification of field data, which cannot be done by NN or SB. The poor performances of the NN method can be also attributed to the well-known overfitting problem, which often results in a sub-optimal prediction capability of the network in the case of noisy data. Another drawback of the NN approach, not discussed here, is related to its local nature. Indeed, in an NN optimization, the weights associated with each neuron are usually randomly initialized and are adjusted using a gradient-based strategy. This demonstrates the importance of a good initial model to prevent convergence towards a local minimum of the error function. To attenuate this limitation, the NN optimization should be repeated many times, starting from different initial models, and finally selecting the best result. Obviously, this procedure is very time consuming in the case of large networks or a large number of input samples.

Figure 14.

Results of facies classification obtained for the in-line section shown in figure 13. (a)–(d) The results for the DA, NN, SB and MC, respectively. For each CMP position the results obtained by classifying the MAP solution of the AVA inversion are shown. The white dashed lines indicate the well position. Black, yellow and red colors refer to shale, brine sand and gas sand, respectively.

Figure 14 shows that the two Bayesian methods identify a possible gas sand interval between cross-line 8420 and 8440 and located around 2.5 s. In this position decreases in I_p and density together with strong increases in I_s were detected by the AVA inversion (see figure 13). Note again that the MC approach ensures a result with a physically plausible vertical succession of litho-fluid facies.

6 Conclusions

We tested four different approaches to classify the P-wave and S-wave impedances estimated from pre-stack seismic data by means of a Bayesian AVA inversion. The methods were applied for facies identification in a clastic, gas-saturated reservoir located in the offshore Nile Delta. The four classification approaches are the quadratic discriminant analysis (DA), the neural network (NN), the standard Bayesian classification (SB) and the Bayesian classification with a 1D Markov chain prior model (MC). The first two methods yield a deterministic classification, whereas the Bayesian approaches cast the classification into a probabilistic framework. In particular, the classification procedure in the two Bayesian approaches is driven and constrained by a set of a priori information that can be derived from available borehole information.

We observed that the deep depth interval in which the reservoir zone is located (2300–2700 m), produces a significant overlap between the elastic properties of each facies. This significant depth interval tends to partially hide the so-called fluid effect, thus making the discrimination between different saturation conditions particularly challenging. However, the example of synthetic data demonstrated that all the methods are potentially able to correctly identify the reservoir zone. In addition, the similar results produced by the DA and NN approaches demonstrated that, in the analyzed case, the discriminant surfaces that separate the different litho-fluid classes in the feature space can be conveniently approximated to quadratic surfaces.

The field data classification was more problematic due to the very low resolution of the seismic data. In this field data test, the quadratic DA and the NN approach are not able to correctly discriminate the fluid effect on the elastic properties related to the different saturations conditions. The quadratic DA tends to overestimate the occurrence of brine saturated sands, while the NN method tends to overestimate the gas sand intervals. In particular, important limitations of the NN method are its local nature, the overfitting problem and the difficulty in designing the optimal NN architecture. For the analyzed case, the consideration of any kind of a priori information in the classification process was revealed to be the crucial aspect to deriving a reliable facies classification. In particular, the SB approach, which simply takes into account the a priori information about the overall proportions of facies around the investigated interval, produces final results in which the actual proportion of facies is closely respected. If a downward 1D Markov chain a prior model is incorporated into the Bayesian classification procedure, the predicted facies profiles are further improved and show a physically plausible vertical succession of litho-fluid classes.

Acknowledgments

This study is a part of a collaborative research project between EDISON RD&I and the Earth Science Department of the University of Pisa. The authors wish to thank EDISON E&P for making the well log data available and for the permission to publish this work and Alfredo Mazzotti for useful discussions and insightful comments.