Abstract:
Conventional migration of short offset seismic data is performed. An interval velocity is obtained using, for example, check shots. An initial model of effective anellipticity parameter as a function of depth is obtained by flattening long offsets within a common image point. From these, interval anellipticity and horizontal velocity are obtained as a function of time. These initial models are used for anisotropic imaging. Any remaining residuals are used to update the interval anellipticity model and the process of migration is repeated until satisfactory results are obtained.

Description:
BACKGROUND OF THE INVENTION  
         [0001]    1. Field of the Invention  
           [0002]    The present invention is a method of obtaining an anisotropic velocity model for proper depth imaging of seismic data.  
           [0003]    2. Related Prior Art  
           [0004]    The search for subsurface hydrocarbon deposits typically involves a multifaceted sequence of data acquisition, analysis, and interpretation procedures. The data acquisition phase involves use of an energy source to generate signals which propagate into the earth and reflect from various subsurface geologic structures. The reflected signals are recorded by a multitude of receivers on or near the surface of the earth, or in an overlying body of water. The received signals, which are often referred to as seismic traces, consist of amplitudes of acoustic energy which vary as a function of time, receiver position, and source position and, most importantly, vary as a function of the physical properties of the structures from which the signals reflect. The data analyst uses these traces along with a geophysical model to develop an image of the subsurface geologic structures.  
           [0005]    Common Mid Point (CMP) stacking, also sometimes referred to as Common Depth Point or Common Reflection Point, CDP or CRP respectively) of seismic field data is well known. See for example U.S. Pat. No. 3,217,828 to Mendenhall et al, and U.S. Pat. No. 2,732,906 to Mayne, which are incorporated herein by reference as a teaching of the CMP technique. The U.S. Pat. No. 3,217,828 teaches two-dimensional (2-D) data processing where dipping earth layers are projected into a two-dimensional plane, perpendicular to the surface of the earth, along a designated line of profile.  
           [0006]    Seismic data migration typically uses diffraction traveltimes from subsurface imaging points to the source and receiver locations to produce an image of the subsurface reflectors. The diffraction traveltimes are the seismic signal propagation times along raypaths from each imaging point to the source and receiver locations. The propagation times, which are usually plotted as diffraction traveltime curves, are used after appropriate preprocessing of the raw seismic data to generate an estimate of the correct location of the reflector.  
           [0007]    The basic concepts of migration are discussed next with reference to FIG. 1. A portion of the surface of the earth is shown at  10 . A sound source S may be offset from a receiver R by an offset distance 2 k with a midpoint at m 0 . A wavefield generated at S, may travel along ray path  12  to a reflecting point  14  on reflector  16  and reflected back along ray path  18  to receiver R. The amplitude of the wave field as a function of time, f(t), may be recorded on a time-scale recording, hereinafter referred to as trace. Seismic data are generally quantized as digital samples. In the processes next to be discussed, each sample of a seismic trace may be operated upon individually. In the interest of brevity, use of the collective term “seismic trace” in conjunction with a process means that every data sample of that trace has been individually subjected to the named process.  
           [0008]    Because of the offset distance 2 k, the arrival time t of a particular wavelet along path S- 14 -R is greater than the travel time t 0  of a wavelet that might have traveled along a direct path m 0 - 14 -m 0 . The time difference is termed normal moveout (NMO). The quantity f(t) may be corrected for NMO by the relation  
               t   0   2     =       t   k   2     -     (     4          k   2     /     V   2         )               (   1   )                               
 
           [0009]    where V represents the NMO velocity of the medium through which the wavelet traveled. The above relation is a small-offset approximation. It should be observed that for zero offset k and constant velocity V the reflection points for f(t 0 ) such as  14  lie along a semicircle such as  19  centered about m 0 . Wavelet amplitude varies according to the inverse square law due to geometric spreading.  
           [0010]    The hyperbolic time-distance relationship given by eq. (1) is illustrated in FIG. 2 where a plurality of traces  21   a ,  21   b , . . .  21   n  are shown with a reflection event depicted by the curve  25 . The general process of migration involves moving portions of seismic data along a diffraction curve such as  25  to an output location and summing the contributions from a large portion of the seismic section. The reflector in the subsurface may be viewed as being made up of a plurality of diffraction points. At the correct position of a diffractor, contributions from the plurality of migrated traces will sum coherently whereas at other locations, the summation will be incoherent and have close to zero amplitude; in this fashion an image of the subsurface can be reconstructed.  
           [0011]    The migration process may be implemented in the time domain by the summation process or it may be implemented in the frequency domain. Such methods would be known to those versed in the art.  
           [0012]    There is a vast array of literature indicating that the velocity of propagation of seismic waves in the subsurface is anisotropic, i.e., it depends upon the direction of propagation. There are numerous causes of such anisotropy. Backus showed that layering of isotropic materials having different shear modulii and with a layer thickness much shorter than the wavelength of seismic waves produces a transverse isotropy with a vertical symmetry axis (referred to hereafter VTI). In addition, shales exhibit transverse isotropy in their elastic properties even at very short seismic wavelengths. This has been attributed to intrinsic anisotropy. The cause of VTI is not pertinent to the present invention in that the present invention is a method for properly migrating seismic data in the presence of VTI regardless of the origin of the VTI.  
           [0013]    A TI medium is characterized by five elastic modulii. These may be denoted by the tensor  
             [           c   11             c   11     -     2        c   66               c   13         0       0       0               c   11     -     2        c   66               c   11           c   13         0       0       0             c   13           c   13           c   33         0       0       0           0       0       0         c   44         0       0           0       0       0       0         c   44         0           0       0       0       0       0         c   66           ]           (   2   )                               
 
           [0014]    where the modulus c 11  defines the velocity of a horizontally propagating P-wave, c 33  defines the velocity of a vertically propagating P-wave, c 44  defines the velocity of a vertically propagating shear wave, and c 66  defines the velocity of a horizontally propagating S-wave (shear wave) with horizontal polarization. These four parameters are determinable by making suitable measurements of P- and S-waves parallel to and perpendicular to the symmetry axis. The velocity is given by the square root of the ratio of the elastic modulus to the density. However, the parameter c 13  is not readily determinable by measurements parallel to and perpendicular to the symmetry axis. Determination of c 13  requires a measurement off the symmetry axis. It also turns out that this parameter c 13  is of particular importance in reflection seismic prospecting.  
           [0015]    Thomsen defined several parameters characterizing a VTI medium. One of these is the quantity ε defined by  
             ɛ   =         c   11     -     c   33         2        c   33                 (   3   )                               
 
           [0016]    This parameter is directly related to the difference between the horizontal and vertical P-wave velocities. Another parameter defined by Thomsen is the quantity δ defined as  
             δ   =           (       c   13     +     c   44       )     2     -       (       c   33     -     c   44       )     2         2          c   33          (       c   33     -     c   44       )                   (   4   )                               
 
           [0017]    Based on the results derived by Thomsen, for a horizontal reflector in a TI medium, eq.(1) is valid for small offsets with 
             t   n   =z/V   z   (5) 
           and 
             V   nmo   =V=V   z {square root}{square root over (1+2δ)}  (6) 
           [0018]    where z is the depth of the reflector, V z  is the vertical velocity. The normal moveout velocity is thus dependent upon the parameter δ and the vertical velocity.  
           [0019]    Alkhalifah &amp; Tsvankin showed that the P-wave reflection moveout at large offsets is largely controlled by just two parameters. One is the V nmo  given by eq. (6) and the second is a parameter η defined by  
             η   =         1   2          (         V   hor   2       V   nmo   2       -   1     )       =       ɛ   -   δ       1   +     2      δ                   (   7   )                               
 
           [0020]    where 
             V   hor   =V   nmo {square root}{square root over (1+2η)}  (8) 
           [0021]    with the moveout equation given by  
                       t   2          (   x   )       =              t   o   2     +       x   2       V   nmo   2       -       2                 η                   x   4           V   nmo   2     [         t   0   2          V   nmo   2       +       (     1   +     2                 η       )          x   2         ]                     =              t   0   2     +       x   2       V   nmo   2       -         (       V   hor   2     -     V   nmo   2       )          x   4           V   nmo   2     (         t   2          V   nmo   4       +       V   ho   2          x   2         )                       (   9   )                               
 
           [0022]    and x=2 k in eq. (1) and FIG. 1.  
           [0023]    As would be known to those versed in the art, conventional velocity analysis of seismic reflection data is based on use of semblance scans along hyperbolic moveout curves. Gretchka &amp; Tsvankin disclose a method for inversion of seismic data using semblance scans along non-hyperbolic moveout curves using a modified version of eq. (9). The modification is based on empirical observations for simulated results wherein an offset range of x/z≅2 was used. The modified equation used by Gretchka &amp; Tsvankin is  
                 t   2          (   x   )       =       t   0   2     +       x   2       V   nmo   2       -         (       V   hor   2     -     V   nmo   2       )          x   4           V   nmo   2     (         t   0   2          V   nmo   4       +       CV   hor   2          x   2         )                 (   10   )                               
 
           [0024]    the modification being the empirical constant C in the denominator of the last term on the right hand side of eq. (10).  
           [0025]    Accuracy of the subsurface image is of great importance in seismic prospecting: interpretation of the structure and lithology of the subsurface is preferably done on depth migrated sections. The accuracy of the image is closely related to the velocity model. Gretchka &amp; Tsvankin do not address the interaction between iterative changes in the velocity model and subsequent changes in the accuracy of the depth image.  
           [0026]    It is therefore desirable to have a method for determination and iterative updating of subsurface velocities in anisotropic earth formations. Such a method should preferably be efficient in terms of computer resources and should also preferably provide more accurate images of the earth formations. The present invention satisfies this need.  
         SUMMARY OF THE INVENTION  
         [0027]    The present invention is a method of seismic prospecting over a region of interest having velocity anisotropy. Conventionally migrated seismic data are obtained using an isotropic velocity model wherein short offsets within a common image point gather have substantially flat moveouts. A model of vertical interval velocity as a function of depth is obtained using checkshots, VSPs, sonic logs, or geologic information. An initial model of an effective anellipticity parameter as a function of depth over the region of interest is obtained. Using the conventionally migrated seismic data, the model of vertical interval velocity as a function of depth and the anellipticity parameter as a function of depth, an anisotropic migrated image of the subsurface is obtained. By reviewing long offsets within common image points of the anisotropic migrated image, a residual anellipticity parameter is determined, the subsurface model is updated, and the migration is repeated. The anisotropic migrated image is preferably obtained using Kirchoff migration. The conventionally migrated seismic data may be obtained by Kirchoff migration, finite difference, or frequency domain methods.  
           [0028]    Determination of the effective anellipticity parameter further comprises using semblance processing with non-hyperbolic moveout.  
           [0029]    The Kirchoff migration is preferably carried out after a depth-to-time conversion of the effective anellipticity model. That is, following a layer stripping, models of the horizontal interval velocity and the anellipticity parameter are obtained as a function of time: this is followed by a conversion back to depth prior to the anisotropic imaging. 
       
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0030]    The file of this patent contains at least one drawing executed in color: Copies of this patent with color drawing(s) will be provided by the Patent and Trademark Office upon request and payment of the necessary fee. The present invention is best understood with reference to the following figures in which like numerals refer to like elements and in which:  
         [0031]    [0031]FIG. 1 illustrates the basic geometry used in seismic prospecting of subsurface formations;  
         [0032]    [0032]FIG. 2 is an illustration of seismic traces;  
         [0033]    [0033]FIG. 3 is a flow chart illustrating some of the steps of the present invention;  
         [0034]    [0034]FIG. 4 is an example of a seismic section;  
         [0035]    [0035]FIG. 5 is an exemplary initial display showing a conventionally migrated seismic section along with diagnostic displays of a common image point gather and an effective anellipticity parameter;  
         [0036]    [0036]FIGS. 6 a ,  6   b  and  6   c  show displays of initial models of the vertical velocity, horizontal velocity and anellipticity parameter as a function of depth;  
         [0037]    [0037]FIG. 7 shows the final display of the seismic section, the common image point gather and the effective residual anellipticity parameter; and  
         [0038]    [0038]FIG. 8 shows a comparison of conventional migration with a migration using the method of the present invention. 
     
    
     DETAILED DESCRIPTION OF THE INVENTION  
       [0039]    Turning now to FIG. 3, a flow chart that includes the preferred embodiment of the present invention is shown. Anisotropic depth migration requires three independent models: a vertical velocity model V z , a horizontal velocity model V x , and an anisotropy parameter model η satisfy this requirement. Seismic data over the region of interest are acquired or obtained  101  using any one of numerous well-known seismic exploration techniques. This gives seismic traces such as those depicted in FIG. 2 wherein one or more seismic reflections (or diffractions) appear on a plurality of seismic traces. Each seismic trace is associated with a source position and a receiver position. Shown in FIG. 4 is an example of a seismic section. It is to be noted that in actual practice, color displays are used for displays of seismic sections and diagnostics: this greatly increases the ease with which desired features can be seen. Black and white displays, as shown in this document, are usually adequate, but color displays are preferable.  
         [0040]    Still referring to FIG. 3, the data acquired at  101  are then subjected to routine seismic processing in well-known ways. The following techniques exemplify such processing. As part of the routine seismic processing, the data may have the source signature removed, frequency filtered, resampled, and then sorted to common midpoint gathers. Next the gathers might be subjected to multichannel filters to remove noise and multiples. This may be done for all seismograms that were recorded.  
         [0041]    Next, an isotropic velocity model (V m ) is built and the processed seismograms are prestack depth migrated yielding a depth migration dataset  105 . Prestack depth migration uses a velocity model of the subsurface to reposition the seismograms that are recorded into time to proper locations in depth. This initial velocity model may be based upon analysis of well logs, stacking velocities, or tomographic velocities. The geophysicist can then interpret the depth data to locate hydrocarbons or similar valuable minerals. There are several known ways to depth migrate the data including Kirchhoff summation, finite difference methods, or frequency finite difference methods. All rely on an accurate velocity model. In a preferred embodiment of the invention, a Kirchhoff summation is used. Given a velocity model a traveltime is calculated from a source position to a position in the subsurface known as a common image point. Another traveltime is calculated from the receiver position to the common image point. The amplitude from the seismogram corresponding to this source-receiver combination at the resulting time is summed into this common image point. Usually a similar offset is summed into one depth dataset. The dataset is then sorted by position on the earth and offset resulting in a common image point gather (CIP gather). An example of a CIP gather is shown in the top right hand panel of FIG. 5. The source-receiver offset increases from left to right. The location corresponding to CIP panel  201  is also indicated on the upper left hand panel. The bottom right panel of FIG. 5 shows an enlarged portion of the display in the top right panel  
         [0042]    These CIP gathers are then analyzed. If the velocity model is correct the CIP gathers will be flat and the gathers can be stacked. If the CIP gathers are not flat the velocity model is refined using well-known techniques such as vertical updating or ray tomography. The data is then prestack depth migrated and the process repeated until the CIP gathers are flat. In an isotropic medium such as sandstone the horizontal and vertical velocities are the same. In VTI media such as shale (and in finely laminated sequences of isotropic materials) the horizontal velocities are faster than the vertical velocities. Using this isotropic velocity model the seismic traces at small offsets in the CIP gathers will be flat using this technique, but the far offsets which result from energy passing through earth formations at angles closer to horizontal will be influenced more by the horizontal velocity. Usually, the far offsets in the CIP will not be flat in the presence of anisotropy. It can be seen in FIG. 5, that the near offsets in the CIP are flat, but the far offsets are not. For convenience, data obtained at this point in the processing (i.e., after  105 ) is referred to as conventionally migrated seismic data. In addition, data at the output of  103  is referred to as conventionally processed seismic data.  
         [0043]    Still referring to FIG. 3, a vertical interval velocity V z  model as a function of depth is created  107 . This model represents interval velocity in the vertical direction as a function of depth. This velocity information can be obtained from well information such as checkshot surveys, sonic logs, or known tops of geologic units in the well. Data may also be obtained from a Vertical Seismic Profile (VSP). Calculation from well data of V z  as a function of depth uses methods well-known in prior art. It is to be pointed out that Gretchka &amp; Tsvankin do not discuss obtaining a model of the vertical interval velocity as a function of depth.  
         [0044]    Next, an η model is created  109  by analyzing and flattening the far offsets in the CIP gathers from  105 . A semblance analysis for η is calculated on the CIP gathers using the moveout equation:  
               t   2     =       t   0   2     +         x   2       V     m   ,   eff     2       [     1   -     κ        (     x   ,       t   0     ;     V     m   ,   eff         ,     η   eff       )         ]               (   11   )                               
 
           where                         κ   ≡       2                   η   eff          x   2             t   0   2          V   m   2       +       (     1   +     2                   η   eff         )          x   2                   (   12   )                               
 
         [0045]    where C is the same constant as in eq. (10), and t=CIP travel time  
         [0046]    t 0 =two way vertical travel time  
         [0047]    x=offset distance  
         [0048]    V m,eff =effective velocity from short-offset CIP analysis (calculated from  105 )  
         [0049]    η eff =effective anisotropy (anellipticity) coefficient.  
         [0050]    The maximum value of semblance is picked. This gives a model of the effective value of η eff  as a function of time that is converted to a model of the effective velocity of η eff  as a function of depth, using the V z  model as a function of depth from  107 .  
         [0051]    An example of η analysis is shown in the bottom left panel of FIG. 5. This corresponds to a region close to the horizon  203 . If there is no anisotropy in the area η eff  will be zero. The semblance curve is picked to yield the best η eff  value that flattens the CIP gather. The η eff  function may vary laterally and vertically. This analysis is done for several layers in the seismic survey and a model η eff  is produced.  
         [0052]    Proceeding to  111 , an effective velocity model V m,eff  is created in time from the isotropic V m -depth model created in  105 . Using the equation:  
                 V     m   ,   eff     2          (     t   0     )       =         ∫   0     t   0                 τ                       V   m   2          (   τ   )               ∫   0     t   0               τ                 (   13   )                               
 
         [0053]    where V m (t 0 )=short-offset interval velocity. The η eff -depth model from  109  is also stretched to time using the isotropic velocities calculated in  105 .  
         [0054]    Next, an interval η-time model and an interval horizontal velocity V x -time model is created  113  using equations from Grechka and Tsvankin: 
           V   x,eff ( t   0 )= V   m,eff ( t   0 ){square root}{square root over (1+2η eff ( t   0 ))}  (14) 
           f ( t   0 )= V   m,eff   2 ( t   0 )[4 V   x,eff   2 ( t   0 )−3 V   m,eff   2 ( t   0 )]  (15) 
         [0055]    [0055]               η                   (     t   0     )       =       1   8          [             t   0          f        (     t   0     )         -       (       t   0     -     δ                   t   0         )          f        (       t   0     -     δ                   t   0         )                 V   m   4          (     t   0     )                     δ                   t   0         -   1     ]               (   16   )                                                 V   x          (     t   0     )       =       1   2            V   m          (     t   0     )                      t   0          f        (     t   0     )         -       (       t   0     -     δ                   t   0         )          f        (       t   0     -     δ                   t   0         )                 V   m   4          (     t   0     )                     δ                   t   0         +   3                 (   17   )                                 
         [0056]    or equivalently 
           V   x ( t   0 )= V   m ( t   0 ){square root}{square root over (1+2η( t   0 ))}  (18) 
         [0057]    where  
         [0058]    V x,eff (t 0 )=effective horizontal velocity  
         [0059]    V x (t 0 )=interval horizontal velocity  
         [0060]    η(t 0 )=interval anisotropy (anellipticity) parameter  
         [0061]    δt 0 =time sample rate of model.  
         [0062]    It should be pointed out that portions of the invention disclosed hereafter are not part of the teachings of Grechka and Tsvankin.  
         [0063]    Proceeding to  115  the interval η-time model and the interval V x -time model are stretched to true depth using the V z  model created in  107 . The models after this are now in the same depth domain. FIGS. 6 a - 6   c  show examples of the three models. FIG. 6 a  is a display of the vertical velocity V z  as a function of depth (vertical axis) and horizontal location along the seismic line; FIG. 6 b  is a similar display of the horizontal velocity while FIG. 6 c  is a display of the anellipticity parameter η. It can be seen that in the shallow portion of the model (which is dominated by shale) the horizontal velocity is faster than the vertical velocity.  
         [0064]    Next, the three models V z , V x , and η are used in depth migration  117 . In a preferred embodiment of the invention, Kirchhoff migration is used, though other methods of migration could also be used. Pre-stack migration techniques are discussed, for example in Johnson (U.S. Pat. No. 4,964,103) having the same assignee as the present application and the contents of which are fully incorporated herein by reference. Kirchoff migration is discussed, for example in Gray &amp; May and in Sena &amp; Toksoz. In the present invention, the depth imaging travel times are calculated using all three models (V z , V x , and η). The travel time calculation used in this method is an eikonal type, but raytracing and other methods could also be used. The Kirchhoff migration results in a set of CIP gathers similar to those in  105 .  
         [0065]    The CIP gathers from  117  are then analyzed  119  by semblance calculation as in  109 . The seismic traces in the CIP gather will be flat at small offsets. The far offsets in the gathers are flat if the η model is correct. If the gathers are flat the method proceeds to  121  where the CIP gathers are then summed and the imaging complete. If the gathers are not flat the method proceeds to  122 . At this point, it is evident that the model is incorrect. The data and models are stretched back to time  122 . The semblances are re-computed and a residual η-time model is created.  
         [0066]    At  125  the η anisotropy model is then updated with the residual η-time model using one of the following two procedures. In one embodiment, use is made of the following equations:  
                   η                  eff   ,   new                       (     t   0     )       =           κ        (       η                  eff   ,   old       )       -   r       1   -     (       κ        (       η                  eff   ,   old       )       -   r     )                    V     m   ,   eff     2          (     t   0     )              t   h   2          (     x   ref     )           2                   x   ref   2                   (   19   )                 t   h   2     ≡       t   0   2     +       x   2         V     m   ,   eff     2          (     t   0     )                   (   20   )                 r        (       Δ                 t     ;   x     )       ≡     Δ                   t        (       Δ                 t     +     2                   t   0         )                V     m   ,   eff     2          (     t   0     )         x   2                 (   21   )                               
 
         [0067]    where Δt is measured directly from data, as shown in FIG. 5, and updated η eff,new  is determined from eqs. (19-21). In eqs. (19-21)  
         [0068]    η eff,new =updated effective anisotropy parameter  
         [0069]    η eff,old =previous effective anisotropy parameter  
         [0070]    x ref =distance from the source to a reference offset  
         [0071]    Alternatively, eq. (22) is fit directly to data to estimate an approximate η eff,res  and η eff,new  from eq. (23).  
                 t   res   2          (   x   )       ≈       t   0   2     -     2                   η     eff   ,   res              x   4         t   0   2          V     m   ,   eff     4                     (   22   )                               
 η eff,res ≡η eff,new −η eff,old   (23) 
         [0072]    where  
         [0073]    η eff,res =residual effective anisotropy parameter  
         [0074]    x=distance from the source to each receiver  
         [0075]    The method then proceeds to  113  where the interval η-time model and V x -time models are calculated. The seismic data are anisotropic depth migrated and the CIP gathers analyzed. Once the gathers are flat the method proceeds to  121  where the CIP gathers are stacked and the imaging is complete. A properly migrated CIP gather is shown in FIG. 7. This is the same CIP gather as shown in FIG. 5 with only the isotropic velocity model. Comparison of the upper and lower right panels in FIG. 7 with the corresponding panels in FIG. 5 shows considerable improvement in the CIP gather. The semblance panel in the lower left comer of FIG. 7 indicates little or no residual correction to the η model is required after this iteration. By properly representing the earth and building exact anisotropic velocity models the data are properly positioned.  
         [0076]    A completed seismic anisotropic depth migrated image is shown in FIG. 8 b  with the corresponding isotropic depth migrated image FIG. 8 a . Note that the reflections and fault are sharper in the anisotropic depth image. The seismic image thus produced yields a more exact representation of the earth&#39;s subsurface, tying wells better, with structures such as faults and channels in their proper position. Even to the untrained eye, the difference in resolution of the major fault (compare  301  with  301 ′) should be apparent. In addition, the shallow channel  303 ′ in the anisotropic migration appears to be better defined than the corresponding image  303  in the isotropic migration. Other differences are seen by the trained seismic interpreter.  
         [0077]    While the foregoing disclosure is directed to the preferred embodiments of the invention, various modifications will be apparent to those skilled in the art. It is intended that all variations within the scope and spirit of the appended claims be embraced by the foregoing disclosure.