Time Domain Waveform Inversion for the Q Model Based on the First-Order Viscoacoustic Wave Equations

Abstract

Propagating seismic waves are dispersed and attenuated in the subsurface due to the conversion of elastic energy into heat. The absorptive property of a medium can be described by the quality factor Q . In this study, the first-order pressure-velocity viscoacoustic wave equations based on the standard linear solid model are used to incorporate the effect of Q . For the Q model inversion, an iterative procedure is then proposed by minimizing an objective function that measures the misfit energy between the observed data and the modeled data. The adjoint method is applied to derive the gradients of the objective function with respect to the model parameters, that is, bulk modulus, density, and Q -related parameter τ . Numerical tests on the crosswell recording geometry indicate the feasibility of the proposed approach for the Q anomaly estimation.

Guowei Zhang 1,2 and Jinghuai Gao 1,2

Academic Editor:Yannis Dimakopoulos

1, Institute of Wave and Information, Xi'an Jiaotong University, Xi'an 710049, China

2, National Engineering Laboratory for Offshore Oil Exploration, Xi'an 710049, China

Received 2 February 2016; Accepted 17 May 2016

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Seismic wave propagation in the earth has been known to be anelastic due to the conversion of elastic energy into heat. The velocity of waves in anelastic medium is dependent on frequency, which distorts the phases of wavelets. Moreover, anelastic medium decreases the amplitudes of propagating waves [1]. Therefore, these anelastic behaviors on kinematics and dynamics of seismic waves can have significant effects on forward modeling, inversion, and seismic attribute [2].

The anelastic effects on propagating waves are well described by the quality factor Q . It is assumed that Q does not change with frequency of the seismic data acquired for oil and gas exploration [3]. Due to the requirement of obtaining higher resolution images and better reservoir characterization in seismic exploration, many algorithms to estimate Q model from seismic data have been studied, such as spectral ratio method [4] and frequency shift method [3, 5, 6]. These methods are based on variations in spectrum amplitude or frequency characteristic of wavelet, provided that scattering, geometrical spreading, and other non Q -related factors have been removed [7].

Full waveform inversion (FWI), originally developed by Tarantola in the 1980s [8-10], is a promising method to extract material parameters in the underground and has been studied extensively in recent years (see [11] and the references therein). It is implemented by minimizing the misfit energy between the observed and modeled data through a gradient optimization approach. Therefore, it requires tremendous amount of floating point calculations and run times to carry out forward modeling of the wave equation during an iteration. For frequency-domain FWI, the attenuation effect can be easily included by replacing the real-valued velocity with a complex-valued velocity using the relationship between complex velocity and Q [12]. However, frequency-domain forward modeling of wave equations needs to implement a large-scale lower and upper triangular decomposition [11, 13].

In time domain, to take Q into account, the superposition of several standard linear solid (SLS) models is introduced to approximate a frequency-independent Q within a predefined frequency band. As a result, the convolutional constitutive law could be circumvented by a set of first-order linear temporal partial differential equations [14, 15]. Thus the first-order viscoacoustic wave equations can be numerically solved by high-order time domain finite difference method on staggered grids to extrapolate seismic wave fields. Based on this, attenuation compensation for least-squares reverse time migration is implemented [16]. In this paper, the first-order viscoacoustic wave equations are used to estimate subsurface attenuation parameters by waveform inversion in time domain.

The rest of this paper is organized as follows. First of all, the first-order time domain viscoacoustic wave equations are briefly discussed. And the inverse problem is solved by minimizing the least-square misfit between the predicted and the observed seismic data. The adjoint method [17, 18] is applied to derive the gradients of the objective function with respect to the model parameters. Then, the inversion algorithm for the Q model is validated by numerical tests. The final section provides conclusions.

2. Methodology

2.1. Forward Problem: Viscoacoustic Wave Equations

