Patent Description:
It is known, especially in oil exploration, to determine the position of oil reservoirs from the results of geophysical measurements carried out from the surface or in well bores.

According to the technology of reflection seismology, these seismic measurements involve emitting a wave (e.g. acoustic waves) into the subsurface and measuring a signal comprising a plurality of echoes of the wave on geological structures being investigated. These structures are typically surfaces separating distinct materials, faults, etc. Other measurements may be carried out at wells bores.

Chrono-stratigraphic analysis (or sequence stratigraphic analysis) is very important to understand basin evolution, predict the sedimentary facies distribution for both hydrocarbon exploration and development. This analysis is based on the fundamental assumption that seismic reflectors are surfaces of chrono-stratigraphic significance. This assumption implies that an individual seismic reflector is a "time-line" through a depositional basin that represents a surface of the same geological age (i.e. an isochronous surface in geological time).

A seismic image (or seismic section) comprises a juxtaposition in a volume of sampled one-dimensional signals referred to as seismic traces. In the seismic image, the value of a pixel (a. voxel for 3D images) is proportional to the seismic amplitude represented by seismic traces.

Computing a chrono-stratigraphic representation of a seismic image often requires determining seismic horizon surfaces of the seismic image, wherein a seismic horizon surface corresponds to an estimated isochronous surface of the geological formation. Such seismic horizon surfaces can be used to determine an RGT image of the geological formation, i.e. an image in which each pixel provides an estimated geological age for the portion of the geological formation represented by said pixel. For example, [WU2015] describes methods for computing a RGT volume by extracting horizon surfaces using control points. The RGT image is referred to as "relative" because the purpose of the RGT image is mainly to be able to compare the estimated geological ages of different pixels, in order to e.g. identify portions of the geological formation that have the same estimated geological age. Also, in practice, it is usually not possible to estimate an absolute geological age of any given portion of the geological formation.

[LOMASK2006] describes a method for determining seismic horizon surfaces based on a seismic image, by computing the local seismic dip at each pixel of the seismic image and searching iteratively for surfaces having local gradients approaching the local seismic dips. This method relies on a regular grid within the seismic image, usually mapped to the size of the pixels of the seismic image (i.e. the seismic resolution), required for using fast Fourier transforms enabling to reduce computational complexity. However, due to the regular grid used, it is not possible in this method to introduce prior knowledge on e.g. where the seismic reflectors are located and/or on where faults are located. For this reason, it is not possible in this method to ensure e.g. that a determined seismic horizon surface corresponds to an actual seismic reflector. The determined seismic horizon surface is any surface deemed to be isochronous in geological time and may as well be a seismic reflector or an isochronous surface located between seismic reflectors. Also, this method will perform poorly in the presence of faults in the geological formation.

The present disclosure aims at improving the situation. In particular, the present disclosure aims at overcoming at least some of the limitations of the prior art discussed above, by proposing a solution for determining RGT surfaces, i.e. isochronous surfaces of the geological formation, by considering some prior knowledge on the geological formation.

According to a first aspect, the present disclosure relates to a computer implemented method for processing a seismic image obtained from seismic measurements performed on a geological formation, said method comprising:.

Also, combining of the vertical positions of the points of the vertical surfaces uses weighting coefficients, wherein the value of a weighting coefficient used for the vertical position of an RGT surface decreases with a distance between the at least one horizontal position of the considered point of the reference RGT surface and the at least one horizontal position of the setpoint associated to the considered RGT surface.

Hence the processing method computes an RGT image based on the seismic image. Computing the RGT image may use any method known to the skilled person, including any method for determining seismic horizon surfaces of the seismic image known to the skilled person.

Then setpoints are obtained, which represent constraints to be satisfied by a reference RGT surface. For instance, the setpoints may be determined based on measurements performed at a well bore in the geological formation or may correspond to points positioned manually in the seismic image and/or in the RGT image by a human interpreter, which are considered by the human interpreter to correspond to different points of a same seismic reflector.

