# Patent

## Title
Method of seismic processing of the type comprising depth migration before addition

## Abstract

A seismic processing method of the migration type, comprising the steps of:
linearly combining downgoing waves generated at a plurality of shotpoints, and also upgoing waves recorded by a plurality of seismic sensors;

    
    
        notionally propagating the composite waves as obtained in this way in order to obtain migrated downgoing and upgoing waves for different depths; and
        determining at a plurality of depths at least one characteristic of the subsoil as a function of the upgoing and downgoing waves propagated in this way;
    
    


the downgoing (resp. upgoing) composite waves being calculated by a linear combination of downgoing (resp. upgoing) waves in which said downgoing (resp. upgoing) waves are weighted by the coefficients of a spatial modulation matrix which is a function of the positions of the shotpoints.

## Background
CROSS-REFERENCE TO RELATED APPLICATION
This application claims priority to French Application No. 0605123 filed Jun. 9, 2006 entitled “A Method of Seismic Processing of the Type Comprising Depth Migration Before Addition”.
The present invention relates to a method of seismic processing of the type comprising depth migration before addition.
More precisely, it relates to a method of this type comprising a prior step of modulating shotpoints (the term “shot-record” can also be used instead of “shotpoints”).
Depth migration prior to addition is a central step in seismic processing. It consists in focusing seismic events recorded in time as reflections indexed by depth. One of the most precise ways of performing this step is migration by shotpoints, which consists in numerically propagating an incident wave representing the seismic source and a reflected wave. The downgoing wave dn(f,x,y,z) is initialized at the surface at depth z=0, and for each frequency f, by a synthesized representation of the source, and the upgoing wave un(f,x,y,z) is initialized by the wave recorded by the seismic sensors. Numerical propagation propagates these waves stepwise through layers of thickness Dz enabling the downgoing and upgoing waves dn(f,x,y,z) and un(f,x,y,z) to be obtained for any depth z of the grid. The result of the migration is reflectivity calculated by summing over all these frequencies and all the shots the cross-correlation of the downgoing and upgoing waves:




  
    
      
        
          r
          ⁡
          
            (
            
              x
              ,
              y
              ,
              z
            
            )
          
        
        =
        
          
            ∑
            
              f
              ,
              n
            
          
          ⁢
          
            
              
                
                  d
                  n
                
                ⁡
                
                  (
                  
                    f
                    ,
                    x
                    ,
                    y
                    ,
                    z
                  
                  )
                
              
              _
            
            ⁢
            
              
                u
                n
              
              ⁡
              
                (
                
                  f
                  ,
                  x
                  ,
                  y
                  ,
                  z
                
                )
              
            
          
        
      
    
    
      
        (
        1
        )
      
    
  




The drawback of that method is that the set of shotpoints in the seismic acquisition un(f,x,y,z) is represented by a single image r(x,y,z). In general, it is desired to have a migration result that is a gather ri(x,y,z) such that:




  
    
      
        
          r
          ⁡
          
            (
            
              x
              ,
              y
              ,
              z
            
            )
          
        
        =
        
          
            ∑
            i
          
          ⁢
          
            
              r
              i
            
            ⁡
            
              (
              
                x
                ,
                y
                ,
                z
              
              )
            
          
        
      
    
    
      
        (
        2
        )
      
    
  




