Method for detecting geological objects in a seismic image

Abstract

The invention is a method applicable to oil and gas exploration and exploitation for automatically detecting geological objects belonging to a given type of geological object in a seismic image, on a basis of a priori probabilities of belonging to a type of geological object assigned to each of samples of the image to be interpreted. The image is transformed into seismic attributes applied beforehand, followed by a classification method. For each of the classes, an a posteriori probability of belonging to a type of geological object is determined for each of the samples of the class according to the a priori probabilities, of the class, of belonging, and according to a parameter α describing a confidence in the a priori probabilities of belonging. Based on the class of the sample, the determined a posteriori probability of belonging to a type of geological object is assigned for the samples of the class. The geological objects belonging to the type of geological object are detected based on determined of the a posteriori probabilities of belonging to the type of geological object for each of the samples of the image to be interpreted.

Claims

1. A computer-implemented method for automatically detecting at least one geological object belonging to a type of geological object in at least one seismic image of a subterranean formation, on a basis of a priori probabilities of belonging to the type of geological object that are assigned to samples of the at least one image, comprising: A. defining an operation for transforming the at least one seismic image into seismic attributes and applying the operation for transforming to determine seismic attributes of the at least one seismic image; B. applying an unsupervised classification operation to the seismic attributes of the at least one seismic image and determining a class of each sample of the at least one seismic image; C. determining for each class an a posteriori probability of each sample of a class belonging to the type of geological object according to the priori probabilities of each sample of the image of the class belonging to the class according to a parameter describing confidence in the a priori probabilities; and D. for each sample of the image, based on the class of the sample, assigning the determined a posteriori probability of belonging to the type of geological object for samples of the class; and wherein the at least one geological object belonging to the type of geological object which is detected based on the determined a posteriori probabilities of belonging to the type of geological object for each sample of the image and a geological model representative of the formation based on at least the detected geological objects, defining a plan for exploiting the hydrocarbons of the formation based on at least the geological model, and exploiting hydrocarbons of the formation according to the defined plan.

2. The method according to claim 1, wherein the a posteriori probability y.sub.i of the samples of the class i belonging to the type of geological object is determined by: $i.\mathrm{if}{p}_{1,i}>\frac{1}{\alpha +1},y_i=s_1.\mathrm{ii}.\mathrm{if}{p}_{1,i}\le \frac{1}{\alpha +1},q_i=\max \left\{k\in 1\text{:}n\mathrm{such}\mathrm{that}\underset{j\in 1:k}{.Math.}{p}_{j,i}\le \frac{1}{\alpha +1}\right\}\mathrm{is}\mathrm{determined},\mathrm{and}$ $a.\mathrm{if}{.Math.}_{j\in 1:q}{p}_{j,i}=\frac{1}{\alpha +1},\mathrm{then}{y}_i=s_q.b.\mathrm{if}{.Math.}_{j\in 1:q}{p}_{j,i}<\frac{1}{\alpha +1},\mathrm{then}$ $\mathrm{if}q=n,y_i=s_n$ $\mathrm{if}q<n,y_i=s_{q+1}.$ where X is the vector containing a value of a priori probability of belonging for each sample of the seismic image; L is the vector containing, for each sample of the seismic image, an identifier of the class to which the sample belongs; S={s.sub.1, s.sub.2, . . . , s.sub.n}={s.sub.k, k∈1:n} with 0=s.sub.1<s.sub.2< . . . <s.sub.n=1, is the set of values taken by elements of the vector X; C.sub.i, i∈L is the identifier class i; and p.sub.k,i=p(s.sub.k|C.sub.i) is the proportion of points of C.sub.i having an a priori probability value of s.sub.k.

3. The method according to claim 1, wherein the seismic attributes are textural attributes.

4. The method according to claim 1, wherein a decrease in dimensions is applied to the determined attributes by the transforming operation.

5. The method according to claim 4, wherein the decrease in dimensions is carried out by a principal component analysis.

6. The method according to claim 4, wherein the decrease in dimensions is carried out by an analysis using an attribute selection method.

7. The method according to step B) of claim 1, comprising applying an unsupervised classification method which is a generative topographic mapping method.

8. The method according to claim 1, wherein at least steps A) to D) are applied to a first seismic image, and wherein at least the steps A) to D) are applied to a second seismic image comprising: 1. transforming the second seismic image into seismic attributes; 2. applying the unsupervised classification operation to the attributes of the second seismic image and determining a class for each sample of the second seismic image; and 3. for each class of the second seismic image and for each sample of the second seismic image belonging to the class, assigning a determined probability of belonging to the type of geological object in the class in the first seismic image.

9. The method according to claim 1, wherein the geological objects belonging to the type of geological object are detected in the at least one seismic image by using a thresholding method applied to a value of the determined a posteriori probabilities of belonging to the type of geological object for each sample of the at least one seismic image.

10. A method for exploiting a subterranean formation comprising hydrocarbons by use of at least one seismic image relating to the formation, comprising: detecting geological objects belonging to at least one type of geological object in the at least one seismic image by use of the method according to claim 1; defining a geological model that is representative of the formation based on at least the detected at least one geological object and defining a plan for exploiting the hydrocarbons of the formation that is based on at least of the geological model; and exploiting the hydrocarbons of the formation according to the plan for exploiting.

