Feature detection in seismic volumes

Methods, systems, and computer-readable media for analyzing a domain are provided. The method includes defining a mask plane that includes a first dimension of a first number of voxels and a second dimension of a second number of voxels, and selecting a plurality of first angles for orientating the mask plane in the domain with respect to a first axis. The method also includes for each one of the plurality of first angles selected populating, using one or more processors, sum cubes associated with each one of a first plurality of subject voxels. The method also includes selecting a plurality of second angles, and for each one of the plurality of second angles selected, calculating a planar sum for each one of a second plurality of subject voxels selected.

BACKGROUND

Seismic modeling and simulation of subterranean formations may employ a process in which collected seismic data may be used to construct a two-dimensional or three-dimensional model of the formation. In such processes, the seismic data may be in the form of properties of signals that are received after being transmitted through an area of interest, for example, around a borehole, vertically and/or angled from the surface (e.g., land, the ocean floor, etc.), or in other locations. The signal characteristics may include direction, time taken to propagate through the area, amplitude, among a variety of others possible. This information may be combined with knowledge of rock properties and/or the subterranean formation to link these properties to the structure of the formation. Accordingly, the properties of a number of signals may be combined to construct an approximation of the structure of the formation.

The recognition of certain features, e.g., faults, horizons, etc., may be of particular interest in such modeling. However, in some cases, the area of interest, number of signals, number of features, or a combination thereof may be large, such that it would be useful to employ computers to automatically detect, enhance, and extract information about the features. One way this is done is by using a windowed Radon transform. Briefly, a system may discretize the area of interest into pixels (two dimensions) or voxels (three dimensions). The system may then sum the signal property (e.g., intensity) of the pixel or voxel with those of a certain number of its neighbors in a range of directions. If the intensity sum for the subject pixel or voxel is higher in one direction, it may be flagged as part of a possible feature. Next, the system may look for the next pixel or voxel that is part of the same feature or may simply move on to the next pixel or voxel in sequence. Such techniques may employ a “brute force” method, such that a line or plane extending from the subject pixel or voxel sweeps across a range of angles, calculating the intensity sum for each angle, and then moving on to the next pixel or voxel and repeating the process.

SUMMARY

Embodiments of the present disclosure may include systems, methods, and computer-readable mediums for analyzing a domain, e.g., detecting features in the domain, particularly a discretized domain. In one embodiment, the method may employ an efficient version of a windowed Radon transform to rapidly calculate a sum of all energies in a mask plane. The method may precalculate intermediate sums in a certain azimuth angle of values associated with each voxel, for example, using running sums to capitalize on redundancies. The method may then calculate the sum of the intermediate sums for each dip angle in a range of dip angles, such that the sum of the intermediate sums for each dip angle may be conceptually equivalent to a sum of voxels in a plane oriented at the azimuth angle and each of the dip angles. The method may also capitalize on redundancies in the planar sum calculation by using running sums. The method may repeat this until each selected azimuth angle and dip angle combination is analyzed for each selected voxel, thereby completing a three-dimensional the windowed Radon transform with improved efficiency.

This summary is provided to introduce some of the subject matter described below and is not to be considered limiting.

DETAILED DESCRIPTION

The following detailed description refers to the accompanying drawings. Wherever convenient, the same reference numbers are used in the drawings and the following description to refer to the same or similar parts. While several embodiments and features of the present disclosure are described herein, modifications, adaptations, and other implementations are possible, without departing from the spirit and scope of the present disclosure. Accordingly, the following detailed description does not limit the present disclosure. Instead, the proper scope of the disclosure is defined by the appended claims.

In the field of seismic analysis, aspects of a geologic environment may be defined as attributes or “features.” In general, seismic attributes may help to condition seismic amplitude data for improved structural interpretation tasks, such as determining the location of lithological terminations and helping to isolate hidden seismic stratigraphic features of a geologic environment. Attribute analysis may assist in defining a trap in exploration or delineating and characterizing a reservoir at the appraisal and development phase, for example. An attribute generation process may rely on a library of various seismic attributes (e.g., for display and use with seismic interpretation and reservoir characterization workflows). In some cases, generation of attributes on the fly for rapid analysis may be conducted. In other cases, attribute generation may be conducted a background process (e.g., a lower priority thread in a multithreaded computing environment), which may allow for execution of one or more foreground processes (e.g., to enable a user to continue using various components).

