Mathematical Problems in Engineering

Volume 2016, Article ID 4750438, 8 pages

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

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

^{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

Academic Editor: Yannis Dimakopoulos

Copyright © 2016 Guowei Zhang and Jinghuai Gao. 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.

#### 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 . 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 . For the 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 -related parameter . Numerical tests on the crosswell recording geometry indicate the feasibility of the proposed approach for the anomaly estimation.

#### 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 . It is assumed that 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 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 -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 [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 into account, the superposition of several standard linear solid (SLS) models is introduced to approximate a frequency-independent 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 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]where is the pressure field, is the particle velocity vector, is the bulk modulus, is acoustic velocity, is density, and represents a point source term at . The symbol represents the time convolution and 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 . In practice, a single relaxation mechanism approximation could yield results that are sufficiently accurate [19]. and are the strain and stress relaxation times of the SLS, respectively. The unitless variable determines the magnitude of . The relationship between and the relaxation times is defined as [20]Here, is the reference frequency usually chosen to be the central frequency of the source wavelet. Since the time derivative of relaxation function is , (1) can be simplified towith memory variable [15, 20] defined as Taking the derivative of with respect to time yieldsCombining (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 [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 which can explain the observed data very well. It can be achieved by minimizing the misfit energy between the observed data and the predicted data :where the interval denotes the time series of interest and 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, . 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 Fréchet 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 satisfies the viscoacoustic wave equations as discussed above, the PDE-constrained optimization problem is obtained by the method of Lagrange multipliers:where is the domain of integration, is the surface of , and , , 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:Here, is the unit vector pointing outward normal to the surface . The wave fields , , and are usually called state variables that are subject to the initial condition:and the radiation boundary condition [27]:The Lagrange multiplier wave fields , , and are constrained by the final condition:and the boundary condition: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 , , and 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 , , and , the adjoint equations are obtained as follows: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 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, according to [17]. Therefore, the gradients of the objective function , with respect to the model parameters, modulus , the density , and attenuation-related parameter , can be written asNow 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 and . And adjoint equation (13) can be solved for the adjoint wave fields , , and 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 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:Here, and are the last and current conjugate gradients; and are the last and current steepest decent directions, respectively. The weighting factor is given by the Polak-Ribière method [28]:For the first iteration, . 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 and forward wave fields;(ii)calculate the residuals ;(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 .(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 model, the quality factor 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 anomaly with crosshole recording geometry in 2D is investigated. The -related parameter is first inverted and then transformed to model through [7]A simple problem used for testing elastic waveform inversion method for reconstructing velocity [29] is modified to adapt to the purpose of inversion. The spherical low 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 anomaly starting with an initial homogeneous model. As shown in Figure 1(a), the test problem has the dimensions m and consists of a homogeneous full space with a velocity of 2000 m/s, a density 1500 kg/m^{3}, and a of 100, respectively. A spherical anomaly with a diameter of 20 m is embedded in the center of model. It has a 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.