11. A computer program product stored on a tangible recording medium is executed by a processor to perform the method of claim 1.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 shows, in a theoretical example, the distribution of the classes determined by use of a classification method, with respect to a geological object to be detected.

(2) FIGS. 2a and 2b show the samples of the class C1 having an a priori probability of belonging to the geological object to be detected of 0 and 0.75, respectively.

(3) FIG. 2c schematically shows the operation of determining the a posteriori probability associated with the class C1.

(4) FIG. 3a shows the distribution of the samples belonging to the class C2, with respect to the geological object to be detected.

(5) FIG. 3b schematically shows the operation of determining the a posteriori probability associated with the class C2.

(6) FIG. 4a shows the distribution of the samples belonging to the classes C4, C5 and C6, with respect to the geological object to be detected.

(7) FIG. 4b schematically shows the operation of determining the a posteriori probability associated with the classes C4, C5 and C6.

(8) FIG. 5a shows the distribution of the samples belonging to the class C3, with respect to the geological object to be detected.

(9) FIG. 5b schematically shows the operation of determining the a posteriori probability associated with the class C3.

(10) FIGS. 6 and 7 show, respectively, three seismic learning sections selected in an overall seismic volume and the a priori probabilities of belonging to an MTD-type geological object for each of the samples of these sections.

(11) FIG. 8 shows the distribution of the samples of the three learning images from FIG. 6 over 49 classes.

(12) FIG. 9 shows the a posteriori probabilities of belonging to an MTD-type geological object for each of the three learning images from FIG. 6, these a posteriori probabilities having been determined according to the method of the invention.

(13) FIG. 10 shows the result of propagating the a posteriori probabilities of belonging to an MTD-type geological object, determined for each of the three learning images from FIG. 6, to 72 sections in the overall seismic volume to be interpreted.

(14) FIG. 11 shows the samples of the overall seismic volume to be interpreted having been detected as belonging to a PTD-type geological object by means of the method according to the invention.

DETAILED DESCRIPTION OF THE INVENTION

(15) In general, one of the subjects of the invention relates to a computer-implemented method for automatically detecting samples of at least one seismic image of a subterranean formation which belong to a predefined type of geological object. The method according to the invention may advantageously be applied to the detection of a type of geological object that is characterized by heterogeneous seismic properties, or in other words a type of geological object that is characterized by heterogeneous seismic facies. Such heterogeneous seismic facies comprise, nonlimitingly, mass-transport deposits (MTDs), fluvial channels with heterogeneous seismic properties or else geological bevels (such as bevelled contact between tilted geological layers and an erosion surface). The method may however also be applied to the detection of geological objects that are known to have homogeneous seismic properties, such as saline bodies, in particular when the internal facies of these objects vary slightly.

(16) Thus, the method according to the invention is implemented for a given type of geological object, and is advantageously repeated if it is desired to detect another type of geological object in the same seismic image.

(17) The seismic image is either a two-dimensional (2D) or three-dimensional (3D) image, in which one of the axes is the time axis (seismic image in the time domain) or else the depth axis (seismic image in the depth domain), and the one or more other axes correspond to the horizontal axes of a geographical space. The method according to the invention may also be implemented on a plurality of 2D and/or 3D seismic images. These one or more seismic images might have been obtained by use of seismic measurements (generally referred to as “seismic acquisition surveys”) and by processing these seismic measurements (generally referred to as “seismic processing”). The steps to be implemented so as to obtain a seismic image of a subterranean formation are known.

(18) The method according to the invention requires the availability of a priori probabilities of belonging to the sought-after type of geological object with these probabilities being assigned to each of the samples of the seismic image under consideration. In other words, each sample (a pixel in the case of a two-dimensional seismic image or a voxel in the case of a three-dimensional seismic image) of the seismic image under consideration must have a value that is representative of the probability of this sample belonging to a geological object of the sought-after type.

(19) According to one nonlimiting embodiment of the invention, the method may comprise the following steps:

(20) 1. Determining a priori probabilities in the seismic image

(21) 2. Determining seismic attributes in the seismic image

(22) 3. Classifying the seismic image in the attribute space

(23) 4. Determining an a posteriori probability of belonging

(24) 5. Propagating the a posteriori probabilities to another seismic image

(25) 6. Automatically detecting geological objects

(26) 7. Exploiting the subterranean formation

(27) However, in general, the method according to the invention comprises at least steps 2 to 4 and step 6.

(28) Thus, the obtaining a priori probabilities of belonging, required for implementing the invention, may be obtained by implementing step 1, which is optional.

(29) According to one advantageous implementation of the invention, when there are large volumes of seismic data to be interpreted, steps 2 to 4 are applied to at least one seismic image, which is a seismic learning image, and the parameters of the method that are determined using this seismic learning image are subsequently applied, as described in step 5, to at least one second seismic image, for example of larger dimension than the seismic learning image. Thus, in the case of a 2D seismic image to be interpreted, the seismic learning image may for example correspond to a window in the overall seismic image, of larger dimension. In the case of a 3D seismic image to be interpreted, the learning operation may be carried out on a series of 2D seismic sections chosen from the 3D volumetric image, or else one or more 3D sub-volumes in this 3D seismic volume.

