Stable and efficient Q-compensated least-squares migration with compressive sensing, sparsity-promoting, and preconditioning

Journal of Applied Geophysics 145 (2017) 84–99

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Stable and efficient Q-compensated least-squares migration with compressive sensing, sparsity-promoting, and preconditioning Xintao Chai a, * , Shangxu Wang b , Genyang Tang b , Xiangcui Meng b a

China University of Geosciences (Wuhan), Institute of Geophysics and Geomatics, Center for Wave Propagation and Imaging (CWp), Hubei Subsurface Multi-scale Imaging Key Laboratory, Wuhan, Hubei, China b China University of Petroleum (Beijing), State Key Laboratory of Petroleum Resources and Prospecting, China National Petroleum Corporation (CNPC) Key Laboratory of Geophysical Exploration, Changping, Beijing, China

A R T I C L E

I N F O

Article history: Received 16 March 2017 Received in revised form 23 May 2017 Accepted 25 July 2017 Available online 12 August 2017 Keywords: Least-squares migration Q-compensation Compressive sensing Sparsity-promoting Preconditioning Stability

A B S T R A C T The anelastic effects of subsurface media decrease the amplitude and distort the phase of propagating wave. These effects, also referred to as the earth’s Q filtering effects, diminish seismic resolution. Ignoring anelastic effects during seismic imaging process generates an image with reduced amplitude and incorrect position of reflectors, especially for highly absorptive media. The numerical instability and the expensive computational cost are major concerns when compensating for anelastic effects during migration. We propose a stable and efficient Q-compensated imaging methodology with compressive sensing, sparsity-promoting, and preconditioning. The stability is achieved by using the Born operator for forward modeling and the adjoint operator for back propagating the residual wavefields. Constructing the attenuation-compensated operators by reversing the sign of attenuation operator is avoided. The method proposed is always stable. To reduce the computational cost that is proportional to the number of wave-equation to be solved (thereby the number of frequencies, source experiments, and iterations), we first subsample over both frequencies and source experiments. We mitigate the artifacts caused by the dimensionality reduction via promoting sparsity of the imaging solutions. We further employ depth- and Q-preconditioning operators to accelerate the convergence rate of iterative migration. We adopt a relatively simple linearized Bregman method to solve the sparsity-promoting imaging problem. Singular value decomposition analysis of the forward operator reveals that attenuation increases the condition number of migration operator, making the imaging problem more ill-conditioned. The visco-acoustic imaging problem converges slower than the acoustic case. The stronger the attenuation, the slower the convergence rate. The preconditioning strategy evidently decreases the condition number of migration operator, which makes the imaging problem less ill-conditioned and significantly expedites the convergence rate. Numerical experiments verify the stability, benefits and robustness of the method proposed and support our analysis and findings. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Due to the fact that real earth media are anelastic, the propagation of seismic waves in field media is in many respects significantly different from propagation in perfectly elastic media (Mittet et al., 1995). The anelastic properties of the subsurface media cause dissipation of wave energy, thus decreasing the amplitude and distorting the frequency content of the wavelet propagating inside the earth

* Corresponding author at: China University of Geosciences (Wuhan), No. 388 Lumo Road, Wuhan 430074, Hubei, China. E-mail addresses: [email protected] (X. Chai), [email protected] (S. Wang), [email protected] (G. Tang), [email protected] (X. Meng).

http://dx.doi.org/10.1016/j.jappgeo.2017.07.015 0926-9851/© 2017 Elsevier B.V. All rights reserved.

(Quan and Harris, 1997). Fluids and especially gas trapped in overburden structures cause strong attenuation and dispersion of the seismic wave, as observed over many oil-gas fields (Zhu et al., 2014). A quality factor Q was introduced to quantify the anelastic attenuation and dispersion effects (Futterman, 1962; Kjartansson, 1979; Knopoff, 1964; Kolsky, 1956; Morozov and Ahmadi, 2015), which are now generally referred to as Q filtering effects (Wang, 2006). The term Q is used to represent the energy lost for a wavefield propagating along one wave length and is often defined as Q = 2pE/DE, where E is the energy of the wavefield and DE is the energy lost in a wavelength of propagation. Ignoring the Q filtering effects during seismic migration/imaging processes, for example, acoustic reverse time migration (RTM) (e.g., Baysal et al., 1983; Chattopadhyay and McMechan, 2008; Levin, 1984)

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

and acoustic least-squares migration (LSM) (e.g., Herrmann and Li, 2012; Kazemi and Sacchi, 2015; Nemeth et al., 1999; Wu et al., 2015), would produce an image with reduced amplitude and incorrect position of reflectors, especially those below highly attenuating geologic environments (Zhu et al., 2014). To extract more detailed information from Q-filtered seismic data, we must account for the Q filtering effects. By taking the anelastic effects into account during migration, the resultant image is expected to have high resolution with true amplitude and correct timing, which are important for both regional geophysics investigation and reservoir exploration. Mittet et al. (1995) showed that even when imaging with a Q model deviating by 10% from the true one, images would still be better focused and of a higher quality than when no Q-compensation was performed (Wang, 2008). There are methods mitigating the Q filtering effects based on the one-dimensional (1D) wave equation (e.g., Chai et al., 2017, 2016a, 2014; Wang, 2006). However, because the anelastic effects on recorded seismic data occur during the wave propagation and are related to wavepaths for 2D/3D surveys, it is more accurate and physically more consistent to address these effects in a way of 2D/3D wave-equation based pre-stack depth migration. Conventional method to compensate for the Q filtering effects is to reverse the sign of attenuation operator. The method is unstable, since it supports waves with exponential growth. The amplitudes of high-frequency components grow rapidly in the wavefield extrapolation, numerical round-off errors tend to be amplified drastically. This causes a huge amount of undesirable noise in the image, even if the input data are noise-free (Wang, 2006, 2008). Fig. 1 illustrates the intrinsic numerical instability associated with inverse Qcompensation. The scale of the amplitude compensation by reversing the sign of the attenuation operator grows rapidly (see Fig. 1b), which introduces a large number of artifacts into the inverse Q filtered seismogram (see Fig. 1c). The numerical instability is a major concern when compensating for the Q filtering effects (Li et al., 2016; Zhu et al., 2014). Mittet et al. (1995) design one-way compensation extrapolators for up- and down-going waves to address the complete amplitude loss from source to receiver. To stabilize the results, they use a high-cut frequency filter on the initial state for up- and downgoing wavefields. Zhang and Wapenaar (2002) suggest limiting the number of extrapolation steps and the maximum angle of migration dip to derive a conditionally stable extrapolation operator. Mittet (2007) compensates for the Q filtering effects in deriving the space-frequency domain depth extrapolation operators by starting from the ones without such compensation, and needs to consider the stability problem. To stabilize migration with inverse-Q filtering (Wang, 2008), an inverse method implemented by adding a small stabilization factor in the inverse of the operator, is required and implicitly introduces a gain limit. Zhang et al. (2013) develop an effective Q method of compensating for absorption and dispersion in migration. Stabilization is achieved by introducing a smooth, maximum-limited gain function that matches the exact amplitude compensation factor when it is less than the user-specified gain limit. Zhu et al. (2014) propose an approach of compensating for attenuation effects in reverse-time migration (Q-RTM), in which the attenuation- and dispersion-compensated operators were constructed by reversing the sign of attenuation operator and leaving the sign of dispersion operator unchanged. They design a low-pass filter to stabilize the compensating procedure. The visco-acoustic migration methods described in Sun et al. (2015), Li et al. (2016) also apply low-pass filters to prevent high wavenumber noises from growing exponentially. Based on a standard linear solid model (Carcione et al., 1988; Zhu et al., 2013), Dutta and Schuster (2014) use the time-domain visco-acoustic wave equation for wavefield extrapolation and use the linearized least-squares inversion method to compensate for the attenuation loss. The method is denoted as Q-LSRTM. They need not