For each setpoint, an RGT surface is retrieved from the RGT image, i.e. estimated isochronous surfaces are retrieved for each setpoint received.

Then the reference RGT surface, i.e. the RGT surface that takes into account the setpoints received, is obtained by combining the RGT surfaces retrieved for the received setpoints. Basically, for a given horizontal position, the vertical position of a point of the reference RGT surface is influenced mainly by the vertical position of the point of the RGT surface retrieved for the setpoint that is closest in the horizontal plane with the considered horizontal position.

Hence, the setpoints, i.e. the constraints / prior knowledge on the geological formation, are not used for determining the seismic horizon surfaces and/or the RGT image but are used afterwards to combine RGT surfaces retrieved from the RGT image for the received setpoints. By taking into account the setpoints when determining the reference RGT surface, better results can be obtained e.g. in the presence of faults in the geological formation.

In specific embodiments, the processing method can further comprise one or more of the following features, considered either alone or in any technically possible combination.

In specific embodiments, obtaining setpoints for the RGT image comprises obtaining a patch of points of the RGT image, said patch of points having a boundary, wherein the setpoints comprise points located on the boundary of the patch of points, and wherein the reference RGT surface between setpoints located on the boundary of the patch of points is obtained from said patch of points.

In such a case, it is possible to consider constraints that correspond to patches of points. The patches of points may be included as such in the reference RGT image, while the RGT surfaces are combined to obtain the points of the reference RGT image which are located e.g. between different patches of points, the setpoints being points on the boundaries of the patches of points.

In specific embodiments, obtaining a patch of points comprises defining a seed on a seismic trace of the seismic image and propagating laterally the seed by searching for pixels of adjacent seismic traces that satisfy a predetermined criterion of similarity with the seed.

In such a case, the seeds may be determined based on measurements performed at a well bore in the geological formation or may correspond to points positioned manually in the seismic image and/or in the RGT image by a human interpreter, which are considered by the human interpreter to correspond to different points of a same seismic reflector. Then these seeds are propagated laterally by e.g. correlation between adjacent seismic traces of the seismic image in order to produce the patches of points. Hence, this method may be used to combine different approaches such as correlation-based methods and seismic dip based methods (such as [LOMASK2006]).

In specific embodiments, the value of a weighting coefficient used for the vertical position of an RGT surface is:.

In specific embodiments: <MAT> wherein: <MAT> <MAT> wherein:.

In specific embodiments, α > <NUM> and/or α < <NUM>.

In specific embodiments, the processing method comprises converting to frequency-domain the reference RGT surface.

In specific embodiments, converting to frequency-domain the reference RGT surface comprises:.

In specific embodiments, the processing method comprises selecting three analysis frequencies and generating a red-green-blue, RGB, image of the reference RGT surface based on the three frequency-domain representations associated respectively to the three selected analysis frequencies.

According to a second aspect, the present disclosure relates to a computer program product comprising instructions which, when executed by at least one processor, configure said at least one processor to carry out a processing method according to any one of the embodiments of the present disclosure.

According to a third aspect, the present disclosure relates to a computer-readable storage medium comprising instructions which, when executed by at least one processor, configure said at least one processor to carry out a processing method according to any one of the embodiments of the present disclosure.

According to a fourth aspect, the present disclosure relates to a computer system for processing a seismic image, said computer system comprising at least one processor configured to carry out a processing method according to any one of the embodiments of the present disclosure.

The invention will be better understood upon reading the following description, given as an example that is in no way limiting, and made in reference to the figures which show:.

In these figures, identical references from one figure to another designate identical or analogous elements. For reasons of clarity, the elements shown are not to scale, unless explicitly stated otherwise.

As discussed above, the present disclosure relates inter alia to a method <NUM> for processing seismic images.

A seismic image represents a picture of the subsoil arising from a seismic exploration survey. The seismic image comprises at least two dimensions which may comprise at least one horizontal dimension (which usually uses a distance scale, expressed e.g. in meters) and one vertical dimension (which usually uses a distance scale or a time scale, expressed e.g. in seconds). Hence, the seismic image may correspond to a 3D seismic image (with two horizontal dimensions and one vertical dimension) or to a 2D seismic image (with one horizontal dimension and one vertical dimension).