(30) Step 7, which is also optional, relates to the exploitation of the subterranean formation for which a seismic image has been interpreted by at least steps 2 to 4.

(31) The steps for nonlimitingly implementing the method according to the invention described above are detailed below with steps 1, 5 and 7 being optional.

(32) 1. Determining a Priori Probabilities in the Seismic Image

(33) This step is optional, since the a priori probabilities of belonging to the type of geological object that it is sought to be detected may have been determined beforehand and are then provided as input to the method according to the invention.

(34) In this step, it is a question of defining regions within the seismic image under study (which may also be a seismic learning image) and of assigning a probability of belonging to the sought-after geological object to these regions, or in other words an a priori probability of belonging to the type of geological object that it is sought to be detected.

(35) According to one implementation of the invention, an automatic method for defining geological objects is used, which is an alternative to the method according to the invention, of, for example, in (i) calculating a seismic attribute chosen to be (at least partially) discriminating for at least some of the sought-after objects; (ii) normalizing the values of the attribute between 0 and 1; (iii) assigning this result to each region of the seismic image as an a priori probability.

(36) According to another implementation of the invention, it is also possible to determine a priori probabilities of belonging to the sought-after type of geological object by using types of information other than seismic information, such as for example well logs from at least one well passing through the formation under study. According to one implementation of the invention, it is for example possible: (i) for each well for which a well log has been made, to determine a geological facies log by processing well logs conventionally and deducing therefrom a proportion of the object to be detected present with depth down the well under study; (ii) to propagate this determined information to the various wells throughout or in part of the space covered by the seismic image to be interpreted, by interpolation or by a geostatistical method; (iii) to assign this result, normalized between 0 and 1, to each region of the seismic image, as an a priori probability.

(37) According to one implementation of the invention, regions in the seismic image under consideration (which may optionally be a learning image) are delimited by use of at least one of manual and automatic interpretation of the geological objects in the seismic image under consideration and a probability of belonging to the sought geological object is assigned to each of these regions. According to one implementation of the invention: a value of a priori probability of belonging of 1 is assigned for the samples of the seismic image belonging to regions for which it is thought that the delimited region corresponds to the sought type of geological object; and a value of a priori probability of belonging of 0 is assigned for the samples of the seismic image belonging to regions for which it is thought that the delimited region does not correspond, with certainty, to the sought-after type of geological object; intermediate probability values (for example values of 0.25, 0.5 or 0.75) are assigned for the other regions, according to the likelihood of them belonging to the sought type of geological object.

(38) 2. Determining Seismic Attributes in the Seismic Image

(39) In this step, an operation of transforming the seismic image (which may also be a learning seismic image) into a plurality of seismic attributes is defined and the transform thus defined is applied to determine a plurality of seismic attributes that are associated with the seismic image under consideration. Advantageously, at least four seismic attributes are determined, which potentially allows 2.sup.4=16 different types of seismic facies to be distinguished, and thus geological objects having particularly heterogeneous seismic properties to be characterized. Advantageously, attributes are chosen so that they represent the variability of the seismic facies in the seismic image under consideration.

(40) Preferably, the determined seismic attributes comprise textural attributes, which are attributes for which the value in each sample characterizes the texture of the image in the vicinity of this sample. The value in each sample is generally a statistical quantity of the distribution of seismic amplitude values in a lateral and vertical vicinity of the sample of the seismic image (for example taking into account second-order statistics in the vicinity of each sample of a seismic image). Thus, this type of attribute is particularly suitable for distinguishing between various types of seismic facies in a seismic image since they provide information on the local organization of the reflectors in an image at a fine level of detail. Among textual attributes, attributes from the grey-level co-occurrence matrix (GLCM) are known. GLCM attributes were initially introduced (see for example the document (Haralick, 1973)) for image processing in general, and now see wide use in seismic interpretation. An attribute is calculated in a greyscale image, for a predefined orientation theta and inter-pixel distance d. This calculation is performed in two steps: first, the co-occurrence matrix M(r,c) is constructed, which indicates, for each grey-level pair (r,c), the number of pairs of pixels in the adjacency window of a pixel (or sample of the seismic image) that are separated by a vector of norm d and of angle theta and have grey levels r and c; then statistical quantities are calculated using this matrix. In general, textural attributes also have the advantage of being calculated according to one and the same calculation procedure, to the extent that in calculating them all, all of the possible patterns are scanned (for example, with the GLCM attributes, an attribute is calculated per orientation, and they are all calculated in the same way).

(41) According to another implementation of the invention, the predetermined seismic attributes comprise spectral attributes, this type of attribute additionally providing information on size scales. According to another implementation of the invention, it is also possible to use seismic attributes that are conventionally used in seismic interpretation, such as coherence or autocorrelation.

