Method for phase reversal correction of diffraction images

Abstract

Examples of methods are disclosed. The methods may include obtaining seismic data regarding a subsurface region of interest, wherein the seismic data is a function of a plurality of coordinates. The methods may also include obtaining a velocity model regarding the subsurface region; forming a plurality of dip-angle gathers; and generating a filtering mask, wherein the filtering mask is a function of the plurality of coordinates. The methods may further include generating, for each dip-angle gather, a transformed gather by applying a transformation to the dip-angle gather; updating the filtering mask based, at least in part, on an energy extremum of the transformed gather; and forming a filtered diffraction image based, at least in part, on the plurality of dip-angle gathers and the updated filtering mask. The methods may still further include determining a drilling target in the subsurface region based on the filtered diffraction image. Examples of systems are also disclosed.

Description

METHOD FOR PHASE REVERSAL CORRECTION OF DIFFRACTION IMAGESBACKGROUND

[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 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 seismic data includes the generation of a seismic velocity model and seismic migration. 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. Another step in processing a seismic data may be the imaging of diffractors associated with local heterogeneities such as faults and fracture zones. 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, from a seismic acquisition system, seismic data regarding a subsurface region of interest, wherein the seismic data is a function of a plurality of coordinates. The method also includes obtaining, using a seismic processing system, a velocity model regarding the subsurface region of interest; forming a plurality of dip-angle gathers based, at least in part, on the seismic data and the velocity model; and  generating, using the seismic processing system, a filtering mask, wherein the filtering mask is a function of the plurality of coordinates. The method further includes, generating, for each dip-angle gather, a transformed gather by applying a transformation to the dip-angle gather; updating the filtering mask based, at least in part, on an energy extremum of the transformed gather; and forming a filtered diffraction image based, at least in part, on the plurality of dip-angle gathers and the updated filtering mask. The method still further includes determining, using a seismic interpretation system, a drilling target in the subsurface region based on the filtered diffraction image.

[0005] In general, in one aspect, embodiments disclosed herein relate to a system. The system includes a seismic acquisition system, a seismic processing system, and a seismic interpretation system. The seismic acquisition system is configured to record seismic data regarding a subsurface region of interest, wherein the seismic data is a function of a plurality of coordinates. The seismic processing system is configured to receive the seismic data; to obtain a velocity model regarding the subsurface region of interest; to form a plurality of dip-angle gathers based, at least in part, on the seismic data and the velocity model; and to generate a filtering mask, wherein the filtering mask is a function of the plurality of coordinates. The seismic processing system is also configured to generate, for each dip-angle gather, a transformed gather by applying a transformation to the dip-angle gather; update the filtering mask based, at least in part, on an energy extremum of the transformed gather; and form a filtered diffraction image based, at least in part, on the plurality of dip-angle gathers and the updated filtering mask. The seismic interpretation system is configured to determine a drilling target in the subsurface region based on the filtered diffraction image.

[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 illustrates a velocity model for a simplified subsurface structure in accordance with one or more embodiments.

[0011] FIG. 4 illustrates a zero-offset section for a simplified subsurface structure in accordance with one or more embodiments.

[0012] FIG. 5 shows an example of a dip-angle gather in accordance with one or more embodiments.

[0013] FIG. 6 shows an example of a diffraction image in accordance with one or more embodiments.

[0014] FIG. 7 shows a drilling system in accordance with one or more embodiments.

[0015] FIG. 8 shows a flowchart in accordance with one or more embodiments.

[0016] FIG. 9 shows a flowchart in accordance with one or more embodiments.

[0017] FIG. 10 illustrates a filtering mask in accordance with one or more embodiments.

[0018] FIG. 11 illustrates a diffraction image in accordance with one or more embodiments.

[0019] FIG. 12 illustrates a filtered diffraction image in accordance with one or more embodiments.

[0020] FIG. 13 illustrates a velocity model, a diffraction image and a filtered diffraction image, in accordance with one or more embodiments.

[0021] FIG. 14 illustrates a block diagram of a computer system in accordance with one or more embodiments.DETAILED DESCRIPTION

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

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

[0024] In the following description of FIGs. 1-14, 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.

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

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

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

[0028] In general, disclosed embodiments include systems and methods to filter diffraction images by correcting phase reversal using flat diffractor gathers. In particular, some embodiments generate diffraction images by migration techniques that retain diffraction events and exclude reflection events. However, when encountering a diffractor, seismic waves undergo a 180-degree phase shift (also called phase reversal) that may significantly affect the accuracy of the diffractors location. Some embodiments implement techniques to mitigate the effects of phase reversal based on the energy of gathers of flat diffractors. For example, dip-angle gathers may be locally transformed to generate corresponding flat diffractor gathers. Each flat diffractor may be affected differently by the phase reversal of the seismic waves. A filtering mask may be implemented to identify and retain only high-energy regions in flat diffractor gathers. Effects of phase reversal on the diffraction image may be reduced by combining the filtering mask with the diffraction image to generate a filtered diffraction image.

[0029] The resulting filtered diffraction image may then be used for further seismic data interpretation, such as in updating the spatial extension 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.

[0030] 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) .

[0031] 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” . Diffracted waves may also occur, when radiated waves are backscattered by an obstacle or a geological discontinuity (112) whose size is of the same order of magnitude as the wavelength of the radiated waves. 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 be also boundaries between faults, fractures, or groups of fractures within a rock. The geological discontinuities (112) may delineate a hydrocarbon reservoir (104) .

[0032] At the surface, refracted seismic waves (110) , reflected seismic waves (114) and diffracted waves 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) .

[0033] 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 delimits the time sample at which the seismic acquisition system (100) activated the seismic source (106) and acquired the seismic data by the seismic receivers (116) .

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

[0035] Seismic data may also refer to data acquired over different time intervals, such as, for example, in cases where seismic surveys are repeated to obtain time-lapse data. Seismic data may also be pre-processed data, e.g., arranging data in a “common shot gather” (CSG) domain, i.e., sorting waveforms as acquired by different receivers and having a single source location. The noted seismic data is not intended as limiting, and any other suitable seismic data is intended to fall within the scope of the present disclosure.

[0036] 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) illustrates the detection of 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 several seismic receivers (116) . 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” .

[0037] In the CSG (204) shown in FIG. 2 the vertical axis indicates the time (206) and the horizontal axis indicates the offset (208) . In some embodiments, direct seismic waves (122) , refracted seismic waves (110) , and reflected seismic waves (114) may be located in the CSG (204) by their arrival times, i.e., the times instants 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” .

