Numerical Integration for Kirchhoff Migration

Résumé — Intégration numérique pour la migration de Kirchhoff — Dans les acquisitions sis-miques 3D, les points milieux d’une collection à déport constant sont distribués de façon irrégulière. Dans une telle situation, un simple «diffraction stack» ne peut réaliser correctement la migration de Kirchhoff car cette technique ne fournit pas une approximation consistante de l’intégrale de Kirchhoff. La solution à ce problème réside dans l’utilisation de véritables formules d’intégration numérique telles que celles basées sur des interpolations polynomiales de type Lagrange ou Hermite. Une intégration numérique basée sur l’interpolation de Lagrange peut parfaitemeent prendre en compte le caractère irrégulier de la distribution spatiale des points milieux, pour peu que ceux-ci soient sufﬁsamment rapprochés. L’interpolation d’Hermite, bien qu’en théorie plus précise que la précédente, ne fournit qu’une approximation très relative de la qualité des images migrées : cette amélioration n’est en fait perceptible qu’au voisinage de l’événement migré. Les deux approches sont très faciles à implémenter : elles reviennent à faire une compensation adéquate (qui dépend du schéma d’interpolation retenu) de l’amplitude des différentes traces sismiques. La mise en place de ce simple prétraitement permet Abstract — Numerical integration for Kirchhoff migration — In 3D seismic surveys, common offset data often involve an irregular distribution of midpoints. In such a situation, common offset Kirchhoff migration cannot be correctly performed by means of a mere discrete diffraction stack algorithm. Such an algorithm indeed corresponds to an inconsistent numerical integration formula. To overcome this difﬁculty, genuine numerical integration formulas (yielding a consistent approximation of the continuous diffraction stack) have to be used, e.g. numerical integration formulas based on polynomials leading to the so-called Lagrange or Hermite interpolations. A numerical integration formula based on Lagrange interpolation can cope with irregularly sampled midpoints provided that the density of midpoints involved in the common offset gather is sufﬁcient. Besides, Hermite interpolation, more accurate in theory than the former, also provides relative improvement in the images at the vicinity of the migrated event,. Both techniques can be implemented by means of a simple preprocessing (adequate scaling) of the data. Thus they are quite easy to implement in any existing diffraction stack code. In addition, they can be used in combination with ﬁlters to prevent aliasing of the migration operator. The additional computation cost is negligible compared with the cost of running the diffraction stack itself.


INTRODUCTION
Seismic imaging of complex geologic structures is still a challenging problem today.These structures are characterized by complicated interface geometries giving rise to strong lateral velocity variations (like those caused by salt bodies in the North Sea and in the Gulf of Mexico).In such cases, classical post-stack imaging techniques do not work.3D prestack depth migration is recognized as the ideal technique provided an accurate velocity model has been determined.Migration velocity analysis (MVA) has accordingly aroused considerable interest for this purpose.This technique relies on a sound interpretation of 3D common offset migrated data.
In 3D seismic surveys, offset is a vector (characterized by a norm and an azimuth) and irregularities in acquisition (due for instance to streamer feathering in marine acquisition) are frequent.These irregularities give rise to offset variations (in norm and/or in azimuth) and to irregularly sampled data (i.e.nonuniform distribution of midpoints 1 ).The influence of offset variations on migrated images has been investigated by [1].This paper examines the influence of the second kind of irregularities, i.e. the influence of a nonuniform distribution of midpoints, on the quality of the migrated image and solutions are proposed to overcome the difficulties met when imaging such data.
For the study we use a simple kinematic Kirchhoff migration algorithm (a mere diffraction stack) without taking (1) In this paper, midpoint means the true midpoint, not to be confused with the center of the associated bin.
account of dynamic effects, i.e. without attempting to preserve seismic reflection amplitudes.The conclusions we draw and the techniques proposed are nonetheless valid and suitable for a true-amplitude migration.Using a test example, we first illustrate the difficulties encountered when imaging irregularly sampled data with a mere diffraction stack.We then analyze the underlying difficulty and show at this stage that a consistent numerical integration formula is required to obtain a correct image.For the sake of simplicity, we first examine the problem of numerical integration in 1D (for application to 2D migration) and then go on to 2D (for application to 3D migration), which is precisely the situation that interests us.

3D KIRCHHOFF MIGRATION: A TEST EXAMPLE
Irregularities in acquisition cause offset variations (in norm and/or azimuth) and an irregular distribution of midpoints.It has been shown [1] that mixing offsets has no drastic influence on the quality of migrated images if the offsets vary slowly with the midpoint coordinate.However, in the case of a non uniform distribution of midpoints, severe artefacts appear, as in Figures 2 and 3, which show common offset migrated synthetic data 2 .For this experiment, the 3D model shown in Figure 1 was used: two layers are separated by an anticline reflector exhibiting a relatively high curvature (2) The offset variations involved in these data being small (the data come from a marine acquisition) we call a "common offset migration", a migration technique that actually does account for the actual locations of sources and receivers (and hence the offset variations).Small offset (0.2 km) migrated synthetic data associated with an actual marine acquisition geometry: -constant Y sections (top); -constant X sections (bottom).
In the depth migrated data volume obtained with the exact migration velocity (v = 3.0 km/s).The amplitude of the imaging artefacts is ten times lower than that of the migrated event corresponding to the reflector of the model.
and dips up to 45 • .Above the interface the velocity is constant v = 3.0 km/s.The acquisition survey (sourcereceiver positions) used to compute the synthetic data comes from an actual marine acquisition in the North Sea.The synthetic data are generated by convolving travel-times from ray tracing [2] with a Ricker wavelet centered on 25 Hz.Our imaging technique is a mere diffraction stack.In the sequel we consider two subsets of data: one associated with small offsets (around 0.2 km) and the other with large offsets (around 3 km).Figures 2 and 3 show, for the two considered offsets, some in-line and cross-line sections extracted from the depth migrated data volumes obtained with the exact velocity model: as expected, they show the migrated event that peaks on the reflector in the model, but they also show severe migration artefacts that strongly degrade the images (these artefacts are not visible at depths less than 2 km above the reflector because, in order to save computing time, the image was not calculated above this limit; this will be the case throughout this paper).Note that the same difficulty appears for small and large offsets.These artefacts are strong enough to harm a migration velocity analysis process.
This paper proposes a way to remove these artefacts.To understand the difficulties met in imaging, we first present a review of some theoretical aspects involved in Kirchhoff migration.

