Interpretation of geologic depositional systems

196039/2

INTERPRETATION OF GEOLOGIC DEPOSITIONAL SYSTEMS 08901/0040*000 196039/2

RELATED APPLICATION DATA

[0001] This application claims the benefit of and priority under 35 U.S. C. § 119(e) toU.S. Patent Application Nos. 60/815,630, filed 21 June 2006, entitled "Algorithm and Process toCreate Geobody Bounding Surfaces," 60/815,625, filed 21 June 2006, entitled "ComputedAided and Automatic Extraction of Depositional Systems," and 60/815,961, filed 21 June 2006,entitled "Stratal-Slice Domain Transformation of a Seismic Volume".

BACKGROUND

[0002] An exemplary embodiment of this invention is in the field of 3-D interpretation, and more particularly to 3-D seismic interpretation. More specifically, an exemplaryembodiment includes a workflow, including two new processes, implemented as software thatis designed to enable automatic or semi-automatic interpretation of paleo-depositional featuresin three-dimensional seismic data for exploration, development and, for example, productionof hydrocarbons.

[0003] The need for computer-aided, semi-automatic and automatic interpretation ofdepositional systems derives from a combination of factors. Energy resources are becomingsteadily more difficult to find and develop. It has been recognized for many years that themajority of new oil and gas reserves are a function of a complex combination of geological,structural and stratigraphic elements. While the problems of exploration and the efficientdevelopment of hydrocarbon reserves have become more difficult, the volume of data to beinterpreted for each project has become orders of magnitude greater over the past 20 years.Simultaneously, both the number of interpreters and the time allowed for interpretation havebeen substantially reduced. This drives the need for more advanced computer-aided processes 1 WO 2008/005690 PCT/US2007/071733 that can support the interpreter by enabling more efficient, precise and effective interpretation of 3-D seismic data volumes.

[0004] Computer-aided structural interpretation of 3-D seismic data volumes has beenembodied in tools in interactive seismic interpretation for a number of years. Since the early1980s, horizon autotracking tools have been available to help increase the speed andconsistency of horizon interpretation in 3-D seismic surveys (Dorn, 1998). More recently,techniques have been developed to provide computer-aided interpretation of faults andautomatic fault interpretation (e.g., Crawford and Medwedeff, 1999, US Patent 5,987,388;Pederson, S. I., 2002, US Patent 7,203,342), as well as techniques beyond event autotrackingto automatically interpret horizons (Dorn, 1999, US Patent 5,894,417; Stark, 1997, US Patent5,671,344).

[0005] Computer tools to aid in stratigraphic interpretation of seismic volumes havedeveloped much more slowly. Elements of depositional systems can most readily beidentified by an interpreter when the morphology of the paleo-depositional system can beviewed. Similarly, it is most likely that a computer algorithm can be written to recognize,image, and extract elements of depositional systems if the computer algorithm is able tooperate on the data in a domain where the paleo-depositional system's morphology is mostreadily imaged. In both of these cases, the optimal environment is the stratal-slice domain,where the slices through the volume of seismic data are close approximations of paleo-depositional surfaces.

[0006] In an undeformed data volume, horizontal slices (planar slices parallel to the (x,y)plane in the volume) may accurately represent depositional surfaces. However, in volumeswith structural deformation, horizontal slices do not represent depositional surfaces for morethan a small portion of the total volume. Faulting, folding, and velocity anomalies prevent thecomplete representation of such a surface by a simple horizontal slice. 2 WO 2008/005690 PCT/US2007/071733 [0007] Horizon-slicing is defined as creating a slice through a 3-D seismic volume in theshape of an interpreted seismic reflection in that volume. Horizon slicing (as opposed tohorizontal slicing) has provided better images of depositional systems since the mid-to-late 1980s.

[0008] A continuous interval is a package of sediments that represent the same span ofgeologic age, but were deposited at different sedimentation rates in different parts of thevolume. The result is an interval that represents that same amount of geologic time, but doesnot exhibit the same thickness. In such an interval, growth is caused by spatially variablerates of sedimentation. If we assume that sedimentation rates between a pair of boundinghorizons are variable only in space (i.e., not vertically variable in a given location), stratalslices may be extracted by interpolating trace values vertically, where the interpolated sampleinterval at each (x,y) location is controlled by the upper and lower bounding surfaces and thenumber of samples desired in the interval on the output trace. This type of stratal slice hasbeen referred to as a proportional slice.

[0009] Proportional slicing or stratal slicing developed in the mid 1990s (Posamentier, et.al., 1996; Zeng, et. al. 1998a,b,c) provides even better imaging of depositional systems, andbetter discrimination between stacked channel systems in the seismic data because thesesurfaces are typically a better approximation of paleo-depositional surfaces than either horizon slices or horizontal slices.

[0010] Zeng, et. al. (1998 a, b, c) describes the first instance of extracting slices based on geologic age. They reasoned that seismic reflectors do not always follow depositional surfaces. Thus, they interpolated seismic slices between surfaces judged to be time-equivalent.

They referred to these interpolated slices as 'stratal' slices. Stark (2004) describes a similarly motivated effort. He used unwrapped phase as a proxy for user-interpreted age horizons.

Slices were extracted from the data volume by drawing data from points of equal unwrapped3 WO 2008/005690 PCT/US2007/071733 phase. Stark's approach assumes that unwrapped phase closely approximates geologic age,but this is an assumption that is often in error.

[0011] Both horizon slicing and proportional slicing generally suffer from substantiallimitations in that they do not accommodate and compensate for generalized 3-D structuraldeformation subsequent to deposition, nor do they properly account for the wide variety ofdepositional environments. Horizon slicing is only appropriate for a conformable sequence ofhorizons in the seismic volume (i.e., a spatially uniform depositional environment over time).Proportional slicing is only appropriate for an interval that exhibits growth (i.e., a spatiallygradational change in depositional thickness over an area, often due to spatially differentialsubsidence). Horizon and proportional slicing do not properly reconstruct paleo-depositionalsurfaces in other depositional environments, nor do they account for post-depositionalstructural changes (particularly faulting) or post-depositional erosion.

[0012] Among the situations that the proportional or stratal slice volume (as defined byZeng, et. al, 1998 a,b,c) does not handle properly are: • Angular unconformities • Non-linear growth in the interval between two horizons • Carbonate platform intervals • Faulting [0013] For example, both proportional slicing and stratal slicing (as defined by Zeng, et.al., 1998 a, b, c) produce volumes that have gaps or undefined zones where the volume is cutby a dipping fault surface. Figure 1 shows, in a 2-D cross-section, the effects of dipping faultson this simple type of proportional slice for a pair of horizons. The output proportional slicevolume is null or indeterminate at all (x,y) positions where one or both horizons is missing(e.g., Null Zone - 1 Horizon in Figure 1). The proportional slice volume is also indeterminatefor (x,y) positions where both horizons are present but on opposite sides of a dipping fault 4 WO 2008/005690 PCT/US2007/071733 surface (e.g., Null Zone - 2 Horizons in Figure 1).

[0014] The situation for more than one pair of horizons is shown in Figure 2. In this case,there are null or indeterminate zones for each pair of horizons. These indeterminate zones arein different (x,y) positions for each pair of horizons.

[0015] In most previous attempts to solve this problem, where this simple form ofproportional slicing is implemented, the indeterminate zones are filled with input seismic datarather than nulls, which can be quite misleading. Lomask et. al. (2006) have developed anapproach that attempts to create a stratal volume without requiring interpreted horizons, faultsor other surfaces to define and constrain the transformation. The lack of interpreted structuralcontrol in their approach produces poor results for seismic volumes that contain anysignificant structural deformation.

[0016] One exemplary embodiment of the Domain Transform method of this inventionexplicitly requires interpreted horizons, faults, and other geologic surfaces as input, and, as aresult, does not suffer the limitations of the method proposed by Lomask.

[0017] Seismic-Wheeler Volumes (e.g., Stark, 2006) represent interpreted depositional systems tracts as well as hiatuses in deposition based on horizon interpretations of systemboundaries in 3-D. This approach requires recognition of the system tract by the interpreter asa starting point, and does not take into account the effects of post-depositional structuraldeformation and faulting. The implementations of Seismic-Wheeler Volumes described byStark (2006) also depend on association of each seismic sample in the volume with a relativegeologic time (Stark, 2004; Stark 2005a, US Patent 6,850,845; Stark 2005b, US Patent6,853,922). This constraint is not present in the process described here.

[0018] By transforming seismic data from the (x,y,time/depth) domain to the (x,y,stratal- slice) domain, data in a deformed volume can be interpreted in stratal-slice view. One exemplary goal is to reconstruct the data volume along stratal surfaces in an undeformed state5 WO 2008/005690 PCT/US2007/071733 using user-interpreted surfaces and user-supplied information on geologic relationships in thevolume as a guide. Seismic data in this undeformed state is more easily and accuratelyinterpreted for stratigraphy, depositional systems, and depositional environments.

[0019] Finally, a lightweight representation of volumetric data is often necessary for real-time rendering, for the segmentation of interpreted data, and for reducing visual clutter. Anew Surface Wrapping technique has also been developed in accordance with an exemplaryembodiment of this invention, and is described herein. For example, it al lows, for example,the user to create a 3-D polygonal mesh that conforms to the exterior boundary of geobodies(such as stream channels) that offers significant improvements over existing techniques.[0020] An inspiration for this Surface Wrapping approach was the Surface Drapingalgorithm (Dorn, 1999, US Patent 5,894,417), which allows a polygonal mesh to be definedthat reflects the geometry of an interpreted horizon. The surface draping algorithm is basedon the metaphor of laying an elastic sheet over a contoured surface: gravity pulls the sheetdown, causing it to conform to the surface beneath it, and the tension of the elastic materialallows the sheet to smoothly cover small gaps in the surface while preserving the important features.

[0021] Dorn’s Surface Draping allows the user to view seismic data and define a series ofpoints slightly above the desired horizon. These points define the initial shape of the 3-Dmesh, which corresponds to the elastic sheet described above. When the user has completedthis stage, the actual mesh is computed, generally using one vertex per voxel. These verticesare then iteratively "dropped" onto the horizon. At each step, the value of the voxel at eachvertex's position is compared to a range that corresponds to the values found in an interpretedhorizon. If the value falls within that range, the vertex is fixed in place.

[0022] The Surface Draping concept would have benefits if adapted to work ongeobodies and other 3-D volumes. Other approaches have been used to define a mesh that 6 WO 2008/005690 PCT/US2007/071733 surrounds and conforms to the shape of a volume. Acosta et. al. (2006a and b; US Patents7,006,085 and 7,098,908) propose a technique where the bounding surface is defined slice-by-slice by a user as a set of spline curves or general polylines that are then connected in 3-D.Kobbelt et. al. (1999) describes a technique based on successive subdivision of an initiallysimple mesh that completely surrounds the volume. The technique described by Koo et. al.(2005) improves on the same idea by allowing the user to define an arbitrarily shaped gridaround a point cloud, allowing holes in the volume to be interpreted properly. Both of theabove algorithms work by moving each vertex to the nearest point in the volume.

SUMMARY

[0023] It is an aspect of the present invention to provide a workflow and automated or semi-automated method and system for identifying and interpreting depositionalenvironments, depositional systems and elements of depositional systems from 3-D seismic volumes.

[0024] It is a further aspect of this invention to provide such a method and system inwhich noise in the seismic volume after acquisition and seismic processing is removed orminimized at each step in the workflow.