(42) According to one implementation of the invention, if the values of the determined seismic attributes are not of the same order of magnitude (for example if there is a ratio of 1 to 10 between two attributes), the values of the attributes thus determined are normalized.

(43) According to one implementation of the invention, if it turns out that the attributes thus determined are correlated, an operation of decreasing the dimension of the determined attributes is performed. For this, a feature (or attribute) extraction method, such as principal component analysis (PCA) or a feature (or attribute) selection method, may be applied to the determined attributes. According to one advantageous implementation of the invention in which the PCA method is applied, only those attributes which are projected into the reference frame of the eigenvectors from the PCA, the sum of the normalized eigenvalues of which is greater than 0.95, are retained. According to one implementation of the invention in which an attribute selection method is applied, a selection of the initial attributes is retained, the selection being defined such that the selected attributes represent the total variance (or a percentage of 0.95 for example) of the scatter plot and exhibit the smallest possible degree of correlation (optimizing the selection of attributes therefore minimizes the correlation of the selected set while maximizing the variance that is represented).

(44) 3. Classifying the Seismic Image in the Attribute Space

(45) In this step, an unsupervised-type classification method is applied to the seismic attributes determined upon completion of step 2 described above, that is potentially to the attributes after dimension decrease. Thus, in this step, the seismic image in the multi-attribute space, and no longer the seismic image expressed in terms of seismic amplitude, is considered.

(46) According to the invention, the attributes from the seismic image are classified by use of a self-organizing map (SOM) algorithm such as described for example in (Kohonen 1989), or else a generative topographic mapping (GTM) algorithm such as described for example in (Bishop, 1998), or else a K-means algorithm such as described for example in the document (Lloyd, 1982). In the case of an SOM or GTM method, the classification operation generates a non-linear representation of the data by of a manifold of small dimension (typically 2), in which a (potentially irregular) grid containing the class centres is located. During the implementation of the SOM or GTM algorithm, the positions of the class centers in the manifold and in the data space are optimized (and in the case of GTM, the manifold itself is parameterized). For this type of approach, the number of classes is predefined in advance. The way in which a number of classes is defined ad hoc so as to implement an unsupervised-type classification method is known.

(47) Upon completion of this classification step, a class is assigned to each sample (a pixel in a 2D image or a voxel in a 3D seismic image) of the seismic image under consideration. In other words, each sample of the seismic image under consideration possesses a class “label”, or else an identifier of the class to which it belongs. Conventionally, the identifier of the class of each sample of the seismic image is a class number.

(48) 4. Determining an a Posteriori Probability of Belonging

(49) According to the invention, it is a question of assigning, to each class obtained upon completion of the preceding step, an a posteriori probability of belonging to the sought type of geological object. This a posteriori probability is estimated according to: a priori probabilities of belonging to the sought geological object, determined for example as described in step 1; the class division determined upon completion of step 3; and a parameter α, which is predefined (in particular by a specialist), that allows the confidence in the quality of the a priori interpretation of the geological objects in the seismic image (which may also be a learning seismic image) to be qualified. According to one implementation of the invention, the parameter α for a given seismic image is between 0 and 1, the value 0 being attributed in the case of a highly doubtful a priori interpretation of the geological objects and, conversely, 1 being attributed in the case of the confidence in the interpretation being high. The parameter α may also be used to afford the detection algorithm according to the invention a certain “flexibility” (0 for a high degree of flexibility and 1 for no flexibility). The parameter α may also make it possible to take into account the potential absence of annotation for certain objects, for example if it is thought that some geological objects, corresponding to the type of object to be detected, have not been interpreted a priori in the (potentially learning) seismic image.
The following notation is used hereinafter: X is the vector containing a value of a priori probability of belonging for each sample of the seismic image under consideration; L is the vector containing, for each sample of the seismic image under consideration, the identifier of the class to which the sample belongs; S={s.sub.1, s.sub.2, . . . , s.sub.n}={s.sub.k, k∈1:n} with 0=s.sub.1<s.sub.2< . . . <s.sub.n=1, is the set of values taken by the elements of X; C.sub.i, i∈L is the number class i of the classification described by L; p.sub.k,i=p(s.sub.k|C.sub.i) is the proportion of points of C.sub.i having an a priori probability value of s.sub.k; Y is the vector containing, for each sample of the image under consideration, the a posteriori probability.

(50) According to one preferred implementation of the invention, the a posteriori probability y.sub.i of belonging to the sought type of geological object is determined for each number class i of the classification described by L as described below:

(51) $\mathrm{if}{p}_{1,i}>\frac{1}{\alpha +1},y_i=s_1.$

(52) $\mathrm{If}{p}_{1,i}\le \frac{1}{\alpha +1},$
the following is determined:

(53) $q_i=\max \left\{k\in 1\text{:}n\mathrm{such}\mathrm{that}\underset{j\in 1:k}{.Math.}{p}_{j,i}\le \frac{1}{\alpha +1}\right\},\mathrm{and}$$a.\mathrm{if}{.Math.}_{j\in 1:q}{p}_{j,i}=\frac{1}{\alpha +1},y_i=s_q.b.\mathrm{if}{.Math.}_{j\in 1:q}{p}_{j,i}<\frac{1}{\alpha +1},\mathrm{if}q=n,y_i=s_n$$\mathrm{if}q<n,y_i=s_{q+1}.$