The system of wave equations for a viscoacoustic medium can be written as [15] [figure omitted; refer to PDF] where p = p ( x , t ) is the pressure field, v = ( _{v x} , _{v y} , _{v z} ) is the particle velocity vector, κ = κ ( x ) = ρ ^{c 2} is the bulk modulus, c = c ( x ) is acoustic velocity, ρ = ρ ( x ) is density, and f = f ( _{x s} , t ) represents a point source term at _{x s} . The symbol [low *] represents the time convolution and H ( t ) is the Heaviside function ensuring the causality. In (1), only one standard linear solid relaxation mechanism is used to describe the attenuating and dispersive effects in viscoacoustic medium through the relaxation function ( 1 + τ ^{e - t / _{τ σ}} ) H ( t ) . In practice, a single relaxation mechanism approximation could yield results that are sufficiently accurate [19]. _{τ [straight epsilon]} and _{τ σ} are the strain and stress relaxation times of the SLS, respectively. The unitless variable τ = _{τ [straight epsilon]} / _{τ σ} - 1 determines the magnitude of Q . The relationship between Q and the relaxation times is defined as [20] [figure omitted; refer to PDF] Here, _{f 0} is the reference frequency usually chosen to be the central frequency of the source wavelet. Since the time derivative of relaxation function is ( 1 + τ ) δ ( t ) - ( τ / _{τ σ} ) ^{e - t / _{τ σ}} H ( t ) , (1) can be simplified to [figure omitted; refer to PDF] with memory variable [15, 20] defined as [figure omitted; refer to PDF] Taking the derivative of r with respect to time yields [figure omitted; refer to PDF] Combining (3) and (5) leads to the first-order wave equations to describe wave propagation in a viscoacoustic medium. It is convenient to apply staggered-grid finite difference scheme to implement forward modeling for this kind of first-order wave equations [21-23]. From (1), the second-order viscoacoustic wave equations can be also derived and then be applied to waveform inversion with finite difference on centered grids [24].

To minimize numerical dispersion in staggered-grid finite difference modeling, ten samples per wavelength are required for the lowest velocity in the model. Meanwhile, the maximum velocity to compute the Courant number should be adjusted to the highest phase velocity _{c m a x} = _{τ [straight epsilon]} κ / _{τ σ} ρ [15] to avoid instability during numerical simulation. And, in the context of the simulation of seismic wave propagation, the computational domain is restricted to only a part of the true physical domain because of limited computational resources. Therefore, the unsplit convolutional perfect matched layer technique [25] is applied to absorb the reflected energy from the artificial boundaries.

2.2. Inverse Problem

In general, the aim of waveform inversion is to find an optimal model m = ^{( κ , ρ , τ ) T} which can explain the observed data very well. It can be achieved by minimizing the misfit energy between the observed data d ( _{x r} , t ) and the predicted data p ( _{x r} , t ) : [figure omitted; refer to PDF] where the interval [ 0 , T ] denotes the time series of interest and _{x r} denotes the receiver locations. Since the relationship between the data and the model is nonlinear, the inversion is performed iteratively through a gradient optimization approach. The update direction is determined by the derivatives of the objective function with respect to the model parameters, that is, ∂ J ( m ) / ∂ m . To simplify the problem, the assumption is made that the material relaxation parameter _{τ σ} is constant. As the number of unknowns is large in 2D or 3D problems, it is not efficient to obtain the gradient from the Frechet derivatives computed by differentiating the wave field with respect to each model parameters. The adjoint method is widely used for inverse problems (e.g., [17, 26]). It allows us to efficiently calculate the gradient for the minimization procedure as follows.

Since the predicted wave field data p ( _{x r} , t ) satisfies the viscoacoustic wave equations as discussed above, the PDE-constrained optimization problem is obtained by the method of Lagrange multipliers: [figure omitted; refer to PDF] where Ω is the domain of integration, ∂ Ω is the surface of Ω , and _{λ 1} , _{λ 2} , _{λ 3} are the Lagrange multipliers that remain to be determined. Taking the variation of (7), the following expression is found after integration by parts and application of the Gauss divergence theorem: [figure omitted; refer to PDF] Here, n is the unit vector pointing outward normal to the surface ∂ Ω . The wave fields p , r , and v are usually called state variables that are subject to the initial condition: [figure omitted; refer to PDF] and the radiation boundary condition [27]: [figure omitted; refer to PDF] The Lagrange multiplier wave fields _{λ 1} , _{λ 2} , and _{λ 3} are constrained by the final condition: [figure omitted; refer to PDF] and the boundary condition: [figure omitted; refer to PDF] Thus the last two terms of (8) vanish. According to the first-order optimality conditions, the Lagrangian is stationary at the optimum. As a consequence, the variation of the Lagrangian with respect to _{λ 1} , _{λ 2} , and _{λ 3} should be zeros that yields the state equations, that is, the viscoacoustic wave equations (3) and (5). Then, taking the variation of the Lagrangian with respect to the state variables p , v , and r , the adjoint equations are obtained as follows: [figure omitted; refer to PDF] Note that the adjoint equations are similar to the viscoacoustic wave equations. However, the adjoint wave fields are subject to final time conditions and are generated by propagating the residual data _{∑ r} ( p ( _{x r} , t ) - d ( _{x r} , t ) ) from the receiver positions backward in time. Provided that the Lagrange multipliers are determined by the adjoint equation (13), the first three terms of (8) vanish. Moreover, ∂ χ / ∂ m = ∂ J / ∂ m according to [17]. Therefore, the gradients of the objective function J , with respect to the model parameters, modulus κ , the density ρ , and attenuation-related parameter τ , can be written as [figure omitted; refer to PDF] Now the gradients can be obtained by computing the forward wave fields and the adjoint wave fields. For given current models and actual sources, a forward propagation is implemented through viscoacoustic wave equations (3) and (5) to obtain the forward wave fields v and p . And adjoint equation (13) can be solved for the adjoint wave fields _{λ 1} , _{λ 2} , and _{λ 3} with data residuals regarded as sources. The residuals at each time step are injected backward in time. Thus, it is commonly described as a backpropagation of data residuals [27].