Attributes may help extract value from seismic and other data, for example, by providing more detail on subtle lithological variations of a geologic environment (e.g., an environment that includes one or more reservoirs). In general, an accurate reconstruction of paleostress may be difficult to achieve for a geologic environment. In particular, stress magnitudes may be difficult to reconstruct based on borehole data (e.g., as acquired over a field grid). Stress magnitudes may be helpful to understand and exploit resources in reserves such as carbonate reserves.

Carbonate rocks, however, may present a challenge for petrophysical analysis. Carbonates may be deposited primarily through biological activity where the resulting rock composition (e.g., of fossil fragments and other grains of widely varying morphology) produces highly complex pore shapes and sizes. Carbonate mineral species may also be comparatively unstable and subjected to multiple stages of dissolution, precipitation, and recrystallization, adding further complexity to the porosity and permeability of the rocks. Further, relationships between depositional attributes, porosity, and permeability may be obscured by such physical, biological, and chemical influences, operating at different scales, during and continuing after deposition.

According to an embodiment of the method disclosed herein, a “windowed” Radon transform may be employed to enhance identification of features such as cracks, faults, discontinuities, etc., which, in turn, may provide for a more comprehensive understanding of a reservoir environment. While various examples are described with respect to analysis of seismic data, techniques may be used for other types of data especially where implementations of a windowed Radon transform may benefit from reduction in requirements for computational resources. In the windowed Radon transform analysis, runtime may depend upon, e.g., image size, mask (or window) dimensions, and the number of angles, combinations thereof and/or others.

Turning to the implementations and aspects described herein for illustrative purposes,FIG. 1depicts a method100for analyzing a domain, e.g., detecting features in the domain, according to an embodiment. The method100may, for example, provide a process for performing a windowed Radon transform with a reduced runtime.

The method100may begin by receiving seismic data corresponding to a domain, as at102. For example, the seismic data may include intensity, amplitude, velocity, etc., of any type of seismic wave (e.g., compression waves, shear waves, etc.). The seismic data may then be employed to generate a model of the domain. The domain may be three-dimensional, having a length, width, and depth, corresponding to a subterranean area of interest, for example. Further, the domain may be discretized into a plurality of voxels, i.e., smaller cubes (or any other suitable three-dimensional shape), as at104, which each represent a portion of the domain. Each voxel may be associated with information, for example, corresponding to one or more values associated with the seismic data for the volume that the voxel occupies in the domain.

With continuing reference toFIG. 1,FIG. 2illustrates a conceptual view of a domain200and a mask plane202for use in performing a windowed Radon transform analysis of the domain200. The domain200may be discretized into voxels; however, for the sake of clearly illustrating the mask plane202, the voxels are not shown inFIG. 2. The method100may include determining first and second dimensions for the mask plane202, as at106(FIG. 1). The first dimension may be measured by number of voxels D1of the mask plane202and another number of voxels D2may define the second dimension. The numbers of voxels D1, D2may depend on the size and complexity of the domain200and/or the features contained therein, along with the level of granularity of the model that is suitable for analysis. The numbers of voxels D1, D2may be the same or different. Further, the number of voxels D1, D2may be chosen so as to minimize the number of operations for the windowed Radon transform analysis, which may increase for higher numbers of voxels D1, D2, while providing statistically relevant information.

