Chapter 2
One-Return Propagators and the Applications in
Modeling and Imaging
Ru-Shan Wu, Xiao-Bi Xie, and Shengwen Jin
This chapter reviews the one-return method for calculating seismic wave propagation in complex acoustic and elastic models. We give a new, intuitive derivation of the one-return approximation using sequential thin-slab transmission/reflection operators. This deriva-tion can reach the same formuladeriva-tion as the De Wolf approximaderiva-tion, and in the same time provide an efficient implementation of the method. The method is based on the multiple-forescattering-single-backscattering approximation. It neglects the internal reverberations (internal multiples) but can take into account all forward scattering phenomena, such as focusing/defocusing, diffraction, refraction, interference, and model primary reflections from heterogeneities. One-return method can be implemented using an iterative march-ing algorithm shuttlmarch-ing between the space and wavenumber domains. For models where reverberation and resonance scattering can be neglected, this method provides an accurate and highly efficient algorithm, especially for large velocity models at high-frequencies. It has been used for reservoir AVO simulations in heterogeneous visco-elastic media with moderate elastic parameter variations (perturbations less than 40%). In seismic imaging, one-return propagators can migrate turning waves and duplex waves (or prism waves) and therefore, image steeply dipping structures and vertical overhangs. For waveform inversion and migration velocity analysis, the one-return propagator is very efficient for calculating sensitivity kernels that involve reflections from interfaces. We also reviewed recent progress in overcoming the limitations and shortcomings of the way, one-return methods. The one-way propagator with super-wide angle interpolates the wave-fronts from two orthogonal one-way propagators to reconstruct the accurate wavefront up to nearly the full angle range, and was used to image overhanging salt flanks. One-return boundary element method was developed to deal with strong-contrast elastic media. It can generate synthetics with only primary reflections from the salt boundaries, and there-fore, be used to study the artifacts caused by salt multiples in subsalt images. However, the current version uses Green’s functions in homogeneous media. One-return modeling in heterogeneous, strong-contrast media is still a challenging problem.
Keywords: Elastic wave modeling, One way propagator, One-return approximation,
2.1
Introduction
One-way approximation for wave propagation has been introduced and widely used as propagators in forward and inverse problems of scalar, acoustic and elastic waves (e.g., Claerbout, 1970, 1976; Landers and Claerbout, 1972; Flatt´e and Tappert, 1975; Corones, 1975; Tappert, 1977; McCoy, 1977; Hudson, 1980; Ma, 1982; Wales and McCoy, 1983; Fishman and McCoy, 1984, 1985; Wales, 1986; McCoy and Frazer, 1986; Collins, 1989, 1993; Collins and Westwood, 1991; Stoffa et al., 1990; Fisk and McCarter, 1991; Wu and Huang, 1992; Zhang, 1993; Ristow and Ruhl, 1994; Wu, 1994, 1996, 2003; Wu and Xie, 1994; Wu and Jin, 1997; Xie and Wu, 1998, 2001, 2005;Grimbergen et al., 1998; Van Stralen et al., 1998; Wild and Hudson, 1998; Jin et al., 1998, 1999, 2002; Jin and Wu, 1999; Huang et al., 1999a,b; Thomson, 1999, 2005; De Hoop et al., 2000; Lee at al., 2000; Wild et al., 2000; Wu et al., 2000a,b; Le Rousseau and de Hoop, 2001; Wu and Wu, 2001; Han and Wu, 2005; Zhang et al., 2005, 2006). The great advantages of one-way propagation methods are the fast speed of computation, often being several orders of magnitudes faster than the full-wave finite difference and finite element methods, and the huge saving in internal memory.
The successful extension and applications of one-way elastic wave propagation meth-ods has stimulated the research interest in developing similar theory and techniques for reflected or backscattered wave calculation. There are several approaches in extending the one-way propagation method to include backscattering and multiple scattering cal-culations, such as the generalized Bremmer series (GBS) approach (Corones, 1975; De Hoop, 1996; Wapenaar, 1996, 1998; van Stralen et al., 1998; Thomson, 1999, 2005; Le Rousseau and de Hoop, 2001) and the generalized screen propagator (GSP) approach (Wu, 1994, 1996, 2003; Wu and Xie, 1993, 1994; Wu et al., 1995; Wild and Hudson, 1998; de Hoop et al., 2000; Wild et al., 2000; Xie et al., 2000; Xie and Wu, 2001, 2005; Wu et al., 2007). The key difference between these approaches is how to define a ref-erence Green’s function in the multiple scattering series. The GBS approach adopts an asymptotic solution of the acoustic or elastic wave equation in the heterogeneous medium as the Green’s function, i.e. the one-way propagator in the preferred direction. The mul-tiple scattering series is based on the interaction of Green’s field (incident field) with the medium heterogeneities. The GSP approach, on the other hand, does not use asymptotic solutions. Instead, the approach uses the multiple-forward-scattering (MFS) corrected one-way propagator as the Green’s function. When the backscattered field, calculated at each thin-slab, is propagated in the backward direction, it uses the same MFS corrected one-way propagator. In Wu et al. (2007), detailed analysis and comparison of different approximations can be found. In this chapter, we will focus on the GSP approach based on the De Wolf approximation (one-return approximation).
In the generalized screen approach, Wu and Huang (1995) introduced a wide-angle modeling method for backscattered acoustic waves using the multiple-forward-scattering approximation and a phase-screen propagator. Xie and Wu (1995; 2001) extended the complex screen method to include the calculation of backscattered elastic waves under the small-angle approximation. Wu (1996) and Wu et al., (2007) derived a more gen-eral theory for acoustic and elastic waves using the De Wolf approximation, and the
the-ory provided two versions of algorithms: the thin-slab method and the complex-screen method. Wu and Wu (1999, 2003a) introduced a fast implementation of the thin-slab method and a second order improvement for the complex-screen method. Wu and Xie (2009) discussed some practical applications of this method. In the rest part of this chap-ter, we review the one-return method and discuss its applications in seismic modeling, imaging and inversion with numerical examples.
2.2
Primary-Only Modeling and One-Return Approximation
“Primaries” are referred to as the primary transmitted or reflected waves, which involve only single interaction at each interface (transmission or reflection) or at each small-scale volumetric heterogeneity (forward or backscattering) inside a complex medium. The con-cept of primary transmission and reflection is closely related to the acquisition (observa-tion) geometry. Mathematically, primary transmission and reflection (T/R) can be defined with respect to certain preferred direction. In surface seismic reflection survey, z-direction is conveniently defined as the preferred direction (globally). In smoothly inhomogeneous media, the preferred direction can be also defined as being along the ray-direction. In this chapter, we will discuss only the globally preferred direction (z-direction) in wave-equation-based methods.
The mathematical basis of primary-only modeling is the multiple-forescattering-single-backscattering (MFSB) approximation, which is also called one-way-one-return (OWOR) approximation. The physics behind this approximation is the neglect of multiple scatter-ing between the forward and backward directions (with respect to the preferred direc-tion), i.e. the reverberations. The OWOR is totally wave-equation-based, not the high-frequency asymptotics. This approach needs the separation of the forward-propagated and backward-propagated waves in the wave equation. The separation can be done by either decomposing the scattered wavefield or factorizing the scattering operator into two parts: the forward part and the backward part. The wavefield decomposition leads to a marching wave equation method which may include the coupling between the forward and back-ward waves (e.g., Weston, 1989), while the operator decomposition (factorization) forms the theoretical basis of way or return modeling. In the latter approach, the one-way-one-return approximation can be formulated as the first order term in a reordered multiple scattering series (the De Wolf series). In this way, the multiple backscattering and reverberations can be calculated by multiple sweeps of the one-return algorithm. In this chapter we limit the treatment to the approach of operator decomposition.
For operator decomposition, the formulation can be based on either PDE (partial differ-ential equation) or the integral equation. For the PDE approach, the operator factorization results in a square-root operator and many approximations for one-way propagators were derived based on the PsDO (pseudo-differential operator) theory (Fishman and McCoy, 1984, 1985; Zhang, 1993; Thomson, 1999, 2005; De Hoop et al., 2000; Zhang et al., 2005, 2006, 2007). The wave equation solution after the factorization can be casted into generalized Bremmer series (De Hoop, 1996; De Hoop et al., 2000). For the integral equation approach, the approximation can be derived from the De Wolf series and its first
approximation (De Wolf approximation or one-return approximation) (De Wolf, 1971, 1985; Wu, 1994, 1996, 2003; Xie and Wu, 2001; Wu et al., 2007). In this chapter, we will follow the integral equation approach and the De Wolf approximation.
Integral equation, and differential equation are different forms of wave equation, and describe mathematically the same physical phenomena: wave motion and propagation. Integral wave equation, such as the Lippmann-Schwinger equation, can be derived from the differential wave equation using the Green’s function. The solution of the Lippmann-Schwinger equation can be casted into a Born series, of which the first order term is called the Born approximation. The Born approximation is among the simplest solutions and widely used in wave scattering problems. The approximation is a weak scattering approximation, and only valid when the scattered field is much smaller than the incident field. This implies that the heterogeneities are weak and the propagation distance is short (or the scattering volume is small). The valid region of the Born approximation for the forward scattering is very different from that for the backscattering. Scattered fields in the Born approximation are simple summation of contributions from all parts of the scat-tering volume, and scatscat-tering at each part is independent of other parts. Therefore, the Born approximation does not obey energy conservation, In the forward direction, scat-tered fields from different sections along the propagation path arrive at the receiving point in phase with the incident field, so that they will be coherently superposed. This leads to the linear increase of the total field, which is a serious divergence problem for forward scattering. However, the total observed field in the backward direction is the sum of all the backscattered fields from all the scatterers since there is no incident wave in the backward direction. In addition, the volume of the scatterers of coherent stacking for backscattered waves is limited due to the travel time differences. For this reason, backscattering does not have the catastrophic divergence and the validity region for the Born approximation is much larger than the case of forward scattering. The De Wolf approximation, which is a multiple-forescattering-single-backscattering (MFSB) approximation, has been in-troduced to overcome the limitation of the Born approximations in long-range forward propagation.
The mathematical basis for the Born series and De Wolf series as the solution of the Lippmann-Schwinger equation (integral wave equation) has been detailed in Wu et al. (2007). Here we adopt a more intuitive approach to show the physics and calculation procedure of the method. Despite of starting from the physical intuition, we can reach the same mathematical formulation if we make the limiting process using the corresponding notations.
For an arbitrary heterogeneous media of a large volume, as shown in Figure 2.1a, we can slice the volume into numerous thin-slabs transversal to the propagation direction (preferred direction). Shown in Figure 2.1b and Figure 2.1c are examples of individual thin-slabs for a scalar medium (velocity is the only parameter), and an elastic medium. Assume each thin-slab is thin enough so that the Born approximation can be used for the scattering calculation. In Figure 2.1a, u0is the incident field, U is the scattered field
and the total field will be denoted by u= u0+ U. The one-return approximation uses the
scattering operator spitting, so the scattered waves by the thin-slab are divided into the forescattered wave Uf and the backscattered wave Ub. At the thin-slab exit, we add the
Fig. 2.1 (a) Cartoon showing the thin-slab decomposition of a velocity model, where a heterogeneous medium
is sliced into an N-stack of thin-slabs; (b) An acoustic (scalar wave) thin-slab, where the incident plus the forescattered waves form the transmitted wave and the back-scattered wave composes the reflected wave; (c) An elastic thin-slab and the different scattered waves generated due to the interaction between incident P and S waves and the thin-slab. The superscripts PP, PS, SS and SP denote the P-to-P, P-to-S, S-to-S and S-to-P scatterings.
new total forward field will be used as the updated incident field for the next thin-slab. For an elastic thin slab shown in Figure 2.1c, the superscripts denote different wave types (P or S) and conversion modes (P-to-P, P-to-S, S-to-S and S-to-P scatterings). Similar to the scalar wave case, after interacting with an elastic thin slab, all forward scattered waves together with the incident waves form the forward P and S waves at the exit side of the slab (the incident waves for the next thin-slab), while at the entrance side, backscattered fields form the reflected P and S waves. We first discuss the scalar wave case to demonstrate the principle of one-return approximation. In the next section, we will summarize the formulation for the elastic wave case.
Figure 2.2 shows a cartoon illustrating the marching algorithm of the one-return ap-proximation. For the lth thin-slab, we define a slab transmission operator Tl and a slab
reflection operator (backscattering operator) Rl. Together they form the split operators in
the forward and backward directions. The transmission operator can be expressed explic-itly by the forescattering operator Fl,
Fig. 2.2 Cartoon illustrating the double-sweep procedure of the one-return modeling. In the downgoing sweep,
after interacting with each slab, the transmitted P- and S-wavefields are used as the input for the next thin-slab. This process updates the incident wavefields slab-by-thin-slab. During the downgoing sweep, all backscattered wavefields at each thin-slab are stored in the memory. In the upgoing sweep, those stored reflections are retrieved and propagated to the surface using the sequential transmission operators similar to the downgoing sweep.
With the split thin-slab operators, we can write out the forescattering updated incident field uf (the transmitted wavefield) and the backscattered field Ubat each slab entrance as
(see Figure 2.2) uf(z0) = u0 Ub(z0) = R1uf(z0) uf(z1) = T1uf(z0) Ub(z1) = R2uf(z1) . . . uf(zi) = Tiuf(zi−1) = i ∏ l=1 Tl u0, Ub(zi) = Ri+1uf(zi) = Ri+1 i ∏ l=1 Tl u0 . . . uf(zN) = TNuf(zN−1), Ub(zN−1) = RNuf(zN−1) (2.2)
propagator as a sequential thin-slab transmission operator, P(zj, zi) = j
∏
l=i Tl, (2.3)we can write the transmitted wave at the bottom of the medium as
ut(zN) = P(zN, z0)u0(z0) (2.4)
and the primary reflected (backscattered) field at the top of the medium as the summation of the single backscattered waves of all the thin-slabs, propagated to the top by one-way propagators defined in (2.3). ur= N−1 ∑ i=0 P(z0, zi)Ub(zi) = N−1 ∑ i=0 P(z0, zi)Ri+1uf(zi) =N∑−1 i=0 P(z0, zi)Ri+1P(zi, z0)u0(z0). (2.5)
In the case of local Born approximation, thin-slab transmission operator is Tl= Gl1 +χf
(2.6) where Glis the local Green’s operator (in the background media) in the l-th thin-slab, “1”
is the identity operator, andχf is scattering potential to the forward half-space. Applying
to an incident field u0, the transmitted field is
uf(zl) = Tlu0(zl−1) = Gl1 +χf u0(zl−1) = u0(zl) + k2 Z thin-slab d3r′Gl x, zl; r′ χ f r′ u0 r′ . (2.7)
The thin-slab reflection operator is
Rl= Glχb (2.8)
whereχbis the scattering potential in the backward half-space. The one-way propagator
defined in (2.3) is the multiple forescattering Green’s function
Gf(zi, z0) = i
∏
l=1 Tl= i∏
l=1 Gl1 +χf(zl) . (2.9)After expanding the product and reordering the series, we reach the De Wolf approxi-mation (one-return approxiapproxi-mation) (see equation 28 of Wu et al., 2007),
Gf(zi, z0) = i
∑
m=0
where G0is the free-space (background) Green’s function between any two thin-slabs G0(zj, zi) = j
∏
l=i Gl(zl).In the same way, we can derive the updated incident field in terms of multiple forward scattering series uf(zi) = i
∑
m=0 [G0χf]mu0 0 6 m 6 i. (2.11)When taking the limit of infinitely thin-slab, i.e.∆z→ 0, we reach the integral form of one-return approximation (De Wolf approximation),
u(x, z) = uf(x, z) + k2
Z
V
d3r′Gf(x, z; r′)χb(r′)uf(r′). (2.12)
From the above derivation, we see that iterative application of thin-slab propagators can derive the De Wolf approximation (MFSB approximation). Furthermore, it provides a very efficient implementation of the MFSB (multiple-forescattering-single-backscattering). In the next section, we will concentrate on the elastic one-return modeling.
2.3
Elastic One-Return Modeling
As seen from the previous section, once we derive the thin-slab transmission and reflec-tion operators, one-return modeling (De Wolf approximareflec-tion) can be implemented by a double-sweep of marching algorithm with iterative application of thin-slab operator. In the following, we summarize the formulations of elastic thin-slab operator. The method belongs to a more general approach of generalized screen propagators (GSP).
In a linear general heterogeneous medium, the elastic dynamic equation can be ex-pressed as (Aki and Richards, 1980)
−ω2ρ(x)u(x) =∇·σ(x), (2.13) where x is the location,ω the frequency, u the displacement vector,σ the stress tensor (dyadic),∇· the divergence operator andρ the density of the medium, and there is no body force in the medium. The stress-displacement relation is
σ(x) = c(x) :ε(x) =1
2c :(∇u+ u∇), (2.14) where c is the elastic constant tensor of the medium,εis the strain field, u∇stands for the transpose of∇u (the gradient operator) and “ : ” stands for double scalar product
of tensors defined through(ab):(cd) = (b · c)(a · d). Equation (2.13) can be written as a wave equation of the displacement field:
−ω2ρ(x)u(x) =∇· 1
2c :(∇u+ u∇)
Using perturbation theory, the elastic parameters and the total wavefield can be decom-posed as
ρ(x) =ρ0+δρ(x), c(x) = c0+δc(x),
u(x) = u0(x) + U(x), (2.16)
whereρ0and c0are the parameters of the background medium,δρ andδc are the
cor-responding perturbations, u0is the incident field and U is the scattered field, then (2.15) can be rewritten as a wave equation about U
−ω2ρ0U−∇· 1 2c0:(∇U+ U∇) = F (2.17) with F=ω2δρu+∇· [δc :ε] (2.18) as an equivalent body force due to scattering.
We choose z-axis as the primary propagation direction and use u= uf, the updated
forward propagated field, as the incident wave for the thin-slab in the current marching step. Substituting the equivalent body force (2.18) and equation (2.16) into equation (2.17), we have the scattered wavefield generated by the interaction between the incident wave and the heterogeneity,
U(xT, z∗) = Z V d3x′δρω2uf(x′T, z) +∇· δc :εf(x′T, z) · G0 xT, z∗; x′T, z , (2.19)
where z∗is the observation depth, G0is the Green’s function in the background model,
x= (xT, z), with xT and z as horizontal and vertical coordinates.
2.3.1
Local Born Approximation
Following the derivations of Wu et al. (2007), we can express the scattered displace-ment field for a thin-slab in the horizontal wavenumber domain using the local Born approximation. Note that the local incident field uf andεf have the footprints of the
heterogeneities in the previous thin-slabs, so the formulation is not a global Born approx-imation. However, we assume the thin-slab is thin enough so the local incident field will not be influenced by the local heterogeneities in the current thin-slab. In this way, we can apply the Born approximation locally within the thin-slab. After Fourier transform along the transversal plane with respect to xTin equation (2.19), we obtain the local Born
expression in the horizontal wavenumber domain:
U(KT, z∗) = Z z1 z′ dz Z Z d2xT δρω2uf(x T, z) +∇· δc :εf(xT, z · G0(KT, z∗; xT, z) , (2.20)
where the thin slab is located between z′and z1, KT is the horizontal wavenumber vector, and G0(z∗, kT; z, xT) = ik2α 2ρ0ω2 ˆkαˆkα 1 γα eikα·r+ ik 2 β 2ρ0ω2 I− ˆkβˆkβ 1 γβeikβ·r, (2.21)
with I as the unit dyadic, and
γα=
q k2
α− KT2,
γβ =qkβ2− KT2, (2.22)
where kα=ω/α0and kβ =ω/β0are P and S wavenumbers withα0andβ0as the P and
S wave background velocities in the thin-slab, respectively, ˆkα = 1 kα(KT,γα) , (2.23) ˆkβ = 1 kβ KT,γβ , (2.24)
and are unit vectors along P and S wave propagation directions. For isotropic media,
δc(x) :ε(x) =δλ(x) |ε|I + 2δ µ(x)ε(x). (2.25) Substituting (2.21) into (2.20), we can derive the dual-domain expressions for scattered displacement fields in isotropic elastic media.
For P to P scattering: UPP(KT, z∗) = ik2α 2γα Z z1 z′ dze ikαz(z∗−z)· ˆkαˆkα·ZZ d2xTe−iKT·xTδρ(xT, z) ρ uαf(xT, z) −ˆkα ZZ d2xTe−iKT·xTδλ (xT, z) λ+ 2µ 1 ikα∇· u f α(xT, z) −ˆkα ˆkαˆkα : ZZ d2xTe−iKT·xTδ µ (xT, z) λ+ 2µ 1 ikαε f α(xT, z) , (2.26)
with kαz = +γαfor forescattering and kαz = −γαfor backscattering, and ˆkα= (KT, kαz)/kα.
Note that we replacedρ0,λ0, µ0 in denominators byρ =ρ0+δρ, λ =λ0+δλ and
µ=µ0+δ µ. This replacement is the result of asymptotic matching between the Born
approximation for large-angle scattering and the h-f asymptotic travel-time (phase) for forward propagation. It is proved (Wu and Wu, 2003a) that with this replacement (asymp-totic matching), the phase-shift in the exact forward direction is accurate and the phase error for small angles is reduced compared with the Born approximation. In the mean-while, the phase error for large angle scattering is much smaller than that of the phase screen approximation.
In (2.26) uαf(xT, z), ∇· uαf(xT, z) andεαf(xT, z) can be calculated by a dual domain
method: uαf(xT, z) = 1 4π ZZ d2K′TeiK′T·xTu0 α(K′T)eiγ ′ α(z−z′) (2.27)
1 ikα∇· u f α(xT, z) = 1 4π ZZ d2K′TeiK′T·xTˆk′ α· u0α(K′T)eiγ ′ α(z−z′) , (2.28) 1 ikαε f α(xT, z) = 1 4π ZZ d2K′TeiK′T·xT1 2[ˆk ′ αu0α(K′T) + ˆkα′U0α(K′T)]eiγ ′ α(z−z′) = 1 4π ZZ d2K′TeiK′T·xTˆk′ αˆkα′u0α(K′T)eiγ ′ α(z−z′) , (2.29) where u0α(K′ T) = u0α(K′T) and ˆkα′ = (K′T,γ′α)/kα. For P to S scattering: UPS(KT, z∗) = ik2β 2γβ Z z1 z′ dze ikβz(z∗−z) I− ˆkβˆkβ · ZZ d2xTe−iKT·xTδρ (xT, z) ρ uαf(xT, z) − I − ˆkβˆkβ · ˆkβ·ZZ d2xTe−iKT·xT2δ µ (xT, z) µ 1 ikβε f α(xT, z) , (2.30) where ˆkβ= (KT, kβz)/kβ. For S to P scattering: USP(KT, z∗) = ikα2 2γα Z z1 z′ dze ikαz(z∗−z) ˆkαˆkα·ZZ d2xTe−iKT·xTδρ (xT, z) ρ uβf(xT, z) − kkα β ˆkα ˆkαˆkα : ZZ d2xTe−iKT·xT2δ µ(xT, z) µ 1 ikβε f β(xT, z) . (2.31) For S to S scattering: USS(KT, z∗) = ik2β 2γβ Z z1 z′ dzeikβz(z∗−z) (I − ˆkβˆkβ) · ZZ d2xTe−iKT·xTδρ (xT, z) ρ uβf(xT, z) − I − ˆkβˆkβ · ˆkβ·ZZ d2xTe−iKT·xT2δ µ (xT, z) µ 1 ikβε f β(xT, z) . (2.32) In equations (2.31) and (2.32), uβf(xT, z) andεβf (xT, z) can be calculated by
uβf(xT, z) = 1 4π2 ZZ d2K′TeiK′T·xTu0 β(K′T)e iγ′β(z−z′) (2.33) 1 ikβε f β(xT, z) = 1 4π2 ZZ d2K′TeiK′T·xT1 2 h ˆk′ βu0β(K′T) + u0β(K′T)ˆkβ′ i eiγ′β(z−z′), (2.34) where ˆk′β= (K′ T,γ′β)/kβ.
2.3.2
The Thin Slab Approximation
From equations (2.26) to (2.34), we see that the leading-order interaction between incident fields and heterogeneities is expressed in three-dimensional volume integrals. Also the
scattered and incident wavenumbers are coupled with each other. So the computation of these equations is still intensive. In this section, parts of the integration over z in the equations are analytically estimated. Assume that the slab for each marching step is thin enough that the parameters (velocity and density) can be approximately taken as invariant along z, the integration with respect to z in equation (2.26) can be calculated as
Z z1 z′ dze ikαz(z∗−z)+iγ′α(z−z′)= ∆zei(γα+γ′α)∆z/2sinchγα−γ′α 2 ∆z i for forescattering(z∗= z1), ∆zei(γα+γ′α)∆z/2sinchγα+γ′α 2 ∆z i for backscattering(z∗= z′). (2.35) We see that the integration over z has been done analytically; however,γα andγ′αare still coupled, which prevents the fast computation of the thin-slab method. To decouple
γαandγ′α, we neglect the angular variation of amplitude factors but keep the phase infor-mation untouched by taking the approxiinfor-mationγα=γ′α= kαfor the amplitude factors in equation (2.35). This assumption is valid for the case where the small-angle scattering is dominant, and therefore the direction of the scattered waves are not far from the incident direction. Under this approximation, equation (2.35) becomes
Z z1 z′ dzeikαz(z∗−z)+iγ′α(z−z′)≈ ∆zei(γα+γ′α)∆z/2 for forescattering(z∗= z1), ∆zei(γα+γ′α)∆z/2sinc(kα∆z) for backscattering(z∗= z′).
(2.36) For the scattered fields P-S, S-P and S-S, similar approximations can be obtained as follows. For P-S or S-P scattering,
Z z1 z′ dze ikβz(z∗−z)+iγ′α(z−z′) ≈
∆zei(γ′α+γβ)∆z/2sinc[(kα− kβ)∆z/2] for forescattering(z∗= z1), ∆zei(γ′α+γβ)∆z/2sinc[(kα+ kβ)∆z/2] for backscattering(z∗= z′).
(2.37) For S-S scattering, Z z1 z′ dzeik β z(z∗−z)+iγ′β(z−z′) ≈ ∆zei(γβ+γ′β)∆z/2 for forescattering(z∗= z1), ∆zei(γβ+γ′β)∆z/2sinc(kβ∆z) for backscattering(z∗= z′). (2.38) After integration over z, the integration over transverse plane xT in equations (2.26)
– (2.34) can be carried out by the FFT. In order to further expedite the computation, we can group the scattered field equations (2.26) to (2.34) into UP(KT, z∗) = UPP(KT, z∗) +
USP(KT, z∗), and US(KT, z∗) = UPS(KT, z∗) + USS(KT, z∗) , i.e. UP(KT, z∗) = ikα2 2γαe iγα∆z/2∆zˆk α ˆkα· ZZ d2xTe−iKT·xTδρ (xT) ρ h ηPPuf α(xT)+ηSPuβf(xT) i − ZZ d2xTe−iKT·xTδλ (xT) λ+ 2µ 1 ikα∇· h ηPPuf α(xT) i − ˆkαˆkα : ZZ d2xTe−iKT·xT 2δ µ(xT) λ+ 2µ 1 ikα h ηPPεf α(xT) +ηSPεf β(xT) i , (2.39) US(KT, z∗) = ik2β 2γβe iγβ∆z/2∆z(I − ˆk βˆkβ) · ZZ d2xTe−iKT·xTδρ (xT) ρ h ηPSuf α(xT) +ηSSuf β(xT) i − ˆkβ· ZZ d2xTe−iKT·xT 2δ µ(xT) µ 1 ikβ h ηPSεf α(xT) +ηSSεβf(xT) i , (2.40)
where z∗(z∗= z′or z∗= z1) indicates the position of the receiver plane. The modulation
factorsηPP,ηSP,ηPSandηSSare
ηPP=
1 for forescattering sin c(kα∆z) for backscattering
, (2.41)
ηPS=ηSP=
sin c[(kα− kβ)∆z/2] for forescattering sin c[(kα+ kβ)∆z/2] for backscattering
, (2.42) ηSS= 1 for forescattering sin c(kβ∆z) for backscattering .
(2.43)
Note that the factors eiγα(z−z′)and eiγβ(z−z′)have been replaced by eiγα∆z/2and eiγβ∆z/2for
calculating the background fields. The phase matching (asymptotic matching) has been applied in equations (2.39) and (2.40).
2.3.3
Small-Angle Approximation and the Screen Propagator
In the thin-slab method (2.39) and (2.40), we need to calculate the incident displacement and strain fields in each marching step. We can further simplify the calculation by using the complex-screen approximation which is based on the small angle approximation. The approximation needs to be made in the wavenumber domain. The wavenumber domain
formulation can be obtained by substituting equations (2.29), (2.34) into equations (2.26) and (2.30)-(2.32). UPP(KT, K′T) = i 2γαk 2 αuP0ˆkα ˆkα· ˆkα′δρρ( ˜k)−λδλ+ 2( ˜k)µ− (ˆkα· ˆkα′)22λδ µ+ 2( ˜k)µ , (2.44) UPS(KT, K′T) = i 2γβk 2 β I− ˆkβˆkβ · uP0 δρ( ˜k) ρ −2 β0 α0 ˆkα· ˆk′ β δ µ( ˜k) µ , (2.45) USP(KT, K′T) = i 2γαk 2 α uS0· ˆkα ˆk α δρ( ˜k) ρ − 2 β0 α0 ˆkβ· ˆk′αδ µµ( ˜k) , (2.46) USS(KT, K′T) = i 2γβk 2 β I− ˆkβˆkβ · uS0 δρ ( ˜k) ρ − (ˆkβ· ˆkβ′)δ µµ( ˜k) −ˆkβ′ uS0· ˆkβ δ µ( ˜k) µ , (2.47)
where uP0 is the spectral field of the incident P wave, andδρ( ˜k),δλ( ˜k) andδ µ( ˜k) are the three-dimensional Fourier transforms of medium perturbations, and ˜k= ˆk − ˆk′is the exchange wavenumber with ˆk′ and ˆk as incident and scattering wavenumber vectors, re-spectively.
Under small angle approximation, both incoming and outgoing wavenumbers have small transversal components KT compared to the longitudinal component γ (vertical
wavenumber). For an isotropic medium, the scattered fields can be expressed as (Xie and Wu, 2001) UPPf (KT, z0) = −ikα∆zˆkαηPPf ZZ dx′Te−iKT·x′TuP 0(x′T, z0) δα(x′ T) α0 , (2.48) UPSf (KT, z0) = −ikβ∆zηPSf ˆkβ× ˆkβ×ZZ dx′Te−iKT·x′TuP 0(x′T, z0) β 0 α0− 1 2 δρ (x′T) ρ0 +2 β 0 α0 δβ(x′ T) β0 , (2.49) USPf (KT, z0) = −ikα∆zηSPf ˆkα ˆkα·ZZ dx′Te−iKT·x′TuS 0(x′T, z0) β 0 α0− 1 2 δρ(x′ T) ρ0 +2 β 0 α0 δβ (x′T) β0 , (2.50) USSf (KT, z0) = −ikβ∆zηSSf ˆkβ× ˆkβ× ZZ dx′Te−iKT·x′TuS 0(x′T, z0)δβ (x′ T) β0 , (2.51)
UPPb (KT, z0) = −ikα∆zˆkαηbPP ZZ dx′Te−iKT·x′TuP 0(x′T, z0) δZα(x′T) Zα0 , (2.52) UPSb (KT, z0) = ikβ∆zηbPSˆkβ× ˆkβ× ZZ dx′Te−iKT·x′TuP 0(x′T, z0) β 0 α0 +1 2 δρ(x′ T) ρ0 +2 β 0 α0 δβ (x′T) β0 , (2.53) USPb (KT, z0) = ikα∆zηbSPˆkα ˆkα·ZZ dx′Te−iKT·x′TuS 0(x′T, z0) β 0 α0 +1 2 δρ (x′T) ρ0 +2 β 0 α0 δβ(x′ T) β0 , (2.54) USSb (KT, z0) = ikβ∆zηbSSˆkβ× ˆkβ×ZZ dx′Te−iKT·x′TuS 0(x′T, z0) δZβ(x′ T) Zβ0 , (2.55)
where Zα=ραand Zβ =ρβ are P and S wave impedances and Zα0=ρ0α0and Zβ0=
ρ0β0are their background values.ηPP,ηPS,ηSPandηSScan be calculated with equations
(2.41)-(2.43). The subscripts f and b denote the forward and backward scattered waves, respectively. From equations (2.48)-(2.55), it can be seen that transmitted waves UPPf and
USSf are controlled by P- and S-wave velocity perturbations and the reflected waves UPPb and USS
b are controlled by impedance perturbations.
In the forward direction, there are both incidence waves and scattered waves. The forward propagated wavefield uf can be obtained as follows:
uf(xT, z1) = 1 4π2 Z dK′T[uPf(K′T, z1) + uSf(K′T, z1)]eiK ′ T·xT, (2.56) where uPf(K′T, z1) = eiγ ′ α|z1−z0|uP 0(K′T, z0) + UPPf (K′T, z0) + USPf (K′T, z0) , (2.57) uSf(K′T, z1) = e iγβ′|z1−z0|uS 0(K′T, z0) + USSf (K′T, z0) + UPSf (K′T, z0) . (2.58)
The backward propagated wavefield ubis composed of pure scattered waves and can
be obtained as follows: ub(xT, z1) = 1 4π2 Z dK′T[uPb(K′T, z1) + uSb(K′T, z1)]eiK ′ T·xT, (2.59) where uPb(K′T, z1) =UPPb (K′T, z0) + USPb (K′T, z0) , (2.60) uSb(K′T, z1) =USSb (K′T, z0) + UPSb (K′T, z0) . (2.61)
In equations (2.56) and (2.59) uPf and uSf are forward propagated P and S waves at the exit side of the slab, and uPb and uSbare backscattered P and S waves at the entrance side of the slab. They are calculated by summing up the incidence waves uP0 and uS0, and
scattered waves UPPf , USPf , USSf , UPSf , UPPb , uSPb , USSb and UPSb using equations (2.57)-(2.58) and (2.60)-(2.61). These scattered fields can be calculated using equations (2.26)-(2.32) if thin slab approximation is adopted. In the latter case, the scattering patternsηcan be calculated using either zero or first order approximations (for details see Wu et al., 2007).
2.3.4
Numerical Implementation
One-return treatment of the propagator
Equations (2.39) and (2.40) (for the thin-slab approximation), or (2.56) and (2.59) (for the screen approximation) provide the solution to calculate the response of a thin slice of the model to the incident wave. To calculate the response of the entire model to the incidence wave, we use the multiple forward scattering and single back scattering approximations described in Section 2.2. A marching method is adopted for updating the wavefield as illustrated in Figure 2.2. A double-sweep method is used to calculate the transmitted and reflected (backscattered) wavefield. We first calculate the down-going wavefield from the source to the bottom of model (down sweep). For each depth step, the interaction of the forward (here it is downward) field with the thin-slab is updated and propagated to the next depth using equations (2.56)-(2.58), and is used as the incidence for the next depth. In a marching way, all multiple forward scatterings are included in the down-going wave. At the same time, at each depth, the back scattered wavefields are calculated using equations (2.59)-(2.61) and recorded in the memory or disk. After the down-going propagation sweeps the entire model, we use the same propagator for the second time to calculate the up-going wave from the bottom to the top of the model (up sweep). At each depth, the up-going waves are updated by the thin-slab transmission operator. In addition, the previously recorded back scattered waves are retrieved at each depth and added to the up-going field. After the up-sweep of the entire model from the bottom to the top, we obtain all the back scattered waves under the one-return approximation. By using the marching algorithm in the double sweeps, all multiple forward scattered and single back scattered waves are included in the results. If we calculate and keep the records of the back scattered fields in the up-sweep too, then we can model the third and higher-order multiples in the subsequent sweeps.
Dual domain implementation
In calculating scattering waves, i.e., equations (2.26)-(2.32) or equations (2.48)-(2.55), interactions between the incident wave and the perturbations of elastic parameters are cal-culated in the space domain. However, the forward propagation of the wavefield in the background velocity is calculated in the wavenumber domain. The inverse Fourier trans-form in equations (2.56) and (2.59) transtrans-form the wavefield from wavenumber domain to the space domain, while the Fourier transform in equations (2.26)-(2.32) or equations (2.48)-(2.55) transfer the wavefield from space domain to the wavenumber domain. In both domains, the calculations are localized, forming a very efficient algorithm.
2.3.5
Elastic, Acoustic and Scalar Cases
The equations presented in the previous sections are for elastic wave case, where both P-and S-waves are multi-component vector waves. Under the thin-slab approximation, the interactions between the medium and waves involve tensor (strain field) operations. Due to the symmetric properties of these tensors, there are only 6 independent components for each tensor. If complex screen approximation is adopted, the interactions involve only vector operations. The equations for acoustic and scalar waves are not included in the current chapter. However, the interested readers are referred to the papers by Wu (1996) and Wu et al. (2007).
2.4
Applications of One-Return Propagators in Modeling, Imaging
and Inversion
2.4.1
Applications to Modeling
First we show the applications of one-return approximation to the calculation of synthetic seismograms in elastic wave modeling. Figure 2.3a shows a 3D French model, where the parameters of the background medium areα0=3.6 km/s,β0=2.08 km/s andρ0=
2.2 g/cm3, and the grey colored structure has a perturbation of−10% for both P- and S-wave velocities. A Ricker wavelet with a dominant frequency of 10 Hz is used. For the 8th-order 3D elastic finite-difference method (Yoon, 1996), the spacing interval is 20 m. The actual grid size used is 250× 250 × 250 including 25 grids of absorbing boundary for each face of the model. Time interval used is 0.001 sec and 2500 time steps are calculated. The calculations took approximately 28 hours. For thin-slab method, the spacing interval used is 20 m in transversal plane, which is the same as that used for the finite-difference method. But a fine grid size of 5 m is used in propagation direction. We did the same calculation on the same machine using the thin-slab propagator. It took 2.7 hours. We see that the two methods agree with each other fairly well, but the thin-slab is about 10 times faster than the finite-difference method.
Shown in Figure 2.4 are primary waves calculated in the 2D acoustic SEG/EAGE salt model using the one-return method. These snapshots are composed of down- and up-going waves. We see that the primary reflections from major sediment/salt boundaries and from the interfaces in the sedimentary layers are properly modeled.
Figure 2.5 shows the comparisons of reflection coefficients calculated by the thin-slab method and the exact solutions (Reflection/Transmission theoretical calculation). The upper panel corresponds to P wave and the lower panel to S wave incidences respectively. The perturbation of the bottom layer is 20% with respect to the top layer for both P and S wave velocities. We can see the good agreement between these two methods.
Figure 2.6 gives one application example of the one-return modeling to the reservoir AVO (amplitude variation with offsets) calculations. It models the combined effects of random scattering and intrinsic attenuation in oil/gas reservoir AVO. In practice, the
geo-Fig. 2.3 Synthetic examples calculated in a 3D velocity model. Illustrated in (a) is the velocity model modified
from the French model (French 1974), and in (b) are comparisons between the synthetic seismograms calculated using the one-return method (blue lines) and the finite-difference method (red lines).
logic models may contain arbitrary spatial variations in P- and S-wave quality factors, as well as density and velocities. The Q values vary from 10 to 150 for different reservoir rocks. The correlation lengths of the random field for perturbing Q and elastic parameters are the same. The rms values are 4% for elastic parameters and 25% for Q. The source and receiver array is located in shale and 1200 m away from the interface. The dotted lines in Figure 2.6 correspond to the homogeneous cases with constant Q’s. The one-return modeling method can be applied to the complicated geological setting and the calculation is quite efficient. We see from the calculated AVO that intrinsic attenuation mainly affects the absolute reflected amplitudes and heterogeneities in Q’s and elastic parameters affect local amplitude fluctuation with offset. AVO responses of the target interfaces have been significantly deformed due to both the heterogeneities and intrinsic attenuation.
Fig. 2.4 Snapshots from the 2D acoustic SEG/EAGE salt model. The source is an 18 Hz Ricker wavelet. The
Fig. 2.5 Comparisons of reflection coefficients calculated by the thin-slab method and the exact solutions. The
upper panel corresponds to P wave and the lower panel to S wave incidences respectively. The perturbation of the bottom layer is 20% with respect to the top layer for both P and S wave velocities.
Fig. 2.6 An example showing the one-return modeling applied to the reservoir AVO (amplitude variation with
offsets) calculation. The upper left panel shows a reservoir model. The formation is anelastic and heterogeneous. The rest three panels show AVO’s for various kinds of interfaces: (a) shale/gas, (b) shale/oil, and (c) shale/brine. Three different background Q′
S(Q = infinity, 150, 50) are given to shale. The sand has background QP= QS=
10. The correlation lengths of the random field for perturbing Q and elastic parameters are the same. The rms values are 4% for elastic parameters and 25% for Q.
2.4.2
One-Return Propagators Used in Migration Imaging
Conventional one-way wave equation migration handles images of primary reflections, but neglects turning waves, multiples and other multiple reflections such as duplex waves. Duplex waves exist in the geologic structures with vertical features such as faults and flanks of salt bodies. In this case, primary reflections of steep reflectors may not be recorded in a limited acquisition aperture. Therefore, one-way wave equation migration cannot produce the image of such events. On the other hand, doubly reflected duplex waves may be recorded and should be taken into consideration in the migration. Based on the concept of multiple fore-scattering and single back-scattering (MFSB) a one-return wave equation migration that extrapolates waves with a single returning point can be ap-plied to the duplex waves. The principle of one-return propagator with double sweeps has been shown in Section 2.2. Followed by a properly designed imaging condition, one-return migration produces the depth image of primary reflections, the same as con-ventional one-way wave equation migration, as well as the image from the contribution of turning waves and duplex waves that one-way wave equation fails to handle (Jin et al., 2006; Xu and Jin, 2006).
A conventional one-way propagator only calculates waves propagated downward (down-going or back propagated up-(down-going waves), while the one-return propagator can handle waves reflected from a interface then propagate upward. Thus, it can generate different source and receiver wave modes including PSD, PSU, PGDand PGU (see Figure 2.7), where subscripts S and G denote source and receiver (geophone) side waves, and superscripts D and U denote down- and up-going waves. Accordingly, the corresponding partial images generated by these wave modes are (Jin et al., 2006)
IDD=
∑
ωP D SPD ∗ G , (2.62)Fig. 2.7 Illustration of imaging conditions for different combinations of source and receiver waves. IDDis the
conventional imaging condition for down-/down-going waves of source and receiver sides. IDUis an imaging
IUU =
∑
ωP U SP U∗ G (2.63) and IDU=∑
ωP D SPU ∗ G +∑
ω P U SPD ∗ G , (2.64)where “∗” denotes applying complex conjugate. In these equations, IDDis the image for
down-/down-going waves for both source and receivers and it is the same as the traditional one-way wave equation migration (refer to Figure 2.7). IDU is the image for
down-/up-and up-/down-going waves for both source down-/up-and receivers, down-/up-and it hdown-/up-andles the duplex waves. IUU is the image for up-/up-going waves. This image condition actually handles multiples
but will not be covered in this chapter. The final image can be obtained by summing up these partial image volumes. Here, we use synthetic examples to demonstrate the unique capability of this approach to image the steep and overhang structures.
The one-return propagator can be used to migrate duplex waves and prism waves. Shown in Figure 2.8a is a simple two-layer model with a vertical fault. The wave path la-beled with “A” is the direct arrival; lala-beled with “B” is the primary reflection and lala-beled with “C” is a doubly bounced reflection against the vertical fault wall. The double bounc-ing reflection is also called duplex wave which is a special kind of prism wave. Illustrated in Figure 2.8b is the wavefield snapshot at 2.4 second simulated using finite-difference forward modeling. Different phases can be clearly seen from this snapshot. Shown in
Fig. 2.8 Image vertical fault wall using one-return propagator. (a) A two-layer velocity model with a vertical
fault wall. The direct arrival, the primary reflection and the doubly bounced reflection against the vertical fault wall are all illustrated in the figure. (b) The wavefield snapshot at 2.4 second. (c) Depth image using the conventional one-way wave-equation-based propagator. (d) Depth image using the one-return propagator. Note that the vertical fault wall has been properly reconstructed.
Figure 2.8c is the depth image which is calculated using the one-way wave equation mi-gration and from a single shot record. The vertical fault structure is totally missing from this image. The depth image in Figure 2.8d is calculated using the one-return wave equa-tion migraequa-tion and from the same shot record. The vertical fault wall is fully reconstructed by duplex waves (Jin et al., 2006).
The one-return propagator can be also used to migrate turning waves that propagate be-yond 90 degrees, thus to image overhang structure. Shown in Figure 2.9a is a typical salt model with a steeply dipping overhang. Note that primary reflections propagate through salt body while strong overturned waves exist in the sediments. Figure 2.9b illustrates the depth image generated using only one-way propagator. Due to the lack of contribu-tions from turning waves, the image amplitude of salt overhang is very weak. Shown in Figure 2.9c is the partial image generated by the contributions from turning waves using one-return wave equation migration. Illustrated in Figure 2.9d is the final depth image by summing up images from both one-way propagator and one-return propagator. In this way, the image amplitude of overhang salt flank is significantly enhanced (Xu and Jin, 2006).
Fig. 2.9 Imaging the overhang salt flank using both one-way and one-return propagators. (a) A typical salt
model with a steeply dipping salt overhang, noticing that primary reflections propagate through salt body while strong overturned waves exist in the sediments; (b) Depth image using one-way wave equation migration; (c) Partial image contributed by turning waves using one-return wave equation migration; (d) Final depth image by summing up both (b) and (c).
The next example is for the BP 2004 Molar benchmark dataset. Illustrated in Figure 2.10a is the central part of the velocity model which consists of a deeply rooted salt body with sedimentary inclusions which present a lot of challenges in depth imaging. Figure 2.10b shows depth image using one-way wave-equation-based propagator. The boundary of the deeply rooted salt body is not imaged. Shown in Figure 2.10c is the partial image calculated using the one-return wave equation migration. By including the contributions from both the prism waves and turning waves, the image of the steeply dipping salt boundary is clearly seen. Finally, shown in Figure 2.10d is the full image of one-return wave equation migration. The boundaries from both the top of the salt crown and the deeply rooted salt body are well reconstructed (Jin et al., 2006; Xu and Jin, 2006, 2007).
Fig. 2.10 One-return wave equation migration for the BP 2004 Molar benchmark dataset. (a) The central part of
the velocity model; (b) Depth image using only the one-way wave-equation-based propagator; (c) Partial image using the one-return wave equation migration; (d) the full image by summing up both (b) and (c). (Synthetic data courtesy of BP).
2.4.3
Calculate Finite-Frequency Sensitivity Kernels Used in Velocity
Inversion
In seismic migration, a correct velocity model plays an important role in obtaining high quality image. The process to update the velocity model based on the migration image is
the migration velocity analysis (MVA), which is a special type of velocity tomography. Comparing to the conventional tomography, where the information regarding the velocity model error is extracted from the data domain (seismograms), the MVA extracts the infor-mation from inconsistency of the depth image. The most commonly used inconsistency is the residual moveout (RMO) in different types of common image gathers. To update the velocity model, the most important part in MVA is converting the observed residual move-out into velocity corrections and back-projecting them into the model space for velocity updating. In the past, this was dominated by the ray-tracing-based tomography method which assumes an infinitely high frequency. Recently, the sensitivity of finite-frequency signals to velocity model has been investigated by researchers working in different fields (Woodward, 1992; Vasco et al., 1995; Dahlen et al., 2000; Zhao, et al., 2000; Skarsoulis and Cornuelle, 2004; Spetzler and Snieder, 2004; Sava and Biondi, 2004; Jocker, et al., 2006; and Buursink and Routh, 2007; Fliedner et al., 2007). Finite-frequency sensitiv-ity kernels have been calculated and used for solving many tomography problems with great success. The major obstacle that prevents this method from being used in migration velocity analysis is that these finite-frequency sensitivity kernels are mostly derived for transmitted waves (e.g., travel time delays or amplitude fluctuations in seismograms). On the contrary, the seismic migration extracts the information regarding the velocity error from the reflectivity image, which is related to the calculation of reflect waves. Based on the Born and Rytov approximations, Xie and Yang (2007, 2008) formulated the sensitiv-ity kernel for the shot-record prestack depth migration. This sensitivsensitiv-ity kernel relates the observed RMO in depth image to the velocity correction in the model. This is a wave-equation-based method which avoids many disadvantages of the ray-based tomography. The one-return method is ideal in calculating sensitivity kernels involving reflection. The single frequency travel time sensitivity kernel is composed of two parts, the down-going leg and the up-going leg.
KDF(r, rS, rI,ω) = imag 2k20GD(r; rS,ω) G (r; rI,ω) GD(rI; rS,ω) , (2.65) KUF(r, rS, rI,ω) = imag 2k20G ∗ U(r; rS,ω) G (r; rI,ω) GU∗(rI; rS,ω) , (2.66)
whereωis frequency, k0=ω/v(r) is the background wavenumber, v(r) is the background
velocity, G is the Green’s function, subscripts U and D are for up- and down-going waves,
r is the space location, rSand rIare the source and image locations, and imag (·) denotes
taking imaginary part. In equations (2.65) and (2.66), the sensitivity kernels are composed by the correlation of two Green’s functions, one radiated from the source and the other from the image point, normalized by the source wavefield at the image location.
The major difference between the up- and down-going legs is that in prestack depth migration, the up-going wave is obtained from the time-reversed reflect wave. The time reversal in the time domain is equivalent to the complex conjugate in the frequency do-main. Thus complex conjugate is applied in the up-going leg. The broadband sensitivity kernel can be obtained by integrating the single frequency sensitivity kernel over
frequen-cies KDB(r, rS, rI) = Z A(ω) ω KDF(r, rS, rI,ω)dω, (2.67) KUB(r, rS, rI) = Z A(ω) ω KUF(r, rS, rI,ω)dω (2.68)
where A(ω) is a factor related to the source spectrum (Xie and Yang, 2008). Finally, the sensitivity kernel for the residual moveout R(rI, rS) can be obtained as follows:
R(rI, rS) = −
v0(rI)
2 cos[θ(rI, rS)]δ
t(rI, rS) (2.69)
Fig. 2.11 The Green’s functions used to construct the sensitivity kernel for migration velocity analysis. (a)
with the factor v0(rI)/2 cos[θ(rI, rS)] being a local time-depth convertor andδt=δtD+ δtU, where δtD(rI, rS) = Z V m r′ KDB r′, rS, rIdv′, (2.70) δtU(rI, rS) = Z V m r′ KB U r′, rS, rIdv′. (2.71)
Equations (2.69)-(2.71) form the velocity inversion system. The cartoon in Figure 2.11 illustrates the three Green’s functions used to calculate the sensitivity kernel for velocity updating, where Figure 2.11a and Figure 2.11b is for down- and up-going waves, and Figure 2.11c is for waves radiated from the image point. The up-going Green’s function in Figure 2.11b is calculated using the one-return method.
Shown in Figure 2.12a is a 10 Hz broadband sensitivity kernel for migration geometry calculated using the one-return method. The source side sensitivity kernel is similar to that for transmitted waves. However, the receiver side kernel is quite different. Near the source location and immediately above the image point, there are two sensitive regions where the velocity model affects the RMO most. As comparisons, shown in Figure 2.12b is the actually measured sensitivity map from the migration process (Xie and Yang, 2008). We see these kernels show excellent consistency.
Figure 2.13 compares sensitivity kernels calculated for selected image points in a five-layer/four-reflector model. Shown in the left column is sensitivity kernels calculated using
Fig. 2.12 The sensitivity kernels for a shot gather image, with (a) theoretical sensitivity kernel, and (b) the
one-return method. As a comparison, the right column shows the actually measured sen-sitivity maps in the same model. The results show that the sensen-sitivity kernels calculated using the one-return method are consistent with the measured sensitivity maps.
Fig. 2.13 Comparison between the theoretically calculated kernels (left column) and actually measured
2.5
Other Development of One-Return Modeling
One-way and one-return propagators have great advantages over the full-wave propaga-tor, especially for the applications to imaging and inversion. The method is efficient and produces less artifacts compared with full-wave methods when applied to backpropaga-tion and imaging. It has the flexibility to select different conversion modes in elastic wave propagation and imaging. It can model multiple scattering in the order of scattering mul-tiplicity by multi-sweeps. It has great potential in applying to different inversion schemes. However, it has also some drawbacks which prevent its use in some applications. One se-rious drawback is the angle limitation. Even the wide-angle one-way propagators become much less accurate when the propagation angles are close to 90◦(relative to the preferred direction, z-axis). Standard one-way propagators cannot go beyond 90◦. Even though the one-return propagator can model the turning wave by double-sweeps, the waves prop-agating nearly along the horizontal direction will have less accuracy in both phase and amplitude. The research in improving accuracy of wide-angle or super-wide angle (be-yond 90◦) is one of the future directions of one-way, and one-return modeling. Another drawback of the one-return modeling is the convergence of the marching algorithm in strong contrast media. The one-return approximation, or the De Wolf series, is based on the perturbation series formulation. The elastic one-return methods using elastic thin-slab propagator or elastic complex-screen propagator have to limit the applications to moderate contrast media (parameter perturbations smaller than 40%) (Wu et al., 2007). Reflectivity method has been incorporated into the one-return method (Wu and Wu, 2003), but is lim-ited only to flat layers. In the following, we summarize the progress in overcoming the two basic drawbacks described above. There are other developments in way, and one-return modeling, such as one-way propagator in anisotropic media (Angus et al, 2004), and improvement of the amplitude accuracy by multi-one-way modeling (Kiyashchenko et al., 2005), which will not be discussed in this chapter.
2.5.1
Super-Wide Angle One-Way Propagator
In order to improve the accuracy of wide-angle waves and overcome the fundamental an-gle limitation of one-way propagators, some work has been done using tilted coordinates or curved coordinate system (Shan et al., 2009; Shragge and Shan, 2010). The other approach is to use two orthogonally propagated one-way propagators to reconstruct the accurate wavefront up to nearly the full angle range (Wu and Jia, 2006; Xu and Jin, 2007; Jia and Wu, 2009a, 2009b). As shown in Figure 2.14, large-angle waves with respect to z-axis become small-angle waves to x-axis. Therefore, wavefront reconstruction method using weighted average of the two orthogonal propagated waves can keep good accuracy to super-wide angle ranges. Figure 2.15 shows the comparison of impulse responses be-tween the regular one-way propagator (blue curve) and the superwide-angle propagator (red curve). The finite difference wavefront (green curve) is also shown as reference. The velocity of the model varies linearly in both lateral and vertical directions. We see that the wavefront from the regular one-way method has been distorted significantly at
Fig. 2.14 Wavefield and wavefront reconstruction by taking a weighted average of two orthogonally propagated
one-way wavefields.
L
Fig. 2.15 Comparison of wavefront accuracy between regular one-way propagator (GSP) (blue curve) and
superwide-angle GSP propagator (red curve). The finite difference wavefront (green curve) is shown as refer-ence. The model has velocity gradients along±x and z direction, v(x,z) = v0+ |x − 6.12|v0/6.12 + z0v0/6.12.
large propagation angles (close to 90◦). On the other hand, the new propagator models the wavefront accurately up to 135◦. The performance is degenerated only for nearly backpropagation angles (close to 180◦). Figure 2.16 compares the prestack images by the regular one-way propagator (left panel) and by the superwide-angle one-way propagator (right panel) for the 2D BP benchmark model. Due to the angle limitation, the regular downward one-way propagator can hardly simulate turning waves and therefore cannot image the overhanging flank. By contrast, these limitations are eliminated in the imaging
Fig. 2.16 Comparison of images obtained by the regular GSP migration (a) and by the superwide-angle GSP
migration (b).
result of the superwide-angle method. The image of overhanging flank is sharp and clear and its location is accurate.
2.5.2
One-Way Boundary Element Method
As we pointed out, the one-return modeling is based on the perturbation method and has difficulty in applying to strong contrast media, such as the media with complicated salt and basalt structures. For salt models, the velocity perturbation could be more than 200%. One-way boundary-element method is proposed to partially solve the problem (He and Wu, 2009). The concepts of one-way and one-return boundary element method can be illustrated by the simple models shown in Figure 2.17a: a layered model, and 2.17b: an inclusion model. The key operation for one-way boundary-element method is to decouple multiple interactions between the upper boundary S1and the lower boundary
S2and therefore eliminate the internal multiples between these two boundaries. It is quite
straightforward to apply this concept to layered media. However, to extend this method to models with an inclusion, extra care must be taken. Here, we take a simple inclusion model (Figure 2.17b) as an example. The inclusion could be a salt dome. We separate the boundary of the salt dome into top part S1and bottom part S2based on the shape of
