From surface seismic data to reservoir elastic parameters using a full-wavefield redatuming approach

Traditionally, overburden target-oriented In this paper, we present a full-waveﬁeld approach, called reservoir-oriented joint migration inversion (JMI- res ), to estimate the high-resolution reservoir elastic parameters from surface seismic data. As a ﬁrst step in JMI- res , we reconstruct the fully redatumed data (local impulse responses) at a suitable depth above the reservoir from the surface seismic data, while correctly accounting for the overburden interal multiples and transmission losses. Next, we apply a localized elastic full waveform inversion on the estimated impulse responses to get the elastic parameters. We show that JMI- res thus provides much more reliable local target impulse responses, thus yielding high-resolution elastic parameters, compared to a standard redatuming procedure based on time reversal of data. Moreover, by using this kind of approach we avoid the need to apply a full elastic full waveform inversion-type process for the whole subsurface, as within JMI- res elastic full waveform inversion is only restricted to the reservoir target domain.

one-way up-/downgoing wavefields at the depth above the reservoir via JMI. These estimated wavefields are equivalent to the up-/downgoing decomposed Green's functions with real sources at the surface and the virtual receivers in the earth subsurface. Then, the estimated up-/downgoing wavefields are used to estimate the accurate local impulse responses via least-squares based multidimensional deconvolution that is called proximity transformation (van der Neut & Herrmann 2013;da Costa et al. 2015;Garg & Verschuur 2016;van der Neut et al. 2017). The estimated local impulse responses represent fully redatumed data as both sources and receiver are inside the earth subsurface. They are free of spurious events related to multiple reflected energy from the overburden as JMI properly accounts for all the complex propagation (transmission) and scattering effects (internal multiples). In addition, JMI also updates the velocity model of the overburden and provides an accurate propagation velocity model for the target area . Finally, the internal multiples imprint-free estimated impulse responses are used to estimate the reservoir elastic parameters using the target-oriented full waveform inversion (Gisolf & van der Berg 2012;Haffinger 2013). There is also another class of redatuming method called Marchenko-based redatuming (Wapenaar et al. 2014;Ravasi 2017) that can construct Green's functions with virtual sources inside the subsurface and with receivers at the surface while correctly accounting for the internal multiples. However, unlike JMI, it cannot update the propagation velocity model. Moreover, as shown by Garg & Verschuur (2017), JMI has the ability to explain the elastic nature of the input P-P reflection data, without explicitly imposing the full elastic wave equation. Therefore, the elastic amplitudes are preserved in the estimated local impulse responses above the reservoir without going 'full elastic' for the whole subsurface.
This paper is organised as follows: We first describes all the steps in JMI-res (Fig. 1a). Then, using a numerical elastic P-P reflection data, we demonstrate the JMI-res process to estimate the reservoir elastic parameters. Here, we will also focus on the effect of overburden internal multiples on the estimated impulse responses at the reservoir level and, consequently, on the reservoir elastic parameters. For this, we will estimate the impulse responses via standard redatuming, based on time reversal of surface seismic data (Beryhill 1984), where we don't account for any internal multiples in the data and apply the same target-oriented inversion scheme to it (Fig. 1b). We will show the higher resolution that we get in the estimated elastic reservoir parameters via JMI-res as compared to the standard redatuming route, courtesy of proper handling of overburden internal multiples in JMI-res. Note that this paper is an elaborate version of the work published by Garg et al. (2018).

R E S E RV O I R -O R I E N T E D J M I
In this section, we will review and explain all the steps of the JMI-res. Fig. 1(a) shows all the steps of the JMI-res workflow using the block diagram, while Fig. 1(b) shows all the steps for standard redatuming route.

