Anti-aliasing and Anti-leakage regularization of high-dimensional seismic data

Abstract

Methods and systems may include receiving input seismic data and determining an input spectrum from the input seismic data. The input spectrum includes an input amplitude at each spectral element. The methods and systems may further include, iteratively or recursively, generating, using a weight function, weight data from the input spectrum, generating a weighted input spectrum by weighting the input spectrum using the weight data, identifying neighboring spectral elements in the weighted input spectrum, updating an output spectrum by replacing an output amplitude at each neighboring spectral element with the input amplitude at each of the neighboring spectral elements in the input spectrum, and updating the input spectrum by removing the input amplitude at each neighboring spectral element. The methods and systems may still further include determining output seismic data from the updated output spectrum and generating a seismic image based, at least in part, on the output seismic data.

Description

ANTI-ALIASING AND ANTI-LEAKAGE REGULARIZATION OF HIGH-DIMENSIONAL SEISMIC DATABACKGROUND

[0001] In the oil and gas industry, a seismic survey may be used to characterize a subterranean region of interest in the search for hydrocarbon reservoirs. During the seismic survey, seismic waves are emitted into the subterranean region of interest and, in a land environment, the vibrations of the earth caused by the seismic waves are detected and recorded as seismic data. The seismic survey may collect irregular, and sometimes sparse, seismic data.

[0002] Following seismic processing of the seismic data, the seismic data may be used to generate a seismic image of the subterranean region of interest that characterizes the subterranean region of interest. The seismic image may then be used to locate hydrocarbon reservoirs within the subterranean region of interest. In turn, a wellbore path may be planned and drilled to penetrate a located hydrocarbon reservoir to ultimately produce hydrocarbons to the surface for use.SUMMARY

[0003] This summary is provided to introduce a selection of concepts that are further described below in the detailed description. This summary is not intended to identify key or essential features of the claimed subject matter, nor is it intended to be used as an aid in limiting the scope of the claimed subject matter.

[0004] In general, in one aspect, embodiments relate to a method. The method includes initializing an output spectrum, receiving, from a seismic acquisition system, input seismic data for a subterranean region of interest, and determining, using a transform, an input spectrum from the input seismic data. The input seismic data includes toxels. The input spectrum includes an input amplitude at each of spectral elements. The method further includes, iteratively or recursively, until a convergence criterion is satisfied, generating,  using a weight function, weight data from the input spectrum, generating a weighted input spectrum by weighting the input spectrum using the weight data, identifying neighboring spectral elements among the spectral elements in the weighted input spectrum, updating the output spectrum by replacing an output amplitude at each of the neighboring spectral elements in the output spectrum with the input amplitude at each of the neighboring spectral elements in the input spectrum, and updating the input spectrum by removing the input amplitude at each of the neighboring spectral elements in the input spectrum. The spectral element among the neighboring spectral elements includes a maximum amplitude. The method still further includes determining, using an inverse transform, output seismic data from the updated output spectrum, generating a seismic image of the subterranean region of interest based, at least in part, on the output seismic data, and identifying a drilling target of a hydrocarbon reservoir within the subterranean region of interest using, at least in part, the seismic image.

[0005] In general, in one aspect, embodiments relate to a system. The system includes a seismic processing system and seismic interpretation workstation. The seismic processing system is configured to initialize an output spectrum, receive, from a seismic acquisition system, input seismic data for a subterranean region of interest, and determine, using a transform, an input spectrum from the input seismic data. The input seismic data includes toxels. The input spectrum includes an input amplitude at each of spectral elements. The seismic processing system is further configured to, iteratively or recursively, until a convergence criterion is satisfied, generate, using a weight function, weight data from the input spectrum, generate a weighted input spectrum by weighting the input spectrum using the weight data, identify neighboring spectral elements among the spectral elements in the weighted input spectrum, update the output spectrum by replacing an output amplitude at each of the neighboring spectral elements in the output spectrum with the input amplitude at each of the neighboring spectral elements in the input spectrum, and update the input spectrum by removing the input amplitude at each of the neighboring spectral elements in the input spectrum. The spectral element among the neighboring spectral elements includes a maximum amplitude. The seismic processing system is still further configured to determine, using an inverse transform, output seismic data from the updated output  spectrum and generate a seismic image of the subterranean region of interest based, at least in part, on the output seismic data. The seismic interpretation workstation is configured to identify a drilling target of a hydrocarbon reservoir within the subterranean region of interest using, at least in part, the seismic image.

[0006] In general, in one aspect, embodiments relate to a non-transitory computer-readable memory having computer-executable instructions stored thereon that, when executed by a computer processor, perform steps. The steps include initializing an output spectrum, receiving, from a seismic acquisition system, input seismic data for a subterranean region of interest, and determining, using a transform, an input spectrum from the input seismic data. The input seismic data includes toxels. The input spectrum includes an input amplitude at each of spectral elements. The steps further include, iteratively or recursively, until a convergence criterion is satisfied, generating, using a weight function, weight data from the input spectrum, generating a weighted input spectrum by weighting the input spectrum using the weight data, identifying neighboring spectral elements among the spectral elements in the weighted input spectrum, updating the output spectrum by replacing an output amplitude at each of the neighboring spectral elements in the output spectrum with the input amplitude at each of the neighboring spectral elements in the input spectrum, and updating the input spectrum by removing the input amplitude at each of the neighboring spectral elements in the input spectrum. The spectral element among the neighboring spectral elements includes a maximum amplitude. The steps still further include determining, using an inverse transform, output seismic data from the updated output spectrum, generating a seismic image of the subterranean region of interest based, at least in part, on the output seismic data, and identifying a drilling target of a hydrocarbon reservoir within the subterranean region of interest using, at least in part, the seismic image.

[0007] Other aspects and advantages of the claimed subject matter will be apparent from the following description and the appended claims.BRIEF DESCRIPTION OF DRAWINGS

[0008] Specific embodiments of the disclosed technology will now be described in detail with reference to the accompanying figures. Like elements in the various figures are denoted by like reference numerals for consistency.

[0009] FIG. 1 illustrates a seismic survey in accordance with one or more embodiments.

[0010] FIG. 2A displays regularly-and densely-sampled input seismic data in accordance with one or more embodiments.

[0011] FIG. 2B displays an associated true input spectrum in accordance with one or more embodiments.

[0012] FIG. 3A displays irregularly-and sparsely-sampled input seismic data in accordance with one or more embodiments.

[0013] FIG. 3B displays an associated aliased and leaky input spectrum in accordance with one or more embodiments.

[0014] FIGs. 4A and 4B display input seismic data in accordance with one or more embodiments.

[0015] FIG. 5A displays an input spectrum in accordance with one or more embodiments.

[0016] FIG. 5B displays a true input spectrum in accordance with one or more embodiments.

[0017] FIG. 6A displays weight data in accordance with one or more embodiments.

[0018] FIG. 6B displays a weighted input spectrum in accordance with one or more embodiments.

[0019] FIGs. 7A and 7B display output seismic data in accordance with one or more embodiments.

[0020] FIGs. 8A and 8B display true seismic data in accordance with one or more embodiments.

[0021] FIGs. 9A and 9B display comparison output seismic data in accordance with one or more embodiments.

[0022] FIG. 10 displays error per iteration in accordance with one or more embodiments.

[0023] FIGs. 11 and 12 describe methods in accordance with one or more embodiments.

[0024] FIG. 13 illustrates a computer system in accordance with one or more embodiments.

[0025] FIG. 14 illustrates a drilling system in accordance with one or more embodiments.

[0026] FIG. 15 illustrates a systems flowchart in accordance with one or more embodiments.DETAILED DESCRIPTION