This makes it possible to analyze velocities by verifying that the seismic arrivals have the same arrival times over all the gathers. It is also possible to perform processing to attenuate multiple arrivals by rejecting events that present curvature in the depth-i plane of the gather.
STATE OF THE ART
Migration by shotpoint makes it possible to produce a gather by shotpoint in which rn(x,y,z) is the reflectivity provided by shotpoint n. However, this is not convenient given the large number of shotpoints (105 to 106) in modern acquisition. In general, it is desired for gathers to comprise several tens of components.
Plane wave migration (Schultz and Claerbout (1978), Rietveld et al. (1992), Duquet et al. (2001)) is a method in which composite shots are constructed from individual shots, each composite shot being obtained by summing the individual shots after applying linear delays proportional to a given slowness value p (where slowness is the inverse of velocity), or in equivalent manner by constructing the composite shot of index m by weighting in the frequency domain the shot of index n on abscissa xn by Cmn(f)=exp(−2jpfpmxn). This linear superposition of downgoing waves and upgoing waves is followed by migration analogous to migration by shotpoint, with the composite shotpoints replacing the individual shotpoints. That technique has two drawbacks: it does not enable exactly the same image to be obtained as is obtained by shotpoint migration, and it does not give a criterion on the number of values of p to be migrated.
Romero et al. (2000) used composite shots obtained from other weightings that do not depend on the position of the source xn (linear phases, random phases, chirp, modified chirp). Those methods suffer from the defects of the preceding methods and in addition they do not make it possible to obtain interpretable gathers.
OBJECTS AND GENERAL PRESENTATION OF THE INVENTION
An object of the invention is to propose a method of processing with migration based on a composition of shots that enables exactly the same image to be obtained as is obtained by migrating shotpoints.
Another object of the invention is to propose a processing method of this type with a minimum number of composite shots to be migrated.
Yet another object of the invention is to propose a processing method with migration that enables gathers to be obtained that are analogous to those obtained when performing plane wave migration.
Specifically, the invention provides a seismic processing method of the migration type, comprising the steps of:

    
    
        linearly combining downgoing waves generated at a plurality of shotpoints, and upgoing waves recorded by a plurality of seismic sensors;
        notionally propagating the composite waves as obtained in this way in order to obtain migrated downgoing and upgoing waves for different depths; and
        determining at a plurality of depths at least one characteristic of the subsoil as a function of the upgoing and downgoing waves propagated in this way;
    
    


the method further comprising the step of calculating the downgoing (resp. upgoing) composite waves by a linear combination of downgoing (resp. upgoing) waves in which said downgoing (resp. upgoing) waves are weighted by the coefficients of a spatial modulation matrix which is a function of the positions of the shotpoints.



DESCRIPTION OF THE FIGURES
FIG. 1 represents in an f-k plane an example of composite shots to be migrated when performing migration by modulating shotpoints in accordance with an embodiment of the invention.
FIG. 2 represents in an f-k plane an example of composite shots to be migrated for plane wave migration.
FIGS. 3 and 4 are images obtained by shot modulation migration (FIG. 3) and for the gather (FIG. 4), applicable to a synthetic data set known as the Sigsbee model and to an embodiment of the invention.
FIG. 5 is a diagram showing the steps of an embodiment of the invention.



DESCRIPTION OF ONE OR MORE EMBODIMENTS OF THE INVENTION
General Principles
Condition for Exactness of Migration by Composite Shotpoint
Since migration by shotpoint is the most exact of migrations, the term “condition of exactness” can be used to designate a condition guaranteeing that the image obtained by some other migration is the same as that obtained by shotpoint migration. Such a condition is obtained below for composite shot migration.
The downgoing and upgoing waves at depth z=0 and corresponding to n individual shots of an acquisition are written dn(f,x,y,z=0) and un(f,x,y,z=0) for n=0,N−1, and migration by composite shotpoints defines a composition matrix C(f) of dimension (M,N) enabling the following composite downgoing and upgoing waves to be obtained

Dm(f,x,y,z=0) and Um(f,x,y,z=0):





  
    
      
        
          
            
              U
              m
            
            ⁡
            
              (
              
                f
                ,
                x
                ,
                y
                ,
                
                  z
                  =
                  0
                
              
              )
            
          
          =
          
            
              ∑
              
                n
                =
                0
              
              
                N
                -
                1
              
            
            ⁢
            
              
                
                  C
                  mn
                
                ⁡
                
                  (
                  f
                  )
                
              
              ⁢
              
                
                  u
                  n
                
                ⁡
                
                  (
                  
                    f
                    ,
                    x
                    ,
                    y
                    ,
                    
                      z
                      =
                      0
                    
                  
                  )
                
              
            
          
        
        ,
        
          

        
        ⁢
        
          m
          =
          0
        
        ,
        
          M
          -
          1
        
      
    
    
      
        (
        3
        )
      
    
  




with the same linear relationship relating Dm to dn. Thereafter, shotpoint migration is performed for these two composite waves. To do that, the two waves are extrapolated in depth and the linearity of the extrapolation guarantees that the linear relationship (3) is valid for every depth z between Um(f,x,y,z) and un(f,x,y,z) and similarly between Dm(f,x,y,z) and dn(f,x,y,z). This can be written in matrix notation, using the vectors:

d(f)=[d1(f,x,y,z), . . . , dN−1(f,x,y,z)]T

and

D(f)=[D1(f,x,y,z), . . . , DM−1(f,x,y,z)]T

and the matrix C(f) of elements Cmn(f):

D(f)=C(f)d(f),U(f)=C(f)u(f)  (4)

The image obtained by composite shot migration is:




  
    
      
        
          
            
              
                R
                ⁡
                
                  (
                  
                    x
                    ,
                    y
                    ,
                    z
                  
                  )
                
              
              =
              
                
                  ∑
                  f
                
                ⁢
                
                  
                    ∑
                    
                      n
                      =
                      0
                    
                    
                      M
                      -
                      1
                    
                  
                  ⁢
                  
                    
                      
                        
                          D
                          m
                        
                        ⁡
                        
                          (
                          
                            f
                            ,
                            x
                            ,
                            y
                            ,
                            z
                          
                          )
                        
                      
                      _
                    
                    ⁢
                    
                      
                        U
                        m
                      
                      ⁡
                      
                        (
                        
                          f
                          ,
                          x
                          ,
                          y
                          ,
                          z
                        
                        )
                      
                    
                  
                
              
            
          
        
        
          
            
              =
              
                
                  ∑
                  f
                
                ⁢
                
                  
                    
                      D
                      *
                    
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                  ⁢
                  
                    U
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                
              
            
          
        
        
          
            
              =
              
                
                  ∑
                  f
                
                ⁢
                
                  
                    
                      d
                      *
                    
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                  ⁢
                  
                    
                      C
                      *
                    
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                  ⁢
                  
                    C
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                  ⁢
                  
                    u
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                
              
            
          
        
      
    
    
      
        (
        5
        )
      
    
  




If the matrix C(f) is unitary for all f, then C*(f)C(f):





  
    
      
        
          
            
              
                R
                ⁡
                
                  (
                  
                    x
                    ,
                    y
                    ,
                    z
                  
                  )
                
              
              =
              
                
                  ∑
                  f
                
                ⁢
                
                  
                    
                      d
                      *
                    
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                  ⁢
                  
                    
                      C
                      *
                    
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                  ⁢
                  
                    C
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                  ⁢
                  
                    u
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                
              
            
          
        
        
          
            
              =
              
                
                  ∑
                  f
                
                ⁢
                
                  
                    
                      d
                      *
                    
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                  ⁢
                  
                    u
                    ⁡
                    
                      (
                      f
                      )
                    
                  
                
              
            
          
        
        
          
            
              =
              
                r
                ⁡
                
                  (
                  
                    x
                    ,
                    y
                    ,
                    z
                  
                  )
                
              
            
          
        
        
          
            
              =
              Id
            
          
        
      
    
    
      
        (
        6
        )
      
    
  




which means that migration by composite shot with a unitary composition matrix produces the same image as migration by shotpoint.

Migration by Shotpoint Modulation

Consideration is given to migration by composite shotpoint where the number of composite shots M is equal to the number of original shots N and where the composition matrix is, for all frequencies, the matrix of the discrete spatial Fourier transform:




  
    
      
        
          
            C
            mn
          
          ⁡
          
            (
            f
            )
          
        
        =
        
          
            C
            mn
          
          =
          
            
              1
              
                N
              
            
            ⁢
            
              exp
              ⁡
              
                (
                
                  2
                  ⁢
                  jπ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    mn
                    N
                  
                
                )
              
            
          
        
      
    
    
      
        (
        7
        )
      
    
  




This matrix (which does not depend on f) is unitary so that the migration produces exactly the same image as migration by shotpoint. The equality of the images stems directly from Parseval's theorem.
As a general rule, the positions of the shots are regularly spaced apart on a grid of pitch Dx. The position of the shot of index n is then given by xn=nDx. Since there are N shots, the total length of the set of shots is xmax=NDx. The pitch is defined in terms of wave number Dk=1/xmax which makes it possible to define a grid of wave numbers km=mDk. This notation gives:




  
    
      
        
          C
          mn
        
        =
        
          
            1
            
              N
            
          
          ⁢
          
            exp
            ⁡
            
              (
              
                2
                ⁢
                jπ
                ⁢
                
                  
                
                ⁢
                
                  k
                  m
                
                ⁢
                
                  x
                  n
                
              
              )
            
          
        
      
    
    
      
        (
        8
        )
      
    
  




