System and method for individual-order multiple-modelling based on multiple imaging

Abstract

A method 300 for generating a final image of a subsurface when seismic data is scarce for a shallow region includes calculating 306 an initial image of the subsurface based on input seismic data D ass

Description

[0013] According to an embodiment, there is a method for generating a final image of a subsurface when seismic data is scarce for a shallow region. The method includes calculating an initial image of the subsurface based on input seismic data D associated with the subsurface, and a velocity field V of the subsurface, generating first-order modelled data P using the initial image, the velocity field V, and a source wavelet, which is injected at a source location, generating second-order modelled data M using the first-order modelled data P, the initial image, and the velocity field V, calculating demultiple data by subtracting the second-order modelled data M from the input seismic data D, and generating a final image of the subsurface based on the demultiple data.

[0014] According to another embodiment, there is a computing system for generating a final image of a subsurface when seismic data is scarce for a shallow region. The computing system includes an interface configured to receive input seismic data D associated with the subsurface and to receive a velocity field V of the subsurface, and a processor connected to the interface. The processor is configured to calculate an initial image of the subsurface based on the input seismic data D and the velocity field V, generate first-order modelled data P using the initial image, the velocity field V, and a source wavelet, which is injected at a source location, generate second-order modelled data M using the first-order modelled data P, the initial image, and the velocity field V, calculate demultiple data by subtracting the second-order modelled data M from the input seismic data D, and generate a final image of the subsurface based on the demultiple data.

[0015] 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 herein. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] For a more complete understanding of the present invention, reference is now made to the following descriptions taken in conjunction with the accompanying drawings, in which:

[0017] FIG. 1 is a schematic diagram of a marine seismic survey acquisition system;

[0018] FIG. 2 illustrates primaries and multiples generated by a source, and a streamer not being positioned to record near-offset waves;

[0019] FIG. 3 is a flow chart of a method for generating a final image of a subsurface when near-offset data for a shallow region is scarce or missing;

[0020] FIG. 4 illustrates a method for modelling a source wavelet for generating a primary;

[0021] FIG. 5A illustrates the input seismic data, FIG. 5B illustrates the reconstruction of primaries and multiples based on a traditional approach, and FIGs. 5C to 5E illustrate the one by one reconstruction of primaries and multiples;

[0022] FIG. 6 is a flow chart of a seismic data processing method; and

[0023] FIG. 7 is a schematic diagram of a computing device that implements one or more of the methods discussed in this disclosure. DETAILED DESCRIPTION OF THE INVENTION

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

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

[0026] According to an embodiment, a method for generating a final image of the subsurface, when near offset data is lacking or input data is coarsely sampled, is now introduced. The method relies on deriving a reflectivity of the subsurface based on the multiple periodicity in the recorded seismic data, and injecting a source wavelet into the reflectivity for first-order modelling, which corresponds to the primaries. The process is repeated with the first-order modelling dataset being injected and generating second-order modelling dataset, and so on. After multiple modelling (i.e., individual-order multiple-modelling), the second-order modelling dataset (which in fact corresponds to the first multiples) and other higher-order modelling datasets, if available, may be subtracted from the input data (recorded seismic data) to generate demultiple data. The demultiple data is then used to generate the final image of the subsurface. While the approach is especially useful when input data is lacking in near offsets or is coarsely sampled, it may also be used when the input sampling contains near offsets (even zero offset) and / or is densely spatially sampled. Details of this approach are now discussed with regard to the figures.

[0027] A flow chart of a method 300 for generating the final image of the subsurface based on demultiple data is illustrated in FIG. 3. Note that this method for generating the final image of a subsurface works when seismic data is scarce for a shallow region. The scarce seismic data may relate to missing near offsets and / or coarsely sampled data. The method 300 starts in step 302 by receiving seismic data D associated with a subsurface. The seismic data D includes at least one of: primaries P, free-surface multiples (i.e., waves that reflect at least once from the water-air interface), internal multiples (i.e., waves that reflect at least three times from internal interfaces), source ghosts (i.e., waves generated by an actual source that first reflect from the water-air interface before being recorded by the receivers), and receiver ghosts (i.e., waves that reflect from an internal interface before reflecting from the free surface and then being recorded by the receivers). Any mixture of these types of seismic waves may be part of the input seismic data D. In one embodiment, the input seismic data includes direct arrivals (i.e., waves that propagate from the source directly to the receiver, without any reflection). In another embodiment, the input seismic data D contains primaries P and free-surface multiples M.