[0027] In the following detailed description of embodiments of the disclosure, numerous specific details are set forth in order to provide a more thorough understanding of the disclosure. However, it will be apparent to one of ordinary skill in the art that the disclosure may be practiced without these specific details. In other instances, well-known features have not been described in detail to avoid unnecessarily complicating the description.

[0028] Throughout the application, ordinal numbers (e g , first, second, third, etc ) may be used as an adjective for an element (i e , any noun in the application) . The use of ordinal numbers is not to imply or create any particular ordering of the elements nor to limit any element to being only a single element unless expressly disclosed, such as using the terms “before, ” “after, ” “single, ” and other such terminology. Rather, the use of ordinal numbers is to distinguish between the elements. By way of an example, a first element is distinct from a second element, and the first element may encompass more than one element and succeed (or precede) the second element in an ordering of elements.

[0029] It is to be understood that the singular forms “a, ” “an, ” and “the” include plural referents unless the context clearly dictates otherwise. Thus, for example, reference to an “element” includes reference to one or more of such elements.

[0030] Terms such as “approximately, ” “substantially, ” etc , mean that the recited characteristic, parameter, or value need not be achieved exactly, but that deviations or variations, including for example, tolerances, measurement error, measurement accuracy  limitations and other factors known to those of skill in the art, may occur in amounts that do not preclude the effect the characteristic was intended to provide.

[0031] It is to be understood that one or more of the steps shown in the flowcharts may be omitted, repeated, and / or performed in a different order than the order shown. Accordingly, the scope disclosed herein should not be considered limited to the specific arrangement of steps shown in the flowcharts.

[0032] Although multiple dependent claims are not introduced, it would be apparent to one of ordinary skill that the subject matter of the dependent claims of one or more embodiments may be combined with other dependent claims.

[0033] In the following description of FIGs. 1-15, any component described regarding a figure, in various embodiments disclosed herein, may be equivalent to one or more like-named components described regarding any other figure. For brevity, descriptions of these components will not be repeated regarding each figure. Thus, each and every embodiment of the components of each figure is incorporated by reference and assumed to be optionally present within every other figure having one or more like-named components. Additionally, in accordance with various embodiments disclosed herein, any description of the components of a figure is to be interpreted as an optional embodiment which may be implemented in addition to, in conjunction with, or in place of the embodiments described regarding a corresponding like-named component in any other figure.

[0034] Systems and methods are disclosed to attenuate aliasing and spectral leakage (collectively “spectral artifacts” ) in high-dimensional seismic data. Hereinafter, the high-dimensional seismic data is at least a function of four dimensions. As such, the high-dimensional seismic data may be referred to as being at least four dimensional though a person of ordinary skill in the art will appreciate that the seismic data includes an amplitude that is a function of the at least four dimensions. Each dimension may be any type of dimension some of which are described below.

[0035] The disclosed methods may be an improvement over existing methods, such as transform-based data regularization methods, convolution-based methods, wavefield continuation operator methods, minimum weighted norm inversion methods, deep learning  methods, and anti-leakage Fourier transform (ALFT) methods, for at least one of the following reasons. One reason is that the disclosed methods may be computationally faster while maintaining or increasing quality relative to existing methods due to the removal and replacement of amplitudes within a collection of neighboring elements rather than removing and replacing an amplitude of one element one at a time. Another reason is that the disclosed methods may lead to improved results relative to existing methods because high-dimensional seismic data (i e , additional data) is used. A third reason is that the disclosed methods may be designed to robustly identify true spectral events unlike existing methods that may erroneously identify aliased events (i e , aliasing) as true spectral events.

[0036] To clearly explain the concepts associated with the disclosed methods, low-dimensional seismic data may be displayed and / or described herein. If two-dimensional seismic data is displayed and / or described, a pixel represents each two-dimensional position within the two-dimensional seismic data. An amplitude is associated with each pixel. Hereinafter, the term “amplitude” may generically denote any type of amplitude. Examples of an amplitude include, but are not limited to, an amplitude, input amplitude, weight amplitude, weighted input amplitude, and output amplitude. If three-dimensional seismic data is displayed and / or described, a voxel represents each three-dimensional position within the three-dimensional seismic data. An amplitude is associated with each voxel. If at least four-dimensional seismic data is described, at least a toxel represents each at least four-dimensional position within the at least four-dimensional seismic data. An amplitude is associated with each at least toxel. Due to the dimensional limitations of visual displays, a toxel may not be displayable. As such, an at least toxel is a dimensional extension of a voxel and a voxel is a dimensional extension of a pixel. Hereinafter, the term “element” generically denotes a pixel, voxel, or toxel. Further, hereinafter, the term “spectral element” generically denotes a spectral pixel, spectral voxel, or spectral toxel.

[0037] To collect high-dimensional seismic data, a seismic survey may be performed over the subterranean region of interest. FIG. 1 illustrates a seismic survey 100 of a subterranean region of interest 105 in accordance with one or more embodiments. While FIG. 1 specifically illustrates a surface seismic survey, a person of ordinary skill in the art will  appreciate that other seismic surveys, such as vertical seismic profile (VSP) surveys, may be alternatively or additionally used without departing from the scope of the disclosure.

[0038] Returning to FIG. 1, the subterranean region of interest 105 is defined based on a coordinate system of three spatial dimensions 110, denoted x, y, and z in FIG. 1. The subterranean region of interest 105 may be made up of layers of rock 115 separated by geological boundaries 120 or other geological structures. The subterranean region of interest 105 may further include a hydrocarbon reservoir 125. The hydrocarbon reservoir 125 may be rock 115 filled with fluid such as oil, gas, water, brine, and / or a combination thereof.

[0039] The seismic survey 100 may be performed using a seismic acquisition system 130. The seismic acquisition system 130 may include a seismic source 135 and seismic receivers 140 positioned on or near the surface of the earth 145.

[0040] The seismic survey 100 may initially rely on the seismic source 135 configured to generate radiated seismic waves 150 (i e , emitted energy, wavefield) . The type of seismic source 135 may depend on the environment in which it is used. For example, on land, the seismic source 135 may be a vibroseis truck or explosive charge. In water, the seismic source 135 may be an airgun. The radiated seismic waves 150 may radiate along and into the subterranean region of interest 105. The radiated seismic waves 150 may return to the surface of the earth 145 as refracted seismic waves (not shown) or may be reflected by geological boundaries 120 or other geological structures and return to the surface of the earth 145 as reflected seismic waves 155. The radiated seismic waves 150 may also propagate along the surface of the earth 145 as Rayleigh waves or Love waves, collectively known as “ground roll” 160. Vibrations associated with ground roll 160 do not penetrate far beneath the surface of the earth 145 and, hence, are not influenced by, nor contain information about, portions of the subterranean region of interest 105 where hydrocarbon reservoirs 125 typically reside. The seismic receivers 140 located on or near the surface of the earth 145 are configured to detect radiated seismic waves 150, reflected seismic waves 155, refracted seismic waves, and ground roll 160.

[0041] Denoting the two-dimensional position of the seismic source 135 as (xs, ys) and the two-dimensional position of each seismic receiver 140 as (xr, yr) , respectively, the seismic trace recorded by each seismic receiver 140 may then be denoted S (xs, ys, xr, yr, t) , where t denotes recording time (i e , the time elapsed after the activation of the seismic source 135) . The collection of all seismic traces acquired during the seismic survey may be described as high-dimensional seismic data. As such, the high-dimensional seismic data may be initially collected in five dimensions.