which shows that the composite shot of index m is obtained by summing individual shots after modulation with the wave number km.

If starting from N individual shots, it is decided to migrate M<N modulated components, then xmax=MDx is defined which becomes a periodization distance. The following are defined in the same manner Dk=1/xmax and km=mDk, and Cmn retains the expression (8). The properties of the Fourier transform ensure, in this case, that the image obtained by migration by shotpoint modulation is the same as the image that would have been obtained by migration in which the downgoing and upgoing waves were periodized:




  
    
      
        
          R
          ⁡
          
            (
            
              x
              ,
              y
              ,
              z
            
            )
          
        
        =
        
          
            ∑
            f
          
          ⁢
          
            
              ∑
              
                n
                =
                0
              
              
                M
                -
                1
              
            
            ⁢
            
              
                
                  (
                  
                    
                      ∑
                      k
                    
                    ⁢
                    
                      
                        d
                        
                          n
                          +
                          kM
                        
                      
                      ⁡
                      
                        (
                        
                          f
                          ,
                          x
                          ,
                          y
                          ,
                          z
                        
                        )
                      
                    
                  
                  )
                
                *
              
              ⁢
              
                (
                
                  
                    ∑
                    k
                  
                  ⁢
                  
                    
                      u
                      
                        n
                        +
                        kM
                      
                    
                    ⁡
                    
                      (
                      
                        f
                        ,
                        x
                        ,
                        y
                        ,
                        z
                      
                      )
                    
                  
                
                )
              
            
          
        
      
    
    
      
        (
        9
        )
      
    
  




Expression (9) shows that the image obtained is the same as that obtained by shotpoint migration providing the downgoing and upgoing waves dn and un are decorrelated for |n−n′|≧M. This applies if M is equal to the number of individual shots N.
Number of Composite Shots to be Migrated
Nevertheless, it is the practice in shotpoint migration to define an illumination distance xill, and to consider that it is pointless to propagate downgoing and upgoing waves outside a range [−xill, xill] centered on the position of the source, since the waves are zero or negligible outside said range.
Under such circumstances, a periodization distance xmax=2xill suffices to guarantee the exactness of migration by shotpoint modulation since, because a downgoing wave dn and an upgoing wave un+M do not overlap, cross-correlation between them does not contribute to the image.
The number of composite shots to be migrated for migration by shotpoint modulation can be further reduced. A maximum slowness parameter pmax can be defined which can be calculated from a reference velocity v and a maximum angle θmax by pmax=sin θmax/v. By default, v is the surface speed and θmax=90°. A number of components to be migrated by frequency Mf is defined such that kmax=(Mf/2)Δk is the smaller of the following two values: the spatial Nyquist interval kNyz=1⁄2Δx and fpmax. This means that components in k lying outside the range [−fpmax, fpmax] are not migrated. This does not degrade the image since it is equivalent to reducing the number of individual shots by resampling them over the largest grid that satisfies the spatial Nyquist criterion for the given frequency. This means that migration by shotpoint modulation can be just as exact while being less expensive than migration by shotpoint.
Obtaining Gathers Indexed by p During Migration by Shotpoint Modulation
The above-described algorithm, which can be referred to as migration by shot modulation, makes it possible to obtain the same migration as when performing shotpoint migration but for smaller cost. There follows a description of a method suitable, during such a migration, for obtaining gathers indexed by slowness p. The range [−pmax, pmax] is subdivided into Np ranges of width Δp. During migration of the composite shot corresponding to the frequency and to the wave number (f,km), the contribution to the image of the composite shot is accumulated over the range in p corresponding to p=km/f. This is equivalent to resampling a regular grid of pitch Δk using a regular grid of pitch fΔp, which can be done in various ways such as linear interpolation, cardinal sine, or band-limited interpolation, etc. . . . . This resampling can be performed using precalculated weightings λi(f,km), i=0, Np−1 such that Σiλi(f,km)=1 which are used during the migration to update the various gathers by:




  
    
      
        
          
            R
            i
          
          ⁡
          
            (
            
              x
              ,
              y
              ,
              z
            
            )
          
        
        =
        
          
            ∑
            
              f
              ,
              m
            
          
          ⁢
          
            
              
                λ
                i
              
              ⁡
              
                (
                
                  f
                  ,
                  
                    k
                    m
                  
                
                )
              
            
            ⁢
            
              
                R
                
                  f
                  ,
                  m
                
              
              ⁡
              
                (
                
                  x
                  ,
                  y
                  ,
                  z
                
                )
              
            
          
        
      
    
    
      
        (
        10
        )
      
    
  




