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

1. IntroductionSeismic 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]

[mathematical expression not reproducible] (1)

where p = p(x, t) is the pressure field, v = {[v.sub.x], [v.sub.y], [v.sub.z]} is the particle velocity vector, [kappa] = [kappa](x) = [pc.sup.2] is the bulk modulus, c = c(x) is acoustic velocity, p = p(x) is density, and f = f([x.sub.s], t) represents a point source term at [x.sub.s]. The symbol * 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 [mathematical expression not reproducible]. In practice, a single relaxation mechanism approximation could yield results that are sufficiently accurate [19]. [[tau].sub.[epsilon]] and [[tau].sub.[sigma]] are the strain and stress relaxation times of the SLS, respectively. The unitless variable [tau] = [[tau].sub.[epsilon]]/[[tau].sub.[sigma]] - 1 determines the magnitude of Q. The relationship between Q and the relaxation times is defined as [20]

[[tau].sub.[sigma]] = [square root of [Q.sup.2] + 1 - 1]/2[pi][f.sub.0]Q [[tau].sub.[epsilon]] = [square root of [Q.sup.2] + 1 + 1]/2[pi][f.sub.0]Q (2)

Here, [f.sub.0] is the reference frequency usually chosen to be the central frequency of the source wavelet. Since the time derivative of relaxation function is [mathematical expression not reproducible], (1) can be simplified to

[partial derivative]p/[partial derivative]t = -[kappa] (1 + [tau]) [nabla] * v - r + f [partial derivative]v/[partial derivative]t = - 1/[rho] [nabla] p, (3)

with memory variable [15, 20] defined as

[mathematical expression not reproducible] (4)

Taking the derivative of with respect to time yields

[partial derivative] = -[kappa] [tau]/[[tau].sub.[sigma]] [nabla] * v = 1/[[tau].sub.[sigma]] r.(5)

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.sub.max] = [square root of ([[tau].sub.[epsilon]][kappa]/[[tau].sub.[sigma]][rho])] [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 = [([kappa], [rho], [tau]).sup.T] which can explain the observed data very well. It can be achieved by minimizing the misfit energy between the observed data d([x.sub.r], t) and the predicted data p([x.sub.r], t):

J (m) = 1/2 [summation over (R)] [[integral].sup.T.sub.0] [[p([x.sub.r], t) - d ([x.sub.r], t)].sup.2] dt, (6)

where the interval [0, T] denotes the time series of interest and [x.sub.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, [partial derivative]J(m)/[partial derivative]m. To simplify the problem, the assumption is made that the material relaxation parameter [[tau].sub.[sigma]] 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.sub.r], t) satisfies the viscoacoustic wave equations as discussed above, the PDE-constrained optimization problem is obtained by the method of Lagrange multipliers:

[mathematical expression not reproducible] (7)