(54) Thus, the method according to the invention allows both an a priori probability of belonging and information on the structure of the dataset, represented by the result of the classification operation, to be employed while taking into account a degree of confidence in the a priori interpretation, which is predefined.

(55) Upon completion of this step, a single a posteriori probability value is obtained per class, which is close to the arithmetic mean of the a priori probabilities of all of the samples of the class (in particular, it is equal to the arithmetic mean in the case that the value of the parameter α is equal to 1) but which may include a certain degree of flexibility, via the parameter α, according to the confidence in the a priori interpretation.

(56) Next, an a posteriori probability of belonging to the sought geological object is assigned to each sample of the seismic image under consideration according to the class (for example identified by a class identifier) to which the sample under consideration belongs.

(57) 5. Propagating the a Posteriori Probabilities to Another Seismic Image

(58) This step is optional. It is applied in the case that the a priori probabilities of belonging to the sought-after geological object have been predefined using a learning seismic image, such as defined above, and that it is desired to automatically detect geological objects of the type under consideration in another seismic image, referred to hereinafter as the “overall seismic image”. It is obvious that the learning seismic image and the overall seismic image should have, preferably and as far as possible, similar seismic parameters. For example, the seismic data on which these two seismic images are based should preferably have been recorded using the same acquisition parameters (in particular the sampling intervals) or have been corrected so as to get them to the same acquisition parameters (this type of correction is conventionally carried out in 4D seismic processing for example). Additionally, the seismic data on which these two seismic images are based should preferably have undergone the same seismic processing or have been corrected so as to get them to the same seismic processing (this type of correction is conventionally carried out in 4D seismic processing for example); specifically, it is important that the seismic amplitudes of the learning images and of the overall image to be interpreted be in similar ranges. These constraints and the obtaining of them are known, using routine seismic processing methods.

(59) Upon completion of step 4, a posteriori probabilities of belonging to the sought geological object, associated with each of the classes determined using the learning seismic image, are obtained in particular. In this step, it is a question of “propagating” these a posteriori probabilities, determined using the learning seismic image, to the overall seismic image.

(60) According to the invention, this propagation operation is performed as follows: I. The seismic attribute transform, such as defined in step 2 and applied to the learning seismic image, is applied to the overall seismic image. Next, a plurality of seismic attributes of the overall seismic image are determined, which attributes are the same (in terms of type and number) as those determined using the learning image. Preferably, the post-processing treatments (normalization, dimension decrease; see step 2 above) applied to the attributes determined for the learning seismic image are advantageously also applied to the attributes of the overall seismic image, using the same post-processing operation parameters as those used for the learning seismic image. II. The unsupervised classification method, such as defined in step 3 and applied to the attributes of the learning seismic image, are applied to the attributes of the overall seismic image. Here, the same classification method as that applied in step 3 to the learning image, parameterized in the same way, and in particular with the same number of classes, is applied. A class identifier for each sample of the overall seismic image is thus obtained, for example. III. For each of the classes of the overall seismic image and for each of the samples of this overall seismic image belonging to this class (for example identified by a class identifier), the a posteriori probability of belonging determined for this class using the learning seismic image is applied. Thus, upon completion of this step, a posteriori probabilities of belonging to the sought geological object are obtained for an overall seismic image on the basis of a posteriori probabilities of belonging to the sought geological object that were determined using a learning seismic image. In particular, this makes possible saving a substantial amount of time when large volumes of seismic data are to be interpreted, since the a posteriori probabilities of belonging are determined finely only for the learning seismic image. This implementation of the invention also allows the size of the computer memory used to determine the a posteriori probabilities of belonging to be limited.

(61) 6. Automatically Detecting Geological Objects

(62) In this step, it is a question of automatically detecting, on the basis of the a posteriori probabilities of belonging to the sought geological object that were determined in step 4 and/or 5, geological objects of the type sought after in the seismic image of interest.

(63) This detection operation may be applied by use of any type of known method for automatically detecting geological objects having homogeneous internal properties, since the a posteriori probability values that have thus been determined may then be used as attributes for the implementation of this type of method. In other words, the method according to the invention makes it possible, on the basis of highly heterogeneous and complex seismic information, to determine an image having more homogeneous properties, making the application of methods for detecting geological objects having homogeneous internal properties possible.

(64) According to one implementation of the invention, it is possible, on the basis of the a posteriori probabilities of belonging for each sample of the seismic image under consideration, to automatically detect the geological objects in this image of the type of that sought by means of the following methods: thresholding is applied to the a posteriori probability values: for example, only those samples of the seismic image having a posteriori probability values that are higher than 0.8 are retained; and/or a region-growing algorithm is used, implemented with a predefined criterion pertaining to the values of a posteriori probabilities of belonging that are expected inside objects of the type of object to be detected. Reference may be made to the document (Adams and Bishof, 1994), which describes this type of method. In general, a region-growing algorithm defines seed points for the regions, then in growing these regions by measuring similarity between samples until a stop criterion is met. It is well known that these methods perform best when the seed points are chosen in homogeneous zones of the image to be interpreted. An implicit or explicit edge detection algorithm is used. Reference may be made to the document (Mumford and Shah, 1989), which describes this type of method in general.

