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

SUMMARY

Traditionally, reservoir elastic parameters inversion suffers from the overburden multiple scattering and transmission imprint in the local input data used for the target-oriented inversion. In this paper, we present a full-wavefield approach, called reservoir-oriented joint migration inversion (JMI-res), to estimate the high-resolution reservoir elastic parameters from surface seismic data. As a first 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.

Image processing, Inverse theory, Tomography, Waveform inversion, Wave scattering and diffraction

1 INTRODUCTION

One of the final goals in exploration seismics, apart from delineating the prospective reservoir structures, is the characterization of reservoir properties (Chopra & Marfurt 2007; Doyen 2007) in terms of reservoir elastic parameters. These reservoir elastic parameters provide much more insight about the prospective reservoir as we get a quantitative image in terms of actual elastic parameters that can be linked to reservoir properties such as porosity and lithology (Eberhart-Phillips et al. 1989; Feng et al. 2017).

Over the years, geophysicists had estimated the reservoir elastic parameters by applying Zoeppritz-based AVO/AVA inversion methods (Castagna & Backus 1993; Buland & Omre 2003) to pre-stack image angle gathers (de Bruin et al. 1991; Sava & Fomel 2003) obtained after migration (Baysal et al. 1983; Bednar 2005). Recently, there has also been a growing interest in the target-oriented full waveform inversion (Gisolf & van der Berg 2012; Biondi et al. 2018; Garg et al. 2018), where first the virtual source–receiver reflection data is estimated at the depth above the potential reservoir via redatuming (Beryhill 1984; Wapenaar et al. 1992; Schuster & Zhou 2006; Garg & Verschuur 2016; Ravasi 2017) and then elastic full waveform inversion (Vigh & Starr 2008; Virieux & Operto 2009; Haffinger 2013) is applied only to this fully redatumed data in a local inversion domain.

The local input data for either the Zoeppritz-based or target-oriented full waveform inversion approach should be accurate in order to get the reliable elastic parameters. First of all, they should have preserved elastic reflection characteristics, which is essential to get the correct elastic properties. Secondly, they should be free of any spurious events related to multiples or transmission energy imprint from the overburden. However, most migration methods (Baysal et al. 1983; Bednar 2005) and model-based or correlation-based redatuming methods (Beryhill 1984; Wapenaar et al. 1992; Schuster & Zhou 2006) fail to correctly explain the internal multiples generated from the overburden. Thus, in the case of a strongly scattering overburden, the events related to multiple (internal) reflection energy from the overburden obscure and overlay the local input data creating both spurious events and affecting the amplitude characteristics. This in turn affects the local inversion (AVA/AVO inversion or target-oriented FWI) and resolution of estimated reservoir elastic parameters.

In this paper, we present the reservoir-oriented joint migration inversion (JMI-res) that first generates the fully redatumed virtual source–receiver data at the reservoir level, free of overburden multiple scattering imprint, using JMI (Berkhout 2014c; Staal 2015; Verschuur et al. 2016) and then applies a target-oriented local inversion process to estimate the reservoir elastic parameters. In JMI-res, we first estimate elastic 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 (Verschuur et al. 2016). 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).

Figure 1.

Flowchart depicting (a) the JMI-res and (b) the standard redatuming route.

Open in new tab Download slide

2 RESERVOIR-ORIENTED JMI

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.