It is emphasized that the expressions "horizontal dimension" and "vertical dimension" are not be interpreted as requiring these dimensions to be respectively strictly horizontal and strictly vertical. These expressions mean that one of the dimensions, referred to as "vertical dimension", is representative of the depth of the geological formation, and that the other dimensions, referred to as "horizontal dimensions" are both orthogonal to the vertical dimension.

The seismic image is composed of pixels which may be 2D in case of a 2D seismic image or 3D (voxels) in case of a 3D seismic image. The pixels are regularly distributed according to a horizontal resolution on each horizontal dimension and a vertical resolution on the vertical dimension. The seismic image comprises, along each horizontal dimension:.

Each pixel is associated with a seismic value which may be a gray value, for instance between <NUM> and <NUM> (or <NUM>). Each seismic value is representative of the amplitude of the seismic signal measured for the portion of the geological formation represented by the corresponding pixel.

In the sequel, a point corresponds to a set of coordinates in the grid of the seismic image, i.e. a set comprising a horizontal position along each horizontal dimension and a vertical position along the vertical dimension. A pixel therefore corresponds to a point with a value associated thereto (i.e. a seismic value in the case of a pixel of the seismic image, an estimated geological age in the case of a pixel of the RGT image, etc.).

<FIG> represents an example of seismic image. The seismic image represented by <FIG> is 2D and may correspond to a 2D seismic image or to a 2D section of a 3D seismic image, in a vertical plane comprising one of its horizontal dimensions. As can be seen in <FIG>, the seismic values highlight the composition of the geological formation, since high amplitude seismic values are usually associated to strong seismic reflectors, which are usually located at the interfaces between geological layers having different acoustic impedances.

<FIG> represents schematically the main steps of an exemplary embodiment of a method <NUM> for processing a seismic image.

The processing method <NUM> is carried out by a computer system (not represented in the figures). In preferred embodiments, the computer system comprises one or more processors (which may belong to a same computer or to different computers) and storage means (magnetic hard disk, optical disk, electronic memory, or any computer readable storage medium) in which a computer program product is stored, in the form of a set of program-code instructions to be executed in order to implement all or part of the steps of the processing method <NUM>. Alternatively, or in combination thereof, the computer system can comprise one or more programmable logic circuits (FPGA, PLD, etc.), and/or one or more specialized integrated circuits (ASIC), etc., adapted for implementing all or part of said steps of the processing method <NUM>. In other words, the computer system comprises a set of means configured by software (specific computer program product) and/or by hardware (processor, FPGA, PLD, ASIC, etc.) to implement the steps of the processing method <NUM>.

As illustrated by <FIG>, the processing method <NUM> comprises a step S30 of determining an RGT image based on the seismic image.

This step S30 of determining an RGT image may use any method known to the skilled person. For instance, determining the RGT image based on the seismic image may rely on determining seismic horizon surfaces. Preferably a seismic horizon surface may be determined for a plurality of local extrema of a plurality of seismic traces of the seismic image, i.e. each seismic horizon surface determined comprises a local extremum of a seismic trace. Determining seismic horizon surfaces may use any method know to the skilled person for determining seismic horizon surfaces, for instance the method described in [LOMASK2006], in the patent applications <CIT> and <CIT>, etc..

Then, the RGT image may be computed based on the seismic horizon surfaces as follows.

We assume that Mh seismic horizon surfaces τn (<NUM> ≤ n ≤ Mh) have been determined, preferably distributed throughout the vertical dimension of the seismic image. Assuming that the seismic image comprises a horizontal dimension x with Nx pixels, a horizontal dimension y with Ny pixels and a vertical dimension t with Nt pixels, then the seismic horizon surface τn of index n is the following set of points of the seismic image {(i,j, τn(i, j)), <NUM> ≤ i ≤ Nx, <NUM> ≤ j ≤ Ny}.