85

to regularize the adjoint wave equations during the receiver-side wavefield extrapolation to compensate for the attenuation loss. Their method is stable. The shortcoming of Q-LSRTM is that this method is expensive (the computational cost per iteration is more than six times that of standard RTM, and the total computational cost is proportional to the number of least-squares iterations). With the sparse-recovery and randomized sampling theory from compressive sensing (CS) (Candès et al., 2006b; Donoho, 2006), and the depth- and Q-preconditioning strategies, we propose a stable and efficient Q-compensated iterative imaging approach in the frequency-domain. The frequency domain is the natural domain for dealing with dispersion and attenuation effects (Song and Williamson, 1995). First, attenuation parameters and dispersive velocity relationships can be incorporated into forward modeling through the use of complex-valued frequency-dependent velocity (Marfurt, 1984; Prieux et al., 2013; Song et al., 1995; Tarantola, 1988). Both phase velocity and attenuation coefficients are allowed to vary with frequency and spatial location. Second, it is flexible for implementing any kind of viscous attenuation mechanism. Besides, for certain geometries, only a few frequency components are required to perform wave-equation based inversion and imaging (Chen et al., 2013; Pratt et al., 1998). During the imaging process, we use the Born operator for forward modeling and the adjoint operator for back propagating the residual wavefields, thus avoiding the construction of attenuation compensation operator. To reduce the computational cost, which is proportional to the number of partial-differential equation (PDE) to be solved (thereby the number of frequencies, source experiments and iterations), we subsample over both frequencies and source experiments. We suppress the artifacts caused by the dimensionality reduction by promoting sparsity of the imaging solutions. Thus we are using a linearized Q inversion method (Q-LSRTM) like Dutta and Schuster (2014) combined with sparsity promoting imaging techniques to compensate for the attenuation loss during migration. Since seismic attenuation may significantly slow down the convergence rate of the least-squares iterative inversion process without proper preconditioning (Sun et al., 2016), we employ a depth-preconditioning operator and a Q-preconditioning operator to accelerate the convergence rate of iterative migration. Following Herrmann et al. (2015), Yang et al. (2016), we adopt a relatively simple linearized Bregman (LB) method (Lorenz et al., 2014b) to solve the sparsity-promoting visco-acoustic imaging problem. Singular value decomposition (SVD) analysis of the forward operator reveals that attenuation increases the condition number of migration operator, making the imaging problem more ill-conditioned. The stronger the attenuation, the slower the convergence rate. The preconditioning operator evidently reduces the condition number of migration operator making the imaging problem less ill-conditioned, and accelerates the convergence rate of LSM. The rest of this paper is organized as follows. After notation declarations, we briefly introduce the basic theory of visco-acoustic seismic modeling and imaging. We then describe the proposed approach in detail. Numerical experiments serve to support our analysis and findings. After discussing the limitations and extensions of the method, we summarize our conclusions. 2. Theory 2.1. Notations Throughout this paper, unless specified otherwise, lowercase boldface symbols (e.g., x) refer to column vectors and uppercase boldface symbols (e.g., X) refer to matrices or operators. The jth element of a vector x is denoted as xj . I is the identity matrix. 0 is the zero vector. log( • ) is the natural logarithm function. log b (x) denotes the logarithm of x to base b. i is the imaginary unit. diag( • ) when

86

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

applied to a vector constructs a diagonal matrix. vec( • ) means to vectorize a matrix. The symbol ‘’ means ‘defined as’. The superscripts T and H denote transposition and (Hermitian) conjugate transpose, respectively. The superscript * denotes the adjoint.  • 2 and  • 1 represent the vector 2 - and 1 -norms, respectively. | • | gets the absolute value. max ( • ) obtains the maximum value. The Kronecker product is denoted by ⊗. 2.2. Wave-equation based modeling The 2-D constant-density visco-acoustic wave equation in the frequency domain reads (Zhang and Wapenaar, 2002)

∂ 2 u(x, z, y) ∂ 2 u(x, z, y) + + my2 a2 (y, Q)u(x, z, y) = s(x, z, y), (1) ∂ x2 ∂ z2 with the complex wavefield u(x, z, y) depending on the space positions (x, z), the angular frequency y, the source term s(x, z, y), and the model parameter m = v12 with v the real-valued velocity. According

to the Kolsky-Futterman Q model (Futterman, 1962; Kolsky, 1956), the attenuation coefficient a(y, Q) is given by      1 y i a(y, Q) = 1 − log , 1− pQ yr 2Q

where yr is the reference frequency. It is well known that constant Q may be a fair approximation for acoustic energy in the seismic frequency band (Mittet, 2007). According to the constant-Q model theory (Kjartansson, 1979; Liu et al., 1976), Q is assumed approximately independent of frequency within the seismic bandwidth. When Q → +∞, a(y, Q) → 1, Eq. (1) is reduced to the acoustic wave equation. Eq. (1) can be solved using the finite-difference method (Chen et al., 2013). Discretizing Eq. (1), the wave equation can be expressed compactly in a matrix form of Hj (m)Uj (m) = Sj ,

(3)

where Uj (m) is the complex wavefield for a single frequency. The subscript j = 1, 2, · · · , Nf indexes the frequencies with Nf being the total number of discretized frequencies. Sj contains the complex source term. Hj (m) is the discretized monochromatic Helmholtz operator, Hj (m) = y2j diag(m)a2 (yj , Q) + L,

Fig. 1. The intrinsic numerical instability associated with inverse Q-compensation. A 50-Hz zero-phase Ricker wavelet is used as the source wavelet. The time sampling interval is 2 ms, and the seismogram consists of a series of wavelets at the time t = 0.1, 0.4, · · · , 1.9 s. (a) Synthetic Q-filtered seismogram with Q = 25 constant with depth, clearly showing the Q filtering effects. Panel (b) plots the value of exp [pft/Q eff ], which describes the amplitude compensation. An effective Q value Q eff represents the Q-filtering effects accumulated from time zero to the current travel time. (c) Results of full and exact inverse Q filtering for the seismogram in panel (a).

(2)

(4)

Fig. 2. Transform coefficients of the medium perturbation dm (i.e., the target image), sorted from large to small, normalized and plotted in percentage of the number of the transform coefficients. Thin line — physical domain (i.e., no transform is applied). We plotted the coefficients for the model in the physical domain by directly sorting the (normalized) amplitudes of the medium perturbation dm from large to small. Bolded line — curvelet domain. (a) For dm of the following graben model example. (b) For dm of the following Sigsbee2B model example.

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