(65) 7. Exploiting the Subterranean Formation

(66) From at least the geological objects detected by implementing at least steps 2 to 4 and 6 described above for at least one type of geological object to be detected, valuable information relating to the subterranean formation which is sought to exploit, for example, hydrocarbons is obtained.

(67) Specifically, a conventional approach is in particular to construct a meshed representation of the formation under study, which representation is known as a geomodel, for the purpose of exploiting this formation. Each mesh of this meshed representation is filled with one or more petrophysical properties, such as porosity, permeability, the type of sedimentary deposit, etc. This meshed representation reflects in particular the geometry of the geological objects encountered in the formation under study. It is thus most frequently structured in layers, with a group of meshes being assigned to each geological layer of the modelled basin. However, it advantageously reflects the presence of smaller-scale geological objects (such as MTDs, fluvial channels, etc.) by use of meshes having petrophysical properties that are distinct from the layers in which these objects are embedded. In general, this model is constructed on the basis of data acquired in seismic surveys, from well logs, from core samples, etc. Thus, the information obtained by implementing the method according to the invention, by improving knowledge of the geological objects that are present in the formation, will contribute to the construction of a more accurate geological model of the formation under study.

(68) From such a meshed representation, in particular zones of interest are selected in the subterranean formation, and in particular geological reservoirs that are likely to be highly productive in hydrocarbons (on the basis of criteria such as high porosity, high permeability, the presence of a cap rock, etc.). From such a meshed representation, a prediction may be made of the movements of fluid through the formation under study and to plan its future development by determining plans for exploiting the formation. In particular, determining a plan for exploiting a subterranean formation comprises defining a number, a geometry and siting (position and spacing) of injection and production wells, determining an enhanced recovery type (water, surfactant injection, etc.), etc. A plan for exploiting a hydrocarbon reservoir should for example allow a high rate of recovery of hydrocarbons trapped in the identified geological reservoir, over a long period of exploitation, requiring a limited number of wells and/or amount of infrastructure. Conventionally, a hydrocarbon exploitation plan is determined using a numerical flow simulator. One example of a flow simulator is the PumaFlow® software (IFP Energies nouvelles, France).

(69) Next, once an exploitation plan has been defined, the hydrocarbons trapped in the reservoir are exploited as per this exploitation plan, in particular by drilling at least one of the injection and production wells according to the determined exploitation plan, and by installing the production infrastructure required for the development of the deposit.

(70) Equipment and Computer Program Product

(71) The method according to the invention is implemented by use of equipment (for example a computer workstation) comprising data processing (a processor) and data storing (a memory, in particular a hard disk), and an input/output interface for inputting data and returning the results of the method.

(72) The data processor is in particular configured to carry out the following steps: a plurality of seismic attributes of the seismic image are determined; an unsupervised classification method is applied to the seismic attributes; an a posteriori probability of belonging to the sought type of geological object is determined for the samples of each of the classes; an a posteriori probability of belonging to the sought type of geological object is determined for each of the samples of the seismic image; and and at least one geological object of the type sought is detected on the basis of the a posteriori probabilities of belonging for each sample of the seismic image.

(73) Furthermore, the invention further relates to a computer program product stored on a tangible recording medium that is readable by computer and executable by a processor, comprising program code instructions for implementing the method such as described above, when the program is executed by the processor.

(74) Variant: Prior Acquisition of a Seismic Image

(75) The method according to the invention is implemented using at least one seismic image of a subterranean formation. According to one variant embodiment of the invention, a seismic image is determined prior to the implementation of the method according to the invention. The acquisition of such a seismic image comprises a seismic measurement step, then a seismic processing step of processing these seismic measurements.

(76) The seismic measurement step is carried out by a seismic measurement device. The seismic measurement device comprises means such as explosives or a vibrator for terrestrial seismology, or an air gun or a water gun for marine seismology which emits one or more seismic waves into the subterranean formation of interest, as well as means (such as acceleration sensors (seismometers), vibration sensors (geophones), pressure sensors (hydrophones), or by a combination of elementary sensors of the preceding types (for example multicomponent sensors)) for recording the waves thus emitted and being at least partially reflected in the formation under study. The equipment for acquiring seismic measurements produces a seismic image of a subterranean formation.

(77) The seismic processing step is applied to the seismic measurements taken of the subterranean formation under study. Specifically, the seismic measurements recorded as described above are very often unusable. Conventionally, the seismic processing operation may comprise a step of correcting the recorded seismic amplitudes, a deconvolution operation, statics corrections, (random or coherent) noise filtering, NMO (normal moveout) correction, stacking and migration (depthwise or in time, before or after stacking). These processing steps, requiring calculations that are often long and complex, are carried out by computer. The resulting seismic data are the seismic image. These seismic images are most frequently depicted on a computer, by means of a mesh or grid. Each mesh corresponds to a lateral and vertical position (the vertical direction corresponding to the time or to the depth at which the processing operation resulted in a time image or a depth image) within the formation under study, and further being characterized by a seismic amplitude. The equipment for applying seismic processing ad hoc to seismic measurements to generate a seismic image of a subterranean formation, which image is intended for seismic interpretation is known.