JMI
JMI (Berkhout 2014c;Staal 2015;Verschuur et al. 2016) is a full waveform type data-fitting inversion process that estimates both reflectivity image and propagation velocity model, while estimating the up-and downgoing wavefields at all depth levels. The backbone of JMI is its forward modelling engine called full wavefield modelling (FWMod; Berkhout 2014a). Full wavefield modelling recursively and sequentially models all orders of scattering (primaries and surface/internal multiples) for the given velocity model and reflectivity image. Unlike finitedifference methods (Virieux 1986;Wang & Schuster 1996), FWMod does not attempt to solve the differential form of the wave equation directly, instead, it uses integral operators to compute wave propagation in the subsurface.
We consider a 2-D subsurface with M ∈ N + depth samples and N x ∈ N + lateral samples. For this subsurface, we consider a 2-D seismic data with N s ∈ N + sources, N r ∈ N + receivers and N ω ∈ N + frequency samples. Let P ± (x r , x s , ω; z n ) and Q ± (x r , x s , ω; z n ) be the 3-D scalar functions that represent the up-/downgoing wavefields at the depth level z n , while W ± (x i , x j , ω; z n ± 1, z n ) be the 3-D scalar function up-/downgoing propagation operators between z n and z n ± 1 in the frequency domain (see Fig. 2a). Here, x r and x s represent source and receiver locations, respectively, while x i and x j represent the lateral locations in the subsurface. The signs, + and -, represent the 'downgoing 'and 'upgoing', respectively. Let R ∪ (z n ), R ∩ (z n ) ∈ R Nr ×Nr and T + (z n ), T − (z n ) ∈ R Nr ×Nr be the matrices that represent the reflection operators and the transmission operators at depth z n , respectively, as shown in Fig. 2(a). Note, the variables after the semi-colon in all the notations indicate the variables that are fixed.
Using the vector-matrix notations introduced by Berkhout (1982), for a wavefield with single frequency ω f , f = 1, 2, ..., N ω and at one shot location x s, ξ , ξ = 1, 2, ..., N s , the wavefields P ± (x r ; x s, ξ , ω f , z n ) can be written as a vector P ± (z n ). Similarly, for the single frequency (ω f ), the propagation operators W ± (x i , x j ; ω f , z n ± 1, z n ) can be represented by the matrices W ± (z n ± 1 , z n ). A monochromatic downgoing wavefield P + (z n ) at depth level z n is reflected whenever there is a sharp discontinuity. This reflected PP wavefield, for the case where we can neglect the wave conversion, can be written using the reflection operator R ∪ (z n ): However, this upgoing wavefield is incomplete in the sense that it only contains the reflected energy and neglects the transmitted wavefield at depth level z n coming from below. Thus, to include the transmitted energy too, the total upgoing wavefield Q − (z n ) (see Fig. 2a) moving away from the depth level z n can be expressed as follows: where P − (z n ) is the upgoing incoming wavefield at depth level z n and T − (z n ) represents the transmission operator at the discontinuity. T ± (z n ) are given by: where the differential transmission operator δT ± (z n ) = 0 in case there is no contrast at z n . Using eq. (3), we can also write eq. (2) as: where the last two terms are the scattering terms at depth level z n . R ∪ (z n ) and δT − (z n ) can also be viewed as the backward and forward scattering operators, respectively (Berkhout 2014a). Thus, the scattering terms in eq. (4), forward (second term) and backward (third term), can be combined into one scattered wavefield δS − (z n ). δS − (z n ) can be viewed as the wavefield from the secondary sources at z n resulting in scattering of the incoming wavefields. Therefore, eq. (4) can be written in compact form as follows: Similarly, the total downgoing wavefield Q + (z n ) that leaves the depth level z n can be written as follows: The wavefields Q ± (z n ) at z n (see Fig. 2a) propagate to the neighbouring depth levels z n ± 1 via wavefield extrapolation: where W ± are the propagation operators defined in the space-frequency domain (see Berkhout 1982). The described scattering (eqs 2 and 6) and propagation processes (eqs 9 and 10) occurs at each depth level. The whole process can be summarized by the following equations (Staal 2015): (i) for the downgoing wavefields m = 1, 2, ..., M: (ii) for the upgoing wavefields m = 0, 1, ..., M − 1: where M and j indicates the maximum numbers of depth levels and the shot number, respectively. S + j (z 0 ) represents the jth downgoing source wavefield at the surface. The hybrid propagation operators W ± in eqs (11) and (12) are described as: The eqs (11) and (12) describe the scattering and propagation processes in FWMod at all depth levels in one round trip. By recursively applying eqs (11) and (12), more orders of scattering are included in the wavefield modelling, thereby making the up-/downgoing wavefields. This kind of forward modelling scheme is similar to Bremmer series, which is defined for the horizontally stratefied and plane wave propagation (Bremmer 1951). Moreover, for the case where wave conversion can be neglected and with small or moderate dipping reflectors, we can describe the scattering process in FWMod using only a single operator R via the following approximations: where R descibes the PP reflectivity. At the same time, W + (z n−1 , z n ) and W − (z n , z n−1 ) between the two consecutive depth levels are the transposed function of each other (see Berkhout 2014c; Staal 2015), thus only one W has to be considered. Thus, the FWMod modelled data is a function of a single scattering (R) and a single propagation operator (W). Here, R encode all the amplitude changes due to the scattering, while W explains the phase changes due to the propagation velocities. JMI (Berkhout 2014c;Staal 2015;Verschuur et al. 2016) iteratively minimizes the difference between the observed and the FWMod modelled data in a least-squares sense. Hence, the output of the JMI are the angle-dependent reflectivity image, free of internal multiple imprint and a propagation velocity model along with the up-/downgoing wavefields modelled at all depth levels. Note, as we are explaining the data sample-by-sample, JMI belongs to the full waveform inversion processes. Fig. 2(b) explains the implementation of JMI using a block diagram. The objective function for JMI is given as follows: where the, P − obs (z 0 ) is the collection of all recorded surface seismic shot records while P − mod (z 0 ) describes the FWMod modelled surface shot records as a function of reflectivity R and propagation velocity v. Moreover, Garg & Verschuur (2017) showed as the reflectivity operator (R) Downloaded from https://academic.oup.com/gji/article/221/1/115/5680486 by Bibliotheek TU Delft user on 08 February 2021 is not fully constrained by the forward modelling engine in JMI. Thus, it has the ability to explain the elastic P-P reflection data without explicitly imposing the elastic wave equation. Therefore, the modelled up-/downgoing wavefields contain the full extended coda with correct primary and internal multiples information and also have the preserved elastic amplitudes. Thus, they can be used to create the elastic amplitude preserved local impulse responses within the subsurface above the reservoir without going 'full elastic' for the whole subsurface.
We refer to Berkhout (2014c) and Staal (2015) for details on gradients estimations for scattering and propagation operators. Also, in the current implementation of JMI, we do not invert for angle-dependent reflectivity and velocity together as doing so leads to the risk of running into the null space due to overparametrization (Qu et al. 2018), but we rather apply JMI under an angle-independent or scalar reflectivity assumption. Thus, our intermediate outputs are the scalar reflectivity image and the propagation model. When we have a good enough velocity model, we apply JMI with fixed propagation velocity model, that is full wavefield migration (FWM, Berkhout 2014b; Davydenko & Verschuur 2017a), in angle-dependent reflectivity mode so as to explain the true elastic amplitudes present in the dataset and get the elastic downand upgoing wavefields. Note, recently Davydenko & Verschuur (2017b) presented an approach on how to estimate both angle-dependent reflectivity and propagation velocity within the JMI framework.