The mask plane202may have an orientation defined by an azimuth angles θ and a dip angle φ. The azimuth angles θ may define the angle between the x-axis and the mask plane202, and the dip angle φ may define the angle between the z-axis and the mask plane202. Accordingly, for a full three-dimensional windowed Radon transform by the “brute force” method, sums from each voxel intersected by the mask plane202in each dip angle φ and each azimuth angles θ combination are provided. To reduce the time complexity associated with such an analysis, the number of mask planes202may be limited, e.g., such that sequential planes associated with a given voxel are separated by at least one voxel at their outer extremities, thereby, for example, providing a unique set of voxels intersected by the mask plane202in each orientation. In one example implementation, D1×D2×16 planes202may be provided for each voxel. To fully calculate the windowed Radon transform for an individual voxel using the “brute force” method, values for each of the D1×D2voxels in each mask plane202are summed. Since there are, e.g., D1×D2×16 mask planes202per voxel, the total number of operations to calculate the sum of each of the values of the voxels intersected by the mask plane202may be D12×D22×16. In some cases, the windowed Radon transform may include calculating this sum for each voxel of the domain200. If the domain has dimensions of N1voxels by N2voxels by N3voxels, the total number of operations may be N1×N2×N3×D12×D22×16 to perform the windowed Radon transform for the domain200.

In various embodiments, the method100may advance over the brute force method. For example, the method100may include calculating sum cubes for at least some of the voxels of the domain200, for example, by separating the analysis for azimuth angles and dip angles for each mask plane202.FIG. 3illustrates the domain in discretized form, labeled as300, according to an embodiment. As shown, the discretized domain300may be discretized into a plurality of voxels302at104(FIG. 1). The voxels302may be aligned in layers (four shown:303-1,303-2,303-3, and303-4). Each layer303-1,303-2,303-3,303-4may extend parallel to a common reference plane, for example, the x-y plane.

Referring again toFIG. 1, the method100may proceed to selecting an azimuth angle θ for the mask plane202(FIG. 2) for analysis, as at108. The method100may, in an embodiment, regard the voxels302individually. Further, in the illustrated example, the selected azimuth angle θ1may be 90°.

For purposes of discussion, the voxel302currently being analyzed is referred to herein as a “subject” voxel. Accordingly, the method100may proceed to selecting a subject voxel from among the plurality of voxels302, disposed in a layer, as at110. For example, as shown inFIG. 3, the method100may include selecting a subject voxel304disposed in layer303-1. Further, in the domain300, the first number of voxels D1determined at106, that is, the length of the mask plane202in the x-y plane (seeFIG. 2), may be three.

With the azimuth angle θ1selected at108and the subject voxel304selected at110, the method100may, in an embodiment, proceed to defining an azimuth line, as at112.FIG. 3further illustrates an azimuth line306extending at the first azimuth angle θ1(90°), with the azimuth line306intersecting the subject voxel304. In an embodiment, the azimuth line306may be centered on the subject voxel304, as shown; however, in other embodiments, the azimuth line306may intersect the subject voxel304at a non-center location.

Further, the azimuth line306may extend the first number D1of voxels302in the same layer303-1as the subject voxel304. Accordingly, in at least one embodiment in which the azimuth line306is centered on the subject voxel304and extends a length of three voxels (D1), the azimuth line306may intersect the subject voxel304and two additional voxels308,310. The voxels that are intersected by the azimuth line306may be referred to as “neighbor” voxels to the subject voxel304. Neighbor voxels may be any number of voxels away from the subject voxel304, for example, with each voxel therebetween also being a neighbor voxel.

The method may then proceed to populating a sum cube with an intermediate sum associated with the subject voxel304and the azimuth line306, as at114. In an embodiment, calculating the intermediate sum for the subject voxel304may proceed by summing values associated with the subject voxel304and the other voxels308,310(e.g., neighbor voxels) that are intersected by the azimuth line306. In one example, the values may be associated with seismic data contained in the voxels304,308,310. Further, this intermediate sum may be stored in the sum cube, i.e., indexed to the location of the subject voxel304and the azimuth angle θ1of the azimuth line306, thereby populating the sum cube for the subject voxel304.

The method100may determine whether to calculate intermediate sums for additional subject voxels302, as at116. If the final subject voxel for calculation has not yet been reached, the method100may loop back to selecting a (next) subject voxel, as at110. In an embodiment, this subsequent selection at110may shift the method100to considering the voxel310as the next subject voxel, as illustrated inFIG. 4.

