Elastic least-squares reverse time migration with velocities and density perturbation

Elastic least-squares reverse time migration with velocities and density perturbation Summary Elastic least-squares reverse time migration (LSRTM) based on the non-density-perturbation assumption can generate false-migrated interfaces caused by density variations. We perform an elastic LSRTM scheme with density variations for multicomponent seismic data to produce high-quality images in Vp, Vs and ρ components. However, the migrated images may suffer from crosstalk artefacts caused by P- and S-waves coupling in elastic LSRTM no matter what model parametrizations used. We have proposed an elastic LSRTM with density variations method based on wave modes separation to reduce these crosstalk artefacts by using P- and S-wave decoupled elastic velocity–stress equations to derive demigration equations and gradient formulae with respect to Vp, Vs and ρ. Numerical experiments with synthetic data demonstrate the capability and superiority of the proposed method. The imaging results suggest that our method promises imaging results with higher quality and has a faster residual convergence rate. Sensitivity analysis of migration velocity, migration density and stochastic noise verifies the robustness of the proposed method for field data. Image processing, Inverse theory, Seismic tomography, Wave propagation 1 INTRODUCTION With the development of multicomponent seismic exploration, conventional seismic imaging methods based on the acoustic approximation appear to be imperfect, because the assumption that only P waves propagate in the earth is not valid theoretically. Therefore, elastic imaging technologies attract more and more attention of exploration geophysicists. Several types of elastic imaging methods are used for multicomponent seismic data, mainly including ray-based migration methods, such as elastic Kirchhoff migration (Kuo & Dai 1984; Dai & Kuo 1986; Hokstad 2000; Du & Hou 2008) and Gaussian beam migration (Protasov & Tcheverda 2012; Han et al.2013), one-way wave equation migration (Wu & Xie 1994; Wu 1994; Xie & Wu 2005) and reverse time migration (RTM; Chang & McMechan 1987, 1994; Yan & Sava 2008; Qu et al.2015a,b; Yu et al.2016). In these methods, elastic RTM (ERTM) is the most promising and popular migration method because it is free from dip-angle limitation and high-frequency assumption. The coupling between P and S waves can generate crosstalk artefacts, resulting in inaccurate imaging results. To reduce crosstalk artefacts, ERTM algorithms use wave modes separation to produce scale images. Scalar RTM algorithm is first used for multicomponent data (Sun et al.2001, 2006). This method separates the multicomponent data into P- and S-component ones by calculating their divergence and curl, and uses them to implement acoustic RTM. The scalar RTM, ignoring the elastic characteristics of the vector wavefields, cannot accurately migrate converted waves. Therefore, vector ERTM methods, separating P and S waves when calculating forward and backward wavefields, have been developed over the past decade to produce PP, PS-, SP- and SS-wave images (Yan & Sava 2008). To obtain higher quality PS and SP images, polarity corrections need to be applied to S waves (Balch & Erdemir 1994; Du et al.2012), which cannot be accurately implemented in complex structures (Du et al.2017). In isotropic media, wave modes separation can be easily implemented by employing Helmholtz decomposition (Aki & Richards 2002; Qu et al.2015c) to the vector wavefields. But the calculation of divergence and curl operator will lead to a phase shift (Sun et al.2001) and amplitude changes (Sun et al. 2011). Recently, ERTM methods based on decoupled wave equations (Ma & Zhu 2003; Li et al.2007) are proposed (Gu et al.2015; Du et al.2017). Conventional migration applies the adjoint of a forward modeling operator to observed data instead of its inverse (Claerbout 1992; Nemeth et al.1999). A key problem with conventional migration is that it can produce a low-quality image with data acquisition footprint because of aliased seismic data, poor acquisition geometry, irregular data sampling, and a narrow source–receiver aperture (Kuehl & Sacchi 1999; Duquet et al.2000; Etgen et al.2009). Compared with the conventional RTM method, least-squares reverse time migration (LSRTM) is a powerful tool to produce a high-quality image with higher spatial resolution, fewer migration artefacts and better balanced reflector amplitudes (Dai & Schuster 2013; Tan & Huang 2014; Qu et al.2015c, 2016a; Hou & Symes 2015, 2016; Xue et al.2014, 2016; Yang et al.2016). With the development of computer technologies, elastic LSRTM (ELSRTM) is being applied to multicomponent seismic data (Duan et al.2016; Xu et al.2016; Feng & Schuster 2017). Ren et al. (2017) analyse the influence of different model parametrizations and objective functions in ELSRTM. Compared with acoustic LSRTM, ELSRTM produces better resolution images. These ELSRTM methods are all based on the non-density-perturbation assumption, resulting in the false-migrated interfaces caused by density variations. To overcome this problem, an ELSRTM scheme with density variations is developed. The approach is called V-ELSRTM. However, the crosstalk artefacts caused by the P and S waves coupling in multiparameter inversion are still a problem no matter what model parametrizations used (Feng & Schuster 2017). Métivier et al. (2014) incorporated the inverse Hessian operator within inversion method to alleviate this difficulty. We propose a new V-ELSRTM method based on P- and S-wavefield separation (V-SELSRTM) to reduce the crosstalk artefacts and enhance imaging quality. 2 ELSRTM FOR VELOCITIES AND DENSITY PERTURBATIONS BASED ON DECOUPLED ELASTIC WAVE EQUATION In the time domain, the velocity–stress wave equation with variable density is written as:   \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_x}}}{{\partial t}} = \frac{{\partial {\tau _{xx}}}}{{\partial x}} + \frac{{\partial {\tau _{xz}}}}{{\partial z}}\\ \rho \frac{{\partial {v_z}}}{{\partial t}} = \frac{{\partial {\tau _{xz}}}}{{\partial x}} + \frac{{\partial {\tau _{zz}}}}{{\partial z}}\\ \frac{{\partial {\tau _{xx}}}}{{\partial t}} = (\lambda + 2\mu )\frac{{\partial {v_x}}}{{\partial x}} + \lambda \frac{{\partial {v_z}}}{{\partial z}} + {f_x}\\ \frac{{\partial {\tau _{zz}}}}{{\partial t}} = (\lambda + 2\mu )\frac{{\partial {v_z}}}{{\partial z}} + \lambda \frac{{\partial {v_x}}}{{\partial x}} + {f_z}\\ \frac{{\partial {\tau _{xz}}}}{{\partial t}} = \mu \left(\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}\right) \end{array} \right.,\end{equation} (1)where vx and vz are the velocities in horizontal and vertical components, respectively, τxx and τzz are the normal stress, τxz is the shear stress, fx and fz are the two components of the source terms, t is time, ρ is the density, λ and μ are the Lamé parameters. The relationships between the Lamé parameters λ and μ, density ρ, P-wave velocity Vp, and S wave velocity Vs are   \begin{equation}\left\{ \begin{array}{@{}l@{}} \lambda = \rho (V_p^2 - 2V_s^2)\\ \mu = \rho V_s^2 \end{array} \right..\end{equation} (2) To reduce crosstalk artefacts of V-ELSRTM caused by the coupling of different wave modes, we develop a V-SELSRTM method for velocities and density perturbations based on the decoupled elastic wave equation. The decoupled elastic velocity–stress wave equations can separate coupled elastic waves into P- and S-mode waves, and are given by (Li et al.2007)   \begin{equation}\left\{ \begin{array}{@{}l@{}} {v_x} = {v_{xp}} + {v_{xs}}\\ {v_z} = {v_{zp}} + {v_{zs}} \end{array} \right.,\end{equation} (3a)  \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_{xp}}}}{{\partial t}} = \frac{{\partial {\tau _{xxp}}}}{{\partial x}}\\ \rho \frac{{\partial {v_{zp}}}}{{\partial t}} = \frac{{\partial {\tau _{zzp}}}}{{\partial z}}\\ \frac{{\partial {\tau _{xxp}}}}{{\partial t}} = (\lambda + 2\mu )\left(\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}\right) + {f_x}\\ \frac{{\partial {\tau _{zzp}}}}{{\partial t}} = (\lambda + 2\mu )\left(\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}\right) + {f_z} \end{array} \right.,\end{equation} (3b)  \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_{xs}}}}{{\partial t}} = \frac{{\partial {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial {\tau _{xzs}}}}{{\partial z}}\\ \rho \frac{{\partial {v_{zs}}}}{{\partial t}} = \frac{{\partial {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial {\tau _{zzs}}}}{{\partial z}}\\ \frac{{\partial {\tau _{xxs}}}}{{\partial t}} = - 2\mu \frac{{\partial {v_z}}}{{\partial z}} + {f_x}\\ \frac{{\partial {\tau _{zzs}}}}{{\partial t}} = - 2\mu \frac{{\partial {v_x}}}{{\partial x}} + {f_z}\\ \frac{{\partial {\tau _{xzs}}}}{{\partial t}} = \mu \left(\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}\right) \end{array} \right.,\end{equation} (3c)where the subscript p and s denote the P-wavefield Up and S-wavefield Us, respectively. Eqs (3a)–(3c) are mixed, P-wave and S-wave equations, respectively. Source terms fx and fz are added in both P- and S-wave equations. Only when shear sources are excited, source term is added only in the fifth equation of S-wave equations (3c). The correctness of the P- and S-wave modes separated first-order velocity–stress wave equation is demonstrated in Appendix A. Least-squares migration can be seen as a method to solve an optimization problem. The purpose is to find a perturbation model in the least-squares sense to fit the observed data by minimizing the following objective function:   \begin{equation}E({{\bf m}}) = \frac{1}{2}\left\| {{{{\bf d}}_{{\rm{cal}}}} - {{{\bf d}}_{{\rm{obs}}}}} \right\|_2^2,\end{equation} (4)where dcal denote the synthetic data and dcal = (dcalp, dcals), and dobs denote the observed data and dobs = (dobsp, dobss); the synthetic data can be given by dcal = Lm, in which L denotes a linear modeling operator and L = (Lp, Ls), m is the reflectivity model. Perturbation wavefields δUp and δUs can be obtained by solving (details see Appendix C)   \begin{equation}\left\{ \begin{array}{@{}l@{}} \delta {v_x} = \delta {v_{xp}} + \delta {v_{xs}}\\ \delta {v_z} = \delta {v_{zp}} + \delta {v_{zs}} \end{array} \right.,\end{equation} (5a)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxp}}}}{{\partial x}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{xp0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{zzp}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{zp0}}}}{{\partial t}}\\ \frac{1}{{V_{p0}^2}}\frac{{\partial \delta {\tau _{xxp}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + \frac{{2{{\bf m}}({V_p})}}{{V_{p0}^2}}\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\\ \frac{1}{{V_{p0}^2}}\frac{{\partial \delta {\tau _{zzp}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + \frac{{2{{\bf m}}({V_p})}}{{V_{p0}^2}}\frac{{\partial {\tau _{zzp0}}}}{{\partial t}} \end{array} \right.,\end{equation} (5b)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{xzs}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{xs0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{zzs}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{zs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{xxs}}}}{{\partial t}} = - 2{\rho _0}\frac{{\partial \delta {v_z}}}{{\partial z}} + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{zzs}}}}{{\partial t}} = - 2{\rho _0}\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{zzs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{xzs}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial z}} + \frac{{\partial \delta {v_z}}}{{\partial x}}\right) + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{xzs0}}}}{{\partial t}} \end{array} \right.,\end{equation} (5c)where Up0 and Us0 are the solutions associated with the background parameters V0 of eq. (3), δUp and δUs denote the perturbation wavefields, in which Up = (vxp, vzp, τxxp, τzzp) and Us= (vxs, vzs, τxxs, τzzs, τxzs), and δV is parameter perturbations and V = (Vp, Vs, ρ). The reflectivity models associated with Vp, Vs and ρ are defined by   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{\bf m}}({V_p}) = \frac{{\delta {V_p}}}{{{V_{p0}}}}\\ {{\bf m}}({V_s}) = \frac{{\delta {V_s}}}{{{V_{s0}}}}\\ {{\bf m}}(\rho ) = \frac{{\delta \rho }}{{{\rho _0}}} \end{array} \right..\end{equation} (6) The P- and S-wave migration equations can also be given by   \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial v_{xp}^R}}{{\partial t}} = \frac{{\partial \tau _{xxp}^R}}{{\partial x}} + \Delta {{{\bf d}}_x}({{{\bf x}}_r},t)\\ \rho \frac{{\partial v_{zp}^R}}{{\partial t}} = \frac{{\partial \tau _{zzp}^R}}{{\partial z}} + \Delta {{{\bf d}}_z}({{{\bf x}}_r},t)\\ \frac{{\partial \tau _{xxp}^R}}{{\partial t}} = \rho V_p^2\left(\frac{{\partial v_x^R}}{{\partial x}} + \frac{{\partial v_z^R}}{{\partial z}}\right)\\ \frac{{\partial \tau _{zzp}^R}}{{\partial t}} = \rho V_p^2\left(\frac{{\partial v_x^R}}{{\partial x}} + \frac{{\partial v_z^R}}{{\partial z}}\right) \end{array} \right.,\end{equation} (7a)  \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial v_{xs}^R}}{{\partial t}} = \frac{{\partial \tau _{xxs}^R}}{{\partial x}} + \frac{{\partial \tau _{xzs}^R}}{{\partial z}} + \Delta {{{\bf d}}_x}({{{\bf x}}_r},t)\\ \rho \frac{{\partial v_{zs}^R}}{{\partial t}} = \frac{{\partial \tau _{xzs}^R}}{{\partial x}} + \frac{{\partial \tau _{zzs}^R}}{{\partial z}} + \Delta {{{\bf d}}_z}({{{\bf x}}_r},t)\\ \frac{{\partial \tau _{xxs}^R}}{{\partial t}} = - 2\rho V_s^2\frac{{\partial v_z^R}}{{\partial z}}\\ \frac{{\partial \tau _{zzs}^R}}{{\partial t}} = - 2\rho V_s^2\frac{{\partial v_x^R}}{{\partial x}}\\ \frac{{\partial \tau _{xzs}^R}}{{\partial t}} = \rho V_s^2\left(\frac{{\partial v_x^R}}{{\partial z}} + \frac{{\partial v_z^R}}{{\partial x}}\right) \end{array} \right.,\end{equation} (7b)where $${{{\bf U}}^R} = ({{\bf U}}_p^R,{{\bf U}}_s^R)$$ denotes the adjoint matrix of U = (Up, Us), xr denotes the location of receivers, Δdx(xr, t) and Δdz(xr, t) are the difference between observed data and synthetic data in horizontal and vertical directions, respectively, which are calculated by   \begin{equation}\left\{ \begin{array}{@{}l@{}} \Delta {{{\bf d}}_x}({{{\bf x}}_r},t) = {{{\bf d}}_{{\rm{cal}}x}}({{{\bf x}}_r},t) - {{{\bf d}}_{{\rm{obs}}x}}({{{\bf x}}_r},t)\\ \Delta {{{\bf d}}_z}({{{\bf x}}_r},t) = {{{\bf d}}_{{\rm{cal}}z}}({{{\bf x}}_r},t) - {{{\bf d}}_{{\rm{obs}}z}}({{{\bf x}}_r},t) \end{array} \right.,\end{equation} (8)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {{{\bf d}}_{{\rm{cal}}x}}({{{\bf x}}_r},t) = \delta {v_x}({{{\bf x}}_r},t)\\ {{{\bf d}}_{{\rm{cal}}z}}({{{\bf x}}_r},t) = \delta {v_z}({{{\bf x}}_r},t) \end{array} \right..\end{equation} (9) Noted that the input observed data are multicomponent shot records dobs (dobsx, dobsz) without separating P and S wavefields. The predicted data are mixed multicomponent shot records dcal (dcalx, dcalz), which are mixed by eq. (3a). And then the residual P and S wavefields are separated during the backpropagation. Then the gradient formulae with respect to Vp, Vs and ρ can be calculated by (details see Appendix C)   \begin{eqnarray}\left\{ \begin{array}{@{}l@{}} {{{\bf g}}_{vp}} = \frac{{\partial {{\bf E}}}}{{\partial {V_p}}} = \frac{1}{{V_{p0}^3}}\int{{\left( {\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\tau _{xxp}^R + \frac{{\partial {\tau _{zzp0}}}}{{\partial t}}\tau _{zzp}^R} \right)}}dt\\ {{{\bf g}}_{vs}} = \frac{{\partial {{\bf E}}}}{{\partial {V_s}}} = \frac{1}{{V_{s0}^3}}\int{{\left( {\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}\tau _{xxs}^R + \frac{{\partial {\tau _{zzs0}}}}{{\partial t}}\tau _{zzs}^R + \frac{{\partial {\tau _{xzs0}}}}{{\partial t}}\tau _{xzs}^R} \right)}}dt\\ {{{\bf g}}_{\rho v}} = \frac{{\partial {{\bf E}}}}{{\partial \rho }} = \int{{\left( {\frac{{\partial {v_{xp0}}}}{{\partial t}}v_{xp}^R + \frac{{\partial {v_{zp0}}}}{{\partial t}}v_{zp}^R + \frac{{\partial {v_{xs0}}}}{{\partial t}}v_{xs}^R + \frac{{\partial {v_{zs0}}}}{{\partial t}}v_{zs}^R} \right)}}dt \end{array} \right.,\nonumber\\ \end{eqnarray} (10)where gvp, gvs and gρv are gradient formulae with respect to Vp, Vs and ρ, respectively. With the relationships between velocities (Vp and Vs) and impedances (P impedance Ip and S impedance Is)   \begin{equation}\left\{ \begin{array}{@{}l@{}} {I_p} = \rho {V_p}\\ {I_s} = \rho {V_s} \end{array} \right.,\end{equation} (11)the reflectivity of Ip and Iscan be written by   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{\bf m}}({I_p}) = \frac{{\delta {I_p}}}{{{I_{p0}}}} = \frac{{\delta (\rho {V_p})}}{{{\rho _0}{V_{p0}}}} = \frac{{\delta {V_p}}}{{{v_{p0}}}} + \frac{{\delta \rho }}{{{\rho _0}}} = {{\bf m}}({V_p}) + {{\bf m}}(\rho )\\ {{\bf m}}({I_s}) = \frac{{\delta {I_s}}}{{{I_{s0}}}} = \frac{{\delta (\rho {V_s})}}{{{\rho _0}{V_{s0}}}} = \frac{{\delta {V_s}}}{{{v_{s0}}}} + \frac{{\delta \rho }}{{{\rho _0}}} = {{\bf m}}({V_s}) + {{\bf m}}(\rho ) \end{array} \right..\end{equation} (12) The gradient formulae with respect to Ip, Is and ρ can be obtained by applying the chain rule (Mora 1987; Crase et al.1990; Köhn et al.2012):   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{{\bf g}}_{Ip}} = \frac{1}{\rho }{{{\bf g}}_{vp}}\\ {{{\bf g}}_{Is}} = \frac{1}{\rho }{{{\bf g}}_{vs}}\\ {{{\bf g}}_{\rho i}} = - \frac{{{V_p}}}{\rho }{{{\bf g}}_{vp}} - \frac{{{V_s}}}{\rho }{{{\bf g}}_{vs}} + {{{\bf g}}_{\rho v}} \end{array} \right..\end{equation} (13) Similarly, when performing V-ELSRTM for Lamé parameters and density perturbations based on decoupled elastic wave equation, the gradients related to λ and μ, ρ can be accessed as (Appendix C)   \begin{eqnarray}\left\{ \begin{array}{@{}l@{}} {{{\bf g}}_\lambda } = \int{{\left( {\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)(\tau _{xxp}^R + \tau _{zzp}^R)} \right)}}dt\\ {{{\bf g}}_\mu } = \int\left( 2\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)(\tau _{xxp}^R + \tau _{zzp}^R) - \frac{{\partial {v_{z0}}}}{{\partial z}}\tau _{xxs}^R - \frac{{\partial {v_{x0}}}}{{\partial x}}\tau _{zzs}^R \right. \\ \qquad \qquad\left.+\left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right)\tau _{xzs}^R \right)dt\\ {{{\bf g}}_\rho } = \int{{\left( {\frac{{\partial {v_{xp0}}}}{{\partial t}}v_{xp}^R + \frac{{\partial {v_{zp0}}}}{{\partial t}}v_{zp}^R + \frac{{\partial {v_{xs0}}}}{{\partial t}}v_{xs}^R + \frac{{\partial {v_{zs0}}}}{{\partial t}}v_{zs}^R} \right)}}dt \end{array} \right.\end{eqnarray} (14) The gradients for velocities and impedances can also be calculated by applying the chain rule. In this paper, we use eq. (10) to update V-SELSRTM images. The pre-conditioned conjugate gradient (Dai et al.2012) method is applied to solve the least-squares inverse problem. The calculation flow of elastic LSRTM with velocities and density perturbation is illustrated in Fig. 1. The left- and right-hand columns show the conventional V-ELSRTM (details see Appendix B) and the V-SELSRTM, respectively. Figure 1. View largeDownload slide The calculation flow of the conventional V-ELSRTM (the left-hand column) and V-SELSRTM (the right-hand column). Figure 1. View largeDownload slide The calculation flow of the conventional V-ELSRTM (the left-hand column) and V-SELSRTM (the right-hand column). 3 NUMERICAL EXAMPLES First, the performances of the P- and S-wave modes separated elastic wave equation and linear modeling operators are demonstrated through a numerical example using a simple two-layer model. Then the application of the proposed V-SELSRTM compared with the conventional V-ELSRTM and ELSRTM based on the non-density-perturbation assumption (N-ELSRTM, Feng & Schuster 2017) is demonstrated with synthetic data examples using two models: (1) a layered elastic model with different Vp and Vs anomalies and (2) the SEG/EAGE salt model. In these synthetic examples, the observed data are produced by linear forward modeling, which has variable density. All our simulations are implemented using a time-domain staggered-grid finite-difference algorithm with the second-order accurate in time and tenth-order-accurate in space scheme. Multisource LSRTM is introduced to compress the seismic data and improve the efficiency of LSRTM (Dai et al.2012). Dynamic encoding technology is used to suppress the crosstalks caused by different sources in one supershot (Qu et al.2016b, 2017), which are different from those caused by coupling of P and S waves. To suppress artificial boundary reflections, complex frequency-shifted perfectly matched layer (CFS-PML) boundary conditions (Drossaert & Giannopoulos 2007) are used in all four edges of the models. 3.1 Wavefield separation test with a simple two-layer model We consider an isotropic elastic two-layer model characterized by the P velocity, S velocity and density illustrated in Figs 2(a)–(c) to demonstrate the effectiveness of the P- and S-wave modes separated elastic wave equation and linear modeling operators. The flat interfaces of the P velocity, S velocity and density models are located at different depths. The physical property parameters are presented in Fig. 2. We choose the two-layer model because the P and S waves have different traveltime, and therefore we can determine their individual amplitudes and phases without mutual interference. The dimensions of the model are 201 × 201 gridpoints in the horizontal and vertical directions, with 8 m grid spacing. An explosive source using the Ricker wavelet with a dominant frequency of 25 Hz is excited at (800 m, 8 m), and the total recording length is 1.5 s, with a time sampling of 0.5 ms. The 201 geophones are spread out over the surface. Fig. 3 shows the snapshots at 240 and 320 ms generated using the P- and S-wave modes separated elastic wave equations (3). A comparison of the mixed- and separated-mode wavefield indicates that the wave modes separation method works well for isotropic media. The reflectivity models of the two-layer model (calculated by eq. 6) are displayed in Fig. 4, which are used to perform linear modeling. Fig. 5 shows mixed synthetic shot records obtained from eq. (B4) (Figs 5a and d) and separated synthetic shot records obtained from eq. (5) (Figs 5b,c,e and f). These results suggest that the primary reflected P and S waves are simulated accurately and separated completely and there are no change in amplitude and phase for both P and S waves after wavefield separation, demonstrating the correctness of eq. (5) and the effectiveness of the linear modeling for V-SELSERM. Figure 2. View largeDownload slide A two-layer model. (a) Vp, (b) Vs and (c) ρ. Figure 2. View largeDownload slide A two-layer model. (a) Vp, (b) Vs and (c) ρ. Figure 3. View largeDownload slide Snapshots at (a)–(f) 240 ms and (g)–(l) 320 ms. (a) and (g) Mixed wavefield in x-direction; (b) and (h) P wavefield in x-direction; (c) and (i) S wavefield record in x-direction; (d) and (j) mixed wavefield in z-direction; (e) and (k) P wavefield in z-direction; (f) and (l) S wavefield in z-direction. Figure 3. View largeDownload slide Snapshots at (a)–(f) 240 ms and (g)–(l) 320 ms. (a) and (g) Mixed wavefield in x-direction; (b) and (h) P wavefield in x-direction; (c) and (i) S wavefield record in x-direction; (d) and (j) mixed wavefield in z-direction; (e) and (k) P wavefield in z-direction; (f) and (l) S wavefield in z-direction. Figure 4. View largeDownload slide The reflectivity of the two-layer model. (a) Vp, (b) Vs and (c) ρ. Figure 4. View largeDownload slide The reflectivity of the two-layer model. (a) Vp, (b) Vs and (c) ρ. Figure 5. View largeDownload slide (a) and (d) Mixed synthetic shot records obtained from eq. (B4) and (b), (c), (e) and (f) separated synthetic shot records obtained from eq. (5). (a)–(c) x-direction and (d)–(f) z-direction. Figure 5. View largeDownload slide (a) and (d) Mixed synthetic shot records obtained from eq. (B4) and (b), (c), (e) and (f) separated synthetic shot records obtained from eq. (5). (a)–(c) x-direction and (d)–(f) z-direction. 3.2 A layered elastic model with different Vp and Vs anomalies The second model considered here is a layered elastic model with different Vp and Vs anomalies (Fig. 6). We use this model to demonstrate the advantages of V-ELSRTM and V-SELSRTM over N-ELSRTM. The velocity models have three flat layers with a distance of 4 km × 2 km in x- and z-directions, discretized on a 401 × 201 grids with a spatial grid interval of 10 m. Two anomalies with Vp and Vs of 4000 and 2309 m s−1 are inserted into different locations of the second layer, the size of which is 1.0 km × 0.4 km in x- and z-directions. The material properties of the elastic model are indicated in Fig. 6. The migration Vp, Vs and density models (Figs 6a, c and e) are smoothed versions of the true ones (Figs 6b,d and f). The wide-aperture survey of generating the synthetic data is made up of 200 explosive P-wave sources spaced every 20 m and located in a depth of 10 m, which emit Ricker wavelets with the dominant frequency of 25 Hz. These shots are compressed into 10 supershots by using dynamic encodes. The synthetic records are received by 401 two-component receivers, spaced every 10 m on the surface, calculated up to 1.6 s with a time sample rate of 0.8 ms. Fig. 7 presents the 101st shot gather with removed direct waves recorded at surface. Figure 6. View largeDownload slide An elastic layered model. (a) True Vp, (b) migration Vp, (c) true Vs, (d) migration Vs, (e) true ρ and (e) migration ρ. Figure 6. View largeDownload slide An elastic layered model. (a) True Vp, (b) migration Vp, (c) true Vs, (d) migration Vs, (e) true ρ and (e) migration ρ. Figure 7. View largeDownload slide A shot gather (shot 101). (a) x-direction and (b) z-direction. Figure 7. View largeDownload slide A shot gather (shot 101). (a) x-direction and (b) z-direction. We first perform the conventional N-ELSRTM on the synthetic data sets generated with variable density. Fig. 8 gives the N-ELSRTM images of the 40th iteration after Laplacian filtering for Vp and Vs components. The ELSRTM image in Vs component has higher resolution than that in Vp component because S wave has shorter wavelength than P wave. However, in the N-ELSRTM images, there are false interfaces caused by the variable density, which are indicated by the red ellipses. Then, we conduct the three V-ELSRTM methods using the same synthetic shot records on the layered model. Figs 9 and 10 display the forward- and backward-propagated wavefields at 400 ms, in which the mixed and separated wavefields are generated by the solutions of eqs (3) and (7), respectively. The results further verify that the P and S waves separated forward modeling method can simulate pure P and S wavefields during seismic wave forward and backward propagation, thus providing a theoretical basis for the proposed V-SELSRTM method. Figure 8. View largeDownload slide Comparisons of (a) the true Vp perturbation, (b) the true Vs perturbation, (c) conventional N-ELSRTM images after 40 iterations for Vp component and (d) conventional N-ELSRTM for Vs component. Figure 8. View largeDownload slide Comparisons of (a) the true Vp perturbation, (b) the true Vs perturbation, (c) conventional N-ELSRTM images after 40 iterations for Vp component and (d) conventional N-ELSRTM for Vs component. Figure 9. View largeDownload slide Forward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 9. View largeDownload slide Forward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 10. View largeDownload slide Backward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 10. View largeDownload slide Backward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Fig. 11 compares the V-ELSRTM and V-SELSRTM images of the 40th iteration after Laplacian filtering. For Vp-, Vs- and ρ-component images, the amplitude is balanced with high resolution and signal-to-noise ratio (SNR). The false interfaces caused by the perturbation of velocities (indicated by the blue ellipses) or density (indicated by the red ellipses) gradually weakened and disappeared against the iteration number, which is demonstrated by the fact that the V-ELSRTM and V-SELSRTM images for Vp, Vs and ρ components have the same pattern as the true Vp (Fig. 11a), Vs (Fig. 11b) and ρ perturbations (Fig. 11c). By the comparisons of the V-ELSRTM and V-SELSRTM images, there are some artefacts in V-ELSRTM images caused by P- and S-wave crosstalks (shown by the red arrows in Figs 11d and e). These crosstalk artefacts are mitigated when using the V-SELSRTM (Figs 11g and h). Fig. 12 presents the convergence curves of model (Fig. 12a) and data residuals (Fig. 12b) over iterations for N-ELSRTM, V-ELSRTM and V-SELSRTM. We track the data residuals using the L2 norm of mixed data residuals as the conventional ELSRTM. It can be observed that in the V-ELSRTM and V-SELSRTM cases, these residuals converge to much smaller values than the N-ELSRTM case. The fluctuations of the residuals are as a result of the dynamic encodes change. In particular, after 10 iterations, the model and data residuals with the N-ELSRTM are clearly larger than those of V-ELSRTM and V-SELSRTM, and the convergence curve of the N-ELSRTM becomes flat after 25 iterations. The curves of V-SELSRTM have faster convergence speeds, and the model and data residuals converge to smaller minimums than V-ELSRTM. Fig. 13 compares the corresponding V-ELSRTM and V-SELSRTM images using seismic impedances (Ip, Is and ρ) with the true impedance perturbations (m(Ip), m(Is) and m(ρ)), in which the false interfaces caused by the perturbation of velocities are indicated by the blue ellipses. The results also suggest that the V-SELSRTM images of impedances show better accordance with the true perturbations of impedances than the V-ELSRTM by reducing P- and S-waves coupling crosstalks (shown by the red arrows). Figure 11. View largeDownload slide (a) Comparisons of the true Vp perturbation, (b) the true Vs perturbation , (c) the true ρ perturbation, (d) V-ELSRTM images for Vp component, (e) Vs component and (f) ρ component, (g) V-SELSRTM images for Vp component, (h) Vs component and ρ component. Figure 11. View largeDownload slide (a) Comparisons of the true Vp perturbation, (b) the true Vs perturbation , (c) the true ρ perturbation, (d) V-ELSRTM images for Vp component, (e) Vs component and (f) ρ component, (g) V-SELSRTM images for Vp component, (h) Vs component and ρ component. Figure 12. View largeDownload slide (a) Normalized model residuals and (b) data residuals. The black lines indicate the residuals obtained from the conventional N-ELSRTM, the red lines indicate the residuals obtained from V-ELSRTM and the blue lines indicate the residuals obtained from V-SELSRTM. Figure 12. View largeDownload slide (a) Normalized model residuals and (b) data residuals. The black lines indicate the residuals obtained from the conventional N-ELSRTM, the red lines indicate the residuals obtained from V-ELSRTM and the blue lines indicate the residuals obtained from V-SELSRTM. Figure 13. View largeDownload slide (a) Comparisons of the true Ip perturbation, (b) the true Is perturbation, (c) the true ρ perturbation, (d) V-ELSRTM images for Ip component, (e) Is component, (f) ρ component, (g) V-SELSRTM images for Ip component, (h) Is component and ρ component. Figure 13. View largeDownload slide (a) Comparisons of the true Ip perturbation, (b) the true Is perturbation, (c) the true ρ perturbation, (d) V-ELSRTM images for Ip component, (e) Is component, (f) ρ component, (g) V-SELSRTM images for Ip component, (h) Is component and ρ component. 3.3 SEG/EAGE salt model To further analyse the capacity of our V-SELSRTM for imaging complex model, we apply the new algorithm to a SEG/EAGE salt model (Fig. 14), discretized on an 645 × 150 grids with a grid spacing of Δx = Δz = 8 m. In the fixed-spread acquisition geometry, we evenly deploy 320 sources on the surface, with the interval between neighbouring sources of 16 m. These shots are compressed into 16 supershots by using dynamic encodes. The explosive point source chooses a Ricker wavelet with the dominant frequency of 20 Hz. There are 645 fixed two-component receivers with a spacing of 8 m. To ensure numerical stability, sample every 0.8 ms across a 3.5 s recording period. There are 50 CFS-PML grids surrounding the SEG/EAGE model boundaries. Figs 15 and 16 illustrate the snapshot comparisons of forward- and backward-propagated wavefields, respectively, after 0.4 s. Figure 14. View largeDownload slide The SEG/EAGE salt model. (a) True Vp, (b) migration Vp, (c) true Vs, (d) migration Vs, (e) true ρ and (f) migration ρ. Figure 14. View largeDownload slide The SEG/EAGE salt model. (a) True Vp, (b) migration Vp, (c) true Vs, (d) migration Vs, (e) true ρ and (f) migration ρ. Figure 15. View largeDownload slide Forward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 15. View largeDownload slide Forward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 16. View largeDownload slide Backpropagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 16. View largeDownload slide Backpropagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. The true reflectivity models are displayed in Fig. 17, which are used only for comparison. Fig. 18 plots the V-ELSRTM and V-SELSRTM images after 40 iterations. Similar to the previous example, the Vs-component images show higher resolution than the Vp-component images. The ρ-component images achieve lower SNR than Vp- and Vs-component images, which have some low-frequency noises. V-ELSRTM images show the same resolution as that of V-SELSRTM. However, the V-SELSRTM produces much fewer image artefacts than V-ELSRTM. To compare the resolution and SNR of the salt and subsalt areas in both images, we illustrate the magnified images in Figs 19 and 20. The red and blue boxes in Fig. 18 show the regions for zoom views. The magnified views suggest that the V-ELSRTM method produces more serious ringing artefacts around the reflecting interfaces and the V-SELSRTM method generates better images of the salt and subsalt structures with fewer ringing and low-frequency artefacts. Fig. 21 gives V-SELSRTM images after Laplacian filtering, showing high-quality images with high resolution, fewer artefacts and balanced amplitude. The vertical profiles of the true reflectivity and migration images at the horizontal distance of 3.36 km are shown in Fig. 22. It is obviously observed that the Vp-, Vs- and ρ-component images obtained from the V-SELSRTM method are more accurate than those obtained from the V-ELSRTM method. Fig. 23 presents the convergence curve of data and model residuals for V-ELSRTM and V-SELSRTM. By comparison, V-ELSRTM has a lower convergence rate and residual of V-ELSRTM converges to a larger value than that of the V-SELSRTM. Figure 17. View largeDownload slide The true reflectivity models. (a) The Vp perturbation, (b) the Vs perturbation and (c) the ρ perturbation. Figure 17. View largeDownload slide The true reflectivity models. (a) The Vp perturbation, (b) the Vs perturbation and (c) the ρ perturbation. Figure 18. View largeDownload slide (a) Comparisons of V-ELSRTM images after 40 iterations for Vp component, (c) Vs component, (e) ρ component, (b) V-SELSRTM images for Vp component, (d) Vs component and (f) ρ component. Figure 18. View largeDownload slide (a) Comparisons of V-ELSRTM images after 40 iterations for Vp component, (c) Vs component, (e) ρ component, (b) V-SELSRTM images for Vp component, (d) Vs component and (f) ρ component. Figure 19. View largeDownload slide Magnified views of Fig. 18 shown by the red boxes. (a) V-ELSRTM images for Vp component, (b) Vs component, (c) ρ component, (d) V-SELSRTM images for Vp component, (e) Vs component and (f) ρ component. Figure 19. View largeDownload slide Magnified views of Fig. 18 shown by the red boxes. (a) V-ELSRTM images for Vp component, (b) Vs component, (c) ρ component, (d) V-SELSRTM images for Vp component, (e) Vs component and (f) ρ component. Figure 20. View largeDownload slide Magnified views of Fig. 18 shown by the blue boxes. (a) V-ELSRTM images for Vp component, (b) Vs component, (c) ρ component, (d) V-SELSRTM images for Vp component, (e) Vs component and (f) ρ component. Figure 20. View largeDownload slide Magnified views of Fig. 18 shown by the blue boxes. (a) V-ELSRTM images for Vp component, (b) Vs component, (c) ρ component, (d) V-SELSRTM images for Vp component, (e) Vs component and (f) ρ component. Figure 21. View largeDownload slide V-SELSRTM images after Laplacian filtering for (a) Vp component, (b) Vs component and (c) ρ component. Figure 21. View largeDownload slide V-SELSRTM images after Laplacian filtering for (a) Vp component, (b) Vs component and (c) ρ component. Figure 22. View largeDownload slide Vertical profile at the horizontal distance of 3.36 km. (a) Vp component, (b) Vs component and (c) ρ component. The solid red lines indicate the true reflectivity model, the dashed black lines denote V-ELSRTM images after 40 iterations and the dashed blue lines denote V-SELSRTM images after 40 iterations. Figure 22. View largeDownload slide Vertical profile at the horizontal distance of 3.36 km. (a) Vp component, (b) Vs component and (c) ρ component. The solid red lines indicate the true reflectivity model, the dashed black lines denote V-ELSRTM images after 40 iterations and the dashed blue lines denote V-SELSRTM images after 40 iterations. Figure 23. View largeDownload slide (a) Normalized residuals, where the black line with circles denote the model residual for V-ELSRTM, the red line with circles denote the model residual for V-SELSRTM, the black line with stars denote the data residual for V-ELSRTM and the red line with stars denote the data residual for V-SELSRTM. (b) Normalized data residuals in log domain. Figure 23. View largeDownload slide (a) Normalized residuals, where the black line with circles denote the model residual for V-ELSRTM, the red line with circles denote the model residual for V-SELSRTM, the black line with stars denote the data residual for V-ELSRTM and the red line with stars denote the data residual for V-SELSRTM. (b) Normalized data residuals in log domain. 3.4 Sensitivity analysis to migration velocity error In real cases, migration velocity errors should be taken into account to produce high-quality images. To study the sensitivity of our V-SELSRTM to migration Vp and Vs errors, we smooth the migration Vp and Vs models using different smoothing parameter values and add small errors to them (shown in Fig. 24), with the assumption that the migration density is correct (Fig. 14f). For quantification, we give the relative velocity error by using the following equation   \begin{equation}{\rm{ERROR}} = \frac{{||{{{\bf m}}_{{\rm{true}}}} - {{{\bf m}}_{{\rm{init}}}}||}}{{||{{{\bf m}}_{{\rm{true}}}}||}}.\end{equation} (15) Figure 24. View largeDownload slide Migration velocities with different error of (a) and (b) 7.1 per cent, (c) and (d) 10.8 per cent, (e) and (f) 13.2 per cent and (g) and (h) 16.5 per cent. (a), (c), (e) and (g) migration Vp. (b), (d), (f) and (h) migration Vs. Figure 24. View largeDownload slide Migration velocities with different error of (a) and (b) 7.1 per cent, (c) and (d) 10.8 per cent, (e) and (f) 13.2 per cent and (g) and (h) 16.5 per cent. (a), (c), (e) and (g) migration Vp. (b), (d), (f) and (h) migration Vs. The relative velocity errors between the migration Vp and Vs models and the original ones are 7.1 per cent (Figs 24a, b), 10.8 per cent (Figs 24c and d), 13.2 per cent (Figs 24e and f) and 16.5 per cent (Figs 24g and h). The V-SELSRTM images obtained with the different wrong migration Vp and Vs models are displayed in Fig. 25. Because of the errors in Vp and Vs models, there are some spurious images and obvious low-frequency noises in the shallow and the salt areas, and the subsalt zone is weakly imaged. The wrong velocities would image the reflection data towrong positions. The ρ-component images are more seriously influenced by the velocity errors. As the errors in Vp and Vs models increase, the quality of the V-SELSRTM images gradually decreases. The convergence curves of the data residual over iterations when using imperfect migration velocities are illustrated in Fig. 26. The convergence curves for the four cases convergence to different values with different convergence rate, which suggest that our V-SELSRTM method can produce acceptable inversion results as long as the migration velocity models are within pre-defined thresholds. Figure 25. View largeDownload slide (a)–(c) The V-SELSRTM images after 40 iterations using different migration velocities with the error of 7.1 per cent, (d)–(f)10.8 per cent, (g)–(i) 13.2 per cent and (j)–(l) 16.5 per cent. (a), (d), (g) and (j) Vp component, (b), (e), (h) and (k) Vs component, (c), (f), (i) and (l) ρ component. Figure 25. View largeDownload slide (a)–(c) The V-SELSRTM images after 40 iterations using different migration velocities with the error of 7.1 per cent, (d)–(f)10.8 per cent, (g)–(i) 13.2 per cent and (j)–(l) 16.5 per cent. (a), (d), (g) and (j) Vp component, (b), (e), (h) and (k) Vs component, (c), (f), (i) and (l) ρ component. Figure 26. View largeDownload slide Normalized data residual using migration velocities with different error. The black, red, blue and pink lines denote migration velocities with the error of 7.1 per cent, 10.8 per cent, 13.2 per cent and 16.5 per cent, respectively. Figure 26. View largeDownload slide Normalized data residual using migration velocities with different error. The black, red, blue and pink lines denote migration velocities with the error of 7.1 per cent, 10.8 per cent, 13.2 per cent and 16.5 per cent, respectively. 3.5 Migration density error sensitivity analysis We then study the sensitivity of our V-SELSRTM to the density error in the migration density model. In real cases, the density model is difficult to obtain accurately. Therefore in this case, we use a constant migration density model (the value is set to be 2411 kg m−3 by averaging the true density values) and the accurate migration velocity models (Figs 14b and d) to conduct the V-SELSRTM. Fig. 27 shows the V-SELSRTM images after 40 iterations using a constant migration densitymodel. The ρ-component image shows serious low-frequency noise and weak-imaging energy beneath subsalt areas, in turn resulting in some artefacts on migrated VP- and VS-component images. Fig. 28 also displays the convergence curve of data residuals. The data residuals of V-SELSRTM decrease over iterations and converge to 10−2. In general, the V-SELSRTM is less sensitive to the errors in the migration density model than the migration velocity models. Figure 27. View largeDownload slide The V-SELSRTM images after 40 iterations using constant migration density model. (a) Vp component, (b) Vs component and (c) ρ component. Figure 27. View largeDownload slide The V-SELSRTM images after 40 iterations using constant migration density model. (a) Vp component, (b) Vs component and (c) ρ component. Figure 28. View largeDownload slide Normalized data residual using constant migration density model. Figure 28. View largeDownload slide Normalized data residual using constant migration density model. 3.6 Noise sensitivity analysis In real cases, the field data are always contaminated by noise, so we add different levels of stochastic noise into the shot records obtained from the SEG/EAGE model. Fig. 29 displays the 160th shot records with an SNR of 30 dB (Fig. 29a), 25 dB (Fig. 29b) and 20 dB (Fig. 29c). The SNR of shot records are calculated by   \begin{equation}{\rm SNR} = \log \left( {\frac{{\sqrt {\frac{{\sum\limits_{j = {j_m} - M/2}^{{j_m} + M/2} {{{\bf d}}({i_m},j)} }}{M}} }}{{\sqrt {\frac{{\sum\limits_{i = 1}^{nx} {\sum\limits_{j = 1}^{nz} {({{{\bf d}}_n}(i,j) - {{\bf d}}(i,j))} } }}{{nx \times nz}}} }}} \right),\end{equation} (16)where d and dn denote the observed shot records without and with noise, respectively, nx and nz are the grid dimensions in x- and z-directions, im and jm are coordinates where the shot amplitude is the largest, and M = λ · Δz, where λ denotes the wavelength. We implement V-SELSRTM using these noisy shot records shown in Fig. 29. In Fig. 30, we display the migrated images after 40 iterations from data with different SNR. It can be observed that all migrated images still have high-resolution and migration artefacts are more and more obvious as the SNR decreases. Fig. 31 shows the convergence curves of the data residual over iterations from data with different SNR. The result shows that when inputting the noisy data set with different SNR, the data residuals of the V-SELSRTM converge to smaller values, suggesting that our V-SELSRTM could work well and be robust for field data. Figure 29. View largeDownload slide The 160th shot records with an SNR of (a) and (d) 30 dB, (b) and (e) 25 dB and (c) and (f) 20 dB. Figure 29. View largeDownload slide The 160th shot records with an SNR of (a) and (d) 30 dB, (b) and (e) 25 dB and (c) and (f) 20 dB. Figure 30. View largeDownload slide The V-SELSRTM images after 40 iterations obtained from shot records with an SNR of (a)–(c) 30 dB, (d)–(f) 25 dB and (g)–(i) 20 dB. (a), (d) and (g) Vp component, (b), (e) and (h) Vs component, and (c), (f) and (i) ρ component. Figure 30. View largeDownload slide The V-SELSRTM images after 40 iterations obtained from shot records with an SNR of (a)–(c) 30 dB, (d)–(f) 25 dB and (g)–(i) 20 dB. (a), (d) and (g) Vp component, (b), (e) and (h) Vs component, and (c), (f) and (i) ρ component. Figure 31. View largeDownload slide Normalized data residual obtained from shot records with an SNR of 30, 25 and 20 dB. The black, red and blue lines denote shot records with an SNR of 30, 25 and 20 dB, respectively. Figure 31. View largeDownload slide Normalized data residual obtained from shot records with an SNR of 30, 25 and 20 dB. The black, red and blue lines denote shot records with an SNR of 30, 25 and 20 dB, respectively. 3.7 Computational cost To compare the computational efficiency of V-SELSRTM to V-ELSRTM, we must consider both the computational cost for each iteration and the convergence rate. All numerical examples are parallelized with MPI. All codes are executed on a Dell T7600 with 32 cores Intel Xeon 2.7 GHz processors and 64 GB of memory. In this test, the iterations will stop if the end conditions are satisfied (the end conditions for different model and data residuals are given in Table 1). The average computational cost of V-SELSRTM for each iteration is 1.87 times that of V-ELSRTM. The reason is that compared with eq. (1), eq. (3) has more variables and computational equations, therefore the computational cost of V-SELSRTM for each iteration is more than that of V-ELSRTM when calculating wavefield forward and backward propagations. However, when the end conditions of model residuals are 0.5 and 0.3, the total computational costs of V-SELSRTM are 0.71 and 0.90 times that of V-ELSRTM, respectively; when the end conditions of data residuals are 10−1.5 and 10−1.8, the total computational costs of V-SELSRTM are 0.52 and 0.43 times that of V-ELSRTM, respectively. V-SELSRTM has faster convergence rate with fewer iterations than V-ELSRTM to satisfy the end conditions because the P and S waves are separated to decrease the crosstalk artefacts between Vp and Vs. Table 1. Computational cost comparison of V-SELSRTM and V-ELSRTM.   Computational cost      Model residual  Data residual  Imaging method  Each iteration (s)  <0.5 (s)  <0.3 (s)  <10−1.5 (s)  <10−1.8 (s)  V-SELSRTM  404  3224  7218  2464  3733  V-ELSRTM  216  4542  7985  4750  8616    Computational cost      Model residual  Data residual  Imaging method  Each iteration (s)  <0.5 (s)  <0.3 (s)  <10−1.5 (s)  <10−1.8 (s)  V-SELSRTM  404  3224  7218  2464  3733  V-ELSRTM  216  4542  7985  4750  8616  View Large 4 CONCLUSIONS Compared with N-ELSRTM based on the non-density-perturbation assumption, V-ELSRTM can produce higher quality images, which is very important for subsequent amplitude versus offset/amplitude versus angle (AVO/AVA) analysis. The numerical experiments using a layered elastic model with different Vp and Vs anomalies demonstrate the advantage of V-ELSRTM over N-ELSRTM. To reduce crosstalk artefacts caused by P- and S-waves coupling in V-ELSRTM, a V-SELSRTM scheme is proposed. The results also suggest that compared with V-ELSRTM, V-SELSRTM can minimize crosstalk artefacts between Vp and Vs because P and S waves are separated to calculate forward- and backward-propagated wavefields and gradient directions. Numerical tests on synthetic data sets from SEG/EAGE salt model further verify that V-SELSRTM can obtain higher resolution and more balanced amplitude images than V-ELSRTM. Tests performed with imperfect migration velocities and constant migration density suggest that the V-SELSRTM can produce acceptable inversion results as long as the migration velocity models are within pre-defined thresholds and it is less sensitive to the errors in the migration density model than the migration velocity models. Sensitivity analysis with data sets with different SNR further verifies the robustness of the V-SELSRTM for field data. Compared with the V-ELSRTM, the V-SELSRTM is more efficient due to its faster convergence rate. V-SELSRTM cannot completely eliminate the crosstalks between Vp and Vs because of the conversion waves, and there still exist obvious crosstalks between velocities and density. Therefore, we can introduce the incorporation of the inverse Hessian operator into this method to mitigate the crosstalks. ACKNOWLEDGEMENTS We would like to thank SWPI group from China University of Petroleum (East China) and Alan Levander from Rice University, and thank the editors and reviewers for their precious suggestions on improving the paper. REFERENCES Aki K., Richards P., 2002. Quantitative Seismology , 2nd edn, University Science Books. Balch A.H., Erdemir C., 1994. Sign-change correction for prestack migration of P-S converted wave reflections1, Geophys. Prospect. , 42, 637– 663. https://doi.org/10.1111/j.1365-2478.1994.tb00233.x Google Scholar CrossRef Search ADS   Chang W., Mcmechan G.A., 1987. Elastic reverse-time migration, Geophysics , 52, 1365– 1375. https://doi.org/10.1190/1.1442249 Google Scholar CrossRef Search ADS   Chang W., Mcmechan G.A., 1994. 3-D elastic prestack, reverse-time depth migration, Geophysics , 59, 597– 609. https://doi.org/10.1190/1.1443620 Google Scholar CrossRef Search ADS   Claerbout J., 1992. Earth Soundings Analysis: Processing Versus Inversion , Blackwell Scientific. Crase E., Pica A., Noble M., Mcdonald J., Tarantola A., 1990. Robust elastic nonlinear waveform inversion: application to real data, Geophysics , 55, 527– 538. https://doi.org/10.1190/1.1442864 Google Scholar CrossRef Search ADS   Dai T., Kuo J.T., 1986. Real data results of Kirchhoff elastic wave migration, Geophysics , 51, 1006– 1011. https://doi.org/10.1190/1.1442139 Google Scholar CrossRef Search ADS   Dai W., Schuster G.T., 2013. Plane-wave least-squares reverse-time migration, Geophysics , 78, S165– S177. https://doi.org/10.1190/geo2012-0377.1 Google Scholar CrossRef Search ADS   Dai W., Fowler P., Schuster G.T., 2012. Multi-source least-squares reverse time migration, Geophys. Prospect. , 60, 681– 695. https://doi.org/10.1111/j.1365-2478.2012.01092.x Google Scholar CrossRef Search ADS   Drossaert F.H., Giannopoulos A., 2007. A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves, Geophysics , 72, T9– T17. https://doi.org/10.1190/1.2424888 Google Scholar CrossRef Search ADS   Du Q., Hou B., 2008. Elastic Kirchhoff migration of vectorial wave-fields, Appl. Geophys. , 5, 284– 293. https://doi.org/10.1007/s11770-008-0045-z Google Scholar CrossRef Search ADS   Du Q., Zhu Y., Ba J., 2012. Polarity reversal correction for elastic reverse time migration, Geophysics , 77, S31– S41. https://doi.org/10.1190/geo2011-0348.1 Google Scholar CrossRef Search ADS   Du Q., Guo C., Zhao Q., Gong X., Wang C., Li X., 2017. Vector-based elastic reverse time migration based on scalar imaging condition, Geophysics , 82, S111– S127. https://doi.org/10.1190/geo2016-0146.1 Google Scholar CrossRef Search ADS   Duan Y, Sava P., Guitton A., 2016. Elastic least-squares reverse time migration, in SEG Technical Program Expanded Abstracts , pp. 4152– 4157. Duquet B., Marfurt K.J., Dellinger J.A., 2000. Kirchhoff modeling, inversion for reflectivity, and subsurface illumination, Geophysics , 65, 1195– 1209. https://doi.org/10.1190/1.1444812 Google Scholar CrossRef Search ADS   Etgen J., Gray S.H., Zhang Y., 2009. An overview of depth imaging in exploration geophysics, Geophysics , 74, WCA5– WCA17. https://doi.org/10.1190/1.3223188 Google Scholar CrossRef Search ADS   Feng Z., Schuster G.T., 2017. Elastic least-squares reverse time migration, Geophysics , 82, S143– S157. https://doi.org/10.1190/geo2016-0254.1 Google Scholar CrossRef Search ADS   Gu B., Li Z., Ma X., Liang G., 2015. Multi-component elastic reverse time migration based on the P- and S-wave separated velocity-stress equations, J. appl. Geophys. , 112, 62– 78. https://doi.org/10.1016/j.jappgeo.2014.11.008 Google Scholar CrossRef Search ADS   Han J., Wang Y., Lu J., 2013. Multi-component Gaussian beam prestack depth migration, J. geophys. Eng. , 10, doi:10.1088/17422132/10/5/055008. https://doi.org/10.1088/1742-2132/10/5/055008 Hokstad K., 2000. Multicomponent Kirchhoff migration, Geophysics , 65, 861– 873. https://doi.org/10.1190/1.1444783 Google Scholar CrossRef Search ADS   Hou J., Symes W.W., 2015. An approximate inverse to the extended born modeling operator, Geophysics , 80, R331– R349. https://doi.org/10.1190/geo2014-0592.1 Google Scholar CrossRef Search ADS   Hou J., Symes W.W., 2016. Accelerating extended least-squares migration with weighted conjugate gradient iteration, Geophysics , 81, S165– S179. https://doi.org/10.1190/geo2015-0499.1 Google Scholar CrossRef Search ADS   Kuehl H., Sacchi M., 1999. Least-squares split-step migration using the Hartley transform, in SEG Technical Program Expanded Abstracts , pp. 1548– 1551. Köhn D., De Nil D., Kurzmann A., Przebindowska A., Bohlen T., 2012. On the influence of model parametrization in elastic full waveform tomography, Geophys. J. Int. , 191, 325– 345. https://doi.org/10.1111/j.1365-246X.2012.05633.x Google Scholar CrossRef Search ADS   Kuo J.T., Dai T., 1984. Kirchhoff elastic wave migration for the case of noncoincident source and receiver, Geophysics , 49, 1223– 1238. https://doi.org/10.1190/1.1441751 Google Scholar CrossRef Search ADS   Li Z., Zhang H., Liu Q., Han W., 2007. Numeric simulation of elastic wavefield separation by staggering grid high-order finite-difference algorithm, Oil Geophys. Prospect. , 42, 510– 515. Ma D., Zhu G., 2003. P- and S-wave separated elastic wave equation numerical modeling, Oil Geophys. Prospect. , 38, 482– 486. Mora P., 1987. Nonlinear two-dimensional elastic inversion of multioffset seismic data, Geophysics , 52, 1211– 1228. https://doi.org/10.1190/1.1442384 Google Scholar CrossRef Search ADS   Métivier L., Brossier R., Operto S., Virieux J., 2014. Multi-parameter FWI – An illustration of the hessian operator role for mitigating trade-off between parameter classes, in 6th EAGE Saint Petersburg International Conference and Exhibition . Nemeth T., Wu C., Schuster G.T., 1999. Least-squares migration of incomplete reflection data, Geophysics , 64, 208– 221. https://doi.org/10.1190/1.1444517 Google Scholar CrossRef Search ADS   Protasov M.I., Tcheverda V.A., 2012. True amplitude elastic Gaussian beam imaging of multicomponent walkaway vertical seismic profiling data, Geophys. Prospect. , 60, 1030– 1042. https://doi.org/10.1111/j.1365-2478.2012.01068.x Google Scholar CrossRef Search ADS   Qu Y., Li Z., Huang J., Li Q., Han Y., 2015a. Pre-stack elastic wave reverse time migration of irregular surface based on layered mapping method, in Near-Surface Asia Pacific Conference , Waikoloa, Hawaii, pp. 104– 107. Qu Y., Li Z., Huang J., Li Q., Yang Y., Li J., 2015b. Elastic wave modeling and wave field separation of irregular free-surface based on multi-block mapping method, in SEG Technical Program Expanded Abstracts , pp. 3724– 3728. Qu Y., Li Z., Huang J., Deng W., Li J., 2015c. The least-squares reverse time migration for viscoacoustic medium based on a stable reverse-time propagator, in SEG Technical Program Expanded Abstracts , pp. 3977– 3980. Qu Y., Li Z., Huang J., Li J., 2016a. ViscoacousticVTI least square reverse time migration, in EAGE Technical Program Expanded Abstracts . Qu Y., Li Z., Huang J., Li J., 2016b. Multi-scale full waveform inversion for areas with irregular surface topography in an auxiliary coordinate system, Explor. Geophys. , doi:10.1071/EG16037. https://doi.org/10.1071/EG16037 Qu Y., Li Z., Huang J., Li J., 2017, Viscoacoustic anisotropic full waveform inversion, J. appl. Geophys. , 136, 484– 497. https://doi.org/10.1016/j.jappgeo.2016.12.001 Google Scholar CrossRef Search ADS   Sun R., Mcmechan G.A., Lee C., Chow J., Chen C., 2006. Prestack scalar reverse-time depth migration of 3D elastic seismic data, Geophysics , 71, S199– S207. https://doi.org/10.1190/1.2227519 Google Scholar CrossRef Search ADS   Sun R., Chow J., Chen K., 2001. Phase correction in separating P-and S-waves in elastic data, Geophysics , 66, 1515– 1518. https://doi.org/10.1190/1.1487097 Google Scholar CrossRef Search ADS   Sun R., Mcmechan G.A., Chuang H., 2011. Amplitude balancing in separating P- and S-waves in 2D and 3D elastic seismic data, Geophysics , 76, S103– S113. https://doi.org/10.1190/1.3555529 Google Scholar CrossRef Search ADS   Tan S., Huang L., 2014. Least-squares reverse-time migration with a wavefield-separation imaging condition and updated source wavefields, Geophysics , 79, S195– S205. https://doi.org/10.1190/geo2014-0020.1 Google Scholar CrossRef Search ADS   Wu R., 1994. Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method, J. geophys. Res. , 99, 751– 766. https://doi.org/10.1029/93JB02518 Google Scholar CrossRef Search ADS   Wu R. S., Xie X. B., 1994. Multi-screen back propagator for fast 3D elastic prestack migration, in Mathematical Methods in Geophysical Imaging II.  pp. 181– 193. International Society for Optical Engineering. Xie X., Wu R., 2005. Multicomponent prestack depth migration using the elastic screen method, Geophysics , 70, S30– S37. https://doi.org/10.1190/1.1852787 Google Scholar CrossRef Search ADS   Xu L., Stanton A., Sacchi M., 2016. Elastic least-squares reverse time migration, in SEG Technical Program Expanded Abstracts , pp. 2289– 2293. Xue Z., Chen Y., Fomel S., Sun J., 2014. Imaging incomplete data and simultaneous-source data using least-squares reverse-time migration with shaping regularization, in SEG Technical Program Expanded Abstracts , pp. 3991– 3996. Xue Z., Chen Y., Fomel S., Sun J., 2016. Seismic imaging of incomplete data and simultaneous-source data using least-squares reverse time migration with shaping regularization, Geophysics , 81, S11– S20. https://doi.org/10.1190/geo2014-0524.1 Google Scholar CrossRef Search ADS   Yan J., Sava P., 2008. Isotropic angle-domain elastic reverse-time migration, Geophysics , 73, S229– S239. https://doi.org/10.1190/1.2981241 Google Scholar CrossRef Search ADS   Yang J., Liu Y., Dong L., 2016. Least-squares reverse time migration in the presence of density variations, Geophysics , 81, S497– S509. https://doi.org/10.1190/geo2016-0075.1 Google Scholar CrossRef Search ADS   Ying-Ming Q., Jian-Ping H., Zhen-Chun L., Qing-Yang L., Jin-Liang Z., Xiu-Zhi L., 2015, Elastic wave modeling and pre-stack reverse time migration of irregular free-surface based on layered mapping method, Chin. J. Geophys. , 58, 544– 560. https://doi.org/10.1002/cjg2.20194 Google Scholar CrossRef Search ADS   Yu P., Geng J., Li X., Wang C., 2016. Acoustic-elastic coupled equation for ocean bottom seismic data elastic reverse time migration, Geophysics , 81, S333– S345. https://doi.org/10.1190/geo2015-0535.1 Google Scholar CrossRef Search ADS   APPENDIX A: VERIFICATION FOR THE P- AND S-WAVE MODES SEPARATED FIRST-ORDER VELOCITY–STRESS EQUATION In this section, we demonstrate the correctness of the P- and S-wave modes separated first-order velocity–stress equation. (1) Verifying that the solutions of eqs (3) satisfy eqs (1). The wavefields are made up of P and S wavefields:   \begin{equation}{{\bf U}} = {{{\bf U}}_p} + {{{\bf U}}_s},\end{equation} (A1)where   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{\bf U}} = \left\{ {{v_x},{v_z},{\tau _{xx}},{\tau _{zz}},{\tau _{xz}}} \right\}\\ {{{\bf U}}_p} = \left\{ {{v_{xp}},{v_{zp}},{\tau _{xxp}},{\tau _{zzp}},{\tau _{xzp}}} \right\}\\ {{{\bf U}}_s} = \left\{ {{v_{xs}},{v_{zs}},{\tau _{xxs}},{\tau _{zzs}},{\tau _{xzs}}} \right\} \end{array} \right..\end{equation} (A2) From the fifth equation in eqs (1), we can obtain that   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\tau _{xzp}} = 0\\ {\tau _{xz}} = {\tau _{xzs}} \end{array} \right..\end{equation} (A3) Adding the first equation in eqs (3b) and (3c) together, and substituting eq. (A3) into them, we can get the first equation in eqs (1). Similarly, we can reach the same conclusion for the second to the fifth equation in eqs (1). So the solutions of eqs (3) satisfy eq. (1). (2) Verifying that eqs (3) can produce pure P and S waves. As is known to us all, P wavefield is curl-free field and S wavefield is divergence-free field. In the first-order velocity–stress equation, only vx and vz are received by receivers though there are five elements. Therefore, we only need to verify that vp(vxp, vzp) is curl-free field and vs(vxs, vzs) is divergence-free fields. We first calculate the divergence of vs under the assumptions of a local homogeneous medium:   \begin{eqnarray} \frac{{{\partial ^2}(\nabla \cdot {{{\bf v}}_s})}}{{\partial {t^2}}} &=& \frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {\frac{{\partial {v_{xs}}}}{{\partial x}} + \frac{{\partial {v_{zs}}}}{{\partial z}}} \right)\nonumber \\ &=& \frac{{{\partial ^2}}}{{\partial x\partial t}}\left( {\frac{{\partial {v_{xs}}}}{{\partial t}}} \right) + \frac{{{\partial ^2}}}{{\partial z\partial t}}\left( {\frac{{\partial {v_{zs}}}}{{\partial t}}} \right)\nonumber \\ &=& \frac{{{\partial ^2}}}{{\partial x\partial t}}\left( {\frac{1}{\rho }\frac{{\partial {\tau _{xxs}}}}{{\partial x}} + \frac{1}{\rho }\frac{{\partial {\tau _{xzs}}}}{{\partial z}}} \right) + \frac{{{\partial ^2}}}{{\partial z\partial t}}\left( {\frac{1}{\rho }\frac{{\partial {\tau _{xzs}}}}{{\partial x}} + \frac{1}{\rho }\frac{{\partial {\tau _{zzs}}}}{{\partial z}}} \right)\nonumber \\ &=& \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial {x^2}}}\left( {\frac{{\partial {\tau _{xxs}}}}{{\partial t}}} \right) + \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial x\partial z}}\left( {\frac{{\partial {\tau _{xzs}}}}{{\partial t}}} \right) + \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial z\partial x}}\left( {\frac{{\partial {\tau _{xzs}}}}{{\partial t}}} \right) + \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial {z^2}}}\left( {\frac{{\partial {\tau _{zzs}}}}{{\partial t}}} \right)\nonumber \\ &=& \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial {x^2}}}\left( { - 2\mu \frac{{\partial {v_z}}}{{\partial z}}} \right) + 2\frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial x\partial z}}\left( {\mu \left(\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}\right)} \right) + \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial {z^2}}}\left( { - 2\mu \frac{{\partial {v_x}}}{{\partial x}}} \right)\nonumber\\ &=& - 2\mu \frac{1}{\rho }\frac{{{\partial ^2}{v_z}}}{{\partial {x^2}\partial z}} + 2\mu \frac{1}{\rho }\frac{{{\partial ^2}{v_x}}}{{\partial x\partial {z^2}}} + 2\mu \frac{1}{\rho }\frac{{{\partial ^2}{v_z}}}{{\partial {x^2}\partial z}} - 2\mu \frac{1}{\rho }\frac{{{\partial ^2}{v_x}}}{{\partial x\partial {z^2}}}\nonumber \\ &=& 0. \end{eqnarray} (A4) From eq. (A4), we conclude that ∇ · vs is a constant or a linear function. According to wave characteristics, ∇ · vs cannot be a linear function, so ∇ · vs can only be a constant. In the initial conditions for solving wave equations, ∇ · vs|t = 0 = 0, so ∇ · vs ≡ 0. Similarly, we can prove that ∇ × vp ≡ 0. Therefore, vp in eq. (3b) is non-curl field, and vs in eq. (3c) is non-divergence field. Eqs (3) can produce pure P and S waves. APPENDIX B: GRADIENTS FOR LAMÉ PARAMETERS AND DENSITY PERTURBATIONS BASED ON COUPLED ELASTIC WAVE EQUATION Based on Born approximation theorem, parameter perturbations   \begin{equation}\delta {{\bf V}}(\lambda ,\mu ,\rho ) = {{\bf V}}(\lambda ,\mu ,\rho ) - {{{\bf V}}_0}(\lambda ,\mu ,\rho ),\end{equation} (B1)indicate wavefield perturbations   \begin{equation}\delta {{\bf U}} = {{\bf U}} - {{{\bf U}}_0},\end{equation} (B2)where U0 and V0 denote the background wavefields and parameters, and δU and δV denote the perturbation wavefields and parameters. U = (vx, vz, τxx, τzz, τxz) and V = (λ, μ, vp, vs, ρ). The background wavefield U0 can be identified by   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial {v_{x0}}}}{{\partial t}} = \frac{{\partial {\tau _{xx0}}}}{{\partial x}} + \frac{{\partial {\tau _{xz0}}}}{{\partial z}}\\ {\rho _0}\frac{{\partial {v_{z0}}}}{{\partial t}} = \frac{{\partial {\tau _{xz0}}}}{{\partial x}} + \frac{{\partial {\tau _{zz0}}}}{{\partial z}}\\ \frac{{\partial {\tau _{xx0}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\frac{{\partial {v_{x0}}}}{{\partial x}} + {\lambda _0}\frac{{\partial {v_{z0}}}}{{\partial z}} + {f_x}\\ \frac{{\partial {\tau _{zz0}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\frac{{\partial {v_{z0}}}}{{\partial z}} + {\lambda _0}\frac{{\partial {v_{x0}}}}{{\partial x}} + {f_z}\\ \frac{{\partial {\tau _{xz0}}}}{{\partial t}} = {\mu _0}\left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right..\end{equation} (B3) Substituting U(x, t) = U0(x, t) + δU(x, t) into eq. (1), subtracting eq. (B3) and neglecting the high-order terms of δV yields:   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_x}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xx}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{xz}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{x0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_z}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xz}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{zz}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{z0}}}}{{\partial t}}\\ \frac{{\partial \delta {\tau _{xx}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\frac{{\partial \delta {v_x}}}{{\partial x}} + {\lambda _0}\frac{{\partial \delta {v_z}}}{{\partial z}} + {{\bf m}}(\lambda ){\lambda _0}\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2{{\bf m}}(\mu ){\mu _0}\frac{{\partial {v_{x0}}}}{{\partial x}}\\ \frac{{\partial \delta {\tau _{zz}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\frac{{\partial \delta {v_z}}}{{\partial z}} + {\lambda _0}\frac{{\partial \delta {v_x}}}{{\partial x}} + {{\bf m}}(\lambda ){\lambda _0} \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2{{\bf m}}(\mu ){\mu _0}\frac{{\partial {v_{z0}}}}{{\partial z}}\\ \frac{{\partial \delta {\tau _{xz}}}}{{\partial t}} = {\mu _0}\left(\frac{{\partial \delta {v_x}}}{{\partial z}} + \frac{{\partial \delta {v_z}}}{{\partial x}}\right) + {{\bf m}}(\mu ){\mu _0}\left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right.,\end{equation} (B4)where the reflectivity models with respect to parameter perturbations are defined as   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{\bf m}}(\lambda ) = \frac{{\delta \lambda }}{{{\lambda _0}}}\\ {{\bf m}}(\mu ) = \frac{{\delta \mu }}{{{\mu _0}}}\\ {{\bf m}}(\rho ) = \frac{{\delta \rho }}{{{\rho _0}}} \end{array} \right..\end{equation} (B5) The perturbation of the objective function δE(m) can be written by   \begin{eqnarray} \delta E({{\bf m}}) = \frac{1}{2}\delta \left\| {{{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}}} \right\|_2^2 &=& \frac{1}{2}\delta \left\langle {{{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}},{{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}}} \right\rangle = \left\langle {\delta ({{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}}),} \right.\nonumber \\ \left. {({{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle &=& \left\langle {\delta ({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}}),({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle = \left\langle {{{\bf B}}\delta {{\bf U}},({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle, \end{eqnarray} (B6)where B denotes a sampling operator to restrict the data in the whole model space onto the receiver locations. From eq. (B4), we get   \begin{equation}\delta {{\bf U}} = {{\bf L}}{{{\bf F}}^\prime},\end{equation} (B7)where F΄ is a new constructed source matrix used to produce wavefield perturbation δU, and given by   \begin{equation}{{{\bf F}}^\prime} = \left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{x0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{z0}}}}{{\partial t}}\\ \delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\\ \delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\\ \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right).\end{equation} (B8) Thus, eq. (B6) can be changed to   \begin{equation}\delta E({{\bf m}}) = \left\langle {{{\bf B}}({{\bf L}}{{{\bf F}}^\prime}),({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle = \left\langle {{{{\bf F}}^\prime},{{{\bf L}}^R}{{{\bf B}}^R}({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle ,\end{equation} (B9)where superscript R denotes the adjoint operator. BR denotes where the data at the receivers extends to the whole model space, LR denotes the wavefield reverse-time backpropagation operator. So δE(λ, μ, ρ) can be written as   \begin{eqnarray} \delta E(\lambda ,\mu ,\rho ) &=& {\left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{x0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{z0}}}}{{\partial t}}\\ \delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\\ \delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\\ \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right)^T}\left( \begin{array}{@{}l@{}} v_x^R\\ v_z^R\\ \tau _{xx}^R\\ \tau _{zz}^R\\ \tau _{xz}^R \end{array} \right)\nonumber\\ &=& - \delta \rho \frac{{\partial {v_{x0}}}{{\partial t}}} v_x^R - \delta \rho \frac{{\partial {v_{z0}}}}{{\partial t}}v_z^R + \left[ {\delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}} \right]\tau _{xx}^R\nonumber\\ &&+\,\, \left[ {\delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}} \right]\tau _{zz}^R + \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right)\tau _{xz}^R.\end{eqnarray} (B10) To minimize the objective function of eq. (1), we obtain the gradient directions of E(m) for different material parameters:   \begin{equation}{{{\bf g}}_\lambda } = \frac{{\partial {{\bf E}}}}{{\partial \lambda }} = \int{{\left( {\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\tau _{xx}^R + \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\tau _{zz}^R} \right)}}dt,\end{equation} (B11)  \begin{equation}{{{\bf g}}_\mu } = \frac{{\partial {{\bf E}}}}{{\partial \mu }} = \int{{\left( {2\frac{{\partial {v_{x0}}}}{{\partial x}}\tau _{xx}^R + 2\frac{{\partial {v_{z0}}}}{{\partial z}}\tau _{zz}^R + \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right)\tau _{xz}^R} \right)}}dt,\end{equation} (B12)  \begin{equation}{{{\bf g}}_\rho } = \frac{{\partial {{\bf E}}}}{{\partial \rho }} = \int{{\left( { - \frac{{\partial {v_{x0}}}}{{\partial t}}v_x^R - \frac{{\partial {v_{z0}}}}{{\partial t}}v_z^R} \right)}}dt.\end{equation} (B13) APPENDIX C: GRADIENTS FOR ELSRTM BASED ON DECOUPLED ELASTIC WAVE EQUATION Rewriting eqs (3b) and (3c) using the relationship $${V_p} = \sqrt {(\lambda + 2\mu )/\rho } $$ and $${V_s} = \sqrt {\mu /\rho } $$ yields   \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_{xp}}}}{{\partial t}} = \frac{{\partial {\tau _{xxp}}}}{{\partial x}}\\ \rho \frac{{\partial {v_{zp}}}}{{\partial t}} = \frac{{\partial {\tau _{zzp}}}}{{\partial z}}\\ \frac{1}{{V_p^2}}\frac{{\partial {\tau _{xxp}}}}{{\partial t}} = \rho \left(\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}\right) + {f_x}\\ \frac{1}{{V_p^2}}\frac{{\partial {\tau _{zzp}}}}{{\partial t}} = \rho \left(\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}\right) + {f_z} \end{array} \right.,\end{equation} (C1)and   \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_{xs}}}}{{\partial t}} = \frac{{\partial {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial {\tau _{xzs}}}}{{\partial z}}\\ \rho \frac{{\partial {v_{zs}}}}{{\partial t}} = \frac{{\partial {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial {\tau _{zzs}}}}{{\partial z}}\\ \frac{1}{{V_s^2}}\frac{{\partial {\tau _{xxs}}}}{{\partial t}} = - 2\rho \frac{{\partial {v_z}}}{{\partial z}} + {f_x}\\ \frac{1}{{V_s^2}}\frac{{\partial {\tau _{zzs}}}}{{\partial t}} = - 2\rho \frac{{\partial {v_x}}}{{\partial x}} + {f_z}\\ \frac{1}{{V_s^2}}\frac{{\partial {\tau _{xzs}}}}{{\partial t}} = \rho \left(\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}\right) \end{array} \right..\end{equation} (C2) Substituting U0 (Up0, Us0) and V0 (Vp0, Vs0, ρ) into eqs (C1) and (C2) to obtain the P- and S-wave background equations. Expanding the Vp and Vs terms according to   \begin{equation}\frac{1}{{V_p^2}} \approx \frac{1}{{V_{p0}^2}} - \frac{{2\delta {V_p}}}{{V_{p0}^3}},\end{equation} (C3)  \begin{equation}\frac{1}{{V_s^2}} \approx \frac{1}{{V_{s0}^2}} - \frac{{2\delta {V_s}}}{{V_{s0}^3}},\end{equation} (C4)and subtracting the P- and S-wave background equations from eqs (C1) and (C2), respectively, yields the equations for δUp and δUs after neglecting the high-order terms yields:   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxp}}}}{{\partial x}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{xp0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{zzp}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{zp0}}}}{{\partial t}}\\ \frac{1}{{V_{p0}^2}}\frac{{\partial \delta {\tau _{xxp}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + \frac{{2{{\bf m}}({V_p})}}{{V_{p0}^2}}\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\\ \frac{1}{{V_{p0}^2}}\frac{{\partial \delta {\tau _{zzp}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + \frac{{2{{\bf m}}({V_p})}}{{V_{p0}^2}}\frac{{\partial {\tau _{zzp0}}}}{{\partial t}} \end{array} \right.,\end{equation} (C5)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{xzs}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{xs0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{zzs}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{zs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{xxs}}}}{{\partial t}} = - 2{\rho _0}\frac{{\partial \delta {v_z}}}{{\partial z}} + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{zzs}}}}{{\partial t}} = - 2{\rho _0}\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{zzs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{xzs}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial z}} + \frac{{\partial \delta {v_z}}}{{\partial x}}\right) + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{xzs0}}}}{{\partial t}} \end{array} \right.,\end{equation} (C6)where the reflectivity models associated with velocities and density perturbations are defined as   \begin{equation}{{\bf m}}({V_p}) = \frac{{\delta {V_p}}}{{{V_{p0}}}},\end{equation} (C7)  \begin{equation}{{\bf m}}({V_s}) = \frac{{\delta {V_s}}}{{{V_{s0}}}},\end{equation} (C8)  \begin{equation}{{\bf m}}(\rho ) = \frac{{\delta \rho }}{\rho }.\end{equation} (C9) For velocities and density perturbations, eq. (4) can be expressed as   \begin{eqnarray} E({V_p},{V_s},\rho ) &=& \left( {{{\left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{xp0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{zp0}}}}{{\partial t}}\\ \frac{{2\delta {V_p}}}{{V_{p0}^3}}\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\\ \frac{{2\delta {V_p}}}{{V_{p0}^3}}\frac{{\partial {\tau _{zzp0}}}}{{\partial t}} \end{array} \right)}^T}\left( \begin{array}{@{}l@{}} v_{xp}^R\\ v_{zp}^R\\ \tau _{xxp}^R\\ \tau _{zzp}^R \end{array} \right) + {{\left( {\begin{array}{@{}*{1}{c}@{}} { - \delta \rho \frac{{\partial {v_{xs0}}}}{{\partial t}}}\\ { - \delta \rho \frac{{\partial {v_{zs0}}}}{{\partial t}}}\\ {\frac{{2\delta {V_s}}}{{V_{s0}^3}}\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}}\\ {\frac{{2\delta {V_s}}}{{V_{s0}^3}}\frac{{\partial {\tau _{zzs0}}}}{{\partial t}}}\\ {\frac{{2\delta {V_s}}}{{V_{s0}^3}}\frac{{\partial {\tau _{xzs0}}}}{{\partial t}}} \end{array}} \right)}^T}\left( \begin{array}{@{}l@{}} v_{xs}^R\\ v_{zs}^R\\ \tau _{xxs}^R\\ \tau _{zzs}^R\\ \tau _{xzs}^R \end{array} \right)} \right)\nonumber\\ &=& \frac{{\delta {V_p}}}{{V_{p0}^3}}\left( {\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\tau _{xxp}^R + \frac{{\partial {\tau _{zzp0}}}}{{\partial t}}\tau _{zzp}^R} \right) + \frac{{\delta {V_s}}}{{V_{s0}^3}}\left( {\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}\tau _{xxs}^R + \frac{{\partial {\tau _{zzs0}}}}{{\partial t}}\tau _{zzs}^R}+\, \frac{{\partial {\tau _{xzs0}}}}{{\partial t}}\tau _{xzs}^R \right) \nonumber\\ &&-\,\, \delta \rho \left( {\frac{{\partial {v_{xp0}}}}{{\partial t}}v_{xp}^R + \frac{{\partial {v_{zp0}}}}{{\partial t}}v_{zp}^R + \frac{{\partial {v_{xs0}}}}{{\partial t}}v_{xs}^R + \frac{{\partial {v_{zs0}}}}{{\partial t}}v_{zs}^R} \right).\end{eqnarray} (C10) The gradient formulae with respect toVp, Vs and ρ can be written as eq. (10) by calculating the partial derivative of Vp, Vs and ρ, respectively. Similarly, when performing ELSRTM for Lamé parameters and density perturbations based on decoupled elastic wave equation, the equations for δUp and δUs are given by   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxp}}}}{{\partial x}} - \delta \rho \frac{{\partial {v_{xp0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{zzp}}}}{{\partial z}} - \delta \rho \frac{{\partial {v_{zp0}}}}{{\partial t}}\\ \frac{{\partial \delta {\tau _{xxp}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\\ \frac{{\partial \delta {\tau _{zzp}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) \end{array} \right.,\end{equation} (C11)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{xzs}}}}{{\partial z}} - \delta \rho \frac{{\partial {v_{xs0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{zzs}}}}{{\partial z}} - \delta \rho \frac{{\partial {v_{zs0}}}}{{\partial t}}\\ \frac{{\partial \delta {\tau _{xxs}}}}{{\partial t}} = - 2{\mu _0}\frac{{\partial \delta {v_z}}}{{\partial z}} - 2\delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\\ \frac{{\partial \delta {\tau _{zzs}}}}{{\partial t}} = - 2{\mu _0}\frac{{\partial \delta {v_x}}}{{\partial x}} - 2\delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\\ \frac{{\partial \delta {\tau _{xzs}}}}{{\partial t}} = {\mu _0}\left(\frac{{\partial \delta {v_x}}}{{\partial z}} + \frac{{\partial \delta {v_z}}}{{\partial x}}\right) + \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right..\end{equation} (C12) For velocities and density perturbations, eq. (4) can be written as   \begin{eqnarray} {{\bf E}}(\lambda ,\mu ,\rho ) &=& {\left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{xp0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{zp0}}}}{{\partial t}}\\ (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\\ (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) \end{array} \right)^T}\left( \begin{array}{@{}l@{}} v_{xp}^R\\ v_{zp}^R\\ \tau _{xxp}^R\\ \tau _{zzp}^R \end{array} \right) + {\left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{xs0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{zs0}}}}{{\partial t}}\\ - \delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\\ - \delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\\ \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right)^T}\left( \begin{array}{@{}l@{}} v_{xs}^R\\ v_{zs}^R\\ \tau _{xxs}^R\\ \tau _{zzs}^R\\ \tau _{xzs}^R \end{array} \right)\nonumber\\ &=& - \delta \rho \frac{{\partial {v_{xp0}}}}{{\partial t}} v_{xp}^R - \delta \rho \frac{{\partial {v_{zp0}}}}{{\partial t}}v_{zp}^R + (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\left(\tau _{xxp}^R + \tau _{zzp}^R\right)\nonumber\\ &&-\,\,\delta \rho \frac{{\partial {v_{xs0}}}}{{\partial t}}v_{xs}^R - \delta \rho \frac{{\partial {v_{zs0}}}}{{\partial t}}v_{zs}^R - \delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\tau _{xxs}^R - \delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\tau _{zzs}^R + \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right)\tau _{xzs}^R. \end{eqnarray} (C13) Therefore, the gradients with respect to the Lamé parameters and density can be accessed as eq. (14) by calculating the partial derivative of λ, μ and ρ, respectively. © Crown copyright 2017. This article contains public sector information licensed under the Open Government Licence v3.0 (http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Elastic least-squares reverse time migration with velocities and density perturbation

