Elastic Full-Waveform Inversion Using Migration-Based Depth Reflector Representation in the Data Domain

Abstract

The depth velocity model is a critical element for providing seismic data processing success, as it is responsible for the times of waves’ propagation and, therefore, prescribes the location of geological objects in the resulting seismic images. Constructing a deep velocity model is the most time-consuming part of the entire seismic data processing, which usually requires interactive human intervention. This article introduces the consistently numerical method for reconstructing a depth velocity model based on the modified version of the elastic Full Waveform Inversion (FWI). The specific feature of this approach to FWI is the decomposition of the space of admissible velocity models into subspaces of propagator (macro velocity) and reflector components. In turn, the latter transforms to the data space reflectivity on the base of migration transformation. Finally, we perform minimisation in two different spaces: (1) Macro velocity as a smooth spatial function; (2) Migration transforms data space reflectivity to the spatial reflectivity. We present numerical experiments confirming less sensitiveness of the modified version of FWI to the lack of the low time frequencies in the data acquired. In our computations, we use synthetic data with valuable time frequencies from 5 Hz.

1. Introduction

The approach for recovery of the wave propagation velocity named Full Waveform Inversion (FWI) proposed back in the 70–80s of the last century [1,2,3]. It reduces the inverse dynamic seismic problem to minimise data misfit functional characterising the difference between the recorded and synthetic data. To compute the data deviation, the people usually use the L₂ norm; that is, the FWI (Full Waveform Inversion) is inherently nothing more than using a nonlinear least-squares technique.

One of the main difficulties in applying FWI is recovering the macro velocity component in the absence of low temporal frequencies in the data.

Traditionally, the resolution of FWI-based inverse dynamic seismic problem has the following two main two stages:

(1)

Determination of a macro velocity model (propagator) that describes the smooth distribution of the wave propagation velocity in the medium and determines the wave travel times.

(2)

Recovery of a sharply changing component (reflectors) corresponding to the medium’s interfaces, geometrical localisation, and reflection/transmission coefficients.

Reconstruction of the macro velocity component in the absence of low temporal frequencies in the data is one of the main stumbling blocks for the successful application of FWI for processing seismic data in a realistic setting ([4,5]). As noted in [6], this is due to the complex structure of the data misfit functional. Numerical experiments proved that this function is not convex concerning the macro velocity component in the absence of low temporal frequencies [7]. One way to overcome this difficulty based on the initial use of the lowest available time frequencies with their gradual increase was proposed and implemented in [6,8]. It follows a more stable recovery of the macro velocity component. However, there are still very significant restrictions on low temporal frequencies in the recorded data spectrum [6,7]. Another drawback of this approach is its implementation in the time domain, which lead to the necessity to perform special filters transforming input data to the ones within the prescribed time-frequency interval. This increase claims additional computational resources.

Therefore, the next step in developing FWI for macro velocity reconstruction for the realistic acquisitions was its implementation for the frequency-time domain made by G.Pratt in [9]. In this domain appear no problems when implementing gradual increasing of the time frequencies used for the FWI. Still, one should possess effective 3D solver for the frequency domain formulation of acoustic/elastic wave equations (see, for example, [10]). Below we also do minimisation in the time-frequency domain and should use this kind of the solver. Our previous study of the FWI in the time-frequency domain developed an effective iterative solver for the inversion of both acoustic and elastic wave fields [11,12]. Its main ingredient is a novel preconditioner, which accelerates the convergence of the iteration essentially. We have developed and justified a method to effectively invert the preconditioner based on the 2D fast Fourier transform and solve a system of linear algebraic equations with a banded matrix. We parallelise the solver using the conventional hybrid parallelisation (MPI in conjunction with OpenMP). On this way, we have got the high scalability for the widespread 3D SEG/EAGE overthrust model. We find that our method has a high potential for low-frequency simulations in land models with moderate lateral variations and vertical variations (see details in [13]).

This article develops an original technique for reconstructing the macro velocity component, based on modifying the objective function. This modification uses the decomposition of the space of the admissible velocity models into two subspaces:

• a subspace of smoothly varying velocity functions that do not introduce any significant deviations in the direction of propagation of seismic energy excited the sources on the acquisition, but determine the wave propagation time;
• a subspace of abruptly changing functions describing the reflecting boundaries’ location within the macro velocity component.