[0025] It is a further aspect of this invention to provide a technique whereby the original 3-D seismic volume is transformed to a volume where every horizontal slice through thevolume represents a paleo-depositional (stratal) surface, such that the effects of structuraldeformation are effectively removed from the volume.

[0026] It is a further aspect of this invention to provide a means of imaging, recognizingand obtaining the bounding surfaces of depositional systems or elements of depositionalsystems in the transformed seismic volume.

[0027] It is a further aspect of this invention to provide a technique whereby the imaging 7 WO 2008/005690 PCT/US2007/071733 or attribute volumes created from the transformed seismic volume may be inverse-transformed to the coordinate space of the original seismic volume.

[0028] It is a further aspect of this invention to provide a technique whereby the boundingsurfaces obtained for depositional systems or elements of depositional systems obtained inthe transformed seismic volume may be inverse-transformed to the coordinate space of theoriginal seismic volume.

[0029] In accordance with an exemplary embodiment of this invention, an approachincluding a unique new workflow that includes a combination of existing and new novelprocesses is presented for computer-aided interpretation of depositional systems in 3-Dseismic volumes. In this discussion, channels are used as the example of a depositionalsystem, but the approach will work for the full range of depositional systems and environments recorded in 3-D seismic data volumes.

[0030] This unique workflow includes the following general steps, illustrated in Figure 3a: • Load (Input) 3-D Seismic Volume • Structural Interpretation • Domain Transformation • Optional Structural Refinement • Stratigraphic Interpretation • Inverse Domain Transformation • Output Stratigraphic Volumes and Bodies [0031] Individual steps and series of steps of this workflow may be applied recursively tothe data volume to improve the results of the overall process.

[0032] Structural Interpretation refers to the interpretation of horizons and faults imagedin the 3-D seismic volume. The original seismic volume and its structural interpretation is 8 WO 2008/005690 PCT/US2007/071733 typically described in an orthogonal cartesian coordinate system indicated by (x,y,z) or (x,y,t),where x and y represent horizontal distance, z represents vertical distance, and t typicallyrepresents vertical composite (also called two-way) reflection time. The proposed workflowcan be applied to volumes that have been processed into either (x,y,z) or (x,y,t) volumes.[0033] Domain Transformation refers to a novel process of changing the coordinate spaceof the seismic volume from (x,y,z or t) to (x,y,s) where s represents "stratal-slice." A stratal-slice is defined as a slice along an approximate paleo-depositional surface, that is, a surfaceupon which at some time in the past, geologic deposition (e.g., sedimentation or erosion) wasoccurring. The Domain Transformation creates a stratal-slice volume - a volume where everyhorizontal slice in the volume represents a stratal-slice or paleo depositional surface. Thisstratal-slice volume, created by the Domain Transformation process, is a volume that issubstantially free of deformation. This Domain Transformation process is unique in that itremoves the effects of deformation that has occurred both during and subsequent to thedeposition, and will properly construct a stratal-slice volume for all types of geologic surfaces and intervals.

[0034] The Domain Transformation not only produces an ideal volume for theinterpretation or extraction of elements of depositional systems, it also provides a valuabletool to highlight errors or omissions in the structural interpretation. Such errors or omissionsare highlighted in the domain-transformed volume. Using the transformed volume to imageproblems in the structural interpretation, Optional Structural Refinement uniquely enables theinterpreter to correct these errors and omissions in either the (x,y,s) volume or the (x,y,z or t)volume and improve both the structural interpretation and the results of the Domain

Transformation.

[0035] Stratigraphic Interpretation, as used here, encompasses both the processing of theDomain Transformed volume to improve the imaging of elements of depositional systems 9 WO 2008/005690 PCT/US2007/071733 (herein referred to as attribute calculation), and the process of extracting the boundingsurfaces of those elements of depositional systems. The bounding surface extraction process(herein referred to as Surface Wrapping) is a unique process that provides numerousadvantages over processes currently practiced by individuals with ordinary skill in the art toobtain the bounding surfaces of elements of depositional systems. Surface Wrapping'sapplicability extends to the extraction of the bounding surfaces of bodies or aspects imaged inany type of volumetric data from any discipline.

[0036] Inverse Domain Transformation refers to a process of changing the coordinatespace of the seismic volume, any attribute volumes, the refined structural interpretation, andbounding surfaces from (x,y,s) to (x,y,z or t).

[0037] As mentioned, individual steps and series of steps may be applied recursively tothe data volume to improve the results of the overall process. For example, initial StructuralInterpretation of key horizons and major faults followed by Stratal-slice DomainTransformation of the Seismic Volume may highlight secondary horizons or smaller faults(additional structural interpretation) that must be interpreted and honored in the DomainTransformation process to achieve higher quality results.

[0038] Numerous exemplary benefits derive from the workflow and processes contained therein. • Domain Transformation creates a stratal-sliced volume for any seismicvolume of any structural complexity. Structural effects can be removed from the volume. • An optimized view of stratigraphic features provides improved recognitionand interpretation of depositional features. • Stratigraphic features that are obscured by structural deformation in the input seismic volume are clearly imaged, recognizable, and interpretable in the10 WO 2008/005690 PCT/US2007/071733 transformed volume. • A unique check and refinement of the structural interpretation is provided bythe transformed volume (stratal-sliced volume), refinement of the structuralinterpretation (e.g., horizons and faults) in the stratal domain, and inverse Domain

Transformation of the refined structural surfaces. • This unique workflow integrates the structural and stratigraphic interpretationsof the data in an internally self-consistent manner, not possible previously, therebyimproving the quality of the interpretation. • A novel Surface Wrapping technique provides a tool to obtain a connected,closed bounding surface for a 3-D object (geobody, stratigraphic feature, or anyother arbitrary 3-D body) even where portions of the body are poorly imaged in the data volume. • This Surface Wrapping technique has wide application to the extraction ofcomplex 3-D bodies from any form of volumetric data. Applications would alsoinclude, but are not limited to, obtaining the bounding surfaces of complex 3-Dsalt bodies and canyons in seismic data, and obtaining the bounding surfaces oftissue structures imaged in 3-D medical imaging volumes (e.g., CT, MRI, MRA,PET volumes, and the like). • Seismic attributes that are determined using the Domain Transformed volumeshow a significant improvement both qualitatively and quantitatively whencompared to the same attributes calculated using the original input seismic volume.The workflow of first transforming the volume, determining the attribute, and theninverse transforming the attribute volume produces markedly improved resultswhen compared to directly determining the attribute on the input seismic volume. • Attributes determined using the Domain Transformed volume improve the 11 WO 2008/005690 PCT/US2007/071733 quality of the imaging of stratigraphy when compared to the same attributesdetermined using the input seismic volume:

Improves the imaging of depositional systems

Improves the correlation of attributes with well data for improved geophysical reservoir characterization. • Since stratigraphic features are better imaged, more complete and more easilyinterpreted in the transformed domain, the workflow improves the efficiency,accuracy, and completeness of the interpretation of depositional systems whencompared to other approaches.

[0039] This Summary of the Invention is neither intended nor should it be construed asbeing representative of the full extent and scope of the present invention. While variousembodiments of the present invention have been described in detail, it is apparent that modifications and alterations of those embodiments will occur to those skilled in the art.

However, it is to be expressly understood that such modifications and alterations are withinthe scope and spirit of the present invention.

[0040] These and other features and advantages of this invention are described in, or areapparent from, the following detailed description of the exemplary embodiments.

BRIEF DESCRIPTION OF THE DRAWINGS

[0041] The exemplary embodiments of the invention will be described in detail, withreference to the following figures. It should be understood that the drawings are notnecessarily shown to scale. In certain instances, details which are not necessary for anunderstanding of the invention or which render other details difficult to perceive may havebeen omitted. It should be understood, of course, that the invention is not necessarily limited to the particular embodiments illustrated herein. 12 WO 2008/005690 PCT/US2007/071733 [0042] Figure 1 shows simple proportional slicing for two horizons and a dipping fault surface: (a) data prior to proportional slicing; (b) data after proportional slicing. The null zones are regions where data is not properly handled by the simple proportional slice algorithm.

[0043] Figure 2 shows simple proportional slicing for three horizons and a dipping faultsurface: (a) data prior to proportional slicing; (b) data after proportional slicing. Note that thenull data zones shift laterally between pairs of horizons due to the dip of the fault surface.[0044] Figure 3a is a flow diagram illustrating the exemplary general workflow andprocesses in performing automatic or semi-automatic interpretation of depositional systems in3-D seismic data according to this invention.

[0045] Figure 3b illustrates an exemplary system capable of performing automatic orsemi-automatic interpretation of depositional systems in 3-D seismic data according to this invention.

[0046] Figure 4 shows cross-sections of the Balmoral 3-D survey cutting across a 1 kmwide deepwater turbidite channel. The section in (a) cuts the channel at an angle of 45°. Thesection in (b) cuts the channel at an angle of 90°. The edges of the channel are indicated bythe vertical yellow arrows.

[0047] Figure 5 is a horizon slice of seismic amplitude from the Balmoral 3-D survey.

The 1 km wide deepwater turbidite channel cut by the channels in Figures 4a and 4b isreadily visible in this slice.

[0048] Figure 6 is a flow diagram illustrating the process of structural interpretation of 3- D seismic data in accordance with an exemplary embodiment of the invention.

[0049] Figure 7 shows a vertical section extracted from a 3-D seismic volume showing two interpreted horizons (labeled 1 and 2), and seven steeply dipping fault surfaces.

[0050] Figure 8 shows a schematic cross-section of four intervals (A - D), two of which13 WO 2008/005690 PCT/US2007/071733 have laterally varying thickness (C and D) in (x,y,z) space in (a), and in the transformeddomain (x,y,s) in (b). The "tick" marks illustrate the repositioning of input samples from theinput volume to the output volume.

[0051] Figure 9 shows the effects of the reefs presence as: (1) 'pull-up' of underlyinghorizons, and (2) truncation of horizontally adjacent horizons (not illustrated).

[0052] Figure 10 shows that the velocity pull-up associated with the reef has been corrected.

[0053] Figure 11 is a schematic cross-section through a simple carbonate reef. A sectionthrough the input (x,y,z) space in shown in (a) including a velocity "pull-up" of a horizon atthe base of the reef.. The domain transformed equivalent is shown in (b) with the horizon atthe base of the reef flattened, and the top reef retained. The shaded area in (b) represents nullpoints in the transformed volume.

[0054] Figure 12 illustrates a schematic cross-section through a canyon. A sectionthrough the input (x,y,z) space is shown in (a). A section through the (x,y,s) domaintransformed space is shown in (b). Note that the canyon fill B and the layers above it areyounger than the country rock A, and the paleo-depositional surfaces through the canyon fillare different that those in the surrounding country rock. The shaded area in (b) represents nullpoints in the transformed volume.

[0055] Figure 13 shows proportional slicing honoring dipping faults in 2-D for twohorizons: (a) data prior to proportional slicing; (b) data after proportional slicing. Byprojecting slices in toward their intersection with the fault surface, the null zones are reducedto a narrower dipping zone centered on the fault surface.

[0056] Figure 14 shows proportional slicing honoring dipping fault in 2-D for threehorizons: (a) prior to proportional slicing; (b) data after proportional slicing. The narrow nulldata zone is continuous between intervals and is centered on the dipping fault surface. 14 WO 2008/005690 PCT/US2007/071733 [0057] Figure 15 shows proportional slices honoring dipping faults in 3-D for twohorizons: (a) data prior to proportional slicing; (b) data after proportional slicing. A final step,a horizontal shift, is added in (b) to account for the horizontal displacement along the fault surface.

