GB2700026A - System and method for joint free-surface and internal multiple imaging - Google Patents

Abstract

A method 300 for imaging geological features in a subsurface includes receiving 302 seismic data associated with the subsurface. The seismic data includes primaries, first and second free-surface multiples, and internal multiples. A first image depth range is selected 304. A first image is generated 308 corresponding to the first image depth range. The first image is based at least on the first free-surface multiples associated with the first image depth range. A second image is generated 312 of the subsurface based jointly on second free-surface multiples and internal multiples. The step of generating a second image of the subsurface involves a downgoing reflection within the first image of the subsurface and may involve modelling the first free-surface multiples based on the first image and generating partial demultiple data by subtracting the modelled multiples from the seismic data. Figure 3 to accompany the abstract

Description

SYSTEM AND METHOD FOR JOINT FREE-SURFACE AND INTERNAL MULTIPLE

IMAGING

BACKGROUND

TECHNICAL FIELD

[0001] Embodiments of the subject matter disclosed herein generally relate to processing seismic data acquired over land or water and, more specifically, to generating an image of the subsurface based on jointly processed free-surface and internal multiples.

DISCUSSION OF THE BACKGROUND

[0002] Subsurface exploration (e.g., hydrocarbon exploration and development, ore exploration, etc.) uses waves (e.g., seismic waves or electromagnetic waves) to explore the structure of underground formations on land and/or at sea (i.e., formations or features in the subsurface). These waves are generated by a seismic source and their reflections and/or refractions are collected during a seismic acquisition campaign, either on land or on water. More specifically, energy generated by the seismic source propagates as seismic waves downward into a geological formation, and part of the energy is reflected and/or refracted (for this description, the term 'reflected' or 'reflection' should be considered to refer to an interaction in the subsurface, which may relate to a reflection, a refraction, a diffraction, etc.) back up to the surface. Characteristics of the reflected/refracted energy detected by seismic sensors are used to produce an image of the earth's reflectivity.

[0003] As schematically illustrated in FIG. 1, when a marine seismic acquisition system 100 is used (the same applies for a land seismic survey), waves 112 emitted by a seismic source 110 into water 105, at a known location, penetrate an unexplored formation 120, which is located below the water bottom 106, and are reflected at plural interfaces 122, 124, 126 that separate the formation's layers, as they have different layer impedances. Sensors 130 (only one shown for simplicity), which are distributed along a streamer 132 towed by a vessel 102 (for the case of a marine seismic survey), detect the reflected waves 140, 150. The detected waves 140, 150 include primary reflections such as wave 140, which travel directly from a formation interface to a sensor, and multiple reflections (called "multiples") such as wave 150, which additionally undergo at least one more downward reflection and one more upward reflection. In the case of wave 150 in FIG. 1, the downward reflection is at the water surface 104 and the upward reflection is at reflector 122 and for this reason this multiple is called a free-surface related multiple.

[0004] Other multiples exist as schematically illustrated in FIGs. 2A to 2D. FIG. 2A shows primaries 140 reflected at two different interfaces, labeled R1 and R2 in the figure. The R1 and R2 label horizons or interfaces 122 to 126 in the figures. The primaries P1 and P2 are recorded by the sensors (not shown), as illustrated in FIG. 2A. The primary wavefield 140 is reflected at locations 210-1 and 210-2 in the horizons R1 and R2, respectively. The primary wavefield 140 may subsequently undergo an upper reflection 212-1, as illustrated in FIG. 2B, at the free-surface 104, followed by a lower reflection 212-2 or 212-3, in the subsurface 202, resulting in free-surface multiple arrivals 260. The free-surface multiple arrivals 260 (also labelled as SM11 and SM12) are called in this document first order free-surface multiples.

[0005] The primary wavefield 140 may have first a deep reflection at point 210-2, as schematically illustrated in FIG. 2C, followed by a reflection at the water surface at point 212-3, followed by a shallow reflection at point 212-4, by horizon R1, resulting in the multiple 262 (also called SM21). A similar multiple 262 that experiences a second deep reflection 212-3 is also shown in the figure as SM22. While FIGs. 2B and 2C show first order multiples (following initial shallow and deep reflections, respectively), FIG. 2D shows an internal multiple 264. The internal multiple 264 experiences an internal reflection at point 214-1 instead of having a reflection at the water surface 104. The internal multiples 264 are also called IM1 and IM2 in this figure. Free-surface multiples 260, 262 and internal multiples 264 both involve an upper reflection (at the free surface 104 or horizon R1) and a lower reflection (horizon R1 or R2). While Figure 1 shows first-order multiples, the process may be repeated to produce second-and higher-order multiples.

[0006] Note that, as used herein, the term "formation" refers to any geophysical structure into which source energy is injected to perform seismic surveying, e.g., land or marine, such that a "formation" may include a water layer when the context is marine seismic surveying. A similar picture may be used to describe the rays for a land seismic survey, except that there is no reflection from the water surface as there is no water.

[0007] The recorded seismic waves (recorded at the sensors 130) include primary arrivals and various types of multiples, e.g., surface-related multiples and internal multiples of various orders. Multiple arrivals contaminate the primary image and in many processing projects, multiples are modelled and subtracted (often adaptively) from the recorded data. There are many approaches in the art for modeling or attenuating multiples. Some of the most known methods are free-surface convolution-driven methods [1-3], free-surface wave-equation methods [4-6], internal convolution-driven methods [7, 8], and internal wave-equation methods [9]. These methods typically remove the multiples prior to imaging.

[0008] Instead of taking this approach, others have tried to use the multiples in the recorded data (typically free-surface multiples) for imaging the subsurface. This approach was taken by the authors in [10-13]. Alternatively, multiple prediction [6] associated with a subsurface may be generated based on free-surface multiple imaging. While the internal multiple modelling approach of [9] involves an upper reflection above a splitting horizon, and a lower reflection below a splitting horizon, the modelling step of the least-squares full-wavefield migration approach of [14] generalizes the internal multiple generation mechanism so that a splitting horizon does not need defining.

[0009] However, as is known in the seismic field, each of the methods mentioned in this section can be ineffective under certain conditions/environments. Therefore, imaging the subsurface with contaminated data remains a subject of continuing research, with new opportunities and challenges occurring as data acquisition systems evolve.

SUMMARY

[0010] According to an embodiment, there is a method for imaging geological features in a subsurface. The method includes receiving seismic data do associated with the subsurface, wherein the seismic data do includes primaries, first and second free-surface multiples, and internal multiples, selecting a first image depth range, generating a first image Ii of the subsurface, corresponding to the first image depth range, the first image la being generated based at least on the first free-surface multiples associated with the first image depth range, and generating a second image /2 of the subsurface based jointly on second free-surface multiples and internal multiples. The step of generating a second image of the subsurface involves a downgoing reflection within the first image /./ of the subsurface.