[0042] The high-dimensional seismic data may be processed to transform the high-dimensional seismic data from one domain to another domain. In some embodiments, one or more transforms and / or seismic processing steps may be performed to transform the initially-collected five-dimensional seismic data to other high-dimensional seismic data. For example, initially-collected five-dimensional seismic data may be transformed to four-dimensional seismic data where the four dimensions include three spatial dimensions that correspond to the three spatial dimensions 110 of the subterranean region of interest 105 that the high-dimensional seismic data aims to characterize and a time dimension t.

[0043] However, in any domain, the high-dimensional seismic data may be irregularly sampled. In the context of this disclosure, the term “irregular” refers to the idea that the seismic data is not collected over equal intervals. Instead, the high-dimensional seismic data is collected over unequal intervals where the high-dimensional seismic data may be undersampled in some areas and oversampled in other areas. The high-dimensional seismic data may be irregularly sampled due to the cost of obtaining regularly-sampled seismic data, complex topography, cable feathering, and the removal of unusable seismic traces. However, use of irregularly-sampled high-dimensional seismic data may lead to an inadequate characterization of the subterranean region of interest 105. Regularization may be performed to determine regularly-sampled high-dimensional seismic data from collected irregularly-sampled high-dimensional seismic data. Further, in some embodiments, the irregularly-sampled high-dimensional seismic data may also be sparsely-sampled.

[0044] Regularization may be performed by suppressing or attenuating the aliasing and spectral leakage within the spectrum associated with the irregularly-sampled high-dimensional seismic data. This type of regularization may be referred to as “anti-aliasing and anti-leakage regularization” or “anti-aliasing and anti-leakage Fourier reconstruction. ” This type of regularization may also be considered a method of interpolation. FIGs. 2A, 2B, 3A, and 3B illustrate aliasing and spectral leakage in the context of two-dimensional seismic data for clarity.

[0045] FIG. 2A displays regularly-and densely-sampled input seismic data 200 in accordance with one or more embodiments. Hereinafter, the terms “dense” and “sparse” are considered antonymous. FIG. 2A displays the amplitude of the regularly-and densely-sampled input seismic data 200 as a function of two dimensions, a time dimension t and a space dimension x. As such, each element 205, specifically each pixel, within the regularly-and densely-sampled input seismic data 200 is associated with a time t, spatial position x, and amplitude. The regularly-and densely-sampled input seismic data 200 clearly displays three linear and crossing true seismic events 210 that may be manifestations of geological structures, such as geological boundaries 120, that physically exist within the subterranean region of interest 105.

[0046] FIG. 2B displays an associated true input spectrum 215 associated with the regularly-and densely-sampled input seismic data 200 of FIG. 2A in accordance with one or more embodiments. The time dimension t has been transformed to a frequency dimension ω. The space dimension x has been transformed to a spatial frequency dimension kx (hereinafter also “wavenumber dimension” ) . Each spectral element 220, specifically each spectral pixel, within the associated true input spectrum 215 is associated with a temporal frequency ω, wavenumber kx, and input amplitude.

[0047] The associated true input spectrum 215 displays the Fourier components of the regularly-and densely-sampled input seismic data 200 as the input amplitudes. Fourier components with high input amplitudes may be referred to as “spectral events. ” Relative to FIG. 2B, the associated true input spectrum 215 displays only true spectral events 225.

[0048] Turning to FIG. 3A, FIG. 3A displays irregularly-and sparsely-sampled input seismic data 300 in accordance with one or more embodiments. Here, 7 / 8ths of the seismic traces within the regularly-and densely-sampled input seismic data 200 displayed in FIG. 2A have been removed in regular intervals and then 30%of those seismic traces have been randomly removed to determine the irregularly-and sparsely-sampled input seismic data 300 displayed in FIG. 3A.

[0049] The associated aliased and leaky input spectrum 305 is displayed in FIG. 3B in accordance with one or more embodiments. The associated aliased and leaky input spectrum 305 continues to display the true spectral events 225 and now also displays aliased events (hereinafter “aliasing” 310) and spectral leakage 315. The aliasing 310 displays as periodic spectral events on either side of the true spectral events 225. Spectral leakage 315 displays as spectral noise between the true spectral events 225 and / or the periodic aliasing 310.

[0050] However, in the context of this disclosure, the input seismic data is at least four dimensional unlike the two-dimensional input seismic data 200 and 300 displayed in FIGs. 2A and 3A. As such, any of the input seismic data 200 and 300 displayed in FIGs. 2A and 3A may be considered a “slice” within the at least four-dimensional input seismic data. Each slice may be any dimensionality that is smaller than the dimensionality of the at least four-dimensional input seismic data. Hereinafter, the at least four-dimensional input seismic data is simply referred to as “input seismic data. ”

[0051] FIGs. 4A and 4B display input seismic data 400 in accordance with one or more embodiments. Because the input seismic data 400 is four dimensional here, FIG. 4A displays a first view of the input seismic data 400 and FIG. 4B displays a second view of the input seismic data 400. Here, the four dimensions include three spatial dimensions x, y, and z and one time dimension t. As such, each element 205, specifically each toxel, within the input seismic data 400 is associated with three spatial positions x, y, and z, time t, and an amplitude. However, FIGs. 4A and 4B can only display an amplitude at each voxel due to the dimensional limitations of visual displays. The input seismic data 400 of FIGs. 4A and 4B are irregularly and sparsely sampled.

[0052] One or more transforms may be applied to the input seismic data 400 to determine an at least four-dimensional input spectrum (hereinafter also simply “input spectrum” ) . The input seismic data 400 may be denoted d (xn, t) where xn are three-dimensional positions [x0, x1, ..., xN-1] and t is time. In some embodiments, a fast Fourier transform (FFT) may be applied to the input seismic data d (xn, t) 400 along the time dimension t to determine an intermediate input spectrum d (xn, ω) . In some embodiments, a nonuniform discrete Fourier transform (NDFT) may then be applied to the intermediate input spectrum d (xn, ω) to determine the input spectrum D (mΔk, ω) where:

[0053] In Equation (1) , Δk=2π /  (xN-1-x0) , m is a wavenumber index where 0≤m≤N-1, and

[0054] FIG. 5A displays an input spectrum 500 in accordance with one or more embodiments. Each spectral element 220, specifically each spectral toxel, within the input spectrum 500 is associated with three wavenumbers kx, ky, and ky, a temporal frequency ω, and an input amplitude. However, FIG. 5A can only display an input amplitude at each spectral voxel due to the dimensional limitations of visual displays. The input spectrum 500 may include true spectral events 225 as well as aliasing 310 and / or spectral leakage 315. For comparison, FIG. 5B displays a true input spectrum 505 in accordance with one or more embodiments.

[0055] Anti-aliasing and anti-leakage regularization may attenuate aliasing 310 and / or spectral leakage 315 within the input spectrum 500. In turn, an output spectrum with attenuated aliasing 310 and / or spectral leakage 315 may be inversely transformed to determine regularly-and densely-sampled output seismic data.

[0056] To attenuate aliasing 310 and / or spectral leakage 315 within the input spectrum 500, a weight function may be used and take the form: W (k′, ω′) =ΣΩ|D (k′ω / ω′, ω) |,      Equation (2)

[0057] where W is weight data, D is the input spectrum 500, ω′is the temporal frequency associated with each spectral element 220 along each radial line Ω, and k′is the  wavenumber vector associated with each spectral element 220 along each radial line Ω where k′= (k′x, k′y, k′z) . However, a person of ordinary skill in the art will appreciate that D, ω′, and k′may represent other dimensions and be other dimensionalities without departing from the scope of the disclosure.

