WO2016005815A2 - Reverse time migration and multiple migration based methods - Google Patents

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 D_p and multiple data D_m; calculating (704) 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 (706) 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 (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 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.

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

[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(x_r;t;x_s), 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 D_p(x_r;t;x_s) and multiple reflections (simply called multiple or multiples herein) D_m(x_r;t;x_s). 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 P_s and P_r. The wave-field P_s (Liu et al., 201 1 ), if forward propagated, can be found by solving equations:

and the receiver wave-field P_R is backward propagated and it can be solved, by reducing time:

where x=(x,y,z) is the subsurface imaging location; x_s and x_r denote the source and receiver locations, respectively; v denotes the velocity; and V and V² 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 P_s .^x> t; x_s)P_r .x; t; x_s)dt dx_s . (3) [0027] To simplify the analysis in this embodiment, it is assumed that a shot record contains primaries D_p and its first and second order multiples D_m ¹ and D_m ², respectively, as follows:

D = D_p + D_m = D_p + D_m ¹ + D_m ². (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 (D_PFD_m + D_mlD_m ) dt dx_s + // Ό_ρ Ό^_Β dtdx_s. (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 d_m represents the simulated (estimated) multiples.

[0032] Note that the 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).

[0033] Denoting the simulated multiples as d_m = M_mr, where M_m represents the reverse time multiple demigration operator (RTMDM) described by Eqs. (6a) and (6b), and the image as r = Mj_nM_m (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 D_m ^■ M_mr dx_sdx_rdt = 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 d_m with the measured multiple data D_m: 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 D_m and the simulated multiple data d_m = M_mr 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 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.

[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 M_mr 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)

where || -|| denotes the L₂ norm. Applying a Taylor expansion to Eq. (9), the gradient satisfies

and the Hessian is given by: dWmdrY 2

[0039] It is assumed in the following calculations that the observed (i.e., measured) seismic data D, which includes multiples D_m, initial reflectivity image r₀ and the initial simulated data d_mo = M_mr_Q 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 = M_m ^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: 5d_mi = M_m<5rj. In step 108, a Newton update is performed on the simulated multiple seismic data using expression d_mi = d_mi_ - a 5d_mi and an update on the inverted image using expression r_t = r_i→ - a Sr_{i 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(x_r;t;x_s), which includes primary reflections D_p(x_r;t;x_s) and surface- related multiple reflections D_m(x_r;t;x_s), the RTMM, in the simplified case of the acoustic, isotropic and constant density wave equation, is described by Eqs. (1 a)-(2b) discussed above.

[0049] Eq. (3) can be written as

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 - ^E_m(r , (13) where η ε [0,1] is a weighting scalar, (14)

and dx_sdx_r (15)

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, and d_p = M_pr and d_m = M_mr 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 M_p to image r and simulated multiple data is calculated by applying reverse time multiple demigration operator M_m 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:

dE_p(r) dE_m(r) \

d E(r, ?7) =— η— dr— (1— η)— dr + E_m(r)— E_p(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:

E_m(r + 5r) -

[0053] Applying a Taylor expansion of <5r to the first order, Eq. (17) can be approximated as:

E_m r + fir) - E_m(r)

Jf D_mdt f (M_mr)²dt / M_mr^■ M_mr dt

M_mr dt (18)

[0054] Denoting the RTMM multiple migration operator in Eqs. (1 )-(3) as M_m the gradient of the multiple data objective function can be expressed as

where the relationship J r · M^TD dx = / Mr^■ D dx_sdx_rdt (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 D_p and multiple events D_m, the initial reflectivity image r₀ and the initial simulated data do,_P=M_pro and do,m=M_mr₀ 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: 5d_ltP = M_p<5s; and 5d_ltm = M_m5s_h

[0057] In step 212 the method updates the simulated data using expression d_ltP = di_:P - a5di_:P for the primary data and expression d _m - d _m - a5d_ltm for the multiple data. The same step may also update the inverted image using relation r, = r_h1 - 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 ηι=ηί-ι - (E_m(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 M_pr (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). The objective function in the CLSRTM method is:

f _Dr · D

-— dtdx_sdx_r . (21)

M_pr^■ M_pr dt / D^■ D dt

[0065] Here, M_p denotes the reverse-time demigration (RTDM) operator to predict the seismic data d=M_pr using the stacked reflectivity image r(x) by forward propagating both p_s and p_R wavefields simultaneously:

where g(t) is the source wavelet with a flat spectrum, v denotes the background velocity and V² denotes the Laplacian operator. Note that the receiver wave-field p_R propagates only waves that are generated by the fictitious sources resulting from the multiplication of the time derivative of the source wave-field p_s and the stacked image r.

[0066] For simplicity, it is assumed in this embodiment that 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.

[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

dtdx_sdx_r . (24) j f M_pr^■ M_pr 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=M_pr (Guitton and Verschuur, 2004). Because the modelled data M_pr 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 r₀ (which also can be calculated) and initial simulated data do=M_pr₀. 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 d_hi 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 M_p ^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: <¾ = M_p5s_h 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 ₁ - 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 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.

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

[0080] Alternatively, 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. In this case, 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. 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 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 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 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.

[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: 85^th 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: 83^rd 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: 83^rd 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: 83^rd 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: 83^rd 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: 75^th EAGE Conference & Exhibition, Expanded Abstracts.

- Zhang, Y., L. Duan, and Y. Xie, 2013b, A stable and practical implementation of least-squares reverse time migration: 84^th 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: 84^th 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

WHAT IS CLAIMED IS:

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 D_p and multiple data D_m; calculating (704) 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 (706) 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 (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 d_m 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 d_m 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 D_p with estimated primary seismic data d_p.

6. The method of Claim 5, further comprising:

calculating the estimated primary seismic data d_p based on (i) a reverse time demigration operator M_p 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 M_p and with the reverse time multiple demigration operator M_m; updating the simulated primary seismic data d_p and the simulated seismic multiple data d_m;

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 d_m, updating the simulated primary seismic data d_p, 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 D_m; and

a processor (902) connected to the interface (908) and configured to, calculate (704) 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 (706) 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 (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 d_m 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 d_m 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 D_p with estimated primary seismic data d_p.

15. The device of Claim 14, wherein the processor is further configured to: calculate the estimated primary seismic data d_p based on (i) a reverse time demigration operator M_p 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 M_p and with the reverse time multiple demigration operator M_m; update the simulated primary seismic data d_p and the simulated seismic multiple data d_m;

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 d_m, update the simulated primary seismic data d_p, 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 D_p and multiple data D_m; calculating (704) 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 (706) 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 (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.