For example, the value of each pixel of the RGT image may correspond to the number of seismic horizon surfaces that comprise said considered pixel or that comprise any pixel located in the same column as the considered pixel, between the considered pixel and a reference pixel in the same column. The reference pixel on the vertical axis is the pixel of index k = Nt or, preferably, the pixel of index k = <NUM>.

For instance, it is possible to compute a stack image STK. The value of each pixel of the stack image STK corresponds to the number of seismic horizon surfaces that comprise said considered pixel. We can define a function Pos(i, j, k, n) which is such that: <MAT>.

Hence, the function Pos(i, j, k, n) indicates whether the seismic horizon surface τn passes by the pixel having the position (i, j, k). Based on the function Pos(i, j, k, n), the stack image STK may be computed as follows: <MAT> for each <NUM> ≤ i ≤ Nx, <NUM> ≤ j ≤ Ny, <NUM> ≤ k ≤ Nt, or limited to the pixels which are located inside a predetermined survey volume in the seismic image.

Then, assuming that the reference pixel is the pixel of index k = <NUM>, the RGT image may be computed as follows: <MAT> for each <NUM> ≤ i ≤ Nx, <NUM> ≤ j ≤ Ny, <NUM> ≤ k ≤ Nt, or limited to the pixels which are located inside the survey volume. For the purpose of chrono-stratigraphic analysis, it is possible, in some embodiments, to normalize the RGT image by a predetermined reference geological age, such that the maximum value of the pixels of the RGT image is equal to the reference geological age. Hence, in the present example, the pixels which represent the deepest portions of the geological formation, at least, will have their values equal to the reference geological age.

As illustrated by <FIG>, the processing method <NUM> comprises also a step S31 of obtaining at least two setpoints for the RGT image.

The setpoints represent constraints to be satisfied by a reference RGT surface to be determined. For instance, the setpoints may be determined based on measurements performed at a well bore in the geological formation. Preferably, the setpoints may correspond to points positioned manually in the seismic image and/or in the RGT image by a human interpreter, in which case said setpoints may correspond to points considered by the human interpreter to belong to a same seismic reflector.

In order to facilitate the positioning of the setpoints by a human interpreter, it is possible, in preferred embodiments, to compute an RGT seismic image based on both the RGT image and the seismic image. Basically, the RGT seismic image merges the RGT image with the seismic image by superimposing the RGT image on the seismic traces of the seismic image. For instance, different colors may be used for different layers having substantially the same estimated geological age in the RGT image, and different shades of the color may be used for different seismic values of the seismic traces. For instance, if the color "blue" is used for a given layer of pixels having substantially the same geological age, a "dark blue" may be used for a high seismic value and a "light blue" may be used for a low seismic value of the seismic traces. Such provisions are advantageous in that they allow the human interpreter to visualize simultaneously the estimated geological age of each pixel and its seismic value from the seismic image, thereby facilitating the detection of seismic reflectors. The purpose of positioning the setpoints may be to apply slight corrections to the RGT image. The human interpreter uses the RGT seismic image to detect seismic reflectors and provides setpoints that are considered to have the same real geological age while their respective estimated geological ages in the RGT image may be different. Indeed, the presence of a fault in the geological formation may introduce errors in the estimated geological ages of some pixels of the RGT image, and the setpoints may be positioned to introduce some knowledge regarding points that should be considered to have the same real geological age, in order to compensate for said errors.

We assume that Nsetpoint ≥ <NUM> setpoints are obtained. The position of the setpoint of index n (<NUM> ≤ n ≤ Nsetpoint) is given by (in, jn, kn), wherein (in, jn) are the horizontal positions of the setpoint and kn its vertical position.