[0028] In this or another embodiment, the input seismic data D includes: up-going waves (for example, combination of hydrophone (P) and accelerometer (Vz) data, P+Vz); down-going waves (for example, combination of P and Vz data, P-Vz); only hydrophone data; only particle velocity (Vz, Vx, Vy) data; particle acceleration data; ocean bottom survey data (i.e., from ocean bottom nodes (OBN), ocean bottom cable (OBC), permanent reservoir monitoring (PRM)); towed streamer data; land survey data; transition zone survey data; distributed acoustic sensing (DAS) data; or vertical seismic profile (VSP) data. The input seismic data D may be acquired with a seismic source that is one of an airgun, marine vibrator, land vibroseis, dynamite, sparker, boomer, or another source used in this field.

[0029] The method 300 continues with a step 304 of receiving a velocity field V representing the surveyed subsurface. The velocity field V may be obtained from a previous seismic survey or may be inferred based on theoretical calculations. In one application, the velocity field V may be obtained by compiling time and depth pairs from well data. Any other known method in the field may be used for obtaining the velocity field V.

[0030] In step 306, the method derives an initial image of the subsurface based on the input seismic data D and the velocity field V (also called velocity model V). Note that this initial image is a first image of the subsurface. In one embodiment, an image of the subsurface is understood as including the various reflections that are present between the plural interfaces between layers and / or formations of the subsurface and having different impedances. A goal of the method is to generate a final image of the subsurface, which is more accurate than this initial image. In one application, the initial image of the subsurface is derived using forward propagation of the input seismic data D. In another application, the initial image of the subsurface is derived using backward propagation of the input seismic data D. In one embodiment, the initial image of the subsurface is derived using primary imaging. In another embodiment, the initial image of the subsurface is derived using multiple imaging. In yet another application, the initial image of the subsurface is derived using free-surface multiple reflections. In another application, the initial image of the subsurface is derived using internal multiple reflections. In yet another application, the image of the subsurface is jointly derived by using primaries and multiples. In one embodiment, one of anisotropy, absorption, or s-wave velocity are additionally used.

[0031] In step 308, the method generates first-order modelled data using the initial image of the subsurface from step 302, the velocity field from step 304, and a source wavelet are now discussed. The first-order modelled data relates to a primary reflection. In one application, the first-order modelled data is generated as illustrated in FIG. 4. The source wavelet 410, which may be a spike, is injected at a source location S, and forward propagated 412 into the water 140 and subsurface 402 to generate first-order forward propagated data FD1. The source wavelet 410 represents a seismic signal emitted by a source. It is characterized by shape, amplitude, phase, and frequency content. The first-order forward propagated data FD1 is reflected at a reflection point 414, using the image (i.e., reflectivity) of the subsurface from step 304, to generate first-order reflecting data RD1. The first-order reflecting data RD1 is forward propagated to generate the first-order modelled data MD1, which corresponds to the primary data P. While FIG. 4 shows the generation of the first-order modelled data MD1 due to a reflection at an interface between the water bottom 140 and a first layer 116 of the subsurface 402, other reflection points between other layers also contribute to this data. Also, FIG. 4 shows a single wave 412 being emitted at the source location S due to the source wavelet 410. However, in reality, plural waves 412 are emitted in all directions. The propagations for the first-order modelled data MD1 involve one-way or two-way wavefield propagation.

[0032] The SRMM method described in [2] involves the injection of the recorded data into a modelling approach that forward propagates the injected data into the subsurface, undergoes a reflection in the subsurface, and forward propagates the reflection to the receivers. This approach changes each arrival in the injected data into its following order of arrival (primaries become first order multiples, first order multiples become second order multiples, etc). The approach conveniently models all orders of multiples, but does rely on an adequately sampled input, in terms of offset aperture as well as trace spacing. In the proposed approach, a source wavelet is injected into the modelling approach which results in the simulation of first-order modelled data (primaries). The primaries are reinjected into the modelling approach for a second time, resulting in second-order modelled data (first order multiples). The approach requires more modelling effort than [2], but avoids the sampling limitations previously discussed. Note that the first-order modelled data is left on the propagation grid, and not extracted at the receiver positions. By doing so, the modelling is not limited by the receiver sampling.

[0033] The first-order modelled data MD1 may be generated using one-way modelling, two-way modelling, Kirchhoff modelling (demigration), or another modelling. In one embodiment, the first-order modelled data MD1 is generated with acoustic or elastic modelling. The first-order modelled data is based on wavepropagation in the shot-domain or the receiver domain. In one embodiment, the first- order modelled data are free-surface multiples, or internal multiples, or a combination of free-surface and internal multiples.