The voxel310may be chosen as the next subject voxel at least because it is adjacent to the previous subject voxel304and is intersected by the azimuth line306. The analysis of the voxel310may be similar to the analysis of the subject voxel304, but may be simplified by capitalizing on redundancy using a “running sum” technique. That is, an azimuth line312may be defined, as at114, centered on the subject voxel310and extending in the layer303-1for the first number D1(e.g., three) of voxels302, as shown. The azimuth line312may extend at the selected first azimuth angle θ1, as the azimuth line306did. Further, the subject voxel310may be displaced from the previous subject voxel304by a number N of voxels302, e.g. one voxel (i.e., voxels304and310are adjacent), in the direction of the azimuth lines306,312. Accordingly, the azimuth lines306and312may be overlapping, intersecting D1-N of the same voxels302. For the illustrated embodiment in which N is one and D1is three, the azimuth line306extends to intersect voxel308, but not voxel316, while azimuth line312intersects voxel316but not308. Thus, the azimuth line306intersects two voxels, voxels304and310, which the azimuth line306also intersects.

Thus, shifting consideration to the subject voxel310after the subject voxel304may yield all the same voxels302intersected by the azimuth lines306,312, except N number of voxels302at the end (or beginning) of the azimuth lines306,312. Thus, rather than summing all D1(three) values associated with the D1intersected voxels304,310,316to arrive at the intermediate sum for the subject voxel310in the azimuth angle θ1of the azimuth line312, calculating the intermediate sum for populating the sum cube at114may employ a “running sum” technique.

For example, the running sum may initially be the sum calculated for the previous subject voxel304. From this sum, the value associated with the voxel308that is intersected by the azimuth line306for the previous subject voxel304, but not intersected by the azimuth line312of the current subject voxel310may be subtracted from the running sum. Further, the value associated with the voxel316that is intersected by the azimuth line312, but not by the azimuth line306may be added to the running sum. The order in which this subtraction and addition occurs may be switched. Further, in embodiments where N is greater than one, a corresponding number of additions and subtractions may also be employed to keep a running sum.

As such, instead of performing the first number D1(e.g., three) of operations to calculate the sums for the first number D1of voxels304,310,314intersected by the azimuth line312, the method100performs 2×N (one subtraction and one addition for each non-overlapping voxel) operations to calculate the same sum. In cases where the first number D1is large (e.g., 10, 100, 1000, 10,000, or more), this may result in a large number of operations being avoided. Further, in cases where many (e.g., millions) of intermediate sums are calculated, even a reduction of one operation per calculation may reduce runtime.

In this manner, proceeding from one subject voxel to another, employing the same first azimuth angle θ1and capitalizing on redundancies using the running sums, the method100may populate the sum cube with intermediate sums for each (or a selected subset) of the voxels302corresponding to sums for the azimuth angle θ1. As such, the method100may iterate through blocks110-116until the sum cubes are populated with intermediate sums for one, some, or a select subset of all the voxels302and using azimuth angle θ1as the direction for the azimuth lines.

The method100may also include determining planar values for each, or a select subset, of the voxels302with the azimuth angle θ1and a range of dip angles φ (i.e., the angle of the mask plane202with respect to the z-axis).FIG. 5illustrates an example case of such determining, according to an embodiment. The range of dip angles φ may be −90° and +90°; however, a smaller subset of this range may be prescribed, for example, employing dip angle φ guidance, using predetermined information to determine what dip angles φ to avoid, etc. Accordingly, the method100may proceed to selecting a dip angle φ1in the range of dip angles φ, as at118. For example, in the illustrated case, the selected dip angle φ1is zero degrees. The method100may then proceed to again selecting a subject voxel for analysis, as at120. In the example case, voxel404is selected as the subject voxel.

As shown, an embodiment of the method100also includes defining a dip line400for the subject voxel402, as at122. The dip line400may extend at the selected dip angle φ1. With the selected dip angle φ being zero in this case, the dip line400is parallel to the z-axis. The dip line400may also extend perpendicular to the layers303-1,303-2,303-3,303-4at each selected dip angle φ1, φ2, . . . φm. Further, the dip line400may extend by the second number D2of voxels302, as defined by the mask plane202at106. In this case, the second number D2may be three.