[0038] 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 in a form that approximates a straight line.

[0039] 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) .

[0040] Seismic data (202) may be processed by a seismic processing system (220) . Processing seismic data (202) may consist of several key groups of functions, each serving a specific purpose in the processing workflow. For example, the steps may include data injection, the loading, sorting and arrangement of raw seismic data (202) , acquired from various sources such as seismographs or land-based sensors, into the processing system. This data may include seismic waveforms, and may also include well logs, and survey information.

[0041] Data quality control is critical in seismic processing. The seismic processing system (220) employs various tools and techniques to identify and correct any artifacts, noise, or errors in the data. This step ensures the accuracy and reliability of subsequent processing steps.

[0042] Further, the raw seismic data (202) may be “conditioned” , i.e., the raw seismic data (202) is pre-processed to enhance its quality and make it suitable for further analysis. This  step may include procedures such as filtering, deconvolution, noise suppression, and signal enhancement.

[0043] In addition, data may be “stacked” . Stacking involves combining multiple seismic traces to improve data quality and increase signal-to-noise ratio. This may enhance the identification of subsurface features and reduces random noise interference.

[0044] Velocity analysis is crucial for accurate imaging and interpretation of subsurface structures. It involves estimating the time-depth relationship of seismic reflections and determining the velocity model of the subsurface. The seismic processing system (220) will provide methods for multiple methods of performing velocity analysis, including normal moveout analysis, iterative Kirchhoff time-and depth-migration, tomography, and full waveform inversion.

[0045] Migration is a key step that transforms the processed seismic data (202) from the time domain to the depth domain, providing a more accurate representation of subsurface structures. It helps in locating and positioning geological features accurately.

[0046] The seismic processing system (220) may provide visualization tools to render the seismic data (202) in a visual format, enabling geoscientists to analyze, interpret, and perform visual quality control more effectively. This can include 2D / 3D seismic displays, depth slices, horizon maps, and virtual reality visualization.

[0047] The final step involves generating reports and documenting the results of the seismic processing workflow. This includes recording the processing parameters, interpretation results, and any uncertainties or limitations associated with the data processing. The seismic processing system (220) is required to perform these groups of steps for even a small commercial seismic survey.

[0048] The seismic processing system (220) may consist of various hardware components that work together to process and analyze seismic data (202) . Seismic processing requires significant computational power and storage capacity. High-performance servers and workstations are used to handle the massive amount of data and perform complex processing algorithms efficiently. Seismic data (202) can be massive, reaching terabytes  or even petabytes in size. Reliable and high-capacity storage systems, such as Network Attached Storage (NAS) or Storage Area Networks (SAN) , are utilized to store and manage the seismic data (202) effectively. In some cases, where processing demands are extremely high, the seismic processing system (220) may utilize cluster systems. Clusters are groups of interconnected computers or servers that work together to distribute the processing workload, enabling parallel processing and faster data analysis. A robust and high-speed network infrastructure is vital for seamless data transfer between different components of the seismic processing system (220) . This ensures efficient communication and data sharing, especially in multi-node or distributed processing environments.

[0049] The seismic processing system (220) may use GPUs for accelerating the computation of seismic processing algorithms. Their parallel processing capabilities significantly speed up tasks such as migration, inversion, and visualization. Despite advances in storage technology, data on tapes is still often used for long-term archiving and backup purposes. Tape systems provide high-capacity, cost-effective, and reliable storage solutions for seismic data (202) . Various peripherals such as monitors, keyboards, mice, network switches, uninterruptible power supply (UPS) , and backup power generators complete the hardware setup of a seismic processing system (220) . These peripherals ensure smooth operation, user interaction, and data integrity.

[0050] The software / firmware are at least as integral a part of the seismic processing system (220) as the hardware components and a seismic processing system (220) equipped with a unique software program is at least as distinctively different from other seismic processing systems without the unique software program as a seismic processing system (220) with GPUs is different from one without GPUs.