2.3. The Inversion Algorithm for the Q Model

After obtaining the gradients of the model parameters, a gradient-based optimization method may be applied to perform an inversion for bulk modulus, density, and τ . It is preferred to use the conjugate gradient method to help to speed up convergence: [figure omitted; refer to PDF] Here, _{c n - 1} and _{c n} are the last and current conjugate gradients; _{g n - 1} and _{g n} are the last and current steepest decent directions, respectively. The weighting factor β is given by the Polak-Ribière method [28]: [figure omitted; refer to PDF] For the first iteration, β = 0 . Finally, the general inversion algorithm is performed as follows:

(a) Calculate the gradient for each source:

(i) solve the forward problem for a given model, that is, the viscoacoustic wave equations (3) and (5) to generate modeled seismic data p ( _{x r} , t ) and forward wave fields;

(ii) calculate the residuals p ( _{x r} , t ) - d ( _{x r} , t ) ;

(iii): solve the adjoint problem, that is, (13), by propagating the residuals as acting source backward in time to generate the adjoint wave fields;

(iv) compute the gradient through (14).

(b) Stack the results from all sources and then calculate the conjugate gradient using (15) and (16).

(c) Calculate the step length μ via line search algorithm.

(d) Update the model _{m n + 1} = _{m n} - μ _{c n} .

(e) Repeat steps (a) to (d), until the residual energy is smaller than a given value or the maximum iteration number is reached.

And during the update process, it is necessary to apply some additional constraints to obtain physically meaningful model parameters. For the inversion for the Q model, the quality factor Q at each point should not be negative. In addition, some appropriate preconditioning operators are also suggested to be applied to the gradient to accelerate convergence [7, 11, 29].

3. Numerical Results

To validate the proposed method, inversion for the low Q anomaly with crosshole recording geometry in 2D is investigated. The Q -related parameter τ is first inverted and then transformed to Q model through [7] [figure omitted; refer to PDF] A simple problem used for testing elastic waveform inversion method for reconstructing velocity [29] is modified to adapt to the purpose of Q inversion. The spherical low Q anomaly is set in the center of model instead of the low velocity anomaly. The crosshole seismic data is acquired from two parallel boreholes. With these data, the proposed inversion method is applied to reconstruct the low Q anomaly starting with an initial homogeneous Q model. As shown in Figure 1(a), the test problem has the dimensions 160 × 260 m and consists of a homogeneous full space with a velocity of 2000 m/s, a density 1500 kg/m³ , and a Q of 100, respectively. A spherical anomaly with a diameter of 20 m is embedded in the center of model. It has a Q of 20 that means strong attenuation while the velocity and density remain the same. In the left-side borehole, 31 sources are located with a signature of Ricker wavelet whose peak frequency is 125 Hz. A string of 181 receivers are placed in the right borehole. They are uniformly distributed in a depth between 40 m and 220 m. The vertical spacing of source is 6 m while the receiver interval is 1 m, and their horizontal offset is 80 m. This true model is used to generate the "observed" seismic data.

Figure 1: (a) The true model: a spherical anomaly with a Q of 20 is embedded in a homogeneous full space with a Q of 100. The red dot and the blue dot in the left figure denote the source and receiver locations, respectively. (b) The starting Q model has a homogeneous Q of 100.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