[0011] According to another embodiment, there is a computing system for imaging geological features in a subsurface, and the computing system includes an interface configured to receive seismic data do associated with the subsurface, wherein the seismic data do includes primaries, first and second free-surface multiples, and internal multiples and a processor connected to the interface and configured to select a first image depth range, generate a first image of the subsurface, corresponding to the first image depth range, the first image li being generated based at least on the first free-surface multiples associated with the first image depth range, and generate a second image /2 of the subsurface based jointly on second free-surface multiples and internal multiples. The step of generating a second image of the subsurface involves a downgoing reflection within the first image // of the subsurface.

[0012] According to yet another embodiment, there is a non-transitory computer readable medium including computer executable instructions, wherein the instructions, when executed by a processor, implement the methods discussed above.

BRIEF DESCRIPTION OF THE DRAWINGS

[0013] The accompanying drawings, which are incorporated in and constitute a pad of the specification, illustrate one or more embodiments and, together with the description, explain these embodiments. In the drawings: [0014] FIG. 1 is a schematic diagram illustrating a marine seismic data acquisition system that acquires primaries and multiples; [0015] FIGs. 2A to 2D are schematic diagrams illustrating primary reflections, first order multiples following a shallow reflection, first order multiples following a deep reflection, and internal multiples, respectively; [0016] FIG. 3 is a flow chart of a method for joint free-surface and internal multiples imaging; [0017] FIG. 4 illustrates plural image depth ranges and corresponding upper-reflectivities for the joint free-surface and internal multiples imaging; [0018] FIGs. 5A, 5B and 5C schematically illustrate the selection of the free-surface multiples and internal multiples for each image depth range; [0019] FIGs. 6A to 6C schematically illustrate a wavefield propagation mechanism for wavefields resulting in free-surface multiple modelling, internal multiple modelling, and internal multiple modelling with muted input data, respectively; [0020] FIG. 7A illustrates free-surface multiple imaging of a region extending from the water bottom up to a given depth, according to an existing method, FIG. 7B illustrates the free-surface multiple imaging of the water bottom, FIG. 7C illustrates free-surface multiple imaging of a sub-region of the region of FIG. 7A using an existing method, FIG. 7D illustrates free-surface multiple imaging of the same sub-region as in FIG. 7C, but using another existing method, FIG. 7E illustrates the joint free-surface and interbed multiple imaging of the sub-region using the method of FIG. 3, and FIG. 7F illustrates the differences between the results in FIGs. 7D and 7E; [0021] FIG. 8A illustrates the primary migration of the seismic data, FIG. 8B illustrates a free-surface water bottom multiple model, FIG. 8C illustrates the migrated seismic data of FIG. 8A from which the free-surface water bottom multiples have been removed, FIG. 8D illustrates a free-surface gas multiple model, FIG. 8E illustrates the migrated seismic data of FIG. 8A from which the free-surface water bottom multiples and the free-surface gas multiples have been removed based on the second image of the method of FIG. 3, FIG. 8F illustrates a water bottom gas interbed multiple model obtained based on first and second images of the method of FIG. 3, and FIG. 8G illustrates the migrated seismic data of FIG. 8A from which the free-surface water bottom multiples, the free-surface gas multiples, and the water bottom gas interbed multiples have been removed; and [0022] FIG. 9 is a schematic diagram of a computing device configured to implement methods for processing data according to the embodiments discussed herein.

DETAILED DESCRIPTION

[0023] The following description of the embodiments refers to the accompanying drawings. The same reference numbers in different drawings identify the same or similar elements. The following detailed description does not limit the invention. Instead, the scope of the invention is defined by the appended claims. The following embodiments are discussed, for simplicity, with regard to a marine seismic data acquisition. However, similar embodiments and methods may be used for a land data acquisition system.

[0024] Reference throughout the specification to "one embodiment" or "an embodiment" means that a particular feature, structure or characteristic described in connection with an embodiment is included in at least one embodiment of the subject matter disclosed. Thus, the appearance of the phrases "in one embodiment' or "in an embodiment" in various places throughout the specification is not necessarily referring to the same embodiment. Further, the particular features, structures or characteristics may be combined in any suitable manner in one or more embodiments.

[0025] According to various embodiments, a process or method of cascaded imaging and multiple attenuation approach jointly using surface and internal multiples is now discussed based on FIG. 3. The method 300 begins at step 302 by receiving initial seismic data (or input data) do. The input data do includes first free-surface multiples 260, second free-surface multiples 262, and internal-multiples 264, which are recorded by receivers 130 following an excitation of the seismic source 110. Note that the first free-surface multiples 260, second free-surface multiples 262, and internal-multiples 264 have been defined above with regard to FIGs. 2B to 2D. The method then optionally selects in step 304 plural image depth ranges 410, as illustrated in FIG. 4, for constructing the images of the subsurface. The depth ranges may be fixed for the whole volume or may vary in space, for example following the water bottom shape. FIG. 4 shows a table 400 in which plural images 402 are associated with corresponding image depth ranges 410, i.e., each image of the subsurface is generated for a limited depth range 410. Table 400 also indicates the upper-reflectivity 412 for each image, i.e., the subsurface ranges that are used for calculating the multiple upper reflection. The image depth ranges 410 may be selected by the operator to form an overall continuous interval, e.g., 0 to 1600 m in this example. However, the operator may also select the image depth ranges 410 to have gaps, e.g., 0-200m, 250-400m, 450-800m, etc., or to overlap, e.g., 0-200m, 150-400m, 350m-800m, etc. [0026] The method generates (using any known method, a couple of which are discussed later) in step 308 a first image 11 of the subsurface, using the at least free-surface multiples dFM. This step uses as the upper reflecting surface the water-air interface 106 (i.e., the free-surface). While FIGs. 2B and 2C show free-surface multiples 260 and 262 associated with two horizons R1 and R2, it is noted that for step 308 only those receiver-side free-surface multiples associated with a first image depth range are imaged (i.e., only the horizons associated with the first image are considered), which is a subset of all the recorded free-surface multiples. This concept is schematically illustrated in FIGs. 5A to 5C. For this discussion, it is assumed that the receiver-side multiples associated with the first horizon RI are associated with the first image (index 1 in FIG. 4, which is associated with the first image depth range), called herein 11," and receiver-side multiples associated with the second horizon R2 are associated with a second image (index 2 in FIG. 4, which is associated with a second image depth range), called herein "12." FIG. 5A shows receiver-side free-surface multiples 510 and 511 associated with horizon RI, which are generated by source S and recorded at receiver R while FIG. 5B shows receiver-side, free-surface multiples 512 and 513 associated with horizon R2. Figure 5C shows internal multiples 514 and 515 associated with horizons RI and R2 and R3 and R4, respectively. As the first image 11 is associated with the first image depth range (e.g., 0-200 m in FIG. 4), only the receiver-side free-surface multiples 510 and 511 (see FIG. 5A) are linked to the first image 11, and not the receiver-side free-surface multiples 512 and 513, or internal multiples 514, or 515. In other words, step 308 models only free-surface multiples 510 and 511, and not free-surface multiples 512 and 513, assuming that only four free-surface multiples exist (see FIGs. 5A and 5B). In reality, there are a large number of multiples associated with each image depth range. The free-surface multiples associated with the first image range, imaged in step 308, are labeled "dFM" herein. This discussion relates to receiver-side free-surface multiple imaging. It should be noted that algorithms for source-side free-surface multiple imaging or joint source-and receiver-side free-surface multiple imaging may be used.