EMBODIMENT EXAMPLES

(78) The features and advantages of the method according to the invention will become more clearly apparent from studying the following application examples.

First Example

(79) In this example it is a question of illustrating, using a purely theoretical example, the principle of the method according to the invention.

(80) FIG. 1 shows a geological object (bold outline) which has been manually delimited in a seismic image, along with the classes (classes C1 to C6) resulting from the application of an applied classification method. The a priori probabilities are discretized as follows [s.sub.1,s.sub.2,s.sub.3,s.sub.4,s.sub.5]=[0, 0.25, 0.5, 0.75, 1]. The probability s.sub.4 has been attributed to the interior of the outlined zone, and the probability s.sub.1 has been attributed to the samples outside this outline. Thus, for this example, only two probability values are shown in the overall distribution of a priori probabilities.

(81) In this example, the risk of having forgotten to interpret a geological object of the sought type in the seismic image is considered to be around 40%. According to the embodiment defined above, this amounts to setting α=0.6 (recall that α=1 corresponds to the case in which the interpreter is certain to have labelled all of the objects).

(82) The a posteriori probabilities per class are examined below.

(83) Class 1: Class 1 corresponds to the pattern C1 from FIG. 1. The probability of the previously produced outline corresponding to points of the sought object has been estimated as 0.75. Outside this surface, the probability is 0. It is therefore now necessary to estimate the proportion of points of class 1 in each of the a priori probability zones, i.e. the proportion p.sub.4.1=p(s.sub.4|C.sub.1) of pixels of C.sub.1 that are in the a priori probability zone s.sub.4=0.75, and the proportion p(s.sub.1|C.sub.1) of pixels of C.sub.1 that are in the a priori probability zone s.sub.1=0. The other a priori probability values are not represented in C.sub.1, hence the proportions p(s.sub.2|C.sub.1), p(s.sub.3|C.sub.1) and p(s.sub.5|C.sub.1) are all zero. For this example, the area of the segmented surface having an a priori probability of 0 is small (see hatched portion in FIG. 2a). The area of the segmented surface having an a priori probability of 0.75 is larger (see hatched portion in FIG. 2b). The calculation of the a posteriori probability associated with class 1, takes this distribution of values into account and is illustrated schematically in FIG. 2c. For each a priori probability value s.sub.k, the proportion of points of class 1 having this value s.sub.k is represented by a bar surmounted by a grey ring. The cumulative value of these proportions (Σ.sup.j∈1:kp.sub.j,1) is represented by a dashed line, and the black ring represents the first point of the cumulative distribution curve that is higher than

(84) $\frac{1}{\alpha +1}.$
Specifically, the value of

(85) $\frac{1}{\alpha +1}$
determines the threshold for the search for q as defined in section 4 above. Since

(86) $p_1\le \frac{1}{\alpha +1},q$
is determined, (in this instance q=3) and the a posteriori probability y.sub.1 of class C.sub.1 is therefore s.sub.q+1=0.75.

(87) Class 2 The implementation of the method according to the invention for class 2 (pattern C2 in FIG. 1) is identical to that presented for class 1. FIG. 3a illustrates the distribution of the samples of the image belonging to class 2 (see pattern C2 in FIG. 1) with respect to the previously interpreted object (delimited by the bold outline). The calculation of y.sub.2, the a posteriori probability associated with class 2, takes this distribution of values into account and is illustrated schematically in FIG. 3b. The threshold value

(88) $\frac{1}{\alpha +1}$
is compared with the cumulative distribution curve (Σ.sup.j∈1:kp.sub.j,1) and the a posteriori probability y.sub.2 of class C.sub.2 is therefore s.sub.q+1=0.75.

(89) Classes 4, 5 and 6

(90) The method according to the invention is applied to classes 4, 5 and 6 (patterns C4, C5 and C6 in FIG. 1). FIG. 4a illustrates the distribution of the samples of the image belonging to these classes (see patterns C4, C5 and C6 in FIG. 1). In these three cases, the intersection between the zones for which the a priori probability of the presence of an object of interest and the zones resulting from the classification step is void. Thus,

(91) $p_1>\frac{1}{\alpha +1},$
which implies that y.sub.4,5,6=0. In other words, since the first point of the cumulative distribution curve is above the threshold

(92) 0$\frac{1}{\alpha +1},$
all of the points of these zones take, as the a posteriori probability, the abscissa of this first point, namely s.sub.1=0, as illustrated in FIG. 4b.

(93) Class 3 Lastly, class 3 (see pattern C3 in FIG. 1) matches only one a priori probability surface, as shown in FIG. 5a, but, unlike classes 4, 5 and 6, all of its points are found at the a priori probability s.sub.4=0.75. The cumulative distribution curve (see FIG. 5b) leaves zero only at 0.75. Consequently, s.sub.q=0.5 and the a posteriori probability retained for class 3 is s.sub.q+1=0.75.

