EP2883087B1  Method for enhancing the determination of a seismic horizon  Google Patents
Method for enhancing the determination of a seismic horizon Download PDFInfo
 Publication number
 EP2883087B1 EP2883087B1 EP13745135.7A EP13745135A EP2883087B1 EP 2883087 B1 EP2883087 B1 EP 2883087B1 EP 13745135 A EP13745135 A EP 13745135A EP 2883087 B1 EP2883087 B1 EP 2883087B1
 Authority
 EP
 European Patent Office
 Prior art keywords
 horizon
 pseudo
 points
 rectangle
 axes
 Prior art date
 Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
 Active
Links
 230000002708 enhancing Effects 0.000 title claims description 7
 230000001131 transforming Effects 0.000 claims description 74
 238000004422 calculation algorithm Methods 0.000 claims description 53
 230000000875 corresponding Effects 0.000 claims description 36
 238000004590 computer program Methods 0.000 claims description 9
 230000004048 modification Effects 0.000 claims description 9
 238000006011 modification reaction Methods 0.000 claims description 9
 238000003860 storage Methods 0.000 claims description 2
 210000003414 Extremities Anatomy 0.000 description 11
 238000000034 method Methods 0.000 description 9
 238000004364 calculation method Methods 0.000 description 8
 238000000844 transformation Methods 0.000 description 5
 239000011159 matrix material Substances 0.000 description 4
 238000009826 distribution Methods 0.000 description 3
 238000005304 joining Methods 0.000 description 3
 238000005457 optimization Methods 0.000 description 3
 230000001965 increased Effects 0.000 description 2
 238000005259 measurement Methods 0.000 description 2
 239000000203 mixture Substances 0.000 description 2
 239000010749 BS 2869 Class C1 Substances 0.000 description 1
 239000010750 BS 2869 Class C2 Substances 0.000 description 1
 230000002596 correlated Effects 0.000 description 1
 230000001419 dependent Effects 0.000 description 1
 239000006185 dispersion Substances 0.000 description 1
 230000000694 effects Effects 0.000 description 1
 238000007429 general method Methods 0.000 description 1
 230000003993 interaction Effects 0.000 description 1
 239000003921 oil Substances 0.000 description 1
 229920000123 polythiophene Polymers 0.000 description 1
 230000000135 prohibitive Effects 0.000 description 1
 229910052710 silicon Inorganic materials 0.000 description 1
 239000010703 silicon Substances 0.000 description 1
 239000007787 solid Substances 0.000 description 1
Images
Classifications

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
 G06T17/05—Geographic models

 G—PHYSICS
 G01—MEASURING; TESTING
 G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
 G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
 G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
 G01V1/30—Analysis

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
 G06T17/10—Constructive solid geometry [CSG] using solid primitives, e.g. cylinders, cubes

 G—PHYSICS
 G01—MEASURING; TESTING
 G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
 G01V2210/00—Details of seismic processing or analysis
 G01V2210/60—Analysis
 G01V2210/64—Geostructures, e.g. in 3D data cubes
 G01V2210/643—Horizon tracking

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T2207/00—Indexing scheme for image analysis or image enhancement
 G06T2207/20—Special algorithmic details
 G06T2207/20036—Morphological image processing
 G06T2207/20041—Distance transform

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T2207/00—Indexing scheme for image analysis or image enhancement
 G06T2207/20—Special algorithmic details
 G06T2207/20048—Transform domain processing

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T2207/00—Indexing scheme for image analysis or image enhancement
 G06T2207/20—Special algorithmic details
 G06T2207/20068—Projection on vertical or horizontal image axis

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T2207/00—Indexing scheme for image analysis or image enhancement
 G06T2207/30—Subject of image; Context of image processing
 G06T2207/30181—Earth observation
Description
 The invention pertains to the field of methods implemented in order to determine seismic horizons.
 The invention more specifically relates to a method that enhances the determination of a seismic horizon without suffering from some of the drawbacks of the prior art.
 Geological surveys involving generators of seismic waves and detectors of their reflections in the ground are often conducted to determine the position of oil reservoirs and/or to get to know the composition and thickness of the many layers that form the underground. Seismic reflection techniques consist in generating a seismic wave that propagates through the ground and reflects at the interfaces thereof. A precise measurement of these echoes and more specifically of their arrival times enables a determination of the shape, depth and composition of the layers that the seismic waves went through.
 In a first phase following the measurement of these data signals, image generation algorithms, wellknown in the art, are used to reconstruct a raw picture of the underground in the form of seismic images, sometimes also referred to as echographic images. These images can be either twodimensional in shape or threedimensional. Such seismic images comprise pixels the intensity of which is correlated to a seismic wave amplitude, dependent on the local impedance variation.
 Geophysicists are used to manipulating such seismic images displaying information relating to amplitude. By merely looking at such seismic images, a geophysicist is capable of identifying areas of the underground having distinct characteristics, and use these to determine the corresponding structure of the underground.
 Automatic techniques for extracting structural information from seismic images are known. These generally involve seismic horizon reconstruction algorithms that analyze amplitude gradients in a seismic image and extract the tangent of the local dip in a direction that is transverse to that gradient. Examples of techniques used for reconstructing a seismic horizon using a seismic image are for example described in the French patent
FR 2 869 693 US 20130083973 .  Sometimes the exact depth of a layer can be known due to other data inputs or because of reliable geological information. Therefore, it is sometimes useful to define fixed related control points on a seismic image which are known to belong to a seismic horizon. It is then useful to compute a seismic horizon by implementing a seismic reconstruction algorithm with imposed conditions on a certain limited number of related control points.
 One method for reconstructing a seismic horizon with imposed conditions on a number of related control points is described in the article "Flattening with geological constraints" in Annual Meeting Expanded Abstracts, Society of Exploration Geophysicists (SEG), 2006, pp. 10531056 by J. Lomask and A. Guitton.
 The method disclosed in this article considers a global approach by solving a twodimensional nonlinear partial derivative equation relied on local dip. The partial derivative equation is solved using a GaussNewton approach by an iterative algorithm whose crucial step is the resolution of a Poisson equation. The approach is global in that it systematically computes a seismic horizon on the entire domain of the seismic image, no matter the number of related control points received as input.
 Even if it provides realistic seismic horizons, the method proposed by Lomask et al. suffers from two major drawbacks: its computational cost is often prohibitive for large data volumes, and it requires solving an iterative algorithm on the entire domain of the seismic image every time a change occurs in the number and/or position of the related control points received as input. Zhong Z et al: "Solution of Poisson and Laplace equations by quadrilateral quadrature element", internation Journal of Solids and Structures, vol. 35, no. 21, 1998, pp. 28052819 discloses a general efficient method for solving the Poisson's equation in a twodimensional domain with the use of pseudorectangles and transformation.
 The high computational cost of the horizon reconstruction algorithm implemented by Lomask is further increased by the computational means for solving the Poisson equation that forms the core step of the iterative algorithm. In general, another iterative algorithm may be used to solve the Poisson equation. The method disclosed by Lomask therefore comprises an iterative algorithm within another iterative algorithm.
 To overcome these drawbacks, an enhancement of the determination of a seismic horizon that optimizes the computational speed of the horizon reconstruction algorithm is sought.
 To achieve such an optimization and thereby overcome the drawbacks of the prior art, the invention provides a method for enhancing the determination, from a seismic image, of at least a portion of a seismic horizon in a threedimensional domain comprising axes
X ,Y ,Z . In this threedimensional domain, the seismic horizon is a function of coordinates along axesX ,Y . The method comprises:  receiving the seismic image, the seismic image having points associated with coordinates along axes
X ,Y ,Z ;  receiving a plurality of related control points associated with coordinates on axes
X ,Y ,Z ;  in a reference plane defined by axes
X andY , defining, for at least one related control point among the plurality of related control points, an associated reference point with coordinates along axesX ,Y , among a plurality of reference points, the reference point having coordinates on axesX andY identical to coordinates on axesX andY of the related control point,  defining pseudorectangles in said reference plane, each pseudorectangle comprising a reference point among a plurality of reference points.
 In a subsequent step, the invention consists in, for each current pseudorectangle among the defined pseudorectangles:
 applying a diffeomorphic transformation F to the current pseudorectangle, the diffeomorphic transformation F:
 being a function of coordinates along X, Y and defining a new domain comprising axes , , ;
 transforming points of the seismic image having coordinates along axes X, Y identical to coordinates along axes X, Y of points in the current pseudorectangle, the points of the seismic image including the related control point associated with the current pseudorectangle;
 transforming the current pseudorectangle into a corresponding rectangle;
 applying a horizon reconstruction algorithm to the transformed points, to determine a part of a transformed horizon, the part of a transformed horizon comprising the transformed related control point, the reconstruction of the seismic horizon comprising solving the Poisson equation Δ(δτ)=div(r), where δτ is an unknown function of coordinates along axes X', Y', Δ denotes the Laplace operator in the new domain, div denotes the divergence vector operator in the new domain and r is a fixed function of coordinates along axes X', Y';
 computing a part of the horizon, the computing of a part of the horizon comprising applying an inverse diffeomorphic transformation F^{1} of the diffeomorphic transformation F to the determined part of a transformed horizon.
 The term pseudorectangle is used to refer to any quadrangle or quadrilateral that has a convex shape, that is to say that each of its inner angles is smaller than 180°. Simple diffeomorphic transformations can be used to transform a convex quadrangle into a rectangle.
 Axes