The dip line400may intersect voxels402,410, and412. Populating the sum cubes at114(e.g., previously) with each of the intermediate sums may include populating the sum cubes with the values of the voxels302intersected by azimuth lines404,406,408, as shown. That is, each of the voxels402,410,412may have been selected as subject voxels at110for the azimuth angle θ1. Thus, summing or otherwise calculating a combination of these intermediate sums, as at124, may be effectively summing all of the voxels302intersected by the mask plane202when it is centered on (or otherwise associated with) the subject voxel402and oriented at dip angle θ1(0°) and the azimuth angle θ1(90°). This may yield a “planar sum” associated with the mask plane202having the prescribed orientation and dimensions and associated with the subject voxel selected at120.

The method100may then determine whether additional subject voxels are to be analyzed, as at126. When it is determined that there is one or more additional subject voxels to be selected and analyzed, the method100may proceed to iterating back to selecting a (next) subject voxel410, as at120.FIG. 6illustrates the analysis for the next subject voxel410, according to an embodiment. The method100may again define a dip line414intersecting the subject voxel410, and calculating the planar sum for the mask plane202centered thereat, or otherwise associated therewith, and oriented at the at the first selected dip angle φ1(0°) and the first selected azimuth angle θ1(90°).

As with calculating the intermediate sums using the azimuth lines described above with respect toFIG. 4, the method100may again capitalize on redundancy by using running sums to calculate the planar sums at124. Accordingly, the voxel410may be chosen as the next subject voxel at120since it is adjacent to the previous subject voxel402and intersected by the dip line400. A dip line414may be defined, extending the second number D2of voxels302(e.g., three) at the first selected dip angle φ1, thus intersecting voxels402and410, as well as a voxel416. As such, the dip line414overlaps the dip line400by one less than the second number D2of voxels302(e.g., the overlap is for two of the three voxels302since, in this case, the second number D2is three). Accordingly, the sum for previous subject voxel402may be subtracted by the intermediate sum associated with the voxel412that is intersected by the dip line400but not intersected by the dip line414, and added to the intermediate sum associated with the voxel416that is intersected by the dip line414and not intersected by the dip line400. In various embodiments, sequential subject voxels may be displaced by more than one voxel, as discussed above for calculating the intermediate sums.

The intermediate sums for each of the voxels402,410,416may have already been calculated at iterations passing though the block114ofFIG. 1with voxels402,410, and416acting as subject voxels in turn and defining azimuth lines406,404,418, respectively. Thus, the sum for the dip line414is the sum for the mask plane202centered at (or otherwise associated with) the subject voxel410and oriented at the dip angle φ1of 0° and the azimuth angle θ1of 90°. This sequence of blocks120-124, e.g., using running sums, may be repeated for each of the selected subject voxels.

After summing one, a select subset, or all of the voxels302in the domain300using the dip angle φ1and the azimuth angle θ1, the method100may determine that no additional subject voxels are to be analyzed with these parameters, at126. The method100may then proceed to determining whether another dip angle φ is to be employed as a parameter for calculating planar sums for subject voxels with the azimuth angle θ1, as at128. When this is determined, the method may loop back to selecting a dip angle φ at118. In this example, a second dip angle φ2may be chosen, for example, according to a preset dip angle interval, dip angle guidance, etc.

FIG. 7illustrates the analysis for the subject voxels using the second selected dip angle φ2, for example, of about 45°, as a parameter. Such calculation may be accomplished in a similar fashion as described for the first selected dip angle φ1. Accordingly, the method100may include selecting a subject voxel422at120and defining a dip line420intersecting, e.g., centered, at the subject voxel422, as at122. The dip line420may extend at the second dip angle φ2by the second number D2of voxels302. The intermediate sums for each of the voxels302(including subject voxel422) intersected by the dip line420may then be summed, to result in the planar sum for the mask plane202centered at or otherwise associated with the subject voxel422and oriented at the first azimuth angle θ1(e.g., 90°) and the second dip angle φ2(e.g., 45°), as at124. The calculation of planar sums using the second selected dip angle φ2(as well as subsequent selected dip angles φn) may also incorporate the running sum technique described above.

This process outlined in blocks118-126may be repeated, for example, each time the method100includes determining that additional dip angles φ are to be employed as a parameter for analysis with the first azimuth angle θ1at128. Eventually, the method100may include determining at128that no further dip angles φ are to be used as a parameter with the first selected azimuth angle θ1.

