Image sharpening and spectrum enhancement based on geometric flow

Abstract

Examples of methods and systems are disclosed. The methods may include obtaining seismic data regarding a subsurface region of interest, wherein the seismic data includes a plurality of time-space waveforms. The methods may also include obtaining a seismic velocity model regarding the subsurface region of interest and determining a seismic image based, at least in part, on the plurality of time-space waveforms and the seismic velocity model. The methods may further include generating a filtered seismic image by applying a geometric flow filter to the seismic image, wherein applying the geometric flow filter comprises simulating a geometric flow using a tension spline function. The methods may still further include determining a geological boundary based, at least in part, on the filtered seismic image. The method may also include determining a drilling target in the subsurface region based, at least in part, on the geological boundary.

Description

IMAGE SHARPENING AND SPECTRUM ENHANCEMENT BASED ON GEOMETRIC FLOWBACKGROUND

[0001] In the oil and gas industry, seismic surveys are conducted over subsurface regions of interest during the search for, and characterization of, hydrocarbon reservoirs. In seismic surveys, a seismic source generates seismic waves that propagate through the subterranean region of interest and are detected by seismic receivers. The seismic receivers detect and may store a time-series of samples of earth motion caused by the seismic waves. The collection of time-series of samples recorded at many receiver locations generated by a seismic source at many source locations constitutes a seismic data set.

[0002] To determine the earth structure, including the presence of hydrocarbons, the seismic data set may be processed. Processing a seismic data set includes a sequence of steps designed to correct for a number of issues, such as near-surface effects, noise, irregularities in the seismic survey geometry, etc. Seismic data may be also processed to generate high resolution images of subsurface regions to assist in identification of fine geological features. A properly processed seismic data set may aid in decisions as to if and where to drill for hydrocarbons.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 disclosed herein relate to a method. The method includes obtaining seismic data regarding a subsurface region of interest, wherein the seismic data includes a plurality of time-space waveforms, obtaining a seismic velocity model regarding the subsurface region of interest, and determining a seismic  image based, at least in part, on the plurality of time-space waveforms and the seismic velocity model. The method further includes generating a filtered seismic image by applying a geometric flow filter to the seismic image, wherein applying the geometric flow filter includes simulating a geometric flow using a tension spline function, determining a geological boundary based, at least in part, on the filtered seismic image, and determining a drilling target in the subsurface region based, at least in part, on the geological boundary.

[0005] In general, in one aspect, embodiments disclosed herein relate to a system. The system includes a seismic acquisition system, a seismic processor, and a seismic interpretation workstation. The seismic acquisition system is configured to record seismic data regarding a subsurface region of interest, wherein the seismic data includes a plurality of time-space waveforms. The seismic processor is configured to obtain the seismic data and a seismic velocity model regarding the subsurface region of interest. The seismic processor is also configured to determine a seismic image based, at least in part, on the plurality of time-space waveforms and the seismic velocity model, and to generate a filtered seismic image by applying a geometric flow filter to the seismic image, wherein applying the geometric flow filter includes simulating a geometric flow using a tension spline function. The seismic interpretation workstation is configured to receive the filtered seismic image and to determine a geological boundary based, at least in part, on the filtered seismic image, and to determine a drilling target in the subsurface region based, at least in part, on the geological boundary.

[0006] It is intended that the subject matter of any of the embodiments described herein may be combined with other embodiments described separately, except where otherwise contradictory.

[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] FIG. 1 shows a seismic acquisition system of a subsurface region of interest, according to one or more embodiments of the present disclosure.

[0009] FIG. 2 shows examples of seismic data produced by a seismic acquisition system in accordance with one or more embodiments.

[0010] FIG. 3 shows an example of a seismic image for subsurface structure according to one or more embodiments.

[0011] FIG. 4 shows a drilling system in accordance with one or more embodiments.

[0012] FIG. 5 shows a flowchart in accordance with one or more embodiments.

[0013] FIG. 6 shows examples of filtered seismic images, according to one or more embodiments.

[0014] FIG. 7 shows examples of profiles of filtered seismic images, according to one or more embodiments.

[0015] FIG. 8 shows examples of spectra of filtered seismic images, according to one or more embodiments.

[0016] FIG. 9 depicts a schematic diagram of a computer system in accordance with one or more embodiments.DETAILED DESCRIPTION

[0017] 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.

[0018] 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.

[0019] In the following description of FIGs. 1-9, any component described regarding a figure, in various embodiments disclosed herein, may be equivalent to one or more like-named components described with regard to 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 with regard to a corresponding like-named component in any other figure.

[0020] 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 “a seismic signal” includes reference to one or more of such seismic signals.

[0021] 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.

[0022] 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.

[0023] In general, disclosed embodiments include systems and methods to filter seismic images for image sharpening using a geometric flow equation. In particular, some embodiments acquire seismic signals during geophysical explorations to map the structure  of a subsurface region. The seismic signals may be composed of diffracted and / or reflected seismic energy, due to the presence of discontinuities, subsurface interfaces, and / or diffractors among other geological features. Seismic imaging techniques aim at recovering the geological features from recorded seismic data.