Kirchhoff Migration: the Continuous Diffraction Stack
For the sake of simplicity, we consider here an offset that does not vary with the midpoint coordinate.The common offset  In the depth migrated data volume obtained with the exact migration velocity (v = 3.0 km/s).The amplitude of the imaging artefacts is ten times lower than that of the migrated event corresponding to the reflector of the model.migrated image m v h at a subsurface point M for an offset h = r − s (where s and r characterize the locations of the source S and the receiver R, respectively) and a velocity field v, is defined as the superposition of elementary migrated images m v el , each of those elementary images being characterized by a specific midpoint with location defined by vector q = s + r 2 : The value of the elementary image m v el at point M is obtained as: where d( q, h, t) is the seismic trace associated with midpoint q and offset h.t v (M; q, h) is the travel-time from source S to subsurface point M to receiver R, for the migration velocity field considered v.
Considering seismic data with a single event localized in the vicinity of the arrival time function t e ( q, h), the use of the stationary phase method ( [3], for instance) yields the location of the support of the common offset migrated event.
It is the set of points M solution of the system: Note that this result relies on different assumptions: high frequency asymptotics, smooth kinematics in the data and a smooth migration velocity field.
The first Equation of system (3), namely the phase condition, gives the geometry of an elementary migrated event: it is localized in the vicinity of the isochron surface, which is the set of possible reflection points M in the subsurface.The second Equation of system (3), namely the stationarity condition, gives the geometry of the migrated event: it is the envelope of the family of isochron surfaces, parameterized by q.Thus, common offset migrated data show reflections imaged as migrated events with the geometry described above regardless of the migration velocity field.This means that even if data are migrated with an erroneous migration velocity field, we can expect the migrated data to show organized events rather than a fuzzy image.These organized events are called migrated events: they show a deformed geometry of the associated reflectors, the image of the geologic interface being all the less deformed as the migration velocity is right.

Numerical Diffraction Stack
A widely used discrete diffraction stack involves a straightforward summation of the elementary migrated images (2): We now examine whether the migrated data obtained with this discrete formula can approximate the data that would be obtained by using the continuous formula (1).For the sake of simplicity, we consider a 2D migration, thus involving, for given M and h, a 1D integration in variable q.Let us consider the function to be integrated, m v el (q), having the shape described in Figure 4.If the midpoint sampling is both regular and fine (small and constant sampling interval ∆q), discrete formula ( 4) is (up to the constant factor ∆q) a consistent approximation of the continuous diffraction stack: it is the classical rectangle formula for numerical integration of continuous functions.However, if the midpoint sampling is irregular, i.e. if the midpoint distribution is not uniform, the discrete sum ( 4) is unlikely to approximate the integral of the continuous function.From the migration standpoint, this means that using the discrete sum (4) yields a mere superposition of elementary migrated events that can hardly be compared with the continuous diffraction stack: the constructive interferences between elementary migrated events that are expected to build the migrated event (the event that peaks on the envelope) only take place to a partial extent so that the migrated event shows up among a multitude of migration smiles (the elementary migrated events).This is observable in Figures 2 and 3.The important point at this stage is to realize that the difficulty lies in the use of the mere diffraction stack (4) in conjunction with a non uniform distribution of midpoints, and this difficulty persists even if the distance is very short between neighboring midpoints.This is illustrated by the following experiment.
Let us again consider the same acquisition survey as before, but with midpoint sampling refined by adding synthetic data to the original acquisition.The technique used for this refinement is explained in detail in the Section 5 (Fig. 22) but at this stage we only need to know that the midpoint distribution is still not uniform and that the mean distance between neighboring midpoints is two times shorter than before.Figure 5 shows the corresponding migrated data obtained for the large offset (3 km) and the exact velocity model: it still shows a multitude of migration smiles and the image quality (in terms of interpretability) is not improved at all as compared with the image shown in Figure 3, which confirms our hypothesis.
This paper proposes an imaging technique that can cope with non uniformly distributed midpoints.To do this, we Sketch of the function m v el to be integrated.Left: the function is sampled with a constant spacing between discretization points q i ; in this situation, the sum of the samples approximates (up to the factor ∆q) the integral of the continuous function (rectangle formula).Right: the function is sampled with a nonconstant spacing between discretization points q i ; in this situation, the sum of the samples can hardly be compared with the integral of the continuous function.have to make use of a consistent (meaning that if the distances between neighboring midpoints go to zero, the result provided by the numerical integration formula converges to the integral of the continuous function) and, if possible, accurate numerical integration formula.The derivation and testing of these formulas are the goals of the subsequent sections.However, before opening this subject, we want to state what can and cannot be expected from a consistent numerical integration formula.In the case of an exact migration velocity 3 : -We expect that, when subsurface point M is below the reflector, the numerical integration formula yields a zero value for the migrated image.This condition is, in fact, almost always fulfilled as elementary migrated images are in general zero below the reflector (however, some (3) These expectations also hold in the case of an erroneous migration velocity provided that we replace "reflector" by "migrated event".
elementary migrated events could cross a strongly tilted envelope).-We expect that, when subsurface point M is on the reflector, even though we do not want to preserve amplitudes (but to approximate the continuous diffraction stack), the numerical integration formula gives an indication of the reflectivity for the considered offset.Looking at the values taken by the different elementary migrated images for such a point M (Fig. 6), we find that the function to be integrated is smooth and that our sampled data give detailed information on this function (note, however, that if the midpoint sampling were too coarse (holes in the acquisition), the samples would not correctly represent the continuous function.This is a hopeless situation: no numerical integration formula can create missing information).All we can expect from a numerical integration formula is that it gives a good approximation of the integral of the continuous function.From the stationary phase approximation, we know that the only samples that con- Values of function m v el at regularly spaced midpoints q with a 50 m spacing both in the x and y directions, for a subsurface point located on the interface.Left: 3D view of the function.Right: section obtained for q x = 2.6 km.Values of by function m v el at regularly spaced midpoints q with a 50 m spacing both in the x and y directions, for a subsurface point located 2.5 km above the interface.Left: 3D view of the function.Right: section obtained for q x = 4.0 km.tribute to the reflectivity evaluation are those in the neighborhood of the point where the phase is stationary (the point with highest positive amplitude) and that the contribution of all the other samples must cancel each other out.-However, looking at the values of the different elementary migrated images for a point M located 2.5 km above the reflector (Fig. 7), we cannot expect an accurate result from the best numerical integration formula: the function to be integrated is not sampled finely enough regarding its very rapid variations.In such a situation, the migration operator turns out to be aliased and, the only alternative to another survey, is to apply filters, as described e.g. by [4] or [5], which in fact make the function to be integrated vary more slowly.

The Underlying Difficulty: Conclusion
This analysis shows that if midpoints are irregularly sampled in an acquisition, Kirchhoff migration cannot be correctly performed by using a mere discrete diffraction stack.A consistent numerical integration is needed to approximate the continuous diffraction stack.Besides, at points M far above the reflector, the function to be integrated varies very rapidly with midpoint coordinate.In this situation, midpoint sampling is usually not fine enough to correctly represent the function to be integrated, thereby entailing the use of appropriate filters to prevent aliasing in the migration operator.The use of these filters in combination with a consistent and, if possible, accurate numerical integration formula is the key for obtaining good migrated images.Designing such filters has been the object of different studies (see for instance [4] or [5]).We ignore this problem in the sequel to focus rather on the design of consistent and accurate numerical integration formulas for use in combination with such filters.We describe the principle of the derivation of such formulas in a 1D situation (i.e. for applications to 2D migration) for information, and then switch to the materially important 2D situation (for application to 3D migration).