where Rf,m(x,y,z) is the contribution to the image of the composite shot corresponding to (f,km);

Rf,m(x,y,z)= Dm(f,x,y,z)Um(f,x,y,z)  (11)

and Ri(x,y,z) is the element of the gather in p corresponding to the index i.

Migration by shotpoint modulation thus makes it possible to obtain gathers in p that are analogous to plane wave migration. However, the method of obtaining them is different. Migration by shot modulation uses shot composition that is independent of frequency and that corresponds to modulation, whereas migrations by plane waves uses shot composition that depends on frequency and corresponds to delays. A step specific to modulating migration, distribution depending on the frequency of the current contribution to the image amongst the gathers, enables modulated migration to obtain gathers that are indexed in p.
The advantage of migration by shotpoint modulation is that the number of shots to be migrated is smaller than with plane wave migration. This is illustrated by FIGS. 1 and 2, where each cross corresponds to a composite shot. Under such circumstances, the cost of migration by shotpoint modulation is 42% the cost of migration by plane wave.
Migration by shotpoint modulation can use a composition matrix that depends on frequency, e.g. by using a periodization distance xmax that depends on frequency xmax(f).
Implementation Example
FIG. 5 shows the main steps (modulation, linear combination, then calculating gathers) in an example implementation of shotpoint migration of the type described above.
Shot modulation migration has two main parameters: the periodization distance xmax and the maximum slowness pmax. For N being the number of shots and Δx the distance between shots, the composition matrix is as follows:




  
    
      
        
          
            C
            mn
          
          =
          
            
              1
              
                N
              
            
            ⁢
            
              exp
              ⁡
              
                (
                
                  2
                  ⁢
                  jπ
                  ⁢
                  
                    
                  
                  ⁢
                  m
                  ⁢
                  
                    
                      n
                      ⁢
                      
                        
                      
                      ⁢
                      Δ
                      ⁢
                      
                        
                      
                      ⁢
                      x
                    
                    
                      x
                      max
                    
                  
                
                )
              
            
          
        
        ,
        
          

        
        ⁢
        
          m
          ∈
          
            [
            
              
                -
                
                  
                    M
                    f
                  
                  2
                
              
              ,
              
                
                  M
                  f
                
                2
              
            
            ]
          
        
        ,
        
          n
          ∈
          
            [
            
              0
              ,
              
                N
                -
                1
              
            
            ]
          
        
      
    
    
      
        (
        12
        )
      
    
  




with





  
    
      
        
          M
          f
        
        =
        
          min
          ⁡
          
            (
            
              
                2
                ⁢
                
                  fp
                  max
                
                ⁢
                
                  x
                  max
                
              
              ,
              
                
                  x
                  max
                
                
                  Δ
                  ⁢
                  
                    
                  
                  ⁢
                  x
                
              
            
            )
          
        
      
    
    
      
        (
        13
        )
      
    
  




where Mf is the number of composite shots to be migrated for the frequency f, the composite shot of index m corresponding to wave number km=m/xmax.

In order to calculate the gather of slownesses, after selecting the interval Δp between gathers:




  
    
      
        
          
            R
            i
          
          ⁡
          
            (
            
              x
              ,
              y
              ,
              z
            
            )
          
        
        =
        
          
            ∑
            
              f
              ,
              m
            
          
          ⁢
          
            
              
                λ
                i
              
              ⁡
              
                (
                
                  f
                  ,
                  
                    k
                    m
                  
                
                )
              
            
            ⁢
            
              
                R
                
                  f
                  ,
                  m
                
              
              ⁡
              
                (
                
                  x
                  ,
                  y
                  ,
                  z
                
                )
              
            
          
        
      
    
    
      
        (
        14
        )
      
    
  




where Rf,m(x,y,z) is the contribution to the composite shot image and λi(f,km) are the coefficients of a filter that resamples a grid of pitch 1/xmax to a grid of pitch fΔp.