[0024] However, seismic images generated from acquired seismic data may not have a sufficiently high resolution to identify geological features. The seismic data may not contain the necessary low and high frequency components necessary to accurately image fine geological structures. Commonly used techniques to improve the resolution of a seismic image include broadband data acquisition, pre-processing of seismic data before imaging, and high-accurately numerical imaging techniques. Another approach is implementing techniques to improve the generated seismic image (also called post-processing techniques) . These latter techniques may be preferable because of their reduced complexity and low-computational cost, however they may introduce numerical errors and noise to the seismic image. Methods and processing techniques to improve the resolution of seismic images with limited introduction of errors and noise may assist in improving the quality of the generated seismic images.

[0025] The resulting seismic image may then be used for seismic data interpretation, such as defining the spatial location and extent of a hydrocarbon reservoir. Thus, the disclosed methods are integrated into the established practical applications for improving seismic images and searching for an extraction of hydrocarbons from subsurface hydrocarbon reservoirs. The disclosed methods represent an improvement over existing methods for at least the reasons of lower cost and increased efficacy.

[0026] FIG. 1 shows a seismic acquisition system (100) of a subsurface region of interest (102) , according to one or more embodiments. In some cases, the subsurface region of interest (102) may lie beneath a lake, sea, or ocean. In other cases, the subsurface region of interest (102) may lie beneath an area of dry land. The subsurface region of interest (102) may contain a hydrocarbon deposit (120) that may form part of a hydrocarbon reservoir (104) . The seismic acquisition system (100) may utilize a seismic source (106) that generates radiated seismic waves (108) . The type of seismic source (106) may depend on  the environment in which it is used. For example, on land the seismic source (106) may be a vibroseis truck or an explosive charge, but in water the seismic source (106) may be an airgun. The radiated seismic waves (108) may return to the surface as refracted seismic waves (110) or reflected seismic waves (114) .

[0027] Refracted seismic waves (110) and reflected seismic waves (114) may occur, for example, due to geological discontinuities (112) that may be also known as “seismic reflectors” . The geological discontinuities (112) may be, for example, planes or surfaces that mark changes in physical or chemical characteristics in a geological structure. The geological discontinuities (112) may also be boundaries between faults, fractures, or groups of fractures within a rock. The geological discontinuities (112) may delineate a hydrocarbon reservoir (104) .

[0028] At the surface, refracted seismic waves (110) and reflected seismic waves (114) may be detected by seismic receivers (116) . Radiated seismic waves (108) that propagate from the seismic source (106) directly to the seismic receivers (116) , known as direct seismic waves (122) , are also detected by the seismic receivers (116) .

[0029] In some embodiments, a seismic source (106) may be positioned at a location denoted (xs, ys) , where x and y represent orthogonal axes on the earth’s surface above the subsurface region of interest (102) . The seismic receivers (116) may be positioned at a plurality of seismic receiver locations denoted (xr, yr) , with the distance between each receiver and the source being termed “the source-receiver offset” , or simply “the offset” . Thus, the direct seismic waves (122) , refracted seismic waves (110) , and reflected seismic waves (114) generated by a single activation of the seismic source (106) may be represented in the axes (xs, ys, xr, yr, t) . The t-axis indicates the recording time between the activation of the seismic acquisition system (100) and the sample time at which the seismic wave is detected by the seismic receivers (116) .

[0030] Seismic processing may reduce five-dimensional seismic data produced by a seismic acquisition system (100) to three-dimensional (x, y, t) seismic data by, for example, correcting the recorded time for the time of travel from the seismic source (106) to the seismic receiver (116) and summing ( “stacking” ) samples over two horizontal space  dimensions. Stacking of samples over a predetermined time interval may be performed as desired, for example, to reduce noise and improve the quality of the signals.

[0031] Seismic data may also refer to data acquired at different time intervals, such as, for example, in cases where seismic surveys are repeated after a period of weeks, months, or years, to obtain time-lapse data. Seismic data may also be pre-processed or partially-processed data, e.g., arranged as “common shot gathers” (CSG) , i.e., sorting waveforms as acquired by different receivers and having a single source location. The type of seismic data is not intended as limiting, and any other suitable seismic data is intended to fall within the scope of the present disclosure.

[0032] FIG. 2 shows examples of seismic data (202) produced by a seismic acquisition system (100) in accordance with one or more embodiments. An example of a CSG (204) depicts direct seismic waves (122) , refracted seismic waves (110) , and reflected seismic waves (114) generated by a single activation of the seismic source (106) and recorded by a plurality of seismic receivers (116) deployed, for example in a line, on the surface of the earth. Each seismic receiver (116) may record a time-series representing the amplitude of ground-motion at a sequence of discrete times. This time-series may be denoted or otherwise referred to as a “waveform” . Seismic data therefore may include a plurality of time-space waveforms (205) associated to the plurality of seismic receivers (116) .

[0033] In the CSG (204) shown in FIG. 2 the vertical axis represents the time scale (206) and the horizontal axis represents the offset (208) . In some embodiments, direct seismic waves (122) , refracted seismic waves (110) , and reflected seismic waves (114) may be identified in the CSG (204) by their arrival times, i.e., the times at which they are first detected by the seismic receivers (116) . The location of a particular type of wave in seismic data (202) acquired in time and space, such as in CSG (204) , may be termed as an “arrival” or as an “event” .