X ,Y ,Z are used to define corresponding coordinates x, y and z for each point in the threedimensional domain.  For the sake of clarity, any point belonging to the reference plane will be referred to using the adjective reference, e.g. a reference center, and the corresponding points on the seismic horizon having the same x and y coordinates will be referred to using the adjective related, e.g. a related central point.
 One advantageous feature of the invention resides in the definition of pseudorectangles that delimit portions of the threedimensional domain. Each of these portions has a pseudorectangular section and comprises points in the vicinity of a related control point. A horizon reconstruction algorithm is applied to the points of these portions of the threedimensional domain. The combined volume of these portions, corresponding to the sum of all the volumes of the portions defined by pseudorectangles, may be smaller than the volume of the domain corresponding to the entire seismic image. This reduction in volume provides a first enhancement of the computational speed of the horizon reconstruction algorithm.
 A second advantageous feature of the invention is that it provides fast means for solving the Poisson equation, the latter generally implementing an iterative algorithm within the horizon reconstruction algorithm. To do so, the invention introduces for each previously defined pseudorectangle, a corresponding diffeomorphic transformation F which transforms each pseudorectangle into a corresponding rectangle in a transformed reference plane defined by axes X' and Y'. The same diffeomorphic transformation F also transforms the points of the corresponding portion of the threedimensional domain into transformed points which are within a transformed portion of the threedimensional domain delimited by the corresponding rectangle. The purpose of this transformation is to meet some conditions in which the Poisson equation can be solved in one step, i.e. using direct calculation techniques that do not rely on an iterative algorithm. It is known, by a man skilled in the art of solving Poisson equations on discrete systems, that at least two conditions can be met to enable such a fast computation:
 the portion of the threedimensional domain on which the equation is solved advantageously has a rectangular or circular section and,
 either
 at least one related control point belongs to the latter portion of the threedimensional domain, this being also associated with specific conditions on the boundaries called Neumann condition, or
 the boundary conditions along the edges of the portion of the threedimensional domain are known, the latter condition being also referred to as Dirichlet boundary condition.
 In the invention, the diffeomorphic transformation of pseudorectangles into rectangles ensures that the first condition is met. The diffeomorphic transformation associated with a pseudorectangle is applied to all the points of the seismic image whose coordinates along axes
X andY match those of points in the pseudorectangle. The coordinates along axisZ are not affected by that transformation. The second condition is met by defining pseudorectangles that comprise reference points associated with related control points received as input and the coordinates of which are known.  Another original feature of the invention resides in the fact that each diffeomorphic transformation is applied to the portion of the seismic image comprising points having the same x and y coordinates as points in the pseudorectangles. Therefore, it may not be necessary to replace the Laplace operator of the Poisson equation by a differential operator with variable coefficients, which would render the resolution of the Poisson equation complex. In the invention, the divergence operator and the fixed function r are the ones that are transformed, thereby enabling the implementation of fast solvers and not necessarily matrix methods.
 Finally, another original feature of this invention is the possibility of choosing pseudorectangles delimiting portions of the threedimensional domain having any section suitable for encompassing the received related control points. This is particularly interesting in situations where the related control points are inhomogeneously scattered in the threedimensional domain, with areas locally having higher concentrations of related control points. In such situations, defining a portion of the threedimensional domain with a rectangular section may prove difficult insofar as it may require defining rectangles with small dimensions, sometimes referred to as degenerated rectangles. Horizon reconstruction algorithms might suffer from an insufficient number of data points in portions delimited by such degenerated rectangles and provide less accurate results. The use of pseudorectangles gives more freedom in choosing shapes adapted to the local distribution of related control points without suffering from the disadvantages that arise when defining portions of the threedimensional domain delimited by rectangles.
 More specifically, it may be advantageous that a pseudorectangle is defined so that the reference point comprised in a pseudorectangle belongs to a current reference edge of said pseudorectangle.
 In this embodiment, the portion of a seismic horizon is determined by first determining the boundaries of the portion of the domain delimited by the current pseudorectangle. Having a reference point on a current reference edge may increase the efficiency of the algorithm by providing means for calculating these boundaries of the seismic horizon. Indeed, when a reference point belongs to a current reference edge of a pseudorectangle, the associated related control point belongs to a related edge of the seismic horizon. It may then be possible to implement a calculation of the boundaries on the sought seismic horizon.
 A further improvement of the method of the invention may consist in choosing advantageous methods for finding boundary conditions in the portion of the threedimensional domain delimited by a pseudorectangle comprising reference points on a current reference edge.
 To this end, prior to applying a diffeomorphic transformation F, the method may comprise applying, for each current pseudorectangle comprising a reference point belonging to a current reference edge of said pseudorectangle among the defined pseudorectangles, for each current reference edge of said current pseudorectangle, a horizon reconstruction algorithm to edge points having coordinates along axes
X ,Y identical to the coordinates along axesX ,Y of reference edge points of said current reference edge.  The horizon reconstruction algorithm implemented to compute these boundary conditions may be a simplified algorithm insofar as its solutions are functions that can be graphically represented in two dimensions as lines. A first current reference edge may advantageously be chosen as being the one comprising the reference point associated with the related control point. A first horizon line comprising said related control point and forming a first related edge associated with reference edge points of the first current reference edge may be determined. The extremities of this first related edge may be used to determine, respectively, a second and third related edge, by implementing horizon reconstruction algorithms in a similar fashion on points of faces of the portion of the threedimensional domain delimited by the current pseudorectangle associated with reference edge points of a second and third current reference edge. Two extremities of the second and third related edge may correspond to extremities of a fourth related edge. Therefore the fourth related edge may be determined by implementing a horizon reconstruction algorithm on edge points of a face associated with a fourth current reference edge, with the condition that the horizon line passes through both extremities of the fourth related edge.
 It may be advantageous to perform the calculation of the boundaries prior to applying a diffeomorphic transformation to each pseudorectangle, insofar as some pseudorectangles and therefore, the portions of the threedimensional domain that is delimited by these pseudorectangles, may share at least a portion of an edge. In this way, it may be possible to reduce the number of calculations that are performed to determine the boundary conditions by using the already calculated boundaries of portions of the threedimensional domain delimited by adjacent pseudorectangles. It may however also be possible to perform these calculations individually for each pseudorectangle in the transformed domain after applying a diffeomorphic transformation F. In this alternative embodiment of the invention, it may be possible to use the corresponding inverse diffeomorphic transformation F^{1} to reuse the portions of boundaries that are identical for the portions of the threedimensional domain delimited by two adjacent pseudorectangles.
 Some techniques for defining pseudorectangles may be particularly advantageous, may further reduce the computation time of the algorithm, and may be easy to implement.
 For instance, it may be possible to define pseudorectangles such that at least one reference corner of each pseudorectangle among the defined pseudorectangles may have coordinates along axes
X ,Y identical to the coordinates along axesX ,Y of a related control point among the plurality of related control points.  In such an embodiment, each pseudorectangle among the defined pseudorectangles may have a reference corner associated with a related control point, thus enabling an easy calculation of the boundary conditions, for example by applying successive horizon reconstruction algorithms to points of faces of the portion of the threedimensional domain comprising a reference edge comprising said corner and axis
Z .  In a particularly advantageous configuration, the received plurality of related control points may comprise at least three related control points, and defining pseudorectangles comprises:
 identifying reference points in the reference plane;
 identifying triangles having a first reference corner, a second reference corner and a third reference corner among the identified reference points using a triangulation, and
 in each of the identified triangles:
 identifying a reference centroid of said triangle,
 identifying a first reference center of the segment defined by the first reference corner and the second reference corner;
 identifying a second reference center of the segment defined by the first reference corner and the third reference corner;
 Such a method of defining pseudorectangles may provide several advantages. First of all, it can be easily implemented by a computer program, no matter the distribution of the related control points. Secondly, this method may optimize the size distribution of the pseudorectangles, since the area of the pseudorectangles that are part of a given triangle is substantially the same. Thirdly, this way of defining pseudorectangles may greatly facilitate the determination of boundary conditions, since a reference corner of each pseudorectangle is associated with a related control point, and the triangles define lines joining reference points. These lines enable an easy calculation of the corresponding horizon line by applying a horizon reconstruction algorithm to points of a plane comprising axis
Z and two of the related control points.  More specifically, when pseudorectangles are defined in this way, the method of the invention may advantageously comprise, for an identified triangle, and prior to applying a diffeomorphic transformation F:
 identifying a first, second and third related control point among the plurality of related control points associated with corresponding first , second and third reference corners of said identified triangle;
 applying a horizon reconstruction algorithm to points of a plane comprising axis
Z and comprising the first and second related control points to determine a first portion of a first local horizon;  identifying a first related central point on the first portion of the first local horizon having coordinates along axes
X andY identical to coordinates along axesX andY of the first reference center;  applying a horizon reconstruction algorithm to points of a plane comprising axis
Z and comprising the first and third related control points to determine a second portion of a second local horizon;  identifying a second related central point on the second portion of the second local horizon having coordinates along axes
X andY identical to coordinates along axesX andY of the second reference center;  computing a coordinate along axis Z of a related middle point having coordinates along axes
X andY identical to coordinates along axesX andY of the reference centroid of said identified triangle, the computation of said coordinate along axisZ being a function of the coordinates of a point on said determined first or second local horizons.  More specifically, the computation of the z coordinate of the related middle point can be a function of any point belonging to the first or second local horizon. For example, it could advantageously be a function of one of the extremities of the first or second local horizons, or either related central point.
 The method described above may benefit from one major advantage: it may be particularly efficient from a computational point of view because many steps are implemented once for a first identified triangle, but can be skipped when applying the method to points associated with adjacent triangles. This more specifically concerns the portions of local horizons joining two related control points associated with two reference corners of a triangle. These portions of local horizons may be shared by two adjacent portions of the threedimensional domain delimited by two adjacent triangles.
 It may be possible to compute the coordinates along axis