[0058] As illustrated in FIGs. 3B and 5A, a first endpoint 320 of each radial line Ω 325 is positioned at the origin 330 of the associated aliased and leaky input spectrum 305 and input spectrum 500. A second endpoint 335 of each radial line 325 is positioned at the spectral element 220 that each weight within the weight data W is ultimately assigned to. A summation (i e , stacking) is performed along each radial line Ω 325 at unique orientations and lengths to generate the weight data W as Equation (2) mathematically describes. For example, FIG. 3B displays two radial lines 325. The first endpoint 320 of each radial line Ω 325 is positioned at the origin 330 of the associated aliased and leaky input spectrum 305. The orientation of each radial line 325 is defined by an angle α 340 relative to a datum 345 and length l 350. As such, in these embodiments, each radial line Ω 325 is a function of the angle α 340 and the length l 350 (i e , Ω (α, l) ) . The angle α 340 and the length l 350 are selected from a set of angles and set of lengths, respectively. Stacking is performed along each radial line Ω 325 where each radial line Ω 325 is associated with each angle α 340 in the set of angles and each angle α 340 is ultimately associated with every length l 350 in the set of lengths. For example, if the set of angles includes three angles and the set of lengths includes four lengths, 12 unique radial lines Ω 325 exist and Equation (2) is determined 12 times (i e , 3·4 radial lines and times) . However, a person of ordinary skill in the art will appreciate that, in practice, the set of angles and set of lengths may include tens to hundreds of angles and lengths, respectively. Further, a person of ordinary skill in the art will appreciate that each radial line Ω 325 may be further defined by more than one angle illustrated using α and β in FIG. 5A. The number of angles may depend on the dimensionality of the input spectrum 500 and / or the dimensionality of slices, if used.

[0059] FIG. 6A illustrates weight data W 600 in accordance with one or more embodiments. Here the weight data W 600 is generated by applying the weight function  of Equation (2) to the associated aliased and leaky input spectrum 305 displayed in FIG. 3B. A person of ordinary skill in the art will appreciate that, in practice, the weight data 600 displayed in FIG. 6A may take a higher dimensional form or may be considered a weight data slice without departing from the scope of the disclosure. As such, in some embodiments, at least four-dimensional weight data 600 may be immediately generated. In other embodiments, the at least four-dimensional weight data 600 may be generated on a slice-by-slice basis to determine weight data slices. As such, in some embodiments, the weight data 600 displayed in FIG. 6A may be considered a weight data slice.

[0060] The associated aliased and leaky input spectrum 305 displayed in FIG. 3B may be weighted by the weight data 600 displayed in FIG. 6A. In some embodiments, the weight data 600 and associated aliased and leaky input spectrum 305 may be multiplied together to generate a weighted input spectrum. FIG. 6B displays a weighted input spectrum 605 in accordance with one or more embodiments. In some embodiments, an at least four-dimensional weighted input spectrum 605 may be immediately generated. In other embodiments, the at least four-dimensional weighted input spectrum 605 may be generated on a slice-by-slice basis to determine weighted input spectrum slices. As such, in some embodiments, the weighted input spectrum 605 displayed in FIG. 6B may be considered a weighted input spectrum slice.

[0061] A spectral element 610 associated with the maximum weighted input amplitude (hereinafter “maximum amplitude” 615) within the weighted input spectrum 605 may then be identified. For example, FIG. 6B illustrates the spectral element 610 associated with the maximum amplitude 615. Neighboring spectral elements 620 that include the spectral element 610 associated with the maximum amplitude 615 may then be identified. The neighboring spectral elements 620 may take any dimensionality less than or equal to the dimensionality of the input seismic data 400. Further, the neighboring spectral elements 620 may take any shape, such as, but not limited to, a square, cube, or hypercube.

[0062] The input spectrum 500 and an output spectrum may then be updated based, at least in part, on the input amplitude associated with each of the neighboring spectral elements 620 within the input spectrum 500. To update the output spectrum, the input amplitude  associated with each of the neighboring spectral elements 620 within the input spectrum 500 replaces the output amplitude associated with each of the neighboring spectral elements 620 within the output spectrum. To update the input spectrum 500, the input amplitude associated with each of the neighboring spectral elements 620 within the input spectrum 500 is removed where:

[0063] Updated weight data 600 may then be generated by applying the weight function of Equation (2) to the updated input spectrum 500. An updated weighted input spectrum 605 may then be generated by weighting the updated input spectrum 500 using the updated weight data 600. The updated spectral element 610 associated with the updated maximum amplitude 615 is then identified within the updated weighted input spectrum 605. The updated neighboring spectral elements 620 are then identified that includes the updated spectral element 610 associated with the updated maximum amplitude 615. The updated output spectrum is further updated by replacing the output amplitude associated with each of the updated neighboring spectral elements 620 within the updated output spectrum with the input amplitude associated with each of the updated neighboring spectral elements 620 within the updated input spectrum 500. The updated input spectrum 500 is further updated by removing the input amplitude associated with each of the updated neighboring spectral elements 620 within the updated input spectrum 500.

[0064] This iterative or recursive updating process continues until a convergence criterion is satisfied. The convergence criterion may be based on a pre-defined residual threshold or a pre-defined number of iterations, in some embodiments.

[0065] Once the convergence criterion is satisfied, the output spectrum may include true spectral events 225 and attenuated aliasing 310 and / or spectral leakage 315. One or more inverse transforms may be applied to the output spectrum to determine output seismic data. In some embodiments, an inverse FFT may be applied to the output spectrum along one dimension, such as a temporal frequency dimension ω, to determine an intermediate output spectrum. In some embodiments, an inverse NDFT may then be applied to the intermediate output spectrum to determine the output seismic data. FIGs. 7A and 7B display output  seismic data 700 after 80 iterations in accordance with one or more embodiments. Because the output seismic data 700 is four dimensional here, FIG. 7A displays a first view of the output seismic data 700 and FIG. 7B displays a second view of the output seismic data 700. Comparison of the output seismic data 700 in FIGs. 7A and 7B to the input seismic data 400 in FIGs. 4A and 4B illustrates that the output seismic data 700 may now be regularly and densely sampled. As such, the output seismic data 700 includes clear true seismic events 210 and attenuated aliasing 310 and / or spectral leakage 315.

[0066] FIGs. 8A and 8B display true seismic data 800 in accordance with one or more embodiments. Comparison of the output seismic data 700 in FIGs. 7A and 7B to the true seismic data 800 in FIGs. 8A and 8B illustrates that the disclosed anti-aliasing and anti-leakage regularization methods may adequately attenuate aliasing 310 and / or spectral leakage 315 to determine regularly-and densely-sampled seismic data.

[0067] FIGs. 9A and 9B display comparison output seismic data 900 in accordance with one or more embodiments. The comparison output seismic data 900 is determined using the disclosed anti-aliasing and anti-leakage regularization methods but where replacement and removal of only the input amplitude associated with the spectral element 610 with the maximum amplitude 615 is used (hereinafter referred to as the “alternative method” ) rather than the input amplitude associated with each of the neighboring spectral elements 620. Note that the comparison output seismic data 900 in FIGs. 9A and 9B are determined following 80 iterations as the output seismic data 700 in FIGs. 7A and 7B are. Comparison of the comparison output seismic data 900 in FIGs. 9A and 9B to the output seismic data 700 in FIGs. 7A and 7B and the true seismic data 800 in FIGs. 8A and 8B qualitatively illustrate that the alternative method may not adequately attenuate spectral leakage 315.

[0068] Further, FIG. 10 displays the error per iteration associated with the disclosed anti-aliasing and anti-leakage regularization methods versus the alternative method. The error is the L2 norm of the difference between each method and the true seismic data 800. FIG. 10 quantitatively illustrates that the disclosed anti-aliasing and anti-leakage regularization methods may converge quicker (i e , be computationally faster) than the alternative method. For example, the disclosed anti-aliasing and anti-leakage regularization methods  completes 40 iterations in 4.2 seconds while the alternative method completes 100 iterations in 6.1 seconds to reach approximately the same error.