Proximity transformation
The elastic up-and downgoing wavefields at any subsurface depth level are equivalent to the up-/down decomposed Green's functions with real sources at the surface and the virtual receivers in the Earth's subsurface. Thus, the up-and downgoing wavefields at the target depth level from shot j can be related as follows: where P + j (z d ; z 0 ) is the downgoing wavefield at target depth level z d , which contains all orders of scattering (primary and internal multiples) and acts as the effective source wavefield for estimating the upgoing wavefield P − j (z d ; z 0 ). Thus, both wavefields contain the overburden coda. X(z d , z d ) represents the impulse responses from the area below depth level z d in the frequency domain. We can already anticipate, via eq. (19), that a more complex overburden above target depth level z d will result in an effective source wavefield containing more scattering, thus, providing a richer illumination in the area below the target depth. The impulse responses [X(z d , z d )] represents fully-redatumed reflection responses with both virtual sources and virtual receivers at depth level z d .
The impulse response can be estimated by solving the eq. (19) in the least-squares sense. However, given the ill-posedness associated with it, we use a non-quadratic mixed l 1 − l 2 norm as sparsity constraint and impose the principle of reciprocity (Knopoff & Gangi 1959) on the impulse responses in the time domain. These constraints help to reduce the ill-posedness of the problem and helps to reduce extrapolation anticausal artifacts. Thus, the constrained objective function is given by: The first term is the data residual term, while the second term is the constraint that promotes sparsity in the time domain impulse responses, where the local impulse responses for one time sample is represented byx matrix. The reciprocity constraint is applied as a hard constraint. Note, these estimated impulse responses (x) are dipole pressure source responses. Thus in strict terms,x is not equal to its transpose (x T ) unless we account for the obliquity factor that transform them to monopole responses, as shown by van der Neut et al. (2017). However, they are approximately equal, especially in kinematic terms. The 1 and σ x are the weight applied to the sparsity constraint and model parameters, respectively. The superscript T represents the transpose operation. We refer to this sparsity promoting least-squares inversion with reciprocity constraint as proximity transformation (da Costa et al. 2015;Garg & Verschuur 2016). This minimization problem is then solved using a gradient descent optimization approach with gradient and descent direction given as follows: where ∇J (k) and X (k) are the gradient and the descent direction for the kth iteration, respectively. The superscript H indicates the conjugate transpose of P + j .