which is square (with dimensions of Nx Nz × Nx Nz ) and extremely sparse. Nx and Nz are the numbers of grid points in the model in the x- and z-directions. L is the discretized Laplacian (van Leeuwen and Herrmann, 2013). The term a(y, Q) is a vector composed of the attenuation coefficients. For simplicity, such solutions to the problem (3) can be represented as Uj (m) = H−1 (m)Sj . j

(5)

Eq. (5) is only representational, as no matrix inverse is actually computed (Hicks and Pratt, 2001), and is often solved using direct matrix factorization methods, such as a lower-upper (LU) factorization. Once LU decomposition is performed to solve Eq. (5), the matrix factors can be re-used to solve rapidly the forward problem for any new source vector. Given the source positions ps = (xs , zs ), the measured data at the jth frequency Dj , at the receiver positions pr = (xr , zr ), is calculated by (m)P∗s Wj , Dj = Pr H−1 j

(6)

87

where Dj ∈ CNr ×Ns , with Nr common-receiver gathers in its rows and Ns common-shot gathers in its columns. Pr is a detection operator that collects data at the receiver locations, and Ps injects the source wavefield. Throughout this paper, the source wavelet is assumed to be the same for all shots. Consequently, Wj ∈ CNs ×Ns is a diagonal matrix, with the diagonal entries being the Fourier spectrum of the source wavelet at the jth frequency wj , i.e., Wj = Iwj . We define the vectorized data at the jth frequency dj as   dj = vec Dj  F (m, Q, yj , Wj ).

(7)

2.3. Wave-equation based imaging 2.3.1. Born approximation Seismic migration aims to reconstruct the earth’s reflectivity image from the recorded waveform data under the single-scattering approximation (i.e., the first-order Born approximation) (Etgen et al., 2009; Nemeth et al., 1999; Tarantola, 1984). For the single-scattering approximation, let m = m0 + dm, where m0 = 12 represents v0

a smooth background model that describes the kinematics of the seismic waves and should not produce significant scattering in the seismic frequency band, and dm = v12 − 12 is the perturbation of the v0

Fig. 3. (a) Medium perturbation dm of a graben model. The overlaid dotted-white line shows the reference trace at the horizontal 1000 m. (b) Variation of the condition number j(A) of the matrix A in Eq. (26) with Q. Thin-circle line — without depthpreconditioning (i.e., Pdepth = I). Bolded-asterisk line — with depth-preconditioning (i.e., Pdepth defined by Eq. (20)). The condition numbers are plotted on a log10 scale.

Fig. 4. (a) Shortest-offset section of the visco-acoustic data. The overlaid white and black lines in panel (a) denote the sample seismograms of the acoustic- and viscoacoustic data, respectively. The thin and bolded lines in panel (b) represent the amplitude spectra of the acoustic- and visco-acoustic data, respectively.

88

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

Fig. 6. After the same number of iterations, (a) acoustic LSM image for the acoustic data, (b) visco-acoustic LSM image for the visco-acoustic data. By visco-acoustic LSM, we mean the proposed sparsity promoting Q-LSRTM method rather than the standard Q-LSRTM method.

subsurface medium that is responsible for the scattering of incoming waves. Then we have F (m, Q, yj , Wj ) ≈ F (m0 , Q, yj , Wj ) + ∇ F (m0 , Q, yj , Wj )dm.

(8)

The symbol ‘≈’ indicates that by the linearized modeling, the higherorder scatterings are neglected. ∇ F (m0 , Q, yj , Wj ) is the first-order derivative of F (m, Q, yj , Wj ) with respect to the model parameter m evaluated at m0 , which is known as the Born modeling operator representing the Jacobian matrix of first derivatives. ∇ F (m0 , Q, yj , Wj ) can be obtained by taking the partial derivative of both sides of Eq. (3) with respect to model parameter m

Fig. 5. For Q-filtered data, (a) migrated image produced by Eq. (15) (i.e., RTM image), (b) acoustic LSM image (i.e., without Q-compensation), (c) visco-acoustic LSM image generated by the method proposed. Following Fig. 3a, the overlaid solid-black line displays the correspondingly migrated trace at the horizontal 1000 m. Panel (d) plots the amplitude spectra for the profiles with and without Q-compensation to see the recovery in the high-frequencies/wavenumbers.

∂ Hj (m) ∂ Uj (m) U (m) + Hj (m) = 0. ∂m j ∂m

(9)

Then, we have 



∂ Uj (m) ∂ Hj (m) = H−1 U (m) . (m) − j ∂m ∂m j

(10)

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

89

Combining Eqs. (4) – (6), and (10), yields  = −Pr H−1 (m0 )y2j diag a2 (yj , Q) diag(dm)H−1 (m0 )P∗s Wj , dDlinear j j j (11) where dDlinear is linearized modeling data. The vectorized linear data j dDlinear is defined as j

 ddlinear = vec dDlinear j j  ∇ F (m0 , Q, yj , Wj )dm.

(12)

Expression (12) relates the medium perturbation dm to the data in a linear manner. The nonlinear data ddnonlinear is defined as j ddnonlinear = F (m, Q, yj , Wj ) − F (m0 , Q, yj , Wj ), j

(13)

which contains the higher-order scatterings. Therefore, ∇ F (m0 , Q, yj , Wj )dm = F (m, Q, yj , Wj ) − F (m0 , Q, yj , Wj ).

(14)

The above expression for multiple frequencies can be obtained by stacking over the rows (i.e., vertical concatenation along the columns). For numerical implementation, a perfectly matched layer (PML) boundary is included at the four edges of the model to reduce edge reflections (Berenger, 1994). We use an optimal nine-point finite-difference scheme for the Helmholtz equation with PML (Chen et al., 2013). 2.3.2. Basic equations of imaging A standard seismic migration operator can be regarded as the adjoint of a modeling operator as used in the iterations of fullwaveform inversion (Tarantola, 1984; Virieux and Operto, 2009). Applying the adjoint operator to the data is equivalent to backpropagating the recorded data or data residuals. The migrated image by applying the adjoint operator is given by N

dm =

f

∇ F ∗ (m0 , Q, yj , Wj )ddj .

(15)

j=1

Although an adjoint operator is a useful approximation to the inverse of the modeling operator, it is not the exact inverse. The accuracy of the migrated image can be significantly improved by LSM, which is an iterative migration method used to find the correct amplitudes by minimizing a least-squares misfit (Nemeth et al., 1999) xLS = arg min x

1

Ax − b 2 , 2 2

(16)

where x = dm, and ⎤ ⎡ ∇ F (m0 , Q, y1 , W1 ) dd1 ⎢ ∇ F (m0 , Q, y2 , W2 ) ⎥ ⎢ dd2 ⎥ ⎢ ⎢ ⎥,b = ⎢ . A=⎢ .. ⎥ ⎢ ⎢ . ⎦ ⎣ ⎣ . . ∇ F (m0 , Q, yNf , WNf ) ddNf ⎡