[0058] Figure 16 shows a schematic drawing of a simple faulted interval: (a) shows theinterval in the input (x,y,z) space; (b) shows the interval in the transformed domain (x,y,s).The triangular regions adjacent to the fault require special handling in the transform.

[0059] Figure 17 shows a schematic 3-D drawing of a faulted geologic volume with threelayers in (x,y,z) space in (a); and the domain transformed version of the volume in (b). Thefault has displacement down the fault surface (dip-slip) as well as a small amount ofhorizontal (strike-slip) displacement. The components of displacement are shown in the insetto (a) where a is the vertical dip-slip, β is the horizontal dip-slip, and γ is the strike-slipcomponent of displacement.

[0060] Figure 18 shows a schematic drawing of a folded structure in a section in (x,y,z)space. The solid lines show the path for a vertical interpolation. The geologic interval may bebetter represented by an interpolation that is normal to the interval bounding surfaces, shownby the dashed lines.

[0061] Figure 19 illustrates a schematic drawing of an angular unconformity (base ofinterval A) in (x,y,z) space in (a), and in the domain transformed volume in (b). The shadedarea in (b) represents null space.

[0062] Figure 20 illustrates a schematic drawing of a section through a salt body in(x,y,z) space in (a), and the corresponding section in the domain transformed space in (b).

The shaded area in (b) represents null space.

[0063] Figure 21 illustrates a flow diagram illustrating a high level view of the exemplaryprocess of domain transformation of the seismic volume from the (x,y,z or t) domain to the 15 WO 2008/005690 PCT/US2007/071733 (x,y,s) domain.

[0064] Figure 22 illustrates a flow diagram illustrating the exemplary process oftransform parameter calculation, a part of the domain transformation of seismic volumes.[0065] Figure 23 illustrates a flow diagram illustrating the exemplary process of forwarddomain transformation, a part of the domain transformation of seismic volumes.

[0066] Figure 24 illustrates an exemplary seismic volume before Domain Transformation(a), interpreted horizons and faults used in the transformation (b), and the DomainTransformed stratal-slice volume (c). The input seismic volume in (a) has deformationassociated with syn and post depositional faulting. The output Domain Transformed volume(c) is substantially free of deformation.

[0067] Figure 25 illustrates an exemplary Domain Transformed volume created from thevolume in Figure 24a, using all of the 24 faults in Figure 24b, but using only two of the fivehorizons (the top-most and bottom-most horizons). With insufficient interpretive control,there is substantial deformation remaining in the volume.

[0068] Figure 26 illustrates an exemplary flow diagram illustrating the process of refiningthe structural interpretation after initial domain transformation of the seismic volume.

[0069] Figure 27 illustrates an exemplary flow diagram illustrating the process of inversetransformation of the refined structural interpretation from the transform domain to thedomain of the original seismic volume. This process is part of the refinement of the structuralinterpretation illustrated in Figure 26.

[0070] Figure 28 illustrates a comparison of vertical and horizontal slices extracted from seismic volumes before and after domain transformation: (a) a vertical section and interpreted horizons and faults from the input seismic volume; (b) the corresponding vertical section from the domain transformed volume; (c) a horizontal slice from the input seismic volume showing a small portion of a stream channel (arrow in the lower right corner); (d) a horizontal16 WO 2008/005690 PCT/US2007/071733 slice from the domain transformed volume showing the full extent of the stream channel.[0071] Figure 29 illustrates an exemplary flow diagram illustrating the process ofstratigraphic interpretation of the domain transformed volume.

[0072] Figure 30 illustrates an exemplary defining of the initial bounding mesh (lightgrey) on a horizontal slice through the seismic volume. A channel-like object is imaged in thevolume in darker grey.

[0073] Figure 31 illustrates an example of a very simple initial bounding surface meshconsisting of two abutting cubes. The exterior faces have been tessellated, while the two interior faces have been discarded.

[0074] Figure 32 illustrates an example of an initial connected mesh defined by thesurface wrapping process. Note that this mesh is from a different example than that used inFigure 30.

[0075] Figure 33 illustrates two 2-D examples demonstrating the Surface Wrappingalgorithm, showing an initial ring of connected vertices collapsing onto: (a) two rectangularobjects; (b) a slice from an MRI (Magnetic Resonance Imaging) volume of a person's head.[0076] Figure 34 illustrates a 2-D diagram of fixed vertex determination in the SurfaceWrapping algorithm.

[0077] Figure 35 illustrates an exemplary method of determining the centralized vertexposition in the Surface Wrapping algorithm.

[0078] Figure 36 illustrates an exemplary method of determining the final vertex positionbased on the projected vertex position and centralized vertex position, using an elasticity factor of 0.8.

[0079] Figure 37 illustrates a demonstration of the progression of several iterations of thesurface wrapping algorithm, beginning with the initial bounding surface (upper left) andcontinuing to the detailed segmentation of the geobody (lower right). 17 WO 2008/005690 PCT/US2007/071733 [0080] Figure 38 illustrates determining the “sharpness” of a vertex for simulation of apermeable surface in the Surface Wrapping algorithm. The magnitude of the sum of thesurface vectors for a sharp vertex (left) is smaller than that of a blunt vertex (right).

[0081] Figure 39 illustrates an exemplary a flow diagram illustrating the process ofsurface wrapping of the elements of depositional systems. This process is part of thestratigraphic interpretation of the domain transformed volume illustrated in Figure 29.

[0082] Figure 40 illustrates an exemplary flow diagram of the process of inversetransformation from (x,y,s) to (x,y,z or t) of the surfaces and attribute volumes created in thestratigraphic interpretation process illustrated in Figures 29 and 39.

DETAILED DESCRIPTION

[0083] The exemplary embodiments of this invention will be described in relation tointerpretation of data. However, it should be appreciated, that in general, the systems andmethods of this invention will work equally well for any type of 3-D data (such as seismicdata) from any environment.

[0084] The exemplary systems and methods of this invention will also be described inrelation to seismic interpretation. However, to avoid unnecessarily obscuring the presentinvention, the following description omits well-known structures and devices that may beshown in block diagram form or otherwise summarized.

[0085] For purposes of explanation, numerous details are set forth in order to provide athorough understanding of the present invention. However, it should be appreciated that thepresent invention may be practiced in a variety of ways beyond the specific details set forth herein.

[0086] Furthermore, while the exemplary embodiments illustrated herein show thevarious components of the system collocated, it is to be appreciated that the various 18 196039/2 components of the system can be located at distant portions of a distributed network, such asa telecommunications network and/or the Internet, or within a dedicated secure, unsecuredand/or encrypted system. Thus, it should be appreciated that the components of the systemcan be combined into one or more devices or collocated on a particular node of a distributednetwork, such as a telecommunications network. As will be appreciated from the followingdescription, and for reasons of computational efficiency, the components of the system can bearranged at any location within a distributed network without affecting the operation of thesystem.

[0087] Furthermore, it should be appreciated that various links can be used to connect the elements and can be wired or wireless links, or any combination thereof, or any otherknown or later developed element(s) that is capable of supplying and/or communicating datato and from the connected elements. The term module as used herein can refer to any knownor later developed hardware, software, firmware, or combination thereof that is capable ofperforming the functionality associated with that element. The terms determine, calculate andcompute, and variations thereof, as used herein are used interchangeably and include any typeof methodology, process, mathematical operation or technique, including those performed bya system, such as an expert system or neural network.

[0089] Figure 3a shows an overall view of the CASI Workflow, the novel workflow proposed in this patent and Fig. 3b an architecture capable of performing the method. Theseismic interpretation system 300 comprises a structural interpretation module 310, astructural refinement module 320, a controller 330, a memory 340, storage 350, afiltering/conditioning module 360, a domain transformation module 370, an inverse domaintransformation module 380 and a stratigraphic interpretation module 390 which includes a 19 WO 2008/005690 PCT/US2007/071733 bounding surface module 392 and an attribute determination module 394. The functions ofthe various components of the seismic interpretation system 300 will be discussed in relationto the following figures.

[0090] In this exemplary embodiment illustrated in Fig. 3a, a processed seismic datavolume is loaded (30) into the computer for processing. This input seismic volume may havecoordinate axes that are (x,y,z) or (x,y,t), where x, y, and z are spatial dimensions (e.g., withunits of distance), where t is the measured one-way or two-way reflection time for therecorded seismic data, or where x, y, and z are simply indices incrementing or decrementingfrom initial values at the position defined as the origin of the volume.

[0091] Philosophically, the approach is based on presenting the data to an interpreter andthe computer processes in a manner that optimizes the imaging of depositional systems. Forexample, an interpreter can best recognize the existence of elements of depositional systemsby looking at slices through the data that closely approximate paleo-depositional surfaces.

The depositional elements are recognized in these slices from their characteristic morphologyor shape, and can be readily recognized even if their presence is difficult or impossible tointerpret from vertical sections of seismic data.

[0092] For example, Figures 4a and 4b show two vertical cross-sections taken from the 3-D survey. The sections cut across a 1 km wide deepwater turbidite channel. The section inFigure 4a is oriented at an angle of about 45 degrees to the direction of the channel. Thesection in Figure 4b is oriented at an angle of about 90 degrees to the channel. Evenexperienced interpreters might miss this large channel on vertical sections.

[0093] Figure 5 is a horizon slice of seismic reflection amplitude through the channel.

The existence, location, and direction of the channel are obvious from the horizon slice (inthis case it is close to a stratal-slice), and it is readily identified based on the morphology(shape) of the channel. 20 WO 2008/005690 PCT/US2007/071733 [0094] Structural Interpretation [0095] Figure 6 depicts a typical structural interpretation workflow (32). The inputseismic volume from 30 is examined by the interpreter to determine if any additional dataconditioning is required to achieve a reliable structural interpretation (decision 45). Ifdecision 45 is "Yes", then the input data volume may be filtered to remove or minimize avariety of types of noise which may improve the structural interpretation (process 46,Structural Data Conditioning). This may include processes to remove random noise, coherentnoise, or any artifacts from the volume that were introduced into or resulted from the seismicacquisition, and any processing steps preceding interpretation.

[0096] Examples of such processes would include, but are not limited to, noise filteringof the data along the z or t axis (1-D filtering), spatial filtering along the (x,y) planes (2-Dfiltering), 3-D filter operators, and any combinations of these processes. Temporal (z or t axis1-D filtering) includes, but is not limited to, high, low and band pass filtering, spectralshaping filters, and other trace filters commonly known to individuals schooled in the art ofseismic processing and interpretation. Spatial (2-D) filters include, but are not limited to,mean and median filters, spatial wavelet filtering (e.g., using a Daubechies wavelet filter),and edge preserving filtering (Al-Dossary, et. al., 2002; Jervis, 2006), and non-lineardiffusion filtering (Imhof, 2003). 2-D spatial filters may operate on the volume alonghorizontal slices, or may be guided by local estimates of structural dip in the volume. Incertain instances, the 2-D spatial filter operators may be extended into 3-D operators,depending on the type of data volume being filtered.

[0097] The above filters are all designed to reduce the level of random noise in the seismic volume. Coherent noise in the volume (e.g., "acquisition footprint" - remnant features associated with the geometry used to acquire and process the seismic data) may also need to be reduced by using a variety of coherent noise filtering techniques commonly know in the21 WO 2008/005690 PCT/US2007/071733 industry.