Target-oriented full waveform inversion (FWI-res)
FWI-res is a target-oriented elastic wave-equation based inversion method developed by Gisolf & van der Berg (2012). In this inversion method, the media parameters V p , V s and ρ are redefined in terms of true elastic properties, notably the compressibility κ = 1 / K (with K being the bulk modulus), the shear compliance M = 1 / μ (with μ being the shear modulus) and ρ. Gisolf (2016) showed that this kind of parametrization, directly in terms of elastic properties (κ, M and ρ), is more closely related to the hydrocarbon saturation and porosity than impedances or velocities. Thus, is better suitable for the reservoir-oriented local inversion. Here, κ and M are related to V p , V s and ρ as follows: In FWI-res, the properties κ, M and ρ are reduced to contrasts against a background, instead of the property differences across the interface, according to: Thus, we solve for the contrasts functions (χ ) for the known space-variant background κ 0 , M 0 and ρ 0 . The 2-D implementation of this approach is done in the linear Radon (τ − p) domain with a locally laterally homogeneous medium assumption (i.e. local 1.5-D assumption) at the target area. Hence the input data, being the estimated local impulse responses via proximity transformation, are first sorted in the common midpoint-offset domain (CMP) and then transformed to the linear Radon domain before the inversion. Also, as we are solving the elastic wave equation, it can account for the non-linearities in the local impulse responses i.e. it can correctly account for the multiples and transmission effects generated within the reservoir area. Haffinger et al. (2016) showed this for the seismic dataset from the Sleipner gas field in the North Sea.
Here, we will briefly discuss the formulation and implementation of FWI-res. For detailed exposition and mathematical formulation, we suggest the readers to refer Gisolf & van der Berg (2012) and Rizutti & Gisolf (2017).
The FWI-res is based on the elastic scattering integral equation that is given as follows: where u and u inc represents the total wavefield and the incident wavefield respectively. G represents the Green function operator defined in the background medium and W[χ ] the matrix that depends on the contrast properties. A detailed explanation on eq. (25) and about its formulation is given in Yang et al. (2008) and Rizutti & Gisolf (2017). The key ingredient in the elastic scattering integral is the decomposition of the total wavefields in terms of the incident and the scattering wavefields where the elastic scattering integral is defined with the help of the object equation and the data equation (Rizutti & Gisolf 2017). The object equation relates the total wavefields with the incident wavefields and the contrasts functions where the incident wavefields are the solution of the elastic wave equation in the smoothly inhomogeneous media, calculated with the help of the WKJB approximation. On the other hand, the data equation relates the actual recorded data in the linear Radon domain to the total wavefields and the contrasts in the object domain. The inversion formulation in the FWI-res refers to simultaneously solving the object and the data equations of the forward elastic scattering formulation for the total wavefields and the contrast functions. In the inversion process, we do not find the exact solution to the object and data equations, rather we update the total wavefields and the contrasts as depicted in Fig. 3. In every iteration, we first estimate the contrasts functions (χ ) using the data equation given the best estimate of the total wavefields. As we assume fixed total wavefields while solving for the contrasts, it is reduced to a linear inversion. However, this linear inversion requires to be regularised using multiplicative regularisation based on minimal total variation (Abubakar et al. 2004). This is what constitutes the inner loop in Fig. 3. The total fields update is done by solving for the object equation together with the estimated contrasts from the inner loop and the previous estimate of the total wavefields using the optimised Neumann series, which is based on the Krylov subspace method (Kleinman & van der Berg 1991). This constitutes the outer loop in Fig. 3.