Fig. 7. (a) Migrated image without compressive sensing and sparsity-promoting (by setting S = I, and threshold l = 0). (b) Migrated image with compressive sensing and sparsity-promoting (via setting S = C, threshold l > 0, C is the curvelet operator). (c) Differences between (a) and (b).

⎤ ⎥ ⎥ ⎥. ⎥ ⎦

(17)

Because of the size of the imaging problem, matrix A cannot be formed explicitly. Wide adoption of these iterative methods is largely hindered by the expensive computational burden (Dutta and Schuster, 2014; Zhang and Ulrych, 2010), which is in proportion of the number of PDE to be solved and the number of least-squares iterations carried out. Each iteration roughly requires 4Ns Nf PDE to be

90

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

Fig. 8. (a) Migrated image without any preconditioning (by setting Pdepth = I, and PQ = I). (b) Migrated image with depth-preconditioning only (via setting Pdepth according to Eq. (20), and PQ = I). (c) Migrated image with both depth-preconditioning and Q-preconditioning (via setting Pdepth and PQ according to Eqs. (20) and (22), respectively). (d) Differences between (b) and (c). Panel (d) illustrates the benefits of employing the newly designed Q-preconditioning operator.

solved. The number of iterations done for LSM is the most important factor contributing to the computational expense. 2.4. Q-compensated compressive imaging To reduce the number of PDE to be solved, which is proportional to the number of frequencies Nf and source experiments Ns , inspired by Herrmann and Li (2012), we subsample over both frequencies and source experiments. We only use a randomized subset of frequencies, denoted by Y, with Nfsub being the number of the subsampled frequencies, Nfsub  Nf . We replace the forward modeling in Eq. (12) by   ddj = vec dDj E = ∇ F (m0 , Q, yj , Wj E)dm = ∇ F (m0 , Q, yj , Wj )dm.

(18)

The underlined quantities refer to subsampled source and receiver wavefields. We subsample over source experiments by randomized source superimpositions via the action of a source-encoding matrix sub E ∈ RNs ×Ns , which has randomized Gaussian-distributed entries. sub Ns is the number of simultaneous sources, Nssub  Ns . However, this dimensionality reduction introduces source crosstalks into the migrated image (Herrmann and Li, 2012). 2.4.1. Imaging with compressive sensing and sparsity-promoting Compressed sensing (CS) theory (Donoho, 2006) has shown that if signals’ energy is concentrated in a few large coefficients, these

signals can be recovered from incomplete and inaccurate measurements by solving sparsity-promoting problems (Candès et al., 2006b; Chai et al., 2017, 2016a, 2014, 2016b; Herrmann, 2010; Herrmann et al., 2008). CS determines conditions for which we can recover sparse   vectors x from measurements Ax = b, where A ∈ Rn×N or Cn×N ,   n×1 n×1 b ∈ R or C , with n  N. When x can be accurately approximated by its k-largest coefficients, stable recovery is possible for certain matrices A as long as n  k log(N/n) (Candès et al., 2006b). CS uses randomization to turn subsampling-related interferences (e.g., aliasing and source crosstalks) into relatively ‘harmless’ noise (Herrmann and Li, 2012). Therefore, the artifacts caused by the dimensionality reduction can be significantly mitigated by promoting sparsity of the medium perturbation dm. To maximally exploit the sparsity, we can employ some transforms (e.g., curvelet, seislet) that afford a sparse representation of the target solution (Candès et al., 2006a; Dutta, 2016; Dutta et al., 2016, 2017, 2015; Dutta and Schuster, 2015; Fomel and Liu, 2010; Herrmann et al., 2008). As an illustration, Fig. 2 shows a comparison between the decay of coefficients in the physical domain and the curvelet domain. A faster decay of the curvelet coefficients is evident, which indicates that the curvelet domain offers a better sparse representation of the medium perturbation dm than the physical domain. We therefore obtain the migrated image with transform-domain sparsity promotion from small subsets of simultaneous-source experiments. Specifically, we solve

 ∇ F m0 , Q, yj , Wj S∗ x = ddj ,

(19)

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

91

Fig. 10. (a) Shortest-offset section of 50% random noise added visco-acoustic data. Panel (b) displays the migrated image for the noisy visco-acoustic data shown in panel (a). Fig. 9. For Q true = 20, the migrated images in cases of (a) Q used = 10 < Q true and (b) Q used = 100 > Q true .

where S∗ denotes the adjoint of the sparse transform operator S. The function of S∗ is to apply the inverse sparse transform. Then, x is the transform domain coefficients of the migrated image dm, i.e., dm = S∗ x. If no sparse transform is employed, S = I.

that partly compensates for attenuation of seismic waves, and is defined by PQ = diag (q) ,

(22)

where 2.4.2. Preconditioning the imaging problem To accelerate the convergence rate of iterative migration, inspired by Herrmann et al. (2009), we apply a depth-preconditioning operator Pdepth , which consists of diagonal scaling that compensates for spherical spreading of seismic waves, and is defined by Pdepth

√  = diag z ⊗ INx ×Nx ,

where zl∈[1 solve

Nz ]

(20)

= lDz, Dz is the depth sampling interval. Then, we

 ∇ F m0 , Q, yj , Wj Pdepth S∗ x = ddj .

(21)

Since we are dealing with the Q filtering effects, we further design a Q-preconditioning operator PQ , which consists of diagonal scaling

q[lz +(mx −1)Nz ] =

lz k=1

1 , Q(k, mx )

(23)

and then,

q [lz +(mx −1)Nz ] + 1, −1 q[lz +(mx −1)Nz ] = pmax Q max(q)

(24)

is a user-specified where lz ∈ [1Nz ], mx ∈ [1Nx ], max(q) = 0, and pmax Q parameter controlling the maximal weight given to the deep reflectors. The Q-preconditioning operator PQ is functionally similar to the depth-preconditioning operator Pdepth , which also gives the lower reflectors more weights. Then, we solve

 ∇ F m0 , Q, yj , Wj Pdepth PQ S∗ x = ddj .

(25)

92

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

Fig. 11. A sedimentary part of the Sigsbee2B model, (a) true velocity model v, (b) the migration used background velocity model v0 , (c) true Q model, (d) the migration used Q model.

To facilitate further problem formulation, we cast the imaging problem in the canonical form of solving a linear system, i.e., Ax = b, where ⎡

⎤ ⎡ .. .. . ⎢ ⎥ ⎢ .

⎥ ⎢ ⎢ ∗ ⎥ dd A=⎢ ⎢ ∇ F m0 , Q, yj , Wj ⎥ Pdepth PQ S , b = ⎢ ⎣ j ⎣ ⎦ .

.. .

..

⎤ ⎥ ⎥ ⎥ , j ∈ Y. ⎦

attention recently because of its efficiency and simplicity in solving

1 -regularized problems. It is especially useful for problems in which the linear measurements are slow and expensive to obtain (Lorenz et al., 2014b). The LB method solves a regularized version of problem (27), i.e.,

(26) min x

lx1 +

1 x22 subject to Ax = b, 2

(28)

The final migrated image is obtained by dm = Pdepth PQ S∗ x. 2.4.3. Solution of Q-compensated imaging problem Sparse solutions can be obtained via 1 -norm minimization (Chen et al., 2001), i.e., min x1 subject to Ax = b. x