[0098] Once the input seismic volume (30) has been conditioned (46), or if noconditioning was necessary (i.e., the answer to decision 45 was "No"), then the interpreterproceeds to the interpretation of horizons and faults in the volume of data. The interpretationof horizons and faults may be conducted manually, accomplished using automatic processes,or by any combination of manual and automatic techniques. The interpretation of horizonsand faults may be conducted by interpreting horizons first, faults first, or by intermingling theinterpretation of horizons and faults. Thus, the process of interpreting horizons and faults areshown in parallel in Figure 6.

[0099] Decision 47 represents the decision by the interpreter to interpret faults ("Yes") ornot interpret faults ("No"). Decision 49 represents the decision by the interpreter to interprethorizons ("Yes") or not interpret horizons ("No"). If both decisions 47 and 49 are "No", thendecision 51 is "Yes" and the entire process and workflow is stopped. The workflow describedhere requires that either faults, or horizons, or both faults and horizons, be interpreted in theinput seismic volume (30) or in the conditioned seismic volume (46), or some combination of the two volumes.

[00100] If decision 47 is "Yes", then faults are interpreted from the seismic volumes (30and/or 46) using any fault interpretation technique of the interpreter's choosing - eithermanual, automatic, or a combination of manual and automatic. If decision 49 is "Yes", thenhorizons are interpreted from the seismic volumes (30 and/or 46) using any horizoninterpretation technique of the interpreter's choosing - either manual, automatic, or acombination of manual and automatic. Figure 7 shows a typical vertical seismic sectionextracted from a 3-D seismic survey with two interpreted horizons (labeled 1 and 2 in Figure7) and seven steeply dipping interpreted fault surfaces.

[00101] Once the interpreter has completed the interpretation of the horizons and faults of22 WO 2008/005690 PCT/US2007/071733 interest in process 32 shown in Figures 3a and 6, then the seismic volume, and the interpretedstructure (horizons and faults) are then passed into the Domain Transformation (process 34,Figure 3 a).

[00102] Domain Transformation [00103] Domain Transformation is a trace-by-trace approach to undoing the deformingeffects of syn-depositional and post-depositional geologic processes. Syn-depositionalprocesses take place at the same time, or as a result of, sediment deposition. Some examplesof deformation resulting from these processes are differential deposition and differentialcompaction. The term “differential” implies variation in the horizontal direction. Post-depositional processes deform the rocks present after deposition of sediments has finished.Examples of these processes are faulting and folding of sediment layers, or the rock layersthat eventually are formed by the aforementioned sediments.

[00104] All types of geologic intervals and surfaces can be accommodated including, but not limited to: • Continuous conformable intervals • Continuous intervals that exhibit growth • Intervals with reefs or carbonate platforms • Unconformities (including angular unconformities) and disconformities • Intervals with differential compaction [00105] The process may also account for post-depositional structural geologicdeformation including, but not limited to: • 3-D fault surfaces and displacement • Folding • Salt tectonics [00106] The domain-transformation algorithm requires several types of data to be input. 23 WO 2008/005690 PCT/US2007/071733

These include the seismic data volume, interpreted horizons and faults, and user inputregarding horizon types and interval types. All transform changes to be performed are storedfor each trace segment in the volume. These stored parameters consist of the starting timeand sampling rate in the original volume, as well as the storage location in the stratal-volumeand the number of sample to be interpolated during the forward transform process (62).[00107] The data volume is broken into several pieces for the purpose of DomainTransformation. There are two subdivisions used. The first is that each pair of user-suppliedhorizons defines an Interval. Each Interval may then contain one or more trace segments pertrace location (an inline and crossline intersection). The trace segments are bounded by auser-supplied horizon and either a fault or another horizon (if no fault is present in thatinterval), or by two faults.

[00108] The Domain Transformation is performed interval-by-interval through the volume.The calculation could proceed through the interpreted intervals in any order. In its preferredimplementation, the calculations proceed from the shallowest interval to the deepest interval.Within each interval, the Domain Transformation interpolates the input seismic datafollowing a set of geometric rules. The geometric rules are a function of the type of geologicinterval on which the Domain Transformation is operating.

[00109] Example Intervals and Geometric Rules [00110] The set of intervals included below is intended as a set of examples and is notinclusive of all the possible intervals that can be handled using the Domain Transformationapproach. This subset is chosen for illustrative purposes. All types of geologic intervals canbe handled using the approach presented for Domain Transformation.

[00111] Proportional Intervals: [00112] Proportional intervals include conformable intervals and growth intervals, with or without post-depositional folding and differential compaction. For continuous unfaulted 24 WO 2008/005690 PCT/US2007/071733 proportional intervals, such as in Figure 8, three steps are performed. First, a global search is performed for all trace segments contained between the two bounding horizons. This search is meant to locate the thickest portion of the interval (Zm). The number of samples in this interval (N) at its maximum thickness is calculated by dividing this maximum time thickness by the sample rate (S) of the input data volume.

N = ZM/S

[00113] In the case of Figure 8, the maximum thickness for interval C is at the right edge,and the maximum thickness for interval B is at the left edge.

[00114] For relatively shallow dips, and vertical re-sampling of the volume, the desiredsample rate for every other trace in the interval is equal to the local thickness (Zl) divided bythe maximum number of samples (N).

Sl = Zl/N

[00115] This resampling of the input seismic volume may be accomplished byinterpolation between the existing samples. The simple equations for determining the localdesired sample rate Sl guarantee that if a volume is forward transformed and then inversetransformed, none of the frequency content of the original volume will have been lost (i.e.,the output sample rate is everywhere at least as frequent as the input sample rate, therebypreserving the frequency content, and avoiding aliasing problems).

[00116] This resampling will result in thinner sections of the interval having a highersample rate in the input (x,y,z) domain than in the thicker sections. The "tick" marks on theright and left side of Figure 8a show schematically the relationship between thickness of theinput interval and the desired relative locations in z and the sample rate to produce the outputvolume with stratal slices shown in Figure 8b. Note that the intervals in Figure 8b areconstant in thickness, and that all samples are evenly spaced vertically in the output stratal slice domain. 25 WO 2008/005690 PCT/US2007/071733 [00117] All other intervals involve a generalization of this process described here for theproportional intervals.

[00118] Carbonate Platform Interval: [00119] Carbonate reefs and the intervals that are immediately overly them require specialhandling. Intervals that contain carbonate reefs represent a two-fold problem. The firstproblem is that they represent a velocity anomaly that results in a velocity 'pull-up' of theunderlying strata. The second problem is that they interrupt the horizontal continuity ofadjacent intervals (Figure 9). The transform algorithm must correct both of these issues.[00120] The first problem of a velocity pull-up is comected by handling the strataimmediately below the reef as if they were continuous flat surfaces. The second problem iscorrected by assuming that the top reef structure should remain constant in thetransformation (i.e., the shape of the top reef structure should be the same in the output stratalsliced volume as it is in the original input volume). The net result of these two corrections(shown in Figure 10) is that the base of the carbonate platform is flattened, eliminatingvelocity pull-ups, and the top structure is unchanged. The truncations of the clasticsedimentary layers against the top of the reef are maintained, and the reflections in the clasticsection around the top reef structure are flattened.

[00121] The manner of data extraction is demonstrated in Figure 11 a. The overlying interval has trace segments extracted from the top down. The actual reef has trace segments extracted from the base upward. For both intervals, the maximum number of samples is calculated from the global maximum thickness (interval A + B). The maximum number of samples for the overlying interval A is equal to the maximum thickness of A divided by the sample rate of the input volume. The maximum number of samples for the reef interval B is calculated in a similar manner, with the additional step of multiplying the initial number of samples by a velocity-contrast correction factor. This velocity-contrast correction factor is26 WO 2008/005690 PCT/US2007/071733 the ratio of the seismic velocity of the carbonate reef divided by the seismic velocity of theoverlying sediment. If this ratio is unknown, an assumed contrast (or no contrast at all) maybe used. The local number of samples for each trace segment (above or below the top reefsurface) is then calculated by multiplying that interval's maximum number of samples by theratio of the local time thickness divided by the maximum time thickness of the interval. Theresulting output section is shown Figure 1 lb. The shaded region represents a combination ofthe null regions output by both intervals, and is included in the Domain Transform outputvolume because the clastic sediments in interval A are geologically younger than the reef structure, interval B.

[00122] The result of reef correction is that continuous stratal slices can be output evenwhen they are 'cut' by a reef. Data within the reef are stretched vertically in order to correctfor the anomalous velocities within the reef (Figures 10 and lib). This technique also worksfor other types of velocity anomalies (such as a gas zone).

[00123] Canyon Intervals: [00124] Like carbonate reefs, intervals that include canyons require special handling.Although there typically is no velocity anomaly associated with the canyon, the sediment fillin the canyon is significantly younger, and belongs to different stratal slices than the "countryrock" around the canyon.

[00125] Figure 12a shows a schematic cross-section through a canyon in a seismic datavolume. Figure 12b shows the same section through a Domain Transformed volume of thecanyon. The stratal slices in the country rock (A) around the canyon and below the top of thecanyon are handled independently of the stratal slices in the fill rock (B) in the canyon andabove the top of the canyon. Because erosion had to occur to create the canyon, followed bylater deposition of (B) the stratal slices are also separated by null data values in the domain transformed volume. 27 WO 2008/005690 PCT/US2007/071733 [00126] The manner of data extraction begins by retaining the shape of the canyonunchanged in the transform. The overlying interval (B) has trace segments extracted from thetop down, including the canyon fill. The country rock through which the canyon was cut byerosion (A) has trace segments extracted from bottom up. For both intervals, the maximumnumber of samples is calculated from the global maximum thickness (interval A + B). Themaximum number of samples for the overlying interval A is equal to the maximum thicknessof A divided by the sample rate of the input volume. The maximum number of samples forthe canyon fill B is calculated in a similar manner. The local number of samples for eachtrace segment (above or below the top reef surface) is calculated by multiplying that interval'smaximum number of samples by the ratio of the local time thickness divided by themaximum time thickness of the interval. The resulting output section is shown Figure 12b.The shaded region represents a combination of the null regions output by both intervals, andis included in the Domain Transform output volume because the clastic sediments in intervalA are geologically older than the sediments in interval B. The result of canyon correction isthat continuous stratal slices can be output even when they are 'cut' by a canyon.

[00127] Faulted Intervals: [00128] A faulted interval can be treated as a continuous interval where both the upper andlower bounding surfaces are present. However, difficulties arise in the vicinity of the faultwhere only one bounding surface is defined on either side of the fault (Figures 1 and 2). Inthese fault zones, it becomes necessary to project the missing horizon inward to fill the faultzone with data points (Figures 13 and 14).

[00129] Conceptually, projection of the missing horizon is achieved by assuming that the time thickness in the fault zone is equal to the time thickness derived from the closest fully bound trace. This procedure takes place in two steps. First, an increasing radius search is performed in the (x,y) plane until the nearest trace is located that is bound by both horizons. 28 WO 2008/005690 PCT/US2007/071733

Next, the time thickness is then calculated for this full trace, and assumed to be the same for the fault zone trace. This results in a projection where the missing horizon is assumed to beequidistant from the existing horizon in the fault zone. Honoring the horizontal component ofdip-slip requires that data traces be shifted laterally in the (x,y) plane (Figure 15).