As illustrated by <FIG>, the processing method <NUM> comprises a step S32 of retrieving an RGT surface of the RGT image for each of the setpoints. If we denote by Sn the RGT surface retrieved from the RGT image for the setpoint of index n (<NUM> ≤ n ≤ Nsetpoint), then it corresponds to the following set of points of the RGT image {(i, j,Sn(i, j)), <NUM> ≤ i ≤ Nx, <NUM> ≤ j ≤ Ny}, wherein <MAT> for each <NUM> ≤ i ≤ Nx, <NUM> ≤ j ≤ Ny, or limited to the pixels which are located inside the survey volume. Sn(i, j) is therefore the vertical position of the point (i, j, Sn(i, j)) of the RGT surface Sn having the horizontal positions (i,j). The RGT surface Sn comprises the setpoint (in, jn, kn) of index n, i.e. Sn(in, jn) = kn.

Of course, the sequence of the steps of the processing method <NUM>, as illustrated by <FIG>, is only illustrative and other sequences are possible. For instance, steps S31 and S32 may be combined. Indeed, it is possible to obtain a first setpoint and to retrieve a first RGT surface associated to said first setpoint. The first RGT surface may then be displayed, e.g. on a seismic image or an RGT seismic image, to enable a human interpreter to position a second setpoint, e.g. in a position where the human interpreter considers that the first RGT surface should be corrected. Then a second RGT surface may be retrieved, etc..

As illustrated by <FIG>, the processing method <NUM> comprises a step S33 of determining a reference RGT surface by combining the RGT surfaces Sn (<NUM> ≤ n ≤ Nsetpoint) determined for the setpoints. Basically, the vertical position Sref(i, j) of a point (i, j, Sref(i, j)) of the reference RGT surface Sref is obtained by combining the vertical positions Sn(i, j) of the points of the RGT surfaces having the same horizontal positions (i,j) as the considered point of the reference RGT surface Sref. In other words: <MAT> wherein wn(i, j) ≥ <NUM> is the weighting coefficient used for the vertical position of the RGT surface Sn for the horizontal positions (i, j).

Also, the value of a weighting coefficient (which is ≥ <NUM>) used for the vertical position Sn(i, j) of an RGT surface Sn decreases with a distance in the horizontal plane between the horizontal positions (i, j) of the considered point (i, j, Sref(i, j)) of the reference RGT surface Sref and the horizontal positions (in, jn) of the setpoint associated to the considered RGT surface Sn. In other words: <MAT> wherein dist((i, j), (in, jn)) is the distance in the horizontal plane between the point having the horizontal positions (i, j) and the horizontal positions (in, jn) of the setpoint of index n (<NUM> ≤ n ≤ Nsetpoint). Any mathematical formula for computing a distance between two points may be used, and in particular any p-norm distance. Preferably, the <NUM>-norm distance (a. Euclidian distance) is used: <MAT>.

Hence, with such weighting coefficients wn(i, j), the contribution to the reference RGT surface Sref of an RGT surface Sn (and its influence thereof) diminishes with the distance to the setpoint associated to said RGT surface Sn. In other words, the RGT surface Sn of index n is considered to be more accurate in the vicinity of its setpoint than any other RGT surface Sm of index m ≠ n.

In preferred embodiments, the value of a weighting coefficient wn(i, j) used for the vertical position Sn(i, j) of an RGT surface Sn is:.

By "substantially equal to one", we mean that <NUM> ≤ wn(in, jn) ≤ <NUM>, or even <NUM> ≤ wn(in, jn) ≤ <NUM>. By "substantially equal to zero", we mean that <NUM> ≤ wn(im, jm) ≤ <NUM>, or even <NUM> ≤ wn(im, jm) ≤ <NUM>, m ≠ n.

According to a non-limitative example, the weighting coefficients wn(i, j) used for weighting the vertical positions Sn(i, j) of an RGT surface Sn are given by the following expressions: <MAT> wherein: <MAT> <MAT> wherein α > <NUM> is a predetermined factor. In preferred embodiments, α > <NUM> and/or α < <NUM>. Basically, the factor α defines how the influence of an RGT surface Sn decreases with the distance to its setpoint, i.e. with dist((i, j), (in, jn)). Hence, the higher the factor α, the quicker the influence of an RGT surface Sn decreases with the distance to its setpoint.