(27)

Some 1 -norm solvers developed in the CS filed (e.g., SPG1 , van den Berg and Friedlander, 2008) are rather complicated that it is difficult to implement them on industry problems. Following Herrmann et al. (2015), Yang et al. (2016), Chai et al. (2016b), we adopt a relatively simple linearized Bregman (LB) method (Cai et al., 2009; Lorenz et al., 2014b) to solve problem (25). The LB method has received a lot of

with regularization parameter l > 0. Lorenz et al. (2014a) proved the convergence of the LB method. Goldstein and Osher (2009) discussed the benefits of the Bregman iteration for 1 -regularized problems. It was shown in Yin (2010) that the solution to problem (28) is also a solution to problem (27) as long as l is chosen large enough. Comparison of the LB method for LSM in Herrmann et al. (2015) shows a performance that rivals and even exceeds the performance of state-of-the-art solver SPG1 . To tackle the presence of noise, we consider the following problem

min lx1 + x

1 x22 subject to Ax − b ≤ s, 2

(29)

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

where s is a user-specified parameter that allows for noise in the data. The choice of the norm of Ax − b is dictated by the noise characteristics, e.g., the 2 -norm is appropriate for Gaussian noise, the ∞ -norm for impulsive noise, and the 1 -norm for uniformly distributed noise (Lorenz et al., 2014a). The projections Ps (Ax − b) are designed to address the presence of noise, which involve orthogonal projections onto the respective norm-balls. For the 2 -norm,

  Ps (Ax − b) = max 0, 1 − s / Ax − b 2 (Ax − b),

(30)

and the ∞ -norm, Ps (Ax − b) = shrink(Ax − b, s),

(31)

where the component-wise soft shrinkage function shrink( • ) is defined as (Lorenz et al., 2014b) shrink(x, s) = max(|x| − s, 0)sign(x),

(32)

and fast algorithms are available for projections onto the 1 -normball (Duchi et al., 2008). Algorithm 1 gives a pseudo-code of the proposed stable and efficient Q-compensated imaging method. Lines 6–7 of Algorithm 1 explain the LB iterations in detail. Algorithm 1. (Stable and efficient Q-compensated imaging with compressive sensing, sparsity-promoting, and preconditioning).

The dynamic stepsize tk (Algorithm 1 line 6) is defined by (Lorenz et al., 2014b)

2

2 tk = Ps (Ak xk − bk ) 2 / A∗k Ps (Ak xk − bk ) 2 .

(33)

Generally, the loop can be stopped if the relative residual Ak xk − bk  /  bk  or the number of iterations k reach the target value. For the threshold l (Algorithm 1 line 7), we can generally sort the transform coefficients in descending order and set l to be the value corresponding to a user-specified percentage (e.g., 5%) of the total number of the transform coefficients. We can also typically set l as a fraction of the maximum value of the solution calculated at the first iteration, as Herrmann et al. (2015) did for sparse acoustic LSM. 2.4.4. Stability of Q-compensated imaging During each iteration, we use the Born operator for forward modeling and the adjoint operator for back propagating the residual

93

wavefields (see Algorithm 1 line 6). The adjoint operator A∗ is to calculate the complex conjugate transpose of the modeling operator A. As the amplitude attenuation is related to the real part of the Q filtering operator, the adjoint operator A∗ does not reverse the sign of the amplitude attenuation term, then the amplitude will not be exponential growth. Constructing the attenuation-compensated operators by reversing the sign of attenuation operator (e.g., Li et al., 2016; Wang, 2008; Zhang et al., 2013; Zhu et al., 2014) is avoided. The proposed method (Q-LSRTM using a sparsity formulation) is always stable. The principle of the proposed method circumnutating the intrinsic instability associated with Q-compensation is consistent with that described in Chai et al. (2017, 2016a, 2014) and Dutta and Schuster (2014). 2.4.5. SVD analysis of imaging operator As the iterative migration is using the pseudo-inverse of the forward operator to improve the accuracy of the migrated image, and the condition number is a measure of how close a matrix is to being singular (not invertible), SVD analysis of the forward operator can reveal some details of the visco-acoustic imaging/inversion as suggested by Kamei and Pratt (2013). The condition number j(A) of a matrix A is defined as j(A) = kmax (A)/kmin (A),

(34)

where kmax (A) and kmin (A) are maximal and minimal singular values of A, respectively. If j(A) is not too much larger than one, A is wellconditioned, which means its inverse can be computed with good accuracy. If j(A) is very large, A is said to be ill-conditioned. We attempt SVD analysis of the imaging operator with a small model (shown in Fig. 3a) for limited frequencies. The size of the model is a 134 × 267 grid, with a grid spacing of dx = dz = 7.5 m. A constant velocity 1500 m/s is used to generate the background model m0 . There is no attenuation in the first layer. We put 267 co-located sources and receivers with a 7.5 m lateral grid spacing at a depth of 15 m. Here “co-located” implies that sources and receivers are at the same grid locations. We use a Ricker wavelet with a peak frequency of 15 Hz, and the peak of the wavelet is at 0.1 s. We record for 4 s, which in the frequency domain yields 156 discretized frequencies up to 40 Hz (roughly the highest frequency that can be used for a 7.5 m grid spacing to avoid numerical dispersion). A total number of 144 frequencies in the range of [3 40] Hz are used for the imaging process. Fig. 3b shows the variations of j(A) with Q, for which the Q value of the rest of the layers is set to be 10, 20, 30, 40, 50, 60, 70, 80, 100, 150, 200, and +∞, respectively. We have two major findings deduced from Fig. 3b. First, attenuation causes j(A) to increase, making the imaging problem more ill-conditioned. It indicates that the visco-acoustic imaging problem converges slower than the acoustic imaging problem. The stronger the attenuation, the slower the convergence rate, the more the iterations needed. Second, the preconditioning operator evidently reduces j(A), making the imaging problem less ill-conditioned. It supports that the preconditioning strategies can accelerate the convergence rate of iterative migration. 3. Numerical experiments The following experiments are designed to check the proposed visco-acoustic imaging method avoiding the intrinsic instability associated with inverse Q-compensation, and to verify the benefits of employing the depth- and Q-preconditioning operators, the findings deduced from SVD analysis of imaging operator, and the robustness of the method proposed. Synthetic frequency-domain data were computed using the frequency-domain modeling method, followed by Fourier synthesis to create time-domain data. For simulating different signal-to-noise ratio (SNR) data, random noise (e.g., SNR = 2)

94

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

Fig. 13. With the same number of iterations, (a) acoustic LSM image for the acoustic data, (b) visco-acoustic LSM image for the visco-acoustic data.

was randomly generated as 50% of the 2 -norm of the noise-free data and then added to the noise-free data. Q true and Qused denote the true and the used Q models in the imaging process, respectively. When we want to highlight one influencing factor, we kept all the other parameters being the same. 3.1. Graben model example

Fig. 12. For visco-acoustic data, (a) acoustic LSM image (i.e., without Qcompensation), (b) visco-acoustic LSM image generated by the proposed method. (c) Differences between (a) and (b). White lines — without Q compensation. Black lines — with Q compensation.