Principle
Numerical integration consists in calculating an approximation of the integral, over some interval [a, b], of a function f (x) from its values f at different points x i (called in the sequel interpolation nodes in the sequel) in [a, b].Let these interpolation nodes be numbered from 1 to I, assuming that x 1 = a and x I = b.A straightforward approach to numerical integration is based on interpolation procedures: we calculate a function f int (x) (int stands for interpolating) that interpolates the values f (x i ) (Lagrange interpolation) and possibly the first order derivatives f (x i ) or even higher order derivatives (Hermite interpolation) and then approxim- f int (x) dx .This procedure assumes that function f int (x) approximates function f (x) which, obviously, requires that the distance between successive points x i be small enough regarding the roughness of function f .This section introduces the basic elements involved in numerical integration in a simple 1D situation, enabling us to be concise in the more complex 2D situation.

Lagrange Interpolation
We begin with the simplest interpolation procedure, namely the Lagrange interpolation: this interpolation only uses the values f (x i ).We also deal with its simplest implementation (known as P 1 finite element), in which the interpolation function is piecewise linear (Fig. 8).
The set of these interpolation functions is the vector space generated by basis functions as described in Figure 9, each of them associated with one interpolation node.These basis functions are the piecewise first order polynomials defined by: where Thus the interpolation function f int (x) can be expressed as: Figure 8 In the simplest 1D Lagrange interpolation procedure we interpolate between points x i by a linear function.The interpolation function has C 0 regularity.
Lagrange basis function e i (x): this function is associated with interpolation node with abscissa x i .

Numerical Integration Formula
The numerical integration formula is obtained by integrating function f int (x) which, from ( 6), gives: with: From Equation (7) we observe that numerical integration appears as a weighted summation of the values f (x i ), the weights being the numbers w i defined in (8).From Figure 8 we also see that this numerical integration is an application of the very classical trapezoidal formula.

Error Estimate
To determine the accuracy of the numerical integration formula, it is important to calculate an error estimate denoted by E L : this estimate gives a bound of the error (i.e. the absolute value of the difference between b a f (x) dx and b a f int (x) dx).More precisely, it allows us to highlight the dependency of the error on the typical distance 4 h between interpolation nodes (this distance can, for instance, be the minimum distance between successive interpolation nodes) and on the smoothness of function.This smoothness is quantified by a norm involving the derivatives of function f (x).The error estimate, whose derivation is recalled in Appendix A, is of the form: where f the subscript and superscript L refer to the technique based on Lagrange interpolation.

Hermite Interpolation
To obtain a more accurate numerical integration procedure (i.e. a procedure yielding an error estimate showing for instance a h 3 dependency) more sophisticated interpolation procedures must be followed including techniques based on Hermite interpolation.In Hermite interpolation, we interpolate, at interpolation nodes x i , the values of the derivatives of function f (x) as well as the values of the function itself (Fig. 10).We shall restrict ourselves to the simplest case where only first order derivatives are taken into account.
(4) Not to be confused with offset vector h.
Hermite interpolation in 1D: the interpolation function It is straightforward to check that the interpolation function can be expressed as: where the basis functions e 0 i (x) and e 1 i (x) are piecewise third degree polynomials defined by: e 1 i x j = 0 and These basis functions are shown in Figure 11.It is straightforward to obtain the analytic expression of the basis functions over some interval x j, x j+1 , j = 1, . . ., I − 1: Equations ( 11) or (12) give the information for determining, in a unique way, the four coefficients defining a third Basis functions e 0 i (x) (left) and e 1 i (x) (right) associated with Hermite interpolation.
degree polynomial.The uniqueness in this determination shows that, for any Hermite interpolation data f (x i ) and f (x i ), i = 1, . . ., I, there exists a unique function which is a third degree polynomial over each interval x j, x j+1 , j = 1, . . ., I − 1, that matches the Hermite interpolation data.It is important to note that the interpolation function f int now has C 1 regularity as we have imposed a given value for its derivative at the boundaries of intervals x j, x j+1 , j = 1, . . ., I − 1.

Numerical Integration Formula
As already stated, the numerical integration formula is obtained by integrating f int (x).We obtain: dx Now the numerical integration formula appears not only as a linear combination of the values of f (x i ) but also of the values of f (x i ).Note that coefficients w 0 i and w 1 i are straightforward to calculate since functions e 0 i and e 1 i are third degree polynomials.

Error Estimate
As shown in Appendix B, the error estimate is of the form where C H f is a constant that now depends on the norm of f (x) in C 4 ([a, b]).H used as subscript or superscript refers to Hermite interpolation.Thus, for a smooth function f (having not too large 4 th order derivative), the numerical integration procedure based on Hermite interpolation is more accurate than that based on Lagrange interpolation (P 1 finite elements): the former yields an error estimate in O(h 4 ) (see Eq. ( 14)) whereas the latter yields an error estimate in O(h 2 ) (see Eq. ( 9)).

Numerical Differentiation
In the problem at hand, the values f (x i ) are not directly available.It is hence normal to try to estimate them from the values f (x i ) by a numerical differentiation procedure such as the finite difference technique.If we wish to preserve the good accuracy associated with Hermite interpolation, the numerical differentiation procedure must also be accurate enough.This section describes a numerical differentiation procedure that gives an error in O(h n ) with n chosen such that the overall error involved in the resulting numerical integration formula is of a higher order than that involved in Lagrange interpolation (Eq.( 9)).
We begin with the classical description of the numerical differentiation procedure (see for instance [6]) and the derivation of the associated error estimates.We wish to estimate the derivative of function f (x) at a given interpolation node, which we place, by an appropriate translation of the origin, at x = 0. We want the derivative to have the form of a linear combination of the values f (x i ) (the derivation operator is linear) plus an error estimate in O(h n ): Here, the subscript k refers to a local (i.e.relative to the node that has been placed at x = 0 which corresponds to k = 0) numbering (Fig. 12) of nodes (as opposed to subscript i associated with the overall node numbering).The summation on k only takes place at interpolation nodes located in the neighborhood of the one at x = 0, the latter being included.The numbers α k (as well as the number of neighbor interpolation nodes that must be included in the summation) are determined so as to yield an error estimate in O(h n ).
The derivation of the error estimate is obtained classically by means of Taylor's expansions.We write for any interpolation node x k , including x k = 0: where ξ k is between 0 and x k .
given by Equation ( 16) we find that if numbers Numerical differentiation in 1D.In this procedure we try to estimate the value of the derivative at the node at x = 0 from the values of the function at neighboring nodes.Note the use of a local numbering (index k) of nodes, the one at x = 0 corresponding to k = 0. satisfy the system: we then have: As shown in Appendix C, the term ) so that finally Equation (15) holds with approximation order n = 3.This accuracy requires numbers α k to satisfy the system of four linear equations ( 17)-( 20): hence four interpolation nodes are taken into account (including the one at x = 0).If we want Equation (15) to hold only with approximation order n = 2 we simply need system (17)-( 19) to be fulfilled, thus having to account for only three interpolation nodes.This yields the simplest implementation for a numerical estimation of derivatives that gives an overall error (see Appendix D) for the resulting numerical integration formula in O(h 3 ).

NUMERICAL INTEGRATION IN 2D
The principle of numerical integration in 2D is the same as in 1D: we construct an interpolation function and then calculate the integral, on the considered area, of the interpolation function, which is theoretically simple since the interpolation function is a piecewise polynomial.