[0051] Seismic data (202) may be processed by a seismic processing system (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. Pressure waves are often referred to as “primary-waves” or “P-waves” . Shear waves are often referred to “secondary waves” or “S-waves” . At any given location, shear waves propagate at a different, slower, velocity than pressure 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 may be 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.

[0052] In some embodiments, seismic data (202) is processed by a seismic processing system (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 “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.

[0053] In one or more embodiments, the seismic image (230) may be a diffraction image. A seismic image (230) may include reflection events (reflectors) and diffraction events (diffractors) . Because reflectors are usually associated to more energetic events than diffractors, the generation of diffraction images, containing only diffractors, may facilitate differentiating between reflection events and diffraction events. Diffraction images provide useful information to reveal fault surfaces, fracture zones and erosional surfaces. Furthermore, diffraction images may provide details of the subsurface region beyond the resolution of images containing reflections.

[0054] Turning to FIGs. 3 and 4. FIG. 3 shows an example of a simplified velocity model, and FIG. 4 illustrates the zero-offset section of the corresponding seismic image. The  velocity model of FIG. 3 is composed of a region (302) of velocity v2 and a smaller region (304) of velocity v1. The vertical axis indicates a vertical spatial coordinate (306) , or depth, and the horizontal axis indicates a horizontal spatial coordinate (308) . The boundary of the smaller region (304) may be considered as a reflector. However, the corner (310) of the smaller region causes the spread of incident seismic waves, and thus, the corner (310) of the smaller region (304) may be considered as a diffractor. In FIG. 3 the corner’s horizontal coordinate (312) is 500 and the sharp edge vertical coordinate (314) is also 500.

[0055] In FIG. 4 the vertical axis indicates the time coordinate (402) and the horizontal axis indicates a horizontal spatial coordinate (404) . At the left end of the zero-offset section a reflector (406) is visible in the form of a flat boundary, and a diffraction (408) appears as an arc traversing the complete zero-offset section in the direction of the horizontal coordinate (404) . FIG. 4 illustrates that the diffraction (408) , or diffraction event, undergoes phase reversal on either side of the apex (414) of the arc. More specifically, the amplitude of the diffraction (408) , for example at point (410) , changes from black (positive) to white (negative) along the time coordinate (402) . At point (412) , on the contrary, the amplitude of the diffraction (408) changes from white (negative) to black (positive) to. In this example, phase reversal takes place at the apex (414) of the diffraction (408) . The apex horizontal coordinate (416) provides the horizontal location of the corner (310) considered as the diffraction (408) . Even though the results of FIG. 4 are expected as diffracted waves exhibit phase reversal when backscattered, the phase reversal may lead to inaccuracies in the location of the corner (310) if not correctly accounted for in forming an image.

[0056] Generating diffraction images may include the use of dip-angle gathers having dip-angle in one axis and migrated depth in another axis. The dip angle is the angle between the vertical axis and the normal direction of a presumed reflector. At constant offset and with the correct velocity, reflections in dip-angle gathers may form concave shapes. In contrast, due to the multiangle scattering, diffractions may appear as flat events. Dip-angle gathers may thus facilitate the separation of reflectors and diffractors. Methods to generate dip-angle gathers may include techniques based on ray theory, wave-equation  migration, phase shift plus interpolation (PSPI) , or any other technique known to those skilled in the art without limiting the scope of the invention.

[0057] FIG. 5 shows an example of a dip-angle gather (500) for a horizontal coordinate of 500 corresponding to the velocity model of FIG. 3. The vertical axis in FIG. 5 indicates the vertical coordinate (502) , or depth, and the horizontal axis indicates the dip angle (504) in degrees. Phase reversal is clearly observed by the positive amplitude region (506) adjacent to the negative amplitude region (508) . As expected, the phase reversal takes place at points with vertical coordinate value (510) of approximately 500, corresponding to the location of the corner vertical coordinate (314) . If phase reversal is not corrected, processing the dip-angle gather to generate a diffraction image without reflectors may lead to inaccurate results.

[0058] FIG. 6 illustrates a diffraction image (600) corresponding to the velocity model of FIG. 3. The diffraction image (600) has been generated without correcting the phase reversal observed in FIGs. 4 and 5. The vertical axis indicates the vertical spatial coordinate (602) or depth, and the horizontal axis indicates the horizontal spatial coordinate (604) . FIG. 6 shows that the reflector has been effectively removed and is absent in the diffraction image (600) . However, even though only manifestations of the diffractor (606) , the corner (310) , remains in the diffraction image (600) the location of the diffractor, expected to have a diffractor vertical coordinate (608) of 500 and a diffractor horizontal coordinate (610) of 500, is not imaged with precision. The unfocused imaged diffractor (606) can be observed more clearly in the magnified portion (612) . The imaged diffraction event is spread around and not precisely located at the position of the corner (310) , indicated by the circle (614) . Thus, processing techniques that take into account the effect of phase reversal may assist in generated diffraction images with improved clarity and accuracy.

[0059] As illustrated in FIGs. 2, 4 and 6, processing of seismic data (202) may generate a seismic image (230) that may reveal the three-dimensional geometry of a subsurface region of interest (102) . In particular, the geological boundaries (236, 238) may delineate a hydrocarbon reservoir (104) . If a seismic image (230) indicates the potential presence  of hydrocarbons in the subsurface region of interest (102) , a drilling system may drill a wellbore (118) to confirm the presence of those hydrocarbons.

[0060] FIG. 7 shows a drilling system (700) in accordance with one or more embodiments. As shown in FIG. 7, a wellbore (118) following a wellbore trajectory (704) may be drilled by a drill bit (706) attached by a drillstring (708) to a drilling rig (710) located on the surface (124) of the earth. The drilling rig (710) may include framework, such as a derrick (714) to hold drilling machinery. A crown block (711) may be mounted at the top of the derrick (714) , and a traveling block (713) may hang down from the crown block (711) by means of a cable (715) or drilling line. One end of the cable (715) may be connected to a drawworks (not shown) , which is a reeling device that may be used to adjust the length of the cable (715) so that the traveling block (713) may move up or down the derrick (714) .

[0061] A top drive (716) provides clockwise torque via the drive shaft (718) to the drillstring (708) in order to drill the wellbore (118) . The drillstring (708) may include a plurality of sections of drillpipe attached at the uphole end to the drive shaft (718) and downhole to a bottomhole assembly ( “BHA” ) (720) . The BHA (720) 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 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 (720) and the drillstring (708) may include other drilling tools known in the art but not specifically shown.

[0062] The wellbore (118) may traverse a plurality of overburden (722) layers and one or more formations (724) to a hydrocarbon reservoir (104) within the subterranean region of interest (728) , and specifically to a drilling target (730) within the hydrocarbon reservoir (104) . The wellbore trajectory (704) may be a curved or a straight trajectory.  All or part of the wellbore trajectory (704) may be vertical, and some parts of the wellbore trajectory (704) may be deviated or have horizontal sections. One or more portions of the wellbore (118) may be cased with casing (732) in accordance with a wellbore plan.

[0063] To start drilling, or “spudding in” the well, the hoisting system lowers the drillstring (708) suspended from the derrick (714) 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 (716) to rotate the drillstring (708) through the drive shaft (718) . The weight of the drillstring (708) combined with the rotational motion enables the drill bit (706) to bore the wellbore (118) .

[0064] The drilling system (700) may be disposed at and communicate with other systems in the well environment, such as the seismic processing system (220) , a seismic interpretation system (740) , and a wellbore planning system (738) . The drilling system (700) 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 (700) 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 (700) , the well-measured data may include mud properties, flow rates, drill volume and penetration rates, rock physical properties, etc.

[0065] A seismic interpretation system (740) is primarily used by geoscientists, seismic interpreters, and exploration teams in the oil and gas industry for analyzing seismic data to understand subsurface geological structures. Seismic interpreters use the workstation to visualize seismic data, including 2D and 3D seismic volumes, cross-sections, time slices, and attribute maps. These visualizations provide insights into subsurface structures, faults, and potential hydrocarbon reservoirs.

[0066] Interpreters may pick and interpret key geological horizons within seismic data to identify stratigraphic layers, boundaries, and structural features. Horizon interpretation tools and workflows allow for the accurate extraction of geological information from seismic volumes. A seismic interpretation system (740) enable interpreters to identify and interpret subsurface faults that may impact hydrocarbon reservoirs. Fault interpretation  tools and visualization techniques help in understanding fault geometry, connectivity, and spatial relationships. Seismic attributes, such as amplitude, frequency, and gradient, provide additional information about subsurface properties and can be analyzed using various algorithms and statistical methods. Attribute analysis tools in the workstation aid in defining reservoir characteristics, identifying anomalies, and highlighting potential hydrocarbon traps.

[0067] Interpreters may use the seismic interpretation system (740) to build 3D geological models by integrating seismic data with well-log data, geological knowledge, and other geophysical information. These models help in estimating reservoir properties, optimizing well locations, and predicting hydrocarbon distribution. Interpreters may analyze and characterize hydrocarbon reservoirs by integrating different data sources, including seismic data, well logs, production data, and seismic inversion results. Workstations provide tools for reservoir property estimation, quantitative analysis, and reservoir performance evaluation.

[0068] The seismic interpretation system (740) may facilitate prospect generation and evaluation, where interpreters identify and assess areas with high hydrocarbon exploration potential. They can perform detailed geological and geophysical analysis, identify drilling targets, and quantify the risk and uncertainty associated with potential prospects. Finally, workstations enable interpreters to collaborate with team members, share interpretation results, and communicate findings effectively. Interpretation software allows for the creation of reports, annotated images, and presentations to communicate geological interpretations to stakeholders.

[0069] The seismic interpretation system (740) are essential tools for geoscientists involved in exploration and production activities, helping them make informed decisions about drilling locations, optimize production strategies, and understand complex subsurface geological structures. The seismic interpretation system (740) may be a specialized computer system used by geoscientists and seismic interpreters for analyzing and interpreting seismic data. The seismic interpretation system (740) may be implemented on a computing device such as that shown in FIG. 14.

[0070] Seismic interpretation involves intensive tasks like data visualization, horizon picking, attribute analysis, and 3D modeling. A high-performance seismic interpretation system (740) with a powerful processor, ample memory, and a high-resolution display is essential to handle these computationally demanding tasks efficiently. Dedicated GPUs may be crucial for real-time rendering of seismic data, enabling smooth and interactive visualization. GPUs with high memory and parallel processing capabilities accelerate tasks like volume rendering and horizon visualization.

[0071] Seismic interpretation often involves working with large and complex datasets. Multiple high-resolution monitors allow interpreters to view seismic data, cross-sections, time slices, attribute maps, and other visualizations simultaneously, enhancing productivity and analysis accuracy. The seismic interpretation system (740) may be equipped with industry-standard software applications tailored for seismic interpretation, such as seismic data processing and visualization tools, horizon and fault interpretation systems, attribute analysis software, and 3D modeling software.

[0072] Seismic interpretation projects generate substantial amounts of data, including seismic volumes, processed data, interpretation results, and velocity models. A high-capacity and fast storage system, such as solid-state drives (SSDs) or RAID arrays, is necessary to store and access this data efficiently. The seismic interpretation system (740) often requires network connectivity to access centralized data repositories, collaborate with colleagues, and share interpretation results. A robust network infrastructure with fast Ethernet or fiber connections ensures smooth data transfer and collaboration capabilities.

[0073] Essential peripherals like keyboards, mice, and graphics tablets enable efficient interaction with data and software interfaces. Additionally, color-calibrated and high-accuracy input devices enhance the precision of interpretation tasks like picking horizons or drawing geological features. The seismic interpretation system (740) should have backup solutions in place to protect valuable data from loss or damage. Automated backup systems, external storage devices, or network-attached storage (NAS) can be utilized to ensure data safety. In some cases, seismic interpreters may need remote access to the seismic interpretation system (740) or collaborate with colleagues remotely. Setting  up remote access capabilities, such as Virtual Private Networks (VPNs) or remote desktop solutions, allows interpreters to work from different locations and share their work effectively. The seismic interpretation system (740) may be customized to meet the needs of interpreters and the specific requirements of projects. The hardware specifications may vary based on factors like the complexity of interpretations, the size of datasets, and the software tools utilized.

[0074] In some embodiments, rock physical properties may be used by the seismic interpretation system (740) to help determine a location of a hydrocarbon reservoir (104) (or other subterranean features) . Knowledge of the existence and location of the hydrocarbon reservoir (104) and other subterranean features may be transferred from the seismic interpretation system (740) to a wellbore planning system (738) . The wellbore planning system (738) may use information regarding the hydrocarbon reservoir (104) location to plan a well, including a wellbore trajectory (704) 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 (704) 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.

[0075] 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 (704) may be transferred to the drilling system (700) described in FIG. 7. The drilling system (700) may drill the wellbore (118) along the planned wellbore trajectory (704) to access the drilling target (730) in the hydrocarbon reservoir (104) .

[0076] The wellbore planning system (738) is used in the oil and gas industry for designing and planning drilling operations. It assists drilling engineers and teams in making strategic decisions related to wellbore placement, casing design, trajectory planning, and well path  optimization. The wellbore planning system (738) allows drilling engineers to visualize and interact with wellbore data in a 3D environment. It provides a graphical representation of the planned well trajectory, existing well paths, geological formations, and potential hazards.

[0077] The wellbore planning system (738) integrates geological models, well logs, seismic data, and other subsurface information to facilitate the creation of accurate and realistic wellbore plans. By incorporating geological models, drilling engineers can optimize well placement in reservoir targets and avoid geohazards. Furthermore, the wellbore planning system (738) may assist in designing optimal well trajectories based on reservoir targets, geologic constraints, and drilling objectives. Engineers can define well paths that maximize drilling efficiency, reach specific targets (horizontal or vertical) , and account for geological formations and structural complexities.

[0078] The wellbore planning system (738) incorporates collision-avoidance algorithms to assess potential collision risks between nearby wells, salt bodies, or other subsurface infrastructure. By considering uncertainties in subsurface data and drilling conditions, the wellbore planning system (738) may assess collision probabilities for planned well paths. This analysis helps in quantifying risks associated with collision potential and improving well placement decisions. The wellbore planning system (738) provides real-time alerts to prevent wellbore collisions and maintain drilling safety.

[0079] The wellbore planning system (738) helps drilling engineers in designing casing strings and selecting appropriate tubulars based on the wellbore conditions, planned drilling operations, and regulatory requirements. It considers factors such as pressure, temperature, well depth, formation properties, and casing load capacity. Furthermore, the wellbore planning system (738) performs torque and drag analysis to evaluate the forces and stresses acting on the drillstring during drilling operations. This analysis helps in identifying potential issues such as differential sticking, buckling, or limitations in the drilling equipment. The wellbore planning system (738) may have the capability to integrate real-time drilling data, such as downhole measurements, drilling parameters, and formation evaluation results. This integration allows engineers to monitor the drilling  progress, make on-the-fly adjustments to the well plan, and optimize drilling efficiency. Furthermore, the wellbore planning system (738) provides tools for generating reports, exporting data, and documenting drilling plans and decisions. These reports can be shared with regulatory agencies, drilling contractors, and other stakeholders to ensure alignment and compliance throughout the drilling lifecycle.

[0080] The wellbore planning system (738) assists drilling engineers in designing optimal well trajectories, minimizing risks, and maximizing drilling efficiency. They integrate various subsurface data sources, perform complex analyses, and provide visualization tools to support informed decision-making in well planning and drilling operations.

[0081] Turning to FIG. 8, FIG. 8 shows a flowchart in accordance with one or more embodiments. Specifically, FIG. 8 describes a general method to filter a diffraction image by generating a filtering mask from flat diffractor gathers. While the various blocks in FIG. 8 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.

[0082] In Block 800, seismic data (202) associated with a subsurface region of interest is obtained, in accordance with one or more embodiments. The seismic data (202) is a function of a plurality of coordinates. 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. Thus, the plurality of coordinates may be time and space coordinates. The seismic data (202) may be acquired using a seismic acquisition system (100) above a subsurface region of interest (102) . 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.

[0083] In Block 810, a velocity model regarding the subsurface region of interest 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.

[0084] An accurate velocity model (219) may be useful for accurate simulations and forward modeling of seismic propagation in the subsurface. A velocity model (219) 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 (219) may be considered an inverse problem, where the applied process must determine the subsurface velocity model (219) that resulted in the recorded seismic data (202) . The various processes and techniques used to process seismic data (202) to form a velocity model (219) 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 (219) from given seismic data (202) .

[0085] In Block 820, a plurality of dip-angle gathers is formed based, at least in part, on the seismic data and the velocity model, in accordance with one or more embodiments. Various types of algorithms may be used in the generation of dip-angle gathers. For example, some techniques use migration based on the ray theory, others use Kirchhoff depth migration with plane wave destructors. In some embodiments, diffractors in dip-angle gathers are imaged using wave-equation-based migration.

[0086] Given the propagation angles of source and receiver rays (denoted by the angles θs and θr, respectively) , the dip angle may be defined as: β=(θr-θs)  / 2   (1)

[0087] where β is the angle between the vertical axis and the normal direction of a presumed reflector. The angles θs and θr are measured from the vertical direction. Dip-angle gathers may be generated by binning and accumulating migrated images according to their dip angles. At the correct depth, specular reflections focus at a particular angle, the dip angle β of the reflector that generates the reflections. In contrast, diffractions, due to the  multiangle scattering, appear as flat events in dip-angle gathers. Summing over the dip-angle variable β produces the conventional migrated images, that include reflections and diffractions.

[0088] An example of a method for generating dip angle gathers is illustrated by the flowchart of FIG. 9, in accordance with one or more embodiments. The remaining steps of the method of FIG. 8 are presented below, after the discussion of FIG. 9. Turning to FIG. 9, in Block 910, a source wavefield and a receiver wavefield may be extrapolated based, at least in part, on the seismic data and the velocity model. A one-way equation may be used to extrapolate the two wavefields from the surface to all depths in the subsurface. The down-going one-way wave equation may be used to extrapolate the source wavefield, and the up-going one-way wave equation may be used to extrapolate the receiver wavefield. In some embodiments, a phase shift plus interpolation (PSPI) method may be used to extrapolate the wavefields.

[0089] In other embodiments, the receiver wavefield and the source wavefield may be extrapolated using reverse time migration (RTM) with the two-way wave equation. In RTM, a source wavefield may be obtained by forward modelling the propagation of a synthetic source function using the velocity model. A receiver wavefield may be generated using the same velocity model 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.

[0090] In some embodiments, dip-angle gathers may be generated by decomposing the incident and scattered waves into local plane waves at every imaging point. In Block 920, a plurality of local plane-wave receiver wavefields may be determined from the receiver wavefield, in accordance with one or more embodiments. The receiver wavefield may be transformed into a transformed receiver wavefield in a ray-parameter, p-domain. Each local plane-wave receiver wavefield Rd (x, y, z, ω; p) is parametrized by a value of the ray-parameter p. The plurality of local plane-wave receiver wavefields may be generated by decomposing an extrapolated receiver wavefield R (x, y, z, ω) into local plane-wave  components. According to one or more embodiments, the local plane-wave receiver wavefield may be generated by applying a windowed Radon transform to the extrapolated receiver wavefield R (x, y, z, ω) as follows:

[0091] where δx is the half window width and x′ represents the shift towards the center of the window. Even though in Equation (2) computations are performed with windows in the x-direction, Equation (2) may be adapted to use windows in other horizontal directions. The parameter p is the plane slope given by

[0092] where v is the local interval velocity. Furthermore, in one or more embodiments, the plane slope may have two horizontal components associated with an azimuthal angle  as follows

[0093] and then the extrapolated receiver wavefield R (x, y, z, ω) may be decomposed into local plane-wave components Rd (x, y, z, ω; px, py) with a 2D Radon Transform.

[0094] In other embodiments, the implementation of a recursive Radon transform may reduce the computational time spent to generate each local plane-wave receiver wavefield. A recursive Radon transform may be obtained by writing Equation (2) in terms of auxiliary variables as follows: Rd (x, y, z, ω; p) =exp (-iωpx) F (x, y, z, ω; p)    (5) C (x′, y, z, ω; p) =R (x′, y, z, ω) exp [-iωpx′]       (7)

[0095] where C and F are the auxiliary variables. Equations (4) - (6) are equivalent to Equation (2) and allow a recursive implementation to compute the Radon transform. Let xj denote the center point of a j-th window containing 2δx+1 points. The Radon transform at the next, (j+1) -th window, may be evaluated using Equation (4) : Rd (xj+1, y, z, ω; p) =exp (-iωpxj+1) F (xj+1, y, z, ω; p) . (8)

[0096] The variable F in Equation (5) may be obtained for the (j+1) -th window using the sum for the j-th window as follows: F (xj+1, y, z, ω; p) =F (xj, y, z, ω; p) -C (xj-δx, y, z, ω; p) +C (xj+δx+1, y, z, ω; p)     (9)

[0097] Thus, the evaluation of F for the (j+1) -th window may performed with only two additional operations, the subtraction of C at the point xj-δx, and the addition of the C at the point xj-δx+1. Furthermore, the phase shifting factor exp (-iωpxj+1) may be multiplied in Equation (7) after performing the subtracting and addition operations. The resulting recursive algorithm may reduce the computational complexity of evaluating the Radon Transform from O (M x N2) to O (M x N) , by avoiding the calculation of shared points between adjacent windows, where O denotes Landau’s symbol describing the asymptotic behavior of a function.

[0098] In some embodiments, the Radon Transform may be applied to a given source wavefield S (x, y, z, ω) to generate local plane-wave source wavefields to be combined with the local plane-wave receiver wavefields. Alternatively, only one incidence angle θs may be associated to the source wavefield S (x, y, z, ω) at each imaging point (x, y, z) . In some embodiments, the incidence angle θs of the source wavefield S (x, y, z, ω) at every point may be determined using a structure tensor method, as shown in Block 930 of FIG. 9. The structure tensor method may be used to determine directional changes in local neighborhoods, and / or to distinguish the type of directionality in local neighborhoods, such as, for example, a unique orientation, no local orientation, or equal orientation in all directions. The structure tensor may be computed, for example, as a linear convolution of a smoothing filter and gradients, where partial derivatives in gradients are approximated by discrete derivative operators.

[0099] The incidence angle θs of the source wavefield S (x, y, z, ω) may be computed for a predetermined frequency, that may be chosen from the high frequency band of the data. For example, the predetermined frequency may be chosen to be 50 Hz. The incidence angle θs may be associated to all frequencies of the source wavefield S (x, y, z, ω) . The source wavefield S (x, y, z, ω) may be generated by extrapolating the source function at the predetermined frequency from the source location to all points (x, y, z) .

[0100] In Block 940, each dip-angle gather may be generated by applying an imaging condition to a local plane-wave receiver wavefield and the source wavefield, in accordance with one or more embodiments. The 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. By transforming the source wavefield S (x, y, z, ω) and the receiver wavefield R (x, y, z, ω) to the frequency domain, the dip-angle gather may be given by the following equation:

[0101] Other frequency-domain migration methods known in the art may be used to generate the dip-angle gathers, such as the split-step Fourier, the Fourier finite-differences, implicit ω-x finite-differences, phase shift, etc.

[0102] Keeping with FIG. 9, in Block 950, amplitude balancing may be applied to the dip-angle gathers, in accordance with one or more embodiments. Amplitude balancing may be used to bring up the intensity of the diffraction energy to a level similar to the intensity of the reflection energy. Reflection energy usually has stronger intensity than diffraction energy, thus amplitude balancing may enhance the imaging of diffractors. Amplitude balancing may be applied using automatic gain control, logarithmic-based operations, statistical indicators, or any other technique known in the art.

[0103] In Block 960, a plurality of dip angle bins may be generated, and the dip-angle gathers may be replaced with the dip-angle bins, in accordance with one or more embodiments. Each dip-angle bin may be a combination of a set of dip angle gathers. The  combination may include one or more operations, such as, for example, binning and accumulating. Binning may facilitate the separation of reflection energy from diffraction energy. Returning to the method of FIG. 8, once the plurality of dip-angle gathers is formed, a filtering mask may be used to enhance the resolution of the diffraction image. Thus, as shown in Block 830, a filtering mask is generated, wherein the filtering mask is a function of the plurality of coordinates, in accordance with one or more embodiments. For example, the filtering mask M may be a function of spatial coordinates (x, y, z) . In an initial step, all entries of the filtering mask M may be set equal to zero. Then the plurality of dip-angle gathers may be processed to identify the regions of diffraction energy and exclude changes in phase that may lead to inaccuracies in the diffraction image, as those shown in FIG. 6 discussed above. The regions retained after processing the dip-angle gather may then be used in the construction of the filtering mask M.

[0104] In Block 840 of FIG. 8, a transformed gather is generated for each dip-angle gather, by applying a transformation to the dip-angle gather, in accordance with one or more embodiments. The transformation may be applied to each single dip-angle gather, or to dip-angle bins if binning has been performed. In some embodiments, the linear Radon Transform may be used to transform the dip-angle gather Ig (x, y, z, β) to a transformed gather It (x, y, z′, sp) , where z=spβ+z′, and sp is a slope parameter. Events of diffraction that appear as flat events in the dip-angle gather may be associated to events with zero slope parameter sp in the transformed gather.

[0105] In Block 850, the filtering mask M is updated based, at least in part, on an energy extremum of the transformed gather It (x, y, z′, sp) , in accordance with one or more embodiments. Each transformed gather It (x, y, z′, sp) may be scanned to identify the most energetic regions, some of which may correspond to diffraction events. An energy extremum may be for example, a maximum energy. One or more positions of the energy extremum may be then determined, as shown in Block 852. Each position may be given in terms of one or more coordinate values (xm, ym, z′m) . For each position (xm, ym, z′m) , a diffraction event may be identified if the absolute value of slope parameter spm is smaller than a predetermined value. The predetermined value may be for example, 2  degrees, corresponding to events that are almost flat. If spm is smaller than the predetermined value, a non-zero value may be assigned at a position of the filtering mask M corresponding to the coordinate values (xm, ym, zm) , as shown in Block 854. In some embodiments, the diffraction event is identified if neighboring slope parameters sp are smaller than the predetermined value. Neighboring slope parameters sp may be located for example, at coordinates (xm, ym, z′m) but with 1 or 2 grid points adjacent to the slope parameter spm. The filtering mask M may be progressively updated as diffraction events and their positions are identified. The filtering mask M may then be used as a filter that retains regions of high energy corresponding to the identified diffraction events.

[0106] Continuing with the method of FIG. 8, in Block 856 a smoothing function may be applied to the updated filtering mask M, in accordance with one or more embodiments. Examples of smoothing functions include a Gaussian function, or average smoothing, or any other smoothing technique known in the art. FIG. 10 shows an example of an updated filtering mask M (1000) , generated for the example of velocity model of FIG. 3 and smoothed with a Gaussian function, in accordance with one or more embodiments. The vertical axis (1002) indicates depth, and the horizontal axis (1004) indicates a horizontal distance. Dark regions of the updated filtering mask M, as those indicated by the circles (1006) , correspond to non-zero valued regions that will be retained when the updated filtering mask is applied.

[0107] In Block 860 of FIG. 8, a filtered diffraction image If (x, y, z) is formed based, at least in part, on the plurality of dip-angle gathers and the updated filtering mask, in accordance with one or more embodiments. Forming the filtered diffraction image If (x, y, z) may include applying a windowed median filter to each of the plurality of dip-angle gathers, as shown in Block 857. The windowed median filter may be applied over the dip angle variable, β, in order to remove reflection energy while preserving the diffraction energy. Because of its flatness, the diffraction energy in the dip-angle gather domain is not removed by the windowed median filter. In some embodiments, the windowed median filter is applied to the absolute value of each of the plurality of dip-angle gathers.

[0108] Further, an initial diffraction image may be generated by merging the plurality of the dip-angle gathers, as shown in Block 858. Reflections and diffractions may have been already located onto to their generators with the migration algorithm used to generate the dip-angle gathers. A diffraction image is directed to locate diffractor-only events, thus reflections may need to be reduced or minimized. A proper domain and gathering may facilitate the reduction of reflections. For example, continuous reflections may be removed by local slant stack or poststack sections. Examples of methods to detect diffractors include signal semblance along diffraction trajectories, focusing-defocusing in prestack shot gathers, binning migrated images according to angles, and plane-wave slope sections, among others. In some embodiments, the initial diffraction image may be obtained by merging, or stacking, the dip-angle gathers. Each dip-angle gather may have been first filtered with the windowed median filter of Block 857. The merging may be expressed, for example, by the following equation:

[0109] where Mβ denotes the median filter operator applied along the dip angle variable β within a window [mi-δβ, mi+δβ] that is centered at mi and with length 2δβ. The parameter nwin denotes the number of windows and ε is a stabilization term.

[0110] The initial diffraction image Id (x, y, z) may then be combined with the updated filtering mask M to generate the filtered diffraction image If (x, y, z) . For example, the filtered diffraction image If (x, y, z) may be determined according to the following product: If (x, y, z) =Id (x, y, z) ·M (x, y, z)     (12)

[0111] The product of initial diffraction image Id (x, y, z) and the updated filtering mask M (x, y, z) may attenuate inaccuracies of the initial diffraction image Id (x, y, z) due to reflections and phase reversal of diffractors, leading to a more accurate location of the diffractors.

[0112] Keeping with FIG. 8, In Block 870, a drilling target in the subsurface region may be determined based on the filtered diffraction image If (x, y, z) , in accordance with one or more embodiments. A drilling target (730) in a wellbore (118) may be based on, for example, an expected presence of gas or another hydrocarbon. Locations in a filtered diffraction image If (x, y, z) may indicate an elevated probability of the presence of a hydrocarbon and may be targeted by well designers. On the other hand, locations in filtered diffraction image If (x, y, z) indicating a low probability of the presence of a hydrocarbon may be avoided by well designers.

[0113] In Block 880, a wellbore trajectory to intersect the drilling target in the subsurface region is planned using the filtered diffraction image If (x, y, z) , in accordance with one or more embodiments. Knowledge of the location of the drilling target (730) and the filtered diffraction image If (x, y, z) may be transferred to a wellbore planning system (738) . Instructions associated with the wellbore planning system (738) may be stored, for example, in the memory (1409) within the computer system (1400) described in FIG. 14 below. The wellbore planning system (738) may use the knowledge of the location of the drilling target (730) and of the filtered diffraction image If (x, y, z) to plan a wellbore trajectory (704) within the subterranean region of interest (728) .

[0114] In Block 890, a wellbore is drilled guided by the planned wellbore trajectory, in accordance with one or more embodiments. The wellbore planning system (738) may transfer the planned wellbore trajectory (704) to the drilling system (700) described in FIG. 7. The drilling system (700) may drill the wellbore (118) along the planned wellbore trajectory (704) to access and produce the hydrocarbon reservoir (104) to the surface (124) .

[0115] The proposed method is illustrated using synthetic examples as is common in the art and is advantageous at least because the solutions are known a priori. The example starts with a numerical representation of a simplified portion of the earth (the "model" ) . Then synthetic data is generated from the numerical solution of the differential equation describing the physics of seismic wave propagation. Embodiments of the disclosed  invention are applied to the synthetic data and the results are compared with the model to assess the performance of the disclosed embodiments.

[0116] FIG. 11 shows an example of an initial diffraction image (1100) corresponding to the velocity model of FIG. 3, according to one or more embodiments. In the initial diffraction image (1100) of FIG. 11 the vertical axis indicates the depth-dimension (1102) , and the horizontal axis indicates one horizontal dimension (1104) . A windowed median filter was applied to the absolute value of each of the dip-angle gathers used to generate the initial diffraction image (1100) , with a window length of 60 grid points. As seen in FIG. 11 diffractors (1106) are visible in the initial diffraction image (1100) , however, even though their location appears in the expected region, is not imaged with precision.

[0117] On the other hand, FIG. 12 shows an example of a filtered diffraction image (1200) corresponding to the velocity model of FIG. 3, according to one or more embodiments. The filtered diffraction image (1200) of this example has been obtained by combining the initial diffraction image (1100) of FIG. 11 and the updated filtering mask (1000) of FIG. 10.In the filtered diffraction image (1200) the vertical axis indicates the depth-dimension (1202) , and the horizontal axis indicates one horizontal dimension (1204) . It can be observed that in the filtered diffraction image (1200) , the diffractor inside the rectangle (1206) , corresponding to the corner (310) shown in FIG. 3, is located with a desirable high accuracy. To better display the edge image, the area inside the rectangle (1206) is zoomed-in and shown in rectangle (1208) . The circle inside rectangle (1208) indicates the theoretical location of the edge (1210) , which appears to match the location provided by the filtered diffraction image (1200) . Such good agreement between the theoretical location of the edge (1210) and its diffractor image demonstrates the effectiveness of the proposed method.

[0118] Another example that showcases the capabilities of the proposed method is shown in FIG. 13. FIG. 13 shows a velocity model (1302) that is used to perform numerical simulation of wave propagation. 75 common-shot gathers with 150 receivers are generated using the 2D acoustic finite-difference method. The source function is defined  in terms of a Ricker wavelet with a central frequency of 16 Hz. The maximum duration of the simulations is 2.5 s with a time-sampling interval of 2 ms. The vertical axis indicates the depth-dimension (1304) and the horizontal axis indicates one horizontal dimension (1306) . Dip angles are considered ranging from -70° to 70° with a sampling rate of 3.5°. The scalebar (1307) indicates the magnitude of the velocity, ranging from 1500 m / s (0.9 miles / s) to 4400 m / s (2.7 miles / s) .

[0119] FIG. 13 also shows the initial diffraction image (1308) generated with the 75 common-shot gathers obtained with the above described numerical simulations. The wave equation was used as the migration technique for the source wavefield and the receivers wavefield. Dip angles ranging from -70° to 70° with a sampling rate of 3.5° were considered in the generation of the initial diffraction image (1308) . It can be observed that the initial diffraction image (1308) is free of reflection energy and the diffractors, as those indicated by the circles (1310) , are imaged in the correct regions but not at their precise locations. Upon filtering the initial diffraction image (1308) with the proposed methods, the filtered diffraction image (1312) is obtained, also shown in FIG. 13. Comparing diffraction images (1308) and (1312) , it can be observed that in the filtered diffraction image (1312) the diffractors, as those indicated by the circles (1314) , are well focused and accurately located at the right positions.

[0120] In some embodiments the wellbore planning system (738) and the seismic processing system (220) may each be implemented within the context of a computer system. FIG. 14 is a block diagram of a computer system (1400) 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 (1400) 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 (1400) 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 (1400) , including digital data, visual, or audio information (or a combination of information) , or a GUI.

[0121] The computer (1400) 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 (1400) is communicably coupled with a network (1402) . In some implementations, one or more components of the computer (1400) may be configured to operate within environments, including cloud-computing-based, local, global, or other environment (or a combination of environments) .

[0122] At a high level, the computer (1400) 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 (1400) 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) .

[0123] The computer (1400) can receive requests over network (1402) from a client application (for example, executing on another computer (1400) ) 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 (1400) 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.

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

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

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

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

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

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

[0130] In some embodiments, the computer (1400) 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) .