[0034] Next, the method 300 in FIG. 3 generates in step 310 the second-order modelled data MD2 (i.e., the first-order multiples Mi in equation (2)) using the first-order modelled data MD1, the initial image, and the velocity field V. The second-order modelled data is generated, similar to the step 308 illustrated in FIG. 4, by forward propagating the first-order modelled data MD1 into the subsurface to generate second-order forward propagated data, reflecting the second-order forward propagated data using the image (reflectivity) of the subsurface from step 304, to generate second-order reflecting data, and then forward propagating the second-order reflecting data to generate the second-order modelled data MD2. Optionally, it is possible to extract the second-order modelled data MD2 at the sampling of the input seismic data D. Propagations for the second-order modelled data involve oneway or two-way wavefield forward propagation. In principle, higher order modelled data may be calculated by repeating these steps, if desired.

[0035] The second-order modelled data may be generated using one-way modelling, two-way modelling, Kirchhoff modelling (demigration), or another modelling. The second-order modelled data is generated with acoustic or elastic modelling. In one embodiment, the second-order modelled data is based on wave-propagation in the shot-domain or the receiver domain. In one embodiment, the second-order modelled data are free-surface multiples, or internal multiples, or a combination of free-surface and internal multiples.

[0036] In step 312, the demultiple data D-M is calculated by subtracting the calculated multiples M (i.e., the second-order modelled data MD2) from the input seismic data D. The subtraction may be a straight subtraction or an adaptive subtraction. In one embodiment, the adaptive subtraction is based on deriving a least-squares filter in the space-time domain or another domain.

[0037] In one embodiment, the subtraction in step 312 additionally involves a second multiple modelling that may come from a second multiple-modelling approach. The second multiple-modelling approach may be a Radon multiple prediction, a Wave-equation propagation as in [2] or [3], Internal multiple prediction (e.g. Weglein et al., 1997), SRME as in [5], Model based multiple prediction as in [6, 7], etc.

[0038] In step 314, a final image of the subsurface is generated based on the demultiple data D-M calculated in step 312. Any known imaging algorithm or procedure may be used for this step. The final image of the subsurface may be indicative of an oil and gas reservoir 144, a fault line, or any other subsurface feature that is characterized by a change in impedance relative to an adjacent formation or volume. Based on this final image, the geologist determines to drill a well, or to place a windmill, or to bury CO2 in a subsurface porous reservoir, or to exploit a given ore or geothermal water (or other resources).

[0039] In the previous embodiment, free-surface multiples are predicted using a first lower reflection in the initial image of the subsurface, an upper reflection at the free-surface, and a second lower reflection in the initial image of the subsurface. In the case of internal multiples, the upper reflection may be within the subsurface (below the free-surface). In one embodiment, internal multiples may be predicted using the initial image of the subsurface as follows: - Forward propagating a source wavelet into the subsurface to generate direct arrival data; - Generating first reflection (first lower reflection) data using the direct arrival data and the initial image of the subsurface; - Forward propagating the first reflection data to generate forward propagated first reflection data; - Generating second reflection data (upper reflection) using the forward propagated first reflection data and the initial image of the subsurface; - Forward propagating the second reflection data to generate forward propagated second reflection data; - Generating third reflection data (second lower reflection) using the forward propagated second reflection data and the initial image of the subsurface; - Forward propagating the third reflection data to generate forward propagated third reflection data; - Extracting multiple prediction data using the forward propagated third reflection data; and - Subtracting multiple prediction data from the input data.

[0040] The reflectivity for the second reflection data may use the negative of the initial image of the subsurface as this reflection is from below the interface rather than from above the interface.

[0041] The method of FIG. 3 is now applied to a synthetic dataset for comparing its output to a traditional method. FIG. 5A shows a synthetic dataset 500 generated based on a constant velocity medium at 1490 ms'1 followed by a reflection at 300 m depth. The reflection coefficient is 0.4. The figure shows the primary reflection P at 400 ms (near-offset timing), followed by a first order multiple M1 at 800 ms, a second order multiple M2 at 1200 ms, and a third order multiple M3 at 1600 ms. The multiples get successively weaker and include polarity changes (shown by the white and black lines) relating to an upper reflection at -1 at the free-surface. Note that the X axis indicates the offset (in meter (m) or km) between the source or the source injection point and the point where the reflected wave is recorded, and the Y axis indicates the travel time (in seconds) for the wave to arrive at the receiver from the source.