Z of the related middle point of the identified triangle by applying a horizon reconstruction algorithm to points of a plane comprising axisZ , and comprising the segment connecting the first reference center with the reference centroid or the segment connecting the second reference center with the reference centroid.  Doing so may increase the precision of the above mentioned method.
 Alternatively, computing a coordinate along axis
Z of the related middle point can also be achieved by calculating the mean value of the coordinates along axisZ of at least the first and second related central points.  This technique may be very quick and provide a good accuracy especially if the size of the triangle is small.
 Several techniques may be foreseen to solve the Poisson equation that is computed in the horizon reconstruction algorithm. Once the conditions required for a onestep direct resolution of the equation are met, it may be advantageous to solve the Poisson equation using a Fourier transform algorithm.
 The latter algorithms are wellknown and easy to implement in a computer program for instance, due to the multitude of existing libraries for performing Fourier transforms on discrete data. Furthermore, Fourier transform algorithms are excellent alternatives to matrix methods, the latter being a lot more complex to compute.
 The method described above can be implemented on portions of the threedimensional domain comprising points having the same x and y coordinates as individualized pseudorectangles.
 However it is possible to define pseudorectangles that map a continuous portion of the reference plane.
 This may increase the computational speed of the method due to the fact that some of the computed data, for example the boundaries, can be reused on portions of the threedimensional domain delimited by neighboring pseudorectangles.
 In a final step, once twodimensional portions of a horizon have been calculated for each of the defined pseudorectangles, the method may further comprise assembling all these portions of horizons to define a finalized portion of a reconstructed horizon.
 To do so, the method may comprise computing a portion of a seismic horizon from at least the computed part of the horizon of each current pseudorectangle among the defined pseudorectangles.
 When pseudorectangles were defined using a triangulation as described above, the method may further comprise computing a portion of a seismic horizon from at least the computed part of the horizon of each current pseudorectangle among the defined pseudorectangles, and after computing a portion of a seismic horizon, the method may comprise:
 receiving modification information relating to the related control points;
 identifying pseudorectangles affected by the received modification information relating to the related control points;
 defining a new set of pseudorectangles in a local area corresponding to the area occupied by the pseudorectangles affected by said received modification information relating to the related control points;
 for each current pseudorectangle among the new set of pseudo rectangles:
 applying a diffeomorphic transformation F, said diffeomorphic transformation F:
 being a function of coordinates along
X ,Y and defining a new domain comprising axesX' ,Y' ,Z ;  transforming points of the seismic image having coordinates along axes
X ,Y identical to coordinates along axesX ,Y of points in said current pseudorectangle, said points of the seismic image including the related control point associated with the current pseudorectangle;  transforming said current pseudorectangle into a corresponding rectangle;
 being a function of coordinates along
 applying a horizon reconstruction algorithm to the transformed points, to determine a part of a transformed horizon, said part of a transformed horizon comprising the transformed related control point , the reconstruction of the seismic horizon comprising solving the Poisson equation Δ(δτ)=div(r), where δτ is an unknown function of coordinates along axes
X' ,Y' , Δ denotes the Laplace operator in the new domain, div denotes the divergence vector operator in the new domain and r is a fixed function of coordinates along axesX' ,Y' ;  computing a part of the horizon, said computing of a part of the horizon comprising applying an inverse diffeomorphic transformation F^{1} to the determined part of a transformed horizon.
 applying a diffeomorphic transformation F, said diffeomorphic transformation F:
 Therefore, whenever new related control points are added, or former related control points are removed, the method can efficiently limit the portion of the threedimensional domain on which new calculations are performed to the portion of the threedimensional domain concerned by the modifications that were performed.
 The invention also pertains to a device for enhancing the determination, from a seismic image, of at least a portion of a seismic horizon in a threedimensional domain comprising axes , , , said seismic horizon being a function of coordinates along axes X, Y in said threedimensional domain, wherein said device comprises:
 an input interface for receiving the seismic image, the seismic image having points associated with coordinates along axes X, Y, Z; and for receiving a plurality of related control points associated with coordinates on axes X, Y, Z;
 a circuit for defining, in a reference plane defined by axes X and Y, for at least one related control point among the plurality of related control points, an associated reference point with coordinates along axes X, Y, among a plurality of reference points, the reference point having coordinates on axes X and Y identical to coordinates on axes X and Y of the related control point,
 a circuit for defining pseudorectangles in the reference plane, each pseudorectangle comprising a reference point among a plurality of reference points;
 a circuit being adapted for, for each current pseudorectangle among the defined pseudorectangles:
 applying a diffeomorphic transformation F to the current pseudorectangle, said diffeomorphic transformation F:
  being a function of coordinates along X, Y and defining a new domain comprising axes , , ;
  transforming points of the seismic image having coordinates along axes X, Y identical to coordinates along axes X, Y of points in said current pseudorectangle, said points of the seismic image including the related control point associated with the current pseudorectangle;
  transforming said current pseudorectangle into a corresponding rectangle;
 applying a horizon reconstruction algorithm to the transformed points, to determine a part of a transformed horizon, said part of a transformed horizon comprising the transformed related control point, the reconstruction of the seismic horizon comprising solving the Poisson equation Δ(δτ)=div(r), where δτ is an unknown function of coordinates along axes X', Y', Δ denotes the Laplace operator in the new domain, div denotes the divergence vector operator in the new domain and r is a fixed function of coordinates along axes X', Y';
 computing a part of the horizon, said computing of a part of the horizon comprising applying an inverse diffeomorphic transformation F^{1} of the diffeomorphic transformation F to the determined part of a transformed horizon.
 applying a diffeomorphic transformation F to the current pseudorectangle, said diffeomorphic transformation F:
 The invention also pertains to a nontransitory computer readable storage medium, having stored thereon a computer program comprising program instructions, the computer program being loadable into a dataprocessing unit and adapted to cause the dataprocessing unit to carry out the sequence of operations of the method described above when the computer program is run by the dataprocessing device.
 The method of the invention will be better understood by reading the detailed description of exemplary embodiments presented below. These embodiments are illustrative and by no means limitative. They are provided with the appended figures and drawings on which:

figure 1 is a schematic representation of a seismic image in a threedimensional domain; and 
figure 2 is a schematic representation of the threedimensional domain offigure 1 comprising related control points and their associated reference points in the reference plane; and 
figure 3 is a schematic representation of the reference plane offigure 2 ; and 
figure 4 is a schematic representation of a plane pointed at onfigure 3 and comprising axisZ , a portion of seismic image, a current reference edge of a pseudorectangle and a related control point associated with a reference point on the current reference edge; and 
figure 5 is a schematic representation of the threedimensional domain offigure 1 comprising one related control point the associated current pseudorectangle and the boundaries of the sought seismic horizon delimited by the current pseudorectangle; and 
figure 6 presents schematic representations (A and B) of the transformation operated by the diffeomorphic transformation F associated with the pseudorectangle offigure 5 ; and 
figure 7 presents schematic representations (A and B) of the transformation operated by the inverse diffeomorphic transformation F^{1} associated with the pseudorectangle offigures 5, 6 element A and 6 element B; and 
figure 8 is a schematic representation of the threedimensional domain offigure 1 comprising related control points and their associated portions of a reconstructed seismic horizon; and 
figure 9 is a schematic representation of the reference plane offigure 2 according to a second embodiment; and 
figure 10 is a schematic representation of the reference plane offigure 9 with three pseudorectangles defined in accordance with the second embodiment; and 
figure 11 is a schematic representation of the reference plane offigure 9 illustrating the pseudorectangles affected by the addition of a related control point; and 
figure 12 is a flow chart illustrating the main steps implemented by the horizon reconstruction method; and 
figure 13 is a possible embodiment for a device that enables the present invention.  For the sake of clarity, the dimensions of features represented on these figures may not necessarily correspond to the realsize proportions of the corresponding elements. Like reference numerals on the figures correspond to similar elements or items.

Figure 1 represents an exemplary seismic image in a threedimensional domain 1 associated with axesX ,Y ,Z . Such an image comprises dark regions 101, 102, 103 alternating with brighter regions 110, 120, 130. From the data contained in the seismic image offigure 1 , geophysicists may extract the tangent of the local dip p associated with every data point of the seismic image. The tangent of the local dip is expressed as a function of class C^{1} of x, y, z coordinates. The aim of a horizon reconstruction method is to find a twodimensional surface in the threedimensional domain 1, that can be numerically represented as a function of class C^{2}:$$\mathrm{\tau}:\left(\mathrm{x},\mathrm{y}\right)\to \mathrm{\tau}\left(\mathrm{x},\mathrm{y}\right)$$ of x, y coordinates and verifying the condition :$$\mathrm{\tau}=\underset{\mathrm{f}\in {\mathrm{C}}^{2}}{\mathrm{arg}\phantom{\rule{1ex}{0ex}}\mathrm{min}}{{\displaystyle \underset{\mathrm{\Omega}}{\int}\Vert \nabla \mathrm{f}\left(\mathrm{x},\mathrm{y}\right)\mathrm{p}\left[\mathrm{x};\mathrm{y};\mathrm{f}\left(\mathrm{x},\mathrm{y}\right)\right]\Vert}}^{2}\mathrm{d\Omega}$$ where   denotes a norm, for example the absolute value, ∇ denotes the gradient operator and Ω the portion of the threedimensional domain 1 on which the seismic horizon is calculated. Iterative horizon reconstruction algorithms to solve the above equation are wellknown from the existing prior art, such as for example from the abovecited article by Lomask et al.  In the process of implementing a horizon reconstruction algorithm, one constraint resides in the fact that any calculated horizon must pass through all the related control points received as input.
 Several key steps are implemented in such an algorithm. Generally, a first horizon corresponding to a function τ = τ_{0} is initialized. Then, a residual term r is calculated. This term r is another function of coordinates x, y, verifying the condition r(x;y) =∇τ(x;y)p[x;y;τ(x;y)], which corresponds to the difference between the tangent of the local dip of the seismic image and the gradient of the horizon.
 While implementing the iterative horizon reconstruction algorithm, the main challenge resides in minimizing this residual term r. This is done by progressively correcting function x, so that after each step k of the horizon reconstruction algorithm, τ_{k+1}= τ_{k}+δτ_{k}. At each step, an update term δτ is computed, the latter verifying:
$$\mathrm{\delta \tau}=\underset{\mathrm{f}\in {\mathrm{C}}^{2}}{\mathrm{arg}\phantom{\rule{1ex}{0ex}}\mathrm{min}}{{\displaystyle \underset{\mathrm{\Omega}}{\int}\Vert \nabla \mathrm{f}\left(\mathrm{x},\mathrm{y}\right)+\mathrm{r}\left(\mathrm{x},\mathrm{y}\right)\Vert}}^{2}\mathrm{d\Omega}$$ 
 As mentioned above, the invention resides in the way this Poisson equation is calculated.
 As illustrated on
figure 2 , the method comprises receiving related control points 201, 202, 203, 204, 205, 206, 207, 208 in the threedimensional domain 1. These related control points 201, 202, 203, 204, 205, 206, 207, 208 may for example be points that are known to belong to a given horizon because of drills realized in the ground or because of reliable geological data. The horizon reconstruction algorithm relies on using the x and y coordinates of the points of the threedimensional domain 1 as input, and calculating a corresponding coordinate along axisZ to determine a reconstructed horizon. The method of the invention involves transformations on these points, that only affect their x and y coordinates, but do not change their z coordinate. To simplify the process of defining pseudorectangles and diffeomorphic transformations that are part of this invention, reference points 210, 220, 230, 240, 250, 260, 270, 280 associated with said related control points are defined in a reference plane 10, this reference plane being defined by axesX andY . The reference points 210, 220, 230, 240, 250, 260, 270, 280 have the same x and y coordinates along axesX andY as the related control points 201, 202, 203, 204, 205, 206, 207, 208 i.e. the point 210 (respectively 220, 230, 240, 250, 260, 270, 280) is a projection of the related control point 201 (respectively 202, 203, 204, 205, 206, 207, 208) on a plane surface (X ,Y ).  As illustrated on
figure 3 , the invention then consists in defining pseudorectangles in the reference plane 10 comprising the reference points 210, 220, 230, 240, 250, 260, 270, 280 associated with related control points. This may be done in many different ways, some of which are illustrated onfigures 3 ,9 and 10 . Onfigure 3 , pseudorectangles with random shapes map a portion of the reference plane 10. Each of these pseudorectangles contains one of the reference points 210, 220, 230, 240, 250, 260, 270, 280. The latter points can be located anywhere on a current pseudorectangle. For example, reference point 280 belongs to a reference corner of a current pseudorectangle, and reference point 220 belongs to a current reference edge of a current pseudorectangle 3220.  The pseudorectangles comprising reference points 210, 220, 230, 240, 250, 260, 270, 280 verify the boundary conditions called Neumann conditions, which state that for a unique point of fixed coordinates on the horizon, the derivative of the update term along the exterior normal
ω to the boundary is assumed to be equal to zero and its mean value fixed to zero. In other words, for any value of coordinates x and y along the edges of the horizon in the portion Ω of the threedimensional domain 1 delimited by the current pseudorectangle, the following scalar product is equal to zero: ∇δτ(x;y).ω (x;y)=0. In such pseudorectangles, it is advantageous to avoid calculating the boundary conditions since these boundaries are not required to rapidly solve the horizon reconstruction algorithm. It may also be advantageous to verify that adjacent calculated portions of a seismic horizon form a continuous surface, and implement corrections to ensure that there is no discontinuity at their shared boundary.  In another embodiment, it may be advantageous to compute the boundary conditions on the edges of the horizon in the portion Ω of the threedimensional domain 1 delimited by the current pseudorectangle, to verify the Dirichlet conditions and in order to be sure that the different determined horizons for each pseudorectangle are continuous. On
figure 3 , a plane 20 defined by axisZ and containing reference point 220 and reference corners 2220, 2210 is represented. This plane 20 comprises the current reference edge 320 of the current pseudorectangle 3220. Onfigure 3 , this plane 20 appears as a line.  On
figure 4 , the same plane 20 is represented with the points from the seismic image having the same coordinates in the threedimensional domain 1 as points from the plane 20, reference point 220, the related control point 202, and the reference corners 2220, 2210. To find the related edge 302 comprising related control point 202 and belonging to the seismic horizon, a horizon reconstruction algorithm can be applied to points of plane 20. This horizon reconstruction algorithm is easier to implement since it resolves the Poisson equation in twodimensions, that is to say, it computes a function τ which can be expressed as a function of one variable and which can be graphically represented in a plane. As can be seen onfigure 4 , the reconstructed horizon line 302 tends to follow the tangent of the dip of the points from the seismic image.  The boundaries of the sought horizon are represented on
figure 5. Figure 5 represents the portion of the threedimensional domain 1 delimited by pseudorectangle 3220. This portion comprises four faces: face 501 appears on the left side, face 504 on the right side, face 502 at the back and face 503 at the front of the illustration onfigure 5 . Knowing a related edge 302, corresponding to a horizon line of the sought horizon, comprised in face 501, it is possible to compute the boundaries 420. The horizon line 302 can be used to compute the other horizon lines along the adjacent faces 502, 503 of the current portion of the threedimensional domain 1 delimited by the current pseudorectangle 3220. To do so, the extremities 2201 and 2202 of the horizon line are used in two horizon reconstruction algorithms to determine a second and third horizon lines. The second horizon line passes through extremity 2202, comprises another extremity 2203 and is comprised in face 502. The third horizon line passes through extremity 2201, comprises another extremity 2204 and is comprised in face 503. The horizon line comprised in the remaining face 504 is determined by applying a horizon reconstruction algorithm to points of the remaining face 504, so that the horizon line passes through extremities 2203 of the second and 2204 third horizon line.  This step by step approach leads to the determination of the boundary conditions in the portion Ω of the threedimensional domain 1 delimited by the current pseudorectangle, thereby fulfilling the Dirichlet boundary conditions.
Figure 5 illustrates the determined boundaries 420 in the current portion Ω of the threedimensional domain 1 delimited by the current pseudorectangle associated with related control point 202.  It is to be noted that although the above description and illustrations describe a way of determining the boundary conditions in the current portion Ω, it is possible to skip this step and proceed with the method described below. Indeed, the method of this invention is also efficient in the case where a single related control point is contained in the current portion Ω. Alternatives such as the configuration in which a related control point has the same x and y coordinates as a reference corner of the current pseudorectangle, as is the case for related control point 208, is also compatible with the invention. As long as any one of the boundary conditions is met, the method of the invention further proceeds by identifying, for a current pseudorectangle, a diffeomorphic transformation F which transforms the current pseudorectangle into a corresponding rectangle. For a current pseudorectangle, such a diffeomorphic transformation F is a function which transforms coordinates (x;y) into corresponding coordinates (x',y') so that:
$$\left[\begin{array}{c}\mathrm{x}\prime \\ \mathrm{y}\prime \end{array}\right]=\mathrm{F}\left(\mathrm{x},\mathrm{y}\right)=\left[\begin{array}{c}{\mathrm{F}}_{\mathrm{x}\prime}\left(\mathrm{x},\mathrm{y}\right)\\ {\mathrm{F}}_{\mathrm{y}\prime}\left(\mathrm{x},\mathrm{y}\right)\end{array}\right]$$ 
Figure 6 (element A) illustrates a current portion Ω of the threedimensional domain 1 delimited by the current pseudorectangle associated with related control point 202, for which the Dirichlet conditions, represented by boundaries 420, have been computed. All the points of this current portion Ω are transformed using diffeomorphic transformation F to obtain the corresponding rectangle and the new domain Ω' delimited by the corresponding rectangle illustrated onfigure 6 element B. The boundary conditions 620 in the new domain as well the transformed related control point 602 are also represented. The new domain is associated with the transformed axesX' ,Y' ,Z . In addition to transforming the current portion Ω into the new domain Ω', the method of the invention also transforms the corresponding portion of the seismic image, to obtain a set of transformed points in the new domain. The gradient field of the function τ is therefore relied on a vector field by a partial differential equation:$$\nabla \mathrm{\tau}\left(\mathrm{x}\prime ;\mathrm{y}\prime \right)=\mathrm{p}\prime \left[\mathrm{x}\prime ;\mathrm{y}\prime ;\mathrm{\tau}\left(\mathrm{x}\prime ;\mathrm{y}\prime \right)\right]$$ where p' is the tangent of the transformed local dip p. It can be expressed as:$$\mathrm{p}\prime ={\mathrm{J}}_{\mathrm{F}}^{1}\mathrm{p}$$ where${\mathrm{J}}_{\mathrm{F}}^{1}$ is the inverse of the transformation Jacobian matrix J_{F} defined by:$${\mathrm{J}}_{\mathrm{F}}=\left[\begin{array}{cc}\frac{\partial \mathrm{x}\prime}{\partial \mathrm{x}}& \frac{\partial \mathrm{y}\prime}{\partial \mathrm{x}}\\ \frac{\partial \mathrm{x}\prime}{\partial \mathrm{y}}& \frac{\partial \mathrm{y}\prime}{\partial \mathrm{y}}\end{array}\right]$$  The diffeomorphic transformation F transforming a current pseudorectangle into a corresponding rectangle is a homography defined by a 3x3 matrix H = └h_{ij}┘. This transformation is given, for any x, y coordinates in the current portion Ω by:
$$\left[\begin{array}{c}\mathrm{x}\prime \\ \mathrm{y}\prime \end{array}\right]==\left[\begin{array}{c}\frac{{\mathrm{h}}_{11}\mathrm{x}+{\mathrm{h}}_{12}\mathrm{y}+{\mathrm{h}}_{13}}{{\mathrm{h}}_{31}\mathrm{x}+{\mathrm{h}}_{32}\mathrm{y}+{\mathrm{h}}_{33}}\\ \frac{{\mathrm{h}}_{21}\mathrm{x}+{\mathrm{h}}_{22}\mathrm{y}+{\mathrm{h}}_{23}}{{\mathrm{h}}_{31}\mathrm{x}+{\mathrm{h}}_{32}\mathrm{y}+{\mathrm{h}}_{33}}\end{array}\right]$$  The four terms of the Jacobian are then defined by:
$$\frac{\partial \mathrm{x}\prime}{\partial \mathrm{x}}\left(\mathrm{x},\mathrm{y}\right)=\frac{\left({\mathrm{h}}_{11}{\mathrm{h}}_{32}{\mathrm{h}}_{31}{\mathrm{h}}_{12}\right)\mathrm{y}+{\mathrm{h}}_{11}{\mathrm{h}}_{33}{\mathrm{h}}_{31}{\mathrm{h}}_{13}}{{\left({\mathrm{h}}_{31}\mathrm{x}+{\mathrm{h}}_{32}\mathrm{y}+{\mathrm{h}}_{33}\right)}^{2}}$$ $$\frac{\partial \mathrm{y}\prime}{\partial \mathrm{x}}\left(\mathrm{x},\mathrm{y}\right)=\frac{\left({\mathrm{h}}_{21}{\mathrm{h}}_{32}{\mathrm{h}}_{31}{\mathrm{h}}_{22}\right)\mathrm{y}+{\mathrm{h}}_{21}{\mathrm{h}}_{33}{\mathrm{h}}_{32}{\mathrm{h}}_{23}}{{\left({\mathrm{h}}_{31}\mathrm{x}+{\mathrm{h}}_{32}\mathrm{y}+{\mathrm{h}}_{33}\right)}^{2}}$$ $$\frac{\partial \mathrm{x}\prime}{\partial \mathrm{y}}\left(\mathrm{x},\mathrm{y}\right)=\frac{\left({\mathrm{h}}_{12}{\mathrm{h}}_{31}{\mathrm{h}}_{31}{\mathrm{h}}_{121}\right)\mathrm{y}+{\mathrm{h}}_{21}{\mathrm{h}}_{33}{\mathrm{h}}_{31}{\mathrm{h}}_{23}}{{\left({\mathrm{h}}_{31}\mathrm{x}+{\mathrm{h}}_{32}\mathrm{y}+{\mathrm{h}}_{33}\right)}^{2}}$$ $$\frac{\partial \mathrm{y}\prime}{\partial \mathrm{y}}\left(\mathrm{x},\mathrm{y}\right)=\frac{\left({\mathrm{h}}_{22}{\mathrm{h}}_{31}{\mathrm{h}}_{32}{\mathrm{h}}_{21}\right)\mathrm{x}+{\mathrm{h}}_{22}{\mathrm{h}}_{33}{\mathrm{h}}_{31}{\mathrm{h}}_{23}}{{\left({\mathrm{h}}_{31}\mathrm{x}+{\mathrm{h}}_{32}\mathrm{y}+{\mathrm{h}}_{33}\right)}^{2}}$$  It is therefore possible to compute, for each point of the new domain, a transformed residual term r and solve the Poisson equation in the transformed domain.
 With the elements obtained so far, two conditions are met to allow a direct and onestep resolution of the Poisson equation: the domain on which a solution is searched corresponds to points having x and y coordinates identical to those of a rectangle, and either at least one related control point is within this new domain, or the boundary conditions of the solution are known.
 The determination of the update term, the solution of the Poisson equation, can be calculated using fast Fourier transform algorithms, for example by solving the equation:
$$\mathrm{\delta \tau}={\mathrm{FT}}^{1}\left[\frac{\mathrm{FT}\left[\mathrm{div}\left(\mathrm{r}\right)\right]}{\mathrm{FT}\left[\mathrm{\Delta}\right]}\right]$$ where FT denotes a Fourier transform and FT^{1} denotes an inverse Fourier transform.  Advantageously, the Fourier transform is a discrete Fourier transform, and even more advantageously a fast Fourier transform. If the size of the new domain can be expressed as a number verifying 2^{a}3^{b}5^{c}7^{d}11^{e}13^{f}, where a, b, c, d, e and f are positive integers and e+f is smaller than 1, then a particularly efficient fast Fourier transform can be implemented to further reduce the computation time of the method of the invention.
 As represented on
figure 7 element A, once the transformed part of a reconstructed horizon 7020 is obtained, the method comprises applying the inverse diffeomorphic transformation F^{1} to the transformed part of a reconstructed horizon to obtain a part of a reconstructed horizon 720, as represented onfigure 7 element B.  Finally, the invention advantageously comprises assembling all the parts of a reconstructed horizon to obtain a reconstructed horizon on a portion of the threedimensional domain 1 as represented on
figure 8 .  Besides the general method described above, the invention may advantageously benefit from substantial optimizations that allow it to be performed faster and be easily programmed to be executed with minimal input from the user.
 To this end,
figure 9 represents a method for defining pseudorectangles that have a substantially similar shape and which allows a fast and reliable calculation of the boundary conditions in each pseudorectangle.  On
figure 9 reference points 210, 220, 230, 240, 250, 260, 270, 280 associated with related control points 201, 202, 203, 204, 205, 206, 207, 208 are represented in the reference plane 10. A triangulation, advantageously a Delaunay triangulation, connecting all these reference points to form triangles is implemented. Then, as represented onfigure 10 , the center of each side of an identified triangle is selected.Figure 10 represents the triangle identified by corners corresponding to reference points 210, 220 and 230. The reference centers 223, 212 and 213 of the sides of this triangle are also used to determine the centroid 2123 of this triangle, the centroid being the point where the median lines of the triangle cross. In this manner, the obtained three pseudorectangles have substantially the same area in each triangle, and the method can systematically be implemented by a computer program.  Other advantages arise from the method of defining pseudorectangles represented on
figures 9 and 10 . The sides of each triangle are lines joining two reference points having the same x and y coordinates as related control points, and boundary conditions can be easily computed in the plane comprising axisZ and comprising two related control points by using a horizon reconstruction algorithm to obtain a horizon line. Since it may occur, as seen onfigure 9 , that several triangles share a common side, the calculation of boundary conditions may not have to be computed for each triangle in the portion of the threedimensional domain 1 delimited by a triangle. Indeed the results obtained in the portion of the threedimensional domain 1 delimited by a previously identified triangle may be reused in the portion of the threedimensional domain 1 delimited by subsequent triangles.  The centroid of each triangle, called reference centroid 2123, shares the same x and y coordinates as a related middle point of the horizon. This related middle point is shared by three portions of horizon in three adjacent portions of the threedimensional domain 1. There are several options for determining the z coordinate of that middle point of the horizon.
 It is for example possible to make realistic approximations that are likely to be valid for triangles having a small area compared to the size of the threedimensional domain 1. One of these consists in calculating the mean value of the z coordinate of related central points of the horizon, associated with reference centers 212, 223, 213 of at least two of the three sides of a current triangle. Another consists in assuming the z coordinate of that related middle point is equal to the z coordinate of any related point of the horizon associated with a reference point of the triangle, for example a reference corner 220, 230, 210 or a reference center 212, 223, 213 of a side of the triangle. Another method consists in applying a horizon reconstruction algorithm to points of the plane comprising axis
Z and comprising one of the segments connecting a reference center 212, 223, 213 of a side of the triangle, and the reference centroid 2123, to obtain a horizon line.  In an alternative embodiment, it is possible to define pseudorectangles by combining the identified triangles two by two. Two adjacent triangles are combined by removing the segment they have in common. This embodiment is advantageous in that it makes it even easier to determine the boundary conditions of the portion Ω of the threedimensional domain 1 delimited by a pseudorectangle, since every reference corner of each pseudorectangle is associated with a related control point. In this embodiment, horizon lines passing through the related control points define the boundary conditions of each pseudorectangle.
 The method of the invention nonetheless also offers another major advantage over the existing prior art. Indeed, it is very efficient for computing portions of a seismic horizon when a related control point is added to or removed from a set of related control points.

