TY - JOUR
AU - Wu,, Di
AB - Abstract Acoustic impedance estimation has a significant importance to seismic exploration. In this paper, we use full-waveform inversion to recover the impedance from seismic data, and analyze the sensitivity of the acoustic impedance with respect to the source-receiver offset of seismic data and to the initial velocity model. We parameterize the acoustic wave equation with velocity and impedance, and demonstrate three key aspects of acoustic impedance inversion. First, short-offset data are most suitable for acoustic impedance inversion. Second, acoustic impedance inversion is more compatible with the data generated by density contrasts than velocity contrasts. Finally, acoustic impedance inversion requires the starting velocity model to be very accurate for achieving a high-quality inversion. Based upon these observations, we propose a workflow for acoustic impedance inversion as: (1) building a background velocity model with travel-time tomography or reflection waveform inversion; (2) recovering the intermediate wavelength components of the velocity model with full-waveform inversion constrained by Gardner’s relation; (3) inverting the high-resolution acoustic impedance model with short-offset data through full-waveform inversion. We verify this workflow by the synthetic tests based on the Marmousi model. acoustic impedance, full-waveform inversion, model parameterization, short offset reflection Introduction Seismic exploration aims to recover seismic properties of the subsurface. Generally, there are two ways to obtain these seismic properties: seismic imaging and seismic inversion. The seismic imaging is achieved through migration, which moves reflection events in the records to the location where the reflection happens and stacks them to form the migration images (Yilmaz 2001). Usually, migration images only show the structure of the subsurface and they are mainly used for geological structural interpretation. As the adjoint operators of the forward modeling operators are used instead of their inverse operators, the mainstream methods for migration are not quantitative procedures (Claerbout 1992) except formed in an inversion framework, e.g. least-squares migration (Nemeth et al1999, Kühl and Sacchi 2003, Liu et al2005, Kaplan et al2010, Dai et al2012, Yao and Wu 2015, Zhang et al2015, Wu et al2016, Yao and Jakubowicz 2016). Seismic inversion, however, extracts the quantitative parameters of subsurface rocks, such as P- and S-wave velocity and impedance. Currently, surface acquisition dominates seismic exploration. With this type of acquisition geometry, reflections are the main source carrying subsurface information. Reflections are generated by impedance contrasts between different materials. Consequently, impedance is a key property in seismic exploration for characterizing the geological formations. Wang (2016, chapter 14) proposed a waveform tomography method for seismic impedance, in which the acoustic wave equation is defined in terms of acoustic impedance rather than the velocity model. There are three categories of inversion methods for estimating impedance. The first one is based on the fact that the amplitude of reflections changes with offsets or reflection angles. This phenomenon is called amplitude versus offset or angle, short for AVO or AVA (Koefoed 1962, Tooley et al1965, Ostrander 1984, Zhang and Li 2015, 2016). The reflectivity for different combinations of rocks can be calculated by using the Zoeppritz equations. Thus, P-wave velocity and S-wave velocity and density can be inverted by using the Zoeppritz equations. This category of inversion is then referred to as AVO inversion. Since the Zoeppritz equations are nonlinear and complicated, geophysicists have made efforts to simplify them with approximations (e.g., Shuey 1985, Fatti et al1994, Aki and Richards 1997, Wang 1999). After velocities and density are recovered, the impedance can be computed easily by multiplying the velocity and density. The second category of methods uses the technique of post-stack inversion to convert the migration image into the impedance profile. The principle of this technique is very simple. Migration images represent reflectivity (Yao and Wu 2015, Zhang et al2015) or impedance perturbation (Dai et al2012). The former is a derivative of the latter (Yao et al2017a). Thus, the post-stack inversion can be done easily. First, computing the impedance log according to the well logs of velocity and density. Second, inverting a matched filter, which can convert the migration image near the well to the impedance log by convolution. Finally, the impedance profile can be calculated by convolving the matched filter with the seismic image (Lancaster and Whitcombe 2000, Latimer et al2000). The final step is based on the following observation: the spectrum of the Earth’s reflection coefficient series with respect to location in a field varies insignificantly (Velzeboer 1981, Walden and Hosken 1985). This post-stack inversion technique is simple and robust, so it is widely used. The third category of methods is impedance inversion using full waveform inversion (FWI). FWI is a geophysical inversion method which aims to find an optimal seismic model to fit the seismic data (Tarantola 1984, Pratt and Shipp 1999, Virieux and Operto 2009, Wang and Rao 2009, da Silva and Yao 2018). If the acoustic wave equation used in FWI is parameterized as velocity and impedance, and the velocity is fixed during the inversion, then acoustic impedance inversion with FWI can be achieved. Since impedance contrast between different layers generates the reflections, reflections are the main data in the impedance inversion. One of the main applications of impedance inversion is to produce high-resolution impedance models, which can be used to interpret lithology directly (Routh et al2016), identify and characterize hydrocarbon reservoirs (Lu 2016) and shallow gas hazards (Bright et al2015). Usually, the frequencies used in the inversion are pushed up to 40 Hz and higher to achieve high-resolution images. Consequently, the computational cost is much higher than FWI for inverting velocity, which normally uses frequencies lower than 10 Hz. Inversion at such high frequencies requires efficient finite-difference algorithms (Zhang and Yao 2012, Liu 2013, Wang et al2014, Yao et al2016, Yao et al2017b); in addition, some other inversion strategies, e.g. spectral shaping (Lazaratos et al2011, Plessix and Li 2013), are used to improve the convergence rate. In this paper, we first review the algorithm for acoustic impedance inversion with FWI. The sensitivity analysis of the impedance inversion will then be carried out. The analysis demonstrates that (1) short-offset data are the most suitable input, (2) the impedance inversion algorithm is more compatible with data generated by density contrasts than velocity contrasts, and (3) the inversion needs an accurate starting velocity model. The numerical tests on the Marmousi model validate these analyses. Finally, we propose a cascading inversion workflow to achieve a high-quality impedance inversion. This inversion workflow is demonstrated on the Marmousi model. Impedance inversion with full-waveform inversion The acoustic wave equation reads |$ begin{eqnarray}\left(\displaystyle \frac{1}{\rho ({\bf{x}}){v}^{2}({\bf{x}})}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{\rm{\nabla }}\cdot \displaystyle \frac{1}{\rho ({\bf{x}})}{\rm{\nabla }}\right)p({\bf{x}},t)=f({\bf{x}},t),\end{equation} $| 1ρ(x)v2(x)∂2∂t2−∇⋅1ρ(x)∇p(x,t)=f(x,t), 1 where p(x, t) is the pressure wavefield, v(x) is the acoustic wave speed, ρ(x) is the density, f(x, t) is the source term, x is the vector of the spatial coordinates, and t indicates time. Seismic exploration with active sources, usually, applies a point source near the surface in each shot, so the source term is commonly expressed as \begin{equation}f({\bf{x}},t)=s(t)\delta ({\bf{x}}-{{\bf{x}}}_{s}),\end{equation} f(x,t)=s(t)δ(x−xs), 2 where s(t) is the source wavelet at the position xs and δ is the Dirac delta function. The initial conditions are normally set as begin{eqnarray}\left\{\begin{array}{c}p({\bf{x}},t=0)=0,\\ \displaystyle \frac{\partial }{\partial t}p({\bf{x}},t=0)=0.\end{array}\right.\end{equation} p(x,t=0)=0,∂∂tp(x,t=0)=0. 3 As a result, all the energy of the pressure wavefield, p, comes from the source term. Except the initial conditions, wave equation (1) needs boundary conditions to ensure a unique solution. Usually, the boundary condition is |$ begin{eqnarray}p({\rm{x}}\in \partial {\rm{\Omega }},t)=0,\end{equation} $| p(x∈∂Ω,t)=0, 4 where ∂Ω denotes the boundary of the domain. In seismic exploration, usually, the top boundary is a free surface, so that equation (4) is explicitly satisfied. In addition, imposing absorbing lateral and bottom boundaries also guarantees that the boundary condition is met. Hence, equation (4) is satisfied along the entire boundary. If we parameterize equation (1) with impedance, |$ begin{eqnarray}I=\rho v,\end{equation} $| I=ρv, 5 and velocity, then equation (1) becomes |$ begin{eqnarray}\left(\displaystyle \frac{1}{I({\bf{x}})v({\bf{x}})}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{I({\bf{x}})}{\rm{\nabla }}\right)p({\bf{x}},t)=f({\bf{x}},t),\end{equation} $| 1I(x)v(x)∂2∂t2−∇⋅v(x)I(x)∇p(x,t)=f(x,t), 6 where the pressure wavefield p(x, t) is a function of both I and v. The wave equation (6) represents a linear relationship between p and s. After discretization in both time and space with numerical methods, e.g. finite difference, equation (6) therefore can be converted into a matrix form as |$ begin{eqnarray}{\bf{A}}{\bf{p}}={\bf{f}},\end{equation} $| Ap=f, 7 where f and p are column vectors that represent the source and wavefield at discrete points in space and time, A is a matrix that represents the discrete numerical implementation of the operator |$ begin{eqnarray}\displaystyle \frac{1}{I({\bf{x}})v({\bf{x}})}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{I({\bf{x}})}{\rm{\nabla }}.\end{equation} $| 1I(x)v(x)∂2∂t2−∇⋅v(x)I(x)∇. 8 More details about the matrix formulation of equation (7) can be found in appendix A. Differentiating equation (7) with respect to impedance gives |$ begin{eqnarray}\displaystyle \frac{\partial {\bf{p}}}{\partial I({\bf{x}})}=-{{\bf{A}}}^{-1}\displaystyle \frac{\partial {\bf{A}}}{\partial I({\bf{x}})}{\bf{p}}.\end{equation} $| ∂p∂I(x)=−A−1∂A∂I(x)p. 9 The objective function for FWI is usually defined as the square of ℓ2-norm of the data misfit, which is |$ begin{eqnarray}\phi =\displaystyle \frac{1}{2}{\parallel {\bf{d}}\,-\,{{\bf{d}}}_{0}\parallel }_{2}^{2}=\displaystyle \frac{1}{2}{({\bf{d}}-{{\bf{d}}}_{0})}^{{\rm{T}}}({\bf{d}}\,-\,{{\bf{d}}}_{0})=\displaystyle \frac{1}{2}\delta {{\bf{d}}}^{{\rm{T}}}\delta {\bf{d}},\end{equation} $| ϕ=12∥d-d0∥22=12(d-d0)T(d-d0)=12δdTδd, 10 where T denotes transpose, d0 is the column vector of the observed data at discrete points in space and time, d is the column vector of the predicted data at the same position of d0, and δd is the residual. The predicted data d is a subset of the full predicted data p in equation (7) and can be mathematically expressed as |$ begin{eqnarray}{\bf{d}}={\bf{R}}{\bf{p}},\end{equation} $| d=Rp, 11 where R is a matrix that extracts the wavefields at the receiver position. The number of columns of R is the same as the number of elements of p while the number of rows of R is identical to the number of elements of d0. Thus, R is a short matrix usually and contains a very small number of non-zero element. R is referred to as the restriction matrix. By differentiating equation (10) with respect to I and using equations (9) and (11) and properties of transpose, we obtain the gradient of the objective function: |$ begin{eqnarray}{g}_{I}({\bf{x}}) & = & \displaystyle \frac{\partial \phi }{\partial I({\bf{x}})}={\left(\displaystyle \frac{\partial {\bf{d}}}{\partial I({\bf{x}})}\right)}^{{\rm{T}}}\delta {\bf{d}}\\ & = & {\left(\displaystyle \frac{\partial {\bf{R}}{\bf{p}}}{\partial I({\bf{x}})}\right)}^{{\rm{T}}}\delta {\bf{d}}\\ & = & {\left({\bf{R}}\displaystyle \frac{\partial {\bf{p}}}{\partial I({\bf{x}})}\right)}^{{\rm{T}}}\delta {\bf{d}}\\ & = & -{{\bf{p}}}^{{\rm{T}}}{\left(\displaystyle \frac{\partial {\bf{A}}}{\partial I({\bf{x}})}\right)}^{{\rm{T}}}{{\bf{A}}}^{-{\rm{T}}}{{\bf{R}}}^{{\rm{T}}}\delta {\bf{d}}\\ & = & -{{\bf{p}}}^{{\rm{T}}}{\left(-\displaystyle \frac{1}{{I}^{2}({\bf{x}})v({\bf{x}})}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}+{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{{I}^{2}({\bf{x}})}{\rm{\nabla }}\right)}^{{\rm{T}}}\\ & & \times {{\bf{A}}}^{-{\rm{T}}}{{\bf{R}}}^{{\rm{T}}}\delta {\bf{d}}\\ & = & \displaystyle \frac{1}{{I}^{2}({\bf{x}})v({\bf{x}})}{\left(\displaystyle \frac{{\partial }^{2}{\bf{p}}}{\partial {t}^{2}}\right)}^{{\rm{T}}}{{\bf{A}}}^{-{\rm{T}}}{{\bf{R}}}^{{\rm{T}}}\delta {\bf{d}}\\ & & +\displaystyle \frac{v({\bf{x}})}{{I}^{2}({\bf{x}})}{\rm{\nabla }}{\bf{p}}\cdot {\rm{\nabla }}({{\bf{A}}}^{-{\rm{T}}}{{\bf{R}}}^{{\rm{T}}}\delta {\bf{d}}),\end{equation} $| gI(x)=∂ϕ∂I(x)=∂d∂I(x)Tδd=∂Rp∂I(x)Tδd=R∂p∂I(x)Tδd=−pT∂A∂I(x)TA−TRTδd=−pT−1I2(x)v(x)∂2∂t2+∇⋅v(x)I2(x)∇T×A−TRTδd=1I2(x)v(x)∂2p∂t2TA−TRTδd+v(x)I2(x)∇p⋅∇(A−TRTδd), 12 where RT is the transpose of the restriction matrix and therefore a tall-skinny matrix, RTδd means to extend δd to the same length as p, and the adjoint operator A−T means (A−1)T. The transpose of A−1 reverses its time-propagation direction, so that A−TRTδd reverse-time propagates the residual into the Earth to generate the residual wavefield (more details about the adjoint operator A−T can be found in appendix A. Equation (12) shows that the gradient of impedance includes two terms: the first term is formed by zero-lag cross-correlating the second time derivative of the source wavefield and the residual wavefield, and then scaling the cross-correlation by a factor of | $\displaystyle \frac{1}{{I}^{2}({\bf{x}})v({\bf{x}})};$ | 1I2(x)v(x); the second term is generated by zero-lag cross-correlating the spatial derivatives of the source wavefield and the residual wavefield, and then scaling the cross-correlation by a factor of | $\displaystyle \frac{v({\bf{x}})}{{I}^{2}({\bf{x}})}.$ | v(x)I2(x). Similarly, the gradient of velocity can be derived as |$ begin{eqnarray}{g}_{v}({\bf{x}}) & = & \displaystyle \frac{\partial \phi }{\partial v({\bf{x}})}={\left(\displaystyle \frac{\partial {\bf{d}}}{\partial v({\bf{x}})}\right)}^{{\rm{T}}}\delta {\bf{d}}\\ & = & -{{\bf{p}}}^{{\rm{T}}}{\left(\displaystyle \frac{\partial {\bf{A}}}{\partial v({\bf{x}})}\right)}^{{\rm{T}}}{{\bf{A}}}^{-{\rm{T}}}{{\bf{R}}}^{{\rm{T}}}\delta {\bf{d}}\\ & = & -{{\bf{p}}}^{{\rm{T}}}{\left(-\displaystyle \frac{1}{I({\bf{x}}){v}^{2}({\bf{x}})}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{\rm{\nabla }}\cdot \displaystyle \frac{1}{I({\bf{x}})}{\rm{\nabla }}\right)}^{{\rm{T}}}\\ & & \,\,\times {{\bf{A}}}^{-{\rm{T}}}{{\bf{R}}}^{{\rm{T}}}\delta {\bf{d}}\\ & = & \displaystyle \frac{1}{I({\bf{x}}){v}^{2}({\bf{x}})}{\left(\displaystyle \frac{{\partial }^{2}{\bf{p}}}{\partial {t}^{2}}\right)}^{{\rm{T}}}{{\bf{A}}}^{-{\rm{T}}}{{\bf{R}}}^{{\rm{T}}}\delta {\bf{d}}\\ & & -\displaystyle \frac{1}{I({\bf{x}})}{\rm{\nabla }}{\bf{p}}\cdot {\rm{\nabla }}({{\bf{A}}}^{-{\rm{T}}}{{\bf{R}}}^{{\rm{T}}}\delta {\bf{d}}).\end{equation} $| gv(x)=∂ϕ∂v(x)=∂d∂v(x)Tδd=−pT∂A∂v(x)TA−TRTδd=−pT−1I(x)v2(x)∂2∂t2−∇⋅1I(x)∇T×A−TRTδd=1I(x)v2(x)∂2p∂t2TA−TRTδd−1I(x)∇p⋅∇(A−TRTδd). 13 According to the impedance gradient of the objective function, i.e. equation (12), the acoustic impedance inversion can be achieved by applying gradient-based methods (Wang 2016), e.g. steep-descent or conjugate-gradient, to minimize the objective function (10). Sensitivity analyses of the acoustic impedance inversion Herein we carry out the sensitivity analysis by deriving radiation patterns for impedance, velocity and density with the two parameterizations, i.e. velocity-impedance and velocity-density. The rationale is based on deriving the response of the gradient of the objective function for a diffracting anomaly over a homogeneous background. This gives the diffraction response, and allows to determine how the sensitivity to model parameters changes with the source-receiver opening angle. We outline this analysis in the frequency domain as it simplifies the derivation of the radiation patterns with Green’s functions. By using Green’s function, the source wavefield ps generated by a source s at xs is computed by solving equation (7), in the frequency domain can be expressed as |$ begin{eqnarray}{p}_{s}({\bf{x}},\omega )={G}_{s}({\bf{x}}| {{\bf{x}}}_{s},\omega )s(\omega ),\end{equation} $| ps(x,ω)=Gs(x∣xs,ω)s(ω), 14 where Gs(x∣xs, ω) is the source Green’s function from the source position xs to x, and s(ω) is the source wavelet in the frequency domain. The source wavefield ps is the same as p in equations (12) and (13). Similarly, the residual wavefield, which is computed by backward propagating of the residual, i.e. A−TRTδd, can be expressed as |$ begin{eqnarray}{p}_{r}({\bf{x}},\omega )={G}_{r}^{* }({\bf{x}}| {{\bf{x}}}_{r},\omega )\delta d({{\bf{x}}}_{r},\omega ),\end{equation} $| pr(x,ω)=Gr*(x∣xr,ω)δd(xr,ω), 15 where Gr(x∣xr, ω) is the receiver Green’s function from the receiver position xr to x, δd(xr, ω) is the data residual in the frequency domain, and * denotes the complex conjugate. In the frequency domain, the complex conjugate of Green’s function Gr acting on the data residual δd is equivalent to time-reversal propagation in the time domain. The residual wavefield pr is generated by the data residual at the single receiver position xr. This can simplify the following derivation. For more receiver positions, there should be an integral over xr in equation (15). With the asymptotic approximation (Červený 2001), the source Green’s function can be expressed as |$ begin{eqnarray}{G}_{s}\approx {A}_{s}{{\rm{e}}}^{{\rm{j}}\omega {\tau }_{s}({\bf{x}})},\end{equation} $| Gs≈Asejωτs(x), 16 and the receiver Green’s function can be expressed As and Ar are the amplitude of the source and residual Green’s functions, respectively, τs and τr are the travel-time from the source and from the receiver to a given point in the model space, respectively. By applying Fourier transform on the gradient of impedance (equation (12)) and inserting equations (15)–(17) into it, one can obtain |$ begin{eqnarray}{g}_{I}({\bf{x}}) & = & \displaystyle \sum _{\omega }\left(-\displaystyle \frac{{\omega }^{2}}{{I}^{2}v}{p}_{s}^{\ast }{p}_{r}+\displaystyle \frac{v}{{I}^{2}}{\rm{\nabla }}{p}_{s}^{\ast }\cdot {\rm{\nabla }}{p}_{r}\right)\\ & = & \displaystyle \sum _{\omega }\left(-\displaystyle \frac{{\omega }^{2}}{{I}^{2}v}{p}_{s}^{\ast }{p}_{r}-\displaystyle \frac{{\omega }^{2}v}{{I}^{2}}{A}_{s}{{\rm{e}}}^{-{\rm{j}}\omega {\tau }_{s}}s{A}_{r}{{\rm{e}}}^{-{\rm{j}}\omega {\tau }_{r}}\right.\\ & & \,\times \delta d{\rm{\nabla }}{\tau }_{s}\cdot {\rm{\nabla }}{\tau }_{r})\\ & = & \displaystyle \sum _{\omega }\left(-\displaystyle \frac{{\omega }^{2}}{{I}^{2}v}{p}_{s}^{\ast }{p}_{r}-\displaystyle \frac{{\omega }^{2}v}{{I}^{2}}{p}_{s}^{\ast }{p}_{r}{\rm{\nabla }}{\tau }_{s}\cdot {\rm{\nabla }}{\tau }_{r}\right),\end{equation} $| gI(x)=∑ω−ω2I2vps⁎pr+vI2∇ps⁎⋅∇pr=∑ω−ω2I2vps⁎pr−ω2vI2Ase−jωτssAre−jωτr×δd∇τs⋅∇τr)=∑ω−ω2I2vps⁎pr−ω2vI2ps⁎pr∇τs⋅∇τr, 18 where * denotes the complex conjugate. For isotropic media, the dot product of the gradient of the travel-times is given by |$ begin{eqnarray}{\rm{\nabla }}{\tau }_{s}({\bf{x}})\cdot {\rm{\nabla }}{\tau }_{r}({\bf{x}})=\displaystyle \frac{\cos \,\theta }{{v}^{2}({\bf{x}})},\end{equation} $| ∇τs(x)⋅∇τr(x)=cosθv2(x), 19 where θ is the opening angle between the source position xs and the receiver position xr to the model element position x (Bleistein et al2005, Virieux and Operto 2009). Plugging equation (19) into equation (18) yields |$ begin{eqnarray}{g}_{I}({\bf{x}},\theta ) & = & \displaystyle \sum _{\omega }\left(-\displaystyle \frac{{\omega }^{2}}{{I}^{2}v}{p}_{s}^{\ast }{p}_{r}-\displaystyle \frac{{\omega }^{2}\,\cos \,\theta }{{I}^{2}v}{p}_{s}^{\ast }{p}_{r}\right)\\ & = & \displaystyle \sum _{\omega }\left(-\displaystyle \frac{{\omega }^{2}(1+\,\cos \,\theta )}{{I}^{2}({\bf{x}})v({\bf{x}})}{p}_{s}^{\ast }({\bf{x}},\omega ){p}_{r}({\bf{x}},\omega )\right).\end{equation} $| gI(x,θ)=∑ω−ω2I2vps⁎pr−ω2cosθI2vps⁎pr=∑ω−ω2(1+cosθ)I2(x)v(x)ps⁎(x,ω)pr(x,ω). 20 Perturbing the impedance model at x, applying the Born approximation to equation (6) and combining the asymptotic analysis above, we can obtain the diffractor response to the impedance perturbation against the opening angle θ as |$ begin{eqnarray}\delta p({{\bf{x}}}_{r},\omega ,\theta ) & = & -\displaystyle \frac{{\omega }^{2}(1+\,\cos \,\theta )}{{I}^{2}({\bf{x}})v({\bf{x}})}G({{\bf{x}}}_{r}| {\bf{x}},\omega )\\ & & \times {p}_{s}({\bf{x}},\omega )\delta I({\bf{x}}),\end{equation} $| δp(xr,ω,θ)=−ω2(1+cosθ)I2(x)v(x)G(xr∣x,ω)×ps(x,ω)δI(x), 21 where θ is the opening angle between the source position xs and the receiver position xr to the model element position x, G(xr∣x, ω) is the Green’s function from x to xr, δI is the impedance perturbation, and δp is the diffractor response of δI. In appendix B we outline a detailed derivation of equation (21). Equation (21) shows how the diffraction changes against the diffraction angle when perturbing the impedance but velocity is fixed. In equations (20) and (21), | $\displaystyle \frac{1+\,\cos \,\theta }{2}$ | 1+cosθ2 shows the normalized radiation pattern, which is also called diffraction pattern (e.g. Zhou et al2015, da Silva et al2016). By comparison, one can see that the gradient shows the same radiation pattern as the diffraction response. Similarly to equation (20), we can derive the gradient of velocity under the asymptotic approximation from equation (13) as |$ begin{eqnarray}{g}_{v}\,({\bf{x}},\theta ) & = & \displaystyle \sum _{\omega }\left(-\displaystyle \frac{{\omega }^{2}}{I{v}^{2}}{p}_{s}^{\ast }{p}_{r}+\displaystyle \frac{{\omega }^{2}\,\cos \,\theta }{I{v}^{2}}{p}_{s}^{\ast }{p}_{r}\right)\\ & = & \displaystyle \sum _{\omega }\left(-\displaystyle \frac{{\omega }^{2}(1-\,\cos \,\theta )}{I({\bf{x}}){v}^{2}({\bf{x}})}{p}_{s}^{\ast }({\bf{x}},\omega ){p}_{r}({\bf{x}},\omega )\right),\end{equation} $| gv(x,θ)=∑ω−ω2Iv2ps⁎pr+ω2cosθIv2ps⁎pr=∑ω−ω2(1−cosθ)I(x)v2(x)ps⁎(x,ω)pr(x,ω), 22 and the diffractor response as |$ begin{eqnarray}\delta p({\bf{x}},\omega ,\theta ) & = & -\displaystyle \frac{{\omega }^{2}(1-\,\cos \,\theta )}{I({\bf{x}}){v}^{2}({\bf{x}})}G({{\bf{x}}}_{r}| {\bf{x}},\omega )\\ & & \times {p}_{s}({\bf{x}},\omega )\delta v({\bf{x}}).\end{equation} $| δp(x,ω,θ)=−ω2(1−cosθ)I(x)v2(x)G(xr∣x,ω)×ps(x,ω)δv(x). 23 Equations (20) and (21) show that the data have a decreasing sensitivity to the impedance with an increase of the opening angle. When θ is 180°, the gradient vanishes and the diffractor response is null. This means that if two events, one of which is from the source wavefield and the other of which is from the residual wavefield, propagate opposite to each other in the inversion, they will not contribute to the gradient at all; in the other words, the inverted impedance model will not fix the mismatch in these events. Physically these events are generated by refractions and head-waves in the record. Thus, they do not contain information on impedance. By contrast, equations (22) and (23) show that the gradient of velocity has an opposite behavior to the gradient of impedance. When θ is 0°, the gradient vanishes. This means that events propagating opposite to each other will be emphasized in the gradient while the events having small opening angles will be suppressed. Generally, events recorded at small opening angles are from short source-receiver offset reflections. The radiation patterns of impedance and velocity are plotted in figure 1. Figure 1. View largeDownload slide The radiation pattern of the gradient of (a) velocity and (b) impedance. Figure 1. View largeDownload slide The radiation pattern of the gradient of (a) velocity and (b) impedance. Similarly, we can apply the same analysis for a parameterization in terms of velocity and density. The gradient of velocity under this parameterization is |$ begin{eqnarray}{g}_{v}({\bf{x}})=\displaystyle \sum _{\omega }\left(-\displaystyle \frac{2{\omega }^{2}}{\rho ({\bf{x}}){v}^{3}({\bf{x}})}{p}_{s}^{\ast }({\bf{x}},\omega ){p}_{r}({\bf{x}},\omega )\right),\end{equation} $| gv(x)=∑ω−2ω2ρ(x)v3(x)ps⁎(x,ω)pr(x,ω), 24 and the gradient of density is |$ begin{eqnarray}{g}_{\rho }({\bf{x}},\theta )=\displaystyle \sum _{\omega }\left(-\displaystyle \frac{{\omega }^{2}(1+\,\cos \,\theta )}{{\rho }^{2}({\bf{x}}){v}^{2}({\bf{x}})}{p}_{s}^{\ast }({\bf{x}},\omega ){p}_{r}({\bf{x}},\omega )\right).\end{equation} $| gρ(x,θ)=∑ω−ω2(1+cosθ)ρ2(x)v2(x)ps⁎(x,ω)pr(x,ω). 25 According to equations (24) and (25), we can plot out the radiation patterns of the gradients of velocity and density, which are shown in figure 2. Figure 2. View largeDownload slide The radiation pattern of the gradient of (a) velocity and (b) density. Figure 2. View largeDownload slide The radiation pattern of the gradient of (a) velocity and (b) density. Comparing figures 1 and 2 and equations (20), (22), (24) and (25), the radiation pattern of impedance is the same as that of density but different from that of velocity in both the velocity-impedance and velocity-density parameterizations, especially at large angles. Thus, impedance inversion becomes more and more incompatible with the data generated by velocity contrasts with the increase of the angle. In addition, under the condition that velocity is fixed during the impedance inversion, a density can be always found to produce a predicted data identical to the updated impedance. These imply Short offset (small opening angle) data are more suitable for acoustic impedance inversion. Acoustic impedance inversion is more compatible with the data generated by density contrasts than velocity contrasts. A simple two-layer model with a one-source-one-receiver geometry is chosen to demonstrate the four gradients from the two pairs of parameterizations. The caption of figure 3 contains detailed information on the models and on the geometry used in this example. Figures 3(a) and (b) show the velocity and impedance gradients with the velocity-impedance parameterization. It can be observed that the gradients include an ellipse-like part and a rabbit-ear-like part. The ellipse-like part is formed from the direct cross-correlation of the incident source and residue wavefields. This part is the migration component of FWI (Yao et al2014, Yao and Warner 2015, Zhou et al2015, Yao and Wu 2017). Connecting the source and the receiver indicated by the crosses to a point on the ellipse forms the opening angle of the source and residual wavefields. We can see the angle increases from the point at the most left or right of the ellipse to the point at the top or bottom of the ellipse. It is obvious that the amplitude of the impedance gradient decreases with increasing the opening angle, while the velocity gradient has an opposite behavior. This observation is reflected in their radiation patterns shown in figure 1. More interestingly, one can observe that the rabbit-ear parts are prominent in the velocity gradient but nearly vanish in the impedance gradient. This phenomenon can be explained by the radiation pattern shown in figure 1 and equations (20) and (22) because this part of the gradient is generated by the first-order reflections of the source and residual wavefields; thereby, they have an opening angle of 180°. The component of the gradient generated by events that have an opening angle of 180° contains information on the long-wavelengths. Thus, it is useful for updating the background model of velocity, also referred to as the tomography component of FWI (Yao et al2014, Wu and Alkhalifah 2015, Yao and Warner 2015, Zhou et al2015, Yao and Wu 2017). This phenomenon implies that the impedance inversion does not fix the long-wavelength background errors. Figure 3. View largeDownload slide The gradient of a two-layer model with a one-source-one-receiver configuration. The velocity and density are 2000 m s−1 and 1000 kg m−3 for the top layer, 3000 m s−1 and 1500 kg m−3 for the bottom layer, respectively. The interface is at the depth of 4 km indicated by the dashed line in each subfigure. The source and the receiver are located at the distance of 1.25 km and 3.75 km, respectively. Both the source and the receiver are set at the depth of 2.25 km. They are indicated by the blue and green crosses in each subfigure. A 10 Hz Ricker is used as the source wavelet. The initial model is also a two-layer model. It has a velocity of 1950 m s−1 and a density of 950 kg m−3 for the top layer, 2925 m s−1 and 1425 kg m−3 for the bottom layer. (a), (b) Show the gradient of velocity and impedance respectively under the velocity-impedance parameterization. (c), (d) Show the gradient of velocity and density respectively under the velocity-density parameterization. All the gradients are normalized and the amplitude is clipped from −1 to 1. Figure 3. View largeDownload slide The gradient of a two-layer model with a one-source-one-receiver configuration. The velocity and density are 2000 m s−1 and 1000 kg m−3 for the top layer, 3000 m s−1 and 1500 kg m−3 for the bottom layer, respectively. The interface is at the depth of 4 km indicated by the dashed line in each subfigure. The source and the receiver are located at the distance of 1.25 km and 3.75 km, respectively. Both the source and the receiver are set at the depth of 2.25 km. They are indicated by the blue and green crosses in each subfigure. A 10 Hz Ricker is used as the source wavelet. The initial model is also a two-layer model. It has a velocity of 1950 m s−1 and a density of 950 kg m−3 for the top layer, 2925 m s−1 and 1425 kg m−3 for the bottom layer. (a), (b) Show the gradient of velocity and impedance respectively under the velocity-impedance parameterization. (c), (d) Show the gradient of velocity and density respectively under the velocity-density parameterization. All the gradients are normalized and the amplitude is clipped from −1 to 1. Figures 3(c) and (d) show the gradients of velocity and density with the velocity-density parameterization. By comparison, it can be seen that the amplitude of the velocity gradient does not decrease nor increase with the opening angle while it decreases gradually for the density gradient. Furthermore, the gradient of density is identical to that of impedance. As a result, the impedance gradient has an identical response as the density gradient for all opening angles but a similar response as the velocity gradient for only the small opening angle, which corresponds to small offsets. If the starting model is accurate enough, such that the region of convergence is located in the basin of attraction of the global minimum, then an optimal solution can be found by minimizing ϕ with localized gradient methods. The first-order optimality conditions are begin{eqnarray}\left\{\begin{array}{c}{g}_{v}=0,\\ {g}_{I}=0.\end{array}\right.\end{equation} gv=0,gI=0. 26 However, the velocity is fixed during acoustic impedance inversion with FWI, which means the inversion only achieves the second optimality condition in equation (26); therefore the first optimality condition is a precondition for the acoustic impedance inversion, which means an accurate starting velocity model is expected for the impedance inversion. As a result, equation (26) implies that Acoustic impedance inversion requires very accurate starting velocity models for achieving a high-quality acoustic impedance inversion. Sensitivity tests on the Marmousi model In this section, we demonstrate the points made above on the sensitivity of impedance inversion with FWI utilizing the Marmousi model. The velocity and density of the Marmousi model are shown in figures 4(a) and (b). Its impedance model, which is the product of velocity and density, is shown in figure 4(c). For generating the surface record and the starting models of inversion, the true models are smoothed by using the SeismicUnix command, ‘smooth2’, with the parameters of 10 and 60. The smoothed version of velocity and density with the parameter of 60 is shown in figures 5(a) and (b), respectively. Figures 5(c) and (d) are the counterpart of figures 5(a) and (b) with the parameter set to 10. One can observe that setting the parameter to 60 gives much smoother models than that of 10. Thus, we consider figures 5(a) and (b) as inaccurate starting models but figures 5(c) and (d) as accurate starting models. Figure 4. View largeDownload slide The (a) velocity and (b) density of Marmousi model. (c) The impedance of Marmousi model, which is the product of velocity and density. Figure 4. View largeDownload slide The (a) velocity and (b) density of Marmousi model. (c) The impedance of Marmousi model, which is the product of velocity and density. Figure 5. View largeDownload slide Two pairs of initial models: (a) and (b) are the velocity and density models, respectively, which are generated by smoothing the true models with the smooth parameter of 60; (c) and (d) are the counterpart of (a) and (b) for the smooth parameter of 10. The smoothing process is achieved with the SeismicUnix command ‘smooth2’. Figure 5. View largeDownload slide Two pairs of initial models: (a) and (b) are the velocity and density models, respectively, which are generated by smoothing the true models with the smooth parameter of 60; (c) and (d) are the counterpart of (a) and (b) for the smooth parameter of 10. The smoothing process is achieved with the SeismicUnix command ‘smooth2’. We generated a data set of 255 shots with a spacing of 50 m. The source of the first shot is located at the distance of 187.5 m (figure 4). The source wavelet is a 10 Hz Ricker wavelet. The data for each shot is recorded with a fixed geometry array of 1051 receivers across the whole model and with a spacing of 12.5 m, which is also the cell size of the finite difference in these experiments. All sources and receivers are fixed at the depth of 25 m. In order to mimic the marine environment, the free surface is applied to generate surface-related multiples in the record. Figure 6(a) shows a shot gather generated with the true velocity and density while the shot in figure 6(b) is generated from the true velocity and the heavily smoothed density model (figure 5(b)). The shot in figure 6(c) is generated from the true density and the heavily smoothed velocity model (figure 5(a)). Thereby, the reflections in figure 6(a) are from the contrast of both velocity and density, the reflections in figure 6(b) are mainly from the contrast of velocity, and the reflections in figure 6(c) are mainly from the contrast of density. Comparing the three shot gathers, it is obvious that the one in figure 6(b) is more similar to the one in figure 6(a), especially the refracted energy in the far-offset region, than figure 6(c), but figure 6(c) has quite strong reflections in short offsets. This explains one reason why far-offset refraction data are favorable for FWI to recover the velocity model (Warner et al2013). Figure 6. View largeDownload slide One shot record gather from (a) the true velocity and density models, (b) the true velocity model and the smoothed density model shown in figures 5(b) and (c) the true density model and the smoothed velocity model shown in figure 5(a). The source is located at the distance of 6.175 km. All the displays are clipped to the same range. Figure 6. View largeDownload slide One shot record gather from (a) the true velocity and density models, (b) the true velocity model and the smoothed density model shown in figures 5(b) and (c) the true density model and the smoothed velocity model shown in figure 5(a). The source is located at the distance of 6.175 km. All the displays are clipped to the same range. We carried out two group tests. The first group test is with a range of offsets, i.e. 7.5, 5 and 3.75 km, and two pairs of initial velocity and density models shown in figure 5. The inversions are carried out with 30 iterations. In the case of heavily smoothed initial models, which are shown in figures 5(a) and (b), the multi-scale strategy (Bunks et al1995) is applied. The inversion is up-scaled in frequency at 3, 4, 6 and 8 Hz sequentially, and the algorithm is iterated 5 times per frequency band. Note that each frequency means a narrow frequency band. In the 20 multi-scale iterations, the inversion has recovered the impedance from low wavenumbers to high wavenumbers. The last 10 iterations use the full bandwidth to finalize the inversion. In the case of slightly smoothed initial models, which are shown in figures 5(c) and (d), the inversion uses the full bandwidth of data for all the 30 iterations because the initial models are relatively accurate. Figure 7 shows the inverted models. According to the results, we can observe four interesting phenomena. First, the artifacts, e.g. vertical strips, are gradually reduced with the decrease of offsets. Secondly, the value of impedance becomes closer to the true value when the offsets used for inversion are decreased. Third, increasing accuracy in the starting models gives much better inversion results. Finally, as the velocity model is fixed during the inversion, the impedance inversion hardly makes any correction to the location of the reflectors. The first two phenomena are caused by the fact that the forward modeling of impedance inversion with FWI is through density since the velocity is fixed during the inversion. As a result, the AVO cannot be correct if the data are generated by velocity contract. The residuals caused by incorrect AVO then will be mapped into vertical strips and wrong impedance value by inversion. The mechanism of the generation of the vertical strips can be explained as: since the AVO is inaccurate especially for the far offsets, and the impedance inversion mainly uses the migration component of FWI, the residual of each traces will be mapped into an ellipse-like shape, and further stack of all of them removes their horizontal parts but remains some of their vertical parts in the model. The reason for the last two phenomena is that the gradient of the impedance of FWI suppresses the tomography component of the FWI gradient and therefore is insensitive to the background velocity errors but it is sensitive to the impedance perturbation, which is directly related to the contrast of velocity and density between different layers. Consequently, as the starting velocity model is fixed, the impedance inversion with FWI will not fix the errors of background velocity. This test demonstrated the first and third statements of the sensitivity of the impedance inversion with FWI: short-offset data and accurate starting velocity models are crucial for a high-quality impedance inversion with FWI. Figure 7. View largeDownload slide Impedance inversion results with different initial models and source-receiver maximum offset. (a), (c) and (e) show the inverted impedance models starting from a smoothed version of the true velocity and density models with the smoothing parameter of 60. (b), (d) and (f) are the same as (a), (c) and (e), respectively, but the smoothing parameter for generating the initial models is 10. The maximum offset of records used for inversion is 7.5 km for (a) and (b), 5 km for (c) and (d), and 3.75 km for (e) and (f). The color scale is the same as figure 4(c). Figure 7. View largeDownload slide Impedance inversion results with different initial models and source-receiver maximum offset. (a), (c) and (e) show the inverted impedance models starting from a smoothed version of the true velocity and density models with the smoothing parameter of 60. (b), (d) and (f) are the same as (a), (c) and (e), respectively, but the smoothing parameter for generating the initial models is 10. The maximum offset of records used for inversion is 7.5 km for (a) and (b), 5 km for (c) and (d), and 3.75 km for (e) and (f). The color scale is the same as figure 4(c). In the second group of tests, the data for inversion are generated either using the true velocity model with heavily smoothed density or heavily smoothed velocity with true density. As a result, the reflection data are generated by either the velocity contrast or the density contrast. To generate the record, we used one true model and one heavily smoothed model; in the inversion, to form the starting models, we use the models generating the record but smooth the true model slightly with the parameter of 10. The inversion is carried out with 30 iterations and full bandwidth is used in each iteration. The inverted impedance model is shown in figure 8. By comparison, we can observe that the inversion leads to a better-reconstructed impedance model with data generated by density contrast than velocity contrast. In addition, the offset range does not affect the result for the data generated by density contrast but the large offset range for the data generated by velocity contrast tends to introduce more artifacts and inaccurate value of impedance than the short-offset range. These phenomena are caused by the fact that the forward modeling of impedance inversion with FWI is through density. They can also be explained by the fact that the radiation pattern of impedance is identical to that of density rather than that of velocity. As a result, data generated by density contrast are more compatible than data generated by velocity contrast with the impedance inversion. In addition, short-offset data from velocity contrast are more compatible with the impedance inversion than its large offsets. This experiment demonstrates the second statement of the sensitivity analysis: acoustic impedance inversion with FWI is more compatible with the data generated by density contrast than velocity contrast. Figure 8. View largeDownload slide Impedance inversion with the data generated with (a), (c) heavily smoothed true velocity and true density or (b), (d) true velocity and heavily smoothed true density. The input data for inversion have a maximum offset of 7.5 km for (a) and (b) but 3.75 km for (c) and (d). The smooth parameter is 60 for generating the heavily smoothed models. In the inversion for (a), (c), the starting velocity model is the heavily smoothed velocity model, thus it is already correct; however, the starting density model is a smoothed true density model with the smooth parameter of 10. On the contrary, the inversion for (b), (d) starts from a correct density model but the starting velocity is a smoothed true velocity model with the smooth parameter of 10. (e) Shows the true impedance model for (a), (c) while (f) shows the true impedance model for (b), (d). The color scale is the same as figure 4(c). Figure 8. View largeDownload slide Impedance inversion with the data generated with (a), (c) heavily smoothed true velocity and true density or (b), (d) true velocity and heavily smoothed true density. The input data for inversion have a maximum offset of 7.5 km for (a) and (b) but 3.75 km for (c) and (d). The smooth parameter is 60 for generating the heavily smoothed models. In the inversion for (a), (c), the starting velocity model is the heavily smoothed velocity model, thus it is already correct; however, the starting density model is a smoothed true density model with the smooth parameter of 10. On the contrary, the inversion for (b), (d) starts from a correct density model but the starting velocity is a smoothed true velocity model with the smooth parameter of 10. (e) Shows the true impedance model for (a), (c) while (f) shows the true impedance model for (b), (d). The color scale is the same as figure 4(c). A practical workflow for the acoustic impedance inversion In the previous sections, we demonstrated the necessary conditions for successful impedance inversion with an FWI algorithm, through numerical experiments. In this section, we introduce a workflow for obtaining high-quality acoustic impedance inversion with FWI in practice. Firstly, travel-time tomography or reflection waveform inversion (RWI) (Yao et al2014, Wu and Alkhalifah 2015, Yao and Warner 2015, Zhou et al2015, Yao and Wu 2017) is applied to build an accurate background velocity model with correct long wavelength; secondly, FWI is used to fix the intermediate error of velocity with Gardner’s density and all-offset data; finally, FWI is employed to invert the acoustic impedance with the short-offset data starting from the FWI velocity and Gardner’s density. We need the first step to fix the background velocity error. The reason is that FWI fixes the background velocity error mainly through the refractions for the surface acquisitions because refractions mainly contribute to the tomography component of FWI, which is explained in the sensitivity analysis section. However, very long offsets are needed to capture the refractions from deep targets. Unfortunately, the maximum offset of the modern surface seismic surveys is usually not long enough to capture the refractions from the deep targets. In the second step, we update the velocity with the velocity gradient of the velocity-density parameterization instead of the velocity gradient of the velocity-impedance parameterization, which is equation (13). This is because Gardner’s density is used in each iteration as prior information, and thus the velocity-density parameterization is a natural choice. Besides, our numerical tests show that the velocity gradient shown in equation (13) with Gardner’s density produces less accurate velocity models as well as more artifacts around edges. In this step, the refractions are used for fixing the remaining long-wavelength background errors whilst reflections aim to correct the intermediate wavelength errors. In the final step, the impedance is updated with the impedance gradient expressed in equation (12). The sensitivity analyses have shown that the impedance gradient does not sense long-wavelength background errors but is effective to fix the short-wavelength errors. In addition, the second step aims to minimize the objective function with respect to velocity while the final step is for minimization with respect to impedance. Thus, the last two steps are complementary to each other. We apply this workflow to invert the Marmousi model. The data are generated from the true velocity and true density models. The geometry is the same as that used in the previous section. Figure 6(a) shows a shot gather generated for this model. We use the heavily smoothed velocity model (figure 5(a)) to mimic the result from travel-time tomography or RWI, which can recover the long-wavelength background. FWI is then applied to refine the velocity model with full-offset data and Gardner’s density, which is defined as |$ begin{eqnarray}\rho =310\,{v}^{0.25},\end{equation} $| ρ=310v0.25, 27 where v is the rock velocity in m s−1 (Gardner et al1974). Gardner’s relation (equation (27)) is only for rocks while the water density is set as 1000 kg m−3 in the inversion. The inversion includes 30 iterations and uses the same multi-scale strategy as the one used in the previous section. The recovered velocity model is shown in figure 9(a). Figure 9(b) shows the corresponding impedance model of figure 9(a) with its Gardner’s density. By comparison with the true impedance model shown in figure 4(c), the impedance model has correct structures and the positions of each layer are correct; however, the fine layers are not recovered. Starting from this result, acoustic impedance inversion with FWI gives the result shown in figure 9(c). In the inversion, full bandwidth and offsets smaller than 3.75 km are used in all the 30 iterations. By comparison, it is obvious that the impedance inversion recovers the fine layers by accurately reconstructing the interfaces. From this test, we can see that in this workflow, travel-time tomography or RWI recovers the long-wavelength component of the model, FWI with Gardner’s density fixes the intermediate-wavelength component, and the acoustic impedance inversion with FWI builds the short-wavelength component. With the cascading inversion, the full spectra of models are recovered. Figure 9. View largeDownload slide The cascading inversion of FWI. (a) The inverted velocity model of FWI started from a smoothed version of the true velocity model. The smooth parameter is 60. In the inversion, the density is computed by using Gardner’s relation. (b) Shows the corresponding impedance of (a), which is the product of (a) and its Gardner’s density. (c) Further impedance inversion result with FWI started from (a) and its Gardner’s density. The color scales are identical to those shown in figure 4. Figure 9. View largeDownload slide The cascading inversion of FWI. (a) The inverted velocity model of FWI started from a smoothed version of the true velocity model. The smooth parameter is 60. In the inversion, the density is computed by using Gardner’s relation. (b) Shows the corresponding impedance of (a), which is the product of (a) and its Gardner’s density. (c) Further impedance inversion result with FWI started from (a) and its Gardner’s density. The color scales are identical to those shown in figure 4. The inversion results using the workflow show well reconstructed models of impedance. The convergence can be observed from the predicted data and the residuals. Comparing the predicted data (figure 10(b)) from the FWI velocity model (figure 9(a)) with the observed data (figure 10(a)), it can be seen that the far-offset data, especially the refractions, is matched very well, but the short-offset data, especially the reflections, does not have a good fit. However, the acoustic impedance inversion with FWI improves the fitness of the short-offset data, which can be seen by comparison of figures 10(a)–(c). The convergence also can be seen from the residuals, which are shown in figure 11. Figure 10. View largeDownload slide (a) One shot record gather from the true velocity and density models. The source is located at the distance of 6.175 km. (b) The corresponding predicted shot gather from the inverted velocity model shown in figure 9(a) and its Gardner density. (c) The corresponding predicted shot gather from the inverted impedance model shown in figure 9(c). All the displays are clipped to the same range. Figure 10. View largeDownload slide (a) One shot record gather from the true velocity and density models. The source is located at the distance of 6.175 km. (b) The corresponding predicted shot gather from the inverted velocity model shown in figure 9(a) and its Gardner density. (c) The corresponding predicted shot gather from the inverted impedance model shown in figure 9(c). All the displays are clipped to the same range. Figure 11. View largeDownload slide The corresponding residual of the predicted shot gather shown in figure 10. (a) Residual from the inverted velocity model and its Gardner’s density; (b) residual from the inverted impedance model. The vertical dashed lines indicate the maximum offset used for the impedance inversion. Figure 11. View largeDownload slide The corresponding residual of the predicted shot gather shown in figure 10. (a) Residual from the inverted velocity model and its Gardner’s density; (b) residual from the inverted impedance model. The vertical dashed lines indicate the maximum offset used for the impedance inversion. Conclusions In this paper, we analyzed the sensitivity of acoustic impedance inversion with FWI. In the inversion, the wave equation is parameterized as velocity and impedance. During the inversion, the velocity is fixed. Consequently, all the residuals are mapped into the impedance model. The radiation pattern shows that the impedance has high sensitivity at small opening angles, which corresponds to short offsets. This implies that data recorded in the short offsets are more suitable for impedance inversion with FWI. Since the velocity is fixed, the forward modeling for the impedance inversion, in fact, changes density to simulate the update of impedance. In addition, the impedance gradient has an identical radiation pattern as the density gradient. As a result, the impedance inversion is more compatible with the data generated from density contrast than velocity contrast. Since the radiation pattern shows that density and velocity have a similar response in the small opening angle, we should choose the short-offset data for the best compatibility of velocity contrast with the impedance inversion. We expect the inversion to optimize the objective function, which is a function of both velocity and impedance. Besides, the velocity is fixed during the impedance inversion. As a result, an accurate starting velocity model is a precondition for a high-quality impedance inversion. In order to optimize the objective function, we propose a cascading inversion workflow: firstly, using tomography or RWI to recover a long-wavelength background velocity model; secondly, FWI with Gardner’s density fixes the remaining long-wavelength background velocity errors as well as the intermediate wavelength velocity errors; finally, acoustic impedance inversion with FWI recovers the short-wavelength impedance interfaces. All these analyses and the workflow have been validated by the numerical experiments. We expect the concluded sensitivity analyses and the cascading workflow will give a guide to the applications to field data. Acknowledgments The authors are grateful to the NSFC (Grant No. 41504106, 41630209), National Basic Research Program of China (973 Program) (No.2015CB250903) and Science Foundation of China University of Petroleum, Beijing (No. 2462018BJC001) for partially supporting this research. The authors would like to thank the editor and Dr Samuel Gray and another anonymous reviewer for their very meaningful comments and suggestions, which help to improve and clarify the manuscript significantly. The authors also thank Adrian Umpleby for the discussion about the adjoint operator of wave equations. Appendix A. Derivation of the adjoint operator of the acoustic wave equation In this appendix, we demonstrate the adjoint operator of the acoustic wave equation with variable density. The procedure is similar to appendix A in da Silva et al (2016). If the acoustic wave equation, |$ begin{eqnarray}\left(\displaystyle \frac{1}{\rho ({\bf{x}}){v}^{2}({\bf{x}})}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{\rm{\nabla }}\cdot \displaystyle \frac{1}{\rho ({\bf{x}})}{\rm{\nabla }}\right)p({\bf{x}},t)=f({\bf{x}},t),\end{equation} $| 1ρ(x)v2(x)∂2∂t2−∇⋅1ρ(x)∇p(x,t)=f(x,t), 1 is solved by using explicitly finite-difference methods with the second-order accuracy for the time derivative, then the semi-discretized form of equation (1) becomes |$ begin{eqnarray}\displaystyle \frac{1}{\rho {v}^{2}{\rm{\Delta }}{t}^{2}}p(t+{\rm{\Delta }}t)+\left(-\displaystyle \frac{2}{\rho {v}^{2}{\rm{\Delta }}{t}^{2}}-{\rm{\nabla }}\cdot \displaystyle \frac{1}{\rho }{\rm{\nabla }}\right)p(t)\\ \,\,+\,\displaystyle \frac{1}{\rho {v}^{2}{\rm{\Delta }}{t}^{2}}p(t-{\rm{\Delta }}t)=f({\bf{x}},t),\end{equation} $| 1ρv2Δt2p(t+Δt)+−2ρv2Δt2−∇⋅1ρ∇p(t)+1ρv2Δt2p(t−Δt)=f(x,t), A.1 where Δt is the time-stepping interval. After discretizing | ${\rm{\nabla }}\cdot \tfrac{1}{\rho }{\rm{\nabla }}$ | ∇⋅1ρ∇ in equation (A.1), the fully-discretized wave equation (1) can be written in a matrix form as \begin{equation}\left[\begin{array}{ccccccc}{\bf{D}} & & & & & & \\ {\bf{B}} & {\bf{D}} & & & & & \\ {\bf{D}} & {\bf{B}} & {\bf{D}} & & & & \\ & \ddots & \ddots & \ddots & & & \\ & & {\bf{D}} & {\bf{B}} & {\bf{D}} & & \\ & & & {\bf{D}} & {\bf{B}} & {\bf{D}} & \\ & & & & {\bf{D}} & {\bf{B}} & {\bf{D}}\end{array}\right]\left[\begin{array}{c}{\bf{p}}({t}_{1})\\ {\bf{p}}({t}_{2})\\ {\bf{p}}({t}_{3})\\ \vdots \\ {\bf{p}}({t}_{nt-2})\\ {\bf{p}}({t}_{nt-1})\\ {\bf{p}}({t}_{nt})\end{array}\right]=\left[\begin{array}{c}{\bf{f}}({t}_{0})\\ {\bf{f}}({t}_{1})\\ {\bf{f}}({t}_{2})\\ \vdots \\ {\bf{f}}({t}_{nt-3})\\ {\bf{f}}({t}_{nt-2})\\ {\bf{f}}({t}_{nt-1})\end{array}\right],\end{equation} DBDDBD⋱⋱⋱DBDDBDDBDp(t1)p(t2)p(t3)⋮p(tnt−2)p(tnt−1)p(tnt)=f(t0)f(t1)f(t2)⋮f(tnt−3)f(tnt−2)f(tnt−1), A.2 where | $nt$ | nt is the total time steps. Note that the matrix on the left-hand side (lhs) of equation (A.2) is never explicitly formed in the numerical implementation. It just denotes how the wavefield evolves with time stepping. The matrix equation can be denoted symbolically as |$ begin{eqnarray}{\bf{A}}{\bf{p}}={\bf{f}}.\end{equation} $| Ap=f. 7 In equation (A.2), D is a diagonal matrix, the value of whose diagonal elements is | $\displaystyle \frac{1}{\rho {v}^{2}{\rm{\Delta }}{t}^{2}},$ | 1ρv2Δt2,B is a symmetric matrix, the value of which is computed by discretizing | $-\displaystyle \frac{2}{\rho {v}^{2}{\rm{\Delta }}{t}^{2}}-{\rm{\nabla }}\cdot \displaystyle \frac{1}{\rho }{\rm{\nabla }}$ | −2ρv2Δt2−∇⋅1ρ∇ with finite difference. Since blocks D and B on the lhs of equation (A.2) handle the spatial part of the large matrix A, so their size is nx × nx, where | $nx$ | nx is the total number of spatial grid points of the discretized model. There are | $nt$ | nt rows blocks in the matrix | ${\bf{A}}.$ | A. Each | ${\bf{p}}({t}_{i})$ | p(ti) and | ${\bf{f}}({t}_{i})$ | f(ti) are column vectors of the pressure wavefield and the source for the time step | $i$ | i at each spatial grid points of the model, and can be explicitly expressed as begin{eqnarray}{\bf{p}}({t}_{i}) & = & \left[\begin{array}{cc}p({{\bf{x}}}_{1},{t}_{i}) & p({{\bf{x}}}_{2},{t}_{i})\,\begin{array}{cc}\cdots & p({{\bf{x}}}_{j},{t}_{i})\end{array}\end{array}\right.\\ & & \begin{array}{cc}\cdots & {p({{\bf{x}}}_{nx},{t}_{i})]}^{{\rm{T}}}\end{array}\end{equation} p(ti)=p(x1,ti)p(x2,ti)⋯p(xj,ti)⋯p(xnx,ti)]T A.3 and begin{eqnarray}{\bf{f}}({t}_{i}) & = & \left[\begin{array}{cc}f({{\bf{x}}}_{1},{t}_{i}) & f({{\bf{x}}}_{2},{t}_{i})\,\begin{array}{cc}\cdots & f({{\bf{x}}}_{j},{t}_{i})\end{array}\end{array}\right.\\ & & {\left.\begin{array}{cc}\cdots & f({{\bf{x}}}_{nx},{t}_{i})\end{array}\right]}^{{\rm{T}}},\end{equation} f(ti)=f(x1,ti)f(x2,ti)⋯f(xj,ti)⋯f(xnx,ti)T, A.4 respectively, where T means transpose and | ${{\bf{x}}}_{j}$ | xj indicates the position of the | $j$ | jth grid points of the discretized model that has | $nx$ | nx grid points in total. As a result, the total length of | ${\bf{p}}$ | p and | ${\bf{f}}$ | f is | $nx\times nt,$ | nx×nt, and their most rapidly-varying index is space. Equation (A.2) encompasses the initial conditions (equation (3)) by setting | ${\bf{p}}({t}_{-1})={\bf{p}}({t}_{0})=0.$ | p(t−1)=p(t0)=0. As can be seen, | ${\bf{A}}$ | A is a lower triangular matrix, so equation (A.2) should be solved from top to bottom; in other words, it should be solved forwards in time, which is called forward propagation. According to equations (7) and (A.2), the adjoint wave equation can be written as |$ begin{eqnarray}{{\bf{A}}}^{{\rm{T}}}\delta {\bf{p}}=\delta {\bf{f}},\end{equation} $| ATδp=δf, A.5 where | $\delta {\bf{f}}$ | δf and | $\delta {\bf{p}}$ | δp are the column vectors of the adjoint source and the adjoint wavefield. Equation (A.5) can be explicitly expressed as begin{eqnarray}\left[\begin{array}{ccccccc}{\bf{D}} & {{\bf{B}}}^{{\rm{T}}} & {\bf{D}} & & & & \\ & {\bf{D}} & {{\bf{B}}}^{{\rm{T}}} & {\bf{D}} & & & \\ & & {\bf{D}} & {{\bf{B}}}^{{\rm{T}}} & {\bf{D}} & & \\ \, & & & \ddots & \ddots & \ddots & \\ & & & & {\bf{D}} & {{\bf{B}}}^{{\rm{T}}} & {\bf{D}}\\ & & & & & {\bf{D}} & {{\bf{B}}}^{{\rm{T}}}\\ & & & & & & {\bf{D}}\end{array}\right]\left[\begin{array}{c}\delta {\bf{p}}({t}_{0})\\ \delta {\bf{p}}({t}_{1})\\ \delta {\bf{p}}({t}_{2})\\ \vdots \\ \delta {\bf{p}}({t}_{nt-3})\\ \delta {\bf{p}}({t}_{nt-2})\\ \delta {\bf{p}}({t}_{nt-1})\end{array}\right]=\left[\begin{array}{c}\delta {\bf{f}}({t}_{1})\\ \delta {\bf{f}}({t}_{2})\\ \delta {\bf{f}}({t}_{3})\\ \vdots \\ \delta {\bf{f}}({t}_{nt-2})\\ \delta {\bf{f}}({t}_{nt-1})\\ \delta {\bf{f}}({t}_{nt})\end{array}\right],\end{equation} DBTDDBTDDBTD⋱⋱⋱DBTDDBTDδp(t0)δp(t1)δp(t2)⋮δp(tnt−3)δp(tnt−2)δp(tnt−1)=δf(t1)δf(t2)δf(t3)⋮δf(tnt−2)δf(tnt−1)δf(tnt), A.6 where the diagonal matrix | ${\bf{D}}$ | D is the same as that in equation (A.2) but | ${{\bf{B}}}^{{\rm{T}}}$ | BT is the transpose of the matrix | ${\bf{B}}.$ | B. Similar to equation (A.2), we used the condition: | ${\bf{p}}({t}_{nt})\,={\bf{p}}({t}_{nt+1})=0$ | p(tnt)=p(tnt+1)=0 for deriving equation (A.6). Since the matrix | ${\bf{B}}$ | B is symmetric, equation (A.6) becomes begin{eqnarray}\left[\begin{array}{ccccccc}{\bf{D}} & {\bf{B}} & {\bf{D}} & & & & \\ & {\bf{D}} & {\bf{B}} & {\bf{D}} & & & \\ & & {\bf{D}} & {\bf{B}} & {\bf{D}} & & \\ \, & & & \ddots & \ddots & \ddots & \\ & & & & {\bf{D}} & {\bf{B}} & {\bf{D}}\\ & & & & & {\bf{D}} & {\bf{B}}\\ & & & & & & {\bf{D}}\end{array}\right]\left[\begin{array}{c}\delta {\bf{p}}({t}_{0})\\ \delta {\bf{p}}({t}_{1})\\ \delta {\bf{p}}({t}_{2})\\ \vdots \\ \delta {\bf{p}}({t}_{nt-3})\\ \delta {\bf{p}}({t}_{nt-2})\\ \delta {\bf{p}}({t}_{nt-1})\end{array}\right]=\left[\begin{array}{c}\delta {\bf{f}}({t}_{1})\\ \delta {\bf{f}}({t}_{2})\\ \delta {\bf{f}}({t}_{3})\\ \vdots \\ \delta {\bf{f}}({t}_{nt-2})\\ \delta {\bf{f}}({t}_{nt-1})\\ \delta {\bf{f}}({t}_{nt})\end{array}\right].\end{equation} DBDDBDDBD⋱⋱⋱DBDDBDδp(t0)δp(t1)δp(t2)⋮δp(tnt−3)δp(tnt−2)δp(tnt−1)=δf(t1)δf(t2)δf(t3)⋮δf(tnt−2)δf(tnt−1)δf(tnt). A.7 As | ${{\bf{A}}}^{{\rm{T}}}$ | AT is an upper triangular matrix, equation (A.7) should be solved backward in time, which is called backward propagation. Comparing equations (A.2) and (A.7), one can observe that since wave equation (1) is formulated to be symmetric (i.e. self-adjoint) in space, the same code for | ${\bf{A}}$ | A(i.e. forward modeling) can be used for its adjoint operator | ${{\bf{A}}}^{{\rm{T}}}$ | AT (i.e. back propagation). Similarly, we can reformulate the elastic wave equation to have a symmetric operator in space, which can be found in equation (4.11) of Guasch (2011). However, when more physical parameters, e.g. anisotropy (equations (A.5) and (A.7) in da Silva et al2016) and attenuation, are introduced into wave equations, it becomes more difficult to design symmetric wave-equation operators in space. Usually, the operators are asymmetric, so the gradients of FWI formulated by using the same operator for both forward and backward propagation are not the true gradients but points to a similar direction. Consequently, updating models with the computed gradients still decreases the value of the objective function. In other words, more iterations still converge the objective function and produce good results. This assertion has been verified by many other authors (e.g. Warner et al2013). Appendix B. Derivation of the diffractor response The acoustic wave equation with impedance and velocity parameterization can be expressed as |$ begin{eqnarray}\left(\displaystyle \frac{1}{I({\bf{x}})v({\bf{x}})}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{I({\bf{x}})}{\rm{\nabla }}\right)p({\bf{x}},t)=f({\bf{x}},t).\end{equation} $| 1I(x)v(x)∂2∂t2−∇⋅v(x)I(x)∇p(x,t)=f(x,t). 6 If the impedance | $I$ | I is decomposed as background impedance | ${I}_{0}$ | I0 and impedance perturbation | $\delta I,$ | δI, and similarly, the wavefield | $p$ | p is decomposed as background wavefield | ${p}_{0}$ | p0 and wavefield perturbation | $\delta p,$ | δp, then using equation (2) and applying temporal Fourier transform on equation (6), it then becomes |$ begin{eqnarray}\left(-\displaystyle \frac{{\omega }^{2}}{({I}_{0}({\bf{x}})+\delta I({\bf{x}}))v({\bf{x}})}-{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{{I}_{0}({\bf{x}})+\delta I({\bf{x}})}{\rm{\nabla }}\right)\\ \,\times \,({p}_{0}({\bf{x}},\omega )+\delta p({\bf{x}},\omega ))=s(t)\delta ({\bf{x}}-{{\bf{x}}}_{s}).\end{equation} $| −ω2(I0(x)+δI(x))v(x)−∇⋅v(x)I0(x)+δI(x)∇×(p0(x,ω)+δp(x,ω))=s(t)δ(x−xs). B.1 Applying the first-order Taylor expansion, one can have |$ begin{eqnarray}\displaystyle \frac{1}{{I}_{0}({\bf{x}})+\delta I({\bf{x}})}\approx \displaystyle \frac{1}{{I}_{0}({\bf{x}})}-\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}.\end{equation} $| 1I0(x)+δI(x)≈1I0(x)−δI(x)I02(x). B.2 Inserting equation (B.2) into equation (B.1) yields an equation system, begin{eqnarray}\left\{\begin{array}{l}\left(-\displaystyle \frac{{\omega }^{2}}{{I}_{0}({\bf{x}})v({\bf{x}})}-{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{{I}_{0}({\bf{x}})}{\rm{\nabla }}\right){p}_{0}({\bf{x}},\omega )=s(t)\delta ({\bf{x}}-{{\bf{x}}}_{s}),\\ \left(-\displaystyle \frac{{\omega }^{2}}{{I}_{0}({\bf{x}})v({\bf{x}})}-{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{{I}_{0}({\bf{x}})}{\rm{\nabla }}\right)\delta p({\bf{x}},\omega )\\ =\left(-{\omega }^{2}\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}\displaystyle \frac{1}{v({\bf{x}})}-{\rm{\nabla }}\cdot \displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}}){\rm{\nabla }}\right)({p}_{0}({\bf{x}},\omega )+\delta p({\bf{x}},\omega )).\end{array}\right.\end{equation} −ω2I0(x)v(x)−∇⋅v(x)I0(x)∇p0(x,ω)=s(t)δ(x−xs),−ω2I0(x)v(x)−∇⋅v(x)I0(x)∇δp(x,ω)=−ω2δI(x)I02(x)1v(x)−∇⋅δI(x)I02(x)v(x)∇(p0(x,ω)+δp(x,ω)). B.3 Applying Born approximation, i.e. ignoring | $\delta p$ | δp at the right-hand side of the second equation of (B.3) gives begin{eqnarray}\left\{\begin{array}{l}\left(-\displaystyle \frac{{\omega }^{2}}{{I}_{0}({\bf{x}})v({\bf{x}})}-{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{{I}_{0}({\bf{x}})}{\rm{\nabla }}\right)({p}_{0}({\bf{x}},\omega ))=s(t)\delta ({\bf{x}}-{{\bf{x}}}_{s}),\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\\ \left(-\displaystyle \frac{{\omega }^{2}}{{I}_{0}({\bf{x}})v({\bf{x}})}-{\rm{\nabla }}\cdot \displaystyle \frac{v({\bf{x}})}{{I}_{0}({\bf{x}})}{\rm{\nabla }}\right)\delta p({\bf{x}},\omega )\\ =\left(-{\omega }^{2}\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}\displaystyle \frac{1}{v({\bf{x}})}-{\rm{\nabla }}\cdot \displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}}){\rm{\nabla }}\right){p}_{0}({\bf{x}},\omega ).\end{array}\right.\end{equation} −ω2I0(x)v(x)−∇⋅v(x)I0(x)∇(p0(x,ω))=s(t)δ(x−xs),−ω2I0(x)v(x)−∇⋅v(x)I0(x)∇δp(x,ω)=−ω2δI(x)I02(x)1v(x)−∇⋅δI(x)I02(x)v(x)∇p0(x,ω). B.4 Equation (B.4) can be expressed in an integral form with Green’s functions as begin{eqnarray}\left\{\begin{array}{l}{p}_{0}({\bf{x}},\omega )={G}_{s}({\bf{x}}| {{\bf{x}}}_{s},\omega )s(t),\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\\ \delta p({{\bf{x}}}_{r},\omega )=\displaystyle \int {G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega )\\ \left(-{\omega }^{2}\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}\displaystyle \frac{1}{v({\bf{x}})}-{\rm{\nabla }}\cdot \displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}}){\rm{\nabla }}\right){p}_{0}({\bf{x}},\omega ){\rm{d}}{\bf{x}},\end{array}\right.\end{equation} p0(x,ω)=Gs(x∣xs,ω)s(t),δp(xr,ω)=∫Gr(xr∣x,ω)−ω2δI(x)I02(x)1v(x)−∇⋅δI(x)I02(x)v(x)∇p0(x,ω)dx, B.5 where | ${G}_{s}({\bf{x}}| {{\bf{x}}}_{s},\omega )$ | Gs(x∣xs,ω) is the Green’s function from the source position | ${{\bf{x}}}_{s}$ | xs to | ${\bf{x}},$ | x, and | ${G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega )$ | Gr(xr∣x,ω) is the Green’s function from | ${\bf{x}}$ | x to the receiver position | ${{\bf{x}}}_{r}.$ | xr. Since |$ begin{eqnarray}a\,{\rm{\nabla }}\cdot {\bf{v}}=-({\rm{\nabla }}a)\cdot {\bf{v}}+{\rm{\nabla }}\cdot (a\,{\bf{v}}),\end{equation} $| a∇⋅v=−(∇a)⋅v+∇⋅(av), B.6 where | $a$ | a is a scalar field and | ${\bf{v}}$ | v is a vector field, the second equation in equation (B.5) can be written as |$ begin{eqnarray}\delta p({{\bf{x}}}_{r},\omega ) & = & \displaystyle \int -{\omega }^{2}\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}\displaystyle \frac{1}{v({\bf{x}})}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega ){p}_{0}({\bf{x}},\omega ){\rm{d}}{\bf{x}}\\ & & +\displaystyle \int \displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}}){\rm{\nabla }}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega )\cdot {\rm{\nabla }}{p}_{0}({\bf{x}},\omega ){\rm{d}}{\bf{x}}\\ & & -\displaystyle \int {\rm{\nabla }}\cdot \left(\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}}){G}_{r}\right.({{\bf{x}}}_{r}| {\bf{x}},\omega ){\rm{\nabla }}{p}_{0}({\bf{x}},\omega )){\rm{d}}{\bf{x}}.\end{equation} $| δp(xr,ω)=∫−ω2δI(x)I02(x)1v(x)Gr(xr∣x,ω)p0(x,ω)dx+∫δI(x)I02(x)v(x)∇Gr(xr∣x,ω)⋅∇p0(x,ω)dx−∫∇⋅δI(x)I02(x)v(x)Gr(xr∣x,ω)∇p0(x,ω))dx. B.7 According to Gauss’s theorem, |$ begin{eqnarray}\displaystyle \int {\rm{\nabla }}\cdot {\bf{v}}{\rm{d}}{\bf{x}}=\displaystyle {\int }_{s}{\bf{v}}\cdot {\rm{d}}{\bf{s}},\end{equation} $| ∫∇⋅vdx=∫sv⋅ds, B.8 the last integral in equation (B.7) can be rewritten as |$ begin{eqnarray}\displaystyle \int {\rm{\nabla }}\cdot \left(\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}})G{}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega ){\rm{\nabla }}{p}_{0}({\bf{x}},\omega )\right){\rm{d}}{\bf{x}}\\ \,=\displaystyle {\int }_{s}\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}}){G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega ){\rm{\nabla }}{p}_{0}({\bf{x}},\omega )\cdot {\rm{d}}{\bf{s}},\end{equation} $| ∫∇⋅δI(x)I02(x)v(x)G(xr∣x,ω)r∇p0(x,ω)dx=∫sδI(x)I02(x)v(x)Gr(xr∣x,ω)∇p0(x,ω)⋅ds, B.9 which is equal to zero. This is because we can choose a bounding integral surface with zero impedance perturbation, which is interpreted as the homogeneous boundary condition by Tarantola (1984). Hence the wavefield perturbation can be expressed as |$ begin{eqnarray}\delta p({{\bf{x}}}_{r},\omega ) & = & \displaystyle \int -{\omega }^{2}\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}\displaystyle \frac{1}{v({\bf{x}})}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega ){p}_{0}({\bf{x}},\omega ){\rm{d}}{\bf{x}}\\ & & +\displaystyle \int \displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}}){\rm{\nabla }}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega )\,\cdot \,{\rm{\nabla }}{p}_{0}({\bf{x}},\omega ){\rm{d}}{\bf{x}}.\end{equation} $| δp(xr,ω)=∫−ω2δI(x)I02(x)1v(x)Gr(xr∣x,ω)p0(x,ω)dx+∫δI(x)I02(x)v(x)∇Gr(xr∣x,ω)⋅∇p0(x,ω)dx. B.10 Plugging the first equation of (B.5) into equation (B.10) gives |$ begin{eqnarray}\delta p({{\bf{x}}}_{r},\omega ) & = & \displaystyle \int -{\omega }^{2}\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}\displaystyle \frac{1}{v({\bf{x}})}{G}_{r}\\ & & \times ({{\bf{x}}}_{r}| {\bf{x}},\omega ){G}_{s}({\bf{x}}| {{\bf{x}}}_{s},\omega )s(\omega ){\rm{d}}{\bf{x}}\\ & & +\displaystyle \int \displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}v({\bf{x}})({\rm{\nabla }}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega )\\ & & \times {\rm{\nabla }}{G}_{s}({\bf{x}}| {{\bf{x}}}_{s},\omega ))s(\omega ){\rm{d}}{\bf{x}}.\end{equation} $| δp(xr,ω)=∫−ω2δI(x)I02(x)1v(x)Gr×(xr∣x,ω)Gs(x∣xs,ω)s(ω)dx+∫δI(x)I02(x)v(x)(∇Gr(xr∣x,ω)×∇Gs(x∣xs,ω))s(ω)dx. B.11 From ray theory, the Green’s function can be approximated as |$ begin{eqnarray}G\approx A{{\rm{e}}}^{-{\rm{i}}\omega \tau }.\end{equation} $| G≈Ae−iωτ. B.12 Thus, its gradient is |$ begin{eqnarray}{\rm{\nabla }}G\approx -{\rm{i}}\omega \,{\rm{\nabla }}\tau \,G.\end{equation} $| ∇G≈−iω∇τG. B.13 The dot product in equation (B.11) becomes |$ begin{eqnarray}{\rm{\nabla }}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega )\cdot {\rm{\nabla }}{G}_{s}({\bf{x}}| {{\bf{x}}}_{s},\omega )=-{\omega }^{2}\\ \,\times ({\rm{\nabla }}{\tau }_{r}({\bf{x}})\cdot {\rm{\nabla }}{\tau }_{s}({\bf{x}})){G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega ){G}_{s}({\bf{x}}| {{\bf{x}}}_{s},\omega ),\end{equation} $| ∇Gr(xr∣x,ω)⋅∇Gs(x∣xs,ω)=−ω2×(∇τr(x)⋅∇τs(x))Gr(xr∣x,ω)Gs(x∣xs,ω), B.14 where | ${\tau }_{s}({\bf{x}})$ | τs(x) is the travel-time function from the source to | ${\bf{x}},$ | x, and | ${\tau }_{r}({\bf{x}})$ | τr(x) is the travel-time function from | ${\bf{x}}$ | x to the receiver position | ${{\bf{x}}}_{r}.$ | xr. For isotropic media, |$ begin{eqnarray}| {\rm{\nabla }}{\tau }_{s}| =| {\rm{\nabla }}{\tau }_{r}| =\displaystyle \frac{1}{v({\bf{x}})},\end{equation} $| ∣∇τs∣=∣∇τr∣=1v(x), B.15 thus, |$ begin{eqnarray}{\rm{\nabla }}{\tau }_{r}({\bf{x}})\cdot {\rm{\nabla }}{\tau }_{s}({\bf{x}})=\displaystyle \frac{1}{{v}^{2}({\bf{x}})}\,\cos \,\theta ,\end{equation} $| ∇τr(x)⋅∇τs(x)=1v2(x)cosθ, B.16 where | $\theta $ | θ is the angle between the two gradients of the travel-time functions. Substituting equations (B.14) and (B.16) into equation (B.11) yields |$ begin{eqnarray}\delta p({{\bf{x}}}_{r},\omega ,\theta ) & = & \displaystyle \int -{\omega }^{2}\displaystyle \frac{\delta I({\bf{x}})}{{I}_{0}^{2}({\bf{x}})}\displaystyle \frac{(1+\,\cos \,\theta )}{v({\bf{x}})}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega ){G}_{s}({\bf{x}}| {{\bf{x}}}_{s},\omega )s(\omega ){\rm{d}}{\bf{x}}\\ & = & \displaystyle \int -\displaystyle \frac{{\omega }^{2}(1+\,\cos \,\theta )}{{I}_{0}^{2}({\bf{x}})v({\bf{x}})}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega ){p}_{0}({\bf{x}},\omega )\delta I({\bf{x}}){\rm{d}}{\bf{x}}.\end{equation} $| δp(xr,ω,θ)=∫−ω2δI(x)I02(x)(1+cosθ)v(x)Gr(xr∣x,ω)Gs(x∣xs,ω)s(ω)dx=∫−ω2(1+cosθ)I02(x)v(x)Gr(xr∣x,ω)p0(x,ω)δI(x)dx. B.17 According to equation (B.17), the diffraction at a position | ${{\bf{x}}}_{r}$ | xr along the diffraction angle | $\theta $ | θ caused by the impedance perturbation | $\delta I$ | δI at | ${{\bf{x}}}_{0}$ | x0 can be expressed as |$ begin{eqnarray}\delta p({{\bf{x}}}_{r},\omega ,\theta )=\displaystyle \int -\displaystyle \frac{{\omega }^{2}(1+\,\cos \,\theta )}{{I}_{0}^{2}({\bf{x}})v({\bf{x}})}{G}_{r}({{\bf{x}}}_{r}| {\bf{x}},\omega ){p}_{0}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,({\bf{x}},\omega )\delta I({\bf{x}})\delta ({\bf{x}}-{{\bf{x}}}_{0}){\rm{d}}{\bf{x}}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=\,-\displaystyle \frac{{\omega }^{2}(1+\,\cos \,\theta )}{{I}_{0}^{2}({{\bf{x}}}_{0})v({{\bf{x}}}_{0})}{G}_{r}({{\bf{x}}}_{r}| {{\bf{x}}}_{0},\omega ){p}_{0}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,({{\bf{x}}}_{0},\omega )\delta I({{\bf{x}}}_{0}).\end{equation} $| δp(xr,ω,θ)=∫−ω2(1+cosθ)I02(x)v(x)Gr(xr∣x,ω)p0(x,ω)δI(x)δ(x−x0)dx=−ω2(1+cosθ)I02(x0)v(x0)Gr(xr∣x0,ω)p0(x0,ω)δI(x0). B.18 which is equivalent to equation (21). Similarly, following the same derivation, the diffraction caused by the velocity perturbation | $\delta v$ | δv can be obtained as |$ begin{eqnarray}\delta p({{\bf{x}}}_{r},\omega ,\theta )=-\displaystyle \frac{{\omega }^{2}(1-\,\cos \,\theta )}{I({{\bf{x}}}_{0}){v}_{0}^{2}({{\bf{x}}}_{0})}{G}_{r}({{\bf{x}}}_{r}| {{\bf{x}}}_{0},\omega )\\ \,\times {p}_{0}({{\bf{x}}}_{0},\omega )\delta v({{\bf{x}}}_{0}),\end{equation} $| δp(xr,ω,θ)=−ω2(1−cosθ)I(x0)v02(x0)Gr(xr∣x0,ω)×p0(x0,ω)δv(x0), B.19 which is equivalent to equation (23). References Aki K , Richards P . , 1997 , Quantitative Seismology-Theory and Method Mill Valley, CA University Science Books Bleistein N , Zhang Y , Xu S , Zhang G , Gray S H . , 2005 Migration/inversion: think image point coordinates, process in acquisition surface coordinates , Inverse Problems , vol. 21 (pg. 1715 - 1744 )1715–44 https://doi.org/10.1088/0266-5611/21/5/013 Google Scholar Crossref Search ADS Bright D , Jones C E , Belhassen Y , Brasil R , Macintyre H , França C . , 2015 Full waveform inversion for shallow hazard identification on a narrow azimuth dataset EAGE Expanded Abstracts: 77th EAGE Conf. and Exhibition https://doi.org/10.3997/2214-4609.201413540 Bunks C , Saleck F , Zaleski S , Chavent G . , 1995 Multiscale seismic waveform inversion , Geophysics , vol. 60 (pg. 1457 - 1473 )1457–73 https://doi.org/10.1190/1.1443880 Google Scholar Crossref Search ADS Červený V . , 2001 , Seismic Ray Theory Cambridge Cambridge University Press Claerbout J F . , 1992 , Earth Soundings Analysis: Processing Versus Inversion Oxford Blackwell Dai W , Fowler P , Schuster G T . , 2012 Multi-source least-squares reverse time migration , Geophys. Prospect. , vol. 60 (pg. 681 - 695 )681–95 https://doi.org/10.1111/j.1365-2478.2012.01092.x Google Scholar Crossref Search ADS da Silva N V , Ratcliffe A , Vinje V , Conroy G . , 2016 A new parameter set for anisotropic multiparameter full-waveform inversion and application to a north sea data set , Geophysics , vol. 81 (pg. U25 - 38 )U25–38 https://doi.org/10.1190/geo2015-0349.1 Google Scholar Crossref Search ADS da Silva N V , Yao G . , 2018 Wavefield reconstruction inversion with a multiplicative cost function , Inverse Problems , vol. 34 015004 https://doi.org/10.1088/1361-6420/aa9830 Fatti J , Smith G , Vail P , Strauss P , Levitt P . , 1994 Detection of gas in sandstone reservoirs using AVO analysis: a 3D seismic case history using the geostack technique , Geophysics , vol. 59 (pg. 1362 - 1376 )1362–76 https://doi.org/10.1190/1.1443695 Google Scholar Crossref Search ADS Gardner G , Gardner L , Gregory A . , 1974 Formation velocity and density—the diagnostic basics for stratigraphic traps , Geophysics , vol. 39 (pg. 770 - 780 )770–80 https://doi.org/10.1190/1.1440465 Google Scholar Crossref Search ADS Guasch L . , 2011 3D elastic full-waveform inversion , PhD Thesis Imperial College London (http://hdl.handle.net/10044/1/9974) Kaplan S T , Routh P S , Sacchi M D . , 2010 Derivation of forward and adjoint operators for least-squares shot-profile split-step migration , Geophysics , vol. 75 (pg. S225 - S235 )S225–35 https://doi.org/10.1190/1.3506146 Google Scholar Crossref Search ADS Koefoed O . , 1962 Reflection and transmission coefficients for plane longitudinal incident waves , Geophys. Prospect. , vol. 10 (pg. 304 - 351 )304–51 https://doi.org/10.1111/j.1365-2478.1962.tb02016.x Google Scholar Crossref Search ADS Kühl H , Sacchi M D . , 2003 Least-squares wave-equation migration for AVP/AVA inversion , Geophysics , vol. 68 (pg. 262 - 273 )262–73 https://doi.org/10.1190/1.1543212 Google Scholar Crossref Search ADS Lancaster S , Whitcombe D . , 2000 Fast-track ‘coloured’ inversion SEG Expanded Abstracts: 70th Annual Int. Meeting (pg. 1572 - 1575 )pp 1572–5 https://doi.org/10.1190/1.1815711 Latimer R , Davidson R , van Riel P . , 2000 An interpreter’s guide to understanding and working with seismic-derived acoustic impedance data , Leading Edge , vol. 19 (pg. 242 - 256 )242–56 https://doi.org/10.1190/1.1438580 Google Scholar Crossref Search ADS Lazaratos S , Chikichev I , Wang K . , 2011 Improving the convergence rate of full wavefield inversion using spectral shaping inversion SEG Expanded Abstracts: 81st Annual Int. Meeting (pg. 2428 - 2432 )pp 2428–32 https://doi.org/10.1190/1.3627696 Liu Y , Sun H , Chang X . , 2005 Least-squares wave-path migration , Geophys. Prospect. , vol. 53 (pg. 811 - 816 )811–6 https://doi.org/10.1111/j.1365-2478.2005.00505.x Google Scholar Crossref Search ADS Liu Y . , 2013 Globally optimal finite-difference schemes based on least squares , Geophysics , vol. 78 (pg. T113 - T132 )T113–32 https://doi.org/10.1190/geo2012-0480.1 Google Scholar Crossref Search ADS Lu R . , 2016 Revealing overburden and reservoir complexity with high-resolution FWI SEG Expanded Abstracts: 86th Annual Int. Meeting (pg. 1242 - 1246 )pp 1242–6 https://doi.org/10.1190/segam2016-13872562.1 Nemeth T , Wu C , Schuster G T . , 1999 Least-squares migration of incomplete reflection data , Geophysics , vol. 64 (pg. 208 - 221 )208–21 https://doi.org/10.1190/1.1444517 Google Scholar Crossref Search ADS Ostrander W . , 1984 Plane-wave reflection coefficients for gas sands at nonnormal angles of incidence , Geophysics , vol. 49 (pg. 1637 - 1648 )1637–48 https://doi.org/10.1190/1.1441571 Google Scholar Crossref Search ADS Plessix R É , Li Y . , 2013 Waveform acoustic impedance inversion with spectral shaping , Geophys. J. Int. , vol. 195 (pg. 301 - 314 )301–14 https://doi.org/10.1093/gji/ggt233 Google Scholar Crossref Search ADS Pratt R G , Shipp R M . , 1999 Seismic waveform inversion in the frequency domain: I. Theory and verification in a physical scale model , Geophysics , vol. 64 (pg. 888 - 901 )888–901 https://doi.org/10.1190/1.1444597 Google Scholar Crossref Search ADS Routh P , et al. , 2016 Impact of high-resolution FWI in the Western Black Sea: revealing overburden and reservoir complexity , Leading Edge , vol. 36 (pg. 60 - 66 )60–6 https://doi.org/10.1190/tle36010060.1 Google Scholar Crossref Search ADS Shuey R . , 1985 A simplification of the Zoeppritz equations , Geophysics , vol. 50 (pg. 609 - 614 )609–14 https://doi.org/10.1190/1.1441936 Google Scholar Crossref Search ADS Tarantola A . , 1984 Inversion of seismic reflection data in the acoustic approximation , Geophysics , vol. 49 (pg. 1259 - 1266 )1259–66 https://doi.org/10.1190/1.1441754 Google Scholar Crossref Search ADS Tooley R , Spencer T , Sagoci H . , 1965 Reflection and transmission of plane compressional waves , Geophysics , vol. 30 (pg. 552 - 570 )552–70 https://doi.org/10.1190/1.1439622 Google Scholar Crossref Search ADS Velzeboer C . , 1981 The theoretical seismic reflection response of sedimentary sequences , Geophysics , vol. 46 (pg. 843 - 853 )843–53 https://doi.org/10.1190/1.1441222 Google Scholar Crossref Search ADS Virieux J , Operto S . , 2009 An overview of full-waveform inversion in exploration geophysics , Geophysics , vol. 74 (pg. WCC1 - CC26 )WCC1–CC26 https://doi.org/10.1190/1.3238367 Google Scholar Crossref Search ADS Walden A T , Hosken J W J . , 1985 An investigation of the spectral properties of primary reflection coefficients , Geophys. Prospect. , vol. 33 (pg. 400 - 435 )400–35 https://doi.org/10.1111/j.1365-2478.1985.tb00443.x Google Scholar Crossref Search ADS Wang Y . , 1999 Approximations to the Zoeppritz equations and their use in AVO analysis , Geophysics , vol. 64 (pg. 1920 - 1927 )1920–7 https://doi.org/10.1190/1.1444698 Google Scholar Crossref Search ADS Wang Y H , Rao Y . , 2009 Reflection seismic waveform tomography , J. Geophys. Res. , vol. 114 B03304 https://doi.org/10.1029/2008JB005916 Wang Y H . , 2016 , Seismic Inversion: Theory and Applications Oxford Wiley-Blackwell Wang Y , Liang W , Nashed Z , Li X , Liang G , Yang C . , 2014 Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method , Geophysics , vol. 79 (pg. T277 - T285 )T277–85 https://doi.org/10.1190/geo2014-0078.1 Google Scholar Crossref Search ADS Warner M , et al. , 2013 Anisotropic 3D full-waveform inversion , Geophysics , vol. 78 (pg. R59 - 80 )R59–80 https://doi.org/10.1190/geo2012-0338.1 Google Scholar Crossref Search ADS Wu D , Yao G , Cao J , Wang Y . , 2016 Least-squares RTM with L1 norm regularisation , J. Geophys. Eng. , vol. 13 (pg. 666 - 673 )666–73 https://doi.org/10.1088/1742-2132/13/5/666 Google Scholar Crossref Search ADS Wu Z , Alkhalifah T . , 2015 Simultaneous inversion of the background velocity and the perturbation in full-waveform inversion , Geophysics , vol. 80 (pg. R317 - R329 )R317–29 https://doi.org/10.1190/geo2014-0365.1 Google Scholar Crossref Search ADS Yao G , Warner M , Silverton A . , 2014 Reflection FWI for both reflectivity and background velocity EAGE Expanded Abstracts: 76th EAGE Conf. and Exhibition https://doi.org/10.3997/2214-4609.20141089 Yao G , Warner M . , 2015 Bootstrapped waveform inversion: long-wavelength velocities from pure reflection data EAGE Expanded Abstracts: 77th EAGE Conf. and Exhibition https://doi.org/10.3997/2214-4609.201413208 Yao G , Wu D . , 2015 Least-squares reverse-time migration for reflectivity imaging , Sci. China Earth Sci. , vol. 58 (pg. 1982 - 1992 )1982–92 https://doi.org/10.1007/s11430-015-5143-1 Google Scholar Crossref Search ADS Yao G , Jakubowicz H . , 2016 Least-squares reverse-time migration in a matrix-based formulation , Geophys. Prospect. , vol. 64 (pg. 611 - 621 )611–21 https://doi.org/10.1111/1365-2478.12305 Google Scholar Crossref Search ADS Yao G , Wu D , Debens H A . , 2016 Adaptive finite difference for seismic wavefield modelling in acoustic media , Sci. Rep. , vol. 6 pg. 30302 https://doi.org/10.1038/srep30302 Google Scholar Crossref Search ADS PubMed Yao G , Wu D . , 2017 Reflection full waveform inversion , Sci. China Earth Sci. , vol. 60 (pg. 1783 - 1794 )1783–94 https://doi.org/10.1007/s11430-016-9091-9 Google Scholar Crossref Search ADS Yao G , da Silva N V , Debens H A , Wu D . , 2017a Accurate seabed modeling using finite difference methods , Comput. Geosci. https://doi.org/10.1007/s10596-017-9705-5 Yao G , da Silva N V , Wu D . , 2017b Forward modelling formulas for least-squares reverse-time migration , Explor. Geophys. https://doi.org/10.1071/EG16157 Yilmaz O . , 2001 , Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data Tulsa, OK Society of Exploration Geophysicists Zhang F , Li X . , 2015 Exact elastic impedance tensor in isotropic media , Sci. China Earth Sci. , vol. 58 (pg. 1350 - 1360 )1350–60 https://doi.org/10.1007/s11430-015-5079-5 Google Scholar Crossref Search ADS Zhang F , Li X . , 2016 Exact elastic impedance matrices for transversely isotropic media , Geophysics , vol. 81 (pg. 1350 - 1360 )1350–60 https://doi.org/10.1190/geo2015-0163.1 Zhang J , Yao Z . , 2012 Optimized finite-difference operator for broadband seismic wave modeling , Geophysics , vol. 78 (pg. A13 - A18 )A13–8 https://doi.org/10.1190/geo2012-0277.1 Google Scholar Crossref Search ADS Zhang Y , Duan L , Xie Y . , 2015 A stable and practical implementation of least-squares reverse time migration , Geophysics , vol. 80 (pg. V23 - 31 )V23–31 https://doi.org/10.1190/geo2013-0461.1 Google Scholar Crossref Search ADS Zhou W , Brossier R , Operto S , Virieux J . , 2015 Full waveform inversion of diving & reflected waves for velocity model building with impedance inversion based on scale separation , Geophys. J. Int. , vol. 202 (pg. 1535 - 1554 )1535–54 https://doi.org/10.1093/gji/ggv228 Google Scholar Crossref Search ADS © 2018 Sinopec Geophysical Research Institute
TI - Sensitivity analyses of acoustic impedance inversion with full-waveform inversion
JO - Journal of Geophysics and Engineering
DO - 10.1088/1742-2140/aaa980
DA - 2018-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/sensitivity-analyses-of-acoustic-impedance-inversion-with-full-CccJubAmpC
SP - 461
VL - 15
IS - 2
DP - DeepDyve
ER -