[0034] The CSG (204) illustrates how the arrivals are detected at later times by the seismic receivers (116) that are farther from the seismic source (106) . In some embodiments, arrivals of direct seismic waves (122) in the CSG (204) may be characterized by a straight line, while arrivals of reflected seismic waves (114) may present a hyperbolic shape, as  seen in FIG. 2. Refracted seismic waves (110) may be characterized by arrivals approximating a straight line in offset-time.

[0035] In one or more embodiments, seismic data (202) acquired by a seismic acquisition system (100) may be arranged in a plurality of CSGs (210) to create a 3D seismic dataset. Alternatively, the seismic data may be represented as a “seismic volume” (212) consisting of a plurality of time-space waveforms with a time axis (214) , a first spatial dimension (216) , and a second spatial dimension (218) , where the first spatial dimension (216) and second spatial dimension (218) are orthogonal and span the Earth’s surface above the subsurface region of interest (102) .

[0036] Seismic data (202) may be processed by a seismic processor (220) to generate a seismic velocity model (219) of the subterranean region of interest (102) . A seismic velocity model (219) is a representation of seismic velocity at a plurality of locations within a subterranean region of interest (102) . Seismic velocity is the speed at which a seismic wave, that may be a pressure-wave or a shear-wave, travel through a medium. Pressures waves are often referred to as “primary-waves” or “P-waves” . Shear waves are often referred to a “secondary-waves” or “S-waves” . Seismic velocities in a seismic velocity model (219) may vary in vertical depth, in one or more horizontal directions, or both. Layers of rock are created from different materials or created under varying conditions. Each layer of rock may have different physical properties from neighboring layers and these different physical properties may include seismic velocity.

[0037] FIG. 2 schematically illustrates that in some embodiments seismic data (202) may be processed by a seismic processor (220) to generate a seismic image (230) of the subterranean region of interest (102) . For example, a time-domain seismic image (232) may be generated using a process called seismic migration (also referred to as simply “migration” herein) using a seismic velocity model (219) . In seismic migration, seismic events (e.g., reflections, refractions) recorded at the surface are relocated in either time or space to the location the event occurred in the subsurface. In some embodiments, migration may transform pre-processed shot gathers from a time-domain to a depth-domain seismic image (234) . In a depth-domain seismic image (234) , seismic events in a migrated shot  gather may represent geological boundaries (236, 238) in the subsurface. Various types of migration algorithms may be used in seismic imaging. For example, one type of migration algorithm corresponds to reverse time migration.

[0038] A seismic image (230) of high resolution may be obtained if densely-recorded data is acquired by using closely-spaced shots and seismic receivers. For example, seismic waves with a bandwidth extending up to 100 Hz or more may resolve thin features. However, acquiring densely-recorded data may be costly and in some cases not possible due to the presence of obstacles in the subsurface region of interest (102) . As a result, the the resolution of the seismic image (230) may not high enough to identify and determine fine geological attributes in the subsurface regions of interest (102) .

[0039] FIG. 3 illustrates a two-dimensional seismic image. The seismic image (300) is an example of seismic images used to determine seismic attributes regarding a subsurface region of interest (102) . A distance from the surface (depth) is indicated on the vertical axis (302) and a horizontal distance is indicated on the horizontal axis (304) . The grayscale indicates the initial amplitude range (306) of the seismic image. It is frequently valuable to determine the presence, and if present the location, of a geological fault from such seismic images. Continuous (unfaulted) geological interfaces produce continuous and smoothly varying bands of high or low reflection amplitudes. The location of the geological fault in the seismic image (300) may be determined by breaks or vertical displacements ( “throws” ) in such continuous, smoothly varying bands, such as the fault (308) . However, the bandlimited nature of the seismic image (300) means the location of these faults cannot always be precisely determined and some small anomalies (309) may be overlooked or misinterpreted.

[0040] As illustrated in FIGs. 2 and 3, processing of an observed seismic dataset (202) may generate a seismic image (230) that may reveal the two or three-dimensional geometry of a subsurface region of interest (102) . In particular, the geological boundaries (236, 238) may delineate a hydrocarbon reservoir (104) . Identifying geological boundaries (236, 238) and other geological objects, such as faults, may be performed using a seismic interpretation workstation. If a seismic image (230) indicates the potential presence of  hydrocarbons in the subsurface region of interest (102) , a wellbore (118) may be planned with a wellbore planning system. Further, a drilling system may drill a wellbore (118) to confirm the presence of those hydrocarbons.

[0041] FIG. 4 shows a drilling system (400) in accordance with one or more embodiments. As shown in FIG. 4, a wellbore (118) following a wellbore trajectory (404) may be drilled by a drill bit (406) attached by a drillstring (408) to a drilling rig (410) located on the surface (124) of the earth. The drilling rig (410) may include framework, such as a derrick (414) to hold drilling machinery. A crown block (411) may be mounted at the top of the derrick (414) , and a traveling block (413) may hang down from the crown block (411) by means of a cable (415) or drilling line. One end of the cable (415) may be connected to a drawworks (not shown) , which is a reeling device that may be used to adjust the length of the cable (415) so that the traveling block (413) may move up or down the derrick (414) .