2.1 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 finite-difference 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 \in \mathbb {N^{+}}$| depth samples and |$N_{x} \in \mathbb {N^{+}}$| lateral samples. For this subsurface, we consider a 2-D seismic data with |$N_{s} \in \mathbb {N^{+}}$| sources, |$N_{r} \in \mathbb {N^{+}}$| receivers and |$N_{\omega } \in \mathbb {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 |$\mathbf {R}^{\cup }(z_n), \mathbf {R}^{\cap }(z_n) \in \mathbb {R}^{N_{r} \times N_{r}}$| and |$\mathbf {T}^{+}(z_n), \mathbf {T}^{-}(z_n) \in \mathbb {R}^{N_{r} \times N_{r}}$| 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.

Figure 2.

(a) Schematic representation of reflection and transmission at a single depth level along with the wavefield extrapolation, to the neighbouring levels. (b) Flowchart depicting the JMI.

Open in new tab Download slide

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 |$\skew{4}\vec{P}^{\pm }(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 |$\skew{4}\vec{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 |$\mathbf {R}^{\cup }(z_n)$|⁠:

$$\begin{eqnarray*} \skew{4}\vec{Q}^{-}(z_n) = \mathbf {R}^{\cup }(z_n)\skew{4}\vec{P}^{+}(z_n). \end{eqnarray*}$$

(1)

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 |$\skew{4}\vec{Q}^{-}(z_n)$| (see Fig. 2a) moving away from the depth level z_n can be expressed as follows:

$$\begin{eqnarray*} \skew{4}\vec{Q}^{-}(z_n) = \mathbf {T}^{-}(z_n)\skew{4}\vec{P}^{-}(z_n) + \mathbf {R}^{\cup }(z_n)\skew{4}\vec{P}^{+}(z_n), \end{eqnarray*}$$

(2)

where |$\skew{4}\vec{P}^{-}(z_n)$| is the upgoing incoming wavefield at depth level z_n and |$\mathbf {T}^{-}(z_n)$| represents the transmission operator at the discontinuity. |$\mathbf {T}^{\pm }(z_n)$| are given by:

$$\begin{eqnarray*} \mathbf {T}^{\pm }(z_n) = \mathbf {I} + \delta \mathbf {T}^{\pm }(z_n), \end{eqnarray*}$$

(3)

where the differential transmission operator |$\delta \mathbf {T}^{\pm }(z_n) = 0$| in case there is no contrast at z_n. Using eq. (3), we can also write eq. (2) as:

$$\begin{eqnarray*} \skew{4}\vec{Q}^{-}(z_n) = \skew{4}\vec{P}^{-}(z_n) + \delta \mathbf {T}^{-}(z_n)\skew{4}\vec{P}^{-}(z_n) + \mathbf {R}^{\cup }(z_n)\skew{4}\vec{P}^{+}(z_n), \end{eqnarray*}$$

(4)

where the last two terms are the scattering terms at depth level z_n. |$\mathbf {R}^{\cup }(z_n)$| and |$\delta \mathbf {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 |$\delta \mathbf {S}^{-}(z_n)$|⁠. |$\delta \mathbf {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:

$$\begin{eqnarray*} \skew{4}\vec{Q}^{-}(z_n) = \skew{4}\vec{P}^{-}(z_n) + \delta \mathbf {S}^{-}(z_n). \end{eqnarray*}$$

(5)

Similarly, the total downgoing wavefield |$\skew{4}\vec{Q}^{+}(z_n)$| that leaves the depth level z_n can be written as follows:

$$\begin{eqnarray*} \skew{4}\vec{Q}^{+}(z_n) &=& \mathbf {T}^{+}(z_n)\skew{4}\vec{P}^{+}(z_n) + \mathbf {R}^{\cap }(z_n)\skew{4}\vec{P}^{-}(z_n) \\ &=& \skew{4}\vec{P}^{+}(z_n) + \delta \mathbf {T}^{+}(z_n)\skew{4}\vec{P}^{+}(z_n) + \mathbf {R}^{\cap }(z_n)\skew{4}\vec{P}^{-}(z_n) \\ &=& \skew{4}\vec{P}^{+}(z_n) + \delta \mathbf {S}^{+}(z_n). \end{eqnarray*}$$

(6),(7),(8)

The wavefields |$\skew{4}\vec{Q}^{\pm }(z_n)$| at z_n (see Fig. 2a) propagate to the neighbouring depth levels z_{n ± 1} via wavefield extrapolation:

$$\begin{eqnarray*} \skew{4}\vec{P}^{-}(z_{n-1}) = \mathbf {W}^{-}(z_{n-1},z_n)\skew{4}\vec{Q}^{-}(z_n), \end{eqnarray*}$$

(9)

$$\begin{eqnarray*} \skew{4}\vec{P}^{+}(z_{n+1}) = \mathbf {W}^{+}(z_{n+1},z_n)\skew{4}\vec{Q}^{+}(z_n), \end{eqnarray*}$$

(10)

where |$\mathbf {W}^{\pm }$| 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):

for the downgoing wavefields m = 1, 2, ..., M:
$$\begin{eqnarray*} \skew{4}\vec{P}_{j}^{+}(z_m;z_0) = \mathbf {\underline {W}}^{+}(z_m,z_0)\skew{4}\vec{S}_{j}^{+}(z_0)+\sum _{n=0}^{m-1}\mathbf {\underline {W}}^{+}(z_m,z_n)\mathbf {R}^{\cap }(z_n)\skew{4}\vec{P}^{-}_{j}(z_n;z_0), \end{eqnarray*}$$
(11)
for the upgoing wavefields m = 0, 1, ..., M − 1:
$$\begin{eqnarray*} \skew{4}\vec{P}_{j}^{-}(z_m;z_0) = \mathbf {\underline {W}}^{-}(z_m,z_M)\skew{4}\vec{P}_{j}^{-}(z_M;z_0)+\sum _{n=m+1}^{M}\mathbf {\underline {W}}^{-}(z_m,z_n)\mathbf {R}^{\cup }(z_n)\skew{4}\vec{P}^{+}_{j}(z_n;z_0), \end{eqnarray*}$$
(12)

where M and j indicates the maximum numbers of depth levels and the shot number, respectively. |$\skew{4}\vec{S}^{+}_{j}(z_0)$| represents the jth downgoing source wavefield at the surface. The hybrid propagation operators |$\mathbf {\underline {W}}^{\pm }$| in eqs (11) and (12) are described as:

$$\begin{eqnarray*} \mathbf {\underline {W}}^{+}(z_m,z_n) = \mathbf {W}^{+}(z_{m},z_{m-1})\prod _{l=m-1}^{n+1}\mathbf {T}^{+}(z_l)\mathbf {W}^{+}(z_{l},z_{l-1}), \end{eqnarray*}$$

(13)

$$\begin{eqnarray*} \mathbf {\underline {W}}^{-}(z_m,z_n) = \mathbf {W}^{-}(z_{m},z_{m+1})\prod _{l=m+1}^{n-1}\mathbf {T}^{-}(z_l)\mathbf {W}^{-}(z_{l},z_{l+1}). \end{eqnarray*}$$

(14)

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 |$\mathbf {R}$| via the following approximations:

$$\begin{eqnarray*} \mathbf {\delta T}^{+} = \mathbf {R}^{\cup }, \end{eqnarray*}$$

(15)

$$\begin{eqnarray*} \mathbf {\delta T}^{-} = \mathbf {R}^{\cap }, \end{eqnarray*}$$

(16)

$$\begin{eqnarray*} \mathbf {R} = \mathbf {R}^{\cup } = -\mathbf {R}^{\cap }, \end{eqnarray*}$$

(17)

where |$\mathbf {R}$| descibes the PP reflectivity. At the same time, |$\mathbf {W}^{+}(z_{n-1},z_n)$| and |$\mathbf {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 |$\mathbf {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, |$\mathbf {R}$| encode all the amplitude changes due to the scattering, while |$\mathbf {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:

$$\begin{eqnarray*} J = \sum _{\omega }\left|\left|\mathbf {P}_{obs}^{-}(z_0)-\mathbf {P}_{mod}^{-}(z_0)\right|\right|_{2}^{2}, \end{eqnarray*}$$

(18)

where the, |${\bf P}_{obs}^{-}(z_0)$| is the collection of all recorded surface seismic shot records while |${\bf 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 (⁠|$\mathbf {R}$|⁠) 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 down- and 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.

2.2 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:

$$\begin{eqnarray*} \skew{4}\vec{P}_{j}^{-}(z_d;z_0) = \mathbf {X}(z_d,z_d)\skew{4}\vec{P}_{j}^{+}(z_d;z_0), \end{eqnarray*}$$

(19)

where |$\skew{4}\vec{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 |$\skew{4}\vec{P}_{j}^{-}(z_d;z_0)$|⁠. Thus, both wavefields contain the overburden coda. |$\mathbf {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 [|$\mathbf {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₁ − l₂ 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:

$$\begin{eqnarray*} J & = & \sum _{j}\sum _{\omega }\left|\left|\skew{4}\vec{P}^{-}_{j}(z_d;z_0)-\mathbf {\hat{X}}(z_d,z_d)\skew{4}\vec{P}_{j}^{+}(z_d;z_0)\right|\right|_{2}^{2} + \epsilon _1^2\sum _{t}\sum _{l,m}\sqrt{1+\frac{\mathbf {\hat{x}}^{2}_{lm}(z_d,z_d)}{{\sigma }_{x}^{2}}},\nonumber \\ &\rm{{ subject \, to}} & \hat{{\mathbf {x}}} = \hat{{\mathbf {x}}}^T. \end{eqnarray*}$$

(20)

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 by |$\mathbf {\hat{x}}$| matrix. The reciprocity constraint is applied as a hard constraint. Note, these estimated impulse responses (⁠|$\mathbf {\hat{x}}$|⁠) are dipole pressure source responses. Thus in strict terms, |$\mathbf {\hat{x}}$| is not equal to its transpose (⁠|$\mathbf {\hat{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 ϵ₁ 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:

$$\begin{eqnarray*} \nabla J^{(k)} & = -&\sum _{j}\sum _{\omega } \left\lbrace \left[\skew{4}\vec{P}^{-}_{j}(z_d;z_0)-\mathbf {\hat{X}}^{(k)}(z_d,z_d)\skew{4}\vec{P}^{+}_{j}(z_d;z_0)\right]\left[\skew{4}\vec{P}^{+}_{j}(z_d;z_0)\right]^{H}\right\rbrace \nonumber \\ & + & \epsilon _1^2\sum _{t}\sum _{l,m}\frac{\mathbf {\hat{x}}^{k}_{lm}(z_d,z_d)}{\sqrt{{\sigma }_{x}^{4}+{\sigma }_{x}^{2}\mathbf {\hat{x}}^{(k)2}_{lm}}}, {{ subject \, to }} \hat{{\mathbf {x}}} = \hat{{\mathbf {x}}}^T, \end{eqnarray*}$$

(21)

$$\begin{eqnarray*} \Delta \mathbf {\hat{X}}^{(k)}(z_d,z_d) = -\nabla J^{(k)}, \end{eqnarray*}$$

(22)

where ∇J^(k) and |$\Delta \mathbf {\hat{X}}^{(k)}$| are the gradient and the descent direction for the kth iteration, respectively. The superscript H indicates the conjugate transpose of |$\skew{4}\vec{P}^{+}_{j}$|⁠.

2.3 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 κ = ¹/_K (with K being the bulk modulus), the shear compliance M = ¹/_μ (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:

$$\begin{eqnarray*} \kappa = \frac{1}{\rho ({V_p}^2 - {}^{4}\!/_{3}{V_s}^2)} \text{, } M = \frac{1}{\rho {V_s}^2}. \end{eqnarray*}$$

(23)

In FWI-res, the properties κ, M and ρ are reduced to contrasts against a background, instead of the property differences across the interface, according to:

$$\begin{eqnarray*} {\chi }_{\kappa } = \frac{\kappa - {\kappa }_{b}}{{\kappa }_{b}}, {\chi }_{M} = \frac{M - {M}_{b}}{M_{b}}, {\chi }_{\rho } = \frac{\rho - {\rho }_{b}}{{\rho }_{b}}. \end{eqnarray*}$$

(24)

Thus, we solve for the contrasts functions (χ) for the known space-variant background κ₀, M₀ and ρ₀. 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:

$$\begin{eqnarray*} \mathbf {u} = \mathbf {u}^{inc} + \mathcal {G}\mathcal {W}[\chi ]\mathbf {u}, \end{eqnarray*}$$

(25)

where |$\mathbf {u}$| and |$\mathbf {u}^{inc}$| represents the total wavefield and the incident wavefield respectively. |$\mathcal {G}$| represents the Green function operator defined in the background medium and |$\mathcal {W}[\chi ]$| 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.

Figure 3.

Flowchart depicting the FWI-res [after Gisolf & van der Berg (2012)].

Open in new tab Download slide

3 NUMERICAL EXAMPLE

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.

Figure 4.

Subsurface (a) V_p (b) ρ (c) V_s models used to generate the elastic FD data.

Open in new tab Download slide

Figure 5.

Surface shot records at three location depicting the internal multiples imprint in the deeper section, as indicated by the yellow arrows.

Open in new tab Download slide

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.

Figure 6.

(a) Initial reflectivity image and (b) intial velocity image used for JMI step.

Open in new tab Download slide

Figure 7.

Estimated (a) structural image and (b) P-wave velocity model in the JMI step.

Open in new tab Download slide

Figure 8.

Image angle gathers at three lateral locations obtained after the JMI step.

Open in new tab Download slide

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 [|$\skew{4}\vec{P}^{+}(z_0,z_d)$|⁠, Fig. 9a], which acts as an extended coda and provides extra illumination in the deeper section. The limited offsets in the upgoing wavefields [|$\skew{4}\vec{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₀, z_d), Fig. 9b].

Figure 9.

Estimated (a) down- and (b) upgoing wavefields for three lateral locations at the target depth of z_d = 680 m in the JMI step.

Open in new tab Download slide

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.

Figure 10.

Comparison of impulse responses (X) at three lateral locations at the target depth of z_d = 680 m estimated in (a) the JMI-res route and (b) via standard redatuming.

Open in new tab Download slide

Figure 11.

Comparison of CMP-τ/p domain impulse responses (X) at three lateral locations at the target depth of z_d = 680 m estimated (a) through the JMI-res process and (b) via standard redatuming.

Open in new tab Download slide

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.

Figure 12.

Diagram depicting the well location at x = 1500 m, target depth and target domain for the localized FWI-res.

Open in new tab Download slide

Figure 13.

Estimated wavelet for (a) the JMI-res and (b) the standard redatumed impulse response at the well location on the target depth, respectively.

Open in new tab Download slide

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.

Figure 14.

True (red), estimated (blue) and background (green) elastic parameters (κ, M) obtained via (a)–(c) the JMI-res route and (d)–(f) via the standard redatuming route. (a), (c) Results at x = 1100 m, (b), (e) at x = 1500 m and (c), (f) at x = 1900 m.

Open in new tab Download slide

Figure 15.

(a) True and (b), (c) estimated elastic parameters (κ, M) at all locations in the target domain (b) via the JMI-res route and (c) via the standard redatuming route.

Open in new tab Download slide

4 DISCUSSION

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.

5 CONCLUSIONS

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.

ACKNOWLEDGEMENTS

The authors would like to thank Siddharth Sharma for all the fruitful discussions. The authors also thank the sponsors of the Delphi consortium for their support.

REFERENCES

Abubakar

van der Berg

P.M.

Habashy

T.M.

Braunisch

2004

A multiplicative regularization approach for deblurring problems

IEEE Trans. Image Process.

(

1524

–

1532

Month:	Total Views:
December 2019	51
January 2020	48
February 2020	200
March 2020	106
April 2020	91
May 2020	47
June 2020	32
July 2020	41
August 2020	21
September 2020	31
October 2020	13
November 2020	11
December 2020	24
January 2021	17
February 2021	19
March 2021	26
April 2021	21
May 2021	14
June 2021	22
July 2021	20
August 2021	19
September 2021	11
October 2021	13
November 2021	12
December 2021	15
January 2022	13
February 2022	25
March 2022	12
April 2022	13
May 2022	8
June 2022	11
July 2022	18
August 2022	27
September 2022	17
October 2022	19
November 2022	14
December 2022	19
January 2023	17
February 2023	25
March 2023	12

Article Contents

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

SUMMARY

1 INTRODUCTION

2 RESERVOIR-ORIENTED JMI

2.1 JMI

2.2 Proximity transformation

2.3 Target-oriented full waveform inversion (FWI-res)

3 NUMERICAL EXAMPLE

4 DISCUSSION

5 CONCLUSIONS

ACKNOWLEDGEMENTS

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only