The Interpolation Problem
In 2D interpolation consists in constructing a function f int (x, y) that assumes specified values at the interpolation nodes and, in the case of Hermite interpolation, whose gradient (or even higher order derivatives) also assumes specified values at these interpolation nodes.The interpolation nodes are, of course, distributed on the area over which the integral is to be calculated.Returning to the physics of the problem, the interpolation nodes are the midpoints involved in the common offset gather: the locations of these midpoints are defined by vectors q i , i = 1, . . ., I. Integration takes place in the midpoint coordinates (see Eq. ( 1)).We shall use the Figure 13 An example of triangulation.Note that triangles can have summits or edges in common but not pieces of edges.classical technique of finite elements (see [7] for instance) to solve this interpolation problem.More specifically, we shall use simplectic finite elements here: we start by splitting the integration domain into non overlapping triangles whose summits are the interpolation nodes (Fig. 13).This splitting is called triangulation.More precisely, we assume that in such a triangulation, no summit lies inside the edge of another triangle: triangles can have summits or edges in common but not pieces of edges.We also assume that the integration domain, denoted by Ω, is the union of the triangles included in the triangulation.

The Technology of Finite Elements
The technology of finite elements is based on discrete functions (i.e.functions depending on a finite number of parameters).These discrete functions are polynomials on each triangle of the triangulation considered, thus allowing easy computation of basic mathematical quantities: the values of the function or of some of its derivatives or its integral on some domain, etc.The nature of the polynomial function depends in particular on the problem to be solved (Lagrange versus Hermite interpolation in our case).The discrete functions considered belong to a (finite dimensional) vector space generated by some basis functions: any discrete function is a linear combination of these basis functions.The above mentioned computation of basic mathematical quantities is then straightforward once the basis functions are exhibited.In the finite elements we consider, the support of these basis functions is a set of triangles at the vicinity of a specific interpolation node (Fig. 14).It is normal to define the basis functions piece by piece, i.e. separately in each of In the different finite elements considered in this paper, the basis functions are associated with a specific vertex of the triangulation whose location is defined by vector q i (and, for Hermite interpolation, with the quantity to be interpolated at this vertex).The basis functions have their support localized on the triangles that surround the vertex considered.

Figure 15
Piecewise construction of some basis functions associated with node q i .To construct the piece of basis function associated with triangle T on the right, we use the affine invertible transformation F −1 that maps it onto the reference triangle on the left.We introduce a local (i.e.attached to the triangle considered) node numbering: node q i in the global numbering becomes q 1 in the local numbering.Note that these nodes q 1 have a specific role: the vertex shared by all the triangles in Figure 14.This affine transformation maps q 1 , q 2 , q 3 into points with coordinates (0,0), (0,1) and (1,0), respectively.The basis function on the right hand triangle is straightforwardly obtained once the basis function on the reference triangle is exhibited.This task is simplified in practice by using barycentric coordinates.
these triangles, and then to assemble the contribution of the different pieces.
It is classic to construct a piece of basis function associated with a specific triangle from a piece of basis function defined on a reference triangle, and to obtain the piece of basis function on the triangle of interest using the affine transformation that maps it into the reference triangle.This construction relies on the following general result.
Theorem: Consider two triangles (Fig. 15), one of them being considered as a reference triangle.Consider the affine invertible transformation F that maps the reference triangle T onto the other triangle T. Consider a general interpolation problem on triangle T: we seek a function within a specific space (or subspace) of polynomials (for instance, the space of third order polynomials that match a given linear equation), this space being called the space of interpolation functions and denoted by P, whose values, or the values of some derivatives, at some points (called interpolation nodes) in triangle T, match some specified values (called interpolation values).Assume that this interpolation problem in triangle T gives rise to a uniquely defined interpolation operator Π: for whatever specified interpolation values, there exists a unique interpolation function in P that matches these interpolation values.Then the associated interpolation problem on triangle T (the problem in which the interpolation values are to be matched at the interpolation nodes defined as the images by F of the interpolation nodes in triangle T) has a unique solution in the space of interpolation functions defined by the composition of polynomials in P with This result is obvious for simple interpolation problems such as the simplest implementation of Lagrange interpolation but, contrary to what may appear at first glance, is not that obvious for more complex problems, such as those associated with Hermite interpolation.The reader can refer to textbooks on the finite elements method such as [8], for a proof of this theorem.
From a practical standpoint, this theorem also gives the construction of a piece of basis function in any triangle once the basis function has been exhibited on the reference triangle.Barycentric coordinates play an important role here since they remain invariant in the affine transformation considered.In general, the use of barycentric coordinates greatly simplifies the expression of basis functions in triangle T. If the locations of the vertices of triangle T are defined by vectors q 1 , q 2 , q 3 , the barycentric coordinates λ 1 , λ 2 , λ 3 of a point whose location is defined by vector q are defined as the solution of: with: It is easy to check that the barycentric coordinates of the vertex at q 1 are (1, 0, 0) and that λ 1 remains constant while moving in a direction parallel to vector q 3 − q 2 .
The nature of the affine transformation F also plays an important role in the derivation of error estimates, which is an important point for our goal of accurate numerical integration (see [8] for a detailed analysis).All this is rather formal, merely the guideline for solving our 2D Lagrange and Hermite interpolation problems.In addressing these problems, we will stop being formal and become very concrete.

First Order Lagrange Interpolation in 2D
(P 1 Finite Elements) In 2D Lagrange interpolation, we calculate an interpolation function that assumes specified values at the interpolation nodes.The simplest choice (first order Lagrange interpolation) consists in finding an interpolation function that is affine on each triangle in the triangulation (P 1 finite elements): this function is defined in a unique way from the values specified at the interpolation nodes.This interpolation function obviously displays C 0 regularity.It is straightforward to realize that the interpolation function can be written as: where basis function e i ( q ) is the piecewise affine function such that: Its support is localized on the triangles around node q i (Fig. 14).To exhibit the piece of basis function associated with a triangle which has q i as summit we use a reference triangle: the affine transformation F −1 that maps the original triangle onto the reference triangle maps q i into the point with coordinates (0, 0) (Fig. 15).The basis function on the reference triangle is quite simple.In terms of barycentric coordinates its expression is λ 1 (1 is the number of the summit with coordinates (0, 0)).It is also the expression of the piece of basis function for the right hand triangle in Figure 15: barycentric coordinates are invariant in the affine transformation.

Numerical Integration
Thanks to Equation (24) integration of the interpolation function is quite simple.The integral is a weighted summation of the values f q i : with: In fact, since basis functions are simple first degree polynomials on each triangle, w i has a rather simple expression: it is one third of the sum of the areas of the triangles surrounding node i.

Error Estimate
We wish to estimate the error E L involved in the numerical integration formula (26), namely the error made in replacing function f (x, y) by function f int (x, y) in the integral.The derivation of the error estimate (see [8] for instance) is fairly technical and we merely give the result: Here h is the characteristic dimension of the triangles in the triangulation, constant C L f depends on the deformations of the triangles in the triangulation (the best triangles are equilateral triangles) and on the norm of function f in the Sobolev space W 2,1+ε (Ω) (ε is an arbitrary strictly positive number).Such a norm is large whenever the second order derivatives of function f are large.Note that, for this error estimate, we find the same order as in the 1D situation.