[0042] A top drive (416) provides clockwise torque via the drive shaft (418) to the drillstring (408) in order to drill the wellbore (118) . The drillstring (408) may comprise a plurality of sections of drillpipe attached at the uphole end to the drive shaft (418) and downhole to a bottomhole assembly ( “BHA” ) (420) . The BHA (420) may be composed of a plurality of sections of heavier drillpipe and one or more measurement-while-drilling ( “MWD” ) tools configured to measure drilling parameters, such as torque, weight-on-bit, drilling direction, temperature, etc., and one or more logging-while-drilling ( “LWD" ) tools configured to measure parameters of the rock surrounding the wellbore (118) , such as electrical resistivity, density, sonic propagation velocities, gamma-ray emission, etc. MWD and logging tools may include sensors and hardware to measure downhole drilling parameters, and these measurements may be transmitted to the surface (124) using any suitable telemetry system known in the art. The BHA (420) and the drillstring (408) may include other drilling tools known in the art but not specifically listed.

[0043] The wellbore (118) may traverse a plurality of overburden (422) layers and one or more seals or cap-rock formations (424) to a hydrocarbon reservoir (104) within the subterranean region (428) , and specifically to a drilling target (430) within the  hydrocarbon reservoir (104) . The wellbore trajectory (404) may be a curved or a straight. All or part of the wellbore trajectory (404) may be vertical, and some portions of the wellbore trajectory (404) may be deviated from the vertical or horizontal. One or more portions of the wellbore (118) may be cased with casing (432) in accordance with a wellbore plan.

[0044] To start drilling, or “spudding in” the well, the hoisting system lowers the drillstring (408) suspended from the derrick (414) towards the planned surface location of the wellbore (118) . An engine, such as an electric motor, may be used to supply power to the top drive (416) to rotate the drillstring (408) through the drive shaft (418) . The weight of the drillstring (408) combined with the rotational motion enables the drill bit (406) to bore the wellbore (118) .

[0045] The drilling system (400) may be disposed at and communicate with other systems in the well environment, such as a seismic processor (220) , an interpretation workstation (440) , and a wellbore planning system (438) . The drilling system (400) 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 drilling system (400) may receive well-measured data from one or more sensors and / or logging tools arranged to measure controllable parameters of the drilling operation. During operation of the drilling system (400) , the well-measured data may include mud properties, flow rates, drill volume and penetration rates, rock physical properties, etc.

[0046] In some embodiments, the rock physical properties may be used by the seismic interpretation workstation (440) to determine a location of a hydrocarbon reservoir (104) (or other subterranean features) , including the drilling target (430) . Knowledge of the existence and location of the hydrocarbon reservoir (104) and other subterranean features may be transferred from the seismic interpretation workstation (440) to a wellbore planning system (438) . The wellbore planning system (438) may use information regarding the hydrocarbon reservoir (104) location to plan a well, including a planned wellbore trajectory (404) from the surface (124) of the earth to penetrate the hydrocarbon reservoir (104) . In addition, to the depth and geographic location of the hydrocarbon  reservoir (104) , the planned wellbore trajectory (404) may be constrained by surface limitations, such as suitable locations for the surface position of the wellhead, i.e., the location of potential or preexisting drilling rigs, drilling ships or from a natural or man-made island.

[0047] Typically, the wellbore plan is generated based on best available information at the time of planning from a geophysical model, geomechanical models encapsulating subterranean stress conditions, the trajectory of any existing wellbores (which it may be desirable to avoid) , and the existence of other drilling hazards, such as shallow gas pockets, over-pressure zones, and active fault planes. Information regarding the planned wellbore trajectory (404) may be transferred to the drilling system (400) described in FIG. 4. The drilling system (400) may drill the wellbore (118) along the planned wellbore trajectory (404) to access the drilling target (430) in the hydrocarbon reservoir (104) .

[0048] Turning to FIG. 5, FIG. 5 shows a flowchart in accordance with one or more embodiments. Specifically, FIG. 5 describes an embodiment of the inventive method to filter a seismic image using a geometric flow filter. While the various blocks in FIG. 5 are presented and described sequentially, one of ordinary skill in the art will appreciate that some or all of the blocks may be executed in different orders, may be combined or omitted, and some or all of the blocks may be executed in parallel. Furthermore, the blocks may be performed actively or passively.

[0049] In Block 500, seismic data (202) associated with a subsurface region of interest (102) may be obtained, in accordance with one or more embodiments. The seismic data (202) may be acquired using a seismic acquisition system (100) above the subsurface region of interest (102) . The seismic data (202) may be processed to attenuate noise and may be organized in one or more spatial dimensions (216, 218) and a time axis (214) to form a plurality of time-space waveforms (205) . In some embodiments, one or more CSGs (210) may be generated with the source position corresponding to the middle of the offset, as illustrated in FIG. 2.

[0050] In Block 510, a seismic velocity model (219) regarding the subsurface region of interest (102) is obtained, in accordance with one or more embodiments. The seismic  velocity model (219) provides an estimate of at least one seismic wave propagation velocity at each location in the depth domain within the subterranean region of interest (102) . Typically, a seismic velocity model (219) is specified by at least one seismic velocity for a particular wave type at a plurality of discrete grid points spanning the subsurface region of interest, but other specifications are possible. For example, the seismic velocity model (219) may be defined by a plurality of continuously varying mathematical functions.

[0051] An accurate velocity model may be useful for accurate simulations of seismic wave propagation in the subsurface. A velocity model may be constructed by processing recorded seismic data (202) obtained using a seismic acquisition system (100) . Processing seismic data (202) to obtain a velocity model may be considered an inverse problem, where the applied process must determine the subsurface velocity model that resulted in the recorded seismic data (202) . The various processes and techniques used to process seismic data (202) to form a velocity model may generally be categorized as either a “data-domain approach” , such as full-waveform inversion (FWI) , or an “image-domain approach” , such as migration velocity analysis. Machine learning models and deep learning frameworks may be also implemented to determine a velocity model from given seismic data (202) .

[0052] In Block 520 of Fig. 5, a seismic image (230) may be determined based on the plurality of time-space waveforms (205) and the seismic velocity model (219) , in accordance with one or more embodiments. Various types of algorithms may be used in the generation of the seismic image (230) . For example, the plurality of time-space waveforms (205) may include a source wavefield and a receiver wavefield. The source wavefield may include one or more source functions at one or more source locations (xs, ys) used during activation of the seismic source (106) . The receiver wavefield may include the plurality of waveforms acquired at the seismic receivers (116) .

[0053] In some embodiments, the seismic image (230) may be a migrated seismic image. The migrated seismic image may be generated with reverse time migration (RTM) using the two-way wave equation. In RTM, the source wavefield may be obtained by forward  modelling the propagation of a synthetic source function using the seismic velocity model (219) . A receiver wavefield may be generated using the same seismic velocity model (219) by backward propagating in time the receiver wavefield. In other words, the receiver wavefield may be first reversed in time and the used as a source function applied at the corresponding seismic receivers (116) to simulate a radiated wavefield.

[0054] The migrated seismic image may then be formed by applying an imaging condition to the receiver wavefield and the source wavefield. In some embodiments, the first imaging condition may be represented by a cross‐correlation between the source wavefield with the receiver wavefield under the basic assumption that the source wavefield represents the down‐going wave‐field and the receiver wave‐field the up‐going wave‐field.

[0055] In some embodiments the seismic image (230) may be obtained by merging, or stacking, different partial migrated seismic images. Each partial migrated seismic image may be generated from seismic data (202) acquired upon activation of one or more seismic sources (106) .

[0056] The seismic image (230) may map the spatial variation of a seismic attribute. Non-limiting examples of seismic attributes mapped by a seismic image (230) include acoustic impedance, wave velocity, density, or any other hydrocarbon indicator.

[0057] Keeping with to FIG. 5, in Block 530, a filtered seismic image is generated by applying a geometric flow filter to the seismic image (230) , where the geometric flow filter includes simulating a geometric flow using a tension spline function, in accordance with one or more embodiments. Application of a post-processing filtering technique may provide an effective and computational efficient method to improve the resolution of the seismic image (230) . For example, in a non-limiting embodiment, filtering may include applying image segmentation to decompose an image domain into approximately homogenous regions, that are separated by sharp boundaries. Image segmentation may be implemented to recover target boundaries even in noisy images or signals. Problems of sharpening of boundaries may be solved by image processing with Partial Differential Equations (PDEs) , i.e., with mathematical models that simulate a kind of flow, such as a  fluid flow, crystal growth, or geometrical flow. The simulations are parametrized by the variable t playing the role of time in the physical context.

[0058] In particular, geometric flow equations are parabolic PDEs that may describe the deformation of surfaces driven by their curvature in various forms. In some embodiments, a geometrical flow filter may include simulations of surfaces that undergo continuous deformations driven by their mean-curvature. Intuitively, surfaces may evolve driven by mean-curvature if the normal component of the velocity at a moving point of the surface is given by the mean curvature of the surface. For example, if evolving based on mean curvature, a sphere would shrink inwardly since its mean curvature vector points inward.

[0059] In one or more embodiments, simulating a geometric flow may include numerical simulations of surface mean-curvature time-evolution, as shown in Block 534. Surface mean-curvature time-evolution may be simulated by solving the following non-linear PDE, known as the Allen-Cahn equation:

[0060] where u represents an updating intensity function starting at an initial image intensity u0, and is a function of both time t and space coordinates (x, y, z) . The constant ε is a diffusion coefficient, that determines how fast the updating intensity function u diffuses in space. The larger the value of ε, the faster the spread and decay of u. A non-limiting example of the constant ε, associated with slow decay and improved accuracy, is ε=0.5. The Allen-Cahn equation was initially developed for applications of interface propagation of binary phase fluid or multi-phase fluids. In applications of image processing, the function u may bring the initial image intensity u0 as close as possible to boundaries, and the boundaries are progressively sharpened with time.

[0061] In some embodiments, the initial image intensity u0 may be the intensity of a normalized image. As shown in Block 532, the normalized image may be obtained by mapping the initial amplitude range of the seismic image to a predetermined amplitude range. An example of a two-dimensional normalized image (610) is shown in FIG. 6, according to one or more embodiments. In all seismic images of FIG. 6 the depth is  indicated on the vertical axis (602) , and the horizontal distance is indicated on the horizontal axis (604) . The normalized image (610) is obtained from mapping the amplitude range (306) of the seismic image (300) into the amplitude range [-1, 1] (606) as shown in the scale bar.

[0062] Eq. (1) is a non-linear parabolic PDE, since the last term on the right hand side is non-linear. PDEs such as Equation (1) may be solved with ordinary numerical methods such as finite-differences or finite-elements, which are focused on discretizing the PDE in its different space dimensions. The discretized PDE results in a system of linear equations, that may be efficiently solved with standard methods of linear algebra. The obtained numerical solution is a solution of the discretized PDE, and an approximation to the solution of the original, continuous, PDE. The difference between the solution of the discretized PDE and the solution of the continuous PDE is known as the discretization error. When the spatial gradients (variation of the function u in space) are high, discretization methods such as finite-differences and finite elements are known for being affected by numerical dispersion, which increases the discretization error. Numerical dispersion results when the velocity of a propagating wave according to the discrete PDE deviates from the corresponding velocity of the continuous PDE.

[0063] In some embodiments, the accuracy of the numerical solution may be increased by approximating the function u as a piece-wise function, composed of a series of splines, where each spline connects two nodes. The simplest spline is a linear function. However, when accuracy of derivatives of different order (and therefore gradients) is desired, a series of cubic polynomials, referred to as cubic splines, may be used to approximate the function u. In accordance with one or more embodiments, the numerical solution of Equation (1) is expressed in terms of a tension spline function. In one-dimensional space, the tension spline function is obtained by minimizing the following variational:

[0064] where σ is the tensor factor and a and b are the limits of integration. The solution to Equation (2) is given by u (x) =a+bx+c sinh (σx) +d cosh (σx)   (3)

[0065] where a, b, c, and d are constants determined by the boundary conditions. In the particular case of σ approaching zero, the variational (2) reduces to:

[0066] the solution of which is the cubic spline. Compared to the cubic spline, the tension spline function can reduce the number of saddle points of the curve as well as the arc length.

[0067] In some embodiments, the tension spline function includes a large tensor factor σ, as shown in Block 535 of FIG. 5. When the tensor factor σ is large, the piece-wise functions approach linear functions. When the tensor factor σ is small, the piece-wise functions approach the cubic spline. As a non-limiting example, the tensor factor σ may take values between 10 and 20, which may be appropriate for preserving faults and sharp structures. Because two resulting tension spline functions on either side of a node have the same first and second derivative values, the obtained numerical solution is affected with less high-frequency noises and less numerical dispersion.

[0068] As stated above, geometric flow simulations may be performed by solving the equation describing the surface mean-curvature time-evolution, given by Equation (1) . The numerical update intensity function u may be obtained in terms of one or more spatial dimensions (x, y, z) and time, t ∈ [0, T] where T is the total duration of the simulation. A simulated seismic image may be obtained when reaching the total duration of the simulation. According to some embodiments Equation (1) may be solved by first discretizing it in space. Discretizing in space may involve the construction of a spatial grid related to the subsurface region of interest (102) . As such, a spatial grid may be two-dimensional or three-dimensional, where grid nodes provide discretized local solutions for the update intensity function u during a simulation. The discretized local solutions for the update intensity function u may be expressed in terms of the tension spline function. The grid data may be obtained from a parameter file, where the parameter file includes grid node numbers and grid space dimensional information. While the grid data may  correspond to a regular grid, irregular grids are also contemplated with nonuniform grids, tree-based data structures, or triangulated meshes.

[0069] In accordance with one or more embodiments, a PDE that has been discretized in space may be discretized in time by, for example, subdividing the simulation interval [0, T] into time steps, Δt. In some embodiments, Δt may be constant and in other embodiments, Δt may vary throughout the simulation. The discretization may result in a system of discrete differential equations that may be solved for any simulation interval [0, T] using an explicit or iterative time integration scheme. The time integration scheme may include approximating the discrete differential equations with finite-differences, as shown in Block 536. Other time-marching or time-discretization methods known in the art may be implemented without departing from the scope of the present disclosure. Further, in some embodiments the simulation interval [0, T] may be predetermined. In other embodiments, the simulation may be performed by iterating in time until a stopping condition in reached, as shown in Block 538. The stopping condition may be a user’s choice and may be, for example, a criterion regarding the resolution of the simulated seismic image, or a criterion regarding the frequency band of the spectrum of the simulated seismic image.

[0070] As shown in Block 539, in some embodiments, the amplitude range of the simulated seismic image may be mapped to the initial amplitude range of the seismic image. The resulting filtered seismic image may then be more easily interpreted based on the initial amplitude range.

[0071] Returning to FIG. 6, FIG. 6 shows examples of filtered seismic images that demonstrate the capabilities of the proposed method, according to one or more embodiments. Seismic images (620) , (630) , (640) , and (650) correspond to the results of filtering the normalized seismic image (610) with the geometric flow filter of Equation (1) , after 41, 81, 121 and 161 time steps, respectively. The amplitude range (607) of filtered seismic images (620) - (650) corresponds to the initial amplitude range (306) of the seismic image (300) . It can be observed that as the number of time steps increases, the discontinuities in the filtered seismic images (620) - (650) remain at the right locations  and are better defined. Major structures and layers are well separated by sharpened boundaries. For example, the extension of the major fault (608) with depth can be more precisely determined. On the other hand, the location of the minor fault (609) along the depth dimension (602) becomes more visible.

[0072] The sharpening effects of the proposed method are further illustrated in FIG. 7, that shows amplitude vertical profiles (611) extracted from seismic images (610) - (650) of Fig. 6. In FIG. 7 the vertical axis indicates the amplitude (702) , and the horizontal axis indicates depth (704) . The vertical profiles show that as the time step of the geometrical flow filter increases, the slopes (710) of the amplitudes at the discontinuities approach infinite values. On other words, the vertical profiles approach the form of a piecewise block function (706) , resembling block structures that are commonly seen in geological formations. Small variations (708) of the vertical profiles are “smoothed out” by the processing with the geometrical flow filter.

[0073] In order to illustrate the effects of the proposed method on the frequency content of the seismic signal, the spectra of the vertical profiles of FIG. 7 are shown in FIG. 8. The vertical axis indicates the spectral amplitude (802) , and the horizontal axis the frequency in Hz (804) . Fig. 8 indicates that amplitudes at both, the low frequency end (806) and high frequency end (808) , of the frequency band are amplified by processing with the geometric flow filter. Therefore, the geometric flow filter appears to be effective in enhancing the spectral content of the seismic images, and thus, improving its resolution.

[0074] Returning to FIG. 5, in Block 540, a geological boundary is determined based on the filtered seismic image, in accordance with one or more embodiments. The seismic processor (220) may transfer the filtered seismic image to a seismic interpretation workstation (440) . The seismic interpretation workstation (440) may use the knowledge of the filtered seismic image to identify geological boundaries (236, 238) and geological structures. Examples of geological structures include, but are not limited to, faults, salt bodies and salt caverns. A filtered seismic image with improved resolution in space may facilitate the description of fine geological structures delimited by geological boundaries (236, 238) . In some embodiments, description of geological structures may assist in  monitoring storage of chemicals, such as for example, CO2 or nitrogen, in deep underground geologic formations. Detailed description of fine geological structures may be important information to be used in the prevention of leakage of chemicals from the underground geologic formations.

[0075] In Block 550, a drilling target in the subsurface region is determined based, at least in part, on the geological boundary, in accordance with one or more embodiments. A drilling target (430) in a wellbore (118) may be determined by the seismic interpretation workstation (440) , and may be based on, for example, an expected presence of gas or another hydrocarbon. Locations in a seismic image (230) may be delimited by geological boundaries (236, 238) and may indicate a probability of the presence of a hydrocarbon. Locations in a seismic image (230) may indicate an elevated probability of the presence of a hydrocarbon and may be targeted by well designers. On the other hand, locations in a seismic image (230) indicating a low probability of the presence of a hydrocarbon may be avoided by well designers.

[0076] In Block 560, a wellbore trajectory to intersect the drilling target is planned using the filtered seismic image, in accordance with one or more embodiments. Knowledge of the location of the drilling target (430) and the filtered seismic image may be transferred by the seismic interpretation workstation (440) to a wellbore planning system (438) . Instructions associated with the wellbore planning system (438) may be stored, for example, in the memory (909) within the computer system (900) described in FIG. 9 below. The wellbore planning system (438) may use the knowledge of the location of the drilling target (430) and of filtered the seismic image to plan a wellbore trajectory (404) within the subterranean region of interest (102) .

[0077] In Block 570, a wellbore is drilled guided by the planned wellbore trajectory, in accordance with one or more embodiments. The wellbore planning system (438) may transfer the planned wellbore trajectory (404) to the drilling system (400) described in FIG. 4. The drilling system (400) may drill a portion of the wellbore (118) along the planned wellbore trajectory (404) to access and produce the hydrocarbon reservoir (104) to the surface (124) .

[0078] In some embodiments the wellbore planning system (438) , the seismic interpretation workstation (440) , and the seismic processor (220) may each be implemented within the context of a computer system. FIG. 9 is a block diagram of a computer system (900) used to provide computational functionalities associated with described algorithms, methods, functions, processes, flows, and procedures as described in the instant disclosure, according to an implementation. The illustrated computer (900) is intended to encompass any computing device such as a high performance computing (HPC) device, 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 (900) may include a computer that includes an input device, such as a keypad, keyboard, touch screen, or other device that can accept user information, and an output device that conveys information associated with the operation of the computer (900) , including digital data, visual, or audio information (or a combination of information) , or a GUI.

[0079] The computer (900) can serve in a role as a client, network component, a server, a database or other persistency, or any other component (or a combination of roles) of a computer system for performing the subject matter described in the instant disclosure. The illustrated computer (900) is communicably coupled with a network (902) . In some implementations, one or more components of the computer (900) may be configured to operate within environments, including cloud-computing-based, local, global, or other environment (or a combination of environments) .

[0080] At a high level, the computer (900) is an electronic computing device operable to receive, transmit, process, store, or manage data and information associated with the described subject matter. According to some implementations, the computer (900) 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) .