[0042] According to a conventional approach, the method receives the input seismic data 500 and derives a reflectivity r of the subsurface based on the multiple periodicity in the data (multiple imaging). Although not shown here, the reflectivity r is in the (x, y, z) domain and represents the subsurface (see, for example, [3] for more information). The reflectivity includes the interface at 300 m depth. With conventional modeling, the recorded data 500 of FIG. 5A is injected into the reflectivity r and the obtained modelled data 502 is shown in FIG. 5B. This modelling step involves forward extrapolating the data 500 into the subsurface, reflecting the forward propagated data at each depth step (e.g., each 2 m for up to 100 Hz frequency spectrum, or 1 m for up to 200 Hz frequency spectrum), and forward extrapolating the reflecting wavefield to the receivers. Due to a restricted recording aperture of the input data 500 (as described previously with respect to FIG. 2), the multiple model at short offsets is reduced in amplitude (see regions 510 and 520 in FIG. 5B). While it is possible to extrapolate near offsets to reduce this effect, extrapolated data is often low frequency, and the resulting conventional modelling is still not accurate. Additionally, data may be interpolated between shots to avoid modelling artifacts, however, this modelling may also be inaccurate.

[0043] However, different from the traditional approach, the approach illustrated in FIG. 3 injects the source wavelet 410 (e.g., a spike) into the reflectivity r for the purpose of generating the first-order modelling dataset MD1 (which corresponds to the primary P), which is shown in FIG. 5C. In this case, the forward propagation, reflection, and forward propagation results in modelling the primary arrivals for the reflection at 300 m depth. Next, method 300 injects the first-order modelling dataset MD1 into the modelling approach, which results in the second-order modelling dataset MD2 (which corresponds to the first order multiple Mi) estimate, which is illustrated in FIG. 5D.

[0044] As the input to this modelling approach was dependent only on the source wavelet 410, and not on the sampling of the receivers 114, the second-order modelling dataset MD2 does not suffer from the near-offset dimming. The second-order modelling dataset MD2 may again be injected into the modelling approach 300 to produce third-order modelling data (corresponding to second order multiples M2), as illustrated in FIG. 5E. This individual-order multiple-modelling may be repeated a couple of times.

[0045] After multiple modelling steps, the method may subtract the second-order modelling result (e.g., the results in FIG. 5D) from the input data D (shown in FIG. 5A) to produce demultiple data (not shown). The subtraction (i.e., step 312 in FIG. 3) may be a numerical ‘straight-subtraction’ or may be an adaptive subtraction. The second-order modelling only creates a targeted multiple model, and many multiples are not present in this multiple prediction. For that reason, the method may jointly subtract one or more other models (e.g., jointly subtract the conventional model, the second-order model, and optionally, the third-order model, or from another modelling approach).

[0046] As generally discussed above, a purpose of seismic exploration is to render the most accurate possible graphic representation of specific portions of the Earth's subsurface geologic structure (also referred to as a GAI). The images produced allow exploration companies to accurately and cost-effectively evaluate a promising target (prospect) for its oil and gas yielding potential (i.e., hydrocarbon deposits 144). FIG. 6 illustrates a flow chart of a general method for seismic exploration (method 600). There are five main steps: a detailed discussion of any one of the process steps would far exceed the scope of this document, but a general overview of the process should aid in understanding where the different aspects of the embodiments can be used. Step 602 of method 600 involves positioning and surveying of the potential site for seismic exploration. In step 604, a determination of what type of seismic energy source should be used, and then causing seismic signals to be transmitted. While method 600 applies equally to both marine and land seismic exploration systems, each will use different types of equipment, especially in generating seismic signals that are used to develop data about the Earth’s subsurface geologic structure. In step 606, data recording occurs. In a first part of this step, receivers 114 receive and most often digitize the data, and in a second part of the step 606, the data is transferred to a recording station. In step 608, data processing occurs. Data processing generally involves enormous amounts of computer processing resources, including the storage of vast amounts of data, multiple processors or computers running in parallel, and include the steps 306 to 314 of FIG. 3. Finally, in step 610, data interpretation occurs and results can be displayed, sometimes in two-dimensional form, more often now in three dimensional form. Four dimensional data presentations (a 3D plot or graph, over time (the fourth dimension) are also possible, when needed to track the effects of other processes, for example.