[0069] FIGs. 11 and 12 describe the disclosed anti-aliasing and anti-leakage regularization methods in accordance with one or more embodiments. In step 1100 in FIG. 11, an output spectrum is initialized. In some embodiments, the output spectrum may be initialized as a null array. However, a person of ordinary skill in the art will appreciate that the output spectrum may be initialized to include any values within an array without departing from the scope of the disclosure.

[0070] In step 1105, input seismic data 400 is received for a subterranean region of interest 105. In some embodiments, the input seismic data 400 may be received from a seismic acquisition system 130. The input seismic data 400 is at least four dimensions. Each dimension may be any type of dimension. In some embodiments, the dimensions of the input seismic data 400 may include the two-dimensional spatial dimensions, xs and ys, of the seismic source 135, the two-dimensional spatial dimensions, xr and yr, of each of the seismic receivers 140, and a time dimension t. In other embodiments, the dimensions of the input seismic data 400 may include three spatial dimensions, x, y, and z, that correspond to the three spatial dimensions 110 of the subterranean region of interest 105 that the input seismic data 400 aims to characterize and a time dimension t. As such, the input seismic data 400 includes at least toxels. Each at least toxel is associated with an amplitude. For reference, FIGs. 2A and 3A display two-dimensional input seismic data 200 and 300 and FIGs. 4A and 4B display four-dimensional input seismic data 400 using two three-dimensional views. In some embodiments, the input seismic data 400 may be irregularly sampled and / or sparsely sampled.

[0071] In step 1110, an input spectrum 500 is determined from the input seismic data 400 using one or more transforms. In some embodiments, the one or more transforms may include an FFT and NDFT used in serial as previously described. In some embodiments, the input spectrum 500 may be immediately determined. In other embodiments, the input spectrum 500 may be determined on a slice-by-slice basis to determine input spectrum slices. The input spectrum 500 includes spectral elements 220. The spectral elements 220  may be spectral pixels, spectral voxels, or spectral toxels. The spectral elements 220 may be spectral pixels or spectral voxels if input spectrum slices are two or three dimensional, respectively. Each spectral element 220 is associated with an input amplitude. For reference, FIGs. 2B and 3B display two-dimensional input spectra 215 and 305 and FIGs. 5A and 5B display a four-dimensional input spectrum 500 using two three-dimensional views. In some embodiments, the input spectrum 500 may include true spectral events 225, aliasing 310, and / or spectral leakage 315.

[0072] In step 1115, an updated output spectrum is determined iteratively or recursively until a convergence criterion is satisfied. The details of step 1115 are described relative to FIG. 12. In some embodiments, the updated output spectrum may include true spectral events 225 and attenuated aliasing 310 and / or spectral leakage 315.

[0073] In step 1120, output seismic data 700 is determined from the updated output spectrum using one or more inverse transforms. In some embodiments, the one or more inverse transforms may include an inverse FFT and inverse NDFT used in serial. In some embodiments, the output seismic data 700 may be immediately determined. In other embodiments, the output seismic data 700 may be determined on a slice-by-slice basis to determine output seismic data slices.

[0074] In step 1125, a seismic image is generated of the subterranean region of interest 105 based, at least in part, on the output seismic data 700. Generating the seismic image may include other seismic processing steps known to a person of ordinary skill in the art. One or more seismic processing steps may be designed to correct for anomalies such as near-surface effects, noise, seismic survey geometry irregularities, migration artifacts, etc.

[0075] In step 1130, a drilling target of a hydrocarbon reservoir 125 within the subterranean region of interest 105 is identified using, at least in part, the seismic image.

[0076] Turning to FIG. 12, steps 1115a-e are performed iteratively or recursively until a convergence criterion is satisfied. In some embodiments, steps 1115a-e may be performed over the entire input spectrum 500. In other embodiments, steps 1115a-e may be performed on a slice-by-slice basis for each input spectrum slice within the input spectrum 500. In these embodiments, each input spectrum slice may be any dimensionality that is less than  the dimensionality of the input seismic data 400. Further, in these embodiments, the convergence criterion may be satisfied for each input spectrum slice individually.

[0077] In step 1115a, weight data 600 is generated from the input spectrum 500 using a weight function. In some embodiments, the weight function may take the form of Equation (2) as previously described. In some embodiments, the weight data 600 may be the same dimensionality as the input spectrum 500. In other embodiments, the weight data 600 may be a weight data slice generated from each input spectrum slice within the input spectrum 500. For reference, FIG. 6A displays two-dimensional weight data 600, which may be considered a weight data slice.

[0078] In step 1115b, a weighted input spectrum 605 is generated by weighting the input spectrum 500 using the weight data 600. In some embodiments, the weighted input spectrum 605 may be the same dimensionality as the input spectrum 500. In other embodiments, the weighted input spectrum 605 may be organized into weighted input spectrum slices of a lower dimensionality or generated by weighting each input spectrum slice using each weight data slice.

[0079] In step 1115c, neighboring spectral elements 620 are identified in the weighted input spectrum 605 or, in some embodiments, each weighted input spectrum slice. A spectral element 610 among the neighboring spectral elements 620 is associated with a maximum amplitude 615 within the weighted input spectrum 605 or each weighted input spectrum slice. The neighboring spectral elements 620 may take any dimensionality less than or equal to the dimensionality of the input seismic data 400. Further, the neighboring spectral elements 620 may take any shape, such as a square, rectangle, cube, cuboid, hypercube, hypercuboid, or any other shape that may be formed from neighboring spectral elements 620. If two or three dimensional weighted input spectrum slices are used, the neighboring spectral elements 620 may be neighboring spectral pixels or neighboring spectral voxels.

[0080] In step 1115d, the output spectrum is updated. To do so, the input amplitude associated with each of the neighboring spectral elements 620 within the input spectrum  500 replaces the output amplitude associated with each of the neighboring spectral elements 620 within the output spectrum or output spectrum slice.

[0081] In step 1115e, the input spectrum is updated. To do so, the input amplitude associated with each of the neighboring spectral elements 620 within the input spectrum 500 is removed from the input spectrum 500 or input spectrum slice.

[0082] Steps 1100, 1105, 1110, 1115 (including 1115a-e) , 1120, and 1125 may be performed using a seismic processing system. Step 1130 may be performed using a seismic interpretation workstation. The seismic processing system and seismic interpretation workstation may each be a specifically-configured computer system.

[0083] FIG. 13 illustrates a generic computer system 1300 in accordance with one or more embodiments. As mentioned, the computer system 1300 may be specifically configured for seismic processing and denoted a “seismic processing system. ” Alternatively, the computer system 1300 may be specifically configured for seismic interpretation and denoted a “seismic interpretation workstation. ” The seismic processing system, seismic interpretation workstation, or generic computer system 1300 may store and be used by the wellbore planning system 1340. While the generic term computer system 1300 may be used to describe each of the parts of a computer system 1300 in the following paragraphs, the terms seismic processing system or seismic interpretation workstation may replace the term computer system 1300 without departing from the scope of the disclosure.

[0084] The computer system 1300 is intended to depict any computing device such as a server, desktop computer, laptop / notebook computer, wireless data port, smart phone, personal data assistant (PDA) , tablet computing device, one or more processors within these devices, or any other suitable processing device, including both physical or virtual instances (or both) of the computing device. Additionally, the computer system 1300 may include an input device, such as a keypad, keyboard, touch screen, or other device that can accept user information, and an output device that displays information, including digital data, visual or audio information (or a combination of both) , or a graphical user interface (GUI) . Specifically, a seismic interpretation workstation may include a robust graphics card for the detailed rendering of the seismic image such that the seismic image may be  displayed and manipulated in a virtual reality system using 3D goggles, a mouse, or a wand to identify a drilling target of a hydrocarbon reservoir 125 within the subterranean region of interest 105.

