WO2016005815A2 - Reverse time migration and multiple migration based methods - Google Patents
Reverse time migration and multiple migration based methods Download PDFInfo
- Publication number
- WO2016005815A2 WO2016005815A2 PCT/IB2015/001367 IB2015001367W WO2016005815A2 WO 2016005815 A2 WO2016005815 A2 WO 2016005815A2 IB 2015001367 W IB2015001367 W IB 2015001367W WO 2016005815 A2 WO2016005815 A2 WO 2016005815A2
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- data
- image
- subsurface
- objective function
- seismic data
- Prior art date
Links
- 238000000034 method Methods 0.000 title claims abstract description 141
- 238000013508 migration Methods 0.000 title claims description 48
- 230000005012 migration Effects 0.000 title claims description 47
- 230000006870 function Effects 0.000 description 46
- 238000003384 imaging method Methods 0.000 description 17
- 238000002310 reflectometry Methods 0.000 description 9
- 230000008569 process Effects 0.000 description 8
- 238000005286 illumination Methods 0.000 description 6
- 238000007493 shaping process Methods 0.000 description 6
- 238000012545 processing Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 230000006872 improvement Effects 0.000 description 4
- 230000003044 adaptive effect Effects 0.000 description 3
- 238000009472 formulation Methods 0.000 description 3
- 238000012067 mathematical method Methods 0.000 description 3
- 239000000203 mixture Substances 0.000 description 3
- 238000007781 pre-processing Methods 0.000 description 3
- 230000002238 attenuated effect Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000002939 conjugate gradient method Methods 0.000 description 2
- 238000005314 correlation function Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 230000000644 propagated effect Effects 0.000 description 2
- 230000001902 propagating effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000001881 scanning electron acoustic microscopy Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000002945 steepest descent method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/56—De-ghosting; Reverberation compensation
Definitions
- Embodiments of the subject matter disclosed herein generally relate to methods and systems for generating an image of a surveyed subsurface and, more particularly, to mechanisms and techniques for processing data related to the subsurface for generating the image with higher quality and enhanced resolution comparative to the available methods.
- Seismic data acquisition and processing may be used to generate a profile (image) of geophysical structures under the ground (subsurface). While this profile does not provide an accurate location for oil and gas reservoirs, it suggests, to those trained in the field, the presence or absence of such reservoirs. Thus, providing a high-resolution image of the subsurface is important, for example, to those who need to determine where oil and gas reservoirs are located.
- RTM True-amplitude reverse time-migration
- RTMM reverse time multiple migration
- LSM least-squares migration
- the method includes receiving measured seismic data D related to the subsurface of the earth, where the measured seismic data D includes primary data D p and multiple data D m ; calculating estimated multiple seismic data d m based on (i) a reverse time multiple demigration operator M m and (ii) an initial image r of the subsurface; generating an objective function E that includes a first term defined by a cross-correlation of the measured multiple data D m with the estimated multiple seismic data d m ; and calculating the optimal image of the subsurface by maximizing or minimizing the objective function E.
- a computing device for generating an optimal image of a subsurface of the earth.
- the computing device includes an interface for receiving measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data D p and multiple data D m ; and a processor connected to the interface.
- the processor is configured to calculate estimated multiple seismic data d m based on (i) a reverse time multiple demigration operator M m and (ii) an initial image r of the subsurface, generate an objective function E that includes a first term defined by a cross-correlation of the measured multiple data D m with the estimated multiple seismic data d m , and calculate the optimal image of the subsurface by maximizing or minimizing the objective function E.
- a computer-readable medium including computer-executable instructions, wherein the instructions, when executed by a processor, implement instructions for generating an image of a subsurface of the earth.
- the instructions implement the method steps discussed above.
- Figure 1 is a flowchart of a correlative least-squares RTMM method
- FIG. 2 is a flowchart of an adaptive joint least-squares RTM method
- Figure 3A illustrates a migration velocity model
- Figure 3B illustrates an image generated by using conventional RTM method
- Figure 3C illustrates an image generated by using conventional RTMM
- Figure 4A illustrates an image obtained using 10 iterations of the
- Figure 5 is a flowchart of a guided correlative LSRTM method
- Figure 6A illustrates an image obtained using the RTM method
- Figure 6B illustrates an image obtained using the GLSRTM method
- Figure 6C illustrates an image obtained using a conventional CLSRTM method
- Figure 7 is a flowchart of a method for calculating an optimal image based on an objective function
- Figure 8 is a flowchart of another method for calculating an optimal image based on an objective function.
- Figure 9 is a schematic diagram of a computing device that implements one or more of the above-methods.
- a correlative least-square reverse time multiple migration (CLSRTMM) method for calculating the image of the subsurface.
- CLSRTMM maximizes a cross-correlation of the simulated and observed multiple data at zero lag.
- an adaptive joint LSRTM (AJLSRTM) method which uses both the primary and multiple reflections as signals and adaptively modify the data matching.
- a guided correlative LSRTM (GLSRTM) method which uses the modelled data from CLSRTM to guide the seismic data to ease the defocusing and incoherency due to inaccurate velocities.
- this method can be used to remove ground-roll residual in land data before maximizing the cross-correlation of the simulated and observed data at zero lag through an iterative inversion scheme. Details about each of these methods are now discussed.
- the CLSRTMM method which is based on CLSRTM, maximizes the cross-correlation of the simulated and observed multiple data at zero lag.
- the concept of zero lag is associated with the cross-correlation function (to be discussed later).
- a cross-correlation function or operator is a measure of similarity of two series as a function of the lag of one relative to the other.
- the cross-correlation may be calculated when the lag between the two series (simulated and observed multiple data) is zero, i.e., zero lag.
- D p x r ;t;x s
- D m x r ;t;x s
- LSM least-squares migration
- the RTDM formulation proposed in Zhang and Duan (2012) is extended herein to predict the seismic multiples from the stacked image, or reflectivity model, r(x), using the following partial differential equations:
- d m represents the simulated (estimated) multiples.
- receiver wave-field P R (x; t; x s ) propagates only waves that are generated by fictitious sources that result from the multiplication of source wave- field P s (x; t;x s ) and reflectivity r(x) on the right hand side of Eq. (6b).
- simulated multiple data d m may be obtained during processing, prior to starting this algorithm.
- Methods for estimating multiple data are known in the art, for example, surface-related multiple elimination (SRME) method and, thus, not repeated herein.
- the RTMM imaging operator also known as the reverse time migration operator.
- step 1 10 the objective function E is checked against a predetermined value and a stopping criterion is estimated.
- the predetermined value may be -1 and the stopping criterion may be whether a difference between the updated E and the predetermined value is smaller than a desired target. If the stopping criterion is met, the process advances to step 1 12, wherein the last calculated image ⁇ is generated and displayed. Otherwise, the process returns to step 104 to calculate a new gradient direction and steps 104 to 108 are repeated for updating again the simulated multiple data and the image.
- the CLSRTMM method is an attractive technique for attenuating crosstalk, improving image resolution and illumination, and suppressing migration artifacts.
- the CLSRTMM method also provides stable solutions. Since RTMDM simulates multiple data with the actual shot record injected as the virtual sources for any acquisition, CLSRTMM method is capable of suppressing the effect of source signature and handling some of the difficult imaging issues caused by acquisition, such as free surface ghosts, acquisition hole and cross-talk caused by simultaneous shooting and multiple migration imaging condition.
- the AJLSRTM method is now discussed. By adaptively maximizing a cross-correlation between the simulated and acquired primary and multiple data, the AJLSRTM method derives a simultaneous steepest descent for both the primary and multiple in seeking the optimal image.
- the novel AJLSRTM method provides high- quality images with balanced amplitudes, improved focusing and enhanced resolution.
- the method is also capable of removing free surface ghosts caused by towed streamer acquisition and filling structures at acquisition gaps.
- the method is more generally formulated than other methods and provides more desirable imaging comparing to LSRTM or least-square reverse time multiple migration (LSRTMM).
- a direct implementation of the conventional amplitude-matching-based Least-Squares RTMM may not be the most efficient way of attenuating the cross-talk.
- the cross-talk generally appears as the peg-leg of the actual event demigrates into the same multiple events on the recorded data.
- the earth is at least a visco-elastic medium with density variations, much more complicated than the models used to propagate acoustic wave-fields in seismic imaging.
- the amplitude matching is never perfect.
- AJLSRTM aids the cross-talk attenuation efficiency using the primary reflections matching while benefits from an improved data matching using a stacked image with improved balanced amplitude.
- a time-domain steepest descent direction is determined and a conjugate gradient simultaneous descent scheme in the image domain is used for generating the image.
- the desired subsurface image g(x) is formed by the cross-correlation between the forward propagation of the (primary or multiple) reflection and the backward propagation of the associated next order of the multiple, while the undesired cross-talk image c(x) is also formed due to the cross-correlation of higher order multiples.
- a cross-correlation based objective function that adaptively matches the closeness between the simulated and observed primaries and between the simulated and observed multiples is given by:
- Eq. (14) describes the cross-correlation objective function part for the primary reflections D p
- Eq. (15) describes the cross-correlation objective function part for the multiple reflections D m
- the simulated primary data is calculated by applying reverse time demigration operator M p to image r and simulated multiple data is calculated by applying reverse time multiple demigration operator M m to image r.
- (13) can be considered as a generalized framework of the LSRTM and LSRTMM methods, wherein LSRTM is performed when ⁇ is 1 , and LSRTMM is performed when ⁇ is 0.
- the scale ⁇ plays a substantial role in efficiently matching different part of the data.
- the first two terms in Eq. 16 are the primary and multiples gradients in the image domain while the last term is used to adjust the weighting scalar.
- the multiple gradient is derived through a small perturbation of the reflectivity ⁇ 5r:
- some of these quantities are determined based on methods known in the art and others are measured.
- step 204 the method calculates the steepest descent direction.
- weighting coefficient /3 ⁇ 4 max(0, J ⁇ 3 ⁇ 4;. ( ⁇ 3 ⁇ 4;. - ) based on the
- the weighting coefficient is part of the conjugate gradients method.
- step 216 the objective function E is checked against a predetermined value and a stopping criterion is estimated.
- the predetermined value may be -1 and the stopping criterion may be whether a difference between the updated E and the predetermined value is smaller than a desired target. If the stopping criterion is met, the process advances to step 218, wherein the last calculated image r is generated and displayed. Otherwise, the process returns to step 204 to calculate a new steepest descent direction and steps 204 to 216 are repeated for updating again the simulated primary data, the simulated multiple data and the inverted image r.
- FIG. 3A-4C An example showing the capabilities of the AJLSRTM method is now presented with regard to Figures 3A-4C.
- the AJLSRTM method is performed on the 2D synthetic data set Sigsbee2A (Paffenholz et al., 2002), which has a complex geology.
- the seismic data are generated using a (fine) stratigraphic velocity model, with 50 m shot spacing, 25 m receiver spacing and about 8,000 maximum offset between the source and the last receiver. Both source and receiver ghosts are recorded at 7.6 m depth.
- the seismic data is migrated using the (smoothed) migration velocity illustrated in Figure 3A.
- the initial RTM image has poor illumination at the subsalt 300, while the initial RTMM image contains quite a large area 302 of cross-talk. After 10 iterations of least-square migration, the
- AJLSRTM method provides the most desired image (see Figure 4C) comparing to the LSRTM method (see Figure 4A) and the LSRTMM method (see Figure 4B). Not only the image illumination at subsalt appears to be more balanced using the AJLSRTM method (box 400 in Figure 4C) compared to using the LSRTM method (see box 402 in Figure 4A), but also the cross-talk in the RTMM method (ellipse 302 in Figure 3C) and the LSRTMM method (ellipse 404 in Figure 4B) are best attenuated using AJLSRTM method (ellipse 406 in Figure 4C).
- Figures 4A-C shows, in the subsalt area, the structures are better imaged and diffractors are better focused using the AJLSRTM method (arrows in Figure 4C) than using LSRTM (arrows in Figure 4A) and LSRTMM (arrows in Figure 4B).
- the AJLSRTM method is an attractive technique for improving image resolution and illumination, and suppressing migration artifacts.
- the AJLSRTM method provides stable solutions.
- the AJLSRTM method is more generalized and better capable of obtaining desired image than using either LSRTM or LSRTMM alone. Since RTMDM can simulate data by injecting data containing ghosts as virtual sources, the AJLSRTM method is capable of handling some of the difficult imaging issues caused by acquisition, such as free surface ghosts for towed streamer and crosstalk in the multiple migration.
- RTM has always been regarded as the state-of-the-art technology for imaging. However, as real acquisition is never perfect, RTM acts as an imperfect inverse operator of the seismic data when seeking the image (Luo and Hale, 2014).
- LSRTM seeks an inverted image which generates simulated data best matching the observed seismic data (Dai et al., 201 1 ).
- LSM provides an attractive imaging improvement (Duan et al., 2014).
- travel-time inconsistency often leads to data-matching-failures and even degraded image (Yousefzadeh and Bancroft, 2012).
- the GLSRTM uses the modelled data from LSRTM to guide the seismic data to ease the defocusing and incoherency due to inaccurate velocities.
- the modelled data can be used to remove ground- roll residual in land data before maximizing the cross-correlation of the simulated and observed data at zero lag through an iterative inversion scheme.
- the proposed GLSRTM method advantageously reveals poorly-imaged structures and improves image coherence and focusing in the presence of velocity imperfectness and ground-roll residual.
- the C LSRTM method proposed by Zhang et al. (2013) seeks an optimal image which minimises an objective function measuring the cross-correlation between simulated primary data M p r (M p is the reverse time demigration operator for primary data) and recorded data D(x r ;t;x s ) (recorded at receiver locations x r due to specified source location x s ).
- M p is the reverse time demigration operator for primary data
- D(x r ;t;x s ) recorded at receiver locations x r due to specified source location x s .
- g(t) is the source wavelet with a flat spectrum
- v denotes the background velocity
- V 2 denotes the Laplacian operator
- the seismic data D(x r ;t;x s ) can be modeled using the acoustic propagation of the wave-field p(x;t;x s ) through an isotropic medium with velocity v and density p, as described by equations: where f(t) is the source signature and V is the gradient operator.
- RTDM system 22a and 22b
- acoustic modelling system 23a and 23b
- the conjugate gradient method may be used, similar to the embodiment illustrated in Figure 2.
- the gradient is part of the conjugate gradient method.
- step 504 the shaping operator S can be calculated using observed data D and the simulated data d h i based on relative time shifts between the two data sets following Hale (2013).
- step 506 the steepest descent direction is calculated first in the time domain and then it is converted to the image domain using the migration operator M p T of the RTM, as derived in Zhang et. al (2013):
- step 508 the weighting coefficient 1 ⁇ 4 is calculated as 1 ⁇ 4
- the objective function E is checked against a predetermined value and a stopping criterion is estimated.
- the predetermined value may be -1 and the stopping criterion may be whether a difference between the updated E and the predetermined value is smaller than a desired target.
- step 518 the process advances to step 518, wherein the last calculated image r is generated and displayed. Otherwise, the process returns to step 504 to calculate a new shaping operator S and steps 504 to 516 are repeated for updating again the simulated data and the inverted image r.
- FIG. 6A An example of real seismic data collected from the central North Sea is now discussed for illustrating the GCLSRTM method.
- the seismic data D are acquired using towed streamer cable and has been bandpass-filtered and regularized. Note that any one of the methods discussed above may be applied not only to marine data, but also to land data.
- the migration/demigration model is TTI (Tilted Transverse Isotropy).
- TTI Transmission Transverse Isotropy
- the image bandwidth is quite limited due to the ghosts effect.
- the image bandwidth is gradually extended using either the CLSRTM method as illustrated in Figure 6C or the GLSRTM method as illustrated in Figure 6B as ghosts modelling are utilized in the proposed RTDM at each iteration.
- a new guided correlative Least-Squares RTM method is introduced by reformulating the conventional CLSRTM to incorporate a shaping operator which matches simulated and shifted recorded data when velocities are inaccurate.
- the GLSRTM method naturally removes ground-roll residual in land data.
- the GLSRTM method can attenuate migration artifacts, enhance image resolution, improve structural continuity and reveal weakly imaged details which are otherwise unachievable using RTM or conventional CLSRTM.
- Figure 7 is a flowchart of a method for generating an optimal image of a subsurface of the earth.
- the method includes a step 702 of receiving measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data D p and multiple data D m ; a step 704 of calculating estimated multiple seismic data d m based on (i) a reverse time multiple demigration operator M m and (ii) an image r of the subsurface; a step 706 of generating an objective function E that includes a first term defined by a cross- correlation of the measured multiple data D m with the estimated multiple seismic data d m ; and a step 708 of calculating the optimal image of the subsurface by maximizing or minimizing the objective function E.
- the cross-correlation may be calculated at zero lag.
- the method may include a step of calculating a gradient direction of the objective function in the image domain, a step of demigrating the gradient direction to calculate the gradient in the shot domain, a step of updating the simulated seismic multiple data d m based on the conjugate gradient in the shot domain; a step of updating an inverted image r based on the conjugate gradient direction in the image domain; and a step of repeating the calculating, demigrating, updating the simulated seismic multiple data d m and updating the inverted image until a criterion is met.
- the optimal image is the last calculated inverted image r when the criterion is met.
- the objective function E includes a second term that is defined as a cross-correlation of the measured primary data D p with estimated primary seismic data d p .
- the method may include a step of calculating the estimated primary seismic data d p based on (i) a reverse time demigration operator M p and (ii) the image r of the subsurface.
- the first term is multiplied by a scalar ⁇ and the second term is multipled by 1 - ⁇ .
- the method may further include one or more of a step of calculating a steepest descent direction of the objective function in the image domain; a step of computing a weighting coefficient based on the steepest descent direction; a step of updating a conjugate gradient direction based on the weighting coefficient and the steepest descent direction; a step of demigrating the updated conjugate gradient migration with the reverse time demigration operator M p and with the reverse time multiple demigration operator M m ; a step of updating the simulated seismic data d p and the simulated seismic multiple data d m ; a step of updating an inverted image r based on the updated conjugate gradient direction; a step of updating a weighting scalar ⁇ ; and a step of repeating the calculating, computing, demigrating, updating the simulated seismic multiple data d m , updating the simulated seismic data d p , and updating the inverted image until a criterion is met.
- the optimal image is the last calculated inverted image r when the
- Figure 8 illustrates another method for generating an optimal image of a subsurface of the earth.
- The includes a step 802 of receiving measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data D p and multiple data D m ; a step 804 of calculating estimated primary seismic data d p based on (i) a reverse time demigration operator M p and (ii) an image r of the subsurface; a step 806 of generating an objective function E that includes a cross-correlation of (i) the measured primary data D p shaped with a shaping filter S, and (ii) the estimated primary seismic data d p ; and a step 808 of calculating the optimal image of the subsurface by maximizing or minimizing the objective function E.
- FIG. 9 Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein.
- the exemplary computing system 900 suitable for performing the activities described in the exemplary embodiments may include a server 901 .
- a server 901 may include 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 processor 902 may communicate with other internal and external components through input/output (I/O) circuitry 908 and bussing 910, to provide control signals and the like.
- the processor 902 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.
- Server 901 may also include one or more data storage devices, including a hard drive 912, CD-ROM drives 914, and other hardware capable of reading and/or storing information such as DVD, etc.
- software for carrying out the above-discussed steps may be stored and distributed on a CD- or DVD-ROM 916, removable memory device 918 or other form of media capable of portably storing information. These storage media may be inserted into, and read by, devices such as the CD-ROM drive 914, the disk drive 912, etc.
- Server 901 may be coupled to a display 920, which may be any type of known display or presentation screen, such as LCD, LED displays, plasma displays, cathode ray tubes (CRT), etc.
- a user input interface 922 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.
- Server 901 may be coupled to other computing devices, such as landline and/or wireless terminals, via a network.
- 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 client devices.
- GAN global area network
- the computing device may be implemented on a vehicle that performs a land seismic survey.
- the disclosed exemplary embodiments provide a system and a method for generating a better resolution image of a surveyed subsurface. By improving the accuracy of the subsurface image, an oil and gas company would have a better sense of where to drill a next well, thus improving the technological process of drilling.
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
A device, medium and method for generating an optimal image of a subsurface of the earth. The method includes receiving (702) measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data Dp and multiple data Dm; calculating (704) estimated multiple seismic data dm based on (i) a reverse time multiple demigration operator Mm and (ii) an initial image r of the subsurface; generating (706) an objective function E that includes a first term defined by a cross-correlation of the measured multiple data Dm with the estimated multiple seismic data dm; and calculating (708) the optimal image of the subsurface by maximizing or minimizing the objective function E.
Description
Reverse Time Migration and Multiple Migration Based Methods
CROSS REFERENCE TO RELATED APPLICATIONS
[0001] The present application claims the benefit of priority under 35 U.S.C. § 1 19(e) to U.S. Provisional Application No. 62/022,219 filed on July 9, 2014, U.S. Provisional Application No. 62/084,213 filed on November 25, 2014, and U.S.
Provisional Application No. 62/1 10,819 filed on February 2, 2015, the entire contents of which are hereby incorporated by reference.
BACKGROUND
TECHNICAL FIELD
[0002] Embodiments of the subject matter disclosed herein generally relate to methods and systems for generating an image of a surveyed subsurface and, more particularly, to mechanisms and techniques for processing data related to the subsurface for generating the image with higher quality and enhanced resolution comparative to the available methods.
DISCUSSION OF THE BACKGROUND
[0003] Seismic data acquisition and processing may be used to generate a profile (image) of geophysical structures under the ground (subsurface). While this profile does not provide an accurate location for oil and gas reservoirs, it suggests, to those trained in the field, the presence or absence of such reservoirs. Thus, providing a high-resolution image of the subsurface is important, for example, to those who need to determine where oil and gas reservoirs are located.
[0004] Various methods exist for processing the seismic data for generating an image of the subsurface. True-amplitude reverse time-migration, RTM, (Zhang et al., 2009 and Xu et al., 201 1 ) is the state-of-the-art technology for imaging and interpreting subtle and complex geologic features using the primaries of the seismic data. An addition to RTM is the reverse time multiple migration (RTMM) method (Liu et al., 201 1 ), which adopts the free-surface multiples as virtual sources. For this method, the seismic data was shown to be able to better illuminate the subsurface areas compared to conventional RTM (Lu, 201 1 ).
[0005] However, it is widely accepted that cross-talk noise is always present in the RTMM image due the cross-correlation of different orders of multiples (Wang et al., 2014). Thus, the goal of many improvements to these methods is to reduce the crosstalks noise.
[0006] Over the years, although the modeling and imaging operators have been improved from Kirchhoff (Shuster, 1993 and Nemeth et al., 1999) and one-way wave equation (Tang, 2008) to recent RTM (Dai et al., 201 1 ), the least-squares migration (LSM), which seeks an inverted image with simulated seismic data to best match the amplitude of the observed data, has been consistently shown to enhance resolution and is capable of removing migration artefacts from imperfect acquisition. These observations make LSM an attractive option for attenuating the unavoidable cross-talk in RTMM.
[0007] Because the real earth model is a lot more complicated than the one used in least-squares RTM (LSRTM), amplitude matching is not an easy task in practice and requires considerable effort in preprocessing the observed and simulated data to make LSRTM work (Dong et al., 2012 and Yao, 2012). These practical issues inspired the more recent correlative LSRTM (CLSRTM) (Zhang et al., 2013), which is based on maximizing the cross-correlation of the simulated and observed data. The CLSRTM method relaxes the amplitude constraints and, thus, it can be applied to real data to seek the optimal image with stable performance.
[0008] However, the CLSRTM only uses primary energy. RTMM uses both primaries and multiples and enhances image quality but suffers from cross-talk. Thus, there is a need to develop alternative methods capable of dealing more efficiently with the cross-talk and to generate even better resolution images.
SUMMARY OF THE INVENTION
[0009] According to an embodiment, there is a method for generating an optimal image of a subsurface of the earth. The method includes receiving measured seismic data D related to the subsurface of the earth, where the measured seismic data D includes primary data Dp and multiple data Dm; calculating estimated multiple seismic data dm based on (i) a reverse time multiple demigration operator Mm and (ii) an initial image r of the subsurface; generating an objective function E that includes a first term defined by a cross-correlation of the measured multiple data Dm with the estimated
multiple seismic data dm; and calculating the optimal image of the subsurface by maximizing or minimizing the objective function E.
[0010] According to another embodiment, there is a computing device for generating an optimal image of a subsurface of the earth. The computing device includes an interface for receiving measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data Dp and multiple data Dm; and a processor connected to the interface. The processor is configured to calculate estimated multiple seismic data dm based on (i) a reverse time multiple demigration operator Mm and (ii) an initial image r of the subsurface, generate an objective function E that includes a first term defined by a cross-correlation of the measured multiple data Dm with the estimated multiple seismic data dm, and calculate the optimal image of the subsurface by maximizing or minimizing the objective function E.
[0011] According to still another embodiment, there is a computer-readable medium including computer-executable instructions, wherein the instructions, when executed by a processor, implement instructions for generating an image of a subsurface of the earth. The instructions implement the method steps discussed above.
BRIEF DESCRIPTION OF THE DRAWINGS
[0012] 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:
[0013] Figure 1 is a flowchart of a correlative least-squares RTMM method;
[0014] Figure 2 is a flowchart of an adaptive joint least-squares RTM method;
[0015] Figure 3A illustrates a migration velocity model, Figure 3B illustrates an image generated by using conventional RTM method and Figure 3C illustrates an image generated by using conventional RTMM;
[0016] Figure 4A illustrates an image obtained using 10 iterations of the
CLSRTM method, Figure 4B illustrates an image obtained using 10 iterations of the CLSRTMM method and Figure 4C illustrates an image obtained using 10 iterations of the AJLSRTM method;
[0017] Figure 5 is a flowchart of a guided correlative LSRTM method;
[0018] Figure 6A illustrates an image obtained using the RTM method, Figure 6B illustrates an image obtained using the GLSRTM method, and Figure 6C illustrates an image obtained using a conventional CLSRTM method;
[0019] Figure 7 is a flowchart of a method for calculating an optimal image based on an objective function;
[0020] Figure 8 is a flowchart of another method for calculating an optimal image based on an objective function; and
[0021] Figure 9 is a schematic diagram of a computing device that implements one or more of the above-methods.
DETAILED DESCRIPTION OF THE INVENTION
[0022] 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 seismic models of the following embodiments are solved, for simplicity, with Newton or conjugate gradients methods for achieving numerical convergence. However, the embodiments to be discussed next are not limited to these mathematical methods, but they may solved with other mathematical methods.
[0023] 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.
[0024] According to an embodiment, it is proposed a correlative least-square reverse time multiple migration (CLSRTMM) method for calculating the image of the subsurface. The CLSRTMM maximizes a cross-correlation of the simulated and observed multiple data at zero lag. According to another embodiment, it is proposed an adaptive joint LSRTM (AJLSRTM) method, which uses both the primary and multiple reflections as signals and adaptively modify the data matching. In still another embodiment, it is proposed a guided correlative LSRTM (GLSRTM) method, which uses the modelled data from CLSRTM to guide the seismic data to ease the
defocusing and incoherency due to inaccurate velocities. In one application, this method can be used to remove ground-roll residual in land data before maximizing the cross-correlation of the simulated and observed data at zero lag through an iterative inversion scheme. Details about each of these methods are now discussed.
CLSRTMM method
[0025] The CLSRTMM method, which is based on CLSRTM, maximizes the cross-correlation of the simulated and observed multiple data at zero lag. The concept of zero lag is associated with the cross-correlation function (to be discussed later). It is known that a cross-correlation function (or operator) is a measure of similarity of two series as a function of the lag of one relative to the other. Thus, the cross-correlation may be calculated when the lag between the two series (simulated and observed multiple data) is zero, i.e., zero lag.
[0026] A shot record D(xr;t;xs), i.e., seismic data recorded with hydrophones, geophones, accelerometers, other sensors or any combination of these sensors, be it on land or marine environment, includes primary reflections Dp(xr;t;xs) and multiple reflections (simply called multiple or multiples herein) Dm(xr;t;xs). For this data, in the simplified case of the acoustic, isotropic and constant density wave equation, the traditional RTMM method can be summarized based on the forward and backward propagation of wave-fields Ps and Pr. The wave-field Ps (Liu et al., 201 1 ), if forward propagated, can be found by solving equations:
and the receiver wave-field PR is backward propagated and it can be solved, by reducing time:
where x=(x,y,z) is the subsurface imaging location; xs and xr denote the source and receiver locations, respectively; v denotes the velocity; and V and V2 denote the gradient and Laplacian operators, respectively. To interpret the subsurface structure, geophysicists utilize the full-stacked image generated using a "cross-correlation" imaging condition (Zhang and Sun, 2009):
*) = JT Ps .x> t; xs)Pr .x; t; xs)dt dxs . (3)
[0027] To simplify the analysis in this embodiment, it is assumed that a shot record contains primaries Dp and its first and second order multiples Dm 1 and Dm 2, respectively, as follows:
D = Dp + Dm = Dp + Dm 1 + Dm 2. (4)
[0028] Using Eq. (4) and denoting the source (receiver) wave-field result from forward (backwards) propagating the event by subscript F (B), the stacked image generated using Eq. (3) can be expressed as follows:
r(*) = f (DPFDm + DmlDm ) dt dxs + // Όρ Ό^Β dtdxs. (5)
[0029] While the first term of Eq. (5) provides the desired stacked image which is better illuminated in some subsurface areas than conventional RTM (Lu, 201 1 ), the second term is widely accepted as the undesired cross-talk noise (Liu et al., 201 1 ).
[0030] It has been shown in many publications (Yao and Jakubowicz, 2012, Zhang et al., 2013 and many references herein) that the migration artefacts resulted from imperfect acquisition using only the primaries in migration, can be attenuated with the least-squares migration (LSM). However, in the RTMM setting, LSM formulation is not an easy extension from LSRTM, as it must be ensured that the reverse time multiple demigration (RTMDM) is adjoint to the proposed RTMM above (Nemeth et al., 1999 and Zhan et al., 2010).
[0031] Thus, according to an embodiment, the RTDM formulation proposed in Zhang and Duan (2012) is extended herein to predict the seismic multiples from the stacked image, or reflectivity model, r(x), using the following partial differential equations:
where dm represents the simulated (estimated) multiples.
[0032] Note that the receiver wave-field PR(x; t; xs) propagates only waves that are generated by fictitious sources that result from the multiplication of source wave- field Ps(x; t;xs) and reflectivity r(x) on the right hand side of Eq. (6b).
[0033] Denoting the simulated multiples as dm = Mmr, where Mm represents the reverse time multiple demigration operator (RTMDM) described by Eqs. (6a) and (6b), and the image as r = MjnMm (where represents the reverse time multiple
migration operator), and following the work of Duan et al. (2014), it can be shown that the proposed RTMDM and RTMM operators are adjoint (also called Hermitian adjoint, Hermitian conjugate or Hermitian transpose) to each other in the cross-correlation setting, i.e.,:
jjj Dm ■ Mmr dxsdxrdt = j M^Mm■ r dx. (7) [0034] Thus, according to an embodiment, it is introduced a cross-correlation based objective function E that cross-correlates the simulated multiple data dm with the measured multiple data Dm: E(r) =— (8)
[0035] Note that while the objective function E's structure is similar to that noted in the CLSRTM method (Zhang et al., 2013), the input to this function is different as noted in Eq. (8).
[0036] The process, as described later, tries to find the optimal image r which maximizes the cross-correlation between the measured seismic multiple data Dm and the simulated multiple data dm = Mmr at zero lag, or equivalently, to minimize the objective function of Eq. (8). In one embodiment, the objective function E is
maximized. Note that simulated multiple data dm may be obtained during processing, prior to starting this algorithm. Methods for estimating multiple data are known in the art, for example, surface-related multiple elimination (SRME) method and, thus, not repeated herein.
[0037] In an advantageous scenario, when the two multiple datasets are identical or with a constant scaling difference, the objective function reaches its minimum, -1. The objective function in Eq. (8) is practical and stable because the value of E is unchanged by rescaling the modeled multiple seismic data Mmr and the exact overall scaling of the seismic data strength can be ignored. The proposed CLSRTMM method can thus be adapted without comprehensive preprocessing operations, which are normally required in the conventional amplitude-matching LSM.
[0038] A numerical solution to Eq. (8) can be found using, for example,
Newton's method. Other methods are known in the computational field for numerically solving Eq. (8). If the approach noted above is adopted, the method calculates the gradient and Hessian of the objective function (8) through a small perturbation 5r of the reflectivity image r, as follows:
E(r + <5r) — E(r)
[0039] It is assumed in the following calculations that the observed (i.e., measured) seismic data D, which includes multiples Dm, initial reflectivity image r0 and the initial simulated data dmo = MmrQ are provided in step 102 of method 100 illustrated in Figure 1 . For various iterations = 1,2, ... , the CLSRTMM calculates in step 104 the gradient direction of the objective function in the image domain for iteration i, i.e., <5rj = Mm T where Mr is, as discussed above,
the RTMM imaging operator, also known as the reverse time migration operator.
[0040] In step 106, the method performs seismic modeling of the gradient direction, i.e., it demigrates the gradient direction to calculate the gradient in the shot domain: 5dmi = Mm<5rj. In step 108, a Newton update is performed on the simulated multiple seismic data using expression dmi = dmi_ - a 5dmi and an update on the inverted image using expression rt = ri→ - a Sri t where the scalar a is found to minimize the objective function (8) using a linear search method.
[0041] In step 1 10, the objective function E is checked against a predetermined value and a stopping criterion is estimated. For example, the predetermined value may be -1 and the stopping criterion may be whether a difference between the updated E and the predetermined value is smaller than a desired target. If the stopping criterion is met, the process advances to step 1 12, wherein the last calculated image η is generated and displayed. Otherwise, the process returns to step 104 to calculate a new gradient direction and steps 104 to 108 are repeated for updating again the simulated multiple data and the image.
[0042] The CLSRTMM method is an attractive technique for attenuating crosstalk, improving image resolution and illumination, and suppressing migration artifacts.
The CLSRTMM method also provides stable solutions. Since RTMDM simulates multiple data with the actual shot record injected as the virtual sources for any acquisition, CLSRTMM method is capable of suppressing the effect of source signature and handling some of the difficult imaging issues caused by acquisition, such as free surface ghosts, acquisition hole and cross-talk caused by simultaneous shooting and multiple migration imaging condition.
AJLSRTM method
[0043] The AJLSRTM method is now discussed. By adaptively maximizing a cross-correlation between the simulated and acquired primary and multiple data, the AJLSRTM method derives a simultaneous steepest descent for both the primary and multiple in seeking the optimal image. The novel AJLSRTM method provides high- quality images with balanced amplitudes, improved focusing and enhanced resolution. The method is also capable of removing free surface ghosts caused by towed streamer acquisition and filling structures at acquisition gaps. The method is more generally formulated than other methods and provides more desirable imaging comparing to LSRTM or least-square reverse time multiple migration (LSRTMM).
[0044] Primary reflections have always been considered as the only useful information for subsurface imaging, while multiple reflections have always been considered as coherent noise which requires to be suppressed in data processing prior to the imaging stage. However, with the RTMM, multiples can be treated as useful signals which can improve illumination of the image, especially the shallow section. On the other hand, the cross-talk due to different orders of surface reflections (i.e., multiples) is unavoidable and difficult to remove in the RTMM image.
[0045] A direct implementation of the conventional amplitude-matching-based Least-Squares RTMM (LSRTMM) may not be the most efficient way of attenuating the cross-talk. The cross-talk generally appears as the peg-leg of the actual event demigrates into the same multiple events on the recorded data. As a result, in most cases, it is only possible to use the amplitude information of the recorded and simulated multiples to attenuate the cross-talk.
[0046] The earth is at least a visco-elastic medium with density variations, much more complicated than the models used to propagate acoustic wave-fields in seismic imaging. As a result, the amplitude matching is never perfect. These practical issues require considerable effort in preprocessing both the observed and simulated data to
correctly use the conventional amplitude-matching-based LSRTM or LSRTMM formulation.
[0047] Thus, according to this embodiment, a new method, which uses both the primary and multiple reflections as signals and adaptively modifies the data matching, is proposed. The new method, AJLSRTM, aids the cross-talk attenuation efficiency using the primary reflections matching while benefits from an improved data matching using a stacked image with improved balanced amplitude. A time-domain steepest descent direction is determined and a conjugate gradient simultaneous descent scheme in the image domain is used for generating the image.
[0048] For a zero-phased and de-signatured shot record with source and receiver ghosts D(xr;t;xs), which includes primary reflections Dp(xr;t;xs) and surface- related multiple reflections Dm(xr;t;xs), the RTMM, in the simplified case of the acoustic, isotropic and constant density wave equation, is described by Eqs. (1 a)-(2b) discussed above.
where the desired subsurface image g(x) is formed by the cross-correlation between the forward propagation of the (primary or multiple) reflection and the backward propagation of the associated next order of the multiple, while the undesired cross-talk image c(x) is also formed due to the cross-correlation of higher order multiples.
[0050] According to an embodiment, a cross-correlation based objective function that adaptively matches the closeness between the simulated and observed primaries and between the simulated and observed multiples is given by:
E(r, ?7) ^E r - (l - ^Em(r , (13) where η ε [0,1] is a weighting scalar, (14)
and dxsdxr (15)
Eq. (14) describes the cross-correlation objective function part for the primary reflections Dp, Eq. (15) describes the cross-correlation objective function part for the multiple reflections Dm, and dp = Mpr and dm = Mmr denote the RTDM and RTMDM
simulated primary and simulated multiple data, respectively (Zhang et. al, 2014). Note that the simulated primary data is calculated by applying reverse time demigration operator Mp to image r and simulated multiple data is calculated by applying reverse time multiple demigration operator Mm to image r. The objective function of Eq. (13) can be considered as a generalized framework of the LSRTM and LSRTMM methods, wherein LSRTM is performed when η is 1 , and LSRTMM is performed when η is 0. The scale η plays a substantial role in efficiently matching different part of the data.
[0051] According to an embodiment, a method for finding the optimal image r, which maximizes the cross-correlation between the observed and the simulated primary and multiples at zero lag, or equivalently, which minimizes the defined objective function (13), is now discussed. A numerical solution for this method can be found using the steepest descent method. Those skilled in the art would know that other methods may be used for obtaining a numerical solution for the image. Applying the total derivative to E, the gradient of the objective function (13) is:
dEp(r) dEm(r) \
d E(r, ?7) =— η— dr— (1— η)— dr + Em(r)— Ep(r)J άη. (16)
[0052] The first two terms in Eq. 16 are the primary and multiples gradients in the image domain while the last term is used to adjust the weighting scalar. The multiple gradient is derived through a small perturbation of the reflectivity <5r:
Em(r + 5r) -
[0053] Applying a Taylor expansion of <5r to the first order, Eq. (17) can be approximated as:
Jf Dmdt f (Mmr)2dt / Mmr■ Mmr dt
Mmr dt (18)
[0054] Denoting the RTMM multiple migration operator in Eqs. (1 )-(3) as Mm the gradient of the multiple data objective function can be expressed as
where the relationship J r · MTD dx = / Mr■ D dxsdxrdt (Zhang et al., 2013b) is applied. Similarly, using Eq. (9) and the gradient in Eq. (17-19), and denoting the RTM migration operator in Zhang et al. (2014) as , the gradient of the primary data
data D, which contains primary events Dp and multiple events Dm, the initial reflectivity image r0 and the initial simulated data do,P=Mpro and do,m=Mmr0 are provided in step 202. In one application, some of these quantities are determined based on methods known in the art and others are measured. Method 200, illustrated in Figure 2, may also include a sub-step of setting an initial conjugate gradient direction variation 5so = 0, and the initial weighting scalar
for the initial iteration. These quantities are then updated as the iteration number / increases. In step 204, the method calculates the steepest descent direction. This may be calculated first in the time domain and then, it may be converted to the image domain using RTM or RTMM operators as derived in Eqs. (16), (19) and (20) to obtain: 5n=dE/dn. In ste 206, the method computes weighting coefficient /¾ =max(0, J <¾;. (<¾;. - ) based on the
steepest descent direction. The weighting coefficient is part of the conjugate gradients method.
[0056] In step 208, the method updates the conjugate gradient direction: <5s, = <5r, + /¼<¾/. and in step 210 the method performs seismic modeling of the conjugate gradient direction, i.e. , demigrates the conjugate gradient direction based on
equations: 5dltP = Mp<5s; and 5dltm = Mm5sh
[0057] In step 212 the method updates the simulated data using expression dltP = di:P - a5di:P for the primary data and expression d m - d m - a5dltm for the multiple data. The same step may also update the inverted image using relation r, = rh1 - ads,, where the scalar a is found to minimize the objective function (13) using the linear searching method.
[0058] In step 214, the method updates the weighting scalar ηι=ηί-ι - (Em(n) - Ep(n)) using Eqs. (14) and (15). In step 216, the objective function E is checked
against a predetermined value and a stopping criterion is estimated. For example, the predetermined value may be -1 and the stopping criterion may be whether a difference between the updated E and the predetermined value is smaller than a desired target. If the stopping criterion is met, the process advances to step 218, wherein the last calculated image r is generated and displayed. Otherwise, the process returns to step 204 to calculate a new steepest descent direction and steps 204 to 216 are repeated for updating again the simulated primary data, the simulated multiple data and the inverted image r.
[0059] An example showing the capabilities of the AJLSRTM method is now presented with regard to Figures 3A-4C. The AJLSRTM method is performed on the 2D synthetic data set Sigsbee2A (Paffenholz et al., 2002), which has a complex geology. The seismic data are generated using a (fine) stratigraphic velocity model, with 50 m shot spacing, 25 m receiver spacing and about 8,000 maximum offset between the source and the last receiver. Both source and receiver ghosts are recorded at 7.6 m depth. The seismic data is migrated using the (smoothed) migration velocity illustrated in Figure 3A. As illustrated in Figure 3B, the initial RTM image has poor illumination at the subsalt 300, while the initial RTMM image contains quite a large area 302 of cross-talk. After 10 iterations of least-square migration, the
AJLSRTM method provides the most desired image (see Figure 4C) comparing to the LSRTM method (see Figure 4A) and the LSRTMM method (see Figure 4B). Not only the image illumination at subsalt appears to be more balanced using the AJLSRTM method (box 400 in Figure 4C) compared to using the LSRTM method (see box 402 in Figure 4A), but also the cross-talk in the RTMM method (ellipse 302 in Figure 3C) and the LSRTMM method (ellipse 404 in Figure 4B) are best attenuated using AJLSRTM method (ellipse 406 in Figure 4C). The comparison of Figures 4A-C shows, in the subsalt area, the structures are better imaged and diffractors are better focused using the AJLSRTM method (arrows in Figure 4C) than using LSRTM (arrows in Figure 4A) and LSRTMM (arrows in Figure 4B).
[0060] Thus, the AJLSRTM method is an attractive technique for improving image resolution and illumination, and suppressing migration artifacts. Using a cross- correlation objective function and adaptively adjusting the matching between the simulated and recorded primary and multiple data, the AJLSRTM method provides stable solutions. The AJLSRTM method is more generalized and better capable of obtaining desired image than using either LSRTM or LSRTMM alone. Since RTMDM
can simulate data by injecting data containing ghosts as virtual sources, the AJLSRTM method is capable of handling some of the difficult imaging issues caused by acquisition, such as free surface ghosts for towed streamer and crosstalk in the multiple migration.
GLSRTM method
[0061] RTM has always been regarded as the state-of-the-art technology for imaging. However, as real acquisition is never perfect, RTM acts as an imperfect inverse operator of the seismic data when seeking the image (Luo and Hale, 2014).
[0062] To overcome this issue and obtain an improved image, LSRTM seeks an inverted image which generates simulated data best matching the observed seismic data (Dai et al., 201 1 ). With an accurate smooth background velocity, LSM provides an attractive imaging improvement (Duan et al., 2014). However, when the velocity is inaccurate, travel-time inconsistency often leads to data-matching-failures and even degraded image (Yousefzadeh and Bancroft, 2012).
[0063] To overcome these problems, the GLSRTM uses the modelled data from LSRTM to guide the seismic data to ease the defocusing and incoherency due to inaccurate velocities. In particular, the modelled data can be used to remove ground- roll residual in land data before maximizing the cross-correlation of the simulated and observed data at zero lag through an iterative inversion scheme. As will be discussed later, the proposed GLSRTM method advantageously reveals poorly-imaged structures and improves image coherence and focusing in the presence of velocity imperfectness and ground-roll residual.
[0064] The C LSRTM method proposed by Zhang et al. (2013) seeks an optimal image which minimises an objective function measuring the cross-correlation between simulated primary data Mpr (Mp is the reverse time demigration operator for primary data) and recorded data D(xr;t;xs) (recorded at receiver locations xr due to specified source location xs). The objective function in the CLSRTM method is:
f Dr · D
-— dtdxsdxr . (21)
Mpr■ Mpr dt / D■ D dt
[0065] Here, Mp denotes the reverse-time demigration (RTDM) operator to predict the seismic data d=Mpr using the stacked reflectivity image r(x) by forward propagating both ps and pR wavefields simultaneously:
where g(t) is the source wavelet with a flat spectrum, v denotes the background velocity and V2 denotes the Laplacian operator. Note that the receiver wave-field pR propagates only waves that are generated by the fictitious sources resulting from the multiplication of the time derivative of the source wave-field ps and the stacked image r.
[0066] For simplicity, it is assumed in this embodiment that the seismic data D(xr;t;xs) can be modeled using the acoustic propagation of the wave-field p(x;t;xs) through an isotropic medium with velocity v and density p, as described by equations:
where f(t) is the source signature and V is the gradient operator.
[0067] In the ideal situation, where the exact velocity and the source information is known (v = v and / = g), the objective function (21 ) is able to ease the amplitude ambiguity and robustly provide image improvement (Zhang et al. , 2013), even though
RTDM system (22a and 22b) and acoustic modelling system (23a and 23b) produce different propagation amplitudes.
[0068] However, when the velocity is inaccurate or imperfect as in most of real applications, a travel-time inconsistency between the simulated and recorded data can be encountered. This inconsistency is not included in objective function (21 ). The lack of accounting for this inconsistency often leads to an undesired LSM result as showed in Yousefzadeh and Bancroft (2012). Thus, according to an embodiment, a novel objective function is introduced which better matches the simulated and recorded data by taking into account the above noted inconsistency by introducing a shaping operator S for this purpose. The novel objective function is as follows:
f Dr · S(r)D
dtdxsdxr . (24) j f Mpr■ Mpr dt / S r)D■ S(r)D dt
[0069] Here, S(r) is the shaping operator and it is designed to provide the minimum of
where σ is a scalar to preserve the amplitude of the
recorded data and C(D) is the matrix representing the convolution with the time series for the data d=Mpr (Guitton and Verschuur, 2004). Because the modelled data Mpr is used to guide the recorded data D, the proposed LSM using objective function (24) is called the guided correlative LSRTM, i.e., GLSRTM.
[0070] To solve the objective function, i.e., to determine the inverted image r, the conjugate gradient method may be used, similar to the embodiment illustrated in Figure 2. Those skilled in the art would know that other mathematical methods may be used to find a numerical solution for the GLSRTM method. Regarding the term "inverted image," it is noted that this term simply means an image determined with an inversion method. Similar to the embodiment discussed with regard to Figure 2, the present method 500, illustrated in Figure 5, includes a step 502 of receiving the observed data D, initial reflectivity image r0 (which also can be calculated) and initial simulated data do=Mpr0. In this step, it is possible to set an initial gradient to zero, i.e., 5so =0. The gradient is part of the conjugate gradient method.
[0071] In step 504, the shaping operator S can be calculated using observed data D and the simulated data dhi based on relative time shifts between the two data sets following Hale (2013). In step 506 the steepest descent direction is calculated first in the time domain and then it is converted to the image domain using the migration operator Mp T of the RTM, as derived in Zhang et. al (2013):
[0072] In step 508, the weighting coefficient ¼ is calculated as ¼
=max(0, / δτ^δτι - Sr^^dx / / drfiridx). In step 510, the conjugate gradient direction is updated based on the previous conjugate gradient and the current inverted image as described by expression: <5s, = 5r, + /¼<¾/- .
[0073] In step 512, the method performs seismic modelling of the conjugate gradient direction, which is demigrating the gradient direction: <¾ = Mp5sh followed by step 514 in which the method updates the simulated data using expression d, = d, - a¾ and also updates the inverted image using expression r, = r 1 - ads,, where the scalar a is found to minimize the objective function (24) using the linear searching method.
[0074] In step 516, the objective function E is checked against a predetermined value and a stopping criterion is estimated. For example, the predetermined value may be -1 and the stopping criterion may be whether a difference between the updated E and the predetermined value is smaller than a desired target. If the stopping criterion is met, the process advances to step 518, wherein the last calculated image r is generated and displayed. Otherwise, the process returns to step 504 to calculate a new shaping operator S and steps 504 to 516 are repeated for updating again the simulated data and the inverted image r.
[0075] An example of real seismic data collected from the central North Sea is now discussed for illustrating the GCLSRTM method. The seismic data D are acquired using towed streamer cable and has been bandpass-filtered and regularized. Note that any one of the methods discussed above may be applied not only to marine data, but also to land data. The migration/demigration model is TTI (Tilted Transverse Isotropy). In the initial RTM image illustrate in Figure 6A, the image bandwidth is quite limited due to the ghosts effect. The image bandwidth is gradually extended using either the CLSRTM method as illustrated in Figure 6C or the GLSRTM method as illustrated in Figure 6B as ghosts modelling are utilized in the proposed RTDM at each iteration. In addition, some of the structure at shallow depth (arrows 600) and above the chalk section (arrows 602) in the initial RTM image (Figure 6A) appears to be weakly imaged and possibly defocused due to velocity inaccuracy. The GLSRTM method is able to improve the focusing and coherence at these places and further provide improved and better-interpretable structures as shown in Figure 6B, while the conventional CLSRTM method, without the guiding step discussed above, fails with a lot of diffraction type artifacts.
[0076] In particular, a large section of faults 604 is weakly imaged in the initial image of Figure 6A. These fault details are retrieved using the GLSRTM (Figure 6B), but not using CLSRTM (Figure 6C). It is also noted further improvements on image illumination and structural continuity compared to the result of the conventional CLSRTM method.
[0077] Thus, according to one or more embodiments, a new guided correlative Least-Squares RTM method is introduced by reformulating the conventional CLSRTM to incorporate a shaping operator which matches simulated and shifted recorded data when velocities are inaccurate. The GLSRTM method naturally removes ground-roll residual in land data. When the velocity is inaccurate, the GLSRTM method can
attenuate migration artifacts, enhance image resolution, improve structural continuity and reveal weakly imaged details which are otherwise unachievable using RTM or conventional CLSRTM.
[0078] The above discussed methods are now briefly summarized with regard to Figures 7 and 8. In one embodiment, Figure 7 is a flowchart of a method for generating an optimal image of a subsurface of the earth. The method includes a step 702 of receiving measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data Dp and multiple data Dm; a step 704 of calculating estimated multiple seismic data dm based on (i) a reverse time multiple demigration operator Mm and (ii) an image r of the subsurface; a step 706 of generating an objective function E that includes a first term defined by a cross- correlation of the measured multiple data Dm with the estimated multiple seismic data dm; and a step 708 of calculating the optimal image of the subsurface by maximizing or minimizing the objective function E.
[0079] The cross-correlation may be calculated at zero lag. The method may include a step of calculating a gradient direction of the objective function in the image domain, a step of demigrating the gradient direction to calculate the gradient in the shot domain, a step of updating the simulated seismic multiple data dm based on the conjugate gradient in the shot domain; a step of updating an inverted image r based on the conjugate gradient direction in the image domain; and a step of repeating the calculating, demigrating, updating the simulated seismic multiple data dm and updating the inverted image until a criterion is met. The optimal image is the last calculated inverted image r when the criterion is met.
[0080] Alternatively, the objective function E includes a second term that is defined as a cross-correlation of the measured primary data Dp with estimated primary seismic data dp. In this case, the method may include a step of calculating the estimated primary seismic data dp based on (i) a reverse time demigration operator Mp and (ii) the image r of the subsurface. In one application, the first term is multiplied by a scalar η and the second term is multipled by 1 - η.
[0081] The method may further include one or more of a step of calculating a steepest descent direction of the objective function in the image domain; a step of computing a weighting coefficient based on the steepest descent direction; a step of updating a conjugate gradient direction based on the weighting coefficient and the steepest descent direction; a step of demigrating the updated conjugate gradient
migration with the reverse time demigration operator Mp and with the reverse time multiple demigration operator Mm; a step of updating the simulated seismic data dp and the simulated seismic multiple data dm; a step of updating an inverted image r based on the updated conjugate gradient direction; a step of updating a weighting scalar η; and a step of repeating the calculating, computing, demigrating, updating the simulated seismic multiple data dm, updating the simulated seismic data dp, and updating the inverted image until a criterion is met. The optimal image is the last calculated inverted image r when the criterion is met.
[0082] Figure 8 illustrates another method for generating an optimal image of a subsurface of the earth. The includes a step 802 of receiving measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data Dp and multiple data Dm; a step 804 of calculating estimated primary seismic data dp based on (i) a reverse time demigration operator Mp and (ii) an image r of the subsurface; a step 806 of generating an objective function E that includes a cross-correlation of (i) the measured primary data Dp shaped with a shaping filter S, and (ii) the estimated primary seismic data dp; and a step 808 of calculating the optimal image of the subsurface by maximizing or minimizing the objective function E.
[0083] The above methods and others may be implemented in a computing system specifically configured to calculate the image of the subsurface. An example of a representative computing system capable of carrying out operations in
accordance with the exemplary embodiments is illustrated in Figure 9. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein.
[0084] The exemplary computing system 900 suitable for performing the activities described in the exemplary embodiments may include a server 901 . Such a server 901 may include 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 processor 902 may communicate with other internal and external components through input/output (I/O) circuitry 908 and bussing 910, to provide control signals and the like. The processor 902 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.
[0085] Server 901 may also include one or more data storage devices, including a hard drive 912, CD-ROM drives 914, 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- or DVD-ROM 916, removable memory device 918 or other form of media capable of portably storing information. These storage media may be inserted into, and read by, devices such as the CD-ROM drive 914, the disk drive 912, etc. Server 901 may be coupled to a display 920, which may be any type of known display or presentation screen, such as LCD, LED displays, plasma displays, cathode ray tubes (CRT), etc. A user input interface 922 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.
[0086] Server 901 may be coupled to other computing devices, such as landline and/or wireless terminals, via a network. 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 client devices. The computing device may be implemented on a vehicle that performs a land seismic survey.
[0087] The disclosed exemplary embodiments provide a system and a method for generating a better resolution image of a surveyed subsurface. By improving the accuracy of the subsurface image, an oil and gas company would have a better sense of where to drill a next well, thus improving the technological process of drilling.
[0088] It should be understood that the above 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.
[0089] 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.
[0090] 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
- Nemeth, T., C. Wu, and G. Schuster, 1999, Least-squares migration of incomplete reflection data: Geophysics, Vol. 64, No. 1 , S208-S221 .
- Dai, W., C. Boonyasiriwat, and G. Schuster, 2010, 3D Multi-source Least-squares Reverse Time Migration: SEG, Expanded Abstracts, 29, 3120-3124.
- Dai, W., X. Wang and G. Schuster, 201 1 , Least-squares migration of multisource data with a deblurring filter: Geophysics, 76, no. 5, 135-146.
- Duan, L., Y. Zhang and L. Peng, 2014, Band-limited impedance perturbation inversion using cross-correlative least-squares RTM: 85th Annual International Meeting, SEG Technical Program Expanded Abstracts, To appear.
- Dong, S., J. Cai, M. Guo, S. Suh, Z. Zhang, B. Wang and Z. Li, 2012, Least-squares reverse time migration: towards true amplitude imaging and improving the resolution: 83rd Annual International Meeting, SEG, Expanded Abstracts, 1425-1429.
- Guitton, A., and D. Verschuur, 2004, Adaptive subtraction of multiples using the L1 - norm: Geophysical Prospecting 52.1 , 27-38.
- Hale, D., 2013, Dynamic warping of seismic images: Geophysics, 78, S105-S1 15.
- Liu, Y., X. Chang, D. Jin, R. He, H. Sun and Y. Zheng, 201 1 , Reverse time migration of multiples for subsalt imaging. Geophysics, 76, No. 5, WB209-WB216.
- Luo, S., and D. Hale, 2014, Least-squares migration in the presence of velocity errors: SEG Technical Program Expanded Abstracts 2014, 3980-3984.
- Lu, S., N. D. Whitmore, A. A. Valenciano and N, Chemingui, 201 1 , Imaging of Primaries and Multiples with 3D SEAM Synthetic: SEG Technical Program Expanded Abstracts, 3217-3221 .
- Schuster, G. T., 1993, Least-squares crosswell migration: 63th Annual International Meeting, SEG Technical Program Expanded Abstracts, 1 10-1 13.
- Tang, Y., 2008, Wave-equation Hessian by phase encoding: 78th Annual
International Meeting, SEG, Expanded Abstracts, 27, 2201-2205.
- Tang, Y., and B. Biodi, 2009, Least-squares migration/inversion of blended data: SEG, Expanded Abstracts, 28, 2859-2863.
- Yousefzadeh, A., and B. John, 2012, Velocity evaluation using least squares prestack migration: SEG Technical Program Expanded Abstracts 2012, 1 -5.
- Wang, Y., X. Chang and H. Hao, Multiple migration without prediction: Geophysics, 79, No. 1 , S1 -S9.
- Xu, S., Y. Zhang and B. Tang, 201 1 , 3D angle gathers from reverse time migration: Geophysics, 76, No. 2, S77-S92.
- Yao, G. and H. Jakubowicz, 2012, Least-squares reverse-time migration: 83rd Annual International Meeting, SEG, Expanded Abstracts, 3406-3410.
- Zhan, G. and G. Schuster, 2010, Skeletonized Least Squares Wave Equation
Migration: SEG, Expanded Abstracts, 29, 3380-3384.
- Zhang, Y. and J. Sun, 2009, Practical issues of reverse time migration: true- amplitude gathers, noise removal and harmonic-source encoding: First Break, 26, 19- 25.
- Zhang, Y. and L. Duan, 2012a, Predicting multiples using a reverse time demigration: 83rd Annual International Meeting, SEG, Expanded Abstracts, 520-524.
- Zhang, Y., G. Roberts, and A. Khalil, 2012b, Compensating for source and receiver ghost effects in reverse time migration: 83rd Annual International Meeting, SEG, Expanded Abstracts, SEG-2012-0626.
- Zhang, Y., L. Duan, and G. Roberts, 2013a, True amplitude reverse time migration- from reflectivity to velocity and impedance perturbations: 75th EAGE Conference & Exhibition, Expanded Abstracts.
- Zhang, Y., L. Duan, and Y. Xie, 2013b, A stable and practical implementation of least-squares reverse time migration: 84th Annual International Meeting, SEG
Technical Program Expanded Abstracts, 3716-3720.
- Zhang, Y., A. Ratcliffe, L. Duan and G. Roberts, 2013c, Velocity and impedance inversion based on true amplitude Reverse Time Migration: 84th Annual International Meeting, SEG Technical Program Expanded Abstracts, 949-953.
- Zhang, Y., L. Duan, and Y. Xie, , 2014, A stable and practical implementation of least-squares reverse time migration: Geophysics, To appear.
Claims
1 . A method for generating an optimal image of a subsurface of the earth, the method comprising:
receiving (702) measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data Dp and multiple data Dm; calculating (704) estimated multiple seismic data dm based on (i) a reverse time multiple demigration operator Mm and (ii) an initial image r of the subsurface;
generating (706) an objective function E that includes a first term defined by a cross-correlation of the measured multiple data Dm with the estimated multiple seismic data dm; and
calculating (708) the optimal image of the subsurface by maximizing or minimizing the objective function E.
2. The method of Claim 1 , wherein the cross-correlation is calculated at zero lag.
3. The method of Claim 1 , further comprising:
calculating a gradient direction of the objective function in an image domain; and
demigrating the gradient direction to calculate a gradient in a shot domain.
4. The method of Claim 3, further comprising:
updating the simulated seismic multiple data dm based on the gradient in the shot domain;
updating an inverted image r based on the gradient direction in the image domain; and
repeating the calculating, demigrating, updating the simulated seismic multiple data dm and updating the inverted image until a criterion is met,
wherein the optimal image is the last calculated inverted image r when the criterion is met.
5. The method of Claim 1 , wherein the objective function E includes a second term that is defined as a cross-correlation of the measured primary data Dp with estimated primary seismic data dp.
6. The method of Claim 5, further comprising:
calculating the estimated primary seismic data dp based on (i) a reverse time demigration operator Mp and (ii) the initial image r of the subsurface.
7. The method of Claim 5, wherein the first term is multiplied by a scalar η and the second term is multipled by 1 - η.
8. The method of Claim 5, further comprising:
calculating a steepest descent direction of the objective function in the image domain;
computing a weighting coefficient based on the steepest descent direction; and updating a conjugate gradient direction based on the weighting coefficient and the steepest descent direction.
9. The method of Claim 8, further comprising:
demigrating the updated conjugate gradient migration with the reverse time demigration operator Mp and with the reverse time multiple demigration operator Mm; updating the simulated primary seismic data dp and the simulated seismic multiple data dm;
updating an inverted image r based on the updated conjugate gradient direction;
updating a weighting scalar η; and
repeating the calculating, computing, demigrating, updating the simulated seismic multiple data dm, updating the simulated primary seismic data dp, and updating the inverted image until a criterion is met,
wherein the optimal image is the last calculated inverted image r when the criterion is met.
10. A computing device (900) for generating an optimal image of a subsurface of the earth, the computing device comprising:
an interface (908) for receiving (702) measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data Dp and multiple data Dm; and
a processor (902) connected to the interface (908) and configured to, calculate (704) estimated multiple seismic data dm based on (i) a reverse time multiple demigration operator Mm and (ii) an initial image r of the subsurface,
generate (706) an objective function E that includes a first term defined by a cross-correlation of the measured multiple data Dm with the estimated multiple seismic data dm, and
calculate (708) the optimal image of the subsurface by maximizing or minimizing the objective function E.
1 1. The device of Claim 10, wherein the cross-correlation is calculated at zero lag.
12. The device of Claim 10, wherein the processor is further configured to: calculate a gradient direction of the objective function in an image domain; and demigrate the gradient direction to calculate a gradient in a shot domain.
13. The device of Claim 12, wherein the processor is further configured to: update the simulated seismic multiple data dm based on the gradient in the shot domain;
update an inverted image r based on the gradient direction in the image domain; and
repeat the calculating, demigrating, updating the simulated seismic multiple data dm and updating the inverted image until a criterion is met,
wherein the optimal image is the last calculated inverted image r when the criterion is met.
14. The device of Claim 10, wherein the objective function E includes a second term that is defined as a cross-correlation of the measured primary data Dp with estimated primary seismic data dp.
15. The device of Claim 14, wherein the processor is further configured to:
calculate the estimated primary seismic data dp based on (i) a reverse time demigration operator Mp and (ii) the initial image r of the subsurface.
16. The device of Claim 14, wherein the first term is multiplied by a scalar η and the second term is multipled by 1 - η.
17. The device of Claim 14, wherein the processor is further configured to: calculate a steepest descent direction of the objective function in the image domain;
compute a weighting coefficient based on the steepest descent direction; and update a conjugate gradient direction based on the weighting coefficient and the steepest descent direction.
18. The device of Claim 17, wherein the processor is further configured to: demigrate the updated conjugate gradient migration with the reverse time demigration operator Mp and with the reverse time multiple demigration operator Mm; update the simulated primary seismic data dp and the simulated seismic multiple data dm;
update an inverted image r based on the updated conjugate gradient direction; update a weighting scalar η; and
repeat the calculate, compute, demigrate, update the simulated seismic multiple data dm, update the simulated primary seismic data dp, and update the inverted image until a criterion is met,
wherein the optimal image is the last calculated inverted image r when the criterion is met.
19. A non-transitory computer readable medium including computer executable instructions, wherein the instructions, when executed by a computer, implement a method for generating an optimal image of a subsurface of the earth, the method comprising:
receiving (702) measured seismic data D related to the subsurface of the earth, wherein the measured seismic data D includes primary data Dp and multiple data Dm; calculating (704) estimated multiple seismic data dm based on (i) a reverse time multiple demigration operator Mm and (ii) an initial image r of the subsurface;
generating (706) an objective function E that includes a first term defined by a cross-correlation of the measured multiple data Dm with the estimated multiple seismic data dm; and
calculating (708) the optimal image of the subsurface by maximizing or minimizing the objective function E.
20. The medium of Claim 19, further comprising:
calculating a gradient direction of the objective function in an image domain; and
demigrating the gradient direction to calculate a gradient in a shot domain.
Applications Claiming Priority (6)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US201462022219P | 2014-07-09 | 2014-07-09 | |
US62/022,219 | 2014-07-09 | ||
US201462084213P | 2014-11-25 | 2014-11-25 | |
US62/084,213 | 2014-11-25 | ||
US201562110819P | 2015-02-02 | 2015-02-02 | |
US62/110,819 | 2015-02-02 |
Publications (2)
Publication Number | Publication Date |
---|---|
WO2016005815A2 true WO2016005815A2 (en) | 2016-01-14 |
WO2016005815A3 WO2016005815A3 (en) | 2016-03-10 |
Family
ID=54329849
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
PCT/IB2015/001367 WO2016005815A2 (en) | 2014-07-09 | 2015-07-03 | Reverse time migration and multiple migration based methods |
Country Status (1)
Country | Link |
---|---|
WO (1) | WO2016005815A2 (en) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107589443A (en) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | Method and system based on elastic wave least square reverse-time migration imaging |
CN110133713A (en) * | 2019-04-24 | 2019-08-16 | 中国石油大学(华东) | A kind of the multiple wave least square reverse-time migration imaging method and system of full propagation path attenuation compensation |
US11733413B2 (en) | 2021-04-30 | 2023-08-22 | Saudi Arabian Oil Company | Method and system for super resolution least-squares reverse time migration |
-
2015
- 2015-07-03 WO PCT/IB2015/001367 patent/WO2016005815A2/en active Application Filing
Non-Patent Citations (25)
Title |
---|
DAI, W.; C. BOONYASIRIWAT; G. SCHUSTER: "3D Multi-source Least-squares Reverse Time Migration", SEG, EXPANDED ABSTRACTS, vol. 29, 2010, pages 3120 - 3124 |
DAI, W.; X. WANG; G. SCHUSTER: "Least-squares migration of multisource data with a deblurring filter", GEOPHYSICS, vol. 76, no. 5, 2011, pages 135 - 146 |
DONG, S.; J. CAI; M. GUO; S. SUH; Z. ZHANG; B. WANG; Z. LI: "Least-squares reverse time migration: towards true amplitude imaging and improving the resolution", 83RD ANNUAL INTERNATIONAL MEETING, SEG, EXPANDED ABSTRACTS, 2012, pages 425 - 1429 |
DUAN, L.; Y. ZHANG; L. PENG: "Band-limited impedance perturbation inversion using cross-correlative least-squares RTM", 85TH ANNUAL INTERNATIONAL MEETING, SEG TECHNICAL PROGRAM EXPANDED ABSTRACTS, TO APPEAR, 2014 |
GUITTON, A.; D. VERSCHUUR: "Adaptive subtraction of multiples using the Ll - norm", GEOPHYSICAL PROSPECTING, vol. 52.1, 2004, pages 27 - 38 |
HALE, D.: "Dynamic warping of seismic images", GEOPHYSICS, vol. 78, 2013, pages S105 - S115 |
LIU, Y.; X. CHANG; D. JIN; R. HE; H. SUN; Y. ZHENG: "Reverse time migration of multiples for subsalt imaging", GEOPHYSICS, vol. 76, no. 5, 2011, pages WB209 - WB216 |
LU, S.; N. D. WHITMORE; A. A. VALENCIANO; N, CHEMINGUI: "Imaging of Primaries and Multiples with 3D SEAM Synthetic", SEG TECHNICAL PROGRAM EXPANDED ABSTRACTS, 2011, pages 3217 - 3221 |
LUO, S.; D. HALE: "Least-squares migration in the presence of velocity errors", SEG TECHNICAL PROGRAM EXPANDED ABSTRACTS, 2014, pages 3980 - 3984 |
NEMETH, T.; C. WU; G. SCHUSTER: "Least-squares migration of incomplete reflection data", GEOPHYSICS, vol. 64, no. 1, 1999, pages S208 - S221 |
SCHUSTER, G. T.: "Least-squares crosswell migration", 63TH ANNUAL INTERNATIONAL MEETING, SEG TECHNICAL PROGRAM EXPANDED ABSTRACTS, 1993, pages 110 - 113 |
TANG, Y.: "Wave-equation Hessian by phase encoding", 78TH ANNUAL INTERNATIONAL MEETING, SEG, EXPANDED ABSTRACTS, vol. 27, 2008, pages 2201 - 2205 |
TANG, Y.; B. BIODI: "Least-squares migration/inversion of blended data", SEG, EXPANDED ABSTRACTS, vol. 28, 2009, pages 2859 - 2863 |
WANG, Y.; X. CHANG; H. HAO: "Multiple migration without prediction", GEOPHYSICS, vol. 79, no. 1, pages S1 - S9 |
XU, S.; Y. ZHANG; B. TANG: "3D angle gathers from reverse time migration", GEOPHYSICS, vol. 76, no. 2, 2011, pages S77 - S92 |
YAO, G.; H. JAKUBOWICZ: "Least-squares reverse-time migration", 83RD ANNUAL INTERNATIONAL MEETING, SEG, EXPANDED ABSTRACTS, 2012, pages 3406 - 3410 |
YOUSEFZADEH, A.; B. JOHN: "Velocity evaluation using least squares prestack migration", SEG TECHNICAL PROGRAM EXPANDED ABSTRACTS, 2012, pages 1 - 5 |
ZHAN, G.; G. SCHUSTER: "Skeletonized Least Squares Wave Equation Migration", SEG, EXPANDED ABSTRACTS, vol. 29, 2010, pages 3380 - 3384 |
ZHANG, Y.; A. RATCLIFFE; L. DUAN; G. ROBERTS: "Velocity and impedance inversion based on true amplitude Reverse Time Migration", 84TH ANNUAL INTERNATIONAL MEETING, SEG TECHNICAL PROGRAM EXPANDED ABSTRACTS, 2013, pages 949 - 953 |
ZHANG, Y.; G. ROBERTS; A. KHALIL: "Compensating for source and receiver ghost effects in reverse time migration", 83RD ANNUAL INTERNATIONAL MEETING, SEG, EXPANDED ABSTRACTS, 2012 |
ZHANG, Y.; J. SUN: "Practical issues of reverse time migration: true-amplitude gathers, noise removal and harmonic-source encoding", FIRST BREAK, vol. 26, 2009, pages 19 - 25 |
ZHANG, Y.; L. DUAN: "Predicting multiples using a reverse time demigration", 83RD ANNUAL INTERNATIONAL MEETING, SEG, EXPANDED ABSTRACTS, 2012, pages 520 - 524 |
ZHANG, Y.; L. DUAN; G. ROBERTS: "True amplitude reverse time migration-from reflectivity to velocity and impedance perturbations", 75TH EAGE CONFERENCE & EXHIBITION, EXPANDED ABSTRACTS, 2013 |
ZHANG, Y.; L. DUAN; Y. XIE: "A stable and practical implementation of least-squares reverse time migration", 84TH ANNUAL INTERNATIONAL MEETING, SEG TECHNICAL PROGRAM EXPANDED ABSTRACTS, 2013, pages 3716 - 3720 |
ZHANG, Y.; L. DUAN; Y. XIE: "A stable and practical implementation of least-squares reverse time migration", GEOPHYSICS, TO APPEAR, 2014 |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107589443A (en) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | Method and system based on elastic wave least square reverse-time migration imaging |
CN110133713A (en) * | 2019-04-24 | 2019-08-16 | 中国石油大学(华东) | A kind of the multiple wave least square reverse-time migration imaging method and system of full propagation path attenuation compensation |
US11733413B2 (en) | 2021-04-30 | 2023-08-22 | Saudi Arabian Oil Company | Method and system for super resolution least-squares reverse time migration |
Also Published As
Publication number | Publication date |
---|---|
WO2016005815A3 (en) | 2016-03-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
EP3168653B1 (en) | Device and method for full waveform inversion | |
Pestana et al. | Time evolution of the wave equation using rapid expansion method | |
RU2582480C2 (en) | Coding of simultaneous sources and separation of sources as practical solution for full wave field inversion | |
US9541661B2 (en) | Device and method for deghosting variable depth streamer data | |
Luo et al. | Least-squares migration in the presence of velocity errors | |
US20170031045A1 (en) | Method and apparatus for modeling and separation of primaries and multiples using multi-order green's function | |
Kim et al. | Frequency-domain reverse-time migration with source estimation | |
Yang et al. | Illumination compensation for image-domain wavefield tomography | |
US20170115418A1 (en) | Seismic data processing using matching filter based cost function optimization | |
Yao et al. | Least-squares reverse-time migration | |
da Costa Filho et al. | Elastic P-and S-wave autofocus imaging with primaries and internal multiples | |
AU2014201896A1 (en) | Device and method for stable least-squares reverse time migration | |
Jia et al. | A practical implementation of subsalt Marchenko imaging with a Gulf of Mexico data set | |
US10976457B2 (en) | Full waveform inversion of seismic data using partial match filtering | |
US8634271B2 (en) | Variable depth streamer SRME | |
Hou et al. | Inversion velocity analysis in the subsurface-offset domain | |
Zeng et al. | A guide to least-squares reverse time migration for subsalt imaging: Challenges and solutions | |
Kim et al. | Estimated source wavelet‐incorporated reverse‐time migration with a virtual source imaging condition | |
EP3014308A1 (en) | Processing survey data containing ghost data | |
Araujo et al. | Symplectic scheme and the Poynting vector in reverse-time migration | |
Liu et al. | Eliminating the redundant source effects from the cross-correlation reverse-time migration using a modified stabilized division | |
Tang et al. | Target-oriented wavefield tomography using synthesized Born data | |
EP2622378A1 (en) | Hybrid method for full waveform inversion using simultaneous and sequential source method | |
WO2016005815A2 (en) | Reverse time migration and multiple migration based methods | |
Guo et al. | Time-domain extended-source full-waveform inversion: Algorithm and practical workflow |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 15781400 Country of ref document: EP Kind code of ref document: A2 |
|
NENP | Non-entry into the national phase |
Ref country code: DE |
|
122 | Ep: pct application non-entry in european phase |
Ref document number: 15781400 Country of ref document: EP Kind code of ref document: A2 |