A Finite Element for 2D Hermite Interpolation
We now wish to find an interpolation function f int ( q ) which takes specified values at the vertices of the triangulation and whose gradient also takes specified values at these vertices.Function f int ( q ) thus has to fulfill the system of equations: where ∇ x and ∇ y are the x and y components of the gradient operator ∇, respectively.Quantities f q i , ∇ x f q i , ∇ y f q i are given data for i = 1, . . ., I. Note that this interpolation function is a well-known finite element that is described for instance in the book of [9] (Chapter 10: bending of plates).We also suggest [8] (Chapter 5, example 2.2) for a complete treatment of the problem.
As usual we shall now seek an interpolation function f int ( q ) which is a superposition of basis functions: where basis function e 0 i ( q ) is defined by: ∇ x e 0 i q j = 0 (34) ∇ y e 0 i q j = 0 (35) q 1 q 3 q 3 q 2 q 2 1 The three kinds of basis functions for Hermite interpolation in 2D for a reference triangle.In the first (left) the function has value one at the origin, zero at the other two vertices and, at all vertices, the derivatives along edges are zero.In the second (middle) the function has value zero at all vertices and, at all vertices, the derivatives along edges are zero except for edge [ q 1 , q 2 ] at the origin, where this derivative has value one.
In the third (right) the function has value zero at all vertices and, at all vertices, the derivatives along edges are zero except for edge [ q 1 , q 3 ] at the origin, where this derivative has value one.
e x i ( q ) by: e x i q j = 0 (36) ∇ x e x i q j = δ ij (37) ∇ y e x i q j = 0 (38) and e y i ( q ) by: e y i q j = 0 (39) ∇ x e y i q j = 0 (40) ∇ y e y i q j = δ ij (41) The above equations must be fulfilled for all i and j with values 1, . . ., I. On each triangle, each basis function has to fulfill 9 conditions: on each triangle, the basis functions are to be sought within the space of third order polynomials (in fact, a subspace of space P 3 of third order polynomials since the definition of these polynomials requires the specification of 10 coefficients).
To exhibit the piece of basis function e 0 i associated with a triangle which has q i as summit we can equivalently impose that at all summits of the triangle, the derivatives along the edges are zero and exhibit this function by using a reference triangle (Fig. 16, left).This function can be expressed in terms of barycentric coordinates (λ k is the barycentric coordinate associated with vertex number k in the numbering associated with the triangle) in closed form: Instead of exhibiting basis functions e x i and e y i , we shall exhibit basis functions for which we impose, at the different vertices, conditions on derivatives along edges instead of derivatives along the x or y directions (it is straightforward to obtain basis functions e x i and e y i once the new basis functions have been exhibited).Again making use of the reference triangle and of the node numbering attached to this triangle, the basis functions e 2 i as described in Figure 16 (middle) can be expressed in terms of barycentric coordinates as: whereas the basis functions e 3 i as described in Figure 16 (right) can be expressed as: It can be verified that the basis functions obtained by assembling the pieces described above have C 0 regularity (instead of C 1 regularity in the 1D situation).

Numerical Integration
Once again, numerical integration turns out to be a weighted summation of the different quantities that have been interpolated: with: w y i = Ω e y i ( q ) dx dy (48)

Error Estimate
The derivation of the error estimate E H involved in the numerical integration formula (45) is technical.In [8] (example V-2-2, page VI-23) we find the following result: Again, h is the characteristic dimension of the triangles in the triangulation; constant C H f depends on the deformations of the triangles in the triangulation (the best triangles are equilateral triangles) and on the norm of function f in the Sobolev space W 3,1+ε (Ω) (ε is an arbitrary strictly positive number).Such a norm is large whenever the third order derivatives of function f are large: thus the Hermite interpolation quadrature formula is more accurate than the one based on Lagrange interpolation only for smooth functions f .Note that we loose one order in h as compared with the 1D situation (Eq.( 14)): this is because the space of finite elements is no longer P 3 but just a subspace of P 3 .However, this error estimate will turn out to be quite compatible with the errors that arise from a numerical evaluation of derivatives (see next section).

Numerical Differentiation
If we wish to use a quadrature formula based on Hermite interpolation with no direct measurements of the derivatives, we must evaluate these derivatives by some numerical procedure.Here we resume for the 2D situation what was previously done for the 1D situation, restricting ourselves to a specific order of accuracy, namely n = 2.The aim is again to obtain an estimation of the gradient of function f ( q ) at a given interpolation node (which we place, by an appropriate translation of the origin, at q = 0) of the form: where first, subscript k refers to a local (i.e.relative to the node placed at q = 0 which corresponds to k = 0) numbering 5 of nodes (the global numbering is defined by subscript i ); secondly, the summation on k takes place only (5) This local node numbering must not be confused with the one attached to a specific triangle that we used previously.
at those interpolation nodes located at the vicinity of the one at q = 0; and thirdly, the vectors α k (as well as the number of neighbor interpolation nodes that have to be involved in the summation) must be determined so as to yield an error estimate in O(h 2 ).
We begin with the Taylor expansion at the vicinity of q = 0: where .denotes the dot product in R 2 and [D 2 f ( 0)] denotes the 2 × 2 matrix of second order derivatives of function f at q = 0.
If the α k 's in Equation ( 50) are chosen so as to satisfy: , a 2×2 symmetric matrix, then we obtain a second order accurate estimate (keep in mind the scaling based argument given in Appendix C) of ∇f ( q = 0) given by Equation ( 50).The total number of conditions to be fulfilled by the α k 's is 12: 2 are given by Equation (52), 4 by Equation (53), and 6 by Equation (54).Hence the α k 's can be computed by solving the linear system (52)-( 54) involving the values of function f at q = 0 and at 5 properly chosen (so that the linear system is not singular) neighbor interpolation nodes with locations q k (typical situations are shown in Figure 17).This technique can be applied for each interpolation node: if f denotes the vector whose components are the values f q i , i = 1, . . ., I (subscript i refers to the global numbering of interpolation nodes), we can construct, from this vector, the vector ∇ h f (here again, subscript h refers to the characteristic dimension of the triangles) whose components are estimates of the gradient of f evaluated at the different interpolation nodes q i , i = 1, . . ., I. Since the evaluation of this gradient involves a linear process (see Eq. ( 50)), the construction of this vector is nothing but the application of a linear operator (a I × I matrix) denoted ∇ h to vector f.The technique described above explains how to construct operator ∇ h .
It is straightforward to use this numerical differentiation technique together with the numerical integration formula based on Hermite interpolation: instead of Equation ( 45) we obtain: where w 0 is the vector that gathers the weights w 0 i defined in (46) (the first dot product involves a summation over Figure 17 Typical situations for the choice of 5 neighboring interpolation nodes, with locations q k , in the vicinity of the interpolation node q 0 , for the numerical differentiation procedure.Left: a simple choice of 5 neighbors q k among those associated with the 6 neighbors of q 0 .Middle: necessity to find one summit, q 5 , in a neighboring triangle.Right: as q 0 is on the boundary of the integration domain, it appears necessary to find two other summits, q 4 and q 5 , in two neighboring triangles. subscript i), and w is the vector (with the same dimensions as ∇ h f) that gathers the weights w x i and w y i defined in ( 47) and (48), respectively (the second dot product involves a summation over all components of the vectors considered, namely over subscript i and over the x and y components).
As shown in Appendix D, the error made when using numerical integration formula (55) is in O h 3 .