This approach’s roots are the MBTT (Migration Based Travel Time) technique introduced by G.Chavent with coauthors in the papers [14,15]. The modified version of FWI is again in the time domain, which brings some technical and computational issues in the implementation due to the gradient’s complicated computations for the macrovelocity component.

We modified this approach to the inversion of acoustic waves in the frequency-time domain in the paper [16] and successfully verified it by applying to Marmousi data for realistic acquisition. We do the next step towards the real seismic data processing based on the modified FWI version by developing it for inversion of elastic data for realistic land seismic acquisition.

2. Method and Theory

2.1. Data Space Reflectivity FWI Reformulation

Let us treat the dynamical inverse problem of the seismic waves’ propagation as some nonlinear operator equation:

$$F\left(m\right)=d,$$

with nonlinear forward map $F:M\to D$ transforming the model m to the data d.

Following acoustic MBTT formulation [14,15] we decompose both P- and S-waves propagation velocities in smooth propagator $p$ and quick changing in space reflector $r$. As we deal with multi-parameter inverse problem searching for both $V_p$ and $V_s$ velocity models, is a vector $p=\left(p{V}_p,p{V}_s\right)$, where $p{V}_p$—P-velocity propagator and $p{V}_s$—S-velocity propagator.

The next step of the data misfit functional modification relies on the fact that with the known macro velocity model the locations of the interfaces within the medium are set due to the migration operator’s action to the data d. These data should undergo some preprocessing procedures like multiple attenuations, amplitude correction and others [17]. In other words, there is a connection between reflectors and some transformation of the data in time/frequency-time domain: r = M (p) < s>. Note that here we introduce the new unknown s— Data Space Reflectivity (DSR). On this base, we represent wave propagation velocities as follows:

$$V_p=p{V}_p+M_p\left(p{V}_p,p{V}_s\right)\langle s\rangle$$

$$V_s=p{V}_s+M_s\left(p{V}_p,p{V}_s\right)\langle s\rangle .$$

where $p{V}_p, p{V}_s$ are P- and S-waves propagators respectively, while $M_p, M_s$are P- and S-waves migration operators.

The critical moment of this decomposition is propagator-reflector interrelation $r=M\left(p\right)\langle s\rangle$ with operator $M\left(p\right)$ being some kind of true-amplitude prestack migration/linearised inversion. In the paper we use like these operators reweighted version of the operators on the base of the adjoint state technique [Plessix]:

$$M\left(pVp,pVs\right)s=W\circ Re\left\{D{F}^{*}\langle s\rangle \right\}$$

here $DF$ is Frechet derivative of the forward modelling operator $F$, * denotes the adjoint operator, and $W$ is some linear reweighting operator providing true amplitudes imaging/migration [18,19,20,21].

Finally, we come to the following modified data misfit functional:

$$E\left(Vp,Vs,s\right)=\parallel F\left(\left(pVp,pVs\right)+M\left(pVp,pVs\right)\langle s\rangle \right)-d{\parallel }_D^2.$$

Minimisation with respect to propagator p and data-space reflectivity s are independent and doing by turn. We start with admitting $s=d$ and search for some intermediate value of propagator $p$. After stabilisation of this process, the search is switched to data-space reflectivity $s$ and so on.

We implemented the search for the propagator using Gauss-Newton method [9]. The dimension of propagator space is small enough concerning the size of the whole model space $M$ and therefore we can calculate partial SVD (only highest singular vectors which correspond the highest singular values) of the approximate Hessian ${\mathcal{H}}_a$ on the first iteration using matrix-free methods (such as Krylov-subspace based methods, for numerical implementation, see [11,12,13]). On the next step, we use r-pseudoinverse of this operator [22,23,24] as a preconditioner in the Conjugate-Gradient method during propagator minimisation:

$$p_{k+1}=p_k+{\mu }_k{S}_k,S_0={\nabla }_0,$$

$$S_k=-P{\nabla }_k-\frac{{\langle P{\nabla }_k,{\nabla }_k-{\nabla }_{k-1}\rangle }_M}{{\langle P{\nabla }_{k-1},{\nabla }_{k-1}\rangle }_M}{S}_{k-1},$$