The experiment setup is mostly consistent with that given in the SVD analysis subsection, except that Q value of the rest of the layers is specified to be 20. Fig. 4a shows the zero-offset section of the generated visco-acoustic data. The amplitude attenuation and phase dispersion effects are evident in Fig. 4a. The high-frequency components of the visco-acoustic data are more attenuated than those of the acoustic data, and the dominant frequency band moves toward to low-frequency (see Fig. 4b).

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

Fig. 14. Migrated image (a) without and (b) with compressive sensing and sparsitypromoting. (c) Differences between (a) and (b).

95

Fig. 15. (a) Migrated image without preconditioning. (b) Migrated image with depthand Q-preconditioning. (c) Differences between (a) and (b).

96

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

Fig. 3a provides the reference result and all the following LSM images are plotted in the same scale as Fig. 3a. The migrated image obtained by applying the adjoint operator (Fig. 5a) does not have the correct amplitude, and the graben structure is unrecovered. Fig. 5b displays the acoustic LSM image for the visco-acoustic data (i.e., no Q-compensation). The first event is correctly recovered because there is no attenuation. The graben structure is migrated with weak amplitude and a slight mispositioning in its location. All the reflectors of the migrated image produced by the proposed method are recovered with accurate amplitudes and positions (see Fig. 5c). There is no noise caused by the intrinsic instability associated with conventional inverse Q-compensation. Comparison of Fig. 5b–d reveals the significance of accomplishing Q-compensation during seismic migration process. With the same number of iterations, Fig. 6a represents the acoustic LSM image for the acoustic data, and Fig. 6b denotes the visco-acoustic LSM image for the visco-acoustic data. For producing the images in Fig. 6, we use 6 randomly selected frequencies and 23 simultaneous sources, without using preconditioning operators. For the acoustic data, the graben reflector of the acoustic LSM image (Fig. 6a) is close to the true solution. However, for the Q-filtered data, the graben reflector of the visco-acoustic LSM image (Fig. 6b) is far from the true one. This verifies that the acoustic imaging problem converges faster than the visco-acoustic case, as attenuation makes j(A) increasing and the imaging problem more ill-conditioned. Therefore, solving the visco-acoustic imaging problem needs more iterations than the acoustic case. Fig. 7 shows the migrated image without and with the compressive sensing and sparsity-promoting strategies. For producing the images in Fig. 7, we use 6 randomly selected frequencies and 1 simultaneous source, with the depth- and Q-preconditioning operators. The dimensionality reduction introduces severe source crosstalks into the migrated image (Fig. 7a). However, those artifacts (shown in Fig. 7c) are significantly mitigated via CS and sparsity-promoting strategies (see Fig. 7b). This underlines the importance of sparsity promotion, which enable us to mitigate the noisy source crosstalks while still benefiting from computational gains related to the reduction of the number of PDE to be solved. Fig. 8 illustrates the advantages of employing the depth- and Q-preconditioning operators. Without the preconditioning strategies, we cannot recover the strongly attenuated weak events with the visco-acoustic LSM method despite after many iterations (see Fig. 8a). However, we obtain better results after adopting the depth preconditioning operator (see Fig. 8b). Then we get even better results after utilizing both depth- and Q-preconditioning strategies (see Fig. 8c). This verifies that the preconditioning techniques make the imaging problem less ill-conditioned and accelerate the convergence rate of LSM. The graben reflector of the migrated image with Q used < Q true (Fig. 9a) is over-compensated. On the contrary, the graben reflector of the migrated image with Q used > Q true (Fig. 9b) is undercompensated. Fig. 10 serves to verify the robustness of the proposed method for noisy data. The migration results mostly resemble those of the true model and remain at acceptable recovery levels. This indicates that the proposed Q-compensated imaging method is satisfactorily stable in the presence of severe noise for this structurally simple model. 3.2. Sigsbee2B model example To evaluate the performance of the proposed methodology for a more realistic setup, we verify the method with a more complex earth model (Fig. 11a). The model is cropped from the Sigsbee2B model that contains abundant sedimentary layers. The background model v0 (Fig. 11b) is obtained by smoothing the true model. The size of the model is a 268 × 401 grid. The grid spacing is 7.5 m. The Q model used for modeling the visco-acoustic data (shown in Fig. 11c) is determined from velocity in km/s by using the following empirical

Fig. 16. (a) Section of 50% random noise added visco-acoustic data at shot number 101 at a distance of 1.5 km. Panel (b) displays the migrated image for the noisy viscoacoustic data shown in panel (a).

relation (Li et al., 2016) Q ≈ 14v2.2 .

(35)

Q used (shown in Fig. 11d) is also calculated with Eq. (35), but using the background velocity v0 . For the overlying water layer, Q → +∞. There are 201 co-located sources and receivers with a 15 m lateral grid spacing at a depth of 7.5 m. We use a Ricker wavelet with a peak frequency of 15 Hz, and the peak of the wavelet is at 0.2 s. A total number of 144 frequencies in the range [3 40] Hz are used for the imaging process. For Q-filtered data, Fig. 12a shows the migrated image without Q-compensation, and Fig. 12b shows the Q-compensated image generated by the proposed method. All the following migrated images are plotted in the same scale. On both edges of the migrated images, the recovered amplitude decreases because of the weak illumination. The reflectors in Fig. 12a are migrated with weak amplitudes (especially for the deep reflectors), because of the uncompensated attenuation and dispersion effects. However, the reflectors in Fig. 12b are recovered with more balanced amplitudes and higher resolution.

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

97

Again, there is no noise caused by the intrinsic instability associated with conventional inverse Q compensation. Fig. 12 further confirms the significance of achieving Q-compensation during imaging. Fig. 13a shows the acoustic LSM image for the acoustic data, and Fig. 13b shows the visco-acoustic LSM image for the visco-acoustic data. Here we use 6 randomly selected frequencies and all sources, without using preconditioning operators to isolate the preconditioning effect. The migration results in Fig. 13a are much better than that shown in Fig. 13b. This confirms that the acoustic LSM converges faster than the visco-acoustic LSM. Therefore, the visco-acoustic LSM needs more iterations than the acoustic LSM. The artifacts caused by the dimensionality reduction are evident in Fig. 14a (see also Fig. 14c). However, those crosstalks shown in Fig. 14c are significantly mitigated by CS and sparsity-promoting strategies (see Fig. 14b). It indicates that sparsity-promoting enable us to achieve a significant speed-up without notable compromise on the image quality with the subsampling strategy. Fig. 15 verifies the benefits of employing the depth- and Q-preconditioning operators. The resolution of the deeper weak events are improved with the preconditioning strategies. Fig. 16 shows the migrated images for Q-filtered data with 50% noise added. The migration results remain at acceptable recovery levels. Note that for all the migrated images, Qused is a smoothed Q model instead of the true Q model, which indicates that relatively imprecise Q estimates can still produce adequate seismic imaging.