At this point, for example, the method100may determine whether additional azimuth angles θ are to be employed for analysis, as at130. When the method100includes determining that additional azimuth angles θ are to be employed for analysis, the method100may return to selecting one of the additional azimuth angles θ for analysis. In at least one embodiment, a range for azimuth angles θ may be 0° to 180° or a subset thereof limited based on any one of a variety of factors, such as a priori knowledge of the subterranean formation. For each azimuth angle θ selected, the portions of the method100outlined in blocks110-128may repeat until the planar sums for the selected subject voxels are determined for each combination of the selected azimuth angle θ1, θ2, . . . , θ, and the selected dip angle φ1, φ2, . . . , φm. Once this is determined, the method100may perform one or more various other analyses or may conclude.

An example of pseudo-code for this embodiment of the method100may be substantially as follows:

Define the first dimension D1;Define the second dimension D2;For each azimuth angle θ:Calculate running horizontal sums of D1 voxels runningthrough the domain in direction θ;Save the intermediate results in a sum cube;For each dip angle φ:Calculate running sums of D2 voxels through the sumcube with dip angle φ, and perpendicular to the direction ofthe horizontal sums.

In another embodiment, the order of analysis for the azimuth angles θ and the dip angles φ may be reversed. As such, the dip angle φmmay be the “first” angle, i.e., the first angle selected in the method100, and the intermediate sums may be calculated by proceeding through the selected (or all) voxels302at the dip angle φm. Thereafter, the azimuth angles θnmay be employed as the “second angle” with the dip angle φmremaining constant, thereby providing for a calculation of the planar sums. In other words, the iterations through the azimuth angles θ1-nmay be determined for each dip angle φ1-m.

Thus, the number of operations needed for calculating the discrete windowed Radon transform with D1×D2×16 local planes for each voxel consistent with an embodiment of the present disclosure may be D1×8+D1×D2×32 per voxel. This may compare well with the “brute force” method discussed above, in which the time complexity may be D12×D22×16 for each voxel. In some cases, this transition to a lower time complexity may enable a three-dimensional windowed Radon transform to be interactive or semi-interactive (i.e., able to be manipulated on demand) rather than impractical or very slow (i.e., creating a lag time of minutes, hours, days, etc.).

In at least one embodiment, the mask plane202may also include a thickness τ, also measured in numbers of voxels. This may promote detection of nearly, but not completely, planar distributions of energy in the domain300. Moreover, the mask plane202with thickness τ may be estimated while minimizing additional computing time relative to a mask plane202without a thickness. For example, prior to populating the sum cubes with the intermediate sums at114, the method100may include calculating the local sums of voxel sub-cubes having a dimension, for example, of τ×D3×D4, where D3and D4may be any suitable number, for example, 1 each or τ each.

Further, the method100may rescale the domain300or any sub-cube thereof to at least partially distribute data uniformly. For example, the method100may include assessing the statistical significance of the sums if the distribution of the input is known (uniform in this example case). This is because sums of random uniformly distributed data may be normally distributed. Thus, the method100may include a data rescale transform that may be performed to generate a rescaled cube that subsequently may be processed as described above for the domain300. In another embodiment, the method100may apply such data rescale transform “on the fly” as a preconditioning step for one or more calculations that are performed.

Additionally, the populating the sum cube114and/or calculating the planar sum at124may or may not proceed as simple summations of values associated with the selected voxels302. For example, such calculations of the intermediate and/or planar sums may proceed by incorporating one or more additional calculations, statistical analyses, etc., to promote accuracy, precision, and/or conservation of computing resources. For example, the calculations may include parallel or simultaneous calculation of several statistics, such as the sum of squared values associated with each of the voxels302. In another example, fault line detection may include choosing preferred line segments based on both semblance (e.g., mean) and normalized variance.

Moreover, the method100may include dip-field estimation, for example, using data known about the domain300. Any suitable dip field estimation technique may be employed to steer the search for faults in orientations deviating significantly from the stratigraphic layers, especially, for example, in areas where variation in dip angles φ is too large to narrow the dip range globally.