[0081] The computer (900) can receive requests over network (902) from a client application (for example, executing on another computer (900) ) and responding to the received requests by processing the said requests in an appropriate software application. In addition, requests may also be sent to the computer (900) 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 computers.

[0082] Each of the components of the computer (900) can communicate using a system bus (903) . In some implementations, any or all of the components of the computer (900) , both hardware or software (or a combination of hardware and software) , may interface with each other or the interface (904) (or a combination of both) over the system bus (903) using an application programming interface (API) (907) or a service layer (908) (or a combination of the API (907) and service layer (908) . The API (907) may include specifications for routines, data structures, and object classes. The API (907) 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 (908) provides software services to the computer (900) or other components (whether or not illustrated) that are communicably coupled to the computer (900) . The functionality of the computer (900) may be accessible for all service consumers using this service layer (908) . Software services, such as those provided by the service layer (908) , 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 other suitable format. While illustrated as an integrated component of the computer (900) , alternative implementations may illustrate the API (907) or the service layer (908) as stand-alone components in relation to other components of the computer (900) or other components (whether or not illustrated) that are communicably coupled to the computer (900) . Moreover, any or all parts of the API (907) or the service layer (908) 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.

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

[0084] The computer (900) includes at least one computer processor (905) . Although illustrated as a single computer processor (905) in FIG. 9, two or more processors may be used according to particular needs, desires, or particular implementations of the computer (900) . Generally, the computer processor (905) executes instructions and manipulates data to perform the operations of the computer (900) and any algorithms, methods, functions, processes, flows, and procedures as described in the instant disclosure.

[0085] The computer (900) also includes a memory (909) that holds data for the computer (900) or other components (or a combination of both) that may be connected to the network (902) . For example, memory (909) may be a database storing data consistent with this disclosure. Although illustrated as a single memory (909) in FIG. 9, two or more memories may be used according to particular needs, desires, or particular implementations of the computer (900) and the described functionality. While memory (909) is illustrated as an integral component of the computer (900) , in alternative implementations, memory (909) may be external to the computer (900) .

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

[0087] There may be any number of computers (900) associated with, or external to, a computer system containing computer (900) , each computer (900) communicating over network (902) . 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 one computer (900) , or that one user may use multiple computers (900) .

[0088] In some embodiments, the computer (900) is implemented as part of a cloud computing system. For example, a cloud computing system may include one or more remote servers along with various other cloud components, such as cloud storage units and edge servers. In particular, a cloud computing system may perform one or more computing operations without direct active management by a user device or local computer system. As such, a cloud computing system may have different functions distributed over multiple locations from a central server, which may be performed using one or more Internet connections. More specifically, cloud computing system may operate according to one or more service models, such as infrastructure as a service (IaaS) , platform as a service (PaaS) , software as a service (SaaS) , mobile "backend" as a service (MBaaS) , serverless computing, artificial intelligence (AI) as a service (AIaaS) , and / or function as a service (FaaS) .

[0089] 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:obtaining seismic data regarding a subsurface region of interest, wherein the seismic data comprises a plurality of time-space waveforms;obtaining a seismic velocity model regarding the subsurface region of interest;determining a seismic image based, at least in part, on the plurality of time-space waveforms and the seismic velocity model;generating a filtered seismic image by applying a geometric flow filter to the seismic image, wherein applying the geometric flow filter comprises simulating a geometric flow using a tension spline function;determining a geological boundary based, at least in part, on the filtered seismic image; anddetermining a drilling target in the subsurface region based, at least in part, on the geological boundary.2.The method of claim 1, further comprising planning, using a wellbore planning system, a planned wellbore trajectory to intersect the drilling target.3.The method of claim 2, further comprising drilling, using a drilling system, a portion of a wellbore guided by the planned wellbore trajectory.4.The method of claim 1, wherein simulating the geometric flow comprises simulating a surface mean-curvature time-evolution.5.The method of claim 1, wherein using a tension spline function comprises using a large tensor factor σ.6.The method of claim 1, wherein generating the filtered seismic image comprises mapping an initial amplitude range of the seismic image to a predetermined amplitude range.7.The method of claim 1, wherein simulating the geometric flow comprises generating a simulated seismic image by iterating in time until a stopping condition is reached.8.The method of claim 7, wherein generating the filtered seismic image comprises mapping an amplitude range of the simulated seismic image to an initial amplitude range of the seismic image.9.The method of claim 1, wherein simulating the geometric flow comprises discretizing a PDE in time using finite-differences.10.The method of claim 1, wherein determining the seismic image comprises performing reverse time migration.11.A system comprising:a seismic acquisition system configured to record seismic data regarding a subsurface region of interest, wherein the seismic data comprises a plurality of time-space waveforms;a seismic processor configured to:obtain the seismic data,obtain a seismic velocity model regarding the subsurface region of interest,determine a seismic image based, at least in part, on the plurality of time-space waveforms and the seismic velocity model, andgenerate a filtered seismic image by applying a geometric flow filter to the seismic image, wherein applying the geometric flow filter comprises simulating a geometric flow using a tension spline function; anda seismic interpretation workstation configured to receive the filtered seismic image and to:determine a geological boundary based, at least in part, on the filtered seismic image, anddetermine a drilling target in the subsurface region based, at least in part, on the geological boundary.12.The system of claim 11, further comprising a wellbore planning system configured to plan a planned wellbore trajectory to intersect the drilling target.13.The system of claim 12, further comprising a drilling system configured to drill a portion of a wellbore guided by the planned wellbore trajectory.14.The system of claim 11, wherein the geometric flow comprises surface mean-curvature evolution.15.The system of claim 11, wherein the seismic processor is further configured to map an initial amplitude range of the seismic image to a predetermined amplitude range.16.The system of claim 11, wherein the seismic processor is further configured to generate a simulated seismic image by iterating in time until a stopping condition is reached.17.The system of claim 16, wherein the seismic processor is further configured to map an amplitude range of the simulated seismic image to an initial amplitude range of the seismic image.18.The system of claim 11, wherein the seismic processor is further configured to discretize a PDE in time using finite-differences.19.The system of claim 11, wherein the seismic processor is further configured to perform reverse time migration.20.The system of claim 11, wherein the tension spline function comprises a large tensor factor σ.

Images