[00130] The manner of data extraction is illustrated in Figure 16 for a normal-offset fault.Other types of faults (e.g., reverse-offset faults, thrust faults, and growth faults) may beprocessed in a similar fashion. Outside of the fault zone, all trace segments may be handledas normal proportional interval trace segments. However, in the fault zone, it is not possibleto calculate the local time thickness in a normal manner as one horizon is missing (or existson the other side of the fault). In these cases, the local vertical thickness is estimated byfinding the closest fully-bounded trace segment on the same side of the fault (solid verticallines between horizons A and B in Figure 16a). The vertical thicknesses of the proximalcomplete trace segments (Zl and Zr) are then used as estimates of the local verticalthicknesses that would be present in the zones a and β were the fault not present. A numberof samples which is less than the maximum number of samples for the interval is then outputfor the local trace segment for traces in the triangular shaped regions (a and β) in Figure 16a.The local number of samples (Nl) is calculated by multiplying the maximum number ofsamples (N) by the ratio of the thickness from the horizon to the fault (Zq and Ζβ), divided by the local estimated vertical thickness. NL = Ν x (Za / ZL) and Nr = Ν χ (Ζβ / ZL) [00131] These samples are interpolated from the top down for the hanging wall of the fault, and from the bottom up for the foot wall. The output interval in Figure 16b has the same total thickness throughout. The vertical dashed lines in Figure 16b indicate the section of the interval that received special fault handling. 29 WO 2008/005690 PCT/US2007/071733 [00132] The estimation of the local vertical thickness that would be present in the zones a and β were the fault not present described above assumes a constant thickness of the intervalAB in the region of the fault. A refinement of this approach is to determine both the verticalthickness of the proximal complete trace segment, and the gradient (rate of change or firstderivative) of this thickness as the interval approaches the fault. Then, instead of projecting aconstant thickness from the proximal trace toward the fault, the estimated thickness in zones a and β would be calculated from the thickness of the proximal trace segment plus a constantgradient of that thickness.

[00133] Real faulting is commonly more complex than 2-D diagrams would imply. Figure17 is a 3-D diagram of a channel being cut by a fault. The fault exhibits both dip-slip (motionperpendicular to the long axis of the fault) and strike-slip (motion parallel to the long axis ofthe fault). Full closure of the fault requires handling both types of motion on the fault.

[00134] The algorithm described for faults above, compensates for the dip-slip componentof fault motion. The strike-slip component of motion is handled by a horizontal adjustment of the voxels in the transformed volume on one side of the fault relative to the other. The amount of the adjustment may be calculated based on a number of criteria. In its simplestform, the strike slip adjustment, if required, is the lateral displacement along the faultrequired to minimize the difference in amplitude across the fault on any given output stratalslice. This type of operation is demonstrated in Figure 17b. After transformation, the channel is continuous and unbroken.

[00135] Interpolation for Steeply Dipping Intervals [00136] The interpolation that has been described above (interpolation of traces vertically) works well for all cases except those that involve steeply dipping intervals. In steeply dipping intervals, extracting data vertically from the input volume is not a sufficient approximation of the geology present in the volume. Correct handling of steeply dipping intervals requires the30 WO 2008/005690 PCT/US2007/071733 use of non-vertical data traces, as can be seen in Figure 18. This figure represents a foldedinterval that contains a flat top and dipping flanks. The solid lines indicate vertical traces.The dashed lines indicate traces that are correctly oriented perpendicularly to the interval'sbounding horizons. The difference between the two methods of data extraction becomesgreater as dip increases on the flanks of the structure. Thus, these types of intervals requireinterpolation through the volume along paths that are non-vertical.

[00137] The path through the volume for any point in the interval A in Figure 18 may bedetermined by several methods. Surface normals may be calculated for either the upper orlower bounding surfaces of interval A. If normals are calculated for the upper surface, thenthese normals would be projected down to their intersection with the lower surface at everypoint. If normals are calculated for the lower surface, then these normals would be projectedup to their intersection with the upper surface at every point. A third, and perhaps betteralternative, is to create a surface that is mid-way between the upper and lower intervalbounding surfaces, calculate the surface normals to that intermediate surface, and extendthose normals in each direction to the upper and lower bounding surfaces at every point.[00138] Whichever method is used to calculate and project the surface normals, theprojected normals define the path of interpolation in 3-D in the interval.

[00139] Unconformities: [00140] Figure 19 demonstrates the modified handling of an interval containing anunconformity. The unconformity is the horizon at the base of interval A in the cross-section.The use of non-vertical traces (perpendicular to the bounding horizons) is shown in the pre-transform section on the left side of the figure by the dashed line. Intervals overlying anunconformity (e.g., Interval A) are treated in the normal proportional manner. Beneath theunconformity, intervals (e.g., B - E) are also handled in the normal proportional manner awayfrom the unconformity surface where both bounding horizons for each interval are present. 31 WO 2008/005690 PCT/US2007/071733 [00141] In the presence of the unconformity, the unconformity interval is handled in amanner similar to the foot wall of a fault. A search is performed to find the closest completetrace segment (both vertical and bed normal thickness are indicated in Figure 19a). This tracesegment is then used to calculate the approximate local thickness. As with fault handling, thelocal number of output samples is calculated by multiplying the maximum number ofsamples by the ratio of the thickness from the horizon to the fault divided by the localapproximate thickness. In the output section (shown on the right side of the figure), thisresults in the tapered configuration of the output interval. The shaded areas represent nullregions, which are not represented in the input volume.

[00142] Salt Boundaries: [00143] Many seismic volumes contain complex 3-D salt bodies. Intervals that arepartially bounded by salt are handled in a manner similar to faults and to reef top boundaries.As with faults, trace segments that are fully bounded by non-salt horizons are handled in thenormal proportional manner. Where a trace segment is bounded by the salt boundary and anon-salt horizon, a search must be made for the closest full non-salt bounded trace segment.As with fault handling, the thickness from this proximal fully-bounded trace segment is usedto determine the number of samples to be output in the local trace segment. This number ofsamples is equal to the maximum number of samples multiplied by the ratio of the thicknessbetween the salt and non-salt horizons divided by the thickness of the fully-bounded tracesegment. The resulting output of this interval is demonstrated in Figure 20. The shape of thesalt horizon is preserved, while removing complexities and growth that exist due to the non-salt horizons.

[00144] Domain Transformation Specification and Determination [00145] The Domain Transformation process described above is implemented as process 34 (Figure 3a) and is shown in detail in Figures 21, 22 and 23. Figure 21 shows an overview32 WO 2008/005690 PCT/US2007/071733 of the Domain Transform process 34. Figure 22 shows the detailed flow of process 60,Transform Parameter Calculation, which is part of process 34. Figure 23 shows the detailedflow of process 62, Forward Domain Transformation, which is also part of process 34.[00146] In Domain Transformation, the structural surfaces and geologic informationprovided by the interpreter regarding the types of geologic surfaces and intervals representedby the data are used to transform the seismic volume of data into a stratal-sliced volume. TheDomain Transformation ideally removes all of the effects of structural deformation of theportion of the earth represented by the seismic volume. This results in a new seismic volumewhere each horizontal slice represents a paleo-depositional surface - a surface upon whichdeposition occurred at some time in the geologic past.

[00147] The inputs to the Domain Transformation process (34) are the interpretedstructure and seismic volume(s) (conditioned or not conditioned) from process 32. There mayor may not be additional input to process 34 from process 36, Refine Structural Interpretation.Process 36 is shown in detail in Figure 26 and is described in detail below, after thedescription of process 34.

[00148] Upon input of data to process 34, decision 57 (Condition Data) is made regarding conditioning of the input data (horizons, faults, and volumes) prior to process 60 (Transform

Parameter Calculation). If decision 57 is "Yes", then the input data volume(s) and surfaces may be filtered to remove or minimize a variety of types of noise, thus improving the results of the Transform Parameter Calculation (process 60) and the results of the Forward Domain

Transform (62). This may include processes to remove random noise, coherent noise, or any artifacts from the volume that were introduced or resulted from the seismic acquisition and processing steps preceding interpretation. Examples of such processes would include, but are not limited to, mean, median or wavelet filtering to the volume, and acquisition footprint removal. It is important to note that the actual techniques used to Condition Data in process33 WO 2008/005690 PCT/US2007/071733 58 may differ from those used in process 46 to condition data for Structural Interpretation.[00149] Once the input volume(s) and surfaces (32) have been conditioned (58), or if noconditioning is necessary (i.e., the answer to decision 57 was "No"), then these data(represented in Figure 21 by 59) are input into process 60, Transform Parameter Calculation,which is described in detail in Figure 22. The Transform Parameter Calculation process (60)also requires geologic information as input from the interpreter in 63. The interpreter shouldprovide information regarding the types of geologic surfaces that are being input to process60 and the types of geologic intervals that exist between the surfaces. This "geologicknowledge" is input into the algorithm in terms commonly used by individualsknowledgeable in the practice of seismic interpretation and geologic modeling. Such surfaceswould include, but are not limited to, horizons, faults, unconformities, angular unconformities,and tops or bases of carbonate platforms. Intervals should include, but are not limited to,conformable intervals, growth intervals, and carbonate intervals.

[00150] Once all of the requisite input has been supplied, the Interval Index is initialized(64), the data for the first interval is obtained from the computer's memory (66, Get IntervalData), and the Maximum Interval Thickness is calculated (68). A Trace Segment is theportion of a seismic trace between the bounding surfaces that define an interval. The TraceSegment Index is initialized (70), and the Trace Segment is obtained from the volume.

[00151 ] All Domain Transformation operations are performed once per trace segment present in the volume. For example, in a 3-D seismic volume with two interpreted horizonsbounding one interval with no faults present, the number of trace segments will be equal tothe number of inlines present in the volume multiplied by the number of cross-lines present inthe volume. If there are three horizons present that define two unique intervals, the numberof trace segments will be twice the single-interval case. Furthermore, if the same volume hadfaults present in the intervals to be Domain Transformed, the number of trace segments 34 WO 2008/005690 PCT/US2007/071733 would increase by one for each fault at each inline and cross-line intersection that has a faultpresent inside a Domain Transformation interval.

[00152] Any seismic Trace Segment in an interval between may be cut into one or moreSub-Segments by faults. Thus, once the Trace Segment has been obtained, decision 73 ismade to determine if the trace is faulted in the interval. If the result of decision 73 is "Yes",then process 74 determines the Number of Sub-Segments into which the trace segment is cutby faults. The Trace Sub-Segment Index is initialized (76), and the Initial Trace Segment isobtained (78). Trace Sub-Segment Transform Parameters are calculated for each Sub-Segment in process 80. Decision 81 is evaluated to determine if there are More Sub-Segments on that Trace Segment in the interval. If the result of decision 81 is "Yes", the Sub-Segment Index is incremented in process 82, the next Trace Sub-Segment is obtained (78),and its Trace Sub-Segment Transform Parameters are calculated (80). This continues until allsub-segments of the trace have been processed.

[00153] If the result of decision 73 is "No", then the Trace Segment is input to process 83, and the Trace Segment Parameters are calculated. The Trace Segment and Trace Sub-Segment Transform Parameters are collected in process 84. These Transform Parametersdefine how each Trace Segment and Trace Sub-Segment must be processed in the DomainTransformation process to properly transform that segment or sub-segment given thedefinitions of the bounding geologic surfaces and the geologic interval containing that sub- segment.