More than one output value may be defined from the method100, such as the azimuth angles θ and dip angle φ of the fault or other feature. Such additional information may be used to guide automated fault extraction, e.g., by indicating in which direction to shift for analyzing subsequent voxels and/or limiting or informing the selection of subsequent dip angles φ and/or subsequent azimuth angles θ. In another embodiment, given fault seed points or a seed fault stick interpretation and an edge attribute seismic volume as input, the method100may provide robust feature (e.g., fault) auto-tracking by determining the approximate azimuth angles θ and dip angle φ of the feature, and subsequently calculating the three-dimensional transform using a narrow range of dip angles φ and azimuth angles θ to inform the selection of subsequent dip angles φ and/or subsequent azimuth angles θ.

Thus, the method100may include estimating the partial derivatives of a fault plane in three dimensions. The calculation of the steering may be efficient, since the transform may be calculated locally and for a narrow range of parameters, e.g., for a small part of the domain300. The method100may also promote the tracking staying on the path with the “highest energy” in the edge cube, and may be robust as to high noise-to-signal ratios if large operator sizes are used.

Accordingly, in various embodiments, the method100may be employed for improved, three-dimensional feature detection. For example, the method100may include finding the maximum mask planar sum or average over all azimuth angles θ and all dip angles φ as the result output. This planar sum may correspond to “energy” in the plane, with high planar energy indicating the presence of a feature. Thus, the algorithm may be employed in a variety of contexts to search for extreme high or low energy in a domain of voxels, depending on the nature of the input data (e.g. seismic edge attribute), and thereby detect features.

Moreover,FIG. 8illustrates one example, among many contemplated, of a system800in which one or more embodiments of the method100may be implemented. The system800includes various management components810to manage various aspects of a geologic environment850, for which an embodiment of the method100may be employed to model and/or detect/track features. For example, the management components810may allow for direct or indirect management of sensing, drilling, injecting, extracting, etc. with respect to the geologic environment850. In turn, further information about the geologic environment850may become available as feedback860(e.g., optionally as input to one or more of the management components810).

In the example illustrated inFIG. 8, the management components810include a seismic data component812, an information component814, a processing component816, a simulation component820, an attribute component830, an analysis/visualization component842, and a workflow component844. In operation, according to an embodiment, seismic data and other information provided per the components812and814may be inputted to the simulation component820, for example, after processing using the processing component816, which may be configured to implement one or more embodiments of the method100. The simulation component820may process information to conform to one or more attributes, for example, as specified by the attribute component830, which may be a library of attributes. Such processing may occur prior to input to the simulation component820(e.g., per the processing component816). Further, the simulation component820may perform operations on input information based on one or more attributes specified by the attribute component830. The simulation component820may construct one or more models of the geologic environment850, which may be relied upon to simulate behavior of the geologic environment850(e.g., response to one or more acts). In the example ofFIG. 8, the analysis/visualization component842may allow for interaction with a model or model-based results. Additionally, output from the simulation component820may be inputted to one or more other workflows, as indicated by the workflow component844.

Further, the management components810may include features of a commercially-available simulation framework such as the PETREL® seismic to simulation software framework (Schlumberger Limited, Houston, Tex.). The PETREL® framework may provide components that allow for optimization of exploration and development operations. Further, the PETREL® framework may include seismic to simulation software components that may output information for use in increasing reservoir performance, for example, by improving asset team productivity. Through the use of such a framework, various professionals (e.g., geophysicists, geologists, and reservoir engineers) may develop collaborative workflows and integrate operations to streamline processes.

The management components810may include features for geology and geological modeling to generate high-resolution geological models of reservoir structure and stratigraphy (e.g., classification and estimation, facies modeling, well correlation, surface imaging, structural and fault analysis, well path design, data analysis, fracture modeling, workflow editing, uncertainty and optimization modeling, petrophysical modeling, etc.). Particular features may allow for performance of rapid three-dimensional seismic interpretation, for example, for integration with geological and engineering tools (e.g., classification and estimation, well-path design, seismic interpretation, seismic attribute analysis, seismic sampling, seismic volume rendering, geobody extraction, domain conversion, etc.). In reservoir engineering, for a generated model, one or more features may allow for simulation workflow to perform streamline simulation, reduce uncertainty and assist in future well planning (e.g., uncertainty analysis and optimization workflow, well path design, advanced gridding and upscaling, history match analysis, etc.). The management components810may include features for drilling workflows including well path design, drilling visualization, and real-time model updates (e.g., via real-time data links).