N U M E R I C A L E X A M P L E
We use a 2-D subsurface model with three high velocity anomalies (V p and V s ) and high densities (ρ) as overburden (Fig. 4) to generate the input elastic data. Since, our focus is to only explain the P-P reflection surface seismic data and not the converted waves, we smooth the  overburden lenses in the V s to minimize the converted waves in the elastic input data. The input P-P reflection data is generated via a finite difference implementation of the elastic wave equation (Thorbecke & Draganov 2011) using a Ricker wavelet of peak frequency 20 Hz and with sources and receivers at 20 m interval. Note, to decrease the converted waves at large offsets, we also apply high-angle filtering to the input data in the linear Radon domain. Fig. 5 shows the input data at three shot locations. The zoomed box in the Fig. 5 shows an overburden internal multiple with the apex at 0.9 s along with some converted waves events that obscure the P-P responses from the deeper section.
We apply JMI, that is, the first step in the JMI-res route (Fig. 1a) as explained in Section 2.1. Fig. 6 shows the inital velocity model and reflectivity image used for this inversion process. Fig. 7 shows the estimated velocity and structural image via JMI. We get a good estimate of the propagation velocity, which is further validated by the flat reflectivity angle gathers in Fig. 8, which were generated for Q.C. purpose. Fig. 9 shows the corresponding estimated down-and upgoing wavefields at the chosen target depth of z d = 680 m (see Fig. 12). We can see the correctly resolved overburden multiples in the downgoing wavefields [ P + (z 0 , z d ), Fig. 9a], which acts as an extended coda and Downloaded from https://academic.oup.com/gji/article/221/1/115/5680486 by Bibliotheek TU Delft user on 08 February 2021  provides extra illumination in the deeper section. The limited offsets in the upgoing wavefields [ P − (z 0 , z d ), Fig. 9b] are due to the limited receiver-aperture for the used acquisition geometry. Moreover, we can also see some unexplained converted waves artifacts in the upgoing wavefields [P − (z 0 , z d ), Fig. 9b].
The estimated down-and upgoing wavefields at target depth (z d = 680 m) are used to estimate the impulse responses, that is the virtual source-receiver data, via proximity transformation, which is the second step in the JMI-res process. In Fig. 10(a), we can see the estimated impulse responses with three clear reflections that correspond to the target area below the datum level. As a comparison, we also estimated the impulse responses via a conventional standard redatuming, based on time reversal of the recorded data. Note that we use the same velocity model estimated via JMI for standard redatuming (Fig. 1b). The impulse responses via standard redatuming (Fig. 10b) have the clear imprint of the unexplained overburden scattering on and around the second reflection event, even though we used the same good velocity model that was estimated in the JMI step. This comparison becomes even clearer when we transform the impulse responses to τ /p-CMP gathers for the localised target-oriented inversion, as shown in Fig. 11, with visible overburden internal multiple imprint. This imprint from the overburden internal multiples in standard redatuming case (Fig. 1b) has the potential to affect the local elastic parameters estimation.  Before applying the last step, that is, FWI-res, for elastic parameters estimation, we require a source wavelet estimate at the target depth. Thus, we assume that we have a well log at the central location of the model (Fig. 12), which is usually the case when using field data. Fig. 13 shows the estimated wavelet for both JMI-res and standard redatuming case at the well location on the target depth. The wavelet in the each case is estimated by least-squares matching of the estimated and the synthetic data modelled via Kennett modelling (Kennett 1983) using the well-log information. Note, we only consider the area between x = 1000 m to x = 2000 m as the target domain (Fig. 12) given it being the central region and will have the maximum illumination. Also, we only invert for κ and M as stable ρ inversion requires high-plane wave angles and P-S impulse responses too. Furthermore, the background κ, M and ρ for the whole target area are taken from the well log. Fig. 14 shows the estimated elastic parameters (κ and M) via FWI-res at 3 locations for both impulse responses estimated via JMI-res and standard redatuming. We clearly see the effect of the unresolved overburden internal multiples in the standard redatuming case with severe effects on the estimated M, especially around the second interface. At the same time, we also make the 2-D section of the inverted local elastic parameters by applying FWI-res at all the locations in the target area. Fig. 15 shows the comparison of the estimated κ and M for the whole target area. We clearly see the high-resolution results we get in elastic parameters estimated via JMI-res compared to the standard redatuming route. We get better elastic parameters inversion values and sharp layer boundaries in the JMI-res case due to the more accurate  input impulse responses, courtesy of correctly explained multiple scattering energy in the JMI-res redatuming step. The unresolved internal multiple (also pointed out in the input data in Fig. 5) in the standard redatuming case affects the elastic parameters estimation as observed from the unrecoved elastic parameter values and no clear layer boundaries demarcations, especially at the second interface for M.