where $p_k=\left(p{V}_p,p{V}_s\right)$—propagator on k-th iteration, $P={\left({\mathcal{H}}_a\right)}_r^{†}$—pseudoinverse of the approximate Hessian. The gradient ${\nabla }_k$ for the propagator is calculated as follows (see Appendix A):

$${\nabla }_k=Re\{D{F}^{*}\delta d_k+{(D^{2}F\left(p_k\right)W^{*}D{F}^{*}\delta d_k)}^{*}s\},$$

here $\delta d_k=F\left(m_k\right)-d$ is data residual on current iteration, $DF$—a first derivative of the forward map calculated at the point $m_k=p_k+M\left(p_k\right)s$ and $D^{2}F$—second derivative calculated at the point $p_k$. Incorporation into inversion process Hessian ${\mathcal{H}}_a$ helps to mitigate cross-talk between $V_p$ and $V_s$ models. Minimisation of the data-space reflectivity $s$ is doing via steepest descent method. The corresponding gradient ${\nabla }_s$ is calculated as follows:

$${\nabla }_s=M^{*}\left(p\right)D{F}^{*}\delta d,$$

where $M^{*}\left(p\right)$—is the adjoint to the migration operator, also known as the demigration.

2.2. The Validation of Propagator/Reflector Decomposition

The first question to be answered is the possibility of representing the depth reflector $r=\left(r{V}_p,r{V}_s\right)$ through the migration operator with current propagator $p=\left(p{V}_p,p{V}_s\right)$. To assess this representation’s feasibility, we consider the vertically inhomogeneous layered medium (see Figure 1). To evaluate the cross-talks between $V_p$ and $V_s$ models during minimisation of modified misfit function $E\left(p,s\right)$ for data-space reflectivity the targets horizons for $V_p$ and $V_s$ we choose models such that they do not match in the upper part. Input data are synthesised for the set of ten uniformly sampled frequencies in the range 5–30 Hz. The acquisition system has 25 vertical force point sources and 100 2C-receivers located at depth 10 m with a lateral spacing of 80 m and 20 m. For the sake of simplicity in the definition of migration operator, we use the identity weight matrix W=Id, that is, we do not use true-amplitude reweighting in this experiment. Propagators $p{V}_p$ and $p{V}_s$ (see Figure 1, red plots) provide the information on macro velocity—we are concentrated only on the reconstruction of space reflectivity.

During the minimisation process, we performed 30 nonlinear Conjugate Gradients iterations for data-space reflectivity unknown. The results of inversion one can see in Figure 1 and Figure 2. As one can observe there are no cross-talk effects between $V_p$ and $V_s$ models and the depth reflectors are sitting in the correct positions.

2.3. Choice of Depth Reflectivity Parametrisation

Another critical point is the parameterisation of the depth reflectivity:

Which parameterisation of elastic models fits better for computation of the spatial reflectivity in the context of the inverse problem?

One of the most widespread approaches to choosing optimal parameterisation is to analyse radiation pattern [25,26], but we prefer justification of choice by SVD-analysis of the corresponding linear operators.

Consider two different parameterisations for the depth reflectivity: $r_1=\left(r{V}_p,r{V}_s\right)$ and $r_2=\left(r{I}_p,r{I}_s\right)$, where $I_p=\rho V_p$, $I_s=\rho V_s$—elastic impedances, $\rho$ is a density. For a fixed acquisition geometry and reference model $m_0$ one may consider two linearised forward maps $L_1=\frac{\partial F}{\partial m_1}$ and $L_2=\frac{\partial F}{\partial m_2}$—first Frechet derivatives of the forward map $F$ concerning the different elastic model parametrisations, evaluated at the point $m_0$. Following [16] the comparative resolution analysis of singular value decomposition (SVD) of two linear maps $L_1$ and $L_2$ will answer the question. The numerically calculated singular spectrums we present in Figure 3. As one can see, parameterisation by elastic velocities allows using the larger base of singular vectors with the same condition number. In other words, it provides more stable results for the same level of input noise.

Linearised inverse problem for parametrisation $\left(V_p,V_s\right)$ is better conditioned (see Figure 3) in comparison with the case of elastic impedances.