where [OMEGA] is the domain of integration, [partial derivative][OMEGA] is the surface of [OMEGA], and [[lambda].sub.1], [[lambda].sub.2], [[lambda].sub.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:

[mathematical expression not reproducible] (8)

Here, n is the unit vector pointing outward normal to the surface [partial derivative][OMEGA]. The wave fields p, r, and v are usually called state variables that are subject to the initial condition:

p (x,0) = 0, V (x, 0) = 0, r (x, 0) = 0 (9)

and the radiation boundary condition [27]:

p (x, t)[|.sub.x[right arrow][infinity]] = 0, v (x, t) [|.sub.x[right arrow][infinity]] = 0. (10)

The Lagrange multiplier wave fields [[lambda].sub.1], [[lambda].sub.2], and [[lambda].sub.3] are constrained by the final condition:

[[lambda].sub.1] (x, T) = 0, [[lambda].sub.2] (x, T) = 0, [[lambda].sub.3] (x, T) = 0 (11)

and the boundary condition:

[[lambda].sub.1] (x, t) [|.sub.x[right arrow][infinity]] = 0, [[lambda].sub.2] (x, t) [|.sub.x[right arrow][infinity]] = 0, [[lambda].sub.3] (x, t) [|.sub.x[right arrow][infinity]] = 0. (12)

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 [[lambda].sub.1], [[lambda].sub.2], and [[lambda].sub.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:

[mathematical expression not reproducible] (13)

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 [[summation].sub.r](p([x.sub.r], t) - d([x.sub.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, [partial derivative]X/[partial derivative]m = [partial derivative]j/[partial derivative]m according to [17]. Therefore, the gradients of the objective function J, with respect to the model parameters, modulus [kappa], the density [rho], and attenuation-related parameter [mu], can be written as

[mathematical expression not reproducible] (14)

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 [[lambda].sub.1], [[lambda].sub.2], and [[lambda].sub.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 maybe applied to perform an inversion for bulk modulus, density, and [tau]. It is preferred to use the conjugate gradient method to help to speed up convergence:

[c.sub.n] = [g.sub.n] + [beta][c.sub.n-1]. (15)

Here, [c.sub.n-1] and [c.sub.n] are the last and current conj'ugate gradients; [g.sub.n-1] and [g.sub.n] are the last and current steepest decent directions, respectively. The weighting factor [beta] is given by the PolakRibiere method [28]:

[beta] = [g.sup.t.sub.n]([g.sub.n] - [g.sub.n-1]/ [g.sup.t.sub.n-1][g.sup.n-1]. (16)

For the first iteration, [beta] = 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.sub.r], t) and forward wave fields;

(ii) calculate the residuals p([x.sub.r], t) - d([x.sub.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 [mu] via line search algorithm.

(d) Update the model [m.sub.n+1] = [m.sub.n] -[mu][c.sub.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 [tau] is first inverted and then transformed to Q model through [7]

Q = [square root of ([(2/[tau] + 1).sup.2] - 1)].(17)

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 x 260 m and consists of a homogeneous full space with a velocity of 2000 m/s, a density 1500 kg/[m.sup.3], 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.

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 r. 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].

After 30 iteration steps, the attenuation anomaly is largely recovered and the artifact is strongly suppressed as shown in Figure 4. And the [L.sub.2] 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.

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.sup.9], the magnitude of density is approximately [10.sup.3], and the magnitude of r is just about [10.sup.-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.

http://dx.doi.org/10.1155/2016/4750438

Competing Interests

The authors declare that they have no competing interests.

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 and P. G. Richards, Quantitative Seismology, W. H. Freeman, San Francisco, Calif, USA, 1980.

[2] Z. Wang, J. Gao, Q. Zhou, K. Li, and 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 and 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, and 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 and D. Yingst, "Q estimation through waveform inversion," in Proceedings of the 75th Annual International Conference and Exhibition (EAGE '13), 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 and 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, and Y. Dou, "Seismic waveform inversion using the finite-difference contrast source inversion method," Journal of Applied Mathematics, vol. 2014, Article ID 532159, 11 pages, 2014.

[14] J. M. Carcione, D. Kosloff, and 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, and W. W. Symes, "Viscoelastic finite-difference modeling," Geophysics, vol. 59, no. 9, pp. 1444-1456, 1994.

[16] G. Dutta and 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 and 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, and 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 and Y. Zhang, "Staggered-grid finite difference method with variable-order accuracy for porous media," Mathematical Problems in Engineering, vol. 2013, Article ID 157071, 10 pages, 2013.

[23] A. S. Serdyukov and A. A. Duchkov, "Hybrid kinematicdynamic approach to seismic wave-equation modeling, imaging, and tomography," Mathematical Problems in Engineering, vol. 2015, Article ID 543540, 8 pages, 2015.

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

[25] R. Martin and 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, and Y.-H. Zeng, "The adjoint method for the inverse problem of option pricing," Mathematical Problems in Engineering, vol. 2014, Article ID 314104, 7 pages, 2014.

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

[28] E. Polak and G. Ribiere, "Note sur la convergence de methodes de directions conjuguees," Revue Frangaise d'Informatique et de Recherche Operationnelle, vol. 3, no. 16, pp. 35-43,1969.

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

Guowei Zhang (1, 2) and Jinghuai Gao (1, 2)

(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

Correspondence should be addressed to Jinghuai Gao; jhgao@mail.xjtu.edu.cn

Received 2 February 2016; Accepted 17 May 2016

Academic Editor: Yannis Dimakopoulos

Caption: 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.

Caption: 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.

Caption: FIGURE 3: Gradients of the Q-related parameter r 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.

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

Caption: 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.

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

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Zhang, Guowei; Gao, Jinghuai |

Publication: | Mathematical Problems in Engineering |

Date: | Jan 1, 2016 |

Words: | 4043 |

Previous Article: | Multiscale Fluctuation Features of the Dynamic Correlation between Bivariate Time Series. |

Next Article: | The Optimization Model of Earthquake Emergency Supplies Collecting with the Limited Period and Double-Level Multihub. |

Topics: |