Loading next page...
 
/lp/ou_press/elastic-least-squares-reverse-time-migration-with-velocities-and-w6Giql0nh4
Publisher
Oxford University Press
Copyright
© Crown copyright 2017.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx468
Publisher site
See Article on Publisher Site

Abstract

Summary Elastic least-squares reverse time migration (LSRTM) based on the non-density-perturbation assumption can generate false-migrated interfaces caused by density variations. We perform an elastic LSRTM scheme with density variations for multicomponent seismic data to produce high-quality images in Vp, Vs and ρ components. However, the migrated images may suffer from crosstalk artefacts caused by P- and S-waves coupling in elastic LSRTM no matter what model parametrizations used. We have proposed an elastic LSRTM with density variations method based on wave modes separation to reduce these crosstalk artefacts by using P- and S-wave decoupled elastic velocity–stress equations to derive demigration equations and gradient formulae with respect to Vp, Vs and ρ. Numerical experiments with synthetic data demonstrate the capability and superiority of the proposed method. The imaging results suggest that our method promises imaging results with higher quality and has a faster residual convergence rate. Sensitivity analysis of migration velocity, migration density and stochastic noise verifies the robustness of the proposed method for field data. Image processing, Inverse theory, Seismic tomography, Wave propagation 1 INTRODUCTION With the development of multicomponent seismic exploration, conventional seismic imaging methods based on the acoustic approximation appear to be imperfect, because the assumption that only P waves propagate in the earth is not valid theoretically. Therefore, elastic imaging technologies attract more and more attention of exploration geophysicists. Several types of elastic imaging methods are used for multicomponent seismic data, mainly including ray-based migration methods, such as elastic Kirchhoff migration (Kuo & Dai 1984; Dai & Kuo 1986; Hokstad 2000; Du & Hou 2008) and Gaussian beam migration (Protasov & Tcheverda 2012; Han et al.2013), one-way wave equation migration (Wu & Xie 1994; Wu 1994; Xie & Wu 2005) and reverse time migration (RTM; Chang & McMechan 1987, 1994; Yan & Sava 2008; Qu et al.2015a,b; Yu et al.2016). In these methods, elastic RTM (ERTM) is the most promising and popular migration method because it is free from dip-angle limitation and high-frequency assumption. The coupling between P and S waves can generate crosstalk artefacts, resulting in inaccurate imaging results. To reduce crosstalk artefacts, ERTM algorithms use wave modes separation to produce scale images. Scalar RTM algorithm is first used for multicomponent data (Sun et al.2001, 2006). This method separates the multicomponent data into P- and S-component ones by calculating their divergence and curl, and uses them to implement acoustic RTM. The scalar RTM, ignoring the elastic characteristics of the vector wavefields, cannot accurately migrate converted waves. Therefore, vector ERTM methods, separating P and S waves when calculating forward and backward wavefields, have been developed over the past decade to produce PP, PS-, SP- and SS-wave images (Yan & Sava 2008). To obtain higher quality PS and SP images, polarity corrections need to be applied to S waves (Balch & Erdemir 1994; Du et al.2012), which cannot be accurately implemented in complex structures (Du et al.2017). In isotropic media, wave modes separation can be easily implemented by employing Helmholtz decomposition (Aki & Richards 2002; Qu et al.2015c) to the vector wavefields. But the calculation of divergence and curl operator will lead to a phase shift (Sun et al.2001) and amplitude changes (Sun et al. 2011). Recently, ERTM methods based on decoupled wave equations (Ma & Zhu 2003; Li et al.2007) are proposed (Gu et al.2015; Du et al.2017). Conventional migration applies the adjoint of a forward modeling operator to observed data instead of its inverse (Claerbout 1992; Nemeth et al.1999). A key problem with conventional migration is that it can produce a low-quality image with data acquisition footprint because of aliased seismic data, poor acquisition geometry, irregular data sampling, and a narrow source–receiver aperture (Kuehl & Sacchi 1999; Duquet et al.2000; Etgen et al.2009). Compared with the conventional RTM method, least-squares reverse time migration (LSRTM) is a powerful tool to produce a high-quality image with higher spatial resolution, fewer migration artefacts and better balanced reflector amplitudes (Dai & Schuster 2013; Tan & Huang 2014; Qu et al.2015c, 2016a; Hou & Symes 2015, 2016; Xue et al.2014, 2016; Yang et al.2016). With the development of computer technologies, elastic LSRTM (ELSRTM) is being applied to multicomponent seismic data (Duan et al.2016; Xu et al.2016; Feng & Schuster 2017). Ren et al. (2017) analyse the influence of different model parametrizations and objective functions in ELSRTM. Compared with acoustic LSRTM, ELSRTM produces better resolution images. These ELSRTM methods are all based on the non-density-perturbation assumption, resulting in the false-migrated interfaces caused by density variations. To overcome this problem, an ELSRTM scheme with density variations is developed. The approach is called V-ELSRTM. However, the crosstalk artefacts caused by the P and S waves coupling in multiparameter inversion are still a problem no matter what model parametrizations used (Feng & Schuster 2017). Métivier et al. (2014) incorporated the inverse Hessian operator within inversion method to alleviate this difficulty. We propose a new V-ELSRTM method based on P- and S-wavefield separation (V-SELSRTM) to reduce the crosstalk artefacts and enhance imaging quality. 2 ELSRTM FOR VELOCITIES AND DENSITY PERTURBATIONS BASED ON DECOUPLED ELASTIC WAVE EQUATION In the time domain, the velocity–stress wave equation with variable density is written as:   \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_x}}}{{\partial t}} = \frac{{\partial {\tau _{xx}}}}{{\partial x}} + \frac{{\partial {\tau _{xz}}}}{{\partial z}}\\ \rho \frac{{\partial {v_z}}}{{\partial t}} = \frac{{\partial {\tau _{xz}}}}{{\partial x}} + \frac{{\partial {\tau _{zz}}}}{{\partial z}}\\ \frac{{\partial {\tau _{xx}}}}{{\partial t}} = (\lambda + 2\mu )\frac{{\partial {v_x}}}{{\partial x}} + \lambda \frac{{\partial {v_z}}}{{\partial z}} + {f_x}\\ \frac{{\partial {\tau _{zz}}}}{{\partial t}} = (\lambda + 2\mu )\frac{{\partial {v_z}}}{{\partial z}} + \lambda \frac{{\partial {v_x}}}{{\partial x}} + {f_z}\\ \frac{{\partial {\tau _{xz}}}}{{\partial t}} = \mu \left(\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}\right) \end{array} \right.,\end{equation} (1)where vx and vz are the velocities in horizontal and vertical components, respectively, τxx and τzz are the normal stress, τxz is the shear stress, fx and fz are the two components of the source terms, t is time, ρ is the density, λ and μ are the Lamé parameters. The relationships between the Lamé parameters λ and μ, density ρ, P-wave velocity Vp, and S wave velocity Vs are   \begin{equation}\left\{ \begin{array}{@{}l@{}} \lambda = \rho (V_p^2 - 2V_s^2)\\ \mu = \rho V_s^2 \end{array} \right..\end{equation} (2) To reduce crosstalk artefacts of V-ELSRTM caused by the coupling of different wave modes, we develop a V-SELSRTM method for velocities and density perturbations based on the decoupled elastic wave equation. The decoupled elastic velocity–stress wave equations can separate coupled elastic waves into P- and S-mode waves, and are given by (Li et al.2007)   \begin{equation}\left\{ \begin{array}{@{}l@{}} {v_x} = {v_{xp}} + {v_{xs}}\\ {v_z} = {v_{zp}} + {v_{zs}} \end{array} \right.,\end{equation} (3a)  \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_{xp}}}}{{\partial t}} = \frac{{\partial {\tau _{xxp}}}}{{\partial x}}\\ \rho \frac{{\partial {v_{zp}}}}{{\partial t}} = \frac{{\partial {\tau _{zzp}}}}{{\partial z}}\\ \frac{{\partial {\tau _{xxp}}}}{{\partial t}} = (\lambda + 2\mu )\left(\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}\right) + {f_x}\\ \frac{{\partial {\tau _{zzp}}}}{{\partial t}} = (\lambda + 2\mu )\left(\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}\right) + {f_z} \end{array} \right.,\end{equation} (3b)  \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_{xs}}}}{{\partial t}} = \frac{{\partial {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial {\tau _{xzs}}}}{{\partial z}}\\ \rho \frac{{\partial {v_{zs}}}}{{\partial t}} = \frac{{\partial {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial {\tau _{zzs}}}}{{\partial z}}\\ \frac{{\partial {\tau _{xxs}}}}{{\partial t}} = - 2\mu \frac{{\partial {v_z}}}{{\partial z}} + {f_x}\\ \frac{{\partial {\tau _{zzs}}}}{{\partial t}} = - 2\mu \frac{{\partial {v_x}}}{{\partial x}} + {f_z}\\ \frac{{\partial {\tau _{xzs}}}}{{\partial t}} = \mu \left(\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}\right) \end{array} \right.,\end{equation} (3c)where the subscript p and s denote the P-wavefield Up and S-wavefield Us, respectively. Eqs (3a)–(3c) are mixed, P-wave and S-wave equations, respectively. Source terms fx and fz are added in both P- and S-wave equations. Only when shear sources are excited, source term is added only in the fifth equation of S-wave equations (3c). The correctness of the P- and S-wave modes separated first-order velocity–stress wave equation is demonstrated in Appendix A. Least-squares migration can be seen as a method to solve an optimization problem. The purpose is to find a perturbation model in the least-squares sense to fit the observed data by minimizing the following objective function:   \begin{equation}E({{\bf m}}) = \frac{1}{2}\left\| {{{{\bf d}}_{{\rm{cal}}}} - {{{\bf d}}_{{\rm{obs}}}}} \right\|_2^2,\end{equation} (4)where dcal denote the synthetic data and dcal = (dcalp, dcals), and dobs denote the observed data and dobs = (dobsp, dobss); the synthetic data can be given by dcal = Lm, in which L denotes a linear modeling operator and L = (Lp, Ls), m is the reflectivity model. Perturbation wavefields δUp and δUs can be obtained by solving (details see Appendix C)   \begin{equation}\left\{ \begin{array}{@{}l@{}} \delta {v_x} = \delta {v_{xp}} + \delta {v_{xs}}\\ \delta {v_z} = \delta {v_{zp}} + \delta {v_{zs}} \end{array} \right.,\end{equation} (5a)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxp}}}}{{\partial x}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{xp0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{zzp}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{zp0}}}}{{\partial t}}\\ \frac{1}{{V_{p0}^2}}\frac{{\partial \delta {\tau _{xxp}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + \frac{{2{{\bf m}}({V_p})}}{{V_{p0}^2}}\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\\ \frac{1}{{V_{p0}^2}}\frac{{\partial \delta {\tau _{zzp}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + \frac{{2{{\bf m}}({V_p})}}{{V_{p0}^2}}\frac{{\partial {\tau _{zzp0}}}}{{\partial t}} \end{array} \right.,\end{equation} (5b)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{xzs}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{xs0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{zzs}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{zs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{xxs}}}}{{\partial t}} = - 2{\rho _0}\frac{{\partial \delta {v_z}}}{{\partial z}} + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{zzs}}}}{{\partial t}} = - 2{\rho _0}\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{zzs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{xzs}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial z}} + \frac{{\partial \delta {v_z}}}{{\partial x}}\right) + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{xzs0}}}}{{\partial t}} \end{array} \right.,\end{equation} (5c)where Up0 and Us0 are the solutions associated with the background parameters V0 of eq. (3), δUp and δUs denote the perturbation wavefields, in which Up = (vxp, vzp, τxxp, τzzp) and Us= (vxs, vzs, τxxs, τzzs, τxzs), and δV is parameter perturbations and V = (Vp, Vs, ρ). The reflectivity models associated with Vp, Vs and ρ are defined by   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{\bf m}}({V_p}) = \frac{{\delta {V_p}}}{{{V_{p0}}}}\\ {{\bf m}}({V_s}) = \frac{{\delta {V_s}}}{{{V_{s0}}}}\\ {{\bf m}}(\rho ) = \frac{{\delta \rho }}{{{\rho _0}}} \end{array} \right..\end{equation} (6) The P- and S-wave migration equations can also be given by   \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial v_{xp}^R}}{{\partial t}} = \frac{{\partial \tau _{xxp}^R}}{{\partial x}} + \Delta {{{\bf d}}_x}({{{\bf x}}_r},t)\\ \rho \frac{{\partial v_{zp}^R}}{{\partial t}} = \frac{{\partial \tau _{zzp}^R}}{{\partial z}} + \Delta {{{\bf d}}_z}({{{\bf x}}_r},t)\\ \frac{{\partial \tau _{xxp}^R}}{{\partial t}} = \rho V_p^2\left(\frac{{\partial v_x^R}}{{\partial x}} + \frac{{\partial v_z^R}}{{\partial z}}\right)\\ \frac{{\partial \tau _{zzp}^R}}{{\partial t}} = \rho V_p^2\left(\frac{{\partial v_x^R}}{{\partial x}} + \frac{{\partial v_z^R}}{{\partial z}}\right) \end{array} \right.,\end{equation} (7a)  \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial v_{xs}^R}}{{\partial t}} = \frac{{\partial \tau _{xxs}^R}}{{\partial x}} + \frac{{\partial \tau _{xzs}^R}}{{\partial z}} + \Delta {{{\bf d}}_x}({{{\bf x}}_r},t)\\ \rho \frac{{\partial v_{zs}^R}}{{\partial t}} = \frac{{\partial \tau _{xzs}^R}}{{\partial x}} + \frac{{\partial \tau _{zzs}^R}}{{\partial z}} + \Delta {{{\bf d}}_z}({{{\bf x}}_r},t)\\ \frac{{\partial \tau _{xxs}^R}}{{\partial t}} = - 2\rho V_s^2\frac{{\partial v_z^R}}{{\partial z}}\\ \frac{{\partial \tau _{zzs}^R}}{{\partial t}} = - 2\rho V_s^2\frac{{\partial v_x^R}}{{\partial x}}\\ \frac{{\partial \tau _{xzs}^R}}{{\partial t}} = \rho V_s^2\left(\frac{{\partial v_x^R}}{{\partial z}} + \frac{{\partial v_z^R}}{{\partial x}}\right) \end{array} \right.,\end{equation} (7b)where $${{{\bf U}}^R} = ({{\bf U}}_p^R,{{\bf U}}_s^R)$$ denotes the adjoint matrix of U = (Up, Us), xr denotes the location of receivers, Δdx(xr, t) and Δdz(xr, t) are the difference between observed data and synthetic data in horizontal and vertical directions, respectively, which are calculated by   \begin{equation}\left\{ \begin{array}{@{}l@{}} \Delta {{{\bf d}}_x}({{{\bf x}}_r},t) = {{{\bf d}}_{{\rm{cal}}x}}({{{\bf x}}_r},t) - {{{\bf d}}_{{\rm{obs}}x}}({{{\bf x}}_r},t)\\ \Delta {{{\bf d}}_z}({{{\bf x}}_r},t) = {{{\bf d}}_{{\rm{cal}}z}}({{{\bf x}}_r},t) - {{{\bf d}}_{{\rm{obs}}z}}({{{\bf x}}_r},t) \end{array} \right.,\end{equation} (8)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {{{\bf d}}_{{\rm{cal}}x}}({{{\bf x}}_r},t) = \delta {v_x}({{{\bf x}}_r},t)\\ {{{\bf d}}_{{\rm{cal}}z}}({{{\bf x}}_r},t) = \delta {v_z}({{{\bf x}}_r},t) \end{array} \right..\end{equation} (9) Noted that the input observed data are multicomponent shot records dobs (dobsx, dobsz) without separating P and S wavefields. The predicted data are mixed multicomponent shot records dcal (dcalx, dcalz), which are mixed by eq. (3a). And then the residual P and S wavefields are separated during the backpropagation. Then the gradient formulae with respect to Vp, Vs and ρ can be calculated by (details see Appendix C)   \begin{eqnarray}\left\{ \begin{array}{@{}l@{}} {{{\bf g}}_{vp}} = \frac{{\partial {{\bf E}}}}{{\partial {V_p}}} = \frac{1}{{V_{p0}^3}}\int{{\left( {\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\tau _{xxp}^R + \frac{{\partial {\tau _{zzp0}}}}{{\partial t}}\tau _{zzp}^R} \right)}}dt\\ {{{\bf g}}_{vs}} = \frac{{\partial {{\bf E}}}}{{\partial {V_s}}} = \frac{1}{{V_{s0}^3}}\int{{\left( {\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}\tau _{xxs}^R + \frac{{\partial {\tau _{zzs0}}}}{{\partial t}}\tau _{zzs}^R + \frac{{\partial {\tau _{xzs0}}}}{{\partial t}}\tau _{xzs}^R} \right)}}dt\\ {{{\bf g}}_{\rho v}} = \frac{{\partial {{\bf E}}}}{{\partial \rho }} = \int{{\left( {\frac{{\partial {v_{xp0}}}}{{\partial t}}v_{xp}^R + \frac{{\partial {v_{zp0}}}}{{\partial t}}v_{zp}^R + \frac{{\partial {v_{xs0}}}}{{\partial t}}v_{xs}^R + \frac{{\partial {v_{zs0}}}}{{\partial t}}v_{zs}^R} \right)}}dt \end{array} \right.,\nonumber\\ \end{eqnarray} (10)where gvp, gvs and gρv are gradient formulae with respect to Vp, Vs and ρ, respectively. With the relationships between velocities (Vp and Vs) and impedances (P impedance Ip and S impedance Is)   \begin{equation}\left\{ \begin{array}{@{}l@{}} {I_p} = \rho {V_p}\\ {I_s} = \rho {V_s} \end{array} \right.,\end{equation} (11)the reflectivity of Ip and Iscan be written by   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{\bf m}}({I_p}) = \frac{{\delta {I_p}}}{{{I_{p0}}}} = \frac{{\delta (\rho {V_p})}}{{{\rho _0}{V_{p0}}}} = \frac{{\delta {V_p}}}{{{v_{p0}}}} + \frac{{\delta \rho }}{{{\rho _0}}} = {{\bf m}}({V_p}) + {{\bf m}}(\rho )\\ {{\bf m}}({I_s}) = \frac{{\delta {I_s}}}{{{I_{s0}}}} = \frac{{\delta (\rho {V_s})}}{{{\rho _0}{V_{s0}}}} = \frac{{\delta {V_s}}}{{{v_{s0}}}} + \frac{{\delta \rho }}{{{\rho _0}}} = {{\bf m}}({V_s}) + {{\bf m}}(\rho ) \end{array} \right..\end{equation} (12) The gradient formulae with respect to Ip, Is and ρ can be obtained by applying the chain rule (Mora 1987; Crase et al.1990; Köhn et al.2012):   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{{\bf g}}_{Ip}} = \frac{1}{\rho }{{{\bf g}}_{vp}}\\ {{{\bf g}}_{Is}} = \frac{1}{\rho }{{{\bf g}}_{vs}}\\ {{{\bf g}}_{\rho i}} = - \frac{{{V_p}}}{\rho }{{{\bf g}}_{vp}} - \frac{{{V_s}}}{\rho }{{{\bf g}}_{vs}} + {{{\bf g}}_{\rho v}} \end{array} \right..\end{equation} (13) Similarly, when performing V-ELSRTM for Lamé parameters and density perturbations based on decoupled elastic wave equation, the gradients related to λ and μ, ρ can be accessed as (Appendix C)   \begin{eqnarray}\left\{ \begin{array}{@{}l@{}} {{{\bf g}}_\lambda } = \int{{\left( {\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)(\tau _{xxp}^R + \tau _{zzp}^R)} \right)}}dt\\ {{{\bf g}}_\mu } = \int\left( 2\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)(\tau _{xxp}^R + \tau _{zzp}^R) - \frac{{\partial {v_{z0}}}}{{\partial z}}\tau _{xxs}^R - \frac{{\partial {v_{x0}}}}{{\partial x}}\tau _{zzs}^R \right. \\ \qquad \qquad\left.+\left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right)\tau _{xzs}^R \right)dt\\ {{{\bf g}}_\rho } = \int{{\left( {\frac{{\partial {v_{xp0}}}}{{\partial t}}v_{xp}^R + \frac{{\partial {v_{zp0}}}}{{\partial t}}v_{zp}^R + \frac{{\partial {v_{xs0}}}}{{\partial t}}v_{xs}^R + \frac{{\partial {v_{zs0}}}}{{\partial t}}v_{zs}^R} \right)}}dt \end{array} \right.\end{eqnarray} (14) The gradients for velocities and impedances can also be calculated by applying the chain rule. In this paper, we use eq. (10) to update V-SELSRTM images. The pre-conditioned conjugate gradient (Dai et al.2012) method is applied to solve the least-squares inverse problem. The calculation flow of elastic LSRTM with velocities and density perturbation is illustrated in Fig. 1. The left- and right-hand columns show the conventional V-ELSRTM (details see Appendix B) and the V-SELSRTM, respectively. Figure 1. View largeDownload slide The calculation flow of the conventional V-ELSRTM (the left-hand column) and V-SELSRTM (the right-hand column). Figure 1. View largeDownload slide The calculation flow of the conventional V-ELSRTM (the left-hand column) and V-SELSRTM (the right-hand column). 3 NUMERICAL EXAMPLES First, the performances of the P- and S-wave modes separated elastic wave equation and linear modeling operators are demonstrated through a numerical example using a simple two-layer model. Then the application of the proposed V-SELSRTM compared with the conventional V-ELSRTM and ELSRTM based on the non-density-perturbation assumption (N-ELSRTM, Feng & Schuster 2017) is demonstrated with synthetic data examples using two models: (1) a layered elastic model with different Vp and Vs anomalies and (2) the SEG/EAGE salt model. In these synthetic examples, the observed data are produced by linear forward modeling, which has variable density. All our simulations are implemented using a time-domain staggered-grid finite-difference algorithm with the second-order accurate in time and tenth-order-accurate in space scheme. Multisource LSRTM is introduced to compress the seismic data and improve the efficiency of LSRTM (Dai et al.2012). Dynamic encoding technology is used to suppress the crosstalks caused by different sources in one supershot (Qu et al.2016b, 2017), which are different from those caused by coupling of P and S waves. To suppress artificial boundary reflections, complex frequency-shifted perfectly matched layer (CFS-PML) boundary conditions (Drossaert & Giannopoulos 2007) are used in all four edges of the models. 3.1 Wavefield separation test with a simple two-layer model We consider an isotropic elastic two-layer model characterized by the P velocity, S velocity and density illustrated in Figs 2(a)–(c) to demonstrate the effectiveness of the P- and S-wave modes separated elastic wave equation and linear modeling operators. The flat interfaces of the P velocity, S velocity and density models are located at different depths. The physical property parameters are presented in Fig. 2. We choose the two-layer model because the P and S waves have different traveltime, and therefore we can determine their individual amplitudes and phases without mutual interference. The dimensions of the model are 201 × 201 gridpoints in the horizontal and vertical directions, with 8 m grid spacing. An explosive source using the Ricker wavelet with a dominant frequency of 25 Hz is excited at (800 m, 8 m), and the total recording length is 1.5 s, with a time sampling of 0.5 ms. The 201 geophones are spread out over the surface. Fig. 3 shows the snapshots at 240 and 320 ms generated using the P- and S-wave modes separated elastic wave equations (3). A comparison of the mixed- and separated-mode wavefield indicates that the wave modes separation method works well for isotropic media. The reflectivity models of the two-layer model (calculated by eq. 6) are displayed in Fig. 4, which are used to perform linear modeling. Fig. 5 shows mixed synthetic shot records obtained from eq. (B4) (Figs 5a and d) and separated synthetic shot records obtained from eq. (5) (Figs 5b,c,e and f). These results suggest that the primary reflected P and S waves are simulated accurately and separated completely and there are no change in amplitude and phase for both P and S waves after wavefield separation, demonstrating the correctness of eq. (5) and the effectiveness of the linear modeling for V-SELSERM. Figure 2. View largeDownload slide A two-layer model. (a) Vp, (b) Vs and (c) ρ. Figure 2. View largeDownload slide A two-layer model. (a) Vp, (b) Vs and (c) ρ. Figure 3. View largeDownload slide Snapshots at (a)–(f) 240 ms and (g)–(l) 320 ms. (a) and (g) Mixed wavefield in x-direction; (b) and (h) P wavefield in x-direction; (c) and (i) S wavefield record in x-direction; (d) and (j) mixed wavefield in z-direction; (e) and (k) P wavefield in z-direction; (f) and (l) S wavefield in z-direction. Figure 3. View largeDownload slide Snapshots at (a)–(f) 240 ms and (g)–(l) 320 ms. (a) and (g) Mixed wavefield in x-direction; (b) and (h) P wavefield in x-direction; (c) and (i) S wavefield record in x-direction; (d) and (j) mixed wavefield in z-direction; (e) and (k) P wavefield in z-direction; (f) and (l) S wavefield in z-direction. Figure 4. View largeDownload slide The reflectivity of the two-layer model. (a) Vp, (b) Vs and (c) ρ. Figure 4. View largeDownload slide The reflectivity of the two-layer model. (a) Vp, (b) Vs and (c) ρ. Figure 5. View largeDownload slide (a) and (d) Mixed synthetic shot records obtained from eq. (B4) and (b), (c), (e) and (f) separated synthetic shot records obtained from eq. (5). (a)–(c) x-direction and (d)–(f) z-direction. Figure 5. View largeDownload slide (a) and (d) Mixed synthetic shot records obtained from eq. (B4) and (b), (c), (e) and (f) separated synthetic shot records obtained from eq. (5). (a)–(c) x-direction and (d)–(f) z-direction. 3.2 A layered elastic model with different Vp and Vs anomalies The second model considered here is a layered elastic model with different Vp and Vs anomalies (Fig. 6). We use this model to demonstrate the advantages of V-ELSRTM and V-SELSRTM over N-ELSRTM. The velocity models have three flat layers with a distance of 4 km × 2 km in x- and z-directions, discretized on a 401 × 201 grids with a spatial grid interval of 10 m. Two anomalies with Vp and Vs of 4000 and 2309 m s−1 are inserted into different locations of the second layer, the size of which is 1.0 km × 0.4 km in x- and z-directions. The material properties of the elastic model are indicated in Fig. 6. The migration Vp, Vs and density models (Figs 6a, c and e) are smoothed versions of the true ones (Figs 6b,d and f). The wide-aperture survey of generating the synthetic data is made up of 200 explosive P-wave sources spaced every 20 m and located in a depth of 10 m, which emit Ricker wavelets with the dominant frequency of 25 Hz. These shots are compressed into 10 supershots by using dynamic encodes. The synthetic records are received by 401 two-component receivers, spaced every 10 m on the surface, calculated up to 1.6 s with a time sample rate of 0.8 ms. Fig. 7 presents the 101st shot gather with removed direct waves recorded at surface. Figure 6. View largeDownload slide An elastic layered model. (a) True Vp, (b) migration Vp, (c) true Vs, (d) migration Vs, (e) true ρ and (e) migration ρ. Figure 6. View largeDownload slide An elastic layered model. (a) True Vp, (b) migration Vp, (c) true Vs, (d) migration Vs, (e) true ρ and (e) migration ρ. Figure 7. View largeDownload slide A shot gather (shot 101). (a) x-direction and (b) z-direction. Figure 7. View largeDownload slide A shot gather (shot 101). (a) x-direction and (b) z-direction. We first perform the conventional N-ELSRTM on the synthetic data sets generated with variable density. Fig. 8 gives the N-ELSRTM images of the 40th iteration after Laplacian filtering for Vp and Vs components. The ELSRTM image in Vs component has higher resolution than that in Vp component because S wave has shorter wavelength than P wave. However, in the N-ELSRTM images, there are false interfaces caused by the variable density, which are indicated by the red ellipses. Then, we conduct the three V-ELSRTM methods using the same synthetic shot records on the layered model. Figs 9 and 10 display the forward- and backward-propagated wavefields at 400 ms, in which the mixed and separated wavefields are generated by the solutions of eqs (3) and (7), respectively. The results further verify that the P and S waves separated forward modeling method can simulate pure P and S wavefields during seismic wave forward and backward propagation, thus providing a theoretical basis for the proposed V-SELSRTM method. Figure 8. View largeDownload slide Comparisons of (a) the true Vp perturbation, (b) the true Vs perturbation, (c) conventional N-ELSRTM images after 40 iterations for Vp component and (d) conventional N-ELSRTM for Vs component. Figure 8. View largeDownload slide Comparisons of (a) the true Vp perturbation, (b) the true Vs perturbation, (c) conventional N-ELSRTM images after 40 iterations for Vp component and (d) conventional N-ELSRTM for Vs component. Figure 9. View largeDownload slide Forward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 9. View largeDownload slide Forward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 10. View largeDownload slide Backward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 10. View largeDownload slide Backward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Fig. 11 compares the V-ELSRTM and V-SELSRTM images of the 40th iteration after Laplacian filtering. For Vp-, Vs- and ρ-component images, the amplitude is balanced with high resolution and signal-to-noise ratio (SNR). The false interfaces caused by the perturbation of velocities (indicated by the blue ellipses) or density (indicated by the red ellipses) gradually weakened and disappeared against the iteration number, which is demonstrated by the fact that the V-ELSRTM and V-SELSRTM images for Vp, Vs and ρ components have the same pattern as the true Vp (Fig. 11a), Vs (Fig. 11b) and ρ perturbations (Fig. 11c). By the comparisons of the V-ELSRTM and V-SELSRTM images, there are some artefacts in V-ELSRTM images caused by P- and S-wave crosstalks (shown by the red arrows in Figs 11d and e). These crosstalk artefacts are mitigated when using the V-SELSRTM (Figs 11g and h). Fig. 12 presents the convergence curves of model (Fig. 12a) and data residuals (Fig. 12b) over iterations for N-ELSRTM, V-ELSRTM and V-SELSRTM. We track the data residuals using the L2 norm of mixed data residuals as the conventional ELSRTM. It can be observed that in the V-ELSRTM and V-SELSRTM cases, these residuals converge to much smaller values than the N-ELSRTM case. The fluctuations of the residuals are as a result of the dynamic encodes change. In particular, after 10 iterations, the model and data residuals with the N-ELSRTM are clearly larger than those of V-ELSRTM and V-SELSRTM, and the convergence curve of the N-ELSRTM becomes flat after 25 iterations. The curves of V-SELSRTM have faster convergence speeds, and the model and data residuals converge to smaller minimums than V-ELSRTM. Fig. 13 compares the corresponding V-ELSRTM and V-SELSRTM images using seismic impedances (Ip, Is and ρ) with the true impedance perturbations (m(Ip), m(Is) and m(ρ)), in which the false interfaces caused by the perturbation of velocities are indicated by the blue ellipses. The results also suggest that the V-SELSRTM images of impedances show better accordance with the true perturbations of impedances than the V-ELSRTM by reducing P- and S-waves coupling crosstalks (shown by the red arrows). Figure 11. View largeDownload slide (a) Comparisons of the true Vp perturbation, (b) the true Vs perturbation , (c) the true ρ perturbation, (d) V-ELSRTM images for Vp component, (e) Vs component and (f) ρ component, (g) V-SELSRTM images for Vp component, (h) Vs component and ρ component. Figure 11. View largeDownload slide (a) Comparisons of the true Vp perturbation, (b) the true Vs perturbation , (c) the true ρ perturbation, (d) V-ELSRTM images for Vp component, (e) Vs component and (f) ρ component, (g) V-SELSRTM images for Vp component, (h) Vs component and ρ component. Figure 12. View largeDownload slide (a) Normalized model residuals and (b) data residuals. The black lines indicate the residuals obtained from the conventional N-ELSRTM, the red lines indicate the residuals obtained from V-ELSRTM and the blue lines indicate the residuals obtained from V-SELSRTM. Figure 12. View largeDownload slide (a) Normalized model residuals and (b) data residuals. The black lines indicate the residuals obtained from the conventional N-ELSRTM, the red lines indicate the residuals obtained from V-ELSRTM and the blue lines indicate the residuals obtained from V-SELSRTM. Figure 13. View largeDownload slide (a) Comparisons of the true Ip perturbation, (b) the true Is perturbation, (c) the true ρ perturbation, (d) V-ELSRTM images for Ip component, (e) Is component, (f) ρ component, (g) V-SELSRTM images for Ip component, (h) Is component and ρ component. Figure 13. View largeDownload slide (a) Comparisons of the true Ip perturbation, (b) the true Is perturbation, (c) the true ρ perturbation, (d) V-ELSRTM images for Ip component, (e) Is component, (f) ρ component, (g) V-SELSRTM images for Ip component, (h) Is component and ρ component. 3.3 SEG/EAGE salt model To further analyse the capacity of our V-SELSRTM for imaging complex model, we apply the new algorithm to a SEG/EAGE salt model (Fig. 14), discretized on an 645 × 150 grids with a grid spacing of Δx = Δz = 8 m. In the fixed-spread acquisition geometry, we evenly deploy 320 sources on the surface, with the interval between neighbouring sources of 16 m. These shots are compressed into 16 supershots by using dynamic encodes. The explosive point source chooses a Ricker wavelet with the dominant frequency of 20 Hz. There are 645 fixed two-component receivers with a spacing of 8 m. To ensure numerical stability, sample every 0.8 ms across a 3.5 s recording period. There are 50 CFS-PML grids surrounding the SEG/EAGE model boundaries. Figs 15 and 16 illustrate the snapshot comparisons of forward- and backward-propagated wavefields, respectively, after 0.4 s. Figure 14. View largeDownload slide The SEG/EAGE salt model. (a) True Vp, (b) migration Vp, (c) true Vs, (d) migration Vs, (e) true ρ and (f) migration ρ. Figure 14. View largeDownload slide The SEG/EAGE salt model. (a) True Vp, (b) migration Vp, (c) true Vs, (d) migration Vs, (e) true ρ and (f) migration ρ. Figure 15. View largeDownload slide Forward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 15. View largeDownload slide Forward-propagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 16. View largeDownload slide Backpropagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. Figure 16. View largeDownload slide Backpropagated wavefields at 400 ms. (a) Mixed wavefield in x-direction; (b) P wavefield in x-direction; (c) S wavefield record in x-direction; (d) mixed wavefield in z-direction; (e) P wavefield in z-direction and (f) S wavefield in z-direction. The true reflectivity models are displayed in Fig. 17, which are used only for comparison. Fig. 18 plots the V-ELSRTM and V-SELSRTM images after 40 iterations. Similar to the previous example, the Vs-component images show higher resolution than the Vp-component images. The ρ-component images achieve lower SNR than Vp- and Vs-component images, which have some low-frequency noises. V-ELSRTM images show the same resolution as that of V-SELSRTM. However, the V-SELSRTM produces much fewer image artefacts than V-ELSRTM. To compare the resolution and SNR of the salt and subsalt areas in both images, we illustrate the magnified images in Figs 19 and 20. The red and blue boxes in Fig. 18 show the regions for zoom views. The magnified views suggest that the V-ELSRTM method produces more serious ringing artefacts around the reflecting interfaces and the V-SELSRTM method generates better images of the salt and subsalt structures with fewer ringing and low-frequency artefacts. Fig. 21 gives V-SELSRTM images after Laplacian filtering, showing high-quality images with high resolution, fewer artefacts and balanced amplitude. The vertical profiles of the true reflectivity and migration images at the horizontal distance of 3.36 km are shown in Fig. 22. It is obviously observed that the Vp-, Vs- and ρ-component images obtained from the V-SELSRTM method are more accurate than those obtained from the V-ELSRTM method. Fig. 23 presents the convergence curve of data and model residuals for V-ELSRTM and V-SELSRTM. By comparison, V-ELSRTM has a lower convergence rate and residual of V-ELSRTM converges to a larger value than that of the V-SELSRTM. Figure 17. View largeDownload slide The true reflectivity models. (a) The Vp perturbation, (b) the Vs perturbation and (c) the ρ perturbation. Figure 17. View largeDownload slide The true reflectivity models. (a) The Vp perturbation, (b) the Vs perturbation and (c) the ρ perturbation. Figure 18. View largeDownload slide (a) Comparisons of V-ELSRTM images after 40 iterations for Vp component, (c) Vs component, (e) ρ component, (b) V-SELSRTM images for Vp component, (d) Vs component and (f) ρ component. Figure 18. View largeDownload slide (a) Comparisons of V-ELSRTM images after 40 iterations for Vp component, (c) Vs component, (e) ρ component, (b) V-SELSRTM images for Vp component, (d) Vs component and (f) ρ component. Figure 19. View largeDownload slide Magnified views of Fig. 18 shown by the red boxes. (a) V-ELSRTM images for Vp component, (b) Vs component, (c) ρ component, (d) V-SELSRTM images for Vp component, (e) Vs component and (f) ρ component. Figure 19. View largeDownload slide Magnified views of Fig. 18 shown by the red boxes. (a) V-ELSRTM images for Vp component, (b) Vs component, (c) ρ component, (d) V-SELSRTM images for Vp component, (e) Vs component and (f) ρ component. Figure 20. View largeDownload slide Magnified views of Fig. 18 shown by the blue boxes. (a) V-ELSRTM images for Vp component, (b) Vs component, (c) ρ component, (d) V-SELSRTM images for Vp component, (e) Vs component and (f) ρ component. Figure 20. View largeDownload slide Magnified views of Fig. 18 shown by the blue boxes. (a) V-ELSRTM images for Vp component, (b) Vs component, (c) ρ component, (d) V-SELSRTM images for Vp component, (e) Vs component and (f) ρ component. Figure 21. View largeDownload slide V-SELSRTM images after Laplacian filtering for (a) Vp component, (b) Vs component and (c) ρ component. Figure 21. View largeDownload slide V-SELSRTM images after Laplacian filtering for (a) Vp component, (b) Vs component and (c) ρ component. Figure 22. View largeDownload slide Vertical profile at the horizontal distance of 3.36 km. (a) Vp component, (b) Vs component and (c) ρ component. The solid red lines indicate the true reflectivity model, the dashed black lines denote V-ELSRTM images after 40 iterations and the dashed blue lines denote V-SELSRTM images after 40 iterations. Figure 22. View largeDownload slide Vertical profile at the horizontal distance of 3.36 km. (a) Vp component, (b) Vs component and (c) ρ component. The solid red lines indicate the true reflectivity model, the dashed black lines denote V-ELSRTM images after 40 iterations and the dashed blue lines denote V-SELSRTM images after 40 iterations. Figure 23. View largeDownload slide (a) Normalized residuals, where the black line with circles denote the model residual for V-ELSRTM, the red line with circles denote the model residual for V-SELSRTM, the black line with stars denote the data residual for V-ELSRTM and the red line with stars denote the data residual for V-SELSRTM. (b) Normalized data residuals in log domain. Figure 23. View largeDownload slide (a) Normalized residuals, where the black line with circles denote the model residual for V-ELSRTM, the red line with circles denote the model residual for V-SELSRTM, the black line with stars denote the data residual for V-ELSRTM and the red line with stars denote the data residual for V-SELSRTM. (b) Normalized data residuals in log domain. 3.4 Sensitivity analysis to migration velocity error In real cases, migration velocity errors should be taken into account to produce high-quality images. To study the sensitivity of our V-SELSRTM to migration Vp and Vs errors, we smooth the migration Vp and Vs models using different smoothing parameter values and add small errors to them (shown in Fig. 24), with the assumption that the migration density is correct (Fig. 14f). For quantification, we give the relative velocity error by using the following equation   \begin{equation}{\rm{ERROR}} = \frac{{||{{{\bf m}}_{{\rm{true}}}} - {{{\bf m}}_{{\rm{init}}}}||}}{{||{{{\bf m}}_{{\rm{true}}}}||}}.\end{equation} (15) Figure 24. View largeDownload slide Migration velocities with different error of (a) and (b) 7.1 per cent, (c) and (d) 10.8 per cent, (e) and (f) 13.2 per cent and (g) and (h) 16.5 per cent. (a), (c), (e) and (g) migration Vp. (b), (d), (f) and (h) migration Vs. Figure 24. View largeDownload slide Migration velocities with different error of (a) and (b) 7.1 per cent, (c) and (d) 10.8 per cent, (e) and (f) 13.2 per cent and (g) and (h) 16.5 per cent. (a), (c), (e) and (g) migration Vp. (b), (d), (f) and (h) migration Vs. The relative velocity errors between the migration Vp and Vs models and the original ones are 7.1 per cent (Figs 24a, b), 10.8 per cent (Figs 24c and d), 13.2 per cent (Figs 24e and f) and 16.5 per cent (Figs 24g and h). The V-SELSRTM images obtained with the different wrong migration Vp and Vs models are displayed in Fig. 25. Because of the errors in Vp and Vs models, there are some spurious images and obvious low-frequency noises in the shallow and the salt areas, and the subsalt zone is weakly imaged. The wrong velocities would image the reflection data towrong positions. The ρ-component images are more seriously influenced by the velocity errors. As the errors in Vp and Vs models increase, the quality of the V-SELSRTM images gradually decreases. The convergence curves of the data residual over iterations when using imperfect migration velocities are illustrated in Fig. 26. The convergence curves for the four cases convergence to different values with different convergence rate, which suggest that our V-SELSRTM method can produce acceptable inversion results as long as the migration velocity models are within pre-defined thresholds. Figure 25. View largeDownload slide (a)–(c) The V-SELSRTM images after 40 iterations using different migration velocities with the error of 7.1 per cent, (d)–(f)10.8 per cent, (g)–(i) 13.2 per cent and (j)–(l) 16.5 per cent. (a), (d), (g) and (j) Vp component, (b), (e), (h) and (k) Vs component, (c), (f), (i) and (l) ρ component. Figure 25. View largeDownload slide (a)–(c) The V-SELSRTM images after 40 iterations using different migration velocities with the error of 7.1 per cent, (d)–(f)10.8 per cent, (g)–(i) 13.2 per cent and (j)–(l) 16.5 per cent. (a), (d), (g) and (j) Vp component, (b), (e), (h) and (k) Vs component, (c), (f), (i) and (l) ρ component. Figure 26. View largeDownload slide Normalized data residual using migration velocities with different error. The black, red, blue and pink lines denote migration velocities with the error of 7.1 per cent, 10.8 per cent, 13.2 per cent and 16.5 per cent, respectively. Figure 26. View largeDownload slide Normalized data residual using migration velocities with different error. The black, red, blue and pink lines denote migration velocities with the error of 7.1 per cent, 10.8 per cent, 13.2 per cent and 16.5 per cent, respectively. 3.5 Migration density error sensitivity analysis We then study the sensitivity of our V-SELSRTM to the density error in the migration density model. In real cases, the density model is difficult to obtain accurately. Therefore in this case, we use a constant migration density model (the value is set to be 2411 kg m−3 by averaging the true density values) and the accurate migration velocity models (Figs 14b and d) to conduct the V-SELSRTM. Fig. 27 shows the V-SELSRTM images after 40 iterations using a constant migration densitymodel. The ρ-component image shows serious low-frequency noise and weak-imaging energy beneath subsalt areas, in turn resulting in some artefacts on migrated VP- and VS-component images. Fig. 28 also displays the convergence curve of data residuals. The data residuals of V-SELSRTM decrease over iterations and converge to 10−2. In general, the V-SELSRTM is less sensitive to the errors in the migration density model than the migration velocity models. Figure 27. View largeDownload slide The V-SELSRTM images after 40 iterations using constant migration density model. (a) Vp component, (b) Vs component and (c) ρ component. Figure 27. View largeDownload slide The V-SELSRTM images after 40 iterations using constant migration density model. (a) Vp component, (b) Vs component and (c) ρ component. Figure 28. View largeDownload slide Normalized data residual using constant migration density model. Figure 28. View largeDownload slide Normalized data residual using constant migration density model. 3.6 Noise sensitivity analysis In real cases, the field data are always contaminated by noise, so we add different levels of stochastic noise into the shot records obtained from the SEG/EAGE model. Fig. 29 displays the 160th shot records with an SNR of 30 dB (Fig. 29a), 25 dB (Fig. 29b) and 20 dB (Fig. 29c). The SNR of shot records are calculated by   \begin{equation}{\rm SNR} = \log \left( {\frac{{\sqrt {\frac{{\sum\limits_{j = {j_m} - M/2}^{{j_m} + M/2} {{{\bf d}}({i_m},j)} }}{M}} }}{{\sqrt {\frac{{\sum\limits_{i = 1}^{nx} {\sum\limits_{j = 1}^{nz} {({{{\bf d}}_n}(i,j) - {{\bf d}}(i,j))} } }}{{nx \times nz}}} }}} \right),\end{equation} (16)where d and dn denote the observed shot records without and with noise, respectively, nx and nz are the grid dimensions in x- and z-directions, im and jm are coordinates where the shot amplitude is the largest, and M = λ · Δz, where λ denotes the wavelength. We implement V-SELSRTM using these noisy shot records shown in Fig. 29. In Fig. 30, we display the migrated images after 40 iterations from data with different SNR. It can be observed that all migrated images still have high-resolution and migration artefacts are more and more obvious as the SNR decreases. Fig. 31 shows the convergence curves of the data residual over iterations from data with different SNR. The result shows that when inputting the noisy data set with different SNR, the data residuals of the V-SELSRTM converge to smaller values, suggesting that our V-SELSRTM could work well and be robust for field data. Figure 29. View largeDownload slide The 160th shot records with an SNR of (a) and (d) 30 dB, (b) and (e) 25 dB and (c) and (f) 20 dB. Figure 29. View largeDownload slide The 160th shot records with an SNR of (a) and (d) 30 dB, (b) and (e) 25 dB and (c) and (f) 20 dB. Figure 30. View largeDownload slide The V-SELSRTM images after 40 iterations obtained from shot records with an SNR of (a)–(c) 30 dB, (d)–(f) 25 dB and (g)–(i) 20 dB. (a), (d) and (g) Vp component, (b), (e) and (h) Vs component, and (c), (f) and (i) ρ component. Figure 30. View largeDownload slide The V-SELSRTM images after 40 iterations obtained from shot records with an SNR of (a)–(c) 30 dB, (d)–(f) 25 dB and (g)–(i) 20 dB. (a), (d) and (g) Vp component, (b), (e) and (h) Vs component, and (c), (f) and (i) ρ component. Figure 31. View largeDownload slide Normalized data residual obtained from shot records with an SNR of 30, 25 and 20 dB. The black, red and blue lines denote shot records with an SNR of 30, 25 and 20 dB, respectively. Figure 31. View largeDownload slide Normalized data residual obtained from shot records with an SNR of 30, 25 and 20 dB. The black, red and blue lines denote shot records with an SNR of 30, 25 and 20 dB, respectively. 3.7 Computational cost To compare the computational efficiency of V-SELSRTM to V-ELSRTM, we must consider both the computational cost for each iteration and the convergence rate. All numerical examples are parallelized with MPI. All codes are executed on a Dell T7600 with 32 cores Intel Xeon 2.7 GHz processors and 64 GB of memory. In this test, the iterations will stop if the end conditions are satisfied (the end conditions for different model and data residuals are given in Table 1). The average computational cost of V-SELSRTM for each iteration is 1.87 times that of V-ELSRTM. The reason is that compared with eq. (1), eq. (3) has more variables and computational equations, therefore the computational cost of V-SELSRTM for each iteration is more than that of V-ELSRTM when calculating wavefield forward and backward propagations. However, when the end conditions of model residuals are 0.5 and 0.3, the total computational costs of V-SELSRTM are 0.71 and 0.90 times that of V-ELSRTM, respectively; when the end conditions of data residuals are 10−1.5 and 10−1.8, the total computational costs of V-SELSRTM are 0.52 and 0.43 times that of V-ELSRTM, respectively. V-SELSRTM has faster convergence rate with fewer iterations than V-ELSRTM to satisfy the end conditions because the P and S waves are separated to decrease the crosstalk artefacts between Vp and Vs. Table 1. Computational cost comparison of V-SELSRTM and V-ELSRTM.   Computational cost      Model residual  Data residual  Imaging method  Each iteration (s)  <0.5 (s)  <0.3 (s)  <10−1.5 (s)  <10−1.8 (s)  V-SELSRTM  404  3224  7218  2464  3733  V-ELSRTM  216  4542  7985  4750  8616    Computational cost      Model residual  Data residual  Imaging method  Each iteration (s)  <0.5 (s)  <0.3 (s)  <10−1.5 (s)  <10−1.8 (s)  V-SELSRTM  404  3224  7218  2464  3733  V-ELSRTM  216  4542  7985  4750  8616  View Large 4 CONCLUSIONS Compared with N-ELSRTM based on the non-density-perturbation assumption, V-ELSRTM can produce higher quality images, which is very important for subsequent amplitude versus offset/amplitude versus angle (AVO/AVA) analysis. The numerical experiments using a layered elastic model with different Vp and Vs anomalies demonstrate the advantage of V-ELSRTM over N-ELSRTM. To reduce crosstalk artefacts caused by P- and S-waves coupling in V-ELSRTM, a V-SELSRTM scheme is proposed. The results also suggest that compared with V-ELSRTM, V-SELSRTM can minimize crosstalk artefacts between Vp and Vs because P and S waves are separated to calculate forward- and backward-propagated wavefields and gradient directions. Numerical tests on synthetic data sets from SEG/EAGE salt model further verify that V-SELSRTM can obtain higher resolution and more balanced amplitude images than V-ELSRTM. Tests performed with imperfect migration velocities and constant migration density suggest that the V-SELSRTM can produce acceptable inversion results as long as the migration velocity models are within pre-defined thresholds and it is less sensitive to the errors in the migration density model than the migration velocity models. Sensitivity analysis with data sets with different SNR further verifies the robustness of the V-SELSRTM for field data. Compared with the V-ELSRTM, the V-SELSRTM is more efficient due to its faster convergence rate. V-SELSRTM cannot completely eliminate the crosstalks between Vp and Vs because of the conversion waves, and there still exist obvious crosstalks between velocities and density. Therefore, we can introduce the incorporation of the inverse Hessian operator into this method to mitigate the crosstalks. ACKNOWLEDGEMENTS We would like to thank SWPI group from China University of Petroleum (East China) and Alan Levander from Rice University, and thank the editors and reviewers for their precious suggestions on improving the paper. REFERENCES Aki K., Richards P., 2002. Quantitative Seismology , 2nd edn, University Science Books. Balch A.H., Erdemir C., 1994. Sign-change correction for prestack migration of P-S converted wave reflections1, Geophys. Prospect. , 42, 637– 663. https://doi.org/10.1111/j.1365-2478.1994.tb00233.x Google Scholar CrossRef Search ADS   Chang W., Mcmechan G.A., 1987. Elastic reverse-time migration, Geophysics , 52, 1365– 1375. https://doi.org/10.1190/1.1442249 Google Scholar CrossRef Search ADS   Chang W., Mcmechan G.A., 1994. 3-D elastic prestack, reverse-time depth migration, Geophysics , 59, 597– 609. https://doi.org/10.1190/1.1443620 Google Scholar CrossRef Search ADS   Claerbout J., 1992. Earth Soundings Analysis: Processing Versus Inversion , Blackwell Scientific. Crase E., Pica A., Noble M., Mcdonald J., Tarantola A., 1990. Robust elastic nonlinear waveform inversion: application to real data, Geophysics , 55, 527– 538. https://doi.org/10.1190/1.1442864 Google Scholar CrossRef Search ADS   Dai T., Kuo J.T., 1986. Real data results of Kirchhoff elastic wave migration, Geophysics , 51, 1006– 1011. https://doi.org/10.1190/1.1442139 Google Scholar CrossRef Search ADS   Dai W., Schuster G.T., 2013. Plane-wave least-squares reverse-time migration, Geophysics , 78, S165– S177. https://doi.org/10.1190/geo2012-0377.1 Google Scholar CrossRef Search ADS   Dai W., Fowler P., Schuster G.T., 2012. Multi-source least-squares reverse time migration, Geophys. Prospect. , 60, 681– 695. https://doi.org/10.1111/j.1365-2478.2012.01092.x Google Scholar CrossRef Search ADS   Drossaert F.H., Giannopoulos A., 2007. A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves, Geophysics , 72, T9– T17. https://doi.org/10.1190/1.2424888 Google Scholar CrossRef Search ADS   Du Q., Hou B., 2008. Elastic Kirchhoff migration of vectorial wave-fields, Appl. Geophys. , 5, 284– 293. https://doi.org/10.1007/s11770-008-0045-z Google Scholar CrossRef Search ADS   Du Q., Zhu Y., Ba J., 2012. Polarity reversal correction for elastic reverse time migration, Geophysics , 77, S31– S41. https://doi.org/10.1190/geo2011-0348.1 Google Scholar CrossRef Search ADS   Du Q., Guo C., Zhao Q., Gong X., Wang C., Li X., 2017. Vector-based elastic reverse time migration based on scalar imaging condition, Geophysics , 82, S111– S127. https://doi.org/10.1190/geo2016-0146.1 Google Scholar CrossRef Search ADS   Duan Y, Sava P., Guitton A., 2016. Elastic least-squares reverse time migration, in SEG Technical Program Expanded Abstracts , pp. 4152– 4157. Duquet B., Marfurt K.J., Dellinger J.A., 2000. Kirchhoff modeling, inversion for reflectivity, and subsurface illumination, Geophysics , 65, 1195– 1209. https://doi.org/10.1190/1.1444812 Google Scholar CrossRef Search ADS   Etgen J., Gray S.H., Zhang Y., 2009. An overview of depth imaging in exploration geophysics, Geophysics , 74, WCA5– WCA17. https://doi.org/10.1190/1.3223188 Google Scholar CrossRef Search ADS   Feng Z., Schuster G.T., 2017. Elastic least-squares reverse time migration, Geophysics , 82, S143– S157. https://doi.org/10.1190/geo2016-0254.1 Google Scholar CrossRef Search ADS   Gu B., Li Z., Ma X., Liang G., 2015. Multi-component elastic reverse time migration based on the P- and S-wave separated velocity-stress equations, J. appl. Geophys. , 112, 62– 78. https://doi.org/10.1016/j.jappgeo.2014.11.008 Google Scholar CrossRef Search ADS   Han J., Wang Y., Lu J., 2013. Multi-component Gaussian beam prestack depth migration, J. geophys. Eng. , 10, doi:10.1088/17422132/10/5/055008. https://doi.org/10.1088/1742-2132/10/5/055008 Hokstad K., 2000. Multicomponent Kirchhoff migration, Geophysics , 65, 861– 873. https://doi.org/10.1190/1.1444783 Google Scholar CrossRef Search ADS   Hou J., Symes W.W., 2015. An approximate inverse to the extended born modeling operator, Geophysics , 80, R331– R349. https://doi.org/10.1190/geo2014-0592.1 Google Scholar CrossRef Search ADS   Hou J., Symes W.W., 2016. Accelerating extended least-squares migration with weighted conjugate gradient iteration, Geophysics , 81, S165– S179. https://doi.org/10.1190/geo2015-0499.1 Google Scholar CrossRef Search ADS   Kuehl H., Sacchi M., 1999. Least-squares split-step migration using the Hartley transform, in SEG Technical Program Expanded Abstracts , pp. 1548– 1551. Köhn D., De Nil D., Kurzmann A., Przebindowska A., Bohlen T., 2012. On the influence of model parametrization in elastic full waveform tomography, Geophys. J. Int. , 191, 325– 345. https://doi.org/10.1111/j.1365-246X.2012.05633.x Google Scholar CrossRef Search ADS   Kuo J.T., Dai T., 1984. Kirchhoff elastic wave migration for the case of noncoincident source and receiver, Geophysics , 49, 1223– 1238. https://doi.org/10.1190/1.1441751 Google Scholar CrossRef Search ADS   Li Z., Zhang H., Liu Q., Han W., 2007. Numeric simulation of elastic wavefield separation by staggering grid high-order finite-difference algorithm, Oil Geophys. Prospect. , 42, 510– 515. Ma D., Zhu G., 2003. P- and S-wave separated elastic wave equation numerical modeling, Oil Geophys. Prospect. , 38, 482– 486. Mora P., 1987. Nonlinear two-dimensional elastic inversion of multioffset seismic data, Geophysics , 52, 1211– 1228. https://doi.org/10.1190/1.1442384 Google Scholar CrossRef Search ADS   Métivier L., Brossier R., Operto S., Virieux J., 2014. Multi-parameter FWI – An illustration of the hessian operator role for mitigating trade-off between parameter classes, in 6th EAGE Saint Petersburg International Conference and Exhibition . Nemeth T., Wu C., Schuster G.T., 1999. Least-squares migration of incomplete reflection data, Geophysics , 64, 208– 221. https://doi.org/10.1190/1.1444517 Google Scholar CrossRef Search ADS   Protasov M.I., Tcheverda V.A., 2012. True amplitude elastic Gaussian beam imaging of multicomponent walkaway vertical seismic profiling data, Geophys. Prospect. , 60, 1030– 1042. https://doi.org/10.1111/j.1365-2478.2012.01068.x Google Scholar CrossRef Search ADS   Qu Y., Li Z., Huang J., Li Q., Han Y., 2015a. Pre-stack elastic wave reverse time migration of irregular surface based on layered mapping method, in Near-Surface Asia Pacific Conference , Waikoloa, Hawaii, pp. 104– 107. Qu Y., Li Z., Huang J., Li Q., Yang Y., Li J., 2015b. Elastic wave modeling and wave field separation of irregular free-surface based on multi-block mapping method, in SEG Technical Program Expanded Abstracts , pp. 3724– 3728. Qu Y., Li Z., Huang J., Deng W., Li J., 2015c. The least-squares reverse time migration for viscoacoustic medium based on a stable reverse-time propagator, in SEG Technical Program Expanded Abstracts , pp. 3977– 3980. Qu Y., Li Z., Huang J., Li J., 2016a. ViscoacousticVTI least square reverse time migration, in EAGE Technical Program Expanded Abstracts . Qu Y., Li Z., Huang J., Li J., 2016b. Multi-scale full waveform inversion for areas with irregular surface topography in an auxiliary coordinate system, Explor. Geophys. , doi:10.1071/EG16037. https://doi.org/10.1071/EG16037 Qu Y., Li Z., Huang J., Li J., 2017, Viscoacoustic anisotropic full waveform inversion, J. appl. Geophys. , 136, 484– 497. https://doi.org/10.1016/j.jappgeo.2016.12.001 Google Scholar CrossRef Search ADS   Sun R., Mcmechan G.A., Lee C., Chow J., Chen C., 2006. Prestack scalar reverse-time depth migration of 3D elastic seismic data, Geophysics , 71, S199– S207. https://doi.org/10.1190/1.2227519 Google Scholar CrossRef Search ADS   Sun R., Chow J., Chen K., 2001. Phase correction in separating P-and S-waves in elastic data, Geophysics , 66, 1515– 1518. https://doi.org/10.1190/1.1487097 Google Scholar CrossRef Search ADS   Sun R., Mcmechan G.A., Chuang H., 2011. Amplitude balancing in separating P- and S-waves in 2D and 3D elastic seismic data, Geophysics , 76, S103– S113. https://doi.org/10.1190/1.3555529 Google Scholar CrossRef Search ADS   Tan S., Huang L., 2014. Least-squares reverse-time migration with a wavefield-separation imaging condition and updated source wavefields, Geophysics , 79, S195– S205. https://doi.org/10.1190/geo2014-0020.1 Google Scholar CrossRef Search ADS   Wu R., 1994. Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method, J. geophys. Res. , 99, 751– 766. https://doi.org/10.1029/93JB02518 Google Scholar CrossRef Search ADS   Wu R. S., Xie X. B., 1994. Multi-screen back propagator for fast 3D elastic prestack migration, in Mathematical Methods in Geophysical Imaging II.  pp. 181– 193. International Society for Optical Engineering. Xie X., Wu R., 2005. Multicomponent prestack depth migration using the elastic screen method, Geophysics , 70, S30– S37. https://doi.org/10.1190/1.1852787 Google Scholar CrossRef Search ADS   Xu L., Stanton A., Sacchi M., 2016. Elastic least-squares reverse time migration, in SEG Technical Program Expanded Abstracts , pp. 2289– 2293. Xue Z., Chen Y., Fomel S., Sun J., 2014. Imaging incomplete data and simultaneous-source data using least-squares reverse-time migration with shaping regularization, in SEG Technical Program Expanded Abstracts , pp. 3991– 3996. Xue Z., Chen Y., Fomel S., Sun J., 2016. Seismic imaging of incomplete data and simultaneous-source data using least-squares reverse time migration with shaping regularization, Geophysics , 81, S11– S20. https://doi.org/10.1190/geo2014-0524.1 Google Scholar CrossRef Search ADS   Yan J., Sava P., 2008. Isotropic angle-domain elastic reverse-time migration, Geophysics , 73, S229– S239. https://doi.org/10.1190/1.2981241 Google Scholar CrossRef Search ADS   Yang J., Liu Y., Dong L., 2016. Least-squares reverse time migration in the presence of density variations, Geophysics , 81, S497– S509. https://doi.org/10.1190/geo2016-0075.1 Google Scholar CrossRef Search ADS   Ying-Ming Q., Jian-Ping H., Zhen-Chun L., Qing-Yang L., Jin-Liang Z., Xiu-Zhi L., 2015, Elastic wave modeling and pre-stack reverse time migration of irregular free-surface based on layered mapping method, Chin. J. Geophys. , 58, 544– 560. https://doi.org/10.1002/cjg2.20194 Google Scholar CrossRef Search ADS   Yu P., Geng J., Li X., Wang C., 2016. Acoustic-elastic coupled equation for ocean bottom seismic data elastic reverse time migration, Geophysics , 81, S333– S345. https://doi.org/10.1190/geo2015-0535.1 Google Scholar CrossRef Search ADS   APPENDIX A: VERIFICATION FOR THE P- AND S-WAVE MODES SEPARATED FIRST-ORDER VELOCITY–STRESS EQUATION In this section, we demonstrate the correctness of the P- and S-wave modes separated first-order velocity–stress equation. (1) Verifying that the solutions of eqs (3) satisfy eqs (1). The wavefields are made up of P and S wavefields:   \begin{equation}{{\bf U}} = {{{\bf U}}_p} + {{{\bf U}}_s},\end{equation} (A1)where   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{\bf U}} = \left\{ {{v_x},{v_z},{\tau _{xx}},{\tau _{zz}},{\tau _{xz}}} \right\}\\ {{{\bf U}}_p} = \left\{ {{v_{xp}},{v_{zp}},{\tau _{xxp}},{\tau _{zzp}},{\tau _{xzp}}} \right\}\\ {{{\bf U}}_s} = \left\{ {{v_{xs}},{v_{zs}},{\tau _{xxs}},{\tau _{zzs}},{\tau _{xzs}}} \right\} \end{array} \right..\end{equation} (A2) From the fifth equation in eqs (1), we can obtain that   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\tau _{xzp}} = 0\\ {\tau _{xz}} = {\tau _{xzs}} \end{array} \right..\end{equation} (A3) Adding the first equation in eqs (3b) and (3c) together, and substituting eq. (A3) into them, we can get the first equation in eqs (1). Similarly, we can reach the same conclusion for the second to the fifth equation in eqs (1). So the solutions of eqs (3) satisfy eq. (1). (2) Verifying that eqs (3) can produce pure P and S waves. As is known to us all, P wavefield is curl-free field and S wavefield is divergence-free field. In the first-order velocity–stress equation, only vx and vz are received by receivers though there are five elements. Therefore, we only need to verify that vp(vxp, vzp) is curl-free field and vs(vxs, vzs) is divergence-free fields. We first calculate the divergence of vs under the assumptions of a local homogeneous medium:   \begin{eqnarray} \frac{{{\partial ^2}(\nabla \cdot {{{\bf v}}_s})}}{{\partial {t^2}}} &=& \frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {\frac{{\partial {v_{xs}}}}{{\partial x}} + \frac{{\partial {v_{zs}}}}{{\partial z}}} \right)\nonumber \\ &=& \frac{{{\partial ^2}}}{{\partial x\partial t}}\left( {\frac{{\partial {v_{xs}}}}{{\partial t}}} \right) + \frac{{{\partial ^2}}}{{\partial z\partial t}}\left( {\frac{{\partial {v_{zs}}}}{{\partial t}}} \right)\nonumber \\ &=& \frac{{{\partial ^2}}}{{\partial x\partial t}}\left( {\frac{1}{\rho }\frac{{\partial {\tau _{xxs}}}}{{\partial x}} + \frac{1}{\rho }\frac{{\partial {\tau _{xzs}}}}{{\partial z}}} \right) + \frac{{{\partial ^2}}}{{\partial z\partial t}}\left( {\frac{1}{\rho }\frac{{\partial {\tau _{xzs}}}}{{\partial x}} + \frac{1}{\rho }\frac{{\partial {\tau _{zzs}}}}{{\partial z}}} \right)\nonumber \\ &=& \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial {x^2}}}\left( {\frac{{\partial {\tau _{xxs}}}}{{\partial t}}} \right) + \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial x\partial z}}\left( {\frac{{\partial {\tau _{xzs}}}}{{\partial t}}} \right) + \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial z\partial x}}\left( {\frac{{\partial {\tau _{xzs}}}}{{\partial t}}} \right) + \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial {z^2}}}\left( {\frac{{\partial {\tau _{zzs}}}}{{\partial t}}} \right)\nonumber \\ &=& \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial {x^2}}}\left( { - 2\mu \frac{{\partial {v_z}}}{{\partial z}}} \right) + 2\frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial x\partial z}}\left( {\mu \left(\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}\right)} \right) + \frac{1}{\rho }\frac{{{\partial ^2}}}{{\partial {z^2}}}\left( { - 2\mu \frac{{\partial {v_x}}}{{\partial x}}} \right)\nonumber\\ &=& - 2\mu \frac{1}{\rho }\frac{{{\partial ^2}{v_z}}}{{\partial {x^2}\partial z}} + 2\mu \frac{1}{\rho }\frac{{{\partial ^2}{v_x}}}{{\partial x\partial {z^2}}} + 2\mu \frac{1}{\rho }\frac{{{\partial ^2}{v_z}}}{{\partial {x^2}\partial z}} - 2\mu \frac{1}{\rho }\frac{{{\partial ^2}{v_x}}}{{\partial x\partial {z^2}}}\nonumber \\ &=& 0. \end{eqnarray} (A4) From eq. (A4), we conclude that ∇ · vs is a constant or a linear function. According to wave characteristics, ∇ · vs cannot be a linear function, so ∇ · vs can only be a constant. In the initial conditions for solving wave equations, ∇ · vs|t = 0 = 0, so ∇ · vs ≡ 0. Similarly, we can prove that ∇ × vp ≡ 0. Therefore, vp in eq. (3b) is non-curl field, and vs in eq. (3c) is non-divergence field. Eqs (3) can produce pure P and S waves. APPENDIX B: GRADIENTS FOR LAMÉ PARAMETERS AND DENSITY PERTURBATIONS BASED ON COUPLED ELASTIC WAVE EQUATION Based on Born approximation theorem, parameter perturbations   \begin{equation}\delta {{\bf V}}(\lambda ,\mu ,\rho ) = {{\bf V}}(\lambda ,\mu ,\rho ) - {{{\bf V}}_0}(\lambda ,\mu ,\rho ),\end{equation} (B1)indicate wavefield perturbations   \begin{equation}\delta {{\bf U}} = {{\bf U}} - {{{\bf U}}_0},\end{equation} (B2)where U0 and V0 denote the background wavefields and parameters, and δU and δV denote the perturbation wavefields and parameters. U = (vx, vz, τxx, τzz, τxz) and V = (λ, μ, vp, vs, ρ). The background wavefield U0 can be identified by   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial {v_{x0}}}}{{\partial t}} = \frac{{\partial {\tau _{xx0}}}}{{\partial x}} + \frac{{\partial {\tau _{xz0}}}}{{\partial z}}\\ {\rho _0}\frac{{\partial {v_{z0}}}}{{\partial t}} = \frac{{\partial {\tau _{xz0}}}}{{\partial x}} + \frac{{\partial {\tau _{zz0}}}}{{\partial z}}\\ \frac{{\partial {\tau _{xx0}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\frac{{\partial {v_{x0}}}}{{\partial x}} + {\lambda _0}\frac{{\partial {v_{z0}}}}{{\partial z}} + {f_x}\\ \frac{{\partial {\tau _{zz0}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\frac{{\partial {v_{z0}}}}{{\partial z}} + {\lambda _0}\frac{{\partial {v_{x0}}}}{{\partial x}} + {f_z}\\ \frac{{\partial {\tau _{xz0}}}}{{\partial t}} = {\mu _0}\left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right..\end{equation} (B3) Substituting U(x, t) = U0(x, t) + δU(x, t) into eq. (1), subtracting eq. (B3) and neglecting the high-order terms of δV yields:   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_x}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xx}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{xz}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{x0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_z}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xz}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{zz}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{z0}}}}{{\partial t}}\\ \frac{{\partial \delta {\tau _{xx}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\frac{{\partial \delta {v_x}}}{{\partial x}} + {\lambda _0}\frac{{\partial \delta {v_z}}}{{\partial z}} + {{\bf m}}(\lambda ){\lambda _0}\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2{{\bf m}}(\mu ){\mu _0}\frac{{\partial {v_{x0}}}}{{\partial x}}\\ \frac{{\partial \delta {\tau _{zz}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\frac{{\partial \delta {v_z}}}{{\partial z}} + {\lambda _0}\frac{{\partial \delta {v_x}}}{{\partial x}} + {{\bf m}}(\lambda ){\lambda _0} \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2{{\bf m}}(\mu ){\mu _0}\frac{{\partial {v_{z0}}}}{{\partial z}}\\ \frac{{\partial \delta {\tau _{xz}}}}{{\partial t}} = {\mu _0}\left(\frac{{\partial \delta {v_x}}}{{\partial z}} + \frac{{\partial \delta {v_z}}}{{\partial x}}\right) + {{\bf m}}(\mu ){\mu _0}\left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right.,\end{equation} (B4)where the reflectivity models with respect to parameter perturbations are defined as   \begin{equation}\left\{ \begin{array}{@{}l@{}} {{\bf m}}(\lambda ) = \frac{{\delta \lambda }}{{{\lambda _0}}}\\ {{\bf m}}(\mu ) = \frac{{\delta \mu }}{{{\mu _0}}}\\ {{\bf m}}(\rho ) = \frac{{\delta \rho }}{{{\rho _0}}} \end{array} \right..\end{equation} (B5) The perturbation of the objective function δE(m) can be written by   \begin{eqnarray} \delta E({{\bf m}}) = \frac{1}{2}\delta \left\| {{{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}}} \right\|_2^2 &=& \frac{1}{2}\delta \left\langle {{{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}},{{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}}} \right\rangle = \left\langle {\delta ({{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}}),} \right.\nonumber \\ \left. {({{\bf Lm}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle &=& \left\langle {\delta ({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}}),({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle = \left\langle {{{\bf B}}\delta {{\bf U}},({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle, \end{eqnarray} (B6)where B denotes a sampling operator to restrict the data in the whole model space onto the receiver locations. From eq. (B4), we get   \begin{equation}\delta {{\bf U}} = {{\bf L}}{{{\bf F}}^\prime},\end{equation} (B7)where F΄ is a new constructed source matrix used to produce wavefield perturbation δU, and given by   \begin{equation}{{{\bf F}}^\prime} = \left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{x0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{z0}}}}{{\partial t}}\\ \delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\\ \delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\\ \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right).\end{equation} (B8) Thus, eq. (B6) can be changed to   \begin{equation}\delta E({{\bf m}}) = \left\langle {{{\bf B}}({{\bf L}}{{{\bf F}}^\prime}),({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle = \left\langle {{{{\bf F}}^\prime},{{{\bf L}}^R}{{{\bf B}}^R}({{\bf BU}} - {{{\bf d}}_{{\rm{obs}}}})} \right\rangle ,\end{equation} (B9)where superscript R denotes the adjoint operator. BR denotes where the data at the receivers extends to the whole model space, LR denotes the wavefield reverse-time backpropagation operator. So δE(λ, μ, ρ) can be written as   \begin{eqnarray} \delta E(\lambda ,\mu ,\rho ) &=& {\left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{x0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{z0}}}}{{\partial t}}\\ \delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\\ \delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\\ \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right)^T}\left( \begin{array}{@{}l@{}} v_x^R\\ v_z^R\\ \tau _{xx}^R\\ \tau _{zz}^R\\ \tau _{xz}^R \end{array} \right)\nonumber\\ &=& - \delta \rho \frac{{\partial {v_{x0}}}{{\partial t}}} v_x^R - \delta \rho \frac{{\partial {v_{z0}}}}{{\partial t}}v_z^R + \left[ {\delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}} \right]\tau _{xx}^R\nonumber\\ &&+\,\, \left[ {\delta \lambda \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) + 2\delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}} \right]\tau _{zz}^R + \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right)\tau _{xz}^R.\end{eqnarray} (B10) To minimize the objective function of eq. (1), we obtain the gradient directions of E(m) for different material parameters:   \begin{equation}{{{\bf g}}_\lambda } = \frac{{\partial {{\bf E}}}}{{\partial \lambda }} = \int{{\left( {\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\tau _{xx}^R + \left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\tau _{zz}^R} \right)}}dt,\end{equation} (B11)  \begin{equation}{{{\bf g}}_\mu } = \frac{{\partial {{\bf E}}}}{{\partial \mu }} = \int{{\left( {2\frac{{\partial {v_{x0}}}}{{\partial x}}\tau _{xx}^R + 2\frac{{\partial {v_{z0}}}}{{\partial z}}\tau _{zz}^R + \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right)\tau _{xz}^R} \right)}}dt,\end{equation} (B12)  \begin{equation}{{{\bf g}}_\rho } = \frac{{\partial {{\bf E}}}}{{\partial \rho }} = \int{{\left( { - \frac{{\partial {v_{x0}}}}{{\partial t}}v_x^R - \frac{{\partial {v_{z0}}}}{{\partial t}}v_z^R} \right)}}dt.\end{equation} (B13) APPENDIX C: GRADIENTS FOR ELSRTM BASED ON DECOUPLED ELASTIC WAVE EQUATION Rewriting eqs (3b) and (3c) using the relationship $${V_p} = \sqrt {(\lambda + 2\mu )/\rho } $$ and $${V_s} = \sqrt {\mu /\rho } $$ yields   \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_{xp}}}}{{\partial t}} = \frac{{\partial {\tau _{xxp}}}}{{\partial x}}\\ \rho \frac{{\partial {v_{zp}}}}{{\partial t}} = \frac{{\partial {\tau _{zzp}}}}{{\partial z}}\\ \frac{1}{{V_p^2}}\frac{{\partial {\tau _{xxp}}}}{{\partial t}} = \rho \left(\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}\right) + {f_x}\\ \frac{1}{{V_p^2}}\frac{{\partial {\tau _{zzp}}}}{{\partial t}} = \rho \left(\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}\right) + {f_z} \end{array} \right.,\end{equation} (C1)and   \begin{equation}\left\{ \begin{array}{@{}l@{}} \rho \frac{{\partial {v_{xs}}}}{{\partial t}} = \frac{{\partial {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial {\tau _{xzs}}}}{{\partial z}}\\ \rho \frac{{\partial {v_{zs}}}}{{\partial t}} = \frac{{\partial {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial {\tau _{zzs}}}}{{\partial z}}\\ \frac{1}{{V_s^2}}\frac{{\partial {\tau _{xxs}}}}{{\partial t}} = - 2\rho \frac{{\partial {v_z}}}{{\partial z}} + {f_x}\\ \frac{1}{{V_s^2}}\frac{{\partial {\tau _{zzs}}}}{{\partial t}} = - 2\rho \frac{{\partial {v_x}}}{{\partial x}} + {f_z}\\ \frac{1}{{V_s^2}}\frac{{\partial {\tau _{xzs}}}}{{\partial t}} = \rho \left(\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}\right) \end{array} \right..\end{equation} (C2) Substituting U0 (Up0, Us0) and V0 (Vp0, Vs0, ρ) into eqs (C1) and (C2) to obtain the P- and S-wave background equations. Expanding the Vp and Vs terms according to   \begin{equation}\frac{1}{{V_p^2}} \approx \frac{1}{{V_{p0}^2}} - \frac{{2\delta {V_p}}}{{V_{p0}^3}},\end{equation} (C3)  \begin{equation}\frac{1}{{V_s^2}} \approx \frac{1}{{V_{s0}^2}} - \frac{{2\delta {V_s}}}{{V_{s0}^3}},\end{equation} (C4)and subtracting the P- and S-wave background equations from eqs (C1) and (C2), respectively, yields the equations for δUp and δUs after neglecting the high-order terms yields:   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxp}}}}{{\partial x}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{xp0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{zzp}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{zp0}}}}{{\partial t}}\\ \frac{1}{{V_{p0}^2}}\frac{{\partial \delta {\tau _{xxp}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + \frac{{2{{\bf m}}({V_p})}}{{V_{p0}^2}}\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\\ \frac{1}{{V_{p0}^2}}\frac{{\partial \delta {\tau _{zzp}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + \frac{{2{{\bf m}}({V_p})}}{{V_{p0}^2}}\frac{{\partial {\tau _{zzp0}}}}{{\partial t}} \end{array} \right.,\end{equation} (C5)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{xzs}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{xs0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{zzs}}}}{{\partial z}} - {{\bf m}}(\rho ){\rho _0}\frac{{\partial {v_{zs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{xxs}}}}{{\partial t}} = - 2{\rho _0}\frac{{\partial \delta {v_z}}}{{\partial z}} + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{zzs}}}}{{\partial t}} = - 2{\rho _0}\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{zzs0}}}}{{\partial t}}\\ \frac{1}{{V_{s0}^2}}\frac{{\partial \delta {\tau _{xzs}}}}{{\partial t}} = {\rho _0}\left(\frac{{\partial \delta {v_x}}}{{\partial z}} + \frac{{\partial \delta {v_z}}}{{\partial x}}\right) + \frac{{2{{\bf m}}({V_s})}}{{V_{s0}^2}}\frac{{\partial {\tau _{xzs0}}}}{{\partial t}} \end{array} \right.,\end{equation} (C6)where the reflectivity models associated with velocities and density perturbations are defined as   \begin{equation}{{\bf m}}({V_p}) = \frac{{\delta {V_p}}}{{{V_{p0}}}},\end{equation} (C7)  \begin{equation}{{\bf m}}({V_s}) = \frac{{\delta {V_s}}}{{{V_{s0}}}},\end{equation} (C8)  \begin{equation}{{\bf m}}(\rho ) = \frac{{\delta \rho }}{\rho }.\end{equation} (C9) For velocities and density perturbations, eq. (4) can be expressed as   \begin{eqnarray} E({V_p},{V_s},\rho ) &=& \left( {{{\left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{xp0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{zp0}}}}{{\partial t}}\\ \frac{{2\delta {V_p}}}{{V_{p0}^3}}\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\\ \frac{{2\delta {V_p}}}{{V_{p0}^3}}\frac{{\partial {\tau _{zzp0}}}}{{\partial t}} \end{array} \right)}^T}\left( \begin{array}{@{}l@{}} v_{xp}^R\\ v_{zp}^R\\ \tau _{xxp}^R\\ \tau _{zzp}^R \end{array} \right) + {{\left( {\begin{array}{@{}*{1}{c}@{}} { - \delta \rho \frac{{\partial {v_{xs0}}}}{{\partial t}}}\\ { - \delta \rho \frac{{\partial {v_{zs0}}}}{{\partial t}}}\\ {\frac{{2\delta {V_s}}}{{V_{s0}^3}}\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}}\\ {\frac{{2\delta {V_s}}}{{V_{s0}^3}}\frac{{\partial {\tau _{zzs0}}}}{{\partial t}}}\\ {\frac{{2\delta {V_s}}}{{V_{s0}^3}}\frac{{\partial {\tau _{xzs0}}}}{{\partial t}}} \end{array}} \right)}^T}\left( \begin{array}{@{}l@{}} v_{xs}^R\\ v_{zs}^R\\ \tau _{xxs}^R\\ \tau _{zzs}^R\\ \tau _{xzs}^R \end{array} \right)} \right)\nonumber\\ &=& \frac{{\delta {V_p}}}{{V_{p0}^3}}\left( {\frac{{\partial {\tau _{xxp0}}}}{{\partial t}}\tau _{xxp}^R + \frac{{\partial {\tau _{zzp0}}}}{{\partial t}}\tau _{zzp}^R} \right) + \frac{{\delta {V_s}}}{{V_{s0}^3}}\left( {\frac{{\partial {\tau _{xxs0}}}}{{\partial t}}\tau _{xxs}^R + \frac{{\partial {\tau _{zzs0}}}}{{\partial t}}\tau _{zzs}^R}+\, \frac{{\partial {\tau _{xzs0}}}}{{\partial t}}\tau _{xzs}^R \right) \nonumber\\ &&-\,\, \delta \rho \left( {\frac{{\partial {v_{xp0}}}}{{\partial t}}v_{xp}^R + \frac{{\partial {v_{zp0}}}}{{\partial t}}v_{zp}^R + \frac{{\partial {v_{xs0}}}}{{\partial t}}v_{xs}^R + \frac{{\partial {v_{zs0}}}}{{\partial t}}v_{zs}^R} \right).\end{eqnarray} (C10) The gradient formulae with respect toVp, Vs and ρ can be written as eq. (10) by calculating the partial derivative of Vp, Vs and ρ, respectively. Similarly, when performing ELSRTM for Lamé parameters and density perturbations based on decoupled elastic wave equation, the equations for δUp and δUs are given by   \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxp}}}}{{\partial x}} - \delta \rho \frac{{\partial {v_{xp0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zp}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{zzp}}}}{{\partial z}} - \delta \rho \frac{{\partial {v_{zp0}}}}{{\partial t}}\\ \frac{{\partial \delta {\tau _{xxp}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\\ \frac{{\partial \delta {\tau _{zzp}}}}{{\partial t}} = ({\lambda _0} + 2{\mu _0})\left(\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_z}}}{{\partial z}}\right) + (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) \end{array} \right.,\end{equation} (C11)  \begin{equation}\left\{ \begin{array}{@{}l@{}} {\rho _0}\frac{{\partial \delta {v_{xs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xxs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{xzs}}}}{{\partial z}} - \delta \rho \frac{{\partial {v_{xs0}}}}{{\partial t}}\\ {\rho _0}\frac{{\partial \delta {v_{zs}}}}{{\partial t}} = \frac{{\partial \delta {\tau _{xzs}}}}{{\partial x}} + \frac{{\partial \delta {\tau _{zzs}}}}{{\partial z}} - \delta \rho \frac{{\partial {v_{zs0}}}}{{\partial t}}\\ \frac{{\partial \delta {\tau _{xxs}}}}{{\partial t}} = - 2{\mu _0}\frac{{\partial \delta {v_z}}}{{\partial z}} - 2\delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\\ \frac{{\partial \delta {\tau _{zzs}}}}{{\partial t}} = - 2{\mu _0}\frac{{\partial \delta {v_x}}}{{\partial x}} - 2\delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\\ \frac{{\partial \delta {\tau _{xzs}}}}{{\partial t}} = {\mu _0}\left(\frac{{\partial \delta {v_x}}}{{\partial z}} + \frac{{\partial \delta {v_z}}}{{\partial x}}\right) + \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right..\end{equation} (C12) For velocities and density perturbations, eq. (4) can be written as   \begin{eqnarray} {{\bf E}}(\lambda ,\mu ,\rho ) &=& {\left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{xp0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{zp0}}}}{{\partial t}}\\ (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\\ (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right) \end{array} \right)^T}\left( \begin{array}{@{}l@{}} v_{xp}^R\\ v_{zp}^R\\ \tau _{xxp}^R\\ \tau _{zzp}^R \end{array} \right) + {\left( \begin{array}{@{}l@{}} - \delta \rho \frac{{\partial {v_{xs0}}}}{{\partial t}}\\ - \delta \rho \frac{{\partial {v_{zs0}}}}{{\partial t}}\\ - \delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\\ - \delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\\ \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right) \end{array} \right)^T}\left( \begin{array}{@{}l@{}} v_{xs}^R\\ v_{zs}^R\\ \tau _{xxs}^R\\ \tau _{zzs}^R\\ \tau _{xzs}^R \end{array} \right)\nonumber\\ &=& - \delta \rho \frac{{\partial {v_{xp0}}}}{{\partial t}} v_{xp}^R - \delta \rho \frac{{\partial {v_{zp0}}}}{{\partial t}}v_{zp}^R + (\delta \lambda + 2\delta \mu )\left(\frac{{\partial {v_{x0}}}}{{\partial x}} + \frac{{\partial {v_{z0}}}}{{\partial z}}\right)\left(\tau _{xxp}^R + \tau _{zzp}^R\right)\nonumber\\ &&-\,\,\delta \rho \frac{{\partial {v_{xs0}}}}{{\partial t}}v_{xs}^R - \delta \rho \frac{{\partial {v_{zs0}}}}{{\partial t}}v_{zs}^R - \delta \mu \frac{{\partial {v_{z0}}}}{{\partial z}}\tau _{xxs}^R - \delta \mu \frac{{\partial {v_{x0}}}}{{\partial x}}\tau _{zzs}^R + \delta \mu \left(\frac{{\partial {v_{x0}}}}{{\partial z}} + \frac{{\partial {v_{z0}}}}{{\partial x}}\right)\tau _{xzs}^R. \end{eqnarray} (C13) Therefore, the gradients with respect to the Lamé parameters and density can be accessed as eq. (14) by calculating the partial derivative of λ, μ and ρ, respectively. © Crown copyright 2017. This article contains public sector information licensed under the Open Government Licence v3.0 (http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/).

Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off