To implement the forward modeling and adjoint modeling, a staggered-grid finite difference scheme of second-order accuracy in time and eighth-order accuracy in space is used with an unsplit convolutional perfectly matched layer. To avoid numerical dispersion and instability, the grid spacing and time step size is chosen as 1.0 m and 0.1 ms, respectively. The perfect matched layer has a thickness of 20 m. Therefore, the spatial model consists of 201 grid-points in horizontal direction and 301 grid-points in vertical direction, respectively. The total time step has a number of 1500.

For the iterative inversion, a homogeneous Q model without the anomaly is regarded as the starting model as shown in Figure 1(b). The synthetic data of the starting model, the observed seismic data of the true model, and their residuals for the representative source at depth 130 m are illustrated in Figure 2. As a consequence of the anomaly, the observed data is subject to strong attenuation that leads to the large initial data residual. Propagating the residuals of each shot backward in time from the receiver positions, the adjoint wave fields are obtained to calculate the gradient of Q -related parameter τ . In Figure 3, gradients of the shots at depth 40 m, 130 m, and 220 m and the total gradient that is computed by the superposition of all 31 sources are shown. The gradients for the individual shot are shaped like Banana-Doughnut sensitivity kernels [18, 29]. Moreover, the gradient of all shots shows the shape of the anomaly, but the "X" shaped artifact is surrounding the anomaly which comes from the superposition of these sensitivity kernels. It is also seen that the gradient is contaminated by some high-amplitude artifacts around the acquisition geometry, particularly at the source locations. These artifacts are due to the large amplitude of forward and adjoint wave field nearby and can be mitigated by applying a taper function as a preconditioning operator at source and receiver locations [29].

Figure 2: (a) The synthetic data of the starting Q model, (b) the "observed" seismic data of the true Q model, and (c) their residuals for the representative source at depth 130 m.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

Figure 3: Gradients of the Q -related parameter τ for the different sources. (a) The source at depth 40 m, (b) the source at depth 130 m, (c) the source at depth 220 m, and (d) the total gradient that is computed by the superposition of all 31 sources.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

After 30 iteration steps, the attenuation anomaly is largely recovered and the artifact is strongly suppressed as shown in Figure 4. And the L₂ norm of the final data residuals is reduced to a value of 0.21% of the initial one. In Figure 5, the synthetic data of the inverted model, the observed seismic data of the true model, and their residuals for the same source at depth 130 m are illustrated. It shows that the inverted Q model can reconstruct the observed data quite well. The decrease of the objective function through the inversion is shown in Figure 6. The convergence is fast within a few iteration steps. The results show the feasibility of the proposed approach to reconstruct the Q anomaly.

Figure 4: Inversion result for the Q model after 30 iterations.

[figure omitted; refer to PDF]

Figure 5: (a) The synthetic data of the inverted Q model, (b) the "observed" seismic data of the true Q model, and (c) their residuals for the representative source at depth 130 m.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

Figure 6: The normalized data misfit energy with respect to the number of iteration.

[figure omitted; refer to PDF]

4. Conclusions

The adjoint method provides a computationally efficient way to calculate the gradients of the objective function. In this study, this method is applied to derive a time domain waveform inversion algorithm to estimate the Q model based on the first-order viscoacoustic wave equations. Numerical tests are implemented to demonstrate the feasibility of the proposed method on the crosswell recording geometry in 2D. The attenuation anomaly is well reconstructed with known bulk modulus and density model. Generally, the magnitude of the modulus is about 10⁹ , the magnitude of density is approximately 10³ , and the magnitude of τ is just about 10^-2 . Thus their contributions to the objective function are different. Simultaneous inversion for all kinds of model parameters underground remains to be studied in the future.

Acknowledgments

The authors greatly appreciate the Major Programs of the National Natural Science Foundation of China under Grants nos. 41390450 and 41390454, the Major Research Plan of the National Natural Science Foundation of China under Grant no. 91330204, and Beijing Center for Mathematics and Information Interdisciplinary Science (BCMIIS) for their financial support.

References

[1] K. Aki, P. G. Richards Quantitative Seismology , W. H. Freeman, San Francisco, Calif, USA, 1980.

[2] Z. Wang, J. Gao, Q. Zhou, K. Li, J. Peng, "A new extension of seismic instantaneous frequency using a fractional time derivative," Journal of Applied Geophysics , vol. 98, pp. 176-181, 2013.

[3] C. Zhang Seismic absorption estimation and compensation [Ph.D. thesis] , The University of British Columbia, Vancouver, Canada, 2008.

[4] R. Tonn, "The determination of the seismic quality factor Q from VSP data: a comparison of different computational methods," Geophysical Prospecting , vol. 39, no. 1, pp. 1-27, 1991.