Another Implementation Example
In a variant, instead of synthesizing modulated shots on the surface, these shots can be synthesized at depth. This procedure, described in the general case by Rietveld et al. (1992) applies to migration by shot modulation. Virtual sources are placed at a depth z0 in a medium L at abscissa points x1, and they are given respective amplitudes equal to the element of the composition matrix Cm1(f),m that corresponds to a given wave number km (the virtual sources may also be placed at different depths). For each m, the rising field generated by those virtual sources at depth is calculated at the positions of the N real sources (the positions of the individual shots, usually at depth z=0 and at abscissa points xn). Wavelets Smn(f) are obtained. The complex conjugate Cmn(f)=Smn(f)* is then applied as the composition matrix for the surface sources. The quasi-unitary properties of propagation matrices mean that the downgoing wave from the composite shot will have the desired form for a sequence of modulated shots once it has been extrapolated to the depth z0. This makes it possible to have gathers in which the complexity of propagation between the surface and the depth z0 has been removed.
The sources may also be placed at different depths.
The 3D Case
The method described can be generalized to 3D processing.
The sources are indexed n=0,N−1 and occupy positions (xn, yn) on a grid of pitch Δx in the x direction and Δy in the y direction.
Definition of parameters xmax and ymax, pmax, Δpx, pymax, and Δpy 
For each frequency, Mx(f) and My(f) are calculated as follows:




  
    
      
        
          
            
              M
              x
            
            ⁡
            
              (
              f
              )
            
          
          =
          
            min
            ⁡
            
              (
              
                
                  2
                  ⁢
                  
                    fp
                    
                      x
                      max
                    
                  
                
                ,
                
                  
                    x
                    max
                  
                  
                    Δ
                    ⁢
                    
                      
                    
                    ⁢
                    x
                  
                
              
              )
            
          
        
        ,
        
          

        
        ⁢
        
          
            
              M
              y
            
            ⁡
            
              (
              f
              )
            
          
          =
          
            min
            ⁡
            
              (
              
                
                  2
                  ⁢
                  
                    fp
                    
                      y
                      max
                    
                  
                
                ,
                
                  
                    y
                    max
                  
                  
                    Δ
                    ⁢
                    
                      
                    
                    ⁢
                    y
                  
                
              
              )
            
          
        
      
    
    
      
        (
        15
        )
      
    
  




Wave numbers kxmx and kymy are calculated indexed by mx, my with