Application to Kirchhoff Migration (Diffraction Stack)
We shall now apply these interpolation procedures to our concrete problem, i.e. the numerical integration for Kirchhoff migration.
For given M and h, the Lagrange interpolation function mv el (M; q, h), that interpolates the different m v el (M; q i , h) defined by Equation ( 2) for i = 1, . . ., I, can be written, according to Equation (24), as: Equation ( 26) then provides the migrated image mv h at a subsurface point M as a weighted summation of the values of m v el at midpoints q i of the acquisition: Note that weights w i do not depend on the point M where the image is to be computed.Note also that, by linearity of the migration operator, w i m v el (M; q i , h) is the elementary migrated image associated with the seismic trace w i d q i , h, t v (M; q i , h) .
Thus, for a given offset gather, we can compute once and for all the weights w i associated with each midpoint q i and scale the corresponding seismic trace by multiplying it by w i .The use of the simple numerical diffraction stack (4) with scaled seismic traces yields the same result as the numerical integration formula (57): this simple preprocessing makes the mere numerical diffraction stack a consistent approximation of the continuous diffraction stack.
The same kind of procedure also applies to Hermite interpolation.The Hermite interpolation function mv el can indeed be written, accordingly to equation (32), as: where the basis functions e 0 i , e x i , e y i are defined by Equations (33) to (41).
If the derivatives are evaluated by the numerical differentiation procedure described above, the corresponding migrated image mv h at a subsurface point M is defined, according to Equation (55), by: where w 0 , w and ∇ h are defined in (55) and m v el (M) denotes the vector whose components are the values m v el (M; q i , h) of the elementary migrated images.
Introducing ∇ T h ,the transposed operator of ∇ h , the quadrature formula can be rewritten as: This formulation is interesting since vector w 0 + ∇ T h w gives the weights to be applied to each elementary migrated image m v el (M; q i , h).As in the Lagrange interpolation, these weights are independent on the subsurface point M, and a preprocessing involving a scaling of each seismic trace in the common offset gather by the corresponding weight makes the mere numerical diffraction stack a consistent and very Figure 18 Migration using a Lagrange numerical integration formula, of small offset (0.2 km) data associated with an actual marine acquisition geometry (Fig. 2).Constant Y sections (top) and constant X sections (bottom) in the depth migrated data volume obtained with the correct migration velocity (v = 3.0 km/s).accurate (error in O(h 3 )) approximation of the continuous diffraction stack.
In conclusion, our numerical integration formulas based on Lagrange and Hermite interpolation can be very easily implemented in any Kirchhoff migration program involving a superposition of elementary migrated images.The additional cost is a scaling of the seismic traces, which is quite negligible compared with the cost of running the migration itself.

Application to the Test Example
We shall now test the diffraction stack based on the Lagrange and Hermite interpolation schemes on the previously studied synthetic data built using a real marine acquisition geometry.
Using the Lagrange numerical integration formula, the corresponding migrated data associated with a small offset (0.2 km) are shown in Figure 18.A substantial improvement is observed in the quality 6 of the migrated images as compared with those obtained with the mere diffraction stack formula (4) (Fig. 2): a considerable number of migration artefacts has been removed.The same kind of result is obtained for the large offset (3 km) (Figure 20, to be compared with Figure 3).For the large offset, however the number of remaining artefacts is much larger than for the smaller offset.The remaining artefacts are caused by local sparsities in the midpoint distribution.In such a situation, as explained in the Section 2, we do not have enough samples (values taken by the elementary migrated images at the image point considered) to correctly represent the function to be integrated.This situation is Figure 19 Migration using a Hermite numerical integration formula, of small offset (0.2 km) data associated with an actual marine acquisition geometry (Fig. 2).Constant Y sections (top) and constant X sections (bottom) in the depth migrated data volume obtained with the correct migration velocity (v = 3.0 km/s).frequently encountered when migrating data associated with large offsets.
Using the Hermite interpolation scheme (which is, theoreretically, more accurate than the Lagrange scheme), the corresponding migrated data associated with a small offset (0.2 km) and with a large offset (3 km) are shown in Figures 19 and 21, respectively.The result is apparently degraded by a kind of random noise as compared with the result obtained by Lagrange interpolation.
In order to understand this surprising result we shall carry out a similar experiment with a survey involving a finer sampling in midpoints (actually the same data than those used in Figure 5).This finer sampling was obtained by adding midpoints to the original acquisition geometry.These additional midpoints are located at the barycenter of each triangle in the triangulation built from the original acquisition (Fig. 22) and the added traces correspond to offset such that h = 3 km with azimuth equal to zero.With this finer sampling interval, our numerical integration formulas can be expected to be more accurate.
The results obtained by means of Hermite interpolation are shown in Figure 23: despite irregularities in the midpoint distribution, they appear as perfect (these results should be compared with those in Figure 5 obtained by migrating the same data but with the mere diffraction stack formula (4)).Both theory and software thus seem to be correct.These good results may serve as a reference (i.e.be considered as correct) and therefore be used to evaluate the error made when computing the results with the original acquisition.The errors (Figs. 24, 25 and 26) evaluated for offset 3 km (which gave rise to the strongest artefacts) are the differences between the reference result (Fig. 23) and the results in Figures 3, 20 and 21, respectively.
Comparing the observed errors we note that: -the use of a consistent numerical integration formula drastically improves the quality of the migrated images as compared with those obtained with the mere diffraction stack; -in the vicinity of the reflector, results obtained by means of Hermite interpolation are more accurate than those obtained by Lagrange interpolation; -at some distance above the reflector, the use of Hermite interpolation reveals many more migration smiles, hence the random noise observed in Figures 19 and 21.
These experimental results are well explained by the theory: -the Hermite based integration formula can be expected to be more accurate than the formula based on Lagrange interpolation when the function to be integrated is smooth, that is, as shown in Figure 6, at the vicinity of the reflector; -however, if the function to be integrated is rough (typically the situation shown in Figure 7), that is, at some distance above the reflector, the Lagrange based integration formula turns out to be more robust.
This relative robustness, however, is still not sufficient to achieve a correct numerical integration when we move far above the reflector: the midpoint distribution is then not dense enough to make the samples correctly represent such a rapidly varying function.