[00154] Decision 85 is evaluated to determine if there are more Trace Segments in theinterval being processed. If the result of decision 85 is "Yes", then the Trace Segment Indexis incremented (86), the next Trace Segment is obtained (72), and decision 73 is evaluated forthis new Trace Segment. If the result of decision 85 is "No", then decision 87 is evaluated todetermine if there are more intervals to process. If the answer to decision 87 is "Yes", then 35 WO 2008/005690 PCT/US2007/071733 the Interval Index is incremented (88), and the Interval Data for the next interval is retrieved (66).

[00155] These processes continue until Transform Parameters have been calculated for allTrace Sub-Segments and Trace Segments in all Intervals. While calculating the tracesegments the Transform Displacement Volume is also created (89). This volume has the samedimensions as the output (re-sampled) Stratal Volume. Whereas the Stratal Volume storesthe Domain Transformed version of the input data volume, the Transform DisplacementVolume stores the x, y, and z coordinates of each data point in the Domain Transformed(x,y,s) volume. With this volume, the position of any interpretation produced from the StratalVolume can be inverse transformed from (x,y,s) to the original (x,y,z) coordinates of thesurvey. Moreover, attribute volumes calculated from the stratal sliced volume can also beinverse transformed back to the original (x,y,z) coordinates as new 3-D attribute volumes.[00156] F orward Domain T ransformation [00157] Once process 60 is complete, the data (including the seismic volume, horizons,faults, Transform Parameters, and Transform Displacement Volume) are passed from process60 to process' 62 (Figure 23). The Transform Parameters stored by process 60 comprise thestarting time and sample rate for each trace segment in the original input volume, as well asthe number of samples to be interpolated and the location to store them in the DomainTransformed (stratal) output volume. These Transform Parameters are used to build theStratal Volume by interpolating data points from the original input volume.

[00158] In process 62 (Figure 23), the Transform Parameters calculated in process 60 are applied to the seismic volume, horizons and faults to transform the seismic volume and the surfaces from the (x,y,z) or (x,y,t) domain into the (x,y,s) domain, where s, the vertical dimension of the transformed data signifies "stratal-slice." The first two steps in process 62 are Initialize the Interval Index (90) and Initialize the Trace Index in the Interval (92). 36 WO 2008/005690 PCT/US2007/071733

Process 94 then retrieves the Local Trace Segment Parameters. Process 96 then calculates theLocal Transformed Trace segment. Decision 97 is then made to determine if there are MoreTraces to be processed in the Interval. If the result of decision 97 is "Yes", then the TraceIndex is incremented (98) and the Transform Parameters are retrieved for the next LocalTrace Segment (94). If the result of decision 97 is "No", then decision 99 is evaluated todetermine if there are more intervals to be processed. If the result of decision 99 is "Yes", theInterval Index is incremented (100), and the Trace Index in the Interval is initialized (92).Process 62 continues in this manner until all trace segments in all intervals have beentransformed to the (x,y,s) domain.

[00159] The type of interpolation performed by process 96 may be one of many established interpolation techniques. These techniques include, but are not limited to, linearinterpolation, spline interpolation, and "sine function" interpolation (also known as “(sinx)/x” interpolation). If interpolation is being performed in the trace (z or t) direction, thepreferred implementation would use sine interpolation. In the case of interpolation along a 3-D path in the volume (e.g., because of steep dips), some combination of techniques may beused with the horizontal and vertical parts of the interpolation operation handled separately.[00160] The Transform Parameter Calculation process (60) outputs all of the inputsurfaces in Stratal-Domain coordinates. All types of surfaces are output as listed previously.Horizons are output as planar features between the intervals that they separate in the StratalDomain. Reef Tops are output with the same shape as they were input, but with theirposition following the interface between the extracted reef values and the null values thatexist above the valid values. Similarly, salt boundaries and erosive surfaces are output alongthe interface of the valid data they bound and the null regions they define. These transformedsurfaces act as a cue for the interpreter, indicating how the Stratal Domain volume relates to the input volume. 37 WO 2008/005690 PCT/US2007/071733 [00161] When sufficient structural control in the form of interpreted horizons, faults, saltboundaries, canyon boundaries, and other possible geological surfaces is used, the DomainTransformation produces a stratal-slice volume that is substantially free of any deformation.This deformation may have been caused by post or syn-depositional folding or faulting,differential compaction and/or differential sedimentation.