mx=[(−Mx(f)/2,Mx(f)/2] and my=[(−My(f)/2,My(f)/2]

as follows:





  
    
      
        
          
            kx
            
              m
              x
            
          
          =
          
            
              m
              x
            
            
              x
              max
            
          
        
        ,
        
          
            ky
            
              m
              y
            
          
          =
          
            
              m
              y
            
            
              y
              max
            
          
        
      
    
    
      
        (
        16
        )
      
    
  




The downgoing and upgoing waves of the composite shots derived from the waves of the individuals shots are calculated as follows:





  
    
      
        
          
            
              D
              
                
                  m
                  x
                
                ,
                
                  m
                  y
                
              
            
            ⁡
            
              (
              
                f
                ,
                x
                ,
                y
                ,
                
                  z
                  =
                  0
                
              
              )
            
          
          =
          
            
              ∑
              
                n
                =
                0
              
              
                N
                -
                1
              
            
            ⁢
            
              
                C
                
                  
                    m
                    x
                  
                  ,
                  
                    m
                    y
                  
                  ,
                  n
                
              
              ⁢
              
                
                  d
                  n
                
                ⁡
                
                  (
                  
                    f
                    ,
                    x
                    ,
                    y
                    ,
                    
                      z
                      =
                      0
                    
                  
                  )
                
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              U
              
                
                  m
                  x
                
                .
                
                  m
                  y
                
              
            
            ⁡
            
              (
              
                f
                ,
                x
                ,
                y
                ,
                
                  z
                  =
                  0
                
              
              )
            
          
          =
          
            
              ∑
              
                n
                =
                0
              
              
                N
                -
                1
              
            
            ⁢
            
              
                C
                
                  
                    m
                    x
                  
                  ,
                  
                    m
                    y
                  
                  ,
                  n
                
              
              ⁢
              
                
                  u
                  n
                
                ⁡
                
                  (
                  
                    f
                    ,
                    x
                    ,
                    y
                    ,
                    
                      z
                      =
                      0
                    
                  
                  )
                
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              m
              x
            
            =
            
              -
              
                
                  
                    M
                    x
                  
                  ⁡
                  
                    (
                    f
                    )
                  
                
                2
              
            
          
          ,
          
            
              
                M
                x
              
              ⁡
              
                (
                f
                )
              
            
            2
          
          ,
          
            
              m
              y
            
            =
            
              -
              
                
                  
                    M
                    y
                  
                  ⁡
                  
                    (
                    f
                    )
                  
                
                2
              
            
          
          ,
          
            
              
                M
                y
              
              ⁡
              
                (
                f
                )
              
            
            2
          
        
      
    
    
      
        (
        17
        )
      
    
  




with:





  
    
      
        
          C
          
            
              m
              x
            
            ,
            
              m
              y
            
            ,
            n
          
        
        =
        
          
            1
            
              N
            
          
          ⁢
          
            exp
            ⁡
            
              [
              
                2
                ⁢
                
                  jπ
                  ⁡
                  
                    (
                    
                      
                        
                          kx
                          
                            m
                            x
                          
                        
                        ⁢
                        
                          x
                          n
                        
                      
                      +
                      
                        
                          ky
                          
                            m
                            y
                          
                        
                        ⁢
                        
                          y
                          n
                        
                      
                    
                    )
                  
                
              
              ]
            
          
        
      
    
    
      
        (
        18
        )
      
    
  




For all mx, my, Δmax,my and Umx,my are extrapolated for all values of z over the grid of pitch Δz.
The contribution of the composite shot mx, my to the image is calculated as follows:

Rf,mx,my(x,y,z)= Dmx,my(f,x,y,z)Umx,my(f,x,y,z)  (19)

Slowness vector gathers p=(px,py) indexed in ix, iy corresponding to px=ixDpx, py=iyDpy are updated as follows:




  
    
      
        
          
            R
            
              
                i
                x
              
              ,
              
                i
                y
              
            
          
          ⁡
          
            (
            
              x
              ,
              y
              ,
              z
            
            )
          
        
        =
        
          
            ∑
            
              f
              ,
              
                m
                x
              
              ,
              
                m
                y
              
            
          
          ⁢
          
            
              
                λ
                
                  
                    i
                    x
                  
                  ,
                  
                    i
                    y
                  
                
              
              ⁡
              
                (
                
                  f
                  ,
                  
                    kx
                    
                      m
                      x
                    
                  
                  ,
                  
                    ky
                    
                      m
                      y
                    
                  
                
                )
              
            
            ⁢
            
              
                R
                
                  f
                  ,
                  
                    m
                    x
                  
                  ,
                  
                    m
                    y
                  
                
              
              ⁡
              
                (
                
                  x
                  ,
                  y
                  ,
                  z
                
                )
              
            
          
        
      
    
    
      
        (
        20
        )
      
    
  




where λix,iy(f,kxmx,kymy) is a filter that interpolates a two-dimensional grid of pitch Δkx, Δky to a grid of pitch fΔpx, fΔpy.

REFERENCES


    P. S. Schultz, J. F. Claerbout (1978): Velocity estimation and downward continuation by wavefront synthesis, Geophysics, 43, 691-714.
    W. E. A. Rietveld, A. J. Berkhout, C. P. A. Wapenaar (1992): Optimum seismic illumination of hydrocarbon reservoirs, Geophysics, 57, 1334-1345.
    C. C. Ober., L. A. Romero, D. C. Ghiglia (2000): Method of migrating seismic records, U.S. Pat. No. 6,021,094, granted Feb. 1, 2000.
    P. Lailly, B. Duquet, A. Ehinger (2000): Méthode pour réaliser en 3D avant sommation, une migration de données sismiques [A method of performing 3D migration of seismic data prior to summing], French patent No. 2 784 195, published on Apr. 7, 2000.
    B. Duquet, P. Lailly, A. Ehinger (2001): 3D plane wave migration of streamer data. SEG 71st Annual International Meeting, Expanded Abstracts 1033-1036.
    L. A. Romero, D. C. Ghiglia, C. C. Ober, S. A. Morton (2000): Phase encoding of shot records in prestack migration, Geophysics, 65, 426-436.

## Claims

1. A seismic processing method of the migration type, comprising the steps of:
linearly combining downgoing seismic waves generated at a plurality of shotpoints by different seismic sources, and upgoing seismic waves recorded by a plurality of seismic sensors;
notionally propagating the composite seismic waves as obtained in this way in order to obtain migrated downgoing and upgoing seismic waves for different depths; and
determining at a plurality of depths at least one characteristic of the subsoil as a function of the upgoing and downgoing seismic waves propagated in this way; and
the method further comprising the step of calculating the downgoing (respectively upgoing) composite waves by a linear combination of downgoing (respectively upgoing) waves in which said downgoing (respectively upgoing) waves are weighted by the coefficients of a spatial modulation matrix which is a function of the positions of the shotpoints.

2. A seismic processing method according to claim 1, in which, for a line of shotpoints, a weighting applied to a downgoing wave or an upgoing wave is:




  exp
  ⁡
  
    (
    
      2
      ⁢
      j
      ⁢
      
        
      
      ⁢
      
        Πk
        m
      
      ⁢
      
        x
        n
      
    
    )
  



where km is a wave number defined by mΔkx, with m being an integer corresponding to the index of the calculated composite wave and





  
    
      Δ
      ⁢
      k
    
    x
  
  =
  
    1
    /
    
      x
      max
    
  



 where xmax is a periodization distance.

3. A method according to claim 1, in which, for a grid of shotpoints in two dimensions x and y, a weighting applied to a downgoing or upgoing wave for shots at positions xn, yn is:




  exp
  ⁡
  
    [
    
      2
      ⁢
      j
      ⁢
      
        
      
      ⁢
      
        Π
        ⁡
        
          (
          
            
              k
              ⁢
              
                
              
              ⁢
              
                x
                mx
              
              ⁢
              
                x
                n
              
            
            +
            
              k
              ⁢
              
                
              
              ⁢
              
                y
                my
              
              ⁢
              
                y
                n
              
            
          
          )
        
      
    
    ]
  



where kmx and kmy are two wave numbers given by:





  
    k
    mx
  
  =
  
    
      m
      ⁢
      
        
      
      ⁢
      x
      ⁢
      
        
      
      ⁢
      
        
          Δ
          ⁢
          k
        
        x
      
      ⁢
      
        
      
      ⁢
      and
      ⁢
      
        
      
      ⁢
      
        k
        my
      
    
    =
    
      m
      ⁢
      
        
      
      ⁢
      y
      ⁢
      
        
      
      ⁢
      
        
          Δ
          ⁢
          k
        
        y
      
    
  



with mx and my being two integers corresponding to the indices of the calculated composite wave, and





  
    
      Δ
      ⁢
      k
    
    x
  
  =
  
    1
    /
    
      x
      max
    
  



where xmax is a periodization distance, and





  
    
      Δ
      ⁢
      k
    
    y
  
  =
  
    1
    /
    
      y
      max
    
  



 where ymax is a periodization distance.

4. A method according to claim 2, or claim 3, in which the periodization distance(s) depend(s) on frequency.

5. A method according to claim 2, in which for each frequency f, a number Mf of downgoing or upgoing composite waves that are calculated and over which the migration processing is applied is Mf such that kmax=(Mf/2)Δk is the smaller of two values constituted firstly by the spatial Nyquist interval kNyq=1⁄2Δx and secondly fpmax, where pmax is a maximum slowness parameter.

6. A method according to claim 3, in which for each frequency f, the indices mx and my lie respectively in the ranges [−Mx(f)/2, Mx(f)/2] and [−My(f)/2, My(f)/2], where Mxf and Myf are such that kxmax=(Mxf/2)Δk is the smaller of two values constituted firstly by the spatial Nyquist intervals kNyq=1⁄2Δx and secondly by fpmax, where pmax is a maximum slowness parameter, and kymax=(Myf/2)Δk is the smaller of two values constituted firstly by the spatial Nyquist interval kNyq=1⁄2Δy and secondly fpmax, where pmax is a maximum slowness parameter.

7. A method according to any preceding claim, in which the migrated gathers are calculated by distributing the contribution to the image of the composite shotpoints over ranges that are a function of slowness parameters.