4. Discussion As demonstrated in Figs. 6 and 13, the visco-acoustic LSM needs more iterations than the acoustic LSM. In conventional visco-acoustic imaging approaches, researchers usually work with all shots and frequencies, but can therefore only afford a few iterations. With the proposed approach, we work with significantly smaller data subsets and can therefore afford more iterations since each iteration is relatively cheaper. Fig. 17 indicates that more iterations with cheaper per-iteration cost lead to better Q-compensated images. As shown in Figs. 8 and 15, the depth- and Q-preconditioning operators significantly accelerate the convergence rate of visco-acoustic LSM. Therefore, other helpful preconditioners should be investigated. For example, the source-side illumination (Plessix and Mulder, 2004) may be used as a diagonal preconditioner. Same as most of wave-equation based imaging techniques (e.g., Herrmann and Li, 2012; Nemeth et al., 1999; Plessix and Mulder, 2004), the proposed approach hinges on the availability of a smooth background velocity model that accurately captures the kinematics of reflected waves. Also, same as almost all of the Q-compensated imaging schemes (e.g., Dutta and Schuster, 2014; Li et al., 2016; Wang, 2008; Wang and Guo, 2004; Zhang et al., 2013; Zhu et al., 2014), the method will depend on how good the preceding estimation of the Q model. The proposed methodology can also be adjusted for the time-domain Q-compensated imaging (Sun et al., 2016; Zhu and Sun, 2017).

5. Conclusions

Fig. 17. The visco-acoustic LSM image generated by the proposed method with (a) all 201 shots and 144 frequencies but 10 iterations, (b) 67 shots (i.e., 13 shot-subsampling) 1 and 6 frequencies (i.e., 24 frequency-subsampling) but 100 iterations. (c) Differences between (a) and (b).

We proposed a stable and efficient Q-compensated imaging method with compressive sensing, sparsity-promoting, and preconditioning. The stability is achieved by using the linearized Born operator for forward modeling and the adjoint operator for back propagating the residual wavefields. Constructing the attenuationcompensated operators by reversing the sign of attenuation operator, which arises the numerical instability, is circumvented. To reduce the computational cost that is proportional to the number of PDE to be solved, we first subsample over both frequencies and source experiments, then we employ a depth-preconditioning

98

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99

operator and a Q-preconditioning operator to accelerate the convergence rate of iterative migration. Numerical experiments verify that the proposed approach can avoid the intrinsic instability associated with inverse Q-compensation, and it is robust in the presence of severe noise. Attenuation causes the condition number of the migration operator to increase, which consequently makes the imaging problem more ill-conditioned. The visco-acoustic imaging problem converges slower than the acoustic case. The stronger the attenuation, the more the iterations needed. The preconditioning strategy evidently reduces the condition number of the migration operator, making the imaging problem less ill-conditioned and accelerating the convergence rate of LSM. Acknowledgments This work was supported by the National Key Basic Research Development Program (grant no. 2013CB228600), the National Natural Science Foundation of China Young Scholar Program (grant no. 41304042, 41604092), and the Hubei Subsurface Multi-scale Imaging Key Laboratory (China University of Geosciences) Program (grant no. SMIL-2017-01). X. Chai greatly thanks all the people at University of British Columbia (UBC)-Seismic Laboratory for Imaging and Modeling (SLIM) for their beneficial discussion and generous help. The authors are grateful to Co-editor-in-chief J. Behura, reviewers T. Bai and G. Dutta for their critical comments and useful suggestions, which significantly improved the quality of this paper.

References Baysal, E., Kosloff, D.D., Sherwood, J.W.C., 1983. Reverse time migration. Geophysics 48, 1514–1524. Berenger, J., 1994. A perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys. 114, 185–200. Cai, J.-F., Osher, S., Shen, Z., 2009. Linearized Bregman iterations for compressed sensing. Math. Comput. 78, 1515–1536. Candès, E., Demanet, L., Donoho, D., Ying, L., 2006. Fast discrete curvelet transforms. Multiscale Model. Simul. 5, 861–899. Candès, E.J., Romberg, J.K., Tao, T., 2006. Stable signal recovery from incomplete and inaccurate measurements. Commun. Pure Appl. Math. 59, 1207–1223. Carcione, J.M., Kosloff, D., Kosloff, R., 1988. Wave propagation simulation in a linear viscoacoustic medium. Geophys. J. Int. 93, 393–401. Chai, X., Wang, S., Tang, G., 2017. Sparse reflectivity inversion for nonstationary seismic data with surface-related multiples: numerical and field-data experiments. Geophysics 82, R199–R217. Chai, X., Wang, S., Wei, J., Li, J., Yin, H., 2016. Reflectivity inversion for attenuated seismic data: physical modeling and field data experiments. Geophysics 81, T11–T24. Chai, X., Wang, S., Yuan, S., Zhao, J., Sun, L., Wei, X., 2014. Sparse reflectivity inversion for nonstationary seismic data. Geophysics 79, V93–V105. Chai, X., Yang, M., Witte, P., Wang, R., Fang, Z., Herrmann, F., 2016. A linearized Bregman method for compressive waveform inversion. SEG Technical Program Expanded Abstracts 2016. pp. 1449–1454. Chattopadhyay, S., McMechan, G.A., 2008. Imaging conditions for prestack reverse time migration. Geophysics 73, S81–S89. Chen, S.S., Donoho, D.L., Saunders, M.A., 2001. Atomic decomposition by basis pursuit. SIAM Rev. 43, 129–159. Chen, Z., Cheng, D., Feng, W., Wu, T., 2013. An optimal 9-point finite difference scheme for the Helmholtz equation with PML. Int. J. Numer. Anal. Model. 10, 389–410. Donoho, D.L., 2006. Compressed sensing. IEEE Trans. Inf. Theory 52, 1289–1306. Duchi, J., Shalev-Shwartz, S., Singer, Y., Chandra, T., 2008. Efficient projections onto the 1 -ball for learning in high dimensions. Proc. 25th Int. Conf. Mach. Learn. 78, 272–279. Dutta, G., 2016. Sparse least-squares reverse time migration using seislets. J. Appl. Geophys. 136, 142–155. Dutta, G., Agut, C., Giboli, M., Williamson, P., 2016. Least-squares reverse time migration with radon preconditioning. SEG Technical Program Expanded Abstracts. pp. 4198–4203. Dutta, G., Giboli, M., Agut, C., Williamson, P., Schuster, G.T., 2017. Least-squares reverse time migration with local radon-based preconditioning. Geophysics 82, S75–S84. Dutta, G., Giboli, M., Williamson, P., Schuster, G.T., 2015. Least-squares reverse time migration with factorization-free priorconditioning. SEG Technical Program Expanded Abstracts. pp. 4270–4275. Dutta, G., Schuster, G.T., 2014. Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation. Geophysics 79, S251–S262.

