Abstract:
A system and method for subsurface characterization including depth and structural uncertainty estimation is disclosed. In one embodiment, the method may include determining a detectability threshold for moveout in a seismic data gather based on the seismic data and computing a depth uncertainty function, wherein the depth uncertainty function represents an error estimate that is used to analyze an interpretation of the seismic data. In another embodiment, the method may include receiving a depth uncertainty volume and at least one interpreted horizon from seismic data, extracting a depth uncertainty cage for each of the interpreted horizons based on the depth uncertainty volume, and simulating multiple realizations for each of the interpreted horizons, constrained by the depth uncertainty cage. The multiple realizations may be used for analyzing changes to geometrical or structural properties of the at least one interpreted horizon. The changes may be plotted as at least one distribution and may be used to make P10, P50 and P90 estimates.

Description:
FIELD OF THE INVENTION 
     The present invention relates generally to methods and systems for characterization of subsurface formations using seismic data and, in particular, methods and systems for subsurface characterization that estimate uncertainty in the depth or structure of the interpreted subsurface based on an anisotropic seismic velocity model. 
     BACKGROUND OF THE INVENTION 
     Oil and gas prospect evaluation and field development require accurate characterization of subsurface features. Seismic acquisition over subsurface structure generally produces time-domain data, which is then migrated to, for example, depth image data. The migration process necessarily involves certain assumptions regarding the propagation velocity of elastic waves through the subsurface materials and structures. Moreover, there is generally some degree of anisotropy in geological formations. That is, while it may be possible to determine vertical velocities using well data, the velocities estimated using multi-offset seismic techniques will necessarily be somewhat different from measured vertical velocities. Finally, because assumptions, based on measurement or estimation, regarding both velocities and degrees of anisotropy may be incorrect, there is some inherent uncertainty in the resulting depth image, both in the depth of imaged events and in the structural interpretation of the events. 
     Models of the subsurface obtained from geophysical measurements are inherently non-unique. Geophysical measurements are finite in resolution, and relate to many orders of magnitude of scale. Uncertainty in the measurements results from a variety of sources, including signal-to-noise ratio, data acquisition parameter selection, processing algorithms, or the above mentioned velocity and anisotropy parameter selection. It is therefore important to understand the degree of that uncertainty when evaluating model results. That is, it is important to quantitatively understand to what degree the models are sensitive to a given change or group of changes in the assumptions regarding velocities, anisotropy or the other factors impacting uncertainty. An understanding of the uncertainty and the range of possible characterizations allows interpreters of the data to make business decisions regarding reserve estimation, well placement and count, development scenarios, secondary recovery strategies and other factors that ultimately impact recovery and project economics. 
     SUMMARY OF THE INVENTION 
     Described herein are implementations of various approaches for a computer-implemented method for subsurface characterization from seismic data. In one embodiment, the method may include receiving effective parameters and properties from an anisotropic seismic velocity model, at least one velocity function, and key descriptors from seismic data; determining a detectability threshold for moveout in a seismic data gather based on the seismic data; and computing a depth uncertainty function based on the detectability threshold, the effective parameters and properties, and the at least one velocity function, wherein the depth uncertainty function represents an error estimate that is used to analyze an interpretation of the seismic data. 
     In another embodiment, the computer-implemented method may include receiving a depth uncertainty volume and at least one interpreted horizon from seismic data; extracting a depth uncertainty cage for each of the at least one interpreted horizons based on the depth uncertainty volume; and simulating multiple realizations for each of the at least one interpreted horizons, constrained by the depth uncertainty cage. The multiple realizations may be used for analyzing changes to geometrical or structural properties of the at least one interpreted horizon. The changes may be plotted as at least one distribution and may be used to make P10, P50 and P90 estimates. 
     Additional embodiments include systems and articles of manufacture configured to implement the methods described herein. 
     The above summary section is provided to introduce a selection of concepts in a simplified form that are further described below in the detailed description section. The summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject matter. Furthermore, the claimed subject matter is not limited to implementations that solve any or all disadvantages noted in any part of this disclosure. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       These and other features of the present invention will become better understood with regard to the following description, claims and accompanying drawings where: 
         FIG. 1  is a flowchart illustrating a method in accordance with an embodiment of the present invention; 
         FIG. 2  is a flowchart illustrating a method for using detectability in accordance with an embodiment of the present invention; 
         FIG. 3  illustrates exemplary moveout curves computed by a method in accordance with the present invention; 
         FIG. 4  illustrates exemplary moveout curves including wavelet thickness computed by a method in accordance with the present invention; 
         FIG. 5  illustrates exemplary velocity curves computed by a method in accordance with the present invention; 
         FIG. 6  illustrates an exemplary depth uncertainty curve computed by a method in accordance with the present invention; 
         FIG. 7  is a flowchart illustrating a method in accordance with an embodiment of the present invention; and 
         FIG. 8  schematically illustrates a system for performing a method in accordance with an embodiment of the invention. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The present invention may be described and implemented in the general context of a system and computer methods to be executed by a computer. Such computer-executable instructions may include programs, routines, objects, components, data structures, and computer software technologies that can be used to perform particular tasks and process abstract data types. Software implementations of the present invention may be coded in different languages for application in a variety of computing platforms and environments. It will be appreciated that the scope and underlying principles of the present invention are not limited to any particular computer software technology. 
     Moreover, those skilled in the art will appreciate that the present invention may be practiced using any one or combination of hardware and software configurations, including but not limited to a system having single and/or multiple processor computers, hand-held devices, programmable consumer electronics, mini-computers, mainframe computers, and the like. The invention may also be practiced in distributed computing environments where tasks are performed by servers or other processing devices that are linked through a one or more data communications network. In a distributed computing environment, program modules may be located in both local and remote computer storage media including memory storage devices. The present invention may also be practiced as part of a down-hole sensor or measuring device or as part of a laboratory measuring device. 
     Also, an article of manufacture for use with a computer processor, such as a CD, pre-recorded disk or other equivalent devices, may include a computer program storage medium and program means recorded thereon for directing the computer processor to facilitate the implementation and practice of the present invention. Such devices and articles of manufacture also fall within the spirit and scope of the present invention. 
     Referring now to the drawings, embodiments of the present invention will be described. The invention can be implemented in numerous ways, including, for example, as a system (including a computer processing system), a method (including a computer implemented method), an apparatus, a computer readable medium, a computer program product, a graphical user interface, a web portal, or a data structure tangibly fixed in a computer readable memory. Several embodiments of the present invention are discussed below. The appended drawings illustrate only typical embodiments of the present invention and therefore are not to be considered limiting of its scope and breadth. 
     The present invention relates to characterizing the subsurface using seismic data by estimating uncertainty in depth or structure based on the uncertainty in the heterogeneous and/or anisotropic velocity model. For the purposes of this description, the term, velocity, refers specifically to propagation velocity of the acoustic wave energy, or p-waves. The present analysis, however, can be applied with equal validity to shear wave energy, or s-waves. In this document, anisotropic velocity refers both to anisotropy and heterogeneity. The uncertainty in the anisotropic velocity model arises from the limit of velocity resolution based on detectability of a small perturbation in the curvature of the arrival time versus source-receiver offset or angle, commonly called moveout. 
     The process of transforming or migrating acquired seismic data in the time domain to the depth domain uses a velocity model. Often, each volume of similar source-receiver offset traces in a seismic survey are migrated together. The volumes of different source-receiver offsets can then be re-sorted to show the continuum of source-receiver offset traces at each output location in the migrated seismic data. In evaluating the velocity model used in a migration, one factor that can be applied to verify that the resulting model is accurate is the existence of flat gathers. That is, the response due to a particular seismic reflector is indicated at the same depth across all source-receiver offsets at the same seismic trace location, indicating no moveout. It should be noted that the method described herein is not limited to offset-domain common image gathers, but may find application in subsurface angle and subsurface angle plus azimuth gathers, offset-domain plus azimuth, and other gather methods. 
     Since it is true that a perfect velocity model would produce flat gathers, it is generally assumed that flat gathers imply that the anisotropic velocity model is correct. However, in practice the production of flat gathers does not necessarily imply that the model is perfect, because it is possible that the perturbation in the moveout is not detectable due to the noise inherent in the data, the finite source-receiver offset ranges, and the width of the seismic wavelets. In other words, flat gathers are a necessary but insufficient condition for an accurate velocity model. Furthermore, for a given set of data, there are likely different velocity models that all produce apparently flat gathers but would also yield different realizations for the subsurface structure. The resulting differences in structure may be such that a hydrocarbon reservoir appears to be larger or smaller and that the target interval for a selected well location can vary in depth. Both factors can result in improper selection of well locations and drill depths impacting the ultimate productivity and economic value of the reservoir. 
     The discussion above implies that the seismic data be processed with a prestack depth migration algorithm. However, those experienced in the art would realize that the seismic data could also be processed with a poststack depth migration algorithm. Although a poststack depth migration is known not to be as accurate as a prestack depth migration, certain features of this invention could be employed using an algorithm of either type. The use of a poststack depth migration algorithm in this invention precludes the analysis of migrated offset (or azimuth) versus depth gathers; however, there are other methodologies, such as image coherence, to judge the range of acceptable solutions in order to practice this invention. 
     To fully characterize the velocity field for a given subsurface volume of interest, the velocity data may need to be known in a variety of directions. The common terminology used includes v o  to indicate the vertical velocity (the seismic velocity vertically in the Earth), v nmo  to indicate the near offset moveout velocity of the seismic energy traveling in the Earth, η to represent a difference between the horizontal velocity of the seismic energy in the Earth and v nmo , δ to represent a difference between the vertical velocity of the seismic energy in the Earth and v nmo , and ε to represent a difference between the vertical and horizontal velocities of the seismic energy in the Earth. The parameters which represent anisotropy, η, δ, and ε, are known to those skilled in the art. Further, those experienced in the art will appreciate that if the symmetry axis of the velocity in the Earth is not vertical, but rather tilted at arbitrary dip and strike angles, the velocity along the symmetry axis may be substituted for the vertical velocity in the above description. The velocity orthogonal to the symmetry axis would then be substituted for the horizontal velocity in the above description. 
       FIG. 1  is a flowchart illustrating a method  100  of the present invention that uses the detectability criterion to determine depth uncertainty from seismic data and an anisotropic velocity model. This embodiment may begin, for example, with an anisotropic velocity model  10 , which may have been obtained from well log data, well core data, seismic data or some other data, and a PSDM seismic data volume  12 , which may have been migrated by wave-equation methods, Kirchhoff type methods, or some other migration method. The seismic data might also be poststack migration data and might be 2D or 3D. 
     Using the anisotropic velocity model, this embodiment of the present invention extracts functions at multiple spatial locations at step  11 , where the locations are commonly identified in terms of x in the 2D case or (x,y) in the 3D case. The functions that may be extracted are, for example, depth functions of the v o  (vertical velocity), η int  (interval η), and δ, which is also an interval function. One skilled in the art will appreciate that these functions might also be obtained in other ways, such as from regional geologic models, and that they might also be constants, i.e. have the same value for all depths. 
     The anisotropic velocity model may also be used to calculate so-called effective parameters and properties at step  13 . These parameters and properties may include η eff , which is the effective η, v nmo , which is the interval NMO velocity, and V nmo , which is the effective NMO velocity. All three of these parameters may vary in three dimensions. 
     The parameters η eff  and V nmo  describe the moveout of a given seismic event as a composite effect of the overlying layers, and can be written in the form: 
                     t   h   2     =       t   0   2     +       X   2       V   nmo   2       -       2   ⁢     η   eff     ⁢     X   4           V   nmo   2     [         t   0   2     ⁢     V   nmo   2       +       (     1   +     2   ⁢     η   eff         )     ⁢     X   2                       Eqn   .           ⁢   1               
where t o  is the zero-offset intercept time, X is the offset, and t h  is the arrival time at that offset. This equation may be called the anisotropic moveout equation. One skilled in the art will appreciate that it is possible to determine a similar relationship for angle gathers.
 
     The effective moveout velocity V nmo , the interval moveout velocity v nmo , and the vertical interval v o  are related as: 
     
       
         
           
             
               
                 
                   
                     
                       v 
                       nmo 
                       2 
                     
                     ⁡ 
                     
                       ( 
                       
                         t 
                         2 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         
                           t 
                           2 
                         
                         ⁢ 
                         
                           
                             V 
                             nmo 
                             2 
                           
                           ⁡ 
                           
                             ( 
                             
                               t 
                               2 
                             
                             ) 
                           
                         
                       
                       - 
                       
                         
                           t 
                           1 
                         
                         ⁢ 
                         
                           
                             V 
                             nmo 
                             2 
                           
                           ⁡ 
                           
                             ( 
                             
                               t 
                               1 
                             
                             ) 
                           
                         
                       
                     
                     
                       
                         t 
                         2 
                       
                       - 
                       
                         t 
                         1 
                       
                     
                   
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   2 
                 
               
             
             
               
                 
                   
                     v 
                     nmo 
                   
                   = 
                   
                     
                       v 
                       0 
                     
                     ⁢ 
                     
                       
                         1 
                         + 
                         
                           2 
                           ⁢ 
                           δ 
                         
                       
                     
                   
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   3 
                 
               
             
           
         
       
     
     Referring again to  FIG. 1 , at step  14  the seismic data is used to estimate key descriptors. These key descriptors are characteristics of the seismic data that can be easily estimated or determined by one skilled in the art. They may include, for example, the minimum and maximum offset, the angle mute, the highest usable frequency in the data at shallow depths and at the maximum depth, the maximum depth, and the general level of signal with respect to noise. Once these key descriptors are determined, the seismic data is not needed for the remainder of the method  100 . 
     The parameters and properties from steps  11 ,  13 , and  14  are passed to step  15  which computes depth uncertainty functions Δz(z) for each spatial location. This computation uses the concept of detectability, which is a measure of whether non-flat moveout is detectable. The detectability computation produces bounding velocity and eta curves, (V nmoLO , η effLO ) and (V nmoHI , η effHI ). Examples of moveout curves generated at a single location (x,y) and arrival time t 0  from these bounding velocity and eta curves can be seen in  FIGS. 3 and 4 , where the curve  30  in  FIG. 3  shows the exact moveout as completely flat, curve  32  shows the moveout at the low (LO) end of detectability, and curve  34  shows the moveout at the high (HI) end of detectability.  FIG. 4  shows the same curves  40 ,  42 , and  44  but has also thickened curve  44  to represent the half-wavelet thickness λ/2, which helps to illustrate why it is difficult to detect non-flat moveouts that fall within the bounds calculated by the detectability criteria. An example of the bounding V nmo  curves in time can be seen in  FIG. 5 , where curve  54  represents the effective V nmo  that produces a flat moveout curve such as curve  30  in  FIG. 3 , curve  50  represents the V nmoHI  that would produce a moveout curve such as curve  34  in  FIG. 3 , and curve  52  represents the V nmoLO  that would produce a moveout curve such as curve  32  in  FIG. 3 . 
     Referring again to  FIG. 1 , step  15  will use the bounding velocity and eta curves to compute the depth uncertainty functions Δz(z). To do this, the effective properties represented by the bounding curves may first be converted into interval properties. This process can be described as:
 
 v   nmo =ƒ( V   nmo ),η int   =g ( v   nmo ,η eff ), and δ= e (η int )  Eqn. 4
 
where ƒ, g, and e are functions that convert the effective properties into interval properties.
 
     For the function ƒ one option is to use a stabilized version of Eqn. 2, which is the Dix equation, applied to convert from effective moveout velocity, V nmo , to interval moveout velocity, v nmo . This example is not meant not be limiting; other velocity-conversion methods are within the scope of the present invention. 
     For the function g, one of several options is selected based on characteristics of the data set under consideration. These options may include using crossplots and curve fitting to derive a version of the relationship g that is specific to the current velocity model, using a relationship derived in a manner similar to the first option based on regional behavior observed and learned from estimation of multiple velocity models, or by stabilization of the generalized Dix-like relationship: 
     
       
         
           
             
               
                 
                   
                     
                       
                         η 
                         int 
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             t 
                             0 
                           
                           , 
                           
                             t 
                             1 
                           
                         
                         ) 
                       
                     
                     = 
                     
                       
                         
                           
                             
                               
                                 
                                   
                                     
                                       ( 
                                       
                                         
                                           ∫ 
                                           0 
                                           
                                             t 
                                             1 
                                           
                                         
                                         ⁢ 
                                         
                                           
                                             v 
                                             nmo 
                                             2 
                                           
                                           ⁢ 
                                           
                                               
                                           
                                           ⁢ 
                                           
                                             ⅆ 
                                             τ 
                                           
                                         
                                       
                                       ) 
                                     
                                     2 
                                   
                                   
                                     t 
                                     1 
                                   
                                 
                                 ⁢ 
                                 
                                   ( 
                                   
                                     
                                       
                                         η 
                                         eff 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           t 
                                           1 
                                         
                                         ) 
                                       
                                     
                                     - 
                                     
                                       
                                         η 
                                         iso 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           t 
                                           1 
                                         
                                         ) 
                                       
                                     
                                   
                                   ) 
                                 
                               
                               - 
                             
                           
                         
                         
                           
                             
                               
                                 
                                   
                                     ( 
                                     
                                       
                                         ∫ 
                                         0 
                                         
                                           t 
                                           0 
                                         
                                       
                                       ⁢ 
                                       
                                         
                                           v 
                                           nmo 
                                           2 
                                         
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         
                                           ⅆ 
                                           τ 
                                         
                                       
                                     
                                     ) 
                                   
                                   2 
                                 
                                 
                                   t 
                                   0 
                                 
                               
                               ⁢ 
                               
                                 ( 
                                 
                                   
                                     
                                       η 
                                       eff 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       
                                         t 
                                         0 
                                       
                                       ) 
                                     
                                   
                                   - 
                                   
                                     
                                       η 
                                       iso 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       
                                         t 
                                         0 
                                       
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                             
                           
                         
                       
                       
                         ( 
                         
                           
                             ∫ 
                             
                               t 
                               0 
                             
                             
                               t 
                               1 
                             
                           
                           ⁢ 
                           
                             
                               v 
                               nmo 
                               4 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               ⅆ 
                               τ 
                             
                           
                         
                         ) 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   where 
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   5 
                 
               
             
             
               
                 
                   
                     
                       n 
                       iso 
                     
                     ⁡ 
                     
                       ( 
                       
                         t 
                         1 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       1 
                       8 
                     
                     ⁡ 
                     
                       [ 
                       
                         
                           
                             
                               t 
                               1 
                             
                             ⁢ 
                             
                               
                                 ∫ 
                                 0 
                                 
                                   t 
                                   1 
                                 
                               
                               ⁢ 
                               
                                 
                                   v 
                                   nmo 
                                   4 
                                 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   ⅆ 
                                   τ 
                                 
                               
                             
                           
                           
                             
                               ( 
                               
                                 
                                   ∫ 
                                   0 
                                   
                                     t 
                                     1 
                                   
                                 
                                 ⁢ 
                                 
                                   
                                     v 
                                     nmo 
                                     2 
                                   
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   
                                     ⅆ 
                                     τ 
                                   
                                 
                               
                               ) 
                             
                             2 
                           
                         
                         - 
                         1 
                       
                       ] 
                     
                   
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   6 
                 
               
             
           
         
       
     
     One skilled in the art may recognize other methods that may be used to derive g. These examples are not meant to be limiting. 
     For the function e, an empirical relationship may be used to derive δ from η int . These types of relationships are commonly used in practice when well data is not incorporated into the anisotropic velocity estimation problem. A specific example of one of these empirical relationships is: 
                   δ   =       η   int     3             Eqn   .           ⁢   7               
but other relationships may be used based on knowledge of the particular rock formations in the subsurface.
 
     When v nmo  and δ have been computed for the high and low bounds of detectability, the vertical velocities v oLo  and v oHI  can be computed using Eqn. 3. Note that v oLo  and v oHI  are functions of time t o . The depth functions z LO (t o ) and z HI (t o ) may be computed from v oLO  and v oHI  by integrating over t o  which allows all the bounding curves to be converted from time to depth. The resultant interval depth models z LO (z) and z HI (z) are functions of depth calculated from the input velocity model, not the depth as calculated from the bounding functions. 
     Once z LO (z) and z HI (z) are computed, the depth uncertainty as a function of depth Δz(z) is computed by:
 
Δ z ( z )= z   HI ( z )− z   LO ( z ).  Eqn. 8
 
     The Δz(z) for each spatial location may be output, either to a storage device or a display device, at step  16  of method  100 . An example of a Δz(z) for a single location can be seen as curve  60  in  FIG. 6 , which shows Δz(z)/2, the +/− error estimate. The Δz(z) may also be mapped into a depth uncertainty volume, Δz(x,y,z), at step  17  which is then output at step  18 . The depth uncertainty volume gives an estimate of the depth uncertainty at any point in the volume from analysis of the seismic data alone. The depth uncertainty volume may be constructed by discrete depths or by using conformable structural layers. 
     The above analysis is described for a specific symmetry of anisotropy, the VTI medium. The analysis is equally valid for TTI (“Tilted, Transverse Isotropic”) media, when equations suitable to TTI are substituted. In TTI media, the symmetry axis becomes distinct from the vertical axis, and the equations relating interval and effective properties are generalized to accommodate this distinction. 
     The detectability computation to which step  15  of  FIG. 1  refers is explained in  FIG. 2  as method  200 . This computation generates bounding velocity and eta curves, (V nmoLO , η effLO ) and (V nmoHI , η effHI ) for all desired zero-offset intercept times t o  and each of the multiple spatial locations. In step  20 , input is received. This input may include the effective parameters from step  13  of  FIG. 1 , the velocity functions from step  11  of  FIG. 1 , and the key descriptors from step  14  of  FIG. 1 . Method  200  then begins at the first desired t o , set at step  21 . For this t o , the maximum usable frequency f max  is computed based on the usable frequency information from the seismic data key descriptors. This may be done based on a linear trend computed from a f max  determined for shallow depths and a f max  determined for the maximum depth or may be selected from a f max  function in time or depth that was supplied as a key descriptor. Step  21  also computes the X max  maximum offset within the angle mute at this t o  or may set X max  to the maximum offset from the key descriptors. If step  13  of  FIG. 1  did not compute the required effective parameters η eff ; v nmo , and V nmo  for this t o , they may also be computed at step  21  using the calculations explained for step  13  of  FIG. 1 . 
     Referring again to  FIG. 2 , the detectability threshold D is computed at step  23 . This threshold may be based on prior knowledge or may be computed based on the signal to noise ratio of the seismic data, which is one of the key descriptors, and the f max  computed at step  22 . For example, the threshold might be set as D=(τ, τ/2, τ/3, or τ/4) where τ=1/f max , depending on whether the data is very noisy, noisy, normal, or clean, respectively. 
     In step  24 , the V nmoLO , η effLO , V nmoHI , and η effHI  are calculated for this t o . This may be done, for example, by the following process:
         calculate a base arrival time t base  for each offset between the minimum offset supplied by the key descriptors and the X max  from step  22 , and calculated from the velocity model;   calculate test arrival times t test  for each offset by scanning over a range of effective NMO velocities V nmo  (such as a range from 0.7*V nmo −1.2*V nmo ) and a range of effective η eff  (such as a range from 0.0 to 0.3) and using Eqn. 1 to compute t test ;   compute the difference between each t test  and t base  and keep the largest negative and positive differences, Δt −   test  and Δt +   test  respectively;   if the total difference between Δt −   test  and Δt +   test  is less than the detectability threshold D, the V nmo , η eff  pair associated with that t test  is within the uncertainty, as its residual moveout is not detectable;   when all V nmo , η eff  pairs within the uncertainty are identified through the scanning process, find the minimum and maximum V nmo  values and call those V nmoLO  and V nmoHI . The η eff  values that are paired with those V nmo  values will be η effLO  and η effHI .       

     At step  25 , the process may move on to the next t o  and repeat steps  22 - 24 . Once the maximum time desired is reached, the process may move on to step  26  where the bounding velocity and eta curves, (V nmoLO , η effLO ) and (V nmoHI , η effHI ) for all t o  are produced to be passed on for further computation in step  15  of  FIG. 1 . 
     The method  200  of  FIG. 2  described above is one aspect of the present invention and is not meant to be limiting. Other methods of setting a detectability threshold and calculating the bounding velocities and etas fall within the scope of this method. 
     The depth uncertainty volume Δz(x,y,z) produced by method  100  of  FIG. 1  may be used directly by seismic analysts and interpreters and may also be used to determine structural uncertainty, as explained by the method  700  in  FIG. 7 . This method uses the depth uncertainty volume and one or more interpreted seismic horizons to produce multiple structural realizations that lay within the depth uncertainty bounds. These multiple realizations can aid in analysis of geometrical and structural properties of the subsurface. 
     The method  700  receives Δz(x,y,z) at step  70  and the interpreted horizons from the seismic data at step  71 . At step  72 , for each horizon, i, interpreted from the imaged seismic data volume, Z hzni (x,y), on which uncertainty analysis is to be performed, extract the depth uncertainty from Δz(x,y,z) at each point on the interpreted surface, giving, ΔZ hzni (x,y). This quantity may be called the “cage” of uncertainty for this horizon i, since it describes the allowable limits that the horizon may be distorted in depth while safely remaining within the limits of uncertainty. 
     At step  74 , the horizon is distorted in depth in multiple ways to create unique realizations, or realistic versions, of the horizon, each of which differs in its depth structure from the horizon as originally interpreted, but remains confined within the constraints of the depth uncertainty resulting from the above analysis. One skilled in the art will appreciate that the distortion of the horizon may be done in many ways. One method is to start with the uncertainty cage for interpreted horizon Z hzn (x,y), ΔZ hzn (x,y), where we drop the subscript, i, and understand that we refer to the i th  interpreted horizon, and where a corresponding analysis can for done for another interpreted horizon. It is possible to create a Deviation Surface appropriate for this horizon, D(x,y), such that for each location, (x,y), the Deviation Surface is constrained by |D(x,y)|≦ΔZ hzn (x,y)/2, and where D(x,y) follows the rules listed below. The Deviation Surface is used to create an alternate realization of the horizon surface according to the relation, R hzn (x,y)=Z hzn (x,y)+D(x,y), where R is the alternate realization. The Deviation Surface may be created mathematically, such as to follow a tilted plane, or it may be created using a geologically-driven model, such as by mimicking the structural trend of an overlying stratum. Multiple alternate realizations are obtained by generating alternate versions of D(x,y). Using the subscript, j, we denote the j th  realization as R hznj (x,y)=Z hzn (x,y)+D j (x,y). The manner in which these horizons are distorted to create unique realizations is constrained to follow specific rules. An exemplary set of three rules is listed below. Other rules can be generated to constrain further the allowable set of realizations.
         Each realization must portray a depth structure for the horizon which is geologically reasonable. Whether a depth structure is geologically reasonable depends on the depositional environment and subsequent geologic history of the subsurface of which one skilled in the art will be aware. This rule may act as a rejection criterion for a realization from a mathematically-created Deviation Surface. It can also be a construction constraint for a geologically-created Deviation Surface, for example that it mimics the isopach contours of an overlying stratigraphic layer.   Spatial variations in structural features added to create each new realization, e.g., undulations, must not exceed the original interpreted horizon in complexity. This rule can also be stated that the distribution of spectral amplitude over the range of spatial wavenumbers of the horizon spectrum, A(k x ,k y ), should be similar in each realization to the spectrum of the original horizon.   Surface discontinuities cannot be created where not present in the original interpretation, and must be maintained where present in the original interpretation.       

     Once the set of horizon realizations are created, they can be output at step  73  to a data storage device, a display device, and/or further processing. For example, they may be passed to step  75  for continuing analysis. 
     Continuing analysis focuses on the variation in some specific geometrical or structural property of the horizon, which may differ according to each realization. For example, when performing hydrocarbon reserve calculations, the container volume within the spill point contour, the column height, and the area of the spill point contour can be calculated for each realization. 
     Other properties can be calculated for each realization. Each property calculated from the multiple realizations of structural surfaces can be portrayed graphically as a distribution at step  77 , from which P10 and P90 estimates of each property can be made at step  76 . These P10 and P90 estimates may also be output to a data storage and/or display device. 
     A system  800  for performing the method is schematically illustrated in  FIG. 8 . The system includes a data source/storage device  82  which may include, among others, a data storage device or computer memory. The device  82  may contain an anisotropic velocity model, seismic data, parameters relating to the velocity model, parameters relating to the seismic data, and/or a depth uncertainty volume. The data from device  82  may be made available to a processor  84 , such as a programmable general purpose computer. The processor  84  is configured to execute a processing module  86  that can perform the depth uncertainty method can interact with the detectability module  85 . The processor  84  may also execute the structural uncertainty module  87 . The system may include interface components such as user interface  88 , and is used to implement the above-described transforms in accordance with embodiments of the invention. The user interface  88  may be used both to display data and processed data products and to allow the user to select among options for implementing aspects of the method. By way of example and not limitation, the depth uncertainty volume and/or multiple structural realizations computed on the processor  84  may be displayed on the user interface  88 , stored on the data storage device or memory  82 , or both displayed and stored. 
     While in the foregoing specification this invention has been described in relation to certain preferred embodiments thereof, and many details have been set forth for purpose of illustration, it will be apparent to those skilled in the art that the invention is susceptible to alteration and that certain other details described herein can vary considerably without departing from the basic principles of the invention. In addition, it should be appreciated that structural features or method steps shown or described in any one embodiment herein can be used in other embodiments as well.