During propagating inside the Earth, seismic wave energies are gradually converted into heat due to the presence of internal friction (Knopoff, 1964;Walcott, 1970). Therefore, seismograms are affected by the attenuation-associated phase dispersion and amplitude dissipation. Measuring such waveform distortions is an important way to constrain intrinsic attenuation in the Earth. To complement velocity structures, seismic attenuation provides important information to study the fluid properties of hydrocarbon reservoirs in exploration seismology (Müller et al., 2010), as well as temperature distributions, partial melting, and water content within the crust and mantle in earthquake seismology (Durek & Ekström, 1996;Romanowicz & Mitchell, 2007;Stachnik et al., 2004;Wilcock et al., 1992;Zhu et al., 2013).
Early estimates of subsurface distributions of velocity and attenuation are based on ray-based tomography (Brzostowski & McMechan, 1992;Haberland & Rietbrock, 2001;Nolet, 1987;Pozgay et al., 2009;Quan & Harris, 1997;Roth et al., 1999;Wilcock et al., 1992;Wu & Toksöz, 1987). Increasing computational capability makes it possible to directly solve the wave equation using numerical methods, such as finite differences, finite elements, and spectral elements. Using these numerical solvers, the inverse problem can be formulated with the adjoint-state method to estimate (an)elastic parameters in the subsurface, which is known as full waveform inversion (FWI) (Liu & Gu, 2012;Pratt, 1999;Pratt et al., 1998;Tarantola, 1984Tarantola, , 1986Tarantola, , 2005;;Tromp et al., 2005). In exploration seismology, FWI has been successfully used to characterize detailed structures for hydrocarbon reservoirs (Ratcliffe et al., 2011;Sears et al., 2010;Sirgue et al., 2010;Virieux & Operto, 2009). In earthquake seismology, by fitting passive-source seismograms with calculated data, FWI has been used to constrain heterogeneous structures within the crust and mantle (Bozdaǧ et al., 2016;Tape et al., 2010; 10.1029/2019JB019129 Zhu et al., 2015), which provide important information to investigate tectonic processes and physical properties of the Earth's materials (Bozdaǧ et al., 2015;Chen et al., 2017;Fichtner & Villaseñor, 2015;Govers & Fichtner, 2016;Tape et al., 2009;Zhu et al., 2012;Zhu & Tromp, 2013).
In the frequency domain, incorporating attenuation into FWI can be implemented using a complex-valued velocity (Aki & Richards, 1980;Kamei & Pratt, 2013;Liao & McMechan, 1996). Solving the Helmholtz equation for the forward and adjoint wavefields requires large computer memory costs, which is still challenging for large 3-D models with current computational capability (Operto et al., 2007). The standard linear solid (SLS) is a typical rheology model for incorporating attenuation into the time domain wave equation (Carcione, 2007), in which case the quality factor Q has to be converted to the stress and strain relaxation times. When multiple attenuation mechanisms are superimposed to approximate a frequency-independent Q within a certain frequency band, it is difficult to directly estimate the corresponding stress and strain relaxation times using FWI; the relaxation times with different reference frequencies might depend on each other, and the FWI results are nonunique. To mitigate these problems in the SLS methods, Fichtner and van Driel (2014) present a new method for modeling frequency-dependent and independent Q in the time domain and computing the corresponding Fréchet derivatives. Alternatively, Zhu and Harris (2014) simplify the fractional time derivative in the constant-Q equation (Kjartansson, 1979) using a fractional Laplacian operator and derive a practical viscoacoustic wave equation. They subsequently apply this equation to compensate attenuation effects in reverse-time migration (RTM) (Zhu, 2014;Zhu et al., 2014), reflectivity inversion (Sun et al., 2016), and FWI (Xue et al., 2018). Based on a modified constant-Q wave equation (Xing & Zhu, 2018), Xing and Zhu (2019) derive the corresponding Fréchet kernels and misfit gradients for traveltime, amplitude, and waveform measurements. Yang et al. (2016) give a systematic review of the viscoelastic wave equations and the corresponding applications in FWI. Using a constant relaxation time in the SLS wave equation, Fabien-Ouellet et al. (2017) define a parameter 휏 to describe the attenuation level in viscoelastic FWI, which is used to monitor CO 2 injection. Karaoglu and Romanowicz (2017) test three misfit functionals for 3-D imaging of shear attenuation in the Earth's upper mantle at a global scale and observe that the waveform misfit is less robust than waveform envelope and spectral amplitude ratio for constraining Q distribution. They subsequently develop a hybrid full-waveform inversion method for mapping global shear wave attenuation (Karaoglu & Romanowicz, 2018).
Using a second-order polynomial to approximate the dispersion term, and a pseudodifferential operator to approximate the dissipation term, Yang and Zhu (2018a) derive a time domain complex-valued wave equation for modeling wave propagation in viscoacoustic media. Because an imaginary term is introduced into the dispersion approximation, this wave equation is complex-valued, and their real and imaginary parts are coupled during wave propagation. In comparison with traditional frequency-domain viscoacoustic wave equation and the time domain SLS equation, the new complex-valued wave equation has the following advantages. First, the dispersion and dissipation effects are separated, which enables us to compensate amplitude loss in migration by flipping the sign of the dissipation term and keeping the dispersion term unchanged as in the constant-Q method (Yang & Zhu, 2018b;Zhu et al., 2014). Second, Q is explicitly incorporated into the wave equation, which makes it easier to derive the misfit gradient with respect to Q in FWI than in the SLS-based methods. Third, it can be numerically solved using finite-difference time marching and Fourier transform, which has lower memory requirements compared with directly solving the Helmholtz equation in the frequency domain.
Based on this time domain, complex-valued wave equation, we develop an FWI framework to simultaneously estimate P wave velocity and attenuation distributions in this study. Using the Lagrange multiplier method, we derive the adjoint wave equation and sensitivity kernels. Numerical analysis for the Gauss-Newton Hessian demonstrates that there is strong interparameter crosstalk between the misfit gradients, especially for the leakage from P wave velocity to Q. To reduce these artifacts in viscoacoustic FWI applications, we utilize an optimized L-BFGS method (Liu & Nocedal, 1989;Zhu et al., 1997) to approximate the full Hessian, in which the diagonal Hessian is used as a preconditioner to accelerate convergence speed. In comparison with the conjugate gradient method, the proposed simultaneous inversion strategy can effectively reduce interparameter crosstalk and accurately estimate P wave velocity and attenuation distributions. In the following sections, we first give a brief review of the time domain, complex-valued viscoacoustic wave equation and then describe the simultaneous inversion method. Several synthetic examples are used 10.1029/2019JB019129 to illustrate the feasibility and robustness of the proposed viscoacoustic FWI framework for simultaneously mapping velocity and attenuation in the subsurface.
In the frequency domain, the dispersion and dissipation effects caused by intrinsic attenuation can be described using a complex-valued velocity (Aki & Richards, 1980):
where x denotes a location within the subsurface region of interest, 휔 is the angular frequency, and v(x) is the velocity at the reference frequency 휔 0 , which is set as 1 Hz in this study. Q(x) is the quality factor, and we assume it is frequency independent. i is the imaginary unit, and sgn is the sign function.
Using a second-order polynomial to approximate the dispersion term and a pseudodifferential operator to approximate the dissipation term, Yang and Zhu (2018a) derived the following viscoacoustic wave equation:
with
where p(x, t) is the pressure wavefield; 휌(x) is the density; f(t) is the source time function; x s is the source location; and ∇, ∇•, and ∇ 2 denote the gradient, divergence, and Laplacian operators, respectively. a, b, and c are the polynomial coefficients used in the dispersion approximation (J. Yang & Zhu, 2018a). They are related with the reference frequency 휔 0 and can be calculated by solving a regression problem within a certain frequency band as
For instance, if the reference frequency is set as 1 Hz and the fitting frequencies range from 1 to 100 Hz, we have a = 5.3343, b = -510.8381, and c = 2.0837e4. The terms including C 1 , C 2 , and C 3 in Equation 2 determine phase dispersion, and the term with C 4 describes amplitude dissipation. Since an imaginary unit is introduced in the dispersion approximation, Equation 2 is complex-valued, and the real and imaginary parts of the wavefields are coupled during propagation.
With the time domain complex-valued wave equation, the adjoint viscoacoustic wave equation and sensitivity kernels can be derived using the adjoint-state method (Liu & Tromp, 2006;Plessix, 2006;Tromp et al., 2005). The adjoint wave equation has the following expression (the detailed derivation is given in Appendix A):
where T denotes the duration time, x r is the receiver location, 훿 is the Kronecker delta function, and p † (x, t) is the adjoint wavefield. C 0 1 , C 0 2 , C 0 3 , and C 0 4 have the same expressions as in Equation 3, except replacing the true model parameters 휌, v, and Q with the background parameters 휌 0 , v 0 , and Q 0 . 푓 † (x, t) is the adjoint source and has the following form
10.1029/2019JB019129
where * denotes complex conjugate and d obs (x r ) and d syn (x r ) are the observed and predicted data, respectively. 휌 0 , v 0 , and Q 0 are the background density, P wave velocity, and quality factor, respectively. As in the forward wave equation, the adjoint wave equation is also complex-valued, and its energy is dissipated during back-propagation. The complex conjugate of the data residual is used as the adjoint source in back-propagation.
With the adjoint wave Equation 5, the P wave velocity and attenuation sensitivity kernels can be written as
where ℜ denotes the real part, and the relative model perturbations, that is, 훿 ln v = 훿v∕v 0 and 훿 ln Q = 훿Q∕Q 0 , are used in the derivation of the sensitivity kernels (see Appendix A). Equation 7illustrates that the velocity and Q sensitivity kernels are computed by convolving the forward and adjoint wavefields, followed by stacking all kernels over sources and receivers to calculate final misfit gradients. Since the energies of both forward and adjoint wavefields are attenuated, Q leads to double-damping in the gradients. In addition, the difference between P wave velocity and Q kernels is the coefficients of p(x, t) in the first term of Equation 7, which indicates that the radiation patterns for velocity and Q are similar. Therefore, it is not easy to separate velocity and Q crosstalk by only using different scattering-angle data.
Interparameter crosstalk is a critical problem in multiple-parameter inversion. Since different elastic model parameters, such as density and P and S wave impedances and velocities, have different scattering radiation patterns (Operto et al., 2013;Virieux & Operto, 2009;Zhou et al., 2015), utilizing different scattering-angle data to constrain these model parameters is a common strategy to mitigate interparameter crosstalk. But Malinowski et al. (2011) demonstrate that the partial derivatives of the pressure wavefield with respect to P wave velocity and Q have similar radiation patterns, and their crosstalk is difficult to separate by only using different scattering-angle information. Therefore, the introduction of the Hessian into viscoacoustic FWI is more important to reduce the crosstalk between velocity and Q than that among density, velocity, and impedance. In this section, we first present the Gauss-Newton Hessian based on the time domain complex-valued wave equation and then analyze the crosstalk between velocity and Q using a simple numerical example. The Gauss-Newton Hessian can be expressed as
with the elements as
where d syn (x s , x r , t) is the synthetic data recorded at the receiver location x r excited at the source location x s . 휕d s푦n ∕휕 ln v and 휕d s푦n ∕휕 ln Q are the partial derivative wavefields (Pratt et al., 1998) from a unit velocity (훿 ln v) and quality factor (훿 ln Q) perturbation, respectively. Physically, each element of H vv (x, x ′ ) represents the derivative of the misfit gradient for P wave velocity at location x ′ caused by a unit velocity perturbation at location x. H QQ (x, x ′ ) has a similar physical meaning except for the quality factor.
10.1029/2019JB019129 mean the derivatives of the misfit gradients for one model parameter caused by the other model parameter perturbations and therefore are related with interparameter crosstalk.
The partial derivative wavefields from velocity and Q perturbations can be derived using the Born approximation (the detailed derivation is given in Appendix B). With these partial derivative wavefields, the elements of the Gauss-Newton Hessian can be written as
with
where G(x, x r , t) is the Green's function from an image point x to the receiver location x r , and p(x s , x, t) is the forward wavefield propagating from the source location x s to x. Note that the Gauss-Newton Hessian describes the single scattering effect from an image point x to another point x ′ in terms of a given acquisition geometry and background models. Numerically, each column of the Gauss-Newton Hessian can be computed using a Born modeling (the last two terms involving G(x, x r , t) and P(x s , x, t) in Equation 10) and an adjoint-based reverse-time migration (the first two terms involving G * (x, x r , t) and P * (x s , x, t) in Equation 10); its computational cost is about 4 times that of a normal forward modeling. Therefore, the total cost for computing the whole Gauss-Newton Hessian is about 4N × N s × C sim , where N is the discrete node number for velocity and Q models, N s is the source number, and C sim is the cost of one normal forward simulation.
To analyze the interparameter crosstalk between P wave velocity (Vp) and Q, we compute the true Gauss-Newton Hessian for a homogeneous model with Vp = 3 km/s, Q = 100, and a constant density of 1.0 g/cm 3 . The model size is 51 × 51 with 20 m grid spacing. Twenty-five shots are evenly distributed along the left boundary, and 51 receivers are used to record the data for each shot (Figure 1a). Another shot-receiver group is deployed along the top and bottom boundaries. The Gauss-Newton Hessian is shown in Figure 1b, and the point-spread-functions (PSFs), which are the derivative of the misfit gradient for a point perturbation (Humphreys & Clayton, 1988;Yu et al., 2002) Figures 2a and 2b show the true perturbations of P wave velocity and Q for an uncorrelated-square model. The positive perturbations for velocity (훿v∕v 0 ) and quality factor (훿Q∕Q 0 ) are 33% and 100%, respectively. The corresponding misfit gradients are presented in Figures 2c and 2d. Although the P wave velocity gradient is blurred (Figure 2c) due to the finite-frequency effect, the square structure is imaged at the correct location.
In contrast, the Q gradient (Figure 2d) is similar to the velocity gradient, except with an opposite polarity. This negative gradient is a crosstalk artifact leaked from the P wave velocity, and the true Q perturbation is submerged in these artifacts. To reduce the crosstalk, the inverse of the truncated Gauss-Newton Hessian is applied to P wave velocity and Q gradients. The Hessian inverse is computed using the singular-value decomposition (SVD). The size of the Hessian matrix is 5,202. The optimized gradients using the truncated Hessian with the rank of 2000, 3000, and 5202 are shown in Figure 3. Using the rank of 2000, the truncated Gauss-Newton Hessian improves the spatial resolution for the P wave velocity gradient, but the Q perturbation is still missing. Increasing the rank to 3000, the finite-frequency artifacts in the P wave velocity gradient are further reduced, and the square structure for Q is recovered with blurred edges. Using the full rank, both P wave velocity and Q perturbations are recovered with high spatial resolution. This experiment illustrates that the crosstalk from velocity to Q models is strong, and it even submerges the true Q perturbations. Therefore, the Gauss-Newton or full Hessian should be introduced into FWI to reduce these interparameter artifacts, which facilitates to produce accurate velocity and Q updates. In addition, this experiment indicates that the two-parameter (Vp and Q) Gauss-Newton Hessian is not low-rank, and only using a few major ranks to compute the Hessian inverse is difficult to effectively reduce the crosstalk artifacts.
Directly computing the Gauss-Newton Hessian and its inverse involves a large number of wavefields simulations and is prohibitive for large-scale problems under current computational capability. In this section, we derive an approximate diagonal Hessian and use it as a preconditioner for the misfit gradients in viscoacoustic FWI. By setting x ′ = x in Equation 10, the Gauss-Newton Hessian is simplified as a diagonal form, whose elements can be written as
Since the diagonal Hessians in Equation 12 require computing the Green's function for every receiver, it is still expensive for dense arrays. Tang (2009) proposed a phase-encoding strategy to compute the receiver wavefield using a composite source
where w(x s , x r ) is the source and receiver mask operator and 휙(x r , 휔) is a random phase-encoding function.
Using this strategy, the Hessians in Equation 12 can be simplified to
where R(x s , x, t) is the receiver wavefield computed using the composite source in Equation 13. As analyzed by Tang (2009), multiple implementations for computing the receiver wavefields are needed to suppress the crosstalk artifacts caused by the phase encoding. Although this is more expensive than one normal wavefield simulation, its computational cost is much lower than that of Equation 12. Ten implementations are used to computed the phased encoding Hessian in subsequent numerical examples.
Under current computational limitations, it is too expensive to directly compute the Hessian and its inverse in FWI applications. Métivier et al. (2013Métivier et al. ( , 2014) ) proposed a truncated Gauss-Newton and Newton method to avoid explicitly computing the Hessian. This method includes a nested loop in numerical implementation: an outer loop for solving a nonlinear problem to update model parameters and an inner loop for solving a linear problem to reduce Hessian effects and optimize gradients. When a large iteration number is used in the inner loop, this approach is still expensive in real applications (Pan et al., 2017). Another popular Hessian-free method is the L-BFGS algorithm (Byrd et al., 1994;Nocedal, 1980;Nocedal & Wright, 2006), which has a good trade-off between computational efficiency and accuracy for large-scale inverse problems.
To accelerate the convergence, we use a preconditioned L-BFGS approach in the simultaneous inversion for P wave velocity and attenuation. The diagonal Hessian in Equation 14 is used as a preconditioner for the misfit gradients. The optimized update direction is equivalent to the solution of the following linear problem
10.1029/2019JB019129 where
T is the model parameter vector with size N × 1, N is the number of the discrete model parameters, 훿m = [훿 ln v, 훿 ln Q] T is the optimized update direction, and 훿 ln v = 훿v∕v 0 and 훿 ln Q = 훿Q∕Q 0 are the relative model perturbations. It is notable that the preconditioned L-BFGS framework approximates H D H(m k ) instead of H as the traditional L-BFGS method. But mathematically, they have the same solution
The preconditioner in Equation 15 can reduce the conditioner number of the Hessian matrix and therefore accelerate the convergence speed of the iterative solver. Because the absolute perturbations of velocity and Q might be quite different, we use a relative update scheme in the multiparameter inversion:
where Diag(m) denotes a N × N matrix with vector m as the diagonal, I is a vector with all entries of one, 훼 is a scalar step length, and k and k + 1 denote the current and next iterations, respectively. The detailed computational steps can be summarized as (1) computing the gradient g = ∇J = [g v , g Q ] T using the adjoint-state method as described in the previous section;
(2) preconditioning the gradient with the inverse of the phase-encoding diagonal Hessian H D = [H vv , H QQ ] T to improve subsurface illumination; (3) computing the update direction 훿m = [훿m v , 훿m Q ] T using the L-BFGS method, with preconditioned gradients and models at the current and the previous l iterations as input, where l is the history length; (4) normalizing the velocity and Q update directions with their respective maximum values; and (5) computing a scalar step length 훼 using a parabolic linear-search algorithm and updating velocity and Q models based on Equation 17. In the L-BFGS algorithm for the numerical examples, the history length l is set as five. Theoretically, the L-BFGS method approximates the full Hessian and therefore takes both first-and high-order scattering effects into account. This makes it better to reduce interparameter crosstalk than the Gauss-Newton method in practical applications, especially for models with large velocity contrasts.
In this section, we use four numerical examples to illustrate the performance of the simultaneous inversion method for estimating P wave velocity and attenuation models. In the numerical implementation, the preconditioner for the conjugate gradient and L-BFGS algorithms is the same, that is, the phase-encoding diagonal Hessian in Equation 12. Ten implementations are used in the phase-encoding diagonal Hessian to suppress the crosstalk artifacts between receivers. The complex-valued forward and adjoint wave equations are solved using the staggered-grid finite-difference and Fourier transform (Yang & Zhu, 2018a).
To illustrate the crosstalk effects in multiparameter inversion, we use an uncorrelated Camembert model (Figure 4) in this example. The discrete model size is a 201 × 201 grid with a 15 m spacing. A Ricker wavelet with the peak frequency of 10 Hz is used as the source time function. The Camembert perturbations for attenuation (훿Q∕Q 0 ) is 50%, which is about 3.8 times those of P wave velocity perturbations (i.e, 훿v∕v 0 =13%). Forty-two sources are evenly deployed along left and top boundaries, and 201 receivers are deployed along right and bottom boundaries. The source-receiver configuration is similar to Figure 1a. The record duration is 2 s, and the time sample increment is 4 ms. Two representative shots are displayed in Figures 4c and 4d.
The imaginary parts of the seismograms are computed using the Hilbert transform (the blue dashed line in Figure 4e).
The estimated P wave velocity and Q models using different FWI schemes are presented in Figure 5. The initial model for the P wave velocity is a constant of 3 km/s, and Q is a constant of 100. Thirty iterations are performed for each inversion scheme. The first scheme is single-parameter inversion. When estimating the P wave velocity model, the true Q model is used as the input and is not updated during iterations. Similarly, the true velocity model is used as the input when estimating the Q model. The preconditioned conjugate gradient (PCG) algorithm is used in the single-parameter inversions. Since there is no interparameter crosstalk problem for the single-parameter inversion, both P wave velocity and attenuation models can be accurately recovered (Figures 5a and 5d). Because the Hessian is not used to optimize the gradients, the edges of the estimated P wave velocity and Q models are blurry. The second inversion scheme is a simultaneous inversion using the same PCG solver. The constant velocity and Q are used as the input, and both of them are updated simultaneously during iterations. Although P wave velocity perturbations are recovered at correct locations, there are weak crosstalk artifacts at the Q perturbation locations (Figure 5b). Strong crosstalk artifacts occur in the estimated Q model at the velocity Camembert locations, and the Q values at the correct locations are not accurately recovered (Figure 5e). The last inversion scheme is a simultaneous inversion using the preconditioned L-BFGS method. Both velocity and Q Camemberts are recovered with correct spatial patterns and physical values (Figures 5c and 5f). The evolution history of the updated velocity and Q models in Figure 6 demonstrates that at early stages, there are strong P wave velocity leakage to the Q model (e.g., the second iteration). This is because the L-BGFS algorithm does not provide a good approximation for the Hessian inverse for a few iterations. With iteration increasing, the approximation accuracy of the L-BFGS solver for the Hessian is improved, which helps to enhance the spatial resolution of P wave velocity and reduce the crosstalk artifacts in the attenuation model. Theoretically, the L-BFGS is a rank 2 updating scheme (Liu & Nocedal, 1989). The good performance for reducing interparameter crosstalk between velocity and Q in this example indicates that approximating the Hessian inverse using a low-rank method might be more effective than directly approximating the Hessian matrix as illustrated in the previous section.
For comparison, we compute the simultaneous inversion results from the L-BFGS method without precondition (Figure 7). At the thirtieth iteration, the P wave velocity Camemberts are resolved. But the attenuation model has strong crosstalk artifacts, which has as large Q values as the true Q Camemberts. Increasing iteration to 100, the crosstalk artifacts are significantly reduced, and the correct Q Camembert structures are recovered. The estimated velocity and attenuation are similar to those computed using the preconditioned L-BFGS method with 30 iterations. This indicates that the diagonal Hessian preconditioner enables us to accelerate the convergence speed and reduce computational costs. In addition, the PCG, classical L-BFGS, and preconditioned L-BGFS produce good velocity constraints with 30 iterations, suggesting velocity inversion is more robust than attenuation. The incorporation of the Hessian into FWI is required to reduce the interparameter crosstalk and produce accurate subsurface attenuation structures.
Next, we use a cross-well example (Figure 8) to test the performance of the proposed simultaneous inversion for transmission data. The model size is 497×225 with a 15 m grid spacing. There is a low-velocity and low-Q body in the middle of the model. Seventy-one sources (the stars in Figure 8a) are evenly distributed along the left boundary. Each source is recorded by 497 receivers deployed along the right boundary (triangles in Figure 8a). The source time function is a Ricker wavelet with the peak frequency of 5 Hz. The record duration is 1 s, and the time sampling increment is 2 ms.
The initial models (Figures 9a and 9d) are generated by smoothing the true velocity and Q models with a 1, 200×450 m Gaussian filter. The simultaneous inversion results using the PCG and preconditioned L-BFGS 10.1029/2019JB019129 Figure 9. Initial models (a, d) and FWI results using the PCG solver (b, e) and the preconditioned L-BFGS method (c, f) for the cross-well example. The top row is for P wave velocity, and the bottom row is for Q model.
solvers after 80 iterations are presented in Figure 9. The estimated P wave velocity models are similar from these two inversion methods (Figures 9b and 9c). The low-velocity body in the middle of the model is accurately recovered, and the layer boundaries are well defined. The accurate recovery for velocity reduces most of the data residuals (Figure 10). Without using the Hessian information, the estimated Q model in Figure 9e has a low resolution, and the recovered Q value in the central low-Q body is larger than its true value. This inaccuracy results in larger data residuals at far offsets than near offsets. In contrast, the preconditioned L-BGFS method produces more accurate attenuation distributions and better data fitting results (Figure 9f). This is because the approximated Hessian can reduce velocity and Q crosstalk and produce more balanced updates for Q than the PCG solver. In terms of the layer interfaces, the estimated Q model is not defined as well as the velocity model. One explanation is that the dispersion and dissipation in seismograms are caused by integrating Q effects along wave paths, and the sharp contrasts in Q distribution produce much weaker scattering than those in the velocity model. Therefore, the macro attenuation models can be effectively estimated using FWI, but the high-wavenumber perturbations are difficult to recover in comparison to velocity models.
Following the work of Wilcock et al. (1992) and Han et al. (2014), we design a mid-ocean ridge model (Figure 11). The layered structures in P wave velocity model represent the oceanic crust strata, including the shallow basalts and deep gabbros. The low-velocity anomaly corresponds to a magma chamber. The attenuation model contains a major axial magma lens (AML) and two small off-axis magma lenses (OAMLs), which have Q value as low as 35 (Figure 11b). We discretize this model with a 270 × 1, 272 grid and a 12.5 m spacing. Sixty-three shots are evenly deployed along the surface with a spacing of 250 m. Each shot is recorded by 636 receivers at the depth of 12.5 m with a horizontal interval of 25 m. The Ricker wavelet with the peak frequency of 12 Hz is used as the source time function. The observed data set (Figure 12a) is simulated by solving the SLS wave equation using the staggered-grid finite-difference with eighth-order accuracy in space and second-order accuracy in time. Three mechanisms with the reference frequencies of 2, 8, and 14 Hz are used in the SLS method to approximate the frequency-independent attenuation effect. For comparison, we compute the viscoacoustic data (Figure 12b) using Equation 2. The waveform consistency (Figures 12c and 12d) between these two data sets verifies the accuracy of Equation 2 for modeling wave propagation in viscoacoustic media.
The initial velocity model (Figure 13a) is generated by smoothing the true model with a Gaussian filter with the width of 375 m × 5km, and a 1-D initial Q model is used in the simultaneous FWI (Figure 13b). A multiscale inversion scheme with five frequency ranges, that is, 1-3, 1-5, 1-7, 1-10, and 1-12 Hz, is used to avoid the cycle-skipping problem. Twenty iterations are performed for each frequency band. The resulting P wave velocity and attenuation models are presented in Figures 13c and 13d. Acoustic FWI is also applied to the same data set, and the result is shown in Figure 13e. To quantitatively evaluate the inversion results, a relative model error is defined as
where m is the true model, m ′ denotes the initial or estimated model, and |||| 2 is the L 2 norm. The relative errors of viscoacoustic and acoustic FWIs are presented in Table 1, and their correlations with the true model are shown in Figure 14. The Hessian approximated with the preconditioned L-BFGS reduces the contamination of interparameter crosstalk and produces good velocity and Q estimations. The fine layers in the velocity model are accurately recovered (Figures 13c and 13f), and the OAMLs and AML are resolved with correct Q values (Figures 13d and 13f). Most data residuals are reduced, and synthetic data computed using the FWI results match observed data well (Figure 15). Compared with P wave velocity, the estimated Q distributions Note. Q-FWI and AC-FWI denote the results computed by viscoacoustic and acoustic FWI, respectively. have a larger model error (Table 1), despite the spatial patterns for the AML and OAMLs are resolved. Their correlations with the true model are not as focused as the estimated velocity model (Figure 14). This might be because seismic data are more sensitive to subsurface variations in velocity than attenuation. Without incorporating attenuation into the wave equation, acoustic FWI cannot simulate physical dispersion and amplitude loss in the propagator. This makes it difficult to match viscoacoustic data (Figures 15e and 15f) and produce inaccurate velocity estimations, especially at large depths (Figure 13e, Table 1, and the green squares in Figure 14a).
Finally, the SEG/EAGE overthrust model (Figure 16) is used to test the robustness of the simultaneous inversion for complex structures. Initial P wave velocity and Q models are presented in Figures 16c and 16d. The model size is 187×801, and the spatial increment is 12.5 m. Eighty sources are uniformly deployed along the surface at the depth of 12.5 m. A fixed spread with 801 receivers at the surface is used to record pressure data. The source time function is a Ricker wavelet with peak frequency of 12 Hz. The record duration is 3 s, and the time interval is 4 ms.
Figures 16e and 16f are the estimated P wave velocity and attenuation models using viscoacoustic FWI. Five frequency ranges are used in the multiscale inversion: 1-3, 1-5, 1-7, 1-9, and 1-12 Hz. Twenty iterations are performed in each frequency range. The diagonal Hessian and gradients at the first iteration of the 1-3 Hz stage are shown in Figure 17. Because of geometric spreading and double-damping, the velocity and Q gradients (Figures 17c and 17d) have unbalanced subsurface illumination and weak amplitudes at great depths.
The precondition with the diagonal Hessian (Figures 17a and 17b) enhances the deep illumination and produces balanced amplitudes (Figures 17e and 17f). The overthrust-associated faults and folds are resolved well in the estimated velocity model (Figure 16e). The low-Q layers are recovered in the estimated attenuation model (Figure 16f), but it has a lower spatial resolution than the estimated velocity model. Comparisons 10.1029/2019JB019129 between observed and predicted data are shown in Figure 18. The initial velocity and attenuation models mainly produce direct and large-offset turning waves, which cannot match observed data in terms of traveltimes and amplitudes (Figures 18a and 18c). In contrast, the seismograms computed using the estimated velocity and Q models fit observed data well for both reflections and refractions. The significant reduction of the data residuals (Figures 18b and 18d) indicates viscoacoustic FWI converges to a good solution.
We demonstrate the possibility of using viscoacoustic FWI to simultaneously map P wave velocity and attenuation distributions in the subsurface. Since attenuation gives rise to phase dispersion and amplitude dissipation in seismograms, incorporating attenuation into FWI produces more accurate velocity models. In addition, estimated Q structures can also be used to delineate specific geological structures, such as unsaturated gas reservoirs, AMLs, and OAMLs beneath mid-ocean ridges. Explicit incorporation of Q in the wave equation and the approximated Hessian used in the proposed method successfully reduce interparameter crosstalk and produce accurate P wave velocity and Q estimates.
Synthetic examples have been used to illustrate the performance of the proposed method. Several critical steps should be considered in practical applications. First, except for attenuation, elasticity and anisotropy are two other important factors that affect wave propagation inside the Earth. Elasticity introduces wave-mode conversions, which produces amplitude differences from acoustic propagation. Anisotropy results in different wave speeds along different directions, which affects both traveltimes and amplitudes of seismograms. Theoretically, these effects should also be included in waveform inversion. But incorporating more parameters expands the search space and increases the nonlinearity of the problem. To mitigate the amplitude effects from elasticity and anisotropy in viscoacoustic FWI applications, one strategy is to use the misfit functions based on the frequency-shift (Dutta & Schuster, 2016;Quan & Harris, 1997) and spectral-amplitude ratio (Karaoglu & Romanowicz, 2017). This can be introduced into the proposed simultaneous inversion by modifying the misfit and associated adjoint sources. Second, cycle skipping is a common problem in FWI, when the initial model is far from the true solution (Virieux & Operto, 2009). Enlarging acquisition apertures and recording low-frequency data (less than 3 Hz) can effectively mitigate this problem but with increasing costs. As analyzed by Malinowski et al. (2011), because attenuation has significant effects after wave propagating for a long time (commonly larger than 3 times the target depth), large-offset data are important to map attenuation distribution in the subsurface. The requirement for low-frequency data can be partially relaxed by introducing novel misfit functions, such as waveform envelopes (Bozdag et al., 2011;Karaoglu & Romanowicz, 2018;Wu et al., 2014) and Wasserstein distances (Engquist et al., 2016;Métivier et al., 2016aMétivier et al., , 2016b;;Yang & Engquist, 2016b). Third, nonuniqueness gives multiple solutions that match the observed data equally well. Combining the inversion results from gravity and electromagnetic data might help to mitigate the nonuniqueness problem.
The amplitude loss caused by seismic attenuation leads to weak signals, which might be contaminated by ambient noises. Directly fitting these noisy data using waveform inversion will result in artifacts in the estimated model parameters (Karaoglu & Romanowicz, 2017). Therefore, careful preprocessing, including muting, trace killing, weak signal enhancement, denoisings, source-receiver calibration, and multiple attenuation, is necessary to ensure the success of viscoacoustic FWI in field data applications (Baeten et al., 2013;Malinowski et al., 2011;Ratcliffe et al., 2011;Prieux et al., 2013;Vigh et al., 2010). Since the forward and adjoint wave equations are solved using the complex-valued finite-difference and Fourier transform (Yang & Zhu, 2018a, 2018b), the computational cost is about 3 times that of conventional acoustic wave equation solvers. Therefore, high-performance clusters and optimized parallel computational algorithms are required for large-scale 3-D experiments.
In the numerical experiments, we mainly use compressional waves to constrain attenuation structures for hydrocarbon reservoir and the uppermost crust. At regional and global scales, strong seismic attenuation has also been observed, especially for shear waves (Cammarano & Romanowicz, 2008;Dalton et al., 2009;De Siena et al., 2010;Nuttli, 1973;Nicolas et al., 1982;Pujades et al., 1997;Romanowicz & Mitchell, 2007;Solomon, 1972;Zhu et al., 2013). The proposed inversion strategies can be extended to regional and global FWI to mitigate velocity and attenuation crosstalk and map attenuation structures in the crust and mantle. Both surface and body waves are commonly used in regional and global FWI, and their sensitivities to 10.1029/2019JB019129 attenuation perturbations might be different from compressional waves. The applications of the proposed attenuation inversion framework to regional and global FWI need further investigations.
We present a viscoacoustic FWI framework for simultaneously mapping the distributions of P wave velocity and attenuation in the subsurface, based on a time domain complex-valued wave equation. Explicit incorporation of Q in the wave equation enables us to directly derive sensitivity kernels with respect to P wave velocity and Q. It does not require transferring Q to stress and strain relaxation times as in the SLS-based method. By numerically analyzing the Gauss-Newton Hessian, we observe strong interparameter crosstalk, especially for the leakage from P wave velocity to Q. To save computational costs, we use a preconditioned L-BGFS method to approximate the Hessian in viscoacoustic FWI. Numerical examples demonstrate that the proposed method can effectively reduce the crosstalk artifacts and produce accurate P wave velocity and Q estimations. Since seismograms are not sensitive to high-wavenumber attenuation structures, the estimated Q distribution has lower spatial resolution than the estimated velocity models.
In this appendix, we give a detailed derivation for the adjoint viscoacoustic wave equation and sensitivity kernels using the Lagrange multiplier method (Liu & Tromp, 2006;Tromp et al., 2005). The augmented misfit function can be written as
where x r is the receiver location, 휆(x, t) is a complex-valued multiplier, superscript * denotes complex conjugate, and ℜ denotes the real part. d syn (x r , t) and d obs (x r , t) are the calculated and observed data, respectively. Because the wavefields are complex-valued, we measure the misfit using the real part of the conjugate multiplication of the data residual. The imaginary parts of the observed data are computed using the Hilbert transform.
In this study, we do not consider the density perturbations. Therefore, taking the variation of the misfit function in Equation A1 and neglecting the high-order terms of the perturbed model parameters and wavefield yield
) p(x, t)
) 훿p(x, t)
where 휌 0 , v 0 and Q 0 are the background density, P wave velocity, and attenuation, respectively. 훿p is the perturbed wavefield, 훿 ln v = 훿v∕v 0 and 훿 ln Q = 훿Q∕Q 0 are the relative velocity and Q perturbations, and 훿(x) is the Kronecker delta function.
According to the property of the complex conjugate, we have the relation as
ℜ{휆 * (x, t)[휅p(x, t)]} = ℜ{휆(x, t)[휅 * p * (x, t)]}, (A3) where 휅 is a complex coefficient of the wavefield p(x, t). Using this relation, Equation A2 can be rewritten as 훿J = ℜ { ∫ T 0 ∫ Ω 훿d s푦n (x) * [ d s푦n (x, t) -d obs (x, t) ] 훿(x -x r ) -[ -2훿 ln v 휌 0 (x)v 2 0 (x) ( C 0 1 (x) 휕 2 휕t 2 -iC 0 2 (x) 휕 휕t + C 0 3 (x) ) p(x, t)휆 * (x, t) + -훿 ln Q 휌 0 (x)v 2 0 (x) ( -2a 휋Q 0 (x) 휕 2 휕t 2 -iC 0 2 (x) 휕 휕t + C 0 3 (x) ) p(x, t)휆 * (x, t) + 1 휌 0 (x)v 2 0 (x) ( C 0 1 (x) 휕 2 휕t 2 + iC 0 2 (x) 휕 휕t + C 0 3 (x) ) 훿p * (x, t)휆(x, t) + (훿 ln v + 훿 ln Q) C 0 4 (x) √ -∇ 2 휕p(x, t) 휕t 휆 * (x, t) -C 4 (x) √ -∇ 2 휕훿p * (x, t) 휕t 휆(x, t) -∇ • ( 1 휌(x) ∇훿p * (x, t) ) 휆(x, t) ] d 3 xdt } . (A4) Using the integration by part, we have the following equality ∫ T 0 휕훿p * (x, t) 휕t 휆(x, t)dt = [훿p * (x, t)휆(x, t)]| T 0 -∫ T 0 휕휆(x, t) 휕t 훿p * (x, t)dt. (A5) Considering the starting condition 훿p(x, t) = 0 and 휕훿p(x, t)∕휕t = 0, as well as the end condition 휆(x, t) = 0 and 휕휆(x, t)∕휕t = 0, Equation A6 reduces to ∫ T 0 휕훿p * (x, t) 휕t 휆(x, t)dt = -∫ T 0 휕휆(x, t) 휕t 훿p * (x, t)dt. (A6) Applying similar relations as in Equation A6, Equation A4 can be simplified as 훿J = ℜ { ∫ T 0 ∫ Ω 훿d s푦n (x) * [ d s푦n (x, t) -d obs (x, t) ] 훿(x -x r )
-훿p * (x, t)
) p(x, t)휆 * (x, t)
Provided that the multiplier 휆(x, t) satisfies the following equation { 1
(A8)