[0085] The computer system 1300 can serve in a role as a client, network component, server, database, or any other component (or a combination of roles) of a computer system 1300 as required for seismic processing and seismic interpretation. The illustrated computer system 1300 is communicably coupled with a network 1305. For example, a seismic processing system and a seismic interpretation workstation may be communicably coupled using a network 1305. In some implementations, one or more components of each computer system 1300 may be configured to operate within environments, including cloud-computing-based, local, global, or other environment (or a combination of environments) .

[0086] At a high level, the computer system 1300 is an electronic computing device operable to receive, transmit, process, store, and / or manage data and information associated with the disclosed methods. According to some implementations, the computer system 1300 may also include or be communicably coupled with an application server, e-mail server, web server, caching server, streaming data server, business intelligence (BI) server, or other server (or a combination of servers) .

[0087] Because seismic processing and seismic interpretation may not be sequential, the computer system 1300 can receive requests over network 1305 from other computer systems 1300 or another client application and respond to the received requests by processing the requests appropriately. In addition, requests may also be sent to the computer system 1300 from internal users (for example, from a command console or by other appropriate access method) , external or third-parties, other automated applications, as well as any other appropriate entities, individuals, systems, or computer systems 1300.

[0088] Each of the components of the computer system 1300 can communicate using a system bus 1310. In some implementations, any or all of the components of each computer system 1300, both hardware or software (or a combination of hardware and software) , may interface with each other or the interface 1315 (or a combination of both) over the system bus 1310 using an application programming interface (API) 1320 or a service layer 1325  (or a combination of the API 1320 and service layer 1325. The API 1320 may include specifications for routines, data structures, and object classes. The API 1320 may be either computer-language independent or dependent and refer to a complete interface, a single function, or even a set of APIs. The service layer 1325 provides software services to each computer system 1300 or other components (whether or not illustrated) that are communicably coupled to each computer system 1300. The functionality of each computer system 1300 may be accessible for all service consumers using this service layer 1325. Software services, such as those provided by the service layer 1325, provide reusable, defined business functionalities through a defined interface. For example, the interface may be software written in JAVA, C++, or other suitable language providing data in extensible markup language (XML) format or another suitable format. While illustrated as an integrated component of each computer system 1300, alternative implementations may illustrate the API 1320 or the service layer 1325 as stand-alone components in relation to other components of each computer system 1300 or other components (whether or not illustrated) that are communicably coupled to each computer system 1300. Moreover, any or all parts of the API 1320 or the service layer 1325 may be implemented as child or sub-modules of another software module, enterprise application, or hardware module without departing from the scope of this disclosure.

[0089] The computer system 1300 includes an interface 1315. Although illustrated as a single interface 1315 in FIG. 13, two or more interfaces 1315 may be used according to particular needs, desires, or particular implementations of each computer system 1300. The interface 1315 is used by each computer system 1300 for communicating with other systems in a distributed environment that are connected to the network 1305. Generally, the interface 1315 includes logic encoded in software or hardware (or a combination of software and hardware) and operable to communicate with the network 1305. More specifically, the interface 1315 may include software supporting one or more communication protocols associated with communications such that the network 1305 or interface’s hardware is operable to communicate physical signals within and outside of the illustrated computer system 1300.

[0090] The computer system 1300 includes at least one computer processor 1330. Generally, a computer processor 1330 executes any instructions, algorithms, methods, functions, processes, flows, and procedures as described above. A computer processor 1330 may be a central processing unit (CPU) and / or a graphics processing unit (GPU) . The input seismic data 400 may be tens to hundreds of terabytes or even petabytes in size. To efficiently process the input seismic data 400, in series and / or in parallel, to generate a seismic image, a seismic processing system may consist of an array of CPUs with one or more subarrays of GPUs attached to each CPU. Further, tape readers or high-capacity hard-drives may be connected to the CPUs using wide-band system buses 1310.

[0091] The computer system 1300 also includes a memory 1335 that stores data and software for the computer system 1300 or other components (or a combination of both) that can be connected to the network 1305. For example, the memory 1335 may store the input seismic data 400 as well as the wellbore planning system 1340 in the form of software as shown in FIG. 13. In some embodiments, the wellbore planning system 1340 may be configured to plan the wellbore path. Although illustrated as a single memory 1335 in FIG. 13, two or more memories may be used according to particular needs, desires, or particular implementations of the computer system 1300 and the described functionality. While memory 1335 is illustrated as an integral component of each computer system 1300, in alternative implementations, memory 1335 can be external to each computer system 1300.

[0092] The application 1345 is an algorithmic software engine providing functionality according to particular needs, desires, or particular implementations of the computer system 1300, particularly with respect to functionality described in this disclosure. For example, application 1345 can serve as one or more components, modules, applications, etc. Further, although illustrated as a single application 1345, the application 1345 may be implemented as multiple applications 1345 on each computer system 1300. In addition, although illustrated as integral to each computer system 1300, in alternative implementations, the application 1345 can be external to each computer system 1300.

[0093] There may be any number of computer systems 1300, such as computer clusters, associated with, or external to, a seismic processing system and a seismic interpretation  workstation, where each computer system 1300 communicates over network 1305. Further, the term “client, ” “user, ” and other appropriate terminology may be used interchangeably as appropriate without departing from the scope of this disclosure. Moreover, this disclosure contemplates that many users may use the computer system 1300, or that one user may use multiple computer systems 1300.

[0094] Returning to the discussion of the wellbore plan, the wellbore plan may be additionally informed by the best available information at the time of planning. This may include models encapsulating stress conditions of the subterranean region of interest 105, the trajectory of any existing wellbores (which may be desirable to avoid) , and the existence of other drilling hazards, such as shallow gas pockets, over-pressure zones, and active fault planes.

[0095] The wellbore path may include a starting surface location of the wellbore, or a subsurface location within an existing wellbore, from which the wellbore may be drilled. The wellbore path may further include a terminal location that may intersect with the drilling target of the hydrocarbon reservoir 125. The wellbore path may further still include wellbore geometry information such as wellbore diameter and inclination angle and when each of these change along the depth of the wellbore. If casing is used, the wellbore plan may include casing type or casing depths. Furthermore, the wellbore plan may consider other engineering constraints such as the maximum wellbore curvature ( “dog-log” ) that a drillstring of a drilling system may tolerate and the maximum torque and drag values that the drilling system may tolerate. The wellbore plan may further define associated drilling parameters, such as the planned depths at which casing will be inserted to support the wellbore to prevent formation fluids entering the wellbore and the drilling mud weights (densities) and types that may be used during drilling of the wellbore.

[0096] The wellbore plan may then be transferred to a drilling system 1400 such that the wellbore 1405 guided by the wellbore path 1410 may be drilled as illustrated in FIG. 14 in accordance with one or more embodiments. Although the drilling system 1400 shown in FIG. 14 is used to drill the wellbore 1405 on land, the drilling system 1400 may be a marine wellbore drilling system. Further, although the drilling system 1400 shown in FIG. 14 is  used to drill a new wellbore 1405, the wellbore 1405 being drilled may be a sidetrack wellbore or an offset wellbore. As such, the example of the drilling system 1400 shown in FIG. 14 is not meant to limit the present disclosure.

[0097] As shown in FIG. 14, the wellbore 1405 may be drilled using a drill rig that may be situated on a land drill site, an offshore platform, such as a jack-up rig, a semi-submersible, or a drill ship. The drill rig may be equipped with a hoisting system, such as a derrick 1415, which can raise or lower the drillstring 1420 and other tools required to drill the wellbore 1405. The drillstring 1420 may include one or more drill pipes connected to form conduit and a bottom hole assembly 1425 (BHA) disposed at the distal end of the drillstring 1420. The BHA 1425 may include a drill bit 1430 to cut into rock 115, including cap rock 115a. The BHA 1425 may further include measurement tools, such as a measurement-while-drilling (MWD) tool and logging-while-drilling (LWD) tool. MWD tools may include sensors and hardware to measure downhole drilling parameters, such as the azimuth and inclination of the drill bit 1430, the weight-on-bit, and the torque. The LWD measurements may include sensors, such as resistivity, gamma ray, and neutron density sensors, to characterize the rock 115 surrounding the wellbore 1405. Both MWD and LWD measurements may be transmitted to the surface of the earth 145 using any suitable telemetry system known in the art, such as a mud-pulse or by wired-drill pipe.

[0098] To start drilling, or “spudding in, ” the wellbore 1405, the hoisting system lowers the drillstring 1420 suspended from the derrick 1415 of the drill rig towards the planned surface location of the wellbore 1405. An engine, such as a diesel engine, may be used to supply power to the top drive 1435 to rotate the drillstring 1420 via the drive shaft 1440. The weight of the drillstring 1420 combined with the rotational motion enables the drill bit 1430 to bore the wellbore 1405.

[0099] The near-surface rock 115 of the subterranean region of interest 105 is typically made up of loose or soft sediment or rock, so large diameter casing 1445 (e g , “base pipe” or “conductor casing” ) is often put in place while drilling to stabilize and isolate the wellbore 1405. At the top of the base pipe is the wellhead, which serves to provide pressure control through a series of spools, valves, or adapters (not shown) . Once near-surface  drilling has begun, water or drill fluid may be used to force the base pipe into place using a pumping system until the wellhead is situated just above the surface of the earth 145.

[0100] Drilling may continue without any casing 1445 once deeper or more compact rock 115 is reached. While drilling, a drilling mud system 1450 may pump drilling mud from a mud tank on the surface of the earth 145 through the drill pipe. Drilling mud serves various purposes, including pressure equalization, removal of rock cuttings, and drill bit cooling and lubrication.

[0101] At planned depth intervals, drilling may be paused and the drillstring 1420 withdrawn from the wellbore 1405. Sections of casing 1445 may be connected, inserted, and cemented into the wellbore 1405. Casing string may be cemented in place by pumping cement and mud, separated by a “cementing plug, ” from the surface of the earth 145 through the drill pipe. The cementing plug and drilling mud force the cement through the drill pipe and into the annular space between the casing 1445 and the wall of the wellbore 1405. Once the cement cures, drilling may recommence. The drilling process is often performed in several stages. Therefore, the drilling and casing cycle may be repeated more than once, depending on the depth of the wellbore 1405 and the pressure on the walls of the wellbore 1405 from surrounding rock 115.

[0102] Due to the high pressures experienced by deep wellbores 1405, a blowout preventer (BOP) may be installed at the wellhead to protect the rig and environment from unplanned oil or gas releases. As the wellbore 1405 becomes deeper, both successively smaller drill bits 1430 and casing 1445 may be used. Drilling deviated or horizontal wellbores 1405 may require specialized drill bits 1430 or drill assemblies.

[0103] The drilling system 1400 may be disposed at and communicate with other systems in the wellbore environment. The drilling system 1400 may control at least a portion of a drilling operation by providing controls to various components of the drilling operation. In one or more embodiments, the system may receive data from one or more sensors arranged to measure controllable parameters of the drilling operation. As a non-limiting example, sensors may be arranged to measure weight-on-bit, drill rotational speed (RPM) , flow rate of the mud pumps (GPM) , and rate of penetration of the drilling operation (ROP) . Each  sensor may be positioned or configured to measure a desired physical stimulus. Drilling may be considered complete when a drilling target with the hydrocarbon reservoir 125 is reached or the presence of hydrocarbons is established.

[0104] Turning to FIG. 15, FIG. 15 illustrates a systems flowchart in accordance with one or more embodiments. A seismic acquisition system 130 may be configured to collect the input seismic data 400 as previously described relative to FIG. 1. The input seismic data 400 may be transferred to a seismic processing system 1300a via a network 1305 and stored on the seismic processing system 1300a as previously described relative to FIG. 13. The seismic processing system 1300a may be configured to perform steps 1100, 1105, 1110, 1115 (including 1115a-e) , 1120, and 1125 as described relative to FIGs. 11 and 12 to ultimately generate a seismic image. Prior, during, and / or after steps 1100, 1105, 1110, 1115 (including 1115a-e) , 1120, and 1125 are performed, the seismic processing system 1300a may apply one or more general seismic processing steps to any of the data and / or spectra. General seismic processing steps may include, but are not limited to, sorting, conditioning, stacking, quality controling, migrating, and velocity analyzing. Sorting may include organizing seismic traces by attribute to determine gathers within any seismic data. Conditioning may include filtering, deconvolving, noise suppressing, and signal enhancing any of the data and / or spectra to enhance quality. Stacking may include combining seismic traces of any of the data and / or spectra to enhance quality and increase signal-to-noise ratio. Quality controling may include identifying and correcting any random artifacts, noise, and / or errors in the data and / or spectra. Migrating may include shifting the manifestation of geological structures, such as geological boundaries 120, within any seismic data to locations associated with the true locations of the structures that physically exist within the subterranean region of interest 105. Velocity analyzing may include estimating a time-depth relationship for true seismic events 210 within any seismic data. Prior, during, and / or after steps 1100, 1105, 1110, 1115 (including 1115a-e) , 1120, and 1125 are performed, the data and / or spectra may be displayed to allow geoscientists to determine how to further process the data and / or spectra.

[0105] Returning to the seismic image, the seismic image may be transferred from the seismic processing system 1300a to a seismic interpretation workstation 1300b via the  network 1305 and stored on the seismic interpretation workstation 1300b as previously described relative to FIG. 13. The seismic interpretation workstation 1300b may be configured to perform step 1130 as described relative to FIG. 11 to identify a drilling target of a hydrocarbon reservoir 125 within the subterranean region of interest 105 using, at least in part, the seismic image. To do so, the seismic interpretation workstation 1300b may display a high-resolution version of the seismic image as a volume, cross sections, time slices, and / or attribute maps such that the hydrocarbon reservoir 125 may be identified and the drilling target within the hydrocarbon reservoir 125 identified.

[0106] The seismic image with identified drilling target may be transferred to a wellbore planning system 1340. In some embodiments, the wellbore planning system 1340 may be stored on a memory 1335 of the generic computer system 1300, seismic processing system 1300a, or seismic interpretation workstation 1300b as illustrated in FIG. 13. The wellbore planning system 1340 may be configured to plan a wellbore path 1410 that penetrates the identified drilling target of the hydrocarbon reservoir 125. The wellbore planning system 1340 may be further configured to design drilling operations associated with the drilling of the planned wellbore path 1410. The wellbore planning system 1340 may visually display the planned wellbore path 1410 with the seismic image and other data, such as geological models and well logs, such that drilling engineers can visualize and interact with the planned wellbore path 1410, seismic image, and other data.

[0107] The planned wellbore path 1410 may be transferred to and stored on a drilling system 1400 as described relative to FIG. 14. The drilling system 1400 may be configured to drill the wellbore 1405 guided by the planned wellbore path 1410. The wellbore 1405 may be a continuation of a previously partially-drilled wellbore 1405, sidetrack wellbore, or offset wellbore. Following drilling and completion of the wellbore 1405, the wellbore 1405 may be configured to produce hydrocarbons from the hydrocarbon reservoir 125 to the surface of the earth 145.

[0108] Although only a few example embodiments have been described in detail above, those skilled in the art will readily appreciate that many modifications are possible in the example embodiments without materially departing from this invention. Accordingly, all  such modifications are intended to be included within the scope of this disclosure as defined in the following claims.

Claims

1.A method comprising;using a seismic processing system;initializing an output spectrum,receiving, from a seismic acquisition system, input seismic data for a subterranean region of interest,wherein the input seismic data comprises a plurality of toxels,determining, using a transform, an input spectrum from the input seismic data,wherein the input spectrum comprises an input amplitude at each of a plurality of spectral elements,iteratively or recursively, until a convergence criterion is satisfied;generating, using a weight function, weight data from the input spectrum;generating a weighted input spectrum by weighting the input spectrum using the weight data;identifying neighboring spectral elements among the plurality of spectral elements in the weighted input spectrum,wherein a spectral element among the neighboring spectral elements comprises a maximum amplitude;updating the output spectrum by replacing an output amplitude at each of the neighboring spectral elements in the output spectrum with the input amplitude at each of the neighboring spectral elements in the input spectrum; andupdating the input spectrum by removing the input amplitude at each of the neighboring spectral elements in the input spectrum,determining, using an inverse transform, output seismic data from the updated output spectrum, andgenerating a seismic image of the subterranean region of interest based, at least in part, on the output seismic data; andusing a seismic interpretation workstation;identifying a drilling target of a hydrocarbon reservoir within the subterranean region of interest using, at least in part, the seismic image.2.The method of claim 1, further comprising planning, using a wellbore planning system, a wellbore path that penetrates the hydrocarbon reservoir.3.The method of claim 2, further comprising drilling, using a drilling system, a wellbore guided by the planned wellbore path.4.The method of claim 1, wherein the input seismic data is irregularly sampled.5.The method of claim 1, wherein the input seismic data is sparsely sampled.6.The method of claim 1, wherein the transform comprises a nonuniform discrete Fourier transform (NDFT) .7.The method of claim 1, wherein using the weight function comprises stacking the input amplitude at each spectral element among the plurality of spectral elements in the input spectrum along a radial line,wherein a first endpoint of the radial line is positioned at an origin of the input spectrum, andwherein a second endpoint of the radial line is positioned at each spectral element.8.The method of claim 1, wherein the neighboring spectral elements comprise a cube of spectral elements.9.The method of claim 1, wherein the input spectrum is organized into input spectrum slices,wherein each input spectrum slice among the input spectrum slices comprises the input amplitude at each of a plurality of spectral voxels.10.The method of claim 9, wherein iteratively or recursively, until the convergence criterion is satisfied, further comprises;generating, using the weight function, a weight data slice from each input spectrum slice;generating a weighted input spectrum slice by weighting each input spectrum slice using the weight data slice;identifying neighboring spectral voxels among the plurality of spectral voxels in the weighted input spectrum slice,wherein a spectral voxel among the neighboring spectral voxels comprises the maximum amplitude;updating the output spectrum by replacing the output amplitude at each of the neighboring spectral voxels in the output spectrum with the input amplitude at each of the neighboring spectral voxels in each input spectrum slice; andupdating each input spectrum slice by removing the input amplitude at each of the neighboring spectral voxels in each input spectrum slice.11.A system comprising;a seismic processing system configured to;initialize an output spectrum,receive, from a seismic acquisition system, input seismic data for a subterranean region of interest,wherein the input seismic data comprises a plurality of toxels,determine, using a transform, an input spectrum from the input seismic data,wherein the input spectrum comprises an input amplitude at each of a plurality of spectral elements,iteratively or recursively, until a convergence criterion is satisfied;generate, using a weight function, weight data from the input spectrum;generate a weighted input spectrum by weighting the input spectrum using the weight data;identify neighboring spectral elements among the plurality of spectral elements in the weighted input spectrum,wherein a spectral element among the neighboring spectral elements comprises a maximum amplitude;update the output spectrum by replacing an output amplitude at each of the neighboring spectral elements in the output spectrum with the input amplitude at each of the neighboring spectral elements in the input spectrum; andupdate the input spectrum by removing the input amplitude at each of the neighboring spectral elements in the input spectrum,determine, using an inverse transform, output seismic data from the updated output spectrum, andgenerate a seismic image of the subterranean region of interest based, at least in part, on the output seismic data; anda seismic interpretation workstation configured to;identify a drilling target of a hydrocarbon reservoir within the subterranean region of interest using, at least in part, the seismic image.12.The system of claim 11, further comprising a wellbore planning system configured to plan a wellbore path that penetrates the hydrocarbon reservoir.13.The system of claim 12, further comprising a drilling system configured to drill a wellbore guided by the planned wellbore path.14.A non-transitory computer-readable memory having computer-executable instructions stored thereon that, when executed by a computer processor, perform steps comprising;initializing an output spectrum;receiving, from a seismic acquisition system, input seismic data for a subterranean region of interest,wherein the input seismic data comprises a plurality of toxels;determining, using a transform, an input spectrum from the input seismic data;wherein the input spectrum comprises an input amplitude at each of a plurality of spectral elements;iteratively or recursively, until a convergence criterion is satisfied;generating, using a weight function, weight data from the input spectrum,generating a weighted input spectrum by weighting the input spectrum using the weight data,identifying neighboring spectral elements among the plurality of spectral elements in the weighted input spectrum;wherein a spectral element among the neighboring spectral elements comprises a maximum amplitude,updating the output spectrum by replacing an output amplitude at each of the neighboring spectral elements in the output spectrum with the input amplitude at each of the neighboring spectral elements in the input spectrum, andupdating the input spectrum by removing the input amplitude at each of the neighboring spectral elements in the input spectrum;determining, using an inverse transform, output seismic data from the updated output spectrum;generating a seismic image of the subterranean region of interest based, at least in part, on the output seismic data; andidentifying a drilling target of a hydrocarbon reservoir within the subterranean region of interest using, at least in part, the seismic image.15.The non-transitory computer-readable memory of claim 14, wherein the non-transitory computer-readable memory further stores instructions for a wellbore planning system that, when executed by the computer processor, perform a step comprising planning a wellbore path that penetrates the hydrocarbon reservoir.16.The non-transitory computer-readable memory of claim 14, wherein the input seismic data is irregularly sampled.17.The non-transitory computer-readable memory of claim 14, wherein using the weight function comprises stacking the input amplitude at each spectral element among the plurality of spectral elements in the input spectrum along a radial line,wherein a first endpoint of the radial line is positioned at an origin of the input spectrum, andwherein a second endpoint of the radial line is positioned at each spectral element.18.The non-transitory computer-readable memory of claim 14, wherein the neighboring spectral elements comprise a hypercube of spectral elements.19.The non-transitory computer-readable memory of claim 14, wherein the input spectrum is organized into input spectrum slices,wherein each input spectrum slice among the input spectrum slices comprises the input amplitude at each of a plurality of spectral voxels.20.The non-transitory computer-readable memory of claim 19, wherein iteratively or recursively, until the convergence criterion is satisfied, further comprises;generating, using the weight function, a weight data slice from each input spectrum slice;generating a weighted input spectrum slice by weighting each input spectrum slice using the weight data slice;identifying neighboring spectral voxels among the plurality of spectral voxels in the weighted input spectrum slice,wherein a spectral voxel among the neighboring spectral voxels comprises the maximum amplitude;updating the output spectrum by replacing the output amplitude at each of the neighboring spectral voxels in the output spectrum with the input amplitude at each of the neighboring spectral voxels in each input spectrum slice; andupdating each input spectrum slice by removing the input amplitude at each of the neighboring spectral voxels in each input spectrum slice.

Images