[0131] 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, from a seismic acquisition system, seismic data regarding a subsurface region of interest, wherein the seismic data is a function of a plurality of coordinates;using a seismic processing system:obtaining a velocity model regarding the subsurface region of interest;forming a plurality of dip-angle gathers based, at least in part, on the seismic data and the velocity model,generating a filtering mask, wherein the filtering mask is a function of the plurality of coordinates,for each dip-angle gather:generating a transformed gather by applying a transformation to the dip-angle gather,updating the filtering mask based, at least in part, on an energy extremum of the transformed gather, andforming a filtered diffraction image based, at least in part, on the plurality of dip-angle gathers and the updated filtering mask; anddetermining, using a seismic interpretation system, a drilling target in the subsurface region based on the filtered diffraction image.2.The method of claim 1, further comprising:planning, using a wellbore planning system, a planned wellbore trajectory to intersect the drilling target in the subsurface region using the filtered diffraction image; anddrilling, using a drilling system, a wellbore guided by the planned wellbore trajectory.3.The method of claim 1, wherein forming the plurality of dip-angle gathers comprises:extrapolating a receiver wavefield and a source wavefield based, at least in part, on the seismic data and the velocity model;determining a plurality of local plane-wave receiver wavefields from the receiver wavefield, wherein each local plane-wave receiver wavefield is parametrized by a value of a ray-parameter; andgenerating each dip-angle gather by applying an imaging condition to a local plane-wave receiver wavefield and the source wavefield.4.The method of claim 3, wherein forming a plurality of dip-angle gathers further comprises determining an incidence angle of the source wavefield by applying a structure tensor.5.The method of claim 1, wherein forming a plurality of dip-angle gathers comprises:generating a plurality of dip-angle bins, wherein each dip-angle bin is a combination of a set of dip angle gathers; andreplacing the dip-angle gathers with the dip-angle bins.6.The method of claim 1, wherein generating a transformed gather comprises applying a Radon Transform.7.The method of claim 1, wherein updating the filtering mask comprises:determining a coordinate value of the energy extremum; andsetting a non-zero value at a position of the filtering mask corresponding to the coordinate value.8.The method of claim 1, wherein forming the filtered diffraction image comprises applying a windowed median filter to an absolute value of each dip-angle gather.9.The method of claim 1, wherein generating a dip-angle gather comprises amplitude balancing.10.The method of claim 9, wherein amplitude balancing is performed by automatic gain control.11.The method of claim 1, wherein forming a filtered diffraction image comprises generating an initial diffraction image by merging the plurality of dip-angle gathers.12.The method of claim 1, further comprising applying a smoothing function to the updated filtering mask.13.The method of claim 12, wherein the smoothing function comprises a Gaussian function.14.A system, comprising:a seismic acquisition system configured to record seismic data regarding a subsurface region of interest, wherein the seismic data is a function of a plurality of coordinates;a seismic processing system configured to receive the seismic data and to:obtain a velocity model regarding the subsurface region of interest,form a plurality of dip-angle gathers based, at least in part, on the seismic data and the velocity model,generate a filtering mask, wherein the filtering mask is a function of the plurality of coordinates,for each dip-angle gather:generate a transformed gather by applying a transformation to the dip-angle gather,update the filtering mask based, at least in part, on an energy extremum of the transformed gather, andform a filtered diffraction image based, at least in part, on the plurality of dip-angle gathers and the updated filtering mask; anda seismic interpretation system configured to determine a drilling target in the subsurface region based on the filtered diffraction image.15.The system of claim 14, further comprising:a wellbore planning system configured to plan a planned wellbore trajectory to intersect the drilling target in the subsurface region using the filtered diffraction image; anda drilling system configured to drill a wellbore guided by the planned wellbore trajectory.16.The system of claim 14, wherein the seismic processing system is further configured to:extrapolate a receiver wavefield and a source wavefield based, at least in part, on the seismic data and the velocity model;determine a plurality of local plane-wave receiver wavefields from the receiver wavefield, wherein each local plane-wave receiver wavefield is parametrized by a value of a ray-parameter; andgenerate each dip-angle gather by applying an imaging condition to a local plane-wave receiver wavefield and the source wavefield.17.The system of claim 16, wherein the seismic processing system is further configured to determine an incidence angle of the source wavefield by applying a structure tensor.18.The system of claim 14, wherein the seismic processing system is further configured to:determine a coordinate value of the energy extremum; andset a non-zero value at a position of the filtering mask corresponding to the coordinate value.19.The system of claim 14, wherein the seismic processing system is further configured to generate an initial diffraction image by merging the plurality of dip-angle gathers.20.The system of claim 14, wherein the seismic processing system is further configured to apply a windowed median filter to an absolute value of each dip-angle gather.

Images