[5] Y. Quan, J. M. Harris, "Seismic attenuation tomography using the frequency shift method," Geophysics , vol. 62, no. 3, pp. 895-905, 1997.

[6] J. Gao, S. Yang, D. Wang, R. Wu, "Estimation of quality factor Q from the instantaneous frequency at the envelope peak of a seismic signal," Journal of Computational Acoustics , vol. 19, no. 2, pp. 155-179, 2011.

[7] J. Bai, D. Yingst, "Q estimation through waveform inversion," in Proceedings of the 75th Annual International Conference and Exhibition (EAGE '13), of Extended Abstracts, Th-10-01, 2013.

[8] A. Tarantola, "Inversion of seismic reflection data in the acoustic approximation," Geophysics , vol. 49, no. 8, pp. 1259-1266, 1984.

[9] A. Tarantola Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation , Elsevier, Amsterdam, The Netherland, 1987.

[10] A. Tarantola, "Theoretical background for the inversion of seismic waveforms including elasticity and attenuation," Pure and Applied Geophysics , vol. 128, no. 1, pp. 365-399, 1988.

[11] J. Virieux, S. Operto, "An overview of full-waveform inversion in exploration geophysics," Geophysics , vol. 74, no. 6, pp. WCC1-WCC26, 2009.

[12] Y. Wang Seismic Inverse Q Filtering , Blackwell, Oxford, UK, 2008.

[13] B. Han, Q. He, Y. Chen, Y. Dou, "Seismic waveform inversion using the finite-difference contrast source inversion method," Journal of Applied Mathematics , vol. 2014, 2014.

[14] J. M. Carcione, D. Kosloff, R. Kosloff, "Wave propagation simulation in a linear viscoelastic medium," Geophysical Journal , vol. 95, no. 3, pp. 597-611, 1988.

[15] J. O. A. Robertsson, J. O. Blanch, W. W. Symes, "Viscoelastic finite-difference modeling," Geophysics , vol. 59, no. 9, pp. 1444-1456, 1994.

[16] G. Dutta, G. T. Schuster, "Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation," Geophysics , vol. 79, no. 6, pp. S251-S262, 2014.

[17] R.-E. Plessix, "A review of the adjoint-state method for computing the gradient of a functional with geophysical applications," Geophysical Journal International , vol. 167, no. 2, pp. 495-503, 2006.

[18] Q. Liu, J. Tromp, "Finite-frequency kernels based on adjoint methods," Bulletin of the Seismological Society of America , vol. 96, no. 6, pp. 2383-2397, 2006.

[19] T. Zhu, J. M. Carcione, J. M. Harris, "Approximating constant-Q seismic propagation in the time domain," Geophysical Prospecting , vol. 61, no. 5, pp. 931-940, 2013.

[20] J. M. Carcione Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media , Elsevier Science, Amsterdam, The Netherlands, 2007.

[21] J. Virieux, "P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method," Geophysics , vol. 51, no. 4, pp. 889-901, 1986.

[22] J. Gao, Y. Zhang, "Staggered-grid finite difference method with variable-order accuracy for porous media," Mathematical Problems in Engineering , vol. 2013, 2013.

[23] A. S. Serdyukov, A. A. Duchkov, "Hybrid kinematic-dynamic approach to seismic wave-equation modeling, imaging, and tomography," Mathematical Problems in Engineering , vol. 2015, 2015.

[24] J. Bai, D. Yingst, R. Bloor, J. Leveille, "Viscoacoustic waveform inversion of velocity structures in the time domain," Geophysics , vol. 79, no. 3, pp. R103-R119, 2014.

[25] R. Martin, D. Komatitsch, "An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation," Geophysical Journal International , vol. 179, no. 1, pp. 333-344, 2009.

[26] S.-L. Wang, Y.-F. Yang, Y.-H. Zeng, "The adjoint method for the inverse problem of option pricing," Mathematical Problems in Engineering , vol. 2014, 2014.

[27] N. Kamath, I. Tsvankin, "Gradient computation for elastic full-waveform inversion in 2D VTI media," CWP Report , no. 752, 2013.

[28] E. Polak, G. Ribière, "Note sur la convergence de methodes de directions conjuguees," Revue Française d'Informatique et de Recherche Operationnelle , vol. 3, no. 16, pp. 35-43, 1969.

[29] D. Köhn Time domain 2D elastic full waveform tomography [Ph.D. thesis] , Christian-Albrechts-Universitat zu Kiel, 2011.