Application to a Real Data Set
We shall now compare the behavior of the mere diffraction stack and the diffraction stack based on the Lagrange interpolation scheme on a real data set.The Lagrange interpolation scheme was chosen because this technique appeared in the previous section as the most robust.
The real data set came from a marine acquisition in the North Sea.This data set covers a narrow salt body (about 3 km wide).The data acquisition system was a single-source double-streamer.In theory, the spacing between both the streamers and the navigation lines was 75 m.The spacing between both the shots and the receivers was 25 m.From the prestack data, we built a cube of near offset data (offset range [80 m, 100 m]). Figure 27 shows the midpoint distribution of a part of the acquisition and the triangulation that was used for the Lagrange interpolation scheme (again the midpoints  We applied a 3D prestack common offset time migration software to the cube of near offset data.This software carries out a mere diffraction stack but takes account of the true locations of the source-receiver pairs.The migration velocity model is a 1D model with constant gradient.Although this model is quite simple, it provides a reasonable approximation of the velocity distribution in the sediments surrounding the salt body.A cross-line section in the cube of migrated data is shown in Figure 28.Many migration smiles can be observed.This is not surprising given the irregularity in the midpoint distribution.Note that the migration smiles could also have been caused by inaccuracies in our migration velocity model.However the effect of these inaccuracies should not be important in the shallower part of the migrated image.Thus we focus our interest on this shallower part namely in depths less than 1 second two-way travel-time.Figure 29 shows the result obtained by applying the same migration software to the data after they have been scaled (preprocessing associated with the Lagrange interpolation scheme) so as to compensate for irregularities in the midpoint distribution.This figure shows that the migration smiles, instead of disappearing, appear to be strengthened: in fact, large triangles have been generated in those areas where the midpoint distribution is sparse (i.e. at the vicinity of acquisition holes), hence giving huge weights to the seismic traces associated with the vertices of these triangles (Figs. 30 and 31) To overcome this difficulty, we set an upper bound, namely 1800, to the weight values.In doing so, the migration result (Fig. 32, top) appears somewhat clearer but basically Cross-line section (the same as in Figure 28) in the cube of real data migrated using the scheme based on Lagrange interpolation.The migration smiles appear strengthened compared with the result in Figure 28.quite similar to the one obtained using the mere (unweighted) diffraction stack (Fig. 28). Figure 32 bottom shows the differences between the two results (after scaling to allow quantitative comparison).The difference basically consists, as expected, of migration smiles caused by the irregular midpoint distribution.On average, the amplitude of the difference turns out to be ten times smaller than the amplitude of the individual results: we can thus understand why the results obtained using the mere diffraction stack and those based on a consistent numerical formula yield appeared so similar!In conclusion, the application of a consistent numerical integration formula to carry out a Kirchhoff migration does not drastically improve the quality of migrated images in practice.Moreover, this approach cannot compensate for large acquisition holes: yet, the weights have to be truncated to avoid spurious migration smiles arising from those areas where the midpoint distribution is sparse.However, once this upper bound has been properly chosen, the method can be applied at negligible additional cost to yield a result at least more satisfactory from a quantitative standpoint.Hence the main context for an application of this approach is quantitative migration (see e.g.[10]) for reservoir characterization.

CONCLUSION
Theoretically speaking, Kirchhoff migration cannot be performed with a mere discrete diffraction stack since, for irregularly sampled data, i.e. for a non uniform midpoint distribution in an acquisition, it yields an inconsistent approximation  of the continuous version of the diffraction stack.In fact, using a mere discrete diffraction stack gives a result which is a superposition of elementary migrated events (migration smiles) and this superposition does not give rise to the construction of the migrated event (the envelope of the elementary migrated events).To carry out a migration with non uniformly distributed midpoints, a consistent approximation of the continuous diffraction stack must be used.This can be done by implementing a numerical integration based on interpolation procedures.Two interpolation schemes have been examined: first order Lagrange interpolation and third order Hermite interpolation.The implementation of these schemes is straightforward: we can use any available code computing a mere discrete diffraction stack after a simple scaling of the input seismic traces, the scaling factor being determined by the selected (Lagrange or Hermite) numerical integration technique.In turn, since this scaling is quite cheap, using the proposed schemes does not really increase the CPU time required to run the migration.
We have shown that, when the correct migration velocity is used, the results obtained with either scheme are very good, provided the midpoint distribution is dense enough to account correctly for the variations of the function to be integrated.This also holds for an erroneous migration velocity.
The technique we have investigated is designed to compensate for irregularities in the midpoint distribution.The difficulty raised by the sparsity of midpoints is another problem.More specifically, variations in the function to be integrated can be extremely rapid and it is therefore quite common, especially when shallow dipping events are to be imaged, that the distance between midpoints is too large to allow correct imaging: the image is then contaminated by migration smiles.This difficulty is often referred to as the aliasing of the migration operator.A number of filters have been designed to remove these artefacts (for instance [4] or [5]) and implemented in industrial migration programs.The counterpart is an unavoidable degradation of the spatial resolution.Since the schemes we propose amount to a data preprocessing, they can, and must, be used in combination with these filters.
From the practical standpoint, the numerical integration formula based on Lagrange interpolation is more robust than the formula based on Hermite interpolation.An upper bound must also be set for the weights scaling the input seismic traces: lack of care at this stage can be harmful.The improvement provided by the method is in fact slight: its main field of application is for the quantitative determination of reservoir heterogeneities.
Of course, if we are interested in true amplitude imaging, a Kirchhoff migration based on a diffraction stack turns out to be simplistic, and dedicated techniques (see e.g.[10]) must be used.These techniques, obviously, also require numerical integration when applied to real, discrete data.The ones presented in this paper also apply to these more sophisticated imaging techniques, provided the function to be integrated is modified accordingly.
Finally, we have concentrated on common offset depth migration, but the techniques proposed are also applicable to time migration or when imaging is carried out in another domain: shot record domain, common angle domain ( [11] for instance), etc.We find that: We introduce the higher order error function ε 1 (x) defined by ε 1 (x) = f (x) − f 1 int (x).This error function can be rewritten as: Because of (A4), ε 1 (x) has at least 3 zeros in [x i , x i+1 ].
Hence, by Rolle's theorem, d 2 ε 1 dx 2 (x) has at least one zero ξ in ]x i , x i+1 [.Hence there exists ξ ∈]x i , x i+1 [ (which, of course, depends on u) such that: and therefore: Noting that d 2 f int dx 2 (ξ) = 0 and that w (ξ) = 2 (see (A3)), we finally obtain ∀u ∈]x i , x i+1 [: The total quadrature error E L is then: Noting that we obtain: Finally, introducing the minimum sampling interval h = min i=1,...,I−1 (x i+1 − x i ) and assuming that the discretization is such that the ratio between the largest and the smallest sampling intervals is bounded by a bound B in O( 1 where e 0 i (x) and e 1 i (x) are third degree polynomials on each interval [x i , x i+1 ] defined by: e 0 i x j = δ ij ; de 0 i dx (x j ) = 0 ∀i, j ∈ [1, I] (B14) e 1 i x j = 0; Our goal is to estimate the error involved in the interpolation procedure and, from this to derive the error involved in the associated quadrature formula.To do this, we use an adaptation of a very classical technique ( [12] for instance).
On interval ]x i , x i+1 [ we introduce an arbitrary additional node with abscissa u and construct another interpolation function f 1 int (x) defined by: with w(x) = (x − x i )(x − x i+1 ).
We find that: and that: and considering that f int is identically zero, we end up with: The total quadrature error E H is then: Noting that:

D OVERALL ERROR ESTIMATE OF HERMITE INTERPOLATION IN 1D INCLUDING NUMERICAL DERIVATION
When a numerical derivation procedure is used to carry out a Hermite interpolation, the interpolation function obtained is f int (x) whereas it would have been f int (x) if the derivatives f (x i ), (for i = 1, . . ., I) had been available.Regarding numerical integration, we compute where f num (x i ) is the numerical evaluation of f (x i ) by finite differences.We have seen that, when this evaluation makes use of two interpolation nodes in the vicinity of x i , the error

Figure 1 Synthetic
Figure 1 Synthetic 3D model.Map of interface (left) and different sections in the model: -top row: constant Y sections, left to right Y = 2 km, Y = 4 km, Y = 6 km; -bottom row: constant X sections, left to right X = 1 km, X = 7 km, X = 17 km.The velocity above the reflector is constant (v = 3.0 km/s)

Figure 3 Large
Figure 3 Large offset (3 km) migrated synthetic data associated with an actual marine acquisition geometry: -constant Y sections (top); -constant X sections (bottom).In the depth migrated data volume obtained with the exact migration velocity (v = 3.0 km/s).The amplitude of the imaging artefacts is ten times lower than that of the migrated event corresponding to the reflector of the model.

Figure 5 Large
Figure 5 Large offset (3 km) migrated synthetic data associated with the refined actual marine acquisition geometry: constant Y sections (top) and constant X sections (bottom) in the depth migrated data volume obtained with the exact migration velocity (v = 3.0 km/s).The amplitude of the imaging artefacts is ten times lower than that of the migrated event corresponding to the reflector of the model.

Figure 20 Migration
Figure 20Migration using a Lagrange numerical integration formula, of large offset (3 km) data associated with an actual marine acquisition geometry (Fig.3).Constant Y sections (top) and constant X sections (bottom) in the depth migrated data volume obtained with the correct migration velocity (v = 3.0 km/s).

Figure 21 Migration
Figure 21Migration using a Hermite numerical integration formula, of large offset (3 km) data associated with an actual marine acquisition geometry (Fig.3).Constant Y sections (top) and constant X sections (bottom) in the depth migrated data volume obtained with the correct migration velocity (v = 3.0 km/s).

Figure 22 Densification
Figure 22Densification of midpoints in an acquisition.The midpoints in the original acquisition are shown by circles: they are the vertices of the triangles in the original triangulation.New midpoints (squares) are added to the original acquisition as the barycenters of the triangles from the original triangulation.They serve to refine the original triangulation: the refined triangulation is obtained by splitting each triangle in the original triangulation into three pieces.

Figure 23 Migration
Figure 23 Migration using a Hermite numerical integration formula, of large offset (3 km) data associated with the refined actual marine acquisition geometry.Constant Y sections (top) and constant X sections (bottom) in the depth migrated data volume obtained with the correct migration velocity (v = 3.0 km/s).

Figure 24 Evaluation
Figure 24Evaluation of the error involved using the mere diffraction stack formula (4): this error was evaluated by computing the difference between the results shown in Figure3(obtained by the mere diffraction stack) and those shown in Figure23(used as reference).Constant Y sections (top) and constant X sections (bottom).

Figure 25 Evaluation
Figure 25Evaluation of the error involved in the migration scheme based on Lagrange interpolation: this error was evaluated by computing the difference between the results shown in Figure20(obtained by numerical integration) and those shown in Figure23(used as reference).Constant Y sections (top) and constant X sections (bottom).

Figure 26 Evaluation
Figure 26Evaluation of the error involved in the migration scheme based on Hermite interpolation: this error was evaluated by computing the difference between the results shown in Figure21(obtained by numerical integration) and those shown in Figure23(used as reference).Constant Y sections (top) and constant X sections (bottom).

Figure 27 Part
Figure 27 Part of the triangulation of another marine acquisition from the North Sea.

Figure 28 Cross
Figure 28Cross-line section in the cube of real data migrated using the mere diffraction stack.
Figure 30 Histogram of the weight distribution (output of the Lagrange interpolation scheme) used to scale the traces in order to compensate for irregularities in the midpoint distribution.Large weight values are associated with midpoints surrounded by large triangles.
Figure 31Color plot of weight values as a function of midpoint location (zoom).

Figure 32 Top
Figure 32 Top: same as the left part of Figure 29 but with weight values truncated above 1800.Bottom: difference between the results in Figures 32 top and 28, after normalization of these results.

12 I− 1 i=1B 3
) (B, which is of course less than (b − a)/h, must be in O(1) otherwise the nodes would not yield a genuine discretization of [a, b]), we obtain:E L ≤ C h 3 ≤ C(b − a) 12 B 3 h 2 (A12)B ERROR ESTIMATE OF HERMITE INTERPOLATION IN 1DWe use a finite element technique to implement a Hermite interpolation of a function f (x) on [a, b].We have discretized interval [a, b] by introducing some nodes with abscissas x i , for i = 1, . . ., I, these nodes being such that x 1 = a and x I = b.The nodes may be irregularly distributed on [a, b].Function f being given in C 4 ([a, b]), we denote by f int (x) the corresponding interpolation function.It is given by: = f (x) for x = x i and x = x i+1 (B19) Besides, ε 1 (x) = f (x) − f 1 int (x) satisfies: ε 1 (x) = ε(x) − w 2 (x) w 2 (u) ε(u) (B20)Now, because of (B18) and (B19), ε 1 (x) has at least one single and 2 double distinct zeros, i.e. at least 5 zeros in [x i , x i+1 ].By recursive application of Rolle's theorem, we determine that d 4 ε 1 dx 4 has at least one zero ξ (which, of course, depends on u) in ]x i , x i+1 [. ξ) = 4! we obtain:

α k x 4 k d 4 f
− x i )5  (B26) Introducing h and B as defined in Appendix A, the same reasoning leads to the result:E H ≤ D 4!30 (b − a)B 5 h 4 (B27)C SCALING ARGUMENTThe α k 's are the solution of system (17)-(20) for given x k 's.If, by scaling of discretization, each x k becomes λx k , each α k becomes α k /λ.Besides, h is the characteristic size of the intervals involved in the discretization: this quantity can be viewed as the result of a scaling by h of an original discretization.This shows that the α k 's are O(h −1 ), so that, for a function f (x) ∈ C 4 ([a, b]), the quantityk dx 4 (ξ k ) is O(h 3 ).

f
int (x) dx instead of b a f int (x) dx.The overall error ĒH is then given by:ĒH = b a f (x) − f int (x) dx ≤ b a f (x) − f int (x) dx + b a f int (x) − f int (x) dx (D28) Then ĒH ≤ E H + E D (D29)where E H is the error evaluated in Appendix B. E D is the error arising from the numerical evaluation of the derivatives, i.e.:E D = b a f int (x) − f int (x) dx = I i=1 f (x i ) − f num (x i ) ). Besides, a straightforward calculation shows that: and h as defined in Appendix A.Hence finally E D is O(h 3 ), and, according to the estimation of E H given in Appendix B, so is ĒH .