Figure 11 represents reference plane 10 containing reference points 210, 220, 230, 240, 250, 260, 270, 280 associated with related control points 201, 202, 203, 204, 205, 206, 207, 208. First, modification information relating to the related control points is received, for example the addition of a related control point. Then, the reference point 1100 in the reference plane 10 associated with the added related control point requires locally redefining pseudorectangles. Nevertheless, the effect is only local as shown onfigure 11 , on which the darkest pseudorectangles correspond to the affected area that is chosen for a recalculation of the local horizon. In general, adding a related control point only affects the pseudorectangle or pseudorectangles to which the added reference point associated with the added related control point belongs. Nevertheless, it is advantageous to identify an affected area by identifying the triangle or triangles to which the reference point belongs. This may enable defining new pseudorectangles having substantially the same size as already defined surrounding pseudorectangles. Since the pseudorectangles comprising the added reference point may share boundaries with neighboring pseudorectangles, two of which may belong to neighboring triangles, it is advantageous to include these neighboring triangles into the affected area and triangulate a new set of pseudorectangles on this affected area. Onfigure 11 , the area affected by the addition of reference point 1100 implies a new triangulation giving rise to twelve new pseudorectangles. Similar conclusions arise when a related control point is removed.  For the above reason, the invention is very efficient in terms of computation time required to determine a horizon, for example when a user decides to add several related control points in a portion of the threedimensional domain 1 which requires a finer resolution in the reconstructed horizon.

Figure 12 is a flowchart schematically illustrating the different steps that are implemented by the method of this invention.  In a first step S1, a seismic image SEISM_IMG 1 is received. The seismic image 1 can for example be received from a raw seismic data treatment program that outputs the data points in the threedimensional domain 1.
 In a second step S2, related control points CTRL._PTs 201, 202, 203, 204, 205, 206, 207, 208 are received. The x, y, z coordinates of these points are fixed and they all belong to the same horizon.
 In a subsequent step S3, pseudorectangles PSEUD._RECT. are defined, in such a way that each pseudorectangle is in a reference plane and comprises at least one reference point 210, 220, 230, 240, 250, 260, 270, 280.
 In step S4, it is possible to apply, for each pseudorectangle PSEUD._RECT. one or several horizon reconstruction algorithms to points of an edge of a portion of the threedimensional domain 1 delimited by the current pseudorectangle, in order to find the boundaries 420.
 In step S50, a diffeomorphic transformation F is identified for each pseudorectangle. An identified diffeomorphic transformation F is applied to a current pseudorectangle to transform it into a corresponding rectangle. By doing so, the method generates conditions in which solving the Poisson equation can be greatly simplified.
 Step S50 also comprises applying said transformation to the points of the seismic image having the same x and y coordinates as points of the pseudorectangle.
 The invention further comprises the horizon reconstruction algorithm per se. It starts with step S51 which comprises identifying a horizon corresponding to an initialization function τ_{k} at k=0 and proceeding iteratively as follows:
 comparing the number of iterations to a preset value N. It is assumed that the calculated horizon converges to a reliable solution typically after a few tens of iterations. In case the number of iterations is smaller than the preset value N, the method proceeds by;
 computing a residual term r_{k} using the horizon τ_{k} and the tangent of the transformed local dip p at step S54;
 applying a horizon reconstruction algorithm using Fourier transforms to solve the Poisson equation in the new domain Ω' at step S54;
 incrementing k by one digit at step S55 and returning to step S52.
 When the number k of iterations reaches the target value N, the method proceeds with step S6 by applying the inverse diffeomorphic transformation F^{1} that can transform the corresponding rectangle into the current pseudorectangle, to the computed horizon τ_{k}.
 Finally, all the portions of a reconstructed horizon obtained for each pseudorectangle can be assembled to form the portion of a reconstructed horizon represented on
figure 8 .  A comparison of the method of the invention and the global optimization method disclosed by Lomask et al. was performed on real seismic data defining a volume of 1750m by 4000m by 1600m. Complex geometries and convergent structures of the treated data resulted in an extremely noisy estimated dip, so a set of twenty seven related control points were sequentially received in critical regions corresponding for example to peaks or basins of the horizon to be reconstructed, starting from an initial set of thirteen related control points.
 The number of iterations in the horizon reconstruction algorithm to reach convergence of both methods was set to thirty. For the method of the invention, each identified triangle is subdivided in three pseudorectangles as described above. The twenty seven related control points then lead to one hundred and twenty six pseudorectangles. For the global optimization method disclosed by Lomask et al. each update term δτ computation through a direction descent approach required three hundred iterations and the algorithm had to be initialized with a function τ_{0} close to the solution. This function τ_{0} was obtained from a horizon reconstructed over the entire domain by assuming that only one particular related control point was known.
 Table 1 resumes the computation time in seconds that was measured using both methods. The time in parentheses corresponds to the time measured for the calculations dedicated to the Fourier transforms.
Table 1 Size of rectangular domain (new domain) Method of the invention Method disclosed by Lomask et al. Normal size Optimal size smallest 3.3 s (1.41 s) 2.7 s (0.561 s) 79.1 s largest 9.98 s (5.47 s) 6.43 s (2.41 s) arithmetic mean 5.82 s (2.9 s) 4.26 s (1.56 s) geometric mean 5.4 s (2.54 s) 3.78 s (1.4 s)  Table 1 shows the time required to do calculations on the portions of the threedimensional domain 1 based on the size of the domain. The column labeled normal size gives the measured time that elapsed during the implementation of the method of the invention on portions of a domain that did not have a size optimized for fast Fourier transforms. The column labeled optimal size gives the same data but measured on portions of a domain that had a size suitable for implementing a fast Fourier transform algorithm. The line labeled smallest corresponds to the smallest defined portions of domains, the line labeled largest corresponds to the largest defined portions of domains, and the arithmetic and geometric means give times calculated based on a mean value of the size of the rectangular domains. It arises from the data of table 1 that the method of the invention enables reducing the computation time by as much as thirty times when compared to global approaches like the one disclosed by Lomask et al..
 Another test was conducted to determine the time that can be saved using the method of the invention when modification instructions regarding the related control points are received. Table 2 summarizes the times in seconds measured for implementing the method of the invention when increasing the number of related control points from thirteen to twentyseven. The time in parentheses corresponds to the time measured for the calculations dedicated to the Fourier transforms. In the column labeled entire reconstruction, the measured times are substantially the same, since the volume on which the computation is implemented is the entire threedimensional domain 1. In the column labeled incremental reconstruction, the method is only applied to the portion of the threedimensional domain 1 which is affected by the addition of new related control points.
Table 2 Number of related control points Entire reconstruction Incremental reconstruction 13 3.8 s (1.4 s) 18 3.73 s (1.4 s) 0.627 s (0.219 s) 23 3.72 s (1.38 s) 0.603 s (0.233 s) 27 3.78 s (1.4 s) 0.5 s (0.184 s)  It appears from table 2 that the selective computation of portions of a horizon on only those parts that are affected by the addition or removal of related control points further enhances the computational speed of the method.

Figure 13 is a possible embodiment for a device that enables the present invention.  In this embodiment, the device 1300 comprises a computer, this computer comprising a memory 1305 to store program instructions loadable into a circuit and adapted to cause circuit 1304 to carry out the steps of the present invention when the program instructions are run by the circuit 1304.
 The memory 1305 may also store data and useful information for carrying the steps of the present invention as described above.
 The circuit 1304 may be for instance:
 a processor or a processing unit adapted to interpret instructions in a computer language, the processor or the processing unit may comprise, may be associated with or be attached to a memory comprising the instructions, or
 the association of a processor / processing unit and a memory, the processor or the processing unit adapted to interpret instructions in a computer language, the memory comprising said instructions, or
 an electronic card wherein the steps of the invention are described within silicon, or
 a programmable electronic chip such as a FPGA chip (for « FieldProgrammable Gate Array »).
 This computer comprises an input interface 1303 for the reception of data used for the above method according to the invention and an output interface 1306 for providing a stacked model.
 To ease the interaction with the computer, a screen 1301 and a keyboard 1302 may be provided and connected to the computer circuit 1304.
 The invention is not limited to the embodiments described above and may encompass equivalent embodiments.
 For example, it is possible to define non quadrangular surfaces in the reference plane. Instead of defining pseudorectangles, it may for example be possible to define surfaces for which diffeomorphic transformations, transforming these surfaces into circles, can be obtained. Indeed, a rapid resolution of the Poisson equation in a domain having a circular section, instead of a rectangular section, is possible.
 It is possible to apply the diffeomorphic transformation F to a current pseudorectangle before calculating boundary conditions associated with the current pseudorectangle.
 It is also possible to define some pseudorectangles which are not associated with any related control point. Although doing so might seem less advantageous from a computational point of view, it may be interesting in the case in which large gaps exist between local concentrations of related control points. Defining pseudorectangles that are not associated with any related control point may allow mapping a continuous portion of the threedimensional domain 1 without having a high dispersion in the size of the pseudorectangles. It is also possible to have pseudorectangles that are not associated with any related control point, but which are adjacent to other pseudorectangles which are. Thereby, it is possible to use the boundary conditions of the neighboring pseudorectangles to meet the conditions enabling a direct resolution of the Poisson equation.
 The method described above may also be implemented in a domain comprising more than three dimensions.
 One may also define quadrangles that are not pseudorectangles, although this may render the calculation of the diffeomorphic transformations more complicated.
Claims (14)
 Method for enhancing the determination, from a seismic image, of at least a portion of a seismic horizon in a threedimensional domain (1) comprising axes
