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/m3, 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.
(a)
(b)
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 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 -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].
(a)
(b)
(c)
(a)
(b)
(c)
(d)
After 30 iteration steps, the attenuation anomaly is largely recovered and the artifact is strongly suppressed as shown in Figure 4. And the L2 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 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 anomaly.
(a)
(b)
(c)
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 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 109, the magnitude of density is approximately 103, 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.
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.