[0047] The methods discussed herein may be applied not only to the field of subsurface exploration, for example, hydrocarbon exploration and development, geothermal exploration and development, and carbon capture and sequestration, or other natural resource exploration and exploitation. They could also be employed for surveying and monitoring for windfarm applications, both onshore and offshore, and also for medical imaging applications.

[0048] The above-discussed procedures and methods may be implemented in a computing device as illustrated in FIG. 7. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein. The computing device 700 is suitable for performing the activities described in the above embodiments and may include a server 701. Such a server 701 may include a central processor (CPU) 702 coupled to a random access memory (RAM) 704 and to a read-only memory (ROM) 706. ROM 706 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. Processor 702 may communicate with other internal and external components through input / output (I / O) circuitry 708 and bussing 710 to provide control signals and the like. Processor 702 carries out a variety of functions as are known in the art, as dictated by software and / or firmware instructions.

[0049] Server 701 may also include one or more data storage devices, including hard drives 712, CD-ROM drives 714 and other hardware capable of reading and / or storing information, such as DVD, etc. In one embodiment, software for carrying out the above-discussed steps may be stored and distributed on a CD-ROM or DVD 716, a USB storage device 718 or other form of media capable of portably storing information. These storage media may be inserted into, and read by, devices such as CD-ROM drive 714, disk drive 712, etc. Server 701 may be coupled to a display 720, which may be any type of known display or presentation screen, such as LCD, plasma display, cathode ray tube (CRT), etc. A user input interface 722 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.

[0050] Server 701 may be coupled to other devices, such as a seismic source, plural receivers of a streamer, etc. The server may be part of a larger network configuration as in a global area network (GAN) such as the Internet 728, which allows ultimate connection to various landline and / or mobile computing devices.

[0051] As described above, the apparatus 700 may be embodied by a computing device. However, in some embodiments, the apparatus may be embodied as a chip or chip set. In other words, the apparatus may comprise one or more physical packages (e.g., chips) including materials, components and / or wires on a structural assembly (e.g., a baseboard). The structural assembly may provide physical strength, conservation of size, and / or limitation of electrical interaction for component circuitry included thereon. The apparatus may therefore, in some cases, be configured to implement an embodiment of the present invention on a single chip or as a single “system on a chip.” As such, in some cases, a chip or chipset may constitute means for performing one or more operations for providing the functionalities described herein.

[0052] In an example embodiment, the processor 702 may be configured to execute instructions stored in the memory device 704 or otherwise accessible to the processor. Alternatively or additionally, the processor may be configured to execute hard coded functionality. As such, whether configured by hardware or software methods, or by a combination thereof, the processor may represent an entity (e.g., physically embodied in circuitry) capable of performing operations according to an embodiment of the present invention while configured accordingly.

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

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

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

[0056] The disclosed embodiments provide a system and method for generating a final image of the subsurface based on individual-order multiplemodelling of multiple imaging. It should be understood that this description is not intended to limit the invention. On the contrary, the 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 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.

[0057] Although the features and elements of the present 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.

[0058] 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. References The entire content of all the publications listed herein is incorporated by reference in this patent application. [1] Whitmore, N.D., Valenciano, A.A., Sollner, W. and Lu, S.

[2010] Imaging of primaries and multiples using a dual-sensor towed streamer. 80th Annual International Meeting, SEG, Expanded Abstracts, 3187-3192. [2] 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 Annual International Meeting, SEG, Expanded Abstracts, 2080-2084. [3] Poole, G.

[2019] , Shallow water surface related multiple attenuation using multisailline 3D deconvolution imaging. 81st EAGE Conference and Exhibition, Extended Abstracts, Tu R1 5. [4] Weglein, A.B., Araujo Gasparotto, F., Carvalho, P.M. and Stolt, R.H., 1997. An inverse- scattering series method for attenuating multiples in seismic reflection data. Geophysics, 62: 1975-1989. [5] Berkhout, A. J. and Verschuur, D.J.

[1997] , Estimation of multiple scattering by iterative inversion, Part I: theoretical consideration. Geophysics, 62, 1586-1595. [6] Wiggins, W.

[1988] . Attenuation of complex water-bottom multiples by wave-equation-based prediction and subtraction. Geophysics, 53, 1527-1539. [7] Wang, P., Jin, H., Xu, S. and Zhang, Y.

[2011] , Model-based water-layer demultiple. 81 st Annual International Meeting, SEG, Expanded Abstracts.

Claims