X ,Y ,Z , said seismic horizon being a function of coordinates along axesX ,Y in said threedimensional domain (1),
wherein said method comprises: receiving (S1) the seismic image, the seismic image having points associated with coordinates along axesX ,Y ,Z ; receiving (S2) a plurality of related control points (201, 202, 203, 204, 205, 206, 207, 208) associated with coordinates on axesX ,Y ,Z ; in a reference plane (10) defined by axesX andY , defining, for at least one related control point among the plurality of related control points (201, 202, 203, 204, 205, 206, 207, 208), an associated reference point with coordinates along axesX ,Y , among a plurality of reference points (210, 220, 230, 240, 250, 260, 270, 280), the reference point having coordinates on axesX andY identical to coordinates on axesX andY of the related control point, defining (S3) pseudorectangles in said reference plane (10), each pseudorectangle comprising a reference point among a plurality of reference points (210, 220, 230, 240, 250, 260, 270, 280); for each current pseudorectangle among the defined pseudorectangles: applying a diffeomorphic transformation F (S50) to the current pseudorectangle, said diffeomorphic transformation F: being a function of coordinates alongX ,Y and defining a new domain comprising axesX' ,Y' ,Z ; transforming points of the seismic image having coordinates along axesX ,Y identical to coordinates along axesX ,Y of points in said current pseudorectangle, said points of the seismic image including the related control point associated with the current pseudorectangle; transforming said current pseudorectangle into a corresponding rectangle; applying (S52, S53, S54, S55) a horizon reconstruction algorithm to the transformed points, to determine a part of a transformed horizon (7020), said part of a transformed horizon (7020) comprising the transformed related control point (602), the reconstruction of the seismic horizon comprising solving (S54) the Poisson equation ∇(δτ) = div(r), where δτ is an unknown function of coordinates along axesX' ,Y' , Δ denotes the Laplace operator in the new domain, div denotes the divergence vector operator in the new domain and r is a fixed function of coordinates along axesX' , Y'; computing a part of the horizon (720), said computing of a part of the horizon (720) comprising applying (S6) an inverse diffeomorphic transformation F^{1} of the diffeomorphic transformation F to the determined part of a transformed horizon (7020).  Method according to claim 1 wherein, a pseudorectangle is defined so that the reference point (220) comprised in a pseudorectangle (3220) belongs to a current reference edge (320) of said pseudorectangle (3220).
 Method according to claim 2, wherein prior to applying a diffeomorphic transformation F (S50), said method comprises applying, for each current pseudorectangle (3220) comprising a reference point (220) belonging to a current reference edge (320) of said pseudorectangle (3220) among the defined pseudorectangles, for each current reference edge of said current pseudorectangle (3220), a horizon reconstruction algorithm to edge points having coordinates along axes
X ,Y identical to the coordinates along axesX ,Y of reference edge points of said current reference edge.  Method according any one of the preceding claims wherein at least one reference corner of each pseudorectangle among the defined pseudorectangles has coordinates along axes
X, Y identical to the coordinates along axesX ,Y of a related control point among the plurality of related control points (201, 202, 203, 204, 205, 206, 207, 208).  Method according to any one of the preceding claims, wherein the received plurality of related control points (201, 202, 203, 204, 205, 206, 207, 208) comprises at least three related control points (201, 202, 203), and wherein defining pseudorectangles comprises: identifying reference points in a reference plane (10); identifying triangles having a first reference corner (210), a second reference corner (220) and a third reference corner (230) among the identified reference points (210, 220, 230, 240, 250, 260, 270, 280) using a triangulation, and in each of the identified triangles:wherein a pseudo rectangle is defined by segments connecting the first reference corner (210) with the first reference center (212), the first reference center (212) with the reference centroid (2123), the reference centroid (2123) with the second reference center (213) and the second reference center (213) with the first reference corner (210). identifying a reference centroid (2123) of said triangle, identifying a first reference center (212) of the segment defined by the first reference corner (210) and the second reference corner (220); identifying a second reference center (213) of the segment defined by the first reference corner (210) and the third reference corner (230);
 Method according to claim 5, wherein prior to applying a diffeomorphic transformation F (S50), the method comprises, for an identified triangle: identifying a first (201), second (202) and third (203) related control point among the plurality of related control points associated with corresponding first (210), second (220) and third (230) reference corners of said identified triangle; applying a horizon reconstruction algorithm to points of a plane comprising axis
Z and comprising the first (201) and second (202) related control points to determine a first portion of a first local horizon; identifying a first related central point on the first portion of the first local horizon having coordinates along axesX andY identical to coordinates along axesX andY of the first reference center (212); applying a horizon reconstruction algorithm to points of a plane comprising axisZ and comprising the first (201) and third (203) related control points to determine a second portion of a second local horizon; identifying a second related central point on the second portion of the second local horizon having coordinates along axesX andY identical to coordinates along axesX andY of the second reference center (213); computing a coordinate along axisZ of a related middle point having coordinates along axesX andY identical to coordinates along axesX andY of the reference centroid (2123) of said identified triangle, the computation of said coordinate along axisZ being a function of the coordinates of a point on said determined first or second local horizons.  Method according to claim 6, wherein computing a coordinate along axis
Z of the related middle point of said identified triangle is achieved by applying a horizon reconstruction algorithm to points of a plane comprising axisZ and comprising the segment connecting the first (212) reference center with the reference centroid or the segment connecting the second (213) reference center with the reference centroid (2123).  Method according to claim 6, wherein computing a coordinate along axis
Z of the related middle point is achieved by calculating the mean value of the coordinates along axisZ of at least the first (212) and second (213) related central points.  Method according to any one of the preceding claims, wherein the Poisson equation is solved (S54) using a Fourier transform algorithm.
 Method according to any one of the preceding claims, wherein the defined pseudorectangles map a continuous portion of the reference plane (10).
 Method according to any one of the preceding claims, wherein the method further comprises computing a portion of a seismic horizon (800) from at least the computed part of the horizon (720) of each current pseudorectangle among the defined pseudorectangles.
 Method according to any one of claims 5 to 8 and any of claims 9 to 10, wherein the method further comprises computing a portion of a seismic horizon (800) from at least the computed part of the horizon (720) of each current pseudorectangle among the defined pseudorectangles, and after computing a portion of a seismic horizon (800), the method comprises: receiving modification information relating to the related control points; identifying pseudorectangles affected by said received modification information relating to the related control points; defining a new set of pseudorectangles in a local area corresponding to the area occupied by said pseudorectangles affected by said received modification information relating to the related control points ; for each current pseudorectangle among the new set of pseudorectangles: applying a diffeomorphic transformation F (S50) to the current pseudorectangle, said diffeomorphic transformation F: being a function of coordinates along
X, Y and defining a new domain comprising axesX' ,Y' ,Z ; transforming points of the seismic image having coordinates along axesX ,Y identical to coordinates along axesX ,Y of points in said current pseudorectangle, said points of the seismic image including the related control point associated with the current pseudorectangle; transforming said current pseudorectangle into a corresponding rectangle; applying (S52, S53, S54, S55) a horizon reconstruction algorithm to the transformed points, to determine a part of a transformed horizon (7020), said part of a transformed horizon (7020) comprising the transformed related control point (602), the reconstruction of the seismic horizon comprising solving (S54) the Poisson equation Δ(δτ) = div(r), where δτ is an unknown function of coordinates along axesX' ,Y' , Δ denotes the Laplace operator in the new domain, div denotes the divergence vector operator in the new domain and r is a fixed function of coordinates along axesX' , Y'; computing a part of the horizon (720), said computing of a part of the horizon (720) comprising applying (S6) an inverse diffeomorphic transformation F^{1} of the diffeomorphic transformation F to the determined part of a transformed horizon (7020).  Device (1300) for enhancing the determination, from a seismic image, of at least a portion of a seismic horizon in a threedimensional domain (1) comprising axes
X ,Y, Z , said seismic horizon being a function of coordinates along axesX ,Y in said threedimensional domain (1),
wherein said device (1300) comprises: an input interface (1303) for receiving (S1) the seismic image, the seismic image having points associated with coordinates along axesX ,Y ,Z ; and for receiving (S2) a plurality of related control points (201, 202, 203, 204, 205, 206, 207, 208) associated with coordinates on axesX ,Y ,Z ; a circuit (1304) for defining, in a reference plane (10) defined by axesX andY , for at least one related control point among the plurality of related control points (201, 202, 203, 204, 205, 206, 207, 208), an associated reference point with coordinates along axesX ,Y , among a plurality of reference points (210, 220, 230, 240, 250, 260, 270, 280), the reference point having coordinates on axesX andY identical to coordinates on axesX andY of the related control point, a circuit (1304) for defining pseudorectangles in the reference plane (10), each pseudorectangle comprising a reference point among a plurality of reference points (210, 220, 230, 240, 250, 260, 270, 280); a circuit (1304) being adapted for, for each current pseudorectangle among the defined pseudorectangles: applying a diffeomorphic transformation F (S50) to the current pseudorectangle, said diffeomorphic transformation F: being a function of coordinates alongX ,Y and defining a new domain comprising axesX' ,Y' ,Z ; transforming points of the seismic image having coordinates along axesX ,Y identical to coordinates along axesX ,Y of points in said current pseudorectangle, said points of the seismic image including the related control point associated with the current pseudorectangle; transforming said current pseudorectangle into a corresponding rectangle; applying (S52, S53, S54, S55) a horizon reconstruction algorithm to the transformed points, to determine a part of a transformed horizon (7020), said part of a transformed horizon (7020) comprising the transformed related control point (602), the reconstruction of the seismic horizon comprising solving (S54) the Poisson equation Δ(δτ) = div(r), where δτ is an unknown function of coordinates along axesX' ,Y' , Δ denotes the Laplace operator in the new domain, div denotes the divergence vector operator in the new domain and r is a fixed function of coordinates along axesX' ,Y '; computing a part of the horizon (720), said computing of a part of the horizon (720) comprising applying (S6) an inverse diffeomorphic transformation F^{1} of the diffeomorphic transformation F to the determined part of a transformed horizon (7020).  A nontransitory computer readable storage medium, having stored thereon a computer program comprising program instructions, the computer program being loadable into a dataprocessing unit and adapted to cause the dataprocessing unit to carry out the steps of any of claims 1 to 12 when the computer program is run by the dataprocessing device.
Priority Applications (2)
Application Number  Priority Date  Filing Date  Title 