According to another non-limitative example, the weighting coefficients wn(i, j) used for weighting the vertical positions Sn(i, j) of an RGT surface Sn are given by the following expressions: <MAT> wherein: <MAT> <MAT>.

<FIG> represents an example of reference RGT surface Sref obtained based on RGT surfaces Sn. The image represented by <FIG> corresponds to a 2D section of a 3D seismic image, in a vertical plane comprising one of its horizontal dimensions. In <FIG>, we consider two setpoints (Nsetpoint = <NUM>). Part a) of <FIG> represents schematically a first setpoint SET<NUM> and its retrieved RGT surface S<NUM>. Part b) of <FIG> represents schematically a second setpoint SET<NUM> and its retrieved RGT surface S<NUM>. Finally, part c) of <FIG> represents schematically the reference RGT surface Sref obtained by combining the RGT surfaces S<NUM> and S<NUM>. As can be seen in part c) of <FIG>, the reference RGT surface Sref, in this example, passes by the two setpoints SET<NUM> and SET<NUM> (i.e. wn(in, jn) = <NUM> and wn(im, jm) = <NUM> ∀m ≠ n).

In preferred embodiments, obtaining setpoints for the RGT image comprises obtaining a patch of points of the RGT image used as a constraint. Basically, such a patch of points corresponds to a bounded surface constraining the reference RGT surface Sref, i.e. that is to be included in the reference RGT surface Sref as a portion of said reference RGT surface Sref. In such a case, the setpoints for which RGT surfaces Sn are to be retrieved correspond to points located on a boundary of the patch of points. Hence, between setpoints located on the boundary of a same patch of points, the reference RGT surface Sref is equal to said patch of points RGT, and the RGT surfaces Sn retrieved for the setpoints are combined to determine the reference RGT surface Sref outside the obtained patches of points. Hence, all that has been described hereinabove still applies for combining the RGT surfaces Sn outside the patches of points.

A patch of points may for instance be derived from measurements performed at a well bore in the geological formation and/or from seeds positioned manually by human interpreters.

<FIG> represents schematically the main steps of a preferred embodiment of the step S31 of obtaining the setpoints.

As illustrated by <FIG>, the step S31 of obtaining the setpoints comprises a step S310 of obtaining a seed on at least one seismic trace. The seed may be e.g. derived from measurements at a well bore or positioned manually by a human interpreter.

As illustrated by <FIG>, the step S31 of obtaining the setpoints comprises also, for each seed, a step S311 of processing the seismic image (by the computer system) to propagate laterally the seed by searching for pixels of adjacent seismic traces that satisfy a predetermined criterion of similarity with the seed. For instance, the lateral search may be performed by correlating an adjacent seismic trace with vertically translated copies of the seismic trace comprising the seed in order to identify the vertical translation that optimizes the similarity (e.g. maximum correlation value and/or correlation value above a predetermined threshold) with the seismic trace comprising the seed.

The pixels around a seed that satisfy the similarity criterion with said seed, and their respective positions in the seismic image, define the patch of points associated this this seed.

As illustrated by <FIG>, the step S31 of obtaining the setpoints comprises also a step S312 of selecting setpoints on the boundary of each patch of points obtained by processing the seismic image.

Similarity based methods are known to yield results that are accurate in the vicinity of the seeds. Hence, if a seed is positioned on a seismic reflector, similarity based methods are able to identify pixels of adjacent seismic traces which correspond to the same seismic reflector with a good accuracy at least in the vicinity of the seed. However, such similarity based methods perform poorly in the presence of e.g. noise, where seismic dip based methods (such as [LOMASK2006]) usually provide better results. This method may be used to combine different approaches, and the advantages thereof, such as similarity based methods for determining patches of points and seismic dip based methods for determining the RGT surfaces.

