Abstract:
This invention relates to a method for processing seismic data comprising a collection of seismic traces with different offsets, comprising the steps of: a) breaking up one or several first trace(s) of the collection of seismic traces into a series of trace segment&#39;s according to a predetermined segmentation interval; b) defining a series of expansion coefficients, each expansion coefficient being associated with a segment of the first trace or traces; c) applying the associated expansion coefficient to each segment of the first trace or traces; d) comparing the first trace or traces thus expanded with a second trace from the collection of seismic traces to evaluate their similarity; e) repeating steps b), and d) with a new series of expansion coefficients, f) determining an optimum series of expansion coefficients that maximizes the similarity between the first expanded trace and the second trace in order to obtain one or several corrected first trace(s).

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     This national phase application is based on PCT/EP2005/056789 filed on Dec. 14, 2005 which claims priority to French Application No. 0413260 filed Dec. 14, 2004 entitled “Correction Method for Processing Seismic Traces.” 
     STATEMENT REGARDING FEDERALLY-SPONSORED RESEARCH OR DEVELOPMENT 
     N/A 
     NAMES OF THE PARTIES TO A JOINT RESEARCH AGREEMENT 
     N/A 
     INCORPORATION-BY-REFERENCE OF MATERIAL SUBMITTED ON A COMPACT DISK 
     N/A 
     BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The invention relates to the field of seismic data processing. More precisely, it relates to the inversion of seismic data. 
     2. Description of the Related Art 
     In general, for seismic exploration, a plurality of sources and receivers are distributed at ground level and at a distance from each other. The seismic sources are activated to generate seismic waves that are propagated in the subsurface. These seismic waves undergo deviations during their propagation. They are refracted, reflected and diffracted at the interfaces of the subsurface. Certain waves propagated in the subsurface are detected by seismic receivers and recorded in real time in the form of signals (called traces). The recorded traces can be processed to obtain an image of the subterranean geological structures. 
     During processing, the summing step (or stacking) consists of adding together the traces corresponding to seismic waves that are reflected at a same point in the subsurface. This step makes it possible to increase the signal-to-noise ratio and the primary-reflections-to-multiple-reflections ratio in the seismic data processed. 
     Starting from the assumption of a subsurface that is horizontally stratified without lateral variation of propagation velocities, it can be shown that the traces having the property of illuminating the same point of the subsurface for variable source-receiver distances (or offset) are those having the same mid-point in common between source and receiver. 
     However, the waves reflected in the subsurface are recorded at variable times according to the offset. Before adding the traces, it is therefore necessary to correct these traces to reduce them to a common reference, the zero offset trace. This correction is carried out during a step called NMO (Normal Move Out) correction. 
     The NMO step requires prior knowledge of a model for propagation velocities of seismic waves in the subsurface. 
     For example, the Dix model is based on the assumption that the subsurface is formed of horizontal layers or strata in which each layer is isotropic and has an associated given propagation velocity (interval velocity). The NMO correction relies on the model thus defined to correct the arrival time of a reflection recorded with a given offset x by bringing it to the theoretical time τ 0  at which it would have been recorded with a zero offset x=0. 
     Given that the velocities are not known a priori, the correction step is carried out on the seismic traces by sweeping a range of velocities. Next, only the velocity optimising the semblance of the traces as a whole is retained. 
     Thus it is possible, in the most favourable cases, to deduce an estimate of the interval velocities between the highest energy events reflected. 
     In general the estimation of interval velocities does not take into account the anisotropy of the subsurface, that is to say the variation of velocity in the layers as a function of the propagation direction. 
     Moreover, neither does the estimate of interval velocities take into account the variation in reflectivity of the subsurface interfaces as a function of the angle of incidence of the wave. 
     BRIEF SUMMARY OF INVENTION 
     An aim of the invention is to propose a method for processing seismic data making it possible to obtain more precise information about the properties of the subsurface than the methods of prior art. 
     This problem is solved in accordance with the present invention by a method for processing seismic data comprising a collection of seismic traces with different offsets, comprising the steps consisting of:
         a) breaking up one or several first trace(s) of the collection of seismic traces into a series of trace segments according to a predetermined segmentation interval;   b) defining a series of expansion coefficients, each expansion coefficient being associated with a segment of the first trace or traces;   c) applying the associated expansion coefficient to each segment of the first trace or traces;   d) comparing the first trace or traces thus expanded with a second trace from the collection of seismic traces to evaluate their similarity;   e) repeating steps b), c) and d) with a new series of expansion coefficients, and   f) determining an optimum series of expansion coefficients that maximizes the similarity between the first expanded trace or traces and the second trace in order to obtain one or several corrected first trace(s).       

     According to the method of the invention, the NMO correction is carried out for each trace, segment by segment, with a predetermined arbitrary segmentation interval, and not on the basis of the highest energy events. 
     The method according to the invention allows a global optimisation of the expansions applied to the traces without privileging certain events. 
     The method according to the invention makes it possible to obtain information about the properties of the subsurface with a higher resolution than methods of prior art. The resolution of the obtained information is directly linked with the chosen segmentation interval. In particular, it makes it possible to deduce the following parameters: 
     local P wave velocity contrasts, 
     density contrasts of the subsurface, 
     S wave velocity contrasts, 
     anisotropy parameters of the subsurface. 
     In an embodiment of the method according to the invention, the predetermined segmentation interval is a multiple of the sampling interval for recording the first seismic trace. 
     Furthermore, the second trace can be a zero offset trace of the collection of seismic traces. 
     The second seismic trace can also be a trace of an offset immediately below the first seismic trace in the seismic traces collection. 
     Furthermore, the second seismic trace can itself be a corrected trace. 
     In an embodiment of the invention, the expansion coefficients series in step b) is defined according to an algorithm with random or pseudo-random selection. 
     In particular, the expansion coefficients series maximizing the similarity between the first expanded trace and the second trace can be determined by a Monte Carlo method. 
     In an embodiment of the invention, the comparison step d) comprises the correlation of the first expanded trace with the second trace. 
     In an embodiment of the invention, step d) comprises the determination of a cost function estimating the similarity between the first expanded trace and the second trace. 
     In an embodiment of the invention, steps a) to f) are applied to each trace or groups of traces from the seismic traces collection, to obtain a corrected collection of traces. 
     In particular, steps a) to f) can be applied by increasing order of trace offset. 
     Furthermore, the method can comprise a step g) consisting of deducing, from the optimum expansion coefficient series associated with each trace, a propagation velocity of a compressional (P) seismic wave in the subsurface as a function of depth. 
     In particular, step g) can advantageously include a ray tracing substep. 
     In an embodiment of the invention, the method includes a step h) consisting of deducing the density contrast data of the subsurface as a function of depth. 
     The method can also include a step i) consisting of deducing the propagation velocity of a shear (S) seismic wave as a function of depth. 
     The method can also include a step j) consisting of deducing the anisotropy parameters of a subsurface from the series of optimum expansion coefficients associated with a plurality of traces. 
     The invention also refers to an inversion method for seismic data comprising a collection of seismic traces with different offsets, comprising the steps of:
         a) breaking up one or several first trace(s) of the collection of seismic traces into a series of trace segments according to a predetermined segmentation interval;   b) defining a series of expansion coefficients, each expansion coefficient being associated with a segment of the first trace or traces;   c) applying the associated expansion coefficient to each segment of the first trace or traces;   d) comparing the first trace or traces thus expanded with a second trace from the collection of seismic traces to evaluate their similarity;   e) repeating steps b), c) and d) with a new series of expansion coefficients;   f) determining an optimum series of expansion coefficients that maximizes the similarity between the first expanded trace and the second trace in order to obtain one or several corrected first traces.   g) deducing the values of seismic wave propagation velocities as a function of depth in the subsurface.       

     Other characteristics and advantages will become even clearer from reading the following description, given as a purely illustrative and non-limiting example, to be read with reference to the attached figures, amongst which; 
    
    
     
       BRIEF DESCRIPTION OF THE FIGURES 
         FIG. 1  shows diagrammatically the paths of seismic waves propagating between source-receiver pairs having the same common mid-point, 
         FIG. 2  shows diagrammatically the path of seismic waves propagating in the subsurface between a source and a receiver, 
         FIG. 3  shows diagrammatically a common mid-point collection of seismic traces, 
         FIG. 4  shows diagrammatically the breaking up of traces from the collection of traces into a series of trace segments according to a predetermined segmentation interval, 
         FIG. 5  shows diagrammatically the traces obtained after correction through application of expansion coefficients, 
         FIG. 6  shows diagrammatically the expansion of a trace segment, 
         FIG. 7  is a flow chart showing diagrammatically the different steps of a method for processing seismic data according to a first embodiment of the invention, 
         FIG. 8  is a flow chart showing diagrammatically the different steps of a method for processing seismic data according to a second embodiment of the invention, 
         FIG. 9  is a flow chart showing diagrammatically the different steps of a method for processing seismic data according to a third embodiment of the invention, 
         FIG. 10  shows diagrammatically the path of a seismic wave in the subsurface, 
         FIG. 11  is a diagram showing the amplitudes of trace segments as a function of the angle of incidence in a given section of the subsurface. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     On  FIG. 1 , a set of sources S 1 , S 2 , S 3  is shown together with a set of receivers R 1 , R 2 , R 3  set at ground level. The source-receiver pairs (S 1 , R 1 ), (S 2 , R 2 ), (S 3 , R 3 ) have a common mid-point O. It is supposed that the subsurface is formed of a series of horizontal layers delimited by interfaces. In this figure, the receivers R 1 , R 2 , R 3  record a same event corresponding to a reflection at the level of a minor point P on a subsurface interface, with different slants (offsets). 
       FIG. 2  shows diagrammatically the path of seismic waves propagating in the subsurface between a source S 1  and a receiver R 1  having mid-point O. At the time of a seismic acquisition, the source S 1  is activated to generate seismic waves that propagate in the subsurface. The seismic waves are reflected at the different interfaces and are recorded by the receiver R 1 . The recording obtained by the receiver R 1  contains signals corresponding to different events. 
     As shown in  FIG. 2 , during its propagation, a seismic wave undergoes deviations due to the variations in reflectivity of the subsurface interfaces as a function of the angle of incidence of the wave. 
       FIG. 3  shows diagrammatically the seismic traces A 0 , A 1 , A 2 , A 3  of amplitudes a 0 , a 1 , a 2 , a 3  recorded over time by receivers located in abscissae  0 , x 1 , x 2  and x 3 . The traces are grouped into a common mid point collection. Each trace is associated with a given offset and consequently with a given slant. 
     The traces A 0 , A 1 , A 2 , A 3  contain signals corresponding to identical events. Nonetheless, these signals are recorded at variable times t as a function of the offset x. 
       FIG. 7  is a diagram illustrating the different steps of a processing method according to a first embodiment of the invention 
     A collection of traces is considered comprising traces A 0 , A 1 , A 2 , . . . A N  ranged by increasing offset 0, x 0 , x 1 , x 2 , . . . x N  (or increasing slant). 
     According to a first step  10  (shown in  FIG. 4 ), the trace A 1  of offset x 1  is broken up into n trace segments A 1   1 , A 1   2 , A 1   3 , . . . A 1   n  of amplitudes of a 1   1 , a 1   2 , a 1   3 , . . . a 1   n  according to a predetermined segmentation interval e. The segmentation interval e can, for example, be chosen as a multiple of the sampling interval for recording seismic traces. Typically, in seismic acquisition, the sampling interval for recording traces is of the order of 2 to 4 milliseconds. The segmentation interval is chosen as a function of the subsurface image resolution one wishes to obtain and taking into consideration the calculation time this implies. 
     Breaking up the seismic traces by segments is equivalent to dividing the subsurface into n parallel horizontal sections or layers, each section having a thickness z corresponding to a propagation time e of the zero offset seismic wave (x=0). 
     According to a second step  20 , a series of n expansion coefficients d 1   1 , d 1   2 , d 1   3 , . . . d 1   n  is defined, each expansion coefficient 
             d   ⁢     k   1           