Next step is resolution analysis. Following [23,24] for a fixed condition number of truncated operators, we construct the projection on stable subspaces for both operators. We built projections for a fixed condition number of truncated operator equal to ${10}^1$ (see Figure 4, Figure 5 and Figure 6). As one may recognise, there is no significant difference in projections for both parametrisations, which means they both are adequate for depth reflector reconstruction. We find it more natural to choose $\left(V_p,V_s\right)$ as elastic model parametrization.

3. Results

3.1. Experiment 1

The realistic example invoked in the study is based on the part of SEG/EAGE overthrust synthetic $V_p$-velocity model [27], $V_s$ is a scaled version of $V_p$ model (by a factor $1/\sqrt{3}$). Input data are synthesised for the set of ten uniformly sampled frequencies in the range 5 ÷ 10 Hz. The acquisition system has 31 vertical point force sources and 400 2C-receivers located at surface z = 0 m with a lateral spacing of 325 m and 25 m. Initial guess for the propagator is a vertically heterogeneous model (see Figure 7). As an initial guess for the data-space reflectivity variable s we used the observed data itself.

The results of elastic Full Waveform Inversion in the modified formulation one can see in Figure 7. To search the propagator, we use the space of 2D B-Splines functions of order 3. In this way, one may guarantee that the search is doing in the smooth propagator space. We perform the minimisation using the projected conjugate gradient method, where orthogonal projector onto the subspace spanned by B-splines is used as a projection onto the feasible set.

3.2. Experiment 2

The second experiment uses the synthetic Gullfaks model [28]. The original model is shown in Figure 8 (on the left). One should notice that the Vp/Vs ratio varies significantly throughout the model. The observed synthetic data are using a Ricker wavelet with a dominant frequency 10Hz. The source type is a vertical force, and we measure two velocity components at each receiver’s position. We use a sparse acquisition system placed at z = 0 m with 100 m spacing between seismic sources and 20 m spacing between receivers for the inversion. We performed the FWI in the time-frequency domain for the frequency range {5:1:20} Hz, where each time-frequency has an equal contribution in the inversion result. The starting velocities (Figure 8, in the middle) is a very poor vertically inhomogeneous approximation of the Gullfax real model. As a result of minimisation, we arrive at the model which is in good agreement with the original one for all three parameters: Vp, Vs, and Vp/Vs ratio (Figure 8, on the right). To demonstrate the data fit, we generated the synthetic seismogram for one source placed at x = 700 m (Figure 9). As one can see, the inverted velocities explain most of the data for the two receiver components.

3.3. Experiment 3

The last numerical example is the complex Marmousi2 realistic synthetic model [29]. The synthetic seismogram consists of 46 shot gathers with a shot interval 200 m and each shot is recorded by 460 receivers positioned at z = 0 having a 20 m receiver interval. We performed DSR FWI in the time-frequency domain for the frequency range {5:1:15} Hz. As in the previous example, each frequency has an equal contribution to the inversion process. We obtain the starting smooth model by horizontal averaging the original Marmousi2 velocities with the consequent gaussian smoothing. We used input data as an approximation of the DSR unknown s.

We show the inversion results in Figure 10. As before, we seek the propagator unknown in the smooth linear subspace spanned onto 2D B-Splines of order 3. In this way, we assure the evenness of the inverted propagator and significantly reduce the basis for the smooth models’ description. As one may observe, the inverted P- and S- velocities demonstrate a good correlation with the original Marmousi2 elastic model. To assess the reconstructed model’s accuracy, we compare the offset-depth domain Common Image Gathers (CIGs) in Figure 11. We can see that CIGs associated with the inverted velocity model are considerably flatter than those from the initial model.

4. Conclusions

We propose, investigate and implement the numerical method for solving the seismic inverse dynamical problem in the elastic isotropic approximation based on a modified least-squares technique. This modification bases on the decomposition of the wave propagation velocities into two subspaces: smoothly varying propagators and sharply variable reflectors. The modified functional has a strong nonlinearity to the macro velocity component. However, the smoothness of macrovelocities reduces the basis for their representation and opens the way to apply effective minimisation techniques. The presented numerical experiments confirm less sensitiveness of the modified version of FWI to the lack of the low time frequencies in the acquired data.