Dutta, G., Schuster, G.T., 2015. Sparse least-squares reverse time migration using seislets. SEG Technical Program Expanded Abstracts. pp. 4232–4237. Etgen, J., Gray, S.H., Zhang, Y., 2009. An overview of depth imaging in exploration geophysics. Geophysics 74, WCA5–WCA17. Fomel, S., Liu, Y., 2010. Seislet transform and seislet frame. Geophysics 75, V25–V38. Futterman, W.I., 1962. Dispersive body waves. J. Geophys. Res. 67, 5279–5291. Goldstein, T., Osher, S., 2009. The split Bregman method for 1 -regularized problems. SIAM J. Imag. Sci. 2, 323–343. Herrmann, F.J., 2010. Randomized sampling and sparsity: getting more information from fewer samples. Geophysics 75, WB173–WB187. Herrmann, F.J., Brown, C.R., Erlangga, Y.A., Moghaddam, P.P., 2009. Curvelet-based migration preconditioning and scaling. Geophysics 74, A41–A46. Herrmann, F.J., Li, X., 2012. Efficient least-squares imaging with sparsity promotion and compressive sensing. Geophys. Prospect. 60, 696–712. Herrmann, F.J., Moghaddam, P., Stolk, C.C., 2008. Sparsity- and continuity-promoting seismic image recovery with curvelet frames. Appl. Comput. Harmon. Anal. 24, 150–173. Herrmann, F.J., Tu, N., Esser, E., 2015. Fast “online” migration with compressive sensing. EAGE Annual Conference Proceedings. Hicks, G.J., Pratt, R.G., 2001. Reflection waveform inversion using local descent methods: estimating attenuation and velocity over a gas-sand deposit. Geophysics 66, 598–612. Kamei, R., Pratt, R.G., 2013. Inversion strategies for visco-acoustic waveform inversion. Geophys. J. Int. 194, 859–884. Kazemi, N., Sacchi, M.D., 2015. Block row recursive least-squares migration. Geophysics 80, A95–A101. Kjartansson, E., 1979. Constant-Q wave propagation and attenuation. J. Geophys. Res. 84, 4737–4748. Knopoff, L., 1964. Q. Rev. Geophys. 2, 625–660. Kolsky, H., 1956. The propagation of stress pulses in viscoelastic solids. Philos. Mag. 1, 693–710. Levin, S.A., 1984. Principle of reverse-time migration. Geophysics 49, 581–583. Li, Q., Zhou, H., Zhang, Q., Chen, H., Sheng, S., 2016. Efficient reverse time migration based on fractional Laplacian viscoacoustic wave equation. Geophys. J. Int. 204, 488–504. Liu, H., Anderson, D.L., Kanamori, H., 1976. Velocity dispersion due to anelasticity; implications for seismology and mantle composition. Geophys. J. R. Astron. Soc. 47, 41–58. Lorenz, D.A., Schöpfer, F., Wenger, S., 2014. The linearized Bregman method via split feasibility problems: analysis and generalizations. SIAM J. Imag. Sci. 7, 1237–1262. Lorenz, D.A., Wenger, S., Schöpfer, F., Magnor, M., 2014. A sparse Kaczmarz solver and a linearized Bregman method for online compressed sensing. IEEE International Conference on Image Processing. pp. 1347–1351. Marfurt, K., 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics 49, 533–549. Mittet, R., 2007. A simple design procedure for depth extrapolation operators that compensate for absorption and dispersion. Geophysics 72, S105–S112. Mittet, R., Sollie, R., Hokstad, K., 1995. Prestack depth migration with compensation for absorption and dispersion. Geophysics 60, 1485–1494. Morozov, I., Ahmadi, A.B., 2015. Taxonomy of Q. Geophysics 80, T41–T49. Nemeth, T., Wu, C., Schuster, G.T., 1999. Least-squares migration of incomplete reflection data. Geophysics 64, 208–221. Plessix, R.-E., Mulder, W.A., 2004. Frequency-domain finite-difference amplitude-preserving migration. Geophys. J. Int. 157, 975–987. Pratt, G., Shin, C., Hicks, 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophys. J. Int. 133, 341–362. Prieux, V., Brossier, R., Operto, S., Virieux, J., 2013. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field. Part 1: imaging compressional wave speed, density and attenuation. Geophys. J. Int. 194, 1640–1664. Quan, Y., Harris, J.M., 1997. Seismic attenuation tomography using the frequency shift method. Geophysics 62, 895–905. Song, Z.-M., Williamson, P.R., 1995. Frequency-domain acoustic-wave modeling and inversion of crosshole data: part 1 – 2.5-D modeling method. Geophysics 60, 784–795. Song, Z.M., Williamson, P.R., Pratt, R.G., 1995. Frequency-domain acoustic-wave modeling and inversion of crosshole data: part II-inversion method, synthetic experiments and real-data results. Geophysics 60, 796–809. Sun, J., Fomel, S., Zhu, T., Hu, J., 2016. Q-compensated least-squares reverse time migration using low-rank one-step wave extrapolation. Geophysics 81, S271–S279. Sun, J., Zhu, T., Fomel, S., 2015. Viscoacoustic modeling and imaging using low-rank approximation. Geophysics 80, A103–A108. Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics 49, 1259–1266. Tarantola, A., 1988. Theoretical background for the inversion of seismic waveforms including elasticity and attenuation. Pure Appl. Geophys. 128, 365–399. van den Berg, E., Friedlander, M.P., 2008. Probing the Pareto frontier for basis pursuit solutions. SIAM J. Sci. Comput. 31, 890–912. van Leeuwen, T., Herrmann, F.J., 2013. Mitigating local minima in full-waveform inversion by expanding the search space. Geophys. J. Int. 195, 661–667. Virieux, J., Operto, S., 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics 74, WCC1–WCC26. Wang, Y., 2006. Inverse Q-filter for seismic resolution enhancement. Geophysics 71, V51–V60.

X. Chai et al. / Journal of Applied Geophysics 145 (2017) 84–99 Wang, Y., 2008. Inverse-Q filtered migration. Geophysics 73, S1–S6. Wang, Y., Guo, J., 2004. Seismic migration with inverse Q filtering. Geophys. Res. Lett. 31, 1–4. Wu, S., Wang, Y., Zheng, Y., Chang, X., 2015. Limited-memory BFGS based least-squares pre-stack kirchhoff depth migration. Geophys. J. Int. 202, 738–747. Yang, M., Witte, P.A., Fang, Z., Herrmann, F.J., 2016. Time-domain sparsity-promoting least-squares migration with source estimation. 86th Annual International Meeting, SEG, Expanded Abstracts. Yin, W., 2010. Analysis and generalizations of the linearized Bregman method. SIAM J. Imag. Sci. 3, 856–877. Zhang, C., Ulrych, T.J., 2010. Refocusing migrated seismic images in absorptive media. Geophysics 75, S103–S110.

99

Zhang, J., Wapenaar, K., 2002. Wavefield extrapolation and prestack depth migration in anelastic inhomogeneous media. Geophys. Prospect. 50, 629–643. Zhang, J., Wu, J., Li, X., 2013. Compensation for absorption and dispersion in prestack migration: an effective Q approach. Geophysics 78, S1–S14. Zhu, T., Carcione, J.M., Harris, J.M., 2013. Approximating constant-Q seismic propagation in the time domain. Geophys. Prospect. 61, 931–940. Zhu, T., Harris, J.M., Biondi, B., 2014. Q-compensated reverse-time migration. Geophysics 79, S77–S87. Zhu, T., Sun, J., 2017. Viscoelastic reverse time migration with attenuation compensation. Geophysics 82, S61–S73.