being associated with a segment
 
             A   ⁢     k   1           
of the first trace A 1 .
 
     The expansion coefficients d 1   1 , d 1   2 , d 1   3 , . . . d 1   n  are, for example, determined by a random or pseudo-random selection algorithm. The expansion coefficients are taken within predetermined ranges corresponding to the orders of expansion generally found. 
     According to a third step  30 , the associated expansion coefficient 
             d   ⁢     k   1           
determined in the preceding step is applied to each segment
 
             A   ⁢     k   1           
of the first trace A 1 . Thus one obtains a first expanded A 1 ′ with amplitude a 1 ′.
 
     According to a fourth step  40 , the first expanded trace A 1 ′ is compared with the reference trace A 0  with zero offset (or zero slant) to evaluate their similarity. 
     For this, the correlation product of traces A 1  and A 0  can be determined. 
     It is also possible to calculate a cost function F of the type: 
     
       
         
           
             
               
                 
                   F 
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         0 
                       
                       N 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         ( 
                         
                           
                             a 
                             1 
                             k 
                           
                           - 
                           
                             a 
                             
                               1 
                               ′ 
                             
                             k 
                           
                         
                         ) 
                       
                       2 
                     
                   
                 
               
               
                 
                   [ 
                   1 
                   ] 
                 
               
             
           
         
       
     
     The correlation product or the cost function is a measure of the similarity between traces A 1 ′ and A 0 . 
     Next, the second, third and fourth steps  20 ,  30 , and  40  are repeated with a new series of expansion coefficients. The expansion coefficients are again determined by the random or pseudo-random selection algorithm. 
     According to a fifth step  50 , on the basis of comparisons carried out with the different series of expansion coefficient, a series of expansion coefficients d 1   1 , d 1   2 , d 1   3 , . . . d 1   n  is determined which maximizes the resemblance between the first expanded trace A 1 ′ and the reference trace A 0 . In order to do this, a Monte Carlo non-linear optimisation method is used such as, for example, a simulated annealing method. 
     The offset correction method just described is applied to each of the traces A 1 , A 2 , . . . A N  of the seismic trace collection according to increasing order of slant (or offset). This method leads to a corrected collection of traces A 1 ′, A 2 ′, . . . A N ′ being obtained. According to a possible embodiment mode of the invention, each trace A j+1  is corrected taking a reference trace the corrected trace A j ′ of an offset immediately below trace A j+1  in the collection of seismic traces. 
     According to this embodiment, the traces of the trace collection are corrected one after the other, which leads to determination of an associated series of expansion coefficients d j   1 , d j   2 , . . . d j   n  for each trace A j . 
       FIG. 8  illustrates a second embodiment mode of the processing method whereby two traces are corrected simultaneously. In this variant, the expansion coefficients are determined for two traces A 1  and A 2  or for a partial sum of two traces. 
     This second embodiment mode makes it possible to take into account the amplitude of the traces and consequently the trace amplitude variation as a function of the offset or slant (AVO). 
     The fourth step  40  is then modified to calculate a cost function F of the type: 
     
       
         
           
             
               
                 
                   F 
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         0 
                       
                       n 
                     
                     ⁢ 
                     
                       
                         ( 
                         
                           
                             ( 
                             
                               
                                 
                                   sin 
                                   2 
                                 
                                 ( 
                                 
                                   ar 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   
                                     cos 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         1 
                                         / 
                                         
                                           d 
                                           1 
                                           k 
                                         
                                       
                                       ) 
                                     
                                   
                                 
                               
                               
                                 
                                   sin 
                                   2 
                                 
                                 ( 
                                 
                                   ar 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   
                                     cos 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         1 
                                         / 
                                         
                                           d 
                                           2 
                                           k 
                                         
                                       
                                       ) 
                                     
                                   
                                 
                               
                             
                             ) 
                           
                           - 
                           
                             ( 
                             
                               
                                 
                                   a 
                                   
                                     1 
                                     ′ 
                                   
                                   k 
                                 
                                 - 
                                 
                                   a 
                                   0 
                                   k 
                                 
                               
                               
                                 
                                   a 
                                   
                                     2 
                                     ′ 
                                   
                                   k 
                                 
                                 - 
                                 
                                   a 
                                   0 
                                   k 
                                 
                               
                             
                             ) 
                           
                         
                         ) 
                       
                       2 
                     
                   
                 
               
               
                 
                   [ 
                   2 
                   ] 
                 
               
             
           
         
       
     
     In so-called AVO analyses, one can benefit from the AVO phenomenon (Amplitude Versus Offset). Knowledge of the amplitude of a reflected wave as a function of the incidence angle of the reflection makes it possible to extract richer information about the elastic properties of rocks on either side of an interface, which normal reflectivity alone does not furnish. 
     According to a possible embodiment mode of the invention, the traces A j+1  and A j+2  are corrected taking as reference trace the corrected trace of A j ′ of an offset immediately below trace A j+1  in the seismic traces collection. 
     The fourth step  40  consists of calculating a cost function F of the type: 
                   F   =       ∑     k   =   0     n     ⁢       (       (         sin   2     ⁡     (     ar   ⁢           ⁢     cos   ⁡     (     1   /     d     j   +   1     k       )         )           sin   2     ⁡     (     ar   ⁢           ⁢     cos   ⁡     (     1   /     d     j   +   2     k       )         )         )     -     (         a     j   +     1   ′       k     -     a     j   ′     k           a     j   +     2   ′       k     -     a     j   ′     k         )       )     2               [   3   ]               