[00162] Figure 24a shows an input seismic volume in (x,y,t) space. The set of horizonsand faults that were interpreted in the volume are shown in Fig. 24b, and the stratal-slicevolume output by the Domain Transformation process in shown in Fig. 24c. A total of fivehorizons and 24 faults were used in this transformation. The input volume shows substantialdeformation from faulting, and from differential sedimentation (note the increasing thicknessbetween interpreted horizons moving from the left edge of the image to the right edge.[00163] The output volume is substantially free of deformation, in that there are nosignificant remnant effects from the faulting or from the differential sedimentation. The reflection events in the stratal-sliced volume are all flat.

[00164] Refine Structural Interpretation [00165] One key feature of the Domain Transformed volume is that if there are any errorsor omissions in the interpretation of horizons or faults in process 32 (Figure 6), those errorsand omissions are highlighted or emphasized in the Domain Transformed volume that isoutput from process 34.

[00166] Figure 25 shows a Domain Transformed volume, created from the input volume shown in Figure 24a using all of the interpreted faults, but only two of the five interpretedhorizons shown in Fig. 24b. This Domain Transformed volume is not a stratal-sliced volume,in that it still retains a substantial amount of deformation from both faulting and from differential sedimentation. This deformation is most evident in the middle of the volume at the points that are farthest from the two bounding horizon surfaces that were used in the 38 WO 2008/005690 PCT/US2007/071733 transformation. This volume clearly requires refinement of the structural interpretation (inclusion of additional interpreted horizons in this case). With insufficient interpretive control, there is substantial deformation remaining in the volume.

[00167] In decision 35 (Figure 3a), the Domain Transformed volume is examined for anyof these errors and/or omissions. If Structural Refinement is needed (if decision 35 is "Yes")then process 36 is used to Refine the Structural Interpretation, thus correcting those errorsand/or omissions. After this refinement process, the Domain Transformation (process 34)must be applied again in order for the refined structural interpretation to be applied in thetransformation process.

[00168] The details of the Refine Structural Interpretation process (36) are shown in

Figure 26. If decision (35) has been evaluated as "Yes", that Structural Refinement is needed,process 36 begins with a decision (101) to determine whether to Refine Structural in theTransformed Volume or in the original seismic volume. If decision 101 is evaluated as "No",then control is passed to process 32 Structural Interpretation, and the structural refinementoccurs in the original seismic volume using process 32. If decision 101 is evaluated as "Yes",the structural refinement will be performed in process 36 on the domain transformed seismic volume.

[00169] Data passed into process 36 from process 34 and decision 35 includes the domaintransformed volume, the domain transformed surfaces (horizons, faults, etc.), and theTransform Displacement Volume. Structural refinement may involve the interpretation ofadditional horizons and/or faults, and may also include editing or changing the horizons orfaults that were initially interpreted in process 32. This interpretation and editing may beconducted manually, accomplished using automatic processes, or by any combination ofmanual and automatic techniques. The interpretation of horizons and faults may be conductedby interpreting horizons first, faults first, or by intermingling the interpretation of horizons 39 WO 2008/005690 PCT/US2007/071733 and faults. Thus, the processes of refining the interpretation of horizons and faults are shownin parallel in Figure 26.

[00170] Decision 103 represents the decision by the interpreter to interpret faults ("Yes")or not interpret faults ("No"). Decision 105 represents the decision by the interpreter tointerpret horizons ("Yes") or not interpret horizons ("No"). If both decisions 103 and 105 are"No", then decision 107 is "Yes" and the refinement process is stopped, and the transformeddata and workflow passes on to process 38, Stratigraphic Interpretation. The workflowdescribed here does not require that the interpretation be refined, even if a need for additionalrefinement of the structural interpretation is indicated by decision 35.

[00171] If decision 103 is "Yes", then faults are interpreted and/or edited in the domaintransformed seismic volumes (from 35) using any fault interpretation technique of theinterpreter's choosing - either manual, automatic, or a combination of manual and automatic.If decision 105 is "Yes", then horizons are interpreted and/or in the domain transformedseismic volumes (from 35) using any horizon interpretation technique of the interpreter'schoosing - either manual, automatic, or a combination of manual and automatic.

[00172] Inverse Transform Refined Structural Interpretation [00173] Once the interpreter has completed the refinement of the interpretation of thehorizons and faults of interest in processes 108 and/or 110 in Figure 26, then the domaintransformed horizons and faults, and the refined structure (horizons and faults) interpreted inthe domain transformed volume are then passed (111) into the Inverse Transform RefinedStructural Interpretation process (process 112).

[00174] The refined structural interpretation, if it is performed in the transformed volume in process 36, must be inverse transformed and merged with the original structural interpretation in process 32. The Inverse Transform of the Refined Structural Interpretation (process 112) is shown in detail in Figure 27. 40 WO 2008/005690 PCT/US2007/071733 [00175] The input to process 112 is the refined interpretation of faults and horizons (111), and the Transform Displacement Volume. Process 114 initializes the Surface Index, and the first Surface is obtained by process 116. Process 118 initializes the Point Index for points on the selected surface, and the first Point on the Surface is obtained by process 120. Process 122 is then used to convert the coordinates of the Point on the Surface from the transformed or stratal slice domain to the corresponding coordinates in the original domain of the seismic volume from 30. Decision 123 is then evaluated to determine whether there are More Points on the Surface. If decision 123 is evaluated as "Yes", then the Point Index is incremented byprocess 124, and the next Point on the Surface is obtained by process 120. If decision 123 isevaluated as "No" (there are no more points on the current surface), then decision 125 is then evaluated to determine whether there are more surfaces. If decision 123 is evaluated as "Yes"then process 126 increments the Surface Index, and the next Surface is obtained by process116. If decision 123 is evaluated as "No", then process 112 is completed and passes controlback to process 34, to repeat the Domain Transformation of the original seismic volumeusing the refined structural interpretation.

[00176] Stratigraphic Interpretation [00177] Once the interpreter decides that no more Structural Refinement is required (i.e., decision 35 in Figure 3a is "No") either because the structural interpretation is complete, or because the interpreter has chosen to proceed with an incomplete structural interpretation, then process 38, Stratigraphic Interpretation, is applied to the Domain Transformed volume.

The purpose of the Stratigraphic Interpretation process is to assist the interpreter in the identification and interpretation of elements of depositional systems, or other seismic stratigraphic features represented in the Domain Transformed volume. The identification of these elements of depositional systems is accomplished by calculation of a variety of seismic attribute volumes within process 38. Once these elements of depositional systems have been41 WO 2008/005690 PCT/US2007/071733 identified in the attribute volumes, the bounding surfaces of these elements are also created inprocess 38.

[00178] Figure 28 shows a comparison of slices extracted from an input seismic volumeand a domain transformed volume. Five horizons and more than twenty faults were used inthe domain transformation process applied to this volume. Figure 28a shows a vertical sectionextracted from the input seismic volume, with the fault surfaces and bounding horizons thatintersect that section. Figure 28b shows the corresponding section extracted from the domaintransformed volume. The arrows on Figure 28a show the points that correspond to the four comers of the section extracted from the domain transformed volume.

[00179] Figures 28c and 28d respectively show 3-D views of horizontal slices takenthrough the input seismic volume and the domain transformed volume. The arrows indicate achannel on the two slices. In the input volume, only a small portion of the channel is visibleon the horizontal slice because of the growth and faulting that is present in the volume. The entire channel is visible on the horizontal slice from the domain transformed volume because the effects of growth and faulting have been removed by the domain transformation process.

[00180] The Domain Transformed seismic volume and Domain Transformed interpreted surfaces are input to process 38 (Stratigraphic Interpretation) from process 34 and decision 35 (Figure 3a). Process 38 is shown in detail in Figure 29. Once data are input into process 38, decision 127 is evaluated to determine if the transformed data (both the transformed seismic volume and the interpreted surfaces) need to be conditioned prior to stratigraphic interpretation. If decision 127 is "Yes", then the input transformed data volume and surfaces may be filtered to remove or minimize a variety of types of noise which may improve the stratigraphic interpretation (process 128, Stratigraphic Data Conditioning). This conditioning may include processes to remove random noise, coherent noise, or any artifacts from the volume that were introduced or resulted from the seismic acquisition and processing steps, or42 WO 2008/005690 PCT/US2007/071733 any noise or artifacts introduced in the Domain Transformation process (34). Examples ofsuch processes would include, but are not limited to, mean, median or wavelet filtering to thevolume, and acquisition footprint removal.

[00181] Once the input domain transformed seismic volume and surfaces have beenconditioned (128), or if no conditioning was necessary (i.e., the answer to decision 127 is"No"), then a number of Stratigraphic Attribute Volumes may be calculated in process 130.The goal of calculating these attribute volumes is to produce through a single volume, orthrough a combination of volumes, a data volume or volumes that provide improved imagingof depositional systems when compared to the domain transformed seismic volume.Individuals practiced in the art of stratigraphic interpretation from 3-D seismic data are familiar with these attribute volumes.

[00182] Attribute imaging of stratigraphy is improved by first transforming the seismicvolume, and by then calculating the attribute volume in the transformed domain. This can beseen when compared to the typical practice of calculating the attributes directly on the inputseismic volume without using domain transformation. Because of this, it is advantageous toprovide the attribute volumes created by this workflow in the transform domain as outputfrom process 38 to be inverse transformed in process 40, as is shown in Figure 29.

[00183] Once the Stratigraphic Attribute Volumes have been calculated by applying process 130 to the transformed seismic volume, decision 131 is evaluated to determinewhether Multi-Attribute Imaging is going to be used to aid in the imaging of depositionalsystems using the attribute volumes. If decision 131 is evaluated as "Yes", then process 132 isapplied to identify optimum combinations of attributes to image the el ements of depositionalsystems in the transformed space. There are a number of techniques in the industry thatwould be familiar to one who is practiced in the art of stratigraphic interpretation. Theseinclude, but are not limited to, neural network and neural network related techniques to 43 WO 2008/005690 PCT/TJS2007/071733 analyze combinations of attributes for clusters that might identify elements of depositionalsystems (e.g., Kohonen Self-Organizing Maps and Growing Neural Gas), direct clusteranalysis techniques (e.g., K-Means clustering), and techniques such as attribute cross-plotmatrices and multi-dimensional attribute crossplot visualization techniques. Any of thesetechniques could be used in process 132 to image and analyze the multiple attribute volumescreated by process 130 to image the depositional systems in the data.

[00184] Surface Wrapping [00185] As elements of depositional systems are imaged in individual attribute volumes, or combinations of attribute volumes, these elements are then "interpreted" by obtaining theirbounding surfaces (i.e., the surface that completely encloses the depositional feature).Decision 135 is evaluated to determine how the bounding surfaces of the depositionalsystems will be interpreted. If decision 135 is evaluated as "No", the interpreter will manuallyinterpret the bounding surface of the depositional system or depositional system elementmanually (136). If decision 135 is evaluated as "Yes", then process 142 (Surface Wrapping)will be used to semi-automatically interpret the bounding surface of the depositional systemor depositional system element.

[00186] Surface Draping (Dorn, 1999) is an effective technique for creating polygonalrepresentations of surfaces that are essentially planar, but it cannot be applied directly to theproblem of finding the bounding surface of a geobody. The Surface-Wrapping algorithmcreates a 3-D polygonal mesh that entirely surrounds 3-D objects.

[00187] Surface Wrapping is a semi-automatic approach for segmentation of a geobodybounding surface within volumetric data. The approach is metaphorically based upon theconcept of collapsing an elastic surface onto a physical object. The desired output of theprocess is a polygonal mesh that may be stored as data, displayed to the end user, or used infurther data processing techniques. This approach has advantages over fully automated 44 WO 2008/005690 PCT/US2007/071733 segmentation algorithms in that it may be applied to data where the volume to be segmentedis not fully imaged, or where a high level of noise is present. This approach is alsosignificantly less time consuming for the human analyst than fully manual segmentationtechniques, in that the user need only define an approximate initial bounding surface prior toapplication of the algorithm which determines a more detailed and accurate bounding surface.[00188] In the preferred embodiment, an interpreter first specifies a range of voxel valuesthat best isolates the voxels that correspond to the boundary of the intended geobody in thevolume. The interpreter then defines an initial 3-D bounding surface which completelyencloses the intended geobody and approximates its contours, isolating the voxels belongingto the geobody boundary from the rest of the volume.

[00189] The initial bounding surface may be constructed using manual, automatic, orsemi-automatic methods, or any combination thereof. In the preferred embodiment, themethod for defining the initial bounding surface is based on a technique described by Kobbeltet. al. (1999). In this method, which is similar to graphical user interfaces that are commonlyfound in fully manual volume segmentation software, one slice of the volume is displayed onthe screen, and the user defines the region contained by the initial bounding surface using avirtual brush to “paint” the region on the screen, as shown in Figure 30. Though similar tomanual segmentation interfaces, the interface used in Surface Wrapping is different in twoways. First, the painted region must fully enclose the boundary of the intended geobody(alternatively it must almost fill the boundary of the intended geobody), but need notprecisely track the contour of the volume. Second, the brush defines the same 2-D region on auser-defined range of slices simultaneously, thus extending the approximate bounding region into 3-D.

[00190] Internally, the painted region is represented as a collection of cubes of equal dimensions, where each cube corresponds to a small portion of the volume that is contained45 WO 2008/005690 PCT/US2007/071733 within the initial bounding surface. Figure 31 shows two such adjacent cubes. The size of thecubes is globally adjustable by the user; the smaller the cube size, the denser the boundingsurface mesh. To construct the initial bounding surface mesh, the hull of the painted region isfound by discarding cube faces that share the same spatial coordinates. The remaining facesare tessellated into two triangles per face, which collectively form the visible polygons of thebounding surface mesh (Figure 31). An example of an initial mesh constructed from a largenumber of cubes is shown in Figure 32.

[00191] In order to allow the Surface Wrapping algorithm to process the bounding surfacemesh as if it were an elastic material, each vertex in the mesh maintains a record of itsneighboring vertices, where a neighboring vertex is defined as any vertex to which it isdirectly connected by an edge of a triangle. Each vertex also maintains a record of alltriangles of which it is a part. Vertex locations correspond to index coordinates relative to thedata volume, and there may be at most one vertex data structure in the mesh at any givenspatial coordinate, thus ensuring connectivity of vertices over the entire mesh.

[00192] Once the initial bounding surface has been defined, the Surface Wrapping processiteratively moves each vertex in the mesh toward the boundary of the intended geobody, asillustrated in 2-D in Figures 33a and 33b, wherein an initial ring of connected verticescollapses onto: (33a) two rectangular objects; (33b) a slice from an MRI (MagneticResonance Imaging) volume of a person's head.

[00193] The process by which the vertices are moved is detailed below.

[00194] Each iteration of the surface wrapping algorithm begins with the calculation of the outward vertex normal vector for the first vertex in the mesh. The vertex normal is calculated as the normalized mean of the adjacent face normals, with a unit length corresponding to the grid spacing of the voxels in the data volume. A face normal, N„ is calculated as the cross product of vectors A and B, which are vectors lying on the local "plane" of the surface. 46 WO 2008/005690 PCT/US2007/071733

Ns = AxB

[00195] The vertex normal Nw is calculated as the ratio:

Nw = N' [00196] where the vector Nv is defined as the sum of face normals NSi divided by the totalnumber of adjacent faces, n: Σν-

Nv= ——n [00197] and where |NV| is the magnitude of vector Nv.

[00198] If the initial mesh has been created everywhere outside the object whose boundaryis sought, a projected location for that vertex is then calculated to be at a point one unit lengthfrom the vertex's current position in the direction opposite to the outward unit normal at thatvertex. If the initial mesh has been created everywhere inside the object whose boundary issought, a projected location for that vertex is then calculated to be at a point one unit lengthfrom the vertex's current position in the direction of the outward unit normal at that vertex.[00199] If the voxel value at the projected location falls within the range specified by theinterpreter as corresponding to the boundary of the body being wrapped, then the vertex isflagged as "fixed" (Figure 34) and the projected location is not recorded. If the voxel value atthe projected location does not fall within the specified range, the projected location is storedin the vertex data structure. This process is repeated for each vertex in the mesh, and is notorder-dependent.

[00200] Following the calculation of the projected location for each non-fixed vertex, asecond location is computed for each non-fixed vertex, referred to here as the centralizedlocation. For a given vertex, the centralized location is determined to be the mean of thecurrent locations of its neighboring vertices, as illustrated in Figure 35. This process is 47 WO 2008/005690 PCT/US2007/071733 repeated for each vertex in the mesh, and is not order-dependent.

[00201] When both the projected and centralized locations have been computed for a non-fixed vertex, its actual updated location corresponds to a point partway along the linesegment between the projected location and centralized location. The proximity of theupdated location to either end of the line segment is determined by a user-adjustable elasticityfactor, which is defined as a percentage of the distance from the centralized location to theprojected location. A higher elasticity factor causes the mesh to be treated as a more pliantmaterial, while a lower elasticity factor simulates the effect of increased surface tension.Figure 36 shows a 2-D illustration of calculating a final vertex location based on a vertex’sprojected location and centralized location, using an elasticity factor of 0.8.

[00202] If at this point the mesh must be processed further, as determined eitherautomatically or by the interpreter, the above process is then repeated until the interpreter issatisfied that the desired bounding surface has been achieved. Figure 37 illustrates theprogression from top left to bottom right of successive iterations of this process as used tocreate a bounding surface mesh of a salt dome. In order to provide the user with real-timefeedback for the result of successive iterations of the Surface Wrapping algorithm, vertexpositions at each iteration may be pre-calculated prior to yielding control of the graphicalinterface to the user. Using, for example, a scroll bar, a user can reveal the results of variousvertex calculations and display them graphically. For example, by adjusting a scroll bar, auser could be presented with the series of images displayed in Fig. 37.

[00203] An additional feature that may be incorporated into the Surface Wrappingalgorithm is the simulation of a semi-permeable surface, which allows outlying voxels to“push through” the mesh while maintaining the overall desired structure of the boundingsurface. In the preferred implementation, this is accomplished by the use of an additionalpredicate immediately prior to the calculation of a vertex’s projected location which 48 WO 2008/005690 PCT/US2007/071733 determines if the vertex is creating a sharp point in the mesh. This predicate sums the surfacenormal vectors of the triangles that connect to the vertex and calculates the magnitude of theresulting vector. This calculation is illustrated in 2-D in Figure 38, which shows thedifference in the magnitude of the summed surface vectors for a sharp vertex versus a bluntvertex. If the magnitude is below a user-definable threshold, the vertex is flagged as non-fixed for the current iteration of the Surface Wrapping algorithm, which has the effect ofsmoothing out the spike in the mesh.

[00204] The Surface Wrapping algorithm may also be applied to a subset of the vertices inthe bounding surface mesh, allowing localized editing operations. A typical mechanism forselecting the affected vertices is any picking operation in a 3-D rendered view of the mesh,but the selection of vertices can be accomplished via any combination of manual or automatictechniques.

[00205] Surface Wrapping (process 142) is shown in detail in Figure 39. In SurfaceWrapping, the interpreter Creates an Initial Bounding Surface Mesh (146), which is a veryapproximate connected three-dimensional mesh which either completely surrounds thedepositonal system or system element of interest (Case A), or is completely included withinthe depositional system or system element of interest (Case B). The Surface Wrappingprocess will shrink the initial mesh inward to the boundary of the depositional system orsystem element in Case A, or will expand the initial mesh outward to the boundary of thedepositional system or system element in Case B.

[00206] Once the Initial Bounding Surface Mesh has been created by the interpreter,process 148 Initializes the Vertex Index, and process 150 gets the initial Mesh Vertex.

Process 152 calculates the Projected Vertex Location based on a movement of the vertex from its initial position along the direction of its unit normal toward the element of the depositional system. Process 154 then calculates an estimate of the Voxel Value at the 49 WO 2008/005690 PCT/US2007/071733

Projected Vertex Location. Decision 155 (Voxel Value in Range) is evaluated to determine ifthe voxel has encountered the boundary of the depositional system or system element. If theVoxel Value is in the interpreter specified range (decision 155 evaluated as "Yes") then thevertex is Flagged as Fixed by process 158. If the Voxel Value is outside the interpreterspecified range (decision 155 evaluated as "No") then process 156 Moves the Vertex to theProjected Location.

[00207] After either process 156 or 158 is complete, decision 159 is evaluated todetermine if there are More Vertices. If decision 159 evaluates as "Yes", then process 160Increments the Vertex Index, and process 150 gets the next Mesh Vertex. If decision 159evaluates as "No", then process 162 re-initializes the Vertex Index, and process 164 Gets the

Mesh Vertex. Decision 165 is then evaluated to determine if the Vertex is Fixed. If decision 165 is evaluated as "Yes" (i.e., the vertex has been flagged as fixed) then process 168Increments the Vertex Index, and process 164 Gets the next Mesh Vertex. If decision 165 isevaluated as "No" (i.e. the vertex has not been flagged as fixed) the process 166 Centralizesthe Vertex Location with respect to the neighboring vertices in the mesh.

[00208] After process 166, decision 167 is evaluated to determine whether there are MoreVertices in the mesh. If decision 167 is evaluated as "Yes" (i.e., there are more vertices), thenprocess 168 Increments the Vertex Index, and process 164 Gets the next Mesh Vertex. Ifdecision 167 is evaluated as "No" (i.e., there are no more vertices), then decision 169 isevaluated to determine whether to Continue Shrinking the Mesh. If decision 169 is evaluatedas "Yes", then process 148 Initializes the Vertex Index for the next shrinking step, andprocess 150 gets the Mesh Vertex. Note that once a vertex has been flagged as fixed inprocess 158, its projected position in process 154 remains fixed. If decision 169 is evaluatedas "No", then the surface mesh represents the bounding surface of the depositional system orsystem element to the conditions set by the interpreter at the beginning of the Surface 50 WO 2008/005690 PCT/US2007/071733

Wrapping process (142).

[00209] Referring back to Figure 29, once the bounding surface or surfaces have beenobtained using process 136 or process 142, then the resulting surface or surfaces are outputfrom process 38 to process 40.

[00210] Inverse Domain Transformation [00211] Referring to Figure 3a, as individual elements of depositional systems areidentified, interpreted and their bounding surfaces are created in process 38, the attributevolumes, and/or the bounding surfaces are transformed back into the spatial domain of theinput seismic volume (30) by applying process 40, the Inverse Domain Transformation.Alternatively, the interpreter may choose not to apply process 40 until several, or all, of theelements of depositional systems have been identified and interpreted. Process 40 requires asinput information from the original seismic volume (30), from the Domain Transformation(34 through 35), and from the Stratigraphic Interpretation process (38).

[00212] Data input into process 40 include the Domain Transformed Volume, attributevolumes calculated from the Domain Transformed Volume, the Domain TransformedSurfaces, all interpreted stratigraphic surfaces (either from manual interpretation or fromSurface Wrapping), and the Transform Displacement Volume.

[00213] The Inverse Transform process (40) is shown in detail in Figure 40. Both attribute volumes and surfaces may be inverse transformed from the stratal slice domain into the spatial domain of the original input seismic volume (30). The workflows for surfaces and volumes are shown in parallel in Figure 40 for process 40. Upon entry into process 40, decisions 171 and 173 are evaluated. If decision 171 (Inverse Transform Surfaces) is evaluated as "Yes", then the workflow for Inverse Transform of Surfaces is invoked. If decision 173 (Inverse Transform Volumes) is evaluated as "Yes", then the workflow for

Inverse Transform Volumes is invoked. If both decision 171 and decision 173 are evaluated51 WO 2008/005690 PCT/US2007/071733 as "No", then all desired surfaces and volumes have been inverse transformed, and control is passed to process 42.

[00214] The Inverse Transform process (40) allows interpretation produced in the StratalDomain to be inverted back to the coordinates of the original volume. It also allowsattributes that are produced with superior quality in the Stratal Domain to be inverted to theoriginal coordinates. In process 38, if decision 171 is evaluated as "Yes" (Inverse TransformSurfaces), each point of the surfaces in the interpretation is inverted by finding its nearestneighbors in the Transform Displacement Volume. The original positions (in spatialcoordinates) of these nearest neighbors are used to invert the position of that point. Inprocess 38, if decision 173 is evaluated as "Yes", (Inverse Transform Volumes), each trace isre-sampled (stretched) back to the original coordinate using similar interpolation schemes to those described for the Forward Transform.

[00215] If decision 171 (Inverse Transform Surfaces) is evaluated as "Yes", then process 174 initializes the Surface Index, and the first Surface is obtained by process 176. Process178 initializes the Point Index for points on the selected surface, and the first Point on theSurface is obtained by process 180. Process 182 is then used to convert the coordinates of thePoint on the Surface from the transformed or stratal slice domain to the correspondingcoordinates in the original domain of the seismic volume from 30. Decision 181 is then evaluated to determine whether there are More Points on the Surface. If decision 181 is evaluated as "Yes", then the Point Index is incremented by process 184, and the next Point onthe Surface is obtained by process 180. If decision 181 is evaluated as "No" (there are nomore points on the current surface), then decision 183 is then evaluated to determine whetherthere are more surfaces. If decision 183 is evaluated as "Yes", then process 186 incrementsthe Surface Index, and the next Surface is obtained by process 176. If decision 183 isevaluated as "No", the workflow for inverse transformation of surfaces is completed. 52 WO 2008/005690 PCT/US2007/071733 [00216] If decision 173 (Inverse Transform Volumes) is evaluated as "Yes", then process 188 initializes the Volume Index, and the first Volume is obtained by process 190. Process192 initializes the Voxel Index for points in the selected volume, and the first Voxel in theVolume is obtained by process 194. Process 196, with additional input from process 35, Getsthe Stratal Voxel Spatial Coordinates. Process 198 Gets the Stratal Voxel Value, and process200, with additional input from process 30, Gets the Spatial Voxel Target Position in theinverse transformed vol ume. The results of processes 196, 198, and 200 are input into process202, which then calculates the Voxel Value at the Target Position.

[00217] Decision 201 is then evaluated to determine whether there are More Voxels in theVolume. If decision 201 is evaluated as "Yes", then the Voxel Index is incremented byprocess 204, and the next Voxel in the Volume is obtained by process 194. If decision 201 isevaluated as "No" (there are no more voxels in the current volume), then decision 203 is then evaluated to determine whether there are more volumes. If decision 203 is evaluated as "Yes"then process 206 increments the Volume Index, and the next Volume is obtained by process190. If decision 203 is evaluated as "No", the workflow for inverse transformation ofvolumes is completed.

[00218] Referring to Figure 3, the output of the workflow (42) may include seismicattribute volumes, which were created in process 38 in the Domain Transformed coordinatesystem, and which have been Inverse Transformed by process 40 into the coordinate space ofthe input seismic volume (30). The output of the workflow (42) may also include thebounding surfaces of the interpreted elements of depositional systems, transformed byprocess 40 from surfaces in the coordinate space of the Domain Transformed volume, to thecoordinate space of the input seismic volume (30).

[00219] While the above-described flowcharts have been discussed in relation to a particular sequence of events, it should be appreciated that changes to this sequence can 53 WO 2008/005690 PCT/US2007/071733 occur without materially effecting the operation of the invention. Additionally, the exact sequence of events need not occur as set forth in the exemplary embodiments. Additionally, the exemplary techniques illustrated herein are not limited to the specifically illustrated embodiments but can also be utilized with the other exemplary embodiments and each described feature is individually and separately claimable.

[00220] The systems, methods and techniques of this invention can be implemented on aspecial purpose computer, a programmed microprocessor or microcontroller and peripheralintegrated circuit element(s), an ASIC or other integrated circuit, a digital signal processor, ahard-wired electronic or logic circuit such as discrete element circuit, a programmable logicdevice such as PLD, PLA, FPGA, PAL, any means, or the like. In general, any devicecapable of implementing a state machine that is in turn capable of implementing themethodology illustrated herein can be used to implement the various communication methodsand techniques according to this invention.

[00221] Furthermore, the disclosed methods may be readily implemented in softwareusing object or object-oriented software development environments that provide portablesource code that can be used on a variety of computer or workstation platforms.

Alternatively, the disclosed system may be implemented partially or fully in hardware usingstandard logic circuits or VLSI design. Whether software or hardware is used to implementthe systems in accordance with this invention is dependent on the speed and/or efficiencyrequirements of the system, the particular function, and the particular software or hardwaresystems or microprocessor or microcomputer systems being utilized. The systems, methodsand techniques illustrated herein can be readily implemented in hardware and/or softwareusing any known or later developed systems or structures, devices and/or software by thoseof ordinary skill in the applicable art from the functional description provided herein and witha general basic knowledge of the computer and geologic arts. 54 196039/2 [00222] Moreover, the disclosed methods may be readily implemented in software that canbe stored on a storage medium, executed on programmed general-purpose computer with thecooperation of a controller and memory, a special purpose computer, a microprocessor, or thelike. The systems and methods of this invention can be implemented as program embedded onpersonal computer such as an applet, JAVA® or CGI script, in C or C++, Fortran, or the like, asa resource residing on a server or computer workstation, as a routine embedded in a dedicatedsystem or system component, or the like. The system can also be implemented by physicallyincorporating the system and/or method into a software and/or hardware system, such as thehardware and software systems of a dedicated seismic interpretation device. 55 , DTixan rwzn ova inia^a pnow pnszn irn πτ qaoa,ρνιη ηχΰπ laoana ma’na nairmaa np’ioa .zruwan rwaa mpnan ρηυ1? oxnmοιηπη iv

<img img-format="tif" img-content="drawing" file="IL196039AD00021.tif" id="idf0001" />

.(moia παΊηπ) cros&amp;amp;'an rwa