[0027] The method then models (using any known method, a couple of which are discussed later) in step 310 the free-surface-related multiples mo, corresponding to the first image 11 of the subsurface, and generates in step 311 the first partial-demultiple data di, (generically denoted as dn in this document, where n is an integer that describes the first, second, etc. image, in this case, n=1, as this step generates the first image 11 of the subsurface). The partial-demultiple data di is obtained by subtracting the modelled free-surface multiples mo associated with the first image 11, from the input data do, i.e., di = do -mo. However, as the method iterates through various data sets, it generates further images, i.e., n increases and corresponds to the image index 402 in FIG. 4. The partial-demultiple data dn corresponding to each image is called "partial" because only the multiples mo_i (e.g., 510 and 511 in FIG. 5A) corresponding to reflections of a previous image (11 if n =2) of the subsurface have been removed, and other multiples (e.g., 512 and 513 in FIG. 5B) reflecting from other depth ranges (e.g., R2 horizon) in the subsurface remain in the partial-demultiple data dn.

[0028] Although the example shown in FIG. 4 suggests that the first image 11 is related to a depth of 0 to 200m, the first image of the subsurface may, for example, relate only to the water bottom or to a shallow depth range, for example up to 500m depth. The maximum depth of the first image may relate to a depth just before the earliest causal crosstalk arrives, e.g. just before twice the waterbottom depth. Any other depth range may be used. The partial-demultiple data di generated in step 311, is then used in step 312 to generate the second image 12 of the subsurface, which is associated with the second image depth range shown in FIG. 4. The second image 12 of the subsurface is generated using combined surface-related multiples 512 and 513 and internal multiples 514, as schematically illustrated in FIGs. 5B and 5C. The internal multiples contributing to the second image 12 will have an upper reflection within the first image 11. Note that the receiver-side multiples 512 and 513 used for generating the second image 12 are related to the second image depth, which includes the second horizon R2, but not the first horizon R1. Also note that an internal multiple 515 that is not reflecting on the horizon R2 (associated with the second image range depth) is not considered for generating the second image 12. However, the internal multiple 515 will play a role when generating the third image 13 as horizons R3 and R4 are associated with the third image range depth in this hypothetical example. Thus, FIGs. 5B and 5C show that the multiples 512, 513, and 514 must be related to R2. The horizon R1 (associated with the first image 11) serves as the upper-reflectivity for imaging interbed reflection 514. More specifically, the first image 11 of the subsurface is used to define the upper-reflectivity for the internal multiple imaging for the second image 12 of the subsurface. In one application, the first image 11 of the subsurface may need to be reversed in polarity to relate to the upper reflection rather than a lower reflection.

[0029] Still with regard to step 312, the primaries have one reflection below the free-surface, the free-surface multiples 260 and 262 have two reflections below the free-surface, and the internal multiples 264 have at least two lower reflections and one upper reflection below the free-surface.

[0030] In step 314, the method verifies whether a desired depth for imaging the subsurface has been reached. If the result of this step is that the desired depth has not been reached, i.e., branch 316, the process models in step 317 free-surface multiples and internal multiples and returns to steps 311 and 312 for generating the next demultiple data dz, followed by generating the third image 13, based on the updated next demultiple data dz. This process continues until the desired depth for imaging is reached. When this happens, the process follows the 318 branch and may display 320 the final image IF of the subsurface. In one application, this step combines together all the previous images 11, 12, etc. for generating the final image IF. In another application, this step presents only the last generated image as the final image.

[0031] The final image IF generated in step 320 may be used for identifying a desired geological feature, for example, a fault that might indicate an oil and gas reservoir, or an interface between layers, body, reservoir, etc. The geological feature may also be a formation feature such as, for example, a formation top (e.g., a boundary of a layer of a geological region through which a well passes).

[0032] The resulting images 11, 12, etc. of the subsurface may be used to predict multiple models for subtraction from the subsurface. The multiple models may be jointly or individually subtracted from the input data do. Alternatively, the images 11, 12, etc. of the subsurface may be combined (e.g., summed), and a combined multiple model may be derived and used for generating the next set of partial demultiple data.

[0033] The reflectivity r of a plane/interface discussed above and used in the above noted multiples determination method, is understood herein to be an "image" of the surveyed subsurface. In this regard, an image of the subsurface is generally understood to be a representation of the reflectivity in the earth, defined in space (x-z for 2D, and x-y-z for 3D). This may also be referred to as a migration, the migrated image, or the reflectivity. Thus, these terms are used interchangeably herein. An image of the subsurface may be the result of a single-step or optimized migration.

[0034] In some embodiments, to form an image, an imaging condition may need to be applied to the data. In one application, the imaging condition is a mathematical function applied to extrapolated down-going and up-going wavefields to form the image. The most common imaging condition is the cross-correlation imaging condition. Other options are the deconvolution imaging conditions, a variety of which are known in the field, for example, smoothing imaging condition, 2D deconvolution imaging condition or multi-dimensional deconvolution imaging condition.

[0035] The method 300 of cascaded, joint free-surface and internal multiple imaging process may also be expressed using pseudo-code language, which may run on a computing system, which is discussed later. More specifically, the system is configured to initialize a depth Zspiit, corresponding to image depth ranges 410 in FIG. 4. Thus, the depth Zspiit has initially a value zso = 0, corresponding to the free surface 104. As the method 300 iterates through steps 310 to 312, the value of Zspiit is increased to describe each image depth range 410. This means that this value limits the seismic data used for generating each image In, where n is an integer starting at 1. For the initial Zspiit value of 0, the input data includes at least the free-surface multiples.

[0036] The computing system then calculates, for example, using a least-squares approach, the imaging to depth zsn, using partial demultiple data dm which results in the image In (also known as rn, i.e., the reflectivities for the first image depth range). For a given n larger than 1, the depth range is zsn_i to zsn. The demultiple data dn includes combined free-surface and internal multiples. As noted above, for the first image (or first iteration) the method performs at least free-surface multiple imaging while for the subsequent iterations, the reflectivity from the previous iteration (rn_i) may be used to define the upper reflectivity (e.g., flip polarity of rn_i). The polarity may be flipped for the upper reflectivity as the seismic signal is reflected from below.

[0037] The computing system then models (step 310) the free-surface and internal multiples corresponding to the various reflectivities rn, resulting in the modeled multiples mn. The computing system subtracts in step 311, the multiples m" from partial demultiple data dn, to create partial demultiple data dn./ for the next iteration. In one application, the subtraction is adaptive. At this step, it is possible to program the computing system to additionally remove the primaries corresponding to the current image, from the partial demultiple data. The computing system then returns to step 310 if the desired final image depth range has not been reached in step 314. Otherwise, the computing system displays the final image of the analyzed subsurface, showing one or more geological features.