US201261681005P true  20120808  20120808  
PCT/EP2013/066492 WO2014023737A2 (en)  20120808  20130806  Method for enhancing the determination of a seismic horizon 
Publications (2)
Publication Number  Publication Date 

EP2883087A2 EP2883087A2 (en)  20150617 
EP2883087B1 true EP2883087B1 (en)  20210623 
Family
ID=48916099
Family Applications (1)
Application Number  Title  Priority Date  Filing Date 

EP13745135.7A Active EP2883087B1 (en)  20120808  20130806  Method for enhancing the determination of a seismic horizon 
Country Status (4)
Country  Link 

US (1)  US9214041B2 (en) 
EP (1)  EP2883087B1 (en) 
JP (1)  JP6019230B2 (en) 
WO (1)  WO2014023737A2 (en) 
Families Citing this family (4)
Publication number  Priority date  Publication date  Assignee  Title 

US9442206B2 (en) *  20121108  20160913  Total Sa  Method and device to process a three dimensional seismic image 
NO342146B1 (en) *  20160518  20180403  Roxar Software Solutions As  Method of constructing a geologic model 
US10571585B2 (en) *  20160831  20200225  Chevron U.S.A. Inc.  System and method for timelapsing seismic imaging 
US10914852B2 (en)  20170316  20210209  International Business Machines Corporation  Unsupervised identification of seismic horizons using swarms of cooperating agents 
Family Cites Families (10)
Publication number  Priority date  Publication date  Assignee  Title 

FR2652180B1 (en) *  19890920  19911227  Mallet Jean Laurent  Method for modeling a surface and device for implementing same. 
US7127100B2 (en) *  20010625  20061024  National Instruments Corporation  System and method for analyzing an image 
US6615158B2 (en) *  20010625  20030902  National Instruments Corporation  System and method for analyzing a surface by mapping sample points onto the surface and sampling the surface at the mapped points 
FR2869693B1 (en)  20040430  20060602  Total France Sa  METHOD AND PROGRAM FOR PROPAGATION OF A SEISMIC MARKER IN A SET OF SEISMIC TRACES 
US20060047429A1 (en) *  20040824  20060302  Adams Steven L  Method of estimating geological formation depths by converting interpreted seismic horizons from the time domain to the depth domain 
FR2923312B1 (en) *  20071106  20091218  Total Sa  METHOD FOR PROCESSING SEISMIC IMAGES OF THE BASEMENT 
EP2281212B1 (en) *  20080522  20190227  Exxonmobil Upstream Research Company  Seismic horizon skeletonization 
AU2009344283B2 (en) *  20090416  20130912  Landmark Graphics Corporation  Seismic imaging systems and methods employing a fast targetoriented illumination calculation 
US8694262B2 (en) *  20110815  20140408  Chevron U.S.A. Inc.  System and method for subsurface characterization including uncertainty estimation 
FR2980854A1 (en)  20111004  20130405  Total Sa  METHOD FOR SCORING SEISMIC HORIZONS DISCONTINUOUS IN SEISMIC IMAGES 

2013
 20130806 EP EP13745135.7A patent/EP2883087B1/en active Active
 20130806 US US14/420,267 patent/US9214041B2/en active Active
 20130806 WO PCT/EP2013/066492 patent/WO2014023737A2/en active Application Filing
 20130806 JP JP2015525866A patent/JP6019230B2/en active Active
Also Published As
Publication number  Publication date 

US20150199845A1 (en)  20150716 
JP6019230B2 (en)  20161102 
WO2014023737A2 (en)  20140213 
JP2015529814A (en)  20151008 
US9214041B2 (en)  20151215 
WO2014023737A3 (en)  20140501 
EP2883087A2 (en)  20150617 
Similar Documents
Publication  Publication Date  Title 

Yang et al.  Application of optimal transport and the quadratic Wasserstein metric to fullwaveform inversion  
EP2883087B1 (en)  Method for enhancing the determination of a seismic horizon  
Leung et al.  Eulerian Gaussian beams for highfrequency wave propagation  
Béréziat et al.  Solving illposed image processing problems using data assimilation  
Riyanti et al.  A new iterative solver for the timeharmonic wave equation  
EP2271952A2 (en)  Visulation of geologic features using data representations thereof  
US9341729B2 (en)  Amplitude contrast seismic attribute  
Luo et al.  Higherorder schemes for 3D firstarrival traveltimes and amplitudes  
Merland et al.  Voronoi grids conforming to 3D structural features  
Bois et al.  A fully optimal anisotropic mesh adaptation method based on a hierarchical error estimator  
CN107407736B (en)  Generate the multistage full wave field inversion processing of the data set without multiple wave  
US9542507B2 (en)  Feature detection in seismic volumes  
Martinez et al.  Denoising of gravity gradient data using an equivalent source technique  
EP3211594B1 (en)  Seismic modeling system providing seismic survey data inpainting based upon suspect region boundary comparisons and related methods  
Zinck et al.  Fast seismic horizon reconstruction based on local dip transformation  
AU2016220145B2 (en)  Black hole boundary conditions  
Motta et al.  A 3D sketchbased formulation to model salt bodies from seismic data  
Martínez et al.  Skeletal representations of orthogonal shapes  
Hjelle et al.  A numerical framework for modeling folds in structural geology  
Chen et al.  Constructing misfit function for full waveform inversion based on sliced Wasserstein distance  
Rosman et al.  Efficient Beltrami image filtering via vector extrapolation methods  
OA17254A (en)  Method for enhancing the determination of a seismic horizon.  
Jain et al.  A comparative study of iterative solvers for image denoising  
Noack et al.  Fast computation of eikonal and transport equations on graphics processing units computer architectures  
Doghraji et al.  Seismic horizon reconstruction on polygonal domains using the SchwarzChristoffel transformation 
Legal Events
Date  Code  Title  Description 

PUAI  Public reference made under article 153(3) epc to a published international application that has entered the european phase 
Free format text: ORIGINAL CODE: 0009012 

AK  Designated contracting states 
Kind code of ref document: A2 Designated state(s): AL AT BE BG CH CY CZ DE DK EE ES FI FR GB GR HR HU IE IS IT LI LT LU LV MC MK MT NL NO PL PT RO RS SE SI SK SM TR 

17P  Request for examination filed 
Effective date: 20150211 

AX  Request for extension of the european patent 
Extension state: BA ME 

DAX  Request for extension of the european patent (deleted)  
RAP1  Party data changed (applicant data changed or rights of an application transferred) 
Owner name: UNIVERSITE DE BORDEAUX Owner name: ECOLE NATIONALE SUPERIEURE DES SCIENCES AGRONOMIQUES DE BORDEAUXAQUITAINE Owner name: INSTITUT POLYTECHNIQUE DE BORDEAUX Owner name: CENTRE NATIONAL DE LA RECHERCHE SCIENTIFIQUE (C.N.R.S.) Owner name: TOTAL SE 

STAA  Information on the status of an ep patent application or granted ep patent 
Free format text: STATUS: GRANT OF PATENT IS INTENDED 

GRAP  Despatch of communication of intention to grant a patent 
Free format text: ORIGINAL CODE: EPIDOSNIGR1 

INTG  Intention to grant announced 
Effective date: 20210120 

GRAS  Grant fee paid 
Free format text: ORIGINAL CODE: EPIDOSNIGR3 

GRAA  (expected) grant 
Free format text: ORIGINAL CODE: 0009210 

STAA  Information on the status of an ep patent application or granted ep patent 
Free format text: STATUS: THE PATENT HAS BEEN GRANTED 

AK  Designated contracting states 
Kind code of ref document: B1 Designated state(s): AL AT BE BG CH CY CZ DE DK EE ES FI FR GB GR HR HU IE IS IT LI LT LU LV MC MK MT NL NO PL PT RO RS SE SI SK SM TR 

REG  Reference to a national code 
Ref country code: GB Ref legal event code: FG4D 

REG  Reference to a national code 
Ref country code: CH Ref legal event code: EP 

REG  Reference to a national code 
Ref country code: AT Ref legal event code: REF Ref document number: 1404772 Country of ref document: AT Kind code of ref document: T Effective date: 20210715 Ref country code: DE Ref legal event code: R096 Ref document number: 602013078059 Country of ref document: DE 

REG  Reference to a national code 
Ref country code: IE Ref legal event code: FG4D 

REG  Reference to a national code 
Ref country code: NL Ref legal event code: FP 

PGFP  Annual fee paid to national office [announced via postgrant information from national office to epo] 
Ref country code: NL Payment date: 20210819 Year of fee payment: 9 

REG  Reference to a national code 
Ref country code: LT Ref legal event code: MG9D 

REG  Reference to a national code 
Ref country code: NO Ref legal event code: T2 Effective date: 20210623 