<FIG> represents an example of reference RGT surface Sref obtained based on a patch of points and based on RGT surfaces Sn. The image represented by <FIG> corresponds to a 2D section of a 3D seismic image, in a vertical plane comprising one of its horizontal dimensions. Part a) of <FIG> represents schematically a patch of points obtained by propagating laterally a seed SEED<NUM>. In the 2D section represented by <FIG>, the boundary of the patch of points comprises two setpoints SET<NUM> and SET<NUM>, and part a) of <FIG> represents an RGT surface S<NUM>,<NUM> obtained by combining the RGT surfaces S<NUM> and S<NUM> retrieved for the setpoints SET<NUM> and SET<NUM>. Part b) of <FIG> represents schematically a third setpoint SET<NUM> and its retrieved RGT surface S<NUM>. Finally, part c) of <FIG> represents schematically the reference RGT surface Sref obtained by combining the patch of points and the RGT surfaces S<NUM>,<NUM> and S<NUM>. As can be seen in part c) of <FIG>, the reference RGT surface Sref passes by all three setpoints SET<NUM>, SET<NUM> and SET<NUM>, and is equal to the patch of points where the patch of points is defined.

<FIG> represents schematically the main steps of a preferred embodiment of the processing method <NUM>. As illustrated by <FIG>, the processing method <NUM> comprises, in addition to the steps described in reference to <FIG>, a step S34 of converting to frequency-domain the reference RGT surface Sref. For instance, and as illustrated in a non-limitative manner by <FIG>, the frequency-domain conversion step S34 may comprise:.

During step S340, the intersection point of the reference RGT surface Sref with each seismic trace of the seismic image is determined. Basically, for the seismic trace (column of the seismic image) having the horizontal positions (i,j), then the intersection point corresponds to the point (i, j, Sref(i, j)).

We denote by IS the seismic image, which corresponds to the set of pixels {IS(i, j, k), <NUM> ≤ i ≤ Nx, <NUM> ≤ j ≤ Ny, <NUM> ≤ k ≤ Nt}. During step S341, intervals of pixels centered on the intersection points are determined for each seismic trace of the seismic image IS. If we assume that each interval of pixels comprises Nint = (<NUM> × Nint/<NUM>) pixels, the interval of pixels for the seismic trace having the horizontal positions (i, j) may correspond to the set of pixels: <MAT>.

<FIG> represents schematically the intervals of pixels centered on a reference RGT surface Sref. Basically, the intervals of pixels lie between two boundary surfaces Strans<NUM> and Strans<NUM> (represented by white discontinuous lines in <FIG>) obtained by translating the reference RGT surface Sref by respectively (-Nint/<NUM> + <NUM>) and Nint/<NUM> pixels along the vertical dimension.

Then during step S342, the values of each interval of pixels are converted into frequency-domain. This frequency-domain conversion may use any method known to the skilled person, for instance a Fourier transform and preferably a Sparse Fourier Transform, SFT. Hence, a frequency-domain representation is obtained for each interval of pixels, over a plurality of analysis frequencies. For a given analysis frequency, the set of frequency-domain values obtained for all the intervals of pixels corresponds to a frequency-domain representation of the reference RGT surface Sref. Hence, the frequency-domain conversion step S342 provides a plurality frequency-domain representations of the reference RGT surface Sref associated to respective analysis frequencies.

Claim 1:
- Computer implemented method (<NUM>) for processing a seismic image obtained from seismic measurements performed on a geological formation, said method comprising:
- (S30) determining a relative geological time, RGT, image based on the seismic image;
- (S31) obtaining setpoints for the RGT image, said setpoints representing constraints to be satisfied by a reference RGT surface to be determined;
- (S32) retrieving an RGT surface of the RGT image for each of the setpoints, wherein each point of an RGT surface is defined by a vertical position and at least one horizontal position;
- (S33) determining the reference RGT surface by combining the RGT surfaces determined for the setpoints, wherein the vertical position of a point of the reference RGT surface is obtained by combining the vertical positions of the points of the RGT surfaces having the same at least one horizontal position as the considered point of the reference RGT surface;
wherein said combining of the vertical positions uses weighting coefficients;
wherein the value of a weighting coefficient used for the vertical position of an RGT surface decreases with a distance between the at least one horizontal position of the considered point of the reference RGT surface and the at least one horizontal position of the setpoint associated to the considered RGT surface.