[0038] The imaging step 312 involves solving a least square (LSQ) problem that simultaneously models the free-surface and internal multiples as follow: u(f = F FD (f, x, y, z)r(z, y, z), (1) where: -r is the subsurface reflectivity to be found by LSQ inversion (e.g. conjugate gradients); - D is the down-going wavefield extrapolated into the subsurface; - F is the forward extrapolation and accumulation operator of reflecting

wavefield;

- F is the extraction of data at receiver positions (it is assumed that there are n receivers). This operator may relate to spatial sinc functions to extract traces at the recording positions from the propagated data. The output data may be at the surface, the seabed, or at other arbitrary positions; and - u is the data recorded at receivers 135.

The downgoing wavefield (D = DE. D1) is a summation of: -DF is the forward propagated data; for free-surface multiple modelling, it is possible to use the formalism in [6], and - is the reverse-propagated and downward reflected data; for internal multiple modelling, it is possible to use the formalism in [15].

[0039] FIGs. 6A to 6C schematically illustrate the propagation mechanism for free-surface and internal multiple modelling. More specifically, FIG. 6A shows the down-going wavefield DF that models free-surface multiples. The down-going wavefield OF is generated by forward propagating the recorded data do into the subsurface. In this case, the recorded data do relates to data recorded by a plurality of receivers mounted on towed streamers, as discussed above with regard to FIG. 1. In another embodiment, the data may relate to ocean bottom node data in the shot or receiver domain. In yet another embodiment, the data may relate to land data in the shot or receiver domain. In the case the receivers are positioned significantly below the free-surface, they may be forward propagated to complete their path to the free-surface, and then be reflected at the free-surface and further propagated into the subsurface. This step simulates the wavefield DF that would have resulted from the wavefield being reflected at the free-surface 104 and continued into the subsurface.

[0040] The downgoing wavefield Di for internal multiple modelling is illustrated in FIG. 6B and is being generated by the following steps: * Reverse propagating 602 the recorded data do into the upper reflectivity (in this case including horizon R1); * Reflecting downwards 604 in the upper-reflectivity; and * Forward propagating 606 into the lower reflectivity.

[0041] It should be noted that the primary reflection P1 from R1 in FIG. 6B is reverse propagated 602 to the reflection point 604 at R1 and then reflected downwards 606, effectively in this case changing the primary into a pseudo-direct arrival. This pseudo-direct arrival would go on to image the horizon R2 along with deeper reflectors. For this reason, it may be advantageous to remove (e.g., mute or subtract) the primary arrival P1 relating to R1 from the input data (and other events in the upper reflectivity for a more complex case) for the interbed-related forward wavefield propagation, as illustrated in FIG. 6C.

[0042] The method illustrated in FIG. 3 for generating an image of geological features in a subsurface may be modified as now discussed. The method receives the input data, which includes first free-surface multiples, second free-surface multiples, and internal-multiples recorded by receivers following excitation of a seismic source. Next, the method generates a first image of a portion of the subsurface using the input data. Based on the first image, the method models first free-surface multiple data. By subtracting the modelled first free-surface multiple data from the input data, the method generates partial free-surface demultiple data. Next, a second image is generated using the partial free-surface demultiple data. These steps may be repeated for generating multiple images. At the end, these multiple images may be combined to generate a final image of the subsurface.

[0043] The modelled first free-surface multiple data may be generated by using one of SRME, MWD, deconvolution, SRMM, targeted wave-equation multiple modelling, or WEDecon, as disclosed in [1-6]. In one application, the first free-surface multiple data is dependent on the first image of the subsurface. The dependency may relate to defining a multiple generator. In another application, the first free-surface multiples have a shorter multiple period than the second free-surface multiples. The partial free-surface demultiple data has a reduced first free-surface multiple content.

[0044] The step of generating the first image of the subsurface may be based on primary imaging, or free-surface multiple imaging, or jointly on primary imaging and free-surface multiple imaging.

[0045] In one embodiment, the step of generating the second image of the subsurface is based jointly on free-surface multiple imaging and internal multiple imaging. This step is dependent on the first image of the subsurface for a downward reflection within the first image of the subsurface.

[0046] In another embodiment, the step of generating the second image of the subsurface is based jointly on primary imaging, free-surface multiple imaging, and internal multiple imaging and the internal multiple imaging is dependent on the first image of the subsurface.

[0047] In one embodiment, the step of generating the first image of the subsurface includes the following sub-steps: 1. Generating source-side forward propagated data by forward propagating the input data.

2. Generating reverse propagated receiver-side data by reverse propagating the input data.

3. Generating the first image of the subsurface using the source-side forward propagated data with the reverse propagated receiver-side data. This step may involve the use of an imaging condition (e.g., cross-correlation, deconvolution, etc.) or an optimized migration.

[0048] In one embodiment, the step of generating the second image of the subsurface involves a joint free-surface and interbed multiple modelling. The joint free-surface and interbed multiple imaging (second image of the subsurface) includes generating forward propagated joint source-side data where the generation of forward propagated joint source-side data involves: * Generating forward propagated free-surface multiple source-side data by forward propagating the input data; * Generating the forward propagated internal multiple source-side data by: 1. Receiving second input data, 2. Generating reverse propagated second input data by reverse propagating the second input data, 3. Generating downward reflecting data by reflecting the reverse propagated second input data at a depth within the first image of the subsurface, and 4. Generating forward propagated internal multiple source-side data by forward propagating the downward reflecting data.

* Generating the joint forward propagated source-side data by combining (e.g., summing) the forward propagated free-surface multiple source-side data and the forward propagated internal multiple source-side data.

* The second input data may be generated by muting the input data.

[0049] The joint free-surface and interbed multiple imaging (second image of the subsurface) involves reverse propagating the partial free-surface demultiple data.

[0050] In another embodiment, the step of generating the second image of the subsurface involves application of an imaging condition (e.g., cross-correlation, deconvolution, etc.) using the joint forward propagated source-side data with reverse propagated partial free-surface demultiple data.

[0051] In yet another embodiment, the step of generating the second image of the subsurface involves optimizing a cost function based on the joint forward propagated source-side data with the reverse propagated partial free-surface demultiple data.

[0052] The second image of the subsurface may optionally be used to model second multiples. The second multiples are free-surface multiples, and/or internal multiples. In one application, the second multiples are subtracted (straight or adaptively) from the input data to generate second partial-demultiple data. In another application, the second multiples are subtracted (straight or adaptively) from the first partial demultiple data to generate second partial-demultiple data.

[0053] In one embodiment, the modelled first free-surface multiple data and the second multiples are jointly subtracted from the input data. The above process may be repeated to produce a third image of the subsurface using the second partial-demultiple data. In one application, demultiple relating to the above approach is combined with another free-surface or internal demultiple approach, for example as discussed in the background. The second partial-demultiple data is used to produce a final image of the subsurface.

[0054] Using one or more of the features discussed above with regard to the method of FIG. 3, has been applied to actual seismic data as now discussed. Towed streamer data was used first and this data includes free-surface multiples, and internal multiples with the water bottom acting as the upper generator. FIG. 7A shows a migrated image of the input data. FIG. 7A shows free-surface multiple imaging of the data for a depth range from the water bottom 702 to a depth 704, which is 700 m in this case, using data before multiple attenuation as input. Water bottom sparseness weights were used to ensure sharp water bottom imaging [16]. The horizontal gas free-surface multiple generator 706 is also imaged at 400 m depth. Interbed multiples generated between the water bottom and the gas are also visible on the display, highlighted by the white arrows. Note that for the imaging of FIG. 7A, the source-side input used data before demultiple and the receiver-side input also used data before demultiple.

[0055] FIG. 7B shows free-surface multiple imaging results for the water bottom reflection only, derived using image domain sparseness weights [16]. The same source-side and receiver-side data has been used as for FIG. 7A. FIG. 7C shows free-surface multiple imaging results using data before demultiple for the depth range 328m to 700m. Strong water bottom cross-talk energy (black arrow) is visible on this display. The free-surface multiple image of the water bottom (FIG. 7B) was used to predict free-surface water bottom multiples using the approach of [4]. The free-surface water bottom multiples in the input data were attenuated using the water bottom multiple prediction, and an adaptive subtraction was used in this case.

[0056] Free-surface multiple imaging of the data after water bottom free-surface multiple attenuation (partial demultiple data) is illustrated in FIG. 7D, and was obtained using an existing method. Compared to FIG. 7C, the water bottom multiple cross-talk has been significantly reduced. The surface-related gas multiples are still present in the data and feature strongly in the image. The internal multiples highlighted in FIG. 7A are still present in the image.

[0057] The data after water bottom free-surface demultiple (partial demultiple data) was taken forward for joint free-surface and internal multiple imaging (referred to as joint multiple imaging). As described previously, the down-going wavefield for joint imaging includes two parts: (1) a free-surface contribution, and (2) an internal contribution. The free-surface contribution (1) uses a forward propagation of input data including all free-surface and internal multiples. The internal contribution (2) includes: * Reverse propagation of input data after water bottom primary muting to the upper reflectivity depth range (mute water bottom in this case); * Multiplication by the downwards-reflecting water bottom reflectivity (FIG. 7B after polarity reversal); and * Forward propagation of downwards reflecting data into the deeper subsurface.

[0058] The up-going wavefield related to the input data after water bottom free-surface multiple attenuation (the partial demultiple data).

[0059] The result of the joint free-surface and internal (based on an upper reflectivity including the water bottom) imaging is given in FIG. 7E (which corresponds to the second image, step 312 in FIG. 3). Note that the image in FIG. 7E used data before demultiple for the source-side free surface, data before demultiple with water bottom mute for the source-side internal multiples, and data after water bottom free-surface demultiple for the receiver-side. The image contains fewer interbed multiples when compared to the free-surface multiple imaging result in FIG. 7D. The difference between free-surface multiple imaging (FIG. 7D) and joint free-surface and internal multiple imaging (FIG. 7E) is given in FIG. 7F. This highlights the interbed multiple contamination that was removed by the joint multiple imaging process of FIG. 3.

[0060] The reflectivity images were used for multiple modelling using the method of [4] for free-surface multiple modelling, and the method of [9] for interbed multiple modelling. FIG. 8A shows the data before demultiple, FIG. 8B shows the free-surface water bottom multiple model (FSWB, based on reflectivity image of FIG. 7B) and FIG. 8C shows the data after adaptive subtraction of the corresponding multiple model from FIG. 8B. Comparing with FIG. 8A, it can be seen that the water bottom multiple has been largely removed, but many multiples from other generators are present in the data. FIG. 8D shows the multiple prediction between the free-surface and the deeper reflectivity, including the gas (FIG. 7E), obtained with the method of FIG. 3, and FIG. 8E shows the demultiple result after removing the free-surface water bottom and free-surface gas multiple models (the second image in FIG. 3). The image in FIG. 8E has significantly less residual multiple than the result from just the water-bottom multiple (see FIG. 8C). Finally, interbed multiples between the water bottom and the deeper section were modelled as illustrated in FIG. 8F (based on images 11 and 12 in FIG. 3) and then subtracted from the image of FIG. 8E, resulting in the image of FIG. 8G.

[0061] Various technical terms and imaging procedures have been used for describing the method of FIG. 3, and also the images illustrated in FIGs. 7A to 8G. Some of these terms and procedures are now defined for avoiding any confusion.

[0062] The combined or joint imaging of primaries and multiples used in this disclosure may also be called full wavefield migration. In this case, an image of the subsurface is generated based on both primaries and multiples. This approach may involve combining primary and multiple images after migration, or by jointly imaging primary and multiple arrivals in a single imaging step.

[0063] The image of the subsurface may relate to a single image, as defined later by r, but may also relate to an extended image, for example, with subsurface offsets for one or more of the x-, y-, or z-directions and/or with non-zero lag arrivals (e.g., as with tau-gathers). In these cases, the image of the subsurface may be described as: * Subsurface offsets: r(o,, oy, Oz, x, y z); example for subsurface offsets in x, y, and z. * Tau-gathers: r(r,x, y, z) * Alternatively, it is possible to have different images for different surface offset ranges.

[0064] An image of the subsurface is a representation of the reflectivity in the earth, defined in space (x-z for 2D, x-y-z for 3D). This may also be referred to as a migrated image or the reflectivity. An image of the subsurface may be the result of a single-step or optimized migration. The image may also be an extended image as defined above.

[0065] The image has been discussed to be determined by applying an imaging condition. The imaging condition is a mathematical function applied to the forward extrapolated and backwards extrapolated wavefields to form an image. The most common imaging condition is the cross-correlation imaging condition. Other options include deconvolution imaging conditions, varieties of which have been described in the literature, e.g., smoothing imaging condition, and multi-dimensional deconvolution imaging condition.

[0066] The term migration is also used to mean "imaging." Migration is a term used to generate an image of the subsurface from sensor recordings following the excitation of a source. Wave equation migration is commonly applied in the shot domain or in the receiver domain. In the shot domain, the method uses a group of traces corresponding to different sensor positions for one shot excitation. In the receiver domain, the method uses a group of traces for one receiver position corresponding to different shot excitations.

[0067] There are different types of migration, e.g., one-way migration, two-way reverse time migration (RTM), Kirchhoff migration, common reflection angle migration (CRAM), Wave equation Kirchhoff migration (WEK), etc. One type of migration is optimized migration. The migration may use an extended image. Migration may be of primaries, multiples, or a combination of primaries and multiples.

[0068] The method of FIG. 3 refers to a step 310 of modelling. In this regard, an image of the subsurface may be used to model primaries and/or model multiples, e.g. Born modelling. This may also be known as demigration. Modelling may be performed with Kirchhoff (e.g., diffraction modelling), one-way propagation, two-way propagation, etc. [0069] The term multiples was used through this disclosure. Primary arrivals correspond to a signal leaving the source, interacting once with the subsurface (e.g. reflection, refractions or diffractions), and being recorded by the receiver. In the case there is more than one interaction, this event is called a multiple. A first order free-surface multiples may relate to two interactions in the subsurface and one reflection at the free-surface (water-air interface). With higher order multiples, these interactions may take place more than once. Internal multiples may include two upward interactions and one downward interaction in the subsurface. Raw seismic data will include a mixture of primaries, surface related multiples, and intemal multiples. One step in typical processing of seismic data involves modelling and attenuating multiples in the raw data, so that the primary arrivals may be imaged.

[0070] Multiple modelling/prediction, also known as surface related multiple modelling (SRMM) was introduced by [4]. With multiple modelling, it is possible to estimate multiples in the space-time domain using recorded data and an image of the subsurface, as illustrated by the following equation: u = [DFr]F* (2) [0071] In this case the forward extrapolated down-going wavefield, DF, reflects from the image of the subsurface, r, to generate a reflecting wavefield, Dpr. The reflecting wavefield is subsequently extrapolated forward to the up-going, u, recording positions. These modelled up-going data will contain only multiples.

[0072] The method of FIG. 3 also employs the concept of multiple imaging.

Multiple imaging involves calculating a reflectivity based on: (1) forward propagation of a first dataset, and (2) backward propagation of a second dataset. The first dataset at least includes primaries and multiples (may optionally include direct arrival). The second dataset at least includes multiples (may optionally include primary).

[0073] An optimized migration defines an inverse problem where the method tries to find an image of the subsurface, which when demigrated, respects the input data. In other words, the goal is to find an image of the subsurface such that when it is used to model data in the space-time domain, it should equal the recorded data.

[0074] This may be the result of a least-squares inversion, based on the following definitions: * DF(t,x, y, z) : Forward extrapolated down-going wavefield.

* uF(t,x, y, z) : Up-going wavefield, reverse extrapolated into the subsurface.

* r (x, y, z) : Image of the subsurface.

[0075] In the following equation, the forward extrapolated wavefield, DE., is transformed into a reflecting wavefield through multiplication by the image of the subsurface, r. Ideally, the method tries to find an image of the subsurface which results in a reflecting wavefield that equals the reverse extrapolated up-going wavefield, up, for every depth step. An exact solution in practice is generally not possible, hence the approximation: up ti D Fr. (3) [0076] This may be expressed as a minimization of the difference between both sides of the previous equation, for example, with an L2 norm (other norms may be used, e.g., L1, Cauchy, etc): r = arg minlluR -DFr. (4) [0077] Alternatively, this may be expressed as a maximization of the similarity between both sides of the equation, for example, with a zero-lag cross-correlative norm (other norms may be used, e.g., including appropriate normalisation of the up-and down-going data, etc.): r = argmaxIluR DprlIT.o * (5) [0078] It is common to pre-multiply both sides of the minimization equation by the adjoint of DF (i.e., the adjoint modelling operator), as shown below. This allows many stable solvers to be used, e.g conjugate gradients, steepest descent, etc. DTuR DTDFr. (6) [0079] The left-hand side of equation (6), DTuR, may relate to a standard migration with cross-correlation imaging condition. One alternative is to use the deconvolution imaging condition when applying the adjoint modelling operator. This may not strictly be a least-squares problem that would pass the dot-product test, but in practice it may offer faster convergence. The deconvolution imaging condition would also be used for the adjoint' modelling operation on the righthand side, i.e., 14.1"uR DP"DFr. (7) [0080] An alternative formalism may involve evaluating the receiver side

wavefield, for example:

u [Dpdp. (8) [0081] In this case, the up-going wavefield, u, may be at the recording datum, and the reflecting wavefield (Dpr) for each depth in the image of the subsurface is forward extrapolated to the surface and accumulated.

[0082] An alternative way of writing equation (8) may be: it FDEr, (9) where the linear operator F extrapolates the reflecting wavefield from each depth in the image of the subsurface and accumulates at the receiver positions.

[0083] Constraints may be added to the problem, for example, image domain sparseness weights or total variance regularisation. Data domain weights may also be used, for example to respect a mute function or to respect recorded data positions, e.g., following [6].

[0084] The optimization problem may be solved with conjugate gradients, steepest descent, an inverse Hessian approach, or another solver. The first iteration of most optimized migrations may relate to a standard migration with an imaging condition. Non-linear optimized migration algorithms may also be used. Optimized migrations may be used for primary imaging, multiple imaging, or combined primary and multiple imaging as discussed elsewhere in this document.

[0085] The recorded data may be pre-processed prior to imaging. Pre-processing may involve noise attenuation (e.g., swell noise), guided wave attenuation, deblending, source signature compensation, source and/or receiver deghosting, wavefield separation, demultiple, data interpolation/regularization (for down-going and/or up-going wavefields), redatuming (e.g., to free-surface) and other approaches known in the art. Deghosting and/or wavefield separation may be on the source side or receiver side.

[0086] The input data or recorded data used in step 302 of the method of FIG. 3 may be data recorded by any sensor type, examples include hydrophones, geophones, particle motion sensors, particle velocity sensors, accelerometers, near-field hydrophones, near-field accelerometers or other sensor configured to detect seismic energy. The data may be from towed streamer, ocean-bottom sensor (node or cable) acquisition, land acquisition, transition zone campaign, or borehole (e.g., vertical seismic profiling (VSP), distributed acoustic sensing (DAS)).

[0087] The data may be recorded at a constant datum, or a variable datum (e.g., variable depth streamer, or varying topography in land, e.g., floating datum). In the case of land data, the geophone recordings may be propagated forwards and backwards to form the multiple image, and the multiple image may be used to predict multiples. In general, the terms OBN (ocean bottom node), OBC (ocean bottom cable), OBS (ocean bottom survey/sensor), and PRM (permanent reservoir monitoring) systems may be used interchangeably. The recorded data may be input to multiple imaging before or after wavefield separation. OBN receiver gather data may correspond to hydrophone, geophone (x, or, y, or, z), receiver-upgoing, receiver-downgoing, etc. [0088] The term reflectivity is also known as the image of the subsurface, the migrated image, or the image. This is another term for the migrated image, it describes how the down-going wavefield may change into the up-going wavefield.

[0089] The seismic source used for generating the input data may be any type of seismic source. Examples of the seismic source includes an airgun, sparker, boomer, marine vibrator, land vibrator, dynamite, etc. The source may be a single source (e.g., single airgun) or an array of sources (e.g., airgun array).

[0090] The seismic data may be acquired according to a timelapse (also known as 4D) survey. In areas where it is necessary to monitor subsurface changes as a function of time, the seismic acquisitions over the same geographical area may be repeated. These different acquisitions are known as 'vintages'. It is customary to designate one of the vintages as a reference, or 'baseline', while the other vintages may be called 'monitor' datasets.

[0091] Up-going and down-going data or wavefields were used in FIGs. 6A to 6C.

In general, the terms up-going and down-going may be used loosely to mean datasets that will be backwards propagated or forwards propagated as part of the imaging process, respectively. Up-going and down-going data may relate to primary and ghost arrivals separated by a deghosting process (e.g., PZsum and PZminus or another deghosting process). Alternatively, the same data may be used for both the up-going and down-going data, for example hydrophone, Vz, or deghosted data. The data may or may not have been redatumed to the free-surface or another datum.

[0092] Wavefield extrapolation, also known as wavefield propagation was also used with regard to FIGs. 6A to 6C. Seismic data extrapolation is a method to simulate recordings of a wavefield at a position other to which it was recorded. For example, in horizontal towed streamer acquisition, it is possible to record hydrophone measurements at a depth of 12 m. Then, it is possible to extrapolate these measurements to a new datum, for example, at 100 m depth. The extrapolation may use an estimate of the subsurface properties, e.g., velocity, anisotropy, absorption, etc. The extrapolation may be with one-way or two-way extrapolations and may be forward or reverse in direction.

[0093] Any of the above discussed methods may be combined with the method of FIG. 3 and implemented in a computing device, whose schematic is shown in FIG. 9. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations of the above-described methods. The computing device 900 may include a server 901, which has a central processor (CPU) 902 coupled to a random access memory (RAM) 904 and to a read-only memory (ROM) 906. The ROM 906 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. The CPU 902 communicates with other internal and external components through input/output (I/O) circuitry 908 and bussing 910, to receive the data acquired using a streamer, OBN, DAS, land nodes, etc. and output one or more of the multiples, the primaries and/or the image of the geological formation. The processor 902 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.

[0094] The server 901 may also include one or more data storage devices, including hard drives 912, solid state drives 914, and other hardware capable of reading and/or storing information. In one embodiment, software for carrying out the above-discussed methods may be stored and distributed on external hard drive 916, a USB storage device 918, or other form of media capable of portably storing information. Server 901 may be coupled to a display 920, which may be any type of known display or presentation screen, such as LCD, plasma display, cathode ray tubes (CRT), etc. The images of the geological formations or models similar to the ones in Figures 7A to 8G may be shown on a display 920. A user input interface 922 is provided and may include one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc. [0095] Server 901 may be coupled to other devices, such as sources, receivers, etc. The server may be part of a larger network configuration as in a global area network (GAN) such as the Internet 928, which allows ultimate connection to various landline and/or mobile computing devices.

[0096] The term "about" is used in this application to mean a variation of up to 20% of the parameter characterized by this term.

[0097] It will be understood that, although the terms first, second, etc. may be used herein to describe various elements, these elements should not be limited by these terms. These terms are only used to distinguish one element from another. For example, a first object or step could be termed a second object or step, and, similarly, a second object or step could be termed a first object or step, without departing from the scope of the present disclosure. The first object or step, and the second object or step, are both, objects or steps, respectively, but they are not to be considered the same object or step.

[0098] The terminology used in the description herein is for the purpose of describing particular embodiments and is not intended to be limiting. As used in this description and the appended claims, the singular forms "a," "an" and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will also be understood that the term "and/or" as used herein refers to and encompasses any possible combinations of one or more of the associated listed items. It will be further understood that the terms "includes," "including," "comprises" and/or "comprising," when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof. Further, as used herein, the term "if' may be construed to mean "when" or "upon" or "in response to determining" or "in response to detecting," depending on the context.

[0099] The disclosed embodiments enable jointly processing free-surface and internal multiples for imaging a subsurface. It should be understood that this description is not intended to limit the invention. On the contrary, the exemplary embodiments are intended to cover alternatives, modifications and equivalents, which are included in the spirit and scope of the invention as defined by the appended claims. Further, in the detailed description of the exemplary embodiments, numerous specific details are set forth in order to provide a comprehensive understanding of the claimed invention. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

[0100] Although the features and elements of the present exemplary embodiments are described in the embodiments in particular combinations, each feature or element can be used alone without the other features and elements of the embodiments or in various combinations with or without other features and elements disclosed herein.

[0101] This written description uses examples of the subject matter disclosed to enable any person skilled in the art to practice the same, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the subject matter is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims.

[0102] The entire content of all the publications listed herein is incorporated by reference in this patent application.

[1] Berkhout, A. J. and Verschuur, D. J. [1997] Estimation of multiple scattering by iterative inversion, Part I: theoretical consideration. Geophysics, 62(5), 1586-1595.

[2] Wang, P., Jin, H., Xu, S. and Zhang, Y. [2011]. Model-based water-layer demultiple. 81st Annual International Meeting, SEG, Expanded Abstracts.

[3] Biersteker, J. [2001] MAGIC: Shell's surface multiple attenuation technique. 71st SEG Annual International Meeting, Expanded Abstracts, 1301-1304.

[4] Pica, A., Poulain, G., David, B., Magesan, M., Baldock, S., Weisser, T., Hugonnet, P. and Herrmann, P. [2005] 3D surface-related multiple modeling, principles and results. 75th SEG Annual International Meeting, Expanded Abstracts, 2080-2083.

[5] Wiggins, W. [1988]. Attenuation of complex water-bottom multiples by waveequation-based prediction and subtraction. Geophysics, 53, 1527-1539.

[6] Poole, G. [2019] Shallow water surface related multiple attenuation using multisailline 3D deconvolution imaging. 81st EAGE Conference and Exhibition, Extended Abstracts, Tu R1 5.

[7] Jacubowicz, H. [1998] Wave equation predication and removal of interbed multiples. EAGE conference and proceedings.

[8] Weglein, A.B., Gasparotto, F.A., Carvalho, P.M., and Stolt, R.H.,1997, An inverse scattering series method for attenuating multiples in seismic reflection data, Geophysics, 62, no. 6, 1975-1989.

[9] Pica, A. and Delmas, L., 2008, Wave equation based internal multiple modeling in 3D: 78th Meeting, SEG, Expanded Abstracts, 2476-2480.

[10] Muijs, R., J. O.A. Robertsson, and K. Holliger, 2007, Prestack migration of primary and surface-related multiple reflections: Part I-Imaging. Geophysics.

Muijs, R., J. O.A. Robertsson, and K. Holliger, 2007, Prestack depth migration of primary and surface-related multiple reflections: Part II-Identification and removal of residual multiples. Geophysics.

[11] Lu, S., D. Whitmore, A. Valenciano, and N. Chemingui, 2015, Separated-wavefield imaging using primary and multiple energy: The Leading Edge, 34, no. 7,770-778, //dotorg/10.1190/ t1e34070770.1.

[12] Lu, S., D. Whitmore, A. Valenciano, N. Chemingui, and G. Ronholt, 2016, A practical crosstalk attenuation method for separated wavefield imaging: 86th Annual International Meeting, SEG, Expanded Abstracts, 4235.

[13] Lu, S., Liu, F., Chemingui, N., Valenciano, A. and Long, A. [2018] Least-squared full-wavefield migration. The Leading Edge, SEG.

[14] Berkhout, A. J., and D. J. Verschuur, 2011, Full wavefield migration, utilizing surface and internal multiple scattering: 81st Annual International Meeting, SEG, Expanded Abstracts, 3212-3216.

[15] Poole, G. and Tickle, J. [2023] Wave-equation deconvolution for angle-dependent reflectivitiy and internal multiple prection, EAGE conference proceedings.

[16] Poole, G., Moore, H., Blaszczak, E., Kerrison, H., Keynejad, S., Taboga, A., and Chappell, M. [2021] Sparse wave-equation deconvolution imaging for improved shallow water demultiple. EAGE conference proceedings.

Claims (20)

1. WHAT IS CLAIMED IS: 1. A method for imaging geological features in a subsurface (202), the method comprising: receiving (302) seismic data do associated with the subsurface (202), wherein the seismic data do includes primaries (140), first and second free-surface multiples (260, 262), and internal multiples (264); selecting (304) a first image depth range (410); generating (308) a first image li of the subsurface (202), corresponding to the first image depth range (410), the first image li being generated based at least on the first free-surface multiples (260) associated with the first image depth range (410); and generating (312) a second image /2 of the subsurface based jointly on second free-surface multiples (262) and internal multiples (264), wherein the step of generating a second image of the subsurface involves a downgoing reflection within the first image of the subsurface.
2. 2. The method of Claim 1, wherein the step of generating the second image comprises: modelling the first free-surface multiples based on the first image; generating partial demultiple data di by subtracting the modelled first free surface multiples from the seismic data; and generating the second image of the subsurface using the partial demultiple data.
3. 3. The method of Claim 1, wherein the step of generating the second image is based jointly on the second free-surface multiples and the internal multiples.
4. 4. The method of Claim 1, wherein the first free-surface multiples have a shorter multiple period than the second free-surface multiples.
5. 5. The method of Claim 2, wherein the partial demultiple data di has a reduced first free-surface multiples content than the seismic data do.
6. 6. The method of Claim 1, wherein the step of generating the second image of the subsurface further comprises: modelling joint free-surface and internal multiples.
7. 7. The method of Claim 6, wherein the step of modelling joint free-surface and internal multiples comprises: forward propagating joint source-side data to generate joint forward propagated source-side data; reverse propagating the partial demultiple data to obtain reverse propagated partial demultiple data; and applying an imaging condition to the joint forward propagated source-side data and the reverse propagated partial demultiple data.
8. 8. The method of Claim 1, further comprising: modelling second multiples based on the second image of the subsurface, wherein the second multiples includes at least one of first and second free-surface multiples, and internal multiples.
9. 9. The method of Claim 1, wherein the second image /2 of the subsurface corresponds to a second image depth range, which is different from the first image depth range.
10. 10. The method of Claim 1, further comprising: generating a final image of the subsurface based on the first and second images.
11. 11. A computing system (900) for imaging geological features in a subsurface (202), the computing system (900) comprising: an interface (908) configured to receive (302) seismic data do associated with the subsurface (202), wherein the seismic data do includes primaries (140), first and second free-surface multiples (260, 262), and internal multiples (264); and a processor (902) connected to the interface (908) and configured to, select (304) a first image depth range (410); generate (308) a first image h of the subsurface (202), corresponding to the first image depth range (410), the first image h being generated based at least on the first free-surface multiples (260) associated with the first image depth range (410); and generate (312) a second image /2 of the subsurface based jointly on second free-surface multiples (262) and internal multiples (264), wherein the step of generating a second image of the subsurface involves a downgoing reflection within the first image of the subsurface.
12. 12. The system of Claim 11, wherein the processor is further configured to: model the first free-surface multiples based on the first image; generate partial demultiple data di by subtracting the modelled first free surface multiples from the seismic data; and generate the second image of the subsurface using the partial demultiple data.
13. 13. The system of Claim 11, wherein the step of generating the second image is based jointly on the second free-surface multiple and the internal multiples.
14. 14. The system of Claim 13, wherein the first free-surface multiples have a shorter multiple period than the second free-surface multiples.
15. 15. The system of Claim 12, wherein the partial demultiple data di has a reduced first free-surface multiples content than the seismic data do.
16. 16. The system of Claim 11, wherein the processor is further configured to: model joint free-surface and internal multiples.
17. 17. The system of Claim 16, wherein the processor is further configured to: forward propagate joint source-side data to generate joint forward propagated source-side data; reverse propagate the partial demultiple data to obtain reverse propagated partial demultiple data; and apply an imaging condition to the joint forward propagated source-side data and the reverse propagated partial demultiple data.
18. 18. The system of Claim 11, wherein the second image /2 of the subsurface corresponds to a second image depth range, which is different from the first image depth range.
19. 19. A non-transitory computer readable medium including computer executable instructions, wherein the instructions, when executed by a processor, implement a method for imaging geological features in a subsurface (202), the medium comprising instructions for: receiving (302) seismic data do associated with the subsurface (202), wherein the seismic data do includes primaries (140), first and second free-surface multiples (260, 262), and internal multiples (264); selecting (304) a first image depth range (410); generating (308) a first image Ii of the subsurface (202), corresponding to the first image depth range (410), the first image II being generated based at least on the first free-surface multiples (260) associated with the first image depth range (410); and generating (312) a second image /2 of the subsurface based jointly on second free-surface multiples (262) and internal multiples (264), wherein the step of generating a second image of the subsurface involves a downgoing reflection within the first image li of the subsurface.
20. 20. The medium of Claim 19, wherein the step of generating the second image comprises: modelling the first free-surface multiples based on the first image; generating partial demultiple data di by subtracting the modelled first free surface multiples from the seismic data; and generating the second image of the subsurface using the partial demultiple data.Intellectual Property Office Application GB2500333.6 Search report under Section 17 of the Patents Act 1977 Date search completed: 23 June 2025 Claims searched: 1-20 International classification Subclass and subgroup Valid from GO 1 V1/30 01/01/2006Field of searchWorldwide search of patent documents classified in the following areas of the IPC: GO1V Databases used in the preparation of this search report: SEARCH-NPL; SEARCH-PATENT Documents considered to be relevant Patent literature Category Relevant claims Document of relevance A US 2017248714 Al (RICKETT), Ir al Property Office is ting name of the Patent Office Alpo Non-patent iterature Category Relevant claims Document of relevance A Geophysics, 80,01/09/2015, Singh et al, Marchenko imaging: Imaging with primaries, internal multiples, and free-surface multiples, S165-S174, -Categories Letter or symbolXDescriptionDocument indicating lack of novelty or inventive step.Y Document indicating lack of inventive step, if combined with another document of the same category.Member of the same patent family.A Document indicating technological background.P Document published on or after the priority date but before the fling date of the present application.E Earlier application published on or after the filing date of the present application.