Various aspects of the management components810may be add-ons or plug-ins that operate according to specifications of a framework environment. For example, a commercially available framework environment marketed as the OCEAN® framework environment (Schlumberger Limited) allows for seamless integration of add-ons (or plug-ins) into a PETREL® framework workflow. AntTracker™ may be one such add-on or plug-in that may provide one or more aspects of edge detection using the three-dimensional, windowed Radon transform techniques as part of embodiments of the method100.

Embodiments of the disclosure may also include one or more processor (i.e., computing) systems for implementing one or more embodiments of the method100.FIG. 9illustrates a schematic view of such a processor system900, according to an embodiment. The processor system900may include one or more processors902of varying core (including multiple core) configurations and clock frequencies. The one or more processors902may be operable to execute instructions, apply logic, etc. It will be appreciated that these functions may be provided by multiple processors or multiple cores on a single chip operating in parallel and/or communicably linked together.

The processor system900may also include a memory system, which may be or include one or more memory devices and/or computer-readable media904of varying physical dimensions, accessibility, storage capacities, etc. such as flash drives, hard drives, disks, random access memory, etc., for storing data, such as images, files, and program instructions for execution by the processor902. In an embodiment, the computer-readable media904may store instructions that, when executed by the processor902, are configured to cause the processor system900to perform operations. For example, execution of such instructions may cause the processor system900to implement one or more portions and/or embodiments of the method100described above.

The processor system900may also include one or more network interfaces908. The network interfaces908may include any hardware, applications, and/or other software. Accordingly, the network interfaces908may include Ethernet adapters, wireless transceivers, PCI interfaces, and/or serial network components, for communicating over wired or wireless media using protocols, such as Ethernet, wireless Ethernet, etc.

The processor system900may further include one or more peripheral interfaces906, for communication with a display screen, projector, keyboards, mice, touchpads, sensors, other types of input and/or output peripherals, and/or the like. In some implementations, the components of processor system900need not be enclosed within a single enclosure or even located in close proximity to one another, but in other implementations, the components and/or others may be provided in a single enclosure.

The memory device904may be physically or logically arranged or configured to store data on one or more storage devices910. The storage device910may include one or more file systems or databases in any suitable format. The storage device910may also include one or more software programs912, which may contain interpretable or executable instructions for performing one or more of the disclosed processes. When requested by the processor902, one or more of the software programs912, or a portion thereof, may be loaded from the storage devices910to the memory devices904for execution by the processor902.

Those skilled in the art will appreciate that the above-described componentry is merely one example of a hardware configuration, as the processor system900may include any type of hardware components, including any necessary accompanying firmware or software, for performing the disclosed implementations. The processor system900may also be implemented in part or in whole by electronic circuit components or processors, such as application-specific integrated circuits (ASICs) or field-programmable gate arrays (FPGAs).

The foregoing description of several possible embodiments has been presented for purposes of illustration only. It is not exhaustive and does not limit the present disclosure to the precise form disclosed. Those skilled in the art will appreciate from the foregoing description that modifications and variations are possible in light of the above teachings or may be acquired from practicing the disclosed embodiments.

For example, the same techniques described herein with reference to the processor system900may be used to execute programs according to instructions received from another program or from another computing system altogether. Similarly, commands may be received, executed, and their output returned entirely within the processing and/or memory of the processor system900. Accordingly, neither a visual interface command terminal nor any terminal at all is strictly necessary for performing the described embodiments.

Likewise, the steps described need not be performed in the same sequence discussed or with the same degree of separation. Various steps may be omitted, repeated, combined, or divided, as necessary to achieve the same or similar objectives or enhancements. Accordingly, the present disclosure is not limited to the above-described embodiments, but instead is defined by the appended claims in light of their full scope of equivalents.