1. A method for generating a final image of a subsurface when seismic data is scarce for a shallow region, the method comprising:calculating (306) an initial image of the subsurface (402) based on input seismic data D associated with the subsurface (402), and a velocity field V of the subsurface (402);generating (308) first-order modelled data P using the initial image, the velocity field V, and a source wavelet, which is injected at a source location;generating (310) second-order modelled data M using the first-order modelled data P, the initial image, and the velocity field V;calculating (312) demultiple data by subtracting the second-order modelled data M from the input seismic data D; andgenerating (314) a final image of the subsurface (402) based on the demultiple data.

2. The method of Claim 1, wherein the step of generating second-order modelled data comprises:injecting the first-order modelled data as an areal source.

3. The method of Claim 1, wherein the step of generating first-order modelled data comprises:injecting the source wavelet, at the source location;forward propagating the source wavelet into the subsurface to generate first-order forward propagated data;reflecting the first-order forward propagated data using the initial image of thesubsurface to generate first-order reflecting data; andforward propagating the first-order reflecting data to generate the first-order modelled data.

4. The method of Claim 3, wherein the forward propagating uses a one-way or two-way wavefield propagation.

5. The method of Claim 1, wherein the step of calculating an initial image of the subsurface comprises:using multiple imaging or using free-surface multiple reflections.

6. The method of Claim 1, wherein the step of calculating an initial image of the subsurface comprises:jointly using primary reflections and free-surface multiple reflections.

7. The method of Claim 1, wherein the first-order modelled data is related to a primary reflection.

8. The method of Claim 1, wherein the second-order modelled data is related to multiples.

9. The method of Claim 1, wherein an upper reflection is at the surface or within the subsurface.

10. The method of Claim 1, wherein the initial image includes plural reflectivities corresponding to plural interfaces between layers or formations inside the subsurface.

11. A computing system (700) for generating a final image of a subsurface when seismic data is scarce for a shallow region, the computing system (700) comprising:an interface (708) configured to receive (302) input seismic data D associated with the subsurface (402) and to receive (304) a velocity field V of the subsurface (402); anda processor (702) connected to the interface (708) and configured to, calculate (306) an initial image of the subsurface (402) based on the input seismic data D and the velocity field V;generate (308) first-order modelled data P using the initial image, the velocity field V, and a source wavelet, which is injected at a source location;generate (310) second-order modelled data M using the first-order modelled data P, the initial image, and the velocity field V;calculate (312) demultiple data by subtracting the second-order modelled data M from the input seismic data D; andgenerate (314) a final image of the subsurface (402) based on the demultiple data.

12. The computing system of Claim 11, wherein the processor is furtherconfigured to:inject the first-order modelled data as an areal source.

13. The computing system of Claim 11, wherein the processor is further configured to:inject the source wavelet, at the source location;forward propagate the source wavelet into the subsurface to generate first-order forward propagated data;reflect the first-order forward propagated data using the initial image of the subsurface to generate first-order reflecting data; andforward propagate the first-order reflecting data to generate the first-order modelled data.

14. The computing system of Claim 13, wherein the forward propagating uses a one-way or two-way wavefield propagation.

15. The computing system of Claim 11, wherein the processor is configured to:use multiple imaging or using free-surface multiple reflections for calculating an initial image of the subsurface.

16. The computing system of Claim 11, wherein the first-order modelled data is related to a primary reflection and the second-order modelled data is related to multiples.

17. The computing system of Claim 11, wherein the near-offset is less than 500 m and the shallow region is less than 1000 m.

18. The computing system of Claim 11, wherein the initial image includes plural reflectivities corresponding to plural interfaces between layers or formations inside the subsurface.

19. A non-transitory computer readable medium including computer executable instructions, wherein the instructions, when executed by a processor, implement a method for generating a final image of a subsurface when seismic data is scarce for a shallow region, the medium comprising instructions for:calculating (306) an initial image of the subsurface (402) based on an input seismic data D and a velocity field V;generating (308) first-order modelled data P using the initial image, the velocity field V, and a source wavelet, which is injected at a source location;generating (310) second-order modelled data M using the first-order modelled data P, the initial image, and the velocity field V;calculating (312) demultiple data by subtracting the second-order modelled data M from the input seismic data D; andgenerating (314) a final image of the subsurface (402) based on the demultiple data.

20. The medium of Claim 19, wherein the step of generating first-order modelled data comprises:injecting the source wavelet, at the source location;forward propagating the source wavelet into the subsurface to generate first-order forward propagated data;reflecting the first-order forward propagated data using the initial image of the subsurface to generate first-order reflecting data; andforward propagating the first-order reflecting data to generate the first-order modelled data.

Images