1000850305 SYSTEM AND METHOD FOR SUBSURFACE CHARACTERIZATION INCLUDING UNCERTAINTY ESTIMATION FIELD OF THE INVENTION [0001] 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. [0002] This application is related to Australian Patent Application No. 2012295567, filed on 18 April 2012, which is incorporated herein in its entirety by reference. BACKGROUND OF THE INVENTION [0003] 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. [0004] 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 1 1000850305 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. [0005] Reference to any prior art in the specification is not an acknowledgment or suggestion that this prior art forms part of the common general knowledge in any jurisdiction or that this prior art could reasonably be expected to be understood, regarded as relevant, and/or combined with other pieces of prior art by a skilled person in the art. SUMMARY OF THE INVENTION [0006] Described herein are implementations of various approaches for a computer implemented method for subsurface characterization from seismic data. The described method may include receiving, at a computer processor, effective parameters and properties from an anisotropic seismic velocity model, at least one velocity function, and key descriptors from seismic data; determining, via a computer processor, a detectability threshold for moveout in a seismic data gather based on the seismic data; and computing, via a computer processor, 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. [0007] In an embodiment, the computer-implemented method may include receiving, at a computer processor, 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. 2 1001352833 [0008] In a further embodiment, a system for subsurface characterization from seismic data, the system comprising: a data source containing data representative of the seismic data, a depth uncertainty volume and a depth uncertainty cage extracted from the depth uncertainty volume; at least one computer processor being configured to communicate with the data source and to execute a structural uncertainty module for producing multiple realizations of at least one interpreted horizon in the seismic data based on the depth uncertainty cage and analyzing the multiple realizations for changes to geometrical or structural properties of the at least one interpreted horizon; and a user interface. [0008A] In a further embodiment, an article of manufacture including a computer readable medium having computer readable code on it, the computer readable code being configured to implement a method for subsurface characterization from seismic data, the method comprising: a. Extracting a depth uncertainty cage for at least one interpreted horizon from the seismic data based on a depth uncertainty volume; and b. Simulating multiple realizations for each of the at least one horizons, constrained by the depth uncertainty cage; and c. Analysing the multiple realizations for changes to geometrical or structural properties of at least one interpreted horizon. [0009] Additional embodiments include systems and articles of manufacture configured to implement the methods described herein. [0010] 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 [0011] These and other features of the present invention will become better understood with regard to the following description, claims and accompanying drawings where: [0012] Figure 1 is a flowchart illustrating a method in accordance with an embodiment of the present invention; 3 1001352833 [00131 Figure 2 is a flowchart illustrating a method for using detectability in accordance with an embodiment of the present invention; [0014] Figure 3 illustrates exemplary moveout curves computed by a method in accordance with the present invention; [0015] Figure 4 illustrates exemplary moveout curves including wavelet thickness computed by a method in accordance with the present invention; 3a 1000850305 [0016] Figure 5 illustrates exemplary velocity curves computed by a method in accordance with the present invention; [0017] Figure 6 illustrates an exemplary depth uncertainty curve computed by a method in accordance with the present invention; [0018] Figure 7 is a flowchart illustrating a method in accordance with an embodiment of the present invention; and [0019] Figure 8 schematically illustrates a system for performing a method in accordance with an embodiment of the invention. DETAILED DESCRIPTION OF THE INVENTION [0020] 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. [0021] 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. 4 1000850305 [0022] 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. [0023] 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. [0024] 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. [0025] 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 5 1000850305 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. [0026] 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. [0027] 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. [0028] 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) to indicate the vertical velocity (the seismic velocity vertically in the Earth), vumo to indicate the near offset moveout velocity of the seismic energy traveling in the Earth, r to represent a difference between the horizontal velocity of the seismic energy in the Earth and 6 1000850305 vumo, 6 to represent a difference between the vertical velocity of the seismic energy in the Earth and vumo, and s to represent a difference between the vertical and horizontal velocities of the seismic energy in the Earth. The parameters which represent anisotropy, r, 6, and s, 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. [0029] Figure 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. [0030] 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 vo (vertical velocity), lint (interval f), and 6, 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. [0031] 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 l, vumo, which is the interval NMO velocity, and Vmo, which is the effective NMO velocity. All three of these parameters may vary in three dimensions. [0032] The parameters eff and Vmo describe the moveout of a given seismic event as a composite effect of the overlying layers, and can be written in the form: 7 1000850305 2 t2 + _ x2___ _____ ____ t= + - vfletfX+ Eqn. 1 V.2." [t2 V2." + (12ff)X2] where to is the zero-offset intercept time, X is the offset, and th 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. [0033] The effective moveout velocity Vmo, the interval moveout velocity vamo, and the vertical interval vo are related as: V'mo(t2 ) = t 2 viot 2 )- tiV12Iy(ti) Eqn. 2 t 2 - tI Vnmo = voV1 + 26 Eqn. 3 [0034] Referring again to Figure 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. [0035] The parameters and properties from steps 11, 13, and 14 are passed to step 15 which computes depth uncertainty functions Az(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, (VnmoLO, leffLO) and (VnmoHI, lefifl). Examples of moveout curves generated at a single location (x,y) and arrival time to from these bounding velocity and eta curves can be seen in Figures 3 and 4, where the curve 30 in Figure 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. Figure 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 8 1000850305 the bounding Vumo curves in time can be seen in Figure 5, where curve 54 represents the effective Vumo that produces a flat moveout curve such as curve 30 in Figure 3, curve 50 represents the VnmoHI that would produce a moveout curve such as curve 34 in Figure 3, and curve 52 represents the VnmoLO that would produce a moveout curve such as curve 32 in Figure 3. [0036] Referring again to Figure 1, step 15 will use the bounding velocity and eta curves to compute the depth uncertainty functions Az(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: Vnmo = f(Vnmo), 77int = g(vnmo,eff), and 6 = e(l7int) Eqn. 4 wheref, g, and e are functions that convert the effective properties into interval properties. [0037] For the function f one option is to use a stabilized version of Eqn. 2, which is the Dix equation, applied to convert from effective moveout velocity, Vumo, to interval moveout velocity, vamo. This example is not meant not be limiting; other velocity-conversion methods are within the scope of the present invention. [0038] 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: (fol vkimo dr) 2 (f o voimocdr) 2 1ftt tl (feff(ti)- iso (ti))- t 0 (eff(to) - iso (to)) 77int (to, t1) =iV4Odrt E qn. 5 where 1 rt 1 1 Amod 1 SsftimoT 2-1 Eqn. 6 9 1000850305 [0039] One skilled in the art may recognize other methods that may be used to derive g. These examples are not meant to be limiting. [0040] For the function e, an empirical relationship may be used to derive 6 from . 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: 6 = Eqn. 7 but other relationships may be used based on knowledge of the particular rock formations in the subsurface. [0041] When vamo and 6 have been computed for the high and low bounds of detectability, the vertical velocities voLo and voHI can be computed using Eqn. 3. Note that voLO and voHI are functions of time to. The depth functions ZLO(to) and ZHI(to) may be computed from v o Lo and voHI by integrating over to which allows all the bounding curves to be converted from time to depth. The resultant interval depth models ZLO(Z) and ZHI(Z) are functions of depth calculated from the input velocity model, not the depth as calculated from the bounding functions. [0042] Once ZLO(Z) and ZHI(Z) are computed, the depth uncertainty as a function of depth Az(z) is computed by: Az(z) = ZHI(Z) - ZLO (Z). Eqn. 8 [0043] The Az(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 Az(z) for a single location can be seen as curve 60 in Figure 6, which shows Az(z)/2, the +/- error estimate. The Az(z) may also be mapped into a depth uncertainty volume, Az(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. 10 1000850305 [0044] 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. [0045] The detectability computation to which step 15 of Figure 1 refers is explained in Figure 2 as method 200. This computation generates bounding velocity and eta curves, (VnmoLO, leffLO) and (VnmoHI, feffHI) for all desired zero-offset intercept times to and each of the multiple spatial locations. In step 20, input is received. This input may include the effective parameters from step 13 of Figure 1, the velocity functions from step 11 of Figure 1, and the key descriptors from step 14 of Figure 1. Method 200 then begins at the first desired to, set at step 21. For this to, the maximum usable frequency fmax 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 fmax determined for shallow depths and a fmax determined for the maximum depth or may be selected from a fmax function in time or depth that was supplied as a key descriptor. Step 21 also computes the Xmax maximum offset within the angle mute at this to or may set Xmax to the maximum offset from the key descriptors. If step 13 of Figure 1 did not compute the required effective parameters leff, vamo, and Vumo for this to, they may also be computed at step 21 using the calculations explained for step 13 of Figure 1. [0046] Referring again to Figure 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 fmax computed at step 22. For example, the threshold might be set as D=(T, T/2, T/3, or T/4) where T=1/fmax, depending on whether the data is very noisy, noisy, normal, or clean, respectively. [0047] In step 24, the VnmoLO, leffLO, VnmoHI, and leffHJ are calculated for this to. This may be done, for example, by the following process: - calculate a base arrival time tbase for each offset between the minimum offset supplied by the key descriptors and the Xmax from step 22, and calculated from the velocity model; 11 1000850305 - calculate test arrival times test for each offset by scanning over a range of effective NMO velocities Vamo (such as a range from 0.
7 *Vmo -1.2*Vamo) and a range of effective eff (such as a range from 0.0 to 0.3) and using Eqn. 1 to compute ttest; - compute the difference between each ttest and tbase and keep the largest negative and positive differences, At-test and At'test respectively; - if the total difference between At-test and At'test is less than the detectability threshold D, the Vumo, eff pair associated with that ttest is within the uncertainty, as its residual moveout is not detectable; - when all Vumo, eff pairs within the uncertainty are identified through the scanning process, find the minimum and maximum Vmo values and call those VnmoLO and VnmoHI. The eff values that are paired with those Vumo values will be leffLO and leffHI [0048] At step 25, the process may move on to the next to 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, (VnmoLO, leffLO) and (VnmoHI, leffHI) for all to are produced to be passed on for further computation in step 15 of Figure 1. [0049] The method 200 of Figure 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. [0050] The depth uncertainty volume Az(x,y,z) produced by method 100 of Figure 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 Figure 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. [0051] The method 700 receives Az(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, Zhzn 1 i(x,y), on which uncertainty analysis is to be performed, extract the depth 12 1000850305 uncertainty from Az(x,y,z) at each point on the interpreted surface, giving, AZhzni(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. [0052] 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 Zhzn(x,y), AZhzn(x,y), where we drop the subscript, i, and understand that we refer to the ih 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 ID(x,y) < AZhzn(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, Rhzn(X,y) = Zhzn(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 jth realization as Rhznj(X,y) = Zhzn(x,y) + Dj(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. 13 1000850305 - 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(kx,ky), 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. [0053] 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. [0054] 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. [0055] 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. [0056] A system 800 for performing the method is schematically illustrated in Figure 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, 14 1000850305 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. [0057] 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. [0058] As used herein, except where the context requires otherwise, the term "comprise" and variations of the term, such as "comprising", "comprises" and "comprised", are not intended to exclude further features, components, integers or steps. 15