wherein d j+1   k  and d j+2   k  are expansion coefficients associated with the k-th segments of traces A j+1  and A j+2 , a k   j ′, a k   j+1′  and a k   j+2 ′ are the amplitudes of the k-th segments of the corrected traces A j ′, A j+1 ′ and A j+2 ′.
 
     The cost function F estimates the similarity between the corrected traces A j+1 ′ and A j+2 ′ with the trace A j ′ being taken as reference trace. 
       FIG. 9  illustrates a third embodiment mode for the processing method whereby three traces are corrected simultaneously. In this variant, the expansion coefficients are determined simultaneously for three traces A 1 , A 2  and A 3 . 
     This third embodiment mode makes it possible to take into account the effect of anisotropy VTI (Vertical Transverse Isotropy) of the subsurface. 
     Calculation in the first section compares the three slanted times to propagation in the section k, T A1 , T A2  and T A3  with the vertical time T A0  and provides the three following parameters. 
     the vertical velocity ν k0 . 
     the Thomsen anisotropy parameters ε k  and δ k . According to a possible embodiment mode of the invention, the traces A j+1 , A j+2  and A j+3  are corrected taking as reference trace the corrected trace A j ′ of an offset immediately below the trace A j+1  in the seismic traces collection. 
     The offset correction method makes it possible to deduce information about the properties of the subsurface with greater resolution than with methods of prior art. In particular, it becomes possible to deduce the following parameters: 
     the local P wave velocity contrasts, 
     the density contrasts in the subsurface, 
     the S wave velocity contrasts, 
     the anisotropy parameters (vertical velocity ν P0  and ε k  and δ k  parameters). 
     1/Determination of Local P Wave Velocity Contrasts 
     As shown on  FIG. 10 , the subsurface is considered as formed of n parallel and isotropic sections, with vertical travel time ke, (k=1 . . . n) corresponding to the segmentation interval for the trace A 0  with zero slant. 
     One notes the angle of incidences i xkn  at the base of the section k for a wave emitted by a source at offset x reflecting at the base of section n. 
     According to a first step, the travel time of a wave is determined between the source at offset x and the base of the section of order n at normal incidence. 
     According to a second step, the average propagation velocity νp of the wave is determined between the source at offset x and the base of the section or order n at normal incidence. 
     Knowing the travel time between the source and each interface for each of the traces, it is possible to calculate the average propagation velocity of the wave for each trace. 
     For each section 1 to n, the angle i xnn  is determined by the relation: 
     
       
         
           
             
               
                 
                   
                     d 
                     xn 
                   
                   = 
                   
                     1 
                     
                       cos 
                       ⁡ 
                       
                         ( 
                         
                           i 
                           xnn 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   [ 
                   4 
                   ] 
                 
               
             
           
         
       
     
     The velocities and the angles of incidence satisfy the relation: 
     
       
         
           
             
               
                 
                   
                     
                       sin 
                       ⁡ 
                       
                         ( 
                         
                           i 
                           xnk 
                         
                         ) 
                       
                     
                     
                       v 
                       Pk 
                     
                   
                   = 
                   
                     
                       
                         sin 
                         ⁡ 
                         
                           ( 
                           
                             i 
                             xnq 
                           
                           ) 
                         
                       
                       
                         v 
                         Pq 
                       
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     for 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     any 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     k 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     and 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     q 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     sections 
                   
                 
               
               
                 
                   [ 
                   5 
                   ] 
                 
               
             
           
         
       
     
     According to a third step, a complete ray tracing is made between each source and the zero offset point with a depth: 
                   z   =     e   ⁡     (       ∑     k   =   0     n     ⁢     v   Pk       )               [   6   ]               
where  ν P k  satisfies the refraction law
 
     
       
         
           
             
               
                 v 
                 Pk 
               
               
                 v 
                 
                   Pk 
                   - 
                   1 
                 
               
             
             = 
             
               
                 
                   sin 
                   ⁡ 
                   
                     ( 
                     
                       i 
                       
                         1 
                         ⁢ 
                         nk 
                       
                     
                     ) 
                   
                 
                 
                   sin 
                   ⁡ 
                   
                     ( 
                     
                       i 
                       
                         
                           1 
                           ⁢ 
                           nk 
                         
                         + 
                         1 
                       
                     
                     ) 
                   
                 
               
               . 
             
           
         
       
     
     Thus, 
               dv   Pk       v   Pk           
is deduced for each section k, k=1 . . . n.
 
     According to a fourth step, the P wave is determined for each section k, k=1 . . . n from the relation: 
     
       
         
           
             
               
                 
                   
                     
                       dv 
                       Pk 
                     
                     
                       v 
                       Pk 
                     
                   
                   = 
                   
                     2 
                     ⁢ 
                     
                       
                         
                           v 
                           Pk 
                         
                         - 
                         
                           v 
                           
                             Pk 
                             - 
                             1 
                           
                         
                       
                       
                         
                           v 
                           Pk 
                         
                         + 
                         
                           v 
                           
                             Pk 
                             + 
                             1 
                           
                         
                       
                     
                   
                 
               
               
                 
                   [ 
                   7 
                   ] 
                 
               
             
           
         
       
     
     2/Determination Subsurface Density Contrasts 
     The density contrasts are obtained by subtracting the velocity contrasts from the velocity coefficients of reflection coefficients with incidence zero. Since it is known how to calibrate the seismic samples, it suffices to subtract the velocity contrasts from the zero incidence trace with a suitable calibration coefficient (F factor). 
     For each section 
             k   ,     k   =     1   ⁢           ⁢   …   ⁢           ⁢   n       ,       dv   Pk       v   Pk             
is known resulting from the sequence of P wave velocities determined in accordance with 1/.
 
     According to a first step, the density contrast 
               dv   ρk       v   ρk           
is determined for each section k, k=1 . . . n.
 
     In order to do this, a sequence of density contrasts 
                 d   ⁢           ⁢   ρ1     ρ1     ,       d   ⁢           ⁢   ρ2     ρ2     ,     …   ⁢           ⁢       d   ⁢           ⁢   ρ   ⁢           ⁢   n       ρ   ⁢           ⁢   n               
is defined and the sum
 
                 dv   Pk       v   Pk       +       d   ρk       ρ   ⁢           ⁢   k             
is calculated. The similarity between
 
                 dv   Pk       v   Pk       +       d   ρk       ρ   ⁢           ⁢   k             
and the zero offset A 0  trace is measured. Further density contrasts are swept to obtain maximum similarity.
 
     For example, the density contrasts 
                 d   ⁢           ⁢   ρ1     ρ1     ,       d   ⁢           ⁢   ρ2     ρ2     ,     …   ⁢           ⁢       d   ⁢           ⁢   ρ   ⁢           ⁢   n       ρ   ⁢           ⁢   n               
are determined by a random or pseudo-random selection algorithm. The density contrasts are taken initially from the predetermined ranges corresponding to orders of magnitude for the density contrasts generally observed. For example, initially it is possible to choose a series of density contrasts such that for any segment
 
     
       
         
           
             k 
             , 
             
               
                 
                   d 
                   ρk 
                 
                 
                   ρ 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   k 
                 
               
               = 
               
                 2 
                 , 
                 7. 
               
             
           
         
       
     
     The similarity is, for example, measured by calculating the correlation product: 
     
       
         
           
             
               
                 
                   
                     a 
                     0 
                     k 
                   
                   ⊗ 
                   
                     ( 
                     
                       
                         
                           dv 
                           Pk 
                         
                         
                           v 
                           Pk 
                         
                       
                       + 
                       
                         
                           d 
                           ρk 
                         
                         
                           ρ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           k 
                         
                       
                     
                     ) 
                   
                 
               
               
                 
                   [ 
                   8 
                   ] 
                 
               
             
           
         
       
     
     The sequence 
                 d   ⁢           ⁢   ρ1     ρ1     ,       d   ⁢           ⁢   ρ2     ρ2     ,     …   ⁢           ⁢       d   ⁢           ⁢   ρ   ⁢           ⁢   n       ρ   ⁢           ⁢   n               
is selected, leading to maximum similarity between
 
                 dv   Pk       v   Pk       +         d   ⁢               ρ   ⁢           ⁢   k         ρ   ⁢           ⁢   k             
for the trace A 0 .
 
     According to a second step, the factor F is educed as: 
                   F   =         (         dv   Pk       v   Pk       +       d   ⁢           ⁢     ρ   k         ρ   k         )     ⁢   moy       2   ×     (     a   0   k     )     ⁢   moy               (   9   )               
where  moy  designates the average value.
 
     The factor F is such that the products of F with the segments a 0   k  of the zero offset trace A 0  are equal to the reflection coefficients in the seismic band. 
     3/Determination of Local S Wave Velocity Contrasts 
     According to a first step, a gradient G of AVO amplitude is determined for the corrected seismic traces A 0 , A 1 , ‘A 2 ’, . . . A N ′ as a function of the angle of incidence of the wave at the base of section k. 
     For traces A 0 , A 1 , ‘A 2 ’, . . . A N ′  FIG. 11  is a diagram showing the amplitude a j   k ′ of the j-th corrected trace A j ′ as a function of sin 2 (i xjkn ), i xjkn  being the angle of incidence at the base of section k of a wave emitted by a source at the offset x j  reflecting on the base of the section n. 
     The amplitude a j   k  is a linear function of sin 2 (i xjkn ). From this, it is possible to deduce an estimate of the amplitude gradient G k  as the slope of a linear function. 
     Moreover, the amplitude gradient G k  in the section k verifies the following relation: 
                     v   Sk     =         v   Pk     2     ×           dv   Pk       v   Pk       -     2   ⁢     G   k             (       2   ⁢       dv   Sk       v   Sk         +         d   ⁢               ρ   ⁢           ⁢   k         ρ   ⁢           ⁢   k         )       1   /   2                   (   10   )               
where ν pk  is the P wave velocity, Pk is the acoustic impedance of the subsurface and ν Sk  is the S wave velocity, in section k.
 
     
       
         
           
             
               
                 dv 
                 Pk 
               
               
                 v 
                 Pk 
               
             
             = 
             
               2 
               ⁢ 
               
                 
                   
                     v 
                     Pk 
                   
                   - 
                   
                     v 
                     
                       Pk 
                       - 
                       1 
                     
                   
                 
                 
                   
                     v 
                     Pk 
                   
                   + 
                   
                     v 
                     
                       Pk 
                       + 
                       1 
                     
                   
                 
               
             
           
         
       
       
         
           
             
               
                 dv 
                 Sk 
               
               
                 v 
                 Sk 
               
             
             = 
             
               2 
               ⁢ 
               
                 
                   
                     v 
                     Sk 
                   
                   - 
                   
                     v 
                     
                       Sk 
                       - 
                       1 
                     
                   
                 
                 
                   
                     v 
                     Sk 
                   
                   + 
                   
                     v 
                     
                       Sk 
                       + 
                       1 
                     
                   
                 
               
             
           
         
       
     
     The relation (10) is equivalent to: 
     
       
         
           
             
               
                 
                   
                     2 
                     ⁢ 
                     Gk 
                   
                   = 
                   
                     
                       
                         dv 
                         Pk 
                       
                       
                         v 
                         Pk 
                       
                     
                     × 
                     2 
                     ⁢ 
                     
                       
                         ( 
                         
                           
                             v 
                             Pk 
                           
                           
                             v 
                             Sk 
                           
                         
                         ) 
                       
                       2 
                     
                     × 
                     
                       ( 
                       
                         
                           2 
                           ⁢ 
                           
                             
                               dv 
                               Sk 
                             
                             
                               v 
                               Sk 
                             
                           
                         
                         + 
                         
                           
                             
                               d 
                               ⁢ 
                               
                                   
                               
                             
                             
                               ρ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               k 
                             
                           
                           
                             ρ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             k 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
           
         
       
     
     According to a second step, knowing the gradients G k , the P wave velocity contrasts 
               dv   Pk       v   Pk           
and the impedance contrasts
 
                 d   ⁢               ρ   ⁢           ⁢   k         ρ   ⁢           ⁢   k           
in all the sections,
 
               dv   Sk       v   Sk           
is deduced for each section k using the relation (11).
 
     In order to do this, a sequence of velocities S ν S1 , ν S2 , . . . ν Sn  is defined and the gradients are determined for k=1 . . . n. The similarity between the gradient G k  obtained by the relation (11) and the gradient G k  estimated during the first step is measured. The S wave velocity sequences are swept to obtain maximum similarity between the gradients. 
     The S wave velocities ν S1 , ν S2 , . . . ν Sn  are, for example, determined by a random or pseudo-random selection algorithm. The S wave velocities are taken initially within predetermined ranges corresponding to the orders of usual S wave velocity magnitudes. For example, initially it is possible to choose a series of S wave velocities such as 
     
       
         
           
             
               v 
               Sk 
             
             = 
             
               
                 
                   v 
                   Pk 
                 
                 2 
               
               . 
             
           
         
       
     
     The similarity is measured, for example, by calculating the correlation product between gradients. 
     The sequence of S wave velocities, ν S1 , ν S2 , . . . ν Sn  is determined, leading to maximum similarity between gradients. 
     4/Determination of Anisotropy Parameters in P Mode 
     Now it is assumed that the subsurface shows axial anisotropy around a vertical or anisotropic VTI (Vertical Transverse Isotropy) axis of symmetry. The angles of incidence i xnk  of the waves no longer result directly from the expansions d 1   1 , d 1   2 , d 1   3 , . . . d 1   n . 
     The path of a seismic ray is not perpendicular to the wave front. The velocity all along a ray ν r  and the normal velocity at the wave front ν h  named the phase velocity can thus be distinguished. 
     In P mode, these velocities depend on the angle r between the ray and the axis of symmetry, or the angle h between the perpendicular to the wave front and the axis of symmetry. Thus one has: 
                     v   r     =       v     r   ⁢           ⁢   0       ×     (     1   +     δ   ⁢           ⁢     sin   2     ⁢   r   ⁢           ⁢     cos   2     ⁢   r     +     ɛ   ⁢           ⁢     sin   4     ⁢   r       )               (   12   )               
where ν r0  is the velocity of a ray at zero incidence.
 
                     v   h     =       v     h   ⁢           ⁢   0       ×     (     1   +     δ   ⁢           ⁢     sin   2     ⁢   r   ⁢           ⁢     cos   2     ⁢   h     +     ɛ   ⁢           ⁢     sin   4     ⁢   h       )               (   13   )               
where ν h0  is the velocity of the perpendicular to the wave front at zero incidence. It is to be noted that, at zero incidence, the velocity all along a ray and the phase velocity are equal (with ν r0 =ν h0 ).
 
     Descartes law applies to the phase velocities ν h . This is why the ray tracing step is divided into sub-steps. 
     According to a first sub-step, ν hn  is determined from ν rn  for section n. 
     According to a second sub-step, ν hn+1  is calculated from Descartes law: 
     
       
         
           
             
               
                 
                   
                     
                       sin 
                       ⁡ 
                       
                         ( 
                         
                           i 
                           nn 
                         
                         ) 
                       
                     
                     
                       v 
                       
                         h 
                         n 
                       
                     
                   
                   = 
                   
                     
                       sin 
                       ⁡ 
                       
                         ( 
                         
                           i 
                           
                             nn 
                             + 
                             1 
                           
                         
                         ) 
                       
                     
                     
                       v 
                       
                         h 
                         
                           n 
                           + 
                           1 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
           
         
       
     
     According to a third sub-step, ν rn+1  is determined from ν hn+1  by using the relation (15). 
     The fourth sub-step consists of determining the time or distance of propagation in the section n+1: 
     
       
         
           
             
               
                 
                   
                     tan 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     h 
                   
                   = 
                   
                     
                       tan 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       r 
                     
                     
                       1 
                       + 
                       
                         2 
                         ⁢ 
                         δ 
                       
                       + 
                       
                         4 
                         ⁢ 
                         
                           ( 
                           
                             ɛ 
                             - 
                             δ 
                           
                           ) 
                         
                         ⁢ 
                         
                           sin 
                           2 
                         
                         ⁢ 
                         r 
                       
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     The ray tracing step is applied simultaneously on at least three traces A k , A k+1 , A k+2  to determine ν Pk  at zero incidence x=0, ε k  and δ k . 
     The angles of incidence i xnk  of the waves no longer depend simply on the expansions d 1   1 , d 1   2 , d 1   3 , . . . d 1   n . One has: 
     
       
         
           
             
               
                 
                   
                     d 
                     j 
                     k 
                   
                   = 
                   
                     
                       v 
                       
                         P 
                         
                           k 
                           - 
                           1 
                         
                       
                     
                     
                       
                         v 
                         
                           P 
                           k 
                         
                       
                       × 
                       
                         ( 
                         
                           
                             cos 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               r 
                               k 
                             
                           
                           + 
                           
                             sin 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               r 
                               k 
                             
                             × 
                             
                               tan 
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     r 
                                     k 
                                   
                                   - 
                                   
                                     h 
                                     k 
                                   
                                 
                                 ) 
                               
                             
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   16 
                   ) 
                 
               
             
           
         
       
     
     The same method for determining the anisotropy parameters can be applied to the case of a subsurface having axial anisotropy around an inclined axis or TTI anisotropy (Tilted Transverse Isotropy). In this case, it is necessary to take into account the dip of the sections.