Second Example

(94) For this example, the geological objects are to be delimited in a 3D seismic image (a 3D seismic cube after stacking, also called a “post-stack” seismic cube) and correspond to MTD (mass-transport deposit)-type sedimentary deposits. These objects can be recognized in the seismic image by an informed interpreter by the various seismic facies that they contain and their arrangement. The seismic image has a vertical sampling interval of 4 ms, and an intertrace distance of 25 m in both (crossline and inline) directions. The MTDs that it is sought to detect have a size of several tens of kilometres laterally, and are at most 100 ms thick. Three 2D sections of the cube are used as learning seismic images (see FIG. 6, which shows these three learning seismic images A, B and C). A pixel of a seismic image therefore has a size of 4 ms×25 m. These three seismic sections have been interpreted by an interpreter, who has delimited zones with a very high a priori probability of belonging P to an MTD-type of geological object (white zones in FIG. 7 for the three learning seismic images A, B and C), a lower probability (dark grey zones in FIG. 7 for the three sections A, B and C) and zones which do not comprise MTD (black zones in FIG. 7 for the three learning seismic images A, B and C).

(95) For this illustrative example, textural attributes of GLCM type, as defined above, are used. In particular, the textural attributes conventionally called “contrast”, “correlation”, “energy” and “homogeneity” are used (see for example the documents (Haralick, 1976; West et al., 2002)). A study window (the vicinity around the pixel, in which the co-occurrence matrix will be calculated) of 11×11 pixels is chosen. Two different scales are used: distances of d=1 and of d=2 pixels. For each scale, all possible orientations theta are used (four possible orientations for scale 1 and eight for scale 2). For the calculation, it is necessary to limit the number of possible levels for the seismic amplitude values (grey levels). This number is often 256, 64, 8 or lower; in the present case only two grey levels are chosen. The calculations are thus shorter, since the time for calculating the GLCM attributes increases with the square of the number of grey levels. To compensate for the loss of amplitude information, an attribute from the envelope of the seismic signal is added to the set of GLCM attributes: for a pixel, it is the result of Gaussian-filtering the image of the seismic envelope, with a 2D Gaussian kernel, close to the window of 11×11 pixels in size. The GLCM attributes, as well as this attribute based on the envelope, are calculated for all of the pixels of the learning image. The attributes are then normalized.

(96) For this implementation of the invention, the dimension of the attribute space is decreased. For this, principal component analysis (PCA) is used, which is applied separately to the attributes from the scales d=1 and d=2. The attribute from the envelope is left outside the PCA. For each of the two samples, only those eigenvectors for which the cumulative sum of the normalized eigenvalues is higher than 0.97 are retained. There are then three new attributes for the scale d=1 and six for the scale d=2. Weightings of 0.5 are imposed on the attributes of scale 2 (to ensure an equal representation of scales 1 and 2), and the attribute from the envelope, with a weighting of 3, is added to this set of nine new attributes.

(97) Thus, the base of the new attribute space consists of these 10 eigenvectors: the first three eigenvectors from the PCA at scale 1, the following six (with a weighting of 0.5) from the PCA at scale 2, and the last from the envelope, with a weighting of 3.

(98) Next, a GTM (generative topographic mapping)-type unsupervised classification algorithm, such as defined above, is applied in its 2D version (the manifold generated will be of dimension 2). Thus, for each pixel of each learning image, a class label from among 49 classes (see FIG. 8, which shows the distribution of the 49 classes in each of the learning seismic images A, B and C) is determined.

(99) Next, the preferred implementation of the method according to the invention is applied, as described in step 4, to assign an a posteriori probability of the presence of the objects to each class. The parameter α is set at 0.8 (degree of confidence in the interpretation). The associated a posteriori probability is then determined for each class. FIG. 9 shows the a posteriori probabilities P′ thus obtained for each of the learning seismic images A, B and C.

(100) Next, the a posteriori probabilities thus determined for each of the learning seismic images are “propagated” to a series of 72 parallel 2D sections extracted from the 3D seismic cube, which are, advantageously, parallel to the learning images. The sections are extracted with a sampling interval of one section every 250 m (one section out of 10), which is relatively fine with respect to the size of the sought objects. FIG. 10 shows a portion of these sections depicted as a posteriori probabilities P′.

(101) For illustrative purposes and to determine a complete seismic volume, each section has been replicated 10 times, to represent the overall 3D seismic volume. An interpolation of the values of the a posteriori probabilities between each section could advantageously have been carried out. Next, a post-processing treatment is applied by thresholding the values of the probabilities. More specifically, for this example, only those pixels for which the a posteriori probability values are equal to 1 are retained, followed by thresholding on the basis of the size of the objects and the distribution of the class numbers inside the objects.

(102) The delimitation of the sought-after MTDs is achieved, as presented in FIG. 11. It may be seen that the method according to the invention has thus made it possible (i) to automatically delimit these objects without omitting portions due to the difference in facies, and (ii) to retain the information relating to these various facies within the delimited object itself, which is subsequently useful when it comes to interpreting the properties of the object.