D I S C U S S I O N
Although the aim of the paper is to explain all the steps associated with reservoir-oriented JMI (JMI-res) and how it provides high-resolution local elastic parameters, there are certain areas that are not explored. We neglected the converted wave modes that partly end up in the estimated impulse responses. Berkhout (2014c) proposed to extend the current JMI process for explaining the PS and SP converted data. Note, however, the estimated impulse responses in JMI-res have the elastic characteristics preserved in the PP reflections and are accurate enough for the target-oriented full waveform inversion (FWI-res). Another aspect that is not discussed here is the 3-D extension of the JMI-res. Recently, Davydenko & Verschuur (2018) proposed a way to estimate the local impulse responses in a full wavefield migration/JMI framework for the 3-D case.
However, we have shown the ability of JMI-res to estimate the accurate reservoir impulse responses and associated reservoir elastic parameters from surface seismic elastic data for a complex overburden scenario, courtesy of correctly resolving the overburden internal multiples in the JMI-based redatuming step. In this sense the output of a JMI-based redatuming step is similar to that of the Marchenko redatuming scheme proposed by Ravasi (2017) as both can handle the overburden effects while estimating the reservoir impulse responses. Note that the Marchenko method is fully data driven, while the JMI method relies on a physical model (via propagation and reflection operators). As such, Marchenko will be sensitive to data sampling, while JMI can handle more sparsely sampled data (Garg & Verschuur 2016). Both methods are unlike any standard redatuming procedure where we don't account for overburden multiple scattering and transmission effects.
As a next step, all the aspects of the proposed JMI-res process need to be verified on field data. In this paper we introduced the methodology and showed its feasibility on synthetic data. We leave the field data application for future research.

C O N C L U S I O N S
In this paper, we presented the reservoir-oriented JMI (JMI-res) approach that first estimates the accurate local impulse responses, free of overburden multiples and transmission imprint, and then uses them as input for a target-oriented full waveform inversion process to estimate the high-resolution local elastic parameters. The comparison of JMI-res results with the standard redatuming route results for a synthetic data case with a strong-reflecting overburden showed the effects of the overburden internal multiples in the estimated local impulse responses and how it ultimately affects the localised inversion process for the estimation of reservoir elastic parameters if they are not properly accounted for. This is besides the propagation velocity estimation being an integral part of JMI-res. Moreover in JMI-res, we estimated the elastic characteristics preserved local impulse responses at the target level without going 'full elastic' for the whole subsurface. Thus, this kind of target-oriented approach avoids applying computationally expensive and non-linear elastic-FWI for the whole subsurface in the single-component PP reflection data case, and allows an elastic-FWI application only restricted to the reservoir target domain.